CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PrimaryVertexAnalyzer4PU.cc
Go to the documentation of this file.
1 #include "Validation/RecoVertex/interface/PrimaryVertexAnalyzer4PU.h"
2 
3 #include "FWCore/ParameterSet/interface/ParameterSet.h"
4 #include "FWCore/MessageLogger/interface/MessageLogger.h"
5 #include "FWCore/Version/interface/GetReleaseVersion.h"
6 #include "MagneticField/Engine/interface/MagneticField.h"
7 
8 // reco track and vertex
9 #include "DataFormats/TrackReco/interface/Track.h"
10 #include "DataFormats/VertexReco/interface/Vertex.h"
11 #include "DataFormats/Math/interface/Point3D.h"
12 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
13 
14 // simulated vertices,..., add <use name=SimDataFormats/Vertex> and <../Track>
15 #include <SimDataFormats/Vertex/interface/SimVertex.h>
16 #include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
17 #include <SimDataFormats/Track/interface/SimTrack.h>
18 #include <SimDataFormats/Track/interface/SimTrackContainer.h>
19 
20 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
21 #include "SimDataFormats/CrossingFrame/interface/CrossingFramePlaybackInfo.h"
22 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
23 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
24 
25 //generator level + CLHEP
26 #include "HepMC/GenEvent.h"
27 #include "HepMC/GenVertex.h"
28 
29 
30 // TrackingParticle
31 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
32 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
33 //associator
34 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
35 
36 
37 // fit
38 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
39 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
40 
41 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
42 
43 // Root
44 #include <TH1.h>
45 #include <TH2.h>
46 #include <TFile.h>
47 #include <TProfile.h>
48 
49 #include <cmath>
50 #include <gsl/gsl_math.h>
51 #include <gsl/gsl_eigen.h>
52 
53 
54 // cluster stufff
55 //#include "DataFormats/TrackRecoTrack.h"
56 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
57 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
58 
59 
60 using namespace edm;
61 using namespace reco;
62 using namespace std;
63 //
64 // constants, enums and typedefs
65 //
66 typedef reco::Vertex::trackRef_iterator trackit_t;
67 //
68 // static data member definitions
69 //
70 
71 //
72 // constructors and destructor
73 //
74 PrimaryVertexAnalyzer4PU::PrimaryVertexAnalyzer4PU(const ParameterSet& iConfig):theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"))
75 {
76  //now do what ever initialization is needed
77  simG4_=iConfig.getParameter<edm::InputTag>( "simG4" );
78  recoTrackProducer_= iConfig.getUntrackedParameter<std::string>("recoTrackProducer");
79  // open output file to store histograms}
80  outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
81 
82  rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE");
83  verbose_= iConfig.getUntrackedParameter<bool>("verbose", false);
84  doMatching_= iConfig.getUntrackedParameter<bool>("matching", false);
85  printXBS_= iConfig.getUntrackedParameter<bool>("XBS", false);
86  simUnit_= 1.0; // starting with CMSSW_1_2_x ??
87  if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
88  simUnit_=0.1; // for use in CMSSW_1_1_1 tutorial
89  }
90 
91  dumpPUcandidates_=iConfig.getUntrackedParameter<bool>("dumpPUcandidates", false);
92 
93  zmatch_=iConfig.getUntrackedParameter<double>("zmatch", 0.0500);
94  cout << "PrimaryVertexAnalyzer4PU: zmatch=" << zmatch_ << endl;
95  eventcounter_=0;
96  dumpcounter_=0;
97  ndump_=10;
98  DEBUG_=false;
99  //DEBUG_=true;
100 }
101 
102 
103 
104 
105 PrimaryVertexAnalyzer4PU::~PrimaryVertexAnalyzer4PU()
106 {
107  // do anything here that needs to be done at desctruction time
108  // (e.g. close files, deallocate resources etc.)
109  delete rootFile_;
110 }
111 
112 
113 
114 //
115 // member functions
116 //
117 
118 
119 std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::bookVertexHistograms(){
120  std::map<std::string, TH1*> h;
121 
122  // release validation histograms used in DoCompare.C
123  h["nbtksinvtx"] = new TH1F("nbtksinvtx","reconstructed tracks in vertex",40,-0.5,39.5);
124  h["nbtksinvtxPU"] = new TH1F("nbtksinvtxPU","reconstructed tracks in vertex",40,-0.5,39.5);
125  h["nbtksinvtxTag"]= new TH1F("nbtksinvtxTag","reconstructed tracks in vertex",40,-0.5,39.5);
126  h["resx"] = new TH1F("resx","residual x",100,-0.04,0.04);
127  h["resy"] = new TH1F("resy","residual y",100,-0.04,0.04);
128  h["resz"] = new TH1F("resz","residual z",100,-0.1,0.1);
129  h["resz10"] = new TH1F("resz10","residual z",100,-1.0,1.);
130  h["pullx"] = new TH1F("pullx","pull x",100,-25.,25.);
131  h["pully"] = new TH1F("pully","pull y",100,-25.,25.);
132  h["pullz"] = new TH1F("pullz","pull z",100,-25.,25.);
133  h["vtxchi2"] = new TH1F("vtxchi2","chi squared",100,0.,100.);
134  h["vtxndf"] = new TH1F("vtxndf","degrees of freedom",500,0.,100.);
135  h["vtxndfc"] = new TH1F("vtxndfc","expected 2nd highest ndof",500,0.,100.);
136 
137  h["vtxndfvsntk"] = new TH2F("vtxndfvsntk","ndof vs #tracks",20,0.,100, 20, 0., 200.);
138  h["vtxndfoverntk"]= new TH1F("vtxndfoverntk","ndof / #tracks",40,0.,2.);
139  h["vtxndf2overntk"]= new TH1F("vtxndf2overntk","(ndof+2) / #tracks",40,0.,2.);
140  h["tklinks"] = new TH1F("tklinks","Usable track links",2,-0.5,1.5);
141  h["nans"] = new TH1F("nans","Illegal values for x,y,z,xx,xy,xz,yy,yz,zz",9,0.5,9.5);
142 
143 
144  // raw
145  add(h, new TH1F("szRecVtx","size of recvtx collection",20, -0.5, 19.5));
146  add(h, new TH1F("isFake","fake vertex",2, -0.5, 1.5));
147  add(h, new TH1F("isFake1","fake vertex or ndof<0",2, -0.5, 1.5));
148  add(h, new TH1F("bunchCrossing","bunchCrossing",4000, 0., 4000.));
149  add(h, new TH2F("bunchCrossingLogNtk","bunchCrossingLogNtk",4000, 0., 4000.,5,0., 5.));
150  add(h, new TH1F("highpurityTrackFraction","fraction of high purity tracks",20, 0., 1.));
151  add(h, new TH2F("trkchi2vsndof","vertices chi2 vs ndof",50, 0., 100., 50, 0., 200.));
152  add(h, new TH1F("trkchi2overndof","vertices chi2 / ndof",50, 0., 5.));
153  // two track vertices
154  add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
155  add(h,new TH1F("2trkmassSS","two-track vertices mass (same sign)",100, 0., 2.));
156  add(h,new TH1F("2trkmassOS","two-track vertices mass (opposite sign)",100, 0., 2.));
157  add(h,new TH1F("2trkdphi","two-track vertices delta-phi",360, 0, 2*M_PI));
158  add(h,new TH1F("2trkseta","two-track vertices sum-eta",50, -2., 2.));
159  add(h,new TH1F("2trkdphicurl","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
160  add(h,new TH1F("2trksetacurl","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
161  add(h,new TH1F("2trkdetaOS","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
162  add(h,new TH1F("2trkdetaSS","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
163  // two track PU vertices
164  add(h,new TH1F("2trkmassSSPU","two-track vertices mass (same sign)",100, 0., 2.));
165  add(h,new TH1F("2trkmassOSPU","two-track vertices mass (opposite sign)",100, 0., 2.));
166  add(h,new TH1F("2trkdphiPU","two-track vertices delta-phi",360, 0, 2*M_PI));
167  add(h,new TH1F("2trksetaPU","two-track vertices sum-eta",50, -2., 2.));
168  add(h,new TH1F("2trkdphicurlPU","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
169  add(h,new TH1F("2trksetacurlPU","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
170  add(h,new TH1F("2trkdetaOSPU","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
171  add(h,new TH1F("2trkdetaSSPU","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
172  // three track vertices
173  add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
174  add(h,new TH2F("3trkchi2vsndof","three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
175  add(h,new TH2F("4trkchi2vsndof","four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
176  add(h,new TH2F("5trkchi2vsndof","five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
177  // same for fakes
178  add(h,new TH2F("fake2trkchi2vsndof","fake two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
179  add(h,new TH2F("fake3trkchi2vsndof","fake three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
180  add(h,new TH2F("fake4trkchi2vsndof","fake four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
181  add(h,new TH2F("fake5trkchi2vsndof","fake five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
182  h["resxr"] = new TH1F("resxr","relative residual x",100,-0.04,0.04);
183  h["resyr"] = new TH1F("resyr","relative residual y",100,-0.04,0.04);
184  h["reszr"] = new TH1F("reszr","relative residual z",100,-0.1,0.1);
185  h["pullxr"] = new TH1F("pullxr","relative pull x",100,-25.,25.);
186  h["pullyr"] = new TH1F("pullyr","relative pull y",100,-25.,25.);
187  h["pullzr"] = new TH1F("pullzr","relative pull z",100,-25.,25.);
188  h["vtxprob"] = new TH1F("vtxprob","chisquared probability",100,0.,1.);
189  h["eff"] = new TH1F("eff","efficiency",2, -0.5, 1.5);
190  h["efftag"] = new TH1F("efftag","efficiency tagged vertex",2, -0.5, 1.5);
191  h["zdistancetag"] = new TH1F("zdistancetag","z-distance between tagged and generated",100, -0.1, 0.1);
192  h["abszdistancetag"] = new TH1F("abszdistancetag","z-distance between tagged and generated",1000, 0., 1.0);
193  h["abszdistancetagcum"] = new TH1F("abszdistancetagcum","z-distance between tagged and generated",1000, 0., 1.0);
194  h["puritytag"] = new TH1F("puritytag","purity of primary vertex tags",2, -0.5, 1.5);
195  h["effvsptsq"] = new TProfile("effvsptsq","efficiency vs ptsq",20, 0., 10000., 0, 1.);
196  h["effvsnsimtrk"] = new TProfile("effvsnsimtrk","efficiency vs # simtracks",50, 0., 50., 0, 1.);
197  h["effvsnrectrk"] = new TProfile("effvsnrectrk","efficiency vs # rectracks",50, 0., 50., 0, 1.);
198  h["effvsnseltrk"] = new TProfile("effvsnseltrk","efficiency vs # selected tracks",50, 0., 50., 0, 1.);
199  h["effvsz"] = new TProfile("effvsz","efficiency vs z",20, -20., 20., 0, 1.);
200  h["effvsz2"] = new TProfile("effvsz2","efficiency vs z (2mm)",20, -20., 20., 0, 1.);
201  h["effvsr"] = new TProfile("effvsr","efficiency vs r",20, 0., 1., 0, 1.);
202  h["xresvsntrk"] = new TProfile("xresvsntrk","xresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
203  h["yresvsntrk"] = new TProfile("yresvsntrk","yresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
204  h["zresvsntrk"] = new TProfile("zresvsntrk","zresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
205  h["cpuvsntrk"] = new TProfile("cpuvsntrk","cpu time vs # of fitted tracks",40, 0., 200., 0, 200.);
206  h["cpucluvsntrk"] = new TProfile("cpucluvsntrk","clustering cpu time # of tracks",40, 0., 200., 0, 10.);
207  h["cpufit"] = new TH1F("cpufit","cpu time for fitting",100, 0., 200.);
208  h["cpuclu"] = new TH1F("cpuclu","cpu time for clustering",100, 0., 200.);
209  h["nbtksinvtx2"] = new TH1F("nbtksinvtx2","reconstructed tracks in vertex",40,0.,200.);
210  h["nbtksinvtxPU2"] = new TH1F("nbtksinvtxPU2","reconstructed tracks in vertex",40,0.,200.);
211  h["nbtksinvtxTag2"] = new TH1F("nbtksinvtxTag2","reconstructed tracks in vertex",40,0.,200.);
212 
213  h["xrec"] = new TH1F("xrec","reconstructed x",100,-0.1,0.1);
214  h["yrec"] = new TH1F("yrec","reconstructed y",100,-0.1,0.1);
215  h["zrec"] = new TH1F("zrec","reconstructed z",100,-20.,20.);
216  h["err1"] = new TH1F("err1","error 1",100,0.,0.1);
217  h["err2"] = new TH1F("err2","error 2",100,0.,0.1);
218  h["errx"] = new TH1F("errx","error x",100,0.,0.1);
219  h["erry"] = new TH1F("erry","error y",100,0.,0.1);
220  h["errz"] = new TH1F("errz","error z",100,0.,2.0);
221  h["errz1"] = new TH1F("errz1","error z",100,0.,0.2);
222 
223  h["zrecNt100"] = new TH1F("zrecNt100","reconstructed z for high multiplicity events",80,-40.,40.);
224  add(h, new TH2F("zrecvsnt","reconstructed z vs number of tracks",100,-50., 50., 20, 0., 200.));
225  add(h, new TH2F("xyrec","reconstructed xy",100, -4., 4., 100, -4., 4.));
226  h["xrecBeam"] = new TH1F("xrecBeam","reconstructed x - beam x",100,-0.1,0.1);
227  h["yrecBeam"] = new TH1F("yrecBeam","reconstructed y - beam y",100,-0.1,0.1);
228  h["zrecBeam"] = new TH1F("zrecBeam","reconstructed z - beam z",100,-20.,20.);
229  h["xrecBeamvsz"] = new TH2F("xrecBeamvsz","reconstructed x - beam x vs z", 20, -20., 20.,100,-0.1,0.1);
230  h["yrecBeamvsz"] = new TH2F("yrecBeamvsz","reconstructed y - beam y vs z", 20, -20., 20.,100,-0.1,0.1);
231  h["xrecBeamvszprof"] = new TProfile("xrecBeamvszprof","reconstructed x - beam x vs z-z0", 20, -20., 20.,-0.1,0.1);
232  h["yrecBeamvszprof"] = new TProfile("yrecBeamvszprof","reconstructed y - beam y vs z-z0", 20, -20., 20.,-0.1,0.1);
233  add(h, new TH2F("xrecBeamvsdxXBS","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1)); // just a test
234  add(h, new TH2F("yrecBeamvsdyXBS","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
235  add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1));
236  add(h, new TH2F("yrecBeamvsdy","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
237  add(h, new TH2F("xrecBeamvsdxR2","reconstructed x - beam x vs resolution",20,0., 0.04, 100, -0.1,0.1));
238  add(h, new TH2F("yrecBeamvsdyR2","reconstructed z - beam z vs resolution",20,0., 0.04, 100, -0.1,0.1));
239  // add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",100,-0.1,0.1, 10, 0., 0.04));
240  // add(h, new TH2F("yrecBeamvsdy","reconstructed y - beam y vs resolution",100,-0.1,0.1, 10, 0., 0.04));
241  h["xrecBeamvsdxprof"] = new TProfile("xrecBeamvsdxprof","reconstructed x - beam x vs resolution",10, 0., 0.04,-0.1,0.1 );
242  h["yrecBeamvsdyprof"] = new TProfile("yrecBeamvsdyprof","reconstructed y - beam y vs resolution",10, 0., 0.04,-0.1,0.1 );
243  add(h, new TProfile("xrecBeam2vsdx2prof","reconstructed x - beam x vs resolution",10,0., 0.002, 0., 0.01));
244  add(h, new TProfile("yrecBeam2vsdy2prof","reconstructed y - beam y vs resolution",10,0., 0.002, 0., 0.01));
245  add(h,new TH2F("xrecBeamvsdx2","reconstructed x - beam x vs resolution",10,0., 0.002, 100, -0.01, 0.01));
246  add(h,new TH2F("yrecBeamvsdy2","reconstructed y - beam y vs resolution",10,0., 0.002, 100, -0.01, 0.01));
247  h["xrecb"] = new TH1F("xrecb","reconstructed x - beam x",100,-0.01,0.01);
248  h["yrecb"] = new TH1F("yrecb","reconstructed y - beam y",100,-0.01,0.01);
249  h["zrecb"] = new TH1F("zrecb","reconstructed z - beam z",100,-20.,20.);
250  h["xrec1"] = new TH1F("xrec1","reconstructed x",100,-4,4);
251  h["yrec1"] = new TH1F("yrec1","reconstructed y",100,-4,4); // should match the sim histos
252  h["zrec1"] = new TH1F("zrec1","reconstructed z",100,-80.,80.);
253  h["xrec2"] = new TH1F("xrec2","reconstructed x",100,-1,1);
254  h["yrec2"] = new TH1F("yrec2","reconstructed y",100,-1,1);
255  h["zrec2"] = new TH1F("zrec2","reconstructed z",100,-40.,40.);
256  h["xrec3"] = new TH1F("xrec3","reconstructed x",100,-0.1,0.1);
257  h["yrec3"] = new TH1F("yrec3","reconstructed y",100,-0.1,0.1);
258  h["zrec3"] = new TH1F("zrec3","reconstructed z",100,-20.,20.);
259  add(h, new TH1F("xrecBeamPull","normalized residuals reconstructed x - beam x",100,-20,20));
260  add(h, new TH1F("yrecBeamPull","normalized residuals reconstructed y - beam y",100,-20,20));
261  add(h, new TH1F("zdiffrec","z-distance between vertices",200,-20., 20.));
262  add(h, new TH1F("zdiffrec2","z-distance between ndof>2 vertices",200,-20., 20.));
263  add(h, new TH1F("zdiffrec4","z-distance between ndof>4 vertices",200,-20., 20.));
264  add(h, new TH1F("zdiffrec8","z-distance between ndof>8 vertices",200,-20., 20.));
265  add(h, new TH2F("zvszrec2","z positions of multiple vertices",200,-20., 20., 200,-20., 20.));
266  add(h, new TH2F("pzvspz2","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
267  add(h, new TH2F("zvszrec4","z positions of multiple vertices",100,-20., 20., 100,-20., 20.));
268  add(h, new TH2F("pzvspz4","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
269  add(h, new TH2F("zdiffvsz","z-distance vs z",100,-10.,10.,30,-15.,15.));
270  add(h, new TH2F("zdiffvsz4","z-distance vs z (ndof>4)",100,-10.,10.,30,-15.,15.));
271  add(h, new TProfile("eff0vsntrec","efficiency vs # reconstructed tracks",50, 0., 50., 0, 1.));
272  add(h, new TProfile("eff0vsntsel","efficiency vs # selected tracks",50, 0., 50., 0, 1.));
273  add(h, new TProfile("eff0ndof0vsntsel","efficiency (ndof>0) vs # selected tracks",50, 0., 50., 0, 1.));
274  add(h, new TProfile("eff0ndof8vsntsel","efficiency (ndof>8) vs # selected tracks",50, 0., 50., 0, 1.));
275  add(h, new TProfile("eff0ndof2vsntsel","efficiency (ndof>2) vs # selected tracks",50, 0., 50., 0, 1.));
276  add(h, new TProfile("eff0ndof4vsntsel","efficiency (ndof>4) vs # selected tracks",50, 0., 50., 0, 1.));
277  add(h, new TH1F("xbeamPUcand","x - beam of PU candidates (ndof>4)",100,-1., 1.));
278  add(h, new TH1F("ybeamPUcand","y - beam of PU candidates (ndof>4)",100,-1., 1.));
279  add(h, new TH1F("xbeamPullPUcand","x - beam pull of PU candidates (ndof>4)",100,-20., 20.));
280  add(h, new TH1F("ybeamPullPUcand","y - beam pull of PU candidates (ndof>4)",100,-20., 20.));
281  add(h, new TH1F("ndofOverNtk","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
282  //add(h, new TH1F("sumwOverNtk","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
283  add(h, new TH1F("ndofOverNtkPUcand","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
284  //add(h, new TH1F("sumwOverNtkPUcand","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
285  add(h, new TH1F("zPUcand","z positions of PU candidates (all)",100,-20., 20.));
286  add(h, new TH1F("zPUcand2","z positions of PU candidates (ndof>2)",100,-20., 20.));
287  add(h, new TH1F("zPUcand4","z positions of PU candidates (ndof>4)",100,-20., 20.));
288  add(h, new TH1F("ndofPUcand","ndof of PU candidates (all)",50,0., 100.));
289  add(h, new TH1F("ndofPUcand2","ndof of PU candidates (ndof>2)",50,0., 100.));
290  add(h, new TH1F("ndofPUcand4","ndof of PU candidates (ndof>4)",50,0., 100.));
291  add(h, new TH1F("ndofnr2","second highest ndof",500,0., 100.));
292  add(h, new TH1F("ndofnr2d1cm","second highest ndof (dz>1cm)",500,0., 100.));
293  add(h, new TH1F("ndofnr2d2cm","second highest ndof (dz>2cm)",500,0., 100.));
294 
295  h["nrecvtx"] = new TH1F("nrecvtx","# of reconstructed vertices", 50, -0.5, 49.5);
296  h["nrecvtx8"] = new TH1F("nrecvtx8","# of reconstructed vertices with ndof>8", 50, -0.5, 49.5);
297  h["nrecvtx2"] = new TH1F("nrecvtx2","# of reconstructed vertices with ndof>2", 50, -0.5, 49.5);
298  h["nrecvtx4"] = new TH1F("nrecvtx4","# of reconstructed vertices with ndof>4", 50, -0.5, 49.5);
299  // h["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
300  h["nrectrk"] = new TH1F("nrectrk","# of reconstructed tracks", 100, -0.5, 99.5);
301  add(h, new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5));
302  add(h, new TH1F("nsimtrkSignal","# of simulated tracks (Signal)", 100, -0.5, 99.5));
303  add(h, new TH1F("nsimtrkPU","# of simulated tracks (PU)", 100, -0.5, 99.5));
304  h["nsimtrk"]->StatOverflows(kTRUE);
305  h["nsimtrkPU"]->StatOverflows(kTRUE);
306  h["nsimtrkSignal"]->StatOverflows(kTRUE);
307  h["xrectag"] = new TH1F("xrectag","reconstructed x, signal vtx",100,-0.05,0.05);
308  h["yrectag"] = new TH1F("yrectag","reconstructed y, signal vtx",100,-0.05,0.05);
309  h["zrectag"] = new TH1F("zrectag","reconstructed z, signal vtx",100,-20.,20.);
310  h["nrectrk0vtx"] = new TH1F("nrectrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
311  h["nseltrk0vtx"] = new TH1F("nseltrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
312  h["nrecsimtrk"] = new TH1F("nrecsimtrk","# rec tracks matched to sim tracks in vertex",100,-0.5, 99.5);
313  h["nrecnosimtrk"] = new TH1F("nrecsimtrk","# rec tracks not matched to sim tracks in vertex",100,-0.5, 99.5);
314  h["trackAssEffvsPt"] = new TProfile("trackAssEffvsPt","track association efficiency vs pt",20, 0., 100., 0, 1.);
315 
316  // cluster stuff
317  h["nseltrk"] = new TH1F("nseltrk","# of reconstructed tracks selected for PV", 100, -0.5, 99.5);
318  h["nclutrkall"] = new TH1F("nclutrkall","# of reconstructed tracks in clusters", 100, -0.5, 99.5);
319  h["nclutrkvtx"] = new TH1F("nclutrkvtx","# of reconstructed tracks in clusters of reconstructed vertices", 100, -0.5, 99.5);
320  h["nclu"] = new TH1F("nclu","# of clusters", 100, -0.5, 99.5);
321  h["nclu0vtx"] = new TH1F("nclu0vtx","# of clusters in events with no PV", 100, -0.5, 99.5);
322  h["zlost1"] = new TH1F("zlost1","z of lost vertices (bad z)", 100, -20., 20.);
323  h["zlost2"] = new TH1F("zlost2","z of lost vertices (no matching cluster)", 100, -20., 20.);
324  h["zlost3"] = new TH1F("zlost3","z of lost vertices (vertex too far from beam)", 100, -20., 20.);
325  h["zlost4"] = new TH1F("zlost4","z of lost vertices (invalid vertex)", 100, -20., 20.);
326  h["selstat"] = new TH1F("selstat","selstat", 5, -2.5, 2.5);
327 
328 
329  // properties of fake vertices (MC only)_
330  add(h, new TH1F("fakeVtxZNdofgt05","z of fake vertices with ndof>0.5", 100, -20., 20.));
331  add(h, new TH1F("fakeVtxZNdofgt2","z of fake vertices with ndof>2", 100, -20., 20.));
332  add(h, new TH1F("fakeVtxZ","z of fake vertices", 100, -20., 20.));
333  add(h, new TH1F("fakeVtxNdof","ndof of fake vertices", 500,0., 100.));
334  add(h,new TH1F("fakeVtxNtrk","number of tracks in fake vertex",20,-0.5, 19.5));
335  add(h,new TH1F("matchedVtxNdof","ndof of matched vertices", 500,0., 100.));
336 
337 
338  // histograms of track quality (Data and MC)
339  string types[] = {"all","sel","bachelor","sellost","wgt05","wlt05","M","tagged","untagged","ndof4","PUcand","PUfake"};
340  for(int t=0; t<12; t++){
341  add(h, new TH1F(("rapidity_"+types[t]).c_str(),"rapidity ",100,-3., 3.));
342  h["z0_"+types[t]] = new TH1F(("z0_"+types[t]).c_str(),"z0 ",80,-40., 40.);
343  h["phi_"+types[t]] = new TH1F(("phi_"+types[t]).c_str(),"phi ",80,-3.14159, 3.14159);
344  h["eta_"+types[t]] = new TH1F(("eta_"+types[t]).c_str(),"eta ",80,-4., 4.);
345  h["pt_"+types[t]] = new TH1F(("pt_"+types[t]).c_str(),"pt ",100,0., 20.);
346  h["found_"+types[t]] = new TH1F(("found_"+types[t]).c_str(),"found hits",20, 0., 20.);
347  h["lost_"+types[t]] = new TH1F(("lost_"+types[t]).c_str(),"lost hits",20, 0., 20.);
348  h["nchi2_"+types[t]] = new TH1F(("nchi2_"+types[t]).c_str(),"normalized track chi2",100, 0., 20.);
349  h["rstart_"+types[t]] = new TH1F(("rstart_"+types[t]).c_str(),"start radius",100, 0., 20.);
350  h["tfom_"+types[t]] = new TH1F(("tfom_"+types[t]).c_str(),"track figure of merit",100, 0., 100.);
351  h["expectedInner_"+types[t]] = new TH1F(("expectedInner_"+types[t]).c_str(),"expected inner hits ",10, 0., 10.);
352  h["expectedOuter_"+types[t]] = new TH1F(("expectedOuter_"+types[t]).c_str(),"expected outer hits ",10, 0., 10.);
353  h["logtresxy_"+types[t]] = new TH1F(("logtresxy_"+types[t]).c_str(),"log10(track r-phi resolution/um)",100, 0., 5.);
354  h["logtresz_"+types[t]] = new TH1F(("logtresz_"+types[t]).c_str(),"log10(track z resolution/um)",100, 0., 5.);
355  h["tpullxy_"+types[t]] = new TH1F(("tpullxy_"+types[t]).c_str(),"track r-phi pull",100, -10., 10.);
356  add(h, new TH2F( ("lvseta_"+types[t]).c_str(),"cluster length vs eta",60,-3., 3., 20, 0., 20));
357  add(h, new TH2F( ("lvstanlambda_"+types[t]).c_str(),"cluster length vs tan lambda",60,-6., 6., 20, 0., 20));
358  add(h, new TH1F( ("restrkz_"+types[t]).c_str(),"z-residuals (track vs vertex)", 200, -5., 5.));
359  add(h, new TH2F( ("restrkzvsphi_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
360  add(h, new TH2F( ("restrkzvseta_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.,3.,200, -5., 5.));
361  add(h, new TH2F( ("pulltrkzvsphi_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
362  add(h, new TH2F( ("pulltrkzvseta_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.,3.,100, -5., 5.));
363  add(h, new TH1F( ("pulltrkz_"+types[t]).c_str(),"normalized z-residuals (track vs vertex)", 100, -5., 5.));
364  add(h, new TH1F( ("sigmatrkz0_"+types[t]).c_str(),"z-resolution (excluding beam)", 100, 0., 5.));
365  add(h, new TH1F( ("sigmatrkz_"+types[t]).c_str(),"z-resolution (including beam)", 100,0., 5.));
366  add(h, new TH1F( ("nbarrelhits_"+types[t]).c_str(),"number of pixel barrel hits", 10, 0., 10.));
367  add(h, new TH1F( ("nbarrelLayers_"+types[t]).c_str(),"number of pixel barrel layers", 10, 0., 10.));
368  add(h, new TH1F( ("nPxLayers_"+types[t]).c_str(),"number of pixel layers (barrel+endcap)", 10, 0., 10.));
369  add(h, new TH1F( ("nSiLayers_"+types[t]).c_str(),"number of Tracker layers ", 20, 0., 20.));
370  add(h, new TH1F( ("trackAlgo_"+types[t]).c_str(),"track algorithm ", 30, 0., 30.));
371  add(h, new TH1F( ("trackQuality_"+types[t]).c_str(),"track quality ", 7, -1., 6.));
372  }
373  add(h, new TH1F("trackWt","track weight in vertex",100,0., 1.));
374 
375 
376  h["nrectrk"]->StatOverflows(kTRUE);
377  h["nrectrk"]->StatOverflows(kTRUE);
378  h["nrectrk0vtx"]->StatOverflows(kTRUE);
379  h["nseltrk0vtx"]->StatOverflows(kTRUE);
380  h["nrecsimtrk"]->StatOverflows(kTRUE);
381  h["nrecnosimtrk"]->StatOverflows(kTRUE);
382  h["nseltrk"]->StatOverflows(kTRUE);
383  h["nbtksinvtx"]->StatOverflows(kTRUE);
384  h["nbtksinvtxPU"]->StatOverflows(kTRUE);
385  h["nbtksinvtxTag"]->StatOverflows(kTRUE);
386  h["nbtksinvtx2"]->StatOverflows(kTRUE);
387  h["nbtksinvtxPU2"]->StatOverflows(kTRUE);
388  h["nbtksinvtxTag2"]->StatOverflows(kTRUE);
389 
390  // pile-up and track assignment related histograms (MC)
391  h["npu0"] =new TH1F("npu0","Number of simulated vertices",40,0.,40.);
392  h["npu1"] =new TH1F("npu1","Number of simulated vertices with >0 track",40,0.,40.);
393  h["npu2"] =new TH1F("npu2","Number of simulated vertices with >1 track",40,0.,40.);
394  h["nrecv"] =new TH1F("nrecv","# of reconstructed vertices", 40, 0, 40);
395  add(h,new TH2F("nrecvsnpu","#rec vertices vs number of sim vertices with >0 tracks", 40, 0., 40, 40, 0, 40));
396  add(h,new TH2F("nrecvsnpu2","#rec vertices vs number of sim vertices with >1 tracks", 40, 0., 40, 40, 0, 40));
397  add(h,new TH1F("sumpt","sumpt of simulated tracks",100,0.,100.));
398  add(h,new TH1F("sumptSignal","sumpt of simulated tracks in Signal events",100,0.,200.));
399  add(h,new TH1F("sumptPU","sumpt of simulated tracks in PU events",100,0.,200.));
400  add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
401  add(h,new TH1F("sumpt2","sumpt2 of simulated tracks",100,0.,100.));
402  add(h,new TH1F("sumpt2Signal","sumpt2 of simulated tracks in Signal events",100,0.,200.));
403  add(h,new TH1F("sumpt2PU","sumpt2 of simulated tracks in PU events",100,0.,200.));
404  add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
405  add(h,new TH1F("sumpt2recSignal","sumpt2 of reconstructed and matched tracks in Signal events",100,0.,200.));
406  add(h,new TH1F("sumpt2recPU","sumpt2 of reconstructed and matched tracks in PU events",100,0.,200.));
407  add(h,new TH1F("nRecTrkInSimVtx","number of reco tracks matched to sim-vertex", 101, 0., 100));
408  add(h,new TH1F("nRecTrkInSimVtxSignal","number of reco tracks matched to signal sim-vertex", 101, 0., 100));
409  add(h,new TH1F("nRecTrkInSimVtxPU","number of reco tracks matched to PU-vertex", 101, 0., 100));
410  add(h,new TH1F("nPrimRecTrkInSimVtx","number of reco primary tracks matched to sim-vertex", 101, 0., 100));
411  add(h,new TH1F("nPrimRecTrkInSimVtxSignal","number of reco primary tracks matched to signal sim-vertex", 101, 0., 100));
412  add(h,new TH1F("nPrimRecTrkInSimVtxPU","number of reco primary tracks matched to PU-vertex", 101, 0., 100));
413  // obsolete, scheduled for removal
414  add(h,new TH1F("recPurity","track purity of reconstructed vertices", 101, 0., 1.01));
415  add(h,new TH1F("recPuritySignal","track purity of reconstructed Signal vertices", 101, 0., 1.01));
416  add(h,new TH1F("recPurityPU","track purity of reconstructed PU vertices", 101, 0., 1.01));
417  // end of obsolete
418  add(h,new TH1F("recmatchPurity","track purity of all vertices", 101, 0., 1.01));
419  add(h,new TH1F("recmatchPurityTag","track purity of tagged vertices", 101, 0., 1.01));
420  add(h,new TH1F("recmatchPurityTagSignal","track purity of tagged Signal vertices", 101, 0., 1.01));
421  add(h,new TH1F("recmatchPurityTagPU","track purity of tagged PU vertices", 101, 0., 1.01));
422  add(h,new TH1F("recmatchPuritynoTag","track purity of untagged vertices", 101, 0., 1.01));
423  add(h,new TH1F("recmatchPuritynoTagSignal","track purity of untagged Signal vertices", 101, 0., 1.01));
424  add(h,new TH1F("recmatchPuritynoTagPU","track purity of untagged PU vertices", 101, 0., 1.01));
425  add(h,new TH1F("recmatchvtxs","number of sim vertices contributing to recvtx", 10, 0., 10.));
426  add(h,new TH1F("recmatchvtxsTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
427  add(h,new TH1F("recmatchvtxsnoTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
428  //
429  add(h,new TH1F("trkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
430  add(h,new TH1F("trkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
431  add(h,new TH1F("trkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
432  add(h,new TH1F("primtrkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
433  add(h,new TH1F("primtrkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
434  add(h,new TH1F("primtrkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
435  add(h,new TH1F("vtxMultiplicity", "number of rec vertices containg tracks from one true vertex", 10, 0., 10.) );
436  add(h,new TH1F("vtxMultiplicitySignal", "number of rec vertices containg tracks from the Signal Vertex", 10, 0., 10.) );
437  add(h,new TH1F("vtxMultiplicityPU", "number of rec vertices containg tracks from a PU Vertex", 10, 0., 10.) );
438 
439  add(h,new TProfile("vtxFindingEfficiencyVsNtrk","finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
440  add(h,new TProfile("vtxFindingEfficiencyVsNtrkSignal","Signal vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
441  add(h,new TProfile("vtxFindingEfficiencyVsNtrkPU","PU vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
442 
443  add(h,new TH1F("TagVtxTrkPurity","TagVtxTrkPurity",100,0.,1.01));
444  add(h,new TH1F("TagVtxTrkEfficiency","TagVtxTrkEfficiency",100,0.,1.01));
445 
446  add(h,new TH1F("matchVtxFraction","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
447  add(h,new TH1F("matchVtxFractionSignal","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
448  add(h,new TH1F("matchVtxFractionPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
449  add(h,new TH1F("matchVtxFractionCum","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
450  add(h,new TH1F("matchVtxFractionCumSignal","fraction of sim vertexs track found in a recvertex",101,0,1.01));
451  add(h,new TH1F("matchVtxFractionCumPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
452  add(h,new TH1F("matchVtxEfficiency","efficiency for finding matching rec vertex",2,-0.5,1.5));
453  add(h,new TH1F("matchVtxEfficiencySignal","efficiency for finding matching rec vertex",2,-0.5,1.5));
454  add(h,new TH1F("matchVtxEfficiencyPU","efficiency for finding matching rec vertex",2,-0.5,1.5));
455  add(h,new TH1F("matchVtxEfficiency2","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
456  add(h,new TH1F("matchVtxEfficiency2Signal","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
457  add(h,new TH1F("matchVtxEfficiency2PU","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
458  add(h,new TH1F("matchVtxEfficiency5","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
459  add(h,new TH1F("matchVtxEfficiency5Signal","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
460  add(h,new TH1F("matchVtxEfficiency5PU","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
461 
462 
463  add(h,new TH1F("matchVtxEfficiencyZ","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
464  add(h,new TH1F("matchVtxEfficiencyZSignal","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
465  add(h,new TH1F("matchVtxEfficiencyZPU","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
466 
467  add(h,new TH1F("matchVtxEfficiencyZ1","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
468  add(h,new TH1F("matchVtxEfficiencyZ1Signal","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
469  add(h,new TH1F("matchVtxEfficiencyZ1PU","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
470 
471  add(h,new TH1F("matchVtxEfficiencyZ2","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
472  add(h,new TH1F("matchVtxEfficiencyZ2Signal","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
473  add(h,new TH1F("matchVtxEfficiencyZ2PU","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
474 
475  add(h,new TH1F("matchVtxZ","z distance to matched recvtx",100, -0.1, 0.1));
476  add(h,new TH1F("matchVtxZPU","z distance to matched recvtx",100, -0.1, 0.1));
477  add(h,new TH1F("matchVtxZSignal","z distance to matched recvtx",100, -0.1, 0.1));
478 
479  add(h,new TH1F("matchVtxZCum","z distance to matched recvtx",1001,0, 1.01));
480  add(h,new TH1F("matchVtxZCumSignal","z distance to matched recvtx",1001,0, 1.01));
481  add(h,new TH1F("matchVtxZCumPU","z distance to matched recvtx",1001,0, 1.01));
482 
483  add(h, new TH1F("unmatchedVtx","number of unmatched rec vertices (fakes)",10,0.,10.));
484  add(h, new TH1F("unmatchedVtxFrac","fraction of unmatched rec vertices (fakes)",1000,0.,1.0));
485  add(h, new TH1F("unmatchedVtxZ","z of unmached rec vertices (fakes)",100,-20., 20.));
486  add(h, new TH1F("unmatchedVtxNdof","ndof of unmatched rec vertices (fakes)", 500,0., 100.));
487  add(h,new TH2F("correctlyassigned","pt and eta of correctly assigned tracks", 60, -3., 3., 100, 0, 10.));
488  add(h,new TH2F("misassigned","pt and eta of mis assigned tracks", 60, -3., 3., 100, 0, 10.));
489 
490  add(h,new TH1F("ptcat","pt of correctly assigned tracks", 100, 0, 10.));
491  add(h,new TH1F("etacat","eta of correctly assigned tracks", 60, -3., 3.));
492  add(h,new TH1F("phicat","phi of correctly assigned tracks", 100, -3.14159, 3.14159));
493  add(h,new TH1F("dzcat","dz of correctly assigned tracks", 100, 0., 1.));
494 
495  add(h,new TH1F("ptmis","pt of mis-assigned tracks", 100, 0, 10.));
496  add(h,new TH1F("etamis","eta of mis-assigned tracks", 60, -3., 3.));
497  add(h,new TH1F("phimis","phi of mis-assigned tracks",100, -3.14159, 3.14159));
498  add(h,new TH1F("dzmis","dz of mis-assigned tracks", 100, 0., 1.));
499 
500 
501  add(h,new TH1F("Tc","Tc computed with Truth matched Reco Tracks",100,0.,20.));
502  add(h,new TH1F("TcSignal","Tc of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
503  add(h,new TH1F("TcPU","Tc of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
504 
505  add(h,new TH1F("logTc","log Tc computed with Truth matched Reco Tracks",100,-2.,8.));
506  add(h,new TH1F("logTcSignal","log Tc of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
507  add(h,new TH1F("logTcPU","log Tc of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
508 
509  add(h,new TH1F("xTc","Tc of merged clusters",100,0.,20.));
510  add(h,new TH1F("xTcSignal","Tc of signal vertices merged with PU",100,0.,20.));
511  add(h,new TH1F("xTcPU","Tc of merged PU vertices",100,0.,20.));
512 
513  add(h,new TH1F("logxTc","log Tc merged vertices",100,-2.,8.));
514  add(h,new TH1F("logxTcSignal","log Tc of signal vertices merged with PU",100,-2.,8.));
515  add(h,new TH1F("logxTcPU","log Tc of merged PU vertices ",100,-2.,8.));
516 
517  add(h,new TH1F("logChisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,-2.,8.));
518  add(h,new TH1F("logChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
519  add(h,new TH1F("logChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
520 
521  add(h,new TH1F("logxChisq","Chisq/ntrk of merged clusters",100,-2.,8.));
522  add(h,new TH1F("logxChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,-2.,8.));
523  add(h,new TH1F("logxChisqPU","Chisq/ntrk of merged PU vertices",100,-2.,8.));
524 
525  add(h,new TH1F("Chisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,0.,20.));
526  add(h,new TH1F("ChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
527  add(h,new TH1F("ChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
528 
529  add(h,new TH1F("xChisq","Chisq/ntrk of merged clusters",100,0.,20.));
530  add(h,new TH1F("xChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,0.,20.));
531  add(h,new TH1F("xChisqPU","Chisq/ntrk of merged PU vertices",100,0.,20.));
532 
533  add(h,new TH1F("dzmax","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
534  add(h,new TH1F("dzmaxSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
535  add(h,new TH1F("dzmaxPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
536 
537  add(h,new TH1F("xdzmax","dzmax of merged clusters",100,0.,2.));
538  add(h,new TH1F("xdzmaxSignal","dzmax of signal vertices merged with PU",100,0.,2.));
539  add(h,new TH1F("xdzmaxPU","dzmax of merged PU vertices",100,0.,2.));
540 
541  add(h,new TH1F("dztrim","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
542  add(h,new TH1F("dztrimSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
543  add(h,new TH1F("dztrimPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
544 
545  add(h,new TH1F("xdztrim","dzmax of merged clusters",100,0.,2.));
546  add(h,new TH1F("xdztrimSignal","dzmax of signal vertices merged with PU",100,0.,2.));
547  add(h,new TH1F("xdztrimPU","dzmax of merged PU vertices",100,0.,2.));
548 
549  add(h,new TH1F("m4m2","m4m2 computed with Truth matched Reco Tracks",100,0.,100.));
550  add(h,new TH1F("m4m2Signal","m4m2 of signal vertices computed with Truth matched Reco Tracks",100,0.,100.));
551  add(h,new TH1F("m4m2PU","m4m2 of PU vertices computed with Truth matched Reco Tracks",100,0.,100.));
552 
553  add(h,new TH1F("xm4m2","m4m2 of merged clusters",100,0.,100.));
554  add(h,new TH1F("xm4m2Signal","m4m2 of signal vertices merged with PU",100,0.,100.));
555  add(h,new TH1F("xm4m2PU","m4m2 of merged PU vertices",100,0.,100.));
556 
557  return h;
558 }
559 
560 
561 //
562 // member functions
563 //
564 void PrimaryVertexAnalyzer4PU::beginJob(){
565  std::cout << " PrimaryVertexAnalyzer4PU::beginJob conversion from sim units to rec units is " << simUnit_ << std::endl;
566 
567 
568  rootFile_->cd();
569 
570  TDirectory *noBS = rootFile_->mkdir("noBS");
571  noBS->cd();
572  hnoBS=bookVertexHistograms();
573  for(std::map<std::string,TH1*>::const_iterator hist=hnoBS.begin(); hist!=hnoBS.end(); hist++){
574  hist->second->SetDirectory(noBS);
575  }
576 
577  TDirectory *withBS = rootFile_->mkdir("BS");
578  withBS->cd();
579  hBS=bookVertexHistograms();
580  for(std::map<std::string,TH1*>::const_iterator hist=hBS.begin(); hist!=hBS.end(); hist++){
581  hist->second->SetDirectory(withBS);
582  }
583 
584  TDirectory *DA = rootFile_->mkdir("DA");
585  DA->cd();
586  hDA=bookVertexHistograms();
587  for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
588  hist->second->SetDirectory(DA);
589  }
590 
591 // TDirectory *PIX = rootFile_->mkdir("PIX");
592 // PIX->cd();
593 // hPIX=bookVertexHistograms();
594 // for(std::map<std::string,TH1*>::const_iterator hist=hPIX.begin(); hist!=hPIX.end(); hist++){
595 // hist->second->SetDirectory(PIX);
596 // }
597 
598 // TDirectory *MVF = rootFile_->mkdir("MVF");
599 // MVF->cd();
600 // hMVF=bookVertexHistograms();
601 // for(std::map<std::string,TH1*>::const_iterator hist=hMVF.begin(); hist!=hMVF.end(); hist++){
602 // hist->second->SetDirectory(MVF);
603 // }
604 
605  rootFile_->cd();
606  hsimPV["rapidity"] = new TH1F("rapidity","rapidity ",100,-10., 10.);
607  hsimPV["chRapidity"] = new TH1F("chRapidity","charged rapidity ",100,-10., 10.);
608  hsimPV["recRapidity"] = new TH1F("recRapidity","reconstructed rapidity ",100,-10., 10.);
609  hsimPV["pt"] = new TH1F("pt","pt ",100,0., 20.);
610 
611  hsimPV["xsim"] = new TH1F("xsim","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
612  hsimPV["ysim"] = new TH1F("ysim","simulated y",100,-0.01,0.01);
613  hsimPV["zsim"] = new TH1F("zsim","simulated z",100,-20.,20.);
614 
615  hsimPV["xsim1"] = new TH1F("xsim1","simulated x",100,-4.,4.);
616  hsimPV["ysim1"] = new TH1F("ysim1","simulated y",100,-4.,4.);
617  hsimPV["zsim1"] = new TH1F("zsim1","simulated z",100,-40.,40.);
618 
619  add(hsimPV, new TH1F("xsim2PU","simulated x (Pile-up)",100,-1.,1.));
620  add(hsimPV, new TH1F("ysim2PU","simulated y (Pile-up)",100,-1.,1.));
621  add(hsimPV, new TH1F("zsim2PU","simulated z (Pile-up)",100,-20.,20.));
622  add(hsimPV, new TH1F("xsim2Signal","simulated x (Signal)",100,-1.,1.));
623  add(hsimPV, new TH1F("ysim2Signal","simulated y (Signal)",100,-1.,1.));
624  add(hsimPV, new TH1F("zsim2Signal","simulated z (Signal)",100,-20.,20.));
625 
626  hsimPV["xsim2"] = new TH1F("xsim2","simulated x",100,-1,1); // 0.01cm = 100 um
627  hsimPV["ysim2"] = new TH1F("ysim2","simulated y",100,-1,1);
628  hsimPV["zsim2"] = new TH1F("zsim2","simulated z",100,-20.,20.);
629  hsimPV["xsim3"] = new TH1F("xsim3","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
630  hsimPV["ysim3"] = new TH1F("ysim3","simulated y",100,-0.1,0.1);
631  hsimPV["zsim3"] = new TH1F("zsim3","simulated z",100,-20.,20.);
632  hsimPV["xsimb"] = new TH1F("xsimb","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
633  hsimPV["ysimb"] = new TH1F("ysimb","simulated y",100,-0.01,0.01);
634  hsimPV["zsimb"] = new TH1F("zsimb","simulated z",100,-20.,20.);
635  hsimPV["xsimb1"] = new TH1F("xsimb1","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
636  hsimPV["ysimb1"] = new TH1F("ysimb1","simulated y",100,-0.1,0.1);
637  hsimPV["zsimb1"] = new TH1F("zsimb1","simulated z",100,-20.,20.);
638  add(hsimPV, new TH1F("xbeam","beamspot x",100,-1.,1.));
639  add(hsimPV, new TH1F("ybeam","beamspot y",100,-1.,1.));
640  add(hsimPV, new TH1F("zbeam","beamspot z",100,-1.,1.));
641  add(hsimPV, new TH1F("wxbeam","beamspot sigma x",100,-1.,1.));
642  add(hsimPV, new TH1F("wybeam","beamspot sigma y",100,-1.,1.));
643  add(hsimPV, new TH1F("wzbeam","beamspot sigma z",100,-1.,1.));
644  hsimPV["xsim2"]->StatOverflows(kTRUE);
645  hsimPV["ysim2"]->StatOverflows(kTRUE);
646  hsimPV["zsim2"]->StatOverflows(kTRUE);
647  hsimPV["xsimb"]->StatOverflows(kTRUE);
648  hsimPV["ysimb"]->StatOverflows(kTRUE);
649  hsimPV["zsimb"]->StatOverflows(kTRUE);
650  hsimPV["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
651  // hsimPV["nsimtrk"] = new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5); // not filled right now, also exists in hBS..
652  // hsimPV["nsimtrk"]->StatOverflows(kTRUE);
653  hsimPV["nbsimtksinvtx"]= new TH1F("nbsimtksinvtx","simulated tracks in vertex",100,-0.5,99.5);
654  hsimPV["nbsimtksinvtx"]->StatOverflows(kTRUE);
655 
656 }
657 
658 
659 void PrimaryVertexAnalyzer4PU::endJob() {
660  std::cout << "this is void PrimaryVertexAnalyzer4PU::endJob() " << std::endl;
661  //cumulate some histos
662  double sumDA=0,sumBS=0,sumnoBS=0, sumPIX=0,sumMVF=0;
663  for(int i=101; i>0; i--){
664  sumDA+=hDA["matchVtxFractionSignal"]->GetBinContent(i)/hDA["matchVtxFractionSignal"]->Integral();
665  hDA["matchVtxFractionCumSignal"]->SetBinContent(i,sumDA);
666  sumBS+=hBS["matchVtxFractionSignal"]->GetBinContent(i)/hBS["matchVtxFractionSignal"]->Integral();
667  hBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumBS);
668  sumnoBS+=hnoBS["matchVtxFractionSignal"]->GetBinContent(i)/hnoBS["matchVtxFractionSignal"]->Integral();
669  hnoBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumnoBS);
670 // sumPIX+=hPIX["matchVtxFractionSignal"]->GetBinContent(i)/hPIX["matchVtxFractionSignal"]->Integral();
671 // hPIX["matchVtxFractionCumSignal"]->SetBinContent(i,sumPIX);
672 // sumMVF+=hMVF["matchVtxFractionSignal"]->GetBinContent(i)/hMVF["matchVtxFractionSignal"]->Integral();
673 // hMVF["matchVtxFractionCumSignal"]->SetBinContent(i,sumMVF);
674  }
675  sumDA=0,sumBS=0,sumnoBS=0,sumPIX=0,sumMVF=0;
676  for(int i=1; i<1001; i++){
677  sumDA+=hDA["abszdistancetag"]->GetBinContent(i);
678  hDA["abszdistancetagcum"]->SetBinContent(i,sumDA/float(hDA["abszdistancetag"]->GetEntries()));
679  sumBS+=hBS["abszdistancetag"]->GetBinContent(i);
680  hBS["abszdistancetagcum"]->SetBinContent(i,sumBS/float(hBS["abszdistancetag"]->GetEntries()));
681  sumnoBS+=hnoBS["abszdistancetag"]->GetBinContent(i);
682  hnoBS["abszdistancetagcum"]->SetBinContent(i,sumnoBS/float(hnoBS["abszdistancetag"]->GetEntries()));
683 // sumPIX+=hPIX["abszdistancetag"]->GetBinContent(i);
684 // hPIX["abszdistancetagcum"]->SetBinContent(i,sumPIX/float(hPIX["abszdistancetag"]->GetEntries()));
685 // sumMVF+=hMVF["abszdistancetag"]->GetBinContent(i);
686 // hMVF["abszdistancetagcum"]->SetBinContent(i,sumMVF/float(hMVF["abszdistancetag"]->GetEntries()));
687  }
688 
689  Cumulate(hBS["matchVtxZCum"]); Cumulate(hBS["matchVtxZCumSignal"]); Cumulate(hBS["matchVtxZCumPU"]);
690  Cumulate(hnoBS["matchVtxZCum"]); Cumulate(hnoBS["matchVtxZCumSignal"]); Cumulate(hnoBS["matchVtxZCumPU"]);
691  Cumulate(hDA["matchVtxZCum"]); Cumulate(hDA["matchVtxZCumSignal"]); Cumulate(hDA["matchVtxZCumPU"]);
692  /*
693  h->ComputeIntegral();
694  Double_t *integral = h->GetIntegral();
695  h->SetContent(integral);
696  */
697 
698  // make a reference for ndofnr2
699  //hDA["vtxndof"]->ComputeIntegral();
700  //Double_t *integral = hDA["vtxndf"]->GetIntegral();
701  // h->SetContent(integral);
702  double p;
703  for(int i=1; i<501; i++){
704  if(hDA["vtxndf"]->GetEntries()>0){
705  p= hDA["vtxndf"]->Integral(i,501)/hDA["vtxndf"]->GetEntries(); hDA["vtxndfc"]->SetBinContent(i,p*hDA["vtxndf"]->GetBinContent(i));
706  }
707  if(hBS["vtxndf"]->GetEntries()>0){
708  p= hBS["vtxndf"]->Integral(i,501)/hBS["vtxndf"]->GetEntries(); hBS["vtxndfc"]->SetBinContent(i,p*hBS["vtxndf"]->GetBinContent(i));
709  }
710  if(hnoBS["vtxndf"]->GetEntries()>0){
711  p=hnoBS["vtxndf"]->Integral(i,501)/hnoBS["vtxndf"]->GetEntries(); hnoBS["vtxndfc"]->SetBinContent(i,p*hnoBS["vtxndf"]->GetBinContent(i));
712  }
713 // if(hPIX["vtxndf"]->GetEntries()>0){
714 // p=hPIX["vtxndf"]->Integral(i,501)/hPIX["vtxndf"]->GetEntries(); hPIX["vtxndfc"]->SetBinContent(i,p*hPIX["vtxndf"]->GetBinContent(i));
715 // }
716  }
717 
718  rootFile_->cd();
719  for(std::map<std::string,TH1*>::const_iterator hist=hsimPV.begin(); hist!=hsimPV.end(); hist++){
720  std::cout << "writing " << hist->first << std::endl;
721  hist->second->Write();
722  }
723  rootFile_->Write();
724  std::cout << "PrimaryVertexAnalyzer4PU::endJob: done" << std::endl;
725 }
726 
727 
728 
729 
730 // helper functions
731 std::vector<PrimaryVertexAnalyzer4PU::SimPart> PrimaryVertexAnalyzer4PU::getSimTrkParameters(
732  edm::Handle<edm::SimTrackContainer> & simTrks,
733  edm::Handle<edm::SimVertexContainer> & simVtcs,
734  double simUnit)
735 {
736  std::vector<SimPart > tsim;
737  if(simVtcs->begin()==simVtcs->end()){
738  if(verbose_){
739  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters no simvtcs" << endl;
740  }
741  return tsim;
742  }
743  if(verbose_){
744  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters simVtcs n=" << simVtcs->size() << endl;
745  cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters 1st position" << setw(8) << setprecision(4) << simVtcs->begin()->position() << endl;
746  }
747  double t0=simVtcs->begin()->position().e();
748 
749  for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
750  t!=simTrks->end(); ++t){
751  if (t->noVertex()){
752  std::cout << "simtrk has no vertex" << std::endl;
753  }else{
754  // get the vertex position
755  //HepLorentzVector v=(*simVtcs)[t->vertIndex()].position();
756  math::XYZTLorentzVectorD v((*simVtcs)[t->vertIndex()].position().x(),
757  (*simVtcs)[t->vertIndex()].position().y(),
758  (*simVtcs)[t->vertIndex()].position().z(),
759  (*simVtcs)[t->vertIndex()].position().e());
760  int pdgCode=t->type();
761 
762  if( pdgCode==-99 ){
763  // such entries cause crashes, no idea what they are
764  std::cout << "funny particle skipped , code=" << pdgCode << std::endl;
765  }else{
766  double Q=0; //double Q=HepPDT::theTable().getParticleData(pdgCode)->charge();
767  if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==-321)||(pdgCode==-3222)){Q=-1;}
768  else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)||(pdgCode==3222)){Q=1;}
769  else {
770  //std::cout << pdgCode << " " <<std::endl;
771  }
772  math::XYZTLorentzVectorD p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
773  if ( (Q != 0) && (p.pt()>0.1) && (fabs(t->momentum().eta())<3.0)
774  && fabs(v.z()*simUnit<20) && (sqrt(v.x()*v.x()+v.y()*v.y())<10.)){
775  double x0=v.x()*simUnit;
776  double y0=v.y()*simUnit;
777  double z0=v.z()*simUnit;
778  double kappa=-Q*0.002998*fBfield_/p.pt();
779  double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
780  double q=sqrt(1.-2.*kappa*D0);
781  double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
782  double s1;
783  if (fabs(kappa*s0)>0.001){
784  s1=asin(kappa*s0)/kappa;
785  }else{
786  double ks02=(kappa*s0)*(kappa*s0);
787  s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
788  }
789  SimPart sp;//ParameterVector par;
790  sp.par[reco::TrackBase::i_qoverp] = Q/p.P();
791  sp.par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
792  sp.par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
793  sp.par[reco::TrackBase::i_dxy] = -2.*D0/(1.+q);
794  sp.par[reco::TrackBase::i_dsz] = z0*sin(p.theta())-s1*cos(p.theta());
795 
796  sp.pdg=pdgCode;
797  if (v.t()-t0<1e-15){
798  sp.type=0; // primary
799  }else{
800  sp.type=1; //secondary
801  }
802 
803  // now get zpca (get perigee wrt beam)
804  double x1=x0-0.033; double y1=y0-0.; // FIXME how do we get the simulated beam position?
805  D0=x1*sin(p.phi())-y1*cos(p.phi())-0.5*kappa*(x1*x1+y1*y1);
806  q=sqrt(1.-2.*kappa*D0);
807  s0=(x1*cos(p.phi())+y1*sin(p.phi()))/q;
808  if (fabs(kappa*s0)>0.001){
809  s1=asin(kappa*s0)/kappa;
810  }else{
811  double ks02=(kappa*s0)*(kappa*s0);
812  s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
813  }
814  sp.ddcap=-2.*D0/(1.+q);
815  sp.zdcap=z0-s1/tan(p.theta());
816  sp.zvtx=z0;
817  sp.xvtx=x0;
818  sp.yvtx=y0;
819 
820  tsim.push_back(sp);
821  }
822  }
823  }// has vertex
824  }//for loop
825  return tsim;
826 }
827 
828 
829 int* PrimaryVertexAnalyzer4PU::supf(std::vector<SimPart>& simtrks, const reco::TrackCollection & trks){
830  int nsim=simtrks.size();
831  int nrec=trks.size();
832  int *rectosim=new int[nrec]; // pointer to associated simtrk
833  double** pij=new double*[nrec];
834  double mu=100.; // initial chi^2 cut-off (5 dofs !)
835  int nmatch=0;
836  int i=0;
837  for(reco::TrackCollection::const_iterator t=trks.begin(); t!=trks.end(); ++t){
838  pij[i]=new double[nsim];
839  rectosim[i]=-1;
840  ParameterVector par = t->parameters();
841  //reco::TrackBase::CovarianceMatrix V = t->covariance();
842  reco::TrackBase::CovarianceMatrix S = t->covariance();
843  S.Invert();
844  for(int j=0; j<nsim; j++){
845  simtrks[j].rec=-1;
846  SimPart s=simtrks[j];
847  double c=0;
848  for(int k=0; k<5; k++){
849  for(int l=0; l<5; l++){
850  c+=(par(k)-s.par[k])*(par(l)-s.par[l])*S(k,l);
851  }
852  }
853  pij[i][j]=exp(-0.5*c);
854 
855 // double c0=pow((par[0]-s.par[0])/t->qoverpError(),2)*0.1
856 // +pow((par[1]-s.par[1])/t->lambdaError(),2)
857 // +pow((par[2]-s.par[2])/t->phiError(),2)
858 // +pow((par[3]-s.par[3])/t->dxyError(),2)*0.1;
859 // +pow((par[4]-s.par[4])/t->dszError(),2)*0.1;
860 // pij[i][j]=exp(-0.5*c0);
861 
862 // if( c0 <100 ){
863 // cout << setw(3) << i << " rec " << setw(6) << par << endl;
864 // cout << setw(3) << j << " sim " << setw(6) << s.par << " ---> C=" << c << endl;
865 // cout << " " << setw(6)
866 // << (par[0]-s.par[0])<< ","
867 // << (par[1]-s.par[1])<< ","
868 // << (par[2]-s.par[2])<< ","
869 // << (par[3]-s.par[3])<< ","
870 // << (par[4]-s.par[4])
871 // << " match=" << match(simtrks[j].par, trks.at(i).parameters())
872 // << endl;
873 // cout << " " << setw(6)
874 // << (par[0]-s.par[0])/t->qoverpError() << ","
875 // << (par[1]-s.par[1])/t->lambdaError() << ","
876 // << (par[2]-s.par[2])/t->phiError() << ","
877 // << (par[3]-s.par[3])/t->dxyError() << ","
878 // << (par[4]-s.par[4])/t->dszError() << " c0=" << c0
879 // << endl <<endl;
880 // }
881 
882  }
883  i++;
884  }
885 
886  for(int k=0; k<nrec; k++){
887  int imatch=-1; int jmatch=-1;
888  double pmatch=0;
889  for(int j=0; j<nsim; j++){
890  if ((simtrks[j].rec)<0){
891  double psum=exp(-0.5*mu); //cutoff
892  for(int i=0; i<nrec; i++){
893  if (rectosim[i]<0){ psum+=pij[i][j];}
894  }
895  for(int i=0; i<nrec; i++){
896  if ((rectosim[i]<0)&&(pij[i][j]/psum>pmatch)){
897  pmatch=pij[i][j]/psum;
898  imatch=i; jmatch=j;
899  }
900  }
901  }
902  }// sim loop
903  if((jmatch>=0)||(imatch>=0)){
904  //std::cout << pmatch << " " << pij[imatch][jmatch] << " match=" <<
905  // match(simtrks[jmatch].par, trks.at(imatch).parameters()) <<std::endl;
906  if (pmatch>0.01){
907  rectosim[imatch]=jmatch;
908  simtrks[jmatch].rec=imatch;
909  nmatch++;
910  }else if (match(simtrks[jmatch].par, trks.at(imatch).parameters())){
911  // accept it anyway if it matches crudely and relax the cut-off
912  rectosim[imatch]=jmatch;
913  simtrks[jmatch].rec=imatch;
914  nmatch++;
915  mu=mu*2;
916  }
917  }
918  }
919 
920 // std::cout << ">>>>>>>>>>>>>>>--------------supf----------------------" << std::endl;
921 // std::cout <<"nsim=" << nsim << " nrec=" << nrec << " nmatch=" << nmatch << std::endl;
922 // std::cout << "rec to sim " << std::endl;
923 // for(int i=0; i<nrec; i++){
924 // std::cout << i << " ---> " << rectosim[i] << std::endl;
925 // }
926 // std::cout << "sim to rec " << std::endl;
927 // for(int j=0; j<nsim; j++){
928 // std::cout << j << " ---> " << simtrks[j].rec << std::endl;
929 // }
930 
931  std::cout << "unmatched sim " << std::endl;
932  for(int j=0; j<nsim; j++){
933  if(simtrks[j].rec<0){
934  double pt= 1./simtrks[j].par[0]/tan(simtrks[j].par[1]);
935  if((fabs(pt))>1.){
936  std::cout << setw(3) << j << setw(8) << simtrks[j].pdg
937  << setw(8) << setprecision(4) << " (" << simtrks[j].xvtx << "," << simtrks[j].yvtx << "," << simtrks[j].zvtx << ")"
938  << " pt= " << pt
939  << " phi=" << simtrks[j].par[2]
940  << " eta= " << -log(tan(0.5*(M_PI/2-simtrks[j].par[1])))
941  << std::endl;
942  }
943  }
944  }
945 // std::cout << "<<<<<<<<<<<<<<<--------------supf----------------------" << std::endl;
946 
947  //delete rectosim; // or return it?
948  for(int i=0; i<nrec; i++){delete pij[i];}
949  delete pij;
950  return rectosim; // caller must delete it
951 }
952 
953 
954 
955 
956 
957 
958 
959 
960 bool PrimaryVertexAnalyzer4PU::match(const ParameterVector &a,
961  const ParameterVector &b){
962  double dqoverp =a(0)-b(0);
963  double dlambda =a(1)-b(1);
964  double dphi =a(2)-b(2);
965  double dsz =a(4)-b(4);
966  if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
967  // return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<0.1) );
968  return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<1.0) );
969 }
970 
971 
972 bool PrimaryVertexAnalyzer4PU::matchVertex(const simPrimaryVertex &vsim,
973  const reco::Vertex &vrec){
974  return (fabs(vsim.z*simUnit_-vrec.z())<zmatch_);
975 }
976 
977 bool PrimaryVertexAnalyzer4PU::isResonance(const HepMC::GenParticle * p){
978  double ctau=(pdt_->particle( abs(p->pdg_id()) ))->lifetime();
979  //std::cout << "isResonance " << p->pdg_id() << " " << ctau << std::endl;
980  return ctau >0 && ctau <1e-6;
981 }
982 
983 bool PrimaryVertexAnalyzer4PU::isFinalstateParticle(const HepMC::GenParticle * p){
984  return ( !p->end_vertex() && p->status()==1 );
985 }
986 
987 
988 bool PrimaryVertexAnalyzer4PU::isCharged(const HepMC::GenParticle * p){
989  const ParticleData * part = pdt_->particle( p->pdg_id() );
990  if (part){
991  return part->charge()!=0;
992  }else{
993  // the new/improved particle table doesn't know anti-particles
994  return pdt_->particle( -p->pdg_id() )!=0;
995  }
996 }
997 
998 
999 
1000 
1001 void PrimaryVertexAnalyzer4PU::fillTrackHistos(std::map<std::string, TH1*> & h, const std::string & ttype, const reco::Track & t, const reco::Vertex * v){
1002  Fill(h,"rapidity_"+ttype,t.eta());
1003  Fill(h,"z0_"+ttype,t.vz());
1004  Fill(h,"phi_"+ttype,t.phi());
1005  Fill(h,"eta_"+ttype,t.eta());
1006  Fill(h,"pt_"+ttype,t.pt());
1007  Fill(h,"found_"+ttype,t.found());
1008  Fill(h,"lost_"+ttype,t.lost());
1009  Fill(h,"nchi2_"+ttype,t.normalizedChi2());
1010  Fill(h,"rstart_"+ttype,(t.innerPosition()).Rho());
1011 
1012  double d0Error=t.d0Error();
1013  double d0=t.dxy(vertexBeamSpot_.position());
1014  if (d0Error>0){
1015  Fill(h,"logtresxy_"+ttype,log(d0Error/0.0001)/log(10.));
1016  Fill(h,"tpullxy_"+ttype,d0/d0Error);
1017  }
1018  //double z0=t.vz();
1019  double dzError=t.dzError();
1020  if(dzError>0){
1021  Fill(h,"logtresz_"+ttype,log(dzError/0.0001)/log(10.));
1022  }
1023 
1024  //
1025  Fill(h,"sigmatrkz_"+ttype,sqrt(pow(t.dzError(),2)+wxy2_/pow(tan(t.theta()),2)));
1026  Fill(h,"sigmatrkz0_"+ttype,t.dzError());
1027 
1028  // track vs vertex
1029  if((! (v==NULL)) && (v->ndof()>10.)) {
1030  // emulate clusterizer input
1031  //const TransientTrack & tt = theB_->build(&t); wrong !!!!
1032  TransientTrack tt = theB_->build(&t); tt.setBeamSpot(vertexBeamSpot_); // need the setBeamSpot !
1033  double z=(tt.stateAtBeamLine().trackStateAtPCA()).position().z();
1034  double tantheta=tan((tt.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1035  double dz2= pow(tt.track().dzError(),2)+wxy2_/pow(tantheta,2);
1036 
1037  Fill(h,"restrkz_"+ttype,z-v->position().z());
1038  Fill(h,"restrkzvsphi_"+ttype,t.phi(), z-v->position().z());
1039  Fill(h,"restrkzvseta_"+ttype,t.eta(), z-v->position().z());
1040  Fill(h,"pulltrkzvsphi_"+ttype,t.phi(), (z-v->position().z())/sqrt(dz2));
1041  Fill(h,"pulltrkzvseta_"+ttype,t.eta(), (z-v->position().z())/sqrt(dz2));
1042 
1043  Fill(h,"pulltrkz_"+ttype,(z-v->position().z())/sqrt(dz2));
1044 
1045  double x1=t.vx()-vertexBeamSpot_.x0(); double y1=t.vy()-vertexBeamSpot_.y0();
1046  double kappa=-0.002998*fBfield_*t.qoverp()/cos(t.theta());
1047  double D0=x1*sin(t.phi())-y1*cos(t.phi())-0.5*kappa*(x1*x1+y1*y1);
1048  double q=sqrt(1.-2.*kappa*D0);
1049  double s0=(x1*cos(t.phi())+y1*sin(t.phi()))/q;
1050  double s1;
1051  if (fabs(kappa*s0)>0.001){
1052  s1=asin(kappa*s0)/kappa;
1053  }else{
1054  double ks02=(kappa*s0)*(kappa*s0);
1055  s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
1056  }
1057  // sp.ddcap=-2.*D0/(1.+q);
1058  //double zdcap=t.vz()-s1/tan(t.theta());
1059 
1060  }
1061  //
1062 
1063  // collect some info on hits and clusters
1064  Fill(h,"nbarrelLayers_"+ttype,t.hitPattern().pixelBarrelLayersWithMeasurement());
1065  Fill(h,"nPxLayers_"+ttype,t.hitPattern().pixelLayersWithMeasurement());
1066  Fill(h,"nSiLayers_"+ttype,t.hitPattern().trackerLayersWithMeasurement());
1067  Fill(h,"expectedInner_"+ttype,t.trackerExpectedHitsInner().numberOfHits());
1068  Fill(h,"expectedOuter_"+ttype,t.trackerExpectedHitsOuter().numberOfHits());
1069  Fill(h,"trackAlgo_"+ttype,t.algo());
1070  Fill(h,"trackQuality_"+ttype,t.qualityMask());
1071 
1072  //
1073  int longesthit=0, nbarrel=0;
1074  for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
1075  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1076  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1077  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1078  if (barrel){
1079  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1080  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1081  if (clust.isNonnull()) {
1082  nbarrel++;
1083  if (clust->sizeY()>longesthit) longesthit=clust->sizeY();
1084  if (clust->sizeY()>20.){
1085  Fill(h,"lvseta_"+ttype,t.eta(), 19.9);
1086  Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), 19.9);
1087  }else{
1088  Fill(h,"lvseta_"+ttype,t.eta(), float(clust->sizeY()));
1089  Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), float(clust->sizeY()));
1090  }
1091  }
1092  }
1093  }
1094  }
1095  Fill(h,"nbarrelhits_"+ttype,float(nbarrel));
1096  //-------------------------------------------------------------------
1097 
1098 }
1099 
1100 void PrimaryVertexAnalyzer4PU::dumpHitInfo(const reco::Track & t){
1101  // collect some info on hits and clusters
1102  int longesthit=0, nbarrel=0;
1103  cout << Form("%5.2f %5.2f : ",t.pt(),t.eta());
1104  for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
1105  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1106  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1107  //bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1108  if (barrel){
1109  nbarrel++;
1110  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1111  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1112  if (clust.isNonnull()) {
1113  cout << Form("%4d",clust->sizeY());
1114  if (clust->sizeY()>longesthit) longesthit=clust->sizeY();
1115  }
1116  }
1117  }
1118  }
1119  //cout << endl;
1120 }
1121 
1122 void PrimaryVertexAnalyzer4PU::printRecVtxs(const Handle<reco::VertexCollection> recVtxs, std::string title){
1123  int ivtx=0;
1124  std::cout << std::endl << title << std::endl;
1125  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1126  v!=recVtxs->end(); ++v){
1127  string vtype=" recvtx ";
1128  if( v->isFake()){
1129  vtype=" fake ";
1130  }else if (v->ndof()==-5){
1131  vtype=" cluster "; // pos=selector[iclu],cputime[iclu],clusterz[iclu]
1132  }else if(v->ndof()==-3){
1133  vtype=" event ";
1134  }
1135  std::cout << "vtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
1136  << vtype
1137  << " #trk " << std::fixed << std::setprecision(4) << std::setw(3) << v->tracksSize()
1138  << " chi2 " << std::fixed << std::setw(4) << std::setprecision(1) << v->chi2()
1139  << " ndof " << std::fixed << std::setw(5) << std::setprecision(2) << v->ndof()
1140  << " x " << std::setw(8) <<std::fixed << std::setprecision(4) << v->x()
1141  << " dx " << std::setw(8) << v->xError()// << std::endl
1142  << " y " << std::setw(8) << v->y()
1143  << " dy " << std::setw(8) << v->yError()//<< std::endl
1144  << " z " << std::setw(8) << v->z()
1145  << " dz " << std::setw(8) << v->zError()
1146  << std::endl;
1147  }
1148 }
1149 
1150 
1151 void PrimaryVertexAnalyzer4PU::printSimVtxs(const Handle<SimVertexContainer> simVtxs){
1152  int i=0;
1153  for(SimVertexContainer::const_iterator vsim=simVtxs->begin();
1154  vsim!=simVtxs->end(); ++vsim){
1155  if ( vsim->position().x()*vsim->position().x()+vsim->position().y()*vsim->position().y() < 1.){
1156  std::cout << i++ << ")" << std::scientific
1157  << " evtid=" << vsim->eventId().event() << "," << vsim->eventId().bunchCrossing()
1158  << " sim x=" << vsim->position().x()*simUnit_
1159  << " sim y=" << vsim->position().y()*simUnit_
1160  << " sim z=" << vsim->position().z()*simUnit_
1161  << " sim t=" << vsim->position().t()
1162  << " parent=" << vsim->parentIndex()
1163  << std::endl;
1164  }
1165  }
1166 }
1167 
1168 
1169 
1170 
1171 
1172 
1173 
1174 void PrimaryVertexAnalyzer4PU::printSimTrks(const Handle<SimTrackContainer> simTrks){
1175  std::cout << " simTrks type, (momentum), vertIndex, genpartIndex" << std::endl;
1176  int i=1;
1177  for(SimTrackContainer::const_iterator t=simTrks->begin();
1178  t!=simTrks->end(); ++t){
1179  //HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
1180  std::cout << i++ << ")"
1181  << t->eventId().event() << "," << t->eventId().bunchCrossing()
1182  << (*t)
1183  << " index="
1184  << (*t).genpartIndex();
1185  //if (gp) {
1186  // HepMC::GenVertex *gv=gp->production_vertex();
1187  // std::cout << " genvertex =" << (*gv);
1188  //}
1189  std::cout << std::endl;
1190  }
1191 }
1192 
1193 
1194 
1195 void PrimaryVertexAnalyzer4PU::printRecTrks(const Handle<reco::TrackCollection> &recTrks ){
1196 
1197  cout << "printRecTrks" << endl;
1198  int i =1;
1199  for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
1200  // reco::TrackBase::ParameterVector par = t->parameters();
1201 
1202  cout << endl;
1203  cout << "Track "<<i << " " ; i++;
1204  //enum TrackQuality { undefQuality=-1, loose=0, tight=1, highPurity=2, confirmed=3, goodIterative=4, qualitySize=5};
1205  if( t->quality(reco::TrackBase::loose)){ cout << "loose ";};
1206  if( t->quality(reco::TrackBase::tight)){ cout << "tight ";};
1207  if( t->quality(reco::TrackBase::highPurity)){ cout << "highPurity ";};
1208  if( t->quality(reco::TrackBase::confirmed)){ cout << "confirmed ";};
1209  if( t->quality(reco::TrackBase::goodIterative)){ cout << "goodIterative ";};
1210  cout << endl;
1211 
1212  TransientTrack tk = theB_->build(&(*t)); tk.setBeamSpot(vertexBeamSpot_);
1213  double ipsig=0;
1214  if (tk.stateAtBeamLine().isValid()){
1215  ipsig= tk.stateAtBeamLine().transverseImpactParameter().significance();
1216  }else{
1217  ipsig=-1;
1218  }
1219 
1220  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;
1221 
1222 
1223  cout << Form(" found=%6d lost=%6d chi2/ndof=%10.3f ",t->found(), t->lost(),t->normalizedChi2())<<endl;
1224  const reco::HitPattern & p= t->hitPattern();
1225  cout << "subdet layers valid lost" << endl;
1226  cout << Form(" barrel %2d %2d %2d",p.pixelBarrelLayersWithMeasurement(),p.numberOfValidPixelBarrelHits(), p.numberOfLostPixelBarrelHits()) << endl;
1227  cout << Form(" fwd %2d %2d %2d",p.pixelEndcapLayersWithMeasurement(),p.numberOfValidPixelEndcapHits(), p.numberOfLostPixelEndcapHits()) << endl;
1228  cout << Form(" pixel %2d %2d %2d",p.pixelLayersWithMeasurement(), p.numberOfValidPixelHits(), p.numberOfLostPixelHits()) << endl;
1229  cout << Form(" tracker %2d %2d %2d",p.trackerLayersWithMeasurement(), p.numberOfValidTrackerHits(), p.numberOfLostTrackerHits()) << endl;
1230  cout << endl;
1231  const reco::HitPattern & pinner= t->trackerExpectedHitsInner();
1232  const reco::HitPattern & pouter= t->trackerExpectedHitsOuter();
1233  int ninner=pinner.numberOfHits();
1234  int nouter=pouter.numberOfHits();
1235 
1236  // double d0Error=t.d0Error();
1237  // double d0=t.dxy(myBeamSpot);
1238 
1239  //
1240  for(trackingRecHit_iterator hit=t->recHitsBegin(); hit!=t->recHitsEnd(); hit++){
1241  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1242  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1243  bool endcap = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1244  if (barrel){
1245  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1246  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1247  if (clust.isNonnull()) {
1248  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;
1249  }
1250  }else if(endcap){
1251  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1252  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1253  if (clust.isNonnull()) {
1254  cout << Form(" endcap cluster size=%2d charge=%6.1f wx=%2d wy=%2d",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY()) << endl;
1255  }
1256  }
1257  }
1258  }
1259  cout << "hitpattern" << endl;
1260  for(int i=0; i<p.numberOfHits(); i++){ p.printHitPattern(i,std::cout); }
1261  cout << "expected inner " << ninner << endl;
1262  for(int i=0; i<pinner.numberOfHits(); i++){ pinner.printHitPattern(i,std::cout); }
1263  cout << "expected outer " << nouter << endl;
1264  for(int i=0; i<pouter.numberOfHits(); i++){ pouter.printHitPattern(i,std::cout); }
1265  }
1266 }
1267 
1268 namespace {
1269 
1270  bool recTrackLessZ(const reco::TransientTrack & tk1,
1271  const reco::TransientTrack & tk2)
1272  {
1273  return tk1.stateAtBeamLine().trackStateAtPCA().position().z() < tk2.stateAtBeamLine().trackStateAtPCA().position().z();
1274  }
1275 
1276 }
1277 
1278 
1279 
1280 
1281 void PrimaryVertexAnalyzer4PU::printPVTrks(const Handle<reco::TrackCollection> &recTrks,
1282  const Handle<reco::VertexCollection> &recVtxs,
1283  std::vector<SimPart>& tsim,
1284  std::vector<SimEvent>& simEvt,
1285  bool selectedOnly){
1286  // make a printout of the tracks selected for PV reconstructions, show matching MC tracks, too
1287 
1288  vector<TransientTrack> selTrks;
1289  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1290  t!=recTrks->end(); ++t){
1291  TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
1292  if( (!selectedOnly) || (selectedOnly && theTrackFilter(tt))){ selTrks.push_back(tt); }
1293  }
1294  if (selTrks.size()==0) return;
1295  stable_sort(selTrks.begin(), selTrks.end(), recTrackLessZ);
1296 
1297  // select tracks like for PV reconstruction and match them to sim tracks
1298  reco::TrackCollection selRecTrks;
1299 
1300  for(unsigned int i=0; i<selTrks.size(); i++){ selRecTrks.push_back(selTrks[i].track());}
1301  int* rectosim=supf(tsim, selRecTrks);
1302 
1303 
1304 
1305  // now dump in a format similar to the clusterizer
1306  cout << "printPVTrks" << endl;
1307  cout << "---- z +/- dz 1bet-m ip +/-dip pt phi eta";
1308  if((tsim.size()>0)||(simEvt.size()>0)) {cout << " type pdg zvtx zdca dca zvtx zdca dsz";}
1309  cout << endl;
1310 
1311  cout.precision(4);
1312  int isel=0;
1313  for(unsigned int i=0; i<selTrks.size(); i++){
1314  if (selectedOnly || (theTrackFilter(selTrks[i]))) {
1315  cout << setw (3)<< isel;
1316  isel++;
1317  }else{
1318  cout << " ";
1319  }
1320 
1321 
1322  // is this track in the tracklist of a recvtx ?
1323  int vmatch=-1;
1324  int iv=0;
1325  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1326  v!=recVtxs->end(); ++v){
1327  if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
1328  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
1329  const reco::Track & RTv=*(tv->get());
1330  if(selTrks[i].track().vz()==RTv.vz()) {vmatch=iv;}
1331  }
1332  iv++;
1333  }
1334 
1335  double tz=(selTrks[i].stateAtBeamLine().trackStateAtPCA()).position().z();
1336  double tantheta=tan((selTrks[i].stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1337  double tdz0= selTrks[i].track().dzError();
1338  double tdz2= pow(selTrks[i].track().dzError(),2)+ (pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2))/pow(tantheta,2);
1339 
1340  if(vmatch>-1){
1341  cout << "["<<setw(2)<<vmatch<<"]";
1342  }else{
1343  //int status=theTrackFilter.status(selTrks[i]);
1344  int status=0;
1345  if(status==0){
1346  cout <<" ";
1347  }else{
1348  if(status&0x1){cout << "i";}else{cout << ".";};
1349  if(status&0x2){cout << "c";}else{cout << ".";};
1350  if(status&0x4){cout << "h";}else{cout << ".";};
1351  if(status&0x8){cout << "q";}else{cout << ".";};
1352  }
1353  }
1354  cout << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tdz2);
1355 
1356 
1357  // track quality and hit information, see DataFormats/TrackReco/interface/HitPattern.h
1358  if(selTrks[i].track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";}
1359  if(selTrks[i].track().hitPattern().hasValidHitInFirstPixelBarrel()){cout <<"+";}else{cout << "-";}
1360  cout << setw(1) << selTrks[i].track().hitPattern().pixelBarrelLayersWithMeasurement();
1361  cout << setw(1) << selTrks[i].track().hitPattern().pixelEndcapLayersWithMeasurement();
1362  cout << setw(1) << hex << selTrks[i].track().hitPattern().trackerLayersWithMeasurement()-selTrks[i].track().hitPattern().pixelLayersWithMeasurement()<<dec;
1363  cout << "-" << setw(1)<<hex <<selTrks[i].track().trackerExpectedHitsOuter().numberOfHits() << dec;
1364 
1365 
1366  Measurement1D IP=selTrks[i].stateAtBeamLine().transverseImpactParameter();
1367  cout << setw (8) << IP.value() << "+/-" << setw (6) << IP.error();
1368  if(selTrks[i].track().ptError()<1){
1369  cout << " " << setw(8) << setprecision(2) << selTrks[i].track().pt()*selTrks[i].track().charge();
1370  }else{
1371  cout << " " << setw(7) << setprecision(1) << selTrks[i].track().pt()*selTrks[i].track().charge() << "*";
1372  //cout << "+/-" << setw(6)<< setprecision(2) << selTrks[i].track().ptError();
1373  }
1374  cout << " " << setw(5) << setprecision(2) << selTrks[i].track().phi()
1375  << " " << setw(5) << setprecision(2) << selTrks[i].track().eta() ;
1376 
1377 
1378 
1379  // print MC info, if available
1380  if(simEvt.size()>0){
1381  reco::Track t=selTrks[i].track();
1382  try{
1383  TrackingParticleRef tpr = z2tp_[t.vz()];
1384  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1385  //double vx=parentVertex->position().x(); // problems with tpr->vz()
1386  //double vy=parentVertex->position().y(); work in progress
1387  double vz=parentVertex->position().z();
1388  cout << " " << tpr->eventId().event();
1389  cout << " " << setw(5) << tpr->pdgId();
1390  cout << " " << setw(8) << setprecision(4) << vz;
1391  }catch (...){
1392  cout << " not matched";
1393  }
1394  }else{
1395  // cout << " " << rectosim[i];
1396  if(rectosim[i]>0){
1397  if(tsim[rectosim[i]].type==0){ cout << " prim " ;}else{cout << " sec ";}
1398  cout << " " << setw(5) << tsim[rectosim[i]].pdg;
1399  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zvtx;
1400  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zdcap;
1401  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].ddcap;
1402  double zvtxpull=(tz-tsim[rectosim[i]].zvtx)/sqrt(tdz2);
1403  cout << setw(6) << setprecision(1) << zvtxpull;
1404  double zdcapull=(tz-tsim[rectosim[i]].zdcap)/tdz0;
1405  cout << setw(6) << setprecision(1) << zdcapull;
1406  double dszpull=(selTrks[i].track().dsz()-tsim[rectosim[i]].par[4])/selTrks[i].track().dszError();
1407  cout << setw(6) << setprecision(1) << dszpull;
1408  }
1409  }
1410  cout << endl;
1411  }
1412  delete rectosim;
1413 }
1414 
1415 
1416 void PrimaryVertexAnalyzer4PU::matchRecTracksToVertex(simPrimaryVertex & pv,
1417  const std::vector<SimPart > & tsim,
1418  const edm::Handle<reco::TrackCollection> & recTrks)
1419 {
1420  // find all recTracks that belong to this simulated vertex (not just the ones that are used in
1421  // matching recVertex)
1422 
1423  std::cout << "dump rec tracks: " << std::endl;
1424  int irec=0;
1425  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1426  t!=recTrks->end(); ++t){
1427  reco::TrackBase::ParameterVector p = t->parameters();
1428  std::cout << irec++ << ") " << p << std::endl;
1429  }
1430 
1431  std::cout << "match sim tracks: " << std::endl;
1432  pv.matchedRecTrackIndex.clear();
1433  pv.nMatchedTracks=0;
1434  int isim=0;
1435  for(std::vector<SimPart>::const_iterator s=tsim.begin();
1436  s!=tsim.end(); ++s){
1437  std::cout << isim++ << " " << s->par;// << std::endl;
1438  int imatch=-1;
1439  int irec=0;
1440  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1441  t!=recTrks->end(); ++t){
1442  reco::TrackBase::ParameterVector p = t->parameters();
1443  if (match(s->par,p)){ imatch=irec; }
1444  irec++;
1445  }
1446  pv.matchedRecTrackIndex.push_back(imatch);
1447  if(imatch>-1){
1448  pv.nMatchedTracks++;
1449  std::cout << " matched to rec trk" << imatch << std::endl;
1450  }else{
1451  std::cout << " not matched" << std::endl;
1452  }
1453  }
1454 }
1455 /********************************************************************************************************/
1456 
1457 
1458 
1459 
1460 
1461 
1462 /********************************************************************************************************/
1463 
1464 void PrimaryVertexAnalyzer4PU::getTc(const std::vector<reco::TransientTrack>& tracks,
1465  double & Tc, double & chsq, double & dzmax, double & dztrim, double & m4m2){
1466  if (tracks.size()<2){ Tc=-1; chsq=-1; dzmax=-1; dztrim=-1; m4m2=-1; return;}
1467 
1468  double sumw=0, sumwz=0, sumww=0,sumwwz=0,sumwwzz=0;
1469  double zmin=1e10, zmin1=1e10, zmax1=-1e10, zmax=-1e10;
1470  double m4=0,m3=0,m2=0,m1=0,m0=0;
1471  for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
1472  double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1473  reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
1474  double z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
1475  double dz2= pow((*it).track().dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
1476  // 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);
1477  double w=1./dz2; // take p=1
1478  sumw += w;
1479  sumwz += w*z;
1480  sumwwz += w*w*z;;
1481  sumwwzz += w*w*z*z;
1482  sumww += w*w;
1483  m0 += w;
1484  m1 += w*z;
1485  m2 += w*z*z;
1486  m3 += w*z*z*z;
1487  m4 += w*z*z*z*z;
1488  if(dz2<pow(0.1,2)){
1489  if(z<zmin1){zmin1=z;} if(z<zmin){zmin1=zmin; zmin=z;}
1490  if(z>zmax1){zmax1=z;} if(z>zmax){zmax1=zmax; zmax=z;}
1491  }
1492  }
1493  double z=sumwz/sumw;
1494  double a=sumwwzz-2*z*sumwwz+z*z*sumww;
1495  double b=sumw;
1496  if(tracks.size()>1){
1497  chsq=(m2-m0*z*z)/(tracks.size()-1);
1498  Tc=2.*a/b;
1499  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));
1500  }else{
1501  chsq=0;
1502  Tc=0;
1503  m4m2=0;
1504  }
1505  dzmax=zmax-zmin;
1506  dztrim=zmax1-zmin1;// truncated
1507 }
1508 /********************************************************************************************************/
1509 
1510 
1511 
1512 
1513 
1514 /********************************************************************************************************/
1515 bool PrimaryVertexAnalyzer4PU::truthMatchedTrack( edm::RefToBase<reco::Track> track, TrackingParticleRef & tpr)
1516 
1517 /********************************************************************************************************/
1518 // for a reco track select the matching tracking particle, always use this function to make sure we
1519 // are consistent
1520 // to get the TrackingParticle form the TrackingParticleRef, use ->get();
1521 {
1522  double f=0;
1523  try{
1524  std::vector<std::pair<TrackingParticleRef, double> > tp = r2s_[track];
1525  for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = tp.begin();
1526  it != tp.end(); ++it) {
1527 
1528  if (it->second>f){
1529  tpr=it->first;
1530  f=it->second;
1531  }
1532  }
1533  } catch (Exception event) {
1534  // silly way of testing whether track is in r2s_
1535  }
1536 
1537  // sanity check on track parameters?
1538 
1539  return (f>0.5);
1540 }
1541 /********************************************************************************************************/
1542 
1543 
1544 
1545 
1546 
1547 
1548 /********************************************************************************************************/
1549 std::vector< edm::RefToBase<reco::Track> > PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks(
1550  const reco::Vertex& v
1551  )
1552 // for vertex v get a list of tracks for which truth matching is available
1553 /********************************************************************************************************/
1554 {
1555  std::vector< edm::RefToBase<reco::Track> > b;
1556  TrackingParticleRef tpr;
1557 
1558  for(trackit_t tv=v.tracks_begin(); tv!=v.tracks_end(); tv++){
1559  edm::RefToBase<reco::Track> track=*tv;
1560  if (truthMatchedTrack(track, tpr)){
1561  b.push_back(*tv);
1562  }
1563  }
1564 
1565 
1566  return b;
1567 }
1568 /********************************************************************************************************/
1569 
1570 
1571 
1572 
1573 
1574 /********************************************************************************************************/
1575 std::vector<PrimaryVertexAnalyzer4PU::SimEvent> PrimaryVertexAnalyzer4PU::getSimEvents
1576 (
1577  // const Event& iEvent, const EventSetup& iSetup,
1578  edm::Handle<TrackingParticleCollection> TPCollectionH,
1579  edm::Handle<TrackingVertexCollection> TVCollectionH,
1580  edm::Handle<View<Track> > trackCollectionH
1581  ){
1582 
1583  const TrackingParticleCollection* simTracks = TPCollectionH.product();
1584  const View<Track> tC = *(trackCollectionH.product());
1585 
1586 
1587  vector<SimEvent> simEvt;
1588  map<EncodedEventId, unsigned int> eventIdToEventMap;
1589  map<EncodedEventId, unsigned int>::iterator id;
1590  bool dumpTP=false;
1591  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1592 
1593  if( fabs(it->parentVertex().get()->position().z())>100.) continue; // skip funny entries @ -13900
1594 
1595  unsigned int event=0; //note, this is no longer the same as eventId().event()
1596  id=eventIdToEventMap.find(it->eventId());
1597  if (id==eventIdToEventMap.end()){
1598 
1599  // new event here
1600  SimEvent e;
1601  e.eventId=it->eventId();
1602  event=simEvt.size();
1603  const TrackingVertex *parentVertex= it->parentVertex().get();
1604  if(DEBUG_){
1605  cout << "getSimEvents: ";
1606  cout << it->eventId().bunchCrossing() << "," << it->eventId().event()
1607  << " z="<< it->vz() << " "
1608  << parentVertex->eventId().bunchCrossing() << "," <<parentVertex->eventId().event()
1609  << " z=" << parentVertex->position().z() << endl;
1610  }
1611  if (it->eventId()==parentVertex->eventId()){
1612  //e.x=it->vx(); e.y=it->vy(); e.z=it->vz();// wrong ???
1613  e.x=parentVertex->position().x();
1614  e.y=parentVertex->position().y();
1615  e.z=parentVertex->position().z();
1616  if(e.z<-100){
1617  dumpTP=true;
1618  }
1619  }else{
1620  e.x=0; e.y=0; e.z=-88.;
1621  }
1622  simEvt.push_back(e);
1623  eventIdToEventMap[e.eventId]=event;
1624  }else{
1625  event=id->second;
1626  }
1627 
1628 
1629  simEvt[event].tp.push_back(&(*it));
1630  if( (abs(it->eta())<2.5) && (it->charge()!=0) ){
1631  simEvt[event].sumpt2+=pow(it->pt(),2); // should keep track of decays ?
1632  simEvt[event].sumpt+=it->pt();
1633  }
1634  }
1635 
1636  if(dumpTP){
1637  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1638  std::cout << *it << std::endl;
1639  }
1640  }
1641 
1642 
1643  for(View<Track>::size_type i=0; i<tC.size(); ++i) {
1644  RefToBase<Track> track(trackCollectionH, i);
1645  TrackingParticleRef tpr;
1646  if( truthMatchedTrack(track,tpr)){
1647 
1648  if( eventIdToEventMap.find(tpr->eventId())==eventIdToEventMap.end() ){ cout << "Bug in getSimEvents" << endl; break; }
1649 
1650  z2tp_[track.get()->vz()]=tpr;
1651 
1652  unsigned int event=eventIdToEventMap[tpr->eventId()];
1653  double ipsig=0,ipdist=0;
1654  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1655  double vx=parentVertex->position().x(); // problems with tpr->vz()
1656  double vy=parentVertex->position().y();
1657  double vz=parentVertex->position().z();
1658  double d=sqrt(pow(simEvt[event].x-vx,2)+pow(simEvt[event].y-vy,2)+pow(simEvt[event].z-vz,2))*1.e4;
1659  ipdist=d;
1660  double dxy=track->dxy(vertexBeamSpot_.position());
1661  ipsig=dxy/track->dxyError();
1662 
1663 
1664  TransientTrack t = theB_->build(tC[i]);
1665  t.setBeamSpot(vertexBeamSpot_);
1666  if (theTrackFilter(t)){
1667  simEvt[event].tk.push_back(t);
1668  if(ipdist<5){simEvt[event].tkprim.push_back(t);}
1669  if(ipsig<5){simEvt[event].tkprimsel.push_back(t);}
1670  }
1671 
1672  }
1673  }
1674 
1675 
1676 
1677  AdaptiveVertexFitter theFitter;
1678  cout << "SimEvents " << simEvt.size() << endl;
1679  for(unsigned int i=0; i<simEvt.size(); i++){
1680 
1681  if(simEvt[i].tkprim.size()>0){
1682 
1683  getTc(simEvt[i].tkprimsel, simEvt[i].Tc, simEvt[i].chisq, simEvt[i].dzmax, simEvt[i].dztrim, simEvt[i].m4m2);
1684  TransientVertex v = theFitter.vertex(simEvt[i].tkprim, vertexBeamSpot_);
1685  if (v.isValid()){
1686  simEvt[i].xfit=v.position().x();
1687  simEvt[i].yfit=v.position().y();
1688  simEvt[i].zfit=v.position().z();
1689  // if (simEvt[i].z<-80.){simEvt[i].z=v.position().z();} // for now
1690  }
1691  }
1692 
1693 
1694  if(DEBUG_){
1695  cout << i <<" ) nTP=" << simEvt[i].tp.size()
1696  << " z=" << simEvt[i].z
1697  << " recTrks=" << simEvt[i].tk.size()
1698  << " recTrksPrim=" << simEvt[i].tkprim.size()
1699  << " zfit=" << simEvt[i].zfit
1700  << endl;
1701  }
1702  }
1703 
1704  return simEvt;
1705 }
1706 
1707 
1708 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1709  const Handle<HepMCProduct> evtMC)
1710 {
1711  if(DEBUG_){std::cout << "getSimPVs HepMC " << std::endl;}
1712 
1713  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1714  const HepMC::GenEvent* evt=evtMC->GetEvent();
1715  if (evt) {
1716 // std::cout << "process id " <<evt->signal_process_id()<<std::endl;
1717 // std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
1718 // evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
1719 // std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
1720 
1721 
1722  //int idx=0;
1723  for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
1724  vitr != evt->vertices_end(); ++vitr )
1725  { // loop for vertex ...
1726 
1727  HepMC::FourVector pos = (*vitr)->position();
1728  // if (pos.t()>0) { continue;} // skip secondary vertices, doesn't work for some samples
1729 
1730  if (fabs(pos.z())>1000) continue; // skip funny junk vertices
1731 
1732  bool hasMotherVertex=false;
1733  //std::cout << "mothers" << std::endl;
1734  for ( HepMC::GenVertex::particle_iterator
1735  mother = (*vitr)->particles_begin(HepMC::parents);
1736  mother != (*vitr)->particles_end(HepMC::parents);
1737  ++mother ) {
1738  HepMC::GenVertex * mv=(*mother)->production_vertex();
1739  if (mv) {hasMotherVertex=true;}
1740  //std::cout << "\t"; (*mother)->print();
1741  }
1742  if(hasMotherVertex) {continue;}
1743 
1744 
1745  // could be a new vertex, check all primaries found so far to avoid multiple entries
1746  const double mm=0.1;
1747  simPrimaryVertex sv(pos.x()*mm,pos.y()*mm,pos.z()*mm);
1748  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
1749  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1750  v0!=simpv.end(); v0++){
1751  if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
1752  vp=&(*v0);
1753  break;
1754  }
1755  }
1756  if(!vp){
1757  // this is a new vertex
1758  //std::cout << "this is a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;
1759  simpv.push_back(sv);
1760  vp=&simpv.back();
1761 // }else{
1762 // std::cout << "this is not a new vertex" << std::endl;
1763  }
1764 
1765 
1766  // store the gen vertex barcode with this simpv
1767  vp->genVertex.push_back((*vitr)->barcode());
1768 
1769 
1770  // collect final state descendants and sum up momenta etc
1771  for ( HepMC::GenVertex::particle_iterator
1772  daughter = (*vitr)->particles_begin(HepMC::descendants);
1773  daughter != (*vitr)->particles_end(HepMC::descendants);
1774  ++daughter ) {
1775  //std::cout << "checking daughter type " << (*daughter)->pdg_id() << " final :" <<isFinalstateParticle(*daughter) << std::endl;
1776  if (isFinalstateParticle(*daughter)){
1777  if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
1778  == vp->finalstateParticles.end()){
1779  vp->finalstateParticles.push_back((*daughter)->barcode());
1780  HepMC::FourVector m=(*daughter)->momentum();
1781  //std::cout << "adding particle to primary " << m.px()<< " " << m.py() << " " << m.pz() << std::endl;
1782  vp->ptot.setPx(vp->ptot.px()+m.px());
1783  vp->ptot.setPy(vp->ptot.py()+m.py());
1784  vp->ptot.setPz(vp->ptot.pz()+m.pz());
1785  vp->ptot.setE(vp->ptot.e()+m.e());
1786  vp->ptsq+=(m.perp())*(m.perp());
1787  // count relevant particles
1788  if ( (m.perp()>0.2) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
1789  vp->nGenTrk++;
1790  }
1791 
1792  hsimPV["rapidity"]->Fill(m.pseudoRapidity());
1793  if( (m.perp()>0.8) && isCharged( *daughter ) ){
1794  hsimPV["chRapidity"]->Fill(m.pseudoRapidity());
1795  }
1796  hsimPV["pt"]->Fill(m.perp());
1797  }//new final state particle for this vertex
1798  }//daughter is a final state particle
1799  }
1800 
1801 
1802  //idx++;
1803  }
1804  }
1805  if(verbose_){
1806  cout << "------- PrimaryVertexAnalyzer4PU simPVs -------" << endl;
1807  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1808  v0!=simpv.end(); v0++){
1809  cout << "z=" << v0->z
1810  << " px=" << v0->ptot.px()
1811  << " py=" << v0->ptot.py()
1812  << " pz=" << v0->ptot.pz()
1813  << " pt2="<< v0->ptsq
1814  << endl;
1815  }
1816  cout << "-----------------------------------------------" << endl;
1817  }
1818  return simpv;
1819 }
1820 
1821 
1822 
1823 
1824 
1825 
1826 
1827 
1828 /* get sim pv from TrackingParticles/TrackingVertex */
1829 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1830  const edm::Handle<TrackingVertexCollection> tVC
1831  )
1832 {
1833  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1834  //std::cout <<"number of vertices " << tVC->size() << std::endl;
1835 
1836  if(DEBUG_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;}
1837 
1838  for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) {
1839 
1840  if(DEBUG_){
1841  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;
1842  for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){
1843  std::cout << *gv << std::endl;
1844  }
1845  std::cout << "----------" << std::endl;
1846  }
1847 
1848  // bool hasMotherVertex=false;
1849  if ((unsigned int) v->eventId().event()<simpv.size()) continue;
1850  //if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work)
1851  if (fabs(v->position().z())>1000) continue; // skip funny junk vertices
1852 
1853  // could be a new vertex, check all primaries found so far to avoid multiple entries
1854  const double mm=1.0; // for tracking vertices
1855  simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm);
1856 
1857  //cout << "sv: " << (v->eventId()).event() << endl;
1858  sv.eventId=v->eventId();
1859 
1860  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){
1861  //cout <<((**iTrack).eventId()).event() << " "; // an iterator of Refs, dereference twice. Cool eyh?
1862  sv.eventId=(**iTrack).eventId();
1863  }
1864  //cout <<endl;
1865  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
1866  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1867  v0!=simpv.end(); v0++){
1868  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)){
1869  vp=&(*v0);
1870  break;
1871  }
1872  }
1873  if(!vp){
1874  // this is a new vertex
1875  if(DEBUG_){std::cout << "this is a new vertex " << sv.eventId.event() << " " << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
1876  simpv.push_back(sv);
1877  vp=&simpv.back();
1878  }else{
1879  if(DEBUG_){std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
1880  }
1881 
1882 
1883  // Loop over daughter track(s)
1884  if(DEBUG_){
1885  for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) {
1886  std::cout << " Daughter momentum: " << (*(*iTP)).momentum();
1887  std::cout << " Daughter type " << (*(*iTP)).pdgId();
1888  std::cout << std::endl;
1889  }
1890  }
1891  }
1892  if(DEBUG_){
1893  cout << "------- PrimaryVertexAnalyzer4PU simPVs from TrackingVertices -------" << endl;
1894  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1895  v0!=simpv.end(); v0++){
1896  cout << "z=" << v0->z << " event=" << v0->eventId.event() << endl;
1897  }
1898  cout << "-----------------------------------------------" << endl;
1899  }
1900  return simpv;
1901 }
1902 
1903 
1904 
1905 
1906 
1907 
1908 // ------------ method called to produce the data ------------
1909 void
1910 PrimaryVertexAnalyzer4PU::analyze(const Event& iEvent, const EventSetup& iSetup)
1911 {
1912 
1913  bool MC=false;
1914  std::vector<simPrimaryVertex> simpv; // a list of primary MC vertices
1915  std::vector<SimPart> tsim;
1916  std::string mcproduct="generator"; // starting with 3_1_0 pre something
1917 
1918  eventcounter_++;
1919  run_ = iEvent.id().run();
1920  luminosityBlock_ = iEvent.luminosityBlock();
1921  event_ = iEvent.id().event();
1922  bunchCrossing_ = iEvent.bunchCrossing();
1923  orbitNumber_ = iEvent.orbitNumber();
1924 
1925  dumpThisEvent_ = false;
1926 
1927 
1928 
1929  if(verbose_){
1930  std::cout << endl
1931  << "PrimaryVertexAnalyzer4PU::analyze event counter=" << eventcounter_
1932  << " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
1933  << " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
1934  //<< " selected = " << good
1935  << std::endl;
1936  }
1937 
1938 
1939  try{
1940  iSetup.getData(pdt_);
1941  }catch(const Exception&){
1942  std::cout << "Some problem occurred with the particle data table. This may not work !" <<std::endl;
1943  }
1944 
1945  Handle<reco::VertexCollection> recVtxs;
1946  bool bnoBS=iEvent.getByLabel("offlinePrimaryVertices", recVtxs);
1947 
1948  Handle<reco::VertexCollection> recVtxsBS;
1949  bool bBS=iEvent.getByLabel("offlinePrimaryVerticesWithBS", recVtxsBS);
1950 
1951  Handle<reco::VertexCollection> recVtxsDA;
1952  bool bDA=iEvent.getByLabel("offlinePrimaryVerticesDA", recVtxsDA);
1953 
1954 // Handle<reco::VertexCollection> recVtxsPIX;
1955 // bool bPIX=iEvent.getByLabel("pixelVertices", recVtxsPIX);
1956 // bPIX=false;
1957 
1958 // Handle<reco::VertexCollection> recVtxsMVF;
1959 // bool bMVF=iEvent.getByLabel("offlinePrimaryVerticesMVF", recVtxsMVF);
1960 
1961  Handle<reco::TrackCollection> recTrks;
1962  iEvent.getByLabel(recoTrackProducer_, recTrks);
1963 
1964 
1965  int nhighpurity=0, ntot=0;
1966  for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
1967  ntot++;
1968  if(t->quality(reco::TrackBase::highPurity)) nhighpurity++;
1969  }
1970  if(ntot>10) hnoBS["highpurityTrackFraction"]->Fill(float(nhighpurity)/float(ntot));
1971  if((recoTrackProducer_=="generalTracks")&&(nhighpurity<0.25*ntot)){
1972  if(verbose_){ cout << "rejected, " << nhighpurity << " highpurity out of " << ntot << " total tracks "<< endl<< endl;}
1973  return;
1974  }
1975 
1976 
1977 
1978 
1979 
1980  if(iEvent.getByType(recoBeamSpotHandle_)){
1981  vertexBeamSpot_= *recoBeamSpotHandle_;
1982  wxy2_=pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2);
1983  Fill(hsimPV, "xbeam",vertexBeamSpot_.x0()); Fill(hsimPV, "wxbeam",vertexBeamSpot_.BeamWidthX());
1984  Fill(hsimPV, "ybeam",vertexBeamSpot_.y0()); Fill(hsimPV, "wybeam",vertexBeamSpot_.BeamWidthY());
1985  Fill(hsimPV, "zbeam",vertexBeamSpot_.z0());// Fill("wzbeam",vertexBeamSpot_.BeamWidthZ());
1986  }else{
1987  cout << " beamspot not found, using dummy " << endl;
1988  vertexBeamSpot_=reco::BeamSpot();// dummy
1989  }
1990 
1991 
1992  if(bnoBS){
1993  if((recVtxs->begin()->isValid())&&(recVtxs->begin()->ndof()>1)&&(recVtxs->begin()->ndof()>(0.0*recVtxs->begin()->tracksSize()))){ // was 5 and 0.2
1994  Fill(hnoBS,"xrecBeamvsdxXBS",recVtxs->begin()->xError(),recVtxs->begin()->x()-vertexBeamSpot_.x0());
1995  Fill(hnoBS,"yrecBeamvsdyXBS",recVtxs->begin()->yError(),recVtxs->begin()->y()-vertexBeamSpot_.y0());
1996 
1997  if(printXBS_) {
1998  cout << Form("XBS %10d %5d %10d %4d %lu %5.2f %8f %8f %8f %8f %8f %8f",
1999  run_,luminosityBlock_,event_,bunchCrossing_,
2000  (unsigned long)(recVtxs->begin()->tracksSize()), recVtxs->begin()->ndof(),
2001  recVtxs->begin()->x(), recVtxs->begin()->xError(),
2002  recVtxs->begin()->y(), recVtxs->begin()->yError(),
2003  recVtxs->begin()->z(), recVtxs->begin()->zError()
2004  ) << endl;
2005  }
2006 
2007  }
2008  }
2009 
2010 
2011  // for the associator
2012  Handle<View<Track> > trackCollectionH;
2013  iEvent.getByLabel(recoTrackProducer_,trackCollectionH);
2014 
2015  Handle<HepMCProduct> evtMC;
2016 
2017  Handle<SimVertexContainer> simVtxs;
2018  iEvent.getByLabel( simG4_, simVtxs);
2019 
2020  Handle<SimTrackContainer> simTrks;
2021  iEvent.getByLabel( simG4_, simTrks);
2022 
2023 
2024 
2025 
2026 
2027  edm::Handle<TrackingParticleCollection> TPCollectionH ;
2028  edm::Handle<TrackingVertexCollection> TVCollectionH ;
2029  bool gotTP=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionH);
2030  bool gotTV=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TVCollectionH);
2031 
2032 
2033  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
2034  fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
2035 
2036 
2037  vector<SimEvent> simEvt;
2038  if (gotTP && gotTV ){
2039 
2040  edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
2041  iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
2042  associatorByHits_ = (TrackAssociatorBase *) theHitsAssociator.product();
2043  r2s_ = associatorByHits_->associateRecoToSim (trackCollectionH,TPCollectionH, &iEvent );
2044  //simEvt=getSimEvents(iEvent, TPCollectionH, TVCollectionH, trackCollectionH);
2045  simEvt=getSimEvents(TPCollectionH, TVCollectionH, trackCollectionH);
2046 
2047  if (simEvt.size()==0){
2048  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2049  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2050  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2051  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2052  cout << " !!!!!!!!!!!!!!!!!!!!!! got TrackingParticles but no simEvt !!!!!!!!!!!!!!!!!" << endl;
2053  cout << " dumping Tracking particles " << endl;
2054  const TrackingParticleCollection* simTracks = TPCollectionH.product();
2055  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
2056  cout << *it << endl;
2057  }
2058  cout << " dumping Tracking Vertices " << endl;
2059  const TrackingVertexCollection* tpVtxs = TVCollectionH.product();
2060  for(TrackingVertexCollection::const_iterator iv=tpVtxs->begin(); iv!=tpVtxs->end(); iv++){
2061  cout << *iv << endl;
2062  }
2063  if(iEvent.getByLabel(mcproduct,evtMC)){
2064  cout << "dumping simTracks" << endl;
2065  for(SimTrackContainer::const_iterator t=simTrks->begin(); t!=simTrks->end(); ++t){
2066  cout << *t << endl;
2067  }
2068  cout << "dumping simVertexes" << endl;
2069  for(SimVertexContainer::const_iterator vv=simVtxs->begin();
2070  vv!=simVtxs->end();
2071  ++vv){
2072  cout << *vv << endl;
2073  }
2074  }else{
2075  cout << "no hepMC" << endl;
2076  }
2077  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2078 
2079  const HepMC::GenEvent* evt=evtMC->GetEvent();
2080  if(evt){
2081  std::cout << "process id " <<evt->signal_process_id()<<std::endl;
2082  std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
2083  evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
2084  std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
2085  evt->print();
2086  }else{
2087  cout << "no event in HepMC" << endl;
2088  }
2089  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2090 
2091  }
2092  }
2093 
2094 
2095 
2096 
2097  if(gotTV){
2098 
2099  MC=true;
2100  if(verbose_){
2101  cout << "Found Tracking Vertices " << endl;
2102  }
2103  simpv=getSimPVs(TVCollectionH);
2104 
2105 
2106  }else if(iEvent.getByLabel(mcproduct,evtMC)){
2107 
2108  MC=true;
2109  simpv=getSimPVs(evtMC);
2110 
2111  if(verbose_){
2112  cout << "Using HepMCProduct " << endl;
2113  std::cout << "simtrks " << simTrks->size() << std::endl;
2114  }
2115  tsim = PrimaryVertexAnalyzer4PU::getSimTrkParameters(simTrks, simVtxs, simUnit_);
2116 
2117  }else{
2118  MC=false;
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  double zmatch=-99; 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; zmatch=ev->z;
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.,isnan(v->position().x())*1.);
3376  Fill(h,"nans",2.,isnan(v->position().y())*1.);
3377  Fill(h,"nans",3.,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., isnan(v->covariance(i, j))*1.);
3384  if (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
ExpressReco_HICollisions_FallBack.e
tuple e
Definition: ExpressReco_HICollisions_FallBack.py:1068
reco::TrackBase::qoverp
double qoverp() const
q/p
Definition: TrackBase.h:114
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:128
type
type
Definition: HCALResponse.h:22
PrimaryVertexAnalyzer4PU::theTrackFilter
TrackFilterForPVFinding theTrackFilter
Definition: PrimaryVertexAnalyzer4PU.h:359
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:44
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:69
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
reco::TrackBase::d0Error
double d0Error() const
error on d0
Definition: TrackBase.h:204
ExpressReco_HICollisions_FallBack.track
tuple track
Definition: ExpressReco_HICollisions_FallBack.py:4641
PrimaryVertexAnalyzer4PU::matchVertex
bool matchVertex(const simPrimaryVertex &vsim, const reco::Vertex &vrec)
Definition: PrimaryVertexAnalyzer4PU.cc:972
PrimaryVertexAnalyzer4PU::getSimTrkParameters
std::vector< SimPart > getSimTrkParameters(edm::Handle< edm::SimTrackContainer > &simTrks, edm::Handle< edm::SimVertexContainer > &simVtcs, double simUnit=1.0)
Definition: PrimaryVertexAnalyzer4PU.cc:731
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:1576
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:1001
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:45
debug_cff.d0
tuple d0
Definition: debug_cff.py:3
reco::TransientTrack::setBeamSpot
void setBeamSpot(const reco::BeamSpot &beamSpot)
Definition: TransientTrack.h:43
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:110
PrimaryVertexAnalyzer4PU::ParameterVector
reco::TrackBase::ParameterVector ParameterVector
Definition: PrimaryVertexAnalyzer4PU.h:79
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
PrimaryVertexAnalyzer4PU::hsimPV
std::map< std::string, TH1 * > hsimPV
Definition: PrimaryVertexAnalyzer4PU.h:353
TrackingVertex::g4v_iterator
std::vector< SimVertex >::const_iterator g4v_iterator
Definition: TrackingVertex.h:34
PrimaryVertexAnalyzer4PU::recoBeamSpotHandle_
edm::Handle< reco::BeamSpot > recoBeamSpotHandle_
Definition: PrimaryVertexAnalyzer4PU.h:362
PrimaryVertexAnalyzer4PU::SimEvent::tk
std::vector< reco::TransientTrack > tk
Definition: PrimaryVertexAnalyzer4PU.h:148
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:116
reco::TrackBase::dxyError
double dxyError() const
error on dxy
Definition: TrackBase.h:202
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
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
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:960
PrimaryVertexAnalyzer4PU::theB_
edm::ESHandle< TransientTrackBuilder > theB_
Definition: PrimaryVertexAnalyzer4PU.h:363
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:138
NULL
#define NULL
Definition: scimark2.h:8
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
theta
Geom::Theta< T > theta() const
Definition: newBasic3DVector.h:115
edm::Event::getByType
bool getByType(Handle< PROD > &result) const
Definition: Event.h:403
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
reco::TrackBase::ParameterVector
math::Vector< dimension >::type ParameterVector
parameter vector
Definition: TrackBase.h:69
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:132
reco::HitPattern::pixelLayersWithMeasurement
int pixelLayersWithMeasurement() const
Definition: HitPattern.cc:865
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:1122
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:860
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:93
reco::TransientTrack::stateAtBeamLine
TrajectoryStateClosestToBeamLine stateAtBeamLine() const
Definition: TransientTrack.h:88
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:67
PrimaryVertexAnalyzer4PU::hnoBS
std::map< std::string, TH1 * > hnoBS
Definition: PrimaryVertexAnalyzer4PU.h:349
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
reco::TrackBase::algo
TrackAlgorithm algo() const
Definition: TrackBase.h:303
PrimaryVertexAnalyzer4PU::Fill
void Fill(std::map< std::string, TH1 * > &h, std::string s, double x)
Definition: PrimaryVertexAnalyzer4PU.h:186
PrimaryVertexAnalyzer4PU::hDA
std::map< std::string, TH1 * > hDA
Definition: PrimaryVertexAnalyzer4PU.h:350
PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP
void analyzeVertexCollectionTP(std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< SimEvent > &simEvt, const std::string message="")
Definition: PrimaryVertexAnalyzer4PU.cc:2449
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:248
PrimaryVertexAnalyzer4PU::isFinalstateParticle
bool isFinalstateParticle(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:983
PrimaryVertexAnalyzer4PU::isCharged
bool isCharged(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:988
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:220
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:140
z
Definition: DDAxes.h:10
PrimaryVertexAnalyzer4PU::truthMatchedTrack
bool truthMatchedTrack(edm::RefToBase< reco::Track >, TrackingParticleRef &)
Definition: PrimaryVertexAnalyzer4PU.cc:1515
ExpressReco_HICollisions_FallBack.y
tuple y
Definition: ExpressReco_HICollisions_FallBack.py:4645
dt_offlineAnalysis_common_cff.reco
tuple reco
Definition: dt_offlineAnalysis_common_cff.py:54
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.
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
EncodedEventId::bunchCrossing
int bunchCrossing() const
get the detector field from this detid
Definition: EncodedEventId.h:24
TransientVertex::position
GlobalPoint position() const
Definition: TransientVertex.h:113
PrimaryVertexAnalyzer4PU::bookVertexHistograms
std::map< std::string, TH1 * > bookVertexHistograms()
Definition: PrimaryVertexAnalyzer4PU.cc:119
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:130
lt
< trclass="colgroup">< tdclass="colgroup"colspan=5 > Muon Digi collections</td >< tr >< td >< ahref="classDTDigi.html"> DTDigi</a ></td >< td >< ahref="DataFormats_DTDigi.html"> MuonDigiCollection & lt
Definition: MuonDigiCollections.doi:5
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:1281
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
ExpressReco_HICollisions_FallBack.beamSpot
tuple beamSpot
Definition: ExpressReco_HICollisions_FallBack.py:592
reco::TrackBase::qualityMask
int qualityMask() const
Definition: TrackBase.h:269
AdaptiveVertexFitter::vertex
virtual CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &) const
Definition: AdaptiveVertexFitter.cc:155
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
edm::View::size_type
unsigned int size_type
Definition: View.h:90
reco::HitPattern::numberOfHits
int numberOfHits() const
Definition: HitPattern.cc:312
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:118
PrimaryVertexAnalyzer4PU::getSimPVs
std::vector< simPrimaryVertex > getSimPVs(const edm::Handle< edm::HepMCProduct > evtMC)
Definition: PrimaryVertexAnalyzer4PU.cc:1708
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:77
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:218
PrimaryVertexAnalyzer4PU::printSimVtxs
void printSimVtxs(const edm::Handle< edm::SimVertexContainer > simVtxs)
Definition: PrimaryVertexAnalyzer4PU.cc:1151
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
ExpressReco_HICollisions_FallBack.id
tuple id
Definition: ExpressReco_HICollisions_FallBack.py:9307
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:216
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
ExpressReco_HICollisions_FallBack.pt
tuple pt
Definition: ExpressReco_HICollisions_FallBack.py:1070
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:10
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:89
TrackAssociatorBase::associateRecoToSim
virtual reco::RecoToSimCollection associateRecoToSim(edm::Handle< edm::View< reco::Track > > &tCH, edm::Handle< TrackingParticleCollection > &tPCH, const edm::Event *event=0, const edm::EventSetup *setup=0) const
compare reco to sim the handle of reco::Track and TrackingParticle collections
Definition: TrackAssociatorBase.h:44
nt
int nt
Definition: AMPTWrapper.h:32
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:136
PrimaryVertexAnalyzer4PU::hBS
std::map< std::string, TH1 * > hBS
Definition: PrimaryVertexAnalyzer4PU.h:348
PrimaryVertexAnalyzer4PU::SimEvent::tkprim
std::vector< reco::TransientTrack > tkprim
Definition: PrimaryVertexAnalyzer4PU.h:149
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:208
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:146
PrimaryVertexAnalyzer4PU::associatorByHits_
TrackAssociatorBase * associatorByHits_
Definition: PrimaryVertexAnalyzer4PU.h:355
DetId
Definition: DetId.h:20
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:78
PrimaryVertexAnalyzer4PU::matchRecTracksToVertex
void matchRecTracksToVertex(simPrimaryVertex &pv, const std::vector< SimPart > &tsim, const edm::Handle< reco::TrackCollection > &recTrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1416
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:977
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
reco::TransientTrack::track
const Track & track() const
Definition: TransientTrack.h:84
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
PrimaryVertexAnalyzer4PU::printRecTrks
void printRecTrks(const edm::Handle< reco::TrackCollection > &recTrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1195
PrimaryVertexAnalyzer4PU::getTc
void getTc(const std::vector< reco::TransientTrack > &, double &, double &, double &, double &, double &)
Definition: PrimaryVertexAnalyzer4PU.cc:1464
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
part
part
Definition: HCALResponse.h:21
TrackingVertexCollection
std::vector< TrackingVertex > TrackingVertexCollection
Definition: TrackingVertexContainer.h:8
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
PrimaryVertexAnalyzer4PU::r2s_
reco::RecoToSimCollection r2s_
Definition: PrimaryVertexAnalyzer4PU.h:356
reco::HitPattern::pixelBarrelLayersWithMeasurement
int pixelBarrelLayersWithMeasurement() const
Definition: HitPattern.cc:877
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:62
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:79
edm::View::size
size_type size() const
b
double b
Definition: hdecay.h:120
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
argparse.message
string message
Definition: argparse.py:126
Measurement1D::value
double value() const
Definition: Measurement1D.h:28
PrimaryVertexAnalyzer4PU::supf
int * supf(std::vector< SimPart > &simtrks, const reco::TrackCollection &trks)
Definition: PrimaryVertexAnalyzer4PU.cc:829
edm::Handle::product
T const * product() const
Definition: Handle.h:74
PrimaryVertexAnalyzer4PU::PrimaryVertexAnalyzer4PU
PrimaryVertexAnalyzer4PU(const edm::ParameterSet &)
Definition: PrimaryVertexAnalyzer4PU.cc:74
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
TrackingVertex::eventId
const EncodedEventId & eventId() const
Definition: TrackingVertex.h:95
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:23
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:144
ExpressReco_HICollisions_FallBack.x
tuple x
Definition: ExpressReco_HICollisions_FallBack.py:4646
PrimaryVertexAnalyzer4PU::printSimTrks
void printSimTrks(const edm::Handle< edm::SimTrackContainer > simVtrks)
Definition: PrimaryVertexAnalyzer4PU.cc:1174
begin
#define begin
Definition: vmac.h:31
PrimaryVertexAnalyzer4PU::add
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
Definition: PrimaryVertexAnalyzer4PU.h:184
a
double a
Definition: hdecay.h:121
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:41
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:112
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:1144
asciidump.s
string s
Definition: asciidump.py:422
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: TrackingParticleFwd.h:8
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:1100
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:239
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
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:120
PV3DBase::x
T x() const
Definition: PV3DBase.h:56
PrimaryVertexAnalyzer4PU::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexAnalyzer4PU.cc:1910
TransientVertex::isValid
bool isValid() const
Definition: TransientVertex.h:133
PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks
std::vector< edm::RefToBase< reco::Track > > getTruthMatchedVertexTracks(const reco::Vertex &)
Definition: PrimaryVertexAnalyzer4PU.cc:1549
edm::RefToBase::get
value_type const * get() const
Definition: RefToBase.h:207
PrimaryVertexAnalyzer4PU::pdt_
edm::ESHandle< ParticleDataTable > pdt_
Definition: PrimaryVertexAnalyzer4PU.h:329
PrimaryVertexAnalyzer4PU::z2tp_
std::map< double, TrackingParticleRef > z2tp_
Definition: PrimaryVertexAnalyzer4PU.h:357
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:322
v
mathSSE::Vec4< T > v
Definition: newBasic3DVector.h:210
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
ExpressReco_HICollisions_FallBack.beamspot
tuple beamspot
Definition: ExpressReco_HICollisions_FallBack.py:1673
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:134
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:142
reco::TrackBase::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:71
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:27
Gflash::par
const double par[8 *NPar][4]
Definition: GflashNameSpace.h:177
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