MadTVAnalysis Class Reference

#include <MadTVAnalysis.h>

Inheritance diagram for MadTVAnalysis:

List of all members.

Public Member Functions
void	CreatePAN (std::string s, const char *bmonpath="")
	MadTVAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
	MadTVAnalysis (JobC *, string, int)
virtual	~MadTVAnalysis ()
Float_t	ClosestPointToLine (const Float_t &x1, const Float_t &y1, const Float_t &x2, const Float_t &y2, const Float_t &xpt, const Float_t &ypt, Float_t &x, Float_t &y)
void	ForceMCBeam (const BeamType::BeamType_t &beam)
const BeamType::BeamType_t &	GetForcedMCBeam ()
void	SetDataUseMCBeam (bool x)
void	FillMCHists (TFile *fout)
void	SetRustemNearPID (string fname)
void	SetRustemFarPID (string fname)
void	SetRustemConfig (string fname)
Static Public Member Functions
Bool_t	PittInFidVol (Float_t x, Float_t y, Float_t z, Float_t u, Float_t v)
Int_t	PittTrkContained (Float_t x, Float_t y, Float_t z, Float_t u, Float_t v)
Int_t	InPartialRegion (UShort_t plane, UShort_t strip)
Int_t	NDRadialFidVol (Float_t x, Float_t y, Float_t z)
int	FarTrkContained (Float_t x, Float_t y, Float_t z)
Public Attributes
MadNsID	nsid
MadDpID	dpid
MadQEID	qeid
MadAbID	abid
bool	openedabpidfile
Protected Member Functions
void	SigInOut (Float_t &sigfull, Float_t &sigpart, Float_t &trkfull, Float_t &trkpart)
void	SigInOut (Int_t trkidx, Float_t &sigfull, Float_t &sigpart, Float_t &trkfull, Float_t &trkpart)
virtual Float_t	RecoMuDCosNeuND (Int_t itr=0, Float_t *vtx=0)
virtual Float_t	RecoMuDCosNeuFD (Int_t itr=0, Float_t *vtx=0)
virtual bool	InFidVol (Int_t evidx)
virtual bool	InFidVol ()
bool	InFidVol (float vx, float vy, float vz)
void	SetBECFile (const char *f)
Registry &	GetBECRegistry ()
Int_t	NPlanesInObj (TObject *rec, TObject *obj, float phcut, float &sumph)
Int_t	NStripsInObj (TObject *rec, TObject *obj, float phcut, float &sumph)
void	GetEvtTimeChar (double &timemin, double &timemax, double &timeln)
void	GetEvtTimeCharPHC (double &timemin, double &timemax, double &timeln)
int	FindBatchNumber (double mintimediff)
float	MuonShowerEnergy (const NtpSREvent *evt, const NtpSRTrack *trk)
void	EarlyActivity (int evidx, float &earlywph, float &earlywphpw, float &earlywphsphere, float &ph1000, float &ph500, float &ph100, float &lateph1000, float &lateph500, float &lateph100, float &lateph50)
void	DupRecoStrips (int evidx, int &nduprs, float &dupphrs)
float	ComputeLowPHRatio ()
Int_t	FDRCBoundary ()
Float_t	GetNDCoilCurrent (const VldContext &vc)
Private Attributes
Registry	fBECConfig
BeamType::BeamType_t	fMCBeam
bool	fDataUseMCBeam
Anp::Interface *	phnti
string	rustempidfile_near
string	rustempidfile_far
string	rustemconfigfile
Static Private Attributes
const Float_t	kVetoEndZ = 1.2154
const Float_t	kTargEndZ = 3.5914
const Float_t	kCalorEndZ = 7.1554
const Float_t	kHalfCalorEndZ = kTargEndZ+(kCalorEndZ-kTargEndZ)*0.5
const Float_t	kSpecEndZ = 16.656
const Float_t	kXcenter = 1.4828
const Float_t	kYcenter = 0.2384
const Double_t	fEarlyActivityWindowSize = 10.
const Double_t	fPMTTimeConstant = 700.
const Int_t	fNPlaneEarlyAct = 30
const Double_t	fEarlyPlaneSphere = 0.5

---

Constructor & Destructor Documentation

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

Definition at line 60 of file MadTVAnalysis.cxx. References MadBase::Clear(), fBECConfig, fDataUseMCBeam, fMCBeam, openedabpidfile, phnti, and Registry::Set(). : MadAnalysis(chainSR, chainMC, chainTH, chainEM) { openedabpidfile=false; fMCBeam=BeamType::kUnknown; fDataUseMCBeam=false; if(!chainSR) { record = 0; emrecord = 0; mcrecord = 0; threcord = 0; Clear(); whichSource = -1; isMC = true; isTH = true; isEM = true; Nentries = -1; return; } static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root"; fBECConfig.Set("beam:flux_file",default_bec_file); fBECConfig.Set("beam:beam_type",BeamType::kLE); fBECConfig.Set("beam:detector",Detector::kNear); fBECConfig.Set("beam:low_energy_limit",0.5); fBECConfig.Set("beam:high_energy_limit",60.0); // InitChain(chainSR,chainMC,chainTH,chainEM); whichSource = 0; phnti = new Anp::Interface(); }

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

Definition at line 96 of file MadTVAnalysis.cxx. References MadBase::Clear(), fBECConfig, fDataUseMCBeam, fMCBeam, openedabpidfile, phnti, and Registry::Set(). : MadAnalysis(j, path, entries) { openedabpidfile=false; fMCBeam=BeamType::kUnknown; fDataUseMCBeam=false; record = 0; emrecord = 0; mcrecord = 0; threcord = 0; Clear(); isMC = true; isTH = true; isEM = true; Nentries = entries; whichSource = 1; jcPath = path; fJC = j; fChain = NULL; static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root"; fBECConfig.Set("beam:flux_file",default_bec_file); fBECConfig.Set("beam:beam_type",BeamType::kLE); fBECConfig.Set("beam:detector",Detector::kNear); fBECConfig.Set("beam:low_energy_limit",0.5); fBECConfig.Set("beam:high_energy_limit",60.0); phnti = new Anp::Interface(); }

MadTVAnalysis::~MadTVAnalysis (   )  [virtual]

Definition at line 129 of file MadTVAnalysis.cxx. References phnti. { if(phnti){ delete phnti; phnti=0; } }

---

Member Function Documentation

Float_t MadTVAnalysis::ClosestPointToLine (  const Float_t &  x1, const Float_t &  y1, const Float_t &  x2, const Float_t &  y2, const Float_t &  xpt, const Float_t &  ypt, Float_t &  x, Float_t &  y )

Definition at line 2227 of file MadTVAnalysis.cxx. { // given a line with end points (x1,y1) (x2,y2) // and a spatial point (xpt,ypt) // report the closest point on the line (x,y) to the spatial point // and compute the distance from (xpt,ypt) to (x,y) Float_t u= (xpt-x1)*(x2-x1)+(ypt-y1)*(y2-y1); u/=sqrt( pow(x2-x1,2) + pow(y2-y1,2)); x=x1+u*(x2-x1); y=y1+u*(y2-y1); Float_t dist =sqrt(pow(xpt-x, 2)+pow(ypt-y,2)); return dist; }

float MadTVAnalysis::ComputeLowPHRatio (   )  [protected]

Definition at line 3134 of file MadTVAnalysis.cxx. References MadBase::GetEventSummary(), MadBase::LoadStrip(), NtpSREventSummary::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, and NtpSRStrip::ph1. Referenced by CreatePAN(). { float lowphsum=0.; float totphsum=0.; for(unsigned int i=0;i<GetEventSummary()->nstrip;i++){ if(!LoadStrip(i)) continue; totphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe; if(ntpStrip->ph1.pe+ntpStrip->ph0.pe<2){ lowphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe; } } float result=-1.; if(totphsum!=0){ result = lowphsum/totphsum; } return result; }

void MadTVAnalysis::CreatePAN (  std::string  s, const char *  bmonpath = "" )

Definition at line 137 of file MadTVAnalysis.cxx. References abid, abs(), Moments::AddPoint(), MadQEID::AltCalcQEID(), ReleaseType::AsString(), base, BeamMonSpill::BeamType(), NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, MadAbID::CalcPID(), MadDpID::CalcPID(), NtpSRFitTrack::chi2, MadDpID::ChoosePDFs(), MadNsID::ChooseWeightFile(), NearbyVertexFinder::ClosestSpatial(), NearbyVertexFinder::ClosestTemporal(), NtpSRShower::clu, CoilTools::CoilCurrent(), NtpSRDetStatus::coilstatus, NtpTHTrack::completeslc, ComputeLowPHRatio(), Anp::Interface::Config(), NtpSRTrack::contained, NtpSRVertex::dcosx, NtpSRVertex::dcosy, det, dpid, DupRecoStrips(), EarlyActivity(), NtpSRTrack::end, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, fBECConfig, BeamMonSpill::fBpmInt, fDataUseMCBeam, FDRCBoundary(), BeamMonSpill::fHadInt, BeamMonSpill::fHornCur, FillMCHists(), Anp::Interface::FillSnarl(), FindBatchNumber(), NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, NtpMCFluxInfo::fluxevtno, NtpMCFluxInfo::fluxrun, BeamMonSpill::fMuInt1, BeamMonSpill::fMuInt2, BeamMonSpill::fMuInt3, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTargProfX, BeamMonSpill::fTargProfY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTr101d, BeamMonSpill::fTrtgtd, EnergyCorrections::FullyCorrectMomentumFromRange(), EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(), BDSpillAccessor::Get(), Registry::Get(), VldContext::GetDetector(), MadBase::GetEntry(), GetEvtTimeChar(), GetEvtTimeCharPHC(), MadNsID::GetPID(), MadAbID::GetRecoCosTheta(), MadAbID::GetRecoE(), MadAbID::GetRecoY(), RecPhysicsHeader::GetRemoteSpillType(), Moments::GetRMS(), RecDataHeader::GetRun(), VldTimeStamp::GetSeconds(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), WeightCalculator::GetStandardConfig(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), MadAbID::GetTrackEMCharge(), MadAbID::GetTrackLikePlanes(), MadAbID::GetTrackPHfrac(), MadAbID::GetTrackPHmean(), MadAbID::GetTrackPlanes(), MadAbID::GetTrackQPsigmaQP(), RecPhysicsHeader::GetTrigSrc(), Anp::Interface::GetVar(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, gSystem(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSRCluster::id, InFidVol(), MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::INu(), MadQuantities::INuNoOsc(), MadQuantities::IResonance(), ReleaseType::IsCedar(), MadQuantities::IsFid_2008(), MadQuantities::IsFidAll(), MadQuantities::IsFidFD_AB(), MadQuantities::IsFidVtxEvt(), DataUtil::IsGoodFDData(), LISieve::IsLI(), CoilTools::IsOK(), CoilTools::IsReverse(), kXcenter, kYcenter, NtpSRShowerPulseHeight::linCCgev, NtpSREventSummary::litime, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), MadBase::LoadSlice(), BDSpillAccessor::LoadSpill(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), Registry::LockValues(), ReleaseType::MakeReleaseType(), NtpSRStripPulseHeight::mip, NtpSRTrack::momentum, MuonShowerEnergy(), NtpMCFluxInfo::mupare, NtpMCFluxInfo::muparpx, NtpMCFluxInfo::muparpy, NtpMCFluxInfo::muparpz, NtpSRPlane::n, NtpSRShower::ncluster, NtpMCFluxInfo::ndecay, NtpSRFitTrack::ndof, NDRadialFidVol(), NtpMCFluxInfo::ndxdz, NtpMCFluxInfo::ndxdzfar, NtpMCFluxInfo::ndxdznear, NtpMCFluxInfo::ndydz, NtpMCFluxInfo::ndydzfar, NtpMCFluxInfo::ndydznear, NtpMCFluxInfo::necm, NtpMCFluxInfo::nenergy, NtpMCFluxInfo::nenergyfar, NtpMCFluxInfo::nenergynear, NtpTHTrack::neumc, NtpSREventSummary::nevent, NtpMCFluxInfo::nimpwt, NtpSRDmxStatus::nonphysicalfail, NtpMCFluxInfo::norig, NPlanesInObj(), NtpMCFluxInfo::npz, NtpSREvent::nshower, nsid, NtpSRTrack::nstrip, NtpSREvent::nstrip, NStripsInObj(), NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpMCFluxInfo::ntype, MadQuantities::Nucleus(), MadQuantities::NumFSNeutron(), MadQuantities::NumFSPion(), MadQuantities::NumFSPiZero(), MadQuantities::NumFSProton(), NtpMCFluxInfo::nwtfar, NtpMCFluxInfo::nwtnear, openedabpidfile, NtpSRFitTrack::pass, MadAnalysis::PassBasicCuts(), MadQEID::PassMEDQECut(), MadQuantities::PassTrackCuts(), NtpMCFluxInfo::pdpx, NtpMCFluxInfo::pdpy, NtpMCFluxInfo::pdpz, NtpSRCluster::ph, NtpSRTrack::ph, NtpSREvent::ph, phnti, PittInFidVol(), PittTrkContained(), NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRCluster::planeview, NtpMCFluxInfo::ppdxdz, NtpMCFluxInfo::ppdydz, NtpMCFluxInfo::ppenergy, NtpMCFluxInfo::ppmedium, NtpMCFluxInfo::pppz, NtpMCFluxInfo::ppvx, NtpMCFluxInfo::ppvy, NtpMCFluxInfo::ppvz, Registry::Print(), NearbyVertexFinder::ProcessEntry(), NtpMCFluxInfo::ptype, MadQuantities::Q2(), qeid, NtpSRMomentum::qp, NtpSRMomentum::range, MadAbID::ReadPDFs(), RecoMuDCosNeuFD(), RecoMuDCosNeuND(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergyNew(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergyNew(), run(), rustemconfigfile, rustempidfile_far, rustempidfile_near, BMSpillAna::SelectSpill(), Registry::Set(), MadDpID::SetPHCorrection(), BMSpillAna::SetSpill(), BMSpillAna::SetTimeDiff(), NtpSRShower::shwph, NtpSRPulseHeight::sigcor, SigInOut(), NtpSREvent::slc, NtpSRCluster::slope, BeamMonSpill::SpillTime(), MadQuantities::Target4Vector(), NtpMCFluxInfo::tgen, MadQuantities::TotFSNeutronEnergy(), MadQuantities::TotFSPionEnergy(), MadQuantities::TotFSPiZeroEnergy(), MadQuantities::TotFSProtonEnergy(), NtpSRCluster::tposvtx, NtpMCFluxInfo::tptype, NtpMCFluxInfo::tpx, NtpMCFluxInfo::tpy, NtpMCFluxInfo::tpz, NtpSREventSummary::trigtime, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpMCFluxInfo::tvx, NtpMCFluxInfo::tvy, NtpMCFluxInfo::tvz, NtpSRVertex::u, BMSpillAna::UseDatabaseCuts(), NtpSRVertex::v, NtpSRDmxStatus::validplanesfail, NtpSRDmxStatus::vertexplanefail, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, NtpMCFluxInfo::vx, NtpMCFluxInfo::vy, NtpMCFluxInfo::vz, MadQuantities::W2(), NtpSRShowerPulseHeight::wtCCgev, NtpSRVertex::x, MadQuantities::X(), NtpMCFluxInfo::xpoint, NtpSRVertex::y, MadQuantities::Y(), NtpMCFluxInfo::ypoint, NtpSRVertex::z, and NtpMCFluxInfo::zpoint. { // is this MC?? this->GetEntry(0); const bool file_is_mc =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC); //set mc and recoversions for pid histogram retrevial // string recover = EnergyCorrections::GetCorrectionAsString(); string mcver = "daikon"; //note this is hard coded and bad, but I dont' care. ReleaseType::Release_t reltype = ReleaseType::MakeReleaseType(strecord->GetTitle(),mcver); string recover = ReleaseType::AsString(reltype); cout<<"Running pans over "<<recover<<" "<<mcver<<endl; //updated ireg to from rustems for his cc-pid selection 12-21-2007 //set up anti numu selection and rustem's cc pid Registry ireg(false); ireg.Set("InterfaceConfigPath", rustemconfigfile.c_str()); if(ntpHeader->GetVldContext().GetDetector() == Detector::kNear) { ireg.Set("FillkNNFilePath", rustempidfile_near.c_str()); } else if(ntpHeader->GetVldContext().GetDetector() == Detector::kFar) { ireg.Set("FillkNNFilePath", rustempidfile_far.c_str()); } // // Do not erase events that fail SelectAntiNeutrino algorithm // ireg.Set("QuietInterface", "yes"); ireg.Set("SelectAntiNeutrinoErase", "no"); ireg.LockValues(); phnti->Config(ireg); //PAN Quantities //Truth: Float_t true_enu; // true neutrino energy (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 Int_t true_pitt_fid; // was event truly in fiducial volue? Float_t true_trk_cmplt; //true track completeness //For Neugen: Float_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 inu; // stdhep inu Int_t inunoosc; // stdhep inu unoscilllated Int_t resnum; // IdHEP%1000 of resonance state Short_t npp; // # final state pi+ Short_t npm; // # fs pi- Short_t npz; // # fs pi0 Short_t np; // # fs p Short_t nn; // # fs n float epp; // energy final state pi+ float epm; // energy fs pi- float epz; // energy fs pi0 float ep; // energy fs p float en; // energy fs n //Reco Quantities Int_t detector; // Near = 1, Far = 2; Int_t run; // run number Int_t snarl; // snarl number Int_t subrun; // subrun number UInt_t trgsrc; // trigger source double trgtime; // trigger time Double_t vld_sec; // seconds field of header vld context Int_t spilltype; // comes from mdSpillTypeEnum // 0=none, 1=reported, 2=predicted // 3=fake, 4=local Int_t nevent; // number of events in snarl Int_t event; // event index Int_t mcevent; // mc event index Int_t ntrack; // number of tracks in event Int_t nshower; // number of showers in event Int_t nstrip; // number of strips in event Int_t is_fid; // pass fid cut. true = 1, false = 0 Int_t is_fid_2007; // pass fid cut. true = 1, false = 0 Int_t is_cev; // fully contained. true = 1, false = 0 Int_t sr_contained; // containment support from reconstruction. 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 Int_t is_pitt_fid; // passes pitt ND vtx fid cut true=1, false=0 Int_t is_trk_fid; // is the trkvtx in fid (pitt fid for ND) // true=1, false=0, notrack = -1 Int_t nd_radial_fid; // a bitfield, bits correspond to r cuts Int_t pitt_evt_class; // CC event class as defined by pitt group // 1 = track is fully contained in the upstream region // 2 = track is partially contained in the upstream region // 3 = track is fully contained in the downstream region // 4 = track is partially contained in the downstream region Int_t pitt_trk_qual; // pitt tracking quality cuts Int_t duvvtx; //delta (Uvtx-Vvtx) Float_t pid0; // pid parameter (usual method) Float_t pid1; // pid parameter (alternative method 1) Float_t pid2; // pid parameter (alternative method 2) Int_t pass_track; // the primary track passes track cuts Int_t trk_fit_pass; // trk.fit.pass Int_t pass; // pass analysis cuts. true = 1, false = 0 Int_t emu_meth; // method used to determine muon energy Float_t reco_enu; // reco neutrino energy Float_t reco_emu; // reco muon energy Float_t reco_mushw; // showers along muon track Float_t reco_eshw; // reco shower energy (largest shower Float_t reco_wtd_eshw; // as above but weighted Float_t reco_vtx_eshw; // reco shower energy (=0 if no vertex shower) Float_t reco_vtx_wtd_eshw;// as above but weighted (=0 if no vertex shower) Float_t reco_eshw_uncorr; // uncorrected shower energy Float_t reco_wtd_eshw_uncorr; // uncorrected weighted shower energy Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino Float_t reco_dircosx; // reco track vtx x-dircos Float_t reco_dircosy; // reco track vtx y-dircos 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 reco_vtxr;// radial position of vertex Float_t reco_trk_vtxx; // reco x track vtx Float_t reco_trk_vtxy; // reco y track vtx Float_t reco_trk_vtxz; // reco z track vtx Float_t reco_trk_vtxr;// radial position of vertex Float_t reco_shw_vtxx; // reco x track vtx Float_t reco_shw_vtxy; // reco y track vtx Float_t reco_shw_vtxz; // reco z track vtx Float_t reco_shw_vtxr;// radial position of vertex Float_t reco_trk_endx; // reco x track end Float_t reco_trk_endy; // reco y track end Float_t reco_trk_endz; // reco z track end Float_t reco_trk_endr;// radial position of vertex Int_t evt_nplanes; //number of planes in event above some ph Int_t trk_nplanes; //number of planes in track above some ph Int_t slc_nplanes; //number of planes in slice above some ph Int_t shw_nplanes_cut; //number of planes in shower above some ph Int_t shw_nplanes_raw; //number of planes in shower Float_t shw_ph_cut, shw_ph_raw; // sigcor in shower, w/ & w/o ph cut Int_t evt_nstrips; //number of strips in event above some ph Int_t trk_nstrips; //number of strips in track above some ph Int_t slc_nstrips; //number of strips in slice above some ph // reconstructed kinematics // formulae given for high-energy CC interactions // Assume nu = Ehad Float_t reco_y;// y = (Enu-Emu)/Enu Float_t reco_Q2;// Q2 = -q^{2} = 2 Enu Emu (1 - cos(theta_mu)) Float_t reco_x;// x = Q2 / (2M* Ehad) {M=nucleon mass} Float_t reco_W2;// W2 = M^{2} + q^{2} + 2M(Enu-Emu) // reco quantities, assuming a QE event Float_t reco_qe_enu; // reco qe neutrino energy Float_t reco_qe_q2; // reco qe q2 Float_t evtlength; // reco event length Float_t shwlength; // reco shower length Float_t trklength; // reco track length (planes) Int_t ntrklike; // number of track-like planes 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 trkchi2; Int_t trkndf; Int_t trkend, trkendu, trkendv; // track end plane, u, v Int_t shwend,shwbeg; // shower begin,end planes Int_t nss[6];// Float_t ess[6];// Int_t ntotss; // variables to tag nearest event in space and time // nvsi = nearest vertex space : index // nvsd = nearest vertex space : distance // nvst = nearest vertex space : distance // nvsp = nearest vertex space : pulse height // nvsevt = nearest vertex space : mc event // nvti = nearest vertex time : index // .... // bvsi = back vertex space : index // the "back" variables hold the same info as the "near" variables // but for that nearest event // I'm interested in this case // // (this vtx)------>(nearest)--->(nearest to it) // // if you are going to try to recombine events // you need to somehow account for this // Int_t nvsi, bvsi, nvti, bvti; Int_t nvsevt, bvsevt, nvtevt, bvtevt; Float_t nvsd, bvsd, nvtd, bvtd; Float_t nvst, bvst, nvtt, bvtt; Float_t nvsp, bvsp, nvtp, bvtp; // this event's pulseheight (in sigcor) Float_t evtph; Float_t evtphgev; // the pulse height of the event and track // in the fully and partially instrumented regions Float_t evtsigfull,trksigfull,evtsigpart,trksigpart; //time structure of event double evttimemin, evttimemax, evttimeln; double evttimeminphc, evttimemaxphc, evttimelnphc; //early activity variables float earlywph, earlywphpw, earlywphsphere; float ph1000, ph500, ph100; float lateph1000, lateph500, lateph100, lateph50; //duplicate reconstructed strip variables int nduprs; //number of reco strips that are duplicated in another event float dupphrs; //pulseheight of reco strips //that are duplicated in another event //time last li event in snarl, -1 if none; double litime; int is_li_sieve=0; //demux failures (for fd) UChar_t dmx_stat; // std::cout<<"Second Print"<<std::endl; // david's FD data quality Int_t dp_fd_quality; //brian's cut on lowph ratio to remove junk events //in FD that pass the dmx_stat variable float lowphrat; //Trees // output file std::string savename = "PAN_" + tag + ".root"; TFile save(savename.c_str(),"RECREATE"); save.cd(); TTree *tree = new TTree("pan","pan"); //Truth tree->Branch("true_enu",&true_enu,"true_enu/F",32000); tree->Branch("true_emu",&true_emu,"true_emu/F",32000); tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000); tree->Branch("true_x",&true_x,"true_x/F",32000); tree->Branch("true_y",&true_y,"true_y/F",32000); tree->Branch("true_q2",&true_q2,"true_q2/F",32000); tree->Branch("true_w2",&true_w2,"true_w2/F",32000); Float_t nu_px, nu_py, nu_pz; Float_t tar_px, tar_py, tar_pz, tar_e; tree->Branch("nu_px", &nu_px, "nu_px/F"); tree->Branch("nu_py", &nu_py, "nu_py/F"); tree->Branch("nu_pz", &nu_pz, "nu_pz/F"); tree->Branch("tar_px", &tar_px, "tar_px/F"); tree->Branch("tar_py", &tar_py, "tar_py/F"); tree->Branch("tar_pz", &tar_pz, "tar_pz/F"); tree->Branch("tar_e", &tar_e, "tar_e/F"); tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000); tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000); tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000); tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000); tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000); tree->Branch("true_pitt_fid",&true_pitt_fid,"true_pitt_fid/I",32000); tree->Branch("true_trk_cmplt",&true_trk_cmplt,"true_trk_cmplt/F",32000); //Neugen tree->Branch("flavour",&flavour,"flavour/F",32000); tree->Branch("process",&process,"process/I",32000); tree->Branch("nucleus",&nucleus,"nucleus/I",32000); tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000); tree->Branch("inu",&inu,"inu/I",32000); tree->Branch("inunoosc",&inunoosc,"inunoosc/I",32000); tree->Branch("initial_state",&initial_state,"initial_state/I",32000); tree->Branch("resnum",&resnum,"resnum/I",32000); tree->Branch("npp",&npp,"npp/S",32000); tree->Branch("npm",&npm,"npm/S",32000); tree->Branch("npz",&npz,"npz/S",32000); tree->Branch("np",&np,"np/S",32000); tree->Branch("nn",&nn,"nn/S",32000); tree->Branch("epp",&epp,"epp/F",32000); tree->Branch("epm",&epm,"epm/F",32000); tree->Branch("epz",&epz,"epz/F",32000); tree->Branch("ep",&ep,"ep/F",32000); tree->Branch("en",&en,"en/F",32000); //Reco tree->Branch("detector",&detector,"detector/I",32000); tree->Branch("run",&run,"run/I",32000); tree->Branch("snarl",&snarl,"snarl/I",32000); tree->Branch("subrun",&subrun,"subrun/I",32000); tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000); tree->Branch("trgtime",&trgtime,"trgtime/D",32000); tree->Branch("vld_sec",&vld_sec,"vld_sec/D",32000); tree->Branch("spilltype",&spilltype,"spiltype/I",32000); tree->Branch("nevent",&nevent,"nevent/I",32000); tree->Branch("event",&event,"event/I",32000); tree->Branch("mcevent",&mcevent,"mcevent/I",32000); tree->Branch("ntrack",&ntrack,"ntrack/I",32000); tree->Branch("nshower",&nshower,"nshower/I",32000); tree->Branch("nstrip",&nstrip,"nstrip/I",32000); tree->Branch("is_fid",&is_fid,"is_fid/I",32000); tree->Branch("is_fid_2007",&is_fid_2007,"is_fid_2007/I",32000); tree->Branch("is_cev",&is_cev,"is_cev/I",32000); tree->Branch("sr_contained",&sr_contained,"sr_contained/I",32000); tree->Branch("is_mc",&is_mc,"is_mc/I",32000); tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000); tree->Branch("pass_track",&pass_track,"pass_track/I",32000); tree->Branch("trk_fit_pass",&trk_fit_pass,"trk_fit_pass/I",32000); tree->Branch("emu_meth",&emu_meth,"emu_meth/I",32000); tree->Branch("is_pitt_fid",&is_pitt_fid,"is_pitt_fid/I",32000); tree->Branch("is_trk_fid",&is_trk_fid,"is_trk_fid/I",32000); tree->Branch("nd_radial_fid",&nd_radial_fid,"nd_radial_fid/I",32000); tree->Branch("pitt_evt_class",&pitt_evt_class,"pitt_evt_class/I",32000); tree->Branch("pitt_trk_qual",&pitt_trk_qual,"pitt_trk_qual/I",32000); tree->Branch("duvvtx",&duvvtx,"duvvtx/I",32000); tree->Branch("pid0",&pid0,"pid0/F",32000); tree->Branch("pid1",&pid1,"pid1/F",32000); tree->Branch("pid2",&pid2,"pid2/F",32000); tree->Branch("pass",&pass,"pass/I",32000); float med_qe_pid; int med_qe_pass; tree->Branch("med_qe_pid",&med_qe_pid,"med_qe_pid/F",32000); tree->Branch("med_qe_pass",&med_qe_pass,"med_qe_pass/I",32000); double niki_cc_pid; tree->Branch("niki_cc_pid",&niki_cc_pid,"niki_cc_pid/D",32000); float dave_cc_pid; tree->Branch("dave_cc_pid",&dave_cc_pid,"dave_cc_pid/F",32000); tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000); tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000); tree->Branch("reco_mushw",&reco_mushw,"reco_mushw/F",32000); tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000); tree->Branch("reco_wtd_eshw",&reco_wtd_eshw,"reco_wtd_eshw/F",32000); tree->Branch("reco_vtx_eshw",&reco_vtx_eshw,"reco_vtx_eshw/F",32000); tree->Branch("reco_vtx_wtd_eshw",&reco_vtx_wtd_eshw,"reco_vtx_wtd_eshw/F",32000); tree->Branch("reco_eshw_uncorr",&reco_eshw_uncorr,"reco_eshw_uncorr/F",32000); tree->Branch("reco_wtd_eshw_uncorr",&reco_wtd_eshw_uncorr,"reco_wtd_eshw_uncorr/F",32000); tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000); tree->Branch("reco_dircosx",&reco_dircosx,"reco_dircosx/F",32000); tree->Branch("reco_dircosy",&reco_dircosy,"reco_dircosy/F",32000); tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000); tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000); tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000); tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000); tree->Branch("reco_vtxr",&reco_vtxr,"reco_vtxr/F",32000); tree->Branch("reco_trk_vtxx",&reco_trk_vtxx,"reco_trk_vtxx/F",32000); tree->Branch("reco_trk_vtxy",&reco_trk_vtxy,"reco_trk_vtxy/F",32000); tree->Branch("reco_trk_vtxz",&reco_trk_vtxz,"reco_trk_vtxz/F",32000); tree->Branch("reco_trk_vtxr",&reco_trk_vtxr,"reco_trk_vtxr/F",32000); tree->Branch("reco_shw_vtxx",&reco_shw_vtxx,"reco_shw_vtxx/F",32000); tree->Branch("reco_shw_vtxy",&reco_shw_vtxy,"reco_shw_vtxy/F",32000); tree->Branch("reco_shw_vtxz",&reco_shw_vtxz,"reco_shw_vtxz/F",32000); // tree->Branch("reco_shw_vtxr",&reco_shw_vtxr,"reco_shw_vtxr/F",32000); tree->Branch("reco_trk_endx",&reco_trk_endx,"reco_trk_endx/F",32000); tree->Branch("reco_trk_endy",&reco_trk_endy,"reco_trk_endy/F",32000); tree->Branch("reco_trk_endz",&reco_trk_endz,"reco_trk_endz/F",32000); tree->Branch("reco_trk_endr",&reco_trk_endr,"reco_trk_endr/F",32000); tree->Branch("evt_nplanes",&evt_nplanes,"evt_nplanes/I",32000); tree->Branch("trk_nplanes",&trk_nplanes,"trk_nplanes/I",32000); tree->Branch("slc_nplanes",&slc_nplanes,"slc_nplanes/I",32000); tree->Branch("shw_nplanes_cut",&shw_nplanes_cut,"shw_nplanes_cut/I",32000); tree->Branch("shw_nplanes_raw",&shw_nplanes_raw,"shw_nplanes_raw/I",32000); tree->Branch("shw_ph_cut",&shw_ph_cut,"shw_ph_cut/F",32000); tree->Branch("shw_ph_raw",&shw_ph_raw,"shw_ph_raw/F",32000); tree->Branch("evt_nstrips",&evt_nstrips,"evt_nstrips/I",32000); tree->Branch("trk_nstrips",&trk_nstrips,"trk_nstrips/I",32000); tree->Branch("slc_nstrips",&slc_nstrips,"slc_nstrips/I",32000); tree->Branch("reco_y",&reco_y,"reco_y/F",32000); tree->Branch("reco_Q2",&reco_Q2,"reco_Q2/F",32000); tree->Branch("reco_x",&reco_x,"reco_x/F",32000); tree->Branch("reco_W2",&reco_W2,"reco_W2/F",32000); tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000); tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000); tree->Branch("evtlength",&evtlength,"evtlength/F",32000); tree->Branch("trklength",&trklength,"trklength/F",32000); tree->Branch("shwlength",&shwlength,"shwlength/F",32000); tree->Branch("ntrklike",&ntrklike,"ntrklike/I",32000); tree->Branch("trkend",&trkend,"trkend/I",32000); tree->Branch("trkendu",&trkendu,"trkendu/I",32000); tree->Branch("trkendv",&trkendv,"trkendv/I",32000); tree->Branch("shwbeg",&shwbeg,"shwbeg/I",32000); tree->Branch("shwend",&shwend,"shwend/I",32000); Int_t trknstp; tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000); tree->Branch("trkqp",&trkqp,"trkqp/F",32000); tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000); tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000); tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000); tree->Branch("trkchi2",&trkchi2,"trkchi2/F",32000); tree->Branch("trkndf",&trkndf,"trkndf/I",32000); tree->Branch("trknstp",&trknstp,"trknstp/I",32000); // vertex track back variables (described above and in code) tree->Branch("nvsi", &nvsi, "nvsi/I"); tree->Branch("nvsevt", &nvsevt, "nvsevt/I"); tree->Branch("nvsd", &nvsd, "nvsd/F"); tree->Branch("nvst", &nvst, "nvst/F"); tree->Branch("nvsp", &nvsp, "nvsp/F"); tree->Branch("nvti", &nvti, "nvti/I"); tree->Branch("nvtevt", &nvtevt, "nvtevt/I"); tree->Branch("nvtd", &nvtd, "nvtd/F"); tree->Branch("nvtt", &nvtt, "nvtt/F"); tree->Branch("nvtp", &nvtp, "nvtp/F"); /* tree->Branch("bvsi", &bvsi, "bvsi/I"); tree->Branch("bvsevt", &bvsevt, "bvsevt/I"); tree->Branch("bvsd", &bvsd, "bvsd/F"); tree->Branch("bvst", &bvst, "bvst/F"); tree->Branch("bvsp", &bvsp, "bvsp/F"); tree->Branch("bvti", &bvti, "bvti/I"); tree->Branch("bvtevt", &bvtevt, "bvtevt/I"); tree->Branch("bvtd", &bvtd, "bvtd/F"); tree->Branch("bvtt", &bvtt, "bvtt/F"); tree->Branch("bvtp", &bvtp, "bvtp/F"); */ tree->Branch("evtph", &evtph, "evtph/F"); tree->Branch("evtphgev", &evtphgev, "evtphgev/F"); tree->Branch("evtsigfull", &evtsigfull, "evtsigfull/F"); tree->Branch("evtsigpart", &evtsigpart, "evtsigpart/F"); tree->Branch("trksigfull", &trksigfull, "trksigfull/F"); tree->Branch("trksigpart", &trksigpart, "trksigpart/F"); tree->Branch("evttimemin",&evttimemin,"evttimemin/D"); tree->Branch("evttimemax",&evttimemax,"evttimemax/D"); tree->Branch("evttimeln",&evttimeln,"evttimeln/D"); tree->Branch("evttimeminphc",&evttimeminphc,"evttimeminphc/D"); tree->Branch("evttimemaxphc",&evttimemaxphc,"evttimemaxphc/D"); tree->Branch("evttimelnphc",&evttimelnphc,"evttimelnphc/D"); tree->Branch("earlywph",&earlywph,"earlywph/F"); tree->Branch("earlywphpw",&earlywphpw,"earlywphpw/F"); tree->Branch("earlywphsphere",&earlywphsphere,"earlywphsphere/F"); tree->Branch("ph1000",&ph1000,"ph1000/F"); tree->Branch("ph500",&ph500,"ph500/F"); tree->Branch("ph100",&ph100,"ph100/F"); tree->Branch("lateph1000",&lateph1000,"lateph1000/F"); tree->Branch("lateph500",&lateph500,"lateph500/F"); tree->Branch("lateph100",&lateph100,"lateph100/F"); tree->Branch("lateph50",&lateph50,"lateph50/F"); Int_t largest_event; tree->Branch("largest_event",&largest_event,"largest_event/I"); tree->Branch("nduprs",&nduprs,"ndubrs/I"); tree->Branch("dupphrs",&dupphrs,"dupphrs/F"); tree->Branch("nss", nss, "nss[6]/I"); tree->Branch("ess", ess, "ess[6]/F"); tree->Branch("ntotss", &ntotss, "ntotss/I"); Float_t etu, etv, elu, elv; tree->Branch("etu", &etu, "etu/F"); // transveres u and v energy tree->Branch("etv", &etv, "etv/F"); // uses subshowers tree->Branch("elu", &elu, "elu/F"); // longitudinal energy tree->Branch("elv", &elv, "elv/F"); Float_t swu,swv,wswu,wswv; tree->Branch("swu",&swu, "swu/F");// shower width in u tree->Branch("swv",&swv, "swv/F");// shower width in u tree->Branch("wswu",&wswu, "wswu/F");// shower width in u tree->Branch("wswv",&wswv, "wswv/F");// shower width in u // beam energy weights // 1/E flux, le, z=050, z=100, z=200, z=250 cm Float_t bec_inve, bec_le, bec_050, bec_100, bec_200, bec_250; tree->Branch("bec_inve", &bec_inve, "bec_inve/F"); tree->Branch("bec_le", &bec_le, "bec_le/F"); tree->Branch("bec_050", &bec_050, "bec_050/F"); tree->Branch("bec_100", &bec_100, "bec_100/F"); tree->Branch("bec_200", &bec_200, "bec_200/F"); tree->Branch("bec_250", &bec_250, "bec_250/F"); //li time tree->Branch("litime",&litime,"litime/D"); Int_t rc_boundary; tree->Branch("rc_boundary",&rc_boundary,"rc_boundary/I"); tree->Branch("is_li_sieve",&is_li_sieve,"is_li_sieve/I"); //demux status tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b"); tree->Branch("lowphrat",&lowphrat,"lowphrat/F"); // David's FD data quality cut tree->Branch("dp_fd_quality",&dp_fd_quality,"dp_fd_quality/I"); //anti numu selection float anumupass; tree->Branch("anumupass",&anumupass,"anumupass/F"); //rustem cc pid float rustem_cc_pid; tree->Branch("rustem_cc_pid",&rustem_cc_pid,"rustem_cc_pid/F"); //andy cc pid float andy_cc_pid; tree->Branch("andy_cc_pid",&andy_cc_pid,"andy_cc_pid/F"); //andy's variables float abid_costheta; float abid_eventenergy; float abid_trackcharge; float abid_trackenergy; float abid_tracklikeplanes; float abid_trackphfrac; float abid_trackphmean; float abid_trackqpsigmaqp; float abid_y; tree->Branch("abid_costheta", &abid_costheta, "abid_costheta/F"); tree->Branch("abid_eventenergy", &abid_eventenergy, "abid_eventenergy/F"); tree->Branch("abid_trackcharge", &abid_trackcharge, "abid_trackcharge/I"); tree->Branch("abid_trackenergy", &abid_trackenergy, "abid_trackenergy/i"); tree->Branch("abid_tracklikeplanes", &abid_tracklikeplanes, "abid_tracklikeplanes/I"); tree->Branch("abid_trackphfrac", &abid_trackphfrac, "abid_trackphfrac/F"); tree->Branch("abid_trackphmean", &abid_trackphmean, "abid_trackphmean/F"); tree->Branch("abid_trackqpsigmaqp", &abid_trackqpsigmaqp, "abid_trackqpsigmaqp/F"); tree->Branch("abid_y", &abid_y, "abid_y/F"); float ro_knn_01; float ro_knn_10; float ro_knn_20; float ro_knn_40; tree->Branch("ro_knn_01",&ro_knn_01,"ro_knn01/F"); tree->Branch("ro_knn_10",&ro_knn_10,"ro_knn10/F"); tree->Branch("ro_knn_20",&ro_knn_20,"ro_knn20/F"); tree->Branch("ro_knn_40",&ro_knn_40,"ro_knn40/F"); // beam monitoring variables Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill,dt_beam_mon; //goodbeam==1 if spill meets criteria defined below //(see the line where goodbeam gets set) //beamcnf is a unsigned int as defined in BeamMonSpill::StatusBits // UInt_t beamcnf; //now beamcnf is defined by BeamType in conventions BeamType::BeamType_t beamcnf; //leaving nuTarZ in until we've phased out the summary ntuples Double_t nuTarZ; //goodbeam==0 if spill doesnt meet these criteria Int_t goodbeam; //these next variables weren't available in the old beam monitoring ntuples UInt_t horn_on, target_in, n_batches; Float_t tor101, tr101d, tortgt, trtgtd; Float_t hadmon, mumon1, mumon2, mumon3; Float_t hposbyspill[6]; Float_t vposbyspill[6]; Float_t bibyspill[6]; int batchnumber; tree->Branch("hbw",&hbw,"hbw/D"); tree->Branch("vbw",&vbw,"vbw/D"); tree->Branch("hpos1",&hpos1,"hpos1/D"); tree->Branch("vpos1",&vpos1,"vpos1/D"); tree->Branch("hpos2",&hpos2,"hpos2/D"); tree->Branch("vpos2",&vpos2,"vpos2/D"); tree->Branch("hposbyspill",hposbyspill,"hposbyspill[6]/F"); tree->Branch("vposbyspill",vposbyspill,"vposbyspill[6]/F"); tree->Branch("bibyspill",bibyspill,"bibyspill[6]/F"); tree->Branch("hornI",&hornI,"hornI/D"); tree->Branch("closestspill",&closestspill,"closestspill/D"); tree->Branch("dt_beam_mon",&dt_beam_mon,"dt_beam_mon/D"); tree->Branch("beamcnf",&beamcnf,"beamcnf/i"); tree->Branch("nuTarZ",&nuTarZ,"nuTarZ/D"); tree->Branch("goodbeam",&goodbeam,"goodbeam/I"); tree->Branch("horn_on",&horn_on,"horn_on/i"); tree->Branch("target_in",&target_in,"target_in/i"); tree->Branch("n_batches",&n_batches,"n_batches/i"); tree->Branch("tor101",&tor101,"tor101/F"); tree->Branch("tr101d",&tr101d,"tr101d/F"); tree->Branch("tortgt",&tortgt,"tortgt/F"); tree->Branch("trtgtd",&trtgtd,"trtgtd/F"); tree->Branch("hadmon",&hadmon,"hadmon/F"); tree->Branch("mumon1",&mumon1,"mumon1/F"); tree->Branch("mumon2",&mumon2,"mumon2/F"); tree->Branch("mumon3",&mumon3,"mumon3/F"); tree->Branch("batchnumber",&batchnumber,"batchnumber/I"); //coil status Short_t coil_is_ok; Short_t coil_is_reverse; Short_t coil_stat; Float_t coil_current; tree->Branch("coil_is_ok",&coil_is_ok,"coil_is_ok/S"); tree->Branch("coil_is_reverse",&coil_is_reverse,"coil_is_reverse/S"); tree->Branch("coil_stat",&coil_stat,"coil_stat/S"); tree->Branch("coil_current",&coil_current,"coil_current/F"); //gnumiflux info Int_t fluxrun; // gnumi flux run number Int_t fluxevtno; // gnumi flux event number Float_t nudxdz; // dx/dz slope, neutrino rndm decay Float_t nudydz; // dy/dz slope, neutrino rndm decay Float_t nupz; // Pz of neutrino (GeV) rndm decay Float_t nuenergy; // E(neutrino) (GeV) rndm decay Float_t nudxdznear; // dx/dz slope, neutrino NearDet center Float_t nudydznear; // dy/dz slope, neutrino NearDet center Float_t nuenergynear; // E(neutrino) (GeV) NearDet center Float_t nwtnear; // Weight of nu NearDet center Float_t nudxdzfar; // dx/dz slope, neutrino FarDet center Float_t nudydzfar; // dy/dz slope, neutrino FarDet center Float_t nuenergyfar; // E(neutrino) (GeV) FarDet center Float_t nwtfar; // Weight of nu FarDet center Int_t norig; // (ignore) Int_t ndecay; // Tag of decay mode Int_t ntype; // Neutrino type (translated to PDG) Float_t vx; // x vertex of hadron (cm) Float_t vy; // y vertex of hadron (cm) Float_t vz; // z vertex of hadron (cm) Float_t pdpx; // nu parent px at decay point Float_t pdpy; // nu parent py at decay point Float_t pdpz; // nu parent pz at decay point Float_t ppdxdz; // nu parent slope at production Float_t ppdydz; // nu parent slope at production Float_t pppz; // nu parent pz at production Float_t ppenergy; // nu parent energy at production Int_t ppmedium; // GEANT medium of nu parent at parent production Int_t ptype; // nu parent type (translated to PDG) Float_t ppvx; // nu parent production vtx x Float_t ppvy; // nu parent production vtx y Float_t ppvz; // nu parent production vtx z Float_t muparpx; // if parent=mu, hadron parent px Float_t muparpy; // if parent=mu, hadron parent py Float_t muparpz; // if parent=mu, hadron parent pz Float_t mupare; // if parent=mu, hadron parent energy Float_t necm; // E(nu) in parent cm Float_t nimpwt; // importance weight Float_t xpoint; // (unused) Float_t ypoint; // (unused) Float_t zpoint; // (unused) Float_t tvx; // target exit point (x) of parent Float_t tvy; // target exit point (y) of parent Float_t tvz; // target exit point (z) of parent Float_t tpx; // parent px at target exit Float_t tpy; // parent py at target exit Float_t tpz; // parent pz at target exit Int_t tptype; // parent particle type (translated to PDG) Int_t tgen; // parent generation in cascade tree->Branch("fluxrun",&fluxrun,"fluxrun/I"); tree->Branch("fluxevtno",&fluxevtno,"fluxevtno/I"); tree->Branch("nudxdz",&nudxdz,"nudxdz/F"); tree->Branch("nudydz",&nudydz,"nudydz/F"); tree->Branch("nupz",&nupz,"nupz/F"); tree->Branch("nuenergy",&nuenergy,"nuenergy/F"); tree->Branch("nudxdznear",&nudxdznear,"nudxdznear/F"); tree->Branch("nudydznear",&nudydznear,"nudydznear/F"); tree->Branch("nuenergynear",&nuenergynear,"nuenergynear/F"); tree->Branch("nwtnear",&nwtnear,"nwtnear/F"); tree->Branch("nudxdzfar",&nudxdzfar,"nudxdzfar/F"); tree->Branch("nudydzfar",&nudydzfar,"nudydzfar/F"); tree->Branch("nuenergyfar",&nuenergyfar,"nuenergyfar/F"); tree->Branch("nwtfar",&nwtfar,"nwtfar/F"); tree->Branch("norig",&norig,"norig/I"); tree->Branch("ndecay",&ndecay,"ndecay/I"); tree->Branch("ntype",&ntype,"ntype/I"); tree->Branch("vx",&vx,"vx/F"); tree->Branch("vy",&vy,"vy/F"); tree->Branch("vz",&vz,"vz/F"); tree->Branch("pdpx",&pdpx,"pdpx/F"); tree->Branch("pdpy",&pdpy,"pdpy/F"); tree->Branch("pdpz",&pdpz,"pdpz/F"); tree->Branch("ppdxdz",&ppdxdz,"ppdxdz/F"); tree->Branch("ppdydz",&ppdydz,"ppdydz/F"); tree->Branch("pppz",&pppz,"pppz/F"); tree->Branch("ppenergy",&ppenergy,"ppenergy/F"); tree->Branch("ppmedium",&ppmedium,"ppmedium/I"); tree->Branch("ptype",&ptype,"ptype/I"); tree->Branch("ppvx",&ppvx,"ppvx/F"); tree->Branch("ppvy",&ppvy,"ppvy/F"); tree->Branch("ppvz",&ppvz,"ppvz/F"); tree->Branch("muparpx",&muparpx,"muparpx/F"); tree->Branch("muparpy",&muparpy,"muparpy/F"); tree->Branch("muparpz",&muparpz,"muparpz/F"); tree->Branch("mupare",&mupare,"mupare/F"); tree->Branch("necm",&necm,"necm/F"); tree->Branch("nimpwt",&nimpwt,"nimpwt/F"); tree->Branch("xpoint",&xpoint,"xpoint/F"); tree->Branch("ypoint",&ypoint,"ypoint/F"); tree->Branch("zpoint",&zpoint,"zpoint/F"); tree->Branch("tvx",&tvx,"tvx/F"); tree->Branch("tvy",&tvy,"tvy/F"); tree->Branch("tvz",&tvz,"tvz/F"); tree->Branch("tpx",&tpx,"tpx/F"); tree->Branch("tpy",&tpy,"tpy/F"); tree->Branch("tpz",&tpz,"tpz/F"); tree->Branch("tptype",&tptype,"tptype/I"); tree->Branch("tgen",&tgen,"tgen/I"); //pot counting float pottor101=0.; float pottortgt=0.; float pot=0.; int NRUNS=0; int NSNARLS=0; TTree *pottree = new TTree("pottree","pottree"); pottree->Branch("pot",&pot,"pot/F"); pottree->Branch("pottor101",&pottor101,"pottor101/F"); pottree->Branch("pottortgt",&pottortgt,"pottortgt/F"); pottree->Branch("nruns",&NRUNS,"nruns/I"); pottree->Branch("nsnarls",&NSNARLS,"nsnarls/I"); //make a BMSpillAna so we can tell if it's goodbeam BMSpillAna bmana; bmana.UseDatabaseCuts(); BeamEnergyCalculator ecalc(&fBECConfig); double bec_high=0; double bec_low=0; ecalc.GetStandardConfig()->Print(); ecalc.GetStandardConfig()->Get("beam:high_energy_limit",bec_high); ecalc.GetStandardConfig()->Get("beam:low_energy_limit",bec_low); // mark's QE id // MadQEID qeid; // qeid.ReadPDFs("/afs/fnal.gov/files/data/minos/d04/kordosky/madqe_files/QEpdfs.root"); // bool beamdb=false; // bool checkeddb=false; int lastrun=-1; for(int ii=0;ii<Nentries;ii++){ if(ii%1000==0) std::cout << float(ii)*100./float(Nentries) << "% done" << std::endl; if(!this->GetEntry(ii)) continue; NSNARLS++; // fill some histograms for mc acceptance if(file_is_mc) FillMCHists(&save); //give the snarl to the PhysicsNtuple interface bool anumuready = phnti->FillSnarl(strecord); if(anumuready){ // cout<<"ITHINK THE STUPID INTERFACE THING WORKED WITH THE SNARL"<<endl; } VldContext vc = ntpHeader->GetVldContext(); nevent = eventSummary->nevent; run = ntpHeader->GetRun(); if(run!=lastrun){ NRUNS++; lastrun = run; } //get detector type: const Detector::Detector_t det = ntpHeader->GetVldContext().GetDetector(); snarl = ntpHeader->GetSnarl(); subrun = ntpHeader->GetSubRun(); vld_sec = ntpHeader->GetVldContext().GetTimeStamp().GetSeconds(); trgsrc = ntpHeader->GetTrigSrc(); spilltype = ntpHeader->GetRemoteSpillType(); trgtime = eventSummary->trigtime; litime = eventSummary->litime; is_li_sieve = LISieve::IsLI(*strecord); coil_stat = detStatus->coilstatus; coil_is_ok = 1; coil_is_reverse = 0; coil_current = 0; if(!file_is_mc) { coil_is_ok = CoilTools::IsOK(vc); coil_is_reverse = CoilTools::IsReverse(vc); coil_current = CoilTools::CoilCurrent(vc).first; } dmx_stat=0; //figure out demux status dmx_stat+=(dmxStatus->nonphysicalfail); dmx_stat+=((dmxStatus->validplanesfail<<1)); dmx_stat+=((dmxStatus->vertexplanefail<<2)); //figure out how much of snarl ph was deposited as low ph lowphrat=0.; lowphrat = ComputeLowPHRatio(); //is this record in a time period defined as "bad" by Dave Petyt? // 1=good, 0=bad if(det==Detector::kFar){ // dp_fd_quality=david_passfail(vld_sec); // dp_fd_quality=passfail(vld_sec); dp_fd_quality=DataUtil::IsGoodFDData(strecord); } else if(det==Detector::kNear){ dp_fd_quality=1; } // look through "events" and characterize how close together // their vertices are // idea is to flag runt events caused by late activity near vertex NearbyVertexFinder nby(nevent); nby.ProcessEntry(this); // run through and tag event with largest signal int lgst_evt_idx=-1; float big_ph=0.0; for(int i=0;i<eventSummary->nevent;i++){ if(!LoadEvent(i)) continue; //no event found if(ntpEvent->ph.sigcor > big_ph){ big_ph=ntpEvent->ph.sigcor; lgst_evt_idx=i; } } detector = -1; if(det==Detector::kFar){ detector = 2; } else if(det==Detector::kNear){ detector = 1; } // deal with beam type as best we can BeamType::BeamType_t current_beam=BeamType::kUnknown; if(file_is_mc){ // for MC set beam to fMCbeam current_beam=fMCBeam; // for data we will try to deduce from beam monitoring } // first thing, if data do beam monitoring hbw=vbw=hpos1=vpos1=hpos2=vpos2=hornI=closestspill=dt_beam_mon=-999.0; for(int scnt=0;scnt<6;scnt++){ hposbyspill[scnt]=0.; vposbyspill[scnt]=0.; bibyspill[scnt]=0.; } beamcnf=BeamType::kUnknown; nuTarZ=0.0; goodbeam=0; horn_on=target_in=n_batches=0; tor101=tr101d=tortgt=trtgtd=0.0; hadmon=mumon1=mumon2=mumon3=0.0; VldTimeStamp vts = vc.GetTimeStamp(); if(!file_is_mc){ // if(beamdb){ BDSpillAccessor &sa = BDSpillAccessor::Get(); const BeamMonSpill *bs = sa.LoadSpill(vts); hbw=bs->fProfWidX*1000.;//make it in mm vbw=bs->fProfWidY*1000.;//make it in mm hpos1=bs->fTargProfX*1000.;//make it in mm vpos1=bs->fTargProfY*1000.;//make it in mm n_batches=0; float btint=0.; hpos2=0.; vpos2=0.; for(int bt=0;bt<6;bt++){ //cout<<"batch "<<bt<<" x "<<bs->fTargBpmX[bt]<<" y "<<bs->fTargBpmY[bt] // <<" int "<<bs->fBpmInt[bt]<<endl; hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]); vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]); hposbyspill[bt]=bs->fTargBpmX[bt]; vposbyspill[bt]=bs->fTargBpmY[bt]; bibyspill[bt]=bs->fBpmInt[bt]; if(bs->fBpmInt[bt]>0.001){ n_batches++; } btint+=bs->fBpmInt[bt]; } if(btint!=0){ hpos2=hpos2*1000./btint;//make it in mm vpos2=vpos2*1000./btint;//make it in mm } else{ hpos2=-999.; vpos2=-999.; } hornI=bs->fHornCur; SpillTimeFinder &sfind= SpillTimeFinder::Instance(); closestspill = sfind.GetTimeToNearestSpill(vc); // beamcnf =bs->GetStatusBits().beam_type; beamcnf =bs->BeamType(); nuTarZ=0.; horn_on = bs->GetStatusBits().horn_on; target_in = bs->GetStatusBits().target_in; tor101=bs->fTor101; tr101d=bs->fTr101d; tortgt=bs->fTortgt; trtgtd=bs->fTrtgtd; hadmon=bs->fHadInt; mumon1=bs->fMuInt1; mumon2=bs->fMuInt2; mumon3=bs->fMuInt3; goodbeam=0; //use Mark D. BMSpillAna to set good beam bmana.SetSpill(*bs); // bmana.SetTimeDiff(closestspill); // Mark D. suggests that we do the following rather // than the line above VldTimeStamp delta_vts=vts-bs->SpillTime(); dt_beam_mon=delta_vts.GetSeconds(); bmana.SetTimeDiff(delta_vts.GetSeconds()); if(bmana.SelectSpill()){ goodbeam=1; } else{ goodbeam=0; } if(fDataUseMCBeam==true){ //allow us to force the data to use the MC beam current_beam=fMCBeam; } else{ // attempt to deduce beam type from beam monitoring // current_beam = BeamType::FromBeamMon(beamcnf); current_beam = beamcnf; } // cout<<"Current beam "<<current_beam<<" "<<BeamType::AsString(current_beam)<<endl; //don't count fake spills if(goodbeam==1&&spilltype!=3){ //don't count spills where ND coil not working if(det==Detector::kFar||coil_is_ok==1){ pot+=tortgt; pottor101+=tor101; pottortgt+=tortgt; } } } // Write a special record for each spill! event=-1; ntrack=-1; nshower=-1; reco_emu=0.0; reco_enu=0.0; reco_eshw=0.0; reco_eshw_uncorr=0.0; true_enu=0.0; true_emu=0.0; true_eshw=0.0; pass=0; is_fid=0; is_fid_2007=0; anumupass=0; rustem_cc_pid=0; ro_knn_01=0.; ro_knn_10=0.; ro_knn_20=0.; ro_knn_40=0.; andy_cc_pid=-1.; abid_costheta =0.; abid_eventenergy =0.; abid_trackcharge =0; abid_trackenergy =0; abid_tracklikeplanes =0; abid_trackphfrac =0.; abid_trackphmean =0.; abid_trackqpsigmaqp =0.; abid_y =0.; tree->Fill(); if(!openedabpidfile){ string base=""; base=getenv("SRT_PRIVATE_CONTEXT"); if(base!=""&&base!="."){ // check if directory exists in SRT_PRIVATE_CONTEXT std::string path = base + "/Mad/data"; void *dir_ptr = gSystem -> OpenDirectory(path.c_str()); if(!dir_ptr){ // if it doesn't exist then use SRT_PUBLIC_CONTEXT base=getenv("SRT_PUBLIC_CONTEXT"); } } else{ base=getenv("SRT_PUBLIC_CONTEXT"); } if(base=="") { cout<<"No SRT_PUBLIC_CONTEXT set"<<std::endl; assert(false); } base+="/Mad/data"; string fname=base; string sdet=""; if(det==Detector::kNear){ sdet="near"; } else if(det==Detector::kFar){ sdet="far"; } else{ cout<<"Could not set ab pid file, defaulting to far"<<endl; sdet="far"; } string rmc=""; if(ReleaseType::IsCedar(reltype)&&mcver=="daikon"){ rmc="cedar_daikon"; } else{ cout<<"Dont know reco/mc versions " <<"defaulting to cedar_daikon"<<endl; rmc="cedar_daikon"; } string sbeam=""; switch (current_beam){ case BeamType::kL010z000i: sbeam="le0"; break; case BeamType::kL010z170i: sbeam="le170"; break; case BeamType::kL010z185i: sbeam="le"; break; case BeamType::kL010z200i: sbeam="le200"; break; case BeamType::kL100z200i: sbeam="pme"; break; case BeamType::kL150z200i: sbeam="pmhe"; break; case BeamType::kL250z200i: sbeam="phe"; break; case BeamType::kLE: sbeam="le"; break; default: cout<<"Don't know beam type " <<" defaulting to LE"<<endl; sbeam="le"; break; } fname+="/ab_pdf_"+sdet+"_"+sbeam+"_"+rmc+".root"; cout<<"Reading AB PDFS from "<<fname<<endl; abid.ReadPDFs(fname.c_str()); openedabpidfile=true; } if(nevent==0) continue; //fill tree once for each reconstructed event for(int i=0;i<eventSummary->nevent;i++){ if(!LoadEvent(i)) continue; //no event found //set event index: event = i; ntrack = ntpEvent->ntrack; nshower = ntpEvent->nshower; nstrip = ntpEvent->nstrip; //anti nu selection anumupass=0; rustem_cc_pid=0; ro_knn_01=0.; ro_knn_10=0.; ro_knn_20=0.; ro_knn_40=0.; andy_cc_pid=-1.; abid_costheta =0.; abid_eventenergy =0.; abid_trackcharge =0; abid_trackenergy =0; abid_tracklikeplanes =0; abid_trackphfrac =0.; abid_trackphmean =0.; abid_trackqpsigmaqp =0.; abid_y =0.; if(anumuready){ // cout<<"Filling event rustem stuff"<<endl; anumupass = phnti->GetVar("numubar", ntpEvent); rustem_cc_pid = phnti->GetVar("knn_pid", ntpEvent); ro_knn_01 = phnti->GetVar("knn_01", ntpEvent); ro_knn_10 = phnti->GetVar("knn_10", ntpEvent); ro_knn_20 = phnti->GetVar("knn_20", ntpEvent); ro_knn_40 = phnti->GetVar("knn_40", ntpEvent); } //zero all tree values true_enu = 0; true_emu = 0; true_eshw = 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; true_pitt_fid=0; true_trk_cmplt=0.; resnum=0; npp=npm=npz=np=nn=0; epp=epm=epz=ep=en=0.0; fluxrun=0; fluxevtno=0; nudxdz=0.; nudydz=0.; nupz=0.; nuenergy=0.; nudxdznear=0.; nudydznear=0.; nuenergynear=0.; nwtnear=0.; nudxdzfar=0.; nudydzfar=0.; nuenergyfar=0.; nwtfar=0.; norig=0; ndecay=0; ntype=0; vx=0.; vy=0.; vz=0.; pdpx=0.; pdpy=0.; pdpz=0.; ppdxdz=0.; ppdydz=0.; pppz=0.; ppenergy=0.; ppmedium=0; ptype=0; ppvx=0.; ppvy=0.; ppvz=0.; muparpx=0.; muparpy=0.; muparpz=0.; mupare=0.; necm=0.; nimpwt=0.; xpoint=0.; ypoint=0.; zpoint=0.; tvx=0.; tvy=0.; tvz=0.; tpx=0.; tpy=0.; tpz=0.; tptype=0; tgen=0; flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0; rc_boundary=0; mcevent = -1; is_fid = 0;is_fid_2007 = 0; is_cev = 0; sr_contained = 0; is_mc = 0; pass_basic = 0; is_pitt_fid=0; pitt_evt_class=0; pitt_trk_qual=0; nd_radial_fid=0; is_trk_fid=0; duvvtx=0; pid0 = 0; pid1 = 0; pid2 = 0; pass = 0; reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_wtd_eshw=0; reco_eshw_uncorr=0; reco_wtd_eshw_uncorr=0; reco_vtx_eshw=0; reco_vtx_wtd_eshw=0; reco_qe_enu = 0; reco_mushw=0; reco_qe_q2 = 0; reco_dircosneu = 0; reco_dircosx = 0;; reco_dircosy = 0;; reco_dircosz = 0; reco_vtxx = reco_vtxy = reco_vtxz = reco_vtxr=0; reco_trk_vtxx = reco_trk_vtxy = reco_trk_vtxz = reco_trk_vtxr=0; reco_shw_vtxx = reco_shw_vtxy = reco_shw_vtxz = reco_shw_vtxr=0; reco_trk_endx = reco_trk_endy = reco_trk_endz = reco_trk_endr=0; evt_nstrips = slc_nstrips = trk_nstrips=0; evt_nplanes = slc_nplanes = trk_nplanes=0; shw_nplanes_cut=shw_nplanes_raw=0; shw_ph_cut=shw_ph_raw=0.0; evtlength = trklength = shwlength=trknstp=0; trkphfrac = trkphplane = 0; ntrklike= trkend=trkendu=trkendv=shwbeg=shwend=0; trkmomrange = 0; trkqp = 0; trkeqp = 0; trkchi2=0.0; trkndf=0; pass_track=0; emu_meth=0; trk_fit_pass =0; // kinematic quantities always positive // convention: -1 indicates unfilled variable (NC events) reco_y = reco_x = reco_Q2 = reco_W2 = -1; // vertex back tracking nvsi=nvti=bvsi=bvti=-1; nvsevt=bvsevt=nvtevt=bvtevt=-1; nvsd=nvtd=bvsd=bvtd=999999; nvst=bvst=nvtt=bvtt=1e6; nvsp=nvtp=bvsp=bvtp=-1.0; evtph=0; evtphgev=0; evtsigfull=evtsigpart=trksigfull=trksigpart=0; evttimemin=evttimemax=evttimeln=0.; evttimeminphc=evttimemaxphc=evttimelnphc=0.; earlywph=earlywphpw=earlywphsphere=0.; ph1000=ph500=ph100=0.; lateph1000=lateph500=lateph100=lateph50=0.; nduprs=0; dupphrs=0.; nu_px=nu_py=nu_pz=tar_px=tar_py=tar_pz=tar_e=0.0; inu=0; inunoosc=0; ntotss=0; for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;} med_qe_pid=-1001.0; med_qe_pass=0; niki_cc_pid=-999; dave_cc_pid=-999; //check for a corresponding mc event // first, looks to match reconstructed track // then, looks to match reconstructed shower if(LoadTruthForRecoTH(i,mcevent)) { // should flag with a "bad truth word" // in case the function above returns false //true info for tree: is_mc = 1; nucleus = Nucleus(mcevent); flavour = Flavour(mcevent); initial_state = Initial_state(mcevent); process = IResonance(mcevent); cc_nc = IAction(mcevent); true_enu = TrueNuEnergy(mcevent); 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); //flux info if(LoadTruth(mcevent)){ fluxrun=ntpTruth->flux.fluxrun; fluxevtno=ntpTruth->flux.fluxevtno; nudxdz=ntpTruth->flux.ndxdz; nudydz=ntpTruth->flux.ndydz; nupz=ntpTruth->flux.npz; nuenergy=ntpTruth->flux.nenergy; nudxdznear=ntpTruth->flux.ndxdznear; nudydznear=ntpTruth->flux.ndydznear; nuenergynear=ntpTruth->flux.nenergynear; nwtnear=ntpTruth->flux.nwtnear; nudxdzfar=ntpTruth->flux.ndxdzfar; nudydzfar=ntpTruth->flux.ndydzfar; nuenergyfar=ntpTruth->flux.nenergyfar; nwtfar=ntpTruth->flux.nwtfar; norig=ntpTruth->flux.norig; ndecay=ntpTruth->flux.ndecay; ntype=ntpTruth->flux.ntype; vx=ntpTruth->flux.vx; vy=ntpTruth->flux.vy; vz=ntpTruth->flux.vz; pdpx=ntpTruth->flux.pdpx; pdpy=ntpTruth->flux.pdpy; pdpz=ntpTruth->flux.pdpz; ppdxdz=ntpTruth->flux.ppdxdz; ppdydz=ntpTruth->flux.ppdydz; pppz=ntpTruth->flux.pppz; ppenergy=ntpTruth->flux.ppenergy; ppmedium=ntpTruth->flux.ppmedium; ptype=ntpTruth->flux.ptype; ppvx=ntpTruth->flux.ppvx; ppvy=ntpTruth->flux.ppvy; ppvz=ntpTruth->flux.ppvz; muparpx=ntpTruth->flux.muparpx; muparpy=ntpTruth->flux.muparpy; muparpz=ntpTruth->flux.muparpz; mupare=ntpTruth->flux.mupare; necm=ntpTruth->flux.necm; nimpwt=ntpTruth->flux.nimpwt; xpoint=ntpTruth->flux.xpoint; ypoint=ntpTruth->flux.ypoint; zpoint=ntpTruth->flux.zpoint; tvx=ntpTruth->flux.tvx; tvy=ntpTruth->flux.tvy; tvz=ntpTruth->flux.tvz; tpx=ntpTruth->flux.tpx; tpy=ntpTruth->flux.tpy; tpz=ntpTruth->flux.tpz; tptype=ntpTruth->flux.tptype; tgen=ntpTruth->flux.tgen; } if(ntpTHTrack){ if(ntpTHTrack->neumc==mcevent){ true_trk_cmplt=ntpTHTrack->completeslc; } } Float_t* true_p = TrueNuMom(mcevent); if(true_p){ nu_px=true_p[0]; nu_py=true_p[1]; nu_pz=true_p[2]; delete[] true_p; true_p=0; } true_p = Target4Vector(mcevent); if(true_p){ tar_px=true_p[0]; tar_py=true_p[1]; tar_pz=true_p[2]; tar_e=true_p[3]; delete[] true_p; true_p=0; } Float_t *vtx = TrueVtx(mcevent); true_vtxx = vtx[0]; true_vtxy = vtx[1]; true_vtxz = vtx[2]; float true_vtxu = (true_vtxx + true_vtxy)/sqrt(2.0); float true_vtxv = (true_vtxy - true_vtxx)/sqrt(2.0); true_pitt_fid=PittInFidVol(true_vtxx, true_vtxy, true_vtxz, true_vtxu, true_vtxv); resnum = HadronicFinalState(mcevent); npp = NumFSPion(mcevent,+1); npm = NumFSPion(mcevent,-1); npz = NumFSPiZero(mcevent); np = NumFSProton(mcevent); nn = NumFSNeutron(mcevent); epp = TotFSPionEnergy(mcevent,+1); epm = TotFSPionEnergy(mcevent,-1); epz = TotFSPiZeroEnergy(mcevent); ep = TotFSProtonEnergy(mcevent); en = TotFSNeutronEnergy(mcevent); // beam energy weights // inu=0; inu = INu(mcevent); inunoosc = INuNoOsc(mcevent); bec_inve=bec_le=bec_050=bec_100=bec_200=bec_250=1.0; /* bec_inve=ecalc.GetWeight(BeamType::kInverseE,det,inu, true_enu, bec_high,bec_low); bec_le=ecalc.GetWeight(BeamType::kLE,det,inu, true_enu, bec_high,bec_low); bec_050=ecalc.GetWeight(BeamType::k050,det,inu, true_enu, bec_high,bec_low); bec_100=ecalc.GetWeight(BeamType::k100,det,inu, true_enu, bec_high,bec_low); bec_200=ecalc.GetWeight(BeamType::k200,det,inu, true_enu, bec_high,bec_low); bec_250=ecalc.GetWeight(BeamType::k250,det,inu, true_enu, bec_high,bec_low); */ delete [] vtx; } // is the event in the fiducial volume // MadQuantities::IsFidVtxEvt() // this is used by Dave/Niki analysis // provisional, will update to using track vertices // if there is a track. is_fid_2007 = IsFidVtxEvt(event); is_fid = IsFid_2008(event); if(det==Detector::kFar){ if(IsFidFD_AB(event)){ is_fid_2007 = 1; } else{ is_fid_2007=0; } } // is this event on a readout crate boundary rc_boundary=0; if(det==Detector::kFar){ // if(FDRCBoundary(eventSummary->planeall.beg,eventSummary->planeall.end)) rc_boundary=1; // if(FDRCBoundary(eventSummary->plane.beg,eventSummary->plane.end)) rc_boundary=1; rc_boundary=FDRCBoundary(); } // is this event the largest (sicor) in the spill if(event==lgst_evt_idx) largest_event=1; else largest_event=0; if(PassBasicCuts(event)) { // this does nothing. //reco info for tree: pass_basic = 1; //index will -1 if no track/shower present in event int track_index = -1; // track spanning largest number of planes LoadLargestTrackFromEvent(i,track_index); if(track_index>-1&&det==Detector::kNear){ //continue to use PRL cuts for ND fid. volume //FD fid vol has been updated and is computed above is_fid_2007=InFidVol(ntpTrack->vtx.x,ntpTrack->vtx.y,ntpTrack->vtx.z); } pass_track = PassTrackCuts(track_index); //methods should all be protected against index = -1 // reco_emu = RecoMKMuEnergy(emu_meth,track_index, // !file_is_mc,det); reco_emu = RecoMuEnergyNew(0,track_index,vc, reltype,EnergyCorrections::kDefault); int shower_index = -1; // LoadLargestShowerFromEvent(i,shower_index); // Now use 0th shower (set by Jim) with cuts to remove brems LoadShower_Jim(i,shower_index,det); // shower with the best match to the vertex // reco_eshw = RecoShwEnergy(shower_index, det,1,!file_is_mc); reco_eshw=RecoShwEnergyNew(shower_index,0,vc, reltype,EnergyCorrections::kDefault); reco_wtd_eshw=RecoShwEnergyNew(shower_index,1,vc, reltype,EnergyCorrections::kDefault); // cout<<"reco_eshw output"<<reco_eshw<<endl; // cout<<"VLD context "<<endl; // vc.Print(); // cout<<endl; if(shower_index>-1) { reco_eshw_uncorr = ntpShower->shwph.linCCgev; reco_wtd_eshw_uncorr = ntpShower->shwph.wtCCgev; } reco_vtx_eshw = reco_eshw; reco_vtx_wtd_eshw = reco_wtd_eshw; // provisional, to be updated later reco_enu = reco_emu + reco_eshw; // event energy according to offline evtphgev = ntpEvent->ph.gev; if(reco_enu>0){ reco_qe_enu = RecoQENuEnergy(track_index); reco_qe_q2 = RecoQEQ2(track_index); // note, NtpSRFiducial isn't filled for ND events // is_cev determines if the track is a stopper or punch-through // used in Dave/Niki analysis is_cev = IsFidAll(track_index); // use Pitt group's fiducial cuts for ND events if(detector==1){ // for now, use event vertex // but maybe should use track vertex... // if vertex finder works well it shouldn't matter... is_pitt_fid = PittInFidVol(ntpEvent->vtx.x, ntpEvent->vtx.y, ntpEvent->vtx.z, ntpEvent->vtx.u, ntpEvent->vtx.v); nd_radial_fid = NDRadialFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y, ntpEvent->vtx.z); } else{ is_pitt_fid = InFidVol(); } // have already checked that ntpEvent exists reco_vtxx = ntpEvent->vtx.x; reco_vtxy = ntpEvent->vtx.y; reco_vtxz = ntpEvent->vtx.z; // check if distance from vertex to shower is too large // if so, set reco_eshw =0 // idea is to handle case of no vertex shower but runt downstream if(shower_index!=-1){ // bool shw_ok=true; const float max_shw_dist=14*0.06; // about 2 interaction lengths float dist_z = fabs(reco_vtxz - ntpShower->vtx.z); float dist_r = sqrt( pow(ntpEvent->vtx.u-ntpShower->vtx.u,2) + pow(ntpEvent->vtx.v-ntpShower->vtx.v,2) ); if( (dist_z + dist_r) > max_shw_dist ){ reco_vtx_eshw=0; reco_vtx_wtd_eshw=0; } } // radial position of vertex if(detector==1){ reco_vtxr = pow(reco_vtxx-kXcenter,2) + pow(reco_vtxy-kYcenter,2); reco_vtxr=sqrt(reco_vtxr); } else reco_vtxr = sqrt (pow(reco_vtxx,2)+pow(reco_vtxy,2)); // compute signal in fully and partially instrumented regions SigInOut(track_index,evtsigfull,evtsigpart,trksigfull,trksigpart); if(detector==2){ //if far det, all signal is in fully instrumented region evtsigfull=ntpEvent->ph.sigcor; evtsigpart=0.; if(track_index!=-1){ trksigfull=ntpTrack->ph.sigcor; } else{ trksigfull=0.; } trksigpart=0.; } // angle reconstruction if(track_index>-1) { reco_dircosx = ntpTrack->vtx.dcosx; reco_dircosy = ntpTrack->vtx.dcosy; } reco_dircosz = RecoMuDCosZ(track_index); // event vertex float vtx_array[3]; vtx_array[0]=ntpEvent->vtx.x; vtx_array[2]=ntpEvent->vtx.y; vtx_array[2]=ntpEvent->vtx.z; if(detector==1) reco_dircosneu = RecoMuDCosNeuND(track_index, vtx_array); else reco_dircosneu = RecoMuDCosNeuFD(track_index, vtx_array); evtlength = ntpEvent->plane.end-ntpEvent->plane.beg; if(track_index!=-1){ //check track is present before using ntpTrack trklength = ntpTrack->plane.end-ntpTrack->plane.beg; ntrklike =ntpTrack->plane.ntrklike; trkend = ntpTrack->plane.end; trkendu = ntpTrack->plane.endu; trkendv = ntpTrack->plane.endv; trknstp = ntpTrack->nstrip; trkmomrange = ntpTrack->momentum.range; trk_fit_pass = ntpTrack->fit.pass; // trkmomrange = EnergyCorrections::CorrectMomentumFromRange(trkmomrange,!file_is_mc,det); //now use FullyCorrect Versions trkmomrange = EnergyCorrections::FullyCorrectMomentumFromRange(trkmomrange,vc,reltype); trkqp = ntpTrack->momentum.qp; if(trkqp!=0){ // trkqp = 1.0/EnergyCorrections::CorrectSignedMomentumFromCurvature(1.0/trkqp, trkqp = 1.0/EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(1.0/trkqp,vc,reltype); } trkeqp = ntpTrack->momentum.eqp; trkphfrac = ntpTrack->ph.sigcor/ntpEvent->ph.sigcor; trkphplane = ntpTrack->ph.mip/ntpTrack->plane.n; trkchi2 = ntpTrack->fit.chi2; trkndf = ntpTrack->fit.ndof; reco_trk_vtxx = ntpTrack->vtx.x; reco_trk_vtxy = ntpTrack->vtx.y; reco_trk_vtxz = ntpTrack->vtx.z; reco_trk_endx = ntpTrack->end.x; reco_trk_endy = ntpTrack->end.y; reco_trk_endz = ntpTrack->end.z; sr_contained = ntpTrack->contained; // compute CC event class as defined by Pitt group pitt_evt_class = PittTrkContained(ntpTrack->end.x, ntpTrack->end.y, ntpTrack->end.z, ntpTrack->end.u, ntpTrack->end.v); duvvtx=abs(ntpTrack->plane.begu - ntpTrack->plane.begv); // based on the event class, go back and recompute // the muon and neutrino energy // modified here, now use is_cev=IsFidAll() int do_meth=0; // emu reco method we should use if( is_cev==1 ) do_meth=2; // stoppers --> range else if( is_cev==0 ) do_meth=1; // punch-through --> curvature // with the advent of trk.fit.pass==0 being accepted // it is possible that we have to modify do_meth to // account for qp==0 // in this case we have to use the range if(trkqp==0) do_meth=2; if(do_meth==1 || do_meth==2){ // reco_emu = RecoMKMuEnergy(do_meth,track_index,!file_is_mc,det); reco_emu = RecoMuEnergyNew(do_meth,track_index, vc,reltype, EnergyCorrections::kDefault); // update reco_enu to account for possible change // in reco_emu reco_enu = reco_eshw + reco_emu; } emu_meth=do_meth; // was forgetting about this //check to see if track vertex is in fid. vol if(detector==1){ is_trk_fid = PittInFidVol(ntpTrack->vtx.x, ntpTrack->vtx.y, ntpTrack->vtx.z, ntpTrack->vtx.u, ntpTrack->vtx.v); } else{ is_trk_fid = InFidVol(ntpTrack->vtx.x,ntpTrack->vtx.y, ntpTrack->vtx.z); } // look for showers along muon track reco_mushw=MuonShowerEnergy(ntpEvent, ntpTrack); // pitt tracking quality cuts Int_t temp_ndof = ntpTrack->fit.ndof; if(temp_ndof==0) temp_ndof=1; if( (ntpTrack->fit.pass>0) && (ntpTrack->fit.chi2/temp_ndof <20 ) && abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 ) pitt_trk_qual=1; const double M=(0.93827 + 0.93957)/2.0; // <nucleon mass> if(reco_emu>0) reco_Q2 = 2*reco_enu*reco_emu*(1.0 - reco_dircosneu); reco_W2 = M*M - reco_Q2 + 2*M*reco_eshw; if(reco_enu>0) reco_y = reco_eshw/reco_enu; if(reco_eshw>0 && reco_Q2>0) reco_x = reco_Q2/(2*M*reco_eshw); const double muM=0.10566; // muon mass const double muP=sqrt(reco_emu*reco_emu -muM*muM); reco_qe_enu = (M*reco_emu - muM*muM/2.0) /(M - reco_emu + muP*reco_dircosneu); reco_qe_q2 = abs(-2.0*reco_qe_enu*(reco_emu-muP*reco_dircosneu)+muM*muM); med_qe_pid=qeid.AltCalcQEID(this,event,reco_enu,reco_W2); med_qe_pass=qeid.PassMEDQECut(reco_enu,med_qe_pid); if(nsid.ChooseWeightFile(det,current_beam)){ if(!nsid.GetPID(this,event,det,niki_cc_pid)) niki_cc_pid=-999; } else niki_cc_pid=-999; // original code made detector dependent fiducial volume cut here // but, it's not necessary since PID has already been // calculated... // must reevaluate efficiency or make same cut downstream if(det==Detector::kFar && !file_is_mc){ if(recover==""||recover.find("birch")<recover.size()|| recover.find("Birch")<recover.size()|| recover.find("BIRCH")<recover.size()|| recover.find("R1.18")<recover.size()){ dpid.SetPHCorrection(1.018); } if(recover.find("cedar")<recover.size()|| recover.find("Cedar")<recover.size()|| recover.find("CEDAR")<recover.size()|| recover.find("R1.24")<recover.size()){ dpid.SetPHCorrection(1.0); } } if(dpid.ChoosePDFs(det,current_beam,recover,mcver)){ dave_cc_pid = dpid.CalcPID(this,event,0); } // if(abid.PassCuts(ntpEvent,strecord)){ andy_cc_pid=abid.CalcPID(ntpEvent,strecord); abid_costheta =abid.GetRecoCosTheta(ntpEvent,strecord); abid_eventenergy =abid.GetRecoE(ntpEvent,strecord); abid_trackcharge =abid.GetTrackEMCharge(ntpEvent,strecord); abid_trackenergy =abid.GetTrackPlanes(ntpEvent,strecord); abid_tracklikeplanes =abid.GetTrackLikePlanes(ntpEvent,strecord); abid_trackphfrac =abid.GetTrackPHfrac(ntpEvent,strecord); abid_trackphmean =abid.GetTrackPHmean(ntpEvent,strecord); abid_trackqpsigmaqp =abid.GetTrackQPsigmaQP(ntpEvent,strecord); abid_y =abid.GetRecoY(ntpEvent,strecord); // } TObject *recobj=0; if(record!=0){ recobj=record; } else{ recobj=strecord; } LoadTrack(track_index); //figure out planes in trk > 1.5 pe float temp_ph; trk_nplanes=NPlanesInObj(recobj,ntpTrack,1.5,temp_ph); //figure out strips in trk > 1.5 pe trk_nstrips=NStripsInObj(recobj,ntpTrack,1.5,temp_ph); } else { trklength = 0; trkmomrange = 0; trkqp = 0; trkeqp = 0; trkphfrac = 0; trkphplane = 0; is_trk_fid = -1; } //this used to load the slice associated with the track //but we really want the slice associated witht the event //tv 10-19-05 LoadSlice(ntpEvent->slc); // cout<<"*********SNARL="<<snarl<<" "<<event // <<"*********************************"<<endl; TObject *recobj=0; if(record!=0){ recobj=record; } else{ recobj=strecord; } float temp_ph; //figure out planes in slice > 1.5 pe slc_nplanes=NPlanesInObj(recobj,ntpSlice,1.5,temp_ph); //figure out planes in evt > 1.5 pe evt_nplanes=NPlanesInObj(recobj,ntpEvent,1.5,temp_ph); //figure out strips in slice > 1.5 pe slc_nstrips=NStripsInObj(recobj,ntpSlice,1.5,temp_ph); //figure out strips in evt > 1.5 pe evt_nstrips=NStripsInObj(recobj,ntpEvent,1.5,temp_ph); if(shower_index!=-1 && reco_eshw>0){ // case of no vertex shower LoadShower(shower_index); shwlength=ntpShower->plane.n; shw_nplanes_raw=NPlanesInObj(recobj,ntpShower,0.0,shw_ph_raw); shw_nplanes_cut=NPlanesInObj(recobj,ntpShower,1.5,shw_ph_cut); shwend=ntpShower->plane.end; shwbeg=ntpShower->plane.beg; reco_shw_vtxx=ntpShower->vtx.x; reco_shw_vtxy=ntpShower->vtx.y; reco_shw_vtxz=ntpShower->vtx.z; int ncl = ntpShower->ncluster; //Float_t zvtx //Float_t tposvtx //Float_t slope //Float_t avgdev /* static TH1F hswu("hswu","",400,4,-4); static TH1F hswu_wtd("hswu_wtd","",400,4,-4); static TH1F hswv("hswv","",400,4,-4); static TH1F hswv_wtd("hswv_wtd","",400,4,-4); hswu.Reset(); hswu_wtd.Reset(); hswv.Reset(); hswv_wtd.Reset(); */ Moments mom_swu; Moments mom_swu_wtd; Moments mom_swv; Moments mom_swv_wtd; etu=etv=elu=elv=0.0; swu=swv=wswu=wswv=0.0; // const float ss_cut=0.150;//GeV const float ss_cut=0.0;//GeV float totsigssu=0.0; // em and hadr subshowers float totsigssv=0.0; // em and hadr subshowers for(int ii=0; ii<ncl; ii++){ if(LoadCluster(ntpShower->clu[ii])){ if(ntpCluster->ph.gev >ss_cut){ nss[ntpCluster->id]+=1; ess[ntpCluster->id]+=ntpCluster->ph.gev; ntotss++; if(ntpCluster->id <2){ // slope = tpos/zpos? float sinang = sin(atan(ntpCluster->slope)); if(ntpCluster->planeview==PlaneView::kU){ totsigssu+=ntpCluster->ph.gev; //hswu.Fill(ntpCluster->tposvtx); mom_swu.AddPoint(ntpCluster->tposvtx); mom_swu_wtd.AddPoint(ntpCluster->tposvtx, ntpCluster->ph.gev); etu += sinang*ntpCluster->ph.gev; elu += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev; } else if(ntpCluster->planeview==PlaneView::kV){ totsigssv+=ntpCluster->ph.gev; //hswv.Fill(ntpCluster->tposvtx); mom_swv.AddPoint(ntpCluster->tposvtx); mom_swv_wtd.AddPoint(ntpCluster->tposvtx, ntpCluster->ph.gev); etv += sinang*ntpCluster->ph.gev; elv += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev; } } } } } /* for(int ii=0; ii<ncl; ii++){ if(LoadCluster(ntpShower->clu[ii])){ if(ntpCluster->ph.gev >ss_cut){ if(ntpCluster->id <2){ if(ntpCluster->planeview==PlaneView::kU){ totsigssu+=ntpCluster->ph.gev; hswu_wtd.Fill(ntpCluster->tposvtx); } else if(ntpCluster->planeview==PlaneView::kV){ totsigssv+=ntpCluster->ph.gev; hswv_wtd.Fill(ntpCluster->tposvtx); } } } } } */ /* swv=hswv.GetRMS(); swu=hswu.GetRMS(); wswv=hswv_wtd.GetRMS(); wswu=hswu_wtd.GetRMS(); */ swv=mom_swv.GetRMS(); swu=mom_swu.GetRMS(); wswv=mom_swv_wtd.GetRMS(); wswu=mom_swu_wtd.GetRMS(); } /* if(ntrack==1 && reco_enu>0.0 && reco_enu<140 && is_pitt_fid && pitt_trk_qual && (is_trk_fid!=0)) pass=1; */ if(ntrack>0 && trk_fit_pass==1 && is_fid==1) pass=1; } // if(reco_enu>0) } // if(PassBasicCuts()) // for both near in time and near in space // index of nearest, dist of nearest, ph of nearest // near_vtx_space_idx, near_vtx_space_dist, near_vtx_space_ph // nvsi, nvsd, nvsp // near_vtx_time_idx, near_vtx_time_dist, near_vtx_time_ph // nvti, nvtd, nvtp // for the nearest event, the index of the event nearest it // ph of event nearest it, distance of event nearest it // maybe use for event rehashing // back_vtx_space_idx, back_vtx_space_dist, back_vtx_space_ph // bvsi, bvsd, bvsp // back_vtx_time_idx, back_vtx_time_dist, back_vtx_time_ph // bvti, bvtd, bvtp nby.ClosestSpatial(i,nvsi,nvsd,nvst,nvsp,nvsevt); nby.ClosestSpatial(nvsi,bvsi,bvsd,bvst,bvsp,bvsevt); nby.ClosestTemporal(i,nvti,nvtd,nvtt,nvtp,nvtevt); nby.ClosestTemporal(nvti,bvti,bvtd,bvtt,bvtp,bvtevt); GetEvtTimeChar(evttimemin,evttimemax,evttimeln); GetEvtTimeCharPHC(evttimeminphc,evttimemaxphc,evttimelnphc); EarlyActivity(i,earlywph,earlywphpw,earlywphsphere, ph1000,ph500,ph100,lateph1000,lateph500, lateph100,lateph50); evtph = ntpEvent->ph.sigcor; DupRecoStrips(i,nduprs,dupphrs); //figure out batch number //note this will only match up with the spillmon info for //sure when there are 6 batches //if there are only 5 batches, we don't know if there //were only 5 batches in the MI, or if the first of those //batches went elsewhere. batchnumber=FindBatchNumber(evttimemin); tree->Fill(); } // loop over events } // loop over entries // delete nuparent; // should I do this? save.cd(); pottree->Fill(); pottree->Write(); tree->Write(); save.Write(); save.Close(); }

void MadTVAnalysis::DupRecoStrips (  int  evidx, int &  nduprs, float &  dupphrs )  [protected]

Definition at line 3156 of file MadTVAnalysis.cxx. References MadBase::LoadEvent(), MadBase::LoadStrip(), NtpSREventSummary::nevent, NtpSREvent::nstrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::sigcor, and NtpSREvent::stp. Referenced by CreatePAN(). { //zero the variables nduprs=0; dupphrs=0.; if(!LoadEvent(evidx)) return; //no event found map<int,float> sindex; for(int i=0;i<ntpEvent->nstrip;i++){ int si = ntpEvent->stp[i]; if(!LoadStrip(si)) continue; float ph = ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor; sindex[ntpEvent->stp[i]]=ph; } for(int i=0;i<eventSummary->nevent;i++){ if(i==evidx){ continue; } LoadEvent(i); //load a new event for(int j=0;j<ntpEvent->nstrip;j++){ map<int,float>::iterator it(sindex.find(ntpEvent->stp[j])); if(it!=sindex.end()){ nduprs++; dupphrs+=it->second; } } } //reload the current event LoadEvent(evidx); }

void MadTVAnalysis::EarlyActivity (  int  evidx, float &  earlywph, float &  earlywphpw, float &  earlywphsphere, float &  ph1000, float &  ph500, float &  ph100, float &  lateph1000, float &  lateph500, float &  lateph100, float &  lateph50 )  [protected]

Definition at line 2800 of file MadTVAnalysis.cxx. References fPMTTimeConstant, MadBase::LoadEvent(), MadBase::LoadStrip(), NtpSREventSummary::nstrip, NtpSREvent::nstrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRVertex::plane, NtpSRStrip::plane, NtpSRStrip::planeview, NtpSRStrip::pmtindex0, NtpSRStrip::pmtindex1, s(), NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRStrip::time0, NtpSRStrip::time1, NtpSRStrip::tpos, NtpSREventSummary::trigtime, NtpSRVertex::u, NtpSRVertex::v, NtpSREvent::vtx, NtpSRStrip::z, and NtpSRVertex::z. Referenced by CreatePAN(). { //zero the variables earlywph=0.; earlywphpw=0.; earlywphsphere=0.; ph1000=0.; ph500=0.; ph100=0.; lateph1000=0.; lateph500=0.; lateph100=0.; lateph50=0.; if(!LoadEvent(evidx)) return; //no event found //find the time of the earliest strip in the event double earliestEventTime = 1.e20; double latestEventTime = 0.; map<int,int> eventPlanes; double trigtime = eventSummary->trigtime; //loop over strips in the event //find the earliest strip time //make a map of the pmt's hit in this event for(int s=0;s<ntpEvent->nstrip;s++){ int stripindex = ntpEvent->stp[s]; if(!LoadStrip(stripindex)) continue; //calculate ph weighted strip time (over the two ends, //should be equivalent to ph1 in ND float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+ ntpStrip->ph0.sigcor*ntpStrip->time0)/ (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor)); double time = 1.e9*(phwt-trigtime); if(time<earliestEventTime){ earliestEventTime = time; } if(time>latestEventTime){ latestEventTime=time; } if(ntpStrip->ph1.sigcor>0){ if(eventPlanes.find(ntpStrip->pmtindex1)==eventPlanes.end()){ eventPlanes[ntpStrip->pmtindex1] = 1; } } if(ntpStrip->ph0.sigcor>0){ if(eventPlanes.find(ntpStrip->pmtindex0)==eventPlanes.end()){ eventPlanes[ntpStrip->pmtindex1] = 1; } } } //now find the weigthed sum of ADC for the early strips - make sure the early //strips come from the same pmts as the event strips. for(unsigned int s=0;s<eventSummary->nstrip;s++){ if(!LoadStrip(s)) continue; //work in ns float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+ ntpStrip->ph0.sigcor*ntpStrip->time0)/ (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor)); double newtime = 1.e9*(phwt-trigtime); double dt=(earliestEventTime-newtime); double dtlate=newtime-latestEventTime; float d=1000.; if((int)ntpStrip->planeview==PlaneView::kU){ d=sqrt((pow(ntpEvent->vtx.u-ntpStrip->tpos,2)+ pow(ntpEvent->vtx.z-ntpStrip->z,2))); } else{ d=sqrt((pow(ntpEvent->vtx.v-ntpStrip->tpos,2)+ pow(ntpEvent->vtx.z-ntpStrip->z,2))); } if(d<fEarlyPlaneSphere&&dtlate>0.){ if(dtlate<1000.){ lateph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } if(dtlate<500.){ lateph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } if(dtlate<100.){ lateph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } if(dtlate<50.){ lateph50+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } } if(d<fEarlyPlaneSphere&&dt>0.){ if(dt<1000.){ ph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } if(dt<500.){ ph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } if(dt<100.){ ph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; } } if(dt>0.&&dt<1000.*fEarlyActivityWindowSize){ double eadc = ((ntpStrip->ph1.sigcor+ntpStrip->ph1.sigcor)* TMath::Exp(-1.*dt/fPMTTimeConstant)); //if this strip has activity in the same pmt as our event, //sum up the early pulse height if(eventPlanes.find(ntpStrip->pmtindex1)!=eventPlanes.end()){ earlywph += eadc; } //if this strip is within some radius of event vertex, //add up the early pulse height if(d<fEarlyPlaneSphere){ earlywphsphere+=eadc; } //if this strip is within the plane window of event vertex, //add up the early pulseheight if(ntpStrip->plane>=ntpEvent->vtx.plane&& ntpStrip->plane<=ntpEvent->vtx.plane+fNPlaneEarlyAct){ earlywphpw+=eadc; } } } }

int MadTVAnalysis::FarTrkContained (  Float_t  x, Float_t  y, Float_t  z )  [static]

Definition at line 3116 of file MadTVAnalysis.cxx. { int is_cont=0; //arguments should correspond to track end //returns: //1 if track is partially contained //2 if track is fully contained const Float_t r = sqrt(x*x+y*y); if(r<=3.5&&r>.5&&((z>=0.5&&z<=14.5)||(z>=16.5&&z<=29.4))){//contained is_cont=2; } else{//not contained is_cont=1; } return is_cont; }

Int_t MadTVAnalysis::FDRCBoundary (   )  [protected]

Definition at line 3189 of file MadTVAnalysis.cxx. References NtpSRPlane::beg, NtpSRPlane::end, NtpSREventSummary::nshower, NtpSREventSummary::ph, NtpSREventSummary::plane, and NtpSRPulseHeight::raw. Referenced by CreatePAN(). { Int_t 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++; } return litag; }

void MadTVAnalysis::FillMCHists (  TFile *  fout  )

Definition at line 2931 of file MadTVAnalysis.cxx. References NtpMCTruth::iaction, NtpMCTruth::inu, NtpMCTruth::iresonance, MadBase::LoadTruth(), NDRadialFidVol(), NtpMCSummary::nmc, NtpMCTruth::p4neu, PittInFidVol(), NtpMCTruth::vtxx, NtpMCTruth::vtxy, and NtpMCTruth::vtxz. Referenced by CreatePAN(). { static bool first=true; static TH1* h_numu_pf_tot = new TH1F("h_numu_pf_tot","CC; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_pf_qe = new TH1F("h_numu_pf_qe","QE; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_pf_res = new TH1F("h_numu_pf_res","RES; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_pf_dis = new TH1F("h_numu_pf_dis","DIS; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_pf_nc = new TH1F("h_numu_pf_nc","NC; true neutrino energy (GeV)", 240,0.0,120.0); if(first){ h_numu_pf_tot->SetDirectory(fout); h_numu_pf_qe->SetDirectory(fout); h_numu_pf_res->SetDirectory(fout); h_numu_pf_dis->SetDirectory(fout); h_numu_pf_nc->SetDirectory(fout); } static TH1* h_numubar_pf_tot = new TH1F("h_numubar_pf_tot","CC; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_pf_qe = new TH1F("h_numubar_pf_qe","QE; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_pf_res = new TH1F("h_numubar_pf_res","RES; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_pf_dis = new TH1F("h_numubar_pf_dis","DIS; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_pf_nc = new TH1F("h_numubar_pf_nc","NC; true neutrino energy (GeV)", 240,0.0,120.0); if(first){ h_numubar_pf_tot->SetDirectory(fout); h_numubar_pf_qe->SetDirectory(fout); h_numubar_pf_res->SetDirectory(fout); h_numubar_pf_dis->SetDirectory(fout); h_numubar_pf_nc->SetDirectory(fout); } static TH1* h_numu_1m_tot = new TH1F("h_numu_1m_tot","CC; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_1m_qe = new TH1F("h_numu_1m_qe","QE; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_1m_res = new TH1F("h_numu_1m_res","RES; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_1m_dis = new TH1F("h_numu_1m_dis","DIS; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numu_1m_nc = new TH1F("h_numu_1m_nc","NC; true neutrino energy (GeV)", 240,0.0,120.0); if(first){ h_numu_1m_tot->SetDirectory(fout); h_numu_1m_qe->SetDirectory(fout); h_numu_1m_res->SetDirectory(fout); h_numu_1m_dis->SetDirectory(fout); h_numu_1m_nc->SetDirectory(fout); } static TH1* h_numubar_1m_tot = new TH1F("h_numubar_1m_tot","CC; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_1m_qe = new TH1F("h_numubar_1m_qe","QE; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_1m_res = new TH1F("h_numubar_1m_res","RES; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_1m_dis = new TH1F("h_numubar_1m_dis","DIS; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_numubar_1m_nc = new TH1F("h_numubar_1m_nc","NC; true neutrino energy (GeV)", 240,0.0,120.0); if(first){ h_numubar_1m_tot->SetDirectory(fout); h_numubar_1m_qe->SetDirectory(fout); h_numubar_1m_res->SetDirectory(fout); h_numubar_1m_dis->SetDirectory(fout); h_numubar_1m_nc->SetDirectory(fout); } static TH1* h_other_pf = new TH1F("h_other_pf", "other; true neutrino energy (GeV)", 240,0.0,120.0); static TH1* h_other_1m = new TH1F("h_other_1m", "other; true neutrino energy (GeV)", 240,0.0,120.0); if(first){ h_other_pf->SetDirectory(fout); h_other_1m->SetDirectory(fout); } if(!mcHeader) { first=true; return; } for(Int_t i=0; i<mcHeader->nmc; i++){ if(!LoadTruth(i)) continue; // no ntpmctruth const Float_t k = 1.0/sqrt(2.0); Float_t x = ntpTruth->vtxx; Float_t y = ntpTruth->vtxy; Float_t z = ntpTruth->vtxz; Float_t u = k*(y+x); Float_t v = k*(y-x); Bool_t inpf = PittInFidVol(x,y,z,u,v); Int_t ndr = NDRadialFidVol(x,y,z); Bool_t inr =kFALSE; if((ndr&0x8)>0) inr=kTRUE; Int_t ccnc = ntpTruth->iaction; Int_t process = ntpTruth->iresonance; Int_t inu = ntpTruth->inu; Float_t E = ntpTruth->p4neu[3]; if(inu==14){ // neutrinos if(ccnc==0){ if(inpf) h_numu_pf_nc->Fill(E); if(inr) h_numu_1m_nc->Fill(E); } else if(ccnc==1){ if(inpf) h_numu_pf_tot->Fill(E); if(inr) h_numu_1m_tot->Fill(E); switch(process){ case 1001: if(inpf) h_numu_pf_qe->Fill(E); if(inr) h_numu_1m_qe->Fill(E); break; case 1002: if(inpf) h_numu_pf_res->Fill(E); if(inr) h_numu_1m_res->Fill(E); break; case 1003: if(inpf) h_numu_pf_dis->Fill(E); if(inr) h_numu_1m_dis->Fill(E); break; default: break; } } } else if(inu==-14){ if(ccnc==0){ if(inpf) h_numubar_pf_nc->Fill(E); if(inr) h_numubar_1m_nc->Fill(E); } else if(ccnc==1){ if(inpf) h_numubar_pf_tot->Fill(E); if(inr) h_numubar_1m_tot->Fill(E); switch(process){ case 1001: if(inpf) h_numubar_pf_qe->Fill(E); if(inr) h_numubar_1m_qe->Fill(E); break; case 1002: if(inpf) h_numubar_pf_res->Fill(E); if(inr) h_numubar_1m_res->Fill(E); break; case 1003: if(inpf) h_numubar_pf_dis->Fill(E); if(inr) h_numubar_1m_dis->Fill(E); break; default: break; } } } else{ if(inpf) h_other_pf->Fill(E); if(inr) h_other_1m->Fill(E); } } first=false; }

int MadTVAnalysis::FindBatchNumber (  double  mintimediff  )  [protected]

Definition at line 3245 of file MadTVAnalysis.cxx. Referenced by CreatePAN(). { //cuts from peter shannahan //Batch 0: 1.7e-6 to 3.0e-6 //Batch 1: 3.3e-6 to 4.5e-6 //Batch 2: 4.9e-6 to 6.2e-6 //Batch 3: 6.5e-6 to 7.8e-6 //Batch 4: 8.1e-6 to 9.4e-6 //Batch 5: 9.7e-6 to 11.0e-6 //mintime already has triggertime subtracted double tdiff = mintime*1.e6; if(tdiff>1.7&&tdiff<3.0) return 0; if(tdiff>3.3&&tdiff<4.5) return 1; if(tdiff>4.9&&tdiff<6.2) return 2; if(tdiff>6.5&&tdiff<7.8) return 3; if(tdiff>8.1&&tdiff<9.4) return 4; if(tdiff>9.7&&tdiff<11.0) return 5; else return -1; }

void MadTVAnalysis::ForceMCBeam (  const BeamType::BeamType_t &  beam  )  [inline]

Definition at line 63 of file MadTVAnalysis.h. References fMCBeam. {fMCBeam=beam;}

Registry& MadTVAnalysis::GetBECRegistry (   )  [inline, protected]

Definition at line 113 of file MadTVAnalysis.h. { return fBECConfig; }

void MadTVAnalysis::GetEvtTimeChar (  double &  timemin, double &  timemax, double &  timeln )  [protected]

Definition at line 2627 of file MadTVAnalysis.cxx. References VldContext::GetDetector(), VldTimeStamp::GetNanoSec(), VldContext::GetTimeStamp(), RecHeader::GetVldContext(), MadBase::LoadStrip(), NtpSREvent::nstrip, NtpSREvent::stp, striptime, NtpSRStrip::time0, and NtpSRStrip::time1. Referenced by CreatePAN(). { //cut and pasted, with minor changes from MadDpAnalysis.cxx // double trgtime=eventSummary->trigtime; double trgtime=ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1.e9; double timemax=0.; double timemin=1.e10; // Find max min time of events by loading strips for ND if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); if(!LoadStrip(stp_index)) continue; Double_t striptime=ntpStrip->time1-trgtime; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } } if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); if(!LoadStrip(stp_index)) continue; Double_t striptime1=ntpStrip->time1-trgtime; Double_t striptime0=ntpStrip->time0-trgtime; // Double_t striptime1=ntpStrip->time1; // Double_t striptime0=ntpStrip->time0; Double_t striptime=(striptime0+striptime1)/2.; // if(striptime1>0 && striptime0<0) striptime=striptime1; // if(striptime0>0 && striptime1<0) striptime=striptime0; // if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; if(striptime1>-1 && striptime0<-1) striptime=striptime1; if(striptime0>-1 && striptime1<-1) striptime=striptime0; if(striptime0>-1 && striptime1>-1) striptime=(striptime0+striptime1)/2.; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } } evttimemax=timemax; evttimemin=timemin; evttimeln=timemax-timemin; }

void MadTVAnalysis::GetEvtTimeCharPHC (  double &  timemin, double &  timemax, double &  timeln )  [protected]

Definition at line 2673 of file MadTVAnalysis.cxx. References VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadStrip(), NtpSREvent::nstrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::sigcor, NtpSREvent::stp, striptime, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSREventSummary::trigtime. Referenced by CreatePAN(). { //same as GetEvtTimeChar but only looking at hits with //pulse height > 200 sig cor (~2pe) in near det //pulse height > 160 sig cor in far det double trgtime=eventSummary->trigtime; double timemax=0.; double timemin=1.e10; // Find max min time of events by loading strips for ND if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); if(!LoadStrip(stp_index)) continue; Double_t striptime=ntpStrip->time1-trgtime; if(ntpStrip->ph1.sigcor>200){ if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } } } if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ Int_t stp_index=((ntpEvent->stp)[evss]); if(!LoadStrip(stp_index)) continue; Double_t striptime1=ntpStrip->time1-trgtime; Double_t striptime0=ntpStrip->time0-trgtime; Double_t striptime=(striptime0+striptime1)/2.; if(striptime1>0 && striptime0<0) striptime=striptime1; if(striptime0>0 && striptime1<0) striptime=striptime0; if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; if(ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor>160){ if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } } } evttimemax=timemax; evttimemin=timemin; evttimeln=timemax-timemin; }

const BeamType::BeamType_t& MadTVAnalysis::GetForcedMCBeam (   )  [inline]

Definition at line 64 of file MadTVAnalysis.h. {return fMCBeam;}

Float_t MadTVAnalysis::GetNDCoilCurrent (  const VldContext &  vc  )  [protected]

Definition at line 3235 of file MadTVAnalysis.cxx. References Dcs_Mag_Near::GetCurrent(), DbiResultPtr< T >::GetNumRows(), and DbiResultPtr< T >::GetRow(). { Float_t result = 0.0; DbiResultPtr<Dcs_Mag_Near> ptr(vc); if(ptr.GetNumRows()){ const Dcs_Mag_Near* row = ptr.GetRow(0); result = row->GetCurrent(); } return result; }

bool MadTVAnalysis::InFidVol (  float  vx, float  vy, float  vz )  [protected]

Definition at line 2068 of file MadTVAnalysis.cxx. References VldContext::GetDetector(), RecHeader::GetVldContext(), kXcenter, and kYcenter. { // assumes event has been loaded // check that vertex is in fiducial volume bool is_fid=false; // FD: from doc-1351-v2 : there is no coil cut! if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ //fiducial criteria on vtx for far detector is_fid = true; if((vz<0.5 || vz>28.0) || //ends (vz<16.2 && vz>14.3) || //between SMs ((vx*vx)+(vy*vy))>14.0) is_fid = false; } else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ //fiducial criteria on vtx for near detector // float beamzerox = 1.4885; // float beamzeroy = 0.1397-reco_vtxz*TMath::Tan(.058); // float r=sqrt(pow((reco_vtxx-beamzerox),2)+pow((reco_vtxy-beamzeroy),2)); is_fid = false; // MAK: Sept 19, 2005: Modified to account for beam slope /* double radius = pow(vx-kXcenter,2) + pow(vy-(kYcenter -vz*TMath::Tan(0.058)),2); radius=sqrt(radius); */ // MAK: removed beam slope. Dave/Niki don't do it. //double radius = sqrt(vx*vx + vy*vy); double radius = pow(vx-kXcenter,2) + pow(vy-kYcenter,2); radius=sqrt(radius); // restrict vertices to target region and 1m radius if(vz>=1.0 && vz<=5.0 && radius<=1.0) is_fid = true; } return is_fid; }

bool MadTVAnalysis::InFidVol (   )  [protected, virtual]

Definition at line 2063 of file MadTVAnalysis.cxx. References NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by CreatePAN(), and InFidVol(). { return InFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,ntpEvent->vtx.z); }

bool MadTVAnalysis::InFidVol (  Int_t  evidx  )  [protected, virtual]

Definition at line 2058 of file MadTVAnalysis.cxx. References InFidVol(), and MadBase::LoadEvent(). { if(!LoadEvent(evidx)) return false; //no event found return InFidVol(); }

Int_t MadTVAnalysis::InPartialRegion (  UShort_t  plane, UShort_t  strip )  [static]

Definition at line 2346 of file MadTVAnalysis.cxx. Referenced by SigInOut(). { // figure out if this (plane,strip) corresponds to something in // the partially instrumented region // // this region is defined as: // v planes: (strip<=4 || strip>=67) // partial u: (strip==0 || strip=63) // full u: (strip<=26 || strip>=88) // // if so, return 1 // if not, return -1 // if error, return 0 // make a lookup ptype to hold the type of each plane // 1 = v partial 2 = u partial // 3 = v full 4 = u full // 0 = uninstrumented static bool first=true; static UShort_t ptype[282]; if(first){ ptype[0]=0; for(int i=1; i<=281; i++){ if(i%2==0) ptype[i]=1; // a v plane else ptype[i]=2; // a u plane if((i-1)%5 == 0) ptype[i]+=2; // fully instrumented else if(i>120) ptype[i]=0; // not instrumented } first=false; } if(plane>281){ // std::cerr<<"InPartialRegion passed plane = "<<plane<<std::endl; return 0; } UShort_t pt = ptype[plane]; Int_t result; switch(pt){ case 1: case 3: if(strip<=4 || strip>=67) result=1; else result = -1; break; case 2: if(strip==0 || strip == 63) result=1; else result = -1; break; case 4: if(strip<=26 || strip>=88) result=1; else result = -1; break; case 0: default: result=0; break; } return result; }

float MadTVAnalysis::MuonShowerEnergy (  const NtpSREvent *  evt, const NtpSRTrack *  trk )  [protected]

Definition at line 2720 of file MadTVAnalysis.cxx. References NtpSRShowerPulseHeight::EMgev, MadBase::LoadShower(), NtpSREvent::nshower, NtpSRTrack::nstrip, NtpSREvent::shw, NtpSRShower::shwph, NtpSRTrack::stpu, NtpSRTrack::stpv, NtpSRTrack::stpz, NtpSRVertex::t, NtpSRVertex::u, NtpSRVertex::v, NtpSRTrack::vtx, NtpSREvent::vtx, NtpSRShower::vtx, and NtpSRVertex::z. Referenced by CreatePAN(). { // Purpose: try to sum up brems along the track static TH2* hr = new TH2F("h_mushw_rph", "; track - shw r dist (m); shower energy (GeV)", 10,0,0.2,40,0,2); static TH2* hz = new TH2F("h_mushw_zph", "; track - shw z dist (m); shower energy (GeV)", 20,-0.2,0.2,40,0,2); static TH2* ht = new TH2F("h_mushw_tph", "; track - shw t dist (m); shower energy (GeV)", 600,-1e-6,10e-6,40,0,2); if(evt==0 || trk==0) return 0; float result=0; NtpSRShower* save_shwr = ntpShower; // loop through all showers in event const float max_vtx_dist=14*0.06; // 2 interaction lengths for(int i=0; i<evt->nshower; i++){ LoadShower(evt->shw[i]); if(!ntpShower) continue; // calculate distance to vertex // don't want to use showers near the vertex float dist_vtx = pow(ntpShower->vtx.z - ntpEvent->vtx.z,2) + pow(ntpShower->vtx.u - ntpEvent->vtx.u,2) + pow(ntpShower->vtx.v - ntpEvent->vtx.v,2); dist_vtx=sqrt(dist_vtx); if(dist_vtx<max_vtx_dist) continue; // loop through track's strip array // calculate distance from strip to shower in time and space // save spatial and time distance of closest approach to track // if shower is close enough add the shower energy to running sum. float min_dist=1e6; float min_r=1e6; float min_z=1e6; float min_t=1e6; float min_ph=0; float add_ph=0; for(int j=0; j<trk->nstrip; j++){ float dist_trk= pow(trk->stpz[j]-ntpShower->vtx.z,2) + pow(trk->stpu[j]-ntpShower->vtx.u,2) + pow(trk->stpv[j]-ntpShower->vtx.v,2); dist_trk=sqrt(dist_trk); float tdist = (trk->vtx.t + (trk->stpz[j]-trk->vtx.z)/3e8) - ntpShower->vtx.t; // = trk->stpt0[j]*trk->stpph0[j]+trk->stp1[j]*trk->stpph1[j]; // if( (trk->stpph0[j]+trk->stpph1[j]) > 0){ // tdist/=(trk->stpph0[j]+trk->stpph1[j]); // } if(fabs(dist_trk)<fabs(min_dist)) { min_ph=ntpShower->shwph.EMgev; min_z = fabs(ntpShower->vtx.z-trk->stpz[j]); min_r = sqrt( pow(ntpShower->vtx.u-trk->stpu[j],2) + pow(ntpShower->vtx.v-trk->stpv[j],2)); min_t = tdist; } const float max_trk_dist=1.76*3.0; // 3 radiation lengths const float max_trk_tdist=40e-9; if(fabs(dist_trk)<max_trk_dist && fabs(tdist)<max_trk_tdist){ // shower correllated with track add_ph=min_ph; } } result+=min_ph; // // fill histograms hr->Fill(min_r,min_ph); hz->Fill(min_z,min_ph); ht->Fill(min_t,min_ph); }// end loop over showers // reset shower pointer ntpShower=save_shwr; // all done return result; }

Int_t MadTVAnalysis::NDRadialFidVol (  Float_t  x, Float_t  y, Float_t  z )  [static]

Definition at line 2262 of file MadTVAnalysis.cxx. References kTargEndZ, kVetoEndZ, kXcenter, and kYcenter. Referenced by CreatePAN(), and FillMCHists(). { Float_t rings[5]={0.1,0.25,0.5,1.0,1.5}; // in meters Int_t result=0; for (unsigned int i=0; i<5; i++){ Float_t radius = pow(x-kXcenter,2) + pow(y-(kYcenter - z*TMath::Tan(0.058)),2); radius=sqrt(radius); // restrict vertices to target region and 1m radius if(z>=kVetoEndZ && z<=kTargEndZ && radius<=rings[i]) result|=(1<<i); } return result; }

Int_t MadTVAnalysis::NPlanesInObj (  TObject *  rec, TObject *  obj, float  phcut, float &  sumph )  [protected]

Definition at line 2464 of file MadTVAnalysis.cxx. References NtpStRecord::evthdr, NtpSRRecord::evthdr, NtpSREventSummary::nstrip, NtpSRShower::nstrip, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRRecord::stp, NtpSRShower::stp, NtpSREvent::stp, NtpSRTrack::stp, and NtpSRSlice::stp. Referenced by CreatePAN(). { sumph=0.0; std::set<UShort_t> planes; NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec); NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec); Int_t *stripindex=0; Int_t nstrip=0; NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj); NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj); NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj); NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj); if(slc!=0){ stripindex=slc->stp; nstrip=slc->nstrip; // cout<<"Counting planes in slice"<<endl; } else if(trk!=0){ stripindex=trk->stp; nstrip=trk->nstrip; // cout<<"Counting planes in trk"<<endl; } else if(evt!=0){ stripindex=evt->stp; nstrip=evt->nstrip; // cout<<"Counting planes in event"<<endl; } else if(shw!=0){ stripindex=shw->stp; nstrip=shw->nstrip; // cout<<"Counting planes in event"<<endl; } TClonesArray *striplist=0; int TOTSTRIPS=0; if(srrec!=0){ striplist = srrec->stp; TOTSTRIPS=srrec->evthdr.nstrip; // cout<<"Found a sr"<<endl; } else{ striplist = strec->stp; TOTSTRIPS=strec->evthdr.nstrip; // cout<<"found a st"<<endl; } if(striplist!=0){ // cout<<"Looping over strip list, nstrip="<<nstrip<<endl; for(int i=0;i<nstrip;i++){ int sindex = stripindex[i]; if(sindex<0) continue; if(sindex<TOTSTRIPS){ NtpSRStrip *strip = dynamic_cast<NtpSRStrip *> ((*striplist)[sindex]); // cout<<"got strip "<<sindex<<" plane "<<strip->plane // <<" strip "<<strip->strip // <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl; if(strip->ph0.pe+strip->ph1.pe>phcut){ planes.insert(strip->plane); sumph+=strip->ph0.sigcor+strip->ph1.sigcor; } } } } else{ // cout<<"Getting number of planes in snarl"<<endl; for(int i=0;i<TOTSTRIPS;i++){ NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]); // cout<<"got strip "<<i<<" plane "<<strip->plane // <<" strip "<<strip->strip // <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl; if(strip->ph0.pe+strip->ph1.pe>phcut){ planes.insert(strip->plane); sumph+=strip->ph0.sigcor+strip->ph1.sigcor; } } } // cout<<"I count "<<planes.size()<<" planes > "<<phcut<<" pe"<<endl; return planes.size(); }

Int_t MadTVAnalysis::NStripsInObj (  TObject *  rec, TObject *  obj, float  phcut, float &  sumph )  [protected]

Definition at line 2551 of file MadTVAnalysis.cxx. References NtpStRecord::evthdr, NtpSRRecord::evthdr, NtpSREventSummary::nstrip, NtpSRShower::nstrip, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRRecord::stp, NtpSRShower::stp, NtpSREvent::stp, NtpSRTrack::stp, and NtpSRSlice::stp. Referenced by CreatePAN(). { sumph=0.0; std::vector<UShort_t> strips; NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec); NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec); Int_t *stripindex=0; Int_t nstrip=0; NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj); NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj); NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj); NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj); if(slc!=0){ stripindex=slc->stp; nstrip=slc->nstrip; } else if(trk!=0){ stripindex=trk->stp; nstrip=trk->nstrip; } else if(evt!=0){ stripindex=evt->stp; nstrip=evt->nstrip; } else if(shw!=0){ stripindex=shw->stp; nstrip=shw->nstrip; } TClonesArray *striplist=0; int TOTSTRIPS; if(srrec!=0){ striplist = srrec->stp; TOTSTRIPS=srrec->evthdr.nstrip; } else{ striplist = strec->stp; TOTSTRIPS=strec->evthdr.nstrip; } // cout<<"CHECKING nstrips="<<nstrip<<endl; if(striplist!=0){ for(int i=0;i<nstrip;i++){ int sindex=stripindex[i]; if(sindex<0) continue; NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[sindex]); if(strip->ph0.pe+strip->ph1.pe>phcut){ strips.push_back(strip->plane); sumph+=strip->ph0.sigcor+strip->ph1.sigcor; } } } else{ if(srrec!=0){ nstrip=srrec->evthdr.nstrip; } else{ nstrip=strec->evthdr.nstrip; } for(int i=0;i<nstrip;i++){ NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]); // cout<<"Checking "<<i<<" stripph: "<<strip->ph0.pe<<" "<<strip->ph1.pe<<endl; if(strip->ph0.pe+strip->ph1.pe>phcut){ strips.push_back(strip->plane); sumph+=strip->ph0.sigcor+strip->ph1.sigcor; } } } // cout<<"returning "<<strips.size()<<endl; return strips.size(); }

Bool_t MadTVAnalysis::PittInFidVol (  Float_t  x, Float_t  y, Float_t  z, Float_t  u, Float_t  v )  [static]

Definition at line 2278 of file MadTVAnalysis.cxx. Referenced by CreatePAN(), and FillMCHists(). { // Is the event vertex, (x,y,u,v,z) inside the fiducial volume // used by the pittsburgh group? if( !( z>0.6 && z<3.56) ) return kFALSE; if( !( u>0.3 && u<1.8) ) return kFALSE; if( !( v>-1.8 && v< -0.3) ) return kFALSE; if( x>=2.4) return kFALSE; const Float_t coil_cut=0.8*0.8; if( (x*x + y*y) <= coil_cut) return kFALSE; return kTRUE; }

Int_t MadTVAnalysis::PittTrkContained (  Float_t  x, Float_t  y, Float_t  z, Float_t  u, Float_t  v )  [static]

Definition at line 2294 of file MadTVAnalysis.cxx. Referenced by CreatePAN(). { // what is the containment status of this track // as defined by the pitt group // // takes the track end point (x,y,u,v,z) // returns: // 1 = track is fully contained in the upstream region // 2 = track is partially contained in the upstream region // 3 = track is fully contained in the downstream region // 4 = track is partially contained in the downstream region // // the rationale is that these categories have different momentum resolution // const Float_t radsq = x*x + y*y; int result=0; if(z<7.0) { // track in upstream region static const Float_t coil_cut=0.8*0.8; // was asking for u>3 && u<1.8 here before -- a bug if( (u>0.3 && u<1.8) && (v>-1.8 && v<-0.3) && (x<2.4) && (radsq>coil_cut) ) result=1; else result=2; } else{ // downstream region // uses an ellipse in x,y // major axis (x) = a = 1.7 m // minor axis (y) = b = 1.4 m // centered on x,y = (0.8,0) // (x-x0)^2/a^2 + (y-y0)^2/b^2 = 1 // should change this to 0.8 // static const Float_t coil_cut=0.5*0.5; // note differs from above //changed from 0.5 to 0.8 by TV on 7-29-05 static const Float_t coil_cut=0.8*0.8; static const Float_t x0=0.8; static const Float_t y0=0.0; static const Float_t a=1.7; static const Float_t b=1.4; const Float_t xsc = (x-x0)/a; // rescale ellipse to unit circle const Float_t ysc = (y-y0)/b; //MAK: cut on z was 16.0 till Aug 2, 2005 if( (sqrt(xsc*xsc + ysc*ysc)<1.0) && (radsq>coil_cut) && (z<15.6)) result = 3; else result = 4; } return result; }

Float_t MadTVAnalysis::RecoMuDCosNeuFD (  Int_t  itr = 0, Float_t *  vtx = 0 )  [protected, virtual]

Definition at line 2187 of file MadTVAnalysis.cxx. References NtpSRVertex::dcosy, NtpSRVertex::dcosz, MadBase::LoadTrack(), and NtpSRTrack::vtx. Referenced by CreatePAN(). { if(!LoadTrack(itr)) return 0.; Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam Float_t bl_y = sqrt(1. - bl_z*bl_z); Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y; return costhbl; }

Float_t MadTVAnalysis::RecoMuDCosNeuND (  Int_t  itr = 0, Float_t *  vtx = 0 )  [protected, virtual]

Definition at line 2197 of file MadTVAnalysis.cxx. References NtpSRVertex::dcosy, NtpSRVertex::dcosz, MadBase::LoadTrack(), and NtpSRTrack::vtx. Referenced by CreatePAN(). { if(!LoadTrack(itr)) return 0.; /* // simple correction based on the vertical position float vtxy=0; if(vtx) vtxy=vtx[1]; // in meters // cosine of the typical neutrino angle in the yz plane // calculated as py / sqrt( py^2 + pz^2) float nu_cos = -5.799E-2; if(vtxy>-2.0 && vtxy<2.0){ //prevents further nuttyness if vtxy is silly nu_cos += vtxy*1.23304E-3 + vtxy*vtxy*1.08212E-5 + vtxy*vtxy*vtxy*(-4.634E-5); } */ float nu_cos = -5.799E-2; float nu_sin = sqrt(1 -nu_cos*nu_cos); float cosz = ntpTrack->vtx.dcosz*nu_sin + ntpTrack->vtx.dcosy*nu_cos; return cosz; }

void MadTVAnalysis::SetBECFile (  const char *  f  )  [protected]

Definition at line 2458 of file MadTVAnalysis.cxx. References fBECConfig, and Registry::Set(). { if(f){ fBECConfig.Set("beam:flux_file",f); } }

void MadTVAnalysis::SetDataUseMCBeam (  bool  x  )  [inline]

Definition at line 65 of file MadTVAnalysis.h. References fDataUseMCBeam. {fDataUseMCBeam=x;}

void MadTVAnalysis::SetRustemConfig (  string  fname  )  [inline]

Definition at line 73 of file MadTVAnalysis.h. References rustemconfigfile. {rustemconfigfile=fname;}

void MadTVAnalysis::SetRustemFarPID (  string  fname  )  [inline]

Definition at line 72 of file MadTVAnalysis.h. References rustempidfile_far. {rustempidfile_far=fname;}

void MadTVAnalysis::SetRustemNearPID (  string  fname  )  [inline]

Definition at line 71 of file MadTVAnalysis.h. References rustempidfile_near. {rustempidfile_near=fname;}

void MadTVAnalysis::SigInOut (  Int_t  trkidx, Float_t &  sigfull, Float_t &  sigpart, Float_t &  trkfull, Float_t &  trkpart )  [protected]

Definition at line 2450 of file MadTVAnalysis.cxx. References MadBase::LoadTrack(), and SigInOut(). { // will try to load the track pointed to by trkidx if(trkidx==-1) ntpTrack=0; else if (!(LoadTrack(trkidx))) ntpTrack=0; SigInOut(sigfull,sigpart,trkfull,trkpart); }

void MadTVAnalysis::SigInOut (  Float_t &  sigfull, Float_t &  sigpart, Float_t &  trkfull, Float_t &  trkpart )  [protected]

Definition at line 2406 of file MadTVAnalysis.cxx. References InPartialRegion(), MadBase::LoadStrip(), NtpSRTrack::nstrip, NtpSREvent::nstrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSRTrack::stp, NtpSREvent::stp, and NtpSRStrip::strip. Referenced by CreatePAN(), and SigInOut(). { // compute the amount of signal in the partially instrumented region // and the amount in the fully instrumented region // // this region is defined as: // v planes: (strip<=4 || strip>=67) // partial u: (strip==0 || strip=63) // full u: (strip<=26 || strip>=88) sigfull=sigpart=0; // loop over all strips in the event // and sum the signals in the partial and full regions for(int i=0; i<ntpEvent->nstrip; i++){ if(!LoadStrip(ntpEvent->stp[i])) continue; Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip); if(pr==1){ sigpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; } else if(pr==-1){ sigfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; } } // loop over those strips in the primary track trkfull=trkpart=0; if(ntpTrack){ for(int i=0; i<ntpTrack->nstrip; i++){ if(!LoadStrip(ntpTrack->stp[i])) continue; Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip); if(pr==1){ trkpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; } else if(pr==-1){ trkfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor; } } } }

---

Member Data Documentation

MadAbID MadTVAnalysis::abid

Definition at line 80 of file MadTVAnalysis.h. Referenced by CreatePAN().

MadDpID MadTVAnalysis::dpid

Definition at line 78 of file MadTVAnalysis.h. Referenced by CreatePAN().

Registry MadTVAnalysis::fBECConfig [private]

Definition at line 167 of file MadTVAnalysis.h. Referenced by CreatePAN(), MadTVAnalysis(), and SetBECFile().

bool MadTVAnalysis::fDataUseMCBeam [private]

Definition at line 169 of file MadTVAnalysis.h. Referenced by CreatePAN(), MadTVAnalysis(), and SetDataUseMCBeam().

const Double_t MadTVAnalysis::fEarlyActivityWindowSize = 10. [static, private]

Definition at line 55 of file MadTVAnalysis.cxx.

const Double_t MadTVAnalysis::fEarlyPlaneSphere = 0.5 [static, private]

Definition at line 58 of file MadTVAnalysis.cxx.

BeamType::BeamType_t MadTVAnalysis::fMCBeam [private]

Definition at line 168 of file MadTVAnalysis.h. Referenced by ForceMCBeam(), and MadTVAnalysis().

const Int_t MadTVAnalysis::fNPlaneEarlyAct = 30 [static, private]

Definition at line 57 of file MadTVAnalysis.cxx.

const Double_t MadTVAnalysis::fPMTTimeConstant = 700. [static, private]

Definition at line 56 of file MadTVAnalysis.cxx. Referenced by EarlyActivity().

const Float_t MadTVAnalysis::kCalorEndZ = 7.1554 [static, private]

Definition at line 44 of file MadTVAnalysis.cxx.

const Float_t MadTVAnalysis::kHalfCalorEndZ = kTargEndZ+(kCalorEndZ-kTargEndZ)*0.5 [static, private]

Definition at line 45 of file MadTVAnalysis.cxx.

const Float_t MadTVAnalysis::kSpecEndZ = 16.656 [static, private]

Definition at line 46 of file MadTVAnalysis.cxx.

const Float_t MadTVAnalysis::kTargEndZ = 3.5914 [static, private]

Definition at line 43 of file MadTVAnalysis.cxx. Referenced by NDRadialFidVol().

const Float_t MadTVAnalysis::kVetoEndZ = 1.2154 [static, private]

Definition at line 42 of file MadTVAnalysis.cxx. Referenced by NDRadialFidVol().

const Float_t MadTVAnalysis::kXcenter = 1.4828 [static, private]

Definition at line 51 of file MadTVAnalysis.cxx. Referenced by CreatePAN(), InFidVol(), and NDRadialFidVol().

const Float_t MadTVAnalysis::kYcenter = 0.2384 [static, private]

Definition at line 52 of file MadTVAnalysis.cxx. Referenced by CreatePAN(), InFidVol(), and NDRadialFidVol().

MadNsID MadTVAnalysis::nsid

Definition at line 77 of file MadTVAnalysis.h. Referenced by CreatePAN().

bool MadTVAnalysis::openedabpidfile

Definition at line 81 of file MadTVAnalysis.h. Referenced by CreatePAN(), and MadTVAnalysis().

Anp::Interface* MadTVAnalysis::phnti [private]

Definition at line 171 of file MadTVAnalysis.h. Referenced by CreatePAN(), MadTVAnalysis(), and ~MadTVAnalysis().

MadQEID MadTVAnalysis::qeid

Definition at line 79 of file MadTVAnalysis.h. Referenced by CreatePAN().

string MadTVAnalysis::rustemconfigfile [private]

Definition at line 174 of file MadTVAnalysis.h. Referenced by CreatePAN(), and SetRustemConfig().

string MadTVAnalysis::rustempidfile_far [private]

Definition at line 173 of file MadTVAnalysis.h. Referenced by CreatePAN(), and SetRustemFarPID().

string MadTVAnalysis::rustempidfile_near [private]

Definition at line 172 of file MadTVAnalysis.h. Referenced by CreatePAN(), and SetRustemNearPID().

---

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

• MadTVAnalysis.h
• MadTVAnalysis.cxx

---