Message commented out
[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 {
80   //default constructor for response function object
81   fDirect=direct;
82   if (np!=0)fNRF = np;
83   else (fNRF=fgNRF);
84   fcharge = new Float_t[fNRF];
85   if (step>0) fDSTEPM1=1./step;
86   else fDSTEPM1 = 1./fgRFDSTEP;
87   fSigma = 0;
88   fGRF = 0;
89   fkNorm = 0.5;
90   fpadWidth = 3.5;
91   forigsigma=0.;
92   fOffset = 0.;
93 }
94
95 AliTPCRF1D::AliTPCRF1D(const AliTPCRF1D &prf):TObject(prf)
96 {
97   
98   //
99   memcpy(this, &prf, sizeof(prf)); 
100   fcharge = new Float_t[fNRF];
101   memcpy(fcharge,prf.fcharge, fNRF);
102   fGRF = new TF1(*(prf.fGRF));
103   //PH Change the name (add 0 to the end)
104   TString s(fGRF->GetName());
105   s+="0";
106   fGRF->SetName(s.Data());
107 }
108
109 AliTPCRF1D & AliTPCRF1D::operator = (const AliTPCRF1D &prf)
110 {
111   //  
112   if (fcharge) delete fcharge;
113   if (fGRF) delete fGRF;
114   memcpy(this, &prf, sizeof(prf)); 
115   fcharge = new Float_t[fNRF];
116   memcpy(fcharge,prf.fcharge, fNRF);
117   fGRF = new TF1(*(prf.fGRF));
118    //PH Change the name (add 0 to the end)
119   TString s(fGRF->GetName());
120   s+="0";
121   fGRF->SetName(s.Data());
122   return (*this);
123 }
124
125
126
127 AliTPCRF1D::~AliTPCRF1D()
128 {
129   //
130   if (fcharge!=0) delete [] fcharge;
131   if (fGRF !=0 ) fGRF->Delete();
132 }
133
134 Float_t AliTPCRF1D::GetRF(Float_t xin)
135 {
136   //function which return response
137   //for the charge in distance xin 
138   //return linear aproximation of RF
139   Float_t x = TMath::Abs((xin-fOffset)*fDSTEPM1)+fNRF/2;
140   Int_t i1=Int_t(x);
141   if (x<0) i1-=1;
142   Float_t res=0;
143   if (i1+1<fNRF)
144     res = fcharge[i1]*(Float_t(i1+1)-x)+fcharge[i1+1]*(x-Float_t(i1));    
145   return res;
146 }
147
148 Float_t  AliTPCRF1D::GetGRF(Float_t xin)
149 {  
150   //function which returnoriginal charge distribution
151   //this function is just normalised for fKnorm
152   if (fGRF != 0 ) 
153     return fkNorm*fGRF->Eval(xin)/fInteg;
154       else
155     return 0.;
156 }
157
158    
159 void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
160                        Float_t kNorm, Float_t sigma)
161 {
162   //adjust parameters of the original charge distribution
163   //and pad size parameters
164    fpadWidth = padwidth;
165    fGRF = GRF;
166    fkNorm = kNorm;
167    if (sigma==0) sigma= fpadWidth/TMath::Sqrt(12.);
168    forigsigma=sigma;
169    fDSTEPM1 = 10/TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
170    sprintf(fType,"User");
171    //   Update();   
172 }
173   
174
175 void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
176                       Float_t kNorm)
177 {
178   // 
179   // set parameters for Gauss generic charge distribution
180   //
181   fpadWidth = padWidth;
182   fkNorm = kNorm;
183   if (fGRF !=0 ) fGRF->Delete();
184   fGRF = new TF1("funGauss",funGauss,-5,5,1);
185   funParam[0]=sigma;
186   forigsigma=sigma;
187   fGRF->SetParameters(funParam);
188    fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
189   //by default I set the step as one tenth of sigma  
190   sprintf(fType,"Gauss");
191 }
192
193 void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
194                      Float_t kNorm)
195 {
196   // 
197   // set parameters for Cosh generic charge distribution
198   //
199   fpadWidth = padWidth;
200   fkNorm = kNorm;
201   if (fGRF !=0 ) fGRF->Delete();
202   fGRF = new TF1("funCosh",     funCosh, -5.,5.,2);   
203   funParam[0]=sigma;
204   fGRF->SetParameters(funParam);
205   forigsigma=sigma;
206   fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12); 
207   //by default I set the step as one tenth of sigma
208   sprintf(fType,"Cosh");
209 }
210
211 void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
212                      Float_t kNorm)
213 {
214   // 
215   // set parameters for Gati generic charge distribution
216   //
217   fpadWidth = padWidth;
218   fkNorm = kNorm;
219   if (fGRF !=0 ) fGRF->Delete();
220   fGRF = new TF1("funGati",     funGati, -5.,5.,2);   
221   funParam[0]=padDistance;
222   funParam[1]=K3;  
223   fGRF->SetParameters(funParam);
224   forigsigma=padDistance;
225   fDSTEPM1 = 10./TMath::Sqrt(padDistance*padDistance+fpadWidth*fpadWidth/12); 
226   //by default I set the step as one tenth of sigma
227   sprintf(fType,"Gati");
228 }
229
230 void AliTPCRF1D::DrawRF(Float_t x1,Float_t x2,Int_t N)
231
232   //
233   //Draw prf in selected region <x1,x2> with nuber of diviision = n
234   //
235   char s[100];
236   TCanvas  * c1 = new TCanvas("canRF","Pad response function",700,900);
237   c1->cd();
238   TPad * pad1 = new TPad("pad1RF","",0.05,0.55,0.95,0.95,21);
239   pad1->Draw();
240   TPad * pad2 = new TPad("pad2RF","",0.05,0.05,0.95,0.45,21);
241   pad2->Draw();
242
243   sprintf(s,"RF response function for %1.2f cm pad width",
244           fpadWidth);  
245   pad1->cd();
246   TH1F * hRFo = new TH1F("hRFo","Original charge distribution",N+1,x1,x2);
247   pad2->cd();
248    gStyle->SetOptFit(1);
249    gStyle->SetOptStat(0); 
250   TH1F * hRFc = new TH1F("hRFc",s,N+1,x1,x2);
251   Float_t x=x1;
252   Float_t y1;
253   Float_t y2;
254
255   for (Float_t i = 0;i<N+1;i++)
256     {
257       x+=(x2-x1)/Float_t(N);
258       y1 = GetRF(x);
259       hRFc->Fill(x,y1);
260       y2 = GetGRF(x);
261       hRFo->Fill(x,y2);      
262     };
263   pad1->cd();
264   hRFo->Fit("gaus");
265   pad2->cd();
266   hRFc->Fit("gaus");
267 }
268
269 void AliTPCRF1D::Update()
270 {
271   //
272   //update fields  with interpolated values for
273   //PRF calculation
274
275   //at the begining initialize to 0
276   for (Int_t i =0; i<fNRF;i++)  fcharge[i] = 0;
277   if ( fGRF == 0 ) return;
278   fInteg  = fGRF->Integral(-5*forigsigma,5*forigsigma,funParam,0.00001);
279   if ( fInteg == 0 ) fInteg = 1; 
280   if (fDirect==kFALSE){
281   //integrate charge over pad for different distance of pad
282   for (Int_t i =0; i<fNRF;i++)
283     {      //x in cm fpadWidth in cm
284       Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
285       Float_t x1=TMath::Max(x-fpadWidth/2,-5*forigsigma);
286       Float_t x2=TMath::Min(x+fpadWidth/2,5*forigsigma);
287       fcharge[i] = 
288         fkNorm*fGRF->Integral(x1,x2,funParam,0.0001)/fInteg;
289     };   
290   }
291   else{
292     for (Int_t i =0; i<fNRF;i++)
293       {      //x in cm fpadWidth in cm
294         Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
295         fcharge[i] = fkNorm*fGRF->Eval(x);
296       };   
297   }  
298   fSigma = 0; 
299   Float_t sum =0;
300   Float_t mean=0;
301   for (Float_t  x =-fNRF/fDSTEPM1; x<fNRF/fDSTEPM1;x+=1/fDSTEPM1)
302     {      //x in cm fpadWidth in cm
303       Float_t weight = GetRF(x+fOffset);
304       fSigma+=x*x*weight; 
305       mean+=x*weight;
306       sum+=weight;
307     };  
308   if (sum>0){
309     mean/=sum;
310     fSigma = TMath::Sqrt(fSigma/sum-mean*mean);   
311   }
312   else fSigma=0; 
313 }
314
315 void AliTPCRF1D::Streamer(TBuffer &R__b)
316 {
317    // Stream an object of class AliTPCRF1D.
318    if (R__b.IsReading()) {
319       AliTPCRF1D::Class()->ReadBuffer(R__b, this);
320       //read functions
321  
322       if (strncmp(fType,"Gauss",3)==0) {delete fGRF; fGRF = new TF1("funGauss",funGauss,-5.,5.,4);}
323       if (strncmp(fType,"Cosh",3)==0)  {delete fGRF; fGRF = new TF1("funCosh",funCosh,-5.,5.,4);}
324       if (strncmp(fType,"Gati",3)==0)  {delete fGRF; fGRF = new TF1("funGati",funGati,-5.,5.,4);}  
325       if (fGRF) fGRF->SetParameters(funParam);     
326
327    } else {
328       AliTPCRF1D::Class()->WriteBuffer(R__b, this);
329    }
330 }
331