Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-07 02:51:00 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ValidationMisalignedTracker
0004 // Class:      ValidationMisalignedTracker
0005 //
0006 /**\class ValidationMisalignedTracker ValidationMisalignedTracker.cc Alignment/OfflineValidation/src/ValidationMisalignedTracker.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Nicola De Filippis
0015 //         Created:  Thu Dec 14 13:13:32 CET 2006
0016 // $Id: ValidationMisalignedTracker.cc,v 1.8 2013/01/07 20:46:23 wmtan Exp $
0017 //
0018 //
0019 
0020 #include "Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.h"
0021 
0022 // user include files
0023 
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Utilities/interface/transform.h"
0028 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0029 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0030 
0031 //
0032 // constructors and destructor
0033 //
0034 ValidationMisalignedTracker::ValidationMisalignedTracker(const edm::ParameterSet& iConfig)
0035     : geomToken_(esConsumes()),
0036       magFieldToken_(esConsumes()),
0037       selection_eff(iConfig.getUntrackedParameter<bool>("selection_eff", false)),
0038       selection_fake(iConfig.getUntrackedParameter<bool>("selection_fake", true)),
0039       ZmassSelection_(iConfig.getUntrackedParameter<bool>("ZmassSelection", false)),
0040       simobject(iConfig.getUntrackedParameter<std::string>("simobject", "g4SimHits")),
0041       trackassociator(iConfig.getUntrackedParameter<std::string>("TrackAssociator", "ByHits")),
0042       associators(iConfig.getParameter<std::vector<std::string> >("associators")),
0043       label(iConfig.getParameter<std::vector<edm::InputTag> >("label")),
0044       label_tp_effic(iConfig.getParameter<edm::InputTag>("label_tp_effic")),
0045       label_tp_fake(iConfig.getParameter<edm::InputTag>("label_tp_fake")),
0046 
0047       rootfile_(iConfig.getUntrackedParameter<std::string>("rootfile", "myroot.root")),
0048       evtToken_(consumes<edm::HepMCProduct>(edm::InputTag("source"))),
0049       tpeffToken_(consumes<TrackingParticleCollection>(label_tp_effic)),
0050       tpfakeToken_(consumes<TrackingParticleCollection>(label_tp_fake)),
0051       trackToken_(consumes<edm::View<reco::Track> >(label[0])),
0052       assocToken_{edm::vector_transform(
0053           label, [this](const edm::InputTag& lab) { return consumes<reco::TrackToTrackingParticleAssociator>(lab); })} {
0054   //now do what ever initialization is needed
0055   mzmu = 0., recmzmu = 0., ptzmu = 0., recptzmu = 0., etazmu = 0., recetazmu = 0., thetazmu = 0., recthetazmu = 0.,
0056   phizmu = 0., recphizmu = 0.;
0057   recenezmu = 0., enezmu = 0., pLzmu = 0., recpLzmu = 0., yzmu = 0., recyzmu = 0., mxptmu = 0., recmxptmu = 0.,
0058   minptmu = 0., recminptmu = 0.;
0059   // mzele=0.,recmzele=0.
0060 
0061   flag = 0, flagrec = 0, count = 0, countrec = 0;
0062   nAssoc = 0;
0063 
0064   for (int i = 0; i < 2; i++) {
0065     countpart[i] = 0;
0066     countpartrec[i] = 0;
0067     for (int j = 0; j < 2; j++) {
0068       ene[i][j] = 0.;
0069       p[i][j] = 0.;
0070       px[i][j] = 0.;
0071       py[i][j] = 0.;
0072       pz[i][j] = 0.;
0073       ptmu[i][j] = 0.;
0074       recene[i][j] = 0.;
0075       recp[i][j] = 0.;
0076       recpx[i][j] = 0.;
0077       recpy[i][j] = 0.;
0078       recpz[i][j] = 0.;
0079       recptmu[i][j] = 0.;
0080     }
0081   }
0082 
0083   eventCount_ = 0;
0084 
0085   file_ = new TFile(rootfile_.c_str(), "RECREATE");
0086 
0087   // initialize the tree
0088   tree_eff = new TTree("EffTracks", "Efficiency Tracks Tree");
0089 
0090   tree_eff->Branch("Run", &irun, "irun/i");
0091   tree_eff->Branch("Event", &ievt, "ievt/i");
0092 
0093   // SimTrack
0094   tree_eff->Branch("TrackID", &trackType, "trackType/i");
0095   tree_eff->Branch("pt", &pt, "pt/F");
0096   tree_eff->Branch("eta", &eta, "eta/F");
0097   tree_eff->Branch("CotTheta", &cottheta, "cottheta/F");
0098   tree_eff->Branch("phi", &phi, "phi/F");
0099   tree_eff->Branch("d0", &d0, "d0/F");
0100   tree_eff->Branch("z0", &z0, "z0/F");
0101   tree_eff->Branch("nhit", &nhit, "nhit/i");
0102 
0103   // RecTrack
0104   tree_eff->Branch("recpt", &recpt, "recpt/F");
0105   tree_eff->Branch("receta", &receta, "receta/F");
0106   tree_eff->Branch("CotRecTheta", &reccottheta, "reccottheta/F");
0107   tree_eff->Branch("recphi", &recphi, "recphi/F");
0108   tree_eff->Branch("recd0", &recd0, "recd0/F");
0109   tree_eff->Branch("recz0", &recz0, "recz0/F");
0110   tree_eff->Branch("nAssoc", &nAssoc, "nAssoc/i");
0111   tree_eff->Branch("recnhit", &recnhit, "recnhit/i");
0112   tree_eff->Branch("CHISQ", &recchiq, "recchiq/F");
0113 
0114   tree_eff->Branch("reseta", &reseta, "reseta/F");
0115   tree_eff->Branch("respt", &respt, "respt/F");
0116   tree_eff->Branch("resd0", &resd0, "resd0/F");
0117   tree_eff->Branch("resz0", &resz0, "resz0/F");
0118   tree_eff->Branch("resphi", &resphi, "resphi/F");
0119   tree_eff->Branch("rescottheta", &rescottheta, "rescottheta/F");
0120   tree_eff->Branch("eff", &eff, "eff/F");
0121 
0122   // Invariant masses, pt of Z
0123   tree_eff->Branch("mzmu", &mzmu, "mzmu/F");
0124   tree_eff->Branch("ptzmu", &ptzmu, "ptzmu/F");
0125   tree_eff->Branch("pLzmu", &pLzmu, "pLzmu/F");
0126   tree_eff->Branch("enezmu", &enezmu, "enezmu/F");
0127   tree_eff->Branch("etazmu", &etazmu, "etazmu/F");
0128   tree_eff->Branch("thetazmu", &thetazmu, "thetazmu/F");
0129   tree_eff->Branch("phizmu", &phizmu, "phizmu/F");
0130   tree_eff->Branch("yzmu", &yzmu, "yzmu/F");
0131   tree_eff->Branch("mxptmu", &mxptmu, "mxptmu/F");
0132   tree_eff->Branch("minptmu", &minptmu, "minptmu/F");
0133 
0134   tree_eff->Branch("recmzmu", &recmzmu, "recmzmu/F");
0135   tree_eff->Branch("recptzmu", &recptzmu, "recptzmu/F");
0136   tree_eff->Branch("recpLzmu", &recpLzmu, "recpLzmu/F");
0137   tree_eff->Branch("recenezmu", &recenezmu, "recenezmu/F");
0138   tree_eff->Branch("recetazmu", &recetazmu, "recetazmu/F");
0139   tree_eff->Branch("recthetazmu", &recthetazmu, "recthetazmu/F");
0140   tree_eff->Branch("recphizmu", &recphizmu, "recphizmu/F");
0141   tree_eff->Branch("recyzmu", &recyzmu, "recyzmu/F");
0142   tree_eff->Branch("recmxptmu", &recmxptmu, "recmxptmu/F");
0143   tree_eff->Branch("recminptmu", &recminptmu, "recminptmu/F");
0144 
0145   //tree->Branch("mzele",&ntmzele,"ntmzele/F");
0146   //tree->Branch("recmzele",&ntmzeleRec,"ntmzeleRec/F");
0147   tree_eff->Branch("chi2Associator", &recchiq, "recchiq/F");
0148 
0149   // Fake
0150 
0151   tree_fake = new TTree("FakeTracks", "Fake Rate Tracks Tree");
0152 
0153   tree_fake->Branch("Run", &irun, "irun/i");
0154   tree_fake->Branch("Event", &ievt, "ievt/i");
0155 
0156   // SimTrack
0157   tree_fake->Branch("fakeTrackID", &faketrackType, "faketrackType/i");
0158   tree_fake->Branch("fakept", &fakept, "fakept/F");
0159   tree_fake->Branch("fakeeta", &fakeeta, "fakeeta/F");
0160   tree_fake->Branch("fakeCotTheta", &fakecottheta, "fakecottheta/F");
0161   tree_fake->Branch("fakephi", &fakephi, "fakephi/F");
0162   tree_fake->Branch("faked0", &faked0, "faked0/F");
0163   tree_fake->Branch("fakez0", &fakez0, "fakez0/F");
0164   tree_fake->Branch("fakenhit", &fakenhit, "fakenhit/i");
0165 
0166   // RecTrack
0167   tree_fake->Branch("fakerecpt", &fakerecpt, "fakerecpt/F");
0168   tree_fake->Branch("fakereceta", &fakereceta, "fakereceta/F");
0169   tree_fake->Branch("fakeCotRecTheta", &fakereccottheta, "fakereccottheta/F");
0170   tree_fake->Branch("fakerecphi", &fakerecphi, "fakerecphi/F");
0171   tree_fake->Branch("fakerecd0", &fakerecd0, "fakerecd0/F");
0172   tree_fake->Branch("fakerecz0", &fakerecz0, "fakerecz0/F");
0173   tree_fake->Branch("fakenAssoc", &fakenAssoc, "fakenAssoc/i");
0174   tree_fake->Branch("fakerecnhit", &fakerecnhit, "fakerecnhit/i");
0175   tree_fake->Branch("fakeCHISQ", &fakerecchiq, "fakerecchiq/F");
0176 
0177   tree_fake->Branch("fakereseta", &fakereseta, "fakereseta/F");
0178   tree_fake->Branch("fakerespt", &fakerespt, "fakerespt/F");
0179   tree_fake->Branch("fakeresd0", &fakeresd0, "fakeresd0/F");
0180   tree_fake->Branch("fakeresz0", &fakeresz0, "fakeresz0/F");
0181   tree_fake->Branch("fakeresphi", &fakeresphi, "fakeresphi/F");
0182   tree_fake->Branch("fakerescottheta", &fakerescottheta, "fakerescottheta/F");
0183   tree_fake->Branch("fake", &fake, "fake/F");
0184 
0185   // Invariant masses, pt of Z
0186   tree_fake->Branch("fakemzmu", &fakemzmu, "fakemzmu/F");
0187   tree_fake->Branch("fakeptzmu", &fakeptzmu, "fakeptzmu/F");
0188   tree_fake->Branch("fakepLzmu", &fakepLzmu, "fakepLzmu/F");
0189   tree_fake->Branch("fakeenezmu", &fakeenezmu, "fakeenezmu/F");
0190   tree_fake->Branch("fakeetazmu", &fakeetazmu, "fakeetazmu/F");
0191   tree_fake->Branch("fakethetazmu", &fakethetazmu, "fakethetazmu/F");
0192   tree_fake->Branch("fakephizmu", &fakephizmu, "fakephizmu/F");
0193   tree_fake->Branch("fakeyzmu", &fakeyzmu, "fakeyzmu/F");
0194   tree_fake->Branch("fakemxptmu", &fakemxptmu, "fakemxptmu/F");
0195   tree_fake->Branch("fakeminptmu", &fakeminptmu, "fakeminptmu/F");
0196 
0197   tree_fake->Branch("fakerecmzmu", &fakerecmzmu, "fakerecmzmu/F");
0198   tree_fake->Branch("fakerecptzmu", &fakerecptzmu, "fakerecptzmu/F");
0199   tree_fake->Branch("fakerecpLzmu", &fakerecpLzmu, "fakerecpLzmu/F");
0200   tree_fake->Branch("fakerecenezmu", &fakerecenezmu, "fakerecenezmu/F");
0201   tree_fake->Branch("fakerecetazmu", &fakerecetazmu, "fakerecetazmu/F");
0202   tree_fake->Branch("fakerecthetazmu", &fakerecthetazmu, "fakerecthetazmu/F");
0203   tree_fake->Branch("fakerecphizmu", &fakerecphizmu, "fakerecphizmu/F");
0204   tree_fake->Branch("fakerecyzmu", &fakerecyzmu, "fakerecyzmu/F");
0205   tree_fake->Branch("fakerecmxptmu", &fakerecmxptmu, "fakerecmxptmu/F");
0206   tree_fake->Branch("fakerecminptmu", &fakerecminptmu, "fakerecminptmu/F");
0207 
0208   tree_fake->Branch("fakechi2Associator", &fakerecchiq, "fakerecchiq/F");
0209 }
0210 
0211 ValidationMisalignedTracker::~ValidationMisalignedTracker() {
0212   edm::LogVerbatim("ValidationMisalignedTracker")
0213       << "ValidationMisalignedTracker::endJob Processed " << eventCount_ << " events";
0214 
0215   // store the tree in the output file
0216   file_->Write();
0217 
0218   // Closing the file deletes the tree.
0219   file_->Close();
0220   tree_eff = nullptr;
0221   tree_fake = nullptr;
0222 }
0223 
0224 //
0225 // member functions
0226 //
0227 void ValidationMisalignedTracker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0228   edm::ParameterSetDescription desc;
0229   desc.add<std::vector<edm::InputTag> >("label", {});
0230   desc.add<edm::InputTag>("label_tp_effic", edm::InputTag("FinalTrackRefitter"));
0231   desc.add<edm::InputTag>("label_tp_fake", edm::InputTag("TrackRefitter"));
0232   desc.add<std::vector<std::string> >("associators", {});
0233   desc.addUntracked<bool>("selection_eff", false);
0234   desc.addUntracked<bool>("selection_fake", true);
0235   desc.addUntracked<bool>("ZmassSelection", false);
0236   desc.addUntracked<std::string>("simobject", "g4SimHits");
0237   desc.addUntracked<std::string>("TrackAssociator", "ByHits");
0238   desc.addUntracked<std::string>("rootfile", "myroot.root");
0239   descriptions.add("validationMisAlignedTracker", desc);
0240 }
0241 
0242 // ------------ method called to for each event  ------------
0243 void ValidationMisalignedTracker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0244   std::vector<const reco::TrackToTrackingParticleAssociator*> associatore;
0245 
0246   {
0247     for (unsigned int w = 0; w < associators.size(); w++) {
0248       const edm::Handle<reco::TrackToTrackingParticleAssociator>& theAssociator = iEvent.getHandle(assocToken_[w]);
0249       associatore.push_back(theAssociator.product());
0250     }
0251   }
0252 
0253   edm::LogInfo("Tracker Misalignment Validation") << "\n Starting!";
0254 
0255   // Monte Carlo Z selection
0256   skip = false;
0257   std::vector<int> indmu;
0258 
0259   if (selection_eff && ZmassSelection_) {
0260     const edm::Handle<edm::HepMCProduct>& evt = iEvent.getHandle(evtToken_);
0261     bool accepted = false;
0262     bool foundmuons = false;
0263     HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0264 
0265     for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0266       if (!accepted && ((*p)->pdg_id() == 23) && (*p)->status() == 3) {
0267         accepted = true;
0268         for (HepMC::GenVertex::particle_iterator aDaughter = (*p)->end_vertex()->particles_begin(HepMC::descendants);
0269              aDaughter != (*p)->end_vertex()->particles_end(HepMC::descendants);
0270              aDaughter++) {
0271           if (abs((*aDaughter)->pdg_id()) == 13) {
0272             foundmuons = true;
0273             if ((*aDaughter)->status() != 1) {
0274               for (HepMC::GenVertex::particle_iterator byaDaughter =
0275                        (*aDaughter)->end_vertex()->particles_begin(HepMC::descendants);
0276                    byaDaughter != (*aDaughter)->end_vertex()->particles_end(HepMC::descendants);
0277                    byaDaughter++) {
0278                 if ((*byaDaughter)->status() == 1 && abs((*byaDaughter)->pdg_id()) == 13) {
0279                   indmu.push_back((*byaDaughter)->barcode());
0280                   edm::LogVerbatim("ValidationMisalignedTracker")
0281                       << "Stable muon from Z with charge " << (*byaDaughter)->pdg_id() << " and index "
0282                       << (*byaDaughter)->barcode();
0283                 }
0284               }
0285             } else {
0286               indmu.push_back((*aDaughter)->barcode());
0287               edm::LogVerbatim("ValidationMisalignedTracker")
0288                   << "Stable muon from Z with charge " << (*aDaughter)->pdg_id() << " and index "
0289                   << (*aDaughter)->barcode();
0290             }
0291           }
0292         }
0293         if (!foundmuons) {
0294           edm::LogVerbatim("ValidationMisalignedTracker") << "No muons from Z ...skip event";
0295           skip = true;
0296         }
0297       }
0298     }
0299     if (!accepted) {
0300       edm::LogVerbatim("ValidationMisalignedTracker") << "No Z particles in the event ...skip event";
0301       skip = true;
0302     }
0303   } else {
0304     skip = false;
0305   }
0306 
0307   //
0308   // Retrieve tracker geometry from event setup
0309   //
0310   const TrackerGeometry* trackerGeometry = &iSetup.getData(geomToken_);
0311   auto testGeomDet = trackerGeometry->detsTOB().front();
0312   edm::LogVerbatim("ValidationMisalignedTracker") << testGeomDet->position();
0313 
0314   //Dump Run and Event
0315   irun = iEvent.id().run();
0316   ievt = iEvent.id().event();
0317 
0318   // Reset tree variables
0319   int countpart[2] = {0, 0}, countpartrec[2] = {0, 0}, flag = 0, flagrec = 0, count = 0, countrec = 0;
0320   //int countsim=0;
0321   float ene[2][2], px[2][2], py[2][2], pz[2][2], ptmu[2][2];
0322   float recene[2][2], recp[2][2], recpx[2][2], recpy[2][2], recpz[2][2], recptmu[2][2];
0323 
0324   for (int i = 0; i < 2; i++) {
0325     for (int j = 0; j < 2; j++) {
0326       ene[i][j] = 0.;
0327       px[i][j] = 0.;
0328       py[i][j] = 0.;
0329       pz[i][j] = 0.;
0330       ptmu[i][j] = 0.;
0331       recene[i][j] = 0.;
0332       recp[i][j] = 0.;
0333       recpx[i][j] = 0.;
0334       recpy[i][j] = 0.;
0335       recpz[i][j] = 0.;
0336       recptmu[i][j] = 0.;
0337     }
0338   }
0339 
0340   const edm::Handle<TrackingParticleCollection>& TPCollectionHeff = iEvent.getHandle(tpeffToken_);
0341   const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
0342 
0343   const edm::Handle<TrackingParticleCollection>& TPCollectionHfake = iEvent.getHandle(tpfakeToken_);
0344   const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
0345 
0346   int w = 0;
0347   for (unsigned int ww = 0; ww < associators.size(); ww++) {
0348     //
0349     //get collections from the event
0350     //
0351 
0352     const edm::Handle<edm::View<reco::Track> >& trackCollection = iEvent.getHandle(trackToken_);
0353     const edm::View<reco::Track> tC = *(trackCollection.product());
0354 
0355     //associate tracks
0356     LogTrace("TrackValidator") << "Calling associateRecoToSim method\n";
0357     reco::RecoToSimCollection recSimColl = associatore[ww]->associateRecoToSim(trackCollection, TPCollectionHfake);
0358 
0359     LogTrace("TrackValidator") << "Calling associateSimToReco method\n";
0360     reco::SimToRecoCollection simRecColl = associatore[ww]->associateSimToReco(trackCollection, TPCollectionHeff);
0361 
0362     //
0363     //compute number of tracks per eta interval
0364     //
0365 
0366     if (selection_eff && !skip) {
0367       edm::LogVerbatim("ValidationMisalignedTracker") << "Computing Efficiency";
0368 
0369       edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles (before cuts): " << tPCeff.size() << "\n";
0370       int ats = 0;
0371       int st = 0;
0372       for (TrackingParticleCollection::size_type i = 0; i < tPCeff.size(); i++) {
0373         // Initialize variables
0374         eta = 0., theta = 0., phi = 0., pt = 0., cottheta = 0., costheta = 0.;
0375         d0 = 0., z0 = 0.;
0376         nhit = 0;
0377         receta = 0., rectheta = 0., recphi = 0., recpt = 0., reccottheta = 0., recd0 = 0., recz0 = 0.;
0378         respt = 0., resd0 = 0., resz0 = 0., reseta = 0., resphi = 0., rescottheta = 0.;
0379         recchiq = 0.;
0380         recnhit = 0;
0381         trackType = 0;
0382         eff = 0;
0383 
0384         // typedef edm::Ref<TrackingParticleCollection> TrackingParticleRef;
0385         TrackingParticleRef tp(TPCollectionHeff, i);
0386         if (tp->charge() == 0)
0387           continue;
0388         st++;
0389         //pt=sqrt(tp->momentum().perp2());
0390         //eta=tp->momentum().eta();
0391         //vpos=tp->vertex().perp2()));
0392 
0393         const SimTrack* simulatedTrack = &(*tp->g4Track_begin());
0394 
0395         const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0396         FreeTrajectoryState ftsAtProduction(
0397             GlobalPoint(tp->vertex().x(), tp->vertex().y(), tp->vertex().z()),
0398             GlobalVector(
0399                 simulatedTrack->momentum().x(), simulatedTrack->momentum().y(), simulatedTrack->momentum().z()),
0400             TrackCharge(tp->charge()),
0401             theMF);
0402         TSCPBuilderNoMaterial tscpBuilder;
0403         TrajectoryStateClosestToPoint tsAtClosestApproach =
0404             tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
0405         GlobalPoint v = tsAtClosestApproach.theState().position();
0406         GlobalVector p = tsAtClosestApproach.theState().momentum();
0407 
0408         //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
0409         //  double lambdaSim = M_PI/2-p.theta();
0410         //  double phiSim    = p.phi();
0411         double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0412         double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0413         d0 = float(-dxySim);
0414         z0 = float(dszSim * p.mag() / p.perp());
0415 
0416         if (abs(simulatedTrack->type()) == 13 && simulatedTrack->genpartIndex() != -1) {
0417           edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA SIM DI MUONI ";
0418           edm::LogVerbatim("ValidationMisalignedTracker") << "Gen part " << simulatedTrack->genpartIndex();
0419           trackType = simulatedTrack->type();
0420           theta = simulatedTrack->momentum().theta();
0421           costheta = cos(theta);
0422           cottheta = 1. / tan(theta);
0423 
0424           eta = simulatedTrack->momentum().eta();
0425           phi = simulatedTrack->momentum().phi();
0426           pt = simulatedTrack->momentum().pt();
0427           nhit = tp->matchedHit();
0428 
0429           edm::LogVerbatim("ValidationMisalignedTracker")
0430               << "3) Before assoc: SimTrack of type = " << simulatedTrack->type() << " ,at eta = " << eta
0431               << " ,with pt at vertex = " << simulatedTrack->momentum().pt() << " GeV/c"
0432               << " ,d0 =" << d0 << " ,z0 =" << z0 << " ,nhit=" << nhit;
0433 
0434           if (ZmassSelection_) {
0435             if (abs(trackType) == 13 &&
0436                 (simulatedTrack->genpartIndex() == indmu[0] || simulatedTrack->genpartIndex() == indmu[1])) {
0437               edm::LogVerbatim("ValidationMisalignedTracker") << " TRACK sim of muons from Z ";
0438               flag = 0;
0439               count = countpart[0];
0440               countpart[0]++;
0441             } else if (abs(trackType) == 11) {
0442               //edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA SIM DI ELETTRONI ";
0443               flag = 1;
0444               count = countpart[1];
0445               countpart[1]++;
0446             }
0447 
0448             px[flag][count] = simulatedTrack->momentum().x();
0449             py[flag][count] = simulatedTrack->momentum().y();
0450             pz[flag][count] = simulatedTrack->momentum().z();
0451             ptmu[flag][count] = simulatedTrack->momentum().pt();
0452             ene[flag][count] = simulatedTrack->momentum().e();
0453           }
0454 
0455           std::vector<std::pair<edm::RefToBase<reco::Track>, double> > rt;
0456           if (simRecColl.find(tp) != simRecColl.end()) {
0457             rt = simRecColl[tp];
0458             if (!rt.empty()) {
0459               edm::RefToBase<reco::Track> t = rt.begin()->first;
0460               ats++;
0461 
0462               // bool flagptused=false;
0463               // for (unsigned int j=0;j<ptused.size();j++){
0464               //   if (fabs(t->pt()-ptused[j])<0.001) {
0465               //     flagptused=true;
0466               //   }
0467               // }
0468 
0469               edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st << " with pt=" << t->pt()
0470                                                  << " associated with quality:" << rt.begin()->second << "\n";
0471               edm::LogVerbatim("ValidationMisalignedTracker") << "Reconstructed Track:" << t->pt();
0472               edm::LogVerbatim("ValidationMisalignedTracker") << "\tpT: " << t->pt();
0473               edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:d0: " << t->d0();
0474               edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:z0: " << t->dz();
0475               edm::LogVerbatim("ValidationMisalignedTracker")
0476                   << "\tAzimuthal angle of point of closest approach:" << t->phi();
0477               edm::LogVerbatim("ValidationMisalignedTracker") << "\tcharge: " << t->charge();
0478               edm::LogVerbatim("ValidationMisalignedTracker") << "\teta: " << t->eta();
0479               edm::LogVerbatim("ValidationMisalignedTracker") << "\tnormalizedChi2: " << t->normalizedChi2();
0480 
0481               recnhit = t->numberOfValidHits();
0482               recchiq = t->normalizedChi2();
0483               rectheta = t->theta();
0484               reccottheta = 1. / tan(rectheta);
0485               //receta=-log(tan(rectheta/2.));
0486               receta = t->momentum().eta();
0487               //       reccostheta=cos(matchedrectrack->momentum().theta());
0488               recphi = t->phi();
0489               recpt = t->pt();
0490               ptused.push_back(recpt);
0491               recd0 = t->d0();
0492               recz0 = t->dz();
0493 
0494               edm::LogVerbatim("ValidationMisalignedTracker")
0495                   << "5) After call to associator: the best match has " << recnhit << " hits, Chi2 = " << recchiq
0496                   << ", pt at vertex = " << recpt << " GeV/c, "
0497                   << ", recd0 = " << recd0 << ", recz0= " << recz0;
0498 
0499               respt = recpt - pt;
0500               resd0 = recd0 - d0;
0501               resz0 = recz0 - z0;
0502               reseta = receta - eta;
0503               resphi = recphi - phi;
0504               rescottheta = reccottheta - cottheta;
0505               eff = 1;
0506 
0507               edm::LogVerbatim("ValidationMisalignedTracker")
0508                   << "6) Transverse momentum residual=" << respt << " ,d0 residual=" << resd0
0509                   << " ,z0 residual=" << resz0 << " with eff=" << eff;
0510 
0511               if (ZmassSelection_) {
0512                 if (abs(trackType) == 13) {
0513                   edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA RECO DI MUONI ";
0514                   flagrec = 0;
0515                   countrec = countpartrec[0];
0516                   countpartrec[0]++;
0517                 } else if (abs(trackType) == 11) {
0518                   edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA RECO DI ELETTRONI ";
0519                   flagrec = 1;
0520                   countrec = countpartrec[1];
0521                   countpartrec[1]++;
0522                 }
0523 
0524                 recp[flagrec][countrec] = sqrt(t->momentum().mag2());
0525                 recpx[flagrec][countrec] = t->momentum().x();
0526                 recpy[flagrec][countrec] = t->momentum().y();
0527                 recpz[flagrec][countrec] = t->momentum().z();
0528                 recptmu[flagrec][countrec] =
0529                     sqrt((t->momentum().x() * t->momentum().x()) + (t->momentum().y() * t->momentum().y()));
0530                 if (abs(trackType) == 13)
0531                   recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.105 * 0.105);
0532                 if (abs(trackType) == 11)
0533                   recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.0005 * 0.0005);
0534               }
0535 
0536               edm::LogVerbatim("ValidationMisalignedTracker") << "7) Transverse momentum reconstructed =" << recpt
0537                                                               << " at  eta= " << receta << " and phi= " << recphi;
0538             }
0539           } else {
0540             edm::LogVerbatim("TrackValidator")
0541                 << "TrackingParticle #" << st << " with pt=" << sqrt(tp->momentum().perp2())
0542                 << " NOT associated to any reco::Track"
0543                 << "\n";
0544             receta = -100.;
0545             recphi = -100.;
0546             recpt = -100.;
0547             recd0 = -100.;
0548             recz0 = -100;
0549             respt = -100.;
0550             resd0 = -100.;
0551             resz0 = -100.;
0552             resphi = -100.;
0553             reseta = -100.;
0554             rescottheta = -100.;
0555             recnhit = 100;
0556             recchiq = -100;
0557             eff = 0;
0558             flagrec = 100;
0559           }
0560 
0561           edm::LogVerbatim("ValidationMisalignedTracker") << "Eff=" << eff;
0562 
0563           // simulated muons
0564 
0565           edm::LogVerbatim("ValidationMisalignedTracker") << "Flag is" << flag;
0566           edm::LogVerbatim("ValidationMisalignedTracker") << "RecFlag is" << flagrec;
0567 
0568           if (countpart[0] == 2 && flag == 0) {
0569             mzmu =
0570                 sqrt((ene[0][0] + ene[0][1]) * (ene[0][0] + ene[0][1]) - (px[0][0] + px[0][1]) * (px[0][0] + px[0][1]) -
0571                      (py[0][0] + py[0][1]) * (py[0][0] + py[0][1]) - (pz[0][0] + pz[0][1]) * (pz[0][0] + pz[0][1]));
0572             edm::LogVerbatim("ValidationMisalignedTracker") << "Mzmu " << mzmu;
0573             ptzmu = sqrt((px[0][0] + px[0][1]) * (px[0][0] + px[0][1]) + (py[0][0] + py[0][1]) * (py[0][0] + py[0][1]));
0574 
0575             pLzmu = pz[0][0] + pz[0][1];
0576             enezmu = ene[0][0] + ene[0][1];
0577             phizmu = atan2((py[0][0] + py[0][1]), (px[0][0] + px[0][1]));
0578             thetazmu = atan2(ptzmu, (pz[0][0] + pz[0][1]));
0579             etazmu = -log(tan(thetazmu * 3.14 / 360.));
0580             yzmu = 0.5 * log((enezmu + pLzmu) / (enezmu - pLzmu));
0581             mxptmu = std::max(ptmu[0][0], ptmu[0][1]);
0582             minptmu = std::min(ptmu[0][0], ptmu[0][1]);
0583           } else {
0584             mzmu = -100.;
0585             ptzmu = -100.;
0586             pLzmu = -100.;
0587             enezmu = -100.;
0588             etazmu = -100.;
0589             phizmu = -100.;
0590             thetazmu = -100.;
0591             yzmu = -100.;
0592             mxptmu = -100.;
0593             minptmu = -100.;
0594           }
0595 
0596           // reconstructed muons
0597           if (countpartrec[0] == 2 && flagrec == 0) {
0598             recmzmu = sqrt((recene[0][0] + recene[0][1]) * (recene[0][0] + recene[0][1]) -
0599                            (recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) -
0600                            (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]) -
0601                            (recpz[0][0] + recpz[0][1]) * (recpz[0][0] + recpz[0][1]));
0602             edm::LogVerbatim("ValidationMisalignedTracker") << "RecMzmu " << recmzmu;
0603             recptzmu = sqrt((recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) +
0604                             (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]));
0605 
0606             recpLzmu = recpz[0][0] + recpz[0][1];
0607             recenezmu = recene[0][0] + recene[0][1];
0608             recphizmu = atan2((recpy[0][0] + recpy[0][1]), (recpx[0][0] + recpx[0][1]));
0609             recthetazmu = atan2(recptzmu, (recpz[0][0] + recpz[0][1]));
0610             recetazmu = -log(tan(recthetazmu * 3.14 / 360.));
0611             recyzmu = 0.5 * log((recenezmu + recpLzmu) / (recenezmu - recpLzmu));
0612             recmxptmu = std::max(recptmu[0][0], recptmu[0][1]);
0613             recminptmu = std::min(recptmu[0][0], recptmu[0][1]);
0614           } else {
0615             recmzmu = -100.;
0616             recptzmu = -100.;
0617             recpLzmu = -100.;
0618             recenezmu = -100.;
0619             recetazmu = -100.;
0620             recphizmu = -100.;
0621             recthetazmu = -100.;
0622             recyzmu = -100.;
0623             recmxptmu = -100;
0624             recminptmu = -100.;
0625           }
0626 
0627           tree_eff->Fill();
0628 
0629         }  // end of loop on muons
0630       }    // end of loop for tracking particle
0631     }      // end of loop for efficiency
0632 
0633     //
0634     // Fake Rate
0635     //
0636     if (selection_fake) {
0637       edm::LogVerbatim("ValidationMisalignedTracker") << "Computing Fake Rate";
0638 
0639       fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0640       faked0 = 0., fakez0 = 0.;
0641       fakenhit = 0;
0642       fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0643       fakerecz0 = 0.;
0644       fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0645       fakerecchiq = 0.;
0646       fakerecnhit = 0;
0647       faketrackType = 0;
0648       fake = 0;
0649 
0650       //      int at=0;
0651       int rT = 0;
0652       for (reco::TrackCollection::size_type i = 0; i < tC.size(); ++i) {
0653         edm::RefToBase<reco::Track> track(trackCollection, i);
0654         rT++;
0655 
0656         fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0657         faked0 = 0., fakez0 = 0.;
0658         fakenhit = 0;
0659         fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0660         fakerecz0 = 0.;
0661         fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0662         fakerecchiq = 0.;
0663         fakerecnhit = 0;
0664         faketrackType = 0;
0665         fake = 0;
0666 
0667         fakerecnhit = track->numberOfValidHits();
0668         fakerecchiq = track->normalizedChi2();
0669         fakerectheta = track->theta();
0670         fakereccottheta = 1. / tan(rectheta);
0671         //fakereceta=-log(tan(rectheta/2.));
0672         fakereceta = track->momentum().eta();
0673         //     fakereccostheta=cos(track->momentum().theta());
0674         fakerecphi = track->phi();
0675         fakerecpt = track->pt();
0676         fakerecd0 = track->d0();
0677         fakerecz0 = track->dz();
0678 
0679         edm::LogVerbatim("ValidationMisalignedTracker") << "1) Before assoc: TkRecTrack at eta = " << fakereceta;
0680         edm::LogVerbatim("ValidationMisalignedTracker") << "Track number " << i;
0681         edm::LogVerbatim("ValidationMisalignedTracker") << "\tPT: " << track->pt();
0682         edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:d0: " << track->d0();
0683         edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:z0: " << track->dz();
0684         edm::LogVerbatim("ValidationMisalignedTracker")
0685             << "\tAzimuthal angle of point of closest approach:" << track->phi();
0686         edm::LogVerbatim("ValidationMisalignedTracker") << "\tcharge: " << track->charge();
0687         edm::LogVerbatim("ValidationMisalignedTracker") << "\teta: " << track->eta();
0688         edm::LogVerbatim("ValidationMisalignedTracker") << "\tnormalizedChi2: " << track->normalizedChi2();
0689 
0690         std::vector<std::pair<TrackingParticleRef, double> > tp;
0691 
0692         //Compute fake rate vs eta
0693         if (recSimColl.find(track) != recSimColl.end()) {
0694           tp = recSimColl[track];
0695           if (!tp.empty()) {
0696             edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0697                                                << " associated with quality:" << tp.begin()->second << "\n";
0698 
0699             TrackingParticleRef tpr = tp.begin()->first;
0700             const SimTrack* fakeassocTrack = &(*tpr->g4Track_begin());
0701 
0702             const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0703             FreeTrajectoryState ftsAtProduction(
0704                 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0705                 GlobalVector(
0706                     fakeassocTrack->momentum().x(), fakeassocTrack->momentum().y(), fakeassocTrack->momentum().z()),
0707                 TrackCharge(tpr->charge()),
0708                 theMF);
0709             TSCPBuilderNoMaterial tscpBuilder;
0710             TrajectoryStateClosestToPoint tsAtClosestApproach =
0711                 tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
0712             GlobalPoint v = tsAtClosestApproach.theState().position();
0713             GlobalVector p = tsAtClosestApproach.theState().momentum();
0714 
0715             //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
0716             //  double lambdaSim = M_PI/2-p.theta();
0717             //  double phiSim    = p.phi();
0718             double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0719             double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0720             faked0 = float(-dxySim);
0721             fakez0 = float(dszSim * p.mag() / p.perp());
0722 
0723             faketrackType = fakeassocTrack->type();
0724             faketheta = fakeassocTrack->momentum().theta();
0725             fakecottheta = 1. / tan(faketheta);
0726             fakeeta = fakeassocTrack->momentum().eta();
0727             fakephi = fakeassocTrack->momentum().phi();
0728             fakept = fakeassocTrack->momentum().pt();
0729             fakenhit = tpr->matchedHit();
0730 
0731             edm::LogVerbatim("ValidationMisalignedTracker")
0732                 << "4) After call to associator: the best SimTrack match is of type" << fakeassocTrack->type()
0733                 << " ,at eta = " << fakeeta << " and phi = " << fakephi << " ,with pt at vertex = " << fakept
0734                 << " GeV/c"
0735                 << " ,d0 global = " << faked0 << " ,z0 = " << fakez0;
0736             fake = 1;
0737 
0738             fakerespt = fakerecpt - fakept;
0739             fakeresd0 = fakerecd0 - faked0;
0740             fakeresz0 = fakerecz0 - fakez0;
0741             fakereseta = -log(tan(fakerectheta / 2.)) - (-log(tan(faketheta / 2.)));
0742             fakeresphi = fakerecphi - fakephi;
0743             fakerescottheta = fakereccottheta - fakecottheta;
0744           }
0745         } else {
0746           edm::LogVerbatim("TrackValidator")
0747               << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any TrackingParticle"
0748               << "\n";
0749 
0750           fakeeta = -100.;
0751           faketheta = -100;
0752           fakephi = -100.;
0753           fakept = -100.;
0754           faked0 = -100.;
0755           fakez0 = -100;
0756           fakerespt = -100.;
0757           fakeresd0 = -100.;
0758           fakeresz0 = -100.;
0759           fakeresphi = -100.;
0760           fakereseta = -100.;
0761           fakerescottheta = -100.;
0762           fakenhit = 100;
0763           fake = 0;
0764         }
0765 
0766         tree_fake->Fill();
0767       }
0768 
0769     }  // End of loop on fakerate
0770 
0771     w++;
0772 
0773   }  // End of loop on associators
0774 }
0775 
0776 // ------------ method called once each job just after ending the event loop  ------------
0777 void ValidationMisalignedTracker::endJob() {
0778   edm::LogVerbatim("ValidationMisalignedTracker") << "\t Misalignment analysis completed \n";
0779 }
0780 
0781 DEFINE_FWK_MODULE(ValidationMisalignedTracker);
[ 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!