MadQuantities Class Reference

#include <MadQuantities.h>

Inheritance diagram for MadQuantities:
MadBase MadAnalysis NuMadAnalysis MadCBSQEAnalysis MadCluAnalysis MadDpAnalysis MadEdAnalysis MadMKAnalysis MadPIDAnalysis MadTestAnalysis MadTVAnalysis MadUserAnalysis

List of all members.

Public Member Functions

 MadQuantities (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadQuantities (JobC *, string, int)
 ~MadQuantities ()
Bool_t IsTrack ()
Bool_t IsSingleTrack ()
Bool_t PassTrackCuts (Int_t itrk)
Bool_t PassCuts (Int_t itrk)
Bool_t PassCuts ()
Bool_t IsFidVtx (Int_t itrk=0)
Bool_t IsFidAll (Int_t itrk=0)
Bool_t IsFidVtxEvt (Int_t ievt=0)
Bool_t IsFidVtxTrk (Int_t itrk=0)
Bool_t IsFidAllEvt (Int_t ievt=0)
Bool_t IsFidFD_AB (Int_t ievt=0)
Bool_t IsFid_2008 (Int_t ievt=0)
Float_t * CCNCSepVars (Int_t itrk=0)
Bool_t PassQuasiCuts (Int_t itrk=0, Float_t frac=0.1)
Int_t INu (Int_t itr=0)
Int_t INuNoOsc (Int_t itr=0)
Int_t Flavour (Int_t itr=0)
Int_t IAction (Int_t itr=0)
Int_t IResonance (Int_t itr=0)
Float_t X (Int_t itr=0)
Float_t Y (Int_t itr=0)
Float_t Q2 (Int_t itr=0)
Float_t W2 (Int_t itr=0)
Int_t Initial_state (Int_t itr=0)
Int_t Nucleus (Int_t itr=0)
Int_t HadronicFinalState (Int_t itr=0)
Float_t * Target4Vector (Int_t itr=0)
Float_t TrueNuEnergy (Int_t itr=0)
Float_t * TrueNuMom (Int_t itr=0)
Float_t * TrueVtx (Int_t itr=0)
Float_t TrueMuEnergy (Int_t itr=0)
Float_t TrueLeptonEnergy (Int_t itr=0)
Float_t TrueMuQP (Int_t itr=0)
Float_t TrueShwEnergy (Int_t itr=0)
Float_t GetTrueShwSC (Int_t itr=0)
Float_t TrueMuDCosZ (Int_t itr=0)
Float_t TrueMuDCosNeu (Int_t itr=0)
Int_t NumFSGeantino (Int_t itr=0)
Float_t GeantinoEnergy (Int_t itr=0)
Int_t NumFSProton (Int_t itr=0)
Float_t TotFSProtonEnergy (Int_t itr=0)
Float_t MaxFSProtonEnergy (Int_t itr=0)
Int_t NumFSNeutron (Int_t itr=0)
Float_t TotFSNeutronEnergy (Int_t itr=0)
Float_t MaxFSNeutronEnergy (Int_t itr=0)
Int_t NumFSPion (Int_t itr=0, Int_t pm=0)
Float_t TotFSPionEnergy (Int_t itr=0, Int_t pm=0)
Float_t MaxFSPionEnergy (Int_t itr=0, Int_t pm=0)
Int_t NumFSKaon (Int_t itr=0)
Float_t TotFSKaonEnergy (Int_t itr=0)
Float_t MaxFSKaonEnergy (Int_t itr=0)
Int_t NumFSPiZero (Int_t itr=0)
Float_t TotFSPiZeroEnergy (Int_t itr=0)
Float_t MaxFSPiZeroEnergy (Int_t itr=0)
Float_t RecoMuEnergy (Int_t opt=0, Int_t itrk=0, Bool_t data=true, Detector::Detector_t det=Detector::kNear)
Float_t RecoMuEnergyNew (Int_t opt, Int_t itrk, VldContext cx, ReleaseType::Release_t reltype, EnergyCorrections::WhichCorrection_t corrver)
Float_t RecoMuQP (Int_t itrk=0)
Float_t GetNonTrkSC (Int_t itrk=0, Int_t npln=0)
Float_t GetShwSC (Int_t itrk=0)
Float_t RecoShwEnergy (Int_t ishw=0, Int_t opt=-1, Int_t det=1, bool isdata=true, int mode=1)
Float_t RecoShwEnergyNew (Int_t ishw, Int_t opt, VldContext cx, ReleaseType::Release_t reltype, EnergyCorrections::WhichCorrection_t corrver)
Float_t RecoShwEnergySqrt (Int_t ievt=0)
Float_t GetSqrtTrkSkimSC (Int_t ievt=0)
Float_t GetSqrtShwSC (Int_t ievt=0)
Int_t GetChargeSign (Int_t itrk=0)
Float_t RecoQENuEnergy (Int_t itrk=0)
Float_t RecoQEQ2 (Int_t itrk=0)
Float_t RecoMuDCosZ (Int_t itr=0)
Float_t RecoMuDCosNeu (Int_t itr=0)
Double_t RecoEMFrac (Int_t shower, Int_t method, Double_t &EMFrac0, Double_t &EMFrac1)
Int_t * GetNShowerStrip (Int_t event, Int_t opt=0)
Float_t * ClosestDistanceToEvent (Int_t event, Float_t timeWeight=1)
void MakeValidationFile (std::string)
void ShowerValidation (char *fileName, Int_t startEnt=0)

Protected Member Functions

void SetStdHepVars (Int_t itr)

Protected Attributes

Int_t fCurStdHepEntry
Int_t fCurTruthIndex
Int_t fNumFSGeantino
Float_t fGeantinoEnergy
Int_t fNumFSProton
Float_t fTotFSProtonEnergy
Float_t fMaxFSProtonEnergy
Int_t fNumFSNeutron
Float_t fTotFSNeutronEnergy
Float_t fMaxFSNeutronEnergy
Int_t fNumFSPiPlus
Float_t fTotFSPiPlusEnergy
Float_t fMaxFSPiPlusEnergy
Int_t fNumFSPiMinus
Float_t fTotFSPiMinusEnergy
Float_t fMaxFSPiMinusEnergy
Int_t fNumFSPiZero
Float_t fTotFSPiZeroEnergy
Float_t fMaxFSPiZeroEnergy
Int_t fNumFSKaon
Float_t fTotFSKaonEnergy
Float_t fMaxFSKaonEnergy

Detailed Description

Definition at line 15 of file MadQuantities.h.


Constructor & Destructor Documentation

MadQuantities::MadQuantities ( TChain *  chainSR = 0,
TChain *  chainMC = 0,
TChain *  chainTH = 0,
TChain *  chainEM = 0 
)

Definition at line 26 of file MadQuantities.cxx.

References MadBase::Clear(), MadBase::emrecord, fCurStdHepEntry, fCurTruthIndex, fGeantinoEnergy, fMaxFSKaonEnergy, fMaxFSNeutronEnergy, fMaxFSPiMinusEnergy, fMaxFSPiPlusEnergy, fMaxFSPiZeroEnergy, fMaxFSProtonEnergy, fNumFSGeantino, fNumFSKaon, fNumFSNeutron, fNumFSPiMinus, fNumFSPiPlus, fNumFSPiZero, fNumFSProton, fTotFSKaonEnergy, fTotFSNeutronEnergy, fTotFSPiMinusEnergy, fTotFSPiPlusEnergy, fTotFSPiZeroEnergy, fTotFSProtonEnergy, MadBase::InitChain(), MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.

00028    : MadBase(chainSR, chainMC, chainTH, chainEM)
00029 {
00030 
00031   if(!chainSR) {
00032     record = 0;
00033     strecord = 0;
00034     emrecord = 0;
00035     mcrecord = 0;
00036     threcord = 0;
00037     Clear();
00038     whichSource = -1;
00039     isMC = true;
00040     isTH = true;
00041     isEM = true;
00042     Nentries = -1;
00043     return;
00044   }
00045 
00046   fNumFSGeantino = 0;
00047   fGeantinoEnergy = 0;
00048   fNumFSProton = 0;
00049   fTotFSProtonEnergy = 0;
00050   fMaxFSProtonEnergy = 0;
00051   fNumFSNeutron = 0;
00052   fTotFSNeutronEnergy = 0;
00053   fMaxFSNeutronEnergy = 0;
00054   fNumFSPiPlus = 0;
00055   fTotFSPiPlusEnergy = 0;
00056   fMaxFSPiPlusEnergy = 0;
00057   fNumFSPiMinus = 0;
00058   fTotFSPiMinusEnergy = 0;
00059   fMaxFSPiMinusEnergy = 0;
00060   fNumFSPiZero = 0;
00061   fTotFSPiZeroEnergy = 0;
00062   fMaxFSPiZeroEnergy = 0;
00063   fNumFSKaon = 0;
00064   fTotFSKaonEnergy = 0;
00065   fMaxFSKaonEnergy = 0;
00066   fCurStdHepEntry = -1;
00067   fCurTruthIndex = -1;
00068 
00069   InitChain(chainSR,chainMC,chainTH,chainEM);
00070   whichSource = 0;
00071 
00072 }

MadQuantities::MadQuantities ( JobC j,
string  path,
int  entries 
)

Definition at line 76 of file MadQuantities.cxx.

References MadBase::Clear(), MadBase::emrecord, MadBase::fChain, fGeantinoEnergy, MadBase::fJC, fMaxFSKaonEnergy, fMaxFSNeutronEnergy, fMaxFSPiMinusEnergy, fMaxFSPiPlusEnergy, fMaxFSPiZeroEnergy, fMaxFSProtonEnergy, fNumFSGeantino, fNumFSKaon, fNumFSNeutron, fNumFSPiMinus, fNumFSPiPlus, fNumFSPiZero, fNumFSProton, fTotFSKaonEnergy, fTotFSNeutronEnergy, fTotFSPiMinusEnergy, fTotFSPiPlusEnergy, fTotFSPiZeroEnergy, fTotFSProtonEnergy, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.

00077    : MadBase(j,path,entries)
00078 {
00079   record = 0;
00080   strecord = 0;
00081   emrecord = 0;
00082   mcrecord = 0;
00083   threcord = 0;
00084   Clear();
00085   isMC = true;
00086   isTH = true;
00087   isEM = true;
00088   Nentries = entries;
00089   whichSource = 1;
00090   jcPath = path;
00091   fJC = j;
00092   fChain = NULL;
00093 
00094   fNumFSGeantino = 0;
00095   fGeantinoEnergy = 0;
00096   fNumFSProton = 0;
00097   fTotFSProtonEnergy = 0;
00098   fMaxFSProtonEnergy = 0;
00099   fNumFSNeutron = 0;
00100   fTotFSNeutronEnergy = 0;
00101   fMaxFSNeutronEnergy = 0;
00102   fNumFSPiPlus = 0;
00103   fTotFSPiPlusEnergy = 0;
00104   fMaxFSPiPlusEnergy = 0;
00105   fNumFSPiMinus = 0;
00106   fTotFSPiMinusEnergy = 0;
00107   fMaxFSPiMinusEnergy = 0;
00108   fNumFSPiZero = 0;
00109   fTotFSPiZeroEnergy = 0;
00110   fMaxFSPiZeroEnergy = 0;
00111   fNumFSKaon = 0;
00112   fTotFSKaonEnergy = 0;
00113   fMaxFSKaonEnergy = 0;
00114 
00115 }

MadQuantities::~MadQuantities (  ) 

Definition at line 119 of file MadQuantities.cxx.

00120 {
00121 
00122 }


Member Function Documentation

Float_t * MadQuantities::CCNCSepVars ( Int_t  itrk = 0  ) 

Definition at line 425 of file MadQuantities.cxx.

References NtpSRTrack::ds, MadBase::isST, MadBase::LoadTrack(), NtpSRTrack::nstrip, MadBase::ntpStrip, MadBase::ntpTrack, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, MadBase::record, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRRecord::stp, NtpSRTrack::stp, and MadBase::strecord.

Referenced by MadEdAnalysis::DataHist(), MadCBSQEAnalysis::LikeliQE(), MadCBSQEAnalysis::MakeQEFile(), MakeValidationFile(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyLikeliQE(), MadEdAnalysis::MyMakeQEFile(), and MadCBSQEAnalysis::TestQEDiscrim().

00425                                              {
00426   
00427   Float_t *sepVars = new Float_t[3];
00428   sepVars[0] = sepVars[1] = sepVars[2] = 0.0;
00429   
00430   if(!LoadTrack(itrk)) return sepVars;
00431 
00432   Int_t numtrkstp = ntpTrack->nstrip;
00433   Int_t *trkstrips = ntpTrack->stp;
00434   Float_t planeCharge[500] = {};
00435   Float_t totTrkCharge = 0;
00436   Int_t firstPlane = 500;
00437   Int_t lastPlane = 0;
00438 
00439   TClonesArray* pointStripArray = NULL;
00440   if(isST) pointStripArray = (strecord->stp);
00441   else pointStripArray = (record->stp);
00442   TClonesArray& stripArray = *pointStripArray;
00443 
00444   for(Int_t i=0;i<numtrkstp;i++){
00445     Int_t index = trkstrips[i];
00446     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
00447     totTrkCharge += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
00448     Int_t tmpPln = ntpStrip->plane;
00449     planeCharge[tmpPln] += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; 
00450     if(ntpStrip->plane<firstPlane) firstPlane = ntpStrip->plane;
00451     if(ntpStrip->plane>lastPlane) lastPlane = ntpStrip->plane;
00452   }
00453   
00454   Int_t numPlanes = (lastPlane-firstPlane+1);
00455   Float_t chargePerPlane = totTrkCharge/ntpTrack->ds;
00456   Float_t firstHalf = 0;
00457   Float_t lastHalf = 0.00001;
00458   
00459   for(Int_t i=firstPlane;i<=lastPlane;i++){
00460     if(numPlanes%2==0){
00461       if(i<firstPlane+Int_t(numPlanes/2)) firstHalf+=planeCharge[i];
00462       else lastHalf+=planeCharge[i];
00463     }
00464     else {
00465       if(i<firstPlane+Int_t(numPlanes/2)) firstHalf+=planeCharge[i];
00466       else if(i==firstPlane+Int_t(numPlanes/2)+1) {
00467         firstHalf+=planeCharge[i];
00468         lastHalf+=planeCharge[i];
00469       }
00470       else lastHalf+=planeCharge[i];
00471     }
00472   }
00473   Float_t halfRatio = firstHalf/lastHalf;
00474   
00475   sepVars[0] = ntpTrack->ds;
00476   sepVars[1] = chargePerPlane;
00477   sepVars[2] = halfRatio;
00478   
00479   return sepVars;
00480 
00481 }

Float_t * MadQuantities::ClosestDistanceToEvent ( Int_t  event,
Float_t  timeWeight = 1 
)

Definition at line 1756 of file MadQuantities.cxx.

References MadBase::eventSummary, MadBase::LoadEvent(), MadBase::LoadStrip(), NtpSREventSummary::nstrip, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpStrip, NtpSRStrip::planeview, NtpSREvent::stp, NtpSRStrip::time0, NtpSRStrip::time1, NtpSRStrip::tpos, and NtpSRStrip::z.

Referenced by ShowerValidation().

01757 {
01758   Float_t *values = new Float_t[3];
01759   values[0] = values[1] = values[2] = -1;
01760   if(!LoadEvent(evt)) return values;
01761   Float_t smallestMetricValue = 0;
01762   Float_t smallestDeltaPos = 0;
01763   Float_t smallestDeltaTime = 0;
01764   for(int i=0;i<ntpEvent->nstrip;i++){
01765     if(!LoadStrip(ntpEvent->stp[i])) continue;
01766     Float_t z_val = ntpStrip->z;
01767     Float_t tpos_val = ntpStrip->tpos;
01768     Double_t time0_val = 1e9*ntpStrip->time0;
01769     Double_t time1_val = 1e9*ntpStrip->time1;    
01770     Double_t time_val = 0;
01771     if(time0_val>0 && time1_val>0) time_val = (time0_val + time1_val)/2.;
01772     else if(time0_val>0) time_val = time0_val;
01773     else if(time1_val>0) time_val = time1_val;
01774     Int_t planeview = ntpStrip->planeview;
01775     Float_t minDeltaStrip = 999999;
01776     Float_t minDeltaStripPos = 999999;
01777     Double_t minDeltaStripTime = 999999;
01778     for(UInt_t j=0;j<eventSummary->nstrip;j++){
01779       Bool_t useStrip = true;
01780       for(int k=0;k<ntpEvent->nstrip;k++){
01781         if(ntpEvent->stp[k]==Int_t(j)) {
01782           useStrip=false;
01783           break;
01784         }
01785       }
01786       if(!useStrip) continue;
01787       if(!LoadStrip(j)) continue;
01788       if(ntpStrip->planeview!=planeview) continue;
01789       Float_t deltaStripPos = TMath::Sqrt(TMath::Power(z_val - ntpStrip->z,2) +
01790                                           TMath::Power(tpos_val = ntpStrip->tpos,2));
01791       Double_t t0_val = 1e9*ntpStrip->time0; 
01792       Double_t t1_val = 1e9*ntpStrip->time1;
01793       Double_t t_val = 0;
01794       if(t0_val>0 && t1_val>0) t_val = (t0_val + t1_val)/2.;
01795       else if(t0_val>0) t_val = t0_val;
01796       else if(t1_val>0) t_val = t1_val;      
01797       Double_t deltaStripTime = TMath::Abs(time_val-t_val);
01798       Double_t metricValue = deltaStripPos + timeWeight*Float_t(deltaStripTime)*0.3;
01799       if(metricValue<minDeltaStrip) {
01800         minDeltaStrip = metricValue;
01801         minDeltaStripTime = deltaStripTime;
01802         minDeltaStripPos = deltaStripPos;
01803       }
01804     }
01805     smallestMetricValue += minDeltaStrip;
01806     smallestDeltaTime += Float_t(minDeltaStripTime);
01807     smallestDeltaPos += minDeltaStripPos;
01808   }
01809   values[0] = smallestMetricValue/Float_t(ntpEvent->nstrip);
01810   values[1] = smallestDeltaPos/Float_t(ntpEvent->nstrip);
01811   values[2] = smallestDeltaTime/Float_t(ntpEvent->nstrip);
01812   return values;
01813 }

Int_t MadQuantities::Flavour ( Int_t  itr = 0  ) 

Definition at line 506 of file MadQuantities.cxx.

References INu().

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadDpID::FillPDFs(), MadTestAnalysis::MakeMyFile(), MadPIDAnalysis::MakeMyFile(), MadDpAnalysis::MakeMyFile(), MadEdAnalysis::MyCreatePAN(), and MadEdAnalysis::MyMakeMyFile().

00506                                      {
00507   Int_t f2=abs(INu(itr));
00508   Int_t flavcode=0;
00509   switch (f2) {
00510   case 12:
00511     flavcode=1; // nu_e
00512     break;
00513   case 14:
00514     flavcode=2; // nu_mu
00515     break;
00516   case 16:
00517     flavcode=3; // nu_tau
00518     break;
00519   default:
00520     flavcode=0;
00521     break;
00522   }
00523   return flavcode;
00524 }

Float_t MadQuantities::GeantinoEnergy ( Int_t  itr = 0  ) 

Definition at line 634 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fGeantinoEnergy, MadBase::lastEntry, and SetStdHepVars().

Referenced by MakeValidationFile().

00634                                               {
00635   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00636     SetStdHepVars(itr);
00637   return fGeantinoEnergy;
00638 }

Int_t MadQuantities::GetChargeSign ( Int_t  itrk = 0  ) 

Definition at line 1643 of file MadQuantities.cxx.

References RecoMuQP().

Referenced by RecoQENuEnergy().

01643                                             {  
01644   if(RecoMuQP(itrk)>0) return 1;
01645   return -1;
01646 }

Float_t MadQuantities::GetNonTrkSC ( Int_t  itrk = 0,
Int_t  npln = 0 
)

Definition at line 1441 of file MadQuantities.cxx.

References MadBase::eventSummary, MadBase::isST, MadBase::LoadStrip(), MadBase::LoadTrack(), NtpSREventSummary::nstrip, NtpSRTrack::nstrip, MadBase::ntpStrip, MadBase::ntpTrack, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRVertex::plane, MadBase::record, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRRecord::stp, NtpSRTrack::stp, MadBase::strecord, and NtpSRTrack::vtx.

01441                                                        {
01442   
01443   Float_t Summed_Stp_SC = 0;
01444   Float_t Summed_Trk_SC = 0;
01445   
01446   bool useAll = false;
01447   if(npln<=0) useAll=true;  //Most basic thing in the world ever: 
01448                             //if it's not in the track it's in the shower
01449   
01450   TClonesArray* pointStripArray = NULL;
01451   if(isST) pointStripArray = (strecord->stp);
01452   else pointStripArray = (record->stp);
01453   TClonesArray& stripArray = *pointStripArray;
01454   if(LoadTrack(itrk)) {
01455     Int_t numtrkstp = ntpTrack->nstrip;
01456     for(Int_t i=0;i<numtrkstp;i++){
01457       Int_t index = ntpTrack->stp[i];
01458       ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01459       if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
01460           && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
01461         Summed_Trk_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01462       }
01463     }
01464   }
01465   
01466   for(Int_t i=0;i<int(eventSummary->nstrip);i++) {
01467     if(!LoadStrip(i)) continue;
01468     if((ntpStrip->plane <= ntpTrack->vtx.plane + npln
01469         && ntpStrip->plane >= ntpTrack->vtx.plane - npln)||useAll){
01470       Summed_Stp_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01471     }
01472   }
01473   return Summed_Stp_SC - Summed_Trk_SC;
01474 }

Int_t * MadQuantities::GetNShowerStrip ( Int_t  event,
Int_t  opt = 0 
)

Definition at line 1735 of file MadQuantities.cxx.

References NtpSRShower::clu, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadShower(), NtpSRShower::ncluster, NtpSRCluster::nstrip, MadBase::ntpCluster, MadBase::ntpShower, and NtpSRCluster::planeview.

Referenced by MadCluAnalysis::CallUserFunctions(), and ShowerValidation().

01736 {
01737   Int_t *NShowerStrip = new Int_t[2];
01738   NShowerStrip[0] = 0; NShowerStrip[1] = 0;
01739   if(!LoadShower(shower)) {
01740     NShowerStrip[0] = -10; NShowerStrip[1] = -10;
01741     return NShowerStrip;
01742   }
01743   Int_t *clusters = ntpShower->clu;  
01744   for(int l=0;l<ntpShower->ncluster;l++){    
01745     Int_t index = clusters[l];
01746     if(!LoadCluster(index)) continue;
01747     if(opt==0||ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3) {
01748       if(ntpCluster->planeview==2) NShowerStrip[0] += ntpCluster->nstrip;
01749       else if(ntpCluster->planeview==3) NShowerStrip[1] += ntpCluster->nstrip;
01750     }
01751   }
01752   return NShowerStrip;
01753 }

Float_t MadQuantities::GetShwSC ( Int_t  itrk = 0  ) 

Definition at line 1478 of file MadQuantities.cxx.

References MadBase::eventSummary, MadBase::isST, MadBase::LoadShower(), MadBase::LoadTrack(), NtpSREventSummary::nshower, NtpSRShower::nstrip, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTrack, NtpSRStrip::ph0, NtpSRStrip::ph1, MadBase::record, NtpSRPulseHeight::sigcor, NtpSRShower::stp, NtpStRecord::stp, NtpSRRecord::stp, MadBase::strecord, NtpSRShower::vtx, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by MakeValidationFile().

01478                                          {
01479 
01480   //Try something a bit smarter:
01481   //Look for a reconstructed shower near to the track vtx and sum the SigCors
01482   if(!LoadTrack(itrk)) return 0.;
01483   
01484   Float_t Summed_Shw_SC = 0;
01485   Int_t bestshower = -1;
01486   Float_t bestdistance = 1000.;
01487   if(eventSummary->nshower!=0){
01488     for(Int_t i=0;i<eventSummary->nshower;i++){
01489       if(!LoadShower(i)) continue;
01490       Float_t distance = sqrt((ntpTrack->vtx.x - ntpShower->vtx.x)
01491                               *(ntpTrack->vtx.x - ntpShower->vtx.x)
01492                               +(ntpTrack->vtx.y - ntpShower->vtx.y)
01493                               *(ntpTrack->vtx.y - ntpShower->vtx.y)
01494                               +(ntpTrack->vtx.z - ntpShower->vtx.z)
01495                               *(ntpTrack->vtx.z - ntpShower->vtx.z));
01496       
01497       if(distance<0.5){
01498         if(bestshower==-1) {
01499           bestshower = i;
01500           bestdistance = distance;
01501         }
01502         else if(distance<bestdistance) bestshower = i;
01503       }
01504     }
01505     
01506     TClonesArray* pointStripArray = NULL;
01507     if(isST) pointStripArray = (strecord->stp);
01508     else pointStripArray = (record->stp);
01509     TClonesArray& stripArray = *pointStripArray;
01510     if(bestshower!=-1){
01511       LoadShower(bestshower);
01512       Int_t *shwstrips = ntpShower->stp;
01513       for(Int_t j=0;j<ntpShower->nstrip;j++){
01514         Int_t index = shwstrips[j];
01515         ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01516         Summed_Shw_SC += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01517       }
01518       return Summed_Shw_SC;
01519     }
01520     return 0;
01521   }
01522   else return 0; //no reconstructed shower
01523 
01524 }

Float_t MadQuantities::GetSqrtShwSC ( Int_t  ievt = 0  ) 

Definition at line 1622 of file MadQuantities.cxx.

References MadBase::isST, MadBase::LoadLargestShowerFromEvent(), NtpSRShower::nstrip, MadBase::ntpShower, MadBase::ntpStrip, NtpSRStrip::ph0, NtpSRStrip::ph1, MadBase::record, NtpSRPulseHeight::sigcor, NtpSRShower::stp, NtpStRecord::stp, NtpSRRecord::stp, and MadBase::strecord.

Referenced by RecoShwEnergySqrt().

01622                                              {
01623 
01624   Int_t shower = -1;
01625   if(!LoadLargestShowerFromEvent(ievt,shower)) return 0;
01626 
01627   Float_t summedSqrtSC = 0;
01628   TClonesArray* pointStripArray = NULL;
01629   if(isST) pointStripArray = (strecord->stp);
01630   else pointStripArray = (record->stp);
01631   TClonesArray& stripArray = *pointStripArray;
01632   Int_t *shwstrips = ntpShower->stp;
01633   for(Int_t j=0;j<ntpShower->nstrip;j++){
01634     Int_t index = shwstrips[j];
01635     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01636     summedSqrtSC += sqrt(ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
01637   }
01638   return summedSqrtSC;  
01639 }

Float_t MadQuantities::GetSqrtTrkSkimSC ( Int_t  ievt = 0  ) 

Definition at line 1594 of file MadQuantities.cxx.

References NtpSRVertex::dcosz, NtpSRStripPulseHeight::gev, MadBase::isST, MadBase::LoadLargestTrackFromEvent(), NtpSRTrack::nstrip, MadBase::ntpStrip, MadBase::ntpTrack, NtpSRTrack::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRVertex::plane, MadBase::record, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRRecord::stp, NtpSRTrack::stp, MadBase::strecord, and NtpSRTrack::vtx.

Referenced by RecoShwEnergySqrt().

01594                                                  {
01595 
01596   Float_t skimmedSqrtSC = 0;
01597   Int_t track = -1;
01598   if(!LoadLargestTrackFromEvent(ievt,track)) return 0;
01599 
01600   if(ntpTrack->ph.gev>0){
01601     TClonesArray* pointStripArray = NULL;
01602     if(isST) pointStripArray = (strecord->stp);
01603     else pointStripArray = (record->stp);
01604     TClonesArray& stripArray = *pointStripArray;
01605     Int_t numtrkstp = ntpTrack->nstrip;
01606     for(Int_t i=0;i<numtrkstp;i++){
01607       Int_t index = ntpTrack->stp[i];
01608       ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray[index]);
01609       if(ntpStrip->plane<ntpTrack->vtx.plane+11){
01610         skimmedSqrtSC += ((ntpStrip->ph0.sigcor + 
01611                            ntpStrip->ph1.sigcor) > 1600/ntpTrack->vtx.dcosz) ? 
01612           sqrt((ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor) - 
01613                (1600/ntpTrack->vtx.dcosz)) : 0;
01614       }
01615     }
01616   }
01617   return skimmedSqrtSC;  
01618 }

Float_t MadQuantities::GetTrueShwSC ( Int_t  itr = 0  ) 

Definition at line 1320 of file MadQuantities.cxx.

References MadBase::LoadTruth().

Referenced by MakeValidationFile().

01320                                             {
01321 
01322   if(!LoadTruth(itr)) return 0.;
01323   //Requires truth digit information which is not present in ntuples  
01324   Float_t Summed_Shw_SC = 0;
01325   
01326   return Summed_Shw_SC;
01327 }

Int_t MadQuantities::HadronicFinalState ( Int_t  itr = 0  ) 

Definition at line 1175 of file MadQuantities.cxx.

References NtpMCStdHep::IdHEP, IResonance(), MadBase::isST, NtpMCStdHep::IstHEP, MadBase::LoadStdHep(), MadBase::LoadTruth(), NtpMCStdHep::mc, MadBase::mcrecord, MadBase::ntpStdHep, NtpStRecord::stdhep, NtpMCRecord::stdhep, and MadBase::strecord.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), and MadEdAnalysis::MyCreatePAN().

01175                                                 {
01176   if(!LoadTruth(itr)) return 0;
01177   Int_t hfs=0;
01178   Int_t proc = IResonance(itr);
01179   TClonesArray* pointStdhepArray = NULL;
01180   if(isST) pointStdhepArray = (strecord->stdhep);
01181   else pointStdhepArray = (mcrecord->stdhep);
01182   TClonesArray& stdhepArray = *pointStdhepArray;
01183   Int_t nStdHep = stdhepArray.GetEntries();
01184   if(proc==1002){
01185     for(int i=0;i<nStdHep;i++){    
01186       LoadStdHep(i);
01187       if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3 &&
01188          !(abs(ntpStdHep->IdHEP)==15)){  //not a tau lepton
01189         hfs = ntpStdHep->IdHEP;
01190         break;
01191       }
01192     }
01193     hfs = hfs%1000;
01194   }
01195   else {
01196     for(int i=0;i<nStdHep;i++){    
01197       LoadStdHep(i);
01198       if(ntpStdHep->mc==itr && ntpStdHep->IstHEP==3){
01199         hfs = ntpStdHep->IdHEP;
01200         break;
01201       }
01202     }
01203     hfs = hfs%1000;
01204   }
01205   return hfs;
01206 }

Int_t MadQuantities::IAction ( Int_t  itr = 0  ) 
Int_t MadQuantities::Initial_state ( Int_t  itr = 0  ) 

Definition at line 1128 of file MadQuantities.cxx.

References NtpMCStdHep::IdHEP, MadBase::isST, NtpMCStdHep::IstHEP, MadBase::LoadStdHep(), MadBase::LoadTruth(), NtpMCStdHep::mc, MadBase::mcrecord, MadBase::ntpStdHep, NtpStRecord::stdhep, NtpMCRecord::stdhep, and MadBase::strecord.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::MyCreatePAN(), and MadAnalysis::RecoMC().

01128                                            {
01129   
01130   if(!LoadTruth(itr)) return 0;
01131   
01132   Int_t initial_state=0;  
01133   TClonesArray* pointStdhepArray = NULL;
01134   if(isST) pointStdhepArray = (strecord->stdhep);
01135   else pointStdhepArray = (mcrecord->stdhep);
01136   TClonesArray& stdhepArray = *pointStdhepArray;
01137   Int_t nStdHep = stdhepArray.GetEntries();
01138   
01139   int protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
01140   int nubarnu = 0;    // +1 = neutrino, -1 = antineutrino
01141   
01142   for(int i=0;i<nStdHep;i++){    
01143     LoadStdHep(i);
01144     if(ntpStdHep->mc==itr){
01145       
01146       if(ntpStdHep->IstHEP==0){  //initial state particle
01147         if(abs(ntpStdHep->IdHEP)==12 || 
01148            abs(ntpStdHep->IdHEP)==14 || 
01149            abs(ntpStdHep->IdHEP)==16){   //(anti)neutrino         
01150           nubarnu = ntpStdHep->IdHEP/abs(ntpStdHep->IdHEP);  //get sign
01151         }
01152       }
01153       if(ntpStdHep->IstHEP==11){    //target nucleon
01154         if(ntpStdHep->IdHEP==2212) protneut = 1;  //proton
01155         else if(ntpStdHep->IdHEP==2112) protneut = 0;  //neutron
01156         else if(abs(ntpStdHep->IdHEP)>1000000000) protneut = 2;  //nucleus
01157         else protneut = 3; //atom - probably never get here since IdHEP A>N?
01158       }
01159     }
01160   }
01161 
01162   if(protneut==1 && nubarnu==1)  initial_state=1;  //p + v
01163   if(protneut==0 && nubarnu==1)  initial_state=2;  //n + v
01164   if(protneut==1 && nubarnu==-1) initial_state=3;  //p + vbar
01165   if(protneut==0 && nubarnu==-1) initial_state=4;  //n + vbar
01166   if(protneut==2 && nubarnu==1)  initial_state=5;  //N + v
01167   if(protneut==3 && nubarnu==1)  initial_state=6;  //A + v
01168   if(protneut==2 && nubarnu==-1) initial_state=7;  //N + vbar
01169   if(protneut==3 && nubarnu==-1) initial_state=8;  //A + vbar
01170 
01171   return initial_state;
01172 }

Int_t MadQuantities::INu ( Int_t  itr = 0  ) 
Int_t MadQuantities::INuNoOsc ( Int_t  itr = 0  ) 

Definition at line 499 of file MadQuantities.cxx.

References NtpMCTruth::inunoosc, MadBase::LoadTruth(), and MadBase::ntpTruth.

Referenced by MadTVAnalysis::CreatePAN().

00499                                       {
00500   if(!LoadTruth(itr)) return 0;
00501   return ntpTruth->inunoosc;
00502 }

Int_t MadQuantities::IResonance ( Int_t  itr = 0  ) 
Bool_t MadQuantities::IsFid_2008 ( Int_t  ievt = 0  ) 

Definition at line 286 of file MadQuantities.cxx.

References choose_infid_set(), infid(), NueStandard::IsInFid(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, and MadBase::strecord.

Referenced by MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), and MadPIDAnalysis::MakeMyFile().

00287 { 
00288 
00289 //first up initialise the cuts to be the cc2008 version
00290 //just need to do this once
00291   static Bool_t infid_initialised=false;
00292   if (!infid_initialised) {
00293      infid_initialised=true;
00294      choose_infid_set("cc2008");
00295   }
00296 
00297   if(!LoadEvent(event)) return false;
00298 
00299   Int_t track = -1;
00300   LoadLargestTrackFromEvent(event,track);
00301 
00302 
00303   bool IsInFid=false;
00304 
00305 //  Then find out if the track is in the fiducial volume:
00306   if(track!=-1){
00307      //if it has a track, use the track vertex
00308      IsInFid=infid(*strecord,*ntpTrack);
00309   }
00310   else{
00311      // for events with no track, use event vertex
00312      IsInFid=infid(*strecord,*ntpEvent);
00313   }
00314 
00315   return IsInFid;
00316 
00317 }

Bool_t MadQuantities::IsFidAll ( Int_t  itrk = 0  ) 

Definition at line 324 of file MadQuantities.cxx.

References NtpSRTrack::contained, NtpSRTrack::end, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kNear, MadBase::LoadTrack(), MadBase::ntpHeader, MadBase::ntpTrack, NtpSRVertex::plane, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::DataHist(), MadCBSQEAnalysis::MakeQEFile(), MakeValidationFile(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyCreatePANData(), MadEdAnalysis::MyMakeQEFile(), RecoMuEnergy(), RecoMuEnergyNew(), and MadCBSQEAnalysis::TestQEDiscrim().

00324                                         {
00325   
00326   if(!LoadTrack(itrk)) return false;
00327 
00328   // 'empirical' definiton of CEV-like events
00329  
00330 
00331   // new definition - coil hole cuts removed for cedar
00332   bool contained=true;
00333   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00334     
00335     if (!(ntpTrack->end.z<15 && 
00336           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00337           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00338           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00339           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00340           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00341           ntpTrack->end.y>ntpTrack->end.x-3.55)) {contained=false;}
00342     
00343     //updated for 2008 analysis
00344     if (ntpTrack->end.x>1.3) contained=ntpTrack->contained;
00345   }
00346   else if(pow(ntpTrack->end.x,2)+pow(ntpTrack->end.y,2)>14 || ntpTrack->end.plane>475 || ntpTrack->end.plane<5) contained=false;
00347 
00348   return contained;
00349 
00350 
00351 
00352 
00353 // old definition - cuts out ND & FD events close to coil in cedar...
00354 
00355 /*
00356  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00357     
00358     if (!(ntpTrack->end.z<16 && sqrt(pow(ntpTrack->end.x,2)+
00359                                      pow(ntpTrack->end.y,2))>0.4 &&
00360           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00361           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00362           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00363           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00364           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00365           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00366      }      
00367      else if(ntpTrack->fidall.dr<0.5||ntpTrack->fidall.dz<0.5) return false;
00368  */
00369 
00370 
00371   return true;
00372 
00373 }

Bool_t MadQuantities::IsFidAllEvt ( Int_t  ievt = 0  ) 

Definition at line 376 of file MadQuantities.cxx.

References NtpSRTrack::end, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kNear, MadBase::LoadLargestTrackFromEvent(), MadBase::ntpHeader, MadBase::ntpTrack, NtpSRVertex::plane, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

00376                                            {
00377   Int_t track = -1;
00378   if(!LoadLargestTrackFromEvent(ievt,track)) return false;
00379   
00380   // 'empirical' definiton of CEV-like events
00381 
00382   // new definition - coil hole cuts removed for cedar
00383 
00384   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00385     
00386     if (!(ntpTrack->end.z<15 && 
00387           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00388           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00389           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00390           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00391           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00392           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00393     
00394 
00395   }
00396   else if(pow(ntpTrack->end.x,2)+pow(ntpTrack->end.y,2)>14 || ntpTrack->end.plane>475 || ntpTrack->end.plane<5) return false;
00397 
00398 
00399 
00400 
00401 // old definition - cuts out ND & FD events close to coil in cedar...
00402 
00403 /*
00404  if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) {//near det
00405     
00406     if (!(ntpTrack->end.z<16 && sqrt(pow(ntpTrack->end.x,2)+
00407                                      pow(ntpTrack->end.y,2))>0.4 &&
00408           ntpTrack->end.x<2.7 && ntpTrack->end.x>-1.65 && 
00409           ntpTrack->end.y<1.65 && ntpTrack->end.y>-1.65 &&
00410           ntpTrack->end.y>(-ntpTrack->end.x)-1.65 && 
00411           ntpTrack->end.y<ntpTrack->end.x+1.65 &&
00412           ntpTrack->end.y<(-ntpTrack->end.x)+3.55 && 
00413           ntpTrack->end.y>ntpTrack->end.x-3.55)) {return false;}
00414      }      
00415      else if(ntpTrack->fidall.dr<0.5||ntpTrack->fidall.dz<0.5) return false;
00416  */
00417 
00418   return true;
00419 
00420 }

Bool_t MadQuantities::IsFidFD_AB ( Int_t  ievt = 0  ) 

Definition at line 246 of file MadQuantities.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpTrack, NtpSRVertex::plane, NtpSREvent::vtx, NtpSRTrack::vtx, NtpSRVertex::x, and NtpSRVertex::y.

Referenced by MadPIDAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00247 { 
00248 
00249   if(!LoadEvent(event)) return false;
00250 
00251   Int_t track = -1;
00252   LoadLargestTrackFromEvent(event,track);
00253 
00254   // only check for FD
00255 
00256   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00257  
00258 
00259     // for events with no track, use event vertex
00260     
00261     if (track==-1 && (ntpEvent->vtx.plane<=4 || 
00262                       ntpEvent->vtx.plane>=466 ||     // ends
00263                       (ntpEvent->vtx.plane<=253 && 
00264                        ntpEvent->vtx.plane>=241) ||   // between SMs
00265                       pow(ntpEvent->vtx.x,2)+         // radial cut
00266                       pow(ntpEvent->vtx.y,2)>14)) return false;
00267     
00268     
00269     // otherwise, use track vertex
00270     
00271     if (track!=-1 && (ntpTrack->vtx.plane<=4 || 
00272                       ntpTrack->vtx.plane>=466 ||   //ends
00273                       (ntpTrack->vtx.plane<=253 && 
00274                        ntpTrack->vtx.plane>=241) ||  //between SMs
00275                       pow(ntpTrack->vtx.x,2)+           //radial cut
00276                       pow(ntpTrack->vtx.y,2)>14)) return false;
00277 
00278   }
00279 
00280   return true; 
00281 
00282 }

Bool_t MadQuantities::IsFidVtx ( Int_t  itrk = 0  ) 

Definition at line 167 of file MadQuantities.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadTrack(), MadBase::ntpHeader, MadBase::ntpTrack, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by MadCBSQEAnalysis::PassAnalysisCuts(), and PassCuts().

00167                                         {
00168   if(!LoadTrack(itrk)) return false;
00169   
00170   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00171     
00172     if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>29.4 ||   //ends
00173        (ntpTrack->vtx.z<16.5&&ntpTrack->vtx.z>14.5) ||  //between SMs
00174        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00175             +(ntpTrack->vtx.y*ntpTrack->vtx.y))>3.5 ||
00176        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00177             +(ntpTrack->vtx.y*ntpTrack->vtx.y))<0.4) return false;
00178     
00179   }
00180   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00181     if(ntpTrack->vtx.z<1 || ntpTrack->vtx.z>5 ||
00182        sqrt(((ntpTrack->vtx.x-1.4828)*(ntpTrack->vtx.x-1.4828)) +
00183             ((ntpTrack->vtx.y-0.2384)*(ntpTrack->vtx.y-0.2384)))>1) return false;
00184     //      sqrt(((ntpTrack->vtx.x-1.4885)*(ntpTrack->vtx.x-1.4885)) +
00185     //      ((ntpTrack->vtx.y-0.1397)*(ntpTrack->vtx.y-0.1397)))>1) return false;
00186   }
00187   
00188   return true;
00189 }

Bool_t MadQuantities::IsFidVtxEvt ( Int_t  ievt = 0  ) 

Definition at line 192 of file MadQuantities.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::ntpEvent, MadBase::ntpHeader, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadCluAnalysis::DataDistributions(), MadCluAnalysis::FOM(), MadCBSQEAnalysis::MakeQEFile(), MadTestAnalysis::PassBasicCuts(), MadTestAnalysis::PassFid(), MadPIDAnalysis::PassFid(), MadDpAnalysis::PassFid(), MadCluAnalysis::Plot(), ShowerValidation(), and MadCBSQEAnalysis::TestQEDiscrim().

00192                                            {
00193   if(!LoadEvent(ievt)) return false;
00194   
00195   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00196     
00197     if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 ||   //ends
00198        (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) ||  //between SMs
00199        sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00200             +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5 ||
00201        sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x)           //radial cut
00202             +(ntpEvent->vtx.y*ntpEvent->vtx.y))<0.4) return false;
00203     
00204   }
00205   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00206     if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 ||
00207        sqrt(((ntpEvent->vtx.x-1.4828)*(ntpEvent->vtx.x-1.4828)) +
00208             ((ntpEvent->vtx.y-0.2384)*(ntpEvent->vtx.y-0.2384)))>1) return 
00209 false;
00210     //       sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) +
00211     //      ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) return 
00212     //false;
00213   }
00214   
00215   return true;
00216 }

Bool_t MadQuantities::IsFidVtxTrk ( Int_t  itrk = 0  ) 

Definition at line 220 of file MadQuantities.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadTrack(), MadBase::ntpHeader, MadBase::ntpTrack, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by ShowerValidation().

00220                                            {
00221   if(!LoadTrack(itrk)) return false;
00222   
00223   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
00224     
00225     if(ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>29.4 ||   //ends
00226        (ntpTrack->vtx.z<16.5&&ntpTrack->vtx.z>14.5) ||  //between SMs
00227        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00228             +(ntpTrack->vtx.y*ntpTrack->vtx.y))>3.5 ||
00229        sqrt((ntpTrack->vtx.x*ntpTrack->vtx.x)           //radial cut
00230             +(ntpTrack->vtx.y*ntpTrack->vtx.y))<0.4) return false;
00231     
00232   }
00233   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00234     if(ntpTrack->vtx.z<1 || ntpTrack->vtx.z>5 ||
00235        sqrt(((ntpTrack->vtx.x-1.4828)*(ntpTrack->vtx.x-1.4828)) +
00236             ((ntpTrack->vtx.y-0.2384)*(ntpTrack->vtx.y-0.2384)))>1) return false;
00237     //       sqrt(((ntpTrack->vtx.x-1.4885)*(ntpTrack->vtx.x-1.4885)) +
00238     //      ((ntpTrack->vtx.y-0.1397)*(ntpTrack->vtx.y-0.1397)))>1) return false;
00239   }
00240   return true;
00241 }

Bool_t MadQuantities::IsSingleTrack (  ) 

Definition at line 133 of file MadQuantities.cxx.

References MadBase::eventSummary, and NtpSREventSummary::ntrack.

00133                                    {
00134   if(eventSummary->ntrack==1) return true;
00135   return false;
00136 }

Bool_t MadQuantities::IsTrack (  ) 

Definition at line 126 of file MadQuantities.cxx.

References MadBase::eventSummary, and NtpSREventSummary::ntrack.

00126                              {
00127   if(eventSummary->ntrack>0) return true;
00128   return false;
00129 }

void MadQuantities::MakeValidationFile ( std::string  tag  ) 

Definition at line 1816 of file MadQuantities.cxx.

References CCNCSepVars(), NtpSRTrack::ds, NtpSRTrack::end, MadBase::eventSummary, GeantinoEnergy(), MadBase::GetEntry(), GetShwSC(), GetTrueShwSC(), IAction(), INu(), IResonance(), IsFidAll(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadTruthForRecoTH(), MadBase::Nentries, NtpSREventSummary::nevent, MadBase::ntpEvent, MadBase::ntpTrack, NumFSGeantino(), PassCuts(), PassQuasiCuts(), RecoMuEnergy(), RecoMuQP(), RecoQENuEnergy(), RecoShwEnergy(), TrueMuEnergy(), TrueMuQP(), TrueNuEnergy(), TrueShwEnergy(), TrueVtx(), NtpSREvent::vtx, NtpSRTrack::vtx, W2(), NtpSRVertex::x, Y(), NtpSRVertex::y, and NtpSRVertex::z.

01816                                                    {
01817 
01818   std::string savename = "RecoHistos_" + tag + ".root";
01819   TFile save(savename.c_str(),"RECREATE"); 
01820   save.cd();
01821   
01822   //#declare lots of histos:
01823   TH1F *TrueEnergy = new TH1F("TrueEnergy","True Neutrino Energy",200,0,50);
01824 
01825   TH1F *RecoEnergy = new TH1F("RecoEnergy","Reco Neutrino Energy",200,0,50);
01826 
01827   TH1F *TrueMuEn = new TH1F("TrueMuEn","True Muon Energy",200,0,50);
01828 
01829   TH1F *RecoMuEn = new TH1F("RecoMuEn","Reco Muon Energy",200,0,50);
01830   
01831   TH1F *TrueShwEn = new TH1F("TrueShwEn","True Shower Energy",200,0,50);
01832 
01833   TH1F *RecoShwEn = new TH1F("RecoShwEn","Reco Shower Energy",200,0,50);
01834   
01835   TH1F *TrueMuQP = new TH1F("TrueMuQP","True Muon q/p",300,-3,3);
01836 
01837   TH1F *RecoMuQP = new TH1F("RecoMuQP","Reconstructed Muon q/p",300,-3,3);
01838   
01839   TH2F *TrueVsReco = new TH2F("TrueVsReco",
01840                               "Reco Vs True Neutrino Energy",
01841                               200,0,50,200,0,50);
01842 
01843   TH2F *TrueVsRecoMu = new TH2F("TrueVsRecoMu",
01844                                 "Reco Vs True Muon Energy",200,0,50,200,0,50);
01845 
01846   TH2F *TrueVsRecoShw = new TH2F("TrueVsRecoShw",
01847                                  "Reco Vs True Shower Energy",
01848                                  200,0,50,200,0,50);
01849   
01850   TH2F *TrueRecoDiff 
01851     = new TH2F("TrueRecoDiff",
01852                "(True - Reco)/True Neutrino Energy vs True Neutrino Energy",
01853                200,0,50,1000,-9,1);
01854 
01855   TH2F *TrueRecoDiffMuCurv 
01856     = new TH2F("TrueRecoDiffMuCurv",
01857                "(True - Reco_{Curv})/True Muon Energy vs True Muon Energy",
01858                200,0,50,1000,-9,1);
01859   
01860   TH2F *TrueRecoDiffMuRange 
01861     = new TH2F("TrueRecoDiffMuRange",
01862                "(True - Reco_{Range})/True Muon Energy vs True Muon Energy",
01863                200,0,50,1000,-9,1);
01864 
01865   TH2F *TrueRecoDiffMuRangeHardCuts 
01866     = new TH2F("TrueRecoDiffMuRangeHardCuts",
01867                "(True - Reco_{Range})/True Muon Energy vs True Muon Energy with Hard Fiducial Cuts",200,0,50,1000,-9,1);
01868 
01869   TH2F *TrueRecoDiffShw 
01870     = new TH2F("TrueRecoDiffShw",
01871                "(True - Reco)/True Shower Energy vs True Shower Energy",
01872                200,0,50,1000,-9,1);
01873   
01874   TH2F *TrueRangeMu 
01875     = new TH2F("TrueRangeMu",
01876                "Reco Muon Energy from Range vs True Muon Energy",
01877                200,0,50,200,0,50); 
01878 
01879   TH2F *TrueCurvMu 
01880     = new TH2F("TrueCurvMu",
01881                "Reco Muon Energy from Curvature vs True Muon Energy",
01882                200,0,50,200,0,50); 
01883 
01884   TH2F *RecoRangeCurvMu 
01885     =new TH2F("RecoRangeCurvMu",
01886               "Reco Muon Energy from Range vs Reco Muon Energy from Curvature",
01887               200,0,50,200,0,50);
01888   
01889   TH2F *ShwSCVsTrueEn = new TH2F("ShwSCVsTrueEn",
01890                                  "Summed Shower SCs vs True Shower Energy",
01891                                  60,0,30,1000,0,100000);
01892 
01893   TH2F *ShwSCVsTrueEnZoom 
01894     = new TH2F("ShwSCVsTrueEnZoom",
01895                "Zoomed Summed Shower SCs vs True Shower Energy",
01896                50,0,5,1000,0,50000);
01897 
01898   TH1F *ShwSCperGeV = new TH1F("ShwSCperGeV",
01899                                "Summed Shower SCs per GeV Energy (from Truth)",
01900                                1000,0,10000);
01901 
01902   TH2F *RecoShwSCVsTrueShwEn 
01903     = new TH2F("RecoShwSCVsTrueShwEn",
01904                "Reconstructed Shower SCs Vs True Shower Energy",
01905                120,0,30,1000,0,100000);
01906     
01907   //Attempting to estimate neutrino energy in QE events using track angle:
01908   TH1F *TrueQENuEn = new TH1F("TrueQENuEn",
01909                               "True Neutrino Energy for QE-like Events",
01910                               200,0,50);
01911 
01912   TH1F *TrueNuEn_QEBackGround 
01913     = new TH1F("TrueNuEn_QEBackGround",
01914                "True Neutrino Energy for Background QE-like NC Events",
01915                200,0,50);
01916 
01917   TH1F *RecoQENuEn = new TH1F("RecoQENuEn","Reconstructed Neutrino Energy for QE-like Events using Muon Energy and Angle",200,0,50);
01918 
01919   TH1F *RecoQEMuEn = new TH1F("RecoQEMuEn",
01920                               "Reconstructed Muon Energy for QE-like Events",
01921                               200,0,50);
01922 
01923   TH2F *TrueRecoQEDiff 
01924     = new TH2F("TrueRecoQEDiff","(True - Reco_{QE})/True Neutrino Energy vs True Neutrino Energy for QE-like Events",200,0,50,1000,-9,1);
01925 
01926   TH2F *TrueRecoQEDiff_TrueQE 
01927     = new TH2F("TrueRecoQEDiff_TrueQE","(True - Reco_{QE})/True Neutrino Energy vs True Neutrino Energy for True QE Events",200,0,50,1000,-9,1);
01928 
01929   TH2F *TrueNuRecoMuDiff 
01930     = new TH2F("TrueNuRecoMuDiff","(TrueNu - RecoMu_{QE})/TrueNu Energy vs True Neutrino Energy for QE-like Events",200,0,50,1000,-9,1);
01931 
01932   TH2F *TrueNuRecoMuDiff_TrueQE 
01933     = new TH2F("TrueNuRecoMuDiff_TrueQE","(TrueNu - RecoMu_{QE})/TrueNu Energy vs True Neutrino Energy for True QE Events",200,0,50,1000,-9,1);
01934   
01935   //#Efficiency histos:  
01936   TH2F *MuEff = new TH2F("MuEff","MuEff",100,0,100,10,0,1);
01937   TH2F *MuEff_all = new TH2F("MuEff_all","MuEff_all",100,0,100,10,0,1);
01938   TH2F *MuEff_fid_all = new TH2F("MuEff_fid_all",
01939                                  "MuEff_fid_all",100,0,100,10,0,1);
01940 
01941   TH2F *NCMisEff = new TH2F("NCMisEff","NCMisEff",100,0,100,10,0,1);
01942   TH2F *NCMisEff_all = new TH2F("NCMisEff_all",
01943                                 "NCMisEff_all",100,0,100,10,0,1);
01944   TH2F *NCMisEff_fid_all = new TH2F("NCMisEff_fid_all",
01945                                     "NCMisEff_fid_all",100,0,100,10,0,1);
01946   
01947   //#Pi background histos:
01948   TH1F *TrkLenCC = new TH1F("TrkLenCC","Track length - NuMu CC",50,0,50);
01949   TH1F *TrkLenNC = new TH1F("TrkLenNC","Track length - NC",50,0,50);
01950   
01951   TH1F *dedxCC = new TH1F("dedxCC","Average dEdx - NuMu CC",1000,0,2000);
01952   TH1F *dedxNC = new TH1F("dedxNC","Average dEdx - NC",1000,0,2000);
01953 
01954   TH1F *HalfRatioCC 
01955     = new TH1F("HalfRatioCC",
01956                "Charge Ratio: First Half/Second Half of Track - NuMu CC",
01957                150,-1,14);
01958   TH1F *HalfRatioNC 
01959     = new TH1F("HalfRatioNC",
01960                "Charge Ratio: First Half/Second Half of Track - NuMu NC",
01961                150,-1,14);
01962   
01963   TH2F *MuRangeCurvDiffVsTrkLen 
01964       = new TH2F("MuRangeCurvDiffVsTrkLen",
01965                  "Diff in Momentum from Range and Curv Vs Track Length (CC)",
01966                  100,0,30,200,-3,17);
01967   TH2F *PiRangeCurvDiffVsTrkLen 
01968     = new TH2F("PiRangeCurvDiffVsTrkLen",
01969                "Diff in Momentum from Range and Curv Vs Track Length (NC)",
01970                100,0,30,200,-3,17);
01971   
01972   //#quasi plots
01973   TH2F *TrueShwSCVsQuasiNuEn
01974     = new TH2F("TrueShwSCVsQuasiNuEn",
01975                "True Shower SCs Vs Neutrino Energy for True QE events",
01976                60,0,30,500,0,10000);
01977   
01978   TH2F *TrueShwSCVsNonQENuEn
01979     = new TH2F("TrueShwSCVsNonQENuEn",
01980                "True Shower SCs Vs Neutrino Energy for True Non-QE events",
01981                60,0,30,500,0,10000);
01982   
01983   TH1F *TrueQENuEn_Pass 
01984     = new TH1F("TrueQENuEn_Pass",
01985                "True Neutrino Energy for QE Events that pass Reco Cuts",
01986                200,0,50);
01987   
01988   TH1F *SRMCVtxDistHist
01989     = new TH1F("SRMCVtxDistHist",
01990                "Distance between SR and MC Event Vertex",
01991                1000,0,5);
01992 
01993   TH1F *geantinoEnergy_CC[3];
01994   geantinoEnergy_CC[0] = new TH1F("geantinoEnergy_QECC",
01995                                   "Geantino Energy in Event",1000,0,5);
01996   geantinoEnergy_CC[1] = new TH1F("geantinoEnergy_RESCC",
01997                                   "Geantino Energy in Event",1000,0,5);
01998   geantinoEnergy_CC[2] = new TH1F("geantinoEnergy_DISCC",
01999                                   "Geantino Energy in Event",1000,0,5);
02000 
02001   TH1F *geantinoEnergy_NC[3];
02002   geantinoEnergy_NC[0] = new TH1F("geantinoEnergy_QENC",
02003                                   "Geantino Energy in Event",1000,0,5);
02004   geantinoEnergy_NC[1] = new TH1F("geantinoEnergy_RESNC",
02005                                   "Geantino Energy in Event",1000,0,5);
02006   geantinoEnergy_NC[2] = new TH1F("geantinoEnergy_DISNC",
02007                                   "Geantino Energy in Event",1000,0,5);
02008 
02009   TH1F *fracGeantinoEnergy_CC[3];
02010   fracGeantinoEnergy_CC[0] = new TH1F("fracGeantinoEnergy_QECC",
02011                                       "Fractional Geantino Energy in Event",
02012                                       1000,0,1);
02013   fracGeantinoEnergy_CC[1] = new TH1F("fracGeantinoEnergy_RESCC",
02014                                       "Fractional Geantino Energy in Event",
02015                                       1000,0,1);
02016   fracGeantinoEnergy_CC[2] = new TH1F("fracGeantinoEnergy_DISCC",
02017                                       "Fractional Geantino Energy in Event",
02018                                       1000,0,1);
02019 
02020   TH1F *fracGeantinoEnergy_NC[3];
02021   fracGeantinoEnergy_NC[0] = new TH1F("fracGeantinoEnergy_QENC",
02022                                       "Fractional Geantino Energy in Event",
02023                                       1000,0,1);
02024   fracGeantinoEnergy_NC[1] = new TH1F("fracGeantinoEnergy_RESNC",
02025                                       "Fractional Geantino Energy in Event",
02026                                       1000,0,1);
02027   fracGeantinoEnergy_NC[2] = new TH1F("fracGeantinoEnergy_DISNC",
02028                                       "Fractional Geantino Energy in Event",
02029                                       1000,0,1);
02030 
02031   TH1F *numGeantino_CC[3];
02032   numGeantino_CC[0] = new TH1F("numGeantino_QECC",
02033                                "Number of Geantino in Event",10,0,10);
02034   numGeantino_CC[1] = new TH1F("numGeantino_RESCC",
02035                                "Number of Geantino in Event",10,0,10);
02036   numGeantino_CC[2] = new TH1F("numGeantino_DISCC",
02037                                "Number of Geantino in Event",10,0,10);
02038 
02039   TH1F *numGeantino_NC[3];
02040   numGeantino_NC[0] = new TH1F("numGeantino_QENC",
02041                                "Number of Geantino in Event",10,0,10);
02042   numGeantino_NC[1] = new TH1F("numGeantino_RESNC",
02043                                "Number of Geantino in Event",10,0,10);
02044   numGeantino_NC[2] = new TH1F("numGeantino_DISNC",
02045                                "Number of Geantino in Event",10,0,10);
02046 
02047 
02048   for(int i=0;i<Nentries;i++){
02049     
02050     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
02051                             << "% done" << std::endl;
02052     
02053     if(!GetEntry(i)) continue;
02054 
02055     for(int j=0;j<eventSummary->nevent;j++){
02056       
02057       if(!LoadEvent(j)) continue;
02058       int mcevent = -1;
02059       if(!LoadTruthForRecoTH(j,mcevent)) continue;
02060 
02061       if(IAction(mcevent)==1) {
02062         Int_t index = IResonance(mcevent)-1001;
02063         if(index==3) index = 2;
02064         if(NumFSGeantino(mcevent)>0){
02065           geantinoEnergy_CC[index]->Fill(GeantinoEnergy(mcevent));
02066           fracGeantinoEnergy_CC[index]->Fill(GeantinoEnergy(mcevent) / 
02067                                              TrueNuEnergy(mcevent));
02068         }       
02069         numGeantino_CC[index]->Fill(NumFSGeantino(mcevent));
02070       }
02071       else {
02072         Int_t index = IResonance(mcevent)-1001;
02073         if(index==3) index = 2;
02074         if(NumFSGeantino(mcevent)>0){
02075           geantinoEnergy_NC[index]->Fill(GeantinoEnergy(mcevent));
02076           fracGeantinoEnergy_NC[index]->Fill(GeantinoEnergy(mcevent) / 
02077                                              TrueNuEnergy(mcevent));
02078         }
02079         numGeantino_NC[index]->Fill(NumFSGeantino(mcevent));
02080       }
02081       
02082       int track = -1;
02083       LoadLargestTrackFromEvent(j,track);
02084       int shower = -1;
02085       LoadLargestShowerFromEvent(j,shower);
02086       
02087       Float_t *vtx = this->TrueVtx(mcevent);
02088       float dist = sqrt(TMath::Power(vtx[0]-ntpEvent->vtx.x,2)+
02089                         TMath::Power(vtx[1]-ntpEvent->vtx.y,2)+
02090                         TMath::Power(vtx[2]-ntpEvent->vtx.z,2));
02091       SRMCVtxDistHist->Fill(dist);
02092       
02093       if(abs(this->INu(mcevent))==14&&this->IAction(mcevent)==1){
02094         
02095         MuEff_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02096         
02097         if(sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1])<=3.5
02098            &&vtx[2]>=0.5){
02099           MuEff_fid_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02100           if(this->PassCuts(track)) 
02101             MuEff->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02102         }
02103         
02104         if(sqrt(this->W2(mcevent))<1.0) {
02105           TrueShwSCVsQuasiNuEn->Fill(this->TrueNuEnergy(mcevent),
02106                                      this->GetTrueShwSC(mcevent));
02107         }
02108         else TrueShwSCVsNonQENuEn->Fill(this->TrueNuEnergy(mcevent),
02109                                         this->GetTrueShwSC(mcevent));
02110         
02111         if(this->PassCuts(track)){
02112           //neutrino
02113           TrueEnergy->Fill(this->TrueNuEnergy(mcevent));
02114           RecoEnergy->Fill(this->RecoMuEnergy(1,track)
02115                            +this->RecoShwEnergy(shower));
02116           TrueVsReco->Fill(this->TrueNuEnergy(mcevent),
02117                            this->RecoMuEnergy(1,track) 
02118                            + this->RecoShwEnergy(shower));
02119           TrueRecoDiff->Fill(this->TrueNuEnergy(mcevent),
02120                              (this->TrueNuEnergy(mcevent)
02121                               -(this->RecoMuEnergy(1,track)
02122                                 +this->RecoShwEnergy(shower)))
02123                              /this->TrueNuEnergy(mcevent));
02124           
02125           //shower        
02126           TrueShwEn->Fill(this->TrueShwEnergy(mcevent));
02127           RecoShwEn->Fill(this->RecoShwEnergy(shower));
02128           TrueVsRecoShw->Fill(this->TrueShwEnergy(mcevent),
02129                               this->RecoShwEnergy(shower));
02130           TrueRecoDiffShw->Fill(this->TrueShwEnergy(mcevent),
02131                                 (this->TrueShwEnergy(mcevent)
02132                                  -this->RecoShwEnergy(shower))
02133                                 /this->TrueShwEnergy(mcevent));
02134           
02135           ShwSCVsTrueEn->Fill(this->TrueShwEnergy(mcevent),
02136                               this->GetTrueShwSC(mcevent));
02137           ShwSCVsTrueEnZoom->Fill(this->TrueShwEnergy(mcevent),
02138                                   this->GetTrueShwSC(mcevent));
02139           ShwSCperGeV->Fill(this->GetTrueShwSC(mcevent)
02140                             /this->TrueShwEnergy(mcevent));
02141           RecoShwSCVsTrueShwEn->Fill(this->TrueShwEnergy(mcevent),
02142                                      this->GetShwSC(track));
02143           
02144           //muon
02145           TrueMuEn->Fill(this->TrueMuEnergy(mcevent));
02146           RecoMuEn->Fill(this->RecoMuEnergy(0,track));
02147           TrueVsRecoMu->Fill(this->TrueMuEnergy(mcevent),
02148                              this->RecoMuEnergy(1,track));
02149             
02150           //curvature
02151           TrueRecoDiffMuCurv->Fill(this->TrueMuEnergy(mcevent),
02152                                    (this->TrueMuEnergy(mcevent)
02153                                     -this->RecoMuEnergy(1,track))
02154                                    /this->TrueMuEnergy(mcevent));
02155           TrueCurvMu->Fill(this->TrueMuEnergy(mcevent),
02156                            this->RecoMuEnergy(1,track));
02157           
02158           TrueMuQP->Fill(this->TrueMuQP(mcevent));
02159           RecoMuQP->Fill(this->RecoMuQP(track));
02160           
02161           //range
02162           if(this->IsFidAll(track)){
02163             TrueRecoDiffMuRange->Fill(this->TrueMuEnergy(mcevent),
02164                                       (this->TrueMuEnergy(mcevent)
02165                                        -this->RecoMuEnergy(2,track))
02166                                       /this->TrueMuEnergy(mcevent));      
02167               
02168             if(sqrt(TMath::Power(ntpTrack->vtx.x,2)
02169                     +TMath::Power(ntpTrack->vtx.y,2))>0.5
02170                &&sqrt(TMath::Power(ntpTrack->end.x,2)
02171                       +TMath::Power(ntpTrack->end.y,2))>0.5
02172                &&ntpTrack->end.z<25.) {
02173               TrueRecoDiffMuRangeHardCuts->Fill(this->TrueMuEnergy(mcevent),
02174                                                 (this->TrueMuEnergy(mcevent)
02175                                                  -this->RecoMuEnergy(2,track))
02176                                                 /this->TrueMuEnergy(mcevent));    
02177             }
02178             
02179             TrueRangeMu->Fill(this->TrueMuEnergy(mcevent),
02180                               this->RecoMuEnergy(2,track));
02181             
02182             RecoRangeCurvMu->Fill(this->RecoMuEnergy(1,track),
02183                                   this->RecoMuEnergy(2,track));
02184             MuRangeCurvDiffVsTrkLen->Fill(ntpTrack->ds,
02185                                           (this->RecoMuEnergy(1,track)
02186                                            -this->RecoMuEnergy(2,track))
02187                                           /this->RecoMuEnergy(1,track));
02188           }
02189           
02190           //QE-like
02191           if(this->PassQuasiCuts(track)){
02192             RecoQENuEn->Fill(this->RecoQENuEnergy(track));
02193             RecoQEMuEn->Fill(this->RecoMuEnergy(0,track));
02194             TrueQENuEn->Fill(this->TrueNuEnergy(mcevent));
02195             TrueRecoQEDiff->Fill(this->TrueNuEnergy(mcevent),
02196                                  (this->TrueNuEnergy(mcevent)
02197                                   -this->RecoQENuEnergy(track))
02198                                  /this->TrueNuEnergy(mcevent));
02199             TrueNuRecoMuDiff->Fill(this->TrueNuEnergy(mcevent),
02200                                    (this->TrueNuEnergy(mcevent)
02201                                     -this->RecoMuEnergy(0,track))
02202                                    /this->TrueNuEnergy(mcevent));
02203             if(this->W2(mcevent)<1.0) {
02204               TrueQENuEn_Pass->Fill(this->TrueNuEnergy(mcevent));
02205               TrueRecoQEDiff_TrueQE->Fill(this->TrueNuEnergy(mcevent),
02206                                           (this->TrueNuEnergy(mcevent)
02207                                            -this->RecoQENuEnergy(track))
02208                                           /this->TrueNuEnergy(mcevent));
02209               TrueNuRecoMuDiff_TrueQE->Fill(this->TrueNuEnergy(mcevent),
02210                                             (this->TrueNuEnergy(mcevent)
02211                                              -this->RecoMuEnergy(0,track))
02212                                             /this->TrueNuEnergy(mcevent));
02213             }
02214           }
02215           
02216           Float_t *CCNCVars = this->CCNCSepVars(track);
02217           TrkLenCC->Fill(CCNCVars[0]);
02218           dedxCC->Fill(CCNCVars[1]/100.);
02219           HalfRatioCC->Fill(CCNCVars[2]);
02220           delete CCNCVars;
02221         }
02222       }
02223       delete vtx;
02224       if(this->IAction(mcevent)==0){
02225         
02226         NCMisEff_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02227         
02228         Float_t *vtx = this->TrueVtx(mcevent);
02229         if(sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1])<=3.5
02230            &&vtx[2]>=0.5){
02231           NCMisEff_fid_all->Fill(this->TrueNuEnergy(mcevent),this->Y(mcevent));
02232           if(this->PassCuts(0)) NCMisEff->Fill(this->TrueNuEnergy(mcevent),
02233                                                this->Y(mcevent));
02234         }
02235         delete vtx;
02236         
02237         if(this->PassCuts(track)) {
02238           Float_t *CCNCVars = this->CCNCSepVars(track);
02239           TrkLenNC->Fill(CCNCVars[0]);
02240           dedxNC->Fill(CCNCVars[1]/100.);
02241           HalfRatioNC->Fill(CCNCVars[2]);
02242           delete CCNCVars;
02243           if(this->IsFidAll(track)){
02244             PiRangeCurvDiffVsTrkLen->Fill(ntpTrack->ds,
02245                                           (this->RecoMuEnergy(1,track)
02246                                            -this->RecoMuEnergy(2,track))
02247                                           /this->RecoMuEnergy(1,track));
02248           }
02249           
02250           if(this->PassQuasiCuts(track)){
02251             TrueNuEn_QEBackGround->Fill(this->TrueNuEnergy(mcevent));
02252           }
02253         }
02254       }
02255     }
02256   }
02257 
02258   //#close file
02259   save.Write();
02260   save.Close();
02261  
02262 }

Float_t MadQuantities::MaxFSKaonEnergy ( Int_t  itr = 0  ) 

Definition at line 719 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fMaxFSKaonEnergy, MadBase::lastEntry, and SetStdHepVars().

00719                                                {
00720   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00721     SetStdHepVars(itr);
00722   return fMaxFSKaonEnergy;
00723 }

Float_t MadQuantities::MaxFSNeutronEnergy ( Int_t  itr = 0  ) 

Definition at line 676 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fMaxFSNeutronEnergy, MadBase::lastEntry, and SetStdHepVars().

00676                                                   {
00677   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00678     SetStdHepVars(itr);
00679   return fMaxFSNeutronEnergy;
00680 }

Float_t MadQuantities::MaxFSPionEnergy ( Int_t  itr = 0,
Int_t  pm = 0 
)

Definition at line 698 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fMaxFSPiMinusEnergy, fMaxFSPiPlusEnergy, MadBase::lastEntry, and SetStdHepVars().

00698                                                         {
00699   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00700     SetStdHepVars(itr);
00701   if(pm==-1) return fMaxFSPiMinusEnergy;
00702   else if(pm==1) return fMaxFSPiPlusEnergy;  
00703   if(fMaxFSPiPlusEnergy>fMaxFSPiMinusEnergy) return fMaxFSPiPlusEnergy;
00704   return fMaxFSPiMinusEnergy;
00705 }

Float_t MadQuantities::MaxFSPiZeroEnergy ( Int_t  itr = 0  ) 

Definition at line 737 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fMaxFSPiZeroEnergy, MadBase::lastEntry, and SetStdHepVars().

00737                                                  {
00738   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00739     SetStdHepVars(itr);
00740   return fMaxFSPiZeroEnergy;
00741 }

Float_t MadQuantities::MaxFSProtonEnergy ( Int_t  itr = 0  ) 

Definition at line 658 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fMaxFSProtonEnergy, MadBase::lastEntry, and SetStdHepVars().

00658                                                  {
00659   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00660     SetStdHepVars(itr);
00661   return fMaxFSProtonEnergy;
00662 }

Int_t MadQuantities::Nucleus ( Int_t  itr = 0  ) 

Definition at line 773 of file MadQuantities.cxx.

References MuELoss::a, NtpMCTruth::a, MadBase::LoadTruth(), MadBase::ntpTruth, and NtpMCTruth::z.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::MyCreatePAN(), and MadAnalysis::RecoMC().

00773                                      {
00774  
00775   if(!LoadTruth(itr)) return 0;
00776 
00777   Int_t z=int(ntpTruth->z);
00778   Int_t a=int(ntpTruth->a);
00779   Int_t nucleus=1;
00780 
00781   switch (z) {
00782     //  case 1:
00783     //nucleus=0;   // free nucleon
00784     //break;
00785   case 1:
00786     switch (a) {
00787     case 1:
00788       nucleus=256;   // hydrogen1
00789       break;
00790     case 2:
00791       nucleus=257;   // hydrogen2
00792       break;
00793     case 3:
00794       nucleus=258;   // hydrogen2
00795       break;
00796     default:
00797       nucleus=256;   // hydrogen1
00798       break;
00799     }
00800     break;
00801   case 6:
00802     switch (a) {
00803     case 11:
00804       nucleus=274; // carbon11   
00805       break;
00806     case 12:
00807       nucleus=275; // carbon12
00808       break;
00809     case 13:
00810       nucleus=276; // carbon13
00811       break;
00812     case 14:
00813       nucleus=277; // carbon14
00814       break;
00815     default:
00816       nucleus=275; // carbon12
00817       break;
00818     }
00819     break;
00820   case 7:
00821     switch (a) {
00822     case 13:
00823       nucleus=278; // nitrogen13   
00824       break;
00825     case 14:
00826       nucleus=279; // nitrogen14
00827       break;
00828     case 15:
00829       nucleus=280; // nitrogen15
00830       break;
00831     case 16:
00832       nucleus=281; // nitrogen16
00833       break;
00834     case 17:
00835       nucleus=282; // nitrogen17
00836       break;
00837     default:
00838       nucleus=279; // nitrogen14
00839       break;
00840     }
00841     break;
00842   case 8:
00843     switch (a) {
00844     case 15:
00845       nucleus=283;   // oxygen15
00846       break;
00847     case 16:
00848       nucleus=284;   // oxygen16
00849       break;
00850     case 17:
00851       nucleus=285;   // oxygen17
00852       break;
00853     case 18:
00854       nucleus=286;   // oxygen18
00855       break;
00856     default:
00857       nucleus=284;   // oxygen16
00858       break;
00859     }
00860     break;
00861   case 13:
00862     switch (a) {
00863     case 26:
00864       nucleus=303;   // aluminium26
00865       break;
00866     case 27:
00867       nucleus=304;   // aluminium27
00868       break;
00869     case 28:
00870       nucleus=305;   // aluminium28
00871       break;
00872     case 29:
00873       nucleus=306;   // aluminium29
00874       break;
00875     default:
00876       nucleus=304;   // aluminium27
00877       break;
00878     }
00879     break;
00880   case 14:
00881     switch (a) {
00882     case 27:
00883       nucleus=307;   // silicon27
00884       break;
00885     case 28:
00886       nucleus=308;   // silicon28
00887       break;
00888     case 29:
00889       nucleus=309;   // silicon29
00890       break;
00891     case 30:
00892       nucleus=310;   // silicon30
00893       break;
00894     default:         
00895       nucleus=308;   // silicon28
00896       break;
00897     }
00898     break;
00899   case 15:
00900     switch (a) {
00901     case 30:
00902       nucleus=311;   //phosphorus30
00903       break;
00904     case 31:
00905       nucleus=312;   //phosphorus31
00906       break;
00907     case 32:
00908       nucleus=313;   //phosphorus32
00909       break;
00910     case 33:
00911       nucleus=314;   //phosphorus33
00912       break;
00913     default:
00914       nucleus=312;   //phosphorus31
00915       break;
00916     }
00917     break;
00918   case 16:
00919     switch (a) {
00920     case 31:
00921       nucleus=315;   //sulphur31
00922       break;
00923     case 32:
00924       nucleus=316;   //sulphur32
00925       break;
00926     case 33:
00927       nucleus=317;   //sulphur33
00928       break;
00929     case 34:
00930       nucleus=318;   //sulphur34
00931       break;
00932     case 35:
00933       nucleus=319;   //sulphur35
00934       break;
00935     case 36:
00936       nucleus=320;   //sulphur36
00937       break;
00938     default:
00939       nucleus=316;   //sulphur32
00940       break;
00941     }
00942     break;
00943   case 22:
00944     switch (a) {
00945     case 45:
00946       nucleus=347;   //titanium45
00947       break;
00948     case 46:
00949       nucleus=348;   //titanium46
00950       break;
00951     case 47:
00952       nucleus=349;   //titanium47
00953       break;
00954     case 48:
00955       nucleus=350;   //titanium48
00956       break;
00957     case 49:
00958       nucleus=351;   //titanium49
00959       break;
00960     case 50:
00961       nucleus=352;   //titanium50
00962       break;
00963     default:
00964       nucleus=350;   //titanium48
00965       break;
00966     }
00967     break;
00968   case 23:
00969     switch (a) {
00970     case 49:
00971       nucleus=353;   //vanadium49
00972       break;
00973     case 50:
00974       nucleus=354;   //vanadium50
00975       break;
00976     case 51:
00977       nucleus=355;   //vanadium51
00978       break;
00979     case 52:
00980       nucleus=356;   //vanadium52
00981       break;
00982     case 53:
00983       nucleus=357;   //vanadium53
00984       break;
00985     default:
00986       nucleus=355;   //vanadium51
00987       break;
00988     }
00989     break;
00990   case 24:
00991     switch (a) {
00992     case 49:
00993       nucleus=358;   //chromium49
00994       break;
00995     case 50:
00996       nucleus=359;   //chromium50
00997       break;
00998     case 51:
00999       nucleus=360;   //chromium51
01000       break;
01001     case 52:
01002       nucleus=361;   //chromium52
01003       break;
01004     case 53:
01005       nucleus=362;   //chromium53
01006       break;
01007     case 54:
01008       nucleus=363;   //chromium54
01009       break;
01010     default:
01011       nucleus=361;   //chromium52
01012       break;
01013     }
01014     break;
01015   case 25:
01016     switch (a) {
01017     case 53:
01018       nucleus=364;   //manganese53
01019       break;
01020     case 54:
01021       nucleus=365;   //manganese54
01022       break;    
01023     case 55:
01024       nucleus=366;   //manganese55
01025       break;    
01026     case 56:
01027       nucleus=367;   //manganese56
01028       break;    
01029     case 57:
01030       nucleus=368;   //manganese57
01031       break;    
01032     default:
01033       nucleus=366;   //manganese55
01034       break;
01035     }
01036     break;
01037   case 26:
01038     switch (a) {
01039     case 53:
01040       nucleus=369;   //iron53
01041       break;
01042     case 54:
01043       nucleus=370;   //iron54
01044       break;
01045     case 55:
01046       nucleus=371;   //iron55
01047       break;
01048     case 56:
01049       nucleus=372;   //iron56
01050       break;
01051     case 57:
01052       nucleus=373;   //iron57
01053       break;
01054     case 58:
01055       nucleus=374;   //iron58
01056       break;
01057     default:
01058       nucleus=372;   //iron56
01059       break;
01060     }
01061     break;
01062   case 28:
01063     switch (a) {
01064     case 57:
01065       nucleus=382;   //nickel57
01066       break;
01067     case 58:
01068       nucleus=383;   //nickel58
01069       break;
01070     case 59:
01071       nucleus=384;   //nickel59
01072       break;
01073     case 60:
01074       nucleus=385;   //nickel60
01075       break;
01076     case 61:
01077       nucleus=386;   //nickel61
01078       break;
01079     case 62:
01080       nucleus=387;   //nickel62
01081       break;
01082     case 63:
01083       nucleus=388;   //nickel63
01084       break;
01085     case 64:
01086       nucleus=389;   //nickel64
01087       break;
01088     default:
01089       nucleus=383;   //nickel58
01090       break;
01091     }
01092     break;
01093   case 29:
01094     switch (a) {
01095     case 62:
01096       nucleus=390;   //copper62
01097       break;
01098     case 63:
01099       nucleus=391;   //copper63
01100       break;
01101     case 64:
01102       nucleus=392;   //copper64
01103       break;
01104     case 65:
01105       nucleus=393;   //copper65
01106       break;
01107     case 66:
01108       nucleus=394;   //copper66
01109       break;
01110     case 67:
01111       nucleus=395;   //copper67
01112       break;
01113     default:
01114       nucleus=392;   //copper64
01115       break;
01116     }
01117     break;
01118   default:
01119     nucleus=1;   // unknown
01120     break;
01121   }
01122 
01123   return nucleus;
01124 }

Int_t MadQuantities::NumFSGeantino ( Int_t  itr = 0  ) 

Definition at line 640 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSGeantino, MadBase::lastEntry, and SetStdHepVars().

Referenced by MakeValidationFile().

00640                                            {
00641   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00642     SetStdHepVars(itr);
00643   return fNumFSGeantino;
00644 }

Int_t MadQuantities::NumFSKaon ( Int_t  itr = 0  ) 

Definition at line 707 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSKaon, MadBase::lastEntry, and SetStdHepVars().

00707                                        {
00708   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00709     SetStdHepVars(itr);
00710   return fNumFSKaon;
00711 }

Int_t MadQuantities::NumFSNeutron ( Int_t  itr = 0  ) 

Definition at line 664 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSNeutron, MadBase::lastEntry, and SetStdHepVars().

Referenced by MadMKAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00664                                           {
00665   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00666     SetStdHepVars(itr);
00667   return fNumFSNeutron;
00668 }

Int_t MadQuantities::NumFSPion ( Int_t  itr = 0,
Int_t  pm = 0 
)

Definition at line 682 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSPiMinus, fNumFSPiPlus, MadBase::lastEntry, and SetStdHepVars().

Referenced by MadMKAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00682                                                 {
00683   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00684     SetStdHepVars(itr);
00685   if(pm==-1) return fNumFSPiMinus;
00686   else if(pm==1) return fNumFSPiPlus;
00687   return fNumFSPiPlus+fNumFSPiMinus;
00688 }

Int_t MadQuantities::NumFSPiZero ( Int_t  itr = 0  ) 

Definition at line 725 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSPiZero, MadBase::lastEntry, and SetStdHepVars().

Referenced by MadMKAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00725                                          {
00726   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00727     SetStdHepVars(itr);
00728   return fNumFSPiZero;
00729 }

Int_t MadQuantities::NumFSProton ( Int_t  itr = 0  ) 

Definition at line 646 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fNumFSProton, MadBase::lastEntry, and SetStdHepVars().

Referenced by MadMKAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00646                                          {
00647   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00648     SetStdHepVars(itr);
00649   return fNumFSProton;
00650 }

Bool_t MadQuantities::PassCuts (  ) 

Definition at line 158 of file MadQuantities.cxx.

References TrueMuEnergy().

Referenced by MadEdAnalysis::MakeEff(), MakeValidationFile(), and PassQuasiCuts().

00158                               {  //specifically for event display
00159   //demand PassCuts for the first track
00160   if(!PassCuts(0)) return false;
00161   //make any cuts you like to define events to be displayed:
00162   if(TrueMuEnergy(0)==0) return false;
00163   return true;
00164 }

Bool_t MadQuantities::PassCuts ( Int_t  itrk  ) 

Definition at line 150 of file MadQuantities.cxx.

References IsFidVtx(), and PassTrackCuts().

00150                                         {
00151   //check basic fit track quality + check for contained vtx:
00152   if(!PassTrackCuts(itrk)) return false;
00153   if(!IsFidVtx(itrk)) return false;
00154   return true;
00155 }

Bool_t MadQuantities::PassQuasiCuts ( Int_t  itrk = 0,
Float_t  frac = 0.1 
)

Definition at line 484 of file MadQuantities.cxx.

References PassCuts(), RecoMuEnergy(), and RecoShwEnergy().

Referenced by MakeValidationFile().

00484                                                           {
00485   if(!PassCuts(itrk)) return false;
00486   if(RecoShwEnergy()>frac*(RecoShwEnergy()+RecoMuEnergy(0,itrk))) return false;
00487   return true;
00488 }

Bool_t MadQuantities::PassTrackCuts ( Int_t  itrk  ) 
Float_t MadQuantities::Q2 ( Int_t  itr = 0  ) 
Double_t MadQuantities::RecoEMFrac ( Int_t  shower,
Int_t  method,
Double_t &  EMFrac0,
Double_t &  EMFrac1 
)

Definition at line 1681 of file MadQuantities.cxx.

References NtpSRShower::clu, NtpSRStripPulseHeight::gev, NtpSRCluster::id, MadBase::LoadCluster(), MadBase::LoadShower(), NtpSRShower::ncluster, MadBase::ntpCluster, MadBase::ntpShower, NtpSRCluster::ph, NtpSRCluster::planeview, and NtpSRCluster::probem.

Referenced by MadAnalysis::CreatePAN(), and MadCluAnalysis::DataDistributions().

01683 {
01684 
01685   if(!LoadShower(shower)) return 0;
01686   
01687   float EMFrac = 0;
01688   float EMFrac_end0 = 0;
01689   float EMFrac_end1 = 0;
01690   float EnTot = 0;
01691   float EnTot_end0 = 0;
01692   float EnTot_end1 = 0;
01693   for(int k=0; k<ntpShower->ncluster; k++){
01694     int index = ntpShower->clu[k];
01695     if(!LoadCluster(index)) continue;
01696     if(ntpCluster->id==0 && 
01697        (method==0 || ntpCluster->probem>0.2)){
01698       EMFrac+=ntpCluster->ph.gev;
01699       EnTot+=ntpCluster->ph.gev;
01700       if(ntpCluster->planeview==2){
01701         EMFrac_end0+=ntpCluster->ph.gev;
01702         EnTot_end0+=ntpCluster->ph.gev;
01703       }
01704       if(ntpCluster->planeview==3){
01705         EMFrac_end1+=ntpCluster->ph.gev;
01706         EnTot_end1+=ntpCluster->ph.gev;
01707       }
01708     }
01709     else {
01710       EnTot+=ntpCluster->ph.gev;
01711       if(ntpCluster->planeview==2){
01712         EnTot_end0+=ntpCluster->ph.gev;
01713       }
01714       if(ntpCluster->planeview==3){
01715         EnTot_end1+=ntpCluster->ph.gev;
01716       }
01717     }
01718   }
01719 
01720   if(EnTot==0) return 0;
01721   if(EnTot_end0==0) EMFrac = EMFrac_end1/EnTot_end1;
01722   else if(EnTot_end1==0) EMFrac = EMFrac_end0/EnTot_end0;
01723   else if(EMFrac_end0/EnTot_end0>EMFrac_end1/EnTot_end1) 
01724     EMFrac = EMFrac_end0/EnTot_end0;
01725   else EMFrac = EMFrac_end1/EnTot_end1;
01726   
01727   if(EnTot_end0>0) EMFrac0 = EMFrac_end0/EnTot_end0;
01728   if(EnTot_end1>1) EMFrac1 = EMFrac_end1/EnTot_end1;
01729   
01730   return EMFrac;
01731 
01732 }

Float_t MadQuantities::RecoMuDCosNeu ( Int_t  itr = 0  ) 

Definition at line 1422 of file MadQuantities.cxx.

References NtpSRVertex::dcosy, NtpSRVertex::dcosz, MadBase::LoadTrack(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::DataHist(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyCreatePANData(), RecoQENuEnergy(), and RecoQEQ2().

01422                                               { //ds_mu/ds_neu
01423   if(!LoadTrack(itrk)) return 0.;
01424   
01425   Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
01426   Float_t bl_y = sqrt(1. - bl_z*bl_z);
01427   Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
01428    
01429   return  costhbl;
01430 }

Float_t MadQuantities::RecoMuDCosZ ( Int_t  itr = 0  ) 
Float_t MadQuantities::RecoMuEnergy ( Int_t  opt = 0,
Int_t  itrk = 0,
Bool_t  data = true,
Detector::Detector_t  det = Detector::kNear 
)

Definition at line 1330 of file MadQuantities.cxx.

References EnergyCorrections::CorrectMomentumFromRange(), EnergyCorrections::CorrectSignedMomentumFromCurvature(), IsFidAll(), MadBase::LoadTrack(), NtpSRTrack::momentum, MadBase::ntpTrack, NtpSRMomentum::qp, and NtpSRMomentum::range.

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadCluAnalysis::DataDistributions(), MadEdAnalysis::DataHist(), MadCBSQEAnalysis::MakeQEFile(), MakeValidationFile(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyCreatePANData(), MadEdAnalysis::MyMakeQEFile(), MadUserAnalysis::PassBasicCuts(), PassQuasiCuts(), MadTestAnalysis::RecoAnalysisEnergy(), MadCluAnalysis::RecoAnalysisEnergy(), MadCBSQEAnalysis::RecoAnalysisEnergy(), MadPIDAnalysis::RecoAnalysisEnergy(), MadDpAnalysis::RecoAnalysisEnergy(), MadUserAnalysis::RecoAnalysisEnergy(), RecoQENuEnergy(), RecoQEQ2(), ShowerValidation(), and MadCBSQEAnalysis::TestQEDiscrim().

01330                                                                                             {
01331   const float mumass=0.10566;
01332 
01333   if(LoadTrack(itrk)){
01334     Float_t mr=ntpTrack->momentum.range;
01335     if(mr>0) mr=CorrectMomentumFromRange(mr,isdata,det);
01336     float mc=0.0; // signed momentum from curvature
01337     if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
01338     mc=CorrectSignedMomentumFromCurvature(mc,isdata,det);
01339     
01340     if(opt==0){  //return the most appropriate measure of momentum
01341       if(IsFidAll(itrk) || ntpTrack->momentum.qp==0) {
01342         return sqrt(mr*mr+ mumass*mumass);
01343       }
01344       else {
01345         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01346         else return sqrt(mc*mc+ mumass*mumass);
01347       }
01348     }
01349     else if(opt==1) {
01350       //return curvature measurement
01351       // in R1.9 the tracker will apparently return (q/p)=0.0
01352       // maybe it's when a track looks perfectly rigid?
01353       // if so, we have to do something
01354       // I don't want to use the range, that could be very wrong
01355       // so, we'll return an obviously ridiculous value of 10 TeV
01356       if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01357       else return sqrt(mc*mc+ mumass*mumass);
01358     }
01359     else if(opt==2) //return range measurement
01360       return sqrt(mr*mr+ mumass*mumass);
01361     else return 0;
01362   }
01363   return 0.;
01364 }

Float_t MadQuantities::RecoMuEnergyNew ( Int_t  opt,
Int_t  itrk,
VldContext  cx,
ReleaseType::Release_t  reltype,
EnergyCorrections::WhichCorrection_t  corrver 
)

Definition at line 1368 of file MadQuantities.cxx.

References EnergyCorrections::FullyCorrectMomentumFromRange(), EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(), IsFidAll(), MadBase::LoadTrack(), NtpSRTrack::momentum, MadBase::ntpTrack, NtpSRMomentum::qp, and NtpSRMomentum::range.

Referenced by MadPIDAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

01368                                                                                                                                                      {
01369 
01370   // using the new version of energy corrections developed for Cedar/Daikon
01371 
01372   const float mumass=0.10566;
01373 
01374   if(LoadTrack(itrk)){
01375 
01376     Float_t mr=ntpTrack->momentum.range;
01377 
01378     if(mr>0) mr=FullyCorrectMomentumFromRange(mr,cx,reltype,corrver);
01379 
01380 
01381     float mc=0.0; // signed momentum from curvature
01382     if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
01383 
01384     mc=FullyCorrectSignedMomentumFromCurvature(mc,cx,reltype,corrver);
01385 
01386    
01387     if(opt==0){  //return the most appropriate measure of momentum
01388       if(IsFidAll(itrk) || ntpTrack->momentum.qp==0) {
01389         return sqrt(mr*mr+ mumass*mumass);
01390       }
01391       else {
01392         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01393         else return sqrt(mc*mc+ mumass*mumass);
01394       }
01395     }
01396     else if(opt==1) {
01397       //return curvature measurement
01398       // in R1.9 the tracker will apparently return (q/p)=0.0
01399       // maybe it's when a track looks perfectly rigid?
01400       // if so, we have to do something
01401       // I don't want to use the range, that could be very wrong
01402       // so, we'll return an obviously ridiculous value of 10 TeV
01403       if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
01404       else return sqrt(mc*mc+ mumass*mumass);
01405     }
01406     else if(opt==2) //return range measurement
01407       return sqrt(mr*mr+ mumass*mumass);
01408     else return 0;
01409   }
01410   return 0.;
01411 }

Float_t MadQuantities::RecoMuQP ( Int_t  itrk = 0  ) 

Definition at line 1434 of file MadQuantities.cxx.

References MadBase::LoadTrack(), NtpSRTrack::momentum, MadBase::ntpTrack, and NtpSRMomentum::qp.

Referenced by GetChargeSign(), and MakeValidationFile().

01434                                          {
01435   if(LoadTrack(itrk)) return ntpTrack->momentum.qp;
01436   return 0.;
01437 }

Float_t MadQuantities::RecoQENuEnergy ( Int_t  itrk = 0  ) 

Definition at line 1650 of file MadQuantities.cxx.

References GetChargeSign(), MadBase::LoadTrack(), RecoMuDCosNeu(), and RecoMuEnergy().

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::DataHist(), MakeValidationFile(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyCreatePANData(), MadCBSQEAnalysis::RecoAnalysisEnergy(), RecoQEQ2(), and MadCBSQEAnalysis::TestQEDiscrim().

01650                                                {
01651   
01652   if(!LoadTrack(itrk)) return 0.;
01653   Float_t nucleonMass = 0.93956563; //mass of neutron by default
01654   if(GetChargeSign(itrk)==1) nucleonMass = 0.93827231; //proton mass for nubar
01655   Float_t muonEnergy = RecoMuEnergy(0,itrk);
01656   Float_t muonMass = 0.10555;
01657   Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
01658   Float_t costhbl = RecoMuDCosNeu(itrk);
01659 
01660   Float_t Eneu = (nucleonMass*muonEnergy - muonMass*muonMass/2.)
01661     /(nucleonMass - muonEnergy + muonMomentum*costhbl);  
01662 
01663   return Eneu;
01664 }

Float_t MadQuantities::RecoQEQ2 ( Int_t  itrk = 0  ) 

Definition at line 1668 of file MadQuantities.cxx.

References MadBase::LoadTrack(), RecoMuDCosNeu(), RecoMuEnergy(), and RecoQENuEnergy().

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::DataHist(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), and MadEdAnalysis::MyCreatePANData().

01668                                          {
01669   
01670   if(!LoadTrack(itrk)) return 0.;
01671   Float_t Eneu = RecoQENuEnergy(itrk);
01672   Float_t muonEnergy = RecoMuEnergy(0,itrk);
01673   Float_t muonMass = 0.10555;
01674   Float_t muonMomentum = sqrt(muonEnergy*muonEnergy - muonMass*muonMass);
01675   Float_t costhbl = RecoMuDCosNeu(itrk);
01676   return -2.*Eneu*(muonEnergy-muonMomentum*costhbl)+muonMass*muonMass;
01677 
01678 }

Float_t MadQuantities::RecoShwEnergy ( Int_t  ishw = 0,
Int_t  opt = -1,
Int_t  det = 1,
bool  isdata = true,
int  mode = 1 
)

Definition at line 1528 of file MadQuantities.cxx.

References EnergyCorrections::CorrectShowerEnergyFar(), EnergyCorrections::CorrectShowerEnergyNear(), NtpSRShowerPulseHeight::EMgev, NtpSRStripPulseHeight::gev, CandShowerHandle::kCC, CandShowerHandle::kNC, CandShowerHandle::kWtCC, CandShowerHandle::kWtNC, NtpSRShowerPulseHeight::linCCgev, NtpSRShowerPulseHeight::linNCgev, MadBase::LoadShower(), MadBase::ntpShower, NtpSRShower::ph, NtpSRShower::shwph, NtpSRShowerPulseHeight::wtCCgev, and NtpSRShowerPulseHeight::wtNCgev.

Referenced by MadAnalysis::CreateANtpPAN(), MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::DataHist(), MadCBSQEAnalysis::MakeQEFile(), MakeValidationFile(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyCreatePANData(), PassQuasiCuts(), MadTestAnalysis::RecoAnalysisEnergy(), MadCBSQEAnalysis::RecoAnalysisEnergy(), MadPIDAnalysis::RecoAnalysisEnergy(), MadDpAnalysis::RecoAnalysisEnergy(), MadUserAnalysis::RecoAnalysisEnergy(), and ShowerValidation().

01528                                                                                            {
01529   // mode is used to deterimine if Niki or Andy's tuning should be used
01530   // they are more-or-less equivalent but we have choosen Niki's
01531   // mode=1 is set in the function declaration, so this choice will
01532   // be automatic for most people.
01533 
01534   Float_t shower_ene=0; 
01535 
01536   if(LoadShower(ishw)) { 
01537   if(det==1){
01538     if(opt==-1) shower_ene = ntpShower->ph.gev;
01539     if(opt==0)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,mode,isdata);
01540     if(opt==1)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,mode,isdata);
01541     if(opt==2)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,mode,isdata);
01542     if(opt==3)  shower_ene = CorrectShowerEnergyNear(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,mode,isdata);
01543     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01544     return shower_ene;
01545   }
01546    if(det==2){
01547     if(opt==-1) shower_ene = ntpShower->ph.gev;
01548     if(opt==0)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,mode,isdata);
01549     if(opt==1)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,mode,isdata);
01550     if(opt==2)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,mode,isdata);
01551     if(opt==3)  shower_ene = CorrectShowerEnergyFar(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,mode,isdata);
01552     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01553     return shower_ene;
01554 
01555   }
01556 } 
01557   return 0;
01558 }

Float_t MadQuantities::RecoShwEnergyNew ( Int_t  ishw,
Int_t  opt,
VldContext  cx,
ReleaseType::Release_t  reltype,
EnergyCorrections::WhichCorrection_t  corrver 
)

Definition at line 1563 of file MadQuantities.cxx.

References NtpSRShowerPulseHeight::EMgev, EnergyCorrections::FullyCorrectShowerEnergy(), NtpSRStripPulseHeight::gev, CandShowerHandle::kCC, CandShowerHandle::kNC, CandShowerHandle::kWtCC, CandShowerHandle::kWtNC, NtpSRShowerPulseHeight::linCCgev, NtpSRShowerPulseHeight::linNCgev, MadBase::LoadShower(), MadBase::ntpShower, NtpSRShower::ph, NtpSRShower::shwph, NtpSRShowerPulseHeight::wtCCgev, and NtpSRShowerPulseHeight::wtNCgev.

Referenced by MadPIDAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

01563                                                                                                                                                       {
01564 
01565 
01566   Float_t shower_ene=0; 
01567 
01568   if(LoadShower(ishw)) { 
01569 
01570     if(opt==-1) shower_ene = ntpShower->ph.gev;
01571     if(opt==0)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.linCCgev,CandShowerHandle::kCC,cx,reltype,corrver);
01572     if(opt==1)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.wtCCgev,CandShowerHandle::kWtCC,cx,reltype,corrver);
01573     if(opt==2)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.linNCgev,CandShowerHandle::kNC,cx,reltype,corrver);
01574     if(opt==3)  shower_ene = FullyCorrectShowerEnergy(ntpShower->shwph.wtNCgev,CandShowerHandle::kWtNC,cx,reltype,corrver);
01575     if(opt==4)  shower_ene = ntpShower->shwph.EMgev;
01576     return shower_ene;
01577   }
01578  
01579   return 0;
01580 }

Float_t MadQuantities::RecoShwEnergySqrt ( Int_t  ievt = 0  ) 

Definition at line 1584 of file MadQuantities.cxx.

References MuELoss::e, GetSqrtShwSC(), and GetSqrtTrkSkimSC().

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), and MadEdAnalysis::MyCreatePAN().

01584                                                   {
01585   Float_t sqrtSC = this->GetSqrtShwSC(ievt)+this->GetSqrtTrkSkimSC(ievt);
01586   if(sqrtSC==0) return 0;
01587   Float_t shwEnergy = (sqrtSC<3400) ? 0.108 + sqrtSC*0.0017 + 
01588     sqrtSC*sqrtSC*2.07e-7 : -3.76 + sqrtSC*0.00354;
01589   return shwEnergy;
01590 }

void MadQuantities::SetStdHepVars ( Int_t  itr  )  [protected]

Definition at line 542 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fGeantinoEnergy, fMaxFSKaonEnergy, fMaxFSNeutronEnergy, fMaxFSPiMinusEnergy, fMaxFSPiPlusEnergy, fMaxFSPiZeroEnergy, fMaxFSProtonEnergy, fNumFSGeantino, fNumFSKaon, fNumFSNeutron, fNumFSPiMinus, fNumFSPiPlus, fNumFSPiZero, fNumFSProton, fTotFSKaonEnergy, fTotFSNeutronEnergy, fTotFSPiMinusEnergy, fTotFSPiPlusEnergy, fTotFSPiZeroEnergy, fTotFSProtonEnergy, NtpMCStdHep::IdHEP, MadBase::isST, NtpMCStdHep::IstHEP, MadBase::lastEntry, MadBase::LoadStdHep(), MadBase::LoadTruth(), NtpMCStdHep::mass, NtpMCStdHep::mc, MadBase::mcrecord, MadBase::ntpStdHep, NtpMCStdHep::p4, NtpStRecord::stdhep, NtpMCRecord::stdhep, and MadBase::strecord.

Referenced by MadAnalysis::CreatePAN(), GeantinoEnergy(), MaxFSKaonEnergy(), MaxFSNeutronEnergy(), MaxFSPionEnergy(), MaxFSPiZeroEnergy(), MaxFSProtonEnergy(), NumFSGeantino(), NumFSKaon(), NumFSNeutron(), NumFSPion(), NumFSPiZero(), NumFSProton(), TotFSKaonEnergy(), TotFSNeutronEnergy(), TotFSPionEnergy(), TotFSPiZeroEnergy(), and TotFSProtonEnergy().

00542                                           {
00543 
00544   fNumFSGeantino = 0;
00545   fGeantinoEnergy = 0;  
00546   fNumFSProton = 0;
00547   fTotFSProtonEnergy = 0;
00548   fMaxFSProtonEnergy = 0;
00549   fNumFSNeutron = 0;
00550   fTotFSNeutronEnergy = 0;
00551   fMaxFSNeutronEnergy = 0;
00552   fNumFSPiPlus = 0;
00553   fTotFSPiPlusEnergy = 0;
00554   fMaxFSPiPlusEnergy = 0;
00555   fNumFSPiMinus = 0;
00556   fTotFSPiMinusEnergy = 0;
00557   fMaxFSPiMinusEnergy = 0;
00558   fNumFSPiZero = 0;
00559   fTotFSPiZeroEnergy = 0;
00560   fMaxFSPiZeroEnergy = 0;
00561   fNumFSKaon = 0;
00562   fTotFSKaonEnergy = 0;
00563   fMaxFSKaonEnergy = 0;
00564 
00565   fCurStdHepEntry = lastEntry;
00566   fCurTruthIndex = itr;
00567   
00568   if(!LoadTruth(itr)) return;
00569 
00570   //otherwise start calculating things
00571 
00572   //First get indices to use:
00573   TClonesArray* pointStdhepArray = NULL;
00574   if(isST) pointStdhepArray = (strecord->stdhep);
00575   else pointStdhepArray = (mcrecord->stdhep);
00576   TClonesArray& stdhepArray = *pointStdhepArray;
00577   Int_t nStdHep = stdhepArray.GetEntries();
00578   Int_t *indicesToUse = new Int_t[nStdHep];
00579   Int_t cnt = 0;
00580   for(int i=0;i<nStdHep;i++){
00581     LoadStdHep(i);
00582     if(ntpStdHep->mc==itr){
00583       indicesToUse[cnt] = i; 
00584       cnt++;
00585     }
00586   }
00587 
00588   //now loop through these and fill variables
00589   for(int i=0;i<cnt;i++){
00590     LoadStdHep(indicesToUse[i]);
00591     if(ntpStdHep->IstHEP==1){ //final state particles:
00592       Float_t mom = TMath::Sqrt(TMath::Abs(ntpStdHep->p4[3]*ntpStdHep->p4[3] - 
00593                                            ntpStdHep->mass*ntpStdHep->mass));
00594       if(ntpStdHep->IdHEP==28) {
00595         fNumFSGeantino+=1;
00596         fGeantinoEnergy+=mom;
00597       }
00598       else if(ntpStdHep->IdHEP==2212) {
00599         fNumFSProton+=1;
00600         fTotFSProtonEnergy+=mom;
00601         if(mom>fMaxFSProtonEnergy) fMaxFSProtonEnergy = mom;
00602       }
00603       else if(ntpStdHep->IdHEP==2112) {
00604         fNumFSNeutron+=1;
00605         fTotFSNeutronEnergy+=mom;
00606         if(mom>fMaxFSNeutronEnergy) fMaxFSNeutronEnergy = mom;
00607       }
00608       else if(ntpStdHep->IdHEP==211) {
00609         fNumFSPiPlus+=1;
00610         fTotFSPiPlusEnergy+=mom;
00611         if(mom>fMaxFSPiPlusEnergy) fMaxFSPiPlusEnergy = mom;
00612       }
00613       else if(ntpStdHep->IdHEP==-211) {
00614         fNumFSPiMinus+=1;
00615         fTotFSPiMinusEnergy+=mom;
00616         if(mom>fMaxFSPiMinusEnergy) fMaxFSPiMinusEnergy = mom;
00617       }
00618       else if(ntpStdHep->IdHEP==111) {
00619         fNumFSPiZero+=1;
00620         fTotFSPiZeroEnergy+=mom;
00621         if(mom>fMaxFSPiZeroEnergy) fMaxFSPiZeroEnergy = mom;
00622       }
00623       else if(abs(ntpStdHep->IdHEP)==321) {
00624         fNumFSKaon+=1;
00625         fTotFSKaonEnergy+=mom;
00626         if(mom>fMaxFSKaonEnergy) fMaxFSPiZeroEnergy = mom;
00627       }
00628     }
00629   }
00630 
00631   delete[] indicesToUse;
00632 }

void MadQuantities::ShowerValidation ( char *  fileName,
Int_t  startEnt = 0 
)

Definition at line 2264 of file MadQuantities.cxx.

References ClosestDistanceToEvent(), NtpSRShower::clu, NtpTHTrack::completeall, NtpTHEvent::completeall, NtpTHShower::completeall, NtpTHShower::completeslc, NtpTHTrack::completeslc, NtpTHEvent::completeslc, MadBase::eventSummary, BeamMonSpill::fBpmInt, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTortgt, BDSpillAccessor::Get(), MadBase::GetEntry(), GetNShowerStrip(), VldContext::GetSimFlag(), BeamMonSpill::GetStatusBits(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpSRCluster::id, NtpSREvent::index, SpillTimeFinder::Instance(), NtpMCTruth::iresonance, IsFidVtxEvt(), IsFidVtxTrk(), Dbi::kDefaultTask, Dbi::kDisabled, SimFlag::kMC, NtpSRShowerPulseHeight::linCCgev, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShowerAtTrackVertex(), BDSpillAccessor::LoadSpill(), MadBase::LoadStrip(), MadBase::LoadTHEvent(), MadBase::LoadTHShower(), MadBase::LoadTHTrack(), MadBase::LoadTruthForRecoTH(), NtpSRPlane::n, NtpSRShower::ncluster, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSREventSummary::nstrip, NtpSRCluster::nstrip, NtpSRShower::nstrip, MadBase::ntpCluster, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTHEvent, MadBase::ntpTHShower, MadBase::ntpTHTrack, MadBase::ntpTrack, MadBase::ntpTruth, NtpSREvent::ntrack, NtpSRShower::nUcluster, NtpSRShower::nVcluster, NtpMCTruth::p4mu1, NtpMCTruth::p4neu, NtpMCTruth::p4shw, NtpSRPulseHeight::pe, NtpSRCluster::ph, NtpSREvent::ph, NtpSRStrip::ph1, NtpSRShower::plane, NtpSRStrip::plane, NtpSRCluster::planeview, NtpTHEvent::purity, NtpTHTrack::purity, NtpTHShower::purity, RecoMuEnergy(), RecoShwEnergy(), NtpSRShower::shwph, NtpSRShower::stp, NtpSRCluster::stp, NtpSRStrip::strip, NtpSRVertex::t, NtpSRStrip::time0, NtpSRStrip::time1, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, NtpMCTruth::y, and NtpSRVertex::z.

02265 {
02266 
02267   this->GetEntry(0);
02268   const bool file_is_mc
02269     =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
02270   Double_t pot = 0;
02271 
02272   char savename[256];
02273   sprintf(savename,"ShowerValidation_%s.root",fileName);
02274   TFile *save = new TFile(savename,"RECREATE");
02275   save->cd();
02276 
02277   //strip level:
02278   TH1F *stripPH = new TH1F("stripPH","stripPH",1000,0,100);
02279   TH1F *stripTime = new TH1F("stripTime","stripTime",1000,-1000,1000);
02280   TH1F *StripPHStripCut = new TH1F("StripPHStripCut","StripPHStripCut",1000,0,100);
02281   TH1F *StripTimeStripCut = new TH1F("StripTimeStripCut","StripTimeStripCut",1000,-1000,1000);
02282   TH1F *StripPHSSCut = new TH1F("StripPHSSCut","StripPHSSCut",1000,0,100);
02283   TH1F *StripTimeSSCut = new TH1F("StripTimeSSCut","StripTimeSSCut",1000,-1000,1000);
02284   TH1F *StripPHBothCut = new TH1F("StripPHBothCut","StripPHBothCut",1000,0,100);
02285   TH1F *StripTimeBothCut = new TH1F("StripTimeBothCut","StripTimeBothCut",1000,-1000,1000);
02286 
02287   TH1F *stripPH_Phys = new TH1F("stripPH_Phys","stripPH_Phys",1000,0,100);
02288   TH1F *stripTime_Phys = new TH1F("stripTime_Phys","stripTime_Phys",1000,-1000,1000);
02289   TH1F *StripPHStripCut_Phys = new TH1F("StripPHStripCut_Phys","StripPHStripCut_Phys",1000,0,100);
02290   TH1F *StripTimeStripCut_Phys = new TH1F("StripTimeStripCut_Phys","StripTimeStripCut_Phys",1000,-1000,1000);
02291   TH1F *StripPHSSCut_Phys = new TH1F("StripPHSSCut_Phys","StripPHSSCut_Phys",1000,0,100);
02292   TH1F *StripTimeSSCut_Phys = new TH1F("StripTimeSSCut_Phys","StripTimeSSCut_Phys",1000,-1000,1000);
02293   TH1F *StripPHBothCut_Phys = new TH1F("StripPHBothCut_Phys","StripPHBothCut_Phys",1000,0,100);
02294   TH1F *StripTimeBothCut_Phys = new TH1F("StripTimeBothCut_Phys","StripTimeBothCut_Phys",1000,-1000,1000);
02295 
02296   TH2F *stripTimeDiff = new TH2F("stripTimeDiff","stripTimeDiff",1000,0,1000,100,0,10);
02297 
02298   TH1F *multiplyHitStripFraction = new TH1F("multiplyHitStripFraction",
02299                                             "multiplyHitStripFraction",
02300                                             100,0,1.01);
02301   TH1F *multiplyHitStripFractionPhys = new TH1F("multiplyHitStripFractionPhys",
02302                                             "multiplyHitStripFractionPhys",
02303                                             100,0,1.01);
02304 
02305   //shower level:  
02306   TH1F *recoEshw = new TH1F("recoEshw","recoEshw",100,0,50);
02307   TH1F *recoEshwStripCut = new TH1F("recoEshwStripCut","recoEshwStripCut",100,0,50);
02308   TH1F *recoEshwSSCut = new TH1F("recoEshwSSCut","recoEshwSSCut",100,0,50);
02309   TH1F *recoEshwBothCut = new TH1F("recoEshwBothCut","recoEshwBothCut",100,0,50);
02310 
02311   TH1F *shwLengthTime = new TH1F("shwLengthTime","shwLengthTime",1000,0,1000);
02312   TH1F *shwLengthPlane = new TH1F("shwLengthPlane","shwLengthPlane",60,0,60);
02313   TH1F *shwStripN = new TH1F("shwStripN","shwStripN",200,0,200);
02314   TH1F *shwPhysStpN = new TH1F("shwPhysStpN","shwPhysStpN",200,0,200);
02315   TH1F *shwClusterN = new TH1F("shwClusterN","shwClusterN",60,0,60);
02316   TH1F *shwPhysCluN = new TH1F("shwPhysCluN","shwPhysCluN",60,0,60);
02317 
02318   //Truth Helper histos:
02319   TH1F *eventCompleteAll = new TH1F("eventCompleteAll",
02320                                     "Event Completeness",100,0,1);
02321   TH1F *eventCompleteSlc = new TH1F("eventCompleteSlc",
02322                                     "Event Completeness in Slice",100,0,1);
02323   TH1F *eventPurity = new TH1F("eventPurity","eventPurity",100,0,1);
02324   TH1F *eventEnergyRes = new TH1F("eventEnergyRes","eventEnergyRes",100,-1,9);
02325   TH1F *eventClosestDistance = new TH1F("eventClosestDistance","eventClosestDistance",
02326                                         1000,0,200);
02327   TH1F *eventLowResClosestDistance = new TH1F("eventLowResClosestDistance",
02328                                               "eventLowResClosestDistance",
02329                                               1000,0,200);
02330   TH2F *eventClosestDistTime = new TH2F("eventClosestDistTime","eventClosestDistTime",
02331                                         200,0,40,200,0,40);
02332   TH2F *eventLRClosestDistTime = new TH2F("eventLRClosestDistTime","eventLRClosestDistTime",
02333                                           200,0,40,200,0,40);
02334   TH2F *eventLRTrkShw = new TH2F("eventLRTrkShw","eventLRTrkShw",4,-0.5,3.5,4,-0.5,3.5);
02335   
02336 
02337   TH1F *trackCompleteAll = new TH1F("trackCompleteAll",
02338                                     "Track Completeness",100,0,1);
02339   TH1F *trackCompleteSlc = new TH1F("trackCompleteSlc",
02340                                     "Track Completeness in Slice",100,0,1);
02341   TH1F *trackPurity = new TH1F("trackPurity","trackPurity",100,0,1);
02342   TH1F *trackEnergyRes = new TH1F("trackEnergyRes","trackEnergyRes",100,-1,9);
02343 
02344   TH1F *showerCompleteAll = new TH1F("showerCompleteAll",
02345                                      "Shower Completeness",100,0,1);
02346   TH1F *showerCompleteSlc = new TH1F("showerCompleteSlc",
02347                                      "Shower Completeness in Slice",100,0,1);
02348   TH1F *showerPurity = new TH1F("showerPurity","showerPurity",100,0,1);
02349   TH1F *showerEnergyRes = new TH1F("showerEnergyRes","showerEnergyRes",100,-1,9);
02350   TH1F *showerVtxDist = new TH1F("showerVtxDist","showerVtxDist",100,0,6);
02351 
02352   //investigate low completeness showers:
02353   TH1F *showerLowCompleteProcess = new TH1F("showerLowCompleteProcess",
02354                                             "showerLowCompleteProcess",4,-0.5,3.5);
02355   TH1F *showerLowCompleteEnergy = new TH1F("showerLowCompleteEnergy",
02356                                            "showerLowCompleteEnergy",100,0,10);
02357   TH1F *showerLowCompleteCluID = new TH1F("showerLowCompleteCluID",
02358                                           "showerLowCompleteCluID",14,-0.5,13.5);
02359   TH1F *showerCompleteSlcSSOnly = new TH1F("showerCompleteSlcSSOnly",
02360                                           "showerCompleteSlcSSOnly",100,0,1);
02361   TH1F *showerLowCompleteVtxDist = new TH1F("showerLowCompleteVtxDist",
02362                                             "showerLowCompleteVtxDist",
02363                                             100,0,6);
02364   TH1F *showerLowCompleteProcessSSOnly = new TH1F("showerLowCompleteProcessSSOnly",
02365                                                   "showerLowCompleteProcessSSOnly",
02366                                                   4,-0.5,3.5);
02367 
02368   bool beamdb=false;
02369   bool checkeddb=false;
02370   for(int i=startEnt;i<Nentries;i++){
02371     if(i%1000==0) std::cout << float(i)*100./float(Nentries) 
02372                             << "% done" << std::endl;
02373     if(!this->GetEntry(i)) continue;
02374     if(!checkeddb){
02375       DbiResultPtr<BeamMonSpill> testbs(ntpHeader->GetVldContext(), 
02376                                         Dbi::kDefaultTask, 
02377                                         Dbi::kDisabled);
02378       checkeddb=true;
02379       if(testbs.GetNumRows()>0){
02380         beamdb=true;
02381         cout<<"Found the Beam Monitoring Database Tables"<<endl;
02382       }
02383       else{     
02384         cout<<"Didn't find the Beam Monitoring Database Tables"<<endl;
02385         cout<<"Resorting to old beamsummary ntuples"<<endl;
02386       }
02387     }
02388 
02389     if(!file_is_mc){ //file is mc
02390       Bool_t goodbeam = false;
02391       VldContext vc = ntpHeader->GetVldContext();
02392       VldTimeStamp vts = vc.GetTimeStamp();
02393       if(beamdb){ //beam db
02394         BDSpillAccessor &sa = BDSpillAccessor::Get();
02395         const BeamMonSpill *bs = sa.LoadSpill(vts);
02396         SpillTimeFinder &sfind= SpillTimeFinder::Instance();
02397         Double_t hpos2=0.;
02398         Double_t vpos2=0.;
02399         Double_t btint=0;
02400         for(int bt=0;bt<6;bt++){
02401           hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
02402           vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);
02403           btint+=bs->fBpmInt[bt];
02404         }
02405         if(btint!=0){
02406           hpos2=hpos2*1000./btint;//make it in mm
02407           vpos2=vpos2*1000./btint;//make it in mm
02408         }
02409         else{
02410           hpos2=-999.;
02411           vpos2=-999.;
02412         }
02413         /*
02414           cout << bs->GetStatusBits().horn_on << " "
02415           << bs->GetStatusBits().target_in << " "
02416           << bs->fTortgt << " "
02417           << bs->fProfWidX*1000. << " "
02418           << bs->fProfWidY*1000. << " "
02419           << hpos2 << " " << vpos2 << " "
02420           << fabs(sfind.GetTimeToNearestSpill(vc)) << endl;
02421         */
02422         if(bs->GetStatusBits().horn_on && 
02423            bs->GetStatusBits().target_in &&
02424            bs->fTortgt > 0.1 &&
02425            bs->fProfWidX*1000.<1.5 &&
02426            bs->fProfWidY*1000.<2.0 &&
02427            hpos2>-2 && hpos2<0 &&    
02428            vpos2>0 && vpos2<2 &&
02429            fabs(sfind.GetTimeToNearestSpill(vc))<2 ) goodbeam=true;
02430         if(goodbeam) pot += bs->fTortgt;
02431       }
02432       if(!goodbeam) continue;
02433     }
02434 
02435     //Still have all good snarls here    
02436 
02437     if(eventSummary->nevent==0) continue;
02438 
02439     //find and store earliest time of hit strip
02440     Double_t stripArrayTime[500][192][2];
02441     for(int ii=0;ii<500;ii++) { 
02442       for(int jj=0;jj<192;jj++) {
02443         stripArrayTime[ii][jj][0] = stripArrayTime[ii][jj][1] = -1.0;
02444       }
02445     }
02446     for(int j=0;j<int(eventSummary->nstrip);j++){
02447       if(!LoadStrip(j)) continue;
02448       Int_t pl = ntpStrip->plane;
02449       Int_t st = ntpStrip->strip;
02450       if(ntpStrip->time1<stripArrayTime[pl][st][1] || 
02451          stripArrayTime[pl][st][1]<=0) {
02452         stripArrayTime[pl][st][1] = ntpStrip->time1;
02453       }
02454     }
02455 
02456     Int_t is_fid = 0; 
02457     for(int j=0;j<eventSummary->nevent;j++){ 
02458       if(!LoadEvent(j)) continue; //no event found
02459 
02460       int track_index   = -1;
02461       int shower_index  = -1;
02462       Bool_t vertex_shower = true;
02463       if(LoadLargestTrackFromEvent(j,track_index)) {
02464         if(!LoadShowerAtTrackVertex(j,track_index,shower_index)) {        
02465           LoadLargestShowerFromEvent(j,shower_index);
02466           vertex_shower = false;
02467         }
02468       }
02469       else LoadLargestShowerFromEvent(j,shower_index);
02470 
02471       is_fid = 1;
02472       if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
02473       if(ntpEvent->ntrack>0&&!IsFidVtxTrk(track_index)) is_fid = 0;
02474       if(is_fid==0) continue;
02475 
02476       //Still have all fiducial events from good snarls here
02477       Float_t phwCluID = 0;
02478       Float_t phClu = 0;
02479       if(ntpEvent->nshower>0) { //is shower?
02480 
02481         recoEshw->Fill(ntpShower->shwph.linCCgev);
02482         shwLengthPlane->Fill(ntpShower->plane.n);
02483         Double_t shwstptimemax = 0;
02484         Double_t shwstptimemin = 0;
02485         for(int k=0;k<ntpShower->nstrip;k++){
02486           if(!LoadStrip(ntpShower->stp[k])) continue;
02487           Double_t t0 = ntpStrip->time0;
02488           Double_t t1 = ntpStrip->time1;
02489           Double_t avet = 0;
02490           if(t0>0&&t1>0) avet = (t0+t1)/2.;
02491           else if(t0>0) avet = t0;
02492           else if(t1>0) avet = t1;
02493           if(k==0) shwstptimemin = shwstptimemax = avet;
02494           else {
02495             if(avet>shwstptimemax) shwstptimemax = avet;
02496             if(avet<shwstptimemin) shwstptimemin = avet;
02497           }
02498         }
02499         shwLengthTime->Fill(1e9*(shwstptimemax-shwstptimemin));
02500 
02501         Int_t nMultiplyHitStrips = 0;
02502         Int_t nPhysClusterU = 0;
02503         Int_t nPhysClusterV = 0;        
02504         Int_t *clusters = ntpShower->clu;
02505         for(int k=0; k<ntpShower->ncluster; k++){
02506           Int_t index = clusters[k];
02507           if(!LoadCluster(index)) continue;
02508           phwCluID += ntpCluster->id*ntpCluster->ph.gev;
02509           phClu += ntpCluster->ph.gev;
02510           for(int l=0;l<ntpCluster->nstrip;l++){
02511             if(!LoadStrip(ntpCluster->stp[l])) continue;
02512             Int_t pl = ntpStrip->plane;
02513             Int_t st = ntpStrip->strip;
02514             if(ntpStrip->time1>stripArrayTime[pl][st][1]) {       
02515               stripTimeDiff->Fill((ntpStrip->time1 - 
02516                                    stripArrayTime[pl][st][1])*1e9,
02517                                   ntpStrip->ph1.pe);
02518               nMultiplyHitStrips+=1;
02519             }
02520           }
02521           if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02522             if(ntpCluster->planeview==2) nPhysClusterU+=1;
02523             if(ntpCluster->planeview==3) nPhysClusterV+=1;
02524           }
02525         }
02526 
02527         if(phClu>0) phwCluID/=phClu;
02528 
02529         shwClusterN->Fill(ntpShower->nUcluster+ntpShower->nVcluster);
02530         shwPhysCluN->Fill(nPhysClusterU+nPhysClusterV);
02531         
02532         Int_t *nss = GetNShowerStrip(shower_index);
02533         Int_t *npss = GetNShowerStrip(shower_index,1);
02534     
02535         shwStripN->Fill(nss[0]+nss[1]);
02536         shwPhysStpN->Fill(npss[0]+npss[1]);
02537 
02538         multiplyHitStripFraction->Fill(Float_t(nMultiplyHitStrips)/
02539                                        Float_t(ntpShower->nstrip));
02540         
02541         if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02542           recoEshwStripCut->Fill(ntpShower->shwph.linCCgev);
02543           if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02544             recoEshwBothCut->Fill(ntpShower->shwph.linCCgev);
02545           }
02546         }
02547         if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02548           recoEshwSSCut->Fill(ntpShower->shwph.linCCgev);
02549         }
02550         
02551         multiplyHitStripFractionPhys->Fill(Float_t(nMultiplyHitStrips)/
02552                                            Float_t(ntpShower->nstrip));
02553         
02554         
02555         for(int k=0; k<ntpShower->ncluster; k++){
02556           Int_t index = clusters[k];
02557           if(!LoadCluster(index)) continue;
02558           for(int l=0;l<ntpCluster->nstrip;l++){
02559             if(!LoadStrip(ntpCluster->stp[l])) continue;
02560             stripPH->Fill(ntpStrip->ph1.pe);
02561             stripTime->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));       
02562             if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02563               stripPH_Phys->Fill(ntpStrip->ph1.pe);
02564               stripTime_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02565             }
02566             
02567             if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02568               StripPHStripCut->Fill(ntpStrip->ph1.pe);
02569               StripTimeStripCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));               
02570               if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02571                 StripPHStripCut_Phys->Fill(ntpStrip->ph1.pe);
02572                 StripTimeStripCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02573               }
02574             }
02575             
02576             if(1.-float(npss[0]+npss[1])/float(nss[0]+nss[1]) < 0.85) {
02577               StripPHSSCut->Fill(ntpStrip->ph1.pe);
02578               StripTimeSSCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));                  
02579               if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02580                 StripPHSSCut_Phys->Fill(ntpStrip->ph1.pe);
02581                 StripTimeSSCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02582               }
02583               
02584               if(Float_t(nMultiplyHitStrips)/Float_t(ntpShower->nstrip)<0.3) {
02585                 StripPHBothCut->Fill(ntpStrip->ph1.pe);
02586                 StripTimeBothCut->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));              
02587                 if(ntpCluster->id==0||ntpCluster->id==1||ntpCluster->id==3){
02588                   StripPHBothCut_Phys->Fill(ntpStrip->ph1.pe);
02589                   StripTimeBothCut_Phys->Fill(1e9*(ntpStrip->time1-ntpEvent->vtx.t));
02590                 }
02591               }
02592             }
02593           }
02594         }
02595         delete [] nss;
02596         delete [] npss;
02597       }
02598       //if this is MC:
02599       Int_t mcevent = -1;
02600       if(file_is_mc&&LoadTruthForRecoTH(j,mcevent)) {
02601         if(LoadTHEvent(j)) {
02602           eventCompleteAll->Fill(ntpTHEvent->completeall);
02603           eventCompleteSlc->Fill(ntpTHEvent->completeslc);
02604           eventPurity->Fill(ntpTHEvent->purity);
02605           Float_t true_en_to_use = ntpTruth->p4neu[3];
02606           if(ntpTruth->iaction==0) true_en_to_use = ntpTruth->p4shw[3];
02607           Float_t reco_en_to_use = RecoMuEnergy(track_index,0) + 
02608             RecoShwEnergy(shower_index,0);
02609           reco_en_to_use = ntpEvent->ph.gev;
02610           eventEnergyRes->Fill((reco_en_to_use - 
02611                                 true_en_to_use)/true_en_to_use);
02612           Float_t *closestDistanceToOtherEvent = 
02613             ClosestDistanceToEvent(ntpEvent->index,0.033); 
02614           eventClosestDistance->Fill(closestDistanceToOtherEvent[0]);
02615           eventClosestDistTime->Fill(closestDistanceToOtherEvent[2],
02616                                      closestDistanceToOtherEvent[1]);
02617           if((reco_en_to_use - true_en_to_use)/true_en_to_use<-0.9) {
02618             eventLowResClosestDistance->Fill(closestDistanceToOtherEvent[0]);
02619             eventLRClosestDistTime->Fill(closestDistanceToOtherEvent[2],
02620                                          closestDistanceToOtherEvent[1]);
02621             eventLRTrkShw->Fill(ntpEvent->ntrack,ntpEvent->nshower);
02622           }
02623           delete [] closestDistanceToOtherEvent;
02624         }
02625         if(LoadTHTrack(track_index)){
02626           trackCompleteAll->Fill(ntpTHTrack->completeall);
02627           trackCompleteSlc->Fill(ntpTHTrack->completeslc);
02628           trackPurity->Fill(ntpTHTrack->purity);
02629           Float_t true_en_to_use = ntpTruth->p4mu1[3];
02630           Float_t reco_en_to_use = RecoMuEnergy(track_index,0);
02631           if(true_en_to_use > 0)
02632             trackEnergyRes->Fill((reco_en_to_use - 
02633                                   true_en_to_use)/true_en_to_use);
02634           else trackEnergyRes->Fill(-1);
02635         }
02636         if(LoadTHShower(shower_index)&&vertex_shower){
02637           showerCompleteAll->Fill(ntpTHShower->completeall);
02638           showerCompleteSlc->Fill(ntpTHShower->completeslc);
02639           showerPurity->Fill(ntpTHShower->purity);
02640           Float_t true_en_to_use = ntpTruth->p4neu[3]*ntpTruth->y;
02641           if(true_en_to_use == 0) true_en_to_use = ntpTruth->p4shw[3];
02642           Float_t reco_en_to_use = RecoShwEnergy(shower_index,0);
02643           showerEnergyRes->Fill((reco_en_to_use - 
02644                                  true_en_to_use)/true_en_to_use);
02645           if(ntpTrack) showerVtxDist->Fill(TMath::Abs(ntpTrack->vtx.z - 
02646                                                       ntpShower->vtx.z));
02647           if(phwCluID>=5) {
02648             showerCompleteSlcSSOnly->Fill(ntpTHShower->completeslc);
02649             if(ntpTHShower->completeslc<0.1){
02650               showerLowCompleteProcessSSOnly->Fill((ntpTruth->iresonance-1000) * 
02651                                                    ntpTruth->iaction);
02652             }
02653           }
02654           if(ntpTHShower->completeslc<0.1){
02655             showerLowCompleteProcess->Fill((ntpTruth->iresonance-1000) * 
02656                                            ntpTruth->iaction);
02657             showerLowCompleteEnergy->Fill(ntpTruth->p4shw[3]);
02658             showerLowCompleteCluID->Fill(phwCluID);
02659             if(ntpTrack) showerLowCompleteVtxDist->Fill(TMath::Abs(ntpTrack->vtx.z - 
02660                                                                    ntpShower->vtx.z));
02661           }
02662         }
02663       }
02664     }
02665   }
02666 
02667   recoEshw->Write();
02668   recoEshwSSCut->Write();
02669   recoEshwStripCut->Write();
02670   recoEshwBothCut->Write();
02671 
02672   stripPH->Write();
02673   stripTime->Write();  
02674   stripPH_Phys->Write();
02675   stripTime_Phys->Write();  
02676 
02677   StripPHStripCut->Write();
02678   StripTimeStripCut->Write();
02679   StripPHStripCut_Phys->Write();
02680   StripTimeStripCut_Phys->Write();
02681 
02682   StripPHSSCut->Write();
02683   StripTimeSSCut->Write();
02684   StripPHSSCut_Phys->Write();
02685   StripTimeSSCut_Phys->Write();
02686 
02687   StripPHBothCut->Write();
02688   StripTimeBothCut->Write();
02689   StripPHBothCut_Phys->Write();
02690   StripTimeBothCut_Phys->Write();
02691 
02692   stripTimeDiff->Write();
02693   multiplyHitStripFraction->Write();
02694   multiplyHitStripFractionPhys->Write();
02695 
02696   shwLengthTime->Write();
02697   shwLengthPlane->Write();
02698   shwStripN->Write();
02699   shwPhysStpN->Write();
02700   shwClusterN->Write();
02701   shwPhysCluN->Write();
02702 
02703   eventCompleteAll->Write();
02704   eventCompleteSlc->Write();
02705   eventPurity->Write();
02706   eventEnergyRes->Write();
02707   eventClosestDistance->Write();
02708   eventLowResClosestDistance->Write();
02709   eventClosestDistTime->Write();
02710   eventLRClosestDistTime->Write();
02711   eventLRTrkShw->Write();
02712 
02713   trackCompleteAll->Write();
02714   trackCompleteSlc->Write();
02715   trackPurity->Write();
02716   trackEnergyRes->Write();
02717 
02718   showerCompleteAll->Write();
02719   showerCompleteSlc->Write();
02720   showerPurity->Write();
02721   showerEnergyRes->Write();
02722   showerVtxDist->Write();
02723 
02724   showerCompleteSlcSSOnly->Write();
02725   showerLowCompleteProcess->Write();
02726   showerLowCompleteEnergy->Write();
02727   showerLowCompleteCluID->Write();
02728   showerLowCompleteVtxDist->Write();
02729   showerLowCompleteProcessSSOnly->Write();
02730 
02731   save->Close();
02732   cout << "POT = " << pot << endl;
02733   return;
02734 }

Float_t * MadQuantities::Target4Vector ( Int_t  itr = 0  ) 

Definition at line 1209 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, and NtpMCTruth::p4tgt.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), and MadEdAnalysis::MyCreatePAN().

01209                                               {
01210   Float_t *vec = new Float_t[4];
01211   vec[0] = 0.; vec[1] = 0.; vec[2] = 0.; vec[3] = 0.;
01212   if(!LoadTruth(itr)) return vec;
01213   vec[0] = ntpTruth->p4tgt[0];
01214   vec[1] = ntpTruth->p4tgt[1];
01215   vec[2] = ntpTruth->p4tgt[2];
01216   vec[3] = ntpTruth->p4tgt[3];
01217   return vec;
01218 }

Float_t MadQuantities::TotFSKaonEnergy ( Int_t  itr = 0  ) 

Definition at line 713 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fTotFSKaonEnergy, MadBase::lastEntry, and SetStdHepVars().

00713                                                {
00714   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00715     SetStdHepVars(itr);
00716   return fTotFSKaonEnergy;
00717 }

Float_t MadQuantities::TotFSNeutronEnergy ( Int_t  itr = 0  ) 
Float_t MadQuantities::TotFSPionEnergy ( Int_t  itr = 0,
Int_t  pm = 0 
)

Definition at line 690 of file MadQuantities.cxx.

References fCurStdHepEntry, fCurTruthIndex, fTotFSPiMinusEnergy, fTotFSPiPlusEnergy, MadBase::lastEntry, and SetStdHepVars().

Referenced by MadMKAnalysis::CreatePAN(), and MadTVAnalysis::CreatePAN().

00690                                                         {
00691   if(fCurStdHepEntry!=lastEntry||fCurTruthIndex!=itr) 
00692     SetStdHepVars(itr);
00693   if(pm==-1) return fTotFSPiMinusEnergy;
00694   else if(pm==1) return fTotFSPiPlusEnergy;
00695   return fTotFSPiPlusEnergy+fTotFSPiMinusEnergy;
00696 }

Float_t MadQuantities::TotFSPiZeroEnergy ( Int_t  itr = 0  ) 
Float_t MadQuantities::TotFSProtonEnergy ( Int_t  itr = 0  ) 
Float_t MadQuantities::TrueLeptonEnergy ( Int_t  itr = 0  ) 

Definition at line 1263 of file MadQuantities.cxx.

References NtpMCTruth::iaction, NtpMCTruth::inu, MadBase::LoadTruth(), MadBase::ntpTruth, NtpMCTruth::p4el1, NtpMCTruth::p4mu1, and NtpMCTruth::p4tau.

Referenced by MadAnalysis::CreatePAN().

01263                                                 {
01264   if(!LoadTruth(itr)) return 0;
01265   if(ntpTruth->iaction==0) return 0;
01266   else if(ntpTruth->inu==12) return fabs(ntpTruth->p4el1[3]);
01267   else if(ntpTruth->inu==14) return fabs(ntpTruth->p4mu1[3]);
01268   else if(ntpTruth->inu==16) return fabs(ntpTruth->p4tau[3]);
01269   return 0;
01270 }

Float_t MadQuantities::TrueMuDCosNeu ( Int_t  itr = 0  ) 

Definition at line 1304 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, NtpMCTruth::p4mu1, and NtpMCTruth::p4neu.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::MCHist(), MadEdAnalysis::MyCreatePAN(), and MadEdAnalysis::MyMakeQEFile().

01304                                              { //ds_mu/ds_nu
01305   if(!LoadTruth(itr)) return 0.;
01306   TVector3 *nuvec = new TVector3(ntpTruth->p4neu[0],
01307                                  ntpTruth->p4neu[1],
01308                                  ntpTruth->p4neu[2]);
01309   TVector3 *muvec = new TVector3(ntpTruth->p4mu1[0],
01310                                  ntpTruth->p4mu1[1],
01311                                  ntpTruth->p4mu1[2]);  
01312   Float_t MuAng = nuvec->Angle(*muvec); //angle in rads
01313   delete nuvec;
01314   delete muvec;
01315   return TMath::Cos(MuAng);
01316 }

Float_t MadQuantities::TrueMuDCosZ ( Int_t  itr = 0  ) 

Definition at line 1293 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, and NtpMCTruth::p4mu1.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MadEdAnalysis::MyCreatePAN(), MadEdAnalysis::MyMakeQEFile(), and MadAnalysis::RecoMC().

01293                                            { //dz/ds
01294   if(!LoadTruth(itr)) return 0.;
01295   Float_t mom = sqrt((ntpTruth->p4mu1[0]*ntpTruth->p4mu1[0]) + 
01296                      (ntpTruth->p4mu1[1]*ntpTruth->p4mu1[1]) + 
01297                      (ntpTruth->p4mu1[2]*ntpTruth->p4mu1[2]));
01298   if(mom==0) return 0;
01299   return ntpTruth->p4mu1[2]/mom;
01300 }

Float_t MadQuantities::TrueMuEnergy ( Int_t  itr = 0  ) 
Float_t MadQuantities::TrueMuQP ( Int_t  itr = 0  ) 

Definition at line 1274 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, and NtpMCTruth::p4mu1.

Referenced by MakeValidationFile().

01274                                         {
01275   if(!LoadTruth(itr)) return 0;
01276   if(ntpTruth->p4mu1[3]!=0) {
01277     Float_t p = 1./sqrt(ntpTruth->p4mu1[3]*ntpTruth->p4mu1[3]-0.10555*0.10555);
01278     Float_t q = ntpTruth->p4mu1[3]/fabs(ntpTruth->p4mu1[3]);
01279     return q*p;
01280   }
01281   return 0;
01282 }

Float_t MadQuantities::TrueNuEnergy ( Int_t  itr = 0  ) 
Float_t * MadQuantities::TrueNuMom ( Int_t  itr = 0  ) 

Definition at line 1228 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, and NtpMCTruth::p4neu.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), and MadEdAnalysis::MyCreatePAN().

01228                                           {
01229   Float_t *mom = new Float_t[3];  
01230   mom[0] = 0.; 
01231   mom[1] = 0.; 
01232   mom[2] = 0.;
01233   if(!LoadTruth(itr)) return mom;
01234   mom[0] = ntpTruth->p4neu[0];
01235   mom[1] = ntpTruth->p4neu[1];  
01236   mom[2] = ntpTruth->p4neu[2];
01237   return mom;
01238 }

Float_t MadQuantities::TrueShwEnergy ( Int_t  itr = 0  ) 
Float_t * MadQuantities::TrueVtx ( Int_t  itr = 0  ) 

Definition at line 1242 of file MadQuantities.cxx.

References MadBase::LoadTruth(), MadBase::ntpTruth, NtpMCTruth::vtxx, NtpMCTruth::vtxy, and NtpMCTruth::vtxz.

Referenced by MadDpAnalysis::CreatePAN(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), MadAnalysis::CreatePAN(), MakeValidationFile(), MadEdAnalysis::MyCreatePAN(), and MadEdAnalysis::MyMakeQEFile().

01242                                         {
01243   Float_t *vtx = new Float_t[3];  
01244   vtx[0] = 0.; 
01245   vtx[1] = 0.; 
01246   vtx[2] = 0.;
01247   if(!LoadTruth(itr)) return vtx;
01248   vtx[0] = ntpTruth->vtxx;
01249   vtx[1] = ntpTruth->vtxy;  
01250   vtx[2] = ntpTruth->vtxz;
01251   return vtx;
01252 }

Float_t MadQuantities::W2 ( Int_t  itr = 0  ) 
Float_t MadQuantities::X ( Int_t  itr = 0  ) 
Float_t MadQuantities::Y ( Int_t  itr = 0  ) 

Member Data Documentation

Int_t MadQuantities::fCurStdHepEntry [protected]
Int_t MadQuantities::fCurTruthIndex [protected]
Float_t MadQuantities::fGeantinoEnergy [protected]
Float_t MadQuantities::fMaxFSKaonEnergy [protected]
Float_t MadQuantities::fMaxFSPiPlusEnergy [protected]
Float_t MadQuantities::fMaxFSPiZeroEnergy [protected]
Float_t MadQuantities::fMaxFSProtonEnergy [protected]
Int_t MadQuantities::fNumFSGeantino [protected]
Int_t MadQuantities::fNumFSKaon [protected]

Definition at line 46 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSKaon(), and SetStdHepVars().

Int_t MadQuantities::fNumFSNeutron [protected]

Definition at line 30 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSNeutron(), and SetStdHepVars().

Int_t MadQuantities::fNumFSPiMinus [protected]

Definition at line 38 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSPion(), and SetStdHepVars().

Int_t MadQuantities::fNumFSPiPlus [protected]

Definition at line 34 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSPion(), and SetStdHepVars().

Int_t MadQuantities::fNumFSPiZero [protected]

Definition at line 42 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSPiZero(), and SetStdHepVars().

Int_t MadQuantities::fNumFSProton [protected]

Definition at line 26 of file MadQuantities.h.

Referenced by MadAnalysis::CreatePAN(), MadQuantities(), NumFSProton(), and SetStdHepVars().

Float_t MadQuantities::fTotFSKaonEnergy [protected]
Float_t MadQuantities::fTotFSPiPlusEnergy [protected]
Float_t MadQuantities::fTotFSPiZeroEnergy [protected]
Float_t MadQuantities::fTotFSProtonEnergy [protected]

The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1