MadTVAnalysis Class Reference

#include <MadTVAnalysis.h>

Inheritance diagram for MadTVAnalysis:
MadAnalysis MadQuantities MadBase

List of all members.

Public Member Functions

void CreatePAN (std::string s, const char *bmonpath="")
 MadTVAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0)
 MadTVAnalysis (JobC *, string, int)
virtual ~MadTVAnalysis ()
Float_t ClosestPointToLine (const Float_t &x1, const Float_t &y1, const Float_t &x2, const Float_t &y2, const Float_t &xpt, const Float_t &ypt, Float_t &x, Float_t &y)
void ForceMCBeam (const BeamType::BeamType_t &beam)
const BeamType::BeamType_tGetForcedMCBeam ()
void SetDataUseMCBeam (bool x)
void FillMCHists (TFile *fout)
void SetRustemNearPID (string fname)
void SetRustemFarPID (string fname)
void SetRustemConfig (string fname)

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
MadAbID abid
bool openedabpidfile

Protected Member Functions

void SigInOut (Float_t &sigfull, Float_t &sigpart, Float_t &trkfull, Float_t &trkpart)
void SigInOut (Int_t trkidx, Float_t &sigfull, Float_t &sigpart, Float_t &trkfull, Float_t &trkpart)
virtual Float_t RecoMuDCosNeuND (Int_t itr=0, Float_t *vtx=0)
virtual Float_t RecoMuDCosNeuFD (Int_t itr=0, Float_t *vtx=0)
virtual bool InFidVol (Int_t evidx)
virtual bool InFidVol ()
bool InFidVol (float vx, float vy, float vz)
void SetBECFile (const char *f)
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)
int FindBatchNumber (double mintimediff)
float MuonShowerEnergy (const NtpSREvent *evt, const NtpSRTrack *trk)
void EarlyActivity (int evidx, float &earlywph, float &earlywphpw, float &earlywphsphere, float &ph1000, float &ph500, float &ph100, float &lateph1000, float &lateph500, float &lateph100, float &lateph50)
void DupRecoStrips (int evidx, int &nduprs, float &dupphrs)
float ComputeLowPHRatio ()
Int_t FDRCBoundary ()
Float_t GetNDCoilCurrent (const VldContext &vc)

Private Attributes

Registry fBECConfig
BeamType::BeamType_t fMCBeam
bool fDataUseMCBeam
Anp::Interfacephnti
string rustempidfile_near
string rustempidfile_far
string rustemconfigfile

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.4828
static const Float_t kYcenter = 0.2384
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 33 of file MadTVAnalysis.h.


Constructor & Destructor Documentation

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

Definition at line 60 of file MadTVAnalysis.cxx.

References MadBase::Clear(), MadBase::emrecord, fBECConfig, fDataUseMCBeam, fMCBeam, MadBase::isEM, MadBase::isMC, MadBase::isTH, BeamType::kLE, Detector::kNear, BeamType::kUnknown, MadBase::mcrecord, MadBase::Nentries, openedabpidfile, phnti, MadBase::record, Registry::Set(), MadBase::threcord, and MadBase::whichSource.

00062   : MadAnalysis(chainSR, chainMC, chainTH, chainEM)
00063 {
00064   openedabpidfile=false;
00065   fMCBeam=BeamType::kUnknown;
00066   fDataUseMCBeam=false;
00067   if(!chainSR) {
00068     record = 0;
00069     emrecord = 0;
00070     mcrecord = 0;
00071     threcord = 0;
00072     Clear();
00073     whichSource = -1;
00074     isMC = true;
00075     isTH = true;
00076     isEM = true;
00077     Nentries = -1;
00078     return;
00079   }
00080   static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
00081   fBECConfig.Set("beam:flux_file",default_bec_file);
00082   fBECConfig.Set("beam:beam_type",BeamType::kLE);
00083   fBECConfig.Set("beam:detector",Detector::kNear);
00084   fBECConfig.Set("beam:low_energy_limit",0.5);
00085   fBECConfig.Set("beam:high_energy_limit",60.0);
00086 
00087   
00088   //  InitChain(chainSR,chainMC,chainTH,chainEM);
00089   whichSource = 0;  
00090 
00091   phnti = new Anp::Interface();
00092 
00093 }

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

Definition at line 96 of file MadTVAnalysis.cxx.

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

00096                                                             : MadAnalysis(j, path, entries)
00097 {
00098   openedabpidfile=false;
00099   fMCBeam=BeamType::kUnknown;
00100   fDataUseMCBeam=false;
00101 
00102   record = 0;
00103   emrecord = 0;
00104   mcrecord = 0;
00105   threcord = 0;
00106   Clear();
00107   isMC = true;
00108   isTH = true;
00109   isEM = true;
00110   Nentries = entries;
00111   whichSource = 1;
00112   jcPath = path;
00113   fJC = j;
00114   fChain = NULL;
00115 
00116   
00117   static const char* default_bec_file="/afs/fnal.gov/files/data/minos/d04/kordosky/bec_files/pbeams_0.5_gnumi_v16.root";
00118   fBECConfig.Set("beam:flux_file",default_bec_file);
00119   fBECConfig.Set("beam:beam_type",BeamType::kLE);
00120   fBECConfig.Set("beam:detector",Detector::kNear);
00121   fBECConfig.Set("beam:low_energy_limit",0.5);
00122   fBECConfig.Set("beam:high_energy_limit",60.0);
00123 
00124   phnti = new Anp::Interface();
00125 
00126 }

MadTVAnalysis::~MadTVAnalysis (  )  [virtual]

Definition at line 129 of file MadTVAnalysis.cxx.

References phnti.

00130 {
00131   if(phnti){
00132     delete phnti;
00133     phnti=0;
00134   }
00135 }


Member Function Documentation

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

Definition at line 2227 of file MadTVAnalysis.cxx.

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

float MadTVAnalysis::ComputeLowPHRatio (  )  [protected]

Definition at line 3134 of file MadTVAnalysis.cxx.

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

Referenced by CreatePAN().

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

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

Definition at line 137 of file MadTVAnalysis.cxx.

References abid, Moments::AddPoint(), MadQEID::AltCalcQEID(), bfld::AsString(), base, BeamMonSpill::BeamType(), NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, MadDpID::CalcPID(), MadAbID::CalcPID(), NtpSRFitTrack::chi2, MadDpID::ChoosePDFs(), MadNsID::ChooseWeightFile(), NearbyVertexFinder::ClosestSpatial(), NearbyVertexFinder::ClosestTemporal(), NtpSRShower::clu, CoilTools::CoilCurrent(), NtpSRDetStatus::coilstatus, NtpTHTrack::completeslc, ComputeLowPHRatio(), Anp::Interface::Config(), NtpSRTrack::contained, NtpSRVertex::dcosx, NtpSRVertex::dcosy, det, MadBase::detStatus, MadBase::dmxStatus, dpid, DupRecoStrips(), EarlyActivity(), NtpSRPlane::end, NtpSRTrack::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, MadBase::eventSummary, fBECConfig, BeamMonSpill::fBpmInt, fDataUseMCBeam, FDRCBoundary(), BeamMonSpill::fHadInt, BeamMonSpill::fHornCur, FillMCHists(), Anp::Interface::FillSnarl(), FindBatchNumber(), NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, NtpMCFluxInfo::fluxevtno, NtpMCFluxInfo::fluxrun, fMCBeam, BeamMonSpill::fMuInt1, BeamMonSpill::fMuInt2, BeamMonSpill::fMuInt3, fname, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTargProfX, BeamMonSpill::fTargProfY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTr101d, BeamMonSpill::fTrtgtd, EnergyCorrections::FullyCorrectMomentumFromRange(), EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(), BDSpillAccessor::Get(), Registry::Get(), VldContext::GetDetector(), MadBase::GetEntry(), GetEvtTimeChar(), GetEvtTimeCharPHC(), MadNsID::GetPID(), MadAbID::GetRecoCosTheta(), MadAbID::GetRecoE(), MadAbID::GetRecoY(), RecPhysicsHeader::GetRemoteSpillType(), Moments::GetRMS(), RecDataHeader::GetRun(), VldTimeStamp::GetSeconds(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), WeightCalculator::GetStandardConfig(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), MadAbID::GetTrackEMCharge(), MadAbID::GetTrackLikePlanes(), MadAbID::GetTrackPHfrac(), MadAbID::GetTrackPHmean(), MadAbID::GetTrackPlanes(), MadAbID::GetTrackQPsigmaQP(), RecPhysicsHeader::GetTrigSrc(), Anp::Interface::GetVar(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, gSystem(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSRCluster::id, InFidVol(), MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::INu(), MadQuantities::INuNoOsc(), MadQuantities::IResonance(), ReleaseType::IsCedar(), MadQuantities::IsFid_2008(), MadQuantities::IsFidAll(), MadQuantities::IsFidFD_AB(), MadQuantities::IsFidVtxEvt(), DataUtil::IsGoodFDData(), LISieve::IsLI(), RunStatus::IsOK(), CoilTools::IsReverse(), EnergyCorrections::kDefault, Detector::kFar, BeamType::kL010z000i, BeamType::kL010z170i, BeamType::kL010z185i, BeamType::kL010z200i, BeamType::kL100z200i, BeamType::kL150z200i, BeamType::kL250z200i, BeamType::kLE, SimFlag::kMC, Detector::kNear, PlaneView::kU, BeamType::kUnknown, PlaneView::kV, kXcenter, kYcenter, NtpSRShowerPulseHeight::linCCgev, NtpSREventSummary::litime, MadBase::LoadCluster(), MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), MadBase::LoadSlice(), BDSpillAccessor::LoadSpill(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), Registry::LockValues(), ReleaseType::MakeReleaseType(), NtpSRStripPulseHeight::mip, NtpSRTrack::momentum, MuonShowerEnergy(), NtpMCFluxInfo::mupare, NtpMCFluxInfo::muparpx, NtpMCFluxInfo::muparpy, NtpMCFluxInfo::muparpz, NtpSRPlane::n, NtpSRShower::ncluster, NtpMCFluxInfo::ndecay, NtpSRFitTrack::ndof, NDRadialFidVol(), NtpMCFluxInfo::ndxdz, NtpMCFluxInfo::ndxdzfar, NtpMCFluxInfo::ndxdznear, NtpMCFluxInfo::ndydz, NtpMCFluxInfo::ndydzfar, NtpMCFluxInfo::ndydznear, NtpMCFluxInfo::necm, NtpMCFluxInfo::nenergy, NtpMCFluxInfo::nenergyfar, NtpMCFluxInfo::nenergynear, 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::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, openedabpidfile, NtpSRFitTrack::pass, MadAnalysis::PassBasicCuts(), MadQEID::PassMEDQECut(), MadQuantities::PassTrackCuts(), NtpMCFluxInfo::pdpx, NtpMCFluxInfo::pdpy, NtpMCFluxInfo::pdpz, NtpSRCluster::ph, NtpSRTrack::ph, NtpSREvent::ph, phnti, PittInFidVol(), PittTrkContained(), NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRCluster::planeview, NtpMCFluxInfo::ppdxdz, NtpMCFluxInfo::ppdydz, NtpMCFluxInfo::ppenergy, NtpMCFluxInfo::ppmedium, NtpMCFluxInfo::pppz, NtpMCFluxInfo::ppvx, NtpMCFluxInfo::ppvy, NtpMCFluxInfo::ppvz, Registry::Print(), NearbyVertexFinder::ProcessEntry(), NtpMCFluxInfo::ptype, MadQuantities::Q2(), qeid, NtpSRMomentum::qp, NtpSRMomentum::range, MadAbID::ReadPDFs(), RecoMuDCosNeuFD(), RecoMuDCosNeuND(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergyNew(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadBase::record, MadQuantities::RecoShwEnergyNew(), run(), rustemconfigfile, rustempidfile_far, rustempidfile_near, BMSpillAna::SelectSpill(), Registry::Set(), MadDpID::SetPHCorrection(), BMSpillAna::SetSpill(), BMSpillAna::SetTimeDiff(), NtpSRShower::shwph, NtpSRPulseHeight::sigcor, SigInOut(), NtpSREvent::slc, NtpSRCluster::slope, BeamMonSpill::SpillTime(), 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, BMSpillAna::UseDatabaseCuts(), 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, MadQuantities::Y(), NtpSRVertex::y, NtpMCFluxInfo::ypoint, NtpSRVertex::z, and NtpMCFluxInfo::zpoint.

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

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

Definition at line 3156 of file MadTVAnalysis.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().

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

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

Definition at line 2800 of file MadTVAnalysis.cxx.

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

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

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

Definition at line 3116 of file MadTVAnalysis.cxx.

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

Int_t MadTVAnalysis::FDRCBoundary (  )  [protected]

Definition at line 3189 of file MadTVAnalysis.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, MadBase::eventSummary, NtpSREventSummary::nshower, NtpSREventSummary::ph, NtpSREventSummary::plane, and NtpSRPulseHeight::raw.

Referenced by CreatePAN().

03189                                  {
03190   Int_t litag=0;
03191   Int_t numshower=eventSummary->nshower;
03192   Float_t allph=eventSummary->ph.raw;
03193   Int_t plbeg=eventSummary->plane.beg;
03194   Int_t plend=eventSummary->plane.end;
03195 
03196   if (numshower) {
03197     if (allph>1e6) litag+=10;
03198 
03199     if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) ||
03200         ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) ||
03201         ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) ||
03202         ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++;
03203     if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) ||
03204         ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) ||
03205         ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) ||
03206         ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++;
03207   }
03208   return litag;
03209 
03210 }

void MadTVAnalysis::FillMCHists ( TFile *  fout  ) 

Definition at line 2931 of file MadTVAnalysis.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().

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

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

Definition at line 3245 of file MadTVAnalysis.cxx.

Referenced by CreatePAN().

03246 {
03247   //cuts from peter shannahan
03248     //Batch 0: 1.7e-6 to 3.0e-6
03249     //Batch 1: 3.3e-6 to 4.5e-6
03250     //Batch 2: 4.9e-6 to 6.2e-6
03251     //Batch 3: 6.5e-6 to 7.8e-6
03252     //Batch 4: 8.1e-6 to 9.4e-6
03253     //Batch 5: 9.7e-6 to 11.0e-6
03254 
03255   //mintime already has triggertime subtracted
03256   double tdiff = mintime*1.e6;
03257   if(tdiff>1.7&&tdiff<3.0) return 0;
03258   if(tdiff>3.3&&tdiff<4.5) return 1;
03259   if(tdiff>4.9&&tdiff<6.2) return 2;
03260   if(tdiff>6.5&&tdiff<7.8) return 3;
03261   if(tdiff>8.1&&tdiff<9.4) return 4;
03262   if(tdiff>9.7&&tdiff<11.0) return 5;
03263   else return -1;
03264 }

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

Definition at line 63 of file MadTVAnalysis.h.

References fMCBeam.

00063 {fMCBeam=beam;}

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

Definition at line 113 of file MadTVAnalysis.h.

References fBECConfig.

00113 { return fBECConfig; }

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

Definition at line 2627 of file MadTVAnalysis.cxx.

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

Referenced by CreatePAN().

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

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

Definition at line 2673 of file MadTVAnalysis.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().

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

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

Definition at line 64 of file MadTVAnalysis.h.

References fMCBeam.

00064 {return fMCBeam;}

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

Definition at line 3235 of file MadTVAnalysis.cxx.

References Dcs_Mag_Near::GetCurrent(), DbiResultPtr< T >::GetNumRows(), and DbiResultPtr< T >::GetRow().

03235                                                            {
03236   Float_t result = 0.0;
03237   DbiResultPtr<Dcs_Mag_Near> ptr(vc);
03238   if(ptr.GetNumRows()){
03239     const Dcs_Mag_Near* row = ptr.GetRow(0);
03240     result = row->GetCurrent();
03241   }
03242   return result;
03243 }

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

Definition at line 2068 of file MadTVAnalysis.cxx.

References VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, kXcenter, kYcenter, and MadBase::ntpHeader.

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

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

Definition at line 2063 of file MadTVAnalysis.cxx.

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

Referenced by CreatePAN(), and InFidVol().

02064  {
02065    return InFidVol(ntpEvent->vtx.x,ntpEvent->vtx.y,ntpEvent->vtx.z);
02066  }

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

Definition at line 2058 of file MadTVAnalysis.cxx.

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

02058                                           {
02059    if(!LoadEvent(evidx)) return false; //no event found
02060    return InFidVol();
02061  }

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

Definition at line 2346 of file MadTVAnalysis.cxx.

Referenced by SigInOut().

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

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

Definition at line 2720 of file MadTVAnalysis.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().

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

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

Definition at line 2262 of file MadTVAnalysis.cxx.

References kTargEndZ, kVetoEndZ, kXcenter, and kYcenter.

Referenced by CreatePAN(), and FillMCHists().

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

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

Definition at line 2464 of file MadTVAnalysis.cxx.

References NtpStRecord::evthdr, NtpSRRecord::evthdr, NtpSREventSummary::nstrip, 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().

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

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

Definition at line 2551 of file MadTVAnalysis.cxx.

References NtpStRecord::evthdr, NtpSRRecord::evthdr, NtpSREventSummary::nstrip, 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().

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

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

Definition at line 2278 of file MadTVAnalysis.cxx.

Referenced by CreatePAN(), and FillMCHists().

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

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

Definition at line 2294 of file MadTVAnalysis.cxx.

References MuELoss::a.

Referenced by CreatePAN().

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

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

Definition at line 2187 of file MadTVAnalysis.cxx.

References NtpSRVertex::dcosy, NtpSRVertex::dcosz, MadBase::LoadTrack(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by CreatePAN().

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

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

Definition at line 2197 of file MadTVAnalysis.cxx.

References NtpSRVertex::dcosy, NtpSRVertex::dcosz, MadBase::LoadTrack(), MadBase::ntpTrack, and NtpSRTrack::vtx.

Referenced by CreatePAN().

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

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

Definition at line 2458 of file MadTVAnalysis.cxx.

References fBECConfig, and Registry::Set().

02458                                            {
02459   if(f){
02460     fBECConfig.Set("beam:flux_file",f);
02461   }
02462 }

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

Definition at line 65 of file MadTVAnalysis.h.

References fDataUseMCBeam.

00065 {fDataUseMCBeam=x;}

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

Definition at line 73 of file MadTVAnalysis.h.

References rustemconfigfile.

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

Definition at line 72 of file MadTVAnalysis.h.

References rustempidfile_far.

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

Definition at line 71 of file MadTVAnalysis.h.

References rustempidfile_near.

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

Definition at line 2450 of file MadTVAnalysis.cxx.

References MadBase::LoadTrack(), MadBase::ntpTrack, and SigInOut().

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

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

Definition at line 2406 of file MadTVAnalysis.cxx.

References InPartialRegion(), MadBase::LoadStrip(), 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().

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


Member Data Documentation

Definition at line 80 of file MadTVAnalysis.h.

Referenced by CreatePAN().

Definition at line 78 of file MadTVAnalysis.h.

Referenced by CreatePAN().

Definition at line 167 of file MadTVAnalysis.h.

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

Definition at line 169 of file MadTVAnalysis.h.

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

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

Definition at line 160 of file MadTVAnalysis.h.

Referenced by EarlyActivity().

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

Definition at line 163 of file MadTVAnalysis.h.

Referenced by EarlyActivity().

Definition at line 168 of file MadTVAnalysis.h.

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

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

Definition at line 162 of file MadTVAnalysis.h.

Referenced by EarlyActivity().

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

Definition at line 161 of file MadTVAnalysis.h.

Referenced by EarlyActivity().

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

Definition at line 152 of file MadTVAnalysis.h.

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

Definition at line 153 of file MadTVAnalysis.h.

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

Definition at line 154 of file MadTVAnalysis.h.

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

Definition at line 151 of file MadTVAnalysis.h.

Referenced by NDRadialFidVol().

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

Definition at line 150 of file MadTVAnalysis.h.

Referenced by NDRadialFidVol().

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

Definition at line 156 of file MadTVAnalysis.h.

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

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

Definition at line 157 of file MadTVAnalysis.h.

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

Definition at line 77 of file MadTVAnalysis.h.

Referenced by CreatePAN().

Definition at line 81 of file MadTVAnalysis.h.

Referenced by CreatePAN(), and MadTVAnalysis().

Definition at line 171 of file MadTVAnalysis.h.

Referenced by CreatePAN(), MadTVAnalysis(), and ~MadTVAnalysis().

Definition at line 79 of file MadTVAnalysis.h.

Referenced by CreatePAN().

Definition at line 174 of file MadTVAnalysis.h.

Referenced by CreatePAN(), and SetRustemConfig().

Definition at line 173 of file MadTVAnalysis.h.

Referenced by CreatePAN(), and SetRustemFarPID().

Definition at line 172 of file MadTVAnalysis.h.

Referenced by CreatePAN(), and SetRustemNearPID().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1