KaliVeda
Toolkit for HIC analysis
KVLightEnergyCsI.cpp
1 #include "KVLightEnergyCsI.h"
2 #include "TMath.h"
3 #include "TH1.h"
4 #include "TFitResultPtr.h"
5 #include "TFitResult.h"
6 
8 
9 
10 //___________________________________________________________________________
11 
12 
13 
25 Double_t KVLightEnergyCsI::CalculLumiere(Double_t* x, Double_t* par)
26 {
27  //Calculate total light output for given Z, A and energy
28  //
29  //~~~~~~~~~~~~~~~~~~
30  // x[0] = energy (MeV)
31  // par[0] = a1 : gain factor
32  // par[1] = a2 : nuclear & recombination quenching term
33  // par[2] = a3 : threshold (MeV/u) for delta-ray production
34  // par[3] = a4 : fractional energy loss removed by delta rays
35  //~~~~~~~~~~~~~~~~~~
36  //
37 
38  Double_t energie = x[0];
39  Double_t c1 = par[0];
40  Double_t c2 = Z * Z * A * par[1];
41  Double_t c3 = A * par[2];
42  Double_t c4 = par[3];
43  Double_t T = 8 * A;
44  Double_t c4_new = c4 / (1. + TMath::Exp((c3 - energie) / T));
45  Double_t c0 = c4 / (1. + TMath::Exp(c3 / T));
46 
47  Double_t lumcalc = c1 * energie;
48  if (c2 > 0.0) {
49  Double_t xm = -c1 * c0 * c2 * TMath::Log(c2 / (c2 + c3));
50  lumcalc = lumcalc - c1 * c2 * TMath::Log(1. + energie / c2) + c1 * c2 * c4_new * TMath::Log((energie + c2) / (c3 + c2)) + xm;
51  }
52 
53  return lumcalc;
54 }
55 
56 
57 
65 
67 {
68  // Inverse function for fitting data of calculated CsI energy vs. measured light output
69  //
70  // x[0] = light
71  //
72  // Same parameters as CalculLumiere, we plug their values into the calibration function
73  // then return its inverted value
74 
76  return do_inversion(*x);
77 }
78 
79 
80 
83 
85 {
86  // \param[in] make_func used by child class constructors to inhibit creation of internal calibration function
87 
88  SetType("LightEnergyCsI");
89  if (make_func) SetCalibFunction(new TF1("fLight_CsI", this, &KVLightEnergyCsI::CalculLumiere, 0., 10000., 4));
91 }
92 
93 
94 
104 
106 {
107  // Calculate the calibrated energy (in MeV) for a given total light output.
108  //
109  // The Z and A of the particle should be given as extra parameters:
110  //
111  //~~~~~~~~~~~~~~~{.cpp}
112  // Compute(light, "Z=3,A=7");
113  //~~~~~~~~~~~~~~~
114  //
115 
116  if (!IsAvailableFor(z_and_a)) {
117  Error("Compute", "[%s] Cannot compute energy for : ", GetName());
118  z_and_a.ls();
119  return -1;
120  }
121  Z = z_and_a.GetIntValue("Z");
122  A = z_and_a.GetIntValue("A");
123  return KVCalibrator::Compute(light, z_and_a);
124 }
125 
126 
127 
139 
141 {
142  // Given the calibrated (or simulated) energy in MeV,
143  // calculate the corresponding total light output according to the
144  // calibration parameters (useful for filtering simulations).
145  //
146  // The Z and A of the particle should be given as extra parameters:
147  //
148  //~~~~~~~~~~~~~~~{.cpp}
149  // Invert(light, "Z=3,A=7");
150  //~~~~~~~~~~~~~~~
151  //
152 
153  if (!IsAvailableFor(z_and_a)) {
154  Error("Compute", "Cannot compute without knowing Z and A of nucleus");
155  return -1;
156  }
157  Z = z_and_a.GetIntValue("Z");
158  A = z_and_a.GetIntValue("A");
159  return KVCalibrator::Invert(energy, z_and_a);
160 }
161 
162 
163 
166 
168 {
169  // \return kFALSE if parameter list does not contain "Z=..." and "A=..."
170  return z_and_a.HasIntParameter("Z") && z_and_a.HasIntParameter("A");
171 }
172 
173 
174 
184 
185 void KVLightEnergyCsI::Fit(TH1 *lite_vs_e, int z, int a, double a3, double a4)
186 {
187  // Allows to easily fit a TProfile of mean total light output versus CsI energy
188  // for a given isotope \f$(Z,A)\f$.
189  //
190  // Reasonable ranges for \f$a_{1}\f$ and \f$a_{2}\f$ parameters are set.
191  //
192  // \f$a_{3}\f$ is fixed to the value given by argument a3 (default=6); if a3<0 we let the parameter vary between 0 and (-a3)
193  //
194  // \f$a_{4}\f$ is fixed to the value given by argument a4 (default=0.385); if argument a4<0 we let it vary by +/-(-a4) around this value
195  Z=z; A=a;
196  GetCalibFunction()->SetParameters(10.,0.5,6.,0.385);
197  GetCalibFunction()->SetParLimits(0,0.,1000.);
198  GetCalibFunction()->SetParLimits(1,0.,1.);
199  if(a3<0)
200  GetCalibFunction()->SetParLimits(2,0,-a3);
201  else
203  if(a4<0)
204  GetCalibFunction()->SetParLimits(3,0.385+a4,0.385-a4);
205  else
207  auto fit_result = lite_vs_e->Fit(GetCalibFunction(),"EMGS");
208  std::cout << "chi**2/ndf=" << fit_result->Chi2()/fit_result->Ndf()
209  << " proba=" << fit_result->Prob()*100 << "%\n";
210  std::cout << GetCalibFunction()->GetParameter(0) << ","
211  << GetCalibFunction()->GetParameter(1) << ","
212  << GetCalibFunction()->GetParameter(2) << ","
213  << GetCalibFunction()->GetParameter(3) << std::endl;
214 }
215 
216 
217 
227 
228 void KVLightEnergyCsI::invFit(TH1*e_vs_lite, int z, int a, double a3, double a4)
229 {
230  // Allows to easily fit a TProfile of mean calculated CsI energy vs. measured total light output
231  // for a given isotope \f$(Z,A)\f$ using the inverted light-energy relation.
232  //
233  // Reasonable ranges for \f$a_{1}\f$ and \f$a_{2}\f$ parameters are set.
234  //
235  // \f$a_{3}\f$ is fixed to the value given by argument a3 (default=6); if a3<0 we let the parameter vary between 0 and (-a3)
236  //
237  // \f$a_{4}\f$ is fixed to the value given by argument a4 (default=0.385); if argument a4<0 we let it vary by +/-(-a4) around this value
238 
239  if(!inv_fit_func)
240  inv_fit_func = std::make_unique<TF1>("fInv_Light_CsI", this, &KVLightEnergyCsI::InvCalc, 0., 100000., 4);
241  Z=z; A=a;
242  inv_fit_func->SetParameters(10.,0.5,6.,0.385);
243  inv_fit_func->SetParLimits(0,0.,1000.);
244  inv_fit_func->SetParLimits(1,0.,1.);
245  if(a3<0)
246  inv_fit_func->SetParLimits(2,0,-a3);
247  else
248  inv_fit_func->FixParameter(2,a3);
249  if(a4<0)
250  inv_fit_func->SetParLimits(3,0.385+a4,0.385-a4);
251  else
252  inv_fit_func->FixParameter(3,a4);
253  auto fit_result = e_vs_lite->Fit(inv_fit_func.get(),"EMGS");
254  std::cout << "chi**2/ndf=" << fit_result->Chi2()/fit_result->Ndf()
255  << " proba=" << fit_result->Prob()*100 << "%\n";
256  std::cout << inv_fit_func->GetParameter(0) << ","
257  << inv_fit_func->GetParameter(1) << ","
258  << inv_fit_func->GetParameter(2) << ","
259  << inv_fit_func->GetParameter(3) << std::endl;
260 }
261 
262 
bool Bool_t
double Double_t
void Error(const char *method, const char *msgfmt,...) const override
Definition: KVBase.cpp:1650
virtual void SetType(const Char_t *str)
Definition: KVBase.h:172
Base class for all detector calibrations.
Definition: KVCalibrator.h:99
TF1 * GetCalibFunction() const
Definition: KVCalibrator.h:118
void SetCalibFunction(TF1 *f)
Definition: KVCalibrator.h:112
void SetUseInverseFunction(Bool_t yes=kTRUE)
Definition: KVCalibrator.h:241
virtual Double_t Invert(Double_t x, const KVNameValueList &) const
Definition: KVCalibrator.h:191
Double_t do_inversion(Double_t x) const
Definition: KVCalibrator.h:122
virtual Double_t Compute(Double_t x, const KVNameValueList &) const
Definition: KVCalibrator.h:185
Light-energy calibration function for CsI detectors using a Fermi-function dependence on energy for d...
KVLightEnergyCsI(Bool_t make_func=kTRUE)
std::unique_ptr< TF1 > inv_fit_func
Double_t CalculLumiere(Double_t *, Double_t *)
Double_t Compute(Double_t chan, const KVNameValueList &z_and_a="") const override
void invFit(TH1 *, int z, int a, double a3=6.0, double a4=0.385)
Double_t Invert(Double_t, const KVNameValueList &z_and_a="") const override
Bool_t IsAvailableFor(const KVNameValueList &z_and_a) const override
Double_t InvCalc(Double_t *, Double_t *)
void Fit(TH1 *, int z, int a, double a3=6.0, double a4=0.385)
Handles lists of named parameters with different types, a list of KVNamedParameter objects.
Int_t GetIntValue(const Char_t *name) const
void ls(Option_t *opt="") const override
Bool_t HasIntParameter(const Char_t *name) const
double Chi2() const
virtual Double_t GetParameter(const TString &name) const
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
virtual void SetParameters(const Double_t *params)
virtual void FixParameter(Int_t ipar, Double_t value)
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
const char * GetName() const override
return c1
Double_t x[n]
return c2
return c3
double T(double x)
c0
Double_t Exp(Double_t x)
Double_t Log(Double_t x)
TArc a
ClassImp(TPyArg)