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