Do not hide variables from the outer scope (Solaris x86)
[u/mrichter/AliRoot.git] / TPC / AliTPCRF1D.cxx
CommitLineData
4c039060 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
88cb7938 16/* $Id$ */
19364939 17
4c039060 18
8c555625 19//-----------------------------------------------------------------------------
73042f01 20//
8c555625 21//
22//
23// Origin: Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk
24//
25// Declaration of class AliTPCRF1D
26//
27//-----------------------------------------------------------------------------
cc80f89e 28
73042f01 29//
30
8c555625 31#include "TMath.h"
32#include "AliTPCRF1D.h"
33#include "TF2.h"
19364939 34#include <Riostream.h>
8c555625 35#include "TCanvas.h"
36#include "TPad.h"
37#include "TStyle.h"
38#include "TH1.h"
2a336e15 39#include <TString.h>
8c555625 40
73042f01 41extern TStyle * gStyle;
42
43Int_t AliTPCRF1D::fgNRF=100; //default number of interpolation points
44Float_t AliTPCRF1D::fgRFDSTEP=0.01; //default step in cm
8c555625 45
46static Double_t funGauss(Double_t *x, Double_t * par)
47{
cc80f89e 48 //Gauss function -needde by the generic function object
8c555625 49 return TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]));
50}
51
52static Double_t funCosh(Double_t *x, Double_t * par)
53{
cc80f89e 54 //Cosh function -needde by the generic function object
8c555625 55 return 1/TMath::CosH(3.14159*x[0]/(2*par[0]));
56}
57
58static Double_t funGati(Double_t *x, Double_t * par)
59{
cc80f89e 60 //Gati function -needde by the generic function object
73042f01 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));
8c555625 65 Float_t l=x[0]/par[0];
73042f01 66 Float_t tan2=TMath::TanH(k2*l);
8c555625 67 tan2*=tan2;
73042f01 68 Float_t res = k1*(1-tan2)/(1+k3*tan2);
8c555625 69 return res;
70}
71
8c555625 72///////////////////////////////////////////////////////////////////////////
73///////////////////////////////////////////////////////////////////////////
74
8c555625 75ClassImp(AliTPCRF1D)
76
77
78AliTPCRF1D::AliTPCRF1D(Bool_t direct,Int_t np,Float_t step)
179c6296 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.)
8c555625 92{
cc80f89e 93 //default constructor for response function object
8c555625 94 fDirect=direct;
73042f01 95 if (np!=0)fNRF = np;
96 else (fNRF=fgNRF);
8c555625 97 fcharge = new Float_t[fNRF];
73042f01 98 if (step>0) fDSTEPM1=1./step;
99 else fDSTEPM1 = 1./fgRFDSTEP;
8c555625 100}
101
179c6296 102AliTPCRF1D::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.)
73042f01 116{
176aff27 117
73042f01 118 //
119 memcpy(this, &prf, sizeof(prf));
120 fcharge = new Float_t[fNRF];
121 memcpy(fcharge,prf.fcharge, fNRF);
2a336e15 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());
73042f01 127}
128
129AliTPCRF1D & 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));
2a336e15 138 //PH Change the name (add 0 to the end)
139 TString s(fGRF->GetName());
140 s+="0";
141 fGRF->SetName(s.Data());
73042f01 142 return (*this);
143}
144
145
8c555625 146
147AliTPCRF1D::~AliTPCRF1D()
148{
73042f01 149 //
cc80f89e 150 if (fcharge!=0) delete [] fcharge;
8c555625 151 if (fGRF !=0 ) fGRF->Delete();
152}
153
154Float_t AliTPCRF1D::GetRF(Float_t xin)
155{
cc80f89e 156 //function which return response
157 //for the charge in distance xin
8c555625 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
168Float_t AliTPCRF1D::GetGRF(Float_t xin)
169{
cc80f89e 170 //function which returnoriginal charge distribution
171 //this function is just normalised for fKnorm
8c555625 172 if (fGRF != 0 )
173 return fkNorm*fGRF->Eval(xin)/fInteg;
174 else
175 return 0.;
176}
177
178
179void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
180 Float_t kNorm, Float_t sigma)
181{
cc80f89e 182 //adjust parameters of the original charge distribution
183 //and pad size parameters
8c555625 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
195void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
196 Float_t kNorm)
197{
cc80f89e 198 //
199 // set parameters for Gauss generic charge distribution
200 //
8c555625 201 fpadWidth = padWidth;
202 fkNorm = kNorm;
203 if (fGRF !=0 ) fGRF->Delete();
2a336e15 204 fGRF = new TF1("funGauss",funGauss,-5,5,1);
8c555625 205 funParam[0]=sigma;
206 forigsigma=sigma;
207 fGRF->SetParameters(funParam);
208 fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
cc80f89e 209 //by default I set the step as one tenth of sigma
8c555625 210 sprintf(fType,"Gauss");
211}
212
213void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
214 Float_t kNorm)
215{
cc80f89e 216 //
217 // set parameters for Cosh generic charge distribution
218 //
8c555625 219 fpadWidth = padWidth;
220 fkNorm = kNorm;
221 if (fGRF !=0 ) fGRF->Delete();
2a336e15 222 fGRF = new TF1("funCosh", funCosh, -5.,5.,2);
8c555625 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
8c555625 228 sprintf(fType,"Cosh");
229}
230
231void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
232 Float_t kNorm)
233{
cc80f89e 234 //
235 // set parameters for Gati generic charge distribution
236 //
8c555625 237 fpadWidth = padWidth;
238 fkNorm = kNorm;
239 if (fGRF !=0 ) fGRF->Delete();
2a336e15 240 fGRF = new TF1("funGati", funGati, -5.,5.,2);
8c555625 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
8c555625 247 sprintf(fType,"Gati");
248}
249
73042f01 250void AliTPCRF1D::DrawRF(Float_t x1,Float_t x2,Int_t N)
8c555625 251{
cc80f89e 252 //
253 //Draw prf in selected region <x1,x2> with nuber of diviision = n
254 //
8c555625 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
289void AliTPCRF1D::Update()
290{
cc80f89e 291 //
292 //update fields with interpolated values for
293 //PRF calculation
294
295 //at the begining initialize to 0
8c555625 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
335void AliTPCRF1D::Streamer(TBuffer &R__b)
336{
cc80f89e 337 // Stream an object of class AliTPCRF1D.
8c555625 338 if (R__b.IsReading()) {
2ab0c725 339 AliTPCRF1D::Class()->ReadBuffer(R__b, this);
8c555625 340 //read functions
8c555625 341
2a336e15 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);}
2ab0c725 345 if (fGRF) fGRF->SetParameters(funParam);
8c555625 346
347 } else {
2ab0c725 348 AliTPCRF1D::Class()->WriteBuffer(R__b, this);
8c555625 349 }
350}
351