#include <MadDpAnalysis.h>
Inheritance diagram for MadDpAnalysis:
Public Member Functions | |
| MadDpAnalysis (TChain *chainSR=0, TChain *chainMC=0, TChain *chainTH=0, TChain *chainEM=0) | |
| MadDpAnalysis (JobC *, string, int) | |
| ~MadDpAnalysis () | |
| void | MakeMyFile (std::string, int, int) |
| void | ReadPIDFile (std::string) |
| void | CreatePAN (std::string, Int_t, Int_t) |
Static Public Member Functions | |
| static bool | InFidVol (const Detector::Detector_t &det, const Float_t &x, const Float_t &y, const Float_t &z) |
Protected Member Functions | |
| Bool_t | PassTruthSignal (Int_t event=0) |
| Bool_t | PassAnalysisCuts (Int_t event=0) |
| Float_t | PID (Int_t event=0, Int_t method=0) |
| Bool_t | PassBasicCuts (Int_t event=0) |
| Bool_t | PassFid (Int_t event=0) |
| Bool_t | PassFidNew (Int_t event=0) |
| Bool_t | PassFidNewN (Int_t event=0) |
| Float_t | RecoAnalysisEnergy (Int_t event=0) |
| AnnInputBlock | AnnVar (Int_t event=0) |
| Double_t | GetAnnPid (AnnInputBlock anninput, Int_t det, Int_t tar, Int_t fid, Int_t prior, Bool_t first_ann) |
| Double_t | Sigmoid (Double_t x) |
| Float_t | SingleTimeFrame (Int_t snarlentry, Int_t nentries) |
Protected Attributes | |
| TFile * | fLikeliFile |
| TH1F * | fLikeliHist [6] |
Private Attributes | |
| Double_t | Ann_weight_vector [226] |
Definition at line 8 of file MadDpAnalysis.h.
| MadDpAnalysis::MadDpAnalysis | ( | TChain * | chainSR = 0, |
|
| TChain * | chainMC = 0, |
|||
| TChain * | chainTH = 0, |
|||
| TChain * | chainEM = 0 | |||
| ) |
Definition at line 27 of file MadDpAnalysis.cxx.
References MadBase::Clear(), MadBase::emrecord, fLikeliFile, fLikeliHist, MadBase::InitChain(), MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.
00029 { 00030 00031 if(!chainSR) { 00032 record = 0; 00033 strecord = 0; 00034 emrecord = 0; 00035 mcrecord = 0; 00036 threcord = 0; 00037 Clear(); 00038 whichSource = -1; 00039 isMC = true; 00040 isTH = true; 00041 isEM = true; 00042 Nentries = -1; 00043 return; 00044 } 00045 00046 InitChain(chainSR,chainMC,chainTH,chainEM); 00047 whichSource = 0; 00048 fLikeliFile = NULL; 00049 for(int i=0;i<6;i++) fLikeliHist[i] = NULL; 00050 00051 }
| MadDpAnalysis::MadDpAnalysis | ( | JobC * | , | |
| string | , | |||
| int | ||||
| ) |
Definition at line 54 of file MadDpAnalysis.cxx.
References MadBase::Clear(), MadBase::emrecord, MadBase::fJC, fLikeliFile, fLikeliHist, MadBase::isEM, MadBase::isMC, MadBase::isTH, MadBase::jcPath, MadBase::mcrecord, MadBase::Nentries, MadBase::record, MadBase::strecord, MadBase::threcord, and MadBase::whichSource.
00055 { 00056 record = 0; 00057 strecord = 0; 00058 emrecord = 0; 00059 mcrecord = 0; 00060 threcord = 0; 00061 Clear(); 00062 isMC = true; 00063 isTH = true; 00064 isEM = true; 00065 Nentries = entries; 00066 whichSource = 1; 00067 jcPath = path; 00068 fJC = j; 00069 fLikeliFile = NULL; 00070 for(int i=0;i<6;i++) fLikeliHist[i] = NULL; 00071 }
| MadDpAnalysis::~MadDpAnalysis | ( | ) |
Definition at line 74 of file MadDpAnalysis.cxx.
References fLikeliFile.
00075 { 00076 if(fLikeliFile) fLikeliFile->Close(); 00077 }
| AnnInputBlock MadDpAnalysis::AnnVar | ( | Int_t | event = 0 |
) | [protected] |
Definition at line 172 of file MadDpAnalysis.cxx.
References AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhcommon, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwdig, AnnInputBlock::aShwph, AnnInputBlock::aShwphper, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwphperpl, AnnInputBlock::aShwphperstp, AnnInputBlock::aShwpl, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aShwstp, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrklen, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkph, AnnInputBlock::aTrkphper, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, AnnInputBlock::aTrkstp, NtpSRPlane::beg, NtpSRPlane::end, NtpSRTrack::fit, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower_Jim(), MadBase::LoadStrip(), NtpSRPlane::n, NtpSRShower::ndigit, NtpSRShower::nstrip, NtpSRTrack::nstrip, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpStrip, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRTrackPlane::ntrklike, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRFitTrack::pass, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRShower::plane, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, NtpSRShower::stp, NtpSREvent::stp, NtpSRStrip::strip, NtpSRStrip::time0, and NtpSRStrip::time1.
Referenced by CreatePAN().
00173 { 00174 AnnInputBlock anninput; 00175 // Initialize 00176 anninput.aTotrk =0.; 00177 anninput.aTotstp =0.; 00178 anninput.aTotph =0.; 00179 anninput.aTotlen =0.; 00180 anninput.aPhperpl =0.; 00181 anninput.aPhperstp =0.; 00182 00183 00184 anninput.aTrkpass =0.; 00185 anninput.aTrkph =0.; 00186 anninput.aTrklen =0.; 00187 anninput.aTrkphperpl =0.; 00188 anninput.aTrkphper =0.; 00189 anninput.aTrkplu =0.; 00190 anninput.aTrkplv =0.; 00191 anninput.aTrkstp =0.; 00192 00193 anninput.aShwph =0.; 00194 anninput.aShwstp =0.; 00195 anninput.aShwdig =0.; 00196 anninput.aShwpl =0.; 00197 anninput.aShwphper =0.; 00198 anninput.aShwphperpl =0.; 00199 anninput.aShwphperdig =0.; 00200 anninput.aShwphperstp =0.; 00201 anninput.aShwplu =0.; 00202 anninput.aShwplv =0.; 00203 anninput.aPh3 =0.; 00204 anninput.aPh6 =0.; 00205 anninput.aPhlast =0.; 00206 anninput.aPhcommon =0.; 00207 anninput.aTimemax =0.; 00208 anninput.aTimemin =0.; 00209 // 00210 00211 // Calculate variables needed for ANN (and a few more) 00212 Int_t detector = 0; 00213 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear) detector=1; 00214 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar) detector=2; 00215 00216 LoadEvent(event) ; 00217 int track_index = -1; 00218 LoadLargestTrackFromEvent(event,track_index); 00219 int shower_index = -1; 00220 // LoadLargestShowerFromEvent(event,shower_index); 00221 LoadShower_Jim(event,shower_index,detector); 00222 00223 anninput.aTotrk =ntpEvent->ntrack; 00224 anninput.aTotstp =ntpEvent->nstrip; 00225 anninput.aTotph =ntpEvent->ph.sigcor; 00226 anninput.aTotlen =ntpEvent->plane.end-ntpEvent->plane.beg+1; // prepei na to alla3w 00227 00228 if (ntpEvent->plane.n>0) anninput.aPhperpl =ntpEvent->ph.sigcor/ntpEvent->plane.n; 00229 if (ntpEvent->nstrip>0) anninput.aPhperstp =ntpEvent->ph.sigcor/ntpEvent->nstrip; 00230 00231 if(track_index!=-1) { 00232 anninput.aTrkpass =ntpTrack->fit.pass; 00233 anninput.aTrkph =ntpTrack->ph.sigcor; 00234 anninput.aTrklen =ntpTrack->plane.ntrklike; 00235 if(ntpTrack->plane.ntrklike>0) anninput.aTrkphperpl =ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike; 00236 if(ntpEvent->ph.sigcor>0) anninput.aTrkphper =ntpTrack->ph.sigcor/ntpEvent->ph.sigcor; 00237 anninput.aTrkplu =ntpTrack->plane.nu; 00238 anninput.aTrkplv =ntpTrack->plane.nv; 00239 anninput.aTrkstp =ntpTrack->nstrip; 00240 } 00241 if(shower_index!=-1) { 00242 anninput.aShwph =ntpShower->ph.sigcor; 00243 anninput.aShwstp =ntpShower->nstrip; 00244 anninput.aShwdig =ntpShower->ndigit; 00245 anninput.aShwpl =ntpShower->plane.end-ntpShower->plane.beg+1; 00246 if (ntpEvent->ph.sigcor>0) anninput.aShwphper =ntpShower->ph.sigcor/ntpEvent->ph.sigcor; 00247 anninput.aShwphperpl =ntpShower->ph.sigcor/(ntpShower->plane.end-ntpShower->plane.beg+1); 00248 00249 if (ntpShower->ndigit>0) anninput.aShwphperdig =ntpShower->ph.sigcor/ntpShower->ndigit; 00250 if (ntpShower->nstrip>0) anninput.aShwphperstp =ntpShower->ph.sigcor/ntpShower->nstrip; 00251 anninput.aShwplu =ntpShower->plane.nu; 00252 anninput.aShwplv =ntpShower->plane.nv; 00253 } 00254 00255 Int_t ssind,ssplind; 00256 Double_t ssphtot; 00257 Bool_t foundshw,foundtrk; 00258 Int_t planes=ntpEvent->plane.beg; 00259 00260 Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere; 00261 ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.; 00262 // loop over strips to compute ph3 ph6 phlast phcommon 00263 Double_t timemin=9.*10e30; 00264 Double_t timemax=-9999999; 00265 00266 for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){ 00267 Int_t stp_index=((ntpEvent->stp)[evss]); 00268 if(stp_index!=-1){ 00269 LoadStrip(stp_index); 00270 // TIMING INFO // 00271 // Find max min time of events by loading strips for ND 00272 if(detector==1){ 00273 Double_t striptime=ntpStrip->time1; 00274 if(striptime<=timemin) timemin=striptime; 00275 if(striptime>=timemax) timemax=striptime; 00276 } 00277 00278 if(detector==2){ 00279 Double_t striptime1=ntpStrip->time1; 00280 Double_t striptime0=ntpStrip->time0; 00281 Double_t striptime=0; 00282 if(striptime1>0 && striptime0<0) striptime=striptime1; 00283 if(striptime0>0 && striptime1<0) striptime=striptime0; 00284 if(striptime0>0 && striptime1>0) striptime=(striptime0+striptime1)/2.; 00285 if(striptime<=timemin) timemin=striptime; 00286 if(striptime>=timemax) timemax=striptime; 00287 00288 } 00289 00290 ssind=ntpStrip->strip; 00291 ssplind=ntpStrip->plane; 00292 ssphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; 00293 00294 foundshw=false; 00295 foundtrk=false; 00296 00297 Double_t phstrips=0; 00298 Double_t phstript=0; 00299 // shower strips 00300 Int_t sshwind,sshwplind; 00301 Double_t sshwphtot; 00302 00303 if(shower_index!=-1) { 00304 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){ 00305 Int_t stp_indexs=((ntpShower->stp)[jj]); 00306 if(stp_indexs!=-1){ 00307 LoadStrip(stp_indexs); 00308 00309 if(!foundshw){ 00310 sshwind=ntpStrip->strip; 00311 sshwplind=ntpStrip->plane; 00312 sshwphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; 00313 if(sshwind==ssind && sshwplind==ssplind) { 00314 foundshw=true; 00315 phstrips=sshwphtot; 00316 } 00317 } 00318 } 00319 } 00320 00321 } 00322 // tracks strips 00323 Int_t strkind,strkplind; 00324 Double_t strkphtot; 00325 if(track_index!=-1) { 00326 for(Int_t jj=0;jj<ntpTrack->nstrip;jj++){ 00327 Int_t stp_indext=((ntpTrack->stp)[jj]); 00328 00329 if(stp_indext!=-1){ 00330 LoadStrip(stp_indext); 00331 if(!foundtrk){ 00332 strkind=ntpStrip->strip; 00333 strkplind=ntpStrip->plane; 00334 strkphtot=ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor; 00335 if(strkind==ssind && strkplind==ssplind) { 00336 foundtrk=true; 00337 phstript=strkphtot; 00338 } 00339 } 00340 } 00341 00342 } 00343 } 00344 00345 if(foundshw && foundtrk) { 00346 hitcommon=hitcommon+1; 00347 phcommon=phcommon+phstrips+phstript; 00348 } 00349 if(!foundshw && ! foundtrk) { 00350 hitnowere=hitnowere+1; 00351 phnowere=phnowere+ssphtot; 00352 } 00353 if(ssplind>=planes && ssplind<=(planes+3)){ 00354 ph3=ph3+ssphtot; 00355 } 00356 else if(ssplind>(planes+3) && ssplind<=(planes+6)){ 00357 ph6=ph6+ssphtot; 00358 } 00359 else { 00360 phlast=phlast+ssphtot; 00361 } 00362 } // end if strip != -1 (!!!) 00363 } // end of looping over event strips.... 00364 00365 anninput.aTimemin = timemin; 00366 anninput.aTimemax = timemax; 00367 00368 // 00369 anninput.aPh3 =ph3; 00370 anninput.aPh6 =ph6; 00371 anninput.aPhlast =phlast; 00372 anninput.aPhcommon =phcommon; 00373 return anninput; 00374 }
| void MadDpAnalysis::CreatePAN | ( | std::string | , | |
| Int_t | , | |||
| Int_t | ||||
| ) |
Definition at line 944 of file MadDpAnalysis.cxx.
References AnnVar(), AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhcommon, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aPhperstp, AnnInputBlock::aShwdig, AnnInputBlock::aShwph, AnnInputBlock::aShwphper, AnnInputBlock::aShwphperdig, AnnInputBlock::aShwphperpl, AnnInputBlock::aShwphperstp, AnnInputBlock::aShwpl, AnnInputBlock::aShwplu, AnnInputBlock::aShwplv, AnnInputBlock::aShwstp, AnnInputBlock::aTimemax, AnnInputBlock::aTimemin, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTotstp, AnnInputBlock::aTrklen, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkph, AnnInputBlock::aTrkphper, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, AnnInputBlock::aTrkplv, AnnInputBlock::aTrkstp, NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, NtpSRMomentum::best, NtpSRFitTrack::chi2, NtpSRTrack::contained, NtpSREventSummary::date, NtpSRDate::day, Munits::day, NtpSRTrack::ds, NtpSRTrack::end, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRMomentum::eqp, MadBase::eventSummary, BeamMonSpill::fHornCur, NtpSRTrack::fit, MadQuantities::Flavour(), NtpMCTruth::flux, BeamMonSpill::fProfWidX, BeamMonSpill::fProfWidY, BeamMonSpill::fTargBpmX, BeamMonSpill::fTargBpmY, BeamMonSpill::fTor101, BeamMonSpill::fTortgt, BeamMonSpill::fTrtgtd, BDSpillAccessor::Get(), GetAnnPid(), ScanList::GetDecision(), VldContext::GetDetector(), SpillInfoBlock::GetHorn(), SpillInfoBlock::GetHpos(), SpillInfoBlock::GetMagnet(), VldTimeStamp::GetNanoSec(), SpillTimeFinder::GetNearestSpill(), SpillInfoBlock::GetPot(), RecDataHeader::GetRun(), VldTimeStamp::GetSec(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), SpillInfo::GetSpillInfo(), BeamMonSpill::GetStatusBits(), RecDataHeader::GetSubRun(), SpillInfoBlock::GetTgtpos(), SpillTimeFinder::GetTimeOfNearestSpill(), SpillTimeND::GetTimeStamp(), VldContext::GetTimeStamp(), SpillTimeFinder::GetTimeToNearestSpill(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), SpillInfoBlock::GetVpos(), MadQuantities::HadronicFinalState(), MadQuantities::IAction(), NtpSREvent::index, MadQuantities::Initial_state(), SpillTimeFinder::Instance(), MadQuantities::IResonance(), MadQuantities::IsFidAll(), MadQuantities::IsFidVtxEvt(), Detector::kFar, Detector::kNear, NtpSREventSummary::litime, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadShower_Jim(), BDSpillAccessor::LoadSpill(), MadBase::LoadTrack(), MadBase::LoadTruth(), MadBase::LoadTruthForRecoTH(), NtpSRTrack::momentum, NtpSRDate::month, month, NtpSRPlane::n, NtpSRShower::ncluster, NtpSRTrack::ndigit, NtpSRShower::ndigit, NtpSREvent::ndigit, NtpSRFitTrack::ndof, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, NtpSREventSummary::nshower, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSREvent::nstrip, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, MadBase::ntpTruth, NtpSREvent::ntrack, NtpSRPlane::nu, MadQuantities::Nucleus(), NtpSRPlane::nv, NtpSRFitTrack::pass, PassAnalysisCuts(), PassBasicCuts(), passfail(), PassFid(), PassFidNew(), PassFidNewN(), NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSREventSummary::ph, PID(), NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, MadQuantities::Q2(), NtpSRMomentum::qp, NtpSRMomentum::range, NtpSRPulseHeight::raw, MadQuantities::RecoMuDCosNeu(), MadQuantities::RecoMuDCosZ(), MadQuantities::RecoMuEnergy(), MadQuantities::RecoQENuEnergy(), MadQuantities::RecoQEQ2(), MadQuantities::RecoShwEnergy(), MadQuantities::RecoShwEnergySqrt(), run(), NtpSRPulseHeight::sigcor, BeamMonSpill::SpillTime(), MadQuantities::Target4Vector(), NtpMCFluxInfo::tptype, NtpMCFluxInfo::tpx, NtpMCFluxInfo::tpy, NtpMCFluxInfo::tpz, MadQuantities::TrueMuDCosNeu(), MadQuantities::TrueMuDCosZ(), MadQuantities::TrueMuEnergy(), MadQuantities::TrueNuEnergy(), MadQuantities::TrueNuMom(), MadQuantities::TrueShwEnergy(), MadQuantities::TrueVtx(), NtpSRVertex::u, NtpSRVertex::v, NtpSRTrack::vtx, NtpSRShower::vtx, NtpSREvent::vtx, MadQuantities::W2(), NtpSRVertex::x, MadQuantities::X(), NtpSRVertex::y, MadQuantities::Y(), NtpSRDate::year, Munits::year, NtpSRVertex::z, and SpillInfoBlock::Zero().
00944 { 00945 00946 // 1 = le-10 2 = pme 3 = phe there is no other way in an MC file to now what target positions it 00947 // was...! 00948 00949 std::string savename = "PAN_" + tag + ".root"; 00950 TFile save(savename.c_str(),"RECREATE"); 00951 save.cd(); 00952 /* 00953 GnumiInterface gnumi; 00954 if(!gnumi.Status()) { 00955 cout << "Error MadAnalysis::CreatePAN() - Flux files not open." 00956 << " Will not be filling NuParent info" 00957 << endl; 00958 cout << "Set environmental variable $GNUMIAUX to point to the " 00959 << "directory containing the gnumi files to fix this!" 00960 << endl; 00961 00962 } 00963 */ 00964 SpillInfo pppinfo; 00965 00966 00967 //ann quantities 00968 00969 Float_t ann_aTotrk =0.; 00970 Float_t ann_aTotstp =0.; 00971 Float_t ann_aTotph =0.; 00972 Float_t ann_aTotlen =0.; 00973 Float_t ann_aPhperpl =0.; 00974 Float_t ann_aPhperstp =0.; 00975 00976 00977 Float_t ann_aTrkpass =0.; 00978 Float_t ann_aTrkph =0.; 00979 Float_t ann_aTrklen =0.; 00980 Float_t ann_aTrkphperpl =0.; 00981 Float_t ann_aTrkphper =0.; 00982 Float_t ann_aTrkplu =0.; 00983 Float_t ann_aTrkplv =0.; 00984 Float_t ann_aTrkstp =0.; 00985 00986 Float_t ann_aShwph =0.; 00987 Float_t ann_aShwstp =0.; 00988 Float_t ann_aShwdig =0.; 00989 Float_t ann_aShwpl =0.; 00990 Float_t ann_aShwphper =0.; 00991 Float_t ann_aShwphperpl =0.; 00992 Float_t ann_aShwphperdig =0.; 00993 Float_t ann_aShwphperstp =0.; 00994 Float_t ann_aShwplu =0.; 00995 Float_t ann_aShwplv =0.; 00996 Float_t ann_aPh3 =0.; 00997 Float_t ann_aPh6 =0.; 00998 Float_t ann_aPhlast =0.; 00999 Float_t ann_aPhcommon =0.; 01000 01001 //PAN Quantities 01002 //Truth: 01003 Float_t true_enu; // true neutrino energy (GeV) 01004 Float_t true_pxnu; // true neutrino momentum-x (GeV) 01005 Float_t true_pynu; // true neutrino momentum-y (GeV) 01006 Float_t true_pznu; // true neutrino momentum-z (GeV) 01007 Float_t true_etgt; // true target energy (GeV) 01008 Float_t true_pxtgt; // true target momentum-x (GeV) 01009 Float_t true_pytgt; // true target momentum-y (GeV) 01010 Float_t true_pztgt; // true target momentum-z (GeV) 01011 Float_t true_emu; // true muon energy 01012 Float_t true_eshw; // true shower energy 01013 Float_t true_x; // true x 01014 Float_t true_y; // true y 01015 Float_t true_q2; // true q2 01016 Float_t true_w2; // true w2 01017 Float_t true_dircosneu; // true muon dircos w.r.t neutrino 01018 Float_t true_dircosz; // track z-direction cosine 01019 Float_t true_vtxx; // true x vtx 01020 Float_t true_vtxy; // true y vtx 01021 Float_t true_vtxz; // true z vtx 01022 01023 //For Neugen: 01024 Int_t flavour; // true flavour: 1-e 2-mu 3-tau 01025 Int_t process; // process: 1001-QEL 1002-RES 1003-DIS 1004-COH 01026 Int_t nucleus; // target nucleus: 274-C 284-O 372-Fe 01027 Int_t cc_nc; // cc/nc flag: 1-cc 2-nc 01028 Int_t initial_state; // initial state: 1-vp 2-vn 3-vbarp 4-vbarn ... 01029 Int_t had_fs; // hadronic final state, number between 200-300 01030 01031 //For Beam Reweighting: 01032 // NuParent *nuparent = new NuParent(); 01033 01034 //Reco Quantities 01035 Int_t detector; // Near = 1, Far = 2; 01036 Int_t run; // run number 01037 Int_t snarl; // snarl number 01038 Int_t nevent; // number of events in snarl 01039 Int_t event; // event index 01040 Int_t subrun; // subrun number 01041 Int_t trigsrc; // trigger source 01042 Int_t mcevent; // mc event index 01043 Int_t ntrack; // number of tracks in event 01044 Int_t nshower; // number of showers in event 01045 01046 Int_t is_fid; // pass fid cut. true = 1, false = 0 01047 01048 Int_t is_fid_dp; // pass dp analysis fid cut. true = 1, false = 0 01049 Int_t is_fid_ns; // pass nd analysis fid cut true=1, false=0 // !! NEW 01050 01051 Int_t is_cev; // fully contained. true = 1, false = 0 01052 Int_t is_cont_trk; // new containemnt support by IsContained() for the largest track 01053 Int_t is_mc; // is a corresponding mc event found 01054 Int_t pass_basic; // Pass basic cuts. true = 1, false = 0 01055 Float_t pid0; // pid parameter (usual method) 01056 Float_t pid1; // pid parameter (alternative method 1) 01057 Float_t pid2; // pid parameter (alternative method 2) 01058 Float_t annpid; // ANN PID 01059 Float_t annpid_f1p1; // ANN PID FAR FIDUCIAL DP PRIOR 1 // UN-OSICLLATED 01060 Float_t annpid_f1p2; // ANN PID FAR FIDUCIAL DP PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95 01061 Float_t annpid_f1p3; // ANN PID FAR FIDUCIAL DP PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95 01062 Float_t annpid_f2p1; // ANN PID FAR FIDUCIAL NS PRIOR 1 UN-OSICLLATED 01063 Float_t annpid_f2p2; // ANN PID FAR FIDUCIAL NS PRIOR 2 // OSCILLATED WITH DM2 0.0025 SIN2THETA 0.95 01064 Float_t annpid_f2p3; // ANN PID FAR FIDUCIAL NS PRIOR 3 // OSCILLATED WITH DM2 0.0020 SIN2THETA 0.95 01065 01066 Int_t pass; // pass analysis cuts. true = 1, false = 0 01067 01068 Float_t reco_enu; // reco neutrino energy 01069 Float_t reco_emu; // reco muon energy 01070 Float_t reco_eshw; // reco shower energy (shw.ph.gev/1.23) 01071 01072 Int_t pass_fd_qualcuts; // pass FD quality cuts (HV trips etc) 01073 Int_t litag; 01074 01075 01076 Float_t reco_eshwcc; // linear cc version 01077 Float_t reco_eshwnc; // linear nc version 01078 Float_t reco_eshwccw; // weighted cc version 01079 Float_t reco_eshwncw; // weighted nc version 01080 01081 Float_t reco_eshw_sqrt; // reco shower energy using sqrt method 01082 Float_t reco_qe_enu; // reco qe neutrino energy 01083 Float_t reco_qe_q2; // reco qe q2 01084 Float_t reco_y; // reco y 01085 Float_t reco_dircosneu; // reco track vtx dircos w.r.t neutrino 01086 Float_t reco_dircosz; // reco track vtx z-dircos 01087 Float_t reco_vtxx; // reco x vtx 01088 Float_t reco_vtxy; // reco y vtx 01089 Float_t reco_vtxz; // reco z vtx 01090 01091 Float_t evtlength; // reco event length 01092 Float_t evtph; // reco event ph 01093 Float_t trklength; // reco track length 01094 Float_t trkmomrange; // reco track momentum from range 01095 Float_t trkqp; // reco track q/p 01096 Float_t trkeqp; // reco track q/p error 01097 Float_t trkmombest; // reco best track momentum based on new containment support 01098 Float_t trkphfrac; // reco pulse height fraction in track 01099 Float_t trkphplane; // reco average track pulse height/plane 01100 Float_t trkds; // track DS 01101 Float_t rawph; // Raw ph 01102 // Additional track info 01103 Float_t trkendz; // track end Z position 01104 Float_t trkendx; // track end X position 01105 Float_t trkendy; // track end Y position 01106 Float_t trkendu; // track end U position 01107 Float_t trkendv; // track end V position 01108 Float_t trkvtxz; // track begin Z position 01109 Float_t trkvtxx; // track begin X position 01110 Float_t trkvtxy; // track begin Y position 01111 Float_t trkvtxu; // track begin U position 01112 Float_t trkvtxv; // track begin V position 01113 Float_t trkplbu; // track begin plane u 01114 Float_t trkplbv; // track begin plane v 01115 Float_t trkpleu; // track end plane u 01116 Float_t trkplev; // track end plane v 01117 Float_t trkple; // track end plane 01118 Float_t trkplb; // track begin plane 01119 Int_t trkfitpass; // 01120 Int_t totdigits,totstrips; 01121 Int_t trkdigits,trkstrips; 01122 Float_t trkph,trkchi2,trkndof; 01123 Int_t trkdiffuv; 01124 // shower variables 01125 Float_t shwph; 01126 Int_t shwdigits, shwstrips; 01127 Float_t shwvtxz; // shower begin Z position 01128 Float_t shwvtxx; // shower begin X position 01129 Float_t shwvtxy; // shower begin Y position 01130 Float_t shwvtxu; // shower begin U position 01131 Float_t shwvtxv; // shower begin V position 01132 01133 // cluster variables 01134 01135 Int_t shwncluster; 01136 01137 // SCAN DECISION // <- define an instance of ScanList class here 01138 Int_t scandecision; 01139 ScanList scan; 01140 01141 // EVENT TIMING 01142 Double_t evttimemin; //MIN STRIP TIME OF EVENT 01143 Double_t evttimemax; //MAX STRIP TIME OF EVENT 01144 Double_t snarlt; 01145 Double_t litime; 01146 01147 // BEAM INFO 01148 Float_t pot,horn,tar,vvpos,hhpos,magn; 01149 // Beam INFO DATABASE 01150 Float_t potdb,horndb,tardb,vvposdb,hhposdb,vvwidthdb,hhwidthdb; 01151 Double_t timedb; 01152 // TIME INFO DATABASE 01153 Double_t neartdb,fartdb,neardifftdb; 01154 01155 // mc flux info for xf,pt weighting 01156 Float_t mcflux_tpx; 01157 Float_t mcflux_tpy; 01158 Float_t mcflux_tpz; 01159 Float_t mcflux_tptype; 01160 01161 // IS SNARL CUT IN PIECES 01162 Float_t singletimeframe; 01163 // Beam Info Block 01164 SpillInfoBlock *spinfo = new SpillInfoBlock(); 01165 spinfo->Zero(); 01166 singletimeframe=-999; 01167 // 01168 Int_t day,month,year; 01169 Int_t pass_beamcuts; // pass default beam quality cuts 01170 //Trees 01171 TTree *tree = new TTree("pan","pan"); 01172 //Truth 01173 tree->Branch("true_enu",&true_enu,"true_enu/F",512000); 01174 tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",512000); 01175 tree->Branch("true_pynu",&true_pynu,"true_pynu/F",512000); 01176 tree->Branch("true_pznu",&true_pznu,"true_pznu/F",512000); 01177 tree->Branch("true_etgt",&true_etgt,"true_etgt/F",512000); 01178 tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",512000); 01179 tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",512000); 01180 tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",512000); 01181 tree->Branch("true_emu",&true_emu,"true_emu/F",512000); 01182 tree->Branch("true_eshw",&true_eshw,"true_eshw/F",512000); 01183 tree->Branch("true_x",&true_x,"true_x/F",512000); 01184 tree->Branch("true_y",&true_y,"true_y/F",512000); 01185 tree->Branch("true_q2",&true_q2,"true_q2/F",512000); 01186 tree->Branch("true_w2",&true_w2,"true_w2/F",512000); 01187 tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",512000); 01188 tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",512000); 01189 tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",512000); 01190 tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",512000); 01191 tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",512000); 01192 //Neugen 01193 tree->Branch("flavour",&flavour,"flavour/I",512000); 01194 tree->Branch("process",&process,"process/I",512000); 01195 tree->Branch("nucleus",&nucleus,"nucleus/I",512000); 01196 tree->Branch("cc_nc",&cc_nc,"cc_nc/I",512000); 01197 tree->Branch("initial_state",&initial_state,"initial_state/I",512000); 01198 tree->Branch("had_fs",&had_fs,"had_fs/I",512000); 01199 //NuParent 01200 //tree->Branch("nuparent","NuParent",&nuparent,8000,1); 01201 //Reco 01202 tree->Branch("detector",&detector,"detector/I",512000); 01203 tree->Branch("run",&run,"run/I",512000); 01204 tree->Branch("snarl",&snarl,"snarl/I",512000); 01205 tree->Branch("event",&event,"event/I",512000); 01206 tree->Branch("mcevent",&mcevent,"mcevent/I",512000); 01207 tree->Branch("ntrack",&ntrack,"ntrack/I",512000); 01208 tree->Branch("nshower",&nshower,"nshower/I",512000); 01209 tree->Branch("subrun",&subrun,"subrun/I",512000); 01210 tree->Branch("trigsrc",&trigsrc,"trigsrc/I",512000); 01211 tree->Branch("litime",&litime,"litime/D",512000); 01212 tree->Branch("is_fid",&is_fid,"is_fid/I",512000); 01213 tree->Branch("is_fid_dp",&is_fid_dp,"is_fid_dp/I",512000); 01214 tree->Branch("is_fid_ns",&is_fid_ns,"is_fid_ns/I",512000); 01215 01216 01217 tree->Branch("pass_fd_qualcuts", &pass_fd_qualcuts, "pass_fd_qualcuts/I"); 01218 tree->Branch("litag", &litag, "litag/I"); 01219 01220 tree->Branch("is_cev",&is_cev,"is_cev/I",512000); 01221 tree->Branch("is_cont_trk",&is_cont_trk,"is_cont_trk/I",512000); 01222 tree->Branch("is_mc",&is_mc,"is_mc/I",512000); 01223 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",512000); 01224 tree->Branch("pid0",&pid0,"pid0/F",512000); 01225 tree->Branch("pid1",&pid1,"pid1/F",512000); 01226 tree->Branch("pid2",&pid2,"pid2/F",512000); 01227 tree->Branch("annpid",&annpid,"annpid/F",512000); 01228 tree->Branch("annpid_f1p1",&annpid_f1p1,"annpid_f1p1/F",512000); 01229 tree->Branch("annpid_f1p2",&annpid_f1p2,"annpid_f1p2/F",512000); 01230 tree->Branch("annpid_f1p3",&annpid_f1p3,"annpid_f1p3/F",512000); 01231 tree->Branch("annpid_f2p1",&annpid_f2p1,"annpid_f2p1/F",512000); 01232 tree->Branch("annpid_f2p2",&annpid_f2p2,"annpid_f2p2/F",512000); 01233 tree->Branch("annpid_f2p3",&annpid_f2p3,"annpid_f2p3/F",512000); 01234 tree->Branch("pass",&pass,"pass/I",512000); 01235 01236 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",512000); 01237 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",512000); 01238 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",512000); 01239 tree->Branch("reco_eshw_sqrt",&reco_eshw_sqrt,"reco_eshw_sqrt/F",512000); 01240 tree->Branch("reco_eshwcc", &reco_eshwcc, "reco_eshwcc/F",512000); 01241 tree->Branch("reco_eshwnc", &reco_eshwnc, "reco_eshwnc/F",512000); 01242 tree->Branch("reco_eshwccw", &reco_eshwccw, "reco_eshwccw/F",512000); 01243 tree->Branch("reco_eshwncw", &reco_eshwncw, "reco_eshwncw/F",512000); 01244 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",512000); 01245 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",512000); 01246 tree->Branch("reco_y",&reco_y,"reco_y/F",512000); 01247 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",512000); 01248 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",512000); 01249 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",512000); 01250 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",512000); 01251 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",512000); 01252 01253 tree->Branch("evtlength",&evtlength,"evtlength/F",512000); 01254 tree->Branch("evtph",&evtph,"evtph/F",512000); 01255 tree->Branch("trklength",&trklength,"trklength/F",512000); 01256 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",512000); 01257 tree->Branch("trkqp",&trkqp,"trkqp/F",512000); 01258 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",512000); 01259 tree->Branch("trkmombest",&trkmombest,"trkmombest/F",512000); 01260 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",512000); 01261 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",512000); 01262 tree->Branch("rawph",&rawph,"rawph/F",512000); 01263 01264 tree->Branch("trkendz",&trkendz,"trkendz/F",512000); 01265 tree->Branch("trkendu",&trkendu,"trkendu/F",512000); 01266 tree->Branch("trkendv",&trkendv,"trkendv/F",512000); 01267 tree->Branch("trkendx",&trkendx,"trkendx/F",512000); 01268 tree->Branch("trkendy",&trkendy,"trkendy/F",512000); 01269 tree->Branch("trkplbu",&trkplbu,"trkplbu/F",512000); 01270 tree->Branch("trkplbv",&trkplbv,"trkplbv/F",512000); 01271 tree->Branch("trkpleu",&trkpleu,"trkpleu/F",512000); 01272 tree->Branch("trkplev",&trkplev,"trkplev/F",512000); 01273 tree->Branch("trkple",&trkple,"trkple/F",512000); 01274 tree->Branch("trkplb",&trkplb,"trkplb/F",512000); 01275 01276 tree->Branch("trkvtxz",&trkvtxz,"trkvtxz/F",512000); 01277 tree->Branch("trkvtxu",&trkvtxu,"trkvtxu/F",512000); 01278 tree->Branch("trkvtxv",&trkvtxv,"trkvtxv/F",512000); 01279 tree->Branch("trkvtxx",&trkvtxx,"trkvtxx/F",512000); 01280 tree->Branch("trkvtxy",&trkvtxy,"trkvtxy/F",512000); 01281 tree->Branch("trkds", &trkds,"trkds/F",512000); 01282 01283 tree->Branch("evttimemin",&evttimemin,"evttimemin/D"); 01284 tree->Branch("evttimemax",&evttimemax,"evttimemax/D"); 01285 01286 tree->Branch("totdigits",&totdigits,"totdigits/I"); 01287 tree->Branch("totstrips",&totstrips,"totstrips/I"); 01288 tree->Branch("trkdigits",&trkdigits,"trkdigits/I"); 01289 tree->Branch("trkstrips",&trkstrips,"trkstrips/I"); 01290 tree->Branch("trkph",&trkph,"trkph/F"); 01291 tree->Branch("trkdiffuv",&trkdiffuv,"trkdiffuv/I",512000); 01292 01293 tree->Branch("trkchi2",&trkchi2,"trkchi2/F"); 01294 tree->Branch("trkndof",&trkndof,"trkndof/F"); 01295 tree->Branch("trkfitpass",&trkfitpass,"trkfitpass/I",512000); 01296 01297 tree->Branch("shwdigits",&shwdigits,"shwdigits/I"); 01298 tree->Branch("shwstrips",&shwstrips,"shwstrips/I"); 01299 tree->Branch("shwph",&shwph,"shwph/F"); 01300 01301 tree->Branch("shwvtxz",&shwvtxz,"shwvtxz/F",512000); 01302 tree->Branch("shwvtxu",&shwvtxu,"shwvtxu/F",512000); 01303 tree->Branch("shwvtxv",&shwvtxv,"shwvtxv/F",512000); 01304 tree->Branch("shwvtxx",&shwvtxx,"shwvtxx/F",512000); 01305 tree->Branch("shwvtxy",&shwvtxy,"shwvtxy/F",512000); 01306 tree->Branch("shwncluster",&shwncluster,"shwncluster/I"); 01307 01308 // SCAN DECISIONS 01309 tree->Branch("scandecision",&scandecision,"scandecision/I"); 01310 01311 //BEAM INFO 01312 tree->Branch("pot", &pot, "pot/F",512000); 01313 tree->Branch("horn", &horn, "horn/F",512000); 01314 tree->Branch("tar", &tar, "tar/F",512000); 01315 tree->Branch("vvpos", &vvpos, "vvpos/F",512000); 01316 tree->Branch("hhpos", &hhpos, "hhpos/F",512000); 01317 tree->Branch("magn", &magn, "magn/F",512000); 01318 tree->Branch("singletimeframe",&singletimeframe,"singletimeframe/F",512000); 01319 01320 // BEAM INFO DB 01321 tree->Branch("potdb", &potdb, "potdb/F",512000); 01322 tree->Branch("horndb", &horndb, "horndb/F",512000); 01323 tree->Branch("tardb", &tardb, "tardb/F",512000); 01324 tree->Branch("vvposdb", &vvposdb, "vvposdb/F",512000); 01325 tree->Branch("hhposdb", &hhposdb, "hhposdb/F",512000); 01326 tree->Branch("vvwidthdb", &vvwidthdb, "vvwidthdb/F",512000); 01327 tree->Branch("hhwidthdb", &hhwidthdb, "hhwidthdb/F",512000); 01328 tree->Branch("timedb", &timedb, "timedb/D"); 01329 01330 // TIME DB 01331 tree->Branch("neartdb", &neartdb, "neartdb/D"); 01332 tree->Branch("fartdb", &fartdb, "fartdb/D"); 01333 tree->Branch("neardifftdb", &neardifftdb,"neardifftdb/D"); 01334 tree->Branch("snarlt", &snarlt, "snarlt/D"); 01335 // mc flux info 01336 01337 tree->Branch("mcflux_tpx", &mcflux_tpx, "mcflux_tpx/F"); 01338 tree->Branch("mcflux_tpy", &mcflux_tpy, "mcflux_tpy/F"); 01339 tree->Branch("mcflux_tpz", &mcflux_tpz, "mcflux_tpz/F"); 01340 tree->Branch("mcflux_tptype", &mcflux_tptype, "mcflux_tptype/F"); 01341 01342 // ann vars 01343 01344 tree->Branch("ann_aTotrk", &ann_aTotrk, "ann_aTotrk/F"); 01345 tree->Branch("ann_aTotstp", &ann_aTotstp, "ann_aTotstp/F"); 01346 tree->Branch("ann_aTotph", &ann_aTotph, "ann_aTotph/F"); 01347 tree->Branch("ann_aTotlen", &ann_aTotlen, "ann_aTotlen/F"); 01348 tree->Branch("ann_aPhperpl", &ann_aPhperpl, "ann_aPhperpl/F"); 01349 tree->Branch("ann_aPhperstp", &ann_aPhperstp, "ann_aPhperstp/F"); 01350 tree->Branch("ann_aTrkpass", &ann_aTrkpass, "ann_aTrkpass/F"); 01351 tree->Branch("ann_aTrkph", &ann_aTrkph, "ann_aTrkph/F"); 01352 tree->Branch("ann_aTrklen", &ann_aTrklen, "ann_aTrklen/F"); 01353 tree->Branch("ann_aTrkphperpl", &ann_aTrkphperpl, "ann_aTrkphperpl/F"); 01354 tree->Branch("ann_aTrkphper", &ann_aTrkphper, "ann_aTrkphper/F"); 01355 tree->Branch("ann_aTrkplu", &ann_aTrkplu, "ann_aTrkplu/F"); 01356 tree->Branch("ann_aTrkplv", &ann_aTrkplv, "ann_aTrkplv/F"); 01357 tree->Branch("ann_aTrkstp", &ann_aTrkstp, "ann_aTrkstp/F"); 01358 tree->Branch("ann_aShwph", &ann_aShwph, "ann_aShwph/F"); 01359 tree->Branch("ann_aShwstp", &ann_aShwstp, "ann_aShwstp/F"); 01360 tree->Branch("ann_aShwdig", &ann_aShwdig, "ann_aShwdig/F"); 01361 tree->Branch("ann_aShwpl", &ann_aShwpl, "ann_aShwpl/F"); 01362 tree->Branch("ann_aShwphper", &ann_aShwphper, "ann_aShwphper/F"); 01363 tree->Branch("ann_aShwphperpl", &ann_aShwphperpl, "ann_aShwphperpl/F"); 01364 tree->Branch("ann_aShwphperdig", &ann_aShwphperdig,"ann_aShwphperdig/F"); 01365 tree->Branch("ann_aShwphperstp", &ann_aShwphperstp,"ann_aShwphperstp/F"); 01366 tree->Branch("ann_aShwplu", &ann_aShwplu, "ann_aShwplu/F"); 01367 tree->Branch("ann_aShwplv", &ann_aShwplv, "ann_aShwplv/F"); 01368 tree->Branch("ann_aPh3", &ann_aPh3, "ann_aPh3/F"); 01369 tree->Branch("ann_aPh6", &ann_aPh6, "ann_aPh6/F"); 01370 tree->Branch("ann_aPhlast", &ann_aPhlast, "ann_aPhlast/F"); 01371 tree->Branch("ann_aPhcommon", &ann_aPhcommon, "ann_aPhcommon/F"); 01372 01373 tree->Branch("day", &day, "day/I"); 01374 tree->Branch("month", &month, "month/I"); 01375 tree->Branch("year", &year, "year/I"); 01376 tree->Branch("pass_beamcuts", &pass_beamcuts, "pass_beamcuts/I"); 01377 Bool_t first_ann = true; 01378 for(int i=0;i<Nentries;i++){ 01379 01380 if(i%400==0) std::cout << float(i)*100./float(Nentries) 01381 << "% done" << std::endl; 01382 01383 if(!this->GetEntry(i)) continue; 01384 01385 nevent = eventSummary->nevent; 01386 if(nevent==0) continue; 01387 01388 run = ntpHeader->GetRun(); 01389 snarl = ntpHeader->GetSnarl(); 01390 subrun = ntpHeader->GetSubRun(); 01391 trigsrc = ntpHeader->GetTrigSrc(); 01392 day = eventSummary->date.day; 01393 month = eventSummary->date.month; 01394 year = eventSummary->date.year; 01395 01396 Double_t snarltime =ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9); 01397 Int_t month = eventSummary->date.month; 01398 litime = eventSummary->litime; 01399 01400 pass_fd_qualcuts =passfail(snarltime); 01401 01402 litag=0; 01403 Int_t numshower=eventSummary->nshower; 01404 Float_t allph=eventSummary->ph.raw; 01405 Int_t plbeg=eventSummary->plane.beg; 01406 Int_t plend=eventSummary->plane.end; 01407 01408 if (numshower) { 01409 if (allph>1e6) litag+=10; 01410 01411 if (((plbeg==1 || plbeg==2) && (plend==63 || plend==64)) || 01412 ((plbeg==65 || plbeg==66) && (plend==127 || plend==128)) || 01413 ((plbeg==129 || plbeg==130) && (plend==191 || plend==192)) || 01414 ((plbeg==193 || plbeg==194) && (plend==247 || plend==248))) litag++; 01415 if (((plbeg==250 || plbeg==251) && (plend==312 || plend==313)) || 01416 ((plbeg==314 || plbeg==315) && (plend==376 || plend==377)) || 01417 ((plbeg==378 || plbeg==379) && (plend==440 || plend==441)) || 01418 ((plbeg==442 || plbeg==443) && (plend==484 || plend==485))) litag++; 01419 } 01420 01421 01422 01423 // SCAN INFO 01424 01425 scandecision=scan.GetDecision(run,snarl); 01426 01427 // only pass through scanned snarls if scanfilter=1 01428 if (scanfilter==0 || scandecision>0) { 01429 01430 // Get Beam Monitorint Info if NOT MC: 01431 01432 if(ntpHeader->GetVldContext().GetSimFlag()!=4){ 01433 pppinfo.GetSpillInfo(month,year,snarltime,*spinfo); 01434 01435 pot =spinfo->GetPot(); 01436 horn =spinfo->GetHorn(); 01437 hhpos =spinfo->GetHpos(); 01438 vvpos =spinfo->GetVpos(); 01439 tar =spinfo->GetTgtpos(); 01440 magn =spinfo->GetMagnet(); 01441 01442 potdb = -999; 01443 horndb = -999; 01444 tardb = -999; 01445 vvposdb = -999; 01446 hhposdb = -999; 01447 vvwidthdb = -999; 01448 hhwidthdb = -999; 01449 timedb = -999; 01450 neartdb = -999; 01451 fartdb = -999; 01452 neardifftdb = -999; 01453 01454 // BEAM INFO DB 01455 Detector::Detector_t detectort; 01456 SimFlag::SimFlag_t simflag; 01457 VldTimeStamp timestamp; 01458 VldContext cx; 01459 cx = ntpHeader->GetVldContext(); 01460 detectort = ntpHeader->GetVldContext().GetDetector(); 01461 simflag = ntpHeader->GetVldContext().GetSimFlag(); 01462 timestamp = ntpHeader->GetVldContext().GetTimeStamp(); 01463 01464 01465 BDSpillAccessor &spillAccess=BDSpillAccessor::Get(); 01466 const BeamMonSpill *spillInfoDB = spillAccess.LoadSpill(timestamp); 01467 if(spillInfoDB) { 01468 if(spillInfoDB->fTortgt !=0.) potdb = spillInfoDB->fTortgt; 01469 else if(spillInfoDB->fTrtgtd !=0.) potdb = spillInfoDB->fTrtgtd; 01470 else if(spillInfoDB->fTor101 !=0.) potdb = spillInfoDB->fTor101; 01471 else potdb = -999; 01472 01473 horndb = spillInfoDB->fHornCur; 01474 tardb = (int)spillInfoDB->GetStatusBits().target_in; 01475 vvposdb = spillInfoDB->fTargBpmY[0]; 01476 hhposdb = spillInfoDB->fTargBpmX[0]; 01477 vvwidthdb = spillInfoDB->fProfWidY; 01478 hhwidthdb = spillInfoDB->fProfWidX; 01479 timedb = spillInfoDB->SpillTime(); 01480 // apply beam quality cuts 01481 01482 if (potdb>0.5 && fabs(snarltime-timedb)<1 && 01483 (hhposdb<0 && hhposdb>-0.002) && 01484 (vvposdb>0 && vvposdb<0.002) && 01485 (hhwidthdb>0.0001 && hhwidthdb<0.0022) && 01486 (vvwidthdb>0.0001 && vvwidthdb<0.005) && 01487 (horndb<-155 && horndb>-200)) pass_beamcuts=1; 01488 01489 } 01490 01491 01492 // TIME DB 01493 SpillTimeFinder& spillfinder = SpillTimeFinder::Instance(); 01494 fartdb = spillfinder.GetTimeOfNearestSpill(cx); 01495 const SpillTimeND& spnd = spillfinder.GetNearestSpill(cx); 01496 neartdb = spnd.GetTimeStamp().GetSec()+(spnd.GetTimeStamp().GetNanoSec()/1e9); 01497 neardifftdb = spillfinder.GetTimeToNearestSpill(cx); 01498 snarlt = ntpHeader->GetVldContext().GetTimeStamp().GetSec()+(ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9); 01499 01500 } 01501 01502 //singletimeframe=SingleTimeFrame(i,Nentries); 01503 01504 Int_t trkcount=0; 01505 Int_t shwcount=0; 01506 Float_t totph=0; 01507 01508 Int_t evtmin=0; 01509 Int_t evtmax=nevent; 01510 01511 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) { 01512 // cout << "MULTIPLE events!!\n"; 01513 01514 Int_t maxevent=0; 01515 Float_t maxph=0; 01516 for(int j=0;j<nevent;j++){ 01517 if(!LoadEvent(j)) continue; 01518 // cout << "evtlength=" << ntpEvent->plane.n << " totph=" 01519 // << ntpEvent->ph.sigcor << "ntrack=" << ntpEvent->ntrack 01520 // << " nshower=" << ntpEvent->nshower << endl; 01521 Float_t evtpheight=ntpEvent->ph.sigcor; 01522 if (evtpheight>maxph) { 01523 maxph=evtpheight; 01524 maxevent=j; 01525 } 01526 } 01527 // cout << "BIGGEST event=" << maxevent << endl; 01528 evtmin=maxevent; 01529 evtmax=maxevent+1; 01530 } 01531 01532 01533 // EVENT LOOP - fill tree once for each reconstructed event 01534 01535 for(int j=evtmin;j<evtmax;j++){ 01536 01537 if(!LoadEvent(j)) continue; //no event found 01538 01539 totph = 0; 01540 01541 //set event index: 01542 event = j; 01543 ntrack = ntpEvent->ntrack; 01544 nshower = ntpEvent->nshower; 01545 totdigits=ntpEvent->ndigit; 01546 totstrips=ntpEvent->nstrip; 01547 01548 // Initialize everything 01549 //nuparent->Zero(); 01550 mcflux_tpx=0; mcflux_tpy=0; mcflux_tpz=0; mcflux_tptype=0; 01551 true_enu = 0; true_emu = 0; true_eshw = 0; 01552 true_pxnu = 0; true_pynu = 0; true_pznu = 0; 01553 true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0; 01554 true_dircosneu = 0; true_dircosz = 0; 01555 true_vtxx = 0; true_vtxy = 0; true_vtxz = 0; 01556 flavour = 0; process = 0; nucleus = 0; cc_nc = 0; 01557 initial_state = 0; had_fs = 0; 01558 true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0; 01559 detector = -1; mcevent = -1; 01560 is_fid = 0; is_cev = 0; is_mc = 0; pass_basic = 0; 01561 pid0 = -999; pid1 = -999; pid2 = -999; pass = 0; 01562 reco_enu = 0; reco_emu = 0; reco_eshw = 0; reco_eshw_sqrt = 0; 01563 reco_qe_enu = 0; reco_qe_q2 = 0; reco_y = 0; reco_dircosneu = 0; 01564 reco_dircosz = 0; 01565 reco_vtxx = -999; reco_vtxy = -999; reco_vtxz = -999; 01566 evtlength = -1; trklength = -1; trkphfrac = -1; trkphplane = -1; 01567 trkmomrange = -999; trkqp = -999; trkeqp = -999; 01568 trkendz=100.; trkendu=100.; trkendv=100.; trkendx=100.; trkendy=100.; 01569 trkendz=999; trkendu=999; trkendv=999; trkendx=999; 01570 trkendy=-999;trkvtxz=-999; trkvtxu=-999; trkvtxv=-999; trkvtxx=-999; trkvtxy=-999; 01571 trkds=-1; evttimemin=-1; evttimemax=-1; 01572 evtph=0.; 01573 annpid=-999.; 01574 annpid_f1p1=-999.;annpid_f1p2=-999.;annpid_f1p3=-999.; 01575 annpid_f2p1=-999.;annpid_f2p2=-999.;annpid_f2p3=-999.; 01576 reco_eshwcc=0.;reco_eshwnc=0.;reco_eshwccw=0.;reco_eshwncw=0.; 01577 is_fid_dp=0;is_fid_ns=0; 01578 is_cont_trk=0; trkmombest=-999.; 01579 01580 trkfitpass=0; 01581 trkdiffuv=0; 01582 trkdigits=0; trkstrips=0; trkph=0; trkndof=0; trkchi2=0; 01583 shwdigits=0; shwstrips=0; shwph=0; 01584 shwvtxu=0; 01585 shwvtxv=0; 01586 shwvtxx=0; 01587 shwvtxy=0; 01588 shwvtxz=0; 01589 shwncluster=0; 01590 01591 01592 ann_aTotrk =0.; 01593 ann_aTotstp =0.; 01594 ann_aTotph =0.; 01595 ann_aTotlen =0.; 01596 ann_aPhperpl =0.; 01597 ann_aPhperstp =0.; 01598 01599 01600 ann_aTrkpass =0.; 01601 ann_aTrkph =0.; 01602 ann_aTrklen =0.; 01603 ann_aTrkphperpl =0.; 01604 ann_aTrkphper =0.; 01605 ann_aTrkplu =0.; 01606 ann_aTrkplv =0.; 01607 ann_aTrkstp =0.; 01608 01609 ann_aShwph =0.; 01610 ann_aShwstp =0.; 01611 ann_aShwdig =0.; 01612 ann_aShwpl =0.; 01613 ann_aShwphper =0.; 01614 ann_aShwphperpl =0.; 01615 ann_aShwphperdig =0.; 01616 ann_aShwphperstp =0.; 01617 ann_aShwplu =0.; 01618 ann_aShwplv =0.; 01619 ann_aPh3 =0.; 01620 ann_aPh6 =0.; 01621 ann_aPhlast =0.; 01622 ann_aPhcommon =0.; 01623 01624 rawph=ntpEvent->ph.raw; 01625 //get detector type: 01626 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ 01627 detector = 2; 01628 totph=eventSummary->ph.sigcor; 01629 //fiducial criteria on vtx for far detector 01630 is_fid = 1; 01631 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0; 01632 } 01633 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ 01634 detector = 1; 01635 01636 //get total charge in trk/shw 01637 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) { 01638 LoadTrack(ii); 01639 totph+=ntpTrack->ph.sigcor; 01640 } 01641 trkcount+=ntrack; 01642 01643 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) { 01644 LoadShower(ii); 01645 totph+=ntpShower->ph.sigcor; 01646 } 01647 shwcount+=nshower; 01648 01649 //fiducial criteria on vtx for near detector 01650 is_fid = 1; 01651 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0; 01652 } 01653 01654 if (detector==2) is_fid_dp=PassFidNew(ntpEvent->index); 01655 else is_fid_dp=PassFid(ntpEvent->index); 01656 01657 if (detector==2) is_fid_ns=PassFidNewN(ntpEvent->index); 01658 else is_fid_ns=PassFid(ntpEvent->index); 01659 01660 //check for a corresponding mc event 01661 if(ntpHeader->GetVldContext().GetSimFlag()==4){ 01662 if(LoadTruthForRecoTH(j,mcevent)) { 01663 Float_t *vtx = TrueVtx(mcevent); 01664 Float_t *nu3mom = TrueNuMom(mcevent); 01665 Float_t *targ4vec = Target4Vector(mcevent); 01666 //true info for tree: 01667 is_mc = 1; 01668 nucleus = Nucleus(mcevent); 01669 flavour = Flavour(mcevent); 01670 initial_state = Initial_state(mcevent); 01671 had_fs = HadronicFinalState(mcevent); 01672 process = IResonance(mcevent); 01673 cc_nc = IAction(mcevent); 01674 true_enu = TrueNuEnergy(mcevent); 01675 true_pxnu = nu3mom[0]; 01676 true_pynu = nu3mom[1]; 01677 true_pznu = nu3mom[2]; 01678 true_emu = TrueMuEnergy(mcevent); 01679 true_eshw = TrueShwEnergy(mcevent); 01680 true_x = X(mcevent); 01681 true_y = Y(mcevent); 01682 true_q2 = Q2(mcevent); 01683 true_w2 = W2(mcevent); 01684 true_dircosz = TrueMuDCosZ(mcevent); 01685 true_dircosneu = TrueMuDCosNeu(mcevent); 01686 true_vtxx = vtx[0]; 01687 true_vtxy = vtx[1]; 01688 true_vtxz = vtx[2]; 01689 true_etgt = targ4vec[3]; 01690 true_pxtgt = targ4vec[0]; 01691 true_pytgt = targ4vec[1]; 01692 true_pztgt = targ4vec[2]; 01693 01694 /* 01695 if(gnumi.Status()) { 01696 if(isST) gnumi.GetParent(strecord,mcevent,*nuparent); 01697 else gnumi.GetParent(mcrecord,mcevent,*nuparent); 01698 } 01699 */ 01700 if (LoadTruth(mcevent)) { 01701 mcflux_tpx=ntpTruth->flux.tpx; 01702 mcflux_tpy=ntpTruth->flux.tpy; 01703 mcflux_tpz=ntpTruth->flux.tpz; 01704 mcflux_tptype=ntpTruth->flux.tptype; 01705 } 01706 01707 delete [] vtx; 01708 delete [] nu3mom; 01709 delete [] targ4vec; 01710 } 01711 } 01712 01713 if(PassBasicCuts(event)) {pass_basic=1;} // make this reflect pass/fail flag 01714 01715 //reco info for tree: 01716 reco_vtxx = ntpEvent->vtx.x; 01717 reco_vtxy = ntpEvent->vtx.y; 01718 reco_vtxz = ntpEvent->vtx.z; 01719 evtlength = ntpEvent->plane.end-ntpEvent->plane.beg+1; 01720 evtph = ntpEvent->ph.sigcor; 01721 01722 // different definition for neardet - accounts for 01723 // uninstrumented planes 01724 01725 AnnInputBlock anninput=AnnVar(event); 01726 01727 ann_aTotrk = anninput.aTotrk; 01728 ann_aTotstp = anninput.aTotstp; 01729 ann_aTotph = anninput.aTotph; 01730 ann_aTotlen = anninput.aTotlen; 01731 ann_aPhperpl = anninput.aPhperpl; 01732 ann_aPhperstp = anninput.aPhperstp; 01733 01734 01735 ann_aTrkpass = anninput.aTrkpass; 01736 ann_aTrkph = anninput.aTrkph; 01737 ann_aTrklen = anninput.aTrklen; 01738 ann_aTrkphperpl = anninput.aTrkphperpl; 01739 ann_aTrkphper = anninput.aTrkphper; 01740 ann_aTrkplu = anninput.aTrkplu; 01741 ann_aTrkplv = anninput.aTrkplv; 01742 ann_aTrkstp = anninput.aTrkstp; 01743 01744 ann_aShwph = anninput.aShwph; 01745 ann_aShwstp = anninput.aShwstp; 01746 ann_aShwdig = anninput.aShwdig; 01747 ann_aShwpl = anninput.aShwpl; 01748 ann_aShwphper = anninput.aShwphper; 01749 ann_aShwphperpl = anninput.aShwphperpl; 01750 ann_aShwphperdig = anninput.aShwphperdig ; 01751 ann_aShwphperstp = anninput.aShwphperstp; 01752 ann_aShwplu = anninput.aShwplu; 01753 ann_aShwplv = anninput.aShwplv; 01754 ann_aPh3 = anninput.aPh3; 01755 ann_aPh6 = anninput.aPh6; 01756 ann_aPhlast = anninput.aPhlast; 01757 ann_aPhcommon = anninput.aPhcommon; 01758 01759 // unused: Double_t trgtime=eventSummary->trigtime; 01760 // evttimemin = anninput.aTimemin-trgtime; 01761 // evttimemax = anninput.aTimemax-trgtime; 01762 01763 // cedar - use snarl time rather than trigger time 01764 evttimemin = anninput.aTimemin-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9; 01765 evttimemax = anninput.aTimemax-ntpHeader->GetVldContext().GetTimeStamp().GetNanoSec()/1e9; 01766 01767 01768 // ANN PID // 01769 01770 01771 Int_t target=1; 01772 if(ntpHeader->GetVldContext().GetSimFlag()!=4){ 01773 //if(tar<5000. && tar>3500.) target=1; 01774 //if(tar>5000. && tar<40000.) target=2; 01775 //if(tar>40000. && tar<100000.) target=3; // Used to work just fine for the pre shutdown data 01776 // after the shutdown this variable was not properly filled from ACNET data. in fact it is not 01777 // recoverable since it is corrupted in ACNET data. 01778 if(mcneartype==1) target=1; 01779 if(mcneartype==2) target=2; 01780 if(mcneartype==3) target=3; 01781 } 01782 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==1){ 01783 if(mcneartype==1) target=1; 01784 if(mcneartype==2) target=2; 01785 if(mcneartype==3) target=3; 01786 } 01787 if(ntpHeader->GetVldContext().GetSimFlag()==4 && detector==2){ 01788 target=0; // FOR THE MOMENT IN THE ABSENCE OF ANY FAR PME PHE MC 01789 } 01790 // SET FIDUCIAL AND TARGET DEFAULT IS DAVIDS FIDUCIAL AND NO OSCILLATION TRAINING 01791 01792 Int_t fid =1; 01793 Int_t prior = 1; 01794 01795 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ 01796 if(!PassFid(event)) annpid=-999; 01797 else { 01798 annpid=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01799 first_ann=false; 01800 } 01801 } 01802 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ 01803 if(!is_fid_ns){ 01804 annpid_f2p1=-999; 01805 annpid_f2p2=-999; 01806 annpid_f2p3=-999; 01807 } 01808 if(is_fid_ns){ 01809 fid=2; 01810 prior=1; 01811 annpid_f2p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01812 prior=2; 01813 annpid_f2p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01814 prior=3; 01815 annpid_f2p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01816 first_ann=false; 01817 } 01818 if(!is_fid_dp){ 01819 annpid_f1p1=-999; 01820 annpid_f1p2=-999; 01821 annpid_f1p3=-999; 01822 } 01823 if(is_fid_dp){ 01824 fid=1; 01825 prior=1; 01826 annpid_f1p1=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01827 prior=2; 01828 annpid_f1p2=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01829 prior=3; 01830 annpid_f1p3=GetAnnPid(anninput,detector,target,fid,prior,first_ann); 01831 first_ann=false; 01832 } 01833 } 01834 // MAK March 8, 2005 01835 // data or MC? 01836 // needed by energy corrections below 01837 bool isdata=true; 01838 if(ntpHeader->GetVldContext().GetSimFlag()==4) isdata=false; 01839 01840 //index will -1 if no track/shower present in event 01841 int track_index = -1; 01842 LoadLargestTrackFromEvent(j,track_index); 01843 int shower_index = -1; 01844 //LoadLargestShowerFromEvent(j,shower_index); // Change 01845 LoadShower_Jim(j,shower_index,detector); 01846 if(shower_index==-1) nshower=0; 01847 if(shower_index!=-1){ //check shower is present before using ntpShower 01848 01849 shwdigits = ntpShower->ndigit; 01850 shwstrips = ntpShower->nstrip; 01851 shwph = ntpShower->ph.sigcor; 01852 shwvtxu = ntpShower->vtx.u; 01853 shwvtxv = ntpShower->vtx.v; 01854 shwvtxx = ntpShower->vtx.x; 01855 shwvtxy = ntpShower->vtx.y; 01856 shwvtxz = ntpShower->vtx.z; 01857 shwncluster = ntpShower->ncluster; 01858 01859 // MAK March 8, 2005 01860 Int_t shwtype=0; 01861 const Int_t shw_correction_mode=1; // Niki's re-tuning 01862 shwtype=-1; 01863 reco_eshw = RecoShwEnergy(shower_index,shwtype, 01864 detector,isdata, 01865 shw_correction_mode); 01866 shwtype=0; 01867 reco_eshwcc = RecoShwEnergy(shower_index,shwtype, 01868 detector,isdata, 01869 shw_correction_mode); 01870 shwtype=2; 01871 reco_eshwnc = RecoShwEnergy(shower_index,shwtype, 01872 detector,isdata, 01873 shw_correction_mode); 01874 shwtype=1; 01875 reco_eshwccw = RecoShwEnergy(shower_index,shwtype, 01876 detector,isdata, 01877 shw_correction_mode); 01878 shwtype=3; 01879 reco_eshwncw = RecoShwEnergy(shower_index,shwtype, 01880 detector,isdata, 01881 shw_correction_mode); 01882 01883 reco_eshw_sqrt = RecoShwEnergySqrt(event); 01884 01885 } 01886 01887 01888 if(track_index!=-1){ //check track is present before using ntpTrack 01889 // MAK March 8, 2005 01890 reco_emu = RecoMuEnergy(0,track_index,isdata, 01891 ntpHeader->GetVldContext().GetDetector()); 01892 trklength = ntpTrack->plane.end-ntpTrack->plane.beg+1; 01893 trkmomrange = ntpTrack->momentum.range; 01894 trkqp = ntpTrack->momentum.qp; 01895 trkeqp = ntpTrack->momentum.eqp; 01896 trkendz = ntpTrack->end.z; 01897 trkendu = ntpTrack->end.u; 01898 trkendv = ntpTrack->end.v; 01899 trkendx = ntpTrack->end.x; 01900 trkendy = ntpTrack->end.y; 01901 01902 trkvtxz = ntpTrack->vtx.z; 01903 trkvtxu = ntpTrack->vtx.u; 01904 trkvtxv = ntpTrack->vtx.v; 01905 trkvtxx = ntpTrack->vtx.x; 01906 trkvtxy = ntpTrack->vtx.y; 01907 trkds = ntpTrack->ds; 01908 trkplbu = ntpTrack->plane.begu; 01909 trkplbv = ntpTrack->plane.begv; 01910 trkpleu = ntpTrack->plane.endu; 01911 trkplev = ntpTrack->plane.endv; 01912 trkple = ntpTrack->plane.end; 01913 trkplb = ntpTrack->plane.beg; 01914 trkfitpass = ntpTrack->fit.pass; 01915 trkdiffuv = abs(ntpTrack->plane.nu-ntpTrack->plane.nv); 01916 trkdigits = ntpTrack->ndigit; 01917 trkstrips = ntpTrack->nstrip; 01918 trkph = ntpTrack->ph.sigcor; 01919 trkndof = ntpTrack->fit.ndof; 01920 trkchi2 = ntpTrack->fit.chi2; 01921 01922 Float_t phtrack =ntpTrack->ph.sigcor; 01923 Float_t phevent =ntpEvent->ph.sigcor; 01924 01925 /* In case of Birch 01926 // apply 1.8% correction to pulse height variables 01927 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 01928 ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData ){ 01929 phtrack=phtrack*1.018; 01930 phevent=phevent*1.018; 01931 } 01932 */ 01933 01934 if (phevent>0) {trkphfrac=phtrack/phevent;} 01935 if(ntpTrack->plane.n>0) { 01936 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n; 01937 } 01938 01939 is_cont_trk = ntpTrack->contained; 01940 trkmombest =ntpTrack->momentum.best; 01941 } 01942 else { 01943 trklength = 0; 01944 trkmomrange = 0; 01945 trkqp = 0; 01946 trkeqp = 0; 01947 trkphfrac = 0; 01948 trkphplane = 0; 01949 } 01950 01951 reco_enu = reco_emu + reco_eshwcc; 01952 reco_qe_enu = RecoQENuEnergy(track_index); 01953 if (reco_enu>0) { 01954 reco_y = reco_eshwcc/reco_enu; 01955 } 01956 reco_qe_q2 = RecoQEQ2(track_index); 01957 reco_dircosz = RecoMuDCosZ(track_index); 01958 reco_dircosneu = RecoMuDCosNeu(track_index); 01959 is_cev = IsFidAll(track_index); 01960 01961 // only calculate likelihood PID for events that pass event 01962 // fiducial cuts and track quality cuts 01963 01964 if(is_fid_dp){ 01965 pid0 = PID(event,0); 01966 pid1 = PID(event,1); 01967 pid2 = PID(event,2); 01968 if(PassAnalysisCuts(event)) pass = 1; 01969 } 01970 else{ 01971 pid0 = -999.; 01972 pid1 = -999.; 01973 pid2 = -999.; 01974 pass = 0; 01975 } 01976 tree->Fill(); 01977 01978 } // end of EVENT LOOP 01979 } // end of SCANINFO loop 01980 } // end of NTUPLE ENTRIES loop 01981 01982 // delete nuparent; 01983 delete spinfo; 01984 01985 save.cd(); 01986 tree->Write(); 01987 save.Write(); 01988 save.Close(); 01989 }
| Double_t MadDpAnalysis::GetAnnPid | ( | AnnInputBlock | anninput, | |
| Int_t | det, | |||
| Int_t | tar, | |||
| Int_t | fid, | |||
| Int_t | prior, | |||
| Bool_t | first_ann | |||
| ) | [protected] |
Definition at line 392 of file MadDpAnalysis.cxx.
References Ann_weight_vector, AnnInputBlock::aPh3, AnnInputBlock::aPh6, AnnInputBlock::aPhlast, AnnInputBlock::aPhperpl, AnnInputBlock::aShwplu, AnnInputBlock::aTotlen, AnnInputBlock::aTotph, AnnInputBlock::aTotrk, AnnInputBlock::aTrkpass, AnnInputBlock::aTrkphperpl, AnnInputBlock::aTrkplu, and Sigmoid().
Referenced by CreatePAN().
00393 { 00394 00395 00396 // DET IS THE DETECTOR 1 = NEAR 2 = FAR 00397 // TAR IS THE ENERGY 1 = LE 2 = PME 3 = PHE 00398 // FID IS THE FIDUCIAL REGION FOR FAR , 2 = DAVIDS 1 = MINE (MORE STRICT) 00399 // PRIOR IS THE FAR TRAINED METHOD , 1 = NO OSCILLATIONS , 2 = DM2 0.0025 SINTHETA2 = 0.95 00400 // 3 = DM2 0.002 SINTHETA2 = 0.95 00401 00402 00403 // int inneuron = 13; 00404 int inneuron = 7; 00405 int hidneuron = 15; 00406 double var; 00407 00408 double out[15]; // input neurons 00409 double rin[15]; // first hidden layer 00410 00411 // double weight1[13][15];// weights first layer 00412 double weight1[7][15]; 00413 double constant1[15]; // constant first layer 00414 00415 double weighto[15];// weights second (output) layer 00416 double constanto[1]; // constant second (output) layer 00417 double prob; 00418 00419 // Initialize 00420 00421 for(Int_t i=0;i<hidneuron;i++) { 00422 out[i]=0.; 00423 rin[i]=0.; 00424 constant1[i]=0.; 00425 } 00426 00427 for(Int_t i=0;i<inneuron;i++) { 00428 for(Int_t j=0;j<hidneuron;j++){ 00429 weight1[i][j]=0.; 00430 } 00431 } 00432 00433 constanto[0]=0.; 00434 for(Int_t i=0;i<hidneuron;i++){ 00435 weighto[i]=0.; 00436 } 00437 Int_t all = inneuron*hidneuron+hidneuron+hidneuron+1; 00438 Int_t all1 = inneuron*hidneuron+hidneuron; 00439 Int_t n1 =0; 00440 Int_t nw1 =0; 00441 Int_t no =0; 00442 Int_t nwo =0; 00443 ifstream weightfile; 00444 00445 // INPUT VARIABLES 00446 if(first_ann){ 00447 00448 if(det==2){ 00449 if(prior==1){ 00450 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00451 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00452 } 00453 if(prior==2){ 00454 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00455 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00456 } 00457 if(prior==3){ 00458 if(fid==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00459 if(fid==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_farle_cedar_daikon7v2.dat"); 00460 } 00461 } 00462 00463 if(det==1){ 00464 if(tar==1) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearle_cedar_daikon7v2.dat"); 00465 if(tar==2) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearpme_cedar_daikon7v2.dat"); 00466 if(tar==3) weightfile.open("/afs/fnal.gov/files/data/minos/d80/niki/Ann_Files/weights_nearphe_cedar_daikon7v2.dat"); 00467 } 00468 for(Int_t k=0;k<all;k++){ 00469 weightfile >> var ; 00470 Ann_weight_vector[k]=var; 00471 } 00472 weightfile.close(); 00473 } 00474 if(det==2){ 00475 00476 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/2000.; 00477 out[1] = (anninput.aPhperpl/10000.); 00478 // out[1] = (anninput.aTotph/50000.); 00479 out[2] = (anninput.aTotlen/40); 00480 out[3] = (anninput.aPh3/(anninput.aTotph+1.)); 00481 out[4] = (anninput.aPh6/(anninput.aTotph+1.)); 00482 out[5] = (anninput.aPhlast/(anninput.aTotph+1.)); 00483 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.; 00484 /* 00485 out[0] = anninput.aTrkphperpl/2000.; 00486 out[1] = anninput.aTotrk*anninput.aTrkpass; 00487 out[2] = (anninput.aTotstp/100.); 00488 out[3] = (anninput.aTotph/50000.); 00489 out[4] = (anninput.aTotlen/40); 00490 out[5] = (anninput.aPhperpl/5000.); 00491 out[6] = (anninput.aPhperstp/1000.); 00492 out[10] = ((anninput.aTrkplv-anninput.aShwplv)+20.)/40.; 00493 out[11] = ((anninput.aTrkplu-anninput.aShwplu)+20.)/40.; 00494 */ 00495 } 00496 00497 00498 if(det==1){ 00499 00500 out[0] = (anninput.aTrkphperpl*anninput.aTotrk*anninput.aTrkpass)/5000.; 00501 out[1] = (anninput.aPhperpl/10000.); 00502 //out[1] = (anninput.aTotph/100000.); 00503 out[2] = (anninput.aTotlen/40); 00504 out[3] = (anninput.aPh3/(anninput.aTotph+1.)); 00505 out[4] = (anninput.aPh6/(anninput.aTotph+1.)); 00506 out[5] = (anninput.aPhlast/(anninput.aTotph+1.)); 00507 out[6] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.; 00508 /* 00509 out[0] = anninput.aTrkphperpl/5000.; 00510 out[1] = anninput.aTotrk*anninput.aTrkpass; 00511 out[2] = (anninput.aTotstp/100.); 00512 out[3] = (anninput.aTotph/100000.); 00513 out[4] = (anninput.aTotlen/40); 00514 out[5] = (anninput.aPhperpl/10000.); 00515 out[6] = (anninput.aPhperstp/2000.); 00516 out[10] = ((anninput.aTrkplv-anninput.aShwplv)+10.)/20.; 00517 out[11] = ((anninput.aTrkplu-anninput.aShwplu)+10.)/20.; 00518 */ 00519 } 00520 /* 00521 out[7] = (anninput.aPh3/(anninput.aTotph+1.)); 00522 out[8] = (anninput.aPh6/(anninput.aTotph+1.)); 00523 out[9] = (anninput.aPhlast/(anninput.aTotph+1.)); 00524 out[12] = (anninput.aShwphperdig/800.); 00525 */ 00526 // INPUT FILE 00527 00528 00529 for(Int_t k=0;k<all;k++){ 00530 if(k<all1){ 00531 if(k%(inneuron+1)==0){ 00532 nw1=0; 00533 constant1[n1]=Ann_weight_vector[k]; 00534 n1=n1+1; 00535 } 00536 else { 00537 weight1[nw1][n1-1]=Ann_weight_vector[k]; 00538 nw1=nw1+1; 00539 } 00540 } 00541 else if(k>=all1){ 00542 if((k-all1)%(hidneuron+1)==0){ 00543 nwo=0; 00544 constanto[no]=Ann_weight_vector[k]; 00545 no=no+1; 00546 } 00547 else { 00548 weighto[nwo]=Ann_weight_vector[k]; 00549 nwo=nwo+1; 00550 } 00551 } 00552 } 00553 00554 00555 // end of reading the weights 00556 00557 // first layer 00558 for(Int_t i=0;i<hidneuron;i++){ 00559 rin[i]=constant1[i]; 00560 for(Int_t j=0;j<inneuron;j++){ 00561 rin[i]=rin[i]+weight1[j][i]*out[j]; 00562 } 00563 } 00564 // output 00565 for(Int_t i=0;i<hidneuron;i++){ 00566 out[i]=Sigmoid(rin[i]); 00567 } 00568 00569 prob=constanto[0]; 00570 for(Int_t i=0;i<hidneuron;i++) prob=prob+weighto[i]*out[i]; 00571 00572 00573 if(anninput.aTotlen>50) prob=0.; 00574 00575 return prob; 00576 }
| bool MadDpAnalysis::InFidVol | ( | const Detector::Detector_t & | det, | |
| const Float_t & | x, | |||
| const Float_t & | y, | |||
| const Float_t & | z | |||
| ) | [static] |
Definition at line 2017 of file MadDpAnalysis.cxx.
References Detector::kFar, and Detector::kNear.
02019 { 02020 // Mike Kordosky: 11 August 2005 02021 // Defines a fiducial volume used to pass events which fill the PDFs 02022 // Transcribed from MakeMyFile() above 02023 02024 bool is_fid=false; 02025 02026 if(det==Detector::kFar){ 02027 02028 02029 is_fid = true; 02030 if(z<0.5 || z>29.4 || //ends 02031 (z<16.5&&z>14.5) || //between SMs 02032 sqrt((x*x) //radial cut 02033 +(y*y))>3.5) is_fid = false; 02034 } 02035 else if(det==Detector::kNear){ 02036 02037 //fiducial criteria on vtx for near detector 02038 is_fid = true; 02039 if(z<1 || z>5 || 02040 sqrt(((x-1.4885)*(x-1.4885)) + 02041 ((y-0.1397)*(y-0.1397)))>1) is_fid = false; 02042 } 02043 02044 return is_fid; 02045 02046 }
| void MadDpAnalysis::MakeMyFile | ( | std::string | , | |
| int | , | |||
| int | ||||
| ) |
Definition at line 657 of file MadDpAnalysis.cxx.
References NtpSRPlane::beg, NtpSRPlane::end, MadBase::eventSummary, MadQuantities::Flavour(), VldContext::GetDetector(), RecHeader::GetVldContext(), MadQuantities::IAction(), Detector::kFar, Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::LoadShower(), MadBase::LoadTHTrack(), MadBase::LoadTrack(), MadBase::LoadTruthForRecoTH(), NtpSRPlane::n, MadBase::Nentries, NtpSREventSummary::nevent, NtpSREvent::nshower, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpShower, MadBase::ntpTrack, NtpSREvent::ntrack, NtpSRShower::ph, NtpSRTrack::ph, NtpSREvent::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSRPulseHeight::sigcor, MadQuantities::TrueNuEnergy(), NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.
00657 { 00658 00659 00660 // pdf binning: 0 - old binning, max 50 bins, 1 - new finer binning 00661 00662 std::string savename = "DPHistos_" + tag + ".root"; 00663 TFile save(savename.c_str(),"RECREATE"); 00664 save.cd(); 00665 00666 00667 00668 00669 //#declare lots of histos: 00670 00671 Int_t evlbins=0; 00672 00673 if (pdfbinning==0) evlbins=60; 00674 else evlbins=300; 00675 00676 00677 TH1F *evlength_nc = new TH1F("Event length - nc", 00678 "Event length - nc", 00679 evlbins,0,600); 00680 evlength_nc->SetXTitle("Event length (planes)"); 00681 TH1F *evlength_cc = new TH1F("Event length - cc", 00682 "Event length - cc", 00683 evlbins,0,600); 00684 evlength_cc->SetXTitle("Event length (planes)"); 00685 00686 00687 00688 TH1F *phfrac_nc = new TH1F("Track ph frac - nc", 00689 "Track ph frac - nc", 00690 50,0,1); 00691 phfrac_nc->SetXTitle("pulse height fraction"); 00692 00693 00694 TH1F *phfrac_cc = new TH1F("Track ph frac - cc", 00695 "Track ph frac - cc", 00696 50,0,1); 00697 phfrac_cc->SetXTitle("pulse height fraction"); 00698 00699 00700 00701 Int_t phplanebins=0; 00702 Float_t phplanemax=0; 00703 00704 if (pdfbinning==0) { 00705 phplanebins=50; 00706 phplanemax=5000; 00707 } 00708 else { 00709 phplanebins=100; 00710 phplanemax=3000; 00711 } 00712 00713 TH1F *phplane_nc = new TH1F("ph per plane (nc)", 00714 "ph per plane (nc)", 00715 phplanebins,0,phplanemax); 00716 phplane_nc->SetXTitle("pulse height per plane"); 00717 TH1F *phplane_cc = new TH1F("ph per plane (cc)", 00718 "ph per plane (cc)", 00719 phplanebins,0,phplanemax); 00720 phplane_cc->SetXTitle("pulse height per plane"); 00721 00722 00723 evlength_nc->SetDirectory(&save); 00724 evlength_cc->SetDirectory(&save); 00725 00726 phfrac_nc->SetDirectory(&save); 00727 phfrac_cc->SetDirectory(&save); 00728 00729 phplane_nc->SetDirectory(&save); 00730 phplane_cc->SetDirectory(&save); 00731 00732 00733 00734 00735 00736 // gDirectory->ls(); 00737 for(int i=0;i<Nentries;i++){ 00738 00739 if(i%1000==0) std::cout << float(i)*100./float(Nentries) 00740 << "% done" << std::endl; 00741 if(!this->GetEntry(i)) continue; 00742 00743 Int_t nevent = eventSummary->nevent; 00744 if(nevent==0) continue; 00745 00746 Int_t trkcount=0; 00747 Int_t shwcount=0; 00748 00749 Int_t evtmin=0; 00750 Int_t evtmax=nevent; 00751 00752 if (nevent>1 && ntpHeader->GetVldContext().GetDetector()==Detector::kFar) { 00753 00754 Int_t maxevent=0; 00755 Float_t maxph=0; 00756 for(int j=0;j<nevent;j++){ 00757 if(!LoadEvent(j)) continue; 00758 Float_t evtpheight=ntpEvent->ph.sigcor; 00759 if (evtpheight>maxph) { 00760 maxph=evtpheight; 00761 maxevent=j; 00762 } 00763 } 00764 evtmin=maxevent; 00765 evtmax=maxevent+1; 00766 } 00767 00768 // EVENT LOOP - fill tree once for each reconstructed event 00769 00770 for(int j=evtmin;j<evtmax;j++){ 00771 00772 00773 if(!LoadEvent(j)) continue; //no event found 00774 00775 Int_t ntrack = 0; 00776 ntrack=ntpEvent->ntrack; 00777 Int_t nshower = 0; 00778 nshower=ntpEvent->nshower; 00779 00780 Float_t totph=0; 00781 00782 00783 for (Int_t ii=trkcount;ii<ntrack+trkcount;ii++) { 00784 LoadTrack(ii); 00785 00786 totph+=ntpTrack->ph.sigcor; 00787 LoadTHTrack(ii); 00788 00789 } 00790 00791 trkcount+=ntrack; 00792 00793 00794 for (Int_t ii=shwcount;ii<nshower+shwcount;ii++) { 00795 LoadShower(ii); 00796 00797 totph+=ntpShower->ph.sigcor; 00798 } 00799 00800 shwcount+=nshower; 00801 00802 00803 00804 Int_t cc_nc = 0; 00805 Int_t flavour = 0; 00806 Int_t detector = -1; 00807 Int_t is_fid = 0; 00808 Double_t neu_e = -1; 00809 Double_t weight = -1; 00810 00811 //get detector type: 00812 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar){ 00813 detector = 2; 00814 //fiducial criteria on vtx for far detector 00815 if (evlength_cc->GetNbinsX()==evlbins) { 00816 if (pdfbinning==0) { 00817 evlength_cc->SetBins(50,0,500); 00818 evlength_nc->SetBins(50,0,500); 00819 } 00820 else { 00821 evlength_cc->SetBins(250,0,500); 00822 evlength_nc->SetBins(250,0,500); 00823 } 00824 } 00825 is_fid = 1; 00826 if(ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>29.4 || //ends 00827 (ntpEvent->vtx.z<16.5&&ntpEvent->vtx.z>14.5) || //between SMs 00828 sqrt((ntpEvent->vtx.x*ntpEvent->vtx.x) //radial cut 00829 +(ntpEvent->vtx.y*ntpEvent->vtx.y))>3.5) is_fid = 0; 00830 } 00831 else if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ 00832 detector = 1; 00833 00834 if (pdfbinning==1 && evlength_cc->GetBinLowEdge(250)>200) { 00835 evlength_cc->SetBins(150,0,300); 00836 evlength_nc->SetBins(150,0,300); 00837 } 00838 if (pdfbinning==0 && evlength_cc->GetBinLowEdge(50)>200) { 00839 evlength_cc->SetBins(60,0,300); 00840 evlength_nc->SetBins(60,0,300); 00841 } 00842 00843 00844 //fiducial criteria on vtx for near detector 00845 is_fid = 1; 00846 if(ntpEvent->vtx.z<1 || ntpEvent->vtx.z>5 || 00847 sqrt(((ntpEvent->vtx.x-1.4885)*(ntpEvent->vtx.x-1.4885)) + 00848 ((ntpEvent->vtx.y-0.1397)*(ntpEvent->vtx.y-0.1397)))>1) is_fid = 0; 00849 } 00850 00851 //check for a corresponding mc event 00852 00853 Int_t mcevent=-1; 00854 00855 if(LoadTruthForRecoTH(j,mcevent)) { 00856 //true info for tree: 00857 00858 cc_nc = IAction(mcevent); 00859 flavour = Flavour(mcevent); 00860 neu_e = TrueNuEnergy(mcevent); 00861 00862 // tarpos 1 = LE 2 = ME 3 = HE 00863 if(tarpos==2) weight=1; 00864 if(tarpos==1) weight=1; 00865 if(tarpos==3) weight=1; 00866 } 00867 00868 int track_index = -1; 00869 LoadLargestTrackFromEvent(j,track_index); 00870 00871 if (track_index!=-1) { 00872 00873 // if (is_fid && trkpass) { 00874 if (is_fid) { 00875 00876 Float_t evlength=0; 00877 Float_t phfrac=0; 00878 Float_t phplane=0; 00879 00880 evlength=ntpEvent->plane.n; 00881 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ 00882 evlength=ntpEvent->plane.end-ntpEvent->plane.beg+1; 00883 } 00884 00885 Float_t phtrack=ntpTrack->ph.sigcor; 00886 Float_t phevent=ntpEvent->ph.sigcor; 00887 if (phevent>0) {phfrac=phtrack/phevent;} 00888 00889 Float_t trkplane=ntpTrack->plane.n; 00890 if (trkplane>0) {phplane=phtrack/trkplane;} 00891 00892 if(flavour==2 && cc_nc==1) { 00893 evlength_cc->Fill(evlength,weight); 00894 phfrac_cc->Fill(phfrac,weight); 00895 phplane_cc->Fill(phplane,weight); 00896 // cout<<"I've filled a histogram, where it is, I don't know"<<endl; 00897 // cout<<phplane_cc->GetEntries()<<endl; 00898 // gDirectory->ls(); 00899 } 00900 00901 else if(cc_nc==0) { 00902 evlength_nc->Fill(evlength,weight); 00903 phfrac_nc->Fill(phfrac,weight); 00904 phplane_nc->Fill(phplane,weight); 00905 } 00906 // std::cout<<"evlenght "<<evlength<<" weight "<<weight<<std::endl; 00907 00908 } 00909 } 00910 } 00911 } 00912 //#close file 00913 // std::cout<<"writing hist "<<endl; 00914 // save.ls(); 00915 save.Write(); 00916 save.Close(); 00917 }
| Bool_t MadDpAnalysis::PassAnalysisCuts | ( | Int_t | event = 0 |
) | [protected, virtual] |
Reimplemented from MadAnalysis.
Definition at line 88 of file MadDpAnalysis.cxx.
References VldContext::GetDetector(), RecHeader::GetVldContext(), MadBase::LoadEvent(), MadBase::ntpHeader, and PID().
Referenced by CreatePAN().
00089 { 00090 if(!LoadEvent(event)) return false; 00091 // if(PassBasicCuts(event)) return true; 00092 Float_t pidcut=-0.4; 00093 if(ntpHeader->GetVldContext().GetDetector()==1) {pidcut=-0.1;} 00094 00095 if(!(PID(event,0)>pidcut)) return false; 00096 return true; 00097 }
| Bool_t MadDpAnalysis::PassBasicCuts | ( | Int_t | event = 0 |
) | [protected, virtual] |
Reimplemented from MadAnalysis.
Definition at line 152 of file MadDpAnalysis.cxx.
References MadBase::LoadEvent(), MadBase::ntpEvent, and NtpSREvent::ntrack.
Referenced by CreatePAN().
00153 { 00154 if(!LoadEvent(event)) return false; 00155 if(ntpEvent->ntrack==0) return false; 00156 // Int_t track = -1; 00157 // LoadLargestTrackFromEvent(event,track); 00158 // if(!PassTrackCuts(track)) return false; 00159 // if(!IsFidVtxEvt(ntpEvent->index)) return false; 00160 // if(!IsFidVtx(track)) return false; 00161 return true; 00162 }
| Bool_t MadDpAnalysis::PassFid | ( | Int_t | event = 0 |
) | [protected] |
Definition at line 163 of file MadDpAnalysis.cxx.
References NtpSREvent::index, MadQuantities::IsFidVtxEvt(), MadBase::LoadEvent(), and MadBase::ntpEvent.
Referenced by CreatePAN().
00164 { 00165 if(!LoadEvent(event)) return false; 00166 if(!IsFidVtxEvt(ntpEvent->index)) return false; 00167 return true; 00168 }
| Bool_t MadDpAnalysis::PassFidNew | ( | Int_t | event = 0 |
) | [protected] |
Definition at line 98 of file MadDpAnalysis.cxx.
References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSRTrack::vtx, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.
Referenced by CreatePAN().
00099 { 00100 00101 if(!LoadEvent(event)) return false; 00102 00103 Int_t track = -1; 00104 LoadLargestTrackFromEvent(event,track); 00105 00106 // for events with no track, use event vertex 00107 00108 if (track==-1 && (ntpEvent->vtx.z<0.5 || ntpEvent->vtx.z>28.0 || //ends 00109 (ntpEvent->vtx.z<16.2 && ntpEvent->vtx.z>14.3) || //between SMs 00110 (pow(ntpEvent->vtx.x,2)+ //radial cut 00111 pow(ntpEvent->vtx.y,2))>14)) return false; 00112 00113 // otherwise, use track vertex 00114 00115 if (track!=-1 && (ntpTrack->vtx.z<0.5 || ntpTrack->vtx.z>28.0 || //ends 00116 (ntpTrack->vtx.z<16.2 && ntpTrack->vtx.z>14.3) || //between SMs 00117 (pow(ntpTrack->vtx.x,2)+ //radial cut 00118 pow(ntpTrack->vtx.y,2))>14)) return false; 00119 00120 return true; 00121 00122 }
| Bool_t MadDpAnalysis::PassFidNewN | ( | Int_t | event = 0 |
) | [protected] |
Definition at line 126 of file MadDpAnalysis.cxx.
References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadBase::ntpEvent, MadBase::ntpTrack, NtpSRTrack::vtx, NtpSREvent::vtx, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.
Referenced by CreatePAN().
00127 { 00128 00129 if(!LoadEvent(event)) return false; 00130 00131 Int_t track = -1; 00132 LoadLargestTrackFromEvent(event,track); 00133 00134 // for events with no track, use event vertex 00135 00136 if (track==-1 && (ntpEvent->vtx.z<1.0 || ntpEvent->vtx.z>27.0 || //ends 00137 (ntpEvent->vtx.z<17.0 && ntpEvent->vtx.z>14.0) || //between SMs 00138 (pow(ntpEvent->vtx.x,2)+ //radial cut 00139 pow(ntpEvent->vtx.y,2))>(3.5*3.5))) return false; 00140 00141 // otherwise, use track vertex 00142 00143 if (track!=-1 && (ntpTrack->vtx.z<1.0 || ntpTrack->vtx.z>27.0 || //ends 00144 (ntpTrack->vtx.z<17. && ntpTrack->vtx.z>14.0) || //between SMs 00145 (pow(ntpTrack->vtx.x,2)+ //radial cut 00146 pow(ntpTrack->vtx.y,2))>(3.5*3.5))) return false; 00147 00148 return true; 00149 00150 }
| Bool_t MadDpAnalysis::PassTruthSignal | ( | Int_t | event = 0 |
) | [protected, virtual] |
Reimplemented from MadAnalysis.
Definition at line 80 of file MadDpAnalysis.cxx.
References NtpMCTruth::iaction, NtpMCTruth::inu, MadBase::LoadTruth(), and MadBase::ntpTruth.
00081 { 00082 if(!LoadTruth(mcevent)) return false; 00083 if(abs(ntpTruth->inu)==14&&ntpTruth->iaction==1) return true; 00084 return false; 00085 }
| Float_t MadDpAnalysis::PID | ( | Int_t | event = 0, |
|
| Int_t | method = 0 | |||
| ) | [protected, virtual] |
Reimplemented from MadAnalysis.
Definition at line 579 of file MadDpAnalysis.cxx.
References NtpSRPlane::beg, NtpSRPlane::end, MadBase::eventSummary, fLikeliHist, VldContext::GetDetector(), RecHeader::GetVldContext(), Detector::kNear, MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), NtpSRPlane::n, MadBase::ntpEvent, MadBase::ntpHeader, MadBase::ntpTrack, NtpSREvent::ph, NtpSRTrack::ph, NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, and NtpSRPulseHeight::sigcor.
Referenced by CreatePAN(), and PassAnalysisCuts().
00580 { 00581 00582 00583 if(!LoadEvent(event)) return -10; 00584 Int_t track = -1; 00585 if(!LoadLargestTrackFromEvent(event,track)) return -10; 00586 if(method!=0) return -10; 00587 00588 Float_t ProbNC = 1.; 00589 Float_t ProbCC = 1.; 00590 Float_t PidCC = 0.; 00591 00592 00593 Float_t var1=eventSummary->plane.n; 00594 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){ 00595 var1=ntpEvent->plane.end-ntpEvent->plane.beg+1; 00596 } 00597 Float_t var2=0; 00598 Float_t var3=0; 00599 00600 00601 Float_t phtrack=ntpTrack->ph.sigcor; 00602 Float_t phevent=ntpEvent->ph.sigcor; 00603 00604 /* In case of Birch 00605 // apply 1.8% correction to pulse height variables 00606 if(ntpHeader->GetVldContext().GetDetector()==Detector::kFar && 00607 ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kData ){ 00608 phtrack=phtrack*1.018; 00609 phevent=phevent*1.018; 00610 } 00611 */ 00612 00613 if (phevent>0) {var2=phtrack/phevent;} 00614 00615 if (ntpTrack->plane.n!=0) var3 = float(phtrack) 00616 /float(ntpTrack->plane.n); 00617 00618 00619 Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1)); 00620 Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2)); 00621 Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3)); 00622 00623 if (prob1==0) {prob1=0.0001;} 00624 if (prob2==0) {prob2=0.0001;} 00625 if (prob3==0) {prob3=0.0001;} 00626 00627 ProbCC=prob1*prob2*prob3; 00628 00629 prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1)); 00630 prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2)); 00631 prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3)); 00632 00633 if (prob1==0) {prob1=0.0001;} 00634 if (prob2==0) {prob2=0.0001;} 00635 if (prob3==0) {prob3=0.0001;} 00636 00637 ProbNC=prob1*prob2*prob3; 00638 00639 PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC))); 00640 00641 return PidCC; 00642 }
| void MadDpAnalysis::ReadPIDFile | ( | std::string | ) |
Definition at line 1993 of file MadDpAnalysis.cxx.
References fLikeliFile, and fLikeliHist.
01993 { 01994 01995 fLikeliFile = new TFile(tag.c_str(),"READ"); 01996 01997 if(fLikeliFile->IsOpen() && !fLikeliFile->IsZombie()) { //if file is found 01998 01999 fLikeliFile->cd(); 02000 02001 fLikeliHist[0] = (TH1F*) fLikeliFile->Get("Event length - cc"); 02002 fLikeliHist[1] = (TH1F*) fLikeliFile->Get("Event length - nc"); 02003 fLikeliHist[2] = (TH1F*) fLikeliFile->Get("Track ph frac - cc"); 02004 fLikeliHist[3] = (TH1F*) fLikeliFile->Get("Track ph frac - nc"); 02005 fLikeliHist[4] = (TH1F*) fLikeliFile->Get("ph per plane (cc)"); 02006 fLikeliHist[5] = (TH1F*) fLikeliFile->Get("ph per plane (nc)"); 02007 02008 for(int i=0;i<6;i++){ 02009 float integ = fLikeliHist[i]->Integral(); 02010 fLikeliHist[i]->Scale(1./integ); 02011 } 02012 } 02013 else fLikeliFile = NULL; 02014 }
| Float_t MadDpAnalysis::RecoAnalysisEnergy | ( | Int_t | event = 0 |
) | [protected, virtual] |
Reimplemented from MadAnalysis.
Definition at line 645 of file MadDpAnalysis.cxx.
References MadBase::LoadEvent(), MadBase::LoadLargestTrackFromEvent(), MadQuantities::RecoMuEnergy(), and MadQuantities::RecoShwEnergy().
00645 { 00646 00647 if(!LoadEvent(event)) return 0; 00648 int largestTrack = -1; 00649 LoadLargestTrackFromEvent(event,largestTrack); 00650 float emu = RecoMuEnergy(0,largestTrack); 00651 float ehad = RecoShwEnergy(event); 00652 float ereco = emu + ehad; 00653 return ereco; 00654 00655 }
| Double_t MadDpAnalysis::Sigmoid | ( | Double_t | x | ) | [protected] |
Definition at line 376 of file MadDpAnalysis.cxx.
Referenced by GetAnnPid().
00377 { 00378 double sig; 00379 if(x>37.){ 00380 sig = 1.; 00381 } 00382 else if(x<-37.){ 00383 sig = 0.; 00384 } 00385 else { 00386 sig = 1./(1.+exp(-x)); 00387 } 00388 return sig; 00389 }
| Float_t MadDpAnalysis::SingleTimeFrame | ( | Int_t | snarlentry, | |
| Int_t | nentries | |||
| ) | [protected] |
Definition at line 919 of file MadDpAnalysis.cxx.
References MadBase::GetEntry(), RecPhysicsHeader::GetTimeFrame(), and MadBase::ntpHeader.
00919 { 00920 00921 Int_t tf=0,tfbef=0,tfaft=0; 00922 00923 00924 if(this->GetEntry(snarlentry)) tf = ntpHeader->GetTimeFrame(); 00925 if(snarlentry<(nentries-1) && this->GetEntry(snarlentry+1)) tfaft = ntpHeader->GetTimeFrame(); 00926 if(snarlentry>0 && this->GetEntry(snarlentry-1)) tfbef = ntpHeader->GetTimeFrame(); 00927 00928 Float_t singletimeframe=0; 00929 if(snarlentry<(nentries-1) && snarlentry>0 ){ 00930 if(tf!=tfaft && tf!=tfbef) singletimeframe=1; 00931 } 00932 if(snarlentry==(nentries-1)) { 00933 if(tf!=tfbef) singletimeframe=1; 00934 } 00935 if(snarlentry==0) { 00936 if(tf!=tfaft) singletimeframe=1; 00937 } 00938 00939 this->GetEntry(snarlentry); 00940 return singletimeframe; 00941 00942 }
Double_t MadDpAnalysis::Ann_weight_vector[226] [private] |
TFile* MadDpAnalysis::fLikeliFile [protected] |
Definition at line 28 of file MadDpAnalysis.h.
Referenced by MadDpAnalysis(), ReadPIDFile(), and ~MadDpAnalysis().
TH1F* MadDpAnalysis::fLikeliHist[6] [protected] |
Definition at line 29 of file MadDpAnalysis.h.
Referenced by MadDpAnalysis(), PID(), and ReadPIDFile().
1.4.7