StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AliTPCRF1D.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 
17 
22 
23 #include <RVersion.h>
24 #include <Riostream.h>
25 #include <TCanvas.h>
26 #include <TClass.h>
27 #include <TF2.h>
28 #include <TH1.h>
29 #include <TMath.h>
30 #include <TPad.h>
31 #include <TString.h>
32 #include <TStyle.h>
33 #include "TBuffer.h"
34 
35 #include "AliTPCRF1D.h"
36 
37 extern TStyle * gStyle;
38 
39 Int_t AliTPCRF1D::fgNRF=100;
40 Float_t AliTPCRF1D::fgRFDSTEP=0.01;
41 
42 static Double_t funGauss(Double_t *x, Double_t * par)
43 {
45 
46  return TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]));
47 }
48 
49 static Double_t funCosh(Double_t *x, Double_t * par)
50 {
52 
53  return 1/TMath::CosH(3.14159*x[0]/(2*par[0]));
54 }
55 
56 static Double_t funGati(Double_t *x, Double_t * par)
57 {
59 
60  Float_t k3=par[1];
61  Float_t k3R=TMath::Sqrt(k3);
62  Float_t k2=(TMath::Pi()/2)*(1-k3R/2.);
63  Float_t k1=k2*k3R/(4*TMath::ATan(k3R));
64  Float_t l=x[0]/par[0];
65  Float_t tan2=TMath::TanH(k2*l);
66  tan2*=tan2;
67  Float_t res = k1*(1-tan2)/(1+k3*tan2);
68  return res;
69 }
70 
73 
75 ClassImp(AliTPCRF1D)
77 
78 
79 AliTPCRF1D::AliTPCRF1D(Bool_t direct,Int_t np,Float_t step)
80  :TObject(),
81  fNRF(0),
82  fDSTEPM1(0.),
83  fcharge(0),
84  forigsigma(0.),
85  fpadWidth(3.5),
86  fkNorm(0.5),
87  fInteg(0.),
88  fGRF(0),
89  fSigma(0.),
90  fOffset(0.),
91  fDirect(kFALSE),
92  fPadDistance(0.)
93 {
94  //default constructor for response function object
95  fDirect=direct;
96  if (np!=0)fNRF = np;
97  else (fNRF=fgNRF);
98  fcharge = new Float_t[fNRF];
99  if (step>0) fDSTEPM1=1./step;
100  else fDSTEPM1 = 1./fgRFDSTEP;
101  for(Int_t i=0;i<5;i++) {
102  funParam[i]=0.;
103  fType[i]=0;
104  }
105 
106 }
107 
108 AliTPCRF1D::AliTPCRF1D(const AliTPCRF1D &prf)
109  :TObject(prf),
110  fNRF(prf.fNRF),
111  fDSTEPM1(prf.fDSTEPM1),
112  fcharge(0),
113  forigsigma(prf.forigsigma),
114  fpadWidth(prf.fpadWidth),
115  fkNorm(prf.fkNorm),
116  fInteg(prf.fInteg),
117  fGRF(new TF1(*(prf.fGRF))),
118  fSigma(prf.fSigma),
119  fOffset(prf.fOffset),
120  fDirect(prf.fDirect),
121  fPadDistance(prf.fPadDistance)
122 {
123  //
124  //
125  for(Int_t i=0;i<5;i++) {
126  funParam[i]=0.;
127  fType[i]=0;
128  }
129  fcharge = new Float_t[fNRF];
130  memcpy(fcharge,prf.fcharge, fNRF*sizeof(Float_t));
131 
132  //PH Change the name (add 0 to the end)
133  TString s(fGRF->GetName());
134  s+="0";
135  fGRF->SetName(s.Data());
136 }
137 
138 AliTPCRF1D & AliTPCRF1D::operator = (const AliTPCRF1D &prf)
139 {
140  if(this!=&prf) {
141  TObject::operator=(prf);
142  fNRF=prf.fNRF;
143  fDSTEPM1=prf.fDSTEPM1;
144  delete [] fcharge;
145  fcharge = new Float_t[fNRF];
146  memcpy(fcharge,prf.fcharge, fNRF*sizeof(Float_t));
147  forigsigma=prf.forigsigma;
148  fpadWidth=prf.fpadWidth;
149  fkNorm=prf.fkNorm;
150  fInteg=prf.fInteg;
151  delete fGRF;
152  fGRF=new TF1(*(prf.fGRF));
153  //PH Change the name (add 0 to the end)
154  TString s(fGRF->GetName());
155  s+="0";
156  fGRF->SetName(s.Data());
157  fSigma=prf.fSigma;
158  fOffset=prf.fOffset;
159  fDirect=prf.fDirect;
160  fPadDistance=prf.fPadDistance;
161  }
162  return *this;
163 }
164 
165 
166 
167 AliTPCRF1D::~AliTPCRF1D()
168 {
169  //
170  delete [] fcharge;
171  delete fGRF;
172 }
173 
174 Float_t AliTPCRF1D::GetRF(Float_t xin)
175 {
176  //function which return response
177  //for the charge in distance xin
178  //return linear aproximation of RF
179  Float_t x = (xin-fOffset)*fDSTEPM1+fNRF/2;
180  Int_t i1=Int_t(x);
181  if (x<0) i1-=1;
182  Float_t res=0;
183  if (i1+1<fNRF &&i1>0)
184  res = fcharge[i1]*(Float_t(i1+1)-x)+fcharge[i1+1]*(x-Float_t(i1));
185  return res;
186 }
187 
188 Float_t AliTPCRF1D::GetGRF(Float_t xin)
189 {
190  //function which returnoriginal charge distribution
191  //this function is just normalised for fKnorm
192  if (fGRF != 0 )
193  return fkNorm*fGRF->Eval(xin)/fInteg;
194  else
195  return 0.;
196 }
197 
198 
199 void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
200  Float_t kNorm, Float_t sigma)
201 {
202  //adjust parameters of the original charge distribution
203  //and pad size parameters
204  fpadWidth = padwidth;
205  fGRF = GRF;
206  fkNorm = kNorm;
207  if (sigma==0) sigma= fpadWidth/TMath::Sqrt(12.);
208  forigsigma=sigma;
209  fDSTEPM1 = 10/TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
210  //sprintf(fType,"User");
211  snprintf(fType,5,"User");
212  // Update();
213 }
214 
215 
216 void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
217  Float_t kNorm)
218 {
219  //
220  // set parameters for Gauss generic charge distribution
221  //
222  fpadWidth = padWidth;
223  fkNorm = kNorm;
224  if (fGRF !=0 ) fGRF->Delete();
225  fGRF = new TF1("funGauss",funGauss,-5,5,1);
226  funParam[0]=sigma;
227  forigsigma=sigma;
228  fGRF->SetParameters(funParam);
229  fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
230  //by default I set the step as one tenth of sigma
231  //sprintf(fType,"Gauss");
232  snprintf(fType,5,"Gauss");
233 }
234 
235 void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
236  Float_t kNorm)
237 {
238  //
239  // set parameters for Cosh generic charge distribution
240  //
241  fpadWidth = padWidth;
242  fkNorm = kNorm;
243  if (fGRF !=0 ) fGRF->Delete();
244  fGRF = new TF1("funCosh", funCosh, -5.,5.,2);
245  funParam[0]=sigma;
246  fGRF->SetParameters(funParam);
247  forigsigma=sigma;
248  fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
249  //by default I set the step as one tenth of sigma
250  //sprintf(fType,"Cosh");
251  snprintf(fType,5,"Cosh");
252 }
253 
254 void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
255  Float_t kNorm)
256 {
257  //
258  // set parameters for Gati generic charge distribution
259  //
260  fpadWidth = padWidth;
261  fkNorm = kNorm;
262  if (fGRF !=0 ) fGRF->Delete();
263  fGRF = new TF1("funGati", funGati, -5.,5.,2);
264  funParam[0]=padDistance;
265  funParam[1]=K3;
266  fGRF->SetParameters(funParam);
267  forigsigma=padDistance;
268  fDSTEPM1 = 10./TMath::Sqrt(padDistance*padDistance+fpadWidth*fpadWidth/12);
269  //by default I set the step as one tenth of sigma
270  //sprintf(fType,"Gati");
271  snprintf(fType,5,"Gati");
272 }
273 
274 
275 
276 void AliTPCRF1D::DrawRF(Float_t x1,Float_t x2,Int_t N)
277 {
278  //
279  //Draw prf in selected region <x1,x2> with nuber of diviision = n
280  //
281  char s[100];
282  TCanvas * c1 = new TCanvas("canRF","Pad response function",700,900);
283  c1->cd();
284  TPad * pad1 = new TPad("pad1RF","",0.05,0.55,0.95,0.95,21);
285  pad1->Draw();
286  TPad * pad2 = new TPad("pad2RF","",0.05,0.05,0.95,0.45,21);
287  pad2->Draw();
288 
289  //sprintf(s,"RF response function for %1.2f cm pad width",
290  // fpadWidth);
291  snprintf(s,60,"RF response function for %1.2f cm pad width",fpadWidth);
292  pad1->cd();
293  TH1F * hRFo = new TH1F("hRFo","Original charge distribution",N+1,x1,x2);
294  pad2->cd();
295  gStyle->SetOptFit(1);
296  gStyle->SetOptStat(0);
297  TH1F * hRFc = new TH1F("hRFc",s,N+1,x1,x2);
298  Float_t x=x1;
299  Float_t y1;
300  Float_t y2;
301 
302  for (Float_t i = 0;i<N+1;i++)
303  {
304  x+=(x2-x1)/Float_t(N);
305  y1 = GetRF(x);
306  hRFc->Fill(x,y1);
307  y2 = GetGRF(x);
308  hRFo->Fill(x,y2);
309  };
310  pad1->cd();
311  hRFo->Fit("gaus");
312  pad2->cd();
313  hRFc->Fit("gaus");
314 }
315 
316 void AliTPCRF1D::Update()
317 {
318  //
319  //update fields with interpolated values for
320  //PRF calculation
321 
322  //at the begining initialize to 0
323  for (Int_t i =0; i<fNRF;i++) fcharge[i] = 0;
324  if ( fGRF == 0 ) return;
325  // This form is no longer available
326 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
327  fInteg = fGRF->Integral(-5*forigsigma,5*forigsigma,funParam,0.00001);
328 #else
329  TArrayD savParam(fGRF->GetNpar(), fGRF->GetParameters());
330  fGRF->SetParameters(funParam);
331  fInteg = fGRF->Integral(-5*forigsigma,5*forigsigma,0.00001);
332 #endif
333  if ( fInteg == 0 ) fInteg = 1;
334  if (fDirect==kFALSE){
335  //integrate charge over pad for different distance of pad
336  for (Int_t i =0; i<fNRF;i++)
337  { //x in cm fpadWidth in cm
338  Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
339  Float_t x1=TMath::Max(x-fpadWidth/2,-5*forigsigma);
340  Float_t x2=TMath::Min(x+fpadWidth/2,5*forigsigma);
341 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
342  fcharge[i] = fkNorm*fGRF->Integral(x1,x2,funParam,0.0001)/fInteg;
343 #else
344  fcharge[i] = fkNorm*fGRF->Integral(x1,x2,0.0001)/fInteg;
345 #endif
346  };
347  }
348  else{
349  for (Int_t i =0; i<fNRF;i++)
350  { //x in cm fpadWidth in cm
351  Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
352  fcharge[i] = fkNorm*fGRF->Eval(x);
353  };
354  }
355  fSigma = 0;
356  Float_t sum =0;
357  Float_t mean=0;
358  for (Float_t x =-fNRF/fDSTEPM1; x<fNRF/fDSTEPM1;x+=1/fDSTEPM1)
359  { //x in cm fpadWidth in cm
360  Float_t weight = GetRF(x+fOffset);
361  fSigma+=x*x*weight;
362  mean+=x*weight;
363  sum+=weight;
364  };
365  if (sum>0){
366  mean/=sum;
367  fSigma = TMath::Sqrt(fSigma/sum-mean*mean);
368  }
369  else fSigma=0;
370 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
371  fGRF->SetParameters(savParam.GetArray());
372 #endif
373 }
374 #if 0
375 void AliTPCRF1D::Streamer(TBuffer &R__b)
376 {
377  // Stream an object of class AliTPCRF1D.
378  if (R__b.IsReading()) {
379  AliTPCRF1D::Class()->ReadBuffer(R__b, this);
380  //read functions
381 
382  if (strncmp(fType,"Gauss",3)==0) {delete fGRF; fGRF = new TF1("funGauss",funGauss,-5.,5.,4);}
383  if (strncmp(fType,"Cosh",3)==0) {delete fGRF; fGRF = new TF1("funCosh",funCosh,-5.,5.,4);}
384  if (strncmp(fType,"Gati",3)==0) {delete fGRF; fGRF = new TF1("funGati",funGati,-5.,5.,4);}
385  if (fGRF) fGRF->SetParameters(funParam);
386 
387  } else {
388  AliTPCRF1D::Class()->WriteBuffer(R__b, this);
389  }
390 }
391 
392 #endif
393 Double_t AliTPCRF1D::Gamma4(Double_t x, Double_t p0, Double_t p1){
394  //
395  // Gamma 4 Time response function of ALTRO
396  // Marian Ivanov: 01/20/2022, Time using in time bins (bin=100 ns)
397  // TF1 f1("f1","AliTPCRF1D::Gamma4((x+1.98094355865156)*0.1,55,0.160)",-5,5)
398  //
399  if (x<0) return 0;
400  Double_t g1 = TMath::Exp(-4.*x/p1);
401  Double_t g2 = TMath::Power(x/p1,4);
402  return p0*g1*g2;
403 }
404 
AliTPCRF1D
Declaration of class AliTPCRF1D.
Definition: AliTPCRF1D.h:19
AliTPCRF1D::Update
void Update()
it&#39;s on user !!!!
Definition: AliTPCRF1D.cxx:316
Generated by   doxygen 1.8.5