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.hitPattern().numberOfHits(HitPattern::MISSING_INNER_HITS));
1118  Fill(h, "expectedOuter_" + ttype, t.hitPattern().numberOfHits(HitPattern::MISSING_OUTER_HITS));
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 
1271  cout << "subdet layers valid lost" << endl;
1272  cout << Form(" barrel %2d %2d %2d",
1273  p.pixelBarrelLayersWithMeasurement(),
1274  p.numberOfValidPixelBarrelHits(),
1275  p.numberOfLostPixelBarrelHits(HitPattern::TRACK_HITS))
1276  << endl;
1277 
1278  cout << Form(" fwd %2d %2d %2d",
1279  p.pixelEndcapLayersWithMeasurement(),
1280  p.numberOfValidPixelEndcapHits(),
1281  p.numberOfLostPixelEndcapHits(HitPattern::TRACK_HITS))
1282  << endl;
1283  cout << Form(" pixel %2d %2d %2d",
1284  p.pixelLayersWithMeasurement(),
1285  p.numberOfValidPixelHits(),
1286  p.numberOfLostPixelHits(HitPattern::TRACK_HITS))
1287  << endl;
1288  cout << Form(" tracker %2d %2d %2d",
1289  p.trackerLayersWithMeasurement(),
1290  p.numberOfValidTrackerHits(),
1291  p.numberOfLostTrackerHits(HitPattern::TRACK_HITS))
1292  << endl;
1293  cout << endl;
1294 
1295  // double d0Error=t.d0Error();
1296  // double d0=t.dxy(myBeamSpot);
1297 
1298  //
1299  for(trackingRecHit_iterator hit=t->recHitsBegin(); hit!=t->recHitsEnd(); hit++){
1300  if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
1301  bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1302  bool endcap = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1303  if (barrel){
1304  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1305  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1306  if (clust.isNonnull()) {
1307  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;
1308  }
1309  }else if(endcap){
1310  const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
1311  edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
1312  if (clust.isNonnull()) {
1313  cout << Form(" endcap cluster size=%2d charge=%6.1f wx=%2d wy=%2d",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY()) << endl;
1314  }
1315  }
1316  }
1317  }
1318 
1319  cout << "hitpattern" << endl;
1320  for(int i = 0; i < p.numberOfHits(HitPattern::TRACK_HITS); i++){
1321  p.printHitPattern(HitPattern::TRACK_HITS, i, std::cout);
1322  }
1323 
1324  cout << "expected inner " << p.numberOfHits(HitPattern::MISSING_INNER_HITS) << endl;
1325  for(int i = 0; i < p.numberOfHits(HitPattern::MISSING_INNER_HITS); i++){
1326  p.printHitPattern(HitPattern::MISSING_INNER_HITS, i, std::cout);
1327  }
1328 
1329  cout << "expected outer " << p.numberOfHits(HitPattern::MISSING_OUTER_HITS) << endl;
1330  for(int i = 0; i < p.numberOfHits(HitPattern::MISSING_OUTER_HITS); i++){
1331  p.printHitPattern(HitPattern::MISSING_OUTER_HITS, i, std::cout);
1332  }
1333 
1334  }
1335 }
1336 
1337 namespace {
1338 
1339  bool recTrackLessZ(const reco::TransientTrack & tk1,
1340  const reco::TransientTrack & tk2)
1341  {
1342  return tk1.stateAtBeamLine().trackStateAtPCA().position().z() < tk2.stateAtBeamLine().trackStateAtPCA().position().z();
1343  }
1344 
1345 }
1346 
1347 
1348 
1349 void PrimaryVertexAnalyzer4PU::printPVTrks(const Handle<reco::TrackCollection> &recTrks,
1350  const Handle<reco::VertexCollection> &recVtxs,
1351  std::vector<SimPart>& tsim,
1352  std::vector<SimEvent>& simEvt,
1353  bool selectedOnly){
1354  // make a printout of the tracks selected for PV reconstructions, show matching MC tracks, too
1355 
1356  vector<TransientTrack> selTrks;
1357  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1358  t!=recTrks->end(); ++t){
1359  TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
1360  if( (!selectedOnly) || (selectedOnly && theTrackFilter(tt))){ selTrks.push_back(tt); }
1361  }
1362  if (selTrks.size()==0) return;
1363  stable_sort(selTrks.begin(), selTrks.end(), recTrackLessZ);
1364 
1365  // select tracks like for PV reconstruction and match them to sim tracks
1366  reco::TrackCollection selRecTrks;
1367 
1368  for(unsigned int i=0; i<selTrks.size(); i++){ selRecTrks.push_back(selTrks[i].track());}
1369  std::vector<int> rectosim=supf(tsim, selRecTrks);
1370 
1371 
1372 
1373  // now dump in a format similar to the clusterizer
1374  cout << "printPVTrks" << endl;
1375  cout << "---- z +/- dz 1bet-m ip +/-dip pt phi eta";
1376  if((tsim.size()>0)||(simEvt.size()>0)) {cout << " type pdg zvtx zdca dca zvtx zdca dsz";}
1377  cout << endl;
1378 
1379  cout.precision(4);
1380  int isel=0;
1381  for(unsigned int i=0; i<selTrks.size(); i++){
1382  if (selectedOnly || (theTrackFilter(selTrks[i]))) {
1383  cout << setw (3)<< isel;
1384  isel++;
1385  }else{
1386  cout << " ";
1387  }
1388 
1389 
1390  // is this track in the tracklist of a recvtx ?
1391  int vmatch=-1;
1392  int iv=0;
1393  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1394  v!=recVtxs->end(); ++v){
1395  if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
1396  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
1397  const reco::Track & RTv=*(tv->get());
1398  if(selTrks[i].track().vz()==RTv.vz()) {vmatch=iv;}
1399  }
1400  iv++;
1401  }
1402 
1403  double tz=(selTrks[i].stateAtBeamLine().trackStateAtPCA()).position().z();
1404  double tantheta=tan((selTrks[i].stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1405  double tdz0= selTrks[i].track().dzError();
1406  double tdz2= pow(selTrks[i].track().dzError(),2)+ (pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2))/pow(tantheta,2);
1407 
1408  if(vmatch>-1){
1409  cout << "["<<setw(2)<<vmatch<<"]";
1410  }else{
1411  //int status=theTrackFilter.status(selTrks[i]);
1412  int status=0;
1413  if(status==0){
1414  cout <<" ";
1415  }else{
1416  if(status&0x1){cout << "i";}else{cout << ".";};
1417  if(status&0x2){cout << "c";}else{cout << ".";};
1418  if(status&0x4){cout << "h";}else{cout << ".";};
1419  if(status&0x8){cout << "q";}else{cout << ".";};
1420  }
1421  }
1422  cout << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tdz2);
1423 
1424 
1425  // track quality and hit information, see DataFormats/TrackReco/interface/HitPattern.h
1426  if(selTrks[i].track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";}
1427  if(selTrks[i].track().hitPattern().hasValidHitInFirstPixelBarrel()){ cout << "+"; } else { cout << "-"; }
1428  cout << setw(1) << selTrks[i].track().hitPattern().pixelBarrelLayersWithMeasurement();
1429  cout << setw(1) << selTrks[i].track().hitPattern().pixelEndcapLayersWithMeasurement();
1430  cout << setw(1) << hex << selTrks[i].track().hitPattern().trackerLayersWithMeasurement()
1431  - selTrks[i].track().hitPattern().pixelLayersWithMeasurement() << dec;
1432  cout << "-" << setw(1) << hex << selTrks[i].track().hitPattern().numberOfHits(HitPattern::MISSING_OUTER_HITS) << dec;
1433 
1434 
1435  Measurement1D IP=selTrks[i].stateAtBeamLine().transverseImpactParameter();
1436  cout << setw (8) << IP.value() << "+/-" << setw (6) << IP.error();
1437  if(selTrks[i].track().ptError()<1){
1438  cout << " " << setw(8) << setprecision(2) << selTrks[i].track().pt()*selTrks[i].track().charge();
1439  }else{
1440  cout << " " << setw(7) << setprecision(1) << selTrks[i].track().pt()*selTrks[i].track().charge() << "*";
1441  //cout << "+/-" << setw(6)<< setprecision(2) << selTrks[i].track().ptError();
1442  }
1443  cout << " " << setw(5) << setprecision(2) << selTrks[i].track().phi()
1444  << " " << setw(5) << setprecision(2) << selTrks[i].track().eta() ;
1445 
1446 
1447 
1448  // print MC info, if available
1449  if(simEvt.size()>0){
1450  reco::Track t=selTrks[i].track();
1451  try{
1452  TrackingParticleRef tpr = z2tp_[t.vz()];
1453  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1454  //double vx=parentVertex->position().x(); // problems with tpr->vz()
1455  //double vy=parentVertex->position().y(); work in progress
1456  double vz=parentVertex->position().z();
1457  cout << " " << tpr->eventId().event();
1458  cout << " " << setw(5) << tpr->pdgId();
1459  cout << " " << setw(8) << setprecision(4) << vz;
1460  }catch (...){
1461  cout << " not matched";
1462  }
1463  }else{
1464  // cout << " " << rectosim[i];
1465  if(rectosim[i]>0){
1466  if(tsim[rectosim[i]].type==0){ cout << " prim " ;}else{cout << " sec ";}
1467  cout << " " << setw(5) << tsim[rectosim[i]].pdg;
1468  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zvtx;
1469  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zdcap;
1470  cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].ddcap;
1471  double zvtxpull=(tz-tsim[rectosim[i]].zvtx)/sqrt(tdz2);
1472  cout << setw(6) << setprecision(1) << zvtxpull;
1473  double zdcapull=(tz-tsim[rectosim[i]].zdcap)/tdz0;
1474  cout << setw(6) << setprecision(1) << zdcapull;
1475  double dszpull=(selTrks[i].track().dsz()-tsim[rectosim[i]].par[4])/selTrks[i].track().dszError();
1476  cout << setw(6) << setprecision(1) << dszpull;
1477  }
1478  }
1479  cout << endl;
1480  }
1481 }
1482 
1483 
1484 void PrimaryVertexAnalyzer4PU::matchRecTracksToVertex(simPrimaryVertex & pv,
1485  const std::vector<SimPart > & tsim,
1486  const edm::Handle<reco::TrackCollection> & recTrks)
1487 {
1488  // find all recTracks that belong to this simulated vertex (not just the ones that are used in
1489  // matching recVertex)
1490 
1491  std::cout << "dump rec tracks: " << std::endl;
1492  int irec=0;
1493  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1494  t!=recTrks->end(); ++t){
1495  reco::TrackBase::ParameterVector p = t->parameters();
1496  std::cout << irec++ << ") " << p << std::endl;
1497  }
1498 
1499  std::cout << "match sim tracks: " << std::endl;
1500  pv.matchedRecTrackIndex.clear();
1501  pv.nMatchedTracks=0;
1502  int isim=0;
1503  for(std::vector<SimPart>::const_iterator s=tsim.begin();
1504  s!=tsim.end(); ++s){
1505  std::cout << isim++ << " " << s->par;// << std::endl;
1506  int imatch=-1;
1507  int irec=0;
1508  for(reco::TrackCollection::const_iterator t=recTrks->begin();
1509  t!=recTrks->end(); ++t){
1510  reco::TrackBase::ParameterVector p = t->parameters();
1511  if (match(s->par,p)){ imatch=irec; }
1512  irec++;
1513  }
1514  pv.matchedRecTrackIndex.push_back(imatch);
1515  if(imatch>-1){
1516  pv.nMatchedTracks++;
1517  std::cout << " matched to rec trk" << imatch << std::endl;
1518  }else{
1519  std::cout << " not matched" << std::endl;
1520  }
1521  }
1522 }
1523 /********************************************************************************************************/
1524 
1525 
1526 
1527 
1528 /********************************************************************************************************/
1529 
1530 void PrimaryVertexAnalyzer4PU::getTc(const std::vector<reco::TransientTrack>& tracks,
1531  double & Tc, double & chsq, double & dzmax, double & dztrim, double & m4m2){
1532  if (tracks.size()<2){ Tc=-1; chsq=-1; dzmax=-1; dztrim=-1; m4m2=-1; return;}
1533 
1534  double sumw=0, sumwz=0, sumww=0,sumwwz=0,sumwwzz=0;
1535  double zmin=1e10, zmin1=1e10, zmax1=-1e10, zmax=-1e10;
1536  double m4=0,m3=0,m2=0,m1=0,m0=0;
1537  for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
1538  double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
1539  reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
1540  double z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
1541  double dz2= pow((*it).track().dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
1542  // 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);
1543  double w=1./dz2; // take p=1
1544  sumw += w;
1545  sumwz += w*z;
1546  sumwwz += w*w*z;;
1547  sumwwzz += w*w*z*z;
1548  sumww += w*w;
1549  m0 += w;
1550  m1 += w*z;
1551  m2 += w*z*z;
1552  m3 += w*z*z*z;
1553  m4 += w*z*z*z*z;
1554  if(dz2<pow(0.1,2)){
1555  if(z<zmin1){zmin1=z;} if(z<zmin){zmin1=zmin; zmin=z;}
1556  if(z>zmax1){zmax1=z;} if(z>zmax){zmax1=zmax; zmax=z;}
1557  }
1558  }
1559  double z=sumwz/sumw;
1560  double a=sumwwzz-2*z*sumwwz+z*z*sumww;
1561  double b=sumw;
1562  if(tracks.size()>1){
1563  chsq=(m2-m0*z*z)/(tracks.size()-1);
1564  Tc=2.*a/b;
1565  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));
1566  }else{
1567  chsq=0;
1568  Tc=0;
1569  m4m2=0;
1570  }
1571  dzmax=zmax-zmin;
1572  dztrim=zmax1-zmin1;// truncated
1573 }
1574 /********************************************************************************************************/
1575 
1576 
1577 
1578 
1579 
1580 /********************************************************************************************************/
1581 bool PrimaryVertexAnalyzer4PU::truthMatchedTrack( edm::RefToBase<reco::Track> track, TrackingParticleRef & tpr)
1582 
1583 /********************************************************************************************************/
1584 // for a reco track select the matching tracking particle, always use this function to make sure we
1585 // are consistent
1586 // to get the TrackingParticle form the TrackingParticleRef, use ->get();
1587 {
1588  double f=0;
1589  try{
1590  std::vector<std::pair<TrackingParticleRef, double> > tp = r2s_[track];
1591  for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = tp.begin();
1592  it != tp.end(); ++it) {
1593 
1594  if (it->second>f){
1595  tpr=it->first;
1596  f=it->second;
1597  }
1598  }
1599  } catch (Exception event) {
1600  // silly way of testing whether track is in r2s_
1601  }
1602 
1603  // sanity check on track parameters?
1604 
1605  return (f>0.5);
1606 }
1607 /********************************************************************************************************/
1608 
1609 
1610 /********************************************************************************************************/
1611 std::vector< edm::RefToBase<reco::Track> > PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks(
1612  const reco::Vertex& v
1613  )
1614 // for vertex v get a list of tracks for which truth matching is available
1615 /********************************************************************************************************/
1616 {
1617  std::vector< edm::RefToBase<reco::Track> > b;
1618  TrackingParticleRef tpr;
1619 
1620  for(trackit_t tv=v.tracks_begin(); tv!=v.tracks_end(); tv++){
1621  edm::RefToBase<reco::Track> track=*tv;
1622  if (truthMatchedTrack(track, tpr)){
1623  b.push_back(*tv);
1624  }
1625  }
1626 
1627 
1628  return b;
1629 }
1630 /********************************************************************************************************/
1631 
1632 
1633 /********************************************************************************************************/
1634 std::vector<PrimaryVertexAnalyzer4PU::SimEvent> PrimaryVertexAnalyzer4PU::getSimEvents
1635 (
1636  // const Event& iEvent, const EventSetup& iSetup,
1637  edm::Handle<TrackingParticleCollection> TPCollectionH,
1638  edm::Handle<TrackingVertexCollection> TVCollectionH,
1639  edm::Handle<View<Track> > trackCollectionH
1640  ){
1641 
1642  const TrackingParticleCollection* simTracks = TPCollectionH.product();
1643  const View<Track> tC = *(trackCollectionH.product());
1644 
1645 
1646  vector<SimEvent> simEvt;
1647  map<EncodedEventId, unsigned int> eventIdToEventMap;
1648  map<EncodedEventId, unsigned int>::iterator id;
1649  bool dumpTP=false;
1650  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1651 
1652  if( fabs(it->parentVertex().get()->position().z())>100.) continue; // skip funny entries @ -13900
1653 
1654  unsigned int event=0; //note, this is no longer the same as eventId().event()
1655  id=eventIdToEventMap.find(it->eventId());
1656  if (id==eventIdToEventMap.end()){
1657 
1658  // new event here
1659  SimEvent e;
1660  e.eventId=it->eventId();
1661  event=simEvt.size();
1662  const TrackingVertex *parentVertex= it->parentVertex().get();
1663  if(DEBUG_){
1664  cout << "getSimEvents: ";
1665  cout << it->eventId().bunchCrossing() << "," << it->eventId().event()
1666  << " z="<< it->vz() << " "
1667  << parentVertex->eventId().bunchCrossing() << "," <<parentVertex->eventId().event()
1668  << " z=" << parentVertex->position().z() << endl;
1669  }
1670  if (it->eventId()==parentVertex->eventId()){
1671  //e.x=it->vx(); e.y=it->vy(); e.z=it->vz();// wrong ???
1672  e.x=parentVertex->position().x();
1673  e.y=parentVertex->position().y();
1674  e.z=parentVertex->position().z();
1675  if(e.z<-100){
1676  dumpTP=true;
1677  }
1678  }else{
1679  e.x=0; e.y=0; e.z=-88.;
1680  }
1681  simEvt.push_back(e);
1682  eventIdToEventMap[e.eventId]=event;
1683  }else{
1684  event=id->second;
1685  }
1686 
1687 
1688  simEvt[event].tp.push_back(&(*it));
1689  if( (abs(it->eta())<2.5) && (it->charge()!=0) ){
1690  simEvt[event].sumpt2+=pow(it->pt(),2); // should keep track of decays ?
1691  simEvt[event].sumpt+=it->pt();
1692  }
1693  }
1694 
1695  if(dumpTP){
1696  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
1697  std::cout << *it << std::endl;
1698  }
1699  }
1700 
1701 
1702  for(View<Track>::size_type i=0; i<tC.size(); ++i) {
1703  RefToBase<Track> track(trackCollectionH, i);
1704  TrackingParticleRef tpr;
1705  if( truthMatchedTrack(track,tpr)){
1706 
1707  if( eventIdToEventMap.find(tpr->eventId())==eventIdToEventMap.end() ){ cout << "Bug in getSimEvents" << endl; break; }
1708 
1709  z2tp_[track.get()->vz()]=tpr;
1710 
1711  unsigned int event=eventIdToEventMap[tpr->eventId()];
1712  double ipsig=0,ipdist=0;
1713  const TrackingVertex *parentVertex= tpr->parentVertex().get();
1714  double vx=parentVertex->position().x(); // problems with tpr->vz()
1715  double vy=parentVertex->position().y();
1716  double vz=parentVertex->position().z();
1717  double d=sqrt(pow(simEvt[event].x-vx,2)+pow(simEvt[event].y-vy,2)+pow(simEvt[event].z-vz,2))*1.e4;
1718  ipdist=d;
1719  double dxy=track->dxy(vertexBeamSpot_.position());
1720  ipsig=dxy/track->dxyError();
1721 
1722 
1723  TransientTrack t = theB_->build(tC[i]);
1724  t.setBeamSpot(vertexBeamSpot_);
1725  if (theTrackFilter(t)){
1726  simEvt[event].tk.push_back(t);
1727  if(ipdist<5){simEvt[event].tkprim.push_back(t);}
1728  if(ipsig<5){simEvt[event].tkprimsel.push_back(t);}
1729  }
1730 
1731  }
1732  }
1733 
1734 
1735 
1736  AdaptiveVertexFitter theFitter;
1737  cout << "SimEvents " << simEvt.size() << endl;
1738  for(unsigned int i=0; i<simEvt.size(); i++){
1739 
1740  if(simEvt[i].tkprim.size()>0){
1741 
1742  getTc(simEvt[i].tkprimsel, simEvt[i].Tc, simEvt[i].chisq, simEvt[i].dzmax, simEvt[i].dztrim, simEvt[i].m4m2);
1743  TransientVertex v = theFitter.vertex(simEvt[i].tkprim, vertexBeamSpot_);
1744  if (v.isValid()){
1745  simEvt[i].xfit=v.position().x();
1746  simEvt[i].yfit=v.position().y();
1747  simEvt[i].zfit=v.position().z();
1748  // if (simEvt[i].z<-80.){simEvt[i].z=v.position().z();} // for now
1749  }
1750  }
1751 
1752 
1753  if(DEBUG_){
1754  cout << i <<" ) nTP=" << simEvt[i].tp.size()
1755  << " z=" << simEvt[i].z
1756  << " recTrks=" << simEvt[i].tk.size()
1757  << " recTrksPrim=" << simEvt[i].tkprim.size()
1758  << " zfit=" << simEvt[i].zfit
1759  << endl;
1760  }
1761  }
1762 
1763  return simEvt;
1764 }
1765 
1766 
1767 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1768  const Handle<HepMCProduct> evtMC)
1769 {
1770  if(DEBUG_){std::cout << "getSimPVs HepMC " << std::endl;}
1771 
1772  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1773  const HepMC::GenEvent* evt=evtMC->GetEvent();
1774  if (evt) {
1775 // std::cout << "process id " <<evt->signal_process_id()<<std::endl;
1776 // std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
1777 // evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
1778 // std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
1779 
1780 
1781  //int idx=0;
1782  for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
1783  vitr != evt->vertices_end(); ++vitr )
1784  { // loop for vertex ...
1785 
1786  HepMC::FourVector pos = (*vitr)->position();
1787  // if (pos.t()>0) { continue;} // skip secondary vertices, doesn't work for some samples
1788 
1789  if (fabs(pos.z())>1000) continue; // skip funny junk vertices
1790 
1791  bool hasMotherVertex=false;
1792  //std::cout << "mothers" << std::endl;
1793  for ( HepMC::GenVertex::particle_iterator
1794  mother = (*vitr)->particles_begin(HepMC::parents);
1795  mother != (*vitr)->particles_end(HepMC::parents);
1796  ++mother ) {
1797  HepMC::GenVertex * mv=(*mother)->production_vertex();
1798  if (mv) {hasMotherVertex=true;}
1799  //std::cout << "\t"; (*mother)->print();
1800  }
1801  if(hasMotherVertex) {continue;}
1802 
1803 
1804  // could be a new vertex, check all primaries found so far to avoid multiple entries
1805  const double mm=0.1;
1806  simPrimaryVertex sv(pos.x()*mm,pos.y()*mm,pos.z()*mm);
1807  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
1808  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1809  v0!=simpv.end(); v0++){
1810  if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
1811  vp=&(*v0);
1812  break;
1813  }
1814  }
1815  if(!vp){
1816  // this is a new vertex
1817  //std::cout << "this is a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;
1818  simpv.push_back(sv);
1819  vp=&simpv.back();
1820 // }else{
1821 // std::cout << "this is not a new vertex" << std::endl;
1822  }
1823 
1824 
1825  // store the gen vertex barcode with this simpv
1826  vp->genVertex.push_back((*vitr)->barcode());
1827 
1828 
1829  // collect final state descendants and sum up momenta etc
1830  for ( HepMC::GenVertex::particle_iterator
1831  daughter = (*vitr)->particles_begin(HepMC::descendants);
1832  daughter != (*vitr)->particles_end(HepMC::descendants);
1833  ++daughter ) {
1834  //std::cout << "checking daughter type " << (*daughter)->pdg_id() << " final :" <<isFinalstateParticle(*daughter) << std::endl;
1835  if (isFinalstateParticle(*daughter)){
1836  if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
1837  == vp->finalstateParticles.end()){
1838  vp->finalstateParticles.push_back((*daughter)->barcode());
1839  HepMC::FourVector m=(*daughter)->momentum();
1840  //std::cout << "adding particle to primary " << m.px()<< " " << m.py() << " " << m.pz() << std::endl;
1841  vp->ptot.setPx(vp->ptot.px()+m.px());
1842  vp->ptot.setPy(vp->ptot.py()+m.py());
1843  vp->ptot.setPz(vp->ptot.pz()+m.pz());
1844  vp->ptot.setE(vp->ptot.e()+m.e());
1845  vp->ptsq+=(m.perp())*(m.perp());
1846  // count relevant particles
1847  if ( (m.perp()>0.2) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
1848  vp->nGenTrk++;
1849  }
1850 
1851  hsimPV["rapidity"]->Fill(m.pseudoRapidity());
1852  if( (m.perp()>0.8) && isCharged( *daughter ) ){
1853  hsimPV["chRapidity"]->Fill(m.pseudoRapidity());
1854  }
1855  hsimPV["pt"]->Fill(m.perp());
1856  }//new final state particle for this vertex
1857  }//daughter is a final state particle
1858  }
1859 
1860 
1861  //idx++;
1862  }
1863  }
1864  if(verbose_){
1865  cout << "------- PrimaryVertexAnalyzer4PU simPVs -------" << endl;
1866  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
1867  v0!=simpv.end(); v0++){
1868  cout << "z=" << v0->z
1869  << " px=" << v0->ptot.px()
1870  << " py=" << v0->ptot.py()
1871  << " pz=" << v0->ptot.pz()
1872  << " pt2="<< v0->ptsq
1873  << endl;
1874  }
1875  cout << "-----------------------------------------------" << endl;
1876  }
1877  return simpv;
1878 }
1879 
1880 
1881 double
1882 PrimaryVertexAnalyzer4PU::getTrueSeparation(float in_z, const vector<float> & simpv)
1883 {
1884 
1885  double sepL_min = 50.;
1886 
1887  //loop over all sim vertices
1888  for(unsigned sv_idx=0; sv_idx<simpv.size(); ++sv_idx){
1889 
1890  float comp_z = simpv[sv_idx];
1891  double dist_z = fabs(comp_z - in_z);
1892 
1893  if ( dist_z==0. ) continue;
1894 
1895  if ( dist_z<sepL_min ) sepL_min = dist_z;
1896 
1897  }
1898 
1899  return sepL_min;
1900 
1901 }
1902 
1903 
1904 double
1905 PrimaryVertexAnalyzer4PU::getTrueSeparation(SimEvent inEv, vector<SimEvent> & simEv)
1906 {
1907 
1908  EncodedEventId in_ev = inEv.eventId;
1909  double in_z = inEv.z;
1910 
1911  int lastEvent = -1;
1912 
1913  double sepL_min = 50.;
1914 
1915  //loop over all sim vertices
1916  for(unsigned se_idx=0; se_idx<simEv.size(); ++se_idx){
1917 
1918  SimEvent compEv = simEv[se_idx];
1919  EncodedEventId comp_ev = compEv.eventId;
1920 
1921  if ( comp_ev.event()==lastEvent ) continue;
1922  lastEvent = comp_ev.event();
1923 
1924  float comp_z = compEv.z;
1925 
1926  if ( ( fabs(comp_z)>24. ) || ( comp_ev.bunchCrossing()!=0 ) || ( in_ev.event()==comp_ev.event() ) ) continue;
1927 
1928  double dist_z = fabs(comp_z - in_z);
1929 
1930  if ( dist_z<sepL_min ) sepL_min = dist_z;
1931 
1932  }
1933 
1934  return sepL_min;
1935 
1936 }
1937 
1938 vector<int>*
1939 PrimaryVertexAnalyzer4PU::vertex_match(float in_z, const edm::Handle<reco::VertexCollection> vCH)
1940 {
1941 
1942  double zmatch_ = 3.;
1943 
1944  vector<int>* matchs = new vector<int>();
1945 
1946  for(unsigned vtx_idx=0; vtx_idx<vCH->size(); ++vtx_idx){
1947 
1948  VertexRef comp_vtxref(vCH, vtx_idx);
1949 
1950  double comp_z = comp_vtxref->z();
1951  double comp_z_err = comp_vtxref->zError();
1952 
1953  double z_dist = comp_z - in_z;
1954  double z_rel = z_dist / comp_z_err;
1955 
1956  if ( fabs(z_rel)<zmatch_ ) {
1957  matchs->push_back(vtx_idx);
1958  }
1959 
1960  }
1961 
1962  return matchs;
1963 }
1964 
1965 
1966 
1967 
1968 /* get sim pv from TrackingParticles/TrackingVertex */
1969 std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs(
1970  const edm::Handle<TrackingVertexCollection> tVC
1971  )
1972 {
1973  std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
1974  //std::cout <<"number of vertices " << tVC->size() << std::endl;
1975 
1976  if(DEBUG_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;}
1977 
1978  for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) {
1979 
1980  if(DEBUG_){
1981  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;
1982  for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){
1983  std::cout << *gv << std::endl;
1984  }
1985  std::cout << "----------" << std::endl;
1986  }
1987 
1988  // bool hasMotherVertex=false;
1989  if ((unsigned int) v->eventId().event()<simpv.size()) continue;
1990  //if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work)
1991  if (fabs(v->position().z())>1000) continue; // skip funny junk vertices
1992 
1993  // could be a new vertex, check all primaries found so far to avoid multiple entries
1994  const double mm=1.0; // for tracking vertices
1995  simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm);
1996 
1997  //cout << "sv: " << (v->eventId()).event() << endl;
1998  sv.eventId=v->eventId();
1999 
2000  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){
2001  //cout <<((**iTrack).eventId()).event() << " "; // an iterator of Refs, dereference twice. Cool eyh?
2002  sv.eventId=(**iTrack).eventId();
2003  }
2004  //cout <<endl;
2005  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
2006  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
2007  v0!=simpv.end(); v0++){
2008  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)){
2009  vp=&(*v0);
2010  break;
2011  }
2012  }
2013  if(!vp){
2014  // this is a new vertex
2015  if(DEBUG_){std::cout << "this is a new vertex " << sv.eventId.event() << " " << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
2016  simpv.push_back(sv);
2017  vp=&simpv.back();
2018  }else{
2019  if(DEBUG_){std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
2020  }
2021 
2022 
2023  // Loop over daughter track(s)
2024  if(DEBUG_){
2025  for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) {
2026  std::cout << " Daughter momentum: " << (*(*iTP)).momentum();
2027  std::cout << " Daughter type " << (*(*iTP)).pdgId();
2028  std::cout << std::endl;
2029  }
2030  }
2031  }
2032  if(DEBUG_){
2033  cout << "------- PrimaryVertexAnalyzer4PU simPVs from TrackingVertices -------" << endl;
2034  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
2035  v0!=simpv.end(); v0++){
2036  cout << "z=" << v0->z << " event=" << v0->eventId.event() << endl;
2037  }
2038  cout << "-----------------------------------------------" << endl;
2039  }
2040  return simpv;
2041 }
2042 
2043 
2044 
2045 
2046 
2047 
2048 // ------------ method called to produce the data ------------
2049 void
2050 PrimaryVertexAnalyzer4PU::analyze(const Event& iEvent, const EventSetup& iSetup)
2051 {
2052 
2053  std::vector<simPrimaryVertex> simpv; // a list of primary MC vertices
2054  std::vector<float> pui_z;
2055  std::vector<SimPart> tsim;
2056 
2057  eventcounter_++;
2058  run_ = iEvent.id().run();
2059  luminosityBlock_ = iEvent.luminosityBlock();
2060  event_ = iEvent.id().event();
2061  bunchCrossing_ = iEvent.bunchCrossing();
2062  orbitNumber_ = iEvent.orbitNumber();
2063 
2064 
2065 
2066 
2067  if(verbose_){
2068  std::cout << endl
2069  << "PrimaryVertexAnalyzer4PU::analyze event counter=" << eventcounter_
2070  << " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
2071  << " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
2072  //<< " selected = " << good
2073  << std::endl;
2074  }
2075 
2076 
2077  try{
2078  iSetup.getData(pdt_);
2079  }catch(const Exception&){
2080  std::cout << "Some problem occurred with the particle data table. This may not work !" <<std::endl;
2081  }
2082 
2083  try{
2084 
2085  //get the pileup information
2086  Handle< vector<PileupSummaryInfo> > puinfoH;
2087  iEvent.getByToken( vecPileupSummaryInfoToken_, puinfoH );
2088  PileupSummaryInfo puinfo;
2089 
2090  for (unsigned int puinfo_ite=0;puinfo_ite<(*puinfoH).size();++puinfo_ite){
2091  if ((*puinfoH)[puinfo_ite].getBunchCrossing()==0){
2092  puinfo=(*puinfoH)[puinfo_ite];
2093  break;
2094  }
2095  }
2096 
2097  pui_z = puinfo.getPU_zpositions();
2098 
2099  } catch( edm::Exception const & ex ) {
2100  if ( !ex.alreadyPrinted() ) {
2101  std::cout << ex.what() << " Maybe data?" << std::endl;
2102  }
2103  }
2104 
2105  Handle<reco::VertexCollection> recVtxs;
2106  bool bnoBS = iEvent.getByToken( recoVertexCollectionToken_, recVtxs );
2107 
2108  Handle<reco::VertexCollection> recVtxsBS;
2109  bool bBS = iEvent.getByToken( recoVertexCollection_BS_Token_, recVtxsBS );
2110 
2111  Handle<reco::VertexCollection> recVtxsDA;
2112  bool bDA = iEvent.getByToken( recoVertexCollection_DA_Token_, recVtxsDA );
2113 
2114  Handle<reco::TrackCollection> recTrks;
2115  iEvent.getByToken( recoTrackCollectionToken_, recTrks );
2116 
2117  int nhighpurity=0, ntot=0;
2118  for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
2119  ntot++;
2120  if(t->quality(reco::TrackBase::highPurity)) nhighpurity++;
2121  }
2122  if(ntot>10) hnoBS["highpurityTrackFraction"]->Fill(float(nhighpurity)/float(ntot));
2123  if((recoTrackProducer_=="generalTracks")&&(nhighpurity<0.25*ntot)){
2124  if(verbose_){ cout << "rejected, " << nhighpurity << " highpurity out of " << ntot << " total tracks "<< endl<< endl;}
2125  return;
2126  }
2127 
2128 
2129 
2130 
2131 
2132  if(iEvent.getByToken(recoBeamSpotToken_, recoBeamSpotHandle_)){
2133  vertexBeamSpot_= *recoBeamSpotHandle_;
2134  wxy2_=pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2);
2135  Fill(hsimPV, "xbeam",vertexBeamSpot_.x0()); Fill(hsimPV, "wxbeam",vertexBeamSpot_.BeamWidthX());
2136  Fill(hsimPV, "ybeam",vertexBeamSpot_.y0()); Fill(hsimPV, "wybeam",vertexBeamSpot_.BeamWidthY());
2137  Fill(hsimPV, "zbeam",vertexBeamSpot_.z0());// Fill("wzbeam",vertexBeamSpot_.BeamWidthZ());
2138  }else{
2139  cout << " beamspot not found, using dummy " << endl;
2140  vertexBeamSpot_=reco::BeamSpot();// dummy
2141  }
2142 
2143 
2144  if(bnoBS){
2145  if((recVtxs->begin()->isValid())&&(recVtxs->begin()->ndof()>1)&&(recVtxs->begin()->ndof()>(0.0*recVtxs->begin()->tracksSize()))){ // was 5 and 0.2
2146  Fill(hnoBS,"xrecBeamvsdxXBS",recVtxs->begin()->xError(),recVtxs->begin()->x()-vertexBeamSpot_.x0());
2147  Fill(hnoBS,"yrecBeamvsdyXBS",recVtxs->begin()->yError(),recVtxs->begin()->y()-vertexBeamSpot_.y0());
2148 
2149  if(printXBS_) {
2150  cout << Form("XBS %10d %5d %10d %4d %lu %5.2f %8f %8f %8f %8f %8f %8f",
2151  run_,luminosityBlock_,event_,bunchCrossing_,
2152  (unsigned long)(recVtxs->begin()->tracksSize()), recVtxs->begin()->ndof(),
2153  recVtxs->begin()->x(), recVtxs->begin()->xError(),
2154  recVtxs->begin()->y(), recVtxs->begin()->yError(),
2155  recVtxs->begin()->z(), recVtxs->begin()->zError()
2156  ) << endl;
2157  }
2158 
2159  }
2160  }
2161 
2162 
2163  // for the associator
2164  Handle<View<Track> > trackCollectionH;
2165  iEvent.getByToken( edmView_recoTrack_Token_, trackCollectionH );
2166 
2167  Handle<HepMCProduct> evtMC;
2168 
2169  Handle<SimVertexContainer> simVtxs;
2170  iEvent.getByToken( edmSimVertexContainerToken_, simVtxs );
2171 
2172  Handle<SimTrackContainer> simTrks;
2173  iEvent.getByToken( edmSimTrackContainerToken_, simTrks );
2174 
2175 
2176 
2177 
2178  edm::Handle<TrackingParticleCollection> TPCollectionH ;
2179  edm::Handle<TrackingVertexCollection> TVCollectionH ;
2180  bool gotTP = iEvent.getByToken( trackingParticleCollectionToken_, TPCollectionH );
2181  bool gotTV = iEvent.getByToken( trackingVertexCollectionToken_, TVCollectionH );
2182 
2183  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
2184  fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
2185 
2186 
2187  vector<SimEvent> simEvt;
2188  if (gotTP && gotTV ){
2189 
2190  edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
2191  iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
2192  associatorByHits_ = (TrackAssociatorBase *) theHitsAssociator.product();
2193  r2s_ = associatorByHits_->associateRecoToSim (trackCollectionH,TPCollectionH, &iEvent, &iSetup );
2194  //simEvt=getSimEvents(iEvent, TPCollectionH, TVCollectionH, trackCollectionH);
2195  simEvt=getSimEvents(TPCollectionH, TVCollectionH, trackCollectionH);
2196 
2197  if (simEvt.size()==0){
2198  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2199  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2200  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2201  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2202  cout << " !!!!!!!!!!!!!!!!!!!!!! got TrackingParticles but no simEvt !!!!!!!!!!!!!!!!!" << endl;
2203  cout << " dumping Tracking particles " << endl;
2204  const TrackingParticleCollection* simTracks = TPCollectionH.product();
2205  for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
2206  cout << *it << endl;
2207  }
2208  cout << " dumping Tracking Vertices " << endl;
2209  const TrackingVertexCollection* tpVtxs = TVCollectionH.product();
2210  for(TrackingVertexCollection::const_iterator iv=tpVtxs->begin(); iv!=tpVtxs->end(); iv++){
2211  cout << *iv << endl;
2212  }
2213  if( iEvent.getByToken( edmHepMCProductToken_, evtMC ) ){
2214  cout << "dumping simTracks" << endl;
2215  for(SimTrackContainer::const_iterator t=simTrks->begin(); t!=simTrks->end(); ++t){
2216  cout << *t << endl;
2217  }
2218  cout << "dumping simVertexes" << endl;
2219  for(SimVertexContainer::const_iterator vv=simVtxs->begin();
2220  vv!=simVtxs->end();
2221  ++vv){
2222  cout << *vv << endl;
2223  }
2224  }else{
2225  cout << "no hepMC" << endl;
2226  }
2227  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2228 
2229  const HepMC::GenEvent* evt=evtMC->GetEvent();
2230  if(evt){
2231  std::cout << "process id " <<evt->signal_process_id()<<std::endl;
2232  std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
2233  evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
2234  std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
2235  evt->print();
2236  }else{
2237  cout << "no event in HepMC" << endl;
2238  }
2239  cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
2240 
2241  }
2242  }
2243 
2244 
2245 
2246 
2247  if(gotTV){
2248 
2249  if(verbose_){
2250  cout << "Found Tracking Vertices " << endl;
2251  }
2252  simpv=getSimPVs(TVCollectionH);
2253 
2254 
2255  }else if( iEvent.getByToken( edmHepMCProductToken_, evtMC ) ) {
2256 
2257  simpv=getSimPVs(evtMC);
2258 
2259  if(verbose_){
2260  cout << "Using HepMCProduct " << endl;
2261  std::cout << "simtrks " << simTrks->size() << std::endl;
2262  }
2263  tsim = PrimaryVertexAnalyzer4PU::getSimTrkParameters(simTrks, simVtxs, simUnit_);
2264 
2265  }else{
2266  // if(verbose_({cout << "No MC info at all" << endl;}
2267  }
2268 
2269 
2270 
2271 
2272  // get the beam spot from the appropriate dummy vertex, if available
2273  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2274  v!=recVtxs->end(); ++v){
2275  if ((v->ndof()==-3) && (v->chi2()==0) ){
2276  myBeamSpot=math::XYZPoint(v->x(), v->y(), v->z());
2277  }
2278  }
2279 
2280 
2281 
2282 
2283  hsimPV["nsimvtx"]->Fill(simpv.size());
2284  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
2285  vsim!=simpv.end(); vsim++){
2286  if(doMatching_){
2287  matchRecTracksToVertex(*vsim, tsim, recTrks); // hepmc, based on track parameters
2288  }
2289 
2290  hsimPV["nbsimtksinvtx"]->Fill(vsim->nGenTrk);
2291  hsimPV["xsim"]->Fill(vsim->x*simUnit_);
2292  hsimPV["ysim"]->Fill(vsim->y*simUnit_);
2293  hsimPV["zsim"]->Fill(vsim->z*simUnit_);
2294  hsimPV["xsim1"]->Fill(vsim->x*simUnit_);
2295  hsimPV["ysim1"]->Fill(vsim->y*simUnit_);
2296  hsimPV["zsim1"]->Fill(vsim->z*simUnit_);
2297  Fill(hsimPV,"xsim2",vsim->x*simUnit_,vsim==simpv.begin());
2298  Fill(hsimPV,"ysim2",vsim->y*simUnit_,vsim==simpv.begin());
2299  Fill(hsimPV,"zsim2",vsim->z*simUnit_,vsim==simpv.begin());
2300  hsimPV["xsim2"]->Fill(vsim->x*simUnit_);
2301  hsimPV["ysim2"]->Fill(vsim->y*simUnit_);
2302  hsimPV["zsim2"]->Fill(vsim->z*simUnit_);
2303  hsimPV["xsim3"]->Fill(vsim->x*simUnit_);
2304  hsimPV["ysim3"]->Fill(vsim->y*simUnit_);
2305  hsimPV["zsim3"]->Fill(vsim->z*simUnit_);
2306  hsimPV["xsimb"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
2307  hsimPV["ysimb"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
2308  hsimPV["zsimb"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
2309  hsimPV["xsimb1"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
2310  hsimPV["ysimb1"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
2311  hsimPV["zsimb1"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
2312  }
2313 
2314 
2315 
2316 
2317 
2318  if(bnoBS){
2319  analyzeVertexCollection(hnoBS, recVtxs, recTrks, simpv, pui_z, "noBS");
2320  analyzeVertexCollectionTP(hnoBS, recVtxs, recTrks, simEvt, r2s_,"noBS");
2321  }
2322  if(bBS){
2323  analyzeVertexCollection(hBS, recVtxsBS, recTrks, simpv, pui_z, "BS");
2324  analyzeVertexCollectionTP(hBS, recVtxsBS, recTrks, simEvt, r2s_,"BS");
2325  }
2326  if(bDA){
2327  analyzeVertexCollection(hDA, recVtxsDA, recTrks, simpv, pui_z, "DA");
2328  analyzeVertexCollectionTP(hDA, recVtxsDA, recTrks, simEvt, r2s_,"DA");
2329  }
2330 // if(bPIX){
2331 // analyzeVertexCollection(hPIX, recVtxsPIX, recTrks, simpv,"PIX");
2332 // analyzeVertexCollectionTP(hPIX, recVtxsPIX, recTrks, simEvt,"PIX");
2333 // }
2334 
2335 
2336 
2337  if((dumpThisEvent_&& (dumpcounter_<100)) ||(verbose_ && (eventcounter_<ndump_))){
2338  cout << endl << "Event dump" << endl
2339  << "event counter=" << eventcounter_
2340  << " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
2341  << " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
2342  << std::endl;
2343  dumpcounter_++;
2344 
2345  //evtMC->GetEvent()->print();
2346  //printRecTrks(recTrks); // very verbose !!
2347 
2348 // if (bPIX) printRecVtxs(recVtxsPIX,"pixel vertices");
2349  if (bnoBS) {printRecVtxs(recVtxs,"Offline without Beamspot");}
2350  if (bnoBS && (!bDA)){ printPVTrks(recTrks, recVtxs, tsim, simEvt, false);}
2351  if (bBS) printRecVtxs(recVtxsBS,"Offline with Beamspot");
2352  if (bDA) {
2353  printRecVtxs(recVtxsDA,"Offline DA");
2354  printPVTrks(recTrks, recVtxsDA, tsim, simEvt, false);
2355  }
2356  if (dumpcounter_<2){cout << "beamspot " << vertexBeamSpot_ << endl;}
2357  }
2358 
2359  if(verbose_){
2360  std::cout << std::endl;
2361  }
2362 }
2363 
2364 namespace {
2365 bool lt(const std::pair<double,unsigned int>& a,const std::pair<double,unsigned int>& b ){
2366  return a.first<b.first;
2367 }
2368 }
2369 
2370 /***************************************************************************************/
2371 void PrimaryVertexAnalyzer4PU::printEventSummary(std::map<std::string, TH1*> & h,
2372  const edm::Handle<reco::VertexCollection> recVtxs,
2373  const edm::Handle<reco::TrackCollection> recTrks,
2374  vector<SimEvent> & simEvt,
2375  const string message){
2376  // make a readable summary of the vertex finding if the TrackingParticles are availabe
2377  if (simEvt.size()==0){return;}
2378 
2379 
2380  // sort vertices in z ... for nicer printout
2381 
2382  vector< pair<double,unsigned int> > zrecv;
2383  for(unsigned int idx=0; idx<recVtxs->size(); idx++){
2384  if ( (recVtxs->at(idx).ndof()<0) || (recVtxs->at(idx).chi2()<=0) ) continue; // skip clusters
2385  zrecv.push_back( make_pair(recVtxs->at(idx).z(),idx) );
2386  }
2387  stable_sort(zrecv.begin(),zrecv.end(),lt);
2388 
2389  // same for simulated vertices
2390  vector< pair<double,unsigned int> > zsimv;
2391  for(unsigned int idx=0; idx<simEvt.size(); idx++){
2392  zsimv.push_back(make_pair(simEvt[idx].z, idx));
2393  }
2394  stable_sort(zsimv.begin(), zsimv.end(),lt);
2395 
2396 
2397 
2398 
2399  cout << "---------------------------" << endl;
2400  cout << "event counter = " << eventcounter_ << " " << message << endl;
2401  cout << "---------------------------" << endl;
2402  cout << " z[cm] rec --> ";
2403  cout.precision(4);
2404  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2405  cout << setw(7) << fixed << itrec->first;
2406  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2407  }
2408  cout << endl;
2409  cout << " ";
2410  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2411  cout << setw(7) << fixed << recVtxs->at(itrec->second).tracksSize();
2412  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2413  }
2414  cout << " rec tracks" << endl;
2415  cout << " ";
2416  map<unsigned int, int> truthMatchedVertexTracks;
2417  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2418  truthMatchedVertexTracks[itrec->second]=getTruthMatchedVertexTracks(recVtxs->at(itrec->second)).size();
2419  cout << setw(7) << fixed << truthMatchedVertexTracks[itrec->second];
2420  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2421  }
2422  cout << " truth matched " << endl;
2423 
2424  cout << "sim ------- trk prim ----" << endl;
2425 
2426 
2427 
2428  map<unsigned int, unsigned int> rvmatch; // reco vertex matched to sim vertex (sim to rec)
2429  map<unsigned int, double > nmatch; // highest number of truth-matched tracks of ev found in a recvtx
2430  map<unsigned int, double > purity; // highest purity of a rec vtx (i.e. highest number of tracks from the same simvtx)
2431  map<unsigned int, double > wpurity; // same for the sum of weights
2432 
2433  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2434  purity[itrec->second]=0.;
2435  wpurity[itrec->second]=0.;
2436  }
2437 
2438  for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
2439  SimEvent* ev =&(simEvt[itsim->second]);
2440 
2441 
2442  cout.precision(4);
2443  if (itsim->second==0){
2444  cout << setw(8) << fixed << ev->z << ")*" << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
2445  }else{
2446  cout << setw(8) << fixed << ev->z << ") " << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
2447  }
2448 
2449  nmatch[itsim->second]=0; // highest number of truth-matched tracks of ev found in a recvtx
2450  double matchpurity=0,matchwpurity=0;
2451 
2452  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2453  const reco::Vertex *v = &(recVtxs->at(itrec->second));
2454 
2455  // count tracks found in both, sim and rec
2456  double n=0,wt=0;
2457  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2458  const reco::Track& RTe=te->track();
2459  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2460  const reco::Track & RTv=*(tv->get());
2461  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2462  }
2463  }
2464  cout << setw(7) << int(n)<< " ";
2465 
2466  if (n > nmatch[itsim->second]){
2467  nmatch[itsim->second]=n;
2468  rvmatch[itsim->second]=itrec->second;
2469  matchpurity=n/truthMatchedVertexTracks[itrec->second];
2470  matchwpurity=wt/truthMatchedVertexTracks[itrec->second];
2471  }
2472 
2473  if(n > purity[itrec->second]){
2474  purity[itrec->second]=n;
2475  }
2476 
2477  if(wt > wpurity[itrec->second]){
2478  wpurity[itrec->second]=wt;
2479  }
2480 
2481  }// end of reco vertex loop
2482 
2483  cout << " | ";
2484  if (nmatch[itsim->second]>0 ){
2485  if(matchpurity>0.5){
2486  cout << "found ";
2487  }else{
2488  cout << "merged ";
2489  }
2490  cout << " max eff. = " << setw(8) << nmatch[itsim->second]/ev->tk.size() << " p=" << matchpurity << " w=" << matchwpurity << endl;
2491  }else{
2492  if(ev->tk.size()==0){
2493  cout << " invisible" << endl;
2494  }else if (ev->tk.size()==1){
2495  cout << "single track " << endl;
2496  }else{
2497  cout << "lost " << endl;
2498  }
2499  }
2500  }
2501  cout << "---------------------------" << endl;
2502 
2503  // the purity of the reconstructed vertex
2504  cout << " purity ";
2505  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2506  cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
2507  if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2508  }
2509  cout << endl;
2510 
2511 // // classification of reconstructed vertex fake/real
2512 // cout << " ";
2513 // for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2514 // cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
2515 // if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
2516 // }
2517 // cout << endl;
2518  cout << "---------------------------" << endl;
2519 
2520 
2521 
2522 
2523  // list problematic tracks
2524  for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
2525  SimEvent* ev =&(simEvt[itsim->second]);
2526 
2527  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2528  const reco::Track& RTe=te->track();
2529 
2530  int ivassign=-1; // will become the index of the vertex to which a track was assigned
2531 
2532  for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
2533  const reco::Vertex *v = &(recVtxs->at(itrec->second));
2534 
2535  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2536  const reco::Track & RTv=*(tv->get());
2537  if(RTe.vz()==RTv.vz()) {ivassign=itrec->second;}
2538  }
2539  }
2540  double tantheta=tan((te->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
2541  reco::BeamSpot beamspot=(te->stateAtBeamLine()).beamSpot();
2542  //double z=(te->stateAtBeamLine().trackStateAtPCA()).position().z();
2543  double dz2= pow(RTe.dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
2544 
2545  if(ivassign==(int)rvmatch[itsim->second]){
2546  Fill(h,"correctlyassigned",RTe.eta(),RTe.pt());
2547  Fill(h,"ptcat",RTe.pt());
2548  Fill(h,"etacat",RTe.eta());
2549  Fill(h,"phicat",RTe.phi());
2550  Fill(h,"dzcat",sqrt(dz2));
2551  }else{
2552  Fill(h,"misassigned",RTe.eta(),RTe.pt());
2553  Fill(h,"ptmis",RTe.pt());
2554  Fill(h,"etamis",RTe.eta());
2555  Fill(h,"phimis",RTe.phi());
2556  Fill(h,"dzmis",sqrt(dz2));
2557  cout << "vertex " << setw(8) << fixed << ev->z;
2558 
2559  if (ivassign<0){
2560  cout << " track lost ";
2561  // for some clusterizers there shouldn't be any lost tracks,
2562  // are there differences in the track selection?
2563  }else{
2564  cout << " track misassigned " << setw(8) << fixed << recVtxs->at(ivassign).z();
2565  }
2566 
2567  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);
2568 
2569  //
2570  //cout << " ztrack=" << te->track().vz();
2571  TrackingParticleRef tpr = z2tp_[te->track().vz()];
2572  double zparent=tpr->parentVertex().get()->position().z();
2573  if(zparent==ev->z) {
2574  cout << " prim";
2575  }else{
2576  cout << " sec ";
2577  }
2578  cout << " id=" << tpr->pdgId();
2579  cout << endl;
2580 
2581  //
2582  }
2583  }// next simvertex-track
2584 
2585  }//next simvertex
2586 
2587  cout << "---------------------------" << endl;
2588 
2589 }
2590 /***************************************************************************************/
2591 
2592 
2593 
2594 
2595 /***************************************************************************************/
2596 void PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP(std::map<std::string, TH1*> & h,
2597  const edm::Handle<reco::VertexCollection> recVtxs,
2598  const edm::Handle<reco::TrackCollection> recTrks,
2599  vector<SimEvent> & simEvt,
2600  reco::RecoToSimCollection rsC,
2601  const string message){
2602 
2603  // cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP size=" << simEvt.size() << endl;
2604  if(simEvt.size()==0)return;
2605 
2606  printEventSummary(h, recVtxs,recTrks,simEvt,message);
2607 
2608  //const int iSignal=0;
2609  EncodedEventId iSignal=simEvt[0].eventId;
2610  Fill(h,"npu0",simEvt.size());
2611 
2612 
2613  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2614  Fill(h,"Tc", ev->Tc, ev==simEvt.begin());
2615  Fill(h,"Chisq", ev->chisq, ev==simEvt.begin());
2616  if(ev->chisq>0)Fill(h,"logChisq", log(ev->chisq), ev==simEvt.begin());
2617  Fill(h,"dzmax", ev->dzmax, ev==simEvt.begin());
2618  Fill(h,"dztrim",ev->dztrim,ev==simEvt.begin());
2619  Fill(h,"m4m2", ev->m4m2, ev==simEvt.begin());
2620  if(ev->Tc>0){ Fill(h,"logTc",log(ev->Tc)/log(10.),ev==simEvt.begin());}
2621 
2622 
2623  for(vector<SimEvent>::iterator ev2=ev+1; ev2!=simEvt.end(); ev2++){
2624  vector<TransientTrack> xt;
2625  if((ev->tkprimsel.size()>0)&&(ev2->tkprimsel.size()>0)&&(ev->tkprimsel.size()+ev2->tkprimsel.size())>1){
2626  xt.insert (xt.end() ,ev->tkprimsel.begin(),ev->tkprimsel.end());
2627  xt.insert (xt.end() ,ev2->tkprimsel.begin(),ev2->tkprimsel.end());
2628  double xTc,xChsq,xDzmax,xDztrim,xm4m2;
2629  getTc(xt, xTc, xChsq, xDzmax, xDztrim,xm4m2);
2630  if(xTc>0){
2631  Fill(h,"xTc",xTc,ev==simEvt.begin());
2632  Fill(h,"logxTc", log(xTc)/log(10),ev==simEvt.begin());
2633  Fill(h,"xChisq", xChsq,ev==simEvt.begin());
2634  if(xChsq>0){Fill(h,"logxChisq", log(xChsq),ev==simEvt.begin());};
2635  Fill(h,"xdzmax", xDzmax,ev==simEvt.begin());
2636  Fill(h,"xdztrim", xDztrim,ev==simEvt.begin());
2637  Fill(h,"xm4m2", xm4m2,ev==simEvt.begin());
2638 
2639  }
2640  }
2641  }
2642  }
2643 
2644  // --------------------------------------- count actual rec vtxs ----------------------
2645  int nrecvtxs=0;//, nrecvtxs1=0, nrecvtxs2=0;
2646  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2647  v!=recVtxs->end(); ++v){
2648  if ( (v->isFake()) || (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
2649  nrecvtxs++;
2650  }
2651 
2652  // --------------------------------------- fill the track assignment matrix ----------------------
2653  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2654  ev->ntInRecVz.clear(); // just in case
2655  ev->zmatch=-99.;
2656  ev->nmatch=0;
2657  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2658  v!=recVtxs->end(); ++v){
2659  double n=0, wt=0;
2660  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2661  const reco::Track& RTe=te->track();
2662  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2663  const reco::Track & RTv=*(tv->get());
2664  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2665  }
2666  }
2667  ev->ntInRecVz[v->z()]=n;
2668  if (n > ev->nmatch){ ev->nmatch=n; ev->zmatch=v->z(); ev->zmatch=v->z(); }
2669  }
2670  }
2671 
2672  // call a vertex a fake if for every sim vertex there is another recvertex containing more tracks
2673  // from that sim vertex than the current recvertex
2674  double nfake=0;
2675  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2676  v!=recVtxs->end(); ++v){
2677  bool matched=false;
2678  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2679  if ((ev->nmatch>0)&&(ev->zmatch==v->z())){
2680  matched=true;
2681  }
2682  }
2683  if(!matched && !v->isFake()) {
2684  nfake++;
2685  cout << " fake rec vertex at z=" << v->z() << endl;
2686  // some histograms of fake vertex properties here
2687  Fill(h,"unmatchedVtxZ",v->z());
2688  Fill(h,"unmatchedVtxNdof",v->ndof());
2689  }
2690  }
2691  if(nrecvtxs>0){
2692  Fill(h,"unmatchedVtx",nfake);
2693  Fill(h,"unmatchedVtxFrac",nfake/nrecvtxs);
2694  }
2695 
2696  // --------------------------------------- match rec to sim ---------------------------------------
2697  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2698  v!=recVtxs->end(); ++v){
2699 
2700  if ( (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
2701  double nmatch=-1; // highest number of tracks in recvtx v found in any event
2702  EncodedEventId evmatch;
2703  bool matchFound=false;
2704  double nmatchvtx=0; // number of simvtcs contributing to recvtx v
2705 
2706  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2707 
2708  double n=0; // number of tracks that are in both, the recvtx v and the event ev
2709  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2710 
2711  const reco::Track& RTe=te->track();
2712  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2713  const reco::Track & RTv=*(tv->get());
2714  if(RTe.vz()==RTv.vz()){ n++;}
2715  }
2716  }
2717 
2718  // find the best match in terms of the highest number of tracks
2719  // from a simvertex in this rec vertex
2720  if (n > nmatch){
2721  nmatch=n;
2722  evmatch=ev->eventId;
2723  matchFound=true;
2724  }
2725  if(n>0){
2726  nmatchvtx++;
2727  }
2728  }
2729 
2730  double nmatchany=getTruthMatchedVertexTracks(*v).size();
2731  if (matchFound && (nmatchany>0)){
2732  // highest number of tracks in recvtx matched to (the same) sim vertex
2733  // purity := -----------------------------------------------------------------
2734  // number of truth matched tracks in this recvtx
2735  double purity =nmatch/nmatchany;
2736  Fill(h,"recmatchPurity",purity);
2737  if(v==recVtxs->begin()){
2738  Fill(h,"recmatchPurityTag",purity, (bool)(evmatch==iSignal));
2739  }else{
2740  Fill(h,"recmatchPuritynoTag",purity,(bool)(evmatch==iSignal));
2741  }
2742  }
2743  Fill(h,"recmatchvtxs",nmatchvtx);
2744  if(v==recVtxs->begin()){
2745  Fill(h,"recmatchvtxsTag",nmatchvtx);
2746  }else{
2747  Fill(h,"recmatchvtxsnoTag",nmatchvtx);
2748  }
2749 
2750 
2751 
2752  } // recvtx loop
2753  Fill(h,"nrecv",nrecvtxs);
2754 
2755 
2756  // --------------------------------------- match sim to rec ---------------------------------------
2757 
2758  int npu1=0, npu2=0;
2759 
2760  vector<int> used_indizesV;
2761  int lastEvent = 999;
2762 
2763  for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
2764 
2765  if(ev->tk.size()>0) npu1++;
2766  if(ev->tk.size()>1) npu2++;
2767 
2768  double sep = getTrueSeparation(*ev, simEvt);
2769 
2770  bool isSignal= ev->eventId==iSignal;
2771 
2772  Fill(h,"nRecTrkInSimVtx",(double) ev->tk.size(),isSignal);
2773  Fill(h,"nPrimRecTrkInSimVtx",(double) ev->tkprim.size(),isSignal);
2774  Fill(h,"sumpt2rec",sqrt(ev->sumpt2rec),isSignal);
2775  Fill(h,"sumpt2",sqrt(ev->sumpt2),isSignal);
2776  Fill(h,"sumpt",sqrt(ev->sumpt),isSignal);
2777 
2778  double nRecVWithTrk=0; // vertices with tracks from this simvertex
2779  double nmatch=0, ntmatch=0, zmatch=-99;
2780 
2781  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
2782  v!=recVtxs->end(); ++v){
2783  if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
2784  // count tracks found in both, sim and rec
2785  double n=0, wt=0;
2786  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2787  const reco::Track& RTe=te->track();
2788  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2789  const reco::Track & RTv=*(tv->get());
2790  if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
2791  }
2792  }
2793 
2794  if(n>0){ nRecVWithTrk++; }
2795  if (n > nmatch){
2796  nmatch=n; ntmatch=v->tracksSize(); zmatch=v->position().z();
2797  }
2798 
2799  }// end of reco vertex loop
2800 
2801 
2802  // nmatch is the highest number of tracks from this sim vertex found in a single reco-vertex
2803  if(ev->tk.size()>0){ Fill(h,"trkAssignmentEfficiency", nmatch/ev->tk.size(), isSignal); };
2804  if(ev->tkprim.size()>0){ Fill(h,"primtrkAssignmentEfficiency", nmatch/ev->tkprim.size(), isSignal); };
2805 
2806  // matched efficiency = efficiency for finding a reco vertex with > 50% of the simvertexs reconstructed tracks
2807 
2808  double ntsim=ev->tk.size(); // may be better to use the number of primary tracks here ?
2809  double matchpurity=nmatch/ntmatch;
2810 
2811  if(ntsim>0){
2812 
2813  Fill(h,"matchVtxFraction",nmatch/ntsim,isSignal);
2814  if(nmatch/ntsim>=0.5){
2815  Fill(h,"matchVtxEfficiency",1.,isSignal);
2816  if(ntsim>1){Fill(h,"matchVtxEfficiency2",1.,isSignal);}
2817  if(matchpurity>0.5){Fill(h,"matchVtxEfficiency5",1.,isSignal);}
2818  }else{
2819  Fill(h,"matchVtxEfficiency",0.,isSignal);
2820  if(ntsim>1){Fill(h,"matchVtxEfficiency2",0.,isSignal);}
2821  Fill(h,"matchVtxEfficiency5",0.,isSignal); // no (matchpurity>5) here !!
2822  if(isSignal){
2823  cout << "Signal vertex not matched " << message << " event=" << eventcounter_ << " nmatch=" << nmatch << " ntsim=" << ntsim << endl;
2824  }
2825  }
2826  } // ntsim >0
2827 
2828 
2829  if(zmatch>-99){
2830  Fill(h,"matchVtxZ",zmatch-ev->z);
2831  Fill(h,"matchVtxZ",zmatch-ev->z,isSignal);
2832  Fill(h,"matchVtxZCum",fabs(zmatch-ev->z));
2833  Fill(h,"matchVtxZCum",fabs(zmatch-ev->z),isSignal);
2834 
2835  }else{
2836  Fill(h,"matchVtxZCum",1.0);
2837  Fill(h,"matchVtxZCum",1.0,isSignal);
2838 
2839  }
2840  if(fabs(zmatch-ev->z)<zmatch_){
2841  Fill(h,"matchVtxEfficiencyZ",1.,isSignal);
2842  }else{
2843  Fill(h,"matchVtxEfficiencyZ",0.,isSignal);
2844  }
2845 
2846  if(ntsim>0) Fill(h, "matchVtxEfficiencyZ1", fabs(zmatch-ev->z)<zmatch_ , isSignal);
2847  if(ntsim>1) Fill(h, "matchVtxEfficiencyZ2", fabs(zmatch-ev->z)<zmatch_ , isSignal);
2848 
2849 
2850  Fill(h,"vtxMultiplicity",nRecVWithTrk,isSignal);
2851 
2852  // efficiency vs number of tracks, use your favorite definition of efficiency here
2853  //if(nmatch>=0.5*ntmatch){ // purity
2854  if(fabs(zmatch-ev->z)<zmatch_){ // zmatch
2855  Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),1.);
2856  if(isSignal){
2857  Fill(h,"vtxFindingEfficiencyVsNtrkSignal",ev->tk.size(),1.);
2858  }else{
2859  Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
2860  }
2861 
2862  }else{
2863  Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),0.);
2864  if(isSignal){
2865  Fill(h,"vtxFindingEfficiencyVsNtrkSignal",(double) ev->tk.size(),1.);
2866  }else{
2867  Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
2868  }
2869  }
2870 
2871  //part of the separation efficiency profiles
2872 
2873  if ( ev->eventId.event()==lastEvent ) continue;
2874  lastEvent = ev->eventId.event();
2875 
2876  if ( ( fabs(ev->z)>24. ) || ( ev->eventId.bunchCrossing()!=0 ) ) continue;
2877 
2878  int max_match = 0;
2879  int best_match = 999;
2880 
2881  for ( unsigned rv_idx=0; rv_idx<recVtxs->size(); rv_idx++ ) {
2882 
2883  VertexRef vtx_ref(recVtxs, rv_idx);
2884 
2885  int match = 0;
2886 
2887  for ( vector<TrackBaseRef>::const_iterator rv_trk_ite=vtx_ref->tracks_begin(); rv_trk_ite!=vtx_ref->tracks_end(); rv_trk_ite++ ) {
2888 
2889  vector<pair<TrackingParticleRef, double> > tp;
2890  if ( rsC.find(*rv_trk_ite)!=rsC.end() ) tp = rsC[*rv_trk_ite];
2891 
2892  for ( unsigned matched_tp_idx=0; matched_tp_idx<tp.size(); matched_tp_idx++ ) {
2893 
2894  TrackingParticle trackpart = *(tp[matched_tp_idx].first);
2895  EncodedEventId tp_ev = trackpart.eventId();
2896 
2897  if ( ( tp_ev.bunchCrossing()==ev->eventId.bunchCrossing() ) && ( tp_ev.event()==ev->eventId.event() ) ) {
2898  match++;
2899  break;
2900  }
2901 
2902  }
2903 
2904  }
2905 
2906 
2907  if ( match > max_match ) {
2908 
2909  max_match = match;
2910  best_match = rv_idx;
2911 
2912  }
2913 
2914  }
2915 
2916  double eff = 0.;
2917  double effwod = 0.;
2918  double dsimed = 0.;
2919 
2920  bool dsflag = false;
2921 
2922  for (unsigned i=0;i<used_indizesV.size(); i++) {
2923  if ( used_indizesV.at(i)==best_match ) {
2924  dsflag = true;
2925  break;
2926  }
2927  }
2928 
2929  if ( best_match!=999 ) eff = 1.;
2930  if ( dsflag ) dsimed = 1.;
2931  if ( ( best_match!=999 ) && ( !dsflag ) ) effwod = 1.;
2932 
2933  if ( best_match!=999 ) {
2934  used_indizesV.push_back(best_match);
2935  }
2936 
2937  Fill(h,"tveffvszsep", sep, eff);
2938  Fill(h,"tveffwodvszsep", sep, effwod);
2939  Fill(h,"tvmergedvszsep", sep, dsimed);
2940 
2941 
2942  }
2943 
2944  Fill(h,"npu1",npu1);
2945  Fill(h,"npu2",npu2);
2946 
2947  Fill(h,"nrecvsnpu",npu1,float(nrecvtxs));
2948  Fill(h,"nrecvsnpu2",npu2,float(nrecvtxs));
2949 
2950  // --------------------------------------- sim-signal vs rec-tag ---------------------------------------
2951  SimEvent* ev=&(simEvt[0]);
2952  const reco::Vertex* v=&(*recVtxs->begin());
2953 
2954  double n=0;
2955  for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
2956  const reco::Track& RTe=te->track();
2957  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
2958  const reco::Track & RTv=*(tv->get());
2959  if(RTe.vz()==RTv.vz()) {n++;}
2960  }
2961  }
2962 
2963  cout << "Number of tracks in reco tagvtx " << v->tracksSize() << endl;
2964  cout << "Number of selected tracks in sim event vtx " << ev->tk.size() << " (prim=" << ev->tkprim.size() << ")"<<endl;
2965  cout << "Number of tracks in both " << n << endl;
2966  double ntruthmatched=getTruthMatchedVertexTracks(*v).size();
2967  if (ntruthmatched>0){
2968  cout << "TrackPurity = "<< n/ntruthmatched <<endl;
2969  Fill(h,"TagVtxTrkPurity",n/ntruthmatched);
2970  }
2971  if (ev->tk.size()>0){
2972  cout << "TrackEfficiency = "<< n/ev->tk.size() <<endl;
2973  Fill(h,"TagVtxTrkEfficiency",n/ev->tk.size());
2974  }
2975 }
2976 
2977 /***************************************************************************************/
2978 
2979 
2980 
2981 
2982 
2983 
2984 /***************************************************************************************/
2985 
2986 void PrimaryVertexAnalyzer4PU::analyzeVertexCollection(std::map<std::string, TH1*> & h,
2987  const Handle<reco::VertexCollection> recVtxs,
2988  const Handle<reco::TrackCollection> recTrks,
2989  std::vector<simPrimaryVertex> & simpv,
2990  const std::vector<float> & pui_z,
2991  const std::string message)
2992 {
2993  //cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollection (HepMC), simpvs=" << simpv.size() << endl;
2994  int nrectrks=recTrks->size();
2995  int nrecvtxs=recVtxs->size();
2996  int nseltrks=-1;
2997  reco::TrackCollection selTrks; // selected tracks
2998  reco::TrackCollection lostTrks; // selected but lost tracks (part of dropped clusters)
2999 
3000  // extract dummy vertices representing clusters
3001  reco::VertexCollection clusters;
3002  reco::Vertex allSelected;
3003  double cpufit=0;
3004  double cpuclu=0;
3005  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3006  v!=recVtxs->end(); ++v){
3007  if ( (fabs(v->ndof()+3.)<0.0001) && (v->chi2()<=0) ){
3008  // this dummy vertex is for the full event
3009  allSelected=(*v);
3010  nseltrks=(allSelected.tracksSize());
3011  nrecvtxs--;
3012  cpuclu=-v->chi2();
3013  continue;
3014  }else if( (fabs(v->ndof()+2.)<0.0001) && (v->chi2()==0) ){
3015  // this is a cluster, not a vertex
3016  clusters.push_back(*v);
3017  Fill(h,"cpuvsntrk",(double) v->tracksSize(),fabs(v->y()));
3018  cpufit+=fabs(v->y());
3019  Fill(h,"nclutrkall",(double) v->tracksSize());
3020  Fill(h,"selstat",v->x());
3021  //Fill(h,"nclutrkvtx",);// see below
3022  nrecvtxs--;
3023  }
3024  }
3025  Fill(h,"cpuclu",cpuclu);
3026  Fill(h,"cpufit",cpufit);
3027  Fill(h,"cpucluvsntrk",nrectrks, cpuclu);
3028 
3029 
3030 
3031  if(simpv.size()>0){//this is mc
3032  double dsimrecx=0.;
3033  double dsimrecy=0.;//0.0011;
3034  double dsimrecz=0.;//0.0012;
3035 
3036  // vertex matching and efficiency bookkeeping
3037  int nsimtrk=0;
3038  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
3039  vsim!=simpv.end(); vsim++){
3040 
3041  nsimtrk+=vsim->nGenTrk;
3042  // look for a matching reconstructed vertex
3043  vsim->recVtx=NULL;
3044  vsim->cluster=-1;
3045 
3046  for(reco::VertexCollection::const_iterator vrec=recVtxs->begin(); vrec!=recVtxs->end(); ++vrec){
3047 
3048  if( vrec->isFake() ) {
3049  continue; // skip fake vertices (=beamspot)
3050  cout << "fake vertex" << endl;
3051  }
3052 
3053  if( vrec->ndof()<0. )continue; // skip dummy clusters, if any
3054  // if ( matchVertex(*vsim,*vrec) ){
3055 
3056  // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
3057  if( ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z-dsimrecz)>fabs(vrec->z()-vsim->z-dsimrecz)))
3058  || (!vsim->recVtx) )
3059  {
3060  vsim->recVtx=&(*vrec);
3061 
3062  // find the corresponding cluster
3063  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3064  if( fabs(clusters[iclu].position().z()-vrec->position().z()) < 0.001 ){
3065  vsim->cluster=iclu;
3066  vsim->nclutrk=clusters[iclu].position().y();
3067  }
3068  }
3069  }
3070 
3071  // the following only works in MC samples without pile-up
3072  if ((simpv.size()==1) && ( fabs(vsim->recVtx->position().z()-vsim->z-dsimrecz)>zmatch_ )){
3073 
3074  // now we have a recvertex without a matching simvertex, I would call this fake
3075  // however, the G4 info does not contain pile-up
3076  Fill(h,"fakeVtxZ",vrec->z());
3077  if (vrec->ndof()>=0.5) Fill(h,"fakeVtxZNdofgt05",vrec->z());
3078  if (vrec->ndof()>=2.0) Fill(h,"fakeVtxZNdofgt2",vrec->z());
3079  Fill(h,"fakeVtxNdof",vrec->ndof());
3080  //Fill(h,"fakeVtxNdof1",vrec->ndof());
3081  Fill(h,"fakeVtxNtrk",vrec->tracksSize());
3082  if(vrec->tracksSize()==2){ Fill(h,"fake2trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
3083  if(vrec->tracksSize()==3){ Fill(h,"fake3trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
3084  if(vrec->tracksSize()==4){ Fill(h,"fake4trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
3085  if(vrec->tracksSize()==5){ Fill(h,"fake5trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
3086  }
3087  }
3088 
3089 
3090  Fill(h,"nsimtrk",float(nsimtrk));
3091  Fill(h,"nsimtrk",float(nsimtrk),vsim==simpv.begin());
3092  Fill(h,"nrecsimtrk",float(vsim->nMatchedTracks));
3093  Fill(h,"nrecnosimtrk",float(nsimtrk-vsim->nMatchedTracks));
3094 
3095  // histogram properties of matched vertices
3096  if (vsim->recVtx && ( fabs(vsim->recVtx->z()-vsim->z*simUnit_)<zmatch_ )){
3097 
3098  if(verbose_){std::cout <<"primary matched " << message << " " << setw(8) << setprecision(4) << vsim->x << " " << vsim->y << " " << vsim->z << std:: endl;}
3099  Fill(h,"matchedVtxNdof", vsim->recVtx->ndof());
3100  // residuals an pulls with respect to simulated vertex
3101  Fill(h,"resx", vsim->recVtx->x()-vsim->x*simUnit_ );
3102  Fill(h,"resy", vsim->recVtx->y()-vsim->y*simUnit_ );
3103  Fill(h,"resz", vsim->recVtx->z()-vsim->z*simUnit_ );
3104  Fill(h,"resz10", vsim->recVtx->z()-vsim->z*simUnit_ );
3105  Fill(h,"pullx", (vsim->recVtx->x()-vsim->x*simUnit_)/vsim->recVtx->xError() );
3106  Fill(h,"pully", (vsim->recVtx->y()-vsim->y*simUnit_)/vsim->recVtx->yError() );
3107  Fill(h,"pullz", (vsim->recVtx->z()-vsim->z*simUnit_)/vsim->recVtx->zError() );
3108  Fill(h,"resxr", vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx);
3109  Fill(h,"resyr", vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy );
3110  Fill(h,"reszr", vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz);
3111  Fill(h,"pullxr", (vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx)/vsim->recVtx->xError() );
3112  Fill(h,"pullyr", (vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy)/vsim->recVtx->yError() );
3113  Fill(h,"pullzr", (vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz)/vsim->recVtx->zError() );
3114 
3115 
3116  // efficiency with zmatch within 500 um (or whatever zmatch is)
3117  Fill(h,"eff", 1.);
3118  if(simpv.size()==1){
3119  if (vsim->recVtx==&(*recVtxs->begin())){
3120  Fill(h,"efftag", 1.);
3121  }else{
3122  Fill(h,"efftag", 0.);
3123  cout << "signal vertex not tagged " << message << " " << eventcounter_ << endl;
3124  // call XXClusterizerInZ.vertices(seltrks,3)
3125  }
3126  }
3127 
3128  Fill(h,"effvsptsq",vsim->ptsq,1.);
3129  Fill(h,"effvsnsimtrk",vsim->nGenTrk,1.);
3130  Fill(h,"effvsnrectrk",nrectrks,1.);
3131  Fill(h,"effvsnseltrk",nseltrks,1.);
3132  Fill(h,"effvsz",vsim->z*simUnit_,1.);
3133  Fill(h,"effvsz2",vsim->z*simUnit_,1.);
3134  Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,1.);
3135 
3136  }else{ // no matching rec vertex found for this simvertex
3137 
3138  bool plapper=verbose_ && vsim->nGenTrk;
3139  if(plapper){
3140  // be quiet about invisble vertices
3141  std::cout << "primary not found " << message << " " << eventcounter_ << " x=" <<vsim->x*simUnit_ << " y=" << vsim->y*simUnit_ << " z=" << vsim->z*simUnit_ << " nGenTrk=" << vsim->nGenTrk << std::endl;
3142  }
3143  int mistype=0;
3144  if (vsim->recVtx){
3145  if(plapper){
3146  std::cout << "nearest recvertex at " << vsim->recVtx->z() << " dz=" << vsim->recVtx->z()-vsim->z*simUnit_ << std::endl;
3147  }
3148 
3149  if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.2 ){
3150  Fill(h,"effvsz2",vsim->z*simUnit_,1.);
3151  }
3152 
3153  if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.5 ){
3154  if(verbose_){std::cout << "type 1, lousy z vertex" << std::endl;}
3155  Fill(h,"zlost1", vsim->z*simUnit_,1.);
3156  mistype=1;
3157  }else{
3158  if(plapper){std::cout << "type 2a no vertex anywhere near" << std::endl;}
3159  mistype=2;
3160  }
3161  }else{// no recVtx at all
3162  mistype=2;
3163  if(plapper){std::cout << "type 2b, no vertex at all" << std::endl;}
3164  }
3165 
3166  if(mistype==2){
3167  int selstat=-3;
3168  // no matching vertex found, is there a cluster?
3169  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3170  if( fabs(clusters[iclu].position().z()-vsim->z*simUnit_) < 0.1 ){
3171  selstat=int(clusters[iclu].position().x()+0.1);
3172  if(verbose_){std::cout << "matching cluster found with selstat=" << clusters[iclu].position().x() << std::endl;}
3173  }
3174  }
3175  if (selstat==0){
3176  if(plapper){std::cout << "vertex rejected (distance to beam)" << std::endl;}
3177  Fill(h,"zlost3", vsim->z*simUnit_,1.);
3178  }else if(selstat==-1){
3179  if(plapper) {std::cout << "vertex invalid" << std::endl;}
3180  Fill(h,"zlost4", vsim->z*simUnit_,1.);
3181  }else if(selstat==1){
3182  if(plapper){std::cout << "vertex accepted, this cannot be right!!!!!!!!!!" << std::endl;}
3183  }else if(selstat==-2){
3184  if(plapper){std::cout << "dont know what this means !!!!!!!!!!" << std::endl;}
3185  }else if(selstat==-3){
3186  if(plapper){std::cout << "no matching cluster found " << std::endl;}
3187  Fill(h,"zlost2", vsim->z*simUnit_,1.);
3188  }else{
3189  if(plapper){std::cout << "dont know what this means either !!!!!!!!!!" << selstat << std::endl;}
3190  }
3191  }//mistype==2
3192 
3193 
3194  Fill(h,"eff", 0.);
3195  if(simpv.size()==1){ Fill(h,"efftag", 0.); }
3196 
3197  Fill(h,"effvsptsq",vsim->ptsq,0.);
3198  Fill(h,"effvsnsimtrk",float(vsim->nGenTrk),0.);
3199  Fill(h,"effvsnrectrk",nrectrks,0.);
3200  Fill(h,"effvsnseltrk",nseltrks,0.);
3201  Fill(h,"effvsz",vsim->z*simUnit_,0.);
3202  Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,0.);
3203 
3204  //part of the efficiency vs separation
3205 
3206 
3207  } // no recvertex for this simvertex
3208 
3209  }// end of sim/rec matching
3210 
3211 
3212  // purity of event vertex tags
3213  if (recVtxs->size()>0){
3214  Double_t dz=(*recVtxs->begin()).z() - (*simpv.begin()).z*simUnit_;
3215  Fill(h,"zdistancetag",dz);
3216  Fill(h,"abszdistancetag",fabs(dz));
3217  if( fabs(dz)<zmatch_){
3218  Fill(h,"puritytag",1.);
3219  }else{
3220  // bad tag: the true primary was more than 500 um (or zmatch) away from the tagged primary
3221  Fill(h,"puritytag",0.);
3222  }
3223  }
3224 
3225  }else{
3226  //if(verbose_) cout << "PrimaryVertexAnalyzer4PU::analyzeVertexCollection: simPV is empty!" << endl;
3227  }
3228 
3229 
3230  //******* the following code does not require MC and will/should work for data **********
3231 
3232 
3233  Fill(h,"bunchCrossing",bunchCrossing_);
3234  if(recTrks->size()>0) Fill(h,"bunchCrossingLogNtk",bunchCrossing_,log(recTrks->size())/log(10.));
3235 
3236  // ----------------- reconstructed tracks ------------------------
3237  // the list of selected tracks can only be correct if the selection parameterset and track collection
3238  // is the same that was used for the reconstruction
3239 
3240  int nt=0;
3241  for(reco::TrackCollection::const_iterator t=recTrks->begin();
3242  t!=recTrks->end(); ++t){
3243  if((recVtxs->size()>0) && (recVtxs->begin()->isValid())){
3244  fillTrackHistos(h,"all",*t,&(*recVtxs->begin()));
3245  }else{
3246  fillTrackHistos(h,"all",*t);
3247  }
3248  if(recTrks->size()>100) fillTrackHistos(h,"M",*t);
3249 
3250 
3251 
3252  TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
3253  if (theTrackFilter(tt)){
3254  selTrks.push_back(*t);
3255  fillTrackHistos(h,"sel",*t);
3256  int foundinvtx=0;
3257  int nvtemp=-1;
3258  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3259  v!=recVtxs->end(); ++v){
3260  nvtemp++;
3261  if(( v->isFake()) || (v->ndof()<-2) ) break;
3262  for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++ ){
3263  if( ((**tv).vz()==t->vz()&&((**tv).phi()==t->phi())) ) {
3264  foundinvtx++;
3265  }
3266  }
3267 
3268  }
3269  if(foundinvtx==0){
3270  fillTrackHistos(h,"sellost",*t);
3271  }else if(foundinvtx>1){
3272  cout << "hmmmm " << foundinvtx << endl;
3273  }
3274  }
3275  nt++;
3276  }
3277 
3278 
3279  if (nseltrks<0){
3280  nseltrks=selTrks.size();
3281  }else if( ! (nseltrks==(int)selTrks.size()) ){
3282  std::cout << "Warning: inconsistent track selection !" << std::endl;
3283  }
3284 
3285 
3286 
3287  // fill track histograms of vertex tracks
3288  int nrec=0, nrec0=0, nrec8=0, nrec2=0, nrec4=0;
3289  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3290  v!=recVtxs->end(); ++v){
3291 
3292  if (! (v->isFake()) && v->ndof()>0 && v->chi2()>0 ){
3293  nrec++;
3294  if (v->ndof()>0) nrec0++;
3295  if (v->ndof()>8) nrec8++;
3296  if (v->ndof()>2) nrec2++;
3297  if (v->ndof()>4) nrec4++;
3298  for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){
3299  if(v==recVtxs->begin()){
3300  fillTrackHistos(h,"tagged",**t, &(*v));
3301  }else{
3302  fillTrackHistos(h,"untagged",**t, &(*v));
3303  }
3304 
3305  Float_t wt=v->trackWeight(*t);
3306  //dumpHitInfo(**t); cout << " w=" << wt << endl;
3307  Fill(h,"trackWt",wt);
3308  if(wt>0.5){
3309  fillTrackHistos(h,"wgt05",**t, &(*v));
3310  if(v->ndof()>4) fillTrackHistos(h,"ndof4",**t, &(*v));
3311  }else{
3312  fillTrackHistos(h,"wlt05",**t, &(*v));
3313  }
3314  }
3315  }
3316  }
3317 
3318 
3319  // bachelor tracks (only available through clusters right now)
3320  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3321  if (clusters[iclu].tracksSize()==1){
3322  for(trackit_t t = clusters[iclu].tracks_begin();
3323  t!=clusters[iclu].tracks_end(); t++){
3324  fillTrackHistos(h,"bachelor",**t);
3325  }
3326  }
3327  }
3328 
3329 
3330  // ----------------- reconstructed vertices ------------------------
3331 
3332  // event
3333  Fill(h,"szRecVtx",recVtxs->size());
3334  Fill(h,"nclu",clusters.size());
3335  Fill(h,"nseltrk",nseltrks);
3336  Fill(h,"nrectrk",nrectrks);
3337  Fill(h,"nrecvtx",nrec);
3338  Fill(h,"nrecvtx2",nrec2);
3339  Fill(h,"nrecvtx4",nrec4);
3340  Fill(h,"nrecvtx8",nrec8);
3341 
3342  if(nrec>0){
3343  Fill(h,"eff0vsntrec",nrectrks,1.);
3344  Fill(h,"eff0vsntsel",nseltrks,1.);
3345  }else{
3346  Fill(h,"eff0vsntrec",nrectrks,0.);
3347  Fill(h,"eff0vsntsel",nseltrks,0.);
3348  if((nseltrks>1)&&(verbose_)){
3349  cout << Form("PrimaryVertexAnalyzer4PU: %s may have lost a vertex %10d %10d %4d / %4d ",message.c_str(),run_, event_, nrectrks,nseltrks) << endl;
3350  dumpThisEvent_=true;
3351  }
3352  }
3353  if(nrec0>0) { Fill(h,"eff0ndof0vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof0vsntsel",nseltrks,0.);}
3354  if(nrec2>0) { Fill(h,"eff0ndof2vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof2vsntsel",nseltrks,0.);}
3355  if(nrec4>0) { Fill(h,"eff0ndof4vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof4vsntsel",nseltrks,0.);}
3356  if(nrec8>0) { Fill(h,"eff0ndof8vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof8vsntsel",nseltrks,0.);}
3357 
3358  if((nrec>1)&&(DEBUG_)) {
3359  cout << "multivertex event" << endl;
3360  dumpThisEvent_=true;
3361  }
3362 
3363  if((nrectrks>10)&&(nseltrks<3)){
3364  cout << "small fraction of selected tracks " << endl;
3365  dumpThisEvent_=true;
3366  }
3367 
3368  // properties of events without a vertex
3369  if((nrec==0)||(recVtxs->begin()->isFake())){
3370  Fill(h,"nrectrk0vtx",nrectrks);
3371  Fill(h,"nseltrk0vtx",nseltrks);
3372  Fill(h,"nclu0vtx",clusters.size());
3373  }
3374 
3375 
3376  // properties of (valid) vertices
3377  double ndof2=-10,ndof1=-10, zndof1=0, zndof2=0;
3378  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
3379  v!=recVtxs->end(); ++v){
3380  if(v->isFake()){ Fill(h,"isFake",1.);}else{ Fill(h,"isFake",0.);}
3381  if(v->isFake()||((v->ndof()<0)&&(v->ndof()>-3))){ Fill(h,"isFake1",1.);}else{ Fill(h,"isFake1",0.);}
3382 
3383  if((v->isFake())||(v->ndof()<0)) continue;
3384  if(v->ndof()>ndof1){ ndof2=ndof1; zndof2=zndof1; ndof1=v->ndof(); zndof1=v->position().z();}
3385  else if(v->ndof()>ndof2){ ndof2=v->ndof(); zndof2=v->position().z();}
3386 
3387 
3388  // some special histogram for two track vertices
3389  if(v->tracksSize()==2){
3390  const TrackBaseRef& t1= *(v->tracks_begin());
3391  const TrackBaseRef& t2=*(v->tracks_begin()+1);
3392  bool os=(t1->charge()*t2->charge()<0);
3393  double dphi=t1->phi()-t2->phi(); if (dphi<0) dphi+=2*M_PI;
3394  double m12=sqrt(pow( sqrt(pow(0.139,2)+pow( t1->p(),2)) +sqrt(pow(0.139,2)+pow( t2->p(),2)) ,2)
3395  -pow(t1->px()+t2->px(),2)
3396  -pow(t1->py()+t2->py(),2)
3397  -pow(t1->pz()+t2->pz(),2)
3398  );
3399  if(os){
3400  Fill(h,"2trkdetaOS",t1->eta()-t2->eta());
3401  Fill(h,"2trkmassOS",m12);
3402  }else{
3403  Fill(h,"2trkdetaSS",t1->eta()-t2->eta());
3404  Fill(h,"2trkmassSS",m12);
3405  }
3406  Fill(h,"2trkdphi",dphi);
3407  Fill(h,"2trkseta",t1->eta()+t2->eta());
3408  if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurl",t1->eta()+t2->eta());
3409  if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurl",dphi);
3410  // fill separately for extra vertices
3411  if(v!=recVtxs->begin()){
3412  if(os){
3413  Fill(h,"2trkdetaOSPU",t1->eta()-t2->eta());
3414  Fill(h,"2trkmassOSPU",m12);
3415  }else{
3416  Fill(h,"2trkdetaSSPU",t1->eta()-t2->eta());
3417  Fill(h,"2trkmassSSPU",m12);
3418  }
3419  Fill(h,"2trkdphiPU",dphi);
3420  Fill(h,"2trksetaPU",t1->eta()+t2->eta());
3421  if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurlPU",t1->eta()+t2->eta());
3422  if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurlPU",dphi);
3423  }
3424  }
3425 
3426 
3427  Fill(h,"trkchi2vsndof",v->ndof(),v->chi2());
3428  if(v->ndof()>0){ Fill(h,"trkchi2overndof",v->chi2()/v->ndof()); }
3429  if(v->tracksSize()==2){ Fill(h,"2trkchi2vsndof",v->ndof(),v->chi2()); }
3430  if(v->tracksSize()==3){ Fill(h,"3trkchi2vsndof",v->ndof(),v->chi2()); }
3431  if(v->tracksSize()==4){ Fill(h,"4trkchi2vsndof",v->ndof(),v->chi2()); }
3432  if(v->tracksSize()==5){ Fill(h,"5trkchi2vsndof",v->ndof(),v->chi2()); }
3433 
3434  Fill(h,"nbtksinvtx",v->tracksSize());
3435  Fill(h,"nbtksinvtx2",v->tracksSize());
3436  Fill(h,"vtxchi2",v->chi2());
3437  Fill(h,"vtxndf",v->ndof());
3438  Fill(h,"vtxprob",ChiSquaredProbability(v->chi2() ,v->ndof()));
3439  Fill(h,"vtxndfvsntk",v->tracksSize(), v->ndof());
3440  Fill(h,"vtxndfoverntk",v->ndof()/v->tracksSize());
3441  Fill(h,"vtxndf2overntk",(v->ndof()+2)/v->tracksSize());
3442  Fill(h,"zrecvsnt",v->position().z(),float(nrectrks));
3443  if(nrectrks>100){
3444  Fill(h,"zrecNt100",v->position().z());
3445  }
3446 
3447  if(v->ndof()>2.0){ // enter only vertices that really contain tracks
3448  Fill(h,"xrec",v->position().x());
3449  Fill(h,"yrec",v->position().y());
3450  Fill(h,"zrec",v->position().z());
3451  Fill(h,"xrec1",v->position().x());
3452  Fill(h,"yrec1",v->position().y());
3453  Fill(h,"zrec1",v->position().z());
3454  Fill(h,"xrec2",v->position().x());
3455  Fill(h,"yrec2",v->position().y());
3456  Fill(h,"zrec2",v->position().z());
3457  Fill(h,"xrec3",v->position().x());
3458  Fill(h,"yrec3",v->position().y());
3459  Fill(h,"zrec3",v->position().z());
3460  Fill(h,"xrecb",v->position().x()-vertexBeamSpot_.x0());
3461  Fill(h,"yrecb",v->position().y()-vertexBeamSpot_.y0());
3462  Fill(h,"zrecb",v->position().z()-vertexBeamSpot_.z0());
3463  Fill(h,"xrecBeam",v->position().x()-vertexBeamSpot_.x0());
3464  Fill(h,"yrecBeam",v->position().y()-vertexBeamSpot_.y0());
3465  Fill(h,"zrecBeam",v->position().z()-vertexBeamSpot_.z0());
3466  Fill(h,"xrecBeamPull",(v->position().x()-vertexBeamSpot_.x0())/sqrt(pow(v->xError(),2)+pow(vertexBeamSpot_.BeamWidthX(),2)));
3467  Fill(h,"yrecBeamPull",(v->position().y()-vertexBeamSpot_.y0())/sqrt(pow(v->yError(),2)+pow(vertexBeamSpot_.BeamWidthY(),2)));
3468 
3469  Fill(h,"xrecBeamvsdx",v->xError(),v->position().x()-vertexBeamSpot_.x0());
3470  Fill(h,"yrecBeamvsdy",v->yError(),v->position().y()-vertexBeamSpot_.y0());
3471  Fill(h,"xrecBeamvsdxR2",v->position().x()-vertexBeamSpot_.x0(),v->xError());
3472  Fill(h,"yrecBeamvsdyR2",v->position().y()-vertexBeamSpot_.y0(),v->yError());
3473  Fill(h,"xrecBeam2vsdx2prof",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
3474  Fill(h,"yrecBeam2vsdy2prof",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
3475  Fill(h,"xrecBeamvsdx2",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
3476  Fill(h,"yrecBeamvsdy2",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
3477  Fill(h,"xrecBeamvsz",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
3478  Fill(h,"yrecBeamvsz",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
3479  Fill(h,"xrecBeamvszprof",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
3480  Fill(h,"yrecBeamvszprof",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
3481  Fill(h,"xrecBeamvsdxprof",v->xError(),v->position().x()-vertexBeamSpot_.x0());
3482  Fill(h,"yrecBeamvsdyprof",v->yError(),v->position().y()-vertexBeamSpot_.y0());
3483 
3484 
3485  Fill(h,"errx",v->xError());
3486  Fill(h,"erry",v->yError());
3487  Fill(h,"errz",v->zError());
3488  double vxx=v->covariance(0,0);
3489  double vyy=v->covariance(1,1);
3490  double vxy=v->covariance(1,0);
3491  double dv=0.25*(vxx+vyy)*(vxx+vyy-(vxx*vyy-vxy*vxy));
3492  if(dv>0){
3493  double l1=0.5*(vxx+vyy)+sqrt(dv);
3494  Fill(h,"err1",sqrt(l1));
3495  double l2=sqrt(0.5*(vxx+vyy)-sqrt(dv));
3496  if(l2>0) Fill(h,"err2",sqrt(l2));
3497  }
3498 
3499 
3500  // look at the tagged vertex separately
3501  if (v==recVtxs->begin()){
3502  Fill(h,"nbtksinvtxTag",v->tracksSize());
3503  Fill(h,"nbtksinvtxTag2",v->tracksSize());
3504  Fill(h,"xrectag",v->position().x());
3505  Fill(h,"yrectag",v->position().y());
3506  Fill(h,"zrectag",v->position().z());
3507  }else{
3508  Fill(h,"nbtksinvtxPU",v->tracksSize());
3509  Fill(h,"nbtksinvtxPU2",v->tracksSize());
3510  }
3511 
3512  // vertex resolution vs number of tracks
3513  Fill(h,"xresvsntrk",v->tracksSize(),v->xError());
3514  Fill(h,"yresvsntrk",v->tracksSize(),v->yError());
3515  Fill(h,"zresvsntrk",v->tracksSize(),v->zError());
3516 
3517  }
3518 
3519  // cluster properties (if available)
3520  for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
3521  if( fabs(clusters[iclu].position().z()-v->position().z()) < 0.0001 ){
3522  Fill(h,"nclutrkvtx",clusters[iclu].tracksSize());
3523  }
3524  }
3525 
3526 
3527 
3528  // properties of (valid) neighbour vertices
3529  reco::VertexCollection::const_iterator v1=v; v1++;
3530  for(; v1!=recVtxs->end(); ++v1){
3531  if((v1->isFake())||(v1->ndof()<0)) continue;
3532  Fill(h,"zdiffrec",v->position().z()-v1->position().z());
3533 // if(fabs(v->position().z()-v1->position().z())>1){
3534 // cout << message << " zdiffrec=" << v->position().z()-v1->position().z() << " " << v->ndof() << " " << v1->ndof() << endl;
3535 // dumpThisEvent_=true;
3536 // }
3537 
3538  double z0=v->position().z()-vertexBeamSpot_.z0();
3539  double z1=v1->position().z()-vertexBeamSpot_.z0();
3540  Fill(h,"zPUcand",z0); Fill(h,"zPUcand",z1);
3541  Fill(h,"ndofPUcand",v->ndof()); Fill(h,"ndofPUcand",v1->ndof());
3542 
3543  Fill(h,"zdiffvsz",z1-z0,0.5*(z1+z0));
3544 
3545  if ((v->ndof()>2) && (v1->ndof()>2)){
3546  Fill(h,"zdiffrec2",v->position().z()-v1->position().z());
3547  Fill(h,"zPUcand2",z0);
3548  Fill(h,"zPUcand2",z1);
3549  Fill(h,"ndofPUcand2",v->ndof());
3550  Fill(h,"ndofPUcand2",v1->ndof());
3551  Fill(h,"zvszrec2",z0, z1);
3552  Fill(h,"pzvspz2",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
3553  }
3554 
3555  if ((v->ndof()>4) && (v1->ndof()>4)){
3556  Fill(h,"zdiffvsz4",z1-z0,0.5*(z1+z0));
3557  Fill(h,"zdiffrec4",v->position().z()-v1->position().z());
3558  Fill(h,"zvszrec4",z0, z1);
3559  Fill(h,"pzvspz4",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
3560  //cout << "ndof4 pu-candidate " << run_ << " " << event_ << endl ;
3561  if(fabs(z0-z1)>1.0){
3562  Fill(h,"xbeamPUcand",v->position().x()-vertexBeamSpot_.x0());
3563  Fill(h,"ybeamPUcand",v->position().y()-vertexBeamSpot_.y0());
3564  Fill(h,"xbeamPullPUcand",(v->position().x()-vertexBeamSpot_.x0())/v->xError());
3565  Fill(h,"ybeamPullPUcand",(v->position().y()-vertexBeamSpot_.y0())/v->yError());
3566  //Fill(h,"sumwOverNtkPUcand",sumw/v->tracksSize());
3567  //Fill(h,"sumwOverNtkPUcand",sumw/v1->tracksSize());
3568  Fill(h,"ndofOverNtkPUcand",v->ndof()/v->tracksSize());
3569  Fill(h,"ndofOverNtkPUcand",v1->ndof()/v1->tracksSize());
3570  Fill(h,"xbeamPUcand",v1->position().x()-vertexBeamSpot_.x0());
3571  Fill(h,"ybeamPUcand",v1->position().y()-vertexBeamSpot_.y0());
3572  Fill(h,"xbeamPullPUcand",(v1->position().x()-vertexBeamSpot_.x0())/v1->xError());
3573  Fill(h,"ybeamPullPUcand",(v1->position().y()-vertexBeamSpot_.y0())/v1->yError());
3574  Fill(h,"zPUcand4",z0);
3575  Fill(h,"zPUcand4",z1);
3576  Fill(h,"ndofPUcand4",v->ndof());
3577  Fill(h,"ndofPUcand4",v1->ndof());
3578  for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){ if(v->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v)); }
3579  for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v1)); }
3580  }
3581  }
3582 
3583  if ((v->ndof()>4) && (v1->ndof()>2) && (v1->ndof()<4)){
3584  for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUfake",**t, &(*v1)); }
3585  }
3586 
3587  if ((v->ndof()>8) && (v1->ndof()>8)){
3588  Fill(h,"zdiffrec8",v->position().z()-v1->position().z());
3589  if(dumpPUcandidates_ && fabs(z0-z1)>1.0){
3590  cout << "PU candidate " << run_ << " " << event_ << " " << message << " zdiff=" <<z0-z1 << endl;
3591 // cerr << "PU candidate " << run_ << " "<< event_ << " " << message << endl;
3592  dumpThisEvent_=true;
3593  }
3594  }
3595 
3596  }
3597 
3598  // test track links, use reconstructed vertices
3599  bool problem = false;
3600  Fill(h,"nans",1.,std::isnan(v->position().x())*1.);
3601  Fill(h,"nans",2.,std::isnan(v->position().y())*1.);
3602  Fill(h,"nans",3.,std::isnan(v->position().z())*1.);
3603 
3604  int index = 3;
3605  for (int i = 0; i != 3; i++) {
3606  for (int j = i; j != 3; j++) {
3607  index++;
3608  Fill(h,"nans",index*1., std::isnan(v->covariance(i, j))*1.);
3609  if (std::isnan(v->covariance(i, j))) problem = true;
3610  // in addition, diagonal element must be positive
3611  if (j == i && v->covariance(i, j) < 0) {
3612  Fill(h,"nans",index*1., 1.);
3613  problem = true;
3614  }
3615  }
3616  }
3617 
3618  try {
3619  for(trackit_t t = v->tracks_begin();
3620  t!=v->tracks_end(); t++) {
3621  // illegal charge
3622  if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
3623  Fill(h,"tklinks",0.);
3624  }
3625  else {
3626  Fill(h,"tklinks",1.);
3627  }
3628  }
3629  } catch (...) {
3630  // exception thrown when trying to use linked track
3631  Fill(h,"tklinks",0.);
3632  }
3633 
3634  if (problem) {
3635  // analyze track parameter covariance definiteness
3636  double data[25];
3637  try {
3638  int itk = 0;
3639  for(trackit_t t = v->tracks_begin();
3640  t!=v->tracks_end(); t++) {
3641  std::cout << "Track " << itk++ << std::endl;
3642  int i2 = 0;
3643  for (int i = 0; i != 5; i++) {
3644  for (int j = 0; j != 5; j++) {
3645  data[i2] = (**t).covariance(i, j);
3646  std::cout << std:: scientific << data[i2] << " ";
3647  i2++;
3648  }
3649  std::cout << std::endl;
3650  }
3651  gsl_matrix_view m
3652  = gsl_matrix_view_array (data, 5, 5);
3653 
3654  gsl_vector *eval = gsl_vector_alloc (5);
3655  gsl_matrix *evec = gsl_matrix_alloc (5, 5);
3656 
3657  gsl_eigen_symmv_workspace * w =
3658  gsl_eigen_symmv_alloc (5);
3659 
3660  gsl_eigen_symmv (&m.matrix, eval, evec, w);
3661 
3662  gsl_eigen_symmv_free (w);
3663 
3664  gsl_eigen_symmv_sort (eval, evec,
3665  GSL_EIGEN_SORT_ABS_ASC);
3666 
3667  // print sorted eigenvalues
3668  {
3669  int i;
3670  for (i = 0; i < 5; i++) {
3671  double eval_i
3672  = gsl_vector_get (eval, i);
3673  //gsl_vector_view evec_i
3674  // = gsl_matrix_column (evec, i);
3675 
3676  printf ("eigenvalue = %g\n", eval_i);
3677  // printf ("eigenvector = \n");
3678  // gsl_vector_fprintf (stdout,
3679  // &evec_i.vector, "%g");
3680  }
3681  }
3682  }
3683  }
3684  catch (...) {
3685  // exception thrown when trying to use linked track
3686  break;
3687  }// catch()
3688  }// if (problem)
3689 
3690 
3691 
3692  } // vertex loop (v)
3693 
3694 
3695  // 2nd highest ndof
3696  if (ndof2>0){
3697  Fill(h,"ndofnr2",ndof2);
3698  if(fabs(zndof1-zndof2)>1) Fill(h,"ndofnr2d1cm",ndof2);
3699  if(fabs(zndof1-zndof2)>2) Fill(h,"ndofnr2d2cm",ndof2);
3700  }
3701 
3702  //part of the separation efficiency profiles
3703 
3704  if ( pui_z.size()>0 ) {
3705 
3706  vector<int> used_indizesV;
3707 
3708  for (unsigned spv_idx=0; spv_idx<pui_z.size(); spv_idx++) {
3709 
3710  float sv = pui_z[spv_idx];
3711 
3712  if ( fabs(sv)>24. ) continue;
3713 
3714  double eff = 0.;
3715  double effwod = 0.;
3716  double dreco = 0.;
3717  double dsimed = 0.;
3718 
3719  vector<int>* matchedV = vertex_match(sv, recVtxs);
3720  unsigned numMatchs = matchedV->size();
3721 
3722  bool dsflag = false;
3723 
3724  for (unsigned i=0;i<used_indizesV.size(); i++) {
3725  for (unsigned j=0;j<numMatchs; j++) {
3726  if ( used_indizesV.at(i)==matchedV->at(j) ) {
3727  dsflag = true;
3728  break;
3729  }
3730  }
3731  }
3732 
3733  if ( numMatchs>0 ) eff = 1.;
3734  if ( numMatchs>1 ) dreco = 1.;
3735  if ( dsflag ) dsimed = 1.;
3736  if ( ( numMatchs>0 ) && ( !dsflag ) ) effwod = 1.;
3737 
3738  for (unsigned i=0;i<numMatchs; i++) {
3739  used_indizesV.push_back(matchedV->at(i));
3740  }
3741 
3742  float sep = getTrueSeparation(sv, pui_z);
3743 
3744  Fill(h,"effvszsep", sep, eff);
3745  Fill(h,"effwodvszsep", sep, effwod);
3746  Fill(h,"mergedvszsep", sep, dsimed);
3747  Fill(h,"splitvszsep", sep, dreco);
3748 
3749  }
3750 
3751  } else {
3752  std::cout << "Strange PileUpSummaryInfo in the event" << std::endl;
3753  }
3754 
3755 
3756 }
3757 
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
reco::TrackBase::qoverp
double qoverp() const
q / p
Definition: TrackBase.h:549
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:591
type
type
Definition: HCALResponse.h:21
PrimaryVertexAnalyzer4PU::theTrackFilter
TrackFilterForPVFinding theTrackFilter
Definition: PrimaryVertexAnalyzer4PU.h:341
edm::RefToBase::get
value_type const * get() const
Definition: RefToBase.h:213
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
edm::View::size_type
unsigned int size_type
Definition: View.h:85
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:778
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
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
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:199
PrimaryVertexAnalyzer4PU::getSimEvents
std::vector< PrimaryVertexAnalyzer4PU::SimEvent > getSimEvents(edm::Handle< TrackingParticleCollection >, edm::Handle< TrackingVertexCollection >, edm::Handle< edm::View< reco::Track > >)
Definition: PrimaryVertexAnalyzer4PU.cc:1635
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
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:537
PrimaryVertexAnalyzer4PU::ParameterVector
reco::TrackBase::ParameterVector ParameterVector
Definition: PrimaryVertexAnalyzer4PU.h:71
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
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:555
reco::TrackBase::dxyError
double dxyError() const
error on dxy
Definition: TrackBase.h:772
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:13
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:917
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
edm::View::size
size_type size() const
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:621
NULL
#define NULL
Definition: scimark2.h:8
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
reco::TrackBase::ParameterVector
math::Vector< dimension >::type ParameterVector
parameter vector
Definition: TrackBase.h:73
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:603
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:2986
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:78
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:55
reco::TrackBase::algo
TrackAlgorithm algo() const
Definition: TrackBase.h:403
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
PrimaryVertexAnalyzer4PU::getTrueSeparation
double getTrueSeparation(float, const std::vector< float > &)
Definition: PrimaryVertexAnalyzer4PU.cc:1882
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::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:627
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:1581
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:912
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:597
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:1349
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:832
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
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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:561
PrimaryVertexAnalyzer4PU::getSimPVs
std::vector< simPrimaryVertex > getSimPVs(const edm::Handle< edm::HepMCProduct > evtMC)
Definition: PrimaryVertexAnalyzer4PU.cc:1767
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
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:1939
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::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
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:104
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:2371
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:615
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:790
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:645
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:1484
Measurement1D::significance
double significance() const
Definition: Measurement1D.h:32
PrimaryVertexAnalyzer4PU::isResonance
bool isResonance(const HepMC::GenParticle *p)
Definition: PrimaryVertexAnalyzer4PU.cc:1027
edm::Handle::product
T const * product() const
Definition: Handle.h:81
reco::TransientTrack::track
const Track & track() const
Definition: TransientTrack.h:114
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:364
SiStripMonitorClusterAlca_cfi.zmin
tuple zmin
Definition: SiStripMonitorClusterAlca_cfi.py:182
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:1530
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:440
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
b
double b
Definition: hdecay.h:120
Measurement1D::value
double value() const
Definition: Measurement1D.h:28
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:639
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:194
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:2596
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:543
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:62
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:72
w
T w() const
Definition: extBasic3DVector.h:245
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
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:567
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
PrimaryVertexAnalyzer4PU::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PrimaryVertexAnalyzer4PU.cc:2050
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:1611
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:609
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:633
reco::TrackBase::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:76
SiStripMonitorClusterAlca_cfi.zmax
tuple zmax
Definition: SiStripMonitorClusterAlca_cfi.py:183
reco::HitPattern::numberOfHits
int numberOfHits(HitCategory category) const
Definition: HitPattern.h:721
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:109
Generated for CMSSW Reference Manual by doxygen 1.8.5