ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc
(Generate patch)

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.13 by jengbou, Fri Dec 4 08:37:40 2009 UTC vs.
Revision 1.22 by yygao, Mon Mar 29 15:27:24 2010 UTC

# Line 12 | Line 12 | Implementation:
12   */
13   //
14   // Original Author:  "Geng-yuan Jeng/UC Riverside"
15 < //                   "Yanyan Gao/Fermilab ygao@fnal.gov"
15 > //                    "Yanyan Gao/Fermilab ygao@fnal.gov"
16   //         Created:  Thu Aug 20 11:55:40 CDT 2009
17   // $Id$
18   //
# Line 25 | Line 25 | Implementation:
25   #include <vector>
26   #include <iostream>
27   #include <sstream>
28 – #include <algorithm>
28  
29   // user include files
30   #include "FWCore/Framework/interface/Frameworkfwd.h"
# Line 39 | Line 38 | Implementation:
38   #include "DataFormats/TrackReco/interface/Track.h"
39   #include "DataFormats/TrackReco/interface/TrackFwd.h"
40   #include "FWCore/ServiceRegistry/interface/Service.h"
41 < #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
41 > #include "CommonTools/UtilAlgos/interface/TFileService.h"
42   #include "UserCode/PVStudy/interface/PVStudy.h"
43   //
44   #include "DataFormats/VertexReco/interface/Vertex.h"
# Line 50 | Line 49 | Implementation:
49   #include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
50   #include <SimDataFormats/Track/interface/SimTrack.h>
51   #include <SimDataFormats/Track/interface/SimTrackContainer.h>
52 + // BeamSpot
53 + #include "DataFormats/BeamSpot/interface/BeamSpot.h"
54  
55   //root
56   #include <TROOT.h>
# Line 60 | Line 61 | Implementation:
61   #include <TPad.h>
62  
63   using namespace std;
64 + typedef math::XYZTLorentzVectorF LorentzVector;
65 + typedef math::XYZPoint Point;
66  
67   PVStudy::PVStudy(const edm::ParameterSet& iConfig)
68   {
# Line 72 | Line 75 | PVStudy::PVStudy(const edm::ParameterSet
75    splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
76    vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
77    splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
78 <  splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");  
76 <  pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");  
78 >  splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
79    verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
80    realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
81    analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
82    saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);  
83 <  outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
83 >  outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");  
84 >  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
85    ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83 –  zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
86    nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87    nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88 <  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
89 <
88 >  zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89 >  fHighPuritycut_            = iConfig.getUntrackedParameter<double>("fHighPuritycut");  
90 >  useHWTrk_                  = iConfig.getUntrackedParameter<bool>("useHWTrk",false);
91 >  bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
92 >  
93    //now do what ever initialization is needed
94    pvinfo.clear();
95  
91 –  
96    // Specify the data mode vector
97    if(realData_) datamode.push_back(0);
98    else {
# Line 97 | Line 101 | PVStudy::PVStudy(const edm::ParameterSet
101      datamode.push_back(2);
102      datamode.push_back(3);
103    }
104 <
101 <  // Create ntuple files if needed
104 > // Create ntuple files if needed
105    if(saventuple_) {
106      file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107      ftree_ = new TTree("mytree","mytree");
# Line 109 | Line 112 | PVStudy::PVStudy(const edm::ParameterSet
112  
113      // pvtxtree_ analyzing the pvtxs ootb
114      pvtxtree_ = new TTree("pvtxtree","pvtxtree");
112 –    //pvtx using all tracks
113 –    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
114 –    pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
115 –    pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");    
116 –    pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
117 –    pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
118 –    pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
119 –    pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
120 –    pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
121 –    pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
122 –    pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
123 –    pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
115  
116 +    // Event information for the data
117 +    pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
118 +    pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
119 +    pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
120 +    pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");  
121 +    
122 +
123 +    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
124 +    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
125 +    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
126 +    // Event level information
127      pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
128      pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
129      pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
130  
131 <    //pvtx using splittracks1
132 <    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
133 <    pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");  
134 <    pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");  
133 <    pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
134 <    pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
135 <    pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
136 <    pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
137 <    pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
138 <    pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
139 <    pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
140 <    pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");  
141 <  
142 <    //pvtx using splittracks2
143 <    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
144 <    pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");    
145 <    pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");  
146 <    pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
147 <    pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
148 <    pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
149 <    pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
150 <    pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
151 <    pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
152 <    pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
153 <    pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");  
154 <    
155 <    //pixeVertices
156 <    pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
157 <    pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");    
158 <    pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");  
159 <    pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
160 <    pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
161 <    pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
162 <    pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
163 <    pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
164 <    pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165 <    pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166 <    pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");  
131 >    //Fill the variables in the twovtx pair (recvtx1, recvtx2)  
132 >    // Information related to the analyzing the two-vertex method
133 >
134 >
135  
168 –    //Fill the variables in the twovtx pair (recvtx1, recvtx2)
136      pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
170 –    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
171 –    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
172 –    pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
173 –    pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
137      pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
138      pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
139      pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
# Line 181 | Line 144 | PVStudy::PVStudy(const edm::ParameterSet
144      pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
145      pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
146      pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
147 +
148 +    // Information for the splitVertexColl1
149 +    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
150 +    pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
151 +    pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
152 +    pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
153 +    pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
154 +    pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
155 +    pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
156 +  
157 +    
158 +    // Information for the splitVertexColl2
159 +    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
160 +    pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
161 +    pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
162 +    pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
163 +    pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
164 +    pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
165 +    pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
166 +  
167      
168      // MC variables
169      if(!realData_) {  
# Line 198 | Line 181 | PVStudy::PVStudy(const edm::ParameterSet
181        pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
182        pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
183        pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
184 <      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
184 >      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");    
185 >      pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
186 >      pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
187 >      pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
188 >      pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
189 >      pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
190 >      pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
191 >      pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
192 >
193        // For pvtxs with splittracks1
194        pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
195        pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
# Line 207 | Line 198 | PVStudy::PVStudy(const edm::ParameterSet
198        pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
199        pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
200        pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
201 +      pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
202 +      pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
203 +      pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
204 +      pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
205 +      pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
206 +      pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
207 +      pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
208 +        
209        // For pvtxs with splittracks1
210        pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
211        pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
# Line 214 | Line 213 | PVStudy::PVStudy(const edm::ParameterSet
213        pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
214        pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
215        pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
216 <      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
216 >      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");  
217 >      pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
218 >      pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
219 >      pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");  
220 >      pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
221 >      pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
222 >      pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
223 >      pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
224 >
225      }
226    }
227 <
227 >  
228    setRootStyle();
229  
230 +  //========================================
231 +  // Booking histograms
232 +  //========================================
233 +  
234    // Create a root file for histograms
235    theFile = new TFile(histoFileName_.c_str(), "RECREATE");
236    // make diretories
# Line 227 | Line 238 | PVStudy::PVStudy(const edm::ParameterSet
238    theFile->cd();
239    theFile->mkdir("Others");
240    theFile->cd();
241 +  if (analyze_) {
242 +    theFile->mkdir("Results");
243 +    theFile->cd();
244 +  }
245  
246    // Book Histograms:
247    h_pvtrk   = new PVHistograms();
233 –  h_pixvtx  = new PVHistograms();
248    h_misc    = new PVHistograms();
249    h_summary = new PVHistograms();
250    h_others  = new PVHistograms();
# Line 240 | Line 254 | PVStudy::PVStudy(const edm::ParameterSet
254      else {
255        stringstream ss;
256        ss << i;
257 <      h_pvtrk->Init("pvTrk",ss.str(),"spl");
257 >      h_pvtrk->Init("pvTrk", ss.str(),"spl");
258      }
259    }
246 –  h_pixvtx->Init("pixVtx");
260    h_misc->Init("misc");
261  
262    // Book MC only plots
# Line 251 | Line 264 | PVStudy::PVStudy(const edm::ParameterSet
264      h_gen = new PVHistograms();
265      h_gen->Init("generator");
266    }
267 +  if (analyze_) h_ana = new PVHistograms();
268  
269    //Book histograms sensitive to data/mc
270    for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
271      string suffix;
272 <    if(verbose_)  cout<<"datamode = "<< *it<<endl;
272 >    edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
273      switch(*it) {
274      case 0: suffix = "";
275        break;
# Line 276 | Line 290 | PVStudy::PVStudy(const edm::ParameterSet
290      }
291  
292      // Book residual and pull histograms only when analyzeOntheFly is enabled
293 <    if(analyze_) {
280 <      h_ana = new PVHistograms();
281 <      h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
282 <    }
293 >    if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
294      suffix.clear();
295    } // End of Book histograms sensitive to data/mc    
296 +  
297  
298   }
299  
300   PVStudy::~PVStudy()
301   {
302 +
303 +  // do anything here that needs to be done at desctruction time
304 +  // (e.g. close files, deallocate resources etc.)
305    theFile->cd();
306    theFile->cd("Summary");
307    h_pvtrk->Save();
293 –  h_pixvtx->Save();
294 –  h_gen->Save();
308    h_misc->Save();
309    h_summary->Save();
310 <  if (analyze_)
310 >  if (!realData_)
311 >    h_gen->Save();
312 >  if (analyze_) {
313 >    theFile->cd();
314 >    theFile->cd("Results");
315      h_ana->Save();
316 +  }
317    theFile->cd();
318    theFile->cd("Others");
319    h_others->Save();
# Line 324 | Line 342 | void PVStudy::setRootStyle() {
342    gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
343    gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
344    gStyle->SetStatFontSize(0.08);      // size for stat. box
345 <  gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
346 <  gStyle->SetTitleFont(32, "");       // same for pad title
347 <  gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
348 <  gStyle->SetStatFont(32);            // same for stat. box
345 >  gStyle->SetTitleFont(42, "XYZ");    // times-bold-italic font (p. 153) for axes
346 >  gStyle->SetTitleFont(42, "");       // same for pad title
347 >  gStyle->SetLabelFont(42, "XYZ");    // same for axis labels
348 >  gStyle->SetStatFont(42);            // same for stat. box
349    gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
350    //
351    return;
# Line 337 | Line 355 | void PVStudy::setRootStyle() {
355   //
356   std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
357   {
358 <  std::vector<PVStudy::simPrimaryVertex> simpv;
358 > std::vector<PVStudy::simPrimaryVertex> simpv;
359    const HepMC::GenEvent* evt=evtMC->GetEvent();
360    if (evt) {
361 <    if(verbose_) {
362 <      cout << "process id " << evt->signal_process_id()<<endl;
363 <      cout << "signal process vertex " << ( evt->signal_process_vertex() ?
364 <                                            evt->signal_process_vertex()->barcode() : 0 )   <<endl;
347 <      cout << "number of vertices " << evt->vertices_size() << endl;
348 <    }
361 >    edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
362 >    edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
363 >                                                                        evt->signal_process_vertex()->barcode() : 0 ) <<endl;
364 >    edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
365  
366      int idx=0; int npv=0;
367      for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
# Line 356 | Line 372 | std::vector<PVStudy::simPrimaryVertex> P
372        // t component of PV:
373        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
374              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
375 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
375 >        //        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
376          HepMC::GenVertex * mv=(*mother)->production_vertex();
377          if( ((*mother)->status() == 3) && (!mv)) {
378 <          //        std::cout << "npv= " << npv << std::endl;
378 >          //        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
379            if (npv == 0) {
380              h_gen->Fill1d("genPart_cT", pos.t()); // mm
381              h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
# Line 370 | Line 386 | std::vector<PVStudy::simPrimaryVertex> P
386        //       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
387          
388        bool hasMotherVertex=false;
389 <      if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
389 >      if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
390        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
391              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
392          HepMC::GenVertex * mv=(*mother)->production_vertex();
393 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
393 >        // if (verbose_) cout << "Status = " << (*mother)->status() << endl;
394          if (mv) {
395            hasMotherVertex=true;
396            if(!verbose_) break; //if verbose_, print all particles of gen vertices
# Line 401 | Line 417 | std::vector<PVStudy::simPrimaryVertex> P
417  
418        if(!vp){
419          // this is a new vertex
420 <        if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
420 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
421          simpv.push_back(sv);
422          vp=&simpv.back();
423        }else{
424 <        if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
424 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
425        }
426        vp->genVertex.push_back((*vitr)->barcode());
427        // collect final state descendants
# Line 516 | Line 532 | PVStudy::analyze(const edm::Event& iEven
532    using namespace edm;
533    using namespace reco;
534    
535 <  if (verbose_)
536 <    cout<<"PVStudy::analyze():"<<endl;
537 <  //=======================================================
535 >  //========================================================================
536 >  // Step 0: Prepare root variables and get information from the Event
537 >  //========================================================================
538 >  
539 >  edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
540    // Initialize Root-tuple variables if needed
523 –  //=======================================================
524 –  //if(saventuple_)
541    SetVarToZero();
542    
543 <  //=======================================================
528 <  // Track accessors
529 <  //=======================================================
530 <  //trackCollection
543 >  // ====== TrackCollection
544    static const reco::TrackCollection s_empty_trackColl;
545    const reco::TrackCollection *trackColl = &s_empty_trackColl;
546    edm::Handle<reco::TrackCollection> trackCollectionHandle;
# Line 535 | Line 548 | PVStudy::analyze(const edm::Event& iEven
548    if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
549      trackColl = trackCollectionHandle.product();
550    } else {
551 <    cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
551 >    edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
552    }
553 <
541 <  //splitTrackCollection1
553 >  // ====== splitTrackCollection1
554    static const reco::TrackCollection s_empty_splitTrackColl1;
555    const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
556    edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
# Line 546 | Line 558 | PVStudy::analyze(const edm::Event& iEven
558    if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
559      splitTrackColl1 = splitTrackCollection1Handle.product();
560    } else {
561 <    cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
561 >    edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
562    }
563 <
552 <  //splitTrackCollection2
563 >  // ====== splitTrackCollection2
564    static const reco::TrackCollection s_empty_splitTrackColl2;
565    const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
566    edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
# Line 557 | Line 568 | PVStudy::analyze(const edm::Event& iEven
568    if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
569      splitTrackColl2 = splitTrackCollection2Handle.product();
570    } else {
571 <    cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
571 >   edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
572    }
573 <
574 <
564 <  //=======================================================
565 <  // Fill trackparameters of the input tracks to pvtx fitter
566 <  //=======================================================
567 <  if(verbose_)
568 <    cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
569 <  //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
570 <  // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
571 <  fillTrackHisto(trackColl, 0);
572 <  fillTrackHisto(splitTrackColl1, 1);
573 <  fillTrackHisto(splitTrackColl2, 2);
574 <  if(verbose_)
575 <    cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
576 <
577 <  //=======================================================
578 <  // PVTX accessors
579 <  //=======================================================
580 <  //vertexCollection
573 >  
574 >  // ======= PrimaryVertexCollection
575    static const reco::VertexCollection s_empty_vertexColl;
576    const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
577    edm::Handle<reco::VertexCollection> vertexCollectionHandle;
# Line 585 | Line 579 | PVStudy::analyze(const edm::Event& iEven
579    if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
580      vertexColl = vertexCollectionHandle.product();
581    } else {
582 <    cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
582 >   edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
583    }
584 <
591 <  //splitVertexCollection1
584 >  // ====== splitVertexCollection1
585    static const reco::VertexCollection s_empty_splitVertexColl1;
586    const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
587    edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
# Line 596 | Line 589 | PVStudy::analyze(const edm::Event& iEven
589    if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
590      splitVertexColl1 = splitVertexCollection1Handle.product();
591    } else {
592 <    cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593 <  }
594 <
602 <  //splitVertexCollection2
592 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593 >  }
594 >  // ====== splitVertexCollection2
595    static const reco::VertexCollection s_empty_splitVertexColl2;
596    const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
597    edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
# Line 607 | Line 599 | PVStudy::analyze(const edm::Event& iEven
599    if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
600      splitVertexColl2 = splitVertexCollection2Handle.product();
601    } else {
602 <    cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
602 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
603    }
604 +
605 +  
606 +  // ======== BeamSpot accessors
607 +  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
608 +  iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
609 +  reco::BeamSpot bs = *recoBeamSpotHandle;    
610 +  const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
611    
612 <  if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
612 >  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
613    
614 <  //=======================================================
616 <  // MC simvtx accessor
617 <  //=======================================================
614 >  // ========== MC simvtx accessor
615    if (!realData_) {
616      edm::Handle<SimVertexContainer> simVtxs;
617      iEvent.getByLabel( simG4_, simVtxs);
# Line 623 | Line 620 | PVStudy::analyze(const edm::Event& iEven
620      iEvent.getByLabel( simG4_, simTrks);
621    }
622  
623 <  //=======================================================
627 <  // GET PDT
628 <  //=======================================================
623 >  // ========== GET PDT
624    try{
625      iSetup.getData(pdt);
626    }catch(const Exception&){
627 <    cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
627 >    edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
628    }
629    
630 <  //=======================================================
631 <  // GET pixelVertices
632 <  //=======================================================
633 <  static const reco::VertexCollection s_empty_pixelVertexColl;
634 <  const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
635 <  edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
636 <  iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
637 <  if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
638 <    pixelVertexColl = pixelVertexCollectionHandle.product();
639 <  } else {
640 <    cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
646 <  }
647 <
648 <  //=======================================================
649 <  // Fill pixelVertices related histograms
650 <  //=======================================================
651 <  nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
652 <  if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
653 <    //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
654 <    fillTrackHistoInPV(pixelVertexColl, 4, false, true);
655 <    h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());    
656 <    // Get the dZ error of the tracks in the leading pixelVertexColl
657 <    for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
658 <        t!= (pixelVertexColl->begin())->tracks_end(); t++) {
659 <      
660 <      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
661 <        //        h_pvtrk->Fill1d("trkPtPV", 0.);
662 <      }
663 <      else
664 <        h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
665 <    }
666 <    // Fill the track->dz() in the PV difference from first pixelpvtx
667 <    for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
668 <        t!= (vertexColl->begin())->tracks_end(); t++) {
669 <      // illegal charge
670 <      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
671 <        //        h_pvtrk->Fill1d("trkPtPV", 0.);
672 <      }
673 <      else {
674 <        if(pixelVertexColl->size()>0) {
675 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
676 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
677 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
678 <          h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
679 <        }
680 <      }
630 >  // ======= Get MC information if needed
631 >  bool MC=false;  
632 >  Handle<HepMCProduct> evtMC;
633 >  if (!realData_) {
634 >    iEvent.getByLabel("generator",evtMC);
635 >    if (!evtMC.isValid()) {
636 >      MC=false;
637 >      edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
638 >    } else {
639 >      edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
640 >      MC=true;
641      }
642    }
683 –
684 –  //=======================================================
685 –  // Fill number of reconstructed vertices
686 –  //=======================================================
643  
644 <  if(verbose_)  {
645 <    cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
644 >  //========================================================================
645 >  // Step 1:  Apply event cleaning for data and MC
646 >  //          WARNING: event selection cut are hard coded!!
647 >  //========================================================================
648 >  // =====Loop over the trackColl to get the fraction of HighPurity tracks
649 >  int nTracks = 0;
650 >  int nHighPurityTracks=0;
651 >  double fHighPurity=0;
652 >
653 >  for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
654 >    TrackRef tkref(trackCollectionHandle,i);
655 >    if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
656 >  }
657 >  if(nTracks>0)
658 >    fHighPurity =  double(nHighPurityTracks)/double(nTracks);
659 >  if(verbose_)
660 >    cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
661 >    
662 >  if( (fHighPurity<fHighPuritycut_ && nTracks>10) || vertexColl->begin()->isFake())  return;
663 >  
664 >  glob_runno_ = iEvent.id().run();
665 >  glob_evtno_ = iEvent.id().event();
666 >  glob_ls_   = iEvent.luminosityBlock();
667 >  glob_bx_  = iEvent.bunchCrossing();  
668 >  
669 >  //========================================================================
670 >  // Step 2: Fill histograms for the splitting consistency checks
671 >  //========================================================================
672 >  
673 >  // === Fill trackparameters of the input tracks to pvtx fitter
674 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
675 >  //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
676 >  // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
677 >  fillTrackHisto(trackColl, 0, beamSpot);
678 >  fillTrackHisto(splitTrackColl1, 1, beamSpot);
679 >  fillTrackHisto(splitTrackColl2, 2, beamSpot);
680 >  edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
681 >  
682 >  
683 >  // ==== Fill number of reconstructed vertices
684 >  edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
685 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
686 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
687 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
688 >  if (verbose_) {
689      printRecVtxs(vertexCollectionHandle);
691 –    cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
690      printRecVtxs(splitVertexCollection1Handle);
693 –    cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
691      printRecVtxs(splitVertexCollection2Handle);
695 –    cout<<"Start filling pvtx parameters."<<endl;
692    }
693 <  
693 >
694    nrecPV_ = int (vertexColl->size());
695    nrecPV1_spl_ = int (splitVertexColl1->size());
696    nrecPV2_spl_ = int (splitVertexColl2->size());
# Line 704 | Line 700 | PVStudy::analyze(const edm::Event& iEven
700    h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
701    h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
702    
707 –  //=================================================================
708 –  // Fill track parameter ntuple/hist for tracks used in recoVertices
709 –  //=================================================================
710 –
711 –  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
712 –  if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
713 –    fillTrackHistoInPV(vertexColl, 0, true, true);
703    
704 <  if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
705 <    fillTrackHistoInPV(splitVertexColl1, 1, false, true);
706 <
707 <  if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
719 <    fillTrackHistoInPV(splitVertexColl2, 2, false, true);
720 <
721 <
722 <  //=======================================================
723 <  // Compare offlinePrimaryVertices with pixelVertices
724 <  //=======================================================
725 <  if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
726 <      && (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {    
727 <    h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
728 <    // difference in reconstructed position of the leading pvtx
729 <    //cout<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;  
730 <    //cout<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
731 <    h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
732 <    h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
733 <    h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
734 <  }
735 <  
736 <
737 <
738 <  //==========================================================
739 <  // Fill secondary/primary min separations for multi-vertices
740 <  //==========================================================
704 >  // ======= Fill track parameter ntuple/hist for tracks used in recoVertices
705 >  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
706 >  if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
707 >    fillTrackHistoInPV(vertexColl, 0, beamSpot);
708    
709 +  // ======= Fill secondary/primary min separations for multi-vertices
710    if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {  
743 –
711      double min_xsep = 9999.0;
712      double min_ysep = 9999.0;
713      double min_zsep = 9999.0;    
# Line 750 | Line 717 | PVStudy::analyze(const edm::Event& iEven
717      double min_ntrksep = 9999.0;
718      double min_sumpt2sep = 9999999.0;
719      
720 <    if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
720 >    edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
721  
722      // Looping through the secondary vertices to find the mininum separation to the primary
723      for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
# Line 797 | Line 764 | PVStudy::analyze(const edm::Event& iEven
764      min_zsep_ = min_zsep;
765      min_ntrksep_ = min_ntrksep;
766      min_sumpt2sep_ = min_sumpt2sep;
800 –
801 –
767    } // End of  if(vertexColl->size()>1) {
768    
769 <  if(verbose_)
805 <    cout<<"End filling pvtx parameters."<<endl;
806 <  
807 <  //=======================================================
808 <  //           PrimaryVertex Matching
809 <  // 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
810 <  // 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
811 <  // Assume the first match is the primary vertex,
812 <  // since vertexColl are sorted in decreasing order of sum(pT)
813 <  //=======================================================
769 >  edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
770  
771 <  if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
771 >  //========================================================================
772 >  // Step 3:  PrimaryVertex Matching
773 >  //   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
774 >  //   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
775 >  //   3. The first match is the primary vertex, which has largest sum(pT^2)
776 >  //========================================================================
777 >  
778 >  edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
779    reco::VertexCollection s_empty_matchedVertexColl1;
780    reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
781    matchedVertexColl1->reserve(splitVertexColl1->size());
# Line 831 | Line 794 | PVStudy::analyze(const edm::Event& iEven
794        if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
795        reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
796        if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
797 <      if(verbose_) {
798 <        cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
799 <        cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
800 <        cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
801 <      }
797 >      edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
798 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
799 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
800 >      
801 >      double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
802 >      h_misc->Fill1d("twovtxzsign", twovtxsig);
803        if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
804 <        if(verbose_) {
805 <          cout<<"The two splitted vertices match in Z. "<<endl;
804 >        edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
805 >
806 >        int nTrkPV1, nTrkPV2;
807 >        if(useHWTrk_) {
808 >          nTrkPV1 = nHWTrkRecVtx(*recvtx1);
809 >          nTrkPV2 = nHWTrkRecVtx(*recvtx2);
810 >        }
811 >        else {
812 >          nTrkPV1 = recvtx1->tracksSize();
813 >          nTrkPV2 = recvtx2->tracksSize();
814          }
843 –        int nTrkPV1 = recvtx1->tracksSize();
844 –        int nTrkPV2 = recvtx2->tracksSize();    
815          double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
816          h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
817          h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
818          if(fabs(ntrkreldiff)<ntrkdiffcut_) {    
819 <          if(verbose_) {
820 <            cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
851 <          }
819 >          edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
820 >
821            matchedVertexColl1->push_back(*recvtx1);
822            matchedVertexColl2->push_back(*recvtx2);
823            stopmatching = true; // stop the matching when the first match is found!
824            break;
825          }
826          else
827 <          cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
827 >          edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
828        }
829        else {
830 <        cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
831 <        cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
832 <        cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;  
830 >        edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
831 >        edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
832 >        edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;  
833        }
834      }
835    }
836 <  if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
837 <  
838 <  //=======================================================
839 <  //  Analyze matchedVertexColls
840 <  //=======================================================
841 <  if(verbose_) {
842 <    cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
843 <  }    
836 >  edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
837 >
838 >  //========================================================================
839 >  // Step 4: Analyze matchedVertexColls
840 >  //========================================================================
841 >  edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
842 >  
843 >  // ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
844 >  if(MC){
845 >    // make a list of primary vertices:
846 >    std::vector<simPrimaryVertex> simpv;
847 >    simpv=getSimPVs(evtMC,"");
848 >    //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
849 >    h_gen->Fill1d("nsimPV", simpv.size());
850 >    nsimPV_ = simpv.size();
851 >    int isimpv = 0;
852 >    for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
853 >        vsim!=simpv.end(); vsim++, isimpv++){
854 >      nsimTrkPV_[isimpv]  =vsim->nGenTrk;
855 >      simx_[isimpv] = vsim->x;
856 >      simy_[isimpv] = vsim->y;
857 >      simz_[isimpv] = vsim->y;
858 >      simptsq_[isimpv] = vsim->ptsq;
859 >      //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
860 >      fillMCHisto(vsim, isimpv, vertexColl, 3);
861 >    }
862 >  }
863    
864    // Compare the reconstructed vertex position and calculate resolution/pulls
865    if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {  
866 <    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
867 <    fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
868 <    fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
866 >    // ==== Analyze the splitted vtx using MC method
867 >    if(MC){
868 >      // make a list of primary vertices:
869 >      std::vector<simPrimaryVertex> simpv;
870 >      simpv=getSimPVs(evtMC,"");
871 >      //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
872 >      int isimpv = 0;
873 >      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
874 >          vsim!=simpv.end(); vsim++, isimpv++){
875 >        //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
876 >        fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
877 >        fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
878 >      }
879 >    }
880 >
881 >    // ==== Analyze res/pull two-vertex method
882 >    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
883 >    fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
884 >    fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
885  
886      reco::VertexCollection::const_iterator v1;
887      reco::VertexCollection::const_iterator v2;
888 +
889      nrecPV_twovtx_ = 0;
890      for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
891          v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
892 <        ++v1, ++v2) {
892 >        ++v1, ++v2) {
893 >      h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
894 >      h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
895        fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
896        fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
897        fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
898        ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
899        ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
900        ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
894 –      fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
901        
902 +      int nTrkPV1, nTrkPV2;
903 +      if(useHWTrk_) {
904 +        nTrkPV1 = nHWTrkRecVtx(*v1);
905 +        nTrkPV2 = nHWTrkRecVtx(*v2);
906 +      }
907 +      else {
908 +        nTrkPV1 = v1->tracksSize();
909 +        nTrkPV2 = v2->tracksSize();
910 +      }
911 +      
912 +      fntrk_ = (nTrkPV1 + nTrkPV2)/2;
913 +
914        h_summary->Fill1d("deltax", fres_[0] );
915        h_summary->Fill1d("deltay", fres_[1] );
916        h_summary->Fill1d("deltaz", fres_[2] );
# Line 909 | Line 927 | PVStudy::analyze(const edm::Event& iEven
927        fillHisto(res(fres_[0], fres_[1], fres_[2]),
928                  error(ferror_[0], ferror_[1], ferror_[2]),
929                  fntrk_,0);  
930 +
931        // Fill the ntuple variables
913 –      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
914 –      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
915 –      sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
916 –      sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
932        nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
933        deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
934        deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
# Line 924 | Line 939 | PVStudy::analyze(const edm::Event& iEven
939        pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
940        pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
941        pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];    
942 +
943 +
944 +      //SplittedVertex
945 +      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = nTrkPV1;
946 +      ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
947 +      normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
948 +      avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
949 +      errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
950 +      erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
951 +      errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
952 +      
953 +      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = nTrkPV2;
954 +      ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
955 +      normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
956 +      avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
957 +      errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
958 +      erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
959 +      errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
960 +  
961        nrecPV_twovtx_++;
962 +      
963 +      // Print some information of the two tracks events
964 +      if(verbose_ && fntrk_ < 4) {
965 +        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
966 +        printRecVtx(*v1);
967 +        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
968 +        printRecVtx(*v2);
969 +      }
970      } // End of analyzing res/pull
971    } // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
972    else
973 <    cout<<"WARNING: Cannot find matching pvtxs"<<endl;
973 >    edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
974    
975 <  if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
934 <
935 <  //=======================================================
936 <  // Fill simulated vertices
937 <  //=======================================================
938 <  if (!realData_) {
939 <    bool MC=false;
940 <    Handle<HepMCProduct> evtMC;
941 <    iEvent.getByLabel("generator",evtMC);
942 <    if (!evtMC.isValid()) {
943 <      MC=false;
944 <      if(verbose_){
945 <        cout << "no HepMCProduct found"<< endl;
946 <      }
947 <    } else {
948 <      if(verbose_){
949 <        cout << "generator HepMCProduct found"<< endl;
950 <      }
951 <      MC=true;
952 <    }
953 <    
954 <    if(MC){
955 <      // make a list of primary vertices:
956 <      std::vector<simPrimaryVertex> simpv;
957 <      simpv=getSimPVs(evtMC,"");
958 <      //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
959 <      h_gen->Fill1d("nsimPV", simpv.size());
960 <      nsimPV_ = simpv.size();
961 <
962 <      int isimpv = 0;
963 <      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
964 <          vsim!=simpv.end(); vsim++, isimpv++){
965 <        nsimTrkPV_[isimpv]  =vsim->nGenTrk;
966 <        simx_[isimpv] = vsim->x;
967 <        simy_[isimpv] = vsim->y;
968 <        simz_[isimpv] = vsim->y;
969 <        simptsq_[isimpv] = vsim->ptsq;
970 <        
971 <        //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
972 <        fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
973 <        fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
974 <        fillMCHisto(vsim, isimpv, vertexColl, 3);
975 <      }
976 <    } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
977 <  } // End of   if (!realData_) {
975 >  edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
976  
977 <  // Finally fill the ftree_
977 >
978 >  //========================================================================
979 >  // Step 5: Fill ntuple if saventuple_ is on
980 >  //========================================================================
981    if(saventuple_) {
982      ftree_->Fill();
983      pvtxtree_->Fill();
984    }
985 +  
986   }
987  
988 +
989 +  
990   // ------------ method called once each job just before starting event loop  ------------
991   void
992   PVStudy::beginJob()
# Line 992 | Line 996 | PVStudy::beginJob()
996   // ------------ method called once each job just after ending the event loop  ------------
997   void
998   PVStudy::endJob() {
999 <  if (verbose_) cout << "Analyzing PV info" << endl;
999 >  edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
1000    // Looping through the datamodes available
1001    for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {  
1002      string suffix;
1003 <    if(verbose_)  cout<<"datamode = "<< *it<<endl;
1003 >    edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1004      switch(*it) {
1005      case 0: suffix = "";
1006        break;
# Line 1007 | Line 1011 | PVStudy::endJob() {
1011      case 3: suffix = "_mct";
1012        break;
1013      }
1014 +    suffix += "_nTrk";
1015      if(analyze_) {
1016        for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1017 <        cout << "ntrk = " << ntrk << endl;
1017 >        edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1018          PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1019 <        if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix+"_nTrk", ntrk,ipv.res_.x(), 1.);
1020 <        if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix+"_nTrk", ntrk,ipv.res_.y(), 1.);
1021 <        if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix+"_nTrk", ntrk,ipv.res_.z(), 1.);
1022 <        if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix+"_nTrk", ntrk,ipv.pull_.x(), 1.);
1023 <        if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix+"_nTrk", ntrk,ipv.pull_.y(), 1.);
1024 <        if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix+"_nTrk", ntrk,ipv.pull_.z(), 1.);
1019 >        if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1020 >        if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1021 >        if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1022 >        if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1023 >        if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1024 >        if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1025        }
1026 <    }
1026 >    }
1027      suffix.clear();
1028    }
1029    if(saventuple_) {
# Line 1027 | Line 1032 | PVStudy::endJob() {
1032      pvtxtree_->Write();
1033      file_->Close();
1034    }
1035 < }
1031 <
1035 >  
1036  
1037 + }
1038  
1039   // Match a recovertex with a simvertex
1040   // ? Should the cut be parameter dependent?
# Line 1037 | Line 1042 | PVStudy::endJob() {
1042   bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1043                            const reco::Vertex & vrec)
1044   {
1045 <  return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1045 >  if(vrec.isFake() || !vrec.isValid()) return false;
1046 >  else
1047 >    return true;
1048 >    //return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1049   }
1050  
1051   // Match two reconstructed vertices
# Line 1046 | Line 1054 | bool PVStudy::matchVertex( const reco::V
1054                             int zsigncut)
1055   {
1056    double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1057 <  if(verbose_)
1050 <    cout<<"The matching criteria in z is "<<cut<<endl;
1057 >  edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1058    return (fabs(vrec1->z()-vrec2->z())<cut);
1059   }
1060  
1061  
1062   bool PVStudy::isResonance(const HepMC::GenParticle * p){
1063    double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1064 <  if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1064 >  edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1065    return  ctau >0 && ctau <1e-6;
1066   }
1067  
# Line 1105 | Line 1112 | void PVStudy::printRecVtxs(const Handle<
1112    }
1113   }
1114  
1115 + void PVStudy::printRecVtx(const reco::Vertex & v){
1116 +
1117 +  cout << "#trk " << std::setw(3) << v.tracksSize()
1118 +         << " chi2 " << std::setw(4) << v.chi2()
1119 +         << " ndof " << std::setw(3) << v.ndof() << endl
1120 +         << " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1121 +         << " dx " << std::setw(8) << v.xError()<< endl
1122 +         << " y "  << std::setw(8) << v.y()
1123 +         << " dy " << std::setw(8) << v.yError()<< endl
1124 +         << " z "  << std::setw(8) << v.z()
1125 +         << " dz " << std::setw(8) << v.zError()
1126 +         << endl;
1127 + }
1128 +
1129   void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1130    cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1131    int i=1;
# Line 1126 | Line 1147 | void PVStudy::printSimTrks(const Handle<
1147   void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1148    stringstream ss;
1149    ss << ntk;
1150 <  string suffix = ss.str();
1151 <  if(datamode == 0 )  suffix = "_" + suffix;
1152 <  if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1153 <  if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1154 <  if(datamode == 3 )  suffix = "_mct_" + suffix;
1150 >  string suffix;
1151 >  if(datamode == 0 )  suffix = "_" + ss.str();
1152 >  if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1153 >  if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1154 >  if(datamode == 3 )  suffix = "_mct_" + ss.str();
1155    
1156    if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1157    // Fill histograms of res and pull of ntk
# Line 1145 | Line 1166 | void PVStudy::fillHisto(res r, error er,
1166    h_others->Fill1d("errPVz"+suffix, er.z());
1167   }
1168  
1169 +
1170   PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1171    map<int,double> data;
1172    data.clear();
# Line 1192 | Line 1214 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1214          data[i] = fpar[0];
1215          break;
1216        }
1217 <      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1217 >      if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1218      }
1219    }
1220    else
# Line 1207 | Line 1229 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1229                              ntk);
1230   }
1231  
1232 +
1233   void PVStudy::fit(TH1 *hdist, double fitPars[]){
1234    TAxis *axis0 = hdist->GetXaxis();
1235    int nbin = axis0->GetLast();
1236    double nOF = hdist->GetBinContent(nbin+1);
1237    double nUF = hdist->GetBinContent(0);
1238    //   double fitRange = axis0->GetBinUpEdge(nbin);
1239 <  double fitRange = axis0->GetXmax();
1239 >  // double fitRange = axis0->GetXmax();
1240 >  double fitRange = 2*hdist->GetRMS();
1241    double sigMax[2] = {0.,0.};
1242  
1243 <  if ( verbose_ ){
1244 <    cout << "Last bin = " << nbin;
1245 <    cout << "; hdist: Overflow Entries = " << nOF;
1246 <    cout << "; Underflow Entries = " << nUF;
1247 <    cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1248 <    cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1225 <  }
1243 >  edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1244 >                           << "; hdist: Overflow Entries = " << nOF
1245 >                           << "; Underflow Entries = " << nUF
1246 >                           << "; hdist->GetEntries() = " << hdist->GetEntries()
1247 >                           << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1248 >
1249    if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1250      for (int bi = 0; bi < nbin; bi++) {
1251        if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
# Line 1234 | Line 1257 | void PVStudy::fit(TH1 *hdist, double fit
1257          if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1258        }
1259      }
1260 <    if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1260 >    //edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1261  
1262      TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1263      fgaus->SetParameter(1, 0.);
1264 <    fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1264 >    fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1265      fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1266      fgaus->SetLineColor(4);
1267 <    hdist->Fit(fgaus,"Q");
1267 >    hdist->Fit(fgaus,"QLRM");
1268      gPad->Update();
1269      TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1270      s->SetOptStat(1111111);
1271      s->SetOptFit(0111);
1272      gPad->Update();
1273 <    if (verbose_) cout << "got function" << endl;
1273 >    //edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1274      fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1275    }
1276    else
1277      fitPars[0] = -999;
1278   }
1279  
1280 < void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1280 >
1281 > void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1282    string suffix;
1283    suffix = ""; // for unsplit tracks
1284    if(datatype == 1) suffix = "1_spl"; // for splittracks 1
# Line 1266 | Line 1290 | void PVStudy::fillTrackHisto(const reco:
1290        ++itTrack) {
1291      h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1292      h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1293 +    h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1294      h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1295      h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1296      h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1297    }
1298   }
1299  
1300 < void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1300 > void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1301    string suffix;
1302    suffix = ""; // for unsplit tracks
1303    if(datatype == 1) suffix = "1_spl"; // for splittracks 1
# Line 1280 | Line 1305 | void PVStudy::fillTrackHistoInPV(const r
1305    int ivtx = 0;
1306    for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1307        v!=vertexColl->end(); ++v, ++ivtx) {  
1308 <    if(fillNtuple) {
1309 <      if(datatype == 0) {
1310 <        nTrkPV_[ivtx] = v->tracksSize();        
1311 <        sumptsq_[ivtx] = sumPtSquared(*v);
1312 <        isValid_[ivtx] = v->isValid();
1313 <        isFake_[ivtx] = v->isFake();    
1314 <        recx_[ivtx] = v->x();
1315 <        recy_[ivtx] = v->y();
1316 <        recz_[ivtx] = v->z();
1317 <        recx_err_[ivtx] = v->xError();
1318 <        recy_err_[ivtx] = v->yError();
1319 <        recz_err_[ivtx] = v->zError();
1320 <      }
1321 <      if(datatype == 1) {
1322 <        nTrkPV1_spl_[ivtx] = v->tracksSize();  
1323 <        sumptsq1_spl_[ivtx] = sumPtSquared(*v);        
1324 <        isValid1_spl_[ivtx] = v->isValid();
1325 <        isFake1_spl_[ivtx] = v->isFake();      
1326 <        recx1_spl_[ivtx] = v->x();
1302 <        recy1_spl_[ivtx] = v->y();
1303 <        recz1_spl_[ivtx] = v->z();
1304 <        recx1_err_spl_[ivtx] = v->xError();
1305 <        recy1_err_spl_[ivtx] = v->yError();
1306 <        recz1_err_spl_[ivtx] = v->zError();
1307 <
1308 <      }
1309 <      if(datatype == 2) {
1310 <        nTrkPV2_spl_[ivtx] = v->tracksSize();  
1311 <        sumptsq2_spl_[ivtx] = sumPtSquared(*v);        
1312 <        isValid2_spl_[ivtx] = v->isValid();
1313 <        isFake2_spl_[ivtx] = v->isFake();
1314 <        recx2_spl_[ivtx] = v->x();
1315 <        recy2_spl_[ivtx] = v->y();
1316 <        recz2_spl_[ivtx] = v->z();
1317 <        recx2_err_spl_[ivtx] = v->xError();
1318 <        recy2_err_spl_[ivtx] = v->yError();
1319 <        recz2_err_spl_[ivtx] = v->zError();
1320 <      }
1321 <      if(datatype == 4) { // for pixelVertices
1322 <        nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();        
1323 <        sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);      
1324 <        isValid_pxlpvtx_[ivtx] = v->isValid();
1325 <        isFake_pxlpvtx_[ivtx] = v->isFake();
1326 <        recx_pxlpvtx_[ivtx] = v->x();
1327 <        recy_pxlpvtx_[ivtx] = v->y();
1328 <        recz_pxlpvtx_[ivtx] = v->z();
1329 <        recx_err_pxlpvtx_[ivtx] = v->xError();
1330 <        recy_err_pxlpvtx_[ivtx] = v->yError();
1331 <        recz_err_pxlpvtx_[ivtx] = v->zError();
1332 <      }
1333 <    }
1334 <  }
1335 <  // For histogram only fill the leading pvtx parameters
1336 <  if(fillHisto) {
1337 <    h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1338 <    try {
1339 <      for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1340 <          t!=vertexColl->begin()->tracks_end(); t++) {
1341 <        // illegal charge
1342 <        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1343 <          //      h_pvtrk->Fill1d("trkPtPV", 0.);
1344 <        }
1345 <        else {
1346 <          h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1347 <          h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1348 <          h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1349 <          h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1350 <          h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1308 >    if(!v->isFake()) {
1309 >      h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1310 >      h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());  
1311 >      h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkRecVtx(*v));
1312 >      try {
1313 >        for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1314 >            t!=vertexColl->begin()->tracks_end(); t++) {
1315 >          // illegal charge
1316 >          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1317 >            //    h_pvtrk->Fill1d("trkPtPV", 0.);
1318 >          }
1319 >          else {
1320 >            h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1321 >            h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());  
1322 >            h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1323 >            h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1324 >            h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1325 >            h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1326 >          }
1327          }
1328        }
1329 <    }
1330 <    catch (...) {
1331 <      // exception thrown when trying to use linked track
1332 <      //          h_pvtrk->Fill1d("trkPtPV", 0.);
1329 >      catch (...) {
1330 >        // exception thrown when trying to use linked track
1331 >        //        h_pvtrk->Fill1d("trkPtPV", 0.);
1332 >      }
1333      }
1334    }
1359 –  
1335   }
1336  
1337   void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
# Line 1375 | Line 1350 | void PVStudy::fillMCHisto(std::vector<si
1350      suffix = "_mct";
1351      nrecPV_mct_ = 0;  
1352    }
1353 <  
1353 >
1354    //========================================================
1355 <  //  look for a matching reconstructed vertex in vertexColl
1356 <  //========================================================        
1357 <
1355 >  //  For each simulated vertex, look for a match in the vertexColl
1356 >  //  If more than 1 recVtx is found, use the one with closest in Z
1357 >  //========================================================  
1358 >    
1359 >  // === look for a matching reconstructed vertex in vertexColl
1360    for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1361        vrec!=vtxColl->end(); ++vrec){
1385 –    int nrectrk = vrec->tracksSize();
1362      vsim->recVtx=NULL;  
1363 <  
1364 <    if(verbose_){
1365 <      cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1366 <      cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1391 <    }
1392 <    if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1363 >    edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1364 >    edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1365 >    
1366 >    if ( matchVertex(*vsim,*vrec)  ) {
1367        vsim->recVtx=&(*vrec);
1368 <      
1395 <      // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1368 >      //if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1369        //if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1370        //|| (!vsim->recVtx) )
1371        //vsim->recVtx=&(*vrec);
1399 –      //}
1400 –      if(verbose_)
1401 –        cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1402 –      
1403 –      double fres_mct[3];
1404 –      double ferror_mct[3];
1405 –      
1406 –      fres_mct[0] = vsim->recVtx->x()-vsim->x;
1407 –      fres_mct[1] = vsim->recVtx->y()-vsim->y;
1408 –      fres_mct[2] = vsim->recVtx->z()-vsim->z;
1409 –      ferror_mct[0] = vsim->recVtx->xError();
1410 –      ferror_mct[1] = vsim->recVtx->yError();
1411 –      ferror_mct[2] = vsim->recVtx->zError();
1412 –      
1413 –      h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1414 –      h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1415 –      h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1416 –      h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1417 –      h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1418 –      h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1419 –      h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1420 –      h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1421 –      h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1422 –      pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1423 –                                       error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1424 –                                       nrectrk));
1425 –      // Fill histo according to its track multiplicity
1426 –      fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1427 –                error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1428 –                nrectrk, datatype);
1429 –      
1430 –      if(saventuple_) {
1431 –        //Fill the values for variables in ftree_  
1432 –        if(datatype == 1) {
1433 –          nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1434 –          deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1435 –          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1436 –          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1437 –          pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];  
1438 –          pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1439 –          pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1440 –          nrecPV_spl1_mct_++;
1441 –        }
1442 –        
1443 –        if(datatype == 2) {
1444 –          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1445 –          deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1446 –          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1447 –          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1448 –          pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];  
1449 –          pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1450 –          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1451 –          nrecPV_spl2_mct_++;
1452 –        }
1453 –        if(datatype == 3) {    
1454 –          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1455 –          deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1456 –          deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1457 –          deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1458 –          pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];  
1459 –          pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1460 –          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];          
1461 –          nrecPV_mct_++;
1462 –        }
1463 –      } // End of  if(saventuple_) {
1464 –    } //  if ( matchVertex(*vsim,*vrec) ) {
1465 –    else {  // no rec vertex found for this simvertex
1466 –      if(verbose_)
1467 –        cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1372      }
1373    }
1374 < }
1374 >  
1375 >  // === If match found fill the residual and pulls
1376 >  if(vsim->recVtx) {
1377 >    edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1378 >    // if study the resolution/pull for all tracks in the PVtx
1379 >    int nrectrk;
1380 >    
1381 >    if(useHWTrk_)
1382 >      nrectrk=nHWTrkRecVtx(*(vsim->recVtx));
1383 >    else
1384 >      nrectrk =  int(vsim->recVtx->tracksSize());
1385  
1386 +    double fres_mct[3];
1387 +    double ferror_mct[3];
1388 +    
1389 +    double ndofPV = vsim->recVtx->ndof();
1390 +    double normchi2PV =  vsim->recVtx->normalizedChi2();
1391 +    double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1392 +    
1393 +    fres_mct[0] = vsim->recVtx->x()-vsim->x;
1394 +    fres_mct[1] = vsim->recVtx->y()-vsim->y;
1395 +    fres_mct[2] = vsim->recVtx->z()-vsim->z;
1396 +    ferror_mct[0] = vsim->recVtx->xError();
1397 +    ferror_mct[1] = vsim->recVtx->yError();
1398 +    ferror_mct[2] = vsim->recVtx->zError();
1399 +    
1400 +    h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1401 +    h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1402 +    h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1403 +    h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1404 +    h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1405 +    h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1406 +    h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1407 +    h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1408 +    h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1409 +    pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1410 +                                     error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1411 +                                     nrectrk));
1412 +    // Fill histo according to its track multiplicity
1413 +    fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1414 +              error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1415 +              nrectrk, datatype);
1416 +    
1417 +    if(saventuple_) {
1418 +      //Fill the values for variables in pvtxtree_  
1419 +      if(datatype == 1) {
1420 +        nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1421 +        ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1422 +        normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1423 +        avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1424 +        deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1425 +        deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1426 +        deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1427 +        pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];  
1428 +        pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1429 +        pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1430 +        errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1431 +        erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1432 +        errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1433 +        nrecPV_spl1_mct_++;
1434 +      }
1435 +      if(datatype == 2) {
1436 +        nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;      
1437 +        ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1438 +        normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1439 +        avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1440 +        deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1441 +        deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1442 +        deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1443 +        pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];  
1444 +        pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1445 +        pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];          
1446 +        errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1447 +        erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1448 +        errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1449 +        nrecPV_spl2_mct_++;
1450 +      }
1451 +      if(datatype == 3) {      
1452 +        nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1453 +        ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1454 +        normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1455 +        avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1456 +        deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1457 +        deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1458 +        deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1459 +        pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];  
1460 +        pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1461 +        pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1462 +        errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1463 +        erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1464 +        errz_mct_[nrecPV_mct_] =  ferror_mct[2];          
1465 +        nrecPV_mct_++;
1466 +      }
1467 +    } // End of  if(saventuple_) {
1468 +  } //  if ( matchVertex(*vsim,*vrec) ) {
1469 +  else {  // no rec vertex found for this simvertex
1470 +    edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1471 +  }
1472 + }
1473  
1474 +                                                                                                                  
1475   void PVStudy::SetVarToZero() {
1474 –  //Greg's variables
1476    fntrk_ = 0;
1477    //pvtx position (x,y,z) residual and error
1478    for(int i = 0; i<3;i++)    {
# Line 1479 | Line 1480 | void PVStudy::SetVarToZero() {
1480      ferror_[i] = 0;
1481    }
1482    
1483 <  //Yanyan's variables
1484 <  // Number of reconstructed vertices  
1483 >  // Event level information
1484 >  glob_runno_ = 0;
1485 >  glob_evtno_ = 0;
1486 >  glob_ls_ = 0;
1487 >  glob_bx_ = 0;
1488    nrecPV_ = 0;
1489    nrecPV1_spl_ = 0;
1490    nrecPV2_spl_ = 0;
# Line 1496 | Line 1500 | void PVStudy::SetVarToZero() {
1500    min_ntrksep_ = 9999.0;
1501    min_sumpt2sep_ = -9999.0;
1502    
1503 <
1503 >  // Variables filled per Vertex
1504    for (int i = 0; i < nMaxPVs_; i++) {
1501 –    // recoVertices with all tracks
1502 –    nTrkPV_[i] = 0; // Number of tracks in the pvtx    
1503 –    sumptsq_[i] = 0;
1504 –    isValid_[i] = -1;
1505 –    isFake_[i] = -1;
1506 –    recx_[i] = 0;
1507 –    recy_[i] = 0;
1508 –    recz_[i] = 0;
1509 –    recx_err_[i] = 0;
1510 –    recy_err_[i] = 0;
1511 –    recz_err_[i] = 0;
1512 –    
1513 –    // recoVertices with splitTrack1
1514 –    nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx    
1515 –    sumptsq1_spl_[i] = 0;
1516 –    isValid1_spl_[i] = -1;
1517 –    isFake1_spl_[i] = -1;
1518 –    recx1_spl_[i] = 0;
1519 –    recy1_spl_[i] = 0;
1520 –    recz1_spl_[i] = 0;
1521 –    recx1_err_spl_[i] = 0;
1522 –    recy1_err_spl_[i] = 0;
1523 –    recz1_err_spl_[i] = 0;
1524 –  
1525 –    // recoVertices with splitTrack2
1526 –    nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx  
1527 –    sumptsq2_spl_[i] = 0;
1528 –    isValid2_spl_[i] = -1;
1529 –    isFake2_spl_[i] = -1;
1530 –    recx2_spl_[i] = 0;
1531 –    recy2_spl_[i] = 0;
1532 –    recz2_spl_[i] = 0;
1533 –    recx2_err_spl_[i] = 0;
1534 –    recy2_err_spl_[i] = 0;
1535 –    recz2_err_spl_[i] = 0;
1536 –    
1537 –    //pixelVertices
1538 –    nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx  
1539 –    sumptsq_pxlpvtx_[i] = 0;
1540 –    isValid_pxlpvtx_[i] = -1;
1541 –    isFake_pxlpvtx_[i] = -1;
1542 –    recx_pxlpvtx_[i] = 0;
1543 –    recy_pxlpvtx_[i] = 0;
1544 –    recz_pxlpvtx_[i] = 0;
1545 –    recx_err_pxlpvtx_[i] = 0;
1546 –    recy_err_pxlpvtx_[i] = 0;
1547 –    recz_err_pxlpvtx_[i] = 0;
1548 –
1505      // matched two-vertices
1506      nTrkPV1_spl_twovtx_[i] = 0;
1507 +    ndofPV1_spl_twovtx_[i] = 0;
1508 +    normchi2PV1_spl_twovtx_[i] = 0;
1509 +    avgPtPV1_spl_twovtx_[i] = 0;
1510 +    errx1_spl_twovtx_[i] = 0;
1511 +    erry1_spl_twovtx_[i] = 0;
1512 +    errz1_spl_twovtx_[i] = 0;
1513 +
1514      nTrkPV2_spl_twovtx_[i] = 0;
1515 <    nTrkPV_twovtx_[i] = 0;  
1516 <    sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1517 <    sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1515 >    ndofPV2_spl_twovtx_[i] = 0;
1516 >    normchi2PV2_spl_twovtx_[i] = 0;
1517 >    avgPtPV2_spl_twovtx_[i] = 0;
1518 >    errx2_spl_twovtx_[i] = 0;
1519 >    erry2_spl_twovtx_[i] = 0;
1520 >    errz2_spl_twovtx_[i] = 0;
1521 >    
1522      deltax_twovtx_[i] = 0;
1523      deltay_twovtx_[i] = 0;
1524      deltaz_twovtx_[i] = 0;
# Line 1576 | Line 1543 | void PVStudy::SetVarToZero() {
1543      pullx_mct_[i] = 0;
1544      pully_mct_[i] = 0;
1545      pullz_mct_[i] = 0;
1546 <        
1546 >    nTrkPV_mct_[i] = 0;
1547 >    ndofPV_mct_[i] = 0;
1548 >    normchi2PV_mct_[i] = 0;
1549 >    avgPtPV_mct_[i] = 0;
1550 >    errx_mct_[i] = 0;
1551 >    erry_mct_[i] = 0;
1552 >    errz_mct_[i] = 0;
1553 >    
1554      deltax_spl1_mct_[i] = 0;
1555      deltay_spl1_mct_[i] = 0;
1556      deltaz_spl1_mct_[i] = 0;
1557      pullx_spl1_mct_[i] = 0;
1558      pully_spl1_mct_[i] = 0;
1559      pullz_spl1_mct_[i] = 0;
1560 <    
1560 >    nTrkPV_spl1_mct_[i] = 0;
1561 >    ndofPV_spl1_mct_[i] = 0;
1562 >    normchi2PV_spl1_mct_[i] = 0;
1563 >    avgPtPV_spl1_mct_[i] = 0;
1564 >    errx_spl1_mct_[i] = 0;
1565 >    erry_spl1_mct_[i] = 0;
1566 >    errz_spl1_mct_[i] = 0;
1567 >
1568      deltax_spl2_mct_[i] = 0;
1569      deltay_spl2_mct_[i] = 0;
1570      deltaz_spl2_mct_[i] = 0;
1571      pullx_spl2_mct_[i] = 0;
1572      pully_spl2_mct_[i] = 0;
1573      pullz_spl2_mct_[i] = 0;
1574 +    nTrkPV_spl2_mct_[i] = 0;
1575 +    ndofPV_spl2_mct_[i] = 0;
1576 +    normchi2PV_spl2_mct_[i] = 0;
1577 +    avgPtPV_spl2_mct_[i] = 0;
1578 +    errx_spl2_mct_[i] = 0;
1579 +    erry_spl2_mct_[i] = 0;
1580 +    errz_spl2_mct_[i] = 0;
1581 +
1582    }
1594 –  
1583   }
1584  
1585   double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1586 <  double sum = 0.;
1587 <  double pT;
1588 <  for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1589 <    pT = (**it).pt();
1590 <    sum += pT*pT;
1586 >  if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1587 >  else {
1588 >    double sum = 0.;
1589 >    double pT;
1590 >    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1591 >      pT = (**it).pt();
1592 >      sum += pT*pT;
1593 >    }
1594 >    return sum;
1595    }
1604 –  return sum;
1596   }
1597  
1598 +
1599 +
1600 + double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1601 +
1602 +  if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1603 +  else {
1604 +    int nHWTrk = 0;
1605 +    double sumpT = 0.;
1606 +    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1607 +      if(v.trackWeight(*it) > 0.5 ) {
1608 +        sumpT +=  (**it).pt();
1609 +        nHWTrk++;
1610 +      }
1611 +  }
1612 +    if(nHWTrk > 0)
1613 +      return sumpT/double(nHWTrk);
1614 +    else
1615 +      return 0;
1616 +  }
1617 + }
1618 +
1619 + int PVStudy::nHWTrkRecVtx(const reco::Vertex & v)  {
1620 +  int nHWTrkPV = 0;
1621 +  for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1622 +    if(v.trackWeight(*it) > 0.5)
1623 +      nHWTrkPV++;
1624 +  }
1625 +  return nHWTrkPV;
1626 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines