FairRoot/PandaRoot
PndFtsSingleStraw.cxx
Go to the documentation of this file.
1 //
2 // ============================================================================
3 //
4 // Straw tube response simulation
5 // Corrections to StrawCharge and Eject for delta electrons (dec 2007)
6 //
7 // version with one threshold only Julich signal
8 //
9 // mandatory call at the beginning
10 // once to set constants
11 //
12 // void TConst(Double_t Radius, Double_t pSTP, Double_t ArP, Double_t CO2P);
13 //
14 // Pysical (Garfield) values:
15 // Radius Press (pSTP) % Ar (ArP) % CO2 (CO2P)
16 // 0.4 1 0.9 0.1
17 // 0.5 1 0.9 0.1
18 // 0.5 2 0.8 0.2
19 //
20 // ------------------------------------------------------------------------------
21 //
22 // Full simulation: calls at each track for MC applications
23 //
24 // define particle and its in-out hits
25 //
26 // void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
27 // Double_t w4, Double_t w5, Double_t w6); define the wire
28 //
29 // Double_t PartToTime(Double_t Mass, Double_t Momentum,
30 // Double_t InOut[] ); straw time
31 //
32 //
33 // Double_t TimnsToDiscm(Double_t time); from time (ns)
34 // to the reconstructed distance (cm)
35 //
36 // Double_t PartToADC(); ADC charge
37 // signal corresponding to the energy loss of StrawCharge
38 // for energy loss simulation
39 //
40 //
41 //-----------------------------------------------------------------------------
42 //
43 // Fast simulation: calls at each track for MC applications
44 //
45 // define particle and its in-out hits
46 //
47 // void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
48 // Double_t w4, Double_t w5, Double_t w6); define the wire
49 //
50 // Double_t FastRec(Double_t TrueDcm, Int_t Flag); fast approximation:
51 // from the true distance (cm)
52 // to the reconstructed one (cm)
53 // Flag=0 standard
54 // Flag=1 Julich expt. curve
55 // (active with press =2 in TCo
56 //
57 // void TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]); constant initialization
58 //
59 // Double_t FastPartToADC(); ADC charge
60 // signal corresponding to the energy loss of StrawCharge
61 // with the Urban model for energy loss simulation
62 //
63 // ----------------------------------------------------------------------------
64 // where:
65 //
66 // time: measured drift time in ns
67 //
68 // pSTP = pressure in bar
69 // Radius = tube radius
70 // ArP and CO2P percentages of Argon and CO2
71 // w1:w6 = coordinates of the extremes of the straw wire (straw axis)
72 // Mass, Momentm = mass and momentum of the MC particle
73 // InOut[0:5] = straw input and output coordinates traversed from the MC particle
74 //
75 //
76 //
77 // For other applications see the routines listed in PndFtsSingleStraw.h
78 //
79 // Author: A. Rotondi (november 2005, revised and extended in August 2006)
80 //
81 // ================================================================================
82 //
83 
84 #include <iostream>
85 #include <string>
86 #include "TRandom.h"
87 #include "TMath.h"
88 #include "PndFtsSingleStraw.h"
89 #include "TImage.h"
90 #include "TVector3.h"
91 #include "TVectorD.h"
92 
93 #include "TMatrixD.h"
94 #include "TMatrixDEigen.h"
95 
96 
97 using namespace std;
98 
99  ClassImp(PndFtsSingleStraw);
100 
101 
102 // ============================================================
103 PndFtsSingleStraw::PndFtsSingleStraw()
104  :CDist(), CDistC(), CNele(), CNeleT(), TeleTime(), AmplSig(),Pulse(), PulseT(), WDist(),
105  Wi(0), ArPerc(0), CO2Perc(0), CH4Perc(0), ArWPerc(0), CO2WPerc(0), CH4WPerc(0), pSTP(0), Radius(0),
106  AAr(0), ZAr(0), RhoAr(0), NclAr(0), EmedAr(0), EminAr(0), EmpAr(0), CsiAr(0), IAr(0), WiAr(0), Ncl(0),
107  Ecl(0), Lcl(0), Ntote(0), GasGain(0), Cutoff(0),
108  EmedCO2(0), EminCO2(0), EmpCO2(0), CsiCO2(0), ACO2(0), ZCO2(0), RhoCO2(0), ICO2(0), WiCO2(0), NclCO2(0),
109  EmedCH4(0), EmpCH4(0), CsiCH4(0), EminCH4(0), ACH4(0), ZCH4(0), RhoCH4(0), ICH4(0), WiCH4(0), NclCH4(0),
110  RhoMixCO2(0), RhoMixCH4(0), PZeta(0), piMass(0), PMass(0), PMom(0), Dx(0), eMass(0), prMass(0), Delta(0),
111  CNumb(0), PEn(0), beta(0), gamma(0), Emed(0), Emin(0), Csi(0), Emax(0), Emp(0), NNClus(0),
112  Xin(0), Yin(0), Zin(0), Xout(0), Yout(0), Zout(0), Rpath(0), NPolya(0), Xmax(0), bPolya(0),
113  Calpha(0), Cbeta(0), Cgamma(0),
114  NUrban(0), SigUrb(0), Eup(0), AvUrb(0), Wx1(0), Wy1(0), Wz1(0), Wx2(0), Wy2(0), Wz2(0),
115  Wp(0), Wq(0), Wr(0), PulseMax(0), PulseTime(0), Thresh1(0), Thresh2(0), Nchann(0), Out1(0), Out2(0), Out3(0)
116 {
117  // class constructor
118 
119  // clear
120  memset(CumClus,0,sizeof(CumClus));
121  memset(CH4Clus,0,sizeof(CH4Clus));
122  memset(PolyaCum,0,sizeof(PolyaCum));
123  memset(Xs,0,sizeof(Xs));
124 
125  CDist.clear();
126  CDistC.clear();
127  CNele.clear();
128  CNeleT.clear();
129  TeleTime.clear();
130  AmplSig.clear();
131  Pulse.clear();
132  PulseT.clear();
133  WDist.clear();
134 }
135 
136 // ----------------------------------------------------------------------
137 void PndFtsSingleStraw::TConst(Double_t R1, Double_t P1, Double_t A1, Double_t C1) {
138 
139  // set constants for the simulation
140  // Input
141  // P1 = tube absolute pressure pSTP
142  // R1 = tube radius Radius
143  // A1 = argon percentage ArPerc
144  // C1 = C=2 percentage C=2Perc
145  //
146 
147  Radius=R1;
148  pSTP=P1;
149  // Input for the media (volume percentages)
150  ArPerc = A1;
151  CO2Perc = C1;
152 
153  // cluster dimensions in Ar and CO2 (experimantal values)
154 
155 // Double_t PClus[20] =
156 // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
157 // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
158 // .00095, .00075, .00063}; // Fischle fit
159 
160  Double_t PClus[20] =
161  {0., .802, .0707, .020, .013, .008, .006, .005, .006,
162  .008, .009, .007, .0050, .0040, .0033, .0029, .0025,
163  .0023, .0022, .002}; // Lapique 1st calculation
164 
165 // Double_t PClus[20] =
166 // {0., .841, .0340, .021, .013, .008, .006, .004, .003,
167 // .008, .013, .008, .0050, .004, .0030, .0028, .0025,
168 // .0023, .0022, .002}; // Lapique 2nd calculation
169 
170 
171 // Double_t PClus[20] =
172 // {0., .656, .150, .064, .035, .0225, .0155, .0105, .0081,
173 // .0061, .0049, .0039, .0030, .0025, .0020, .0016, .0012,
174 // .00080, .00059, .00045}; // Fischle empirical
175 
176 // PDouble_t Clus[20] =
177 // {0., .656, .148, .0649, .0337, .0244, .0141, .0078, .0095,
178 // .0063, .0062, .0042, .0028, .0018, .0023, .0017, .0014,
179 // .00060, .00050, .00063}; // Fischle exp
180 
181  Double_t CO2Clus[20] =
182  {0., .730, .162, .038, .020, .0110, .0147, .0060, .0084,
183  .0052, .0020, .0042, .0021, .0025, .0038, .0021, .0009,
184  .00013, .00064, .00048}; // Fischle exp
185 
186 // Double_t CH4Clus[20] =
187 // {0., .786, .120, .032, .013, .0098, .0055, .0057, .0027,
188 // .0029, .0020, .0016, .0013, .0010, .0012, .0006, .0005,
189 // .00042, .00037, .00033}; // Fischle exp
190 
191 
192  CH4Perc = 0.07;
193 
194  // -----------------------------------------------------
195  // gain of the avalanche
196  // Ar/CO2 90/10 1 bar (NTP) MAGY GARFIELD 20-7-2006
197  GasGain=100000.;
198 
199  // argon ----------------------------------------------------
200  AAr = 39.948; // Argon (39.948)
201  ZAr= 18.0; // Argon (18)
202  RhoAr = pSTP*1.78*1.e-03; // g/cm3 (1.78 mg/cm3)
203  IAr = 188*1.e-09; // ionization potential (GeV) (188 eV)
204  WiAr =27.0; // energy to create an ion pair (standard 26.7 eV)
205  NclAr = 25.; // cluster/cm in Argon
206  // CO2 -----------------------------------------------------
207  ACO2 = 44; // CO2
208  ZCO2 = 22.; // CO2
209  RhoCO2 = pSTP*1.98*1.e-03; // g/cm3 CO2 (1.98 mg/cm3)
210  ICO2 = 95.8*1.e-09; // ionization potential (GeV) (96 eV)
211  WiCO2 = 33.0; // energy to create an ion pair (33 eV)
212  NclCO2 = 35.5; // clusters/cm CO2 35.5
213  // Methane CH4 ---------------------------------------------------------
214  ACH4 = 16; // CO2 (39.948)
215  ZCH4 = 10.; // CO2 (18)
216  RhoCH4 = pSTP*0.71*1.e-03; // g/cm3 CO2 (0.71 mg/cm3)
217  ICH4 = 40.6*1.e-09; // ionization potential (GeV) (45 eV)
218  WiCH4 = 28.0; // energy to create an ion pair
219  NclCH4 = 25.0;
220  // Input for the media (weight percentages) ----------------------------
221  ArWPerc = ArPerc *AAr /(ArPerc*AAr + CO2Perc*ACO2);
222  CO2WPerc = CO2Perc*ACO2/(ArPerc*AAr + CO2Perc*ACO2);
223 
224  // mixture densiies ----------------------------------------------------
225  RhoMixCO2 = 1./((ArWPerc/RhoAr) + (CO2WPerc/RhoCO2));
226  RhoMixCH4 = 1./((ArWPerc/RhoAr) + (CH4WPerc/RhoCH4));
227 
228  //----------------------------------------------------------------------
229  // particles (Gev, energy losses in Kev)
230 
231  PZeta = 1; // projectile charge
232  piMass = 0.139; // particle mass (GeV)
233  eMass = 0.511/1000.; // electron mass (GeV) (0.511 MeV)
234  prMass = 0.93827; // proton mass (GeV)
235 
236  // ---------------------------------------------------------------------
237  // thresholds for the straw tubes (default values) see TInit for current values
238  Thresh1=10;
239  Thresh2=30;
240  // channels for the signal
241  Nchann = 500;
242 
243  // ---------------------------------------------Emin------------------------
244 
245  NPolya= 100; // steps for the calculation of the Polya distributions
246  Xmax=0.; // Polya istribution is calculated between o and Xmax (see Polya)
247  bPolya = 0.5;
248  Polya(bPolya); // cumulative of the Polya distribution
249  // -----------------------------------------------------------------------
250 
251  // cumulative for the number of electron per cluster
252 
253  Double_t Wnorm = (ArPerc*NclAr + CO2Perc*NclCO2);
254  CumClus[0]=(ArPerc*NclAr*PClus[0] + CO2Perc*NclCO2*CO2Clus[0])/Wnorm;
255 
256  for(Int_t i=1; i<=20; i++) {
257  CumClus[i]=(ArPerc*NclAr*PClus[i] + CO2Perc*NclCO2*CO2Clus[i])/Wnorm
258  + CumClus[i-1];
259  }
260  CumClus[20]=1.;
261 
262  Double_t sum=0.;
263  for(Int_t i=0; i<=19; i++) {
264  sum += PClus[i];
265  // cout<<" PClus["<<i<<"] = "<<PClus[i]<<endl;
266  }
267  for(Int_t i=0;i<=20;i++) { //cout<<"CumClus["<<i<<"] = "<<CumClus[i]<<endl;
268 }
269  // cout<<" Sum of Probabilities = "<<sum<<endl;
270 
271 }
272 
273 // ==============================================================
274 void PndFtsSingleStraw::PutWireXYZ(Double_t w1, Double_t w2, Double_t w3,
275  Double_t w4, Double_t w5, Double_t w6){
276  // get wire coordinates
277 
278  Wx1=w1;
279  Wy1=w2;
280  Wz1=w3;
281 
282  Wx2=w4;
283  Wy2=w5;
284  Wz2=w6;
285 
286 }
287 // =============================================================
288 
289 
290 void PndFtsSingleStraw::TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
291 
292  // initialization of the constants for each track and each straw
293 
294  // transfer of data
295  // track geometrical quantities
296  Xin=InOut[0]; Yin=InOut[1]; Zin=InOut[2];
297  Xout=InOut[3]; Yout=InOut[4]; Zout=InOut[5];
298  PMass = xPMass;
299  PMom = xPMom;
300 
301  // path into the straw
302  Dx = sqrt((Xout-Xin)*(Xout-Xin) +
303  (Yout-Yin)*(Yout-Yin) +
304  (Zout-Zin)*(Zout-Zin));
305 
306  // clear
307  CNumb=0;
308  PulseMax=0;
309  PulseTime=0;
310  CDist.clear();
311  CDistC.clear();
312  CNele.clear();
313  CNeleT.clear();
314  WDist.clear();
315  TeleTime.clear();
316  AmplSig.clear();
317  Pulse.clear();
318  PulseT.clear();
319 
320  // ---------------------------------------------------------------------
321 
322  // initialization of the constants
323  // maximum energy transfer Emax (GeV) and related quantities
324 
325  PEn = sqrt(PMom*PMom + PMass*PMass); // particle energy GeV
326  beta = PMom/PEn;
327  gamma= 1/sqrt(1.-beta*beta);
328  Double_t begam = beta*gamma; // local
329  Double_t mratio= eMass/PMass; // local
330 
331  // calculation of the polarization Sternheimer factor for Argon
332  Double_t ms = 2.80;
333  Double_t Xst = log10(begam*sqrt(pSTP));
334  Double_t Cs =-11.92;
335  Double_t X1 = 4;
336  Double_t Xo = 1.96;
337  Double_t as = 0.389;
338  Delta=0;
339  if(Xo <= Xst && Xst <= X1) Delta = 4.6051702*Xst +Cs +as*pow(X1-Xst,ms);
340  else if(X1< Xst) Delta = 4.6051702*Xst +Cs;
341  if(Delta < 0) Delta=0;
342 
343  // calculation of other typical quantites
344 
345  Emax = (1.022*begam*begam/(1.+2.*gamma*mratio+mratio*mratio))/1000.; // GeV
346 
347  CsiAr = 0.5*0.3071*PZeta*PZeta*ZAr*RhoAr/(beta*beta*AAr)*1e-03; //GeV
348  Double_t fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(IAr*IAr);
349  EmedAr = 2.*CsiAr*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
350  EminAr = pSTP*2.70*1.e-06; // GeV/cm (2.5 keV)
351  // most prob energy (GeV)
352  EmpAr = EmedAr + CsiAr*(0.422784 + beta*beta + log(CsiAr/Emax));
353 
354  CsiCO2 = 0.5*0.3071*PZeta*PZeta*ZCO2*RhoCO2/(beta*beta*ACO2)*1e-03; //GeV
355  fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(ICO2*ICO2);
356  EmedCO2 = 2.*CsiCO2*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
357  EminCO2 = pSTP*3.60*1.e-06; // GeV/cm (2.5 keV)
358  // most prob energy (GeV)
359  EmpCO2 = EmedCO2 + CsiCO2*(0.422784 + beta*beta + log(CsiCO2/Emax));
360 
361  Csi = RhoMixCO2 * ((ArWPerc*CsiAr/RhoAr) + (CO2WPerc*CsiCO2/RhoCO2));;
362  Emed = RhoMixCO2 * ((ArWPerc*EmedAr/RhoAr) + (CO2WPerc*EmedCO2/RhoCO2));
363  Emin = RhoMixCO2 * ((ArWPerc*EminAr/RhoAr) + (CO2WPerc*EminCO2/RhoCO2));
364  Emp = RhoMixCO2 * ((ArWPerc*EmpAr/RhoAr) + (CO2WPerc*EmpCO2/RhoCO2));
365 
366  // mean weighted interaction
367  Wi = (ArPerc*NclAr*WiAr + CO2Perc*NclCO2*WiCO2)
368  /(ArPerc*NclAr + CO2Perc*NclCO2);
369  // number of clusters
370  Double_t nAr = ArWPerc * RhoMixCO2/AAr;
371  Double_t nCO2 = CO2WPerc * RhoMixCO2/ACO2;
372  //Ncl=pSTP*((ArPerc*NclAr)+(CO2Perc*NclCO2))*Emed/Emin; // <--- Cluster/cm
373  Ncl=pSTP*((nAr*NclAr+nCO2*NclCO2)/(nAr+nCO2))*Emed/Emin; // <--- Cluster/cm
374 
375  Lcl=1./Ncl; // mean free path between clusters (cm)
376  Ecl = 2.8; // mean number of electrons per clusters (def 2.8)
377  Cutoff = Ncl*Ecl*25; // limit to the number of primary electrons (*20)
378  Ntote = Ecl*Ncl; // total mean number of electrons per cm
379 
380  // ---------------------------------------------------------------------
381  // thresholds for the straw tubes (current values)
382  // total number of electrons * scaling factor (max of signal x electron=1)
383 
384  Thresh1=Ncl*Ecl* 0.10;
385  Thresh2=Ncl*Ecl* 0.15;
386 
387  // -----------------------------------------------------------------------
388 
389  TDirCos(); // director cosine of the track
390 
391  // control prints
392 
393 // cout<<" Dx "<<Dx<<" Csi "<<Csi*1.e+09<<" Emax MeV "<<
394 // Emax*1.e+03<<" PMass "<<PMass<<" gamma "<<
395 // gamma<<endl;
396 // cout<<" RRise "<<RRise(gamma)<<endl;
397 // cout<<" Thresh1 = "<<Thresh1<<" Thresh2 = "<<Thresh2<<endl;
398 
399 }
400 
401 
402 // =============================================================
403 
404 // energy loss in GeV
405 Double_t PndFtsSingleStraw::STEloss() {
406 
407 
408  return gRandom->Landau(Emp*Dx,Csi*Dx); // in GeV
409 
410 }
411 // =============================================================
412 
413 
414 Double_t PndFtsSingleStraw::STUrban() {
415 
416 //
417 // energy loss in GeV from the Urbam distribution NIM A 362(1995)416
418 // Author: A. Rotondi December 2007
419 //
420 
421  Double_t Sig1Ar, Sig2Ar, Sig3Ar;
422  Double_t Sig1Mix, Sig2Mix, Sig3Mix;
423  Double_t Sig1CO2, Sig2CO2, Sig3CO2, IMix, nAr, nCO2;
424  Double_t f1, f2, e1, e2, e2f, ru, Cu, Cuc;
425  Double_t N1Mix, N2Mix, N3Mix, E1Mix, E2Mix, E3Mix, ZMix;
426  Double_t Sig21, Sig22, Sig23, Medd;
427 
428  Double_t eVGeV = 1.e-09;
429 
430  // ru= 0.6;
431  ru= 0.55;
432 
433  Cuc = 2.*PMass*beta*beta*gamma*gamma;
434 
435  // calculation for Argon
436 
437  Cu = EmedAr;
438  f2 = 2./ZAr;
439  f1 = 1. -f2;
440  e2 = 10.*ZAr*ZAr*eVGeV;
441  e2f = TMath::Power(e2,f2);
442  e1 = TMath::Power(IAr/e2f,1./f1);
443  Sig1Ar = (1.-ru)*Cu*(f1/e1)*(TMath::Log(Cuc/e1) - beta*beta)/
444  (TMath::Log(Cuc/IAr) - beta*beta);
445  Sig2Ar = (1.-ru)*Cu*(f2/e2)*(TMath::Log(Cuc/e2) - beta*beta)/
446  (TMath::Log(Cuc/IAr) - beta*beta);
447  Sig3Ar = Cu*Emax*ru / (IAr*(Emax+IAr)*TMath::Log((Emax+IAr)/IAr));
448 
449  // Calculation for CO2
450 
451  Cu = EmedCO2;
452  f2 = 2./ZCO2;
453  f1 = 1. -f2;
454  e2 = 10.*ZCO2*ZCO2*eVGeV;
455  e2f = TMath::Power(e2,f2);
456  e1 = TMath::Power(ICO2/e2f,1./f1);
457 
458  Sig1CO2 = (1.-ru)*Cu*(f1/e1)*(TMath::Log(Cuc/e1) - beta*beta)/
459  (TMath::Log(Cuc/ICO2) - beta*beta);
460  Sig2CO2 = (1.-ru)*Cu*(f2/e2)*(TMath::Log(Cuc/e2) - beta*beta)/
461  (TMath::Log(Cuc/ICO2) - beta*beta);
462  Sig3CO2 = Cu*Emax*ru / (ICO2*(Emax+ICO2)*TMath::Log((Emax+ICO2)/ICO2));
463 
464  // calculation for the mixture
465 
466 
467  nAr = ArWPerc * RhoMixCO2/AAr;
468  nCO2 = CO2WPerc * RhoMixCO2/ACO2;
469  IMix = TMath::Exp((nAr*ZAr*TMath::Log(IAr) + nCO2*ZCO2*TMath::Log(ICO2))/
470  (nAr*ZAr + nCO2*ZCO2));
471 
472  ZMix = ArWPerc*ZAr + CO2WPerc*ZCO2;
473  f2 = 2./ZMix;
474  f1 = 1. -f2;
475  E2Mix = 10.*ZMix*ZMix*eVGeV;
476  e2f = TMath::Power(E2Mix,f2);
477  E1Mix = TMath::Power(IMix/e2f,1./f1);
478 
479  Sig1Mix = Dx*RhoMixCO2*(ArWPerc*Sig1Ar/RhoAr + CO2WPerc*Sig1CO2/RhoCO2);
480  Sig2Mix = Dx*RhoMixCO2*(ArWPerc*Sig2Ar/RhoAr + CO2WPerc*Sig2CO2/RhoCO2);
481  Sig3Mix = Dx*RhoMixCO2*(ArWPerc*Sig3Ar/RhoAr + CO2WPerc*Sig3CO2/RhoCO2);
482 
483  NUrban = Sig1Mix+Sig2Mix+Sig3Mix; // total number of collisions (mean value)
484 
485  // Urban distribution calculation
486 
487  N1Mix = gRandom->Poisson(Sig1Mix);
488  N2Mix = gRandom->Poisson(Sig2Mix);
489  N3Mix = gRandom->Poisson(Sig3Mix);
490 
491  Int_t N3 = TMath::Nint(N3Mix+gRandom->Uniform());
492 
493  E3Mix = 0.;
494  for(Int_t j=1;j<=N3;j++){
495  E3Mix += ((Emax+IMix)*IMix)/(IMix + Emax*(1.-gRandom->Uniform()));
496  }
497 
498  // variance calculation
499 
500  Eup = IMix/(1.-0.98*Emax/(Emax+IMix)); // 98% of the delta electron area (GeV)
501  Sig21 = IMix*(Emax+IMix)/Emax;
502  Medd = Sig21*TMath::Log(Eup/IMix);
503  Sig22 = Sig21*Eup - Medd*Medd;
504  Sig23 = Sig1Mix*E1Mix*E1Mix + Sig2Mix*E2Mix*E2Mix
505  + Sig3Mix*Medd*Medd + Sig3Mix*Sig22;
506 
507  AvUrb = Sig1Mix*E1Mix + Sig2Mix*E2Mix + Sig3Mix*Medd;
508 
509  SigUrb = TMath::Sqrt(Sig23);
510 
511  // return the energy lost in keV
512 
513 
514 
515  return N1Mix*E1Mix + N2Mix*E2Mix + E3Mix; // in GeV
516 
517 
518 }
519 // =============================================================
520 
521 
522 Double_t PndFtsSingleStraw::StrawCharge() {
523 
524  // sampling of the ionization energy loss, number of
525  // primary electrons and their wire distances for each track
526  // energy loss in GeV
527 
528  // clear
529  Double_t TotEnEle = 0., ClusEner=0., Ecurr;
530  CNeleT.clear();
531  WDist.clear();
532 
533  // threshold for delta rays (>1.5 KeV from NIM A301(1991)202)
534  // total energy released by the electrons
535  //Number of cluster in the length
536  NNClus = Cluster(); // set vectors CDist CDistC e CNele
537  // calculation of total energy released by the electrons
538  for(Int_t j=1;j<=NNClus;j++){
539  ClusEner = 0.;
540  for(Int_t jj=1; jj<=(Int_t)(CNele.at(j-1)); jj++){
541  Ecurr = Wi;
542  TotEnEle += Ecurr;
543  ClusEner += Ecurr;
544 
545  }
546 
547  // effective number of electrons per cluster
548 
549  CNeleT.push_back(int(ClusEner/Wi));
550  }
551 
552  // record the distance from the wire of any electron
553  // adding the diffusion effects
554  Int_t owfl =0;
555  for(Int_t k=1; k<=NNClus; k++) {
556  for (Int_t kk=0; kk<(Int_t)CNeleT.at(k-1); kk++) {
557  if(owfl > (Int_t) Cutoff) break;
558  owfl++;
559  WDist.push_back(WDistCalc(CDistC.at(k-1)));
560  }
561  }
562 
563  return TotEnEle*1.e-09; // in Gev
564 }
565 
566 // =============================================================
567 
568 
569 Int_t PndFtsSingleStraw::Cluster() {
570 
571  // calculate the number of clusters CNumb
572  // their distance CDist and
573  // their number of electrons CNele
574 
575  CNumb=0;
576  CDist.clear();
577  CDistC.clear();
578  CNele.clear();
579 
580  Double_t DisTot = 0;
581 
582  for(Int_t k=1; k<=1000; k++) {
583  // distance
584 
585  Double_t path = -Lcl*log(1-gRandom->Uniform()+0.000001);
586  DisTot+= path;
587  if(DisTot>Dx) break;
588  CNumb=k;
589  CDist.push_back(path);
590  CDistC.push_back(DisTot);
591  CNele.push_back((Double_t)Eject() );
592  }
593 
594  return CNumb;
595 }
596 
597 
598 // =============================================================
599 
600 
601 Int_t PndFtsSingleStraw::Eject() {
602  // find the number of electrons in a cluster
603  Int_t Inelect;
604  Double_t nelect=0.;
605 
606  nelect= TMath::BinarySearch(21,CumClus,gRandom->Uniform());
607  if(nelect<19) {
608  return Inelect = (Int_t) (nelect) +1;
609  }
610  else {
611  // long delta electron tail
612  return Inelect = (Int_t)(20./(1.03-gRandom->Uniform()));
613  }
614 }
615 
616 
617 // =============================================================
618 
619 
620 Double_t PndFtsSingleStraw::RRise(Double_t gamma2) {
621 
622  // interpolate the relativisic rise of the
623  // number of cluster per cm, starting from the one
624  // measured at the ionization minimum
625 
626  Double_t Rise;
627  Double_t lg = log10(gamma2);
628 
629  if(1.<=gamma2 && gamma2 <= 2.2){
630  Rise = -2.159*lg +1.7;
631  }
632  else if(2.2<=gamma2 && gamma2 <= 6.){
633  Rise = 1.;
634  }
635  else if(6.<=gamma2 && gamma2 <= 200.){
636  Rise = 0.302*lg + 0.765;
637  }
638  else if(200.<=gamma2 && gamma2 <= 1000.){
639  Rise = 0.1431*lg + 1.131;
640  }
641  else if(1000.<=gamma2){
642  Rise = 1.54;
643  }
644  else{
645  Rise= 1.7;
646  }
647 
648  return Rise;
649 }
650 
651 // ==========================================================================
652 
653 void PndFtsSingleStraw::Polya(Double_t bpar) {
654 
655  // calculate the cumulative of the Polya distribution
656  // for samplig the gain fluctuations
657 
658  Double_t eps = 0.0001;
659  Double_t Dxx = 0.05, xx=0., x1,x2;
660  Double_t PMax;
661 
662  Double_t k=1./bpar -1.;
663 
664  // find Xmax
665 
666  PMax = eps * pow(k,k)*exp(-k);
667  Double_t value = 1.e+06;
668  Xmax =2*k;
669  while(value > PMax){
670  Xmax +=Dxx;
671  value = pow(Xmax,k)*exp(-Xmax);
672  }
673  Xmax += -0.5*Dxx;
674 
675 
676  // calculate the cumulative
677 
678  Double_t dx= Xmax/NPolya;
679  Xs[0]=0;
680  for(Int_t i=1; i<NPolya; i++){
681  x1 = xx;
682  xx += dx;
683  x2= xx;
684  Xs[i]=x2;
685  if(i>1)
686  PolyaCum[i] = PolyaCum[i-1] +
687  0.5*dx*(pow(x1,k)*exp(-x1) + pow(x2,k)*exp(-x2))/TMath::Gamma(k+1);
688  else
689  PolyaCum[i] = PolyaCum[i-1] +
690  0.5*dx* pow(x2,k)*exp(-x2)/TMath::Gamma(k+1);
691  }
692 
693  // adjust the normalization
694 
695  for(Int_t ii=0; ii<NPolya; ii++) PolyaCum[ii] /= PolyaCum[NPolya-1];
696 
697 }
698 
699 
700 // =====================================================================
701 
702 Double_t PndFtsSingleStraw::PolyaSamp() {
703 
704  // sampling a wire gain fluctuation in the gas
705 
706  Double_t xr = gRandom->Uniform();
707  Int_t n = TMath::BinarySearch(NPolya,PolyaCum,xr);
708  Double_t xsamp = Xmax;
709  if(n<NPolya-1){
710  xsamp= Xs[n] +
711  ((xr-PolyaCum[n])/(PolyaCum[n+1] - PolyaCum[n])) *(Xs[n+1]-Xs[n]);
712  }
713  return xsamp;
714 }
715 
716 // =====================================================================
717 
718 TVector3 PndFtsSingleStraw::WDistCalc(Double_t DisTot) {
719 
720  // calculation of the distance of the cluster from the wire
721  // DisTot is the cluster distance from the track entrance
722  // diffusion effects are considered
723 
724  // wire director cosines
725 
726  Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
727  (Wy2-Wy1)*(Wy2-Wy1) +
728  (Wz2-Wz1)*(Wz2-Wz1) );
729  Wp = (Wx2-Wx1)/Wlength;
730  Wq = (Wy2-Wy1)/Wlength;
731  Wr = (Wz2-Wz1)/Wlength;
732 
733  // current track coordinates
734  Double_t xcor = Calpha * DisTot + Xin;
735  Double_t ycor = Cbeta * DisTot + Yin;
736  Double_t zcor = Cgamma * DisTot + Zin;
737 
738  // cross product VV1= (wire x track) for the distance
739  Double_t XX1 = Wr*(ycor-Wy1) - Wq*(zcor-Wz1);
740  Double_t YY1 = Wp*(zcor-Wz1) - Wr*(xcor-Wx1);
741  Double_t ZZ1 = Wq*(xcor-Wx1) - Wp*(ycor-Wy1);
742 
743  // vector for the 3-D distance from the wire (from wire x VV1)
744  Double_t XX = Wq*ZZ1 - Wr*YY1;
745  Double_t YY = Wr*XX1 - Wp*ZZ1;
746  Double_t ZZ = Wp*YY1 - Wq*XX1;
747  //cout<<" XYZ "<<XX<<" "<<YY<<" "<<ZZ<<endl;
748  Double_t DDistcm = sqrt(XX*XX + YY*YY +ZZ*ZZ);
749 
750  //director cosines of the distance vector
751  Double_t cosx = XX/DDistcm;
752  Double_t cosy = YY/DDistcm;
753  Double_t cosz = ZZ/DDistcm;
754 
755  // sampling of the diffusion
756  //Double_t DDist = DDistcm*10.; // distance in mm
757 
758  // longitudinal coefficient of diffusion (GARFIELD) cm --> micron
759  // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
760  Double_t SigL = DiffLong(DDistcm);
761  // ... in cm
762  SigL *= 1.e-04; // in cm
763 
764  // tranverse coefficient (GARFIELD) cm--> micron
765  // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
766  Double_t SigT = DiffTran(DDistcm);
767  SigT *= 1.e-04; // in cm
768 
769  // sampling of Longitudinal and Transverse diffusion
770  Double_t difL = gRandom->Gaus(0.,SigL);
771  Double_t difT = gRandom->Gaus(0.,SigT);
772 
773  // vector addition to the distance
774  // the transverse component has the same dir cos of the wire
775  XX += difL*cosx + difT*Wp;
776  YY += difL*cosy + difT*Wq;
777  ZZ += difL*cosz + difT*Wr;
778  //cout<<" XYZ+ dif "<<XX<<" "<<YY<<" "<<ZZ<<endl;
779  //cout<<" --------------------------------------------------- "<<endl;
780  TVector3 TDist(XX,YY,ZZ);
781 
782  return TDist;
783 
784 }
785 
786 
787 // =====================================================================
788 
789 void PndFtsSingleStraw::TDirCos() {
790 
791  // director cosines of the track (called for each track)
792 
793  //path into the straw
794  Rpath = sqrt((Xout-Xin)*(Xout-Xin) +
795  (Yout-Yin)*(Yout-Yin) +
796  (Zout-Zin)*(Zout-Zin));
797 
798  //director cosines
799  Calpha = (Xout-Xin)/Rpath;
800  Cbeta = (Yout-Yin)/Rpath;
801  Cgamma = (Zout-Zin)/Rpath;
802 }
803 
804 
805 // =====================================================================
806 
807 Int_t PndFtsSingleStraw::TimeEle() {
808 
809  TeleTime.clear();
810 
811  // calculate the arrival times (ns) for each electron in TeleTime
812  // from cm to ns <---------------------------------
813  // return the number of electrons
814  // Author:A. Rotondi 20-7-2006
815 
816  Int_t ido = WDist.size();
817  Double_t etime;
818 
819  // cutoff the total charge to 20 times the average
820  if(ido > (Int_t) Cutoff) ido = (Int_t) Cutoff;
821 
822  for(Int_t k=0; k<ido; k++) {
823  Double_t dst = WDist.at(k).Mag(); // distance in cm
824 
825  // space to time calculation from GARFIELD
826 
827  if(pSTP <1.9){
828  if(Radius<0.5){
829  // V=1600V diameter 4 mm 90/10
830  etime= -1.624e-05 + 0.1258*dst + 0.8079*pow(dst,2) -2.918*pow(dst,3)
831  + 10.33*pow(dst,4) -10.84*pow(dst,5);
832  }
833  else{
834  // V=1600V diameter 5 mm 90/10
835  etime= -6.763e-05 + 0.1471*dst +0.3625*pow(dst,2) +0.3876*pow(dst,3)
836  +1.04*pow(dst,4) -1.693*pow(dst,5);
837  }
838  }
839  else{
840  // 2 bar 2000V, 5 mm, 80/20
841  etime= -0.0001014 + 0.1463*dst -0.1694*pow(dst,2) +2.4248*pow(dst,3)
842  -1.793*pow(dst,4);
843  }
844 
845  etime*=1000.; // nano seconds
846  TeleTime.push_back(etime);
847  }
848 
849  return TeleTime.size();
850 }
851 
852 // =====================================================================
853 
854 Double_t PndFtsSingleStraw::Signal(Double_t t, Double_t t0) {
855 
856  // electric signal at time t of a cluster arriving at t0
857  Double_t elesig;
858  //Double_t A = 1.03e-03; //[R.K. 01/2017] unused variable?
859  //Double_t B = 3.95; //[R.K. 01/2017] unused variable?
860  //Double_t C = 0.228; //[R.K. 01/2017] unused variable?
861  //Double_t D = -3.839e-02; //[R.K. 01/2017] unused variable?
862  //Double_t E = 1.148e-03; //[R.K. 01/2017] unused variable?
863 
864  Double_t x= t - t0;
865  //if(x>0) elesig = A*exp(B*log(x)-C*x)*(1+D*x+E*x*x); // Sokolov Signal
866  if(x>0) elesig =pow(2.*x/10.,2)*exp(-2.*x/10.); // Wirtz signal
867  else elesig = 0.;
868 
869  return elesig;
870 
871 }
872 
873 // =====================================================================
874 
875 Int_t PndFtsSingleStraw::StrawSignal(Int_t nsteps) {
876 
877  // creation of nstep values of
878  // the straw global Pulse (sum on all clusters)
879  // return the number of primary electrons
880 
881  PulseMax=0;
882  Pulse.clear();
883  PulseT.clear();
884  AmplSig.clear();
885 
886  Int_t neltot = TimeEle(); // creation and size of TeleTime (electron times)
887 
888  Double_t Tmax = 1.e-25;
889  for(Int_t k=0; k< neltot; k++) if(Tmax<TeleTime.at(k)) Tmax=TeleTime.at(k);
890  Tmax += 100.;
891 
892  Double_t Dt = Tmax/nsteps; // number of steps of the signal
893 
894  //AmplSig is the amplitude of each electron
895 
896  for(Int_t j=0; j< neltot; j++){
897  AmplSig.push_back(bPolya*PolyaSamp());
898  }
899 
900  // creation of the signal Pulse(PulseT) time in ns
901  for(Int_t j=0; j< nsteps; j++){
902  Double_t te = j*Dt;
903  Double_t sumele=0.;
904  for(Int_t jj=0; jj< neltot; jj++){
905  Double_t te0 = TeleTime.at(jj);
906  sumele += AmplSig.at(jj)*Signal(te,te0);
907  }
908 
909  Pulse.push_back(sumele);
910  PulseT.push_back(te);
911  }
912 
913  // add a random noise (3% of maximum) plus
914  // a 1% of 200 ns periodic signal
915 
916  Double_t Pmax = 1.e-25;
917  for(Int_t k=0; k<(Int_t)Pulse.size(); k++)
918  if(Pmax<Pulse.at(k)) Pmax=Pulse.at(k);
919 
920  for(Int_t k=0; k<(Int_t)Pulse.size(); k++){
921  Pulse.at(k) += 0.03*Pmax*gRandom->Uniform()
922  + 0.01*Pmax*(sin(6.28*PulseT.at(k)/120.));
923  // if(Pulse[k]<0) Pulse[k] *= -1.;
924  }
925 
926  PulseMax=Pmax;
927 
928  // set variable threshold for the signals (constant fraction)
929  // Thresh1=0.05*Pmax;
930  // Thresh2=0.15*Pmax;
931 
932  return neltot;
933 }
934 
935 // =====================================================================
936 
937 Int_t PndFtsSingleStraw::StrawTime() {
938 
939  // simulate the discrimination of the straw signal
940  // and give the time
941  // discriminator technique: set 2 threshold, select the first one
942  // (the low one) only if the second one is fired (FINUDA system)
943 
944 
945  PulseTime=0.;
946 
947  Int_t ind=0;
948  Int_t flag1=0;
949  Int_t flag2=1;
950 
951  for(Int_t k=0; k<(Int_t) Pulse.size(); k++){
952  if(flag1==0 && Pulse[k]>Thresh1) {
953  flag1=1;
954  ind=k;
955  }
956 
957 // if(Pulse[k]<Thresh1) {flag1=0; ind=0;} // reset if signal decreases
958 
959 // if(flag1==1 && Pulse[k]>Thresh2) {
960 // flag2=1;
961 // break;
962 // }
963 
964  }
965  if(flag1==1 && flag2==1) PulseTime=PulseT.at(ind);
966 
967  return ind;
968 }
969 
970 
971 
972 
973 // =====================================================================
974 
975 Double_t PndFtsSingleStraw::TrueDist(Double_t Point[]) {
976  // service routine that finds the distance in cm from the wire
977  // by knowing the wire coordinates (class variables)
978  // and the input-output points Point[6]
979 
980 
981  Double_t truedist = 0;
982 
983  // wire director cosines
984  Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) +
985  (Wy2-Wy1)*(Wy2-Wy1) +
986  (Wz2-Wz1)*(Wz2-Wz1) );
987  Wp = (Wx2-Wx1)/Wlength;
988  Wq = (Wy2-Wy1)/Wlength;
989  Wr = (Wz2-Wz1)/Wlength;
990 
991  // director cosines of the given track
992  Double_t Modu = sqrt( (Point[3]*Point[0])*(Point[3]*Point[0]) +
993  (Point[4]*Point[1])*(Point[4]*Point[1]) +
994  (Point[5]*Point[2])*(Point[5]*Point[2]) );
995  Double_t dcx = (Point[3]-Point[0])/Modu;
996  Double_t dcy = (Point[4]-Point[1])/Modu;
997  Double_t dcz = (Point[5]-Point[2])/Modu;
998 
999  //distance formula
1000  Double_t p1 = (Point[0]-Wx1)*(dcy*Wr-dcz*Wq);
1001  Double_t p2 =-(Point[1]-Wy1)*(dcx*Wr-dcz*Wp);
1002  Double_t p3 = (Point[2]-Wz1)*(dcx*Wq-dcy*Wp);
1003  Double_t Det = p1+p2+p3;
1004  Double_t Disc = sqrt( (dcy*Wr-dcz*Wq)*(dcy*Wr-dcz*Wq) +
1005  (dcz*Wp-dcx*Wr)*(dcz*Wp-dcx*Wr) +
1006  (dcx*Wq-dcy*Wp)*(dcx*Wq-dcy*Wp) );
1007  if(Disc >0) truedist = TMath::Abs(Det/Disc);
1008 
1009  return truedist; // distance in cm
1010 }
1011 
1012 // =====================================================================
1013 
1014 Double_t PndFtsSingleStraw::TimnsToDiscm(Double_t time) {
1015 
1016  // distance in cm from time in ns for pSTP =1 and pSTP=2 atm
1017  // utility routine for the track reconstruction
1018  //last update: A. Rotondi 3-3-2007
1019 
1020  Double_t drift;
1021 
1022 
1023  // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
1024  // ns --> mm
1025 
1026  if(pSTP < 1.9) {
1027  if(Radius < 0.5){
1028 
1029  drift = -0.0140 -1.37281e-01
1030  +5.13978e-02*time
1031  +7.65443e-05*pow(time,2)
1032  -9.53479e-06*pow(time,3)
1033  +1.19432e-07*pow(time,4)
1034  -6.19861e-10*pow(time,5)
1035  +1.35458e-12*pow(time,6)
1036  -1.10933e-15*pow(time,7);
1037  }
1038 
1039 
1040  else{
1041  // 1600 V 5 mm
1042  if(time < 120.){
1043 
1044  drift = 0.0300 -1.07377e-01
1045  +3.65134e-02*time
1046  +1.20909e-03*pow(time,2)
1047  -4.56678e-05*pow(time,3)
1048  +6.70207e-07*pow(time,4)
1049  -4.99204e-09*pow(time,5)
1050  +2.19079e-11*pow(time,6)
1051  -8.01791e-14*pow(time,7)
1052  +2.16778e-16*pow(time,8);
1053  }
1054  else{
1055 
1056  drift = 0.0300 -8.91701e-01
1057  +4.68487e-02*time
1058  +1.00902e-03*pow(time,2)
1059  -4.00359e-05*pow(time,3)
1060  +6.23768e-07*pow(time,4)
1061  -5.20556e-09*pow(time,5)
1062  +2.41502e-11*pow(time,6)
1063  -5.85450e-14*pow(time,7)
1064  +5.77250e-17*pow(time,8);
1065  }
1066  }
1067  }
1068 
1069  else {
1070  // 2 absolute atm 5 mm 80/20 (1 bar overpressure)
1071  if(time <= 50.) {
1072 
1073  // on thresh static 10%
1074 
1075  drift = -0.0300 +1.28551e-02
1076  +1.44029e-02*time
1077  -3.67834e-03*pow(time,2)
1078  +3.32034e-04*pow(time,3)
1079  -6.36592e-06*pow(time,4)
1080  -7.82907e-08*pow(time,5)
1081  +3.58931e-09*pow(time,6)
1082  -2.93491e-11*pow(time,7) ;
1083 
1084 
1085  // one thresh static 3%
1086 // drift = +5.35238e-02
1087 // -4.25959e-02 *time
1088 // +7.59448e-03 *pow(time,2)
1089 // -1.44009e-04 *pow(time,3)
1090 // -1.76365e-06 *pow(time,4)
1091 // +3.29531e-08 *pow(time,5)
1092 // +1.12115e-09 *pow(time,6)
1093 // -1.72919e-11 *pow(time,7) ;
1094 
1095 
1096 
1097  }
1098  else if(50. < time && time < 130.) {
1099 
1100 
1101  // on thresh static 10%
1102 
1103  drift = -0.0190 +4.40993e-01
1104  -2.91442e-02*time
1105  +3.06237e-03*pow(time,2)
1106  -6.07870e-05*pow(time,3)
1107  +5.97431e-07*pow(time,4)
1108  -3.09238e-09*pow(time,5)
1109  +7.70537e-12*pow(time,6)
1110  -6.49086e-15*pow(time,7) ;
1111 
1112 
1113 
1114  // one thresh static 3%
1115 // drift = +2.25702e-02
1116 // +5.17806e-02 *time
1117 // +2.53060e-04 *pow(time,2)
1118 // -1.60338e-05 *pow(time,3)
1119 // +2.14805e-07 *pow(time,4)
1120 // -1.30249e-09 *pow(time,5)
1121 // +3.38791e-12 *pow(time,6)
1122 // -2.10503e-15 *pow(time,7) ;
1123 
1124 
1125  }
1126 
1127  else{
1128 
1129  // on thresh static 10%
1130  drift = -0.0100 +4.28757e-01
1131  -1.95413e-02*time
1132  +3.02333e-03*pow(time,2)
1133  -6.13920e-05*pow(time,3)
1134  +5.93656e-07*pow(time,4)
1135  -3.05271e-09*pow(time,5)
1136  +8.05446e-12*pow(time,6)
1137  -8.59626e-15*pow(time,7) ;
1138 
1139  // one thresh static 3%
1140 // drift = +1.57217e-01
1141 // +5.79365e-02*time
1142 // +2.15636e-04*pow(time,2)
1143 // -1.66405e-05*pow(time,3)
1144 // +2.13726e-07*pow(time,4)
1145 // -1.26888e-09 *pow(time,5)
1146 // +3.68533e-12 *pow(time,6)
1147 // -4.24032e-15 *pow(time,7) ;
1148 
1149 
1150  }
1151  }
1152 
1153 
1154  drift = 0.1*drift;
1155  if(drift < 0.) drift=0.;
1156  return drift;
1157 
1158 }
1159 
1160 // --------------------------------------------------------------------------------
1161 
1162 Double_t PndFtsSingleStraw::PartToTime(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {
1163 
1164  // find the time of a particle of mass xPmass, momentum xPMom, with
1165  // input-output coordinate InOut[6]
1166  // Useful for MC pplication after a call to PutWireXYZ
1167 
1168  TInit(xPMass, xPMom, InOut); // start the event
1169  Out1 = StrawCharge(); // energy loss (GeV) to generate charge
1170  Out2 = StrawSignal(Nchann); // generate the straw signal
1171  Out3 = StrawTime(); // find the straw drift time PulseTime
1172 
1173  return PulseTime;
1174 
1175 }
1176 
1177 // --------------------------------------------------------------------------------
1178 Double_t PndFtsSingleStraw::PartToADC() {
1179 
1180  // return the energy loss in the tube as charge signal
1181  // taking into account the Polya fluctuations
1182  Double_t ADCsignal=0.;
1183 
1184  for(Int_t j=1; j<= NNClus; j++){
1185 
1186  for(Int_t jc=1; jc<=(Int_t) CNeleT[j-1]; jc++)
1187  ADCsignal += bPolya * GasGain * PolyaSamp();
1188  }
1189 
1190  return ADCsignal;
1191 }
1192 
1193 // --------------------------------------------------------------------------------
1194 
1195 Double_t PndFtsSingleStraw::FastPartToADC() {
1196 
1197  // return the energy loss (from the Urban distribution)
1198  // in the tube as charge signal
1199  // taking into account the Polya fluctuations
1200 
1201  Double_t ADCsignal=0.;
1202 
1203  // number of elecrons. Wi is the mean energy lost per free
1204  // electrn in eV
1205  Int_t NtotEle = (Int_t) (1.e+09 * STUrban()/Wi);
1206 
1207  for(Int_t j=1; j<= NtotEle; j++){
1208 
1209  ADCsignal += bPolya * GasGain * PolyaSamp();
1210 
1211  }
1212 
1213  return ADCsignal;
1214 }
1215 
1216 
1217 // ----------------------------------------------------------------------------------
1218 Double_t PndFtsSingleStraw::FastRec(Double_t TrueDcm, Int_t Flag) {
1219 
1220  // having the true distance TrueDcm (cm) as input,
1221  // return the reconstructed distance (cm), in a fast
1222  // and approximated way
1223  // by sampling on the simulated reconstruction curve
1224  // When Press =2 and Flag=1 one uses the julich experimental data
1225  // A. Rotondi March 2007
1226 
1227  Double_t resmic;
1228 
1229  // 1 atm pressure
1230  if(pSTP < 1.9) {
1231  if(Radius < 0.45){
1232  if(TrueDcm < 0.38){
1233  resmic = 1.24506e+02 -1.80117e+02*TrueDcm
1234  +3.76905e+03*pow(TrueDcm,2) -4.63251e+04*pow(TrueDcm,3)
1235  +1.80068e+05*pow(TrueDcm,4) -2.21094e+05*pow(TrueDcm,5);
1236  }
1237  else resmic = 57.;
1238  }
1239  // radius > 0.4 cm
1240  else{
1241  if(TrueDcm < 0.48){
1242  resmic = 1.53656e+02
1243  -5.07749e+03*TrueDcm
1244  +1.73707e+05*pow(TrueDcm,2)
1245  -2.72285e+06*pow(TrueDcm,3)
1246  +2.28719e+07*pow(TrueDcm,4)
1247  -1.12921e+08*pow(TrueDcm,5)
1248  +3.39427e+08*pow(TrueDcm,6)
1249  -6.12741e+08 *pow(TrueDcm,7)
1250  +6.12041e+08 *pow(TrueDcm,8)
1251  -2.60444e+08*pow(TrueDcm,9) ;
1252  }
1253  else resmic = 72.;
1254 
1255  }
1256  }
1257  // 2 atm pressure radius 5 cm
1258 
1259  else {
1260  if(Flag==0) {
1261  if(TrueDcm < 0.48){
1262  resmic = +1.06966e+02
1263  -4.03073e+03 *TrueDcm
1264  +1.60851e+05 *pow(TrueDcm,2)
1265  -2.87722e+06 *pow(TrueDcm,3)
1266  +2.67581e+07 *pow(TrueDcm,4)
1267  -1.43397e+08 *pow(TrueDcm,5)
1268  +4.61046e+08 *pow(TrueDcm,6)
1269  -8.79170e+08 *pow(TrueDcm,7)
1270  +9.17095e+08 *pow(TrueDcm,8)
1271  -4.03253e+08 *pow(TrueDcm,9) ;
1272  }
1273  else resmic=30.;
1274  }
1275  else{
1276  // data from julich
1277  if(TrueDcm < 0.48){
1278  resmic = 20. +1.48048e+02
1279  -3.35951e+02*TrueDcm
1280  -1.87575e+03*pow(TrueDcm,2)
1281  +1.92910e+04*pow(TrueDcm,3)
1282  -6.90036e+04*pow(TrueDcm,4)
1283  +1.07960e+05*pow(TrueDcm,5)
1284  -5.90064e+04*pow(TrueDcm,6) ;
1285  }
1286  else resmic=65.;
1287  }
1288  }
1289 
1290  //real distance in cm
1291  Double_t rsim = gRandom->Gaus(TrueDcm, resmic*0.0001);
1292  if (rsim<0.) rsim = TrueDcm - TrueDcm*gRandom->Uniform(0.,1.);
1293  else if (rsim>0.5) rsim = TrueDcm + (0.5-TrueDcm)*gRandom->Uniform(0.,1.);
1294 
1295  return rsim;
1296 }
1297 // --------------------------------------------------------------------------------
1298 Double_t PndFtsSingleStraw::DiffLong(Double_t Distcm) {
1299 
1300  // return the longitudinal diffusion
1301  // from cm to microns
1302  // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
1303  //
1304  Double_t DiffMic;
1305 
1306  if(pSTP<1.9){
1307  if(Radius < 0.5){
1308  // 4 mm 1600 V 90/10
1309  DiffMic = 0.896 + 1387.*Distcm - 1.888e+04*pow(Distcm,2)
1310  + 1.799e+05*pow(Distcm,3) - 9.848e+05*pow(Distcm,4)
1311  + 3.009e+06*pow(Distcm,5) - 4.777e+06*pow(Distcm,6)
1312  + 3.074e+06*pow(Distcm,7);
1313  }
1314  else{
1315  // 5 mm 1600 V
1316  DiffMic = 1.537 + 1246.*Distcm - 1.357e+04*pow(Distcm,2)
1317  + 1.049e+05*pow(Distcm,3) - 4.755e+05*pow(Distcm,4)
1318  + 1.211e+06*pow(Distcm,5) - 1.6e+06*pow(Distcm,6)
1319  + 8.533e+05*pow(Distcm,7);
1320  }
1321  }
1322 
1323  else{
1324  // 2000V 5 mm 2 bar 80/20
1325  DiffMic = 2.135 +818.*Distcm - 1.044e+04*pow(Distcm,2)
1326  + 8.31e+04*pow(Distcm,3) - 3.492e+05*pow(Distcm,4)
1327  + 7.959e+05*pow(Distcm,5) - 9.378e+05*pow(Distcm,6)
1328  + 4.492e+05*pow(Distcm,7);
1329  }
1330  return DiffMic;
1331 }
1332 // --------------------------------------------------------------------------------
1333 Double_t PndFtsSingleStraw::DiffTran(Double_t Distcm) {
1334 
1335  // return the transverse diffusion in microns
1336  // from cm to microns
1337  // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
1338  //
1339  Double_t DiffMic;
1340 
1341  if(pSTP<1.9){
1342  if(Radius < 0.5){
1343  // 4 mm 1600 V 90/10
1344 // DiffMic = + 0.8513 + 1648.*Distcm - 1.085e+04*pow(Distcm,2)
1345 // + 7.38e+04*pow(Distcm,3) - 3.025e+05*pow(Distcm,4)
1346 // + 6.067e+05*pow(Distcm,5) - 4.643e+04*pow(Distcm,6);
1347  DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
1348  + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
1349  + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
1350  }
1351  else{
1352  // 5 mm 1600 V 90/10
1353  DiffMic = + 1.482 + 1529.*Distcm - 6755.*pow(Distcm,2)
1354  + 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
1355  + 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
1356  }
1357  }
1358  else{
1359 
1360  // 5 mm 2000 V 2 bar 80/20
1361  DiffMic = +2.094 + 1138.*Distcm - 7557.*pow(Distcm,2)
1362  + 2.968e+04*pow(Distcm,3) - 6.577e+04*pow(Distcm,4)
1363  + 7.581e+04*pow(Distcm,5) - 3.497e+04*pow(Distcm,6);
1364  }
1365  return DiffMic;
1366 }
1367 // --------------------------------------------------------
1368 Double_t PndFtsSingleStraw::DistEle(Double_t time) {
1369 
1370  // dist in cm from time in ns for pSTP =1 and pSTP=2
1371  // utility routine for the >SINGLE< electron reconstruction
1372  //last update: A. Rotondi 20-7-2006
1373 
1374  Double_t drift;
1375 
1376  // 1 absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
1377  // ns --> cm
1378 
1379  time *= 0.001; // time in micro s
1380 
1381  if(pSTP < 2.) {
1382  if(Radius <0.5){
1383  // 1600 V 4 mm 90/10
1384  drift = 0.001629 + 6.194*time
1385  - 56.55*pow(time,2) + 355.8*pow(time,3)
1386  - 903.2*pow(time,4);
1387  }
1388  else{
1389  // 1600 V 5 mm 90/10
1390  drift = 0.003365 + 5.734*time
1391  - 41.88*pow(time,2) + 191.2*pow(time,3)
1392  - 333.4*pow(time,4) ;
1393  }
1394  }
1395 
1396  else {
1397 
1398  // 2 absolute atm 2000 V 5 mm 80/20 (1 bar overpressure)
1399 
1400  drift = 0.003365 + 7.056*time
1401  - 62.28*pow(time,2) + 306.1*pow(time,3)
1402  - 558.7*pow(time,4);
1403 
1404  }
1405  return drift; // in cm
1406 
1407 }
1408 
PndFtsSingleStraw::TeleTime
std::vector< Double_t > TeleTime
Definition: PndFtsSingleStraw.h:152
PndFtsSingleStraw::DiffTran
Double_t DiffTran(Double_t distcm)
Definition: PndFtsSingleStraw.cxx:1333
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: P4_F32vec4.h:109
i
Int_t i
Definition: run_full.C:25
f1
TF1 * f1
Definition: reco_analys2.C:50
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndFtsSingleStraw::AmplSig
std::vector< Double_t > AmplSig
Definition: PndFtsSingleStraw.h:152
PndFtsSingleStraw::CDist
std::vector< Double_t > CDist
Definition: PndFtsSingleStraw.h:152
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndFtsSingleStraw::Pulse
std::vector< Double_t > Pulse
Definition: PndFtsSingleStraw.h:153
PndFtsSingleStraw::WDist
std::vector< TVector3 > WDist
Definition: PndFtsSingleStraw.h:154
n
int n
Definition: HypStatDecay.cxx:21
log
friend F32vec4 log(const F32vec4 &a)
Definition: P4_F32vec4.h:110
PndFtsSingleStraw::TInit
void TInit(Double_t Mass, Double_t Momentum, Double_t InOut[])
Definition: PndFtsSingleStraw.cxx:290
lg
TGeoVolume * lg
Definition: createROOTEnDrcGeo.C:126
PndFtsSingleStraw::PulseT
std::vector< Double_t > PulseT
Definition: PndFtsSingleStraw.h:153
jj
Int_t jj
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:434
nsteps
double nsteps
Definition: dedx_bands.C:15
PndFtsSingleStraw::TimnsToDiscm
Double_t TimnsToDiscm(Double_t time)
Definition: PndFtsSingleStraw.cxx:1014
PndFtsSingleStraw::RRise
Double_t RRise(Double_t gamma)
Definition: PndFtsSingleStraw.cxx:620
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndFtsSingleStraw::TConst
void TConst(Double_t Radius, Double_t pSTP, Double_t ArP, Double_t CO2P)
Definition: PndFtsSingleStraw.cxx:137
PndFtsSingleStraw::TrueDist
Double_t TrueDist(Double_t Point[])
Definition: PndFtsSingleStraw.cxx:975
t0
Int_t t0
Definition: hist-t7.C:106
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
eps
double eps(TVector3 v1, TVector3 v2)
Definition: softtrigger_kin5.C:163
Cu
#define Cu
Definition: createSTTwithPassive.C:100
PndFtsSingleStraw::Signal
Double_t Signal(Double_t t, Double_t t0)
Definition: PndFtsSingleStraw.cxx:854
PndFtsSingleStraw::DistEle
Double_t DistEle(Double_t tns)
Definition: PndFtsSingleStraw.cxx:1368
PndFtsSingleStraw::WDistCalc
TVector3 WDistCalc(Double_t d)
Definition: PndFtsSingleStraw.cxx:718
p2
TPad * p2
Definition: hist-t7.C:117
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
PndFtsSingleStraw::PartToTime
Double_t PartToTime(Double_t Mass, Double_t Momentum, Double_t InOut[])
Definition: PndFtsSingleStraw.cxx:1162
PndFtsSingleStraw::PutWireXYZ
void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3, Double_t w4, Double_t w5, Double_t w6)
Definition: PndFtsSingleStraw.cxx:274
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndFtsSingleStraw::CNele
std::vector< Double_t > CNele
Definition: PndFtsSingleStraw.h:152
CAMath::Log
static T Log(const T &x)
Definition: PndCAMath.h:40
PndFtsSingleStraw::StrawSignal
Int_t StrawSignal(Int_t nsteps)
Definition: PndFtsSingleStraw.cxx:875
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
PndFtsSingleStraw::FastRec
Double_t FastRec(Double_t TrueDcm, Int_t Flag)
Definition: PndFtsSingleStraw.cxx:1218
p1
TPad * p1
Definition: hist-t7.C:116
PndFtsSingleStraw::DiffLong
Double_t DiffLong(Double_t distcm)
Definition: PndFtsSingleStraw.cxx:1298
t
TTree * t
Definition: bump_analys.C:13
PndFtsSingleStraw::Polya
void Polya(Double_t bpar)
Definition: PndFtsSingleStraw.cxx:653
PndFtsSingleStraw::CDistC
std::vector< Double_t > CDistC
Definition: PndFtsSingleStraw.h:152
PndFtsSingleStraw::PolyaCum
Double_t PolyaCum[100]
Definition: PndFtsSingleStraw.h:246
ClassImp
ClassImp(PndFtsSingleStraw)
PndFtsSingleStraw::CNeleT
std::vector< Double_t > CNeleT
Definition: PndFtsSingleStraw.h:152
CAMath::Nint
int Nint(float x)
Definition: PndCAMath.h:117
Generated on Sun Aug 18 2019 03:01:27 for FairRoot/PandaRoot by   doxygen 1.8.8