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