Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "Validation/RecoVertex/interface/BSvsPVHistogramMaker.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0009 #include "TH1F.h"
0010 #include "TH2F.h"
0011 #include "TProfile.h"
0012 
0013 BSvsPVHistogramMaker::BSvsPVHistogramMaker(edm::ConsumesCollector&& iC)
0014     : _currdir(nullptr),
0015       m_maxLS(100),
0016       useSlope_(true),
0017       _trueOnly(true),
0018       _runHisto(true),
0019       _runHistoProfile(true),
0020       _runHistoBXProfile(true),
0021       _runHistoBX2D(false),
0022       _histoParameters(),
0023       _rhm(iC) {}
0024 
0025 BSvsPVHistogramMaker::BSvsPVHistogramMaker(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0026     : _currdir(nullptr),
0027       m_maxLS(iConfig.getParameter<unsigned int>("maxLSBeforeRebin")),
0028       useSlope_(iConfig.getParameter<bool>("useSlope")),
0029       _trueOnly(iConfig.getUntrackedParameter<bool>("trueOnly", true)),
0030       _runHisto(iConfig.getUntrackedParameter<bool>("runHisto", true)),
0031       _runHistoProfile(iConfig.getUntrackedParameter<bool>("runHistoProfile", true)),
0032       _runHistoBXProfile(iConfig.getUntrackedParameter<bool>("runHistoBXProfile", true)),
0033       _runHistoBX2D(iConfig.getUntrackedParameter<bool>("runHistoBX2D", false)),
0034       _histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters", edm::ParameterSet())),
0035       _rhm(iC) {}
0036 
0037 BSvsPVHistogramMaker::~BSvsPVHistogramMaker() { delete _currdir; }
0038 
0039 void BSvsPVHistogramMaker::book(const std::string dirname) {
0040   edm::Service<TFileService> tfserv;
0041   TFileDirectory* currdir = &(tfserv->tFileDirectory());
0042 
0043   if (!dirname.empty()) {
0044     currdir = new TFileDirectory(tfserv->mkdir(dirname));
0045     _currdir = currdir;
0046   }
0047 
0048   edm::LogInfo("HistogramBooking") << "Vertex histogram booking in directory " << dirname;
0049 
0050   _hdeltax = currdir->make<TH1F>("deltax",
0051                                  "(PV-BS) X position",
0052                                  _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0053                                  _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0054                                  _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0055   _hdeltax->GetXaxis()->SetTitle("#Delta(X) [cm]");
0056   _hdeltax->GetYaxis()->SetTitle("Vertices");
0057 
0058   _hdeltay = currdir->make<TH1F>("deltay",
0059                                  "(PV-BS) Y position",
0060                                  _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0061                                  _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0062                                  _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0063   _hdeltay->GetXaxis()->SetTitle("#Delta(Y) [cm]");
0064   _hdeltay->GetYaxis()->SetTitle("Vertices");
0065 
0066   _hdeltaz = currdir->make<TH1F>("deltaz",
0067                                  "(PV-BS) Z position",
0068                                  _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0069                                  _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0070                                  _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0071   _hdeltaz->GetXaxis()->SetTitle("#Delta(Z) [cm]");
0072   _hdeltaz->GetYaxis()->SetTitle("Vertices");
0073 
0074   _hdeltaxvsz = currdir->make<TProfile>("deltaxvsz",
0075                                         "(PV-BS) X position vs Z",
0076                                         _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0077                                         _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0078                                         _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0079   _hdeltaxvsz->GetXaxis()->SetTitle("Z [cm]");
0080   _hdeltaxvsz->GetYaxis()->SetTitle("#Delta(X) [cm]");
0081 
0082   _hdeltayvsz = currdir->make<TProfile>("deltayvsz",
0083                                         "(PV-BS) Y position vs Z",
0084                                         _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0085                                         _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0086                                         _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0087   _hdeltayvsz->GetXaxis()->SetTitle("Z [cm]");
0088   _hdeltayvsz->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0089 
0090   if (_runHisto) {
0091     _hdeltaxrun = _rhm.makeTH1F("deltaxrun",
0092                                 "(PV-BS) X position",
0093                                 _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0094                                 _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0095                                 _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0096 
0097     _hdeltayrun = _rhm.makeTH1F("deltayrun",
0098                                 "(PV-BS) Y position",
0099                                 _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0100                                 _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0101                                 _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0102 
0103     _hdeltazrun = _rhm.makeTH1F("deltazrun",
0104                                 "(PV-BS) Z position",
0105                                 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0106                                 _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0107                                 _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0108 
0109     _hdeltaxvszrun = _rhm.makeTProfile("deltaxvszrun",
0110                                        "(PV-BS) X position vs Z",
0111                                        _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0112                                        _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0113                                        _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0114 
0115     _hdeltayvszrun = _rhm.makeTProfile("deltayvszrun",
0116                                        "(PV-BS) Y position vs Z",
0117                                        _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0118                                        _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0119                                        _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0120 
0121     if (_runHistoProfile) {
0122       _hdeltaxvsorbrun = _rhm.makeTProfile(
0123           "deltaxvsorbrun", "(PV-BS) X position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0124       _hdeltayvsorbrun = _rhm.makeTProfile(
0125           "deltayvsorbrun", "(PV-BS) Y position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0126       _hdeltazvsorbrun = _rhm.makeTProfile(
0127           "deltazvsorbrun", "(PV-BS) Z position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0128     }
0129     if (_runHistoBXProfile) {
0130       _hdeltaxvsbxrun = _rhm.makeTProfile("deltaxvsbxrun", "(PV-BS) X position vs BX number", 3564, -0.5, 3563.5);
0131       _hdeltayvsbxrun = _rhm.makeTProfile("deltayvsbxrun", "(PV-BS) Y position vs BX number", 3564, -0.5, 3563.5);
0132       _hdeltazvsbxrun = _rhm.makeTProfile("deltazvsbxrun", "(PV-BS) Z position vs BX number", 3564, -0.5, 3563.5);
0133       if (_runHistoBX2D) {
0134         _hdeltaxvsbx2drun = _rhm.makeTH2F("deltaxvsbx2drun",
0135                                           "(PV-BS) X position vs BX number",
0136                                           3564,
0137                                           -0.5,
0138                                           3563.5,
0139                                           _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0140                                           _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0141                                           _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0142         _hdeltayvsbx2drun = _rhm.makeTH2F("deltayvsbx2drun",
0143                                           "(PV-BS) Y position vs BX number",
0144                                           3564,
0145                                           -0.5,
0146                                           3563.5,
0147                                           _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0148                                           _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0149                                           _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0150         _hdeltazvsbx2drun = _rhm.makeTH2F("deltazvsbx2drun",
0151                                           "(PV-BS) Z position vs BX number",
0152                                           3564,
0153                                           -0.5,
0154                                           3563.5,
0155                                           _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0156                                           _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0157                                           _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0158       }
0159     }
0160   }
0161 }
0162 
0163 void BSvsPVHistogramMaker::beginRun(const unsigned int nrun) {
0164   char runname[100];
0165   sprintf(runname, "run_%d", nrun);
0166 
0167   TFileDirectory* currdir = _currdir;
0168   if (currdir == nullptr) {
0169     edm::Service<TFileService> tfserv;
0170     currdir = &(tfserv->tFileDirectory());
0171   }
0172 
0173   _rhm.beginRun(nrun, *currdir);
0174 
0175   if (_runHisto) {
0176     (*_hdeltaxrun)->GetXaxis()->SetTitle("#Delta(X) [cm]");
0177     (*_hdeltaxrun)->GetYaxis()->SetTitle("Vertices");
0178     (*_hdeltayrun)->GetXaxis()->SetTitle("#Delta(Y) [cm]");
0179     (*_hdeltayrun)->GetYaxis()->SetTitle("Vertices");
0180     (*_hdeltazrun)->GetXaxis()->SetTitle("#Delta(Z) [cm]");
0181     (*_hdeltazrun)->GetYaxis()->SetTitle("Vertices");
0182     (*_hdeltaxvszrun)->GetXaxis()->SetTitle("Z [cm]");
0183     (*_hdeltaxvszrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0184     (*_hdeltayvszrun)->GetXaxis()->SetTitle("Z [cm]");
0185     (*_hdeltayvszrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0186 
0187     if (_runHistoProfile) {
0188       (*_hdeltaxvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0189       (*_hdeltaxvsorbrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0190       (*_hdeltaxvsorbrun)->SetCanExtend(TH1::kAllAxes);
0191       (*_hdeltayvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0192       (*_hdeltayvsorbrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0193       (*_hdeltayvsorbrun)->SetCanExtend(TH1::kAllAxes);
0194       (*_hdeltazvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0195       (*_hdeltazvsorbrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0196       (*_hdeltazvsorbrun)->SetCanExtend(TH1::kAllAxes);
0197     }
0198     if (_runHistoBXProfile) {
0199       (*_hdeltaxvsbxrun)->GetXaxis()->SetTitle("BX");
0200       (*_hdeltaxvsbxrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0201       (*_hdeltayvsbxrun)->GetXaxis()->SetTitle("BX");
0202       (*_hdeltayvsbxrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0203       (*_hdeltazvsbxrun)->GetXaxis()->SetTitle("BX");
0204       (*_hdeltazvsbxrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0205       if (_runHistoBX2D) {
0206         (*_hdeltaxvsbx2drun)->GetXaxis()->SetTitle("BX");
0207         (*_hdeltaxvsbx2drun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0208         (*_hdeltayvsbx2drun)->GetXaxis()->SetTitle("BX");
0209         (*_hdeltayvsbx2drun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0210         (*_hdeltazvsbx2drun)->GetXaxis()->SetTitle("BX");
0211         (*_hdeltazvsbx2drun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0212       }
0213     }
0214   }
0215 }
0216 
0217 void BSvsPVHistogramMaker::fill(const unsigned int orbit,
0218                                 const int bx,
0219                                 const reco::VertexCollection& vertices,
0220                                 const reco::BeamSpot& bs) {
0221   for (reco::VertexCollection::const_iterator vtx = vertices.begin(); vtx != vertices.end(); ++vtx) {
0222     if (!(_trueOnly && vtx->isFake())) {
0223       /*
0224       double deltax = vtx->x()-bs.x0();
0225       double deltay = vtx->y()-bs.y0();
0226       double deltaz = vtx->z()-bs.z0();
0227       */
0228       double deltax = vtx->x() - x(bs, vtx->z());
0229       double deltay = vtx->y() - y(bs, vtx->z());
0230       double deltaz = vtx->z() - bs.z0();
0231 
0232       _hdeltax->Fill(deltax);
0233       _hdeltay->Fill(deltay);
0234       _hdeltaz->Fill(deltaz);
0235       _hdeltaxvsz->Fill(vtx->z(), deltax);
0236       _hdeltayvsz->Fill(vtx->z(), deltay);
0237 
0238       if (_runHisto) {
0239         if (_hdeltaxrun && *_hdeltaxrun)
0240           (*_hdeltaxrun)->Fill(deltax);
0241         if (_hdeltayrun && *_hdeltayrun)
0242           (*_hdeltayrun)->Fill(deltay);
0243         if (_hdeltazrun && *_hdeltazrun)
0244           (*_hdeltazrun)->Fill(deltaz);
0245         if (_hdeltaxvszrun && *_hdeltaxvszrun)
0246           (*_hdeltaxvszrun)->Fill(vtx->z(), deltax);
0247         if (_hdeltayvszrun && *_hdeltayvszrun)
0248           (*_hdeltayvszrun)->Fill(vtx->z(), deltay);
0249         if (_runHistoProfile) {
0250           if (_hdeltaxvsorbrun && *_hdeltaxvsorbrun)
0251             (*_hdeltaxvsorbrun)->Fill(orbit, deltax);
0252           if (_hdeltayvsorbrun && *_hdeltayvsorbrun)
0253             (*_hdeltayvsorbrun)->Fill(orbit, deltay);
0254           if (_hdeltazvsorbrun && *_hdeltazvsorbrun)
0255             (*_hdeltazvsorbrun)->Fill(orbit, deltaz);
0256         }
0257         if (_runHistoBXProfile) {
0258           if (_hdeltaxvsbxrun && *_hdeltaxvsbxrun)
0259             (*_hdeltaxvsbxrun)->Fill(bx % 3564, deltax);
0260           if (_hdeltayvsbxrun && *_hdeltayvsbxrun)
0261             (*_hdeltayvsbxrun)->Fill(bx % 3564, deltay);
0262           if (_hdeltazvsbxrun && *_hdeltazvsbxrun)
0263             (*_hdeltazvsbxrun)->Fill(bx % 3564, deltaz);
0264           if (_runHistoBX2D) {
0265             if (_hdeltaxvsbx2drun && *_hdeltaxvsbx2drun)
0266               (*_hdeltaxvsbx2drun)->Fill(bx % 3564, deltax);
0267             if (_hdeltayvsbx2drun && *_hdeltayvsbx2drun)
0268               (*_hdeltayvsbx2drun)->Fill(bx % 3564, deltay);
0269             if (_hdeltazvsbx2drun && *_hdeltazvsbx2drun)
0270               (*_hdeltazvsbx2drun)->Fill(bx % 3564, deltaz);
0271           }
0272         }
0273       }
0274     }
0275   }
0276 }
0277 
0278 void BSvsPVHistogramMaker::fill(const edm::Event& iEvent,
0279                                 const reco::VertexCollection& vertices,
0280                                 const reco::BeamSpot& bs) {
0281   fill(iEvent.orbitNumber(), iEvent.bunchCrossing(), vertices, bs);
0282 }
0283 
0284 double BSvsPVHistogramMaker::x(const reco::BeamSpot& bs, const double z) const {
0285   double x = bs.x0();
0286 
0287   //  if(useSlope_) x += bs.dxdz()*z;
0288   if (useSlope_)
0289     x += bs.dxdz() * (z - bs.z0());
0290 
0291   //  if(useSlope_) x = bs.x(z);
0292 
0293   return x;
0294 }
0295 
0296 double BSvsPVHistogramMaker::y(const reco::BeamSpot& bs, const double z) const {
0297   double y = bs.y0();
0298 
0299   //  if(useSlope_) y += bs.dydz()*z;
0300   if (useSlope_)
0301     y += bs.dydz() * (z - bs.z0());
0302 
0303   //  if(useSlope_) y = bs.y(z);
0304 
0305   return y;
0306 }