MadMKAnalysis Class Reference

#include <MadMKAnalysis.h>

Inheritance diagram for MadMKAnalysis:

List of all members.

Public Member Functions
void	CreatePAN (std::string s, const char *bmonpath="")
	MadMKAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
	MadMKAnalysis (JobC *, string, int)
virtual	~MadMKAnalysis ()
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	OnlyWriteFidEvents (bool x)
void	FillMCHists (TFile *fout)
void	DoInukeReweight (bool t)
void	SetNRandomSets (int i)
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
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)
Float_t	RecoShwEnergy (int ishw, const Detector::Detector_t &det, int mode=1, bool isdata=true)
Float_t	RecoMKMuEnergy (Int_t &opt, Int_t itrk, bool isdata=true, Detector::Detector_t det=Detector::kNear)
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)
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 ()
Private Attributes
Registry	fBECConfig
BeamType::BeamType_t	fMCBeam
bool	fDataUseMCBeam
bool	fDoInukeReweight
bool	fOnlyWriteFidEvents
int	fNRandomSets
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.4885
const Float_t	kYcenter = 0.1397
const Double_t	fEarlyActivityWindowSize = 10.
const Double_t	fPMTTimeConstant = 700.
const Int_t	fNPlaneEarlyAct = 30
const Double_t	fEarlyPlaneSphere = 0.5

---

Constructor & Destructor Documentation

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

Definition at line 192 of file MadMKAnalysis.cxx. References MadBase::Clear(), fBECConfig, fDataUseMCBeam, fDoInukeReweight, fMCBeam, fNRandomSets, fOnlyWriteFidEvents, and Registry::Set(). : MadAnalysis(chainSR, chainMC, chainTH, chainEM) { fMCBeam=BeamType::kUnknown; fDataUseMCBeam=false; fDoInukeReweight=false; fNRandomSets=0; fOnlyWriteFidEvents=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; }

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

Definition at line 228 of file MadMKAnalysis.cxx. References MadBase::Clear(), fBECConfig, fDataUseMCBeam, fDoInukeReweight, fMCBeam, fNRandomSets, fOnlyWriteFidEvents, and Registry::Set(). : MadAnalysis(j, path, entries) { fMCBeam=BeamType::kUnknown; fDataUseMCBeam=false; fDoInukeReweight=false; fNRandomSets=0; fOnlyWriteFidEvents=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); }

MadMKAnalysis::~MadMKAnalysis (   )  [virtual]

Definition at line 260 of file MadMKAnalysis.cxx. { }

---

Member Function Documentation

Float_t MadMKAnalysis::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 2231 of file MadMKAnalysis.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 MadMKAnalysis::ComputeLowPHRatio (   )  [protected]

Definition at line 3136 of file MadMKAnalysis.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 MadMKAnalysis::CreatePAN (  std::string  tag, const char *  bmonpath = "" )

BEGIN: intranuke reweighting init. END: intranuke reweighting init. fill event-by-event trees here Definition at line 265 of file MadMKAnalysis.cxx. References abs(), inuke_reweight::delta_fate::abs, Moments::AddPoint(), MadQEID::AltCalcQEID(), NtpSRPlane::beg, neugen_wrapper::begin_generation(), NtpSRPlane::begu, NtpSRPlane::begv, inuke_reweight::calc_weights(), MadDpID::CalcPID(), inuke_reweight::delta_fate::cex, NtpSRFitTrack::chi2, MadDpID::ChoosePDFs(), MadNsID::ChooseWeightFile(), NearbyVertexFinder::ClosestSpatial(), NearbyVertexFinder::ClosestTemporal(), NtpSRShower::clu, NtpSRDetStatus::coilstatus, NtpTHShower::completeslc, NtpTHTrack::completeslc, ComputeLowPHRatio(), CorrectMomentumFromRange(), CorrectSignedMomentumFromCurvature(), det, dpid, DupRecoStrips(), EarlyActivity(), inuke_reweight::delta_fate::elas, NtpSRTrack::end, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, FarTrkContained(), fBECConfig, BeamMonSpill::fBpmInt, fDataUseMCBeam, fDoInukeReweight, BeamMonSpill::fHadInt, BeamMonSpill::fHornCur, FillMCHists(), NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, NtpMCFluxInfo::fluxevtno, NtpMCFluxInfo::fluxrun, BeamMonSpill::fMuInt1, BeamMonSpill::fMuInt2, BeamMonSpill::fMuInt3, fNRandomSets, fOnlyWriteFidEvents, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamType::FromBeamMon(), inuke_reweight::delta_scale::ft, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTargProfX, BeamMonSpill::fTargProfY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTr101d, BeamMonSpill::fTrtgtd, inuke_reweight::generate_1sigma_shifts(), inuke_reweight::generate_uncor_shifts(), BDSpillAccessor::Get(), Registry::Get(), VldContext::GetDetector(), MadBase::GetEntry(), GetEvtTimeChar(), GetEvtTimeCharPHC(), GnumiInterface::GetParent(), MadNsID::GetPID(), RecPhysicsHeader::GetRemoteSpillType(), Moments::GetRMS(), RecDataHeader::GetRun(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), WeightCalculator::GetStandardConfig(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, gSystem(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSRCluster::id, inuke_reweight::delta_fate::inel, InFidVol(), MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::INu(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), kXcenter, kYcenter, NtpSREventSummary::litime, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), MadBase::LoadSlice(), BDSpillAccessor::LoadSpill(), MadBase::LoadTHShower(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), inuke_reweight::parameter_limits::lower, 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, NtpSRFitTrack::pass, MadAnalysis::PassBasicCuts(), MadQEID::PassMEDQECut(), MadQuantities::PassTrackCuts(), NtpMCFluxInfo::pdpx, NtpMCFluxInfo::pdpy, NtpMCFluxInfo::pdpz, NtpSRCluster::ph, NtpSRTrack::ph, NtpSREvent::ph, inuke_reweight::parameter_set::pi_fate, inuke_reweight::parameter_set::pi_scale, inuke_reweight::delta_fate::piprod, PittInFidVol(), PittTrkContained(), NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRCluster::planeview, inuke_reweight::parameter_set::pn_fate, inuke_reweight::parameter_set::pn_scale, NtpMCFluxInfo::ppdxdz, NtpMCFluxInfo::ppdydz, NtpMCFluxInfo::ppenergy, NtpMCFluxInfo::ppmedium, NtpMCFluxInfo::pppz, NtpMCFluxInfo::ppvx, NtpMCFluxInfo::ppvy, NtpMCFluxInfo::ppvz, Registry::Print(), inuke_reweight::parameter_limits::print(), inuke_reweight::parameter_set::print(), neugen_wrapper::print_configuration(), NearbyVertexFinder::ProcessEntry(), NtpMCFluxInfo::ptype, NtpTHShower::purity, MadQuantities::Q2(), qeid, NtpSRMomentum::qp, NtpSRMomentum::range, RecoMKMuEnergy(), RecoMuDCosNeuFD(), RecoMuDCosNeuND(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), RecoShwEnergy(), run(), BMSpillAna::SelectSpill(), inuke_reweight::delta_fate::set_abs_fates(), neugen_config::set_config_no(), MadDpID::SetPHCorrection(), BMSpillAna::SetSpill(), BMSpillAna::SetTimeDiff(), NtpSRShower::shwph, NtpSRPulseHeight::sigcor, SigInOut(), NtpSREvent::slc, NtpSRCluster::slope, GnumiInterface::Status(), 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, inuke_reweight::parameter_limits::upper, NtpSRVertex::v, NtpSRDmxStatus::validplanesfail, NtpSRDmxStatus::vertexplanefail, NtpSRTrack::vtx, NtpSRShower::vtx, NtpSREvent::vtx, NtpMCFluxInfo::vx, NtpMCFluxInfo::vy, NtpMCFluxInfo::vz, MadQuantities::W2(), NtpSRShowerPulseHeight::wtCCgev, NtpSRVertex::x, MadQuantities::X(), NtpMCFluxInfo::xpoint, inuke_reweight::delta_scale::xsec, 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); // stuff for bmpt reweighting... from tricia // only do this if mc GnumiInterface* gnumi=0; if(file_is_mc){ gnumi=new GnumiInterface(); if(!gnumi->Status()) { cout << "Warning: MadMKAnalysis::CreatePAN() - Flux files not open." << " Will not be filling NuParent info" << endl; cout << "Set environmental variable $GNUMIAUX to point to the " << "directory containing the gnumi files to fix this!" << endl; } else{ const char* gnumidir= getenv("GNUMIAUX"); cout<<"Read gnumi files from "<<gnumidir<<endl; } } //For Beam Reweighting: // must have this here, because it's in ntuple NuParent *nuparent = new NuParent(); //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 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 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_cev; // fully contained. true = 1, false = 0 Int_t is_mc; // is a corresponding mc event found Int_t pass_basic; // Pass basic cuts. true = 1, false = 0 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 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_dircosneu; // reco track vtx dircos w.r.t neutrino Float_t reco_dircosz; // reco track vtx z-dircos Float_t reco_vtxx; // reco x vtx Float_t reco_vtxy; // reco y vtx Float_t reco_vtxz; // reco z vtx Float_t 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; //demux failures (for fd) UChar_t dmx_stat; // std::cout<<"Second Print"<<std::endl; //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("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("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_cev",&is_cev,"is_cev/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("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); Float_t shw_pur, shw_comp; tree->Branch("shw_pur",&shw_pur,"shw_pur/F",32000); tree->Branch("shw_comp",&shw_comp,"shw_comp/F",32000); tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/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"); 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"); //demux status tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b"); tree->Branch("lowphrat",&lowphrat,"lowphrat/F"); // beam monitoring variables Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill; //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; //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; 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("hornI",&hornI,"hornI/D"); tree->Branch("closestspill",&closestspill,"closestspill/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"); //coil status Short_t coil_stat; tree->Branch("coil_stat",&coil_stat,"coil_stat/S"); //for bmpt reweighting tree->Branch("NuParent","NuParent",&nuparent,8000,1); //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; TTree* inukerw = 0; int n_inuke_wts=0; float* inuke_wts=0; double nucrad=-1.0; double wrad=1.0; std::vector<inuke_reweight::parameter_set> psets; if(file_is_mc && fDoInukeReweight){ inuke_reweight::parameter_set sigmas; neugen_config nc("MODBYRS"); nc.set_config_no(4); neugen_wrapper nw(&nc); nw.print_configuration(); gSystem->Setenv("NEUGEN3_XSEC","MODBYRS4_xsec.dat"); nw.begin_generation(); /* sigmas.pi_fate.inel=0.2; sigmas.pi_fate.abs=0.25; sigmas.pi_scale.ft=0.5; sigmas.pi_scale.xsec=0.1; // sigmas.pn_scale.ft=0.5; sigmas.pn_scale.xsec=0.1; */ sigmas.pi_fate.elas=0.1; sigmas.pi_fate.inel=0.4; sigmas.pi_fate.abs=0.3; sigmas.pi_fate.cex=0.49; sigmas.pi_fate.piprod=0.2; sigmas.pi_scale.xsec=0.1; sigmas.pn_fate.abs=0.2; sigmas.pn_fate.piprod=0.2; sigmas.pn_scale.xsec=0.1; sigmas.pi_scale.ft=0.5; // apply to both pi and p std::cout<<"Inuke Reweight sigmas: "<<std::endl; sigmas.print(std::cout); if(fNRandomSets<=0) { inuke_reweight::generate_1sigma_shifts(sigmas,psets); for (unsigned int iset=0; iset<psets.size(); iset++){ psets[iset].pi_fate.set_abs_fates(); psets[iset].pn_fate.set_abs_fates(); psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime // correlate change in pn xsec with change in pn elastic fate if(psets[iset].pn_scale.xsec!=0) { psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec; } } } else{ inuke_reweight::parameter_limits limits(sigmas,2.0); limits.upper.pi_scale.ft=0.5; // special case, 1 sigma, but uniform dist limits.lower.pi_scale.ft=-0.5; std::cout<<"Inuke Reweight Limits: "<<std::endl; limits.print(std::cout); gRandom->SetSeed(1); for(int i=0; i<fNRandomSets; i++){ inuke_reweight::parameter_set p; inuke_reweight::generate_uncor_shifts(sigmas,limits,p); psets.push_back(p); } // massage parameters for(int iset=0; iset<fNRandomSets; iset++){ psets[iset].pi_fate.set_abs_fates(); psets[iset].pn_fate.set_abs_fates(); psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime // correlate change in pn xsec with change in pn elastic fate if(psets[iset].pn_scale.xsec!=0) { psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec; } } } n_inuke_wts=psets.size(); inuke_wts = new float[n_inuke_wts]; // inukerw = new TTree("inukerw","intranuke weights"); inukerw=tree; inukerw->Branch("r",&run,"r/I"); inukerw->Branch("s",&snarl,"s/I"); inukerw->Branch("e",&event,"e/I"); inukerw->Branch("nwts",&n_inuke_wts,"nwts/I"); inukerw->Branch("wts",inuke_wts,"wts[nwts]/F"); inukerw->Branch("nucrad",&nucrad,"nucrad/D"); inukerw->Branch("wrad",&wrad,"wrad/D"); // save psets // 9*2 fates + 3*2 scales = 24 for(int ipar=0; ipar<24; ipar++){ TVectorD* vptr = new TVectorD(n_inuke_wts); inukerw->GetUserInfo()->Add(vptr); } for(int iset=0; iset<n_inuke_wts; iset++){ for(int ipar=0; ipar<24; ipar++){ TVectorD* vptr = dynamic_cast<TVectorD*>(inukerw->GetUserInfo()->At(ipar)); (*vptr)[iset]=psets[iset].get_par(ipar); } } } // map<VldTimeStamp, BeamMonTV>* beam_mon_map=0; // if(!file_is_mc) beam_mon_map = new BeamMonMap::MakeBeamMonMap(bmonpath); // map<VldTimeStamp, BeamMonTV> beam_mon_map // =BeamMonMap::MakeBeamMonMap(bmonpath); // std::cout<<"map file: "<<bmonpath<<std::endl; // std::cout<<"map size: "<<beam_mon_map.size()<<std::endl; 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); /* if(!checkeddb){ // ntpHeader->GetVldContext().Print(); // cout<<"IM HERE!!!!"<<endl; DbiResultPtr<BeamMonSpill> testbs(ntpHeader->GetVldContext(), Dbi::kDefaultTask, Dbi::kDisabled); // cout<<"tried to get the dbiresult ptr"<<endl; checkeddb=true; // cout<<"****************checking for beam monitoring db:"<<endl; if(testbs.GetNumRows()>0){ beamdb=true; cout<<"Found the Beam Monitoring Database Tables"<<endl; } else{ cout<<"Didn't find the Beam Monitoring Database Tables"<<endl; cout<<"Resorting to old beamsummary ntuples"<<endl; } } */ nevent = eventSummary->nevent; run = ntpHeader->GetRun(); if(run!=lastrun){ NRUNS++; lastrun = run; } snarl = ntpHeader->GetSnarl(); subrun = ntpHeader->GetSubRun(); trgsrc = ntpHeader->GetTrigSrc(); spilltype = ntpHeader->GetRemoteSpillType(); trgtime = eventSummary->trigtime; litime = eventSummary->litime; coil_stat = detStatus->coilstatus; 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(); // 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); detector = -1; //get detector type: const Detector::Detector_t det = ntpHeader->GetVldContext().GetDetector(); 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=0.0; beamcnf=0; nuTarZ=0.0; goodbeam=0; horn_on=target_in=n_batches=0; tor101=tr101d=tortgt=trtgtd=0.0; hadmon=mumon1=mumon2=mumon3=0.0; if(!file_is_mc){ VldContext vc = ntpHeader->GetVldContext(); VldTimeStamp vts = vc.GetTimeStamp(); // 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]); 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; 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; // } /* else{ // could the following return a reference? BeamMonTV btv = BeamMonMap::FindClosestSpill(beam_mon_map,vts); hbw=btv.hbw; vbw=btv.vbw; hpos2=btv.hpos2; vpos2=btv.vpos2; hornI=btv.hornI; closestspill=btv.closestspill; nuTarZ=btv.nuTarZ; //approximately assign beamcnf based on nuTarZ if(nuTarZ<20000&&nuTarZ>=0){ beamcnf=1; target_in=1; } else if(nuTarZ>20000&&nuTarZ<60000){ beamcnf=4; target_in=1; } else{ beamcnf=5; target_in=1; } if(fabs(hornI)>10){ horn_on=1; } n_batches=0; tor101=0.; tr101d=0.; tortgt=btv.bI; trtgtd=0.; hadmon=0.; mumon1=0.; mumon2=0.; mumon3=0.; } */ //set goodbeam==1 if it passes some basic beam monitoring cuts //changing widths to deal with rms instead of fits //getting rid of horn_on and target_in since they sometime //go haywire // if(horn_on&&target_in&&tortgt>0.1&&hbw<1.5&&vbw<2&& // hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(closestspill)<2){ /* if(tor101>0.1&&tortgt>0.1&&hbw<2.9&&vbw<2.9&& hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(hornI)>50&& fabs(closestspill)<2){ goodbeam=1; } */ //use Mark D. BMSpillAna to set good beam bmana.SetSpill(*bs); bmana.SetTimeDiff(closestspill); 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); } if(goodbeam==1&&spilltype!=3){//don't count fake spills 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; true_enu=0.0; true_emu=0.0; true_eshw=0.0; pass=0; is_fid=0; tree->Fill(); for(int iwt=0; iwt<n_inuke_wts; iwt++) inuke_wts[iwt]=-1.0; if(inukerw!=0 && inukerw!=tree) inukerw->Fill(); 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; //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; mcevent = -1; is_fid = 0; is_cev = 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_vtx_eshw=0; reco_vtx_wtd_eshw=0; reco_qe_enu = 0; reco_mushw=0; reco_qe_q2 = 0; reco_dircosneu = 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; // kinematic quantities always positive // convention: -1 indicates unfilled variable (NC events) reco_y = reco_x = reco_Q2 = reco_W2 = -1; shw_pur = shw_comp = -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; 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; is_fid = InFidVol(); //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); /* 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; } if(PassBasicCuts(event)) { //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); pass_track = PassTrackCuts(track_index); //methods should all be protected against index = -1 reco_emu = RecoMKMuEnergy(emu_meth,track_index); int shower_index = -1; // load the largest shower /* if(track_index!=-1) { LoadShowerAtTrackVertex(i,track_index,shower_index); } // shower with the best match to the vertex else LoadLargestShowerFromEvent(i,shower_index); */ /* if(nshower>0){ shower_index=ntpEvent->shw[0]; LoadShower(shower_index); } */ // use the jim shower! LoadShower_Jim(i,shower_index,det); reco_eshw = RecoShwEnergy(shower_index, det,1,!file_is_mc); if(shower_index>-1) reco_wtd_eshw = ntpShower->shwph.wtCCgev; reco_vtx_eshw = reco_eshw; reco_vtx_wtd_eshw = reco_wtd_eshw; // figure shower truth business out if(is_mc && shower_index>-1){ LoadTHShower(shower_index); assert(ntpTHShower); shw_pur=ntpTHShower->purity; shw_comp=ntpTHShower->completeslc; } // 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 is only useful for FD events 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 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; trkmomrange = CorrectMomentumFromRange(trkmomrange,!file_is_mc,det); trkqp = ntpTrack->momentum.qp; if(trkqp!=0){ trkqp = 1.0/CorrectSignedMomentumFromCurvature(1.0/trkqp, !file_is_mc, det); } 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; // 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 if(detector==1){//if it's the near detector if(pitt_evt_class>0){ // if it was classified int do_meth=0; // emu reco method we should use if( (pitt_evt_class == 1) || (pitt_evt_class == 3) ) do_meth=2; // stoppers --> range else if( (pitt_evt_class == 2) || (pitt_evt_class == 4) ) do_meth=1; // punch-through --> curvature if(do_meth==1 || do_meth==2){ reco_emu = RecoMKMuEnergy(do_meth,track_index, !file_is_mc,det); reco_enu = reco_emu+reco_eshw; } emu_meth=do_meth; // was forgetting about this } } else if(detector==2){//if it's the near detector int do_meth=FarTrkContained(ntpTrack->end.x, ntpTrack->end.y, ntpTrack->end.z); if(do_meth==1 || do_meth==2){ reco_emu = RecoMKMuEnergy(do_meth,track_index); reco_enu = reco_emu+reco_eshw; } 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); // CEDAR apparently has 5 plane long tracks, with ndf==0 Double_t tmp_ndof=ntpTrack->fit.ndof; if(tmp_ndof<=0.0) tmp_ndof=1.0; // pitt tracking quality cuts if( (ntpTrack->fit.pass>0) && (ntpTrack->fit.chi2/tmp_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){ dpid.SetPHCorrection(1.018); } if(dpid.ChoosePDFs(det,current_beam)){ dave_cc_pid = dpid.CalcPID(this,event,0); } 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(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); if(file_is_mc){ if(gnumi->Status()) { if(isST) gnumi->GetParent(strecord,mcevent,*nuparent); else gnumi->GetParent(mcrecord,mcevent,*nuparent); } } if(fDoInukeReweight && file_is_mc && inukerw!=0) { if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents){ std::vector<double> weights(n_inuke_wts,1.0); // std::cout<<"snarl, event: "<<snarl<<", "<<event<<std::endl; inuke_reweight::calc_weights(strecord,ntpTruth,psets,weights,nucrad,wrad,false); for (unsigned int iwt=0; iwt<weights.size(); iwt++) inuke_wts[iwt]=weights[iwt]; } // nuclear radius and weighting factor } if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents) { tree->Fill(); if(inukerw!=0 && inukerw!=tree) inukerw->Fill(); } } // loop over events } // loop over entries // delete nuparent; // should I do this? save.cd(); pottree->Fill(); pottree->Write(); tree->Write();if(inukerw!=0 && inukerw!=tree) inukerw->Write(); save.Write(); save.Close(); }

void MadMKAnalysis::DoInukeReweight (  bool  t  )  [inline]

Definition at line 150 of file MadMKAnalysis.h. References fDoInukeReweight. {fDoInukeReweight=t;}

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

Definition at line 3158 of file MadMKAnalysis.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 MadMKAnalysis::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 2802 of file MadMKAnalysis.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 MadMKAnalysis::FarTrkContained (  Float_t  x, Float_t  y, Float_t  z )  [static]

Definition at line 3118 of file MadMKAnalysis.cxx. Referenced by CreatePAN(), MadScanDisplay::DrawTextBox(), and MadEvDisplay::DrawTextBox(). { 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; }

void MadMKAnalysis::FillMCHists (  TFile *  fout  )

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

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

Definition at line 142 of file MadMKAnalysis.h. References fMCBeam. {fMCBeam=beam;}

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

Definition at line 187 of file MadMKAnalysis.h. { return fBECConfig; }

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

Definition at line 2636 of file MadMKAnalysis.cxx. References VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadStrip(), NtpSREvent::nstrip, NtpSREvent::stp, striptime, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSREventSummary::trigtime. Referenced by CreatePAN(). { //cut and pasted, with minor changes from MadDpAnalysis.cxx 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(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 = 0; if(striptime1>0 && striptime0<0) striptime=striptime1; if(striptime0>0 && striptime1<0) striptime=striptime0; if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; if(striptime<=timemin) timemin=striptime; if(striptime>=timemax) timemax=striptime; } } evttimemax=timemax; evttimemin=timemin; evttimeln=timemax-timemin; }

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

Definition at line 2675 of file MadMKAnalysis.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 = 0; if(striptime1>0 && striptime0<0) striptime=striptime1; if(striptime0>0 && striptime1<0) striptime=striptime0; if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; if(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& MadMKAnalysis::GetForcedMCBeam (   )  [inline]

Definition at line 143 of file MadMKAnalysis.h. {return fMCBeam;}

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

Definition at line 2072 of file MadMKAnalysis.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); double radius = pow(vx-kXcenter,2) + pow(vy-kYcenter,2); radius=sqrt(radius); // restrict vertices to target region and 1m radius /* if(vz>=kVetoEndZ && vz<=kTargEndZ && radius<=0.5) is_fid = true; */ // 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 MadMKAnalysis::InFidVol (   )  [protected, virtual]

Definition at line 2067 of file MadMKAnalysis.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 MadMKAnalysis::InFidVol (  Int_t  evidx  )  [protected, virtual]

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

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

Definition at line 2350 of file MadMKAnalysis.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 MadMKAnalysis::MuonShowerEnergy (  const NtpSREvent *  evt, const NtpSRTrack *  trk )  [protected]

Definition at line 2722 of file MadMKAnalysis.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 MadMKAnalysis::NDRadialFidVol (  Float_t  x, Float_t  y, Float_t  z )  [static]

Definition at line 2266 of file MadMKAnalysis.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 MadMKAnalysis::NPlanesInObj (  TObject *  rec, TObject *  obj, float  phcut, float &  sumph )  [protected]

Definition at line 2468 of file MadMKAnalysis.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]); if(!strip) continue; // 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]); if(!strip) continue; // 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 MadMKAnalysis::NStripsInObj (  TObject *  rec, TObject *  obj, float  phcut, float &  sumph )  [protected]

Definition at line 2557 of file MadMKAnalysis.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) continue; 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) continue; 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(); }

void MadMKAnalysis::OnlyWriteFidEvents (  bool  x  )  [inline]

Definition at line 145 of file MadMKAnalysis.h. References fOnlyWriteFidEvents. {fOnlyWriteFidEvents=x;}

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

Definition at line 2282 of file MadMKAnalysis.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 MadMKAnalysis::PittTrkContained (  Float_t  x, Float_t  y, Float_t  z, Float_t  u, Float_t  v )  [static]

Definition at line 2298 of file MadMKAnalysis.cxx. Referenced by CreatePAN(), MadScanDisplay::DrawTextBox(), and MadEvDisplay::DrawTextBox(). { // 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 MadMKAnalysis::RecoMKMuEnergy (  Int_t &  opt, Int_t  itrk, bool  isdata = true, Detector::Detector_t  det = Detector::kNear )  [protected]

Definition at line 2113 of file MadMKAnalysis.cxx. References CorrectMomentumFromRange(), CorrectSignedMomentumFromCurvature(), det, NtpSRFiducial::dr, NtpSRFiducial::dz, NtpSRTrack::fidall, MadBase::LoadTrack(), NtpSRTrack::momentum, NtpSRMomentum::qp, and NtpSRMomentum::range. Referenced by CreatePAN(). { const float mumass=0.10566; if(LoadTrack(itrk)){ float mr=ntpTrack->momentum.range; mr=CorrectMomentumFromRange(mr,isdata,det); float mc=0.0; // signed momentum from curvature if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp); mc=CorrectSignedMomentumFromCurvature(mc,isdata,det); if(opt==0){ //return the most appropriate measure of momentum // assign opt based on our choice if(ntpTrack->fidall.dr>0.5&&ntpTrack->fidall.dz>0.5) { opt=2; return sqrt(mr*mr+ mumass*mumass); } else { opt=1; // in R1.9 the tracker will apparently return (q/p)=0.0 // maybe it's when a track looks perfectly rigid? // if so, we have to do something // I don't want to use the range, that could be very wrong // but wrong in a more subtle way ... // so, we'll return an obviously ridiculous value of 10 TeV if(ntpTrack->momentum.qp == 0.0) return 10000.0; else return sqrt(mc*mc+mumass*mumass); } } else if(opt==1) { //return curvature measurement if(ntpTrack->momentum.qp == 0.0) return 10000.0; else return sqrt(mc*mc+mumass*mumass); } else if(opt==2) //return range measurement return sqrt(mr*mr + mumass*mumass); else return 0; } return 0.; }

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

Definition at line 2191 of file MadMKAnalysis.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 MadMKAnalysis::RecoMuDCosNeuND (  Int_t  itr = 0, Float_t *  vtx = 0 )  [protected, virtual]

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

Float_t MadMKAnalysis::RecoShwEnergy (  int  ishw, const Detector::Detector_t &  det, int  mode = 1, bool  isdata = true )  [protected]

Definition at line 2152 of file MadMKAnalysis.cxx. References CorrectShowerEnergy(), det, NtpSRStripPulseHeight::gev, NtpSRShowerPulseHeight::linCCgev, MadBase::LoadShower(), NtpSRShower::ph, and NtpSRShower::shwph. Referenced by CreatePAN(). { Float_t result = 0; // Return Jim's calculation Bool_t ok = LoadShower(opt); if(!ok) return 0; // test for shwph.linCCGeV<=0.0 // if so, assume that this is an R1.16 object // and use ph.gev result = ntpShower->shwph.linCCgev; if(result<=0.0) result=(ntpShower->ph.gev/1.23); else{ result = CorrectShowerEnergy(result,det,CandShowerHandle::kCC,mode,isdata); } return result; }

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

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

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

Definition at line 144 of file MadMKAnalysis.h. References fDataUseMCBeam. {fDataUseMCBeam=x;}

void MadMKAnalysis::SetNRandomSets (  int  i  )  [inline]

Definition at line 151 of file MadMKAnalysis.h. References fNRandomSets. {fNRandomSets=i;}

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

Definition at line 2454 of file MadMKAnalysis.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 MadMKAnalysis::SigInOut (  Float_t &  sigfull, Float_t &  sigpart, Float_t &  trkfull, Float_t &  trkpart )  [protected]

Definition at line 2410 of file MadMKAnalysis.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

MadDpID MadMKAnalysis::dpid

Definition at line 154 of file MadMKAnalysis.h. Referenced by CreatePAN().

Registry MadMKAnalysis::fBECConfig [private]

Definition at line 235 of file MadMKAnalysis.h. Referenced by CreatePAN(), MadMKAnalysis(), and SetBECFile().

bool MadMKAnalysis::fDataUseMCBeam [private]

Definition at line 237 of file MadMKAnalysis.h. Referenced by CreatePAN(), MadMKAnalysis(), and SetDataUseMCBeam().

bool MadMKAnalysis::fDoInukeReweight [private]

Definition at line 238 of file MadMKAnalysis.h. Referenced by CreatePAN(), DoInukeReweight(), and MadMKAnalysis().

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

Definition at line 187 of file MadMKAnalysis.cxx.

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

Definition at line 190 of file MadMKAnalysis.cxx.

BeamType::BeamType_t MadMKAnalysis::fMCBeam [private]

Definition at line 236 of file MadMKAnalysis.h. Referenced by ForceMCBeam(), and MadMKAnalysis().

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

Definition at line 189 of file MadMKAnalysis.cxx.

int MadMKAnalysis::fNRandomSets [private]

Definition at line 240 of file MadMKAnalysis.h. Referenced by CreatePAN(), MadMKAnalysis(), and SetNRandomSets().

bool MadMKAnalysis::fOnlyWriteFidEvents [private]

Definition at line 239 of file MadMKAnalysis.h. Referenced by CreatePAN(), MadMKAnalysis(), and OnlyWriteFidEvents().

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

Definition at line 188 of file MadMKAnalysis.cxx. Referenced by EarlyActivity().

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

Definition at line 181 of file MadMKAnalysis.cxx.

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

Definition at line 182 of file MadMKAnalysis.cxx.

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

Definition at line 183 of file MadMKAnalysis.cxx.

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

Definition at line 180 of file MadMKAnalysis.cxx. Referenced by NDRadialFidVol().

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

Definition at line 179 of file MadMKAnalysis.cxx. Referenced by NDRadialFidVol().

const Float_t MadMKAnalysis::kXcenter = 1.4885 [static, private]

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

const Float_t MadMKAnalysis::kYcenter = 0.1397 [static, private]

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

MadNsID MadMKAnalysis::nsid

Definition at line 153 of file MadMKAnalysis.h. Referenced by CreatePAN().

MadQEID MadMKAnalysis::qeid

Definition at line 155 of file MadMKAnalysis.h. Referenced by CreatePAN().

---

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

• MadMKAnalysis.h
• MadMKAnalysis.cxx

---