Formatting changes.
[u/mrichter/AliRoot.git] / TPC / 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 /* $Id$ */
17
18
19 //-----------------------------------------------------------------------------
20 //
21 //
22 //
23 //  Origin:   Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk
24 //
25 //  Declaration of class AliTPCRF1D
26 //
27 //-----------------------------------------------------------------------------
28
29 //
30
31 #include "TMath.h"
32 #include "AliTPCRF1D.h"
33 #include "TF2.h"
34 #include <Riostream.h>
35 #include "TCanvas.h"
36 #include "TPad.h"
37 #include "TStyle.h"
38 #include "TH1.h"
39 #include <TString.h>
40
41 extern TStyle * gStyle; 
42
43 Int_t   AliTPCRF1D::fgNRF=100;  //default  number of interpolation points
44 Float_t AliTPCRF1D::fgRFDSTEP=0.01; //default step in cm
45
46 static Double_t funGauss(Double_t *x, Double_t * par)
47 {
48   //Gauss function  -needde by the generic function object 
49   return TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]));
50 }
51
52 static Double_t funCosh(Double_t *x, Double_t * par)
53 {
54   //Cosh function  -needde by the generic function object 
55   return 1/TMath::CosH(3.14159*x[0]/(2*par[0]));  
56 }    
57
58 static Double_t funGati(Double_t *x, Double_t * par)
59 {
60   //Gati function  -needde by the generic function object 
61   Float_t k3=par[1];
62   Float_t k3R=TMath::Sqrt(k3);
63   Float_t k2=(TMath::Pi()/2)*(1-k3R/2.);
64   Float_t k1=k2*k3R/(4*TMath::ATan(k3R));
65   Float_t l=x[0]/par[0];
66   Float_t tan2=TMath::TanH(k2*l);
67   tan2*=tan2;
68   Float_t res = k1*(1-tan2)/(1+k3*tan2);  
69   return res;  
70 }    
71
72 ///////////////////////////////////////////////////////////////////////////
73 ///////////////////////////////////////////////////////////////////////////
74
75 ClassImp(AliTPCRF1D)
76
77
78 AliTPCRF1D::AliTPCRF1D(Bool_t direct,Int_t np,Float_t step)
79            :TObject(),
80             fNRF(0),
81             fDSTEPM1(0.),
82             fcharge(0),
83             forigsigma(0.),
84             fpadWidth(3.5),
85             fkNorm(0.5),
86             fInteg(0.),
87             fGRF(0),
88             fSigma(0.),
89             fOffset(0.),
90             fDirect(kFALSE),
91             fPadDistance(0.)
92 {
93   //default constructor for response function object
94   fDirect=direct;
95   if (np!=0)fNRF = np;
96   else (fNRF=fgNRF);
97   fcharge = new Float_t[fNRF];
98   if (step>0) fDSTEPM1=1./step;
99   else fDSTEPM1 = 1./fgRFDSTEP;
100 }
101
102 AliTPCRF1D::AliTPCRF1D(const AliTPCRF1D &prf)
103            :TObject(prf),
104             fNRF(0),
105             fDSTEPM1(0.),
106             fcharge(0),
107             forigsigma(0.),
108             fpadWidth(3.5),
109             fkNorm(0.5),
110             fInteg(0.),
111             fGRF(0),
112             fSigma(0.),
113             fOffset(0.),
114             fDirect(kFALSE),
115             fPadDistance(0.)
116 {
117   
118   //
119   memcpy(this, &prf, sizeof(prf)); 
120   fcharge = new Float_t[fNRF];
121   memcpy(fcharge,prf.fcharge, fNRF);
122   fGRF = new TF1(*(prf.fGRF));
123   //PH Change the name (add 0 to the end)
124   TString s(fGRF->GetName());
125   s+="0";
126   fGRF->SetName(s.Data());
127 }
128
129 AliTPCRF1D & AliTPCRF1D::operator = (const AliTPCRF1D &prf)
130 {
131   //  
132   if (fcharge) delete fcharge;
133   if (fGRF) delete fGRF;
134   memcpy(this, &prf, sizeof(prf)); 
135   fcharge = new Float_t[fNRF];
136   memcpy(fcharge,prf.fcharge, fNRF);
137   fGRF = new TF1(*(prf.fGRF));
138    //PH Change the name (add 0 to the end)
139   TString s(fGRF->GetName());
140   s+="0";
141   fGRF->SetName(s.Data());
142   return (*this);
143 }
144
145
146
147 AliTPCRF1D::~AliTPCRF1D()
148 {
149   //
150   if (fcharge!=0) delete [] fcharge;
151   if (fGRF !=0 ) fGRF->Delete();
152 }
153
154 Float_t AliTPCRF1D::GetRF(Float_t xin)
155 {
156   //function which return response
157   //for the charge in distance xin 
158   //return linear aproximation of RF
159   Float_t x = TMath::Abs((xin-fOffset)*fDSTEPM1)+fNRF/2;
160   Int_t i1=Int_t(x);
161   if (x<0) i1-=1;
162   Float_t res=0;
163   if (i1+1<fNRF)
164     res = fcharge[i1]*(Float_t(i1+1)-x)+fcharge[i1+1]*(x-Float_t(i1));    
165   return res;
166 }
167
168 Float_t  AliTPCRF1D::GetGRF(Float_t xin)
169 {  
170   //function which returnoriginal charge distribution
171   //this function is just normalised for fKnorm
172   if (fGRF != 0 ) 
173     return fkNorm*fGRF->Eval(xin)/fInteg;
174       else
175     return 0.;
176 }
177
178    
179 void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
180                        Float_t kNorm, Float_t sigma)
181 {
182   //adjust parameters of the original charge distribution
183   //and pad size parameters
184    fpadWidth = padwidth;
185    fGRF = GRF;
186    fkNorm = kNorm;
187    if (sigma==0) sigma= fpadWidth/TMath::Sqrt(12.);
188    forigsigma=sigma;
189    fDSTEPM1 = 10/TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
190    sprintf(fType,"User");
191    //   Update();   
192 }
193   
194
195 void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
196                       Float_t kNorm)
197 {
198   // 
199   // set parameters for Gauss generic charge distribution
200   //
201   fpadWidth = padWidth;
202   fkNorm = kNorm;
203   if (fGRF !=0 ) fGRF->Delete();
204   fGRF = new TF1("funGauss",funGauss,-5,5,1);
205   funParam[0]=sigma;
206   forigsigma=sigma;
207   fGRF->SetParameters(funParam);
208    fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
209   //by default I set the step as one tenth of sigma  
210   sprintf(fType,"Gauss");
211 }
212
213 void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
214                      Float_t kNorm)
215 {
216   // 
217   // set parameters for Cosh generic charge distribution
218   //
219   fpadWidth = padWidth;
220   fkNorm = kNorm;
221   if (fGRF !=0 ) fGRF->Delete();
222   fGRF = new TF1("funCosh",     funCosh, -5.,5.,2);   
223   funParam[0]=sigma;
224   fGRF->SetParameters(funParam);
225   forigsigma=sigma;
226   fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
227   //by default I set the step as one tenth of sigma
228   sprintf(fType,"Cosh");
229 }
230
231 void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
232                      Float_t kNorm)
233 {
234   // 
235   // set parameters for Gati generic charge distribution
236   //
237   fpadWidth = padWidth;
238   fkNorm = kNorm;
239   if (fGRF !=0 ) fGRF->Delete();
240   fGRF = new TF1("funGati",     funGati, -5.,5.,2);   
241   funParam[0]=padDistance;
242   funParam[1]=K3;  
243   fGRF->SetParameters(funParam);
244   forigsigma=padDistance;
245   fDSTEPM1 = 10./TMath::Sqrt(padDistance*padDistance+fpadWidth*fpadWidth/12); 
246   //by default I set the step as one tenth of sigma
247   sprintf(fType,"Gati");
248 }
249
250 void AliTPCRF1D::DrawRF(Float_t x1,Float_t x2,Int_t N)
251
252   //
253   //Draw prf in selected region <x1,x2> with nuber of diviision = n
254   //
255   char s[100];
256   TCanvas  * c1 = new TCanvas("canRF","Pad response function",700,900);
257   c1->cd();
258   TPad * pad1 = new TPad("pad1RF","",0.05,0.55,0.95,0.95,21);
259   pad1->Draw();
260   TPad * pad2 = new TPad("pad2RF","",0.05,0.05,0.95,0.45,21);
261   pad2->Draw();
262
263   sprintf(s,"RF response function for %1.2f cm pad width",
264           fpadWidth);  
265   pad1->cd();
266   TH1F * hRFo = new TH1F("hRFo","Original charge distribution",N+1,x1,x2);
267   pad2->cd();
268    gStyle->SetOptFit(1);
269    gStyle->SetOptStat(0); 
270   TH1F * hRFc = new TH1F("hRFc",s,N+1,x1,x2);
271   Float_t x=x1;
272   Float_t y1;
273   Float_t y2;
274
275   for (Float_t i = 0;i<N+1;i++)
276     {
277       x+=(x2-x1)/Float_t(N);
278       y1 = GetRF(x);
279       hRFc->Fill(x,y1);
280       y2 = GetGRF(x);
281       hRFo->Fill(x,y2);      
282     };
283   pad1->cd();
284   hRFo->Fit("gaus");
285   pad2->cd();
286   hRFc->Fit("gaus");
287 }
288
289 void AliTPCRF1D::Update()
290 {
291   //
292   //update fields  with interpolated values for
293   //PRF calculation
294
295   //at the begining initialize to 0
296   for (Int_t i =0; i<fNRF;i++)  fcharge[i] = 0;
297   if ( fGRF == 0 ) return;
298   fInteg  = fGRF->Integral(-5*forigsigma,5*forigsigma,funParam,0.00001);
299   if ( fInteg == 0 ) fInteg = 1; 
300   if (fDirect==kFALSE){
301   //integrate charge over pad for different distance of pad
302   for (Int_t i =0; i<fNRF;i++)
303     {      //x in cm fpadWidth in cm
304       Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
305       Float_t x1=TMath::Max(x-fpadWidth/2,-5*forigsigma);
306       Float_t x2=TMath::Min(x+fpadWidth/2,5*forigsigma);
307       fcharge[i] = 
308         fkNorm*fGRF->Integral(x1,x2,funParam,0.0001)/fInteg;
309     };   
310   }
311   else{
312     for (Int_t i =0; i<fNRF;i++)
313       {      //x in cm fpadWidth in cm
314         Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
315         fcharge[i] = fkNorm*fGRF->Eval(x);
316       };   
317   }  
318   fSigma = 0; 
319   Float_t sum =0;
320   Float_t mean=0;
321   for (Float_t  x =-fNRF/fDSTEPM1; x<fNRF/fDSTEPM1;x+=1/fDSTEPM1)
322     {      //x in cm fpadWidth in cm
323       Float_t weight = GetRF(x+fOffset);
324       fSigma+=x*x*weight; 
325       mean+=x*weight;
326       sum+=weight;
327     };  
328   if (sum>0){
329     mean/=sum;
330     fSigma = TMath::Sqrt(fSigma/sum-mean*mean);   
331   }
332   else fSigma=0; 
333 }
334
335 void AliTPCRF1D::Streamer(TBuffer &R__b)
336 {
337    // Stream an object of class AliTPCRF1D.
338    if (R__b.IsReading()) {
339       AliTPCRF1D::Class()->ReadBuffer(R__b, this);
340       //read functions
341  
342       if (strncmp(fType,"Gauss",3)==0) {delete fGRF; fGRF = new TF1("funGauss",funGauss,-5.,5.,4);}
343       if (strncmp(fType,"Cosh",3)==0)  {delete fGRF; fGRF = new TF1("funCosh",funCosh,-5.,5.,4);}
344       if (strncmp(fType,"Gati",3)==0)  {delete fGRF; fGRF = new TF1("funGati",funGati,-5.,5.,4);}  
345       if (fGRF) fGRF->SetParameters(funParam);     
346
347    } else {
348       AliTPCRF1D::Class()->WriteBuffer(R__b, this);
349    }
350 }
351