Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:33:09

0001 #ifndef Validation_RecoMuon_Histograms_H
0002 #define Validation_RecoMuon_Histograms_H
0003 
0004 /** \class Histograms
0005  *  No description available.
0006  *
0007  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0008  */
0009 
0010 #include "TH1F.h"
0011 #include "TH2F.h"
0012 #include "TString.h"
0013 #include "TFile.h"
0014 
0015 #include "DQMServices/Core/interface/DQMStore.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 
0018 #include "DataFormats/GeometryVector/interface/Pi.h"
0019 #include <iostream>
0020 #include <vector>
0021 #include <cmath>
0022 
0023 class HTrackVariables {
0024 public:
0025   typedef dqm::legacy::DQMStore DQMStore;
0026   typedef dqm::legacy::MonitorElement MonitorElement;
0027 
0028   HTrackVariables(DQMStore::IBooker &ibooker, std::string dirName_, std::string name, std::string whereIs = "")
0029       : theName(name), where(whereIs) {
0030     ibooker.cd();
0031     const std::string &dirName = dirName_;
0032     //dirName+="/";
0033     //dirName+=name.c_str();
0034 
0035     ibooker.setCurrentFolder(dirName);
0036 
0037     hEta = ibooker.book1D(theName + "_Eta_" + where, "#eta at " + where, 120, -3., 3.);
0038     hPhi = ibooker.book1D(theName + "_Phi_" + where, "#phi at " + where, 100, -Geom::pi(), Geom::pi());
0039     hP = ibooker.book1D(theName + "_P_" + where, "p_{t} at " + where, 1000, 0, 2000);
0040     hPt = ibooker.book1D(theName + "_Pt_" + where, "p_{t} at " + where, 1000, 0, 2000);
0041     hCharge = ibooker.book1D(theName + "_charge_" + where, "Charge at " + where, 4, -2., 2.);
0042 
0043     hEtaVsGen = ibooker.book1D(theName + "_EtaVsGen_" + where, "#eta at " + where, 120, -3., 3.);
0044     hPhiVsGen = ibooker.book1D(theName + "_PhiVsGen_" + where, "#phi at " + where, 100, -Geom::pi(), Geom::pi());
0045     hPtVsGen = ibooker.book1D(theName + "_PtVsGen_" + where, "p_{t} at " + where, 1000, 0, 2000);
0046 
0047     hDeltaR = ibooker.book1D(theName + "_DeltaR_" + where, "Delta R w.r.t. sim track for " + where, 1000, 0, 20);
0048 
0049     theEntries = 0;
0050   }
0051 
0052   /*
0053   HTrackVariables(std::string name, TFile* file):theName(name.c_str()){ 
0054     
0055     hEta = dynamic_cast<TH1F*>( file->Get(theName+"_Eta_"+where));
0056     hPhi = dynamic_cast<TH1F*>( file->Get(theName+"_Phi_"+where));
0057     hP   = dynamic_cast<TH1F*>( file->Get(theName+"_P_"+where));
0058     hPt  = dynamic_cast<TH1F*>( file->Get(theName+"_Pt_"+where));
0059     hCharge = dynamic_cast<TH1F*>( file->Get(theName+"_charge_"+where)); 
0060   }
0061   */
0062 
0063   ~HTrackVariables() {}
0064 
0065   MonitorElement *eta() { return hEta; }
0066   MonitorElement *phi() { return hPhi; }
0067   MonitorElement *p() { return hP; }
0068   MonitorElement *pt() { return hPt; }
0069   MonitorElement *charge() { return hCharge; }
0070   int entries() { return theEntries; }
0071 
0072   void Fill(double p, double pt, double eta, double phi, double charge) {
0073     hEta->Fill(eta);
0074     hPhi->Fill(phi);
0075     hP->Fill(p);
0076     hPt->Fill(pt);
0077     hCharge->Fill(charge);
0078     ++theEntries;
0079   }
0080 
0081   void Fill(double pt, double eta, double phi) {
0082     hEtaVsGen->Fill(eta);
0083     hPhiVsGen->Fill(phi);
0084     hPtVsGen->Fill(pt);
0085   }
0086 
0087   void FillDeltaR(double deltaR) { hDeltaR->Fill(deltaR); }
0088 
0089   double computeEfficiency(HTrackVariables *sim, DQMStore::IBooker &ibooker) {
0090     efficiencyHistos.push_back(computeEfficiency(hEtaVsGen, sim->eta(), ibooker));
0091     efficiencyHistos.push_back(computeEfficiency(hPhiVsGen, sim->phi(), ibooker));
0092     efficiencyHistos.push_back(computeEfficiency(hPtVsGen, sim->pt(), ibooker));
0093 
0094     double efficiency = 100 * entries() / sim->entries();
0095     return efficiency;
0096   }
0097 
0098   MonitorElement *computeEfficiency(MonitorElement *reco, MonitorElement *sim, DQMStore::IBooker &ibooker) {
0099     TH1F *hReco = reco->getTH1F();
0100     TH1F *hSim = sim->getTH1F();
0101 
0102     std::string name = hReco->GetName();
0103     std::string title = hReco->GetTitle();
0104 
0105     MonitorElement *me = ibooker.book1D("Eff_" + name,
0106                                         "Efficiecny as a function of " + title,
0107                                         hSim->GetNbinsX(),
0108                                         hSim->GetXaxis()->GetXmin(),
0109                                         hSim->GetXaxis()->GetXmax());
0110 
0111     me->getTH1F()->Divide(hReco, hSim, 1., 1., "b");
0112 
0113     // Set the error accordingly to binomial statistics
0114     int nBinsEta = me->getTH1F()->GetNbinsX();
0115     for (int bin = 1; bin <= nBinsEta; bin++) {
0116       float nSimHit = hSim->GetBinContent(bin);
0117       float eff = me->getTH1F()->GetBinContent(bin);
0118       float error = 0;
0119       if (nSimHit != 0 && eff <= 1) {
0120         error = sqrt(eff * (1 - eff) / nSimHit);
0121       }
0122       me->getTH1F()->SetBinError(bin, error);
0123     }
0124 
0125     return me;
0126   }
0127 
0128 private:
0129   std::string theName;
0130   std::string where;
0131 
0132   int theEntries;
0133 
0134   MonitorElement *hEta;
0135   MonitorElement *hPhi;
0136   MonitorElement *hP;
0137   MonitorElement *hPt;
0138   MonitorElement *hCharge;
0139 
0140   MonitorElement *hEtaVsGen;
0141   MonitorElement *hPhiVsGen;
0142   MonitorElement *hPtVsGen;
0143 
0144   MonitorElement *hDeltaR;
0145 
0146   std::vector<MonitorElement *> efficiencyHistos;
0147 };
0148 
0149 class HResolution {
0150 public:
0151   typedef dqm::legacy::DQMStore DQMStore;
0152   typedef dqm::legacy::MonitorElement MonitorElement;
0153 
0154   HResolution(DQMStore::IBooker &ibooker, std::string dirName_, std::string name, std::string whereIs)
0155       : theName(name), where(whereIs) {
0156     ibooker.cd();
0157     const std::string &dirName = dirName_;
0158     //dirName+="/";
0159     //dirName+=name.c_str();
0160 
0161     ibooker.setCurrentFolder(dirName);
0162 
0163     double eta = 15.;
0164     int neta = 800;
0165     double phi = 12.;
0166     int nphi = 400;
0167     double pt = 60.;
0168     int npt = 2000;
0169 
0170     hEta = ibooker.book1D(theName + "_Eta_" + where, "#eta " + theName, neta, -eta, eta);  // 400
0171     hPhi = ibooker.book1D(theName + "_Phi_" + where, "#phi " + theName, nphi, -phi, phi);  // 100
0172 
0173     hP = ibooker.book1D(theName + "_P_" + where, "P " + theName, 400, -4, 4);          // 200
0174     hPt = ibooker.book1D(theName + "_Pt_" + where, "P_{t} " + theName, npt, -pt, pt);  // 200
0175 
0176     hCharge = ibooker.book1D(theName + "_charge_" + where, "Charge " + theName, 4, -2., 2.);
0177 
0178     h2Eta = ibooker.book2D(
0179         theName + "_Eta_vs_Eta" + where, "#eta " + theName + " as a function of #eta", 200, -2.5, 2.5, neta, -eta, eta);
0180     h2Phi = ibooker.book2D(theName + "_Phi_vs_Phi" + where,
0181                            "#phi " + theName + " as a function of #phi",
0182                            100,
0183                            -Geom::pi(),
0184                            Geom::pi(),
0185                            nphi,
0186                            -phi,
0187                            phi);
0188 
0189     h2P =
0190         ibooker.book2D(theName + "_P_vs_P" + where, "P " + theName + " as a function of P", 1000, 0, 2000, 400, -4, 4);
0191     h2Pt = ibooker.book2D(
0192         theName + "_Pt_vs_Pt" + where, "P_{t} " + theName + " as a function of P_{t}", 1000, 0, 2000, npt, -pt, pt);
0193 
0194     h2PtVsEta = ibooker.book2D(
0195         theName + "_Pt_vs_Eta" + where, "P_{t} " + theName + " as a function of #eta", 200, -2.5, 2.5, npt, -pt, pt);
0196     h2PtVsPhi = ibooker.book2D(theName + "_Pt_vs_Phi" + where,
0197                                "P_{t} " + theName + " as a function of #phi",
0198                                100,
0199                                -Geom::pi(),
0200                                Geom::pi(),
0201                                npt,
0202                                -pt,
0203                                pt);
0204 
0205     h2EtaVsPt = ibooker.book2D(
0206         theName + "_Eta_vs_Pt" + where, "#eta " + theName + " as a function of P_{t}", 1000, 0, 2000, neta, -eta, eta);
0207     h2EtaVsPhi = ibooker.book2D(theName + "_Eta_vs_Phi" + where,
0208                                 "#eta " + theName + " as a function of #phi",
0209                                 100,
0210                                 -Geom::pi(),
0211                                 Geom::pi(),
0212                                 neta,
0213                                 -eta,
0214                                 eta);
0215 
0216     h2PhiVsPt = ibooker.book2D(
0217         theName + "_Phi_vs_Pt" + where, "#phi " + theName + " as a function of P_{t}", 1000, 0, 2000, nphi, -phi, phi);
0218     h2PhiVsEta = ibooker.book2D(
0219         theName + "_Phi_vs_Eta" + where, "#phi " + theName + " as a function of #eta", 200, -2.5, 2.5, nphi, -phi, phi);
0220   }
0221 
0222   HResolution(DQMStore::IBooker &ibooker, std::string name, TFile *file) : theName(name) {
0223     //    dynamic_cast<TH1F*>( file->Get(theName+"") );
0224   }
0225 
0226   ~HResolution() {}
0227 
0228   void Fill(
0229       double p, double pt, double eta, double phi, double rp, double rpt, double reta, double rphi, double rcharge) {
0230     Fill(rp, rpt, reta, rphi, rcharge);
0231 
0232     h2Eta->Fill(eta, reta);
0233     h2Phi->Fill(phi, rphi);
0234 
0235     h2P->Fill(p, rp);
0236     h2Pt->Fill(pt, rpt);
0237 
0238     h2PtVsEta->Fill(eta, rpt);
0239     h2PtVsPhi->Fill(phi, rpt);
0240 
0241     h2EtaVsPt->Fill(pt, reta);
0242     h2EtaVsPhi->Fill(phi, reta);
0243 
0244     h2PhiVsPt->Fill(pt, rphi);
0245     h2PhiVsEta->Fill(eta, rphi);
0246   }
0247 
0248   void Fill(double p, double pt, double eta, double phi, double rp, double rpt) {
0249     hP->Fill(rp);
0250     hPt->Fill(rpt);
0251 
0252     h2P->Fill(p, rp);
0253     // h2PVsEta->Fill(eta,rp);
0254     // h2PVsPhi->Fill(phi,rp);
0255 
0256     h2Pt->Fill(pt, rpt);
0257     h2PtVsEta->Fill(eta, rpt);
0258     h2PtVsPhi->Fill(phi, rpt);
0259   }
0260 
0261   void Fill(double rp, double rpt, double reta, double rphi, double rcharge) {
0262     hEta->Fill(reta);
0263     hPhi->Fill(rphi);
0264 
0265     hP->Fill(rp);
0266     hPt->Fill(rpt);
0267 
0268     hCharge->Fill(rcharge);
0269   }
0270 
0271 private:
0272   std::string theName;
0273   std::string where;
0274 
0275   MonitorElement *hEta;
0276   MonitorElement *hPhi;
0277 
0278   MonitorElement *hP;
0279   MonitorElement *hPt;
0280 
0281   MonitorElement *hCharge;
0282 
0283   MonitorElement *h2Eta;
0284   MonitorElement *h2Phi;
0285 
0286   MonitorElement *h2P;
0287   MonitorElement *h2Pt;
0288 
0289   MonitorElement *h2PtVsEta;
0290   MonitorElement *h2PtVsPhi;
0291 
0292   MonitorElement *h2EtaVsPt;
0293   MonitorElement *h2EtaVsPhi;
0294 
0295   MonitorElement *h2PhiVsPt;
0296   MonitorElement *h2PhiVsEta;
0297 };
0298 
0299 class HResolution1DRecHit {
0300 public:
0301   typedef dqm::legacy::DQMStore DQMStore;
0302   typedef dqm::legacy::MonitorElement MonitorElement;
0303 
0304   HResolution1DRecHit(DQMStore::IBooker &ibooker, std::string name) : theName(name) {
0305     // Position, sigma, residual, pull
0306     hResX = ibooker.book1D(theName + "_X_Res", "X residual", 5000, -0.5, 0.5);
0307     hResY = ibooker.book1D(theName + "_Y_Res", "Y residual", 5000, -1., 1.);
0308     hResZ = ibooker.book1D(theName + "_Z_Res", "Z residual", 5000, -1.5, 1.5);
0309 
0310     hResXVsEta = ibooker.book2D(theName + "_XResVsEta", "X residual vs eta", 200, -2.5, 2.5, 5000, -1.5, 1.5);
0311     hResYVsEta = ibooker.book2D(theName + "_YResVsEta", "Y residual vs eta", 200, -2.5, 2.5, 5000, -1., 1.);
0312     hResZVsEta = ibooker.book2D(theName + "_ZResVsEta", "Z residual vs eta", 200, -2.5, 2.5, 5000, -1.5, 1.5);
0313 
0314     hXResVsPhi =
0315         ibooker.book2D(theName + "_XResVsPhi", "X residual vs phi", 100, -Geom::pi(), Geom::pi(), 5000, -0.5, 0.5);
0316     hYResVsPhi =
0317         ibooker.book2D(theName + "_YResVsPhi", "Y residual vs phi", 100, -Geom::pi(), Geom::pi(), 5000, -1., 1.);
0318     hZResVsPhi =
0319         ibooker.book2D(theName + "_ZResVsPhi", "Z residual vs phi", 100, -Geom::pi(), Geom::pi(), 5000, -1.5, 1.5);
0320 
0321     hXResVsPos = ibooker.book2D(theName + "_XResVsPos", "X residual vs position", 10000, -750, 750, 5000, -0.5, 0.5);
0322     hYResVsPos = ibooker.book2D(theName + "_YResVsPos", "Y residual vs position", 10000, -740, 740, 5000, -1., 1.);
0323     hZResVsPos = ibooker.book2D(theName + "_ZResVsPos", "Z residual vs position", 10000, -1100, 1100, 5000, -1.5, 1.5);
0324 
0325     hXPull = ibooker.book1D(theName + "_XPull", "X pull", 600, -2, 2);
0326     hYPull = ibooker.book1D(theName + "_YPull", "Y pull", 600, -2, 2);
0327     hZPull = ibooker.book1D(theName + "_ZPull", "Z pull", 600, -2, 2);
0328 
0329     hXPullVsPos = ibooker.book2D(theName + "_XPullVsPos", "X pull vs position", 10000, -750, 750, 600, -2, 2);
0330     hYPullVsPos = ibooker.book2D(theName + "_YPullVsPos", "Y pull vs position", 10000, -740, 740, 600, -2, 2);
0331     hZPullVsPos = ibooker.book2D(theName + "_ZPullVsPos", "Z pull vs position", 10000, -1100, 1100, 600, -2, 2);
0332   }
0333 
0334   HResolution1DRecHit(DQMStore::IBooker &ibooker, const TString &name_, TFile *file) {}
0335 
0336   ~HResolution1DRecHit() {}
0337 
0338   void Fill(double x,
0339             double y,
0340             double z,
0341             double dx,
0342             double dy,
0343             double dz,
0344             double errx,
0345             double erry,
0346             double errz,
0347             double eta,
0348             double phi) {
0349     double rx = dx / x, ry = dy / y, rz = dz / z;
0350 
0351     hResX->Fill(rx);
0352     hResY->Fill(ry);
0353     hResZ->Fill(rz);
0354 
0355     hResXVsEta->Fill(eta, rx);
0356     hResYVsEta->Fill(eta, ry);
0357     hResZVsEta->Fill(eta, rz);
0358 
0359     hXResVsPhi->Fill(phi, rx);
0360     hYResVsPhi->Fill(phi, ry);
0361     hZResVsPhi->Fill(phi, rz);
0362 
0363     hXResVsPos->Fill(x, rx);
0364     hYResVsPos->Fill(y, ry);
0365     hZResVsPos->Fill(z, rz);
0366 
0367     if (errx < 1e-6)
0368       std::cout << "NO proper error set for X: " << errx << std::endl;
0369     else {
0370       hXPull->Fill(dx / errx);
0371       hXPullVsPos->Fill(x, dx / errx);
0372     }
0373 
0374     if (erry < 1e-6)
0375       std::cout << "NO proper error set for Y: " << erry << std::endl;
0376     else {
0377       hYPull->Fill(dy / erry);
0378       hYPullVsPos->Fill(y, dy / erry);
0379     }
0380     if (errz < 1e-6)
0381       std::cout << "NO proper error set for Z: " << errz << std::endl;
0382     else {
0383       hZPull->Fill(dz / errz);
0384       hZPullVsPos->Fill(z, dz / errz);
0385     }
0386   }
0387 
0388 public:
0389   std::string theName;
0390 
0391   MonitorElement *hResX;
0392   MonitorElement *hResY;
0393   MonitorElement *hResZ;
0394 
0395   MonitorElement *hResXVsEta;
0396   MonitorElement *hResYVsEta;
0397   MonitorElement *hResZVsEta;
0398 
0399   MonitorElement *hXResVsPhi;
0400   MonitorElement *hYResVsPhi;
0401   MonitorElement *hZResVsPhi;
0402 
0403   MonitorElement *hXResVsPos;
0404   MonitorElement *hYResVsPos;
0405   MonitorElement *hZResVsPos;
0406 
0407   MonitorElement *hXPull;
0408   MonitorElement *hYPull;
0409   MonitorElement *hZPull;
0410 
0411   MonitorElement *hXPullVsPos;
0412   MonitorElement *hYPullVsPos;
0413   MonitorElement *hZPullVsPos;
0414 };
0415 #endif