New TPC files from M.Kowalski.
[u/mrichter/AliRoot.git] / TPC / AliTPCRF1D.cxx
CommitLineData
8c555625 1//-----------------------------------------------------------------------------
2// $Header$
3//
4//
5// Origin: Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk
6//
7// Declaration of class AliTPCRF1D
8//
9//-----------------------------------------------------------------------------
10#include "TMath.h"
11#include "AliTPCRF1D.h"
12#include "TF2.h"
13#include <iostream.h>
14#include "TCanvas.h"
15#include "TPad.h"
16#include "TStyle.h"
17#include "TH1.h"
18
19extern TStyle * gStyle;
20
21static Double_t funGauss(Double_t *x, Double_t * par)
22{
23 return TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]));
24}
25
26static Double_t funCosh(Double_t *x, Double_t * par)
27{
28 return 1/TMath::CosH(3.14159*x[0]/(2*par[0]));
29}
30
31static Double_t funGati(Double_t *x, Double_t * par)
32{
33 //par[1] = is equal to k3
34 //par[0] is equal to pad wire distance
35 Float_t K3=par[1];
36 Float_t K3R=TMath::Sqrt(K3);
37 Float_t K2=(TMath::Pi()/2)*(1-K3R/2.);
38 Float_t K1=K2*K3R/(4*TMath::ATan(K3R));
39 Float_t l=x[0]/par[0];
40 Float_t tan2=TMath::TanH(K2*l);
41 tan2*=tan2;
42 Float_t res = K1*(1-tan2)/(1+K3*tan2);
43 return res;
44}
45
46
47
48
49///////////////////////////////////////////////////////////////////////////
50///////////////////////////////////////////////////////////////////////////
51///////////////////////////////////////////////////////////////////////////
52///////////////////////////////////////////////////////////////////////////
53
54AliTPCRF1D * gRF1D;
55ClassImp(AliTPCRF1D)
56
57
58AliTPCRF1D::AliTPCRF1D(Bool_t direct,Int_t np,Float_t step)
59{
60 fDirect=direct;
61 fNRF = np;
62 fcharge = new Float_t[fNRF];
63 fDSTEPM1=1./step;
64 fSigma = 0;
65 gRF1D = this;
66 fGRF = 0;
67 fkNorm = 0.5;
68 fpadWidth = 3.5;
69 forigsigma=0.;
70 fOffset = 0.;
71}
72
73
74AliTPCRF1D::~AliTPCRF1D()
75{
76 if (fcharge!=0) delete fcharge;
77 if (fGRF !=0 ) fGRF->Delete();
78}
79
80Float_t AliTPCRF1D::GetRF(Float_t xin)
81{
82 //x xin DSTEP unit
83 //return linear aproximation of RF
84 Float_t x = TMath::Abs((xin-fOffset)*fDSTEPM1)+fNRF/2;
85 Int_t i1=Int_t(x);
86 if (x<0) i1-=1;
87 Float_t res=0;
88 if (i1+1<fNRF)
89 res = fcharge[i1]*(Float_t(i1+1)-x)+fcharge[i1+1]*(x-Float_t(i1));
90 return res;
91}
92
93Float_t AliTPCRF1D::GetGRF(Float_t xin)
94{
95 if (fGRF != 0 )
96 return fkNorm*fGRF->Eval(xin)/fInteg;
97 else
98 return 0.;
99}
100
101
102void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
103 Float_t kNorm, Float_t sigma)
104{
105 fpadWidth = padwidth;
106 fGRF = GRF;
107 fkNorm = kNorm;
108 if (sigma==0) sigma= fpadWidth/TMath::Sqrt(12.);
109 forigsigma=sigma;
110 fDSTEPM1 = 10/TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
111 sprintf(fType,"User");
112 // Update();
113}
114
115
116void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
117 Float_t kNorm)
118{
119 // char s[120];
120 fpadWidth = padWidth;
121 fkNorm = kNorm;
122 if (fGRF !=0 ) fGRF->Delete();
123 fGRF = new TF1("fun",funGauss,-5,5,2);
124 funParam[0]=sigma;
125 forigsigma=sigma;
126 fGRF->SetParameters(funParam);
127 fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
128 //by default I set the step as one tenth of sigma
129 // Update();
130 sprintf(fType,"Gauss");
131}
132
133void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
134 Float_t kNorm)
135{
136 // char s[120];
137 fpadWidth = padWidth;
138 fkNorm = kNorm;
139 if (fGRF !=0 ) fGRF->Delete();
140 fGRF = new TF1("fun", funCosh, -5.,5.,2);
141 funParam[0]=sigma;
142 fGRF->SetParameters(funParam);
143 forigsigma=sigma;
144 fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
145 //by default I set the step as one tenth of sigma
146 // Update();
147 sprintf(fType,"Cosh");
148}
149
150void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
151 Float_t kNorm)
152{
153 // char s[120];
154 fpadWidth = padWidth;
155 fkNorm = kNorm;
156 if (fGRF !=0 ) fGRF->Delete();
157 fGRF = new TF1("fun", funGati, -5.,5.,2);
158 funParam[0]=padDistance;
159 funParam[1]=K3;
160 fGRF->SetParameters(funParam);
161 forigsigma=padDistance;
162 fDSTEPM1 = 10./TMath::Sqrt(padDistance*padDistance+fpadWidth*fpadWidth/12);
163 //by default I set the step as one tenth of sigma
164 // Update();
165 sprintf(fType,"Gati");
166}
167
168void AliTPCRF1D::Draw(Float_t x1,Float_t x2,Int_t N)
169{
170 char s[100];
171 TCanvas * c1 = new TCanvas("canRF","Pad response function",700,900);
172 c1->cd();
173 TPad * pad1 = new TPad("pad1RF","",0.05,0.55,0.95,0.95,21);
174 pad1->Draw();
175 TPad * pad2 = new TPad("pad2RF","",0.05,0.05,0.95,0.45,21);
176 pad2->Draw();
177
178 sprintf(s,"RF response function for %1.2f cm pad width",
179 fpadWidth);
180 pad1->cd();
181 TH1F * hRFo = new TH1F("hRFo","Original charge distribution",N+1,x1,x2);
182 pad2->cd();
183 gStyle->SetOptFit(1);
184 gStyle->SetOptStat(0);
185 TH1F * hRFc = new TH1F("hRFc",s,N+1,x1,x2);
186 Float_t x=x1;
187 Float_t y1;
188 Float_t y2;
189
190 for (Float_t i = 0;i<N+1;i++)
191 {
192 x+=(x2-x1)/Float_t(N);
193 y1 = GetRF(x);
194 hRFc->Fill(x,y1);
195 y2 = GetGRF(x);
196 hRFo->Fill(x,y2);
197 };
198 pad1->cd();
199 hRFo->Fit("gaus");
200 pad2->cd();
201 hRFc->Fit("gaus");
202}
203
204void AliTPCRF1D::Update()
205{
206 //initialize to 0
207 for (Int_t i =0; i<fNRF;i++) fcharge[i] = 0;
208 if ( fGRF == 0 ) return;
209 fInteg = fGRF->Integral(-5*forigsigma,5*forigsigma,funParam,0.00001);
210 if ( fInteg == 0 ) fInteg = 1;
211 if (fDirect==kFALSE){
212 //integrate charge over pad for different distance of pad
213 for (Int_t i =0; i<fNRF;i++)
214 { //x in cm fpadWidth in cm
215 Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
216 Float_t x1=TMath::Max(x-fpadWidth/2,-5*forigsigma);
217 Float_t x2=TMath::Min(x+fpadWidth/2,5*forigsigma);
218 fcharge[i] =
219 fkNorm*fGRF->Integral(x1,x2,funParam,0.0001)/fInteg;
220 };
221 }
222 else{
223 for (Int_t i =0; i<fNRF;i++)
224 { //x in cm fpadWidth in cm
225 Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
226 fcharge[i] = fkNorm*fGRF->Eval(x);
227 };
228 }
229 fSigma = 0;
230 Float_t sum =0;
231 Float_t mean=0;
232 for (Float_t x =-fNRF/fDSTEPM1; x<fNRF/fDSTEPM1;x+=1/fDSTEPM1)
233 { //x in cm fpadWidth in cm
234 Float_t weight = GetRF(x+fOffset);
235 fSigma+=x*x*weight;
236 mean+=x*weight;
237 sum+=weight;
238 };
239 if (sum>0){
240 mean/=sum;
241 fSigma = TMath::Sqrt(fSigma/sum-mean*mean);
242 }
243 else fSigma=0;
244}
245
246void AliTPCRF1D::Streamer(TBuffer &R__b)
247{
248 // Stream an object of class AliTPC.
249
250 if (R__b.IsReading()) {
251 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
252 TObject::Streamer(R__b);
253 //read pad parameters
254 R__b >> fpadWidth;
255 //read charge parameters
256 R__b >> fType[0];
257 R__b >> fType[1];
258 R__b >> fType[2];
259 R__b >> fType[3];
260 R__b >> fType[4];
261 R__b >> forigsigma;
262 R__b >> fkNorm;
263 R__b >> fK3X;
264 R__b >> fPadDistance;
265 R__b >> fInteg;
266 R__b >> fOffset;
267 //read functions
268 if (fGRF!=0) {
269 delete fGRF;
270 fGRF=0;
271 }
272 if (strncmp(fType,"User",3)==0){
273 fGRF= new TF1;
274 R__b>>fGRF;
275 }
276
277 if (strncmp(fType,"Gauss",3)==0)
278 fGRF = new TF1("fun",funGauss,-5.,5.,4);
279 if (strncmp(fType,"Cosh",3)==0)
280 fGRF = new TF1("fun",funCosh,-5.,5.,4);
281 if (strncmp(fType,"Gati",3)==0)
282 fGRF = new TF1("fun",funGati,-5.,5.,4);
283 R__b >>fDSTEPM1;
284 R__b >>fNRF;
285 R__b.ReadFastArray(fcharge,fNRF);
286 R__b.ReadFastArray(funParam,5);
287 if (fGRF!=0) fGRF->SetParameters(funParam);
288
289 } else {
290 R__b.WriteVersion(AliTPCRF1D::IsA());
291 TObject::Streamer(R__b);
292 //write pad width
293 R__b << fpadWidth;
294 //write charge parameters
295 R__b << fType[0];
296 R__b << fType[1];
297 R__b << fType[2];
298 R__b << fType[3];
299 R__b << fType[4];
300 R__b << forigsigma;
301 R__b << fkNorm;
302 R__b << fK3X;
303 R__b << fPadDistance;
304 R__b << fInteg;
305 R__b << fOffset;
306 //write interpolation parameters
307 if (strncmp(fType,"User",3)==0) R__b <<fGRF;
308 R__b <<fDSTEPM1;
309 R__b <<fNRF;
310 R__b.WriteFastArray(fcharge,fNRF);
311 R__b.WriteFastArray(funParam,5);
312
313
314
315 }
316}
317