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.10 by yygao, Sun Nov 8 02:56:57 2009 UTC vs.
Revision 1.21 by yygao, Sun Feb 28 00:07:09 2010 UTC

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

Diff Legend

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