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.17 by yygao, Tue Jan 26 12:06:49 2010 UTC vs.
Revision 1.19 by yygao, Sat Jan 30 18:27:09 2010 UTC

# Line 51 | Line 51 | Implementation:
51   #include <SimDataFormats/Track/interface/SimTrackContainer.h>
52   // BeamSpot
53   #include "DataFormats/BeamSpot/interface/BeamSpot.h"
54 // Trigger
55 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
56 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
57 #include "CondFormats/L1TObjects/interface/L1GtPrescaleFactors.h"
58 #include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsAlgoTrigRcd.h"
59 #include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsTechTrigRcd.h"
60 #include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
61 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
62 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
63 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoAlgoTrigRcd.h"
64 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoTechTrigRcd.h"
65 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
66 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
67 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
68
54  
55   //root
56   #include <TROOT.h>
# Line 95 | Line 80 | PVStudy::PVStudy(const edm::ParameterSet
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");
86    nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87    nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88    zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89    bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
104  applyEvtClean_            = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
105  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
90    
91    //now do what ever initialization is needed
92    pvinfo.clear();
93  
110
94    // Specify the data mode vector
95    if(realData_) datamode.push_back(0);
96    else {
# Line 127 | Line 110 | PVStudy::PVStudy(const edm::ParameterSet
110  
111      // pvtxtree_ analyzing the pvtxs ootb
112      pvtxtree_ = new TTree("pvtxtree","pvtxtree");
130    //pvtx using all tracks
131    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
132    pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
133    pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");    
134    pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
135    pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
136    pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
137    pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
138    pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
139    pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
140    pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
141    pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
113  
114 +    // Event information for the data
115 +    pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
116 +    pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
117 +    pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
118 +    pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");  
119 +    
120 +
121 +    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
122 +    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
123 +    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
124 +    // Event level information
125      pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126      pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127      pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128  
129 <    //pvtx using splittracks1
130 <    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
131 <    pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");  
132 <    pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");  
151 <    pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
152 <    pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
153 <    pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
154 <    pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
155 <    pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
156 <    pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
157 <    pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
158 <    pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");  
159 <  
160 <    //pvtx using splittracks2
161 <    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
162 <    pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");    
163 <    pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");  
164 <    pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
165 <    pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
166 <    pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
167 <    pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
168 <    pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
169 <    pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
170 <    pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
171 <    pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");  
172 <    
173 <    //pixeVertices
174 <    pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
175 <    pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");    
176 <    pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");  
177 <    pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
178 <    pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
179 <    pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
180 <    pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
181 <    pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
182 <    pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
183 <    pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
184 <    pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");  
129 >    //Fill the variables in the twovtx pair (recvtx1, recvtx2)  
130 >    // Information related to the analyzing the two-vertex method
131 >
132 >
133  
186    //Fill the variables in the twovtx pair (recvtx1, recvtx2)
134      pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
188    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
189    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
190    pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
191    pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
135      pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
136      pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
137      pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
# Line 199 | Line 142 | PVStudy::PVStudy(const edm::ParameterSet
142      pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
143      pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
144      pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
145 +
146 +    // Information for the splitVertexColl1
147 +    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
148 +    pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
149 +    pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
150 +    pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
151 +    pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
152 +    pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
153 +    pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
154 +  
155 +    
156 +    // Information for the splitVertexColl2
157 +    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
158 +    pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
159 +    pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
160 +    pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
161 +    pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
162 +    pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
163 +    pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
164 +  
165      
166      // MC variables
167      if(!realData_) {  
# Line 216 | Line 179 | PVStudy::PVStudy(const edm::ParameterSet
179        pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
180        pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
181        pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
182 <      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
182 >      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");    
183 >      pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
184 >      pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
185 >      pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
186 >      pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
187 >      pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
188 >      pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
189 >      pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
190 >
191        // For pvtxs with splittracks1
192        pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
193        pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
# Line 225 | Line 196 | PVStudy::PVStudy(const edm::ParameterSet
196        pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
197        pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
198        pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
199 +      pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
200 +      pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
201 +      pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
202 +      pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
203 +      pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
204 +      pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
205 +      pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
206 +        
207        // For pvtxs with splittracks1
208        pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
209        pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
# Line 232 | Line 211 | PVStudy::PVStudy(const edm::ParameterSet
211        pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
212        pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
213        pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
214 <      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
214 >      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");  
215 >      pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
216 >      pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
217 >      pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");  
218 >      pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
219 >      pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
220 >      pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
221 >      pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
222 >
223      }
224    }
225    
# Line 256 | Line 243 | PVStudy::PVStudy(const edm::ParameterSet
243  
244    // Book Histograms:
245    h_pvtrk   = new PVHistograms();
259  h_pixvtx  = new PVHistograms();
246    h_misc    = new PVHistograms();
247    h_summary = new PVHistograms();
248    h_others  = new PVHistograms();
# Line 269 | Line 255 | PVStudy::PVStudy(const edm::ParameterSet
255        h_pvtrk->Init("pvTrk", ss.str(),"spl");
256      }
257    }
272  h_pixvtx->Init("pixVtx");
258    h_misc->Init("misc");
259  
260    // Book MC only plots
# Line 318 | Line 303 | PVStudy::~PVStudy()
303    theFile->cd();
304    theFile->cd("Summary");
305    h_pvtrk->Save();
321  h_pixvtx->Save();
306    h_misc->Save();
307    h_summary->Save();
308    if (!realData_)
# Line 356 | Line 340 | void PVStudy::setRootStyle() {
340    gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
341    gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
342    gStyle->SetStatFontSize(0.08);      // size for stat. box
343 <  gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
344 <  gStyle->SetTitleFont(32, "");       // same for pad title
345 <  gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
346 <  gStyle->SetStatFont(32);            // same for stat. box
343 >  gStyle->SetTitleFont(42, "XYZ");    // times-bold-italic font (p. 153) for axes
344 >  gStyle->SetTitleFont(42, "");       // same for pad title
345 >  gStyle->SetLabelFont(42, "XYZ");    // same for axis labels
346 >  gStyle->SetStatFont(42);            // same for stat. box
347    gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
348    //
349    return;
# Line 546 | Line 530 | PVStudy::analyze(const edm::Event& iEven
530    using namespace edm;
531    using namespace reco;
532    
533 +  //========================================================================
534 +  // Step 0: Prepare root variables and get information from the Event
535 +  //========================================================================
536 +  
537    edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
550  //=======================================================
538    // Initialize Root-tuple variables if needed
539 <  //=======================================================
540 <  //if(saventuple_)
541 <    SetVarToZero();
555 <  
556 <  //=======================================================
557 <  // Track accessors
558 <  //=======================================================
559 <  //trackCollection
539 >  SetVarToZero();
540 >  
541 >  // ====== TrackCollection
542    static const reco::TrackCollection s_empty_trackColl;
543    const reco::TrackCollection *trackColl = &s_empty_trackColl;
544    edm::Handle<reco::TrackCollection> trackCollectionHandle;
# Line 566 | Line 548 | PVStudy::analyze(const edm::Event& iEven
548    } else {
549      edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
550    }
551 <  //splitTrackCollection1
551 >  // ====== splitTrackCollection1
552    static const reco::TrackCollection s_empty_splitTrackColl1;
553    const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
554    edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
# Line 576 | Line 558 | PVStudy::analyze(const edm::Event& iEven
558    } else {
559      edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
560    }
561 <  //splitTrackCollection2
561 >  // ====== splitTrackCollection2
562    static const reco::TrackCollection s_empty_splitTrackColl2;
563    const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
564    edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
# Line 587 | Line 569 | PVStudy::analyze(const edm::Event& iEven
569     edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
570    }
571    
572 <  //=======================================================
591 <  // PVTX accessors
592 <  //=======================================================
593 <  //vertexCollection
572 >  // ======= PrimaryVertexCollection
573    static const reco::VertexCollection s_empty_vertexColl;
574    const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
575    edm::Handle<reco::VertexCollection> vertexCollectionHandle;
# Line 600 | Line 579 | PVStudy::analyze(const edm::Event& iEven
579    } else {
580     edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
581    }
582 <  //splitVertexCollection1
582 >  // ====== splitVertexCollection1
583    static const reco::VertexCollection s_empty_splitVertexColl1;
584    const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
585    edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
# Line 610 | Line 589 | PVStudy::analyze(const edm::Event& iEven
589    } else {
590      edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
591    }
592 <  //splitVertexCollection2
592 >  // ====== splitVertexCollection2
593    static const reco::VertexCollection s_empty_splitVertexColl2;
594    const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
595    edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
# Line 622 | Line 601 | PVStudy::analyze(const edm::Event& iEven
601    }
602  
603    
604 <  //=======================================================
626 <  // GET pixelVertices
627 <  //=======================================================
628 <  static const reco::VertexCollection s_empty_pixelVertexColl;
629 <  const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
630 <  edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
631 <  iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
632 <  if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
633 <    pixelVertexColl = pixelVertexCollectionHandle.product();
634 <  } else {
635 <    edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
636 <  }
637 <
638 <  //=======================================================
639 <  // BeamSpot accessors
640 <  //=======================================================
641 <
604 >  // ======== BeamSpot accessors
605    edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
606    iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
607    reco::BeamSpot bs = *recoBeamSpotHandle;    
608    const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
609    
610 <  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
610 >  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
611    
612 <  //=======================================================
650 <  // MC simvtx accessor
651 <  //=======================================================
612 >  // ========== MC simvtx accessor
613    if (!realData_) {
614      edm::Handle<SimVertexContainer> simVtxs;
615      iEvent.getByLabel( simG4_, simVtxs);
# Line 657 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618      iEvent.getByLabel( simG4_, simTrks);
619    }
620  
621 <  //=======================================================
661 <  // GET PDT
662 <  //=======================================================
621 >  // ========== GET PDT
622    try{
623      iSetup.getData(pdt);
624    }catch(const Exception&){
625      edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626    }
627 <
628 <  //=======================================================
629 <  // Apply event cleaning for firstcoll data and MC
630 <  //=======================================================
631 <  if(applyEvtClean_) {
632 <    // =====Select on the trigger bits
633 <    edm::ESHandle<L1GtTriggerMenu> menuRcd;
634 <    iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
635 <    const L1GtTriggerMenu* menu = menuRcd.product();
636 <    edm::Handle< L1GlobalTriggerReadoutRecord > gtRecord;
637 <    iEvent.getByLabel( edm::InputTag("gtDigis"), gtRecord);
638 <    
680 <    bool isTechBit0 = false;
681 <    bool isTechBit40 = false;
682 <    bool isBeamHalo = false;
683 <
684 <    TechnicalTriggerWord tw = gtRecord->technicalTriggerWord();
685 <    if ( ! tw.empty() ) {
686 <      // loop over dec. bit to get total rate (no overlap)
687 <      for ( int i = 0; i < 64; ++i ) {
688 <        if ( tw[i] ) {
689 <          //cout<<"technical number "<<i<<"  "<<endl;
690 <          if (i == 0) isTechBit0  = true;
691 <          if (i < 40 && i > 35) isBeamHalo = true;   // The beamHalo bits are 36-39
692 <          if (i == 40 || i == 41) isTechBit40 = true;
693 <        }
694 <      }
695 <    } // =====End select on the trigger bits
696 <    // =====Loop over the trackColl to tet the fraction of HighPurity tracks
697 <    int nTracks = 0;
698 <    int nHighPurityTracks=0;
699 <    double fHighPurity=0;
700 <
701 <    for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
702 <      TrackRef tkref(trackCollectionHandle,i);
703 <      if( (tkref->qualityMask() & 4 ) == 4) ++nHighPurityTracks;
704 <    }
705 <    if(nTracks>0)
706 <      fHighPurity =  double(nHighPurityTracks)/double(nTracks);
707 <    if(verbose_)
708 <      cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
709 <    
710 <    // Apply the event selection for MC events
711 <    if(realData_) {
712 <      if(!isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) {
713 <        glob_runno_ = iEvent.id().run();
714 <        glob_evtno_ = iEvent.id().event();
715 <        glob_ls_   = iEvent.luminosityBlock();
716 <        glob_bx_  = iEvent.bunchCrossing();  
717 <        return;
718 <      }
719 <      else {
720 <        // TechBit 0, beamHalo are not simulted in MC
721 <        if ( !isTechBit0 || !isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) return;
722 <      }
627 >  
628 >  // ======= Get MC information if needed
629 >  bool MC=false;  
630 >  Handle<HepMCProduct> evtMC;
631 >  if (!realData_) {
632 >    iEvent.getByLabel("generator",evtMC);
633 >    if (!evtMC.isValid()) {
634 >      MC=false;
635 >      edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
636 >    } else {
637 >      edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
638 >      MC=true;
639      }
640 <  } // End of Apply event cleaning cuts for firstcoll data
725 <
640 >  }
641  
642 +  //========================================================================
643 +  // Step 1:  Apply event cleaning for data and MC
644 +  //          WARNING: event selection cut are hard coded!!
645 +  //========================================================================
646 +  // =====Loop over the trackColl to get the fraction of HighPurity tracks
647 +  int nTracks = 0;
648 +  int nHighPurityTracks=0;
649 +  double fHighPurity=0;
650 +
651 +  for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
652 +    TrackRef tkref(trackCollectionHandle,i);
653 +    if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
654 +  }
655 +  if(nTracks>0)
656 +    fHighPurity =  double(nHighPurityTracks)/double(nTracks);
657 +  if(verbose_)
658 +    cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
659 +    
660 +  if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
661 +    glob_runno_ = iEvent.id().run();
662 +    glob_evtno_ = iEvent.id().event();
663 +    glob_ls_   = iEvent.luminosityBlock();
664 +    glob_bx_  = iEvent.bunchCrossing();  
665 +    return;
666 +  }
667  
668 <  //=======================================================
669 <  // Fill trackparameters of the input tracks to pvtx fitter
670 <  //=======================================================
668 >  //========================================================================
669 >  // Step 2: Fill histograms for the splitting consistency checks
670 >  //========================================================================
671 >  
672 >  // === Fill trackparameters of the input tracks to pvtx fitter
673    edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
674    //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
675    // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
# Line 735 | Line 677 | PVStudy::analyze(const edm::Event& iEven
677    fillTrackHisto(splitTrackColl1, 1, beamSpot);
678    fillTrackHisto(splitTrackColl2, 2, beamSpot);
679    edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
680 <
681 <
682 <  //=======================================================
741 <  // Fill pixelVertices related histograms
742 <  //=======================================================
743 <  nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
744 <  if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
745 <    //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
746 <    fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
747 <    h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());    
748 <    // Get the dZ error of the tracks in the leading pixelVertexColl
749 <    for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
750 <        t!= (pixelVertexColl->begin())->tracks_end(); t++) {
751 <      
752 <      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
753 <        //        h_pvtrk->Fill1d("trkPtPV", 0.);
754 <      }
755 <      else
756 <        h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
757 <    }
758 <    // Fill the track->dz() in the PV difference from first pixelpvtx
759 <    for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
760 <        t!= (vertexColl->begin())->tracks_end(); t++) {
761 <      // illegal charge
762 <      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
763 <        //        h_pvtrk->Fill1d("trkPtPV", 0.);
764 <      }
765 <      else {
766 <        if(pixelVertexColl->size()>0) {
767 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
768 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
769 <          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
770 <          h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
771 <        }
772 <      }
773 <    }
774 <  }
775 <
776 <  //=======================================================
777 <  // Fill number of reconstructed vertices
778 <  //=======================================================
779 <
680 >  
681 >  
682 >  // ==== Fill number of reconstructed vertices
683    edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
684    edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
685    edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
# Line 786 | Line 689 | PVStudy::analyze(const edm::Event& iEven
689      printRecVtxs(splitVertexCollection1Handle);
690      printRecVtxs(splitVertexCollection2Handle);
691    }
692 +
693    nrecPV_ = int (vertexColl->size());
694    nrecPV1_spl_ = int (splitVertexColl1->size());
695    nrecPV2_spl_ = int (splitVertexColl2->size());
# Line 796 | Line 700 | PVStudy::analyze(const edm::Event& iEven
700    h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
701    
702    
703 < //=================================================================
704 <  // Fill track parameter ntuple/hist for tracks used in recoVertices
801 <  //=================================================================
802 <
803 <  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
703 >  // ======= Fill track parameter ntuple/hist for tracks used in recoVertices
704 >  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
705    if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
706 <    fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
806 <  
807 <  if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
808 <    fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
809 <
810 <  if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
811 <    fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
812 <
813 <  //=======================================================
814 <  // Compare offlinePrimaryVertices with pixelVertices
815 <  //=======================================================
816 <  if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
817 <      && (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {    
818 <    h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
819 <    // difference in reconstructed position of the leading pvtx
820 <    //edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;  
821 <    //edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
822 <    h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
823 <    h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
824 <    h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
825 <  }
826 <  
827 <  
828 < //==========================================================
829 <  // Fill secondary/primary min separations for multi-vertices
830 <  //==========================================================
706 >    fillTrackHistoInPV(vertexColl, 0, beamSpot);
707    
708 +  // ======= Fill secondary/primary min separations for multi-vertices
709    if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {  
833
710      double min_xsep = 9999.0;
711      double min_ysep = 9999.0;
712      double min_zsep = 9999.0;    
# Line 887 | Line 763 | PVStudy::analyze(const edm::Event& iEven
763      min_zsep_ = min_zsep;
764      min_ntrksep_ = min_ntrksep;
765      min_sumpt2sep_ = min_sumpt2sep;
890
891
766    } // End of  if(vertexColl->size()>1) {
767    
768    edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
769  
770 <  
771 <    
772 <  //=======================================================
773 <  //           PrimaryVertex Matching
774 <  // 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
775 <  // 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
902 <  // Assume the first match is the primary vertex,
903 <  // since vertexColl are sorted in decreasing order of sum(pT)
904 <  //=======================================================
770 >  //========================================================================
771 >  // Step 3:  PrimaryVertex Matching
772 >  //   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
773 >  //   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
774 >  //   3. The first match is the primary vertex, which has largest sum(pT^2)
775 >  //========================================================================
776    
777    edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
778    reco::VertexCollection s_empty_matchedVertexColl1;
# Line 955 | Line 826 | PVStudy::analyze(const edm::Event& iEven
826      }
827    }
828    edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
829 <  
830 <  //=======================================================
831 <  //  Analyze matchedVertexColls
832 <  //=======================================================
829 >
830 >  //========================================================================
831 >  // Step 4: Analyze matchedVertexColls
832 >  //========================================================================
833    edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
834    
835 +  // ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
836 +  if(MC){
837 +    // make a list of primary vertices:
838 +    std::vector<simPrimaryVertex> simpv;
839 +    simpv=getSimPVs(evtMC,"");
840 +    //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
841 +    h_gen->Fill1d("nsimPV", simpv.size());
842 +    nsimPV_ = simpv.size();
843 +    int isimpv = 0;
844 +    for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
845 +        vsim!=simpv.end(); vsim++, isimpv++){
846 +      nsimTrkPV_[isimpv]  =vsim->nGenTrk;
847 +      simx_[isimpv] = vsim->x;
848 +      simy_[isimpv] = vsim->y;
849 +      simz_[isimpv] = vsim->y;
850 +      simptsq_[isimpv] = vsim->ptsq;
851 +      //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
852 +      fillMCHisto(vsim, isimpv, vertexColl, 3);
853 +    }
854 +  }
855 +  
856    // Compare the reconstructed vertex position and calculate resolution/pulls
857    if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {  
858 <    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
859 <    fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
860 <    fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
858 >    // ==== Analyze the splitted vtx using MC method
859 >    if(MC){
860 >      // make a list of primary vertices:
861 >      std::vector<simPrimaryVertex> simpv;
862 >      simpv=getSimPVs(evtMC,"");
863 >      //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
864 >      int isimpv = 0;
865 >      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
866 >          vsim!=simpv.end(); vsim++, isimpv++){
867 >        //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
868 >        fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
869 >        fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
870 >      }
871 >    }
872 >
873 >    // ==== Analyze res/pull two-vertex method
874 >    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
875 >    fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
876 >    fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
877  
878      reco::VertexCollection::const_iterator v1;
879      reco::VertexCollection::const_iterator v2;
880 +
881      nrecPV_twovtx_ = 0;
882      for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
883          v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
# Line 999 | Line 908 | PVStudy::analyze(const edm::Event& iEven
908        fillHisto(res(fres_[0], fres_[1], fres_[2]),
909                  error(ferror_[0], ferror_[1], ferror_[2]),
910                  fntrk_,0);  
911 +
912        // Fill the ntuple variables
1003      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
1004      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
1005      sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
1006      sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
913        nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
914        deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
915        deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
# Line 1014 | Line 920 | PVStudy::analyze(const edm::Event& iEven
920        pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
921        pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
922        pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];    
1017      nrecPV_twovtx_++;
1018    } // End of analyzing res/pull
923  
924 <    //=======================================================
925 <    // Fill simulated vertices
1022 <    //=======================================================
1023 <    if (!realData_) {
1024 <      bool MC=false;
1025 <      Handle<HepMCProduct> evtMC;
1026 <      iEvent.getByLabel("generator",evtMC);
1027 <      if (!evtMC.isValid()) {
1028 <        MC=false;
1029 <        edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
1030 <      } else {
1031 <        edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
1032 <        MC=true;
1033 <      }
924 >
925 >      //SplittedVertex
926        
927 <      if(MC){
928 <        // make a list of primary vertices:
929 <        std::vector<simPrimaryVertex> simpv;
930 <        simpv=getSimPVs(evtMC,"");
931 <        //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
932 <        h_gen->Fill1d("nsimPV", simpv.size());
933 <        nsimPV_ = simpv.size();
934 <        
935 <        int isimpv = 0;
936 <        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
937 <            vsim!=simpv.end(); vsim++, isimpv++){
938 <          nsimTrkPV_[isimpv]  =vsim->nGenTrk;
939 <          simx_[isimpv] = vsim->x;
940 <          simy_[isimpv] = vsim->y;
941 <          simz_[isimpv] = vsim->y;
942 <          simptsq_[isimpv] = vsim->ptsq;
943 <          
1052 <          //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1053 <          fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1054 <          fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1055 <          fillMCHisto(vsim, isimpv, vertexColl, 3);
1056 <        }
1057 <      } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1058 <    } // End of   if (!realData_) {
927 >      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
928 >      ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
929 >      normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
930 >      avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
931 >      errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
932 >      erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
933 >      errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
934 >  
935 >      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
936 >      ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
937 >      normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
938 >      avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
939 >      errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
940 >      erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
941 >      errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
942 >  
943 >      nrecPV_twovtx_++;
944  
945 +
946 +      
947 +      // Print some information of the two tracks events
948 +      if(verbose_ && fntrk_ < 4) {
949 +        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
950 +        printRecVtx(*v1);
951 +        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
952 +        printRecVtx(*v2);
953 +      }
954 +    } // End of analyzing res/pull
955    } // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
956    else
957      edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
# Line 1064 | Line 959 | PVStudy::analyze(const edm::Event& iEven
959    edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
960  
961  
962 <
963 <  // Finally fill the ftree_
962 >  //========================================================================
963 >  // Step 5: Fill ntuple if saventuple_ is on
964 >  //========================================================================
965    if(saventuple_) {
966      ftree_->Fill();
967      pvtxtree_->Fill();
# Line 1132 | Line 1028 | bool PVStudy::matchVertex(const PVStudy:
1028   {
1029    if(vrec.isFake() || !vrec.isValid()) return false;
1030    else
1031 <    return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1031 >    return true;
1032 >    //return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1033   }
1034  
1035   // Match two reconstructed vertices
# Line 1199 | Line 1096 | void PVStudy::printRecVtxs(const Handle<
1096    }
1097   }
1098  
1099 + void PVStudy::printRecVtx(const reco::Vertex & v){
1100 +
1101 +  cout << "#trk " << std::setw(3) << v.tracksSize()
1102 +         << " chi2 " << std::setw(4) << v.chi2()
1103 +         << " ndof " << std::setw(3) << v.ndof() << endl
1104 +         << " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1105 +         << " dx " << std::setw(8) << v.xError()<< endl
1106 +         << " y "  << std::setw(8) << v.y()
1107 +         << " dy " << std::setw(8) << v.yError()<< endl
1108 +         << " z "  << std::setw(8) << v.z()
1109 +         << " dz " << std::setw(8) << v.zError()
1110 +         << endl;
1111 + }
1112 +
1113   void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1114    cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1115    int i=1;
# Line 1370 | Line 1281 | void PVStudy::fillTrackHisto(const reco:
1281    }
1282   }
1283  
1284 < void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1284 > void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1285    string suffix;
1286    suffix = ""; // for unsplit tracks
1287    if(datatype == 1) suffix = "1_spl"; // for splittracks 1
# Line 1378 | Line 1289 | void PVStudy::fillTrackHistoInPV(const r
1289    int ivtx = 0;
1290    for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1291        v!=vertexColl->end(); ++v, ++ivtx) {  
1292 <    if(v->isFake()) continue;
1293 <    if(fillNtuple) {
1294 <      if(datatype == 0) {
1295 <        nTrkPV_[ivtx] = v->tracksSize();        
1296 <        sumptsq_[ivtx] = sumPtSquared(*v);
1297 <        isValid_[ivtx] = v->isValid();
1298 <        isFake_[ivtx] = v->isFake();    
1299 <        recx_[ivtx] = v->x();
1300 <        recy_[ivtx] = v->y();
1301 <        recz_[ivtx] = v->z();
1302 <        recx_err_[ivtx] = v->xError();
1303 <        recy_err_[ivtx] = v->yError();
1304 <        recz_err_[ivtx] = v->zError();
1305 <      }
1306 <      if(datatype == 1) {
1307 <        nTrkPV1_spl_[ivtx] = v->tracksSize();  
1308 <        sumptsq1_spl_[ivtx] = sumPtSquared(*v);        
1309 <        isValid1_spl_[ivtx] = v->isValid();
1310 <        isFake1_spl_[ivtx] = v->isFake();      
1311 <        recx1_spl_[ivtx] = v->x();
1312 <        recy1_spl_[ivtx] = v->y();
1402 <        recz1_spl_[ivtx] = v->z();
1403 <        recx1_err_spl_[ivtx] = v->xError();
1404 <        recy1_err_spl_[ivtx] = v->yError();
1405 <        recz1_err_spl_[ivtx] = v->zError();
1406 <
1407 <      }
1408 <      if(datatype == 2) {
1409 <        nTrkPV2_spl_[ivtx] = v->tracksSize();  
1410 <        sumptsq2_spl_[ivtx] = sumPtSquared(*v);        
1411 <        isValid2_spl_[ivtx] = v->isValid();
1412 <        isFake2_spl_[ivtx] = v->isFake();
1413 <        recx2_spl_[ivtx] = v->x();
1414 <        recy2_spl_[ivtx] = v->y();
1415 <        recz2_spl_[ivtx] = v->z();
1416 <        recx2_err_spl_[ivtx] = v->xError();
1417 <        recy2_err_spl_[ivtx] = v->yError();
1418 <        recz2_err_spl_[ivtx] = v->zError();
1419 <      }
1420 <      if(datatype == 4) { // for pixelVertices
1421 <        nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();        
1422 <        sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);      
1423 <        isValid_pxlpvtx_[ivtx] = v->isValid();
1424 <        isFake_pxlpvtx_[ivtx] = v->isFake();
1425 <        recx_pxlpvtx_[ivtx] = v->x();
1426 <        recy_pxlpvtx_[ivtx] = v->y();
1427 <        recz_pxlpvtx_[ivtx] = v->z();
1428 <        recx_err_pxlpvtx_[ivtx] = v->xError();
1429 <        recy_err_pxlpvtx_[ivtx] = v->yError();
1430 <        recz_err_pxlpvtx_[ivtx] = v->zError();
1431 <      }
1432 <    }
1433 <  }
1434 <  // For histogram only fill the leading pvtx parameters
1435 <  if(fillHisto) {
1436 <    h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1437 <    h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1438 <    int nHWTrkPV_ = 0;
1439 <    
1440 <    try {
1441 <      for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1442 <          t!=vertexColl->begin()->tracks_end(); t++) {
1443 <        // illegal charge
1444 <        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1445 <          //      h_pvtrk->Fill1d("trkPtPV", 0.);
1446 <        }
1447 <        else {
1448 <          h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1449 <          h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());  
1450 <          h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1451 <          h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1452 <          h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1453 <          h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1454 <          
1455 <          if(vertexColl->begin()->trackWeight(*t) > 0.5)
1456 <            nHWTrkPV_++;
1292 >    if(!v->isFake()) {
1293 >      h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1294 >      h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1295 >      int nHWTrkPV_ = 0;
1296 >      try {
1297 >        for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1298 >            t!=vertexColl->begin()->tracks_end(); t++) {
1299 >          // illegal charge
1300 >          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1301 >            //    h_pvtrk->Fill1d("trkPtPV", 0.);
1302 >          }
1303 >          else {
1304 >            h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1305 >            h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());  
1306 >            h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1307 >            h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1308 >            h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1309 >            h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1310 >            if(vertexColl->begin()->trackWeight(*t) > 0.5)
1311 >              nHWTrkPV_++;
1312 >          }
1313          }
1314        }
1315 +      catch (...) {
1316 +        // exception thrown when trying to use linked track
1317 +        //        h_pvtrk->Fill1d("trkPtPV", 0.);
1318 +      }
1319 +      h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1320      }
1460    catch (...) {
1461      // exception thrown when trying to use linked track
1462      //          h_pvtrk->Fill1d("trkPtPV", 0.);
1463    }
1464    h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1321    }
1466  
1467  
1322   }
1323  
1324   void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
# Line 1490 | Line 1344 | void PVStudy::fillMCHisto(std::vector<si
1344  
1345    for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1346        vrec!=vtxColl->end(); ++vrec){
1493    int nrectrk = vrec->tracksSize();
1347      vsim->recVtx=NULL;  
1495  
1348      edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1349      edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1350  
# Line 1508 | Line 1360 | void PVStudy::fillMCHisto(std::vector<si
1360        
1361        double fres_mct[3];
1362        double ferror_mct[3];
1363 <      
1363 >      int nrectrk =  int(vsim->recVtx->tracksSize());
1364 >      double ndofPV = vsim->recVtx->ndof();
1365 >      double normchi2PV =  vsim->recVtx->normalizedChi2();
1366 >      double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1367 >
1368        fres_mct[0] = vsim->recVtx->x()-vsim->x;
1369        fres_mct[1] = vsim->recVtx->y()-vsim->y;
1370        fres_mct[2] = vsim->recVtx->z()-vsim->z;
# Line 1534 | Line 1390 | void PVStudy::fillMCHisto(std::vector<si
1390                  nrectrk, datatype);
1391        
1392        if(saventuple_) {
1393 <        //Fill the values for variables in ftree_  
1393 >        //Fill the values for variables in pvtxtree_  
1394          if(datatype == 1) {
1395            nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1396 +          ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1397 +          normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1398 +          avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1399            deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1400 <          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1401 <          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1400 >          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1401 >          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1402            pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];  
1403            pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1404            pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1405 +          errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1406 +          erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1407 +          errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1408            nrecPV_spl1_mct_++;
1409          }
1548        
1410          if(datatype == 2) {
1411 <          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1411 >          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;  
1412 >          ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1413 >          normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1414 >          avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1415            deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1416 <          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1417 <          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1416 >          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1417 >          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1418            pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];  
1419            pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1420 <          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1420 >          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];        
1421 >          errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1422 >          erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1423 >          errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1424            nrecPV_spl2_mct_++;
1425          }
1426          if(datatype == 3) {    
1427 <          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1427 >          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1428 >          ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1429 >          normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1430 >          avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1431            deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1432            deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1433            deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1434            pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];  
1435            pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1436 <          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];          
1436 >          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1437 >          errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1438 >          erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1439 >          errz_mct_[nrecPV_mct_] =  ferror_mct[2];        
1440            nrecPV_mct_++;
1441          }
1442        } // End of  if(saventuple_) {
# Line 1575 | Line 1448 | void PVStudy::fillMCHisto(std::vector<si
1448   }
1449  
1450   void PVStudy::SetVarToZero() {
1578  //Greg's variables
1451    fntrk_ = 0;
1452    //pvtx position (x,y,z) residual and error
1453    for(int i = 0; i<3;i++)    {
# Line 1583 | Line 1455 | void PVStudy::SetVarToZero() {
1455      ferror_[i] = 0;
1456    }
1457    
1458 <  //Yanyan's variables
1459 <  // Number of reconstructed vertices  
1458 >  // Event level information
1459 >  glob_runno_ = 0;
1460 >  glob_evtno_ = 0;
1461 >  glob_ls_ = 0;
1462 >  glob_bx_ = 0;
1463    nrecPV_ = 0;
1464    nrecPV1_spl_ = 0;
1465    nrecPV2_spl_ = 0;
# Line 1600 | Line 1475 | void PVStudy::SetVarToZero() {
1475    min_ntrksep_ = 9999.0;
1476    min_sumpt2sep_ = -9999.0;
1477    
1478 <
1478 >  // Variables filled per Vertex
1479    for (int i = 0; i < nMaxPVs_; i++) {
1605    // recoVertices with all tracks
1606    nTrkPV_[i] = 0; // Number of tracks in the pvtx    
1607    sumptsq_[i] = 0;
1608    isValid_[i] = -1;
1609    isFake_[i] = -1;
1610    recx_[i] = 0;
1611    recy_[i] = 0;
1612    recz_[i] = 0;
1613    recx_err_[i] = 0;
1614    recy_err_[i] = 0;
1615    recz_err_[i] = 0;
1616    
1617    // recoVertices with splitTrack1
1618    nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx    
1619    sumptsq1_spl_[i] = 0;
1620    isValid1_spl_[i] = -1;
1621    isFake1_spl_[i] = -1;
1622    recx1_spl_[i] = 0;
1623    recy1_spl_[i] = 0;
1624    recz1_spl_[i] = 0;
1625    recx1_err_spl_[i] = 0;
1626    recy1_err_spl_[i] = 0;
1627    recz1_err_spl_[i] = 0;
1628  
1629    // recoVertices with splitTrack2
1630    nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx  
1631    sumptsq2_spl_[i] = 0;
1632    isValid2_spl_[i] = -1;
1633    isFake2_spl_[i] = -1;
1634    recx2_spl_[i] = 0;
1635    recy2_spl_[i] = 0;
1636    recz2_spl_[i] = 0;
1637    recx2_err_spl_[i] = 0;
1638    recy2_err_spl_[i] = 0;
1639    recz2_err_spl_[i] = 0;
1640    
1641    //pixelVertices
1642    nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx  
1643    sumptsq_pxlpvtx_[i] = 0;
1644    isValid_pxlpvtx_[i] = -1;
1645    isFake_pxlpvtx_[i] = -1;
1646    recx_pxlpvtx_[i] = 0;
1647    recy_pxlpvtx_[i] = 0;
1648    recz_pxlpvtx_[i] = 0;
1649    recx_err_pxlpvtx_[i] = 0;
1650    recy_err_pxlpvtx_[i] = 0;
1651    recz_err_pxlpvtx_[i] = 0;
1652
1480      // matched two-vertices
1481      nTrkPV1_spl_twovtx_[i] = 0;
1482 +    ndofPV1_spl_twovtx_[i] = 0;
1483 +    normchi2PV1_spl_twovtx_[i] = 0;
1484 +    avgPtPV1_spl_twovtx_[i] = 0;
1485 +    errx1_spl_twovtx_[i] = 0;
1486 +    erry1_spl_twovtx_[i] = 0;
1487 +    errz1_spl_twovtx_[i] = 0;
1488 +
1489      nTrkPV2_spl_twovtx_[i] = 0;
1490 <    nTrkPV_twovtx_[i] = 0;  
1491 <    sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1492 <    sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1490 >    ndofPV2_spl_twovtx_[i] = 0;
1491 >    normchi2PV2_spl_twovtx_[i] = 0;
1492 >    avgPtPV2_spl_twovtx_[i] = 0;
1493 >    errx2_spl_twovtx_[i] = 0;
1494 >    erry2_spl_twovtx_[i] = 0;
1495 >    errz2_spl_twovtx_[i] = 0;
1496 >    
1497      deltax_twovtx_[i] = 0;
1498      deltay_twovtx_[i] = 0;
1499      deltaz_twovtx_[i] = 0;
# Line 1680 | Line 1518 | void PVStudy::SetVarToZero() {
1518      pullx_mct_[i] = 0;
1519      pully_mct_[i] = 0;
1520      pullz_mct_[i] = 0;
1521 <        
1521 >    nTrkPV_mct_[i] = 0;
1522 >    ndofPV_mct_[i] = 0;
1523 >    normchi2PV_mct_[i] = 0;
1524 >    avgPtPV_mct_[i] = 0;
1525 >    errx_mct_[i] = 0;
1526 >    erry_mct_[i] = 0;
1527 >    errz_mct_[i] = 0;
1528 >    
1529      deltax_spl1_mct_[i] = 0;
1530      deltay_spl1_mct_[i] = 0;
1531      deltaz_spl1_mct_[i] = 0;
1532      pullx_spl1_mct_[i] = 0;
1533      pully_spl1_mct_[i] = 0;
1534      pullz_spl1_mct_[i] = 0;
1535 <    
1535 >    nTrkPV_spl1_mct_[i] = 0;
1536 >    ndofPV_spl1_mct_[i] = 0;
1537 >    normchi2PV_spl1_mct_[i] = 0;
1538 >    avgPtPV_spl1_mct_[i] = 0;
1539 >    errx_spl1_mct_[i] = 0;
1540 >    erry_spl1_mct_[i] = 0;
1541 >    errz_spl1_mct_[i] = 0;
1542 >
1543      deltax_spl2_mct_[i] = 0;
1544      deltay_spl2_mct_[i] = 0;
1545      deltaz_spl2_mct_[i] = 0;
1546      pullx_spl2_mct_[i] = 0;
1547      pully_spl2_mct_[i] = 0;
1548      pullz_spl2_mct_[i] = 0;
1549 +    nTrkPV_spl2_mct_[i] = 0;
1550 +    ndofPV_spl2_mct_[i] = 0;
1551 +    normchi2PV_spl2_mct_[i] = 0;
1552 +    avgPtPV_spl2_mct_[i] = 0;
1553 +    errx_spl2_mct_[i] = 0;
1554 +    erry_spl2_mct_[i] = 0;
1555 +    errz_spl2_mct_[i] = 0;
1556 +
1557    }
1558   }
1559  
1560   double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1561 <  double sum = 0.;
1562 <  double pT;
1563 <  for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1564 <    pT = (**it).pt();
1565 <    sum += pT*pT;
1561 >  if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1562 >  else {
1563 >    double sum = 0.;
1564 >    double pT;
1565 >    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1566 >      pT = (**it).pt();
1567 >      sum += pT*pT;
1568 >    }
1569 >    return sum;
1570 >  }
1571 > }
1572 >
1573 >
1574 >
1575 > double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1576 >
1577 >  if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1578 >  else {
1579 >    double sumpT = 0.;
1580 >    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1581 >      sumpT +=  (**it).pt();
1582 >  }
1583 >    return sumpT/double(v.tracksSize());
1584    }
1707  return sum;
1585   }
1586  
1587  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines