TPC_MC_anchoring: Added per region relative gas gain to simulate the change in the voltage settings of GEMs

Author: tubagundem

Detectors/TPC/base/include/TPCBase/ParameterGEM.h

@@ -54,6 +54,7 @@ struct ParameterGEM : public o2::conf::ConfigurableParamHelper<ParameterGEM> {
   float AbsoluteGain[4] = {14.f, 8.f, 53.f, 240.f};              ///< Absolute gain
   float CollectionEfficiency[4] = {1.f, 0.2f, 0.25f, 1.f};       ///< Collection efficiency
   float ExtractionEfficiency[4] = {0.65f, 0.55f, 0.12f, 0.6f};   ///< Extraction efficiency
+  float RelativeGainStack[4] = {1.f, 1.f, 1.f, 1.f};             ///< Relative gain of the stack per region (IROC, OROC1, OROC2, OROC3) for the EffectiveMode
   float TotalGainStack = 2000.f;                                 ///< Total gain of the stack for the EffectiveMode
   float KappaStack = 1.205f;                                     ///< Variable steering the energy resolution of the full stack for the EffectiveMode
   float EfficiencyStack = 0.528f;                                ///< Variable steering the single electron efficiency of the full stack for the EffectiveMode

Detectors/TPC/simulation/include/TPCSimulation/GEMAmplification.h

@@ -55,9 +55,10 @@ class GEMAmplification
   int getStackAmplification(int nElectrons = 1);
 
   /// Compute the number of electrons after amplification in an effective single-stage amplification
+  /// \param region Readout region (0:IROC, 1:OROC1, 2:OROC2, 3:OROC3)
   /// \param nElectrons Number of electrons arriving at the first amplification stage (GEM1)
   /// \return Number of electrons after amplification in an  effective single-stage amplification
-  int getEffectiveStackAmplification(int nElectrons = 1);
+  int getEffectiveStackAmplification(int region, int nElectrons = 1);
 
   /// Compute the number of electrons after amplification in a full stack of four GEM foils
   /// taking into account local variations of the electron amplification
@@ -99,8 +100,8 @@ class GEMAmplification
   math_utils::RandomRing<> mRandomFlat;
   /// Container with random Polya distributions, one for each GEM in the stack
   std::array<math_utils::RandomRing<>, 4> mGain;
-  /// Container with random Polya distributions for the full stack amplification
-  math_utils::RandomRing<> mGainFullStack;
+  /// Container with random Polya distributions for the full stack amplification per region (IROC, OROC1, OROC2, OROC3)
+  std::array<math_utils::RandomRing<>, 4> mGainFullStack;
 
   const ParameterGEM* mGEMParam; ///< Caching of the parameter class to avoid multiple CDB calls
   const ParameterGas* mGasParam; ///< Caching of the parameter class to avoid multiple CDB calls
@@ -118,7 +119,7 @@ inline int GEMAmplification::getStackAmplification(const CRU& cru, const PadPos&
       break;
     }
     case AmplificationMode::EffectiveMode: {
-      return static_cast<int>(static_cast<float>(getEffectiveStackAmplification(nElectrons)) *
+      return static_cast<int>(static_cast<float>(getEffectiveStackAmplification(static_cast<int>(cru.gemStack()), nElectrons)) *
                               mGainMap->getValue(cru, pos.getRow(), pos.getPad()));
       break;
     }

Detectors/TPC/simulation/macro/toyCluster.C

@@ -446,7 +446,7 @@ void createDigitsFromSim(const char* inpFileSim = "o2sim_HitsTPC.root", const st
 
             // convert electrons to ADC signal
             const GlobalPadNumber globalPad = mapper.globalPadNumber(digiPadPos.getGlobalPadPos());
-            const float adcsignal = sampaProcessing.getADCvalue(static_cast<float>((nEleGEM == -1) ? static_cast<int>(gemAmplification.getEffectiveStackAmplification()) : nEleGEM));
+            const float adcsignal = sampaProcessing.getADCvalue(static_cast<float>((nEleGEM == -1) ? static_cast<int>(gemAmplification.getEffectiveStackAmplification(static_cast<int>(cru.gemStack()))) : nEleGEM));
             sampaProcessing.getShapedSignal(adcsignal, driftTime, signalArray);
 
             // set up MC label

Detectors/TPC/simulation/src/GEMAmplification.cxx

@@ -73,28 +73,31 @@ GEMAmplification::GEMAmplification()
     delete polyaDistribution;
   }
 
-  /// Set the polya distribution for the full stack
+  /// Set the polya distribution for the full stack per region (IROC, OROC1, OROC2, OROC3)
   const float gainStack = mGEMParam->TotalGainStack;
   const float kappaStack = mGEMParam->KappaStack;
   const float effStack = mGEMParam->EfficiencyStack;
   // Correct for electron losses occurring when changing kappa and the efficiency
-  const float sStack = gainStack / (kappaStack * effStack);
-  polya % kappaStack % sStack % sStack % (kappaStack - 1) % sStack;
-  std::string name = polya.str();
-  o2::math_utils::CachingTF1* polyaDistribution = nullptr;
-  if (!cacheexists) {
-    polyaDistribution = new o2::math_utils::CachingTF1("polya", name.c_str(), 0, 25.f * gainStack);
-    polyaDistribution->SetNpx(50000);
-  } else {
-    polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get("polyaStack");
-  }
-  mGainFullStack.initialize(*polyaDistribution);
+  for (int i = 0; i < 4; ++i) {
+    const float gainStackPerRegion = mGEMParam->RelativeGainStack[i] * gainStack;
+    LOG(info) << "TPC: Region: " << i << " Gain: " << gainStackPerRegion;
+    const float sStack = gainStackPerRegion / (kappaStack * effStack);
+    polya % kappaStack % sStack % sStack % (kappaStack - 1) % sStack;
+    std::string name = polya.str();
+    o2::math_utils::CachingTF1* polyaDistribution = nullptr;
+    if (!cacheexists) {
+      polyaDistribution = new o2::math_utils::CachingTF1("polya", name.c_str(), 0, 25.f * gainStackPerRegion);
+      polyaDistribution->SetNpx(50000);
+    } else {
+      polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get(TString::Format("polyaStack%d", i).Data());
+    }
+    mGainFullStack[i].initialize(*polyaDistribution);
 
-  if (!cacheexists) {
-    outfile->WriteTObject(polyaDistribution, "polyaStack");
+    if (!cacheexists) {
+      outfile->WriteTObject(polyaDistribution, TString::Format("polyaStack%d", i).Data());
+    }
+    delete polyaDistribution;
   }
-  delete polyaDistribution;
-
   if (outfile) {
     outfile->Close();
   }
@@ -122,17 +125,17 @@ int GEMAmplification::getStackAmplification(int nElectrons)
   return nElectronsGEM4;
 }
 
-int GEMAmplification::getEffectiveStackAmplification(int nElectrons)
+int GEMAmplification::getEffectiveStackAmplification(int region, int nElectrons)
 {
   /// We start with an arbitrary number of electrons given to the first amplification stage
   /// The amplification in the GEM stack is handled for each electron individually and the amplification
-  /// in the stack is handled in an effective manner
+  /// in the stack is handled in an effective manner for different regions (IROC, OROC1, OROC2, OROC3)
   int nElectronsGEM = 0;
   for (int i = 0; i < nElectrons; ++i) {
     if (mRandomFlat.getNextValue() > mGEMParam->EfficiencyStack) {
       continue;
     }
-    nElectronsGEM += mGainFullStack.getNextValue();
+    nElectronsGEM += mGainFullStack[region].getNextValue();
   }
   return nElectronsGEM;
 }

Detectors/TPC/simulation/test/testTPCGEMAmplification.cxx

@@ -73,7 +73,7 @@ BOOST_AUTO_TEST_CASE(GEMamplification_effective_test)
   const int nEleIn = 158; /// Number of electrons liberated in Ne-CO2-N2 by an incident Fe-55 photon
 
   for (int i = 0; i < 500000; ++i) {
-    hTest.Fill(gemStack.getEffectiveStackAmplification(nEleIn));
+    hTest.Fill(gemStack.getEffectiveStackAmplification(0, nEleIn));
   }
 
   hTest.Fit("gaus", "Q0");

Detectors/TPC/spacecharge/macro/createSCHistosFromHits.C

@@ -366,7 +366,7 @@ void createSCHistosFromHits(const int ionDriftTime = 200, const int nEvIon = 1,
             const auto pad = static_cast<size_t>(padPos.getPad());
             const CRU cru = digiPadPos.getCRU();
 
-            const int gain = static_cast<int>(gemAmplification.getEffectiveStackAmplification() * gainMap.getValue(cru, row, pad));
+            const int gain = static_cast<int>(gemAmplification.getEffectiveStackAmplification(static_cast<int>(cru.gemStack())) * gainMap.getValue(cru, row, pad));
             if (gain == 0) {
               continue;
             }
@@ -399,9 +399,9 @@ void createSCHistosFromHits(const int ionDriftTime = 200, const int nEvIon = 1,
             const float chargeEpsilon = vecCalCharge[zbinIDC].getValue(cru, row, pad) + epsilon;
             ((CalPadArr&)(vecCalCharge[zbinIDC].getCalArray(static_cast<size_t>(cru.roc().getRoc())))).setValue(rowRoc, pad, chargeEpsilon);
           } // electron loop
-        }   // hit loop
-      }     // track loop
-    }       // sector loop
+        } // hit loop
+      } // track loop
+    } // sector loop
 
     for (int isec = 0; isec < ::Sector::MAXSECTOR; ++isec) {
       delete arrSectors[isec];