Using get_E_hit function in TFCSLateralShapeParametrizationFluctChain and...

Using get_E_hit function in TFCSLateralShapeParametrizationFluctChain and TFCSLateralShapeParametrizationHitChain

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h

@@ -21,12 +21,19 @@ public:
   bool Initialize(TH1* hist);
   
   virtual void Print(Option_t *option = "") const override;
+  virtual void setMinWeight(float minWeight) {m_minWeight=minWeight;}
+  virtual void setMaxWeight(float maxWeight) {m_maxWeight=maxWeight;}
+  virtual float getMinWeight() const override;
+  virtual float getMaxWeight() const override;
+  
 protected:
   /// Histogram to be used for the shape simulation
   /// The histogram x-axis should be in dR^2=deta^2+dphi^2
   TH1* m_hist{nullptr};
+  float m_minWeight{-1.};
+  float m_maxWeight{-1.};
 
-  ClassDefOverride(TFCSHistoLateralShapeWeight,1)  //TFCSHistoLateralShapeWeight
+  ClassDefOverride(TFCSHistoLateralShapeWeight,2)  //TFCSHistoLateralShapeWeight
 };
 
 #if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationFluctChain.h

@@ -15,6 +15,9 @@ public:
   // set and get the amount of hit energy fluctation around E/n for n hits
   void set_hit_RMS(float RMS);
   float get_hit_RMS() const {return m_RMS;};
+  
+  ///Get hit energy from layer energy and number of hits
+  virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
 
   virtual FCSReturnCode simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
 

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h

@@ -18,6 +18,13 @@ public:
 
   ///Call get_number_of_hits() only once per shower simulation, as it could be calculated with random numbers and give different results each time. Return a value of -1 if this instance can't determine
   virtual int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+  
+  ///Get hit energy from layer energy and number of hits
+  virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+  
+  ///Get minimum and maximum value of weight for hit energy reweighting
+  virtual float getMinWeight() const;
+  virtual float getMaxWeight() const;
 
   class Hit
   {

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitChain.h

@@ -33,9 +33,16 @@ public:
   
   /// Call get_number_of_hits() only once, as it could contain a random number
   virtual int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+  
+  ///Get hit energy from layer energy and number of hits
+  virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
 
   ///Give the effective size sigma^2 of the fluctuations that should be generated.
   virtual float get_sigma2_fluctuation(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+  
+  ///Get minimum and maximum value of weight for hit energy reweighting
+  virtual float getMinWeight() const;
+  virtual float getMaxWeight() const;
 
   static constexpr float s_max_sigma2_fluctuation=1000;//! Do not persistify!
 

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx

@@ -18,14 +18,23 @@
 
 TFCSHistoLateralShapeWeight::TFCSHistoLateralShapeWeight(const char* name, const char* title) :
   TFCSLateralShapeParametrizationHitBase(name,title)
-{
-}
+{}
 
 TFCSHistoLateralShapeWeight::~TFCSHistoLateralShapeWeight()
 {
   if(m_hist) delete m_hist;
 }
 
+float TFCSHistoLateralShapeWeight::getMinWeight() const
+{
+  return m_minWeight;
+}
+
+float TFCSHistoLateralShapeWeight::getMaxWeight() const
+{
+  return m_maxWeight;
+}
+
 FCSReturnCode TFCSHistoLateralShapeWeight::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
 {
   if (!simulstate.randomEngine()) {
@@ -82,7 +91,7 @@ void TFCSHistoLateralShapeWeight::Print(Option_t *option) const
   TFCSLateralShapeParametrizationHitBase::Print(option);
 
   if(longprint) {
-    if(m_hist) ATH_MSG_INFO(optprint <<"  Histogram: "<<m_hist->GetNbinsX()<<" bins ["<<m_hist->GetXaxis()->GetXmin()<<","<<m_hist->GetXaxis()->GetXmax()<<"]");
+    if(m_hist) ATH_MSG_INFO(optprint <<"  Histogram: "<<m_hist->GetNbinsX()<<" bins ["<<m_hist->GetXaxis()->GetXmin()<<","<<m_hist->GetXaxis()->GetXmax()<<"]" << " min weight: " << m_minWeight << " max weight: " << m_maxWeight);
      else ATH_MSG_INFO(optprint <<"  no Histogram");
   }  
 }

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx

@@ -55,20 +55,37 @@ FCSReturnCode TFCSHistoLateralShapeWeightHitAndMiss::simulate_hit(Hit& hit,TFCSS
   if(RMS>0) {
     weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), meanweight, RMS);
   }  
-  if(meanweight<=1) {
-    //if meanweight<=1, give lower energy to hit. 
-    //TFCSLateralShapeParametrizationHitChain needs to be able to generate more hits in this case
-    hit.E()*=weight;
-  } else {
-    //if meanweight>1, accept only 1/meanweight events, but give them a higher energy increased by weight. 
-    //This leads to larger fluctuations, while keeping the shape unchanged.
-    float prob=1.0/meanweight;
+
+  /* -------------------------------------------------------------------
+   * Weight is used to scale hit energy.
+   * 
+   * if (meanweight > m_minWeight and meanweight < m_maxWeight)
+   * 	Hit is accecpted with probability of m_minWeight/meanweight.
+   * 	If not accepted, weight is set to zero (this is
+   *    equivalent to not accept the hit). 
+   * 
+   * if (meanweight<=m_minWeight)
+   * 	Give lower energy to hit with probability 1.  
+   * 	TFCSLateralShapeParametrizationHitChain needs to be able 
+   * 	to generate more hits in this case
+   * 
+   * if (meanweight >= m_maxWeight) 
+   * 	meanweight is above upper threshold. It is set to m_maxWeight.
+   * -------------------------------------------------------------------
+  */
+  if (meanweight >= m_maxWeight) {
+    weight = m_maxWeight;
+  }
+  else if (meanweight > m_minWeight){
+    float prob=m_minWeight/meanweight;
     float rnd=CLHEP::RandFlat::shoot(simulstate.randomEngine());
-    if(rnd<prob) hit.E()*=weight;
-     else hit.E()=0;
+    if(rnd>=prob) weight = 0.;
   }  
 
+
+  hit.E()*=weight;
   ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" meanweight="<<meanweight<<" weight="<<weight);
+  
   return FCSSuccess;
 }
 

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx

@@ -21,6 +21,18 @@ TFCSLateralShapeParametrizationFluctChain::TFCSLateralShapeParametrizationFluctC
 {
 }
 
+float TFCSLateralShapeParametrizationFluctChain::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{
+  const float sigma2 = get_sigma2_fluctuation(simulstate,truth,extrapol);
+  const int sample = calosample();
+  if(sigma2<=0 || sample<0) return -1.;
+  const float maxWeight = getMaxWeight();
+  
+  if(maxWeight>0) return simulstate.E(sample)*sigma2/maxWeight;
+  else return simulstate.E(sample)*sigma2;
+  
+}
+
 FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
   // Call get_number_of_hits() only once, as it could contain a random number
@@ -34,10 +46,9 @@ FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulation
   if(sigma2<1e-8) sigma2=1e-8; 
 
   const float Elayer=simulstate.E(calosample());
-  const int nhit = 1.0 / sigma2;
 
   //Make a good guess of the needed hit energy, assuming all hits would have the same energy
-  const float Eavghit=Elayer/nhit;
+  const float Eavghit=get_E_hit(simulstate,truth,extrapol);
   const float Eavghit_tenth=Eavghit/10;
   float sumEhit=0;
   float error2_sumEhit=0;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitBase.cxx

@@ -26,6 +26,24 @@ int TFCSLateralShapeParametrizationHitBase::get_number_of_hits(TFCSSimulationSta
   return -1;
 }
 
+float TFCSLateralShapeParametrizationHitBase::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{  
+  const int nhits = get_number_of_hits(simulstate,truth,extrapol);
+  const int sample = calosample();
+  if(nhits<=0 || sample<0) return -1.;
+  else return simulstate.E(sample)/nhits;
+}
+
+float TFCSLateralShapeParametrizationHitBase::getMinWeight() const
+{
+  return -1.;
+}
+  
+float TFCSLateralShapeParametrizationHitBase::getMaxWeight() const
+{
+  return -1.;
+}
+
 FCSReturnCode TFCSLateralShapeParametrizationHitBase::simulate_hit(Hit& hit,TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
 {
   int cs=calosample();

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx

@@ -43,6 +43,36 @@ int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationSt
   return 1;
 }
 
+float TFCSLateralShapeParametrizationHitChain::getMinWeight() const
+{
+  for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
+    float weight=hitsim->getMinWeight();
+    if(weight>0.) return weight;
+  } 
+  return -1.;
+}
+  
+float TFCSLateralShapeParametrizationHitChain::getMaxWeight() const
+{
+  for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
+    float weight=hitsim->getMaxWeight();
+    if(weight>0.) return weight;
+  }
+  return -1.;
+}
+
+
+float TFCSLateralShapeParametrizationHitChain::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{  
+  const int nhits = get_number_of_hits(simulstate,truth,extrapol);
+  const int sample = calosample();
+  if(nhits<=0 || sample<0) return -1.;
+  const float maxWeight = getMaxWeight();// maxWeight = -1 if  shapeWeight class is not in m_chain  
+  
+  if(maxWeight>0) return simulstate.E(sample)/(maxWeight*nhits); // maxWeight is used only if shapeWeight class is in m_chain
+  else return simulstate.E(sample)/nhits; // Otherwise, old definition of E_hit is used
+}
+
 float TFCSLateralShapeParametrizationHitChain::get_sigma2_fluctuation(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
   if(m_number_of_hits_simul) {
@@ -67,7 +97,7 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
   }
 
   const float Elayer=simulstate.E(calosample());
-  const float Ehit=Elayer/nhit;
+  const float Ehit=get_E_hit(simulstate,truth,extrapol);
   float sumEhit=0;
 
   const bool debug = msgLvl(MSG::DEBUG);
@@ -107,8 +137,12 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
     }
     sumEhit+=hit.E();
     ++ihit;
-    if(ihit>100*nhit && ihit>1000) {
-      ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): aborting hit chain, iterated " << 100*nhit << " times, expected " << nhit<<" times. Deposited E("<<calosample()<<")="<<sumEhit<<" expected E="<<Elayer);
+    
+    if( (ihit==10*nhit) || (ihit==100*nhit) ) {
+      ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Iterated " << ihit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit);
+    }                                                                                                                         
+    if(ihit>1000*nhit && ihit>1000) {
+      ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Aborting hit chain, iterated " << 1000*nhit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit<<" expected E="<<Elayer);
       break;
     }  
   } while (sumEhit<Elayer);