File indexing completed on 2023-03-17 11:10:43
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020
0021
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
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
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
0060
0061
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
0107
0108
0109
0110
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
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
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
0143
0144 std::vector<GlobalPoint> AllTracks;
0145 std::vector<GlobalPoint> AllTracks1;
0146
0147
0148
0149
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
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
0162 if (tracksCI->momentum().Pt() < 0.7) {
0163
0164 continue;
0165 }
0166
0167
0168 int vertexIndex = tracksCI->vertIndex();
0169
0170
0171 SimVertex vertex(math::XYZVectorD(0., 0., 0.), 0);
0172 if (vertexIndex >= 0)
0173 vertex = (*simVertices)[vertexIndex];
0174
0175
0176
0177
0178
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
0186
0187
0188
0189
0190
0191
0192 TrackDetMatchInfo info =
0193 trackAssociator_.associate(iEvent,
0194 iSetup,
0195 trackAssociator_.getFreeTrajectoryState(
0196 &iSetup.getData(trackAssociatorParameters_.bFieldToken), *tracksCI, vertex),
0197 trackAssociatorParameters_);
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
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
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235 }
0236
0237
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
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
0285 DEFINE_FWK_MODULE(TestIsoSimTracks);