Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:32:36

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