CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PrimaryVertexAnalyzer4PU.cc
Go to the documentation of this file.
1 #include "Validation/RecoVertex/interface/PrimaryVertexAnalyzer4PU.h"
2 
3 #include "FWCore/ParameterSet/interface/ParameterSet.h"
4 #include "FWCore/MessageLogger/interface/MessageLogger.h"
5 #include "FWCore/Version/interface/GetReleaseVersion.h"
6 #include "MagneticField/Engine/interface/MagneticField.h"
7 
8 // reco track and vertex
9 #include "DataFormats/TrackReco/interface/Track.h"
10 #include "DataFormats/VertexReco/interface/Vertex.h"
11 #include "DataFormats/Math/interface/Point3D.h"
12 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
13 
14 // simulated vertices,..., add <use name=SimDataFormats/Vertex> and <../Track>
15 #include <SimDataFormats/Vertex/interface/SimVertex.h>
16 #include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
17 #include <SimDataFormats/Track/interface/SimTrack.h>
18 #include <SimDataFormats/Track/interface/SimTrackContainer.h>
19 
20 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
21 #include "SimDataFormats/CrossingFrame/interface/CrossingFramePlaybackInfo.h"
22 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
23 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
24 
25 //generator level + CLHEP
26 #include "HepMC/GenEvent.h"
27 #include "HepMC/GenVertex.h"
28 
29 
30 // TrackingParticle
31 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
32 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
33 //associator
34 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
35 
36 
37 // fit
38 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
39 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
40 
41 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
42 
43 // Root
44 #include <TH1.h>
45 #include <TH2.h>
46 #include <TFile.h>
47 #include <TProfile.h>
48 
49 #include <cmath>
50 #include <gsl/gsl_math.h>
51 #include <gsl/gsl_eigen.h>
52 
53 
54 // cluster stufff
55 //#include "DataFormats/TrackRecoTrack.h"
56 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
57 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
58 
59 
60 using namespace edm;
61 using namespace reco;
62 using namespace std;
63 //
64 // constants, enums and typedefs
65 //
66 typedef reco::Vertex::trackRef_iterator trackit_t;
67 //
68 // static data member definitions
69 //
70 
71 //
72 // constructors and destructor
73 //
74 PrimaryVertexAnalyzer4PU::PrimaryVertexAnalyzer4PU(const ParameterSet& iConfig):theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"))
75 {
76  //now do what ever initialization is needed
77  simG4_=iConfig.getParameter<edm::InputTag>( "simG4" );
78  recoTrackProducer_= iConfig.getUntrackedParameter<std::string>("recoTrackProducer");
79  // open output file to store histograms}
80  outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
81 
82  rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE");
83  verbose_= iConfig.getUntrackedParameter<bool>("verbose", false);
84  doMatching_= iConfig.getUntrackedParameter<bool>("matching", false);
85  printXBS_= iConfig.getUntrackedParameter<bool>("XBS", false);
86  simUnit_= 1.0; // starting with CMSSW_1_2_x ??
87  if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
88  simUnit_=0.1; // for use in CMSSW_1_1_1 tutorial
89  }
90 
91  dumpPUcandidates_=iConfig.getUntrackedParameter<bool>("dumpPUcandidates", false);
92 
93  zmatch_=iConfig.getUntrackedParameter<double>("zmatch", 0.0500);
94  cout << "PrimaryVertexAnalyzer4PU: zmatch=" << zmatch_ << endl;
95  eventcounter_=0;
96  dumpcounter_=0;
97  ndump_=10;
98  DEBUG_=false;
99  //DEBUG_=true;
100 }
101 
102 
103 
104 
105 PrimaryVertexAnalyzer4PU::~PrimaryVertexAnalyzer4PU()
106 {
107  // do anything here that needs to be done at desctruction time
108  // (e.g. close files, deallocate resources etc.)
109  delete rootFile_;
110 }
111 
112 
113 
114 //
115 // member functions
116 //
117 
118 
119 std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::bookVertexHistograms(){
120  std::map<std::string, TH1*> h;
121 
122  // release validation histograms used in DoCompare.C
123  h["nbtksinvtx"] = new TH1F("nbtksinvtx","reconstructed tracks in vertex",40,-0.5,39.5);
124  h["nbtksinvtxPU"] = new TH1F("nbtksinvtxPU","reconstructed tracks in vertex",40,-0.5,39.5);
125  h["nbtksinvtxTag"]= new TH1F("nbtksinvtxTag","reconstructed tracks in vertex",40,-0.5,39.5);
126  h["resx"] = new TH1F("resx","residual x",100,-0.04,0.04);
127  h["resy"] = new TH1F("resy","residual y",100,-0.04,0.04);
128  h["resz"] = new TH1F("resz","residual z",100,-0.1,0.1);
129  h["resz10"] = new TH1F("resz10","residual z",100,-1.0,1.);
130  h["pullx"] = new TH1F("pullx","pull x",100,-25.,25.);
131  h["pully"] = new TH1F("pully","pull y",100,-25.,25.);
132  h["pullz"] = new TH1F("pullz","pull z",100,-25.,25.);
133  h["vtxchi2"] = new TH1F("vtxchi2","chi squared",100,0.,100.);
134  h["vtxndf"] = new TH1F("vtxndf","degrees of freedom",500,0.,100.);
135  h["vtxndfc"] = new TH1F("vtxndfc","expected 2nd highest ndof",500,0.,100.);
136 
137  h["vtxndfvsntk"] = new TH2F("vtxndfvsntk","ndof vs #tracks",20,0.,100, 20, 0., 200.);
138  h["vtxndfoverntk"]= new TH1F("vtxndfoverntk","ndof / #tracks",40,0.,2.);
139  h["vtxndf2overntk"]= new TH1F("vtxndf2overntk","(ndof+2) / #tracks",40,0.,2.);
140  h["tklinks"] = new TH1F("tklinks","Usable track links",2,-0.5,1.5);
141  h["nans"] = new TH1F("nans","Illegal values for x,y,z,xx,xy,xz,yy,yz,zz",9,0.5,9.5);
142 
143 
144  // raw
145  add(h, new TH1F("szRecVtx","size of recvtx collection",20, -0.5, 19.5));
146  add(h, new TH1F("isFake","fake vertex",2, -0.5, 1.5));
147  add(h, new TH1F("isFake1","fake vertex or ndof<0",2, -0.5, 1.5));
148  add(h, new TH1F("bunchCrossing","bunchCrossing",4000, 0., 4000.));
149  add(h, new TH2F("bunchCrossingLogNtk","bunchCrossingLogNtk",4000, 0., 4000.,5,0., 5.));
150  add(h, new TH1F("highpurityTrackFraction","fraction of high purity tracks",20, 0., 1.));
151  add(h, new TH2F("trkchi2vsndof","vertices chi2 vs ndof",50, 0., 100., 50, 0., 200.));
152  add(h, new TH1F("trkchi2overndof","vertices chi2 / ndof",50, 0., 5.));
153  // two track vertices
154  add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
155  add(h,new TH1F("2trkmassSS","two-track vertices mass (same sign)",100, 0., 2.));
156  add(h,new TH1F("2trkmassOS","two-track vertices mass (opposite sign)",100, 0., 2.));
157  add(h,new TH1F("2trkdphi","two-track vertices delta-phi",360, 0, 2*M_PI));
158  add(h,new TH1F("2trkseta","two-track vertices sum-eta",50, -2., 2.));
159  add(h,new TH1F("2trkdphicurl","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
160  add(h,new TH1F("2trksetacurl","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
161  add(h,new TH1F("2trkdetaOS","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
162  add(h,new TH1F("2trkdetaSS","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
163  // two track PU vertices
164  add(h,new TH1F("2trkmassSSPU","two-track vertices mass (same sign)",100, 0., 2.));
165  add(h,new TH1F("2trkmassOSPU","two-track vertices mass (opposite sign)",100, 0., 2.));
166  add(h,new TH1F("2trkdphiPU","two-track vertices delta-phi",360, 0, 2*M_PI));
167  add(h,new TH1F("2trksetaPU","two-track vertices sum-eta",50, -2., 2.));
168  add(h,new TH1F("2trkdphicurlPU","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
169  add(h,new TH1F("2trksetacurlPU","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
170  add(h,new TH1F("2trkdetaOSPU","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
171  add(h,new TH1F("2trkdetaSSPU","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
172  // three track vertices
173  add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
174  add(h,new TH2F("3trkchi2vsndof","three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
175  add(h,new TH2F("4trkchi2vsndof","four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
176  add(h,new TH2F("5trkchi2vsndof","five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
177  // same for fakes
178  add(h,new TH2F("fake2trkchi2vsndof","fake two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
179  add(h,new TH2F("fake3trkchi2vsndof","fake three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
180  add(h,new TH2F("fake4trkchi2vsndof","fake four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
181  add(h,new TH2F("fake5trkchi2vsndof","fake five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
182  h["resxr"] = new TH1F("resxr","relative residual x",100,-0.04,0.04);
183  h["resyr"] = new TH1F("resyr","relative residual y",100,-0.04,0.04);
184  h["reszr"] = new TH1F("reszr","relative residual z",100,-0.1,0.1);
185  h["pullxr"] = new TH1F("pullxr","relative pull x",100,-25.,25.);
186  h["pullyr"] = new TH1F("pullyr","relative pull y",100,-25.,25.);
187  h["pullzr"] = new TH1F("pullzr","relative pull z",100,-25.,25.);
188  h["vtxprob"] = new TH1F("vtxprob","chisquared probability",100,0.,1.);
189  h["eff"] = new TH1F("eff","efficiency",2, -0.5, 1.5);
190  h["efftag"] = new TH1F("efftag","efficiency tagged vertex",2, -0.5, 1.5);
191  h["zdistancetag"] = new TH1F("zdistancetag","z-distance between tagged and generated",100, -0.1, 0.1);
192  h["abszdistancetag"] = new TH1F("abszdistancetag","z-distance between tagged and generated",1000, 0., 1.0);
193  h["abszdistancetagcum"] = new TH1F("abszdistancetagcum","z-distance between tagged and generated",1000, 0., 1.0);
194  h["puritytag"] = new TH1F("puritytag","purity of primary vertex tags",2, -0.5, 1.5);
195  h["effvsptsq"] = new TProfile("effvsptsq","efficiency vs ptsq",20, 0., 10000., 0, 1.);
196  h["effvsnsimtrk"] = new TProfile("effvsnsimtrk","efficiency vs # simtracks",50, 0., 50., 0, 1.);
197  h["effvsnrectrk"] = new TProfile("effvsnrectrk","efficiency vs # rectracks",50, 0., 50., 0, 1.);
198  h["effvsnseltrk"] = new TProfile("effvsnseltrk","efficiency vs # selected tracks",50, 0., 50., 0, 1.);
199  h["effvsz"] = new TProfile("effvsz","efficiency vs z",20, -20., 20., 0, 1.);
200  h["effvsz2"] = new TProfile("effvsz2","efficiency vs z (2mm)",20, -20., 20., 0, 1.);
201  h["effvsr"] = new TProfile("effvsr","efficiency vs r",20, 0., 1., 0, 1.);
202  h["xresvsntrk"] = new TProfile("xresvsntrk","xresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
203  h["yresvsntrk"] = new TProfile("yresvsntrk","yresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
204  h["zresvsntrk"] = new TProfile("zresvsntrk","zresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
205  h["cpuvsntrk"] = new TProfile("cpuvsntrk","cpu time vs # of fitted tracks",40, 0., 200., 0, 200.);
206  h["cpucluvsntrk"] = new TProfile("cpucluvsntrk","clustering cpu time # of tracks",40, 0., 200., 0, 10.);
207  h["cpufit"] = new TH1F("cpufit","cpu time for fitting",100, 0., 200.);
208  h["cpuclu"] = new TH1F("cpuclu","cpu time for clustering",100, 0., 200.);
209  h["nbtksinvtx2"] = new TH1F("nbtksinvtx2","reconstructed tracks in vertex",40,0.,200.);
210  h["nbtksinvtxPU2"] = new TH1F("nbtksinvtxPU2","reconstructed tracks in vertex",40,0.,200.);
211  h["nbtksinvtxTag2"] = new TH1F("nbtksinvtxTag2","reconstructed tracks in vertex",40,0.,200.);
212 
213  h["xrec"] = new TH1F("xrec","reconstructed x",100,-0.1,0.1);
214  h["yrec"] = new TH1F("yrec","reconstructed y",100,-0.1,0.1);
215  h["zrec"] = new TH1F("zrec","reconstructed z",100,-20.,20.);
216  h["err1"] = new TH1F("err1","error 1",100,0.,0.1);
217  h["err2"] = new TH1F("err2","error 2",100,0.,0.1);
218  h["errx"] = new TH1F("errx","error x",100,0.,0.1);
219  h["erry"] = new TH1F("erry","error y",100,0.,0.1);
220  h["errz"] = new TH1F("errz","error z",100,0.,2.0);
221  h["errz1"] = new TH1F("errz1","error z",100,0.,0.2);
222 
223  h["zrecNt100"] = new TH1F("zrecNt100","reconstructed z for high multiplicity events",80,-40.,40.);
224  add(h, new TH2F("zrecvsnt","reconstructed z vs number of tracks",100,-50., 50., 20, 0., 200.));
225  add(h, new TH2F("xyrec","reconstructed xy",100, -4., 4., 100, -4., 4.));
226  h["xrecBeam"] = new TH1F("xrecBeam","reconstructed x - beam x",100,-0.1,0.1);
227  h["yrecBeam"] = new TH1F("yrecBeam","reconstructed y - beam y",100,-0.1,0.1);
228  h["zrecBeam"] = new TH1F("zrecBeam","reconstructed z - beam z",100,-20.,20.);
229  h["xrecBeamvsz"] = new TH2F("xrecBeamvsz","reconstructed x - beam x vs z", 20, -20., 20.,100,-0.1,0.1);
230  h["yrecBeamvsz"] = new TH2F("yrecBeamvsz","reconstructed y - beam y vs z", 20, -20., 20.,100,-0.1,0.1);
231  h["xrecBeamvszprof"] = new TProfile("xrecBeamvszprof","reconstructed x - beam x vs z-z0", 20, -20., 20.,-0.1,0.1);
232  h["yrecBeamvszprof"] = new TProfile("yrecBeamvszprof","reconstructed y - beam y vs z-z0", 20, -20., 20.,-0.1,0.1);
233  add(h, new TH2F("xrecBeamvsdxXBS","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1)); // just a test
234  add(h, new TH2F("yrecBeamvsdyXBS","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
235  add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1));
236  add(h, new TH2F("yrecBeamvsdy","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
237  add(h, new TH2F("xrecBeamvsdxR2","reconstructed x - beam x vs resolution",20,0., 0.04, 100, -0.1,0.1));
238  add(h, new TH2F("yrecBeamvsdyR2","reconstructed z - beam z vs resolution",20,0., 0.04, 100, -0.1,0.1));
239  // add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",100,-0.1,0.1, 10, 0., 0.04));
240  // add(h, new TH2F("yrecBeamvsdy","reconstructed y - beam y vs resolution",100,-0.1,0.1, 10, 0., 0.04));
241  h["xrecBeamvsdxprof"] = new TProfile("xrecBeamvsdxprof","reconstructed x - beam x vs resolution",10, 0., 0.04,-0.1,0.1 );
242  h["yrecBeamvsdyprof"] = new TProfile("yrecBeamvsdyprof","reconstructed y - beam y vs resolution",10, 0., 0.04,-0.1,0.1 );
243  add(h, new TProfile("xrecBeam2vsdx2prof","reconstructed x - beam x vs resolution",10,0., 0.002, 0., 0.01));
244  add(h, new TProfile("yrecBeam2vsdy2prof","reconstructed y - beam y vs resolution",10,0., 0.002, 0., 0.01));
245  add(h,new TH2F("xrecBeamvsdx2","reconstructed x - beam x vs resolution",10,0., 0.002, 100, -0.01, 0.01));
246  add(h,new TH2F("yrecBeamvsdy2","reconstructed y - beam y vs resolution",10,0., 0.002, 100, -0.01, 0.01));
247  h["xrecb"] = new TH1F("xrecb","reconstructed x - beam x",100,-0.01,0.01);
248  h["yrecb"] = new TH1F("yrecb","reconstructed y - beam y",100,-0.01,0.01);
249  h["zrecb"] = new TH1F("zrecb","reconstructed z - beam z",100,-20.,20.);
250  h["xrec1"] = new TH1F("xrec1","reconstructed x",100,-4,4);
251  h["yrec1"] = new TH1F("yrec1","reconstructed y",100,-4,4); // should match the sim histos
252  h["zrec1"] = new TH1F("zrec1","reconstructed z",100,-80.,80.);
253  h["xrec2"] = new TH1F("xrec2","reconstructed x",100,-1,1);
254  h["yrec2"] = new TH1F("yrec2","reconstructed y",100,-1,1);
255  h["zrec2"] = new TH1F("zrec2","reconstructed z",100,-40.,40.);
256  h["xrec3"] = new TH1F("xrec3","reconstructed x",100,-0.1,0.1);
257  h["yrec3"] = new TH1F("yrec3","reconstructed y",100,-0.1,0.1);
258  h["zrec3"] = new TH1F("zrec3","reconstructed z",100,-20.,20.);
259  add(h, new TH1F("xrecBeamPull","normalized residuals reconstructed x - beam x",100,-20,20));
260  add(h, new TH1F("yrecBeamPull","normalized residuals reconstructed y - beam y",100,-20,20));
261  add(h, new TH1F("zdiffrec","z-distance between vertices",200,-20., 20.));
262  add(h, new TH1F("zdiffrec2","z-distance between ndof>2 vertices",200,-20., 20.));
263  add(h, new TH1F("zdiffrec4","z-distance between ndof>4 vertices",200,-20., 20.));
264  add(h, new TH1F("zdiffrec8","z-distance between ndof>8 vertices",200,-20., 20.));
265  add(h, new TH2F("zvszrec2","z positions of multiple vertices",200,-20., 20., 200,-20., 20.));
266  add(h, new TH2F("pzvspz2","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
267  add(h, new TH2F("zvszrec4","z positions of multiple vertices",100,-20., 20., 100,-20., 20.));
268  add(h, new TH2F("pzvspz4","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
269  add(h, new TH2F("zdiffvsz","z-distance vs z",100,-10.,10.,30,-15.,15.));
270  add(h, new TH2F("zdiffvsz4","z-distance vs z (ndof>4)",100,-10.,10.,30,-15.,15.));
271  add(h, new TProfile("eff0vsntrec","efficiency vs # reconstructed tracks",50, 0., 50., 0, 1.));
272  add(h, new TProfile("eff0vsntsel","efficiency vs # selected tracks",50, 0., 50., 0, 1.));
273  add(h, new TProfile("eff0ndof0vsntsel","efficiency (ndof>0) vs # selected tracks",50, 0., 50., 0, 1.));
274  add(h, new TProfile("eff0ndof8vsntsel","efficiency (ndof>8) vs # selected tracks",50, 0., 50., 0, 1.));
275  add(h, new TProfile("eff0ndof2vsntsel","efficiency (ndof>2) vs # selected tracks",50, 0., 50., 0, 1.));
276  add(h, new TProfile("eff0ndof4vsntsel","efficiency (ndof>4) vs # selected tracks",50, 0., 50., 0, 1.));
277  add(h, new TH1F("xbeamPUcand","x - beam of PU candidates (ndof>4)",100,-1., 1.));
278  add(h, new TH1F("ybeamPUcand","y - beam of PU candidates (ndof>4)",100,-1., 1.));
279  add(h, new TH1F("xbeamPullPUcand","x - beam pull of PU candidates (ndof>4)",100,-20., 20.));
280  add(h, new TH1F("ybeamPullPUcand","y - beam pull of PU candidates (ndof>4)",100,-20., 20.));
281  add(h, new TH1F("ndofOverNtk","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
282  //add(h, new TH1F("sumwOverNtk","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
283  add(h, new TH1F("ndofOverNtkPUcand","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
284  //add(h, new TH1F("sumwOverNtkPUcand","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
285  add(h, new TH1F("zPUcand","z positions of PU candidates (all)",100,-20., 20.));
286  add(h, new TH1F("zPUcand2","z positions of PU candidates (ndof>2)",100,-20., 20.));
287  add(h, new TH1F("zPUcand4","z positions of PU candidates (ndof>4)",100,-20., 20.));
288  add(h, new TH1F("ndofPUcand","ndof of PU candidates (all)",50,0., 100.));
289  add(h, new TH1F("ndofPUcand2","ndof of PU candidates (ndof>2)",50,0., 100.));
290  add(h, new TH1F("ndofPUcand4","ndof of PU candidates (ndof>4)",50,0., 100.));
291  add(h, new TH1F("ndofnr2","second highest ndof",500,0., 100.));
292  add(h, new TH1F("ndofnr2d1cm","second highest ndof (dz>1cm)",500,0., 100.));
293  add(h, new TH1F("ndofnr2d2cm","second highest ndof (dz>2cm)",500,0., 100.));
294 
295  h["nrecvtx"] = new TH1F("nrecvtx","# of reconstructed vertices", 50, -0.5, 49.5);
296  h["nrecvtx8"] = new TH1F("nrecvtx8","# of reconstructed vertices with ndof>8", 50, -0.5, 49.5);
297  h["nrecvtx2"] = new TH1F("nrecvtx2","# of reconstructed vertices with ndof>2", 50, -0.5, 49.5);
298  h["nrecvtx4"] = new TH1F("nrecvtx4","# of reconstructed vertices with ndof>4", 50, -0.5, 49.5);
299  // h["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
300  h["nrectrk"] = new TH1F("nrectrk","# of reconstructed tracks", 100, -0.5, 99.5);
301  add(h, new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5));
302  add(h, new TH1F("nsimtrkSignal","# of simulated tracks (Signal)", 100, -0.5, 99.5));
303  add(h, new TH1F("nsimtrkPU","# of simulated tracks (PU)", 100, -0.5, 99.5));
304  h["nsimtrk"]->StatOverflows(kTRUE);
305  h["nsimtrkPU"]->StatOverflows(kTRUE);
306  h["nsimtrkSignal"]->StatOverflows(kTRUE);
307  h["xrectag"] = new TH1F("xrectag","reconstructed x, signal vtx",100,-0.05,0.05);
308  h["yrectag"] = new TH1F("yrectag","reconstructed y, signal vtx",100,-0.05,0.05);
309  h["zrectag"] = new TH1F("zrectag","reconstructed z, signal vtx",100,-20.,20.);
310  h["nrectrk0vtx"] = new TH1F("nrectrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
311  h["nseltrk0vtx"] = new TH1F("nseltrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
312  h["nrecsimtrk"] = new TH1F("nrecsimtrk","# rec tracks matched to sim tracks in vertex",100,-0.5, 99.5);
313  h["nrecnosimtrk"] = new TH1F("nrecsimtrk","# rec tracks not matched to sim tracks in vertex",100,-0.5, 99.5);
314  h["trackAssEffvsPt"] = new TProfile("trackAssEffvsPt","track association efficiency vs pt",20, 0., 100., 0, 1.);
315 
316  // cluster stuff
317  h["nseltrk"] = new TH1F("nseltrk","# of reconstructed tracks selected for PV", 100, -0.5, 99.5);
318  h["nclutrkall"] = new TH1F("nclutrkall","# of reconstructed tracks in clusters", 100, -0.5, 99.5);
319  h["nclutrkvtx"] = new TH1F("nclutrkvtx","# of reconstructed tracks in clusters of reconstructed vertices", 100, -0.5, 99.5);
320  h["nclu"] = new TH1F("nclu","# of clusters", 100, -0.5, 99.5);
321  h["nclu0vtx"] = new TH1F("nclu0vtx","# of clusters in events with no PV", 100, -0.5, 99.5);
322  h["zlost1"] = new TH1F("zlost1","z of lost vertices (bad z)", 100, -20., 20.);
323  h["zlost2"] = new TH1F("zlost2","z of lost vertices (no matching cluster)", 100, -20., 20.);
324  h["zlost3"] = new TH1F("zlost3","z of lost vertices (vertex too far from beam)", 100, -20., 20.);
325  h["zlost4"] = new TH1F("zlost4","z of lost vertices (invalid vertex)", 100, -20., 20.);
326  h["selstat"] = new TH1F("selstat","selstat", 5, -2.5, 2.5);
327 
328 
329  // properties of fake vertices (MC only)_
330  add(h, new TH1F("fakeVtxZNdofgt05","z of fake vertices with ndof>0.5", 100, -20., 20.));
331  add(h, new TH1F("fakeVtxZNdofgt2","z of fake vertices with ndof>2", 100, -20., 20.));
332  add(h, new TH1F("fakeVtxZ","z of fake vertices", 100, -20., 20.));
333  add(h, new TH1F("fakeVtxNdof","ndof of fake vertices", 500,0., 100.));
334  add(h,new TH1F("fakeVtxNtrk","number of tracks in fake vertex",20,-0.5, 19.5));
335  add(h,new TH1F("matchedVtxNdof","ndof of matched vertices", 500,0., 100.));
336 
337 
338  // histograms of track quality (Data and MC)
339  string types[] = {"all","sel","bachelor","sellost","wgt05","wlt05","M","tagged","untagged","ndof4","PUcand","PUfake"};
340  for(int t=0; t<12; t++){
341  add(h, new TH1F(("rapidity_"+types[t]).c_str(),"rapidity ",100,-3., 3.));
342  h["z0_"+types[t]] = new TH1F(("z0_"+types[t]).c_str(),"z0 ",80,-40., 40.);
343  h["phi_"+types[t]] = new TH1F(("phi_"+types[t]).c_str(),"phi ",80,-3.14159, 3.14159);
344  h["eta_"+types[t]] = new TH1F(("eta_"+types[t]).c_str(),"eta ",80,-4., 4.);
345  h["pt_"+types[t]] = new TH1F(("pt_"+types[t]).c_str(),"pt ",100,0., 20.);
346  h["found_"+types[t]] = new TH1F(("found_"+types[t]).c_str(),"found hits",20, 0., 20.);
347  h["lost_"+types[t]] = new TH1F(("lost_"+types[t]).c_str(),"lost hits",20, 0., 20.);
348  h["nchi2_"+types[t]] = new TH1F(("nchi2_"+types[t]).c_str(),"normalized track chi2",100, 0., 20.);
349  h["rstart_"+types[t]] = new TH1F(("rstart_"+types[t]).c_str(),"start radius",100, 0., 20.);
350  h["tfom_"+types[t]] = new TH1F(("tfom_"+types[t]).c_str(),"track figure of merit",100, 0., 100.);
351  h["expectedInner_"+types[t]] = new TH1F(("expectedInner_"+types[t]).c_str(),"expected inner hits ",10, 0., 10.);
352  h["expectedOuter_"+types[t]] = new TH1F(("expectedOuter_"+types[t]).c_str(),"expected outer hits ",10, 0., 10.);
353  h["logtresxy_"+types[t]] = new TH1F(("logtresxy_"+types[t]).c_str(),"log10(track r-phi resolution/um)",100, 0., 5.);
354  h["logtresz_"+types[t]] = new TH1F(("logtresz_"+types[t]).c_str(),"log10(track z resolution/um)",100, 0., 5.);
355  h["tpullxy_"+types[t]] = new TH1F(("tpullxy_"+types[t]).c_str(),"track r-phi pull",100, -10., 10.);
356  add(h, new TH2F( ("lvseta_"+types[t]).c_str(),"cluster length vs eta",60,-3., 3., 20, 0., 20));
357  add(h, new TH2F( ("lvstanlambda_"+types[t]).c_str(),"cluster length vs tan lambda",60,-6., 6., 20, 0., 20));
358  add(h, new TH1F( ("restrkz_"+types[t]).c_str(),"z-residuals (track vs vertex)", 200, -5., 5.));
359  add(h, new TH2F( ("restrkzvsphi_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
360  add(h, new TH2F( ("restrkzvseta_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.,3.,200, -5., 5.));
361  add(h, new TH2F( ("pulltrkzvsphi_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
362  add(h, new TH2F( ("pulltrkzvseta_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.,3.,100, -5., 5.));
363  add(h, new TH1F( ("pulltrkz_"+types[t]).c_str(),"normalized z-residuals (track vs vertex)", 100, -5., 5.));
364  add(h, new TH1F( ("sigmatrkz0_"+types[t]).c_str(),"z-resolution (excluding beam)", 100, 0., 5.));
365  add(h, new TH1F( ("sigmatrkz_"+types[t]).c_str(),"z-resolution (including beam)", 100,0., 5.));
366  add(h, new TH1F( ("nbarrelhits_"+types[t]).c_str(),"number of pixel barrel hits", 10, 0., 10.));
367  add(h, new TH1F( ("nbarrelLayers_"+types[t]).c_str(),"number of pixel barrel layers", 10, 0., 10.));
368  add(h, new TH1F( ("nPxLayers_"+types[t]).c_str(),"number of pixel layers (barrel+endcap)", 10, 0., 10.));
369  add(h, new TH1F( ("nSiLayers_"+types[t]).c_str(),"number of Tracker layers ", 20, 0., 20.));
370  add(h, new TH1F( ("trackAlgo_"+types[t]).c_str(),"track algorithm ", 30, 0., 30.));
371  add(h, new TH1F( ("trackQuality_"+types[t]).c_str(),"track quality ", 7, -1., 6.));
372  }
373  add(h, new TH1F("trackWt","track weight in vertex",100,0., 1.));
374 
375 
376  h["nrectrk"]->StatOverflows(kTRUE);
377  h["nrectrk"]->StatOverflows(kTRUE);
378  h["nrectrk0vtx"]->StatOverflows(kTRUE);
379  h["nseltrk0vtx"]->StatOverflows(kTRUE);
380  h["nrecsimtrk"]->StatOverflows(kTRUE);
381  h["nrecnosimtrk"]->StatOverflows(kTRUE);
382  h["nseltrk"]->StatOverflows(kTRUE);
383  h["nbtksinvtx"]->StatOverflows(kTRUE);
384  h["nbtksinvtxPU"]->StatOverflows(kTRUE);
385  h["nbtksinvtxTag"]->StatOverflows(kTRUE);
386  h["nbtksinvtx2"]->StatOverflows(kTRUE);
387  h["nbtksinvtxPU2"]->StatOverflows(kTRUE);
388  h["nbtksinvtxTag2"]->StatOverflows(kTRUE);
389 
390  // pile-up and track assignment related histograms (MC)
391  h["npu0"] =new TH1F("npu0","Number of simulated vertices",40,0.,40.);
392  h["npu1"] =new TH1F("npu1","Number of simulated vertices with >0 track",40,0.,40.);
393  h["npu2"] =new TH1F("npu2","Number of simulated vertices with >1 track",40,0.,40.);
394  h["nrecv"] =new TH1F("nrecv","# of reconstructed vertices", 40, 0, 40);
395  add(h,new TH2F("nrecvsnpu","#rec vertices vs number of sim vertices with >0 tracks", 40, 0., 40, 40, 0, 40));
396  add(h,new TH2F("nrecvsnpu2","#rec vertices vs number of sim vertices with >1 tracks", 40, 0., 40, 40, 0, 40));
397  add(h,new TH1F("sumpt","sumpt of simulated tracks",100,0.,100.));
398  add(h,new TH1F("sumptSignal","sumpt of simulated tracks in Signal events",100,0.,200.));
399  add(h,new TH1F("sumptPU","sumpt of simulated tracks in PU events",100,0.,200.));
400  add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
401  add(h,new TH1F("sumpt2","sumpt2 of simulated tracks",100,0.,100.));
402  add(h,new TH1F("sumpt2Signal","sumpt2 of simulated tracks in Signal events",100,0.,200.));
403  add(h,new TH1F("sumpt2PU","sumpt2 of simulated tracks in PU events",100,0.,200.));
404  add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
405  add(h,new TH1F("sumpt2recSignal","sumpt2 of reconstructed and matched tracks in Signal events",100,0.,200.));
406  add(h,new TH1F("sumpt2recPU","sumpt2 of reconstructed and matched tracks in PU events",100,0.,200.));
407  add(h,new TH1F("nRecTrkInSimVtx","number of reco tracks matched to sim-vertex", 101, 0., 100));
408  add(h,new TH1F("nRecTrkInSimVtxSignal","number of reco tracks matched to signal sim-vertex", 101, 0., 100));
409  add(h,new TH1F("nRecTrkInSimVtxPU","number of reco tracks matched to PU-vertex", 101, 0., 100));
410  add(h,new TH1F("nPrimRecTrkInSimVtx","number of reco primary tracks matched to sim-vertex", 101, 0., 100));
411  add(h,new TH1F("nPrimRecTrkInSimVtxSignal","number of reco primary tracks matched to signal sim-vertex", 101, 0., 100));
412  add(h,new TH1F("nPrimRecTrkInSimVtxPU","number of reco primary tracks matched to PU-vertex", 101, 0., 100));
413  // obsolete, scheduled for removal
414  add(h,new TH1F("recPurity","track purity of reconstructed vertices", 101, 0., 1.01));
415  add(h,new TH1F("recPuritySignal","track purity of reconstructed Signal vertices", 101, 0., 1.01));
416  add(h,new TH1F("recPurityPU","track purity of reconstructed PU vertices", 101, 0., 1.01));
417  // end of obsolete
418  add(h,new TH1F("recmatchPurity","track purity of all vertices", 101, 0., 1.01));
419  add(h,new TH1F("recmatchPurityTag","track purity of tagged vertices", 101, 0., 1.01));
420  add(h,new TH1F("recmatchPurityTagSignal","track purity of tagged Signal vertices", 101, 0., 1.01));
421  add(h,new TH1F("recmatchPurityTagPU","track purity of tagged PU vertices", 101, 0., 1.01));
422  add(h,new TH1F("recmatchPuritynoTag","track purity of untagged vertices", 101, 0., 1.01));
423  add(h,new TH1F("recmatchPuritynoTagSignal","track purity of untagged Signal vertices", 101, 0., 1.01));
424  add(h,new TH1F("recmatchPuritynoTagPU","track purity of untagged PU vertices", 101, 0., 1.01));
425  add(h,new TH1F("recmatchvtxs","number of sim vertices contributing to recvtx", 10, 0., 10.));
426  add(h,new TH1F("recmatchvtxsTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
427  add(h,new TH1F("recmatchvtxsnoTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
428  //
429  add(h,new TH1F("trkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
430  add(h,new TH1F("trkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
431  add(h,new TH1F("trkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
432  add(h,new TH1F("primtrkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
433  add(h,new TH1F("primtrkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
434  add(h,new TH1F("primtrkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
435  add(h,new TH1F("vtxMultiplicity", "number of rec vertices containg tracks from one true vertex", 10, 0., 10.) );
436  add(h,new TH1F("vtxMultiplicitySignal", "number of rec vertices containg tracks from the Signal Vertex", 10, 0., 10.) );
437  add(h,new TH1F("vtxMultiplicityPU", "number of rec vertices containg tracks from a PU Vertex", 10, 0., 10.) );
438 
439  add(h,new TProfile("vtxFindingEfficiencyVsNtrk","finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
440  add(h,new TProfile("vtxFindingEfficiencyVsNtrkSignal","Signal vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
441  add(h,new TProfile("vtxFindingEfficiencyVsNtrkPU","PU vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
442 
443  add(h,new TH1F("TagVtxTrkPurity","TagVtxTrkPurity",100,0.,1.01));
444  add(h,new TH1F("TagVtxTrkEfficiency","TagVtxTrkEfficiency",100,0.,1.01));
445 
446  add(h,new TH1F("matchVtxFraction","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
447  add(h,new TH1F("matchVtxFractionSignal","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
448  add(h,new TH1F("matchVtxFractionPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
449  add(h,new TH1F("matchVtxFractionCum","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
450  add(h,new TH1F("matchVtxFractionCumSignal","fraction of sim vertexs track found in a recvertex",101,0,1.01));
451  add(h,new TH1F("matchVtxFractionCumPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
452  add(h,new TH1F("matchVtxEfficiency","efficiency for finding matching rec vertex",2,-0.5,1.5));
453  add(h,new TH1F("matchVtxEfficiencySignal","efficiency for finding matching rec vertex",2,-0.5,1.5));
454  add(h,new TH1F("matchVtxEfficiencyPU","efficiency for finding matching rec vertex",2,-0.5,1.5));
455  add(h,new TH1F("matchVtxEfficiency2","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
456  add(h,new TH1F("matchVtxEfficiency2Signal","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
457  add(h,new TH1F("matchVtxEfficiency2PU","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
458  add(h,new TH1F("matchVtxEfficiency5","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
459  add(h,new TH1F("matchVtxEfficiency5Signal","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
460  add(h,new TH1F("matchVtxEfficiency5PU","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
461 
462 
463  add(h,new TH1F("matchVtxEfficiencyZ","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
464  add(h,new TH1F("matchVtxEfficiencyZSignal","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
465  add(h,new TH1F("matchVtxEfficiencyZPU","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
466 
467  add(h,new TH1F("matchVtxEfficiencyZ1","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
468  add(h,new TH1F("matchVtxEfficiencyZ1Signal","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
469  add(h,new TH1F("matchVtxEfficiencyZ1PU","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
470 
471  add(h,new TH1F("matchVtxEfficiencyZ2","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
472  add(h,new TH1F("matchVtxEfficiencyZ2Signal","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
473  add(h,new TH1F("matchVtxEfficiencyZ2PU","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
474 
475  add(h,new TH1F("matchVtxZ","z distance to matched recvtx",100, -0.1, 0.1));
476  add(h,new TH1F("matchVtxZPU","z distance to matched recvtx",100, -0.1, 0.1));
477  add(h,new TH1F("matchVtxZSignal","z distance to matched recvtx",100, -0.1, 0.1));
478 
479  add(h,new TH1F("matchVtxZCum","z distance to matched recvtx",1001,0, 1.01));
480  add(h,new TH1F("matchVtxZCumSignal","z distance to matched recvtx",1001,0, 1.01));
481  add(h,new TH1F("matchVtxZCumPU","z distance to matched recvtx",1001,0, 1.01));
482 
483  add(h, new TH1F("unmatchedVtx","number of unmatched rec vertices (fakes)",10,0.,10.));
484  add(h, new TH1F("unmatchedVtxFrac","fraction of unmatched rec vertices (fakes)",1000,0.,1.0));
485  add(h, new TH1F("unmatchedVtxZ","z of unmached rec vertices (fakes)",100,-20., 20.));
486  add(h, new TH1F("unmatchedVtxNdof","ndof of unmatched rec vertices (fakes)", 500,0., 100.));
487  add(h,new TH2F("correctlyassigned","pt and eta of correctly assigned tracks", 60, -3., 3., 100, 0, 10.));
488  add(h,new TH2F("misassigned","pt and eta of mis assigned tracks", 60, -3., 3., 100, 0, 10.));
489 
490  add(h,new TH1F("ptcat","pt of correctly assigned tracks", 100, 0, 10.));
491  add(h,new TH1F("etacat","eta of correctly assigned tracks", 60, -3., 3.));
492  add(h,new TH1F("phicat","phi of correctly assigned tracks", 100, -3.14159, 3.14159));
493  add(h,new TH1F("dzcat","dz of correctly assigned tracks", 100, 0., 1.));
494 
495  add(h,new TH1F("ptmis","pt of mis-assigned tracks", 100, 0, 10.));
496  add(h,new TH1F("etamis","eta of mis-assigned tracks", 60, -3., 3.));
497  add(h,new TH1F("phimis","phi of mis-assigned tracks",100, -3.14159, 3.14159));
498  add(h,new TH1F("dzmis","dz of mis-assigned tracks", 100, 0., 1.));
499 
500 
501  add(h,new TH1F("Tc","Tc computed with Truth matched Reco Tracks",100,0.,20.));
502  add(h,new TH1F("TcSignal","Tc of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
503  add(h,new TH1F("TcPU","Tc of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
504 
505  add(h,new TH1F("logTc","log Tc computed with Truth matched Reco Tracks",100,-2.,8.));
506  add(h,new TH1F("logTcSignal","log Tc of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
507  add(h,new TH1F("logTcPU","log Tc of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
508 
509  add(h,new TH1F("xTc","Tc of merged clusters",100,0.,20.));
510  add(h,new TH1F("xTcSignal","Tc of signal vertices merged with PU",100,0.,20.));
511  add(h,new TH1F("xTcPU","Tc of merged PU vertices",100,0.,20.));
512 
513  add(h,new TH1F("logxTc","log Tc merged vertices",100,-2.,8.));
514  add(h,new TH1F("logxTcSignal","log Tc of signal vertices merged with PU",100,-2.,8.));
515  add(h,new TH1F("logxTcPU","log Tc of merged PU vertices ",100,-2.,8.));
516 
517  add(h,new TH1F("logChisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,-2.,8.));
518  add(h,new TH1F("logChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
519  add(h,new TH1F("logChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
520 
521  add(h,new TH1F("logxChisq","Chisq/ntrk of merged clusters",100,-2.,8.));
522  add(h,new TH1F("logxChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,-2.,8.));
523  add(h,new TH1F("logxChisqPU","Chisq/ntrk of merged PU vertices",100,-2.,8.));
524 
525  add(h,new TH1F("Chisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,0.,20.));
526  add(h,new TH1F("ChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
527  add(h,new TH1F("ChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
528 
529  add(h,new TH1F("xChisq","Chisq/ntrk of merged clusters",100,0.,20.));
530  add(h,new TH1F("xChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,0.,20.));
531  add(h,new TH1F("xChisqPU","Chisq/ntrk of merged PU vertices",100,0.,20.));
532 
533  add(h,new TH1F("dzmax","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
534  add(h,new TH1F("dzmaxSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
535  add(h,new TH1F("dzmaxPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
536 
537  add(h,new TH1F("xdzmax","dzmax of merged clusters",100,0.,2.));
538  add(h,new TH1F("xdzmaxSignal","dzmax of signal vertices merged with PU",100,0.,2.));
539  add(h,new TH1F("xdzmaxPU","dzmax of merged PU vertices",100,0.,2.));
540 
541  add(h,new TH1F("dztrim","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
542  add(h,new TH1F("dztrimSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
543  add(h,new TH1F("dztrimPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
544 
545  add(h,new TH1F("xdztrim","dzmax of merged clusters",100,0.,2.));
546  add(h,new TH1F("xdztrimSignal","dzmax of signal vertices merged with PU",100,0.,2.));
547  add(h,new TH1F("xdztrimPU","dzmax of merged PU vertices",100,0.,2.));
548 
549  add(h,new TH1F("m4m2","m4m2 computed with Truth matched Reco Tracks",100,0.,100.));
550  add(h,new TH1F("m4m2Signal","m4m2 of signal vertices computed with Truth matched Reco Tracks",100,0.,100.));
551  add(h,new TH1F("m4m2PU","m4m2 of PU vertices computed with Truth matched Reco Tracks",100,0.,100.));
552 
553  add(h,new TH1F("xm4m2","m4m2 of merged clusters",100,0.,100.));
554  add(h,new TH1F("xm4m2Signal","m4m2 of signal vertices merged with PU",100,0.,100.));
555  add(h,new TH1F("xm4m2PU","m4m2 of merged PU vertices",100,0.,100.));
556 
557  return h;
558 }
559 
560 
561 //
562 // member functions
563 //
564 void PrimaryVertexAnalyzer4PU::beginJob(){
565  std::cout << " PrimaryVertexAnalyzer4PU::beginJob conversion from sim units to rec units is " << simUnit_ << std::endl;
566 
567 
568  rootFile_->cd();
569 
570  TDirectory *noBS = rootFile_->mkdir("noBS");
571  noBS->cd();
572  hnoBS=bookVertexHistograms();
573  for(std::map<std::string,TH1*>::const_iterator hist=hnoBS.begin(); hist!=hnoBS.end(); hist++){
574  hist->second->SetDirectory(noBS);
575  }
576 
577  TDirectory *withBS = rootFile_->mkdir("BS");
578  withBS->cd();
579  hBS=bookVertexHistograms();
580  for(std::map<std::string,TH1*>::const_iterator hist=hBS.begin(); hist!=hBS.end(); hist++){
581  hist->second->SetDirectory(withBS);
582  }
583 
584  TDirectory *DA = rootFile_->mkdir("DA");
585  DA->cd();
586  hDA=bookVertexHistograms();
587  for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
588  hist->second->SetDirectory(DA);
589  }
590 
591 // TDirectory *PIX = rootFile_->mkdir("PIX");
592 // PIX->cd();
593 // hPIX=bookVertexHistograms();
594 // for(std::map<std::string,TH1*>::const_iterator hist=hPIX.begin(); hist!=hPIX.end(); hist++){
595 // hist->second->SetDirectory(PIX);
596 // }
597 
598 // TDirectory *MVF = rootFile_->mkdir("MVF");
599 // MVF->cd();
600 // hMVF=bookVertexHistograms();
601 // for(std::map<std::string,TH1*>::const_iterator hist=hMVF.begin(); hist!=hMVF.end(); hist++){
602 // hist->second->SetDirectory(MVF);
603 // }
604 
605  rootFile_->cd();
606  hsimPV["rapidity"] = new TH1F("rapidity","rapidity ",100,-10., 10.);
607  hsimPV["chRapidity"] = new TH1F("chRapidity","charged rapidity ",100,-10., 10.);
608  hsimPV["recRapidity"] = new TH1F("recRapidity","reconstructed rapidity ",100,-10., 10.);
609  hsimPV["pt"] = new TH1F("pt","pt ",100,0., 20.);
610 
611  hsimPV["xsim"] = new TH1F("xsim","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
612  hsimPV["ysim"] = new TH1F("ysim","simulated y",100,-0.01,0.01);
613  hsimPV["zsim"] = new TH1F("zsim","simulated z",100,-20.,20.);
614 
615  hsimPV["xsim1"] = new TH1F("xsim1","simulated x",100,-4.,4.);
616  hsimPV["ysim1"] = new TH1F("ysim1","simulated y",100,-4.,4.);
617  hsimPV["zsim1"] = new TH1F("zsim1","simulated z",100,-40.,40.);
618 
619  add(hsimPV, new TH1F("xsim2PU","simulated x (Pile-up)",100,-1.,1.));
620  add(hsimPV, new TH1F("ysim2PU","simulated y (Pile-up)",100,-1.,1.));
621  add(hsimPV, new TH1F("zsim2PU","simulated z (Pile-up)",100,-20.,20.));
622  add(hsimPV, new TH1F("xsim2Signal","simulated x (Signal)",100,-1.,1.));
623  add(hsimPV, new TH1F("ysim2Signal","simulated y (Signal)",100,-1.,1.));
624  add(hsimPV, new TH1F("zsim2Signal","simulated z (Signal)",100,-20.,20.));
625 
626  hsimPV["xsim2"] = new TH1F("xsim2","simulated x",100,-1,1); // 0.01cm = 100 um
627  hsimPV["ysim2"] = new TH1F("ysim2","simulated y",100,-1,1);
628  hsimPV["zsim2"] = new TH1F("zsim2","simulated z",100,-20.,20.);
629  hsimPV["xsim3"] = new TH1F("xsim3","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
630  hsimPV["ysim3"] = new TH1F("ysim3","simulated y",100,-0.1,0.1);
631  hsimPV["zsim3"] = new TH1F("zsim3","simulated z",100,-20.,20.);
632  hsimPV["xsimb"] = new TH1F("xsimb","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
633  hsimPV["ysimb"] = new TH1F("ysimb","simulated y",100,-0.01,0.01);
634  hsimPV["zsimb"] = new TH1F("zsimb","simulated z",100,-20.,20.);
635  hsimPV["xsimb1"] = new TH1F("xsimb1","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
636  hsimPV["ysimb1"] = new TH1F("ysimb1","simulated y",100,-0.1,0.1);
637  hsimPV["zsimb1"] = new TH1F("zsimb1","simulated z",100,-20.,20.);
638  add(hsimPV, new TH1F("xbeam","beamspot x",100,-1.,1.));
639  add(hsimPV, new TH1F("ybeam","beamspot y",100,-1.,1.));
640  add(hsimPV, new TH1F("zbeam","beamspot z",100,-1.,1.));
641  add(hsimPV, new TH1F("wxbeam","beamspot sigma x",100,-1.,1.));
642  add(hsimPV, new TH1F("wybeam","beamspot sigma y",100,-1.,1.));
643  add(hsimPV, new TH1F("wzbeam","beamspot sigma z",100,-1.,1.));
644  hsimPV["xsim2"]->StatOverflows(kTRUE);
645  hsimPV["ysim2"]->StatOverflows(kTRUE);
646  hsimPV["zsim2"]->StatOverflows(kTRUE);
647  hsimPV["xsimb"]->StatOverflows(kTRUE);
648  hsimPV["ysimb"]->StatOverflows(kTRUE);
649  hsimPV["zsimb"]->StatOverflows(kTRUE);
650  hsimPV["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
651  // hsimPV["nsimtrk"] = new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5); // not filled right now, also exists in hBS..
652  // hsimPV["nsimtrk"]->StatOverflows(kTRUE);
653  hsimPV["nbsimtksinvtx"]= new TH1F("nbsimtksinvtx","simulated tracks in vertex",100,-0.5,99.5);
654  hsimPV["nbsimtksinvtx"]->StatOverflows(kTRUE);
655 
656 }
657 
658 
659 void PrimaryVertexAnalyzer4PU::endJob() {
660  std::cout << "this is void PrimaryVertexAnalyzer4PU::endJob() " << std::endl;
661  //cumulate some histos
662  double sumDA=0,sumBS=0,sumnoBS=0;
663  // double sumPIX=0,sumMVF=0;
664  for(int i=101; i>0; i--){
665  sumDA+=hDA["matchVtxFractionSignal"]->GetBinContent(i)/hDA["matchVtxFractionSignal"]->Integral();
666  hDA["matchVtxFractionCumSignal"]->SetBinContent(i,sumDA);
667  sumBS+=hBS["matchVtxFractionSignal"]->GetBinContent(i)/hBS["matchVtxFractionSignal"]->Integral();
668  hBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumBS);
669  sumnoBS+=hnoBS["matchVtxFractionSignal"]->GetBinContent(i)/hnoBS["matchVtxFractionSignal"]->Integral();
670  hnoBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumnoBS);
671 // sumPIX+=hPIX["matchVtxFractionSignal"]->GetBinContent(i)/hPIX["matchVtxFractionSignal"]->Integral();
672 // hPIX["matchVtxFractionCumSignal"]->SetBinContent(i,sumPIX);
673 // sumMVF+=hMVF["matchVtxFractionSignal"]->GetBinContent(i)/hMVF["matchVtxFractionSignal"]->Integral();
674 // hMVF["matchVtxFractionCumSignal"]->SetBinContent(i,sumMVF);
675  }
676  sumDA=0,sumBS=0,sumnoBS=0;
677  //sumPIX=0,sumMVF=0;
678  for(int i=1; i<1001; i++){
679  sumDA+=hDA["abszdistancetag"]->GetBinContent(i);
680  hDA["abszdistancetagcum"]->SetBinContent(i,sumDA/float(hDA["abszdistancetag"]->GetEntries()));
681  sumBS+=hBS["abszdistancetag"]->GetBinContent(i);
682  hBS["abszdistancetagcum"]->SetBinContent(i,sumBS/float(hBS["abszdistancetag"]->GetEntries()));
683  sumnoBS+=hnoBS["abszdistancetag"]->GetBinContent(i);
684  hnoBS["abszdistancetagcum"]->SetBinContent(i,sumnoBS/float(hnoBS["abszdistancetag"]->GetEntries()));
685 // sumPIX+=hPIX["abszdistancetag"]->GetBinContent(i);
686 // hPIX["abszdistancetagcum"]->SetBinContent(i,sumPIX/float(hPIX["abszdistancetag"]->GetEntries()));
687 // sumMVF+=hMVF["abszdistancetag"]->GetBinContent(i);
688 // hMVF["abszdistancetagcum"]->SetBinContent(i,sumMVF/float(hMVF["abszdistancetag"]->GetEntries()));
689  }
690 
691  Cumulate(hBS["matchVtxZCum"]); Cumulate(hBS["matchVtxZCumSignal"]); Cumulate(hBS["matchVtxZCumPU"]);
692  Cumulate(hnoBS["matchVtxZCum"]); Cumulate(hnoBS["matchVtxZCumSignal"]); Cumulate(hnoBS["matchVtxZCumPU"]);
693  Cumulate(hDA["matchVtxZCum"]); Cumulate(hDA["matchVtxZCumSignal"]); Cumulate(hDA["matchVtxZCumPU"]);
694  /*
695  h->ComputeIntegral();
696  Double_t *integral = h->GetIntegral();
697  h->SetContent(integral);
698  */
699 
700  // make a reference for ndofnr2
701  //hDA["vtxndof"]->ComputeIntegral();
702  //Double_t *integral = hDA["vtxndf"]->GetIntegral();
703  // h->SetContent(integral);
704  double p;
705  for(int i=1; i<501; i++){
706  if(hDA["vtxndf"]->GetEntries()>0){
707  p= hDA["vtxndf"]->Integral(i,501)/hDA["vtxndf"]->GetEntries(); hDA["vtxndfc"]->SetBinContent(i,p*hDA["vtxndf"]->GetBinContent(i));
708  }
709  if(hBS["vtxndf"]->GetEntries()>0){
710  p= hBS["vtxndf"]->Integral(i,501)/hBS["vtxndf"]->GetEntries(); hBS["vtxndfc"]->SetBinContent(i,p*hBS["vtxndf"]->GetBinContent(i));
711  }
712  if(hnoBS["vtxndf"]->GetEntries()>0){
713  p=hnoBS["vtxndf"]->Integral(i,501)/hnoBS["vtxndf"]->GetEntries(); hnoBS["vtxndfc"]->SetBinContent(i,p*hnoBS["vtxndf"]->GetBinContent(i));
714  }
715 // if(hPIX["vtxndf"]->GetEntries()>0){
716 // p=hPIX["vtxndf"]->Integral(i,501)/hPIX["vtxndf"]->GetEntries(); hPIX["vtxndfc"]->SetBinContent(i,p*hPIX["vtxndf"]->GetBinContent(i));
717 // }
718  }
719 
720  rootFile_->cd();
721  for(std::map<std::string,TH1*>::const_iterator hist=hsimPV.begin(); hist!=hsimPV.end(); hist++){
722  std::cout << "writing " << hist->first << std::endl;
723  hist->second->Write();
724  }
725  rootFile_->Write();
726  std::cout << "PrimaryVertexAnalyzer4PU::endJob: done" << std::endl;
727 }
728 
729 
730 
731 
732 // helper functions
733 std::vector<PrimaryVertexAnalyzer4PU::SimPart> PrimaryVertexAnalyzer4PU::getSimTrkParameters(
734  edm::Handle<edm::SimTrackContainer> & simTrks,
735  edm::Handle<edm::SimVertexContainer> & simVtcs,
736  double simUnit)
737 {
738  std::vector<SimPart > tsim;
739  if(simVtcs->begin()==simVtcs->end()){
740  if(verbose_){
741  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters no simvtcs" << endl;
742  }
743  return tsim;
744  }
745  if(verbose_){
746  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters simVtcs n=" << simVtcs->size() << endl;
747  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters 1st position" << setw(8) << setprecision(4) << simVtcs->begin()->position() << endl;
748  }
749  double t0=simVtcs->begin()->position().e();
750 
751  for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
752  t!=simTrks->end(); ++t){
753  if (t->noVertex()){
754  std::cout << "simtrk has no vertex" << std::endl;
755  }else{
756  // get the vertex position
757  //HepLorentzVector v=(*simVtcs)[t->vertIndex()].position();
758  math::XYZTLorentzVectorD v((*simVtcs)[t->vertIndex()].position().x(),
759  (*simVtcs)[t->vertIndex()].position().y(),
760  (*simVtcs)[t->vertIndex()].position().z(),
761  (*simVtcs)[t->vertIndex()].position().e());
762  int pdgCode=t->type();
763 
764  if( pdgCode==-99 ){
765  // such entries cause crashes, no idea what they are
766  std::cout << "funny particle skipped , code=" << pdgCode << std::endl;
767  }else{
768  double Q=0; //double Q=HepPDT::theTable().getParticleData(pdgCode)->charge();
769  if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==-321)||(pdgCode==-3222)){Q=-1;}
770  else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)||(pdgCode==3222)){Q=1;}
771  else {
772  //std::cout << pdgCode << " " <<std::endl;
773  }
774  math::XYZTLorentzVectorD p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
775  if ( (Q != 0) && (p.pt()>0.1) && (fabs(t->momentum().eta())<3.0)
776  && fabs(v.z()*simUnit<20) && (sqrt(v.x()*v.x()+v.y()*v.y())<10.)){
777  double x0=v.x()*simUnit;
778  double y0=v.y()*simUnit;
779  double z0=v.z()*simUnit;
780  double kappa=-Q*0.002998*fBfield_/p.pt();
781  double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
782  double q=sqrt(1.-2.*kappa*D0);
783  double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
784  double s1;
785  if (fabs(kappa*s0)>0.001){
786  s1=asin(kappa*s0)/kappa;
787  }else{
788  double ks02=(kappa*s0)*(kappa*s0);
789  s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
790  }
791  SimPart sp;//ParameterVector par;
792  sp.par[reco::TrackBase::i_qoverp] = Q/p.P();
793  sp.par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
794  sp.par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
795  sp.par[reco::TrackBase::i_dxy] = -2.*D0/(1.+q);
796  sp.par[reco::TrackBase::i_dsz] = z0*sin(p.theta())-s1*cos(p.theta());
797 
798  sp.pdg=pdgCode;
799  if (v.t()-t0<1e-15){
800  sp.type=0; // primary
801  }else{
802  sp.type=1; //secondary
803  }
804 
805  // now get zpca (get perigee wrt beam)
806  double x1=x0-0.033; double y1=y0-0.; // FIXME how do we get the simulated beam position?
807  D0=x1*sin(p.phi())-y1*cos(p.phi())-0.5*kappa*(x1*x1+y1*y1);
808  q=sqrt(1.-2.*kappa*D0);
809  s0=(x1*cos(p.phi())+y1*sin(p.phi()))/q;
810  if (fabs(kappa*s0)>0.001){
811  s1=asin(kappa*s0)/kappa;
812  }else{
813  double ks02=(kappa*s0)*(kappa*s0);
814  s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
815  }
816  sp.ddcap=-2.*D0/(1.+q);
817  sp.zdcap=z0-s1/tan(p.theta());
818  sp.zvtx=z0;
819  sp.xvtx=x0;
820  sp.yvtx=y0;
821 
822  tsim.push_back(sp);
823  }
824  }
825  }// has vertex
826  }//for loop
827  return tsim;
828 }
829 
830 
831 int* PrimaryVertexAnalyzer4PU::supf(std::vector<SimPart>& simtrks, const reco::TrackCollection & trks){
832  int nsim=simtrks.size();
833  int nrec=trks.size();
834  int *rectosim=new int[nrec]; // pointer to associated simtrk
835  double** pij=new double*[nrec];
836  double mu=100.; // initial chi^2 cut-off (5 dofs !)
837  int nmatch=0;
838  int i=0;
839  for(reco::TrackCollection::const_iterator t=trks.begin(); t!=trks.end(); ++t){
840  pij[i]=new double[nsim];
841  rectosim[i]=-1;
842  ParameterVector par = t->parameters();
843  //reco::TrackBase::CovarianceMatrix V = t->covariance();
844  reco::TrackBase::CovarianceMatrix S = t->covariance();
845  S.Invert();
846  for(int j=0; j<nsim; j++){
847  simtrks[j].rec=-1;
848  SimPart s=simtrks[j];
849  double c=0;
850  for(int k=0; k<5; k++){
851  for(int l=0; l<5; l++){
852  c+=(par(k)-s.par[k])*(par(l)-s.par[l])*S(k,l);
853  }
854  }
855  pij[i][j]=exp(-0.5*c);
856 
857 // double c0=pow((par[0]-s.par[0])/t->qoverpError(),2)*0.1
858 // +pow((par[1]-s.par[1])/t->lambdaError(),2)
859 // +pow((par[2]-s.par[2])/t->phiError(),2)
860 // +pow((par[3]-s.par[3])/t->dxyError(),2)*0.1;
861 // +pow((par[4]-s.par[4])/t->dszError(),2)*0.1;
862 // pij[i][j]=exp(-0.5*c0);
863 
864 // if( c0 <100 ){
865 // cout << setw(3) << i << " rec " << setw(6) << par << endl;
866 // cout << setw(3) << j << " sim " << setw(6) << s.par << " ---> C=" << c << endl;
867 // cout << " " << setw(6)
868 // << (par[0]-s.par[0])<< ","
869 // << (par[1]-s.par[1])<< ","
870 // << (par[2]-s.par[2])<< ","
871 // << (par[3]-s.par[3])<< ","
872 // << (par[4]-s.par[4])
873 // << " match=" << match(simtrks[j].par, trks.at(i).parameters())
874 // << endl;
875 // cout << " " << setw(6)
876 // << (par[0]-s.par[0])/t->qoverpError() << ","
877 // << (par[1]-s.par[1])/t->lambdaError() << ","
878 // << (par[2]-s.par[2])/t->phiError() << ","
879 // << (par[3]-s.par[3])/t->dxyError() << ","
880 // << (par[4]-s.par[4])/t->dszError() << " c0=" << c0
881 // << endl <<endl;
882 // }
883 
884  }
885  i++;
886  }
887 
888  for(int k=0; k<nrec; k++){
889  int imatch=-1; int jmatch=-1;
890  double pmatch=0;
891  for(int j=0; j<nsim; j++){
892  if ((simtrks[j].rec)<0){
893  double psum=exp(-0.5*mu); //cutoff
894  for(int i=0; i<nrec; i++){
895  if (rectosim[i]<0){ psum+=pij[i][j];}
896  }
897  for(int i=0; i<nrec; i++){
898  if ((rectosim[i]<0)&&(pij[i][j]/psum>pmatch)){
899  pmatch=pij[i][j]/psum;
900  imatch=i; jmatch=j;
901  }
902  }
903  }
904  }// sim loop
905  if((jmatch>=0)||(imatch>=0)){
906  //std::cout << pmatch << " " << pij[imatch][jmatch] << " match=" <<
907  // match(simtrks[jmatch].par, trks.at(imatch).parameters()) <<std::endl;
908  if (pmatch>0.01){
909  rectosim[imatch]=jmatch;
910  simtrks[jmatch].rec=imatch;
911  nmatch++;
912  }else if (match(simtrks[jmatch].par, trks.at(imatch).parameters())){
913  // accept it anyway if it matches crudely and relax the cut-off
914  rectosim[imatch]=jmatch;
915  simtrks[jmatch].rec=imatch;
916  nmatch++;
917  mu=mu*2;
918  }
919  }
920  }
921 
922 // std::cout << ">>>>>>>>>>>>>>>--------------supf----------------------" << std::endl;
923 // std::cout <<"nsim=" << nsim << " nrec=" << nrec << " nmatch=" << nmatch << std::endl;
924 // std::cout << "rec to sim " << std::endl;
925 // for(int i=0; i<nrec; i++){
926 // std::cout << i << " ---> " << rectosim[i] << std::endl;
927 // }
928 // std::cout << "sim to rec " << std::endl;
929 // for(int j=0; j<nsim; j++){
930 // std::cout << j << " ---> " << simtrks[j].rec << std::endl;
931 // }
932 
933  std::cout << "unmatched sim " << std::endl;
934  for(int j=0; j<nsim; j++){
935  if(simtrks[j].rec<0){
936  double pt= 1./simtrks[j].par[0]/tan(simtrks[j].par[1]);
937  if((fabs(pt))>1.){
938  std::cout << setw(3) << j << setw(8) << simtrks[j].pdg
939  << setw(8) << setprecision(4) << " (" << simtrks[j].xvtx << "," << simtrks[j].yvtx << "," << simtrks[j].zvtx << ")"
940  << " pt= " << pt
941  << " phi=" << simtrks[j].par[2]
942  << " eta= " << -log(tan(0.5*(M_PI/2-simtrks[j].par[1])))
943  << std::endl;
944  }
945  }
946  }
947 // std::cout << "<<<<<<<<<<<<<<<--------------supf----------------------" << std::endl;
948 
949  //delete rectosim; // or return it?
950  for(int i=0; i<nrec; i++){delete pij[i];}
951  delete pij;
952  return rectosim; // caller must delete it
953 }
954 
955 
956 
957 
958 
959 
960 
961 
962 bool PrimaryVertexAnalyzer4PU::match(const ParameterVector &a,
963  const ParameterVector &b){
964  double dqoverp =a(0)-b(0);
965  double dlambda =a(1)-b(1);
966  double dphi =a(2)-b(2);
967  double dsz =a(4)-b(4);
968  if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
969  // return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<0.1) );
970  return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<1.0) );
971 }
972 
973 
974 bool PrimaryVertexAnalyzer4PU::matchVertex(const simPrimaryVertex &vsim,
975  const reco::Vertex &vrec){
976  return (fabs(vsim.z*simUnit_-vrec.z())<zmatch_);
977 }
978 
979 bool PrimaryVertexAnalyzer4PU::isResonance(const HepMC::GenParticle * p){
980  double ctau=(pdt_->particle( abs(p->pdg_id()) ))->lifetime();
981  //std::cout << "isResonance " << p->pdg_id() << " " << ctau << std::endl;
982  return ctau >0 && ctau <1e-6;
983 }
984 
985 bool PrimaryVertexAnalyzer4PU::isFinalstateParticle(const HepMC::GenParticle * p){
986  return ( !p->end_vertex() && p->status()==1 );
987 }
988 
989 
990 bool PrimaryVertexAnalyzer4PU::isCharged(const HepMC::GenParticle * p){
991  const ParticleData * part = pdt_->particle( p->pdg_id() );
992  if (part){
993  return part->charge()!=0;
994  }else{
995  // the new/improved particle table doesn't know anti-particles
996  return pdt_->particle( -p->pdg_id() )!=0;
997  }
998 }
999 
1000 
1001 
1002 
1003 void PrimaryVertexAnalyzer4PU::fillTrackHistos(std::map<std::string, TH1*> & h, const std::string & ttype, const reco::Track & t, const reco::Vertex * v){
1004  Fill(h,"rapidity_"+ttype,t.eta());
1005  Fill(h,"z0_"+ttype,t.vz());
1006  Fill(h,"phi_"+ttype,t.phi());
1007  Fill(h,"eta_"+ttype,t.eta());
1008  Fill(h,"pt_"+ttype,t.pt());
1009  Fill(h,"found_"+ttype,t.found());
1010  Fill(h,"lost_"+ttype,t.lost());
1011  Fill(h,"nchi2_"+ttype,t.normalizedChi2());
1012  Fill(h,"rstart_"+ttype,(t.innerPosition()).Rho());
1013 
1014  double d0Error=t.d0Error();
1015  double d0=t.dxy(vertexBeamSpot_.position());
1016  if (d0Error>0){
1017  Fill(h,"logtresxy_"+ttype,log(d0Error/0.0001)/log(10.));
1018  Fill(h,"tpullxy_"+ttype,d0/d0Error);
1019  }
1020  //double z0=t.vz();
1021  double dzError=t.dzError();
1022  if(dzError>0){
1023  Fill(h,"logtresz_"+ttype,log(dzError/0.0001)/log(10.));
1024  }
1025 
1026  //
1027  Fill(h,"sigmatrkz_"+ttype,sqrt(pow(t.dzError(),2)+wxy2_/pow(tan(t.theta()),2)));
1028  Fill(h,"sigmatrkz0_"+ttype,t.dzError());
1029 
1030  // track vs vertex
1031  if((! (v==NULL)) && (v->ndof()>10.)) {
1032  // emulate clusterizer input
1033  //const TransientTrack & tt = theB_->build(&t); wrong !!!!
1034  TransientTrack tt = theB_->build(&t); tt.setBeamSpot(vertexBeamSpot_); // need the setBeamSpot !
1035  double z=(tt.stateAtBeamLine().trackStateAtPCA()).position().z();
1036  double tantheta=tan((tt.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1037  double dz2= pow(tt.track().dzError(),2)+wxy2_/pow(tantheta,2);
1038 
1039  Fill(h,"restrkz_"+ttype,z-v->position().z());
1040  Fill(h,"restrkzvsphi_"+ttype,t.phi(), z-v->position().z());
1041  Fill(h,"restrkzvseta_"+ttype,t.eta(), z-v->position().z());
1042  Fill(h,"pulltrkzvsphi_"+ttype,t.phi(), (z-v->position().z())/sqrt(dz2));
1043  Fill(h,"pulltrkzvseta_"+ttype,t.eta(), (z-v->position().z())/sqrt(dz2));
1044 
1045  Fill(h,"pulltrkz_"+ttype,(z-v->position().z())/sqrt(dz2));
1046 
1047  double x1=t.vx()-vertexBeamSpot_.x0(); double y1=t.vy()-vertexBeamSpot_.y0();
1048  double kappa=-0.002998*fBfield_*t.qoverp()/cos(t.theta());
1049  double D0=x1*sin(t.phi())-y1*cos(t.phi())-0.5*kappa*(x1*x1+y1*y1);
1050  double q=sqrt(1.-2.*kappa*D0);
1051  double s0=(x1*cos(t.phi())+y1*sin(t.phi()))/q;
1052  // double s1;
1053  if (fabs(kappa*s0)>0.001){
1054  //s1=asin(kappa*s0)/kappa;
1055  }else{
1056  //double ks02=(kappa*s0)*(kappa*s0);
1057  //s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
1058  }
1059  // sp.ddcap=-2.*D0/(1.+q);
1060  //double zdcap=t.vz()-s1/tan(t.theta());
1061 
1062  }
1063  //
1064 
1065  // collect some info on hits and clusters
1066  Fill(h,"nbarrelLayers_"+ttype,t.hitPattern().pixelBarrelLayersWithMeasurement());
1067  Fill(h,"nPxLayers_"+ttype,t.hitPattern().pixelLayersWithMeasurement());
1068  Fill(h,"nSiLayers_"+ttype,t.hitPattern().trackerLayersWithMeasurement());
1069  Fill(h,"expectedInner_"+ttype,t.trackerExpectedHitsInner().numberOfHits());
1070  Fill(h,"expectedOuter_"+ttype,t.trackerExpectedHitsOuter().numberOfHits());
1071  Fill(h,"trackAlgo_"+ttype,t.algo());
1072  Fill(h,"trackQuality_"+ttype,t.qualityMask());
1073 
1074  //
1075  int longesthit=0, nbarrel=0;
1076  for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
1077  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1078  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1079  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1080  if (barrel){
1081  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1082  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1083  if (clust.isNonnull()) {
1084  nbarrel++;
1085  if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
1086  if (clust->sizeY()>20.){
1087  Fill(h,"lvseta_"+ttype,t.eta(), 19.9);
1088  Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), 19.9);
1089  }else{
1090  Fill(h,"lvseta_"+ttype,t.eta(), float(clust->sizeY()));
1091  Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), float(clust->sizeY()));
1092  }
1093  }
1094  }
1095  }
1096  }
1097  Fill(h,"nbarrelhits_"+ttype,float(nbarrel));
1098  //-------------------------------------------------------------------
1099 
1100 }
1101 
1102 void PrimaryVertexAnalyzer4PU::dumpHitInfo(const reco::Track & t){
1103  // collect some info on hits and clusters
1104  int longesthit=0, nbarrel=0;
1105  cout << Form("%5.2f %5.2f : ",t.pt(),t.eta());
1106  for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
1107  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1108  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1109  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1110  if (barrel){
1111  nbarrel++;
1112  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1113  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1114  if (clust.isNonnull()) {
1115  cout << Form("%4d",clust->sizeY());
1116  if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
1117  }
1118  }
1119  }
1120  }
1121  //cout << endl;
1122 }
1123 
1124 void PrimaryVertexAnalyzer4PU::printRecVtxs(const Handle<reco::VertexCollection> recVtxs, std::string title){
1125  int ivtx=0;
1126  std::cout << std::endl << title << std::endl;
1127  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1128  v!=recVtxs->end(); ++v){
1129  string vtype=" recvtx ";
1130  if( v->isFake()){
1131  vtype=" fake ";
1132  }else if (v->ndof()==-5){
1133  vtype=" cluster "; // pos=selector[iclu],cputime[iclu],clusterz[iclu]
1134  }else if(v->ndof()==-3){
1135  vtype=" event ";
1136  }
1137  std::cout << "vtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
1138  << vtype
1139  << " #trk " << std::fixed << std::setprecision(4) << std::setw(3) << v->tracksSize()
1140  << " chi2 " << std::fixed << std::setw(4) << std::setprecision(1) << v->chi2()
1141  << " ndof " << std::fixed << std::setw(5) << std::setprecision(2) << v->ndof()
1142  << " x " << std::setw(8) <<std::fixed << std::setprecision(4) << v->x()
1143  << " dx " << std::setw(8) << v->xError()// << std::endl
1144  << " y " << std::setw(8) << v->y()
1145  << " dy " << std::setw(8) << v->yError()//<< std::endl
1146  << " z " << std::setw(8) << v->z()
1147  << " dz " << std::setw(8) << v->zError()
1148  << std::endl;
1149  }
1150 }
1151 
1152 
1153 void PrimaryVertexAnalyzer4PU::printSimVtxs(const Handle<SimVertexContainer> simVtxs){
1154  int i=0;
1155  for(SimVertexContainer::const_iterator vsim=simVtxs->begin();
1156  vsim!=simVtxs->end(); ++vsim){
1157  if ( vsim->position().x()*vsim->position().x()+vsim->position().y()*vsim->position().y() < 1.){
1158  std::cout << i++ << ")" << std::scientific
1159  << " evtid=" << vsim->eventId().event() << "," << vsim->eventId().bunchCrossing()
1160  << " sim x=" << vsim->position().x()*simUnit_
1161  << " sim y=" << vsim->position().y()*simUnit_
1162  << " sim z=" << vsim->position().z()*simUnit_
1163  << " sim t=" << vsim->position().t()
1164  << " parent=" << vsim->parentIndex()
1165  << std::endl;
1166  }
1167  }
1168 }
1169 
1170 
1171 
1172 
1173 
1174 
1175 
1176 void PrimaryVertexAnalyzer4PU::printSimTrks(const Handle<SimTrackContainer> simTrks){
1177  std::cout << " simTrks type, (momentum), vertIndex, genpartIndex" << std::endl;
1178  int i=1;
1179  for(SimTrackContainer::const_iterator t=simTrks->begin();
1180  t!=simTrks->end(); ++t){
1181  //HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
1182  std::cout << i++ << ")"
1183  << t->eventId().event() << "," << t->eventId().bunchCrossing()
1184  << (*t)
1185  << " index="
1186  << (*t).genpartIndex();
1187  //if (gp) {
1188  // HepMC::GenVertex *gv=gp->production_vertex();
1189  // std::cout << " genvertex =" << (*gv);
1190  //}
1191  std::cout << std::endl;
1192  }
1193 }
1194 
1195 
1196 
1197 void PrimaryVertexAnalyzer4PU::printRecTrks(const Handle<reco::TrackCollection> &recTrks ){
1198 
1199  cout << "printRecTrks" << endl;
1200  int i =1;
1201  for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
1202  // reco::TrackBase::ParameterVector par = t->parameters();
1203 
1204  cout << endl;
1205  cout << "Track "<<i << " " ; i++;
1206  //enum TrackQuality { undefQuality=-1, loose=0, tight=1, highPurity=2, confirmed=3, goodIterative=4, qualitySize=5};
1207  if( t->quality(reco::TrackBase::loose)){ cout << "loose ";};
1208  if( t->quality(reco::TrackBase::tight)){ cout << "tight ";};
1209  if( t->quality(reco::TrackBase::highPurity)){ cout << "highPurity ";};
1210  if( t->quality(reco::TrackBase::confirmed)){ cout << "confirmed ";};
1211  if( t->quality(reco::TrackBase::goodIterative)){ cout << "goodIterative ";};
1212  cout << endl;
1213 
1214  TransientTrack tk = theB_->build(&(*t)); tk.setBeamSpot(vertexBeamSpot_);
1215  double ipsig=0;
1216  if (tk.stateAtBeamLine().isValid()){
1217  ipsig= tk.stateAtBeamLine().transverseImpactParameter().significance();
1218  }else{
1219  ipsig=-1;
1220  }
1221 
1222  cout << Form("pt=%8.3f phi=%6.3f eta=%6.3f z=%8.4f dz=%6.4f, ipsig=%6.1f",t->pt(), t->phi(), t->eta(), t->vz(), t->dzError(),ipsig) << endl;
1223 
1224 
1225  cout << Form(" found=%6d lost=%6d chi2/ndof=%10.3f ",t->found(), t->lost(),t->normalizedChi2())<<endl;
1226  const reco::HitPattern & p= t->hitPattern();
1227  cout << "subdet layers valid lost" << endl;
1228  cout << Form(" barrel %2d %2d %2d",p.pixelBarrelLayersWithMeasurement(),p.numberOfValidPixelBarrelHits(), p.numberOfLostPixelBarrelHits()) << endl;
1229  cout << Form(" fwd %2d %2d %2d",p.pixelEndcapLayersWithMeasurement(),p.numberOfValidPixelEndcapHits(), p.numberOfLostPixelEndcapHits()) << endl;
1230  cout << Form(" pixel %2d %2d %2d",p.pixelLayersWithMeasurement(), p.numberOfValidPixelHits(), p.numberOfLostPixelHits()) << endl;
1231  cout << Form(" tracker %2d %2d %2d",p.trackerLayersWithMeasurement(), p.numberOfValidTrackerHits(), p.numberOfLostTrackerHits()) << endl;
1232  cout << endl;
1233  const reco::HitPattern & pinner= t->trackerExpectedHitsInner();
1234  const reco::HitPattern & pouter= t->trackerExpectedHitsOuter();
1235  int ninner=pinner.numberOfHits();
1236  int nouter=pouter.numberOfHits();
1237 
1238  // double d0Error=t.d0Error();
1239  // double d0=t.dxy(myBeamSpot);
1240 
1241  //
1242  for(trackingRecHit_iterator hit=t->recHitsBegin(); hit!=t->recHitsEnd(); hit++){
1243  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1244  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1245  bool endcap = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1246  if (barrel){
1247  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1248  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1249  if (clust.isNonnull()) {
1250  cout << Form(" barrel cluster size=%2d charge=%6.1f wx=%2d wy=%2d, expected=%3.1f",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY(),1.+2./fabs(tan(t->theta()))) << endl;
1251  }
1252  }else if(endcap){
1253  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1254  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1255  if (clust.isNonnull()) {
1256  cout << Form(" endcap cluster size=%2d charge=%6.1f wx=%2d wy=%2d",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY()) << endl;
1257  }
1258  }
1259  }
1260  }
1261  cout << "hitpattern" << endl;
1262  for(int i=0; i<p.numberOfHits(); i++){ p.printHitPattern(i,std::cout); }
1263  cout << "expected inner " << ninner << endl;
1264  for(int i=0; i<pinner.numberOfHits(); i++){ pinner.printHitPattern(i,std::cout); }
1265  cout << "expected outer " << nouter << endl;
1266  for(int i=0; i<pouter.numberOfHits(); i++){ pouter.printHitPattern(i,std::cout); }
1267  }
1268 }
1269 
1270 namespace {
1271 
1272  bool recTrackLessZ(const reco::TransientTrack & tk1,
1273  const reco::TransientTrack & tk2)
1274  {
1275  return tk1.stateAtBeamLine().trackStateAtPCA().position().z() < tk2.stateAtBeamLine().trackStateAtPCA().position().z();
1276  }
1277 
1278 }
1279 
1280 
1281 
1282 
1283 void PrimaryVertexAnalyzer4PU::printPVTrks(const Handle<reco::TrackCollection> &recTrks,
1284  const Handle<reco::VertexCollection> &recVtxs,
1285  std::vector<SimPart>& tsim,
1286  std::vector<SimEvent>& simEvt,
1287  bool selectedOnly){
1288  // make a printout of the tracks selected for PV reconstructions, show matching MC tracks, too
1289 
1290  vector<TransientTrack> selTrks;
1291  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1292  t!=recTrks->end(); ++t){
1293  TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
1294  if( (!selectedOnly) || (selectedOnly && theTrackFilter(tt))){ selTrks.push_back(tt); }
1295  }
1296  if (selTrks.size()==0) return;
1297  stable_sort(selTrks.begin(), selTrks.end(), recTrackLessZ);
1298 
1299  // select tracks like for PV reconstruction and match them to sim tracks
1300  reco::TrackCollection selRecTrks;
1301 
1302  for(unsigned int i=0; i<selTrks.size(); i++){ selRecTrks.push_back(selTrks[i].track());}
1303  int* rectosim=supf(tsim, selRecTrks);
1304 
1305 
1306 
1307  // now dump in a format similar to the clusterizer
1308  cout << "printPVTrks" << endl;
1309  cout << "---- z +/- dz 1bet-m ip +/-dip pt phi eta";
1310  if((tsim.size()>0)||(simEvt.size()>0)) {cout << " type pdg zvtx zdca dca zvtx zdca dsz";}
1311  cout << endl;
1312 
1313  cout.precision(4);
1314  int isel=0;
1315  for(unsigned int i=0; i<selTrks.size(); i++){
1316  if (selectedOnly || (theTrackFilter(selTrks[i]))) {
1317  cout << setw (3)<< isel;
1318  isel++;
1319  }else{
1320  cout << " ";
1321  }
1322 
1323 
1324  // is this track in the tracklist of a recvtx ?
1325  int vmatch=-1;
1326  int iv=0;
1327  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1328  v!=recVtxs->end(); ++v){
1329  if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
1330  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
1331  const reco::Track & RTv=*(tv->get());
1332  if(selTrks[i].track().vz()==RTv.vz()) {vmatch=iv;}
1333  }
1334  iv++;
1335  }
1336 
1337  double tz=(selTrks[i].stateAtBeamLine().trackStateAtPCA()).position().z();
1338  double tantheta=tan((selTrks[i].stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1339  double tdz0= selTrks[i].track().dzError();
1340  double tdz2= pow(selTrks[i].track().dzError(),2)+ (pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2))/pow(tantheta,2);
1341 
1342  if(vmatch>-1){
1343  cout << "["<<setw(2)<<vmatch<<"]";
1344  }else{
1345  //int status=theTrackFilter.status(selTrks[i]);
1346  int status=0;
1347  if(status==0){
1348  cout <<" ";
1349  }else{
1350  if(status&0x1){cout << "i";}else{cout << ".";};
1351  if(status&0x2){cout << "c";}else{cout << ".";};
1352  if(status&0x4){cout << "h";}else{cout << ".";};
1353  if(status&0x8){cout << "q";}else{cout << ".";};
1354  }
1355  }
1356  cout << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tdz2);
1357 
1358 
1359  // track quality and hit information, see DataFormats/TrackReco/interface/HitPattern.h
1360  if(selTrks[i].track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";}
1361  if(selTrks[i].track().hitPattern().hasValidHitInFirstPixelBarrel()){cout <<"+";}else{cout << "-";}
1362  cout << setw(1) << selTrks[i].track().hitPattern().pixelBarrelLayersWithMeasurement();
1363  cout << setw(1) << selTrks[i].track().hitPattern().pixelEndcapLayersWithMeasurement();
1364  cout << setw(1) << hex << selTrks[i].track().hitPattern().trackerLayersWithMeasurement()-selTrks[i].track().hitPattern().pixelLayersWithMeasurement()<<dec;
1365  cout << "-" << setw(1)<<hex <<selTrks[i].track().trackerExpectedHitsOuter().numberOfHits() << dec;
1366 
1367 
1368  Measurement1D IP=selTrks[i].stateAtBeamLine().transverseImpactParameter();
1369  cout << setw (8) << IP.value() << "+/-" << setw (6) << IP.error();
1370  if(selTrks[i].track().ptError()<1){
1371  cout << " " << setw(8) << setprecision(2) << selTrks[i].track().pt()*selTrks[i].track().charge();
1372  }else{
1373  cout << " " << setw(7) << setprecision(1) << selTrks[i].track().pt()*selTrks[i].track().charge() << "*";
1374  //cout << "+/-" << setw(6)<< setprecision(2) << selTrks[i].track().ptError();
1375  }
1376  cout << " " << setw(5) << setprecision(2) << selTrks[i].track().phi()
1377  << " " << setw(5) << setprecision(2) << selTrks[i].track().eta() ;
1378 
1379 
1380 
1381  // print MC info, if available
1382  if(simEvt.size()>0){
1383  reco::Track t=selTrks[i].track();
1384  try{
1385  TrackingParticleRef tpr = z2tp_[t.vz()];
1386  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1387  //double vx=parentVertex->position().x(); // problems with tpr->vz()
1388  //double vy=parentVertex->position().y(); work in progress
1389  double vz=parentVertex->position().z();
1390  cout << " " << tpr->eventId().event();
1391  cout << " " << setw(5) << tpr->pdgId();
1392  cout << " " << setw(8) << setprecision(4) << vz;
1393  }catch (...){
1394  cout << " not matched";
1395  }
1396  }else{
1397  // cout << " " << rectosim[i];
1398  if(rectosim[i]>0){
1399  if(tsim[rectosim[i]].type==0){ cout << " prim " ;}else{cout << " sec ";}
1400  cout << " " << setw(5) << tsim[rectosim[i]].pdg;
1401  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zvtx;
1402  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zdcap;
1403  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].ddcap;
1404  double zvtxpull=(tz-tsim[rectosim[i]].zvtx)/sqrt(tdz2);
1405  cout << setw(6) << setprecision(1) << zvtxpull;
1406  double zdcapull=(tz-tsim[rectosim[i]].zdcap)/tdz0;
1407  cout << setw(6) << setprecision(1) << zdcapull;
1408  double dszpull=(selTrks[i].track().dsz()-tsim[rectosim[i]].par[4])/selTrks[i].track().dszError();
1409  cout << setw(6) << setprecision(1) << dszpull;
1410  }
1411  }
1412  cout << endl;
1413  }
1414  delete rectosim;
1415 }
1416 
1417 
1418 void PrimaryVertexAnalyzer4PU::matchRecTracksToVertex(simPrimaryVertex & pv,
1419  const std::vector<SimPart > & tsim,
1420  const edm::Handle<reco::TrackCollection> & recTrks)
1421 {
1422  // find all recTracks that belong to this simulated vertex (not just the ones that are used in
1423  // matching recVertex)
1424 
1425  std::cout << "dump rec tracks: " << std::endl;
1426  int irec=0;
1427  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1428  t!=recTrks->end(); ++t){
1429  reco::TrackBase::ParameterVector p = t->parameters();
1430  std::cout << irec++ << ") " << p << std::endl;
1431  }
1432 
1433  std::cout << "match sim tracks: " << std::endl;
1434  pv.matchedRecTrackIndex.clear();
1435  pv.nMatchedTracks=0;
1436  int isim=0;
1437  for(std::vector<SimPart>::const_iterator s=tsim.begin();
1438  s!=tsim.end(); ++s){
1439  std::cout << isim++ << " " << s->par;// << std::endl;
1440  int imatch=-1;
1441  int irec=0;
1442  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1443  t!=recTrks->end(); ++t){
1444  reco::TrackBase::ParameterVector p = t->parameters();
1445  if (match(s->par,p)){ imatch=irec; }
1446  irec++;
1447  }
1448  pv.matchedRecTrackIndex.push_back(imatch);
1449  if(imatch>-1){
1450  pv.nMatchedTracks++;
1451  std::cout << " matched to rec trk" << imatch << std::endl;
1452  }else{
1453  std::cout << " not matched" << std::endl;
1454  }
1455  }
1456 }
1457 /********************************************************************************************************/
1458 
1459 
1460 
1461 
1462 
1463 
1464 /********************************************************************************************************/
1465 
1466 void PrimaryVertexAnalyzer4PU::getTc(const std::vector<reco::TransientTrack>& tracks,
1467  double & Tc, double & chsq, double & dzmax, double & dztrim, double & m4m2){
1468  if (tracks.size()<2){ Tc=-1; chsq=-1; dzmax=-1; dztrim=-1; m4m2=-1; return;}
1469 
1470  double sumw=0, sumwz=0, sumww=0,sumwwz=0,sumwwzz=0;
1471  double zmin=1e10, zmin1=1e10, zmax1=-1e10, zmax=-1e10;
1472  double m4=0,m3=0,m2=0,m1=0,m0=0;
1473  for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
1474  double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1475  reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
1476  double z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
1477  double dz2= pow((*it).track().dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
1478  // t.dz2= pow((*it).track().dzError(),2) + pow(wxy0/tantheta,2) + 1./(1.+exp(pow(t.ip/t.dip,2)-pow(2.)))*pow(ip/tantheta,2);
1479  double w=1./dz2; // take p=1
1480  sumw += w;
1481  sumwz += w*z;
1482  sumwwz += w*w*z;;
1483  sumwwzz += w*w*z*z;
1484  sumww += w*w;
1485  m0 += w;
1486  m1 += w*z;
1487  m2 += w*z*z;
1488  m3 += w*z*z*z;
1489  m4 += w*z*z*z*z;
1490  if(dz2<pow(0.1,2)){
1491  if(z<zmin1){zmin1=z;} if(z<zmin){zmin1=zmin; zmin=z;}
1492  if(z>zmax1){zmax1=z;} if(z>zmax){zmax1=zmax; zmax=z;}
1493  }
1494  }
1495  double z=sumwz/sumw;
1496  double a=sumwwzz-2*z*sumwwz+z*z*sumww;
1497  double b=sumw;
1498  if(tracks.size()>1){
1499  chsq=(m2-m0*z*z)/(tracks.size()-1);
1500  Tc=2.*a/b;
1501  m4m2=sqrt((m4-4*m3*z+6*m2*z*z-3*m1*z*z*z+m0*z*z*z*z)/(m2-2*m1*z+z*z*m0));
1502  }else{
1503  chsq=0;
1504  Tc=0;
1505  m4m2=0;
1506  }
1507  dzmax=zmax-zmin;
1508  dztrim=zmax1-zmin1;// truncated
1509 }
1510 /********************************************************************************************************/
1511 
1512 
1513 
1514 
1515 
1516 /********************************************************************************************************/
1517 bool PrimaryVertexAnalyzer4PU::truthMatchedTrack( edm::RefToBase<reco::Track> track, TrackingParticleRef & tpr)
1518 
1519 /********************************************************************************************************/
1520 // for a reco track select the matching tracking particle, always use this function to make sure we
1521 // are consistent
1522 // to get the TrackingParticle form the TrackingParticleRef, use ->get();
1523 {
1524  double f=0;
1525  try{
1526  std::vector<std::pair<TrackingParticleRef, double> > tp = r2s_[track];
1527  for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = tp.begin();
1528  it != tp.end(); ++it) {
1529 
1530  if (it->second>f){
1531  tpr=it->first;
1532  f=it->second;
1533  }
1534  }
1535  } catch (Exception event) {
1536  // silly way of testing whether track is in r2s_
1537  }
1538 
1539  // sanity check on track parameters?
1540 
1541  return (f>0.5);
1542 }
1543 /********************************************************************************************************/
1544 
1545 
1546 
1547 
1548 
1549 
1550 /********************************************************************************************************/
1551 std::vector< edm::RefToBase<reco::Track> > PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks(
1552  const reco::Vertex& v
1553  )
1554 // for vertex v get a list of tracks for which truth matching is available
1555 /********************************************************************************************************/
1556 {
1557  std::vector< edm::RefToBase<reco::Track> > b;
1558  TrackingParticleRef tpr;
1559 
1560  for(trackit_t tv=v.tracks_begin(); tv!=v.tracks_end(); tv++){
1561  edm::RefToBase<reco::Track> track=*tv;
1562  if (truthMatchedTrack(track, tpr)){
1563  b.push_back(*tv);
1564  }
1565  }
1566 
1567 
1568  return b;
1569 }
1570 /********************************************************************************************************/
1571 
1572 
1573 
1574 
1575 
1576 /********************************************************************************************************/
1577 std::vector<PrimaryVertexAnalyzer4PU::SimEvent> PrimaryVertexAnalyzer4PU::getSimEvents
1578 (
1579  // const Event& iEvent, const EventSetup& iSetup,
1580  edm::Handle<TrackingParticleCollection> TPCollectionH,
1581  edm::Handle<TrackingVertexCollection> TVCollectionH,
1582  edm::Handle<View<Track> > trackCollectionH
1583  ){
1584 
1585  const TrackingParticleCollection* simTracks = TPCollectionH.product();
1586  const View<Track> tC = *(trackCollectionH.product());
1587 
1588 
1589  vector<SimEvent> simEvt;
1590  map<EncodedEventId, unsigned int> eventIdToEventMap;
1591  map<EncodedEventId, unsigned int>::iterator id;
1592  bool dumpTP=false;
1593  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1594 
1595  if( fabs(it->parentVertex().get()->position().z())>100.) continue; // skip funny entries @ -13900
1596 
1597  unsigned int event=0; //note, this is no longer the same as eventId().event()
1598  id=eventIdToEventMap.find(it->eventId());
1599  if (id==eventIdToEventMap.end()){
1600 
1601  // new event here
1602  SimEvent e;
1603  e.eventId=it->eventId();
1604  event=simEvt.size();
1605  const TrackingVertex *parentVertex= it->parentVertex().get();
1606  if(DEBUG_){
1607  cout << "getSimEvents: ";
1608  cout << it->eventId().bunchCrossing() << "," << it->eventId().event()
1609  << " z="<< it->vz() << " "
1610  << parentVertex->eventId().bunchCrossing() << "," <<parentVertex->eventId().event()
1611  << " z=" << parentVertex->position().z() << endl;
1612  }
1613  if (it->eventId()==parentVertex->eventId()){
1614  //e.x=it->vx(); e.y=it->vy(); e.z=it->vz();// wrong ???
1615  e.x=parentVertex->position().x();
1616  e.y=parentVertex->position().y();
1617  e.z=parentVertex->position().z();
1618  if(e.z<-100){
1619  dumpTP=true;
1620  }
1621  }else{
1622  e.x=0; e.y=0; e.z=-88.;
1623  }
1624  simEvt.push_back(e);
1625  eventIdToEventMap[e.eventId]=event;
1626  }else{
1627  event=id->second;
1628  }
1629 
1630 
1631  simEvt[event].tp.push_back(&(*it));
1632  if( (abs(it->eta())<2.5) && (it->charge()!=0) ){
1633  simEvt[event].sumpt2+=pow(it->pt(),2); // should keep track of decays ?
1634  simEvt[event].sumpt+=it->pt();
1635  }
1636  }
1637 
1638  if(dumpTP){
1639  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1640  std::cout << *it << std::endl;
1641  }
1642  }
1643 
1644 
1645  for(View<Track>::size_type i=0; i<tC.size(); ++i) {
1646  RefToBase<Track> track(trackCollectionH, i);
1647  TrackingParticleRef tpr;
1648  if( truthMatchedTrack(track,tpr)){
1649 
1650  if( eventIdToEventMap.find(tpr->eventId())==eventIdToEventMap.end() ){ cout << "Bug in getSimEvents" << endl; break; }
1651 
1652  z2tp_[track.get()->vz()]=tpr;
1653 
1654  unsigned int event=eventIdToEventMap[tpr->eventId()];
1655  double ipsig=0,ipdist=0;
1656  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1657  double vx=parentVertex->position().x(); // problems with tpr->vz()
1658  double vy=parentVertex->position().y();
1659  double vz=parentVertex->position().z();
1660  double d=sqrt(pow(simEvt[event].x-vx,2)+pow(simEvt[event].y-vy,2)+pow(simEvt[event].z-vz,2))*1.e4;
1661  ipdist=d;
1662  double dxy=track->dxy(vertexBeamSpot_.position());
1663  ipsig=dxy/track->dxyError();
1664 
1665 
1666  TransientTrack t = theB_->build(tC[i]);
1667  t.setBeamSpot(vertexBeamSpot_);
1668  if (theTrackFilter(t)){
1669  simEvt[event].tk.push_back(t);
1670  if(ipdist<5){simEvt[event].tkprim.push_back(t);}
1671  if(ipsig<5){simEvt[event].tkprimsel.push_back(t);}
1672  }
1673 
1674  }
1675  }
1676 
1677 
1678 
1679  AdaptiveVertexFitter theFitter;
1680  cout << "SimEvents " << simEvt.size() << endl;
1681  for(unsigned int i=0; i<simEvt.size(); i++){
1682 
1683  if(simEvt[i].tkprim.size()>0){
1684 
1685  getTc(simEvt[i].tkprimsel, simEvt[i].Tc, simEvt[i].chisq, simEvt[i].dzmax, simEvt[i].dztrim, simEvt[i].m4m2);
1686  TransientVertex v = theFitter.vertex(simEvt[i].tkprim, vertexBeamSpot_);
1687  if (v.isValid()){
1688  simEvt[i].xfit=v.position().x();
1689  simEvt[i].yfit=v.position().y();
1690  simEvt[i].zfit=v.position().z();
1691  // if (simEvt[i].z<-80.){simEvt[i].z=v.position().z();} // for now
1692  }
1693  }
1694 
1695 
1696  if(DEBUG_){
1697  cout << i <<" ) nTP=" << simEvt[i].tp.size()
1698  << " z=" << simEvt[i].z
1699  << " recTrks=" << simEvt[i].tk.size()
1700  << " recTrksPrim=" << simEvt[i].tkprim.size()
1701  << " zfit=" << simEvt[i].zfit
1702  << endl;
1703  }
1704  }
1705 
1706  return simEvt;
1707 }
1708 
1709 
1710 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1711  const Handle<HepMCProduct> evtMC)
1712 {
1713  if(DEBUG_){std::cout << "getSimPVs HepMC " << std::endl;}
1714 
1715  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1716  const HepMC::GenEvent* evt=evtMC->GetEvent();
1717  if (evt) {
1718 // std::cout << "process id " <<evt->signal_process_id()<<std::endl;
1719 // std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
1720 // evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
1721 // std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
1722 
1723 
1724  //int idx=0;
1725  for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
1726  vitr != evt->vertices_end(); ++vitr )
1727  { // loop for vertex ...
1728 
1729  HepMC::FourVector pos = (*vitr)->position();
1730  // if (pos.t()>0) { continue;} // skip secondary vertices, doesn't work for some samples
1731 
1732  if (fabs(pos.z())>1000) continue; // skip funny junk vertices
1733 
1734  bool hasMotherVertex=false;
1735  //std::cout << "mothers" << std::endl;
1736  for ( HepMC::GenVertex::particle_iterator
1737  mother = (*vitr)->particles_begin(HepMC::parents);
1738  mother != (*vitr)->particles_end(HepMC::parents);
1739  ++mother ) {
1740  HepMC::GenVertex * mv=(*mother)->production_vertex();
1741  if (mv) {hasMotherVertex=true;}
1742  //std::cout << "\t"; (*mother)->print();
1743  }
1744  if(hasMotherVertex) {continue;}
1745 
1746 
1747  // could be a new vertex, check all primaries found so far to avoid multiple entries
1748  const double mm=0.1;
1749  simPrimaryVertex sv(pos.x()*mm,pos.y()*mm,pos.z()*mm);
1750  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
1751  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1752  v0!=simpv.end(); v0++){
1753  if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
1754  vp=&(*v0);
1755  break;
1756  }
1757  }
1758  if(!vp){
1759  // this is a new vertex
1760  //std::cout << "this is a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;
1761  simpv.push_back(sv);
1762  vp=&simpv.back();
1763 // }else{
1764 // std::cout << "this is not a new vertex" << std::endl;
1765  }
1766 
1767 
1768  // store the gen vertex barcode with this simpv
1769  vp->genVertex.push_back((*vitr)->barcode());
1770 
1771 
1772  // collect final state descendants and sum up momenta etc
1773  for ( HepMC::GenVertex::particle_iterator
1774  daughter = (*vitr)->particles_begin(HepMC::descendants);
1775  daughter != (*vitr)->particles_end(HepMC::descendants);
1776  ++daughter ) {
1777  //std::cout << "checking daughter type " << (*daughter)->pdg_id() << " final :" <<isFinalstateParticle(*daughter) << std::endl;
1778  if (isFinalstateParticle(*daughter)){
1779  if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
1780  == vp->finalstateParticles.end()){
1781  vp->finalstateParticles.push_back((*daughter)->barcode());
1782  HepMC::FourVector m=(*daughter)->momentum();
1783  //std::cout << "adding particle to primary " << m.px()<< " " << m.py() << " " << m.pz() << std::endl;
1784  vp->ptot.setPx(vp->ptot.px()+m.px());
1785  vp->ptot.setPy(vp->ptot.py()+m.py());
1786  vp->ptot.setPz(vp->ptot.pz()+m.pz());
1787  vp->ptot.setE(vp->ptot.e()+m.e());
1788  vp->ptsq+=(m.perp())*(m.perp());
1789  // count relevant particles
1790  if ( (m.perp()>0.2) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
1791  vp->nGenTrk++;
1792  }
1793 
1794  hsimPV["rapidity"]->Fill(m.pseudoRapidity());
1795  if( (m.perp()>0.8) && isCharged( *daughter ) ){
1796  hsimPV["chRapidity"]->Fill(m.pseudoRapidity());
1797  }
1798  hsimPV["pt"]->Fill(m.perp());
1799  }//new final state particle for this vertex
1800  }//daughter is a final state particle
1801  }
1802 
1803 
1804  //idx++;
1805  }
1806  }
1807  if(verbose_){
1808  cout << "------- PrimaryVertexAnalyzer4PU simPVs -------" << endl;
1809  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1810  v0!=simpv.end(); v0++){
1811  cout << "z=" << v0->z
1812  << " px=" << v0->ptot.px()
1813  << " py=" << v0->ptot.py()
1814  << " pz=" << v0->ptot.pz()
1815  << " pt2="<< v0->ptsq
1816  << endl;
1817  }
1818  cout << "-----------------------------------------------" << endl;
1819  }
1820  return simpv;
1821 }
1822 
1823 
1824 
1825 
1826 
1827 
1828 
1829 
1830 /* get sim pv from TrackingParticles/TrackingVertex */
1831 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1832  const edm::Handle<TrackingVertexCollection> tVC
1833  )
1834 {
1835  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1836  //std::cout <<"number of vertices " << tVC->size() << std::endl;
1837 
1838  if(DEBUG_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;}
1839 
1840  for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) {
1841 
1842  if(DEBUG_){
1843  std::cout << (v->eventId()).event() << v -> position() << v->g4Vertices().size() <<" " <<v->genVertices().size() << " t=" <<v->position().t()*1.e12 <<" ==0:" <<(v->position().t()>0) << std::endl;
1844  for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){
1845  std::cout << *gv << std::endl;
1846  }
1847  std::cout << "----------" << std::endl;
1848  }
1849 
1850  // bool hasMotherVertex=false;
1851  if ((unsigned int) v->eventId().event()<simpv.size()) continue;
1852  //if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work)
1853  if (fabs(v->position().z())>1000) continue; // skip funny junk vertices
1854 
1855  // could be a new vertex, check all primaries found so far to avoid multiple entries
1856  const double mm=1.0; // for tracking vertices
1857  simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm);
1858 
1859  //cout << "sv: " << (v->eventId()).event() << endl;
1860  sv.eventId=v->eventId();
1861 
1862  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){
1863  //cout <<((**iTrack).eventId()).event() << " "; // an iterator of Refs, dereference twice. Cool eyh?
1864  sv.eventId=(**iTrack).eventId();
1865  }
1866  //cout <<endl;
1867  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
1868  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1869  v0!=simpv.end(); v0++){
1870  if( (sv.eventId==v0->eventId) && (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
1871  vp=&(*v0);
1872  break;
1873  }
1874  }
1875  if(!vp){
1876  // this is a new vertex
1877  if(DEBUG_){std::cout << "this is a new vertex " << sv.eventId.event() << " " << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
1878  simpv.push_back(sv);
1879  vp=&simpv.back();
1880  }else{
1881  if(DEBUG_){std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
1882  }
1883 
1884 
1885  // Loop over daughter track(s)
1886  if(DEBUG_){
1887  for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) {
1888  std::cout << " Daughter momentum: " << (*(*iTP)).momentum();
1889  std::cout << " Daughter type " << (*(*iTP)).pdgId();
1890  std::cout << std::endl;
1891  }
1892  }
1893  }
1894  if(DEBUG_){
1895  cout << "------- PrimaryVertexAnalyzer4PU simPVs from TrackingVertices -------" << endl;
1896  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1897  v0!=simpv.end(); v0++){
1898  cout << "z=" << v0->z << " event=" << v0->eventId.event() << endl;
1899  }
1900  cout << "-----------------------------------------------" << endl;
1901  }
1902  return simpv;
1903 }
1904 
1905 
1906 
1907 
1908 
1909 
1910 // ------------ method called to produce the data ------------
1911 void
1912 PrimaryVertexAnalyzer4PU::analyze(const Event& iEvent, const EventSetup& iSetup)
1913 {
1914 
1915  std::vector<simPrimaryVertex> simpv; // a list of primary MC vertices
1916  std::vector<SimPart> tsim;
1917  std::string mcproduct="generator"; // starting with 3_1_0 pre something
1918 
1919  eventcounter_++;
1920  run_ = iEvent.id().run();
1921  luminosityBlock_ = iEvent.luminosityBlock();
1922  event_ = iEvent.id().event();
1923  bunchCrossing_ = iEvent.bunchCrossing();
1924  orbitNumber_ = iEvent.orbitNumber();
1925 
1926  dumpThisEvent_ = false;
1927 
1928 
1929 
1930  if(verbose_){
1931  std::cout << endl
1932  << "PrimaryVertexAnalyzer4PU::analyze event counter=" << eventcounter_
1933  << " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
1934  << " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
1935  //<< " selected = " << good
1936  << std::endl;
1937  }
1938 
1939 
1940  try{
1941  iSetup.getData(pdt_);
1942  }catch(const Exception&){
1943  std::cout << "Some problem occurred with the particle data table. This may not work !" <<std::endl;
1944  }
1945 
1946  Handle<reco::VertexCollection> recVtxs;
1947  bool bnoBS=iEvent.getByLabel("offlinePrimaryVertices", recVtxs);
1948 
1949  Handle<reco::VertexCollection> recVtxsBS;
1950  bool bBS=iEvent.getByLabel("offlinePrimaryVerticesWithBS", recVtxsBS);
1951 
1952  Handle<reco::VertexCollection> recVtxsDA;
1953  bool bDA=iEvent.getByLabel("offlinePrimaryVerticesDA", recVtxsDA);
1954 
1955 // Handle<reco::VertexCollection> recVtxsPIX;
1956 // bool bPIX=iEvent.getByLabel("pixelVertices", recVtxsPIX);
1957 // bPIX=false;
1958 
1959 // Handle<reco::VertexCollection> recVtxsMVF;
1960 // bool bMVF=iEvent.getByLabel("offlinePrimaryVerticesMVF", recVtxsMVF);
1961 
1962  Handle<reco::TrackCollection> recTrks;
1963  iEvent.getByLabel(recoTrackProducer_, recTrks);
1964 
1965 
1966  int nhighpurity=0, ntot=0;
1967  for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
1968  ntot++;
1969  if(t->quality(reco::TrackBase::highPurity)) nhighpurity++;
1970  }
1971  if(ntot>10) hnoBS["highpurityTrackFraction"]->Fill(float(nhighpurity)/float(ntot));
1972  if((recoTrackProducer_=="generalTracks")&&(nhighpurity<0.25*ntot)){
1973  if(verbose_){ cout << "rejected, " << nhighpurity << " highpurity out of " << ntot << " total tracks "<< endl<< endl;}
1974  return;
1975  }
1976 
1977 
1978 
1979 
1980 
1981  if(iEvent.getByType(recoBeamSpotHandle_)){
1982  vertexBeamSpot_= *recoBeamSpotHandle_;
1983  wxy2_=pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2);
1984  Fill(hsimPV, "xbeam",vertexBeamSpot_.x0()); Fill(hsimPV, "wxbeam",vertexBeamSpot_.BeamWidthX());
1985  Fill(hsimPV, "ybeam",vertexBeamSpot_.y0()); Fill(hsimPV, "wybeam",vertexBeamSpot_.BeamWidthY());
1986  Fill(hsimPV, "zbeam",vertexBeamSpot_.z0());// Fill("wzbeam",vertexBeamSpot_.BeamWidthZ());
1987  }else{
1988  cout << " beamspot not found, using dummy " << endl;
1989  vertexBeamSpot_=reco::BeamSpot();// dummy
1990  }
1991 
1992 
1993  if(bnoBS){
1994  if((recVtxs->begin()->isValid())&&(recVtxs->begin()->ndof()>1)&&(recVtxs->begin()->ndof()>(0.0*recVtxs->begin()->tracksSize()))){ // was 5 and 0.2
1995  Fill(hnoBS,"xrecBeamvsdxXBS",recVtxs->begin()->xError(),recVtxs->begin()->x()-vertexBeamSpot_.x0());
1996  Fill(hnoBS,"yrecBeamvsdyXBS",recVtxs->begin()->yError(),recVtxs->begin()->y()-vertexBeamSpot_.y0());
1997 
1998  if(printXBS_) {
1999  cout << Form("XBS %10d %5d %10d %4d %lu %5.2f %8f %8f %8f %8f %8f %8f",
2000  run_,luminosityBlock_,event_,bunchCrossing_,
2001  (unsigned long)(recVtxs->begin()->tracksSize()), recVtxs->begin()->ndof(),
2002  recVtxs->begin()->x(), recVtxs->begin()->xError(),
2003  recVtxs->begin()->y(), recVtxs->begin()->yError(),
2004  recVtxs->begin()->z(), recVtxs->begin()->zError()
2005  ) << endl;
2006  }
2007 
2008  }
2009  }
2010 
2011 
2012  // for the associator
2013  Handle<View<Track> > trackCollectionH;
2014  iEvent.getByLabel(recoTrackProducer_,trackCollectionH);
2015 
2016  Handle<HepMCProduct> evtMC;
2017 
2018  Handle<SimVertexContainer> simVtxs;
2019  iEvent.getByLabel( simG4_, simVtxs);
2020 
2021  Handle<SimTrackContainer> simTrks;
2022  iEvent.getByLabel( simG4_, simTrks);
2023 
2024 
2025 
2026 
2027 
2028  edm::Handle<TrackingParticleCollection> TPCollectionH ;
2029  edm::Handle<TrackingVertexCollection> TVCollectionH ;
2030  bool gotTP=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionH);
2031  bool gotTV=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TVCollectionH);
2032 
2033 
2034  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
2035  fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
2036 
2037 
2038  vector<SimEvent> simEvt;
2039  if (gotTP && gotTV ){
2040 
2041  edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
2042  iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
2043  associatorByHits_ = (TrackAssociatorBase *) theHitsAssociator.product();
2044  r2s_ = associatorByHits_->associateRecoToSim (trackCollectionH,TPCollectionH, &iEvent );
2045  //simEvt=getSimEvents(iEvent, TPCollectionH, TVCollectionH, trackCollectionH);
2046  simEvt=getSimEvents(TPCollectionH, TVCollectionH, trackCollectionH);
2047 
2048  if (simEvt.size()==0){
2049  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2050  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2051  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2052  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2053  cout << " !!!!!!!!!!!!!!!!!!!!!! got TrackingParticles but no simEvt !!!!!!!!!!!!!!!!!" << endl;
2054  cout << " dumping Tracking particles " << endl;
2055  const TrackingParticleCollection* simTracks = TPCollectionH.product();
2056  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
2057  cout << *it << endl;
2058  }
2059  cout << " dumping Tracking Vertices " << endl;
2060  const TrackingVertexCollection* tpVtxs = TVCollectionH.product();
2061  for(TrackingVertexCollection::const_iterator iv=tpVtxs->begin(); iv!=tpVtxs->end(); iv++){
2062  cout << *iv << endl;
2063  }
2064  if(iEvent.getByLabel(mcproduct,evtMC)){
2065  cout << "dumping simTracks" << endl;
2066  for(SimTrackContainer::const_iterator t=simTrks->begin(); t!=simTrks->end(); ++t){
2067  cout << *t << endl;
2068  }
2069  cout << "dumping simVertexes" << endl;
2070  for(SimVertexContainer::const_iterator vv=simVtxs->begin();
2071  vv!=simVtxs->end();
2072  ++vv){
2073  cout << *vv << endl;
2074  }
2075  }else{
2076  cout << "no hepMC" << endl;
2077  }
2078  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2079 
2080  const HepMC::GenEvent* evt=evtMC->GetEvent();
2081  if(evt){
2082  std::cout << "process id " <<evt->signal_process_id()<<std::endl;
2083  std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
2084  evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
2085  std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
2086  evt->print();
2087  }else{
2088  cout << "no event in HepMC" << endl;
2089  }
2090  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2091 
2092  }
2093  }
2094 
2095 
2096 
2097 
2098  if(gotTV){
2099 
2100  if(verbose_){
2101  cout << "Found Tracking Vertices " << endl;
2102  }
2103  simpv=getSimPVs(TVCollectionH);
2104 
2105 
2106  }else if(iEvent.getByLabel(mcproduct,evtMC)){
2107 
2108  simpv=getSimPVs(evtMC);
2109 
2110  if(verbose_){
2111  cout << "Using HepMCProduct " << endl;
2112  std::cout << "simtrks " << simTrks->size() << std::endl;
2113  }
2114  tsim = PrimaryVertexAnalyzer4PU::getSimTrkParameters(simTrks, simVtxs, simUnit_);
2115 
2116  }else{
2117  // if(verbose_({cout << "No MC info at all" << endl;}
2118  }
2119 
2120 
2121 
2122 
2123  // get the beam spot from the appropriate dummy vertex, if available
2124  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2125  v!=recVtxs->end(); ++v){
2126  if ((v->ndof()==-3) && (v->chi2()==0) ){
2127  myBeamSpot=math::XYZPoint(v->x(), v->y(), v->z());
2128  }
2129  }
2130 
2131 
2132 
2133 
2134  hsimPV["nsimvtx"]->Fill(simpv.size());
2135  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
2136  vsim!=simpv.end(); vsim++){
2137  if(doMatching_){
2138  matchRecTracksToVertex(*vsim, tsim, recTrks); // hepmc, based on track parameters
2139  }
2140 
2141  hsimPV["nbsimtksinvtx"]->Fill(vsim->nGenTrk);
2142  hsimPV["xsim"]->Fill(vsim->x*simUnit_);
2143  hsimPV["ysim"]->Fill(vsim->y*simUnit_);
2144  hsimPV["zsim"]->Fill(vsim->z*simUnit_);
2145  hsimPV["xsim1"]->Fill(vsim->x*simUnit_);
2146  hsimPV["ysim1"]->Fill(vsim->y*simUnit_);
2147  hsimPV["zsim1"]->Fill(vsim->z*simUnit_);
2148  Fill(hsimPV,"xsim2",vsim->x*simUnit_,vsim==simpv.begin());
2149  Fill(hsimPV,"ysim2",vsim->y*simUnit_,vsim==simpv.begin());
2150  Fill(hsimPV,"zsim2",vsim->z*simUnit_,vsim==simpv.begin());
2151  hsimPV["xsim2"]->Fill(vsim->x*simUnit_);
2152  hsimPV["ysim2"]->Fill(vsim->y*simUnit_);
2153  hsimPV["zsim2"]->Fill(vsim->z*simUnit_);
2154  hsimPV["xsim3"]->Fill(vsim->x*simUnit_);
2155  hsimPV["ysim3"]->Fill(vsim->y*simUnit_);
2156  hsimPV["zsim3"]->Fill(vsim->z*simUnit_);
2157  hsimPV["xsimb"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
2158  hsimPV["ysimb"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
2159  hsimPV["zsimb"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
2160  hsimPV["xsimb1"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
2161  hsimPV["ysimb1"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
2162  hsimPV["zsimb1"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
2163  }
2164 
2165 
2166 
2167 
2168 
2169  if(bnoBS){
2170  analyzeVertexCollection(hnoBS, recVtxs, recTrks, simpv,"noBS");
2171  analyzeVertexCollectionTP(hnoBS, recVtxs, recTrks, simEvt,"noBS");
2172  }
2173  if(bBS){
2174  analyzeVertexCollection(hBS, recVtxsBS, recTrks, simpv,"BS");
2175  analyzeVertexCollectionTP(hBS, recVtxsBS, recTrks, simEvt,"BS");
2176  }
2177  if(bDA){
2178  analyzeVertexCollection(hDA, recVtxsDA, recTrks, simpv,"DA");
2179  analyzeVertexCollectionTP(hDA, recVtxsDA, recTrks, simEvt,"DA");
2180  }
2181 // if(bPIX){
2182 // analyzeVertexCollection(hPIX, recVtxsPIX, recTrks, simpv,"PIX");
2183 // analyzeVertexCollectionTP(hPIX, recVtxsPIX, recTrks, simEvt,"PIX");
2184 // }
2185 
2186 
2187 
2188  if((dumpThisEvent_&& (dumpcounter_<100)) ||(verbose_ && (eventcounter_<ndump_))){
2189  cout << endl << "Event dump" << endl
2190  << "event counter=" << eventcounter_
2191  << " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
2192  << " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
2193  << std::endl;
2194  dumpcounter_++;
2195 
2196  //evtMC->GetEvent()->print();
2197  //printRecTrks(recTrks); // very verbose !!
2198 
2199 // if (bPIX) printRecVtxs(recVtxsPIX,"pixel vertices");
2200  if (bnoBS) {printRecVtxs(recVtxs,"Offline without Beamspot");}
2201  if (bnoBS && (!bDA)){ printPVTrks(recTrks, recVtxs, tsim, simEvt, false);}
2202  if (bBS) printRecVtxs(recVtxsBS,"Offline with Beamspot");
2203  if (bDA) {
2204  printRecVtxs(recVtxsDA,"Offline DA");
2205  printPVTrks(recTrks, recVtxsDA, tsim, simEvt, false);
2206  }
2207  if (dumpcounter_<2){cout << "beamspot " << vertexBeamSpot_ << endl;}
2208  }
2209 
2210  if(verbose_){
2211  std::cout << std::endl;
2212  }
2213 }
2214 
2215 namespace {
2216 bool lt(const std::pair<double,unsigned int>& a,const std::pair<double,unsigned int>& b ){
2217  return a.first<b.first;
2218 }
2219 }
2220 
2221 /***************************************************************************************/
2222 void PrimaryVertexAnalyzer4PU::printEventSummary(std::map<std::string, TH1*> & h,
2223  const edm::Handle<reco::VertexCollection> recVtxs,
2224  const edm::Handle<reco::TrackCollection> recTrks,
2225  vector<SimEvent> & simEvt,
2226  const string message){
2227  // make a readable summary of the vertex finding if the TrackingParticles are availabe
2228  if (simEvt.size()==0){return;}
2229 
2230 
2231  // sort vertices in z ... for nicer printout
2232 
2233  vector< pair<double,unsigned int> > zrecv;
2234  for(unsigned int idx=0; idx<recVtxs->size(); idx++){
2235  if ( (recVtxs->at(idx).ndof()<0) || (recVtxs->at(idx).chi2()<=0) ) continue; // skip clusters
2236  zrecv.push_back( make_pair(recVtxs->at(idx).z(),idx) );
2237  }
2238  stable_sort(zrecv.begin(),zrecv.end(),lt);
2239 
2240  // same for simulated vertices
2241  vector< pair<double,unsigned int> > zsimv;
2242  for(unsigned int idx=0; idx<simEvt.size(); idx++){
2243  zsimv.push_back(make_pair(simEvt[idx].z, idx));
2244  }
2245  stable_sort(zsimv.begin(), zsimv.end(),lt);
2246 
2247 
2248 
2249 
2250  cout << "---------------------------" << endl;
2251  cout << "event counter = " << eventcounter_ << " " << message << endl;
2252  cout << "---------------------------" << endl;
2253  cout << " z[cm] rec --> ";
2254  cout.precision(4);
2255  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2256  cout << setw(7) << fixed << itrec->first;
2257  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2258  }
2259  cout << endl;
2260  cout << " ";
2261  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2262  cout << setw(7) << fixed << recVtxs->at(itrec->second).tracksSize();
2263  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2264  }
2265  cout << " rec tracks" << endl;
2266  cout << " ";
2267  map<unsigned int, int> truthMatchedVertexTracks;
2268  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2269  truthMatchedVertexTracks[itrec->second]=getTruthMatchedVertexTracks(recVtxs->at(itrec->second)).size();
2270  cout << setw(7) << fixed << truthMatchedVertexTracks[itrec->second];
2271  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2272  }
2273  cout << " truth matched " << endl;
2274 
2275  cout << "sim ------- trk prim ----" << endl;
2276 
2277 
2278 
2279  map<unsigned int, unsigned int> rvmatch; // reco vertex matched to sim vertex (sim to rec)
2280  map<unsigned int, double > nmatch; // highest number of truth-matched tracks of ev found in a recvtx
2281  map<unsigned int, double > purity; // highest purity of a rec vtx (i.e. highest number of tracks from the same simvtx)
2282  map<unsigned int, double > wpurity; // same for the sum of weights
2283 
2284  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2285  purity[itrec->second]=0.;
2286  wpurity[itrec->second]=0.;
2287  }
2288 
2289  for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
2290  SimEvent* ev =&(simEvt[itsim->second]);
2291 
2292 
2293  cout.precision(4);
2294  if (itsim->second==0){
2295  cout << setw(8) << fixed << ev->z << ")*" << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
2296  }else{
2297  cout << setw(8) << fixed << ev->z << ") " << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
2298  }
2299 
2300  nmatch[itsim->second]=0; // highest number of truth-matched tracks of ev found in a recvtx
2301  double matchpurity=0,matchwpurity=0;
2302 
2303  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2304  const reco::Vertex *v = &(recVtxs->at(itrec->second));
2305 
2306  // count tracks found in both, sim and rec
2307  double n=0,wt=0;
2308  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2309  const reco::Track& RTe=te->track();
2310  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2311  const reco::Track & RTv=*(tv->get());
2312  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2313  }
2314  }
2315  cout << setw(7) << int(n)<< " ";
2316 
2317  if (n > nmatch[itsim->second]){
2318  nmatch[itsim->second]=n;
2319  rvmatch[itsim->second]=itrec->second;
2320  matchpurity=n/truthMatchedVertexTracks[itrec->second];
2321  matchwpurity=wt/truthMatchedVertexTracks[itrec->second];
2322  }
2323 
2324  if(n > purity[itrec->second]){
2325  purity[itrec->second]=n;
2326  }
2327 
2328  if(wt > wpurity[itrec->second]){
2329  wpurity[itrec->second]=wt;
2330  }
2331 
2332  }// end of reco vertex loop
2333 
2334  cout << " | ";
2335  if (nmatch[itsim->second]>0 ){
2336  if(matchpurity>0.5){
2337  cout << "found ";
2338  }else{
2339  cout << "merged ";
2340  }
2341  cout << " max eff. = " << setw(8) << nmatch[itsim->second]/ev->tk.size() << " p=" << matchpurity << " w=" << matchwpurity << endl;
2342  }else{
2343  if(ev->tk.size()==0){
2344  cout << " invisible" << endl;
2345  }else if (ev->tk.size()==1){
2346  cout << "single track " << endl;
2347  }else{
2348  cout << "lost " << endl;
2349  }
2350  }
2351  }
2352  cout << "---------------------------" << endl;
2353 
2354  // the purity of the reconstructed vertex
2355  cout << " purity ";
2356  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2357  cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
2358  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2359  }
2360  cout << endl;
2361 
2362 // // classification of reconstructed vertex fake/real
2363 // cout << " ";
2364 // for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2365 // cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
2366 // if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2367 // }
2368 // cout << endl;
2369  cout << "---------------------------" << endl;
2370 
2371 
2372 
2373 
2374  // list problematic tracks
2375  for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
2376  SimEvent* ev =&(simEvt[itsim->second]);
2377 
2378  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2379  const reco::Track& RTe=te->track();
2380 
2381  int ivassign=-1; // will become the index of the vertex to which a track was assigned
2382 
2383  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2384  const reco::Vertex *v = &(recVtxs->at(itrec->second));
2385 
2386  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2387  const reco::Track & RTv=*(tv->get());
2388  if(RTe.vz()==RTv.vz()) {ivassign=itrec->second;}
2389  }
2390  }
2391  double tantheta=tan((te->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
2392  reco::BeamSpot beamspot=(te->stateAtBeamLine()).beamSpot();
2393  //double z=(te->stateAtBeamLine().trackStateAtPCA()).position().z();
2394  double dz2= pow(RTe.dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
2395 
2396  if(ivassign==(int)rvmatch[itsim->second]){
2397  Fill(h,"correctlyassigned",RTe.eta(),RTe.pt());
2398  Fill(h,"ptcat",RTe.pt());
2399  Fill(h,"etacat",RTe.eta());
2400  Fill(h,"phicat",RTe.phi());
2401  Fill(h,"dzcat",sqrt(dz2));
2402  }else{
2403  Fill(h,"misassigned",RTe.eta(),RTe.pt());
2404  Fill(h,"ptmis",RTe.pt());
2405  Fill(h,"etamis",RTe.eta());
2406  Fill(h,"phimis",RTe.phi());
2407  Fill(h,"dzmis",sqrt(dz2));
2408  cout << "vertex " << setw(8) << fixed << ev->z;
2409 
2410  if (ivassign<0){
2411  cout << " track lost ";
2412  // for some clusterizers there shouldn't be any lost tracks,
2413  // are there differences in the track selection?
2414  }else{
2415  cout << " track misassigned " << setw(8) << fixed << recVtxs->at(ivassign).z();
2416  }
2417 
2418  cout << " track z=" << setw(8) << fixed << RTe.vz() << "+/-" << RTe.dzError() << " pt=" << setw(8) << fixed<< RTe.pt() << " eta=" << setw(8) << fixed << RTe.eta()<< " sel=" <<theTrackFilter(*te);
2419 
2420  //
2421  //cout << " ztrack=" << te->track().vz();
2422  TrackingParticleRef tpr = z2tp_[te->track().vz()];
2423  double zparent=tpr->parentVertex().get()->position().z();
2424  if(zparent==ev->z) {
2425  cout << " prim";
2426  }else{
2427  cout << " sec ";
2428  }
2429  cout << " id=" << tpr->pdgId();
2430  cout << endl;
2431 
2432  //
2433  }
2434  }// next simvertex-track
2435 
2436  }//next simvertex
2437 
2438  cout << "---------------------------" << endl;
2439 
2440 }
2441 /***************************************************************************************/
2442 
2443 
2444 
2445 
2446 /***************************************************************************************/
2447 void PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP(std::map<std::string, TH1*> & h,
2448  const edm::Handle<reco::VertexCollection> recVtxs,
2449  const edm::Handle<reco::TrackCollection> recTrks,
2450  vector<SimEvent> & simEvt,
2451  const string message){
2452 
2453  // cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP size=" << simEvt.size() << endl;
2454  if(simEvt.size()==0)return;
2455 
2456  printEventSummary(h, recVtxs,recTrks,simEvt,message);
2457 
2458  //const int iSignal=0;
2459  EncodedEventId iSignal=simEvt[0].eventId;
2460  Fill(h,"npu0",simEvt.size());
2461 
2462 
2463  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2464  Fill(h,"Tc", ev->Tc, ev==simEvt.begin());
2465  Fill(h,"Chisq", ev->chisq, ev==simEvt.begin());
2466  if(ev->chisq>0)Fill(h,"logChisq", log(ev->chisq), ev==simEvt.begin());
2467  Fill(h,"dzmax", ev->dzmax, ev==simEvt.begin());
2468  Fill(h,"dztrim",ev->dztrim,ev==simEvt.begin());
2469  Fill(h,"m4m2", ev->m4m2, ev==simEvt.begin());
2470  if(ev->Tc>0){ Fill(h,"logTc",log(ev->Tc)/log(10.),ev==simEvt.begin());}
2471 
2472 
2473  for(vector<SimEvent>::iterator ev2=ev+1; ev2!=simEvt.end(); ev2++){
2474  vector<TransientTrack> xt;
2475  if((ev->tkprimsel.size()>0)&&(ev2->tkprimsel.size()>0)&&(ev->tkprimsel.size()+ev2->tkprimsel.size())>1){
2476  xt.insert (xt.end() ,ev->tkprimsel.begin(),ev->tkprimsel.end());
2477  xt.insert (xt.end() ,ev2->tkprimsel.begin(),ev2->tkprimsel.end());
2478  double xTc,xChsq,xDzmax,xDztrim,xm4m2;
2479  getTc(xt, xTc, xChsq, xDzmax, xDztrim,xm4m2);
2480  if(xTc>0){
2481  Fill(h,"xTc",xTc,ev==simEvt.begin());
2482  Fill(h,"logxTc", log(xTc)/log(10),ev==simEvt.begin());
2483  Fill(h,"xChisq", xChsq,ev==simEvt.begin());
2484  if(xChsq>0){Fill(h,"logxChisq", log(xChsq),ev==simEvt.begin());};
2485  Fill(h,"xdzmax", xDzmax,ev==simEvt.begin());
2486  Fill(h,"xdztrim", xDztrim,ev==simEvt.begin());
2487  Fill(h,"xm4m2", xm4m2,ev==simEvt.begin());
2488 
2489  }
2490  }
2491  }
2492  }
2493 
2494  // --------------------------------------- count actual rec vtxs ----------------------
2495  int nrecvtxs=0;//, nrecvtxs1=0, nrecvtxs2=0;
2496  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2497  v!=recVtxs->end(); ++v){
2498  if ( (v->isFake()) || (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
2499  nrecvtxs++;
2500  }
2501 
2502  // --------------------------------------- fill the track assignment matrix ----------------------
2503  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2504  ev->ntInRecVz.clear(); // just in case
2505  ev->zmatch=-99.;
2506  ev->nmatch=0;
2507  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2508  v!=recVtxs->end(); ++v){
2509  double n=0, wt=0;
2510  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2511  const reco::Track& RTe=te->track();
2512  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2513  const reco::Track & RTv=*(tv->get());
2514  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2515  }
2516  }
2517  ev->ntInRecVz[v->z()]=n;
2518  if (n > ev->nmatch){ ev->nmatch=n; ev->zmatch=v->z(); ev->zmatch=v->z(); }
2519  }
2520  }
2521 
2522  // call a vertex a fake if for every sim vertex there is another recvertex containing more tracks
2523  // from that sim vertex than the current recvertex
2524  double nfake=0;
2525  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2526  v!=recVtxs->end(); ++v){
2527  bool matched=false;
2528  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2529  if ((ev->nmatch>0)&&(ev->zmatch==v->z())){
2530  matched=true;
2531  }
2532  }
2533  if(!matched && !v->isFake()) {
2534  nfake++;
2535  cout << " fake rec vertex at z=" << v->z() << endl;
2536  // some histograms of fake vertex properties here
2537  Fill(h,"unmatchedVtxZ",v->z());
2538  Fill(h,"unmatchedVtxNdof",v->ndof());
2539  }
2540  }
2541  if(nrecvtxs>0){
2542  Fill(h,"unmatchedVtx",nfake);
2543  Fill(h,"unmatchedVtxFrac",nfake/nrecvtxs);
2544  }
2545 
2546  // --------------------------------------- match rec to sim ---------------------------------------
2547  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2548  v!=recVtxs->end(); ++v){
2549 
2550  if ( (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
2551  double nmatch=-1; // highest number of tracks in recvtx v found in any event
2552  EncodedEventId evmatch;
2553  bool matchFound=false;
2554  double nmatchvtx=0; // number of simvtcs contributing to recvtx v
2555 
2556  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2557 
2558  double n=0; // number of tracks that are in both, the recvtx v and the event ev
2559  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2560 
2561  const reco::Track& RTe=te->track();
2562  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2563  const reco::Track & RTv=*(tv->get());
2564  if(RTe.vz()==RTv.vz()){ n++;}
2565  }
2566  }
2567 
2568  // find the best match in terms of the highest number of tracks
2569  // from a simvertex in this rec vertex
2570  if (n > nmatch){
2571  nmatch=n;
2572  evmatch=ev->eventId;
2573  matchFound=true;
2574  }
2575  if(n>0){
2576  nmatchvtx++;
2577  }
2578  }
2579 
2580  double nmatchany=getTruthMatchedVertexTracks(*v).size();
2581  if (matchFound && (nmatchany>0)){
2582  // highest number of tracks in recvtx matched to (the same) sim vertex
2583  // purity := -----------------------------------------------------------------
2584  // number of truth matched tracks in this recvtx
2585  double purity =nmatch/nmatchany;
2586  Fill(h,"recmatchPurity",purity);
2587  if(v==recVtxs->begin()){
2588  Fill(h,"recmatchPurityTag",purity, (bool)(evmatch==iSignal));
2589  }else{
2590  Fill(h,"recmatchPuritynoTag",purity,(bool)(evmatch==iSignal));
2591  }
2592  }
2593  Fill(h,"recmatchvtxs",nmatchvtx);
2594  if(v==recVtxs->begin()){
2595  Fill(h,"recmatchvtxsTag",nmatchvtx);
2596  }else{
2597  Fill(h,"recmatchvtxsnoTag",nmatchvtx);
2598  }
2599 
2600 
2601 
2602  } // recvtx loop
2603  Fill(h,"nrecv",nrecvtxs);
2604 
2605 
2606  // --------------------------------------- match sim to rec ---------------------------------------
2607 
2608  int npu1=0, npu2=0;
2609 
2610  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2611 
2612  if(ev->tk.size()>0) npu1++;
2613  if(ev->tk.size()>1) npu2++;
2614 
2615  bool isSignal= ev->eventId==iSignal;
2616 
2617  Fill(h,"nRecTrkInSimVtx",(double) ev->tk.size(),isSignal);
2618  Fill(h,"nPrimRecTrkInSimVtx",(double) ev->tkprim.size(),isSignal);
2619  Fill(h,"sumpt2rec",sqrt(ev->sumpt2rec),isSignal);
2620  Fill(h,"sumpt2",sqrt(ev->sumpt2),isSignal);
2621  Fill(h,"sumpt",sqrt(ev->sumpt),isSignal);
2622 
2623  double nRecVWithTrk=0; // vertices with tracks from this simvertex
2624  double nmatch=0, ntmatch=0, zmatch=-99;
2625 
2626  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2627  v!=recVtxs->end(); ++v){
2628  if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
2629  // count tracks found in both, sim and rec
2630  double n=0, wt=0;
2631  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2632  const reco::Track& RTe=te->track();
2633  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2634  const reco::Track & RTv=*(tv->get());
2635  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2636  }
2637  }
2638 
2639  if(n>0){ nRecVWithTrk++; }
2640  if (n > nmatch){
2641  nmatch=n; ntmatch=v->tracksSize(); zmatch=v->position().z();
2642  }
2643 
2644  }// end of reco vertex loop
2645 
2646 
2647  // nmatch is the highest number of tracks from this sim vertex found in a single reco-vertex
2648  if(ev->tk.size()>0){ Fill(h,"trkAssignmentEfficiency", nmatch/ev->tk.size(), isSignal); };
2649  if(ev->tkprim.size()>0){ Fill(h,"primtrkAssignmentEfficiency", nmatch/ev->tkprim.size(), isSignal); };
2650 
2651  // matched efficiency = efficiency for finding a reco vertex with > 50% of the simvertexs reconstructed tracks
2652 
2653  double ntsim=ev->tk.size(); // may be better to use the number of primary tracks here ?
2654  double matchpurity=nmatch/ntmatch;
2655 
2656  if(ntsim>0){
2657 
2658  Fill(h,"matchVtxFraction",nmatch/ntsim,isSignal);
2659  if(nmatch/ntsim>=0.5){
2660  Fill(h,"matchVtxEfficiency",1.,isSignal);
2661  if(ntsim>1){Fill(h,"matchVtxEfficiency2",1.,isSignal);}
2662  if(matchpurity>0.5){Fill(h,"matchVtxEfficiency5",1.,isSignal);}
2663  }else{
2664  Fill(h,"matchVtxEfficiency",0.,isSignal);
2665  if(ntsim>1){Fill(h,"matchVtxEfficiency2",0.,isSignal);}
2666  Fill(h,"matchVtxEfficiency5",0.,isSignal); // no (matchpurity>5) here !!
2667  if(isSignal){
2668  cout << "Signal vertex not matched " << message << " event=" << eventcounter_ << " nmatch=" << nmatch << " ntsim=" << ntsim << endl;
2669  }
2670  }
2671  } // ntsim >0
2672 
2673 
2674  if(zmatch>-99){
2675  Fill(h,"matchVtxZ",zmatch-ev->z);
2676  Fill(h,"matchVtxZ",zmatch-ev->z,isSignal);
2677  Fill(h,"matchVtxZCum",fabs(zmatch-ev->z));
2678  Fill(h,"matchVtxZCum",fabs(zmatch-ev->z),isSignal);
2679  }else{
2680  Fill(h,"matchVtxZCum",1.0);
2681  Fill(h,"matchVtxZCum",1.0,isSignal);
2682  }
2683  if(fabs(zmatch-ev->z)<zmatch_){
2684  Fill(h,"matchVtxEfficiencyZ",1.,isSignal);
2685  }else{
2686  Fill(h,"matchVtxEfficiencyZ",0.,isSignal);
2687  }
2688 
2689  if(ntsim>0) Fill(h, "matchVtxEfficiencyZ1", fabs(zmatch-ev->z)<zmatch_ , isSignal);
2690  if(ntsim>1) Fill(h, "matchVtxEfficiencyZ2", fabs(zmatch-ev->z)<zmatch_ , isSignal);
2691 
2692 
2693  Fill(h,"vtxMultiplicity",nRecVWithTrk,isSignal);
2694 
2695  // efficiency vs number of tracks, use your favorite definition of efficiency here
2696  //if(nmatch>=0.5*ntmatch){ // purity
2697  if(fabs(zmatch-ev->z)<zmatch_){ // zmatch
2698  Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),1.);
2699  if(isSignal){
2700  Fill(h,"vtxFindingEfficiencyVsNtrkSignal",ev->tk.size(),1.);
2701  }else{
2702  Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
2703  }
2704  }else{
2705  Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),0.);
2706  if(isSignal){
2707  Fill(h,"vtxFindingEfficiencyVsNtrkSignal",(double) ev->tk.size(),1.);
2708  }else{
2709  Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
2710  }
2711  }
2712 
2713 
2714  }
2715 
2716  Fill(h,"npu1",npu1);
2717  Fill(h,"npu2",npu2);
2718 
2719  Fill(h,"nrecvsnpu",npu1,float(nrecvtxs));
2720  Fill(h,"nrecvsnpu2",npu2,float(nrecvtxs));
2721 
2722  // --------------------------------------- sim-signal vs rec-tag ---------------------------------------
2723  SimEvent* ev=&(simEvt[0]);
2724  const reco::Vertex* v=&(*recVtxs->begin());
2725 
2726  double n=0;
2727  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2728  const reco::Track& RTe=te->track();
2729  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2730  const reco::Track & RTv=*(tv->get());
2731  if(RTe.vz()==RTv.vz()) {n++;}
2732  }
2733  }
2734 
2735  cout << "Number of tracks in reco tagvtx " << v->tracksSize() << endl;
2736  cout << "Number of selected tracks in sim event vtx " << ev->tk.size() << " (prim=" << ev->tkprim.size() << ")"<<endl;
2737  cout << "Number of tracks in both " << n << endl;
2738  double ntruthmatched=getTruthMatchedVertexTracks(*v).size();
2739  if (ntruthmatched>0){
2740  cout << "TrackPurity = "<< n/ntruthmatched <<endl;
2741  Fill(h,"TagVtxTrkPurity",n/ntruthmatched);
2742  }
2743  if (ev->tk.size()>0){
2744  cout << "TrackEfficiency = "<< n/ev->tk.size() <<endl;
2745  Fill(h,"TagVtxTrkEfficiency",n/ev->tk.size());
2746  }
2747 }
2748 
2749 /***************************************************************************************/
2750 
2751 
2752 
2753 
2754 
2755 
2756 /***************************************************************************************/
2757 
2758 void PrimaryVertexAnalyzer4PU::analyzeVertexCollection(std::map<std::string, TH1*> & h,
2759  const Handle<reco::VertexCollection> recVtxs,
2760  const Handle<reco::TrackCollection> recTrks,
2761  std::vector<simPrimaryVertex> & simpv,
2762  const std::string message)
2763 {
2764  //cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollection (HepMC), simpvs=" << simpv.size() << endl;
2765  int nrectrks=recTrks->size();
2766  int nrecvtxs=recVtxs->size();
2767  int nseltrks=-1;
2768  reco::TrackCollection selTrks; // selected tracks
2769  reco::TrackCollection lostTrks; // selected but lost tracks (part of dropped clusters)
2770 
2771  // extract dummy vertices representing clusters
2772  reco::VertexCollection clusters;
2773  reco::Vertex allSelected;
2774  double cpufit=0;
2775  double cpuclu=0;
2776  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2777  v!=recVtxs->end(); ++v){
2778  if ( (fabs(v->ndof()+3.)<0.0001) && (v->chi2()<=0) ){
2779  // this dummy vertex is for the full event
2780  allSelected=(*v);
2781  nseltrks=(allSelected.tracksSize());
2782  nrecvtxs--;
2783  cpuclu=-v->chi2();
2784  continue;
2785  }else if( (fabs(v->ndof()+2.)<0.0001) && (v->chi2()==0) ){
2786  // this is a cluster, not a vertex
2787  clusters.push_back(*v);
2788  Fill(h,"cpuvsntrk",(double) v->tracksSize(),fabs(v->y()));
2789  cpufit+=fabs(v->y());
2790  Fill(h,"nclutrkall",(double) v->tracksSize());
2791  Fill(h,"selstat",v->x());
2792  //Fill(h,"nclutrkvtx",);// see below
2793  nrecvtxs--;
2794  }
2795  }
2796  Fill(h,"cpuclu",cpuclu);
2797  Fill(h,"cpufit",cpufit);
2798  Fill(h,"cpucluvsntrk",nrectrks, cpuclu);
2799 
2800 
2801 
2802  if(simpv.size()>0){//this is mc
2803  double dsimrecx=0.;
2804  double dsimrecy=0.;//0.0011;
2805  double dsimrecz=0.;//0.0012;
2806 
2807  // vertex matching and efficiency bookkeeping
2808  int nsimtrk=0;
2809  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
2810  vsim!=simpv.end(); vsim++){
2811 
2812 
2813  nsimtrk+=vsim->nGenTrk;
2814  // look for a matching reconstructed vertex
2815  vsim->recVtx=NULL;
2816  vsim->cluster=-1;
2817 
2818  for(reco::VertexCollection::const_iterator vrec=recVtxs->begin(); vrec!=recVtxs->end(); ++vrec){
2819 
2820  if( vrec->isFake() ) {
2821  continue; // skip fake vertices (=beamspot)
2822  cout << "fake vertex" << endl;
2823  }
2824 
2825  if( vrec->ndof()<0. )continue; // skip dummy clusters, if any
2826  // if ( matchVertex(*vsim,*vrec) ){
2827 
2828  // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
2829  if( ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z-dsimrecz)>fabs(vrec->z()-vsim->z-dsimrecz)))
2830  || (!vsim->recVtx) )
2831  {
2832  vsim->recVtx=&(*vrec);
2833 
2834  // find the corresponding cluster
2835  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
2836  if( fabs(clusters[iclu].position().z()-vrec->position().z()) < 0.001 ){
2837  vsim->cluster=iclu;
2838  vsim->nclutrk=clusters[iclu].position().y();
2839  }
2840  }
2841  }
2842 
2843  // the following only works in MC samples without pile-up
2844  if ((simpv.size()==1) && ( fabs(vsim->recVtx->position().z()-vsim->z-dsimrecz)>zmatch_ )){
2845  // now we have a recvertex without a matching simvertex, I would call this fake
2846  // however, the G4 info does not contain pile-up
2847  Fill(h,"fakeVtxZ",vrec->z());
2848  if (vrec->ndof()>=0.5) Fill(h,"fakeVtxZNdofgt05",vrec->z());
2849  if (vrec->ndof()>=2.0) Fill(h,"fakeVtxZNdofgt2",vrec->z());
2850  Fill(h,"fakeVtxNdof",vrec->ndof());
2851  //Fill(h,"fakeVtxNdof1",vrec->ndof());
2852  Fill(h,"fakeVtxNtrk",vrec->tracksSize());
2853  if(vrec->tracksSize()==2){ Fill(h,"fake2trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
2854  if(vrec->tracksSize()==3){ Fill(h,"fake3trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
2855  if(vrec->tracksSize()==4){ Fill(h,"fake4trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
2856  if(vrec->tracksSize()==5){ Fill(h,"fake5trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
2857  }
2858  }
2859 
2860 
2861  Fill(h,"nsimtrk",float(nsimtrk));
2862  Fill(h,"nsimtrk",float(nsimtrk),vsim==simpv.begin());
2863  Fill(h,"nrecsimtrk",float(vsim->nMatchedTracks));
2864  Fill(h,"nrecnosimtrk",float(nsimtrk-vsim->nMatchedTracks));
2865 
2866  // histogram properties of matched vertices
2867  if (vsim->recVtx && ( fabs(vsim->recVtx->z()-vsim->z*simUnit_)<zmatch_ )){
2868 
2869  if(verbose_){std::cout <<"primary matched " << message << " " << setw(8) << setprecision(4) << vsim->x << " " << vsim->y << " " << vsim->z << std:: endl;}
2870  Fill(h,"matchedVtxNdof", vsim->recVtx->ndof());
2871  // residuals an pulls with respect to simulated vertex
2872  Fill(h,"resx", vsim->recVtx->x()-vsim->x*simUnit_ );
2873  Fill(h,"resy", vsim->recVtx->y()-vsim->y*simUnit_ );
2874  Fill(h,"resz", vsim->recVtx->z()-vsim->z*simUnit_ );
2875  Fill(h,"resz10", vsim->recVtx->z()-vsim->z*simUnit_ );
2876  Fill(h,"pullx", (vsim->recVtx->x()-vsim->x*simUnit_)/vsim->recVtx->xError() );
2877  Fill(h,"pully", (vsim->recVtx->y()-vsim->y*simUnit_)/vsim->recVtx->yError() );
2878  Fill(h,"pullz", (vsim->recVtx->z()-vsim->z*simUnit_)/vsim->recVtx->zError() );
2879  Fill(h,"resxr", vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx);
2880  Fill(h,"resyr", vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy );
2881  Fill(h,"reszr", vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz);
2882  Fill(h,"pullxr", (vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx)/vsim->recVtx->xError() );
2883  Fill(h,"pullyr", (vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy)/vsim->recVtx->yError() );
2884  Fill(h,"pullzr", (vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz)/vsim->recVtx->zError() );
2885 
2886 
2887 
2888  // efficiency with zmatch within 500 um (or whatever zmatch is)
2889  Fill(h,"eff", 1.);
2890  if(simpv.size()==1){
2891  if (vsim->recVtx==&(*recVtxs->begin())){
2892  Fill(h,"efftag", 1.);
2893  }else{
2894  Fill(h,"efftag", 0.);
2895  cout << "signal vertex not tagged " << message << " " << eventcounter_ << endl;
2896  // call XXClusterizerInZ.vertices(seltrks,3)
2897 
2898  }
2899  }
2900 
2901  Fill(h,"effvsptsq",vsim->ptsq,1.);
2902  Fill(h,"effvsnsimtrk",vsim->nGenTrk,1.);
2903  Fill(h,"effvsnrectrk",nrectrks,1.);
2904  Fill(h,"effvsnseltrk",nseltrks,1.);
2905  Fill(h,"effvsz",vsim->z*simUnit_,1.);
2906  Fill(h,"effvsz2",vsim->z*simUnit_,1.);
2907  Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,1.);
2908 
2909 
2910  }else{ // no matching rec vertex found for this simvertex
2911 
2912  bool plapper=verbose_ && vsim->nGenTrk;
2913  if(plapper){
2914  // be quiet about invisble vertices
2915  std::cout << "primary not found " << message << " " << eventcounter_ << " x=" <<vsim->x*simUnit_ << " y=" << vsim->y*simUnit_ << " z=" << vsim->z*simUnit_ << " nGenTrk=" << vsim->nGenTrk << std::endl;
2916  }
2917  int mistype=0;
2918  if (vsim->recVtx){
2919  if(plapper){
2920  std::cout << "nearest recvertex at " << vsim->recVtx->z() << " dz=" << vsim->recVtx->z()-vsim->z*simUnit_ << std::endl;
2921  }
2922 
2923  if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.2 ){
2924  Fill(h,"effvsz2",vsim->z*simUnit_,1.);
2925  }
2926 
2927  if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.5 ){
2928  if(verbose_){std::cout << "type 1, lousy z vertex" << std::endl;}
2929  Fill(h,"zlost1", vsim->z*simUnit_,1.);
2930  mistype=1;
2931  }else{
2932  if(plapper){std::cout << "type 2a no vertex anywhere near" << std::endl;}
2933  mistype=2;
2934  }
2935  }else{// no recVtx at all
2936  mistype=2;
2937  if(plapper){std::cout << "type 2b, no vertex at all" << std::endl;}
2938  }
2939 
2940  if(mistype==2){
2941  int selstat=-3;
2942  // no matching vertex found, is there a cluster?
2943  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
2944  if( fabs(clusters[iclu].position().z()-vsim->z*simUnit_) < 0.1 ){
2945  selstat=int(clusters[iclu].position().x()+0.1);
2946  if(verbose_){std::cout << "matching cluster found with selstat=" << clusters[iclu].position().x() << std::endl;}
2947  }
2948  }
2949  if (selstat==0){
2950  if(plapper){std::cout << "vertex rejected (distance to beam)" << std::endl;}
2951  Fill(h,"zlost3", vsim->z*simUnit_,1.);
2952  }else if(selstat==-1){
2953  if(plapper) {std::cout << "vertex invalid" << std::endl;}
2954  Fill(h,"zlost4", vsim->z*simUnit_,1.);
2955  }else if(selstat==1){
2956  if(plapper){std::cout << "vertex accepted, this cannot be right!!!!!!!!!!" << std::endl;}
2957  }else if(selstat==-2){
2958  if(plapper){std::cout << "dont know what this means !!!!!!!!!!" << std::endl;}
2959  }else if(selstat==-3){
2960  if(plapper){std::cout << "no matching cluster found " << std::endl;}
2961  Fill(h,"zlost2", vsim->z*simUnit_,1.);
2962  }else{
2963  if(plapper){std::cout << "dont know what this means either !!!!!!!!!!" << selstat << std::endl;}
2964  }
2965  }//
2966 
2967 
2968  Fill(h,"eff", 0.);
2969  if(simpv.size()==1){ Fill(h,"efftag", 0.); }
2970 
2971  Fill(h,"effvsptsq",vsim->ptsq,0.);
2972  Fill(h,"effvsnsimtrk",float(vsim->nGenTrk),0.);
2973  Fill(h,"effvsnrectrk",nrectrks,0.);
2974  Fill(h,"effvsnseltrk",nseltrks,0.);
2975  Fill(h,"effvsz",vsim->z*simUnit_,0.);
2976  Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,0.);
2977 
2978  } // no recvertex for this simvertex
2979 
2980  }
2981 
2982  // end of sim/rec matching
2983 
2984 
2985  // purity of event vertex tags
2986  if (recVtxs->size()>0){
2987  Double_t dz=(*recVtxs->begin()).z() - (*simpv.begin()).z*simUnit_;
2988  Fill(h,"zdistancetag",dz);
2989  Fill(h,"abszdistancetag",fabs(dz));
2990  if( fabs(dz)<zmatch_){
2991  Fill(h,"puritytag",1.);
2992  }else{
2993  // bad tag: the true primary was more than 500 um (or zmatch) away from the tagged primary
2994  Fill(h,"puritytag",0.);
2995  }
2996  }
2997 
2998  }else{
2999  //if(verbose_) cout << "PrimaryVertexAnalyzer4PU::analyzeVertexCollection: simPV is empty!" << endl;
3000  }
3001 
3002 
3003  //******* the following code does not require MC and will/should work for data **********
3004 
3005 
3006  Fill(h,"bunchCrossing",bunchCrossing_);
3007  if(recTrks->size()>0) Fill(h,"bunchCrossingLogNtk",bunchCrossing_,log(recTrks->size())/log(10.));
3008 
3009  // ----------------- reconstructed tracks ------------------------
3010  // the list of selected tracks can only be correct if the selection parameterset and track collection
3011  // is the same that was used for the reconstruction
3012 
3013  int nt=0;
3014  for(reco::TrackCollection::const_iterator t=recTrks->begin();
3015  t!=recTrks->end(); ++t){
3016  if((recVtxs->size()>0) && (recVtxs->begin()->isValid())){
3017  fillTrackHistos(h,"all",*t,&(*recVtxs->begin()));
3018  }else{
3019  fillTrackHistos(h,"all",*t);
3020  }
3021  if(recTrks->size()>100) fillTrackHistos(h,"M",*t);
3022 
3023 
3024 
3025  TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
3026  if (theTrackFilter(tt)){
3027  selTrks.push_back(*t);
3028  fillTrackHistos(h,"sel",*t);
3029  int foundinvtx=0;
3030  int nvtemp=-1;
3031  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3032  v!=recVtxs->end(); ++v){
3033  nvtemp++;
3034  if(( v->isFake()) || (v->ndof()<-2) ) break;
3035  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++ ){
3036  if( ((**tv).vz()==t->vz()&&((**tv).phi()==t->phi())) ) {
3037  foundinvtx++;
3038  }
3039  }
3040 
3041  }
3042  if(foundinvtx==0){
3043  fillTrackHistos(h,"sellost",*t);
3044  }else if(foundinvtx>1){
3045  cout << "hmmmm " << foundinvtx << endl;
3046  }
3047  }
3048  nt++;
3049  }
3050 
3051 
3052  if (nseltrks<0){
3053  nseltrks=selTrks.size();
3054  }else if( ! (nseltrks==(int)selTrks.size()) ){
3055  std::cout << "Warning: inconsistent track selection !" << std::endl;
3056  }
3057 
3058 
3059 
3060  // fill track histograms of vertex tracks
3061  int nrec=0, nrec0=0, nrec8=0, nrec2=0, nrec4=0;
3062  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3063  v!=recVtxs->end(); ++v){
3064 
3065  if (! (v->isFake()) && v->ndof()>0 && v->chi2()>0 ){
3066  nrec++;
3067  if (v->ndof()>0) nrec0++;
3068  if (v->ndof()>8) nrec8++;
3069  if (v->ndof()>2) nrec2++;
3070  if (v->ndof()>4) nrec4++;
3071  for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){
3072  if(v==recVtxs->begin()){
3073  fillTrackHistos(h,"tagged",**t, &(*v));
3074  }else{
3075  fillTrackHistos(h,"untagged",**t, &(*v));
3076  }
3077 
3078  Float_t wt=v->trackWeight(*t);
3079  //dumpHitInfo(**t); cout << " w=" << wt << endl;
3080  Fill(h,"trackWt",wt);
3081  if(wt>0.5){
3082  fillTrackHistos(h,"wgt05",**t, &(*v));
3083  if(v->ndof()>4) fillTrackHistos(h,"ndof4",**t, &(*v));
3084  }else{
3085  fillTrackHistos(h,"wlt05",**t, &(*v));
3086  }
3087  }
3088  }
3089  }
3090 
3091 
3092  // bachelor tracks (only available through clusters right now)
3093  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3094  if (clusters[iclu].tracksSize()==1){
3095  for(trackit_t t = clusters[iclu].tracks_begin();
3096  t!=clusters[iclu].tracks_end(); t++){
3097  fillTrackHistos(h,"bachelor",**t);
3098  }
3099  }
3100  }
3101 
3102 
3103  // ----------------- reconstructed vertices ------------------------
3104 
3105  // event
3106  Fill(h,"szRecVtx",recVtxs->size());
3107  Fill(h,"nclu",clusters.size());
3108  Fill(h,"nseltrk",nseltrks);
3109  Fill(h,"nrectrk",nrectrks);
3110  Fill(h,"nrecvtx",nrec);
3111  Fill(h,"nrecvtx2",nrec2);
3112  Fill(h,"nrecvtx4",nrec4);
3113  Fill(h,"nrecvtx8",nrec8);
3114 
3115  if(nrec>0){
3116  Fill(h,"eff0vsntrec",nrectrks,1.);
3117  Fill(h,"eff0vsntsel",nseltrks,1.);
3118  }else{
3119  Fill(h,"eff0vsntrec",nrectrks,0.);
3120  Fill(h,"eff0vsntsel",nseltrks,0.);
3121  if((nseltrks>1)&&(verbose_)){
3122  cout << Form("PrimaryVertexAnalyzer4PU: %s may have lost a vertex %10d %10d %4d / %4d ",message.c_str(),run_, event_, nrectrks,nseltrks) << endl;
3123  dumpThisEvent_=true;
3124  }
3125  }
3126  if(nrec0>0) { Fill(h,"eff0ndof0vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof0vsntsel",nseltrks,0.);}
3127  if(nrec2>0) { Fill(h,"eff0ndof2vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof2vsntsel",nseltrks,0.);}
3128  if(nrec4>0) { Fill(h,"eff0ndof4vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof4vsntsel",nseltrks,0.);}
3129  if(nrec8>0) { Fill(h,"eff0ndof8vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof8vsntsel",nseltrks,0.);}
3130 
3131  if((nrec>1)&&(DEBUG_)) {
3132  cout << "multivertex event" << endl;
3133  dumpThisEvent_=true;
3134  }
3135 
3136  if((nrectrks>10)&&(nseltrks<3)){
3137  cout << "small fraction of selected tracks " << endl;
3138  dumpThisEvent_=true;
3139  }
3140 
3141  // properties of events without a vertex
3142  if((nrec==0)||(recVtxs->begin()->isFake())){
3143  Fill(h,"nrectrk0vtx",nrectrks);
3144  Fill(h,"nseltrk0vtx",nseltrks);
3145  Fill(h,"nclu0vtx",clusters.size());
3146  }
3147 
3148 
3149  // properties of (valid) vertices
3150  double ndof2=-10,ndof1=-10, zndof1=0, zndof2=0;
3151  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3152  v!=recVtxs->end(); ++v){
3153  if(v->isFake()){ Fill(h,"isFake",1.);}else{ Fill(h,"isFake",0.);}
3154  if(v->isFake()||((v->ndof()<0)&&(v->ndof()>-3))){ Fill(h,"isFake1",1.);}else{ Fill(h,"isFake1",0.);}
3155 
3156  if((v->isFake())||(v->ndof()<0)) continue;
3157  if(v->ndof()>ndof1){ ndof2=ndof1; zndof2=zndof1; ndof1=v->ndof(); zndof1=v->position().z();}
3158  else if(v->ndof()>ndof2){ ndof2=v->ndof(); zndof2=v->position().z();}
3159 
3160 
3161  // some special histogram for two track vertices
3162  if(v->tracksSize()==2){
3163  const TrackBaseRef& t1= *(v->tracks_begin());
3164  const TrackBaseRef& t2=*(v->tracks_begin()+1);
3165  bool os=(t1->charge()*t2->charge()<0);
3166  double dphi=t1->phi()-t2->phi(); if (dphi<0) dphi+=2*M_PI;
3167  double m12=sqrt(pow( sqrt(pow(0.139,2)+pow( t1->p(),2)) +sqrt(pow(0.139,2)+pow( t2->p(),2)) ,2)
3168  -pow(t1->px()+t2->px(),2)
3169  -pow(t1->py()+t2->py(),2)
3170  -pow(t1->pz()+t2->pz(),2)
3171  );
3172  if(os){
3173  Fill(h,"2trkdetaOS",t1->eta()-t2->eta());
3174  Fill(h,"2trkmassOS",m12);
3175  }else{
3176  Fill(h,"2trkdetaSS",t1->eta()-t2->eta());
3177  Fill(h,"2trkmassSS",m12);
3178  }
3179  Fill(h,"2trkdphi",dphi);
3180  Fill(h,"2trkseta",t1->eta()+t2->eta());
3181  if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurl",t1->eta()+t2->eta());
3182  if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurl",dphi);
3183  // fill separately for extra vertices
3184  if(v!=recVtxs->begin()){
3185  if(os){
3186  Fill(h,"2trkdetaOSPU",t1->eta()-t2->eta());
3187  Fill(h,"2trkmassOSPU",m12);
3188  }else{
3189  Fill(h,"2trkdetaSSPU",t1->eta()-t2->eta());
3190  Fill(h,"2trkmassSSPU",m12);
3191  }
3192  Fill(h,"2trkdphiPU",dphi);
3193  Fill(h,"2trksetaPU",t1->eta()+t2->eta());
3194  if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurlPU",t1->eta()+t2->eta());
3195  if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurlPU",dphi);
3196  }
3197  }
3198 
3199 
3200  Fill(h,"trkchi2vsndof",v->ndof(),v->chi2());
3201  if(v->ndof()>0){ Fill(h,"trkchi2overndof",v->chi2()/v->ndof()); }
3202  if(v->tracksSize()==2){ Fill(h,"2trkchi2vsndof",v->ndof(),v->chi2()); }
3203  if(v->tracksSize()==3){ Fill(h,"3trkchi2vsndof",v->ndof(),v->chi2()); }
3204  if(v->tracksSize()==4){ Fill(h,"4trkchi2vsndof",v->ndof(),v->chi2()); }
3205  if(v->tracksSize()==5){ Fill(h,"5trkchi2vsndof",v->ndof(),v->chi2()); }
3206 
3207  Fill(h,"nbtksinvtx",v->tracksSize());
3208  Fill(h,"nbtksinvtx2",v->tracksSize());
3209  Fill(h,"vtxchi2",v->chi2());
3210  Fill(h,"vtxndf",v->ndof());
3211  Fill(h,"vtxprob",ChiSquaredProbability(v->chi2() ,v->ndof()));
3212  Fill(h,"vtxndfvsntk",v->tracksSize(), v->ndof());
3213  Fill(h,"vtxndfoverntk",v->ndof()/v->tracksSize());
3214  Fill(h,"vtxndf2overntk",(v->ndof()+2)/v->tracksSize());
3215  Fill(h,"zrecvsnt",v->position().z(),float(nrectrks));
3216  if(nrectrks>100){
3217  Fill(h,"zrecNt100",v->position().z());
3218  }
3219 
3220  if(v->ndof()>2.0){ // enter only vertices that really contain tracks
3221  Fill(h,"xrec",v->position().x());
3222  Fill(h,"yrec",v->position().y());
3223  Fill(h,"zrec",v->position().z());
3224  Fill(h,"xrec1",v->position().x());
3225  Fill(h,"yrec1",v->position().y());
3226  Fill(h,"zrec1",v->position().z());
3227  Fill(h,"xrec2",v->position().x());
3228  Fill(h,"yrec2",v->position().y());
3229  Fill(h,"zrec2",v->position().z());
3230  Fill(h,"xrec3",v->position().x());
3231  Fill(h,"yrec3",v->position().y());
3232  Fill(h,"zrec3",v->position().z());
3233  Fill(h,"xrecb",v->position().x()-vertexBeamSpot_.x0());
3234  Fill(h,"yrecb",v->position().y()-vertexBeamSpot_.y0());
3235  Fill(h,"zrecb",v->position().z()-vertexBeamSpot_.z0());
3236  Fill(h,"xrecBeam",v->position().x()-vertexBeamSpot_.x0());
3237  Fill(h,"yrecBeam",v->position().y()-vertexBeamSpot_.y0());
3238  Fill(h,"zrecBeam",v->position().z()-vertexBeamSpot_.z0());
3239  Fill(h,"xrecBeamPull",(v->position().x()-vertexBeamSpot_.x0())/sqrt(pow(v->xError(),2)+pow(vertexBeamSpot_.BeamWidthX(),2)));
3240  Fill(h,"yrecBeamPull",(v->position().y()-vertexBeamSpot_.y0())/sqrt(pow(v->yError(),2)+pow(vertexBeamSpot_.BeamWidthY(),2)));
3241 
3242  Fill(h,"xrecBeamvsdx",v->xError(),v->position().x()-vertexBeamSpot_.x0());
3243  Fill(h,"yrecBeamvsdy",v->yError(),v->position().y()-vertexBeamSpot_.y0());
3244  Fill(h,"xrecBeamvsdxR2",v->position().x()-vertexBeamSpot_.x0(),v->xError());
3245  Fill(h,"yrecBeamvsdyR2",v->position().y()-vertexBeamSpot_.y0(),v->yError());
3246  Fill(h,"xrecBeam2vsdx2prof",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
3247  Fill(h,"yrecBeam2vsdy2prof",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
3248  Fill(h,"xrecBeamvsdx2",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
3249  Fill(h,"yrecBeamvsdy2",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
3250  Fill(h,"xrecBeamvsz",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
3251  Fill(h,"yrecBeamvsz",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
3252  Fill(h,"xrecBeamvszprof",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
3253  Fill(h,"yrecBeamvszprof",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
3254  Fill(h,"xrecBeamvsdxprof",v->xError(),v->position().x()-vertexBeamSpot_.x0());
3255  Fill(h,"yrecBeamvsdyprof",v->yError(),v->position().y()-vertexBeamSpot_.y0());
3256 
3257 
3258  Fill(h,"errx",v->xError());
3259  Fill(h,"erry",v->yError());
3260  Fill(h,"errz",v->zError());
3261  double vxx=v->covariance(0,0);
3262  double vyy=v->covariance(1,1);
3263  double vxy=v->covariance(1,0);
3264  double dv=0.25*(vxx+vyy)*(vxx+vyy-(vxx*vyy-vxy*vxy));
3265  if(dv>0){
3266  double l1=0.5*(vxx+vyy)+sqrt(dv);
3267  Fill(h,"err1",sqrt(l1));
3268  double l2=sqrt(0.5*(vxx+vyy)-sqrt(dv));
3269  if(l2>0) Fill(h,"err2",sqrt(l2));
3270  }
3271 
3272 
3273  // look at the tagged vertex separately
3274  if (v==recVtxs->begin()){
3275  Fill(h,"nbtksinvtxTag",v->tracksSize());
3276  Fill(h,"nbtksinvtxTag2",v->tracksSize());
3277  Fill(h,"xrectag",v->position().x());
3278  Fill(h,"yrectag",v->position().y());
3279  Fill(h,"zrectag",v->position().z());
3280  }else{
3281  Fill(h,"nbtksinvtxPU",v->tracksSize());
3282  Fill(h,"nbtksinvtxPU2",v->tracksSize());
3283  }
3284 
3285  // vertex resolution vs number of tracks
3286  Fill(h,"xresvsntrk",v->tracksSize(),v->xError());
3287  Fill(h,"yresvsntrk",v->tracksSize(),v->yError());
3288  Fill(h,"zresvsntrk",v->tracksSize(),v->zError());
3289 
3290  }
3291 
3292  // cluster properties (if available)
3293  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3294  if( fabs(clusters[iclu].position().z()-v->position().z()) < 0.0001 ){
3295  Fill(h,"nclutrkvtx",clusters[iclu].tracksSize());
3296  }
3297  }
3298 
3299 
3300 
3301  // properties of (valid) neighbour vertices
3302  reco::VertexCollection::const_iterator v1=v; v1++;
3303  for(; v1!=recVtxs->end(); ++v1){
3304  if((v1->isFake())||(v1->ndof()<0)) continue;
3305  Fill(h,"zdiffrec",v->position().z()-v1->position().z());
3306 // if(fabs(v->position().z()-v1->position().z())>1){
3307 // cout << message << " zdiffrec=" << v->position().z()-v1->position().z() << " " << v->ndof() << " " << v1->ndof() << endl;
3308 // dumpThisEvent_=true;
3309 // }
3310 
3311  double z0=v->position().z()-vertexBeamSpot_.z0();
3312  double z1=v1->position().z()-vertexBeamSpot_.z0();
3313  Fill(h,"zPUcand",z0); Fill(h,"zPUcand",z1);
3314  Fill(h,"ndofPUcand",v->ndof()); Fill(h,"ndofPUcand",v1->ndof());
3315 
3316  Fill(h,"zdiffvsz",z1-z0,0.5*(z1+z0));
3317 
3318  if ((v->ndof()>2) && (v1->ndof()>2)){
3319  Fill(h,"zdiffrec2",v->position().z()-v1->position().z());
3320  Fill(h,"zPUcand2",z0);
3321  Fill(h,"zPUcand2",z1);
3322  Fill(h,"ndofPUcand2",v->ndof());
3323  Fill(h,"ndofPUcand2",v1->ndof());
3324  Fill(h,"zvszrec2",z0, z1);
3325  Fill(h,"pzvspz2",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
3326  }
3327 
3328  if ((v->ndof()>4) && (v1->ndof()>4)){
3329  Fill(h,"zdiffvsz4",z1-z0,0.5*(z1+z0));
3330  Fill(h,"zdiffrec4",v->position().z()-v1->position().z());
3331  Fill(h,"zvszrec4",z0, z1);
3332  Fill(h,"pzvspz4",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
3333  //cout << "ndof4 pu-candidate " << run_ << " " << event_ << endl ;
3334  if(fabs(z0-z1)>1.0){
3335  Fill(h,"xbeamPUcand",v->position().x()-vertexBeamSpot_.x0());
3336  Fill(h,"ybeamPUcand",v->position().y()-vertexBeamSpot_.y0());
3337  Fill(h,"xbeamPullPUcand",(v->position().x()-vertexBeamSpot_.x0())/v->xError());
3338  Fill(h,"ybeamPullPUcand",(v->position().y()-vertexBeamSpot_.y0())/v->yError());
3339  //Fill(h,"sumwOverNtkPUcand",sumw/v->tracksSize());
3340  //Fill(h,"sumwOverNtkPUcand",sumw/v1->tracksSize());
3341  Fill(h,"ndofOverNtkPUcand",v->ndof()/v->tracksSize());
3342  Fill(h,"ndofOverNtkPUcand",v1->ndof()/v1->tracksSize());
3343  Fill(h,"xbeamPUcand",v1->position().x()-vertexBeamSpot_.x0());
3344  Fill(h,"ybeamPUcand",v1->position().y()-vertexBeamSpot_.y0());
3345  Fill(h,"xbeamPullPUcand",(v1->position().x()-vertexBeamSpot_.x0())/v1->xError());
3346  Fill(h,"ybeamPullPUcand",(v1->position().y()-vertexBeamSpot_.y0())/v1->yError());
3347  Fill(h,"zPUcand4",z0);
3348  Fill(h,"zPUcand4",z1);
3349  Fill(h,"ndofPUcand4",v->ndof());
3350  Fill(h,"ndofPUcand4",v1->ndof());
3351  for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){ if(v->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v)); }
3352  for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v1)); }
3353  }
3354  }
3355 
3356  if ((v->ndof()>4) && (v1->ndof()>2) && (v1->ndof()<4)){
3357  for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUfake",**t, &(*v1)); }
3358  }
3359 
3360  if ((v->ndof()>8) && (v1->ndof()>8)){
3361  Fill(h,"zdiffrec8",v->position().z()-v1->position().z());
3362  if(dumpPUcandidates_ && fabs(z0-z1)>1.0){
3363  cout << "PU candidate " << run_ << " " << event_ << " " << message << " zdiff=" <<z0-z1 << endl;
3364 // cerr << "PU candidate " << run_ << " "<< event_ << " " << message << endl;
3365  dumpThisEvent_=true;
3366  }
3367  }
3368 
3369  }
3370 
3371  // test track links, use reconstructed vertices
3372  bool problem = false;
3373  Fill(h,"nans",1.,std::isnan(v->position().x())*1.);
3374  Fill(h,"nans",2.,std::isnan(v->position().y())*1.);
3375  Fill(h,"nans",3.,std::isnan(v->position().z())*1.);
3376 
3377  int index = 3;
3378  for (int i = 0; i != 3; i++) {
3379  for (int j = i; j != 3; j++) {
3380  index++;
3381  Fill(h,"nans",index*1., std::isnan(v->covariance(i, j))*1.);
3382  if (std::isnan(v->covariance(i, j))) problem = true;
3383  // in addition, diagonal element must be positive
3384  if (j == i && v->covariance(i, j) < 0) {
3385  Fill(h,"nans",index*1., 1.);
3386  problem = true;
3387  }
3388  }
3389  }
3390 
3391  try {
3392  for(trackit_t t = v->tracks_begin();
3393  t!=v->tracks_end(); t++) {
3394  // illegal charge
3395  if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
3396  Fill(h,"tklinks",0.);
3397  }
3398  else {
3399  Fill(h,"tklinks",1.);
3400  }
3401  }
3402  } catch (...) {
3403  // exception thrown when trying to use linked track
3404  Fill(h,"tklinks",0.);
3405  }
3406 
3407  if (problem) {
3408  // analyze track parameter covariance definiteness
3409  double data[25];
3410  try {
3411  int itk = 0;
3412  for(trackit_t t = v->tracks_begin();
3413  t!=v->tracks_end(); t++) {
3414  std::cout << "Track " << itk++ << std::endl;
3415  int i2 = 0;
3416  for (int i = 0; i != 5; i++) {
3417  for (int j = 0; j != 5; j++) {
3418  data[i2] = (**t).covariance(i, j);
3419  std::cout << std:: scientific << data[i2] << " ";
3420  i2++;
3421  }
3422  std::cout << std::endl;
3423  }
3424  gsl_matrix_view m
3425  = gsl_matrix_view_array (data, 5, 5);
3426 
3427  gsl_vector *eval = gsl_vector_alloc (5);
3428  gsl_matrix *evec = gsl_matrix_alloc (5, 5);
3429 
3430  gsl_eigen_symmv_workspace * w =
3431  gsl_eigen_symmv_alloc (5);
3432 
3433  gsl_eigen_symmv (&m.matrix, eval, evec, w);
3434 
3435  gsl_eigen_symmv_free (w);
3436 
3437  gsl_eigen_symmv_sort (eval, evec,
3438  GSL_EIGEN_SORT_ABS_ASC);
3439 
3440  // print sorted eigenvalues
3441  {
3442  int i;
3443  for (i = 0; i < 5; i++) {
3444  double eval_i
3445  = gsl_vector_get (eval, i);
3446  //gsl_vector_view evec_i
3447  // = gsl_matrix_column (evec, i);
3448 
3449  printf ("eigenvalue = %g\n", eval_i);
3450  // printf ("eigenvector = \n");
3451  // gsl_vector_fprintf (stdout,
3452  // &evec_i.vector, "%g");
3453  }
3454  }
3455  }
3456  }
3457  catch (...) {
3458  // exception thrown when trying to use linked track
3459  break;
3460  }// catch()
3461  }// if (problem)
3462 
3463 
3464 
3465  } // vertex loop (v)
3466 
3467 
3468  // 2nd highest ndof
3469  if (ndof2>0){
3470  Fill(h,"ndofnr2",ndof2);
3471  if(fabs(zndof1-zndof2)>1) Fill(h,"ndofnr2d1cm",ndof2);
3472  if(fabs(zndof1-zndof2)>2) Fill(h,"ndofnr2d2cm",ndof2);
3473  }
3474 
3475 
3476 }
3477 
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:42
reco::TrackBase::qoverp
double qoverp() const
q/p
Definition: TrackBase.h:115
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:129
type
type
Definition: HCALResponse.h:22
PrimaryVertexAnalyzer4PU::theTrackFilter
TrackFilterForPVFinding theTrackFilter
Definition: PrimaryVertexAnalyzer4PU.h:359
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:44
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:69
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
reco::TrackBase::d0Error
double d0Error() const
error on d0
Definition: TrackBase.h:211
PrimaryVertexAnalyzer4PU::matchVertex
bool matchVertex(const simPrimaryVertex &vsim, const reco::Vertex &vrec)
Definition: PrimaryVertexAnalyzer4PU.cc:974
SiPixelRawToDigiRegional_cfi.beamSpot
tuple beamSpot
Definition: SiPixelRawToDigiRegional_cfi.py:11
PrimaryVertexAnalyzer4PU::getSimTrkParameters
std::vector< SimPart > getSimTrkParameters(edm::Handle< edm::SimTrackContainer > &simTrks, edm::Handle< edm::SimVertexContainer > &simVtcs, double simUnit=1.0)
Definition: PrimaryVertexAnalyzer4PU.cc:733
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1093
EncodedEventId::event
int event() const
get the contents of the subdetector field (should be protected?)
Definition: EncodedEventId.h:26
parents
TPRegexp parents
Definition: eve_filter.cc:24
reco::Track::lost
unsigned short lost() const
Number of lost (=invalid) hits on track.
Definition: Track.h:102
PrimaryVertexAnalyzer4PU::getSimEvents
std::vector< PrimaryVertexAnalyzer4PU::SimEvent > getSimEvents(edm::Handle< TrackingParticleCollection >, edm::Handle< TrackingVertexCollection >, edm::Handle< edm::View< reco::Track > >)
Definition: PrimaryVertexAnalyzer4PU.cc:1578
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:15
PrimaryVertexAnalyzer4PU::fillTrackHistos
void fillTrackHistos(std::map< std::string, TH1 * > &h, const std::string &ttype, const reco::Track &t, const reco::Vertex *v=NULL)
Definition: PrimaryVertexAnalyzer4PU.cc:1003
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:116
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:45
reco::TransientTrack::setBeamSpot
void setBeamSpot(const reco::BeamSpot &beamSpot)
Definition: TransientTrack.h:73
reco::TrackBase::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:111
PrimaryVertexAnalyzer4PU::ParameterVector
reco::TrackBase::ParameterVector ParameterVector
Definition: PrimaryVertexAnalyzer4PU.h:79
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
PrimaryVertexAnalyzer4PU::hsimPV
std::map< std::string, TH1 * > hsimPV
Definition: PrimaryVertexAnalyzer4PU.h:353
TrackingVertex::g4v_iterator
std::vector< SimVertex >::const_iterator g4v_iterator
Definition: TrackingVertex.h:34
PrimaryVertexAnalyzer4PU::recoBeamSpotHandle_
edm::Handle< reco::BeamSpot > recoBeamSpotHandle_
Definition: PrimaryVertexAnalyzer4PU.h:362
PrimaryVertexAnalyzer4PU::SimEvent::tk
std::vector< reco::TransientTrack > tk
Definition: PrimaryVertexAnalyzer4PU.h:148
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:117
reco::TrackBase::dxyError
double dxyError() const
error on dxy
Definition: TrackBase.h:209
funct::D0
Divides< arg, void > D0
Definition: Factorize.h:143
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:165
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
reco::HitPattern::pixelLayersWithMeasurement
int pixelLayersWithMeasurement() const
Definition: HitPattern.h:710
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:62
Measurement1D::error
double error() const
Definition: Measurement1D.h:30
abs
#define abs(x)
Definition: mlp_lapack.h:159
edm::EventBase::luminosityBlock
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
PrimaryVertexAnalyzer4PU::match
static bool match(const ParameterVector &a, const ParameterVector &b)
Definition: PrimaryVertexAnalyzer4PU.cc:962
PrimaryVertexAnalyzer4PU::theB_
edm::ESHandle< TransientTrackBuilder > theB_
Definition: PrimaryVertexAnalyzer4PU.h:363
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:139
NULL
#define NULL
Definition: scimark2.h:8
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
edm::Event::getByType
bool getByType(Handle< PROD > &result) const
Definition: Event.h:398
reco::TrackBase::ParameterVector
math::Vector< dimension >::type ParameterVector
parameter vector
Definition: TrackBase.h:70
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:133
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
PrimaryVertexAnalyzer4PU::printRecVtxs
void printRecVtxs(const edm::Handle< reco::VertexCollection > recVtxs, std::string title="Reconstructed Vertices")
Definition: PrimaryVertexAnalyzer4PU.cc:1124
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:93
reco::TransientTrack::stateAtBeamLine
TrajectoryStateClosestToBeamLine stateAtBeamLine() const
Definition: TransientTrack.h:118
detailsBasic3DVector::z
double double double z
Definition: newBasic3DVector.h:17
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:67
PrimaryVertexAnalyzer4PU::hnoBS
std::map< std::string, TH1 * > hnoBS
Definition: PrimaryVertexAnalyzer4PU.h:349
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
reco::TrackBase::algo
TrackAlgorithm algo() const
Definition: TrackBase.h:332
PrimaryVertexAnalyzer4PU::Fill
void Fill(std::map< std::string, TH1 * > &h, std::string s, double x)
Definition: PrimaryVertexAnalyzer4PU.h:186
PrimaryVertexAnalyzer4PU::hDA
std::map< std::string, TH1 * > hDA
Definition: PrimaryVertexAnalyzer4PU.h:350
PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP
void analyzeVertexCollectionTP(std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< SimEvent > &simEvt, const std::string message="")
Definition: PrimaryVertexAnalyzer4PU.cc:2447
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
errorMatrix2Lands_multiChannel.id
tuple id
Definition: errorMatrix2Lands_multiChannel.py:151
PrimaryVertexAnalyzer4PU::isFinalstateParticle
bool isFinalstateParticle(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:985
PrimaryVertexAnalyzer4PU::isCharged
bool isCharged(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:990
iEvent
int iEvent
Definition: GenABIO.cc:243
reco::TrackBase::trackerExpectedHitsOuter
const HitPattern & trackerExpectedHitsOuter() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers after the last...
Definition: TrackBase.h:227
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:141
PrimaryVertexAnalyzer4PU::truthMatchedTrack
bool truthMatchedTrack(edm::RefToBase< reco::Track >, TrackingParticleRef &)
Definition: PrimaryVertexAnalyzer4PU.cc:1517
edm::detail::isnan
bool isnan(float x)
Definition: math.h:13
reco::Vertex::trackWeight
float trackWeight(const TrackBaseRef &r) const
returns the weight with which a Track has contributed to the vertex-fit.
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.h:705
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
EncodedEventId::bunchCrossing
int bunchCrossing() const
get the detector field from this detid
Definition: EncodedEventId.h:24
TransientVertex::position
GlobalPoint position() const
Definition: TransientVertex.h:113
PrimaryVertexAnalyzer4PU::bookVertexHistograms
std::map< std::string, TH1 * > bookVertexHistograms()
Definition: PrimaryVertexAnalyzer4PU.cc:119
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:131
PrimaryVertexAnalyzer4PU::printPVTrks
void printPVTrks(const edm::Handle< reco::TrackCollection > &recTrks, const edm::Handle< reco::VertexCollection > &recVtxs, std::vector< SimPart > &tsim, std::vector< SimEvent > &simEvt, const bool selectedOnly=true)
Definition: PrimaryVertexAnalyzer4PU.cc:1283
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::TrackBase::qualityMask
int qualityMask() const
Definition: TrackBase.h:286
AdaptiveVertexFitter::vertex
virtual CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &) const
Definition: AdaptiveVertexFitter.cc:161
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
edm::View::size_type
unsigned int size_type
Definition: View.h:85
reco::HitPattern::numberOfHits
int numberOfHits() const
Definition: HitPattern.cc:213
j
int j
Definition: DBlmapReader.cc:9
reco::Vertex::z
double z() const
y coordinate
Definition: Vertex.h:99
reco::TrackBase::lambda
double lambda() const
Lambda angle.
Definition: TrackBase.h:119
PrimaryVertexAnalyzer4PU::getSimPVs
std::vector< simPrimaryVertex > getSimPVs(const edm::Handle< edm::HepMCProduct > evtMC)
Definition: PrimaryVertexAnalyzer4PU.cc:1710
f
double f[11][100]
Definition: MuScleFitUtils.cc:80
ChiSquaredProbability
float ChiSquaredProbability(double chiSquared, double nrDOF)
Definition: ChiSquaredProbability.cc:13
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:87
reco::TrackBase::trackerExpectedHitsInner
const HitPattern & trackerExpectedHitsInner() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers before the fir...
Definition: TrackBase.h:225
PrimaryVertexAnalyzer4PU::printSimVtxs
void printSimVtxs(const edm::Handle< edm::SimVertexContainer > simVtxs)
Definition: PrimaryVertexAnalyzer4PU.cc:1153
RPCpg::mu
const int mu
Definition: Constants.h:23
end
#define end
Definition: vmac.h:38
edm::EventBase::orbitNumber
int orbitNumber() const
Definition: EventBase.h:63
reco::Track::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:63
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:223
PrimaryVertexAnalyzer4PU::printEventSummary
void printEventSummary(std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< SimEvent > &simEvt, const std::string message)
Definition: PrimaryVertexAnalyzer4PU.cc:2222
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:76
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:10
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:89
TrackAssociatorBase::associateRecoToSim
virtual reco::RecoToSimCollection associateRecoToSim(edm::Handle< edm::View< reco::Track > > &tCH, edm::Handle< TrackingParticleCollection > &tPCH, const edm::Event *event=0, const edm::EventSetup *setup=0) const
compare reco to sim the handle of reco::Track and TrackingParticle collections
Definition: TrackAssociatorBase.h:44
nt
int nt
Definition: AMPTWrapper.h:32
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:137
PrimaryVertexAnalyzer4PU::hBS
std::map< std::string, TH1 * > hBS
Definition: PrimaryVertexAnalyzer4PU.h:348
PrimaryVertexAnalyzer4PU::SimEvent::tkprim
std::vector< reco::TransientTrack > tkprim
Definition: PrimaryVertexAnalyzer4PU.h:149
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:215
trackit_t
reco::Vertex::trackRef_iterator trackit_t
Definition: PrimaryVertexAnalyzer4PU.cc:66
edm::getReleaseVersion
std::string getReleaseVersion()
Definition: GetReleaseVersion.cc:7
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:147
PrimaryVertexAnalyzer4PU::associatorByHits_
TrackAssociatorBase * associatorByHits_
Definition: PrimaryVertexAnalyzer4PU.h:355
DetId
Definition: DetId.h:20
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:75
PrimaryVertexAnalyzer4PU::matchRecTracksToVertex
void matchRecTracksToVertex(simPrimaryVertex &pv, const std::vector< SimPart > &tsim, const edm::Handle< reco::TrackCollection > &recTrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1418
Measurement1D::significance
double significance() const
Definition: Measurement1D.h:32
M_PI
#define M_PI
Definition: BFit3D.cc:3
PrimaryVertexAnalyzer4PU::isResonance
bool isResonance(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:979
reco::TransientTrack::track
const Track & track() const
Definition: TransientTrack.h:114
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
PrimaryVertexAnalyzer4PU::printRecTrks
void printRecTrks(const edm::Handle< reco::TrackCollection > &recTrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1197
PrimaryVertexAnalyzer4PU::getTc
void getTc(const std::vector< reco::TransientTrack > &, double &, double &, double &, double &, double &)
Definition: PrimaryVertexAnalyzer4PU.cc:1466
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
part
part
Definition: HCALResponse.h:21
TrackingVertexCollection
std::vector< TrackingVertex > TrackingVertexCollection
Definition: TrackingVertexContainer.h:8
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
PrimaryVertexAnalyzer4PU::r2s_
reco::RecoToSimCollection r2s_
Definition: PrimaryVertexAnalyzer4PU.h:356
reco::HitPattern::pixelBarrelLayersWithMeasurement
int pixelBarrelLayersWithMeasurement() const
Definition: HitPattern.cc:337
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:62
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:89
edm::View::size
size_type size() const
b
double b
Definition: hdecay.h:120
Measurement1D::value
double value() const
Definition: Measurement1D.h:28
PrimaryVertexAnalyzer4PU::supf
int * supf(std::vector< SimPart > &simtrks, const reco::TrackCollection &trks)
Definition: PrimaryVertexAnalyzer4PU.cc:831
edm::Handle::product
T const * product() const
Definition: Handle.h:74
PrimaryVertexAnalyzer4PU::PrimaryVertexAnalyzer4PU
PrimaryVertexAnalyzer4PU(const edm::ParameterSet &)
Definition: PrimaryVertexAnalyzer4PU.cc:74
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
TrackingVertex::eventId
const EncodedEventId & eventId() const
Definition: TrackingVertex.h:95
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:23
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:145
PrimaryVertexAnalyzer4PU::printSimTrks
void printSimTrks(const edm::Handle< edm::SimTrackContainer > simVtrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1176
begin
#define begin
Definition: vmac.h:31
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
PrimaryVertexAnalyzer4PU::add
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
Definition: PrimaryVertexAnalyzer4PU.h:184
a
double a
Definition: hdecay.h:121
lt
< trclass="colgroup">< tdclass="colgroup"colspan=5 > DT local reconstruction</td ></tr >< tr >< td >< ahref="classDTRecHit1DPair.html"> DTRecHit1DPair</a ></td >< td >< ahref="DataFormats_DTRecHit.html"> edm::RangeMap & lt
Definition: RecoLocalMuon_RECO.doi:7
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
reco::Track::found
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:100
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:38
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:67
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:113
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:63
reco::HitPattern::printHitPattern
void printHitPattern(int position, std::ostream &stream) const
Definition: HitPattern.cc:556
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:40
benchmark_cfg.pdgId
tuple pdgId
Definition: benchmark_cfg.py:18
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
vdt::x
x
Definition: VDTMath.h:216
PrimaryVertexAnalyzer4PU::dumpHitInfo
void dumpHitInfo(const reco::Track &t)
Definition: PrimaryVertexAnalyzer4PU.cc:1102
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
reco::TrackBase::dxy
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:121
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
PrimaryVertexAnalyzer4PU::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexAnalyzer4PU.cc:1912
TransientVertex::isValid
bool isValid() const
Definition: TransientVertex.h:133
PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks
std::vector< edm::RefToBase< reco::Track > > getTruthMatchedVertexTracks(const reco::Vertex &)
Definition: PrimaryVertexAnalyzer4PU.cc:1551
edm::RefToBase::get
value_type const * get() const
Definition: RefToBase.h:212
PrimaryVertexAnalyzer4PU::pdt_
edm::ESHandle< ParticleDataTable > pdt_
Definition: PrimaryVertexAnalyzer4PU.h:329
PrimaryVertexAnalyzer4PU::z2tp_
std::map< double, TrackingParticleRef > z2tp_
Definition: PrimaryVertexAnalyzer4PU.h:357
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
v
mathSSE::Vec4< T > v
Definition: newBasic3DVector.h:354
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
reco::Vertex::tracksSize
size_t tracksSize() const
number of tracks
Definition: Vertex.cc:35
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:135
PrimaryVertexAnalyzer4PU::analyzeVertexCollection
void analyzeVertexCollection(std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< simPrimaryVertex > &simpv, const std::string message="")
Definition: PrimaryVertexAnalyzer4PU.cc:2758
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:143
reco::TrackBase::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:72
w
T w() const
Definition: newBasic3DVector.h:234
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:22
TrackingVertex::position
const LorentzVector & position() const
Definition: TrackingVertex.h:91
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:65
reco::Track::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:65
Generated for CMSSW Reference Manual by doxygen 1.8.5