CMS 3D CMS Logo

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