MoqQuantities.cxx

Go to the documentation of this file.

#ifndef moqquantities_cxx
#define moqquantities_cxx
#include <iostream>
#include <cmath>
#include "TClonesArray.h"
#include "TMath.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "Swimmer/SwimSwimmer.h"
#include "Swimmer/SwimParticle.h"
#include "Swimmer/SwimZCondition.h"
#include "Conventions/Munits.h"
#include "OfflineMonitor/MoqQuantities.h"

//********************************************************

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

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


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

}

//********************************************************

MoqQuantities::MoqQuantities(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;
  fJC = j;
  fChain = NULL;


}

//********************************************************

MoqQuantities::~MoqQuantities()
{

}

//********************************************************

Bool_t MoqQuantities::IsTrack(){
  if(eventSummary->ntrack>0) return true;
  return false;
}

//********************************************************

Bool_t MoqQuantities::IsSingleTrack(){
  if(eventSummary->ntrack==1) return true;
  return false;
}

//********************************************************

Bool_t MoqQuantities::PassTrackCuts(Int_t itrk){
  if(eventSummary->nevent!=1) return false; //is an event
  if(!IsTrack()) return false; //is a track
  LoadTrack(itrk); //load track
  if(ntpTrack->fit.pass==0) return false; //track fit passes
  if(ntpTrack->momentum.qp!=ntpTrack->momentum.qp) return false; //no nans
  if(fabs(ntpTrack->momentum.qp)<0.0005) return false; //no moms>2000
  return true;
}

//********************************************************

Bool_t MoqQuantities::PassCuts(Int_t itrk){
  //basic fit track cuts:
  if(!PassTrackCuts(itrk)) return false;
  //check for contained vtx:
  //  if(ntpTrack->fidvtx.dr<0.5||ntpTrack->fidvtx.dz<0.5) return false;
  return true;
}


//********************************************************

Bool_t MoqQuantities::IsFidAll(Int_t itrk){
  if(!LoadTrack(itrk)) return false;
  if(ntpTrack->fidall.dr<0.5||ntpTrack->fidall.dz<0.5) return false;
  return true;
}


//********************************************************

Int_t MoqQuantities::INu(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  return ntpTruth->inu;
}

//********************************************************

Int_t MoqQuantities::IAction(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  return ntpTruth->iaction;
}

//********************************************************

Int_t MoqQuantities::IResonance(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  return ntpTruth->iresonance;
}




//********************************************************

Float_t MoqQuantities::Y(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  return ntpTruth->y;
}

//********************************************************

Float_t MoqQuantities::W2(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  return ntpTruth->w2;
}

//********************************************************

Float_t MoqQuantities::TrueNuEnergy(Int_t itr){
  if(!LoadTruth(itr)) return 0.;
  return ntpTruth->p4neu[3];
}

//********************************************************

Float_t *MoqQuantities::TrueVtx(Int_t itr){
  Float_t *vtx = new Float_t[3];  
  vtx[0] = 0.; 
  vtx[1] = 0.; 
  vtx[2] = 0.;
  if(!LoadTruth(itr)) return vtx;
  vtx[0] = ntpTruth->vtxx;
  vtx[1] = ntpTruth->vtxy;  
  vtx[2] = ntpTruth->vtxz;
  return vtx;
}

//********************************************************

Float_t MoqQuantities::TrueMuEnergy(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  if(ntpTruth->p4mu1[3]!=0) return fabs(ntpTruth->p4mu1[3]);
  return 0;
}

//********************************************************

Float_t MoqQuantities::TrueMuQP(Int_t itr){
  if(!LoadTruth(itr)) return 0;
  if(ntpTruth->p4mu1[3]!=0) {
    Float_t p = 1./sqrt(ntpTruth->p4mu1[3]*ntpTruth->p4mu1[3]-0.10555*0.10555);
    Float_t q = ntpTruth->p4mu1[3]/fabs(ntpTruth->p4mu1[3]);
    return q*p;
  }
  return 0;
}

//********************************************************

Float_t MoqQuantities::TrueShwEnergy(Int_t itr){
  if(!LoadTruth(itr)) return 0.;
  return ntpTruth->p4shw[3];
}

//********************************************************

Float_t MoqQuantities::GetTrueShwSC(Int_t itr){

  if(!LoadTruth(itr)) return 0.;
  
  Float_t Summed_Shw_SC = 0;
  
  return Summed_Shw_SC;
}

//********************************************************

Double_t MoqQuantities::TrueVxB(Int_t itr){
  
  if(!LoadTruth(itr)) return -1; //make sure there is an mc event
  if(ntpTruth->p4mu1[3]==0) return -1; //make sure there is a muon
  
  float truvtxu = ntpTruth->vtxx/sqrt(2.) + ntpTruth->vtxy/sqrt(2.);
  float truvtxv = - ntpTruth->vtxx/sqrt(2.) + ntpTruth->vtxy/sqrt(2.);
  float truvtxz = ntpTruth->vtxz;
  TVector3 trupos(truvtxu/100.,truvtxv/100.,truvtxz/100.);
  
  TVector3 trumom(ntpTruth->p4mu1[0]/sqrt(2.)+ntpTruth->p4mu1[1]/sqrt(2.),
                  -ntpTruth->p4mu1[0]/sqrt(2.)+ntpTruth->p4mu1[1]/sqrt(2.),
                  ntpTruth->p4mu1[2]);

  //assume that reco z is the same as true z for now:  
  return SwimToGetVxB(trupos,trumom,
                      int(ntpTruth->p4mu1[3]/fabs(ntpTruth->p4mu1[3])),
                      double(ntpTrack->end.z));
}


//********************************************************

Float_t MoqQuantities::RecoMuEnergy(Int_t opt,Int_t itrk){

  if(LoadTrack(itrk)){
    if(opt==0){  //return the most appropriate measure of momentum
      if(ntpTrack->fidall.dr>0.5&&ntpTrack->fidall.dz>0.5) {
        return sqrt(ntpTrack->momentum.range*ntpTrack->momentum.range
                    + 0.10555*0.10555);
      }
      else {
        return sqrt(1./(ntpTrack->momentum.qp*ntpTrack->momentum.qp)
                    + 0.10555*0.10555);
      }
    }
    else if(opt==1) //return curvature measurement
      return sqrt(1./(ntpTrack->momentum.qp*ntpTrack->momentum.qp)
                  + 0.10555*0.10555);
    else if(opt==2) //return range measurement
      return sqrt(ntpTrack->momentum.range*ntpTrack->momentum.range
                  + 0.10555*0.10555);
    else return 0;
  }
  return 0.;
}

//********************************************************

Float_t MoqQuantities::RecoMuQP(Int_t itrk){
  if(LoadTrack(itrk)) return ntpTrack->momentum.qp;
  return 0.;
}

//********************************************************

Float_t MoqQuantities::GetNonTrkSC(Int_t itrk,Int_t npln){
  
  Float_t Summed_Stp_SC = 0;
  Float_t Summed_Trk_SC = 0;
  
  bool useAll = false;
  if(npln<=0) useAll=true;  //Most basic thing in the world ever: 
                            //if it's not in the track it's in the shower
  
  TClonesArray& stripArray = *(record->stp);
  if(LoadTrack(itrk)) {
    Int_t numtrkstp = ntpTrack->nstrip;
    for(Int_t i=0;i<numtrkstp;i++){
      Int_t index = ntpTrack->stp[i];
      ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
      if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
          && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
        Summed_Trk_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
      }
    }
  }
  
  for(Int_t i=0;i<int(eventSummary->nstrip);i++) {
    if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
        && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
      Summed_Stp_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
    }
  }
  return Summed_Stp_SC - Summed_Trk_SC;
}

//********************************************************

Float_t MoqQuantities::GetShwSC(Int_t itrk){

  //Try something a bit smarter:
  //Look for a reconstructed shower near to the track vtx and sum the SigCors
  if(!LoadTrack(itrk)) return 0.;
  
  Float_t Summed_Shw_SC = 0;
  Int_t bestshower = -1;
  Float_t bestdistance = 1000.;
  if(eventSummary->nshower!=0){
    for(Int_t i=0;i<eventSummary->nshower;i++){
      LoadShower(i);
      Float_t distance = sqrt((ntpTrack->vtx.x - ntpShower->vtx.x)
                              *(ntpTrack->vtx.x - ntpShower->vtx.x)
                              +(ntpTrack->vtx.y - ntpShower->vtx.y)
                              *(ntpTrack->vtx.y - ntpShower->vtx.y)
                              +(ntpTrack->vtx.z - ntpShower->vtx.z)
                              *(ntpTrack->vtx.z - ntpShower->vtx.z));
      
      if(distance<0.5){
        if(bestshower==-1) {
          bestshower = i;
          bestdistance = distance;
        }
        else if(distance<bestdistance) bestshower = i;
      }
    }
    
    TClonesArray& stripArray = *(record->stp);
    if(bestshower!=-1){
      LoadShower(bestshower);
      Int_t *shwstrips = ntpShower->stp;
      for(Int_t j=0;j<ntpShower->nstrip;j++){
        Int_t index = shwstrips[j];
        ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
        Summed_Shw_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
      }
      return Summed_Shw_SC;
    }
    return 0;
  }
  else return 0; //no reconstructed shower

}

//********************************************************

Float_t MoqQuantities::RecoShwEnergy(Int_t opt){
  
  Float_t theGeV = 0;
  //  if(opt==0) theGeV = GetShwSC(0)/10000.0; //MDC
  //else theGeV = GetNonTrkSC()/10000.0;

   // Return Jim's calculation
  if(LoadShower(opt)) theGeV=ntpShower->ph.gev;
  else theGeV=0.;

  return theGeV;
  
}

//********************************************************

Int_t MoqQuantities::GetChargeSign(Int_t itrk){  
  if(RecoMuQP(itrk)>0) return 1;
  return -1;
}

//********************************************************

Float_t MoqQuantities::RecoQENuEnergy(Int_t itrk){
  
  if(!LoadTrack(itrk)) return 0.;
  Float_t nucleonMass = 0.93956563; //mass of neutron by default
  if(GetChargeSign(itrk)==1) nucleonMass = 0.93827231; //proton mass for nubar
  Float_t muonEnergy = RecoMuEnergy(0,itrk);
  Float_t muonMass = 0.10555;
  Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
  Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
  Float_t bl_y = sqrt(1. - bl_z*bl_z);
  Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
  
  return (nucleonMass*muonEnergy - muonMass*muonMass/2.)
    /(nucleonMass - muonEnergy + muonMomentum*costhbl);  

}

//********************************************************

Double_t MoqQuantities::VxB(Int_t itrk){
  
  double mom = RecoMuQP(itrk); //this loads ntpTrack itrk
  if(fabs(mom)<0.0001) return -1;  //check for no track, q/p = 0
  mom = 1./fabs(mom);
  TVector3 position(ntpTrack->vtx.u,
                    ntpTrack->vtx.v,
                    ntpTrack->vtx.z);
  TVector3 momentum(mom*ntpTrack->vtx.dcosu,
                    mom*ntpTrack->vtx.dcosv,
                    mom*ntpTrack->vtx.dcosz);
  
  if(mom>0.05) return SwimToGetVxB(position,momentum,
                                   GetChargeSign(itrk),
                                   double(ntpTrack->end.z));
  return -1;

}

//*******************************************************

Double_t MoqQuantities::SwimToGetVxB(TVector3 pos,TVector3 mom,
                                    Int_t chg,Double_t z){
  
  SwimSwimmer ss(ntpHeader->GetVldContext());
  double muMass = 0.105658357*Munits::GeV;
  SwimParticle swimMu(pos,mom,muMass,chg);
  SwimZCondition zc(z);
  Bool_t done = ss.SwimForward(swimMu,zc);
  if(done) return swimMu.GetVxB();

  return swimMu.GetVxB(); //might not want to do this later

}

//*******************************************************
void MoqQuantities::MakeHistoFile(std::string tag){

  std::string savename =  tag + ".ed.root";
  TFile save(savename.c_str(),"RECREATE"); 
  save.cd();
  
    //#declare lots of histos:
  TH1F *TrueEnergy = new TH1F("TrueEnergy","True Neutrino Energy",200,0,50);

  TH1F *RecoEnergy = new TH1F("RecoEnergy","Reco Neutrino Energy",200,0,50);

  TH1F *TrueMuEn = new TH1F("TrueMuEn","True Muon Energy",200,0,50);

  TH1F *RecoMuEn = new TH1F("RecoMuEn","Reco Muon Energy",200,0,50);
  
  TH1F *TrueShwEn = new TH1F("TrueShwEn","True Shower Energy",200,0,50);

  TH1F *RecoShwEn = new TH1F("RecoShwEn","Reco Shower Energy",200,0,50);
  
  TH1F *TrueMuQP = new TH1F("TrueMuQP","True Muon q/p",300,-3,3);

  TH1F *RecoMuQP = new TH1F("RecoMuQP","Reconstructed Muon q/p",300,-3,3);
  
  TH2F *TrueVsReco = new TH2F("TrueVsReco",
                              "Reco Vs True Neutrino Energy",
                              200,0,50,200,0,50);

  TH2F *TrueVsRecoMu = new TH2F("TrueVsRecoMu",
                                "Reco Vs True Muon Energy",200,0,50,200,0,50);

  TH2F *TrueVsRecoShw = new TH2F("TrueVsRecoShw",
                                 "Reco Vs True Shower Energy",
                                 200,0,50,200,0,50);
  
  TH2F *TrueRecoDiff 
    = new TH2F("TrueRecoDiff",
               "(True - Reco)/True Neutrino Energy vs True Neutrino Energy",
               200,0,50,1000,-9,1);

  //  TH2F *TrueRecoDiffMuCurv 
  //  = new TH2F("TrueRecoDiffMuCurv",
  //       "(True - Reco_{Curv})/True Muon Energy vs True Muon Energy",
  //       200,0,50,1000,-9,1);
  


  TH2F *TrueRecoDiffShw 
    = new TH2F("TrueRecoDiffShw",
               "(True - Reco)/True Shower Energy vs True Shower Energy",
               200,0,50,1000,-9,1);
  

  //  TH2F *TrueCurvMu 
  //= new TH2F("TrueCurvMu",
  //       "Reco Muon Energy from Curvature vs True Muon Energy",
  //       200,0,50,200,0,50); 


  for(int i=0;i<Nentries;i++){
    
    if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
                            << "% done" << std::endl;

    if(!GetEntry(i)) continue;
    
  
    
      if(this->PassCuts(0)){
          //neutrino
        TrueEnergy->Fill(this->TrueNuEnergy());
        RecoEnergy->Fill(this->RecoMuEnergy(1,0) + this->RecoShwEnergy());
        TrueVsReco->Fill(this->TrueNuEnergy(),this->RecoMuEnergy(1,0) 
                         + this->RecoShwEnergy());
        TrueRecoDiff->Fill(this->TrueNuEnergy(),
                           (this->TrueNuEnergy()-(this->RecoMuEnergy(1,0)
                                                  +this->RecoShwEnergy()))
                           /this->TrueNuEnergy());
        

        //shower          
        TrueShwEn->Fill(this->TrueShwEnergy());
        RecoShwEn->Fill(this->RecoShwEnergy());
        TrueVsRecoShw->Fill(this->TrueShwEnergy(),this->RecoShwEnergy());
        TrueRecoDiffShw->Fill(this->TrueShwEnergy(),
                              (this->TrueShwEnergy()
                               -this->RecoShwEnergy())
                              /this->TrueShwEnergy());


        //muon
        TrueMuEn->Fill(this->TrueMuEnergy());
        RecoMuEn->Fill(this->RecoMuEnergy(0,0));
        TrueVsRecoMu->Fill(this->TrueMuEnergy(),this->RecoMuEnergy(1,0));
        
        //curvature
        //TrueRecoDiffMuCurv->Fill(this->TrueMuEnergy(),
        //                       (this->TrueMuEnergy()
        //                -this->RecoMuEnergy(1,0))
        //               /this->TrueMuEnergy());
        //TrueCurvMu->Fill(this->TrueMuEnergy(),this->RecoMuEnergy(1,0));
        
        TrueMuQP->Fill(this->TrueMuQP());
        RecoMuQP->Fill(this->RecoMuQP(0));
        
      }


  }

 //#close file
  save.Write();
  save.Close();
}
    


#endif // #ifdef moqquantities_cxx

---