00001
00002
00003
00004
00005
00007
00008 #include "TClonesArray.h"
00009 #include "TObject.h"
00010
00011 #include "MessageService/MsgFormat.h"
00012 #include "MessageService/MsgService.h"
00013 #include "MCNtuple/NtpMCStdHep.h"
00014 #include "MCNtuple/NtpMCTruth.h"
00015 #include "CandNtupleSR/NtpSREvent.h"
00016 #include "CandNtupleSR/NtpSRShower.h"
00017 #include "CandNtupleSR/NtpSRTrack.h"
00018 #include "StandardNtuple/NtpStRecord.h"
00019 #include "TruthHelperNtuple/NtpTHEvent.h"
00020 #include "TruthHelperNtuple/NtpTHShower.h"
00021 #include "TruthHelperNtuple/NtpTHTrack.h"
00022
00023 #include "NtupleUtils/NuCounter.h"
00024 #include "NtupleUtils/NuEvent.h"
00025 #include "NtupleUtils/NuGeneral.h"
00026 #include "NtupleUtils/NuLibrary.h"
00027
00028 using std::endl;
00029 using std::cout;
00030 using std::map;
00031
00032 CVSID("$Id: NuCounter.cxx,v 1.7 2009/09/04 17:42:34 lefeuvre Exp $");
00033
00034
00035
00036 NuCounter::NuCounter()
00037 {
00038 MSG("NuCounter",Msg::kDebug)
00039 <<"Running NuCounter Constructor..."<<endl;
00040
00041 entry=-1;
00042
00043 goodCoilCurrentSpillCounter=0;
00044 goodSpillCounter=0;
00045 badSpillCounter=0;
00046 totalPot=0;
00047 totalPotDB=0;
00048 totalPotBad=0;
00049 evtCounter=0;
00050 evtNotlitime=0;
00051 evtNotIsLI=0;
00052 goodBeamDetPOTCountingStage=0;
00053 evtInFidVolCounter=0;
00054 evtWithTrkCounter=0;
00055 goodBestTrkCounter=0;
00056 goodTrkPassCounter=0;
00057 goodDirectionCosineCounter=0;
00058 goodDataQualityCounter=0;
00059 goodTimeToNearestSpillCounter=0;
00060 goodBeamInfoDBCounter=0;
00061 goodRecoEnCounter=0;
00062 goodUVDiffCounter=0;
00063 trkInFidVolCounter=0;
00064 goodFitSigQPCounter=0;
00065 goodFitProbCounter=0;
00066 goodFitChi2Counter=0;
00067 goodDeltaTCounter=0;
00068 goodEvtsPerSlcCounter=0;
00069 goodPIDCounter=0;
00070 goodBeyondTrkEndCounter=0;
00071 goodShwFractCounter=0;
00072 nuNQCounter=0;
00073 nuPQCounter=0;
00074 nuNoQCounter=0;
00075 nuNoTrkCounter=0;
00076
00077 MSG("NuCounter",Msg::kDebug)
00078 <<"Finished NuCounter Constructor"<<endl;
00079 }
00080
00081
00082
00083 NuCounter::~NuCounter()
00084 {
00085 MSG("NuCounter",Msg::kDebug)
00086 <<"Running NuCounter Destructor..."<<endl;
00087
00088
00089 MSG("NuCounter",Msg::kDebug)
00090 <<"Finished NuCounter Destructor"<<endl;
00091 }
00092
00093
00094
00095 void NuCounter::Print() const
00096 {
00097
00098 NuLibrary& lib=NuLibrary::Instance();
00099
00100 Float_t totalPotHist=0;
00101 Float_t totalPotBadHist=0;
00102 lib.cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
00103 Float_t pcPotRejected=0;
00104 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
00105
00106 cout<<"Number of good spills="<<goodSpillCounter<<endl
00107 <<"Number of bad spills ="<<badSpillCounter<<endl
00108 <<"Sum of good+bad spill="<<badSpillCounter+goodSpillCounter<<endl
00109 <<endl
00110 <<"Total POT good ="<<totalPot*1e12<<endl
00111 <<"Total POT bad ="<<totalPotBad*1e12<<endl
00112 <<"Total POT good (hist)="<<totalPotHist*1e12<<endl
00113 <<"Total POT bad (hist)="<<totalPotBadHist*1e12
00114 <<" ("<<pcPotRejected<<"%)"<<endl;
00115
00116 MSG("NuAnalysis",Msg::kInfo)
00117 <<endl<<"*************************************"<<endl
00118 <<"Total Events ="<<evtCounter<<endl
00119 <<"Evt not litime ="<<evtNotlitime<<endl
00120 <<"Evt not IsLI ="<<evtNotIsLI<<endl
00121 <<"Evts in Fid ="<<evtInFidVolCounter<<endl
00122 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
00123 <<"Good Best Trk ="<<goodBestTrkCounter<<endl
00124 <<"Track Pass ="<<goodTrkPassCounter<<endl
00125 <<"Good Dir. Cos. ="<<goodDirectionCosineCounter<<endl
00126 <<"Good UVDiff ="<<goodUVDiffCounter<<endl
00127 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
00128 <<"Good reco energy="<<goodRecoEnCounter<<endl
00129 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
00130 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
00131 <<"Good Fit Chi2 ="<<goodFitChi2Counter<<endl
00132 <<"Good Delta T ="<<goodDeltaTCounter<<endl
00133 <<"Good EvtPerSlc ="<<goodEvtsPerSlcCounter<<endl
00134 <<"Good PID ="<<goodPIDCounter<<endl
00135 <<"Good beyond tEnd="<<goodBeyondTrkEndCounter<<endl
00136 <<"Good shw fract ="<<goodShwFractCounter<<endl
00137 <<"*************************************"<<endl;
00138
00139 Float_t pcPQ=0;
00140 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
00141 (nuPQCounter+nuNQCounter);
00142 MSG("NuAnalysis",Msg::kInfo)
00143 <<endl
00144 <<"Neutrinos with NQ="<<nuNQCounter<<endl
00145 <<"Neutrinos with PQ="<<nuPQCounter
00146 <<" ("<<pcPQ<<"%)"<<endl
00147 <<"Neutrinos with undefined charge (q=0)="<<nuNoQCounter<<endl
00148 <<"Events with no track="<<nuNoTrkCounter<<endl
00149 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
00150 <<"Sum NQ+PQ+Q0+NoTrk="<<nuNQCounter+nuPQCounter+nuNoQCounter+nuNoTrkCounter<<endl
00151 <<endl;
00152
00153 }
00154
00155
00156
00157 void NuCounter::PrintNMB() const
00158 {
00159
00160 const NuLibrary& lib=NuLibrary::Instance();
00161
00162 Float_t totalPotHist=0;
00163 Float_t totalPotBadHist=0;
00164 lib.cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
00165 Float_t pcPotRejected=0;
00166 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
00167
00168 MSG("NuAnalysis",Msg::kInfo)
00169 <<endl<<"*************************************"<<endl
00170 <<"POT good beam ="<<totalPot*1e12<<endl
00171 <<"POT bad beam ="<<totalPotBad*1e12<<endl
00172 <<"POT good (hist)="<<totalPotHist*1e12<<endl
00173 <<"POT bad (hist)="<<totalPotBadHist*1e12
00174 <<" ("<<pcPotRejected<<"%)"<<endl;
00175
00176 MSG("NuAnalysis",Msg::kInfo)
00177 <<endl<<"*************************************"<<endl
00178 <<"Snarls processed ="<<entry+1<<endl
00179 <<"Good coil current ="<<goodCoilCurrentSpillCounter<<endl
00180 <<"Good beam spills ="<<goodSpillCounter<<endl
00181 <<"Bad beam spills ="<<badSpillCounter<<endl
00182 <<"Sum good+bad spills="<<badSpillCounter+goodSpillCounter
00183 <<endl;
00184
00185 MSG("NuAnalysis",Msg::kInfo)
00186 <<endl<<"*************************************"<<endl
00187 <<"Total Events ="<<evtCounter<<endl
00188 <<"Evt not litime ="<<evtNotlitime<<endl
00189 <<"Good BeamDet@POT="<<goodBeamDetPOTCountingStage<<endl
00190 <<"Evt Vtx in Fid ="<<evtInFidVolCounter<<endl
00191 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
00192 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
00193 <<"Evt not IsLI ="<<evtNotIsLI<<endl
00194 <<"Good Data Qual. ="<<goodDataQualityCounter<<endl
00195 <<"Good Time2Spill ="<<goodTimeToNearestSpillCounter<<endl
00196 <<"Good Beam DBInfo="<<goodBeamInfoDBCounter<<endl
00197 <<"Good TrkFitPass ="<<goodTrkPassCounter<<endl
00198 <<"Good Dir. Cos. ="<<goodDirectionCosineCounter<<endl
00199 <<"Good PID ="<<goodPIDCounter<<endl
00200 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
00201 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
00202 <<"*************************************"<<endl;
00203
00204 Float_t pcPQ=0;
00205 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
00206 (nuPQCounter+nuNQCounter);
00207 MSG("NuAnalysis",Msg::kInfo)
00208 <<endl
00209 <<"Neutrinos with NQ="<<nuNQCounter<<endl
00210 <<"Neutrinos with PQ="<<nuPQCounter
00211 <<" ("<<pcPQ<<"%)"<<endl
00212 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
00213 <<endl;
00214
00215 }
00216
00217
00218
00219 void NuCounter::PrintFullDST() const
00220 {
00221
00222 const NuLibrary& lib=NuLibrary::Instance();
00223
00224 Float_t totalPotHist=0;
00225 Float_t totalPotBadHist=0;
00226 lib.cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
00227 Float_t pcPotRejected=0;
00228 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
00229
00230 MSG("NuAnalysis",Msg::kInfo)
00231 <<endl<<"*************************************"<<endl
00232 <<"POT good beam ="<<totalPot*1e12<<endl
00233 <<"POT bad beam ="<<totalPotBad*1e12<<endl
00234 <<"POT good (hist)="<<totalPotHist*1e12<<endl
00235 <<"POT bad (hist)="<<totalPotBadHist*1e12
00236 <<" ("<<pcPotRejected<<"%)"<<endl;
00237
00238 MSG("NuAnalysis",Msg::kInfo)
00239 <<endl<<"*************************************"<<endl
00240 <<"Snarls processed ="<<entry+1<<endl
00241 <<"Good coil current ="<<goodCoilCurrentSpillCounter<<endl
00242 <<"Good beam spills ="<<goodSpillCounter<<endl
00243 <<"Bad beam spills ="<<badSpillCounter<<endl
00244 <<"Sum good+bad spills="<<badSpillCounter+goodSpillCounter
00245 <<endl;
00246
00247 MSG("NuAnalysis",Msg::kInfo)
00248 <<endl<<"*************************************"<<endl
00249 <<"Total Events ="<<evtCounter<<endl
00250 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
00251 <<"Evt not litime ="<<evtNotlitime<<endl
00252 <<"Evt Vtx in Fid ="<<evtInFidVolCounter<<endl
00253 <<"*************************************"<<endl;
00254
00255 Float_t pcPQ=0;
00256 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
00257 (nuPQCounter+nuNQCounter);
00258 MSG("NuAnalysis",Msg::kInfo)
00259 <<endl
00260 <<"Neutrinos with NQ="<<nuNQCounter<<endl
00261 <<"Neutrinos with PQ="<<nuPQCounter
00262 <<" ("<<pcPQ<<"%)"<<endl
00263 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
00264 <<endl;
00265 }
00266
00267
00268
00269 void NuCounter::PrintMicroDST() const
00270 {
00271
00272 const NuLibrary& lib=NuLibrary::Instance();
00273
00274 Float_t totalPotHist=0;
00275 Float_t totalPotBadHist=0;
00276 lib.cuts.CalcTotalPot(totalPotHist,totalPotBadHist);
00277 Float_t pcPotRejected=0;
00278 if (totalPotHist) pcPotRejected=100.*totalPotBadHist/totalPotHist;
00279
00280 MSG("NuAnalysis",Msg::kInfo)
00281 <<endl<<"*************************************"<<endl
00282 <<"POT good beam ="<<totalPot*1e12<<endl
00283 <<"POT bad beam ="<<totalPotBad*1e12<<endl
00284 <<"POT good (hist)="<<totalPotHist*1e12<<endl
00285 <<"POT bad (hist)="<<totalPotBadHist*1e12<<endl
00286 <<"POT good (DB) ="<<totalPotDB*1e12<<endl
00287 <<" ("<<pcPotRejected<<"%)"<<endl;
00288
00289 MSG("NuAnalysis",Msg::kInfo)
00290 <<endl<<"*************************************"<<endl
00291 <<"Total Events ="<<evtCounter<<endl
00292 <<"Good BeamDet@POT="<<goodBeamDetPOTCountingStage<<endl
00293 <<"Evt Vtx in Fid ="<<evtInFidVolCounter<<endl
00294 <<"Evts w/ Trk ="<<evtWithTrkCounter<<endl
00295 <<"Trk Vtx in Fid ="<<trkInFidVolCounter<<endl
00296 <<"Evt not IsLI ="<<evtNotIsLI<<endl
00297 <<"Good Data Qual. ="<<goodDataQualityCounter<<endl
00298 <<"Good Time2Spill ="<<goodTimeToNearestSpillCounter<<endl
00299 <<"Good Beam DBInfo="<<goodBeamInfoDBCounter<<endl
00300 <<"Good TrkFitPass ="<<goodTrkPassCounter<<endl
00301 <<"Good Dir. Cos. ="<<goodDirectionCosineCounter<<endl
00302 <<"Good PID ="<<goodPIDCounter<<endl
00303 <<"Good Fit SigQP ="<<goodFitSigQPCounter<<endl
00304 <<"Good Fit Prob ="<<goodFitProbCounter<<endl
00305 <<"*************************************"<<endl;
00306
00307 Float_t pcPQ=0;
00308 if (nuPQCounter+nuNQCounter>0) pcPQ=100.*nuPQCounter/
00309 (nuPQCounter+nuNQCounter);
00310 MSG("NuAnalysis",Msg::kInfo)
00311 <<endl
00312 <<"Neutrinos with NQ="<<nuNQCounter<<endl
00313 <<"Neutrinos with PQ="<<nuPQCounter
00314 <<" ("<<pcPQ<<"%)"<<endl
00315 <<"Sum NQ+PQ ="<<nuNQCounter+nuPQCounter<<endl
00316 <<endl;
00317 }
00318
00319
00320
00321 void NuCounter::CountTrkStdhepId(const NtpStRecord& ntp,
00322 const NtpSREvent& evt,
00323 const NtpSRTrack& trk,
00324 const NuEvent& nu)
00325 {
00326 TClonesArray& thtrkTca=(*ntp.thtrk);
00327 const Int_t numthtrks=thtrkTca.GetEntriesFast();
00328 if (numthtrks<=0) {
00329 return;
00330 }
00331
00332
00333 NuLibrary& lib=NuLibrary::Instance();
00334
00335 Int_t evtmc=-1;
00336 Int_t trkmc=-1;
00337 Int_t shwmc=-1;
00338
00339 TClonesArray& thevtTca=(*ntp.thevt);
00340 const Int_t numthevts=thevtTca.GetEntriesFast();
00341 if (numthevts>0) {
00342 const NtpTHEvent& thevt=
00343 *dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
00344 evtmc=thevt.neumc;
00345 }
00346
00347 TClonesArray& thshwTca=(*ntp.thshw);
00348 const Int_t numthshws=thshwTca.GetEntriesFast();
00349 if (numthshws>0) {
00350 Int_t bestShower=lib.reco.GetBestShower(nu);
00351 const NtpSRShower* pshw=lib.reco.GetShowerWithIndexX(ntp,evt,
00352 bestShower-1);
00353 if (pshw) {
00354 const NtpSRShower& shw=*pshw;
00355 const NtpTHShower& thshw=
00356 *dynamic_cast<NtpTHShower*>(thshwTca[shw.index]);
00357 shwmc=thshw.neumc;
00358 }
00359 }
00360
00361 const NtpTHTrack& thtrk=
00362 *dynamic_cast<NtpTHTrack*>(thtrkTca[trk.index]);
00363 trkmc=thtrk.neumc;
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 TClonesArray& hepTca=(*ntp.stdhep);
00377
00378
00379 const NtpMCStdHep& stdhep=
00380 *dynamic_cast<NtpMCStdHep*>(hepTca[thtrk.trkstdhep]);
00381
00382 Int_t trkstdhepId=stdhep.IdHEP;
00383
00384 if (abs(trkstdhepId)!=13) {
00385 MSG("NuPlots",Msg::kDebug)
00386 <<"Found THTrack, index="<<thtrk.index
00387 <<", neumc="<<thtrk.neumc
00388 <<", neustdhep="<<thtrk.neustdhep
00389 <<", trkstdhep="<<thtrk.trkstdhep
00390 <<", trkstdhepId="<<trkstdhepId
00391 <<", purity="<<thtrk.purity
00392 <<", completeall="<<thtrk.completeall
00393 <<", completeslc="<<thtrk.completeslc
00394 <<endl;
00395 }
00396
00397
00398
00399 if (lib.reco.IsTrkWithDimuon(ntp,trk)) {
00400 if (trkstdhepId>0) trkstdhepId+=100000;
00401 else trkstdhepId+=-100000;
00402 }
00403 if (lib.reco.IsTrkWithMuonFromPionDecay(ntp,trk)) {
00404 if (trkstdhepId>0) trkstdhepId+=900000;
00405 else trkstdhepId+=-900000;
00406 }
00407 if (lib.reco.IsTrkWithMuonFromKaonDecay(ntp,trk)) {
00408 if (trkstdhepId>0) trkstdhepId+=9000000;
00409 else trkstdhepId+=-9000000;
00410 }
00411 if (lib.reco.IsTrkWithMuonFromNotNuNotPiNotKaonNotCharm(ntp,
00412 trk)) {
00413 if (trkstdhepId>0) trkstdhepId+=90000000;
00414 else trkstdhepId+=-90000000;
00415 }
00416
00417 Bool_t print=false;
00418
00419
00420 if (abs(trkstdhepId)==22 && print) {
00421 MSG("NuPlots",Msg::kInfo)
00422 <<"Found THTrack, index="<<thtrk.index
00423 <<", neumc="<<thtrk.neumc
00424 <<", neustdhep="<<thtrk.neustdhep
00425 <<", trkstdhep="<<thtrk.trkstdhep
00426 <<", trkstdhepId="<<trkstdhepId
00427 <<", purity="<<thtrk.purity
00428 <<", completeall="<<thtrk.completeall
00429 <<", completeslc="<<thtrk.completeslc
00430 <<endl;
00431 TClonesArray& mcTca=*ntp.mc;
00432 const NtpMCTruth& mc=
00433 *(dynamic_cast<NtpMCTruth*>(mcTca[thtrk.neumc]));
00434 lib.reco.PrintStdhep(ntp,mc);
00435 cout<<nu.run<<" "<<nu.snarl<<" "<<nu.evt<<endl;
00436 }
00437
00438
00439 mtrkstdhepId[trkstdhepId]++;
00440 }
00441
00442
00443
00444 void NuCounter::PrintStdhepCounter() const
00445 {
00446 const NuGeneral general;
00447
00448 typedef map<Int_t,Int_t>::const_iterator mIter;
00449
00450 Int_t totalEvents=0;
00451
00452 MsgFormat pcfmt("%3.2f");
00453 MsgFormat errfmt("%.1f");
00454
00455
00456 for (mIter it=mtrkstdhepId.begin();it!=mtrkstdhepId.end();++it) {
00457 totalEvents+=it->second;
00458 }
00459 if (totalEvents==0) return;
00460
00461 multimap<Int_t,Int_t> mOpp;
00462
00463 cout<<endl<<"Particle associated with track:"<<endl;
00464
00465
00466 for (mIter it=mtrkstdhepId.begin();it!=mtrkstdhepId.end();++it) {
00467 mOpp.insert(pair<Int_t,Int_t>(it->second,it->first));
00468
00469 string special="";
00470 Int_t trkstdhepId=it->first;
00471 Int_t N=it->second;
00472
00473 if (abs(trkstdhepId)>90000000){
00474 if (trkstdhepId>0) trkstdhepId-=90000000;
00475 else trkstdhepId+=90000000;
00476 special="What is this???";
00477 }
00478 else if (abs(trkstdhepId)>9000000){
00479 if (trkstdhepId>0) trkstdhepId-=9000000;
00480 else trkstdhepId+=9000000;
00481 special="k->mu";
00482 }
00483 else if (abs(trkstdhepId)>900000){
00484 if (trkstdhepId>0) trkstdhepId-=900000;
00485 else trkstdhepId+=900000;
00486 special="pi->mu";
00487 }
00488 else if (abs(trkstdhepId)>100000){
00489 if (trkstdhepId>0) trkstdhepId-=100000;
00490 else trkstdhepId+=100000;
00491 special="dimuon";
00492 }
00493
00494 string idhep=general.IdHEPAsString(trkstdhepId);
00495
00496 Float_t err=1;
00497 if (N!=0) err=sqrt(static_cast<Float_t>(N));
00498
00499 Float_t pc=100.*N/totalEvents;
00500 cout<<idhep<<" (id="<<trkstdhepId<<"), N="<<N
00501 <<" +/- "<<errfmt(err)
00502 <<" ("<<pcfmt(pc)<<"%)"<<" "<<special<<endl;
00503 }
00504
00505
00506 cout<<endl<<"Particle associated with track:"<<endl;
00507
00508 for (mIter it=(--(mOpp.end()));it!=(--(mOpp.begin()));--it) {
00509 string special="";
00510
00511
00512 Int_t trkstdhepId=it->second;
00513 Int_t N=it->first;
00514
00515 if (abs(trkstdhepId)>90000000){
00516 if (trkstdhepId>0) trkstdhepId-=90000000;
00517 else trkstdhepId+=90000000;
00518 special="What is this???";
00519 }
00520 else if (abs(trkstdhepId)>9000000){
00521 if (trkstdhepId>0) trkstdhepId-=9000000;
00522 else trkstdhepId+=9000000;
00523 special="k->mu";
00524 }
00525 else if (abs(trkstdhepId)>900000){
00526 if (trkstdhepId>0) trkstdhepId-=900000;
00527 else trkstdhepId+=900000;
00528 special="pi->mu";
00529 }
00530 else if (abs(trkstdhepId)>100000){
00531 if (trkstdhepId>0) trkstdhepId-=100000;
00532 else trkstdhepId+=100000;
00533 special="dimuon";
00534 }
00535
00536 string idhep=general.IdHEPAsString(trkstdhepId);
00537
00538 Float_t err=1;
00539 if (N!=0) err=sqrt(static_cast<Float_t>(N));
00540
00541 Float_t pc=100.*N/totalEvents;
00542 cout<<idhep<<" (id="<<trkstdhepId<<"), N="<<N
00543 <<" +/- "<<errfmt(err)
00544 <<" ("<<pcfmt(pc)<<"%)"<<" "<<special<<endl;
00545 }
00546 }
00547
00548
00549
