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