Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:34

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/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::EDAnalyzer {
0095 public:
0096   explicit TestIsoTracks(const edm::ParameterSet&);
0097   virtual ~TestIsoTracks(){};
0098 
0099   void setPrimaryVertex(const reco::Vertex& a) { theRecVertex = a; }
0100   void setTracksFromPrimaryVertex(vector<reco::Track>& a) { theTrack = a; }
0101   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0102   void endJob(void);
0103 
0104 private:
0105   TFile* m_Hfile;
0106   struct {
0107     TH1F* eta;
0108     TH1F* phi;
0109     TH1F* p;
0110     TH1F* pt;
0111     TH1F* isomult;
0112   } IsoHists;
0113 
0114   edm::EDGetTokenT<reco::VertexCollection> mInputPVfCTF;
0115   edm::EDGetTokenT<reco::TrackCollection> m_inputTrackToken;
0116 
0117   double theRvert;
0118   reco::Vertex theRecVertex;
0119   vector<reco::Track> theTrack;
0120   TrackDetectorAssociator trackAssociator_;
0121   TrackAssociatorParameters trackAssociatorParameters_;
0122 
0123   edm::EDGetTokenT<edm::SimTrackContainer> simTracksToken_;
0124   edm::EDGetTokenT<edm::SimVertexContainer> simVerticesToken_;
0125 };
0126 
0127 TestIsoTracks::TestIsoTracks(const edm::ParameterSet& iConfig)
0128     : mInputPVfCTF(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("src3"))),
0129       m_inputTrackToken(consumes<reco::TrackCollection>(
0130           edm::InputTag(iConfig.getUntrackedParameter<std::string>("inputTrackLabel", "ctfWithMaterialTracks")))),
0131       theRvert(iConfig.getParameter<double>("rvert")),
0132       simTracksToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simTracksTag"))),
0133       simVerticesToken_(consumes<edm::SimVertexContainer>(iConfig.getParameter<edm::InputTag>("simVerticesTag"))) {
0134   m_Hfile = new TFile("IsoHists.root", "RECREATE");
0135   IsoHists.eta = new TH1F("Eta", "Track eta", 50, -2.5, 2.5);
0136   IsoHists.phi = new TH1F("Phi", "Track phi", 50, -3.2, 3.2);
0137   IsoHists.p = new TH1F("Momentum", "Track momentum", 100, 0., 20.);
0138   IsoHists.pt = new TH1F("pt", "Track pt", 100, 0., 10.);
0139   IsoHists.isomult = new TH1F("IsoMult", "Iso Mult", 10, -0.5, 9.5);
0140 
0141   // Load TrackDetectorAssociator parameters
0142   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0143   edm::ConsumesCollector iC = consumesCollector();
0144   trackAssociatorParameters_.loadParameters(parameters, iC);
0145   trackAssociator_.useDefaultPropagator();
0146 }
0147 
0148 void TestIsoTracks::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0149   using namespace edm;
0150 
0151   // mine! b
0152 
0153   std::vector<GlobalPoint> AllTracks;
0154   std::vector<GlobalPoint> AllTracks1;
0155 
0156   // mine! e
0157 
0158   // Take Reco Vertex collection
0159   edm::Handle<reco::VertexCollection> primary_vertices;  //Define Inputs (vertices)
0160   iEvent.getByToken(mInputPVfCTF, primary_vertices);     //Get Inputs    (vertices)
0161 
0162   // Take Reco Track Collection
0163   edm::Handle<reco::TrackCollection> trackCollection;
0164   iEvent.getByToken(m_inputTrackToken, trackCollection);
0165   const reco::TrackCollection tC = *(trackCollection.product());
0166 
0167   // get list of tracks and their vertices
0168   Handle<SimTrackContainer> simTracks;
0169   iEvent.getByToken(simTracksToken_, simTracks);
0170 
0171   Handle<SimVertexContainer> simVertices;
0172   iEvent.getByToken(simVerticesToken_, simVertices);
0173   if (!simVertices.isValid())
0174     throw cms::Exception("FatalError") << "No vertices found\n";
0175 
0176   // loop over tracks
0177 
0178   std::cout << "Number of tracks found in the event: " << trackCollection->size() << std::endl;
0179 
0180   //   for(SimTrackContainer::const_iterator tracksCI = simTracks->begin();
0181   //       tracksCI != simTracks->end(); tracksCI++){
0182 
0183   for (reco::TrackCollection::const_iterator tracksCI = tC.begin(); tracksCI != tC.end(); tracksCI++) {
0184     double mome_pt =
0185         sqrt(tracksCI->momentum().x() * tracksCI->momentum().x() + tracksCI->momentum().y() * tracksCI->momentum().y());
0186 
0187     // skip low Pt tracks
0188     if (mome_pt < 0.5) {
0189       std::cout << "Skipped low Pt track (Pt: " << mome_pt << ")" << std::endl;
0190       continue;
0191     }
0192 
0193     // get vertex
0194     //      int vertexIndex = tracksCI->vertIndex();
0195     // uint trackIndex = tracksCI->genpartIndex();
0196 
0197     //      SimVertex vertex(Hep3Vector(0.,0.,0.),0);
0198     //      if (vertexIndex >= 0) vertex = (*simVertices)[vertexIndex];
0199 
0200     // skip tracks originated away from the IP
0201     //      if (vertex.position().rho() > 50) {
0202     //   std::cout << "Skipped track originated away from IP: " <<vertex.position().rho()<<std::endl;
0203     //   continue;
0204     //      }
0205 
0206     std::cout << "\n-------------------------------------------------------\n Track (pt,eta,phi): " << mome_pt << " , "
0207               << tracksCI->momentum().eta() << " , " << tracksCI->momentum().phi() << std::endl;
0208 
0209     // Simply get ECAL energy of the crossed crystals
0210     //      std::cout << "ECAL energy of crossed crystals: " <<
0211     //  trackAssociator_.getEcalEnergy(iEvent, iSetup,
0212     //                     trackAssociator_.getFreeTrajectoryState(iSetup, *tracksCI, vertex) )
0213     //    << " GeV" << std::endl;
0214 
0215     //      std::cout << "Details:\n" <<std::endl;
0216     TrackDetMatchInfo info = trackAssociator_.associate(
0217         iEvent,
0218         iSetup,
0219         trackAssociator_.getFreeTrajectoryState(&iSetup.getData(trackAssociatorParameters_.bFieldToken), *tracksCI),
0220         trackAssociatorParameters_);
0221     //      std::cout << "ECAL, if track reach ECAL:     " << info.isGoodEcal << std::endl;
0222     //      std::cout << "ECAL, number of crossed cells: " << info.crossedEcalRecHits.size() << std::endl;
0223     //      std::cout << "ECAL, energy of crossed cells: " << info.ecalEnergy() << " GeV" << std::endl;
0224     //      std::cout << "ECAL, number of cells in the cone: " << info.ecalRecHits.size() << std::endl;
0225     //      std::cout << "ECAL, energy in the cone: " << info.ecalConeEnergy() << " GeV" << std::endl;
0226     //      std::cout << "ECAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtEcal.z() << ", "
0227     //  << info.trkGlobPosAtEcal.R() << " , "   << info.trkGlobPosAtEcal.eta() << " , "
0228     //  << info.trkGlobPosAtEcal.phi()<< std::endl;
0229 
0230     // mine! b
0231 
0232     double rfa =
0233         sqrt(info.trkGlobPosAtEcal.x() * info.trkGlobPosAtEcal.x() +
0234              info.trkGlobPosAtEcal.y() * info.trkGlobPosAtEcal.y() +
0235              info.trkGlobPosAtEcal.z() * info.trkGlobPosAtEcal.z()) /
0236         sqrt(tracksCI->momentum().x() * tracksCI->momentum().x() + tracksCI->momentum().y() * tracksCI->momentum().y() +
0237              tracksCI->momentum().z() * tracksCI->momentum().z());
0238 
0239     if (info.isGoodEcal == 1 && fabs(info.trkGlobPosAtEcal.eta()) < 2.6) {
0240       AllTracks.push_back(GlobalPoint(
0241           info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0242       if (mome_pt > 2. && fabs(info.trkGlobPosAtEcal.eta()) < 2.1)
0243         if (fabs(info.trkGlobPosAtEcal.eta()) < 2.1) {
0244           AllTracks1.push_back(GlobalPoint(
0245               info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0246         }
0247     }
0248 
0249     // mine! e
0250 
0251     //      std::cout << "HCAL, if track reach HCAL:      " << info.isGoodHcal << std::endl;
0252     //      std::cout << "HCAL, number of crossed towers: " << info.crossedTowers.size() << std::endl;
0253     //      std::cout << "HCAL, energy of crossed towers: " << info.hcalEnergy() << " GeV" << std::endl;
0254     //      std::cout << "HCAL, number of towers in the cone: " << info.towers.size() << std::endl;
0255     //      std::cout << "HCAL, energy in the cone: " << info.hcalConeEnergy() << " GeV" << std::endl;
0256     //      std::cout << "HCAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtHcal.z() << ", "
0257     //  << info.trkGlobPosAtHcal.R() << " , "   << info.trkGlobPosAtHcal.eta() << " , "
0258     //  << info.trkGlobPosAtHcal.phi()<< std::endl;
0259   }
0260 
0261   // mine! b
0262 
0263   std::cout << " NUMBER of tracks  " << AllTracks.size() << "  and candidates for iso tracks  " << AllTracks1.size()
0264             << std::endl;
0265 
0266   double imult = 0.;
0267 
0268   for (unsigned int ia1 = 0; ia1 < AllTracks1.size(); ia1++) {
0269     double delta_min = 3.141592;
0270 
0271     for (unsigned int ia = 0; ia < AllTracks.size(); ia++) {
0272       double delta_phi = fabs(AllTracks1[ia1].phi() - AllTracks[ia].phi());
0273       if (delta_phi > 3.141592)
0274         delta_phi = 6.283184 - delta_phi;
0275       double delta_eta = fabs(AllTracks1[ia1].eta() - AllTracks[ia].eta());
0276       double delta_actual = sqrt(delta_phi * delta_phi + delta_eta * delta_eta);
0277 
0278       if (delta_actual < delta_min && delta_actual != 0.)
0279         delta_min = delta_actual;
0280     }
0281 
0282     if (delta_min > 0.5) {
0283       std::cout << "FIND ISOLATED TRACK " << AllTracks1[ia1].mag() << "  " << AllTracks1[ia1].eta() << "  "
0284                 << AllTracks1[ia1].phi() << std::endl;
0285 
0286       IsoHists.eta->Fill(AllTracks1[ia1].eta());
0287       IsoHists.phi->Fill(AllTracks1[ia1].phi());
0288       IsoHists.p->Fill(AllTracks1[ia1].mag());
0289       IsoHists.pt->Fill(AllTracks1[ia1].perp());
0290 
0291       imult = imult + 1.;
0292     }
0293   }
0294 
0295   IsoHists.isomult->Fill(imult);
0296 
0297   // mine! e
0298 }
0299 
0300 void TestIsoTracks::endJob(void) {
0301   m_Hfile->cd();
0302   IsoHists.eta->Write();
0303   IsoHists.phi->Write();
0304   IsoHists.p->Write();
0305   IsoHists.pt->Write();
0306   IsoHists.isomult->Write();
0307   m_Hfile->Close();
0308 }
0309 
0310 //define this as a plug-in
0311 DEFINE_FWK_MODULE(TestIsoTracks);