File indexing completed on 2023-03-17 11:27:25
0001 #include "Validation/Geometry/interface/MaterialBudgetCastorHistos.h"
0002
0003 #include "DetectorDescription/Core/interface/DDFilter.h"
0004 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0005 #include "DetectorDescription/Core/interface/DDSplit.h"
0006 #include "DetectorDescription/Core/interface/DDValue.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ServiceRegistry/interface/Service.h"
0010 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0011
0012 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0013 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0014
0015 #include <string>
0016
0017 MaterialBudgetCastorHistos::MaterialBudgetCastorHistos(const edm::ParameterSet& p) {
0018 binEta = p.getUntrackedParameter<int>("NBinEta", 100);
0019 binPhi = p.getUntrackedParameter<int>("NBinPhi", 180);
0020 etaLow = p.getUntrackedParameter<double>("EtaLow", 5.0);
0021 etaHigh = p.getUntrackedParameter<double>("EtaHigh", 7.0);
0022 fillHistos = p.getUntrackedParameter<bool>("FillHisto", true);
0023 printSum = p.getUntrackedParameter<bool>("PrintSummary", false);
0024 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: FillHisto : " << fillHistos << " PrintSummary "
0025 << printSum << " == Eta plot: NX " << binEta << " Range " << etaLow << " (" << -etaHigh
0026 << ") : " << etaHigh << " (" << -etaLow << ") Phi plot: "
0027 << "NX " << binPhi << " Range " << -pi << ":" << pi << " (Eta limit " << etaLow << ":"
0028 << etaHigh << ")";
0029 if (fillHistos)
0030 book();
0031 }
0032
0033 MaterialBudgetCastorHistos::~MaterialBudgetCastorHistos() {
0034 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Save user "
0035 << "histos ===";
0036 }
0037
0038 void MaterialBudgetCastorHistos::fillStartTrack(const G4Track* aTrack) {
0039 id1 = id2 = steps = 0;
0040 radLen = intLen = stepLen = 0;
0041
0042 const G4ThreeVector& dir = aTrack->GetMomentum();
0043 if (dir.theta() != 0) {
0044 eta_ = dir.eta();
0045 } else {
0046 eta_ = -99;
0047 }
0048 phi_ = dir.phi();
0049 double theEnergy = aTrack->GetTotalEnergy();
0050 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0051
0052 if (printSum) {
0053 matList.clear();
0054 stepLength.clear();
0055 radLength.clear();
0056 intLength.clear();
0057 }
0058
0059 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Track " << aTrack->GetTrackID() << " Code " << theID
0060 << " Energy " << theEnergy / GeV << " GeV; Eta " << eta_ << " Phi " << phi_ / deg
0061 << " PT " << dir.perp() / GeV << " GeV *****";
0062 }
0063
0064 void MaterialBudgetCastorHistos::fillPerStep(const G4Step* aStep) {
0065 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0066 double step = aStep->GetStepLength();
0067 double radl = material->GetRadlen();
0068 double intl = material->GetNuclearInterLength();
0069 double density = material->GetDensity() / (g / cm3);
0070
0071 int id1Old = id1;
0072 int id2Old = id2;
0073 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0074 std::string name = touch->GetVolume(0)->GetName();
0075 const std::string& matName = material->GetName();
0076 if (printSum) {
0077 bool found = false;
0078 for (unsigned int ii = 0; ii < matList.size(); ii++) {
0079 if (matList[ii] == matName) {
0080 stepLength[ii] += step;
0081 radLength[ii] += (step / radl);
0082 intLength[ii] += (step / intl);
0083 found = true;
0084 break;
0085 }
0086 }
0087 if (!found) {
0088 matList.push_back(matName);
0089 stepLength.push_back(step);
0090 radLength.push_back(step / radl);
0091 intLength.push_back(step / intl);
0092 }
0093 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: " << name << " " << step << " " << matName << " "
0094 << stepLen << " " << step / radl << " " << radLen << " " << step / intl << " "
0095 << intLen;
0096 } else {
0097 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Step at " << name << " Length " << step << " in "
0098 << matName << " of density " << density << " g/cc; Radiation Length " << radl
0099 << " mm;"
0100 << " Interaction Length " << intl << " mm\n"
0101 << " Position " << aStep->GetPreStepPoint()->GetPosition()
0102 << " Cylindrical R " << aStep->GetPreStepPoint()->GetPosition().perp()
0103 << " Length (so far) " << stepLen << " L/X0 " << step / radl << "/" << radLen
0104 << " L/Lambda " << step / intl << "/" << intLen;
0105 }
0106
0107 int level = ((touch->GetHistoryDepth()) + 1);
0108 std::string name1 = "XXXX", name2 = "XXXX";
0109 if (level > 3)
0110 name1 = touch->GetVolume(level - 4)->GetName();
0111 if (level > 4)
0112 name2 = touch->GetVolume(level - 5)->GetName();
0113 if (name1 == "CAST") {
0114 id1 = 1;
0115 if (name2 == "CAEC")
0116 id2 = 2;
0117 else if (name2 == "CAHC")
0118 id2 = 3;
0119 else if (name2 == "CEDC")
0120 id2 = 4;
0121 else if (name2 == "CHDC")
0122 id2 = 5;
0123 else
0124 id2 = 0;
0125 } else {
0126 id1 = id2 = 0;
0127 }
0128 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos: Level " << level << " Volume " << name1 << " and " << name2
0129 << " ID1 " << id1 << " (" << id1Old << ") ID2 " << id2 << " (" << id2Old << ")";
0130
0131 if (fillHistos) {
0132 if (id1 != id1Old) {
0133 if (id1 == 0)
0134 fillHisto(id1Old, 1);
0135 else
0136 fillHisto(id1, 0);
0137 }
0138 if (id2 != id2Old) {
0139 if (id2 == 0)
0140 fillHisto(id2Old, 1);
0141 else
0142 fillHisto(id2, 0);
0143 }
0144 }
0145
0146 stepLen += step;
0147 radLen += step / radl;
0148 intLen += step / intl;
0149 }
0150
0151 void MaterialBudgetCastorHistos::fillEndTrack() {
0152 if (fillHistos) {
0153 if (id1 != 0)
0154 fillHisto(id1, 1);
0155 if (id2 != 0)
0156 fillHisto(id2, 1);
0157 }
0158 if (printSum) {
0159 for (unsigned int ii = 0; ii < matList.size(); ii++) {
0160 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: " << matList[ii] << "\t" << stepLength[ii] << "\t"
0161 << radLength[ii] << "\t" << intLength[ii];
0162 }
0163 }
0164 }
0165
0166 void MaterialBudgetCastorHistos::book() {
0167
0168 edm::Service<TFileService> tfile;
0169
0170 if (!tfile.isAvailable())
0171 throw cms::Exception("BadConfig") << "MaterialBudgetCastorHistos: TFileService unavailable: "
0172 << "please add it to config file";
0173
0174 double maxPhi = pi;
0175 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
0176 << "histos === with " << binEta << " bins "
0177 << "in eta from " << etaLow << " to " << etaHigh << " and " << binPhi << " bins "
0178 << "in phi from " << -maxPhi << " to " << maxPhi;
0179
0180 std::string tag1, tag2;
0181
0182 for (int i = 0; i < maxSet; i++) {
0183 double minEta = etaLow;
0184 double maxEta = etaHigh;
0185 int ireg = i;
0186 if (i > 9) {
0187 minEta = -etaHigh;
0188 maxEta = -etaLow;
0189 ireg -= 10;
0190 tag2 = " (-ve Eta Side)";
0191 } else {
0192 tag2 = " (+ve Eta Side)";
0193 }
0194 if ((i % 2) == 0) {
0195 ireg /= 2;
0196 tag1 = " == Start";
0197 } else {
0198 ireg = (ireg - 1) / 2;
0199 tag1 = " == End";
0200 }
0201 std::string title = std::to_string(ireg) + tag1 + tag2;
0202 me100[i] = tfile->make<TProfile>(
0203 std::to_string(i + 100).c_str(), ("MB(X0) prof Eta in region " + title).c_str(), binEta, minEta, maxEta);
0204 me200[i] = tfile->make<TProfile>(
0205 std::to_string(i + 200).c_str(), ("MB(L0) prof Eta in region " + title).c_str(), binEta, minEta, maxEta);
0206 me300[i] = tfile->make<TProfile>(
0207 std::to_string(i + 300).c_str(), ("MB(Step) prof Eta in region " + title).c_str(), binEta, minEta, maxEta);
0208 me400[i] =
0209 tfile->make<TH1F>(std::to_string(i + 400).c_str(), ("Eta in region " + title).c_str(), binEta, minEta, maxEta);
0210 me500[i] = tfile->make<TProfile>(
0211 std::to_string(i + 500).c_str(), ("MB(X0) prof Ph in region " + title).c_str(), binPhi, -maxPhi, maxPhi);
0212 me600[i] = tfile->make<TProfile>(
0213 std::to_string(i + 600).c_str(), ("MB(L0) prof Ph in region " + title).c_str(), binPhi, -maxPhi, maxPhi);
0214 me700[i] = tfile->make<TProfile>(
0215 std::to_string(i + 700).c_str(), ("MB(Step) prof Ph in region " + title).c_str(), binPhi, -maxPhi, maxPhi);
0216 me800[i] =
0217 tfile->make<TH1F>(std::to_string(i + 800).c_str(), ("Phi in region " + title).c_str(), binPhi, -maxPhi, maxPhi);
0218 me900[i] = tfile->make<TProfile2D>(std::to_string(i + 900).c_str(),
0219 ("MB(X0) prof Eta Phi in region " + title).c_str(),
0220 binEta / 2,
0221 minEta,
0222 maxEta,
0223 binPhi / 2,
0224 -maxPhi,
0225 maxPhi);
0226 me1000[i] = tfile->make<TProfile2D>(std::to_string(i + 1000).c_str(),
0227 ("MB(L0) prof Eta Phi in region " + title).c_str(),
0228 binEta / 2,
0229 minEta,
0230 maxEta,
0231 binPhi / 2,
0232 -maxPhi,
0233 maxPhi);
0234 me1100[i] = tfile->make<TProfile2D>(std::to_string(i + 1100).c_str(),
0235 ("MB(Step) prof Eta Phi in region " + title).c_str(),
0236 binEta / 2,
0237 minEta,
0238 maxEta,
0239 binPhi / 2,
0240 -maxPhi,
0241 maxPhi);
0242 me1200[i] = tfile->make<TH2F>(std::to_string(i + 1200).c_str(),
0243 ("Eta vs Phi in region " + title).c_str(),
0244 binEta / 2,
0245 minEta,
0246 maxEta,
0247 binPhi / 2,
0248 -maxPhi,
0249 maxPhi);
0250 }
0251
0252 edm::LogInfo("MaterialBudget") << "MaterialBudgetCastorHistos: Booking user "
0253 << "histos done ===";
0254 }
0255
0256 void MaterialBudgetCastorHistos::fillHisto(int id, int ix) {
0257 int ii = 2 * (id - 1) + ix;
0258 double etaAbs = eta_;
0259 if (eta_ < 0) {
0260 etaAbs = -eta_;
0261 ii += 10;
0262 }
0263 LogDebug("MaterialBudget") << "MaterialBudgetCastorHistos:FillHisto "
0264 << "called with index " << id << ":" << ix << ":" << ii << " eta " << etaAbs << " ("
0265 << eta_ << ") integrated step " << stepLen << " X0 " << radLen << " Lamda " << intLen;
0266
0267 me100[ii]->Fill(eta_, radLen);
0268 me200[ii]->Fill(eta_, intLen);
0269 me300[ii]->Fill(eta_, stepLen);
0270 me400[ii]->Fill(eta_);
0271
0272 if (etaAbs >= etaLow && etaAbs <= etaHigh) {
0273 me500[ii]->Fill(phi_, radLen);
0274 me600[ii]->Fill(phi_, intLen);
0275 me700[ii]->Fill(phi_, stepLen);
0276 me800[ii]->Fill(phi_);
0277 }
0278
0279 me900[ii]->Fill(eta_, phi_, radLen);
0280 me1000[ii]->Fill(eta_, phi_, intLen);
0281 me1100[ii]->Fill(eta_, phi_, stepLen);
0282 me1200[ii]->Fill(eta_, phi_);
0283 }
