Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:45:58 

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