Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestIsoSimTracks
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 class TestIsoSimTracks : public edm::EDAnalyzer {
0081 public:
0082   explicit TestIsoSimTracks(const edm::ParameterSet&);
0083   virtual ~TestIsoSimTracks(){};
0084 
0085   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0086   void endJob(void);
0087 
0088 private:
0089   TFile* m_Hfile;
0090   struct {
0091     TH1F* eta;
0092     TH1F* phi;
0093     TH1F* p;
0094     TH1F* pt;
0095     TH1F* isomult;
0096   } IsoHists;
0097   TrackDetectorAssociator trackAssociator_;
0098   TrackAssociatorParameters trackAssociatorParameters_;
0099 
0100   edm::EDGetTokenT<edm::SimTrackContainer> simTracksToken_;
0101   edm::EDGetTokenT<edm::SimVertexContainer> simVerticesToken_;
0102 };
0103 
0104 TestIsoSimTracks::TestIsoSimTracks(const edm::ParameterSet& iConfig)
0105     : simTracksToken_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simTracksTag"))),
0106       simVerticesToken_(consumes<edm::SimVertexContainer>(iConfig.getParameter<edm::InputTag>("simVerticesTag"))) {
0107   // Fill data labels
0108   //std::vector<std::string> labels = iConfig.getParameter<std::vector<std::string> >("labels");
0109   //boost::regex regExp1 ("([^\\s,]+)[\\s,]+([^\\s,]+)$");
0110   //boost::regex regExp2 ("([^\\s,]+)[\\s,]+([^\\s,]+)[\\s,]+([^\\s,]+)$");
0111   //boost::smatch matches;
0112 
0113   m_Hfile = new TFile("IsoHists.root", "RECREATE");
0114   IsoHists.eta = new TH1F("Eta", "Track eta", 100, -5., 5.);
0115   IsoHists.phi = new TH1F("Phi", "Track phi", 100, -3.5, 3.5);
0116   IsoHists.p = new TH1F("Momentum", "Track momentum", 100, 0., 20.);
0117   IsoHists.pt = new TH1F("pt", "Track pt", 100, 0., 10.);
0118   IsoHists.isomult = new TH1F("IsoMult", "Iso Mult", 10, -0.5, 9.5);
0119 
0120   //for(std::vector<std::string>::const_iterator label = labels.begin(); label != labels.end(); label++) {
0121   //   if (boost::regex_match(*label,matches,regExp1))
0122   //    trackAssociator_.addDataLabels(matches[1],matches[2]);
0123   //     else if (boost::regex_match(*label,matches,regExp2))
0124   //    trackAssociator_.addDataLabels(matches[1],matches[2],matches[3]);
0125   //     else
0126   //    edm::LogError("ConfigurationError") << "Failed to parse label:\n" << *label << "Skipped.\n";
0127   //  }
0128 
0129   // trackAssociator_.addDataLabels("EBRecHitCollection","ecalrechit","EcalRecHitsEB");
0130   // trackAssociator_.addDataLabels("CaloTowerCollection","towermaker");
0131   // trackAssociator_.addDataLabels("DTRecSegment4DCollection","recseg4dbuilder");
0132 
0133   // Load TrackDetectorAssociator parameters
0134   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0135   edm::ConsumesCollector iC = consumesCollector();
0136   trackAssociatorParameters_.loadParameters(parameters, iC);
0137   trackAssociator_.useDefaultPropagator();
0138 }
0139 
0140 void TestIsoSimTracks::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0141   using namespace edm;
0142 
0143   // mine! b
0144 
0145   std::vector<GlobalPoint> AllTracks;
0146   std::vector<GlobalPoint> AllTracks1;
0147 
0148   // mine! e
0149 
0150   // get list of tracks and their vertices
0151   Handle<SimTrackContainer> simTracks;
0152   iEvent.getByToken(simTracksToken_, simTracks);
0153 
0154   Handle<SimVertexContainer> simVertices;
0155   iEvent.getByToken(simVerticesToken_, simVertices);
0156   if (!simVertices.isValid())
0157     throw cms::Exception("FatalError") << "No vertices found\n";
0158 
0159   // loop over simulated tracks
0160   std::cout << "Number of simulated tracks found in the event: " << simTracks->size() << std::endl;
0161   for (SimTrackContainer::const_iterator tracksCI = simTracks->begin(); tracksCI != simTracks->end(); tracksCI++) {
0162     // skip low Pt tracks
0163     if (tracksCI->momentum().Pt() < 0.7) {
0164       //     std::cout << "Skipped low Pt track (Pt: " << tracksCI->momentum().perp() << ")" <<std::endl;
0165       continue;
0166     }
0167 
0168     // get vertex
0169     int vertexIndex = tracksCI->vertIndex();
0170     // uint trackIndex = tracksCI->genpartIndex();
0171 
0172     SimVertex vertex(math::XYZVectorD(0., 0., 0.), 0);
0173     if (vertexIndex >= 0)
0174       vertex = (*simVertices)[vertexIndex];
0175 
0176     // skip tracks originated away from the IP
0177     //      if (vertex.position().rho() > 50) {
0178     //   std::cout << "Skipped track originated away from IP: " <<vertex.position().rho()<<std::endl;
0179     //   continue;
0180     //      }
0181 
0182     std::cout << "\n-------------------------------------------------------\n Track (pt,eta,phi): "
0183               << tracksCI->momentum().Pt() << " , " << tracksCI->momentum().eta() << " , " << tracksCI->momentum().phi()
0184               << std::endl;
0185 
0186     // Simply get ECAL energy of the crossed crystals
0187     //      std::cout << "ECAL energy of crossed crystals: " <<
0188     //  trackAssociator_.getEcalEnergy(iEvent, iSetup,
0189     //                     trackAssociator_.getFreeTrajectoryState(iSetup, *tracksCI, vertex) )
0190     //    << " GeV" << std::endl;
0191 
0192     //      std::cout << "Details:\n" <<std::endl;
0193     TrackDetMatchInfo info =
0194         trackAssociator_.associate(iEvent,
0195                                    iSetup,
0196                                    trackAssociator_.getFreeTrajectoryState(
0197                                        &iSetup.getData(trackAssociatorParameters_.bFieldToken), *tracksCI, vertex),
0198                                    trackAssociatorParameters_);
0199     //      std::cout << "ECAL, if track reach ECAL:     " << info.isGoodEcal << std::endl;
0200     //      std::cout << "ECAL, number of crossed cells: " << info.crossedEcalRecHits.size() << std::endl;
0201     //      std::cout << "ECAL, energy of crossed cells: " << info.ecalEnergy() << " GeV" << std::endl;
0202     //      std::cout << "ECAL, number of cells in the cone: " << info.ecalRecHits.size() << std::endl;
0203     //      std::cout << "ECAL, energy in the cone: " << info.ecalConeEnergy() << " GeV" << std::endl;
0204     //      std::cout << "ECAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtEcal.z() << ", "
0205     //  << info.trkGlobPosAtEcal.R() << " , "   << info.trkGlobPosAtEcal.eta() << " , "
0206     //  << info.trkGlobPosAtEcal.phi()<< std::endl;
0207 
0208     // mine! b
0209 
0210     double rfa =
0211         sqrt(info.trkGlobPosAtEcal.x() * info.trkGlobPosAtEcal.x() +
0212              info.trkGlobPosAtEcal.y() * info.trkGlobPosAtEcal.y() +
0213              info.trkGlobPosAtEcal.z() * info.trkGlobPosAtEcal.z()) /
0214         sqrt(tracksCI->momentum().x() * tracksCI->momentum().x() + tracksCI->momentum().y() * tracksCI->momentum().y() +
0215              tracksCI->momentum().z() * tracksCI->momentum().z());
0216 
0217     if (info.isGoodEcal == 1 && fabs(info.trkGlobPosAtEcal.eta()) < 2.6) {
0218       AllTracks.push_back(GlobalPoint(
0219           info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0220       if (tracksCI->momentum().Pt() > 2. && fabs(info.trkGlobPosAtEcal.eta()) < 2.1) {
0221         AllTracks1.push_back(GlobalPoint(
0222             info.trkGlobPosAtEcal.x() / rfa, info.trkGlobPosAtEcal.y() / rfa, info.trkGlobPosAtEcal.z() / rfa));
0223       }
0224     }
0225 
0226     // mine! e
0227 
0228     //      std::cout << "HCAL, if track reach HCAL:      " << info.isGoodHcal << std::endl;
0229     //      std::cout << "HCAL, number of crossed towers: " << info.crossedTowers.size() << std::endl;
0230     //      std::cout << "HCAL, energy of crossed towers: " << info.hcalEnergy() << " GeV" << std::endl;
0231     //      std::cout << "HCAL, number of towers in the cone: " << info.towers.size() << std::endl;
0232     //      std::cout << "HCAL, energy in the cone: " << info.hcalConeEnergy() << " GeV" << std::endl;
0233     //      std::cout << "HCAL, trajectory point (z,R,eta,phi): " << info.trkGlobPosAtHcal.z() << ", "
0234     //  << info.trkGlobPosAtHcal.R() << " , "   << info.trkGlobPosAtHcal.eta() << " , "
0235     //  << info.trkGlobPosAtHcal.phi()<< std::endl;
0236   }
0237 
0238   // mine! b
0239 
0240   std::cout << " NUMBER of tracks  " << AllTracks.size() << "  and candidates for iso tracks  " << AllTracks1.size()
0241             << std::endl;
0242 
0243   double imult = 0.;
0244 
0245   for (unsigned int ia1 = 0; ia1 < AllTracks1.size(); ia1++) {
0246     double delta_min = 3.141592;
0247 
0248     for (unsigned int ia = 0; ia < AllTracks.size(); ia++) {
0249       double delta_phi = fabs(AllTracks1[ia1].phi() - AllTracks[ia].phi());
0250       if (delta_phi > 3.141592)
0251         delta_phi = 6.283184 - delta_phi;
0252       double delta_eta = fabs(AllTracks1[ia1].eta() - AllTracks[ia].eta());
0253       double delta_actual = sqrt(delta_phi * delta_phi + delta_eta * delta_eta);
0254 
0255       if (delta_actual < delta_min && delta_actual != 0.)
0256         delta_min = delta_actual;
0257     }
0258 
0259     if (delta_min > 0.5) {
0260       std::cout << "FIND ISOLATED TRACK " << AllTracks1[ia1].mag() << "  " << AllTracks1[ia1].eta() << "  "
0261                 << AllTracks1[ia1].phi() << std::endl;
0262 
0263       IsoHists.eta->Fill(AllTracks1[ia1].eta());
0264       IsoHists.phi->Fill(AllTracks1[ia1].phi());
0265       IsoHists.p->Fill(AllTracks1[ia1].mag());
0266       IsoHists.pt->Fill(AllTracks1[ia1].perp());
0267       imult = imult + 1.;
0268     }
0269   }
0270   IsoHists.isomult->Fill(imult);
0271 
0272   // mine! e
0273 }
0274 
0275 void TestIsoSimTracks::endJob(void) {
0276   m_Hfile->cd();
0277   IsoHists.eta->Write();
0278   IsoHists.phi->Write();
0279   IsoHists.p->Write();
0280   IsoHists.pt->Write();
0281   IsoHists.isomult->Write();
0282   m_Hfile->Close();
0283 }
0284 
0285 //define this as a plug-in
0286 DEFINE_FWK_MODULE(TestIsoSimTracks);