MoqBase.cxx

Go to the documentation of this file.

#ifndef moqbase_cxx
#define moqbase_cxx
#include <iostream>
//root includes
#include "TClonesArray.h"
#include "OfflineMonitor/MoqBase.h"

//ClassImp(MoqBase)

MoqBase::MoqBase(TChain *chainSR,TChain *chainMC,
                 TChain *chainTH,TChain *chainEM)
{

  fCurRun = -1;
  lastEntry = -1;
  if(!chainSR) {
    record = 0;
    emrecord = 0;
    mcrecord = 0;
    threcord = 0;
    Clear();
    whichSource = -1;
    isMC = true;
    isTH = true;
    isEM = true;
    Nentries = -1;
    fCurRun = -1;
    return;
  }

  InitChain(chainSR,chainMC,chainTH,chainEM);
  whichSource = 0;

}

MoqBase::MoqBase(JobC *j,string path,int entries)
{
  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  Clear();
  isMC = true;
  isTH = true;
  isEM = true;
  Nentries = entries;
  whichSource = 1;
  jcPath = path;
  fCurRun = -1;
  lastEntry = -1;
  fJC = j;
  fChain = NULL;
}

MoqBase::~MoqBase()
{
  if(fChain) delete fChain;
}


Int_t MoqBase::GetEntry(Int_t entry)
{
  //zero everything before next snarl
  Clear();
  
  //get records
  int status = 0;
  if(whichSource==-1) return 0;
  else if(whichSource==0) {
    if(entry!=lastEntry) {
      status = fChain->Get(entry);
      lastEntry = entry;
    }
    else status = 1;
    record = fChain->Record;
    mcrecord = fChain->mcRecord;
    threcord = fChain->thRecord;
    emrecord = fChain->emRecord;
  }
  else if(whichSource==1) {
    //allow for possibility that reco path can change event by event:
    isMC = isTH = isEM = true; 
    //Get entry:
    //cout << "Getting Cand entry " << entry << " from " 
    // << jcPath.c_str() << endl;
    if(entry!=lastEntry) {
      if(fCurRun!=-1) fJC->Input.GoTo(fCurRun,entry);
      fJC->Path(jcPath.c_str()).RunNin(1);
      lastEntry = entry;
    }
    //Look for records in mom:
    record   = (NtpSRRecord*) fJC->Mom.GetFragment("NtpSRRecord");
    mcrecord = (NtpMCRecord*) fJC->Mom.GetFragment("NtpMCRecord");
    threcord = (NtpTHRecord*) fJC->Mom.GetFragment("NtpTHRecord");
    emrecord = (NtpEMRecord*) fJC->Mom.GetFragment("NtpEMRecord");
    //if no SR record quit
    if(record==0) return 0;
    //otherwise check for other records:
    status = 1;
    if(mcrecord==0) isMC = false;
    if(threcord==0) isTH = false;
    if(emrecord==0) isEM = false;
  }
  else return 0;
  
  //set other things:
  ntpHeader = &(record->GetHeader());
  fCurRun = ntpHeader->GetRun();
  shieldSummary = &(record->vetohdr);
  eventSummary = &(record->evthdr);
  LoadSlice(0);  
  LoadEvent(0);
  LoadTrack(0);
  LoadShower(0);
  LoadShieldStrip(0);

  if(isEM) {
    ntpEMSummary = &(emrecord->emhdr);
    LoadEMShower(0);
  }

  if(isMC) {
    LoadTruth(0);
    LoadStdHep(0);
  }

  if(isTH) {
    LoadTHSlice(0);
    LoadTHTrack(0);
    LoadTHShower(0);
  }
  return status;

}

void MoqBase::InitChain(TChain *chainSR,TChain *chainMC,
                        TChain *chainTH,TChain *chainEM)
{

  record = 0;
  emrecord = 0;
  mcrecord = 0;
  threcord = 0;
  
  fChain = new MoqChain(chainSR,chainMC,chainTH,chainEM);
  isMC = fChain->isMC;
  isTH = fChain->isTH;
  isEM = fChain->isEM;
  Nentries = fChain->Nentries;  

  if(isMC) std::cout << "MC Present" << std::endl;
  if(isTH) std::cout << "TH Present" << std::endl;
  if(isEM) std::cout << "EM Present" << std::endl;
  
}

void MoqBase::Clear(){

  ntpHeader = 0;
  shieldSummary = 0;
  eventSummary = 0;
  ntpSlice=0;
  ntpEvent=0;
  ntpTrack=0;
  ntpShower=0;
  ntpStrip=0;
  ntpShieldStrip=0;

  ntpTruth=0;
  ntpStdHep=0;

  ntpTHSlice=0;
  ntpTHTrack=0;
  ntpTHShower=0;

  ntpEMSummary = 0;
  ntpEMShower = 0;

}


Bool_t MoqBase::LoadTrack(Int_t itrk){

  Bool_t status = false;
  TClonesArray& trackArray = *(record->trk);
  if(itrk>=0&&itrk<trackArray.GetEntries()) {
    ntpTrack = dynamic_cast<NtpSRTrack *>(trackArray[itrk]);
    status = true;
  }
  else ntpTrack = 0;
  
  return status;
}

Bool_t MoqBase::LoadTHTrack(Int_t itrk){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray& trackArray = *(threcord->thtrk);
  if(itrk>=0&&itrk<trackArray.GetEntries()) {
    ntpTHTrack = dynamic_cast<NtpTHTrack *>(trackArray[itrk]);
    status = true;
  }
  else ntpTHTrack = 0;
  
  return status;
}

Bool_t MoqBase::LoadShower(Int_t ishw){

  Bool_t status = false;
  TClonesArray& showerArray = *(record->shw);
  if(ishw>=0&&ishw<showerArray.GetEntries()) {
    ntpShower = dynamic_cast<NtpSRShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpShower = 0;
  
  return status;
}

Bool_t MoqBase::LoadTHShower(Int_t ishw){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray& showerArray = *(threcord->thshw);
  if(ishw>=0&&ishw<showerArray.GetEntries()) {
    ntpTHShower = dynamic_cast<NtpTHShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpTHShower = 0;
  
  return status;
}

Bool_t MoqBase::LoadSlice(Int_t islc){

  Bool_t status = false;
  TClonesArray& sliceArray = *(record->slc);
  if(islc>=0&&islc<sliceArray.GetEntries()) {
    ntpSlice = dynamic_cast<NtpSRSlice *>(sliceArray[islc]);
    status = true;
  }
  else ntpSlice = 0;
  
  return status;
}

Bool_t MoqBase::LoadTHSlice(Int_t islc){

  Bool_t status = false;
  if(!isTH) return status; 
  TClonesArray& sliceArray = *(threcord->thslc);
  if(islc>=0&&islc<sliceArray.GetEntries()) {
    ntpTHSlice = dynamic_cast<NtpTHSlice *>(sliceArray[islc]);
    status = true;
  }
  else ntpTHSlice = 0;
  
  return status;
}

Bool_t MoqBase::LoadEvent(Int_t ievt){

  Bool_t status = false;
  TClonesArray& eventArray = *(record->evt);
  if(ievt>=0&&ievt<eventArray.GetEntries()) {
    ntpEvent = dynamic_cast<NtpSREvent *>(eventArray[ievt]);
    status = true;
  }
  else ntpEvent = 0;
  
  return status;
}

Bool_t MoqBase::LoadTruth(Int_t itr){
  
  Bool_t status = false;
  if(!isMC) return status;
  TClonesArray& mcArray = *(mcrecord->mc);
  if(itr>=0&&itr<mcArray.GetEntries()) {
    ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[itr]);
    status = true;
  }
  else ntpTruth = 0;
  
  return status;
}

Bool_t MoqBase::LoadStdHep(Int_t ish){
  
  Bool_t status = false;
  if(!isMC) return status;
  TClonesArray& stdhepArray = *(mcrecord->stdhep);
  if(ish>=0&&ish<stdhepArray.GetEntries()) {
    ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[ish]);
    status = true;
  }
  else ntpStdHep = 0;
  
  return status;
}

Bool_t MoqBase::LoadTruthForReco(Int_t ievt,Int_t &itruth){

  Bool_t status = false;  
  ntpTruth = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;

  TClonesArray& mcArray = *(mcrecord->mc);
  if(mcArray.GetEntries()==1) { //if only one MC event is present just load it
    LoadTruth(0);
    return true;
  }

  float mindist = 999;
  int itr = 0;
  //loop over all MC events and look for closest vertex to reco vertex
  for(int i=0;i<mcArray.GetEntries();i++){
    ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[i]);
    float dist = sqrt(TMath::Power(ntpTruth->vtxx/100.-ntpEvent->vtx.x,2)
                      +TMath::Power(ntpTruth->vtxy/100.-ntpEvent->vtx.y,2)
                      +TMath::Power(ntpTruth->vtxz/100.-ntpEvent->vtx.z,2));
    if(dist<mindist) { 
      itr = i;
      mindist = dist;
    }
  }

  if(mindist<0.5) {  //Need to place a fairly loose cut to pick up showers
    status = true;   //From looking at FarDet events, a cut of 0.5m 
    LoadTruth(itr);  //correctly pairs up 95% of events (LE beam,unoscillated) 
    itruth = itr;
  }
  return status;
}

Bool_t MoqBase::LoadSliceForRecoTH(Int_t ievt,Int_t &islice){

  Bool_t status = false;
  ntpSlice = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;
  if(!isTH) return status;

  Int_t itru = 0;
  if(LoadTruthForRecoTH(ievt,itru)){
    Int_t count = 0;
    while(!status){
      if(!LoadTHSlice(count)) break;
      if(ntpTHSlice->neumc==itru) {
        if(LoadSlice(count)){
          islice = count;
          status = true;
        }
      }
      count += 1;
    }
  }
  return status;
}


Bool_t MoqBase::LoadTruthForRecoTH(Int_t ievt,Int_t &itruth){

  Bool_t status = false;
  ntpTruth = 0;
  if(!LoadEvent(ievt)) return status;
  if(!isMC) return status;
  if(!isTH) return status;
  
  if(ntpEvent->ntrack>0) {
    int *tracks = ntpEvent->trk;
    if(LoadTHTrack(int(tracks[0]))){
      int ind = ntpTHTrack->neumc;
      if(LoadTruth(ind)) {
        itruth = ind;
        status = true;
      }
    }
  }
  if(!status&&ntpEvent->nshower>0) {
    int *showers = ntpEvent->shw;    
    if(LoadTHShower(int(showers[0]))){
      int ind = ntpTHShower->neumc;
      if(LoadTruth(ind)) {
        itruth = ind;
        status = true;
      }
    }
  }
  return status;

}

Bool_t MoqBase::LoadStrip(Int_t istp){

  Bool_t status = false;
  TClonesArray& stripArray = *(record->stp);
  if(istp>=0&&istp<stripArray.GetEntries()) {
    ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[istp]);
    status = true;
  }
  else ntpStrip = 0;
  
  return status;
}

Bool_t MoqBase::LoadShieldStrip(Int_t istp){

  Bool_t status = false;
  TClonesArray& stripArray = *(record->vetostp);
  if(istp>=0&&istp<stripArray.GetEntries()) {
    ntpShieldStrip = dynamic_cast<NtpSRShieldStrip *>(stripArray[istp]);
    status = true;
  }
  else ntpShieldStrip = 0;
  
  return status;
}

Bool_t MoqBase::LoadEMShower(Int_t ishw){

  Bool_t status = false;
  if(!isEM) return status; 
  TClonesArray& showerArray = *(emrecord->emshw);
  if(ishw>=0&&ishw<showerArray.GetEntries()) {
    ntpEMShower = dynamic_cast<NtpEMShower *>(showerArray[ishw]);
    status = true;
  }
  else ntpEMShower = 0;
  
  return status;
}

#endif // #ifdef moqbase_cxx

---