File indexing completed on 2023-03-17 11:24:59
0001 #include "SimG4Core/MagneticField/test/FieldStepWatcher.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0005
0006 #include "G4LogicalVolume.hh"
0007 #include "G4Step.hh"
0008 #include "G4TransportationManager.hh"
0009 #include "G4VPhysicalVolume.hh"
0010
0011 FieldStepWatcher::FieldStepWatcher(const edm::ParameterSet &p) {
0012 level = p.getUntrackedParameter<int>("DepthLevel", 2);
0013 outFile = p.getUntrackedParameter<std::string>("OutputFile", "field.root");
0014
0015 edm::LogInfo("FieldStepWatcher") << "FieldStepWatcher initialised to monitor"
0016 << " level " << level << " with o/p on " << outFile;
0017 dbe_ = edm::Service<DQMStore>().operator->();
0018 }
0019
0020 FieldStepWatcher::~FieldStepWatcher() {
0021 if (dbe_ && !outFile.empty())
0022 dbe_->save(outFile);
0023 }
0024
0025 void FieldStepWatcher::update(const BeginOfRun *) {
0026 G4VPhysicalVolume *pv =
0027 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0028 findTouch(pv, 0);
0029 lvnames.push_back("Not Found");
0030
0031 edm::LogInfo("FieldStepWatcher") << "FieldStepWatcher: Finds " << lvnames.size() << " different volumes"
0032 << " at level " << level;
0033 unsigned num = lvnames.size();
0034 steps.push_back(0);
0035 for (unsigned int i = 0; i < num; i++) {
0036 edm::LogInfo("FieldStepWatcher") << "FieldStepWatcher: lvnames[" << i << "] = " << lvnames[i];
0037 steps.push_back(0);
0038 }
0039
0040 if (dbe_) {
0041 char titl[60], name[20];
0042 for (unsigned int i = 0; i <= lvnames.size(); i++) {
0043 std::string lvname = "All";
0044 if (i != 0)
0045 lvname = lvnames[i - 1];
0046 sprintf(name, "Step%d", i);
0047 sprintf(titl, "Step Length in Volume %s", lvname.c_str());
0048 MonitorElement *m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0049 meStep.push_back(m1);
0050 sprintf(name, "Call%d", i);
0051 sprintf(titl, "Number of steps in Volume %s", lvname.c_str());
0052 MonitorElement *m2 = dbe_->book1D(name, titl, 50000, 0., 5000000.);
0053 meCall.push_back(m2);
0054 sprintf(name, "StepE%d", i);
0055 sprintf(titl, "Step Length for Electrons in Volume %s", lvname.c_str());
0056 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0057 meStepE.push_back(m1);
0058 sprintf(name, "StepG%d", i);
0059 sprintf(titl, "Step Length for Photons in Volume %s", lvname.c_str());
0060 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0061 meStepG.push_back(m1);
0062 sprintf(name, "StepMu%d", i);
0063 sprintf(titl, "Step Length for Muons in Volume %s", lvname.c_str());
0064 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0065 meStepMu.push_back(m1);
0066 sprintf(name, "StepNu%d", i);
0067 sprintf(titl, "Step Length for Neutrinos in Volume %s", lvname.c_str());
0068 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0069 meStepNu.push_back(m1);
0070 sprintf(name, "StepCH%d", i);
0071 sprintf(titl, "Step Length for Charged Hadrons in Volume %s", lvname.c_str());
0072 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0073 meStepCH.push_back(m1);
0074 sprintf(name, "StepNH%d", i);
0075 sprintf(titl, "Step Length for Neutral Hadrons in Volume %s", lvname.c_str());
0076 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0077 meStepNH.push_back(m1);
0078 sprintf(name, "StepC%d", i);
0079 sprintf(titl, "Step Length for Charged Particles in Volume %s", lvname.c_str());
0080 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0081 meStepC.push_back(m1);
0082 sprintf(name, "StepN%d", i);
0083 sprintf(titl, "Step Length for Neutral Particles in Volume %s", lvname.c_str());
0084 m1 = dbe_->book1D(name, titl, 5000, 0., 10000.);
0085 meStepN.push_back(m1);
0086 }
0087 }
0088 }
0089
0090 void FieldStepWatcher::update(const BeginOfEvent *) {
0091 for (unsigned int i = 0; i < steps.size(); i++)
0092 steps[i] = 0;
0093 }
0094
0095 void FieldStepWatcher::update(const EndOfEvent *) {
0096 if (dbe_) {
0097 for (unsigned int i = 0; i < steps.size(); i++)
0098 meCall[i]->Fill(steps[i]);
0099 }
0100 }
0101
0102 void FieldStepWatcher::update(const G4Step *aStep) {
0103 if (aStep) {
0104 G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0105 const G4VTouchable *pre_touch = preStepPoint->GetTouchable();
0106 int pre_level = ((pre_touch->GetHistoryDepth()) + 1);
0107 std::string name = "Not Found";
0108 if (pre_level > 0 && pre_level >= level) {
0109 int ii = pre_level - level;
0110 name = pre_touch->GetVolume(ii)->GetName();
0111 }
0112 int indx = findName(name);
0113 double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
0114 int code = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0115 LogDebug("FieldStepWatcher") << "FieldStepWatcher:: Step at Level " << pre_level << " with " << name << " at"
0116 << " level " << level << " corresponding"
0117 << " to index " << indx << " due to "
0118 << "particle " << code << " of charge " << charge;
0119 steps[0]++;
0120 if (indx >= 0) {
0121 int i = indx + 1;
0122 steps[indx + 1]++;
0123 if (dbe_) {
0124 meStep[0]->Fill(aStep->GetStepLength());
0125 meStep[i]->Fill(aStep->GetStepLength());
0126 if (code == 11 || code == -11) {
0127 meStepE[0]->Fill(aStep->GetStepLength());
0128 meStepE[i]->Fill(aStep->GetStepLength());
0129 } else if (code == 22) {
0130 meStepG[0]->Fill(aStep->GetStepLength());
0131 meStepG[i]->Fill(aStep->GetStepLength());
0132 } else if (code == 13 || code == -13) {
0133 meStepMu[0]->Fill(aStep->GetStepLength());
0134 meStepMu[i]->Fill(aStep->GetStepLength());
0135 } else if (code == 12 || code == -12 || code == 14 || code == -14 || code == 16 || code == -16) {
0136 meStepNu[0]->Fill(aStep->GetStepLength());
0137 meStepNu[i]->Fill(aStep->GetStepLength());
0138 } else if (charge == 0) {
0139 meStepNH[0]->Fill(aStep->GetStepLength());
0140 meStepNH[i]->Fill(aStep->GetStepLength());
0141 } else {
0142 meStepCH[0]->Fill(aStep->GetStepLength());
0143 meStepCH[i]->Fill(aStep->GetStepLength());
0144 }
0145 if (charge == 0) {
0146 meStepN[0]->Fill(aStep->GetStepLength());
0147 meStepN[i]->Fill(aStep->GetStepLength());
0148 } else {
0149 meStepC[0]->Fill(aStep->GetStepLength());
0150 meStepC[i]->Fill(aStep->GetStepLength());
0151 }
0152 }
0153 }
0154 }
0155 }
0156
0157 void FieldStepWatcher::findTouch(G4VPhysicalVolume *pv, int leafDepth) {
0158 if (leafDepth == 0)
0159 fHistory.SetFirstEntry(pv);
0160 else
0161 fHistory.NewLevel(pv, kNormal, pv->GetCopyNo());
0162
0163 G4LogicalVolume *lv = pv->GetLogicalVolume();
0164 LogDebug("FieldStepWatcher") << "FieldStepWatcher::find Touch " << lv->GetName() << " at level " << leafDepth;
0165 if (leafDepth == level - 1) {
0166 const std::string &lvname = lv->GetName();
0167 if (findName(lvname) < 0)
0168 lvnames.push_back(lvname);
0169 } else if (leafDepth < level - 1) {
0170 int noDaughters = lv->GetNoDaughters();
0171 while ((noDaughters--) > 0) {
0172 G4VPhysicalVolume *pvD = lv->GetDaughter(noDaughters);
0173 if (!pvD->IsReplicated())
0174 findTouch(pvD, leafDepth + 1);
0175 }
0176 }
0177 if (leafDepth > 0)
0178 fHistory.BackLevel();
0179 }
0180
0181 int FieldStepWatcher::findName(std::string name) {
0182 for (unsigned int i = 0; i < lvnames.size(); i++)
0183 if (name == lvnames[i])
0184 return i;
0185 return -1;
0186 }
0187
0188
0189 #include "FWCore/PluginManager/interface/ModuleDef.h"
0190 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0191
0192 DEFINE_SIMWATCHER(FieldStepWatcher);
