MadTestAnalysis Class Reference

#include <MadTestAnalysis.h>

Inheritance diagram for MadTestAnalysis:

List of all members.

Public Member Functions
	MadTestAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
	MadTestAnalysis (JobC *, string, int)
	~MadTestAnalysis ()
void	MakeMyFile (std::string, int)
void	ReadPIDFile (std::string)
void	CreatePAN (std::string, 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)
Public Attributes
MadNsID	nsid
MadDpID	dpid
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)
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)
Double_t	Sigmoid (Double_t x)
Float_t	SingleTimeFrame (Int_t snarlentry, Int_t nentries)
Protected Attributes
TFile *	fLikeliFile
TH1F *	fLikeliHist [6]

---

Constructor & Destructor Documentation

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

Definition at line 25 of file MadTestAnalysis.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; }

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

Definition at line 52 of file MadTestAnalysis.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; }

MadTestAnalysis::~MadTestAnalysis (   )

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

---

Member Function Documentation

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

Definition at line 118 of file MadTestAnalysis.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::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, MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadStrip(), MadBase::LoadTrack(), 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, and NtpSRStrip::strip. 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.; // // Calculate variables needed for ANN (and a few more) LoadEvent(event) ; int track_index = -1; LoadLargestTrackFromEvent(event,track_index); int shower_index = -1; LoadLargestShowerFromEvent(event,shower_index); LoadTrack(track_index); LoadShower(shower_index); 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 anninput.aPhperpl =ntpEvent->ph.sigcor/ntpEvent->plane.n; 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.n; anninput.aShwphper =ntpShower->ph.sigcor/ntpEvent->ph.sigcor; anninput.aShwphperpl =ntpShower->ph.sigcor/ntpShower->plane.n; anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit; 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 for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); LoadStrip(stp_index); 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]); 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]); 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 of looping over event strips.... // anninput.aPh3 =ph3; anninput.aPh6 =ph6; anninput.aPhlast =phlast; anninput.aPhcommon =phcommon; return anninput; }

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

Definition at line 724 of file MadTestAnalysis.cxx. References AnnVar(), NtpSRPlane::beg, MadDpID::CalcPID(), NtpSREventSummary::date, det, dpid, NtpSRTrack::ds, NtpSRTrack::end, NtpSRPlane::end, NtpSRMomentum::eqp, BeamMonSpill::fHornCur, MadQuantities::Flavour(), 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(), GnumiInterface::GetParent(), MadNsID::GetPID(), 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(), MadBase::LoadEvent(), MadBase::LoadLargestShowerFromEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), BDSpillAccessor::LoadSpill(), MadBase::LoadStrip(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), NtpSRTrack::momentum, NtpSRDate::month, month, NtpSRPlane::n, NtpSREventSummary::nevent, NtpSREvent::nshower, nsid, NtpSREvent::nstrip, NtpSREvent::ntrack, MadQuantities::Nucleus(), PassAnalysisCuts(), PassBasicCuts(), PassFid(), NtpSRShower::ph, NtpSRTrack::ph, NtpSREventSummary::ph, NtpSREvent::ph, PID(), NtpSRTrack::plane, NtpSREvent::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, SingleTimeFrame(), BeamMonSpill::SpillTime(), GnumiInterface::Status(), NtpSREvent::stp, striptime, MadQuantities::Target4Vector(), NtpSRStrip::time0, NtpSRStrip::time1, NtpSREventSummary::trigtime, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSRVertex::u, NtpSRVertex::v, NtpSRTrack::vtx, NtpSREvent::vtx, MadQuantities::W2(), NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::Y(), NtpSRDate::year, NtpSRVertex::z, NuParent::Zero(), and SpillInfoBlock::Zero(). { 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; //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_cev; // fully contained. true = 1, false = 0 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 Int_t pass; // pass analysis cuts. true = 1, false = 0 Float_t npid0; // pid parameter (usual method) Float_t npid1; // pid parameter (alternative method 1) Float_t npid2; // pid parameter (alternative method 2) Double_t nannpid; // ANN PID 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) 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 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 // 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; // WEIGHT FOR THE PME BEAM Double_t w_le_me; // 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; // IS SNARL CUT IN PIECES Float_t singletimeframe; // Beam Info Block SpillInfoBlock *spinfo = new SpillInfoBlock(); spinfo->Zero(); singletimeframe=-999; //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("is_fid",&is_fid,"is_fid/I",512000); tree->Branch("is_cev",&is_cev,"is_cev/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("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_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("trkphfrac",&trkphfrac,"trkphfrac/F",512000); tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000); tree->Branch("rawph",&rawph,"rawph/F",512000); // <- !! added this tree->Branch("w_le_me",&w_le_me,"w_le_me/D",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("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"); // 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"); // new pid variables tree->Branch("npid0",&npid0,"npid0/F",512000); tree->Branch("npid1",&npid1,"npid1/F",512000); tree->Branch("npid2",&npid2,"npid2/F",512000); tree->Branch("nannpid",&nannpid,"nannpid/D",512000); // MAIN LOOP (over tree entries) for(int i=0;i<Nentries;i++){ if(i%100==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(); Double_t snarltime=ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9); Int_t month =eventSummary->date.month; Int_t year =eventSummary->date.year; // 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(); } // 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; // EVENT LOOP - fill tree once for each reconstructed event for(int j=0;j<nevent;j++){ if(!LoadEvent(j)) continue; //no event found totph = 0; //set event index: event = j; ntrack = ntpEvent->ntrack; nshower = ntpEvent->nshower; // Initialize everything nuparent->Zero(); 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.; w_le_me=-1; 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.; npid0 = -999; npid1 = -999; npid2 = -999; nannpid=-999; 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; } //check for a corresponding mc event 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); } 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 // ANN PID // AnnInputBlock anninput=AnnVar(event); // AnnInputBlock anninput2=anninput; // anninput2.aTotrk=1000; // if(!nsid.CompareAnnBlocks(this,event,anninput2)) assert(false); Int_t target=0; if(ntpHeader->GetVldContext().GetSimFlag()!=4){ if(tar<900) target=1; if(tar>900 && tar<40000.) target=2; if(tar>40000 && tar <100000) target=3; } if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){ Int_t mod=(int)(run/1e6); if(mod==14) target=1; if(mod==16) target=2; if(mod==18) 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 =2; Int_t prior = 1; if(!PassFid(event)) { annpid=-999; nannpid=-999; } else { Detector::Detector_t det= ntpHeader->GetVldContext().GetDetector(); // ANN PID // AnnInputBlock anninput=AnnVar(event); if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ annpid=GetAnnPid(anninput,detector,target,fid,prior); if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999; if(!((reco_vtxz>1. &&reco_vtxz<5.)&&sqrt((reco_vtxx-1.4885)*(reco_vtxx-1.4885)+(reco_vtxy-0.1397)*(reco_vtxy-0.1397))<1)) { annpid=-999.; nannpid=-999; } } if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ annpid=GetAnnPid(anninput,detector,target,fid,prior); if(!nsid.GetPID(this,event,det,nannpid)) nannpid=-999; if(fid==1) if(!(((reco_vtxz>1. &&reco_vtxz<14.)||(reco_vtxz>17. &&reco_vtxz<27.))&& sqrt(reco_vtxx*reco_vtxx+reco_vtxy*reco_vtxy)<3.5)) { annpid=-999.; nannpid=-999; } } } Double_t timemin=9.*10e30; Double_t timemax=-9999999; Double_t trgtime=eventSummary->trigtime; // Find max min time of events by loading strips for ND if(detector==1){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); LoadStrip(stp_index); Double_t striptime=ntpStrip->time1; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } evttimemax=timemax-trgtime; evttimemin=timemin-trgtime; } if(detector==2){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); LoadStrip(stp_index); 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; } evttimemax=timemax-trgtime; evttimemin=timemin-trgtime; } //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); //methods should all be protected against index = -1 reco_emu = RecoMuEnergy(0,track_index); reco_eshw = RecoShwEnergy(shower_index); reco_eshw_sqrt = RecoShwEnergySqrt(event); reco_enu = reco_emu + reco_eshw; reco_qe_enu = RecoQENuEnergy(track_index); if (reco_enu>0) { reco_y = reco_eshw/reco_enu; } reco_qe_q2 = RecoQEQ2(track_index); reco_dircosz = RecoMuDCosZ(track_index); reco_dircosneu = RecoMuDCosNeu(track_index); is_cev = IsFidAll(track_index); if(track_index!=-1){ //check track is present before using ntpTrack 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; Float_t phtrack=ntpTrack->ph.sigcor; Float_t phevent=ntpEvent->ph.sigcor; if (phevent>0) {trkphfrac=phtrack/phevent;} if(ntpTrack->plane.n>0) { trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n; } } else { trklength = 0; trkmomrange = 0; trkqp = 0; trkeqp = 0; trkphfrac = 0; trkphplane = 0; } // only calculate likelihood PID for events that pass event // fiducial cuts and track quality cuts if(PassBasicCuts(event) && PassFid(event)){ pid0 = PID(event,0); pid1 = PID(event,1); pid2 = PID(event,2); npid0 = dpid.CalcPID(this,event,0); npid1 = dpid.CalcPID(this,event,1); npid2 = dpid.CalcPID(this,event,2); if(PassAnalysisCuts(event)) pass = 1; } else{ pid0 = -999.; pid1 = -999.; pid2 = -999.; pass = 0; npid0 = -999.; npid1 = -999.; npid2 = -999.; } 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 MadTestAnalysis::GetAnnPid (  AnnInputBlock  anninput, Int_t  det, Int_t  tar, Int_t  fid, Int_t  prior )  [protected]

Definition at line 297 of file MadTestAnalysis.cxx. References AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, 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 hidneuron = 15; double var; double out[15]; // input neurons double rin[15]; // first hidden layer double weight1[13][15];// weights first layer 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(det==2){ if(prior==1){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnew12.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farnooscilnewfid12.dat"); } if(prior==2){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1new12.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc1newfid12.dat"); } if(prior==3){ if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2new12.dat"); if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farosc2newfid12.dat"); } 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.; } else if(det==1){ if(tar==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle16.dat"); if(tar==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme16.dat"); if(tar==3) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe16.dat"); 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++){ weightfile >> var ; if(k<all1){ if(k%(inneuron+1)==0){ nw1=0; constant1[n1]=var; n1=n1+1; } else { weight1[nw1][n1-1]=var; nw1=nw1+1; } } else if(k>=all1){ if((k-all1)%(hidneuron+1)==0){ nwo=0; constanto[no]=var; no=no+1; } else { weighto[nwo]=var; nwo=nwo+1; } } } weightfile.close(); // 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>40) prob=0.; return prob; }

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

Definition at line 1422 of file MadTestAnalysis.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){ //fiducial criteria on vtx for far detector /* if (evlength_cc->GetNbinsX()==60) { evlength_cc->SetBins(50,0,500); evlength_nc->SetBins(50,0,500); } */ 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){ /* if (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 = 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 MadTestAnalysis::MakeMyFile (  std::string  , int  )

Definition at line 510 of file MadTestAnalysis.cxx. References NtpSRPlane::beg, NtpSRPlane::end, NtpSRTrack::fit, 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, NtpSRFitTrack::pass, NtpSREvent::ph, NtpSRShower::ph, NtpSRTrack::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, MadQuantities::TrueNuEnergy(), NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. { std::string savename = "DPHistos_" + tag + ".root"; TFile save(savename.c_str(),"RECREATE"); save.cd(); //#declare lots of histos: TH1F *evlength_nc = new TH1F("Event length - nc", "Event length - nc", 60,0,600); evlength_nc->SetXTitle("Event length (planes)"); TH1F *evlength_cc = new TH1F("Event length - cc", "Event length - cc", 60,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"); TH1F *phplane_nc = new TH1F("ph per plane (nc)", "ph per plane (nc)", 50,0,5000); phplane_nc->SetXTitle("pulse height per plane"); TH1F *phplane_cc = new TH1F("ph per plane (cc)", "ph per plane (cc)", 50,0,5000); phplane_cc->SetXTitle("pulse height per plane"); 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; //fill tree once for each reconstructed event for(int j=0;j<eventSummary->nevent;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()==60) { evlength_cc->SetBins(50,0,500); evlength_nc->SetBins(50,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 (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) { Int_t trkpass=ntpTrack->fit.pass; if (is_fid && trkpass) { 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); } else if(cc_nc==0) { evlength_nc->Fill(evlength,weight); phfrac_nc->Fill(phfrac,weight); phplane_nc->Fill(phplane,weight); } } } } } //#close file save.Write(); save.Close(); }

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

Reimplemented from MadAnalysis. Definition at line 86 of file MadTestAnalysis.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 MadTestAnalysis::PassBasicCuts (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 98 of file MadTestAnalysis.cxx. References NtpSREvent::index, MadQuantities::IsFidVtxEvt(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSREvent::ntrack, and MadQuantities::PassTrackCuts(). 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 MadTestAnalysis::PassFid (  Int_t  event = 0  )  [protected]

Definition at line 109 of file MadTestAnalysis.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 MadTestAnalysis::PassTruthSignal (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 78 of file MadTestAnalysis.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 MadTestAnalysis::PID (  Int_t  event = 0, Int_t  method = 0 )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 441 of file MadTestAnalysis.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; if (phevent>0) {var2=phtrack/phevent;} if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor) /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 MadTestAnalysis::ReadPIDFile (  std::string   )

Definition at line 1398 of file MadTestAnalysis.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 MadTestAnalysis::RecoAnalysisEnergy (  Int_t  event = 0  )  [protected, virtual]

Reimplemented from MadAnalysis. Definition at line 498 of file MadTestAnalysis.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 MadTestAnalysis::Sigmoid (  Double_t  x  )  [protected]

Definition at line 281 of file MadTestAnalysis.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 MadTestAnalysis::SingleTimeFrame (  Int_t  snarlentry, Int_t  nentries )  [protected]

Definition at line 699 of file MadTestAnalysis.cxx. References MadBase::GetEntry(), and RecPhysicsHeader::GetTimeFrame(). Referenced by CreatePAN(). { 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

MadDpID MadTestAnalysis::dpid

Definition at line 50 of file MadTestAnalysis.h. Referenced by CreatePAN().

TFile* MadTestAnalysis::fLikeliFile [protected]

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

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

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

MadNsID MadTestAnalysis::nsid

Definition at line 49 of file MadTestAnalysis.h. Referenced by CreatePAN().

---

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

• MadTestAnalysis.h
• MadTestAnalysis.cxx

---