Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:16

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   // Book histograms
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   // total X0
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 }