MadDpAnalysis Class Reference

#include <MadDpAnalysis.h>

Inheritance diagram for MadDpAnalysis:

List of all members.

Public Member Functions
	MadDpAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
	MadDpAnalysis (JobC *, string, int)
	~MadDpAnalysis ()
void	MakeMyFile (std::string, int, int)
void	ReadPIDFile (std::string)
void	CreatePAN (std::string, Int_t, Int_t)
Static Public Member Functions
bool	InFidVol (const Detector::Detector_t &det, const Float_t &x, const Float_t &y, const Float_t &z)
Protected Member Functions
Bool_t	PassTruthSignal (Int_t event=0)
Bool_t	PassAnalysisCuts (Int_t event=0)
Float_t	PID (Int_t event=0, Int_t method=0)
Bool_t	PassBasicCuts (Int_t event=0)
Bool_t	PassFid (Int_t event=0)
Bool_t	PassFidNew (Int_t event=0)
Bool_t	PassFidNewN (Int_t event=0)
Float_t	RecoAnalysisEnergy (Int_t event=0)
AnnInputBlock	AnnVar (Int_t event=0)
Double_t	GetAnnPid (AnnInputBlock anninput, Int_t det, Int_t tar, Int_t fid, Int_t prior, Bool_t first_ann)
Double_t	Sigmoid (Double_t x)
Float_t	SingleTimeFrame (Int_t snarlentry, Int_t nentries)
Protected Attributes
TFile *	fLikeliFile
TH1F *	fLikeliHist [6]
Private Attributes
Double_t	Ann_weight_vector [226]

---

Constructor & Destructor Documentation

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

Definition at line 27 of file MadDpAnalysis.cxx. References MadBase::Clear(), fLikeliFile, fLikeliHist, and MadBase::InitChain(). { if(!chainSR) { record = 0; strecord = 0; emrecord = 0; mcrecord = 0; threcord = 0; Clear(); whichSource = -1; isMC = true; isTH = true; isEM = true; Nentries = -1; return; } InitChain(chainSR,chainMC,chainTH,chainEM); whichSource = 0; fLikeliFile = NULL; for(int i=0;i<6;i++) fLikeliHist[i] = NULL; }

MadDpAnalysis::MadDpAnalysis (  JobC *  , string  , int  )

Definition at line 54 of file MadDpAnalysis.cxx. References MadBase::Clear(), fLikeliFile, and fLikeliHist. { record = 0; strecord = 0; emrecord = 0; mcrecord = 0; threcord = 0; Clear(); isMC = true; isTH = true; isEM = true; Nentries = entries; whichSource = 1; jcPath = path; fJC = j; fLikeliFile = NULL; for(int i=0;i<6;i++) fLikeliHist[i] = NULL; }

MadDpAnalysis::~MadDpAnalysis (   )

Definition at line 74 of file MadDpAnalysis.cxx. References fLikeliFile. { if(fLikeliFile) fLikeliFile->Close(); }

---

Member Function Documentation

AnnInputBlock MadDpAnalysis::AnnVar (  Int_t  event = 0  )  [protected]

Definition at line 172 of file MadDpAnalysis.cxx. References AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhcommon, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwdig, AnnInputBlock::aShwph, AnnInputBlock::aShwphper, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwphperpl, AnnInputBlock::aShwphperstp, AnnInputBlock::aShwpl, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aShwstp, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrklen, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkph, AnnInputBlock::aTrkphper, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, AnnInputBlock::aTrkstp, NtpSRPlane::beg, NtpSRPlane::end, NtpSRTrack::fit, VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower_Jim(), MadBase::LoadStrip(), NtpSRPlane::n, NtpSRShower::ndigit, NtpSRShower::nstrip, NtpSRTrack::nstrip, NtpSREvent::nstrip, NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRFitTrack::pass, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, NtpSRTrack::stp, NtpSRShower::stp, NtpSREvent::stp, NtpSRStrip::strip, striptime, NtpSRStrip::time0, and NtpSRStrip::time1. Referenced by CreatePAN(). { AnnInputBlock anninput; // Initialize anninput.aTotrk =0.; anninput.aTotstp =0.; anninput.aTotph =0.; anninput.aTotlen =0.; anninput.aPhperpl =0.; anninput.aPhperstp =0.; anninput.aTrkpass =0.; anninput.aTrkph =0.; anninput.aTrklen =0.; anninput.aTrkphperpl =0.; anninput.aTrkphper =0.; anninput.aTrkplu =0.; anninput.aTrkplv =0.; anninput.aTrkstp =0.; anninput.aShwph =0.; anninput.aShwstp =0.; anninput.aShwdig =0.; anninput.aShwpl =0.; anninput.aShwphper =0.; anninput.aShwphperpl =0.; anninput.aShwphperdig =0.; anninput.aShwphperstp =0.; anninput.aShwplu =0.; anninput.aShwplv =0.; anninput.aPh3 =0.; anninput.aPh6 =0.; anninput.aPhlast =0.; anninput.aPhcommon =0.; anninput.aTimemax =0.; anninput.aTimemin =0.; // // Calculate variables needed for ANN (and a few more) Int_t detector = 0; if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) detector=1; if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar) detector=2; LoadEvent(event) ; int track_index = -1; LoadLargestTrackFromEvent(event,track_index); int shower_index = -1; // LoadLargestShowerFromEvent(event,shower_index); LoadShower_Jim(event,shower_index,detector); anninput.aTotrk =ntpEvent->ntrack; anninput.aTotstp =ntpEvent->nstrip; anninput.aTotph =ntpEvent->ph.sigcor; anninput.aTotlen =ntpEvent->plane.end-ntpEvent->plane.beg+1; // prepei na to alla3w if (ntpEvent->plane.n>0) anninput.aPhperpl =ntpEvent->ph.sigcor/ntpEvent->plane.n; if (ntpEvent->nstrip>0) anninput.aPhperstp =ntpEvent->ph.sigcor/ntpEvent->nstrip; if(track_index!=-1) { anninput.aTrkpass =ntpTrack->fit.pass; anninput.aTrkph =ntpTrack->ph.sigcor; anninput.aTrklen =ntpTrack->plane.ntrklike; if(ntpTrack->plane.ntrklike>0) anninput.aTrkphperpl =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike; if(ntpEvent->ph.sigcor>0) anninput.aTrkphper =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor; anninput.aTrkplu =ntpTrack->plane.nu; anninput.aTrkplv =ntpTrack->plane.nv; anninput.aTrkstp =ntpTrack->nstrip; } if(shower_index!=-1) { anninput.aShwph =ntpShower->ph.sigcor; anninput.aShwstp =ntpShower->nstrip; anninput.aShwdig =ntpShower->ndigit; anninput.aShwpl =ntpShower->plane.end-ntpShower->plane.beg+1; if (ntpEvent->ph.sigcor>0) anninput.aShwphper =ntpShower->ph.sigcor/ntpEvent->ph.sigcor; anninput.aShwphperpl =ntpShower->ph.sigcor/(ntpShower->plane.end-ntpShower->plane.beg+1); if (ntpShower->ndigit>0) anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit; if (ntpShower->nstrip>0) anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip; anninput.aShwplu =ntpShower->plane.nu; anninput.aShwplv =ntpShower->plane.nv; } Int_t ssind,ssplind; Double_t ssphtot; Bool_t foundshw,foundtrk; Int_t planes=ntpEvent->plane.beg; Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere; ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.; // loop over strips to compute ph3 ph6 phlast phcommon Double_t timemin=9.*10e30; Double_t timemax=-9999999; for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); if(stp_index!=-1){ LoadStrip(stp_index); // TIMING INFO // // Find max min time of events by loading strips for ND if(detector==1){ Double_t striptime=ntpStrip->time1; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } if(detector==2){ Double_t striptime1=ntpStrip->time1; Double_t striptime0=ntpStrip->time0; Double_t striptime=0; if(striptime1>0 && striptime0<0) striptime=striptime1; if(striptime0>0 && striptime1<0) striptime=striptime0; if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } ssind=ntpStrip->strip; ssplind=ntpStrip->plane; ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; foundshw=false; foundtrk=false; Double_t phstrips=0; Double_t phstript=0; // shower strips Int_t sshwind,sshwplind; Double_t sshwphtot; if(shower_index!=-1) { for(Int_t jj=0;jj<ntpShower->nstrip;jj++){ Int_t stp_indexs=((ntpShower->stp)[jj]); if(stp_indexs!=-1){ LoadStrip(stp_indexs); if(!foundshw){ sshwind=ntpStrip->strip; sshwplind=ntpStrip->plane; sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; if(sshwind==ssind && sshwplind==ssplind) { foundshw=true; phstrips=sshwphtot; } } } } } // tracks strips Int_t strkind,strkplind; Double_t strkphtot; if(track_index!=-1) { for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){ Int_t stp_indext=((ntpTrack->stp)[jj]); if(stp_indext!=-1){ LoadStrip(stp_indext); if(!foundtrk){ strkind=ntpStrip->strip; strkplind=ntpStrip->plane; strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; if(strkind==ssind && strkplind==ssplind) { foundtrk=true; phstript=strkphtot; } } } } } if(foundshw && foundtrk) { hitcommon=hitcommon+1; phcommon=phcommon+phstrips+phstript; } if(!foundshw && ! foundtrk) { hitnowere=hitnowere+1; phnowere=phnowere+ssphtot; } if(ssplind>=planes && ssplind<=(planes+3)){ ph3=ph3+ssphtot; } else if(ssplind>(planes+3) && ssplind<=(planes+6)){ ph6=ph6+ssphtot; } else { phlast=phlast+ssphtot; } } // end if strip != -1 (!!!) } // end of looping over event strips.... anninput.aTimemin = timemin; anninput.aTimemax = timemax; // anninput.aPh3 =ph3; anninput.aPh6 =ph6; anninput.aPhlast =phlast; anninput.aPhcommon =phcommon; return anninput; }

void MadDpAnalysis::CreatePAN (  std::string  , Int_t  , Int_t  )

Definition at line 944 of file MadDpAnalysis.cxx. References abs(), AnnVar(), AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhcommon, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwdig, AnnInputBlock::aShwph, AnnInputBlock::aShwphper, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwphperpl, AnnInputBlock::aShwphperstp, AnnInputBlock::aShwpl, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aShwstp, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrklen, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkph, AnnInputBlock::aTrkphper, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, AnnInputBlock::aTrkstp, NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, NtpSRMomentum::best, NtpSRFitTrack::chi2, NtpSRTrack::contained, NtpSREventSummary::date, NtpSRDate::day, NtpSRTrack::ds, NtpSRTrack::end, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, BeamMonSpill::fHornCur, NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTrtgtd, BDSpillAccessor::Get(), GetAnnPid(), ScanList::GetDecision(), VldContext::GetDetector(), MadBase::GetEntry(), SpillInfoBlock::GetHorn(), SpillInfoBlock::GetHpos(), SpillInfoBlock::GetMagnet(), VldTimeStamp::GetNanoSec(), SpillTimeFinder::GetNearestSpill(), SpillInfoBlock::GetPot(), RecDataHeader::GetRun(), VldTimeStamp::GetSec(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), SpillInfo::GetSpillInfo(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), SpillInfoBlock::GetTgtpos(), SpillTimeFinder::GetTimeOfNearestSpill(), SpillTimeND::GetTimeStamp(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), SpillInfoBlock::GetVpos(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSREvent::index, MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), MadQuantities::IsFidVtxEvt(), NtpSREventSummary::litime, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), BDSpillAccessor::LoadSpill(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), NtpSRTrack::momentum, NtpSRDate::month, month, NtpSRPlane::n, NtpSRShower::ncluster, NtpSRTrack::ndigit, NtpSRShower::ndigit, NtpSREvent::ndigit, NtpSRFitTrack::ndof, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSREventSummary::nshower, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSREvent::nstrip, NtpSREvent::ntrack, NtpSRPlane::nu, MadQuantities::Nucleus(), NtpSRPlane::nv, NtpSRFitTrack::pass, PassAnalysisCuts(), PassBasicCuts(), passfail(), PassFid(), PassFidNew(), PassFidNewN(), NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSREventSummary::ph, PID(), NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, MadQuantities::Q2(), NtpSRMomentum::qp, NtpSRMomentum::range, NtpSRPulseHeight::raw, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), MadQuantities::RecoShwEnergySqrt(), run(), NtpSRPulseHeight::sigcor, BeamMonSpill::SpillTime(), MadQuantities::Target4Vector(), NtpMCFluxInfo::tptype, NtpMCFluxInfo::tpx, NtpMCFluxInfo::tpy, NtpMCFluxInfo::tpz, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSRVertex::u, NtpSRVertex::v, NtpSRTrack::vtx, NtpSRShower::vtx, NtpSREvent::vtx, MadQuantities::W2(), NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::Y(), NtpSRDate::year, NtpSRVertex::z, and SpillInfoBlock::Zero(). { // 1 = le-10 2 = pme 3 = phe there is no other way in an MC file to now what target positions it // was...! std::string savename = "PAN_" + tag + ".root"; TFile save(savename.c_str(),"RECREATE"); save.cd(); /* GnumiInterface gnumi; if(!gnumi.Status()) { cout << "Error MadAnalysis::CreatePAN() - Flux files not open." << " Will not be filling NuParent info" << endl; cout << "Set environmental variable $GNUMIAUX to point to the " << "directory containing the gnumi files to fix this!" << endl; } */ SpillInfo pppinfo; //ann quantities Float_t ann_aTotrk =0.; Float_t ann_aTotstp =0.; Float_t ann_aTotph =0.; Float_t ann_aTotlen =0.; Float_t ann_aPhperpl =0.; Float_t ann_aPhperstp =0.; Float_t ann_aTrkpass =0.; Float_t ann_aTrkph =0.; Float_t ann_aTrklen =0.; Float_t ann_aTrkphperpl =0.; Float_t ann_aTrkphper =0.; Float_t ann_aTrkplu =0.; Float_t ann_aTrkplv =0.; Float_t ann_aTrkstp =0.; Float_t ann_aShwph =0.; Float_t ann_aShwstp =0.; Float_t ann_aShwdig =0.; Float_t ann_aShwpl =0.; Float_t ann_aShwphper =0.; Float_t ann_aShwphperpl =0.; Float_t ann_aShwphperdig =0.; Float_t ann_aShwphperstp =0.; Float_t ann_aShwplu =0.; Float_t ann_aShwplv =0.; Float_t ann_aPh3 =0.; Float_t ann_aPh6 =0.; Float_t ann_aPhlast =0.; Float_t ann_aPhcommon =0.; //PAN Quantities //Truth: Float_t true_enu; // true neutrino energy (GeV) Float_t true_pxnu; // true neutrino momentum-x (GeV) Float_t true_pynu; // true neutrino momentum-y (GeV) Float_t true_pznu; // true neutrino momentum-z (GeV) Float_t true_etgt; // true target energy (GeV) Float_t true_pxtgt; // true target momentum-x (GeV) Float_t true_pytgt; // true target momentum-y (GeV) Float_t true_pztgt; // true target momentum-z (GeV) Float_t true_emu; // true muon energy Float_t true_eshw; // true shower energy Float_t true_x; // true x Float_t true_y; // true y Float_t true_q2; // true q2 Float_t true_w2; // true w2 Float_t true_dircosneu; // true muon dircos w.r.t neutrino Float_t true_dircosz; // track z-direction cosine Float_t true_vtxx; // true x vtx Float_t true_vtxy; // true y vtx Float_t true_vtxz; // true z vtx //For Neugen: Int_t flavour; // true flavour: 1-e 2-mu 3-tau Int_t process; // process: 1001-QEL 1002-RES 1003-DIS 1004-COH Int_t nucleus; // target nucleus: 274-C 284-O 372-Fe Int_t cc_nc; // cc/nc flag: 1-cc 2-nc Int_t initial_state; // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ... Int_t had_fs; // hadronic final state, number between 200-300 //For Beam Reweighting: // NuParent *nuparent = new NuParent(); //Reco Quantities Int_t detector; // Near = 1, Far = 2; Int_t run; // run number Int_t snarl; // snarl number Int_t nevent; // number of events in snarl Int_t event; // event index Int_t subrun; // subrun number Int_t trigsrc; // trigger source Int_t mcevent; // mc event index Int_t ntrack; // number of tracks in event Int_t nshower; // number of showers in event Int_t is_fid; // pass fid cut. true = 1, false = 0 Int_t is_fid_dp; // pass dp analysis fid cut. true = 1, false = 0 Int_t is_fid_ns; // pass nd analysis fid cut true=1, false=0 // !! NEW Int_t is_cev; // fully contained. true = 1, false = 0 Int_t is_cont_trk; // new containemnt support by IsContained() for the largest track Int_t is_mc; // is a corresponding mc event found Int_t pass_basic; // Pass basic cuts. true = 1, false = 0 Float_t pid0; // pid parameter (usual method) Float_t pid1; // pid parameter (alternative method 1) Float_t pid2; // pid parameter (alternative method 2) Float_t annpid; // ANN PID Float_t annpid_f1p1; // ANN PID FAR FIDUCIAL DP PRIOR 1 // UN-OSICLLATED Float_t annpid_f1p2; // ANN PID FAR FIDUCIAL DP PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95 Float_t annpid_f1p3; // ANN PID FAR FIDUCIAL DP PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95 Float_t annpid_f2p1; // ANN PID FAR FIDUCIAL NS PRIOR 1 UN-OSICLLATED Float_t annpid_f2p2; // ANN PID FAR FIDUCIAL NS PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95 Float_t annpid_f2p3; // ANN PID FAR FIDUCIAL NS PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95 Int_t pass; // pass analysis cuts. true = 1, false = 0 Float_t reco_enu; // reco neutrino energy Float_t reco_emu; // reco muon energy Float_t reco_eshw; // reco shower energy (shw.ph.gev/1.23) Int_t pass_fd_qualcuts; // pass FD quality cuts (HV trips etc) Int_t litag; Float_t reco_eshwcc; // linear cc version Float_t reco_eshwnc; // linear nc version Float_t reco_eshwccw; // weighted cc version Float_t reco_eshwncw; // weighted nc version Float_t reco_eshw_sqrt; // reco shower energy using sqrt method Float_t reco_qe_enu; // reco qe neutrino energy Float_t reco_qe_q2; // reco qe q2 Float_t reco_y; // reco y Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino Float_t reco_dircosz; // reco track vtx z-dircos Float_t reco_vtxx; // reco x vtx Float_t reco_vtxy; // reco y vtx Float_t reco_vtxz; // reco z vtx Float_t evtlength; // reco event length Float_t evtph; // reco event ph Float_t trklength; // reco track length Float_t trkmomrange; // reco track momentum from range Float_t trkqp; // reco track q/p Float_t trkeqp; // reco track q/p error Float_t trkmombest; // reco best track momentum based on new containment support Float_t trkphfrac; // reco pulse height fraction in track Float_t trkphplane; // reco average track pulse height/plane Float_t trkds; // track DS Float_t rawph; // Raw ph // Additional track info Float_t trkendz; // track end Z position Float_t trkendx; // track end X position Float_t trkendy; // track end Y position Float_t trkendu; // track end U position Float_t trkendv; // track end V position Float_t trkvtxz; // track begin Z position Float_t trkvtxx; // track begin X position Float_t trkvtxy; // track begin Y position Float_t trkvtxu; // track begin U position Float_t trkvtxv; // track begin V position Float_t trkplbu; // track begin plane u Float_t trkplbv; // track begin plane v Float_t trkpleu; // track end plane u Float_t trkplev; // track end plane v Float_t trkple; // track end plane Float_t trkplb; // track begin plane Int_t trkfitpass; // Int_t totdigits,totstrips; Int_t trkdigits,trkstrips; Float_t trkph,trkchi2,trkndof; Int_t trkdiffuv; // shower variables Float_t shwph; Int_t shwdigits, shwstrips; Float_t shwvtxz; // shower begin Z position Float_t shwvtxx; // shower begin X position Float_t shwvtxy; // shower begin Y position Float_t shwvtxu; // shower begin U position Float_t shwvtxv; // shower begin V position // cluster variables Int_t shwncluster; // SCAN DECISION // <- define an instance of ScanList class here Int_t scandecision; ScanList scan; // EVENT TIMING Double_t evttimemin; //MIN STRIP TIME OF EVENT Double_t evttimemax; //MAX STRIP TIME OF EVENT Double_t snarlt; Double_t litime; // BEAM INFO Float_t pot,horn,tar,vvpos,hhpos,magn; // Beam INFO DATABASE Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb; Double_t timedb; // TIME INFO DATABASE Double_t neartdb,fartdb,neardifftdb; // mc flux info for xf,pt weighting Float_t mcflux_tpx; Float_t mcflux_tpy; Float_t mcflux_tpz; Float_t mcflux_tptype; // IS SNARL CUT IN PIECES Float_t singletimeframe; // Beam Info Block SpillInfoBlock *spinfo = new SpillInfoBlock(); spinfo->Zero(); singletimeframe=-999; // Int_t day,month,year; Int_t pass_beamcuts; // pass default beam quality cuts //Trees TTree *tree = new TTree("pan","pan"); //Truth tree->Branch("true_enu",&true_enu,"true_enu/F",512000); tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000); tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000); tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000); tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000); tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000); tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000); tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000); tree->Branch("true_emu",&true_emu,"true_emu/F",512000); tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000); tree->Branch("true_x",&true_x,"true_x/F",512000); tree->Branch("true_y",&true_y,"true_y/F",512000); tree->Branch("true_q2",&true_q2,"true_q2/F",512000); tree->Branch("true_w2",&true_w2,"true_w2/F",512000); tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000); tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000); tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000); tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000); tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000); //Neugen tree->Branch("flavour",&flavour,"flavour/I",512000); tree->Branch("process",&process,"process/I",512000); tree->Branch("nucleus",&nucleus,"nucleus/I",512000); tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000); tree->Branch("initial_state",&initial_state,"initial_state/I",512000); tree->Branch("had_fs",&had_fs,"had_fs/I",512000); //NuParent //tree->Branch("nuparent","NuParent",&nuparent,8000,1); //Reco tree->Branch("detector",&detector,"detector/I",512000); tree->Branch("run",&run,"run/I",512000); tree->Branch("snarl",&snarl,"snarl/I",512000); tree->Branch("event",&event,"event/I",512000); tree->Branch("mcevent",&mcevent,"mcevent/I",512000); tree->Branch("ntrack",&ntrack,"ntrack/I",512000); tree->Branch("nshower",&nshower,"nshower/I",512000); tree->Branch("subrun",&subrun,"subrun/I",512000); tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000); tree->Branch("litime",&litime,"litime/D",512000); tree->Branch("is_fid",&is_fid,"is_fid/I",512000); tree->Branch("is_fid_dp",&is_fid_dp,"is_fid_dp/I",512000); tree->Branch("is_fid_ns",&is_fid_ns,"is_fid_ns/I",512000); tree->Branch("pass_fd_qualcuts", &pass_fd_qualcuts, "pass_fd_qualcuts/I"); tree->Branch("litag", &litag, "litag/I"); tree->Branch("is_cev",&is_cev,"is_cev/I",512000); tree->Branch("is_cont_trk",&is_cont_trk,"is_cont_trk/I",512000); tree->Branch("is_mc",&is_mc,"is_mc/I",512000); tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000); tree->Branch("pid0",&pid0,"pid0/F",512000); tree->Branch("pid1",&pid1,"pid1/F",512000); tree->Branch("pid2",&pid2,"pid2/F",512000); tree->Branch("annpid",&annpid,"annpid/F",512000); tree->Branch("annpid_f1p1",&annpid_f1p1,"annpid_f1p1/F",512000); tree->Branch("annpid_f1p2",&annpid_f1p2,"annpid_f1p2/F",512000); tree->Branch("annpid_f1p3",&annpid_f1p3,"annpid_f1p3/F",512000); tree->Branch("annpid_f2p1",&annpid_f2p1,"annpid_f2p1/F",512000); tree->Branch("annpid_f2p2",&annpid_f2p2,"annpid_f2p2/F",512000); tree->Branch("annpid_f2p3",&annpid_f2p3,"annpid_f2p3/F",512000); tree->Branch("pass",&pass,"pass/I",512000); tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000); tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000); tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000); tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000); tree->Branch("reco_eshwcc", &reco_eshwcc, "reco_eshwcc/F",512000); tree->Branch("reco_eshwnc", &reco_eshwnc, "reco_eshwnc/F",512000); tree->Branch("reco_eshwccw", &reco_eshwccw, "reco_eshwccw/F",512000); tree->Branch("reco_eshwncw", &reco_eshwncw, "reco_eshwncw/F",512000); tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000); tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000); tree->Branch("reco_y",&reco_y,"reco_y/F",512000); tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000); tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000); tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000); tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000); tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000); tree->Branch("evtlength",&evtlength,"evtlength/F",512000); tree->Branch("evtph",&evtph,"evtph/F",512000); tree->Branch("trklength",&trklength,"trklength/F",512000); tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000); tree->Branch("trkqp",&trkqp,"trkqp/F",512000); tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000); tree->Branch("trkmombest",&trkmombest,"trkmombest/F",512000); tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000); tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000); tree->Branch("rawph",&rawph,"rawph/F",512000); tree->Branch("trkendz",&trkendz,"trkendz/F",512000); tree->Branch("trkendu",&trkendu,"trkendu/F",512000); tree->Branch("trkendv",&trkendv,"trkendv/F",512000); tree->Branch("trkendx",&trkendx,"trkendx/F",512000); tree->Branch("trkendy",&trkendy,"trkendy/F",512000); tree->Branch("trkplbu",&trkplbu,"trkplbu/F",512000); tree->Branch("trkplbv",&trkplbv,"trkplbv/F",512000); tree->Branch("trkpleu",&trkpleu,"trkpleu/F",512000); tree->Branch("trkplev",&trkplev,"trkplev/F",512000); tree->Branch("trkple",&trkple,"trkple/F",512000); tree->Branch("trkplb",&trkplb,"trkplb/F",512000); tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000); tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000); tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000); tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000); tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000); tree->Branch("trkds", &trkds,"trkds/F",512000); tree->Branch("evttimemin",&evttimemin,"evttimemin/D"); tree->Branch("evttimemax",&evttimemax,"evttimemax/D"); tree->Branch("totdigits",&totdigits,"totdigits/I"); tree->Branch("totstrips",&totstrips,"totstrips/I"); tree->Branch("trkdigits",&trkdigits,"trkdigits/I"); tree->Branch("trkstrips",&trkstrips,"trkstrips/I"); tree->Branch("trkph",&trkph,"trkph/F"); tree->Branch("trkdiffuv",&trkdiffuv,"trkdiffuv/I",512000); tree->Branch("trkchi2",&trkchi2,"trkchi2/F"); tree->Branch("trkndof",&trkndof,"trkndof/F"); tree->Branch("trkfitpass",&trkfitpass,"trkfitpass/I",512000); tree->Branch("shwdigits",&shwdigits,"shwdigits/I"); tree->Branch("shwstrips",&shwstrips,"shwstrips/I"); tree->Branch("shwph",&shwph,"shwph/F"); tree->Branch("shwvtxz",&shwvtxz,"shwvtxz/F",512000); tree->Branch("shwvtxu",&shwvtxu,"shwvtxu/F",512000); tree->Branch("shwvtxv",&shwvtxv,"shwvtxv/F",512000); tree->Branch("shwvtxx",&shwvtxx,"shwvtxx/F",512000); tree->Branch("shwvtxy",&shwvtxy,"shwvtxy/F",512000); tree->Branch("shwncluster",&shwncluster,"shwncluster/I"); // SCAN DECISIONS tree->Branch("scandecision",&scandecision,"scandecision/I"); //BEAM INFO tree->Branch("pot", &pot, "pot/F",512000); tree->Branch("horn", &horn, "horn/F",512000); tree->Branch("tar", &tar, "tar/F",512000); tree->Branch("vvpos", &vvpos, "vvpos/F",512000); tree->Branch("hhpos", &hhpos, "hhpos/F",512000); tree->Branch("magn", &magn, "magn/F",512000); tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000); // BEAM INFO DB tree->Branch("potdb", &potdb, "potdb/F",512000); tree->Branch("horndb", &horndb, "horndb/F",512000); tree->Branch("tardb", &tardb, "tardb/F",512000); tree->Branch("vvposdb", &vvposdb, "vvposdb/F",512000); tree->Branch("hhposdb", &hhposdb, "hhposdb/F",512000); tree->Branch("vvwidthdb", &vvwidthdb, "vvwidthdb/F",512000); tree->Branch("hhwidthdb", &hhwidthdb, "hhwidthdb/F",512000); tree->Branch("timedb", &timedb, "timedb/D"); // TIME DB tree->Branch("neartdb", &neartdb, "neartdb/D"); tree->Branch("fartdb", &fartdb, "fartdb/D"); tree->Branch("neardifftdb", &neardifftdb,"neardifftdb/D"); tree->Branch("snarlt", &snarlt, "snarlt/D"); // mc flux info tree->Branch("mcflux_tpx", &mcflux_tpx, "mcflux_tpx/F"); tree->Branch("mcflux_tpy", &mcflux_tpy, "mcflux_tpy/F"); tree->Branch("mcflux_tpz", &mcflux_tpz, "mcflux_tpz/F"); tree->Branch("mcflux_tptype", &mcflux_tptype, "mcflux_tptype/F"); // ann vars tree->Branch("ann_aTotrk", &ann_aTotrk, "ann_aTotrk/F"); tree->Branch("ann_aTotstp", &ann_aTotstp, "ann_aTotstp/F"); tree->Branch("ann_aTotph", &ann_aTotph, "ann_aTotph/F"); tree->Branch("ann_aTotlen", &ann_aTotlen, "ann_aTotlen/F"); tree->Branch("ann_aPhperpl", &ann_aPhperpl, "ann_aPhperpl/F"); tree->Branch("ann_aPhperstp", &ann_aPhperstp, "ann_aPhperstp/F"); tree->Branch("ann_aTrkpass", &ann_aTrkpass, "ann_aTrkpass/F"); tree->Branch("ann_aTrkph", &ann_aTrkph, "ann_aTrkph/F"); tree->Branch("ann_aTrklen", &ann_aTrklen, "ann_aTrklen/F"); tree->Branch("ann_aTrkphperpl", &ann_aTrkphperpl, "ann_aTrkphperpl/F"); tree->Branch("ann_aTrkphper", &ann_aTrkphper, "ann_aTrkphper/F"); tree->Branch("ann_aTrkplu", &ann_aTrkplu, "ann_aTrkplu/F"); tree->Branch("ann_aTrkplv", &ann_aTrkplv, "ann_aTrkplv/F"); tree->Branch("ann_aTrkstp", &ann_aTrkstp, "ann_aTrkstp/F"); tree->Branch("ann_aShwph", &ann_aShwph, "ann_aShwph/F"); tree->Branch("ann_aShwstp", &ann_aShwstp, "ann_aShwstp/F"); tree->Branch("ann_aShwdig", &ann_aShwdig, "ann_aShwdig/F"); tree->Branch("ann_aShwpl", &ann_aShwpl, "ann_aShwpl/F"); tree->Branch("ann_aShwphper", &ann_aShwphper, "ann_aShwphper/F"); tree->Branch("ann_aShwphperpl", &ann_aShwphperpl, "ann_aShwphperpl/F"); tree->Branch("ann_aShwphperdig", &ann_aShwphperdig,"ann_aShwphperdig/F"); tree->Branch("ann_aShwphperstp", &ann_aShwphperstp,"ann_aShwphperstp/F"); tree->Branch("ann_aShwplu", &ann_aShwplu, "ann_aShwplu/F"); tree->Branch("ann_aShwplv", &ann_aShwplv, "ann_aShwplv/F"); tree->Branch("ann_aPh3", &ann_aPh3, "ann_aPh3/F"); tree->Branch("ann_aPh6", &ann_aPh6, "ann_aPh6/F"); tree->Branch("ann_aPhlast", &ann_aPhlast, "ann_aPhlast/F"); tree->Branch("ann_aPhcommon", &ann_aPhcommon, "ann_aPhcommon/F"); tree->Branch("day", &day, "day/I"); tree->Branch("month", &month, "month/I"); tree->Branch("year", &year, "year/I"); tree->Branch("pass_beamcuts", &pass_beamcuts, "pass_beamcuts/I"); Bool_t first_ann = true; for(int i=0;i<Nentries;i++){ if(i%400==0) std::cout << float(i)*100./float(Nentries) << "% done" << std::endl; if(!this->GetEntry(i)) continue; nevent = eventSummary->nevent; if(nevent==0) continue; run = ntpHeader->GetRun(); snarl = ntpHeader->GetSnarl(); subrun = ntpHeader->GetSubRun(); trigsrc = ntpHeader->GetTrigSrc(); day = eventSummary->date.day; month = eventSummary->date.month; year = eventSummary->date.year; Double_t snarltime =ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9); Int_t month = eventSummary->date.month; litime = eventSummary->litime; pass_fd_qualcuts =passfail(snarltime); litag=0; Int_t numshower=eventSummary->nshower; Float_t allph=eventSummary->ph.raw; Int_t plbeg=eventSummary->plane.beg; Int_t plend=eventSummary->plane.end; if (numshower) { if (allph>1e6) litag+=10; if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) || ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) || ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) || ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++; if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) || ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) || ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) || ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++; } // SCAN INFO scandecision=scan.GetDecision(run,snarl); // only pass through scanned snarls if scanfilter=1 if (scanfilter==0 || scandecision>0) { // Get Beam Monitorint Info if NOT MC: if(ntpHeader->GetVldContext().GetSimFlag()!=4){ pppinfo.GetSpillInfo(month,year,snarltime,*spinfo); pot =spinfo->GetPot(); horn =spinfo->GetHorn(); hhpos =spinfo->GetHpos(); vvpos =spinfo->GetVpos(); tar =spinfo->GetTgtpos(); magn =spinfo->GetMagnet(); potdb = -999; horndb = -999; tardb = -999; vvposdb = -999; hhposdb = -999; vvwidthdb = -999; hhwidthdb = -999; timedb = -999; neartdb = -999; fartdb = -999; neardifftdb = -999; // BEAM INFO DB Detector::Detector_t detectort; SimFlag::SimFlag_t simflag; VldTimeStamp timestamp; VldContext cx; cx = ntpHeader->GetVldContext(); detectort = ntpHeader->GetVldContext().GetDetector(); simflag = ntpHeader->GetVldContext().GetSimFlag(); timestamp = ntpHeader->GetVldContext().GetTimeStamp(); BDSpillAccessor &spillAccess=BDSpillAccessor::Get(); const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp); if(spillInfoDB) { if(spillInfoDB->fTortgt !=0.) potdb = spillInfoDB->fTortgt; else if(spillInfoDB->fTrtgtd !=0.) potdb = spillInfoDB->fTrtgtd; else if(spillInfoDB->fTor101 !=0.) potdb = spillInfoDB->fTor101; else potdb = -999; horndb = spillInfoDB->fHornCur; tardb = (int)spillInfoDB->GetStatusBits().target_in; vvposdb = spillInfoDB->fTargBpmY[0]; hhposdb = spillInfoDB->fTargBpmX[0]; vvwidthdb = spillInfoDB->fProfWidY; hhwidthdb = spillInfoDB->fProfWidX; timedb = spillInfoDB->SpillTime(); // apply beam quality cuts if (potdb>0.5 && fabs(snarltime-timedb)<1 && (hhposdb<0 && hhposdb>-0.002) && (vvposdb>0 && vvposdb<0.002) && (hhwidthdb>0.0001 && hhwidthdb<0.0022) && (vvwidthdb>0.0001 && vvwidthdb<0.005) && (horndb<-155 && horndb>-200)) pass_beamcuts=1; } // TIME DB SpillTimeFinder& spillfinder = SpillTimeFinder::Instance(); fartdb = spillfinder.GetTimeOfNearestSpill(cx); const SpillTimeND& spnd = spillfinder.GetNearestSpill(cx); neartdb = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9); neardifftdb = spillfinder.GetTimeToNearestSpill(cx); snarlt = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9); } //singletimeframe=SingleTimeFrame(i,Nentries); Int_t trkcount=0; Int_t shwcount=0; Float_t totph=0; Int_t evtmin=0; Int_t evtmax=nevent; if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) { // cout << "MULTIPLE events!!\n"; Int_t maxevent=0; Float_t maxph=0; for(int j=0;j<nevent;j++){ if(!LoadEvent(j)) continue; // cout << "evtlength=" << ntpEvent->plane.n << " totph=" // << ntpEvent->ph.sigcor << "ntrack=" << ntpEvent->ntrack // << " nshower=" << ntpEvent->nshower << endl; Float_t evtpheight=ntpEvent->ph.sigcor; if (evtpheight>maxph) { maxph=evtpheight; maxevent=j; } } // cout << "BIGGEST event=" << maxevent << endl; evtmin=maxevent; evtmax=maxevent+1; } // EVENT LOOP - fill tree once for each reconstructed event for(int j=evtmin;j<evtmax;j++){ if(!LoadEvent(j)) continue; //no event found totph = 0; //set event index: event = j; ntrack = ntpEvent->ntrack; nshower = ntpEvent->nshower; totdigits=ntpEvent->ndigit; totstrips=ntpEvent->nstrip; // Initialize everything //nuparent->Zero(); mcflux_tpx=0; mcflux_tpy=0; mcflux_tpz=0; mcflux_tptype=0; true_enu = 0; true_emu = 0; true_eshw = 0; true_pxnu = 0; true_pynu = 0; true_pznu = 0; true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0; true_dircosneu = 0; true_dircosz = 0; true_vtxx = 0; true_vtxy = 0; true_vtxz = 0; flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0; had_fs = 0; true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0; detector = -1; mcevent = -1; is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; pid0 = -999; pid1 = -999; pid2 = -999; pass = 0; reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0; reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0; reco_dircosz = 0; reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999; evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1; trkmomrange = -999; trkqp = -999; trkeqp = -999; trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.; trkendz=999; trkendu=999; trkendv=999; trkendx=999; trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999; trkds=-1; evttimemin=-1; evttimemax=-1; evtph=0.; annpid=-999.; annpid_f1p1=-999.;annpid_f1p2=-999.;annpid_f1p3=-999.; annpid_f2p1=-999.;annpid_f2p2=-999.;annpid_f2p3=-999.; reco_eshwcc=0.;reco_eshwnc=0.;reco_eshwccw=0.;reco_eshwncw=0.; is_fid_dp=0;is_fid_ns=0; is_cont_trk=0; trkmombest=-999.; trkfitpass=0; trkdiffuv=0; trkdigits=0; trkstrips=0; trkph=0; trkndof=0; trkchi2=0; shwdigits=0; shwstrips=0; shwph=0; shwvtxu=0; shwvtxv=0; shwvtxx=0; shwvtxy=0; shwvtxz=0; shwncluster=0; ann_aTotrk =0.; ann_aTotstp =0.; ann_aTotph =0.; ann_aTotlen =0.; ann_aPhperpl =0.; ann_aPhperstp =0.; ann_aTrkpass =0.; ann_aTrkph =0.; ann_aTrklen =0.; ann_aTrkphperpl =0.; ann_aTrkphper =0.; ann_aTrkplu =0.; ann_aTrkplv =0.; ann_aTrkstp =0.; ann_aShwph =0.; ann_aShwstp =0.; ann_aShwdig =0.; ann_aShwpl =0.; ann_aShwphper =0.; ann_aShwphperpl =0.; ann_aShwphperdig =0.; ann_aShwphperstp =0.; ann_aShwplu =0.; ann_aShwplv =0.; ann_aPh3 =0.; ann_aPh6 =0.; ann_aPhlast =0.; ann_aPhcommon =0.; rawph=ntpEvent->ph.raw; //get detector type: if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ detector = 2; totph=eventSummary->ph.sigcor; //fiducial criteria on vtx for far detector is_fid = 1; if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0; } else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ detector = 1; //get total charge in trk/shw for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) { LoadTrack(ii); totph+=ntpTrack->ph.sigcor; } trkcount+=ntrack; for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) { LoadShower(ii); totph+=ntpShower->ph.sigcor; } shwcount+=nshower; //fiducial criteria on vtx for near detector is_fid = 1; if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0; } if (detector==2) is_fid_dp=PassFidNew(ntpEvent->index); else is_fid_dp=PassFid(ntpEvent->index); if (detector==2) is_fid_ns=PassFidNewN(ntpEvent->index); else is_fid_ns=PassFid(ntpEvent->index); //check for a corresponding mc event if(ntpHeader->GetVldContext().GetSimFlag()==4){ if(LoadTruthForRecoTH(j,mcevent)) { Float_t *vtx = TrueVtx(mcevent); Float_t *nu3mom = TrueNuMom(mcevent); Float_t *targ4vec = Target4Vector(mcevent); //true info for tree: is_mc = 1; nucleus = Nucleus(mcevent); flavour = Flavour(mcevent); initial_state = Initial_state(mcevent); had_fs = HadronicFinalState(mcevent); process = IResonance(mcevent); cc_nc = IAction(mcevent); true_enu = TrueNuEnergy(mcevent); true_pxnu = nu3mom[0]; true_pynu = nu3mom[1]; true_pznu = nu3mom[2]; true_emu = TrueMuEnergy(mcevent); true_eshw = TrueShwEnergy(mcevent); true_x = X(mcevent); true_y = Y(mcevent); true_q2 = Q2(mcevent); true_w2 = W2(mcevent); true_dircosz = TrueMuDCosZ(mcevent); true_dircosneu = TrueMuDCosNeu(mcevent); true_vtxx = vtx[0]; true_vtxy = vtx[1]; true_vtxz = vtx[2]; true_etgt = targ4vec[3]; true_pxtgt = targ4vec[0]; true_pytgt = targ4vec[1]; true_pztgt = targ4vec[2]; /* if(gnumi.Status()) { if(isST) gnumi.GetParent(strecord,mcevent,*nuparent); else gnumi.GetParent(mcrecord,mcevent,*nuparent); } */ if (LoadTruth(mcevent)) { mcflux_tpx=ntpTruth->flux.tpx; mcflux_tpy=ntpTruth->flux.tpy; mcflux_tpz=ntpTruth->flux.tpz; mcflux_tptype=ntpTruth->flux.tptype; } delete [] vtx; delete [] nu3mom; delete [] targ4vec; } } if(PassBasicCuts(event)) {pass_basic=1;} // make this reflect pass/fail flag //reco info for tree: reco_vtxx = ntpEvent->vtx.x; reco_vtxy = ntpEvent->vtx.y; reco_vtxz = ntpEvent->vtx.z; evtlength = ntpEvent->plane.end-ntpEvent->plane.beg+1; evtph = ntpEvent->ph.sigcor; // different definition for neardet - accounts for // uninstrumented planes AnnInputBlock anninput=AnnVar(event); ann_aTotrk = anninput.aTotrk; ann_aTotstp = anninput.aTotstp; ann_aTotph = anninput.aTotph; ann_aTotlen = anninput.aTotlen; ann_aPhperpl = anninput.aPhperpl; ann_aPhperstp = anninput.aPhperstp; ann_aTrkpass = anninput.aTrkpass; ann_aTrkph = anninput.aTrkph; ann_aTrklen = anninput.aTrklen; ann_aTrkphperpl = anninput.aTrkphperpl; ann_aTrkphper = anninput.aTrkphper; ann_aTrkplu = anninput.aTrkplu; ann_aTrkplv = anninput.aTrkplv; ann_aTrkstp = anninput.aTrkstp; ann_aShwph = anninput.aShwph; ann_aShwstp = anninput.aShwstp; ann_aShwdig = anninput.aShwdig; ann_aShwpl = anninput.aShwpl; ann_aShwphper = anninput.aShwphper; ann_aShwphperpl = anninput.aShwphperpl; ann_aShwphperdig = anninput.aShwphperdig ; ann_aShwphperstp = anninput.aShwphperstp; ann_aShwplu = anninput.aShwplu; ann_aShwplv = anninput.aShwplv; ann_aPh3 = anninput.aPh3; ann_aPh6 = anninput.aPh6; ann_aPhlast = anninput.aPhlast; ann_aPhcommon = anninput.aPhcommon; // unused: Double_t trgtime=eventSummary->trigtime; // evttimemin = anninput.aTimemin-trgtime; // evttimemax = anninput.aTimemax-trgtime; // cedar - use snarl time rather than trigger time evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9; evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9; // ANN PID // Int_t target=1; if(ntpHeader->GetVldContext().GetSimFlag()!=4){ //if(tar<5000. && tar>3500.) target=1; //if(tar>5000. && tar<40000.) target=2; //if(tar>40000. && tar<100000.) target=3; // Used to work just fine for the pre shutdown data // after the shutdown this variable was not properly filled from ACNET data. in fact it is not // recoverable since it is corrupted in ACNET data. if(mcneartype==1) target=1; if(mcneartype==2) target=2; if(mcneartype==3) target=3; } if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){ if(mcneartype==1) target=1; if(mcneartype==2) target=2; if(mcneartype==3) target=3; } if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){ target=0; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC } // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING Int_t fid =1; Int_t prior = 1; if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ if(!PassFid(event)) annpid=-999; else { annpid=GetAnnPid(anninput,detector,target,fid,prior,first_ann); first_ann=false; } } if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ if(!is_fid_ns){ annpid_f2p1=-999; annpid_f2p2=-999; annpid_f2p3=-999; } if(is_fid_ns){ fid=2; prior=1; annpid_f2p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann); prior=2; annpid_f2p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann); prior=3; annpid_f2p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann); first_ann=false; } if(!is_fid_dp){ annpid_f1p1=-999; annpid_f1p2=-999; annpid_f1p3=-999; } if(is_fid_dp){ fid=1; prior=1; annpid_f1p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann); prior=2; annpid_f1p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann); prior=3; annpid_f1p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann); first_ann=false; } } // MAK March 8, 2005 // data or MC? // needed by energy corrections below bool isdata=true; if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false; //index will -1 if no track/shower present in event int track_index = -1; LoadLargestTrackFromEvent(j,track_index); int shower_index = -1; //LoadLargestShowerFromEvent(j,shower_index); // Change LoadShower_Jim(j,shower_index,detector); if(shower_index==-1) nshower=0; if(shower_index!=-1){ //check shower is present before using ntpShower shwdigits = ntpShower->ndigit; shwstrips = ntpShower->nstrip; shwph = ntpShower->ph.sigcor; shwvtxu = ntpShower->vtx.u; shwvtxv = ntpShower->vtx.v; shwvtxx = ntpShower->vtx.x; shwvtxy = ntpShower->vtx.y; shwvtxz = ntpShower->vtx.z; shwncluster = ntpShower->ncluster; // MAK March 8, 2005 Int_t shwtype=0; const Int_t shw_correction_mode=1; // Niki's re-tuning shwtype=-1; reco_eshw = RecoShwEnergy(shower_index,shwtype, detector,isdata, shw_correction_mode); shwtype=0; reco_eshwcc = RecoShwEnergy(shower_index,shwtype, detector,isdata, shw_correction_mode); shwtype=2; reco_eshwnc = RecoShwEnergy(shower_index,shwtype, detector,isdata, shw_correction_mode); shwtype=1; reco_eshwccw = RecoShwEnergy(shower_index,shwtype, detector,isdata, shw_correction_mode); shwtype=3; reco_eshwncw = RecoShwEnergy(shower_index,shwtype, detector,isdata, shw_correction_mode); reco_eshw_sqrt = RecoShwEnergySqrt(event); } if(track_index!=-1){ //check track is present before using ntpTrack // MAK March 8, 2005 reco_emu = RecoMuEnergy(0,track_index,isdata, ntpHeader->GetVldContext().GetDetector()); trklength = ntpTrack->plane.end-ntpTrack->plane.beg+1; trkmomrange = ntpTrack->momentum.range; trkqp = ntpTrack->momentum.qp; trkeqp = ntpTrack->momentum.eqp; trkendz = ntpTrack->end.z; trkendu = ntpTrack->end.u; trkendv = ntpTrack->end.v; trkendx = ntpTrack->end.x; trkendy = ntpTrack->end.y; trkvtxz = ntpTrack->vtx.z; trkvtxu = ntpTrack->vtx.u; trkvtxv = ntpTrack->vtx.v; trkvtxx = ntpTrack->vtx.x; trkvtxy = ntpTrack->vtx.y; trkds = ntpTrack->ds; trkplbu = ntpTrack->plane.begu; trkplbv = ntpTrack->plane.begv; trkpleu = ntpTrack->plane.endu; trkplev = ntpTrack->plane.endv; trkple = ntpTrack->plane.end; trkplb = ntpTrack->plane.beg; trkfitpass = ntpTrack->fit.pass; trkdiffuv = abs(ntpTrack->plane.nu-ntpTrack->plane.nv); trkdigits = ntpTrack->ndigit; trkstrips = ntpTrack->nstrip; trkph = ntpTrack->ph.sigcor; trkndof = ntpTrack->fit.ndof; trkchi2 = ntpTrack->fit.chi2; Float_t phtrack =ntpTrack->ph.sigcor; Float_t phevent =ntpEvent->ph.sigcor; /* In case of Birch // apply 1.8% correction to pulse height variables if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData ){ phtrack=phtrack*1.018; phevent=phevent*1.018; } */ if (phevent>0) {trkphfrac=phtrack/phevent;} if(ntpTrack->plane.n>0) { trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n; } is_cont_trk = ntpTrack->contained; trkmombest =ntpTrack->momentum.best; } else { trklength = 0; trkmomrange = 0; trkqp = 0; trkeqp = 0; trkphfrac = 0; trkphplane = 0; } reco_enu = reco_emu + reco_eshwcc; reco_qe_enu = RecoQENuEnergy(track_index); if (reco_enu>0) { reco_y = reco_eshwcc/reco_enu; } reco_qe_q2 = RecoQEQ2(track_index); reco_dircosz = RecoMuDCosZ(track_index); reco_dircosneu = RecoMuDCosNeu(track_index); is_cev = IsFidAll(track_index); // only calculate likelihood PID for events that pass event // fiducial cuts and track quality cuts if(is_fid_dp){ pid0 = PID(event,0); pid1 = PID(event,1); pid2 = PID(event,2); if(PassAnalysisCuts(event)) pass = 1; } else{ pid0 = -999.; pid1 = -999.; pid2 = -999.; pass = 0; } tree->Fill(); } // end of EVENT LOOP } // end of SCANINFO loop } // end of NTUPLE ENTRIES loop // delete nuparent; delete spinfo; save.cd(); tree->Write(); save.Write(); save.Close(); }

Double_t MadDpAnalysis::GetAnnPid (  AnnInputBlock  anninput, Int_t  det, Int_t  tar, Int_t  fid, Int_t  prior, Bool_t  first_ann )  [protected]

Definition at line 392 of file MadDpAnalysis.cxx. References Ann_weight_vector, AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aShwplu, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, det, and Sigmoid(). Referenced by CreatePAN(). { // DET IS THE DETECTOR 1 = NEAR 2 = FAR // TAR IS THE ENERGY 1 = LE 2 = PME 3 = PHE // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS 1 = MINE (MORE STRICT) // PRIOR IS THE FAR TRAINED METHOD , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 // 3 = DM2 0.002 SINTHETA2 = 0.95 // int inneuron = 13; int inneuron = 7; int hidneuron = 15; double var; double out[15]; // input neurons double rin[15]; // first hidden layer // double weight1[13][15];// weights first layer double weight1[7][15]; double constant1[15]; // constant first layer double weighto[15];// weights second (output) layer double constanto[1]; // constant second (output) layer double prob; // Initialize for(Int_t i=0;i<hidneuron;i++) { out[i]=0.; rin[i]=0.; constant1[i]=0.; } for(Int_t i=0;i<inneuron;i++) { for(Int_t j=0;j<hidneuron;j++){ weight1[i][j]=0.; } } constanto[0]=0.; for(Int_t i=0;i<hidneuron;i++){ weighto[i]=0.; } Int_t all = inneuron*hidneuron+hidneuron+hidneuron+1; Int_t all1 = inneuron*hidneuron+hidneuron; Int_t n1 =0; Int_t nw1 =0; Int_t no =0; Int_t nwo =0; ifstream weightfile; // INPUT VARIABLES if(first_ann){ if(det==2){ if(prior==1){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); } if(prior==2){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); } if(prior==3){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); } } if(det==1){ if(tar==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle_cedar_daikon7v2.dat"); if(tar==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme_cedar_daikon7v2.dat"); if(tar==3) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe_cedar_daikon7v2.dat"); } for(Int_t k=0;k<all;k++){ weightfile >> var ; Ann_weight_vector[k]=var; } weightfile.close(); } if(det==2){ out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.; out[1] = (anninput.aPhperpl/10000.); // out[1] = (anninput.aTotph/50000.); out[2] = (anninput.aTotlen/40); out[3] = (anninput.aPh3/(anninput.aTotph+1.)); out[4] = (anninput.aPh6/(anninput.aTotph+1.)); out[5] = (anninput.aPhlast/(anninput.aTotph+1.)); out[6] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.; /* out[0] = anninput.aTrkphperpl/2000.; out[1] = anninput.aTotrk*anninput.aTrkpass; out[2] = (anninput.aTotstp/100.); out[3] = (anninput.aTotph/50000.); out[4] = (anninput.aTotlen/40); out[5] = (anninput.aPhperpl/5000.); out[6] = (anninput.aPhperstp/1000.); out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.; out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.; */ } if(det==1){ out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.; out[1] = (anninput.aPhperpl/10000.); //out[1] = (anninput.aTotph/100000.); out[2] = (anninput.aTotlen/40); out[3] = (anninput.aPh3/(anninput.aTotph+1.)); out[4] = (anninput.aPh6/(anninput.aTotph+1.)); out[5] = (anninput.aPhlast/(anninput.aTotph+1.)); out[6] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.; /* out[0] = anninput.aTrkphperpl/5000.; out[1] = anninput.aTotrk*anninput.aTrkpass; out[2] = (anninput.aTotstp/100.); out[3] = (anninput.aTotph/100000.); out[4] = (anninput.aTotlen/40); out[5] = (anninput.aPhperpl/10000.); out[6] = (anninput.aPhperstp/2000.); out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.; out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.; */ } /* out[7] = (anninput.aPh3/(anninput.aTotph+1.)); out[8] = (anninput.aPh6/(anninput.aTotph+1.)); out[9] = (anninput.aPhlast/(anninput.aTotph+1.)); out[12] = (anninput.aShwphperdig/800.); */ // INPUT FILE for(Int_t k=0;k<all;k++){ if(k<all1){ if(k%(inneuron+1)==0){ nw1=0; constant1[n1]=Ann_weight_vector[k]; n1=n1+1; } else { weight1[nw1][n1-1]=Ann_weight_vector[k]; nw1=nw1+1; } } else if(k>=all1){ if((k-all1)%(hidneuron+1)==0){ nwo=0; constanto[no]=Ann_weight_vector[k]; no=no+1; } else { weighto[nwo]=Ann_weight_vector[k]; nwo=nwo+1; } } } // end of reading the weights // first layer for(Int_t i=0;i<hidneuron;i++){ rin[i]=constant1[i]; for(Int_t j=0;j<inneuron;j++){ rin[i]=rin[i]+weight1[j][i]*out[j]; } } // output for(Int_t i=0;i<hidneuron;i++){ out[i]=Sigmoid(rin[i]); } prob=constanto[0]; for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i]; if(anninput.aTotlen>50) prob=0.; return prob; }

bool MadDpAnalysis::InFidVol (  const Detector::Detector_t &  det, const Float_t &  x, const Float_t &  y, const Float_t &  z )  [static]

Definition at line 2017 of file MadDpAnalysis.cxx. References det. { // Mike Kordosky: 11 August 2005 // Defines a fiducial volume used to pass events which fill the PDFs // Transcribed from MakeMyFile() above bool is_fid=false; if(det==Detector::kFar){ is_fid = true; if(z<0.5 || z>29.4 || //ends (z<16.5&&z>14.5) || //between SMs sqrt((x*x) //radial cut +(y*y))>3.5) is_fid = false; } else if(det==Detector::kNear){ //fiducial criteria on vtx for near detector is_fid = true; if(z<1 || z>5 || sqrt(((x-1.4885)*(x-1.4885)) + ((y-0.1397)*(y-0.1397)))>1) is_fid = false; } return is_fid; }

void MadDpAnalysis::MakeMyFile (  std::string  , int  , int  )

Definition at line 657 of file MadDpAnalysis.cxx. References NtpSRPlane::beg, NtpSRPlane::end, MadQuantities::Flavour(), VldContext::GetDetector(), MadBase::GetEntry(), RecHeader::GetVldContext(), MadQuantities::IAction(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadTHTrack(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), NtpSRPlane::n, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSREvent::ntrack, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, MadQuantities::TrueNuEnergy(), NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. { // pdf binning: 0 - old binning, max 50 bins, 1 - new finer binning std::string savename = "DPHistos_" + tag + ".root"; TFile save(savename.c_str(),"RECREATE"); save.cd(); //#declare lots of histos: Int_t evlbins=0; if (pdfbinning==0) evlbins=60; else evlbins=300; TH1F *evlength_nc = new TH1F("Event length - nc", "Event length - nc", evlbins,0,600); evlength_nc->SetXTitle("Event length (planes)"); TH1F *evlength_cc = new TH1F("Event length - cc", "Event length - cc", evlbins,0,600); evlength_cc->SetXTitle("Event length (planes)"); TH1F *phfrac_nc = new TH1F("Track ph frac - nc", "Track ph frac - nc", 50,0,1); phfrac_nc->SetXTitle("pulse height fraction"); TH1F *phfrac_cc = new TH1F("Track ph frac - cc", "Track ph frac - cc", 50,0,1); phfrac_cc->SetXTitle("pulse height fraction"); Int_t phplanebins=0; Float_t phplanemax=0; if (pdfbinning==0) { phplanebins=50; phplanemax=5000; } else { phplanebins=100; phplanemax=3000; } TH1F *phplane_nc = new TH1F("ph per plane (nc)", "ph per plane (nc)", phplanebins,0,phplanemax); phplane_nc->SetXTitle("pulse height per plane"); TH1F *phplane_cc = new TH1F("ph per plane (cc)", "ph per plane (cc)", phplanebins,0,phplanemax); phplane_cc->SetXTitle("pulse height per plane"); evlength_nc->SetDirectory(&save); evlength_cc->SetDirectory(&save); phfrac_nc->SetDirectory(&save); phfrac_cc->SetDirectory(&save); phplane_nc->SetDirectory(&save); phplane_cc->SetDirectory(&save); // gDirectory->ls(); for(int i=0;i<Nentries;i++){ if(i%1000==0) std::cout << float(i)*100./float(Nentries) << "% done" << std::endl; if(!this->GetEntry(i)) continue; Int_t nevent = eventSummary->nevent; if(nevent==0) continue; Int_t trkcount=0; Int_t shwcount=0; Int_t evtmin=0; Int_t evtmax=nevent; if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) { Int_t maxevent=0; Float_t maxph=0; for(int j=0;j<nevent;j++){ if(!LoadEvent(j)) continue; Float_t evtpheight=ntpEvent->ph.sigcor; if (evtpheight>maxph) { maxph=evtpheight; maxevent=j; } } evtmin=maxevent; evtmax=maxevent+1; } // EVENT LOOP - fill tree once for each reconstructed event for(int j=evtmin;j<evtmax;j++){ if(!LoadEvent(j)) continue; //no event found Int_t ntrack = 0; ntrack=ntpEvent->ntrack; Int_t nshower = 0; nshower=ntpEvent->nshower; Float_t totph=0; for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) { LoadTrack(ii); totph+=ntpTrack->ph.sigcor; LoadTHTrack(ii); } trkcount+=ntrack; for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) { LoadShower(ii); totph+=ntpShower->ph.sigcor; } shwcount+=nshower; Int_t cc_nc = 0; Int_t flavour = 0; Int_t detector = -1; Int_t is_fid = 0; Double_t neu_e = -1; Double_t weight = -1; //get detector type: if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ detector = 2; //fiducial criteria on vtx for far detector if (evlength_cc->GetNbinsX()==evlbins) { if (pdfbinning==0) { evlength_cc->SetBins(50,0,500); evlength_nc->SetBins(50,0,500); } else { evlength_cc->SetBins(250,0,500); evlength_nc->SetBins(250,0,500); } } is_fid = 1; if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 || //ends (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) || //between SMs sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x) //radial cut +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0; } else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ detector = 1; if (pdfbinning==1 && evlength_cc->GetBinLowEdge(250)>200) { evlength_cc->SetBins(150,0,300); evlength_nc->SetBins(150,0,300); } if (pdfbinning==0 && evlength_cc->GetBinLowEdge(50)>200) { evlength_cc->SetBins(60,0,300); evlength_nc->SetBins(60,0,300); } //fiducial criteria on vtx for near detector is_fid = 1; if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 || sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) + ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) is_fid = 0; } //check for a corresponding mc event Int_t mcevent=-1; if(LoadTruthForRecoTH(j,mcevent)) { //true info for tree: cc_nc = IAction(mcevent); flavour = Flavour(mcevent); neu_e = TrueNuEnergy(mcevent); // tarpos 1 = LE 2 = ME 3 = HE if(tarpos==2) weight=1; if(tarpos==1) weight=1; if(tarpos==3) weight=1; } int track_index = -1; LoadLargestTrackFromEvent(j,track_index); if (track_index!=-1) { // if (is_fid && trkpass) { if (is_fid) { Float_t evlength=0; Float_t phfrac=0; Float_t phplane=0; evlength=ntpEvent->plane.n; if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1; } Float_t phtrack=ntpTrack->ph.sigcor; Float_t phevent=ntpEvent->ph.sigcor; if (phevent>0) {phfrac=phtrack/phevent;} Float_t trkplane=ntpTrack->plane.n; if (trkplane>0) {phplane=phtrack/trkplane;} if(flavour==2 && cc_nc==1) { evlength_cc->Fill(evlength,weight); phfrac_cc->Fill(phfrac,weight); phplane_cc->Fill(phplane,weight); // cout<<"I've filled a histogram, where it is, I don't know"<<endl; // cout<<phplane_cc->GetEntries()<<endl; // gDirectory->ls(); } else if(cc_nc==0) { evlength_nc->Fill(evlength,weight); phfrac_nc->Fill(phfrac,weight); phplane_nc->Fill(phplane,weight); } // std::cout<<"evlenght "<<evlength<<" weight "<<weight<<std::endl; } } } } //#close file // std::cout<<"writing hist "<<endl; // save.ls(); save.Write(); save.Close(); }

Bool_t MadDpAnalysis::PassAnalysisCuts (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 88 of file MadDpAnalysis.cxx. References VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadEvent(), and PID(). Referenced by CreatePAN(). { if(!LoadEvent(event)) return false; // if(PassBasicCuts(event)) return true; Float_t pidcut=-0.4; if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;} if(!(PID(event,0)>pidcut)) return false; return true; }

Bool_t MadDpAnalysis::PassBasicCuts (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 152 of file MadDpAnalysis.cxx. References MadBase::LoadEvent(), and NtpSREvent::ntrack. Referenced by CreatePAN(). { if(!LoadEvent(event)) return false; if(ntpEvent->ntrack==0) return false; // Int_t track = -1; // LoadLargestTrackFromEvent(event,track); // if(!PassTrackCuts(track)) return false; // if(!IsFidVtxEvt(ntpEvent->index)) return false; // if(!IsFidVtx(track)) return false; return true; }

Bool_t MadDpAnalysis::PassFid (  Int_t  event = 0  )  [protected]

Definition at line 163 of file MadDpAnalysis.cxx. References NtpSREvent::index, MadQuantities::IsFidVtxEvt(), and MadBase::LoadEvent(). Referenced by CreatePAN(). { if(!LoadEvent(event)) return false; if(!IsFidVtxEvt(ntpEvent->index)) return false; return true; }

Bool_t MadDpAnalysis::PassFidNew (  Int_t  event = 0  )  [protected]

Definition at line 98 of file MadDpAnalysis.cxx. References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRTrack::vtx, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by CreatePAN(). { if(!LoadEvent(event)) return false; Int_t track = -1; LoadLargestTrackFromEvent(event,track); // for events with no track, use event vertex if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 || //ends (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) || //between SMs (pow(ntpEvent->vtx.x,2)+ //radial cut pow(ntpEvent->vtx.y,2))>14)) return false; // otherwise, use track vertex if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 || //ends (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) || //between SMs (pow(ntpTrack->vtx.x,2)+ //radial cut pow(ntpTrack->vtx.y,2))>14)) return false; return true; }

Bool_t MadDpAnalysis::PassFidNewN (  Int_t  event = 0  )  [protected]

Definition at line 126 of file MadDpAnalysis.cxx. References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRTrack::vtx, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by CreatePAN(). { if(!LoadEvent(event)) return false; Int_t track = -1; LoadLargestTrackFromEvent(event,track); // for events with no track, use event vertex if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 || //ends (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) || //between SMs (pow(ntpEvent->vtx.x,2)+ //radial cut pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false; // otherwise, use track vertex if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 || //ends (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) || //between SMs (pow(ntpTrack->vtx.x,2)+ //radial cut pow(ntpTrack->vtx.y,2))>(3.5*3.5))) return false; return true; }

Bool_t MadDpAnalysis::PassTruthSignal (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 80 of file MadDpAnalysis.cxx. References abs(), NtpMCTruth::iaction, NtpMCTruth::inu, and MadBase::LoadTruth(). { if(!LoadTruth(mcevent)) return false; if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true; return false; }

Float_t MadDpAnalysis::PID (  Int_t  event = 0, Int_t  method = 0 )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 579 of file MadDpAnalysis.cxx. References NtpSRPlane::beg, NtpSRPlane::end, fLikeliHist, VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRPlane::n, NtpSREvent::ph, NtpSRTrack::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, and NtpSRPulseHeight::sigcor. Referenced by CreatePAN(), and PassAnalysisCuts(). { if(!LoadEvent(event)) return -10; Int_t track = -1; if(!LoadLargestTrackFromEvent(event,track)) return -10; if(method!=0) return -10; Float_t ProbNC = 1.; Float_t ProbCC = 1.; Float_t PidCC = 0.; Float_t var1=eventSummary->plane.n; if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ var1=ntpEvent->plane.end-ntpEvent->plane.beg+1; } Float_t var2=0; Float_t var3=0; Float_t phtrack=ntpTrack->ph.sigcor; Float_t phevent=ntpEvent->ph.sigcor; /* In case of Birch // apply 1.8% correction to pulse height variables if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData ){ phtrack=phtrack*1.018; phevent=phevent*1.018; } */ if (phevent>0) {var2=phtrack/phevent;} if (ntpTrack->plane.n!=0) var3 = float(phtrack) /float(ntpTrack->plane.n); Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1)); Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2)); Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3)); if (prob1==0) {prob1=0.0001;} if (prob2==0) {prob2=0.0001;} if (prob3==0) {prob3=0.0001;} ProbCC=prob1*prob2*prob3; prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1)); prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2)); prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3)); if (prob1==0) {prob1=0.0001;} if (prob2==0) {prob2=0.0001;} if (prob3==0) {prob3=0.0001;} ProbNC=prob1*prob2*prob3; PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC))); return PidCC; }

void MadDpAnalysis::ReadPIDFile (  std::string   )

Definition at line 1993 of file MadDpAnalysis.cxx. References fLikeliFile, and fLikeliHist. { fLikeliFile = new TFile(tag.c_str(),"READ"); if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found fLikeliFile->cd(); fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc"); fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc"); fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc"); fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc"); fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)"); fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)"); for(int i=0;i<6;i++){ float integ = fLikeliHist[i]->Integral(); fLikeliHist[i]->Scale(1./integ); } } else fLikeliFile = NULL; }

Float_t MadDpAnalysis::RecoAnalysisEnergy (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 645 of file MadDpAnalysis.cxx. References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadQuantities::RecoMuEnergy(), and MadQuantities::RecoShwEnergy(). { if(!LoadEvent(event)) return 0; int largestTrack = -1; LoadLargestTrackFromEvent(event,largestTrack); float emu = RecoMuEnergy(0,largestTrack); float ehad = RecoShwEnergy(event); float ereco = emu + ehad; return ereco; }

Double_t MadDpAnalysis::Sigmoid (  Double_t  x  )  [protected]

Definition at line 376 of file MadDpAnalysis.cxx. Referenced by GetAnnPid(). { double sig; if(x>37.){ sig = 1.; } else if(x<-37.){ sig = 0.; } else { sig = 1./(1.+exp(-x)); } return sig; }

Float_t MadDpAnalysis::SingleTimeFrame (  Int_t  snarlentry, Int_t  nentries )  [protected]

Definition at line 919 of file MadDpAnalysis.cxx. References MadBase::GetEntry(), and RecPhysicsHeader::GetTimeFrame(). { Int_t tf=0,tfbef=0,tfaft=0; if(this->GetEntry(snarlentry)) tf = ntpHeader->GetTimeFrame(); if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1)) tfaft = ntpHeader->GetTimeFrame(); if(snarlentry>0 && this->GetEntry(snarlentry-1)) tfbef = ntpHeader->GetTimeFrame(); Float_t singletimeframe=0; if(snarlentry<(nentries-1) && snarlentry>0 ){ if(tf!=tfaft && tf!=tfbef) singletimeframe=1; } if(snarlentry==(nentries-1)) { if(tf!=tfbef) singletimeframe=1; } if(snarlentry==0) { if(tf!=tfaft) singletimeframe=1; } this->GetEntry(snarlentry); return singletimeframe; }

---

Member Data Documentation

Double_t MadDpAnalysis::Ann_weight_vector[226] [private]

Definition at line 33 of file MadDpAnalysis.h. Referenced by GetAnnPid().

TFile* MadDpAnalysis::fLikeliFile [protected]

Definition at line 28 of file MadDpAnalysis.h. Referenced by MadDpAnalysis(), ReadPIDFile(), and ~MadDpAnalysis().

TH1F* MadDpAnalysis::fLikeliHist[6] [protected]

Definition at line 29 of file MadDpAnalysis.h. Referenced by MadDpAnalysis(), PID(), and ReadPIDFile().

---

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

• MadDpAnalysis.h
• MadDpAnalysis.cxx

---