CMS 3D CMS Logo

Histograms.h

Go to the documentation of this file.
00001 #ifndef Validation_RecoMuon_Histograms_H
00002 #define Validation_RecoMuon_Histograms_H
00003 
00012 #include "TH1F.h"
00013 #include "TH2F.h"
00014 #include "TString.h"
00015 #include "TFile.h"
00016 
00017 #include "DQMServices/Core/interface/DQMStore.h"
00018 #include "DQMServices/Core/interface/MonitorElement.h"
00019 #include "FWCore/ServiceRegistry/interface/Service.h"
00020 
00021 #include "DataFormats/GeometryVector/interface/Pi.h"
00022 #include <iostream>
00023 #include <vector>
00024 #include <math.h>
00025 
00026 class HTrackVariables{
00027 public:
00028   
00029   HTrackVariables(std::string dirName_, std::string name,std::string whereIs =""):theName(name.c_str()),where(whereIs.c_str()){
00030     dbe_ = edm::Service<DQMStore>().operator->();
00031     dbe_->cd();
00032     std::string dirName=dirName_;
00033     //dirName+="/";
00034     //dirName+=name.c_str();
00035     
00036     dbe_->setCurrentFolder(dirName.c_str());
00037     
00038     hEta = dbe_->book1D(theName+"_Eta_"+where,"#eta at "+where,120,-3.,3.);
00039     hPhi = dbe_->book1D(theName+"_Phi_"+where,"#phi at "+where,100,-Geom::pi(),Geom::pi());
00040     hP   = dbe_->book1D(theName+"_P_"+where,"p_{t} at "+where,1000,0,2000);
00041     hPt  = dbe_->book1D(theName+"_Pt_"+where,"p_{t} at "+where,1000,0,2000);
00042     hCharge = dbe_->book1D(theName+"_charge_"+where,"Charge at "+where,4,-2.,2.);
00043 
00044     hEtaVsGen = dbe_->book1D(theName+"_EtaVsGen_"+where,"#eta at "+where,120,-3.,3.);
00045     hPhiVsGen = dbe_->book1D(theName+"_PhiVsGen_"+where,"#phi at "+where,100,-Geom::pi(),Geom::pi());
00046     hPtVsGen  = dbe_->book1D(theName+"_PtVsGen_"+where,"p_{t} at "+where,1000,0,2000);
00047 
00048     hDeltaR = dbe_->book1D(theName+"_DeltaR_"+where,"Delta R w.r.t. sim track for "+where,1000,0,20);
00049 
00050     theEntries = 0;
00051   }
00052   
00053   /*
00054   HTrackVariables(std::string name, TFile* file):theName(name.c_str()){ 
00055     
00056     hEta = dynamic_cast<TH1F*>( file->Get(theName+"_Eta_"+where));
00057     hPhi = dynamic_cast<TH1F*>( file->Get(theName+"_Phi_"+where));
00058     hP   = dynamic_cast<TH1F*>( file->Get(theName+"_P_"+where));
00059     hPt  = dynamic_cast<TH1F*>( file->Get(theName+"_Pt_"+where));
00060     hCharge = dynamic_cast<TH1F*>( file->Get(theName+"_charge_"+where)); 
00061   }
00062   */
00063   
00064   ~HTrackVariables(){}
00065 
00066   MonitorElement *eta() {return hEta;}
00067   MonitorElement *phi() {return hPhi;}
00068   MonitorElement *p() {return hP;}
00069   MonitorElement *pt() {return hPt;}
00070   MonitorElement *charge() {return hCharge;}
00071   int entries() {return theEntries;}
00072 
00073   void Fill(double p, double pt, double eta, double phi, double charge){
00074     hEta->Fill(eta);
00075     hPhi->Fill(phi);
00076     hP->Fill(p);
00077     hPt->Fill(pt);
00078     hCharge->Fill(charge);
00079     ++theEntries;
00080   }
00081 
00082   void Fill(double pt, double eta, double phi){
00083     hEtaVsGen->Fill(eta);
00084     hPhiVsGen->Fill(phi);
00085     hPtVsGen->Fill(pt);
00086   }
00087 
00088   void FillDeltaR(double deltaR){
00089     hDeltaR->Fill(deltaR);
00090   }
00091   
00092   double computeEfficiency(HTrackVariables *sim){
00093     
00094     efficiencyHistos.push_back(computeEfficiency(hEtaVsGen,sim->eta()));
00095     efficiencyHistos.push_back(computeEfficiency(hPhiVsGen,sim->phi()));
00096     //    efficiencyHistos.push_back(computeEfficiency(p(),sim->p()));
00097     efficiencyHistos.push_back(computeEfficiency(hPtVsGen,sim->pt()));
00098     //    efficiencyHistos.push_back(computeEfficiency(charge(),sim->charge()));
00099 
00100     double efficiency = 100*entries()/sim->entries();
00101     return efficiency;
00102   }
00103 
00104   MonitorElement* computeEfficiency(MonitorElement *reco, MonitorElement *sim){
00105     
00106     TH1F* hReco = reco->getTH1F();
00107     TH1F* hSim  = sim->getTH1F();
00108 
00109     std::string name = hReco->GetName();
00110     std::string title = hReco->GetTitle();
00111     
00112     MonitorElement * me = dbe_->book1D(
00113                                        "Eff_"+name,
00114                                        "Efficiecny as a function of "+title,
00115                                        hSim->GetNbinsX(),
00116                                        hSim->GetXaxis()->GetXmin(),
00117                                        hSim->GetXaxis()->GetXmax()
00118                                        );
00119     
00120     me->getTH1F()->Divide(hReco,hSim,1.,1.,"b");
00121     
00122     // Set the error accordingly to binomial statistics
00123     int nBinsEta = me->getTH1F()->GetNbinsX();
00124     for(int bin = 1; bin <=  nBinsEta; bin++) {
00125       float nSimHit = hSim->GetBinContent(bin);
00126       float eff = me->getTH1F()->GetBinContent(bin);
00127       float error = 0;
00128       if(nSimHit != 0 && eff <= 1) {
00129         error = sqrt(eff*(1-eff)/nSimHit);
00130       }
00131       me->getTH1F()->SetBinError(bin, error);
00132     }
00133     
00134     return me;
00135   }
00136   
00137   
00138  private:
00139   DQMStore* dbe_;
00140 
00141   std::string theName;
00142   std::string where;
00143   
00144   int theEntries;
00145 
00146   MonitorElement *hEta;
00147   MonitorElement *hPhi;
00148   MonitorElement *hP;
00149   MonitorElement *hPt;
00150   MonitorElement *hCharge;
00151 
00152   MonitorElement *hEtaVsGen;
00153   MonitorElement *hPhiVsGen;
00154   MonitorElement *hPtVsGen;
00155   
00156   MonitorElement* hDeltaR;
00157   
00158   std::vector<MonitorElement*> efficiencyHistos;
00159 
00160   
00161 };
00162 
00163 
00164 class HResolution {
00165 public:
00166   
00167   HResolution(std::string dirName_,std::string name,std::string whereIs):theName(name.c_str()),where(whereIs.c_str()){
00168     
00169     dbe_ = edm::Service<DQMStore>().operator->();
00170     dbe_->cd();
00171     std::string dirName=dirName_;
00172     //dirName+="/";
00173     //dirName+=name.c_str();
00174     
00175     dbe_->setCurrentFolder(dirName.c_str());
00176     
00177     double eta = 15.; int neta = 800;
00178     double phi = 12.; int nphi = 400;
00179     double pt = 60.; int npt = 2000;
00180 
00181     hEta = dbe_->book1D(theName+"_Eta_"+where,"#eta "+theName,neta,-eta,eta); // 400
00182     hPhi = dbe_->book1D(theName+"_Phi_"+where,"#phi "+theName,nphi,-phi,phi); // 100
00183 
00184     hP = dbe_->book1D(theName+"_P_"+where,"P "+theName,400,-4,4);  // 200
00185     hPt = dbe_->book1D(theName+"_Pt_"+where,"P_{t} "+theName,npt,-pt,pt); // 200
00186 
00187     hCharge = dbe_->book1D(theName+"_charge_"+where,"Charge "+theName,4,-2.,2.);
00188 
00189 
00190     h2Eta = dbe_->book2D(theName+"_Eta_vs_Eta"+where,"#eta "+theName+" as a function of #eta",200,-2.5,2.5,neta,-eta,eta);
00191     h2Phi = dbe_->book2D(theName+"_Phi_vs_Phi"+where,"#phi "+theName+" as a function of #phi",100,-Geom::pi(),Geom::pi(),nphi,-phi,phi);
00192     
00193     h2P = dbe_->book2D(theName+"_P_vs_P"+where,"P "+theName+" as a function of P",1000,0,2000,400,-4,4);
00194     h2Pt = dbe_->book2D(theName+"_Pt_vs_Pt"+where,"P_{t} "+theName+" as a function of P_{t}",1000,0,2000,npt,-pt,pt);
00195     
00196     h2PtVsEta = dbe_->book2D(theName+"_Pt_vs_Eta"+where,"P_{t} "+theName+" as a function of #eta",200,-2.5,2.5,npt,-pt,pt);
00197     h2PtVsPhi = dbe_->book2D(theName+"_Pt_vs_Phi"+where,"P_{t} "+theName+" as a function of #phi",100,-Geom::pi(),Geom::pi(),npt,-pt,pt);
00198 
00199     h2EtaVsPt = dbe_->book2D(theName+"_Eta_vs_Pt"+where,"#eta "+theName+" as a function of P_{t}",1000,0,2000,neta,-eta,eta);
00200     h2EtaVsPhi = dbe_->book2D(theName+"_Eta_vs_Phi"+where,"#eta "+theName+" as a function of #phi",100,-Geom::pi(),Geom::pi(),neta,-eta,eta);
00201     
00202     h2PhiVsPt = dbe_->book2D(theName+"_Phi_vs_Pt"+where,"#phi "+theName+" as a function of P_{t}",1000,0,2000,nphi,-phi,phi);
00203     h2PhiVsEta = dbe_->book2D(theName+"_Phi_vs_Eta"+where,"#phi "+theName+" as a function of #eta",200,-2.5,2.5,nphi,-phi,phi);
00204   }
00205   
00206   HResolution(std::string name, TFile* file):theName(name.c_str()){ 
00207     //    dynamic_cast<TH1F*>( file->Get(theName+"") );
00208   }
00209   
00210   ~HResolution(){}
00211 
00212 
00213   void Fill(double p, double pt, double eta, double phi,
00214             double rp, double rpt,double reta, double rphi, double rcharge){
00215    
00216     Fill(rp, rpt, reta, rphi, rcharge);
00217 
00218     
00219     h2Eta->Fill(eta,reta);
00220     h2Phi->Fill(phi,rphi);
00221     
00222     h2P->Fill(p,rp); 
00223     h2Pt->Fill(pt,rpt);
00224 
00225     h2PtVsEta->Fill(eta,rpt);
00226     h2PtVsPhi->Fill(phi,rpt);
00227     
00228     h2EtaVsPt ->Fill(pt,reta);
00229     h2EtaVsPhi->Fill(phi,reta);
00230     
00231     h2PhiVsPt ->Fill(pt,rphi);
00232     h2PhiVsEta->Fill(eta,rphi);
00233   }
00234 
00235 
00236   void Fill(double p, double pt, double eta, double phi,
00237             double rp, double rpt){
00238    
00239     hP->Fill(rp); 
00240     hPt->Fill(rpt);
00241     
00242     h2P->Fill(p,rp); 
00243     // h2PVsEta->Fill(eta,rp);
00244     // h2PVsPhi->Fill(phi,rp);
00245 
00246     h2Pt->Fill(pt,rpt);
00247     h2PtVsEta->Fill(eta,rpt);
00248     h2PtVsPhi->Fill(phi,rpt);
00249   }
00250     
00251   void Fill(double rp, double rpt, 
00252             double reta, double rphi, double rcharge){
00253     
00254     hEta->Fill(reta);
00255     hPhi->Fill(rphi);
00256     
00257     hP->Fill(rp); 
00258     hPt->Fill(rpt);
00259     
00260     hCharge->Fill(rcharge);
00261   }
00262 
00263 
00264 
00265 private:
00266   DQMStore* dbe_;
00267 
00268   std::string theName;
00269   std::string where;
00270 
00271   MonitorElement *hEta;
00272   MonitorElement *hPhi;
00273 
00274   MonitorElement *hP; 
00275   MonitorElement *hPt; 
00276 
00277   MonitorElement *hCharge;
00278 
00279   MonitorElement *h2Eta;
00280   MonitorElement *h2Phi;
00281 
00282   MonitorElement *h2P; 
00283   MonitorElement *h2Pt; 
00284 
00285   MonitorElement *h2PtVsEta;
00286   MonitorElement *h2PtVsPhi;
00287 
00288   MonitorElement *h2EtaVsPt;
00289   MonitorElement *h2EtaVsPhi;
00290   
00291   MonitorElement *h2PhiVsPt; 
00292   MonitorElement *h2PhiVsEta;
00293 
00294 };
00295 
00296 
00297 class HResolution1DRecHit{
00298  public:
00299   HResolution1DRecHit(std::string name):theName(name.c_str()){
00300 
00301     // Position, sigma, residual, pull
00302     hResX        = dbe_->book1D (theName+"_X_Res", "X residual", 5000, -0.5,0.5);
00303     hResY        = dbe_->book1D (theName+"_Y_Res", "Y residual", 5000, -1.,1.);
00304     hResZ        = dbe_->book1D (theName+"_Z_Res", "Z residual", 5000, -1.5,1.5);
00305 
00306     hResXVsEta   = dbe_->book2D(theName+"_XResVsEta", "X residual vs eta",
00307                             200, -2.5,2.5, 5000, -1.5,1.5);
00308     hResYVsEta   = dbe_->book2D(theName+"_YResVsEta", "Y residual vs eta",
00309                             200, -2.5,2.5, 5000, -1.,1.);
00310     hResZVsEta   = dbe_->book2D(theName+"_ZResVsEta", "Z residual vs eta",
00311                             200, -2.5,2.5, 5000, -1.5,1.5);
00312     
00313     hXResVsPhi   = dbe_->book2D(theName+"_XResVsPhi", "X residual vs phi",
00314                             100,-Geom::pi(),Geom::pi(), 5000, -0.5,0.5);
00315     hYResVsPhi   = dbe_->book2D(theName+"_YResVsPhi", "Y residual vs phi",
00316                             100,-Geom::pi(),Geom::pi(), 5000, -1.,1.);
00317     hZResVsPhi   = dbe_->book2D(theName+"_ZResVsPhi", "Z residual vs phi",
00318                             100,-Geom::pi(),Geom::pi(), 5000, -1.5,1.5);
00319     
00320     hXResVsPos   = dbe_->book2D(theName+"_XResVsPos", "X residual vs position",
00321                             10000, -750,750, 5000, -0.5,0.5);    
00322     hYResVsPos   = dbe_->book2D(theName+"_YResVsPos", "Y residual vs position",
00323                             10000, -740,740, 5000, -1.,1.);    
00324     hZResVsPos   = dbe_->book2D(theName+"_ZResVsPos", "Z residual vs position",
00325                             10000, -1100,1100, 5000, -1.5,1.5);   
00326     
00327     hXPull       = dbe_->book1D (theName+"_XPull", "X pull", 600, -2,2);
00328     hYPull       = dbe_->book1D (theName+"_YPull", "Y pull", 600, -2,2);
00329     hZPull       = dbe_->book1D (theName+"_ZPull", "Z pull", 600, -2,2);
00330 
00331     hXPullVsPos  = dbe_->book2D (theName+"_XPullVsPos", "X pull vs position",10000, -750,750, 600, -2,2);
00332     hYPullVsPos  = dbe_->book2D (theName+"_YPullVsPos", "Y pull vs position",10000, -740,740, 600, -2,2);
00333     hZPullVsPos  = dbe_->book2D (theName+"_ZPullVsPos", "Z pull vs position",10000, -1100,1100, 600, -2,2);
00334   }
00335   
00336   HResolution1DRecHit(TString name_, TFile* file){}
00337 
00338   ~HResolution1DRecHit(){}
00339 
00340   void Fill(double x, double y, double z,
00341             double dx, double dy, double dz,
00342             double errx, double erry, double errz,
00343             double eta, double phi) {
00344     
00345     double rx = dx/x, ry = dy/y, rz = dz/z;
00346 
00347     hResX->Fill(rx);
00348     hResY->Fill(ry);
00349     hResZ->Fill(rz);
00350     
00351     hResXVsEta->Fill(eta,rx);
00352     hResYVsEta->Fill(eta,ry);
00353     hResZVsEta->Fill(eta,rz);
00354     
00355     hXResVsPhi->Fill(phi,rx);
00356     hYResVsPhi->Fill(phi,ry);
00357     hZResVsPhi->Fill(phi,rz);
00358     
00359     hXResVsPos->Fill(x,rx);
00360     hYResVsPos->Fill(y,ry);
00361     hZResVsPos->Fill(z,rz);
00362     
00363     if(errx < 1e-6)
00364       std::cout << "NO proper error set for X: "<<errx<<std::endl;
00365     else{
00366       hXPull->Fill(dx/errx);
00367       hXPullVsPos->Fill(x,dx/errx);
00368     }
00369     
00370     if(erry < 1e-6) 
00371       std::cout << "NO proper error set for Y: "<<erry<<std::endl;
00372     else{
00373       hYPull->Fill(dy/erry);
00374       hYPullVsPos->Fill(y,dy/erry);
00375     }
00376     if(errz < 1e-6)
00377       std::cout << "NO proper error set for Z: "<<errz<<std::endl; 
00378     else{
00379       hZPull->Fill(dz/errz);
00380       hZPullVsPos->Fill(z,dz/errz);
00381     }
00382   }
00383   
00384 
00385   
00386  public:
00387   std::string theName;
00388   
00389   MonitorElement *hResX;
00390   MonitorElement *hResY;
00391   MonitorElement *hResZ;
00392   
00393   MonitorElement *hResXVsEta;
00394   MonitorElement *hResYVsEta;
00395   MonitorElement *hResZVsEta;
00396   
00397   MonitorElement *hXResVsPhi;
00398   MonitorElement *hYResVsPhi;
00399   MonitorElement *hZResVsPhi;
00400   
00401   MonitorElement *hXResVsPos;
00402   MonitorElement *hYResVsPos;
00403   MonitorElement *hZResVsPos;
00404   
00405   MonitorElement *hXPull;
00406   MonitorElement *hYPull;
00407   MonitorElement *hZPull;
00408   
00409   MonitorElement *hXPullVsPos;
00410   MonitorElement *hYPullVsPos;
00411   MonitorElement *hZPullVsPos;
00412 
00413  private:
00414   DQMStore* dbe_;
00415 };
00416 #endif
00417 

Generated on Tue Jun 9 17:49:04 2009 for CMSSW by  doxygen 1.5.4