Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:06 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!