Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:16 

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0015 
0016 #include "DataFormats/Math/interface/deltaR.h"
0017 #include "HepMC/GenEvent.h"
0018 #include "HepMC/GenVertex.h"
0019 #include "SimDataFormats/Track/interface/SimTrack.h"
0020 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0021 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0022 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0023 
0024 #include "TFile.h"
0025 #include "TObjArray.h"
0026 #include <TH1.h>
0027 #include <cmath>
0028 #include <iostream>
0029 #include <vector>
0030 
0031 template <class T>
0032 T sqr(T t) {
0033   return t * t;
0034 }
0035 
0036 class PixelTrackVal : public edm::one::EDAnalyzer<> {
0037 public:
0038   explicit PixelTrackVal(const edm::ParameterSet &conf);
0039   ~PixelTrackVal() override;
0040   void beginJob() override;
0041   void analyze(const edm::Event &ev, const edm::EventSetup &es) override;
0042   void endJob() override;
0043 
0044 private:
0045   int verbose_;
0046   std::string file_;
0047   TObjArray hList;
0048   edm::EDGetTokenT<reco::TrackCollection> trackCollectionToken_;
0049   edm::EDGetTokenT<edm::SimTrackContainer> simTrackContainerToken_;
0050   edm::EDGetTokenT<edm::SimVertexContainer> simVertexContainerToken_;
0051 };
0052 
0053 PixelTrackVal::PixelTrackVal(const edm::ParameterSet &conf)
0054     : verbose_(conf.getUntrackedParameter<unsigned int>("Verbosity",
0055                                                         0))  // How noisy?
0056       ,
0057       file_(conf.getUntrackedParameter<std::string>("HistoFile", "pixelTrackHistos.root")),
0058       hList(0),
0059       trackCollectionToken_(
0060           consumes<reco::TrackCollection>(edm::InputTag(conf.getParameter<std::string>("TrackCollection")))),
0061       simTrackContainerToken_(consumes<edm::SimTrackContainer>(conf.getParameter<edm::InputTag>("simG4"))),
0062       simVertexContainerToken_(consumes<edm::SimVertexContainer>(conf.getParameter<edm::InputTag>("simG4"))) {
0063   edm::LogInfo("PixelTrackVal") << " CTOR";
0064 }
0065 
0066 PixelTrackVal::~PixelTrackVal() { edm::LogInfo("PixelTrackVal") << " DTOR"; }
0067 
0068 void PixelTrackVal::beginJob() {
0069   hList.Add(new TH1F("h_Pt", "h_Pt", 31, -2., 1.2));
0070   hList.Add(new TH1F("h_dR", "h_dR", 30, 0., 0.06));
0071   hList.Add(new TH1F("h_TIP", "h_TIP", 100, -0.1, 0.1));
0072   hList.Add(new TH1F("h_VtxZ", "h_VtxZ", 100, -0.1, 0.1));
0073   hList.Add(new TH1F("h_VtxZ_Pull", "h_VtxZ_Pull", 80, 0., 8));
0074   hList.Add(new TH1F("h_Nan", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz", 9, 0.5, 9.5));
0075   hList.SetOwner();
0076 }
0077 
0078 void PixelTrackVal::analyze(const edm::Event &ev, const edm::EventSetup &es) {
0079   LogTrace("PixelTrackVal") << "*** PixelTrackVal, analyze event: " << ev.id() << std::endl;
0080 
0081   //------------------------ simulated tracks
0082   edm::Handle<reco::TrackCollection> trackCollection;
0083   ev.getByToken(trackCollectionToken_, trackCollection);
0084   const reco::TrackCollection tracks = *(trackCollection.product());
0085 
0086   typedef reco::TrackCollection::const_iterator IT;
0087 
0088   if (verbose_ > 0) {
0089     //    std::cout << *(trackCollection.provenance()) << std::endl;
0090     edm::LogInfo("PixelTrackVal") << "Reconstructed " << tracks.size() << " tracks" << std::endl;
0091   }
0092 
0093   for (unsigned int idx = 0; idx < tracks.size(); idx++) {
0094     const reco::Track *it = &tracks[idx];
0095     TH1 *h = static_cast<TH1 *>(hList.FindObject("h_Nan"));
0096     h->Fill(1., edm::isNotFinite(it->momentum().x()) * 1.);
0097     h->Fill(2., edm::isNotFinite(it->momentum().y()) * 1.);
0098     h->Fill(3., edm::isNotFinite(it->momentum().z()) * 1.);
0099 
0100     bool problem = false;
0101     int index = 3;
0102     for (int i = 0; i != 3; i++) {
0103       for (int j = i; j != 3; j++) {
0104         index++;
0105         static_cast<TH1 *>(hList.FindObject("h_Nan"))->Fill(index * 1., edm::isNotFinite(it->covariance(i, j)) * 1.);
0106         if (edm::isNotFinite(it->covariance(i, j)))
0107           problem = true;
0108         // in addition, diagonal element must be positive
0109         if (j == i && it->covariance(i, j) < 0) {
0110           static_cast<TH1 *>(hList.FindObject("h_Nan"))->Fill(index * 1., 1.);
0111           problem = true;
0112         }
0113       }
0114     }
0115     if (problem)
0116       edm::LogInfo("PixelTrackVal") << " *** PROBLEM **" << std::endl;
0117 
0118     if (verbose_ > 0) {
0119       edm::LogInfo("PixelTrackVal") << "\tmomentum: " << tracks[idx].momentum() << "\tPT: " << tracks[idx].pt()
0120                                     << std::endl;
0121       edm::LogInfo("PixelTrackVal") << "\tvertex: " << tracks[idx].vertex() << "\tTIP: " << tracks[idx].d0() << " +- "
0122                                     << tracks[idx].d0Error() << "\tZ0: " << tracks[idx].dz() << " +- "
0123                                     << tracks[idx].dzError() << std::endl;
0124       edm::LogInfo("PixelTrackVal") << "\tcharge: " << tracks[idx].charge() << std::endl;
0125     }
0126   }
0127 
0128   //------------------------ simulated vertices and tracks
0129 
0130   edm::Handle<edm::SimVertexContainer> simVtcs;
0131   ev.getByToken(simVertexContainerToken_, simVtcs);
0132 
0133   //   edm::LogInfo("PixelTrackVal") << "SimVertex " << simVtcs->size() << std::endl;
0134   //   for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
0135   //       v!=simVtcs->end(); ++v){
0136   //     edm::LogInfo("PixelTrackVal") << "simvtx " << std::setw(10) << std::setprecision(3)
0137   //         << v->position().x() << " " << v->position().y() << " " <<
0138   //         v->position().z() << " "
0139   //         << v->parentIndex() << " " << v->noParent() << " " << std::endl; }
0140 
0141   edm::Handle<edm::SimTrackContainer> simTrks;
0142   ev.getByToken(simTrackContainerToken_, simTrks);
0143   edm::LogInfo("PixelTrackVal") << "simtrks " << simTrks->size() << std::endl;
0144 
0145   //-------------- association
0146   // matching cuts from Marcin
0147   float detaMax = 0.012;
0148   float dRMax = 0.025;
0149   typedef edm::SimTrackContainer::const_iterator IP;
0150   for (IP p = simTrks->begin(); p != simTrks->end(); p++) {
0151     if ((*p).noVertex())
0152       continue;
0153     if ((*p).type() == -99)
0154       continue;
0155     if ((*p).vertIndex() != 0)
0156       continue;
0157 
0158     math::XYZVector mom_gen((*p).momentum().x(), (*p).momentum().y(), (*p).momentum().z());
0159     float phi_gen = (*p).momentum().phi();
0160     //    float pt_gen = (*p).momentum().Pt();
0161     float pt_gen = sqrt((*p).momentum().x() * (*p).momentum().x() + (*p).momentum().y() * (*p).momentum().y());
0162     float eta_gen = (*p).momentum().eta();
0163     math::XYZTLorentzVectorD vtx((*simVtcs)[p->vertIndex()].position().x(),
0164                                  (*simVtcs)[p->vertIndex()].position().y(),
0165                                  (*simVtcs)[p->vertIndex()].position().z(),
0166                                  (*simVtcs)[p->vertIndex()].position().e());
0167     float z_gen = vtx.z();
0168 
0169     //     cout << "\tmomentum: " <<  (*p).momentum()
0170     //          <<" vtx: "<<(*p).vertIndex()<<" type: "<<(*p).type()
0171     //          << endl;
0172 
0173     typedef reco::TrackCollection::const_iterator IT;
0174     for (IT it = tracks.begin(); it != tracks.end(); it++) {
0175       math::XYZVector mom_rec = (*it).momentum();
0176       float phi_rec = (*it).momentum().phi();
0177       float pt_rec = (*it).pt();
0178       float z_rec = (*it).vertex().z();
0179       float eta_rec = (*it).momentum().eta();
0180       //    float chi2   = (*it).chi2();
0181       float dphi = phi_gen - phi_rec;
0182       while (dphi > M_PI)
0183         dphi -= 2 * M_PI;
0184       while (dphi < -M_PI)
0185         dphi += 2 * M_PI;
0186       float deta = eta_gen - eta_rec;
0187       float dz = z_gen - z_rec;
0188       double dR = deltaR(mom_gen, mom_rec);
0189       //
0190       // matched track
0191       //
0192       if (fabs(deta) < 0.3 && fabs(dphi) < 0.3)
0193         static_cast<TH1 *>(hList.FindObject("h_dR"))->Fill(dR);
0194       if (fabs(deta) < detaMax && dR < dRMax) {
0195         static_cast<TH1 *>(hList.FindObject("h_Pt"))->Fill((pt_gen - pt_rec) / pt_gen);
0196         static_cast<TH1 *>(hList.FindObject("h_TIP"))->Fill(it->d0());
0197         static_cast<TH1 *>(hList.FindObject("h_VtxZ"))->Fill(dz);
0198         static_cast<TH1 *>(hList.FindObject("h_VtxZ_Pull"))->Fill(fabs(dz / it->dzError()));
0199       }
0200     }
0201   }
0202 }
0203 
0204 void PixelTrackVal::endJob() {
0205   // Make my little tree
0206   TFile f(file_.c_str(), "RECREATE");
0207   hList.Write();
0208   f.Close();
0209 }
0210 
0211 // float PixelTrackVal::deltaRR(const  math::XYZVector & m1, const
0212 // math::XYZVector & m2) const
0213 //{
0214 //  float dphi = m1.phi()-m2.phi();
0215 //  while (dphi > 2*M_PI) dphi-=2*M_PI;
0216 //  while (dphi < -2*M_PI) dphi+=2*M_PI;
0217 //  float deta = m1.eta() - m2.eta();
0218 //  float dr = sqrt( sqr(dphi) + sqr(deta));
0219 //  return dr;
0220 //}
0221 
0222 DEFINE_FWK_MODULE(PixelTrackVal);
[ 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!