11 #include <CLHEP/Vector/LorentzVector.h> 23 #include "TProfile2D.h" 25 #include "TGraphErrors.h" 30 #include "TLorentzVector.h" 47 histoDir_( outputFile->GetDirectory(name) )
65 virtual void Fill(
const CLHEP::HepLorentzVector & momentum1,
const CLHEP::HepLorentzVector & momentum2 ) {};
66 virtual void Fill(
const CLHEP::HepLorentzVector & momentum1,
const CLHEP::HepLorentzVector & momentum2,
const int charge,
const double &
weight = 1.) {};
74 virtual void Fill(
const CLHEP::HepLorentzVector & momentum,
const int charge,
const double &
weight =1.){};
79 virtual void Fill(
const CLHEP::HepLorentzVector &
p,
const double & likeValue ) {};
80 virtual void Fill(
const unsigned int number ) {};
85 const double & recoMass,
const double & genMass ) {};
90 const double & recoMass,
const double & genMass ) {};
95 virtual void Fill(
const double &
x,
const double &
y ) {};
96 virtual void Fill(
const double &
x,
const double &
y,
const double &
a,
const double &
b ) {};
100 const double &
weight = 1.) {};
101 virtual void Fill(
const CLHEP::HepLorentzVector & momentum1,
102 const CLHEP::HepLorentzVector & momentum2,
103 const CLHEP::HepLorentzVector & momentumRes,
104 const double &
weight = 1.) {};
109 virtual void Write() = 0;
110 virtual void Clear() = 0;
135 const int xBins,
const double &
xMin,
const double &
xMax,
136 const int yBins,
const double &
yMin,
const double &
yMax ) : Histograms(outputFile, dirName),
137 tH2d_( new TH2D(name, title, xBins, xMin, xMax, yBins, yMin, yMax) ),
138 tProfile_( new TProfile(name+
"Prof", title+
" profile", xBins, xMin, xMax, yMin, yMax) ) {}
143 void Fill(
const double &
x,
const double &
y )
override {
145 tProfile_->Fill(x,y);
148 if(histoDir_ !=
nullptr) histoDir_->cd();
157 tH2d_->GetXaxis()->SetTitle(title);
158 tProfile_->GetXaxis()->SetTitle(title);
161 tH2d_->GetYaxis()->SetTitle(title);
162 tProfile_->GetYaxis()->SetTitle(title);
176 const int xBins,
const double &
xMin,
const double &
xMax ) : Histograms(outputFile, name),
177 tH1D_( new TH1D(name, title, xBins, xMin, xMax) ) {}
181 void Fill(
const double &
x,
const double &
y )
override {
185 if(histoDir_ !=
nullptr) histoDir_->cd();
192 tH1D_->GetXaxis()->SetTitle(title);
195 tH1D_->GetYaxis()->SetTitle(title);
207 const int xBins,
const double &
xMin,
const double &
xMax,
208 const double &
yMin,
const double &
yMax ) : Histograms(outputFile, name),
209 tProfile_( new TProfile(name+
"Prof", title+
" profile", xBins, xMin, xMax, yMin, yMax) ) {}
213 void Fill(
const double &
x,
const double &
y )
override {
214 tProfile_->Fill(x,y);
217 if(histoDir_ !=
nullptr) histoDir_->cd();
224 tProfile_->GetXaxis()->SetTitle(title);
227 tProfile_->GetYaxis()->SetTitle(title);
242 hPt_( new TH1F (name+
"_Pt",
"transverse momentum", 100, 0,
maxPt) ),
243 hPtVsEta_( new TH2F (name+
"_PtVsEta",
"transverse momentum vs #eta", 100, 0,
maxPt, 100, -3.0, 3.0) ),
245 hCurvVsEtaNeg_( new TProfile(name+
"_CurvVsEtaNeg",
"q/pT vs #eta neg.", 64, -3.2, 3.2, -1., 0.) ),
246 hCurvVsEtaPos_( new TProfile(name+
"_CurvVsEtaPos",
"q/pT vs #eta pos.", 64, -3.2, 3.2, 0., 1.) ),
247 hCurvVsPhiNeg_( new TProfile(name+
"_CurvVsPhiNeg",
"q/pT vs #phi neg.", 32, -3.2, 3.2, -1., 0.) ),
248 hCurvVsPhiPos_( new TProfile(name+
"_CurvVsPhiPos",
"q/pT vs #phi pos.", 32, -3.2, 3.2, 0., 1.) ),
250 hPtVsPhiNeg_( new TProfile(name+
"_PtVsPhiNeg",
"pT vs #phi neg.", 32, -3.2, 3.2, 0.,100) ),
251 hPtVsPhiPos_( new TProfile(name+
"_PtVsPhiPos",
"pT vs #phi pos.", 32, -3.2, 3.2, 0.,100) ),
254 hEta_( new TH1F (name+
"_Eta",
"pseudorapidity", 64, -3.2, 3.2) ),
255 hPhi_( new TH1F (name+
"_Phi",
"phi angle", 64, -3.2, 3.2) ),
256 hMass_( new TH1F (name+
"_Mass",
"mass", 10000,
minMass,
maxMass) ),
257 hNumber_( new TH1F (name+
"_Number",
"number", 20, -0.5, 19.5) )
262 Histograms(outputFile, name)
265 hPt_ =
new TH1F (name+
"_Pt",
"transverse momentum", 100, 0,
maxPt);
266 hPtVsEta_ =
new TH2F (name+
"_PtVsEta",
"transverse momentum vs #eta", 100, 0,
maxPt, 100, -3.0, 3.0);
268 hPtVsEta_ =
new TH2F (name+
"_PtVsEta",
"transverse momentum vs #eta", 100, 0,
maxPt, 100, -3.0, 3.0);
270 hCurvVsEtaNeg_ =
new TProfile(name+
"_CurvVsEtaNeg",
"q/pT vs #eta neg.", 100, -3.0, 3.0, -1. ,0.);
271 hCurvVsEtaPos_ =
new TProfile(name+
"_CurvVsEtaPos",
"q/pT vs #eta pos.", 100, -3.0, 3.0, 0., 1.);
272 hCurvVsPhiNeg_ =
new TProfile(name+
"_CurvVsPhiNeg",
"q/pT vs #phi neg.", 32, -3.2, 3.2, -1. ,0.);
273 hCurvVsPhiPos_ =
new TProfile(name+
"_CurvVsPhiPos",
"q/pT vs #phi pos.", 32, -3.2, 3.2, 0., 1.);
275 hPtVsPhiNeg_ =
new TProfile(name+
"_PtVsPhiNeg",
"pT vs #phi neg.", 32, -3.2, 3.2, 0.,100);
276 hPtVsPhiPos_ =
new TProfile(name+
"_PtVsPhiPos",
"pT vs #phi pos.", 32, -3.2, 3.2, 0.,100);
281 hEta_ =
new TH1F (name+
"_Eta",
"pseudorapidity", 64, -3.2, 3.2);
282 hPhi_ =
new TH1F (name+
"_Phi",
"phi angle", 64, -3.2, 3.2);
283 hMass_ =
new TH1F (name+
"_Mass",
"mass", 40000,
minMass,
maxMass);
284 hNumber_ =
new TH1F (name+
"_Number",
"number", 20, -0.5, 19.5);
289 hPt_( (TH1F *) file->Get(name_+
"_Pt") ),
290 hPtVsEta_( (TH2F *) file->Get(name_+
"_PtVsEta") ),
293 hCurvVsEtaNeg_( (TProfile *) file->Get(name_+
"_CurvVsEtaNeg") ),
294 hCurvVsEtaPos_( (TProfile *) file->Get(name_+
"_CurvVsEtaPos") ),
295 hCurvVsPhiNeg_( (TProfile *) file->Get(name_+
"_CurvVsPhiNeg") ),
296 hCurvVsPhiPos_( (TProfile *) file->Get(name_+
"_CurvVsPhiPos") ),
298 hPtVsPhiNeg_( (TProfile *) file->Get(name_+
"_PtVsPhiNeg") ),
299 hPtVsPhiPos_( (TProfile *) file->Get(name_+
"_PtVsPhiPos") ),
301 hEta_( (TH1F *) file->Get(name_+
"_Eta") ),
302 hPhi_( (TH1F *) file->Get(name_+
"_Phi") ),
303 hMass_( (TH1F *) file->Get(name_+
"_Mass") ),
305 hNumber_( (TH1F *) file->Get(name_+
"_Number") )
313 delete hCurvVsEtaNeg_;
314 delete hCurvVsEtaPos_;
315 delete hCurvVsPhiNeg_;
316 delete hCurvVsPhiPos_;
330 Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()),charge,
weight);
333 void Fill(
const CLHEP::HepLorentzVector & momentum,
const int charge,
const double &
weight =1.)
override 335 hPt_->Fill(momentum.perp(),
weight);
336 hPtVsEta_->Fill(momentum.perp(), momentum.eta(),
weight);
339 if(charge<0)hCurvVsEtaNeg_->Fill( momentum.eta(),charge/(momentum.perp()),
weight);
340 if(charge>0)hCurvVsEtaPos_->Fill( momentum.eta(),charge/(momentum.perp()),
weight);
341 if(charge<0)hCurvVsPhiNeg_->Fill( momentum.phi(),charge/(momentum.perp()),
weight);
342 if(charge>0)hCurvVsPhiPos_->Fill( momentum.phi(),charge/(momentum.perp()),
weight);
344 if(charge<0)hPtVsPhiNeg_->Fill( momentum.phi(),momentum.perp(),
weight);
345 if(charge>0)hPtVsPhiPos_->Fill( momentum.phi(),momentum.perp(),
weight);
347 hEta_->Fill(momentum.eta(),
weight);
348 hPhi_->Fill(momentum.phi(),
weight);
349 hMass_->Fill(momentum.m(),
weight);
355 void Fill(
unsigned int number )
override 357 hNumber_->Fill(number);
362 if(histoDir_ !=
nullptr) histoDir_->cd();
366 hCurvVsEtaNeg_->Write();
367 hCurvVsEtaPos_->Write();
368 hCurvVsPhiNeg_->Write();
369 hCurvVsPhiPos_->Write();
371 hPtVsPhiNeg_->Write();
372 hPtVsPhiPos_->Write();
386 hCurvVsEtaNeg_->Clear();
387 hCurvVsEtaPos_->Clear();
388 hCurvVsPhiNeg_->Clear();
389 hCurvVsPhiPos_->Clear();
391 hPtVsPhiNeg_->Clear();
392 hPtVsPhiPos_->Clear();
430 hEta_( new TH1F (name+
"_DeltaEta",
"#Delta#eta", 100, 0, 6) ),
431 hEtaSign_( new TH1F (name+
"_DeltaEtaSign",
"#Delta#eta with sign", 100, -6, 6) ),
432 hPhi_( new TH1F (name+
"_DeltaPhi",
"#Delta#phi", 100,0,3.2) ),
433 hTheta_( new TH1F (name+
"_DeltaTheta",
"#Delta#theta", 100,-3.2,3.2) ),
434 hCotgTheta_( new TH1F (name+
"_DeltaCotgTheta",
"#Delta Cotg(#theta )", 100,-3.2,3.2) ),
435 hDeltaR_( new TH1F (name+
"_DeltaR",
"#Delta R", 400, 0, 4 ) )
439 Histograms(outputFile, name),
442 hEta_( new TH1F (name+
"_DeltaEta",
"#Delta#eta", 100, 0, 6) ),
443 hEtaSign_( new TH1F (name+
"_DeltaEtaSign",
"#Delta#eta with sign", 100, -6, 6) ),
444 hPhi_( new TH1F (name+
"_DeltaPhi",
"#Delta#phi", 100,0,3.2) ),
445 hTheta_( new TH1F (name+
"_DeltaTheta",
"#Delta#theta", 100,-3.2,3.2) ),
446 hCotgTheta_( new TH1F (name+
"_DeltaCotgTheta",
"#Delta Cotg(#theta )", 100,-3.2,3.2) ),
447 hDeltaR_( new TH1F (name+
"_DeltaR",
"#DeltaR", 400, 0, 4 ) )
452 hEta_ = (TH1F *) file->Get(name+
"_DeltaEta");
453 hEtaSign_ = (TH1F *) file->Get(name+
"_DeltaEtaSign");
454 hPhi_ = (TH1F *) file->Get(name+
"_DeltaPhi");
455 hTheta_ = (TH1F *) file->Get(name+
"_DeltaTheta");
456 hCotgTheta_ = (TH1F *) file->Get(name+
"_DeltaCotgTheta");
457 hDeltaR_ = (TH1F *) file->Get(name+
"_DeltaR");
470 Fill (CLHEP::HepLorentzVector(p1.x(),p1.y(),p1.z(),p1.t()),
471 CLHEP::HepLorentzVector(p2.x(),p2.y(),p2.z(),p2.t()));
475 Fill (p1,CLHEP::HepLorentzVector(p2.x(),p2.y(),p2.z(),p2.t()));
478 void Fill (
const CLHEP::HepLorentzVector & momentum1,
const CLHEP::HepLorentzVector & momentum2)
override {
479 hEta_->Fill(fabs( momentum1.eta()-momentum2.eta() ));
480 hEtaSign_->Fill(momentum1.eta()-momentum2.eta());
482 hTheta_->Fill(momentum1.theta()-momentum2.theta());
484 double theta1 = momentum1.theta();
485 double theta2 = momentum2.theta();
486 hCotgTheta_->Fill(TMath::Cos(theta1)/TMath::Sin(theta1) - TMath::Cos(theta2)/TMath::Sin(theta2));
487 hDeltaR_->Fill(
sqrt((momentum1.eta()-momentum2.eta())*(momentum1.eta()-momentum2.eta()) +
493 if(histoDir_ !=
nullptr) histoDir_->cd();
499 hCotgTheta_->Write();
509 hCotgTheta_->Clear();
530 hPtVSEta_ =
new TH2F( name+
"_PtVSEta",
"transverse momentum vs pseudorapidity",
531 32, -3.2, 3.2, 200, 0,
maxPt );
532 hMassVSEta_ =
new TH2F( name+
"_MassVSEta",
"mass vs pseudorapidity",
536 hPtVSEta_prof_ =
new TProfile( name+
"_PtVSEta_prof",
"mass vs pseudorapidity",
537 32, -3.2, 3.2, 0,
maxPt );
538 hMassVSEta_prof_ =
new TProfile( name+
"_MassVSEta_prof",
"mass vs pseudorapidity",
540 hCurvVSEta_prof_ =
new TProfile( name+
"_CurvVSEta_prof",
"curvature vs pseudorapidity",
541 32, -3.2, 3.2, 0, 1. );
547 delete hPtVSEta_prof_;
548 delete hMassVSEta_prof_;
549 delete hCurvVSEta_prof_;
553 Fill (CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()),
weight);
556 void Fill (
const CLHEP::HepLorentzVector & momentum,
const double &
weight = 1.)
override {
557 hPtVSEta_->Fill(momentum.eta(),momentum.perp(),
weight);
558 hPtVSEta_prof_->Fill(momentum.eta(),momentum.perp(),
weight);
560 hMassVSEta_->Fill(momentum.eta(),momentum.m(),
weight);
561 hMassVSEta_prof_->Fill(momentum.eta(),momentum.m(),
weight);
562 hCurvVSEta_prof_->Fill(momentum.eta(),1/momentum.perp(),
weight);
567 hPtVSEta_prof_->Write();
568 hCurvVSEta_prof_->Write();
569 hMassVSEta_->Write();
570 hMassVSEta_prof_->Write();
580 hPtVSEta_prof_->Clear();
581 hCurvVSEta_prof_->Clear();
582 hMassVSEta_->Clear();
583 hMassVSEta_prof_->Clear();
624 hMassVSPhi_prof_ =
new TProfile (name+
"_MassVSPhi_prof",
"mass vs phi angle",
625 16, -3.2, 3.2, 70, 110);
626 hPtVSPhi_prof_ =
new TProfile (name+
"_PtVSPhi_prof",
"pt vs phi angle",
627 16, -3.2, 3.2, 0, 200);
633 delete hMassVSPhi_prof_;
634 delete hPtVSPhi_prof_;
646 Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()),
weight);
649 void Fill(
const CLHEP::HepLorentzVector & momentum,
const double &
weight= 1.)
override {
650 hPtVSPhi_->Fill(momentum.phi(),momentum.perp(),
weight);
651 hMassVSPhi_->Fill(momentum.phi(),momentum.m(),
weight);
652 hMassVSPhi_prof_->Fill(momentum.phi(),momentum.m(),
weight);
653 hPtVSPhi_prof_->Fill(momentum.phi(),momentum.perp(),
weight);
666 hMassVSPhi_->Write();
667 hMassVSPhi_prof_->Write();
668 hPtVSPhi_prof_->Write();
714 hMassVSPhi_->Clear();
715 hPtVSPhi_prof_->Clear();
716 hMassVSPhi_prof_->Clear();
749 hMassVSPt_ =
new TH2F( name+
"_MassVSPt",
"mass vs transverse momentum",
750 12, -6, 6, 40, 70, 110 );
752 hMassVSPt_prof_ =
new TProfile( name+
"_MassVSPt_prof",
"mass vs transverse momentum",
753 12, -3, 3, 86, 116 );
758 delete hMassVSPt_prof_;
762 Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()),
weight);
765 void Fill(
const CLHEP::HepLorentzVector & momentum,
const double &
weight = 1.)
override {
766 hMassVSPt_->Fill(momentum.eta(),momentum.m(),
weight);
767 hMassVSPt_prof_->Fill(momentum.eta(),momentum.m(),
weight);
772 hMassVSPt_prof_->Write();
782 hMassVSPt_prof_->Clear();
800 hMassVSPt_ =
new TH2F( name+
"_MassVSPt",
"resonance mass vs muon transverse momentum", 200, 0.,
maxPt, 6000,
minMass,
maxMass );
802 hMassVSEta_ =
new TH2F( name+
"_MassVSEta",
"resonance mass vs muon pseudorapidity", 64, -6.4, 6.4, 6000,
minMass,
maxMass );
803 hMassVSEtaPlus_ =
new TH2F( name+
"_MassVSEtaPlus",
"resonance mass vs muon+ pseudorapidity", 64, -6.4, 6.4, 6000,
minMass,
maxMass );
804 hMassVSEtaMinus_ =
new TH2F( name+
"_MassVSEtaMinus",
"resonance mass vs muon- pseudorapidity", 64, -6.4, 6.4, 6000,
minMass,
maxMass );
806 hMassVSPhiPlus_ =
new TH2F( name+
"_MassVSPhiPlus",
"resonance mass vs muon+ phi angle", 64, -3.2, 3.2, 6000,
minMass,
maxMass );
807 hMassVSPhiMinus_ =
new TH2F( name+
"_MassVSPhiMinus",
"resonance mass vs muon- phi angle", 64, -3.2, 3.2, 6000,
minMass,
maxMass );
815 hMassVSEtaPhiPlus_ =
new TH3F( name+
"_MassVSEtaPhiPlus",
"resonance mass vs muon+ phi/pseudorapidity", 16, -3.2, 3.2, 20, -2.4, 2.4, 300,
minMass,
maxMass );
816 hMassVSEtaPhiMinus_ =
new TH3F( name+
"_MassVSEtaPhiMinus",
"resonance mass vs muon- phi/pseudorapidity", 16, -3.2, 3.2, 20, -2.4, 2.4, 300,
minMass,
maxMass );
819 hMassVSCosThetaCS_ =
new TH2F( name+
"_MassVSCosThetaCS",
"resonance mass vs cos(theta) (CS frame)", 40, -1., 1., 6000,
minMass,
maxMass );
820 hMassVSPhiCS_ =
new TH2F( name+
"_MassVSPhiCS",
"resonance mass vs phi (CS frame)", 64, -3.2, 3.2, 6000,
minMass,
maxMass );
823 hMassVSPhiPlusPhiMinus_ =
new TH3F( name+
"_MassVSPhiPlusPhiMinus",
"resonance mass vs muon+ phi/muon- phi",16, -3.2, 3.2,16, -3.2, 3.2, 6000,
minMass,
maxMass );
824 hMassVSEtaPlusEtaMinus_ =
new TH3F( name+
"_MassVSEtaPlusEtaMinus",
"resonance mass vs muon+ eta/muon- eta",16, -3.2, 3.2,16, -3.2, 3.2, 6000,
minMass,
maxMass );
827 hMassVSPhiPlusMinusDiff_ =
new TH2F( name+
"_MassVSPhiPlusMinusDiff",
"resonance mass vs delta phi between mu+/mu-", 64, -6.4, 6.4, 6000,
minMass,
maxMass );
828 hMassVSEtaPlusMinusDiff_ =
new TH2F( name+
"_MassVSEtaPlusMinusDiff",
"resonance mass vs delta pseudorapidity between mu+/mu-", 32, -4.4, 4.4, 6000,
minMass,
maxMass );
829 hMassVSCosThetaCS_prof =
new TProfile (name+
"_MassVScosTheta_prof",
"resonance mass vs cosTheta", 40, -1., 1., 85., 95.);
839 hMassVSPt_ = (TH2F *) file->Get(name+
"_MassVSPt");
840 hMassVSEta_ = (TH2F *) file->Get(name+
"_MassVSEta");
842 hMassVSEtaPhiPlus_ = (TH3F *) file->Get(name+
"_MassVSEtaPlus");
843 hMassVSEtaPhiMinus_ = (TH3F *) file->Get(name+
"_MassVSEtaMinus");
844 hMassVSEtaPlus_ = (TH2F *) file->Get(name+
"_MassVSEtaPlus");
845 hMassVSEtaMinus_ = (TH2F *) file->Get(name+
"_MassVSEtaMinus");
847 hMassVSPhiPlusMinusDiff_ = (TH2F *) file->Get(name+
"_MassVSPhiPlusMinusDiff");
848 hMassVSEtaPlusMinusDiff_ = (TH2F *) file->Get(name+
"_MassVSEtaPlusMinusDiff");
850 hMassVSPhiPlus_ = (TH2F *) file->Get(name+
"_MassVSPhiPlus");
851 hMassVSPhiMinus_ = (TH2F *) file->Get(name+
"_MassVSPhiMinus");
853 hMassVSCosThetaCS_prof = (TProfile *) file->Get(name+
"_MassVScosTheta_prof");
863 delete hMassVSPhiPlus_;
864 delete hMassVSPhiMinus_;
865 delete hMassVSEtaPhiPlus_;
866 delete hMassVSEtaPhiMinus_;
867 delete hMassVSEtaPlus_;
868 delete hMassVSEtaMinus_;
869 delete hMassVSPhiPlusPhiMinus_;
870 delete hMassVSEtaPlusEtaMinus_;
871 delete hMassVSCosThetaCS_;
872 delete hMassVSPhiCS_;
873 delete hMassVSPhiPlusMinusDiff_;
874 delete hMassVSEtaPlusMinusDiff_;
875 delete hMassVSCosThetaCS_prof;
880 Fill(CLHEP::HepLorentzVector(p41.x(),p41.y(),p41.z(),p41.t()),
881 CLHEP::HepLorentzVector(p42.x(),p42.y(),p42.z(),p42.t()), charge,
weight);
887 const double &
weight = 1.)
override 889 Fill(CLHEP::HepLorentzVector(p41.x(),p41.y(),p41.z(),p41.t()),
890 CLHEP::HepLorentzVector(p42.x(),p42.y(),p42.z(),p42.t()),
891 CLHEP::HepLorentzVector(p4Res.x(),p4Res.y(),p4Res.z(),p4Res.t()),
896 void Fill(
const CLHEP::HepLorentzVector & momentum1,
897 const CLHEP::HepLorentzVector & momentum2,
898 const CLHEP::HepLorentzVector & momentumRes,
899 const double &
weight = 1.)
override 909 double costhetaCS, phiCS;
911 const CLHEP::HepLorentzVector&
mu= momentum1;
912 const CLHEP::HepLorentzVector&
mubar= momentum2;
913 CLHEP::HepLorentzVector
Q(mu+mubar);
914 double muplus = 1.0/
sqrt(2.0) * (mu.e() + mu.z());
915 double muminus = 1.0/
sqrt(2.0) * (mu.e() - mu.z());
916 double mubarplus = 1.0/
sqrt(2.0) * (mubar.e() + mubar.z());
917 double mubarminus = 1.0/
sqrt(2.0) * (mubar.e() - mubar.z());
919 costhetaCS = 2.0 / Q.mag() /
sqrt(
pow(Q.mag(), 2) +
pow(Q.perp(), 2)) * (muplus * mubarminus - muminus * mubarplus);
920 if (momentumRes.rapidity()<0) costhetaCS = -costhetaCS;
932 CLHEP::HepLorentzVector Pbeam(0.,0.,3500.,3500.);
933 CLHEP::Hep3Vector
R = Pbeam.vect().cross(Q.vect());
934 CLHEP::Hep3Vector Runit = R.unit();
938 CLHEP::Hep3Vector Qt = Q.vect(); Qt.setZ(0);
939 CLHEP::Hep3Vector Qtunit = Qt.unit();
942 CLHEP::HepLorentzVector
D(mu-mubar);
943 CLHEP::Hep3Vector Dt = D.vect(); Dt.setZ(0);
944 double tanphi =
sqrt(
pow(Q.mag(), 2) +
pow(Q.perp(), 2)) / Q.mag() * Dt.dot(Runit) / Dt.dot(Qtunit);
945 if (momentumRes.rapidity()<0) tanphi = -tanphi;
946 phiCS = atan(tanphi);
948 hMassVSPhiCS_->Fill(phiCS,momentumRes.m(),
weight);
949 hMassVSCosThetaCS_->Fill(costhetaCS,momentumRes.m(),
weight);
950 hMassVSCosThetaCS_prof ->Fill(costhetaCS,momentumRes.m());
954 hMassVSPhiPlusPhiMinus_->Fill(momentum1.phi(), momentum2.phi(), momentumRes.m(),
weight);
955 hMassVSEtaPlusEtaMinus_->Fill(momentum1.eta(), momentum2.eta(), momentumRes.m(),
weight);
957 hMassVSPhiPlusMinusDiff_->Fill( (momentum1.phi()-momentum2.phi()), momentumRes.m(),
weight);
958 hMassVSEtaPlusMinusDiff_->Fill( (momentum1.eta()-momentum2.eta()), momentumRes.m(),
weight);
961 void Fill(
const CLHEP::HepLorentzVector & momentum1,
const CLHEP::HepLorentzVector & momentum2,
const int charge,
const double &
weight = 1.)
override 963 hMassVSPt_->Fill(momentum1.perp(),momentum2.m(),
weight);
967 hMassVSEta_->Fill(momentum1.eta(),momentum2.m(),
weight);
971 hMassVSPhiPlus_->Fill(momentum1.phi(),momentum2.m(),
weight);
972 hMassVSEtaPlus_->Fill(momentum1.eta(),momentum2.m(),
weight);
973 hMassVSEtaPhiPlus_->Fill(momentum1.phi(),momentum1.eta(),momentum2.m(),
weight);
976 hMassVSPhiMinus_->Fill(momentum1.phi(),momentum2.m(),
weight);
977 hMassVSEtaMinus_->Fill(momentum1.eta(),momentum2.m(),
weight);
978 hMassVSEtaPhiMinus_->Fill(momentum1.phi(),momentum1.eta(),momentum2.m(),
weight);
981 LogDebug(
"Histograms") <<
"HMassVSPart: wrong charge value = " << charge << std::endl;
988 hMassVSEta_->Write();
989 hMassVSPhiPlus_->Write();
990 hMassVSPhiMinus_->Write();
992 hMassVSEtaPhiPlus_->Write();
993 hMassVSEtaPhiMinus_->Write();
994 hMassVSEtaPlus_->Write();
995 hMassVSEtaMinus_->Write();
997 hMassVSPhiPlusPhiMinus_->Write();
998 hMassVSEtaPlusEtaMinus_->Write();
999 hMassVSCosThetaCS_->Write();
1000 hMassVSPhiCS_->Write();
1002 hMassVSPhiPlusMinusDiff_->Write();
1003 hMassVSEtaPlusMinusDiff_->Write();
1004 hMassVSCosThetaCS_prof->Write();
1013 hMassVSPt_->Clear();
1014 hMassVSEta_->Clear();
1015 hMassVSPhiPlus_->Clear();
1016 hMassVSPhiMinus_->Clear();
1018 hMassVSEtaPhiPlus_->Clear();
1019 hMassVSEtaPhiMinus_->Clear();
1020 hMassVSEtaPlus_->Clear();
1021 hMassVSEtaMinus_->Clear();
1023 hMassVSPhiPlusPhiMinus_->Clear();
1024 hMassVSEtaPlusEtaMinus_->Clear();
1025 hMassVSCosThetaCS_->Clear();
1026 hMassVSPhiCS_->Clear();
1027 hMassVSPhiPlusMinusDiff_->Clear();
1028 hMassVSEtaPlusMinusDiff_->Clear();
1029 hMassVSCosThetaCS_prof->Clear();
1075 hMassVSPt_ =
new TProfile2D( name+
"_MassVSPt",
"resonance mass vs muon transverse momentum", 200, 0.,
maxPt, 6000,
minMass,
maxMass, 0., 100. );
1076 hMassVSEta_ =
new TProfile2D( name+
"_MassVSEta",
"resonance mass vs muon pseudorapidity", 64, -6.4, 6.4, 6000,
minMass,
maxMass, 0., 100. );
1077 hMassVSPhiPlus_ =
new TProfile2D( name+
"_MassVSPhiPlus",
"resonance mass vs muon+ phi angle", 64, -3.2, 3.2, 6000,
minMass,
maxMass, 0., 100. );
1078 hMassVSPhiMinus_ =
new TProfile2D( name+
"_MassVSPhiMinus",
"resonance mass vs muon- phi angle", 64, -3.2, 3.2, 6000,
minMass,
maxMass, 0., 100. );
1083 hMassVSPt_ = (TProfile2D *) file->Get(name+
"_MassVSPt");
1084 hMassVSEta_ = (TProfile2D *) file->Get(name+
"_MassVSEta");
1085 hMassVSPhiPlus_ = (TProfile2D *) file->Get(name+
"_MassVSPhiPlus");
1086 hMassVSPhiMinus_ = (TProfile2D *) file->Get(name+
"_MassVSPhiMinus");
1092 delete hMassVSPhiPlus_;
1093 delete hMassVSPhiMinus_;
1097 Fill(CLHEP::HepLorentzVector(p41.x(),p41.y(),p41.z(),p41.t()),
1098 CLHEP::HepLorentzVector(p42.x(),p42.y(),p42.z(),p42.t()), charge,
weight);
1101 void Fill(
const CLHEP::HepLorentzVector & momentum1,
const CLHEP::HepLorentzVector & momentum2,
const int charge,
const double &
weight = 1.)
override {
1102 hMassVSPt_->Fill(momentum1.perp(),momentum2.m(),
weight);
1103 hMassVSEta_->Fill(momentum1.eta(),momentum2.m(),
weight);
1105 hMassVSPhiPlus_->Fill(momentum1.phi(), momentum2.m(),
weight);
1108 hMassVSPhiMinus_->Fill(momentum1.phi(), momentum2.m(),
weight);
1111 LogDebug(
"Histograms") <<
"HMassVSPartProfile: wrong charge value = " << charge << std::endl;
1117 hMassVSPt_->Write();
1118 hMassVSEta_->Write();
1119 hMassVSPhiPlus_->Write();
1120 hMassVSPhiMinus_->Write();
1124 hMassVSPt_->Clear();
1125 hMassVSEta_->Clear();
1126 hMassVSPhiPlus_->Clear();
1127 hMassVSPhiMinus_->Clear();
1143 const double & yMinEta = 0.,
const double & yMaxEta = 2.,
1144 const double & yMinPt = 0.,
const double & yMaxPt = 2.,
1145 const bool doProfiles =
false) : Histograms(outputFile, name), doProfiles_(doProfiles)
1164 hReso_ =
new TH1F( name+
"_Reso",
"resolution", 200, -1, 1 );
1165 hResoVSPtEta_ =
new TH2F( name+
"_ResoVSPtEta",
"resolution VS pt and #eta", 200, 0,
maxPt, 60, -3, 3 );
1166 hResoVSPt_ =
new TH2F( name+
"_ResoVSPt",
"resolution VS pt", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1167 hResoVSPt_Bar_ =
new TH2F( name+
"_ResoVSPt_Bar",
"resolution VS pt Barrel", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1168 hResoVSPt_Endc_17_ =
new TH2F( name+
"_ResoVSPt_Endc_1.7",
"resolution VS pt Endcap (1.4<eta<1.7)", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1169 hResoVSPt_Endc_20_ =
new TH2F( name+
"_ResoVSPt_Endc_2.0",
"resolution VS pt Endcap (1.7<eta<2.0)", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1170 hResoVSPt_Endc_24_ =
new TH2F( name+
"_ResoVSPt_Endc_2.4",
"resolution VS pt Endcap (2.0<eta<2.4)", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1171 hResoVSPt_Ovlap_ =
new TH2F( name+
"_ResoVSPt_Ovlap",
"resolution VS pt overlap", 200, 0,
maxPt, 200, yMinPt, yMaxPt );
1172 hResoVSEta_ =
new TH2F( name+
"_ResoVSEta",
"resolution VS eta", 200, -3, 3, 200, yMinEta, yMaxEta );
1173 hResoVSTheta_ =
new TH2F( name+
"_ResoVSTheta",
"resolution VS theta", 30, 0,
TMath::Pi(), 200, yMinEta, yMaxEta );
1174 hResoVSPhiPlus_ =
new TH2F( name+
"_ResoVSPhiPlus",
"resolution VS phi mu+", 16, -3.2, 3.2, 200, -1, 1 );
1175 hResoVSPhiMinus_ =
new TH2F( name+
"_ResoVSPhiMinus",
"resolution VS phi mu-", 16, -3.2, 3.2, 200, -1, 1 );
1176 hAbsReso_ =
new TH1F( name+
"_AbsReso",
"resolution", 100, 0, 1 );
1177 hAbsResoVSPt_ =
new TH2F( name+
"_AbsResoVSPt",
"Abs resolution VS pt", 200, 0,
maxPt, 100, 0, 1 );
1178 hAbsResoVSEta_ =
new TH2F( name+
"_AbsResoVSEta",
"Abs resolution VS eta", 64, -3.2, 3.2, 100, 0, 1 );
1179 hAbsResoVSPhi_ =
new TH2F( name+
"_AbsResoVSPhi",
"Abs resolution VS phi", 16, -3.2, 3.2, 100, 0, 1 );
1181 hResoVSPt_prof_ =
new TProfile(name+
"_ResoVSPt_prof",
"resolution VS pt", 100, 0,
maxPt, yMinPt, yMaxPt );
1182 hResoVSPt_Bar_prof_ =
new TProfile(name+
"_ResoVSPt_Bar_prof",
"resolution VS pt Barrel", 100, 0,
maxPt, yMinPt, yMaxPt );
1183 hResoVSPt_Endc_17_prof_ =
new TProfile(name+
"_ResoVSPt_Endc_1.7_prof",
"resolution VS pt Endcap (1.4<eta<1.7)", 100, 0,
maxPt, yMinPt, yMaxPt );
1184 hResoVSPt_Endc_20_prof_ =
new TProfile(name+
"_ResoVSPt_Endc_2.0_prof",
"resolution VS pt Endcap (1.7<eta<2.0)", 100, 0,
maxPt, yMinPt, yMaxPt );
1185 hResoVSPt_Endc_24_prof_ =
new TProfile(name+
"_ResoVSPt_Endc_2.4_prof",
"resolution VS pt Endcap (2.0<eta<2.4)", 100, 0,
maxPt, yMinPt, yMaxPt );
1186 hResoVSPt_Ovlap_prof_ =
new TProfile(name+
"_ResoVSPt_Ovlap_prof",
"resolution VS pt Overlap", 100, 0,
maxPt, yMinPt, yMaxPt );
1187 hResoVSEta_prof_ =
new TProfile(name+
"_ResoVSEta_prof",
"resolution VS eta", 200, -3.0, 3.0, yMinEta, yMaxEta );
1188 hResoVSPhi_prof_ =
new TProfile(name+
"_ResoVSPhi_prof",
"resolution VS phi", 16, -3.2, 3.2, -1, 1 );
1194 hReso_ = (TH1F *) file->Get(name+
"_Reso");
1195 hResoVSPtEta_ = (TH2F *) file->Get(name+
"_ResoVSPtEta");
1196 hResoVSPt_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1197 hResoVSPt_Bar_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1198 hResoVSPt_Endc_17_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1199 hResoVSPt_Endc_20_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1200 hResoVSPt_Endc_24_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1201 hResoVSPt_Ovlap_ = (TH2F *) file->Get(name+
"_ResoVSPt");
1202 hResoVSEta_ = (TH2F *) file->Get(name+
"_ResoVSEta");
1203 hResoVSTheta_ = (TH2F *) file->Get(name+
"_ResoVSTheta");
1204 hResoVSPhiPlus_ = (TH2F *) file->Get(name+
"_ResoVSPhiPlus");
1205 hResoVSPhiMinus_ = (TH2F *) file->Get(name+
"_ResoVSPhiMinus");
1206 hAbsReso_ = (TH1F *) file->Get(name+
"_AbsReso");
1207 hAbsResoVSPt_ = (TH2F *) file->Get(name+
"_AbsResoVSPt");
1208 hAbsResoVSEta_ = (TH2F *) file->Get(name+
"_AbsResoVSEta");
1209 hAbsResoVSPhi_ = (TH2F *) file->Get(name+
"_AbsResoVSPhi");
1211 hResoVSPt_prof_ = (TProfile *) file->Get(name+
"_ResoVSPt_prof");
1212 hResoVSPt_Bar_prof_ = (TProfile *) file->Get(name+
"_ResoVSPt_prof");
1213 hResoVSPt_Endc_17_prof_ = (TProfile *) file->Get(name+
"_ResoVSPt_1.7_prof");
1214 hResoVSPt_Endc_20_prof_ = (TProfile *) file->Get(name+
"_ResoVSPt_2.0_prof");
1215 hResoVSPt_Endc_24_prof_ = (TProfile *) file->Get(name+
"_ResoVSPt_2.4_prof");
1216 hResoVSPt_Ovlap_prof_= (TProfile *) file->Get(name+
"_ResoVSPt_prof");
1217 hResoVSEta_prof_ = (TProfile *) file->Get(name+
"_ResoVSEta_prof");
1218 hResoVSPhi_prof_ = (TProfile *) file->Get(name+
"_ResoVSPhi_prof");
1224 delete hResoVSPtEta_;
1226 delete hResoVSPt_Bar_;
1227 delete hResoVSPt_Endc_17_;
1228 delete hResoVSPt_Endc_20_;
1229 delete hResoVSPt_Endc_24_;
1230 delete hResoVSPt_Ovlap_;
1232 delete hResoVSTheta_;
1233 delete hResoVSPhiMinus_;
1234 delete hResoVSPhiPlus_;
1236 delete hAbsResoVSPt_;
1237 delete hAbsResoVSEta_;
1238 delete hAbsResoVSPhi_;
1240 delete hResoVSPt_prof_;
1241 delete hResoVSPt_Bar_prof_;
1242 delete hResoVSPt_Endc_17_prof_;
1243 delete hResoVSPt_Endc_20_prof_;
1244 delete hResoVSPt_Endc_24_prof_;
1245 delete hResoVSPt_Ovlap_prof_;
1246 delete hResoVSEta_prof_;
1247 delete hResoVSPhi_prof_;
1252 double pt = p4.Pt();
1253 double eta = p4.Eta();
1254 hReso_->Fill(resValue);
1255 hResoVSPtEta_->Fill(pt, eta,resValue);
1256 hResoVSPt_->Fill(pt,resValue);
1258 hResoVSPt_Bar_->Fill(pt,resValue);
1259 else if(fabs(eta)>0.9 && fabs(eta)<=1.4)
1260 hResoVSPt_Ovlap_->Fill(pt,resValue);
1261 else if(fabs(eta)>1.4 && fabs(eta)<=1.7)
1262 hResoVSPt_Endc_17_->Fill(pt,resValue);
1263 else if(fabs(eta)>1.7 && fabs(eta)<=2.0)
1264 hResoVSPt_Endc_20_->Fill(pt,resValue);
1266 hResoVSPt_Endc_24_->Fill(pt,resValue);
1268 hResoVSEta_->Fill(eta,resValue);
1269 hResoVSTheta_->Fill(p4.Theta(),resValue);
1271 hResoVSPhiPlus_->Fill(p4.Phi(),resValue);
1273 hResoVSPhiMinus_->Fill(p4.Phi(),resValue);
1274 hAbsReso_->Fill(fabs(resValue));
1275 hAbsResoVSPt_->Fill(pt,fabs(resValue));
1276 hAbsResoVSEta_->Fill(eta,fabs(resValue));
1277 hAbsResoVSPhi_->Fill(p4.Phi(),fabs(resValue));
1279 hResoVSPt_prof_->Fill(p4.Pt(),resValue);
1281 hResoVSPt_Bar_prof_->Fill(p4.Pt(),resValue);
1282 else if(fabs(eta)>0.9 && fabs(eta)<=1.4)
1283 hResoVSPt_Ovlap_prof_->Fill(pt,resValue);
1284 else if(fabs(eta)>1.4 && fabs(eta)<=1.7)
1285 hResoVSPt_Endc_17_prof_->Fill(pt,resValue);
1286 else if(fabs(eta)>1.7 && fabs(eta)<=2.0)
1287 hResoVSPt_Endc_20_prof_->Fill(pt,resValue);
1289 hResoVSPt_Endc_24_prof_->Fill(pt,resValue);
1291 hResoVSEta_prof_->Fill(p4.Eta(),resValue);
1292 hResoVSPhi_prof_->Fill(p4.Phi(),resValue);
1297 if(histoDir_ !=
nullptr) histoDir_->cd();
1300 hResoVSPtEta_->Write();
1301 hResoVSPt_->Write();
1302 hResoVSPt_Bar_->Write();
1303 hResoVSPt_Endc_17_->Write();
1304 hResoVSPt_Endc_20_->Write();
1305 hResoVSPt_Endc_24_->Write();
1306 hResoVSPt_Ovlap_->Write();
1307 hResoVSEta_->Write();
1308 hResoVSTheta_->Write();
1309 hResoVSPhiMinus_->Write();
1310 hResoVSPhiPlus_->Write();
1312 hAbsResoVSPt_->Write();
1313 hAbsResoVSEta_->Write();
1314 hAbsResoVSPhi_->Write();
1316 hResoVSPt_prof_->Write();
1317 hResoVSPt_Bar_prof_->Write();
1318 hResoVSPt_Endc_17_prof_->Write();
1319 hResoVSPt_Endc_20_prof_->Write();
1320 hResoVSPt_Endc_24_prof_->Write();
1321 hResoVSPt_Ovlap_prof_->Write();
1322 hResoVSEta_prof_->Write();
1323 hResoVSPhi_prof_->Write();
1329 hResoVSPtEta_->Clear();
1330 hResoVSPt_->Clear();
1331 hResoVSPt_Bar_->Clear();
1332 hResoVSPt_Endc_17_->Clear();
1333 hResoVSPt_Endc_20_->Clear();
1334 hResoVSPt_Endc_24_->Clear();
1335 hResoVSPt_Ovlap_->Clear();
1336 hResoVSEta_->Clear();
1337 hResoVSTheta_->Clear();
1338 hResoVSPhiPlus_->Clear();
1339 hResoVSPhiMinus_->Clear();
1341 hAbsResoVSPt_->Clear();
1342 hAbsResoVSEta_->Clear();
1343 hAbsResoVSPhi_->Clear();
1345 hResoVSPt_prof_->Clear();
1346 hResoVSPt_Bar_prof_->Clear();
1347 hResoVSPt_Endc_17_prof_->Clear();
1348 hResoVSPt_Endc_20_prof_->Clear();
1349 hResoVSPt_Endc_24_prof_->Clear();
1350 hResoVSPt_Ovlap_prof_->Clear();
1351 hResoVSEta_prof_->Clear();
1352 hResoVSPhi_prof_->Clear();
1395 hLikeVSPt_ =
new TH2F (name+
"_LikelihoodVSPt",
"likelihood vs muon transverse momentum", 100, 0., 100., 100, -100., 100.);
1396 hLikeVSEta_ =
new TH2F (name+
"_LikelihoodVSEta",
"likelihood vs muon pseudorapidity", 100, -4.,4., 100, -100., 100.);
1397 hLikeVSPhi_ =
new TH2F (name+
"_LikelihoodVSPhi",
"likelihood vs muon phi angle", 100, -3.2, 3.2, 100, -100., 100.);
1398 hLikeVSPt_prof_ =
new TProfile (name+
"_LikelihoodVSPt_prof",
"likelihood vs muon transverse momentum", 40, 0., 100., -1000., 1000. );
1399 hLikeVSEta_prof_ =
new TProfile (name+
"_LikelihoodVSEta_prof",
"likelihood vs muon pseudorapidity", 40, -4.,4., -1000., 1000. );
1400 hLikeVSPhi_prof_ =
new TProfile (name+
"_LikelihoodVSPhi_prof",
"likelihood vs muon phi angle", 40, -3.2, 3.2, -1000., 1000.);
1405 hLikeVSPt_ = (TH2F *) file->Get(name+
"_LikelihoodVSPt");
1406 hLikeVSEta_ = (TH2F *) file->Get(name+
"_LikelihoodVSEta");
1407 hLikeVSPhi_ = (TH2F *) file->Get(name+
"_LikelihoodVSPhi");
1408 hLikeVSPt_prof_ = (TProfile *) file->Get(name+
"_LikelihoodVSPt_prof");
1409 hLikeVSEta_prof_ = (TProfile *) file->Get(name+
"_LikelihoodVSEta_prof");
1410 hLikeVSPhi_prof_ = (TProfile *) file->Get(name+
"_LikelihoodVSPhi_prof");
1417 delete hLikeVSPt_prof_;
1418 delete hLikeVSEta_prof_;
1419 delete hLikeVSPhi_prof_;
1423 Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()), likeValue);
1426 void Fill(
const CLHEP::HepLorentzVector& momentum,
const double& likeValue)
override {
1427 hLikeVSPt_->Fill(momentum.perp(),likeValue);
1428 hLikeVSEta_->Fill(momentum.eta(),likeValue);
1429 hLikeVSPhi_->Fill(momentum.phi(),likeValue);
1430 hLikeVSPt_prof_->Fill(momentum.perp(),likeValue);
1431 hLikeVSEta_prof_->Fill(momentum.eta(),likeValue);
1432 hLikeVSPhi_prof_->Fill(momentum.phi(),likeValue);
1436 hLikeVSPt_->Write();
1437 hLikeVSEta_->Write();
1438 hLikeVSPhi_->Write();
1439 hLikeVSPt_prof_->Write();
1440 hLikeVSEta_prof_->Write();
1441 hLikeVSPhi_prof_->Write();
1445 hLikeVSPt_->Reset(
"ICE");
1446 hLikeVSEta_->Reset(
"ICE");
1447 hLikeVSPhi_->Reset(
"ICE");
1448 hLikeVSPt_prof_->Reset(
"ICE");
1449 hLikeVSEta_prof_->Reset(
"ICE");
1450 hLikeVSPhi_prof_->Reset(
"ICE");
1477 totBinsY_ = totBinsY;
1482 deltaX_ = xMax - xMin_;
1483 deltaY_ = yMax - yMin_;
1484 hReso_ =
new TH1F( name+
"_Reso",
"resolution", 1000, 0, 1 );
1485 hResoVSPtEta_ =
new TH2F( name+
"_ResoVSPtEta",
"resolution vs pt and #eta", totBinsX_, xMin_, xMax, totBinsY_, yMin_, yMax );
1487 resoVsPtEta_ =
new double*[totBinsX_];
1488 resoCount_ =
new int*[totBinsX_];
1489 for(
int i=0;
i<totBinsX_; ++
i ) {
1490 resoVsPtEta_[
i] =
new double[totBinsY_];
1491 resoCount_[
i] =
new int[totBinsY_];
1492 for(
int j=0; j<totBinsY_; ++j ) {
1493 resoVsPtEta_[
i][j] = 0;
1494 resoCount_[
i][j] = 0;
1497 hResoVSPt_prof_ =
new TProfile( name+
"_ResoVSPt_prof",
"resolution VS pt", totBinsX_, xMin_, xMax, yMin_, yMax);
1498 hResoVSPt_Bar_prof_ =
new TProfile( name+
"_ResoVSPt_Bar_prof",
"resolution VS pt Barrel", totBinsX_, xMin_, xMax, yMin_, yMax);
1499 hResoVSPt_Endc_17_prof_ =
new TProfile( name+
"_ResoVSPt_Endc_1.7_prof",
"resolution VS pt Endcap (1.4<eta<1.7)", totBinsX_, xMin_, xMax, yMin_, yMax);
1500 hResoVSPt_Endc_20_prof_ =
new TProfile( name+
"_ResoVSPt_Endc_2.0_prof",
"resolution VS pt Endcap (1.7<eta<2.0)", totBinsX_, xMin_, xMax, yMin_, yMax);
1501 hResoVSPt_Endc_24_prof_ =
new TProfile( name+
"_ResoVSPt_Endc_2.4_prof",
"resolution VS pt Endcap (2.0<eta<2.4)", totBinsX_, xMin_, xMax, yMin_, yMax);
1502 hResoVSPt_Ovlap_prof_ =
new TProfile( name+
"_ResoVSPt_Ovlap_prof",
"resolution VS pt Overlap", totBinsX_, xMin_, xMax, yMin_, yMax);
1503 hResoVSEta_prof_ =
new TProfile( name+
"_ResoVSEta_prof",
"resolution VS eta", totBinsY_, yMin_, yMax, 0, 1);
1505 hResoVSPhiPlus_prof_ =
new TProfile( name+
"_ResoVSPhiPlus_prof",
"resolution VS phi mu+", 16, -3.2, 3.2, 0, 1);
1506 hResoVSPhiMinus_prof_ =
new TProfile( name+
"_ResoVSPhiMinus_prof",
"resolution VS phi mu-", 16, -3.2, 3.2, 0, 1);
1507 hResoVSPhi_prof_ =
new TProfile( name+
"_ResoVSPhi_prof",
"resolution VS phi", 16, -3.2, 3.2, -1, 1);
1511 delete hResoVSPtEta_;
1513 for(
int i=0;
i<totBinsX_; ++
i ) {
1514 delete[] resoVsPtEta_[
i];
1515 delete[] resoCount_[
i];
1517 delete[] resoVsPtEta_;
1518 delete[] resoCount_;
1520 delete hResoVSPt_prof_;
1521 delete hResoVSPt_Bar_prof_;
1522 delete hResoVSPt_Endc_17_prof_;
1523 delete hResoVSPt_Endc_20_prof_;
1524 delete hResoVSPt_Endc_24_prof_;
1525 delete hResoVSPt_Ovlap_prof_;
1526 delete hResoVSEta_prof_;
1528 delete hResoVSPhiPlus_prof_;
1529 delete hResoVSPhiMinus_prof_;
1530 delete hResoVSPhi_prof_;
1533 if( resValue != resValue )
return;
1534 hReso_->Fill(resValue);
1537 int xIndex = getXindex(p4.Pt());
1538 int yIndex = getYindex(p4.Eta());
1539 if ( 0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_ ) {
1540 resoVsPtEta_[xIndex][yIndex] += resValue;
1549 resoCount_[xIndex][yIndex] += 1;
1552 hResoVSPt_prof_->Fill(p4.Pt(),resValue);
1553 if(fabs(p4.Eta())<=0.9)
1554 hResoVSPt_Bar_prof_->Fill(p4.Pt(),resValue);
1555 else if(fabs(p4.Eta())>0.9 && fabs(p4.Eta())<=1.4 )
1556 hResoVSPt_Ovlap_prof_->Fill(p4.Pt(),resValue);
1557 else if(fabs(p4.Eta())>1.4 && fabs(p4.Eta())<=1.7 )
1558 hResoVSPt_Endc_17_prof_->Fill(p4.Pt(),resValue);
1559 else if(fabs(p4.Eta())>1.7 && fabs(p4.Eta())<=2.0 )
1560 hResoVSPt_Endc_20_prof_->Fill(p4.Pt(),resValue);
1562 hResoVSPt_Endc_24_prof_->Fill(p4.Pt(),resValue);
1563 hResoVSEta_prof_->Fill(p4.Eta(),resValue);
1566 hResoVSPhiPlus_prof_->Fill(p4.Phi(),resValue);
1568 hResoVSPhiMinus_prof_->Fill(p4.Phi(),resValue);
1569 hResoVSPhi_prof_->Fill(p4.Phi(),resValue);
1574 if(histoDir_ !=
nullptr) histoDir_->cd();
1578 for(
int i=0;
i<totBinsX_; ++
i ) {
1579 for(
int j=0; j<totBinsY_; ++j ) {
1580 int N = resoCount_[
i][j];
1582 if( N != 0 ) hResoVSPtEta_->SetBinContent(
i+1, j+1, resoVsPtEta_[
i][j]/N );
1583 else hResoVSPtEta_->SetBinContent(
i+1, j+1, 0 );
1586 hResoVSPtEta_->Write();
1588 hResoVSPt_prof_->Write();
1589 hResoVSPt_Bar_prof_->Write();
1590 hResoVSPt_Endc_17_prof_->Write();
1591 hResoVSPt_Endc_20_prof_->Write();
1592 hResoVSPt_Endc_24_prof_->Write();
1593 hResoVSPt_Ovlap_prof_->Write();
1594 hResoVSEta_prof_->Write();
1596 hResoVSPhiMinus_prof_->Write();
1597 hResoVSPhiPlus_prof_->Write();
1598 hResoVSPhi_prof_->Write();
1600 TCanvas
canvas(TString(hResoVSPtEta_->GetName())+
"_canvas", TString(hResoVSPtEta_->GetTitle())+
" canvas", 1000, 800);
1603 hResoVSPtEta_->Draw(
"lego");
1605 hResoVSPtEta_->Draw(
"surf5");
1607 hResoVSPtEta_->Write();
1614 hResoVSPtEta_->Clear();
1615 hResoVSPt_prof_->Clear();
1616 hResoVSPt_Bar_prof_->Clear();
1617 hResoVSPt_Endc_17_prof_->Clear();
1618 hResoVSPt_Endc_20_prof_->Clear();
1619 hResoVSPt_Endc_24_prof_->Clear();
1620 hResoVSPt_Ovlap_prof_->Clear();
1621 hResoVSEta_prof_->Clear();
1623 hResoVSPhiPlus_prof_->Clear();
1624 hResoVSPhiMinus_prof_->Clear();
1625 hResoVSPhi_prof_->Clear();
1630 return int((x-xMin_)/deltaX_*totBinsX_);
1633 return int((y-yMin_)/deltaY_*totBinsY_);
1661 histoVarianceCheck_ =
new TH1D**[totBinsX_];
1662 for(
int i=0;
i<totBinsX_; ++
i ) {
1663 histoVarianceCheck_[
i] =
new TH1D*[totBinsY_];
1664 for(
int j=0; j<totBinsY_; ++j ) {
1665 std::stringstream namei;
1666 std::stringstream namej;
1669 histoVarianceCheck_[
i][j] =
new TH1D(name+
"_"+namei.str()+
"_"+namej.str(),
name, 100, 0., 1.);
1674 for(
int i=0;
i<totBinsX_; ++
i ) {
1675 for(
int j=0; j<totBinsY_; ++j ) {
1676 delete histoVarianceCheck_[
i][j];
1678 delete[] histoVarianceCheck_;
1680 delete[] histoVarianceCheck_;
1683 if( resValue != resValue )
return;
1685 int xIndex = getXindex(p4.Pt());
1686 int yIndex = getYindex(p4.Eta());
1688 if ( 0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_ ) {
1689 histoVarianceCheck_[xIndex][yIndex]->Fill(resValue);
1695 if(histoDir_ !=
nullptr) histoDir_->cd();
1696 for(
int xBin=0; xBin<totBinsX_; ++xBin ) {
1698 histoVarianceCheck_[xBin][
yBin]->Write();
1718 const int totBins,
const double &
xMin,
const double &
xMax,
1719 const double &
yMin,
const double &
yMax, TDirectory *
dir =
nullptr) :
1724 if(
dir_ !=
nullptr ) {
1725 dir2D_ = (TDirectory*)
dir_->Get(
"2D");
1726 if(dir2D_ ==
nullptr) dir2D_ =
dir_->mkdir(
"2D");
1727 diffDir_ = (TDirectory*)
dir_->Get(
"deltaXoverX");
1728 if(diffDir_ ==
nullptr) diffDir_ =
dir->mkdir(
"deltaXoverX");
1730 diffHisto_ =
new TProfile(name+
"_prof", title+
" profile", totBins, xMin, xMax, yMin, yMax);
1731 histo2D_ =
new TH2F(name+
"2D", title, totBins, xMin, xMax, 4000, yMin, yMax);
1732 resoHisto_ =
new TH1F(name, title, totBins, xMin, xMax);
1739 Int_t
Fill( Double_t
x, Double_t
y )
override {
1740 diffHisto_->Fill(x, y);
1741 histo2D_->Fill(x, y);
1751 unsigned int totBins = diffHisto_->GetNbinsX();
1753 for(
unsigned int iBin=1; iBin<=totBins; ++iBin ) {
1755 resoHisto_->SetBinContent( iBin, diffHisto_->GetBinError(iBin)*
sqrt(diffHisto_->GetBinEntries(iBin)) );
1758 resoHisto_->Write();
1759 if( diffDir_ !=
nullptr ) diffDir_->cd();
1760 diffHisto_->Write();
1761 if( dir2D_ !=
nullptr ) dir2D_->cd();
1791 void fill(
const double &
x,
const double &
y) {
1799 double meanX = sumX_/N_;
1800 double meanY = sumY_/N_;
1802 return (productXY_/N_ - meanX*meanY);
1822 const int totBinsX,
const double &
xMin,
const double &
xMax,
1823 const int totBinsY,
const double &
yMin,
const double &
yMax,
1824 TDirectory *
dir =
nullptr,
bool varianceCheck =
false ) :
1825 totBinsX_(totBinsX), totBinsY_(totBinsY),
1826 xMin_(xMin), deltaX_(xMax-xMin), yMin_(yMin), deltaY_(yMax-yMin),
1828 varianceCheck_(varianceCheck)
1832 histoCovariance_ =
new TH2D(name+
"Covariance", title+
" covariance", totBinsX, xMin, xMax, totBinsY, yMin, yMax);
1835 for(
int i=0;
i<totBinsX; ++
i ) {
1838 if( varianceCheck_ ) {
1839 histoVarianceCheck_ =
new TH1D**[totBinsX_];
1840 for(
int i=0;
i<totBinsX_; ++
i ) {
1841 histoVarianceCheck_[
i] =
new TH1D*[totBinsY_];
1842 for(
int j=0; j<totBinsY_; ++j ) {
1843 std::stringstream namei;
1844 std::stringstream namej;
1847 histoVarianceCheck_[
i][j] =
new TH1D(name+
"_"+namei.str()+
"_"+namej.str(),
name, 10000, -1, 1);
1856 histoDir_ = (TDirectory*)(inputFile->Get(dirName.Data()));
1857 if( histoDir_ ==
nullptr ) {
1858 std::cout <<
"Error: directory not found" << std::endl;
1861 histoCovariance_ = (TH2D*)(histoDir_->Get(name));
1862 totBinsX_ = histoCovariance_->GetNbinsX();
1863 xMin_ = histoCovariance_->GetXaxis()->GetBinLowEdge(1);
1864 deltaX_ = histoCovariance_->GetXaxis()->GetBinUpEdge(totBinsX_) - xMin_;
1865 totBinsY_ = histoCovariance_->GetNbinsY();
1866 yMin_ = histoCovariance_->GetYaxis()->GetBinLowEdge(1);
1867 deltaY_ = histoCovariance_->GetYaxis()->GetBinUpEdge(totBinsY_) - yMin_;
1871 delete histoCovariance_;
1873 for(
int i=0;
i<totBinsX_; ++
i) {
1874 delete[] covariances_[
i];
1876 delete[] covariances_;
1878 if( varianceCheck_ ) {
1879 for(
int i=0;
i<totBinsX_; ++
i ) {
1880 for(
int j=0; j<totBinsY_; ++j ) {
1881 delete histoVarianceCheck_[
i][j];
1883 delete[] histoVarianceCheck_[
i];
1885 delete[] histoVarianceCheck_;
1893 void Fill(
const double &
x,
const double &
y,
const double &
a,
const double &
b )
override {
1895 int xIndex = getXindex(x);
1896 int yIndex = getYindex(y);
1898 if ( 0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_ ) {
1900 covariances_[xIndex][yIndex].fill(a,b);
1902 if( varianceCheck_ ) histoVarianceCheck_[xIndex][yIndex]->Fill(a);
1906 using Histograms::Get;
1907 double Get(
const double &
x,
const double &
y )
const {
1909 int xIndex = getXindex(x);
1910 int yIndex = getYindex(y);
1912 if ( xIndex < 0 ) xIndex = 0;
1913 if ( xIndex >= totBinsX_ ) xIndex = totBinsX_-1;
1914 if ( yIndex < 0 ) yIndex = 0;
1915 if ( yIndex >= totBinsY_ ) yIndex = totBinsY_-1;
1916 return histoCovariance_->GetBinContent(xIndex+1, yIndex+1);
1921 std::cout <<
"writing: " << histoCovariance_->GetName() << std::endl;
1922 for(
int xBin=0; xBin<totBinsX_; ++xBin ) {
1924 double covariance = covariances_[xBin][
yBin].covariance();
1927 histoCovariance_->SetBinContent(xBin+1,
yBin+1, covariance);
1930 if( histoDir_ !=
nullptr ) histoDir_->cd();
1931 TCanvas
canvas(TString(histoCovariance_->GetName())+
"_canvas", TString(histoCovariance_->GetTitle())+
" canvas", 1000, 800);
1934 histoCovariance_->Draw(
"lego");
1936 histoCovariance_->Draw(
"surf5");
1938 histoCovariance_->Write();
1940 TDirectory * binsDir =
nullptr;
1941 if( varianceCheck_ ) {
1942 if ( histoDir_ !=
nullptr ) {
1944 if( binsDir ==
nullptr ) binsDir = histoDir_->mkdir(name_+
"Bins");
1947 for(
int xBin=0; xBin<totBinsX_; ++xBin ) {
1949 histoVarianceCheck_[xBin][
yBin]->Write();
1956 histoCovariance_->Clear();
1957 if( varianceCheck_ ) {
1958 for(
int i=0;
i<totBinsX_; ++
i ) {
1959 for(
int j=0; j<totBinsY_; ++j ) {
1960 histoVarianceCheck_[
i][j]->Clear();
1967 return int((x-xMin_)/deltaX_*totBinsX_);
1970 return int((y-yMin_)/deltaY_*totBinsY_);
1998 mapHisto_[name+
"Pt"] =
new HCovarianceVSxy(name+
"Pt_",
"Pt", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_,
true);
1999 mapHisto_[name+
"CotgTheta"] =
new HCovarianceVSxy(name+
"CotgTheta_",
"CotgTheta", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_,
true);
2000 mapHisto_[name+
"Phi"] =
new HCovarianceVSxy(name+
"Phi_",
"Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_,
true);
2002 mapHisto_[name+
"Pt-CotgTheta"] =
new HCovarianceVSxy(name+
"Pt_CotgTheta_",
"Pt-CotgTheta", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2003 mapHisto_[name+
"Pt-Phi"] =
new HCovarianceVSxy(name+
"Pt_Phi_",
"Pt-Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2004 mapHisto_[name+
"CotgTheta-Phi"] =
new HCovarianceVSxy(name+
"CotgTheta_Phi_",
"CotgTheta-Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2005 mapHisto_[name+
"Pt1-Pt2"] =
new HCovarianceVSxy(name+
"Pt1_Pt2_",
"Pt1-Pt2", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2006 mapHisto_[name+
"CotgTheta1-CotgTheta2"] =
new HCovarianceVSxy(name+
"CotgTheta1_CotgTheta2_",
"CotgTheta1-CotgTheta2", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2007 mapHisto_[name+
"Phi1-Phi2"] =
new HCovarianceVSxy(name+
"Phi1_Phi2_",
"Phi1-Phi2", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2008 mapHisto_[name+
"Pt12-CotgTheta21"] =
new HCovarianceVSxy(name+
"Pt12_CotgTheta21_",
"Pt12-CotgTheta21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2009 mapHisto_[name+
"Pt12-Phi21"] =
new HCovarianceVSxy(name+
"Pt12_Phi21_",
"Pt12-Phi21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2010 mapHisto_[name+
"CotgTheta12-Phi21"] =
new HCovarianceVSxy(name+
"CotgTheta12_Phi21_",
"CotgTheta12-Phi21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2017 TFile *
inputFile =
new TFile(inputFileName,
"READ");
2021 mapHisto_[name_+
"Pt"] =
new HCovarianceVSxy(inputFile, name_+
"Pt_"+name_, name_);
2022 mapHisto_[name_+
"CotgTheta"] =
new HCovarianceVSxy(inputFile, name_+
"CotgTheta_"+name_, name_);
2023 mapHisto_[name_+
"Phi"] =
new HCovarianceVSxy(inputFile, name_+
"Phi_"+name_, name_);
2025 mapHisto_[name_+
"Pt-CotgTheta"] =
new HCovarianceVSxy(inputFile, name_+
"Pt_CotgTheta_"+name_, name_);
2026 mapHisto_[name_+
"Pt-Phi"] =
new HCovarianceVSxy(inputFile, name_+
"Pt_Phi_"+name_, name_);
2027 mapHisto_[name_+
"CotgTheta-Phi"] =
new HCovarianceVSxy(inputFile, name_+
"CotgTheta_Phi_"+name_, name_);
2028 mapHisto_[name_+
"Pt1-Pt2"] =
new HCovarianceVSxy(inputFile, name_+
"Pt1_Pt2_"+name_, name_);
2029 mapHisto_[name_+
"CotgTheta1-CotgTheta2"] =
new HCovarianceVSxy(inputFile, name_+
"CotgTheta1_CotgTheta2_"+name_, name_);
2030 mapHisto_[name_+
"Phi1-Phi2"] =
new HCovarianceVSxy(inputFile, name_+
"Phi1_Phi2_"+name_, name_);
2031 mapHisto_[name_+
"Pt12-CotgTheta21"] =
new HCovarianceVSxy(inputFile, name_+
"Pt12_CotgTheta21_"+name_, name_);
2032 mapHisto_[name_+
"Pt12-Phi21"] =
new HCovarianceVSxy(inputFile, name_+
"Pt12_Phi21_"+name_, name_);
2033 mapHisto_[name_+
"CotgTheta12-Phi21"] =
new HCovarianceVSxy(inputFile, name_+
"CotgTheta12_Phi21_"+name_, name_);
2037 for (std::map<TString, HCovarianceVSxy*>::const_iterator
histo=mapHisto_.begin();
2039 delete (*histo).second;
2044 return mapHisto_[name_+covarianceName]->Get(recoP1.pt(), recoP1.eta());
2052 double pt1 = recoP1.pt();
2053 double eta1 = recoP1.eta();
2054 double pt2 = recoP2.pt();
2055 double eta2 = recoP2.eta();
2057 double diffPt1 = (pt1 - genP1.pt())/genP1.pt();
2058 double diffPt2 = (pt2 - genP2.pt())/genP2.pt();
2060 double genTheta1 = genP1.theta();
2061 double genTheta2 = genP2.theta();
2062 double recoTheta1 = recoP1.theta();
2063 double recoTheta2 = recoP2.theta();
2065 double genCotgTheta1 = TMath::Cos(genTheta1)/(TMath::Sin(genTheta1));
2066 double genCotgTheta2 = TMath::Cos(genTheta2)/(TMath::Sin(genTheta2));
2067 double recoCotgTheta1 = TMath::Cos(recoTheta1)/(TMath::Sin(recoTheta1));
2068 double recoCotgTheta2 = TMath::Cos(recoTheta2)/(TMath::Sin(recoTheta2));
2072 double diffCotgTheta1 = recoCotgTheta1 - genCotgTheta1;
2073 double diffCotgTheta2 = recoCotgTheta2 - genCotgTheta2;
2081 mapHisto_[name_+
"Pt"]->Fill(pt1, eta1, diffPt1, diffPt1);
2082 mapHisto_[name_+
"Pt"]->Fill(pt2, eta2, diffPt2, diffPt2);
2083 mapHisto_[name_+
"CotgTheta"]->Fill(pt1, eta1, diffCotgTheta1, diffCotgTheta1);
2084 mapHisto_[name_+
"CotgTheta"]->Fill(pt2, eta2, diffCotgTheta2, diffCotgTheta2);
2085 mapHisto_[name_+
"Phi"]->Fill(pt1, eta1, diffPhi1, diffPhi1);
2086 mapHisto_[name_+
"Phi"]->Fill(pt2, eta2, diffPhi2, diffPhi2);
2089 mapHisto_[name_+
"Pt-CotgTheta"]->Fill(pt1, eta1, diffPt1, diffCotgTheta1 );
2090 mapHisto_[name_+
"Pt-CotgTheta"]->Fill(pt2, eta2, diffPt2, diffCotgTheta2 );
2091 mapHisto_[name_+
"Pt-Phi"]->Fill(pt1, eta1, diffPt1, diffPhi1);
2092 mapHisto_[name_+
"Pt-Phi"]->Fill(pt2, eta2, diffPt2, diffPhi2);
2093 mapHisto_[name_+
"CotgTheta-Phi"]->Fill(pt1, eta1, diffCotgTheta1, diffPhi1);
2094 mapHisto_[name_+
"CotgTheta-Phi"]->Fill(pt2, eta2, diffCotgTheta2, diffPhi2);
2099 mapHisto_[name_+
"Pt1-Pt2"]->Fill(pt1, eta1, diffPt1, diffPt2);
2100 mapHisto_[name_+
"Pt1-Pt2"]->Fill(pt2, eta2, diffPt1, diffPt2);
2101 mapHisto_[name_+
"CotgTheta1-CotgTheta2"]->Fill(pt1, eta1, diffCotgTheta1, diffCotgTheta2);
2102 mapHisto_[name_+
"CotgTheta1-CotgTheta2"]->Fill(pt2, eta2, diffCotgTheta1, diffCotgTheta2);
2103 mapHisto_[name_+
"Phi1-Phi2"]->Fill(pt1, eta1, diffPhi1, diffPhi2);
2104 mapHisto_[name_+
"Phi1-Phi2"]->Fill(pt2, eta2, diffPhi1, diffPhi2);
2109 mapHisto_[name_+
"Pt12-CotgTheta21"]->Fill(pt1, eta1, diffPt1, diffCotgTheta2);
2110 mapHisto_[name_+
"Pt12-CotgTheta21"]->Fill(pt2, eta2, diffPt2, diffCotgTheta1);
2111 mapHisto_[name_+
"Pt12-Phi21"]->Fill(pt1, eta1, diffPt1, diffPhi2);
2112 mapHisto_[name_+
"Pt12-Phi21"]->Fill(pt2, eta2, diffPt2, diffPhi1);
2113 mapHisto_[name_+
"CotgTheta12-Phi21"]->Fill(pt1, eta1, diffCotgTheta1, diffPhi2);
2114 mapHisto_[name_+
"CotgTheta12-Phi21"]->Fill(pt2, eta2, diffCotgTheta2, diffPhi1);
2119 for (std::map<TString, HCovarianceVSxy*>::const_iterator
histo=mapHisto_.begin();
2121 (*histo).second->Write();
2126 for (std::map<TString, HCovarianceVSxy*>::const_iterator
histo=mapHisto_.begin();
2128 (*histo).second->Clear();
2150 nameSuffix_[0] =
"Plus";
2151 nameSuffix_[1] =
"Minus";
2152 TString titleSuffix[] = {
" for mu+",
" for mu-"};
2154 mapHisto_[
name] =
new TH1F (name,
"#Delta M/M", 4000, -1, 1);
2155 mapHisto_[name+
"VSPairPt"] =
new HResolution (name+
"VSPairPt",
"resolution VS pt of the pair", 100, 0, 200, -1, 1, histoDir_);
2156 mapHisto_[name+
"VSPairDeltaEta"] =
new HResolution (name+
"VSPairDeltaEta",
"resolution VS #Delta#eta of the pair", 100, -0.1, 6.2, -1, 1, histoDir_);
2157 mapHisto_[name+
"VSPairDeltaPhi"] =
new HResolution (name+
"VSPairDeltaPhi",
"resolution VS #Delta#phi of the pair", 100, -0.1, 3.2, -1, 1, histoDir_);
2159 for(
int i=0;
i<2; ++
i ) {
2160 mapHisto_[name+
"VSPt"+nameSuffix_[
i]] =
new HResolution (name+
"VSPt"+nameSuffix_[
i],
"resolution VS pt"+titleSuffix[i], 100, 0, 200, -1, 1, histoDir_);
2161 mapHisto_[name+
"VSEta"+nameSuffix_[
i]] =
new HResolution (name+
"VSEta"+nameSuffix_[i],
"resolution VS #eta"+titleSuffix[i], 100, -3, 3, -1, 1, histoDir_);
2162 mapHisto_[name+
"VSPhi"+nameSuffix_[
i]] =
new HResolution (name+
"VSPhi"+nameSuffix_[i],
"resolution VS #phi"+titleSuffix[i], 100, -3.2, 3.2, -1, 1, histoDir_);
2166 muMinus.reset(
new HDelta(
"muMinus") );
2167 muPlus.reset(
new HDelta(
"muPlus") );
2171 for (std::map<TString, TH1*>::const_iterator
histo=mapHisto_.begin();
2173 delete (*histo).second;
2181 const double & recoMass,
const double & genMass )
override {
2182 muMinus->Fill(recoP1, genP1);
2183 muPlus->Fill(recoP2, genP2);
2185 Fill( recoP1, charge1, recoP2, charge2, recoMass, genMass );
2192 const double & recoMass,
const double & genMass )
override {
2194 if ( charge1 == charge2 )
std::cout <<
"Error: must get two opposite charge particles" << std::endl;
2196 double massRes = (recoMass - genMass)/genMass;
2199 double pairPt = recoPair.Pt();
2201 double recoPt[2] = {recoP1.Pt(), recoP2.Pt()};
2202 double recoEta[2] = {recoP1.Eta(), recoP2.Eta()};
2203 double recoPhi[2] = {recoP1.Phi(), recoP2.Phi()};
2212 if ( charge1 > 0 ) {
2217 double pairDeltaEta = fabs(recoEta[0] - recoEta[1]);
2220 mapHisto_[name_]->Fill(massRes);
2221 mapHisto_[name_+
"VSPairPt"]->Fill(pairPt, massRes);
2222 mapHisto_[name_+
"VSPairDeltaEta"]->Fill(pairDeltaEta, massRes);
2223 mapHisto_[name_+
"VSPairDeltaPhi"]->Fill(pairDeltaPhi, massRes);
2228 for(
int i=0;
i<2; ++
i ) {
2230 mapHisto_[name_+
"VSPt"+nameSuffix_[
i]]->Fill(recoPt[index], massRes);
2231 mapHisto_[name_+
"VSEta"+nameSuffix_[
i]]->Fill(recoEta[index], massRes);
2232 mapHisto_[name_+
"VSPhi"+nameSuffix_[
i]]->Fill(recoPhi[index], massRes);
2238 for (std::map<TString, TH1*>::const_iterator
histo=mapHisto_.begin();
2240 (*histo).second->Write();
2243 (histoDir_->mkdir(
"singleMuonsVSgen"))->
cd();
2249 for (std::map<TString, TH1*>::const_iterator
histo=mapHisto_.begin();
2251 (*histo).second->Clear();
2273 TString nameSuffix_[2];
static double deltaPhiNoFabs(const double &phi1, const double &phi2)
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance c...
int getXindex(const double &x) const
HTH1D(TFile *outputFile, const TString &name, const TString &title, const int xBins, const double &xMin, const double &xMax)
void Fill(const reco::Particle::LorentzVector &recoP1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const reco::Particle::LorentzVector &genP2) override
TProfile * hMassVSPhi_prof_
virtual void Fill(const CLHEP::HepLorentzVector &momentum, const int charge, const double &weight=1.)
HLikelihoodVSPart(const TString &name)
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
~HFunctionResolution() override
HPartVSEta(const TString &name, const double &minMass=0., const double &maxMass=100., const double &maxPt=100.)
Histograms(TFile *outputFile, const TString &name)
virtual void SetYTitle(const TString &title)
HParticle(TFile *outputFile, const TString &name, const double &minMass=0., const double &maxMass=200., const double &maxPt=100.)
Constructor that puts the histograms inside a TDirectory.
HPartVSPt(const TString &name)
std::unique_ptr< HDelta > muPlus
~HCovarianceVSParts() override
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
TProfile * hPtVSEta_prof_
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const CLHEP::HepLorentzVector &momentumRes, const double &weight=1.)
void Fill(const reco::Particle::LorentzVector &p4, const int charge, const double &weight=1.) override
virtual double Get(const reco::Particle::LorentzVector &recoP1, const TString &covarianceName)
HMassResolutionVSPart(TFile *outputFile, const TString &name)
~HLikelihoodVSPart() override
double Get(const double &x, const double &y) const
HLikelihoodVSPart(const TString &name, TFile *file)
~HMassVSPartProfile() override
virtual void SetXTitle(const TString &title)
virtual void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const reco::Particle::LorentzVector &p4Res, const double &weight=1.)
std::unique_ptr< HDelta > muMinus
TH3F * hMassVSEtaPhiMinus_
void Fill(const double &x, const double &y, const double &a, const double &b) override
TH2F * hResoVSPt_Endc_17_
TProfile * hMassVSPt_prof_
A wrapper for the TH1D histogram to allow it to be put inside the same map as all the other classes i...
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.)
TH3F * hMassVSPhiPlusPhiMinus_
~HCovarianceVSxy() override
void Fill(const reco::Particle::LorentzVector &recoP1, const int charge1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const int charge2, const reco::Particle::LorentzVector &genP2, const double &recoMass, const double &genMass) override
void Fill(unsigned int number) override
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
TProfile * hResoVSPhi_prof_
TProfile2D * hMassVSPhiPlus_
TProfile * hCurvVSEta_prof_
int getYindex(const double &y) const
void Fill(const CLHEP::HepLorentzVector &p1, const reco::Particle::LorentzVector &p2) override
virtual void Fill(const reco::Particle::LorentzVector &recoP1, const int charge1, const reco::Particle::LorentzVector &recoP2, const int charge2, const double &recoMass, const double &genMass)
TProfile * hLikeVSEta_prof_
TH1D *** histoVarianceCheck_
~HResolutionVSPart() override
TProfile * hPtVSPhi_prof_
A wrapper for the TProfile histogram to allow it to be put inside the same map as all the other class...
HParticle(const TString &name, TFile *file)
TProfile * hMassVSCosThetaCS_prof
virtual void SetYTitle(const TString &title)
HCovarianceVSxy(const TString &name, const TString &title, const int totBinsX, const double &xMin, const double &xMax, const int totBinsY, const double &yMin, const double &yMax, TDirectory *dir=0, bool varianceCheck=false)
TProfile * hResoVSPt_Endc_17_prof_
virtual void Fill(const double &x, const double &y)
TProfile * hResoVSPt_Endc_20_prof_
TH3F * hMassVSEtaPlusEtaMinus_
std::map< TString, TH1 * > mapHisto_
void Fill(const double &x, const double &y) override
virtual void SetYTitle(const TString &title)
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
TProfile * hResoVSPt_Ovlap_prof_
TProfile * hCurvVsEtaPos_
TProfile * hResoVSPt_Endc_17_prof_
void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2) override
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
TProfile * hResoVSPt_Bar_prof_
HDelta(const TString &name, TFile *file)
TProfile * hResoVSPt_prof_
virtual void SetXTitle(const TString &title)
std::map< TString, HCovarianceVSxy * > mapHisto_
void fill(const double &x, const double &y)
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
TProfile * hCurvVsEtaNeg_
void Fill(const double &x, const double &y) override
HMassVSPartProfile(const TString &name, const double &minMass=0., const double &maxMass=150., const double maxPt=100.)
Int_t Fill(Double_t x, Double_t y) override
virtual void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2)
void Fill(const CLHEP::HepLorentzVector &momentum, const int charge, const double &weight=1.) override
HFunctionResolutionVarianceCheck(TFile *outputFile, const TString &name, const double ptMax=200)
TProfile * hLikeVSPhi_prof_
HParticle(const TString &name, const double &minMass=0., const double &maxMass=200., const double &maxPt=100.)
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const int charge, const double &weight=1.) override
TProfile * hResoVSPt_Ovlap_prof_
TProfile * hCurvVsPhiPos_
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.) override
virtual void Fill(const unsigned int number)
HMassVSPart(const TString &name, TFile *file)
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &genValue, const double recValue, const int charge)
TProfile * hResoVSPt_Endc_24_prof_
TProfile * hCurvVsPhiNeg_
virtual void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2, const int charge, const double &weight=1.)
virtual void SetXTitle(const TString &title)
virtual void Fill(const reco::Particle::LorentzVector &p4, const int charge, const double &weight=1.)
TH2F * hMassVSCosThetaCS_
void Fill(const reco::Particle::LorentzVector &p4, const double &likeValue) override
~HMassResolutionVSPart() override
TH2F * hResoVSPt_Endc_24_
Covariance ** covariances_
TH3F * hMassVSEtaPhiPlus_
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
HCovarianceVSParts(TFile *outputFile, const TString &name, const double &ptMax)
TProfile * hResoVSPhi_prof_
virtual void Fill(const reco::Particle::LorentzVector &recoP1, const int charge1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const int charge2, const reco::Particle::LorentzVector &genP2, const double &recoMass, const double &genMass)
TProfile * hResoVSPt_prof_
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const reco::Particle::LorentzVector &p4Res, const double &weight=1.) override
HCovarianceVSxy(TFile *inputFile, const TString &name, const TString &dirName)
Contructor to read histograms from file.
TProfile * hMassVSEta_prof_
DecomposeProduct< arg, typename Div::arg > D
TProfile * hResoVSPhiMinus_prof_
virtual void Fill(const reco::Particle::LorentzVector &recoP1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const reco::Particle::LorentzVector &genP2)
virtual void SetWeight(double weight)
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge)
TH2F * hResoVSPt_Endc_20_
TH2F * hMassVSEtaPlusMinusDiff_
virtual void Fill(const CLHEP::HepLorentzVector &p, const double &likeValue)
int getXindex(const double &x) const
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
TProfile * hResoVSPt_Endc_24_prof_
~HFunctionResolutionVarianceCheck() override
HTH2D(TFile *outputFile, const TString &name, const TString &title, const TString &dirName, const int xBins, const double &xMin, const double &xMax, const int yBins, const double &yMin, const double &yMax)
A wrapper for the TH2D histogram to allow it to be put inside the same map as all the other classes i...
HMassVSPartProfile(const TString &name, TFile *file)
HDelta(TFile *outputFile, const TString &name)
virtual TString GetName()
void Fill(const double &x, const double &y) override
int getYindex(const double &y) const
HTProfile(TFile *outputFile, const TString &name, const TString &title, const int xBins, const double &xMin, const double &xMax, const double &yMin, const double &yMax)
TH1D *** histoVarianceCheck_
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.) override
TProfile * hResoVSEta_prof_
TH2F * hMassVSPhiPlusMinusDiff_
HDelta(const TString &name)
void Fill(const CLHEP::HepLorentzVector &momentum, const double &likeValue) override
TProfile * hResoVSPt_Bar_prof_
TProfile2D * hMassVSPhiMinus_
TProfile * hLikeVSPt_prof_
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
static double deltaPhi(const double &phi1, const double &phi2)
TProfile * hResoVSPt_Endc_20_prof_
HFunctionResolution(TFile *outputFile, const TString &name, const double &ptMax=100, const int totBinsY=300)
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2)
HPartVSPhi(const TString &name)
void Fill(const reco::Particle::LorentzVector &recoP1, const int charge1, const reco::Particle::LorentzVector &recoP2, const int charge2, const double &recoMass, const double &genMass) override
HResolutionVSPart(TFile *outputFile, const TString &name, const double maxPt=100, const double &yMinEta=0., const double &yMaxEta=2., const double &yMinPt=0., const double &yMaxPt=2., const bool doProfiles=false)
TProfile * hResoVSEta_prof_
TProfile * hResoVSPhiPlus_prof_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2) override
A set of histograms for resolution.
Histograms(const TString &name)
HCovarianceVSParts(const TString &inputFileName, const TString &name)
Constructor reading the histograms from file.
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const int charge, const double &weight=1.) override
double Get(const reco::Particle::LorentzVector &recoP1, const TString &covarianceName) override
virtual void Fill(const CLHEP::HepLorentzVector &p1, const reco::Particle::LorentzVector &p2)
Power< A, B >::type pow(const A &a, const B &b)
HResolutionVSPart(const TString &name, TFile *file)
HResolution(const TString &name, const TString &title, const int totBins, const double &xMin, const double &xMax, const double &yMin, const double &yMax, TDirectory *dir=0)
HMassVSPart(const TString &name, const double &minMass=0., const double &maxMass=150., const double maxPt=100.)
virtual void Fill(const double &x, const double &y, const double &a, const double &b)
Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0) override
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const CLHEP::HepLorentzVector &momentumRes, const double &weight=1.) override
Used to fill 2D histograms for comparison of opposite charge muons quantities.
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.)