IsolatedParticles/test/TestIsoTracks.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestIsoTracks
0004 // Class:      IsolatedParticles
0005 //
0006 /*
0007
0008
0009  Description: <one line class summary>
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  Sergey Petrushanko
0016 //
0017
0018 // system include files
0019 #include <memory>
0020
0021 // user include files
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0024
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "DataFormats/Common/interface/Handle.h"
0029 #include "DataFormats/Common/interface/OrphanHandle.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033
0034 #include "DataFormats/TrackReco/interface/Track.h"
0035 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0036 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0037 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0038 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0039 //#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0040 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0041
0042 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0043 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0044 #include "SimDataFormats/Track/interface/SimTrack.h"
0045 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0046 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0047 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0049 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0050
0051 // calorimeter info
0052 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0053 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0054 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0055 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0056 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0057 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0058
0059 //#include "Geometry/DTGeometry/interface/DTLayer.h"
0060 //#include "Geometry/DTGeometry/interface/DTGeometry.h"
0061 //#include "Geometry/Records/interface/MuonGeometryRecord.h"
0062
0063 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0064 #include "DataFormats/GeometrySurface/interface/Plane.h"
0065
0066 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0067
0068 #include "MagneticField/Engine/interface/MagneticField.h"
0069 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0070
0071 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0072
0073 #include <boost/regex.hpp>
0074
0075 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0076
0077 #include "TH1F.h"
0078 #include <TFile.h>
0079
0080 #include "RecoTracker/TrackProducer/interface/TrackProducerBase.h"
0081 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0082 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0083 #include "DataFormats/VertexReco/interface/Vertex.h"
0084 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0085
0086 #include <map>
0087 #include <string>
0088 #include <vector>
0089 #include "CLHEP/Vector/LorentzVector.h"
0090
0091 using namespace std;
0092 using namespace reco;
0093
0094 class TestIsoTracks : public edm::one::EDAnalyzer<> {
0095 public:
0096   explicit TestIsoTracks(const edm::ParameterSet&);
0097
0098   void setPrimaryVertex(const reco::Vertex& a) { theRecVertex = a; }
0099   void setTracksFromPrimaryVertex(vector<reco::Track>& a) { theTrack = a; }
0100   void analyze(const edm::Event&, const edm::EventSetup&) override;
0101   void endJob() override;
0102
0103 private:
0104   TFile* m_Hfile;
0105   struct {
0106     TH1F* eta;
0107     TH1F* phi;
0108     TH1F* p;
0109     TH1F* pt;
0110     TH1F* isomult;
0111   } IsoHists;
0112
0113   edm::EDGetTokenT<reco::VertexCollection> mInputPVfCTF;
0114   edm::EDGetTokenT<reco::TrackCollection> m_inputTrackToken;
0115
0116   double theRvert;
0117   reco::Vertex theRecVertex;
0118   vector<reco::Track> theTrack;
0119   TrackDetectorAssociator trackAssociator_;
0120   TrackAssociatorParameters trackAssociatorParameters_;
0121
0122   edm::EDGetTokenT<edm::SimTrackContainer> simTracksToken_;
0123   edm::EDGetTokenT<edm::SimVertexContainer> simVerticesToken_;
0124 };
0125
0126 TestIsoTracks::TestIsoTracks(const edm::ParameterSet& iConfig)
0127     : mInputPVfCTF(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("src3"))),
0128       m_inputTrackToken(consumes<reco::TrackCollection>(
0129           edm::InputTag(iConfig.getUntrackedParameter<std::string>("inputTrackLabel", "ctfWithMaterialTracks")))),
0130       theRvert(iConfig.getParameter<double>("rvert")),
0131       simTracksToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simTracksTag"))),
0132       simVerticesToken_(consumes<edm::SimVertexContainer>(iConfig.getParameter<edm::InputTag>("simVerticesTag"))) {
0133   m_Hfile = new TFile("IsoHists.root", "RECREATE");
0134   IsoHists.eta = new TH1F("Eta", "Track eta", 50, -2.5, 2.5);
0135   IsoHists.phi = new TH1F("Phi", "Track phi", 50, -3.2, 3.2);
0136   IsoHists.p = new TH1F("Momentum", "Track momentum", 100, 0., 20.);
0137   IsoHists.pt = new TH1F("pt", "Track pt", 100, 0., 10.);
0138   IsoHists.isomult = new TH1F("IsoMult", "Iso Mult", 10, -0.5, 9.5);
0139
0140   // Load TrackDetectorAssociator parameters
0141   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0142   edm::ConsumesCollector iC = consumesCollector();
0143   trackAssociatorParameters_.loadParameters(parameters, iC);
0144   trackAssociator_.useDefaultPropagator();
0145 }
0146
0147 void TestIsoTracks::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0148   using namespace edm;
0149
0150   // mine! b
0151
0152   std::vector<GlobalPoint> AllTracks;
0153   std::vector<GlobalPoint> AllTracks1;
0154
0155   // mine! e
0156
0157   // Take Reco Vertex collection
0158   edm::Handle<reco::VertexCollection> primary_vertices;  //Define Inputs (vertices)
0159   iEvent.getByToken(mInputPVfCTF, primary_vertices);     //Get Inputs    (vertices)
0160
0161   // Take Reco Track Collection
0162   edm::Handle<reco::TrackCollection> trackCollection;
0163   iEvent.getByToken(m_inputTrackToken, trackCollection);
0164   const reco::TrackCollection tC = *(trackCollection.product());
0165
0166   // get list of tracks and their vertices
0167   Handle<SimTrackContainer> simTracks;
0168   iEvent.getByToken(simTracksToken_, simTracks);
0169
0170   Handle<SimVertexContainer> simVertices;
0171   iEvent.getByToken(simVerticesToken_, simVertices);
0172   if (!simVertices.isValid())
0173     throw cms::Exception("FatalError") << "No vertices found\n";
0174
0175   // loop over tracks
0176
0177   std::cout << "Number of tracks found in the event: " << trackCollection->size() << std::endl;
0178
0179   //   for(SimTrackContainer::const_iterator tracksCI = simTracks->begin();
0180   //       tracksCI != simTracks->end(); tracksCI++){
0181
0182   for (reco::TrackCollection::const_iterator tracksCI = tC.begin(); tracksCI != tC.end(); tracksCI++) {
0183     double mome_pt =
0184         sqrt(tracksCI->momentum().x() * tracksCI->momentum().x() + tracksCI->momentum().y() * tracksCI->momentum().y());
0185
0186     // skip low Pt tracks
0187     if (mome_pt < 0.5) {
0188       std::cout << "Skipped low Pt track (Pt: " << mome_pt << ")" << std::endl;
0189       continue;
0190     }
0191
0192     // get vertex
0193     //      int vertexIndex = tracksCI->vertIndex();
0194     // uint trackIndex = tracksCI->genpartIndex();
0195
0196     //      SimVertex vertex(Hep3Vector(0.,0.,0.),0);
0197     //      if (vertexIndex >= 0) vertex = (*simVertices)[vertexIndex];
0198
0199     // skip tracks originated away from the IP
0200     //      if (vertex.position().rho() > 50) {
0201     //   std::cout << "Skipped track originated away from IP: " <<vertex.position().rho()<<std::endl;
0202     //   continue;
0203     //      }
0204
0205     std::cout << "\n-------------------------------------------------------\n Track (pt,eta,phi): " << mome_pt << " , "
0206               << tracksCI->momentum().eta() << " , " << tracksCI->momentum().phi() << std::endl;
0207
0208     // Simply get ECAL energy of the crossed crystals
0209     //      std::cout << "ECAL energy of crossed crystals: " <<
0210     //  trackAssociator_.getEcalEnergy(iEvent, iSetup,
0211     //                     trackAssociator_.getFreeTrajectoryState(iSetup, *tracksCI, vertex) )
0212     //    << " GeV" << std::endl;
0213
0214     //      std::cout << "Details:\n" <<std::endl;
0215     TrackDetMatchInfo info = trackAssociator_.associate(
0216         iEvent,
0217         iSetup,
0218         trackAssociator_.getFreeTrajectoryState(&iSetup.getData(trackAssociatorParameters_.bFieldToken), *tracksCI),
0219         trackAssociatorParameters_);
0220     //      std::cout << "ECAL, if track reach ECAL:     " << info.isGoodEcal << std::endl;
0221     //      std::cout << "ECAL, number of crossed cells: " << info.crossedEcalRecHits.size() << std::endl;
0222     //      std::cout << "ECAL, energy of crossed cells: " << info.ecalEnergy() << " GeV" << std::endl;
0223     //      std::cout << "ECAL, number of cells in the cone: " << info.ecalRecHits.size() << std::endl;
0224     //      std::cout << "ECAL, energy in the cone: " << info.ecalConeEnergy() << " GeV" << std::endl;
0225     //      std::cout << "ECAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtEcal.z() << ", "
0226     //  << info.trkGlobPosAtEcal.R() << " , "   << info.trkGlobPosAtEcal.eta() << " , "
0227     //  << info.trkGlobPosAtEcal.phi()<< std::endl;
0228
0229     // mine! b
0230
0231     double rfa =
0232         sqrt(info.trkGlobPosAtEcal.x() * info.trkGlobPosAtEcal.x() +
0233              info.trkGlobPosAtEcal.y() * info.trkGlobPosAtEcal.y() +
0234              info.trkGlobPosAtEcal.z() * info.trkGlobPosAtEcal.z()) /
0235         sqrt(tracksCI->momentum().x() * tracksCI->momentum().x() + tracksCI->momentum().y() * tracksCI->momentum().y() +
0236              tracksCI->momentum().z() * tracksCI->momentum().z());
0237
0238     if (info.isGoodEcal == 1 && fabs(info.trkGlobPosAtEcal.eta()) < 2.6) {
0239       AllTracks.push_back(GlobalPoint(
0240           info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0241       if (mome_pt > 2. && fabs(info.trkGlobPosAtEcal.eta()) < 2.1)
0242         if (fabs(info.trkGlobPosAtEcal.eta()) < 2.1) {
0243           AllTracks1.push_back(GlobalPoint(
0244               info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0245         }
0246     }
0247
0248     // mine! e
0249
0250     //      std::cout << "HCAL, if track reach HCAL:      " << info.isGoodHcal << std::endl;
0251     //      std::cout << "HCAL, number of crossed towers: " << info.crossedTowers.size() << std::endl;
0252     //      std::cout << "HCAL, energy of crossed towers: " << info.hcalEnergy() << " GeV" << std::endl;
0253     //      std::cout << "HCAL, number of towers in the cone: " << info.towers.size() << std::endl;
0254     //      std::cout << "HCAL, energy in the cone: " << info.hcalConeEnergy() << " GeV" << std::endl;
0255     //      std::cout << "HCAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtHcal.z() << ", "
0256     //  << info.trkGlobPosAtHcal.R() << " , "   << info.trkGlobPosAtHcal.eta() << " , "
0257     //  << info.trkGlobPosAtHcal.phi()<< std::endl;
0258   }
0259
0260   // mine! b
0261
0262   std::cout << " NUMBER of tracks  " << AllTracks.size() << "  and candidates for iso tracks  " << AllTracks1.size()
0263             << std::endl;
0264
0265   double imult = 0.;
0266
0267   for (unsigned int ia1 = 0; ia1 < AllTracks1.size(); ia1++) {
0268     double delta_min = 3.141592;
0269
0270     for (unsigned int ia = 0; ia < AllTracks.size(); ia++) {
0271       double delta_phi = fabs(AllTracks1[ia1].phi() - AllTracks[ia].phi());
0272       if (delta_phi > 3.141592)
0273         delta_phi = 6.283184 - delta_phi;
0274       double delta_eta = fabs(AllTracks1[ia1].eta() - AllTracks[ia].eta());
0275       double delta_actual = sqrt(delta_phi * delta_phi + delta_eta * delta_eta);
0276
0277       if (delta_actual < delta_min && delta_actual != 0.)
0278         delta_min = delta_actual;
0279     }
0280
0281     if (delta_min > 0.5) {
0282       std::cout << "FIND ISOLATED TRACK " << AllTracks1[ia1].mag() << "  " << AllTracks1[ia1].eta() << "  "
0283                 << AllTracks1[ia1].phi() << std::endl;
0284
0285       IsoHists.eta->Fill(AllTracks1[ia1].eta());
0286       IsoHists.phi->Fill(AllTracks1[ia1].phi());
0287       IsoHists.p->Fill(AllTracks1[ia1].mag());
0288       IsoHists.pt->Fill(AllTracks1[ia1].perp());
0289
0290       imult = imult + 1.;
0291     }
0292   }
0293
0294   IsoHists.isomult->Fill(imult);
0295
0296   // mine! e
0297 }
0298
0299 void TestIsoTracks::endJob(void) {
0300   m_Hfile->cd();
0301   IsoHists.eta->Write();
0302   IsoHists.phi->Write();
0303   IsoHists.p->Write();
0304   IsoHists.pt->Write();
0305   IsoHists.isomult->Write();
0306   m_Hfile->Close();
0307 }
0308
0309 //define this as a plug-in
0310 DEFINE_FWK_MODULE(TestIsoTracks);