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.18 by jengbou, Fri Jan 29 21:36:22 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  
56   //root
57   #include <TROOT.h>
# Line 59 | Line 62 | Implementation:
62   #include <TPad.h>
63  
64   using namespace std;
65 + typedef math::XYZTLorentzVectorF LorentzVector;
66 + typedef math::XYZPoint Point;
67  
68   PVStudy::PVStudy(const edm::ParameterSet& iConfig)
69   {
# Line 80 | Line 85 | PVStudy::PVStudy(const edm::ParameterSet
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 > //   applyEvtClean_             = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
91 >  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
92 >  
93    //now do what ever initialization is needed
94    pvinfo.clear();
95  
87 –
96    // Specify the data mode vector
97    if(realData_) datamode.push_back(0);
98    else {
# Line 93 | Line 101 | PVStudy::PVStudy(const edm::ParameterSet
101      datamode.push_back(2);
102      datamode.push_back(3);
103    }
104 <
97 <  // Create ntuple files if needed
104 > // Create ntuple files if needed
105    if(saventuple_) {
106      file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107      ftree_ = new TTree("mytree","mytree");
# Line 102 | Line 109 | PVStudy::PVStudy(const edm::ParameterSet
109      ftree_->Branch("residual",&fres_,"fres_[3]/D");
110      ftree_->Branch("error",&ferror_,"ferror_[3]/D");
111      ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
112 <    ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
113 <    ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
114 <    ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
115 <  
116 <    // If it is MC, get the resA, resB
117 <    if(!realData_) {
118 <      ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
119 <      ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
120 <      ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
121 <      ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
122 <      ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
123 <      ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
124 <      ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
125 <      ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
126 <      ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
112 >
113 >    // pvtxtree_ analyzing the pvtxs ootb
114 >    pvtxtree_ = new TTree("pvtxtree","pvtxtree");
115 >    //pvtx using all tracks
116 >    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
117 >    pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
118 >    pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");    
119 >    pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
120 >    pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
121 >    pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
122 >    pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
123 >    pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
124 >    pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
125 >    pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
126 >    pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
127 >
128 >    pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
129 >    pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
130 >    pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
131 >
132 >    //pvtx using splittracks1
133 >    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
134 >    pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");  
135 >    pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");  
136 >    pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
137 >    pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
138 >    pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
139 >    pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
140 >    pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
141 >    pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
142 >    pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
143 >    pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");  
144 >  
145 >    //pvtx using splittracks2
146 >    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
147 >    pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");    
148 >    pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");  
149 >    pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
150 >    pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
151 >    pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
152 >    pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
153 >    pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
154 >    pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
155 >    pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
156 >    pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");  
157 >    
158 >    //pixeVertices
159 >    pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
160 >    pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");    
161 >    pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");  
162 >    pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
163 >    pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
164 >    pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
165 >    pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
166 >    pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
167 >    pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
168 >    pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
169 >    pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");  
170 >
171 >    //Fill the variables in the twovtx pair (recvtx1, recvtx2)
172 >    pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
173 >    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
174 >    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
175 >    pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
176 >    pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
177 >    pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
178 >    pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
179 >    pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
180 >    pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
181 >    pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
182 >    pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
183 >    pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
184 >    pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
185 >    pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
186 >    pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
187 >    
188 >    // MC variables
189 >    if(!realData_) {  
190 >      pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
191 >      pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
192 >      pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
193 >      pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
194 >      pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
195 >      pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
196 >
197 >      // For pvtxs with all the tracks
198 >      pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
199 >      pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
200 >      pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
201 >      pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
202 >      pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
203 >      pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
204 >      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
205 >      // For pvtxs with splittracks1
206 >      pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
207 >      pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
208 >      pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
209 >      pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
210 >      pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
211 >      pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
212 >      pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
213 >      // For pvtxs with splittracks1
214 >      pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
215 >      pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
216 >      pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
217 >      pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
218 >      pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
219 >      pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
220 >      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
221      }
222    }
223 <
123 <  edm::Service<TFileService> fs;
124 <  TFileDirectory subDir = fs->mkdir( "Summary" );
125 <  TFileDirectory subDir1 = fs->mkdir( "Others" );
126 <  //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
127 <
223 >  
224    setRootStyle();
225 +
226 +  //========================================
227 +  // Booking histograms
228 +  //========================================
229    
230 <  //Book histograms for track and pvtx parameters for each splitted pvtx
231 <  //i=0 : x (as in unsplit track collection)
232 <  //i=1 : x1_spl_
233 <  //i=2 : x2_spl_
234 <  
235 <  for(int i=0;i<3;i++)
236 <    {  
237 <      string suffix;
238 <      if(i == 0) suffix = "";
239 <      else {
240 <        stringstream ss;
241 <        ss << i;
242 <        suffix =ss.str()+"_spl";  
243 <      }
244 <      h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
245 <      h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
246 <      h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
247 <      h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
248 <      h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);  
249 <      h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);  
250 <      
251 <      h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
252 <      h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
253 <      h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
254 <      h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
155 <      h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);  
156 <      h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
157 <      h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);    
158 <      suffix.clear();
230 >  // Create a root file for histograms
231 >  theFile = new TFile(histoFileName_.c_str(), "RECREATE");
232 >  // make diretories
233 >  theFile->mkdir("Summary");
234 >  theFile->cd();
235 >  theFile->mkdir("Others");
236 >  theFile->cd();
237 >  if (analyze_) {
238 >    theFile->mkdir("Results");
239 >    theFile->cd();
240 >  }
241 >
242 >  // Book Histograms:
243 >  h_pvtrk   = new PVHistograms();
244 >  h_pixvtx  = new PVHistograms();
245 >  h_misc    = new PVHistograms();
246 >  h_summary = new PVHistograms();
247 >  h_others  = new PVHistograms();
248 >
249 >  for(int i=0;i<3;i++) {  
250 >    if(i == 0) h_pvtrk->Init("pvTrk");
251 >    else {
252 >      stringstream ss;
253 >      ss << i;
254 >      h_pvtrk->Init("pvTrk", ss.str(),"spl");
255      }
256 <  h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
257 <  h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
258 <  h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
259 <  
256 >  }
257 >  h_pixvtx->Init("pixVtx");
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_pixvtx->Save();
307 +  h_misc->Save();
308 +  h_summary->Save();
309 +  if (!realData_)
310 +    h_gen->Save();
311 +  if (analyze_) {
312 +    theFile->cd();
313 +    theFile->cd("Results");
314 +    h_ana->Save();
315 +  }
316 +  theFile->cd();
317 +  theFile->cd("Others");
318 +  h_others->Save();
319  
320 +  theFile->Close();
321   }
322  
323   void PVStudy::setRootStyle() {
# Line 292 | Line 354 | void PVStudy::setRootStyle() {
354   //
355   std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
356   {
357 <  std::vector<PVStudy::simPrimaryVertex> simpv;
357 > std::vector<PVStudy::simPrimaryVertex> simpv;
358    const HepMC::GenEvent* evt=evtMC->GetEvent();
359    if (evt) {
360 <    if(verbose_) {
361 <      cout << "process id " << evt->signal_process_id()<<endl;
362 <      cout << "signal process vertex " << ( evt->signal_process_vertex() ?
363 <                                            evt->signal_process_vertex()->barcode() : 0 )   <<endl;
302 <      cout << "number of vertices " << evt->vertices_size() << endl;
303 <    }
360 >    edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
361 >    edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
362 >                                                                        evt->signal_process_vertex()->barcode() : 0 ) <<endl;
363 >    edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
364  
365      int idx=0; int npv=0;
366      for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
# Line 311 | Line 371 | std::vector<PVStudy::simPrimaryVertex> P
371        // t component of PV:
372        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
373              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
374 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
374 >        //        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
375          HepMC::GenVertex * mv=(*mother)->production_vertex();
376          if( ((*mother)->status() == 3) && (!mv)) {
377 <          //        std::cout << "npv= " << npv << std::endl;
377 >          //        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
378            if (npv == 0) {
379 <            h["genPart_cT"]->Fill(pos.t()); // mm
380 <            h["genPart_T"]->Fill(pos.t()/299.792458); // ns
379 >            h_gen->Fill1d("genPart_cT", pos.t()); // mm
380 >            h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
381            }
382            npv++;
383          }
# Line 325 | Line 385 | std::vector<PVStudy::simPrimaryVertex> P
385        //       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
386          
387        bool hasMotherVertex=false;
388 <      if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
388 >      if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
389        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
390              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
391          HepMC::GenVertex * mv=(*mother)->production_vertex();
392 <        //        std::cout << "Status = " << (*mother)->status() << std::endl;
392 >        // if (verbose_) cout << "Status = " << (*mother)->status() << endl;
393          if (mv) {
394            hasMotherVertex=true;
395            if(!verbose_) break; //if verbose_, print all particles of gen vertices
# Line 356 | Line 416 | std::vector<PVStudy::simPrimaryVertex> P
416  
417        if(!vp){
418          // this is a new vertex
419 <        if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
419 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
420          simpv.push_back(sv);
421          vp=&simpv.back();
422        }else{
423 <        if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
423 >        edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
424        }
425        vp->genVertex.push_back((*vitr)->barcode());
426        // collect final state descendants
# Line 384 | Line 444 | std::vector<PVStudy::simPrimaryVertex> P
444              }
445            }
446          }
447 <      }
448 <    }
447 >      }//loop MC vertices daughters
448 >    }//loop MC vertices
449    }
450    return simpv;
451   }
# Line 418 | Line 478 | std::vector<PVStudy::simPrimaryVertex> P
478                ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
479              }
480            }
481 <          primary =  ( gv->position().t()==0);
481 >          //primary =  ( gv->position().t()==0 );
482 >          primary = true;
483 >          h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
484 >          h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
485 >
486            //resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
487            // no more mother pointer in the improved HepMC GenParticle
488            resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
# Line 467 | Line 531 | PVStudy::analyze(const edm::Event& iEven
531    using namespace edm;
532    using namespace reco;
533    
534 <  if (verbose_)
471 <    cout<<"PVStudy::analyze():"<<endl;
534 >  edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
535    //=======================================================
536    // Initialize Root-tuple variables if needed
537    //=======================================================
# Line 486 | Line 549 | PVStudy::analyze(const edm::Event& iEven
549    if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
550      trackColl = trackCollectionHandle.product();
551    } else {
552 <    cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
552 >    edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
553    }
554    //splitTrackCollection1
555    static const reco::TrackCollection s_empty_splitTrackColl1;
# Line 496 | Line 559 | PVStudy::analyze(const edm::Event& iEven
559    if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
560      splitTrackColl1 = splitTrackCollection1Handle.product();
561    } else {
562 <    cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
562 >    edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
563    }
564    //splitTrackCollection2
565    static const reco::TrackCollection s_empty_splitTrackColl2;
# Line 506 | Line 569 | PVStudy::analyze(const edm::Event& iEven
569    if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
570      splitTrackColl2 = splitTrackCollection2Handle.product();
571    } else {
572 <    cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
572 >   edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
573    }
574    
575    //=======================================================
# Line 520 | Line 583 | PVStudy::analyze(const edm::Event& iEven
583    if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
584      vertexColl = vertexCollectionHandle.product();
585    } else {
586 <    cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
586 >   edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
587    }
588    //splitVertexCollection1
589    static const reco::VertexCollection s_empty_splitVertexColl1;
# Line 530 | Line 593 | PVStudy::analyze(const edm::Event& iEven
593    if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
594      splitVertexColl1 = splitVertexCollection1Handle.product();
595    } else {
596 <    cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
596 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
597    }
598    //splitVertexCollection2
599    static const reco::VertexCollection s_empty_splitVertexColl2;
# Line 540 | Line 603 | PVStudy::analyze(const edm::Event& iEven
603    if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
604      splitVertexColl2 = splitVertexCollection2Handle.product();
605    } else {
606 <    cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
606 >    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
607    }
608  
609 <  if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
609 >  
610 >  //=======================================================
611 >  // GET pixelVertices
612 >  //=======================================================
613 >  static const reco::VertexCollection s_empty_pixelVertexColl;
614 >  const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
615 >  edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
616 >  iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
617 >  if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
618 >    pixelVertexColl = pixelVertexCollectionHandle.product();
619 >  } else {
620 >    edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
621 >  }
622 >
623 >  //=======================================================
624 >  // BeamSpot accessors
625 >  //=======================================================
626 >
627 >  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
628 >  iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
629 >  reco::BeamSpot bs = *recoBeamSpotHandle;    
630 >  const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
631 >  
632 >  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
633    
634    //=======================================================
635    // MC simvtx accessor
# Line 562 | Line 648 | PVStudy::analyze(const edm::Event& iEven
648    try{
649      iSetup.getData(pdt);
650    }catch(const Exception&){
651 <    cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
651 >    edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
652    }
653  
654    //=======================================================
655 <  // Fill trackparameters of the input tracks to pvtx fitter
655 >  // Apply event cleaning for firstcoll data and MC
656    //=======================================================
657  
658 <  // GeneralTracks  
659 <  h["nTrk"]->Fill(trackColl->size());
660 <  for(TrackCollection::const_iterator itTrack = trackColl->begin();
661 <      itTrack != trackColl->end();                      
662 <      ++itTrack) {
663 <    h["trkPt"]->Fill(itTrack->pt());
664 <    h["trkDxy"]->Fill(itTrack->dxy());
665 <    h["trkDz"]->Fill(itTrack->dz());
580 <    h["trkEta"]->Fill(itTrack->eta());
581 <    h["trkPhi"]->Fill(itTrack->phi());
582 <  }
583 <  
584 <  //SplittedTracks1
585 <  h["nTrk1_spl"]->Fill(splitTrackColl1->size());
586 <  for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
587 <      itTrack != splitTrackColl1->end();                      
588 <      ++itTrack) {
589 <    h["trkPt1_spl"]->Fill(itTrack->pt());
590 <    h["trkDxy1_spl"]->Fill(itTrack->dxy());
591 <    h["trkDz1_spl"]->Fill(itTrack->dz());
592 <    h["trkEta1_spl"]->Fill(itTrack->eta());
593 <    h["trkPhi1_spl"]->Fill(itTrack->phi());
594 <  }
595 <  //SplittedTracks2  
596 <  h["nTrk2_spl"]->Fill(splitTrackColl2->size());
597 <  for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
598 <      itTrack != splitTrackColl2->end();                      
599 <      ++itTrack) {
600 <    h["trkPt2_spl"]->Fill(itTrack->pt());
601 <    h["trkDxy2_spl"]->Fill(itTrack->dxy());
602 <    h["trkDz2_spl"]->Fill(itTrack->dz());
603 <    h["trkEta2_spl"]->Fill(itTrack->eta());
604 <    h["trkPhi2_spl"]->Fill(itTrack->phi());
658 >  // =====Loop over the trackColl to tet the fraction of HighPurity tracks
659 >  int nTracks = 0;
660 >  int nHighPurityTracks=0;
661 >  double fHighPurity=0;
662 >
663 >  for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
664 >    TrackRef tkref(trackCollectionHandle,i);
665 >    if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
666    }
667 <  
667 >  if(nTracks>0)
668 >    fHighPurity =  double(nHighPurityTracks)/double(nTracks);
669    if(verbose_)
670 <    cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
670 >    cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
671 >    
672 >  if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
673 >    //  glob_runno_ = iEvent.id().run();
674 >    //  glob_evtno_ = iEvent.id().event();
675 >    //  glob_ls_   = iEvent.luminosityBlock();
676 >    //  glob_bx_  = iEvent.bunchCrossing();  
677 >    return;
678 >  }
679 >
680 >  //=======================================================
681 >  // Fill trackparameters of the input tracks to pvtx fitter
682 >  //=======================================================
683 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
684 >  //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
685 >  // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
686 >  fillTrackHisto(trackColl, 0, beamSpot);
687 >  fillTrackHisto(splitTrackColl1, 1, beamSpot);
688 >  fillTrackHisto(splitTrackColl2, 2, beamSpot);
689 >  edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
690 >
691 >
692 >  //=======================================================
693 >  // Fill pixelVertices related histograms
694 >  //=======================================================
695 >  nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
696 >  if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
697 >    //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
698 >    fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
699 >    h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());    
700 >    // Get the dZ error of the tracks in the leading pixelVertexColl
701 >    for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
702 >        t!= (pixelVertexColl->begin())->tracks_end(); t++) {
703 >      
704 >      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
705 >        //        h_pvtrk->Fill1d("trkPtPV", 0.);
706 >      }
707 >      else
708 >        h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
709 >    }
710 >    // Fill the track->dz() in the PV difference from first pixelpvtx
711 >    for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
712 >        t!= (vertexColl->begin())->tracks_end(); t++) {
713 >      // illegal charge
714 >      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
715 >        //        h_pvtrk->Fill1d("trkPtPV", 0.);
716 >      }
717 >      else {
718 >        if(pixelVertexColl->size()>0) {
719 >          h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
720 >          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
721 >          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
722 >          h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
723 >        }
724 >      }
725 >    }
726 >  }
727  
728    //=======================================================
729 <  // Fill reconstructed vertices
729 >  // Fill number of reconstructed vertices
730    //=======================================================
731 <  if(verbose_)  {
732 <    cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
731 >
732 >  edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
733 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
734 >  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
735 >  edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
736 >  if (verbose_) {
737      printRecVtxs(vertexCollectionHandle);
616 –    cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
738      printRecVtxs(splitVertexCollection1Handle);
618 –    cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
739      printRecVtxs(splitVertexCollection2Handle);
740    }
741 +  nrecPV_ = int (vertexColl->size());
742 +  nrecPV1_spl_ = int (splitVertexColl1->size());
743 +  nrecPV2_spl_ = int (splitVertexColl2->size());
744 +
745 +  h_pvtrk->Fill1d("nrecPV", nrecPV_);
746 +  h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
747 +  h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
748 +  h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
749    
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_));
750    
751 <  // Fill the track paramters for the vertexColl
752 <  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
753 <      v!=vertexColl->end(); ++v){
754 <    h["nTrkPV"]->Fill(v->tracksSize());
755 <    try {
756 <      for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
757 <          t!=v->tracks_end(); t++) {
758 <        // illegal charge
759 <        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
760 <          //      h["trkPtPV"]->Fill(0.);
761 <        }
762 <        else {
763 <          h["trkPtPV"]->Fill((**t).pt());
764 <          h["trkDxyPV"]->Fill((**t).dxy());
765 <          h["trkDzPV"]->Fill((**t).dz());
766 <          h["trkEtaPV"]->Fill((**t).eta());
767 <          h["trkPhiPV"]->Fill((**t).phi());
768 <        }
751 > //=================================================================
752 >  // Fill track parameter ntuple/hist for tracks used in recoVertices
753 >  //=================================================================
754 >
755 >  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
756 >  if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
757 >    fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
758 >  
759 >  if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
760 >    fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
761 >
762 >  if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
763 >    fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
764 >
765 >  //=======================================================
766 >  // Compare offlinePrimaryVertices with pixelVertices
767 >  //=======================================================
768 >  if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
769 >      && (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {    
770 >    h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
771 >    // difference in reconstructed position of the leading pvtx
772 >    //edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;  
773 >    //edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
774 >    h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
775 >    h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
776 >    h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
777 >  }
778 >  
779 >  
780 > //==========================================================
781 >  // Fill secondary/primary min separations for multi-vertices
782 >  //==========================================================
783 >  
784 >  if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {  
785 >
786 >    double min_xsep = 9999.0;
787 >    double min_ysep = 9999.0;
788 >    double min_zsep = 9999.0;    
789 >    double min_xsep_sign = 9999.0;
790 >    double min_ysep_sign = 9999.0;
791 >    double min_zsep_sign = 9999.0;
792 >    double min_ntrksep = 9999.0;
793 >    double min_sumpt2sep = 9999999.0;
794 >    
795 >    edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
796 >
797 >    // Looping through the secondary vertices to find the mininum separation to the primary
798 >    for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
799 >        v!=vertexColl->end(); ++v) {  
800 >      if(v->isValid() && ! v->isFake()) {
801 >        // xsep
802 >        if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
803 >          min_xsep = fabs(v->x()- vertexColl->begin()->x());
804 >        // ysep
805 >        if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
806 >          min_ysep = fabs(v->y()- vertexColl->begin()->y());
807 >        // zsep
808 >        if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
809 >          min_zsep = fabs(v->z()- vertexColl->begin()->z());
810 >        // xsep_sign
811 >        double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
812 >        if(xsep_sign < min_xsep_sign)
813 >          min_xsep_sign = xsep_sign;
814 >        // ysep_sign
815 >        double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
816 >        if(ysep_sign < min_ysep_sign)
817 >          min_ysep_sign = ysep_sign;
818 >        // zsep_sign
819 >        double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
820 >        if(zsep_sign < min_zsep_sign)
821 >          min_zsep_sign = zsep_sign;
822 >        // ntrksep
823 >        if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
824 >          min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
825 >        // sumpt2sep
826 >        if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
827 >          min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
828        }
829      }
830 <    catch (...) {
831 <      // exception thrown when trying to use linked track
832 <      //          h["trkPtPV"]->Fill(0.);
833 <    }
834 <  }
830 >    h_misc->Fill1d("min_xsep", min_xsep);
831 >    h_misc->Fill1d("min_ysep", min_ysep);
832 >    h_misc->Fill1d("min_zsep", min_zsep);
833 >    h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
834 >    h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
835 >    h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
836 >    h_misc->Fill1d("min_ntrksep", min_ntrksep);
837 >    h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
838 >
839 >    min_zsep_ = min_zsep;
840 >    min_ntrksep_ = min_ntrksep;
841 >    min_sumpt2sep_ = min_sumpt2sep;
842 >
843 >
844 >  } // End of  if(vertexColl->size()>1) {
845 >  
846 >  edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
847 >
848 >  
849 >    
850 >  //=======================================================
851 >  //           PrimaryVertex Matching
852 >  // 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
853 >  // 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
854 >  // Assume the first match is the primary vertex,
855 >  // since vertexColl are sorted in decreasing order of sum(pT)
856 >  //=======================================================
857    
858 +  edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
859 +  reco::VertexCollection s_empty_matchedVertexColl1;
860 +  reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
861 +  matchedVertexColl1->reserve(splitVertexColl1->size());
862 +  reco::VertexCollection s_empty_matchedVertexColl2;
863 +  reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
864 +  matchedVertexColl2->reserve(splitVertexColl2->size());
865    
866 <  // == ignore multi-vertices (for now)
867 <  //take single pvt collection and compare
868 <  if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
869 <    const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
870 <    const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
871 <    if (matchVertex(recvtx1, recvtx2))  {
872 <      int nTrkPV1 = recvtx1.tracksSize();
873 <      int nTrkPV2 = recvtx2.tracksSize();    
874 <      h["nTrkPV1_spl"]->Fill(nTrkPV1);
875 <      h["nTrkPV2_spl"]->Fill(nTrkPV2);
876 <      double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
877 <      h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
878 <      h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
866 >  bool stopmatching = false;
867 >  for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
868 >    reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
869 >    if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
870 >    for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {    
871 >      //------------------------------------------------------------------------
872 >      //== If only considering matching the first vertex of the splitVertexColl
873 >      //------------------------------------------------------------------------
874 >      if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
875 >      reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
876 >      if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
877 >      edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
878 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
879 >      edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
880        
881 <      // Fill the track paramters for the splitVertexColl1
882 <      try {
883 <        for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
884 <            t!=recvtx1.tracks_end(); t++) {
885 <          // illegal charge
886 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
887 <            //    h["trkPtPV"]->Fill(0.);
888 <          }
889 <          else {
890 <            h["trkPtPV1_spl"]->Fill((**t).pt());
891 <            h["trkDxyPV1_spl"]->Fill((**t).dxy());
892 <            h["trkDzPV1_spl"]->Fill((**t).dz());
893 <            h["trkEtaPV1_spl"]->Fill((**t).eta());
894 <            h["trkPhiPV1_spl"]->Fill((**t).phi());
895 <          }
881 >      double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
882 >      h_misc->Fill1d("twovtxzsign", twovtxsig);
883 >      if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
884 >        edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
885 >
886 >        int nTrkPV1 = recvtx1->tracksSize();
887 >        int nTrkPV2 = recvtx2->tracksSize();    
888 >        double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
889 >        h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
890 >        h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
891 >        if(fabs(ntrkreldiff)<ntrkdiffcut_) {    
892 >          edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
893 >
894 >          matchedVertexColl1->push_back(*recvtx1);
895 >          matchedVertexColl2->push_back(*recvtx2);
896 >          stopmatching = true; // stop the matching when the first match is found!
897 >          break;
898          }
899 +        else
900 +          edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
901        }
902 <      catch (...) {
903 <        // exception thrown when trying to use linked track
904 <        //        h["trkPtPV"]->Fill(0.);
902 >      else {
903 >        edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
904 >        edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
905 >        edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;  
906        }
907 +    }
908 +  }
909 +  edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
910 +  
911 +  //=======================================================
912 +  //  Analyze matchedVertexColls
913 +  //=======================================================
914 +  edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
915 +  
916 +  // Compare the reconstructed vertex position and calculate resolution/pulls
917 +  if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {  
918 +    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
919 +    fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
920 +    fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
921 +
922 +    reco::VertexCollection::const_iterator v1;
923 +    reco::VertexCollection::const_iterator v2;
924 +    nrecPV_twovtx_ = 0;
925 +    for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
926 +        v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
927 +        ++v1, ++v2) {
928 +      h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
929 +      h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
930 +      fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
931 +      fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
932 +      fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
933 +      ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
934 +      ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
935 +      ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
936 +      fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
937        
938 <      // Fill the track paramters for the splitVertexColl2
939 <      try {
940 <        for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
941 <            t!=recvtx2.tracks_end(); t++) {
942 <          // illegal charge
943 <          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
944 <            //    h["trkPtPV"]->Fill(0.);
945 <          }
946 <          else {
947 <            h["trkPtPV2_spl"]->Fill((**t).pt());
948 <            h["trkDxyPV2_spl"]->Fill((**t).dxy());
949 <            h["trkDzPV2_spl"]->Fill((**t).dz());
950 <            h["trkEtaPV2_spl"]->Fill((**t).eta());
951 <            h["trkPhiPV2_spl"]->Fill((**t).phi());
952 <          }
953 <        }
954 <      }
955 <      catch (...) {
956 <        // exception thrown when trying to use linked track
957 <        //        h["trkPtPV"]->Fill(0.);
938 >      h_summary->Fill1d("deltax", fres_[0] );
939 >      h_summary->Fill1d("deltay", fres_[1] );
940 >      h_summary->Fill1d("deltaz", fres_[2] );
941 >      h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
942 >      h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
943 >      h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
944 >      h_summary->Fill1d("errPVx", ferror_[0] );
945 >      h_summary->Fill1d("errPVy", ferror_[1] );
946 >      h_summary->Fill1d("errPVz", ferror_[2] );
947 >      pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
948 >                                       error(ferror_[0], ferror_[1], ferror_[2]),
949 >                                       fntrk_));
950 >      // Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
951 >      fillHisto(res(fres_[0], fres_[1], fres_[2]),
952 >                error(ferror_[0], ferror_[1], ferror_[2]),
953 >                fntrk_,0);  
954 >      // Fill the ntuple variables
955 >      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
956 >      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
957 >      sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
958 >      sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
959 >      nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
960 >      deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
961 >      deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
962 >      deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
963 >      errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
964 >      erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
965 >      errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
966 >      pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
967 >      pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
968 >      pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];    
969 >      nrecPV_twovtx_++;
970 >    } // End of analyzing res/pull
971 >
972 >    //=======================================================
973 >    // Fill simulated vertices
974 >    //=======================================================
975 >    if (!realData_) {
976 >      bool MC=false;
977 >      Handle<HepMCProduct> evtMC;
978 >      iEvent.getByLabel("generator",evtMC);
979 >      if (!evtMC.isValid()) {
980 >        MC=false;
981 >        edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
982 >      } else {
983 >        edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
984 >        MC=true;
985        }
986        
987 <      if(verbose_)
988 <        cout<<"Done filling track parameters in reconstruced vertices."<<endl;
989 <
990 <      if(fabs(ntrkreldiff)<ntrkdiffcut_) {
991 <        fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
992 <        fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
993 <        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;
987 >      if(MC){
988 >        // make a list of primary vertices:
989 >        std::vector<simPrimaryVertex> simpv;
990 >        simpv=getSimPVs(evtMC,"");
991 >        //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
992 >        h_gen->Fill1d("nsimPV", simpv.size());
993 >        nsimPV_ = simpv.size();
994          
995 <        h["deltax"]->Fill( fres_[0] );
996 <        h["deltay"]->Fill( fres_[1] );
997 <        h["deltaz"]->Fill( fres_[2] );
998 <        h["pullx"]->Fill( fres_[0]/ferror_[0]);
999 <        h["pully"]->Fill( fres_[1]/ferror_[1]);
1000 <        h["pullz"]->Fill( fres_[2]/ferror_[2]);
1001 <        h["errPVx"]->Fill( ferror_[0] );
1002 <        h["errPVy"]->Fill( ferror_[1] );
1003 <        h["errPVz"]->Fill( ferror_[2] );
1004 <        pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
1005 <                                         error(ferror_[0], ferror_[1], ferror_[2]),
1006 <                                         fntrk_));
1007 <        // Fill histo according to its track multiplicity, set datamode = 0
1008 <        fillHisto(res(fres_[0], fres_[1], fres_[2]),
1009 <                  error(ferror_[0], ferror_[1], ferror_[2]),
1010 <                  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;
752 <    }
753 <  } // End of fill histograms for two-vertices Method
995 >        int isimpv = 0;
996 >        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
997 >            vsim!=simpv.end(); vsim++, isimpv++){
998 >          nsimTrkPV_[isimpv]  =vsim->nGenTrk;
999 >          simx_[isimpv] = vsim->x;
1000 >          simy_[isimpv] = vsim->y;
1001 >          simz_[isimpv] = vsim->y;
1002 >          simptsq_[isimpv] = vsim->ptsq;
1003 >          
1004 >          //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1005 >          fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1006 >          fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1007 >          fillMCHisto(vsim, isimpv, vertexColl, 3);
1008 >        }
1009 >      } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1010 >    } // End of   if (!realData_) {
1011  
1012 <  if(verbose_)
1013 <    cout<<"Done filling res/pull with two-vertices method"<<endl;
1012 >  } // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1013 >  else
1014 >    edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1015 >  
1016 >  edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1017  
1018 <  //=======================================================
1019 <  // Fill simulated vertices
1020 <  //=======================================================
1021 <  if (!realData_) {
1022 <    bool MC=false;
1023 <    Handle<HepMCProduct> evtMC;
1024 <    iEvent.getByLabel("generator",evtMC);
1025 <    if (!evtMC.isValid()) {
766 <      MC=false;
767 <      if(verbose_){
768 <        cout << "no HepMCProduct found"<< endl;
769 <      }
770 <    } else {
771 <      if(verbose_){
772 <        cout << "generator HepMCProduct found"<< endl;
773 <      }
774 <      MC=true;
775 <    }
776 <    
777 <    if(MC){
778 <      // make a list of primary vertices:
779 <      std::vector<simPrimaryVertex> simpv;
780 <      simpv=getSimPVs(evtMC,"");
781 <      //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
782 <      h["nsimPV"]->Fill(simpv.size());
783 <      int nsimtrk=0;
784 <      int nrectrk=0;
785 <      
786 <      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
787 <          vsim!=simpv.end(); vsim++){
788 <        nsimtrk+=vsim->nGenTrk;
789 <        // Loop over the datamodes, for MC, fill only 1, 2, 3
790 <        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
883 <      }
884 <    } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
885 <  } // End of Filling MC variables
1018 >
1019 >
1020 >  // Finally fill the ftree_
1021 >  if(saventuple_) {
1022 >    ftree_->Fill();
1023 >    pvtxtree_->Fill();
1024 >  }
1025 >  
1026   }
1027  
1028 +
1029 +  
1030   // ------------ method called once each job just before starting event loop  ------------
1031   void
1032   PVStudy::beginJob()
# Line 894 | Line 1036 | PVStudy::beginJob()
1036   // ------------ method called once each job just after ending the event loop  ------------
1037   void
1038   PVStudy::endJob() {
1039 <  if (verbose_) cout << "Analyzing PV info" << endl;
1039 >  edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
1040    // Looping through the datamodes available
1041    for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {  
1042      string suffix;
1043 <    if(verbose_)  cout<<"datamode = "<< *it<<endl;
1043 >    edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1044      switch(*it) {
1045      case 0: suffix = "";
1046        break;
# Line 909 | Line 1051 | PVStudy::endJob() {
1051      case 3: suffix = "_mct";
1052        break;
1053      }
1054 <    if (analyze_) {
1054 >    suffix += "_nTrk";
1055 >    if(analyze_) {
1056        for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1057 +        edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1058          PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1059 <        if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1060 <        if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1061 <        if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1062 <        if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1063 <        if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1064 <        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.);
1059 >        if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1060 >        if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1061 >        if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1062 >        if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1063 >        if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1064 >        if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1065        }
1066      }
1067      suffix.clear();
# Line 926 | Line 1069 | PVStudy::endJob() {
1069    if(saventuple_) {
1070      file_->cd();
1071      ftree_->Write();
1072 +    pvtxtree_->Write();
1073      file_->Close();
1074    }
1075 < }
932 <
1075 >  
1076  
1077 + }
1078  
1079   // Match a recovertex with a simvertex
1080   // ? Should the cut be parameter dependent?
# Line 938 | Line 1082 | PVStudy::endJob() {
1082   bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1083                            const reco::Vertex & vrec)
1084   {
1085 <  return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1085 >  if(vrec.isFake() || !vrec.isValid()) return false;
1086 >  else
1087 >    return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1088   }
1089  
1090   // Match two reconstructed vertices
1091 < bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1092 <                           const reco::Vertex & vrec2)
1091 > bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1092 >                           const reco::VertexRef & vrec2,
1093 >                           int zsigncut)
1094   {
1095 <  return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1095 >  double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1096 >  edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1097 >  return (fabs(vrec1->z()-vrec2->z())<cut);
1098   }
1099  
1100  
1101   bool PVStudy::isResonance(const HepMC::GenParticle * p){
1102    double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1103 <  if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1103 >  edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1104    return  ctau >0 && ctau <1e-6;
1105   }
1106  
# Line 1023 | Line 1172 | void PVStudy::printSimTrks(const Handle<
1172   void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1173    stringstream ss;
1174    ss << ntk;
1175 <  string suffix = ss.str();
1176 <  if(datamode == 0 )  suffix = "_" + suffix;
1177 <  if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1178 <  if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1179 <  if(datamode == 3 )  suffix = "_mct_" + suffix;
1175 >  string suffix;
1176 >  if(datamode == 0 )  suffix = "_" + ss.str();
1177 >  if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1178 >  if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1179 >  if(datamode == 3 )  suffix = "_mct_" + ss.str();
1180    
1181    if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1182    // Fill histograms of res and pull of ntk
1183 <  h["deltax"+suffix]->Fill(r.x());
1184 <  h["deltay"+suffix]->Fill(r.y());
1185 <  h["deltaz"+suffix]->Fill(r.z());
1186 <  h["pullx"+suffix]->Fill(r.x()/er.x());
1187 <  h["pully"+suffix]->Fill(r.y()/er.y());
1188 <  h["pullz"+suffix]->Fill(r.z()/er.z());  
1189 <  h["errPVx"+suffix]->Fill(er.x());
1190 <  h["errPVy"+suffix]->Fill(er.y());
1191 <  h["errPVz"+suffix]->Fill(er.z());
1183 >  h_others->Fill1d("deltax"+suffix, r.x());
1184 >  h_others->Fill1d("deltay"+suffix, r.y());
1185 >  h_others->Fill1d("deltaz"+suffix, r.z());
1186 >  h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1187 >  h_others->Fill1d("pully"+suffix, r.y()/er.y());
1188 >  h_others->Fill1d("pullz"+suffix, r.z()/er.z());  
1189 >  h_others->Fill1d("errPVx"+suffix, er.x());
1190 >  h_others->Fill1d("errPVy"+suffix, er.y());
1191 >  h_others->Fill1d("errPVz"+suffix, er.z());
1192   }
1193  
1194 +
1195   PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1196    map<int,double> data;
1197    data.clear();
# Line 1059 | Line 1209 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1209      for ( int i = 0; i < 6; ++i) {
1210        switch (i) {
1211        case 0:
1212 <        fit(h["deltax"+suffix], fpar);
1212 >        if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1213 >        fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1214          data[i] = fpar[0]*cm2um;
1215          break;
1216        case 1:
1217 <        fit(h["deltay"+suffix], fpar);
1217 >        if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1218 >        fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1219          data[i] = fpar[0]*cm2um;
1220          break;
1221        case 2:
1222 <        fit(h["deltaz"+suffix], fpar);
1222 >        if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1223 >        fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1224          data[i] = fpar[0]*cm2um;
1225          break;
1226        case 3:
1227 <        fit(h["pullx"+suffix], fpar);
1227 >        if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1228 >        fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1229          data[i] = fpar[0];
1230          break;
1231        case 4:
1232 <        fit(h["pully"+suffix], fpar);
1232 >        if (!h_others->ReadHisto1D("pully"+suffix)) break;
1233 >        fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1234          data[i] = fpar[0];
1235          break;
1236        case 5:
1237 <        fit(h["pullz"+suffix], fpar);
1237 >        if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1238 >        fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1239          data[i] = fpar[0];
1240          break;
1241        }
1242 <      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1242 >      if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1243      }
1244    }
1245    else
# Line 1098 | Line 1254 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1254                              ntk);
1255   }
1256  
1257 < void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1257 >
1258 > void PVStudy::fit(TH1 *hdist, double fitPars[]){
1259    TAxis *axis0 = hdist->GetXaxis();
1260    int nbin = axis0->GetLast();
1261    double nOF = hdist->GetBinContent(nbin+1);
1262    double nUF = hdist->GetBinContent(0);
1263    //   double fitRange = axis0->GetBinUpEdge(nbin);
1264 <  double fitRange = axis0->GetXmax();
1264 >  // double fitRange = axis0->GetXmax();
1265 >  double fitRange = 2*hdist->GetRMS();
1266    double sigMax[2] = {0.,0.};
1267  
1268 <  if ( verbose_ ){
1269 <    cout << "Last bin = " << nbin;
1270 <    cout << "; hdist: Overflow Entries = " << nOF;
1271 <    cout << "; Underflow Entries = " << nUF;
1272 <    cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1273 <    cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1116 <  }
1268 >  edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1269 >                           << "; hdist: Overflow Entries = " << nOF
1270 >                           << "; Underflow Entries = " << nUF
1271 >                           << "; hdist->GetEntries() = " << hdist->GetEntries()
1272 >                           << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1273 >
1274    if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1275      for (int bi = 0; bi < nbin; bi++) {
1276        if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
# Line 1125 | Line 1282 | void PVStudy::fit(TH1 *&hdist, double fi
1282          if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1283        }
1284      }
1285 <    if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1285 >    //edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1286  
1287      TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1288      fgaus->SetParameter(1, 0.);
1289 <    fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1289 >    fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1290      fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1291      fgaus->SetLineColor(4);
1292 <    hdist->Fit(fgaus,"Q");
1292 >    hdist->Fit(fgaus,"QLRM");
1293      gPad->Update();
1294      TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1295      s->SetOptStat(1111111);
1296      s->SetOptFit(0111);
1297      gPad->Update();
1298 <    if (verbose_) cout << "got function" << endl;
1298 >    //edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1299      fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1300    }
1301    else
1302      fitPars[0] = -999;
1303   }
1304  
1305 +
1306 + void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1307 +  string suffix;
1308 +  suffix = ""; // for unsplit tracks
1309 +  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1310 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1311 +  
1312 +  h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1313 +  for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1314 +      itTrack != trackColl->end();                      
1315 +      ++itTrack) {
1316 +    h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1317 +    h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1318 +    h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1319 +    h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1320 +    h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1321 +    h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1322 +  }
1323 + }
1324 +
1325 + void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1326 +  string suffix;
1327 +  suffix = ""; // for unsplit tracks
1328 +  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1329 +  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1330 +  int ivtx = 0;
1331 +  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1332 +      v!=vertexColl->end(); ++v, ++ivtx) {  
1333 +    if(v->isFake()) continue;
1334 +    if(fillNtuple) {
1335 +      if(datatype == 0) {
1336 +        nTrkPV_[ivtx] = v->tracksSize();        
1337 +        sumptsq_[ivtx] = sumPtSquared(*v);
1338 +        isValid_[ivtx] = v->isValid();
1339 +        isFake_[ivtx] = v->isFake();    
1340 +        recx_[ivtx] = v->x();
1341 +        recy_[ivtx] = v->y();
1342 +        recz_[ivtx] = v->z();
1343 +        recx_err_[ivtx] = v->xError();
1344 +        recy_err_[ivtx] = v->yError();
1345 +        recz_err_[ivtx] = v->zError();
1346 +      }
1347 +      if(datatype == 1) {
1348 +        nTrkPV1_spl_[ivtx] = v->tracksSize();  
1349 +        sumptsq1_spl_[ivtx] = sumPtSquared(*v);        
1350 +        isValid1_spl_[ivtx] = v->isValid();
1351 +        isFake1_spl_[ivtx] = v->isFake();      
1352 +        recx1_spl_[ivtx] = v->x();
1353 +        recy1_spl_[ivtx] = v->y();
1354 +        recz1_spl_[ivtx] = v->z();
1355 +        recx1_err_spl_[ivtx] = v->xError();
1356 +        recy1_err_spl_[ivtx] = v->yError();
1357 +        recz1_err_spl_[ivtx] = v->zError();
1358 +
1359 +      }
1360 +      if(datatype == 2) {
1361 +        nTrkPV2_spl_[ivtx] = v->tracksSize();  
1362 +        sumptsq2_spl_[ivtx] = sumPtSquared(*v);        
1363 +        isValid2_spl_[ivtx] = v->isValid();
1364 +        isFake2_spl_[ivtx] = v->isFake();
1365 +        recx2_spl_[ivtx] = v->x();
1366 +        recy2_spl_[ivtx] = v->y();
1367 +        recz2_spl_[ivtx] = v->z();
1368 +        recx2_err_spl_[ivtx] = v->xError();
1369 +        recy2_err_spl_[ivtx] = v->yError();
1370 +        recz2_err_spl_[ivtx] = v->zError();
1371 +      }
1372 +      if(datatype == 4) { // for pixelVertices
1373 +        nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();        
1374 +        sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);      
1375 +        isValid_pxlpvtx_[ivtx] = v->isValid();
1376 +        isFake_pxlpvtx_[ivtx] = v->isFake();
1377 +        recx_pxlpvtx_[ivtx] = v->x();
1378 +        recy_pxlpvtx_[ivtx] = v->y();
1379 +        recz_pxlpvtx_[ivtx] = v->z();
1380 +        recx_err_pxlpvtx_[ivtx] = v->xError();
1381 +        recy_err_pxlpvtx_[ivtx] = v->yError();
1382 +        recz_err_pxlpvtx_[ivtx] = v->zError();
1383 +      }
1384 +    }
1385 +  }
1386 +  // For histogram only fill the leading pvtx parameters
1387 +  if(fillHisto) {
1388 +    h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1389 +    h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1390 +    int nHWTrkPV_ = 0;
1391 +    
1392 +    try {
1393 +      for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1394 +          t!=vertexColl->begin()->tracks_end(); t++) {
1395 +        // illegal charge
1396 +        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1397 +          //      h_pvtrk->Fill1d("trkPtPV", 0.);
1398 +        }
1399 +        else {
1400 +          h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1401 +          h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());  
1402 +          h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1403 +          h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1404 +          h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1405 +          h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1406 +          
1407 +          if(vertexColl->begin()->trackWeight(*t) > 0.5)
1408 +            nHWTrkPV_++;
1409 +        }
1410 +      }
1411 +    }
1412 +    catch (...) {
1413 +      // exception thrown when trying to use linked track
1414 +      //          h_pvtrk->Fill1d("trkPtPV", 0.);
1415 +    }
1416 +    h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1417 +  }
1418 +  
1419 +  
1420 + }
1421 +
1422 + void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1423 + {
1424 +  std::string suffix;
1425 +  // Set the vertexColl and histogram suffix according to datamode
1426 +  if (datatype == 1) {
1427 +    suffix = "_spl1_mct";
1428 +    nrecPV_spl1_mct_ = 0;  
1429 +  }
1430 +  if (datatype == 2) {
1431 +    suffix = "_spl2_mct";
1432 +    nrecPV_spl2_mct_ = 0;    
1433 +  }
1434 +  if (datatype == 3) {
1435 +    suffix = "_mct";
1436 +    nrecPV_mct_ = 0;  
1437 +  }
1438 +  
1439 +  //========================================================
1440 +  //  look for a matching reconstructed vertex in vertexColl
1441 +  //========================================================        
1442 +
1443 +  for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1444 +      vrec!=vtxColl->end(); ++vrec){
1445 +    int nrectrk = vrec->tracksSize();
1446 +    vsim->recVtx=NULL;  
1447 +  
1448 +    edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1449 +    edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1450 +
1451 +    if ( matchVertex(*vsim,*vrec)  ) {
1452 +      vsim->recVtx=&(*vrec);
1453 +      
1454 +      // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1455 +      //if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1456 +      //|| (!vsim->recVtx) )
1457 +      //vsim->recVtx=&(*vrec);
1458 +      //}
1459 +      edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1460 +      
1461 +      double fres_mct[3];
1462 +      double ferror_mct[3];
1463 +      
1464 +      fres_mct[0] = vsim->recVtx->x()-vsim->x;
1465 +      fres_mct[1] = vsim->recVtx->y()-vsim->y;
1466 +      fres_mct[2] = vsim->recVtx->z()-vsim->z;
1467 +      ferror_mct[0] = vsim->recVtx->xError();
1468 +      ferror_mct[1] = vsim->recVtx->yError();
1469 +      ferror_mct[2] = vsim->recVtx->zError();
1470 +      
1471 +      h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1472 +      h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1473 +      h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1474 +      h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1475 +      h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1476 +      h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1477 +      h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1478 +      h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1479 +      h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1480 +      pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1481 +                                       error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1482 +                                       nrectrk));
1483 +      // Fill histo according to its track multiplicity
1484 +      fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1485 +                error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1486 +                nrectrk, datatype);
1487 +      
1488 +      if(saventuple_) {
1489 +        //Fill the values for variables in ftree_  
1490 +        if(datatype == 1) {
1491 +          nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1492 +          deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1493 +          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1494 +          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1495 +          pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];  
1496 +          pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1497 +          pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1498 +          nrecPV_spl1_mct_++;
1499 +        }
1500 +        
1501 +        if(datatype == 2) {
1502 +          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1503 +          deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1504 +          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1505 +          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1506 +          pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];  
1507 +          pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1508 +          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1509 +          nrecPV_spl2_mct_++;
1510 +        }
1511 +        if(datatype == 3) {    
1512 +          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1513 +          deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1514 +          deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1515 +          deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1516 +          pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];  
1517 +          pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1518 +          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];          
1519 +          nrecPV_mct_++;
1520 +        }
1521 +      } // End of  if(saventuple_) {
1522 +    } //  if ( matchVertex(*vsim,*vrec) ) {
1523 +    else {  // no rec vertex found for this simvertex
1524 +      edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1525 +    }
1526 +  }
1527 + }
1528 +
1529   void PVStudy::SetVarToZero() {
1530 +  //Greg's variables
1531 +  fntrk_ = 0;
1532 +  //pvtx position (x,y,z) residual and error
1533 +  for(int i = 0; i<3;i++)    {
1534 +    fres_[i] = 0;
1535 +    ferror_[i] = 0;
1536 +  }
1537    
1538 <  // number of reconstructed vertices
1538 >  //Yanyan's variables
1539 >  // Number of reconstructed vertices  
1540    nrecPV_ = 0;
1541 <  nrecPV1_spl_ = 0;
1541 >  nrecPV1_spl_ = 0;
1542    nrecPV2_spl_ = 0;
1543 <  // number of tracks in the vertex
1544 <  fntrk_ = 0;
1545 <  fntrk_spl1_mct_ = 0;
1546 <  fntrk_spl2_mct_ = 0;
1547 <  fntrk_mct_ = 0;
1548 <  //pvtx position (x,y,z) residual and error
1549 <  for(int i = 0; i<3;i++)
1550 <    {
1551 <      fres_[i] = 0;
1552 <      ferror_[i] = 0;
1553 <      fres_spl1_mct_[i] = 0;
1554 <      ferror_spl1_mct_[i] = 0;
1555 <      fres_spl2_mct_[i] = 0;
1556 <      ferror_spl2_mct_[i] = 0;
1557 <      fres_mct_[i] = 0;
1558 <      ferror_mct_[i] = 0;
1559 <    }
1543 >  nrecPV_twovtx_ = 0;
1544 >  nsimPV_ = 0;
1545 >
1546 >  nrecPV_mct_ = 0;
1547 >  nrecPV_spl1_mct_ = 0;
1548 >  nrecPV_spl2_mct_ = 0;
1549 >
1550 >  // Mininum separation between the secondary pvtxes and leading pvtx
1551 >  min_zsep_ = 9999.0;
1552 >  min_ntrksep_ = 9999.0;
1553 >  min_sumpt2sep_ = -9999.0;
1554 >  
1555 >
1556 >  for (int i = 0; i < nMaxPVs_; i++) {
1557 >    // recoVertices with all tracks
1558 >    nTrkPV_[i] = 0; // Number of tracks in the pvtx    
1559 >    sumptsq_[i] = 0;
1560 >    isValid_[i] = -1;
1561 >    isFake_[i] = -1;
1562 >    recx_[i] = 0;
1563 >    recy_[i] = 0;
1564 >    recz_[i] = 0;
1565 >    recx_err_[i] = 0;
1566 >    recy_err_[i] = 0;
1567 >    recz_err_[i] = 0;
1568 >    
1569 >    // recoVertices with splitTrack1
1570 >    nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx    
1571 >    sumptsq1_spl_[i] = 0;
1572 >    isValid1_spl_[i] = -1;
1573 >    isFake1_spl_[i] = -1;
1574 >    recx1_spl_[i] = 0;
1575 >    recy1_spl_[i] = 0;
1576 >    recz1_spl_[i] = 0;
1577 >    recx1_err_spl_[i] = 0;
1578 >    recy1_err_spl_[i] = 0;
1579 >    recz1_err_spl_[i] = 0;
1580 >  
1581 >    // recoVertices with splitTrack2
1582 >    nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx  
1583 >    sumptsq2_spl_[i] = 0;
1584 >    isValid2_spl_[i] = -1;
1585 >    isFake2_spl_[i] = -1;
1586 >    recx2_spl_[i] = 0;
1587 >    recy2_spl_[i] = 0;
1588 >    recz2_spl_[i] = 0;
1589 >    recx2_err_spl_[i] = 0;
1590 >    recy2_err_spl_[i] = 0;
1591 >    recz2_err_spl_[i] = 0;
1592 >    
1593 >    //pixelVertices
1594 >    nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx  
1595 >    sumptsq_pxlpvtx_[i] = 0;
1596 >    isValid_pxlpvtx_[i] = -1;
1597 >    isFake_pxlpvtx_[i] = -1;
1598 >    recx_pxlpvtx_[i] = 0;
1599 >    recy_pxlpvtx_[i] = 0;
1600 >    recz_pxlpvtx_[i] = 0;
1601 >    recx_err_pxlpvtx_[i] = 0;
1602 >    recy_err_pxlpvtx_[i] = 0;
1603 >    recz_err_pxlpvtx_[i] = 0;
1604 >
1605 >    // matched two-vertices
1606 >    nTrkPV1_spl_twovtx_[i] = 0;
1607 >    nTrkPV2_spl_twovtx_[i] = 0;
1608 >    nTrkPV_twovtx_[i] = 0;  
1609 >    sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1610 >    sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1611 >    deltax_twovtx_[i] = 0;
1612 >    deltay_twovtx_[i] = 0;
1613 >    deltaz_twovtx_[i] = 0;
1614 >    errx_twovtx_[i] = 0;
1615 >    erry_twovtx_[i] = 0;
1616 >    errz_twovtx_[i] = 0;
1617 >    pullx_twovtx_[i] = 0;
1618 >    pully_twovtx_[i] = 0;
1619 >    pullz_twovtx_[i] = 0;
1620 >    
1621 >    //simpv
1622 >    nsimTrkPV_[i] = 0;
1623 >    simx_[i] = 0;
1624 >    simy_[i] = 0;
1625 >    simz_[i] = 0;
1626 >    simptsq_[i] = 0;
1627 >    
1628 >    // residual and pulls with mc truth required
1629 >    deltax_mct_[i] = 0;
1630 >    deltay_mct_[i] = 0;
1631 >    deltaz_mct_[i] = 0;
1632 >    pullx_mct_[i] = 0;
1633 >    pully_mct_[i] = 0;
1634 >    pullz_mct_[i] = 0;
1635 >        
1636 >    deltax_spl1_mct_[i] = 0;
1637 >    deltay_spl1_mct_[i] = 0;
1638 >    deltaz_spl1_mct_[i] = 0;
1639 >    pullx_spl1_mct_[i] = 0;
1640 >    pully_spl1_mct_[i] = 0;
1641 >    pullz_spl1_mct_[i] = 0;
1642 >    
1643 >    deltax_spl2_mct_[i] = 0;
1644 >    deltay_spl2_mct_[i] = 0;
1645 >    deltaz_spl2_mct_[i] = 0;
1646 >    pullx_spl2_mct_[i] = 0;
1647 >    pully_spl2_mct_[i] = 0;
1648 >    pullz_spl2_mct_[i] = 0;
1649 >  }
1650   }
1651  
1652 + double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1653 +  double sum = 0.;
1654 +  double pT;
1655 +  for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1656 +    pT = (**it).pt();
1657 +    sum += pT*pT;
1658 +  }
1659 +  return sum;
1660 + }
1661 +
1662 +
1663  

Diff Legend

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