Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:39 

0001 #ifndef Histograms_H
0002 #define Histograms_H
0003 
0004 /** \class Histograms
0005  *  Collection of histograms for GLB muon analysis
0006  *
0007  *  \author S. Bolognesi - INFN Torino / T.Dorigo - INFN Padova
0008  * revised S. Casasso, E. Migliore - UniTo & INFN Torino 
0009 */
0010 
0011 #include <CLHEP/Vector/LorentzVector.h>
0012 #include "DataFormats/Math/interface/LorentzVector.h"
0013 #include "DataFormats/MuonReco/interface/Muon.h"
0014 #include "MuScleFitUtils.h"
0015 
0016 #include "TH1D.h"
0017 #include "TH1F.h"
0018 #include "TH2F.h"
0019 #include "TH3F.h"
0020 #include "TFile.h"
0021 #include "TString.h"
0022 #include "TProfile.h"
0023 #include "TProfile2D.h"
0024 #include "TF1.h"
0025 #include "TGraphErrors.h"
0026 #include "TFile.h"
0027 #include "TSystem.h"
0028 #include "TCanvas.h"
0029 
0030 #include "TLorentzVector.h"
0031 #include <memory>
0032 
0033 #include "TMath.h"
0034 #include <iostream>
0035 #include <string>
0036 #include <vector>
0037 
0038 class Histograms {
0039 public:
0040   // Constructor
0041   // -----------
0042   Histograms() : theWeight_(1), histoDir_(nullptr){};
0043   Histograms(const TString& name) : theWeight_(1), name_(name), histoDir_(nullptr) {}
0044   Histograms(TFile* outputFile, const TString& name)
0045       : theWeight_(1), name_(name), outputFile_(outputFile), histoDir_(outputFile->GetDirectory(name)) {
0046     if (histoDir_ == nullptr) {
0047       histoDir_ = outputFile->mkdir(name);
0048     }
0049     histoDir_->cd();
0050   }
0051 
0052   // Destructor
0053   // ----------
0054   virtual ~Histograms(){};
0055 
0056   // Operations
0057   // ----------
0058   //   virtual void Fill( const reco::Particle::LorentzVector & p4 ) {};
0059   //   virtual void Fill( const CLHEP::HepLorentzVector & momentum ) {};
0060   virtual void Fill(const reco::Particle::LorentzVector& p1, const reco::Particle::LorentzVector& p2){};
0061   virtual void Fill(const reco::Particle::LorentzVector& p1,
0062                     const reco::Particle::LorentzVector& p2,
0063                     const int charge,
0064                     const double& weight = 1.){};
0065   virtual void Fill(const CLHEP::HepLorentzVector& momentum1, const CLHEP::HepLorentzVector& momentum2){};
0066   virtual void Fill(const CLHEP::HepLorentzVector& momentum1,
0067                     const CLHEP::HepLorentzVector& momentum2,
0068                     const int charge,
0069                     const double& weight = 1.){};
0070   virtual void Fill(const CLHEP::HepLorentzVector& p1, const reco::Particle::LorentzVector& p2){};
0071   virtual void Fill(const reco::Particle::LorentzVector& p4, const double& weight = 1.){};
0072 
0073   //virtual void Fill( const CLHEP::HepLorentzVector & momentum, const double & weight ) {};
0074   //------
0075   virtual void Fill(const reco::Particle::LorentzVector& p4, const int charge, const double& weight = 1.){};
0076 
0077   virtual void Fill(const CLHEP::HepLorentzVector& momentum, const int charge, const double& weight = 1.){};
0078   //------
0079   // virtual void Fill( const reco::Particle::LorentzVector & p4, const double & likeValue ) {};
0080   virtual void Fill(const reco::Particle::LorentzVector& p4, const double& resValue, const int charge){};
0081   virtual void Fill(const reco::Particle::LorentzVector& p4,
0082                     const double& genValue,
0083                     const double recValue,
0084                     const int charge){};
0085   virtual void Fill(const CLHEP::HepLorentzVector& p, const double& likeValue){};
0086   virtual void Fill(const unsigned int number){};
0087   virtual void Fill(const reco::Particle::LorentzVector& recoP1,
0088                     const int charge1,
0089                     const reco::Particle::LorentzVector& genP1,
0090                     const reco::Particle::LorentzVector& recoP2,
0091                     const int charge2,
0092                     const reco::Particle::LorentzVector& genP2,
0093                     const double& recoMass,
0094                     const double& genMass){};
0095   virtual void Fill(const reco::Particle::LorentzVector& recoP1,
0096                     const int charge1,
0097                     // const reco::Particle::LorentzVector & genP1,
0098                     const reco::Particle::LorentzVector& recoP2,
0099                     const int charge2,
0100                     // const reco::Particle::LorentzVector & genP2,
0101                     const double& recoMass,
0102                     const double& genMass){};
0103   virtual void Fill(const reco::Particle::LorentzVector& recoP1,
0104                     const reco::Particle::LorentzVector& genP1,
0105                     const reco::Particle::LorentzVector& recoP2,
0106                     const reco::Particle::LorentzVector& genP2){};
0107   virtual void Fill(const double& x, const double& y){};
0108   virtual void Fill(const double& x, const double& y, const double& a, const double& b){};
0109   virtual void Fill(const reco::Particle::LorentzVector& p41,
0110                     const reco::Particle::LorentzVector& p42,
0111                     const reco::Particle::LorentzVector& p4Res,
0112                     const double& weight = 1.){};
0113   virtual void Fill(const CLHEP::HepLorentzVector& momentum1,
0114                     const CLHEP::HepLorentzVector& momentum2,
0115                     const CLHEP::HepLorentzVector& momentumRes,
0116                     const double& weight = 1.){};
0117 
0118   virtual double Get(const reco::Particle::LorentzVector& recoP1, const TString& covarianceName) { return 0.; };
0119 
0120   virtual void Write() = 0;
0121   virtual void Clear() = 0;
0122 
0123   virtual void SetWeight(double weight) { theWeight_ = weight; }
0124 
0125   virtual TString GetName() { return name_; }
0126 
0127 protected:
0128   double theWeight_;
0129   TString name_;
0130   TFile* outputFile_;
0131   TDirectory* histoDir_;
0132 
0133 private:
0134 };
0135 
0136 /// A wrapper for the TH2D histogram to allow it to be put inside the same map as all the other classes in this file
0137 class HTH2D : public Histograms {
0138 public:
0139   HTH2D(TFile* outputFile,
0140         const TString& name,
0141         const TString& title,
0142         const TString& dirName,
0143         const int xBins,
0144         const double& xMin,
0145         const double& xMax,
0146         const int yBins,
0147         const double& yMin,
0148         const double& yMax)
0149       : Histograms(outputFile, dirName),
0150         tH2d_(new TH2D(name, title, xBins, xMin, xMax, yBins, yMin, yMax)),
0151         tProfile_(new TProfile(name + "Prof", title + " profile", xBins, xMin, xMax, yMin, yMax)) {}
0152   ~HTH2D() override {
0153     delete tH2d_;
0154     delete tProfile_;
0155   }
0156   void Fill(const double& x, const double& y) override {
0157     tH2d_->Fill(x, y);
0158     tProfile_->Fill(x, y);
0159   }
0160   void Write() override {
0161     if (histoDir_ != nullptr)
0162       histoDir_->cd();
0163     tH2d_->Write();
0164     tProfile_->Write();
0165   }
0166   void Clear() override {
0167     tH2d_->Clear();
0168     tProfile_->Clear();
0169   }
0170   virtual void SetXTitle(const TString& title) {
0171     tH2d_->GetXaxis()->SetTitle(title);
0172     tProfile_->GetXaxis()->SetTitle(title);
0173   }
0174   virtual void SetYTitle(const TString& title) {
0175     tH2d_->GetYaxis()->SetTitle(title);
0176     tProfile_->GetYaxis()->SetTitle(title);
0177   }
0178   TH2D* operator->() { return tH2d_; }
0179   TProfile* getProfile() { return tProfile_; }
0180 
0181 protected:
0182   TH2D* tH2d_;
0183   TProfile* tProfile_;
0184 };
0185 
0186 /// A wrapper for the TH1D histogram to allow it to be put inside the same map as all the other classes in this file
0187 class HTH1D : public Histograms {
0188 public:
0189   HTH1D(TFile* outputFile,
0190         const TString& name,
0191         const TString& title,
0192         const int xBins,
0193         const double& xMin,
0194         const double& xMax)
0195       : Histograms(outputFile, name), tH1D_(new TH1D(name, title, xBins, xMin, xMax)) {}
0196   ~HTH1D() override { delete tH1D_; }
0197   void Fill(const double& x, const double& y) override { tH1D_->Fill(x, y); }
0198   void Write() override {
0199     if (histoDir_ != nullptr)
0200       histoDir_->cd();
0201     tH1D_->Write();
0202   }
0203   void Clear() override { tH1D_->Clear(); }
0204   virtual void SetXTitle(const TString& title) { tH1D_->GetXaxis()->SetTitle(title); }
0205   virtual void SetYTitle(const TString& title) { tH1D_->GetYaxis()->SetTitle(title); }
0206   TH1D* operator->() { return tH1D_; }
0207 
0208 protected:
0209   TH1D* tH1D_;
0210 };
0211 
0212 /// A wrapper for the TProfile histogram to allow it to be put inside the same map as all the other classes in this file
0213 class HTProfile : public Histograms {
0214 public:
0215   HTProfile(TFile* outputFile,
0216             const TString& name,
0217             const TString& title,
0218             const int xBins,
0219             const double& xMin,
0220             const double& xMax,
0221             const double& yMin,
0222             const double& yMax)
0223       : Histograms(outputFile, name),
0224         tProfile_(new TProfile(name + "Prof", title + " profile", xBins, xMin, xMax, yMin, yMax)) {}
0225   ~HTProfile() override { delete tProfile_; }
0226   void Fill(const double& x, const double& y) override { tProfile_->Fill(x, y); }
0227   void Write() override {
0228     if (histoDir_ != nullptr)
0229       histoDir_->cd();
0230     tProfile_->Write();
0231   }
0232   void Clear() override { tProfile_->Clear(); }
0233   virtual void SetXTitle(const TString& title) { tProfile_->GetXaxis()->SetTitle(title); }
0234   virtual void SetYTitle(const TString& title) { tProfile_->GetYaxis()->SetTitle(title); }
0235   TProfile* operator->() { return tProfile_; }
0236 
0237 protected:
0238   TProfile* tProfile_;
0239 };
0240 
0241 // -----------------------------------------------------
0242 // A set of histograms of particle kinematical variables
0243 // -----------------------------------------------------
0244 class HParticle : public Histograms {
0245 public:
0246   HParticle(const TString& name, const double& minMass = 0., const double& maxMass = 200., const double& maxPt = 100.)
0247       : Histograms(name),
0248         // Kinematical variables
0249         hPt_(new TH1F(name + "_Pt", "transverse momentum", 100, 0, maxPt)),
0250         hPtVsEta_(new TH2F(name + "_PtVsEta", "transverse momentum vs #eta", 100, 0, maxPt, 100, -3.0, 3.0)),
0251 
0252         hCurvVsEtaNeg_(new TProfile(name + "_CurvVsEtaNeg", "q/pT vs #eta neg.", 64, -3.2, 3.2, -1., 0.)),
0253         hCurvVsEtaPos_(new TProfile(name + "_CurvVsEtaPos", "q/pT vs #eta pos.", 64, -3.2, 3.2, 0., 1.)),
0254         hCurvVsPhiNeg_(new TProfile(name + "_CurvVsPhiNeg", "q/pT vs #phi neg.", 32, -3.2, 3.2, -1., 0.)),
0255         hCurvVsPhiPos_(new TProfile(name + "_CurvVsPhiPos", "q/pT vs #phi pos.", 32, -3.2, 3.2, 0., 1.)),
0256 
0257         hPtVsPhiNeg_(new TProfile(name + "_PtVsPhiNeg", "pT vs #phi neg.", 32, -3.2, 3.2, 0., 100)),
0258         hPtVsPhiPos_(new TProfile(name + "_PtVsPhiPos", "pT vs #phi pos.", 32, -3.2, 3.2, 0., 100)),
0259 
0260         hEta_(new TH1F(name + "_Eta", "pseudorapidity", 64, -3.2, 3.2)),
0261         hPhi_(new TH1F(name + "_Phi", "phi angle", 64, -3.2, 3.2)),
0262         hMass_(new TH1F(name + "_Mass", "mass", 10000, minMass, maxMass)),
0263         hNumber_(new TH1F(name + "_Number", "number", 20, -0.5, 19.5)) {}
0264 
0265   /// Constructor that puts the histograms inside a TDirectory
0266   HParticle(TFile* outputFile,
0267             const TString& name,
0268             const double& minMass = 0.,
0269             const double& maxMass = 200.,
0270             const double& maxPt = 100.)
0271       : Histograms(outputFile, name) {
0272     // Kinematical variables
0273     hPt_ = new TH1F(name + "_Pt", "transverse momentum", 100, 0, maxPt);
0274     hPtVsEta_ = new TH2F(name + "_PtVsEta", "transverse momentum vs #eta", 100, 0, maxPt, 100, -3.0, 3.0);
0275 
0276     hPtVsEta_ = new TH2F(name + "_PtVsEta", "transverse momentum vs #eta", 100, 0, maxPt, 100, -3.0, 3.0);
0277 
0278     hCurvVsEtaNeg_ = new TProfile(name + "_CurvVsEtaNeg", "q/pT vs #eta neg.", 100, -3.0, 3.0, -1., 0.);
0279     hCurvVsEtaPos_ = new TProfile(name + "_CurvVsEtaPos", "q/pT vs #eta pos.", 100, -3.0, 3.0, 0., 1.);
0280     hCurvVsPhiNeg_ = new TProfile(name + "_CurvVsPhiNeg", "q/pT vs #phi neg.", 32, -3.2, 3.2, -1., 0.);
0281     hCurvVsPhiPos_ = new TProfile(name + "_CurvVsPhiPos", "q/pT vs #phi pos.", 32, -3.2, 3.2, 0., 1.);
0282 
0283     hPtVsPhiNeg_ = new TProfile(name + "_PtVsPhiNeg", "pT vs #phi neg.", 32, -3.2, 3.2, 0., 100);
0284     hPtVsPhiPos_ = new TProfile(name + "_PtVsPhiPos", "pT vs #phi pos.", 32, -3.2, 3.2, 0., 100);
0285 
0286     //hPtVSPhi_prof_ = new TProfile (name+"_PtVSPhi_prof", "pt vs phi angle",12, -3.2, 3.2, 0, 200);
0287 
0288     hEta_ = new TH1F(name + "_Eta", "pseudorapidity", 64, -3.2, 3.2);
0289     hPhi_ = new TH1F(name + "_Phi", "phi angle", 64, -3.2, 3.2);
0290     hMass_ = new TH1F(name + "_Mass", "mass", 40000, minMass, maxMass);
0291     hNumber_ = new TH1F(name + "_Number", "number", 20, -0.5, 19.5);
0292   }
0293 
0294   HParticle(const TString& name, TFile* file)
0295       : Histograms(name),
0296         hPt_((TH1F*)file->Get(name_ + "_Pt")),
0297         hPtVsEta_((TH2F*)file->Get(name_ + "_PtVsEta")),
0298 
0299         hCurvVsEtaNeg_((TProfile*)file->Get(name_ + "_CurvVsEtaNeg")),
0300         hCurvVsEtaPos_((TProfile*)file->Get(name_ + "_CurvVsEtaPos")),
0301         hCurvVsPhiNeg_((TProfile*)file->Get(name_ + "_CurvVsPhiNeg")),
0302         hCurvVsPhiPos_((TProfile*)file->Get(name_ + "_CurvVsPhiPos")),
0303 
0304         hPtVsPhiNeg_((TProfile*)file->Get(name_ + "_PtVsPhiNeg")),
0305         hPtVsPhiPos_((TProfile*)file->Get(name_ + "_PtVsPhiPos")),
0306 
0307         hEta_((TH1F*)file->Get(name_ + "_Eta")),
0308         hPhi_((TH1F*)file->Get(name_ + "_Phi")),
0309         hMass_((TH1F*)file->Get(name_ + "_Mass")),
0310         //hMass_fine_ = (TH1F *) file->Get(name_+"_Mass_fine");
0311         hNumber_((TH1F*)file->Get(name_ + "_Number")) {}
0312 
0313   ~HParticle() override {
0314     delete hPt_;
0315     delete hPtVsEta_;
0316 
0317     delete hCurvVsEtaNeg_;
0318     delete hCurvVsEtaPos_;
0319     delete hCurvVsPhiNeg_;
0320     delete hCurvVsPhiPos_;
0321 
0322     delete hPtVsPhiNeg_;
0323     delete hPtVsPhiPos_;
0324 
0325     delete hEta_;
0326     delete hPhi_;
0327     delete hMass_;
0328     // delete hMass_fine_;
0329     delete hNumber_;
0330   }
0331 
0332   void Fill(const reco::Particle::LorentzVector& p4, const int charge, const double& weight = 1.) override {
0333     Fill(CLHEP::HepLorentzVector(p4.x(), p4.y(), p4.z(), p4.t()), charge, weight);
0334   }
0335 
0336   void Fill(const CLHEP::HepLorentzVector& momentum, const int charge, const double& weight = 1.) override {
0337     hPt_->Fill(momentum.perp(), weight);
0338     hPtVsEta_->Fill(momentum.perp(), momentum.eta(), weight);
0339 
0340     //   std::cout<< "charge-> " <<charge<<std::endl;
0341     if (charge < 0)
0342       hCurvVsEtaNeg_->Fill(momentum.eta(), charge / (momentum.perp()), weight);
0343     if (charge > 0)
0344       hCurvVsEtaPos_->Fill(momentum.eta(), charge / (momentum.perp()), weight);
0345     if (charge < 0)
0346       hCurvVsPhiNeg_->Fill(momentum.phi(), charge / (momentum.perp()), weight);
0347     if (charge > 0)
0348       hCurvVsPhiPos_->Fill(momentum.phi(), charge / (momentum.perp()), weight);
0349 
0350     if (charge < 0)
0351       hPtVsPhiNeg_->Fill(momentum.phi(), momentum.perp(), weight);
0352     if (charge > 0)
0353       hPtVsPhiPos_->Fill(momentum.phi(), momentum.perp(), weight);
0354 
0355     hEta_->Fill(momentum.eta(), weight);
0356     hPhi_->Fill(momentum.phi(), weight);
0357     hMass_->Fill(momentum.m(), weight);
0358     //hMass_fine_->Fill(momentum.m(), weight);
0359   }
0360 
0361   void Fill(unsigned int number) override { hNumber_->Fill(number); }
0362 
0363   void Write() override {
0364     if (histoDir_ != nullptr)
0365       histoDir_->cd();
0366     hPt_->Write();
0367     hPtVsEta_->Write();
0368 
0369     hCurvVsEtaNeg_->Write();
0370     hCurvVsEtaPos_->Write();
0371     hCurvVsPhiNeg_->Write();
0372     hCurvVsPhiPos_->Write();
0373 
0374     hPtVsPhiNeg_->Write();
0375     hPtVsPhiPos_->Write();
0376 
0377     hEta_->Write();
0378     hPhi_->Write();
0379     hMass_->Write();
0380     //hMass_fine_->Write();
0381     hNumber_->Write();
0382   }
0383 
0384   void Clear() override {
0385     hPt_->Clear();
0386     hPtVsEta_->Clear();
0387 
0388     hCurvVsEtaNeg_->Clear();
0389     hCurvVsEtaPos_->Clear();
0390     hCurvVsPhiNeg_->Clear();
0391     hCurvVsPhiPos_->Clear();
0392 
0393     hPtVsPhiNeg_->Clear();
0394     hPtVsPhiPos_->Clear();
0395 
0396     hEta_->Clear();
0397     hPhi_->Clear();
0398     hMass_->Clear();
0399     //hMass_fine_->Clear();
0400     hNumber_->Clear();
0401   }
0402 
0403 protected:
0404   TH1F* hPt_;
0405   TH2F* hPtVsEta_;
0406 
0407   TProfile* hCurvVsEtaNeg_;
0408   TProfile* hCurvVsEtaPos_;
0409   TProfile* hCurvVsPhiNeg_;
0410   TProfile* hCurvVsPhiPos_;
0411 
0412   TProfile* hPtVsPhiNeg_;
0413   TProfile* hPtVsPhiPos_;
0414 
0415   TH1F* hEta_;
0416   TH1F* hPhi_;
0417   TH1F* hMass_;
0418   //TH1F* hMass_fine_;
0419   TH1F* hNumber_;
0420 };
0421 
0422 // ---------------------------------------------------
0423 // A set of histograms for distances between particles
0424 // ---------------------------------------------------
0425 class HDelta : public Histograms {
0426 public:
0427   HDelta(const TString& name)
0428       : Histograms(name),
0429         // Kinematical variables
0430         // ---------------------
0431         hEta_(new TH1F(name + "_DeltaEta", "#Delta#eta", 100, 0, 6)),
0432         hEtaSign_(new TH1F(name + "_DeltaEtaSign", "#Delta#eta with sign", 100, -6, 6)),
0433         hPhi_(new TH1F(name + "_DeltaPhi", "#Delta#phi", 100, 0, 3.2)),
0434         hTheta_(new TH1F(name + "_DeltaTheta", "#Delta#theta", 100, -3.2, 3.2)),
0435         hCotgTheta_(new TH1F(name + "_DeltaCotgTheta", "#Delta Cotg(#theta )", 100, -3.2, 3.2)),
0436         hDeltaR_(new TH1F(name + "_DeltaR", "#Delta R", 400, 0, 4)) {}
0437 
0438   HDelta(TFile* outputFile, const TString& name)
0439       : Histograms(outputFile, name),
0440         // Kinematical variables
0441         // ---------------------
0442         hEta_(new TH1F(name + "_DeltaEta", "#Delta#eta", 100, 0, 6)),
0443         hEtaSign_(new TH1F(name + "_DeltaEtaSign", "#Delta#eta with sign", 100, -6, 6)),
0444         hPhi_(new TH1F(name + "_DeltaPhi", "#Delta#phi", 100, 0, 3.2)),
0445         hTheta_(new TH1F(name + "_DeltaTheta", "#Delta#theta", 100, -3.2, 3.2)),
0446         hCotgTheta_(new TH1F(name + "_DeltaCotgTheta", "#Delta Cotg(#theta )", 100, -3.2, 3.2)),
0447         hDeltaR_(new TH1F(name + "_DeltaR", "#DeltaR", 400, 0, 4)) {}
0448 
0449   HDelta(const TString& name, TFile* file) {
0450     name_ = name;
0451     hEta_ = (TH1F*)file->Get(name + "_DeltaEta");
0452     hEtaSign_ = (TH1F*)file->Get(name + "_DeltaEtaSign");
0453     hPhi_ = (TH1F*)file->Get(name + "_DeltaPhi");
0454     hTheta_ = (TH1F*)file->Get(name + "_DeltaTheta");
0455     hCotgTheta_ = (TH1F*)file->Get(name + "_DeltaCotgTheta");
0456     hDeltaR_ = (TH1F*)file->Get(name + "_DeltaR");
0457   }
0458 
0459   ~HDelta() override {
0460     delete hEta_;
0461     delete hEtaSign_;
0462     delete hPhi_;
0463     delete hTheta_;
0464     delete hCotgTheta_;
0465     delete hDeltaR_;
0466   }
0467 
0468   void Fill(const reco::Particle::LorentzVector& p1, const reco::Particle::LorentzVector& p2) override {
0469     Fill(CLHEP::HepLorentzVector(p1.x(), p1.y(), p1.z(), p1.t()),
0470          CLHEP::HepLorentzVector(p2.x(), p2.y(), p2.z(), p2.t()));
0471   }
0472 
0473   void Fill(const CLHEP::HepLorentzVector& p1, const reco::Particle::LorentzVector& p2) override {
0474     Fill(p1, CLHEP::HepLorentzVector(p2.x(), p2.y(), p2.z(), p2.t()));
0475   }
0476 
0477   void Fill(const CLHEP::HepLorentzVector& momentum1, const CLHEP::HepLorentzVector& momentum2) override {
0478     hEta_->Fill(fabs(momentum1.eta() - momentum2.eta()));
0479     hEtaSign_->Fill(momentum1.eta() - momentum2.eta());
0480     hPhi_->Fill(MuScleFitUtils::deltaPhi(momentum1.phi(), momentum2.phi()));
0481     hTheta_->Fill(momentum1.theta() - momentum2.theta());
0482     // hCotgTheta->Fill(1/(TMath::Tan(momentum1.theta()))-1/(TMath::Tan(momentum2.theta())));
0483     double theta1 = momentum1.theta();
0484     double theta2 = momentum2.theta();
0485     hCotgTheta_->Fill(TMath::Cos(theta1) / TMath::Sin(theta1) - TMath::Cos(theta2) / TMath::Sin(theta2));
0486     hDeltaR_->Fill(sqrt((momentum1.eta() - momentum2.eta()) * (momentum1.eta() - momentum2.eta()) +
0487                         (MuScleFitUtils::deltaPhi(momentum1.phi(), momentum2.phi())) *
0488                             (MuScleFitUtils::deltaPhi(momentum1.phi(), momentum2.phi()))));
0489   }
0490 
0491   void Write() override {
0492     if (histoDir_ != nullptr)
0493       histoDir_->cd();
0494 
0495     hEta_->Write();
0496     hEtaSign_->Write();
0497     hPhi_->Write();
0498     hTheta_->Write();
0499     hCotgTheta_->Write();
0500     hDeltaR_->Write();
0501   }
0502 
0503   void Clear() override {
0504     hEta_->Clear();
0505     hEtaSign_->Clear();
0506     hPhi_->Clear();
0507     hTheta_->Clear();
0508     hDeltaR_->Clear();
0509     hCotgTheta_->Clear();
0510   }
0511 
0512 public:
0513   TH1F* hEta_;
0514   TH1F* hEtaSign_;
0515   TH1F* hPhi_;
0516   TH1F* hTheta_;
0517   TH1F* hCotgTheta_;
0518   TH1F* hDeltaR_;
0519 };
0520 
0521 // ------------------------------------------------------------
0522 // A set of histograms of particle kinematical variables vs eta
0523 // ------------------------------------------------------------
0524 class HPartVSEta : public Histograms {
0525 public:
0526   HPartVSEta(const TString& name, const double& minMass = 0., const double& maxMass = 100., const double& maxPt = 100.) {
0527     name_ = name;
0528     hPtVSEta_ = new TH2F(name + "_PtVSEta", "transverse momentum vs pseudorapidity", 32, -3.2, 3.2, 200, 0, maxPt);
0529     hMassVSEta_ = new TH2F(name + "_MassVSEta", "mass vs pseudorapidity", 32, -3.2, 3.2, 40, minMass, maxMass);
0530     // TD profile histograms
0531     // ---------------------
0532     hPtVSEta_prof_ = new TProfile(name + "_PtVSEta_prof", "mass vs pseudorapidity", 32, -3.2, 3.2, 0, maxPt);
0533     hMassVSEta_prof_ =
0534         new TProfile(name + "_MassVSEta_prof", "mass vs pseudorapidity", 32, -3.2, 3.2, minMass, maxMass);
0535     hCurvVSEta_prof_ = new TProfile(name + "_CurvVSEta_prof", "curvature vs pseudorapidity", 32, -3.2, 3.2, 0, 1.);
0536   }
0537 
0538   ~HPartVSEta() override {
0539     delete hPtVSEta_;
0540     delete hMassVSEta_;
0541     delete hPtVSEta_prof_;
0542     delete hMassVSEta_prof_;
0543     delete hCurvVSEta_prof_;
0544   }
0545 
0546   void Fill(const reco::Particle::LorentzVector& p4, const double& weight = 1.) override {
0547     Fill(CLHEP::HepLorentzVector(p4.x(), p4.y(), p4.z(), p4.t()), weight);
0548   }
0549 
0550   void Fill(const CLHEP::HepLorentzVector& momentum, const double& weight = 1.) override {
0551     hPtVSEta_->Fill(momentum.eta(), momentum.perp(), weight);
0552     hPtVSEta_prof_->Fill(momentum.eta(), momentum.perp(), weight);
0553 
0554     hMassVSEta_->Fill(momentum.eta(), momentum.m(), weight);
0555     hMassVSEta_prof_->Fill(momentum.eta(), momentum.m(), weight);
0556     hCurvVSEta_prof_->Fill(momentum.eta(), 1 / momentum.perp(), weight);
0557   }
0558 
0559   void Write() override {
0560     hPtVSEta_->Write();
0561     hPtVSEta_prof_->Write();
0562     hCurvVSEta_prof_->Write();
0563     hMassVSEta_->Write();
0564     hMassVSEta_prof_->Write();
0565 
0566     //     std::vector<TGraphErrors*> graphs( (MuScleFitUtils::fitMass(hMassVSEta_)) );
0567     //     for (std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++) {
0568     //       (*graph)->Write();
0569     //     }
0570   }
0571 
0572   void Clear() override {
0573     hPtVSEta_->Clear();
0574     hPtVSEta_prof_->Clear();
0575     hCurvVSEta_prof_->Clear();
0576     hMassVSEta_->Clear();
0577     hMassVSEta_prof_->Clear();
0578   }
0579 
0580 public:
0581   TH2F* hPtVSEta_;
0582   TH2F* hMassVSEta_;
0583   TProfile* hMassVSEta_prof_;
0584   TProfile* hPtVSEta_prof_;
0585   TProfile* hCurvVSEta_prof_;
0586 };
0587 
0588 //---------------------------------------------------------------------------
0589 // A set of histograms of particle kinematical variables vs phi (in eta bins)
0590 // --------------------------------------------------------------------------
0591 class HPartVSPhi : public Histograms {
0592 public:
0593   HPartVSPhi(const TString& name) {
0594     name_ = name;
0595     //     hPtVSPhi_ = new TH2F (name+"_PtVSPhi", "transverse momentum vs phi angle",
0596     //                           12, -3.2, 3.2, 200, 0, 200);
0597     //     hMassVSPhi_ = new TH2F (name+"_MassVSPhi", "mass vs phi angle",
0598     //                             7, -3.2, 3.2, 40, 70, 110);
0599     //     hMassVSPhiF_ = new TH2F (name+"_MassVSPhiF", "mass vs phi F",
0600     //                              7, -3.2, 3.2, 40, 70, 110);
0601     //     hMassVSPhiWp2_ = new TH2F (name+"_MassVSPhiWp2", "mass vs phi Wp2",
0602     //                                7, -3.2, 3.2, 40, 70, 110);
0603     //     hMassVSPhiWp1_ = new TH2F (name+"_MassVSPhiWp1", "mass vs phi Wp1",
0604     //                                7, -3.2, 3.2, 40, 70, 110);
0605     //     hMassVSPhiW0_ = new TH2F (name+"_MassVSPhiW0", "mass vs phi W0",
0606     //                               7, -3.2, 3.2, 40, 70, 110);
0607     //     hMassVSPhiWm1_ = new TH2F (name+"_MassVSPhiWm1", "mass vs phi Wm1",
0608     //                                7, -3.2, 3.2, 40, 70, 110);
0609     //     hMassVSPhiWm2_ = new TH2F (name+"_MassVSPhiWm2", "mass vs phi Wm2",
0610     //                                7, -3.2, 3.2, 40, 70, 110);
0611     //     hMassVSPhiB_ = new TH2F (name+"_MassVSPhiB", "mass vs phi B",
0612     //                              7, -3.2, 3.2, 40, 70, 110);
0613 
0614     // TD profile histograms
0615     hMassVSPhi_prof_ = new TProfile(name + "_MassVSPhi_prof", "mass vs phi angle", 16, -3.2, 3.2, 70, 110);
0616     hPtVSPhi_prof_ = new TProfile(name + "_PtVSPhi_prof", "pt vs phi angle", 16, -3.2, 3.2, 0, 200);
0617   }
0618 
0619   ~HPartVSPhi() override {
0620     delete hPtVSPhi_;
0621     delete hMassVSPhi_;
0622     delete hMassVSPhi_prof_;
0623     delete hPtVSPhi_prof_;
0624 
0625     //    delete hMassVSPhiB_;
0626     //     delete hMassVSPhiWm2_;
0627     //     delete hMassVSPhiWm1_;
0628     //     delete hMassVSPhiW0_;
0629     //     delete hMassVSPhiWp1_;
0630     //     delete hMassVSPhiWp2_;
0631     //     delete hMassVSPhiF_;
0632   }
0633 
0634   void Fill(const reco::Particle::LorentzVector& p4, const double& weight = 1.) override {
0635     Fill(CLHEP::HepLorentzVector(p4.x(), p4.y(), p4.z(), p4.t()), weight);
0636   }
0637 
0638   void Fill(const CLHEP::HepLorentzVector& momentum, const double& weight = 1.) override {
0639     hPtVSPhi_->Fill(momentum.phi(), momentum.perp(), weight);
0640     hMassVSPhi_->Fill(momentum.phi(), momentum.m(), weight);
0641     hMassVSPhi_prof_->Fill(momentum.phi(), momentum.m(), weight);
0642     hPtVSPhi_prof_->Fill(momentum.phi(), momentum.perp(), weight);
0643 
0644     //     if (momentum.eta()<-1.2)                             hMassVSPhiB_->Fill(momentum.phi(),momentum.m(), weight);
0645     //     else if (momentum.eta()<-0.8 && momentum.eta()>-1.2) hMassVSPhiWm2_->Fill(momentum.phi(),momentum.m(), weight);
0646     //     else if (momentum.eta()<-0.3 && momentum.eta()>-0.8) hMassVSPhiWm1_->Fill(momentum.phi(),momentum.m(), weight);
0647     //     else if ((fabs(momentum.eta())) < 0.3)               hMassVSPhiW0_->Fill(momentum.phi(),momentum.m(), weight);
0648     //     else if (momentum.eta()>0.3 && momentum.eta()<0.8)   hMassVSPhiWp1_->Fill(momentum.phi(),momentum.m(), weight);
0649     //     else if (momentum.eta()>0.8 && momentum.eta()<1.2)   hMassVSPhiWp2_->Fill(momentum.phi(),momentum.m(), weight);
0650     //     else if (momentum.eta()>1.2)                         hMassVSPhiF_->Fill(momentum.phi(),momentum.m(), weight);
0651   }
0652 
0653   void Write() override {
0654     hPtVSPhi_->Write();
0655     hMassVSPhi_->Write();
0656     hMassVSPhi_prof_->Write();
0657     hPtVSPhi_prof_->Write();
0658 
0659     //    hMassVSPhiB_->Write();
0660     //     hMassVSPhiWm2_->Write();
0661     //     hMassVSPhiWm1_->Write();
0662     //     hMassVSPhiW0_->Write();
0663     //     hMassVSPhiWp1_->Write();
0664     //     hMassVSPhiWp2_->Write();
0665     //     hMassVSPhiF_->Write();
0666 
0667     //     std::vector<TGraphErrors*> graphs ((MuScleFitUtils::fitMass(hMassVSPhi_)));
0668     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++){
0669     //       (*graph)->Write();
0670     //     }
0671     //     std::vector<TGraphErrors*> graphsB ((MuScleFitUtils::fitMass(hMassVSPhiB_)));
0672     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsB.begin(); graph != graphsB.end(); graph++){
0673     //       (*graph)->Write();
0674     //     }
0675     //     std::vector<TGraphErrors*> graphsWm2 ((MuScleFitUtils::fitMass(hMassVSPhiWm2_)));
0676     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsWm2.begin(); graph != graphsWm2.end(); graph++){
0677     //       (*graph)->Write();
0678     //     }
0679     //     std::vector<TGraphErrors*> graphsWm1 ((MuScleFitUtils::fitMass(hMassVSPhiWm1_)));
0680     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsWm1.begin(); graph != graphsWm1.end(); graph++){
0681     //       (*graph)->Write();
0682     //     }
0683     //     std::vector<TGraphErrors*> graphsW0 ((MuScleFitUtils::fitMass(hMassVSPhiW0_)));
0684     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsW0.begin(); graph != graphsW0.end(); graph++){
0685     //       (*graph)->Write();
0686     //     }
0687     //     std::vector<TGraphErrors*> graphsWp1 ((MuScleFitUtils::fitMass(hMassVSPhiWp1_)));
0688     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsWp1.begin(); graph != graphsWp1.end(); graph++){
0689     //       (*graph)->Write();
0690     //     }
0691     //     std::vector<TGraphErrors*> graphsWp2 ((MuScleFitUtils::fitMass(hMassVSPhiWp2_)));
0692     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsWp2.begin(); graph != graphsWp2.end(); graph++){
0693     //       (*graph)->Write();
0694     //     }
0695     //     std::vector<TGraphErrors*> graphsF ((MuScleFitUtils::fitMass(hMassVSPhiF_)));
0696     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphsF.begin(); graph != graphsF.end(); graph++){
0697     //       (*graph)->Write();
0698     //     }
0699   }
0700 
0701   void Clear() override {
0702     hPtVSPhi_->Clear();
0703     hMassVSPhi_->Clear();
0704     hPtVSPhi_prof_->Clear();
0705     hMassVSPhi_prof_->Clear();
0706 
0707     //     hMassVSPhiB_->Clear();
0708     //     hMassVSPhiWm2_->Clear();
0709     //     hMassVSPhiWm1_->Clear();
0710     //     hMassVSPhiW0_->Clear();
0711     //     hMassVSPhiWp1_->Clear();
0712     //     hMassVSPhiWp2_->Clear();
0713     //     hMassVSPhiF_->Clear();
0714   }
0715 
0716 public:
0717   TH2F* hPtVSPhi_;
0718   TH2F* hMassVSPhi_;
0719   TProfile* hMassVSPhi_prof_;
0720   TProfile* hPtVSPhi_prof_;
0721 
0722   //   TH2F *hMassVSPhiB_;
0723   //   TH2F *hMassVSPhiWm2_;
0724   //   TH2F *hMassVSPhiWm1_;
0725   //   TH2F *hMassVSPhiW0_;
0726   //   TH2F *hMassVSPhiWp1_;
0727   //   TH2F *hMassVSPhiWp2_;
0728   //   TH2F *hMassVSPhiF_;
0729 };
0730 
0731 //---------------------------------------------------------------------------------------
0732 // A set of histograms of particle VS pt
0733 class HPartVSPt : public Histograms {
0734 public:
0735   HPartVSPt(const TString& name) {
0736     name_ = name;
0737     hMassVSPt_ = new TH2F(name + "_MassVSPt", "mass vs transverse momentum", 12, -6, 6, 40, 70, 110);
0738     // TD profile histograms
0739     hMassVSPt_prof_ = new TProfile(name + "_MassVSPt_prof", "mass vs transverse momentum", 12, -3, 3, 86, 116);
0740   }
0741 
0742   ~HPartVSPt() override {
0743     delete hMassVSPt_;
0744     delete hMassVSPt_prof_;
0745   }
0746 
0747   void Fill(const reco::Particle::LorentzVector& p4, const double& weight = 1.) override {
0748     Fill(CLHEP::HepLorentzVector(p4.x(), p4.y(), p4.z(), p4.t()), weight);
0749   }
0750 
0751   void Fill(const CLHEP::HepLorentzVector& momentum, const double& weight = 1.) override {
0752     hMassVSPt_->Fill(momentum.eta(), momentum.m(), weight);
0753     hMassVSPt_prof_->Fill(momentum.eta(), momentum.m(), weight);
0754   }
0755 
0756   void Write() override {
0757     hMassVSPt_->Write();
0758     hMassVSPt_prof_->Write();
0759 
0760     //     std::vector<TGraphErrors*> graphs( (MuScleFitUtils::fitMass(hMassVSPt_)) );
0761     //     for(std::vector<TGraphErrors*>::const_iterator graph = graphs.begin(); graph != graphs.end(); graph++){
0762     //       (*graph)->Write();
0763     //     }
0764   }
0765 
0766   void Clear() override {
0767     hMassVSPt_->Clear();
0768     hMassVSPt_prof_->Clear();
0769   }
0770 
0771 public:
0772   TH2F* hMassVSPt_;
0773   TProfile* hMassVSPt_prof_;
0774 };
0775 
0776 // ---------------------------------------------------
0777 // A set of histograms of resonance mass versus muon variables
0778 class HMassVSPart : public Histograms {
0779 public:
0780   HMassVSPart(const TString& name, const double& minMass = 0., const double& maxMass = 150., const double maxPt = 100.) {
0781     name_ = name;
0782 
0783     // Kinematical variables
0784     // ---------------------
0785     hMassVSPt_ = new TH2F(
0786         name + "_MassVSPt", "resonance mass vs muon transverse momentum", 200, 0., maxPt, 6000, minMass, maxMass);
0787 
0788     hMassVSEta_ =
0789         new TH2F(name + "_MassVSEta", "resonance mass vs muon pseudorapidity", 64, -6.4, 6.4, 6000, minMass, maxMass);
0790     hMassVSEtaPlus_ = new TH2F(
0791         name + "_MassVSEtaPlus", "resonance mass vs muon+ pseudorapidity", 64, -6.4, 6.4, 6000, minMass, maxMass);
0792     hMassVSEtaMinus_ = new TH2F(
0793         name + "_MassVSEtaMinus", "resonance mass vs muon- pseudorapidity", 64, -6.4, 6.4, 6000, minMass, maxMass);
0794 
0795     hMassVSPhiPlus_ =
0796         new TH2F(name + "_MassVSPhiPlus", "resonance mass vs muon+ phi angle", 64, -3.2, 3.2, 6000, minMass, maxMass);
0797     hMassVSPhiMinus_ =
0798         new TH2F(name + "_MassVSPhiMinus", "resonance mass vs muon- phi angle", 64, -3.2, 3.2, 6000, minMass, maxMass);
0799 
0800     // J/Psi mass -----
0801     //     hMassVSEtaPhiPlus_      = new TH3F( name+"_MassVSEtaPhiPlus", "resonance mass vs muon+ phi/pseudorapidity",6, -3.2, 3.2, 20, -2.5, 2.5, 6000, minMass, maxMass );
0802     //     hMassVSEtaPhiMinus_      = new TH3F( name+"_MassVSEtaPhiMinus", "resonance mass vs muon- phi/pseudorapidity", 6, -3.2, 3.2, 20, -2.5, 2.5, 6000, minMass, maxMass );
0803 
0804     //Z mass -----------
0805     hMassVSEtaPhiPlus_ = new TH3F(name + "_MassVSEtaPhiPlus",
0806                                   "resonance mass vs muon+ phi/pseudorapidity",
0807                                   16,
0808                                   -3.2,
0809                                   3.2,
0810                                   20,
0811                                   -2.4,
0812                                   2.4,
0813                                   300,
0814                                   minMass,
0815                                   maxMass);
0816     hMassVSEtaPhiMinus_ = new TH3F(name + "_MassVSEtaPhiMinus",
0817                                    "resonance mass vs muon- phi/pseudorapidity",
0818                                    16,
0819                                    -3.2,
0820                                    3.2,
0821                                    20,
0822                                    -2.4,
0823                                    2.4,
0824                                    300,
0825                                    minMass,
0826                                    maxMass);
0827 
0828     hMassVSCosThetaCS_ = new TH2F(
0829         name + "_MassVSCosThetaCS", "resonance mass vs cos(theta) (CS frame)", 40, -1., 1., 6000, minMass, maxMass);
0830     hMassVSPhiCS_ =
0831         new TH2F(name + "_MassVSPhiCS", "resonance mass vs phi (CS frame)", 64, -3.2, 3.2, 6000, minMass, maxMass);
0832 
0833     hMassVSPhiPlusPhiMinus_ = new TH3F(name + "_MassVSPhiPlusPhiMinus",
0834                                        "resonance mass vs muon+ phi/muon- phi",
0835                                        16,
0836                                        -3.2,
0837                                        3.2,
0838                                        16,
0839                                        -3.2,
0840                                        3.2,
0841                                        6000,
0842                                        minMass,
0843                                        maxMass);
0844     hMassVSEtaPlusEtaMinus_ = new TH3F(name + "_MassVSEtaPlusEtaMinus",
0845                                        "resonance mass vs muon+ eta/muon- eta",
0846                                        16,
0847                                        -3.2,
0848                                        3.2,
0849                                        16,
0850                                        -3.2,
0851                                        3.2,
0852                                        6000,
0853                                        minMass,
0854                                        maxMass);
0855 
0856     hMassVSPhiPlusMinusDiff_ = new TH2F(name + "_MassVSPhiPlusMinusDiff",
0857                                         "resonance mass vs delta phi between mu+/mu-",
0858                                         64,
0859                                         -6.4,
0860                                         6.4,
0861                                         6000,
0862                                         minMass,
0863                                         maxMass);
0864     hMassVSEtaPlusMinusDiff_ = new TH2F(name + "_MassVSEtaPlusMinusDiff",
0865                                         "resonance mass vs delta pseudorapidity between mu+/mu-",
0866                                         32,
0867                                         -4.4,
0868                                         4.4,
0869                                         6000,
0870                                         minMass,
0871                                         maxMass);
0872     hMassVSCosThetaCS_prof =
0873         new TProfile(name + "_MassVScosTheta_prof", "resonance mass vs cosTheta", 40, -1., 1., 85., 95.);
0874 
0875     //hMassVSPt_prof       = new TProfile (name+"_MassVSPt_prof", "resonance mass vs muon transverse momentum", 100, 0., 200., minMass, maxMass);
0876     //hMassVSEta_prof      = new TProfile (name+"_MassVSEta_prof", "resonance mass vs muon pseudorapidity", 30, -6., 6., minMass, maxMass);
0877     //hMassVSPhiPlus_prof  = new TProfile (name+"_MassVSPhiPlus_prof", "resonance mass vs muon+ phi angle", 32, -3.2, 3.2, minMass, maxMass);
0878     //hMassVSPhiMinus_prof = new TProfile (name+"_MassVSPhiMinus_prof", "resonance mass vs muon- phi angle", 32, -3.2, 3.2, minMass, maxMass);
0879   }
0880 
0881   HMassVSPart(const TString& name, TFile* file) {
0882     name_ = name;
0883     hMassVSPt_ = (TH2F*)file->Get(name + "_MassVSPt");
0884     hMassVSEta_ = (TH2F*)file->Get(name + "_MassVSEta");
0885 
0886     hMassVSEtaPhiPlus_ = (TH3F*)file->Get(name + "_MassVSEtaPlus");
0887     hMassVSEtaPhiMinus_ = (TH3F*)file->Get(name + "_MassVSEtaMinus");
0888     hMassVSEtaPlus_ = (TH2F*)file->Get(name + "_MassVSEtaPlus");
0889     hMassVSEtaMinus_ = (TH2F*)file->Get(name + "_MassVSEtaMinus");
0890 
0891     hMassVSPhiPlusMinusDiff_ = (TH2F*)file->Get(name + "_MassVSPhiPlusMinusDiff");
0892     hMassVSEtaPlusMinusDiff_ = (TH2F*)file->Get(name + "_MassVSEtaPlusMinusDiff");
0893 
0894     hMassVSPhiPlus_ = (TH2F*)file->Get(name + "_MassVSPhiPlus");
0895     hMassVSPhiMinus_ = (TH2F*)file->Get(name + "_MassVSPhiMinus");
0896 
0897     hMassVSCosThetaCS_prof = (TProfile*)file->Get(name + "_MassVScosTheta_prof");
0898     //hMassVSPt_prof       = (TProfile *) file->Get(name+"_MassVSPt_prof");
0899     //hMassVSEta_prof      = (TProfile *) file->Get(name+"_MassVSEta_prof");
0900     //hMassVSPhiPlus_prof  = (TProfile *) file->Get(name+"_MassVSPhiPlus_prof");
0901     //hMassVSPhiMinus_prof = (TProfile *) file->Get(name+"_MassVSPhiMinus_prof");
0902   }
0903 
0904   ~HMassVSPart() override {
0905     delete hMassVSPt_;
0906     delete hMassVSEta_;
0907     delete hMassVSPhiPlus_;
0908     delete hMassVSPhiMinus_;
0909     delete hMassVSEtaPhiPlus_;
0910     delete hMassVSEtaPhiMinus_;
0911     delete hMassVSEtaPlus_;
0912     delete hMassVSEtaMinus_;
0913     delete hMassVSPhiPlusPhiMinus_;
0914     delete hMassVSEtaPlusEtaMinus_;
0915     delete hMassVSCosThetaCS_;
0916     delete hMassVSPhiCS_;
0917     delete hMassVSPhiPlusMinusDiff_;
0918     delete hMassVSEtaPlusMinusDiff_;
0919     delete hMassVSCosThetaCS_prof;
0920   }
0921 
0922   void Fill(const reco::Particle::LorentzVector& p41,
0923             const reco::Particle::LorentzVector& p42,
0924             const int charge,
0925             const double& weight = 1.) override {
0926     Fill(CLHEP::HepLorentzVector(p41.x(), p41.y(), p41.z(), p41.t()),
0927          CLHEP::HepLorentzVector(p42.x(), p42.y(), p42.z(), p42.t()),
0928          charge,
0929          weight);
0930   }
0931 
0932   void Fill(const reco::Particle::LorentzVector& p41,
0933             const reco::Particle::LorentzVector& p42,
0934             const reco::Particle::LorentzVector& p4Res,
0935             const double& weight = 1.) override {
0936     Fill(CLHEP::HepLorentzVector(p41.x(), p41.y(), p41.z(), p41.t()),
0937          CLHEP::HepLorentzVector(p42.x(), p42.y(), p42.z(), p42.t()),
0938          CLHEP::HepLorentzVector(p4Res.x(), p4Res.y(), p4Res.z(), p4Res.t()),
0939          weight);
0940   }
0941 
0942   /// Used to fill 2D histograms for comparison of opposite charge muons quantities
0943   void Fill(const CLHEP::HepLorentzVector& momentum1,
0944             const CLHEP::HepLorentzVector& momentum2,
0945             const CLHEP::HepLorentzVector& momentumRes,
0946             const double& weight = 1.) override {
0947     /************************************************************************
0948      *
0949      * Observable: cos(theta) = 2 Q^-1 (Q^2+Qt^2)^-(1/2) (mu^+ mubar^- - mu^- mubar^+)
0950      * (computed in Collins-Soper frame)
0951      *
0952      ************************************************************************/
0953 
0954     double costhetaCS, phiCS;
0955 
0956     const CLHEP::HepLorentzVector& mu = momentum1;
0957     const CLHEP::HepLorentzVector& mubar = momentum2;
0958     CLHEP::HepLorentzVector Q(mu + mubar);
0959     double muplus = 1.0 / sqrt(2.0) * (mu.e() + mu.z());
0960     double muminus = 1.0 / sqrt(2.0) * (mu.e() - mu.z());
0961     double mubarplus = 1.0 / sqrt(2.0) * (mubar.e() + mubar.z());
0962     double mubarminus = 1.0 / sqrt(2.0) * (mubar.e() - mubar.z());
0963     //double costheta = 2.0 / Q.Mag() / sqrt(pow(Q.Mag(), 2) + pow(Q.Pt(), 2)) * (muplus * mubarminus - muminus * mubarplus);
0964     costhetaCS = 2.0 / Q.mag() / sqrt(pow(Q.mag(), 2) + pow(Q.perp(), 2)) * (muplus * mubarminus - muminus * mubarplus);
0965     if (momentumRes.rapidity() < 0)
0966       costhetaCS = -costhetaCS;
0967 
0968     /************************************************************************
0969    *
0970    * 3) tanphi = (Q^2 + Qt^2)^1/2 / Q (Dt dot R unit) /(Dt dot Qt unit)
0971    *
0972    ************************************************************************/
0973 
0974     // unit vector on R direction
0975     CLHEP::HepLorentzVector Pbeam(0., 0., 3500., 3500.);
0976     CLHEP::Hep3Vector R = Pbeam.vect().cross(Q.vect());
0977     CLHEP::Hep3Vector Runit = R.unit();
0978 
0979     // unit vector on Qt
0980     CLHEP::Hep3Vector Qt = Q.vect();
0981     Qt.setZ(0);
0982     CLHEP::Hep3Vector Qtunit = Qt.unit();
0983 
0984     CLHEP::HepLorentzVector D(mu - mubar);
0985     CLHEP::Hep3Vector Dt = D.vect();
0986     Dt.setZ(0);
0987     double tanphi = sqrt(pow(Q.mag(), 2) + pow(Q.perp(), 2)) / Q.mag() * Dt.dot(Runit) / Dt.dot(Qtunit);
0988     if (momentumRes.rapidity() < 0)
0989       tanphi = -tanphi;
0990     phiCS = atan(tanphi);
0991 
0992     hMassVSPhiCS_->Fill(phiCS, momentumRes.m(), weight);
0993     hMassVSCosThetaCS_->Fill(costhetaCS, momentumRes.m(), weight);
0994     hMassVSCosThetaCS_prof->Fill(costhetaCS, momentumRes.m());
0995     /*************************************************************************
0996      *************************************************************************/
0997 
0998     hMassVSPhiPlusPhiMinus_->Fill(momentum1.phi(), momentum2.phi(), momentumRes.m(), weight);
0999     hMassVSEtaPlusEtaMinus_->Fill(momentum1.eta(), momentum2.eta(), momentumRes.m(), weight);
1000 
1001     hMassVSPhiPlusMinusDiff_->Fill((momentum1.phi() - momentum2.phi()), momentumRes.m(), weight);
1002     hMassVSEtaPlusMinusDiff_->Fill((momentum1.eta() - momentum2.eta()), momentumRes.m(), weight);
1003   }
1004 
1005   void Fill(const CLHEP::HepLorentzVector& momentum1,
1006             const CLHEP::HepLorentzVector& momentum2,
1007             const int charge,
1008             const double& weight = 1.) override {
1009     hMassVSPt_->Fill(momentum1.perp(), momentum2.m(), weight);
1010     //hMassVSPt_prof_->Fill(momentum1.perp(),momentum2.m());
1011 
1012     hMassVSEta_->Fill(momentum1.eta(), momentum2.m(), weight);
1013     //hMassVSEta_prof_->Fill(momentum1.eta(),momentum2.m());
1014 
1015     if (charge > 0) {
1016       hMassVSPhiPlus_->Fill(momentum1.phi(), momentum2.m(), weight);
1017       hMassVSEtaPlus_->Fill(momentum1.eta(), momentum2.m(), weight);
1018       hMassVSEtaPhiPlus_->Fill(momentum1.phi(), momentum1.eta(), momentum2.m(), weight);
1019     } else if (charge < 0) {
1020       hMassVSPhiMinus_->Fill(momentum1.phi(), momentum2.m(), weight);
1021       hMassVSEtaMinus_->Fill(momentum1.eta(), momentum2.m(), weight);
1022       hMassVSEtaPhiMinus_->Fill(momentum1.phi(), momentum1.eta(), momentum2.m(), weight);
1023     } else {
1024       LogDebug("Histograms") << "HMassVSPart: wrong charge value = " << charge << std::endl;
1025       abort();
1026     }
1027   }
1028 
1029   void Write() override {
1030     hMassVSPt_->Write();
1031     hMassVSEta_->Write();
1032     hMassVSPhiPlus_->Write();
1033     hMassVSPhiMinus_->Write();
1034 
1035     hMassVSEtaPhiPlus_->Write();
1036     hMassVSEtaPhiMinus_->Write();
1037     hMassVSEtaPlus_->Write();
1038     hMassVSEtaMinus_->Write();
1039 
1040     hMassVSPhiPlusPhiMinus_->Write();
1041     hMassVSEtaPlusEtaMinus_->Write();
1042     hMassVSCosThetaCS_->Write();
1043     hMassVSPhiCS_->Write();
1044 
1045     hMassVSPhiPlusMinusDiff_->Write();
1046     hMassVSEtaPlusMinusDiff_->Write();
1047     hMassVSCosThetaCS_prof->Write();
1048 
1049     //hMassVSPt_prof_->Write();
1050     //hMassVSEta_prof_->Write();
1051     //hMassVSPhiPlus_prof_->Write();
1052     //hMassVSPhiMinus_prof_->Write();
1053   }
1054 
1055   void Clear() override {
1056     hMassVSPt_->Clear();
1057     hMassVSEta_->Clear();
1058     hMassVSPhiPlus_->Clear();
1059     hMassVSPhiMinus_->Clear();
1060 
1061     hMassVSEtaPhiPlus_->Clear();
1062     hMassVSEtaPhiMinus_->Clear();
1063     hMassVSEtaPlus_->Clear();
1064     hMassVSEtaMinus_->Clear();
1065 
1066     hMassVSPhiPlusPhiMinus_->Clear();
1067     hMassVSEtaPlusEtaMinus_->Clear();
1068     hMassVSCosThetaCS_->Clear();
1069     hMassVSPhiCS_->Clear();
1070     hMassVSPhiPlusMinusDiff_->Clear();
1071     hMassVSEtaPlusMinusDiff_->Clear();
1072     hMassVSCosThetaCS_prof->Clear();
1073 
1074     //hMassVSPt_prof_->Clear();
1075     //hMassVSEta_prof_->Clear();
1076     //hMassVSPhiPlus_prof_->Clear();
1077     //hMassVSPhiMinus_prof_->Clear();
1078   }
1079 
1080 protected:
1081   TH2F* hMassVSPt_;
1082   TH2F* hMassVSEta_;
1083   TH2F* hMassVSPhiPlus_;
1084   TH2F* hMassVSPhiMinus_;
1085   TH2F* hMassVSCosThetaCS_;
1086   TH2F* hMassVSPhiCS_;
1087 
1088   TH3F* hMassVSEtaPhiPlus_;
1089   TH3F* hMassVSEtaPhiMinus_;
1090   TH2F* hMassVSEtaPlus_;
1091   TH2F* hMassVSEtaMinus_;
1092 
1093   TH2F* hMassVSPhiPlusMinusDiff_;
1094   TH2F* hMassVSEtaPlusMinusDiff_;
1095 
1096   TH3F* hMassVSPhiPlusPhiMinus_;
1097   TH3F* hMassVSEtaPlusEtaMinus_;
1098 
1099   TProfile* hMassVSCosThetaCS_prof;
1100 
1101   //TProfile* hMassVSPt_prof_;
1102   //TProfile* hMassVSEta_prof_;
1103   //TProfile* hMassVSPhiPlus_prof_;
1104   //TProfile* hMassVSPhiMinus_prof_;
1105 };
1106 
1107 // ---------------------------------------------------
1108 // A set of histograms of Z mass versus muon variables
1109 class HMassVSPartProfile : public Histograms {
1110 public:
1111   HMassVSPartProfile(const TString& name,
1112                      const double& minMass = 0.,
1113                      const double& maxMass = 150.,
1114                      const double maxPt = 100.) {
1115     name_ = name;
1116 
1117     // Kinematical variables
1118     // ---------------------
1119     hMassVSPt_ = new TProfile2D(name + "_MassVSPt",
1120                                 "resonance mass vs muon transverse momentum",
1121                                 200,
1122                                 0.,
1123                                 maxPt,
1124                                 6000,
1125                                 minMass,
1126                                 maxMass,
1127                                 0.,
1128                                 100.);
1129     hMassVSEta_ = new TProfile2D(
1130         name + "_MassVSEta", "resonance mass vs muon pseudorapidity", 64, -6.4, 6.4, 6000, minMass, maxMass, 0., 100.);
1131     hMassVSPhiPlus_ = new TProfile2D(
1132         name + "_MassVSPhiPlus", "resonance mass vs muon+ phi angle", 64, -3.2, 3.2, 6000, minMass, maxMass, 0., 100.);
1133     hMassVSPhiMinus_ = new TProfile2D(
1134         name + "_MassVSPhiMinus", "resonance mass vs muon- phi angle", 64, -3.2, 3.2, 6000, minMass, maxMass, 0., 100.);
1135   }
1136 
1137   HMassVSPartProfile(const TString& name, TFile* file) {
1138     name_ = name;
1139     hMassVSPt_ = (TProfile2D*)file->Get(name + "_MassVSPt");
1140     hMassVSEta_ = (TProfile2D*)file->Get(name + "_MassVSEta");
1141     hMassVSPhiPlus_ = (TProfile2D*)file->Get(name + "_MassVSPhiPlus");
1142     hMassVSPhiMinus_ = (TProfile2D*)file->Get(name + "_MassVSPhiMinus");
1143   }
1144 
1145   ~HMassVSPartProfile() override {
1146     delete hMassVSPt_;
1147     delete hMassVSEta_;
1148     delete hMassVSPhiPlus_;
1149     delete hMassVSPhiMinus_;
1150   }
1151 
1152   void Fill(const reco::Particle::LorentzVector& p41,
1153             const reco::Particle::LorentzVector& p42,
1154             const int charge,
1155             const double& weight = 1.) override {
1156     Fill(CLHEP::HepLorentzVector(p41.x(), p41.y(), p41.z(), p41.t()),
1157          CLHEP::HepLorentzVector(p42.x(), p42.y(), p42.z(), p42.t()),
1158          charge,
1159          weight);
1160   }
1161 
1162   void Fill(const CLHEP::HepLorentzVector& momentum1,
1163             const CLHEP::HepLorentzVector& momentum2,
1164             const int charge,
1165             const double& weight = 1.) override {
1166     hMassVSPt_->Fill(momentum1.perp(), momentum2.m(), weight);
1167     hMassVSEta_->Fill(momentum1.eta(), momentum2.m(), weight);
1168     if (charge > 0) {
1169       hMassVSPhiPlus_->Fill(momentum1.phi(), momentum2.m(), weight);
1170     } else if (charge < 0) {
1171       hMassVSPhiMinus_->Fill(momentum1.phi(), momentum2.m(), weight);
1172     } else {
1173       LogDebug("Histograms") << "HMassVSPartProfile: wrong charge value = " << charge << std::endl;
1174       abort();
1175     }
1176   }
1177 
1178   void Write() override {
1179     hMassVSPt_->Write();
1180     hMassVSEta_->Write();
1181     hMassVSPhiPlus_->Write();
1182     hMassVSPhiMinus_->Write();
1183   }
1184 
1185   void Clear() override {
1186     hMassVSPt_->Clear();
1187     hMassVSEta_->Clear();
1188     hMassVSPhiPlus_->Clear();
1189     hMassVSPhiMinus_->Clear();
1190   }
1191 
1192 protected:
1193   TProfile2D* hMassVSPt_;
1194   TProfile2D* hMassVSEta_;
1195   TProfile2D* hMassVSPhiPlus_;
1196   TProfile2D* hMassVSPhiMinus_;
1197 };
1198 
1199 //---------------------------------------------------------------------------------------
1200 /// A set of histograms for resolution
1201 class HResolutionVSPart : public Histograms {
1202 public:
1203   HResolutionVSPart(TFile* outputFile,
1204                     const TString& name,
1205                     const double maxPt = 100,
1206                     const double& yMinEta = 0.,
1207                     const double& yMaxEta = 2.,
1208                     const double& yMinPt = 0.,
1209                     const double& yMaxPt = 2.,
1210                     const bool doProfiles = false)
1211       : Histograms(outputFile, name), doProfiles_(doProfiles) {
1212     // Kinematical variables
1213 
1214     // hReso           = new TH1F (name+"_Reso", "resolution", 4000, -1, 1);
1215     // hResoVSPtEta    = new TH2F (name+"_ResoVSPtEta", "resolution VS pt and #eta", 200, 0, 200, 60, -3, 3);
1216     // hResoVSPt       = new TH2F (name+"_ResoVSPt", "resolution VS pt", 200, 0, 200, 8000, -1, 1);
1217     // //hResoVSPt_prof  = new TProfile (name+"_ResoVSPt_prof", "resolution VS pt", 100, 0, 200, -1, 1);
1218     // hResoVSEta      = new TH2F (name+"_ResoVSEta", "resolution VS eta", 60, -3, 3, 8000, yMinEta, yMaxEta);
1219     // hResoVSTheta    = new TH2F (name+"_ResoVSTheta", "resolution VS theta", 30, 0, TMath::Pi(), 8000, -1, 1);
1220     // //hResoVSEta_prof = new TProfile (name+"_ResoVSEta_prof", "resolution VS eta", 10, -2.5, 2.5, -1, 1);
1221     // hResoVSPhiPlus  = new TH2F (name+"_ResoVSPhiPlus", "resolution VS phi mu+", 14, -3.2, 3.2, 8000, -1, 1);
1222     // hResoVSPhiMinus = new TH2F (name+"_ResoVSPhiMinus", "resolution VS phi mu-", 14, -3.2, 3.2, 8000, -1, 1);
1223     // //hResoVSPhi_prof = new TProfile (name+"_ResoVSPhi_prof", "resolution VS phi", 14, -3.2, 3.2, -1, 1);
1224     // hAbsReso        = new TH1F (name+"_AbsReso", "resolution", 100, 0, 1);
1225     // hAbsResoVSPt    = new TH2F (name+"_AbsResoVSPt", "Abs resolution VS pt", 200, 0, 500, 100, 0, 1);
1226     // hAbsResoVSEta   = new TH2F (name+"_AbsResoVSEta", "Abs resolution VS eta", 60, -3, 3, 100, 0, 1);
1227     // hAbsResoVSPhi   = new TH2F (name+"_AbsResoVSPhi", "Abs resolution VS phi", 14, -3.2, 3.2, 100, 0, 1);
1228 
1229     hReso_ = new TH1F(name + "_Reso", "resolution", 200, -1, 1);
1230     hResoVSPtEta_ = new TH2F(name + "_ResoVSPtEta", "resolution VS pt and #eta", 200, 0, maxPt, 60, -3, 3);
1231     hResoVSPt_ = new TH2F(name + "_ResoVSPt", "resolution VS pt", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1232     hResoVSPt_Bar_ = new TH2F(name + "_ResoVSPt_Bar", "resolution VS pt Barrel", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1233     hResoVSPt_Endc_17_ = new TH2F(
1234         name + "_ResoVSPt_Endc_1.7", "resolution VS pt Endcap (1.4<eta<1.7)", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1235     hResoVSPt_Endc_20_ = new TH2F(
1236         name + "_ResoVSPt_Endc_2.0", "resolution VS pt Endcap (1.7<eta<2.0)", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1237     hResoVSPt_Endc_24_ = new TH2F(
1238         name + "_ResoVSPt_Endc_2.4", "resolution VS pt Endcap (2.0<eta<2.4)", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1239     hResoVSPt_Ovlap_ =
1240         new TH2F(name + "_ResoVSPt_Ovlap", "resolution VS pt overlap", 200, 0, maxPt, 200, yMinPt, yMaxPt);
1241     hResoVSEta_ = new TH2F(name + "_ResoVSEta", "resolution VS eta", 200, -3, 3, 200, yMinEta, yMaxEta);
1242     hResoVSTheta_ = new TH2F(name + "_ResoVSTheta", "resolution VS theta", 30, 0, TMath::Pi(), 200, yMinEta, yMaxEta);
1243     hResoVSPhiPlus_ = new TH2F(name + "_ResoVSPhiPlus", "resolution VS phi mu+", 16, -3.2, 3.2, 200, -1, 1);
1244     hResoVSPhiMinus_ = new TH2F(name + "_ResoVSPhiMinus", "resolution VS phi mu-", 16, -3.2, 3.2, 200, -1, 1);
1245     hAbsReso_ = new TH1F(name + "_AbsReso", "resolution", 100, 0, 1);
1246     hAbsResoVSPt_ = new TH2F(name + "_AbsResoVSPt", "Abs resolution VS pt", 200, 0, maxPt, 100, 0, 1);
1247     hAbsResoVSEta_ = new TH2F(name + "_AbsResoVSEta", "Abs resolution VS eta", 64, -3.2, 3.2, 100, 0, 1);
1248     hAbsResoVSPhi_ = new TH2F(name + "_AbsResoVSPhi", "Abs resolution VS phi", 16, -3.2, 3.2, 100, 0, 1);
1249     if (doProfiles_) {
1250       hResoVSPt_prof_ = new TProfile(name + "_ResoVSPt_prof", "resolution VS pt", 100, 0, maxPt, yMinPt, yMaxPt);
1251       hResoVSPt_Bar_prof_ =
1252           new TProfile(name + "_ResoVSPt_Bar_prof", "resolution VS pt Barrel", 100, 0, maxPt, yMinPt, yMaxPt);
1253       hResoVSPt_Endc_17_prof_ = new TProfile(
1254           name + "_ResoVSPt_Endc_1.7_prof", "resolution VS pt Endcap (1.4<eta<1.7)", 100, 0, maxPt, yMinPt, yMaxPt);
1255       hResoVSPt_Endc_20_prof_ = new TProfile(
1256           name + "_ResoVSPt_Endc_2.0_prof", "resolution VS pt Endcap (1.7<eta<2.0)", 100, 0, maxPt, yMinPt, yMaxPt);
1257       hResoVSPt_Endc_24_prof_ = new TProfile(
1258           name + "_ResoVSPt_Endc_2.4_prof", "resolution VS pt Endcap (2.0<eta<2.4)", 100, 0, maxPt, yMinPt, yMaxPt);
1259       hResoVSPt_Ovlap_prof_ =
1260           new TProfile(name + "_ResoVSPt_Ovlap_prof", "resolution VS pt Overlap", 100, 0, maxPt, yMinPt, yMaxPt);
1261       hResoVSEta_prof_ = new TProfile(name + "_ResoVSEta_prof", "resolution VS eta", 200, -3.0, 3.0, yMinEta, yMaxEta);
1262       hResoVSPhi_prof_ = new TProfile(name + "_ResoVSPhi_prof", "resolution VS phi", 16, -3.2, 3.2, -1, 1);
1263     }
1264   }
1265 
1266   HResolutionVSPart(const TString& name, TFile* file) {
1267     name_ = name;
1268     hReso_ = (TH1F*)file->Get(name + "_Reso");
1269     hResoVSPtEta_ = (TH2F*)file->Get(name + "_ResoVSPtEta");
1270     hResoVSPt_ = (TH2F*)file->Get(name + "_ResoVSPt");
1271     hResoVSPt_Bar_ = (TH2F*)file->Get(name + "_ResoVSPt");
1272     hResoVSPt_Endc_17_ = (TH2F*)file->Get(name + "_ResoVSPt");
1273     hResoVSPt_Endc_20_ = (TH2F*)file->Get(name + "_ResoVSPt");
1274     hResoVSPt_Endc_24_ = (TH2F*)file->Get(name + "_ResoVSPt");
1275     hResoVSPt_Ovlap_ = (TH2F*)file->Get(name + "_ResoVSPt");
1276     hResoVSEta_ = (TH2F*)file->Get(name + "_ResoVSEta");
1277     hResoVSTheta_ = (TH2F*)file->Get(name + "_ResoVSTheta");
1278     hResoVSPhiPlus_ = (TH2F*)file->Get(name + "_ResoVSPhiPlus");
1279     hResoVSPhiMinus_ = (TH2F*)file->Get(name + "_ResoVSPhiMinus");
1280     hAbsReso_ = (TH1F*)file->Get(name + "_AbsReso");
1281     hAbsResoVSPt_ = (TH2F*)file->Get(name + "_AbsResoVSPt");
1282     hAbsResoVSEta_ = (TH2F*)file->Get(name + "_AbsResoVSEta");
1283     hAbsResoVSPhi_ = (TH2F*)file->Get(name + "_AbsResoVSPhi");
1284     if (doProfiles_) {
1285       hResoVSPt_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_prof");
1286       hResoVSPt_Bar_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_prof");
1287       hResoVSPt_Endc_17_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_1.7_prof");
1288       hResoVSPt_Endc_20_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_2.0_prof");
1289       hResoVSPt_Endc_24_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_2.4_prof");
1290       hResoVSPt_Ovlap_prof_ = (TProfile*)file->Get(name + "_ResoVSPt_prof");
1291       hResoVSEta_prof_ = (TProfile*)file->Get(name + "_ResoVSEta_prof");
1292       hResoVSPhi_prof_ = (TProfile*)file->Get(name + "_ResoVSPhi_prof");
1293     }
1294   }
1295 
1296   ~HResolutionVSPart() override {
1297     delete hReso_;
1298     delete hResoVSPtEta_;
1299     delete hResoVSPt_;
1300     delete hResoVSPt_Bar_;
1301     delete hResoVSPt_Endc_17_;
1302     delete hResoVSPt_Endc_20_;
1303     delete hResoVSPt_Endc_24_;
1304     delete hResoVSPt_Ovlap_;
1305     delete hResoVSEta_;
1306     delete hResoVSTheta_;
1307     delete hResoVSPhiMinus_;
1308     delete hResoVSPhiPlus_;
1309     delete hAbsReso_;
1310     delete hAbsResoVSPt_;
1311     delete hAbsResoVSEta_;
1312     delete hAbsResoVSPhi_;
1313     if (doProfiles_) {
1314       delete hResoVSPt_prof_;
1315       delete hResoVSPt_Bar_prof_;
1316       delete hResoVSPt_Endc_17_prof_;
1317       delete hResoVSPt_Endc_20_prof_;
1318       delete hResoVSPt_Endc_24_prof_;
1319       delete hResoVSPt_Ovlap_prof_;
1320       delete hResoVSEta_prof_;
1321       delete hResoVSPhi_prof_;
1322     }
1323   }
1324 
1325   void Fill(const reco::Particle::LorentzVector& p4, const double& resValue, const int charge) override {
1326     double pt = p4.Pt();
1327     double eta = p4.Eta();
1328     hReso_->Fill(resValue);
1329     hResoVSPtEta_->Fill(pt, eta, resValue);
1330     hResoVSPt_->Fill(pt, resValue);
1331     if (fabs(eta) <= 0.9)
1332       hResoVSPt_Bar_->Fill(pt, resValue);
1333     else if (fabs(eta) > 0.9 && fabs(eta) <= 1.4)
1334       hResoVSPt_Ovlap_->Fill(pt, resValue);
1335     else if (fabs(eta) > 1.4 && fabs(eta) <= 1.7)
1336       hResoVSPt_Endc_17_->Fill(pt, resValue);
1337     else if (fabs(eta) > 1.7 && fabs(eta) <= 2.0)
1338       hResoVSPt_Endc_20_->Fill(pt, resValue);
1339     else
1340       hResoVSPt_Endc_24_->Fill(pt, resValue);
1341 
1342     hResoVSEta_->Fill(eta, resValue);
1343     hResoVSTheta_->Fill(p4.Theta(), resValue);
1344     if (charge > 0)
1345       hResoVSPhiPlus_->Fill(p4.Phi(), resValue);
1346     else if (charge < 0)
1347       hResoVSPhiMinus_->Fill(p4.Phi(), resValue);
1348     hAbsReso_->Fill(fabs(resValue));
1349     hAbsResoVSPt_->Fill(pt, fabs(resValue));
1350     hAbsResoVSEta_->Fill(eta, fabs(resValue));
1351     hAbsResoVSPhi_->Fill(p4.Phi(), fabs(resValue));
1352     if (doProfiles_) {
1353       hResoVSPt_prof_->Fill(p4.Pt(), resValue);
1354       if (fabs(eta) <= 0.9)
1355         hResoVSPt_Bar_prof_->Fill(p4.Pt(), resValue);
1356       else if (fabs(eta) > 0.9 && fabs(eta) <= 1.4)
1357         hResoVSPt_Ovlap_prof_->Fill(pt, resValue);
1358       else if (fabs(eta) > 1.4 && fabs(eta) <= 1.7)
1359         hResoVSPt_Endc_17_prof_->Fill(pt, resValue);
1360       else if (fabs(eta) > 1.7 && fabs(eta) <= 2.0)
1361         hResoVSPt_Endc_20_prof_->Fill(pt, resValue);
1362       else
1363         hResoVSPt_Endc_24_prof_->Fill(pt, resValue);
1364 
1365       hResoVSEta_prof_->Fill(p4.Eta(), resValue);
1366       hResoVSPhi_prof_->Fill(p4.Phi(), resValue);
1367     }
1368   }
1369 
1370   void Write() override {
1371     if (histoDir_ != nullptr)
1372       histoDir_->cd();
1373 
1374     hReso_->Write();
1375     hResoVSPtEta_->Write();
1376     hResoVSPt_->Write();
1377     hResoVSPt_Bar_->Write();
1378     hResoVSPt_Endc_17_->Write();
1379     hResoVSPt_Endc_20_->Write();
1380     hResoVSPt_Endc_24_->Write();
1381     hResoVSPt_Ovlap_->Write();
1382     hResoVSEta_->Write();
1383     hResoVSTheta_->Write();
1384     hResoVSPhiMinus_->Write();
1385     hResoVSPhiPlus_->Write();
1386     hAbsReso_->Write();
1387     hAbsResoVSPt_->Write();
1388     hAbsResoVSEta_->Write();
1389     hAbsResoVSPhi_->Write();
1390     if (doProfiles_) {
1391       hResoVSPt_prof_->Write();
1392       hResoVSPt_Bar_prof_->Write();
1393       hResoVSPt_Endc_17_prof_->Write();
1394       hResoVSPt_Endc_20_prof_->Write();
1395       hResoVSPt_Endc_24_prof_->Write();
1396       hResoVSPt_Ovlap_prof_->Write();
1397       hResoVSEta_prof_->Write();
1398       hResoVSPhi_prof_->Write();
1399     }
1400   }
1401 
1402   void Clear() override {
1403     hReso_->Clear();
1404     hResoVSPtEta_->Clear();
1405     hResoVSPt_->Clear();
1406     hResoVSPt_Bar_->Clear();
1407     hResoVSPt_Endc_17_->Clear();
1408     hResoVSPt_Endc_20_->Clear();
1409     hResoVSPt_Endc_24_->Clear();
1410     hResoVSPt_Ovlap_->Clear();
1411     hResoVSEta_->Clear();
1412     hResoVSTheta_->Clear();
1413     hResoVSPhiPlus_->Clear();
1414     hResoVSPhiMinus_->Clear();
1415     hAbsReso_->Clear();
1416     hAbsResoVSPt_->Clear();
1417     hAbsResoVSEta_->Clear();
1418     hAbsResoVSPhi_->Clear();
1419     if (doProfiles_) {
1420       hResoVSPt_prof_->Clear();
1421       hResoVSPt_Bar_prof_->Clear();
1422       hResoVSPt_Endc_17_prof_->Clear();
1423       hResoVSPt_Endc_20_prof_->Clear();
1424       hResoVSPt_Endc_24_prof_->Clear();
1425       hResoVSPt_Ovlap_prof_->Clear();
1426       hResoVSEta_prof_->Clear();
1427       hResoVSPhi_prof_->Clear();
1428     }
1429   }
1430 
1431 public:
1432   TH1F* hReso_;
1433   TH2F* hResoVSPtEta_;
1434   TH2F* hResoVSPt_;
1435   TH2F* hResoVSPt_Bar_;
1436   TH2F* hResoVSPt_Endc_17_;
1437   TH2F* hResoVSPt_Endc_20_;
1438   TH2F* hResoVSPt_Endc_24_;
1439   TH2F* hResoVSPt_Ovlap_;
1440   TProfile* hResoVSPt_prof_;
1441   TProfile* hResoVSPt_Bar_prof_;
1442   TProfile* hResoVSPt_Endc_17_prof_;
1443   TProfile* hResoVSPt_Endc_20_prof_;
1444   TProfile* hResoVSPt_Endc_24_prof_;
1445   TProfile* hResoVSPt_Ovlap_prof_;
1446   TH2F* hResoVSEta_;
1447   TH2F* hResoVSTheta_;
1448   TProfile* hResoVSEta_prof_;
1449   TH2F* hResoVSPhiMinus_;
1450   TH2F* hResoVSPhiPlus_;
1451   TProfile* hResoVSPhi_prof_;
1452   TH1F* hAbsReso_;
1453   TH2F* hAbsResoVSPt_;
1454   TH2F* hAbsResoVSEta_;
1455   TH2F* hAbsResoVSPhi_;
1456   bool doProfiles_;
1457 };
1458 
1459 // -------------------------------------------------------------
1460 // A set of histograms of likelihood value versus muon variables
1461 // -------------------------------------------------------------
1462 class HLikelihoodVSPart : public Histograms {
1463 public:
1464   HLikelihoodVSPart(const TString& name) {
1465     name_ = name;
1466 
1467     // Kinematical variables
1468     // ---------------------
1469     hLikeVSPt_ =
1470         new TH2F(name + "_LikelihoodVSPt", "likelihood vs muon transverse momentum", 100, 0., 100., 100, -100., 100.);
1471     hLikeVSEta_ =
1472         new TH2F(name + "_LikelihoodVSEta", "likelihood vs muon pseudorapidity", 100, -4., 4., 100, -100., 100.);
1473     hLikeVSPhi_ = new TH2F(name + "_LikelihoodVSPhi", "likelihood vs muon phi angle", 100, -3.2, 3.2, 100, -100., 100.);
1474     hLikeVSPt_prof_ = new TProfile(
1475         name + "_LikelihoodVSPt_prof", "likelihood vs muon transverse momentum", 40, 0., 100., -1000., 1000.);
1476     hLikeVSEta_prof_ =
1477         new TProfile(name + "_LikelihoodVSEta_prof", "likelihood vs muon pseudorapidity", 40, -4., 4., -1000., 1000.);
1478     hLikeVSPhi_prof_ =
1479         new TProfile(name + "_LikelihoodVSPhi_prof", "likelihood vs muon phi angle", 40, -3.2, 3.2, -1000., 1000.);
1480   }
1481 
1482   HLikelihoodVSPart(const TString& name, TFile* file) {
1483     name_ = name;
1484     hLikeVSPt_ = (TH2F*)file->Get(name + "_LikelihoodVSPt");
1485     hLikeVSEta_ = (TH2F*)file->Get(name + "_LikelihoodVSEta");
1486     hLikeVSPhi_ = (TH2F*)file->Get(name + "_LikelihoodVSPhi");
1487     hLikeVSPt_prof_ = (TProfile*)file->Get(name + "_LikelihoodVSPt_prof");
1488     hLikeVSEta_prof_ = (TProfile*)file->Get(name + "_LikelihoodVSEta_prof");
1489     hLikeVSPhi_prof_ = (TProfile*)file->Get(name + "_LikelihoodVSPhi_prof");
1490   }
1491 
1492   ~HLikelihoodVSPart() override {
1493     delete hLikeVSPt_;
1494     delete hLikeVSEta_;
1495     delete hLikeVSPhi_;
1496     delete hLikeVSPt_prof_;
1497     delete hLikeVSEta_prof_;
1498     delete hLikeVSPhi_prof_;
1499   }
1500 
1501   void Fill(const reco::Particle::LorentzVector& p4, const double& likeValue) override {
1502     Fill(CLHEP::HepLorentzVector(p4.x(), p4.y(), p4.z(), p4.t()), likeValue);
1503   }
1504 
1505   void Fill(const CLHEP::HepLorentzVector& momentum, const double& likeValue) override {
1506     hLikeVSPt_->Fill(momentum.perp(), likeValue);
1507     hLikeVSEta_->Fill(momentum.eta(), likeValue);
1508     hLikeVSPhi_->Fill(momentum.phi(), likeValue);
1509     hLikeVSPt_prof_->Fill(momentum.perp(), likeValue);
1510     hLikeVSEta_prof_->Fill(momentum.eta(), likeValue);
1511     hLikeVSPhi_prof_->Fill(momentum.phi(), likeValue);
1512   }
1513 
1514   void Write() override {
1515     hLikeVSPt_->Write();
1516     hLikeVSEta_->Write();
1517     hLikeVSPhi_->Write();
1518     hLikeVSPt_prof_->Write();
1519     hLikeVSEta_prof_->Write();
1520     hLikeVSPhi_prof_->Write();
1521   }
1522 
1523   void Clear() override {
1524     hLikeVSPt_->Reset("ICE");
1525     hLikeVSEta_->Reset("ICE");
1526     hLikeVSPhi_->Reset("ICE");
1527     hLikeVSPt_prof_->Reset("ICE");
1528     hLikeVSEta_prof_->Reset("ICE");
1529     hLikeVSPhi_prof_->Reset("ICE");
1530   }
1531 
1532 public:
1533   TH2F* hLikeVSPt_;
1534   TH2F* hLikeVSEta_;
1535   TH2F* hLikeVSPhi_;
1536   TProfile* hLikeVSPt_prof_;
1537   TProfile* hLikeVSEta_prof_;
1538   TProfile* hLikeVSPhi_prof_;
1539 };
1540 
1541 /**
1542  * This histogram class fills a TProfile with the resolution evaluated from the resolution
1543  * functions for single muon quantities. The resolution functions are used by MuScleFit to
1544  * evaluate the mass resolution, which is the value seen by minuit and through it,
1545  * corrections are evaluated. <br>
1546  * In the end we will compare the histograms filled by this class (from the resolution
1547  * function, reflecting the parameters changes done by minuit) with those filled comparing
1548  * recoMuons with genMuons (the real resolutions).
1549  */
1550 class HFunctionResolution : public Histograms {
1551 public:
1552   HFunctionResolution(TFile* outputFile, const TString& name, const double& ptMax = 100, const int totBinsY = 300)
1553       : Histograms(outputFile, name) {
1554     name_ = name;
1555     totBinsX_ = 300;
1556     totBinsY_ = totBinsY;
1557     xMin_ = 0.;
1558     yMin_ = -3.0;
1559     double xMax = ptMax;
1560     double yMax = 3.0;
1561     deltaX_ = xMax - xMin_;
1562     deltaY_ = yMax - yMin_;
1563     hReso_ = new TH1F(name + "_Reso", "resolution", 1000, 0, 1);
1564     hResoVSPtEta_ =
1565         new TH2F(name + "_ResoVSPtEta", "resolution vs pt and #eta", totBinsX_, xMin_, xMax, totBinsY_, yMin_, yMax);
1566     // Create and initialize the resolution arrays
1567     resoVsPtEta_ = new double*[totBinsX_];
1568     resoCount_ = new int*[totBinsX_];
1569     for (int i = 0; i < totBinsX_; ++i) {
1570       resoVsPtEta_[i] = new double[totBinsY_];
1571       resoCount_[i] = new int[totBinsY_];
1572       for (int j = 0; j < totBinsY_; ++j) {
1573         resoVsPtEta_[i][j] = 0;
1574         resoCount_[i][j] = 0;
1575       }
1576     }
1577     hResoVSPt_prof_ = new TProfile(name + "_ResoVSPt_prof", "resolution VS pt", totBinsX_, xMin_, xMax, yMin_, yMax);
1578     hResoVSPt_Bar_prof_ =
1579         new TProfile(name + "_ResoVSPt_Bar_prof", "resolution VS pt Barrel", totBinsX_, xMin_, xMax, yMin_, yMax);
1580     hResoVSPt_Endc_17_prof_ = new TProfile(
1581         name + "_ResoVSPt_Endc_1.7_prof", "resolution VS pt Endcap (1.4<eta<1.7)", totBinsX_, xMin_, xMax, yMin_, yMax);
1582     hResoVSPt_Endc_20_prof_ = new TProfile(
1583         name + "_ResoVSPt_Endc_2.0_prof", "resolution VS pt Endcap (1.7<eta<2.0)", totBinsX_, xMin_, xMax, yMin_, yMax);
1584     hResoVSPt_Endc_24_prof_ = new TProfile(
1585         name + "_ResoVSPt_Endc_2.4_prof", "resolution VS pt Endcap (2.0<eta<2.4)", totBinsX_, xMin_, xMax, yMin_, yMax);
1586     hResoVSPt_Ovlap_prof_ =
1587         new TProfile(name + "_ResoVSPt_Ovlap_prof", "resolution VS pt Overlap", totBinsX_, xMin_, xMax, yMin_, yMax);
1588     hResoVSEta_prof_ = new TProfile(name + "_ResoVSEta_prof", "resolution VS eta", totBinsY_, yMin_, yMax, 0, 1);
1589     //hResoVSTheta_prof_    = new TProfile( name+"_ResoVSTheta_prof", "resolution VS theta", 30, 0, TMath::Pi(), 0, 1);
1590     hResoVSPhiPlus_prof_ = new TProfile(name + "_ResoVSPhiPlus_prof", "resolution VS phi mu+", 16, -3.2, 3.2, 0, 1);
1591     hResoVSPhiMinus_prof_ = new TProfile(name + "_ResoVSPhiMinus_prof", "resolution VS phi mu-", 16, -3.2, 3.2, 0, 1);
1592     hResoVSPhi_prof_ = new TProfile(name + "_ResoVSPhi_prof", "resolution VS phi", 16, -3.2, 3.2, -1, 1);
1593   }
1594   ~HFunctionResolution() override {
1595     delete hReso_;
1596     delete hResoVSPtEta_;
1597     // Free the resolution arrays
1598     for (int i = 0; i < totBinsX_; ++i) {
1599       delete[] resoVsPtEta_[i];
1600       delete[] resoCount_[i];
1601     }
1602     delete[] resoVsPtEta_;
1603     delete[] resoCount_;
1604     // Free the profile histograms
1605     delete hResoVSPt_prof_;
1606     delete hResoVSPt_Bar_prof_;
1607     delete hResoVSPt_Endc_17_prof_;
1608     delete hResoVSPt_Endc_20_prof_;
1609     delete hResoVSPt_Endc_24_prof_;
1610     delete hResoVSPt_Ovlap_prof_;
1611     delete hResoVSEta_prof_;
1612     //delete hResoVSTheta_prof_;
1613     delete hResoVSPhiPlus_prof_;
1614     delete hResoVSPhiMinus_prof_;
1615     delete hResoVSPhi_prof_;
1616   }
1617   void Fill(const reco::Particle::LorentzVector& p4, const double& resValue, const int charge) override {
1618     if (resValue != resValue)
1619       return;
1620     hReso_->Fill(resValue);
1621 
1622     // Fill the arrays with the resolution value and count
1623     int xIndex = getXindex(p4.Pt());
1624     int yIndex = getYindex(p4.Eta());
1625     if (0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_) {
1626       resoVsPtEta_[xIndex][yIndex] += resValue;
1627       // ATTENTION: we count only for positive muons because we are considering di-muon resonances
1628       // and we use this counter to compute the mean in the end. The resoVsPtEta value is filled with each muon,
1629       // but they must be considered independently (analogous to a factor 2) so in the end we would have
1630       // to divide by N/2, that is why we only increase the counter for half the muons.
1631       // if( charge > 0 )
1632       // No more. Changing it here influences also other uses of this class. The macro FunctionTerms.cc
1633       // multiplies the terms by the 2 factor.
1634 
1635       resoCount_[xIndex][yIndex] += 1;
1636 
1637       // hResoVSPtEta->Fill(p4.Pt(), p4.Eta(), resValue);
1638       hResoVSPt_prof_->Fill(p4.Pt(), resValue);
1639       if (fabs(p4.Eta()) <= 0.9)
1640         hResoVSPt_Bar_prof_->Fill(p4.Pt(), resValue);
1641       else if (fabs(p4.Eta()) > 0.9 && fabs(p4.Eta()) <= 1.4)
1642         hResoVSPt_Ovlap_prof_->Fill(p4.Pt(), resValue);
1643       else if (fabs(p4.Eta()) > 1.4 && fabs(p4.Eta()) <= 1.7)
1644         hResoVSPt_Endc_17_prof_->Fill(p4.Pt(), resValue);
1645       else if (fabs(p4.Eta()) > 1.7 && fabs(p4.Eta()) <= 2.0)
1646         hResoVSPt_Endc_20_prof_->Fill(p4.Pt(), resValue);
1647       else
1648         hResoVSPt_Endc_24_prof_->Fill(p4.Pt(), resValue);
1649       hResoVSEta_prof_->Fill(p4.Eta(), resValue);
1650       //hResoVSTheta_prof_->Fill(p4.Theta(),resValue);
1651       if (charge > 0)
1652         hResoVSPhiPlus_prof_->Fill(p4.Phi(), resValue);
1653       else if (charge < 0)
1654         hResoVSPhiMinus_prof_->Fill(p4.Phi(), resValue);
1655       hResoVSPhi_prof_->Fill(p4.Phi(), resValue);
1656     }
1657   }
1658 
1659   void Write() override {
1660     if (histoDir_ != nullptr)
1661       histoDir_->cd();
1662 
1663     hReso_->Write();
1664 
1665     for (int i = 0; i < totBinsX_; ++i) {
1666       for (int j = 0; j < totBinsY_; ++j) {
1667         int N = resoCount_[i][j];
1668         // Fill with the mean value
1669         if (N != 0)
1670           hResoVSPtEta_->SetBinContent(i + 1, j + 1, resoVsPtEta_[i][j] / N);
1671         else
1672           hResoVSPtEta_->SetBinContent(i + 1, j + 1, 0);
1673       }
1674     }
1675     hResoVSPtEta_->Write();
1676 
1677     hResoVSPt_prof_->Write();
1678     hResoVSPt_Bar_prof_->Write();
1679     hResoVSPt_Endc_17_prof_->Write();
1680     hResoVSPt_Endc_20_prof_->Write();
1681     hResoVSPt_Endc_24_prof_->Write();
1682     hResoVSPt_Ovlap_prof_->Write();
1683     hResoVSEta_prof_->Write();
1684     //hResoVSTheta_prof_->Write();
1685     hResoVSPhiMinus_prof_->Write();
1686     hResoVSPhiPlus_prof_->Write();
1687     hResoVSPhi_prof_->Write();
1688 
1689     TCanvas canvas(
1690         TString(hResoVSPtEta_->GetName()) + "_canvas", TString(hResoVSPtEta_->GetTitle()) + " canvas", 1000, 800);
1691     canvas.Divide(2);
1692     canvas.cd(1);
1693     hResoVSPtEta_->Draw("lego");
1694     canvas.cd(2);
1695     hResoVSPtEta_->Draw("surf5");
1696     canvas.Write();
1697     hResoVSPtEta_->Write();
1698 
1699     outputFile_->cd();
1700   }
1701 
1702   void Clear() override {
1703     hReso_->Clear();
1704     hResoVSPtEta_->Clear();
1705     hResoVSPt_prof_->Clear();
1706     hResoVSPt_Bar_prof_->Clear();
1707     hResoVSPt_Endc_17_prof_->Clear();
1708     hResoVSPt_Endc_20_prof_->Clear();
1709     hResoVSPt_Endc_24_prof_->Clear();
1710     hResoVSPt_Ovlap_prof_->Clear();
1711     hResoVSEta_prof_->Clear();
1712     //hResoVSTheta_prof_->Clear();
1713     hResoVSPhiPlus_prof_->Clear();
1714     hResoVSPhiMinus_prof_->Clear();
1715     hResoVSPhi_prof_->Clear();
1716   }
1717 
1718 protected:
1719   int getXindex(const double& x) const { return int((x - xMin_) / deltaX_ * totBinsX_); }
1720   int getYindex(const double& y) const { return int((y - yMin_) / deltaY_ * totBinsY_); }
1721   TH1F* hReso_;
1722   TH2F* hResoVSPtEta_;
1723   double** resoVsPtEta_;
1724   int** resoCount_;
1725   TProfile* hResoVSPt_prof_;
1726   TProfile* hResoVSPt_Bar_prof_;
1727   TProfile* hResoVSPt_Endc_17_prof_;
1728   TProfile* hResoVSPt_Endc_20_prof_;
1729   TProfile* hResoVSPt_Endc_24_prof_;
1730   TProfile* hResoVSPt_Ovlap_prof_;
1731   TProfile* hResoVSEta_prof_;
1732   //TProfile* hResoVSTheta_prof_;
1733   TProfile* hResoVSPhiMinus_prof_;
1734   TProfile* hResoVSPhiPlus_prof_;
1735   TProfile* hResoVSPhi_prof_;
1736   int totBinsX_, totBinsY_;
1737   double xMin_, yMin_;
1738   double deltaX_, deltaY_;
1739 };
1740 
1741 class HFunctionResolutionVarianceCheck : public HFunctionResolution {
1742 public:
1743   HFunctionResolutionVarianceCheck(TFile* outputFile, const TString& name, const double ptMax = 200)
1744       : HFunctionResolution(outputFile, name, ptMax) {
1745     histoVarianceCheck_ = new TH1D**[totBinsX_];
1746     for (int i = 0; i < totBinsX_; ++i) {
1747       histoVarianceCheck_[i] = new TH1D*[totBinsY_];
1748       for (int j = 0; j < totBinsY_; ++j) {
1749         std::stringstream namei;
1750         std::stringstream namej;
1751         namei << i;
1752         namej << j;
1753         histoVarianceCheck_[i][j] = new TH1D(name + "_" + namei.str() + "_" + namej.str(), name, 100, 0., 1.);
1754       }
1755     }
1756   }
1757   ~HFunctionResolutionVarianceCheck() override {
1758     for (int i = 0; i < totBinsX_; ++i) {
1759       for (int j = 0; j < totBinsY_; ++j) {
1760         delete histoVarianceCheck_[i][j];
1761       }
1762       delete[] histoVarianceCheck_;
1763     }
1764     delete[] histoVarianceCheck_;
1765   }
1766   void Fill(const reco::Particle::LorentzVector& p4, const double& resValue, const int charge) override {
1767     if (resValue != resValue)
1768       return;
1769     // Need to convert the (x,y) values to the array indeces
1770     int xIndex = getXindex(p4.Pt());
1771     int yIndex = getYindex(p4.Eta());
1772     // Only fill values if they are in the selected range
1773     if (0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_) {
1774       histoVarianceCheck_[xIndex][yIndex]->Fill(resValue);
1775     }
1776     // Call also the fill of the base class
1777     HFunctionResolution::Fill(p4, resValue, charge);
1778   }
1779   void Write() override {
1780     if (histoDir_ != nullptr)
1781       histoDir_->cd();
1782     for (int xBin = 0; xBin < totBinsX_; ++xBin) {
1783       for (int yBin = 0; yBin < totBinsY_; ++yBin) {
1784         histoVarianceCheck_[xBin][yBin]->Write();
1785       }
1786     }
1787     HFunctionResolution::Write();
1788   }
1789 
1790 protected:
1791   TH1D*** histoVarianceCheck_;
1792 };
1793 
1794 /**
1795  * This histogram class can be used to evaluate the resolution of a variable.</br>
1796  * It has a TProfile, a TH2F and a TH1F. The TProfile is used to compute the rms of
1797  * the distribution which is filled in the TH1F (the resolution histogram) in the
1798  * Write method.</br>
1799  * If a TDirectory is passed to the constructor, the different histograms are
1800  * placed in subdirectories.
1801  */
1802 class HResolution : public TH1F {
1803 public:
1804   HResolution(const TString& name,
1805               const TString& title,
1806               const int totBins,
1807               const double& xMin,
1808               const double& xMax,
1809               const double& yMin,
1810               const double& yMax,
1811               TDirectory* dir = nullptr)
1812       : dir_(dir), dir2D_(nullptr), diffDir_(nullptr) {
1813     if (dir_ != nullptr) {
1814       dir2D_ = (TDirectory*)dir_->Get("2D");
1815       if (dir2D_ == nullptr)
1816         dir2D_ = dir_->mkdir("2D");
1817       diffDir_ = (TDirectory*)dir_->Get("deltaXoverX");
1818       if (diffDir_ == nullptr)
1819         diffDir_ = dir->mkdir("deltaXoverX");
1820     }
1821     diffHisto_ = new TProfile(name + "_prof", title + " profile", totBins, xMin, xMax, yMin, yMax);
1822     histo2D_ = new TH2F(name + "2D", title, totBins, xMin, xMax, 4000, yMin, yMax);
1823     resoHisto_ = new TH1F(name, title, totBins, xMin, xMax);
1824   }
1825   ~HResolution() override {
1826     delete diffHisto_;
1827     delete histo2D_;
1828     delete resoHisto_;
1829   }
1830   Int_t Fill(Double_t x, Double_t y) override {
1831     diffHisto_->Fill(x, y);
1832     histo2D_->Fill(x, y);
1833     return 0;
1834   }
1835   Int_t Write(const char* name = nullptr, Int_t option = 0, Int_t bufsize = 0) override {
1836     // Loop on all the bins and take the rms.
1837     // The TProfile bin error is by default the standard error on the mean, that is
1838     // rms/sqrt(N). If it is created with the "S" option (as we did NOT do), it would
1839     // already be the rms. Thus we take the error and multiply it by the sqrt of the
1840     // bin entries to get the rms.
1841     // bin 0 is the underflow, bin totBins+1 is the overflow.
1842     unsigned int totBins = diffHisto_->GetNbinsX();
1843     // std::cout << "totBins = " << totBins << std::endl;
1844     for (unsigned int iBin = 1; iBin <= totBins; ++iBin) {
1845       // std::cout << "iBin = " << iBin << ", " << diffHisto_->GetBinError(iBin)*sqrt(diffHisto_->GetBinEntries(iBin)) << std::endl;
1846       resoHisto_->SetBinContent(iBin, diffHisto_->GetBinError(iBin) * sqrt(diffHisto_->GetBinEntries(iBin)));
1847     }
1848     if (dir_ != nullptr)
1849       dir_->cd();
1850     resoHisto_->Write();
1851     if (diffDir_ != nullptr)
1852       diffDir_->cd();
1853     diffHisto_->Write();
1854     if (dir2D_ != nullptr)
1855       dir2D_->cd();
1856     histo2D_->Write();
1857 
1858     return 0;
1859   }
1860 
1861 protected:
1862   TDirectory* dir_;
1863   TDirectory* dir2D_;
1864   TDirectory* diffDir_;
1865   TProfile* diffHisto_;
1866   TH2F* histo2D_;
1867   TH1F* resoHisto_;
1868 };
1869 
1870 /**
1871  * This class can be used to compute the covariance between two input variables.
1872  * The Fill method needs the two input variables. </br>
1873  * In the end the covariance method computes the covariance as:
1874  * cov(x,y) = Sum_i(x_i*y_i)/N - x_mean*y_mean. </br>
1875  * Of course passing the same variable for x and y gives the variance of that variable.
1876  */
1877 class Covariance {
1878 public:
1879   Covariance() : productXY_(0), sumX_(0), sumY_(0), N_(0) {}
1880   void fill(const double& x, const double& y) {
1881     productXY_ += x * y;
1882     sumX_ += x;
1883     sumY_ += y;
1884     ++N_;
1885   }
1886   double covariance() {
1887     if (N_ != 0) {
1888       double meanX = sumX_ / N_;
1889       double meanY = sumY_ / N_;
1890       // std::cout << "meanX*meanY = "<<meanX<<"*"<<meanY<< " = " << meanX*meanY << std::endl;
1891       return (productXY_ / N_ - meanX * meanY);
1892     }
1893     return 0.;
1894   }
1895   double getN() { return N_; }
1896 
1897 protected:
1898   double productXY_;
1899   double sumX_;
1900   double sumY_;
1901   int N_;
1902 };
1903 
1904 /**
1905  * This class can be used to compute the covariance of two variables with respect to other two variables
1906  * (to see e.g. how does the covariance of ptVSphi vary with respect to (pt,eta).
1907  */
1908 class HCovarianceVSxy : public Histograms {
1909 public:
1910   HCovarianceVSxy(const TString& name,
1911                   const TString& title,
1912                   const int totBinsX,
1913                   const double& xMin,
1914                   const double& xMax,
1915                   const int totBinsY,
1916                   const double& yMin,
1917                   const double& yMax,
1918                   TDirectory* dir = nullptr,
1919                   bool varianceCheck = false)
1920       : totBinsX_(totBinsX),
1921         totBinsY_(totBinsY),
1922         xMin_(xMin),
1923         deltaX_(xMax - xMin),
1924         yMin_(yMin),
1925         deltaY_(yMax - yMin),
1926         readMode_(false),
1927         varianceCheck_(varianceCheck) {
1928     name_ = name;
1929     histoDir_ = dir;
1930     histoCovariance_ = new TH2D(name + "Covariance", title + " covariance", totBinsX, xMin, xMax, totBinsY, yMin, yMax);
1931 
1932     covariances_ = new Covariance*[totBinsX];
1933     for (int i = 0; i < totBinsX; ++i) {
1934       covariances_[i] = new Covariance[totBinsY];
1935     }
1936     if (varianceCheck_) {
1937       histoVarianceCheck_ = new TH1D**[totBinsX_];
1938       for (int i = 0; i < totBinsX_; ++i) {
1939         histoVarianceCheck_[i] = new TH1D*[totBinsY_];
1940         for (int j = 0; j < totBinsY_; ++j) {
1941           std::stringstream namei;
1942           std::stringstream namej;
1943           namei << i;
1944           namej << j;
1945           histoVarianceCheck_[i][j] = new TH1D(name + "_" + namei.str() + "_" + namej.str(), name, 10000, -1, 1);
1946         }
1947       }
1948     }
1949   }
1950   /// Contructor to read histograms from file
1951   HCovarianceVSxy(TFile* inputFile, const TString& name, const TString& dirName) : readMode_(true) {
1952     histoDir_ = (TDirectory*)(inputFile->Get(dirName.Data()));
1953     if (histoDir_ == nullptr) {
1954       std::cout << "Error: directory not found" << std::endl;
1955       exit(0);
1956     }
1957     histoCovariance_ = (TH2D*)(histoDir_->Get(name));
1958     totBinsX_ = histoCovariance_->GetNbinsX();
1959     xMin_ = histoCovariance_->GetXaxis()->GetBinLowEdge(1);
1960     deltaX_ = histoCovariance_->GetXaxis()->GetBinUpEdge(totBinsX_) - xMin_;
1961     totBinsY_ = histoCovariance_->GetNbinsY();
1962     yMin_ = histoCovariance_->GetYaxis()->GetBinLowEdge(1);
1963     deltaY_ = histoCovariance_->GetYaxis()->GetBinUpEdge(totBinsY_) - yMin_;
1964   }
1965 
1966   ~HCovarianceVSxy() override {
1967     delete histoCovariance_;
1968     // Free covariances
1969     for (int i = 0; i < totBinsX_; ++i) {
1970       delete[] covariances_[i];
1971     }
1972     delete[] covariances_;
1973     // Free variance check histograms
1974     if (varianceCheck_) {
1975       for (int i = 0; i < totBinsX_; ++i) {
1976         for (int j = 0; j < totBinsY_; ++j) {
1977           delete histoVarianceCheck_[i][j];
1978         }
1979         delete[] histoVarianceCheck_[i];
1980       }
1981       delete[] histoVarianceCheck_;
1982     }
1983   }
1984 
1985   /**
1986    * x and y should be the variables VS which we are computing the covariance (pt and eta)
1987    * a and b should be the variables OF which we are computing the covariance </br>
1988    */
1989   void Fill(const double& x, const double& y, const double& a, const double& b) override {
1990     // Need to convert the (x,y) values to the array indeces
1991     int xIndex = getXindex(x);
1992     int yIndex = getYindex(y);
1993     // Only fill values if they are in the selected range
1994     if (0 <= xIndex && xIndex < totBinsX_ && 0 <= yIndex && yIndex < totBinsY_) {
1995       // if( TString(histoCovariance_->GetName()).Contains("CovarianceCotgTheta_Covariance") )
1996       covariances_[xIndex][yIndex].fill(a, b);
1997       // Should be used with the variance, in which case a==b
1998       if (varianceCheck_)
1999         histoVarianceCheck_[xIndex][yIndex]->Fill(a);
2000     }
2001   }
2002 
2003   using Histograms::Get;
2004   double Get(const double& x, const double& y) const {
2005     // Need to convert the (x,y) values to the array indeces
2006     int xIndex = getXindex(x);
2007     int yIndex = getYindex(y);
2008     // If the values exceed the limits of the histogram, return the border values
2009     if (xIndex < 0)
2010       xIndex = 0;
2011     if (xIndex >= totBinsX_)
2012       xIndex = totBinsX_ - 1;
2013     if (yIndex < 0)
2014       yIndex = 0;
2015     if (yIndex >= totBinsY_)
2016       yIndex = totBinsY_ - 1;
2017     return histoCovariance_->GetBinContent(xIndex + 1, yIndex + 1);
2018   }
2019 
2020   void Write() override {
2021     if (!readMode_) {
2022       std::cout << "writing: " << histoCovariance_->GetName() << std::endl;
2023       for (int xBin = 0; xBin < totBinsX_; ++xBin) {
2024         for (int yBin = 0; yBin < totBinsY_; ++yBin) {
2025           double covariance = covariances_[xBin][yBin].covariance();
2026           // Histogram bins start from 1
2027           // std::cout << "covariance["<<xBin<<"]["<<yBin<<"] with N = "<<covariances_[xBin][yBin].getN()<<" is: " << covariance << std::endl;
2028           histoCovariance_->SetBinContent(xBin + 1, yBin + 1, covariance);
2029         }
2030       }
2031       if (histoDir_ != nullptr)
2032         histoDir_->cd();
2033       TCanvas canvas(TString(histoCovariance_->GetName()) + "_canvas",
2034                      TString(histoCovariance_->GetTitle()) + " canvas",
2035                      1000,
2036                      800);
2037       canvas.Divide(2);
2038       canvas.cd(1);
2039       histoCovariance_->Draw("lego");
2040       canvas.cd(2);
2041       histoCovariance_->Draw("surf5");
2042       canvas.Write();
2043       histoCovariance_->Write();
2044 
2045       TDirectory* binsDir = nullptr;
2046       if (varianceCheck_) {
2047         if (histoDir_ != nullptr) {
2048           histoDir_->cd();
2049           if (binsDir == nullptr)
2050             binsDir = histoDir_->mkdir(name_ + "Bins");
2051           binsDir->cd();
2052         }
2053         for (int xBin = 0; xBin < totBinsX_; ++xBin) {
2054           for (int yBin = 0; yBin < totBinsY_; ++yBin) {
2055             histoVarianceCheck_[xBin][yBin]->Write();
2056           }
2057         }
2058       }
2059     }
2060   }
2061   void Clear() override {
2062     histoCovariance_->Clear();
2063     if (varianceCheck_) {
2064       for (int i = 0; i < totBinsX_; ++i) {
2065         for (int j = 0; j < totBinsY_; ++j) {
2066           histoVarianceCheck_[i][j]->Clear();
2067         }
2068       }
2069     }
2070   }
2071 
2072 protected:
2073   int getXindex(const double& x) const { return int((x - xMin_) / deltaX_ * totBinsX_); }
2074   int getYindex(const double& y) const { return int((y - yMin_) / deltaY_ * totBinsY_); }
2075   TH2D* histoCovariance_;
2076   Covariance** covariances_;
2077   int totBinsX_, totBinsY_, totBinsZ_;
2078   double xMin_, deltaX_, yMin_, deltaY_;
2079   bool readMode_;
2080   bool varianceCheck_;
2081   TH1D*** histoVarianceCheck_;
2082 };
2083 
2084 /**
2085  * This class uses the HCovariance histograms to compute the covariances between the two input muons kinematic quantities. </br>
2086  * The covariances are computed against pt and eta.
2087  */
2088 class HCovarianceVSParts : public Histograms {
2089 public:
2090   HCovarianceVSParts(TFile* outputFile, const TString& name, const double& ptMax) : Histograms(outputFile, name) {
2091     int totBinsX = 40;
2092     int totBinsY = 40;
2093     double etaMin = -3.;
2094     double etaMax = 3.;
2095     double ptMin = 0.;
2096 
2097     readMode_ = false;
2098 
2099     // Variances
2100     mapHisto_[name + "Pt"] =
2101         new HCovarianceVSxy(name + "Pt_", "Pt", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_, true);
2102     mapHisto_[name + "CotgTheta"] = new HCovarianceVSxy(
2103         name + "CotgTheta_", "CotgTheta", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_, true);
2104     mapHisto_[name + "Phi"] =
2105         new HCovarianceVSxy(name + "Phi_", "Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_, true);
2106     // Covariances
2107     mapHisto_[name + "Pt-CotgTheta"] = new HCovarianceVSxy(
2108         name + "Pt_CotgTheta_", "Pt-CotgTheta", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2109     mapHisto_[name + "Pt-Phi"] =
2110         new HCovarianceVSxy(name + "Pt_Phi_", "Pt-Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2111     mapHisto_[name + "CotgTheta-Phi"] = new HCovarianceVSxy(
2112         name + "CotgTheta_Phi_", "CotgTheta-Phi", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2113     mapHisto_[name + "Pt1-Pt2"] =
2114         new HCovarianceVSxy(name + "Pt1_Pt2_", "Pt1-Pt2", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2115     mapHisto_[name + "CotgTheta1-CotgTheta2"] = new HCovarianceVSxy(name + "CotgTheta1_CotgTheta2_",
2116                                                                     "CotgTheta1-CotgTheta2",
2117                                                                     totBinsX,
2118                                                                     ptMin,
2119                                                                     ptMax,
2120                                                                     totBinsY,
2121                                                                     etaMin,
2122                                                                     etaMax,
2123                                                                     histoDir_);
2124     mapHisto_[name + "Phi1-Phi2"] = new HCovarianceVSxy(
2125         name + "Phi1_Phi2_", "Phi1-Phi2", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2126     mapHisto_[name + "Pt12-CotgTheta21"] = new HCovarianceVSxy(
2127         name + "Pt12_CotgTheta21_", "Pt12-CotgTheta21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2128     mapHisto_[name + "Pt12-Phi21"] = new HCovarianceVSxy(
2129         name + "Pt12_Phi21_", "Pt12-Phi21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2130     mapHisto_[name + "CotgTheta12-Phi21"] = new HCovarianceVSxy(
2131         name + "CotgTheta12_Phi21_", "CotgTheta12-Phi21", totBinsX, ptMin, ptMax, totBinsY, etaMin, etaMax, histoDir_);
2132   }
2133 
2134   /// Constructor reading the histograms from file
2135   HCovarianceVSParts(const TString& inputFileName, const TString& name) {
2136     name_ = name;
2137     TFile* inputFile = new TFile(inputFileName, "READ");
2138     readMode_ = true;
2139 
2140     // Variances
2141     mapHisto_[name_ + "Pt"] = new HCovarianceVSxy(inputFile, name_ + "Pt_" + name_, name_);
2142     mapHisto_[name_ + "CotgTheta"] = new HCovarianceVSxy(inputFile, name_ + "CotgTheta_" + name_, name_);
2143     mapHisto_[name_ + "Phi"] = new HCovarianceVSxy(inputFile, name_ + "Phi_" + name_, name_);
2144     // Covariances
2145     mapHisto_[name_ + "Pt-CotgTheta"] = new HCovarianceVSxy(inputFile, name_ + "Pt_CotgTheta_" + name_, name_);
2146     mapHisto_[name_ + "Pt-Phi"] = new HCovarianceVSxy(inputFile, name_ + "Pt_Phi_" + name_, name_);
2147     mapHisto_[name_ + "CotgTheta-Phi"] = new HCovarianceVSxy(inputFile, name_ + "CotgTheta_Phi_" + name_, name_);
2148     mapHisto_[name_ + "Pt1-Pt2"] = new HCovarianceVSxy(inputFile, name_ + "Pt1_Pt2_" + name_, name_);
2149     mapHisto_[name_ + "CotgTheta1-CotgTheta2"] =
2150         new HCovarianceVSxy(inputFile, name_ + "CotgTheta1_CotgTheta2_" + name_, name_);
2151     mapHisto_[name_ + "Phi1-Phi2"] = new HCovarianceVSxy(inputFile, name_ + "Phi1_Phi2_" + name_, name_);
2152     mapHisto_[name_ + "Pt12-CotgTheta21"] = new HCovarianceVSxy(inputFile, name_ + "Pt12_CotgTheta21_" + name_, name_);
2153     mapHisto_[name_ + "Pt12-Phi21"] = new HCovarianceVSxy(inputFile, name_ + "Pt12_Phi21_" + name_, name_);
2154     mapHisto_[name_ + "CotgTheta12-Phi21"] =
2155         new HCovarianceVSxy(inputFile, name_ + "CotgTheta12_Phi21_" + name_, name_);
2156   }
2157 
2158   ~HCovarianceVSParts() override {
2159     for (std::map<TString, HCovarianceVSxy*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
2160          histo++) {
2161       delete (*histo).second;
2162     }
2163   }
2164 
2165   double Get(const reco::Particle::LorentzVector& recoP1, const TString& covarianceName) override {
2166     return mapHisto_[name_ + covarianceName]->Get(recoP1.pt(), recoP1.eta());
2167   }
2168 
2169   void Fill(const reco::Particle::LorentzVector& recoP1,
2170             const reco::Particle::LorentzVector& genP1,
2171             const reco::Particle::LorentzVector& recoP2,
2172             const reco::Particle::LorentzVector& genP2) override {
2173     double pt1 = recoP1.pt();
2174     double eta1 = recoP1.eta();
2175     double pt2 = recoP2.pt();
2176     double eta2 = recoP2.eta();
2177 
2178     double diffPt1 = (pt1 - genP1.pt()) / genP1.pt();
2179     double diffPt2 = (pt2 - genP2.pt()) / genP2.pt();
2180 
2181     double genTheta1 = genP1.theta();
2182     double genTheta2 = genP2.theta();
2183     double recoTheta1 = recoP1.theta();
2184     double recoTheta2 = recoP2.theta();
2185 
2186     double genCotgTheta1 = TMath::Cos(genTheta1) / (TMath::Sin(genTheta1));
2187     double genCotgTheta2 = TMath::Cos(genTheta2) / (TMath::Sin(genTheta2));
2188     double recoCotgTheta1 = TMath::Cos(recoTheta1) / (TMath::Sin(recoTheta1));
2189     double recoCotgTheta2 = TMath::Cos(recoTheta2) / (TMath::Sin(recoTheta2));
2190 
2191     // double diffCotgTheta1 = (recoCotgTheta1 - genCotgTheta1)/genCotgTheta1;
2192     // double diffCotgTheta2 = (recoCotgTheta2 - genCotgTheta2)/genCotgTheta2;
2193     double diffCotgTheta1 = recoCotgTheta1 - genCotgTheta1;
2194     double diffCotgTheta2 = recoCotgTheta2 - genCotgTheta2;
2195 
2196     // double diffPhi1 = (recoP1.phi() - genP1.phi())/genP1.phi();
2197     // double diffPhi2 = (recoP2.phi() - genP2.phi())/genP2.phi();
2198     double diffPhi1 = MuScleFitUtils::deltaPhiNoFabs(recoP1.phi(), genP1.phi());
2199     double diffPhi2 = MuScleFitUtils::deltaPhiNoFabs(recoP2.phi(), genP2.phi());
2200 
2201     // Fill the variances
2202     mapHisto_[name_ + "Pt"]->Fill(pt1, eta1, diffPt1, diffPt1);
2203     mapHisto_[name_ + "Pt"]->Fill(pt2, eta2, diffPt2, diffPt2);
2204     mapHisto_[name_ + "CotgTheta"]->Fill(pt1, eta1, diffCotgTheta1, diffCotgTheta1);
2205     mapHisto_[name_ + "CotgTheta"]->Fill(pt2, eta2, diffCotgTheta2, diffCotgTheta2);
2206     mapHisto_[name_ + "Phi"]->Fill(pt1, eta1, diffPhi1, diffPhi1);
2207     mapHisto_[name_ + "Phi"]->Fill(pt2, eta2, diffPhi2, diffPhi2);
2208 
2209     // Fill these histograms with both muons
2210     mapHisto_[name_ + "Pt-CotgTheta"]->Fill(pt1, eta1, diffPt1, diffCotgTheta1);
2211     mapHisto_[name_ + "Pt-CotgTheta"]->Fill(pt2, eta2, diffPt2, diffCotgTheta2);
2212     mapHisto_[name_ + "Pt-Phi"]->Fill(pt1, eta1, diffPt1, diffPhi1);
2213     mapHisto_[name_ + "Pt-Phi"]->Fill(pt2, eta2, diffPt2, diffPhi2);
2214     mapHisto_[name_ + "CotgTheta-Phi"]->Fill(pt1, eta1, diffCotgTheta1, diffPhi1);
2215     mapHisto_[name_ + "CotgTheta-Phi"]->Fill(pt2, eta2, diffCotgTheta2, diffPhi2);
2216 
2217     // We fill two (pt, eta) bins for each pair of values. The bin of the
2218     // first and of the second muon. This should take account for the
2219     // assumed symmetry between the exchange of the first with the second muon.
2220     mapHisto_[name_ + "Pt1-Pt2"]->Fill(pt1, eta1, diffPt1, diffPt2);
2221     mapHisto_[name_ + "Pt1-Pt2"]->Fill(pt2, eta2, diffPt1, diffPt2);
2222     mapHisto_[name_ + "CotgTheta1-CotgTheta2"]->Fill(pt1, eta1, diffCotgTheta1, diffCotgTheta2);
2223     mapHisto_[name_ + "CotgTheta1-CotgTheta2"]->Fill(pt2, eta2, diffCotgTheta1, diffCotgTheta2);
2224     mapHisto_[name_ + "Phi1-Phi2"]->Fill(pt1, eta1, diffPhi1, diffPhi2);
2225     mapHisto_[name_ + "Phi1-Phi2"]->Fill(pt2, eta2, diffPhi1, diffPhi2);
2226 
2227     // Fill the following histograms again for each muon (pt, eta) bin. Same
2228     // reason as in the previous case. If the symmetry is true, it does not
2229     // make any difference the order by which we fill the pt and cotgTheta combinations.
2230     mapHisto_[name_ + "Pt12-CotgTheta21"]->Fill(pt1, eta1, diffPt1, diffCotgTheta2);
2231     mapHisto_[name_ + "Pt12-CotgTheta21"]->Fill(pt2, eta2, diffPt2, diffCotgTheta1);
2232     mapHisto_[name_ + "Pt12-Phi21"]->Fill(pt1, eta1, diffPt1, diffPhi2);
2233     mapHisto_[name_ + "Pt12-Phi21"]->Fill(pt2, eta2, diffPt2, diffPhi1);
2234     mapHisto_[name_ + "CotgTheta12-Phi21"]->Fill(pt1, eta1, diffCotgTheta1, diffPhi2);
2235     mapHisto_[name_ + "CotgTheta12-Phi21"]->Fill(pt2, eta2, diffCotgTheta2, diffPhi1);
2236   }
2237   void Write() override {
2238     if (!readMode_) {
2239       histoDir_->cd();
2240       for (std::map<TString, HCovarianceVSxy*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
2241            histo++) {
2242         (*histo).second->Write();
2243       }
2244     }
2245   }
2246   void Clear() override {
2247     for (std::map<TString, HCovarianceVSxy*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
2248          histo++) {
2249       (*histo).second->Clear();
2250     }
2251   }
2252 
2253 protected:
2254   std::map<TString, HCovarianceVSxy*> mapHisto_;
2255   bool readMode_;
2256 };
2257 
2258 /**
2259  * A set of histograms for resolution.</br>
2260  * The fill method requires the two selected muons, their charges and the reconstructed and generated masses.
2261  * It evaluates the resolution on the mass vs:</br>
2262  * Pt of the pair</br>
2263  * DeltaEta of the pair</br>
2264  * DeltaPhi of the pair</br>
2265  * pt, eta and phi of the plus and minus muon separately</br>
2266  */
2267 class HMassResolutionVSPart : public Histograms {
2268 public:
2269   HMassResolutionVSPart(TFile* outputFile, const TString& name) : Histograms(outputFile, name) {
2270     // Kinematical variables
2271     nameSuffix_[0] = "Plus";
2272     nameSuffix_[1] = "Minus";
2273     TString titleSuffix[] = {" for mu+", " for mu-"};
2274 
2275     mapHisto_[name] = new TH1F(name, "#Delta M/M", 4000, -1, 1);
2276     mapHisto_[name + "VSPairPt"] =
2277         new HResolution(name + "VSPairPt", "resolution VS pt of the pair", 100, 0, 200, -1, 1, histoDir_);
2278     mapHisto_[name + "VSPairDeltaEta"] = new HResolution(
2279         name + "VSPairDeltaEta", "resolution VS #Delta#eta of the pair", 100, -0.1, 6.2, -1, 1, histoDir_);
2280     mapHisto_[name + "VSPairDeltaPhi"] = new HResolution(
2281         name + "VSPairDeltaPhi", "resolution VS #Delta#phi of the pair", 100, -0.1, 3.2, -1, 1, histoDir_);
2282 
2283     for (int i = 0; i < 2; ++i) {
2284       mapHisto_[name + "VSPt" + nameSuffix_[i]] = new HResolution(
2285           name + "VSPt" + nameSuffix_[i], "resolution VS pt" + titleSuffix[i], 100, 0, 200, -1, 1, histoDir_);
2286       mapHisto_[name + "VSEta" + nameSuffix_[i]] = new HResolution(
2287           name + "VSEta" + nameSuffix_[i], "resolution VS #eta" + titleSuffix[i], 100, -3, 3, -1, 1, histoDir_);
2288       mapHisto_[name + "VSPhi" + nameSuffix_[i]] = new HResolution(
2289           name + "VSPhi" + nameSuffix_[i], "resolution VS #phi" + titleSuffix[i], 100, -3.2, 3.2, -1, 1, histoDir_);
2290     }
2291 
2292     // single particles histograms
2293     muMinus = std::make_unique<HDelta>("muMinus");
2294     muPlus = std::make_unique<HDelta>("muPlus");
2295   }
2296 
2297   ~HMassResolutionVSPart() override {
2298     for (std::map<TString, TH1*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end(); histo++) {
2299       delete (*histo).second;
2300     }
2301   }
2302 
2303   void Fill(const reco::Particle::LorentzVector& recoP1,
2304             const int charge1,
2305             const reco::Particle::LorentzVector& genP1,
2306             const reco::Particle::LorentzVector& recoP2,
2307             const int charge2,
2308             const reco::Particle::LorentzVector& genP2,
2309             const double& recoMass,
2310             const double& genMass) override {
2311     muMinus->Fill(recoP1, genP1);
2312     muPlus->Fill(recoP2, genP2);
2313 
2314     Fill(recoP1, charge1, recoP2, charge2, recoMass, genMass);
2315   }
2316 
2317   void Fill(const reco::Particle::LorentzVector& recoP1,
2318             const int charge1,
2319             // const reco::Particle::LorentzVector & genP1,
2320             const reco::Particle::LorentzVector& recoP2,
2321             const int charge2,
2322             // const reco::Particle::LorentzVector & genP2,
2323             const double& recoMass,
2324             const double& genMass) override {
2325     if (charge1 == charge2)
2326       std::cout << "Error: must get two opposite charge particles" << std::endl;
2327 
2328     double massRes = (recoMass - genMass) / genMass;
2329 
2330     reco::Particle::LorentzVector recoPair(recoP1 + recoP2);
2331     double pairPt = recoPair.Pt();
2332 
2333     double recoPt[2] = {recoP1.Pt(), recoP2.Pt()};
2334     double recoEta[2] = {recoP1.Eta(), recoP2.Eta()};
2335     double recoPhi[2] = {recoP1.Phi(), recoP2.Phi()};
2336 
2337     // std::cout << "pairPt = " << pairPt << ", massRes = ("<<recoMass<<" - "<<genMass<<")/"<<genMass<<" = " << massRes
2338     //      << ", recoPt[0] = " << recoPt[0] << ", recoPt[1] = " << recoPt[1] << std::endl;
2339 
2340     // Index of the histogram. If the muons have charge1 = -1 and charge2 = 1, they already have the
2341     // correct histogram indeces. Otherwise, swap the indeces.
2342     // In any case the negative muon has index = 0  and the positive muon has index = 1.
2343     int id[2] = {0, 1};
2344     if (charge1 > 0) {
2345       id[0] = 1;
2346       id[1] = 0;
2347     }
2348 
2349     double pairDeltaEta = fabs(recoEta[0] - recoEta[1]);
2350     double pairDeltaPhi = MuScleFitUtils::deltaPhi(recoPhi[0], recoPhi[1]);
2351 
2352     mapHisto_[name_]->Fill(massRes);
2353     mapHisto_[name_ + "VSPairPt"]->Fill(pairPt, massRes);
2354     mapHisto_[name_ + "VSPairDeltaEta"]->Fill(pairDeltaEta, massRes);
2355     mapHisto_[name_ + "VSPairDeltaPhi"]->Fill(pairDeltaPhi, massRes);
2356 
2357     // Resolution vs single muon quantities
2358     // ------------------------------------
2359     int index = 0;
2360     for (int i = 0; i < 2; ++i) {
2361       index = id[i];
2362       mapHisto_[name_ + "VSPt" + nameSuffix_[i]]->Fill(recoPt[index], massRes);  // EM [index] or [i]???
2363       mapHisto_[name_ + "VSEta" + nameSuffix_[i]]->Fill(recoEta[index], massRes);
2364       mapHisto_[name_ + "VSPhi" + nameSuffix_[i]]->Fill(recoPhi[index], massRes);
2365     }
2366   }
2367 
2368   void Write() override {
2369     histoDir_->cd();
2370     for (std::map<TString, TH1*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end(); histo++) {
2371       (*histo).second->Write();
2372     }
2373     // Create the new dir and cd into it
2374     (histoDir_->mkdir("singleMuonsVSgen"))->cd();
2375     muMinus->Write();
2376     muPlus->Write();
2377   }
2378 
2379   void Clear() override {
2380     for (std::map<TString, TH1*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end(); histo++) {
2381       (*histo).second->Clear();
2382     }
2383     muMinus->Clear();
2384     muPlus->Clear();
2385   }
2386 
2387   //   HMassResolutionVSPart(const TString & name, TFile* file){
2388   //     TString nameSuffix[] = {"Plus", "Minus"};
2389   //     name_ = name;
2390   //     hReso                    = (TH1F *)     file->Get(name+"_Reso");
2391   //     hResoVSPairPt            = (TH2F *)     file->Get(name+"_ResoVSPairPt");
2392   //     hResoVSPairDeltaEta      = (TH2F *)     file->Get(name+"_ResoVSPairDeltaEta");
2393   //     hResoVSPairDeltaPhi      = (TH2F *)     file->Get(name+"_ResoVSPairDeltaPhi");
2394   //     for( int i=0; i<2; ++i ) {
2395   //       hResoVSPt[i]           = (TH2F *)     file->Get(name+"_ResoVSPt"+nameSuffix[i]);
2396   //       hResoVSEta[i]          = (TH2F *)     file->Get(name+"_ResoVSEta"+nameSuffix[i]);
2397   //       hResoVSPhi[i]          = (TH2F *)     file->Get(name+"_ResoVSPhi"+nameSuffix[i]);
2398   //     }
2399   //   }
2400 
2401 protected:
2402   std::map<TString, TH1*> mapHisto_;
2403   TString nameSuffix_[2];
2404   std::unique_ptr<HDelta> muMinus;
2405   std::unique_ptr<HDelta> muPlus;
2406 };
2407 
2408 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!