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.19 by yygao, Sat Jan 30 18:27:09 2010 UTC

# 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 >  bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
90 >  
91    //now do what ever initialization is needed
92    pvinfo.clear();
93  
87 –
94    // Specify the data mode vector
95    if(realData_) datamode.push_back(0);
96    else {
# Line 93 | Line 99 | PVStudy::PVStudy(const edm::ParameterSet
99      datamode.push_back(2);
100      datamode.push_back(3);
101    }
102 <
97 <  // Create ntuple files if needed
102 > // Create ntuple files if needed
103    if(saventuple_) {
104      file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
105      ftree_ = new TTree("mytree","mytree");
# Line 102 | Line 107 | PVStudy::PVStudy(const edm::ParameterSet
107      ftree_->Branch("residual",&fres_,"fres_[3]/D");
108      ftree_->Branch("error",&ferror_,"ferror_[3]/D");
109      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 –  }
110  
111 <  edm::Service<TFileService> fs;
112 <  TFileDirectory subDir = fs->mkdir( "Summary" );
125 <  TFileDirectory subDir1 = fs->mkdir( "Others" );
126 <  //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
111 >    // pvtxtree_ analyzing the pvtxs ootb
112 >    pvtxtree_ = new TTree("pvtxtree","pvtxtree");
113  
114 <  setRootStyle();
114 >    // Event information for the data
115 >    pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
116 >    pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
117 >    pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
118 >    pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");  
119 >    
120 >
121 >    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
122 >    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
123 >    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
124 >    // Event level information
125 >    pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126 >    pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127 >    pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128 >
129 >    //Fill the variables in the twovtx pair (recvtx1, recvtx2)  
130 >    // Information related to the analyzing the two-vertex method
131 >
132 >
133 >
134 >    pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
135 >    pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
136 >    pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
137 >    pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
138 >    pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
139 >    pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
140 >    pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
141 >    pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
142 >    pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
143 >    pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
144 >    pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
145 >
146 >    // Information for the splitVertexColl1
147 >    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
148 >    pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
149 >    pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
150 >    pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
151 >    pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
152 >    pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
153 >    pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
154    
155 <  //Book histograms for track and pvtx parameters for each splitted pvtx
156 <  //i=0 : x (as in unsplit track collection)
157 <  //i=1 : x1_spl_
158 <  //i=2 : x2_spl_
159 <  
160 <  for(int i=0;i<3;i++)
161 <    {  
162 <      string suffix;
163 <      if(i == 0) suffix = "";
164 <      else {
165 <        stringstream ss;
166 <        ss << i;
167 <        suffix =ss.str()+"_spl";  
168 <      }
169 <      h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
170 <      h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
171 <      h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
172 <      h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
173 <      h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);  
174 <      h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);  
175 <      
176 <      h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
177 <      h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
178 <      h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
179 <      h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
180 <      h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);  
181 <      h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
182 <      h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);    
183 <      suffix.clear();
155 >    
156 >    // Information for the splitVertexColl2
157 >    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
158 >    pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
159 >    pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
160 >    pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
161 >    pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
162 >    pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
163 >    pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
164 >  
165 >    
166 >    // MC variables
167 >    if(!realData_) {  
168 >      pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
169 >      pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
170 >      pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
171 >      pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
172 >      pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
173 >      pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
174 >
175 >      // For pvtxs with all the tracks
176 >      pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
177 >      pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
178 >      pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
179 >      pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
180 >      pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
181 >      pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
182 >      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");    
183 >      pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
184 >      pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
185 >      pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
186 >      pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
187 >      pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
188 >      pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
189 >      pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
190 >
191 >      // For pvtxs with splittracks1
192 >      pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
193 >      pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
194 >      pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
195 >      pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
196 >      pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
197 >      pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
198 >      pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
199 >      pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
200 >      pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
201 >      pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
202 >      pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
203 >      pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
204 >      pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
205 >      pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
206 >        
207 >      // For pvtxs with splittracks1
208 >      pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
209 >      pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
210 >      pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
211 >      pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
212 >      pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
213 >      pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
214 >      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");  
215 >      pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
216 >      pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
217 >      pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");  
218 >      pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
219 >      pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
220 >      pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
221 >      pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
222 >
223      }
224 <  h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
225 <  h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
226 <  h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
224 >  }
225 >  
226 >  setRootStyle();
227 >
228 >  //========================================
229 >  // Booking histograms
230 >  //========================================
231    
232 +  // Create a root file for histograms
233 +  theFile = new TFile(histoFileName_.c_str(), "RECREATE");
234 +  // make diretories
235 +  theFile->mkdir("Summary");
236 +  theFile->cd();
237 +  theFile->mkdir("Others");
238 +  theFile->cd();
239 +  if (analyze_) {
240 +    theFile->mkdir("Results");
241 +    theFile->cd();
242 +  }
243 +
244 +  // Book Histograms:
245 +  h_pvtrk   = new PVHistograms();
246 +  h_misc    = new PVHistograms();
247 +  h_summary = new PVHistograms();
248 +  h_others  = new PVHistograms();
249 +
250 +  for(int i=0;i<3;i++) {  
251 +    if(i == 0) h_pvtrk->Init("pvTrk");
252 +    else {
253 +      stringstream ss;
254 +      ss << i;
255 +      h_pvtrk->Init("pvTrk", ss.str(),"spl");
256 +    }
257 +  }
258 +  h_misc->Init("misc");
259 +
260 +  // Book MC only plots
261 +  if (!realData_) {
262 +    h_gen = new PVHistograms();
263 +    h_gen->Init("generator");
264 +  }
265 +  if (analyze_) h_ana = new PVHistograms();
266 +
267    //Book histograms sensitive to data/mc
268    for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
269      string suffix;
270 <    if(verbose_)  cout<<"datamode = "<< *it<<endl;
270 >    edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
271      switch(*it) {
272      case 0: suffix = "";
273        break;
# Line 175 | Line 278 | PVStudy::PVStudy(const edm::ParameterSet
278      case 3: suffix = "_mct";
279        break;
280      }
281 <    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);
281 >    h_summary->Init("summary",suffix);
282      
283      // Book residual, error and pull histograms for each nTrk bin
284      for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
285        stringstream ss;
286        ss << ntrk;
287 <      string suffix2 = suffix+"_"+ss.str();
288 <      h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
289 <      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 <    
287 >      h_others->Init("others",suffix,ss.str());
288 >    }
289 >
290      // Book residual and pull histograms only when analyzeOntheFly is enabled
291 <    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_) {
291 >    if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292      suffix.clear();
293 <  } // End of Book histograms sensitive to data/mc  
294 <  // Book MC only plots
295 <  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 <  }
293 >  } // End of Book histograms sensitive to data/mc    
294 >  
295 >
296   }
297  
298   PVStudy::~PVStudy()
# Line 255 | Line 300 | PVStudy::~PVStudy()
300  
301    // do anything here that needs to be done at desctruction time
302    // (e.g. close files, deallocate resources etc.)
303 +  theFile->cd();
304 +  theFile->cd("Summary");
305 +  h_pvtrk->Save();
306 +  h_misc->Save();
307 +  h_summary->Save();
308 +  if (!realData_)
309 +    h_gen->Save();
310 +  if (analyze_) {
311 +    theFile->cd();
312 +    theFile->cd("Results");
313 +    h_ana->Save();
314 +  }
315 +  theFile->cd();
316 +  theFile->cd("Others");
317 +  h_others->Save();
318  
319 +  theFile->Close();
320   }
321  
322   void PVStudy::setRootStyle() {
# Line 279 | Line 340 | void PVStudy::setRootStyle() {
340    gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
341    gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
342    gStyle->SetStatFontSize(0.08);      // size for stat. box
343 <  gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
344 <  gStyle->SetTitleFont(32, "");       // same for pad title
345 <  gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
346 <  gStyle->SetStatFont(32);            // same for stat. box
343 >  gStyle->SetTitleFont(42, "XYZ");    // times-bold-italic font (p. 153) for axes
344 >  gStyle->SetTitleFont(42, "");       // same for pad title
345 >  gStyle->SetLabelFont(42, "XYZ");    // same for axis labels
346 >  gStyle->SetStatFont(42);            // same for stat. box
347    gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
348    //
349    return;
# Line 292 | Line 353 | void PVStudy::setRootStyle() {
353   //
354   std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
355   {
356 <  std::vector<PVStudy::simPrimaryVertex> simpv;
356 > std::vector<PVStudy::simPrimaryVertex> simpv;
357    const HepMC::GenEvent* evt=evtMC->GetEvent();
358    if (evt) {
359 <    if(verbose_) {
360 <      cout << "process id " << evt->signal_process_id()<<endl;
361 <      cout << "signal process vertex " << ( evt->signal_process_vertex() ?
362 <                                            evt->signal_process_vertex()->barcode() : 0 )   <<endl;
302 <      cout << "number of vertices " << evt->vertices_size() << endl;
303 <    }
359 >    edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
360 >    edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
361 >                                                                        evt->signal_process_vertex()->barcode() : 0 ) <<endl;
362 >    edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
363  
364      int idx=0; int npv=0;
365      for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
# Line 311 | Line 370 | std::vector<PVStudy::simPrimaryVertex> P
370        // t component of PV:
371        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
372              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
373 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
373 >        //        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
374          HepMC::GenVertex * mv=(*mother)->production_vertex();
375          if( ((*mother)->status() == 3) && (!mv)) {
376 <          //        std::cout << "npv= " << npv << std::endl;
376 >          //        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
377            if (npv == 0) {
378 <            h["genPart_cT"]->Fill(pos.t()); // mm
379 <            h["genPart_T"]->Fill(pos.t()/299.792458); // ns
378 >            h_gen->Fill1d("genPart_cT", pos.t()); // mm
379 >            h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
380            }
381            npv++;
382          }
# Line 325 | Line 384 | std::vector<PVStudy::simPrimaryVertex> P
384        //       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
385          
386        bool hasMotherVertex=false;
387 <      if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
387 >      if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
388        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
389              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
390          HepMC::GenVertex * mv=(*mother)->production_vertex();
391 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
391 >        // if (verbose_) cout << "Status = " << (*mother)->status() << endl;
392          if (mv) {
393            hasMotherVertex=true;
394            if(!verbose_) break; //if verbose_, print all particles of gen vertices
# Line 356 | Line 415 | std::vector<PVStudy::simPrimaryVertex> P
415  
416        if(!vp){
417          // this is a new vertex
418 <        if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
418 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
419          simpv.push_back(sv);
420          vp=&simpv.back();
421        }else{
422 <        if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
422 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
423        }
424        vp->genVertex.push_back((*vitr)->barcode());
425        // collect final state descendants
# Line 384 | Line 443 | std::vector<PVStudy::simPrimaryVertex> P
443              }
444            }
445          }
446 <      }
447 <    }
446 >      }//loop MC vertices daughters
447 >    }//loop MC vertices
448    }
449    return simpv;
450   }
# Line 418 | Line 477 | std::vector<PVStudy::simPrimaryVertex> P
477                ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
478              }
479            }
480 <          primary =  ( gv->position().t()==0);
480 >          //primary =  ( gv->position().t()==0 );
481 >          primary = true;
482 >          h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
483 >          h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
484 >
485            //resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
486            // no more mother pointer in the improved HepMC GenParticle
487            resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
# Line 467 | Line 530 | PVStudy::analyze(const edm::Event& iEven
530    using namespace edm;
531    using namespace reco;
532    
533 <  if (verbose_)
534 <    cout<<"PVStudy::analyze():"<<endl;
535 <  //=======================================================
533 >  //========================================================================
534 >  // Step 0: Prepare root variables and get information from the Event
535 >  //========================================================================
536 >  
537 >  edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
538    // Initialize Root-tuple variables if needed
539 <  //=======================================================
540 <  //if(saventuple_)
541 <    SetVarToZero();
477 <  
478 <  //=======================================================
479 <  // Track accessors
480 <  //=======================================================
481 <  //trackCollection
539 >  SetVarToZero();
540 >  
541 >  // ====== TrackCollection
542    static const reco::TrackCollection s_empty_trackColl;
543    const reco::TrackCollection *trackColl = &s_empty_trackColl;
544    edm::Handle<reco::TrackCollection> trackCollectionHandle;
# Line 486 | Line 546 | PVStudy::analyze(const edm::Event& iEven
546    if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
547      trackColl = trackCollectionHandle.product();
548    } else {
549 <    cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
549 >    edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
550    }
551 <  //splitTrackCollection1
551 >  // ====== splitTrackCollection1
552    static const reco::TrackCollection s_empty_splitTrackColl1;
553    const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
554    edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
# Line 496 | Line 556 | PVStudy::analyze(const edm::Event& iEven
556    if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
557      splitTrackColl1 = splitTrackCollection1Handle.product();
558    } else {
559 <    cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
559 >    edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
560    }
561 <  //splitTrackCollection2
561 >  // ====== splitTrackCollection2
562    static const reco::TrackCollection s_empty_splitTrackColl2;
563    const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
564    edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
# Line 506 | Line 566 | PVStudy::analyze(const edm::Event& iEven
566    if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
567      splitTrackColl2 = splitTrackCollection2Handle.product();
568    } else {
569 <    cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
569 >   edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
570    }
571    
572 <  //=======================================================
513 <  // PVTX accessors
514 <  //=======================================================
515 <  //vertexCollection
572 >  // ======= PrimaryVertexCollection
573    static const reco::VertexCollection s_empty_vertexColl;
574    const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
575    edm::Handle<reco::VertexCollection> vertexCollectionHandle;
# Line 520 | Line 577 | PVStudy::analyze(const edm::Event& iEven
577    if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
578      vertexColl = vertexCollectionHandle.product();
579    } else {
580 <    cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
580 >   edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
581    }
582 <  //splitVertexCollection1
582 >  // ====== splitVertexCollection1
583    static const reco::VertexCollection s_empty_splitVertexColl1;
584    const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
585    edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
# Line 530 | Line 587 | PVStudy::analyze(const edm::Event& iEven
587    if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
588      splitVertexColl1 = splitVertexCollection1Handle.product();
589    } else {
590 <    cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
590 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
591    }
592 <  //splitVertexCollection2
592 >  // ====== splitVertexCollection2
593    static const reco::VertexCollection s_empty_splitVertexColl2;
594    const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
595    edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
# Line 540 | Line 597 | PVStudy::analyze(const edm::Event& iEven
597    if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
598      splitVertexColl2 = splitVertexCollection2Handle.product();
599    } else {
600 <    cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
600 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
601    }
602  
546 –  if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
603    
604 <  //=======================================================
605 <  // MC simvtx accessor
606 <  //=======================================================
604 >  // ======== BeamSpot accessors
605 >  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
606 >  iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
607 >  reco::BeamSpot bs = *recoBeamSpotHandle;    
608 >  const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
609 >  
610 >  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
611 >  
612 >  // ========== MC simvtx accessor
613    if (!realData_) {
614      edm::Handle<SimVertexContainer> simVtxs;
615      iEvent.getByLabel( simG4_, simVtxs);
# Line 556 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618      iEvent.getByLabel( simG4_, simTrks);
619    }
620  
621 <  //=======================================================
560 <  // GET PDT
561 <  //=======================================================
621 >  // ========== GET PDT
622    try{
623      iSetup.getData(pdt);
624    }catch(const Exception&){
625 <    cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
625 >    edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626 >  }
627 >  
628 >  // ======= Get MC information if needed
629 >  bool MC=false;  
630 >  Handle<HepMCProduct> evtMC;
631 >  if (!realData_) {
632 >    iEvent.getByLabel("generator",evtMC);
633 >    if (!evtMC.isValid()) {
634 >      MC=false;
635 >      edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
636 >    } else {
637 >      edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
638 >      MC=true;
639 >    }
640    }
641  
642 <  //=======================================================
643 <  // Fill trackparameters of the input tracks to pvtx fitter
644 <  //=======================================================
645 <
646 <  // GeneralTracks  
647 <  h["nTrk"]->Fill(trackColl->size());
648 <  for(TrackCollection::const_iterator itTrack = trackColl->begin();
649 <      itTrack != trackColl->end();                      
650 <      ++itTrack) {
651 <    h["trkPt"]->Fill(itTrack->pt());
652 <    h["trkDxy"]->Fill(itTrack->dxy());
653 <    h["trkDz"]->Fill(itTrack->dz());
580 <    h["trkEta"]->Fill(itTrack->eta());
581 <    h["trkPhi"]->Fill(itTrack->phi());
642 >  //========================================================================
643 >  // Step 1:  Apply event cleaning for data and MC
644 >  //          WARNING: event selection cut are hard coded!!
645 >  //========================================================================
646 >  // =====Loop over the trackColl to get the fraction of HighPurity tracks
647 >  int nTracks = 0;
648 >  int nHighPurityTracks=0;
649 >  double fHighPurity=0;
650 >
651 >  for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
652 >    TrackRef tkref(trackCollectionHandle,i);
653 >    if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
654    }
655 <  
656 <  //SplittedTracks1
657 <  h["nTrk1_spl"]->Fill(splitTrackColl1->size());
658 <  for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
659 <      itTrack != splitTrackColl1->end();                      
660 <      ++itTrack) {
661 <    h["trkPt1_spl"]->Fill(itTrack->pt());
662 <    h["trkDxy1_spl"]->Fill(itTrack->dxy());
663 <    h["trkDz1_spl"]->Fill(itTrack->dz());
664 <    h["trkEta1_spl"]->Fill(itTrack->eta());
665 <    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());
655 >  if(nTracks>0)
656 >    fHighPurity =  double(nHighPurityTracks)/double(nTracks);
657 >  if(verbose_)
658 >    cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
659 >    
660 >  if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
661 >    glob_runno_ = iEvent.id().run();
662 >    glob_evtno_ = iEvent.id().event();
663 >    glob_ls_   = iEvent.luminosityBlock();
664 >    glob_bx_  = iEvent.bunchCrossing();  
665 >    return;
666    }
667 +
668 +  //========================================================================
669 +  // Step 2: Fill histograms for the splitting consistency checks
670 +  //========================================================================
671    
672 <  if(verbose_)
673 <    cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
674 <
675 <  //=======================================================
676 <  // Fill reconstructed vertices
677 <  //=======================================================
678 <  if(verbose_)  {
679 <    cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
672 >  // === Fill trackparameters of the input tracks to pvtx fitter
673 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
674 >  //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
675 >  // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
676 >  fillTrackHisto(trackColl, 0, beamSpot);
677 >  fillTrackHisto(splitTrackColl1, 1, beamSpot);
678 >  fillTrackHisto(splitTrackColl2, 2, beamSpot);
679 >  edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
680 >  
681 >  
682 >  // ==== Fill number of reconstructed vertices
683 >  edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
684 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
685 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
686 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
687 >  if (verbose_) {
688      printRecVtxs(vertexCollectionHandle);
616 –    cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
689      printRecVtxs(splitVertexCollection1Handle);
618 –    cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
690      printRecVtxs(splitVertexCollection2Handle);
691    }
692 +
693 +  nrecPV_ = int (vertexColl->size());
694 +  nrecPV1_spl_ = int (splitVertexColl1->size());
695 +  nrecPV2_spl_ = int (splitVertexColl2->size());
696 +
697 +  h_pvtrk->Fill1d("nrecPV", nrecPV_);
698 +  h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
699 +  h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
700 +  h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
701    
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_));
702    
703 <  // Fill the track paramters for the vertexColl
704 <  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
705 <      v!=vertexColl->end(); ++v){
706 <    h["nTrkPV"]->Fill(v->tracksSize());
707 <    try {
708 <      for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
709 <          t!=v->tracks_end(); t++) {
710 <        // illegal charge
711 <        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
712 <          //      h["trkPtPV"]->Fill(0.);
713 <        }
714 <        else {
715 <          h["trkPtPV"]->Fill((**t).pt());
716 <          h["trkDxyPV"]->Fill((**t).dxy());
717 <          h["trkDzPV"]->Fill((**t).dz());
718 <          h["trkEtaPV"]->Fill((**t).eta());
719 <          h["trkPhiPV"]->Fill((**t).phi());
720 <        }
703 >  // ======= Fill track parameter ntuple/hist for tracks used in recoVertices
704 >  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
705 >  if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
706 >    fillTrackHistoInPV(vertexColl, 0, beamSpot);
707 >  
708 >  // ======= Fill secondary/primary min separations for multi-vertices
709 >  if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {  
710 >    double min_xsep = 9999.0;
711 >    double min_ysep = 9999.0;
712 >    double min_zsep = 9999.0;    
713 >    double min_xsep_sign = 9999.0;
714 >    double min_ysep_sign = 9999.0;
715 >    double min_zsep_sign = 9999.0;
716 >    double min_ntrksep = 9999.0;
717 >    double min_sumpt2sep = 9999999.0;
718 >    
719 >    edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
720 >
721 >    // Looping through the secondary vertices to find the mininum separation to the primary
722 >    for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
723 >        v!=vertexColl->end(); ++v) {  
724 >      if(v->isValid() && ! v->isFake()) {
725 >        // xsep
726 >        if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
727 >          min_xsep = fabs(v->x()- vertexColl->begin()->x());
728 >        // ysep
729 >        if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
730 >          min_ysep = fabs(v->y()- vertexColl->begin()->y());
731 >        // zsep
732 >        if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
733 >          min_zsep = fabs(v->z()- vertexColl->begin()->z());
734 >        // xsep_sign
735 >        double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
736 >        if(xsep_sign < min_xsep_sign)
737 >          min_xsep_sign = xsep_sign;
738 >        // ysep_sign
739 >        double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
740 >        if(ysep_sign < min_ysep_sign)
741 >          min_ysep_sign = ysep_sign;
742 >        // zsep_sign
743 >        double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
744 >        if(zsep_sign < min_zsep_sign)
745 >          min_zsep_sign = zsep_sign;
746 >        // ntrksep
747 >        if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
748 >          min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
749 >        // sumpt2sep
750 >        if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
751 >          min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
752        }
753      }
754 <    catch (...) {
755 <      // exception thrown when trying to use linked track
756 <      //          h["trkPtPV"]->Fill(0.);
757 <    }
758 <  }
754 >    h_misc->Fill1d("min_xsep", min_xsep);
755 >    h_misc->Fill1d("min_ysep", min_ysep);
756 >    h_misc->Fill1d("min_zsep", min_zsep);
757 >    h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
758 >    h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
759 >    h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
760 >    h_misc->Fill1d("min_ntrksep", min_ntrksep);
761 >    h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
762 >
763 >    min_zsep_ = min_zsep;
764 >    min_ntrksep_ = min_ntrksep;
765 >    min_sumpt2sep_ = min_sumpt2sep;
766 >  } // End of  if(vertexColl->size()>1) {
767    
768 +  edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
769 +
770 +  //========================================================================
771 +  // Step 3:  PrimaryVertex Matching
772 +  //   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
773 +  //   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
774 +  //   3. The first match is the primary vertex, which has largest sum(pT^2)
775 +  //========================================================================
776    
777 <  // == ignore multi-vertices (for now)
778 <  //take single pvt collection and compare
779 <  if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
780 <    const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
781 <    const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
782 <    if (matchVertex(recvtx1, recvtx2))  {
783 <      int nTrkPV1 = recvtx1.tracksSize();
784 <      int nTrkPV2 = recvtx2.tracksSize();    
785 <      h["nTrkPV1_spl"]->Fill(nTrkPV1);
786 <      h["nTrkPV2_spl"]->Fill(nTrkPV2);
787 <      double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
788 <      h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
789 <      h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
790 <      
791 <      // Fill the track paramters for the splitVertexColl1
792 <      try {
793 <        for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
794 <            t!=recvtx1.tracks_end(); t++) {
795 <          // illegal charge
796 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
797 <            //    h["trkPtPV"]->Fill(0.);
798 <          }
681 <          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 <          }
688 <        }
689 <      }
690 <      catch (...) {
691 <        // exception thrown when trying to use linked track
692 <        //        h["trkPtPV"]->Fill(0.);
693 <      }
777 >  edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
778 >  reco::VertexCollection s_empty_matchedVertexColl1;
779 >  reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
780 >  matchedVertexColl1->reserve(splitVertexColl1->size());
781 >  reco::VertexCollection s_empty_matchedVertexColl2;
782 >  reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
783 >  matchedVertexColl2->reserve(splitVertexColl2->size());
784 >  
785 >  bool stopmatching = false;
786 >  for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
787 >    reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
788 >    if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
789 >    for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {    
790 >      //------------------------------------------------------------------------
791 >      //== If only considering matching the first vertex of the splitVertexColl
792 >      //------------------------------------------------------------------------
793 >      if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
794 >      reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
795 >      if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
796 >      edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
797 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
798 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
799        
800 <      // Fill the track paramters for the splitVertexColl2
801 <      try {
802 <        for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
803 <            t!=recvtx2.tracks_end(); t++) {
804 <          // illegal charge
805 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
806 <            //    h["trkPtPV"]->Fill(0.);
807 <          }
808 <          else {
809 <            h["trkPtPV2_spl"]->Fill((**t).pt());
810 <            h["trkDxyPV2_spl"]->Fill((**t).dxy());
811 <            h["trkDzPV2_spl"]->Fill((**t).dz());
812 <            h["trkEtaPV2_spl"]->Fill((**t).eta());
813 <            h["trkPhiPV2_spl"]->Fill((**t).phi());
814 <          }
800 >      double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
801 >      h_misc->Fill1d("twovtxzsign", twovtxsig);
802 >      if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
803 >        edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
804 >
805 >        int nTrkPV1 = recvtx1->tracksSize();
806 >        int nTrkPV2 = recvtx2->tracksSize();    
807 >        double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
808 >        h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
809 >        h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
810 >        if(fabs(ntrkreldiff)<ntrkdiffcut_) {    
811 >          edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
812 >
813 >          matchedVertexColl1->push_back(*recvtx1);
814 >          matchedVertexColl2->push_back(*recvtx2);
815 >          stopmatching = true; // stop the matching when the first match is found!
816 >          break;
817          }
818 +        else
819 +          edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
820        }
821 <      catch (...) {
822 <        // exception thrown when trying to use linked track
823 <        //        h["trkPtPV"]->Fill(0.);
821 >      else {
822 >        edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
823 >        edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
824 >        edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;  
825        }
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;
826      }
827 <  } // End of fill histograms for two-vertices Method
827 >  }
828 >  edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
829  
830 <  if(verbose_)
831 <    cout<<"Done filling res/pull with two-vertices method"<<endl;
832 <
833 <  //=======================================================
834 <  // Fill simulated vertices
835 <  //=======================================================
836 <  if (!realData_) {
837 <    bool MC=false;
838 <    Handle<HepMCProduct> evtMC;
839 <    iEvent.getByLabel("generator",evtMC);
840 <    if (!evtMC.isValid()) {
841 <      MC=false;
842 <      if(verbose_){
843 <        cout << "no HepMCProduct found"<< endl;
844 <      }
845 <    } else {
846 <      if(verbose_){
847 <        cout << "generator HepMCProduct found"<< endl;
848 <      }
849 <      MC=true;
830 >  //========================================================================
831 >  // Step 4: Analyze matchedVertexColls
832 >  //========================================================================
833 >  edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
834 >  
835 >  // ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
836 >  if(MC){
837 >    // make a list of primary vertices:
838 >    std::vector<simPrimaryVertex> simpv;
839 >    simpv=getSimPVs(evtMC,"");
840 >    //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
841 >    h_gen->Fill1d("nsimPV", simpv.size());
842 >    nsimPV_ = simpv.size();
843 >    int isimpv = 0;
844 >    for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
845 >        vsim!=simpv.end(); vsim++, isimpv++){
846 >      nsimTrkPV_[isimpv]  =vsim->nGenTrk;
847 >      simx_[isimpv] = vsim->x;
848 >      simy_[isimpv] = vsim->y;
849 >      simz_[isimpv] = vsim->y;
850 >      simptsq_[isimpv] = vsim->ptsq;
851 >      //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
852 >      fillMCHisto(vsim, isimpv, vertexColl, 3);
853      }
854 <    
854 >  }
855 >  
856 >  // Compare the reconstructed vertex position and calculate resolution/pulls
857 >  if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {  
858 >    // ==== Analyze the splitted vtx using MC method
859      if(MC){
860        // make a list of primary vertices:
861        std::vector<simPrimaryVertex> simpv;
862        simpv=getSimPVs(evtMC,"");
863        //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
864 <      h["nsimPV"]->Fill(simpv.size());
783 <      int nsimtrk=0;
784 <      int nrectrk=0;
785 <      
864 >      int isimpv = 0;
865        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
866 <          vsim!=simpv.end(); vsim++){
867 <        nsimtrk+=vsim->nGenTrk;
868 <        // Loop over the datamodes, for MC, fill only 1, 2, 3
869 <        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
866 >          vsim!=simpv.end(); vsim++, isimpv++){
867 >        //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
868 >        fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
869 >        fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
870        }
871 <    } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
872 <  } // End of Filling MC variables
871 >    }
872 >
873 >    // ==== Analyze res/pull two-vertex method
874 >    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
875 >    fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
876 >    fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
877 >
878 >    reco::VertexCollection::const_iterator v1;
879 >    reco::VertexCollection::const_iterator v2;
880 >
881 >    nrecPV_twovtx_ = 0;
882 >    for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
883 >        v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
884 >        ++v1, ++v2) {
885 >      h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
886 >      h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
887 >      fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
888 >      fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
889 >      fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
890 >      ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
891 >      ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
892 >      ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
893 >      fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
894 >      
895 >      h_summary->Fill1d("deltax", fres_[0] );
896 >      h_summary->Fill1d("deltay", fres_[1] );
897 >      h_summary->Fill1d("deltaz", fres_[2] );
898 >      h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
899 >      h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
900 >      h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
901 >      h_summary->Fill1d("errPVx", ferror_[0] );
902 >      h_summary->Fill1d("errPVy", ferror_[1] );
903 >      h_summary->Fill1d("errPVz", ferror_[2] );
904 >      pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
905 >                                       error(ferror_[0], ferror_[1], ferror_[2]),
906 >                                       fntrk_));
907 >      // Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
908 >      fillHisto(res(fres_[0], fres_[1], fres_[2]),
909 >                error(ferror_[0], ferror_[1], ferror_[2]),
910 >                fntrk_,0);  
911 >
912 >      // Fill the ntuple variables
913 >      nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
914 >      deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
915 >      deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
916 >      deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
917 >      errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
918 >      erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
919 >      errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
920 >      pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
921 >      pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
922 >      pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];    
923 >
924 >
925 >      //SplittedVertex
926 >      
927 >      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
928 >      ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
929 >      normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
930 >      avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
931 >      errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
932 >      erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
933 >      errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
934 >  
935 >      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
936 >      ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
937 >      normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
938 >      avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
939 >      errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
940 >      erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
941 >      errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
942 >  
943 >      nrecPV_twovtx_++;
944 >
945 >
946 >      
947 >      // Print some information of the two tracks events
948 >      if(verbose_ && fntrk_ < 4) {
949 >        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
950 >        printRecVtx(*v1);
951 >        cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
952 >        printRecVtx(*v2);
953 >      }
954 >    } // End of analyzing res/pull
955 >  } // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
956 >  else
957 >    edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
958 >  
959 >  edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
960 >
961 >
962 >  //========================================================================
963 >  // Step 5: Fill ntuple if saventuple_ is on
964 >  //========================================================================
965 >  if(saventuple_) {
966 >    ftree_->Fill();
967 >    pvtxtree_->Fill();
968 >  }
969 >  
970   }
971  
972 +
973 +  
974   // ------------ method called once each job just before starting event loop  ------------
975   void
976   PVStudy::beginJob()
# Line 894 | Line 980 | PVStudy::beginJob()
980   // ------------ method called once each job just after ending the event loop  ------------
981   void
982   PVStudy::endJob() {
983 <  if (verbose_) cout << "Analyzing PV info" << endl;
983 >  edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
984    // Looping through the datamodes available
985    for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {  
986      string suffix;
987 <    if(verbose_)  cout<<"datamode = "<< *it<<endl;
987 >    edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
988      switch(*it) {
989      case 0: suffix = "";
990        break;
# Line 909 | Line 995 | PVStudy::endJob() {
995      case 3: suffix = "_mct";
996        break;
997      }
998 <    if (analyze_) {
998 >    suffix += "_nTrk";
999 >    if(analyze_) {
1000        for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1001 +        edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1002          PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1003 <        if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1004 <        if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1005 <        if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1006 <        if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1007 <        if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1008 <        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.);
1003 >        if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1004 >        if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1005 >        if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1006 >        if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1007 >        if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1008 >        if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1009        }
1010      }
1011      suffix.clear();
# Line 926 | Line 1013 | PVStudy::endJob() {
1013    if(saventuple_) {
1014      file_->cd();
1015      ftree_->Write();
1016 +    pvtxtree_->Write();
1017      file_->Close();
1018    }
1019 < }
932 <
1019 >  
1020  
1021 + }
1022  
1023   // Match a recovertex with a simvertex
1024   // ? Should the cut be parameter dependent?
# Line 938 | Line 1026 | PVStudy::endJob() {
1026   bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1027                            const reco::Vertex & vrec)
1028   {
1029 <  return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1029 >  if(vrec.isFake() || !vrec.isValid()) return false;
1030 >  else
1031 >    return true;
1032 >    //return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1033   }
1034  
1035   // Match two reconstructed vertices
1036 < bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1037 <                           const reco::Vertex & vrec2)
1036 > bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1037 >                           const reco::VertexRef & vrec2,
1038 >                           int zsigncut)
1039   {
1040 <  return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1040 >  double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1041 >  edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1042 >  return (fabs(vrec1->z()-vrec2->z())<cut);
1043   }
1044  
1045  
1046   bool PVStudy::isResonance(const HepMC::GenParticle * p){
1047    double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1048 <  if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1048 >  edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1049    return  ctau >0 && ctau <1e-6;
1050   }
1051  
# Line 1002 | Line 1096 | void PVStudy::printRecVtxs(const Handle<
1096    }
1097   }
1098  
1099 + void PVStudy::printRecVtx(const reco::Vertex & v){
1100 +
1101 +  cout << "#trk " << std::setw(3) << v.tracksSize()
1102 +         << " chi2 " << std::setw(4) << v.chi2()
1103 +         << " ndof " << std::setw(3) << v.ndof() << endl
1104 +         << " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1105 +         << " dx " << std::setw(8) << v.xError()<< endl
1106 +         << " y "  << std::setw(8) << v.y()
1107 +         << " dy " << std::setw(8) << v.yError()<< endl
1108 +         << " z "  << std::setw(8) << v.z()
1109 +         << " dz " << std::setw(8) << v.zError()
1110 +         << endl;
1111 + }
1112 +
1113   void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1114    cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1115    int i=1;
# Line 1023 | Line 1131 | void PVStudy::printSimTrks(const Handle<
1131   void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1132    stringstream ss;
1133    ss << ntk;
1134 <  string suffix = ss.str();
1135 <  if(datamode == 0 )  suffix = "_" + suffix;
1136 <  if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1137 <  if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1138 <  if(datamode == 3 )  suffix = "_mct_" + suffix;
1134 >  string suffix;
1135 >  if(datamode == 0 )  suffix = "_" + ss.str();
1136 >  if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1137 >  if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1138 >  if(datamode == 3 )  suffix = "_mct_" + ss.str();
1139    
1140    if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1141    // Fill histograms of res and pull of ntk
1142 <  h["deltax"+suffix]->Fill(r.x());
1143 <  h["deltay"+suffix]->Fill(r.y());
1144 <  h["deltaz"+suffix]->Fill(r.z());
1145 <  h["pullx"+suffix]->Fill(r.x()/er.x());
1146 <  h["pully"+suffix]->Fill(r.y()/er.y());
1147 <  h["pullz"+suffix]->Fill(r.z()/er.z());  
1148 <  h["errPVx"+suffix]->Fill(er.x());
1149 <  h["errPVy"+suffix]->Fill(er.y());
1150 <  h["errPVz"+suffix]->Fill(er.z());
1142 >  h_others->Fill1d("deltax"+suffix, r.x());
1143 >  h_others->Fill1d("deltay"+suffix, r.y());
1144 >  h_others->Fill1d("deltaz"+suffix, r.z());
1145 >  h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1146 >  h_others->Fill1d("pully"+suffix, r.y()/er.y());
1147 >  h_others->Fill1d("pullz"+suffix, r.z()/er.z());  
1148 >  h_others->Fill1d("errPVx"+suffix, er.x());
1149 >  h_others->Fill1d("errPVy"+suffix, er.y());
1150 >  h_others->Fill1d("errPVz"+suffix, er.z());
1151   }
1152  
1153 +
1154   PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1155    map<int,double> data;
1156    data.clear();
# Line 1059 | Line 1168 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1168      for ( int i = 0; i < 6; ++i) {
1169        switch (i) {
1170        case 0:
1171 <        fit(h["deltax"+suffix], fpar);
1171 >        if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1172 >        fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1173          data[i] = fpar[0]*cm2um;
1174          break;
1175        case 1:
1176 <        fit(h["deltay"+suffix], fpar);
1176 >        if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1177 >        fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1178          data[i] = fpar[0]*cm2um;
1179          break;
1180        case 2:
1181 <        fit(h["deltaz"+suffix], fpar);
1181 >        if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1182 >        fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1183          data[i] = fpar[0]*cm2um;
1184          break;
1185        case 3:
1186 <        fit(h["pullx"+suffix], fpar);
1186 >        if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1187 >        fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1188          data[i] = fpar[0];
1189          break;
1190        case 4:
1191 <        fit(h["pully"+suffix], fpar);
1191 >        if (!h_others->ReadHisto1D("pully"+suffix)) break;
1192 >        fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1193          data[i] = fpar[0];
1194          break;
1195        case 5:
1196 <        fit(h["pullz"+suffix], fpar);
1196 >        if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1197 >        fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1198          data[i] = fpar[0];
1199          break;
1200        }
1201 <      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1201 >      if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1202      }
1203    }
1204    else
# Line 1098 | Line 1213 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1213                              ntk);
1214   }
1215  
1216 < void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1216 >
1217 > void PVStudy::fit(TH1 *hdist, double fitPars[]){
1218    TAxis *axis0 = hdist->GetXaxis();
1219    int nbin = axis0->GetLast();
1220    double nOF = hdist->GetBinContent(nbin+1);
1221    double nUF = hdist->GetBinContent(0);
1222    //   double fitRange = axis0->GetBinUpEdge(nbin);
1223 <  double fitRange = axis0->GetXmax();
1223 >  // double fitRange = axis0->GetXmax();
1224 >  double fitRange = 2*hdist->GetRMS();
1225    double sigMax[2] = {0.,0.};
1226  
1227 <  if ( verbose_ ){
1228 <    cout << "Last bin = " << nbin;
1229 <    cout << "; hdist: Overflow Entries = " << nOF;
1230 <    cout << "; Underflow Entries = " << nUF;
1231 <    cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1232 <    cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1116 <  }
1227 >  edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1228 >                           << "; hdist: Overflow Entries = " << nOF
1229 >                           << "; Underflow Entries = " << nUF
1230 >                           << "; hdist->GetEntries() = " << hdist->GetEntries()
1231 >                           << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1232 >
1233    if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1234      for (int bi = 0; bi < nbin; bi++) {
1235        if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
# Line 1125 | Line 1241 | void PVStudy::fit(TH1 *&hdist, double fi
1241          if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1242        }
1243      }
1244 <    if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1244 >    //edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1245  
1246      TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1247      fgaus->SetParameter(1, 0.);
1248 <    fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1248 >    fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1249      fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1250      fgaus->SetLineColor(4);
1251 <    hdist->Fit(fgaus,"Q");
1251 >    hdist->Fit(fgaus,"QLRM");
1252      gPad->Update();
1253      TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1254      s->SetOptStat(1111111);
1255      s->SetOptFit(0111);
1256      gPad->Update();
1257 <    if (verbose_) cout << "got function" << endl;
1257 >    //edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1258      fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1259    }
1260    else
1261      fitPars[0] = -999;
1262   }
1263  
1264 +
1265 + void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1266 +  string suffix;
1267 +  suffix = ""; // for unsplit tracks
1268 +  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1269 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1270 +  
1271 +  h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1272 +  for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1273 +      itTrack != trackColl->end();                      
1274 +      ++itTrack) {
1275 +    h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1276 +    h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1277 +    h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1278 +    h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1279 +    h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1280 +    h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1281 +  }
1282 + }
1283 +
1284 + void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1285 +  string suffix;
1286 +  suffix = ""; // for unsplit tracks
1287 +  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1288 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1289 +  int ivtx = 0;
1290 +  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1291 +      v!=vertexColl->end(); ++v, ++ivtx) {  
1292 +    if(!v->isFake()) {
1293 +      h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1294 +      h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1295 +      int nHWTrkPV_ = 0;
1296 +      try {
1297 +        for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1298 +            t!=vertexColl->begin()->tracks_end(); t++) {
1299 +          // illegal charge
1300 +          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1301 +            //    h_pvtrk->Fill1d("trkPtPV", 0.);
1302 +          }
1303 +          else {
1304 +            h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1305 +            h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());  
1306 +            h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1307 +            h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1308 +            h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1309 +            h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1310 +            if(vertexColl->begin()->trackWeight(*t) > 0.5)
1311 +              nHWTrkPV_++;
1312 +          }
1313 +        }
1314 +      }
1315 +      catch (...) {
1316 +        // exception thrown when trying to use linked track
1317 +        //        h_pvtrk->Fill1d("trkPtPV", 0.);
1318 +      }
1319 +      h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1320 +    }
1321 +  }
1322 + }
1323 +
1324 + void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1325 + {
1326 +  std::string suffix;
1327 +  // Set the vertexColl and histogram suffix according to datamode
1328 +  if (datatype == 1) {
1329 +    suffix = "_spl1_mct";
1330 +    nrecPV_spl1_mct_ = 0;  
1331 +  }
1332 +  if (datatype == 2) {
1333 +    suffix = "_spl2_mct";
1334 +    nrecPV_spl2_mct_ = 0;    
1335 +  }
1336 +  if (datatype == 3) {
1337 +    suffix = "_mct";
1338 +    nrecPV_mct_ = 0;  
1339 +  }
1340 +  
1341 +  //========================================================
1342 +  //  look for a matching reconstructed vertex in vertexColl
1343 +  //========================================================        
1344 +
1345 +  for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1346 +      vrec!=vtxColl->end(); ++vrec){
1347 +    vsim->recVtx=NULL;  
1348 +    edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1349 +    edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1350 +
1351 +    if ( matchVertex(*vsim,*vrec)  ) {
1352 +      vsim->recVtx=&(*vrec);
1353 +      
1354 +      // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1355 +      //if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1356 +      //|| (!vsim->recVtx) )
1357 +      //vsim->recVtx=&(*vrec);
1358 +      //}
1359 +      edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1360 +      
1361 +      double fres_mct[3];
1362 +      double ferror_mct[3];
1363 +      int nrectrk =  int(vsim->recVtx->tracksSize());
1364 +      double ndofPV = vsim->recVtx->ndof();
1365 +      double normchi2PV =  vsim->recVtx->normalizedChi2();
1366 +      double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1367 +
1368 +      fres_mct[0] = vsim->recVtx->x()-vsim->x;
1369 +      fres_mct[1] = vsim->recVtx->y()-vsim->y;
1370 +      fres_mct[2] = vsim->recVtx->z()-vsim->z;
1371 +      ferror_mct[0] = vsim->recVtx->xError();
1372 +      ferror_mct[1] = vsim->recVtx->yError();
1373 +      ferror_mct[2] = vsim->recVtx->zError();
1374 +      
1375 +      h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1376 +      h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1377 +      h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1378 +      h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1379 +      h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1380 +      h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1381 +      h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1382 +      h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1383 +      h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1384 +      pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1385 +                                       error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1386 +                                       nrectrk));
1387 +      // Fill histo according to its track multiplicity
1388 +      fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1389 +                error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1390 +                nrectrk, datatype);
1391 +      
1392 +      if(saventuple_) {
1393 +        //Fill the values for variables in pvtxtree_  
1394 +        if(datatype == 1) {
1395 +          nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1396 +          ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1397 +          normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1398 +          avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1399 +          deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1400 +          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1401 +          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1402 +          pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];  
1403 +          pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1404 +          pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1405 +          errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1406 +          erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1407 +          errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1408 +          nrecPV_spl1_mct_++;
1409 +        }
1410 +        if(datatype == 2) {
1411 +          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;  
1412 +          ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1413 +          normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1414 +          avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1415 +          deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1416 +          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1417 +          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1418 +          pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];  
1419 +          pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1420 +          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];        
1421 +          errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1422 +          erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1423 +          errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1424 +          nrecPV_spl2_mct_++;
1425 +        }
1426 +        if(datatype == 3) {    
1427 +          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1428 +          ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1429 +          normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1430 +          avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1431 +          deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1432 +          deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1433 +          deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1434 +          pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];  
1435 +          pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1436 +          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1437 +          errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1438 +          erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1439 +          errz_mct_[nrecPV_mct_] =  ferror_mct[2];        
1440 +          nrecPV_mct_++;
1441 +        }
1442 +      } // End of  if(saventuple_) {
1443 +    } //  if ( matchVertex(*vsim,*vrec) ) {
1444 +    else {  // no rec vertex found for this simvertex
1445 +      edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1446 +    }
1447 +  }
1448 + }
1449 +
1450   void PVStudy::SetVarToZero() {
1451 +  fntrk_ = 0;
1452 +  //pvtx position (x,y,z) residual and error
1453 +  for(int i = 0; i<3;i++)    {
1454 +    fres_[i] = 0;
1455 +    ferror_[i] = 0;
1456 +  }
1457    
1458 <  // number of reconstructed vertices
1458 >  // Event level information
1459 >  glob_runno_ = 0;
1460 >  glob_evtno_ = 0;
1461 >  glob_ls_ = 0;
1462 >  glob_bx_ = 0;
1463    nrecPV_ = 0;
1464 <  nrecPV1_spl_ = 0;
1464 >  nrecPV1_spl_ = 0;
1465    nrecPV2_spl_ = 0;
1466 <  // number of tracks in the vertex
1467 <  fntrk_ = 0;
1468 <  fntrk_spl1_mct_ = 0;
1469 <  fntrk_spl2_mct_ = 0;
1470 <  fntrk_mct_ = 0;
1471 <  //pvtx position (x,y,z) residual and error
1472 <  for(int i = 0; i<3;i++)
1473 <    {
1474 <      fres_[i] = 0;
1475 <      ferror_[i] = 0;
1476 <      fres_spl1_mct_[i] = 0;
1477 <      ferror_spl1_mct_[i] = 0;
1478 <      fres_spl2_mct_[i] = 0;
1479 <      ferror_spl2_mct_[i] = 0;
1480 <      fres_mct_[i] = 0;
1481 <      ferror_mct_[i] = 0;
1466 >  nrecPV_twovtx_ = 0;
1467 >  nsimPV_ = 0;
1468 >
1469 >  nrecPV_mct_ = 0;
1470 >  nrecPV_spl1_mct_ = 0;
1471 >  nrecPV_spl2_mct_ = 0;
1472 >
1473 >  // Mininum separation between the secondary pvtxes and leading pvtx
1474 >  min_zsep_ = 9999.0;
1475 >  min_ntrksep_ = 9999.0;
1476 >  min_sumpt2sep_ = -9999.0;
1477 >  
1478 >  // Variables filled per Vertex
1479 >  for (int i = 0; i < nMaxPVs_; i++) {
1480 >    // matched two-vertices
1481 >    nTrkPV1_spl_twovtx_[i] = 0;
1482 >    ndofPV1_spl_twovtx_[i] = 0;
1483 >    normchi2PV1_spl_twovtx_[i] = 0;
1484 >    avgPtPV1_spl_twovtx_[i] = 0;
1485 >    errx1_spl_twovtx_[i] = 0;
1486 >    erry1_spl_twovtx_[i] = 0;
1487 >    errz1_spl_twovtx_[i] = 0;
1488 >
1489 >    nTrkPV2_spl_twovtx_[i] = 0;
1490 >    ndofPV2_spl_twovtx_[i] = 0;
1491 >    normchi2PV2_spl_twovtx_[i] = 0;
1492 >    avgPtPV2_spl_twovtx_[i] = 0;
1493 >    errx2_spl_twovtx_[i] = 0;
1494 >    erry2_spl_twovtx_[i] = 0;
1495 >    errz2_spl_twovtx_[i] = 0;
1496 >    
1497 >    deltax_twovtx_[i] = 0;
1498 >    deltay_twovtx_[i] = 0;
1499 >    deltaz_twovtx_[i] = 0;
1500 >    errx_twovtx_[i] = 0;
1501 >    erry_twovtx_[i] = 0;
1502 >    errz_twovtx_[i] = 0;
1503 >    pullx_twovtx_[i] = 0;
1504 >    pully_twovtx_[i] = 0;
1505 >    pullz_twovtx_[i] = 0;
1506 >    
1507 >    //simpv
1508 >    nsimTrkPV_[i] = 0;
1509 >    simx_[i] = 0;
1510 >    simy_[i] = 0;
1511 >    simz_[i] = 0;
1512 >    simptsq_[i] = 0;
1513 >    
1514 >    // residual and pulls with mc truth required
1515 >    deltax_mct_[i] = 0;
1516 >    deltay_mct_[i] = 0;
1517 >    deltaz_mct_[i] = 0;
1518 >    pullx_mct_[i] = 0;
1519 >    pully_mct_[i] = 0;
1520 >    pullz_mct_[i] = 0;
1521 >    nTrkPV_mct_[i] = 0;
1522 >    ndofPV_mct_[i] = 0;
1523 >    normchi2PV_mct_[i] = 0;
1524 >    avgPtPV_mct_[i] = 0;
1525 >    errx_mct_[i] = 0;
1526 >    erry_mct_[i] = 0;
1527 >    errz_mct_[i] = 0;
1528 >    
1529 >    deltax_spl1_mct_[i] = 0;
1530 >    deltay_spl1_mct_[i] = 0;
1531 >    deltaz_spl1_mct_[i] = 0;
1532 >    pullx_spl1_mct_[i] = 0;
1533 >    pully_spl1_mct_[i] = 0;
1534 >    pullz_spl1_mct_[i] = 0;
1535 >    nTrkPV_spl1_mct_[i] = 0;
1536 >    ndofPV_spl1_mct_[i] = 0;
1537 >    normchi2PV_spl1_mct_[i] = 0;
1538 >    avgPtPV_spl1_mct_[i] = 0;
1539 >    errx_spl1_mct_[i] = 0;
1540 >    erry_spl1_mct_[i] = 0;
1541 >    errz_spl1_mct_[i] = 0;
1542 >
1543 >    deltax_spl2_mct_[i] = 0;
1544 >    deltay_spl2_mct_[i] = 0;
1545 >    deltaz_spl2_mct_[i] = 0;
1546 >    pullx_spl2_mct_[i] = 0;
1547 >    pully_spl2_mct_[i] = 0;
1548 >    pullz_spl2_mct_[i] = 0;
1549 >    nTrkPV_spl2_mct_[i] = 0;
1550 >    ndofPV_spl2_mct_[i] = 0;
1551 >    normchi2PV_spl2_mct_[i] = 0;
1552 >    avgPtPV_spl2_mct_[i] = 0;
1553 >    errx_spl2_mct_[i] = 0;
1554 >    erry_spl2_mct_[i] = 0;
1555 >    errz_spl2_mct_[i] = 0;
1556 >
1557 >  }
1558 > }
1559 >
1560 > double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1561 >  if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1562 >  else {
1563 >    double sum = 0.;
1564 >    double pT;
1565 >    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1566 >      pT = (**it).pt();
1567 >      sum += pT*pT;
1568      }
1569 +    return sum;
1570 +  }
1571 + }
1572 +
1573 +
1574 +
1575 + double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1576 +
1577 +  if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1578 +  else {
1579 +    double sumpT = 0.;
1580 +    for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1581 +      sumpT +=  (**it).pt();
1582 +  }
1583 +    return sumpT/double(v.tracksSize());
1584 +  }
1585   }
1586  
1587  
1588 +

Diff Legend

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