GnumiInterface.cxx

Go to the documentation of this file.

#ifndef gnumiinterface_cxx
#define gnumiinterface_cxx
#include "StandardNtuple/NtpStRecord.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "MCReweight/GnumiInterface.h"
#include "MCReweight/NuParent.h"
#include "TClonesArray.h"
#include <iostream>

using namespace std;

//****************************
GnumiInterface::GnumiInterface()
{
  fok = false;
  char *dir = getenv("GNUMIAUX");
  if(dir!=NULL) fdir = dir;
  else fdir = ".";
  if(OpenFiles()) SetAddresses();
}

//****************************
GnumiInterface::~GnumiInterface()
{
  CloseFiles();
}

//****************************
Bool_t GnumiInterface::SetFileDir(std::string dir)
{
  fdir = dir;
  if(OpenFiles()) {
    SetAddresses();
    return true;
  }
  return false;
} 

//****************************
Bool_t GnumiInterface::OpenFiles()
{
  map<Int_t,TFile*>::iterator End = fFile.end();
  for(int i=1;i<=20;i++){
    map<Int_t,TFile*>::iterator Cur = fFile.find(i);
    if(Cur!=End){
      TFile *f1 = Cur->second;
      f1->Close();     //close file 
      delete f1;       //delete TFile
      fFile.erase(Cur->first);  //erase file map entry
      fTree.erase(Cur->first);  //erase tree map entry
    }
    char name[256];
    sprintf(name,"%s/job%i/le_%i.root",fdir.c_str(),i,i);
    cout << name << endl;
    TFile *f = new TFile(name,"READ");
    if(f->IsOpen() && !f->IsZombie()){
      fFile[i] = f;
      TTree *tree = (TTree*) f->Get("h10");
      fTree[i] = tree;
    }
    else if(i==1) return false;
    else break;
  }
  fok = true;
  return true;
}

//****************************
void GnumiInterface::SetAddresses()
{
  if(!fok) return;
  map<Int_t,TTree*>::iterator Beg = fTree.begin();
  map<Int_t,TTree*>::iterator End = fTree.end();
  while(Beg!=End){
    TTree *tree = Beg->second;
    tree->SetBranchAddress("run",&fRunNum);
    tree->SetBranchAddress("evtno",&fEvtNum);
    tree->SetBranchAddress("tvx",&fX);
    tree->SetBranchAddress("tvy",&fY);
    tree->SetBranchAddress("tvz",&fZ);
    tree->SetBranchAddress("tpx",&fPX);
    tree->SetBranchAddress("tpy",&fPY);
    tree->SetBranchAddress("tpz",&fPZ);
    tree->SetBranchAddress("tptype",&fPID);
    tree->SetBranchAddress("tgen",&fGen);
    Beg++;
  }

}

//****************************
void GnumiInterface::CloseFiles()
{
  fTree.clear();
  map<Int_t,TFile*>::iterator Beg = fFile.begin();
  map<Int_t,TFile*>::iterator End = fFile.end();
  while(Beg!=End){
    TFile *f = Beg->second;
    f->Close();
    delete f;
    fFile.erase(Beg->first);
    Beg++;
  }
  fok = false;
}

//****************************
void GnumiInterface::GetParent(Int_t fileNum,Int_t entryNum,
                               Int_t runNum,Int_t evtNum,
                               NuParent &par)
{
  par.Zero();
  //NOTES:
  //1) fileNum is the order that a shell would list the files:
  //i.e. 1,10,11,12,..19,2,20,3,4...
  //this code puts the files into a map using the
  //filename# (le_#.root) as the key - this is the same as 
  //the runNum so set fileNum = runNum to get the appropriate file
  fileNum = runNum;

  //2) There are entries in the stdhep array with 
  //runNum = 99, corresponding to flux file job9/le_9.root 
  //if fileNum is 99, change it to 9 in order to pick up the right file
  //(runNum is indeed 99 in le_9.root, so don't change that!) 
  if(fileNum==99) fileNum = 9;

  //3) entryNum is in fortran numbering i.e. 1,2,3....
  //subtract 1 to get to C numbering:
  entryNum -= 1;

  map<Int_t,TTree*>::iterator End = fTree.end();  
  map<Int_t,TTree*>::iterator Cur = fTree.find(fileNum);
  if(Cur==End) return;
  Cur->second->GetEntry(entryNum);
  if(fRunNum==runNum && fEvtNum==evtNum) {
    par.SetX(fX);
    par.SetY(fY);
    par.SetZ(fZ);
    par.SetPx(fPX);
    par.SetPy(fPY);
    par.SetPz(fPZ);
    par.SetPID(fPID);
    par.SetGen(fGen);
  }
  return;
}

//****************************
void GnumiInterface::GetParent(NtpMCRecord *rec,int mcevent,
                               NuParent &par)
{

  Int_t filenum = -1;
  Int_t entrynum = -1;
  Int_t runnum = -1;
  Int_t evtnum = -1;

  TClonesArray& heparray = *(rec->stdhep);
  Int_t nhep = heparray.GetEntries();  
  for(int i=0;i<nhep;i++){
    NtpMCStdHep *sh = static_cast<NtpMCStdHep *>(heparray[i]);
    if(sh->mc==mcevent){
      if(sh->IstHEP==999){ //beam file info
        runnum = Int_t(sh->p4[0]);
        evtnum = Int_t(sh->p4[1]);
        filenum = Int_t(sh->p4[2]);
        entrynum = Int_t(sh->p4[3]);
      }
    }
  }
  this->GetParent(filenum,entrynum,runnum,evtnum,par);
}

//****************************
void GnumiInterface::GetParent(NtpStRecord *rec,int mcevent,
                               NuParent &par)
{

  Int_t filenum = -1;
  Int_t entrynum = -1;
  Int_t runnum = -1;
  Int_t evtnum = -1;

  TClonesArray& heparray = *(rec->stdhep);
  Int_t nhep = heparray.GetEntries();  
  for(int i=0;i<nhep;i++){
    NtpMCStdHep *sh = static_cast<NtpMCStdHep *>(heparray[i]);
    if(sh->mc==mcevent){
      if(sh->IstHEP==999){ //beam file info
        runnum = Int_t(sh->p4[0]);
        evtnum = Int_t(sh->p4[1]);
        filenum = Int_t(sh->p4[2]);
        entrynum = Int_t(sh->p4[3]);
      }
    }
  }
  this->GetParent(filenum,entrynum,runnum,evtnum,par);
}

//****************************
void GnumiInterface::FillANtpTruth(Int_t fileNum,
                                   Int_t entryNum,
                                   Int_t runNum,
                                   Int_t evtNum,
                                   ANtpTruthInfoBeam &ntpBeam)
{
  NuParent parent;
  this->GetParent(fileNum,entryNum,runNum,evtNum,parent);
  ntpBeam.parentX   = parent.GetX();
  ntpBeam.parentY   = parent.GetY();
  ntpBeam.parentZ   = parent.GetZ();
  ntpBeam.parentPX  = parent.GetPx();
  ntpBeam.parentPY  = parent.GetPy();
  ntpBeam.parentPZ  = parent.GetPz();
  ntpBeam.parentPID = parent.GetPID();
  ntpBeam.parentGen = parent.GetGen();
}

//****************************
void GnumiInterface::FillANtpTruth(NtpMCRecord *rec,int mcevent,
                                   ANtpTruthInfoBeam &ntpBeam)
{
  NuParent parent;
  this->GetParent(rec,mcevent,parent);
  ntpBeam.parentX   = parent.GetX();
  ntpBeam.parentY   = parent.GetY();
  ntpBeam.parentZ   = parent.GetZ();
  ntpBeam.parentPX  = parent.GetPx();
  ntpBeam.parentPY  = parent.GetPy();
  ntpBeam.parentPZ  = parent.GetPz();
  ntpBeam.parentPID = parent.GetPID();
  ntpBeam.parentGen = parent.GetGen();
}

//****************************
void GnumiInterface::FillANtpTruth(NtpStRecord *rec,int mcevent,
                                   ANtpTruthInfoBeam &ntpBeam)
{
  NuParent parent;
  this->GetParent(rec,mcevent,parent);
  ntpBeam.parentX   = parent.GetX();
  ntpBeam.parentY   = parent.GetY();
  ntpBeam.parentZ   = parent.GetZ();
  ntpBeam.parentPX  = parent.GetPx();
  ntpBeam.parentPY  = parent.GetPy();
  ntpBeam.parentPZ  = parent.GetPz();
  ntpBeam.parentPID = parent.GetPID();
  ntpBeam.parentGen = parent.GetGen();
}

#endif

---