CMS 3D CMS Logo

Histograms.h
Go to the documentation of this file.
1 #ifndef Histograms_H
2 #define Histograms_H
3 
11 #include <CLHEP/Vector/LorentzVector.h>
14 #include "MuScleFitUtils.h"
15 
16 #include "TH1D.h"
17 #include "TH1F.h"
18 #include "TH2F.h"
19 #include "TH3F.h"
20 #include "TFile.h"
21 #include "TString.h"
22 #include "TProfile.h"
23 #include "TProfile2D.h"
24 #include "TF1.h"
25 #include "TGraphErrors.h"
26 #include "TFile.h"
27 #include "TSystem.h"
28 #include "TCanvas.h"
29 
30 #include "TLorentzVector.h"
31 #include <vector>
32 #include <string>
33 #include <iostream>
34 #include "TMath.h"
35 
36 class Histograms {
37 public:
38 
39  // Constructor
40  // -----------
42  Histograms( const TString & name ) : theWeight_(1), name_(name), histoDir_(nullptr) {}
43  Histograms( TFile * outputFile, const TString & name ) :
44  theWeight_(1),
45  name_(name),
46  outputFile_(outputFile),
47  histoDir_( outputFile->GetDirectory(name) )
48  {
49  if( histoDir_ == nullptr ) {
50  histoDir_ = outputFile->mkdir(name);
51  }
52  histoDir_->cd();
53  }
54 
55  // Destructor
56  // ----------
57  virtual ~Histograms() {};
58 
59  // Operations
60  // ----------
61  // virtual void Fill( const reco::Particle::LorentzVector & p4 ) {};
62  // virtual void Fill( const CLHEP::HepLorentzVector & momentum ) {};
64  virtual void Fill( const reco::Particle::LorentzVector & p1, const reco::Particle::LorentzVector & p2, const int charge, const double & weight = 1.) {};
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.) {};
67  virtual void Fill( const CLHEP::HepLorentzVector & p1, const reco::Particle::LorentzVector & p2 ) {};
68  virtual void Fill( const reco::Particle::LorentzVector & p4, const double & weight = 1. ) {};
69 
70  //virtual void Fill( const CLHEP::HepLorentzVector & momentum, const double & weight ) {};
71  //------
72  virtual void Fill( const reco::Particle::LorentzVector & p4, const int charge, const double & weight = 1. ) {};
73 
74  virtual void Fill( const CLHEP::HepLorentzVector & momentum, const int charge, const double & weight =1.){};
75  //------
76  // virtual void Fill( const reco::Particle::LorentzVector & p4, const double & likeValue ) {};
77  virtual void Fill( const reco::Particle::LorentzVector & p4, const double & resValue, const int charge ) {};
78  virtual void Fill( const reco::Particle::LorentzVector & p4, const double & genValue, const double recValue, const int charge ) {};
79  virtual void Fill( const CLHEP::HepLorentzVector & p, const double & likeValue ) {};
80  virtual void Fill( const unsigned int number ) {};
81  virtual void Fill( const reco::Particle::LorentzVector & recoP1, const int charge1,
82  const reco::Particle::LorentzVector & genP1,
83  const reco::Particle::LorentzVector & recoP2, const int charge2,
84  const reco::Particle::LorentzVector & genP2,
85  const double & recoMass, const double & genMass ) {};
86  virtual void Fill( const reco::Particle::LorentzVector & recoP1, const int charge1,
87  // const reco::Particle::LorentzVector & genP1,
88  const reco::Particle::LorentzVector & recoP2, const int charge2,
89  // const reco::Particle::LorentzVector & genP2,
90  const double & recoMass, const double & genMass ) {};
91  virtual void Fill( const reco::Particle::LorentzVector & recoP1,
92  const reco::Particle::LorentzVector & genP1,
93  const reco::Particle::LorentzVector & recoP2,
94  const reco::Particle::LorentzVector & genP2 ) {};
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 ) {};
97  virtual void Fill(const reco::Particle::LorentzVector & p41,
99  const reco::Particle::LorentzVector & p4Res,
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.) {};
105 
106 
107  virtual double Get( const reco::Particle::LorentzVector & recoP1, const TString & covarianceName ) { return 0.; };
108 
109  virtual void Write() = 0;
110  virtual void Clear() = 0;
111 
112  virtual void SetWeight (double weight) {
113  theWeight_ = weight;
114  }
115 
116  virtual TString GetName() {
117  return name_;
118  }
119 
120 protected:
121  double theWeight_;
122  TString name_;
123  TFile * outputFile_;
124  TDirectory * histoDir_;
125 
126 private:
127 
128 };
129 
131 class HTH2D : public Histograms
132 {
133 public:
134  HTH2D( TFile * outputFile, const TString & name, const TString & title, const TString & dirName,
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) ) {}
139  ~HTH2D() override {
140  delete tH2d_;
141  delete tProfile_;
142  }
143  void Fill( const double & x, const double & y ) override {
144  tH2d_->Fill(x,y);
145  tProfile_->Fill(x,y);
146  }
147  void Write() override {
148  if(histoDir_ != nullptr) histoDir_->cd();
149  tH2d_->Write();
150  tProfile_->Write();
151  }
152  void Clear() override {
153  tH2d_->Clear();
154  tProfile_->Clear();
155  }
156  virtual void SetXTitle(const TString & title) {
157  tH2d_->GetXaxis()->SetTitle(title);
158  tProfile_->GetXaxis()->SetTitle(title);
159  }
160  virtual void SetYTitle(const TString & title) {
161  tH2d_->GetYaxis()->SetTitle(title);
162  tProfile_->GetYaxis()->SetTitle(title);
163  }
164  TH2D * operator->() { return tH2d_; }
165  TProfile * getProfile() { return tProfile_; }
166 protected:
167  TH2D * tH2d_;
168  TProfile * tProfile_;
169 };
170 
172 class HTH1D : public Histograms
173 {
174 public:
175  HTH1D( TFile * outputFile, const TString & name, const TString & title,
176  const int xBins, const double & xMin, const double & xMax ) : Histograms(outputFile, name),
177  tH1D_( new TH1D(name, title, xBins, xMin, xMax) ) {}
178  ~HTH1D() override {
179  delete tH1D_;
180  }
181  void Fill( const double & x, const double & y ) override {
182  tH1D_->Fill(x, y);
183  }
184  void Write() override {
185  if(histoDir_ != nullptr) histoDir_->cd();
186  tH1D_->Write();
187  }
188  void Clear() override {
189  tH1D_->Clear();
190  }
191  virtual void SetXTitle(const TString & title) {
192  tH1D_->GetXaxis()->SetTitle(title);
193  }
194  virtual void SetYTitle(const TString & title) {
195  tH1D_->GetYaxis()->SetTitle(title);
196  }
197  TH1D * operator->() { return tH1D_; }
198 protected:
199  TH1D * tH1D_;
200 };
201 
203 class HTProfile : public Histograms
204 {
205 public:
206  HTProfile( TFile * outputFile, const TString & name, const TString & 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) ) {}
210  ~HTProfile() override {
211  delete tProfile_;
212  }
213  void Fill( const double & x, const double & y ) override {
214  tProfile_->Fill(x,y);
215  }
216  void Write() override {
217  if(histoDir_ != nullptr) histoDir_->cd();
218  tProfile_->Write();
219  }
220  void Clear() override {
221  tProfile_->Clear();
222  }
223  virtual void SetXTitle(const TString & title) {
224  tProfile_->GetXaxis()->SetTitle(title);
225  }
226  virtual void SetYTitle(const TString & title) {
227  tProfile_->GetYaxis()->SetTitle(title);
228  }
229  TProfile * operator->() { return tProfile_; }
230 protected:
231  TProfile * tProfile_;
232 };
233 
234 // -----------------------------------------------------
235 // A set of histograms of particle kinematical variables
236 // -----------------------------------------------------
237 class HParticle : public Histograms {
238  public:
239  HParticle( const TString & name, const double & minMass = 0., const double & maxMass = 200., const double & maxPt = 100. ) :
240  Histograms(name),
241  // Kinematical variables
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) ),
244 
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.) ),
249 
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) ),
252 
253 
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) )
258  {}
259 
261  HParticle( TFile* outputFile, const TString & name, const double & minMass = 0., const double & maxMass = 200., const double & maxPt = 100. ) :
262  Histograms(outputFile, name)
263  {
264  // Kinematical variables
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);
267 
268  hPtVsEta_ = new TH2F (name+"_PtVsEta", "transverse momentum vs #eta", 100, 0, maxPt, 100, -3.0, 3.0);
269 
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.);
274 
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);
277 
278 
279  //hPtVSPhi_prof_ = new TProfile (name+"_PtVSPhi_prof", "pt vs phi angle",12, -3.2, 3.2, 0, 200);
280 
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);
285  }
286 
287  HParticle( const TString & name, TFile* file ) :
288  Histograms(name),
289  hPt_( (TH1F *) file->Get(name_+"_Pt") ),
290  hPtVsEta_( (TH2F *) file->Get(name_+"_PtVsEta") ),
291 
292 
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") ),
297 
298  hPtVsPhiNeg_( (TProfile *) file->Get(name_+"_PtVsPhiNeg") ),
299  hPtVsPhiPos_( (TProfile *) file->Get(name_+"_PtVsPhiPos") ),
300 
301  hEta_( (TH1F *) file->Get(name_+"_Eta") ),
302  hPhi_( (TH1F *) file->Get(name_+"_Phi") ),
303  hMass_( (TH1F *) file->Get(name_+"_Mass") ),
304  //hMass_fine_ = (TH1F *) file->Get(name_+"_Mass_fine");
305  hNumber_( (TH1F *) file->Get(name_+"_Number") )
306  {}
307 
308  ~HParticle() override
309  {
310  delete hPt_;
311  delete hPtVsEta_;
312 
313  delete hCurvVsEtaNeg_;
314  delete hCurvVsEtaPos_;
315  delete hCurvVsPhiNeg_;
316  delete hCurvVsPhiPos_;
317 
318  delete hPtVsPhiNeg_;
319  delete hPtVsPhiPos_;
320 
321  delete hEta_;
322  delete hPhi_;
323  delete hMass_;
324  // delete hMass_fine_;
325  delete hNumber_;
326  }
327 
328  void Fill( const reco::Particle::LorentzVector & p4, const int charge, const double & weight = 1. ) override
329  {
330  Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()),charge, weight);
331  }
332 
333  void Fill( const CLHEP::HepLorentzVector & momentum, const int charge, const double & weight =1.) override
334  {
335  hPt_->Fill(momentum.perp(), weight);
336  hPtVsEta_->Fill(momentum.perp(), momentum.eta(), weight);
337 
338  // std::cout<< "charge-> " <<charge<<std::endl;
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);
343 
344  if(charge<0)hPtVsPhiNeg_->Fill( momentum.phi(),momentum.perp(),weight);
345  if(charge>0)hPtVsPhiPos_->Fill( momentum.phi(),momentum.perp(),weight);
346 
347  hEta_->Fill(momentum.eta(), weight);
348  hPhi_->Fill(momentum.phi(), weight);
349  hMass_->Fill(momentum.m(), weight);
350  //hMass_fine_->Fill(momentum.m(), weight);
351  }
352 
353 
354 
355  void Fill( unsigned int number ) override
356  {
357  hNumber_->Fill(number);
358  }
359 
360  void Write() override
361  {
362  if(histoDir_ != nullptr) histoDir_->cd();
363  hPt_->Write();
364  hPtVsEta_->Write();
365 
366  hCurvVsEtaNeg_->Write();
367  hCurvVsEtaPos_->Write();
368  hCurvVsPhiNeg_->Write();
369  hCurvVsPhiPos_->Write();
370 
371  hPtVsPhiNeg_->Write();
372  hPtVsPhiPos_->Write();
373 
374  hEta_->Write();
375  hPhi_->Write();
376  hMass_->Write();
377  //hMass_fine_->Write();
378  hNumber_->Write();
379  }
380 
381  void Clear() override
382  {
383  hPt_->Clear();
384  hPtVsEta_->Clear();
385 
386  hCurvVsEtaNeg_->Clear();
387  hCurvVsEtaPos_->Clear();
388  hCurvVsPhiNeg_->Clear();
389  hCurvVsPhiPos_->Clear();
390 
391  hPtVsPhiNeg_->Clear();
392  hPtVsPhiPos_->Clear();
393 
394  hEta_->Clear();
395  hPhi_->Clear();
396  hMass_->Clear();
397  //hMass_fine_->Clear();
398  hNumber_->Clear();
399  }
400 
401  protected:
402  TH1F* hPt_;
403  TH2F* hPtVsEta_;
404 
405  TProfile* hCurvVsEtaNeg_;
406  TProfile* hCurvVsEtaPos_;
407  TProfile* hCurvVsPhiNeg_;
408  TProfile* hCurvVsPhiPos_;
409 
410  TProfile* hPtVsPhiNeg_;
411  TProfile* hPtVsPhiPos_;
412 
413  TH1F* hEta_;
414  TH1F* hPhi_;
415  TH1F* hMass_;
416  //TH1F* hMass_fine_;
417  TH1F* hNumber_;
418 };
419 
420 // ---------------------------------------------------
421 // A set of histograms for distances between particles
422 // ---------------------------------------------------
423 class HDelta : public Histograms
424 {
425  public:
426  HDelta (const TString & name) :
427  Histograms(name),
428  // Kinematical variables
429  // ---------------------
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 ) )
436  {}
437 
438  HDelta (TFile* outputFile, const TString & name) :
439  Histograms(outputFile, name),
440  // Kinematical variables
441  // ---------------------
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 ) )
448  {}
449 
450  HDelta (const TString & name, TFile* file) {
451  name_ = name;
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");
458  }
459 
460  ~HDelta() override {
461  delete hEta_;
462  delete hEtaSign_;
463  delete hPhi_;
464  delete hTheta_;
465  delete hCotgTheta_;
466  delete hDeltaR_;
467  }
468 
470  Fill (CLHEP::HepLorentzVector(p1.x(),p1.y(),p1.z(),p1.t()),
471  CLHEP::HepLorentzVector(p2.x(),p2.y(),p2.z(),p2.t()));
472  }
473 
474  void Fill (const CLHEP::HepLorentzVector & p1, const reco::Particle::LorentzVector & p2) override {
475  Fill (p1,CLHEP::HepLorentzVector(p2.x(),p2.y(),p2.z(),p2.t()));
476  }
477 
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());
481  hPhi_->Fill(MuScleFitUtils::deltaPhi(momentum1.phi(),momentum2.phi()));
482  hTheta_->Fill(momentum1.theta()-momentum2.theta());
483  // hCotgTheta->Fill(1/(TMath::Tan(momentum1.theta()))-1/(TMath::Tan(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()) +
488  (MuScleFitUtils::deltaPhi(momentum1.phi(),momentum2.phi()))*
489  (MuScleFitUtils::deltaPhi(momentum1.phi(),momentum2.phi()))));
490  }
491 
492  void Write() override {
493  if(histoDir_ != nullptr) histoDir_->cd();
494 
495  hEta_->Write();
496  hEtaSign_->Write();
497  hPhi_->Write();
498  hTheta_->Write();
499  hCotgTheta_->Write();
500  hDeltaR_->Write();
501  }
502 
503  void Clear() override {
504  hEta_->Clear();
505  hEtaSign_->Clear();
506  hPhi_->Clear();
507  hTheta_->Clear();
508  hDeltaR_->Clear();
509  hCotgTheta_->Clear();
510  }
511 
512  public:
513  TH1F* hEta_;
514  TH1F* hEtaSign_;
515  TH1F* hPhi_;
516  TH1F* hTheta_;
517  TH1F* hCotgTheta_;
518  TH1F* hDeltaR_;
519 };
520 
521 // ------------------------------------------------------------
522 // A set of histograms of particle kinematical variables vs eta
523 // ------------------------------------------------------------
524 class HPartVSEta : public Histograms
525 {
526  public:
527  HPartVSEta(const TString & name, const double & minMass = 0., const double & maxMass = 100., const double & maxPt = 100.)
528  {
529  name_ = name;
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",
533  32, -3.2, 3.2, 40, minMass, maxMass );
534  // TD profile histograms
535  // ---------------------
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",
539  32, -3.2, 3.2, minMass, maxMass );
540  hCurvVSEta_prof_ = new TProfile( name+"_CurvVSEta_prof", "curvature vs pseudorapidity",
541  32, -3.2, 3.2, 0, 1. );
542  }
543 
544  ~HPartVSEta() override {
545  delete hPtVSEta_;
546  delete hMassVSEta_;
547  delete hPtVSEta_prof_;
548  delete hMassVSEta_prof_;
549  delete hCurvVSEta_prof_;
550  }
551 
552  void Fill (const reco::Particle::LorentzVector & p4, const double & weight = 1.) override {
553  Fill (CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()), weight);
554  }
555 
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);
559 
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);
563  }
564 
565  void Write() override {
566  hPtVSEta_->Write();
567  hPtVSEta_prof_->Write();
568  hCurvVSEta_prof_->Write();
569  hMassVSEta_->Write();
570  hMassVSEta_prof_->Write();
571 
572  // std::vector<TGraphErrors*> graphs( (MuScleFitUtils::fitMass(hMassVSEta_)) );
573  // for (std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++) {
574  // (*graph)->Write();
575  // }
576  }
577 
578  void Clear() override {
579  hPtVSEta_->Clear();
580  hPtVSEta_prof_->Clear();
581  hCurvVSEta_prof_->Clear();
582  hMassVSEta_->Clear();
583  hMassVSEta_prof_->Clear();
584 
585  }
586 
587  public:
588 
589  TH2F *hPtVSEta_;
590  TH2F *hMassVSEta_;
591  TProfile *hMassVSEta_prof_;
592  TProfile *hPtVSEta_prof_;
593  TProfile *hCurvVSEta_prof_;
594 };
595 
596 //---------------------------------------------------------------------------
597 // A set of histograms of particle kinematical variables vs phi (in eta bins)
598 // --------------------------------------------------------------------------
599 class HPartVSPhi : public Histograms
600 {
601  public:
602  HPartVSPhi(const TString & name){
603  name_ = name;
604 // hPtVSPhi_ = new TH2F (name+"_PtVSPhi", "transverse momentum vs phi angle",
605 // 12, -3.2, 3.2, 200, 0, 200);
606 // hMassVSPhi_ = new TH2F (name+"_MassVSPhi", "mass vs phi angle",
607 // 7, -3.2, 3.2, 40, 70, 110);
608 // hMassVSPhiF_ = new TH2F (name+"_MassVSPhiF", "mass vs phi F",
609 // 7, -3.2, 3.2, 40, 70, 110);
610 // hMassVSPhiWp2_ = new TH2F (name+"_MassVSPhiWp2", "mass vs phi Wp2",
611 // 7, -3.2, 3.2, 40, 70, 110);
612 // hMassVSPhiWp1_ = new TH2F (name+"_MassVSPhiWp1", "mass vs phi Wp1",
613 // 7, -3.2, 3.2, 40, 70, 110);
614 // hMassVSPhiW0_ = new TH2F (name+"_MassVSPhiW0", "mass vs phi W0",
615 // 7, -3.2, 3.2, 40, 70, 110);
616 // hMassVSPhiWm1_ = new TH2F (name+"_MassVSPhiWm1", "mass vs phi Wm1",
617 // 7, -3.2, 3.2, 40, 70, 110);
618 // hMassVSPhiWm2_ = new TH2F (name+"_MassVSPhiWm2", "mass vs phi Wm2",
619 // 7, -3.2, 3.2, 40, 70, 110);
620 // hMassVSPhiB_ = new TH2F (name+"_MassVSPhiB", "mass vs phi B",
621 // 7, -3.2, 3.2, 40, 70, 110);
622 
623  // TD profile histograms
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);
628  }
629 
630  ~HPartVSPhi() override {
631  delete hPtVSPhi_;
632  delete hMassVSPhi_;
633  delete hMassVSPhi_prof_;
634  delete hPtVSPhi_prof_;
635 
636  // delete hMassVSPhiB_;
637 // delete hMassVSPhiWm2_;
638 // delete hMassVSPhiWm1_;
639 // delete hMassVSPhiW0_;
640 // delete hMassVSPhiWp1_;
641 // delete hMassVSPhiWp2_;
642 // delete hMassVSPhiF_;
643  }
644 
645  void Fill(const reco::Particle::LorentzVector & p4, const double & weight =1.) override {
646  Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()), weight);
647  }
648 
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);
654 
655 // if (momentum.eta()<-1.2) hMassVSPhiB_->Fill(momentum.phi(),momentum.m(), weight);
656 // else if (momentum.eta()<-0.8 && momentum.eta()>-1.2) hMassVSPhiWm2_->Fill(momentum.phi(),momentum.m(), weight);
657 // else if (momentum.eta()<-0.3 && momentum.eta()>-0.8) hMassVSPhiWm1_->Fill(momentum.phi(),momentum.m(), weight);
658 // else if ((fabs(momentum.eta())) < 0.3) hMassVSPhiW0_->Fill(momentum.phi(),momentum.m(), weight);
659 // else if (momentum.eta()>0.3 && momentum.eta()<0.8) hMassVSPhiWp1_->Fill(momentum.phi(),momentum.m(), weight);
660 // else if (momentum.eta()>0.8 && momentum.eta()<1.2) hMassVSPhiWp2_->Fill(momentum.phi(),momentum.m(), weight);
661 // else if (momentum.eta()>1.2) hMassVSPhiF_->Fill(momentum.phi(),momentum.m(), weight);
662  }
663 
664  void Write() override {
665  hPtVSPhi_->Write();
666  hMassVSPhi_->Write();
667  hMassVSPhi_prof_->Write();
668  hPtVSPhi_prof_->Write();
669 
670  // hMassVSPhiB_->Write();
671 // hMassVSPhiWm2_->Write();
672 // hMassVSPhiWm1_->Write();
673 // hMassVSPhiW0_->Write();
674 // hMassVSPhiWp1_->Write();
675 // hMassVSPhiWp2_->Write();
676 // hMassVSPhiF_->Write();
677 
678 // std::vector<TGraphErrors*> graphs ((MuScleFitUtils::fitMass(hMassVSPhi_)));
679 // for(std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++){
680 // (*graph)->Write();
681 // }
682 // std::vector<TGraphErrors*> graphsB ((MuScleFitUtils::fitMass(hMassVSPhiB_)));
683 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsB.begin(); graph != graphsB.end(); graph++){
684 // (*graph)->Write();
685 // }
686 // std::vector<TGraphErrors*> graphsWm2 ((MuScleFitUtils::fitMass(hMassVSPhiWm2_)));
687 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsWm2.begin(); graph != graphsWm2.end(); graph++){
688 // (*graph)->Write();
689 // }
690 // std::vector<TGraphErrors*> graphsWm1 ((MuScleFitUtils::fitMass(hMassVSPhiWm1_)));
691 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsWm1.begin(); graph != graphsWm1.end(); graph++){
692 // (*graph)->Write();
693 // }
694 // std::vector<TGraphErrors*> graphsW0 ((MuScleFitUtils::fitMass(hMassVSPhiW0_)));
695 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsW0.begin(); graph != graphsW0.end(); graph++){
696 // (*graph)->Write();
697 // }
698 // std::vector<TGraphErrors*> graphsWp1 ((MuScleFitUtils::fitMass(hMassVSPhiWp1_)));
699 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsWp1.begin(); graph != graphsWp1.end(); graph++){
700 // (*graph)->Write();
701 // }
702 // std::vector<TGraphErrors*> graphsWp2 ((MuScleFitUtils::fitMass(hMassVSPhiWp2_)));
703 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsWp2.begin(); graph != graphsWp2.end(); graph++){
704 // (*graph)->Write();
705 // }
706 // std::vector<TGraphErrors*> graphsF ((MuScleFitUtils::fitMass(hMassVSPhiF_)));
707 // for(std::vector<TGraphErrors*>::const_iterator graph = graphsF.begin(); graph != graphsF.end(); graph++){
708 // (*graph)->Write();
709 // }
710  }
711 
712  void Clear() override {
713  hPtVSPhi_->Clear();
714  hMassVSPhi_->Clear();
715  hPtVSPhi_prof_->Clear();
716  hMassVSPhi_prof_->Clear();
717 
718 // hMassVSPhiB_->Clear();
719 // hMassVSPhiWm2_->Clear();
720 // hMassVSPhiWm1_->Clear();
721 // hMassVSPhiW0_->Clear();
722 // hMassVSPhiWp1_->Clear();
723 // hMassVSPhiWp2_->Clear();
724 // hMassVSPhiF_->Clear();
725  }
726 
727  public:
728  TH2F *hPtVSPhi_;
729  TH2F *hMassVSPhi_;
730  TProfile *hMassVSPhi_prof_;
731  TProfile *hPtVSPhi_prof_;
732 
733 // TH2F *hMassVSPhiB_;
734 // TH2F *hMassVSPhiWm2_;
735 // TH2F *hMassVSPhiWm1_;
736 // TH2F *hMassVSPhiW0_;
737 // TH2F *hMassVSPhiWp1_;
738 // TH2F *hMassVSPhiWp2_;
739 // TH2F *hMassVSPhiF_;
740 };
741 
742 //---------------------------------------------------------------------------------------
743 // A set of histograms of particle VS pt
744 class HPartVSPt : public Histograms
745 {
746  public:
747  HPartVSPt(const TString & name){
748  name_ = name;
749  hMassVSPt_ = new TH2F( name+"_MassVSPt", "mass vs transverse momentum",
750  12, -6, 6, 40, 70, 110 );
751  // TD profile histograms
752  hMassVSPt_prof_ = new TProfile( name+"_MassVSPt_prof", "mass vs transverse momentum",
753  12, -3, 3, 86, 116 );
754  }
755 
756  ~HPartVSPt() override{
757  delete hMassVSPt_;
758  delete hMassVSPt_prof_;
759  }
760 
761  void Fill(const reco::Particle::LorentzVector & p4, const double & weight = 1.) override {
762  Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()), weight);
763  }
764 
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);
768  }
769 
770  void Write() override {
771  hMassVSPt_->Write();
772  hMassVSPt_prof_->Write();
773 
774 // std::vector<TGraphErrors*> graphs( (MuScleFitUtils::fitMass(hMassVSPt_)) );
775 // for(std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++){
776 // (*graph)->Write();
777 // }
778  }
779 
780  void Clear() override {
781  hMassVSPt_->Clear();
782  hMassVSPt_prof_->Clear();
783  }
784 
785  public:
786  TH2F *hMassVSPt_;
787  TProfile *hMassVSPt_prof_;
788 };
789 
790 // ---------------------------------------------------
791 // A set of histograms of resonance mass versus muon variables
792 class HMassVSPart : public Histograms
793 {
794  public:
795  HMassVSPart( const TString & name, const double & minMass = 0., const double & maxMass = 150., const double maxPt = 100. ) {
796  name_ = name;
797 
798  // Kinematical variables
799  // ---------------------
800  hMassVSPt_ = new TH2F( name+"_MassVSPt", "resonance mass vs muon transverse momentum", 200, 0., maxPt, 6000, minMass, maxMass );
801 
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 );
805 
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 );
808 
809 
810  // J/Psi mass -----
811 // hMassVSEtaPhiPlus_ = new TH3F( name+"_MassVSEtaPhiPlus", "resonance mass vs muon+ phi/pseudorapidity",6, -3.2, 3.2, 20, -2.5, 2.5, 6000, minMass, maxMass );
812 // hMassVSEtaPhiMinus_ = new TH3F( name+"_MassVSEtaPhiMinus", "resonance mass vs muon- phi/pseudorapidity", 6, -3.2, 3.2, 20, -2.5, 2.5, 6000, minMass, maxMass );
813 
814  //Z mass -----------
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 );
817 
818 
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 );
821 
822 
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 );
825 
826 
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.);
830 
831  //hMassVSPt_prof = new TProfile (name+"_MassVSPt_prof", "resonance mass vs muon transverse momentum", 100, 0., 200., minMass, maxMass);
832  //hMassVSEta_prof = new TProfile (name+"_MassVSEta_prof", "resonance mass vs muon pseudorapidity", 30, -6., 6., minMass, maxMass);
833  //hMassVSPhiPlus_prof = new TProfile (name+"_MassVSPhiPlus_prof", "resonance mass vs muon+ phi angle", 32, -3.2, 3.2, minMass, maxMass);
834  //hMassVSPhiMinus_prof = new TProfile (name+"_MassVSPhiMinus_prof", "resonance mass vs muon- phi angle", 32, -3.2, 3.2, minMass, maxMass);
835  }
836 
837  HMassVSPart(const TString & name, TFile* file){
838  name_=name;
839  hMassVSPt_ = (TH2F *) file->Get(name+"_MassVSPt");
840  hMassVSEta_ = (TH2F *) file->Get(name+"_MassVSEta");
841 
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");
846 
847  hMassVSPhiPlusMinusDiff_ = (TH2F *) file->Get(name+"_MassVSPhiPlusMinusDiff");
848  hMassVSEtaPlusMinusDiff_ = (TH2F *) file->Get(name+"_MassVSEtaPlusMinusDiff");
849 
850  hMassVSPhiPlus_ = (TH2F *) file->Get(name+"_MassVSPhiPlus");
851  hMassVSPhiMinus_ = (TH2F *) file->Get(name+"_MassVSPhiMinus");
852 
853  hMassVSCosThetaCS_prof = (TProfile *) file->Get(name+"_MassVScosTheta_prof");
854  //hMassVSPt_prof = (TProfile *) file->Get(name+"_MassVSPt_prof");
855  //hMassVSEta_prof = (TProfile *) file->Get(name+"_MassVSEta_prof");
856  //hMassVSPhiPlus_prof = (TProfile *) file->Get(name+"_MassVSPhiPlus_prof");
857  //hMassVSPhiMinus_prof = (TProfile *) file->Get(name+"_MassVSPhiMinus_prof");
858  }
859 
860  ~HMassVSPart() override{
861  delete hMassVSPt_;
862  delete hMassVSEta_;
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;
876  }
877 
878  void Fill(const reco::Particle::LorentzVector & p41, const reco::Particle::LorentzVector & p42, const int charge, const double & weight = 1.) override
879  {
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);
882  }
883 
885  const reco::Particle::LorentzVector & p42,
886  const reco::Particle::LorentzVector & p4Res,
887  const double & weight = 1.) override
888  {
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()),
892  weight);
893  }
894 
896  void Fill(const CLHEP::HepLorentzVector & momentum1,
897  const CLHEP::HepLorentzVector & momentum2,
898  const CLHEP::HepLorentzVector & momentumRes,
899  const double & weight = 1.) override
900  {
901 
902  /************************************************************************
903  *
904  * Observable: cos(theta) = 2 Q^-1 (Q^2+Qt^2)^-(1/2) (mu^+ mubar^- - mu^- mubar^+)
905  * (computed in Collins-Soper frame)
906  *
907  ************************************************************************/
908 
909  double costhetaCS, phiCS;
910 
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());
918  //double costheta = 2.0 / Q.Mag() / sqrt(pow(Q.Mag(), 2) + pow(Q.Pt(), 2)) * (muplus * mubarminus - muminus * mubarplus);
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;
921 
922 
923 
924 
925  /************************************************************************
926  *
927  * 3) tanphi = (Q^2 + Qt^2)^1/2 / Q (Dt dot R unit) /(Dt dot Qt unit)
928  *
929  ************************************************************************/
930 
931  // unit vector on R direction
932  CLHEP::HepLorentzVector Pbeam(0.,0.,3500.,3500.);
933  CLHEP::Hep3Vector R = Pbeam.vect().cross(Q.vect());
934  CLHEP::Hep3Vector Runit = R.unit();
935 
936 
937  // unit vector on Qt
938  CLHEP::Hep3Vector Qt = Q.vect(); Qt.setZ(0);
939  CLHEP::Hep3Vector Qtunit = Qt.unit();
940 
941 
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);
947 
948  hMassVSPhiCS_->Fill(phiCS,momentumRes.m(), weight);
949  hMassVSCosThetaCS_->Fill(costhetaCS,momentumRes.m(), weight);
950  hMassVSCosThetaCS_prof ->Fill(costhetaCS,momentumRes.m());
951  /*************************************************************************
952  *************************************************************************/
953 
954  hMassVSPhiPlusPhiMinus_->Fill(momentum1.phi(), momentum2.phi(), momentumRes.m(), weight);
955  hMassVSEtaPlusEtaMinus_->Fill(momentum1.eta(), momentum2.eta(), momentumRes.m(), weight);
956 
957  hMassVSPhiPlusMinusDiff_->Fill( (momentum1.phi()-momentum2.phi()), momentumRes.m(), weight);
958  hMassVSEtaPlusMinusDiff_->Fill( (momentum1.eta()-momentum2.eta()), momentumRes.m(), weight);
959  }
960 
961  void Fill(const CLHEP::HepLorentzVector & momentum1, const CLHEP::HepLorentzVector & momentum2, const int charge, const double & weight = 1.) override
962  {
963  hMassVSPt_->Fill(momentum1.perp(),momentum2.m(), weight);
964  //hMassVSPt_prof_->Fill(momentum1.perp(),momentum2.m());
965 
966 
967  hMassVSEta_->Fill(momentum1.eta(),momentum2.m(), weight);
968  //hMassVSEta_prof_->Fill(momentum1.eta(),momentum2.m());
969 
970  if(charge>0){
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);
974  }
975  else if(charge<0){
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);
979  }
980  else {
981  LogDebug("Histograms") << "HMassVSPart: wrong charge value = " << charge << std::endl;
982  abort();
983  }
984  }
985 
986  void Write() override {
987  hMassVSPt_->Write();
988  hMassVSEta_->Write();
989  hMassVSPhiPlus_->Write();
990  hMassVSPhiMinus_->Write();
991 
992  hMassVSEtaPhiPlus_->Write();
993  hMassVSEtaPhiMinus_->Write();
994  hMassVSEtaPlus_->Write();
995  hMassVSEtaMinus_->Write();
996 
997  hMassVSPhiPlusPhiMinus_->Write();
998  hMassVSEtaPlusEtaMinus_->Write();
999  hMassVSCosThetaCS_->Write();
1000  hMassVSPhiCS_->Write();
1001 
1002  hMassVSPhiPlusMinusDiff_->Write();
1003  hMassVSEtaPlusMinusDiff_->Write();
1004  hMassVSCosThetaCS_prof->Write();
1005 
1006  //hMassVSPt_prof_->Write();
1007  //hMassVSEta_prof_->Write();
1008  //hMassVSPhiPlus_prof_->Write();
1009  //hMassVSPhiMinus_prof_->Write();
1010  }
1011 
1012  void Clear() override {
1013  hMassVSPt_->Clear();
1014  hMassVSEta_->Clear();
1015  hMassVSPhiPlus_->Clear();
1016  hMassVSPhiMinus_->Clear();
1017 
1018  hMassVSEtaPhiPlus_->Clear();
1019  hMassVSEtaPhiMinus_->Clear();
1020  hMassVSEtaPlus_->Clear();
1021  hMassVSEtaMinus_->Clear();
1022 
1023  hMassVSPhiPlusPhiMinus_->Clear();
1024  hMassVSEtaPlusEtaMinus_->Clear();
1025  hMassVSCosThetaCS_->Clear();
1026  hMassVSPhiCS_->Clear();
1027  hMassVSPhiPlusMinusDiff_->Clear();
1028  hMassVSEtaPlusMinusDiff_->Clear();
1029  hMassVSCosThetaCS_prof->Clear();
1030 
1031  //hMassVSPt_prof_->Clear();
1032  //hMassVSEta_prof_->Clear();
1033  //hMassVSPhiPlus_prof_->Clear();
1034  //hMassVSPhiMinus_prof_->Clear();
1035  }
1036 
1037  protected:
1038  TH2F* hMassVSPt_;
1044 
1049 
1052 
1055 
1057 
1058  //TProfile* hMassVSPt_prof_;
1059  //TProfile* hMassVSEta_prof_;
1060  //TProfile* hMassVSPhiPlus_prof_;
1061  //TProfile* hMassVSPhiMinus_prof_;
1062 };
1063 
1064 
1065 // ---------------------------------------------------
1066 // A set of histograms of Z mass versus muon variables
1067 class HMassVSPartProfile : public Histograms
1068 {
1069  public:
1070  HMassVSPartProfile( const TString & name, const double & minMass = 0., const double & maxMass = 150., const double maxPt = 100. ) {
1071  name_ = name;
1072 
1073  // Kinematical variables
1074  // ---------------------
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. );
1079  }
1080 
1081  HMassVSPartProfile(const TString & name, TFile* file){
1082  name_=name;
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");
1087  }
1088 
1090  delete hMassVSPt_;
1091  delete hMassVSEta_;
1092  delete hMassVSPhiPlus_;
1093  delete hMassVSPhiMinus_;
1094  }
1095 
1096  void Fill(const reco::Particle::LorentzVector & p41, const reco::Particle::LorentzVector & p42, const int charge, const double & weight = 1.) override {
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);
1099  }
1100 
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);
1104  if(charge>0){
1105  hMassVSPhiPlus_->Fill(momentum1.phi(), momentum2.m(), weight);
1106  }
1107  else if(charge<0){
1108  hMassVSPhiMinus_->Fill(momentum1.phi(), momentum2.m(), weight);
1109  }
1110  else {
1111  LogDebug("Histograms") << "HMassVSPartProfile: wrong charge value = " << charge << std::endl;
1112  abort();
1113  }
1114  }
1115 
1116  void Write() override {
1117  hMassVSPt_->Write();
1118  hMassVSEta_->Write();
1119  hMassVSPhiPlus_->Write();
1120  hMassVSPhiMinus_->Write();
1121  }
1122 
1123  void Clear() override {
1124  hMassVSPt_->Clear();
1125  hMassVSEta_->Clear();
1126  hMassVSPhiPlus_->Clear();
1127  hMassVSPhiMinus_->Clear();
1128  }
1129 
1130  protected:
1131  TProfile2D* hMassVSPt_;
1132  TProfile2D* hMassVSEta_;
1133  TProfile2D* hMassVSPhiPlus_;
1134  TProfile2D* hMassVSPhiMinus_;
1135 };
1136 
1137 //---------------------------------------------------------------------------------------
1139 class HResolutionVSPart : public Histograms
1140 {
1141  public:
1142  HResolutionVSPart(TFile * outputFile, const TString & name, const double maxPt = 100,
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)
1146  {
1147  // Kinematical variables
1148 
1149  // hReso = new TH1F (name+"_Reso", "resolution", 4000, -1, 1);
1150  // hResoVSPtEta = new TH2F (name+"_ResoVSPtEta", "resolution VS pt and #eta", 200, 0, 200, 60, -3, 3);
1151  // hResoVSPt = new TH2F (name+"_ResoVSPt", "resolution VS pt", 200, 0, 200, 8000, -1, 1);
1152  // //hResoVSPt_prof = new TProfile (name+"_ResoVSPt_prof", "resolution VS pt", 100, 0, 200, -1, 1);
1153  // hResoVSEta = new TH2F (name+"_ResoVSEta", "resolution VS eta", 60, -3, 3, 8000, yMinEta, yMaxEta);
1154  // hResoVSTheta = new TH2F (name+"_ResoVSTheta", "resolution VS theta", 30, 0, TMath::Pi(), 8000, -1, 1);
1155  // //hResoVSEta_prof = new TProfile (name+"_ResoVSEta_prof", "resolution VS eta", 10, -2.5, 2.5, -1, 1);
1156  // hResoVSPhiPlus = new TH2F (name+"_ResoVSPhiPlus", "resolution VS phi mu+", 14, -3.2, 3.2, 8000, -1, 1);
1157  // hResoVSPhiMinus = new TH2F (name+"_ResoVSPhiMinus", "resolution VS phi mu-", 14, -3.2, 3.2, 8000, -1, 1);
1158  // //hResoVSPhi_prof = new TProfile (name+"_ResoVSPhi_prof", "resolution VS phi", 14, -3.2, 3.2, -1, 1);
1159  // hAbsReso = new TH1F (name+"_AbsReso", "resolution", 100, 0, 1);
1160  // hAbsResoVSPt = new TH2F (name+"_AbsResoVSPt", "Abs resolution VS pt", 200, 0, 500, 100, 0, 1);
1161  // hAbsResoVSEta = new TH2F (name+"_AbsResoVSEta", "Abs resolution VS eta", 60, -3, 3, 100, 0, 1);
1162  // hAbsResoVSPhi = new TH2F (name+"_AbsResoVSPhi", "Abs resolution VS phi", 14, -3.2, 3.2, 100, 0, 1);
1163 
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 );
1180  if( doProfiles_ ) {
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 );
1189  }
1190  }
1191 
1192  HResolutionVSPart(const TString & name, TFile* file) {
1193  name_=name;
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");
1210  if( doProfiles_ ) {
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");
1219  }
1220  }
1221 
1222  ~HResolutionVSPart() override {
1223  delete hReso_;
1224  delete hResoVSPtEta_;
1225  delete hResoVSPt_;
1226  delete hResoVSPt_Bar_;
1227  delete hResoVSPt_Endc_17_;
1228  delete hResoVSPt_Endc_20_;
1229  delete hResoVSPt_Endc_24_;
1230  delete hResoVSPt_Ovlap_;
1231  delete hResoVSEta_;
1232  delete hResoVSTheta_;
1233  delete hResoVSPhiMinus_;
1234  delete hResoVSPhiPlus_;
1235  delete hAbsReso_;
1236  delete hAbsResoVSPt_;
1237  delete hAbsResoVSEta_;
1238  delete hAbsResoVSPhi_;
1239  if( doProfiles_ ) {
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_;
1248  }
1249  }
1250 
1251  void Fill(const reco::Particle::LorentzVector & p4, const double & resValue, const int charge) override {
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);
1257  if(fabs(eta)<=0.9)
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);
1265  else
1266  hResoVSPt_Endc_24_->Fill(pt,resValue);
1267 
1268  hResoVSEta_->Fill(eta,resValue);
1269  hResoVSTheta_->Fill(p4.Theta(),resValue);
1270  if(charge>0)
1271  hResoVSPhiPlus_->Fill(p4.Phi(),resValue);
1272  else if(charge<0)
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));
1278  if( doProfiles_ ) {
1279  hResoVSPt_prof_->Fill(p4.Pt(),resValue);
1280  if(fabs(eta)<=0.9)
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);
1288  else
1289  hResoVSPt_Endc_24_prof_->Fill(pt,resValue);
1290 
1291  hResoVSEta_prof_->Fill(p4.Eta(),resValue);
1292  hResoVSPhi_prof_->Fill(p4.Phi(),resValue);
1293  }
1294  }
1295 
1296  void Write() override {
1297  if(histoDir_ != nullptr) histoDir_->cd();
1298 
1299  hReso_->Write();
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();
1311  hAbsReso_->Write();
1312  hAbsResoVSPt_->Write();
1313  hAbsResoVSEta_->Write();
1314  hAbsResoVSPhi_->Write();
1315  if( doProfiles_ ) {
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();
1324  }
1325  }
1326 
1327  void Clear() override {
1328  hReso_->Clear();
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();
1340  hAbsReso_->Clear();
1341  hAbsResoVSPt_->Clear();
1342  hAbsResoVSEta_->Clear();
1343  hAbsResoVSPhi_->Clear();
1344  if( doProfiles_ ) {
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();
1353  }
1354  }
1355 
1356  public:
1357  TH1F* hReso_;
1359  TH2F* hResoVSPt_;
1365  TProfile* hResoVSPt_prof_;
1373  TProfile* hResoVSEta_prof_;
1376  TProfile* hResoVSPhi_prof_;
1377  TH1F* hAbsReso_;
1382 };
1383 
1384 // -------------------------------------------------------------
1385 // A set of histograms of likelihood value versus muon variables
1386 // -------------------------------------------------------------
1387 class HLikelihoodVSPart : public Histograms
1388 {
1389  public:
1390  HLikelihoodVSPart(const TString & name){
1391  name_ = name;
1392 
1393  // Kinematical variables
1394  // ---------------------
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.);
1401  }
1402 
1403  HLikelihoodVSPart(const TString & name, TFile* file){
1404  name_ = name;
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");
1411  }
1412 
1414  delete hLikeVSPt_;
1415  delete hLikeVSEta_;
1416  delete hLikeVSPhi_;
1417  delete hLikeVSPt_prof_;
1418  delete hLikeVSEta_prof_;
1419  delete hLikeVSPhi_prof_;
1420  }
1421 
1422  void Fill(const reco::Particle::LorentzVector & p4, const double & likeValue) override {
1423  Fill(CLHEP::HepLorentzVector(p4.x(),p4.y(),p4.z(),p4.t()), likeValue);
1424  }
1425 
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);
1433  }
1434 
1435  void Write() override {
1436  hLikeVSPt_->Write();
1437  hLikeVSEta_->Write();
1438  hLikeVSPhi_->Write();
1439  hLikeVSPt_prof_->Write();
1440  hLikeVSEta_prof_->Write();
1441  hLikeVSPhi_prof_->Write();
1442  }
1443 
1444  void Clear() override {
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");
1451  }
1452 
1453  public:
1454  TH2F* hLikeVSPt_;
1457  TProfile* hLikeVSPt_prof_;
1458  TProfile* hLikeVSEta_prof_;
1459  TProfile* hLikeVSPhi_prof_;
1460 };
1461 
1471 class HFunctionResolution : public Histograms
1472 {
1473  public:
1474  HFunctionResolution(TFile * outputFile, const TString & name, const double & ptMax = 100, const int totBinsY = 300) : Histograms(outputFile, name) {
1475  name_ = name;
1476  totBinsX_ = 300;
1477  totBinsY_ = totBinsY;
1478  xMin_ = 0.;
1479  yMin_ = -3.0;
1480  double xMax = ptMax;
1481  double yMax = 3.0;
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 );
1486  // Create and initialize the resolution arrays
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;
1495  }
1496  }
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);
1504  //hResoVSTheta_prof_ = new TProfile( name+"_ResoVSTheta_prof", "resolution VS theta", 30, 0, TMath::Pi(), 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);
1508  }
1510  delete hReso_;
1511  delete hResoVSPtEta_;
1512  // Free the resolution arrays
1513  for( int i=0; i<totBinsX_; ++i ) {
1514  delete[] resoVsPtEta_[i];
1515  delete[] resoCount_[i];
1516  }
1517  delete[] resoVsPtEta_;
1518  delete[] resoCount_;
1519  // Free the profile histograms
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_;
1527  //delete hResoVSTheta_prof_;
1528  delete hResoVSPhiPlus_prof_;
1529  delete hResoVSPhiMinus_prof_;
1530  delete hResoVSPhi_prof_;
1531  }
1532  void Fill(const reco::Particle::LorentzVector & p4, const double & resValue, const int charge) override {
1533  if( resValue != resValue ) return;
1534  hReso_->Fill(resValue);
1535 
1536  // Fill the arrays with the resolution value and count
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;
1541  // ATTENTION: we count only for positive muons because we are considering di-muon resonances
1542  // and we use this counter to compute the mean in the end. The resoVsPtEta value is filled with each muon,
1543  // but they must be considered independently (analogous to a factor 2) so in the end we would have
1544  // to divide by N/2, that is why we only increase the counter for half the muons.
1545  // if( charge > 0 )
1546  // No more. Changing it here influences also other uses of this class. The macro FunctionTerms.cc
1547  // multiplies the terms by the 2 factor.
1548 
1549  resoCount_[xIndex][yIndex] += 1;
1550 
1551  // hResoVSPtEta->Fill(p4.Pt(), p4.Eta(), resValue);
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);
1561  else
1562  hResoVSPt_Endc_24_prof_->Fill(p4.Pt(),resValue);
1563  hResoVSEta_prof_->Fill(p4.Eta(),resValue);
1564  //hResoVSTheta_prof_->Fill(p4.Theta(),resValue);
1565  if(charge>0)
1566  hResoVSPhiPlus_prof_->Fill(p4.Phi(),resValue);
1567  else if(charge<0)
1568  hResoVSPhiMinus_prof_->Fill(p4.Phi(),resValue);
1569  hResoVSPhi_prof_->Fill(p4.Phi(),resValue);
1570  }
1571  }
1572 
1573  void Write() override {
1574  if(histoDir_ != nullptr) histoDir_->cd();
1575 
1576  hReso_->Write();
1577 
1578  for( int i=0; i<totBinsX_; ++i ) {
1579  for( int j=0; j<totBinsY_; ++j ) {
1580  int N = resoCount_[i][j];
1581  // Fill with the mean value
1582  if( N != 0 ) hResoVSPtEta_->SetBinContent( i+1, j+1, resoVsPtEta_[i][j]/N );
1583  else hResoVSPtEta_->SetBinContent( i+1, j+1, 0 );
1584  }
1585  }
1586  hResoVSPtEta_->Write();
1587 
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();
1595  //hResoVSTheta_prof_->Write();
1596  hResoVSPhiMinus_prof_->Write();
1597  hResoVSPhiPlus_prof_->Write();
1598  hResoVSPhi_prof_->Write();
1599 
1600  TCanvas canvas(TString(hResoVSPtEta_->GetName())+"_canvas", TString(hResoVSPtEta_->GetTitle())+" canvas", 1000, 800);
1601  canvas.Divide(2);
1602  canvas.cd(1);
1603  hResoVSPtEta_->Draw("lego");
1604  canvas.cd(2);
1605  hResoVSPtEta_->Draw("surf5");
1606  canvas.Write();
1607  hResoVSPtEta_->Write();
1608 
1609  outputFile_->cd();
1610  }
1611 
1612  void Clear() override {
1613  hReso_->Clear();
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();
1622  //hResoVSTheta_prof_->Clear();
1623  hResoVSPhiPlus_prof_->Clear();
1624  hResoVSPhiMinus_prof_->Clear();
1625  hResoVSPhi_prof_->Clear();
1626  }
1627 
1628  protected:
1629  int getXindex(const double & x) const {
1630  return int((x-xMin_)/deltaX_*totBinsX_);
1631  }
1632  int getYindex(const double & y) const {
1633  return int((y-yMin_)/deltaY_*totBinsY_);
1634  }
1635  TH1F* hReso_;
1637  double ** resoVsPtEta_;
1638  int ** resoCount_;
1639  TProfile* hResoVSPt_prof_;
1645  TProfile* hResoVSEta_prof_;
1646  //TProfile* hResoVSTheta_prof_;
1649  TProfile* hResoVSPhi_prof_;
1650  int totBinsX_, totBinsY_;
1651  double xMin_, yMin_;
1652  double deltaX_, deltaY_;
1653 };
1654 
1656 {
1657 public:
1658  HFunctionResolutionVarianceCheck(TFile * outputFile, const TString & name, const double ptMax = 200) :
1659  HFunctionResolution(outputFile, name, ptMax)
1660  {
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;
1667  namei << i;
1668  namej << j;
1669  histoVarianceCheck_[i][j] = new TH1D(name+"_"+namei.str()+"_"+namej.str(), name, 100, 0., 1.);
1670  }
1671  }
1672  }
1674  for( int i=0; i<totBinsX_; ++i ) {
1675  for( int j=0; j<totBinsY_; ++j ) {
1676  delete histoVarianceCheck_[i][j];
1677  }
1678  delete[] histoVarianceCheck_;
1679  }
1680  delete[] histoVarianceCheck_;
1681  }
1682  void Fill(const reco::Particle::LorentzVector & p4, const double & resValue, const int charge) override {
1683  if( resValue != resValue ) return;
1684  // Need to convert the (x,y) values to the array indeces
1685  int xIndex = getXindex(p4.Pt());
1686  int yIndex = getYindex(p4.Eta());
1687  // Only fill values if they are in the selected range
1688  if ( 0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_ ) {
1689  histoVarianceCheck_[xIndex][yIndex]->Fill(resValue);
1690  }
1691  // Call also the fill of the base class
1692  HFunctionResolution::Fill( p4, resValue, charge );
1693  }
1694  void Write() override {
1695  if(histoDir_ != nullptr) histoDir_->cd();
1696  for( int xBin=0; xBin<totBinsX_; ++xBin ) {
1697  for( int yBin=0; yBin<totBinsY_; ++yBin ) {
1698  histoVarianceCheck_[xBin][yBin]->Write();
1699  }
1700  }
1702  }
1703 protected:
1705 };
1706 
1715 class HResolution : public TH1F {
1716 public:
1717  HResolution( const TString & name, const TString & title,
1718  const int totBins, const double & xMin, const double & xMax,
1719  const double & yMin, const double & yMax, TDirectory * dir = nullptr) :
1720  dir_(dir),
1721  dir2D_(nullptr),
1722  diffDir_(nullptr)
1723  {
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");
1729  }
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);
1733  }
1734  ~HResolution() override {
1735  delete diffHisto_;
1736  delete histo2D_;
1737  delete resoHisto_;
1738  }
1739  Int_t Fill( Double_t x, Double_t y ) override {
1740  diffHisto_->Fill(x, y);
1741  histo2D_->Fill(x, y);
1742  return 0;
1743  }
1744  Int_t Write(const char* name = nullptr, Int_t option = 0, Int_t bufsize = 0) override {
1745  // Loop on all the bins and take the rms.
1746  // The TProfile bin error is by default the standard error on the mean, that is
1747  // rms/sqrt(N). If it is created with the "S" option (as we did NOT do), it would
1748  // already be the rms. Thus we take the error and multiply it by the sqrt of the
1749  // bin entries to get the rms.
1750  // bin 0 is the underflow, bin totBins+1 is the overflow.
1751  unsigned int totBins = diffHisto_->GetNbinsX();
1752  // std::cout << "totBins = " << totBins << std::endl;
1753  for( unsigned int iBin=1; iBin<=totBins; ++iBin ) {
1754  // std::cout << "iBin = " << iBin << ", " << diffHisto_->GetBinError(iBin)*sqrt(diffHisto_->GetBinEntries(iBin)) << std::endl;
1755  resoHisto_->SetBinContent( iBin, diffHisto_->GetBinError(iBin)*sqrt(diffHisto_->GetBinEntries(iBin)) );
1756  }
1757  if( dir_ != nullptr ) dir_->cd();
1758  resoHisto_->Write();
1759  if( diffDir_ != nullptr ) diffDir_->cd();
1760  diffHisto_->Write();
1761  if( dir2D_ != nullptr ) dir2D_->cd();
1762  histo2D_->Write();
1763 
1764  return 0;
1765  }
1766  protected:
1767  TDirectory * dir_;
1768  TDirectory * dir2D_;
1769  TDirectory * diffDir_;
1770  TProfile * diffHisto_;
1771  TH2F * histo2D_;
1772  TH1F * resoHisto_;
1773 };
1774 
1783 {
1784  public:
1786  productXY_(0),
1787  sumX_(0),
1788  sumY_(0),
1789  N_(0)
1790  {}
1791  void fill(const double & x, const double & y) {
1792  productXY_ += x*y;
1793  sumX_ += x;
1794  sumY_ += y;
1795  ++N_;
1796  }
1797  double covariance() {
1798  if( N_ != 0 ) {
1799  double meanX = sumX_/N_;
1800  double meanY = sumY_/N_;
1801  // std::cout << "meanX*meanY = "<<meanX<<"*"<<meanY<< " = " << meanX*meanY << std::endl;
1802  return (productXY_/N_ - meanX*meanY);
1803  }
1804  return 0.;
1805  }
1806  double getN() {return N_;}
1807  protected:
1808  double productXY_;
1809  double sumX_;
1810  double sumY_;
1811  int N_;
1812 };
1813 
1818 class HCovarianceVSxy : public Histograms
1819 {
1820  public:
1821  HCovarianceVSxy( const TString & name, const TString & title,
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),
1827  readMode_(false),
1828  varianceCheck_(varianceCheck)
1829  {
1830  name_ = name;
1831  histoDir_ = dir;
1832  histoCovariance_ = new TH2D(name+"Covariance", title+" covariance", totBinsX, xMin, xMax, totBinsY, yMin, yMax);
1833 
1834  covariances_ = new Covariance*[totBinsX];
1835  for( int i=0; i<totBinsX; ++i ) {
1836  covariances_[i] = new Covariance[totBinsY];
1837  }
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;
1845  namei << i;
1846  namej << j;
1847  histoVarianceCheck_[i][j] = new TH1D(name+"_"+namei.str()+"_"+namej.str(), name, 10000, -1, 1);
1848  }
1849  }
1850  }
1851  }
1853  HCovarianceVSxy( TFile * inputFile, const TString & name, const TString & dirName ) :
1854  readMode_(true)
1855  {
1856  histoDir_ = (TDirectory*)(inputFile->Get(dirName.Data()));
1857  if( histoDir_ == nullptr ) {
1858  std::cout << "Error: directory not found" << std::endl;
1859  exit(0);
1860  }
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_;
1868  }
1869 
1870  ~HCovarianceVSxy() override {
1871  delete histoCovariance_;
1872  // Free covariances
1873  for(int i=0; i<totBinsX_; ++i) {
1874  delete[] covariances_[i];
1875  }
1876  delete[] covariances_;
1877  // Free variance check histograms
1878  if( varianceCheck_ ) {
1879  for( int i=0; i<totBinsX_; ++i ) {
1880  for( int j=0; j<totBinsY_; ++j ) {
1881  delete histoVarianceCheck_[i][j];
1882  }
1883  delete[] histoVarianceCheck_[i];
1884  }
1885  delete[] histoVarianceCheck_;
1886  }
1887  }
1888 
1893  void Fill( const double & x, const double & y, const double & a, const double & b ) override {
1894  // Need to convert the (x,y) values to the array indeces
1895  int xIndex = getXindex(x);
1896  int yIndex = getYindex(y);
1897  // Only fill values if they are in the selected range
1898  if ( 0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_ ) {
1899  // if( TString(histoCovariance_->GetName()).Contains("CovarianceCotgTheta_Covariance") )
1900  covariances_[xIndex][yIndex].fill(a,b);
1901  // Should be used with the variance, in which case a==b
1902  if( varianceCheck_ ) histoVarianceCheck_[xIndex][yIndex]->Fill(a);
1903  }
1904  }
1905 
1906  using Histograms::Get;
1907  double Get( const double & x, const double & y ) const {
1908  // Need to convert the (x,y) values to the array indeces
1909  int xIndex = getXindex(x);
1910  int yIndex = getYindex(y);
1911  // If the values exceed the limits of the histogram, return the border values
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);
1917  }
1918 
1919  void Write() override {
1920  if( !readMode_ ) {
1921  std::cout << "writing: " << histoCovariance_->GetName() << std::endl;
1922  for( int xBin=0; xBin<totBinsX_; ++xBin ) {
1923  for( int yBin=0; yBin<totBinsY_; ++yBin ) {
1924  double covariance = covariances_[xBin][yBin].covariance();
1925  // Histogram bins start from 1
1926  // std::cout << "covariance["<<xBin<<"]["<<yBin<<"] with N = "<<covariances_[xBin][yBin].getN()<<" is: " << covariance << std::endl;
1927  histoCovariance_->SetBinContent(xBin+1, yBin+1, covariance);
1928  }
1929  }
1930  if( histoDir_ != nullptr ) histoDir_->cd();
1931  TCanvas canvas(TString(histoCovariance_->GetName())+"_canvas", TString(histoCovariance_->GetTitle())+" canvas", 1000, 800);
1932  canvas.Divide(2);
1933  canvas.cd(1);
1934  histoCovariance_->Draw("lego");
1935  canvas.cd(2);
1936  histoCovariance_->Draw("surf5");
1937  canvas.Write();
1938  histoCovariance_->Write();
1939 
1940  TDirectory * binsDir = nullptr;
1941  if( varianceCheck_ ) {
1942  if ( histoDir_ != nullptr ) {
1943  histoDir_->cd();
1944  if( binsDir == nullptr ) binsDir = histoDir_->mkdir(name_+"Bins");
1945  binsDir->cd();
1946  }
1947  for( int xBin=0; xBin<totBinsX_; ++xBin ) {
1948  for( int yBin=0; yBin<totBinsY_; ++yBin ) {
1949  histoVarianceCheck_[xBin][yBin]->Write();
1950  }
1951  }
1952  }
1953  }
1954  }
1955  void Clear() override {
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();
1961  }
1962  }
1963  }
1964  }
1965  protected:
1966  int getXindex(const double & x) const {
1967  return int((x-xMin_)/deltaX_*totBinsX_);
1968  }
1969  int getYindex(const double & y) const {
1970  return int((y-yMin_)/deltaY_*totBinsY_);
1971  }
1974  int totBinsX_, totBinsY_, totBinsZ_;
1975  double xMin_, deltaX_, yMin_, deltaY_;
1979 };
1980 
1985 class HCovarianceVSParts : public Histograms
1986 {
1987  public:
1988  HCovarianceVSParts(TFile * outputFile, const TString & name, const double & ptMax ) : Histograms( outputFile, name ) {
1989  int totBinsX = 40;
1990  int totBinsY = 40;
1991  double etaMin = -3.;
1992  double etaMax = 3.;
1993  double ptMin = 0.;
1994 
1995  readMode_ = false;
1996 
1997  // Variances
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);
2001  // Covariances
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_);
2011  }
2012 
2014  HCovarianceVSParts( const TString & inputFileName, const TString & name )
2015  {
2016  name_ = name;
2017  TFile * inputFile = new TFile(inputFileName, "READ");
2018  readMode_ = true;
2019 
2020  // Variances
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_);
2024  // Covariances
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_);
2034  }
2035 
2036  ~HCovarianceVSParts() override {
2037  for (std::map<TString, HCovarianceVSxy*>::const_iterator histo=mapHisto_.begin();
2038  histo!=mapHisto_.end(); histo++) {
2039  delete (*histo).second;
2040  }
2041  }
2042 
2043  double Get( const reco::Particle::LorentzVector & recoP1, const TString & covarianceName ) override {
2044  return mapHisto_[name_+covarianceName]->Get(recoP1.pt(), recoP1.eta());
2045  }
2046 
2047  void Fill( const reco::Particle::LorentzVector & recoP1,
2048  const reco::Particle::LorentzVector & genP1,
2049  const reco::Particle::LorentzVector & recoP2,
2050  const reco::Particle::LorentzVector & genP2 ) override {
2051 
2052  double pt1 = recoP1.pt();
2053  double eta1 = recoP1.eta();
2054  double pt2 = recoP2.pt();
2055  double eta2 = recoP2.eta();
2056 
2057  double diffPt1 = (pt1 - genP1.pt())/genP1.pt();
2058  double diffPt2 = (pt2 - genP2.pt())/genP2.pt();
2059 
2060  double genTheta1 = genP1.theta();
2061  double genTheta2 = genP2.theta();
2062  double recoTheta1 = recoP1.theta();
2063  double recoTheta2 = recoP2.theta();
2064 
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));
2069 
2070  // double diffCotgTheta1 = (recoCotgTheta1 - genCotgTheta1)/genCotgTheta1;
2071  // double diffCotgTheta2 = (recoCotgTheta2 - genCotgTheta2)/genCotgTheta2;
2072  double diffCotgTheta1 = recoCotgTheta1 - genCotgTheta1;
2073  double diffCotgTheta2 = recoCotgTheta2 - genCotgTheta2;
2074 
2075  // double diffPhi1 = (recoP1.phi() - genP1.phi())/genP1.phi();
2076  // double diffPhi2 = (recoP2.phi() - genP2.phi())/genP2.phi();
2077  double diffPhi1 = MuScleFitUtils::deltaPhiNoFabs(recoP1.phi(), genP1.phi());
2078  double diffPhi2 = MuScleFitUtils::deltaPhiNoFabs(recoP2.phi(), genP2.phi());
2079 
2080  // Fill the variances
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);
2087 
2088  // Fill these histograms with both muons
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);
2095 
2096  // We fill two (pt, eta) bins for each pair of values. The bin of the
2097  // first and of the second muon. This should take account for the
2098  // assumed symmetry between the exchange of the first with the second muon.
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);
2105 
2106  // Fill the following histograms again for each muon (pt, eta) bin. Same
2107  // reason as in the previous case. If the symmetry is true, it does not
2108  // make any difference the order by which we fill the pt and cotgTheta combinations.
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);
2115  }
2116  void Write() override {
2117  if( !readMode_ ) {
2118  histoDir_->cd();
2119  for (std::map<TString, HCovarianceVSxy*>::const_iterator histo=mapHisto_.begin();
2120  histo!=mapHisto_.end(); histo++) {
2121  (*histo).second->Write();
2122  }
2123  }
2124  }
2125  void Clear() override {
2126  for (std::map<TString, HCovarianceVSxy*>::const_iterator histo=mapHisto_.begin();
2127  histo!=mapHisto_.end(); histo++) {
2128  (*histo).second->Clear();
2129  }
2130  }
2131  protected:
2132  std::map<TString, HCovarianceVSxy*> mapHisto_;
2134 };
2135 
2145 class HMassResolutionVSPart : public Histograms
2146 {
2147  public:
2148  HMassResolutionVSPart(TFile * outputFile, const TString & name) : Histograms( outputFile, name ) {
2149  // Kinematical variables
2150  nameSuffix_[0] = "Plus";
2151  nameSuffix_[1] = "Minus";
2152  TString titleSuffix[] = {" for mu+", " for mu-"};
2153 
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_);
2158 
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_);
2163  }
2164 
2165  // single particles histograms
2166  muMinus.reset( new HDelta("muMinus") );
2167  muPlus.reset( new HDelta("muPlus") );
2168  }
2169 
2171  for (std::map<TString, TH1*>::const_iterator histo=mapHisto_.begin();
2172  histo!=mapHisto_.end(); histo++) {
2173  delete (*histo).second;
2174  }
2175  }
2176 
2177  void Fill( const reco::Particle::LorentzVector & recoP1, const int charge1,
2178  const reco::Particle::LorentzVector & genP1,
2179  const reco::Particle::LorentzVector & recoP2, const int charge2,
2180  const reco::Particle::LorentzVector & genP2,
2181  const double & recoMass, const double & genMass ) override {
2182  muMinus->Fill(recoP1, genP1);
2183  muPlus->Fill(recoP2, genP2);
2184 
2185  Fill( recoP1, charge1, recoP2, charge2, recoMass, genMass );
2186  }
2187 
2188  void Fill( const reco::Particle::LorentzVector & recoP1, const int charge1,
2189  // const reco::Particle::LorentzVector & genP1,
2190  const reco::Particle::LorentzVector & recoP2, const int charge2,
2191  // const reco::Particle::LorentzVector & genP2,
2192  const double & recoMass, const double & genMass ) override {
2193 
2194  if ( charge1 == charge2 ) std::cout << "Error: must get two opposite charge particles" << std::endl;
2195 
2196  double massRes = (recoMass - genMass)/genMass;
2197 
2198  reco::Particle::LorentzVector recoPair( recoP1 + recoP2 );
2199  double pairPt = recoPair.Pt();
2200 
2201  double recoPt[2] = {recoP1.Pt(), recoP2.Pt()};
2202  double recoEta[2] = {recoP1.Eta(), recoP2.Eta()};
2203  double recoPhi[2] = {recoP1.Phi(), recoP2.Phi()};
2204 
2205  // std::cout << "pairPt = " << pairPt << ", massRes = ("<<recoMass<<" - "<<genMass<<")/"<<genMass<<" = " << massRes
2206  // << ", recoPt[0] = " << recoPt[0] << ", recoPt[1] = " << recoPt[1] << std::endl;
2207 
2208  // Index of the histogram. If the muons have charge1 = -1 and charge2 = 1, they already have the
2209  // correct histogram indeces. Otherwise, swap the indeces.
2210  // In any case the negative muon has index = 0 and the positive muon has index = 1.
2211  int id[2] = {0,1};
2212  if ( charge1 > 0 ) {
2213  id[0] = 1;
2214  id[1] = 0;
2215  }
2216 
2217  double pairDeltaEta = fabs(recoEta[0] - recoEta[1]);
2218  double pairDeltaPhi = MuScleFitUtils::deltaPhi(recoPhi[0], recoPhi[1]);
2219 
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);
2224 
2225  // Resolution vs single muon quantities
2226  // ------------------------------------
2227  int index = 0;
2228  for( int i=0; i<2; ++i ) {
2229  index = id[i];
2230  mapHisto_[name_+"VSPt"+nameSuffix_[i]]->Fill(recoPt[index], massRes); // EM [index] or [i]???
2231  mapHisto_[name_+"VSEta"+nameSuffix_[i]]->Fill(recoEta[index], massRes);
2232  mapHisto_[name_+"VSPhi"+nameSuffix_[i]]->Fill(recoPhi[index], massRes);
2233  }
2234  }
2235 
2236  void Write() override {
2237  histoDir_->cd();
2238  for (std::map<TString, TH1*>::const_iterator histo=mapHisto_.begin();
2239  histo!=mapHisto_.end(); histo++) {
2240  (*histo).second->Write();
2241  }
2242  // Create the new dir and cd into it
2243  (histoDir_->mkdir("singleMuonsVSgen"))->cd();
2244  muMinus->Write();
2245  muPlus->Write();
2246  }
2247 
2248  void Clear() override {
2249  for (std::map<TString, TH1*>::const_iterator histo=mapHisto_.begin();
2250  histo!=mapHisto_.end(); histo++) {
2251  (*histo).second->Clear();
2252  }
2253  muMinus->Clear();
2254  muPlus->Clear();
2255  }
2256 
2257 // HMassResolutionVSPart(const TString & name, TFile* file){
2258 // TString nameSuffix[] = {"Plus", "Minus"};
2259 // name_ = name;
2260 // hReso = (TH1F *) file->Get(name+"_Reso");
2261 // hResoVSPairPt = (TH2F *) file->Get(name+"_ResoVSPairPt");
2262 // hResoVSPairDeltaEta = (TH2F *) file->Get(name+"_ResoVSPairDeltaEta");
2263 // hResoVSPairDeltaPhi = (TH2F *) file->Get(name+"_ResoVSPairDeltaPhi");
2264 // for( int i=0; i<2; ++i ) {
2265 // hResoVSPt[i] = (TH2F *) file->Get(name+"_ResoVSPt"+nameSuffix[i]);
2266 // hResoVSEta[i] = (TH2F *) file->Get(name+"_ResoVSEta"+nameSuffix[i]);
2267 // hResoVSPhi[i] = (TH2F *) file->Get(name+"_ResoVSPhi"+nameSuffix[i]);
2268 // }
2269 // }
2270 
2271  protected:
2272  std::map<TString, TH1*> mapHisto_;
2273  TString nameSuffix_[2];
2274  std::unique_ptr<HDelta> muMinus;
2275  std::unique_ptr<HDelta> muPlus;
2276 };
2277 
2278 #endif
#define LogDebug(id)
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...
TH1F * hCotgTheta_
Definition: Histograms.h:517
TProfile * getProfile()
Definition: Histograms.h:165
int getXindex(const double &x) const
Definition: Histograms.h:1966
HTH1D(TFile *outputFile, const TString &name, const TString &title, const int xBins, const double &xMin, const double &xMax)
Definition: Histograms.h:175
const double Pi
void Write() override
Definition: Histograms.h:565
void Clear() override
Definition: Histograms.h:2125
void Fill(const reco::Particle::LorentzVector &recoP1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const reco::Particle::LorentzVector &genP2) override
Definition: Histograms.h:2047
TProfile * hMassVSPhi_prof_
Definition: Histograms.h:730
virtual void Fill(const CLHEP::HepLorentzVector &momentum, const int charge, const double &weight=1.)
Definition: Histograms.h:74
HLikelihoodVSPart(const TString &name)
Definition: Histograms.h:1390
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
Definition: Histograms.h:649
~HFunctionResolution() override
Definition: Histograms.h:1509
TFile * outputFile_
Definition: Histograms.h:123
HPartVSEta(const TString &name, const double &minMass=0., const double &maxMass=100., const double &maxPt=100.)
Definition: Histograms.h:527
void Write() override
Definition: Histograms.h:1435
Histograms(TFile *outputFile, const TString &name)
Definition: Histograms.h:43
virtual void SetYTitle(const TString &title)
Definition: Histograms.h:194
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.
Definition: Histograms.h:261
HPartVSPt(const TString &name)
Definition: Histograms.h:747
std::unique_ptr< HDelta > muPlus
Definition: Histograms.h:2275
~HCovarianceVSParts() override
Definition: Histograms.h:2036
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
Definition: Histograms.h:1251
void Clear() override
Definition: Histograms.h:578
TProfile2D * hMassVSPt_
Definition: Histograms.h:1131
static const char dir_[]
TProfile * hPtVSEta_prof_
Definition: Histograms.h:592
TH2F * hPtVSPhi_
Definition: Histograms.h:728
void Write() override
Definition: Histograms.h:184
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const CLHEP::HepLorentzVector &momentumRes, const double &weight=1.)
Definition: Histograms.h:101
TProfile * hPtVsPhiNeg_
Definition: Histograms.h:410
void Fill(const reco::Particle::LorentzVector &p4, const int charge, const double &weight=1.) override
Definition: Histograms.h:328
TH2F * hMassVSEtaMinus_
Definition: Histograms.h:1048
virtual double Get(const reco::Particle::LorentzVector &recoP1, const TString &covarianceName)
Definition: Histograms.h:107
HMassResolutionVSPart(TFile *outputFile, const TString &name)
Definition: Histograms.h:2148
~HLikelihoodVSPart() override
Definition: Histograms.h:1413
TDirectory * dir2D_
Definition: Histograms.h:1768
void Clear() override
Definition: Histograms.h:2248
double Get(const double &x, const double &y) const
Definition: Histograms.h:1907
HLikelihoodVSPart(const TString &name, TFile *file)
Definition: Histograms.h:1403
~HMassVSPartProfile() override
Definition: Histograms.h:1089
virtual void SetXTitle(const TString &title)
Definition: Histograms.h:191
virtual void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const reco::Particle::LorentzVector &p4Res, const double &weight=1.)
Definition: Histograms.h:97
std::unique_ptr< HDelta > muMinus
Definition: Histograms.h:2274
void Clear() override
Definition: Histograms.h:152
~HDelta() override
Definition: Histograms.h:460
TH3F * hMassVSEtaPhiMinus_
Definition: Histograms.h:1046
#define nullptr
void Fill(const double &x, const double &y, const double &a, const double &b) override
Definition: Histograms.h:1893
TH2F * hResoVSPt_Endc_17_
Definition: Histograms.h:1361
TH2F * hMassVSEta_
Definition: Histograms.h:1039
void Write() override
Definition: Histograms.h:2116
double covariance()
Definition: Histograms.h:1797
void Clear() override
Definition: Histograms.h:1012
Definition: weight.py:1
TH1F * hEtaSign_
Definition: Histograms.h:514
TProfile * hMassVSPt_prof_
Definition: Histograms.h:787
A wrapper for the TH1D histogram to allow it to be put inside the same map as all the other classes i...
Definition: Histograms.h:172
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.)
Definition: Histograms.h:66
TH2D * histoCovariance_
Definition: Histograms.h:1972
void Write() override
Definition: Histograms.h:1296
void Clear() override
Definition: Histograms.h:1955
TH3F * hMassVSPhiPlusPhiMinus_
Definition: Histograms.h:1053
~HParticle() override
Definition: Histograms.h:308
~HCovarianceVSxy() override
Definition: Histograms.h:1870
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
Definition: Histograms.h:2177
void Fill(unsigned int number) override
Definition: Histograms.h:355
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
Definition: Histograms.h:1532
TProfile * hResoVSPhi_prof_
Definition: Histograms.h:1376
TProfile2D * hMassVSPhiPlus_
Definition: Histograms.h:1133
TProfile * hCurvVSEta_prof_
Definition: Histograms.h:593
int getYindex(const double &y) const
Definition: Histograms.h:1632
void Fill(const CLHEP::HepLorentzVector &p1, const reco::Particle::LorentzVector &p2) override
Definition: Histograms.h:474
TH2F * hResoVSPhiMinus_
Definition: Histograms.h:1374
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)
Definition: Histograms.h:86
TProfile * hLikeVSEta_prof_
Definition: Histograms.h:1458
TH1D * operator->()
Definition: Histograms.h:197
~HResolutionVSPart() override
Definition: Histograms.h:1222
TProfile * hPtVSPhi_prof_
Definition: Histograms.h:731
A wrapper for the TProfile histogram to allow it to be put inside the same map as all the other class...
Definition: Histograms.h:203
TH2F * hResoVSPhiPlus_
Definition: Histograms.h:1375
TH1F * hPt_
Definition: Histograms.h:402
TH2F * hPtVsEta_
Definition: Histograms.h:403
virtual ~Histograms()
Definition: Histograms.h:57
HParticle(const TString &name, TFile *file)
Definition: Histograms.h:287
TH1F * hPhi_
Definition: Histograms.h:414
virtual void Clear()=0
TString name_
Definition: Histograms.h:122
~HPartVSPt() override
Definition: Histograms.h:756
TProfile * hMassVSCosThetaCS_prof
Definition: Histograms.h:1056
virtual void SetYTitle(const TString &title)
Definition: Histograms.h:226
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)
Definition: Histograms.h:1821
TProfile * hResoVSPt_Endc_17_prof_
Definition: Histograms.h:1367
virtual void Fill(const double &x, const double &y)
Definition: Histograms.h:95
TProfile * hResoVSPt_Endc_20_prof_
Definition: Histograms.h:1642
TH3F * hMassVSEtaPlusEtaMinus_
Definition: Histograms.h:1054
TH1F * hEta_
Definition: Histograms.h:513
std::map< TString, TH1 * > mapHisto_
Definition: Histograms.h:2272
void Fill(const double &x, const double &y) override
Definition: Histograms.h:143
virtual void SetYTitle(const TString &title)
Definition: Histograms.h:160
void Clear() override
Definition: Histograms.h:381
TH1D * tH1D_
Definition: Histograms.h:199
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
Definition: Histograms.h:765
TProfile * hResoVSPt_Ovlap_prof_
Definition: Histograms.h:1370
void Clear() override
Definition: Histograms.h:1123
TH2F * hMassVSPhiCS_
Definition: Histograms.h:1043
TProfile * hCurvVsEtaPos_
Definition: Histograms.h:406
TProfile * hResoVSPt_Endc_17_prof_
Definition: Histograms.h:1641
TH2F * hPtVSEta_
Definition: Histograms.h:589
void Write() override
Definition: Histograms.h:1919
void Clear() override
Definition: Histograms.h:1612
void Write() override
Definition: Histograms.h:770
TH2F * hMassVSEtaPlus_
Definition: Histograms.h:1047
void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2) override
Definition: Histograms.h:469
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
TProfile * hResoVSPt_Bar_prof_
Definition: Histograms.h:1640
void Clear() override
Definition: Histograms.h:780
T sqrt(T t)
Definition: SSEVec.h:18
HDelta(const TString &name, TFile *file)
Definition: Histograms.h:450
double p4[4]
Definition: TauolaWrapper.h:92
TH2D * operator->()
Definition: Histograms.h:164
TProfile * hResoVSPt_prof_
Definition: Histograms.h:1639
virtual void SetXTitle(const TString &title)
Definition: Histograms.h:156
TDirectory * dir_
Definition: Histograms.h:1767
std::map< TString, HCovarianceVSxy * > mapHisto_
Definition: Histograms.h:2132
TDirectory * histoDir_
Definition: Histograms.h:124
void fill(const double &x, const double &y)
Definition: Histograms.h:1791
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
Definition: Histograms.h:761
TH1F * hDeltaR_
Definition: Histograms.h:518
TProfile * hCurvVsEtaNeg_
Definition: Histograms.h:405
void Fill(const double &x, const double &y) override
Definition: Histograms.h:213
HMassVSPartProfile(const TString &name, const double &minMass=0., const double &maxMass=150., const double maxPt=100.)
Definition: Histograms.h:1070
~HTH1D() override
Definition: Histograms.h:178
TProfile * hPtVsPhiPos_
Definition: Histograms.h:411
void Clear() override
Definition: Histograms.h:220
Int_t Fill(Double_t x, Double_t y) override
Definition: Histograms.h:1739
TProfile * tProfile_
Definition: Histograms.h:231
virtual void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2)
Definition: Histograms.h:63
double getN()
Definition: Histograms.h:1806
TDirectory * diffDir_
Definition: Histograms.h:1769
void Fill(const CLHEP::HepLorentzVector &momentum, const int charge, const double &weight=1.) override
Definition: Histograms.h:333
void Write() override
Definition: Histograms.h:216
HFunctionResolutionVarianceCheck(TFile *outputFile, const TString &name, const double ptMax=200)
Definition: Histograms.h:1658
TProfile * hLikeVSPhi_prof_
Definition: Histograms.h:1459
HParticle(const TString &name, const double &minMass=0., const double &maxMass=200., const double &maxPt=100.)
Definition: Histograms.h:239
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const int charge, const double &weight=1.) override
Definition: Histograms.h:1096
TH2F * hMassVSPhi_
Definition: Histograms.h:729
double productXY_
Definition: Histograms.h:1808
TProfile * hResoVSPt_Ovlap_prof_
Definition: Histograms.h:1644
const int mu
Definition: Constants.h:22
void Write() override
Definition: Histograms.h:360
TProfile * hCurvVsPhiPos_
Definition: Histograms.h:408
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.) override
Definition: Histograms.h:961
virtual void Fill(const unsigned int number)
Definition: Histograms.h:80
TProfile * tProfile_
Definition: Histograms.h:168
void Write() override
Definition: Histograms.h:2236
void Clear() override
Definition: Histograms.h:503
HMassVSPart(const TString &name, TFile *file)
Definition: Histograms.h:837
double p2[4]
Definition: TauolaWrapper.h:90
TH2F * hResoVSPt_Ovlap_
Definition: Histograms.h:1364
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &genValue, const double recValue, const int charge)
Definition: Histograms.h:78
const int mubar
Definition: Constants.h:23
TH1F * hNumber_
Definition: Histograms.h:417
TProfile * hResoVSPt_Endc_24_prof_
Definition: Histograms.h:1643
TProfile * hCurvVsPhiNeg_
Definition: Histograms.h:407
TH2F * hMassVSPt_
Definition: Histograms.h:786
virtual void Fill(const reco::Particle::LorentzVector &p1, const reco::Particle::LorentzVector &p2, const int charge, const double &weight=1.)
Definition: Histograms.h:64
virtual void SetXTitle(const TString &title)
Definition: Histograms.h:223
virtual void Fill(const reco::Particle::LorentzVector &p4, const int charge, const double &weight=1.)
Definition: Histograms.h:72
TH1F * resoHisto_
Definition: Histograms.h:1772
TH2F * hMassVSCosThetaCS_
Definition: Histograms.h:1042
void Fill(const reco::Particle::LorentzVector &p4, const double &likeValue) override
Definition: Histograms.h:1422
~HMassResolutionVSPart() override
Definition: Histograms.h:2170
void Clear() override
Definition: Histograms.h:712
TH2F * hResoVSPt_Endc_24_
Definition: Histograms.h:1363
Covariance ** covariances_
Definition: Histograms.h:1973
TH3F * hMassVSEtaPhiPlus_
Definition: Histograms.h:1045
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
Definition: Histograms.h:552
HCovarianceVSParts(TFile *outputFile, const TString &name, const double &ptMax)
Definition: Histograms.h:1988
TProfile * hResoVSPhi_prof_
Definition: Histograms.h:1649
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)
Definition: Histograms.h:81
TProfile * hResoVSPt_prof_
Definition: Histograms.h:1365
void Write() override
Definition: Histograms.h:492
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const reco::Particle::LorentzVector &p4Res, const double &weight=1.) override
Definition: Histograms.h:884
HCovarianceVSxy(TFile *inputFile, const TString &name, const TString &dirName)
Contructor to read histograms from file.
Definition: Histograms.h:1853
TProfile * hMassVSEta_prof_
Definition: Histograms.h:591
#define N
Definition: blowfish.cc:9
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:152
void Write() override
Definition: Histograms.h:1116
TProfile * hResoVSPhiMinus_prof_
Definition: Histograms.h:1647
~HMassVSPart() override
Definition: Histograms.h:860
virtual void Fill(const reco::Particle::LorentzVector &recoP1, const reco::Particle::LorentzVector &genP1, const reco::Particle::LorentzVector &recoP2, const reco::Particle::LorentzVector &genP2)
Definition: Histograms.h:91
TProfile2D * hMassVSEta_
Definition: Histograms.h:1132
virtual void SetWeight(double weight)
Definition: Histograms.h:112
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge)
Definition: Histograms.h:77
TH2F * hResoVSPt_Endc_20_
Definition: Histograms.h:1362
TH2F * hMassVSPhiMinus_
Definition: Histograms.h:1041
TH2F * hMassVSEtaPlusMinusDiff_
Definition: Histograms.h:1051
double b
Definition: hdecay.h:120
virtual void Fill(const CLHEP::HepLorentzVector &p, const double &likeValue)
Definition: Histograms.h:79
void Write() override
Definition: Histograms.h:664
TH2F * hMassVSPhiPlus_
Definition: Histograms.h:1040
int getXindex(const double &x) const
Definition: Histograms.h:1629
void Fill(const CLHEP::HepLorentzVector &momentum, const double &weight=1.) override
Definition: Histograms.h:556
void Clear() override
Definition: Histograms.h:188
void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.) override
Definition: Histograms.h:645
TProfile * hResoVSPt_Endc_24_prof_
Definition: Histograms.h:1369
TH1F * hTheta_
Definition: Histograms.h:516
void Write() override
Definition: Histograms.h:147
virtual void Write()=0
~HFunctionResolutionVarianceCheck() override
Definition: Histograms.h:1673
void Clear() override
Definition: Histograms.h:1327
~HTProfile() override
Definition: Histograms.h:210
double sumY_
Definition: Histograms.h:1810
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)
Definition: Histograms.h:134
A wrapper for the TH2D histogram to allow it to be put inside the same map as all the other classes i...
Definition: Histograms.h:131
double p1[4]
Definition: TauolaWrapper.h:89
HMassVSPartProfile(const TString &name, TFile *file)
Definition: Histograms.h:1081
HDelta(TFile *outputFile, const TString &name)
Definition: Histograms.h:438
virtual TString GetName()
Definition: Histograms.h:116
TH2F * hMassVSPt_
Definition: Histograms.h:1038
void Fill(const double &x, const double &y) override
Definition: Histograms.h:181
int getYindex(const double &y) const
Definition: Histograms.h:1969
HTProfile(TFile *outputFile, const TString &name, const TString &title, const int xBins, const double &xMin, const double &xMax, const double &yMin, const double &yMax)
Definition: Histograms.h:206
void Write() override
Definition: Histograms.h:986
TH1D *** histoVarianceCheck_
Definition: Histograms.h:1978
double a
Definition: hdecay.h:121
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2, const int charge, const double &weight=1.) override
Definition: Histograms.h:1101
TProfile * hResoVSEta_prof_
Definition: Histograms.h:1645
TH2F * hMassVSPhiPlusMinusDiff_
Definition: Histograms.h:1050
HDelta(const TString &name)
Definition: Histograms.h:426
void Fill(const CLHEP::HepLorentzVector &momentum, const double &likeValue) override
Definition: Histograms.h:1426
def canvas(sub, attr)
Definition: svgfig.py:482
TProfile * hResoVSPt_Bar_prof_
Definition: Histograms.h:1366
TProfile2D * hMassVSPhiMinus_
Definition: Histograms.h:1134
TProfile * hLikeVSPt_prof_
Definition: Histograms.h:1457
void Write() override
Definition: Histograms.h:1573
TH1F * hMass_
Definition: Histograms.h:415
void Fill(const reco::Particle::LorentzVector &p4, const double &resValue, const int charge) override
Definition: Histograms.h:1682
static double deltaPhi(const double &phi1, const double &phi2)
TProfile * hResoVSPt_Endc_20_prof_
Definition: Histograms.h:1368
HFunctionResolution(TFile *outputFile, const TString &name, const double &ptMax=100, const int totBinsY=300)
Definition: Histograms.h:1474
virtual void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2)
Definition: Histograms.h:65
HPartVSPhi(const TString &name)
Definition: Histograms.h:602
dbl *** dir
Definition: mlp_gen.cc:35
~HPartVSPhi() override
Definition: Histograms.h:630
TProfile * operator->()
Definition: Histograms.h:229
TH1F * hEta_
Definition: Histograms.h:413
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
Definition: Histograms.h:2188
TH2F * hMassVSEta_
Definition: Histograms.h:590
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)
Definition: Histograms.h:1142
TProfile * hResoVSEta_prof_
Definition: Histograms.h:1373
TH2F * histo2D_
Definition: Histograms.h:1771
TProfile * hResoVSPhiPlus_prof_
Definition: Histograms.h:1648
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
void Fill(const CLHEP::HepLorentzVector &momentum1, const CLHEP::HepLorentzVector &momentum2) override
Definition: Histograms.h:478
A set of histograms for resolution.
Definition: Histograms.h:1139
Histograms(const TString &name)
Definition: Histograms.h:42
HCovarianceVSParts(const TString &inputFileName, const TString &name)
Constructor reading the histograms from file.
Definition: Histograms.h:2014
void Fill(const reco::Particle::LorentzVector &p41, const reco::Particle::LorentzVector &p42, const int charge, const double &weight=1.) override
Definition: Histograms.h:878
~HResolution() override
Definition: Histograms.h:1734
double Get(const reco::Particle::LorentzVector &recoP1, const TString &covarianceName) override
Definition: Histograms.h:2043
void Clear() override
Definition: Histograms.h:1444
virtual void Fill(const CLHEP::HepLorentzVector &p1, const reco::Particle::LorentzVector &p2)
Definition: Histograms.h:67
yBin
Definition: cuy.py:893
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
~HTH2D() override
Definition: Histograms.h:139
double ** resoVsPtEta_
Definition: Histograms.h:1637
double theWeight_
Definition: Histograms.h:121
HResolutionVSPart(const TString &name, TFile *file)
Definition: Histograms.h:1192
~HPartVSEta() override
Definition: Histograms.h:544
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)
Definition: Histograms.h:1717
HMassVSPart(const TString &name, const double &minMass=0., const double &maxMass=150., const double maxPt=100.)
Definition: Histograms.h:795
virtual void Fill(const double &x, const double &y, const double &a, const double &b)
Definition: Histograms.h:96
TH2D * tH2d_
Definition: Histograms.h:167
TH1F * hPhi_
Definition: Histograms.h:515
Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0) override
Definition: Histograms.h:1744
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.
Definition: Histograms.h:896
double sumX_
Definition: Histograms.h:1809
TProfile * diffHisto_
Definition: Histograms.h:1770
virtual void Fill(const reco::Particle::LorentzVector &p4, const double &weight=1.)
Definition: Histograms.h:68