00001
00002
00003
00004
00005
00007
00008 #include <set>
00009 #include <cmath>
00010
00011 #include "TH2F.h"
00012 #include "TH3F.h"
00013 #include "TNtuple.h"
00014 #include "TProfile.h"
00015 #include "TSpline.h"
00016
00017 #include "BeamDataUtil/BeamMonSpill.h"
00018 #include "Conventions/BeamType.h"
00019 #include "BeamDataUtil/BMSpillAna.h"
00020 #include "DataUtil/infid_sr_interface.h"
00021 #include "JobControl/JobCEnv.h"
00022 #include "Mad/MadDpID.h"
00023 #include "Mad/MadAbID.h"
00024 #include "MessageService/MsgFormat.h"
00025 #include "MessageService/MsgService.h"
00026 #include "BeamDataNtuple/NtpBDLiteRecord.h"
00027 #include "CandFitTrackSA/Ntp/NtpFitSARecord.h"
00028 #include "CandNtupleSR/NtpSREvent.h"
00029 #include "CandNtupleSR/NtpSRTrack.h"
00030 #include "CandNtupleSR/NtpSRShower.h"
00031 #include "CandNtupleSR/NtpSRStrip.h"
00032 #include "StandardNtuple/NtpStRecord.h"
00033 #include "MCNtuple/NtpMCTruth.h"
00034 #include "TruthHelperNtuple/NtpTHEvent.h"
00035 #include "TruthHelperNtuple/NtpTHShower.h"
00036 #include "TruthHelperNtuple/NtpTHTrack.h"
00037 #include "Conventions/ReleaseType.h"
00038 #include "MCReweight/SKZPWeightCalculator.h"
00039 #include "UgliGeometry/UgliGeomHandle.h"
00040
00041 #include "NtupleUtils/LISieve.h"
00042 #include "NtupleUtils/NuCounter.h"
00043 #include "NtupleUtils/NuCuts.h"
00044 #include "NtupleUtils/NuEvent.h"
00045 #include "NtupleUtils/NuExtraction.h"
00046 #include "NtupleUtils/NuGeneral.h"
00047 #include "NtupleUtils/NuHistos.h"
00048 #include "NtupleUtils/NuInputEvents.h"
00049 #include "NtupleUtils/NuReco.h"
00050 #include "NtupleUtils/NuMCEvent.h"
00051 #include "NtupleUtils/NuZBeamReweight.h"
00052
00053 #include "NtupleUtils/NuLibrary.h"
00054
00055 #include "NuMuBar/NuAnalysis.h"
00056 #include "NuMuBar/NuBeam.h"
00057 #include "NuMuBar/NuConfig.h"
00058 #include "NuMuBar/NuMadAnalysis.h"
00059 #include "NuMuBar/NuOutputWriter.h"
00060 #include "NuMuBar/NuPlots.h"
00061 #include "NuMuBar/NuPIDInterface.h"
00062 #include "NuMuBar/NuTime.h"
00063
00064 using std::endl;
00065 using std::cout;
00066 using std::map;
00067 using std::vector;
00068
00069 CVSID("$Id: NuAnalysis.cxx,v 1.39 2008/03/12 15:47:55 hartnell Exp $");
00070
00071
00072
00073 NuAnalysis::NuAnalysis()
00074 {
00075 MSG("NuAnalysis",Msg::kDebug)
00076 <<"Running NuAnalysis Constructor..."<<endl;
00077
00078
00079 MSG("NuAnalysis",Msg::kDebug)
00080 <<"Finished NuAnalysis Constructor"<<endl;
00081 }
00082
00083
00084
00085 NuAnalysis::~NuAnalysis()
00086 {
00087 MSG("NuAnalysis",Msg::kDebug)
00088 <<"Running NuAnalysis Destructor..."<<endl;
00089
00090
00091 MSG("NuAnalysis",Msg::kDebug)
00092 <<"Finished NuAnalysis Destructor"<<endl;
00093 }
00094
00095
00096
00097 void NuAnalysis::DemoInfidSRInterface(const NtpStRecord& ntp,
00098 const NtpSREvent& evt,
00099 const NuEvent& nu)
00100 {
00101 MAXMSG("NuAnalysis",Msg::kInfo,3)
00102 <<"Running DemoInfidSRInterface..."<<endl;
00103
00104
00105
00106
00107
00108
00109 const NuLibrary& lib=NuLibrary::Instance();
00110
00111
00112 Int_t bestTrack=lib.reco.GetBestTrack(nu);
00113 if (bestTrack<1) return;
00114 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
00115 (ntp,evt,bestTrack-1);
00116 const NtpSRTrack& trk=*ptrk;
00117
00118
00119 Int_t bestShower=lib.reco.GetBestShower(nu);
00120 if (bestShower<1) return;
00121 const NtpSRShower* pshw=lib.reco.GetShowerWithIndexX
00122 (ntp,evt,bestShower-1);
00123 const NtpSRShower& shw=*pshw;
00124
00128
00129
00130
00131 static Bool_t infid_initialised=false;
00132 if (!infid_initialised) {
00133 infid_initialised=true;
00134 choose_infid_set("cc2008");
00135
00136
00137 for (Float_t z=0;z<30;z+=0.005) {
00138
00139
00140 Float_t x=1.5*Munits::m;
00141 Float_t y=0.6*Munits::m;
00142
00143
00144 const RecCandHeader& rec=ntp.GetHeader();
00145 Bool_t inFid=infid(rec.GetVldContext().GetDetector(),
00146 rec.GetVldContext().GetSimFlag(),
00147 x,y,z);
00148 MSG("NuAnalysis",Msg::kInfo)
00149 <<"z="<<z<<", inFid = "<<inFid<<endl;
00150 }
00151 }
00152
00153
00154
00155 Bool_t evtIsInFid=infid(ntp,evt);
00156 Bool_t trkIsInFid=infid(ntp,trk);
00157 Bool_t shwIsInFid=infid(ntp,shw);
00158
00159 MAXMSG("NuAnalysis",Msg::kInfo,100)
00160 <<"Evt/Trk/Shw IsInFid = "
00161 <<evtIsInFid<<"/"<<trkIsInFid<<"/"<<shwIsInFid<<endl
00162 <<"evt(x,y,z) = ("<<evt.vtx.x<<","<<evt.vtx.y<<","<<evt.vtx.z<<")"
00163 <<endl
00164 <<"trk(x,y,z) = ("<<trk.vtx.x<<","<<trk.vtx.y<<","<<trk.vtx.z<<")"
00165 <<endl
00166 <<"shw(x,y,z) = ("<<shw.vtx.x<<","<<shw.vtx.y<<","<<shw.vtx.z<<")"
00167 <<endl;
00168
00169
00170
00171 if (evtIsInFid!=trkIsInFid && evt.vtx.plane==trk.vtx.plane-1) {
00172 MAXMSG("NuAnalysis",Msg::kInfo,200)
00173 <<"*******************************************************"<<endl
00174 <<"*******************************************************"<<endl
00175 <<"*******************************************************"<<endl
00176 <<"Evt/Trk/Shw IsInFid = "
00177 <<evtIsInFid<<"/"<<trkIsInFid<<"/"<<shwIsInFid<<endl
00178 <<"evt(x,y,z) = ("<<evt.vtx.x<<","<<evt.vtx.y<<","<<evt.vtx.z<<")"
00179 <<endl
00180 <<"trk(x,y,z) = ("<<trk.vtx.x<<","<<trk.vtx.y<<","<<trk.vtx.z<<")"
00181 <<endl
00182 <<"shw(x,y,z) = ("<<shw.vtx.x<<","<<shw.vtx.y<<","<<shw.vtx.z<<")"
00183 <<endl
00184 <<"*******************************************************"<<endl
00185 <<"*******************************************************"<<endl
00186 <<"*******************************************************"<<endl
00187 <<endl;
00188 }
00189 }
00190
00191
00192
00193 void NuAnalysis::BasicPlots()
00194 {
00195 MSG("NuAnalysis",Msg::kInfo)
00196 <<" ** Running BasicPlots method... **"<<endl;
00197
00198
00199 fOutFile=this->OpenFile(100,"NuBasicPlots");
00200
00201 TH1F* hNuInt=new TH1F("hNuInt","hNuInt",10000,-100,1000);
00202 hNuInt->SetTitle("Number of neutrino interactions per spill");
00203 hNuInt->GetXaxis()->SetTitle("# interactions");
00204 hNuInt->GetXaxis()->CenterTitle();
00205 hNuInt->GetYaxis()->SetTitle("");
00206 hNuInt->GetYaxis()->CenterTitle();
00207 hNuInt->SetFillColor(0);
00208 hNuInt->SetLineColor(1);
00209
00210
00211 TH1F* hNuMuBarEn=new TH1F("hNuMuBarEn","hNuMuBarEn",4*352,-32,320);
00212 hNuMuBarEn->GetXaxis()->SetTitle("Energy (GeV)");
00213 hNuMuBarEn->GetXaxis()->CenterTitle();
00214 hNuMuBarEn->GetYaxis()->SetTitle("");
00215 hNuMuBarEn->GetYaxis()->CenterTitle();
00216 hNuMuBarEn->SetFillColor(0);
00217 hNuMuBarEn->SetLineColor(2);
00218
00219
00220 TH1F* hNuMuEn=new TH1F("hNuMuEn","hNuMuEn",4*352,-32,320);
00221 hNuMuEn->GetXaxis()->SetTitle("Energy (GeV)");
00222 hNuMuEn->GetXaxis()->CenterTitle();
00223 hNuMuEn->GetYaxis()->SetTitle("");
00224 hNuMuEn->GetYaxis()->CenterTitle();
00225 hNuMuEn->SetFillColor(0);
00226 hNuMuEn->SetLineColor(1);
00227 hNuMuEn->SetLineWidth(2);
00228
00229
00230 TH1F* hNueBarEn=new TH1F("hNueBarEn","hNueBarEn",4*352,-32,320);
00231 hNueBarEn->GetXaxis()->SetTitle("Energy (GeV)");
00232 hNueBarEn->GetXaxis()->CenterTitle();
00233 hNueBarEn->GetYaxis()->SetTitle("");
00234 hNueBarEn->GetYaxis()->CenterTitle();
00235 hNueBarEn->SetFillColor(0);
00236 hNueBarEn->SetLineColor(2);
00237 hNueBarEn->SetLineWidth(2);
00238 hNueBarEn->SetLineStyle(2);
00239
00240
00241 TH1F* hNueEn=new TH1F("hNueEn","hNueEn",4*352,-32,320);
00242 hNueEn->GetXaxis()->SetTitle("Energy (GeV)");
00243 hNueEn->GetXaxis()->CenterTitle();
00244 hNueEn->GetYaxis()->SetTitle("");
00245 hNueEn->GetYaxis()->CenterTitle();
00246 hNueEn->SetFillColor(0);
00247 hNueEn->SetLineColor(1);
00248 hNueEn->SetLineWidth(1);
00249 hNueEn->SetLineStyle(2);
00250
00251
00255
00256 this->InitialiseLoopVariables();
00257
00258
00259 for(Int_t entry=0;entry<this->GetEntries();entry++){
00260
00261 this->SetLoopVariables(entry);
00262
00263
00264 const NtpStRecord& ntp=(*this->GetNtpStRecord());
00265
00266 TClonesArray& mcTca=*ntp.mc;
00267 Int_t numInt=mcTca.GetEntries();
00268 MAXMSG("NuAnalysis",Msg::kInfo,20)
00269 <<"Number of entries in NtpMCTruth="
00270 <<mcTca.GetEntries()<<endl;
00271
00272
00273 hNuInt->Fill(mcTca.GetEntries());
00274
00275 for (Int_t i=0;i<numInt;i++){
00276 const NtpMCTruth& mc=*(dynamic_cast<NtpMCTruth*>(mcTca[i]));
00277
00278 string iaction="NC";
00279 if (mc.iaction==1) iaction="CC";
00280 string iresonance="QE ";
00281 if (mc.iresonance==1002) iresonance="RES";
00282 else if (mc.iresonance==1003) iresonance="DIS";
00283 else if (mc.iresonance==1004) iresonance=
00284 "Coherent pion production";
00285 string inu="??";
00286 if (mc.inu==14) inu="NuMu ";
00287 else if (mc.inu==-14) inu="NuMuBar";
00288 else if (mc.inu==12) inu="NuE ";
00289 else if (mc.inu==-12) inu="NuEBar ";
00290
00291 MAXMSG("NuAnalysis",Msg::kInfo,100)
00292 <<"Interaction #"<<i<<" is "<<inu
00293 <<" "<<iaction
00294 <<" "<<iresonance<<endl;
00295
00296 if (mc.inu==-14) hNuMuBarEn->Fill(mc.p4neu[3]);
00297 else if (mc.inu==14) hNuMuEn->Fill(mc.p4neu[3]);
00298 else if (mc.inu==-12) hNueBarEn->Fill(mc.p4neu[3]);
00299 else if (mc.inu==12) hNueEn->Fill(mc.p4neu[3]);
00300 else {
00301 MAXMSG("NuAnalysis",Msg::kWarning,100)
00302 <<"particle not defined for filling histo="<<mc.inu<<endl;
00303 }
00304 }
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 }
00323
00327
00328 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
00329
00330 MSG("NuAnalysis",Msg::kInfo)
00331 <<" ** Finished BasicPlots method **"<<endl;
00332 }
00333
00334
00335
00336 void NuAnalysis::EfficienciesOld()
00337 {
00338 MSG("NuAnalysis",Msg::kInfo)
00339 <<" ** Running Efficiencies method... **"<<endl;
00340
00341
00342 fOutFile=this->OpenFile(100,"NuEfficiencies");
00343
00344 TProfile* pNuMuBarQEffVsEn=new TProfile("pNuMuBarQEffVsEn","pNuMuBarQEffVsEn",50,0,50);
00345 pNuMuBarQEffVsEn->SetTitle("Charge Sign Efficiency vs. True Neutrino Energy");
00346 pNuMuBarQEffVsEn->GetXaxis()->SetTitle("Energy (GeV)");
00347 pNuMuBarQEffVsEn->GetXaxis()->CenterTitle();
00348 pNuMuBarQEffVsEn->GetYaxis()->SetTitle("Efficiency");
00349 pNuMuBarQEffVsEn->GetYaxis()->CenterTitle();
00350 pNuMuBarQEffVsEn->SetLineColor(1);
00351 pNuMuBarQEffVsEn->SetFillColor(0);
00352
00353
00354 TProfile* pNuMuQEffVsEn=new TProfile("pNuMuQEffVsEn","pNuMuQEffVsEn",50,0,50);
00355 pNuMuQEffVsEn->SetTitle("Charge Sign Efficiency vs. True Neutrino Energy");
00356 pNuMuQEffVsEn->GetXaxis()->SetTitle("Energy (GeV)");
00357 pNuMuQEffVsEn->GetXaxis()->CenterTitle();
00358 pNuMuQEffVsEn->GetYaxis()->SetTitle("Efficiency");
00359 pNuMuQEffVsEn->GetYaxis()->CenterTitle();
00360 pNuMuQEffVsEn->SetLineColor(1);
00361 pNuMuQEffVsEn->SetFillColor(0);
00362
00363
00364 TProfile* pNuMuBarFitPassVsEn=new TProfile("pNuMuBarFitPassVsEn","pNuMuBarFitPassVsEn",50,0,50);
00365 pNuMuBarFitPassVsEn->SetTitle("Fit Pass Fraction vs. True Neutrino Energy");
00366 pNuMuBarFitPassVsEn->GetXaxis()->SetTitle("Energy (GeV)");
00367 pNuMuBarFitPassVsEn->GetXaxis()->CenterTitle();
00368 pNuMuBarFitPassVsEn->GetYaxis()->SetTitle("Pass Fraction");
00369 pNuMuBarFitPassVsEn->GetYaxis()->CenterTitle();
00370 pNuMuBarFitPassVsEn->SetLineColor(1);
00371 pNuMuBarFitPassVsEn->SetFillColor(0);
00372
00373
00374 TProfile* pNuMuFitPassVsEn=new TProfile("pNuMuFitPassVsEn","pNuMuFitPassVsEn",50,0,50);
00375 pNuMuFitPassVsEn->SetTitle("Fit Pass Fraction vs. True Neutrino Energy");
00376 pNuMuFitPassVsEn->GetXaxis()->SetTitle("Energy (GeV)");
00377 pNuMuFitPassVsEn->GetXaxis()->CenterTitle();
00378 pNuMuFitPassVsEn->GetYaxis()->SetTitle("Pass Fraction");
00379 pNuMuFitPassVsEn->GetYaxis()->CenterTitle();
00380 pNuMuFitPassVsEn->SetLineColor(1);
00381 pNuMuFitPassVsEn->SetFillColor(0);
00382
00383
00387
00388 this->InitialiseLoopVariables();
00389
00390
00391 for(Int_t entry=0;entry<this->GetEntries();entry++){
00392
00393 this->SetLoopVariables(entry);
00394
00395
00396 const NtpStRecord& ntp=(*this->GetNtpStRecord());
00397
00398 TClonesArray& mcTca=*ntp.mc;
00399 Int_t numInt=mcTca.GetEntries();
00400 MAXMSG("NuAnalysis",Msg::kInfo,20)
00401 <<"Number of entries in NtpMCTruth="
00402 <<mcTca.GetEntries()<<endl;
00403
00404
00405 TClonesArray& evtTca=(*ntp.evt);
00406 const Int_t numEvts=evtTca.GetEntriesFast();
00407
00408 TClonesArray& trkTca=(*ntp.trk);
00409 const Int_t numTrks=trkTca.GetEntriesFast();
00410
00411 TClonesArray& thtrkTca=(*ntp.thtrk);
00412 const Int_t numthtrks=thtrkTca.GetEntriesFast();
00413
00414 TClonesArray& thshwTca=(*ntp.thshw);
00415 const Int_t numthshws=thshwTca.GetEntriesFast();
00416
00417 TClonesArray& shwTca=(*ntp.shw);
00418 const Int_t numShws=shwTca.GetEntriesFast();
00419
00420
00421
00422
00423
00424
00425
00426 TClonesArray& hepTca=(*ntp.stdhep);
00427 const Int_t numHeps=hepTca.GetEntriesFast();
00428 MAXMSG("NuAnalysis",Msg::kInfo,1000)
00429 <<"Num stdhep entries="<<numHeps<<endl;
00430
00431 Int_t charmEvent=-1;
00432
00433
00434 for (Int_t ihep=0;ihep<numHeps;++ihep) {
00435 const NtpMCStdHep& stdhep=
00436 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
00437
00438 MAXMSG("NuAnalysis",Msg::kInfo,3000)
00439 <<"ihep="<<ihep
00440 <<", mc="<<stdhep.mc
00441 <<", id="<<stdhep.IdHEP
00442 <<", Ist="<<stdhep.IstHEP
00443 <<", parent=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
00444 <<", child=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
00445 <<", m="<<stdhep.mass
00446 <<", E="<<stdhep.p4[3]
00447 <<endl;
00448
00449 if (abs(stdhep.IdHEP)==411 || abs(stdhep.IdHEP)==421 ||
00450 abs(stdhep.IdHEP)==431 || abs(stdhep.IdHEP)==4122){
00451 charmEvent=stdhep.mc;
00452 MAXMSG("NuAnalysis",Msg::kInfo,3000)
00453 <<"Found Charm event, mc="<<charmEvent
00454 <<", id="<<stdhep.IdHEP<<endl;
00455 }
00456 }
00457
00458 if (charmEvent!=-1){
00459 MAXMSG("NuAnalysis",Msg::kInfo,3000)
00460 <<"CHARM EVENT, mc="<<charmEvent
00461 <<", entry="<<entry<<endl;
00462
00463 for (Int_t ihep=0;ihep<numHeps;++ihep) {
00464 const NtpMCStdHep& stdhep=
00465 *dynamic_cast<NtpMCStdHep*>(hepTca[ihep]);
00466
00467 if (stdhep.mc!=charmEvent) continue;
00468
00469 MAXMSG("NuAnalysis",Msg::kInfo,3000)
00470 <<"CH: i="<<ihep
00471 <<", mc="<<stdhep.mc
00472 <<", id="<<stdhep.IdHEP
00473 <<", Ist="<<stdhep.IstHEP
00474 <<", par=["<<stdhep.parent[0]<<","<<stdhep.parent[1]<<"]"
00475 <<", chi=["<<stdhep.child[0]<<","<<stdhep.child[1]<<"]"
00476 <<", m="<<stdhep.mass
00477 <<", E="<<stdhep.p4[3]
00478 <<endl;
00479 }
00480 }
00481
00482
00483 for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
00484 const NtpSREvent& evt=
00485 *dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
00486
00487
00488 if (evt.ntrack!=1) continue;
00489
00490 const NtpSRTrack& trk=
00491 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
00492
00493
00494 Float_t totalShwEn=0;
00495 for (Int_t ishw=0;ishw<evt.nshower;++ishw) {
00496 const NtpSRShower& shw=
00497 *dynamic_cast<NtpSRShower*>(shwTca[evt.shw[ishw]]);
00498 totalShwEn+=shw.ph.gev;
00499 }
00500
00501 MAXMSG("NuAnalysis",Msg::kInfo,100)
00502 <<"Entry "<<entry<<", event "<<ntpevt
00503 <<" has tracks="<<evt.ntrack
00504 <<", evts="<<numEvts
00505 <<", trks="<<numTrks<<", thtrks="<<numthtrks
00506 <<", shws="<<numShws<<", thshws="<<numthshws
00507 <<endl;
00508
00509
00510
00511
00512 const NtpTHTrack& thtrk=
00513 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
00514 MAXMSG("NuAnalysis",Msg::kInfo,100)
00515 <<"numInt="<<numInt<<", thtrk.neumc="<<thtrk.neumc<<endl;
00516
00517
00518
00519 const NtpMCTruth& mc=
00520 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
00521
00522 string iaction="NC";
00523 if (mc.iaction==1) iaction="CC";
00524 string iresonance="QE ";
00525 if (mc.iresonance==1002) iresonance="RES";
00526 else if (mc.iresonance==1003) iresonance="DIS";
00527 else if (mc.iresonance==1004) iresonance=
00528 "Coherent pion production";
00529 string inu="??";
00530 if (mc.inu==14) inu="NuMu ";
00531 else if (mc.inu==-14) inu="NuMuBar";
00532 else if (mc.inu==12) inu="NuE ";
00533 else if (mc.inu==-12) inu="NuEBar ";
00534
00535 Float_t trkCurveP=0;
00536 if (trk.momentum.qp!=0) trkCurveP=1./trk.momentum.qp;
00537
00538 MAXMSG("NuAnalysis",Msg::kInfo,100)
00539 <<""<<inu
00540 <<" "<<iaction
00541 <<" "<<iresonance
00542 <<", E="<<mc.p4neu[3]<<", mu="<<mc.p4mu1[3]<<" GeV"
00543 <<", Reco: trk="<<trk.momentum.range
00544 <<" (curv="<<trkCurveP<<")"
00545 <<", shw="<<totalShwEn
00546 <<endl;
00547
00548
00549 if (mc.iaction==0) continue;
00550
00551 Int_t chargeMC=0;
00552 if (mc.p4mu1[3]<0) chargeMC=-1;
00553 else if (mc.p4mu1[3]>0) chargeMC=+1;
00554
00555 Int_t charge=0;
00556 if (trk.momentum.qp<0) charge=-1;
00557 else if (trk.momentum.qp>0) charge=+1;
00558
00559
00560 if (mc.inu==14) {
00561 pNuMuFitPassVsEn->Fill(mc.p4neu[3],trk.fit.pass);
00562 }
00563 else if (mc.inu==-14) {
00564 pNuMuBarFitPassVsEn->Fill(mc.p4neu[3],trk.fit.pass);
00565 }
00566
00567 if (trk.fit.pass){
00568 if (mc.inu==14) {
00569 pNuMuQEffVsEn->Fill(mc.p4neu[3],charge==chargeMC);
00570 }
00571 else if (mc.inu==-14) {
00572 pNuMuBarQEffVsEn->Fill(mc.p4neu[3],charge==chargeMC);
00573 }
00574 }
00575
00576
00577
00578
00579
00580
00581
00582 }
00583 }
00584
00588
00589 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
00590
00591 MSG("NuAnalysis",Msg::kInfo)
00592 <<" ** Finished Efficiencies method **"<<endl;
00593 }
00594
00595
00596
00597 TSpline3* NuAnalysis::CreateTSpline3(std::string splName,
00598 std::string varNames,
00599 std::string fileName)
00600 {
00601 TNtuple nt(splName.c_str(),splName.c_str(),varNames.c_str());
00602 nt.ReadFile(fileName.c_str(),varNames.c_str());
00603
00604
00605 nt.Print();
00606
00607
00608
00609
00610
00611 nt.Draw("energy:xsec","","GOFF");
00612
00613 unsigned int n = (int) nt.GetSelectedRows();
00614 int * idx = new int[n];
00615 double * xx = new double[n];
00616 double * yy = new double[n];
00617
00618 TMath::Sort(n,nt.GetV1(),idx,false);
00619
00620 for(unsigned int i=0; i<n; i++) {
00621 int ii = idx[i];
00622 xx[i] = (nt.GetV1())[ii];
00623 yy[i] = (nt.GetV2())[ii];
00624
00625
00626 }
00627
00628 TSpline3 * spl=new TSpline3(("spl"+splName).c_str(),xx,yy,n,"0");
00629
00630 delete [] idx;
00631 delete [] xx;
00632 delete [] yy;
00633
00634 return spl;
00635 }
00636
00637
00638
00639 void NuAnalysis::CrossSections()
00640 {
00641 MSG("NuAnalysis",Msg::kInfo)
00642 <<" ** Running CrossSections method... **"<<endl;
00643
00644
00645 fOutFile=this->OpenFile(100,"NuCrossSections");
00646
00647 cout<<"here"<<endl;
00648 TSpline3* splnumupCC=CreateTSpline3
00649 ("numupCC","energy:xsec",
00650 "/home/hartnell/mytest/NuMuBar/numu-p-CC.txt");
00651 cout<<"here"<<endl;
00652 TSpline3* splnumunCC=CreateTSpline3
00653 ("numunCC","energy:xsec",
00654 "/home/hartnell/mytest/NuMuBar/numu-n-CC.txt");
00655 cout<<"here"<<endl;
00656 TSpline3* splnumubarpCC=CreateTSpline3
00657 ("numubarpCC","energy:xsec",
00658 "/home/hartnell/mytest/NuMuBar/numubar-p-CC.txt");
00659 cout<<"here"<<endl;
00660 TSpline3* splnumubarnCC=CreateTSpline3
00661 ("numubarnCC","energy:xsec",
00662 "/home/hartnell/mytest/NuMuBar/numubar-n-CC.txt");
00663 cout<<"hereend"<<endl;
00664
00665 TH1F* hXsecNuMuPCC=new TH1F("hXsecNuMuPCC","hXsecNuMuPCC",
00666 1000,0,100);
00667 hXsecNuMuPCC->SetTitle("NuMu Cross Section on Protons");
00668 hXsecNuMuPCC->GetXaxis()->SetTitle("Energy (GeV)");
00669 hXsecNuMuPCC->GetXaxis()->CenterTitle();
00670 hXsecNuMuPCC->GetYaxis()->SetTitle("#sigma (cm^{2})");
00671 hXsecNuMuPCC->GetYaxis()->CenterTitle();
00672
00673 TH1F* hXsecNuMuNCC=new TH1F("hXsecNuMuNCC","hXsecNuMuNCC",
00674 1000,0,100);
00675 hXsecNuMuNCC->SetTitle("NuMu Cross Section on Neutrons");
00676 hXsecNuMuNCC->GetXaxis()->SetTitle("Energy (GeV)");
00677 hXsecNuMuNCC->GetXaxis()->CenterTitle();
00678 hXsecNuMuNCC->GetYaxis()->SetTitle("#sigma (cm^{2})");
00679 hXsecNuMuNCC->GetYaxis()->CenterTitle();
00680
00681 TH1F* hXsecNuMuBarPCC=new TH1F("hXsecNuMuBarPCC","hXsecNuMuBarPCC",
00682 1000,0,100);
00683 hXsecNuMuBarPCC->SetTitle("NuMuBar Cross Section on Protons");
00684 hXsecNuMuBarPCC->GetXaxis()->SetTitle("Energy (GeV)");
00685 hXsecNuMuBarPCC->GetXaxis()->CenterTitle();
00686 hXsecNuMuBarPCC->GetYaxis()->SetTitle("#sigma (cm{^2})");
00687 hXsecNuMuBarPCC->GetYaxis()->CenterTitle();
00688
00689 TH1F* hXsecNuMuBarNCC=new TH1F("hXsecNuMuBarNCC","hXsecNuMuBarNCC",
00690 1000,0,100);
00691 hXsecNuMuBarNCC->SetTitle("NuMuBar Cross Section on Neutrons");
00692 hXsecNuMuBarNCC->GetXaxis()->SetTitle("Energy (GeV)");
00693 hXsecNuMuBarNCC->GetXaxis()->CenterTitle();
00694 hXsecNuMuBarNCC->GetYaxis()->SetTitle("#sigma (cm^{2})");
00695 hXsecNuMuBarNCC->GetYaxis()->CenterTitle();
00696
00697
00698 TH1F* hXsecNuMuFeCC=new TH1F("hXsecNuMuFeCC","hXsecNuMuFeCC",
00699 1000,0,100);
00700 hXsecNuMuFeCC->SetTitle("NuMu Cross Section on Fe");
00701 hXsecNuMuFeCC->GetXaxis()->SetTitle("Energy (GeV)");
00702 hXsecNuMuFeCC->GetXaxis()->CenterTitle();
00703 hXsecNuMuFeCC->GetYaxis()->SetTitle("#sigma (cm^{2})");
00704 hXsecNuMuFeCC->GetYaxis()->CenterTitle();
00705
00706 TH1F* hXsecNuMuBarFeCC=new TH1F("hXsecNuMuBarFeCC","hXsecNuMuBarFeCC",
00707 1000,0,100);
00708 hXsecNuMuBarFeCC->SetTitle("NuMuBar Cross Section on Fe");
00709 hXsecNuMuBarFeCC->GetXaxis()->SetTitle("Energy (GeV)");
00710 hXsecNuMuBarFeCC->GetXaxis()->CenterTitle();
00711 hXsecNuMuBarFeCC->GetYaxis()->SetTitle("#sigma (cm^{2})");
00712 hXsecNuMuBarFeCC->GetYaxis()->CenterTitle();
00713
00714
00715 for (Float_t i=0.05;i<100;i+=0.1){
00716 MAXMSG("NuAnalysis",Msg::kInfo,20)
00717 <<"E="<<i<<" GeV"<<endl
00718 <<" numu p CC xsec="<<splnumupCC->Eval(i)<<"e-38 cm^2"
00719 <<", numu n CC xsec="<<splnumunCC->Eval(i)<<"e-38 cm^2"<<endl
00720 <<" numubar p CC xsec="<<splnumubarpCC->Eval(i)<<"e-38 cm^2"
00721 <<", numubar n CC xsec="<<splnumubarnCC->Eval(i)<<"e-38 cm^2"
00722 <<endl;
00723
00724 Float_t xsecNuMuPCC=splnumupCC->Eval(i)<0 ? 0:splnumupCC->Eval(i);
00725 Float_t xsecNuMuNCC=splnumunCC->Eval(i)<0 ? 0:splnumunCC->Eval(i);
00726 Float_t xsecNuMuBarPCC=
00727 splnumubarpCC->Eval(i)<0 ? 0:splnumubarpCC->Eval(i);
00728 Float_t xsecNuMuBarNCC=
00729 splnumubarnCC->Eval(i)<0 ? 0:splnumubarnCC->Eval(i);
00730
00731 hXsecNuMuPCC->Fill(i,xsecNuMuPCC);
00732 hXsecNuMuNCC->Fill(i,xsecNuMuNCC);
00733 hXsecNuMuBarPCC->Fill(i,xsecNuMuBarPCC);
00734 hXsecNuMuBarNCC->Fill(i,xsecNuMuBarNCC);
00735
00736 Float_t FeA=55.845;
00737 Float_t FeZ=26;
00738 Float_t FeP=FeZ;
00739 Float_t FeN=FeA-FeZ;
00740
00741 Float_t xsecNuMuBarFeCC=xsecNuMuBarPCC*FeP+xsecNuMuBarNCC*FeN;
00742 Float_t xsecNuMuFeCC=xsecNuMuPCC*FeP+xsecNuMuNCC*FeN;
00743 hXsecNuMuFeCC->Fill(i,xsecNuMuFeCC);
00744 hXsecNuMuBarFeCC->Fill(i,xsecNuMuBarFeCC);
00745
00746
00747 }
00748
00749 MSG("NuAnalysis",Msg::kInfo)
00750 <<" ** Finished CrossSections method **"<<endl;
00751 }
00752
00753
00754
00755 void NuAnalysis::EnergySpectMC()
00756 {
00757 MSG("NuAnalysis",Msg::kInfo)
00758 <<" ** Running EnergySpectMC method... **"<<endl;
00759
00760
00761 fOutFile=this->OpenFile(100,"NuEnergySpectMC");
00762
00763 TH1F* hNuInt=new TH1F("hNuInt","hNuInt",10000,-100,1000);
00764 hNuInt->SetTitle("Number of neutrino interactions per spill");
00765 hNuInt->GetXaxis()->SetTitle("# interactions");
00766 hNuInt->GetXaxis()->CenterTitle();
00767 hNuInt->GetYaxis()->SetTitle("");
00768 hNuInt->GetYaxis()->CenterTitle();
00769 hNuInt->SetFillColor(0);
00770 hNuInt->SetLineColor(1);
00771
00772
00773 TH1F* hinuCC=new TH1F("hinuCC","hinuCC",80,-20,20);
00774 hinuCC->SetTitle("MC Number of Neutrino");
00775 hinuCC->GetXaxis()->SetTitle("MC Number");
00776 hinuCC->GetXaxis()->CenterTitle();
00777 hinuCC->GetYaxis()->SetTitle("");
00778 hinuCC->GetYaxis()->CenterTitle();
00779 hinuCC->SetFillColor(0);
00780 hinuCC->SetLineColor(1);
00781
00782
00783 TH1F* hinuNC=new TH1F("hinuNC","hinuNC",80,-20,20);
00784 hinuNC->SetTitle("MC Number of Neutrino");
00785 hinuNC->GetXaxis()->SetTitle("MC Number");
00786 hinuNC->GetXaxis()->CenterTitle();
00787 hinuNC->GetYaxis()->SetTitle("");
00788 hinuNC->GetYaxis()->CenterTitle();
00789 hinuNC->SetFillColor(0);
00790 hinuNC->SetLineColor(3);
00791 hinuNC->SetLineWidth(2);
00792
00793
00794
00795 TH1F* hNuMuBarEnCC=new TH1F("hNuMuBarEnCC","hNuMuBarEnCC",
00796 4*352,-32,320);
00797 hNuMuBarEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00798 hNuMuBarEnCC->GetXaxis()->CenterTitle();
00799 hNuMuBarEnCC->GetYaxis()->SetTitle("");
00800 hNuMuBarEnCC->GetYaxis()->CenterTitle();
00801 hNuMuBarEnCC->SetFillColor(0);
00802 hNuMuBarEnCC->SetLineColor(2);
00803
00804
00805 TH1F* hNuMuBarFidEnCC=new TH1F("hNuMuBarFidEnCC","hNuMuBarFidEnCC",
00806 4*352,-32,320);
00807 hNuMuBarFidEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00808 hNuMuBarFidEnCC->GetXaxis()->CenterTitle();
00809 hNuMuBarFidEnCC->GetYaxis()->SetTitle("");
00810 hNuMuBarFidEnCC->GetYaxis()->CenterTitle();
00811 hNuMuBarFidEnCC->SetFillColor(0);
00812 hNuMuBarFidEnCC->SetLineColor(2);
00813
00814
00815 TH1F* hNuMuBarLike=new TH1F("hNuMuBarLike","hNuMuBarLike",
00816 4*352,-32,320);
00817 hNuMuBarLike->GetXaxis()->SetTitle("Energy (GeV)");
00818 hNuMuBarLike->GetXaxis()->CenterTitle();
00819 hNuMuBarLike->GetYaxis()->SetTitle("");
00820 hNuMuBarLike->GetYaxis()->CenterTitle();
00821 hNuMuBarLike->SetFillColor(0);
00822 hNuMuBarLike->SetLineColor(2);
00823
00824
00825 TH1F* hNuMuEnCC=new TH1F("hNuMuEnCC","hNuMuEnCC",4*352,-32,320);
00826 hNuMuEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00827 hNuMuEnCC->GetXaxis()->CenterTitle();
00828 hNuMuEnCC->GetYaxis()->SetTitle("");
00829 hNuMuEnCC->GetYaxis()->CenterTitle();
00830 hNuMuEnCC->SetFillColor(0);
00831 hNuMuEnCC->SetLineColor(1);
00832 hNuMuEnCC->SetLineWidth(2);
00833
00834
00835 TH1F* hNuMuFidEnCC=new TH1F("hNuMuFidEnCC","hNuMuFidEnCC",
00836 4*352,-32,320);
00837 hNuMuFidEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00838 hNuMuFidEnCC->GetXaxis()->CenterTitle();
00839 hNuMuFidEnCC->GetYaxis()->SetTitle("");
00840 hNuMuFidEnCC->GetYaxis()->CenterTitle();
00841 hNuMuFidEnCC->SetFillColor(0);
00842 hNuMuFidEnCC->SetLineColor(1);
00843 hNuMuFidEnCC->SetLineWidth(2);
00844
00845
00846 TH1F* hNuMuLike=new TH1F("hNuMuLike","hNuMuLike",4*352,-32,320);
00847 hNuMuLike->GetXaxis()->SetTitle("Energy (GeV)");
00848 hNuMuLike->GetXaxis()->CenterTitle();
00849 hNuMuLike->GetYaxis()->SetTitle("");
00850 hNuMuLike->GetYaxis()->CenterTitle();
00851 hNuMuLike->SetFillColor(0);
00852 hNuMuLike->SetLineColor(1);
00853 hNuMuLike->SetLineWidth(2);
00854
00855
00856 TH1F* hNuMuBarEn4000CC=new TH1F("hNuMuBarEn4000CC","hNuMuBarEn4000CC",
00857 4000,-10,200);
00858 hNuMuBarEn4000CC->GetXaxis()->SetTitle("Energy (GeV)");
00859 hNuMuBarEn4000CC->GetXaxis()->CenterTitle();
00860 hNuMuBarEn4000CC->GetYaxis()->SetTitle("");
00861 hNuMuBarEn4000CC->GetYaxis()->CenterTitle();
00862 hNuMuBarEn4000CC->SetFillColor(0);
00863 hNuMuBarEn4000CC->SetLineColor(2);
00864
00865
00866 TH1F* hNuMuEn4000CC=new TH1F("hNuMuEn4000CC","hNuMuEn4000CC",
00867 4000,-10,200);
00868 hNuMuEn4000CC->GetXaxis()->SetTitle("Energy (GeV)");
00869 hNuMuEn4000CC->GetXaxis()->CenterTitle();
00870 hNuMuEn4000CC->GetYaxis()->SetTitle("");
00871 hNuMuEn4000CC->GetYaxis()->CenterTitle();
00872 hNuMuEn4000CC->SetFillColor(0);
00873 hNuMuEn4000CC->SetLineColor(1);
00874 hNuMuEn4000CC->SetLineWidth(2);
00875
00876
00877 TH1F* hNueBarEnCC=new TH1F("hNueBarEnCC","hNueBarEnCC",4*352,-32,320);
00878 hNueBarEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00879 hNueBarEnCC->GetXaxis()->CenterTitle();
00880 hNueBarEnCC->GetYaxis()->SetTitle("");
00881 hNueBarEnCC->GetYaxis()->CenterTitle();
00882 hNueBarEnCC->SetFillColor(0);
00883 hNueBarEnCC->SetLineColor(2);
00884 hNueBarEnCC->SetLineWidth(2);
00885 hNueBarEnCC->SetLineStyle(2);
00886
00887
00888 TH1F* hNueEnCC=new TH1F("hNueEnCC","hNueEnCC",4*352,-32,320);
00889 hNueEnCC->GetXaxis()->SetTitle("Energy (GeV)");
00890 hNueEnCC->GetXaxis()->CenterTitle();
00891 hNueEnCC->GetYaxis()->SetTitle("");
00892 hNueEnCC->GetYaxis()->CenterTitle();
00893 hNueEnCC->SetFillColor(0);
00894 hNueEnCC->SetLineColor(1);
00895 hNueEnCC->SetLineWidth(1);
00896 hNueEnCC->SetLineStyle(2);
00897
00898
00899
00900 TH1F* hNuMuBarEnNC=new TH1F("hNuMuBarEnNC","hNuMuBarEnNC",4*352,-32,320);
00901 hNuMuBarEnNC->GetXaxis()->SetTitle("Energy (GeV)");
00902 hNuMuBarEnNC->GetXaxis()->CenterTitle();
00903 hNuMuBarEnNC->GetYaxis()->SetTitle("");
00904 hNuMuBarEnNC->GetYaxis()->CenterTitle();
00905 hNuMuBarEnNC->SetFillColor(0);
00906 hNuMuBarEnNC->SetLineColor(4);
00907
00908
00909 TH1F* hNuMuEnNC=new TH1F("hNuMuEnNC","hNuMuEnNC",4*352,-32,320);
00910 hNuMuEnNC->GetXaxis()->SetTitle("Energy (GeV)");
00911 hNuMuEnNC->GetXaxis()->CenterTitle();
00912 hNuMuEnNC->GetYaxis()->SetTitle("");
00913 hNuMuEnNC->GetYaxis()->CenterTitle();
00914 hNuMuEnNC->SetFillColor(0);
00915 hNuMuEnNC->SetLineColor(1);
00916 hNuMuEnNC->SetLineWidth(3);
00917
00918
00919 TH1F* hNueBarEnNC=new TH1F("hNueBarEnNC","hNueBarEnNC",4*352,-32,320);
00920 hNueBarEnNC->GetXaxis()->SetTitle("Energy (GeV)");
00921 hNueBarEnNC->GetXaxis()->CenterTitle();
00922 hNueBarEnNC->GetYaxis()->SetTitle("");
00923 hNueBarEnNC->GetYaxis()->CenterTitle();
00924 hNueBarEnNC->SetFillColor(0);
00925 hNueBarEnNC->SetLineColor(4);
00926 hNueBarEnNC->SetLineWidth(2);
00927 hNueBarEnNC->SetLineStyle(2);
00928
00929
00930 TH1F* hNueEnNC=new TH1F("hNueEnNC","hNueEnNC",4*352,-32,320);
00931 hNueEnNC->GetXaxis()->SetTitle("Energy (GeV)");
00932 hNueEnNC->GetXaxis()->CenterTitle();
00933 hNueEnNC->GetYaxis()->SetTitle("");
00934 hNueEnNC->GetYaxis()->CenterTitle();
00935 hNueEnNC->SetFillColor(0);
00936 hNueEnNC->SetLineColor(3);
00937 hNueEnNC->SetLineWidth(1);
00938 hNueEnNC->SetLineStyle(2);
00939
00940
00941
00942 TH1F* hNuMuBarEnDepNC=new TH1F("hNuMuBarEnDepNC","hNuMuBarEnDepNC",4*352,-32,320);
00943 hNuMuBarEnDepNC->GetXaxis()->SetTitle("Energy (GeV)");
00944 hNuMuBarEnDepNC->GetXaxis()->CenterTitle();
00945 hNuMuBarEnDepNC->GetYaxis()->SetTitle("");
00946 hNuMuBarEnDepNC->GetYaxis()->CenterTitle();
00947 hNuMuBarEnDepNC->SetFillColor(0);
00948 hNuMuBarEnDepNC->SetLineColor(4);
00949
00950
00951 TH1F* hNuMuEnDepNC=new TH1F("hNuMuEnDepNC","hNuMuEnDepNC",4*352,-32,320);
00952 hNuMuEnDepNC->GetXaxis()->SetTitle("Energy (GeV)");
00953 hNuMuEnDepNC->GetXaxis()->CenterTitle();
00954 hNuMuEnDepNC->GetYaxis()->SetTitle("");
00955 hNuMuEnDepNC->GetYaxis()->CenterTitle();
00956 hNuMuEnDepNC->SetFillColor(0);
00957 hNuMuEnDepNC->SetLineColor(1);
00958 hNuMuEnDepNC->SetLineWidth(3);
00959
00960
00961 TH1F* hNueBarEnDepNC=new TH1F("hNueBarEnDepNC","hNueBarEnDepNC",4*352,-32,320);
00962 hNueBarEnDepNC->GetXaxis()->SetTitle("Energy (GeV)");
00963 hNueBarEnDepNC->GetXaxis()->CenterTitle();
00964 hNueBarEnDepNC->GetYaxis()->SetTitle("");
00965 hNueBarEnDepNC->GetYaxis()->CenterTitle();
00966 hNueBarEnDepNC->SetFillColor(0);
00967 hNueBarEnDepNC->SetLineColor(4);
00968 hNueBarEnDepNC->SetLineWidth(2);
00969 hNueBarEnDepNC->SetLineStyle(2);
00970
00971
00972 TH1F* hNueEnDepNC=new TH1F("hNueEnDepNC","hNueEnDepNC",4*352,-32,320);
00973 hNueEnDepNC->GetXaxis()->SetTitle("Energy (GeV)");
00974 hNueEnDepNC->GetXaxis()->CenterTitle();
00975 hNueEnDepNC->GetYaxis()->SetTitle("");
00976 hNueEnDepNC->GetYaxis()->CenterTitle();
00977 hNueEnDepNC->SetFillColor(0);
00978 hNueEnDepNC->SetLineColor(3);
00979 hNueEnDepNC->SetLineWidth(1);
00980 hNueEnDepNC->SetLineStyle(2);
00981
00982
00983
00984 TH1F* hNuYCC=new TH1F("hNuYCC","hNuYCC",100,-0.1,1.1);
00985 hNuYCC->GetXaxis()->SetTitle("Neutrino y");
00986 hNuYCC->GetXaxis()->CenterTitle();
00987 hNuYCC->GetYaxis()->SetTitle("");
00988 hNuYCC->GetYaxis()->CenterTitle();
00989 hNuYCC->SetFillColor(0);
00990 hNuYCC->SetLineColor(1);
00991 hNuYCC->SetLineWidth(1);
00992 hNuYCC->SetLineStyle(1);
00993
00994
00995 TH1F* hNuYNC=new TH1F("hNuYNC","hNuYNC",100,-0.1,1.1);
00996 hNuYNC->GetXaxis()->SetTitle("Neutrino y");
00997 hNuYNC->GetXaxis()->CenterTitle();
00998 hNuYNC->GetYaxis()->SetTitle("");
00999 hNuYNC->GetYaxis()->CenterTitle();
01000 hNuYNC->SetFillColor(0);
01001 hNuYNC->SetLineColor(3);
01002 hNuYNC->SetLineWidth(2);
01003 hNuYNC->SetLineStyle(1);
01004
01005
01006 TH1F* hNuYFidNuMuBarCC=new TH1F("hNuYFidNuMuBarCC","hNuYFidNuMuBarCC",
01007 100,-0.1,1.1);
01008 hNuYFidNuMuBarCC->GetXaxis()->SetTitle("Neutrino y");
01009 hNuYFidNuMuBarCC->GetXaxis()->CenterTitle();
01010 hNuYFidNuMuBarCC->GetYaxis()->SetTitle("");
01011 hNuYFidNuMuBarCC->GetYaxis()->CenterTitle();
01012 hNuYFidNuMuBarCC->SetFillColor(0);
01013 hNuYFidNuMuBarCC->SetLineColor(1);
01014 hNuYFidNuMuBarCC->SetLineWidth(1);
01015 hNuYFidNuMuBarCC->SetLineStyle(1);
01016
01017
01018 TH1F* hNuYFidNuMuBarNC=new TH1F("hNuYFidNuMuBarNC","hNuYFidNuMuBarNC",
01019 100,-0.1,1.1);
01020 hNuYFidNuMuBarNC->GetXaxis()->SetTitle("Neutrino y");
01021 hNuYFidNuMuBarNC->GetXaxis()->CenterTitle();
01022 hNuYFidNuMuBarNC->GetYaxis()->SetTitle("");
01023 hNuYFidNuMuBarNC->GetYaxis()->CenterTitle();
01024 hNuYFidNuMuBarNC->SetFillColor(0);
01025 hNuYFidNuMuBarNC->SetLineColor(3);
01026 hNuYFidNuMuBarNC->SetLineWidth(2);
01027 hNuYFidNuMuBarNC->SetLineStyle(1);
01028
01029
01030 TH1F* hNuYFidNuMuCC=new TH1F("hNuYFidNuMuCC","hNuYFidNuMuCC",
01031 100,-0.1,1.1);
01032 hNuYFidNuMuCC->GetXaxis()->SetTitle("Neutrino y");
01033 hNuYFidNuMuCC->GetXaxis()->CenterTitle();
01034 hNuYFidNuMuCC->GetYaxis()->SetTitle("");
01035 hNuYFidNuMuCC->GetYaxis()->CenterTitle();
01036 hNuYFidNuMuCC->SetFillColor(0);
01037 hNuYFidNuMuCC->SetLineColor(1);
01038 hNuYFidNuMuCC->SetLineWidth(1);
01039 hNuYFidNuMuCC->SetLineStyle(1);
01040
01041
01042 TH1F* hNuYFidNuMuNC=new TH1F("hNuYFidNuMuNC","hNuYFidNuMuNC",
01043 100,-0.1,1.1);
01044 hNuYFidNuMuNC->GetXaxis()->SetTitle("Neutrino y");
01045 hNuYFidNuMuNC->GetXaxis()->CenterTitle();
01046 hNuYFidNuMuNC->GetYaxis()->SetTitle("");
01047 hNuYFidNuMuNC->GetYaxis()->CenterTitle();
01048 hNuYFidNuMuNC->SetFillColor(0);
01049 hNuYFidNuMuNC->SetLineColor(3);
01050 hNuYFidNuMuNC->SetLineWidth(2);
01051 hNuYFidNuMuNC->SetLineStyle(1);
01052
01053
01054
01055 Int_t nbinsForY=60;
01056 vector<TH1F*> vNuYNuMuBarCC(nbinsForY);
01057 vector<TH1F*> vNuYNuMuCC(nbinsForY);
01058 vector<TH1F*> vNuYNuMuBarNC(nbinsForY);
01059 vector<TH1F*> vNuYNuMuNC(nbinsForY);
01060 for (Int_t i=0;i<nbinsForY;i++){
01061 Int_t en=(i+1)*5;
01062 string sEn=Form("%d",en);
01063 string s="hNuYNuMuBarCC"+sEn;
01064
01065 vNuYNuMuBarCC[i]=new TH1F(s.c_str(),s.c_str(),100,-0.1,1.1);
01066 vNuYNuMuBarCC[i]->GetXaxis()->SetTitle("Neutrino y");
01067 vNuYNuMuBarCC[i]->GetXaxis()->CenterTitle();
01068 vNuYNuMuBarCC[i]->GetYaxis()->SetTitle("");
01069 vNuYNuMuBarCC[i]->GetYaxis()->CenterTitle();
01070 vNuYNuMuBarCC[i]->SetFillColor(0);
01071 vNuYNuMuBarCC[i]->SetLineColor(2);
01072 vNuYNuMuBarCC[i]->SetLineWidth(2);
01073 vNuYNuMuBarCC[i]->SetLineStyle(2);
01074
01075
01076 s="hNuYNuMuCC"+sEn;
01077
01078 vNuYNuMuCC[i]=new TH1F(s.c_str(),s.c_str(),100,-0.1,1.1);
01079 vNuYNuMuCC[i]->GetXaxis()->SetTitle("Neutrino y");
01080 vNuYNuMuCC[i]->GetXaxis()->CenterTitle();
01081 vNuYNuMuCC[i]->GetYaxis()->SetTitle("");
01082 vNuYNuMuCC[i]->GetYaxis()->CenterTitle();
01083 vNuYNuMuCC[i]->SetFillColor(0);
01084 vNuYNuMuCC[i]->SetLineColor(1);
01085 vNuYNuMuCC[i]->SetLineWidth(2);
01086 vNuYNuMuCC[i]->SetLineStyle(1);
01087
01088
01089 s="hNuYNuMuBarNC"+sEn;
01090
01091 vNuYNuMuBarNC[i]=new TH1F(s.c_str(),s.c_str(),100,-0.1,1.1);
01092 vNuYNuMuBarNC[i]->GetXaxis()->SetTitle("Neutrino y");
01093 vNuYNuMuBarNC[i]->GetXaxis()->CenterTitle();
01094 vNuYNuMuBarNC[i]->GetYaxis()->SetTitle("");
01095 vNuYNuMuBarNC[i]->GetYaxis()->CenterTitle();
01096 vNuYNuMuBarNC[i]->SetFillColor(0);
01097 vNuYNuMuBarNC[i]->SetLineColor(2);
01098 vNuYNuMuBarNC[i]->SetLineWidth(2);
01099 vNuYNuMuBarNC[i]->SetLineStyle(2);
01100
01101
01102 s="hNuYNuMuNC"+sEn;
01103
01104 vNuYNuMuNC[i]=new TH1F(s.c_str(),s.c_str(),100,-0.1,1.1);
01105 vNuYNuMuNC[i]->GetXaxis()->SetTitle("Neutrino y");
01106 vNuYNuMuNC[i]->GetXaxis()->CenterTitle();
01107 vNuYNuMuNC[i]->GetYaxis()->SetTitle("");
01108 vNuYNuMuNC[i]->GetYaxis()->CenterTitle();
01109 vNuYNuMuNC[i]->SetFillColor(0);
01110 vNuYNuMuNC[i]->SetLineColor(1);
01111 vNuYNuMuNC[i]->SetLineWidth(2);
01112 vNuYNuMuNC[i]->SetLineStyle(1);
01113
01114 }
01115
01116 TH1F* hZCC=new TH1F("hZCC","hZCC",1000,-10,100);
01117 hZCC->GetXaxis()->SetTitle("Z (m)");
01118 hZCC->GetXaxis()->CenterTitle();
01119 hZCC->GetYaxis()->SetTitle("");
01120 hZCC->GetYaxis()->CenterTitle();
01121 hZCC->SetFillColor(0);
01122 hZCC->SetLineColor(1);
01123 hZCC->SetLineWidth(1);
01124 hZCC->SetLineStyle(1);
01125
01126
01127 TH1F* hZNC=new TH1F("hZNC","hZNC",1000,-10,100);
01128 hZNC->GetXaxis()->SetTitle("Z (m)");
01129 hZNC->GetXaxis()->CenterTitle();
01130 hZNC->GetYaxis()->SetTitle("");
01131 hZNC->GetYaxis()->CenterTitle();
01132 hZNC->SetFillColor(0);
01133 hZNC->SetLineColor(3);
01134 hZNC->SetLineWidth(2);
01135 hZNC->SetLineStyle(2);
01136
01137
01138 TH2F* hYvsXCC=new TH2F("hYvsXCC","hYvsXCC",
01139 100,-5,5,100,-5,5);
01140 hYvsXCC->SetTitle("Y vs X");
01141 hYvsXCC->GetXaxis()->SetTitle("X");
01142 hYvsXCC->GetXaxis()->CenterTitle();
01143 hYvsXCC->GetYaxis()->SetTitle("Y");
01144 hYvsXCC->GetYaxis()->CenterTitle();
01145 hYvsXCC->SetFillColor(0);
01146
01147
01148 TH2F* hYvsXNC=new TH2F("hYvsXNC","hYvsXNC",
01149 100,-5,5,100,-5,5);
01150 hYvsXNC->SetTitle("Y vs X");
01151 hYvsXNC->GetXaxis()->SetTitle("X");
01152 hYvsXNC->GetXaxis()->CenterTitle();
01153 hYvsXNC->GetYaxis()->SetTitle("Y");
01154 hYvsXNC->GetYaxis()->CenterTitle();
01155 hYvsXNC->SetFillColor(0);
01156
01157
01158 TH2F* hYvsXFidCC=new TH2F("hYvsXFidCC","hYvsXFidCC",
01159 100,-5,5,100,-5,5);
01160 hYvsXFidCC->SetTitle("Y vs X");
01161 hYvsXFidCC->GetXaxis()->SetTitle("X");
01162 hYvsXFidCC->GetXaxis()->CenterTitle();
01163 hYvsXFidCC->GetYaxis()->SetTitle("Y");
01164 hYvsXFidCC->GetYaxis()->CenterTitle();
01165 hYvsXFidCC->SetFillColor(0);
01166
01167
01168 TH2F* hYvsXFidNC=new TH2F("hYvsXFidNC","hYvsXFidNC",
01169 100,-5,5,100,-5,5);
01170 hYvsXFidNC->SetTitle("Y vs X");
01171 hYvsXFidNC->GetXaxis()->SetTitle("X");
01172 hYvsXFidNC->GetXaxis()->CenterTitle();
01173 hYvsXFidNC->GetYaxis()->SetTitle("Y");
01174 hYvsXFidNC->GetYaxis()->CenterTitle();
01175 hYvsXFidNC->SetFillColor(0);
01176
01177
01178 TH2F* hYvsX400CC=new TH2F("hYvsX400CC","hYvsX400CC",
01179 400,-20,20,400,-20,20);
01180 hYvsX400CC->SetTitle("Y vs X");
01181 hYvsX400CC->GetXaxis()->SetTitle("X");
01182 hYvsX400CC->GetXaxis()->CenterTitle();
01183 hYvsX400CC->GetYaxis()->SetTitle("Y");
01184 hYvsX400CC->GetYaxis()->CenterTitle();
01185 hYvsX400CC->SetFillColor(0);
01186
01187
01188 TH2F* hYvsX400NC=new TH2F("hYvsX400NC","hYvsX400NC",
01189 400,-20,20,400,-20,20);
01190 hYvsX400NC->SetTitle("Y vs X");
01191 hYvsX400NC->GetXaxis()->SetTitle("X");
01192 hYvsX400NC->GetXaxis()->CenterTitle();
01193 hYvsX400NC->GetYaxis()->SetTitle("Y");
01194 hYvsX400NC->GetYaxis()->CenterTitle();
01195 hYvsX400NC->SetFillColor(0);
01196
01197
01198
01199 TH2F* hYvsZCC=new TH2F("hYvsZCC","hYvsZCC",
01200 1000,-200,50,400,-20,20);
01201 hYvsZCC->SetTitle("Y vs Z");
01202 hYvsZCC->GetXaxis()->SetTitle("Z");
01203 hYvsZCC->GetXaxis()->CenterTitle();
01204 hYvsZCC->GetYaxis()->SetTitle("Y");
01205 hYvsZCC->GetYaxis()->CenterTitle();
01206 hYvsZCC->SetFillColor(0);
01207
01208
01209 TH2F* hYvsZNC=new TH2F("hYvsZNC","hYvsZNC",
01210 1000,-200,50,400,-20,20);
01211 hYvsZNC->SetTitle("Y vs Z");
01212 hYvsZNC->GetXaxis()->SetTitle("Z");
01213 hYvsZNC->GetXaxis()->CenterTitle();
01214 hYvsZNC->GetYaxis()->SetTitle("Y");
01215 hYvsZNC->GetYaxis()->CenterTitle();
01216 hYvsZNC->SetFillColor(0);
01217
01218
01219
01220 TH2F* hXvsZCC=new TH2F("hXvsZCC","hXvsZCC",
01221 1000,-200,50,400,-20,20);
01222 hXvsZCC->SetTitle("X vs Z");
01223 hXvsZCC->GetXaxis()->SetTitle("Z");
01224 hXvsZCC->GetXaxis()->CenterTitle();
01225 hXvsZCC->GetYaxis()->SetTitle("X");
01226 hXvsZCC->GetYaxis()->CenterTitle();
01227 hXvsZCC->SetFillColor(0);
01228
01229
01230 TH2F* hXvsZNC=new TH2F("hXvsZNC","hXvsZNC",
01231 1000,-200,50,400,-20,20);
01232 hXvsZNC->SetTitle("X vs Z");
01233 hXvsZNC->GetXaxis()->SetTitle("Z");
01234 hXvsZNC->GetXaxis()->CenterTitle();
01235 hXvsZNC->GetYaxis()->SetTitle("X");
01236 hXvsZNC->GetYaxis()->CenterTitle();
01237 hXvsZNC->SetFillColor(0);
01238
01239
01240 TH2F* hXvsZZoomCC=new TH2F("hXvsZZoomCC","hXvsZZoomCC",
01241 20000,-2,25,200,-3,4);
01242 hXvsZZoomCC->SetTitle("X vs Z");
01243 hXvsZZoomCC->GetXaxis()->SetTitle("Z");
01244 hXvsZZoomCC->GetXaxis()->CenterTitle();
01245 hXvsZZoomCC->GetYaxis()->SetTitle("X");
01246 hXvsZZoomCC->GetYaxis()->CenterTitle();
01247 hXvsZZoomCC->SetFillColor(0);
01248
01249
01250 TH2F* hXvsZZoomNC=new TH2F("hXvsZZoomNC","hXvsZZoomNC",
01251 20000,-2,25,200,-3,4);
01252 hXvsZZoomNC->SetTitle("X vs Z");
01253 hXvsZZoomNC->GetXaxis()->SetTitle("Z");
01254 hXvsZZoomNC->GetXaxis()->CenterTitle();
01255 hXvsZZoomNC->GetYaxis()->SetTitle("X");
01256 hXvsZZoomNC->GetYaxis()->CenterTitle();
01257 hXvsZZoomNC->SetFillColor(0);
01258
01259
01260
01261 TH3F* hXYZCC=new TH3F("hXYZCC","hXYZCC",
01262 100,-8,8, 400,-150,25, 200,-8,15);
01263 hXYZCC->SetTitle("XYZ");
01264 hXYZCC->GetXaxis()->SetTitle("X");
01265 hXYZCC->GetXaxis()->CenterTitle();
01266 hXYZCC->GetYaxis()->SetTitle("Y");
01267 hXYZCC->GetYaxis()->CenterTitle();
01268 hXYZCC->SetFillColor(0);
01269
01270
01271 TH3F* hXYZNC=new TH3F("hXYZNC","hXYZNC",
01272 100,-8,8, 400,-150,25, 200,-8,15);
01273 hXYZNC->SetTitle("XYZ");
01274 hXYZNC->GetXaxis()->SetTitle("X");
01275 hXYZNC->GetXaxis()->CenterTitle();
01276 hXYZNC->GetYaxis()->SetTitle("Y");
01277 hXYZNC->GetYaxis()->CenterTitle();
01278 hXYZNC->SetFillColor(0);
01279
01280
01281 TH1F* hTrksPerEvt=new TH1F("hTrksPerEvt","hTrksPerEvt",10,-1,9);
01282 hTrksPerEvt->SetTitle("Tracks per Reco'd Event");
01283 hTrksPerEvt->GetXaxis()->SetTitle("# tracks");
01284 hTrksPerEvt->GetXaxis()->CenterTitle();
01285 hTrksPerEvt->GetYaxis()->SetTitle("");
01286 hTrksPerEvt->GetYaxis()->CenterTitle();
01287 hTrksPerEvt->SetFillColor(0);
01288 hTrksPerEvt->SetLineColor(1);
01289 hTrksPerEvt->SetLineWidth(1);
01290 hTrksPerEvt->SetLineStyle(1);
01291
01292
01293 TH1F* hEvtsPerNu=new TH1F("hEvtsPerNu","hEvtsPerNu",10,-1,9);
01294 hEvtsPerNu->SetTitle("Events reco'd per neutrino interaction");
01295 hEvtsPerNu->GetXaxis()->SetTitle("Events reco'd");
01296 hEvtsPerNu->GetXaxis()->CenterTitle();
01297 hEvtsPerNu->GetYaxis()->SetTitle("");
01298 hEvtsPerNu->GetYaxis()->CenterTitle();
01299 hEvtsPerNu->SetFillColor(0);
01300 hEvtsPerNu->SetLineColor(1);
01301 hEvtsPerNu->SetLineWidth(1);
01302 hEvtsPerNu->SetLineStyle(1);
01303
01304
01305 TH1F* hEvtsPerNuCC=new TH1F("hEvtsPerNuCC","hEvtsPerNuCC",10,-1,9);
01306 hEvtsPerNuCC->SetTitle("Events reco'd per neutrino interaction");
01307 hEvtsPerNuCC->GetXaxis()->SetTitle("Events reco'd");
01308 hEvtsPerNuCC->GetXaxis()->CenterTitle();
01309 hEvtsPerNuCC->GetYaxis()->SetTitle("");
01310 hEvtsPerNuCC->GetYaxis()->CenterTitle();
01311 hEvtsPerNuCC->SetFillColor(0);
01312 hEvtsPerNuCC->SetLineColor(1);
01313 hEvtsPerNuCC->SetLineWidth(1);
01314 hEvtsPerNuCC->SetLineStyle(1);
01315
01316
01317 TH1F* hEvtsPerNuNC=new TH1F("hEvtsPerNuNC","hEvtsPerNuNC",10,-1,9);
01318 hEvtsPerNuNC->SetTitle("Events reco'd per neutrino interaction");
01319 hEvtsPerNuNC->GetXaxis()->SetTitle("Events reco'd");
01320 hEvtsPerNuNC->GetXaxis()->CenterTitle();
01321 hEvtsPerNuNC->GetYaxis()->SetTitle("");
01322 hEvtsPerNuNC->GetYaxis()->CenterTitle();
01323 hEvtsPerNuNC->SetFillColor(0);
01324 hEvtsPerNuNC->SetLineColor(3);
01325 hEvtsPerNuNC->SetLineWidth(2);
01326 hEvtsPerNuNC->SetLineStyle(1);
01327
01328
01329 TH1F* hEvtsPerSlc=new TH1F("hEvtsPerSlc","hEvtsPerSlc",10,-1,9);
01330 hEvtsPerSlc->SetTitle("Events reco'd per Slice");
01331 hEvtsPerSlc->GetXaxis()->SetTitle("Events reco'd");
01332 hEvtsPerSlc->GetXaxis()->CenterTitle();
01333 hEvtsPerSlc->GetYaxis()->SetTitle("");
01334 hEvtsPerSlc->GetYaxis()->CenterTitle();
01335 hEvtsPerSlc->SetFillColor(0);
01336 hEvtsPerSlc->SetLineColor(1);
01337 hEvtsPerSlc->SetLineWidth(2);
01338 hEvtsPerSlc->SetLineStyle(1);
01339
01340
01341 TH1F* hEvtsPerSlcg1=new TH1F("hEvtsPerSlcg1","hEvtsPerSlcg1",10,-1,9);
01342 hEvtsPerSlcg1->SetTitle("Events reco'd per Slice");
01343 hEvtsPerSlcg1->GetXaxis()->SetTitle("Events reco'd");
01344 hEvtsPerSlcg1->GetXaxis()->CenterTitle();
01345 hEvtsPerSlcg1->GetYaxis()->SetTitle("");
01346 hEvtsPerSlcg1->GetYaxis()->CenterTitle();
01347 hEvtsPerSlcg1->SetFillColor(0);
01348 hEvtsPerSlcg1->SetLineColor(2);
01349 hEvtsPerSlcg1->SetLineWidth(2);
01350 hEvtsPerSlcg1->SetLineStyle(1);
01351
01352
01353 TH1F* hEvtsPerSlcg1T=new TH1F("hEvtsPerSlcg1T","hEvtsPerSlcg1T",10,-1,9);
01354 hEvtsPerSlcg1T->SetTitle("Events reco'd per Slice");
01355 hEvtsPerSlcg1T->GetXaxis()->SetTitle("Events reco'd");
01356 hEvtsPerSlcg1T->GetXaxis()->CenterTitle();
01357 hEvtsPerSlcg1T->GetYaxis()->SetTitle("");
01358 hEvtsPerSlcg1T->GetYaxis()->CenterTitle();
01359 hEvtsPerSlcg1T->SetFillColor(0);
01360 hEvtsPerSlcg1T->SetLineColor(4);
01361 hEvtsPerSlcg1T->SetLineWidth(2);
01362 hEvtsPerSlcg1T->SetLineStyle(1);
01363
01364
01365 TH2F* hEvtsVsNuEnCC=new TH2F("hEvtsVsNuEnCC","hEvtsVsNuEnCC",
01366 110,-2,20,8,-1,7);
01367 hEvtsVsNuEnCC->SetTitle("Events Reco'd per Neutrino vs True Neutrino Energy");
01368 hEvtsVsNuEnCC->GetXaxis()->SetTitle("True Neutrino Energy (GeV)");
01369 hEvtsVsNuEnCC->GetXaxis()->CenterTitle();
01370 hEvtsVsNuEnCC->GetYaxis()->SetTitle("Events Reco'd");
01371 hEvtsVsNuEnCC->GetYaxis()->CenterTitle();
01372 hEvtsVsNuEnCC->SetFillColor(0);
01373
01374
01375 TH2F* hEvtsVsNuEnNC=new TH2F("hEvtsVsNuEnNC","hEvtsVsNuEnNC",
01376 110,-2,20,8,-1,7);
01377 hEvtsVsNuEnNC->SetTitle("Events Reco'd per Neutrino vs True Neutrino Energy");
01378 hEvtsVsNuEnNC->GetXaxis()->SetTitle("True Neutrino Energy (GeV)");
01379 hEvtsVsNuEnNC->GetXaxis()->CenterTitle();
01380 hEvtsVsNuEnNC->GetYaxis()->SetTitle("Events Reco'd");
01381 hEvtsVsNuEnNC->GetYaxis()->CenterTitle();
01382 hEvtsVsNuEnNC->SetFillColor(0);
01383
01384
01385
01386
01387 TProfile* pNuMuBarQEffVsEn=new TProfile("pNuMuBarQEffVsEn","pNuMuBarQEffVsEn",50,0,50);
01388 pNuMuBarQEffVsEn->SetTitle("Charge Sign Efficiency vs. True Neutrino Energy");
01389 pNuMuBarQEffVsEn->GetXaxis()->SetTitle("Energy (GeV)");
01390 pNuMuBarQEffVsEn->GetXaxis()->CenterTitle();
01391 pNuMuBarQEffVsEn->GetYaxis()->SetTitle("Efficiency");
01392 pNuMuBarQEffVsEn->GetYaxis()->CenterTitle();
01393 pNuMuBarQEffVsEn->SetLineColor(2);
01394 pNuMuBarQEffVsEn->SetFillColor(0);
01395
01396
01397 TProfile* pNuMuQEffVsEn=new TProfile("pNuMuQEffVsEn","pNuMuQEffVsEn",50,0,50);
01398 pNuMuQEffVsEn->SetTitle("Charge Sign Efficiency vs. True Neutrino Energy");
01399 pNuMuQEffVsEn->GetXaxis()->SetTitle("Energy (GeV)");
01400 pNuMuQEffVsEn->GetXaxis()->CenterTitle();
01401 pNuMuQEffVsEn->GetYaxis()->SetTitle("Efficiency");
01402 pNuMuQEffVsEn->GetYaxis()->CenterTitle();
01403 pNuMuQEffVsEn->SetLineColor(1);
01404 pNuMuQEffVsEn->SetFillColor(0);
01405
01406
01407 TProfile* pNuMuBarFitPassVsEn=new TProfile("pNuMuBarFitPassVsEn","pNuMuBarFitPassVsEn",50,0,50);
01408 pNuMuBarFitPassVsEn->SetTitle("Fit Pass Fraction vs. True Neutrino Energy");
01409 pNuMuBarFitPassVsEn->GetXaxis()->SetTitle("Energy (GeV)");
01410 pNuMuBarFitPassVsEn->GetXaxis()->CenterTitle();
01411 pNuMuBarFitPassVsEn->GetYaxis()->SetTitle("Pass Fraction");
01412 pNuMuBarFitPassVsEn->GetYaxis()->CenterTitle();
01413 pNuMuBarFitPassVsEn->SetLineColor(2);
01414 pNuMuBarFitPassVsEn->SetFillColor(0);
01415
01416
01417 TProfile* pNuMuFitPassVsEn=new TProfile("pNuMuFitPassVsEn","pNuMuFitPassVsEn",50,0,50);
01418 pNuMuFitPassVsEn->SetTitle("Fit Pass Fraction vs. True Neutrino Energy");
01419 pNuMuFitPassVsEn->GetXaxis()->SetTitle("Energy (GeV)");
01420 pNuMuFitPassVsEn->GetXaxis()->CenterTitle();
01421 pNuMuFitPassVsEn->GetYaxis()->SetTitle("Pass Fraction");
01422 pNuMuFitPassVsEn->GetYaxis()->CenterTitle();
01423 pNuMuFitPassVsEn->SetLineColor(1);
01424 pNuMuFitPassVsEn->SetFillColor(0);
01425
01426
01427
01428
01429 TH1F* hNuMuBarEnGoodQ=new TH1F("hNuMuBarEnGoodQ","hNuMuBarEnGoodQ",
01430 4*352,-32,320);
01431 hNuMuBarEnGoodQ->GetXaxis()->SetTitle("Energy (GeV)");
01432 hNuMuBarEnGoodQ->GetXaxis()->CenterTitle();
01433 hNuMuBarEnGoodQ->GetYaxis()->SetTitle("");
01434 hNuMuBarEnGoodQ->GetYaxis()->CenterTitle();
01435 hNuMuBarEnGoodQ->SetFillColor(0);
01436 hNuMuBarEnGoodQ->SetLineColor(2);
01437
01438
01439 TH1F* hNuMuEnGoodQ=new TH1F("hNuMuEnGoodQ","hNuMuEnGoodQ",4*352,-32,320);
01440 hNuMuEnGoodQ->GetXaxis()->SetTitle("Energy (GeV)");
01441 hNuMuEnGoodQ->GetXaxis()->CenterTitle();
01442 hNuMuEnGoodQ->GetYaxis()->SetTitle("");
01443 hNuMuEnGoodQ->GetYaxis()->CenterTitle();
01444 hNuMuEnGoodQ->SetFillColor(0);
01445 hNuMuEnGoodQ->SetLineColor(1);
01446 hNuMuEnGoodQ->SetLineWidth(2);
01447
01448
01449 TH1F* hNuMuBarEnBadQ=new TH1F("hNuMuBarEnBadQ","hNuMuBarEnBadQ",
01450 4*352,-32,320);
01451 hNuMuBarEnBadQ->GetXaxis()->SetTitle("Energy (GeV)");
01452 hNuMuBarEnBadQ->GetXaxis()->CenterTitle();
01453 hNuMuBarEnBadQ->GetYaxis()->SetTitle("");
01454 hNuMuBarEnBadQ->GetYaxis()->CenterTitle();
01455 hNuMuBarEnBadQ->SetFillColor(0);
01456 hNuMuBarEnBadQ->SetLineColor(2);
01457
01458
01459 TH1F* hNuMuEnBadQ=new TH1F("hNuMuEnBadQ","hNuMuEnBadQ",4*352,-32,320);
01460 hNuMuEnBadQ->GetXaxis()->SetTitle("Energy (GeV)");
01461 hNuMuEnBadQ->GetXaxis()->CenterTitle();
01462 hNuMuEnBadQ->GetYaxis()->SetTitle("");
01463 hNuMuEnBadQ->GetYaxis()->CenterTitle();
01464 hNuMuEnBadQ->SetFillColor(0);
01465 hNuMuEnBadQ->SetLineColor(1);
01466 hNuMuEnBadQ->SetLineWidth(2);
01467
01468
01469
01470
01471 TH1F* hNuMuBarRecoEnGoodQ=new TH1F("hNuMuBarRecoEnGoodQ","hNuMuBarRecoEnGoodQ",
01472 4*352,-32,320);
01473 hNuMuBarRecoEnGoodQ->GetXaxis()->SetTitle("Energy (GeV)");
01474 hNuMuBarRecoEnGoodQ->GetXaxis()->CenterTitle();
01475 hNuMuBarRecoEnGoodQ->GetYaxis()->SetTitle("");
01476 hNuMuBarRecoEnGoodQ->GetYaxis()->CenterTitle();
01477 hNuMuBarRecoEnGoodQ->SetFillColor(0);
01478 hNuMuBarRecoEnGoodQ->SetLineColor(2);
01479
01480
01481 TH1F* hNuMuRecoEnGoodQ=new TH1F("hNuMuRecoEnGoodQ","hNuMuRecoEnGoodQ",4*352,-32,320);
01482 hNuMuRecoEnGoodQ->GetXaxis()->SetTitle("Energy (GeV)");
01483 hNuMuRecoEnGoodQ->GetXaxis()->CenterTitle();
01484 hNuMuRecoEnGoodQ->GetYaxis()->SetTitle("");
01485 hNuMuRecoEnGoodQ->GetYaxis()->CenterTitle();
01486 hNuMuRecoEnGoodQ->SetFillColor(0);
01487 hNuMuRecoEnGoodQ->SetLineColor(1);
01488 hNuMuRecoEnGoodQ->SetLineWidth(2);
01489
01490
01491 TH1F* hNuMuBarRecoEnBadQ=new TH1F("hNuMuBarRecoEnBadQ","hNuMuBarRecoEnBadQ",
01492 4*352,-32,320);
01493 hNuMuBarRecoEnBadQ->GetXaxis()->SetTitle("Energy (GeV)");
01494 hNuMuBarRecoEnBadQ->GetXaxis()->CenterTitle();
01495 hNuMuBarRecoEnBadQ->GetYaxis()->SetTitle("");
01496 hNuMuBarRecoEnBadQ->GetYaxis()->CenterTitle();
01497 hNuMuBarRecoEnBadQ->SetFillColor(0);
01498 hNuMuBarRecoEnBadQ->SetLineColor(2);
01499
01500
01501 TH1F* hNuMuRecoEnBadQ=new TH1F("hNuMuRecoEnBadQ","hNuMuRecoEnBadQ",4*352,-32,320);
01502 hNuMuRecoEnBadQ->GetXaxis()->SetTitle("Energy (GeV)");
01503 hNuMuRecoEnBadQ->GetXaxis()->CenterTitle();
01504 hNuMuRecoEnBadQ->GetYaxis()->SetTitle("");
01505 hNuMuRecoEnBadQ->GetYaxis()->CenterTitle();
01506 hNuMuRecoEnBadQ->SetFillColor(0);
01507 hNuMuRecoEnBadQ->SetLineColor(1);
01508 hNuMuRecoEnBadQ->SetLineWidth(2);
01509
01510
01511
01512
01513 TH1F* hNuMuBarEnTrk0=new TH1F("hNuMuBarEnTrk0","hNuMuBarEnTrk0",
01514 4*352,-32,320);
01515 hNuMuBarEnTrk0->GetXaxis()->SetTitle("Energy (GeV)");
01516 hNuMuBarEnTrk0->GetXaxis()->CenterTitle();
01517 hNuMuBarEnTrk0->GetYaxis()->SetTitle("");
01518 hNuMuBarEnTrk0->GetYaxis()->CenterTitle();
01519 hNuMuBarEnTrk0->SetFillColor(0);
01520 hNuMuBarEnTrk0->SetLineColor(2);
01521
01522
01523 TH1F* hNuMuEnTrk0=new TH1F("hNuMuEnTrk0","hNuMuEnTrk0",4*352,-32,320);
01524 hNuMuEnTrk0->GetXaxis()->SetTitle("Energy (GeV)");
01525 hNuMuEnTrk0->GetXaxis()->CenterTitle();
01526 hNuMuEnTrk0->GetYaxis()->SetTitle("");
01527 hNuMuEnTrk0->GetYaxis()->CenterTitle();
01528 hNuMuEnTrk0->SetFillColor(0);
01529 hNuMuEnTrk0->SetLineColor(1);
01530 hNuMuEnTrk0->SetLineWidth(2);
01531
01532
01533
01534
01535 TH1F* hNuMuBarChi2GoodQ=new TH1F("hNuMuBarChi2GoodQ",
01536 "hNuMuBarChi2GoodQ",
01537 10*176,-16,160);
01538 hNuMuBarChi2GoodQ->GetXaxis()->SetTitle("Chi2/ndof");
01539 hNuMuBarChi2GoodQ->GetXaxis()->CenterTitle();
01540 hNuMuBarChi2GoodQ->GetYaxis()->SetTitle("");
01541 hNuMuBarChi2GoodQ->GetYaxis()->CenterTitle();
01542 hNuMuBarChi2GoodQ->SetFillColor(0);
01543 hNuMuBarChi2GoodQ->SetLineColor(2);
01544
01545
01546 TH1F* hNuMuChi2GoodQ=new TH1F("hNuMuChi2GoodQ","hNuMuChi2GoodQ",
01547 10*176,-16,160);
01548 hNuMuChi2GoodQ->GetXaxis()->SetTitle("Chi2/ndof");
01549 hNuMuChi2GoodQ->GetXaxis()->CenterTitle();
01550 hNuMuChi2GoodQ->GetYaxis()->SetTitle("");
01551 hNuMuChi2GoodQ->GetYaxis()->CenterTitle();
01552 hNuMuChi2GoodQ->SetFillColor(0);
01553 hNuMuChi2GoodQ->SetLineColor(1);
01554 hNuMuChi2GoodQ->SetLineWidth(2);
01555
01556
01557 TH1F* hNuMuBarChi2BadQ=new TH1F("hNuMuBarChi2BadQ",
01558 "hNuMuBarChi2BadQ",
01559 10*176,-16,160);
01560 hNuMuBarChi2BadQ->GetXaxis()->SetTitle("Chi2/ndof");
01561 hNuMuBarChi2BadQ->GetXaxis()->CenterTitle();
01562 hNuMuBarChi2BadQ->GetYaxis()->SetTitle("");
01563 hNuMuBarChi2BadQ->GetYaxis()->CenterTitle();
01564 hNuMuBarChi2BadQ->SetFillColor(0);
01565 hNuMuBarChi2BadQ->SetLineColor(2);
01566 hNuMuBarChi2BadQ->SetLineStyle(2);
01567
01568
01569 TH1F* hNuMuChi2BadQ=new TH1F("hNuMuChi2BadQ","hNuMuChi2BadQ",
01570 10*176,-16,160);
01571 hNuMuChi2BadQ->GetXaxis()->SetTitle("Chi2/ndof");
01572 hNuMuChi2BadQ->GetXaxis()->CenterTitle();
01573 hNuMuChi2BadQ->GetYaxis()->SetTitle("");
01574 hNuMuChi2BadQ->GetYaxis()->CenterTitle();
01575 hNuMuChi2BadQ->SetFillColor(0);
01576 hNuMuChi2BadQ->SetLineColor(1);
01577 hNuMuChi2BadQ->SetLineWidth(2);
01578 hNuMuChi2BadQ->SetLineStyle(2);
01579
01580
01581
01582
01583 TH1F* hNuMuBarSigqp_qpGoodQ=new TH1F("hNuMuBarSigqp_qpGoodQ",
01584 "hNuMuBarSigqp_qpGoodQ",
01585 5000,-100,100);
01586 hNuMuBarSigqp_qpGoodQ->GetXaxis()->SetTitle("#sigma(Q/P)/(Q/P)");
01587 hNuMuBarSigqp_qpGoodQ->GetXaxis()->CenterTitle();
01588 hNuMuBarSigqp_qpGoodQ->GetYaxis()->SetTitle("");
01589 hNuMuBarSigqp_qpGoodQ->GetYaxis()->CenterTitle();
01590 hNuMuBarSigqp_qpGoodQ->SetFillColor(0);
01591 hNuMuBarSigqp_qpGoodQ->SetLineColor(2);
01592
01593
01594 TH1F* hNuMuSigqp_qpGoodQ=new TH1F("hNuMuSigqp_qpGoodQ",
01595 "hNuMuSigqp_qpGoodQ",
01596 5000,-100,100);
01597 hNuMuSigqp_qpGoodQ->GetXaxis()->SetTitle("#sigma(Q/P)/(Q/P)");
01598 hNuMuSigqp_qpGoodQ->GetXaxis()->CenterTitle();
01599 hNuMuSigqp_qpGoodQ->GetYaxis()->SetTitle("");
01600 hNuMuSigqp_qpGoodQ->GetYaxis()->CenterTitle();
01601 hNuMuSigqp_qpGoodQ->SetFillColor(0);
01602 hNuMuSigqp_qpGoodQ->SetLineColor(1);
01603 hNuMuSigqp_qpGoodQ->SetLineWidth(2);
01604
01605
01606 TH1F* hNuMuBarSigqp_qpBadQ=new TH1F("hNuMuBarSigqp_qpBadQ",
01607 "hNuMuBarSigqp_qpBadQ",
01608 5000,-100,100);
01609 hNuMuBarSigqp_qpBadQ->GetXaxis()->SetTitle("#sigma(Q/P)/(Q/P)");
01610 hNuMuBarSigqp_qpBadQ->GetXaxis()->CenterTitle();
01611 hNuMuBarSigqp_qpBadQ->GetYaxis()->SetTitle("");
01612 hNuMuBarSigqp_qpBadQ->GetYaxis()->CenterTitle();
01613 hNuMuBarSigqp_qpBadQ->SetFillColor(0);
01614 hNuMuBarSigqp_qpBadQ->SetLineColor(2);
01615 hNuMuBarSigqp_qpBadQ->SetLineStyle(2);
01616
01617
01618 TH1F* hNuMuSigqp_qpBadQ=new TH1F("hNuMuSigqp_qpBadQ",
01619 "hNuMuSigqp_qpBadQ",
01620 5000,-100,100);
01621 hNuMuSigqp_qpBadQ->GetXaxis()->SetTitle("#sigma(Q/P)/(Q/P)");
01622 hNuMuSigqp_qpBadQ->GetXaxis()->CenterTitle();
01623 hNuMuSigqp_qpBadQ->GetYaxis()->SetTitle("");
01624 hNuMuSigqp_qpBadQ->GetYaxis()->CenterTitle();
01625 hNuMuSigqp_qpBadQ->SetFillColor(0);
01626 hNuMuSigqp_qpBadQ->SetLineColor(1);
01627 hNuMuSigqp_qpBadQ->SetLineWidth(2);
01628 hNuMuSigqp_qpBadQ->SetLineStyle(2);
01629
01630
01631 TH1F* hDeltaTime=new TH1F("hDeltaTime","hDeltaTime",10100,-100,10000);
01632 hDeltaTime->GetXaxis()->SetTitle("{#Delta}T (ns)");
01633 hDeltaTime->GetXaxis()->CenterTitle();
01634 hDeltaTime->GetYaxis()->SetTitle("");
01635 hDeltaTime->GetYaxis()->CenterTitle();
01636 hDeltaTime->SetFillColor(0);
01637 hDeltaTime->SetLineColor(1);
01638 hDeltaTime->SetLineWidth(2);
01639 hDeltaTime->SetLineStyle(2);
01640
01641
01642 TH1F* hDeltaT=new TH1F("hDeltaT","hDeltaT",10100,-100,10000);
01643 hDeltaT->GetXaxis()->SetTitle("#DeltaT (ns)");
01644 hDeltaT->GetXaxis()->CenterTitle();
01645 hDeltaT->GetYaxis()->SetTitle("");
01646 hDeltaT->GetYaxis()->CenterTitle();
01647 hDeltaT->SetFillColor(0);
01648 hDeltaT->SetLineColor(1);
01649 hDeltaT->SetLineWidth(2);
01650 hDeltaT->SetLineStyle(1);
01651
01652
01653 TH1F* hDeltaTg1=new TH1F("hDeltaTg1","hDeltaTg1",10100,-100,10000);
01654 hDeltaTg1->GetXaxis()->SetTitle("#DeltaT (ns)");
01655 hDeltaTg1->GetXaxis()->CenterTitle();
01656 hDeltaTg1->GetYaxis()->SetTitle("");
01657 hDeltaTg1->GetYaxis()->CenterTitle();
01658 hDeltaTg1->SetFillColor(0);
01659 hDeltaTg1->SetLineColor(2);
01660 hDeltaTg1->SetLineWidth(2);
01661 hDeltaTg1->SetLineStyle(1);
01662
01663
01664 Int_t totalSnarlsCounter=0;
01665 Int_t ccInFidCounter=0;
01666 Int_t recoEvtCounter=0;
01667 Int_t nuIntCounter=0;
01668 Int_t oneRecoEvtCounter=0;
01669 Int_t zeroTrackCounter=0;
01670 Int_t zeroTrackNuMuCounter=0;
01671 Int_t zeroTrackNuMuBarCounter=0;
01672 Int_t twoTrackCounter=0;
01673 Int_t twoTrackNuMuCounter=0;
01674 Int_t twoTrackNuMuBarCounter=0;
01675 Int_t threeTrackCounter=0;
01676 Int_t gr3TrackCounter=0;
01677 Int_t oneTrackCounter=0;
01678 Int_t passFitCounter=0;
01679 Int_t badNuMuCounter=0;
01680 Int_t goodNuMuCounter=0;
01681 Int_t badNuMuBarCounter=0;
01682 Int_t goodNuMuBarCounter=0;
01683
01684 const NuCuts cuts;
01685 const NuReco reco;
01686 const NuExtraction ext;
01687
01688 NuConfig config;
01689 this->ExtractConfig(config);
01690
01691 const NuLibrary& lib=NuLibrary::Instance();
01692
01696
01697 this->InitialiseLoopVariables();
01698
01699
01700 for(Int_t entry=0;entry<this->GetEntries();entry++){
01701
01702 this->SetLoopVariables(entry);
01703
01704 set<Int_t> setNuInFid;
01705 set<Int_t> setNuInFidCC;
01706 map<Int_t,Double_t> deltaTs;
01707 map<Int_t,Int_t> evtsPerSlc;
01708
01709
01710 const NtpStRecord& ntp=(*this->GetNtpStRecord());
01711
01712
01713 const RecCandHeader& rec=ntp.GetHeader();
01714 MAXMSG("NuAnalysis",Msg::kInfo,5)
01715 <<"Found: run="<<rec.GetRun()
01716 <<", subrun="<<rec.GetSubRun()
01717 <<", detector="<<rec.GetVldContext().GetDetector()
01718 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
01719 <<endl;
01720
01721 VldTimeStamp vldts;
01722
01723 VldContext vc(rec.GetVldContext().GetDetector(),
01724 rec.GetVldContext().GetSimFlag(),vldts);
01725
01726 UgliGeomHandle ugh(vc);
01727
01728 TClonesArray& mcTca=*ntp.mc;
01729 Int_t numInt=mcTca.GetEntries();
01730 MAXMSG("NuAnalysis",Msg::kInfo,20)
01731 <<"Number of entries in NtpMCTruth="
01732 <<mcTca.GetEntries()<<endl;
01733
01734
01735 hNuInt->Fill(mcTca.GetEntries());
01736
01737 for (Int_t i=0;i<numInt;i++){
01738 const NtpMCTruth& mc=*(dynamic_cast<NtpMCTruth*>(mcTca[i]));
01739
01740 TVector3 xyz(mc.vtxx,mc.vtxy,mc.vtxz);
01741
01742 TVector3 uvz=ugh.xyz2uvz(xyz);
01743
01744
01745 NuEvent nu;
01746
01747
01748 ext.ExtractGeneralInfo(ntp,nu);
01749
01750 Bool_t isInFidVol=cuts.IsInFidVol(mc.vtxx,mc.vtxy,mc.vtxz,
01751 uvz.X(),uvz.Y(),
01752 0,0,
01753 nu.detector,nu.anaVersion,
01754 nu.releaseType,nu.simFlag);
01755
01756 if (isInFidVol) {
01757 deltaTs.clear();
01758 reco.GetEvtDeltaTs(ntp,deltaTs);
01759 evtsPerSlc.clear();
01760 reco.GetEvtsPerSlc(ntp,evtsPerSlc);
01761 MAXMSG("NuAnalysis",Msg::kInfo,100)
01762 <<"Contained mc="<<mc.index<<endl;
01763 setNuInFid.insert(mc.index);
01764 if (mc.iaction==1 && (mc.inu==-14 || mc.inu==14)) {
01765 setNuInFidCC.insert(mc.index);
01766 }
01767 }
01768
01769 MAXMSG("NuAnalysis",Msg::kDebug,20)
01770 <<"isInFidVol="<<isInFidVol
01771 <<", x="<<mc.vtxx<<", y="<<mc.vtxy<<", z="<<mc.vtxz
01772 <<", u="<<uvz.X()<<", v="<<uvz.Y()<<endl;
01773
01774 if (isInFidVol){
01775 MAXMSG("NuAnalysis",Msg::kInfo,20)
01776 <<"isInFidVol="<<isInFidVol
01777 <<", x="<<mc.vtxx<<", y="<<mc.vtxy<<", z="<<mc.vtxz
01778 <<", u="<<uvz.X()<<", v="<<uvz.Y()<<endl;
01779 }
01780
01781
01782 Int_t enIndex=static_cast<Int_t>(mc.p4neu[3]*2.0+1);
01783
01784 if (mc.iaction==1){
01785 hNuYCC->Fill(mc.y);
01786 if (isInFidVol) {
01787 hYvsXFidCC->Fill(mc.vtxx,mc.vtxy);
01788 }
01789 hYvsXCC->Fill(mc.vtxx,mc.vtxy);
01790 hYvsX400CC->Fill(mc.vtxx,mc.vtxy);
01791 hYvsZCC->Fill(mc.vtxz,mc.vtxy);
01792 hXvsZCC->Fill(mc.vtxz,mc.vtxx);
01793 if (mc.vtxy>-2.5 && mc.vtxy<3) hXvsZZoomCC->
01794 Fill(mc.vtxz,mc.vtxx);
01795 hXYZCC->Fill(mc.vtxx,mc.vtxz,mc.vtxy);
01796 hZCC->Fill(mc.vtxz);
01797 hinuCC->Fill(mc.inu);
01798 if (mc.inu==-14){
01799 if (enIndex<nbinsForY) vNuYNuMuBarCC[enIndex]->Fill(mc.y);
01800 hNuMuBarEnCC->Fill(mc.p4neu[3]);
01801 hNuMuBarEn4000CC->Fill(mc.p4neu[3]);
01802 if (isInFidVol) {
01803 hNuYFidNuMuBarCC->Fill(mc.y);
01804 MAXMSG("NuAnalysis",Msg::kInfo,200)
01805 <<"enIndex="<<enIndex<<", E="<<mc.p4neu[3]<<endl;
01806 hNuMuBarFidEnCC->Fill(mc.p4neu[3]);
01807 hNuMuBarLike->Fill(mc.p4neu[3]);
01808 }
01809 }
01810 else if (mc.inu==14) {
01811 if (enIndex<nbinsForY) vNuYNuMuCC[enIndex]->Fill(mc.y);
01812 hNuMuEnCC->Fill(mc.p4neu[3]);
01813 hNuMuEn4000CC->Fill(mc.p4neu[3]);
01814 if (isInFidVol) {
01815 hNuYFidNuMuCC->Fill(mc.y);
01816 hNuMuFidEnCC->Fill(mc.p4neu[3]);
01817 hNuMuLike->Fill(mc.p4neu[3]);
01818 }
01819 }
01820 else if (mc.inu==-12) hNueBarEnCC->Fill(mc.p4neu[3]);
01821 else if (mc.inu==12) hNueEnCC->Fill(mc.p4neu[3]);
01822 else {
01823 MAXMSG("NuAnalysis",Msg::kWarning,100)
01824 <<"particle not defined for filling histo="<<mc.inu<<endl;
01825 }
01826 }
01827 else if (mc.iaction==0){
01828 hNuYNC->Fill(mc.y);
01829 if (isInFidVol) {
01830 hYvsXFidNC->Fill(mc.vtxx,mc.vtxy);
01831 }
01832 hYvsXNC->Fill(mc.vtxx,mc.vtxy);
01833 hYvsX400NC->Fill(mc.vtxx,mc.vtxy);
01834 hYvsZNC->Fill(mc.vtxz,mc.vtxy);
01835 hXvsZNC->Fill(mc.vtxz,mc.vtxx);
01836
01837 if (mc.vtxy>-2.5 && mc.vtxy<3) hXvsZZoomNC->
01838 Fill(mc.vtxz,mc.vtxx);
01839 hXYZNC->Fill(mc.vtxx,mc.vtxz,mc.vtxy);
01840 hZNC->Fill(mc.vtxz);
01841 hinuNC->Fill(mc.inu);
01842 if (mc.inu==-14) {
01843 if (enIndex<nbinsForY) vNuYNuMuBarNC[enIndex]->Fill(mc.y);
01844 hNuMuBarEnNC->Fill(mc.p4neu[3]);
01845 hNuMuBarEnDepNC->Fill(mc.y*mc.p4neu[3]);
01846 if (isInFidVol) {
01847 hNuYFidNuMuBarNC->Fill(mc.y);
01848 }
01849 }
01850 else if (mc.inu==14) {
01851 if (enIndex<nbinsForY) vNuYNuMuNC[enIndex]->Fill(mc.y);
01852 hNuMuEnNC->Fill(mc.p4neu[3]);
01853 hNuMuEnDepNC->Fill(mc.y*mc.p4neu[3]);
01854 if (isInFidVol) {
01855 hNuYFidNuMuNC->Fill(mc.y);
01856 }
01857 }
01858 else if (mc.inu==-12) {
01859 hNueBarEnNC->Fill(mc.p4neu[3]);
01860 hNueBarEnDepNC->Fill(mc.y*mc.p4neu[3]);
01861 }
01862 else if (mc.inu==12) {
01863 hNueEnNC->Fill(mc.p4neu[3]);
01864 hNueEnDepNC->Fill(mc.y*mc.p4neu[3]);
01865 }
01866 else {
01867 MAXMSG("NuAnalysis",Msg::kWarning,100)
01868 <<"particle not defined for filling histo="<<mc.inu<<endl;
01869 }
01870 }
01871 else cout<<"Ahhh, bad iaction"<<endl;
01872 }
01873
01874 MAXMSG("NuAnalysis",Msg::kInfo,100)
01875 <<"Number of fiducial volume interactions in spill="
01876 <<setNuInFid.size()<<endl;
01877
01878 TClonesArray& thevtTca=(*ntp.thevt);
01879 const Int_t numthevts=thevtTca.GetEntriesFast();
01880 TClonesArray& evtTca=(*ntp.evt);
01881
01882
01883
01884 set<Int_t> evtsInFid;
01885 set<Int_t> evtsInFidCC;
01886
01887 set<Int_t>::iterator setNuInFidEnd=setNuInFid.end();
01888 map<Int_t,Int_t> nuEvtCounter;
01889 map<Int_t,Int_t> nuEvtCounterCC;
01890 map<Int_t,multiset<Double_t> > nuIntEvtTimesCC;
01891 for (Int_t ithevt=0;ithevt<numthevts;ithevt++){
01892 const NtpTHEvent& thevt=
01893 *dynamic_cast<NtpTHEvent*>(thevtTca[ithevt]);
01894 const NtpSREvent& evt=
01895 *dynamic_cast<NtpSREvent*>(evtTca[ithevt]);
01896
01897 if (setNuInFid.find(thevt.neumc)!=setNuInFidEnd){
01898 MAXMSG("NuAnalysis",Msg::kInfo,100)
01899 <<"Found mc="<<thevt.neumc<<endl;
01900
01901
01902 nuEvtCounter[thevt.neumc]++;
01903
01904
01905 evtsInFid.insert(ithevt);
01906
01907 const NtpMCTruth& mc=
01908 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
01909 if (mc.iaction==1 && (mc.inu==-14 || mc.inu==14)) {
01910 nuEvtCounterCC[thevt.neumc]++;
01911 evtsInFidCC.insert(ithevt);
01912
01913
01914
01915 Double_t time=reco.GetEvtMedianTime(ntp,evt);
01916 nuIntEvtTimesCC[thevt.neumc].insert(time);
01917 }
01918 }
01919 }
01920
01921
01922 typedef set<Int_t>::iterator setNuInFidIt;
01923 for (setNuInFidIt it=setNuInFid.begin();
01924 it!=setNuInFid.end();++it){
01925
01926 MAXMSG("NuAnalysis",Msg::kInfo,100)
01927 <<"Looking in evt counter for mc.index="<<*it
01928 <<", evtCo size="<<nuEvtCounter.size()
01929 <<", setNu size="<<setNuInFid.size()
01930 <<endl;
01931 if (nuEvtCounter.find(*it)!=nuEvtCounter.end()){
01932 MAXMSG("NuAnalysis",Msg::kInfo,100)
01933 <<" Found mc.index="<<*it<<endl;
01934 }
01935 else {
01936 MAXMSG("NuAnalysis",Msg::kInfo,30)
01937 <<" NOT found mc.index="<<*it
01938 <<", evtCo size="<<nuEvtCounter.size()
01939 <<", setNu size="<<setNuInFid.size()
01940 <<endl;
01941
01942 nuEvtCounter[*it]=0;
01943 }
01944 }
01945
01946
01947 typedef map<Int_t,Int_t>::iterator evCounterIt;
01948 for (evCounterIt it=nuEvtCounter.begin();
01949 it!=nuEvtCounter.end();++it){
01950
01951 const NtpMCTruth& mc=
01952 *(dynamic_cast<NtpMCTruth*>(mcTca[it->first]));
01953
01954
01955 if (!(mc.inu==-14 || mc.inu==14)) continue;
01956
01957 hEvtsPerNu->Fill(it->second);
01958 if (mc.iaction==1){
01959 hEvtsPerNuCC->Fill(it->second);
01960 hEvtsVsNuEnCC->Fill(mc.p4neu[3],it->second);
01961 }
01962 else if (mc.iaction==0){
01963 hEvtsPerNuNC->Fill(it->second);
01964 hEvtsVsNuEnNC->Fill(mc.p4neu[3],it->second);
01965 }
01966 }
01967
01971
01972
01973 totalSnarlsCounter++;
01974
01975 if (setNuInFidCC.size()<1) continue;
01976 ccInFidCounter++;
01977
01978 if (evtsInFidCC.size()<1) continue;
01979 recoEvtCounter++;
01980
01981
01982 TClonesArray& trkTca=(*ntp.trk);
01983
01984
01985
01986
01987 typedef set<Int_t>::iterator setEvtsInFidIt;
01988 for (setEvtsInFidIt it=evtsInFidCC.begin();
01989 it!=evtsInFidCC.end();++it) {
01990 const NtpTHEvent& thevt=
01991 *dynamic_cast<NtpTHEvent*>(thevtTca[*it]);
01992 const NtpSREvent& evt=
01993 *dynamic_cast<NtpSREvent*>(evtTca[*it]);
01994 const NtpMCTruth& mc=
01995 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
01996
01997
01998 nuIntCounter++;
01999
02000
02001 if (mc.iaction==0) cout<<"Ahhh NC"<<endl;
02002
02003 if (!(mc.inu==-14 || mc.inu==14)) cout<<"Ahhh nue"<<endl;
02004
02005
02006
02007 Int_t numEvtsPerNuInt=-1;
02008 map<Int_t,Int_t>::iterator countIt=
02009 nuEvtCounterCC.find(thevt.neumc);
02010 if (countIt!=nuEvtCounterCC.end()){
02011 numEvtsPerNuInt=countIt->second;
02012
02013
02014 Int_t times=nuIntEvtTimesCC.find(thevt.neumc)->second.size();
02015 if (times>1 || numEvtsPerNuInt>1){
02016 MAXMSG("NuAnalysis",Msg::kInfo,200)
02017 <<"numEvtsPerNuInt="<<numEvtsPerNuInt
02018 <<", times="<<times<<endl;
02019
02020 typedef map<Int_t,multiset<Double_t> >::iterator nuIntIt;
02021 nuIntIt nuInt=nuIntEvtTimesCC.find(thevt.neumc);
02022 typedef multiset<Double_t>::iterator timeIt;
02023 timeIt beg=nuInt->second.begin();
02024 timeIt end=nuInt->second.end();
02025 timeIt lastTimeIt=nuInt->second.end();
02026 for (timeIt it=beg;it!=end;++it){
02027 MsgFormat ffmt("%9.9f");
02028
02029 MAXMSG("NuAnalysis",Msg::kDebug,200)
02030 <<" in multiset time="<<ffmt(*it)
02031 <<", evt.slc="<<evt.slc<<endl;
02032 if (lastTimeIt!=end){
02033 Double_t deltaTime=(*it)-(*lastTimeIt);
02034 MAXMSG("NuAnalysis",Msg::kDebug,200)
02035 <<" deltaTime="<<ffmt(deltaTime)<<endl;
02036 hDeltaTime->Fill(deltaTime/(Munits::ns));
02037 }
02038
02039
02040 lastTimeIt=it;
02041 }
02042
02043 }
02044
02045
02046
02047
02048
02049
02050
02051 }
02052 else cout<<"Can't find mc="<<thevt.neumc<<endl;
02053
02054 if (numEvtsPerNuInt>1) {
02055 const NtpSREvent& evt=
02056 *dynamic_cast<NtpSREvent*>(evtTca[*it]);
02057 const NtpMCTruth& mc=
02058 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
02059 MAXMSG("NuAnalysis",Msg::kDebug,200)
02060 <<"entry="<<entry<<", NuInt: mc="<<thevt.neumc
02061 <<", slc="<<evt.slc
02062 <<", evt="<<*it<<" ("<<numEvtsPerNuInt<<")"
02063 <<", ntrk="<<evt.ntrack
02064 <<", nshw="<<evt.nshower
02065 <<", Nu: E="<<mc.p4neu[3]
02066 <<", y="<<mc.y
02067 <<endl;
02068
02069
02070 if (deltaTs.find(evt.index)!=deltaTs.end()){
02071 hDeltaTg1->Fill(deltaTs[evt.index]/(Munits::ns));
02072
02073 if (deltaTs[evt.index]/(Munits::ns)>50){
02074 hEvtsPerSlcg1T->Fill(evtsPerSlc[evt.slc]);
02075 }
02076 }
02077 else cout<<"No delta T for evt="<<evt.index<<endl;
02078 hEvtsPerSlcg1->Fill(evtsPerSlc[evt.slc]);
02079 }
02080 else {
02081 if (deltaTs.find(evt.index)!=deltaTs.end()){
02082 hDeltaT->Fill(deltaTs[evt.index]/(Munits::ns));
02083 }
02084 else cout<<"No delta T for evt="<<evt.index<<endl;
02085 hEvtsPerSlc->Fill(evtsPerSlc[evt.slc]);
02086 }
02087
02088
02089 if (numEvtsPerNuInt!=1) continue;
02090
02091
02092 oneRecoEvtCounter++;
02093
02094
02095
02096
02097
02098
02099 hTrksPerEvt->Fill(evt.ntrack);
02100
02102
02104 if (evt.ntrack==0) {
02105 zeroTrackCounter++;
02106 if (mc.inu==14) {
02107 zeroTrackNuMuCounter++;
02108 hNuMuEnTrk0->Fill(mc.p4neu[3]);
02109
02110
02111 }
02112 else if (mc.inu==-14) {
02113 zeroTrackNuMuBarCounter++;
02114 hNuMuBarEnTrk0->Fill(mc.p4neu[3]);
02115
02116
02117 hNuMuLike->Fill(mc.p4neu[3]);
02118 }
02119 }
02120 if (evt.ntrack==1) {
02121 oneTrackCounter++;
02122 }
02123 if (evt.ntrack==2) {
02124 twoTrackCounter++;
02125 if (mc.inu==14) {
02126 twoTrackNuMuCounter++;
02127 }
02128 else if (mc.inu==-14) {
02129 twoTrackNuMuBarCounter++;
02130 }
02131 }
02132 if (evt.ntrack==3) {
02133 threeTrackCounter++;
02134 }
02135 if (evt.ntrack>3) {
02136 gr3TrackCounter++;
02137 }
02138
02139
02140
02141 Int_t trkToUse=0;
02142 Float_t trkLength=0;
02143 if (evt.ntrack>1){
02144 for(Int_t itrk=0;itrk<evt.ntrack;itrk++){
02145 const NtpSRTrack& trk=
02146 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[itrk]]);
02147
02148 MAXMSG("NuAnalysis",Msg::kDebug,200)
02149 <<"trk="<<itrk<<"/"<<evt.ntrack
02150 <<", pass="<<trk.fit.pass
02151 <<", ntrklike="<<trk.plane.ntrklike
02152 <<endl;
02153
02154 if (!trk.fit.pass) continue;
02155
02156
02157 if (trk.plane.ntrklike>trkLength) {
02158 trkLength=trk.plane.ntrklike;
02159 trkToUse=itrk;
02160 }
02161 }
02162 MAXMSG("NuAnalysis",Msg::kDebug,200)
02163 <<"selected trkToUse="<<trkToUse<<endl;
02164 }
02165
02166
02167 NuEvent nu;
02168
02169
02170 lib.ext.ExtractGeneralInfo(ntp,nu);
02171
02172 lib.ext.ExtractEvtInfo(evt,nu);
02173
02174 lib.ext.ExtractTrkInfo(ntp,evt,nu);
02175
02176 lib.ext.ExtractShwInfo(ntp,evt,nu);
02177
02179
02180 lib.reco.GetEvtEnergy(nu);
02181
02182 Float_t recoEn=nu.energy;
02183 MAXMSG("NuAnalysis",Msg::kInfo,500)
02184 <<"Reco'd event energy="<<recoEn<<endl;
02185
02186
02187 if (evt.ntrack<1) continue;
02188
02189
02190 const NtpSRTrack& trk=
02191 *dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[trkToUse]]);
02192
02193 Int_t chargeMC=0;
02194 if (mc.p4mu1[3]<0) chargeMC=-1;
02195 else if (mc.p4mu1[3]>0) chargeMC=+1;
02196
02197 Int_t charge=0;
02198 if (trk.momentum.qp<0) charge=-1;
02199 else if (trk.momentum.qp>0) charge=+1;
02200
02201
02202 if (mc.inu==14) {
02203 pNuMuFitPassVsEn->Fill(mc.p4neu[3],trk.fit.pass);
02204 }
02205 else if (mc.inu==-14) {
02206 pNuMuBarFitPassVsEn->Fill(mc.p4neu[3],trk.fit.pass);
02207 }
02208
02209 if (trk.fit.pass){
02210 Float_t chi2PerNdof=-1;
02211 if (trk.fit.ndof) chi2PerNdof=trk.fit.chi2/trk.fit.ndof;
02212 Float_t sigqp_qp=0;
02213 if (trk.momentum.qp) sigqp_qp=trk.momentum.eqp/trk.momentum.qp;
02214
02215 passFitCounter++;
02216 if (mc.inu==14) {
02217 pNuMuQEffVsEn->Fill(mc.p4neu[3],charge==chargeMC);
02218 if (charge!=chargeMC){
02219 badNuMuCounter++;
02220 hNuMuEnBadQ->Fill(mc.p4neu[3]);
02221 hNuMuRecoEnBadQ->Fill(recoEn);
02222
02223 hNuMuBarLike->Fill(mc.p4neu[3]);
02224 hNuMuChi2BadQ->Fill(chi2PerNdof);
02225 hNuMuSigqp_qpBadQ->Fill(sigqp_qp);
02226 }
02227 else{
02228 goodNuMuCounter++;
02229 hNuMuEnGoodQ->Fill(mc.p4neu[3]);
02230 hNuMuRecoEnGoodQ->Fill(recoEn);
02231 hNuMuChi2GoodQ->Fill(chi2PerNdof);
02232 hNuMuSigqp_qpGoodQ->Fill(sigqp_qp);
02233 }
02234 }
02235 else if (mc.inu==-14) {
02236 pNuMuBarQEffVsEn->Fill(mc.p4neu[3],charge==chargeMC);
02237 if (charge!=chargeMC){
02238 badNuMuBarCounter++;
02239 hNuMuBarEnBadQ->Fill(mc.p4neu[3]);
02240 hNuMuBarRecoEnBadQ->Fill(recoEn);
02241
02242 hNuMuLike->Fill(mc.p4neu[3]);
02243 hNuMuBarChi2BadQ->Fill(chi2PerNdof);
02244 hNuMuBarSigqp_qpBadQ->Fill(sigqp_qp);
02245 }
02246 else{
02247 goodNuMuBarCounter++;
02248 hNuMuBarEnGoodQ->Fill(mc.p4neu[3]);
02249 hNuMuBarRecoEnGoodQ->Fill(recoEn);
02250 hNuMuBarChi2GoodQ->Fill(chi2PerNdof);
02251 hNuMuBarSigqp_qpGoodQ->Fill(sigqp_qp);
02252 }
02253 }
02254 }
02255 }
02256 }
02257
02261
02262 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
02263
02264 MSG("NuAnalysis",Msg::kInfo)
02265 <<endl<<" *** Summary ***"<<endl
02266 <<"Total snarls ="<<totalSnarlsCounter<<endl
02267 <<"Snarls w/ CC NuInt in Fid="<<ccInFidCounter<<endl
02268 <<"Snarls w/ reco'd fid evt ="<<recoEvtCounter<<endl
02269 <<endl
02270 <<"Total neutrino interactions="<<nuIntCounter<<endl
02271 <<"Only 1 reco'd evt per nu ="<<oneRecoEvtCounter<<endl
02272 <<"Only 1 track ="<<oneTrackCounter<<endl
02273 <<"Pass track fit ="<<passFitCounter<<endl
02274 <<" bad Q NuMu ="<<badNuMuCounter<<endl
02275 <<" good Q NuMu ="<<goodNuMuCounter<<endl
02276 <<" bad Q NuMuBar ="<<badNuMuBarCounter<<endl
02277 <<" good Q NuMuBar ="<<goodNuMuBarCounter<<endl;
02278
02279 MSG("NuAnalysis",Msg::kInfo)
02280 <<endl<<" *** CC NuMu Event Track Breakdown ***"<<endl
02281 <<"zeroTrackCounter ="<<zeroTrackCounter<<endl
02282 <<" NuMu ="<<zeroTrackNuMuCounter<<endl
02283 <<" NuMuBar="<<zeroTrackNuMuBarCounter<<endl
02284 <<"oneTrackCounter ="<<oneTrackCounter<<endl
02285 <<"twoTrackCounter ="<<twoTrackCounter<<endl
02286 <<" NuMu ="<<twoTrackNuMuCounter<<endl
02287 <<" NuMuBar="<<twoTrackNuMuBarCounter<<endl
02288 <<"threeTrackCounter="<<threeTrackCounter<<endl
02289 <<"g3TrackCounter ="<<gr3TrackCounter<<endl;
02290
02291 MSG("NuAnalysis",Msg::kInfo)
02292 <<" ** Finished EnergySpectMC method **"<<endl;
02293 }
02294
02295
02296
02297 void NuAnalysis::N_1()
02298 {
02299 MSG("NuAnalysis",Msg::kInfo)
02300 <<" ** Running N_1 method... **"<<endl;
02301
02302
02303 fOutFile=this->OpenFile(100,"NuMuBarN_1");
02304
02305 TH1F* hDpID=new TH1F("hDpID","hDpID",4*160,-1.6,1.6);
02306 hDpID->GetXaxis()->SetTitle("PID (from MadDpID)");
02307 hDpID->GetXaxis()->CenterTitle();
02308 hDpID->GetYaxis()->SetTitle("");
02309 hDpID->GetYaxis()->CenterTitle();
02310 hDpID->SetFillColor(0);
02311 hDpID->SetLineColor(1);
02312 hDpID->SetLineWidth(3);
02313
02314
02315 TH1F* hDpIDN_1=new TH1F("hDpIDN_1","hDpIDN_1",4*160,-1.6,1.6);
02316 hDpIDN_1->GetXaxis()->SetTitle("PID (from MadDpID)");
02317 hDpIDN_1->GetXaxis()->CenterTitle();
02318 hDpIDN_1->GetYaxis()->SetTitle("");
02319 hDpIDN_1->GetYaxis()->CenterTitle();
02320 hDpIDN_1->SetFillColor(0);
02321 hDpIDN_1->SetLineColor(2);
02322 hDpIDN_1->SetLineWidth(3);
02323
02324
02325 TH1F* hDpIDN=new TH1F("hDpIDN","hDpIDN",4*160,-1.6,1.6);
02326 hDpIDN->GetXaxis()->SetTitle("PID (from MadDpID)");
02327 hDpIDN->GetXaxis()->CenterTitle();
02328 hDpIDN->GetYaxis()->SetTitle("");
02329 hDpIDN->GetYaxis()->CenterTitle();
02330 hDpIDN->SetFillColor(0);
02331 hDpIDN->SetLineColor(4);
02332 hDpIDN->SetLineWidth(3);
02333
02334
02335
02336
02337 TH1F* hSigmaqp_qp=new TH1F("hSigmaqp_qp","hSigmaqp_qp",200*160,-160,160);
02338 hSigmaqp_qp->GetXaxis()->SetTitle("Sigmaqp_qp");
02339 hSigmaqp_qp->GetXaxis()->CenterTitle();
02340 hSigmaqp_qp->GetYaxis()->SetTitle("");
02341 hSigmaqp_qp->GetYaxis()->CenterTitle();
02342 hSigmaqp_qp->SetFillColor(0);
02343 hSigmaqp_qp->SetLineColor(1);
02344 hSigmaqp_qp->SetLineWidth(3);
02345
02346
02347 TH1F* hSigmaqp_qpN_1=new TH1F("hSigmaqp_qpN_1","hSigmaqp_qpN_1",200*160,-160,160);
02348 hSigmaqp_qpN_1->GetXaxis()->SetTitle("Sigmaqp_qp");
02349 hSigmaqp_qpN_1->GetXaxis()->CenterTitle();
02350 hSigmaqp_qpN_1->GetYaxis()->SetTitle("");
02351 hSigmaqp_qpN_1->GetYaxis()->CenterTitle();
02352 hSigmaqp_qpN_1->SetFillColor(0);
02353 hSigmaqp_qpN_1->SetLineColor(2);
02354 hSigmaqp_qpN_1->SetLineWidth(3);
02355
02356
02357 TH1F* hSigmaqp_qpN=new TH1F("hSigmaqp_qpN","hSigmaqp_qpN",200*160,-160,160);
02358 hSigmaqp_qpN->GetXaxis()->SetTitle("Sigmaqp_qp");
02359 hSigmaqp_qpN->GetXaxis()->CenterTitle();
02360 hSigmaqp_qpN->GetYaxis()->SetTitle("");
02361 hSigmaqp_qpN->GetYaxis()->CenterTitle();
02362 hSigmaqp_qpN->SetFillColor(0);
02363 hSigmaqp_qpN->SetLineColor(4);
02364 hSigmaqp_qpN->SetLineWidth(3);
02365
02366
02367
02368
02369 TH1F* hChi2=new TH1F("hChi2","hChi2",10*176,-16,160);
02370 hChi2->GetXaxis()->SetTitle("Chi2");
02371 hChi2->GetXaxis()->CenterTitle();
02372 hChi2->GetYaxis()->SetTitle("");
02373 hChi2->GetYaxis()->CenterTitle();
02374 hChi2->SetFillColor(0);
02375 hChi2->SetLineColor(1);
02376 hChi2->SetLineWidth(3);
02377
02378
02379 TH1F* hChi2N_1=new TH1F("hChi2N_1","hChi2N_1",10*176,-16,160);
02380 hChi2N_1->GetXaxis()->SetTitle("Chi2");
02381 hChi2N_1->GetXaxis()->CenterTitle();
02382 hChi2N_1->GetYaxis()->SetTitle("");
02383 hChi2N_1->GetYaxis()->CenterTitle();
02384 hChi2N_1->SetFillColor(0);
02385 hChi2N_1->SetLineColor(2);
02386 hChi2N_1->SetLineWidth(3);
02387
02388
02389 TH1F* hChi2N=new TH1F("hChi2N","hChi2N",10*176,-16,160);
02390 hChi2N->GetXaxis()->SetTitle("Chi2");
02391 hChi2N->GetXaxis()->CenterTitle();
02392 hChi2N->GetYaxis()->SetTitle("");
02393 hChi2N->GetYaxis()->CenterTitle();
02394 hChi2N->SetFillColor(0);
02395 hChi2N->SetLineColor(4);
02396 hChi2N->SetLineWidth(3);
02397
02398
02399
02400
02401 TH1F* hFitProb=new TH1F("hFitProb","hFitProb",10*160,-1.6,1.6);
02402 hFitProb->GetXaxis()->SetTitle("FitProb");
02403 hFitProb->GetXaxis()->CenterTitle();
02404 hFitProb->GetYaxis()->SetTitle("");
02405 hFitProb->GetYaxis()->CenterTitle();
02406 hFitProb->SetFillColor(0);
02407 hFitProb->SetLineColor(1);
02408 hFitProb->SetLineWidth(3);
02409
02410
02411 TH1F* hFitProbN_1=new TH1F("hFitProbN_1","hFitProbN_1",10*160,-1.6,1.6);
02412 hFitProbN_1->GetXaxis()->SetTitle("FitProb");
02413 hFitProbN_1->GetXaxis()->CenterTitle();
02414 hFitProbN_1->GetYaxis()->SetTitle("");
02415 hFitProbN_1->GetYaxis()->CenterTitle();
02416 hFitProbN_1->SetFillColor(0);
02417 hFitProbN_1->SetLineColor(2);
02418 hFitProbN_1->SetLineWidth(3);
02419
02420
02421 TH1F* hFitProbN=new TH1F("hFitProbN","hFitProbN",10*160,-1.6,1.6);
02422 hFitProbN->GetXaxis()->SetTitle("FitProb");
02423 hFitProbN->GetXaxis()->CenterTitle();
02424 hFitProbN->GetYaxis()->SetTitle("");
02425 hFitProbN->GetYaxis()->CenterTitle();
02426 hFitProbN->SetFillColor(0);
02427 hFitProbN->SetLineColor(4);
02428 hFitProbN->SetLineWidth(3);
02429
02430
02431
02432
02433 TH1F* hDeltaT=new TH1F("hDeltaT","hDeltaT",2010,-10,2000);
02434 hDeltaT->GetXaxis()->SetTitle("Delta T (ns)");
02435 hDeltaT->GetXaxis()->CenterTitle();
02436 hDeltaT->GetYaxis()->SetTitle("");
02437 hDeltaT->GetYaxis()->CenterTitle();
02438 hDeltaT->SetFillColor(0);
02439 hDeltaT->SetLineColor(1);
02440 hDeltaT->SetLineWidth(3);
02441
02442
02443 TH1F* hDeltaTN_1=new TH1F("hDeltaTN_1","hDeltaTN_1",2010,-10,2000);
02444 hDeltaTN_1->GetXaxis()->SetTitle("Delta T (ns)");
02445 hDeltaTN_1->GetXaxis()->CenterTitle();
02446 hDeltaTN_1->GetYaxis()->SetTitle("");
02447 hDeltaTN_1->GetYaxis()->CenterTitle();
02448 hDeltaTN_1->SetFillColor(0);
02449 hDeltaTN_1->SetLineColor(2);
02450 hDeltaTN_1->SetLineWidth(3);
02451
02452
02453 TH1F* hDeltaTN=new TH1F("hDeltaTN","hDeltaTN",2010,-10,2000);
02454 hDeltaTN->GetXaxis()->SetTitle("Delta T (ns)");
02455 hDeltaTN->GetXaxis()->CenterTitle();
02456 hDeltaTN->GetYaxis()->SetTitle("");
02457 hDeltaTN->GetYaxis()->CenterTitle();
02458 hDeltaTN->SetFillColor(0);
02459 hDeltaTN->SetLineColor(4);
02460 hDeltaTN->SetLineWidth(3);
02461
02462
02463
02464
02465 TH1F* hEvtsPerSlc=new TH1F("hEvtsPerSlc","hEvtsPerSlc",60,-10,50);
02466 hEvtsPerSlc->GetXaxis()->SetTitle("Evts Per Slc");
02467 hEvtsPerSlc->GetXaxis()->CenterTitle();
02468 hEvtsPerSlc->GetYaxis()->SetTitle("");
02469 hEvtsPerSlc->GetYaxis()->CenterTitle();
02470 hEvtsPerSlc->SetFillColor(0);
02471 hEvtsPerSlc->SetLineColor(1);
02472 hEvtsPerSlc->SetLineWidth(3);
02473
02474
02475 TH1F* hEvtsPerSlcN_1=new TH1F("hEvtsPerSlcN_1","hEvtsPerSlcN_1",
02476 60,-10,50);
02477 hEvtsPerSlcN_1->GetXaxis()->SetTitle("Evts Per Slc");
02478 hEvtsPerSlcN_1->GetXaxis()->CenterTitle();
02479 hEvtsPerSlcN_1->GetYaxis()->SetTitle("");
02480 hEvtsPerSlcN_1->GetYaxis()->CenterTitle();
02481 hEvtsPerSlcN_1->SetFillColor(0);
02482 hEvtsPerSlcN_1->SetLineColor(2);
02483 hEvtsPerSlcN_1->SetLineWidth(3);
02484
02485
02486 TH1F* hEvtsPerSlcN=new TH1F("hEvtsPerSlcN","hEvtsPerSlcN",60,-10,50);
02487 hEvtsPerSlcN->GetXaxis()->SetTitle("Evts Per Slc");
02488 hEvtsPerSlcN->GetXaxis()->CenterTitle();
02489 hEvtsPerSlcN->GetYaxis()->SetTitle("");
02490 hEvtsPerSlcN->GetYaxis()->CenterTitle();
02491 hEvtsPerSlcN->SetFillColor(0);
02492 hEvtsPerSlcN->SetLineColor(4);
02493 hEvtsPerSlcN->SetLineWidth(3);
02494
02495
02496
02497 const NuCuts cuts;
02498
02499 const NuPlots plots;
02500 const NuReco reco;
02501 const NuExtraction ext;
02502 const NuBeam beam;
02503
02504 NuConfig config;
02505 this->ExtractConfig(config);
02506
02507
02508
02509 NuMadAnalysis& mad=*new NuMadAnalysis();
02510 MadDpID& dp=*new MadDpID();
02511
02512 BeamType::BeamType_t beamType=static_cast
02513 <BeamType::BeamType_t>(config.beamType);
02514
02515
02516 map<Int_t,Double_t> deltaTs;
02517 map<Int_t,Int_t> evtsPerSlc;
02518
02519 Int_t goodSpillCounter=0;
02520 Int_t badSpillCounter=0;
02521 Float_t totalPot=0;
02522 Float_t totalPotBad=0;
02523 Int_t evtCounter=0;
02524 Int_t evtInFidVolCounter=0;
02525 Int_t evtWithTrkCounter=0;
02526 Int_t goodBestTrkCounter=0;
02527 Int_t goodTrkPassCounter=0;
02528 Int_t goodRecoEnCounter=0;
02529 Int_t goodUVDiffCounter=0;
02530
02531
02532
02533
02534
02535
02536 const NuLibrary& lib=NuLibrary::Instance();
02537
02541
02542 this->InitialiseLoopVariables();
02543
02544
02545 for(Int_t entry=0;entry<this->GetEntries();entry++){
02546
02547 this->SetLoopVariables(entry);
02548
02549
02550 const NtpStRecord& ntp=(*this->GetNtpStRecord());
02551
02552 const RecCandHeader& rec=ntp.GetHeader();
02553 MAXMSG("NuAnalysis",Msg::kInfo,5)
02554 <<"Found: run="<<rec.GetRun()
02555 <<", subrun="<<rec.GetSubRun()
02556 <<", detector="<<rec.GetVldContext().GetDetector()
02557 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
02558 <<endl;
02559
02560
02561 plots.FillGeneralHistos(ntp);
02562
02563
02564 NuEvent nuSnarl;
02565
02566
02567 ext.ExtractGeneralInfo(ntp,nuSnarl);
02568
02569 if (!beam.IsGoodSpillAndFillPot
02570 (this->GetNtpBDLiteRecord(),config,nuSnarl)) {
02571 badSpillCounter++;
02572 totalPotBad+=nuSnarl.pot;
02573 continue;
02574 }
02575 goodSpillCounter++;
02576 totalPot+=nuSnarl.pot;
02577
02578 dp.ChoosePDFs(rec.GetVldContext().GetDetector(),beamType);
02579 mad.SetEntry(&ntp);
02580
02581 TClonesArray& evtTca=(*ntp.evt);
02582 const Int_t numEvts=evtTca.GetEntriesFast();
02583
02584
02585 for (Int_t ievt=0;ievt<numEvts;ievt++){
02586 const NtpSREvent& evt=
02587 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
02588 evtCounter++;
02589
02590
02591 NuEvent nu;
02592
02593
02594 ext.ExtractGeneralInfo(ntp,nu);
02595
02596 lib.ext.ExtractEvtInfo(evt,nu);
02597
02598 lib.ext.ExtractTrkInfo(ntp,evt,nu);
02599
02600 lib.ext.ExtractShwInfo(ntp,evt,nu);
02601
02602 Bool_t isInFidVol=cuts.IsInFidVol(evt.vtx.x,evt.vtx.y,
02603 evt.vtx.z,
02604 evt.vtx.u,evt.vtx.v,
02605 0,0,
02606 nu.detector,nu.anaVersion,
02607 nu.releaseType,nu.simFlag);
02608
02609 if (!isInFidVol) continue;
02610 evtInFidVolCounter++;
02611
02612 nu.dpID=dp.CalcPID(&mad,ievt,0);
02613 MAXMSG("NuAnalysis",Msg::kDebug,50)
02614 <<"dpID="<<nu.dpID<<endl;
02615
02616
02617
02618 if (evt.ntrack!=1) continue;
02619 evtWithTrkCounter++;
02620
02621
02622 Int_t bestTrack=lib.reco.GetBestTrack(nu);
02623 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
02624 (ntp,evt,bestTrack-1);
02625 const NtpSRTrack& trk=*ptrk;
02626 goodBestTrkCounter++;
02627
02628
02629 if (!trk.fit.pass) continue;
02630 goodTrkPassCounter++;
02631
02632
02634
02635 lib.reco.GetEvtEnergy(nu);
02636 goodRecoEnCounter++;
02637
02638 if (!cuts.IsGoodUVVtx(nu)){
02639 continue;
02640 }
02641 goodUVDiffCounter++;
02642
02643
02645
02646 Bool_t goodSigmaQP=cuts.IsGoodSigmaQP_QP(nu);
02647 Bool_t goodFitProb=cuts.IsGoodFitProb(nu);
02648 Bool_t goodChi2=cuts.IsGoodFitChi2PerNdof(nu);
02649 Bool_t goodDpID=cuts.IsGoodDpID(nu);
02650
02651
02652 deltaTs.clear();
02653 reco.GetEvtDeltaTs(ntp,deltaTs);
02654 evtsPerSlc.clear();
02655 reco.GetEvtsPerSlc(ntp,evtsPerSlc);
02656 Bool_t goodDeltaT=deltaTs[ievt]>=50*Munits::ns;
02657 Bool_t goodEvtsPerSlc=evtsPerSlc[evt.slc]==1;
02658
02659
02660 if (nu.sigqp_qp<0) continue;
02661
02662
02663 hSigmaqp_qp->Fill(nu.sigqp_qp);
02664 hFitProb->Fill(nu.prob);
02665 hChi2->Fill(nu.chi2PerNdof);
02666 hDpID->Fill(nu.dpID);
02667 hDeltaT->Fill(deltaTs[ievt]/Munits::ns);
02668 hEvtsPerSlc->Fill(evtsPerSlc[evt.slc]);
02669
02671
02673 if (true &&
02674 goodFitProb &&
02675 goodChi2 &&
02676 goodDpID &&
02677 goodDeltaT &&
02678 goodEvtsPerSlc) {
02679 hSigmaqp_qpN_1->Fill(nu.sigqp_qp);
02680 }
02681
02682 if (goodSigmaQP &&
02683 true &&
02684 goodChi2 &&
02685 goodDpID &&
02686 goodDeltaT &&
02687 goodEvtsPerSlc){
02688 hFitProbN_1->Fill(nu.prob);
02689 }
02690
02691 if (goodSigmaQP &&
02692 goodFitProb &&
02693 true &&
02694 goodDpID &&
02695 goodDeltaT &&
02696 goodEvtsPerSlc){
02697 hChi2N_1->Fill(nu.chi2PerNdof);
02698 }
02699
02700 if (goodSigmaQP &&
02701 goodFitProb &&
02702 goodChi2 &&
02703 true &&
02704 goodDeltaT &&
02705 goodEvtsPerSlc){
02706 hDpIDN_1->Fill(nu.dpID);
02707 }
02708
02709 if (goodSigmaQP &&
02710 goodFitProb &&
02711 goodChi2 &&
02712 goodDpID &&
02713 true &&
02714 goodEvtsPerSlc){
02715 hDeltaTN_1->Fill(deltaTs[ievt]/Munits::ns);
02716 }
02717
02718 if (goodSigmaQP &&
02719 goodFitProb &&
02720 goodChi2 &&
02721 goodDpID &&
02722 goodDeltaT &&
02723 true){
02724 hEvtsPerSlcN_1->Fill(evtsPerSlc[evt.slc]);
02725 }
02726
02728
02730 if (goodSigmaQP &&
02731 goodFitProb &&
02732 goodChi2 &&
02733 goodDpID &&
02734 goodDeltaT &&
02735 goodEvtsPerSlc){
02736 hSigmaqp_qpN->Fill(nu.sigqp_qp);
02737 hFitProbN->Fill(nu.prob);
02738 hChi2N->Fill(nu.chi2PerNdof);
02739 hDpIDN->Fill(nu.dpID);
02740 hDeltaTN->Fill(deltaTs[ievt]/Munits::ns);
02741 hEvtsPerSlcN->Fill(evtsPerSlc[evt.slc]);
02742 }
02743 }
02744 }
02745
02749
02750 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
02751
02752
02753 MSG("NuAnalysis",Msg::kInfo)
02754 <<" ** Finished N_1 method **"<<endl;
02755 }
02756
02757
02758
02759 void NuAnalysis::OldDetermineBeamType(NuConfig& config)
02760 {
02761
02762
02763
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786
02787 if (config.beamTypeSntp==0) {
02788 config.targetPos=-1;
02789 config.sTargetPos="UknownTargetPos";
02790 }
02791 else if (config.beamTypeSntp==1) {
02792 config.targetPos=10;
02793 config.sTargetPos="010";
02794 }
02795 else if (config.beamTypeSntp==2) {
02796 config.targetPos=100;
02797 config.sTargetPos="100";
02798 }
02799 else if (config.beamTypeSntp==3) {
02800 config.targetPos=250;
02801 config.sTargetPos="250";
02802 }
02803 else if (config.beamTypeSntp==4) {
02804 config.targetPos=100;
02805 config.sTargetPos="100";
02806 }
02807 else if (config.beamTypeSntp==5) {
02808 config.targetPos=250;
02809 config.sTargetPos="250";
02810 }
02811 else {
02812 MAXMSG("NuAnalysis",Msg::kError,10)
02813 <<"Ahhh, don't recognise beamTypeSntp="
02814 <<config.beamTypeSntp<<endl;
02815 assert(false);
02816 }
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829 if (fabs(config.hornCurrent)>200-10-2){
02830 config.sHornCurrent="200";
02831 }
02832 else if (fabs(config.hornCurrent)>185-10-2) {
02833 config.sHornCurrent="185";
02834 }
02835 else if (fabs(config.hornCurrent)>170-10-2) {
02836 config.sHornCurrent="170";
02837 }
02838 else if (fabs(config.hornCurrent)>=0) {
02839 config.sHornCurrent="000";
02840 }
02841 else {
02842 MAXMSG("NuAnalysis",Msg::kError,10)
02843 <<"Ahh, no case for hornCurrent="<<config.hornCurrent<<endl;
02844 }
02845
02846
02847 config.sBeamType="L"+config.sTargetPos+config.sHornCurrent;
02848
02849
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878 this->DetermineBeamType(config);
02879 }
02880
02881
02882
02883 void NuAnalysis::DetermineBeamType(NuConfig& config) const
02884 {
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926 if (config.sTargetPos=="000") {
02927 MAXMSG("NuAnalysis",Msg::kInfo,5)
02928 <<"Identified case for sTargetPos="<<config.sTargetPos
02929 <<", now searching for horn current case..."<<endl;
02930 if (config.sHornCurrent=="200") {
02931 config.beamType=BeamType::kL000z200i;
02932 }
02933 else {
02934 MAXMSG("NuAnalysis",Msg::kError,5)
02935 <<"Ahhh, no horn current case for hornCurrent="
02936 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
02937 <<endl;
02938 config.beamType=BeamType::kUnknown;
02939 assert(false);
02940 }
02941 }
02942 else if (config.sTargetPos=="010") {
02943 MAXMSG("NuAnalysis",Msg::kInfo,5)
02944 <<"Identified case for sTargetPos="<<config.sTargetPos
02945 <<", now searching for horn current..."<<endl;
02946 if (config.sHornCurrent=="185") {
02947 config.beamType=BeamType::kL010z185i;
02948 }
02949 else if (config.sHornCurrent=="200") {
02950 config.beamType=BeamType::kL010z200i;
02951 }
02952 else if (config.sHornCurrent=="170") {
02953 config.beamType=BeamType::kL010z170i;
02954 }
02955 else if (config.sHornCurrent=="000") {
02956 config.beamType=BeamType::kL010z000i;
02957 }
02958 else {
02959 MAXMSG("NuAnalysis",Msg::kError,5)
02960 <<"Ahhh, no horn current case for hornCurrent="
02961 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
02962 <<endl;
02963 config.beamType=BeamType::kUnknown;
02964 assert(false);
02965 }
02966 }
02967 else if (config.sTargetPos=="50") {
02968 MAXMSG("NuAnalysis",Msg::kInfo,5)
02969 <<"Identified case for sTargetPos="<<config.sTargetPos
02970 <<", now searching for horn current..."<<endl;
02971 if (config.sHornCurrent=="200") {
02972 config.beamType=BeamType::kL050z200i;
02973 }
02974 else {
02975 MAXMSG("NuAnalysis",Msg::kError,5)
02976 <<"Ahhh, no horn current case for hornCurrent="
02977 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
02978 <<endl;
02979 config.beamType=BeamType::kUnknown;
02980 assert(false);
02981 }
02982 }
02983 else if (config.sTargetPos=="100") {
02984 MAXMSG("NuAnalysis",Msg::kInfo,5)
02985 <<"Identified case for sTargetPos="<<config.sTargetPos
02986 <<", now searching for horn current..."<<endl;
02987 if (config.sHornCurrent=="200") {
02988 config.beamType=BeamType::kL100z200i;
02989 }
02990 else {
02991 MAXMSG("NuAnalysis",Msg::kError,5)
02992 <<"Ahhh, no horn current case for hornCurrent="
02993 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
02994 <<endl;
02995 config.beamType=BeamType::kUnknown;
02996 assert(false);
02997 }
02998 }
02999 else if (config.sTargetPos=="150") {
03000 MAXMSG("NuAnalysis",Msg::kInfo,5)
03001 <<"Identified case for sTargetPos="<<config.sTargetPos
03002 <<", now searching for horn current..."<<endl;
03003 if (config.sHornCurrent=="200") {
03004 config.beamType=BeamType::kL150z200i;
03005 }
03006 else {
03007 MAXMSG("NuAnalysis",Msg::kError,5)
03008 <<"Ahhh, no horn current case for hornCurrent="
03009 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
03010 <<endl;
03011 config.beamType=BeamType::kUnknown;
03012 assert(false);
03013 }
03014 }
03015 else if (config.sTargetPos=="200") {
03016 MAXMSG("NuAnalysis",Msg::kInfo,5)
03017 <<"Identified case for sTargetPos="<<config.sTargetPos
03018 <<", now searching for horn current..."<<endl;
03019 if (config.sHornCurrent=="200") {
03020 config.beamType=BeamType::kL200z200i;
03021 }
03022 else {
03023 MAXMSG("NuAnalysis",Msg::kError,5)
03024 <<"Ahhh, no horn current case for hornCurrent="
03025 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
03026 <<endl;
03027 config.beamType=BeamType::kUnknown;
03028 assert(false);
03029 }
03030 }
03031 else if (config.sTargetPos=="250") {
03032 MAXMSG("NuAnalysis",Msg::kInfo,5)
03033 <<"Identified case for sTargetPos="<<config.sTargetPos
03034 <<", now searching for horn current..."<<endl;
03035 if (config.sHornCurrent=="200") {
03036 config.beamType=BeamType::kL250z200i;
03037 }
03038 else {
03039 MAXMSG("NuAnalysis",Msg::kError,5)
03040 <<"Ahhh, no horn current case for hornCurrent="
03041 <<config.sHornCurrent<<" for sTargetPos="<<config.sTargetPos
03042 <<endl;
03043 config.beamType=BeamType::kUnknown;
03044 assert(false);
03045 }
03046 }
03047 else {
03048 config.beamType=BeamType::kUnknown;
03049 MAXMSG("NuAnalysis",Msg::kError,10)
03050 <<"Ahhhhh, target position="<<config.sTargetPos
03051 <<" case not known..."<<endl;
03052
03053 }
03054 }
03055
03056
03057
03058 void NuAnalysis::ExtractConfig(NuConfig& config)
03059 {
03061
03062 this->ExtractConfig(this->GetNtpStRecord(),
03063 this->GetNtpBDLiteRecord(),config);
03064 }
03065
03066
03067
03068 Bool_t NuAnalysis::ExtractConfig(const NtpStRecord* pntp,
03069 const NtpBDLiteRecord* pntpBD,
03070 NuConfig& config) const
03071 {
03072
03073 const NtpStRecord& ntp=*pntp;
03074
03075 const RecCandHeader& rec=ntp.GetHeader();
03076 MAXMSG("NuAnalysis",Msg::kInfo,5)
03077 <<"ExtractConfig: run="<<rec.GetRun()
03078 <<", subrun="<<rec.GetSubRun()
03079 <<", detector="<<ntp.GetVldContext()->GetDetector()
03080 <<", simFlag="<<ntp.GetVldContext()->GetSimFlag()
03081 <<endl;
03082
03083
03084 config.run=rec.GetRun();
03085 config.subRun=rec.GetSubRun();
03086 config.snarl=rec.GetSnarl();
03087 config.timeSec=ntp.GetVldContext()->GetTimeStamp().GetSec();
03088 config.timeNanoSec=ntp.GetVldContext()->GetTimeStamp().GetNanoSec();
03089 config.detector=ntp.GetVldContext()->GetDetector();
03090 config.simFlag=ntp.GetVldContext()->GetSimFlag();
03091 config.trigSrc=rec.GetTrigSrc();
03092
03093
03094 MAXMSG("NuAnalysis",Msg::kInfo,50)
03095 <<"Using ntp.GetTitle()="<<ntp.GetTitle()<<endl
03096 <<"Using ntp.mchdr.geninfo.codename="
03097 <<ntp.mchdr.geninfo.codename<<endl;
03098
03099
03100 config.releaseType=ReleaseType::MakeReleaseType
03101 (ntp.GetTitle(),ntp.mchdr.geninfo.codename);
03102
03103 config.mcVersion=ReleaseType::GetMCInfo(config.releaseType);
03104
03105 config.recoVersion=ReleaseType::GetRecoInfo(config.releaseType);
03106
03107 MAXMSG("NuAnalysis",Msg::kInfo,50)
03108 <<"Extracted:"<<endl
03109 <<" ReleaseType="<<config.releaseType<<" "
03110 <<ReleaseType::AsString(config.releaseType)<<endl
03111 <<" mcVersion="<<config.mcVersion<<" "
03112 <<ReleaseType::AsString(config.mcVersion)<<endl
03113 <<" recoVersion="<<config.recoVersion<<" "
03114 <<ReleaseType::AsString(config.recoVersion)<<endl;
03115
03117
03119 if (ntp.GetVldContext()->GetSimFlag()==SimFlag::kData &&
03120 (!pntpBD || config.overrideBeamDataConfigExtraction)) {
03121 if (!pntpBD) {
03122 MSG("NuAnalysis",Msg::kInfo)
03123 <<"ExtractConfig: no NtpBDLiteRecord available"
03124 <<", making up settings..."<<endl;
03125 }
03126 else if (config.overrideBeamDataConfigExtraction) {
03127 MSG("NuAnalysis",Msg::kInfo)
03128 <<"ExtractConfig: overriding use of Beam Data"
03129 <<" to extract config, making up settings..."<<endl;
03130 }
03131 config.beamType=BeamType::kL010z185i;
03132 config.hornCurrent=-185;
03133 string btAsString=BeamType::AsString
03134 (static_cast<BeamType::BeamType_t>(config.beamType));
03135 config.sTargetPos=btAsString.substr(1,3);
03136 config.targetPos=-1.*atoi(config.sTargetPos.c_str());
03137 config.sHornCurrent=btAsString.substr(5,7);
03138
03139 config.sBeamType="L"+config.sTargetPos+config.sHornCurrent;
03140
03141 config.beamTypeSntp=1;
03142 }
03144
03146 else if (ntp.GetVldContext()->GetSimFlag()==SimFlag::kData){
03147 MSG("NuAnalysis",Msg::kInfo)
03148 <<"ExtractConfig: running on DATA..."<<endl;
03149 const NtpBDLiteRecord& ntpBD=*pntpBD;
03150
03151
03152
03153 BeamMonSpill spill;
03154 BMSpillAna bmsa;
03155
03156
03157
03158
03159 bmsa.UseDatabaseCuts();
03160
03161
03162 bmsa.SetSpill(ntpBD,spill);
03163
03164 MAXMSG("NuAnalysis",Msg::kInfo,10)
03165 <<"ExtractConfig: after SetSpill & UseDatabaseCuts:"<<endl;
03166 bmsa.Print();
03167 cout<<endl;
03168
03169
03170 if (!bmsa.SelectSpill()) {
03171 MSG("NuBeam",Msg::kWarning)
03172 <<"ExtractConfig: BAD SPILL"
03173 <<", tortgt="<<ntpBD.tortgt
03174 <<", timeDiff="<<ntpBD.GetHeader().GetTimeDiffStreamSpill()
03175 <<", hornCur="<<ntpBD.horncur
03176 <<", bt="<<BeamType::AsString(static_cast<BeamType::BeamType_t>
03177 (spill.BeamType()))
03178 <<" ("<<spill.BeamType()<<")"<<endl;
03179 MSG("NuBeam",Msg::kWarning)
03180 <<endl<<endl
03181 <<"Occasionally the first spill is bad so have to skip to next"
03182 <<" snarl in order to ExtractConfig... return..."
03183 <<endl<<endl<<endl;
03184 return false;
03185 }
03186
03187
03188
03189
03190 config.beamType=spill.BeamType();
03191
03192
03193 config.hornCurrent=ntpBD.horncur;
03194
03195 string btAsString=BeamType::AsString
03196 (static_cast<BeamType::BeamType_t>(config.beamType));
03197
03198 config.sTargetPos=btAsString.substr(1,3);
03199 config.targetPos=-1.*atoi(config.sTargetPos.c_str());
03200 config.sHornCurrent=btAsString.substr(5,7);
03201
03202
03203 config.sBeamType="L"+config.sTargetPos+config.sHornCurrent;
03204
03205
03206
03207
03208 config.beamTypeSntp=spill.GetStatusBits().beam_type;
03209 }
03211
03213 else if (ntp.GetVldContext()->GetSimFlag()==SimFlag::kMC){
03214 MSG("NuAnalysis",Msg::kInfo)
03215 <<"ExtractConfig: running on MC..."<<endl;
03216
03217 string fileName="n13021160_0000_L010185.sntp.R1_18_2.root";
03218
03219 if (this->GetNtpStRecord()) {
03220 fileName=this->GetFirstFileName();
03221 }
03222 else {
03223 fileName=this->GetInputFileName();
03224 }
03225
03226 Int_t lastSlash=fileName.find_last_of("/");
03227 fileName.erase(0,lastSlash+1);
03228 MAXMSG("NuAnalysis",Msg::kInfo,5)
03229 <<"Extracting MC beam info from run filename="<<fileName<<endl;
03230
03231
03232
03233 UInt_t strPosStartPre=0;
03234 UInt_t strLengthPre=6;
03235 if (strPosStartPre+strLengthPre<=fileName.size()) {
03236 string preDiakon=fileName.substr(0,6);
03237 MAXMSG("NuAnalysis",Msg::kDebug,50)
03238 <<"Checking for preDiakon name, fileName.substr(0,6)="
03239 <<preDiakon<<endl;
03240 if (preDiakon=="Cedar-") {
03241 fileName.erase(0,6);
03242 }
03243 }
03244
03245
03246 UInt_t strPosStartTarget=16;
03247 UInt_t strLengthTarget=3;
03248 if (strPosStartTarget+strLengthTarget<=fileName.size()) {
03249 config.sTargetPos=fileName.substr(strPosStartTarget,
03250 strLengthTarget);
03251 }
03252 else {
03253 MSG("NuAnalysis",Msg::kError)
03254 <<"Assuming targetPos=010"<<endl;
03255 config.sTargetPos="010";
03256 }
03257
03258
03259 UInt_t strPosStartHornCur=19;
03260 UInt_t strLengthHornCur=3;
03261 if (strPosStartHornCur+strLengthHornCur<=fileName.size()) {
03262 config.sHornCurrent=fileName.substr(19,3);
03263 }
03264 else {
03265 MSG("NuAnalysis",Msg::kError)
03266 <<"Assuming hornCurrent=185"<<endl;
03267 config.sHornCurrent="185";
03268 }
03269
03270
03271 config.sBeamType="L"+config.sTargetPos+config.sHornCurrent;
03272
03273
03274
03275 MAXMSG("NuAnalysis",Msg::kInfo,5)
03276 <<"Extracted strings: sTargetPos="<<config.sTargetPos
03277 <<", sHornCurrent="<<config.sHornCurrent
03278 <<", mcVersion="<<config.mcVersion
03279 <<endl;
03280
03281
03282 config.hornCurrent=-1.*atoi(config.sHornCurrent.c_str());
03283 config.targetPos=-1.*atoi(config.sTargetPos.c_str());
03284
03285
03286 this->DetermineBeamType(config);
03287 }
03288
03289 MSG("NuAnalysis",Msg::kInfo)
03290 <<"ExtractConfig: final config is beamType="
03291 <<BeamType::AsString(static_cast<BeamType::BeamType_t>
03292 (config.beamType))
03293 <<" ("<<config.beamType<<")"
03294 <<", btSntp="<<config.beamTypeSntp
03295 <<endl
03296 <<" sHornCurrent="<<config.sHornCurrent
03297 <<" ("<<config.hornCurrent<<")"
03298 <<", sTargetPos="<<config.sTargetPos
03299 <<" ("<<config.targetPos<<")"
03300 <<", sBeamType="<<config.sBeamType
03301 <<endl;
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311
03312
03313
03314
03315
03317
03318
03319 return true;
03320 }
03321
03322
03323
03324 void NuAnalysis::NuMuBarAppearance()
03325 {
03326 MSG("NuAnalysis",Msg::kInfo)
03327 <<" ** Running NuMuBarAppearance method... **"<<endl;
03328
03329 NuConfig config;
03330 this->ExtractConfig(config);
03331
03332
03333 string sFilePrefix="NMBApp"+config.sBeamType;
03334 fOutFile=this->OpenFile(config,sFilePrefix.c_str());
03335
03336 TH1F* hRecoEnNMBFull=new TH1F("hRecoEnNMBFull","hRecoEnNMBFull",
03337 4*352,-32,320);
03338 hRecoEnNMBFull->GetXaxis()->SetTitle("True Energy (GeV)");
03339 hRecoEnNMBFull->GetXaxis()->CenterTitle();
03340 hRecoEnNMBFull->GetYaxis()->SetTitle("");
03341 hRecoEnNMBFull->GetYaxis()->CenterTitle();
03342
03343 hRecoEnNMBFull->Sumw2();
03344
03345 TH1F* hRecoEnNMFull=new TH1F("hRecoEnNMFull","hRecoEnNMFull",
03346 4*352,-32,320);
03347 hRecoEnNMFull->GetXaxis()->SetTitle("True Energy (GeV)");
03348 hRecoEnNMFull->GetXaxis()->CenterTitle();
03349 hRecoEnNMFull->GetYaxis()->SetTitle("");
03350 hRecoEnNMFull->GetYaxis()->CenterTitle();
03351
03352 hRecoEnNMFull->Sumw2();
03353
03354 TH1F* hRecoEnNMBOsc=new TH1F("hRecoEnNMBOsc","hRecoEnNMBOsc",
03355 4*352,-32,320);
03356 hRecoEnNMBOsc->GetXaxis()->SetTitle("True Energy (GeV)");
03357 hRecoEnNMBOsc->GetXaxis()->CenterTitle();
03358 hRecoEnNMBOsc->GetYaxis()->SetTitle("");
03359 hRecoEnNMBOsc->GetYaxis()->CenterTitle();
03360
03361 hRecoEnNMBOsc->Sumw2();
03362
03363 TH1F* hRecoEnNMOsc=new TH1F("hRecoEnNMOsc","hRecoEnNMOsc",
03364 4*352,-32,320);
03365 hRecoEnNMOsc->GetXaxis()->SetTitle("True Energy (GeV)");
03366 hRecoEnNMOsc->GetXaxis()->CenterTitle();
03367 hRecoEnNMOsc->GetYaxis()->SetTitle("");
03368 hRecoEnNMOsc->GetYaxis()->CenterTitle();
03369
03370 hRecoEnNMOsc->Sumw2();
03371
03372
03373 Int_t entries=static_cast<Int_t>(this->GetEntries());
03374 Int_t startTimeSecs=2000000000;
03375 Int_t endTimeSecs=1;
03376 vector<Double_t> vPot;
03377 vPot.reserve(entries);
03378 vector<Double_t> vTime;
03379 vTime.reserve(entries);
03380 vector<Double_t> vPotBad;
03381 vPotBad.reserve(entries/10);
03382 vector<Double_t> vTimeBad;
03383 vTimeBad.reserve(entries/10);
03384
03385 vector<Double_t> vTimeNuMuEvt;
03386 vTimeNuMuEvt.reserve(entries/5);
03387 vector<Double_t> vTimeNuMuBarEvt;
03388 vTimeNuMuBarEvt.reserve(entries/5);
03389
03390
03391 map<Int_t,Double_t> deltaTs;
03392 map<Int_t,Int_t> evtsPerSlc;
03393
03394
03395 Int_t goodSpillCounter=0;
03396 Int_t badSpillCounter=0;
03397 Float_t totalPot=0;
03398 Float_t totalPotBad=0;
03399 Int_t evtCounter=0;
03400 Int_t evtNotlitime=0;
03401 Int_t evtNotIsLI=0;
03402 Int_t evtInFidVolCounter=0;
03403 Int_t evtWithTrkCounter=0;
03404 Int_t goodBestTrkCounter=0;
03405 Int_t goodTrkPassCounter=0;
03406 Int_t goodRecoEnCounter=0;
03407 Int_t goodUVDiffCounter=0;
03408 Int_t trkInFidVolCounter=0;
03409 Int_t goodFitSigQPCounter=0;
03410 Int_t goodFitProbCounter=0;
03411 Int_t goodFitChi2Counter=0;
03412 Int_t goodDeltaTCounter=0;
03413 Int_t goodEvtsPerSlcCounter=0;
03414 Int_t goodPIDCounter=0;
03415 Int_t goodBeyondTrkEndCounter=0;
03416 Int_t goodShwFractCounter=0;
03417 Int_t nuNQCounter=0;
03418 Int_t nuPQCounter=0;
03419
03420
03421
03422
03423
03424
03425 const NuCuts cuts;
03426 const NuGeneral general;
03427 const NuPlots plots;
03428 const NuReco reco;
03429 const NuExtraction ext;
03430 const NuBeam beam;
03431
03432
03433
03434 NuMadAnalysis& mad=*new NuMadAnalysis();
03435 MadDpID& dp=*new MadDpID();
03436
03437 BeamType::BeamType_t beamType=static_cast
03438 <BeamType::BeamType_t>(config.beamType);
03439
03443
03444 this->InitialiseLoopVariables();
03445
03446
03447 for(Int_t entry=0;entry<this->GetEntries();entry++){
03448
03449 this->SetLoopVariables(entry);
03450
03451
03452 const NtpStRecord& ntp=(*this->GetNtpStRecord());
03453
03454 const RecCandHeader& rec=ntp.GetHeader();
03455 MAXMSG("NuAnalysis",Msg::kInfo,5)
03456 <<"Found: run="<<rec.GetRun()
03457 <<", subrun="<<rec.GetSubRun()
03458 <<", detector="<<rec.GetVldContext().GetDetector()
03459 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
03460 <<endl;
03461
03462
03463 plots.FillGeneralHistos(ntp);
03464
03465
03466 NuEvent nuSnarl;
03467
03468
03469 ext.ExtractGeneralInfo(ntp,nuSnarl);
03470
03471 if (!beam.IsGoodSpillAndFillPot
03472 (this->GetNtpBDLiteRecord(),config,nuSnarl)) {
03473 badSpillCounter++;
03474 totalPotBad+=nuSnarl.pot;
03475 continue;
03476 }
03477 goodSpillCounter++;
03478 totalPot+=nuSnarl.pot;
03479
03480 dp.ChoosePDFs(rec.GetVldContext().GetDetector(),beamType);
03481 mad.SetEntry(&ntp);
03482
03483 set<Int_t> setNuMuBarInTrueFidCC;
03484
03485 VldTimeStamp vldts;
03486
03487 VldContext vc(rec.GetVldContext().GetDetector(),
03488 rec.GetVldContext().GetSimFlag(),vldts);
03489
03490 UgliGeomHandle ugh(vc);
03491
03492 TClonesArray& mcTca=(*ntp.mc);
03493 Int_t numInt=mcTca.GetEntries();
03494 MAXMSG("NuAnalysis",Msg::kInfo,10)
03495 <<"Number of entries in NtpMCTruth="
03496 <<mcTca.GetEntries()<<endl;
03497
03498
03499 NuEvent nu;
03500
03501
03502 ext.ExtractGeneralInfo(ntp,nu);
03503
03504
03506
03508 for (Int_t i=0;i<numInt;i++){
03509 const NtpMCTruth& mc=*(dynamic_cast<NtpMCTruth*>(mcTca[i]));
03510
03511 TVector3 xyz(mc.vtxx,mc.vtxy,mc.vtxz);
03512
03513 TVector3 uvz=ugh.xyz2uvz(xyz);
03514
03515
03516
03517 if (mc.iaction!=1) continue;
03518 if (abs(mc.inu)!=14) continue;
03519
03520 Bool_t isInFidVol=cuts.IsInFidVol(mc.vtxx,mc.vtxy,mc.vtxz,
03521 uvz.X(),uvz.Y(),0,0,
03522 nu.detector,nu.anaVersion,
03523 nu.releaseType,nu.simFlag);
03524
03525 if (!isInFidVol) continue;
03526 setNuMuBarInTrueFidCC.insert(mc.index);
03527
03528 Double_t oscWeight=general.OscWeight(3e-3,1,fabs(mc.p4neu[3]));
03529
03530 MAXMSG("NuAnalysis",Msg::kInfo,200)
03531 <<"inu="<<mc.inu<<", iaction="<<mc.iaction
03532 <<", E="<<mc.p4neu[3]<<", oscWeight="<<oscWeight<<endl;
03533
03534 if (mc.inu==-14) {
03535 hRecoEnNMBFull->Fill(fabs(mc.p4neu[3]));
03536 hRecoEnNMBOsc->Fill(fabs(mc.p4neu[3]),oscWeight);
03537 }
03538 else if (mc.inu==14) {
03539 hRecoEnNMFull->Fill(fabs(mc.p4neu[3]));
03540 hRecoEnNMOsc->Fill(fabs(mc.p4neu[3]),oscWeight);
03541 }
03542 else {
03543 cout<<"Ahhh, bad inu"<<endl;
03544 }
03545
03546 }
03547
03548 MAXMSG("NuAnalysis",Msg::kInfo,20)
03549 <<"Number of NuMuBar fiducial volume interactions in spill="
03550 <<setNuMuBarInTrueFidCC.size()<<endl;
03551 }
03552
03556
03557 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
03558
03559
03560 Int_t nf=this->GetNumFilesAddedToChain();
03561 TH1F* hRun=dynamic_cast<TH1F*>(gROOT->FindObject("hRun"));
03562 Float_t nr=-1;
03563 if (hRun) nr=hRun->GetEntries();
03564 TH1F* hNumFiles=new TH1F("hNumFiles","hNumFiles",15,-5,10);
03565 hNumFiles->Fill(1,nf);
03566 MSG("NuAnalysis",Msg::kInfo)
03567 <<"Files added to chain="<<hNumFiles->GetMaximum()
03568 <<", runs found="<<nr<<endl;
03569
03570 MSG("NuAnalysis",Msg::kInfo)
03571 <<"Found start time and end time: "
03572 <<startTimeSecs<<" -> "<<endTimeSecs<<endl;
03573 general.EpochTo1995(endTimeSecs);
03574 general.EpochTo1995(startTimeSecs);
03575
03576 MSG("NuAnalysis",Msg::kInfo)<<"Filling vs. Time plots..."<<endl;
03577 TH1F* hSpillsVsTime=new TH1F("hSpillsVsTime","hSpillsVsTime",
03578 100,startTimeSecs,endTimeSecs);
03579 general.TH1FFill(hSpillsVsTime,vTime);
03580 general.SetGraphAxis(hSpillsVsTime->GetXaxis(),
03581 startTimeSecs,endTimeSecs);
03582 hSpillsVsTime->GetXaxis()->CenterTitle();
03583 hSpillsVsTime->GetYaxis()->SetTitle("Spills");
03584 hSpillsVsTime->GetYaxis()->CenterTitle();
03585
03586 TH1F* hSpillsBadVsTime=new TH1F("hSpillsBadVsTime","hSpillsBadVsTime",
03587 100,startTimeSecs,endTimeSecs);
03588 general.TH1FFill(hSpillsBadVsTime,vTimeBad);
03589 general.SetGraphAxis(hSpillsBadVsTime->GetXaxis(),
03590 startTimeSecs,endTimeSecs);
03591 hSpillsBadVsTime->GetXaxis()->CenterTitle();
03592 hSpillsBadVsTime->GetYaxis()->SetTitle("Spills");
03593 hSpillsBadVsTime->GetYaxis()->CenterTitle();
03594
03595 TH1F* hNuMuVsTime=new TH1F("hNuMuVsTime","hNuMuVsTime",
03596 100,startTimeSecs,endTimeSecs);
03597 general.TH1FFill(hNuMuVsTime,vTimeNuMuEvt);
03598 general.SetGraphAxis(hNuMuVsTime->GetXaxis(),
03599 startTimeSecs,endTimeSecs);
03600 hNuMuVsTime->GetXaxis()->CenterTitle();
03601 hNuMuVsTime->GetYaxis()->SetTitle("#nu_#mu events");
03602 hNuMuVsTime->GetYaxis()->CenterTitle();
03603
03604 TH1F* hNuMuBarVsTime=new TH1F("hNuMuBarVsTime","hNuMuBarVsTime",
03605 100,startTimeSecs,endTimeSecs);
03606 general.TH1FFill(hNuMuBarVsTime,vTimeNuMuBarEvt);
03607 general.SetGraphAxis(hNuMuBarVsTime->GetXaxis(),
03608 startTimeSecs,endTimeSecs);
03609 hNuMuBarVsTime->GetXaxis()->CenterTitle();
03610 hNuMuBarVsTime->GetYaxis()->SetTitle("#bar{#nu_#mu} events");
03611 hNuMuBarVsTime->GetYaxis()->CenterTitle();
03612
03613 TH1F* hPotVsTime=new TH1F("hPotVsTime","hPotVsTime",
03614 100,startTimeSecs,endTimeSecs);
03615 general.TH1FFill(hPotVsTime,vTime,vPot);
03616 general.SetGraphAxis(hPotVsTime->GetXaxis(),
03617 startTimeSecs,endTimeSecs);
03618 hPotVsTime->GetXaxis()->CenterTitle();
03619 hPotVsTime->GetYaxis()->SetTitle("POT");
03620 hPotVsTime->GetYaxis()->CenterTitle();
03621
03622 TH1F* hPotBadVsTime=new TH1F("hPotBadVsTime","hPotBadVsTime",
03623 100,startTimeSecs,endTimeSecs);
03624 general.TH1FFill(hPotBadVsTime,vTimeBad,vPotBad);
03625 general.SetGraphAxis(hPotBadVsTime->GetXaxis(),
03626 startTimeSecs,endTimeSecs);
03627 hPotBadVsTime->GetXaxis()->CenterTitle();
03628 hPotBadVsTime->GetYaxis()->SetTitle("POT");
03629 hPotBadVsTime->GetYaxis()->CenterTitle();
03630
03631 Float_t totalPotHist=0;
03632 Float_t totalPotBadHist=0;
03633 cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
03634 Float_t pcPotRejected=0;
03635 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
03636
03637
03638 cout<<"Number of spills ="<<this->GetEntries()<<endl
03639 <<"Number of good spills="<<goodSpillCounter<<endl
03640 <<"Number of bad spills ="<<badSpillCounter<<endl
03641 <<"Sum of good+bad spill="<<badSpillCounter+goodSpillCounter<<endl
03642 <<endl
03643 <<"Total POT good ="<<totalPot*1e12<<endl
03644 <<"Total POT bad ="<<totalPotBad*1e12<<endl
03645 <<"Total POT good (hist)="<<totalPotHist*1e12<<endl
03646 <<"Total POT bad (hist)="<<totalPotBadHist*1e12
03647 <<" ("<<pcPotRejected<<"%)"<<endl;
03648
03649 MSG("NuAnalysis",Msg::kInfo)
03650 <<endl<<"*************************************"<<endl
03651 <<"Total Events ="<<evtCounter<<endl
03652 <<"Evt not litime ="<<evtNotlitime<<endl
03653 <<"Evt not IsLI ="<<evtNotIsLI<<endl
03654 <<"Evts in Fid ="<<evtInFidVolCounter<<endl
03655 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
03656 <<"Good Best Trk ="<<goodBestTrkCounter<<endl
03657 <<"Track Pass ="<<goodTrkPassCounter<<endl
03658 <<"Good reco energy="<<goodRecoEnCounter<<endl
03659 <<"Good UVDiff ="<<goodUVDiffCounter<<endl
03660 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
03661 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
03662 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
03663 <<"Good Fit Chi2 ="<<goodFitChi2Counter<<endl
03664 <<"Good Delta T ="<<goodDeltaTCounter<<endl
03665 <<"Good EvtPerSlc ="<<goodEvtsPerSlcCounter<<endl
03666 <<"Good DP ID ="<<goodPIDCounter<<endl
03667 <<"Good beyond tEnd="<<goodBeyondTrkEndCounter<<endl
03668 <<"Good shw fract ="<<goodShwFractCounter<<endl
03669 <<"*************************************"<<endl;
03670
03671 Float_t pcPQ=0;
03672 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
03673 (nuPQCounter+nuNQCounter);
03674 MSG("NuAnalysis",Msg::kInfo)
03675 <<endl
03676 <<"Neutrinos with NQ="<<nuNQCounter<<endl
03677 <<"Neutrinos with PQ="<<nuPQCounter
03678 <<" ("<<pcPQ<<"%)"<<endl
03679 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
03680 <<endl;
03681
03682 MSG("NuAnalysis",Msg::kInfo)
03683 <<" ** Finished NuMuBarAppearance method **"<<endl;
03684 }
03685
03686
03687
03688 void NuAnalysis::Efficiencies()
03689 {
03690 MSG("NuAnalysis",Msg::kInfo)
03691 <<" ** Running Efficiencies method... **"<<endl;
03692
03693 NuConfig config;
03694 this->ExtractConfig(config);
03695
03696
03697 string sFilePrefix="NMBEff"+config.sBeamType;
03698 fOutFile=this->OpenFile(config,sFilePrefix.c_str());
03699
03700 TH1F* hRecoEnNMBFull=new TH1F("hRecoEnNMBFull","hRecoEnNMBFull",
03701 4*352,-32,320);
03702 hRecoEnNMBFull->GetXaxis()->SetTitle("True Energy (GeV)");
03703 hRecoEnNMBFull->GetXaxis()->CenterTitle();
03704 hRecoEnNMBFull->GetYaxis()->SetTitle("");
03705 hRecoEnNMBFull->GetYaxis()->CenterTitle();
03706
03707
03708 TH1F* hRecoEnNMBReco=new TH1F("hRecoEnNMBReco","hRecoEnNMBReco",
03709 4*352,-32,320);
03710 hRecoEnNMBReco->GetXaxis()->SetTitle("True Energy (GeV)");
03711 hRecoEnNMBReco->GetXaxis()->CenterTitle();
03712 hRecoEnNMBReco->GetYaxis()->SetTitle("");
03713 hRecoEnNMBReco->GetYaxis()->CenterTitle();
03714
03715
03716 TH1F* hRecoEnNMBFid=new TH1F("hRecoEnNMBFid","hRecoEnNMBFid",
03717 4*352,-32,320);
03718 hRecoEnNMBFid->GetXaxis()->SetTitle("True Energy (GeV)");
03719 hRecoEnNMBFid->GetXaxis()->CenterTitle();
03720 hRecoEnNMBFid->GetYaxis()->SetTitle("");
03721 hRecoEnNMBFid->GetYaxis()->CenterTitle();
03722
03723
03724 TH1F* hRecoEnNMBTrack=new TH1F("hRecoEnNMBTrack","hRecoEnNMBTrack",
03725 4*352,-32,320);
03726 hRecoEnNMBTrack->GetXaxis()->SetTitle("True Energy (GeV)");
03727 hRecoEnNMBTrack->GetXaxis()->CenterTitle();
03728 hRecoEnNMBTrack->GetYaxis()->SetTitle("");
03729 hRecoEnNMBTrack->GetYaxis()->CenterTitle();
03730
03731
03732 TH1F* hRecoEnNMBPass=new TH1F("hRecoEnNMBPass","hRecoEnNMBPass",
03733 4*352,-32,320);
03734 hRecoEnNMBPass->GetXaxis()->SetTitle("True Energy (GeV)");
03735 hRecoEnNMBPass->GetXaxis()->CenterTitle();
03736 hRecoEnNMBPass->GetYaxis()->SetTitle("");
03737 hRecoEnNMBPass->GetYaxis()->CenterTitle();
03738
03739
03740 TH1F* hRecoEnNMBQCut=new TH1F("hRecoEnNMBQCut","hRecoEnNMBQCut",
03741 4*352,-32,320);
03742 hRecoEnNMBQCut->GetXaxis()->SetTitle("True Energy (GeV)");
03743 hRecoEnNMBQCut->GetXaxis()->CenterTitle();
03744 hRecoEnNMBQCut->GetYaxis()->SetTitle("");
03745 hRecoEnNMBQCut->GetYaxis()->CenterTitle();
03746
03747
03748 TH1F* hRecoEnNMBSigQP=new TH1F("hRecoEnNMBSigQP","hRecoEnNMBSigQP",
03749 4*352,-32,320);
03750 hRecoEnNMBSigQP->GetXaxis()->SetTitle("True Energy (GeV)");
03751 hRecoEnNMBSigQP->GetXaxis()->CenterTitle();
03752 hRecoEnNMBSigQP->GetYaxis()->SetTitle("");
03753 hRecoEnNMBSigQP->GetYaxis()->CenterTitle();
03754
03755
03756 TH1F* hRecoEnNMBProb=new TH1F("hRecoEnNMBProb","hRecoEnNMBProb",
03757 4*352,-32,320);
03758 hRecoEnNMBProb->GetXaxis()->SetTitle("True Energy (GeV)");
03759 hRecoEnNMBProb->GetXaxis()->CenterTitle();
03760 hRecoEnNMBProb->GetYaxis()->SetTitle("");
03761 hRecoEnNMBProb->GetYaxis()->CenterTitle();
03762
03763
03764 TH1F* hRecoEnNMBProb05=new TH1F("hRecoEnNMBProb05","hRecoEnNMBProb05",
03765 4*352,-32,320);
03766 hRecoEnNMBProb05->GetXaxis()->SetTitle("True Energy (GeV)");
03767 hRecoEnNMBProb05->GetXaxis()->CenterTitle();
03768 hRecoEnNMBProb05->GetYaxis()->SetTitle("");
03769 hRecoEnNMBProb05->GetYaxis()->CenterTitle();
03770
03771
03772 TH1F* hRecoEnNMBDpID01=new TH1F("hRecoEnNMBDpID01","hRecoEnNMBDpID01",
03773 4*352,-32,320);
03774 hRecoEnNMBDpID01->GetXaxis()->SetTitle("True Energy (GeV)");
03775 hRecoEnNMBDpID01->GetXaxis()->CenterTitle();
03776 hRecoEnNMBDpID01->GetYaxis()->SetTitle("");
03777 hRecoEnNMBDpID01->GetYaxis()->CenterTitle();
03778
03779
03780 TH1F* hRecoEnNMBDpID04=new TH1F("hRecoEnNMBDpID04","hRecoEnNMBDpID04",
03781 4*352,-32,320);
03782 hRecoEnNMBDpID04->GetXaxis()->SetTitle("True Energy (GeV)");
03783 hRecoEnNMBDpID04->GetXaxis()->CenterTitle();
03784 hRecoEnNMBDpID04->GetYaxis()->SetTitle("");
03785 hRecoEnNMBDpID04->GetYaxis()->CenterTitle();
03786
03787
03788
03789
03790 TH1F* hRecoEnNMBLeakIn=new TH1F("hRecoEnNMBLeakIn","hRecoEnNMBLeakIn",
03791 4*352,-32,320);
03792 hRecoEnNMBLeakIn->GetXaxis()->SetTitle("True Energy (GeV)");
03793 hRecoEnNMBLeakIn->GetXaxis()->CenterTitle();
03794 hRecoEnNMBLeakIn->GetYaxis()->SetTitle("");
03795 hRecoEnNMBLeakIn->GetYaxis()->CenterTitle();
03796
03797
03798 TH1F* hRecoEnNMBLeakOut=new TH1F("hRecoEnNMBLeakOut",
03799 "hRecoEnNMBLeakOut",
03800 4*352,-32,320);
03801 hRecoEnNMBLeakOut->GetXaxis()->SetTitle("True Energy (GeV)");
03802 hRecoEnNMBLeakOut->GetXaxis()->CenterTitle();
03803 hRecoEnNMBLeakOut->GetYaxis()->SetTitle("");
03804 hRecoEnNMBLeakOut->GetYaxis()->CenterTitle();
03805
03806
03807
03808
03809 TH2F* hYvsXTrueLeakIn=new TH2F("hYvsXTrueLeakIn","hYvsXTrueLeakIn",
03810 100,-5,5,100,-5,5);
03811 hYvsXTrueLeakIn->SetTitle("Y vs X");
03812 hYvsXTrueLeakIn->GetXaxis()->SetTitle("X");
03813 hYvsXTrueLeakIn->GetXaxis()->CenterTitle();
03814 hYvsXTrueLeakIn->GetYaxis()->SetTitle("Y");
03815 hYvsXTrueLeakIn->GetYaxis()->CenterTitle();
03816
03817
03818 TH2F* hYvsXRecoLeakIn=new TH2F("hYvsXRecoLeakIn","hYvsXRecoLeakIn",
03819 100,-5,5,100,-5,5);
03820 hYvsXRecoLeakIn->SetTitle("Y vs X");
03821 hYvsXRecoLeakIn->GetXaxis()->SetTitle("X");
03822 hYvsXRecoLeakIn->GetXaxis()->CenterTitle();
03823 hYvsXRecoLeakIn->GetYaxis()->SetTitle("Y");
03824 hYvsXRecoLeakIn->GetYaxis()->CenterTitle();
03825
03826
03827 TH2F* hYvsXTrueLeakOut=new TH2F("hYvsXTrueLeakOut","hYvsXTrueLeakOut",
03828 100,-5,5,100,-5,5);
03829 hYvsXTrueLeakOut->SetTitle("Y vs X");
03830 hYvsXTrueLeakOut->GetXaxis()->SetTitle("X");
03831 hYvsXTrueLeakOut->GetXaxis()->CenterTitle();
03832 hYvsXTrueLeakOut->GetYaxis()->SetTitle("Y");
03833 hYvsXTrueLeakOut->GetYaxis()->CenterTitle();
03834
03835
03836 TH2F* hYvsXRecoLeakOut=new TH2F("hYvsXRecoLeakOut","hYvsXRecoLeakOut",
03837 100,-5,5,100,-5,5);
03838 hYvsXRecoLeakOut->SetTitle("Y vs X");
03839 hYvsXRecoLeakOut->GetXaxis()->SetTitle("X");
03840 hYvsXRecoLeakOut->GetXaxis()->CenterTitle();
03841 hYvsXRecoLeakOut->GetYaxis()->SetTitle("Y");
03842 hYvsXRecoLeakOut->GetYaxis()->CenterTitle();
03843
03844
03845
03846 TH1F* hRecoEnNMBTrkLeakIn=new TH1F("hRecoEnNMBTrkLeakIn",
03847 "hRecoEnNMBTrkLeakIn",
03848 4*352,-32,320);
03849 hRecoEnNMBTrkLeakIn->GetXaxis()->SetTitle("True Energy (GeV)");
03850 hRecoEnNMBTrkLeakIn->GetXaxis()->CenterTitle();
03851 hRecoEnNMBTrkLeakIn->GetYaxis()->SetTitle("");
03852 hRecoEnNMBTrkLeakIn->GetYaxis()->CenterTitle();
03853
03854
03855 TH1F* hRecoEnNMBTrkLeakOut=new TH1F("hRecoEnNMBTrkLeakOut",
03856 "hRecoEnNMBTrkLeakOut",
03857 4*352,-32,320);
03858 hRecoEnNMBTrkLeakOut->GetXaxis()->SetTitle("True Energy (GeV)");
03859 hRecoEnNMBTrkLeakOut->GetXaxis()->CenterTitle();
03860 hRecoEnNMBTrkLeakOut->GetYaxis()->SetTitle("");
03861 hRecoEnNMBTrkLeakOut->GetYaxis()->CenterTitle();
03862
03863
03864
03865
03866 TH2F* hYvsXTrueTrkLeakIn=new TH2F("hYvsXTrueTrkLeakIn",
03867 "hYvsXTrueTrkLeakIn",
03868 100,-5,5,100,-5,5);
03869 hYvsXTrueTrkLeakIn->SetTitle("Y vs X");
03870 hYvsXTrueTrkLeakIn->GetXaxis()->SetTitle("X");
03871 hYvsXTrueTrkLeakIn->GetXaxis()->CenterTitle();
03872 hYvsXTrueTrkLeakIn->GetYaxis()->SetTitle("Y");
03873 hYvsXTrueTrkLeakIn->GetYaxis()->CenterTitle();
03874
03875
03876 TH2F* hYvsXRecoTrkLeakIn=new TH2F("hYvsXRecoTrkLeakIn",
03877 "hYvsXRecoTrkLeakIn",
03878 100,-5,5,100,-5,5);
03879 hYvsXRecoTrkLeakIn->SetTitle("Y vs X");
03880 hYvsXRecoTrkLeakIn->GetXaxis()->SetTitle("X");
03881 hYvsXRecoTrkLeakIn->GetXaxis()->CenterTitle();
03882 hYvsXRecoTrkLeakIn->GetYaxis()->SetTitle("Y");
03883 hYvsXRecoTrkLeakIn->GetYaxis()->CenterTitle();
03884
03885
03886 TH2F* hYvsXTrueTrkLeakOut=new TH2F("hYvsXTrueTrkLeakOut",
03887 "hYvsXTrueTrkLeakOut",
03888 100,-5,5,100,-5,5);
03889 hYvsXTrueTrkLeakOut->SetTitle("Y vs X");
03890 hYvsXTrueTrkLeakOut->GetXaxis()->SetTitle("X");
03891 hYvsXTrueTrkLeakOut->GetXaxis()->CenterTitle();
03892 hYvsXTrueTrkLeakOut->GetYaxis()->SetTitle("Y");
03893 hYvsXTrueTrkLeakOut->GetYaxis()->CenterTitle();
03894
03895
03896 TH2F* hYvsXRecoTrkLeakOut=new TH2F("hYvsXRecoTrkLeakOut",
03897 "hYvsXRecoTrkLeakOut",
03898 100,-5,5,100,-5,5);
03899 hYvsXRecoTrkLeakOut->SetTitle("Y vs X");
03900 hYvsXRecoTrkLeakOut->GetXaxis()->SetTitle("X");
03901 hYvsXRecoTrkLeakOut->GetXaxis()->CenterTitle();
03902 hYvsXRecoTrkLeakOut->GetYaxis()->SetTitle("Y");
03903 hYvsXRecoTrkLeakOut->GetYaxis()->CenterTitle();
03904
03905
03906
03907 Int_t entries=static_cast<Int_t>(this->GetEntries());
03908 Int_t startTimeSecs=2000000000;
03909 Int_t endTimeSecs=1;
03910 vector<Double_t> vPot;
03911 vPot.reserve(entries);
03912 vector<Double_t> vTime;
03913 vTime.reserve(entries);
03914 vector<Double_t> vPotBad;
03915 vPotBad.reserve(entries/10);
03916 vector<Double_t> vTimeBad;
03917 vTimeBad.reserve(entries/10);
03918
03919 vector<Double_t> vTimeNuMuEvt;
03920 vTimeNuMuEvt.reserve(entries/5);
03921 vector<Double_t> vTimeNuMuBarEvt;
03922 vTimeNuMuBarEvt.reserve(entries/5);
03923
03924
03925 map<Int_t,Double_t> deltaTs;
03926 map<Int_t,Int_t> evtsPerSlc;
03927
03928
03929 Int_t goodSpillCounter=0;
03930 Int_t badSpillCounter=0;
03931 Float_t totalPot=0;
03932 Float_t totalPotBad=0;
03933 Int_t evtCounter=0;
03934 Int_t evtNotlitime=0;
03935 Int_t evtNotIsLI=0;
03936 Int_t evtInFidVolCounter=0;
03937 Int_t evtWithTrkCounter=0;
03938 Int_t goodBestTrkCounter=0;
03939 Int_t goodTrkPassCounter=0;
03940 Int_t goodRecoEnCounter=0;
03941 Int_t goodUVDiffCounter=0;
03942 Int_t trkInFidVolCounter=0;
03943 Int_t goodFitSigQPCounter=0;
03944 Int_t goodFitProbCounter=0;
03945 Int_t goodFitChi2Counter=0;
03946 Int_t goodDeltaTCounter=0;
03947 Int_t goodEvtsPerSlcCounter=0;
03948 Int_t goodPIDCounter=0;
03949 Int_t goodBeyondTrkEndCounter=0;
03950 Int_t goodShwFractCounter=0;
03951 Int_t nuNQCounter=0;
03952 Int_t nuPQCounter=0;
03953
03954
03955
03956
03957
03958
03959 const NuCuts cuts;
03960 const NuGeneral general;
03961 const NuPlots plots;
03962 const NuReco reco;
03963 const NuExtraction ext;
03964 const NuBeam beam;
03965
03966
03967
03968 NuMadAnalysis& mad=*new NuMadAnalysis();
03969 MadDpID& dp=*new MadDpID();
03970
03971 BeamType::BeamType_t beamType=static_cast
03972 <BeamType::BeamType_t>(config.beamType);
03973
03974
03975 const NuLibrary& lib=NuLibrary::Instance();
03976
03980
03981 this->InitialiseLoopVariables();
03982
03983
03984 for(Int_t entry=0;entry<this->GetEntries();entry++){
03985
03986 this->SetLoopVariables(entry);
03987
03988
03989 const NtpStRecord& ntp=(*this->GetNtpStRecord());
03990
03991 const RecCandHeader& rec=ntp.GetHeader();
03992 MAXMSG("NuAnalysis",Msg::kInfo,5)
03993 <<"Found: run="<<rec.GetRun()
03994 <<", subrun="<<rec.GetSubRun()
03995 <<", detector="<<rec.GetVldContext().GetDetector()
03996 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
03997 <<endl;
03998
03999
04000 plots.FillGeneralHistos(ntp);
04001
04002
04003 NuEvent nuSnarl;
04004
04005
04006 ext.ExtractGeneralInfo(ntp,nuSnarl);
04007
04008 if (!beam.IsGoodSpillAndFillPot
04009 (this->GetNtpBDLiteRecord(),config,nuSnarl)) {
04010 badSpillCounter++;
04011 totalPotBad+=nuSnarl.pot;
04012 continue;
04013 }
04014 goodSpillCounter++;
04015 totalPot+=nuSnarl.pot;
04016
04017 dp.ChoosePDFs(rec.GetVldContext().GetDetector(),beamType);
04018 mad.SetEntry(&ntp);
04019
04020 set<Int_t> setNuMuBarInTrueFidCC;
04021
04022 VldTimeStamp vldts;
04023
04024 VldContext vc(rec.GetVldContext().GetDetector(),
04025 rec.GetVldContext().GetSimFlag(),vldts);
04026
04027 UgliGeomHandle ugh(vc);
04028
04029 TClonesArray& mcTca=(*ntp.mc);
04030 Int_t numInt=mcTca.GetEntries();
04031 MAXMSG("NuAnalysis",Msg::kInfo,10)
04032 <<"Number of entries in NtpMCTruth="
04033 <<mcTca.GetEntries()<<endl;
04034
04035
04036 NuEvent nu;
04037
04038
04039 ext.ExtractGeneralInfo(ntp,nu);
04040
04041
04043
04045 for (Int_t i=0;i<numInt;i++){
04046 const NtpMCTruth& mc=*(dynamic_cast<NtpMCTruth*>(mcTca[i]));
04047
04048 TVector3 xyz(mc.vtxx,mc.vtxy,mc.vtxz);
04049
04050 TVector3 uvz=ugh.xyz2uvz(xyz);
04051
04052 if (!(mc.iaction==1 && mc.inu==-14)) continue;
04053
04054 Bool_t isInFidVol=cuts.IsInFidVol(mc.vtxx,mc.vtxy,mc.vtxz,
04055 uvz.X(),uvz.Y(),0,0,
04056 nu.detector,nu.anaVersion,
04057 nu.releaseType,nu.simFlag);
04058
04059 if (!isInFidVol) continue;
04060 setNuMuBarInTrueFidCC.insert(mc.index);
04061
04062 hRecoEnNMBFull->Fill(fabs(mc.p4neu[3]));
04063 }
04064
04065 MAXMSG("NuAnalysis",Msg::kInfo,20)
04066 <<"Number of NuMuBar fiducial volume interactions in spill="
04067 <<setNuMuBarInTrueFidCC.size()<<endl;
04068
04069
04070 TClonesArray& thevtTca=(*ntp.thevt);
04071 const Int_t numthevts=thevtTca.GetEntriesFast();
04072 if (numthevts<=0) {
04073 MSG("NuAnalysis",Msg::kWarning)
04074 <<"No THEvents..."<<endl;
04075 continue;
04076 }
04077
04078 set<Int_t> setNuMuBarInRecoFidCC;
04079
04080 TClonesArray& evtTca=(*ntp.evt);
04081 const Int_t numEvts=evtTca.GetEntriesFast();
04082
04084
04086 for (Int_t ievt=0;ievt<numEvts;ievt++){
04087 const NtpSREvent& evt=
04088 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
04089 evtCounter++;
04090
04091 const NtpTHEvent& thevt=
04092 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
04093
04094
04095
04096 const NtpMCTruth& mc=
04097 *(dynamic_cast<NtpMCTruth*>(mcTca[thevt.neumc]));
04098
04099
04100 if (!(mc.iaction==1 && mc.inu==-14)) continue;
04101
04102
04103 NuEvent nu;
04104
04105
04106 ext.ExtractGeneralInfo(ntp,nu);
04107
04108 lib.ext.ExtractEvtInfo(evt,nu);
04109
04110 lib.ext.ExtractTrkInfo(ntp,evt,nu);
04111
04112 lib.ext.ExtractShwInfo(ntp,evt,nu);
04113
04114 Bool_t isInFidVol=cuts.IsInFidVol(evt.vtx.x,evt.vtx.y,
04115 evt.vtx.z,
04116 evt.vtx.u,evt.vtx.v,0,0,
04117 nu.detector,nu.anaVersion,
04118 nu.releaseType,nu.simFlag);
04119
04120
04121 set<Int_t>::iterator itT=setNuMuBarInRecoFidCC.find(thevt.neumc);
04122 if (itT==setNuMuBarInRecoFidCC.end()) {
04123 MAXMSG("NuAnalysis",Msg::kInfo,20)
04124 <<"Found reco event for mc="<<thevt.neumc<<endl;
04125 setNuMuBarInRecoFidCC.insert(thevt.neumc);
04126 }
04127 else {
04128 MAXMSG("NuAnalysis",Msg::kInfo,20)
04129 <<"Already found reco event for mc="<<thevt.neumc<<endl;
04130 continue;
04131 }
04132
04133
04134 Int_t bestTrack=lib.reco.GetBestTrack(nu);
04135 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
04136 (ntp,evt,bestTrack-1);
04137 Bool_t trkIsInFidVol=false;
04138 if (ptrk!=0) {
04139 const NtpSRTrack& trk=*ptrk;
04140 trkIsInFidVol=cuts.IsInFidVol(trk.vtx.x,trk.vtx.y,
04141 trk.vtx.z,
04142 trk.vtx.u,trk.vtx.v,0,0,
04143 nu.detector,nu.anaVersion,
04144 nu.releaseType,nu.simFlag);
04145 trkIsInFidVol=true;
04146 }
04147
04148
04149 set<Int_t>::iterator itR=setNuMuBarInTrueFidCC.find(thevt.neumc);
04150 if (itR==setNuMuBarInTrueFidCC.end()) {
04151 MAXMSG("NuAnalysis",Msg::kInfo,20)
04152 <<"ievt="<<ievt<<" has no true event in fid"<<endl;
04153 if (isInFidVol){
04154 hRecoEnNMBLeakIn->Fill(fabs(mc.p4neu[3]));
04155 hYvsXRecoLeakIn->Fill(evt.vtx.x,evt.vtx.y);
04156 hYvsXTrueLeakIn->Fill(mc.vtxx,mc.vtxy);
04157 if (trkIsInFidVol){
04158 const NtpSRTrack& trk=*ptrk;
04159 hRecoEnNMBTrkLeakIn->Fill(fabs(mc.p4neu[3]));
04160 hYvsXRecoTrkLeakIn->Fill(trk.vtx.x,trk.vtx.y);
04161 hYvsXTrueTrkLeakIn->Fill(mc.vtxx,mc.vtxy);
04162 }
04163 }
04164 continue;
04165 }
04166
04167
04168 hRecoEnNMBReco->Fill(fabs(mc.p4neu[3]));
04169
04170
04171
04172
04173 if (!trkIsInFidVol) {
04174
04175 if (ptrk){
04176 const NtpSRTrack& trk=*ptrk;
04177 hRecoEnNMBTrkLeakOut->Fill(fabs(mc.p4neu[3]));
04178 hYvsXRecoTrkLeakOut->Fill(trk.vtx.x,trk.vtx.y);
04179 hYvsXTrueTrkLeakOut->Fill(mc.vtxx,mc.vtxy);
04180 }
04181 else{
04182 MAXMSG("NuAnalysis",Msg::kInfo,20)
04183 <<"No track, using evt vtx"<<endl;
04184 hRecoEnNMBTrkLeakOut->Fill(fabs(mc.p4neu[3]));
04185 hYvsXRecoTrkLeakOut->Fill(evt.vtx.x,evt.vtx.y);
04186 hYvsXTrueTrkLeakOut->Fill(mc.vtxx,mc.vtxy);
04187 }
04188 }
04189 if (!isInFidVol) {
04190
04191 hRecoEnNMBLeakOut->Fill(fabs(mc.p4neu[3]));
04192 hYvsXRecoLeakOut->Fill(evt.vtx.x,evt.vtx.y);
04193 hYvsXTrueLeakOut->Fill(mc.vtxx,mc.vtxy);
04194 continue;
04195 }
04196 evtInFidVolCounter++;
04197 hRecoEnNMBFid->Fill(fabs(mc.p4neu[3]));
04198
04199
04200 if (evt.ntrack<1) continue;
04201
04202 evtWithTrkCounter++;
04203 hRecoEnNMBTrack->Fill(fabs(mc.p4neu[3]));
04204
04205
04206 if (ptrk==0) continue;
04207 const NtpSRTrack& trk=*ptrk;
04208 goodBestTrkCounter++;
04209
04210
04211 if (!trk.fit.pass) continue;
04212 goodTrkPassCounter++;
04213 hRecoEnNMBPass->Fill(fabs(mc.p4neu[3]));
04214
04215
04216 if (trk.momentum.qp<=0) continue;
04217 hRecoEnNMBQCut->Fill(fabs(mc.p4neu[3]));
04218
04220
04221 lib.reco.GetEvtEnergy(nu);
04222 goodRecoEnCounter++;
04223
04224 reco.GetTruthInfo(ntp,evt,nu);
04225
04226 nu.charge=trk.momentum.qp<0 ? -1 : 1;
04227 if (trk.momentum.qp==0) nu.charge=0;
04228
04229
04230 if (!cuts.IsGoodSigmaQP_QP(nu)) continue;
04231 goodFitSigQPCounter++;
04232 hRecoEnNMBSigQP->Fill(fabs(mc.p4neu[3]));
04233
04234
04235 if (!cuts.IsGoodFitProb(nu)) continue;
04236 goodFitProbCounter++;
04237 hRecoEnNMBProb->Fill(fabs(mc.p4neu[3]));
04238
04239 nu.dpID=dp.CalcPID(&mad,ievt,0);
04240 if (nu.dpID>=-0.1) {
04241 hRecoEnNMBDpID01->Fill(fabs(mc.p4neu[3]));
04242 }
04243
04244 if (nu.dpID>=0.4) {
04245 hRecoEnNMBDpID04->Fill(fabs(mc.p4neu[3]));
04246 }
04247 }
04248 }
04249
04253
04254 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
04255
04256
04257 Int_t nf=this->GetNumFilesAddedToChain();
04258 TH1F* hRun=dynamic_cast<TH1F*>(gROOT->FindObject("hRun"));
04259 Float_t nr=-1;
04260 if (hRun) nr=hRun->GetEntries();
04261 TH1F* hNumFiles=new TH1F("hNumFiles","hNumFiles",15,-5,10);
04262 hNumFiles->Fill(1,nf);
04263 MSG("NuAnalysis",Msg::kInfo)
04264 <<"Files added to chain="<<hNumFiles->GetMaximum()
04265 <<", runs found="<<nr<<endl;
04266
04267 MSG("NuAnalysis",Msg::kInfo)
04268 <<"Found start time and end time: "
04269 <<startTimeSecs<<" -> "<<endTimeSecs<<endl;
04270 general.EpochTo1995(endTimeSecs);
04271 general.EpochTo1995(startTimeSecs);
04272
04273 MSG("NuAnalysis",Msg::kInfo)<<"Filling vs. Time plots..."<<endl;
04274 TH1F* hSpillsVsTime=new TH1F("hSpillsVsTime","hSpillsVsTime",
04275 100,startTimeSecs,endTimeSecs);
04276 general.TH1FFill(hSpillsVsTime,vTime);
04277 general.SetGraphAxis(hSpillsVsTime->GetXaxis(),
04278 startTimeSecs,endTimeSecs);
04279 hSpillsVsTime->GetXaxis()->CenterTitle();
04280 hSpillsVsTime->GetYaxis()->SetTitle("Spills");
04281 hSpillsVsTime->GetYaxis()->CenterTitle();
04282
04283 TH1F* hSpillsBadVsTime=new TH1F("hSpillsBadVsTime","hSpillsBadVsTime",
04284 100,startTimeSecs,endTimeSecs);
04285 general.TH1FFill(hSpillsBadVsTime,vTimeBad);
04286 general.SetGraphAxis(hSpillsBadVsTime->GetXaxis(),
04287 startTimeSecs,endTimeSecs);
04288 hSpillsBadVsTime->GetXaxis()->CenterTitle();
04289 hSpillsBadVsTime->GetYaxis()->SetTitle("Spills");
04290 hSpillsBadVsTime->GetYaxis()->CenterTitle();
04291
04292 TH1F* hNuMuVsTime=new TH1F("hNuMuVsTime","hNuMuVsTime",
04293 100,startTimeSecs,endTimeSecs);
04294 general.TH1FFill(hNuMuVsTime,vTimeNuMuEvt);
04295 general.SetGraphAxis(hNuMuVsTime->GetXaxis(),
04296 startTimeSecs,endTimeSecs);
04297 hNuMuVsTime->GetXaxis()->CenterTitle();
04298 hNuMuVsTime->GetYaxis()->SetTitle("#nu_#mu events");
04299 hNuMuVsTime->GetYaxis()->CenterTitle();
04300
04301 TH1F* hNuMuBarVsTime=new TH1F("hNuMuBarVsTime","hNuMuBarVsTime",
04302 100,startTimeSecs,endTimeSecs);
04303 general.TH1FFill(hNuMuBarVsTime,vTimeNuMuBarEvt);
04304 general.SetGraphAxis(hNuMuBarVsTime->GetXaxis(),
04305 startTimeSecs,endTimeSecs);
04306 hNuMuBarVsTime->GetXaxis()->CenterTitle();
04307 hNuMuBarVsTime->GetYaxis()->SetTitle("#bar{#nu_#mu} events");
04308 hNuMuBarVsTime->GetYaxis()->CenterTitle();
04309
04310 TH1F* hPotVsTime=new TH1F("hPotVsTime","hPotVsTime",
04311 100,startTimeSecs,endTimeSecs);
04312 general.TH1FFill(hPotVsTime,vTime,vPot);
04313 general.SetGraphAxis(hPotVsTime->GetXaxis(),
04314 startTimeSecs,endTimeSecs);
04315 hPotVsTime->GetXaxis()->CenterTitle();
04316 hPotVsTime->GetYaxis()->SetTitle("POT");
04317 hPotVsTime->GetYaxis()->CenterTitle();
04318
04319 TH1F* hPotBadVsTime=new TH1F("hPotBadVsTime","hPotBadVsTime",
04320 100,startTimeSecs,endTimeSecs);
04321 general.TH1FFill(hPotBadVsTime,vTimeBad,vPotBad);
04322 general.SetGraphAxis(hPotBadVsTime->GetXaxis(),
04323 startTimeSecs,endTimeSecs);
04324 hPotBadVsTime->GetXaxis()->CenterTitle();
04325 hPotBadVsTime->GetYaxis()->SetTitle("POT");
04326 hPotBadVsTime->GetYaxis()->CenterTitle();
04327
04328 Float_t totalPotHist=0;
04329 Float_t totalPotBadHist=0;
04330 cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
04331 Float_t pcPotRejected=0;
04332 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
04333
04334
04335 cout<<"Number of spills ="<<this->GetEntries()<<endl
04336 <<"Number of good spills="<<goodSpillCounter<<endl
04337 <<"Number of bad spills ="<<badSpillCounter<<endl
04338 <<"Sum of good+bad spill="<<badSpillCounter+goodSpillCounter<<endl
04339 <<endl
04340 <<"Total POT good ="<<totalPot*1e12<<endl
04341 <<"Total POT bad ="<<totalPotBad*1e12<<endl
04342 <<"Total POT good (hist)="<<totalPotHist*1e12<<endl
04343 <<"Total POT bad (hist)="<<totalPotBadHist*1e12
04344 <<" ("<<pcPotRejected<<"%)"<<endl;
04345
04346 MSG("NuAnalysis",Msg::kInfo)
04347 <<endl<<"*************************************"<<endl
04348 <<"Total Events ="<<evtCounter<<endl
04349 <<"Evt not litime ="<<evtNotlitime<<endl
04350 <<"Evt not IsLI ="<<evtNotIsLI<<endl
04351 <<"Evts in Fid ="<<evtInFidVolCounter<<endl
04352 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
04353 <<"Good Best Trk ="<<goodBestTrkCounter<<endl
04354 <<"Track Pass ="<<goodTrkPassCounter<<endl
04355 <<"Good reco energy="<<goodRecoEnCounter<<endl
04356 <<"Good UVDiff ="<<goodUVDiffCounter<<endl
04357 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
04358 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
04359 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
04360 <<"Good Fit Chi2 ="<<goodFitChi2Counter<<endl
04361 <<"Good Delta T ="<<goodDeltaTCounter<<endl
04362 <<"Good EvtPerSlc ="<<goodEvtsPerSlcCounter<<endl
04363 <<"Good DP ID ="<<goodPIDCounter<<endl
04364 <<"Good beyond tEnd="<<goodBeyondTrkEndCounter<<endl
04365 <<"Good shw fract ="<<goodShwFractCounter<<endl
04366 <<"*************************************"<<endl;
04367
04368 Float_t pcPQ=0;
04369 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
04370 (nuPQCounter+nuNQCounter);
04371 MSG("NuAnalysis",Msg::kInfo)
04372 <<endl
04373 <<"Neutrinos with NQ="<<nuNQCounter<<endl
04374 <<"Neutrinos with PQ="<<nuPQCounter
04375 <<" ("<<pcPQ<<"%)"<<endl
04376 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
04377 <<endl;
04378
04379 MSG("NuAnalysis",Msg::kInfo)
04380 <<" ** Finished Efficiencies method **"<<endl;
04381 }
04382
04383
04384
04385 void NuAnalysis::ChargeSignCut()
04386 {
04387 MSG("NuAnalysis",Msg::kInfo)
04388 <<" ** Running ChargeSignCut method... **"<<endl;
04389
04390 NuConfig config;
04391 this->ExtractConfig(config);
04392
04393
04394 string sFilePrefix="NMBCutQ"+config.sBeamType;
04395 fOutFile=this->OpenFile(config,sFilePrefix.c_str());
04396
04397 TH1F* hNuMuBarRecoEn=new TH1F("hNuMuBarRecoEn","hNuMuBarRecoEn",
04398 4*352,-32,320);
04399 hNuMuBarRecoEn->GetXaxis()->SetTitle("Energy (GeV)");
04400 hNuMuBarRecoEn->GetXaxis()->CenterTitle();
04401 hNuMuBarRecoEn->GetYaxis()->SetTitle("");
04402 hNuMuBarRecoEn->GetYaxis()->CenterTitle();
04403 hNuMuBarRecoEn->SetFillColor(0);
04404 hNuMuBarRecoEn->SetLineColor(2);
04405
04406
04407 TH1F* hNuMuRecoEn=new TH1F("hNuMuRecoEn","hNuMuRecoEn",4*352,-32,320);
04408 hNuMuRecoEn->GetXaxis()->SetTitle("Energy (GeV)");
04409 hNuMuRecoEn->GetXaxis()->CenterTitle();
04410 hNuMuRecoEn->GetYaxis()->SetTitle("");
04411 hNuMuRecoEn->GetYaxis()->CenterTitle();
04412 hNuMuRecoEn->SetFillColor(0);
04413 hNuMuRecoEn->SetLineColor(1);
04414 hNuMuRecoEn->SetLineWidth(2);
04415
04416
04417 TH1F* hNuMuBarRecoY=new TH1F("hNuMuBarRecoY","hNuMuBarRecoY",
04418 4*144,-0.16,1.28);
04419 hNuMuBarRecoY->GetXaxis()->SetTitle("y");
04420 hNuMuBarRecoY->GetXaxis()->CenterTitle();
04421 hNuMuBarRecoY->GetYaxis()->SetTitle("");
04422 hNuMuBarRecoY->GetYaxis()->CenterTitle();
04423 hNuMuBarRecoY->SetFillColor(0);
04424 hNuMuBarRecoY->SetLineColor(2);
04425
04426
04427 TH1F* hNuMuRecoY=new TH1F("hNuMuRecoY","hNuMuRecoY",4*144,-0.16,1.28);
04428 hNuMuRecoY->GetXaxis()->SetTitle("y");
04429 hNuMuRecoY->GetXaxis()->CenterTitle();
04430 hNuMuRecoY->GetYaxis()->SetTitle("");
04431 hNuMuRecoY->GetYaxis()->CenterTitle();
04432 hNuMuRecoY->SetFillColor(0);
04433 hNuMuRecoY->SetLineColor(1);
04434 hNuMuRecoY->SetLineWidth(2);
04435
04436
04437 TH1F* hDpIDAll=new TH1F("hDpIDAll","hDpIDAll",4*160,-1.6,1.6);
04438 hDpIDAll->GetXaxis()->SetTitle("PID (from MadDpID)");
04439 hDpIDAll->GetXaxis()->CenterTitle();
04440 hDpIDAll->GetYaxis()->SetTitle("");
04441 hDpIDAll->GetYaxis()->CenterTitle();
04442 hDpIDAll->SetFillColor(0);
04443 hDpIDAll->SetLineColor(1);
04444 hDpIDAll->SetLineWidth(2);
04445
04446
04447 TH1F* hDpIDPQ=new TH1F("hDpIDPQ","hDpIDPQ",4*160,-1.6,1.6);
04448 hDpIDPQ->GetXaxis()->SetTitle("PID (from MadDpID)");
04449 hDpIDPQ->GetXaxis()->CenterTitle();
04450 hDpIDPQ->GetYaxis()->SetTitle("");
04451 hDpIDPQ->GetYaxis()->CenterTitle();
04452 hDpIDPQ->SetFillColor(0);
04453 hDpIDPQ->SetLineColor(1);
04454 hDpIDPQ->SetLineWidth(2);
04455
04456
04457 TH1F* hDpIDNQ=new TH1F("hDpIDNQ","hDpIDNQ",4*160,-1.6,1.6);
04458 hDpIDNQ->GetXaxis()->SetTitle("PID (from MadDpID)");
04459 hDpIDNQ->GetXaxis()->CenterTitle();
04460 hDpIDNQ->GetYaxis()->SetTitle("");
04461 hDpIDNQ->GetYaxis()->CenterTitle();
04462 hDpIDNQ->SetFillColor(0);
04463 hDpIDNQ->SetLineColor(1);
04464 hDpIDNQ->SetLineWidth(2);
04465
04466
04467 TH1F* hDpIDPQN_1=new TH1F("hDpIDPQN_1","hDpIDPQN_1",4*160,-1.6,1.6);
04468 hDpIDPQN_1->GetXaxis()->SetTitle("PID (from MadDpID)");
04469 hDpIDPQN_1->GetXaxis()->CenterTitle();
04470 hDpIDPQN_1->GetYaxis()->SetTitle("");
04471 hDpIDPQN_1->GetYaxis()->CenterTitle();
04472 hDpIDPQN_1->SetFillColor(0);
04473 hDpIDPQN_1->SetLineColor(1);
04474 hDpIDPQN_1->SetLineWidth(2);
04475
04476
04477 TH1F* hDpIDNQN_1=new TH1F("hDpIDNQN_1","hDpIDNQN_1",4*160,-1.6,1.6);
04478 hDpIDNQN_1->GetXaxis()->SetTitle("PID (from MadDpID)");
04479 hDpIDNQN_1->GetXaxis()->CenterTitle();
04480 hDpIDNQN_1->GetYaxis()->SetTitle("");
04481 hDpIDNQN_1->GetYaxis()->CenterTitle();
04482 hDpIDNQN_1->SetFillColor(0);
04483 hDpIDNQN_1->SetLineColor(1);
04484 hDpIDNQN_1->SetLineWidth(2);
04485
04486
04487
04488 Int_t entries=static_cast<Int_t>(this->GetEntries());
04489 Int_t startTimeSecs=2000000000;
04490 Int_t endTimeSecs=1;
04491 vector<Double_t> vPot;
04492 vPot.reserve(entries);
04493 vector<Double_t> vTime;
04494 vTime.reserve(entries);
04495 vector<Double_t> vPotBad;
04496 vPotBad.reserve(entries/10);
04497 vector<Double_t> vTimeBad;
04498 vTimeBad.reserve(entries/10);
04499
04500 vector<Double_t> vTimeNuMuEvt;
04501 vTimeNuMuEvt.reserve(entries/5);
04502 vector<Double_t> vTimeNuMuBarEvt;
04503 vTimeNuMuBarEvt.reserve(entries/5);
04504
04505
04506 map<Int_t,Double_t> deltaTs;
04507 map<Int_t,Int_t> evtsPerSlc;
04508
04509
04510 Int_t goodSpillCounter=0;
04511 Int_t badSpillCounter=0;
04512 Float_t totalPot=0;
04513 Float_t totalPotBad=0;
04514 Int_t evtCounter=0;
04515 Int_t evtNotlitime=0;
04516 Int_t evtNotIsLI=0;
04517 Int_t evtInFidVolCounter=0;
04518 Int_t evtWithTrkCounter=0;
04519 Int_t goodBestTrkCounter=0;
04520 Int_t goodTrkPassCounter=0;
04521 Int_t goodRecoEnCounter=0;
04522 Int_t goodUVDiffCounter=0;
04523 Int_t trkInFidVolCounter=0;
04524 Int_t goodFitSigQPCounter=0;
04525 Int_t goodFitProbCounter=0;
04526 Int_t goodFitChi2Counter=0;
04527 Int_t goodDeltaTCounter=0;
04528 Int_t goodEvtsPerSlcCounter=0;
04529 Int_t goodPIDCounter=0;
04530 Int_t goodBeyondTrkEndCounter=0;
04531 Int_t goodShwFractCounter=0;
04532 Int_t nuNQCounter=0;
04533 Int_t nuPQCounter=0;
04534
04535 string sTxtFilePrefix="nmbQCut"+config.sBeamType;
04536 ofstream& nmbTxt=*(this->OpenTxtFile(config,
04537 sTxtFilePrefix.c_str()));
04538 nmbTxt<<"# All NuMuBar like events (with standard CC cuts)"<<endl;
04539 nmbTxt<<"#"<<endl;
04540
04541 sTxtFilePrefix="nmbQCut05"+config.sBeamType;
04542 ofstream& nmb05Txt=*(this->OpenTxtFile(config,
04543 sTxtFilePrefix.c_str()));
04544 nmb05Txt<<"# NuMuBar like events below 5 GeV (with standard CC cuts)"
04545 <<endl;
04546 nmb05Txt<<"#"<<endl;
04547
04548 sTxtFilePrefix="nmbQCut02"+config.sBeamType;
04549 ofstream& nmb02Txt=*(this->OpenTxtFile(config,
04550 sTxtFilePrefix.c_str()));
04551 nmb02Txt<<"# NuMuBar like events below 2 GeV (with standard CC cuts)"
04552 <<endl;
04553 nmb02Txt<<"#"<<endl;
04554
04555
04556 const NuCuts cuts;
04557 const NuGeneral general;
04558 const NuPlots plots;
04559 const NuReco reco;
04560 const NuExtraction ext;
04561 const NuBeam beam;
04562
04563
04564
04565 NuMadAnalysis& mad=*new NuMadAnalysis();
04566 MadDpID& dp=*new MadDpID();
04567
04568 BeamType::BeamType_t beamType=static_cast
04569 <BeamType::BeamType_t>(config.beamType);
04570
04571
04572 const NuLibrary& lib=NuLibrary::Instance();
04573
04577
04578 this->InitialiseLoopVariables();
04579
04580
04581 for(Int_t entry=0;entry<this->GetEntries();entry++){
04582
04583 this->SetLoopVariables(entry);
04584
04585
04586 const NtpStRecord& ntp=(*this->GetNtpStRecord());
04587
04588 const RecCandHeader& rec=ntp.GetHeader();
04589 MAXMSG("NuAnalysis",Msg::kInfo,5)
04590 <<"Found: run="<<rec.GetRun()
04591 <<", subrun="<<rec.GetSubRun()
04592 <<", detector="<<rec.GetVldContext().GetDetector()
04593 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
04594 <<endl;
04595
04596
04597 plots.FillGeneralHistos(ntp);
04598
04599
04600 NuEvent nuSnarl;
04601
04602
04603 ext.ExtractGeneralInfo(ntp,nuSnarl);
04604
04605 if (!beam.IsGoodSpillAndFillPot
04606 (this->GetNtpBDLiteRecord(),config,nuSnarl)) {
04607 badSpillCounter++;
04608 totalPotBad+=nuSnarl.pot;
04609 continue;
04610 }
04611 goodSpillCounter++;
04612 totalPot+=nuSnarl.pot;
04613
04614 dp.ChoosePDFs(rec.GetVldContext().GetDetector(),beamType);
04615 mad.SetEntry(&ntp);
04616
04617 TClonesArray& evtTca=(*ntp.evt);
04618 const Int_t numEvts=evtTca.GetEntriesFast();
04619
04621
04623 for (Int_t ievt=0;ievt<numEvts;ievt++){
04624 const NtpSREvent& evt=
04625 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
04626 evtCounter++;
04627
04628 const NtpSREventSummary& evthdr=ntp.evthdr;
04629
04630 if (evthdr.litime!=-1) {
04631 MAXMSG("NuAnalysis",Msg::kInfo,500)
04632 <<"skipping event: litime="<<evthdr.litime<<endl;
04633 continue;
04634 }
04635 evtNotlitime++;
04636
04637 Bool_t isLI=LISieve::IsLI(ntp);
04638 if (isLI){
04639 MAXMSG("NuAnalysis",Msg::kInfo,5)
04640 <<endl<<endl<<endl
04641 <<"Found LI"
04642 <<", run="<<config.run
04643 <<", entry="<<entry
04644 <<", event="<<ievt
04645 <<", litime="<<evthdr.litime
04646 <<endl<<endl<<endl;
04647 }
04648 else {
04649 MAXMSG("NuAnalysis",Msg::kDebug,50)
04650 <<"Not LI"
04651 <<", run="<<config.run
04652 <<", entry="<<entry
04653 <<", event="<<ievt
04654 <<", litime="<<evthdr.litime
04655 <<endl;
04656 }
04657 if (isLI) continue;
04658 evtNotIsLI++;
04659
04660
04661 NuEvent nu;
04662
04663
04664 ext.ExtractGeneralInfo(ntp,nu);
04665
04666 Bool_t isInFidVol=cuts.IsInFidVol(evt.vtx.x,evt.vtx.y,
04667 evt.vtx.z,
04668 evt.vtx.u,evt.vtx.v,0,0,
04669 nu.detector,nu.anaVersion,
04670 nu.releaseType,nu.simFlag);
04671
04672 if (!isInFidVol) continue;
04673 evtInFidVolCounter++;
04674
04675 nu.dpID=dp.CalcPID(&mad,ievt,0);
04676 MAXMSG("NuAnalysis",Msg::kDebug,50)
04677 <<"dpID="<<nu.dpID<<endl;
04678 hDpIDAll->Fill(nu.dpID);
04679
04680
04681 if (evt.ntrack<1) continue;
04682
04683 evtWithTrkCounter++;
04684
04685
04686 Int_t bestTrack=lib.reco.GetBestTrack(nu);
04687 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
04688 (ntp,evt,bestTrack-1);
04689 if (ptrk==0) continue;
04690 const NtpSRTrack& trk=*ptrk;
04691 goodBestTrkCounter++;
04692
04693
04694 if (!trk.fit.pass) continue;
04695 goodTrkPassCounter++;
04696
04697
04699
04700 lib.reco.GetEvtEnergy(nu);
04701 goodRecoEnCounter++;
04702
04703
04704 reco.GetTruthInfo(ntp,evt,nu);
04705
04706
04707 nu.charge=trk.momentum.qp<0 ? -1 : +1;
04708 if (trk.momentum.qp==0) nu.charge=0;
04709
04710
04711
04712
04713
04714
04715
04716
04717
04718
04719
04720
04721
04722
04723
04724
04725
04726
04727
04728 if (nu.dpID<=-0.1) continue;
04729 goodPIDCounter++;
04730
04731 if (trk.momentum.qp<0) {
04732 hNuMuRecoEn->Fill(nu.energy);
04733 hNuMuRecoY->Fill(nu.y);
04734 hDpIDNQ->Fill(nu.dpID);
04735 nuNQCounter++;
04736 }
04737 else if (trk.momentum.qp>0) {
04738 string siaction="CC";
04739 if (nu.iaction==0) siaction="NC";
04740
04741 nmbTxt<<"#"<<endl;
04742 nmbTxt<<"# recoEn="<<nu.energy
04743 <<", trkEn="<<nu.trkEn
04744 <<", shwEn="<<nu.shwEn
04745 <<", recoy="<<nu.y
04746 <<", true inu="<<nu.inu
04747 <<" "<<siaction
04748 <<endl;
04749 nmbTxt<<nu.run<<" "<<nu.snarl<<" "<<nu.evt<<endl;
04750 if (nu.energy<5){
04751 nmb05Txt<<"#"<<endl;
04752 nmb05Txt<<"# recoEn="<<nu.energy
04753 <<", trkEn="<<nu.trkEn
04754 <<", shwEn="<<nu.shwEn
04755 <<", recoy="<<nu.y
04756 <<", true inu="<<nu.inu
04757 <<" "<<siaction
04758 <<endl;
04759 nmb05Txt<<nu.run<<" "<<nu.snarl<<" "<<nu.evt<<endl;
04760 }
04761
04762 if (nu.energy<2){
04763 nmb02Txt<<"#"<<endl;
04764 nmb02Txt<<"# recoEn="<<nu.energy
04765 <<", trkEn="<<nu.trkEn
04766 <<", shwEn="<<nu.shwEn
04767 <<", recoy="<<nu.y
04768 <<", true inu="<<nu.inu
04769 <<" "<<siaction
04770 <<endl;
04771 nmb02Txt<<nu.run<<" "<<nu.snarl<<" "<<nu.evt<<endl;
04772 }
04773
04774 hNuMuBarRecoEn->Fill(nu.energy);
04775 hNuMuBarRecoY->Fill(nu.y);
04776 hDpIDPQ->Fill(nu.dpID);
04777 nuPQCounter++;
04778 }
04779 else cout<<"ahhh, zero qp"<<endl;
04780
04781 plots.FillXYZHistos(nu);
04782 plots.FillContainmentHistos(nu);
04783 plots.FillTruePIDHistos(nu);
04784 plots.FillTrackResponseHistos(ntp,trk,nu);
04785 plots.FillRangeCurvCompHistos(nu);
04786 }
04787 }
04788
04792
04793 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
04794
04795
04796 Int_t nf=this->GetNumFilesAddedToChain();
04797 TH1F* hRun=dynamic_cast<TH1F*>(gROOT->FindObject("hRun"));
04798 Float_t nr=-1;
04799 if (hRun) nr=hRun->GetEntries();
04800 TH1F* hNumFiles=new TH1F("hNumFiles","hNumFiles",15,-5,10);
04801 hNumFiles->Fill(1,nf);
04802 MSG("NuAnalysis",Msg::kInfo)
04803 <<"Files added to chain="<<hNumFiles->GetMaximum()
04804 <<", runs found="<<nr<<endl;
04805
04806 MSG("NuAnalysis",Msg::kInfo)
04807 <<"Found start time and end time: "
04808 <<startTimeSecs<<" -> "<<endTimeSecs<<endl;
04809 general.EpochTo1995(endTimeSecs);
04810 general.EpochTo1995(startTimeSecs);
04811
04812 MSG("NuAnalysis",Msg::kInfo)<<"Filling vs. Time plots..."<<endl;
04813 TH1F* hSpillsVsTime=new TH1F("hSpillsVsTime","hSpillsVsTime",
04814 100,startTimeSecs,endTimeSecs);
04815 general.TH1FFill(hSpillsVsTime,vTime);
04816 general.SetGraphAxis(hSpillsVsTime->GetXaxis(),
04817 startTimeSecs,endTimeSecs);
04818 hSpillsVsTime->GetXaxis()->CenterTitle();
04819 hSpillsVsTime->GetYaxis()->SetTitle("Spills");
04820 hSpillsVsTime->GetYaxis()->CenterTitle();
04821
04822 TH1F* hSpillsBadVsTime=new TH1F("hSpillsBadVsTime","hSpillsBadVsTime",
04823 100,startTimeSecs,endTimeSecs);
04824 general.TH1FFill(hSpillsBadVsTime,vTimeBad);
04825 general.SetGraphAxis(hSpillsBadVsTime->GetXaxis(),
04826 startTimeSecs,endTimeSecs);
04827 hSpillsBadVsTime->GetXaxis()->CenterTitle();
04828 hSpillsBadVsTime->GetYaxis()->SetTitle("Spills");
04829 hSpillsBadVsTime->GetYaxis()->CenterTitle();
04830
04831 TH1F* hNuMuVsTime=new TH1F("hNuMuVsTime","hNuMuVsTime",
04832 100,startTimeSecs,endTimeSecs);
04833 general.TH1FFill(hNuMuVsTime,vTimeNuMuEvt);
04834 general.SetGraphAxis(hNuMuVsTime->GetXaxis(),
04835 startTimeSecs,endTimeSecs);
04836 hNuMuVsTime->GetXaxis()->CenterTitle();
04837 hNuMuVsTime->GetYaxis()->SetTitle("#nu_#mu events");
04838 hNuMuVsTime->GetYaxis()->CenterTitle();
04839
04840 TH1F* hNuMuBarVsTime=new TH1F("hNuMuBarVsTime","hNuMuBarVsTime",
04841 100,startTimeSecs,endTimeSecs);
04842 general.TH1FFill(hNuMuBarVsTime,vTimeNuMuBarEvt);
04843 general.SetGraphAxis(hNuMuBarVsTime->GetXaxis(),
04844 startTimeSecs,endTimeSecs);
04845 hNuMuBarVsTime->GetXaxis()->CenterTitle();
04846 hNuMuBarVsTime->GetYaxis()->SetTitle("#bar{#nu_#mu} events");
04847 hNuMuBarVsTime->GetYaxis()->CenterTitle();
04848
04849 TH1F* hPotVsTime=new TH1F("hPotVsTime","hPotVsTime",
04850 100,startTimeSecs,endTimeSecs);
04851 general.TH1FFill(hPotVsTime,vTime,vPot);
04852 general.SetGraphAxis(hPotVsTime->GetXaxis(),
04853 startTimeSecs,endTimeSecs);
04854 hPotVsTime->GetXaxis()->CenterTitle();
04855 hPotVsTime->GetYaxis()->SetTitle("POT");
04856 hPotVsTime->GetYaxis()->CenterTitle();
04857
04858 TH1F* hPotBadVsTime=new TH1F("hPotBadVsTime","hPotBadVsTime",
04859 100,startTimeSecs,endTimeSecs);
04860 general.TH1FFill(hPotBadVsTime,vTimeBad,vPotBad);
04861 general.SetGraphAxis(hPotBadVsTime->GetXaxis(),
04862 startTimeSecs,endTimeSecs);
04863 hPotBadVsTime->GetXaxis()->CenterTitle();
04864 hPotBadVsTime->GetYaxis()->SetTitle("POT");
04865 hPotBadVsTime->GetYaxis()->CenterTitle();
04866
04867 Float_t totalPotHist=0;
04868 Float_t totalPotBadHist=0;
04869 cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
04870 Float_t pcPotRejected=0;
04871 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
04872
04873
04874 cout<<"Number of spills ="<<this->GetEntries()<<endl
04875 <<"Number of good spills="<<goodSpillCounter<<endl
04876 <<"Number of bad spills ="<<badSpillCounter<<endl
04877 <<"Sum of good+bad spill="<<badSpillCounter+goodSpillCounter<<endl
04878 <<endl
04879 <<"Total POT good ="<<totalPot*1e12<<endl
04880 <<"Total POT bad ="<<totalPotBad*1e12<<endl
04881 <<"Total POT good (hist)="<<totalPotHist*1e12<<endl
04882 <<"Total POT bad (hist)="<<totalPotBadHist*1e12
04883 <<" ("<<pcPotRejected<<"%)"<<endl;
04884
04885 MSG("NuAnalysis",Msg::kInfo)
04886 <<endl<<"*************************************"<<endl
04887 <<"Total Events ="<<evtCounter<<endl
04888 <<"Evt not litime ="<<evtNotlitime<<endl
04889 <<"Evt not IsLI ="<<evtNotIsLI<<endl
04890 <<"Evts in Fid ="<<evtInFidVolCounter<<endl
04891 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
04892 <<"Good Best Trk ="<<goodBestTrkCounter<<endl
04893 <<"Track Pass ="<<goodTrkPassCounter<<endl
04894 <<"Good reco energy="<<goodRecoEnCounter<<endl
04895 <<"Good UVDiff ="<<goodUVDiffCounter<<endl
04896 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
04897 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
04898 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
04899 <<"Good Fit Chi2 ="<<goodFitChi2Counter<<endl
04900 <<"Good Delta T ="<<goodDeltaTCounter<<endl
04901 <<"Good EvtPerSlc ="<<goodEvtsPerSlcCounter<<endl
04902 <<"Good DP ID ="<<goodPIDCounter<<endl
04903 <<"Good beyond tEnd="<<goodBeyondTrkEndCounter<<endl
04904 <<"Good shw fract ="<<goodShwFractCounter<<endl
04905 <<"*************************************"<<endl;
04906
04907 Float_t pcPQ=0;
04908 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
04909 (nuPQCounter+nuNQCounter);
04910 MSG("NuAnalysis",Msg::kInfo)
04911 <<endl
04912 <<"Neutrinos with NQ="<<nuNQCounter<<endl
04913 <<"Neutrinos with PQ="<<nuPQCounter
04914 <<" ("<<pcPQ<<"%)"<<endl
04915 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
04916 <<endl;
04917
04918 MSG("NuAnalysis",Msg::kInfo)
04919 <<" ** Finished ChargeSignCut method **"<<endl;
04920 }
04921
04922
04923
04924 void NuAnalysis::EnergySpect()
04925 {
04926 MSG("NuAnalysis",Msg::kInfo)
04927 <<" ** Running EnergySpect method... **"<<endl;
04928
04929 NuConfig config;
04930 this->ExtractConfig(config);
04931
04932
04933 string sFilePrefix="NMBEnSp"+config.sBeamType;
04934 fOutFile=this->OpenFile(config,sFilePrefix.c_str());
04935
04936 TH1F* hNuMuBarRecoEn=new TH1F("hNuMuBarRecoEn","hNuMuBarRecoEn",
04937 4*352,-32,320);
04938 hNuMuBarRecoEn->GetXaxis()->SetTitle("Energy (GeV)");
04939 hNuMuBarRecoEn->GetXaxis()->CenterTitle();
04940 hNuMuBarRecoEn->GetYaxis()->SetTitle("");
04941 hNuMuBarRecoEn->GetYaxis()->CenterTitle();
04942 hNuMuBarRecoEn->SetFillColor(0);
04943 hNuMuBarRecoEn->SetLineColor(2);
04944
04945
04946 TH1F* hNuMuRecoEn=new TH1F("hNuMuRecoEn","hNuMuRecoEn",4*352,-32,320);
04947 hNuMuRecoEn->GetXaxis()->SetTitle("Energy (GeV)");
04948 hNuMuRecoEn->GetXaxis()->CenterTitle();
04949 hNuMuRecoEn->GetYaxis()->SetTitle("");
04950 hNuMuRecoEn->GetYaxis()->CenterTitle();
04951 hNuMuRecoEn->SetFillColor(0);
04952 hNuMuRecoEn->SetLineColor(1);
04953 hNuMuRecoEn->SetLineWidth(2);
04954
04955
04956 TH1F* hNuMuBarRecoY=new TH1F("hNuMuBarRecoY","hNuMuBarRecoY",
04957 4*144,-0.16,1.28);
04958 hNuMuBarRecoY->GetXaxis()->SetTitle("y");
04959 hNuMuBarRecoY->GetXaxis()->CenterTitle();
04960 hNuMuBarRecoY->GetYaxis()->SetTitle("");
04961 hNuMuBarRecoY->GetYaxis()->CenterTitle();
04962 hNuMuBarRecoY->SetFillColor(0);
04963 hNuMuBarRecoY->SetLineColor(2);
04964
04965
04966 TH1F* hNuMuRecoY=new TH1F("hNuMuRecoY","hNuMuRecoY",4*144,-0.16,1.28);
04967 hNuMuRecoY->GetXaxis()->SetTitle("y");
04968 hNuMuRecoY->GetXaxis()->CenterTitle();
04969 hNuMuRecoY->GetYaxis()->SetTitle("");
04970 hNuMuRecoY->GetYaxis()->CenterTitle();
04971 hNuMuRecoY->SetFillColor(0);
04972 hNuMuRecoY->SetLineColor(1);
04973 hNuMuRecoY->SetLineWidth(2);
04974
04975
04976 TH1F* hDpIDAll=new TH1F("hDpIDAll","hDpIDAll",4*160,-1.6,1.6);
04977 hDpIDAll->GetXaxis()->SetTitle("PID (from MadDpID)");
04978 hDpIDAll->GetXaxis()->CenterTitle();
04979 hDpIDAll->GetYaxis()->SetTitle("");
04980 hDpIDAll->GetYaxis()->CenterTitle();
04981 hDpIDAll->SetFillColor(0);
04982 hDpIDAll->SetLineColor(1);
04983 hDpIDAll->SetLineWidth(2);
04984
04985
04986 TH1F* hDpIDPQ=new TH1F("hDpIDPQ","hDpIDPQ",4*160,-1.6,1.6);
04987 hDpIDPQ->GetXaxis()->SetTitle("PID (from MadDpID)");
04988 hDpIDPQ->GetXaxis()->CenterTitle();
04989 hDpIDPQ->GetYaxis()->SetTitle("");
04990 hDpIDPQ->GetYaxis()->CenterTitle();
04991 hDpIDPQ->SetFillColor(0);
04992 hDpIDPQ->SetLineColor(1);
04993 hDpIDPQ->SetLineWidth(2);
04994
04995
04996 TH1F* hDpIDNQ=new TH1F("hDpIDNQ","hDpIDNQ",4*160,-1.6,1.6);
04997 hDpIDNQ->GetXaxis()->SetTitle("PID (from MadDpID)");
04998 hDpIDNQ->GetXaxis()->CenterTitle();
04999 hDpIDNQ->GetYaxis()->SetTitle("");
05000 hDpIDNQ->GetYaxis()->CenterTitle();
05001 hDpIDNQ->SetFillColor(0);
05002 hDpIDNQ->SetLineColor(1);
05003 hDpIDNQ->SetLineWidth(2);
05004
05005
05006 TH1F* hDpIDPQN_1=new TH1F("hDpIDPQN_1","hDpIDPQN_1",4*160,-1.6,1.6);
05007 hDpIDPQN_1->GetXaxis()->SetTitle("PID (from MadDpID)");
05008 hDpIDPQN_1->GetXaxis()->CenterTitle();
05009 hDpIDPQN_1->GetYaxis()->SetTitle("");
05010 hDpIDPQN_1->GetYaxis()->CenterTitle();
05011 hDpIDPQN_1->SetFillColor(0);
05012 hDpIDPQN_1->SetLineColor(1);
05013 hDpIDPQN_1->SetLineWidth(2);
05014
05015
05016 TH1F* hDpIDNQN_1=new TH1F("hDpIDNQN_1","hDpIDNQN_1",4*160,-1.6,1.6);
05017 hDpIDNQN_1->GetXaxis()->SetTitle("PID (from MadDpID)");
05018 hDpIDNQN_1->GetXaxis()->CenterTitle();
05019 hDpIDNQN_1->GetYaxis()->SetTitle("");
05020 hDpIDNQN_1->GetYaxis()->CenterTitle();
05021 hDpIDNQN_1->SetFillColor(0);
05022 hDpIDNQN_1->SetLineColor(1);
05023 hDpIDNQN_1->SetLineWidth(2);
05024
05025
05026
05027 Int_t entries=static_cast<Int_t>(this->GetEntries());
05028 Int_t startTimeSecs=2000000000;
05029 Int_t endTimeSecs=1;
05030 vector<Double_t> vPot;
05031 vPot.reserve(entries);
05032 vector<Double_t> vTime;
05033 vTime.reserve(entries);
05034 vector<Double_t> vPotBad;
05035 vPotBad.reserve(entries/10);
05036 vector<Double_t> vTimeBad;
05037 vTimeBad.reserve(entries/10);
05038
05039 vector<Double_t> vTimeNuMuEvt;
05040 vTimeNuMuEvt.reserve(entries/5);
05041 vector<Double_t> vTimeNuMuBarEvt;
05042 vTimeNuMuBarEvt.reserve(entries/5);
05043
05044
05045 map<Int_t,Double_t> deltaTs;
05046 map<Int_t,Int_t> evtsPerSlc;
05047
05048
05049 Int_t goodSpillCounter=0;
05050 Int_t badSpillCounter=0;
05051 Float_t totalPot=0;
05052 Float_t totalPotBad=0;
05053 Int_t evtCounter=0;
05054 Int_t evtNotlitime=0;
05055 Int_t evtNotIsLI=0;
05056 Int_t evtInFidVolCounter=0;
05057 Int_t evtWithTrkCounter=0;
05058 Int_t goodBestTrkCounter=0;
05059 Int_t goodTrkPassCounter=0;
05060 Int_t goodRecoEnCounter=0;
05061 Int_t goodUVDiffCounter=0;
05062 Int_t trkInFidVolCounter=0;
05063 Int_t goodFitSigQPCounter=0;
05064 Int_t goodFitProbCounter=0;
05065 Int_t goodFitChi2Counter=0;
05066 Int_t goodDeltaTCounter=0;
05067 Int_t goodEvtsPerSlcCounter=0;
05068 Int_t goodPIDCounter=0;
05069 Int_t goodBeyondTrkEndCounter=0;
05070 Int_t goodShwFractCounter=0;
05071 Int_t nuNQCounter=0;
05072 Int_t nuPQCounter=0;
05073
05074 string sTxtFilePrefix="nmb"+config.sBeamType;
05075 ofstream& nmbTxt=*(this->OpenTxtFile(config,
05076 sTxtFilePrefix.c_str()));
05077
05078
05079 const NuCuts cuts;
05080 const NuGeneral general;
05081 const NuPlots plots;
05082 const NuReco reco;
05083 const NuExtraction ext;
05084 const NuBeam beam;
05085
05086
05087
05088 NuMadAnalysis& mad=*new NuMadAnalysis();
05089 MadDpID& dp=*new MadDpID();
05090
05091 BeamType::BeamType_t beamType=static_cast
05092 <BeamType::BeamType_t>(config.beamType);
05093
05094
05095 NuCuts::NuAnaVersion_t anaVersion=NuCuts::kJJH1;
05096
05097
05098 const NuLibrary& lib=NuLibrary::Instance();
05099
05103
05104 cout<<endl
05105 <<"************************************************"<<endl
05106 <<"*** Starting main loop over snarls ***"<<endl
05107 <<"************************************************"<<endl;
05108
05109 this->InitialiseLoopVariables();
05110
05111
05112 for(Int_t entry=0;entry<this->GetEntries();entry++){
05113
05114 this->SetLoopVariables(entry);
05115
05116
05117 const NtpStRecord& ntp=(*this->GetNtpStRecord());
05118
05119 const RecCandHeader& rec=ntp.GetHeader();
05120 MAXMSG("NuAnalysis",Msg::kInfo,5)
05121 <<"Found: run="<<rec.GetRun()
05122 <<", subrun="<<rec.GetSubRun()
05123 <<", detector="<<rec.GetVldContext().GetDetector()
05124 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
05125 <<endl;
05126
05127
05128 plots.FillGeneralHistos(ntp);
05129
05130
05131 NuEvent nuSnarl;
05132
05133 nuSnarl.anaVersion=anaVersion;
05134
05135 ext.ExtractGeneralInfo(ntp,nuSnarl);
05136
05137
05138 plots.FillTrueFidEnergySpect(nuSnarl);
05139
05140 Int_t evTime=rec.GetVldContext().GetTimeStamp().GetSec();
05141
05142 if (evTime>endTimeSecs) endTimeSecs=evTime;
05143 if (evTime<startTimeSecs && evTime>0) startTimeSecs=evTime;
05144 if (evTime<=0) cout<<"Bad time="<<evTime<<endl;
05145
05146 general.EpochTo1995(evTime);
05147
05148 if (!beam.IsGoodSpillAndFillPot
05149 (this->GetNtpBDLiteRecord(),config,nuSnarl)) {
05150 badSpillCounter++;
05151 totalPotBad+=nuSnarl.pot;
05152 vPotBad.push_back(nuSnarl.pot*1e12);
05153 vTimeBad.push_back(evTime);
05154 continue;
05155 }
05156 goodSpillCounter++;
05157 totalPot+=nuSnarl.pot;
05158 vPot.push_back(nuSnarl.pot*1e12);
05159 vTime.push_back(evTime);
05160
05161 dp.ChoosePDFs(rec.GetVldContext().GetDetector(),beamType,
05162 ReleaseType::AsString(nuSnarl.releaseType),
05163 ReleaseType::AsString(nuSnarl.releaseType));
05164 mad.SetEntry(&ntp);
05165
05166 TClonesArray& evtTca=(*ntp.evt);
05167 const Int_t numEvts=evtTca.GetEntriesFast();
05168
05170
05172 for (Int_t ievt=0;ievt<numEvts;ievt++){
05173 const NtpSREvent& evt=
05174 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
05175 evtCounter++;
05176
05177 const NtpSREventSummary& evthdr=ntp.evthdr;
05178
05179 if (evthdr.litime!=-1) {
05180 MAXMSG("NuAnalysis",Msg::kInfo,500)
05181 <<"skipping event: litime="<<evthdr.litime<<endl;
05182 continue;
05183 }
05184 evtNotlitime++;
05185
05186 Bool_t isLI=LISieve::IsLI(ntp);
05187 if (isLI){
05188 MAXMSG("NuAnalysis",Msg::kInfo,50)
05189 <<"Found LI"
05190 <<", run="<<config.run
05191 <<", entry="<<entry
05192 <<", event="<<ievt
05193 <<", litime="<<evthdr.litime
05194 <<endl;
05195 }
05196 else {
05197 MAXMSG("NuAnalysis",Msg::kInfo,50)
05198 <<"Not LI"
05199 <<", run="<<config.run
05200 <<", entry="<<entry
05201 <<", event="<<ievt
05202 <<", litime="<<evthdr.litime
05203 <<endl;
05204 }
05205 if (isLI) continue;
05206 evtNotIsLI++;
05207
05208
05209 NuEvent nu;
05210
05211
05212 ext.ExtractGeneralInfo(ntp,nu);
05213
05214 lib.ext.ExtractEvtInfo(evt,nu);
05215
05216 lib.ext.ExtractTrkInfo(ntp,evt,nu);
05217
05218 lib.ext.ExtractShwInfo(ntp,evt,nu);
05219
05220 Bool_t isInFidVol=cuts.IsInFidVol(evt.vtx.x,evt.vtx.y,
05221 evt.vtx.z,
05222 evt.vtx.u,evt.vtx.v,0,0,
05223 nu.detector,nu.anaVersion,
05224 nu.releaseType,nu.simFlag);
05225
05226 if (!isInFidVol) continue;
05227 evtInFidVolCounter++;
05228
05229 nu.dpID=dp.CalcPID(&mad,ievt,0);
05230 MAXMSG("NuAnalysis",Msg::kDebug,50)
05231 <<"dpID="<<nu.dpID<<endl;
05232 hDpIDAll->Fill(nu.dpID);
05233
05234
05235
05236 if (evt.ntrack!=1) continue;
05237 evtWithTrkCounter++;
05238
05239
05240 Int_t bestTrack=lib.reco.GetBestTrack(nu);
05241 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
05242 (ntp,evt,bestTrack-1);
05243 if (ptrk==0) continue;
05244 const NtpSRTrack& trk=*ptrk;
05245 goodBestTrkCounter++;
05246
05247
05248 if (!trk.fit.pass) continue;
05249 goodTrkPassCounter++;
05250
05252
05253 lib.reco.GetEvtEnergy(nu);
05254 goodRecoEnCounter++;
05255
05256
05257 reco.GetTruthInfo(ntp,evt,nu);
05258
05259 if (!cuts.IsGoodUVVtx(nu)){
05260 continue;
05261 }
05262 goodUVDiffCounter++;
05263
05264
05265 Bool_t isTrkInFidVol=cuts.IsInFidVol(trk.vtx.x,trk.vtx.y,
05266 trk.vtx.z,
05267 trk.vtx.u,trk.vtx.v,0,0,
05268 nu.detector,nu.anaVersion,
05269 nu.releaseType,nu.simFlag);
05270 if (!isTrkInFidVol) continue;
05271 trkInFidVolCounter++;
05272
05273 plots.FillSigmaQPPlots(nu);
05274 plots.FillTrueDpIDHistos(nu);
05275
05276
05277 nu.charge=trk.momentum.qp<0 ? -1 : 1;
05278 if (trk.momentum.qp==0) nu.charge=0;
05279
05280
05281 if (!cuts.IsGoodSigmaQP_QP(nu)) continue;
05282 goodFitSigQPCounter++;
05283
05284
05285 if (!cuts.IsGoodFitProb(nu)) continue;
05286 goodFitProbCounter++;
05287
05288 Float_t chi2PerNdof=999999;
05289 if (!cuts.IsGoodFitChi2PerNdof(nu)){
05290 MAXMSG("NuAnalysis",Msg::kVerbose,100)
05291 <<"Cutting out bad chi2PerNdof="<<chi2PerNdof<<endl;
05292 continue;
05293 }
05294 goodFitChi2Counter++;
05295
05296
05297 deltaTs.clear();
05298 reco.GetEvtDeltaTs(ntp,deltaTs);
05299 evtsPerSlc.clear();
05300 reco.GetEvtsPerSlc(ntp,evtsPerSlc);
05301 MsgFormat ffmt("%9.9f");
05302 MAXMSG("NuReco",Msg::kDebug,100)
05303 <<"smallestDelta="<<ffmt(deltaTs[ievt])
05304 <<", evtsPerSlc="<<evtsPerSlc[evt.slc]<<endl;
05305
05306
05307 if (config.detector==Detector::kNear &&
05308 deltaTs[ievt]<50*Munits::ns) continue;
05309 goodDeltaTCounter++;
05310
05311
05312 if (config.detector==Detector::kNear &&
05313 evtsPerSlc[evt.slc]!=1) continue;
05314 goodEvtsPerSlcCounter++;
05315
05316 plots.FillTrueDpIDHistosPQNQ(nu);
05317 if (trk.momentum.qp<0) {
05318 hDpIDNQN_1->Fill(nu.dpID);
05319 }
05320 else if (trk.momentum.qp>0) {
05321 hDpIDPQN_1->Fill(nu.dpID);
05322 }
05323
05324
05325 if (!cuts.IsGoodDpID(nu)) continue;
05326 goodPIDCounter++;
05327
05328 Float_t sigCorBeyond=0;
05329 goodBeyondTrkEndCounter++;
05330
05331 Float_t fractFull=0;
05332
05333 goodShwFractCounter++;
05334
05335 Float_t evtTrkVtxDiff=evt.plane.beg-trk.plane.beg;
05336
05337
05338
05339 if (trk.momentum.qp<0) {
05340 hNuMuRecoEn->Fill(nu.energy);
05341 hNuMuRecoY->Fill(nu.y);
05342 hDpIDNQ->Fill(nu.dpID);
05343 vTimeNuMuEvt.push_back(evTime);
05344 nuNQCounter++;
05345 }
05346 else if (trk.momentum.qp>0) {
05347
05348 nmbTxt<<nu.run<<" "<<nu.snarl<<" "<<nu.evt<<endl;
05349 string sStop="PC";
05350 if (nu.containmentFlag==1 ||
05351 nu.containmentFlag==3) sStop="FC";
05352
05353 Float_t dRng=999999;
05354 if (nu.trkEnMC) dRng=(nu.trkEnRange-nu.trkEnMC)/nu.trkEnMC;
05355 nmbTxt<<"# "
05356 <<sStop<<nu.containmentFlag
05357 <<", dRng="<<dRng
05358 <<", beyond="<<sigCorBeyond
05359 <<", frFull="<<fractFull
05360 <<", e-tVtx="<<evtTrkVtxDiff
05361 <<endl;
05362 nmbTxt<<"# "<<endl;
05363
05364 hNuMuBarRecoEn->Fill(nu.energy);
05365 hNuMuBarRecoY->Fill(nu.y);
05366 hDpIDPQ->Fill(nu.dpID);
05367 vTimeNuMuBarEvt.push_back(evTime);
05368 nuPQCounter++;
05369 }
05370 else cout<<"ahhh, zero qp"<<endl;
05371
05372 plots.FillRecoEnYHistosN(nu);
05373 plots.FillDPIdSigmaQPPlotsN(nu);
05374
05375 plots.FillShwHistos(ntp,evt,nu);
05376 plots.FillXYZHistos(nu);
05377 plots.FillContainmentHistos(nu);
05378 plots.FillTruePIDHistos(nu);
05379 plots.FillTrackResponseHistos(ntp,trk,nu);
05380 }
05381 }
05382
05386
05387 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
05388
05389 MSG("NuAnalysis",Msg::kInfo)
05390 <<"Found start time and end time: "
05391 <<startTimeSecs<<" -> "<<endTimeSecs<<endl;
05392 general.EpochTo1995(endTimeSecs);
05393 general.EpochTo1995(startTimeSecs);
05394
05395 MSG("NuAnalysis",Msg::kInfo)<<"Filling vs. Time plots..."<<endl;
05396 TH1F* hSpillsVsTime=new TH1F("hSpillsVsTime","hSpillsVsTime",
05397 100,startTimeSecs,endTimeSecs);
05398 general.TH1FFill(hSpillsVsTime,vTime);
05399 general.SetGraphAxis(hSpillsVsTime->GetXaxis(),
05400 startTimeSecs,endTimeSecs);
05401 hSpillsVsTime->GetXaxis()->CenterTitle();
05402 hSpillsVsTime->GetYaxis()->SetTitle("Spills");
05403 hSpillsVsTime->GetYaxis()->CenterTitle();
05404
05405 TH1F* hSpillsBadVsTime=new TH1F("hSpillsBadVsTime","hSpillsBadVsTime",
05406 100,startTimeSecs,endTimeSecs);
05407 general.TH1FFill(hSpillsBadVsTime,vTimeBad);
05408 general.SetGraphAxis(hSpillsBadVsTime->GetXaxis(),
05409 startTimeSecs,endTimeSecs);
05410 hSpillsBadVsTime->GetXaxis()->CenterTitle();
05411 hSpillsBadVsTime->GetYaxis()->SetTitle("Spills");
05412 hSpillsBadVsTime->GetYaxis()->CenterTitle();
05413
05414 TH1F* hNuMuVsTime=new TH1F("hNuMuVsTime","hNuMuVsTime",
05415 100,startTimeSecs,endTimeSecs);
05416 general.TH1FFill(hNuMuVsTime,vTimeNuMuEvt);
05417 general.SetGraphAxis(hNuMuVsTime->GetXaxis(),
05418 startTimeSecs,endTimeSecs);
05419 hNuMuVsTime->GetXaxis()->CenterTitle();
05420 hNuMuVsTime->GetYaxis()->SetTitle("#nu_#mu events");
05421 hNuMuVsTime->GetYaxis()->CenterTitle();
05422
05423 TH1F* hNuMuBarVsTime=new TH1F("hNuMuBarVsTime","hNuMuBarVsTime",
05424 100,startTimeSecs,endTimeSecs);
05425 general.TH1FFill(hNuMuBarVsTime,vTimeNuMuBarEvt);
05426 general.SetGraphAxis(hNuMuBarVsTime->GetXaxis(),
05427 startTimeSecs,endTimeSecs);
05428 hNuMuBarVsTime->GetXaxis()->CenterTitle();
05429 hNuMuBarVsTime->GetYaxis()->SetTitle("#bar{#nu_#mu} events");
05430 hNuMuBarVsTime->GetYaxis()->CenterTitle();
05431
05432 TH1F* hPotVsTime=new TH1F("hPotVsTime","hPotVsTime",
05433 100,startTimeSecs,endTimeSecs);
05434 general.TH1FFill(hPotVsTime,vTime,vPot);
05435 general.SetGraphAxis(hPotVsTime->GetXaxis(),
05436 startTimeSecs,endTimeSecs);
05437 hPotVsTime->GetXaxis()->CenterTitle();
05438 hPotVsTime->GetYaxis()->SetTitle("POT");
05439 hPotVsTime->GetYaxis()->CenterTitle();
05440
05441 TH1F* hPotBadVsTime=new TH1F("hPotBadVsTime","hPotBadVsTime",
05442 100,startTimeSecs,endTimeSecs);
05443 general.TH1FFill(hPotBadVsTime,vTimeBad,vPotBad);
05444 general.SetGraphAxis(hPotBadVsTime->GetXaxis(),
05445 startTimeSecs,endTimeSecs);
05446 hPotBadVsTime->GetXaxis()->CenterTitle();
05447 hPotBadVsTime->GetYaxis()->SetTitle("POT");
05448 hPotBadVsTime->GetYaxis()->CenterTitle();
05449
05450 Float_t totalPotHist=0;
05451 Float_t totalPotBadHist=0;
05452 cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
05453 Float_t pcPotRejected=0;
05454 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
05455
05456
05457 cout<<"Number of spills ="<<this->GetEntries()<<endl
05458 <<"Number of good spills="<<goodSpillCounter<<endl
05459 <<"Number of bad spills ="<<badSpillCounter<<endl
05460 <<"Sum of good+bad spill="<<badSpillCounter+goodSpillCounter<<endl
05461 <<endl
05462 <<"Total POT good ="<<totalPot*1e12<<endl
05463 <<"Total POT bad ="<<totalPotBad*1e12<<endl
05464 <<"Total POT good (hist)="<<totalPotHist*1e12<<endl
05465 <<"Total POT bad (hist)="<<totalPotBadHist*1e12
05466 <<" ("<<pcPotRejected<<"%)"<<endl;
05467
05468 MSG("NuAnalysis",Msg::kInfo)
05469 <<endl<<"*************************************"<<endl
05470 <<"Total Events ="<<evtCounter<<endl
05471 <<"Evt not litime ="<<evtNotlitime<<endl
05472 <<"Evt not IsLI ="<<evtNotIsLI<<endl
05473 <<"Evts in Fid ="<<evtInFidVolCounter<<endl
05474 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
05475 <<"Good Best Trk ="<<goodBestTrkCounter<<endl
05476 <<"Track Pass ="<<goodTrkPassCounter<<endl
05477 <<"Good reco energy="<<goodRecoEnCounter<<endl
05478 <<"Good UVDiff ="<<goodUVDiffCounter<<endl
05479 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
05480 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
05481 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
05482 <<"Good Fit Chi2 ="<<goodFitChi2Counter<<endl
05483 <<"Good Delta T ="<<goodDeltaTCounter<<endl
05484 <<"Good EvtPerSlc ="<<goodEvtsPerSlcCounter<<endl
05485 <<"Good DP ID ="<<goodPIDCounter<<endl
05486 <<"Good beyond tEnd="<<goodBeyondTrkEndCounter<<endl
05487 <<"Good shw fract ="<<goodShwFractCounter<<endl
05488 <<"*************************************"<<endl;
05489
05490 Float_t pcPQ=0;
05491 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
05492 (nuPQCounter+nuNQCounter);
05493 MSG("NuAnalysis",Msg::kInfo)
05494 <<endl
05495 <<"Neutrinos with NQ="<<nuNQCounter<<endl
05496 <<"Neutrinos with PQ="<<nuPQCounter
05497 <<" ("<<pcPQ<<"%)"<<endl
05498 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
05499 <<endl;
05500
05501 MSG("NuAnalysis",Msg::kInfo)
05502 <<" ** Finished EnergySpect method **"<<endl;
05503 }
05504
05505
05506
05507 void NuAnalysis::CopyAcrossHistos(TDirectory* dirInput,
05508 TDirectory* dirOutput) const
05509 {
05510
05511 vector<string> vObjectNames;
05512 vObjectNames.push_back("hDetector");
05513 vObjectNames.push_back("hSimFlag");
05514 vObjectNames.push_back("hTrigSrc");
05515 vObjectNames.push_back("hRun");
05516 vObjectNames.push_back("hPottortgt");
05517 vObjectNames.push_back("hPotBadtortgt");
05518 vObjectNames.push_back("hPottrtgtd");
05519 vObjectNames.push_back("hPotBadtrtgtd");
05520 vObjectNames.push_back("hPottor101");
05521 vObjectNames.push_back("hPotBadtor101");
05522 vObjectNames.push_back("hPottr101d");
05523 vObjectNames.push_back("hPotBadtr101d");
05524 vObjectNames.push_back("hSpillsPerFile");
05525
05526 NuGeneral general;
05527 general.CopyAcrossObjects(dirInput,dirOutput,vObjectNames);
05528 }
05529
05530
05531
05532 void NuAnalysis::NMBSummaryTreeAna()
05533 {
05534 string sFilePrefix="NMBSumAna";
05535 fOutFile=this->OpenFile(100,sFilePrefix.c_str());
05536
05537 TDirectory* dirOutput=gDirectory;
05538 cout<<"After opening output file:"<<endl;
05539 dirOutput->Print();
05540
05541 string sBase="/home/hartnell/mytest/NuMuBar/";
05542 string inputFileName=
05543 "NMBChgSepL010185iN00009059_2.root";
05544
05545 inputFileName=sBase+inputFileName;
05546
05547 NuInputEvents* fpInput=new NuInputEvents();
05548 NuInputEvents& input=*fpInput;
05549 input.InputFileName(inputFileName);
05550 input.InitialiseChains();
05551 input.InitialiseNuEventBranch();
05552 TDirectory* dirInput=input.OpenInputFile();
05553
05554
05555 this->CopyAcrossHistos(dirInput,dirOutput);
05556
05557
05558
05559
05560
05561
05562
05563
05564 input.ResetNuEventLoopPositionToStart();
05565
05566 static NuGeneral general;
05567 static NuPlots plots;
05568 static NuCuts cuts;
05569
05570 for (Int_t i=0;i<input.GetEntriesNuEvent();++i) {
05571
05572 this->PrintLoopProgress(i,input.GetEntriesNuEvent(),1);
05573
05574 NuEvent& nu=const_cast<NuEvent&>(input.GetNextNuEvent(Msg::kDebug));
05575
05576
05577 plots.FillTrueDpIDHistos(nu);
05578 plots.FillTrueDpIDHistosPQNQ(nu);
05579 plots.FillDPIdSigmaQPPassPreSelCutPlots(nu);
05580 plots.FillSigmaQPPlots(nu);
05581 plots.FillTruePoIDHistos(nu);
05582 plots.FillTruePoIDHistosPQNQ(nu);
05583 plots.FillTrueAbIDHistos(nu);
05584 plots.FillTrueAbIDHistosPQNQ(nu);
05585
05586
05587 nu.anaVersion=NuCuts::kJJH1;
05588
05589
05590
05591 if (!cuts.IsGoodAbID(nu)) {
05592 plots.FillDPIdSigmaQPFailDpIDCutPlots(nu);
05593 continue;
05594 }
05595 plots.FillDPIdSigmaQPPassDpIDCutPlots(nu);
05596
05597 plots.FillEnergyBinHistos(nu);
05598
05599
05600 if (!cuts.IsGoodSigmaQP_QP(nu)) {
05601 plots.FillDPIdSigmaQPFailSigQPCutPlots(nu);
05602 continue;
05603 }
05604 plots.FillDPIdSigmaQPPassSigQPCutPlots(nu);
05605
05606
05607 if (!cuts.IsGoodFitChi2PerNdof(nu)) continue;
05608
05609
05610
05611
05612
05613 if (!cuts.IsGoodFitProb(nu)) {
05614 plots.FillDPIdSigmaQPFailProbCutPlots(nu);
05615 continue;
05616 }
05617
05618
05619
05620
05621
05622
05623
05624
05625
05626
05627
05628
05629 Bool_t weightForOsc=false;
05630 if (weightForOsc){
05631
05632 Double_t weight=nu.rw;
05633
05634 if (nu.iaction==1) {
05635 Float_t dm2=2.7e-3;
05636 Double_t oscWeight=general.OscWeight(dm2,1,
05637 nu.energyMC);
05638 weight*=oscWeight;
05639 }
05640 nu.rw=weight;
05641 }
05642
05643
05644
05645
05646
05647
05649
05650
05651
05652
05653
05654
05655
05656
05657
05658
05659
05660
05661
05662
05663
05664 plots.FillRecoEnYHistosN(nu);
05665 plots.FillDPIdSigmaQPPlotsN(nu);
05666
05667
05668 plots.FillXYZHistos(nu);
05669 plots.FillContainmentHistos(nu);
05670 plots.FillTruePIDHistos(nu);
05671
05672 plots.FillRangeCurvCompHistos(nu);
05673 plots.FillKinematicsHistos(nu);
05674
05675 MAXMSG("NuMinuit",Msg::kInfo,5)
05676 <<"NMBSummaryTreeAna: NuEvent: index="<<nu.index
05677 <<", energy="<<nu.energy<<", energyMC="<<nu.energyMC<<endl;
05678 }
05679
05680 }
05681
05682
05683
05684 void NuAnalysis::LIRejectionTest()
05685 {
05686 MSG("NuAnalysis",Msg::kInfo)
05687 <<" ** Running LIRejectionTest method... **"<<endl;
05688
05689 NuConfig config;
05690 this->ExtractConfig(config);
05691
05692
05693 string sFilePrefix="LIRejectionTest"+config.sBeamType;
05694 fOutFile=this->OpenFile(config,sFilePrefix.c_str());
05695
05696 string sTxt="litag"+config.sBeamType;
05697 ofstream& litagTxt=*(this->OpenTxtFile(config,sTxt.c_str()));
05698 sTxt="liSieve"+config.sBeamType;
05699 ofstream& liSieveTxt=*(this->OpenTxtFile(config,sTxt.c_str()));
05700
05701
05702 const NuCuts cuts;
05703 const NuGeneral general;
05704 const NuPlots plots;
05705 const NuReco reco;
05706 const NuExtraction ext;
05707 const NuBeam beam;
05708
05709
05710 NuCounter cnt;
05711 NuTime time;
05712
05713
05714 NuPIDInterface pid;
05715
05716
05717 const NuZBeamReweight zBeamReweight;
05718
05719 config.anaVersion=NuCuts::kCC0250Std;
05720 config.useGeneratorReweight=false;
05721 config.reweightVersion=SKZPWeightCalculator::kPiMinus;
05722 Int_t mcVersionOveride=ReleaseType::kDaikon;
05723
05724 MSG("NuAnalysis",Msg::kInfo)
05725 <<"Setting:"<<endl
05726 <<" anaVersion="
05727 <<cuts.AsString(static_cast<NuCuts::NuAnaVersion_t>
05728 (config.anaVersion))
05729 <<" ("<<config.anaVersion<<")"<<endl
05730 <<" useGeneratorReweight="<<config.useGeneratorReweight<<endl
05731 <<" reweightVersion="<<config.reweightVersion<<endl;
05732
05733
05734
05735 if (config.mcVersion==ReleaseType::kData) {
05736 MSG("NuAnalysis",Msg::kInfo)
05737 <<"Overiding mcVersion:"<<endl
05738 <<" new mcVersion="
05739 <<ReleaseType::AsString(mcVersionOveride)
05740 <<" ("<<mcVersionOveride<<")"
05741 <<", old mcVersion="<<ReleaseType::AsString(config.mcVersion)
05742 <<" ("<<config.mcVersion<<")"
05743 <<endl;
05744 config.mcVersion=mcVersionOveride;
05745 }
05746
05747
05748 const NuLibrary& lib=NuLibrary::Instance();
05749
05753
05754 this->InitialiseLoopVariables();
05755
05756
05757 for(Int_t entry=0;entry<this->GetEntries();entry++){
05758
05759 this->SetLoopVariables(entry);
05760
05761
05762 const NtpStRecord& ntp=(*this->GetNtpStRecord());
05763
05764 const RecCandHeader& rec=ntp.GetHeader();
05765 MAXMSG("NuAnalysis",Msg::kInfo,5)
05766 <<"Found: run="<<rec.GetRun()
05767 <<", subrun="<<rec.GetSubRun()
05768 <<", detector="<<rec.GetVldContext().GetDetector()
05769 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
05770 <<", snarl="<<rec.GetSnarl()
05771 <<endl;
05772
05773
05774 plots.FillGeneralHistos(ntp);
05775
05776
05777 NuEvent nuSnarl;
05778 nuSnarl.anaVersion=config.anaVersion;
05779 nuSnarl.mcVersion=config.mcVersion;
05780 nuSnarl.recoVersion=config.recoVersion;
05781 ext.ExtractGeneralInfo(ntp,nuSnarl);
05782
05783
05784 plots.FillTrueFidEnergySpect(nuSnarl);
05785
05786 Int_t evTime=rec.GetVldContext().GetTimeStamp().GetSec();
05787 Bool_t goodSpill=beam.IsGoodSpillAndFillPot
05788 (this->GetNtpBDLiteRecord(),config,nuSnarl);
05789 time.FillTime(evTime,nuSnarl.pot,goodSpill);
05790 if (!goodSpill) {
05791 cnt.badSpillCounter++;
05792 cnt.totalPotBad+=nuSnarl.pot;
05793 continue;
05794 }
05795 cnt.goodSpillCounter++;
05796 cnt.totalPot+=nuSnarl.pot;
05797
05798 TClonesArray& evtTca=(*ntp.evt);
05799 const Int_t numEvts=evtTca.GetEntriesFast();
05800
05802
05804 for (Int_t ievt=0;ievt<numEvts;ievt++){
05805 const NtpSREvent& evt=
05806 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
05807 cnt.evtCounter++;
05808
05809 MAXMSG("NuAnalysis",Msg::kInfo,5)
05810 <<"Found event, "<<ievt+1<<"/"<<numEvts<<endl;
05811
05812 const NtpSREventSummary& evthdr=ntp.evthdr;
05813
05814 if (evthdr.litime!=-1) {
05815 MAXMSG("NuAnalysis",Msg::kInfo,500)
05816 <<"skipping event: litime="<<evthdr.litime<<endl;
05817 continue;
05818 }
05819 cnt.evtNotlitime++;
05820
05821 Bool_t isLI=LISieve::IsLI(ntp);
05822 if (isLI){
05823 MAXMSG("NuAnalysis",Msg::kInfo,5)
05824 <<endl<<endl<<endl
05825 <<"Found LI"
05826 <<", run="<<config.run
05827 <<", entry="<<entry
05828 <<", event="<<ievt
05829 <<", litime="<<evthdr.litime
05830 <<endl<<endl<<endl;
05831 }
05832 else {
05833 MAXMSG("NuAnalysis",Msg::kInfo,50)
05834 <<"Not LI"
05835 <<", run="<<config.run
05836 <<", entry="<<entry
05837 <<", event="<<ievt
05838 <<", litime="<<evthdr.litime
05839 <<endl;
05840 }
05841
05842 if (!isLI) cnt.evtNotIsLI++;
05843
05844
05845 NuEvent nu;
05846 nu.entry=entry;
05847 this->CopyConfig(config,nu);
05848
05849
05850 ext.ExtractGeneralInfo(ntp,nu);
05851
05852 lib.ext.ExtractEvtInfo(evt,nu);
05853
05854 lib.ext.ExtractTrkInfo(ntp,evt,nu);
05855
05856 lib.ext.ExtractShwInfo(ntp,evt,nu);
05857
05858
05859 Int_t litag=reco.FDRCBoundary(nu);
05860 if (litag) {
05861 MAXMSG("NuAnalysis",Msg::kInfo,50)
05862 <<endl<<endl
05863 <<"litag="<<litag
05864 <<", run="<<config.run
05865 <<", entry="<<entry
05866 <<", event="<<ievt
05867 <<", litime="<<evthdr.litime
05868 <<endl
05869 <<", nu.nshw="<<nu.nshw
05870 <<", nu.rawPhEvt="<<nu.rawPhEvt
05871 <<", nu.planeEvtBeg="<<nu.planeEvtBeg
05872 <<", nu.planeEvtEnd="<<nu.planeEvtEnd
05873 <<endl<<endl<<endl;
05874 }
05875 if (litag>10) {
05876 MAXMSG("NuAnalysis",Msg::kInfo,50)
05877 <<endl<<endl
05878 <<"litag="<<litag
05879 <<", run="<<config.run
05880 <<", entry="<<entry
05881 <<", event="<<ievt
05882 <<", litime="<<evthdr.litime
05883 <<endl
05884 <<", nu.nshw="<<nu.nshw
05885 <<", nu.rawPhEvt="<<nu.rawPhEvt
05886 <<", nu.planeEvtBeg="<<nu.planeEvtBeg
05887 <<", nu.planeEvtEnd="<<nu.planeEvtEnd
05888 <<endl<<endl<<endl;
05889 }
05890
05891
05892 if (!(litag || isLI)) continue;
05893
05894
05895 reco.GetEvtEnergy(nu);
05896 cnt.goodRecoEnCounter++;
05897
05898
05899 pid.GetDpID(ntp,evt,nu);
05900 pid.GetAbID(ntp,evt,nu);
05901
05902 MAXMSG("NuAnalysis",Msg::kInfo,10)
05903 <<"After reco:"<<endl;
05904 MAXMSG("NuAnalysis",Msg::kInfo,10)
05905 <<endl<<plots.PrintEventInfo(cout,nu)<<endl;
05906
05907
05908 reco.GetTruthInfo(ntp,evt,nu);
05909
05910 if (litag) plots.PrintEventInfo(litagTxt,nu);
05911 if (isLI) plots.PrintEventInfo(liSieveTxt,nu);
05912
05913
05914 if (isLI) continue;
05915
05916
05917 if (!cuts.IsGoodNumberOfTracks(nu)) continue;
05918 cnt.evtWithTrkCounter++;
05919
05920
05921 zBeamReweight.ExtractZBeamReweight(ntp,evt,nu);
05922 reco.ApplyReweights(nu);
05923
05924 MAXMSG("NuAnalysis",Msg::kInfo,10)
05925 <<"PID: dp="<<nu.dpID<<", ab="<<nu.abID
05926 <<", ro="<<nu.roID<<", po="<<nu.poID
05927 <<", roNMB="<<nu.roIDNuMuBar<<endl;
05928
05929
05930 if (!cuts.IsInFidVolTrk(nu)) continue;
05931 cnt.trkInFidVolCounter++;
05932
05933
05934 if (!cuts.IsGoodTrackFitPass(nu)) continue;
05935 cnt.goodTrkPassCounter++;
05936
05937 if (!cuts.IsGoodDirCos(nu)) {
05938 MAXMSG("NuAnalysis",Msg::kInfo,50)
05939 <<"Bad DirCos, trkvtxdcosz="<<nu.trkvtxdcosz
05940 <<", dirCosNu="<<nu.dirCosNu<<endl;
05941 continue;
05942 }
05943 cnt.goodDirectionCosineCounter++;
05944
05945
05946 ext.ExtractDataQuality(ntp,nu);
05947 if (!cuts.IsGoodDataQuality(nu)) continue;
05948 cnt.goodDataQualityCounter++;
05949
05950
05951 ext.ExtractMinMaxEvtTimes(ntp,evt,nu);
05952 ext.ExtractTimeToNearestSpill(nu);
05953 MAXMSG("NuAnalysis",Msg::kInfo,1000)
05954 <<"Time to spill="
05955 <<(nu.timeEvtMin-nu.timeToNearestSpill)/(Munits::microsecond)
05956 <<" us"<<endl;
05957 if (!cuts.IsGoodTimeToNearestSpill(nu)) continue;
05958 cnt.goodTimeToNearestSpillCounter++;
05959
05960
05961 if (!cuts.IsGoodPID(nu)) {
05962 plots.FillDPIdSigmaQPFailDpIDCutPlots(nu);
05963 continue;
05964 }
05965 cnt.goodPIDCounter++;
05966 plots.FillDPIdSigmaQPPassDpIDCutPlots(nu);
05967
05968
05969 if (!cuts.IsGoodSigmaQP_QP(nu)) {
05970 plots.FillDPIdSigmaQPFailSigQPCutPlots(nu);
05971 continue;
05972 }
05973 cnt.goodFitSigQPCounter++;
05974 plots.FillDPIdSigmaQPPassSigQPCutPlots(nu);
05975
05976
05977 if (!cuts.IsGoodFitProb(nu)) {
05978 plots.FillDPIdSigmaQPFailProbCutPlots(nu);
05979 continue;
05980 }
05981 cnt.goodFitProbCounter++;
05982
05983 if (nu.charge==-1) cnt.nuNQCounter++;
05984 else if (nu.charge==+1) cnt.nuPQCounter++;
05985 else cout<<"ahhh, bad charge"<<endl;
05986
05987 plots.FillRecoEnYHistosN(nu);
05988 plots.FillDPIdSigmaQPPlotsN(nu);
05989 plots.FillEnergyBinHistos(nu);
05990
05991 plots.FillShwHistos(ntp,evt,nu);
05992 plots.FillXYZHistos(nu);
05993 plots.FillContainmentHistos(nu);
05994 plots.FillTruePIDHistos(nu);
05995
05996 plots.FillRangeCurvCompHistos(nu);
05997 plots.FillKinematicsHistos(nu);
05998
05999 }
06000 }
06001
06005
06006 MSG("NuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
06007
06008 cnt.PrintNMB();
06009
06010 MSG("NuAnalysis",Msg::kInfo)
06011 <<" ** Finished LIRejectionTest method **"<<endl;
06012 }
06013
06014
06015
06016 void NuAnalysis::CopyConfig(const NuConfig& config,NuEvent& nu)
06017 {
06018 nu.anaVersion=config.anaVersion;
06019
06020 nu.useGeneratorReweight=config.useGeneratorReweight;
06021 nu.sGeneratorConfigName=config.sGeneratorConfigName;
06022 nu.generatorConfigNo=config.generatorConfigNo;
06023 nu.reweightVersion=config.reweightVersion;
06024
06025 nu.mcVersion=config.mcVersion;
06026 nu.recoVersion=config.recoVersion;
06027 nu.releaseType=config.releaseType;
06028 nu.beamType=config.beamType;
06029
06030 nu.useDBForDataQuality=config.useDBForDataQuality;
06031 nu.useDBForSpillTiming=config.useDBForSpillTiming;
06032 nu.useDBForBeamInfo=config.useDBForBeamInfo;
06033
06034 nu.cutOnDataQuality=config.cutOnDataQuality;
06035 nu.cutOnSpillTiming=config.cutOnSpillTiming;
06036 nu.cutOnBeamInfo=config.cutOnBeamInfo;
06037
06038 nu.calcMajCurv=config.calcMajCurv;
06039 nu.calcRoID=config.calcRoID;
06040 }
06041
06042
06043
06044 void NuAnalysis::CopySnarlInfo(const NuEvent& nuSnarl,NuEvent& nu)
06045 {
06046
06047 nu.goodBeamSntp=nuSnarl.goodBeamSntp;
06048 nu.pot=nuSnarl.pot;
06049 nu.coilIsOk=nuSnarl.coilIsOk;
06050 nu.coilIsReverse=nuSnarl.coilIsReverse;
06051 nu.coilCurrent=nuSnarl.coilCurrent;
06052 }
06053
06054
06055
06056 void NuAnalysis::LoopOverTruthInfo(const NtpStRecord& ntp,
06057 NuOutputWriter& output,
06058 const NuEvent& nuSnarl) const
06059 {
06060 TClonesArray& mcTca=(*ntp.mc);
06061 const Int_t nummcs=mcTca.GetEntriesFast();
06062 static NuPlots plots;
06063 if (nummcs<=0) {
06064 MAXMSG("NuAnalysis",Msg::kDebug,1)
06065 <<"LoopOverTruthInfo: no MCTruth so return..."<<endl;
06066
06067 plots.FillTrueFidEnergySpect(nuSnarl);
06068 return;
06069 }
06070
06071
06072 const NuLibrary& lib=NuLibrary::Instance();
06073
06074
06075 for (Int_t imc=0;imc<nummcs;++imc){
06076 const NtpMCTruth& mc=
06077 *(dynamic_cast<NtpMCTruth*>(mcTca[imc]));
06078
06079
06080
06081 if (nuSnarl.detector==Detector::kNear &&
06082 (mc.vtxz<-0.1*(Munits::m) || mc.vtxz>7.3*(Munits::m))) {
06083 MAXMSG("NuAnalysis",Msg::kInfo,1)
06084 <<"Not looping over true events occuring in front or behind"
06085 <<" the ND calorimeter (-0.1<Z<7.3m)"
06086 <<" (this mc.vtxz="<<mc.vtxz<<")"<<endl;
06087 continue;
06088 }
06089
06090
06091
06092 NuEvent nu=nuSnarl;
06093
06094
06095 nu.mc=imc;
06096
06097
06098 lib.reco.GetTruthInfo(ntp,mc,nu);
06099
06100
06101
06102 nu.shwEn=nu.shwEnMC;
06103 nu.trkEn=fabs(nu.trkEnMC);
06104
06105
06106 lib.zBeamReweight.ExtractZBeamReweight(nu);
06107 lib.reco.ApplyReweights(nu);
06108
06109
06110 if (!lib.cuts.IsInFidVolLoose
06111 (nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,nu.detector)) continue;
06112
06113
06114 NuMCEvent& numc=output.GetNuMCEventToFill();
06115 lib.ext.NuMCEventFromNuEvent(nu,numc);
06116 numc.mc=imc;
06117
06118
06119 output.FillNuMCEventTree();
06120
06121
06122
06123
06124 Bool_t isInFidVol=lib.cuts.IsInFidVol
06125 (nu.vtxxMC,nu.vtxyMC,nu.vtxzMC,
06126 nu.vtxuMC,nu.vtxvMC,
06127 nu.planeTrkVtxMC,nu.rTrkVtxMC,
06128 nu.detector,nu.anaVersion,nu.releaseType,nu.simFlag);
06129 if (!isInFidVol) continue;
06130
06131
06132 plots.FillTrueFidEnergySpect(nu);
06133 }
06134 }
06135
06136
06137
06138 void NuAnalysis::ChargeSeparation()
06139 {
06140 MSG("NuAnalysis",Msg::kInfo)
06141 <<" ** Running ChargeSeparation method... **"<<endl;
06142
06146
06147 cout<<endl
06148 <<"************************************************"<<endl
06149 <<"*** Starting main loop over snarls ***"<<endl
06150 <<"************************************************"<<endl;
06151
06152 this->InitialiseLoopVariables();
06153
06154 for(Int_t entry=0;entry<this->GetEntries();entry++){
06155
06156 this->SetLoopVariables(entry);
06157
06158
06159 this->ChargeSeparationOneSnarl(this->GetNtpStRecord(),
06160 this->GetNtpBDLiteRecord(),0);
06161 }
06162
06163
06164 this->StoreOrFinishTree(NULL,NULL,1);
06165
06166 MSG("NuAnalysis",Msg::kInfo)
06167 <<" ** Finished ChargeSeparation method **"<<endl;
06168 }
06169
06170
06171
06172 void NuAnalysis::StoreOrFinishTree(NuOutputWriter* pOutput,
06173 NuTime* pTime,
06174 Bool_t finish) const
06175 {
06176
06177 static NuTime& time=*pTime;
06178 static NuOutputWriter& output=*pOutput;
06179 static Bool_t firstTime=true;
06180
06181 if (firstTime) {
06182 MSG("NuPIDInterface",Msg::kInfo)
06183 <<"StoreOrFinishTree: TDirectory for NuOutputWriter is:"<<endl;
06184 if (pOutput) output.GetDirectory()->pwd();
06185 else cout<<"Ahhh, no NuOutputWriter"<<endl;
06186
06187 MSG("NuPIDInterface",Msg::kInfo)
06188 <<"StoreOrFinishTree: gDirectory is:"<<endl;
06189 gDirectory->pwd();
06190 }
06191
06192
06193 if (finish){
06194 if (!firstTime) {
06195
06196 const NuLibrary& lib=NuLibrary::Instance();
06197
06198 MSG("NuPIDInterface",Msg::kInfo)
06199 <<"StoreOrFinishTree: gDirectory is:"<<endl;
06200 gDirectory->pwd();
06201 MSG("NuPIDInterface",Msg::kInfo)
06202 <<"After reseting gDirectory to NuOutputWriter location:"<<endl;
06203 gDirectory=output.GetDirectory();
06204 gDirectory->pwd();
06205
06206
06207 time.FillTimeHistos();
06208
06209
06210 lib.cnt.PrintFullDST();
06211 lib.cnt.PrintStdhepCounter();
06212
06213 lib.hist.CalcPOTsFromHistos();
06214
06215
06216 output.Finish();
06217 }
06218 else {
06219 MSG("NuAnalysis",Msg::kWarning)
06220 <<"StoreOrFinishTree was not run in store mode"<<endl;
06221 }
06222 }
06223
06224
06225 firstTime=false;
06226 }
06227
06228
06229
06230 void NuAnalysis::SetAnaFlags(NuConfig& config)
06231 {
06233
06234
06235 const NuLibrary& lib=NuLibrary::Instance();
06236
06237
06238 config.anaVersion=NuCuts::kFullDST;
06239
06240
06241
06242
06243 config.reweightVersion=SKZPWeightCalculator::kDetXs;
06244 config.useGeneratorReweight=false;
06245
06246
06247
06248
06249
06250
06251
06252
06253
06254
06255
06256
06257
06258
06259
06260 Int_t mcVersionOveride=ReleaseType::kDaikon;
06261
06262
06263 MSG("NuAnalysis",Msg::kInfo)
06264 <<"Setting:"<<endl
06265 <<" anaVersion="
06266 <<lib.cuts.AsString(static_cast<NuCuts::NuAnaVersion_t>
06267 (config.anaVersion))
06268 <<" ("<<config.anaVersion<<")"<<endl
06269 <<" useGeneratorReweight="<<config.useGeneratorReweight<<endl
06270 <<" reweightVersion="<<config.reweightVersion<<endl;
06271
06272
06273
06274 if (config.mcVersion==ReleaseType::kData) {
06275 MSG("NuAnalysis",Msg::kInfo)
06276 <<"Overiding mcVersion for PDF selection:"<<endl
06277 <<" new mcVersion="
06278 <<ReleaseType::AsString(mcVersionOveride)
06279 <<" ("<<mcVersionOveride<<")"
06280 <<", old mcVersion="<<ReleaseType::AsString(config.mcVersion)
06281 <<" ("<<config.mcVersion<<")"
06282 <<endl;
06283 config.mcVersion=mcVersionOveride;
06284 }
06285 }
06286
06287
06288
06289 void NuAnalysis::DoExtractions(const NtpStRecord& ntp,
06290 const NtpSREvent& evt,
06291 const NtpFitSARecord* pntpSA,
06292 NuPIDInterface& pid,NuEvent& nu) const
06293 {
06294
06295 const NuLibrary& lib=NuLibrary::Instance();
06296
06297
06298 lib.ext.ExtractAuxiliaryInfo(ntp,evt,pntpSA,nu);
06299
06300 lib.reco.GetTruthInfo(ntp,evt,nu);
06301
06302
06303 this->ExtractPIDsAndWeights(ntp,evt,pid,nu);
06304 }
06305
06306
06307
06308 void NuAnalysis::ExtractPIDsAndWeights(const NtpStRecord& ntp,
06309 const NtpSREvent& evt,
06310 NuPIDInterface& pid,
06311 NuEvent& nu) const
06312 {
06313
06314 const NuLibrary& lib=NuLibrary::Instance();
06315
06316
06317 lib.zBeamReweight.ExtractZBeamReweight(nu);
06318 lib.zBeamReweight.CalcGeneratorReweight(ntp,evt,nu);
06319 lib.reco.ApplyReweights(nu);
06320
06321
06322 pid.GetDpID(ntp,evt,nu);
06323 pid.GetAbID(ntp,evt,nu);
06324 if (nu.calcRoID) pid.GetRoID(ntp,evt,nu);
06325 else pid.GetRoIDNuMuBar(ntp,evt,nu);
06326 pid.GetPoID(nu);
06327 lib.ext.ExtractMajorityCurvature(ntp,evt,nu);
06328 }
06329
06330
06331
06332 void NuAnalysis::EndJob()
06333 {
06334 this->StoreOrFinishTree(NULL,NULL,1);
06335
06336
06337
06338 cout<<"fOutFile="<<fOutFile<<endl;
06339 this->WriteOutHistos();
06340 }
06341
06342
06343
06344 void NuAnalysis::ChargeSeparationOneSnarl(const NtpStRecord* pntp,
06345 const NtpBDLiteRecord* pntpBD,
06346 const NtpFitSARecord* pntpSA)
06347 {
06348 MAXMSG("NuAnalysis",Msg::kError,100)
06349 <<"ChargeSeparationOneSnarl is deprecated. Use MakeFullDst instead"
06350 <<endl;
06351
06354
06355
06356 NuLibrary& lib=NuLibrary::Instance();
06357
06358
06359 lib.cnt.entry++;
06360
06364
06365
06366
06367
06368
06369
06370 static NuConfig config;
06371 static Bool_t configNotExtracted=true;
06372 if (configNotExtracted) {
06373 configNotExtracted=!this->ExtractConfig(pntp,pntpBD,config);
06374 if (configNotExtracted) return;
06375 }
06376
06377 static Bool_t firstTime=true;
06378 static NuOutputWriter output;
06379
06380 static ofstream* pnmbTxt=0;
06381 static ofstream* pnmTxt=0;
06382 static ofstream* pnmbTxtProb=0;
06383 static ofstream* pnmbTxtProb001=0;
06384 static ofstream* pnmbTxtProb003=0;
06385
06386
06387 static const NuPlots plots;
06388 static const NuBeam beam;
06389 static NuTime time;
06390
06391
06392 static NuPIDInterface pid;
06393
06394 if (firstTime) {
06395 firstTime=false;
06396
06397
06398 string sFilePrefix="NuDST"+config.sBeamType;
06399
06400 output.SetupFile(config,sFilePrefix);
06401
06402 string sTxt="nmb"+config.sBeamType;
06403 pnmbTxt=this->OpenTxtFile(config,sTxt.c_str());
06404 string sTxtNM="nm"+config.sBeamType;
06405 pnmTxt=this->OpenTxtFile(config,sTxtNM.c_str());
06406 string sTxtProb="nmbProb"+config.sBeamType;
06407 pnmbTxtProb=this->OpenTxtFile(config,sTxtProb.c_str());
06408 string sTxtProb001="nmbProb001"+config.sBeamType;
06409 pnmbTxtProb001=this->OpenTxtFile(config,sTxtProb001.c_str());
06410 string sTxtProb003="nmbProb003"+config.sBeamType;
06411 pnmbTxtProb003=this->OpenTxtFile(config,sTxtProb003.c_str());
06412
06413
06414 this->SetAnaFlags(config);
06415 config.anaVersion=NuCuts::kCC0250Std;
06416
06417
06418
06419 this->StoreOrFinishTree(&output,&time,false);
06420
06421 cout<<endl
06422 <<"************************************************"<<endl
06423 <<"*** Starting main loop over snarls ***"<<endl
06424 <<"************************************************"<<endl;
06425 }
06426
06427
06428 static ofstream& nmbTxt=*pnmbTxt;
06429 static ofstream& nmTxt=*pnmTxt;
06430 static ofstream& nmbTxtProb=*pnmbTxtProb;
06431 static ofstream& nmbTxtProb001=*pnmbTxtProb001;
06432 static ofstream& nmbTxtProb003=*pnmbTxtProb003;
06433
06434
06435
06436 const NtpStRecord& ntp=*pntp;
06437
06441
06442 const RecCandHeader& rec=ntp.GetHeader();
06443 MAXMSG("NuAnalysis",Msg::kInfo,5)
06444 <<"Found: run="<<rec.GetRun()
06445 <<", subrun="<<rec.GetSubRun()
06446 <<", detector="<<rec.GetVldContext().GetDetector()
06447 <<", simFlag="<<rec.GetVldContext().GetSimFlag()
06448 <<", snarl="<<rec.GetSnarl()
06449 <<endl;
06450
06451
06452 plots.FillGeneralHistos(ntp);
06453
06454
06455 NuEvent nuSnarl;
06456 this->CopyConfig(config,nuSnarl);
06457
06458
06459 lib.ext.ExtractGeneralInfo(ntp,nuSnarl);
06460
06461
06462 this->LoopOverTruthInfo(ntp,output,nuSnarl);
06463
06464
06465
06466 lib.ext.ExtractCoilInfo(nuSnarl);
06467 plots.FillCoilCurrentHistos(nuSnarl);
06468 if (!nuSnarl.coilIsOk) return;
06469 lib.cnt.goodCoilCurrentSpillCounter++;
06470
06471
06472 beam.IsGoodSpillAndFillPot(pntpBD,config,nuSnarl);
06473 time.FillTime(nuSnarl.timeSec,nuSnarl.pot,nuSnarl.goodBeamSntp);
06474 if (!nuSnarl.goodBeamSntp) {
06475 lib.cnt.badSpillCounter++;
06476 lib.cnt.totalPotBad+=nuSnarl.pot;
06477 return;
06478 }
06479 lib.cnt.goodSpillCounter++;
06480 lib.cnt.totalPot+=nuSnarl.pot;
06481
06482
06483 static Float_t potSinceLastEvt=0;
06484 static Float_t potSinceLastEvtGood=0;
06485 static Float_t potSinceLastEvtBad=0;
06486 potSinceLastEvt+=nuSnarl.pot;
06487 if (nuSnarl.goodBeamSntp &&
06488 nuSnarl.coilIsOk) potSinceLastEvtGood+=nuSnarl.pot;
06489 else potSinceLastEvtBad+=nuSnarl.pot;
06490
06491
06492 TClonesArray& evtTca=(*ntp.evt);
06493 const Int_t numEvts=evtTca.GetEntriesFast();
06494
06496
06498 for (Int_t ievt=0;ievt<numEvts;ievt++){
06499 const NtpSREvent& evt=
06500 *dynamic_cast<NtpSREvent*>(evtTca[ievt]);
06501 lib.cnt.evtCounter++;
06502
06503 MAXMSG("NuAnalysis",Msg::kInfo,5)
06504 <<"Found event, "<<ievt+1<<"/"<<numEvts<<endl;
06505
06506
06507 if (evt.ntrack<1) continue;
06508 lib.cnt.evtWithTrkCounter++;
06509
06510
06511 const NtpSREventSummary& evthdr=ntp.evthdr;
06512 if (evthdr.litime!=-1) continue;
06513 lib.cnt.evtNotlitime++;
06514
06515
06516 NuEvent& nu=output.GetNuEventToFill();
06517 nu.entry=lib.cnt.entry;
06518 nu.nevt=numEvts;
06519
06520 this->CopySnarlInfo(nuSnarl,nu);
06521 this->CopyConfig(config,nu);
06522
06523
06524 lib.ext.ExtractBasicInfo(ntp,evt,nu);
06525
06526
06527 if (!lib.cuts.IsInFidVolLoose(nu)) continue;
06528 lib.cnt.evtInFidVolCounter++;
06529
06530
06531 lib.reco.GetEvtEnergy(nu);
06532
06533 MAXMSG("NuAnalysis",Msg::kDebug,10)
06534 <<endl<<"After reco:"<<endl
06535 <<plots.PrintEventInfo(cout,nu)<<endl;
06536
06537
06538 this->DoExtractions(ntp,evt,pntpSA,pid,nu);
06539
06540 MAXMSG("NuAnalysis",Msg::kInfo,10)
06541 <<"PID: dp="<<nu.dpID<<", ab="<<nu.abID<<", ro="<<nu.roID
06542 <<", po="<<nu.poID<<", roNMB="<<nu.roIDNuMuBar<<endl;
06543
06544
06545
06546 plots.FillNtupleEarliestLatestTime(nu);
06547
06548
06549
06550
06551 nu.potSinceLastEvt=potSinceLastEvt;
06552 nu.potSinceLastEvtGood=potSinceLastEvtGood;
06553 nu.potSinceLastEvtBad=potSinceLastEvtBad;
06554 potSinceLastEvt=0;
06555 potSinceLastEvtGood=0;
06556 potSinceLastEvtBad=0;
06557
06559
06560 output.FillNuEventTree();
06561
06562
06563 Int_t bestTrack=lib.reco.GetBestTrack(nu);
06564 const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX
06565 (ntp,evt,bestTrack-1);
06566 const NtpSRTrack& trk=*ptrk;
06567
06568
06569 if (!lib.cuts.IsGoodNumberOfTracks(nu)) continue;
06570 lib.cnt.evtWithTrkCounter++;
06571
06572
06573 if (!lib.cuts.IsInFidVolTrk(nu)) continue;
06574 lib.cnt.trkInFidVolCounter++;
06575
06576
06577 if (lib.cuts.IsLI(nu)) continue;
06578 lib.cnt.evtNotIsLI++;
06579
06580
06581 if (!lib.cuts.IsGoodDataQuality(nu)) continue;
06582 lib.cnt.goodDataQualityCounter++;
06583
06584
06585 plots.FillEvtAndSpillTimingPlots(nu);
06586 if (!lib.cuts.IsGoodTimeToNearestSpill(nu)) continue;
06587 lib.cnt.goodTimeToNearestSpillCounter++;
06588
06589
06590 if (!lib.