bsunanda:Run2-TB11 Adding provision for estimating material budget in geometry description

SimG4CMS/HGCalTestBeam/interface/HGCalTBMB.h

@@ -0,0 +1,47 @@
+#ifndef SimG4CMS_HGCalTestBeam_HGCalTBMB_h
+#define SimG4CMS_HGCalTestBeam_HGCalTBMB_h
+
+#include "SimG4Core/Watcher/interface/SimWatcher.h"
+#include "SimG4Core/Notification/interface/Observer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+class BeginOfTrack;
+class G4Step;
+class EndOfTrack;
+class G4VTouchable;
+
+#include <TH1F.h>
+#include <string>
+#include <vector>
+
+class HGCalTBMB : public SimWatcher, 
+                  public Observer<const BeginOfTrack*>,
+                  public Observer<const G4Step*>,
+                  public Observer<const EndOfTrack*> {
+
+public:
+
+  HGCalTBMB(const edm::ParameterSet&);
+  virtual ~HGCalTBMB();
+
+private:
+
+  HGCalTBMB(const HGCalTBMB&);          // stop default
+  const HGCalTBMB& operator=(const HGCalTBMB&); // ...
+
+  void update(const BeginOfTrack*);
+  void update(const G4Step*);
+  void update(const EndOfTrack*);
+
+  bool stopAfter(const G4Step*);
+  int  findVolume(const G4VTouchable* touch, bool stop) const;
+
+  std::vector<std::string>      listNames_;
+  std::string                   stopName_;
+  double                        stopZ_;
+  unsigned int                  nList_;
+  std::vector<double>           radLen_, intLen_, stepLen_;
+  std::vector<TH1D*>            me100_, me200_, me300_;
+};
+
+#endif

SimG4CMS/HGCalTestBeam/plugins/BuildFile.xml

@@ -12,10 +12,13 @@
 <use   name="SimDataFormats/Track"/>
 <use   name="SimDataFormats/Vertex"/>
 <use   name="SimG4CMS/HGCalTestBeam"/>
+<use   name="SimG4Core/Notification"/>
+<use   name="SimG4Core/Watcher"/>
 <use   name="clhep"/>
 <use   name="rootphysics"/>
 <use   name="root"/>
 <use   name="hepmc"/>
+<flags   EDM_PLUGIN="1"/>
 <library   file="*.cc" name="SimG4CMSTestBeamPlugins">
-  <flags   EDM_PLUGIN="1"/>
 </library>
+

SimG4CMS/HGCalTestBeam/plugins/module.cc

@@ -1,7 +1,11 @@
 #include "SimG4CMS/HGCalTestBeam/interface/HGCalTB16SD01.h"
+#include "SimG4CMS/HGCalTestBeam/interface/HGCalTBMB.h"
 #include "SimG4Core/SensitiveDetector/interface/SensitiveDetectorPluginFactory.h"
+#include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
 #include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
 
 typedef HGCalTB16SD01 HGCalTB1601SensitiveDetector;
 DEFINE_SENSITIVEDETECTOR(HGCalTB1601SensitiveDetector);
+DEFINE_SIMWATCHER (HGCalTBMB);
 

SimG4CMS/HGCalTestBeam/python/HGCalTB161Module8V2XML_cfi.py

@@ -0,0 +1,18 @@
+import FWCore.ParameterSet.Config as cms
+
+XMLIdealGeometryESSource = cms.ESSource("XMLIdealGeometryESSource",
+    geomXMLFiles = cms.vstring('Geometry/CMSCommonData/data/materials.xml',
+        'Geometry/CMSCommonData/data/rotations.xml',
+        'Geometry/HGCalCommonData/data/TB161/cms.xml',
+        'Geometry/HGCalCommonData/data/TB161/hgcal.xml',
+        'Geometry/HGCalCommonData/data/TB161/8ModuleV2/hgcalEE.xml',
+        'Geometry/HGCalCommonData/data/v7/hgcalwafer.xml',
+        'Geometry/HGCalCommonData/data/TB161/hgcalBeam.xml',
+        'Geometry/HGCalCommonData/data/TB161/hgcalsense.xml',
+        'Geometry/HGCalCommonData/data/TB161/hgcProdCuts.xml',
+        'Geometry/HGCalCommonData/data/TB161/8Module/hgcalCons.xml'
+                               ),
+    rootNodeName = cms.string('cms:OCMS')
+)
+
+

SimG4CMS/HGCalTestBeam/python/hgcalTBMBCERN_cfi.py

@@ -0,0 +1,21 @@
+import FWCore.ParameterSet.Config as cms
+from SimG4Core.Configuration.SimG4Core_cff import *
+
+g4SimHits.Watchers = cms.VPSet(cms.PSet(
+        HGCalTBMB = cms.PSet(
+            DetectorNames = cms.vstring(
+                'HGCalBeamWChamb',  
+                'HGCalBeamS1',
+                'HGCalBeamS2',
+                'HGCalBeamS3',
+                'HGCalBeamS4',
+                'HGCalBeamHaloCounter',
+                'HGCalBeamCK3',
+                'HGCalBeamMuonCounter',
+                ),
+            MaximumZ = cms.double(9500.),
+            StopName = cms.string("HGCal"),
+            ),
+        type = cms.string('HGCalTBMB')
+        )
+                               )

SimG4CMS/HGCalTestBeam/python/hgcalTBMBFNAL_cfi.py

@@ -0,0 +1,20 @@
+import FWCore.ParameterSet.Config as cms
+from SimG4Core.Configuration.SimG4Core_cff import *
+
+g4SimHits.Watchers = cms.VPSet(cms.PSet(
+        HGCalTBMB = cms.PSet(
+            DetectorNames = cms.vstring(
+                'HGCMTS6SC3b',  
+                'HGCHeTube',
+                'HGCFeChamber',
+                'HGCScint1',
+                'HGCScint2',
+                'HGCFSiTrack',
+                'HGCAlPlate',
+                ),
+            MaximumZ = cms.double(200.),
+            StopName = cms.string("HGCal"),
+            ),
+        type = cms.string('HGCalTBMB')
+        )
+                               )

SimG4CMS/HGCalTestBeam/src/HGCalTBMB.cc

@@ -0,0 +1,163 @@
+#include "SimG4CMS/HGCalTestBeam/interface/HGCalTBMB.h"
+
+#include "FWCore/Utilities/interface/Exception.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+#include "SimG4Core/Notification/interface/BeginOfTrack.h"
+#include "SimG4Core/Notification/interface/EndOfTrack.h"
+
+#include "G4LogicalVolumeStore.hh"
+#include "G4Step.hh"
+#include "G4Track.hh"
+
+#include <iostream>
+
+#define DebugLog
+
+
+HGCalTBMB::HGCalTBMB(const edm::ParameterSet& p) {
+
+  edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("HGCalTBMB");
+  listNames_   = m_p.getParameter<std::vector<std::string> >("DetectorNames");
+  stopName_    = m_p.getParameter<std::string>("StopName");
+  stopZ_       = m_p.getParameter<double>("MaximumZ");
+  nList_       = listNames_.size();
+  edm::LogInfo("HGCSim") << "HGCalTBMB initialized for " << nList_ <<" volumes\n";
+  for (unsigned int k=0; k<nList_; ++k)
+    edm::LogInfo("HGCSim") << " [" << k << "] " << listNames_[k] << std::endl;
+  edm::LogInfo("HGCSim") << "Stop after " << stopZ_ << " or reaching volume "
+			 << stopName_ << std::endl;
+
+  edm::Service<TFileService> tfile;
+  if ( !tfile.isAvailable() )
+    throw cms::Exception("BadConfig") << "TFileService unavailable: "
+                                      << "please add it to config file";
+  char  name[20], title[80];
+  TH1D* hist;
+  for (unsigned int i=0; i<=nList_; i++) {
+    std::string named = (i == nList_) ? "Total" : listNames_[i];
+    sprintf(name, "RadL%d", i);
+    sprintf(title, "MB(X0) for (%s)", named.c_str());
+    hist = tfile->make<TH1D>(name,title,1000,0.0,100.0);
+    hist->Sumw2(true);
+    me100_.push_back(hist);
+    sprintf(name, "IntL%d", i);
+    sprintf(title, "MB(L0) for (%s)", named.c_str());
+    hist = tfile->make<TH1D>(name,title,1000,0.0,10.0);
+    hist->Sumw2(true);
+    me200_.push_back(hist);
+    sprintf(name, "StepL%d", i);
+    sprintf(title, "MB(Step) for (%s)", named.c_str());
+    hist = tfile->make<TH1D>(name,title,1000,0.0,50000.0);
+    hist->Sumw2(true);
+    me300_.push_back(hist);
+  }
+  edm::LogInfo("HGCSim") << "HGCalTBMB: Booking user histos done ===";
+}
+
+HGCalTBMB::~HGCalTBMB() { }
+
+
+void HGCalTBMB::update(const BeginOfTrack* trk) {
+
+  radLen_  = std::vector<double>(nList_+1,0);
+  intLen_  = std::vector<double>(nList_+1,0);
+  stepLen_ = std::vector<double>(nList_+1,0);
+
+#ifdef DebugLog
+  const G4Track *   aTrack = (*trk)(); // recover G4 pointer if wanted
+  const G4ThreeVector& mom = aTrack->GetMomentum() ;
+  double         theEnergy = aTrack->GetTotalEnergy();
+  int            theID     = (int)(aTrack->GetDefinition()->GetPDGEncoding());
+  std::cout << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID()
+	    << " Code " << theID << " Energy " <<theEnergy/CLHEP::GeV
+	    << " GeV; Momentum " << mom << std::endl;
+#endif
+}
+
+void HGCalTBMB::update(const G4Step* aStep) {
+
+  G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
+  double step    = aStep->GetStepLength();
+  double radl    = material->GetRadlen();
+  double intl    = material->GetNuclearInterLength();
+
+  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
+  int indx = findVolume(touch,false);
+
+  if (indx >= 0) {
+    stepLen_[indx] += step;
+    radLen_[indx]  += (step/radl);
+    intLen_[indx]  += (step/intl);
+  }
+  stepLen_[nList_] += step;
+  radLen_[nList_]  += (step/radl);
+  intLen_[nList_]  += (step/intl);
+#ifdef DebugLog
+  std::cout << "HGCalTBMB::Step in "
+	    << touch->GetVolume(0)->GetLogicalVolume()->GetName()
+	    << " Index " << indx <<" Step " << step << " RadL " << step/radl
+	    << " IntL " << step/intl << std::endl;
+#endif
+
+  if (stopAfter(aStep)) {
+    G4Track* track = aStep->GetTrack();
+    track->SetTrackStatus(fStopAndKill);
+  }
+}
+
+void HGCalTBMB::update(const EndOfTrack* trk) {
+
+  for (unsigned int ii=0; ii<=nList_; ++ii) {
+    me100_[ii]->Fill(radLen_[ii]);
+    me200_[ii]->Fill(intLen_[ii]);
+    me300_[ii]->Fill(stepLen_[ii]);
+#ifdef DebugLog
+    std::string name("Total");
+    if (ii < nList_) name = listNames_[ii];
+    std::cout << "HGCalTBMB::Volume[" << ii << "]: " << name  << " == Step "
+	      << stepLen_[ii] << " RadL " << radLen_[ii] << " IntL "
+	      << intLen_[ii] << std::endl;
+#endif
+  }
+}
+
+bool HGCalTBMB::stopAfter(const G4Step* aStep) {
+
+  bool   flag(false);
+  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
+  G4ThreeVector hitPoint    = aStep->GetPreStepPoint()->GetPosition();
+  if (aStep->GetPostStepPoint() != 0) 
+    hitPoint = aStep->GetPostStepPoint()->GetPosition();
+  double zz    = hitPoint.z();
+
+  if ((findVolume(touch,true) == 0) || (zz > stopZ_)) flag = true;
+#ifdef DebugLog
+  std::cout << " HGCalTBMB::Name " << touch->GetVolume(0)->GetName() << " z "
+	    << zz << " Flag" << flag << std::endl;
+#endif
+  return flag;
+}
+
+int HGCalTBMB::findVolume(const G4VTouchable* touch, bool stop) const {
+
+  int ivol  =-1;
+  int level = (touch->GetHistoryDepth())+1;
+  for (int ii = 0; ii < level; ii++) {
+    std::string name = touch->GetVolume(ii)->GetName();
+    if (stop) {
+      if (strcmp(name.c_str(),stopName_.c_str()) == 0) ivol = 0;
+    } else {
+      for (unsigned int k=0; k<nList_; ++k) {
+	if (strcmp(name.c_str(),listNames_[k].c_str()) == 0) {
+	  ivol = k; break;
+	}
+      }
+    }
+    if (ivol >= 0) break;
+  }
+  return ivol;
+}
+