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&) override;
47  virtual void endJob() override ;
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, consumesCollector())
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  consumes<reco::SecondaryVertexTagInfoCollection>(svTagInfoProducer_);
114 
115  // Name of the traking pariticle collection
116  trackingTruth_ = config.getUntrackedParameter<edm::InputTag> ( "trackingTruth" );
117  consumes<TrackingVertexCollection>(trackingTruth_);
118 
119  // Number of track categories
120  numberVertexClassifier_ = VertexCategories::Unknown+1;
121 
122  // Define histogram for counting categories
123  TH1Index_["VertexClassifier"] = fs_->make<TH1D>(
124  "VertexClassifier",
125  "Frequency for the different track categories",
126  numberVertexClassifier_,
127  -0.5,
128  numberVertexClassifier_ - 0.5
129  );
130 
131  //--- RecoToSim
132  TH1Index_["rs_All_MatchQuality"]= fs_->make<TH1D>( "rs_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01 );
133  TH1Index_["rs_All_FlightDistance2d"]= fs_->make<TH1D>( "rs_All_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
134  TH1Index_["rs_SecondaryVertex_FlightDistance2d"]= fs_->make<TH1D>( "rs_SecondaryVertex_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
135  TH1Index_["rs_BSV_FlightDistance2d"]= fs_->make<TH1D>( "rs_BSV_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
136  TH1Index_["rs_BWeakDecay_FlightDistance2d"]= fs_->make<TH1D>( "rs_BWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
137  TH1Index_["rs_CWeakDecay_FlightDistance2d"]= fs_->make<TH1D>( "rs_CWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
138  TH1Index_["rs_Light_FlightDistance2d"]= fs_->make<TH1D>( "rs_Light_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5 );
139 
140  TH1Index_["rs_All_nRecVtx"]= fs_->make<TH1D>( "rs_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
141  TH1Index_["rs_SecondaryVertex_nRecVtx"]= fs_->make<TH1D>( "rs_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
142  TH1Index_["rs_BSV_nRecVtx"]= fs_->make<TH1D>( "rs_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
143  TH1Index_["rs_BWeakDecay_nRecVtx"]= fs_->make<TH1D>( "rs_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
144  TH1Index_["rs_CWeakDecay_nRecVtx"]= fs_->make<TH1D>( "rs_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
145  TH1Index_["rs_Light_nRecVtx"]= fs_->make<TH1D>( "rs_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
146 
147 
148  //--- SimToReco
149  TH1Index_["sr_All_MatchQuality"]= fs_->make<TH1D>( "sr_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01);
150  TH1Index_["sr_All_nRecVtx"]= fs_->make<TH1D>( "sr_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
151  TH1Index_["sr_SecondaryVertex_nRecVtx"]= fs_->make<TH1D>( "sr_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
152  TH1Index_["sr_BSV_nRecVtx"]= fs_->make<TH1D>( "sr_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
153  TH1Index_["sr_BWeakDecay_nRecVtx"]= fs_->make<TH1D>( "sr_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
154  TH1Index_["sr_CWeakDecay_nRecVtx"]= fs_->make<TH1D>( "sr_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
155  TH1Index_["sr_Light_nRecVtx"]= fs_->make<TH1D>( "sr_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5 );
156 
157 
158  // Set the proper categories names
159  for (Int_t i = 0; i < numberVertexClassifier_; ++i)
160  TH1Index_["VertexClassifier"]->GetXaxis()->SetBinLabel(i+1, VertexCategories::Names[i]);
161 
162  // book histograms
163  bookRecoToSim("rs_All");
164  bookRecoToSim("rs_SecondaryVertex");
165  bookRecoToSim("rs_BSV");
166  bookRecoToSim("rs_BWeakDecay");
167  bookRecoToSim("rs_CWeakDecay");
168  bookRecoToSim("rs_Light");
169 
170  bookSimToReco("sr_All");
171  bookSimToReco("sr_SecondaryVertex");
172  bookSimToReco("sr_BSV");
173  bookSimToReco("sr_BWeakDecay");
174  bookSimToReco("sr_CWeakDecay");
175  bookSimToReco("sr_Light");
176 }
177 
178 
179 void SVTagInfoValidationAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup)
180 {
181  ++n_event;
182 
183  std::cout << "*** Analyzing " << event.id() << " n_event = " << n_event << std::endl << std::endl;
184 
185  // Set the classifier for a new event
186  classifier_.newEvent(event, setup);
187 
188 
189  // Vertex collection
190  edm::Handle<reco::SecondaryVertexTagInfoCollection> svTagInfoCollection;
191  event.getByLabel(svTagInfoProducer_, svTagInfoCollection);
192 
193  // Get a constant reference to the track history associated to the classifier
194  VertexHistory const & tracer = classifier_.history();
195 
196  cout << "* Event " << n_event << " ; svTagInfoCollection->size() = " << svTagInfoCollection->size() << endl;
197 
198  int rs_nall = 0;
199  int rs_nsv = 0;
200  int rs_nbv = 0;
201  int rs_nbsv = 0;
202  int rs_ncv = 0;
203  int rs_nlv = 0;
204  int nfake = 0;
205 
206  int sr_nall = 0;
207  int sr_nsv = 0;
208  int sr_nbv = 0;
209  int sr_nbsv = 0;
210  int sr_ncv = 0;
211  int sr_nlv = 0;
212  int nmiss = 0;
213 
214  // Loop over the svTagInfo collection.
215  for (std::size_t index = 0; index < svTagInfoCollection->size(); ++index){
216 
217  reco::SecondaryVertexTagInfoRef svTagInfo(svTagInfoCollection, index);
218 
219  // Loop over the vertexes in svTagInfo
220  for ( std::size_t vindex = 0; vindex < svTagInfo->nVertices(); ++vindex ){
221 
222 
223  // Classify the vertices
224  classifier_.evaluate(svTagInfo, vindex);
225 
226  //quality of the match
227  double rs_quality = tracer.quality();
228 
229 
230  // Fill the histogram with the categories
231  for (Int_t i = 0; i != numberVertexClassifier_; ++i) {
232 
233  if ( classifier_.is( (VertexCategories::Category) i ) ) {
234 
235  TH1Index_["VertexClassifier"]->Fill(i);
236  }
237  }
238  if ( !classifier_.is(VertexCategories::Fake) ) {
239 
240  cout << "R2S: MatchQuality = " << rs_quality << endl;
241 
242  TH1Index_["rs_All_MatchQuality"]->Fill( rs_quality );
243  fillRecoToSim("rs_All", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
244  TH1Index_["rs_All_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
245  rs_nall++;
246 
247  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
248 
249  //cout << " R2S VertexCategories::SecondaryVertex" << endl;
250  fillRecoToSim("rs_SecondaryVertex", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
251  TH1Index_["rs_SecondaryVertex_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
252  rs_nsv++;
253  }
254 
255  if ( classifier_.is(VertexCategories::BWeakDecay) ) {
256 
257  //cout << " R2S VertexCategories::BWeakDecay" << endl;
258  fillRecoToSim("rs_BWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
259  TH1Index_["rs_BWeakDecay_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
260  rs_nbv++;
261 
262  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
263  //cout << " R2S VertexCategories::BWeakDecay SecondaryVertex" << endl;
264  fillRecoToSim("rs_BSV", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
265  TH1Index_["rs_BSV_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
266  rs_nbsv++;
267 
268  }
269  }//BWeakDecay
270 
271  else if ( classifier_.is(VertexCategories::CWeakDecay) ) {
272 
273  //cout << " R2S VertexCategories::CWeakDecay" << endl;
274  fillRecoToSim("rs_CWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
275  TH1Index_["rs_CWeakDecay_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
276  rs_ncv++;
277 
278  }
279  else {
280  //cout << " R2S Light (rest of categories)" << endl;
281  fillRecoToSim("rs_Light", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
282  TH1Index_["rs_Light_FlightDistance2d"]->Fill( svTagInfo->flightDistance( vindex, true ).value() );
283  rs_nlv++;
284  }
285  }//end if classifier
286 
287  else {
288  cout << " VertexCategories::Fake!!" << endl;
289  nfake++;
290  }
291 
292  }//end loop over vertices in svTagInfo
293 
294  }//loop over svTagInfo
295 
296  TH1Index_["rs_All_nRecVtx"]->Fill( rs_nall );
297  TH1Index_["rs_SecondaryVertex_nRecVtx"]->Fill( rs_nsv );
298  TH1Index_["rs_BWeakDecay_nRecVtx"]->Fill( rs_nbv );
299  TH1Index_["rs_BSV_nRecVtx"]->Fill( rs_nbsv );
300  TH1Index_["rs_CWeakDecay_nRecVtx"]->Fill( rs_ncv );
301  TH1Index_["rs_Light_nRecVtx"]->Fill( rs_nlv );
302  cout << endl;
303 
304  //----------------------------------------------------------------
305  // SIM TO RECO!
306 
307  // Vertex collection
308  edm::Handle<TrackingVertexCollection> TVCollection;
309  event.getByLabel(trackingTruth_, TVCollection);
310 
311  // Loop over the TV collection.
312  for (std::size_t index = 0; index < TVCollection->size(); ++index){
313 
314  TrackingVertexRef trackingVertex(TVCollection, index);
315 
316  classifier_.evaluate(trackingVertex);
317 
318  double sr_quality = tracer.quality();
319 
320  if ( classifier_.is(VertexCategories::Reconstructed) ) {
321 
322  cout << "S2R: MatchQuality = " << sr_quality << endl;
323 
324  //cout << " S2R VertexCategories::Reconstructed" << endl;
325  TH1Index_["sr_All_MatchQuality"]->Fill( sr_quality );
326  fillSimToReco("sr_All", tracer.recoVertex(), trackingVertex);
327  sr_nall++;
328 
329  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
330 
331  //cout << " S2R VertexCategories::Reconstructed::SecondaryVertex" << endl;
332  fillSimToReco("sr_SecondaryVertex", tracer.recoVertex(), trackingVertex);
333  sr_nsv++;
334  }
335 
336  if ( classifier_.is(VertexCategories::BWeakDecay) ) {
337 
338  //cout << " S2R VertexCategories::Reconstructed::BWeakDecay" << endl;
339  fillSimToReco("sr_BWeakDecay", tracer.recoVertex(), trackingVertex);
340  sr_nbv++;
341 
342  if ( classifier_.is(VertexCategories::SecondaryVertex) ) {
343 
344  //cout << " S2R VertexCategories::Reconstructed::BWeakDecay SecondaryVertex" << endl;
345  fillSimToReco("sr_BSV", tracer.recoVertex(), trackingVertex);
346  sr_nbsv++;
347  }
348 
349  }//BWeakDecay
350 
351  else if ( classifier_.is(VertexCategories::CWeakDecay) ) {
352 
353  //cout << " S2R VertexCategories::CWeakDecay" << endl;
354  fillSimToReco("sr_CWeakDecay", tracer.recoVertex(), trackingVertex);
355  sr_ncv++;
356  }
357 
358  else {
359 
360  //cout << " S2R Light (rest of categories)" << endl;
361  fillSimToReco("sr_Light", tracer.recoVertex(), trackingVertex);
362  sr_nlv++;
363  }
364 
365  }//Reconstructed
366  else {
367  //cout << "##### Not reconstructed!" << endl;
368  nmiss++;
369  }
370 
371  }//TVCollection.size()
372 
373  TH1Index_["sr_All_nRecVtx"]->Fill( sr_nall );
374  TH1Index_["sr_SecondaryVertex_nRecVtx"]->Fill( sr_nsv );
375  TH1Index_["sr_BWeakDecay_nRecVtx"]->Fill( sr_nbv );
376  TH1Index_["sr_BSV_nRecVtx"]->Fill( sr_nbsv );
377  TH1Index_["sr_CWeakDecay_nRecVtx"]->Fill( sr_ncv );
378  TH1Index_["sr_Light_nRecVtx"]->Fill( rs_nlv );
379 
380  rs_total_nall += rs_nall;
381  rs_total_nsv += rs_nsv;
382  rs_total_nbv += rs_nbv;
383  rs_total_nbsv += rs_nbsv;
384  rs_total_ncv += rs_ncv ;
385  rs_total_nlv += rs_nlv;
386  total_nfake += nfake;
387 
388  sr_total_nall += sr_nall;
389  sr_total_nsv += sr_nsv;
390  sr_total_nbv += sr_nbv;
391  sr_total_nbsv += sr_nbsv;
392  sr_total_ncv += sr_ncv;
393  sr_total_nlv += sr_nlv;
394  total_nmiss += nmiss;
395 
396 
397 }
398 
399 
400 void SVTagInfoValidationAnalyzer::bookRecoToSim(std::string const & prefix){
401  // Book pull histograms
402 
403  std::string name = prefix + "_Pullx";
404  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
405  name = prefix + "_Pully";
406  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
407  name = prefix + "_Pullz";
408  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
409 
410  name = prefix + "_Resx";
411  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
412  name = prefix + "_Resy";
413  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
414  name = prefix + "_Resz";
415  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
416 
417  name = prefix + "_Chi2Norm";
418  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0, 10.);
419  name = prefix + "_Chi2Prob";
420  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0., 1.);
421 
422  name = prefix + "_nRecTrks";
423  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
424 
425  name = prefix + "_AverageTrackWeight";
426  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.1, 1.1);
427 
428  name = prefix + "_Mass";
429  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 65, 0., 6.5);
430 
431  name = prefix + "_RecPt";
432  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
433 
434  name = prefix + "_RecEta";
435  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
436 
437  name = prefix + "_RecCharge";
438  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
439 
440  name = prefix + "_RecTrackPt";
441  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
442 
443  name = prefix + "_RecTrackEta";
444  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
445 
446  name = prefix + "_nSimTrks";
447  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
448 
449  name = prefix + "_SimPt";
450  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
451 
452  name = prefix + "_SimEta";
453  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
454 
455  name = prefix + "_SimCharge";
456  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
457 
458  name = prefix + "_SimTrackPt";
459  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 500, 0., 500.);
460 
461  name = prefix + "_SimTrackEta";
462  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
463 
464 
465 }
466 
467 
468 void SVTagInfoValidationAnalyzer::bookSimToReco(std::string const & prefix){
469  // Book pull histograms
470 
471  std::string name = prefix + "_Pullx";
472  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
473  name = prefix + "_Pully";
474  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
475  name = prefix + "_Pullz";
476  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -10., 10.);
477 
478  name = prefix + "_Resx";
479  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
480  name = prefix + "_Resy";
481  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
482  name = prefix + "_Resz";
483  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.05, 0.05);
484 
485  name = prefix + "_Chi2Norm";
486  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0, 10.);
487  name = prefix + "_Chi2Prob";
488  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, 0., 1.);
489 
490  name = prefix + "_nRecTrks";
491  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
492 
493  name = prefix + "_AverageTrackWeight";
494  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 100, -0.1, 1.1);
495 
496  name = prefix + "_Mass";
497  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 65, 0., 6.5);
498 
499  name = prefix + "_RecPt";
500  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
501 
502  name = prefix + "_RecEta";
503  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
504 
505  name = prefix + "_RecCharge";
506  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
507 
508  name = prefix + "_RecTrackPt";
509  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
510 
511  name = prefix + "_RecTrackEta";
512  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
513 
514  name = prefix + "_nSimTrks";
515  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 501, -0.5, 500.5);
516 
517  name = prefix + "_SimPt";
518  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 2000, 0., 1000.);
519 
520  name = prefix + "_SimEta";
521  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
522 
523  name = prefix + "_SimCharge";
524  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 21, -0.5, 20.5);
525 
526  name = prefix + "_SimTrackPt";
527  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 500, 0., 500.);
528 
529  name = prefix + "_SimTrackEta";
530  TH1Index_[name] = fs_->make<TH1D>(name.c_str(), name.c_str(), 200, -3., 3.);
531 
532 
533 }
534 
535 
536 void SVTagInfoValidationAnalyzer::fillRecoToSim(std::string const & prefix, reco::Vertex const & vertex, TrackingVertexRef const & simVertex)
537 {
538 
539  double pullx = (vertex.x() - simVertex->position().x())/vertex.xError();
540  double pully = (vertex.y() - simVertex->position().y())/vertex.yError();
541  double pullz = (vertex.z() - simVertex->position().z())/vertex.zError();
542 
543  double resx = vertex.x() - simVertex->position().x();
544  double resy = vertex.y() - simVertex->position().y();
545  double resz = vertex.z() - simVertex->position().z();
546 
547  double chi2norm = vertex.normalizedChi2();
548  double chi2prob = ChiSquaredProbability( vertex.chi2(), vertex.ndof() );
549 
550  double sum_weight = 0.;
551  double weight = 0.;
552  double tracksize = vertex.tracksSize();
553  math::XYZVector momentum;
554  math::XYZTLorentzVector sum;
555  int charge = 0;
556  double thePiMass = 0.13957;
557  for ( reco::Vertex::trackRef_iterator recDaughter = vertex.tracks_begin() ; recDaughter != vertex.tracks_end(); ++recDaughter) {
558 
559  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
560 
561  vec.SetPx( (**recDaughter).px() );
562  vec.SetPy( (**recDaughter).py() );
563  vec.SetPz( (**recDaughter).pz() );
564  vec.SetM(thePiMass);
565 
566  sum += vec;
567 
568  weight = vertex.trackWeight(*recDaughter);
569  sum_weight += weight;
570 
571  math::XYZVector p;
572  p = (**recDaughter).momentum();
573  momentum += p;
574 
575  charge += (*recDaughter)->charge();
576 
577  TH1Index_[prefix + "_RecTrackPt"]->Fill( (*recDaughter)->pt() );
578  TH1Index_[prefix + "_RecTrackEta"]->Fill( (*recDaughter)->eta() );
579  }//end loop to recDaughters
580  //cout << " average sum of weights = " << sum_weight/tracksize << endl;
581 
582  double mass = sum.M();
583  double pt = sqrt( momentum.Perp2() );
584  double eta = momentum.Eta();
585 
586  math::XYZVector simmomentum;
587  int simcharge = 0;
588  for ( TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin(); simDaughter != simVertex->daughterTracks_end(); ++simDaughter ){
589 
590  math::XYZVector p;
591  p = (**simDaughter).momentum();
592  simmomentum += p;
593 
594  simcharge += (*simDaughter)->charge();
595 
596  TH1Index_[prefix + "_SimTrackPt"]->Fill( (*simDaughter)->pt() );
597  TH1Index_[prefix + "_SimTrackEta"]->Fill( (*simDaughter)->eta() );
598  }
599 
600  double simpt = sqrt( simmomentum.Perp2() );
601  double simeta = simmomentum.Eta();
602 
603  //cout << "[fillRecoToSim] vertex.tracksSize() = " << vertex.tracksSize() << " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() << endl;
604 
605  TH1Index_[prefix + "_nRecTrks"]->Fill( vertex.tracksSize() );
606  TH1Index_[prefix + "_nSimTrks"]->Fill( simVertex->nDaughterTracks() );
607  TH1Index_[prefix + "_Pullx"]->Fill(pullx);
608  TH1Index_[prefix + "_Pully"]->Fill(pully);
609  TH1Index_[prefix + "_Pullz"]->Fill(pullz);
610  TH1Index_[prefix + "_Resx"]->Fill(resx);
611  TH1Index_[prefix + "_Resy"]->Fill(resy);
612  TH1Index_[prefix + "_Resz"]->Fill(resz);
613  TH1Index_[prefix + "_AverageTrackWeight"]->Fill( sum_weight/tracksize );
614  TH1Index_[prefix + "_Chi2Norm"]->Fill(chi2norm);
615  TH1Index_[prefix + "_Chi2Prob"]->Fill(chi2prob);
616  TH1Index_[prefix + "_RecPt"]->Fill(pt);
617  TH1Index_[prefix + "_RecEta"]->Fill(eta);
618  TH1Index_[prefix + "_RecCharge"]->Fill(charge);
619  TH1Index_[prefix + "_Mass"]->Fill(mass);
620  TH1Index_[prefix + "_SimPt"]->Fill(simpt);
621  TH1Index_[prefix + "_SimEta"]->Fill(simeta);
622  TH1Index_[prefix + "_SimCharge"]->Fill(simcharge);
623 }
624 
625 
626 void SVTagInfoValidationAnalyzer::fillSimToReco(std::string const & prefix, reco::VertexBaseRef const & vertex, TrackingVertexRef const & simVertex)
627 {
628 
629  double pullx = (vertex->x() - simVertex->position().x())/vertex->xError();
630  double pully = (vertex->y() - simVertex->position().y())/vertex->yError();
631  double pullz = (vertex->z() - simVertex->position().z())/vertex->zError();
632 
633  double resx = vertex->x() - simVertex->position().x();
634  double resy = vertex->y() - simVertex->position().y();
635  double resz = vertex->z() - simVertex->position().z();
636 
637  double chi2norm = vertex->normalizedChi2();
638  double chi2prob = ChiSquaredProbability(vertex->chi2(), vertex->ndof());
639 
640  double sum_weight = 0.;
641  double weight = 0.;
642  double tracksize = vertex->tracksSize();
643  math::XYZVector momentum;
644  math::XYZTLorentzVector sum;
645  int charge = 0;
646  double thePiMass = 0.13957;
647  for ( reco::Vertex::trackRef_iterator recDaughter = vertex->tracks_begin() ; recDaughter != vertex->tracks_end(); ++recDaughter ) {
648 
649  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > vec;
650 
651  vec.SetPx( (**recDaughter).px() );
652  vec.SetPy( (**recDaughter).py() );
653  vec.SetPz( (**recDaughter).pz() );
654  vec.SetM(thePiMass);
655 
656  sum += vec;
657 
658  weight = vertex->trackWeight(*recDaughter);
659  sum_weight += weight;
660 
661  math::XYZVector p;
662  p = (**recDaughter).momentum();
663  momentum += p;
664 
665  charge += (*recDaughter)->charge();
666 
667  TH1Index_[prefix + "_RecTrackPt"]->Fill( (*recDaughter)->pt() );
668  TH1Index_[prefix + "_RecTrackEta"]->Fill( (*recDaughter)->eta() );
669 
670  }
671  //cout << " average sum of weights = " << sum_weight/tracksize << endl;
672 
673  double mass = sum.M();
674  double pt = sqrt( momentum.Perp2() );
675  double eta = momentum.Eta();
676 
677  math::XYZVector simmomentum;
678  int simcharge = 0;
679  for ( TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin(); simDaughter != simVertex->daughterTracks_end(); ++simDaughter ){
680 
681  math::XYZVector p;
682  p = (**simDaughter).momentum();
683  simmomentum += p;
684 
685  simcharge += (*simDaughter)->charge();
686 
687  TH1Index_[prefix + "_SimTrackPt"]->Fill( (*simDaughter)->pt() );
688  TH1Index_[prefix + "_SimTrackEta"]->Fill( (*simDaughter)->eta() );
689  }
690 
691  double simpt = sqrt( simmomentum.Perp2() );
692  double simeta = simmomentum.Eta();
693 
694  //cout << "[fillSimToReco] vertex->tracksSize() = " << vertex->tracksSize() << " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() << endl;
695 
696  TH1Index_[prefix + "_nRecTrks"]->Fill( vertex->tracksSize() );
697  TH1Index_[prefix + "_nSimTrks"]->Fill( simVertex->nDaughterTracks() );
698  TH1Index_[prefix + "_Pullx"]->Fill(pullx);
699  TH1Index_[prefix + "_Pully"]->Fill(pully);
700  TH1Index_[prefix + "_Pullz"]->Fill(pullz);
701  TH1Index_[prefix + "_Resx"]->Fill(resx);
702  TH1Index_[prefix + "_Resy"]->Fill(resy);
703  TH1Index_[prefix + "_Resz"]->Fill(resz);
704  TH1Index_[prefix + "_AverageTrackWeight"]->Fill( sum_weight/tracksize );
705  TH1Index_[prefix + "_Chi2Norm"]->Fill(chi2norm);
706  TH1Index_[prefix + "_Chi2Prob"]->Fill(chi2prob);
707  TH1Index_[prefix + "_RecPt"]->Fill(pt);
708  TH1Index_[prefix + "_RecEta"]->Fill(eta);
709  TH1Index_[prefix + "_RecCharge"]->Fill(charge);
710  TH1Index_[prefix + "_Mass"]->Fill(mass);
711  TH1Index_[prefix + "_SimPt"]->Fill(simpt);
712  TH1Index_[prefix + "_SimEta"]->Fill(simeta);
713  TH1Index_[prefix + "_SimCharge"]->Fill(simcharge);
714 
715 }
716 
717 
718 void
719 SVTagInfoValidationAnalyzer::endJob() {
720 
721  std::cout << std::endl;
722  std::cout << " ====== Total Number of analyzed events: " << n_event << " ====== " << std::endl;
723  std::cout << " ====== Total Number of R2S All: " << rs_total_nall << " ====== " << std::endl;
724  std::cout << " ====== Total Number of R2S SecondaryVertex: " << rs_total_nsv << " ====== " << std::endl;
725  std::cout << " ====== Total Number of R2S BWeakDecay: " << rs_total_nbv << " ====== " << std::endl;
726  std::cout << " ====== Total Number of R2S BWeakDecay::SecondaryVertex: " << rs_total_nbsv << " ====== " << std::endl;
727  std::cout << " ====== Total Number of R2S CWeakDecay: " << rs_total_ncv << " ====== " << std::endl;
728  std::cout << " ====== Total Number of R2S Light: " << rs_total_nlv << " ====== " << std::endl;
729  std::cout << std::endl;
730  std::cout << " ====== Total Number of S2R All: " << sr_total_nall << " ====== " << std::endl;
731  std::cout << " ====== Total Number of S2R SecondaryVertex: " << sr_total_nsv << " ====== " << std::endl;
732  std::cout << " ====== Total Number of S2R BWeakDecay: " << sr_total_nbv << " ====== " << std::endl;
733  std::cout << " ====== Total Number of S2R BWeakDecay::SecondaryVertex: " << sr_total_nbsv << " ====== " << std::endl;
734  std::cout << " ====== Total Number of S2R CWeakDecay: " << sr_total_ncv << " ====== " << std::endl;
735  std::cout << " ====== Total Number of S2R Light: " << sr_total_nlv << " ====== " << std::endl;
736  std::cout << std::endl;
737  std::cout << " ====== Total Number of Fake Vertices: " << total_nfake << " ====== " << std::endl;
738 
739 }
740 
741 
742 DEFINE_FWK_MODULE(SVTagInfoValidationAnalyzer);
VertexHistory
This class traces the simulated and generated history of a given track.
Definition: VertexHistory.h:18
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:48
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:44
reco::Vertex::zError
double zError() const
error on z
Definition: Vertex.h:118
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
reco::Vertex::y
double y() const
y coordinate
Definition: Vertex.h:110
puppiForMET_cff.weight
tuple weight
Definition: puppiForMET_cff.py:35
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
VertexHistory::recoVertex
const reco::VertexBaseRef & recoVertex() const
Return a reference to the reconstructed track.
Definition: VertexHistory.h:64
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:106
SVTagInfoValidationAnalyzer::fillRecoToSim
void fillRecoToSim(std::string const &, reco::Vertex const &, TrackingVertexRef const &)
Definition: SVTagInfoValidationAnalyzer.cc:536
SVTagInfoValidationAnalyzer::fs_
edm::Service< TFileService > fs_
Definition: SVTagInfoValidationAnalyzer.cc:57
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
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:48
SVTagInfoValidationAnalyzer::bookSimToReco
void bookSimToReco(std::string const &)
Definition: SVTagInfoValidationAnalyzer.cc:468
SVTagInfoValidationAnalyzer::TH1Index_
std::map< std::string, TH1D * > TH1Index_
Definition: SVTagInfoValidationAnalyzer.cc:86
reco::Vertex::chi2
double chi2() const
chi-squares
Definition: Vertex.h:95
reco::Vertex::z
double z() const
y coordinate
Definition: Vertex.h:112
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:102
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:400
reco::Vertex::x
double x() const
x coordinate
Definition: Vertex.h:108
reco::Vertex::xError
double xError() const
error on x
Definition: Vertex.h:114
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
VertexClassifierByProxy::evaluate
VertexClassifierByProxy< Collection > const & evaluate(TrackingVertexRef const &vertex)
Classify the TrackingVertex in categories.
Definition: VertexClassifierByProxy.h:36
SVTagInfoValidationAnalyzer::fillSimToReco
void fillSimToReco(std::string const &, reco::VertexBaseRef const &, TrackingVertexRef const &)
Definition: SVTagInfoValidationAnalyzer.cc:626
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:27
reco::Vertex::trackRef_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
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:104
VertexHistory::quality
double quality() const
Return the quality of the match.
Definition: VertexHistory.h:70
VertexCategories::Names
static const char *const Names[]
Name of the different categories.
Definition: VertexCategories.h:55
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:39
SVTagInfoValidationAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: SVTagInfoValidationAnalyzer.cc:179
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:116
reco::Vertex::tracksSize
size_t tracksSize() const
number of tracks
Definition: Vertex.cc:34
Generated for CMSSW Reference Manual by doxygen 1.8.5