00001 #include "Validation/RecoVertex/interface/BSvsPVHistogramMaker.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/ServiceRegistry/interface/Service.h"
00006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00007 #include "DataFormats/VertexReco/interface/Vertex.h"
00008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00009 #include "TH1F.h"
00010 #include "TH2F.h"
00011 #include "TProfile.h"
00012
00013
00014 BSvsPVHistogramMaker::BSvsPVHistogramMaker():
00015 _currdir(0), m_maxLS(100), useSlope_(true), _trueOnly(true),
00016 _runHisto(true), _runHistoProfile(true), _runHistoBXProfile(true), _runHistoBX2D(false), _histoParameters() { }
00017
00018 BSvsPVHistogramMaker::BSvsPVHistogramMaker(const edm::ParameterSet& iConfig):
00019 _currdir(0),
00020 m_maxLS(iConfig.getParameter<unsigned int>("maxLSBeforeRebin")),
00021 useSlope_(iConfig.getParameter<bool>("useSlope")),
00022 _trueOnly(iConfig.getUntrackedParameter<bool>("trueOnly",true)),
00023 _runHisto(iConfig.getUntrackedParameter<bool>("runHisto",true)),
00024 _runHistoProfile(iConfig.getUntrackedParameter<bool>("runHistoProfile",true)),
00025 _runHistoBXProfile(iConfig.getUntrackedParameter<bool>("runHistoBXProfile",true)),
00026 _runHistoBX2D(iConfig.getUntrackedParameter<bool>("runHistoBX2D",false)),
00027 _histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters",edm::ParameterSet())),
00028 _rhm()
00029 { }
00030
00031
00032 BSvsPVHistogramMaker::~BSvsPVHistogramMaker() {
00033
00034 delete _currdir;
00035
00036 }
00037
00038
00039 void BSvsPVHistogramMaker::book(const std::string dirname) {
00040
00041 edm::Service<TFileService> tfserv;
00042 TFileDirectory* currdir = &(*tfserv);
00043
00044 if(dirname!="") {
00045 currdir = new TFileDirectory(tfserv->mkdir(dirname));
00046 _currdir = currdir;
00047 }
00048
00049 edm::LogInfo("HistogramBooking") << "Vertex histogram booking in directory " << dirname;
00050
00051 _hdeltax = currdir->make<TH1F>("deltax","(PV-BS) X position",
00052 _histoParameters.getUntrackedParameter<unsigned int>("nBinX",200),
00053 _histoParameters.getUntrackedParameter<double>("xMin",-1.),
00054 _histoParameters.getUntrackedParameter<double>("xMax",1.)
00055 );
00056 _hdeltax->GetXaxis()->SetTitle("#Delta(X) [cm]"); _hdeltax->GetYaxis()->SetTitle("Vertices");
00057
00058 _hdeltay = currdir->make<TH1F>("deltay","(PV-BS) Y position",
00059 _histoParameters.getUntrackedParameter<unsigned int>("nBinY",200),
00060 _histoParameters.getUntrackedParameter<double>("yMin",-1.),
00061 _histoParameters.getUntrackedParameter<double>("yMax",1.)
00062 );
00063 _hdeltay->GetXaxis()->SetTitle("#Delta(Y) [cm]"); _hdeltay->GetYaxis()->SetTitle("Vertices");
00064
00065 _hdeltaz = currdir->make<TH1F>("deltaz","(PV-BS) Z position",
00066 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ",200),
00067 _histoParameters.getUntrackedParameter<double>("zMin",-20.),
00068 _histoParameters.getUntrackedParameter<double>("zMax",20.)
00069 );
00070 _hdeltaz->GetXaxis()->SetTitle("#Delta(Z) [cm]"); _hdeltaz->GetYaxis()->SetTitle("Vertices");
00071
00072 _hdeltaxvsz = currdir->make<TProfile>("deltaxvsz","(PV-BS) X position vs Z",
00073 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile",40),
00074 _histoParameters.getUntrackedParameter<double>("zMinProfile",-20.),
00075 _histoParameters.getUntrackedParameter<double>("zMaxProfile",20.)
00076 );
00077 _hdeltaxvsz->GetXaxis()->SetTitle("Z [cm]"); _hdeltaxvsz->GetYaxis()->SetTitle("#Delta(X) [cm]");
00078
00079 _hdeltayvsz = currdir->make<TProfile>("deltayvsz","(PV-BS) Y position vs Z",
00080 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile",40),
00081 _histoParameters.getUntrackedParameter<double>("zMinProfile",-20.),
00082 _histoParameters.getUntrackedParameter<double>("zMaxProfile",20.)
00083 );
00084 _hdeltayvsz->GetXaxis()->SetTitle("Z [cm]"); _hdeltayvsz->GetYaxis()->SetTitle("#Delta(Y) [cm]");
00085
00086
00087
00088
00089 if(_runHisto) {
00090 _hdeltaxrun = _rhm.makeTH1F("deltaxrun","(PV-BS) X position",
00091 _histoParameters.getUntrackedParameter<unsigned int>("nBinX",200),
00092 _histoParameters.getUntrackedParameter<double>("xMin",-1.),
00093 _histoParameters.getUntrackedParameter<double>("xMax",1.));
00094
00095 _hdeltayrun = _rhm.makeTH1F("deltayrun","(PV-BS) Y position",
00096 _histoParameters.getUntrackedParameter<unsigned int>("nBinY",200),
00097 _histoParameters.getUntrackedParameter<double>("yMin",-1.),
00098 _histoParameters.getUntrackedParameter<double>("yMax",1.));
00099
00100 _hdeltazrun = _rhm.makeTH1F("deltazrun","(PV-BS) Z position",
00101 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ",200),
00102 _histoParameters.getUntrackedParameter<double>("zMin",-20.),
00103 _histoParameters.getUntrackedParameter<double>("zMax",20.));
00104
00105 _hdeltaxvszrun = _rhm.makeTProfile("deltaxvszrun","(PV-BS) X position vs Z",
00106 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile",40),
00107 _histoParameters.getUntrackedParameter<double>("zMinProfile",-20.),
00108 _histoParameters.getUntrackedParameter<double>("zMaxProfile",20.)
00109 );
00110
00111 _hdeltayvszrun = _rhm.makeTProfile("deltayvszrun","(PV-BS) Y position vs Z",
00112 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile",40),
00113 _histoParameters.getUntrackedParameter<double>("zMinProfile",-20.),
00114 _histoParameters.getUntrackedParameter<double>("zMaxProfile",20.)
00115 );
00116
00117 if(_runHistoProfile) {
00118 _hdeltaxvsorbrun = _rhm.makeTProfile("deltaxvsorbrun","(PV-BS) X position vs orbit number",4*m_maxLS,0.5,m_maxLS*262144+0.5);
00119 _hdeltayvsorbrun = _rhm.makeTProfile("deltayvsorbrun","(PV-BS) Y position vs orbit number",4*m_maxLS,0.5,m_maxLS*262144+0.5);
00120 _hdeltazvsorbrun = _rhm.makeTProfile("deltazvsorbrun","(PV-BS) Z position vs orbit number",4*m_maxLS,0.5,m_maxLS*262144+0.5);
00121 }
00122 if(_runHistoBXProfile) {
00123 _hdeltaxvsbxrun = _rhm.makeTProfile("deltaxvsbxrun","(PV-BS) X position vs BX number",3564,-0.5,3563.5);
00124 _hdeltayvsbxrun = _rhm.makeTProfile("deltayvsbxrun","(PV-BS) Y position vs BX number",3564,-0.5,3563.5);
00125 _hdeltazvsbxrun = _rhm.makeTProfile("deltazvsbxrun","(PV-BS) Z position vs BX number",3564,-0.5,3563.5);
00126 if(_runHistoBX2D) {
00127 _hdeltaxvsbx2drun = _rhm.makeTH2F("deltaxvsbx2drun","(PV-BS) X position vs BX number",3564,-0.5,3563.5,
00128 _histoParameters.getUntrackedParameter<unsigned int>("nBinX",200),
00129 _histoParameters.getUntrackedParameter<double>("xMin",-1.),
00130 _histoParameters.getUntrackedParameter<double>("xMax",1.));
00131 _hdeltayvsbx2drun = _rhm.makeTH2F("deltayvsbx2drun","(PV-BS) Y position vs BX number",3564,-0.5,3563.5,
00132 _histoParameters.getUntrackedParameter<unsigned int>("nBinY",200),
00133 _histoParameters.getUntrackedParameter<double>("yMin",-1.),
00134 _histoParameters.getUntrackedParameter<double>("yMax",1.));
00135 _hdeltazvsbx2drun = _rhm.makeTH2F("deltazvsbx2drun","(PV-BS) Z position vs BX number",3564,-0.5,3563.5,
00136 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ",200),
00137 _histoParameters.getUntrackedParameter<double>("zMin",-20.),
00138 _histoParameters.getUntrackedParameter<double>("zMax",20.));
00139 }
00140 }
00141
00142 }
00143 }
00144
00145 void BSvsPVHistogramMaker::beginRun(const unsigned int nrun) {
00146
00147 char runname[100];
00148 sprintf(runname,"run_%d",nrun);
00149
00150 TFileDirectory* currdir = _currdir;
00151 if(currdir==0) {
00152 edm::Service<TFileService> tfserv;
00153 currdir = &(*tfserv);
00154 }
00155
00156 _rhm.beginRun(nrun,*currdir);
00157
00158 if(_runHisto) {
00159 (*_hdeltaxrun)->GetXaxis()->SetTitle("#Delta(X) [cm]"); (*_hdeltaxrun)->GetYaxis()->SetTitle("Vertices");
00160 (*_hdeltayrun)->GetXaxis()->SetTitle("#Delta(Y) [cm]"); (*_hdeltayrun)->GetYaxis()->SetTitle("Vertices");
00161 (*_hdeltazrun)->GetXaxis()->SetTitle("#Delta(Z) [cm]"); (*_hdeltazrun)->GetYaxis()->SetTitle("Vertices");
00162 (*_hdeltaxvszrun)->GetXaxis()->SetTitle("Z [cm]"); (*_hdeltaxvszrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
00163 (*_hdeltayvszrun)->GetXaxis()->SetTitle("Z [cm]"); (*_hdeltayvszrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
00164
00165 if(_runHistoProfile) {
00166 (*_hdeltaxvsorbrun)->GetXaxis()->SetTitle("time [orbit#]"); (*_hdeltaxvsorbrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
00167 (*_hdeltaxvsorbrun)->SetBit(TH1::kCanRebin);
00168 (*_hdeltayvsorbrun)->GetXaxis()->SetTitle("time [orbit#]"); (*_hdeltayvsorbrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
00169 (*_hdeltayvsorbrun)->SetBit(TH1::kCanRebin);
00170 (*_hdeltazvsorbrun)->GetXaxis()->SetTitle("time [orbit#]"); (*_hdeltazvsorbrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
00171 (*_hdeltazvsorbrun)->SetBit(TH1::kCanRebin);
00172 }
00173 if(_runHistoBXProfile) {
00174 (*_hdeltaxvsbxrun)->GetXaxis()->SetTitle("BX"); (*_hdeltaxvsbxrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
00175 (*_hdeltayvsbxrun)->GetXaxis()->SetTitle("BX"); (*_hdeltayvsbxrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
00176 (*_hdeltazvsbxrun)->GetXaxis()->SetTitle("BX"); (*_hdeltazvsbxrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
00177 if(_runHistoBX2D) {
00178 (*_hdeltaxvsbx2drun)->GetXaxis()->SetTitle("BX"); (*_hdeltaxvsbx2drun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
00179 (*_hdeltayvsbx2drun)->GetXaxis()->SetTitle("BX"); (*_hdeltayvsbx2drun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
00180 (*_hdeltazvsbx2drun)->GetXaxis()->SetTitle("BX"); (*_hdeltazvsbx2drun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
00181 }
00182 }
00183
00184 }
00185 }
00186
00187 void BSvsPVHistogramMaker::fill(const unsigned int orbit, const int bx, const reco::VertexCollection& vertices, const reco::BeamSpot& bs) {
00188
00189 for(reco::VertexCollection::const_iterator vtx=vertices.begin();vtx!=vertices.end();++vtx) {
00190
00191 if(!(_trueOnly && vtx->isFake())) {
00192
00193
00194
00195
00196
00197
00198 double deltax = vtx->x()-x(bs,vtx->z());
00199 double deltay = vtx->y()-y(bs,vtx->z());
00200 double deltaz = vtx->z()-bs.z0();
00201
00202 _hdeltax->Fill(deltax);
00203 _hdeltay->Fill(deltay);
00204 _hdeltaz->Fill(deltaz);
00205 _hdeltaxvsz->Fill(vtx->z(),deltax);
00206 _hdeltayvsz->Fill(vtx->z(),deltay);
00207
00208 if(_runHisto) {
00209 if(_hdeltaxrun && *_hdeltaxrun ) (*_hdeltaxrun)->Fill(deltax);
00210 if(_hdeltayrun && *_hdeltayrun ) (*_hdeltayrun)->Fill(deltay);
00211 if(_hdeltazrun && *_hdeltazrun ) (*_hdeltazrun)->Fill(deltaz);
00212 if(_hdeltaxvszrun && *_hdeltaxvszrun ) (*_hdeltaxvszrun)->Fill(vtx->z(),deltax);
00213 if(_hdeltayvszrun && *_hdeltayvszrun ) (*_hdeltayvszrun)->Fill(vtx->z(),deltay);
00214 if(_runHistoProfile) {
00215 if(_hdeltaxvsorbrun && *_hdeltaxvsorbrun ) (*_hdeltaxvsorbrun)->Fill(orbit,deltax);
00216 if(_hdeltayvsorbrun && *_hdeltayvsorbrun ) (*_hdeltayvsorbrun)->Fill(orbit,deltay);
00217 if(_hdeltazvsorbrun && *_hdeltazvsorbrun ) (*_hdeltazvsorbrun)->Fill(orbit,deltaz);
00218 }
00219 if(_runHistoBXProfile) {
00220 if(_hdeltaxvsbxrun && *_hdeltaxvsbxrun ) (*_hdeltaxvsbxrun)->Fill(bx,deltax);
00221 if(_hdeltayvsbxrun && *_hdeltayvsbxrun ) (*_hdeltayvsbxrun)->Fill(bx,deltay);
00222 if(_hdeltazvsbxrun && *_hdeltazvsbxrun ) (*_hdeltazvsbxrun)->Fill(bx,deltaz);
00223 if(_runHistoBX2D) {
00224 if(_hdeltaxvsbx2drun && *_hdeltaxvsbx2drun ) (*_hdeltaxvsbx2drun)->Fill(bx,deltax);
00225 if(_hdeltayvsbx2drun && *_hdeltayvsbx2drun ) (*_hdeltayvsbx2drun)->Fill(bx,deltay);
00226 if(_hdeltazvsbx2drun && *_hdeltazvsbx2drun ) (*_hdeltazvsbx2drun)->Fill(bx,deltaz);
00227 }
00228 }
00229 }
00230 }
00231 }
00232 }
00233
00234 void BSvsPVHistogramMaker::fill(const edm::Event& iEvent, const reco::VertexCollection& vertices, const reco::BeamSpot& bs) {
00235
00236 fill(iEvent.orbitNumber(),iEvent.bunchCrossing(),vertices,bs);
00237
00238 }
00239
00240 double BSvsPVHistogramMaker::x(const reco::BeamSpot& bs, const double z) const {
00241
00242 double x = bs.x0();
00243
00244
00245 if(useSlope_) x += bs.dxdz()*(z-bs.z0());
00246
00247
00248
00249 return x;
00250
00251 }
00252
00253 double BSvsPVHistogramMaker::y(const reco::BeamSpot& bs, const double z) const {
00254
00255 double y = bs.y0();
00256
00257
00258 if(useSlope_) y += bs.dydz()*(z-bs.z0());
00259
00260
00261
00262 return y;
00263
00264 }
00265