NuBase.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Dec/2005 onwards
//                                                                    
// Contact: j.j.hartnell@rl.ac.uk                         

#include "TChain.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1F.h"
#include "TObjString.h"
#include "TMath.h"
#include "TRegexp.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TSystem.h"

#include "MessageService/MsgService.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRSlice.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "StandardNtuple/NtpStRecord.h"

#include "NuMuBar/NuBase.h"
#include "NuMuBar/NuConfig.h"

using std::endl;
using std::cout;
using std::map;
using std::vector;

string NuBase::fInputFileName="";//use hard coded place by default
Bool_t NuBase::fLoadTrees=true;//load trees by default

//ClassImp(NuBase)

CVSID("$Id: NuBase.cxx,v 1.5 2007/12/14 23:36:49 ahimmel Exp $");

//......................................................................

NuBase::NuBase()
{
  MSG("NuBase",Msg::kInfo)
    <<"Running NuBase Constructor..."<<endl;

  fChain=0;
  fChainBD=0;
  fOutFile=0;
  fRec=0;//sr
  fRecBD=0;//sr beam data
  
  fEntries=0;
  fEntriesBD=0;
  fRun=100;//default
  fSubrun=0;

  fNumFiles=0;//number of files added to chain
  
  fFirstFileName="";
  fS="";//general purpose string


  if (fLoadTrees) {
    MSG("NuBase",Msg::kInfo)<<"Loading trees"<<endl;
    this->MakeChain();
    this->SetRootFileObjects();
    //this->InitialiseHitInfoVariables();
    
    //get basic run number stuff
    if (fRec && fChain){//sr
      //get the first event
      fChain->GetEvent(0); 
      NtpStRecord& ntp=(*fRec);
      const RecCandHeader& rec=ntp.GetHeader();
      MSG("NuBase",Msg::kInfo)
        <<"Found: run="<<rec.GetRun()
        <<", subrun="<<rec.GetSubRun()<<endl;
      fRun=rec.GetRun();
      fSubrun=rec.GetSubRun();
    }
    else cout<<"Ahhh, no root files"<<endl;
  }
  else  MSG("NuBase",Msg::kInfo)<<"Not loading trees"<<endl;

  //setup some root stuff
  //include the under and overflow counts
  gStyle->SetOptStat(1111111);
  gStyle->SetOptFit(1111);
  gStyle->SetPalette(1);
  
  MSG("NuBase",Msg::kInfo)
    <<"Finished NuBase Constructor"<<endl;
}

//......................................................................

NuBase::~NuBase()
{
  MSG("NuBase",Msg::kInfo)
    <<"Running NuBase Destructor..."<<endl;
  
  if (fOutFile){
    //this makes histos disappear from the canvases 
    //so do it at the very end
    fOutFile->Close();
  }
  
  MSG("NuBase",Msg::kDebug)
    <<"Finished NuBase Destructor"<<endl;
}

//......................................................................

void NuBase::InputFileName(string f)
{
  if (f!=""){
    MAXMSG("NuBase",Msg::kInfo,1)
      <<"Running with input file name="<<f<<endl;
    fInputFileName=f;
  }
}

//......................................................................

void NuBase::LoadTrees(Bool_t loadTrees)
{
  fLoadTrees=loadTrees;
}

//......................................................................

void NuBase::WriteOutHistos()
{
  //write out the histos to the file, if it's open
  if (fOutFile){
    if (fOutFile->IsWritable()) {
      fOutFile->cd();

      //store the number of files added to the chain      
      Int_t nf=this->GetNumFilesAddedToChain();
      TH1F* hNumFiles=new TH1F("hNumFiles","hNumFiles",15,-5,10);
      hNumFiles->Fill(1,nf);//store as weight, access with GetMaximum

      //get the run number to print it out below
      TH1F* hRun=dynamic_cast<TH1F*>(gROOT->FindObject("hRun"));
      Float_t nr=-1;
      if (hRun) nr=hRun->GetEntries();
      MSG("NuAnalysis",Msg::kInfo)
        <<"Files added to chain="<<hNumFiles->GetMaximum()
        <<", runs found="<<nr<<endl;

      //create a tree and set branches      
      //TTree runInfo("runInfo","runInfo");
      //have to create local variables here
      //I think it's something to do with the fact that they are
      //read in on other trees in this class as well???
      //No, actually not because fSigCorsPerMip is not in a tree
      //don't understand this crazyness!
      //Int_t lSimFlag=0;
      //runInfo.Branch("SimFlag",&lSimFlag,"SimFlag/I");


      //fill and write      
      //runInfo.Fill();
      //runInfo.Write();

      MSG("NuBase",Msg::kInfo)
        <<"Writing histos to: "<<fOutFile->GetName()<<" ..."<<endl;
      fOutFile->Write();
      //fOutFile->Close();//this makes histos disappear from canvases
      //so do it in the destructor (need to make NuBase on heap)
    }
    else {
      MSG("NuBase",Msg::kWarning)<<"File not writable!"<<endl;
    }
  }
}

//......................................................................

void NuBase::OpenOutFile(TString outName, TString opt)
{
        TH1::AddDirectory(true);
  if (fOutFile) {
                fOutFile->Close();
                delete fOutFile;
        }
        fOutFile = new TFile(outName, opt);
        if (!fOutFile->IsOpen()) {
      MSG("NuBase",Msg::kError)<< "File " << outName << " did not open!"<<endl;         
        }
}

//......................................................................

void NuBase::CloseOutFile()
{
  if (fOutFile) fOutFile->Close();
}

//......................................................................

vector<string> NuBase::MakeFileList()
{

  vector<string> fileList;
  
  if (fInputFileName!=""){

    Int_t findGives=fInputFileName.find(".root");
    MSG("NuBase",Msg::kInfo)
      <<"Checking input string for \"*.root\", position in string="
      <<findGives<<endl;

    if (findGives>0){
      //add the file direct to the list
      fileList.push_back(fInputFileName);
    }
    else{//treat the file as a list of root files 
      ifstream inputFile(fInputFileName.c_str());
      
      //check if file exists
      if (inputFile){
        //variables to hold input from file
        string file="";
        
        //read in from the text file and fill objects
        while(inputFile>>file) {
          MSG("NuBase",Msg::kDebug)
            <<"Found input file name="<<file<<endl;
          fileList.push_back(file);
        }
        MSG("NuBase",Msg::kDebug)
          <<"Files names found in txt file="<<fileList.size()<<endl;
      }
      else{
        MSG("NuBase",Msg::kFatal)
          <<endl<<endl
          <<"**********************************************************"
          <<endl<<"Input txt file of file names does not exist!"<<endl
          <<"InputFileName="<<fInputFileName<<endl
          <<"**********************************************************"
          <<endl<<endl<<"Program will exit here"<<endl;
        exit(0);
      }
    }
  }
  //return the fileList if files were found
  if (fileList.size()>0) return fileList;
  
  //Check the env variable to find files
  char* envVariable=getenv("NUDATA");
  if (envVariable==NULL){
    MSG("NuBase",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"Environmental variable NUDATA not set!"<<endl
      <<"Please set NUDATA to the directory containing the"
      <<" ntuple root files"<<endl
      <<"Note: If more than one file is found they will be"
      <<" concatenated and treated as one"<<endl
      <<"*************************************************************"
      <<endl<<endl<<"Program will exit here"<<endl;
    exit(0);
  }

  string sEnv="";
  sEnv=envVariable;
  sEnv+="/*.root";
  MSG("NuBase",Msg::kInfo)
    <<"Looking for *.root files using the env variable"<<endl
    <<"NUDATA="<<sEnv<<endl;    
  fileList.push_back(sEnv);

  return fileList;
}


//......................................................................

Bool_t NuBase::ObjectExistsInFile(TFile* file,
                                  std::string objectName) const
{
  if (!file) {
    MSG("NuBase",Msg::kWarning)
      <<"File does not exist so can't test if object exists!"<<endl;
    return false;
  }

  TList* listOfKeys=file->GetListOfKeys();

  static TFile* lastFile=0;
  if (lastFile!=file){
    MSG("NuBase",Msg::kInfo)
      <<"File contains these keys:"<<endl;
    listOfKeys->Print();
  }
  lastFile=file;

  TObject* ob=listOfKeys->FindObject(objectName.c_str());
  if (ob) {
    MSG("NuBase",Msg::kInfo)
      <<"Requested object with name="<<objectName
      <<" exists in file:"<<endl;
    ob->Print();
    return true;
  }
  else {
    MSG("NuBase",Msg::kInfo)
      <<"Requested object with name="<<objectName
      <<" does NOT exist in file"<<endl;
    return false;
  }
}

//......................................................................

std::string NuBase::GetFirstFileName(std::string wildcardString) const
{
  //check if there is actually a wildcard
  if (!TString(wildcardString.c_str()).MaybeWildcard()) {
    return wildcardString;
  }
  
  //this was ripped off from ROOT's TChain::AddFile

  vector<string> fileList;

  //wildcarding used in name
  TString basename(wildcardString.c_str());

  Int_t dotslashpos = basename.Index(".root/");
  TString behind_dot_root;
  if (dotslashpos>=0) {
    // Copy the tree name specification
    behind_dot_root = basename(dotslashpos+6,basename.Length()-dotslashpos+6);
    // and remove it from basename
    basename.Remove(dotslashpos+5);
  }

  Int_t slashpos = basename.Last('/');
  TString directory;
  if (slashpos>=0) {
    directory = basename(0,slashpos); // Copy the directory name
    basename.Remove(0,slashpos+1);      // and remove it from basename
  } else {
    directory = gSystem->WorkingDirectory();
  }

  const char *file;
  void *dir = gSystem->OpenDirectory(gSystem->ExpandPathName(directory.Data()));

  if (dir) {
    //create a TList to store the file names (not yet sorted)
    TList l;
    TRegexp re(basename,kTRUE);
    while ((file = gSystem->GetDirEntry(dir))) {
      if (!strcmp(file,".") || !strcmp(file,"..")) continue;
      TString s = file;
      if ( (basename!=file) && s.Index(re) == kNPOS) continue;
      l.Add(new TObjString(file));
    }
    gSystem->FreeDirectory(dir);
    //sort the files in alphanumeric order
    l.Sort();
    TIter next(&l);
    TObjString *obj;
    while ((obj = (TObjString*)next())) {
      file = obj->GetName();
      if (behind_dot_root.Length()!=0){
        string fileName=Form("%s/%s/%s",directory.Data(),
                             file,behind_dot_root.Data());
        fileList.push_back(fileName);
        //nf += AddFile(Form("%s/%s/%s",directory.Data(),file,behind_dot_root.Data()),kBigNumber);
      }
      else {
        string fileName=Form("%s/%s",directory.Data(),file);
        fileList.push_back(fileName);
        //nf += AddFile(Form("%s/%s",directory.Data(),file),kBigNumber);
      }
    }
    l.Delete();
  }

  //check if any files were found
  if (fileList.begin()!=fileList.end()){
    cout<<"Used wildcard expansion to find first file name="
        <<*fileList.begin()<<endl;
    return *fileList.begin();
  }
  else {
    return "";//null string
  }
}

//......................................................................

void NuBase::MakeChain()
{
  //get the files to open
  vector<string> fileList=this->MakeFileList();

  TFile* firstFile=0;
  // case with one single file
  if (!TString(*fileList.begin()).MaybeWildcard()) {
    fFirstFileName=(*fileList.begin());//store first filename
    firstFile=TFile::Open((*fileList.begin()).c_str(),"READ");
  }
  else {
    MSG("NuBase",Msg::kInfo)
      <<"Found a wildcard present in list of file names: "
      <<*fileList.begin()
      <<endl<<"Looking for first filename using wildcard..."<<endl;
    string fileName=this->GetFirstFileName(*fileList.begin());
    fFirstFileName=fileName;//store first filename
    if (fileName!="") {
      firstFile=TFile::Open(fileName.c_str(),"READ");
    }
    else {
      MSG("NuBase",Msg::kInfo)
        <<"No files found matching:"<<endl
        <<*fileList.begin()<<", will exit here..."<<endl;
      exit(1);
    }
  }
  
  //create a chain
  //check if the SR ntuple exists
  if (this->ObjectExistsInFile(firstFile,"NtpSt")) {
    fChain=new TChain("NtpSt");
  }
  else {
    MSG("NuBase",Msg::kInfo)
      <<"Not creating NtpSt TChain because it does not exist in file"
      <<endl;
  }
  
  //check in BD ntuple exists
  if (this->ObjectExistsInFile(firstFile,"NtpBDLite")) {
    fChainBD=new TChain("NtpBDLite");//have to load libraries too
  }
  else {
    MSG("NuBase",Msg::kInfo)
      <<"Not creating NtpBDLite TChain because it does not exist"
      <<" in file"<<endl;
  }
  
  firstFile->Close();
  if (firstFile) delete firstFile;
  firstFile=0;

  Int_t nf=0;
  Int_t nfBD=0;
  //add the files to the chain
  for (vector<string>::iterator file=fileList.begin();
       file!=fileList.end();++file){
    
    //test if file exists
    Bool_t openOk=false;
    //can't test dcache files are open so assume it's fine
    if ((*file).find("dcache")!=string::npos) {//check for null position
      MAXMSG("NuBase",Msg::kInfo,3)
        <<"Filename contains \"dcache\" in path"
        <<" so not testing its existance"<<endl;
      openOk=true;
    }
    else {
      ifstream fileOpenOk((*file).c_str()); 
      openOk=fileOpenOk;
    }

    //check if a wildcard was used because ifstream can't open wildcards
    Int_t stars=(*file).find("*");
    Int_t quest=(*file).find("?");

    //check if file existed
    if (!openOk && !(quest>=0 || stars>=0)){
      MSG("NuBase",Msg::kInfo)
        <<endl<<endl
        <<"***********************************************************"
        <<endl<<"Can't find file="<<*file<<endl
        <<"Note: you can't use '~/'. It has to be the full path"<<endl
        <<"***********************************************************"
        <<endl<<endl
        <<"Exiting here!"<<endl;
      exit(0);
    }
    
    MSG("NuBase",Msg::kInfo)<<"Adding file(s)="<<*file<<endl;
    nf+=fChain->Add((*file).c_str());
    if (fChainBD) nfBD+=fChainBD->Add((*file).c_str());
    //fChain->Print();//very verbose
  }

  if(nf==0){
    MSG("NuBase",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"No *.root files found"<<endl
      <<"Please set NUDATA to the directory containing the"
      <<" *.root files"<<endl
      <<"Or give the txt file containing the files to be input"<<endl
      <<"Note: If more than one file is found they will be"
      <<" concatenated in a TChain and treated as one"<<endl
      <<"*************************************************************"
      <<endl<<endl<<"Program will exit here"<<endl;
    exit(0);
  }
  
  MSG("NuBase",Msg::kInfo)
    <<endl<<"Added "<<nf<<" file(s) to TChain"<<endl;
  MSG("NuBase",Msg::kInfo)
    <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain"
    <<endl;
  
  if (fChain) {
    MSG("NuBase",Msg::kInfo) 
      <<"Ntuple information:"<<endl;
    fChain->Show(0);
    
    if (fChainBD) {
      MSG("NuBase",Msg::kInfo) 
        <<"NtpBDLite information:"<<endl;
      fChainBD->Show(0);
    }
  }
  
  MSG("NuBase",Msg::kInfo)
    <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain."
    <<endl;
  MSG("NuBase",Msg::kInfo)
    <<endl<<"Analysing "<<nf<<" file(s). Reading from disk..."<<endl;

  //set the class data member
  fNumFiles=nf;
}

//......................................................................

void NuBase::SetRootFileObjects()
{
  MSG("NuBase",Msg::kInfo)
    <<"Running the SetRootFileObjects method..."<<endl;

  //set up tree
  fChain->SetBranchAddress("NtpStRecord",&fRec);
  if (fChainBD) fChainBD->SetBranchAddress("NtpBDLiteRecord",&fRecBD);
  //fBDtortgt
    
  //set the number of events in the tree
  fEntries=static_cast<Int_t>(fChain->GetEntries());
  if (fChainBD) fEntriesBD=static_cast<Int_t>(fChainBD->GetEntries());
  MSG("NuBase",Msg::kInfo)
    <<"NtpSt tree has "<<fEntries<<" entries"<<endl;
  MSG("NuBase",Msg::kInfo)
    <<"NtpBDLite tree has "<<fEntriesBD<<" entries"<<endl;
    
  if (fEntries>0){
    //get the snarl/event
    fChain->GetEvent(0); 
    //MSG("NuBase",Msg::kInfo)
    //<<"First run number="<<fRec->fHeader.fRun<<endl;
    MSG("NuBase",Msg::kInfo)
      <<"SetRootFileObjects: first snarl nslice="
      <<fRec->evthdr.nslice<<endl;
  }

  if (fEntriesBD>0){
    //get the snarl/event
    fChainBD->GetEvent(0); 
    MSG("NuBase",Msg::kInfo)
      <<"First tortgt="<<fRecBD->tortgt<<" E12 POT"<<endl;
  }
  MSG("NuBase",Msg::kDebug)
    <<"Finished the SetRootFileObjects method"<<endl;
}

//......................................................................

TFile* NuBase::OpenFile(const NuConfig& config,std::string prefix) const
{
  //call this static function to ensure any histograms created
  //after the file is opened are added
  TH1::AddDirectory(true);

  //create the tfile pointer to be returned
  TFile* outputFile=0;
  
  //get the environmental variable
  char *anaDir=getenv("NUANA_DIR");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("NuBase",Msg::kInfo) 
      <<"Environmental variable $NUANA_DIR not set." 
      <<" Writing file(s) to current directory"<<endl;
    sAnaDir=".";
  }
  
  //convert variables to string
  string sRunNumber=Form("%d",config.run);
  string sZeros="";
  if (config.run>=0 && config.run<10) sZeros="00000000";
  else if (config.run>=10 && config.run<100) sZeros="000000";
  else if (config.run>=100 && config.run<1000) sZeros="00000";
  else if (config.run>=1000 && config.run<10000) sZeros="0000";
  else if (config.run>=10000 && config.run<100000) sZeros="000";
  else if (config.run>=100000 && config.run<1000000) sZeros="00";
  else if (config.run>=1000000 && config.run<10000000) sZeros="0";
  else if (config.run>=10000000 && config.run<100000000) sZeros="";
  sRunNumber=sZeros+sRunNumber;

  string sDetector="UnknownDet";
  if (config.detector==Detector::kNear) {
    sDetector="N";
    if (config.simFlag==SimFlag::kMC) sDetector="n";
  }
  else if (config.detector==Detector::kFar) {
    sDetector="F";
    if (config.simFlag==SimFlag::kMC) sDetector="f";
  }
  else if (config.detector==Detector::kCalDet) {
    sDetector="C";
    if (config.simFlag==SimFlag::kMC) sDetector="c";
  }
  else cout<<"Ahhh, don't know detector="<<config.detector<<endl;
  
  string sPrefix="";//default
  if (prefix!="") sPrefix+=prefix;
  string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
  string sFileName=sBase+".root";
  
  //test if file already exists
  ifstream Test(sFileName.c_str());
  
  //open the appropriate file
  if(!Test){
    outputFile=new TFile(sFileName.c_str(),"RECREATE");
  }
  else {
    //Need new filename
    Int_t fred=1;
    while(Test) {
      Test.close();
      string sAppendage=Form("%d",fred);
      sFileName=sBase+"_"+sAppendage+".root";
      Test.open(sFileName.c_str());
      fred++;
    }
    outputFile=new TFile(sFileName.c_str(),"NEW");
    outputFile->SetCompressionLevel(9);
  }
  
  string sTmp="No File!";
  if (outputFile) sTmp=outputFile->GetName();

  //this is a simple histogram for passing out an error code
  //to the output file if it occured
  //it is meant to be used only to flag errors and to make the
  //user go and look in the log file
  TH1F* hError=new TH1F("hError","hError",110,-10,100);
  hError->GetXaxis()->SetTitle("Error code");
  hError->GetXaxis()->CenterTitle();

  MSG("NuBase",Msg::kInfo) 
    <<"Output file opened: "<<sTmp<<endl;
  return outputFile;
}

//......................................................................

TFile* NuBase::OpenFile(Int_t runNumber,std::string prefix) const
{
  //call this static function to ensure any histograms created
  //after the file is opened are added
  TH1::AddDirectory(true);

  //create the tfile pointer to be returned
  TFile* outputFile=0;
  
  //get the environmental variable
  char *anaDir=getenv("NUANA_DIR");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("NuBase",Msg::kInfo) 
      <<"Environmental variable $NUANA_DIR not set." 
      <<" Writing file(s) to current directory"<<endl;
    sAnaDir=".";
  }
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  string sZeros="";
  if (runNumber>=0 && runNumber<10) sZeros="00000000";
  else if (runNumber>=10 && runNumber<100) sZeros="000000";
  else if (runNumber>=100 && runNumber<1000) sZeros="00000";
  else if (runNumber>=1000 && runNumber<10000) sZeros="0000";
  else if (runNumber>=10000 && runNumber<100000) sZeros="000";
  else if (runNumber>=100000 && runNumber<1000000) sZeros="00";
  else if (runNumber>=1000000 && runNumber<10000000) sZeros="0";
  else if (runNumber>=10000000 && runNumber<100000000) sZeros="";
  sRunNumber=sZeros+sRunNumber;

  //string sDetectorType="C";
  string sDetector="";
  //string sPrefix="h";//default
  string sPrefix="";//default
  if (prefix!="") sPrefix+=prefix;
  string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
  string sFileName=sBase+".root";
  
  //test if file already exists
  ifstream Test(sFileName.c_str());
  
  //open the appropriate file
  if(!Test){
    outputFile=new TFile(sFileName.c_str(),"RECREATE");
  }
  else {
    //Need new filename
    Int_t fred=1;
    while(Test) {
      Test.close();
      string sAppendage=Form("%d",fred);
      sFileName=sBase+"_"+sAppendage+".root";
      Test.open(sFileName.c_str());
      fred++;
    }
    outputFile=new TFile(sFileName.c_str(),"NEW");
    outputFile->SetCompressionLevel(9);
  }
  
  string sTmp="No File!";
  if (outputFile) sTmp=outputFile->GetName();

  MSG("NuBase",Msg::kInfo) 
    <<"Output file opened: "<<sTmp<<endl;
  return outputFile;
}

//......................................................................

TFile* NuBase::OpenFileRECREATE(std::string sFileName) const
{
  //call this static function to ensure any histograms created
  //after the file is opened are added
  TH1::AddDirectory(true);
  
  //create the tfile pointer to be returned
  TFile* outputFile=0;
  
  outputFile=new TFile(sFileName.c_str(),"RECREATE");
  
  MSG("NuBase",Msg::kInfo) 
    <<"Output file opened: "<<outputFile->GetName()<<endl;
  return outputFile;
}

//......................................................................

ofstream* NuBase::OpenTxtFile(const NuConfig& config,
                              std::string prefix) const
{
  //create the tfile pointer to be returned
  ofstream* outputFile=0;
  
  //convert variables to string
  string sRunNumber=Form("%d",config.run);
  string sZeros="";
  if (config.run>=0 && config.run<10) sZeros="00000000";
  else if (config.run>=10 && config.run<100) sZeros="000000";
  else if (config.run>=100 && config.run<1000) sZeros="00000";
  else if (config.run>=1000 && config.run<10000) sZeros="0000";
  else if (config.run>=10000 && config.run<100000) sZeros="000";
  else if (config.run>=100000 && config.run<1000000) sZeros="00";
  else if (config.run>=1000000 && config.run<10000000) sZeros="0";
  else if (config.run>=10000000 && config.run<100000000) sZeros="";
  sRunNumber=sZeros+sRunNumber;

  string sDetector="UnknownDet";
  if (config.detector==Detector::kNear) {
    sDetector="N";
    if (config.simFlag==SimFlag::kMC) sDetector="n";
  }
  else if (config.detector==Detector::kFar) {
    sDetector="F";
    if (config.simFlag==SimFlag::kMC) sDetector="f";
  }
  else if (config.detector==Detector::kCalDet) {
    sDetector="C";
    if (config.simFlag==SimFlag::kMC) sDetector="c";
  }
  else cout<<"Ahhh, don't know detector="<<config.detector<<endl;
  
  string sPrefix="";//default
  string sAnaDir=".";
  if (prefix!="") sPrefix+=prefix;
  string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
  string sFileName=sBase+".txt";

  outputFile=new ofstream(sFileName.c_str());
  
  if (outputFile){
    MSG("NuBase",Msg::kInfo) 
      <<"Output txt file opened: "<<sFileName<<endl;
  }
  else{
    MSG("NuBase",Msg::kInfo)
      <<"Txt file NOT opened: "<<sFileName<<endl;
  }

  return outputFile;
}

//......................................................................

ofstream* NuBase::OpenTxtFile(Int_t runNumber,
                              std::string prefix) const
{
  //create the tfile pointer to be returned
  ofstream* outputFile=0;
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  string sZeros="";
  if (runNumber>=0 && runNumber<10) sZeros="00000000";
  else if (runNumber>=10 && runNumber<100) sZeros="000000";
  else if (runNumber>=100 && runNumber<1000) sZeros="00000";
  else if (runNumber>=1000 && runNumber<10000) sZeros="0000";
  else if (runNumber>=10000 && runNumber<100000) sZeros="000";
  else if (runNumber>=100000 && runNumber<1000000) sZeros="00";
  else if (runNumber>=1000000 && runNumber<10000000) sZeros="0";
  else if (runNumber>=10000000 && runNumber<100000000) sZeros="";
  sRunNumber=sZeros+sRunNumber;

  string sDetectorType="Run";//was F
  string sPrefix="";//default
  string sAnaDir=".";
  if (prefix!="") sPrefix+=prefix;
  string sBase=sAnaDir+"/"+sPrefix+sDetectorType+sRunNumber;
  string sFileName=sBase+".txt";

  outputFile=new ofstream(sFileName.c_str());
  
  if (outputFile){
    MSG("NuBase",Msg::kInfo) 
      <<"Output txt file opened: "<<sFileName<<endl;
  }
  else{
    MSG("NuBase",Msg::kInfo)
      <<"Txt file NOT opened: "<<sFileName<<endl;
  }

  return outputFile;
}

//......................................................................

void NuBase::SetLoopVariables(Int_t entry,Int_t printMode)
{
  this->PrintLoopProgress(entry,fEntries,printMode);
  
  //get the snarl/entry
  fChain->GetEvent(entry);
  if (fChainBD) fChainBD->GetEvent(entry); 
}

//......................................................................

void NuBase::PrintLoopProgress(Int_t entry,Float_t nEntries,
                               Int_t printMode) const
{
  if (printMode==1){
    Float_t fract=ceil(nEntries/20.);  
    if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
      MSG("STAnalysis",Msg::kInfo) 
        <<"Fraction of loop complete: "<<entry 
        <<"/"<<nEntries<<"  ("
        <<(Int_t)(100.*entry/nEntries)<<"%)"<<endl;
    }
  }
}

//......................................................................

void NuBase::InitialiseLoopVariables()
{
  MSG("NuBase",Msg::kDebug)<<"Initialising loop variables..."<<endl;

  MSG("NuBase",Msg::kDebug)<<"Initialisation complete"<<endl;
}

//......................................................................

void NuBase::Test()
{
  MSG("NuBase",Msg::kInfo) 
    <<" ** Running Test method... **"<<endl;

  
  //this->InitialiseLoopVariables();  
  
  //  for(Int_t event=0;event<fEntries;event++){
    
  //this->SetLoopVariables(event);

  for (int i=0;i<5;i++) {
    fChain->GetEntry(i);
    NtpStRecord& ntpst=(*fRec);   // Make a reference instead of a pointer
    
    TClonesArray& tcaTk=(*ntpst.trk);
    Int_t numTrks=tcaTk.GetEntries();
   
    MSG("NuBase",Msg::kInfo)
      <<"numTrks="<<numTrks
      <<", ndigit="<<fRec->evthdr.ndigit
      <<", nTrack="<<fRec->evthdr.ntrack<<endl;

    //Loop over tracks
    for (Int_t itrk=0;itrk<numTrks;itrk++){
      const NtpSRTrack* ptrk=
        dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
      const NtpSRTrack& trk=(*ptrk);
      cout<<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds
          <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:"
          <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl;


      TClonesArray& tcaStp=(*ntpst.stp);
      Int_t numStps=tcaStp.GetEntries();
      MSG("NuBase",Msg::kInfo)
        <<"numStps="<<numStps<<endl;
      
      for (Int_t i=0;i<trk.nstrip;i++){
        
        //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
        const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);

        const NtpSRStrip& stp=(*pstp);
        MSG("NuBase",Msg::kInfo)
          <<"Strip="<<stp.strip<<", pl="<<stp.plane
          <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
      }

    }


    TClonesArray& tcaStp=(*ntpst.stp);
    Int_t numStps=tcaStp.GetEntries();
    MSG("NuBase",Msg::kInfo)
      <<endl<<"2nd: numStps="<<numStps<<endl;
    
    for (Int_t i=0;i<numStps;i++){      
      const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
      const NtpSRStrip& stp=(*pstp);
      MSG("NuBase",Msg::kInfo)
        <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl;
    }
    
    //    for (Int_t nsloop=0;nsloop<trk.nstrip;nsloop++) {
    //const NtpSRStrip* pstp=
    //dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[ns]]);
    //const NtpSRStrip& stp=(*pstp);

    //      if (stp.plane) ;   // Avoid compiler message ifndef PRINTLOTS
    //float du = trk.stpu[ns]-uprev;
    //float dv = trk.stpv[ns]-vprev;
    //float dz = trk.stpz[ns]-zprev;
    //cout<<" "<<ns<<" "<<du<<" "<<dv<<" "<<dz
    //    <<" uv:"<<trk.stpu[ns] <<" "<<trk.stpv[ns]
    //    <<" xyz:"<<trk.stpx[ns] <<" "<<trk.stpy[ns]<<" "<<trk.stpz[ns]
    //    <<" ds:"<<trk.stpds[ns]<<" stp: "<<trk.stp[ns]
    //    <<" "<<stp.plane<<endl;
    //  printf("%2d%8.4f%8.4f%8.4f%8.4f%2d uv:%8.4f%8.4f"
    //       " xyz:%8.4f%8.4f%8.4f%8.4f"
    //       " stp: %4d%4d\n"
    //       ,ns,du,dv,dz,stp.tpos,stp.planeview
    //       ,trk.stpu[ns],trk.stpv[ns]
    //       ,trk.stpx[ns],trk.stpy[ns],trk.stpz[ns]
    //       ,trk.stpds[ns],trk.stp[ns],stp.plane);
      
    //uprev = trk.stpu[ns];
    //vprev = trk.stpv[ns];
    //zprev = trk.stpz[ns];
    //}// stp loop

  }//end of for                                       
  

  MSG("NuBase",Msg::kInfo)<<"Finished main loop"<<endl;

  MSG("NuBase",Msg::kInfo) 
    <<" ** Finished Test method **"<<endl;
}

//......................................................................

void NuBase::TestEventLoop()
{
  MSG("NuBase",Msg::kInfo) 
    <<" ** Running TestEventLoop method... **"<<endl;

  
  this->InitialiseLoopVariables();  
  
  for(Int_t entry=0;entry<fEntries;entry++) {
  //for(Int_t entry=0;entry<5000;entry++) {
    
    this->SetLoopVariables(entry);

    NtpStRecord& ntp=(*fRec);

    TClonesArray& evtTca=(*ntp.evt);
    const Int_t numEvts=evtTca.GetEntriesFast();
    
    TClonesArray& trkTca=(*ntp.trk);
    const Int_t numTrks=trkTca.GetEntriesFast();

    TClonesArray& shwTca=(*ntp.shw);
    const Int_t numShws=shwTca.GetEntriesFast();

    TClonesArray& stpTca=(*ntp.stp);
    const Int_t numStps=stpTca.GetEntries();

    TClonesArray& slcTca=(*ntp.slc);
    const Int_t numSlcs=slcTca.GetEntries();
    
    Float_t totalSnarlAdc=0;
    Float_t totalEvtAdc=0;
    Float_t totalSlcAdc=0;
    Float_t thisSlcAdc=0;
    
    for (Int_t i=0;i<numStps;i++) {
      const NtpSRStrip* pstp=
        dynamic_cast<NtpSRStrip*>(stpTca[i]);
      const NtpSRStrip& stp=(*pstp);
      totalSnarlAdc+=stp.ph0.raw+stp.ph1.raw;
    }
    
    if (numEvts==0 && numTrks>0) {
      cout<<"Ahhh, num events in tca="<<numEvts
          <<", trks="<<numTrks<<endl;
    }
    
    if (numEvts>1) {
      MAXMSG("NuBase",Msg::kInfo,10)
        <<"High num evts="<<numEvts<<", trks="<<numTrks<<endl;
    }

    if (numSlcs!=1) {
      MAXMSG("NuBase",Msg::kInfo,10)
        <<"High num slices in tca="<<numSlcs
        <<", trks="<<numTrks<<endl;
      //continue;
    }
    
    MAXMSG("NuBase",Msg::kInfo,100)
      <<"Num events in tca="<<numEvts
      <<", slcs="<<numSlcs<<", trks="<<numTrks<<endl;
    
    Int_t slcStps=0;
    
    //loop over the slices
    for (Int_t islc=0;islc<numSlcs;++islc) {
      const NtpSRSlice* pslc=
        dynamic_cast<NtpSRSlice*>(slcTca[islc]);
      const NtpSRSlice& slc=(*pslc);

      MAXMSG("NuBase",Msg::kDebug,200)
        <<"  slice "<<islc+1<<" of "<<numSlcs<<endl;
      
      slcStps+=slc.nstrip;

      //loop over strips in slc
      for (Int_t i=0;i<slc.nstrip;++i) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        totalSlcAdc+=stp.ph0.raw+stp.ph1.raw;
      }

      //loop over tracks in slc  ---- can't do this!!!
      //for(int i=0;i<slc.ntrack;i++){    
      //const NtpSRTrack* ptrk=
      //  dynamic_cast<NtpSRTrack*>(trkTca[slc.trk[i]]);
      //const NtpSRTrack& trk=(*ptrk);

      //trkStps=trk.nstrip;
      //MAXMSG("NuBase",Msg::kInfo,300)
      //  <<"    track "<<i+1<<" of "<<slc.ntrack
      //  <<", trkStps="<<trkStps<<endl;
      //}
    }

    map<Int_t,Int_t> evtPerSlc;
    
    //loop over the events
    for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
      const NtpSREvent* pevt=
        dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
      const NtpSREvent& evt=(*pevt);
      
      MAXMSG("NuBase",Msg::kInfo,100)
        <<"Entry "<<entry<<" has tracks="<<evt.ntrack
        <<", shws="<<numShws<<", slcs="<<numSlcs<<endl;

      //get the slice associated with this event
      const NtpSRSlice* pslc=
        dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
      const NtpSRSlice& slc=(*pslc);

      evtPerSlc[evt.slc]++;

      //check that the filter works
      Int_t evts=-1;
      Bool_t goodEvt=true;
      if (goodEvt){
        MAXMSG("NuBase",Msg::kInfo,100)
          <<"TEST::entry="<<entry
          <<", evt="<<ntpevt<<" is good, evtsPerSlc="<<evts
          <<" (slc="<<evt.slc<<")"<<endl;
      }
      else{
        MAXMSG("NuBase",Msg::kInfo,100)
          <<"TEST::entry="<<entry
          <<", evt="<<ntpevt<<" is bad, evtsPerSlc="<<evts
          <<" (slc="<<evt.slc<<")"<<endl;
      }
      
      if (evt.ntrack!=1 || evt.nshower>1) continue;

      Int_t trkStps=0;
      Int_t shwStps=0;
      
      //loop over strips in slc
      for (Int_t i=0;i<slc.nstrip;++i) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        thisSlcAdc+=stp.ph0.raw+stp.ph1.raw;
      }
      
      //loop over tracks in evt
      for(int i=0;i<evt.ntrack;i++) {    
        const NtpSRTrack* ptrk=
          dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[i]]);
        const NtpSRTrack& trk=(*ptrk);

        trkStps=trk.nstrip;
        
        for (Int_t i=0;i<trk.nstrip;i++) {
          const NtpSRStrip* pstp=
            dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
          const NtpSRStrip& stp=(*pstp);
            MAXMSG("NuBase",Msg::kDebug,100)
              <<"strip time="<<stp.time0<<endl;
        }
      }

      //loop over showers in evt
      for(int i=0;i<evt.nshower;i++) {    
        const NtpSRShower* pshw=
          dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
        const NtpSRShower& shw=(*pshw);
        
        shwStps=shw.nstrip;
      }

      //loop over strips in evt
      for (Int_t i=0;i<evt.nstrip;i++) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        totalEvtAdc+=stp.ph0.raw+stp.ph1.raw;
      }
      MAXMSG("NuBase",Msg::kInfo,100)
        <<"e="<<entry<<", evt="<<ntpevt
        <<", ADCs: snl="<<totalSnarlAdc
        <<", evt="<<totalEvtAdc
        <<", slc="<<totalSlcAdc
        <<", tslc="<<thisSlcAdc
        <<", #st="<<numStps<<","<<evt.nstrip<<","<<slcStps
        //<<", trk="<<trkStps<<", shw="<<shwStps
        <<endl;
    }

    //make this into a function
    //Bool_t FilterMultiEvtPerSlc(ntp,slc,entry)
    typedef map<Int_t,Int_t>::const_iterator slcIt;
    MAXMSG("NuBase",Msg::kInfo,30)
      <<"MAIN::Events per slice for entry="<<entry<<endl;
    for (slcIt it=evtPerSlc.begin();it!=evtPerSlc.end();++it) {
      MAXMSG("NuBase",Msg::kInfo,200)
        <<"  slice "<<it->first<<" has "<<it->second<<" event(s)"<<endl;
    }
    
  }//end of for                                       
    
  
  MSG("NuBase",Msg::kInfo)<<"Finished main loop"<<endl;

  MSG("NuBase",Msg::kInfo) 
    <<" ** Finished TestEventLoop method **"<<endl;
}

//......................................................................

---