Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:01 

0001 #include <memory>
0002 
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0011 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0012 
0013 #include "MagneticField/Engine/interface/MagneticField.h"
0014 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0015 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0016 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0017 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0018 
0019 #include <iostream>
0020 #include <string>
0021 
0022 #include <TH1.h>
0023 #include <TH2.h>
0024 #include <TROOT.h>
0025 #include <TFile.h>
0026 #include <TCanvas.h>
0027 
0028 using namespace edm;
0029 using namespace std;
0030 using namespace reco;
0031 
0032 class TrackValidator : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0033 public:
0034   TrackValidator(const edm::ParameterSet& pset)
0035       : sim(pset.getParameter<string>("sim")),
0036         label(pset.getParameter<vector<string> >("label")),
0037         out(pset.getParameter<string>("out")),
0038         open(pset.getParameter<string>("open")),
0039         min(pset.getParameter<double>("min")),
0040         max(pset.getParameter<double>("max")),
0041         nint(pset.getParameter<int>("nint")),
0042         partId(pset.getParameter<int>("partId")),
0043         theMFToken_(esConsumes()),
0044         theBToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) {
0045     hFile = new TFile(out.c_str(), open.c_str());
0046   }
0047 
0048   ~TrackValidator() {
0049     if (hFile != 0) {
0050       hFile->Close();
0051       delete hFile;
0052     }
0053   }
0054 
0055   void beginRun(edm::Run const& run, const edm::EventSetup& setup) override {
0056     for (unsigned int j = 0; j < label.size(); j++) {
0057       vector<double> etaintervalsv;
0058       vector<int> totSIMv, totRECv;
0059       vector<TH1F*> ptdistribv;
0060       vector<TH1F*> etadistribv;
0061 
0062       double step = (max - min) / nint;
0063       ostringstream title, name;
0064       etaintervalsv.push_back(0);
0065       for (double d = min; d < max; d = d + step) {
0066         etaintervalsv.push_back(d + step);
0067         totSIMv.push_back(0);
0068         totRECv.push_back(0);
0069         name.str("");
0070         title.str("");
0071         name << "pt[" << d << "," << d + step << "]";
0072         title << "p_{t} residue " << d << "<#eta<" << d + step;
0073         ptdistribv.push_back(new TH1F(name.str().c_str(), title.str().c_str(), 200, -2, 2));
0074         name.str("");
0075         title.str("");
0076         name << "eta[" << d << "," << d + step << "]";
0077         title << "eta residue " << d << "<#eta<" << d + step;
0078         etadistribv.push_back(new TH1F(name.str().c_str(), title.str().c_str(), 200, -0.2, 0.2));
0079       }
0080       etaintervals.push_back(etaintervalsv);
0081       totSIM.push_back(totSIMv);
0082       totREC.push_back(totRECv);
0083       ptdistrib.push_back(ptdistribv);
0084       etadistrib.push_back(etadistribv);
0085 
0086       h_ptSIM.push_back(new TH1F("ptSIM", "generated p_{t}", 5500, 0, 110));
0087       h_etaSIM.push_back(new TH1F("etaSIM", "generated pseudorapidity", 500, 0, 5));
0088       h_tracksSIM.push_back(new TH1F("tracksSIM", "number of simluated tracks", 100, -0.5, 99.5));
0089       h_vertposSIM.push_back(new TH1F("vertposSIM", "Transverse position of sim vertices", 1000, -0.5, 10000.5));
0090 
0091       //     h_pt     = new TH1F("pt", "p_{t} residue", 2000, -500, 500 );
0092       h_pt.push_back(new TH1F("pullPt", "pull of p_{t}", 100, -10, 10));
0093       h_pt2.push_back(new TH1F("pt2", "p_{t} residue (#tracks>1)", 300, -15, 15));
0094       h_eta.push_back(new TH1F("eta", "pseudorapidity residue", 1000, -0.1, 0.1));
0095       h_tracks.push_back(new TH1F("tracks", "number of reconstructed tracks", 10, -0.5, 9.5));
0096       h_nchi2.push_back(new TH1F("nchi2", "normalized chi2", 200, 0, 20));
0097       h_hits.push_back(new TH1F("hits", "number of hits per track", 30, -0.5, 29.5));
0098       h_effic.push_back(new TH1F("effic", "efficiency vs #eta", nint, &etaintervals[j][0]));
0099       h_ptrmsh.push_back(new TH1F("PtRMS", "PtRMS vs #eta", nint, &etaintervals[j][0]));
0100       h_deltaeta.push_back(new TH1F("etaRMS", "etaRMS vs #eta", nint, &etaintervals[j][0]));
0101       h_charge.push_back(new TH1F("charge", "charge", 3, -1.5, 1.5));
0102 
0103       h_pullTheta.push_back(new TH1F("pullTheta", "pull of theta parameter", 100, -10, 10));
0104       h_pullPhi0.push_back(new TH1F("pullPhi0", "pull of phi0 parameter", 100, -10, 10));
0105       h_pullD0.push_back(new TH1F("pullD0", "pull of d0 parameter", 100, -10, 10));
0106       h_pullDz.push_back(new TH1F("pullDz", "pull of dz parameter", 100, -10, 10));
0107       h_pullK.push_back(new TH1F("pullK", "pull of k parameter", 100, -10, 10));
0108 
0109       chi2_vs_nhits.push_back(new TH2F("chi2_vs_nhits", "chi2 vs nhits", 25, 0, 25, 100, 0, 10));
0110       chi2_vs_eta.push_back(new TH2F("chi2_vs_eta", "chi2 vs eta", nint, min, max, 100, 0, 10));
0111       nhits_vs_eta.push_back(new TH2F("nhits_vs_eta", "nhits vs eta", nint, min, max, 25, 0, 25));
0112       ptres_vs_eta.push_back(new TH2F("ptres_vs_eta", "ptresidue vs eta", nint, min, max, 200, -2, 2));
0113       etares_vs_eta.push_back(new TH2F("etares_vs_eta", "etaresidue vs eta", nint, min, max, 200, -0.1, 0.1));
0114     }
0115   }
0116 
0117   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override {
0118     std::cout << "In TrackValidator\n";
0119     edm::ESHandle<MagneticField> theMF = setup.getHandle(theMFToken_);
0120     edm::ESHandle<TransientTrackBuilder> theB = setup.getHandle(theBToken_);
0121     for (unsigned int w = 0; w < label.size(); w++) {
0122       //
0123       //get collections from the event
0124       //
0125       edm::Handle<SimTrackContainer> simTrackCollection;
0126       event.getByLabel(sim, simTrackCollection);
0127       const SimTrackContainer simTC = *(simTrackCollection.product());
0128 
0129       edm::Handle<SimVertexContainer> simVertexCollection;
0130       event.getByLabel(sim, simVertexCollection);
0131       const SimVertexContainer simVC = *(simVertexCollection.product());
0132 
0133       edm::Handle<reco::TrackCollection> trackCollection;
0134       event.getByLabel(label[w], trackCollection);
0135       const reco::TrackCollection tC = *(trackCollection.product());
0136 
0137       vector<TransientTrack> t_tks = (*theB).build(trackCollection);
0138       cout << "Found: " << t_tks.size() << " reconstructed tracks"
0139            << "\n";
0140 
0141       //
0142       //fill simulation histograms
0143       //
0144       int st = 0;
0145       for (SimTrackContainer::const_iterator simTrack = simTC.begin(); simTrack != simTC.end(); simTrack++) {
0146         if (abs(simTrack->momentum().eta()) > max || abs(simTrack->momentum().eta()) < min)
0147           continue;
0148         st++;
0149         h_ptSIM[w]->Fill(simTrack->momentum().pt());
0150         h_etaSIM[w]->Fill(simTrack->momentum().eta());
0151 
0152         h_vertposSIM[w]->Fill(simVC[simTrack->vertIndex()].position().pt());
0153 
0154         if (simTrack->type() != partId)
0155           continue;
0156         //compute number of tracks per eta interval
0157         int i = 0;
0158         for (vector<double>::iterator h = etaintervals[w].begin(); h != etaintervals[w].end() - 1; h++) {
0159           if (abs(simTrack->momentum().eta()) > etaintervals[w][i] &&
0160               abs(simTrack->momentum().eta()) < etaintervals[w][i + 1]) {
0161             totSIM[w][i]++;
0162             bool doit = false;
0163             for (reco::TrackCollection::const_iterator track = tC.begin(); track != tC.end(); track++) {
0164               if (abs(track->pt() - simTrack->momentum().pt()) < (simTrack->momentum().pt() * 0.1))
0165                 doit = true;
0166             }
0167             if (doit)
0168               totREC[w][i]++;
0169           }
0170           i++;
0171         }
0172       }
0173       if (st != 0)
0174         h_tracksSIM[w]->Fill(st);
0175 
0176       //       for (SimVertexContainer::const_iterator simVertex=simVC.begin();simVertex!=simVC.end();simVertex++){
0177       //    h_vertposSIM[w]->Fill(simVertex->position().perp());
0178       //       }
0179 
0180       //
0181       //fill reconstructed track histograms
0182       //
0183       int rt = 0;
0184       for (vector<TransientTrack>::const_iterator track = t_tks.begin(); track != t_tks.end(); track++) {
0185         TrajectoryStateClosestToPoint tscp = track->impactPointTSCP();
0186         cout << tscp.perigeeParameters().vector() << tscp.perigeeError().covarianceMatrix() << endl;
0187 
0188         if (abs(track->track().eta()) > max || abs(track->track().eta()) < min)
0189           continue;
0190 
0191         rt++;
0192 
0193         //nchi2 and hits global distributions
0194         h_nchi2[w]->Fill(track->normalizedChi2());
0195         h_hits[w]->Fill(track->numberOfValidHits());
0196         chi2_vs_nhits[w]->Fill(track->numberOfValidHits(), track->normalizedChi2());
0197         chi2_vs_eta[w]->Fill(track->track().eta(), track->normalizedChi2());
0198         nhits_vs_eta[w]->Fill(track->track().eta(), track->numberOfValidHits());
0199         h_charge[w]->Fill(track->charge());
0200 
0201         //pt, eta residue, theta, phi0, d0, dz pull
0202         double ptres = 1000;
0203         double etares = 1000;
0204         double thetares = 1000;
0205         double phi0res = 1000;
0206         double d0res = 1000;
0207         double dzres = 1000;
0208         double kres = 1000;
0209         for (SimTrackContainer::const_iterator simTrack = simTC.begin(); simTrack != simTC.end(); simTrack++) {
0210           if (simTrack->type() != partId)
0211             continue;
0212           double tmp = track->track().pt() - simTrack->momentum().pt();
0213           if (tC.size() > 1)
0214             h_pt2[w]->Fill(tmp);
0215           if (abs(tmp) < abs(ptres)) {
0216             ptres = tmp;
0217 
0218             etares = track->initialFreeState().momentum().eta() - simTrack->momentum().eta();
0219             thetares =
0220                 (tscp.perigeeParameters().theta() - simTrack->momentum().theta()) / tscp.perigeeError().thetaError();
0221             phi0res = (tscp.perigeeParameters().phi() - simTrack->momentum().phi()) / tscp.perigeeError().phiError();
0222             d0res =
0223                 (tscp.perigeeParameters().transverseImpactParameter() - simVC[simTrack->vertIndex()].position().pt()) /
0224                 tscp.perigeeError().transverseImpactParameterError();
0225             dzres =
0226                 (tscp.perigeeParameters().longitudinalImpactParameter() - simVC[simTrack->vertIndex()].position().z()) /
0227                 tscp.perigeeError().longitudinalImpactParameterError();
0228 
0229             const math::XYZTLorentzVectorD& vertexPosition = simVC[simTrack->vertIndex()].position();
0230             GlobalVector magField =
0231                 theMF->inTesla(GlobalPoint(vertexPosition.x(), vertexPosition.y(), vertexPosition.z()));
0232             kres = (tscp.perigeeParameters().transverseCurvature() -
0233                     (-track->charge() * 2.99792458e-3 * magField.z() / simTrack->momentum().pt())) /
0234                    tscp.perigeeError().transverseCurvatureError();
0235 
0236             //      cout << "track->d0(): " << track->d0() << endl;
0237             //      cout << "simVC[simTrack->vertIndex()].position().perp(): " << simVC[simTrack->vertIndex()].position().perp() << endl;
0238             //      cout << "track->dz(): " << track->dz() << endl;
0239             //      cout << "simVC[simTrack->vertIndex()].position().z(): " << simVC[simTrack->vertIndex()].position().z() << endl;
0240             //      cout << "track->transverseCurvature(): " << track->transverseCurvature() << endl;
0241 
0242             //      cout << "-track->charge()*2.99792458e-3 * 4./simTrack->momentum().perp(): " << -track->charge()*2.99792458e-3 * 4./simTrack->momentum().perp() << endl;
0243             //      cout << "-track->charge()*2.99792458e-3 * 4./simTrack->momentum().perp(): " << -track->charge()*2.99792458e-3 * magField.z()/simTrack->momentum().perp() << endl;
0244           }
0245         }
0246         cout << etares << endl;
0247         h_pt[w]->Fill(
0248             ptres / (tscp.perigeeError().transverseCurvatureError() / tscp.perigeeParameters().transverseCurvature()));
0249         h_eta[w]->Fill(etares);
0250         ptres_vs_eta[w]->Fill(track->track().eta(), ptres);
0251         etares_vs_eta[w]->Fill(track->track().eta(), etares);
0252         h_pullTheta[w]->Fill(thetares);
0253         h_pullPhi0[w]->Fill(phi0res);
0254         h_pullD0[w]->Fill(d0res);
0255         h_pullDz[w]->Fill(dzres);
0256         h_pullK[w]->Fill(kres);
0257 
0258         //pt residue distribution per eta interval
0259         int i = 0;
0260         for (vector<TH1F*>::iterator h = ptdistrib[w].begin(); h != ptdistrib[w].end(); h++) {
0261           for (SimTrackContainer::const_iterator simTrack = simTC.begin(); simTrack != simTC.end(); simTrack++) {
0262             if (simTrack->type() != partId)
0263               continue;
0264             ptres = 1000;
0265             if (abs(simTrack->momentum().eta()) > etaintervals[w][i] &&
0266                 abs(simTrack->momentum().eta()) < etaintervals[w][i + 1]) {
0267               double tmp = track->track().pt() - simTrack->momentum().pt();
0268               if (abs(tmp) < abs(ptres))
0269                 ptres = tmp;
0270             }
0271           }
0272           (*h)->Fill(ptres);
0273           i++;
0274         }
0275         //eta residue distribution per eta interval
0276         i = 0;
0277         for (vector<TH1F*>::iterator h = etadistrib[w].begin(); h != etadistrib[w].end(); h++) {
0278           for (SimTrackContainer::const_iterator simTrack = simTC.begin(); simTrack != simTC.end(); simTrack++) {
0279             if (simTrack->type() != partId)
0280               continue;
0281             etares = 1000;
0282             ptres = 1000;
0283             if (abs(simTrack->momentum().eta()) > etaintervals[w][i] &&
0284                 abs(simTrack->momentum().eta()) < etaintervals[w][i + 1]) {
0285               double tmp = track->track().pt() - simTrack->momentum().pt();
0286               if (abs(tmp) < abs(ptres))
0287                 etares = track->track().eta() - simTrack->momentum().eta();
0288             }
0289           }
0290           (*h)->Fill(etares);
0291           i++;
0292         }
0293       }
0294       if (rt != 0)
0295         h_tracks[w]->Fill(rt);
0296     }
0297   }
0298 
0299   void endJob() override {
0300     for (unsigned int w = 0; w < label.size(); w++) {
0301       TDirectory* p = hFile->mkdir(label[w].c_str());
0302 
0303       //write simulation histos
0304       TDirectory* simD = p->mkdir("simulation");
0305       simD->cd();
0306       h_ptSIM[w]->Write();
0307       h_etaSIM[w]->Write();
0308       h_tracksSIM[w]->Write();
0309       h_vertposSIM[w]->Write();
0310 
0311       //fill pt rms plot versus eta and write pt residue distribution per eta interval histo
0312       TDirectory* ptD = p->mkdir("ptdistribution");
0313       ptD->cd();
0314       int i = 0;
0315       for (vector<TH1F*>::iterator h = ptdistrib[w].begin(); h != ptdistrib[w].end(); h++) {
0316         (*h)->Write();
0317         h_ptrmsh[w]->Fill(etaintervals[w][i + 1] - 0.00001, (*h)->GetRMS());
0318         i++;
0319       }
0320 
0321       //fill eta rms plot versus eta and write eta residue distribution per eta interval histo
0322       TDirectory* etaD = p->mkdir("etadistribution");
0323       etaD->cd();
0324       i = 0;
0325       for (vector<TH1F*>::iterator h = etadistrib[w].begin(); h != etadistrib[w].end(); h++) {
0326         (*h)->Write();
0327         h_deltaeta[w]->Fill(etaintervals[w][i + 1] - 0.00001, (*h)->GetRMS());
0328         i++;
0329       }
0330 
0331       //write the other histos
0332       p->cd();
0333       int j = 0;
0334       for (vector<int>::iterator h = totSIM[w].begin(); h != totSIM[w].end(); h++) {
0335         if (totSIM[w][j])
0336           h_effic[w]->Fill(etaintervals[w][j + 1] - 0.00001, ((double)totREC[w][j]) / ((double)totSIM[w][j]));
0337         else
0338           h_effic[w]->Fill(etaintervals[w][j + 1] - 0.00001, 0);
0339         j++;
0340       }
0341 
0342       h_pt[w]->Write();
0343       h_pt2[w]->Write();
0344       h_eta[w]->Write();
0345       h_tracks[w]->Write();
0346       h_nchi2[w]->Write();
0347       h_hits[w]->Write();
0348       h_effic[w]->Write();
0349       h_ptrmsh[w]->Write();
0350       h_deltaeta[w]->Write();
0351       chi2_vs_nhits[w]->Write();
0352       chi2_vs_eta[w]->Write();
0353       nhits_vs_eta[w]->Write();
0354       ptres_vs_eta[w]->Write();
0355       etares_vs_eta[w]->Write();
0356       h_charge[w]->Write();
0357 
0358       h_pullTheta[w]->Write();
0359       h_pullPhi0[w]->Write();
0360       h_pullD0[w]->Write();
0361       h_pullDz[w]->Write();
0362       h_pullK[w]->Write();
0363     }
0364 
0365     hFile->Close();
0366   }
0367 
0368   void endRun(edm::Run const& run, const edm::EventSetup& setup) override {}
0369 
0370 private:
0371   string sim;
0372   vector<string> label;
0373   string out, open;
0374   double min, max;
0375   int nint, partId;
0376 
0377   vector<TH1F*> h_ptSIM, h_etaSIM, h_tracksSIM, h_vertposSIM;
0378   vector<TH1F*> h_tracks, h_nchi2, h_hits, h_effic, h_ptrmsh, h_deltaeta, h_charge;
0379   vector<TH1F*> h_pt, h_eta, h_pullTheta, h_pullPhi0, h_pullD0, h_pullDz, h_pullK, h_pt2;
0380   vector<TH2F*> chi2_vs_nhits, chi2_vs_eta, nhits_vs_eta, ptres_vs_eta, etares_vs_eta;
0381 
0382   vector<vector<double> > etaintervals;
0383   vector<vector<int> > totSIM, totREC;
0384 
0385   vector<vector<TH1F*> > ptdistrib;
0386   vector<vector<TH1F*> > etadistrib;
0387   TFile* hFile;
0388 
0389   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken_;
0390   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> theBToken_;
0391 };
0392 
0393 DEFINE_FWK_MODULE(TrackValidator);
[ 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!