CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SVTagInfoValidationAnalyzer.cc
Go to the documentation of this file.
1 #include <map>
2 #include <string>
3 #include <cmath>
4 #include <Math/Functions.h>
5 #include <Math/SVector.h>
6 #include <Math/SMatrix.h>
7 #include "TH1F.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TMath.h"
11 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
12 
13 #include "FWCore/Framework/interface/Frameworkfwd.h"
14 #include "FWCore/Framework/interface/EDAnalyzer.h"
15 #include "FWCore/Framework/interface/MakerMacros.h"
16 #include "FWCore/ParameterSet/interface/ParameterSet.h"
17 #include "FWCore/ServiceRegistry/interface/Service.h"
18 
19 #include "CommonTools/UtilAlgos/interface/TFileService.h"
20 
21 #include "SimTracker/TrackHistory/interface/VertexClassifierByProxy.h"
22 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
23 #include "DataFormats/Math/interface/Vector.h"
24 #include "TROOT.h"
25 #include "Math/VectorUtil.h"
26 #include <TVector3.h>
27 #include <Math/GenVector/PxPyPzE4D.h>
28 #include <Math/GenVector/PxPyPzM4D.h>
29 #include "DataFormats/Math/interface/LorentzVector.h"
30 //
31 // class decleration
32 //
33 using namespace reco;
34 using namespace std;
35 using namespace edm;
36 
37 class SVTagInfoValidationAnalyzer : public edm::EDAnalyzer
38 {
39 
40 public:
41 
42  explicit SVTagInfoValidationAnalyzer(const edm::ParameterSet&);
43 
44 private:
45 
46  virtual void analyze(const edm::Event&, const edm::EventSetup&);
47  virtual void endJob() ;
48  // Member data
49 
50  VertexClassifierByProxy<reco::SecondaryVertexTagInfoCollection> classifier_;
51 
52  Int_t numberVertexClassifier_;
53  edm::InputTag trackingTruth_;
54  edm::InputTag svTagInfoProducer_;
55 
56  // Get the file service
57  edm::Service<TFileService> fs_;
58 
59  Int_t n_event;
60  Int_t rs_total_nall ;
61  Int_t rs_total_nsv ;
62  Int_t rs_total_nbv ;
63  Int_t rs_total_nbsv ;
64  Int_t rs_total_ncv ;
65  Int_t rs_total_nlv ;
66  Int_t total_nfake ;
67 
68  Int_t sr_total_nall ;
69  Int_t sr_total_nsv ;
70  Int_t sr_total_nbv ;
71  Int_t sr_total_nbsv ;
72  Int_t sr_total_ncv ;
73  Int_t sr_total_nlv ;
74  Int_t total_nmiss ;
75 
76 
77  // Bookeeping of all the histograms per category
78  void bookRecoToSim(std::string const &);
79  void bookSimToReco(std::string const &);
80 
81  // Fill all histogram per category
82  void fillRecoToSim(std::string const &, reco::Vertex const &, TrackingVertexRef const &);
83  void fillSimToReco(std::string const &, reco::VertexBaseRef const &, TrackingVertexRef const &);
84 
85  // Histogram handlers
86  std::map<std::string, TH1D *> TH1Index_;
87 
88 };
89 
90 
91 SVTagInfoValidationAnalyzer::SVTagInfoValidationAnalyzer(const edm::ParameterSet& config) : classifier_(config)
92 {
93  //Initialize counters
94  n_event = 0;
95  rs_total_nall = 0;
96  rs_total_nsv = 0;
97  rs_total_nbv = 0;
98  rs_total_nbsv = 0;
99  rs_total_ncv = 0;
100  rs_total_nlv = 0;
101  total_nfake = 0;
102 
103  sr_total_nall = 0;
104  sr_total_nsv = 0;
105  sr_total_nbv = 0;
106  sr_total_nbsv = 0;
107  sr_total_ncv = 0;
108  sr_total_nlv = 0;
109  total_nmiss = 0;
110 
111  // Get the track collection
112  svTagInfoProducer_ = config.getUntrackedParameter<edm::InputTag> ( "svTagInfoProducer" );
113 
114  // Name of the traking pariticle collection
115  trackingTruth_ = config.getUntrackedParameter<edm::InputTag> ( "trackingTruth" );
116 
117  // Number of track categories
118  numberVertexClassifier_ = VertexCategories::Unknown+1;
119 
120  // Define histogram for counting categories
121  TH1Index_["VertexClassifier"] = fs_->make<TH1D>(
122  "VertexClassifier",
123  "Frequency for the different track categories",
124  numberVertexClassifier_,
125  -0.5,
126  numberVertexClassifier_ - 0.5
127  );
128 
129  //--- RecoToSim
130  TH1Index_["rs_All_MatchQuality"]= fs_->make<TH1D>( "rs_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01 );
131  TH1Index_["rs_All_FlightDistance2d"]= fs_->make<TH1D>( "rs_All_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
132  TH1Index_["rs_SecondaryVertex_FlightDistance2d"]= fs_->make<TH1D>( "rs_SecondaryVertex_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
133  TH1Index_["rs_BSV_FlightDistance2d"]= fs_->make<TH1D>( "rs_BSV_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
134  TH1Index_["rs_BWeakDecay_FlightDistance2d"]= fs_->make<TH1D>( "rs_BWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
135  TH1Index_["rs_CWeakDecay_FlightDistance2d"]= fs_->make<TH1D>( "rs_CWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
136  TH1Index_["rs_Light_FlightDistance2d"]= fs_->make<TH1D>( "rs_Light_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
137 
138  TH1Index_["rs_All_nRecVtx"]= fs_->make<TH1D>( "rs_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
139  TH1Index_["rs_SecondaryVertex_nRecVtx"]= fs_->make<TH1D>( "rs_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
140  TH1Index_["rs_BSV_nRecVtx"]= fs_->make<TH1D>( "rs_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
141  TH1Index_["rs_BWeakDecay_nRecVtx"]= fs_->make<TH1D>( "rs_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
142  TH1Index_["rs_CWeakDecay_nRecVtx"]= fs_->make<TH1D>( "rs_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
143  TH1Index_["rs_Light_nRecVtx"]= fs_->make<TH1D>( "rs_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
144 
145 
146  //--- SimToReco
147  TH1Index_["sr_All_MatchQuality"]= fs_->make<TH1D>( "sr_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01);
148  TH1Index_["sr_All_nRecVtx"]= fs_->make<TH1D>( "sr_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
149  TH1Index_["sr_SecondaryVertex_nRecVtx"]= fs_->make<TH1D>( "sr_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
150  TH1Index_["sr_BSV_nRecVtx"]= fs_->make<TH1D>( "sr_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
151  TH1Index_["sr_BWeakDecay_nRecVtx"]= fs_->make<TH1D>( "sr_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
152  TH1Index_["sr_CWeakDecay_nRecVtx"]= fs_->make<TH1D>( "sr_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
153  TH1Index_["sr_Light_nRecVtx"]= fs_->make<TH1D>( "sr_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
154 
155 
156  // Set the proper categories names
157  for (Int_t i = 0; i < numberVertexClassifier_; ++i)
158  TH1Index_["VertexClassifier"]->GetXaxis()->SetBinLabel(i+1, VertexCategories::Names[i]);
159 
160  // book histograms
161  bookRecoToSim("rs_All");
162  bookRecoToSim("rs_SecondaryVertex");
163  bookRecoToSim("rs_BSV");
164  bookRecoToSim("rs_BWeakDecay");
165  bookRecoToSim("rs_CWeakDecay");
166  bookRecoToSim("rs_Light");
167 
168  bookSimToReco("sr_All");
169  bookSimToReco("sr_SecondaryVertex");
170  bookSimToReco("sr_BSV");
171  bookSimToReco("sr_BWeakDecay");
172  bookSimToReco("sr_CWeakDecay");
173  bookSimToReco("sr_Light");
174 }
175 
176 
177 void SVTagInfoValidationAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup)
178 {
179  ++n_event;
180 
181  std::cout << "*** Analyzing " << event.id() << " n_event = " << n_event << std::endl << std::endl;
182 
183  // Set the classifier for a new event
184  classifier_.newEvent(event, setup);
185 
186 
187  // Vertex collection
188  edm::Handle<reco::SecondaryVertexTagInfoCollection> svTagInfoCollection;
189  event.getByLabel(svTagInfoProducer_, svTagInfoCollection);
190 
191  // Get a constant reference to the track history associated to the classifier
192  VertexHistory const & tracer = classifier_.history();
193 
194  cout << "* Event " << n_event << " ; svTagInfoCollection->size() = " << svTagInfoCollection->size() << endl;
195 
196  int rs_nall = 0;
197  int rs_nsv = 0;
198  int rs_nbv = 0;
199  int rs_nbsv = 0;
200  int rs_ncv = 0;
201  int rs_nlv = 0;
202  int nfake = 0;
203 
204  int sr_nall = 0;
205  int sr_nsv = 0;
206  int sr_nbv = 0;
207  int sr_nbsv = 0;
208  int sr_ncv = 0;
209  int sr_nlv = 0;
210  int nmiss = 0;
211 
212  // Loop over the svTagInfo collection.
213  for (std::size_t index = 0; index < svTagInfoCollection->size(); ++index){
214 
215  reco::SecondaryVertexTagInfoRef svTagInfo(svTagInfoCollection, index);
216 
217  // Loop over the vertexes in svTagInfo
218  for ( std::size_t vindex = 0; vindex < svTagInfo->nVertices(); ++vindex ){
219 
220 
221  // Classify the vertices
222  classifier_.evaluate(svTagInfo, vindex);
223 
224  //quality of the match
225  double rs_quality = tracer.quality();
226 
227 
228  // Fill the histogram with the categories
229  for (Int_t i = 0; i != numberVertexClassifier_; ++i) {
230 
231  if ( classifier_.is( (VertexCategories::Category) i ) ) {
232 
233  TH1Index_["VertexClassifier"]->Fill(i);
234  }
235  }
236  if ( !classifier_.is(VertexCategories::Fake) ) {
237 
238  cout << "R2S: MatchQuality = " << rs_quality << endl;
239 
240  TH1Index_["rs_All_MatchQuality"]->Fill( rs_quality );
241  fillRecoToSim("rs_All", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
242  TH1Index_["rs_All_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
243  rs_nall++;
244 
245  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
246 
247  //cout << " R2S VertexCategories::SecondaryVertex" << endl;
248  fillRecoToSim("rs_SecondaryVertex", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
249  TH1Index_["rs_SecondaryVertex_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
250  rs_nsv++;
251  }
252 
253  if ( classifier_.is(VertexCategories::BWeakDecay) ) {
254 
255  //cout << " R2S VertexCategories::BWeakDecay" << endl;
256  fillRecoToSim("rs_BWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
257  TH1Index_["rs_BWeakDecay_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
258  rs_nbv++;
259 
260  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
261  //cout << " R2S VertexCategories::BWeakDecay SecondaryVertex" << endl;
262  fillRecoToSim("rs_BSV", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
263  TH1Index_["rs_BSV_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
264  rs_nbsv++;
265 
266  }
267  }//BWeakDecay
268 
269  else if ( classifier_.is(VertexCategories::CWeakDecay) ) {
270 
271  //cout << " R2S VertexCategories::CWeakDecay" << endl;
272  fillRecoToSim("rs_CWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
273  TH1Index_["rs_CWeakDecay_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
274  rs_ncv++;
275 
276  }
277  else {
278  //cout << " R2S Light (rest of categories)" << endl;
279  fillRecoToSim("rs_Light", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
280  TH1Index_["rs_Light_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
281  rs_nlv++;
282  }
283  }//end if classifier
284 
285  else {
286  cout << " VertexCategories::Fake!!" << endl;
287  nfake++;
288  }
289 
290  }//end loop over vertices in svTagInfo
291 
292  }//loop over svTagInfo
293 
294  TH1Index_["rs_All_nRecVtx"]->Fill( rs_nall );
295  TH1Index_["rs_SecondaryVertex_nRecVtx"]->Fill( rs_nsv );
296  TH1Index_["rs_BWeakDecay_nRecVtx"]->Fill( rs_nbv );
297  TH1Index_["rs_BSV_nRecVtx"]->Fill( rs_nbsv );
298  TH1Index_["rs_CWeakDecay_nRecVtx"]->Fill( rs_ncv );
299  TH1Index_["rs_Light_nRecVtx"]->Fill( rs_nlv );
300  cout << endl;
301 
302  //----------------------------------------------------------------
303  // SIM TO RECO!
304 
305  // Vertex collection
306  edm::Handle<TrackingVertexCollection> TVCollection;
307  event.getByLabel(trackingTruth_, TVCollection);
308 
309  // Loop over the TV collection.
310  for (std::size_t index = 0; index < TVCollection->size(); ++index){
311 
312  TrackingVertexRef trackingVertex(TVCollection, index);
313 
314  classifier_.evaluate(trackingVertex);
315 
316  double sr_quality = tracer.quality();
317 
318  if ( classifier_.is(VertexCategories::Reconstructed) ) {
319 
320  cout << "S2R: MatchQuality = " << sr_quality << endl;
321 
322  //cout << " S2R VertexCategories::Reconstructed" << endl;
323  TH1Index_["sr_All_MatchQuality"]->Fill( sr_quality );
324  fillSimToReco("sr_All", tracer.recoVertex(), trackingVertex);
325  sr_nall++;
326 
327  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
328 
329  //cout << " S2R VertexCategories::Reconstructed::SecondaryVertex" << endl;
330  fillSimToReco("sr_SecondaryVertex", tracer.recoVertex(), trackingVertex);
331  sr_nsv++;
332  }
333 
334  if ( classifier_.is(VertexCategories::BWeakDecay) ) {
335 
336  //cout << " S2R VertexCategories::Reconstructed::BWeakDecay" << endl;
337  fillSimToReco("sr_BWeakDecay", tracer.recoVertex(), trackingVertex);
338  sr_nbv++;
339 
340  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
341 
342  //cout << " S2R VertexCategories::Reconstructed::BWeakDecay SecondaryVertex" << endl;
343  fillSimToReco("sr_BSV", tracer.recoVertex(), trackingVertex);
344  sr_nbsv++;
345  }
346 
347  }//BWeakDecay
348 
349  else if ( classifier_.is(VertexCategories::CWeakDecay) ) {
350 
351  //cout << " S2R VertexCategories::CWeakDecay" << endl;
352  fillSimToReco("sr_CWeakDecay", tracer.recoVertex(), trackingVertex);
353  sr_ncv++;
354  }
355 
356  else {
357 
358  //cout << " S2R Light (rest of categories)" << endl;
359  fillSimToReco("sr_Light", tracer.recoVertex(), trackingVertex);
360  sr_nlv++;
361  }
362 
363  }//Reconstructed
364  else {
365  //cout << "##### Not reconstructed!" << endl;
366  nmiss++;
367  }
368 
369  }//TVCollection.size()
370 
371  TH1Index_["sr_All_nRecVtx"]->Fill( sr_nall );
372  TH1Index_["sr_SecondaryVertex_nRecVtx"]->Fill( sr_nsv );
373  TH1Index_["sr_BWeakDecay_nRecVtx"]->Fill( sr_nbv );
374  TH1Index_["sr_BSV_nRecVtx"]->Fill( sr_nbsv );
375  TH1Index_["sr_CWeakDecay_nRecVtx"]->Fill( sr_ncv );
376  TH1Index_["sr_Light_nRecVtx"]->Fill( rs_nlv );
377 
378  rs_total_nall += rs_nall;
379  rs_total_nsv += rs_nsv;
380  rs_total_nbv += rs_nbv;
381  rs_total_nbsv += rs_nbsv;
382  rs_total_ncv += rs_ncv ;
383  rs_total_nlv += rs_nlv;
384  total_nfake += nfake;
385 
386  sr_total_nall += sr_nall;
387  sr_total_nsv += sr_nsv;
388  sr_total_nbv += sr_nbv;
389  sr_total_nbsv += sr_nbsv;
390  sr_total_ncv += sr_ncv;
391  sr_total_nlv += sr_nlv;
392  total_nmiss += nmiss;
393 
394 
395 }
396 
397 
398 void SVTagInfoValidationAnalyzer::bookRecoToSim(std::string const & prefix){
399  // Book pull histograms
400 
401  std::string name = prefix + "_Pullx";
402  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
403  name = prefix + "_Pully";
404  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
405  name = prefix + "_Pullz";
406  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
407 
408  name = prefix + "_Resx";
409  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
410  name = prefix + "_Resy";
411  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
412  name = prefix + "_Resz";
413  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
414 
415  name = prefix + "_Chi2Norm";
416  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0, 10.);
417  name = prefix + "_Chi2Prob";
418  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0., 1.);
419 
420  name = prefix + "_nRecTrks";
421  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
422 
423  name = prefix + "_AverageTrackWeight";
424  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.1, 1.1);
425 
426  name = prefix + "_Mass";
427  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 65, 0., 6.5);
428 
429  name = prefix + "_RecPt";
430  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
431 
432  name = prefix + "_RecEta";
433  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
434 
435  name = prefix + "_RecCharge";
436  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
437 
438  name = prefix + "_RecTrackPt";
439  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
440 
441  name = prefix + "_RecTrackEta";
442  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
443 
444  name = prefix + "_nSimTrks";
445  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
446 
447  name = prefix + "_SimPt";
448  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
449 
450  name = prefix + "_SimEta";
451  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
452 
453  name = prefix + "_SimCharge";
454  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
455 
456  name = prefix + "_SimTrackPt";
457  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 500, 0., 500.);
458 
459  name = prefix + "_SimTrackEta";
460  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
461 
462 
463 }
464 
465 
466 void SVTagInfoValidationAnalyzer::bookSimToReco(std::string const & prefix){
467  // Book pull histograms
468 
469  std::string name = prefix + "_Pullx";
470  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
471  name = prefix + "_Pully";
472  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
473  name = prefix + "_Pullz";
474  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
475 
476  name = prefix + "_Resx";
477  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
478  name = prefix + "_Resy";
479  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
480  name = prefix + "_Resz";
481  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
482 
483  name = prefix + "_Chi2Norm";
484  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0, 10.);
485  name = prefix + "_Chi2Prob";
486  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0., 1.);
487 
488  name = prefix + "_nRecTrks";
489  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
490 
491  name = prefix + "_AverageTrackWeight";
492  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.1, 1.1);
493 
494  name = prefix + "_Mass";
495  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 65, 0., 6.5);
496 
497  name = prefix + "_RecPt";
498  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
499 
500  name = prefix + "_RecEta";
501  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
502 
503  name = prefix + "_RecCharge";
504  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
505 
506  name = prefix + "_RecTrackPt";
507  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
508 
509  name = prefix + "_RecTrackEta";
510  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
511 
512  name = prefix + "_nSimTrks";
513  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
514 
515  name = prefix + "_SimPt";
516  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
517 
518  name = prefix + "_SimEta";
519  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
520 
521  name = prefix + "_SimCharge";
522  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
523 
524  name = prefix + "_SimTrackPt";
525  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 500, 0., 500.);
526 
527  name = prefix + "_SimTrackEta";
528  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
529 
530 
531 }
532 
533 
534 void SVTagInfoValidationAnalyzer::fillRecoToSim(std::string const & prefix, reco::Vertex const & vertex, TrackingVertexRef const & simVertex)
535 {
536 
537  double pullx = (vertex.x() - simVertex->position().x())/vertex.xError();
538  double pully = (vertex.y() - simVertex->position().y())/vertex.yError();
539  double pullz = (vertex.z() - simVertex->position().z())/vertex.zError();
540 
541  double resx = vertex.x() - simVertex->position().x();
542  double resy = vertex.y() - simVertex->position().y();
543  double resz = vertex.z() - simVertex->position().z();
544 
545  double chi2norm = vertex.normalizedChi2();
546  double chi2prob = ChiSquaredProbability( vertex.chi2(), vertex.ndof() );
547 
548  double sum_weight = 0.;
549  double weight = 0.;
550  double tracksize = vertex.tracksSize();
551  math::XYZVector momentum;
552  math::XYZTLorentzVector sum;
553  int charge = 0;
554  double thePiMass = 0.13957;
555  for ( reco::Vertex::trackRef_iterator recDaughter = vertex.tracks_begin() ; recDaughter != vertex.tracks_end(); ++recDaughter) {
556 
557  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
558 
559  vec.SetPx( (**recDaughter).px() );
560  vec.SetPy( (**recDaughter).py() );
561  vec.SetPz( (**recDaughter).pz() );
562  vec.SetM(thePiMass);
563 
564  sum += vec;
565 
566  weight = vertex.trackWeight(*recDaughter);
567  sum_weight += weight;
568 
569  math::XYZVector p;
570  p = (**recDaughter).momentum();
571  momentum += p;
572 
573  charge += (*recDaughter)->charge();
574 
575  TH1Index_[prefix + "_RecTrackPt"]->Fill( (*recDaughter)->pt() );
576  TH1Index_[prefix + "_RecTrackEta"]->Fill( (*recDaughter)->eta() );
577  }//end loop to recDaughters
578  //cout << " average sum of weights = " << sum_weight/tracksize << endl;
579 
580  double mass = sum.M();
581  double pt = sqrt( momentum.Perp2() );
582  double eta = momentum.Eta();
583 
584  math::XYZVector simmomentum;
585  int simcharge = 0;
586  for ( TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin(); simDaughter != simVertex->daughterTracks_end(); ++simDaughter ){
587 
588  math::XYZVector p;
589  p = (**simDaughter).momentum();
590  simmomentum += p;
591 
592  simcharge += (*simDaughter)->charge();
593 
594  TH1Index_[prefix + "_SimTrackPt"]->Fill( (*simDaughter)->pt() );
595  TH1Index_[prefix + "_SimTrackEta"]->Fill( (*simDaughter)->eta() );
596  }
597 
598  double simpt = sqrt( simmomentum.Perp2() );
599  double simeta = simmomentum.Eta();
600 
601  //cout << "[fillRecoToSim] vertex.tracksSize() = " << vertex.tracksSize() << " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() << endl;
602 
603  TH1Index_[prefix + "_nRecTrks"]->Fill( vertex.tracksSize() );
604  TH1Index_[prefix + "_nSimTrks"]->Fill( simVertex->nDaughterTracks() );
605  TH1Index_[prefix + "_Pullx"]->Fill(pullx);
606  TH1Index_[prefix + "_Pully"]->Fill(pully);
607  TH1Index_[prefix + "_Pullz"]->Fill(pullz);
608  TH1Index_[prefix + "_Resx"]->Fill(resx);
609  TH1Index_[prefix + "_Resy"]->Fill(resy);
610  TH1Index_[prefix + "_Resz"]->Fill(resz);
611  TH1Index_[prefix + "_AverageTrackWeight"]->Fill( sum_weight/tracksize );
612  TH1Index_[prefix + "_Chi2Norm"]->Fill(chi2norm);
613  TH1Index_[prefix + "_Chi2Prob"]->Fill(chi2prob);
614  TH1Index_[prefix + "_RecPt"]->Fill(pt);
615  TH1Index_[prefix + "_RecEta"]->Fill(eta);
616  TH1Index_[prefix + "_RecCharge"]->Fill(charge);
617  TH1Index_[prefix + "_Mass"]->Fill(mass);
618  TH1Index_[prefix + "_SimPt"]->Fill(simpt);
619  TH1Index_[prefix + "_SimEta"]->Fill(simeta);
620  TH1Index_[prefix + "_SimCharge"]->Fill(simcharge);
621 }
622 
623 
624 void SVTagInfoValidationAnalyzer::fillSimToReco(std::string const & prefix, reco::VertexBaseRef const & vertex, TrackingVertexRef const & simVertex)
625 {
626 
627  double pullx = (vertex->x() - simVertex->position().x())/vertex->xError();
628  double pully = (vertex->y() - simVertex->position().y())/vertex->yError();
629  double pullz = (vertex->z() - simVertex->position().z())/vertex->zError();
630 
631  double resx = vertex->x() - simVertex->position().x();
632  double resy = vertex->y() - simVertex->position().y();
633  double resz = vertex->z() - simVertex->position().z();
634 
635  double chi2norm = vertex->normalizedChi2();
636  double chi2prob = ChiSquaredProbability(vertex->chi2(), vertex->ndof());
637 
638  double sum_weight = 0.;
639  double weight = 0.;
640  double tracksize = vertex->tracksSize();
641  math::XYZVector momentum;
642  math::XYZTLorentzVector sum;
643  int charge = 0;
644  double thePiMass = 0.13957;
645  for ( reco::Vertex::trackRef_iterator recDaughter = vertex->tracks_begin() ; recDaughter != vertex->tracks_end(); ++recDaughter ) {
646 
647  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
648 
649  vec.SetPx( (**recDaughter).px() );
650  vec.SetPy( (**recDaughter).py() );
651  vec.SetPz( (**recDaughter).pz() );
652  vec.SetM(thePiMass);
653 
654  sum += vec;
655 
656  weight = vertex->trackWeight(*recDaughter);
657  sum_weight += weight;
658 
659  math::XYZVector p;
660  p = (**recDaughter).momentum();
661  momentum += p;
662 
663  charge += (*recDaughter)->charge();
664 
665  TH1Index_[prefix + "_RecTrackPt"]->Fill( (*recDaughter)->pt() );
666  TH1Index_[prefix + "_RecTrackEta"]->Fill( (*recDaughter)->eta() );
667 
668  }
669  //cout << " average sum of weights = " << sum_weight/tracksize << endl;
670 
671  double mass = sum.M();
672  double pt = sqrt( momentum.Perp2() );
673  double eta = momentum.Eta();
674 
675  math::XYZVector simmomentum;
676  int simcharge = 0;
677  for ( TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin(); simDaughter != simVertex->daughterTracks_end(); ++simDaughter ){
678 
679  math::XYZVector p;
680  p = (**simDaughter).momentum();
681  simmomentum += p;
682 
683  simcharge += (*simDaughter)->charge();
684 
685  TH1Index_[prefix + "_SimTrackPt"]->Fill( (*simDaughter)->pt() );
686  TH1Index_[prefix + "_SimTrackEta"]->Fill( (*simDaughter)->eta() );
687  }
688 
689  double simpt = sqrt( simmomentum.Perp2() );
690  double simeta = simmomentum.Eta();
691 
692  //cout << "[fillSimToReco] vertex->tracksSize() = " << vertex->tracksSize() << " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() << endl;
693 
694  TH1Index_[prefix + "_nRecTrks"]->Fill( vertex->tracksSize() );
695  TH1Index_[prefix + "_nSimTrks"]->Fill( simVertex->nDaughterTracks() );
696  TH1Index_[prefix + "_Pullx"]->Fill(pullx);
697  TH1Index_[prefix + "_Pully"]->Fill(pully);
698  TH1Index_[prefix + "_Pullz"]->Fill(pullz);
699  TH1Index_[prefix + "_Resx"]->Fill(resx);
700  TH1Index_[prefix + "_Resy"]->Fill(resy);
701  TH1Index_[prefix + "_Resz"]->Fill(resz);
702  TH1Index_[prefix + "_AverageTrackWeight"]->Fill( sum_weight/tracksize );
703  TH1Index_[prefix + "_Chi2Norm"]->Fill(chi2norm);
704  TH1Index_[prefix + "_Chi2Prob"]->Fill(chi2prob);
705  TH1Index_[prefix + "_RecPt"]->Fill(pt);
706  TH1Index_[prefix + "_RecEta"]->Fill(eta);
707  TH1Index_[prefix + "_RecCharge"]->Fill(charge);
708  TH1Index_[prefix + "_Mass"]->Fill(mass);
709  TH1Index_[prefix + "_SimPt"]->Fill(simpt);
710  TH1Index_[prefix + "_SimEta"]->Fill(simeta);
711  TH1Index_[prefix + "_SimCharge"]->Fill(simcharge);
712 
713 }
714 
715 
716 void
717 SVTagInfoValidationAnalyzer::endJob() {
718 
719  std::cout << std::endl;
720  std::cout << " ====== Total Number of analyzed events: " << n_event << " ====== " << std::endl;
721  std::cout << " ====== Total Number of R2S All: " << rs_total_nall << " ====== " << std::endl;
722  std::cout << " ====== Total Number of R2S SecondaryVertex: " << rs_total_nsv << " ====== " << std::endl;
723  std::cout << " ====== Total Number of R2S BWeakDecay: " << rs_total_nbv << " ====== " << std::endl;
724  std::cout << " ====== Total Number of R2S BWeakDecay::SecondaryVertex: " << rs_total_nbsv << " ====== " << std::endl;
725  std::cout << " ====== Total Number of R2S CWeakDecay: " << rs_total_ncv << " ====== " << std::endl;
726  std::cout << " ====== Total Number of R2S Light: " << rs_total_nlv << " ====== " << std::endl;
727  std::cout << std::endl;
728  std::cout << " ====== Total Number of S2R All: " << sr_total_nall << " ====== " << std::endl;
729  std::cout << " ====== Total Number of S2R SecondaryVertex: " << sr_total_nsv << " ====== " << std::endl;
730  std::cout << " ====== Total Number of S2R BWeakDecay: " << sr_total_nbv << " ====== " << std::endl;
731  std::cout << " ====== Total Number of S2R BWeakDecay::SecondaryVertex: " << sr_total_nbsv << " ====== " << std::endl;
732  std::cout << " ====== Total Number of S2R CWeakDecay: " << sr_total_ncv << " ====== " << std::endl;
733  std::cout << " ====== Total Number of S2R Light: " << sr_total_nlv << " ====== " << std::endl;
734  std::cout << std::endl;
735  std::cout << " ====== Total Number of Fake Vertices: " << total_nfake << " ====== " << std::endl;
736 
737 }
738 
739 
740 DEFINE_FWK_MODULE(SVTagInfoValidationAnalyzer);
VertexHistory
This class traces the simulated and generated history of a given track.
Definition: VertexHistory.h:17
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
VertexClassifier::history
VertexHistory const & history() const
Returns a reference to the vertex history used in the classification.
Definition: VertexClassifier.h:45
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:45
reco::Vertex::zError
double zError() const
error on z
Definition: Vertex.h:105
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
reco::Vertex::y
double y() const
y coordinate
Definition: Vertex.h:97
VertexHistory::recoVertex
const reco::VertexBaseRef & recoVertex() const
Return a reference to the reconstructed track.
Definition: VertexHistory.h:62
eta
T eta() const
Definition: Basic3DVectorLD.h:172
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
SVTagInfoValidationAnalyzer::fillRecoToSim
void fillRecoToSim(std::string const &, reco::Vertex const &, TrackingVertexRef const &)
Definition: SVTagInfoValidationAnalyzer.cc:534
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
SVTagInfoValidationAnalyzer::fs_
edm::Service< TFileService > fs_
Definition: SVTagInfoValidationAnalyzer.cc:57
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
SVTagInfoValidationAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: SVTagInfoValidationAnalyzer.cc:177
reco::Vertex::trackWeight
float trackWeight(const TrackBaseRef &r) const
returns the weight with which a Track has contributed to the vertex-fit.
VertexCategories::Category
Category
Categories available to vertexes.
Definition: VertexCategories.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
SVTagInfoValidationAnalyzer::bookSimToReco
void bookSimToReco(std::string const &)
Definition: SVTagInfoValidationAnalyzer.cc:466
SVTagInfoValidationAnalyzer::TH1Index_
std::map< std::string, TH1D * > TH1Index_
Definition: SVTagInfoValidationAnalyzer.cc:86
reco::Vertex::chi2
double chi2() const
chi-squares
Definition: Vertex.h:82
reco::Vertex::z
double z() const
y coordinate
Definition: Vertex.h:99
ChiSquaredProbability
float ChiSquaredProbability(double chiSquared, double nrDOF)
Definition: ChiSquaredProbability.cc:13
VertexCategories::is
bool is(Category category) const
Returns track flag for a given category.
Definition: VertexCategories.h:67
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:89
HistoryBase::simVertex
const TrackingVertexRef & simVertex() const
Return the initial tracking vertex from the history.
Definition: HistoryBase.h:83
SVTagInfoValidationAnalyzer::bookRecoToSim
void bookRecoToSim(std::string const &)
Definition: SVTagInfoValidationAnalyzer.cc:398
reco::Vertex::x
double x() const
x coordinate
Definition: Vertex.h:95
reco::Vertex::xError
double xError() const
error on x
Definition: Vertex.h:101
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
VertexClassifierByProxy::evaluate
VertexClassifierByProxy< Collection > const & evaluate(TrackingVertexRef const &vertex)
Classify the TrackingVertex in categories.
Definition: VertexClassifierByProxy.h:33
SVTagInfoValidationAnalyzer::fillSimToReco
void fillSimToReco(std::string const &, reco::VertexBaseRef const &, TrackingVertexRef const &)
Definition: SVTagInfoValidationAnalyzer.cc:624
SVTagInfoValidationAnalyzer::classifier_
VertexClassifierByProxy< reco::SecondaryVertexTagInfoCollection > classifier_
Definition: SVTagInfoValidationAnalyzer.cc:50
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
VertexClassifierByProxy::newEvent
virtual void newEvent(edm::Event const &event, edm::EventSetup const &config)
Pre-process event information (for accessing reconstraction information).
Definition: VertexClassifierByProxy.h:24
scaleCards.mass
tuple mass
Definition: scaleCards.py:27
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:38
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:55
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
reco::Vertex::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f.
Definition: Vertex.h:91
VertexCategories::Names
static const char * Names[]
Name of the different categories.
Definition: VertexCategories.h:55
VertexHistory::quality
double quality() const
Return the quality of the match.
Definition: VertexHistory.h:68
SVTagInfoValidationAnalyzer::SVTagInfoValidationAnalyzer
SVTagInfoValidationAnalyzer(const edm::ParameterSet &)
Definition: SVTagInfoValidationAnalyzer.cc:91
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:40
HcalObjRepresent::setup
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Definition: HcalObjRepresent.h:628
reco::Vertex::yError
double yError() const
error on y
Definition: Vertex.h:103
reco::Vertex::tracksSize
size_t tracksSize() const
number of tracks
Definition: Vertex.cc:35
Generated for CMSSW Reference Manual by doxygen 1.8.5