Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:33:26

0001 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
0002 #include "TH1D.h"
0003 #include "TH1F.h"
0004 #include "TH2D.h"
0005 #include "TMath.h"
0006 #include <Math/Functions.h>
0007 #include <Math/SMatrix.h>
0008 #include <Math/SVector.h>
0009 #include <cmath>
0010 #include <map>
0011 #include <string>
0012 
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 #include "FWCore/Framework/interface/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 
0020 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0021 #include "DataFormats/Math/interface/LorentzVector.h"
0022 #include "DataFormats/Math/interface/Vector.h"
0023 #include "Math/VectorUtil.h"
0024 #include "SimTracker/TrackHistory/interface/VertexClassifierByProxy.h"
0025 #include "TROOT.h"
0026 #include <Math/GenVector/PxPyPzE4D.h>
0027 #include <Math/GenVector/PxPyPzM4D.h>
0028 #include <TVector3.h>
0029 //
0030 // class decleration
0031 //
0032 using namespace reco;
0033 using namespace std;
0034 using namespace edm;
0035 
0036 class recoBSVTagInfoValidationAnalyzer : public edm::EDAnalyzer {
0037 public:
0038   typedef dqm::legacy::DQMStore DQMStore;
0039   typedef dqm::legacy::MonitorElement MonitorElement;
0040 
0041   explicit recoBSVTagInfoValidationAnalyzer(const edm::ParameterSet &);
0042 
0043 private:
0044   void analyze(const edm::Event &, const edm::EventSetup &) override;
0045   void endJob() override;
0046   // Member data
0047 
0048   VertexClassifierByProxy<reco::SecondaryVertexTagInfoCollection> classifier_;
0049 
0050   Int_t numberVertexClassifier_;
0051   edm::InputTag trackingTruth_;
0052   edm::InputTag svTagInfoProducer_;
0053 
0054   DQMStore *dqmStore_;
0055   std::string dqmLabel;
0056 
0057   Int_t n_event;
0058   Int_t rs_total_nall;
0059   Int_t rs_total_nsv;
0060   Int_t rs_total_nbv;
0061   Int_t rs_total_nbsv;
0062   Int_t rs_total_ncv;
0063   Int_t rs_total_nlv;
0064   Int_t total_nfake;
0065 
0066   Int_t sr_total_nall;
0067   Int_t sr_total_nsv;
0068   Int_t sr_total_nbv;
0069   Int_t sr_total_nbsv;
0070   Int_t sr_total_ncv;
0071   Int_t sr_total_nlv;
0072   Int_t total_nmiss;
0073 
0074   // Bookeeping of all the histograms per category
0075   void bookRecoToSim(std::string const &);
0076   void bookSimToReco(std::string const &);
0077 
0078   // Fill all histogram per category
0079   void fillRecoToSim(std::string const &, reco::Vertex const &, TrackingVertexRef const &);
0080   void fillSimToReco(std::string const &, reco::VertexBaseRef const &, TrackingVertexRef const &);
0081 
0082   // Histogram handlers
0083   std::map<std::string, MonitorElement *> HistIndex_;
0084 
0085   // consumes
0086   edm::EDGetTokenT<reco::SecondaryVertexTagInfoCollection> svInfoToken;
0087   edm::EDGetTokenT<TrackingVertexCollection> tvToken;
0088 };
0089 
0090 recoBSVTagInfoValidationAnalyzer::recoBSVTagInfoValidationAnalyzer(const edm::ParameterSet &config)
0091     : classifier_(config, consumesCollector()) {
0092   // Initialize counters
0093   n_event = 0;
0094   rs_total_nall = 0;
0095   rs_total_nsv = 0;
0096   rs_total_nbv = 0;
0097   rs_total_nbsv = 0;
0098   rs_total_ncv = 0;
0099   rs_total_nlv = 0;
0100   total_nfake = 0;
0101 
0102   sr_total_nall = 0;
0103   sr_total_nsv = 0;
0104   sr_total_nbv = 0;
0105   sr_total_nbsv = 0;
0106   sr_total_ncv = 0;
0107   sr_total_nlv = 0;
0108   total_nmiss = 0;
0109 
0110   //  get the store
0111   dqmStore_ = edm::Service<DQMStore>().operator->();
0112   dqmLabel = "SVValidation/";
0113   dqmStore_->setCurrentFolder(dqmLabel);
0114 
0115   // Get the track collection
0116   svInfoToken = consumes<reco::SecondaryVertexTagInfoCollection>(config.getParameter<InputTag>("svTagInfoProducer"));
0117   // Name of the traking pariticle collection
0118   tvToken = consumes<TrackingVertexCollection>(config.getParameter<InputTag>("trackingTruth"));
0119   // Number of track categories
0120   numberVertexClassifier_ = VertexCategories::Unknown + 1;
0121 
0122   // Define histogram for counting categories
0123   HistIndex_["VertexClassifier"] = dqmStore_->book1D("VertexClassifier",
0124                                                      "Frequency for the different track categories",
0125                                                      numberVertexClassifier_,
0126                                                      -0.5,
0127                                                      numberVertexClassifier_ - 0.5);
0128 
0129   //--- RecoToSim
0130   HistIndex_["rs_All_MatchQuality"] = dqmStore_->book1D("rs_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01);
0131   HistIndex_["rs_All_FlightDistance2d"] =
0132       dqmStore_->book1D("rs_All_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0133   HistIndex_["rs_SecondaryVertex_FlightDistance2d"] =
0134       dqmStore_->book1D("rs_SecondaryVertex_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0135   HistIndex_["rs_BSV_FlightDistance2d"] =
0136       dqmStore_->book1D("rs_BSV_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0137   HistIndex_["rs_BWeakDecay_FlightDistance2d"] =
0138       dqmStore_->book1D("rs_BWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0139   HistIndex_["rs_CWeakDecay_FlightDistance2d"] =
0140       dqmStore_->book1D("rs_CWeakDecay_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0141   HistIndex_["rs_Light_FlightDistance2d"] =
0142       dqmStore_->book1D("rs_Light_FlightDistance2d", "Transverse flight distance [cm]", 100, 0, 5);
0143 
0144   HistIndex_["rs_All_nRecVtx"] = dqmStore_->book1D("rs_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0145   HistIndex_["rs_SecondaryVertex_nRecVtx"] =
0146       dqmStore_->book1D("rs_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0147   HistIndex_["rs_BSV_nRecVtx"] = dqmStore_->book1D("rs_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0148   HistIndex_["rs_BWeakDecay_nRecVtx"] =
0149       dqmStore_->book1D("rs_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0150   HistIndex_["rs_CWeakDecay_nRecVtx"] =
0151       dqmStore_->book1D("rs_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0152   HistIndex_["rs_Light_nRecVtx"] =
0153       dqmStore_->book1D("rs_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0154 
0155   //--- SimToReco
0156   HistIndex_["sr_All_MatchQuality"] = dqmStore_->book1D("sr_All_MatchQuality", "Quality of Match", 51, -0.01, 1.01);
0157   HistIndex_["sr_All_nRecVtx"] = dqmStore_->book1D("sr_All_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0158   HistIndex_["sr_SecondaryVertex_nRecVtx"] =
0159       dqmStore_->book1D("sr_SecondaryVertex_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0160   HistIndex_["sr_BSV_nRecVtx"] = dqmStore_->book1D("sr_BSV_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0161   HistIndex_["sr_BWeakDecay_nRecVtx"] =
0162       dqmStore_->book1D("sr_BWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0163   HistIndex_["sr_CWeakDecay_nRecVtx"] =
0164       dqmStore_->book1D("sr_CWeakDecay_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0165   HistIndex_["sr_Light_nRecVtx"] =
0166       dqmStore_->book1D("sr_Light_nRecVtx", "Number of Vertices per event", 11, -0.5, 10.5);
0167 
0168   // Set the proper categories names
0169   for (Int_t i = 0; i < numberVertexClassifier_; ++i)
0170     HistIndex_["VertexClassifier"]->setBinLabel(i + 1, VertexCategories::Names[i]);
0171 
0172   // book histograms
0173   bookRecoToSim("rs_All");
0174   bookRecoToSim("rs_SecondaryVertex");
0175   bookRecoToSim("rs_BSV");
0176   bookRecoToSim("rs_BWeakDecay");
0177   bookRecoToSim("rs_CWeakDecay");
0178   bookRecoToSim("rs_Light");
0179 
0180   bookSimToReco("sr_All");
0181   bookSimToReco("sr_SecondaryVertex");
0182   bookSimToReco("sr_BSV");
0183   bookSimToReco("sr_BWeakDecay");
0184   bookSimToReco("sr_CWeakDecay");
0185   bookSimToReco("sr_Light");
0186 }
0187 
0188 void recoBSVTagInfoValidationAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0189   ++n_event;
0190 
0191   std::cout << "*** Analyzing " << event.id() << " n_event = " << n_event << std::endl << std::endl;
0192 
0193   // Set the classifier for a new event
0194   classifier_.newEvent(event, setup);
0195 
0196   // Vertex collection
0197   edm::Handle<reco::SecondaryVertexTagInfoCollection> svTagInfoCollection;
0198   event.getByToken(svInfoToken, svTagInfoCollection);
0199   // Get a constant reference to the track history associated to the classifier
0200   VertexHistory const &tracer = classifier_.history();
0201 
0202   cout << "* Event " << n_event << " ; svTagInfoCollection->size() = " << svTagInfoCollection->size() << endl;
0203 
0204   int rs_nall = 0;
0205   int rs_nsv = 0;
0206   int rs_nbv = 0;
0207   int rs_nbsv = 0;
0208   int rs_ncv = 0;
0209   int rs_nlv = 0;
0210   int nfake = 0;
0211 
0212   int sr_nall = 0;
0213   int sr_nsv = 0;
0214   int sr_nbv = 0;
0215   int sr_nbsv = 0;
0216   int sr_ncv = 0;
0217   int sr_nlv = 0;
0218   int nmiss = 0;
0219 
0220   // Loop over the svTagInfo collection.
0221   for (std::size_t index = 0; index < svTagInfoCollection->size(); ++index) {
0222     reco::SecondaryVertexTagInfoRef svTagInfo(svTagInfoCollection, index);
0223 
0224     // Loop over the vertexes in svTagInfo
0225     for (std::size_t vindex = 0; vindex < svTagInfo->nVertices(); ++vindex) {
0226       // Classify the vertices
0227       classifier_.evaluate(svTagInfo, vindex);
0228 
0229       // quality of the match
0230       double rs_quality = tracer.quality();
0231 
0232       // Fill the histogram with the categories
0233       for (Int_t i = 0; i != numberVertexClassifier_; ++i) {
0234         if (classifier_.is((VertexCategories::Category)i)) {
0235           HistIndex_["VertexClassifier"]->Fill(i);
0236         }
0237       }
0238       if (!classifier_.is(VertexCategories::Fake)) {
0239         HistIndex_["rs_All_MatchQuality"]->Fill(rs_quality);
0240         fillRecoToSim("rs_All", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0241         HistIndex_["rs_All_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0242         rs_nall++;
0243 
0244         if (classifier_.is(VertexCategories::SecondaryVertex)) {
0245           fillRecoToSim("rs_SecondaryVertex", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0246           HistIndex_["rs_SecondaryVertex_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0247           rs_nsv++;
0248         }
0249 
0250         if (classifier_.is(VertexCategories::BWeakDecay)) {
0251           fillRecoToSim("rs_BWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0252           HistIndex_["rs_BWeakDecay_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0253           rs_nbv++;
0254 
0255           if (classifier_.is(VertexCategories::SecondaryVertex)) {
0256             fillRecoToSim("rs_BSV", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0257             HistIndex_["rs_BSV_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0258             rs_nbsv++;
0259           }
0260         }  // BWeakDecay
0261 
0262         else if (classifier_.is(VertexCategories::CWeakDecay)) {
0263           fillRecoToSim("rs_CWeakDecay", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0264           HistIndex_["rs_CWeakDecay_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0265           rs_ncv++;
0266 
0267         } else {
0268           fillRecoToSim("rs_Light", svTagInfo->secondaryVertex(vindex), tracer.simVertex());
0269           HistIndex_["rs_Light_FlightDistance2d"]->Fill(svTagInfo->flightDistance(vindex, true).value());
0270           rs_nlv++;
0271         }
0272       }  // end if classifier
0273 
0274       else {
0275         cout << "    VertexCategories::Fake!!" << endl;
0276         nfake++;
0277       }
0278 
0279     }  // end loop over vertices in svTagInfo
0280 
0281   }  // loop over svTagInfo
0282 
0283   HistIndex_["rs_All_nRecVtx"]->Fill(rs_nall);
0284   HistIndex_["rs_SecondaryVertex_nRecVtx"]->Fill(rs_nsv);
0285   HistIndex_["rs_BWeakDecay_nRecVtx"]->Fill(rs_nbv);
0286   HistIndex_["rs_BSV_nRecVtx"]->Fill(rs_nbsv);
0287   HistIndex_["rs_CWeakDecay_nRecVtx"]->Fill(rs_ncv);
0288   HistIndex_["rs_Light_nRecVtx"]->Fill(rs_nlv);
0289   cout << endl;
0290 
0291   //----------------------------------------------------------------
0292   // SIM TO RECO!
0293 
0294   // Vertex collection
0295   edm::Handle<TrackingVertexCollection> TVCollection;
0296   event.getByToken(tvToken, TVCollection);
0297   // Loop over the TV collection.
0298   for (std::size_t index = 0; index < TVCollection->size(); ++index) {
0299     TrackingVertexRef trackingVertex(TVCollection, index);
0300 
0301     classifier_.evaluate(trackingVertex);
0302 
0303     double sr_quality = tracer.quality();
0304 
0305     if (classifier_.is(VertexCategories::Reconstructed)) {
0306       HistIndex_["sr_All_MatchQuality"]->Fill(sr_quality);
0307       fillSimToReco("sr_All", tracer.recoVertex(), trackingVertex);
0308       sr_nall++;
0309 
0310       if (classifier_.is(VertexCategories::SecondaryVertex)) {
0311         fillSimToReco("sr_SecondaryVertex", tracer.recoVertex(), trackingVertex);
0312         sr_nsv++;
0313       }
0314 
0315       if (classifier_.is(VertexCategories::BWeakDecay)) {
0316         fillSimToReco("sr_BWeakDecay", tracer.recoVertex(), trackingVertex);
0317         sr_nbv++;
0318 
0319         if (classifier_.is(VertexCategories::SecondaryVertex)) {
0320           fillSimToReco("sr_BSV", tracer.recoVertex(), trackingVertex);
0321           sr_nbsv++;
0322         }
0323 
0324       }  // BWeakDecay
0325 
0326       else if (classifier_.is(VertexCategories::CWeakDecay)) {
0327         fillSimToReco("sr_CWeakDecay", tracer.recoVertex(), trackingVertex);
0328         sr_ncv++;
0329       }
0330 
0331       else {
0332         fillSimToReco("sr_Light", tracer.recoVertex(), trackingVertex);
0333         sr_nlv++;
0334       }
0335 
0336     }  // Reconstructed
0337     else {
0338       // cout << "##### Not reconstructed!" << endl;
0339       nmiss++;
0340     }
0341 
0342   }  // TVCollection.size()
0343 
0344   HistIndex_["sr_All_nRecVtx"]->Fill(sr_nall);
0345   HistIndex_["sr_SecondaryVertex_nRecVtx"]->Fill(sr_nsv);
0346   HistIndex_["sr_BWeakDecay_nRecVtx"]->Fill(sr_nbv);
0347   HistIndex_["sr_BSV_nRecVtx"]->Fill(sr_nbsv);
0348   HistIndex_["sr_CWeakDecay_nRecVtx"]->Fill(sr_ncv);
0349   HistIndex_["sr_Light_nRecVtx"]->Fill(rs_nlv);
0350 
0351   rs_total_nall += rs_nall;
0352   rs_total_nsv += rs_nsv;
0353   rs_total_nbv += rs_nbv;
0354   rs_total_nbsv += rs_nbsv;
0355   rs_total_ncv += rs_ncv;
0356   rs_total_nlv += rs_nlv;
0357   total_nfake += nfake;
0358 
0359   sr_total_nall += sr_nall;
0360   sr_total_nsv += sr_nsv;
0361   sr_total_nbv += sr_nbv;
0362   sr_total_nbsv += sr_nbsv;
0363   sr_total_ncv += sr_ncv;
0364   sr_total_nlv += sr_nlv;
0365   total_nmiss += nmiss;
0366 }
0367 
0368 void recoBSVTagInfoValidationAnalyzer::bookRecoToSim(std::string const &prefix) {
0369   // Book pull histograms
0370 
0371   std::string name = prefix + "_Pullx";
0372   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0373   name = prefix + "_Pully";
0374   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0375   name = prefix + "_Pullz";
0376   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0377 
0378   name = prefix + "_Resx";
0379   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0380   name = prefix + "_Resy";
0381   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0382   name = prefix + "_Resz";
0383   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0384 
0385   name = prefix + "_Chi2Norm";
0386   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, 0, 10.);
0387   name = prefix + "_Chi2Prob";
0388   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, 0., 1.);
0389 
0390   name = prefix + "_nRecTrks";
0391   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 501, -0.5, 500.5);
0392 
0393   name = prefix + "_AverageTrackWeight";
0394   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.1, 1.1);
0395 
0396   name = prefix + "_Mass";
0397   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 65, 0., 6.5);
0398 
0399   name = prefix + "_RecPt";
0400   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0401 
0402   name = prefix + "_RecEta";
0403   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0404 
0405   name = prefix + "_RecCharge";
0406   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 21, -0.5, 20.5);
0407 
0408   name = prefix + "_RecTrackPt";
0409   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0410 
0411   name = prefix + "_RecTrackEta";
0412   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0413 
0414   name = prefix + "_nSimTrks";
0415   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 501, -0.5, 500.5);
0416 
0417   name = prefix + "_SimPt";
0418   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0419 
0420   name = prefix + "_SimEta";
0421   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0422 
0423   name = prefix + "_SimCharge";
0424   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 21, -0.5, 20.5);
0425 
0426   name = prefix + "_SimTrackPt";
0427   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 500, 0., 500.);
0428 
0429   name = prefix + "_SimTrackEta";
0430   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0431 }
0432 
0433 void recoBSVTagInfoValidationAnalyzer::bookSimToReco(std::string const &prefix) {
0434   // Book pull histograms
0435 
0436   std::string name = prefix + "_Pullx";
0437   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0438   name = prefix + "_Pully";
0439   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0440   name = prefix + "_Pullz";
0441   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -10., 10.);
0442 
0443   name = prefix + "_Resx";
0444   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0445   name = prefix + "_Resy";
0446   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0447   name = prefix + "_Resz";
0448   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.05, 0.05);
0449 
0450   name = prefix + "_Chi2Norm";
0451   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, 0, 10.);
0452   name = prefix + "_Chi2Prob";
0453   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, 0., 1.);
0454 
0455   name = prefix + "_nRecTrks";
0456   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 501, -0.5, 500.5);
0457 
0458   name = prefix + "_AverageTrackWeight";
0459   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 100, -0.1, 1.1);
0460 
0461   name = prefix + "_Mass";
0462   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 65, 0., 6.5);
0463 
0464   name = prefix + "_RecPt";
0465   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0466 
0467   name = prefix + "_RecEta";
0468   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0469 
0470   name = prefix + "_RecCharge";
0471   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 21, -0.5, 20.5);
0472 
0473   name = prefix + "_RecTrackPt";
0474   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0475 
0476   name = prefix + "_RecTrackEta";
0477   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0478 
0479   name = prefix + "_nSimTrks";
0480   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 501, -0.5, 500.5);
0481 
0482   name = prefix + "_SimPt";
0483   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 2000, 0., 1000.);
0484 
0485   name = prefix + "_SimEta";
0486   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0487 
0488   name = prefix + "_SimCharge";
0489   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 21, -0.5, 20.5);
0490 
0491   name = prefix + "_SimTrackPt";
0492   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 500, 0., 500.);
0493 
0494   name = prefix + "_SimTrackEta";
0495   HistIndex_[name] = dqmStore_->book1D(name.c_str(), name.c_str(), 200, -3., 3.);
0496 }
0497 
0498 void recoBSVTagInfoValidationAnalyzer::fillRecoToSim(std::string const &prefix,
0499                                                      reco::Vertex const &vertex,
0500                                                      TrackingVertexRef const &simVertex) {
0501   double pullx = (vertex.x() - simVertex->position().x()) / vertex.xError();
0502   double pully = (vertex.y() - simVertex->position().y()) / vertex.yError();
0503   double pullz = (vertex.z() - simVertex->position().z()) / vertex.zError();
0504 
0505   double resx = vertex.x() - simVertex->position().x();
0506   double resy = vertex.y() - simVertex->position().y();
0507   double resz = vertex.z() - simVertex->position().z();
0508 
0509   double chi2norm = vertex.normalizedChi2();
0510   double chi2prob = ChiSquaredProbability(vertex.chi2(), vertex.ndof());
0511 
0512   double sum_weight = 0.;
0513   double weight = 0.;
0514   double tracksize = vertex.tracksSize();
0515   math::XYZVector momentum;
0516   math::XYZTLorentzVector sum;
0517   int charge = 0;
0518   double thePiMass = 0.13957;
0519   for (reco::Vertex::trackRef_iterator recDaughter = vertex.tracks_begin(); recDaughter != vertex.tracks_end();
0520        ++recDaughter) {
0521     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> vec;
0522 
0523     vec.SetPx((**recDaughter).px());
0524     vec.SetPy((**recDaughter).py());
0525     vec.SetPz((**recDaughter).pz());
0526     vec.SetM(thePiMass);
0527 
0528     sum += vec;
0529 
0530     weight = vertex.trackWeight(*recDaughter);
0531     sum_weight += weight;
0532 
0533     math::XYZVector p;
0534     p = (**recDaughter).momentum();
0535     momentum += p;
0536 
0537     charge += (*recDaughter)->charge();
0538 
0539     HistIndex_[prefix + "_RecTrackPt"]->Fill((*recDaughter)->pt());
0540     HistIndex_[prefix + "_RecTrackEta"]->Fill((*recDaughter)->eta());
0541   }  // end loop to recDaughters
0542   // cout << "                   average sum of weights = " <<
0543   // sum_weight/tracksize << endl;
0544 
0545   double mass = sum.M();
0546   double pt = sqrt(momentum.Perp2());
0547   double eta = momentum.Eta();
0548 
0549   math::XYZVector simmomentum;
0550   int simcharge = 0;
0551   for (TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin();
0552        simDaughter != simVertex->daughterTracks_end();
0553        ++simDaughter) {
0554     math::XYZVector p;
0555     p = (**simDaughter).momentum();
0556     simmomentum += p;
0557 
0558     simcharge += (*simDaughter)->charge();
0559 
0560     HistIndex_[prefix + "_SimTrackPt"]->Fill((*simDaughter)->pt());
0561     HistIndex_[prefix + "_SimTrackEta"]->Fill((*simDaughter)->eta());
0562   }
0563 
0564   double simpt = sqrt(simmomentum.Perp2());
0565   double simeta = simmomentum.Eta();
0566 
0567   // cout << "[fillRecoToSim]  vertex.tracksSize() = " << vertex.tracksSize() <<
0568   // " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() <<
0569   // endl;
0570 
0571   HistIndex_[prefix + "_nRecTrks"]->Fill(vertex.tracksSize());
0572   HistIndex_[prefix + "_nSimTrks"]->Fill(simVertex->nDaughterTracks());
0573   HistIndex_[prefix + "_Pullx"]->Fill(pullx);
0574   HistIndex_[prefix + "_Pully"]->Fill(pully);
0575   HistIndex_[prefix + "_Pullz"]->Fill(pullz);
0576   HistIndex_[prefix + "_Resx"]->Fill(resx);
0577   HistIndex_[prefix + "_Resy"]->Fill(resy);
0578   HistIndex_[prefix + "_Resz"]->Fill(resz);
0579   HistIndex_[prefix + "_AverageTrackWeight"]->Fill(sum_weight / tracksize);
0580   HistIndex_[prefix + "_Chi2Norm"]->Fill(chi2norm);
0581   HistIndex_[prefix + "_Chi2Prob"]->Fill(chi2prob);
0582   HistIndex_[prefix + "_RecPt"]->Fill(pt);
0583   HistIndex_[prefix + "_RecEta"]->Fill(eta);
0584   HistIndex_[prefix + "_RecCharge"]->Fill(charge);
0585   HistIndex_[prefix + "_Mass"]->Fill(mass);
0586   HistIndex_[prefix + "_SimPt"]->Fill(simpt);
0587   HistIndex_[prefix + "_SimEta"]->Fill(simeta);
0588   HistIndex_[prefix + "_SimCharge"]->Fill(simcharge);
0589 }
0590 
0591 void recoBSVTagInfoValidationAnalyzer::fillSimToReco(std::string const &prefix,
0592                                                      reco::VertexBaseRef const &vertex,
0593                                                      TrackingVertexRef const &simVertex) {
0594   double pullx = (vertex->x() - simVertex->position().x()) / vertex->xError();
0595   double pully = (vertex->y() - simVertex->position().y()) / vertex->yError();
0596   double pullz = (vertex->z() - simVertex->position().z()) / vertex->zError();
0597 
0598   double resx = vertex->x() - simVertex->position().x();
0599   double resy = vertex->y() - simVertex->position().y();
0600   double resz = vertex->z() - simVertex->position().z();
0601 
0602   double chi2norm = vertex->normalizedChi2();
0603   double chi2prob = ChiSquaredProbability(vertex->chi2(), vertex->ndof());
0604 
0605   double sum_weight = 0.;
0606   double weight = 0.;
0607   double tracksize = vertex->tracksSize();
0608   math::XYZVector momentum;
0609   math::XYZTLorentzVector sum;
0610   int charge = 0;
0611   double thePiMass = 0.13957;
0612   for (reco::Vertex::trackRef_iterator recDaughter = vertex->tracks_begin(); recDaughter != vertex->tracks_end();
0613        ++recDaughter) {
0614     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> vec;
0615 
0616     vec.SetPx((**recDaughter).px());
0617     vec.SetPy((**recDaughter).py());
0618     vec.SetPz((**recDaughter).pz());
0619     vec.SetM(thePiMass);
0620 
0621     sum += vec;
0622 
0623     weight = vertex->trackWeight(*recDaughter);
0624     sum_weight += weight;
0625 
0626     math::XYZVector p;
0627     p = (**recDaughter).momentum();
0628     momentum += p;
0629 
0630     charge += (*recDaughter)->charge();
0631 
0632     HistIndex_[prefix + "_RecTrackPt"]->Fill((*recDaughter)->pt());
0633     HistIndex_[prefix + "_RecTrackEta"]->Fill((*recDaughter)->eta());
0634   }
0635   // cout << "                   average sum of weights = " <<
0636   // sum_weight/tracksize << endl;
0637 
0638   double mass = sum.M();
0639   double pt = sqrt(momentum.Perp2());
0640   double eta = momentum.Eta();
0641 
0642   math::XYZVector simmomentum;
0643   int simcharge = 0;
0644   for (TrackingVertex::tp_iterator simDaughter = simVertex->daughterTracks_begin();
0645        simDaughter != simVertex->daughterTracks_end();
0646        ++simDaughter) {
0647     math::XYZVector p;
0648     p = (**simDaughter).momentum();
0649     simmomentum += p;
0650 
0651     simcharge += (*simDaughter)->charge();
0652 
0653     HistIndex_[prefix + "_SimTrackPt"]->Fill((*simDaughter)->pt());
0654     HistIndex_[prefix + "_SimTrackEta"]->Fill((*simDaughter)->eta());
0655   }
0656 
0657   double simpt = sqrt(simmomentum.Perp2());
0658   double simeta = simmomentum.Eta();
0659 
0660   // cout << "[fillSimToReco]  vertex->tracksSize() = " << vertex->tracksSize()
0661   // << " ; simVertex->nDaughterTracks() = " << simVertex->nDaughterTracks() <<
0662   // endl;
0663 
0664   HistIndex_[prefix + "_nRecTrks"]->Fill(vertex->tracksSize());
0665   HistIndex_[prefix + "_nSimTrks"]->Fill(simVertex->nDaughterTracks());
0666   HistIndex_[prefix + "_Pullx"]->Fill(pullx);
0667   HistIndex_[prefix + "_Pully"]->Fill(pully);
0668   HistIndex_[prefix + "_Pullz"]->Fill(pullz);
0669   HistIndex_[prefix + "_Resx"]->Fill(resx);
0670   HistIndex_[prefix + "_Resy"]->Fill(resy);
0671   HistIndex_[prefix + "_Resz"]->Fill(resz);
0672   HistIndex_[prefix + "_AverageTrackWeight"]->Fill(sum_weight / tracksize);
0673   HistIndex_[prefix + "_Chi2Norm"]->Fill(chi2norm);
0674   HistIndex_[prefix + "_Chi2Prob"]->Fill(chi2prob);
0675   HistIndex_[prefix + "_RecPt"]->Fill(pt);
0676   HistIndex_[prefix + "_RecEta"]->Fill(eta);
0677   HistIndex_[prefix + "_RecCharge"]->Fill(charge);
0678   HistIndex_[prefix + "_Mass"]->Fill(mass);
0679   HistIndex_[prefix + "_SimPt"]->Fill(simpt);
0680   HistIndex_[prefix + "_SimEta"]->Fill(simeta);
0681   HistIndex_[prefix + "_SimCharge"]->Fill(simcharge);
0682 }
0683 
0684 void recoBSVTagInfoValidationAnalyzer::endJob() {
0685   std::cout << std::endl;
0686   std::cout << " ====== Total Number of analyzed events: " << n_event << " ====== " << std::endl;
0687   std::cout << " ====== Total Number of R2S All:                         " << rs_total_nall << " ====== " << std::endl;
0688   std::cout << " ====== Total Number of R2S SecondaryVertex:             " << rs_total_nsv << " ====== " << std::endl;
0689   std::cout << " ====== Total Number of R2S BWeakDecay:                  " << rs_total_nbv << " ====== " << std::endl;
0690   std::cout << " ====== Total Number of R2S BWeakDecay::SecondaryVertex: " << rs_total_nbsv << " ====== " << std::endl;
0691   std::cout << " ====== Total Number of R2S CWeakDecay:                  " << rs_total_ncv << " ====== " << std::endl;
0692   std::cout << " ====== Total Number of R2S Light:                       " << rs_total_nlv << " ====== " << std::endl;
0693   std::cout << std::endl;
0694   std::cout << " ====== Total Number of S2R All:                         " << sr_total_nall << " ====== " << std::endl;
0695   std::cout << " ====== Total Number of S2R SecondaryVertex:             " << sr_total_nsv << " ====== " << std::endl;
0696   std::cout << " ====== Total Number of S2R BWeakDecay:                  " << sr_total_nbv << " ====== " << std::endl;
0697   std::cout << " ====== Total Number of S2R BWeakDecay::SecondaryVertex: " << sr_total_nbsv << " ====== " << std::endl;
0698   std::cout << " ====== Total Number of S2R CWeakDecay:                  " << sr_total_ncv << " ====== " << std::endl;
0699   std::cout << " ====== Total Number of S2R Light:                       " << sr_total_nlv << " ====== " << std::endl;
0700   std::cout << std::endl;
0701   std::cout << " ====== Total Number of Fake Vertices:              " << total_nfake << " ====== " << std::endl;
0702 }
0703 
0704 DEFINE_FWK_MODULE(recoBSVTagInfoValidationAnalyzer);