Skip to content

Commit 5d4427f

Browse files
lobislobis
authored
Merge pull request #114 from rest-for-physics/mariajmz_diffusion
TRestDetectorElectronDiffusionProcess taking into account the fano factor
2 parents c1a601b + 11af6fb commit 5d4427f

File tree

2 files changed

+52
-32
lines changed

2 files changed

+52
-32
lines changed

Diff for: ‎inc/TRestDetectorElectronDiffusionProcess.h

+6-4
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,15 @@ class TRestDetectorElectronDiffusionProcess : public TRestEventProcess {
3939
Double_t fAttachment;
4040
Double_t fGasPressure;
4141
Double_t fWValue;
42+
Double_t fFanoFactor;
4243
Double_t fLongitudinalDiffusionCoefficient;
4344
Double_t fTransversalDiffusionCoefficient;
44-
Bool_t fPoissonElectronExcitation;
4545
Bool_t fUnitElectronEnergy;
4646
UInt_t fMaxHits;
47-
Double_t fSeed = 0;
47+
UInt_t fSeed = 0;
4848
Bool_t fCheckIsInside = true;
49+
Bool_t fUseFanoFactor = false;
50+
Bool_t fPoissonElectronExcitation = false;
4951

5052
public:
5153
RESTValue GetInputEvent() const override { return fInputHitsEvent; }
@@ -61,7 +63,7 @@ class TRestDetectorElectronDiffusionProcess : public TRestEventProcess {
6163
BeginPrintProcess();
6264

6365
RESTMetadata << " eField : " << fElectricField * units("V/cm") << " V/cm" << RESTendl;
64-
RESTMetadata << " attachment coeficient : " << fAttachment << " V/cm" << RESTendl;
66+
RESTMetadata << " attachment coefficient : " << fAttachment << " V/cm" << RESTendl;
6567
RESTMetadata << " gas pressure : " << fGasPressure << " atm" << RESTendl;
6668
RESTMetadata << " longitudinal diffusion coefficient : " << fLongitudinalDiffusionCoefficient
6769
<< " cm^1/2" << RESTendl;
@@ -90,7 +92,7 @@ class TRestDetectorElectronDiffusionProcess : public TRestEventProcess {
9092
// Destructor
9193
~TRestDetectorElectronDiffusionProcess();
9294

93-
ClassDefOverride(TRestDetectorElectronDiffusionProcess, 4); // Template for a REST "event process" class
95+
ClassDefOverride(TRestDetectorElectronDiffusionProcess, 5); // Template for a REST "event process" class
9496
// inherited from TRestEventProcess
9597
};
9698
#endif

Diff for: ‎src/TRestDetectorElectronDiffusionProcess.cxx

+46-28
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ void TRestDetectorElectronDiffusionProcess::Initialize() {
4949
fTransversalDiffusionCoefficient = 0;
5050
fLongitudinalDiffusionCoefficient = 0;
5151
fWValue = 0;
52+
fFanoFactor = 0;
5253

5354
fOutputHitsEvent = new TRestDetectorHitsEvent();
5455
fInputHitsEvent = nullptr;
@@ -96,14 +97,25 @@ void TRestDetectorElectronDiffusionProcess::InitProcess() {
9697
<< RESTendl;
9798
exit(1);
9899
#endif
99-
if (fGasPressure <= 0) fGasPressure = fGas->GetPressure();
100-
if (fElectricField <= 0) fElectricField = fGas->GetElectricField();
101-
if (fWValue <= 0) fWValue = fGas->GetWvalue();
102-
100+
if (fGasPressure <= 0) {
101+
fGasPressure = fGas->GetPressure();
102+
}
103103
fGas->SetPressure(fGasPressure);
104+
105+
if (fElectricField <= 0) {
106+
fElectricField = fGas->GetElectricField();
107+
}
104108
fGas->SetElectricField(fElectricField);
109+
110+
if (fWValue <= 0) {
111+
fWValue = fGas->GetWvalue();
112+
}
105113
fGas->SetW(fWValue);
106114

115+
if (fFanoFactor <= 0) {
116+
fFanoFactor = fGas->GetGasFanoFactor();
117+
}
118+
107119
if (fLongitudinalDiffusionCoefficient <= 0) {
108120
fLongitudinalDiffusionCoefficient = fGas->GetLongitudinalDiffusion();
109121
} // (cm)^1/2
@@ -124,8 +136,8 @@ TRestEvent* TRestDetectorElectronDiffusionProcess::ProcessEvent(TRestEvent* inpu
124136
fInputHitsEvent = (TRestDetectorHitsEvent*)inputEvent;
125137
fOutputHitsEvent->SetEventInfo(fInputHitsEvent);
126138

127-
set<unsigned int> hitsToProcess; // indices of the hits to process (we do not want to process veto hits)
128-
for (unsigned int n = 0; n < fInputHitsEvent->GetNumberOfHits(); n++) {
139+
set<int> hitsToProcess; // indices of the hits to process (we do not want to process veto hits)
140+
for (int n = 0; n < static_cast<int>(fInputHitsEvent->GetNumberOfHits()); n++) {
129141
if (fInputHitsEvent->GetType(n) == REST_HitType::VETO) {
130142
// keep unprocessed hits as they are
131143
fOutputHitsEvent->AddHit(fInputHitsEvent->GetX(n), fInputHitsEvent->GetY(n),
@@ -148,7 +160,7 @@ TRestEvent* TRestDetectorElectronDiffusionProcess::ProcessEvent(TRestEvent* inpu
148160
bool isAttached;
149161

150162
Double_t wValue = fWValue;
151-
const unsigned int totalElectrons = totalEnergy * REST_Units::eV / wValue;
163+
const auto totalElectrons = static_cast<unsigned int>(totalEnergy * REST_Units::eV / wValue);
152164

153165
// TODO: double check this
154166
if (fMaxHits > 0 && totalElectrons > fMaxHits) {
@@ -185,15 +197,21 @@ TRestEvent* TRestDetectorElectronDiffusionProcess::ProcessEvent(TRestEvent* inpu
185197

186198
Double_t driftDistance = plane->GetDistanceTo({x, y, z});
187199

188-
Int_t numberOfElectrons;
189-
if (fPoissonElectronExcitation) {
190-
numberOfElectrons = fRandom->Poisson(energy * REST_Units::eV / fWValue);
191-
if (wValue != fWValue) {
192-
// reduce the number of electrons to improve speed
193-
numberOfElectrons = numberOfElectrons * fWValue / wValue;
200+
double numberOfElectrons = energy * REST_Units::eV / wValue;
201+
if (fUseFanoFactor) {
202+
if (fFanoFactor <= 0) {
203+
throw runtime_error("Fano factor not valid");
194204
}
195-
} else {
196-
numberOfElectrons = energy * REST_Units::eV / wValue;
205+
const auto sigma = TMath::Sqrt(fFanoFactor * numberOfElectrons);
206+
numberOfElectrons = fRandom->Gaus(numberOfElectrons, sigma);
207+
} else if (fPoissonElectronExcitation) {
208+
// this is probably a bad idea, we only keep it for backwards compatibility
209+
numberOfElectrons = fRandom->Poisson(energy * REST_Units::eV / fWValue);
210+
}
211+
212+
if (wValue != fWValue) {
213+
// reduce the number of electrons to improve speed
214+
numberOfElectrons = static_cast<unsigned int>(numberOfElectrons * fWValue / wValue);
197215
}
198216

199217
if (numberOfElectrons <= 0) {
@@ -202,14 +220,12 @@ TRestEvent* TRestDetectorElectronDiffusionProcess::ProcessEvent(TRestEvent* inpu
202220

203221
const Double_t energyPerElectron = energy * REST_Units::eV / numberOfElectrons;
204222

205-
while (numberOfElectrons > 0) {
206-
numberOfElectrons--;
207-
208-
Double_t longitudinalDiffusion =
209-
10. * TMath::Sqrt(driftDistance / 10.) * fLongitudinalDiffusionCoefficient; // mm
210-
Double_t transversalDiffusion =
211-
10. * TMath::Sqrt(driftDistance / 10.) * fTransversalDiffusionCoefficient; // mm
223+
Double_t longitudinalDiffusion =
224+
10. * TMath::Sqrt(driftDistance / 10.) * fLongitudinalDiffusionCoefficient; // mm
225+
Double_t transversalDiffusion =
226+
10. * TMath::Sqrt(driftDistance / 10.) * fTransversalDiffusionCoefficient; // mm
212227

228+
for (unsigned int i = 0; i < numberOfElectrons; i++) {
213229
if (fAttachment > 0) {
214230
// TODO: where is this formula from?
215231
isAttached = (fRandom->Uniform(0, 1) > pow(1 - fAttachment, driftDistance / 10.));
@@ -240,7 +256,7 @@ TRestEvent* TRestDetectorElectronDiffusionProcess::ProcessEvent(TRestEvent* inpu
240256
1E-6 * plane->GetNormal().Z()); // add a delta to make sure readout finds it
241257
}
242258

243-
if (!fCheckIsInside && !plane->IsInside(positionAfterDiffusion)) {
259+
if (fCheckIsInside && !plane->IsInside(positionAfterDiffusion)) {
244260
// electron has been moved outside the readout plane
245261
continue;
246262
}
@@ -280,13 +296,14 @@ void TRestDetectorElectronDiffusionProcess::InitFromConfigFile() {
280296
fGasPressure = GetDblParameterWithUnits("gasPressure", -1.);
281297
fElectricField = GetDblParameterWithUnits("electricField", -1.);
282298
fWValue = GetDblParameterWithUnits("WValue", 0.0) * REST_Units::eV;
299+
fFanoFactor = GetDblParameterWithUnits("fanoFactor", 0.0);
283300
fAttachment = StringToDouble(GetParameter("attachment", "0"));
284301
fLongitudinalDiffusionCoefficient =
285302
StringToDouble(GetParameter("longitudinalDiffusionCoefficient", "-1"));
286303
if (fLongitudinalDiffusionCoefficient == -1)
287304
fLongitudinalDiffusionCoefficient = StringToDouble(GetParameter("longDiff", "-1"));
288305
else {
289-
RESTWarning << "longitudinalDiffusionCoefficient is now OBSOLETE! It will soon dissapear."
306+
RESTWarning << "longitudinalDiffusionCoefficient is now OBSOLETE! It will soon disappear."
290307
<< RESTendl;
291308
RESTWarning << " Please use the shorter form of this parameter : longDiff" << RESTendl;
292309
}
@@ -295,12 +312,13 @@ void TRestDetectorElectronDiffusionProcess::InitFromConfigFile() {
295312
if (fTransversalDiffusionCoefficient == -1)
296313
fTransversalDiffusionCoefficient = StringToDouble(GetParameter("transDiff", "-1"));
297314
else {
298-
RESTWarning << "transversalDiffusionCoefficient is now OBSOLETE! It will soon dissapear." << RESTendl;
315+
RESTWarning << "transversalDiffusionCoefficient is now OBSOLETE! It will soon disappear." << RESTendl;
299316
RESTWarning << " Please use the shorter form of this parameter : transDiff" << RESTendl;
300317
}
301318
fMaxHits = StringToInteger(GetParameter("maxHits", "1000"));
302-
fSeed = StringToDouble(GetParameter("seed", "0"));
303-
fPoissonElectronExcitation = StringToBool(GetParameter("poissonElectronExcitation", "false"));
319+
fSeed = static_cast<UInt_t>(StringToInteger(GetParameter("seed", "0")));
320+
fPoissonElectronExcitation = StringToBool(GetParameter("poissonElectronExcitation", "true"));
304321
fUnitElectronEnergy = StringToBool(GetParameter("unitElectronEnergy", "false"));
305-
fCheckIsInside = StringToBool(GetParameter("checkIsInside", "true"));
322+
fCheckIsInside = StringToBool(GetParameter("checkIsInside", "false"));
323+
fUseFanoFactor = StringToBool(GetParameter("useFano", "false"));
306324
}

0 commit comments

Comments
 (0)