67 std::vector<edm::EDGetTokenT<edm::PCaloHitContainer> >
toks_calo_;
87 maxEnergy_(ps.getUntrackedParameter<double>(
"maxEnergy", 250.0)),
88 maxTime_(ps.getUntrackedParameter<double>(
"maxTime", 1000.0)),
89 tMax_(ps.getUntrackedParameter<double>(
"timeCut", 100.0)),
90 tScale_(ps.getUntrackedParameter<double>(
"timeScale", 1.0)),
91 tCut_(ps.getUntrackedParameter<double>(
"timeThreshold", 15.0)),
98 for (
unsigned int i = 0;
i <
hitLab_.size();
i++)
103 for (
unsigned int i = 0; i <
hitLab_.size(); i++)
113 <<
"please add it to config file";
116 for (
int i = 0; i <
nCalo_; i++) {
117 sprintf(name,
"Hit%d", i);
118 sprintf(title,
"Number of hits in %s", dets[i].c_str());
120 h_hit_[
i]->GetXaxis()->SetTitle(title);
121 h_hit_[
i]->GetYaxis()->SetTitle(
"Events");
122 sprintf(name,
"Time%d", i);
123 sprintf(title,
"Time of the hit (ns) in %s", dets[i].c_str());
125 h_time_[
i]->GetXaxis()->SetTitle(title);
126 h_time_[
i]->GetYaxis()->SetTitle(
"Hits");
127 sprintf(name,
"TimeT%d", i);
128 sprintf(title,
"Time of each hit (ns) in %s", dets[i].c_str());
130 h_timeT_[
i]->GetXaxis()->SetTitle(title);
131 h_timeT_[
i]->GetYaxis()->SetTitle(
"Hits");
132 double ymax = (i > 1) ? 0.01 : 0.1;
133 double ymx0 = (i > 1) ? 0.0025 : 0.025;
134 sprintf(name,
"Edep%d", i);
135 sprintf(title,
"Energy deposit (GeV) in %s", dets[i].c_str());
137 h_edep_[
i]->GetXaxis()->SetTitle(title);
138 h_edep_[
i]->GetYaxis()->SetTitle(
"Hits");
139 sprintf(name,
"EdepT%d", i);
140 sprintf(title,
"Energy deposit (GeV) of each hit in %s", dets[i].c_str());
142 h_edepT_[
i]->GetXaxis()->SetTitle(title);
143 h_edepT_[
i]->GetYaxis()->SetTitle(
"Hits");
144 sprintf(name,
"EdepEM%d", i);
145 sprintf(title,
"Energy deposit (GeV) by EM particles in %s", dets[i].c_str());
149 sprintf(name,
"EdepHad%d", i);
150 sprintf(title,
"Energy deposit (GeV) by hadrons in %s", dets[i].c_str());
154 sprintf(name,
"Edep15%d", i);
155 sprintf(title,
"Energy deposit (GeV) for T > %4.0f ns in %s",
tCut_, dets[i].c_str());
157 h_edep1_[
i]->GetXaxis()->SetTitle(title);
158 h_edep1_[
i]->GetYaxis()->SetTitle(
"Hits");
159 sprintf(name,
"EdepT15%d", i);
160 sprintf(title,
"Energy deposit (GeV) of each hit for T > %4.0f in %s",
tCut_, dets[i].c_str());
165 sprintf(name,
"Etot%d", i);
166 sprintf(title,
"Total energy deposit (GeV) in %s", dets[i].c_str());
168 h_etot_[
i]->GetXaxis()->SetTitle(title);
169 h_etot_[
i]->GetYaxis()->SetTitle(
"Events");
170 sprintf(name,
"EtotG%d", i);
171 sprintf(title,
"Total energy deposit (GeV) in %s (t < 100 ns)", dets[i].c_str());
173 h_etotg_[
i]->GetXaxis()->SetTitle(title);
174 h_etotg_[
i]->GetYaxis()->SetTitle(
"Events");
175 sprintf(name,
"rr%d", i);
176 sprintf(title,
"R of hit point (cm) in %s", dets[i].c_str());
178 h_rr_[
i]->GetXaxis()->SetTitle(title);
179 h_rr_[
i]->GetYaxis()->SetTitle(
"Hits");
180 sprintf(name,
"zz%d", i);
181 sprintf(title,
"z of hit point (cm) in %s", dets[i].c_str());
183 h_zz_[
i]->GetXaxis()->SetTitle(title);
184 h_zz_[
i]->GetYaxis()->SetTitle(
"Hits");
185 sprintf(name,
"eta%d", i);
186 sprintf(title,
"#eta of hit point in %s", dets[i].c_str());
188 h_eta_[
i]->GetXaxis()->SetTitle(title);
189 h_eta_[
i]->GetYaxis()->SetTitle(
"Hits");
190 sprintf(name,
"phi%d", i);
191 sprintf(title,
"#phi of hit point in %s", dets[i].c_str());
193 h_phi_[
i]->GetXaxis()->SetTitle(title);
194 h_phi_[
i]->GetYaxis()->SetTitle(
"Hits");
196 sprintf(title,
"R vs Z of hit point");
197 h_rz_ = tfile->
make<TH2F>(
"rz",
title, 120, 0., 600., 100, 0., 250.);
198 h_rz_->GetXaxis()->SetTitle(
"z (cm)");
199 h_rz_->GetYaxis()->SetTitle(
"R (cm)");
200 sprintf(title,
"R vs Z of hit point for hits with T > %4.0f ns",
tCut_);
201 h_rz1_ = tfile->
make<TH2F>(
"rz2",
title, 120, 0., 600., 100, 0., 250.);
202 h_rz1_->GetXaxis()->SetTitle(
"z (cm)");
203 h_rz1_->GetYaxis()->SetTitle(
"R (cm)");
204 sprintf(title,
"#phi vs #eta of hit point");
210 h_hitp_ = tfile->
make<TH1F>(
"hitp",
"All Steps", 100, 0.0, 20000.0);
211 h_hitp_->GetXaxis()->SetTitle(
"Hits");
212 h_hitp_->GetYaxis()->SetTitle(
"Events");
213 h_trackp_ = tfile->
make<TH1F>(
"trackp",
"All Steps", 100, 0.0, 200000.0);
214 h_hitp_->GetXaxis()->SetTitle(
"Tracks");
215 h_hitp_->GetYaxis()->SetTitle(
"Events");
216 h_edepp_ = tfile->
make<TH1F>(
"edepp",
"All Steps", 100, 0.0, 50.0);
217 h_edepp_->GetXaxis()->SetTitle(
"Energy Deposit (MeV)");
218 h_edepp_->GetYaxis()->SetTitle(
"Hits");
219 h_timep_ = tfile->
make<TH1F>(
"timep",
"All Steps", 100, 0.0, 100.0);
220 h_timep_->GetXaxis()->SetTitle(
"Hits");
221 h_timep_->GetYaxis()->SetTitle(
"Hit Time (ns)");
222 h_stepp_ = tfile->
make<TH1F>(
"stepp",
"All Steps", 1000, 0.0, 100.0);
223 h_stepp_->GetXaxis()->SetTitle(
"Hits");
224 h_stepp_->GetYaxis()->SetTitle(
"Step length (cm)");
225 for (
unsigned int k = 0;
k < detNames_.size(); ++
k) {
226 sprintf(name,
"edept%d",
k);
227 sprintf(title,
"Energy Deposit (MeV) in %s", detNames_[
k].c_str());
228 h_edepTk_.emplace_back(tfile->
make<TH1F>(name, title, 100, 0.0, 1.0));
229 h_edepTk_.back()->GetYaxis()->SetTitle(
"Hits");
230 h_edepTk_.back()->GetXaxis()->SetTitle(title);
231 sprintf(name,
"timet%d",
k);
232 sprintf(title,
"Hit Time (ns) in %s", detNames_[
k].c_str());
233 h_timeTk_.emplace_back(tfile->
make<TH1F>(name, title, 100, 0.0, 100.0));
234 h_timeTk_.back()->GetYaxis()->SetTitle(
"Hits");
235 h_timeTk_.back()->GetXaxis()->SetTitle(title);
242 std::vector<std::string>
labels = {
"EcalHitsEB1",
"EcalHitsEE1",
"HcalHits1"};
243 std::vector<double> times = {1, 1, 1};
245 desc.
add<std::vector<std::string> >(
"hitCollection",
labels);
246 desc.
add<std::vector<double> >(
"timeSliceUnit", times);
254 std::vector<std::string>
names = {
"PixelBarrel",
"PixelForward",
"TIB",
"TID",
"TOB",
"TEC"};
256 descriptions.
add(
"caloSimHitAnalysis", desc);
270 bool getHits = (hitsCalo.
isValid());
271 edm::LogVerbatim(
"HitStudy") <<
"CaloSimHitAnalysis: Input flags Hits[" << i <<
"]: " << getHits;
274 std::vector<PCaloHit> caloHits;
275 caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
276 edm::LogVerbatim(
"HitStudy") <<
"CaloSimHitAnalysis: Hit buffer [" << i <<
"] " << caloHits.size();
284 bool getHits = (hitsPassive.
isValid());
288 std::vector<PassiveHit> passiveHits;
289 passiveHits.insert(passiveHits.end(), hitsPassive->begin(), hitsPassive->end());
290 edm::LogVerbatim(
"HitStudy") <<
"CaloSimHitAnalysis: Passive Hit buffer " << passiveHits.size();
297 int nHit = hits.size();
298 int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nEB = 0, nEE = 0, nBad = 0, iHit = 0;
299 std::map<CaloHitID, std::pair<double, double> > hitMap;
301 for (
unsigned int k = 0;
k <
nCalo_; ++
k) {
302 etot[
k] = etotG[
k] = 0;
304 for (
const auto&
hit : hits) {
305 double edep =
hit.energy();
307 uint32_t
id =
hit.
id();
308 int itra =
hit.geantTrackId();
309 double edepEM =
hit.energyEM();
310 double edepHad =
hit.energyHad();
330 }
else if (subdet == static_cast<int>(
HcalEndcap)) {
333 }
else if (subdet == static_cast<int>(
HcalOuter)) {
336 }
else if (subdet == static_cast<int>(
HcalForward)) {
342 edm::LogVerbatim(
"HitStudy") <<
"Hit[" << iHit <<
":" << nHit <<
":" << idx <<
"] E " << edep <<
":" << edepEM
343 <<
":" << edepHad <<
" T " << time <<
" itra " << itra <<
" ID " << std::hex <<
id 349 auto itr = hitMap.find(hid);
350 if (itr == hitMap.end())
351 hitMap[hid] = std::make_pair(time, edep);
353 ((itr->second).second) += edep;
368 for (
auto itr = hitMap.begin(); itr != hitMap.end(); ++itr) {
371 DetId id((itr->first).unitID());
374 << std::hex <<
id.rawId() <<
std::dec;
380 int subdet =
id.subdetId();
383 }
else if (subdet == static_cast<int>(
HcalEndcap)) {
385 }
else if (subdet == static_cast<int>(
HcalOuter)) {
387 }
else if (subdet == static_cast<int>(
HcalForward)) {
392 double edep = (itr->second).
second;
395 edm::LogVerbatim(
"HitStudy") <<
"Index " << idx <<
":" << nCalo_ <<
" Point " << point <<
" E " << edep <<
" T " 418 h_etot_[indx]->Fill(etot[indx]);
421 h_hit_[indx]->Fill(
double(nEB));
423 h_hit_[indx]->Fill(
double(nEE));
425 h_hit_[2]->Fill(
double(nHB));
426 h_hit_[3]->Fill(
double(nHE));
427 h_hit_[4]->Fill(
double(nHO));
428 h_hit_[5]->Fill(
double(nHF));
435 edm::LogVerbatim(
"HitStudy") <<
"CaloSimHitAnalysis::analyzeHits: EB " << nEB <<
" EE " << nEE <<
" HB " << nHB
436 <<
" HE " << nHE <<
" HO " << nHO <<
" HF " << nHF <<
" Bad " << nBad <<
" All " << nHit
437 <<
" Reduced " << hitMap.size();
443 std::map<std::pair<std::string, uint32_t>,
int> hitx;
444 std::map<int, int>
tracks;
445 for (
auto&
hit : hits) {
447 std::pair<std::string, uint32_t> volume = std::make_pair(name, (
hit.
id() % 1000000));
448 auto itr = hitx.find(volume);
449 if (itr == hitx.end())
453 auto ktr = tracks.find(
hit.trackId());
454 if (ktr == tracks.end())
455 tracks[
hit.trackId()] = 1;
461 if ((name.find(active) != std::string::npos) || (name.find(sensor) != std::string::npos)) {
464 if (name.find(
detNames_[
k]) != std::string::npos) {
477 edm::LogVerbatim(
"HitStudy") <<
"CaloSimHitAnalysis::analyzPassiveHits: Total " << hits.size() <<
" Cells " 478 << hitx.size() <<
" Tracks " << tracks.size();
void analyzeHits(std::vector< PCaloHit > &, int)
static const std::string kSharedResource
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
EventNumber_t event() const
edm::EDGetTokenT< edm::PassiveHitContainer > tok_passive_
HcalSubdetector subdet() const
get the subdetector
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
void analyzePassiveHits(std::vector< PassiveHit > &hits)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
const std::vector< std::string > hitLab_
CaloSimHitAnalysis(const edm::ParameterSet &ps)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
std::vector< TH1F * > h_edepTk_
Geom::Phi< T > phi() const
constexpr uint32_t rawId() const
get the raw id
T * make(const Args &...args) const
make new ROOT object
std::vector< TH1F * > h_timeTk_
const std::string names[nVars_]
const std::string g4Label_
const std::vector< std::string > detNames_
U second(std::pair< T, U > const &p)
std::vector< edm::EDGetTokenT< edm::PCaloHitContainer > > toks_calo_
void endRun(edm::Run const &, edm::EventSetup const &) override
#define DEFINE_FWK_MODULE(type)
TH1F * h_edepHad_[nCalo_]
GlobalPoint getPosition(const DetId &id) const
Get the position of a given detector id.
Abs< T >::type abs(const T &t)
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
const HcalGeometry * hcalGeom_
ParameterDescriptionBase * add(U const &iLabel, T const &value)
GlobalPoint getPosition(const DetId &id) const
const CaloGeometry * caloGeometry_
void beginRun(edm::Run const &, edm::EventSetup const &) override
const std::vector< double > timeSliceUnit_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
void analyze(edm::Event const &, edm::EventSetup const &) override
T const * product() const
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
~CaloSimHitAnalysis() override