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.8 by jengbou, Fri Oct 16 00:11:30 2009 UTC vs.
Revision 1.22 by yygao, Mon Mar 29 15:27:24 2010 UTC

# Line 38 | 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 49 | 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 59 | 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 76 | 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 <
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  
87 –
96    // Specify the data mode vector
97    if(realData_) datamode.push_back(0);
98    else {
# Line 93 | Line 101 | PVStudy::PVStudy(const edm::ParameterSet
101      datamode.push_back(2);
102      datamode.push_back(3);
103    }
104 <
97 <  // 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 102 | Line 109 | PVStudy::PVStudy(const edm::ParameterSet
109      ftree_->Branch("residual",&fres_,"fres_[3]/D");
110      ftree_->Branch("error",&ferror_,"ferror_[3]/D");
111      ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
105 –    ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
106 –    ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
107 –    ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
108 –  
109 –    // If it is MC, get the resA, resB
110 –    if(!realData_) {
111 –      ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
112 –      ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
113 –      ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
114 –      ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
115 –      ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
116 –      ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
117 –      ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
118 –      ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
119 –      ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
120 –    }
121 –  }
112  
113 <  edm::Service<TFileService> fs;
114 <  TFileDirectory subDir = fs->mkdir( "Summary" );
125 <  TFileDirectory subDir1 = fs->mkdir( "Others" );
126 <  //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
113 >    // pvtxtree_ analyzing the pvtxs ootb
114 >    pvtxtree_ = new TTree("pvtxtree","pvtxtree");
115  
116 <  setRootStyle();
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 >    //Fill the variables in the twovtx pair (recvtx1, recvtx2)  
132 >    // Information related to the analyzing the two-vertex method
133 >
134 >
135 >
136 >    pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
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");
140 >    pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
141 >    pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
142 >    pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
143 >    pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
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 <  //Book histograms for track and pvtx parameters for each splitted pvtx
158 <  //i=0 : x (as in unsplit track collection)
159 <  //i=1 : x1_spl_
160 <  //i=2 : x2_spl_
161 <  
162 <  for(int i=0;i<3;i++)
163 <    {  
164 <      string suffix;
165 <      if(i == 0) suffix = "";
166 <      else {
167 <        stringstream ss;
168 <        ss << i;
169 <        suffix =ss.str()+"_spl";  
170 <      }
171 <      h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
172 <      h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
173 <      h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
174 <      h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
175 <      h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);  
176 <      h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);  
177 <      
178 <      h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
179 <      h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
180 <      h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
181 <      h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
182 <      h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);  
183 <      h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
184 <      h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);    
185 <      suffix.clear();
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_) {  
170 >      pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
171 >      pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
172 >      pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
173 >      pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
174 >      pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
175 >      pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
176 >
177 >      // For pvtxs with all the tracks
178 >      pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
179 >      pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
180 >      pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
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");    
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");
196 >      pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
197 >      pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
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");
212 >      pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
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");  
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 <  h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
161 <  h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
162 <  h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
226 >  }
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
237 +  theFile->mkdir("Summary");
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();
248 +  h_misc    = new PVHistograms();
249 +  h_summary = new PVHistograms();
250 +  h_others  = new PVHistograms();
251 +
252 +  for(int i=0;i<3;i++) {  
253 +    if(i == 0) h_pvtrk->Init("pvTrk");
254 +    else {
255 +      stringstream ss;
256 +      ss << i;
257 +      h_pvtrk->Init("pvTrk", ss.str(),"spl");
258 +    }
259 +  }
260 +  h_misc->Init("misc");
261 +
262 +  // Book MC only plots
263 +  if (!realData_) {
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 175 | Line 280 | PVStudy::PVStudy(const edm::ParameterSet
280      case 3: suffix = "_mct";
281        break;
282      }
283 <    h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
179 <    h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
180 <    h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
181 <    h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
182 <    h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
183 <    h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
184 <    h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
185 <    h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
186 <    h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
283 >    h_summary->Init("summary",suffix);
284      
285      // Book residual, error and pull histograms for each nTrk bin
286      for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
287        stringstream ss;
288        ss << ntrk;
289 <      string suffix2 = suffix+"_"+ss.str();
290 <      h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
291 <      h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
195 <      h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
196 <      h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
197 <      h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
198 <      h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
199 <      h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
200 <      h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
201 <      h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
202 <      h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
203 <      h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
204 <      h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
205 <      suffix2.clear();
206 <    } // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
207 <    
289 >      h_others->Init("others",suffix,ss.str());
290 >    }
291 >
292      // Book residual and pull histograms only when analyzeOntheFly is enabled
293 <    if(analyze_) {
210 <      h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
211 <      h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
212 <      h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
213 <      h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
214 <      h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
215 <      h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
216 <      h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
217 <      h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
218 <      h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
219 <      h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
220 <      h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
221 <      h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
222 <      h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
223 <      h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
224 <      h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
225 <      h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
226 <      h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
227 <      h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
228 <      h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
229 <      h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
230 <      h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
231 <      h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
232 <      h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
233 <      h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
234 <      h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
235 <      h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
236 <      h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
237 <      h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
238 <      h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
239 <      h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
240 <    } // End of   if(analyze_) {
293 >    if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
294      suffix.clear();
295 <  } // End of Book histograms sensitive to data/mc  
296 <  // Book MC only plots
297 <  if (!realData_) {  
245 <    h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
246 <    h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
247 <    h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
248 <    h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");    
249 <    h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);  
250 <  }
295 >  } // End of Book histograms sensitive to data/mc    
296 >  
297 >
298   }
299  
300   PVStudy::~PVStudy()
# Line 255 | Line 302 | PVStudy::~PVStudy()
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();
308 +  h_misc->Save();
309 +  h_summary->Save();
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();
320  
321 +  theFile->Close();
322   }
323  
324   void PVStudy::setRootStyle() {
# Line 279 | 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 292 | 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;
302 <      cout << "number of vertices " << evt->vertices_size() << endl;
303 <    }
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 311 | 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["genPart_cT"]->Fill(pos.t()); // mm
381 <            h["genPart_T"]->Fill(pos.t()/299.792458); // ns
380 >            h_gen->Fill1d("genPart_cT", pos.t()); // mm
381 >            h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
382            }
383            npv++;
384          }
# Line 325 | 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 356 | 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 384 | Line 445 | std::vector<PVStudy::simPrimaryVertex> P
445              }
446            }
447          }
448 <      }
449 <    }
448 >      }//loop MC vertices daughters
449 >    }//loop MC vertices
450    }
451    return simpv;
452   }
# Line 418 | Line 479 | std::vector<PVStudy::simPrimaryVertex> P
479                ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
480              }
481            }
482 <          primary =  ( gv->position().t()==0);
482 >          //primary =  ( gv->position().t()==0 );
483 >          primary = true;
484 >          h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
485 >          h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
486 >
487            //resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
488            // no more mother pointer in the improved HepMC GenParticle
489            resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
# Line 467 | 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
541 <  //=======================================================
542 <  //if(saventuple_)
543 <    SetVarToZero();
477 <  
478 <  //=======================================================
479 <  // Track accessors
480 <  //=======================================================
481 <  //trackCollection
541 >  SetVarToZero();
542 >  
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 486 | 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 <  //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 496 | 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 <  //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 506 | 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 <  //=======================================================
513 <  // PVTX accessors
514 <  //=======================================================
515 <  //vertexCollection
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 520 | 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 <  //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 530 | 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;
592 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593    }
594 <  //splitVertexCollection2
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 540 | 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  
546 –  if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
605    
606 <  //=======================================================
607 <  // MC simvtx accessor
608 <  //=======================================================
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 >  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
613 >  
614 >  // ========== MC simvtx accessor
615    if (!realData_) {
616      edm::Handle<SimVertexContainer> simVtxs;
617      iEvent.getByLabel( simG4_, simVtxs);
# Line 556 | Line 620 | PVStudy::analyze(const edm::Event& iEven
620      iEvent.getByLabel( simG4_, simTrks);
621    }
622  
623 <  //=======================================================
560 <  // GET PDT
561 <  //=======================================================
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 >  // ======= 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    }
643  
644 <  //=======================================================
645 <  // Fill trackparameters of the input tracks to pvtx fitter
646 <  //=======================================================
647 <
648 <  // GeneralTracks  
649 <  h["nTrk"]->Fill(trackColl->size());
650 <  for(TrackCollection::const_iterator itTrack = trackColl->begin();
651 <      itTrack != trackColl->end();                      
652 <      ++itTrack) {
653 <    h["trkPt"]->Fill(itTrack->pt());
654 <    h["trkDxy"]->Fill(itTrack->dxy());
655 <    h["trkDz"]->Fill(itTrack->dz());
580 <    h["trkEta"]->Fill(itTrack->eta());
581 <    h["trkPhi"]->Fill(itTrack->phi());
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 <  //SplittedTracks1
665 <  h["nTrk1_spl"]->Fill(splitTrackColl1->size());
666 <  for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
667 <      itTrack != splitTrackColl1->end();                      
588 <      ++itTrack) {
589 <    h["trkPt1_spl"]->Fill(itTrack->pt());
590 <    h["trkDxy1_spl"]->Fill(itTrack->dxy());
591 <    h["trkDz1_spl"]->Fill(itTrack->dz());
592 <    h["trkEta1_spl"]->Fill(itTrack->eta());
593 <    h["trkPhi1_spl"]->Fill(itTrack->phi());
594 <  }
595 <  //SplittedTracks2  
596 <  h["nTrk2_spl"]->Fill(splitTrackColl2->size());
597 <  for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
598 <      itTrack != splitTrackColl2->end();                      
599 <      ++itTrack) {
600 <    h["trkPt2_spl"]->Fill(itTrack->pt());
601 <    h["trkDxy2_spl"]->Fill(itTrack->dxy());
602 <    h["trkDz2_spl"]->Fill(itTrack->dz());
603 <    h["trkEta2_spl"]->Fill(itTrack->eta());
604 <    h["trkPhi2_spl"]->Fill(itTrack->phi());
605 <  }
664 >  glob_runno_ = iEvent.id().run();
665 >  glob_evtno_ = iEvent.id().event();
666 >  glob_ls_   = iEvent.luminosityBlock();
667 >  glob_bx_  = iEvent.bunchCrossing();  
668    
669 <  if(verbose_)
670 <    cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
671 <
672 <  //=======================================================
673 <  // Fill reconstructed vertices
674 <  //=======================================================
675 <  if(verbose_)  {
676 <    cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
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);
616 –    cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
690      printRecVtxs(splitVertexCollection1Handle);
618 –    cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
691      printRecVtxs(splitVertexCollection2Handle);
692    }
693 +
694 +  nrecPV_ = int (vertexColl->size());
695 +  nrecPV1_spl_ = int (splitVertexColl1->size());
696 +  nrecPV2_spl_ = int (splitVertexColl2->size());
697 +
698 +  h_pvtrk->Fill1d("nrecPV", nrecPV_);
699 +  h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
700 +  h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
701 +  h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
702    
622 –  //Fill the tree variables
623 –  nrecPV_ = vertexColl->size();
624 –  nrecPV1_spl_ = splitVertexColl1->size();
625 –  nrecPV2_spl_ = splitVertexColl2->size();
626 –
627 –  h["nrecPV"]->Fill(nrecPV_);
628 –  h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
629 –  h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
630 –  h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
703    
704 <  // Fill the track paramters for the vertexColl
705 <  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
706 <      v!=vertexColl->end(); ++v){
707 <    h["nTrkPV"]->Fill(v->tracksSize());
708 <    try {
709 <      for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
710 <          t!=v->tracks_end(); t++) {
711 <        // illegal charge
712 <        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
713 <          //      h["trkPtPV"]->Fill(0.);
714 <        }
715 <        else {
716 <          h["trkPtPV"]->Fill((**t).pt());
717 <          h["trkDxyPV"]->Fill((**t).dxy());
718 <          h["trkDzPV"]->Fill((**t).dz());
719 <          h["trkEtaPV"]->Fill((**t).eta());
720 <          h["trkPhiPV"]->Fill((**t).phi());
721 <        }
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()) ) {  
711 >    double min_xsep = 9999.0;
712 >    double min_ysep = 9999.0;
713 >    double min_zsep = 9999.0;    
714 >    double min_xsep_sign = 9999.0;
715 >    double min_ysep_sign = 9999.0;
716 >    double min_zsep_sign = 9999.0;
717 >    double min_ntrksep = 9999.0;
718 >    double min_sumpt2sep = 9999999.0;
719 >    
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 ;
724 >        v!=vertexColl->end(); ++v) {  
725 >      if(v->isValid() && ! v->isFake()) {
726 >        // xsep
727 >        if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
728 >          min_xsep = fabs(v->x()- vertexColl->begin()->x());
729 >        // ysep
730 >        if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
731 >          min_ysep = fabs(v->y()- vertexColl->begin()->y());
732 >        // zsep
733 >        if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
734 >          min_zsep = fabs(v->z()- vertexColl->begin()->z());
735 >        // xsep_sign
736 >        double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
737 >        if(xsep_sign < min_xsep_sign)
738 >          min_xsep_sign = xsep_sign;
739 >        // ysep_sign
740 >        double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
741 >        if(ysep_sign < min_ysep_sign)
742 >          min_ysep_sign = ysep_sign;
743 >        // zsep_sign
744 >        double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
745 >        if(zsep_sign < min_zsep_sign)
746 >          min_zsep_sign = zsep_sign;
747 >        // ntrksep
748 >        if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
749 >          min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
750 >        // sumpt2sep
751 >        if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
752 >          min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
753        }
754      }
755 <    catch (...) {
756 <      // exception thrown when trying to use linked track
757 <      //          h["trkPtPV"]->Fill(0.);
758 <    }
759 <  }
755 >    h_misc->Fill1d("min_xsep", min_xsep);
756 >    h_misc->Fill1d("min_ysep", min_ysep);
757 >    h_misc->Fill1d("min_zsep", min_zsep);
758 >    h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
759 >    h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
760 >    h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
761 >    h_misc->Fill1d("min_ntrksep", min_ntrksep);
762 >    h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
763 >
764 >    min_zsep_ = min_zsep;
765 >    min_ntrksep_ = min_ntrksep;
766 >    min_sumpt2sep_ = min_sumpt2sep;
767 >  } // End of  if(vertexColl->size()>1) {
768 >  
769 >  edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
770 >
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());
782 +  reco::VertexCollection s_empty_matchedVertexColl2;
783 +  reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
784 +  matchedVertexColl2->reserve(splitVertexColl2->size());
785    
786 <  // == ignore multi-vertices (for now)
787 <  //take single pvt collection and compare
788 <  if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
789 <    const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
790 <    const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
791 <    if (matchVertex(recvtx1, recvtx2))  {
792 <      int nTrkPV1 = recvtx1.tracksSize();
793 <      int nTrkPV2 = recvtx2.tracksSize();    
794 <      h["nTrkPV1_spl"]->Fill(nTrkPV1);
795 <      h["nTrkPV2_spl"]->Fill(nTrkPV2);
796 <      double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
797 <      h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
798 <      h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
786 >  bool stopmatching = false;
787 >  for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
788 >    reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
789 >    if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
790 >    for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {    
791 >      //------------------------------------------------------------------------
792 >      //== If only considering matching the first vertex of the splitVertexColl
793 >      //------------------------------------------------------------------------
794 >      if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
795 >      reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
796 >      if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
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 <      // Fill the track paramters for the splitVertexColl1
802 <      try {
803 <        for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
804 <            t!=recvtx1.tracks_end(); t++) {
805 <          // illegal charge
806 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
807 <            //    h["trkPtPV"]->Fill(0.);
808 <          }
809 <          else {
682 <            h["trkPtPV1_spl"]->Fill((**t).pt());
683 <            h["trkDxyPV1_spl"]->Fill((**t).dxy());
684 <            h["trkDzPV1_spl"]->Fill((**t).dz());
685 <            h["trkEtaPV1_spl"]->Fill((**t).eta());
686 <            h["trkPhiPV1_spl"]->Fill((**t).phi());
687 <          }
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 >        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 <      }
812 <      catch (...) {
813 <        // exception thrown when trying to use linked track
814 <        //        h["trkPtPV"]->Fill(0.);
815 <      }
816 <      
817 <      // Fill the track paramters for the splitVertexColl2
818 <      try {
819 <        for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
820 <            t!=recvtx2.tracks_end(); t++) {
821 <          // illegal charge
822 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
823 <            //    h["trkPtPV"]->Fill(0.);
824 <          }
703 <          else {
704 <            h["trkPtPV2_spl"]->Fill((**t).pt());
705 <            h["trkDxyPV2_spl"]->Fill((**t).dxy());
706 <            h["trkDzPV2_spl"]->Fill((**t).dz());
707 <            h["trkEtaPV2_spl"]->Fill((**t).eta());
708 <            h["trkPhiPV2_spl"]->Fill((**t).phi());
709 <          }
811 >        else {
812 >          nTrkPV1 = recvtx1->tracksSize();
813 >          nTrkPV2 = recvtx2->tracksSize();
814 >        }
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 >          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 +          edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
828        }
829 <      catch (...) {
830 <        // exception thrown when trying to use linked track
831 <        //        h["trkPtPV"]->Fill(0.);
829 >      else {
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        }
716 –      
717 –      if(verbose_)
718 –        cout<<"Done filling track parameters in reconstruced vertices."<<endl;
719 –
720 –      if(fabs(ntrkreldiff)<ntrkdiffcut_) {
721 –        fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
722 –        fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
723 –        fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);      
724 –        ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
725 –        ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
726 –        ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
727 –        fntrk_ =  (nTrkPV1+nTrkPV2)/2;
728 –        
729 –        h["deltax"]->Fill( fres_[0] );
730 –        h["deltay"]->Fill( fres_[1] );
731 –        h["deltaz"]->Fill( fres_[2] );
732 –        h["pullx"]->Fill( fres_[0]/ferror_[0]);
733 –        h["pully"]->Fill( fres_[1]/ferror_[1]);
734 –        h["pullz"]->Fill( fres_[2]/ferror_[2]);
735 –        h["errPVx"]->Fill( ferror_[0] );
736 –        h["errPVy"]->Fill( ferror_[1] );
737 –        h["errPVz"]->Fill( ferror_[2] );
738 –        pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
739 –                                         error(ferror_[0], ferror_[1], ferror_[2]),
740 –                                         fntrk_));
741 –        // Fill histo according to its track multiplicity, set datamode = 0
742 –        fillHisto(res(fres_[0], fres_[1], fres_[2]),
743 –                  error(ferror_[0], ferror_[1], ferror_[2]),
744 –                  fntrk_,0);  
745 –
746 –        if (saventuple_) ftree_->Fill();        
747 –      } // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
748 –    } // End of  if (matchVertex(recvtx1, recvtx2))  {
749 –    else   {
750 –      cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
751 –      cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
834      }
835 <  } // End of fill histograms for two-vertices Method
835 >  }
836 >  edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
837  
838 <  if(verbose_)
839 <    cout<<"Done filling res/pull with two-vertices method"<<endl;
840 <
841 <  //=======================================================
842 <  // Fill simulated vertices
843 <  //=======================================================
844 <  if (!realData_) {
845 <    bool MC=false;
846 <    Handle<HepMCProduct> evtMC;
847 <    iEvent.getByLabel("generator",evtMC);
848 <    if (!evtMC.isValid()) {
849 <      MC=false;
850 <      if(verbose_){
851 <        cout << "no HepMCProduct found"<< endl;
852 <      }
853 <    } else {
854 <      if(verbose_){
855 <        cout << "generator HepMCProduct found"<< endl;
856 <      }
857 <      MC=true;
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 <    
862 >  }
863 >  
864 >  // Compare the reconstructed vertex position and calculate resolution/pulls
865 >  if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {  
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 <      h["nsimPV"]->Fill(simpv.size());
783 <      int nsimtrk=0;
784 <      int nrectrk=0;
785 <      
872 >      int isimpv = 0;
873        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
874 <          vsim!=simpv.end(); vsim++){
875 <        nsimtrk+=vsim->nGenTrk;
876 <        // Loop over the datamodes, for MC, fill only 1, 2, 3
877 <        for (int i=1;i<=3;i++) {
791 <          static const reco::VertexCollection s_empty_vtxColl;
792 <          const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
793 <          std::string suffix;
794 <          // Set the vertexColl and histogram suffix according to datamode
795 <          if (i == 1) {
796 <            vtxColl = splitVertexColl1;
797 <            suffix = "_spl1_mct";
798 <          }
799 <          if (i == 2) {
800 <            vtxColl =  splitVertexColl2;
801 <            suffix = "_spl2_mct";
802 <          }
803 <          if (i == 3) {
804 <            vtxColl =  vertexColl;
805 <            suffix = "_mct";
806 <          }
807 <          //========================================================
808 <          //  look for a matching reconstructed vertex in vertexColl
809 <          //========================================================    
810 <          vsim->recVtx=NULL;
811 <          for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
812 <              vrec!=vtxColl->end(); ++vrec){
813 <            nrectrk=vrec->tracksSize();
814 <            if(verbose_){
815 <              cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
816 <              cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
817 <            }
818 <            if ( matchVertex(*vsim,*vrec) ){
819 <              // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
820 <              if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
821 <                     || (!vsim->recVtx) )
822 <                vsim->recVtx=&(*vrec);
823 <            }// End of finding the matched reco_vertex
824 <            
825 <            if (vsim->recVtx){
826 <              if(verbose_) {
827 <                cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
828 <              }
829 <              double fres_mct[3];
830 <              double ferror_mct[3];
831 <              int fntrk_mct;
832 <              
833 <              fres_mct[0] = vsim->recVtx->x()-vsim->x;
834 <              fres_mct[1] = vsim->recVtx->y()-vsim->y;
835 <              fres_mct[2] = vsim->recVtx->z()-vsim->z;
836 <              ferror_mct[0] = vsim->recVtx->xError();
837 <              ferror_mct[1] = vsim->recVtx->yError();
838 <              ferror_mct[2] = vsim->recVtx->zError();
839 <              fntrk_mct = nrectrk;
840 <              
841 <              //Fill the values for variables in ftree_
842 <              if(i == 1) {
843 <                for(int j = 0;j<3;j++)
844 <                  fres_spl1_mct_[j] = fres_mct[j];
845 <                fntrk_spl1_mct_ =  fntrk_mct;
846 <              }
847 <              if(i == 2) {
848 <                for(int j = 0;j<3;j++)
849 <                  fres_spl2_mct_[j] = fres_mct[j];
850 <                fntrk_spl2_mct_ =  fntrk_mct;
851 <              }
852 <              if(i == 3) {
853 <                for(int j=0;j<3;j++)
854 <                  fres_mct_[j] = fres_mct[j];
855 <                fntrk_mct_ =  fntrk_mct;
856 <              }
857 <              
858 <              h["deltax"+suffix]->Fill( fres_mct[0] );
859 <              h["deltay"+suffix]->Fill( fres_mct[1] );
860 <              h["deltaz"+suffix]->Fill( fres_mct[2] );
861 <              h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
862 <              h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
863 <              h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
864 <              h["errPVx"+suffix]->Fill( ferror_mct[0] );
865 <              h["errPVy"+suffix]->Fill( ferror_mct[1] );
866 <              h["errPVz"+suffix]->Fill( ferror_mct[2] );
867 <              pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
868 <                                               error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
869 <                                               fntrk_mct));
870 <              // Fill histo according to its track multiplicity
871 <              fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
872 <                        error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
873 <                        fntrk_mct, i);
874 <              if (saventuple_) ftree_->Fill();
875 <              suffix.clear();
876 <            }
877 <            else {  // no rec vertex found for this simvertex
878 <              if(verbose_)
879 <                cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
880 <            }
881 <          }
882 <        } // === End of Looping through vertexColl
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 <    } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
880 <  } // End of Filling MC variables
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) {
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);
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] );
917 >      h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
918 >      h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
919 >      h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
920 >      h_summary->Fill1d("errPVx", ferror_[0] );
921 >      h_summary->Fill1d("errPVy", ferror_[1] );
922 >      h_summary->Fill1d("errPVz", ferror_[2] );
923 >      pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
924 >                                       error(ferror_[0], ferror_[1], ferror_[2]),
925 >                                       fntrk_));
926 >      // Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
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
932 >      nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
933 >      deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
934 >      deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
935 >      deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
936 >      errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
937 >      erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
938 >      errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
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 >    edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
974 >  
975 >  edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
976 >
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 894 | 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 909 | Line 1011 | PVStudy::endJob() {
1011      case 3: suffix = "_mct";
1012        break;
1013      }
1014 <    if (analyze_) {
1014 >    suffix += "_nTrk";
1015 >    if(analyze_) {
1016        for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1017 +        edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1018          PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1019 <        if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1020 <        if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1021 <        if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1022 <        if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1023 <        if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1024 <        if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
921 <        if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(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      }
1027      suffix.clear();
# Line 926 | Line 1029 | PVStudy::endJob() {
1029    if(saventuple_) {
1030      file_->cd();
1031      ftree_->Write();
1032 +    pvtxtree_->Write();
1033      file_->Close();
1034    }
1035 < }
932 <
1035 >  
1036  
1037 + }
1038  
1039   // Match a recovertex with a simvertex
1040   // ? Should the cut be parameter dependent?
# Line 938 | 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
1052 < bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1053 <                           const reco::Vertex & vrec2)
1052 > bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1053 >                           const reco::VertexRef & vrec2,
1054 >                           int zsigncut)
1055   {
1056 <  return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1056 >  double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
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 1002 | 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 1023 | 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
1158 <  h["deltax"+suffix]->Fill(r.x());
1159 <  h["deltay"+suffix]->Fill(r.y());
1160 <  h["deltaz"+suffix]->Fill(r.z());
1161 <  h["pullx"+suffix]->Fill(r.x()/er.x());
1162 <  h["pully"+suffix]->Fill(r.y()/er.y());
1163 <  h["pullz"+suffix]->Fill(r.z()/er.z());  
1164 <  h["errPVx"+suffix]->Fill(er.x());
1165 <  h["errPVy"+suffix]->Fill(er.y());
1166 <  h["errPVz"+suffix]->Fill(er.z());
1158 >  h_others->Fill1d("deltax"+suffix, r.x());
1159 >  h_others->Fill1d("deltay"+suffix, r.y());
1160 >  h_others->Fill1d("deltaz"+suffix, r.z());
1161 >  h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1162 >  h_others->Fill1d("pully"+suffix, r.y()/er.y());
1163 >  h_others->Fill1d("pullz"+suffix, r.z()/er.z());  
1164 >  h_others->Fill1d("errPVx"+suffix, er.x());
1165 >  h_others->Fill1d("errPVy"+suffix, er.y());
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 1059 | Line 1184 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1184      for ( int i = 0; i < 6; ++i) {
1185        switch (i) {
1186        case 0:
1187 <        fit(h["deltax"+suffix], fpar);
1187 >        if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1188 >        fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1189          data[i] = fpar[0]*cm2um;
1190          break;
1191        case 1:
1192 <        fit(h["deltay"+suffix], fpar);
1192 >        if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1193 >        fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1194          data[i] = fpar[0]*cm2um;
1195          break;
1196        case 2:
1197 <        fit(h["deltaz"+suffix], fpar);
1197 >        if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1198 >        fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1199          data[i] = fpar[0]*cm2um;
1200          break;
1201        case 3:
1202 <        fit(h["pullx"+suffix], fpar);
1202 >        if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1203 >        fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1204          data[i] = fpar[0];
1205          break;
1206        case 4:
1207 <        fit(h["pully"+suffix], fpar);
1207 >        if (!h_others->ReadHisto1D("pully"+suffix)) break;
1208 >        fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1209          data[i] = fpar[0];
1210          break;
1211        case 5:
1212 <        fit(h["pullz"+suffix], fpar);
1212 >        if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1213 >        fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
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 1098 | Line 1229 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1229                              ntk);
1230   }
1231  
1232 < void PVStudy::fit(TH1 *&hdist, double fitPars[]){
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;
1116 <  }
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 1125 | Line 1257 | void PVStudy::fit(TH1 *&hdist, double fi
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 +
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
1285 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1286 +  
1287 +  h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1288 +  for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1289 +      itTrack != trackColl->end();                      
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, const Point & bs) {
1301 +  string suffix;
1302 +  suffix = ""; // for unsplit tracks
1303 +  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1304 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1305 +  int ivtx = 0;
1306 +  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1307 +      v!=vertexColl->end(); ++v, ++ivtx) {  
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 +      catch (...) {
1330 +        // exception thrown when trying to use linked track
1331 +        //        h_pvtrk->Fill1d("trkPtPV", 0.);
1332 +      }
1333 +    }
1334 +  }
1335 + }
1336 +
1337 + void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1338 + {
1339 +  std::string suffix;
1340 +  // Set the vertexColl and histogram suffix according to datamode
1341 +  if (datatype == 1) {
1342 +    suffix = "_spl1_mct";
1343 +    nrecPV_spl1_mct_ = 0;  
1344 +  }
1345 +  if (datatype == 2) {
1346 +    suffix = "_spl2_mct";
1347 +    nrecPV_spl2_mct_ = 0;    
1348 +  }
1349 +  if (datatype == 3) {
1350 +    suffix = "_mct";
1351 +    nrecPV_mct_ = 0;  
1352 +  }
1353 +
1354 +  //========================================================
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){
1362 +    vsim->recVtx=NULL;  
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 +      //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);
1372 +    }
1373 +  }
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() {
1476 +  fntrk_ = 0;
1477 +  //pvtx position (x,y,z) residual and error
1478 +  for(int i = 0; i<3;i++)    {
1479 +    fres_[i] = 0;
1480 +    ferror_[i] = 0;
1481 +  }
1482    
1483 <  // 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;
1489 >  nrecPV1_spl_ = 0;
1490    nrecPV2_spl_ = 0;
1491 <  // number of tracks in the vertex
1492 <  fntrk_ = 0;
1493 <  fntrk_spl1_mct_ = 0;
1494 <  fntrk_spl2_mct_ = 0;
1495 <  fntrk_mct_ = 0;
1496 <  //pvtx position (x,y,z) residual and error
1497 <  for(int i = 0; i<3;i++)
1498 <    {
1499 <      fres_[i] = 0;
1500 <      ferror_[i] = 0;
1501 <      fres_spl1_mct_[i] = 0;
1502 <      ferror_spl1_mct_[i] = 0;
1503 <      fres_spl2_mct_[i] = 0;
1504 <      ferror_spl2_mct_[i] = 0;
1505 <      fres_mct_[i] = 0;
1506 <      ferror_mct_[i] = 0;
1491 >  nrecPV_twovtx_ = 0;
1492 >  nsimPV_ = 0;
1493 >
1494 >  nrecPV_mct_ = 0;
1495 >  nrecPV_spl1_mct_ = 0;
1496 >  nrecPV_spl2_mct_ = 0;
1497 >
1498 >  // Mininum separation between the secondary pvtxes and leading pvtx
1499 >  min_zsep_ = 9999.0;
1500 >  min_ntrksep_ = 9999.0;
1501 >  min_sumpt2sep_ = -9999.0;
1502 >  
1503 >  // Variables filled per Vertex
1504 >  for (int i = 0; i < nMaxPVs_; i++) {
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 >    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;
1525 >    errx_twovtx_[i] = 0;
1526 >    erry_twovtx_[i] = 0;
1527 >    errz_twovtx_[i] = 0;
1528 >    pullx_twovtx_[i] = 0;
1529 >    pully_twovtx_[i] = 0;
1530 >    pullz_twovtx_[i] = 0;
1531 >    
1532 >    //simpv
1533 >    nsimTrkPV_[i] = 0;
1534 >    simx_[i] = 0;
1535 >    simy_[i] = 0;
1536 >    simz_[i] = 0;
1537 >    simptsq_[i] = 0;
1538 >    
1539 >    // residual and pulls with mc truth required
1540 >    deltax_mct_[i] = 0;
1541 >    deltay_mct_[i] = 0;
1542 >    deltaz_mct_[i] = 0;
1543 >    pullx_mct_[i] = 0;
1544 >    pully_mct_[i] = 0;
1545 >    pullz_mct_[i] = 0;
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 >    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 >  }
1583 > }
1584 >
1585 > double PVStudy::sumPtSquared(const reco::Vertex & v)  {
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 +  }
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