File indexing completed on 2023-03-17 11:27:26
0001 #include "Validation/Geometry/interface/MaterialBudgetHcalHistos.h"
0002
0003 #include "DataFormats/Math/interface/GeantUnits.h"
0004 #include "DetectorDescription/Core/interface/DDFilter.h"
0005 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0006 #include "DetectorDescription/Core/interface/DDSplit.h"
0007 #include "DetectorDescription/Core/interface/DDValue.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0012
0013 #include "DD4hep/Filter.h"
0014
0015 #include <string>
0016 #include <vector>
0017
0018 using namespace geant_units::operators;
0019
0020 MaterialBudgetHcalHistos::MaterialBudgetHcalHistos(const edm::ParameterSet& p) {
0021 binEta_ = p.getUntrackedParameter<int>("NBinEta", 260);
0022 binPhi_ = p.getUntrackedParameter<int>("NBinPhi", 180);
0023 maxEta_ = p.getUntrackedParameter<double>("MaxEta", 5.2);
0024 etaLow_ = p.getUntrackedParameter<double>("EtaLow", -5.2);
0025 etaHigh_ = p.getUntrackedParameter<double>("EtaHigh", 5.2);
0026 fillHistos_ = p.getUntrackedParameter<bool>("FillHisto", true);
0027 printSum_ = p.getUntrackedParameter<bool>("PrintSummary", false);
0028 fromdd4hep_ = p.getUntrackedParameter<bool>("Fromdd4hep", false);
0029 etaMinP_ = p.getUntrackedParameter<double>("EtaMinP", 5.2);
0030 etaMaxP_ = p.getUntrackedParameter<double>("EtaMaxP", 0.0);
0031 etaLowMin_ = p.getUntrackedParameter<double>("EtaLowMin", 0.783);
0032 etaLowMax_ = p.getUntrackedParameter<double>("EtaLowMax", 0.870);
0033 etaMidMin_ = p.getUntrackedParameter<double>("EtaMidMin", 2.650);
0034 etaMidMax_ = p.getUntrackedParameter<double>("EtaMidMax", 2.868);
0035 etaHighMin_ = p.getUntrackedParameter<double>("EtaHighMin", 2.868);
0036 etaHighMax_ = p.getUntrackedParameter<double>("EtaHighMax", 3.000);
0037 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: FillHisto : " << fillHistos_ << " PrintSummary "
0038 << printSum_ << " == Eta plot: NX " << binEta_ << " Range " << -maxEta_ << ":"
0039 << maxEta_ << " Phi plot: NX " << binPhi_ << " Range " << -1._pi << ":" << 1._pi
0040 << " (Eta limit " << etaLow_ << ":" << etaHigh_ << ")"
0041 << " Eta range (" << etaLowMin_ << ":" << etaLowMax_ << "), (" << etaMidMin_ << ":"
0042 << etaMidMax_ << "), (" << etaHighMin_ << ":" << etaHighMax_
0043 << ") Debug for eta range " << etaMinP_ << ":" << etaMaxP_ << " FromDD4hep "
0044 << fromdd4hep_;
0045 if (fillHistos_)
0046 book();
0047 }
0048
0049 void MaterialBudgetHcalHistos::fillBeginJob(const DDCompactView& cpv) {
0050 constexpr int32_t addLevel = 0;
0051 if (fillHistos_) {
0052 std::string attribute = "ReadOutName";
0053 std::string value = "HcalHits";
0054 DDSpecificsMatchesValueFilter filter1{DDValue(attribute, value, 0)};
0055 DDFilteredView fv1(cpv, filter1);
0056 std::vector<std::string> names = getNames(fv1);
0057 for (auto& name : names) {
0058 std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0059 if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
0060 sensitives_.emplace_back(namx);
0061 }
0062 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
0063 << value << " has " << sensitives_.size() << " elements";
0064 for (unsigned int i = 0; i < sensitives_.size(); i++)
0065 edm::LogVerbatim("MaterialBudgetFull")
0066 << "MaterialBudgetHcalHistos: sensitives[" << i << "] = " << sensitives_[i];
0067 attribute = "Volume";
0068 value = "HF";
0069 DDSpecificsMatchesValueFilter filter2{DDValue(attribute, value, 0)};
0070 DDFilteredView fv2(cpv, filter2);
0071 hfNames_ = getNames(fv2);
0072 fv2.firstChild();
0073 DDsvalues_type sv(fv2.mergedSpecifics());
0074 std::vector<double> temp = getDDDArray("Levels", sv);
0075 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
0076 << value << " has " << hfNames_.size() << " elements";
0077 for (unsigned int i = 0; i < hfNames_.size(); i++) {
0078 int level = static_cast<int>(temp[i]);
0079 hfLevels_.push_back(level + addLevel);
0080 edm::LogVerbatim("MaterialBudgetFull")
0081 << "MaterialBudgetHcalHistos: HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
0082 }
0083
0084 const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
0085 attribute = "ReadOutName";
0086 for (int k = 0; k < 2; k++) {
0087 value = ecalRO[k];
0088 DDSpecificsMatchesValueFilter filter3{DDValue(attribute, value, 0)};
0089 DDFilteredView fv3(cpv, filter3);
0090 std::vector<std::string> senstmp = getNames(fv3);
0091 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute
0092 << " = " << value << " has " << senstmp.size() << " elements";
0093 for (unsigned int i = 0; i < senstmp.size(); i++) {
0094 std::string name = senstmp[i].substr(0, 4);
0095 if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
0096 sensitiveEC_.push_back(name);
0097 }
0098 }
0099 for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
0100 edm::LogVerbatim("MaterialBudgetFull")
0101 << "MaterialBudgetHcalHistos:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
0102 }
0103 }
0104
0105 void MaterialBudgetHcalHistos::fillBeginJob(const cms::DDCompactView& cpv) {
0106 constexpr int32_t addLevel = 1;
0107 if (fillHistos_) {
0108 std::string attribute = "ReadOutName";
0109 std::string value = "HcalHits";
0110 const cms::DDFilter filter1(attribute, value);
0111 cms::DDFilteredView fv1(cpv, filter1);
0112 std::vector<std::string> names = getNames(fv1);
0113 for (auto& name : names) {
0114 std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0115 if (std::find(sensitives_.begin(), sensitives_.end(), namx) == sensitives_.end())
0116 sensitives_.emplace_back(namx);
0117 }
0118 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
0119 << value << " has " << sensitives_.size() << " elements";
0120 for (unsigned int i = 0; i < sensitives_.size(); i++)
0121 edm::LogVerbatim("MaterialBudgetFull")
0122 << "MaterialBudgetHcalHistos: sensitives[" << i << "] = " << sensitives_[i];
0123 attribute = "Volume";
0124 value = "HF";
0125 const cms::DDFilter filter2(attribute, value);
0126 cms::DDFilteredView fv2(cpv, filter2);
0127 std::vector<int> temp = fv2.get<std::vector<int> >("hf", "Levels");
0128 hfNames_ = getNames(fv2);
0129 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute << " = "
0130 << value << " has " << hfNames_.size() << " elements";
0131 for (unsigned int i = 0; i < hfNames_.size(); i++) {
0132 hfLevels_.push_back(temp[i] + addLevel);
0133 edm::LogVerbatim("MaterialBudgetFull")
0134 << "MaterialBudgetHcalHistos: HF[" << i << "] = " << hfNames_[i] << " at level " << hfLevels_[i];
0135 }
0136
0137 const std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
0138 attribute = "ReadOutName";
0139 for (int k = 0; k < 2; k++) {
0140 value = ecalRO[k];
0141 const cms::DDFilter filter(attribute, value);
0142 cms::DDFilteredView fv(cpv, filter);
0143 std::vector<std::string> senstmp = getNames(fv);
0144 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Names to be tested for " << attribute
0145 << " = " << value << " has " << senstmp.size() << " elements";
0146 for (unsigned int i = 0; i < senstmp.size(); i++) {
0147 std::string name = senstmp[i].substr(0, 4);
0148 if (std::find(sensitiveEC_.begin(), sensitiveEC_.end(), name) == sensitiveEC_.end())
0149 sensitiveEC_.push_back(name);
0150 }
0151 }
0152 for (unsigned int i = 0; i < sensitiveEC_.size(); i++)
0153 edm::LogVerbatim("MaterialBudgetFull")
0154 << "MaterialBudgetHcalHistos:sensitiveEC[" << i << "] = " << sensitiveEC_[i];
0155 }
0156 }
0157
0158 void MaterialBudgetHcalHistos::fillStartTrack(const G4Track* aTrack) {
0159 id_ = layer_ = steps_ = 0;
0160 radLen_ = intLen_ = stepLen_ = 0;
0161 nlayHB_ = nlayHE_ = nlayHF_ = nlayHO_ = 0;
0162
0163 const G4ThreeVector& dir = aTrack->GetMomentum();
0164 if (dir.theta() != 0) {
0165 eta_ = dir.eta();
0166 } else {
0167 eta_ = -99;
0168 }
0169 phi_ = dir.phi();
0170 double theEnergy = aTrack->GetTotalEnergy();
0171 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0172
0173 if (printSum_) {
0174 matList_.clear();
0175 stepLength_.clear();
0176 radLength_.clear();
0177 intLength_.clear();
0178 }
0179
0180 if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0181 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code "
0182 << theID << " Energy " << convertUnitsTo(1._GeV, theEnergy) << " GeV; Eta "
0183 << eta_ << " Phi " << convertRadToDeg(phi_) << " PT "
0184 << convertUnitsTo(1._GeV, dir.perp()) << " GeV *****";
0185 }
0186
0187 void MaterialBudgetHcalHistos::fillPerStep(const G4Step* aStep) {
0188 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0189 double step = aStep->GetStepLength();
0190 double radl = material->GetRadlen();
0191 double intl = material->GetNuclearInterLength();
0192 double density = convertUnitsTo(1._g_per_cm3, material->GetDensity());
0193
0194 int idOld = id_;
0195 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0196 std::string name = (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(0)->GetName())));
0197 std::string matName = (static_cast<std::string>(dd4hep::dd::noNamespace(material->GetName())));
0198 if (printSum_) {
0199 bool found = false;
0200 for (unsigned int ii = 0; ii < matList_.size(); ii++) {
0201 if (matList_[ii] == matName) {
0202 stepLength_[ii] += step;
0203 radLength_[ii] += (step / radl);
0204 intLength_[ii] += (step / intl);
0205 found = true;
0206 break;
0207 }
0208 }
0209 if (!found) {
0210 matList_.push_back(matName);
0211 stepLength_.push_back(step);
0212 radLength_.push_back(step / radl);
0213 intLength_.push_back(step / intl);
0214 }
0215 if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0216 edm::LogVerbatim("MaterialBudget") << "Volume " << name << " id " << id_ << ":" << idOld << " Step " << step
0217 << " Material " << matName << " Old Length " << stepLen_ << " X0 "
0218 << step / radl << ":" << radLen_ << " Lambda " << step / intl << ":"
0219 << intLen_;
0220 } else {
0221 if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0222 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Step at " << name << " id " << id_ << ":"
0223 << idOld << " Length " << step << " in " << matName << " of density "
0224 << density << " g/cc; Radiation Length " << radl << " mm; Interaction Length "
0225 << intl << " mm\n Position "
0226 << aStep->GetPreStepPoint()->GetPosition() << " Cylindrical R "
0227 << aStep->GetPreStepPoint()->GetPosition().perp() << " Length (so far) "
0228 << stepLen_ << " L/X0 " << step / radl << "/" << radLen_ << " L/Lambda "
0229 << step / intl << "/" << intLen_;
0230 }
0231
0232 int det = 0, lay = 0;
0233 double abseta = std::abs(eta_);
0234 if (fillHistos_) {
0235 edm::LogVerbatim("MaterialBudgetFull")
0236 << "Volume " << name << ":" << matName << " EC:Sensitive:HF " << isItEC(name) << ":" << isSensitive(name) << ":"
0237 << isItHF(touch) << " Eta " << abseta << " HC " << ((touch->GetReplicaNumber(1)) / 1000) << ":"
0238 << ((touch->GetReplicaNumber(0) / 10) % 100 + 3) << " X0 " << (radLen_ + (step / radl)) << " Lambda "
0239 << (intLen_ + (step / intl));
0240 ;
0241 if (isItEC(name)) {
0242 det = 1;
0243 lay = 1;
0244 } else {
0245 if (isSensitive(name)) {
0246 if (isItHF(touch)) {
0247 det = 5;
0248 lay = 21;
0249 if (lay != layer_)
0250 ++nlayHF_;
0251 } else {
0252 det = (touch->GetReplicaNumber(1)) / 1000;
0253 lay = (touch->GetReplicaNumber(0) / 10) % 100 + 3;
0254 if (det == 4) {
0255 if (abseta < 1.479)
0256 lay = layer_ + 1;
0257 else
0258 lay--;
0259 if (lay < 3)
0260 lay = 3;
0261 if (lay == layer_)
0262 lay++;
0263 if (lay > 20)
0264 lay = 20;
0265 if (lay != layer_)
0266 ++nlayHE_;
0267 } else if (lay != layer_) {
0268 if (lay >= 20)
0269 ++nlayHO_;
0270 else
0271 ++nlayHB_;
0272 }
0273 }
0274 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Det " << det << " Layer " << lay << " Eta "
0275 << eta_ << " Phi " << convertRadToDeg(phi_);
0276 } else if (layer_ == 1) {
0277 det = -1;
0278 lay = 2;
0279 }
0280 }
0281 if (det != 0) {
0282 if (lay != layer_) {
0283 id_ = lay;
0284 layer_ = lay;
0285 }
0286 }
0287
0288 if (id_ > idOld) {
0289 if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
0290 edm::LogVerbatim("MaterialBudget")
0291 << "MaterialBudgetHcalHistos: Step at " << name << " calls filHisto with " << (id_ - 1);
0292 fillHisto(id_ - 1);
0293 }
0294 }
0295
0296 stepLen_ += step;
0297 radLen_ += (step / radl);
0298 intLen_ += (step / intl);
0299 if (fillHistos_) {
0300 if (id_ == 21) {
0301 if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
0302 if ((abseta >= etaMinP_) && (abseta <= etaMaxP_))
0303 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: After HF in " << name << ":"
0304 << static_cast<std::string>(dd4hep::dd::noNamespace(
0305 aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName()))
0306 << " calls fillHisto with " << id_;
0307 fillHisto(idOld);
0308 ++id_;
0309 layer_ = 0;
0310 }
0311 }
0312 }
0313 }
0314
0315 void MaterialBudgetHcalHistos::fillEndTrack() {
0316 if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0317 edm::LogVerbatim("MaterialBudget") << "Number of layers hit in HB:" << nlayHB_ << " HE:" << nlayHE_
0318 << " HO:" << nlayHO_ << " HF:" << nlayHF_;
0319 if (fillHistos_) {
0320 fillHisto(maxSet_ - 1);
0321 fillLayer();
0322 }
0323 if (printSum_) {
0324 for (unsigned int ii = 0; ii < matList_.size(); ii++) {
0325 edm::LogVerbatim("MaterialBudget") << matList_[ii] << "\t" << stepLength_[ii] << "\t" << radLength_[ii] << "\t"
0326 << intLength_[ii];
0327 }
0328 }
0329 }
0330
0331 void MaterialBudgetHcalHistos::book() {
0332
0333 edm::Service<TFileService> tfile;
0334
0335 if (!tfile.isAvailable())
0336 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0337 << "please add it to config file";
0338
0339 double maxPhi = 1._pi;
0340 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos: Booking user histos === with " << binEta_
0341 << " bins in eta from " << -maxEta_ << " to " << maxEta_ << " and " << binPhi_
0342 << " bins in phi from " << -maxPhi << " to " << maxPhi;
0343
0344 std::string iter;
0345 std::string range0 = "(" + std::to_string(etaMidMin_) + ":" + std::to_string(etaMidMax_) + ") ";
0346 std::string range1 = "(" + std::to_string(etaHighMin_) + ":" + std::to_string(etaHighMax_) + ") ";
0347 std::string range2 = "(" + std::to_string(etaLowMin_) + ":" + std::to_string(etaLowMax_) + ") ";
0348
0349 for (int i = 0; i < maxSet_; i++) {
0350 iter = std::to_string(i);
0351 me100[i] = tfile->make<TProfile>(
0352 std::to_string(i + 100).c_str(), ("MB(X0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0353 me200[i] = tfile->make<TProfile>(
0354 std::to_string(i + 200).c_str(), ("MB(L0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0355 me300[i] = tfile->make<TProfile>(
0356 std::to_string(i + 300).c_str(), ("MB(Step) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0357 me400[i] = tfile->make<TH1F>(
0358 std::to_string(i + 400).c_str(), ("Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0359 me500[i] = tfile->make<TProfile>(
0360 std::to_string(i + 500).c_str(), ("MB(X0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0361 me600[i] = tfile->make<TProfile>(
0362 std::to_string(i + 600).c_str(), ("MB(L0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0363 me700[i] = tfile->make<TProfile>(
0364 std::to_string(i + 700).c_str(), ("MB(Step) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0365 me800[i] =
0366 tfile->make<TH1F>(std::to_string(i + 800).c_str(), ("Phi in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0367 me900[i] = tfile->make<TProfile2D>(std::to_string(i + 900).c_str(),
0368 ("MB(X0) prof Eta Phi in region " + iter).c_str(),
0369 binEta_ / 2,
0370 -maxEta_,
0371 maxEta_,
0372 binPhi_ / 2,
0373 -maxPhi,
0374 maxPhi);
0375 me1000[i] = tfile->make<TProfile2D>(std::to_string(i + 1000).c_str(),
0376 ("MB(L0) prof Eta Phi in region " + iter).c_str(),
0377 binEta_ / 2,
0378 -maxEta_,
0379 maxEta_,
0380 binPhi_ / 2,
0381 -maxPhi,
0382 maxPhi);
0383 me1100[i] = tfile->make<TProfile2D>(std::to_string(i + 1100).c_str(),
0384 ("MB(Step) prof Eta Phi in region " + iter).c_str(),
0385 binEta_ / 2,
0386 -maxEta_,
0387 maxEta_,
0388 binPhi_ / 2,
0389 -maxPhi,
0390 maxPhi);
0391 me1200[i] = tfile->make<TH2F>(std::to_string(i + 1200).c_str(),
0392 ("Eta vs Phi in region " + iter).c_str(),
0393 binEta_ / 2,
0394 -maxEta_,
0395 maxEta_,
0396 binPhi_ / 2,
0397 -maxPhi,
0398 maxPhi);
0399 me1600[i] = tfile->make<TProfile>(std::to_string(i + 1600).c_str(),
0400 ("MB(X0) prof Ph in region " + range0 + iter).c_str(),
0401 binPhi_,
0402 -maxPhi,
0403 maxPhi);
0404 me1700[i] = tfile->make<TProfile>(std::to_string(i + 1700).c_str(),
0405 ("MB(L0) prof Ph in region " + range0 + iter).c_str(),
0406 binPhi_,
0407 -maxPhi,
0408 maxPhi);
0409 me1800[i] = tfile->make<TProfile>(std::to_string(i + 1800).c_str(),
0410 ("MB(Step) prof Ph in region " + range0 + iter).c_str(),
0411 binPhi_,
0412 -maxPhi,
0413 maxPhi);
0414 me1900[i] = tfile->make<TProfile>(std::to_string(i + 1900).c_str(),
0415 ("MB(X0) prof Ph in region " + range1 + iter).c_str(),
0416 binPhi_,
0417 -maxPhi,
0418 maxPhi);
0419 me2000[i] = tfile->make<TProfile>(std::to_string(i + 2000).c_str(),
0420 ("MB(L0) prof Ph in region " + range1 + iter).c_str(),
0421 binPhi_,
0422 -maxPhi,
0423 maxPhi);
0424 me2100[i] = tfile->make<TProfile>(std::to_string(i + 2100).c_str(),
0425 ("MB(Step) prof Ph in region " + range1 + iter).c_str(),
0426 binPhi_,
0427 -maxPhi,
0428 maxPhi);
0429 me2200[i] = tfile->make<TProfile>(std::to_string(i + 2200).c_str(),
0430 ("MB(X0) prof Ph in region " + range2 + iter).c_str(),
0431 binPhi_,
0432 -maxPhi,
0433 maxPhi);
0434 me2300[i] = tfile->make<TProfile>(std::to_string(i + 2300).c_str(),
0435 ("MB(L0) prof Ph in region " + range2 + iter).c_str(),
0436 binPhi_,
0437 -maxPhi,
0438 maxPhi);
0439 me2400[i] = tfile->make<TProfile>(std::to_string(i + 2400).c_str(),
0440 ("MB(Step) prof Ph in region " + range2 + iter).c_str(),
0441 binPhi_,
0442 -maxPhi,
0443 maxPhi);
0444 }
0445 for (int i = 0; i < maxSet2_; i++) {
0446 iter = std::to_string(i);
0447 me1300[i] = tfile->make<TH1F>(std::to_string(i + 1300).c_str(),
0448 ("Events with layers Hit (0 all, 1 HB, ..) for " + iter).c_str(),
0449 binEta_,
0450 -maxEta_,
0451 maxEta_);
0452 me1400[i] = tfile->make<TH2F>(std::to_string(i + 1400).c_str(),
0453 ("Eta vs Phi for layers hit in " + iter).c_str(),
0454 binEta_ / 2,
0455 -maxEta_,
0456 maxEta_,
0457 binPhi_ / 2,
0458 -maxPhi,
0459 maxPhi);
0460 me1500[i] = tfile->make<TProfile>(std::to_string(i + 1500).c_str(),
0461 ("Number of layers crossed (0 all, 1 HB, ..) for " + iter).c_str(),
0462 binEta_,
0463 -maxEta_,
0464 maxEta_);
0465 }
0466
0467 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user histos done ===";
0468 }
0469
0470 void MaterialBudgetHcalHistos::fillHisto(int ii) {
0471 if ((std::abs(eta_) >= etaMinP_) && (std::abs(eta_) <= etaMaxP_))
0472 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos:FillHisto called with index " << ii
0473 << " integrated step " << stepLen_ << " X0 " << radLen_ << " Lamda " << intLen_;
0474
0475 if (ii >= 0 && ii < maxSet_) {
0476 me100[ii]->Fill(eta_, radLen_);
0477 me200[ii]->Fill(eta_, intLen_);
0478 me300[ii]->Fill(eta_, stepLen_);
0479 me400[ii]->Fill(eta_);
0480
0481 if (eta_ >= etaLow_ && eta_ <= etaHigh_) {
0482 me500[ii]->Fill(phi_, radLen_);
0483 me600[ii]->Fill(phi_, intLen_);
0484 me700[ii]->Fill(phi_, stepLen_);
0485 me800[ii]->Fill(phi_);
0486 }
0487
0488 me900[ii]->Fill(eta_, phi_, radLen_);
0489 me1000[ii]->Fill(eta_, phi_, intLen_);
0490 me1100[ii]->Fill(eta_, phi_, stepLen_);
0491 me1200[ii]->Fill(eta_, phi_);
0492
0493 if ((std::abs(eta_) >= etaMidMin_) && (std::abs(eta_) <= etaMidMax_)) {
0494 me1600[ii]->Fill(phi_, radLen_);
0495 me1700[ii]->Fill(phi_, intLen_);
0496 me1800[ii]->Fill(phi_, stepLen_);
0497 }
0498
0499 if ((std::abs(eta_) >= etaHighMin_) && (std::abs(eta_) <= etaHighMax_)) {
0500 me1900[ii]->Fill(phi_, radLen_);
0501 me2000[ii]->Fill(phi_, intLen_);
0502 me2100[ii]->Fill(phi_, stepLen_);
0503 }
0504
0505 if ((std::abs(eta_) >= etaLowMin_) && (std::abs(eta_) <= etaLowMax_)) {
0506 me2200[ii]->Fill(phi_, radLen_);
0507 me2300[ii]->Fill(phi_, intLen_);
0508 me2400[ii]->Fill(phi_, stepLen_);
0509 }
0510 }
0511 }
0512
0513 void MaterialBudgetHcalHistos::fillLayer() {
0514 me1300[0]->Fill(eta_);
0515 me1400[0]->Fill(eta_, phi_);
0516 if (nlayHB_ > 0) {
0517 me1300[1]->Fill(eta_);
0518 me1400[1]->Fill(eta_, phi_);
0519 }
0520 if (nlayHB_ >= 16) {
0521 me1300[2]->Fill(eta_);
0522 me1400[2]->Fill(eta_, phi_);
0523 }
0524 if (nlayHE_ > 0) {
0525 me1300[3]->Fill(eta_);
0526 me1400[3]->Fill(eta_, phi_);
0527 }
0528 if (nlayHE_ >= 16) {
0529 me1300[4]->Fill(eta_);
0530 me1400[4]->Fill(eta_, phi_);
0531 }
0532 if (nlayHO_ > 0) {
0533 me1300[5]->Fill(eta_);
0534 me1400[5]->Fill(eta_, phi_);
0535 }
0536 if (nlayHO_ >= 2) {
0537 me1300[6]->Fill(eta_);
0538 me1400[6]->Fill(eta_, phi_);
0539 }
0540 if (nlayHF_ > 0) {
0541 me1300[7]->Fill(eta_);
0542 me1400[7]->Fill(eta_, phi_);
0543 }
0544 if (nlayHB_ > 0 || nlayHE_ > 0 || (nlayHF_ > 0 && std::abs(eta_) > 3.0)) {
0545 me1300[8]->Fill(eta_);
0546 me1400[8]->Fill(eta_, phi_);
0547 }
0548 me1500[0]->Fill(eta_, (double)(nlayHB_ + nlayHO_ + nlayHE_ + nlayHF_));
0549 me1500[1]->Fill(eta_, (double)(nlayHB_));
0550 me1500[2]->Fill(eta_, (double)(nlayHE_));
0551 me1500[4]->Fill(eta_, (double)(nlayHF_));
0552 }
0553
0554 void MaterialBudgetHcalHistos::hend() {
0555 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalHistos: Save user histos ===";
0556 }
0557
0558 std::vector<std::string> MaterialBudgetHcalHistos::getNames(DDFilteredView& fv) {
0559 std::vector<std::string> tmp;
0560 bool dodet = fv.firstChild();
0561 while (dodet) {
0562 const DDLogicalPart& log = fv.logicalPart();
0563 std::string namx = log.name().name();
0564 if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
0565 tmp.push_back(namx);
0566 dodet = fv.next();
0567 }
0568 return tmp;
0569 }
0570
0571 std::vector<std::string> MaterialBudgetHcalHistos::getNames(cms::DDFilteredView& fv) {
0572 std::vector<std::string> tmp;
0573 const std::vector<std::string> notIn = {
0574 "CALO", "HCal", "MBBTL", "MBBTR", "MBBTC", "MBAT", "MBBT_R1M", "MBBT_R1P", "MBBT_R1MX", "MBBT_R1PX", "VCAL"};
0575 while (fv.firstChild()) {
0576 const std::string n{fv.name().data(), fv.name().size()};
0577 if (std::find(notIn.begin(), notIn.end(), n) == notIn.end()) {
0578 std::string::size_type pos = n.find(':');
0579 const std::string namx = (pos == std::string::npos) ? n : std::string(n, pos + 1, n.size() - 1);
0580 if (std::find(tmp.begin(), tmp.end(), namx) == tmp.end())
0581 tmp.push_back(namx);
0582 }
0583 }
0584 return tmp;
0585 }
0586
0587 std::vector<double> MaterialBudgetHcalHistos::getDDDArray(const std::string& str, const DDsvalues_type& sv) {
0588 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalHistos:getDDDArray called for " << str;
0589 DDValue value(str);
0590 if (DDfetch(&sv, value)) {
0591 edm::LogVerbatim("MaterialBudgetFull") << value;
0592 const std::vector<double>& fvec = value.doubles();
0593 int nval = fvec.size();
0594 if (nval < 1) {
0595 throw cms::Exception("MaterialBudgetHcalHistos") << "nval = " << nval << " < 1 for array " << str << "\n";
0596 }
0597
0598 return fvec;
0599 } else {
0600 throw cms::Exception("MaterialBudgetHcalHistos") << "cannot get array " << str << "\n";
0601 }
0602 }
0603
0604 bool MaterialBudgetHcalHistos::isSensitive(const std::string& name) {
0605 std::vector<std::string>::const_iterator it = sensitives_.begin();
0606 std::vector<std::string>::const_iterator itEnd = sensitives_.end();
0607 std::string namx = (name.find('_') == std::string::npos) ? name : name.substr(0, name.find('_'));
0608 for (; it != itEnd; ++it)
0609 if (namx == *it)
0610 return true;
0611 return false;
0612 }
0613
0614 bool MaterialBudgetHcalHistos::isItHF(const G4VTouchable* touch) {
0615 int levels = ((touch->GetHistoryDepth()) + 1);
0616 for (unsigned int it = 0; it < hfNames_.size(); it++) {
0617 if (levels >= hfLevels_[it]) {
0618 std::string name =
0619 (static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(levels - hfLevels_[it])->GetName())))
0620 .substr(0, 4);
0621 if (name == hfNames_[it]) {
0622 return true;
0623 }
0624 }
0625 }
0626 return false;
0627 }
0628
0629 bool MaterialBudgetHcalHistos::isItEC(const std::string& name) {
0630 std::vector<std::string>::const_iterator it = sensitiveEC_.begin();
0631 std::vector<std::string>::const_iterator itEnd = sensitiveEC_.end();
0632 for (; it != itEnd; ++it)
0633 if (name.substr(0, 4) == *it)
0634 return true;
0635 return false;
0636 }
