er  dev
ERSupport.cxx
1 #include "ERSupport.h"
2 
3 #include <stdio.h>
4 #include <string.h>
5 #include <iostream>
6 
7 #include "G4IonTable.hh"
8 #include "G4EmCalculator.hh"
9 #include "G4NistManager.hh"
10 #include "G4ParticleDefinition.hh"
11 
12 #include "FairLogger.h"
13 
14 //--------------------------------------------------------------------------------------------------
15 TMatrixD* ReadCalFile(TString fileName) {
16  std::ifstream in;
17  in.open(fileName);
18  if (!in.is_open()){
19  LOG(FATAL) << "Can`t open calibration file " << fileName << FairLogger::FairLogger::endl;
20  return NULL;
21  }
22 
23  Int_t nRows = -1, nCols = -1;
24  in >> nCols;
25  in >> nRows;
26  if (nCols <= 0 || nRows <= 0){
27  LOG(FATAL) << "Can`t read rows or cols from calibration file " << fileName << FairLogger::FairLogger::endl;
28  return NULL;
29  }
30 
31  TMatrixD* calTable = new TMatrixD(nRows,nCols);
32  Int_t i = 0;
33 
34  while (!in.eof()){
35  if (i >= nRows){
36  break;
37  }
38  in >> (*calTable)[i][0] >> (*calTable)[i][1];
39  i++;
40  }
41 
42  if (i < nRows) {
43  LOG(FATAL) << "Wrong file format in " << fileName << FairLogger::FairLogger::endl;
44  }
45 
46  return calTable;
47 }
48 //-----------------------------------------------------------------------------
49 double EiEo(double tableER[][105],double Tp,double Rp){
50  if(Tp<0.1||Tp>1000.)
51  {printf("Energy is out of range\n"); return -1.;}
52 
53  int look=0;
54  while (Tp>tableER[0][look]) {look++;}
55 
56  double R1=tableER[2][look-1]+
57  +tableER[3][look-1]*Tp+
58  +tableER[4][look-1]*Tp*Tp+
59  +tableER[5][look-1]*Tp*Tp*Tp;
60 
61  double R2 = R1 - Rp;
62 
63  if(Rp>0.)
64  {if(R2<tableER[1][0]) {return 0.;}}
65  else
66  {if(R2>=tableER[1][104]) {return -2.;
67  printf("table is out of range\n");}}
68 
69  look=0;
70  while (R2>tableER[1][look]) {look++;}
71 
72  double E1=tableER[6][look-1]+
73  +tableER[7][look-1]*R2+
74  +tableER[8][look-1]*R2*R2+
75  +tableER[9][look-1]*R2*R2*R2;
76 
77  return E1;
78 };
79 //-----------------------------------------------------------------------------
80 int ReadRint(char* Fname,double Ranges[][105]){
81  char DummyA[256];
82  char termA[]="---";
83  char unR[10],un3[10],un4[10],dee[10],b[10],cc[10];
84  char unE[10]="begin";
85  double Energy,Rng,e,r,den;
86  double x[4],y[4],c[4];
87  int i;
88  double a[4];
89 
90  FILE *F2 = fopen(Fname,"r");
91  if(F2==NULL) {printf("ReadRint: File %s has not been read\n",Fname);
92  return 1;}
93  else
94  {
95  printf("ReadRint: File %s has been read\n",Fname);
96  while (fgets(DummyA,256,F2)) {if(strstr(DummyA,termA)) break;}
97  i=0;
98  while (strcmp(unE,"GeV"))
99  {
100  int res;
101  //printf("%i\n",i);
102  res = fscanf(F2,"%lf %s %s %lf %lf %s %s %s %s %s\n",&e,unE,dee,
103  &den,&r,unR,b,un3,cc,un4);
104  if(!strcmp(unE,"keV")) {Energy=e/1000.0;}
105  else if(!strcmp(unE,"MeV")) {Energy=e;}
106  else if(!strcmp(unE,"GeV")) {Energy=e*1000.0;}
107  else {printf("Error of reading in ReadRint a: %s %s\n",unE,unR);}
108 
109  if(!strcmp(unR,"A")) {Rng=r*1.0E-8;}
110  else if(!strcmp(unR,"um")) {Rng=r/10000.0;}
111  else if(!strcmp(unR,"mm")) {Rng=r/10.0;}
112  else if(!strcmp(unR,"cm")) {Rng=r;}
113  else if(!strcmp(unR,"m")) {Rng=r*100.0;}
114  else if(!strcmp(unR,"km")) {Rng=r*100000.0;}
115  else {printf("Error of reading in ReadRint b: %s %s\n",unE,unR);}
116 
117  Ranges[0][i]=Energy;
118  Ranges[1][i]=Rng;
119  i++;
120  }
121  fclose(F2);
122  }
123  for (i=0;i<4;++i)
124  {
125  x[i]=Ranges[0][i];
126  y[i]=Ranges[1][i];
127  }
128  if(intrp4(x,y,c))
129  {
130  printf("1 Interpolation error 1 intrp4=%i\n",intrp4(x,y,c));
131  return 1;
132  }
133 
134  Ranges[2][0]=c[0];
135  Ranges[3][0]=c[1];
136  Ranges[4][0]=c[2];
137  Ranges[5][0]=c[3];
138 
139  for (i=0;i<4;++i)
140  {
141  x[i]=Ranges[1][i];
142  y[i]=Ranges[0][i];
143  }
144 
145  if(intrp4(x,y,c))
146  {
147  printf("2 Interpolation error 2 intrp4=%i\n",intrp4(x,y,c));return 1;
148  }
149 
150  Ranges[6][0]= c[0];
151  Ranges[7][0]= c[1];
152  Ranges[8][0]= c[2];
153  Ranges[9][0]= c[3];
154 
155  for(i=1;i<103;++i)
156  {
157  x[0] = Ranges[0][i-1];
158  x[1] = Ranges[0][i];
159  x[2] = Ranges[0][i+1];
160  x[3] = Ranges[0][i+2];
161  y[0] = Ranges[1][i-1];
162  y[1] = Ranges[1][i];
163  y[2] = Ranges[1][i+1];
164  y[3] = Ranges[1][i+2];
165 
166  if(intrp4(x,y,c))
167  {printf("3 Interpolation error 3 intrp4=%i\n",intrp4(x,y,c));return 1;}
168  Ranges[2][i] = c[0];
169  Ranges[3][i] = c[1];
170  Ranges[4][i] = c[2];
171  Ranges[5][i] = c[3];
172  }
173  Ranges[2][i] = c[0];
174  Ranges[3][i] = c[1];
175  Ranges[4][i] = c[2];
176  Ranges[5][i] = c[3];
177 
178  Ranges[2][i+1] = c[0];
179  Ranges[3][i+1] = c[1];
180  Ranges[4][i+1] = c[2];
181  Ranges[5][i+1] = c[3];
182 
183  for(i=1;i<103;++i)
184  {
185  y[0] = Ranges[0][i-1];
186  y[1] = Ranges[0][i];
187  y[2] = Ranges[0][i+1];
188  y[3] = Ranges[0][i+2];
189 
190  x[0] = Ranges[1][i-1];
191  x[1] = Ranges[1][i];
192  x[2] = Ranges[1][i+1];
193  x[3] = Ranges[1][i+2];
194 
195  if(intrp4(x,y,c))
196  {
197  printf("4 Interpolation error 4 intrp4=%i\n",intrp4(x,y,c));return 1;
198  }
199  Ranges[6][i] = c[0];
200  Ranges[7][i] = c[1];
201  Ranges[8][i] = c[2];
202  Ranges[9][i] = c[3];
203  }
204 
205  Ranges[6][i] = c[0];
206  Ranges[7][i] = c[1];
207  Ranges[8][i] = c[2];
208  Ranges[9][i] = c[3];
209 
210  Ranges[6][i+1] = c[0];
211  Ranges[7][i+1] = c[1];
212  Ranges[8][i+1] = c[2];
213  Ranges[9][i+1] = c[3];
214 
215  return 0;
216 }
217 
218 //-----------------------------------------------------------------------------
219 int intrp4(double* x,double* y, double* c)
220 {
221  //______________________________________________________________________
222  // returns c0,c1,c2,c3 coeff. of y= c0 + c1*x + c2*x^2 + c3*x^3 function
223  // passing through 4 points:
224  // x1,y1; x2,y2; x3,y3; x4,y4
225  //______________________________________________________________________|
226  double d0,d1,d2,d3;
227  double x12,x13,x22,x23,x32,x33,x42,x43;
228  int rp4;
229 
230  x12 = x[0]*x[0];
231  x13 = x12*x[0];
232  x22 = x[1]*x[1];
233  x23 = x22*x[1];
234  x32 = x[2]*x[2];
235  x33 = x32*x[2];
236  x42 = x[3]*x[3];
237  x43 = x42*x[3];
238 
239  d3=(x13*x22*x[2]-x12*x23*x[2]-x13*x[1]*x32+x[0]*x23*x32+x12*x[1]*x33-
240  x[0]*x22*x33-x13*x22*x[3]+x12*x23*x[3]+x13*x32*x[3]-x23*x32*x[3]-
241  x12*x33*x[3]+x22*x33*x[3]+x13*x[1]*x42-x[0]*x23*x42-x13*x[2]*x42+
242  x23*x[2]*x42+x[0]*x33*x42-x[1]*x33*x42-x12*x[1]*x43+x[0]*x22*x43+
243  x12*x[2]*x43-x22*x[2]*x43-x[0]*x32*x43+x[1]*x32*x43);
244 
245  c[3]=(x22*x[2]*y[0]-x[1]*x32*y[0]-x22*x[3]*y[0]+
246  x32*x[3]*y[0]+x[1]*x42*y[0]-x[2]*x42*y[0]-
247  x12*x[2]*y[1]+x[0]*x32*y[1]+x12*x[3]*y[1]-
248  x32*x[3]*y[1]-x[0]*x42*y[1]+x[2]*x42*y[1]+
249  x12*x[1]*y[2]-x[0]*x22*y[2]-x12*x[3]*y[2]+
250  x22*x[3]*y[2]+x[0]*x42*y[2]-x[1]*x42*y[2]-
251  x12*x[1]*y[3]+x[0]*x22*y[3]+x12*x[2]*y[3]-
252  x22*x[2]*y[3]-x[0]*x32*y[3]+x[1]*x32*y[3])/d3;
253 
254  d2=(x13*x22*x[2]-x12*x23*x[2]-x13*x[1]*x32+x[0]*x23*x32+x12*x[1]*x33-
255  x[0]*x22*x33-x13*x22*x[3]+x12*x23*x[3]+x13*x32*x[3]-x23*x32*x[3]-
256  x12*x33*x[3]+x22*x33*x[3]+x13*x[1]*x42-x[0]*x23*x42-x13*x[2]*x42+
257  x23*x[2]*x42+x[0]*x33*x42-x[1]*x33*x42-x12*x[1]*x43+x[0]*x22*x43+
258  x12*x[2]*x43-x22*x[2]*x43-x[0]*x32*x43+x[1]*x32*x43);
259 
260  c[2]=(-(x23*x[2]*y[0])+x[1]*x33*y[0]+x23*x[3]*y[0]-x33*x[3]*y[0]-x[1]*x43*y[0]+
261  x[2]*x43*y[0]+x13*x[2]*y[1]-x[0]*x33*y[1]-
262  x13*x[3]*y[1]+x33*x[3]*y[1]+x[0]*x43*y[1]-
263  x[2]*x43*y[1]-x13*x[1]*y[2]+x[0]*x23*y[2]+
264  x13*x[3]*y[2]-x23*x[3]*y[2]-x[0]*x43*y[2]+
265  x[1]*x43*y[2]+x13*x[1]*y[3]-x[0]*x23*y[3]-
266  x13*x[2]*y[3]+x23*x[2]*y[3]+x[0]*x33*y[3]-
267  x[1]*x33*y[3])/d2;
268 
269  d1=(x13*x22*x[2]-x12*x23*x[2]-x13*x[1]*x32+x[0]*x23*x32+x12*x[1]*x33-
270  x[0]*x22*x33-x13*x22*x[3]+x12*x23*x[3]+x13*x32*x[3]-x23*x32*x[3]-
271  x12*x33*x[3]+x22*x33*x[3]+x13*x[1]*x42-x[0]*x23*x42-x13*x[2]*x42+
272  x23*x[2]*x42+x[0]*x33*x42-x[1]*x33*x42-x12*x[1]*x43+x[0]*x22*x43+
273  x12*x[2]*x43-x22*x[2]*x43-x[0]*x32*x43+x[1]*x32*x43);
274 
275  c[1]=(x23*x32*y[0]-x22*x33*y[0]-x23*x42*y[0]+x33*x42*y[0]+x22*x43*y[0]-
276  x32*x43*y[0]-x13*x32*y[1]+x12*x33*y[1]+x13*x42*y[1]-x33*x42*y[1]-
277  x12*x43*y[1]+x32*x43*y[1]+x13*x22*y[2]-x12*x23*y[2]-x13*x42*y[2]+
278  x23*x42*y[2]+x12*x43*y[2]-x22*x43*y[2]-x13*x22*y[3]+x12*x23*y[3]+
279  x13*x32*y[3]-x23*x32*y[3]-x12*x33*y[3]+x22*x33*y[3])/d1;
280 
281  d0=(x13*x22*x[2]-x12*x23*x[2]-x13*x[1]*x32+x[0]*x23*x32+x12*x[1]*x33-
282  x[0]*x22*x33-x13*x22*x[3]+x12*x23*x[3]+x13*x32*x[3]-x23*x32*x[3]-
283  x12*x33*x[3]+x22*x33*x[3]+x13*x[1]*x42-x[0]*x23*x42-x13*x[2]*x42+
284  x23*x[2]*x42+x[0]*x33*x42-x[1]*x33*x42-x12*x[1]*x43+x[0]*x22*x43+
285  x12*x[2]*x43-x22*x[2]*x43-x[0]*x32*x43+x[1]*x32*x43);
286 
287  c[0]=(-(x23*x32*x[3]*y[0])+x22*x33*x[3]*y[0]+x23*x[2]*x42*y[0]-x[1]*x33*x42*y[0]-
288  x22*x[2]*x43*y[0]+x[1]*x32*x43*y[0]+x13*x32*x[3]*y[1]-x12*x33*x[3]*y[1]-
289  x13*x[2]*x42*y[1]+x[0]*x33*x42*y[1]+x12*x[2]*x43*y[1]-x[0]*x32*x43*y[1]-
290  x13*x22*x[3]*y[2]+x12*x23*x[3]*y[2]+x13*x[1]*x42*y[2]-x[0]*x23*x42*y[2]-
291  x12*x[1]*x43*y[2]+x[0]*x22*x43*y[2]+x13*x22*x[2]*y[3]-x12*x23*x[2]*y[3]-
292  x13*x[1]*x32*y[3]+x[0]*x23*x32*y[3]+x12*x[1]*x33*y[3]-x[0]*x22*x33*y[3])/d0;
293 
294  if(d0*d1*d2*d3!=0.0) {rp4=0;}
295  else {rp4=-1;}
296 
297  return rp4;
298 }
299 
300 
301 //Double_t CalcElossIntegralVolStep(Double_t T, const G4ParticleDefinition& ion,
302 // const G4Material& mat, const Double_t range,
303 // const Double_t intStep/*=1e-4 [1mk]*/) {
304 // if (mat.GetName() == "vacuum") {
305 // return 0.;
306 // }
307 // if (range <= 0.)
308 // return 0;
309 // Double_t integralEloss = 0.;
310 // Double_t curStep = 0.;
311 // const Double_t tempStep = range / 10.;
312 // G4EmCalculator* calc = new G4EmCalculator();
313 // while (curStep < range) {
314 // Double_t eloss = calc->ComputeDEDX(T, &ion, "hIoni", &mat) * tempStep * 10 /*cm to mm*/;
315 // integralEloss += eloss;
316 // T += eloss;
317 // curStep += tempStep;
318 // }
319 // return integralEloss;
320 //}
321 
322 
323 Double_t CalcElossIntegralVolStep (Double_t kineticEnergy, const G4ParticleDefinition& particle,
324  const G4Material& material, const Double_t range, const Double_t /*intStep=1e-4*/) {
325  if (range <= 0.)
326  return 0;
327  Double_t integralEloss = 0.;
328  const Double_t intStep = range / 20;
329  Double_t curStep = 0.;
330  G4EmCalculator* calc = new G4EmCalculator();
331  while (curStep < range) {
332  Double_t eloss = calc->GetDEDX(kineticEnergy, &particle, &material) * intStep * 10 ; // MeV
333  integralEloss += eloss;
334  kineticEnergy += eloss;
335  curStep += intStep;
336  }
337  return integralEloss;
338 }
339 
340 Double_t ComputeElossIntegralVolStep(Double_t kineticEnergy, const G4ParticleDefinition& particle,
341  const G4Material& material, const TString& process_name,
342  const Double_t range) {
343  if (range <= 0.)
344  return 0;
345  Double_t integralEloss = 0.;
346  const Double_t intStep = range / 20;
347  Double_t curStep = 0.;
348  G4EmCalculator* calc = new G4EmCalculator();
349  while (curStep < range) {
350  Double_t eloss = calc->ComputeDEDX(kineticEnergy, &particle,
351  process_name.Data(), &material) * intStep * 10 ; // MeV
352  integralEloss += eloss;
353  kineticEnergy += eloss;
354  curStep += intStep;
355  }
356  return integralEloss;
357 }
358 
359 Double_t ElossCalculator::CalcDeDx(const TString& material_name, const int pdg, const double kineticEnergy) {
360  const auto* particle = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
361  const auto* material = G4NistManager::Instance()->FindOrBuildMaterial(material_name.Data());
362  G4EmCalculator calc;
363  return calc.GetDEDX(kineticEnergy, particle, material);
364 }
365 
366 Double_t ElossCalculator::CalcDeDx_long(const TString& material_name, const int pdg, double kineticEnergy, const Double_t range) {
367  const auto* particle = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
368  const auto* material = G4NistManager::Instance()->FindOrBuildMaterial(material_name.Data());
369  G4EmCalculator calc;
370 
371  Double_t integralEloss = 0.;
372  const Double_t intStep = range / 20;
373  Double_t curStep = 0.;
374 
375  while (curStep < range) {
376  Double_t eloss = calc.GetDEDX(kineticEnergy, particle, material) * intStep * 10 ; // MeV
377  integralEloss += eloss;
378  kineticEnergy += eloss;
379  curStep += intStep;
380  }
381  return integralEloss;
382  //return calc.GetDEDX(kineticEnergy, particle, material);
383 }
384 
385 ClassImp(ElossCalculator)
Generated on Tue Dec 7 2021 13:58:45 for er by   doxygen 1.8.11