Merge pull request #222 from fucd/birks

SIM: Birks option to apply energy correction for scintillator

Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py

@@ -88,6 +88,10 @@ detsimalg.AnaElems = [
 ]
 detsimalg.RootDetElem = "WorldDetElemTool"
 
+from Configurables import CalorimeterSensDetTool
+cal_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+cal_sensdettool.CalNamesApplyBirks = ["EcalBarrel"]
+
 # output
 from Configurables import PodioOutput
 out = PodioOutput("outputalg")

Detector/DetCRD/src/Calorimeter/CRDEcal_v01.cpp

@@ -182,7 +182,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
   	 DetElement sd(ECAL, _toString(i,"trap%3d"), detid);
   	 sd.setPlacement(plv);
   }
-  
+
   sens.setType("calorimeter");
   ECAL.addExtension< LayeredCalorimeterData >( caloData ) ; 
 

Simulation/DetSimSD/include/DetSimSD/CaloSensitiveDetector.h

@@ -23,6 +23,7 @@ public:
     virtual void Initialize(G4HCofThisEvent* HCE);
     virtual G4bool ProcessHits(G4Step* step,G4TouchableHistory* history);
     virtual void EndOfEvent(G4HCofThisEvent* HCE);
+    void ApplyBirksLaw(){m_applyBirksLaw = true;};
 
 protected:
     CalorimeterHit* find(const HitCollection*, const dd4hep::sim::HitCompare<CalorimeterHit>&);
@@ -31,7 +32,8 @@ protected:
 
     HitCollection* m_hc;
     std::map<unsigned long, CalorimeterHit*> m_hitMap;
-    bool                                     m_isMergeEnabled;
+    bool                                     m_isMergeEnabled = false;
+    bool                                     m_applyBirksLaw  = false;
 };
 
 

Simulation/DetSimSD/src/CaloSensitiveDetector.cpp

@@ -33,6 +33,7 @@ CaloSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
     // std::cout << "CaloSensitiveDetector::ProcessHits" << std::endl;
 
     dd4hep::sim::Geant4StepHandler h(step);
+    if(m_applyBirksLaw) h.doApplyBirksLaw();
     dd4hep::Position pos = 0.5 * (h.prePos() + h.postPos());
     HitContribution contrib = dd4hep::sim::Geant4Hit::extractContribution(step);
     const std::string& name = GetName();
@@ -62,9 +63,9 @@ CaloSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
         m_hc->insert(hit);
     }
     hit->truth.push_back(contrib);
-    hit->energyDeposit += contrib.deposit;
-
-
+    //hit->energyDeposit += contrib.deposit;
+    hit->energyDeposit += h.totalEnergy();
+    //std::cout << "Apply Birk law: before = " << contrib.deposit << " after = " << h.totalEnergy() << std::endl;
     
     return true;
 }

Simulation/DetSimSD/src/CalorimeterSensDetTool.cpp

@@ -42,9 +42,19 @@ CalorimeterSensDetTool::createSD(const std::string& name) {
 	break;
       }
     }
-    G4VSensitiveDetector* sd = new CaloSensitiveDetector(name, *dd4hep_geo, is_merge_enabled);
+
+    CaloSensitiveDetector* sd = new CaloSensitiveDetector(name, *dd4hep_geo, is_merge_enabled);
     warning() << name << " set to merge true/false = " << is_merge_enabled << endmsg;
 
+    for(auto cal_name : m_listCalsApplyBirks){
+      if(cal_name==name){
+	info() << name << " will apply Birks law" << endmsg;
+	sd->ApplyBirksLaw();
+	break;
+      }
+    }
+
+
     return sd;
 }
 

Simulation/DetSimSD/src/CalorimeterSensDetTool.h

@@ -30,6 +30,7 @@ private:
     SmartIF<IGeomSvc> m_geosvc;
 
     Gaudi::Property<std::vector<std::string> > m_listCalsMergeDisable{this, "CalNamesMergeDisable", {}};
+    Gaudi::Property<std::vector<std::string> > m_listCalsApplyBirks{this, "CalNamesApplyBirks", {}};
 };
 
 #endif