MadMKAnalysis Class Reference

#include <MadMKAnalysis.h>

Inheritance diagram for MadMKAnalysis:
MadAnalysis MadQuantities MadBase

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_tGetForcedMCBeam ()
void SetDataUseMCBeam (bool x)
void OnlyWriteFidEvents (bool x)
void FillMCHists (TFile *fout)
void DoInukeReweight (bool t)
void SetNRandomSets (int i)

Static Public Member Functions

static Bool_t PittInFidVol (Float_t x, Float_t y, Float_t z, Float_t u, Float_t v)
static Int_t PittTrkContained (Float_t x, Float_t y, Float_t z, Float_t u, Float_t v)
static Int_t InPartialRegion (UShort_t plane, UShort_t strip)
static Int_t NDRadialFidVol (Float_t x, Float_t y, Float_t z)
static 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)
RegistryGetBECRegistry ()
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

static const Float_t kVetoEndZ = 1.2154
static const Float_t kTargEndZ = 3.5914
static const Float_t kCalorEndZ = 7.1554
static const Float_t kHalfCalorEndZ = kTargEndZ+(kCalorEndZ-kTargEndZ)*0.5
static const Float_t kSpecEndZ = 16.656
static const Float_t kXcenter = 1.4885
static const Float_t kYcenter = 0.1397
static const Double_t fEarlyActivityWindowSize = 10.
static const Double_t fPMTTimeConstant = 700.
static const Int_t fNPlaneEarlyAct = 30
static const Double_t fEarlyPlaneSphere = 0.5

Detailed Description

Definition at line 112 of file MadMKAnalysis.h.


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(), MadBase::emrecord, fBECConfig, fDataUseMCBeam, fDoInukeReweight, fMCBeam, fNRandomSets, fOnlyWriteFidEvents, MadBase::isEM, MadBase::isMC, MadBase::isTH, BeamType::kLE, Detector::kNear, BeamType::kUnknown, MadBase::mcrecord, MadBase::Nentries, MadBase::record, Registry::Set(), MadBase::threcord, and MadBase::whichSource.

00194   : MadAnalysis(chainSR, chainMC, chainTH, chainEM)
00195 {
00196   fMCBeam=BeamType::kUnknown;
00197   fDataUseMCBeam=false;
00198   fDoInukeReweight=false;
00199   fNRandomSets=0;
00200   fOnlyWriteFidEvents=false;
00201   if(!chainSR) {
00202     record = 0;
00203     emrecord = 0;
00204     mcrecord = 0;
00205     threcord = 0;
00206     Clear();
00207     whichSource = -1;
00208     isMC = true;
00209     isTH = true;
00210     isEM = true;
00211     Nentries = -1;
00212     return;
00213   }
00214   static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
00215   fBECConfig.Set("beam:flux_file",default_bec_file);
00216   fBECConfig.Set("beam:beam_type",BeamType::kLE);
00217   fBECConfig.Set("beam:detector",Detector::kNear);
00218   fBECConfig.Set("beam:low_energy_limit",0.5);
00219   fBECConfig.Set("beam:high_energy_limit",60.0);
00220 
00221   
00222   //  InitChain(chainSR,chainMC,chainTH,chainEM);
00223   whichSource = 0;  
00224 
00225 }

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

Definition at line 228 of file MadMKAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, fBECConfig, MadBase::fChain, fDataUseMCBeam, fDoInukeReweight, MadBase::fJC, fMCBeam, fNRandomSets, fOnlyWriteFidEvents, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, BeamType::kLE, Detector::kNear, BeamType::kUnknown, MadBase::mcrecord, MadBase::Nentries, MadBase::record, Registry::Set(), MadBase::threcord, and MadBase::whichSource.

00228                                                             : MadAnalysis(j, path, entries)
00229 {
00230   fMCBeam=BeamType::kUnknown;
00231   fDataUseMCBeam=false;
00232   fDoInukeReweight=false;
00233   fNRandomSets=0;
00234   fOnlyWriteFidEvents=false;
00235   record = 0;
00236   emrecord = 0;
00237   mcrecord = 0;
00238   threcord = 0;
00239   Clear();
00240   isMC = true;
00241   isTH = true;
00242   isEM = true;
00243   Nentries = entries;
00244   whichSource = 1;
00245   jcPath = path;
00246   fJC = j;
00247   fChain = NULL;
00248 
00249 
00250   static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
00251   fBECConfig.Set("beam:flux_file",default_bec_file);
00252   fBECConfig.Set("beam:beam_type",BeamType::kLE);
00253   fBECConfig.Set("beam:detector",Detector::kNear);
00254   fBECConfig.Set("beam:low_energy_limit",0.5);
00255   fBECConfig.Set("beam:high_energy_limit",60.0);
00256 
00257 }

MadMKAnalysis::~MadMKAnalysis (  )  [virtual]

Definition at line 260 of file MadMKAnalysis.cxx.

00261 {
00262 
00263 }


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.

02237                                                                  {
02238   // given a line with end points (x1,y1) (x2,y2) 
02239   // and a spatial point (xpt,ypt)
02240   // report the closest point on the line (x,y) to the spatial point
02241   // and compute the distance from (xpt,ypt) to (x,y)
02242   Float_t u= (xpt-x1)*(x2-x1)+(ypt-y1)*(y2-y1);
02243   u/=sqrt( pow(x2-x1,2) + pow(y2-y1,2));
02244   x=x1+u*(x2-x1);
02245   y=y1+u*(y2-y1);
02246   Float_t dist =sqrt(pow(xpt-x, 2)+pow(ypt-y,2));
02247   return dist;
02248 }

float MadMKAnalysis::ComputeLowPHRatio (  )  [protected]

Definition at line 3136 of file MadMKAnalysis.cxx.

References MadBase::GetEventSummary(), MadBase::LoadStrip(), NtpSREventSummary::nstrip, MadBase::ntpStrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, and NtpSRStrip::ph1.

Referenced by CreatePAN().

03137 {
03138   float lowphsum=0.;
03139   float totphsum=0.;
03140 
03141   for(unsigned int i=0;i<GetEventSummary()->nstrip;i++){
03142     if(!LoadStrip(i)) continue;
03143     totphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe;
03144     if(ntpStrip->ph1.pe+ntpStrip->ph0.pe<2){
03145       lowphsum+=ntpStrip->ph1.pe+ntpStrip->ph0.pe;
03146     }
03147   }
03148 
03149   float result=-1.;
03150   if(totphsum!=0){
03151     result = lowphsum/totphsum;
03152   }
03153 
03154   return result;
03155 
03156 }

void MadMKAnalysis::CreatePAN ( std::string  s,
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 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(), EnergyCorrections::CorrectMomentumFromRange(), EnergyCorrections::CorrectSignedMomentumFromCurvature(), det, MadBase::detStatus, MadBase::dmxStatus, dpid, DupRecoStrips(), EarlyActivity(), inuke_reweight::delta_fate::elas, NtpSRPlane::end, NtpSRTrack::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, MadBase::eventSummary, FarTrkContained(), fBECConfig, BeamMonSpill::fBpmInt, fDataUseMCBeam, fDoInukeReweight, BeamMonSpill::fHadInt, BeamMonSpill::fHornCur, FillMCHists(), NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, NtpMCFluxInfo::fluxevtno, NtpMCFluxInfo::fluxrun, fMCBeam, 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(), MadBase::isST, Detector::kFar, SimFlag::kMC, Detector::kNear, PlaneView::kU, BeamType::kUnknown, PlaneView::kV, 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, MadBase::mcrecord, 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, MadBase::Nentries, NtpTHTrack::neumc, NtpSREventSummary::nevent, NtpMCFluxInfo::nimpwt, NtpSRDmxStatus::nonphysicalfail, NtpMCFluxInfo::norig, NPlanesInObj(), NtpMCFluxInfo::npz, NtpSREvent::nshower, nsid, NtpSREvent::nstrip, NtpSRTrack::nstrip, NStripsInObj(), MadBase::ntpCluster, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpSlice, MadBase::ntpTHShower, MadBase::ntpTHTrack, MadBase::ntpTrack, MadBase::ntpTruth, 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_set::print(), inuke_reweight::parameter_limits::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(), MadBase::record, RecoShwEnergy(), run(), BMSpillAna::SelectSpill(), neugen_config::set_config_no(), MadDpID::SetPHCorrection(), BMSpillAna::SetSpill(), BMSpillAna::SetTimeDiff(), NtpSRShower::shwph, NtpSRPulseHeight::sigcor, SigInOut(), NtpSREvent::slc, NtpSRCluster::slope, GnumiInterface::Status(), MadBase::strecord, 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, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, NtpMCFluxInfo::vx, NtpMCFluxInfo::vy, NtpMCFluxInfo::vz, MadQuantities::W2(), NtpSRShowerPulseHeight::wtCCgev, MadQuantities::X(), NtpSRVertex::x, NtpMCFluxInfo::xpoint, inuke_reweight::delta_scale::xsec, MadQuantities::Y(), NtpSRVertex::y, NtpMCFluxInfo::ypoint, NtpSRVertex::z, and NtpMCFluxInfo::zpoint.

00266 {
00267   
00268   // is this MC??
00269   this->GetEntry(0);
00270   const bool file_is_mc
00271     =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
00272   // stuff for bmpt reweighting... from tricia
00273   // only do this if mc
00274 
00275   GnumiInterface* gnumi=0;  
00276   if(file_is_mc){
00277     gnumi=new GnumiInterface();
00278     
00279     if(!gnumi->Status()) {
00280       cout << "Warning: MadMKAnalysis::CreatePAN() - Flux files not open."
00281            << " Will not be filling NuParent info"
00282            << endl;
00283       cout << "Set environmental variable $GNUMIAUX to point to the "
00284            << "directory containing the gnumi files to fix this!"
00285            << endl;    
00286     }
00287     else{
00288       const char* gnumidir= getenv("GNUMIAUX");
00289       cout<<"Read gnumi files from "<<gnumidir<<endl;
00290     }
00291   }
00292   //For Beam Reweighting:
00293   // must have this here, because it's in ntuple
00294   NuParent *nuparent = new NuParent();
00295 
00296   
00297   //PAN Quantities 
00298   //Truth:
00299   Float_t true_enu;       // true neutrino energy (GeV)
00300   Float_t true_emu;       // true muon energy
00301   Float_t true_eshw;      // true shower energy
00302   Float_t true_x;         // true x
00303   Float_t true_y;         // true y
00304   Float_t true_q2;        // true q2
00305   Float_t true_w2;        // true w2
00306   Float_t true_dircosneu; // true muon dircos w.r.t neutrino
00307   Float_t true_dircosz;   // track z-direction cosine
00308   Float_t true_vtxx;      // true x vtx
00309   Float_t true_vtxy;      // true y vtx
00310   Float_t true_vtxz;      // true z vtx
00311   Int_t true_pitt_fid; // was event truly in fiducial volue?
00312   Float_t true_trk_cmplt;  //true track completeness
00313 
00314   //For Neugen:
00315   Float_t flavour;        // true flavour: 1-e 2-mu 3-tau
00316   Int_t process;          // process: 1001-QEL 1002-RES 1003-DIS 1004-COH
00317   Int_t nucleus;          // target nucleus: 274-C 284-O 372-Fe
00318   Int_t cc_nc;            // cc/nc flag: 1-cc 2-nc
00319   Int_t initial_state;    // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ...
00320   Int_t inu;              // stdhep inu
00321   Int_t resnum; // IdHEP%1000 of resonance state
00322 
00323   Short_t npp; // # final state pi+
00324   Short_t npm; // # fs pi-
00325   Short_t npz; // # fs pi0
00326   Short_t np; // # fs p
00327   Short_t nn; // # fs n
00328   float epp; // energy final state pi+
00329   float epm; // energy fs pi-
00330   float epz; // energy fs pi0
00331   float ep; // energy fs p
00332   float en; // energy fs n
00333 
00334   
00335   //Reco Quantities
00336   Int_t detector;         // Near = 1, Far = 2;
00337   Int_t run;              // run number
00338   Int_t snarl;            // snarl number
00339   Int_t subrun;           // subrun number
00340   UInt_t trgsrc;          // trigger source
00341   double trgtime;         // trigger time
00342   Int_t spilltype;        // comes from mdSpillTypeEnum
00343                           // 0=none, 1=reported, 2=predicted
00344                           // 3=fake, 4=local
00345 
00346   Int_t nevent;           // number of events in snarl
00347   Int_t event;            // event index
00348   Int_t mcevent;          // mc event index
00349   Int_t ntrack;           // number of tracks in event
00350   Int_t nshower;          // number of showers in event
00351   Int_t nstrip;           // number of strips in event
00352 
00353   Int_t is_fid;           // pass fid cut. true = 1, false = 0
00354   Int_t is_cev;           // fully contained. true = 1, false = 0 
00355   Int_t is_mc;            // is a corresponding mc event found
00356   Int_t pass_basic;       // Pass basic cuts. true = 1, false = 0
00357   Int_t is_pitt_fid;      // passes pitt ND vtx fid cut true=1, false=0
00358   Int_t is_trk_fid;       // is the trkvtx in fid (pitt fid for ND)
00359                           // true=1, false=0, notrack = -1
00360   Int_t nd_radial_fid;    // a bitfield, bits correspond to r cuts
00361   Int_t pitt_evt_class;   // CC event class as defined by pitt group
00362   // 1 = track is fully contained in the upstream region
00363   // 2 = track is partially contained in the upstream region
00364   // 3 = track is fully contained in the downstream region
00365   // 4 = track is partially contained in the downstream region
00366   Int_t pitt_trk_qual;   // pitt tracking quality cuts
00367   Int_t duvvtx;          //delta (Uvtx-Vvtx)
00368 
00369   Float_t pid0;           // pid parameter (usual method)
00370   Float_t pid1;           // pid parameter (alternative method 1)
00371   Float_t pid2;           // pid parameter (alternative method 2)
00372   Int_t pass_track;       // the primary track passes track cuts
00373   Int_t pass;             // pass analysis cuts. true = 1, false = 0
00374 
00375   Int_t emu_meth;         // method used to determine muon energy
00376   Float_t reco_enu;       // reco neutrino energy
00377   Float_t reco_emu;       // reco muon energy
00378   Float_t reco_mushw;     // showers along muon track
00379   Float_t reco_eshw;      // reco shower energy (largest shower
00380   Float_t reco_wtd_eshw;    // as above but weighted 
00381   Float_t reco_vtx_eshw;    // reco shower energy (=0 if no vertex shower)
00382   Float_t reco_vtx_wtd_eshw;// as above but weighted (=0 if no vertex shower)
00383   
00384   Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino
00385   Float_t reco_dircosz;   // reco track vtx z-dircos
00386 
00387   Float_t reco_vtxx;      // reco x vtx
00388   Float_t reco_vtxy;      // reco y vtx
00389   Float_t reco_vtxz;      // reco z vtx
00390   Float_t reco_vtxr;// radial position of vertex
00391 
00392   Float_t reco_trk_vtxx;      // reco x track vtx
00393   Float_t reco_trk_vtxy;      // reco y track vtx
00394   Float_t reco_trk_vtxz;      // reco z track vtx
00395   Float_t reco_trk_vtxr;// radial position of vertex
00396 
00397   Float_t reco_shw_vtxx;      // reco x track vtx
00398   Float_t reco_shw_vtxy;      // reco y track vtx
00399   Float_t reco_shw_vtxz;      // reco z track vtx
00400   Float_t reco_shw_vtxr;// radial position of vertex
00401 
00402 
00403   Float_t reco_trk_endx;      // reco x track end
00404   Float_t reco_trk_endy;      // reco y track end
00405   Float_t reco_trk_endz;      // reco z track end
00406   Float_t reco_trk_endr;// radial position of vertex
00407 
00408   Int_t evt_nplanes;        //number of planes in event above some ph
00409   Int_t trk_nplanes;        //number of planes in track above some ph
00410   Int_t slc_nplanes;        //number of planes in slice above some ph
00411   Int_t shw_nplanes_cut;        //number of planes in shower above some ph
00412   Int_t shw_nplanes_raw;        //number of planes in shower
00413   Float_t shw_ph_cut, shw_ph_raw; // sigcor in shower, w/ & w/o ph cut
00414   
00415   Int_t evt_nstrips;        //number of strips in event above some ph
00416   Int_t trk_nstrips;        //number of strips in track above some ph
00417   Int_t slc_nstrips;        //number of strips in slice above some ph
00418 
00419 
00420   // reconstructed kinematics
00421   // formulae given for high-energy CC interactions
00422   // Assume nu = Ehad
00423   Float_t reco_y;// y = (Enu-Emu)/Enu
00424   Float_t reco_Q2;// Q2 = -q^{2} = 2 Enu Emu (1 - cos(theta_mu))
00425   Float_t reco_x;// x = Q2 / (2M* Ehad)  {M=nucleon mass}
00426   Float_t reco_W2;// W2 = M^{2} + q^{2} + 2M(Enu-Emu)
00427 
00428   // reco quantities, assuming a QE event
00429   Float_t reco_qe_enu;    // reco qe neutrino energy
00430   Float_t reco_qe_q2;     // reco qe q2
00431 
00432   
00433   Float_t evtlength;      // reco event length
00434   Float_t shwlength;      // reco shower length
00435   Float_t trklength;      // reco track length (planes)
00436   Int_t ntrklike;      // number of track-like planes
00437   Float_t trkmomrange;    // reco track momentum from range
00438   Float_t trkqp;          // reco track q/p
00439   Float_t trkeqp;         // reco track q/p error
00440   Float_t trkphfrac;      // reco pulse height fraction in track
00441   Float_t trkphplane;     // reco average track pulse height/plane
00442   Float_t trkchi2;
00443   Int_t trkndf;
00444 
00445 
00446   Int_t trkend, trkendu, trkendv; // track end plane, u, v
00447   Int_t shwend,shwbeg; // shower begin,end planes
00448 
00449   Int_t nss[6];//
00450   Float_t ess[6];//
00451   Int_t ntotss;
00452 
00453 
00454 
00455   // variables to tag nearest event  in space and time
00456   // nvsi = nearest vertex space : index
00457   // nvsd = nearest vertex space : distance
00458   // nvst = nearest vertex space : distance
00459   // nvsp = nearest vertex space : pulse height
00460   // nvsevt = nearest vertex space : mc event
00461   // nvti = nearest vertex time : index
00462   // ....  
00463   // bvsi = back vertex space : index
00464   // the "back" variables hold the same info as the "near" variables
00465   // but for that nearest event
00466   // I'm interested in this case
00467   //  
00468   //  (this vtx)------>(nearest)--->(nearest to it)
00469   //
00470   // if you are going to try to recombine events
00471   // you need to somehow account for this
00472   //
00473 
00474   Int_t nvsi, bvsi, nvti, bvti;
00475   Int_t nvsevt, bvsevt, nvtevt, bvtevt;
00476   Float_t nvsd, bvsd, nvtd, bvtd;
00477   Float_t nvst, bvst, nvtt, bvtt;
00478   Float_t nvsp, bvsp, nvtp, bvtp;
00479 
00480   // this event's pulseheight (in sigcor)
00481   Float_t evtph;
00482   Float_t evtphgev;
00483 
00484   // the pulse height of the event and track
00485   // in the fully and partially instrumented regions
00486   Float_t evtsigfull,trksigfull,evtsigpart,trksigpart;
00487 
00488 
00489   //time structure of event
00490   double evttimemin, evttimemax, evttimeln;
00491   double evttimeminphc, evttimemaxphc, evttimelnphc;
00492 
00493   //early activity variables
00494   float earlywph, earlywphpw, earlywphsphere;
00495   float ph1000, ph500, ph100;
00496   float lateph1000, lateph500, lateph100, lateph50;
00497 
00498   //duplicate reconstructed strip variables
00499   int nduprs; //number of reco strips that are duplicated in another event
00500   float dupphrs; //pulseheight of reco strips 
00501                  //that are duplicated in another event
00502 
00503   //time last li event in snarl, -1 if none;
00504   double litime;
00505 
00506   //demux failures (for fd)
00507   UChar_t dmx_stat;
00508   //  std::cout<<"Second Print"<<std::endl;
00509 
00510   //brian's cut on lowph ratio to remove junk events 
00511   //in FD that pass the dmx_stat variable
00512   float lowphrat;
00513 
00514   //Trees
00515   // output file
00516   std::string savename = "PAN_" + tag + ".root";
00517   TFile save(savename.c_str(),"RECREATE"); 
00518   save.cd();
00519   TTree *tree = new TTree("pan","pan");
00520   //Truth
00521   tree->Branch("true_enu",&true_enu,"true_enu/F",32000);
00522   tree->Branch("true_emu",&true_emu,"true_emu/F",32000);
00523   tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000);
00524   tree->Branch("true_x",&true_x,"true_x/F",32000);
00525   tree->Branch("true_y",&true_y,"true_y/F",32000);
00526   tree->Branch("true_q2",&true_q2,"true_q2/F",32000);
00527   tree->Branch("true_w2",&true_w2,"true_w2/F",32000);
00528   Float_t nu_px, nu_py, nu_pz;
00529   Float_t tar_px, tar_py, tar_pz, tar_e;
00530   tree->Branch("nu_px", &nu_px, "nu_px/F");
00531   tree->Branch("nu_py", &nu_py, "nu_py/F");
00532   tree->Branch("nu_pz", &nu_pz, "nu_pz/F");
00533   tree->Branch("tar_px", &tar_px, "tar_px/F");
00534   tree->Branch("tar_py", &tar_py, "tar_py/F");
00535   tree->Branch("tar_pz", &tar_pz, "tar_pz/F");
00536   tree->Branch("tar_e", &tar_e, "tar_e/F");
00537 
00538 
00539   tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000);
00540   tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000);
00541   tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000);
00542   tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000);
00543   tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000);
00544   tree->Branch("true_pitt_fid",&true_pitt_fid,"true_pitt_fid/I",32000);
00545   tree->Branch("true_trk_cmplt",&true_trk_cmplt,"true_trk_cmplt/F",32000);
00546 
00547   //Neugen
00548   tree->Branch("flavour",&flavour,"flavour/F",32000);
00549   tree->Branch("process",&process,"process/I",32000);
00550   tree->Branch("nucleus",&nucleus,"nucleus/I",32000);
00551   tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000);
00552   tree->Branch("inu",&inu,"inu/I",32000);
00553 
00554   tree->Branch("initial_state",&initial_state,"initial_state/I",32000);
00555   tree->Branch("resnum",&resnum,"resnum/I",32000);
00556   tree->Branch("npp",&npp,"npp/S",32000);
00557   tree->Branch("npm",&npm,"npm/S",32000);
00558   tree->Branch("npz",&npz,"npz/S",32000);
00559   tree->Branch("np",&np,"np/S",32000);
00560   tree->Branch("nn",&nn,"nn/S",32000);
00561 
00562   tree->Branch("epp",&epp,"epp/F",32000);
00563   tree->Branch("epm",&epm,"epm/F",32000);
00564   tree->Branch("epz",&epz,"epz/F",32000);
00565   tree->Branch("ep",&ep,"ep/F",32000);
00566   tree->Branch("en",&en,"en/F",32000);
00567 
00568   //Reco
00569   tree->Branch("detector",&detector,"detector/I",32000);
00570   tree->Branch("run",&run,"run/I",32000);
00571   tree->Branch("snarl",&snarl,"snarl/I",32000);
00572   tree->Branch("subrun",&subrun,"subrun/I",32000);
00573   tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000);
00574   tree->Branch("trgtime",&trgtime,"trgtime/D",32000);
00575   tree->Branch("spilltype",&spilltype,"spiltype/I",32000);
00576   tree->Branch("nevent",&nevent,"nevent/I",32000);
00577   tree->Branch("event",&event,"event/I",32000);
00578   tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
00579   tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
00580   tree->Branch("nshower",&nshower,"nshower/I",32000);
00581   tree->Branch("nstrip",&nstrip,"nstrip/I",32000);
00582   
00583   tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
00584   tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
00585   tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
00586   tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
00587   tree->Branch("pass_track",&pass_track,"pass_track/I",32000);
00588   tree->Branch("emu_meth",&emu_meth,"emu_meth/I",32000);
00589   tree->Branch("is_pitt_fid",&is_pitt_fid,"is_pitt_fid/I",32000);
00590   tree->Branch("is_trk_fid",&is_trk_fid,"is_trk_fid/I",32000);
00591 
00592   tree->Branch("nd_radial_fid",&nd_radial_fid,"nd_radial_fid/I",32000);
00593   tree->Branch("pitt_evt_class",&pitt_evt_class,"pitt_evt_class/I",32000);
00594   tree->Branch("pitt_trk_qual",&pitt_trk_qual,"pitt_trk_qual/I",32000);
00595   tree->Branch("duvvtx",&duvvtx,"duvvtx/I",32000);
00596 
00597   tree->Branch("pid0",&pid0,"pid0/F",32000);
00598   tree->Branch("pid1",&pid1,"pid1/F",32000);
00599   tree->Branch("pid2",&pid2,"pid2/F",32000);
00600   tree->Branch("pass",&pass,"pass/I",32000);
00601 
00602   float med_qe_pid;
00603   int med_qe_pass;
00604   tree->Branch("med_qe_pid",&med_qe_pid,"med_qe_pid/F",32000);
00605   tree->Branch("med_qe_pass",&med_qe_pass,"med_qe_pass/I",32000);
00606 
00607   double niki_cc_pid;
00608   tree->Branch("niki_cc_pid",&niki_cc_pid,"niki_cc_pid/D",32000);
00609   float dave_cc_pid;
00610   tree->Branch("dave_cc_pid",&dave_cc_pid,"dave_cc_pid/F",32000);
00611 
00612 
00613   tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
00614   tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
00615   tree->Branch("reco_mushw",&reco_mushw,"reco_mushw/F",32000);
00616   tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
00617   tree->Branch("reco_wtd_eshw",&reco_wtd_eshw,"reco_wtd_eshw/F",32000);
00618   tree->Branch("reco_vtx_eshw",&reco_vtx_eshw,"reco_vtx_eshw/F",32000);
00619   tree->Branch("reco_vtx_wtd_eshw",&reco_vtx_wtd_eshw,"reco_vtx_wtd_eshw/F",32000);
00620   Float_t shw_pur, shw_comp;
00621   tree->Branch("shw_pur",&shw_pur,"shw_pur/F",32000);
00622   tree->Branch("shw_comp",&shw_comp,"shw_comp/F",32000);
00623 
00624 
00625 
00626   tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
00627   tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);
00628 
00629   tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
00630   tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
00631   tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
00632   tree->Branch("reco_vtxr",&reco_vtxr,"reco_vtxr/F",32000);
00633 
00634   tree->Branch("reco_trk_vtxx",&reco_trk_vtxx,"reco_trk_vtxx/F",32000);
00635   tree->Branch("reco_trk_vtxy",&reco_trk_vtxy,"reco_trk_vtxy/F",32000);
00636   tree->Branch("reco_trk_vtxz",&reco_trk_vtxz,"reco_trk_vtxz/F",32000);
00637   tree->Branch("reco_trk_vtxr",&reco_trk_vtxr,"reco_trk_vtxr/F",32000);
00638 
00639   tree->Branch("reco_shw_vtxx",&reco_shw_vtxx,"reco_shw_vtxx/F",32000);
00640   tree->Branch("reco_shw_vtxy",&reco_shw_vtxy,"reco_shw_vtxy/F",32000);
00641   tree->Branch("reco_shw_vtxz",&reco_shw_vtxz,"reco_shw_vtxz/F",32000);
00642   //  tree->Branch("reco_shw_vtxr",&reco_shw_vtxr,"reco_shw_vtxr/F",32000);
00643 
00644   tree->Branch("reco_trk_endx",&reco_trk_endx,"reco_trk_endx/F",32000);
00645   tree->Branch("reco_trk_endy",&reco_trk_endy,"reco_trk_endy/F",32000);
00646   tree->Branch("reco_trk_endz",&reco_trk_endz,"reco_trk_endz/F",32000);
00647   tree->Branch("reco_trk_endr",&reco_trk_endr,"reco_trk_endr/F",32000);
00648 
00649   tree->Branch("evt_nplanes",&evt_nplanes,"evt_nplanes/I",32000);
00650   tree->Branch("trk_nplanes",&trk_nplanes,"trk_nplanes/I",32000);
00651   tree->Branch("slc_nplanes",&slc_nplanes,"slc_nplanes/I",32000);
00652   tree->Branch("shw_nplanes_cut",&shw_nplanes_cut,"shw_nplanes_cut/I",32000);
00653   tree->Branch("shw_nplanes_raw",&shw_nplanes_raw,"shw_nplanes_raw/I",32000);
00654 
00655   tree->Branch("shw_ph_cut",&shw_ph_cut,"shw_ph_cut/F",32000);
00656   tree->Branch("shw_ph_raw",&shw_ph_raw,"shw_ph_raw/F",32000);
00657 
00658 
00659   tree->Branch("evt_nstrips",&evt_nstrips,"evt_nstrips/I",32000);
00660   tree->Branch("trk_nstrips",&trk_nstrips,"trk_nstrips/I",32000);
00661   tree->Branch("slc_nstrips",&slc_nstrips,"slc_nstrips/I",32000);
00662 
00663   tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
00664   tree->Branch("reco_Q2",&reco_Q2,"reco_Q2/F",32000);
00665   tree->Branch("reco_x",&reco_x,"reco_x/F",32000);
00666   tree->Branch("reco_W2",&reco_W2,"reco_W2/F",32000);
00667 
00668   tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
00669   tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);
00670 
00671   tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
00672   tree->Branch("trklength",&trklength,"trklength/F",32000);
00673   tree->Branch("shwlength",&shwlength,"shwlength/F",32000);
00674   tree->Branch("ntrklike",&ntrklike,"ntrklike/I",32000);
00675   tree->Branch("trkend",&trkend,"trkend/I",32000);
00676   tree->Branch("trkendu",&trkendu,"trkendu/I",32000);
00677   tree->Branch("trkendv",&trkendv,"trkendv/I",32000);
00678   tree->Branch("shwbeg",&shwbeg,"shwbeg/I",32000);
00679   tree->Branch("shwend",&shwend,"shwend/I",32000);
00680 
00681 
00682   Int_t trknstp;
00683   tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
00684   tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
00685   tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
00686   tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
00687   tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
00688   tree->Branch("trkchi2",&trkchi2,"trkchi2/F",32000);
00689   tree->Branch("trkndf",&trkndf,"trkndf/I",32000);
00690   tree->Branch("trknstp",&trknstp,"trknstp/I",32000);
00691 
00692 
00693   // vertex track back variables (described above and in code)
00694   tree->Branch("nvsi", &nvsi, "nvsi/I");
00695   tree->Branch("nvsevt", &nvsevt, "nvsevt/I");
00696   tree->Branch("nvsd", &nvsd, "nvsd/F");
00697   tree->Branch("nvst", &nvst, "nvst/F");
00698   tree->Branch("nvsp", &nvsp, "nvsp/F");
00699   tree->Branch("nvti", &nvti, "nvti/I");
00700   tree->Branch("nvtevt", &nvtevt, "nvtevt/I");
00701   tree->Branch("nvtd", &nvtd, "nvtd/F");
00702   tree->Branch("nvtt", &nvtt, "nvtt/F");
00703   tree->Branch("nvtp", &nvtp, "nvtp/F");
00704   /*
00705   tree->Branch("bvsi", &bvsi, "bvsi/I");
00706   tree->Branch("bvsevt", &bvsevt, "bvsevt/I");
00707   tree->Branch("bvsd", &bvsd, "bvsd/F");
00708   tree->Branch("bvst", &bvst, "bvst/F");
00709   tree->Branch("bvsp", &bvsp, "bvsp/F");
00710   tree->Branch("bvti", &bvti, "bvti/I");
00711   tree->Branch("bvtevt", &bvtevt, "bvtevt/I");
00712   tree->Branch("bvtd", &bvtd, "bvtd/F");
00713   tree->Branch("bvtt", &bvtt, "bvtt/F");
00714   tree->Branch("bvtp", &bvtp, "bvtp/F");
00715   */
00716 
00717   tree->Branch("evtph", &evtph, "evtph/F");
00718   tree->Branch("evtphgev", &evtphgev, "evtphgev/F");
00719   
00720   tree->Branch("evtsigfull", &evtsigfull, "evtsigfull/F");
00721   tree->Branch("evtsigpart", &evtsigpart, "evtsigpart/F");
00722 
00723   tree->Branch("trksigfull", &trksigfull, "trksigfull/F");
00724   tree->Branch("trksigpart", &trksigpart, "trksigpart/F");
00725 
00726   tree->Branch("evttimemin",&evttimemin,"evttimemin/D");
00727   tree->Branch("evttimemax",&evttimemax,"evttimemax/D");
00728   tree->Branch("evttimeln",&evttimeln,"evttimeln/D");
00729 
00730   tree->Branch("evttimeminphc",&evttimeminphc,"evttimeminphc/D");
00731   tree->Branch("evttimemaxphc",&evttimemaxphc,"evttimemaxphc/D");
00732   tree->Branch("evttimelnphc",&evttimelnphc,"evttimelnphc/D");
00733 
00734   tree->Branch("earlywph",&earlywph,"earlywph/F");
00735   tree->Branch("earlywphpw",&earlywphpw,"earlywphpw/F");
00736   tree->Branch("earlywphsphere",&earlywphsphere,"earlywphsphere/F");
00737   tree->Branch("ph1000",&ph1000,"ph1000/F");
00738   tree->Branch("ph500",&ph500,"ph500/F");
00739   tree->Branch("ph100",&ph100,"ph100/F");
00740   tree->Branch("lateph1000",&lateph1000,"lateph1000/F");
00741   tree->Branch("lateph500",&lateph500,"lateph500/F");
00742   tree->Branch("lateph100",&lateph100,"lateph100/F");
00743   tree->Branch("lateph50",&lateph50,"lateph50/F");
00744 
00745   tree->Branch("nduprs",&nduprs,"ndubrs/I");
00746   tree->Branch("dupphrs",&dupphrs,"dupphrs/F");
00747 
00748   tree->Branch("nss", nss, "nss[6]/I");
00749   tree->Branch("ess", ess, "ess[6]/F");
00750   tree->Branch("ntotss", &ntotss, "ntotss/I");
00751 
00752   Float_t etu, etv, elu, elv;
00753   tree->Branch("etu", &etu, "etu/F"); // transveres u and v energy
00754   tree->Branch("etv", &etv, "etv/F"); // uses subshowers
00755   tree->Branch("elu", &elu, "elu/F"); // longitudinal energy
00756   tree->Branch("elv", &elv, "elv/F");
00757 
00758   Float_t swu,swv,wswu,wswv; 
00759   tree->Branch("swu",&swu, "swu/F");// shower width in u
00760   tree->Branch("swv",&swv, "swv/F");// shower width in u
00761   tree->Branch("wswu",&wswu, "wswu/F");// shower width in u
00762   tree->Branch("wswv",&wswv, "wswv/F");// shower width in u
00763 
00764   // beam energy weights
00765   // 1/E flux, le, z=050, z=100, z=200, z=250 cm
00766   Float_t bec_inve, bec_le, bec_050, bec_100, bec_200, bec_250;
00767   
00768   tree->Branch("bec_inve", &bec_inve, "bec_inve/F");
00769   tree->Branch("bec_le", &bec_le, "bec_le/F");
00770   tree->Branch("bec_050", &bec_050, "bec_050/F");
00771   tree->Branch("bec_100", &bec_100, "bec_100/F");
00772   tree->Branch("bec_200", &bec_200, "bec_200/F");
00773   tree->Branch("bec_250", &bec_250, "bec_250/F");
00774 
00775   //li time
00776   tree->Branch("litime",&litime,"litime/D");
00777 
00778   //demux status
00779   tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b");
00780   tree->Branch("lowphrat",&lowphrat,"lowphrat/F");
00781 
00782   // beam monitoring variables
00783   Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill;
00784   //goodbeam==1 if spill meets criteria defined below 
00785   //(see the line where goodbeam gets set)
00786   //beamcnf is a unsigned int as defined in BeamMonSpill::StatusBits
00787   UInt_t beamcnf;
00788   //leaving nuTarZ in until we've phased out the summary ntuples
00789   Double_t nuTarZ;
00790   //goodbeam==0 if spill doesnt meet these criteria
00791   Int_t goodbeam;
00792 
00793   //these next variables weren't available in the old beam monitoring ntuples
00794   UInt_t horn_on, target_in, n_batches;
00795   Float_t tor101, tr101d, tortgt, trtgtd;
00796   Float_t hadmon, mumon1, mumon2, mumon3;
00797 
00798   tree->Branch("hbw",&hbw,"hbw/D");
00799   tree->Branch("vbw",&vbw,"vbw/D");
00800   tree->Branch("hpos1",&hpos1,"hpos1/D");
00801   tree->Branch("vpos1",&vpos1,"vpos1/D");
00802   tree->Branch("hpos2",&hpos2,"hpos2/D");
00803   tree->Branch("vpos2",&vpos2,"vpos2/D");
00804   tree->Branch("hornI",&hornI,"hornI/D");
00805   tree->Branch("closestspill",&closestspill,"closestspill/D");
00806   tree->Branch("beamcnf",&beamcnf,"beamcnf/i");
00807   tree->Branch("nuTarZ",&nuTarZ,"nuTarZ/D");
00808   tree->Branch("goodbeam",&goodbeam,"goodbeam/I");
00809   tree->Branch("horn_on",&horn_on,"horn_on/i");
00810   tree->Branch("target_in",&target_in,"target_in/i");
00811   tree->Branch("n_batches",&n_batches,"n_batches/i");
00812   tree->Branch("tor101",&tor101,"tor101/F");
00813   tree->Branch("tr101d",&tr101d,"tr101d/F");
00814   tree->Branch("tortgt",&tortgt,"tortgt/F");
00815   tree->Branch("trtgtd",&trtgtd,"trtgtd/F");
00816   tree->Branch("hadmon",&hadmon,"hadmon/F");
00817   tree->Branch("mumon1",&mumon1,"mumon1/F");
00818   tree->Branch("mumon2",&mumon2,"mumon2/F");
00819   tree->Branch("mumon3",&mumon3,"mumon3/F");
00820   
00821   //coil status
00822   Short_t coil_stat;
00823   tree->Branch("coil_stat",&coil_stat,"coil_stat/S");
00824   
00825 
00826   //for bmpt reweighting
00827   tree->Branch("NuParent","NuParent",&nuparent,8000,1);
00828 
00829   //gnumiflux info
00830   Int_t      fluxrun;     // gnumi flux run number
00831   Int_t      fluxevtno;   // gnumi flux event number
00832   Float_t    nudxdz;       // dx/dz slope, neutrino rndm decay
00833   Float_t    nudydz;       // dy/dz slope, neutrino rndm decay
00834   Float_t    nupz;         // Pz of neutrino (GeV)  rndm decay
00835   Float_t    nuenergy;     // E(neutrino)    (GeV)  rndm decay
00836   Float_t    nudxdznear;   // dx/dz slope, neutrino NearDet center
00837   Float_t    nudydznear;   // dy/dz slope, neutrino NearDet center
00838   Float_t    nuenergynear;    // E(neutrino)    (GeV)  NearDet center
00839   Float_t    nwtnear;     // Weight of nu          NearDet center
00840   Float_t    nudxdzfar;    // dx/dz slope, neutrino FarDet center
00841   Float_t    nudydzfar;    // dy/dz slope, neutrino FarDet center
00842   Float_t    nuenergyfar;    // E(neutrino)    (GeV)  FarDet center
00843   Float_t    nwtfar;      // Weight of nu          FarDet center
00844   Int_t      norig;       // (ignore)
00845   Int_t      ndecay;      // Tag of decay mode
00846   Int_t      ntype;       // Neutrino type (translated to PDG)
00847   Float_t    vx;          // x vertex of hadron (cm)
00848   Float_t    vy;          // y vertex of hadron (cm)
00849   Float_t    vz;          // z vertex of hadron (cm)
00850   Float_t    pdpx;        // nu parent px at decay point
00851   Float_t    pdpy;        // nu parent py at decay point
00852   Float_t    pdpz;        // nu parent pz at decay point
00853   Float_t    ppdxdz;      // nu parent slope at production
00854   Float_t    ppdydz;      // nu parent slope at production
00855   Float_t    pppz;        // nu parent pz at production
00856   Float_t    ppenergy;    // nu parent energy at production
00857   Int_t      ppmedium;    // GEANT medium of nu parent at parent production
00858   Int_t      ptype;       // nu parent type (translated to PDG)
00859   Float_t    ppvx;        // nu parent production vtx x
00860   Float_t    ppvy;        // nu parent production vtx y
00861   Float_t    ppvz;        // nu parent production vtx z
00862   Float_t    muparpx;     // if parent=mu, hadron parent px
00863   Float_t    muparpy;     // if parent=mu, hadron parent py
00864   Float_t    muparpz;     // if parent=mu, hadron parent pz
00865   Float_t    mupare;      // if parent=mu, hadron parent energy
00866   Float_t    necm;        // E(nu) in parent cm
00867   Float_t    nimpwt;      // importance weight
00868   Float_t    xpoint;      // (unused)
00869   Float_t    ypoint;      // (unused)
00870   Float_t    zpoint;      // (unused)
00871   Float_t    tvx;         // target exit point (x) of parent
00872   Float_t    tvy;         // target exit point (y) of parent
00873   Float_t    tvz;         // target exit point (z) of parent
00874   Float_t    tpx;         // parent px at target exit
00875   Float_t    tpy;         // parent py at target exit
00876   Float_t    tpz;         // parent pz at target exit
00877   Int_t      tptype;      // parent particle type (translated to PDG)
00878   Int_t      tgen;        // parent generation in cascade
00879 
00880   tree->Branch("fluxrun",&fluxrun,"fluxrun/I");
00881   tree->Branch("fluxevtno",&fluxevtno,"fluxevtno/I");
00882   tree->Branch("nudxdz",&nudxdz,"nudxdz/F");
00883   tree->Branch("nudydz",&nudydz,"nudydz/F");
00884   tree->Branch("nupz",&nupz,"nupz/F");
00885   tree->Branch("nuenergy",&nuenergy,"nuenergy/F");
00886   tree->Branch("nudxdznear",&nudxdznear,"nudxdznear/F");
00887   tree->Branch("nudydznear",&nudydznear,"nudydznear/F");
00888   tree->Branch("nuenergynear",&nuenergynear,"nuenergynear/F");
00889   tree->Branch("nwtnear",&nwtnear,"nwtnear/F");
00890   tree->Branch("nudxdzfar",&nudxdzfar,"nudxdzfar/F");
00891   tree->Branch("nudydzfar",&nudydzfar,"nudydzfar/F");
00892   tree->Branch("nuenergyfar",&nuenergyfar,"nuenergyfar/F");
00893   tree->Branch("nwtfar",&nwtfar,"nwtfar/F");
00894   tree->Branch("norig",&norig,"norig/I");
00895   tree->Branch("ndecay",&ndecay,"ndecay/I");
00896   tree->Branch("ntype",&ntype,"ntype/I");
00897   tree->Branch("vx",&vx,"vx/F");
00898   tree->Branch("vy",&vy,"vy/F");
00899   tree->Branch("vz",&vz,"vz/F");
00900   tree->Branch("pdpx",&pdpx,"pdpx/F");
00901   tree->Branch("pdpy",&pdpy,"pdpy/F");
00902   tree->Branch("pdpz",&pdpz,"pdpz/F");
00903   tree->Branch("ppdxdz",&ppdxdz,"ppdxdz/F");
00904   tree->Branch("ppdydz",&ppdydz,"ppdydz/F");
00905   tree->Branch("pppz",&pppz,"pppz/F");
00906   tree->Branch("ppenergy",&ppenergy,"ppenergy/F");
00907   tree->Branch("ppmedium",&ppmedium,"ppmedium/I");
00908   tree->Branch("ptype",&ptype,"ptype/I");
00909   tree->Branch("ppvx",&ppvx,"ppvx/F");
00910   tree->Branch("ppvy",&ppvy,"ppvy/F");
00911   tree->Branch("ppvz",&ppvz,"ppvz/F");
00912   tree->Branch("muparpx",&muparpx,"muparpx/F");
00913   tree->Branch("muparpy",&muparpy,"muparpy/F");
00914   tree->Branch("muparpz",&muparpz,"muparpz/F");
00915   tree->Branch("mupare",&mupare,"mupare/F");
00916   tree->Branch("necm",&necm,"necm/F");
00917   tree->Branch("nimpwt",&nimpwt,"nimpwt/F");
00918   tree->Branch("xpoint",&xpoint,"xpoint/F");
00919   tree->Branch("ypoint",&ypoint,"ypoint/F");
00920   tree->Branch("zpoint",&zpoint,"zpoint/F");
00921   tree->Branch("tvx",&tvx,"tvx/F");
00922   tree->Branch("tvy",&tvy,"tvy/F");
00923   tree->Branch("tvz",&tvz,"tvz/F");
00924   tree->Branch("tpx",&tpx,"tpx/F");
00925   tree->Branch("tpy",&tpy,"tpy/F");
00926   tree->Branch("tpz",&tpz,"tpz/F");
00927   tree->Branch("tptype",&tptype,"tptype/I");
00928   tree->Branch("tgen",&tgen,"tgen/I");
00929   
00930   //pot counting
00931   float pottor101=0.;
00932   float pottortgt=0.;
00933   float pot=0.;
00934   int NRUNS=0;
00935   int NSNARLS=0;
00936   TTree *pottree = new TTree("pottree","pottree");
00937   pottree->Branch("pot",&pot,"pot/F");
00938   pottree->Branch("pottor101",&pottor101,"pottor101/F");
00939   pottree->Branch("pottortgt",&pottortgt,"pottortgt/F");
00940   pottree->Branch("nruns",&NRUNS,"nruns/I");
00941   pottree->Branch("nsnarls",&NSNARLS,"nsnarls/I");
00942 
00943   //make a BMSpillAna so we can tell if it's goodbeam
00944   BMSpillAna bmana;
00945 
00946 
00950   TTree* inukerw = 0;
00951   int n_inuke_wts=0;
00952   float* inuke_wts=0;
00953   double nucrad=-1.0;
00954   double wrad=1.0;
00955   std::vector<inuke_reweight::parameter_set> psets;
00956   
00957   if(file_is_mc && fDoInukeReweight){
00958     
00959     inuke_reweight::parameter_set sigmas;
00960 
00961     neugen_config nc("MODBYRS");
00962     nc.set_config_no(4);
00963     neugen_wrapper nw(&nc);
00964     nw.print_configuration();
00965     gSystem->Setenv("NEUGEN3_XSEC","MODBYRS4_xsec.dat");
00966     nw.begin_generation();
00967 
00968     /*    
00969     sigmas.pi_fate.inel=0.2;
00970     sigmas.pi_fate.abs=0.25;
00971     sigmas.pi_scale.ft=0.5;
00972     sigmas.pi_scale.xsec=0.1;    
00973     //    sigmas.pn_scale.ft=0.5;
00974     sigmas.pn_scale.xsec=0.1;
00975     */
00976 
00977     sigmas.pi_fate.elas=0.1;
00978     sigmas.pi_fate.inel=0.4;
00979     sigmas.pi_fate.abs=0.3;
00980     sigmas.pi_fate.cex=0.49;
00981     sigmas.pi_fate.piprod=0.2;
00982 
00983     sigmas.pi_scale.xsec=0.1;
00984 
00985     sigmas.pn_fate.abs=0.2;
00986     sigmas.pn_fate.piprod=0.2;
00987 
00988     sigmas.pn_scale.xsec=0.1;
00989 
00990     sigmas.pi_scale.ft=0.5; // apply to both pi and p
00991 
00992     std::cout<<"Inuke Reweight sigmas: "<<std::endl;
00993     sigmas.print(std::cout);
00994 
00995     if(fNRandomSets<=0) {
00996       inuke_reweight::generate_1sigma_shifts(sigmas,psets);    
00997       for (unsigned int iset=0; iset<psets.size(); iset++){
00998         psets[iset].pi_fate.set_abs_fates();
00999         psets[iset].pn_fate.set_abs_fates();
01000         psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime
01001         // correlate change in pn xsec with change in pn elastic fate
01002         if(psets[iset].pn_scale.xsec!=0) {
01003           psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec;
01004         }
01005       }
01006     }
01007     else{
01008       inuke_reweight::parameter_limits limits(sigmas,2.0);
01009       limits.upper.pi_scale.ft=0.5;  // special case, 1 sigma, but uniform dist
01010       limits.lower.pi_scale.ft=-0.5;
01011       std::cout<<"Inuke Reweight Limits: "<<std::endl;
01012       limits.print(std::cout);
01013       gRandom->SetSeed(1);
01014       for(int i=0; i<fNRandomSets; i++){
01015         inuke_reweight::parameter_set p;
01016         inuke_reweight::generate_uncor_shifts(sigmas,limits,p);
01017         psets.push_back(p);
01018       }
01019       // massage parameters
01020       for(int iset=0; iset<fNRandomSets; iset++){
01021         psets[iset].pi_fate.set_abs_fates();
01022         psets[iset].pn_fate.set_abs_fates();
01023         psets[iset].pn_scale.ft=psets[iset].pi_scale.ft; // pi,pn same Ftime
01024         // correlate change in pn xsec with change in pn elastic fate
01025         if(psets[iset].pn_scale.xsec!=0) {
01026           psets[iset].pn_fate.elas = 2.0*psets[iset].pn_scale.xsec;
01027         }
01028       }
01029 
01030 
01031     }
01032 
01033     n_inuke_wts=psets.size();
01034     inuke_wts = new float[n_inuke_wts];
01035     
01036     //    inukerw = new TTree("inukerw","intranuke weights");
01037     inukerw=tree;
01038     inukerw->Branch("r",&run,"r/I");
01039     inukerw->Branch("s",&snarl,"s/I");
01040     inukerw->Branch("e",&event,"e/I");
01041     inukerw->Branch("nwts",&n_inuke_wts,"nwts/I");
01042     inukerw->Branch("wts",inuke_wts,"wts[nwts]/F");          
01043     inukerw->Branch("nucrad",&nucrad,"nucrad/D");
01044     inukerw->Branch("wrad",&wrad,"wrad/D");
01045 
01046     // save psets
01047     // 9*2 fates + 3*2 scales = 24
01048     for(int ipar=0; ipar<24; ipar++){
01049       TVectorD* vptr = new TVectorD(n_inuke_wts);
01050       inukerw->GetUserInfo()->Add(vptr);
01051     }
01052     for(int iset=0; iset<n_inuke_wts; iset++){
01053       for(int ipar=0; ipar<24; ipar++){ 
01054         TVectorD* vptr = 
01055           dynamic_cast<TVectorD*>(inukerw->GetUserInfo()->At(ipar));    
01056         (*vptr)[iset]=psets[iset].get_par(ipar);
01057       }      
01058     }
01059 
01060 
01061   }
01062 
01066 
01067 
01068   //  map<VldTimeStamp, BeamMonTV>* beam_mon_map=0;
01069   //  if(!file_is_mc) beam_mon_map = new BeamMonMap::MakeBeamMonMap(bmonpath);
01070   //  map<VldTimeStamp, BeamMonTV> beam_mon_map
01071   //    =BeamMonMap::MakeBeamMonMap(bmonpath);
01072   //  std::cout<<"map file: "<<bmonpath<<std::endl;
01073   //  std::cout<<"map size: "<<beam_mon_map.size()<<std::endl;
01074 
01080   BeamEnergyCalculator ecalc(&fBECConfig);
01081   double bec_high=0; double bec_low=0;
01082   ecalc.GetStandardConfig()->Print();
01083   ecalc.GetStandardConfig()->Get("beam:high_energy_limit",bec_high);
01084   ecalc.GetStandardConfig()->Get("beam:low_energy_limit",bec_low);
01085 
01086   // mark's QE id
01087   //  MadQEID qeid;
01088   //  qeid.ReadPDFs("/afs/fnal.gov/files/data/minos/d04/kordosky/madqe_files/QEpdfs.root");
01089 
01090   
01092   //  bool beamdb=false;
01093   //  bool checkeddb=false;
01094   int lastrun=-1;
01095   for(int ii=0;ii<Nentries;ii++){   
01096     if(ii%1000==0) std::cout << float(ii)*100./float(Nentries) 
01097                             << "% done" << std::endl;
01098 
01099     if(!this->GetEntry(ii)) continue;
01100     NSNARLS++;
01101     // fill some histograms for mc acceptance
01102     if(file_is_mc) FillMCHists(&save);
01103 
01104     /*
01105     if(!checkeddb){
01106       //      ntpHeader->GetVldContext().Print();
01107       //      cout<<"IM HERE!!!!"<<endl;
01108       DbiResultPtr<BeamMonSpill> testbs(ntpHeader->GetVldContext(), 
01109                                         Dbi::kDefaultTask, 
01110                                         Dbi::kDisabled);
01111       //      cout<<"tried to get the dbiresult ptr"<<endl;
01112       checkeddb=true;
01113       //      cout<<"****************checking for beam monitoring db:"<<endl;
01114       if(testbs.GetNumRows()>0){
01115         beamdb=true;
01116         cout<<"Found the Beam Monitoring Database Tables"<<endl;
01117       }
01118       else{     
01119         cout<<"Didn't find the Beam Monitoring Database Tables"<<endl;
01120         cout<<"Resorting to old beamsummary ntuples"<<endl;
01121       }
01122     }
01123     */
01124 
01125     nevent = eventSummary->nevent;
01126     run = ntpHeader->GetRun();
01127     if(run!=lastrun){
01128       NRUNS++;
01129       lastrun = run;
01130     }
01131     snarl = ntpHeader->GetSnarl();
01132     subrun = ntpHeader->GetSubRun();
01133     trgsrc = ntpHeader->GetTrigSrc();
01134     spilltype = ntpHeader->GetRemoteSpillType();
01135     trgtime = eventSummary->trigtime;
01136     litime = eventSummary->litime;
01137     coil_stat = detStatus->coilstatus;
01138     dmx_stat=0;
01139 
01140     //figure out demux status
01141     dmx_stat+=(dmxStatus->nonphysicalfail);
01142     dmx_stat+=((dmxStatus->validplanesfail<<1));
01143     dmx_stat+=((dmxStatus->vertexplanefail<<2));
01144 
01145     //figure out how much of snarl ph was deposited as low ph
01146     lowphrat=0.;
01147     lowphrat = ComputeLowPHRatio();
01148 
01150     // look through "events" and characterize how close together
01151     // their vertices are
01152     // idea is to flag runt events caused by late activity near vertex
01154     NearbyVertexFinder nby(nevent);
01155     nby.ProcessEntry(this);
01156 
01157 
01158     detector = -1; 
01159     //get detector type:
01160     const Detector::Detector_t det 
01161       = ntpHeader->GetVldContext().GetDetector();
01162     if(det==Detector::kFar){
01163       detector = 2;
01164     }
01165     else if(det==Detector::kNear){
01166       detector = 1;     
01167     }
01168     // deal with beam type as best we can
01169     BeamType::BeamType_t current_beam=BeamType::kUnknown;
01170     if(file_is_mc){
01171       // for MC set beam to fMCbeam
01172       current_beam=fMCBeam;
01173       // for data we will try to deduce from beam monitoring
01174     }
01175 
01176 
01178     // first thing, if data do beam monitoring
01180     hbw=vbw=hpos1=vpos1=hpos2=vpos2=hornI=closestspill=0.0;
01181     beamcnf=0;
01182     nuTarZ=0.0;
01183     goodbeam=0;
01184     horn_on=target_in=n_batches=0;
01185     tor101=tr101d=tortgt=trtgtd=0.0;
01186     hadmon=mumon1=mumon2=mumon3=0.0;  
01187 
01188     if(!file_is_mc){
01189       VldContext vc = ntpHeader->GetVldContext();
01190       VldTimeStamp vts = vc.GetTimeStamp();
01191       //      if(beamdb){
01192       BDSpillAccessor &sa = BDSpillAccessor::Get();
01193       const BeamMonSpill *bs = sa.LoadSpill(vts);
01194       hbw=bs->fProfWidX*1000.;//make it in mm
01195       vbw=bs->fProfWidY*1000.;//make it in mm
01196       hpos1=bs->fTargProfX*1000.;//make it in mm
01197       vpos1=bs->fTargProfY*1000.;//make it in mm
01198       n_batches=0;
01199       float btint=0.;
01200       hpos2=0.;
01201       vpos2=0.;
01202       for(int bt=0;bt<6;bt++){
01203         //cout<<"batch "<<bt<<" x "<<bs->fTargBpmX[bt]<<" y "<<bs->fTargBpmY[bt]
01204         //              <<" int "<<bs->fBpmInt[bt]<<endl;
01205         hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
01206         vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);
01207         if(bs->fBpmInt[bt]>0.001){
01208           n_batches++;
01209         }
01210         btint+=bs->fBpmInt[bt];
01211       }
01212       if(btint!=0){
01213         hpos2=hpos2*1000./btint;//make it in mm
01214         vpos2=vpos2*1000./btint;//make it in mm
01215       }
01216       else{
01217         hpos2=-999.;
01218         vpos2=-999.;
01219       }
01220       
01221       hornI=bs->fHornCur;       
01222       SpillTimeFinder &sfind= SpillTimeFinder::Instance();
01223       closestspill = 
01224         sfind.GetTimeToNearestSpill(vc);
01225       
01226       beamcnf =bs->GetStatusBits().beam_type;
01227       nuTarZ=0.;
01228       horn_on = bs->GetStatusBits().horn_on;
01229       target_in = bs->GetStatusBits().target_in;
01230       
01231       
01232       tor101=bs->fTor101;
01233       tr101d=bs->fTr101d;
01234       tortgt=bs->fTortgt;
01235       trtgtd=bs->fTrtgtd;
01236       hadmon=bs->fHadInt;
01237       mumon1=bs->fMuInt1;
01238       mumon2=bs->fMuInt2;
01239       mumon3=bs->fMuInt3;
01240       
01241       goodbeam=0;
01242       //      }
01243       /*
01244       else{
01245         // could the following return a reference?
01246         BeamMonTV btv = BeamMonMap::FindClosestSpill(beam_mon_map,vts);
01247         hbw=btv.hbw;
01248         vbw=btv.vbw;
01249         hpos2=btv.hpos2;
01250         vpos2=btv.vpos2;
01251         hornI=btv.hornI;
01252         closestspill=btv.closestspill;
01253         
01254         nuTarZ=btv.nuTarZ;
01255         //approximately assign beamcnf based on nuTarZ
01256         if(nuTarZ<20000&&nuTarZ>=0){
01257           beamcnf=1;
01258           target_in=1;
01259         }
01260         else if(nuTarZ>20000&&nuTarZ<60000){
01261           beamcnf=4;
01262           target_in=1;
01263         }
01264         else{
01265           beamcnf=5;
01266           target_in=1;
01267         }
01268         if(fabs(hornI)>10){
01269           horn_on=1;
01270         }
01271         n_batches=0;
01272         tor101=0.;
01273         tr101d=0.;        
01274         tortgt=btv.bI;
01275         trtgtd=0.;
01276         hadmon=0.;
01277         mumon1=0.;
01278         mumon2=0.;
01279         mumon3=0.;
01280       }
01281       */
01282       //set goodbeam==1 if it passes some basic beam monitoring cuts
01283       //changing widths to deal with rms instead of fits
01284       //getting rid of horn_on and target_in since they sometime
01285       //go haywire
01286       //        if(horn_on&&target_in&&tortgt>0.1&&hbw<1.5&&vbw<2&&
01287       //           hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(closestspill)<2){
01288       /*
01289       if(tor101>0.1&&tortgt>0.1&&hbw<2.9&&vbw<2.9&&
01290          hpos2>-2&&hpos2<0&&vpos2>0&&vpos2<2&&fabs(hornI)>50&&
01291          fabs(closestspill)<2){
01292         goodbeam=1;
01293       }
01294       */
01295       //use Mark D. BMSpillAna to set good beam
01296       bmana.SetSpill(*bs);
01297       bmana.SetTimeDiff(closestspill);
01298       if(bmana.SelectSpill()){
01299         goodbeam=1;
01300       }
01301       else{
01302         goodbeam=0;
01303       }
01304 
01305 
01306       if(fDataUseMCBeam==true){
01307         //allow us to force the data to use the MC beam
01308         current_beam=fMCBeam;
01309       }
01310       else{
01311         // attempt to deduce beam type from beam monitoring
01312         current_beam = BeamType::FromBeamMon(beamcnf);
01313       }
01314       
01315       if(goodbeam==1&&spilltype!=3){//don't count fake spills
01316         pot+=tortgt;
01317         pottor101+=tor101;
01318         pottortgt+=tortgt;
01319       }
01320     }
01321 
01322     // Write a special record for each spill!
01323     event=-1; ntrack=-1; nshower=-1; reco_emu=0.0; reco_enu=0.0; reco_eshw=0;
01324     true_enu=0.0; true_emu=0.0; true_eshw=0.0;
01325     pass=0; is_fid=0;
01326     tree->Fill();
01327     for(int iwt=0; iwt<n_inuke_wts; iwt++) inuke_wts[iwt]=-1.0;
01328     if(inukerw!=0 && inukerw!=tree) inukerw->Fill();
01329 
01330     if(nevent==0) continue;    
01331     //fill tree once for each reconstructed event
01332     for(int i=0;i<eventSummary->nevent;i++){ 
01333       if(!LoadEvent(i)) continue; //no event found
01334 
01335       //set event index:
01336       event = i;
01337       ntrack = ntpEvent->ntrack;
01338       nshower = ntpEvent->nshower;
01339       nstrip = ntpEvent->nstrip;
01340 
01342       //zero all tree values
01343       true_enu = 0; true_emu = 0; true_eshw = 0; true_x = 0; true_y = 0; 
01344       true_q2 = 0; true_w2 = 0; true_dircosneu = 0; true_dircosz = 0;
01345       true_vtxx = 0; true_vtxy = 0; true_vtxz = 0; true_pitt_fid=0;
01346       true_trk_cmplt=0.;
01347       resnum=0;
01348       npp=npm=npz=np=nn=0;
01349       epp=epm=epz=ep=en=0.0;
01350 
01351       fluxrun=0; fluxevtno=0; nudxdz=0.; nudydz=0.; nupz=0.; nuenergy=0.;
01352       nudxdznear=0.; nudydznear=0.; nuenergynear=0.; nwtnear=0.; nudxdzfar=0.;
01353       nudydzfar=0.; nuenergyfar=0.; nwtfar=0.; norig=0; ndecay=0; ntype=0;
01354       vx=0.; vy=0.; vz=0.; pdpx=0.; pdpy=0.; pdpz=0.; ppdxdz=0.; ppdydz=0.; pppz=0.;
01355       ppenergy=0.; ppmedium=0; ptype=0; ppvx=0.; ppvy=0.; ppvz=0.;
01356       muparpx=0.; muparpy=0.; muparpz=0.; mupare=0.; necm=0.; nimpwt=0.;
01357       xpoint=0.; ypoint=0.; zpoint=0.; tvx=0.; tvy=0.; tvz=0.; tpx=0.;
01358       tpy=0.; tpz=0.; tptype=0; tgen=0;
01359       
01360       flavour = 0; process = 0; nucleus = 0; cc_nc = 0; initial_state = 0;
01361 
01362       mcevent = -1;
01363       is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; 
01364       is_pitt_fid=0; pitt_evt_class=0; pitt_trk_qual=0; nd_radial_fid=0;
01365       is_trk_fid=0;
01366       duvvtx=0;
01367       pid0 = 0; pid1 = 0; pid2 = 0; pass = 0;
01368       reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_wtd_eshw=0; 
01369       reco_vtx_eshw=0; reco_vtx_wtd_eshw=0;
01370       reco_qe_enu = 0;  reco_mushw=0;
01371       reco_qe_q2 = 0; reco_dircosneu = 0; reco_dircosz = 0;
01372       reco_vtxx = reco_vtxy = reco_vtxz = reco_vtxr=0; 
01373       reco_trk_vtxx = reco_trk_vtxy = reco_trk_vtxz = reco_trk_vtxr=0; 
01374       reco_shw_vtxx = reco_shw_vtxy = reco_shw_vtxz = reco_shw_vtxr=0; 
01375       reco_trk_endx = reco_trk_endy = reco_trk_endz = reco_trk_endr=0; 
01376       evt_nstrips = slc_nstrips = trk_nstrips=0;
01377       evt_nplanes = slc_nplanes = trk_nplanes=0;
01378       shw_nplanes_cut=shw_nplanes_raw=0;
01379       shw_ph_cut=shw_ph_raw=0.0;
01380 
01381       evtlength = trklength = shwlength=trknstp=0; 
01382       trkphfrac = trkphplane = 0;
01383       ntrklike= trkend=trkendu=trkendv=shwbeg=shwend=0;
01384       trkmomrange = 0; trkqp = 0; trkeqp = 0;
01385       trkchi2=0.0; trkndf=0;
01386       pass_track=0; emu_meth=0;
01387       // kinematic quantities always positive
01388       // convention: -1 indicates unfilled variable (NC events)
01389       reco_y = reco_x = reco_Q2 = reco_W2 = -1;
01390       
01391       shw_pur = shw_comp = -1;
01392 
01393       // vertex back tracking
01394 
01395       nvsi=nvti=bvsi=bvti=-1;
01396       nvsevt=bvsevt=nvtevt=bvtevt=-1;
01397       nvsd=nvtd=bvsd=bvtd=999999;
01398       nvst=bvst=nvtt=bvtt=1e6;
01399       nvsp=nvtp=bvsp=bvtp=-1.0;
01400 
01401       evtph=0; evtphgev=0;
01402       evtsigfull=evtsigpart=trksigfull=trksigpart=0;            
01403 
01404       evttimemin=evttimemax=evttimeln=0.;
01405       evttimeminphc=evttimemaxphc=evttimelnphc=0.;
01406 
01407       earlywph=earlywphpw=earlywphsphere=0.;
01408       ph1000=ph500=ph100=0.;
01409       lateph1000=lateph500=lateph100=lateph50=0.;
01410 
01411       nduprs=0;
01412       dupphrs=0.;
01413       
01414       nu_px=nu_py=nu_pz=tar_px=tar_py=tar_pz=tar_e=0.0;
01415       inu=0;
01416       ntotss=0;
01417       for(int iss=0; iss<6; iss++){nss[iss]=0; ess[iss]=0;}
01418 
01419       med_qe_pid=-1001.0;
01420       med_qe_pass=0;
01421       niki_cc_pid=-999;
01422       dave_cc_pid=-999;
01423 
01424       is_fid = InFidVol();
01425             
01426       //check for a corresponding mc event
01427       // first, looks to match reconstructed track
01428       // then, looks to match reconstructed shower
01429       if(LoadTruthForRecoTH(i,mcevent)) {
01430 
01431         // should flag with a "bad truth word"
01432         // in case the function above returns false
01433 
01434         //true info for tree:
01435         is_mc             = 1;
01436         nucleus           = Nucleus(mcevent);
01437         flavour           = Flavour(mcevent);
01438         initial_state     = Initial_state(mcevent);
01439         process           = IResonance(mcevent); 
01440         cc_nc             = IAction(mcevent);
01441         true_enu          = TrueNuEnergy(mcevent); 
01442         true_emu          = TrueMuEnergy(mcevent); 
01443         true_eshw         = TrueShwEnergy(mcevent); 
01444         true_x            = X(mcevent);
01445         true_y            = Y(mcevent);
01446         true_q2           = Q2(mcevent);
01447         true_w2           = W2(mcevent);
01448         true_dircosz      = TrueMuDCosZ(mcevent);
01449         true_dircosneu    = TrueMuDCosNeu(mcevent);
01450         
01451 
01452         //flux info
01453         if(LoadTruth(mcevent)){
01454           fluxrun=ntpTruth->flux.fluxrun;
01455           fluxevtno=ntpTruth->flux.fluxevtno;
01456           nudxdz=ntpTruth->flux.ndxdz;
01457           nudydz=ntpTruth->flux.ndydz;
01458           nupz=ntpTruth->flux.npz;
01459           nuenergy=ntpTruth->flux.nenergy;
01460           nudxdznear=ntpTruth->flux.ndxdznear;
01461           nudydznear=ntpTruth->flux.ndydznear;
01462           nuenergynear=ntpTruth->flux.nenergynear;
01463           nwtnear=ntpTruth->flux.nwtnear;
01464           nudxdzfar=ntpTruth->flux.ndxdzfar;
01465           nudydzfar=ntpTruth->flux.ndydzfar;
01466           nuenergyfar=ntpTruth->flux.nenergyfar;
01467           nwtfar=ntpTruth->flux.nwtfar;
01468           norig=ntpTruth->flux.norig;
01469           ndecay=ntpTruth->flux.ndecay;
01470           ntype=ntpTruth->flux.ntype;
01471           vx=ntpTruth->flux.vx;
01472           vy=ntpTruth->flux.vy;
01473           vz=ntpTruth->flux.vz;
01474           pdpx=ntpTruth->flux.pdpx;
01475           pdpy=ntpTruth->flux.pdpy;
01476           pdpz=ntpTruth->flux.pdpz;
01477           ppdxdz=ntpTruth->flux.ppdxdz;
01478           ppdydz=ntpTruth->flux.ppdydz;
01479           pppz=ntpTruth->flux.pppz;
01480           ppenergy=ntpTruth->flux.ppenergy;
01481           ppmedium=ntpTruth->flux.ppmedium;
01482           ptype=ntpTruth->flux.ptype;
01483           ppvx=ntpTruth->flux.ppvx;
01484           ppvy=ntpTruth->flux.ppvy;
01485           ppvz=ntpTruth->flux.ppvz;
01486           muparpx=ntpTruth->flux.muparpx;
01487           muparpy=ntpTruth->flux.muparpy;
01488           muparpz=ntpTruth->flux.muparpz;
01489           mupare=ntpTruth->flux.mupare;
01490           necm=ntpTruth->flux.necm;
01491           nimpwt=ntpTruth->flux.nimpwt;
01492           xpoint=ntpTruth->flux.xpoint;
01493           ypoint=ntpTruth->flux.ypoint;
01494           zpoint=ntpTruth->flux.zpoint;
01495           tvx=ntpTruth->flux.tvx;
01496           tvy=ntpTruth->flux.tvy;
01497           tvz=ntpTruth->flux.tvz;
01498           tpx=ntpTruth->flux.tpx;
01499           tpy=ntpTruth->flux.tpy;
01500           tpz=ntpTruth->flux.tpz;
01501           tptype=ntpTruth->flux.tptype;
01502           tgen=ntpTruth->flux.tgen;       
01503         }
01504 
01505         if(ntpTHTrack){
01506           if(ntpTHTrack->neumc==mcevent){
01507             true_trk_cmplt=ntpTHTrack->completeslc;
01508           }
01509         }
01510         
01511         Float_t* true_p = TrueNuMom(mcevent);
01512         if(true_p){
01513           nu_px=true_p[0]; nu_py=true_p[1]; nu_pz=true_p[2];
01514           delete[] true_p; true_p=0;
01515         }
01516         true_p = Target4Vector(mcevent);
01517         if(true_p){
01518           tar_px=true_p[0]; tar_py=true_p[1]; tar_pz=true_p[2]; 
01519           tar_e=true_p[3];
01520           delete[] true_p; true_p=0;
01521         }
01522 
01523         Float_t *vtx = TrueVtx(mcevent);
01524         true_vtxx         = vtx[0];
01525         true_vtxy         = vtx[1];
01526         true_vtxz         = vtx[2];
01527         float true_vtxu = (true_vtxx + true_vtxy)/sqrt(2.0);
01528         float true_vtxv = (true_vtxy - true_vtxx)/sqrt(2.0);
01529         true_pitt_fid=PittInFidVol(true_vtxx, true_vtxy, 
01530                                    true_vtxz,
01531                                    true_vtxu, true_vtxv);
01532 
01533         resnum = HadronicFinalState(mcevent);
01534         npp = NumFSPion(mcevent,+1);
01535         npm = NumFSPion(mcevent,-1);
01536         npz = NumFSPiZero(mcevent);
01537         np = NumFSProton(mcevent);
01538         nn = NumFSNeutron(mcevent);
01539         
01540         epp = TotFSPionEnergy(mcevent,+1);
01541         epm = TotFSPionEnergy(mcevent,-1);
01542         epz = TotFSPiZeroEnergy(mcevent);
01543         ep = TotFSProtonEnergy(mcevent);
01544         en = TotFSNeutronEnergy(mcevent);
01545         
01546         // beam energy weights //
01547         inu=0;
01548         inu = INu(mcevent);
01549         /*
01550         bec_inve=bec_le=bec_050=bec_100=bec_200=bec_250=1.0;
01551         bec_inve=ecalc.GetWeight(BeamType::kInverseE,det,inu,
01552                                  true_enu, bec_high,bec_low);
01553         bec_le=ecalc.GetWeight(BeamType::kLE,det,inu,
01554                                true_enu, bec_high,bec_low);
01555         bec_050=ecalc.GetWeight(BeamType::k050,det,inu,
01556                                  true_enu, bec_high,bec_low);
01557         bec_100=ecalc.GetWeight(BeamType::k100,det,inu,
01558                                  true_enu, bec_high,bec_low);
01559         bec_200=ecalc.GetWeight(BeamType::k200,det,inu,
01560                                  true_enu, bec_high,bec_low);
01561         bec_250=ecalc.GetWeight(BeamType::k250,det,inu,
01562                                  true_enu, bec_high,bec_low);
01563         */
01565         delete [] vtx;
01566       }
01567 
01568       if(PassBasicCuts(event)) {
01569         //reco info for tree:
01570         pass_basic        = 1;
01571         
01572         
01573         //index will -1 if no track/shower present in event
01574         int track_index   = -1;
01575         // track spanning largest number of planes
01576         LoadLargestTrackFromEvent(i,track_index);
01577 
01578         pass_track = PassTrackCuts(track_index);
01579         //methods should all be protected against index = -1
01580         
01581         reco_emu          = RecoMKMuEnergy(emu_meth,track_index);
01582 
01583         int shower_index  = -1;
01584         
01585         // load the largest shower
01586         /*
01587         if(track_index!=-1) {
01588           LoadShowerAtTrackVertex(i,track_index,shower_index);
01589         }
01590         // shower with the best match to the vertex     
01591         else LoadLargestShowerFromEvent(i,shower_index);
01592         */
01593         /*
01594         if(nshower>0){
01595           shower_index=ntpEvent->shw[0];
01596           LoadShower(shower_index);
01597         }
01598         */
01599 
01600         // use the jim shower!
01601         LoadShower_Jim(i,shower_index,det);
01602         reco_eshw         = RecoShwEnergy(shower_index, det,1,!file_is_mc);
01603         if(shower_index>-1) reco_wtd_eshw = ntpShower->shwph.wtCCgev;
01604         reco_vtx_eshw = reco_eshw;
01605         reco_vtx_wtd_eshw = reco_wtd_eshw;
01606 
01607         // figure shower truth business out
01608         if(is_mc && shower_index>-1){
01609           LoadTHShower(shower_index);
01610           assert(ntpTHShower);
01611           shw_pur=ntpTHShower->purity;
01612           shw_comp=ntpTHShower->completeslc;
01613         }
01614         
01615         // provisional, to be updated later
01616         reco_enu = reco_emu + reco_eshw;
01617         // event energy according to offline
01618         evtphgev = ntpEvent->ph.gev;
01619         if(reco_enu>0){
01620           reco_qe_enu       = RecoQENuEnergy(track_index);
01621           reco_qe_q2        = RecoQEQ2(track_index);
01622           // note, NtpSRFiducial isn't filled for ND events
01623           // is_cev is only useful for FD events
01624           is_cev            = IsFidAll(track_index);
01625           // use Pitt group's fiducial cuts for ND events
01626           if(detector==1){
01627             // for now, use event vertex
01628             // but maybe should use track vertex...
01629             // if vertex finder works well it shouldn't matter...
01630             is_pitt_fid = PittInFidVol(ntpEvent->vtx.x, ntpEvent->vtx.y, 
01631                                        ntpEvent->vtx.z,
01632                                        ntpEvent->vtx.u, ntpEvent->vtx.v);
01633             nd_radial_fid = NDRadialFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,
01634                                            ntpEvent->vtx.z);
01635           }
01636           else{
01637             is_pitt_fid = InFidVol();
01638           }
01639           // have already checked that ntpEvent exists
01640           reco_vtxx         = ntpEvent->vtx.x;
01641           reco_vtxy         = ntpEvent->vtx.y;
01642           reco_vtxz         = ntpEvent->vtx.z;
01643           
01644           // check if distance from vertex to shower is too large
01645           // if so, set reco_eshw =0 
01646           // idea is to handle case of no vertex shower but runt downstream
01647           if(shower_index!=-1){
01648             //      bool shw_ok=true;
01649             const float max_shw_dist=14*0.06; // about 2 interaction lengths
01650             float dist_z = fabs(reco_vtxz - ntpShower->vtx.z);
01651             float dist_r = sqrt( pow(ntpEvent->vtx.u-ntpShower->vtx.u,2) +
01652                                  pow(ntpEvent->vtx.v-ntpShower->vtx.v,2) );
01653             if( (dist_z + dist_r) > max_shw_dist ){
01654               reco_vtx_eshw=0;
01655               reco_vtx_wtd_eshw=0;
01656             }
01657             
01658           }
01659           // radial position of vertex
01660           if(detector==1){
01661             reco_vtxr = pow(reco_vtxx-kXcenter,2) 
01662               + pow(reco_vtxy-kYcenter,2);
01663             reco_vtxr=sqrt(reco_vtxr);
01664           }
01665           else reco_vtxr = sqrt (pow(reco_vtxx,2)+pow(reco_vtxy,2));
01666           
01667           // compute signal in fully and partially instrumented regions
01668           SigInOut(track_index,evtsigfull,evtsigpart,trksigfull,trksigpart);
01669            if(detector==2){
01670             //if far det, all signal is in fully instrumented region
01671             evtsigfull=ntpEvent->ph.sigcor;
01672             evtsigpart=0.;
01673             if(track_index!=-1){
01674               trksigfull=ntpTrack->ph.sigcor;
01675             }
01676             else{
01677               trksigfull=0.;
01678             }
01679             trksigpart=0.;
01680           }
01681                  
01682           // angle reconstruction
01683           reco_dircosz      = RecoMuDCosZ(track_index);
01684           
01685           // event vertex
01686           float vtx_array[3]; 
01687           vtx_array[0]=ntpEvent->vtx.x; vtx_array[2]=ntpEvent->vtx.y;
01688           vtx_array[2]=ntpEvent->vtx.z;   
01689           if(detector==1) 
01690             reco_dircosneu = RecoMuDCosNeuND(track_index, vtx_array);
01691           else reco_dircosneu = RecoMuDCosNeuFD(track_index, vtx_array);
01692         
01693           evtlength         = ntpEvent->plane.end-ntpEvent->plane.beg;
01694           if(track_index!=-1){ //check track is present before using ntpTrack
01695 
01696             trklength         = ntpTrack->plane.end-ntpTrack->plane.beg;
01697             ntrklike =ntpTrack->plane.ntrklike;
01698             trkend         = ntpTrack->plane.end;
01699             trkendu         = ntpTrack->plane.endu;
01700             trkendv         = ntpTrack->plane.endv;
01701             trknstp       = ntpTrack->nstrip;
01702             trkmomrange       = ntpTrack->momentum.range;
01703             trkmomrange = CorrectMomentumFromRange(trkmomrange,!file_is_mc,det);
01704             trkqp             = ntpTrack->momentum.qp;
01705             if(trkqp!=0){
01706               trkqp = 1.0/CorrectSignedMomentumFromCurvature(1.0/trkqp,
01707                                                              !file_is_mc, det);
01708             }
01709 
01710             trkeqp            = ntpTrack->momentum.eqp;
01711             trkphfrac         = ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
01712             trkphplane        = ntpTrack->ph.mip/ntpTrack->plane.n;
01713             trkchi2           = ntpTrack->fit.chi2;
01714             trkndf           = ntpTrack->fit.ndof;
01715             reco_trk_vtxx         = ntpTrack->vtx.x;
01716             reco_trk_vtxy         = ntpTrack->vtx.y;
01717             reco_trk_vtxz         = ntpTrack->vtx.z;
01718 
01719             reco_trk_endx         = ntpTrack->end.x;
01720             reco_trk_endy         = ntpTrack->end.y;
01721             reco_trk_endz         = ntpTrack->end.z;
01722             // compute CC event class as defined by Pitt group
01723             pitt_evt_class = PittTrkContained(ntpTrack->end.x,
01724                                               ntpTrack->end.y,
01725                                               ntpTrack->end.z,
01726                                               ntpTrack->end.u,
01727                                               ntpTrack->end.v);
01728             duvvtx=abs(ntpTrack->plane.begu - ntpTrack->plane.begv);
01729             // based on the event class, go back and recompute
01730             // the muon and neutrino energy
01731             if(detector==1){//if it's the near detector
01732               if(pitt_evt_class>0){ // if it was classified
01733                 int do_meth=0; // emu reco method we should use
01734                 if( (pitt_evt_class == 1) || (pitt_evt_class == 3) ) 
01735                   do_meth=2; // stoppers --> range
01736                 else if( (pitt_evt_class == 2) || (pitt_evt_class == 4) ) 
01737                   do_meth=1; // punch-through --> curvature
01738                 if(do_meth==1 || do_meth==2){
01739                    reco_emu = RecoMKMuEnergy(do_meth,track_index,
01740                                              !file_is_mc,det);
01741                    reco_enu = reco_emu+reco_eshw;
01742                 }
01743                 emu_meth=do_meth; // was forgetting about this
01744               }
01745             }
01746             else if(detector==2){//if it's the near detector
01747               int do_meth=FarTrkContained(ntpTrack->end.x,
01748                                           ntpTrack->end.y,
01749                                           ntpTrack->end.z);
01750               if(do_meth==1 || do_meth==2){
01751                 reco_emu = RecoMKMuEnergy(do_meth,track_index);
01752                 reco_enu = reco_emu+reco_eshw;
01753               }
01754               emu_meth=do_meth; // was forgetting about this
01755             }                               
01756 
01757             //check to see if track vertex is in fid. vol
01758             if(detector==1){
01759               is_trk_fid = PittInFidVol(ntpTrack->vtx.x, ntpTrack->vtx.y, 
01760                                         ntpTrack->vtx.z,
01761                                         ntpTrack->vtx.u, ntpTrack->vtx.v);
01762             }
01763             else{
01764               is_trk_fid = InFidVol(ntpTrack->vtx.x,ntpTrack->vtx.y,
01765                                     ntpTrack->vtx.z);
01766             }
01767             
01768             // look for showers along muon track
01769             reco_mushw=MuonShowerEnergy(ntpEvent, ntpTrack);
01770             
01771             // CEDAR apparently has 5 plane long tracks, with ndf==0
01772             Double_t tmp_ndof=ntpTrack->fit.ndof;
01773             if(tmp_ndof<=0.0) tmp_ndof=1.0;
01774             // pitt tracking quality cuts
01775             if( (ntpTrack->fit.pass>0) 
01776                 && (ntpTrack->fit.chi2/tmp_ndof <20 )
01777                 && abs(ntpTrack->plane.begu - ntpTrack->plane.begv)<6 )
01778               pitt_trk_qual=1;
01779             
01781             const double M=(0.93827 + 0.93957)/2.0; // <nucleon mass>
01782             if(reco_emu>0) reco_Q2 = 
01783                              2*reco_enu*reco_emu*(1.0 - reco_dircosneu);
01784             reco_W2 = M*M - reco_Q2 + 2*M*reco_eshw;      
01785             if(reco_enu>0) reco_y = reco_eshw/reco_enu;
01786             if(reco_eshw>0 && reco_Q2>0) reco_x =  reco_Q2/(2*M*reco_eshw);
01788 
01790             const double muM=0.10566; // muon mass
01791             const double muP=sqrt(reco_emu*reco_emu -muM*muM);
01792             reco_qe_enu = (M*reco_emu - muM*muM/2.0)
01793               /(M - reco_emu + muP*reco_dircosneu);
01794             reco_qe_q2 = abs(-2.0*reco_qe_enu*(reco_emu-muP*reco_dircosneu)+muM*muM);
01796             
01798 
01799             med_qe_pid=qeid.AltCalcQEID(this,event,reco_enu,reco_W2);
01800             med_qe_pass=qeid.PassMEDQECut(reco_enu,med_qe_pid);
01802 
01804             if(nsid.ChooseWeightFile(det,current_beam)){
01805               if(!nsid.GetPID(this,event,det,niki_cc_pid)) niki_cc_pid=-999;
01806             }
01807             else niki_cc_pid=-999;
01808             // original code made detector dependent fiducial volume cut here
01809             // but, it's not necessary since PID has already been
01810             // calculated...
01811             // must reevaluate efficiency or make same cut downstream       
01813 
01815             if(det==Detector::kFar && !file_is_mc){
01816               dpid.SetPHCorrection(1.018);
01817             }
01818             if(dpid.ChoosePDFs(det,current_beam)){
01819               dave_cc_pid = dpid.CalcPID(this,event,0);
01820             }
01822 
01823             TObject *recobj=0;
01824             if(record!=0){
01825               recobj=record;
01826             }
01827             else{
01828               recobj=strecord;
01829             }
01830               
01831     
01832             LoadTrack(track_index);
01833             //figure out planes in trk > 1.5 pe
01834             float temp_ph;
01835             trk_nplanes=NPlanesInObj(recobj,ntpTrack,1.5,temp_ph);
01836             //figure out strips in trk > 1.5 pe
01837             trk_nstrips=NStripsInObj(recobj,ntpTrack,1.5,temp_ph);
01838 
01839           }
01840           else {
01841             trklength         = 0;
01842             trkmomrange       = 0;
01843             trkqp             = 0;
01844             trkeqp            = 0;
01845             trkphfrac         = 0;
01846             trkphplane        = 0;
01847             is_trk_fid        = -1;
01848           }
01849 
01850           //this used to load the slice associated with the track
01851           //but we really want the slice associated witht the event
01852           //tv 10-19-05
01853           LoadSlice(ntpEvent->slc);
01854           //        cout<<"*********SNARL="<<snarl<<" "<<event
01855           //            <<"*********************************"<<endl;
01856           TObject *recobj=0;
01857           if(record!=0){
01858             recobj=record;
01859           }
01860           else{
01861             recobj=strecord;
01862           }
01863 
01864           float temp_ph;
01865           //figure out planes in slice > 1.5 pe
01866           slc_nplanes=NPlanesInObj(recobj,ntpSlice,1.5,temp_ph);
01867           //figure out planes in evt > 1.5 pe
01868           evt_nplanes=NPlanesInObj(recobj,ntpEvent,1.5,temp_ph);
01869           //figure out strips in slice > 1.5 pe
01870           slc_nstrips=NStripsInObj(recobj,ntpSlice,1.5,temp_ph);
01871           //figure out strips in evt > 1.5 pe
01872           evt_nstrips=NStripsInObj(recobj,ntpEvent,1.5,temp_ph);
01873           
01874           
01875           if(shower_index!=-1 && reco_eshw>0){ // case of no vertex shower
01876             LoadShower(shower_index);
01877             shwlength=ntpShower->plane.n;
01878             shw_nplanes_raw=NPlanesInObj(recobj,ntpShower,0.0,shw_ph_raw);
01879             shw_nplanes_cut=NPlanesInObj(recobj,ntpShower,1.5,shw_ph_cut);
01880             
01881             shwend=ntpShower->plane.end;
01882             shwbeg=ntpShower->plane.beg;
01883             reco_shw_vtxx=ntpShower->vtx.x;
01884             reco_shw_vtxy=ntpShower->vtx.y;
01885             reco_shw_vtxz=ntpShower->vtx.z;
01886             
01887             
01889             int ncl = ntpShower->ncluster;
01890 
01891             //Float_t   zvtx
01892             //Float_t   tposvtx
01893             //Float_t   slope
01894             //Float_t   avgdev
01895             /*
01896             static TH1F hswu("hswu","",400,4,-4);
01897             static TH1F hswu_wtd("hswu_wtd","",400,4,-4);
01898             static TH1F hswv("hswv","",400,4,-4);
01899             static TH1F hswv_wtd("hswv_wtd","",400,4,-4);
01900             hswu.Reset();
01901             hswu_wtd.Reset();
01902             hswv.Reset();
01903             hswv_wtd.Reset();
01904             */
01905             Moments mom_swu;
01906             Moments mom_swu_wtd;
01907             Moments mom_swv;
01908             Moments mom_swv_wtd;
01909 
01910             etu=etv=elu=elv=0.0; 
01911             swu=swv=wswu=wswv=0.0;
01912 
01913             //      const float ss_cut=0.150;//GeV
01914             const float ss_cut=0.0;//GeV
01915             float totsigssu=0.0; // em and hadr subshowers
01916             float totsigssv=0.0; // em and hadr subshowers
01917             for(int ii=0; ii<ncl; ii++){
01918               if(LoadCluster(ntpShower->clu[ii])){
01919                 if(ntpCluster->ph.gev >ss_cut){
01920                   nss[ntpCluster->id]+=1;
01921                   ess[ntpCluster->id]+=ntpCluster->ph.gev;
01922                   ntotss++;
01923                   if(ntpCluster->id <2){
01924                     // slope = tpos/zpos?
01925                     float sinang = sin(atan(ntpCluster->slope));
01926                     if(ntpCluster->planeview==PlaneView::kU){
01927                       totsigssu+=ntpCluster->ph.gev;
01928                       //hswu.Fill(ntpCluster->tposvtx);
01929                       mom_swu.AddPoint(ntpCluster->tposvtx);
01930                       mom_swu_wtd.AddPoint(ntpCluster->tposvtx,
01931                                            ntpCluster->ph.gev);
01932                       
01933                       etu += sinang*ntpCluster->ph.gev;
01934                       elu += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev;
01935                     }
01936                     else if(ntpCluster->planeview==PlaneView::kV){
01937                       totsigssv+=ntpCluster->ph.gev;
01938                       //hswv.Fill(ntpCluster->tposvtx);
01939                       mom_swv.AddPoint(ntpCluster->tposvtx);
01940                       mom_swv_wtd.AddPoint(ntpCluster->tposvtx,
01941                                            ntpCluster->ph.gev);
01942 
01943                       etv += sinang*ntpCluster->ph.gev;
01944                       elv += sqrt(fabs(1-sinang*sinang))*ntpCluster->ph.gev;
01945                     }
01946                   }
01947                 }
01948               }
01949             }
01950             /*
01951             for(int ii=0; ii<ncl; ii++){
01952               if(LoadCluster(ntpShower->clu[ii])){
01953                 if(ntpCluster->ph.gev >ss_cut){
01954                   if(ntpCluster->id <2){
01955                     if(ntpCluster->planeview==PlaneView::kU){
01956                       totsigssu+=ntpCluster->ph.gev;
01957                       hswu_wtd.Fill(ntpCluster->tposvtx);
01958                     }
01959                     else if(ntpCluster->planeview==PlaneView::kV){
01960                       totsigssv+=ntpCluster->ph.gev;
01961                       hswv_wtd.Fill(ntpCluster->tposvtx);
01962                     }
01963                   }
01964                 }
01965               }
01966             }
01967             */
01968             /*
01969             swv=hswv.GetRMS();
01970             swu=hswu.GetRMS();
01971             wswv=hswv_wtd.GetRMS();
01972             wswu=hswu_wtd.GetRMS();
01973             */
01974             swv=mom_swv.GetRMS();
01975             swu=mom_swu.GetRMS();
01976             wswv=mom_swv_wtd.GetRMS();
01977             wswu=mom_swu_wtd.GetRMS();
01978             
01980           }
01981           if(ntrack==1 && reco_enu>0.0 && reco_enu<140
01982              && is_pitt_fid && pitt_trk_qual && (is_trk_fid!=0)) pass=1;
01983 
01984         } // if(reco_enu>0)
01985       } // if(PassBasicCuts())
01986       
01987       // for both near in time and near in space
01988       // index of nearest, dist of nearest, ph of nearest
01989       // near_vtx_space_idx, near_vtx_space_dist, near_vtx_space_ph 
01990       // nvsi, nvsd, nvsp
01991       // near_vtx_time_idx, near_vtx_time_dist, near_vtx_time_ph 
01992       // nvti, nvtd, nvtp
01993       // for the nearest event, the index of the event nearest it
01994       // ph of event nearest it, distance of event nearest it
01995       //  maybe use for event rehashing
01996       // back_vtx_space_idx, back_vtx_space_dist, back_vtx_space_ph 
01997       // bvsi, bvsd, bvsp
01998       // back_vtx_time_idx, back_vtx_time_dist, back_vtx_time_ph 
01999       // bvti, bvtd, bvtp
02000 
02001       nby.ClosestSpatial(i,nvsi,nvsd,nvst,nvsp,nvsevt);
02002       nby.ClosestSpatial(nvsi,bvsi,bvsd,bvst,bvsp,bvsevt);
02003       nby.ClosestTemporal(i,nvti,nvtd,nvtt,nvtp,nvtevt);
02004       nby.ClosestTemporal(nvti,bvti,bvtd,bvtt,bvtp,bvtevt);
02005 
02006       GetEvtTimeChar(evttimemin,evttimemax,evttimeln);
02007       GetEvtTimeCharPHC(evttimeminphc,evttimemaxphc,evttimelnphc);
02008       EarlyActivity(i,earlywph,earlywphpw,earlywphsphere,
02009                     ph1000,ph500,ph100,lateph1000,lateph500,
02010                     lateph100,lateph50);
02011       evtph = ntpEvent->ph.sigcor;
02012 
02013       DupRecoStrips(i,nduprs,dupphrs);
02014 
02016       if(file_is_mc){
02017         if(gnumi->Status()) {
02018           if(isST) gnumi->GetParent(strecord,mcevent,*nuparent);
02019           else gnumi->GetParent(mcrecord,mcevent,*nuparent);
02020         }               
02021       }
02023 
02027       if(fDoInukeReweight && file_is_mc && inukerw!=0) {
02028         
02029         if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents){
02030           std::vector<double> weights(n_inuke_wts,1.0);
02031           //    std::cout<<"snarl, event: "<<snarl<<", "<<event<<std::endl;
02032           inuke_reweight::calc_weights(strecord,ntpTruth,psets,weights,nucrad,wrad,false);
02033           for (unsigned int iwt=0; iwt<weights.size(); iwt++) inuke_wts[iwt]=weights[iwt];
02034         }
02035         // nuclear radius and weighting factor
02036 
02037       }
02041 
02043       if((fOnlyWriteFidEvents && is_pitt_fid==1) || !fOnlyWriteFidEvents) {
02044         tree->Fill(); if(inukerw!=0 && inukerw!=tree) inukerw->Fill();
02045       }
02046 
02047 
02048     } // loop over events
02049   } // loop over entries
02050   
02051   // delete nuparent; // should I do this?
02052   save.cd();
02053   pottree->Fill();
02054   pottree->Write();
02055   tree->Write();if(inukerw!=0 && inukerw!=tree) inukerw->Write();
02056 
02057   save.Write();
02058   save.Close();
02059 
02060 }

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

Definition at line 150 of file MadMKAnalysis.h.

References fDoInukeReweight.

00150 {fDoInukeReweight=t;}

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

Definition at line 3158 of file MadMKAnalysis.cxx.

References MadBase::eventSummary, it, MadBase::LoadEvent(), MadBase::LoadStrip(), NtpSREventSummary::nevent, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpStrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::sigcor, and NtpSREvent::stp.

Referenced by CreatePAN().

03159 {
03160   //zero the variables
03161   nduprs=0;
03162   dupphrs=0.;
03163 
03164   if(!LoadEvent(evidx)) return; //no event found
03165   
03166   map<int,float> sindex;
03167   for(int i=0;i<ntpEvent->nstrip;i++){
03168     int si = ntpEvent->stp[i];
03169     if(!LoadStrip(si)) continue;
03170     float ph = ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor;
03171     sindex[ntpEvent->stp[i]]=ph;
03172   }
03173 
03174   for(int i=0;i<eventSummary->nevent;i++){
03175     if(i==evidx){
03176       continue;
03177     }
03178     LoadEvent(i);  //load a new event
03179     for(int j=0;j<ntpEvent->nstrip;j++){
03180       map<int,float>::iterator it(sindex.find(ntpEvent->stp[j]));
03181       if(it!=sindex.end()){
03182         nduprs++;
03183         dupphrs+=it->second;
03184       }
03185     }
03186   }
03187 
03188   //reload the current event
03189   LoadEvent(evidx);
03190 }

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 MadBase::eventSummary, fEarlyActivityWindowSize, fEarlyPlaneSphere, fNPlaneEarlyAct, fPMTTimeConstant, PlaneView::kU, MadBase::LoadEvent(), MadBase::LoadStrip(), NtpSREventSummary::nstrip, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpStrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRVertex::plane, NtpSRStrip::planeview, NtpSRStrip::pmtindex0, NtpSRStrip::pmtindex1, 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().

02812 {
02813   //zero the variables
02814   earlywph=0.;
02815   earlywphpw=0.; 
02816   earlywphsphere=0.;
02817   ph1000=0.;
02818   ph500=0.;
02819   ph100=0.;
02820   lateph1000=0.;
02821   lateph500=0.;
02822   lateph100=0.;
02823   lateph50=0.;
02824 
02825   if(!LoadEvent(evidx)) return; //no event found
02826   
02827   //find the time of the earliest strip in the event
02828   double earliestEventTime = 1.e20;
02829   double latestEventTime = 0.;
02830   map<int,int> eventPlanes;
02831   double trigtime = eventSummary->trigtime;
02832 
02833   //loop over strips in the event
02834   //find the earliest strip time 
02835   //make a map of the pmt's hit in this event
02836   for(int s=0;s<ntpEvent->nstrip;s++){
02837     int stripindex = ntpEvent->stp[s];
02838     if(!LoadStrip(stripindex)) continue;
02839     //calculate ph weighted strip time (over the two ends,
02840     //should be equivalent to ph1 in ND
02841     float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+
02842                    ntpStrip->ph0.sigcor*ntpStrip->time0)/
02843                   (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor));
02844       
02845     double time = 1.e9*(phwt-trigtime);
02846     if(time<earliestEventTime){
02847       earliestEventTime = time;
02848     }
02849     if(time>latestEventTime){
02850       latestEventTime=time;
02851     }
02852     if(ntpStrip->ph1.sigcor>0){
02853       if(eventPlanes.find(ntpStrip->pmtindex1)==eventPlanes.end()){
02854         eventPlanes[ntpStrip->pmtindex1] = 1;
02855       }
02856     }
02857     if(ntpStrip->ph0.sigcor>0){
02858       if(eventPlanes.find(ntpStrip->pmtindex0)==eventPlanes.end()){
02859         eventPlanes[ntpStrip->pmtindex1] = 1;
02860       }
02861     }
02862   }
02863 
02864   //now find the weigthed sum of ADC for the early strips - make sure the early
02865   //strips come from the same pmts as the event strips.
02866   for(unsigned int s=0;s<eventSummary->nstrip;s++){
02867     if(!LoadStrip(s)) continue;
02868     //work in ns
02869     float phwt = ((ntpStrip->ph1.sigcor*ntpStrip->time1+
02870                    ntpStrip->ph0.sigcor*ntpStrip->time0)/
02871                   (ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor));
02872     double newtime = 1.e9*(phwt-trigtime);
02873     double dt=(earliestEventTime-newtime);
02874     double dtlate=newtime-latestEventTime;
02875     float d=1000.;
02876     if((int)ntpStrip->planeview==PlaneView::kU){
02877       d=sqrt((pow(ntpEvent->vtx.u-ntpStrip->tpos,2)+
02878               pow(ntpEvent->vtx.z-ntpStrip->z,2)));
02879     }
02880     else{
02881       d=sqrt((pow(ntpEvent->vtx.v-ntpStrip->tpos,2)+
02882               pow(ntpEvent->vtx.z-ntpStrip->z,2)));
02883     }
02884     if(d<fEarlyPlaneSphere&&dtlate>0.){
02885       if(dtlate<1000.){
02886         lateph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02887       }
02888       if(dtlate<500.){
02889         lateph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02890       }
02891       if(dtlate<100.){
02892         lateph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02893       }
02894       if(dtlate<50.){
02895         lateph50+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02896       }
02897     }   
02898     if(d<fEarlyPlaneSphere&&dt>0.){
02899       if(dt<1000.){
02900         ph1000+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02901       }
02902       if(dt<500.){
02903         ph500+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02904       }
02905       if(dt<100.){
02906         ph100+=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
02907       }
02908     }
02909     if(dt>0.&&dt<1000.*fEarlyActivityWindowSize){
02910       double eadc = ((ntpStrip->ph1.sigcor+ntpStrip->ph1.sigcor)*
02911                      TMath::Exp(-1.*dt/fPMTTimeConstant));
02912       //if this strip has activity in the same pmt as our event,
02913       //sum up the early pulse height
02914       if(eventPlanes.find(ntpStrip->pmtindex1)!=eventPlanes.end()){
02915         earlywph += eadc;
02916       }      
02917       //if this strip is within some radius of event vertex, 
02918       //add up the early pulse height
02919       if(d<fEarlyPlaneSphere){
02920         earlywphsphere+=eadc;
02921       }
02922       //if this strip is within the plane window of event vertex,
02923       //add up the early pulseheight
02924       if(ntpStrip->plane>=ntpEvent->vtx.plane&&
02925          ntpStrip->plane<=ntpEvent->vtx.plane+fNPlaneEarlyAct){
02926         earlywphpw+=eadc;
02927       }
02928     }
02929   }
02930 }

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

Definition at line 3118 of file MadMKAnalysis.cxx.

Referenced by CreatePAN(), MadEvDisplay::DrawTextBox(), and MadScanDisplay::DrawTextBox().

03119 {
03120   int is_cont=0;
03121   //arguments should correspond to track end
03122   //returns:
03123   //1 if track is partially contained
03124   //2 if  track is fully contained
03125   const Float_t r = sqrt(x*x+y*y);
03126 
03127   if(r<=3.5&&r>.5&&((z>=0.5&&z<=14.5)||(z>=16.5&&z<=29.4))){//contained
03128     is_cont=2;
03129   }
03130   else{//not contained
03131     is_cont=1;
03132   }
03133   return is_cont;
03134 }

void MadMKAnalysis::FillMCHists ( TFile *  fout  ) 

Definition at line 2933 of file MadMKAnalysis.cxx.

References NtpMCTruth::iaction, NtpMCTruth::inu, NtpMCTruth::iresonance, MadBase::LoadTruth(), MadBase::mcHeader, NDRadialFidVol(), NtpMCSummary::nmc, MadBase::ntpTruth, NtpMCTruth::p4neu, PittInFidVol(), NtpMCTruth::vtxx, NtpMCTruth::vtxy, and NtpMCTruth::vtxz.

Referenced by CreatePAN().

02933                                           {
02934 
02935   static bool first=true;
02936 
02937   static TH1* h_numu_pf_tot = new TH1F("h_numu_pf_tot","CC; true neutrino energy (GeV)",
02938                                  240,0.0,120.0);
02939 
02940   static TH1* h_numu_pf_qe = new TH1F("h_numu_pf_qe","QE; true neutrino energy (GeV)",
02941                                  240,0.0,120.0);
02942 
02943   static TH1* h_numu_pf_res = new TH1F("h_numu_pf_res","RES; true neutrino energy (GeV)",
02944                                  240,0.0,120.0);
02945 
02946   static TH1* h_numu_pf_dis = new TH1F("h_numu_pf_dis","DIS; true neutrino energy (GeV)",
02947                                  240,0.0,120.0);
02948   static TH1* h_numu_pf_nc = new TH1F("h_numu_pf_nc","NC; true neutrino energy (GeV)",
02949                                  240,0.0,120.0);
02950   if(first){
02951     h_numu_pf_tot->SetDirectory(fout);
02952     h_numu_pf_qe->SetDirectory(fout);
02953     h_numu_pf_res->SetDirectory(fout);
02954     h_numu_pf_dis->SetDirectory(fout);
02955     h_numu_pf_nc->SetDirectory(fout);
02956   }
02957 
02958   static TH1* h_numubar_pf_tot = new TH1F("h_numubar_pf_tot","CC; true neutrino energy (GeV)",
02959                                  240,0.0,120.0);
02960 
02961   static TH1* h_numubar_pf_qe = new TH1F("h_numubar_pf_qe","QE; true neutrino energy (GeV)",
02962                                  240,0.0,120.0);
02963 
02964   static TH1* h_numubar_pf_res = new TH1F("h_numubar_pf_res","RES; true neutrino energy (GeV)",
02965                                  240,0.0,120.0);
02966 
02967   static TH1* h_numubar_pf_dis = new TH1F("h_numubar_pf_dis","DIS; true neutrino energy (GeV)",
02968                                  240,0.0,120.0);
02969   static TH1* h_numubar_pf_nc = new TH1F("h_numubar_pf_nc","NC; true neutrino energy (GeV)",
02970                                  240,0.0,120.0);
02971   if(first){
02972     h_numubar_pf_tot->SetDirectory(fout);
02973     h_numubar_pf_qe->SetDirectory(fout);
02974     h_numubar_pf_res->SetDirectory(fout);
02975     h_numubar_pf_dis->SetDirectory(fout);
02976     h_numubar_pf_nc->SetDirectory(fout);
02977   }
02978 
02979 
02980 
02981   static TH1* h_numu_1m_tot = new TH1F("h_numu_1m_tot","CC; true neutrino energy (GeV)",
02982                                  240,0.0,120.0);
02983 
02984   static TH1* h_numu_1m_qe = new TH1F("h_numu_1m_qe","QE; true neutrino energy (GeV)",
02985                                  240,0.0,120.0);
02986 
02987   static TH1* h_numu_1m_res = new TH1F("h_numu_1m_res","RES; true neutrino energy (GeV)",
02988                                  240,0.0,120.0);
02989 
02990   static TH1* h_numu_1m_dis = new TH1F("h_numu_1m_dis","DIS; true neutrino energy (GeV)",
02991                                  240,0.0,120.0);
02992   static TH1* h_numu_1m_nc = new TH1F("h_numu_1m_nc","NC; true neutrino energy (GeV)",
02993                                  240,0.0,120.0);
02994   if(first){
02995     h_numu_1m_tot->SetDirectory(fout);
02996     h_numu_1m_qe->SetDirectory(fout);
02997     h_numu_1m_res->SetDirectory(fout);
02998     h_numu_1m_dis->SetDirectory(fout);
02999     h_numu_1m_nc->SetDirectory(fout);
03000   }
03001 
03002   static TH1* h_numubar_1m_tot = new TH1F("h_numubar_1m_tot","CC; true neutrino energy (GeV)",
03003                                  240,0.0,120.0);
03004 
03005   static TH1* h_numubar_1m_qe = new TH1F("h_numubar_1m_qe","QE; true neutrino energy (GeV)",
03006                                  240,0.0,120.0);
03007 
03008   static TH1* h_numubar_1m_res = new TH1F("h_numubar_1m_res","RES; true neutrino energy (GeV)",
03009                                  240,0.0,120.0);
03010 
03011   static TH1* h_numubar_1m_dis = new TH1F("h_numubar_1m_dis","DIS; true neutrino energy (GeV)",
03012                                  240,0.0,120.0);
03013   static TH1* h_numubar_1m_nc = new TH1F("h_numubar_1m_nc","NC; true neutrino energy (GeV)",
03014                                  240,0.0,120.0);
03015   if(first){
03016     h_numubar_1m_tot->SetDirectory(fout);
03017     h_numubar_1m_qe->SetDirectory(fout);
03018     h_numubar_1m_res->SetDirectory(fout);
03019     h_numubar_1m_dis->SetDirectory(fout);
03020     h_numubar_1m_nc->SetDirectory(fout);
03021   }
03022 
03023 
03024   static TH1* h_other_pf = new TH1F("h_other_pf",
03025                                     "other; true neutrino energy (GeV)",
03026                                     240,0.0,120.0);
03027 
03028   static TH1* h_other_1m = new TH1F("h_other_1m",
03029                                     "other; true neutrino energy (GeV)",
03030                                     240,0.0,120.0);
03031   if(first){
03032     h_other_pf->SetDirectory(fout);
03033     h_other_1m->SetDirectory(fout);
03034   }
03035   if(!mcHeader) {
03036     first=true;
03037     return;
03038   }
03039   for(Int_t i=0; i<mcHeader->nmc; i++){
03040     if(!LoadTruth(i)) continue; // no ntpmctruth
03041     const Float_t k = 1.0/sqrt(2.0);
03042     Float_t x = ntpTruth->vtxx;
03043     Float_t y = ntpTruth->vtxy;
03044     Float_t z = ntpTruth->vtxz;
03045     Float_t u = k*(y+x);
03046     Float_t v = k*(y-x);
03047     Bool_t inpf = PittInFidVol(x,y,z,u,v);    
03048     Int_t ndr = NDRadialFidVol(x,y,z);
03049     Bool_t inr =kFALSE;
03050     if((ndr&0x8)>0) inr=kTRUE;
03051     Int_t ccnc = ntpTruth->iaction;
03052     Int_t process = ntpTruth->iresonance;
03053     Int_t inu = ntpTruth->inu;
03054     Float_t E = ntpTruth->p4neu[3];
03055 
03056     
03057     if(inu==14){ // neutrinos
03058       if(ccnc==0){
03059         if(inpf) h_numu_pf_nc->Fill(E);
03060         if(inr) h_numu_1m_nc->Fill(E);
03061       }
03062       else if(ccnc==1){
03063         if(inpf) h_numu_pf_tot->Fill(E);
03064         if(inr) h_numu_1m_tot->Fill(E);
03065         switch(process){
03066         case 1001:
03067           if(inpf) h_numu_pf_qe->Fill(E);
03068           if(inr) h_numu_1m_qe->Fill(E);
03069           break;
03070         case 1002:
03071           if(inpf) h_numu_pf_res->Fill(E);
03072           if(inr) h_numu_1m_res->Fill(E);
03073           break;
03074         case 1003:
03075           if(inpf) h_numu_pf_dis->Fill(E);
03076           if(inr) h_numu_1m_dis->Fill(E);
03077           break;
03078         default:
03079           break;
03080         }
03081       }
03082     }
03083     else if(inu==-14){
03084       if(ccnc==0){
03085         if(inpf) h_numubar_pf_nc->Fill(E);
03086         if(inr) h_numubar_1m_nc->Fill(E);
03087       }
03088       else if(ccnc==1){
03089         if(inpf) h_numubar_pf_tot->Fill(E);
03090         if(inr) h_numubar_1m_tot->Fill(E);
03091         switch(process){
03092         case 1001:
03093           if(inpf) h_numubar_pf_qe->Fill(E);
03094           if(inr) h_numubar_1m_qe->Fill(E);
03095           break;
03096         case 1002:
03097           if(inpf) h_numubar_pf_res->Fill(E);
03098           if(inr) h_numubar_1m_res->Fill(E);
03099           break;
03100         case 1003:
03101           if(inpf) h_numubar_pf_dis->Fill(E);
03102           if(inr) h_numubar_1m_dis->Fill(E);
03103           break;
03104         default:
03105           break;
03106         }
03107       }
03108     }
03109     else{
03110       if(inpf) h_other_pf->Fill(E);
03111       if(inr) h_other_1m->Fill(E);
03112     }
03113   }
03114   
03115   first=false;
03116 }

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

Definition at line 142 of file MadMKAnalysis.h.

References fMCBeam.

00142 {fMCBeam=beam;}

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

Definition at line 187 of file MadMKAnalysis.h.

References fBECConfig.

00187 { return fBECConfig; }

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

Definition at line 2636 of file MadMKAnalysis.cxx.

References MadBase::eventSummary, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadStrip(), NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpStrip, NtpSREvent::stp, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSREventSummary::trigtime.

Referenced by CreatePAN().

02637 {
02638   //cut and pasted, with minor changes from MadDpAnalysis.cxx
02639   double trgtime=eventSummary->trigtime;
02640   double timemax=0.;
02641   double timemin=1.e10;
02642 
02643   //  Find max min time of events by loading strips for ND 
02644   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02645     for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
02646       Int_t stp_index=((ntpEvent->stp)[evss]);
02647       if(!LoadStrip(stp_index)) continue;
02648       Double_t striptime=ntpStrip->time1-trgtime;
02649       if(striptime<=timemin) timemin=striptime;
02650       if(striptime>=timemax) timemax=striptime;
02651     }
02652   }      
02653   
02654   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02655     for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
02656       Int_t stp_index=((ntpEvent->stp)[evss]);
02657       if(!LoadStrip(stp_index)) continue;
02658       Double_t striptime1=ntpStrip->time1-trgtime;
02659       Double_t striptime0=ntpStrip->time0-trgtime;
02660       Double_t striptime = 0;
02661       if(striptime1>0 && striptime0<0)  striptime=striptime1;
02662       if(striptime0>0 && striptime1<0)  striptime=striptime0;
02663       if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
02664       if(striptime<=timemin) timemin=striptime;
02665       if(striptime>=timemax) timemax=striptime;
02666     }
02667   }
02668 
02669     evttimemax=timemax;
02670     evttimemin=timemin;
02671     evttimeln=timemax-timemin;
02672 
02673 }

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

Definition at line 2675 of file MadMKAnalysis.cxx.

References MadBase::eventSummary, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadStrip(), NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpStrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRStrip::time0, NtpSRStrip::time1, and NtpSREventSummary::trigtime.

Referenced by CreatePAN().

02676 {
02677   //same as GetEvtTimeChar but only looking at hits with 
02678   //pulse height > 200 sig cor (~2pe) in near det
02679   //pulse height > 160 sig cor in far det
02680   double trgtime=eventSummary->trigtime;
02681   double timemax=0.;
02682   double timemin=1.e10;
02683 
02684   //  Find max min time of events by loading strips for ND 
02685   if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02686     for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
02687       Int_t stp_index=((ntpEvent->stp)[evss]);
02688       if(!LoadStrip(stp_index)) continue;
02689       Double_t striptime=ntpStrip->time1-trgtime;
02690       if(ntpStrip->ph1.sigcor>200){
02691         if(striptime<=timemin) timemin=striptime;
02692         if(striptime>=timemax) timemax=striptime;
02693       }
02694     }
02695   }      
02696   
02697   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02698     for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
02699       Int_t stp_index=((ntpEvent->stp)[evss]);
02700       if(!LoadStrip(stp_index)) continue;
02701       Double_t striptime1=ntpStrip->time1-trgtime;
02702       Double_t striptime0=ntpStrip->time0-trgtime;
02703       Double_t striptime = 0;
02704       if(striptime1>0 && striptime0<0)  striptime=striptime1;
02705       if(striptime0>0 && striptime1<0)  striptime=striptime0;
02706       if(striptime0>0 && striptime1>0)  striptime=(striptime0+striptime1)/2.;
02707       if(ntpStrip->ph0.sigcor+ntpStrip->ph1.sigcor>160){
02708         if(striptime<=timemin) timemin=striptime;
02709         if(striptime>=timemax) timemax=striptime;
02710       }
02711     }
02712   }
02713   
02714   evttimemax=timemax;
02715   evttimemin=timemin;
02716   evttimeln=timemax-timemin;
02717 
02718 }

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

Definition at line 143 of file MadMKAnalysis.h.

References fMCBeam.

00143 {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(), Detector::kFar, Detector::kNear, kXcenter, kYcenter, and MadBase::ntpHeader.

02072                                                         {
02073   // assumes event has been loaded
02074   // check that vertex is in fiducial volume
02075   bool is_fid=false;
02076 
02077   // FD: from doc-1351-v2 : there is no coil cut!
02078   if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){
02079         //fiducial criteria on vtx for far detector
02080     is_fid = true;
02081     if(vz<0.5 || vz>28.0 ||   //ends
02082        (vz<16.2&&vz>14.3) ||  //between SMs
02083        ((vx*vx)+(vy*vy))>14.0) is_fid = false;
02084   }
02085   else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
02086     //fiducial criteria on vtx for near detector
02087 
02088     //    float beamzerox = 1.4885;
02089     //    float beamzeroy = 0.1397-reco_vtxz*TMath::Tan(.058);
02090 
02091     //    float r=sqrt(pow((reco_vtxx-beamzerox),2)+pow((reco_vtxy-beamzeroy),2));
02092 
02093 
02094     is_fid = false;
02095     // MAK: Sept 19, 2005: Modified to account for beam slope
02096     //    double radius = pow(vx-kXcenter,2) 
02097     //      + pow(vy-(kYcenter -vz*TMath::Tan(0.058)),2);
02098     double radius = pow(vx-kXcenter,2) + pow(vy-kYcenter,2);
02099 
02100     radius=sqrt(radius);
02101     // restrict vertices to target region and 1m radius
02102     /*
02103     if(vz>=kVetoEndZ && 
02104        vz<=kTargEndZ &&
02105        radius<=0.5) is_fid = true;    
02106     */
02107     // restrict vertices to target region and 1m radius
02108     if(vz>=1.0 && vz<=5.0 && radius<=1.0) is_fid = true;
02109   }
02110   return is_fid;
02111 }

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

Definition at line 2067 of file MadMKAnalysis.cxx.

References MadBase::ntpEvent, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by CreatePAN(), and InFidVol().

02068 {
02069   return InFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,ntpEvent->vtx.z);
02070 }

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

Definition at line 2062 of file MadMKAnalysis.cxx.

References InFidVol(), and MadBase::LoadEvent().

02062                                        {
02063   if(!LoadEvent(evidx)) return false; //no event found
02064   return InFidVol();
02065 }

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

Definition at line 2350 of file MadMKAnalysis.cxx.

Referenced by SigInOut().

02350                                                                   {
02351   // figure out if this (plane,strip) corresponds to something in
02352   // the partially instrumented region
02353   //
02354   // this region is defined as:
02355   // v planes: (strip<=4 || strip>=67)
02356   // partial u: (strip==0 || strip=63)
02357   // full u: (strip<=26 || strip>=88) 
02358   //
02359   // if so, return 1
02360   // if not, return -1
02361   // if error, return 0
02362 
02363 
02364   // make a lookup ptype to hold the type of each plane
02365   // 1 = v partial   2 = u partial
02366   // 3 = v full   4 = u full
02367   // 0 = uninstrumented
02368   static bool first=true;
02369   static UShort_t ptype[282];
02370   if(first){
02371     ptype[0]=0;
02372     for(int i=1; i<=281; i++){
02373       if(i%2==0) ptype[i]=1; // a v plane
02374       else ptype[i]=2; // a u plane
02375       if((i-1)%5 == 0) ptype[i]+=2; // fully instrumented
02376       else if(i>120) ptype[i]=0; // not instrumented
02377     }
02378     first=false;
02379   }
02380   if(plane>281){
02381     //    std::cerr<<"InPartialRegion passed plane = "<<plane<<std::endl;
02382     return 0;
02383   }
02384   UShort_t pt = ptype[plane];
02385   
02386   Int_t result;
02387   switch(pt){
02388   case 1:
02389   case 3:
02390     if(strip<=4 || strip>=67) result=1;
02391     else result = -1;
02392     break;
02393   case 2:
02394     if(strip==0 || strip == 63) result=1;
02395     else result = -1;
02396     break;
02397   case 4:
02398     if(strip<=26 || strip>=88) result=1;
02399     else result = -1;
02400     break;
02401   case 0:
02402   default:
02403     result=0;
02404     break;
02405   }
02406   return result;
02407 
02408 }

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

Definition at line 2722 of file MadMKAnalysis.cxx.

References MuELoss::e, NtpSRShowerPulseHeight::EMgev, MadBase::LoadShower(), NtpSREvent::nshower, NtpSRTrack::nstrip, MadBase::ntpEvent, MadBase::ntpShower, NtpSREvent::shw, NtpSRShower::shwph, NtpSRTrack::stpu, NtpSRTrack::stpv, NtpSRTrack::stpz, NtpSRVertex::t, NtpSRVertex::u, NtpSRVertex::v, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, and NtpSRVertex::z.

Referenced by CreatePAN().

02722                                                                                 {  
02723   // Purpose: try to sum up brems along the track
02724 
02725   static TH2* hr = new TH2F("h_mushw_rph",
02726                             "; track - shw r dist (m); shower energy (GeV)", 
02727                             10,0,0.2,40,0,2);
02728 
02729   static TH2* hz = new TH2F("h_mushw_zph",
02730                             "; track - shw z dist (m); shower energy (GeV)", 
02731                             20,-0.2,0.2,40,0,2);
02732 
02733   static TH2* ht = new TH2F("h_mushw_tph",
02734                             "; track - shw t dist (m); shower energy (GeV)", 
02735                             600,-1e-6,10e-6,40,0,2);
02736   
02737 
02738   if(evt==0 || trk==0) return 0;
02739   float result=0;
02740   NtpSRShower* save_shwr = ntpShower;  
02741   // loop through all showers in event
02742   const float max_vtx_dist=14*0.06; // 2 interaction lengths
02743   for(int i=0; i<evt->nshower; i++){
02744     LoadShower(evt->shw[i]);
02745     if(!ntpShower) continue;
02746     // calculate distance to vertex
02747     // don't want to use showers near the vertex
02748     float dist_vtx = pow(ntpShower->vtx.z - ntpEvent->vtx.z,2)
02749       + pow(ntpShower->vtx.u - ntpEvent->vtx.u,2)
02750       + pow(ntpShower->vtx.v - ntpEvent->vtx.v,2);
02751     dist_vtx=sqrt(dist_vtx);
02752     if(dist_vtx<max_vtx_dist) continue;
02753     // loop through track's strip array
02754     // calculate distance from strip to shower in time and space
02755     // save spatial and time distance of closest approach to track
02756     // if shower is close enough add the shower energy to running sum.
02757     float min_dist=1e6; float min_r=1e6; float min_z=1e6; float min_t=1e6;
02758     float min_ph=0; float add_ph=0;
02759     for(int j=0; j<trk->nstrip; j++){
02760       float dist_trk= pow(trk->stpz[j]-ntpShower->vtx.z,2)
02761         + pow(trk->stpu[j]-ntpShower->vtx.u,2)
02762         + pow(trk->stpv[j]-ntpShower->vtx.v,2);
02763       dist_trk=sqrt(dist_trk);
02764       float tdist = (trk->vtx.t + (trk->stpz[j]-trk->vtx.z)/3e8) 
02765                          - ntpShower->vtx.t;
02766       
02767         // = trk->stpt0[j]*trk->stpph0[j]+trk->stp1[j]*trk->stpph1[j];
02768         //      if( (trk->stpph0[j]+trk->stpph1[j]) > 0){
02769         //      tdist/=(trk->stpph0[j]+trk->stpph1[j]);
02770         //      }
02771       
02772       if(fabs(dist_trk)<fabs(min_dist)) {
02773         min_ph=ntpShower->shwph.EMgev;
02774         min_z = fabs(ntpShower->vtx.z-trk->stpz[j]);
02775         min_r = sqrt( pow(ntpShower->vtx.u-trk->stpu[j],2)
02776                       + pow(ntpShower->vtx.v-trk->stpv[j],2));
02777         min_t = tdist;
02778       }
02779 
02780       const float max_trk_dist=1.76*3.0; // 3 radiation lengths
02781       const float max_trk_tdist=40e-9;
02782       if(fabs(dist_trk)<max_trk_dist && fabs(tdist)<max_trk_tdist){
02783         // shower correllated with track
02784         add_ph=min_ph;
02785       }
02786     }
02787     result+=min_ph; // 
02788     
02789     // fill histograms
02790     hr->Fill(min_r,min_ph);
02791     hz->Fill(min_z,min_ph);
02792     ht->Fill(min_t,min_ph);
02793 
02794   }// end loop over showers
02795 
02796   // reset shower pointer
02797   ntpShower=save_shwr;
02798   // all done
02799   return result;
02800 }

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().

02266                                                                  {
02267 
02268   Float_t rings[5]={0.1,0.25,0.5,1.0,1.5}; // in meters
02269   Int_t result=0;
02270   for (unsigned int i=0; i<5; i++){
02271     Float_t radius = pow(x-kXcenter,2) 
02272       + pow(y-(kYcenter - z*TMath::Tan(0.058)),2);
02273     radius=sqrt(radius);
02274     // restrict vertices to target region and 1m radius
02275     if(z>=kVetoEndZ && 
02276        z<=kTargEndZ &&
02277        radius<=rings[i]) result|=(1<<i);
02278   }
02279   return result;
02280 }

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, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRShower::stp, NtpStRecord::stp, NtpSRSlice::stp, NtpSRRecord::stp, and NtpSRTrack::stp.

Referenced by CreatePAN().

02469 {
02470   sumph=0.0;
02471 
02472   std::set<UShort_t> planes;
02473 
02474   NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec);
02475   NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec);
02476 
02477   Int_t *stripindex=0;
02478   Int_t nstrip=0;
02479   NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj);
02480   NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj);
02481   NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj);
02482   NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj);
02483 
02484   if(slc!=0){
02485     stripindex=slc->stp;
02486     nstrip=slc->nstrip;
02487     //    cout<<"Counting planes in slice"<<endl;
02488   }
02489   else if(trk!=0){
02490     stripindex=trk->stp;
02491     nstrip=trk->nstrip;
02492     //    cout<<"Counting planes in trk"<<endl;
02493   }
02494   else if(evt!=0){
02495     stripindex=evt->stp;
02496     nstrip=evt->nstrip;
02497     //    cout<<"Counting planes in event"<<endl;
02498   }
02499   else if(shw!=0){
02500     stripindex=shw->stp;
02501     nstrip=shw->nstrip;
02502     //    cout<<"Counting planes in event"<<endl;
02503   }
02504   
02505   TClonesArray *striplist=0; 
02506   int TOTSTRIPS=0;
02507   if(srrec!=0){
02508     striplist = srrec->stp;
02509     TOTSTRIPS=srrec->evthdr.nstrip;
02510     //    cout<<"Found a sr"<<endl;
02511   }
02512   else{
02513     striplist = strec->stp;
02514     TOTSTRIPS=strec->evthdr.nstrip;
02515     //    cout<<"found a st"<<endl;
02516   }
02517   
02518   if(striplist!=0){
02519     //    cout<<"Looping over strip list, nstrip="<<nstrip<<endl;
02520     for(int i=0;i<nstrip;i++){
02521       int sindex = stripindex[i];
02522       if(sindex<0) continue;
02523       if(sindex<TOTSTRIPS){       
02524         NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>
02525           ((*striplist)[sindex]);
02526         if(!strip) continue;
02527         //      cout<<"got strip "<<sindex<<" plane "<<strip->plane
02528         //          <<" strip "<<strip->strip
02529         //          <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl;
02530         if(strip->ph0.pe+strip->ph1.pe>phcut){
02531           planes.insert(strip->plane);
02532           sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
02533         }     
02534       }
02535     }
02536   }
02537   else{
02538     //    cout<<"Getting number of planes in snarl"<<endl;
02539     for(int i=0;i<TOTSTRIPS;i++){
02540       NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]);
02541       if(!strip) continue;
02542       //      cout<<"got strip "<<i<<" plane "<<strip->plane
02543       //          <<" strip "<<strip->strip
02544       //          <<" ph "<<strip->ph0.pe+strip->ph1.pe<<endl;
02545       if(strip->ph0.pe+strip->ph1.pe>phcut){
02546         planes.insert(strip->plane);
02547         sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
02548       }     
02549     }
02550   }
02551     
02552   //  cout<<"I count "<<planes.size()<<" planes > "<<phcut<<" pe"<<endl;
02553   return planes.size();
02554 }

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, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRShower::stp, NtpStRecord::stp, NtpSRSlice::stp, NtpSRRecord::stp, and NtpSRTrack::stp.

Referenced by CreatePAN().

02558 {
02559   sumph=0.0;
02560   std::vector<UShort_t> strips;
02561 
02562   NtpSRRecord *srrec = dynamic_cast<NtpSRRecord *>(rec);
02563   NtpStRecord *strec = dynamic_cast<NtpStRecord *>(rec);
02564 
02565   Int_t *stripindex=0;
02566   Int_t nstrip=0;
02567   NtpSRSlice *slc = dynamic_cast<NtpSRSlice *>(obj);
02568   NtpSRTrack *trk = dynamic_cast<NtpSRTrack *>(obj);
02569   NtpSREvent *evt = dynamic_cast<NtpSREvent *>(obj);
02570   NtpSRShower *shw = dynamic_cast<NtpSRShower *>(obj);
02571 
02572   if(slc!=0){
02573     stripindex=slc->stp;
02574     nstrip=slc->nstrip;
02575   }
02576   else if(trk!=0){
02577     stripindex=trk->stp;
02578     nstrip=trk->nstrip;
02579   }
02580   else if(evt!=0){
02581     stripindex=evt->stp;
02582     nstrip=evt->nstrip;
02583   }
02584   else if(shw!=0){
02585     stripindex=shw->stp;
02586     nstrip=shw->nstrip;
02587   }
02588   
02589   TClonesArray *striplist=0; 
02590   int TOTSTRIPS;
02591   if(srrec!=0){
02592     striplist = srrec->stp;
02593     TOTSTRIPS=srrec->evthdr.nstrip;
02594   }
02595   else{
02596     striplist = strec->stp;
02597     TOTSTRIPS=strec->evthdr.nstrip;
02598   }
02599   
02600   //  cout<<"CHECKING nstrips="<<nstrip<<endl;
02601   if(striplist!=0){
02602     for(int i=0;i<nstrip;i++){
02603       int sindex=stripindex[i];
02604       if(sindex<0) continue;
02605       NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[sindex]);
02606       if(!strip) continue;
02607       if(strip->ph0.pe+strip->ph1.pe>phcut){
02608         strips.push_back(strip->plane);
02609         sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
02610       }     
02611     }
02612   }
02613   else{
02614     if(srrec!=0){
02615       nstrip=srrec->evthdr.nstrip;
02616     }
02617     else{
02618       nstrip=strec->evthdr.nstrip;
02619     }
02620     for(int i=0;i<nstrip;i++){
02621 
02622       NtpSRStrip *strip = dynamic_cast<NtpSRStrip *>((*striplist)[i]);
02623       //      cout<<"Checking "<<i<<" stripph: "<<strip->ph0.pe<<" "<<strip->ph1.pe<<endl;
02624       if(!strip) continue;
02625       if(strip->ph0.pe+strip->ph1.pe>phcut){
02626         strips.push_back(strip->plane);
02627         sumph+=strip->ph0.sigcor+strip->ph1.sigcor;
02628       }     
02629     }
02630   }
02631     
02632   //  cout<<"returning "<<strips.size()<<endl;
02633   return strips.size();
02634 }

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

Definition at line 145 of file MadMKAnalysis.h.

References fOnlyWriteFidEvents.

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().

02282                                                                                        {
02283 
02284   
02285   // Is the event vertex, (x,y,u,v,z) inside the fiducial volume
02286   // used by the pittsburgh group?
02287 
02288   if( !( z>0.6 && z<3.56) ) return kFALSE;
02289   if( !( u>0.3 && u<1.8) ) return kFALSE;
02290   if( !( v>-1.8 && v< -0.3) ) return kFALSE;
02291   if( x>=2.4) return kFALSE;
02292   const Float_t coil_cut=0.8*0.8;
02293   if( (x*x + y*y) <= coil_cut) return kFALSE;
02294   return kTRUE;
02295 
02296 }

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.

References MuELoss::a.

Referenced by CreatePAN(), MadEvDisplay::DrawTextBox(), and MadScanDisplay::DrawTextBox().

02299                                                            {
02300   // what is the containment status of this track 
02301   // as defined by the pitt group
02302   //
02303  // takes the track end point (x,y,u,v,z)
02304   // returns:
02305   // 1 = track is fully contained in the upstream region
02306   // 2 = track is partially contained in the upstream region
02307   // 3 = track is fully contained in the downstream region
02308   // 4 = track is partially contained in the downstream region
02309   // 
02310   // the rationale is that these categories have different momentum resolution
02311   //
02312 
02313   const Float_t radsq = x*x + y*y;
02314   int result=0;
02315   if(z<7.0) {
02316     // track in upstream region
02317     static const Float_t coil_cut=0.8*0.8;
02318     // was asking for u>3 && u<1.8 here before -- a bug
02319     if( (u>0.3 && u<1.8) && (v>-1.8 && v<-0.3) && (x<2.4) && (radsq>coil_cut) )
02320       result=1;
02321     else result=2;
02322   }
02323   else{
02324     // downstream region
02325     // uses an ellipse in x,y
02326     // major axis (x) = a = 1.7 m
02327     // minor axis (y) = b = 1.4 m
02328     // centered on x,y = (0.8,0)
02329     // (x-x0)^2/a^2 + (y-y0)^2/b^2 = 1
02330     // should change this to 0.8
02331     //    static const Float_t coil_cut=0.5*0.5; // note differs from above
02332 
02333     //changed from 0.5 to 0.8 by TV on 7-29-05
02334     static const Float_t coil_cut=0.8*0.8;
02335     static const Float_t x0=0.8;
02336     static const Float_t y0=0.0;
02337     static const Float_t a=1.7;
02338     static const Float_t b=1.4;
02339     const Float_t xsc = (x-x0)/a; // rescale ellipse to unit circle
02340     const Float_t ysc = (y-y0)/b;
02341     //MAK: cut on z was 16.0 till Aug 2, 2005
02342     if( (sqrt(xsc*xsc + ysc*ysc)<1.0) && (radsq>coil_cut) && (z<15.6)) 
02343       result = 3;
02344     else result = 4;
02345   }
02346   return result;
02347   
02348 }

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 EnergyCorrections::CorrectMomentumFromRange(), EnergyCorrections::CorrectSignedMomentumFromCurvature(), NtpSRFiducial::dr, NtpSRFiducial::dz, NtpSRTrack::fidall, MadBase::LoadTrack(), NtpSRTrack::momentum, MadBase::ntpTrack, NtpSRMomentum::qp, and NtpSRMomentum::range.

Referenced by CreatePAN().

02113                                                                                                      {
02114   const float mumass=0.10566;
02115   if(LoadTrack(itrk)){
02116      float mr=ntpTrack->momentum.range;
02117      mr=CorrectMomentumFromRange(mr,isdata,det);
02118      float mc=0.0; // signed momentum from curvature
02119      if(ntpTrack->momentum.qp!=0) mc = 1.0/(ntpTrack->momentum.qp);
02120      mc=CorrectSignedMomentumFromCurvature(mc,isdata,det);
02121 
02122     if(opt==0){  
02123       //return the most appropriate measure of momentum
02124       // assign opt based on our choice
02125       if(ntpTrack->fidall.dr>0.5&&ntpTrack->fidall.dz>0.5) {
02126         opt=2;
02127         return sqrt(mr*mr+ mumass*mumass);
02128       }
02129       else {
02130         opt=1;
02131         // in R1.9 the tracker will apparently return (q/p)=0.0
02132         // maybe it's when a track looks perfectly rigid?
02133         // if so, we have to do something
02134         // I don't want to use the range, that could be very wrong
02135         // but wrong in a more subtle way ...
02136         // so, we'll return an obviously ridiculous value of 10 TeV
02137         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
02138         else return sqrt(mc*mc+mumass*mumass);
02139       }
02140     }
02141     else if(opt==1) { //return curvature measurement
02142         if(ntpTrack->momentum.qp == 0.0) return 10000.0; 
02143         else return sqrt(mc*mc+mumass*mumass);
02144     }
02145     else if(opt==2) //return range measurement
02146       return sqrt(mr*mr + mumass*mumass);
02147     else return 0;
02148   }
02149   return 0.;
02150 }

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(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by CreatePAN().

02191                                                                    {
02192   if(!LoadTrack(itr)) return 0.;
02193 
02194   Float_t bl_z = TMath::Cos(TMath::Pi()*3./180.); //3degree beam
02195   Float_t bl_y = sqrt(1. - bl_z*bl_z);
02196   Float_t costhbl = ntpTrack->vtx.dcosz*bl_z + ntpTrack->vtx.dcosy*bl_y;
02197   
02198   return  costhbl;
02199 }

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(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by CreatePAN().

02201                                                                    {
02202   if(!LoadTrack(itr)) return 0.;
02203   /*
02204    // simple correction based on the vertical position
02205   float vtxy=0;
02206   if(vtx) vtxy=vtx[1]; // in meters
02207   // cosine of the typical neutrino angle in the yz plane
02208   // calculated as py / sqrt( py^2 + pz^2)
02209   float nu_cos = -5.799E-2;
02210   if(vtxy>-2.0 && vtxy<2.0){ //prevents further nuttyness if vtxy is silly
02211     nu_cos += vtxy*1.23304E-3 
02212       + vtxy*vtxy*1.08212E-5 
02213       + vtxy*vtxy*vtxy*(-4.634E-5);
02214   }
02215   */
02216   float nu_cos = -5.799E-2;
02217   float nu_sin = sqrt(1 -nu_cos*nu_cos);
02218   float cosz = ntpTrack->vtx.dcosz*nu_sin + ntpTrack->vtx.dcosy*nu_cos;
02219 
02220   return cosz;
02221 
02222 }

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 EnergyCorrections::CorrectShowerEnergy(), NtpSRStripPulseHeight::gev, CandShowerHandle::kCC, NtpSRShowerPulseHeight::linCCgev, MadBase::LoadShower(), MadBase::ntpShower, NtpSRShower::ph, and NtpSRShower::shwph.

Referenced by CreatePAN().

02152                                                                                                  {
02153   
02154   Float_t result = 0;
02155   // Return Jim's calculation
02156   Bool_t ok = LoadShower(opt);
02157   if(!ok) return 0;
02158 
02159   // test for shwph.linCCGeV<=0.0
02160   // if so, assume that this is an R1.16 object
02161   // and use ph.gev
02162   result = ntpShower->shwph.linCCgev;
02163   if(result<=0.0) result=(ntpShower->ph.gev/1.23);
02164   else{
02165      result = CorrectShowerEnergy(result,det,CandShowerHandle::kCC,mode,isdata);
02166   }
02167   return result;
02168   
02169 }

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

Definition at line 2462 of file MadMKAnalysis.cxx.

References fBECConfig, and Registry::Set().

02462                                            {
02463   if(f){
02464     fBECConfig.Set("beam:flux_file",f);
02465   }
02466 }

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

Definition at line 144 of file MadMKAnalysis.h.

References fDataUseMCBeam.

00144 {fDataUseMCBeam=x;}

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

Definition at line 151 of file MadMKAnalysis.h.

References fNRandomSets.

00151 {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(), MadBase::ntpTrack, and SigInOut().

02455                                                                 {
02456   // will try to load the track pointed to by trkidx
02457   if(trkidx==-1) ntpTrack=0;
02458   else if (!(LoadTrack(trkidx))) ntpTrack=0;
02459   SigInOut(sigfull,sigpart,trkfull,trkpart);
02460 }

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(), NtpSREvent::nstrip, NtpSRTrack::nstrip, MadBase::ntpEvent, MadBase::ntpStrip, MadBase::ntpTrack, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRTrack::stp, and NtpSRStrip::strip.

Referenced by CreatePAN(), and SigInOut().

02411                                                                 {
02412   // compute the amount of signal in the partially instrumented region
02413   // and the amount in the fully instrumented region
02414   //
02415   // this region is defined as:
02416   // v planes: (strip<=4 || strip>=67)
02417   // partial u: (strip==0 || strip=63)
02418   // full u: (strip<=26 || strip>=88) 
02419 
02420   sigfull=sigpart=0;
02421  
02422   // loop over all strips in the event
02423   // and sum the signals in the partial and full regions
02424   for(int i=0; i<ntpEvent->nstrip; i++){
02425     if(!LoadStrip(ntpEvent->stp[i])) continue;
02426     Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip);
02427     if(pr==1){
02428       sigpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
02429     }
02430     else if(pr==-1){
02431       sigfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
02432     }
02433   }
02434   
02435   // loop over those strips in the primary track
02436 
02437   trkfull=trkpart=0;
02438 
02439   if(ntpTrack){
02440     for(int i=0; i<ntpTrack->nstrip; i++){
02441       if(!LoadStrip(ntpTrack->stp[i])) continue;
02442       Int_t pr = InPartialRegion(ntpStrip->plane, ntpStrip->strip);
02443       if(pr==1){
02444         trkpart+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
02445       }
02446       else if(pr==-1){
02447         trkfull+=ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
02448       }
02449     }
02450   }
02451   
02452 }


Member Data Documentation

Definition at line 154 of file MadMKAnalysis.h.

Referenced by CreatePAN().

Definition at line 235 of file MadMKAnalysis.h.

Referenced by CreatePAN(), GetBECRegistry(), MadMKAnalysis(), and SetBECFile().

Definition at line 237 of file MadMKAnalysis.h.

Referenced by CreatePAN(), MadMKAnalysis(), and SetDataUseMCBeam().

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 228 of file MadMKAnalysis.h.

Referenced by EarlyActivity().

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

Definition at line 231 of file MadMKAnalysis.h.

Referenced by EarlyActivity().

Definition at line 236 of file MadMKAnalysis.h.

Referenced by CreatePAN(), ForceMCBeam(), GetForcedMCBeam(), and MadMKAnalysis().

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

Definition at line 230 of file MadMKAnalysis.h.

Referenced by EarlyActivity().

Definition at line 240 of file MadMKAnalysis.h.

Referenced by CreatePAN(), MadMKAnalysis(), and SetNRandomSets().

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 229 of file MadMKAnalysis.h.

Referenced by EarlyActivity().

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

Definition at line 220 of file MadMKAnalysis.h.

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

Definition at line 221 of file MadMKAnalysis.h.

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

Definition at line 222 of file MadMKAnalysis.h.

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

Definition at line 219 of file MadMKAnalysis.h.

Referenced by NDRadialFidVol().

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

Definition at line 218 of file MadMKAnalysis.h.

Referenced by NDRadialFidVol().

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

Definition at line 224 of file MadMKAnalysis.h.

Referenced by CreatePAN(), InFidVol(), and NDRadialFidVol().

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

Definition at line 225 of file MadMKAnalysis.h.

Referenced by CreatePAN(), InFidVol(), and NDRadialFidVol().

Definition at line 153 of file MadMKAnalysis.h.

Referenced by CreatePAN().

Definition at line 155 of file MadMKAnalysis.h.

Referenced by CreatePAN().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1