CMS 3D CMS Logo

ElectronStudy.cc
Go to the documentation of this file.
3 
16 
24 
29 
30 #include <TH1F.h>
31 
32 #include <algorithm>
33 #include <fstream>
34 #include <iostream>
35 #include <memory>
36 #include <string>
37 #include <vector>
38 
39 class ElectronStudy : public edm::one::EDAnalyzer<edm::one::SharedResources> {
40 public:
41  explicit ElectronStudy(const edm::ParameterSet& ps);
42  ~ElectronStudy() override {}
43 
44  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
45 
46 private:
47  void analyze(edm::Event const&, edm::EventSetup const&) override;
48  void beginJob() override {}
49 
50  static const int NEtaBins_ = 3;
51  static const int NPBins_ = 8;
52  double pBins_[NPBins_ + 1], etaBins_[NEtaBins_ + 1];
53 
58 
60  bool histos_;
61  TH1F* histoR1[NPBins_ + 1][NEtaBins_ + 1];
62  TH1F* histoR2[NPBins_ + 1][NEtaBins_ + 1];
63  TH1F* histoR3[NPBins_ + 1][NEtaBins_ + 1];
64  TH1F* histoE1x1[NPBins_ + 1][NEtaBins_ + 1];
65  TH1F* histoE3x3[NPBins_ + 1][NEtaBins_ + 1];
66  TH1F* histoE5x5[NPBins_ + 1][NEtaBins_ + 1];
67 };
68 
70  usesResource("TFileService");
71 
72  std::string g4Label = ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits");
73  std::string hitLabEB = ps.getUntrackedParameter<std::string>("EBCollection", "EcalHitsEB");
74  std::string hitLabEE = ps.getUntrackedParameter<std::string>("EECollection", "EcalHitsEE");
75 
76  tok_EBhit_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLabEB));
77  tok_EEhit_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLabEE));
78  tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag(g4Label));
79  tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag(g4Label));
80 
81  hotZone_ = ps.getUntrackedParameter<int>("HotZone", 0);
82  verbose_ = ps.getUntrackedParameter<int>("Verbosity", 0);
83  edm::LogInfo("ElectronStudy") << "Module Label: " << g4Label << " Hits: " << hitLabEB << ", " << hitLabEE;
84 
85  double tempP[NPBins_ + 1] = {0.0, 10.0, 20.0, 40.0, 60.0, 100.0, 500.0, 1000.0, 10000.0};
86  double tempEta[NEtaBins_ + 1] = {0.0, 1.2, 1.6, 3.0};
87 
88  for (int i = 0; i < NPBins_ + 1; i++)
89  pBins_[i] = tempP[i];
90  for (int i = 0; i < NEtaBins_ + 1; i++)
91  etaBins_[i] = tempEta[i];
92 
94  if (!tfile.isAvailable()) {
95  edm::LogInfo("ElectronStudy") << "TFileService unavailable: no histograms";
96  histos_ = false;
97  } else {
98  char name[20], title[200], cpbin[30], cebin[30];
99  histos_ = true;
100  for (unsigned int i = 0; i < NPBins_ + 1; ++i) {
101  if (i == 0)
102  sprintf(cpbin, " All p");
103  else
104  sprintf(cpbin, " p (%6.0f:%6.0f)", pBins_[i - 1], pBins_[i]);
105  for (unsigned int j = 0; j < NEtaBins_ + 1; ++j) {
106  if (j == 0)
107  sprintf(cebin, " All #eta");
108  else
109  sprintf(cebin, " #eta (%4.1f:%4.1f)", etaBins_[j - 1], etaBins_[j]);
110  sprintf(name, "R1%d%d", i, j);
111  sprintf(title, "E1/E9 for %s%s", cpbin, cebin);
112  histoR1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
113  histoR1[i][j]->GetXaxis()->SetTitle(title);
114  histoR1[i][j]->GetYaxis()->SetTitle("Tracks");
115  sprintf(name, "R2%d%d", i, j);
116  sprintf(title, "E1/E25 for %s%s", cpbin, cebin);
117  histoR2[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
118  histoR2[i][j]->GetXaxis()->SetTitle(title);
119  histoR2[i][j]->GetYaxis()->SetTitle("Tracks");
120  sprintf(name, "R3%d%d", i, j);
121  sprintf(title, "E9/E25 for %s%s", cpbin, cebin);
122  histoR3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
123  histoR3[i][j]->GetXaxis()->SetTitle(title);
124  histoR3[i][j]->GetYaxis()->SetTitle("Tracks");
125  sprintf(name, "E1x1%d%d", i, j);
126  sprintf(title, "E1/P for %s%s", cpbin, cebin);
127  histoE1x1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
128  histoE1x1[i][j]->GetXaxis()->SetTitle(title);
129  histoE1x1[i][j]->GetYaxis()->SetTitle("Tracks");
130  sprintf(name, "E3x3%d%d", i, j);
131  sprintf(title, "E9/P for %s%s", cpbin, cebin);
132  histoE3x3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
133  histoE3x3[i][j]->GetXaxis()->SetTitle(title);
134  histoE3x3[i][j]->GetYaxis()->SetTitle("Tracks");
135  sprintf(name, "E5x5%d%d", i, j);
136  sprintf(title, "E25/P for %s%s", cpbin, cebin);
137  histoE5x5[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
138  histoE5x5[i][j]->GetXaxis()->SetTitle(title);
139  histoE5x5[i][j]->GetYaxis()->SetTitle("Tracks");
140  }
141  }
142  }
143 }
144 
147  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
148  desc.addUntracked<std::string>("EBCollection", "EcalHitsEB");
149  desc.addUntracked<std::string>("EECollection", "EcalHitsEE");
150  desc.addUntracked<int>("Verbosity", 0);
151  descriptions.add("electronStudy", desc);
152 }
153 
155  if (verbose_ > 1)
156  edm::LogVerbatim("IsoTrack") << "Run = " << iEvent.id().run() << " Event = " << iEvent.id().event();
157 
158  // get Geometry, B-field, Topology
160  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
161  const MagneticField* bField = bFieldH.product();
162 
164  iSetup.get<CaloGeometryRecord>().get(pG);
165  const CaloGeometry* geo = pG.product();
166 
167  edm::ESHandle<CaloTopology> theCaloTopology;
168  iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
169  const CaloTopology* caloTopology = theCaloTopology.product();
170 
171  // get PCaloHits for ecal barrel
173  iEvent.getByToken(tok_EBhit_, caloHitEB);
174 
175  // get PCaloHits for ecal endcap
177  iEvent.getByToken(tok_EEhit_, caloHitEE);
178 
179  // get sim tracks
181  iEvent.getByToken(tok_simTk_, SimTk);
182 
183  // get sim vertices
185  iEvent.getByToken(tok_simVtx_, SimVtx);
186 
187  if (verbose_ > 0)
188  edm::LogVerbatim("IsoTrack") << "ElectronStudy: hits valid[EB]: " << caloHitEB.isValid()
189  << " valid[EE]: " << caloHitEE.isValid();
190 
191  if (caloHitEB.isValid() && caloHitEE.isValid()) {
192  unsigned int indx;
193  if (verbose_ > 2) {
194  edm::PCaloHitContainer::const_iterator ihit;
195  for (ihit = caloHitEB->begin(), indx = 0; ihit != caloHitEB->end(); ihit++, indx++) {
196  EBDetId id = ihit->id();
197  edm::LogVerbatim("IsoTrack") << "Hit[" << indx << "] " << id << " E " << ihit->energy() << " T "
198  << ihit->time();
199  }
200  for (ihit = caloHitEE->begin(), indx = 0; ihit != caloHitEE->end(); ihit++, indx++) {
201  EEDetId id = ihit->id();
202  edm::LogVerbatim("IsoTrack") << "Hit[" << indx << "] " << id << " E " << ihit->energy() << " T "
203  << ihit->time();
204  }
205  }
206  edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin();
207  for (indx = 0; simTrkItr != SimTk->end(); simTrkItr++, indx++) {
208  if (verbose_ > 0)
209  edm::LogVerbatim("IsoTrack") << "ElectronStudy: Track[" << indx << "] ID " << simTrkItr->trackId() << " type "
210  << simTrkItr->type() << " charge " << simTrkItr->charge() << " p "
211  << simTrkItr->momentum() << " Generator Index " << simTrkItr->genpartIndex()
212  << " vertex " << simTrkItr->vertIndex();
213  if (std::abs(simTrkItr->type()) == 11 && simTrkItr->vertIndex() != -1) {
214  int thisTrk = simTrkItr->trackId();
215  spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, (verbose_ > 1));
216  if (trkD.okECAL) {
217  const DetId isoCell = trkD.detIdECAL;
218  DetId hotCell = isoCell;
219  if (hotZone_ > 0)
220  hotCell = spr::hotCrystal(
221  isoCell, caloHitEB, caloHitEE, geo, caloTopology, hotZone_, hotZone_, -500.0, 500.0, (verbose_ > 1));
222  double e1x1 = spr::eECALmatrix(
223  hotCell, caloHitEB, caloHitEE, geo, caloTopology, 0, 0, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
224  double e3x3 = spr::eECALmatrix(
225  hotCell, caloHitEB, caloHitEE, geo, caloTopology, 1, 1, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
226  double e5x5 = spr::eECALmatrix(
227  hotCell, caloHitEB, caloHitEE, geo, caloTopology, 2, 2, -100.0, -100.0, -500.0, 500.0, (verbose_ > 2));
228  double p = simTrkItr->momentum().P();
229  double eta = std::abs(simTrkItr->momentum().eta());
230  int etaBin = -1, momBin = -1;
231  for (int ieta = 0; ieta < NEtaBins_; ieta++) {
232  if (eta > etaBins_[ieta] && eta < etaBins_[ieta + 1])
233  etaBin = ieta + 1;
234  }
235  for (int ipt = 0; ipt < NPBins_; ipt++) {
236  if (p > pBins_[ipt] && p < pBins_[ipt + 1])
237  momBin = ipt + 1;
238  }
239  double r1 = -1, r2 = -1, r3 = -1;
240  if (e3x3 > 0)
241  r1 = e1x1 / e3x3;
242  if (e5x5 > 0) {
243  r2 = e1x1 / e5x5;
244  r3 = e3x3 / e5x5;
245  }
246  if (verbose_ > 0) {
247  edm::LogVerbatim("IsoTrack") << "ElectronStudy: p " << p << " [" << momBin << "] eta " << eta << " ["
248  << etaBin << "] Cell 0x" << std::hex << isoCell() << std::dec;
249  if (isoCell.subdetId() == EcalBarrel) {
250  edm::LogVerbatim("IsoTrack") << EBDetId(isoCell);
251  } else if (isoCell.subdetId() == EcalEndcap) {
252  edm::LogVerbatim("IsoTrack") << EEDetId(isoCell);
253  }
254  edm::LogVerbatim("IsoTrack") << " e1x1 " << e1x1 << "|" << r1 << "|" << r2 << " e3x3 " << e3x3 << "|" << r3
255  << " e5x5 " << e5x5;
256  }
257  if (histos_) {
258  histoR1[0][0]->Fill(r1);
259  histoR2[0][0]->Fill(r2);
260  histoR3[0][0]->Fill(r3);
261  histoE1x1[0][0]->Fill(e1x1 / p);
262  histoE3x3[0][0]->Fill(e3x3 / p);
263  histoE5x5[0][0]->Fill(e5x5 / p);
264  if (momBin > 0) {
265  histoR1[momBin][0]->Fill(r1);
266  histoR2[momBin][0]->Fill(r2);
267  histoR3[momBin][0]->Fill(r3);
268  histoE1x1[momBin][0]->Fill(e1x1 / p);
269  histoE3x3[momBin][0]->Fill(e3x3 / p);
270  histoE5x5[momBin][0]->Fill(e5x5 / p);
271  }
272  if (etaBin > 0) {
273  histoR1[0][etaBin]->Fill(r1);
274  histoR2[0][etaBin]->Fill(r2);
275  histoR3[0][etaBin]->Fill(r3);
276  histoE1x1[0][etaBin]->Fill(e1x1 / p);
277  histoE3x3[0][etaBin]->Fill(e3x3 / p);
278  histoE5x5[0][etaBin]->Fill(e5x5 / p);
279  if (momBin > 0) {
280  histoR1[momBin][etaBin]->Fill(r1);
281  histoR2[momBin][etaBin]->Fill(r2);
282  histoR3[momBin][etaBin]->Fill(r3);
283  histoE1x1[momBin][etaBin]->Fill(e1x1 / p);
284  histoE3x3[momBin][etaBin]->Fill(e3x3 / p);
285  histoE5x5[momBin][etaBin]->Fill(e5x5 / p);
286  }
287  }
288  }
289  }
290  }
291  }
292  }
293 }
294 
295 //define this as a plug-in
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
EDAnalyzer.h
mps_fire.i
i
Definition: mps_fire.py:428
spr::hotCrystal
DetId hotCrystal(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double tMin=-500, double tMax=500, bool debug=false)
MessageLogger.h
ESHandle.h
etaBin
int etaBin(const l1t::HGCalMulticluster *cl)
Definition: L1EGammaEEProducer.cc:19
edm::EDGetTokenT< edm::PCaloHitContainer >
CaloGeometryRecord
Definition: CaloGeometryRecord.h:30
EBDetId
Definition: EBDetId.h:17
ElectronStudy::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: ElectronStudy.cc:145
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
ElectronStudy::ElectronStudy
ElectronStudy(const edm::ParameterSet &ps)
Definition: ElectronStudy.cc:69
ElectronStudy::tok_simVtx_
edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
Definition: ElectronStudy.cc:57
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
CaloTopologyRecord
Definition: CaloTopologyRecord.h:10
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
edm::one::EDAnalyzer
Definition: EDAnalyzer.h:30
edm::Handle
Definition: AssociativeIterator.h:50
ElectronStudy::histoE5x5
TH1F * histoE5x5[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:66
CaloTopology
Definition: CaloTopology.h:19
EcalBarrel
Definition: EcalSubdetector.h:10
ElectronStudy::NPBins_
static const int NPBins_
Definition: ElectronStudy.cc:51
IdealMagneticFieldRecord
Definition: IdealMagneticFieldRecord.h:11
DetId
Definition: DetId.h:17
MakerMacros.h
CaloGeometry
Definition: CaloGeometry.h:21
edm::EventSetup::get
T get() const
Definition: EventSetup.h:80
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
spr::eECALmatrix
double eECALmatrix(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
ElectronStudy::pBins_
double pBins_[NPBins_+1]
Definition: ElectronStudy.cc:52
spr::propagatedTrackDirection::okECAL
bool okECAL
Definition: CaloPropagateTrack.h:59
Service.h
PVValHelper::eta
Definition: PVValidationHelpers.h:69
ElectronStudy::hotZone_
int hotZone_
Definition: ElectronStudy.cc:59
ElectronStudy::histos_
bool histos_
Definition: ElectronStudy.cc:60
tfile
Definition: tfile.py:1
IdealMagneticFieldRecord.h
edm::ESHandle< MagneticField >
ElectronStudy::verbose_
int verbose_
Definition: ElectronStudy.cc:59
ElectronStudy::tok_simTk_
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
Definition: ElectronStudy.cc:56
eECALMatrix.h
ElectronStudy::etaBins_
double etaBins_[NEtaBins_+1]
Definition: ElectronStudy.cc:52
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum)
Definition: DetId.h:48
EEDetId
Definition: EEDetId.h:14
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
EcalEndcap
Definition: EcalSubdetector.h:10
TFileService.h
ElectronStudy::NEtaBins_
static const int NEtaBins_
Definition: ElectronStudy.cc:50
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
CaloSubdetectorGeometry.h
ElectronStudy
Definition: ElectronStudy.cc:39
edm::ParameterSet
Definition: ParameterSet.h:47
ElectronStudy::histoR2
TH1F * histoR2[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:62
Event.h
ElectronStudy::histoR3
TH1F * histoR3[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:63
PCaloHit.h
CaloTopologyRecord.h
diffTwoXMLs.r2
r2
Definition: diffTwoXMLs.py:73
edm::Service< TFileService >
iEvent
int iEvent
Definition: GenABIO.cc:224
MagneticField.h
edm::EventSetup
Definition: EventSetup.h:57
get
#define get
InputTag.h
compare.tfile
tfile
Definition: compare.py:325
CaloTopology.h
spr::propagateCALO
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
Definition: CaloPropagateTrack.cc:141
CaloSubdetectorTopology.h
ElectronStudy::~ElectronStudy
~ElectronStudy() override
Definition: ElectronStudy.cc:42
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
ElectronStudy::histoE3x3
TH1F * histoE3x3[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:65
Calorimetry_cff.bField
bField
Definition: Calorimetry_cff.py:292
Frameworkfwd.h
diffTwoXMLs.r1
r1
Definition: diffTwoXMLs.py:53
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
CaloPropagateTrack.h
ElectronStudy::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: ElectronStudy.cc:154
ElectronStudy::tok_EEhit_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_EEhit_
Definition: ElectronStudy.cc:55
CaloGeometry.h
spr::propagatedTrackDirection::detIdECAL
DetId detIdECAL
Definition: CaloPropagateTrack.h:60
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
ElectronStudy::tok_EBhit_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_EBhit_
Definition: ElectronStudy.cc:54
Exception.h
PCaloHitContainer.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
edm::Event
Definition: Event.h:73
MagneticField
Definition: MagneticField.h:19
ElectronStudy::histoE1x1
TH1F * histoE1x1[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:64
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
SimTrackContainer.h
edm::InputTag
Definition: InputTag.h:15
ElectronStudy::histoR1
TH1F * histoR1[NPBins_+1][NEtaBins_+1]
Definition: ElectronStudy.cc:61
spr::propagatedTrackDirection
Definition: CaloPropagateTrack.h:53
SimVertexContainer.h
ElectronStudy::beginJob
void beginJob() override
Definition: ElectronStudy.cc:48