RecoVertex/src/TrackParameterAnalyzer.cc

0001 #include "Validation/RecoVertex/interface/TrackParameterAnalyzer.h"
0002
0003 //system includes
0004 #include <memory>
0005 #include <vector>
0006
0007 // core framework
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Version/interface/GetReleaseVersion.h"
0012
0013 // Hep MC stuff from CLHEP
0014 #include "CLHEP/Vector/LorentzVector.h"
0015
0016 // track
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018
0019 // Root
0020 #include <TH1.h>
0021 #include <TH2.h>
0022 #include <TFile.h>
0023
0024 //
0025 //
0026 // constants, enums and typedefs
0027 //
0028
0029 //
0030 // static data member definitions
0031 //
0032
0033 //
0034 // constructors and destructor
0035 //
0036 TrackParameterAnalyzer::TrackParameterAnalyzer(const edm::ParameterSet& iConfig)
0037     : edmSimVertexContainerToken_(consumes<edm::SimVertexContainer>(iConfig.getParameter<edm::InputTag>("simG4"))),
0038       edmSimTrackContainerToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simG4"))),
0039       recoTrackCollectionToken_(consumes<reco::TrackCollection>(
0040           edm::InputTag(iConfig.getUntrackedParameter<std::string>("recoTrackProducer")))),
0041       outputFile_(iConfig.getUntrackedParameter<std::string>("outputFile")),
0042       simUnit_(1.0)  //  starting from  CMSSW_1_2_x, I think
0043       ,
0044       verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)) {
0045   //now do whatever initialization is needed
0046   // open output file to store histograms}
0047   auto tversion = edm::getReleaseVersion();
0048   tversion = tversion.erase(tversion.size() - 1, 1).erase(0, 1);
0049   outputFile_ = tversion + "_" + outputFile_;
0050   rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
0051   if ((edm::getReleaseVersion()).find("CMSSW_1_1_", 0) != std::string::npos) {
0052     simUnit_ = 0.1;  // for use in  CMSSW_1_1_1 tutorial
0053   }
0054 }
0055
0056 TrackParameterAnalyzer::~TrackParameterAnalyzer() {
0057   // do anything here that needs to be done at destruction time
0058   // (e.g. close files, deallocate resources etc.)
0059   delete rootFile_;
0060 }
0061
0062 //
0063 // member functions
0064 //
0065 void TrackParameterAnalyzer::beginJob() {
0066   std::cout << " TrackParameterAnalyzer::beginJob  conversion from sim units to rec units is " << simUnit_ << std::endl;
0067
0068   rootFile_->cd();
0069   h1_pull0_ = new TH1F("pull0", "pull q/p", 100, -25., 25.);
0070   h1_pull1_ = new TH1F("pull1", "pull lambda", 100, -25., 25.);
0071   h1_pull2_ = new TH1F("pull2", "pull phi  ", 100, -25., 25.);
0072   h1_pull3_ = new TH1F("pull3", "pull dca  ", 100, -25., 25.);
0073   h1_pull4_ = new TH1F("pull4", "pull zdca ", 100, -25., 25.);
0074
0075   h1_res0_ = new TH1F("res0", "res q/p", 100, -0.1, 0.1);
0076   h1_res1_ = new TH1F("res1", "res lambda", 100, -0.1, 0.1);
0077   h1_res2_ = new TH1F("res2", "res phi  ", 100, -0.1, 0.1);
0078   h1_res3_ = new TH1F("res3", "res dca  ", 100, -0.1, 0.1);
0079   h1_res4_ = new TH1F("res4", "res zdca ", 100, -0.1, 0.1);
0080
0081   h1_Beff_ = new TH1F("Beff", "Beff", 2000, -10., 10.);
0082   h2_dvsphi_ = new TH2F("dvsphi", "dvsphi", 360, -M_PI, M_PI, 100, -0.1, 0.1);
0083   h1_par0_ = new TH1F("par0", "q/p", 100, -0.1, 0.1);
0084   h1_par1_ = new TH1F("par1", "lambda", 100, -M_PI / 2., M_PI / 2.);
0085   h1_par2_ = new TH1F("par2", "phi  ", 100, -M_PI, M_PI);
0086   h1_par3_ = new TH1F("par3", "dca  ", 100, -0.1, 0.1);
0087   h1_par4_ = new TH1F("par4", "zdca ", 1000, -10., 10.);
0088 }
0089
0090 void TrackParameterAnalyzer::endJob() {
0091   rootFile_->cd();
0092   h1_pull0_->Write();
0093   h1_pull1_->Write();
0094   h1_pull2_->Write();
0095   h1_pull3_->Write();
0096   h1_pull4_->Write();
0097
0098   h1_res0_->Write();
0099   h1_res1_->Write();
0100   h1_res2_->Write();
0101   h1_res3_->Write();
0102   h1_res4_->Write();
0103
0104   h1_Beff_->Write();
0105   h2_dvsphi_->Write();
0106   h1_par0_->Write();
0107   h1_par1_->Write();
0108   h1_par2_->Write();
0109   h1_par3_->Write();
0110   h1_par4_->Write();
0111 }
0112
0113 // helper function
0114 bool TrackParameterAnalyzer::match(const ParameterVector& a, const ParameterVector& b) {
0115   double dtheta = a(1) - b(1);
0116   double dphi = a(2) - b(2);
0117   if (dphi > M_PI) {
0118     dphi -= M_2_PI;
0119   } else if (dphi < -M_PI) {
0120     dphi += M_2_PI;
0121   }
0122   return ((fabs(dtheta) < 0.02) && (fabs(dphi) < 0.04));
0123 }
0124
0125 // ------------ method called to produce the data  ------------
0126 void TrackParameterAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0127   using CLHEP::HepLorentzVector;
0128
0129   const double fBfield = 3.8;
0130
0131   edm::Handle<edm::SimVertexContainer> simVtcs;
0132   iEvent.getByToken(edmSimVertexContainerToken_, simVtcs);
0133   if (verbose_) {
0134     std::cout << "SimVertex " << simVtcs->size() << std::endl;
0135     for (edm::SimVertexContainer::const_iterator v = simVtcs->begin(); v != simVtcs->end(); ++v) {
0136       std::cout << "simvtx " << std::setw(10) << std::setprecision(4) << v->position().x() << " " << v->position().y()
0137                 << " " << v->position().z() << " " << v->parentIndex() << " " << v->noParent() << " " << std::endl;
0138     }
0139   }
0140
0141   // get the simulated tracks, extract perigee parameters
0142   edm::Handle<edm::SimTrackContainer> simTrks;
0143   iEvent.getByToken(edmSimTrackContainerToken_, simTrks);
0144
0145   if (verbose_) {
0146     std::cout << "simtrks " << simTrks->size() << std::endl;
0147   }
0148   std::vector<ParameterVector> tsim;
0149   for (edm::SimTrackContainer::const_iterator t = simTrks->begin(); t != simTrks->end(); ++t) {
0150     if (t->noVertex()) {
0151       std::cout << "simtrk  has no vertex" << std::endl;
0152       return;
0153     } else {
0154       // get the vertex position
0155       HepLorentzVector v((*simVtcs)[t->vertIndex()].position().x(),
0156                          (*simVtcs)[t->vertIndex()].position().y(),
0157                          (*simVtcs)[t->vertIndex()].position().z(),
0158                          (*simVtcs)[t->vertIndex()].position().e());
0159       int pdgCode = t->type();
0160
0161       if (pdgCode == -99) {
0162         // such entries cause crashes, no idea what they are
0163         std::cout << "funny particle skipped  , code=" << pdgCode << std::endl;
0164       } else {
0165         double Q = 0;
0166         if ((pdgCode == 11) || (pdgCode == 13) || (pdgCode == 15) || (pdgCode == -211) || (pdgCode == -2212) ||
0167             (pdgCode == 321)) {
0168           Q = -1;
0169         } else if ((pdgCode == -11) || (pdgCode == -13) || (pdgCode == -15) || (pdgCode == 211) || (pdgCode == 2212) ||
0170                    (pdgCode == 321)) {
0171           Q = 1;
0172         } else {
0173           std::cout << pdgCode << " " << std::endl;
0174         }
0175         HepLorentzVector p(t->momentum().x(), t->momentum().y(), t->momentum().z(), t->momentum().e());
0176         if (verbose_) {
0177           std::cout << "simtrk "
0178                     << " gen=" << std::setw(4) << t->genpartIndex() << " vtx=" << std::setw(4) << t->vertIndex()
0179                     << " pdg=" << std::setw(5) << t->type() << " Q=" << std::setw(3) << Q << " pt=" << std::setw(6)
0180                     << p.perp() << " vx=" << std::setw(6) << v.x() << " vy=" << std::setw(6) << v.y()
0181                     << " vz=" << std::setw(6) << v.z() << std::endl;
0182         }
0183         if ((Q != 0) && (p.perp() > 0.1)) {
0184           double x0 = v.x() * simUnit_;
0185           double y0 = v.y() * simUnit_;
0186           double z0 = v.z() * simUnit_;
0187           double kappa = -Q * 0.002998 * fBfield / p.perp();
0188           double D0 = x0 * sin(p.phi()) - y0 * cos(p.phi()) - 0.5 * kappa * (x0 * x0 + y0 * y0);
0189           double q = sqrt(1. - 2. * kappa * D0);
0190           double s0 = (x0 * cos(p.phi()) + y0 * sin(p.phi())) / q;
0191           double s1;
0192           if (fabs(kappa * s0) > 0.001) {
0193             s1 = asin(kappa * s0) / kappa;
0194           } else {
0195             double ks02 = (kappa * s0) * (kappa * s0);
0196             s1 = s0 * (1. + ks02 / 6. + 3. / 40. * ks02 * ks02 + 5. / 112. * pow(ks02, 3));
0197           }
0198           ParameterVector par;
0199           par[reco::TrackBase::i_qoverp] = Q / sqrt(p.perp2() + p.pz() * p.pz());
0200           par[reco::TrackBase::i_lambda] = M_PI / 2. - p.theta();
0201           par[reco::TrackBase::i_phi] = p.phi() - asin(kappa * s0);
0202           par[reco::TrackBase::i_dxy] = 2. * D0 / (1. + q);
0203           par[reco::TrackBase::i_dsz] = z0 * sin(p.theta()) - s1 * cos(p.theta());
0204           tsim.push_back(par);
0205         }
0206       }
0207     }  // has vertex
0208   }  //for loop
0209
0210   // simtrack parameters are in now tsim
0211   // loop over tracks and try to match them to simulated tracks
0212
0213   edm::Handle<reco::TrackCollection> recTracks;
0214   iEvent.getByToken(recoTrackCollectionToken_, recTracks);
0215
0216   for (reco::TrackCollection::const_iterator t = recTracks->begin(); t != recTracks->end(); ++t) {
0217     reco::TrackBase::ParameterVector p = t->parameters();
0218     reco::TrackBase::CovarianceMatrix c = t->covariance();
0219     if (verbose_) {
0220       std::cout << "reco pars= " << p << std::endl;
0221     }
0222     for (std::vector<ParameterVector>::const_iterator s = tsim.begin(); s != tsim.end(); ++s) {
0223       if (match(*s, p)) {
0224         h1_pull0_->Fill((p(0) - (*s)(0)) / sqrt(c(0, 0)));
0225         h1_pull1_->Fill((p(1) - (*s)(1)) / sqrt(c(1, 1)));
0226         h1_pull2_->Fill((p(2) - (*s)(2)) / sqrt(c(2, 2)));
0227         h1_pull3_->Fill((p(3) - (*s)(3)) / sqrt(c(3, 3)));
0228         h1_pull4_->Fill((p(4) - (*s)(4)) / sqrt(c(4, 4)));
0229
0230         h1_res0_->Fill(p(0) - (*s)(0));
0231         h1_res1_->Fill(p(1) - (*s)(1));
0232         h1_res2_->Fill(p(2) - (*s)(2));
0233         h1_res3_->Fill(p(3) - (*s)(3));
0234         h1_res4_->Fill(p(4) - (*s)(4));
0235
0236         h1_Beff_->Fill(p(0) / (*s)(0) * fBfield);
0237         h2_dvsphi_->Fill(p(2), p(3));
0238         h1_par0_->Fill(p(0));
0239         h1_par1_->Fill(p(1));
0240         h1_par2_->Fill(p(2));
0241         h1_par3_->Fill(p(3));
0242         h1_par4_->Fill(p(4));
0243       }
0244     }
0245   }
0246 }