Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57:11 

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   for (unsigned int ww = 0; ww < associators.size(); ww++) {
0347     //
0348     //get collections from the event
0349     //
0350 
0351     const edm::Handle<edm::View<reco::Track> >& trackCollection = iEvent.getHandle(trackToken_);
0352     const edm::View<reco::Track> tC = *(trackCollection.product());
0353 
0354     //associate tracks
0355     LogTrace("TrackValidator") << "Calling associateRecoToSim method\n";
0356     reco::RecoToSimCollection recSimColl = associatore[ww]->associateRecoToSim(trackCollection, TPCollectionHfake);
0357 
0358     LogTrace("TrackValidator") << "Calling associateSimToReco method\n";
0359     reco::SimToRecoCollection simRecColl = associatore[ww]->associateSimToReco(trackCollection, TPCollectionHeff);
0360 
0361     //
0362     //compute number of tracks per eta interval
0363     //
0364 
0365     if (selection_eff && !skip) {
0366       edm::LogVerbatim("ValidationMisalignedTracker") << "Computing Efficiency";
0367 
0368       edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles (before cuts): " << tPCeff.size() << "\n";
0369       int st = 0;
0370       for (TrackingParticleCollection::size_type i = 0; i < tPCeff.size(); i++) {
0371         // Initialize variables
0372         eta = 0., theta = 0., phi = 0., pt = 0., cottheta = 0., costheta = 0.;
0373         d0 = 0., z0 = 0.;
0374         nhit = 0;
0375         receta = 0., rectheta = 0., recphi = 0., recpt = 0., reccottheta = 0., recd0 = 0., recz0 = 0.;
0376         respt = 0., resd0 = 0., resz0 = 0., reseta = 0., resphi = 0., rescottheta = 0.;
0377         recchiq = 0.;
0378         recnhit = 0;
0379         trackType = 0;
0380         eff = 0;
0381 
0382         // typedef edm::Ref<TrackingParticleCollection> TrackingParticleRef;
0383         TrackingParticleRef tp(TPCollectionHeff, i);
0384         if (tp->charge() == 0)
0385           continue;
0386         st++;
0387         //pt=sqrt(tp->momentum().perp2());
0388         //eta=tp->momentum().eta();
0389         //vpos=tp->vertex().perp2()));
0390 
0391         const SimTrack* simulatedTrack = &(*tp->g4Track_begin());
0392 
0393         const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0394         FreeTrajectoryState ftsAtProduction(
0395             GlobalPoint(tp->vertex().x(), tp->vertex().y(), tp->vertex().z()),
0396             GlobalVector(
0397                 simulatedTrack->momentum().x(), simulatedTrack->momentum().y(), simulatedTrack->momentum().z()),
0398             TrackCharge(tp->charge()),
0399             theMF);
0400         TSCPBuilderNoMaterial tscpBuilder;
0401         TrajectoryStateClosestToPoint tsAtClosestApproach =
0402             tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
0403         GlobalPoint v = tsAtClosestApproach.theState().position();
0404         GlobalVector p = tsAtClosestApproach.theState().momentum();
0405 
0406         //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
0407         //  double lambdaSim = M_PI/2-p.theta();
0408         //  double phiSim    = p.phi();
0409         double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0410         double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0411         d0 = float(-dxySim);
0412         z0 = float(dszSim * p.mag() / p.perp());
0413 
0414         if (abs(simulatedTrack->type()) == 13 && simulatedTrack->genpartIndex() != -1) {
0415           edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA SIM DI MUONI ";
0416           edm::LogVerbatim("ValidationMisalignedTracker") << "Gen part " << simulatedTrack->genpartIndex();
0417           trackType = simulatedTrack->type();
0418           theta = simulatedTrack->momentum().theta();
0419           costheta = cos(theta);
0420           cottheta = 1. / tan(theta);
0421 
0422           eta = simulatedTrack->momentum().eta();
0423           phi = simulatedTrack->momentum().phi();
0424           pt = simulatedTrack->momentum().pt();
0425           nhit = tp->matchedHit();
0426 
0427           edm::LogVerbatim("ValidationMisalignedTracker")
0428               << "3) Before assoc: SimTrack of type = " << simulatedTrack->type() << " ,at eta = " << eta
0429               << " ,with pt at vertex = " << simulatedTrack->momentum().pt() << " GeV/c"
0430               << " ,d0 =" << d0 << " ,z0 =" << z0 << " ,nhit=" << nhit;
0431 
0432           if (ZmassSelection_) {
0433             if (abs(trackType) == 13 &&
0434                 (simulatedTrack->genpartIndex() == indmu[0] || simulatedTrack->genpartIndex() == indmu[1])) {
0435               edm::LogVerbatim("ValidationMisalignedTracker") << " TRACK sim of muons from Z ";
0436               flag = 0;
0437               count = countpart[0];
0438               countpart[0]++;
0439             } else if (abs(trackType) == 11) {
0440               //edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA SIM DI ELETTRONI ";
0441               flag = 1;
0442               count = countpart[1];
0443               countpart[1]++;
0444             }
0445 
0446             px[flag][count] = simulatedTrack->momentum().x();
0447             py[flag][count] = simulatedTrack->momentum().y();
0448             pz[flag][count] = simulatedTrack->momentum().z();
0449             ptmu[flag][count] = simulatedTrack->momentum().pt();
0450             ene[flag][count] = simulatedTrack->momentum().e();
0451           }
0452 
0453           std::vector<std::pair<edm::RefToBase<reco::Track>, double> > rt;
0454           if (simRecColl.find(tp) != simRecColl.end()) {
0455             rt = simRecColl[tp];
0456             if (!rt.empty()) {
0457               edm::RefToBase<reco::Track> t = rt.begin()->first;
0458 
0459               // bool flagptused=false;
0460               // for (unsigned int j=0;j<ptused.size();j++){
0461               //   if (fabs(t->pt()-ptused[j])<0.001) {
0462               //     flagptused=true;
0463               //   }
0464               // }
0465 
0466               edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st << " with pt=" << t->pt()
0467                                                  << " associated with quality:" << rt.begin()->second << "\n";
0468               edm::LogVerbatim("ValidationMisalignedTracker") << "Reconstructed Track:" << t->pt();
0469               edm::LogVerbatim("ValidationMisalignedTracker") << "\tpT: " << t->pt();
0470               edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:d0: " << t->d0();
0471               edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:z0: " << t->dz();
0472               edm::LogVerbatim("ValidationMisalignedTracker")
0473                   << "\tAzimuthal angle of point of closest approach:" << t->phi();
0474               edm::LogVerbatim("ValidationMisalignedTracker") << "\tcharge: " << t->charge();
0475               edm::LogVerbatim("ValidationMisalignedTracker") << "\teta: " << t->eta();
0476               edm::LogVerbatim("ValidationMisalignedTracker") << "\tnormalizedChi2: " << t->normalizedChi2();
0477 
0478               recnhit = t->numberOfValidHits();
0479               recchiq = t->normalizedChi2();
0480               rectheta = t->theta();
0481               reccottheta = 1. / tan(rectheta);
0482               //receta=-log(tan(rectheta/2.));
0483               receta = t->momentum().eta();
0484               //       reccostheta=cos(matchedrectrack->momentum().theta());
0485               recphi = t->phi();
0486               recpt = t->pt();
0487               ptused.push_back(recpt);
0488               recd0 = t->d0();
0489               recz0 = t->dz();
0490 
0491               edm::LogVerbatim("ValidationMisalignedTracker")
0492                   << "5) After call to associator: the best match has " << recnhit << " hits, Chi2 = " << recchiq
0493                   << ", pt at vertex = " << recpt << " GeV/c, "
0494                   << ", recd0 = " << recd0 << ", recz0= " << recz0;
0495 
0496               respt = recpt - pt;
0497               resd0 = recd0 - d0;
0498               resz0 = recz0 - z0;
0499               reseta = receta - eta;
0500               resphi = recphi - phi;
0501               rescottheta = reccottheta - cottheta;
0502               eff = 1;
0503 
0504               edm::LogVerbatim("ValidationMisalignedTracker")
0505                   << "6) Transverse momentum residual=" << respt << " ,d0 residual=" << resd0
0506                   << " ,z0 residual=" << resz0 << " with eff=" << eff;
0507 
0508               if (ZmassSelection_) {
0509                 if (abs(trackType) == 13) {
0510                   edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA RECO DI MUONI ";
0511                   flagrec = 0;
0512                   countrec = countpartrec[0];
0513                   countpartrec[0]++;
0514                 } else if (abs(trackType) == 11) {
0515                   edm::LogVerbatim("ValidationMisalignedTracker") << " TRACCIA RECO DI ELETTRONI ";
0516                   flagrec = 1;
0517                   countrec = countpartrec[1];
0518                   countpartrec[1]++;
0519                 }
0520 
0521                 recp[flagrec][countrec] = sqrt(t->momentum().mag2());
0522                 recpx[flagrec][countrec] = t->momentum().x();
0523                 recpy[flagrec][countrec] = t->momentum().y();
0524                 recpz[flagrec][countrec] = t->momentum().z();
0525                 recptmu[flagrec][countrec] =
0526                     sqrt((t->momentum().x() * t->momentum().x()) + (t->momentum().y() * t->momentum().y()));
0527                 if (abs(trackType) == 13)
0528                   recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.105 * 0.105);
0529                 if (abs(trackType) == 11)
0530                   recene[flagrec][countrec] = sqrt(recp[flagrec][countrec] * recp[flagrec][countrec] + 0.0005 * 0.0005);
0531               }
0532 
0533               edm::LogVerbatim("ValidationMisalignedTracker") << "7) Transverse momentum reconstructed =" << recpt
0534                                                               << " at  eta= " << receta << " and phi= " << recphi;
0535             }
0536           } else {
0537             edm::LogVerbatim("TrackValidator")
0538                 << "TrackingParticle #" << st << " with pt=" << sqrt(tp->momentum().perp2())
0539                 << " NOT associated to any reco::Track"
0540                 << "\n";
0541             receta = -100.;
0542             recphi = -100.;
0543             recpt = -100.;
0544             recd0 = -100.;
0545             recz0 = -100;
0546             respt = -100.;
0547             resd0 = -100.;
0548             resz0 = -100.;
0549             resphi = -100.;
0550             reseta = -100.;
0551             rescottheta = -100.;
0552             recnhit = 100;
0553             recchiq = -100;
0554             eff = 0;
0555             flagrec = 100;
0556           }
0557 
0558           edm::LogVerbatim("ValidationMisalignedTracker") << "Eff=" << eff;
0559 
0560           // simulated muons
0561 
0562           edm::LogVerbatim("ValidationMisalignedTracker") << "Flag is" << flag;
0563           edm::LogVerbatim("ValidationMisalignedTracker") << "RecFlag is" << flagrec;
0564 
0565           if (countpart[0] == 2 && flag == 0) {
0566             mzmu =
0567                 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]) -
0568                      (py[0][0] + py[0][1]) * (py[0][0] + py[0][1]) - (pz[0][0] + pz[0][1]) * (pz[0][0] + pz[0][1]));
0569             edm::LogVerbatim("ValidationMisalignedTracker") << "Mzmu " << mzmu;
0570             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]));
0571 
0572             pLzmu = pz[0][0] + pz[0][1];
0573             enezmu = ene[0][0] + ene[0][1];
0574             phizmu = atan2((py[0][0] + py[0][1]), (px[0][0] + px[0][1]));
0575             thetazmu = atan2(ptzmu, (pz[0][0] + pz[0][1]));
0576             etazmu = -log(tan(thetazmu * 3.14 / 360.));
0577             yzmu = 0.5 * log((enezmu + pLzmu) / (enezmu - pLzmu));
0578             mxptmu = std::max(ptmu[0][0], ptmu[0][1]);
0579             minptmu = std::min(ptmu[0][0], ptmu[0][1]);
0580           } else {
0581             mzmu = -100.;
0582             ptzmu = -100.;
0583             pLzmu = -100.;
0584             enezmu = -100.;
0585             etazmu = -100.;
0586             phizmu = -100.;
0587             thetazmu = -100.;
0588             yzmu = -100.;
0589             mxptmu = -100.;
0590             minptmu = -100.;
0591           }
0592 
0593           // reconstructed muons
0594           if (countpartrec[0] == 2 && flagrec == 0) {
0595             recmzmu = sqrt((recene[0][0] + recene[0][1]) * (recene[0][0] + recene[0][1]) -
0596                            (recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) -
0597                            (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]) -
0598                            (recpz[0][0] + recpz[0][1]) * (recpz[0][0] + recpz[0][1]));
0599             edm::LogVerbatim("ValidationMisalignedTracker") << "RecMzmu " << recmzmu;
0600             recptzmu = sqrt((recpx[0][0] + recpx[0][1]) * (recpx[0][0] + recpx[0][1]) +
0601                             (recpy[0][0] + recpy[0][1]) * (recpy[0][0] + recpy[0][1]));
0602 
0603             recpLzmu = recpz[0][0] + recpz[0][1];
0604             recenezmu = recene[0][0] + recene[0][1];
0605             recphizmu = atan2((recpy[0][0] + recpy[0][1]), (recpx[0][0] + recpx[0][1]));
0606             recthetazmu = atan2(recptzmu, (recpz[0][0] + recpz[0][1]));
0607             recetazmu = -log(tan(recthetazmu * 3.14 / 360.));
0608             recyzmu = 0.5 * log((recenezmu + recpLzmu) / (recenezmu - recpLzmu));
0609             recmxptmu = std::max(recptmu[0][0], recptmu[0][1]);
0610             recminptmu = std::min(recptmu[0][0], recptmu[0][1]);
0611           } else {
0612             recmzmu = -100.;
0613             recptzmu = -100.;
0614             recpLzmu = -100.;
0615             recenezmu = -100.;
0616             recetazmu = -100.;
0617             recphizmu = -100.;
0618             recthetazmu = -100.;
0619             recyzmu = -100.;
0620             recmxptmu = -100;
0621             recminptmu = -100.;
0622           }
0623 
0624           tree_eff->Fill();
0625 
0626         }  // end of loop on muons
0627       }    // end of loop for tracking particle
0628     }      // end of loop for efficiency
0629 
0630     //
0631     // Fake Rate
0632     //
0633     if (selection_fake) {
0634       edm::LogVerbatim("ValidationMisalignedTracker") << "Computing Fake Rate";
0635 
0636       fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0637       faked0 = 0., fakez0 = 0.;
0638       fakenhit = 0;
0639       fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0640       fakerecz0 = 0.;
0641       fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0642       fakerecchiq = 0.;
0643       fakerecnhit = 0;
0644       faketrackType = 0;
0645       fake = 0;
0646 
0647       //      int at=0;
0648       int rT = 0;
0649       for (reco::TrackCollection::size_type i = 0; i < tC.size(); ++i) {
0650         edm::RefToBase<reco::Track> track(trackCollection, i);
0651         rT++;
0652 
0653         fakeeta = 0., faketheta = 0., fakephi = 0., fakept = 0., fakecottheta = 0., fakecostheta = 0.;
0654         faked0 = 0., fakez0 = 0.;
0655         fakenhit = 0;
0656         fakereceta = 0., fakerectheta = 0., fakerecphi = 0., fakerecpt = 0., fakereccottheta = 0., fakerecd0 = 0.,
0657         fakerecz0 = 0.;
0658         fakerespt = 0., fakeresd0 = 0., fakeresz0 = 0., fakereseta = 0., fakeresphi = 0., fakerescottheta = 0.;
0659         fakerecchiq = 0.;
0660         fakerecnhit = 0;
0661         faketrackType = 0;
0662         fake = 0;
0663 
0664         fakerecnhit = track->numberOfValidHits();
0665         fakerecchiq = track->normalizedChi2();
0666         fakerectheta = track->theta();
0667         fakereccottheta = 1. / tan(rectheta);
0668         //fakereceta=-log(tan(rectheta/2.));
0669         fakereceta = track->momentum().eta();
0670         //     fakereccostheta=cos(track->momentum().theta());
0671         fakerecphi = track->phi();
0672         fakerecpt = track->pt();
0673         fakerecd0 = track->d0();
0674         fakerecz0 = track->dz();
0675 
0676         edm::LogVerbatim("ValidationMisalignedTracker") << "1) Before assoc: TkRecTrack at eta = " << fakereceta;
0677         edm::LogVerbatim("ValidationMisalignedTracker") << "Track number " << i;
0678         edm::LogVerbatim("ValidationMisalignedTracker") << "\tPT: " << track->pt();
0679         edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:d0: " << track->d0();
0680         edm::LogVerbatim("ValidationMisalignedTracker") << "\timpact parameter:z0: " << track->dz();
0681         edm::LogVerbatim("ValidationMisalignedTracker")
0682             << "\tAzimuthal angle of point of closest approach:" << track->phi();
0683         edm::LogVerbatim("ValidationMisalignedTracker") << "\tcharge: " << track->charge();
0684         edm::LogVerbatim("ValidationMisalignedTracker") << "\teta: " << track->eta();
0685         edm::LogVerbatim("ValidationMisalignedTracker") << "\tnormalizedChi2: " << track->normalizedChi2();
0686 
0687         std::vector<std::pair<TrackingParticleRef, double> > tp;
0688 
0689         //Compute fake rate vs eta
0690         if (recSimColl.find(track) != recSimColl.end()) {
0691           tp = recSimColl[track];
0692           if (!tp.empty()) {
0693             edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0694                                                << " associated with quality:" << tp.begin()->second << "\n";
0695 
0696             TrackingParticleRef tpr = tp.begin()->first;
0697             const SimTrack* fakeassocTrack = &(*tpr->g4Track_begin());
0698 
0699             const MagneticField* theMF = &iSetup.getData(magFieldToken_);
0700             FreeTrajectoryState ftsAtProduction(
0701                 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0702                 GlobalVector(
0703                     fakeassocTrack->momentum().x(), fakeassocTrack->momentum().y(), fakeassocTrack->momentum().z()),
0704                 TrackCharge(tpr->charge()),
0705                 theMF);
0706             TSCPBuilderNoMaterial tscpBuilder;
0707             TrajectoryStateClosestToPoint tsAtClosestApproach =
0708                 tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
0709             GlobalPoint v = tsAtClosestApproach.theState().position();
0710             GlobalVector p = tsAtClosestApproach.theState().momentum();
0711 
0712             //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
0713             //  double lambdaSim = M_PI/2-p.theta();
0714             //  double phiSim    = p.phi();
0715             double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0716             double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0717             faked0 = float(-dxySim);
0718             fakez0 = float(dszSim * p.mag() / p.perp());
0719 
0720             faketrackType = fakeassocTrack->type();
0721             faketheta = fakeassocTrack->momentum().theta();
0722             fakecottheta = 1. / tan(faketheta);
0723             fakeeta = fakeassocTrack->momentum().eta();
0724             fakephi = fakeassocTrack->momentum().phi();
0725             fakept = fakeassocTrack->momentum().pt();
0726             fakenhit = tpr->matchedHit();
0727 
0728             edm::LogVerbatim("ValidationMisalignedTracker")
0729                 << "4) After call to associator: the best SimTrack match is of type" << fakeassocTrack->type()
0730                 << " ,at eta = " << fakeeta << " and phi = " << fakephi << " ,with pt at vertex = " << fakept
0731                 << " GeV/c"
0732                 << " ,d0 global = " << faked0 << " ,z0 = " << fakez0;
0733             fake = 1;
0734 
0735             fakerespt = fakerecpt - fakept;
0736             fakeresd0 = fakerecd0 - faked0;
0737             fakeresz0 = fakerecz0 - fakez0;
0738             fakereseta = -log(tan(fakerectheta / 2.)) - (-log(tan(faketheta / 2.)));
0739             fakeresphi = fakerecphi - fakephi;
0740             fakerescottheta = fakereccottheta - fakecottheta;
0741           }
0742         } else {
0743           edm::LogVerbatim("TrackValidator")
0744               << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any TrackingParticle"
0745               << "\n";
0746 
0747           fakeeta = -100.;
0748           faketheta = -100;
0749           fakephi = -100.;
0750           fakept = -100.;
0751           faked0 = -100.;
0752           fakez0 = -100;
0753           fakerespt = -100.;
0754           fakeresd0 = -100.;
0755           fakeresz0 = -100.;
0756           fakeresphi = -100.;
0757           fakereseta = -100.;
0758           fakerescottheta = -100.;
0759           fakenhit = 100;
0760           fake = 0;
0761         }
0762 
0763         tree_fake->Fill();
0764       }
0765 
0766     }  // End of loop on fakerate
0767 
0768   }  // End of loop on associators
0769 }
0770 
0771 // ------------ method called once each job just after ending the event loop  ------------
0772 void ValidationMisalignedTracker::endJob() {
0773   edm::LogVerbatim("ValidationMisalignedTracker") << "\t Misalignment analysis completed \n";
0774 }
0775 
0776 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!