IsolatedParticles/test/TestIsoSimTracks.cc

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