Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:49

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 // define this as a plug-in
0189 #include "FWCore/PluginManager/interface/ModuleDef.h"
0190 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0191 
0192 DEFINE_SIMWATCHER(FieldStepWatcher);