Merge pull request #125 from jzennamo/feature/zennamo_ellipsoidrecomb

Adding An Ellipsoidal Model for Recombination
LArSoft · Nov 8, 2023 · 19c24e1 · 19c24e1
2 parents 6c2fee1 + 32ab92a
commit 19c24e1
Show file tree

Hide file tree

Showing 7 changed files with 167 additions and 53 deletions.
diff --git a/larsim/IonizationScintillation/ISCalcCorrelated.cxx b/larsim/IonizationScintillation/ISCalcCorrelated.cxx
@@ -51,7 +51,11 @@ namespace larg4 {
     fRecombk = LArG4PropHandle->Recombk() / detProp.Density(detProp.Temperature());
     fModBoxA = LArG4PropHandle->ModBoxA();
     fModBoxB = LArG4PropHandle->ModBoxB() / detProp.Density(detProp.Temperature());
+    fEllipsModBoxA = LArG4PropHandle->EllipsModBoxA();
+    fEllipsModBoxB = LArG4PropHandle->EllipsModBoxB() / detProp.Density(detProp.Temperature());
+    fEllipsModBoxR = LArG4PropHandle->EllipsModBoxR();
     fUseModBoxRecomb = (bool)LArG4PropHandle->UseModBoxRecomb();
+    fUseEllipsModBoxRecomb = (bool)LArG4PropHandle->UseEllipsModBoxRecomb();
     fUseModLarqlRecomb = (bool)LArG4PropHandle->UseModLarqlRecomb();
     fUseBinomialFlucts = (bool)LArG4PropHandle->UseBinomialFlucts();
     fLarqlChi0A = LArG4PropHandle->LarqlChi0A();
@@ -96,6 +100,22 @@ namespace larg4 {
         double Xi = fModBoxB * dEdx / EFieldStep;
         recomb = std::log(fModBoxA + Xi) / Xi;
       }
+      else if (fUseEllipsModBoxRecomb) {
+
+        double phi = AngleToEFieldAtStep(detProp.Efield(), edep);
+
+        if (std::isnan(phi)) {
+          double Xi = fModBoxB * dEdx / EFieldStep;
+          recomb = std::log(fModBoxA + Xi) / Xi;
+        }
+        else {
+          double B_ellips =
+            fEllipsModBoxB * dEdx /
+            (EFieldStep * std::hypot(std::sin(phi), std::cos(phi) / fEllipsModBoxR));
+
+          recomb = std::log(fEllipsModBoxA + B_ellips) / B_ellips;
+        }
+      }
       // ... or using Birks/Doke
       else {
         recomb = fRecombA / (1. + dEdx * fRecombk / EFieldStep);
@@ -119,7 +139,7 @@ namespace larg4 {
       recomb = 1.;
     }
 
-    // using this recombination, calculate number of ionization electrons
+    // using this recombination, calculate number energy_deposit of ionization electrons
     if (num_ions > 0.)
       num_electrons =
         (fUseBinomialFlucts) ? fBinomialGen.fire(num_ions, recomb) : (num_ions * recomb);
@@ -152,6 +172,45 @@ namespace larg4 {
     auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
     return efield * std::hypot(1 + eFieldOffsets.X(), eFieldOffsets.Y(), eFieldOffsets.Z());
   }
+  //----------------------------------------------------------------------------
+  double ISCalcCorrelated::AngleToEFieldAtStep(double efield, sim::SimEnergyDeposit const& edep)
+  {
+
+    // electric field outside active volume set to zero
+    if (!fISTPC.isScintInActiveVolume(edep.MidPoint())) return 0.;
+
+    TVector3 stepvec(
+      edep.StartX() - edep.EndX(), edep.StartY() - edep.EndY(), edep.StartZ() - edep.EndZ());
+
+    TVector3 elecvec;
+
+    art::ServiceHandle<geo::Geometry const> fGeometry;
+    geo::TPCID tpcid = fGeometry->PositionToTPCID(edep.MidPoint());
+    const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
+
+    if (tpcGeo.DetectDriftDirection() == 1) elecvec.SetXYZ(1, 0, 0);
+    if (tpcGeo.DetectDriftDirection() == -1) elecvec.SetXYZ(-1, 0, 0);
+    if (tpcGeo.DetectDriftDirection() == 2) elecvec.SetXYZ(0, 1, 0);
+    if (tpcGeo.DetectDriftDirection() == -2) elecvec.SetXYZ(0, -1, 0);
+    if (tpcGeo.DetectDriftDirection() == 3) elecvec.SetXYZ(0, 0, 1);
+    if (tpcGeo.DetectDriftDirection() == -3) elecvec.SetXYZ(0, 0, -1);
+
+    elecvec *= efield;
+
+    // electric field inside active volume
+    if (fSCE->EnableSimEfieldSCE()) {
+      auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
+      TVector3 scevec(
+        efield * eFieldOffsets.X(), efield * eFieldOffsets.Y(), efield * eFieldOffsets.Z());
+      elecvec += scevec;
+    }
+
+    double angle = std::acos(stepvec.Dot(elecvec) / (stepvec.Mag() * elecvec.Mag()));
+
+    if (angle > TMath::PiOver2()) { angle = abs(TMath::Pi() - angle); }
+
+    return angle;
+  }
 
   //----------------------------------------------------------------------------
   // LArQL chi0 function = fraction of escaping electrons

diff --git a/larsim/IonizationScintillation/ISCalcCorrelated.h b/larsim/IonizationScintillation/ISCalcCorrelated.h
@@ -38,6 +38,7 @@ namespace larg4 {
     double EFieldAtStep(double efield,
                         sim::SimEnergyDeposit const& edep)
       override; //value of field with any corrections for this step
+    double AngleToEFieldAtStep(double efield, sim::SimEnergyDeposit const& edep);
     ISCalcData CalcIonAndScint(detinfo::DetectorPropertiesData const& detProp,
                                sim::SimEnergyDeposit const& edep) override;
 
@@ -46,24 +47,28 @@ namespace larg4 {
     const spacecharge::SpaceCharge* fSCE;
     CLHEP::RandBinomial fBinomialGen;
 
-    double fGeVToElectrons;  ///< from LArG4Parameters service
-    double fWion;            ///< W_ion (23.6 eV) == 1/fGeVToElectrons
-    double fWph;             ///< from LArG4Parameters service
-    double fScintPreScale;   ///< scintillation pre-scaling factor from LArProperties service
-    double fRecombA;         ///< from LArG4Parameters service
-    double fRecombk;         ///< from LArG4Parameters service
-    double fModBoxA;         ///< from LArG4Parameters service
-    double fModBoxB;         ///< from LArG4Parameters service
-    double fLarqlChi0A;      ///< from LArG4Parameters service
-    double fLarqlChi0B;      ///< from LArG4Parameters service
-    double fLarqlChi0C;      ///< from LArG4Parameters service
-    double fLarqlChi0D;      ///< from LArG4Parameters service
-    double fLarqlAlpha;      ///< from LArG4Parameters service
-    double fLarqlBeta;       ///< from LArG4Parameters service
-    double fQAlpha;          ///< from LArG4Parameters service
-    bool fUseModBoxRecomb;   ///< from LArG4Parameters service
-    bool fUseModLarqlRecomb; ///< from LArG4Parameters service
-    bool fUseBinomialFlucts; ///< from LArG4Parameters service
+    double fGeVToElectrons;      ///< from LArG4Parameters service
+    double fWion;                ///< W_ion (23.6 eV) == 1/fGeVToElectrons
+    double fWph;                 ///< from LArG4Parameters service
+    double fScintPreScale;       ///< scintillation pre-scaling factor from LArProperties service
+    double fRecombA;             ///< from LArG4Parameters service
+    double fRecombk;             ///< from LArG4Parameters service
+    double fModBoxA;             ///< from LArG4Parameters service
+    double fModBoxB;             ///< from LArG4Parameters service
+    double fEllipsModBoxA;       ///< from LArG4Parameters service
+    double fEllipsModBoxB;       ///< from LArG4Parameters service
+    double fEllipsModBoxR;       ///< from LArG4Parameters service
+    double fLarqlChi0A;          ///< from LArG4Parameters service
+    double fLarqlChi0B;          ///< from LArG4Parameters service
+    double fLarqlChi0C;          ///< from LArG4Parameters service
+    double fLarqlChi0D;          ///< from LArG4Parameters service
+    double fLarqlAlpha;          ///< from LArG4Parameters service
+    double fLarqlBeta;           ///< from LArG4Parameters service
+    double fQAlpha;              ///< from LArG4Parameters service
+    bool fUseModBoxRecomb;       ///< from LArG4Parameters service
+    bool fUseEllipsModBoxRecomb; ///< from LArG4Parameters service
+    bool fUseModLarqlRecomb;     ///< from LArG4Parameters service
+    bool fUseBinomialFlucts;     ///< from LArG4Parameters service
 
     double EscapingEFraction(
       double const dEdx) const; //LArQL chi0 function = fraction of escaping electrons

diff --git a/larsim/LegacyLArG4/ISCalculationCorrelated.cxx b/larsim/LegacyLArG4/ISCalculationCorrelated.cxx
@@ -60,13 +60,18 @@ namespace larg4 {
     fRecombk = lgpHandle->Recombk() / density;
     fModBoxA = lgpHandle->ModBoxA();
     fModBoxB = lgpHandle->ModBoxB() / density;
+    fEllipsModBoxA = lgpHandle->EllipsModBoxA();
+    fEllipsModBoxB = lgpHandle->EllipsModBoxB() / detProp.Density(detProp.Temperature());
+    fEllipsModBoxR = lgpHandle->EllipsModBoxR();
     fLarqlChi0A = lgpHandle->LarqlChi0A();
     fLarqlChi0B = lgpHandle->LarqlChi0B();
     fLarqlChi0C = lgpHandle->LarqlChi0C();
     fLarqlChi0D = lgpHandle->LarqlChi0D();
     fLarqlAlpha = lgpHandle->LarqlAlpha();
     fLarqlBeta = lgpHandle->LarqlBeta();
     fUseModBoxRecomb = lgpHandle->UseModBoxRecomb();
+    fUseEllipsModBoxRecomb = (bool)lgpHandle->UseEllipsModBoxRecomb();
+
     fUseModLarqlRecomb = lgpHandle->UseModLarqlRecomb();
 
     // determine the step size using the voxel sizes
@@ -127,6 +132,25 @@ namespace larg4 {
       else
         recomb = 0;
     }
+    else if (fUseEllipsModBoxRecomb) {
+
+      double phi = std::acos(abs(step->GetPreStepPoint()->GetPosition().x() -
+                                 step->GetPostStepPoint()->GetPosition().x()) /
+                             step->GetStepLength());
+
+      if (phi > std::atan(1) * 2) { phi = std::atan(1) * 4 - phi; }
+
+      if (phi != phi) {
+        double Xi = fModBoxB * dEdx / EFieldStep;
+        recomb = std::log(fModBoxA + Xi) / Xi;
+      }
+      else {
+        double B_ellips = fEllipsModBoxB * dEdx /
+                          (EFieldStep * std::hypot(std::sin(phi), std::cos(phi) / fEllipsModBoxR));
+
+        recomb = std::log(fEllipsModBoxA + B_ellips) / B_ellips;
+      }
+    }
     else {
       recomb = fRecombA / (1. + dEdx * fRecombk / EFieldStep);
     }

diff --git a/larsim/LegacyLArG4/ISCalculationCorrelated.h b/larsim/LegacyLArG4/ISCalculationCorrelated.h
@@ -36,23 +36,27 @@ namespace larg4 {
     double StepSizeLimit() const { return fStepSize; }
 
   private:
-    double fStepSize;        ///< maximum step to take
-    double fEfield;          ///< value of electric field from LArProperties service
-    double fWion;            ///< W_ion (23.6 eV) == 1/fGeVToElectrons
-    double fWph;             ///< W_ph (19.5 eV)
-    double fScintPreScale;   ///< scintillation pre-scale from LArProperties service
-    double fRecombA;         ///< from LArG4Parameters service
-    double fRecombk;         ///< from LArG4Parameters service
-    double fModBoxA;         ///< from LArG4Parameters service
-    double fModBoxB;         ///< from LArG4Parameters service
-    double fLarqlChi0A;      ///< from LArG4Parameters service
-    double fLarqlChi0B;      ///< from LArG4Parameters service
-    double fLarqlChi0C;      ///< from LArG4Parameters service
-    double fLarqlChi0D;      ///< from LArG4Parameters service
-    double fLarqlAlpha;      ///< from LArG4Parameters service
-    double fLarqlBeta;       ///< from LArG4Parameters service
-    bool fUseModBoxRecomb;   ///< from LArG4Parameters service
-    bool fUseModLarqlRecomb; ///< from LArG4Parameters service
+    double fStepSize;            ///< maximum step to take
+    double fEfield;              ///< value of electric field from LArProperties service
+    double fWion;                ///< W_ion (23.6 eV) == 1/fGeVToElectrons
+    double fWph;                 ///< W_ph (19.5 eV)
+    double fScintPreScale;       ///< scintillation pre-scale from LArProperties service
+    double fRecombA;             ///< from LArG4Parameters service
+    double fRecombk;             ///< from LArG4Parameters service
+    double fModBoxA;             ///< from LArG4Parameters service
+    double fModBoxB;             ///< from LArG4Parameters service
+    double fLarqlChi0A;          ///< from LArG4Parameters service
+    double fLarqlChi0B;          ///< from LArG4Parameters service
+    double fLarqlChi0C;          ///< from LArG4Parameters service
+    double fLarqlChi0D;          ///< from LArG4Parameters service
+    double fLarqlAlpha;          ///< from LArG4Parameters service
+    double fLarqlBeta;           ///< from LArG4Parameters service
+    double fEllipsModBoxA;       ///< from LArG4Parameters service
+    double fEllipsModBoxB;       ///< from LArG4Parameters service
+    double fEllipsModBoxR;       ///< from LArG4Parameters service
+    bool fUseEllipsModBoxRecomb; ///< from LArG4Parameters service
+    bool fUseModBoxRecomb;       ///< from LArG4Parameters service
+    bool fUseModLarqlRecomb;     ///< from LArG4Parameters service
 
     double EscapingEFraction(
       double const dEdx); //LArQL chi0 function = fraction of escaping electrons

diff --git a/larsim/Simulation/LArG4Parameters.cc b/larsim/Simulation/LArG4Parameters.cc
@@ -46,6 +46,9 @@ namespace sim {
     , fRecombk{pset.get<double>("Recombk", util::kRecombk)}
     , fModBoxA{pset.get<double>("ModBoxA", util::kModBoxA)}
     , fModBoxB{pset.get<double>("ModBoxB", util::kModBoxB)}
+    , fEllipsModBoxA{pset.get<double>("EllipsModBoxA")}
+    , fEllipsModBoxB{pset.get<double>("EllipsModBoxB")}
+    , fEllipsModBoxR{pset.get<double>("EllipsModBoxR")}
     , fLarqlChi0A{pset.get<double>("LarqlChi0A")}
     , fLarqlChi0B{pset.get<double>("LarqlChi0B")}
     , fLarqlChi0C{pset.get<double>("LarqlChi0C")}
@@ -55,6 +58,7 @@ namespace sim {
     , fWph{pset.get<double>("Wph")}
     , fQAlpha{pset.get<double>("QAlpha")}
     , fUseModBoxRecomb{pset.get<bool>("UseModBoxRecomb")}
+    , fUseEllipsModBoxRecomb{pset.get<bool>("UseEllipsModBoxRecomb")}
     , fUseModLarqlRecomb{pset.get<bool>("UseModLarqlRecomb")}
     , fUseBinomialFlucts{pset.get<bool>("UseBinomialFlucts", true)}
     , fIonAndScintCalculator{pset.get<std::string>("IonAndScintCalculator", "Separate")}

diff --git a/larsim/Simulation/LArG4Parameters.h b/larsim/Simulation/LArG4Parameters.h
@@ -43,6 +43,9 @@ namespace sim {
     double Recombk() const { return fRecombk; }
     double ModBoxA() const { return fModBoxA; }
     double ModBoxB() const { return fModBoxB; }
+    double EllipsModBoxA() const { return fEllipsModBoxA; }
+    double EllipsModBoxB() const { return fEllipsModBoxB; }
+    double EllipsModBoxR() const { return fEllipsModBoxR; }
     double LarqlChi0A() const { return fLarqlChi0A; }
     double LarqlChi0B() const { return fLarqlChi0B; }
     double LarqlChi0C() const { return fLarqlChi0C; }
@@ -52,6 +55,7 @@ namespace sim {
     double Wph() const { return fWph; }
     double QAlpha() const { return fQAlpha; }
     bool UseModBoxRecomb() const { return fUseModBoxRecomb; }
+    bool UseEllipsModBoxRecomb() const { return fUseEllipsModBoxRecomb; }
     bool UseModLarqlRecomb() const { return fUseModLarqlRecomb; }
     bool UseBinomialFlucts() const { return fUseBinomialFlucts; }
     double GeVToElectrons() const { return util::kGeVToElectrons; }
@@ -113,20 +117,25 @@ namespace sim {
     bool const fDisableWireplanes;     ///< Turn of LAr sensitivity and remove charge
                                        ///< drift simulation - use for running pure optical sims
     std::vector<unsigned short int> const
-      fSkipWireSignalInTPCs;       ///< selective disabling of drift simulation
-    double const fRecombA;         ///< Possibly override the RecombA parameter
-    double const fRecombk;         ///< Possibly override the Recombk parameter
-    double const fModBoxA;         ///< Possibly override the ModBoxA parameter
-    double const fModBoxB;         ///< Possibly override the ModBoxB parameter
-    double const fLarqlChi0A;      ///< Possibly override the LarqlChi0A parameter
-    double const fLarqlChi0B;      ///< Possibly override the LarqlChi0B parameter
-    double const fLarqlChi0C;      ///< Possibly override the LarqlChi0C parameter
-    double const fLarqlChi0D;      ///< Possibly override the LarqlChi0D parameter
-    double const fLarqlAlpha;      ///< Possibly override the LarqlAlpha parameter
-    double const fLarqlBeta;       ///< Possibly override the LarqlBeta parameter
-    double const fWph;             ///< Possibly override the Wph parameter
-    double const fQAlpha;          ///< Possibly override the QAlpha parameter
-    bool const fUseModBoxRecomb;   ///< Use Modified Box model recombination instead of Birks
+      fSkipWireSignalInTPCs;     ///< selective disabling of drift simulation
+    double const fRecombA;       ///< Possibly override the RecombA parameter
+    double const fRecombk;       ///< Possibly override the Recombk parameter
+    double const fModBoxA;       ///< Possibly override the ModBoxA parameter
+    double const fModBoxB;       ///< Possibly override the ModBoxB parameter
+    double const fEllipsModBoxA; ///< Possibly override the EllipsModBoxA parameter
+    double const fEllipsModBoxB; ///< Possibly override the EllipsModBoxB parameter
+    double const fEllipsModBoxR; ///< Possibly override the EllipsModBoxR parameter
+    double const fLarqlChi0A;    ///< Possibly override the LarqlChi0A parameter
+    double const fLarqlChi0B;    ///< Possibly override the LarqlChi0B parameter
+    double const fLarqlChi0C;    ///< Possibly override the LarqlChi0C parameter
+    double const fLarqlChi0D;    ///< Possibly override the LarqlChi0D parameter
+    double const fLarqlAlpha;    ///< Possibly override the LarqlAlpha parameter
+    double const fLarqlBeta;     ///< Possibly override the LarqlBeta parameter
+    double const fWph;           ///< Possibly override the Wph parameter
+    double const fQAlpha;        ///< Possibly override the QAlpha parameter
+    bool const fUseModBoxRecomb; ///< Use Modified Box model recombination instead of Birks
+    bool const
+      fUseEllipsModBoxRecomb; ///< Use Ellipsoid Modified Box model recombination instead of Birks
     bool const fUseModLarqlRecomb; ///< Use LArQL model recombination correction (dependence on EF)
     bool const fUseBinomialFlucts; ///< Use binomial fluctuations in correlated method
     std::string const

diff --git a/larsim/Simulation/simulationservices.fcl b/larsim/Simulation/simulationservices.fcl
@@ -23,9 +23,10 @@ standard_largeantparameters:
  CosmogenicXSMNBiasOn:     0 # 0 is off. 1 works. 2 still in development.
  CosmogenicXSMNBiasFactor: 1 # Not more than 5-ish cuz of numerical instabilities.
  DisableWireplanes:        false #if set true, charge drift simulation does not run - used for optical sim jobs OR just when you don't wanna drift the e's.
- SkipWireSignalInTPCs:     []     # put here TPC id's which should not receive ionization electrons - used to simulate TPC geom volumes which are actually dead LAr volumes in protoDUNE
- UseModBoxRecomb:          true   # use Modified Box recombination instead of Birks
- UseModLarqlRecomb:        false   # use LArQL recombination corrections (dependence on EF)
+ SkipWireSignalInTPCs:     []    # put here TPC id's which should not receive ionization electrons - used to simulate TPC geom volumes which are actually dead LAr volumes in protoDUNE
+ UseModBoxRecomb:          true  # use Modified Box recombination instead of Birks
+ UseEllipsModBoxRecomb:    false # use Ellipsoid Modified Box recomination instead of Birks
+ UseModLarqlRecomb:        false # use LArQL recombination corrections (dependence on EF)
 
  #* Recombination factor coefficients come from Nucl.Instrum.Meth.A523:275-286,2004
  #* * @f$ dE/dx @f$ is given by the voxel energy deposition, but have to convert it to MeV/cm from GeV/voxel width
@@ -43,6 +44,14 @@ standard_largeantparameters:
  ModBoxA: 0.930
  ModBoxB: 0.212
 
+ #* Recombination factor coefficients come from Gray Putnam 
+ #* * @f$ dE/dx @f$ is given by the voxel energy deposition, but have to convert it to MeV/cm from GeV/voxel width
+ #* * electric field: @f$ E @f$ in kV/cm
+ #* * `EllipsModBoxB` needs to be scaled with the electric field.
+ EllipsModBoxA: 0.906 # +/- 0.008 
+ EllipsModBoxB: 0.203 # +/- 0.008
+ EllipsModBoxR: 1.25  # +/- 0.02
+
  #* Recombination factor coefficients for LArQL
  #* https://doi.org/10.1088/1748-0221/17/07/C07009
  #* * @f$ dE/dx @f$ is given by the energy deposition in MeV/cm