19 binEta_ =
p.getUntrackedParameter<
int>(
"NBinEta", 260);
20 binPhi_ =
p.getUntrackedParameter<
int>(
"NBinPhi", 180);
21 maxEta_ =
p.getUntrackedParameter<
double>(
"MaxEta", 5.2);
22 etaLow_ =
p.getUntrackedParameter<
double>(
"EtaLow", -5.2);
23 etaHigh_ =
p.getUntrackedParameter<
double>(
"EtaHigh", 5.2);
24 fillHistos_ =
p.getUntrackedParameter<
bool>(
"FillHisto",
true);
25 printSum_ =
p.getUntrackedParameter<
bool>(
"PrintSummary",
false);
26 etaMinP_ =
p.getUntrackedParameter<
double>(
"EtaMinP", 5.2);
27 etaMaxP_ =
p.getUntrackedParameter<
double>(
"EtaMaxP", 0.0);
28 etaLowMin_ =
p.getUntrackedParameter<
double>(
"EtaLowMin", 0.783);
29 etaLowMax_ =
p.getUntrackedParameter<
double>(
"EtaLowMax", 0.870);
30 etaMidMin_ =
p.getUntrackedParameter<
double>(
"EtaMidMin", 2.650);
31 etaMidMax_ =
p.getUntrackedParameter<
double>(
"EtaMidMax", 2.868);
32 etaHighMin_ =
p.getUntrackedParameter<
double>(
"EtaHighMin", 2.868);
33 etaHighMax_ =
p.getUntrackedParameter<
double>(
"EtaHighMax", 3.000);
34 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos: FillHisto : " << fillHistos_ <<
" PrintSummary "
35 << printSum_ <<
" == Eta plot: NX " << binEta_ <<
" Range " << -maxEta_ <<
":"
36 << maxEta_ <<
" Phi plot: NX " << binPhi_ <<
" Range " << -1._pi <<
":" << 1._pi
37 <<
" (Eta limit " << etaLow_ <<
":" << etaHigh_ <<
")"
38 <<
" Eta range (" << etaLowMin_ <<
":" << etaLowMax_ <<
"), (" << etaMidMin_ <<
":"
39 << etaMidMax_ <<
"), (" << etaHighMin_ <<
":" << etaHighMax_
40 <<
") Debug for eta range " << etaMinP_ <<
":" << etaMaxP_;
52 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute <<
" = "
53 <<
value <<
" has " << sensitives_.size() <<
" elements";
54 for (
unsigned int i = 0;
i < sensitives_.size();
i++)
56 <<
"MaterialBudgetHcalHistos: sensitives[" <<
i <<
"] = " << sensitives_[
i];
64 std::vector<double>
temp = getDDDArray(
"Levels",
sv);
65 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute <<
" = "
66 <<
value <<
" has " << hfNames_.size() <<
" elements";
67 for (
unsigned int i = 0;
i < hfNames_.size();
i++) {
69 hfLevels_.push_back(
level);
71 <<
"MaterialBudgetHcalHistos: HF[" <<
i <<
"] = " << hfNames_[
i] <<
" at level " << hfLevels_[
i];
74 const std::string ecalRO[2] = {
"EcalHitsEB",
"EcalHitsEE"};
75 attribute =
"ReadOutName";
76 for (
int k = 0;
k < 2;
k++) {
80 std::vector<std::string> senstmp =
getNames(fv3);
81 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute
82 <<
" = " <<
value <<
" has " << senstmp.size() <<
" elements";
83 for (
unsigned int i = 0;
i < senstmp.size();
i++)
84 sensitiveEC_.push_back(senstmp[
i]);
86 for (
unsigned int i = 0;
i < sensitiveEC_.size();
i++)
88 <<
"MaterialBudgetHcalHistos:sensitiveEC[" <<
i <<
"] = " << sensitiveEC_[
i];
99 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute <<
" = "
100 <<
value <<
" has " << sensitives_.size() <<
" elements";
101 for (
unsigned int i = 0;
i < sensitives_.size();
i++)
103 <<
"MaterialBudgetHcalHistos: sensitives[" <<
i <<
"] = " << sensitives_[
i];
104 attribute =
"Volume";
108 std::vector<int>
temp = fv2.
get<std::vector<int> >(
"hf",
"Levels");
110 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute <<
" = "
111 <<
value <<
" has " << hfNames_.size() <<
" elements";
112 for (
unsigned int i = 0;
i < hfNames_.size();
i++) {
113 hfLevels_.push_back(
temp[
i]);
115 <<
"MaterialBudgetHcalHistos: HF[" <<
i <<
"] = " << hfNames_[
i] <<
" at level " << hfLevels_[
i];
118 const std::string ecalRO[2] = {
"EcalHitsEB",
"EcalHitsEE"};
119 attribute =
"ReadOutName";
120 for (
int k = 0;
k < 2;
k++) {
124 std::vector<std::string> senstmp =
getNames(fv);
125 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Names to be tested for " << attribute
126 <<
" = " <<
value <<
" has " << senstmp.size() <<
" elements";
127 for (
unsigned int i = 0;
i < senstmp.size();
i++)
128 sensitiveEC_.push_back(senstmp[
i]);
130 for (
unsigned int i = 0;
i < sensitiveEC_.size();
i++)
132 <<
"MaterialBudgetHcalHistos:sensitiveEC[" <<
i <<
"] = " << sensitiveEC_[
i];
137 id_ =
layer_ = steps_ = 0;
138 radLen_ = intLen_ = stepLen_ = 0;
139 nlayHB_ = nlayHE_ = nlayHF_ = nlayHO_ = 0;
141 const G4ThreeVector&
dir = aTrack->GetMomentum();
142 if (
dir.theta() != 0) {
148 double theEnergy = aTrack->GetTotalEnergy();
149 int theID = (
int)(aTrack->GetDefinition()->GetPDGEncoding());
159 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() <<
" Code "
160 << theID <<
" Energy " <<
convertUnitsTo(1._GeV, theEnergy) <<
" GeV; Eta "
166 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
167 double step = aStep->GetStepLength();
168 double radl = material->GetRadlen();
169 double intl = material->GetNuclearInterLength();
173 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
178 for (
unsigned int ii = 0;
ii < matList_.size();
ii++) {
179 if (matList_[
ii] == matName) {
181 radLength_[
ii] += (
step / radl);
182 intLength_[
ii] += (
step / intl);
188 matList_.push_back(matName);
189 stepLength_.push_back(
step);
190 radLength_.push_back(
step / radl);
191 intLength_.push_back(
step / intl);
195 <<
" Material " << matName <<
" Old Length " << stepLen_ <<
" X0 "
196 <<
step / radl <<
":" << radLen_ <<
" Lambda " <<
step / intl <<
":"
200 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos: Step at " <<
name <<
" id " << id_ <<
":"
201 << idOld <<
" Length " <<
step <<
" in " << matName <<
" of density "
202 <<
density <<
" g/cc; Radiation Length " << radl <<
" mm; Interaction Length "
203 << intl <<
" mm\n Position "
204 << aStep->GetPreStepPoint()->GetPosition() <<
" Cylindrical R "
205 << aStep->GetPreStepPoint()->GetPosition().perp() <<
" Length (so far) "
206 << stepLen_ <<
" L/X0 " <<
step / radl <<
"/" << radLen_ <<
" L/Lambda "
207 <<
step / intl <<
"/" << intLen_;
210 int det = 0, lay = 0;
217 if (isSensitive(
name)) {
224 det = (touch->GetReplicaNumber(1)) / 1000;
225 lay = (touch->GetReplicaNumber(0) / 10) % 100 + 3;
239 }
else if (lay !=
layer_) {
246 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Det " << det <<
" Layer " << lay <<
" Eta "
261 if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
263 <<
"MaterialBudgetHcalHistos: Step at " <<
name <<
" calls filHisto with " << (id_ - 1);
269 radLen_ += (
step / radl);
270 intLen_ += (
step / intl);
273 if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
274 if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
276 <<
"MaterialBudgetHcalHistos: After HF in " <<
name <<
":"
277 << aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName() <<
" calls fillHisto with " << id_;
288 edm::LogVerbatim(
"MaterialBudget") <<
"Number of layers hit in HB:" << nlayHB_ <<
" HE:" << nlayHE_
289 <<
" HO:" << nlayHO_ <<
" HF:" << nlayHF_;
291 fillHisto(maxSet_ - 1);
295 for (
unsigned int ii = 0;
ii < matList_.size();
ii++) {
296 edm::LogVerbatim(
"MaterialBudget") << matList_[
ii] <<
"\t" << stepLength_[
ii] <<
"\t" << radLength_[
ii] <<
"\t"
306 if (!
tfile.isAvailable())
308 <<
"please add it to config file";
311 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos: Booking user histos === with " << binEta_
312 <<
" bins in eta from " << -maxEta_ <<
" to " << maxEta_ <<
" and " << binPhi_
316 std::string range0 =
"(" + std::to_string(etaMidMin_) +
":" + std::to_string(etaMidMax_) +
") ";
317 std::string range1 =
"(" + std::to_string(etaHighMin_) +
":" + std::to_string(etaHighMax_) +
") ";
318 std::string range2 =
"(" + std::to_string(etaLowMin_) +
":" + std::to_string(etaLowMax_) +
") ";
320 for (
int i = 0;
i < maxSet_;
i++) {
321 iter = std::to_string(
i);
322 me100[
i] =
tfile->make<TProfile>(
323 std::to_string(
i + 100).c_str(), (
"MB(X0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
324 me200[
i] =
tfile->make<TProfile>(
325 std::to_string(
i + 200).c_str(), (
"MB(L0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
326 me300[
i] =
tfile->make<TProfile>(
327 std::to_string(
i + 300).c_str(), (
"MB(Step) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
328 me400[
i] =
tfile->make<TH1F>(
329 std::to_string(
i + 400).c_str(), (
"Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
330 me500[
i] =
tfile->make<TProfile>(
331 std::to_string(
i + 500).c_str(), (
"MB(X0) prof Ph in region " + iter).c_str(), binPhi_, -
maxPhi,
maxPhi);
332 me600[
i] =
tfile->make<TProfile>(
333 std::to_string(
i + 600).c_str(), (
"MB(L0) prof Ph in region " + iter).c_str(), binPhi_, -
maxPhi,
maxPhi);
334 me700[
i] =
tfile->make<TProfile>(
335 std::to_string(
i + 700).c_str(), (
"MB(Step) prof Ph in region " + iter).c_str(), binPhi_, -
maxPhi,
maxPhi);
337 tfile->make<TH1F>(std::to_string(
i + 800).c_str(), (
"Phi in region " + iter).c_str(), binPhi_, -
maxPhi,
maxPhi);
338 me900[
i] =
tfile->make<TProfile2D>(std::to_string(
i + 900).c_str(),
339 (
"MB(X0) prof Eta Phi in region " + iter).c_str(),
346 me1000[
i] =
tfile->make<TProfile2D>(std::to_string(
i + 1000).c_str(),
347 (
"MB(L0) prof Eta Phi in region " + iter).c_str(),
354 me1100[
i] =
tfile->make<TProfile2D>(std::to_string(
i + 1100).c_str(),
355 (
"MB(Step) prof Eta Phi in region " + iter).c_str(),
362 me1200[
i] =
tfile->make<TH2F>(std::to_string(
i + 1200).c_str(),
363 (
"Eta vs Phi in region " + iter).c_str(),
370 me1600[
i] =
tfile->make<TProfile>(std::to_string(
i + 1600).c_str(),
371 (
"MB(X0) prof Ph in region " + range0 + iter).c_str(),
375 me1700[
i] =
tfile->make<TProfile>(std::to_string(
i + 1700).c_str(),
376 (
"MB(L0) prof Ph in region " + range0 + iter).c_str(),
380 me1800[
i] =
tfile->make<TProfile>(std::to_string(
i + 1800).c_str(),
381 (
"MB(Step) prof Ph in region " + range0 + iter).c_str(),
385 me1900[
i] =
tfile->make<TProfile>(std::to_string(
i + 1900).c_str(),
386 (
"MB(X0) prof Ph in region " + range1 + iter).c_str(),
390 me2000[
i] =
tfile->make<TProfile>(std::to_string(
i + 2000).c_str(),
391 (
"MB(L0) prof Ph in region " + range1 + iter).c_str(),
395 me2100[
i] =
tfile->make<TProfile>(std::to_string(
i + 2100).c_str(),
396 (
"MB(Step) prof Ph in region " + range1 + iter).c_str(),
400 me2200[
i] =
tfile->make<TProfile>(std::to_string(
i + 2200).c_str(),
401 (
"MB(X0) prof Ph in region " + range2 + iter).c_str(),
405 me2300[
i] =
tfile->make<TProfile>(std::to_string(
i + 2300).c_str(),
406 (
"MB(L0) prof Ph in region " + range2 + iter).c_str(),
410 me2400[
i] =
tfile->make<TProfile>(std::to_string(
i + 2400).c_str(),
411 (
"MB(Step) prof Ph in region " + range2 + iter).c_str(),
416 for (
int i = 0;
i < maxSet2_;
i++) {
417 iter = std::to_string(
i);
418 me1300[
i] =
tfile->make<TH1F>(std::to_string(
i + 1300).c_str(),
419 (
"Events with layers Hit (0 all, 1 HB, ..) for " + iter).c_str(),
423 me1400[
i] =
tfile->make<TH2F>(std::to_string(
i + 1400).c_str(),
424 (
"Eta vs Phi for layers hit in " + iter).c_str(),
431 me1500[
i] =
tfile->make<TProfile>(std::to_string(
i + 1500).c_str(),
432 (
"Number of layers crossed (0 all, 1 HB, ..) for " + iter).c_str(),
438 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos: Booking user histos done ===";
443 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos:FillHisto called with index " <<
ii
444 <<
" integrated step " << stepLen_ <<
" X0 " << radLen_ <<
" Lamda " << intLen_;
446 if (
ii >= 0 &&
ii < maxSet_) {
447 me100[
ii]->Fill(eta_, radLen_);
448 me200[
ii]->Fill(eta_, intLen_);
449 me300[
ii]->Fill(eta_, stepLen_);
450 me400[
ii]->Fill(eta_);
452 if (eta_ >= etaLow_ && eta_ <= etaHigh_) {
453 me500[
ii]->Fill(phi_, radLen_);
454 me600[
ii]->Fill(phi_, intLen_);
455 me700[
ii]->Fill(phi_, stepLen_);
456 me800[
ii]->Fill(phi_);
459 me900[
ii]->Fill(eta_, phi_, radLen_);
460 me1000[
ii]->Fill(eta_, phi_, intLen_);
461 me1100[
ii]->Fill(eta_, phi_, stepLen_);
462 me1200[
ii]->Fill(eta_, phi_);
465 me1600[
ii]->Fill(phi_, radLen_);
466 me1700[
ii]->Fill(phi_, intLen_);
467 me1800[
ii]->Fill(phi_, stepLen_);
471 me1900[
ii]->Fill(phi_, radLen_);
472 me2000[
ii]->Fill(phi_, intLen_);
473 me2100[
ii]->Fill(phi_, stepLen_);
477 me2200[
ii]->Fill(phi_, radLen_);
478 me2300[
ii]->Fill(phi_, intLen_);
479 me2400[
ii]->Fill(phi_, stepLen_);
485 me1300[0]->Fill(eta_);
486 me1400[0]->Fill(eta_, phi_);
488 me1300[1]->Fill(eta_);
489 me1400[1]->Fill(eta_, phi_);
492 me1300[2]->Fill(eta_);
493 me1400[2]->Fill(eta_, phi_);
496 me1300[3]->Fill(eta_);
497 me1400[3]->Fill(eta_, phi_);
500 me1300[4]->Fill(eta_);
501 me1400[4]->Fill(eta_, phi_);
504 me1300[5]->Fill(eta_);
505 me1400[5]->Fill(eta_, phi_);
508 me1300[6]->Fill(eta_);
509 me1400[6]->Fill(eta_, phi_);
512 me1300[7]->Fill(eta_);
513 me1400[7]->Fill(eta_, phi_);
515 if (nlayHB_ > 0 || nlayHE_ > 0 || (nlayHF_ > 0 &&
std::abs(eta_) > 3.0)) {
516 me1300[8]->Fill(eta_);
517 me1400[8]->Fill(eta_, phi_);
519 me1500[0]->Fill(eta_, (
double)(nlayHB_ + nlayHO_ + nlayHE_ + nlayHF_));
520 me1500[1]->Fill(eta_, (
double)(nlayHB_));
521 me1500[2]->Fill(eta_, (
double)(nlayHE_));
522 me1500[4]->Fill(eta_, (
double)(nlayHF_));
526 edm::LogVerbatim(
"MaterialBudget") <<
"MaterialBudgetHcalHistos: Save user histos ===";
530 std::vector<std::string>
tmp;
543 std::vector<std::string>
tmp;
544 const std::vector<std::string> notIn = {
545 "CALO",
"HCal",
"MBBTL",
"MBBTR",
"MBBTC",
"MBAT",
"MBBT_R1M",
"MBBT_R1P",
"VCAL",
"HVQF"};
548 if (
std::find(notIn.begin(), notIn.end(),
n) == notIn.end()) {
559 edm::LogVerbatim(
"MaterialBudgetFull") <<
"MaterialBudgetHcalHistos:getDDDArray called for " <<
str;
563 const std::vector<double>& fvec =
value.doubles();
564 int nval = fvec.size();
566 throw cms::Exception(
"MaterialBudgetHcalHistos") <<
"nval = " << nval <<
" < 1 for array " <<
str <<
"\n";
571 throw cms::Exception(
"MaterialBudgetHcalHistos") <<
"cannot get array " <<
str <<
"\n";
576 std::vector<std::string>::const_iterator it = sensitives_.begin();
577 std::vector<std::string>::const_iterator itEnd = sensitives_.end();
578 for (; it != itEnd; ++it)
585 int levels = ((touch->GetHistoryDepth()) + 1);
586 for (
unsigned int it = 0; it < hfNames_.size(); it++) {
587 if (
levels >= hfLevels_[it]) {
589 if (
name == hfNames_[it]) {
598 std::vector<std::string>::const_iterator it = sensitiveEC_.begin();
599 std::vector<std::string>::const_iterator itEnd = sensitiveEC_.end();
600 for (; it != itEnd; ++it)