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