]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCPRF2D.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Updated flags for low flux case (A. Dainese)
[u/mrichter/AliRoot.git] / TPC / AliTPCPRF2D.cxx
... / ...
CommitLineData
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// AliTPCPRF2D - //
20// Pad response function object in two dimesions //
21// This class contains the basic functions for the //
22// calculation of PRF according generic charge distribution //
23// In Update function object calculate table of response function //
24// in discrete x and y position //
25// This table is used for interpolation od response function in any position //
26// (function GetPRF) //
27// //
28// Origin: Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk //
29// //
30///////////////////////////////////////////////////////////////////////////////
31
32#include <Riostream.h>
33#include <TCanvas.h>
34#include <TClass.h>
35#include <TF2.h>
36#include <TH1.h>
37#include <TMath.h>
38#include <TPad.h>
39#include <TPaveText.h>
40#include <TStyle.h>
41#include <TText.h>
42#include <string.h>
43
44#include "AliH2F.h"
45#include "AliTPCPRF2D.h"
46
47
48extern TStyle * gStyle;
49
50const Double_t AliTPCPRF2D::fgkDegtoRad = 0.01745329251994;
51const Double_t AliTPCPRF2D::fgkSQRT12=3.464101;
52const Int_t AliTPCPRF2D::fgkNPRF = 100;
53
54
55static Double_t funGauss2D(Double_t *x, Double_t * par)
56{
57//Gauss function -needde by the generic function object
58 return ( TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]))*
59 TMath::Exp(-(x[1]*x[1])/(2*par[1]*par[1])));
60
61}
62
63static Double_t funCosh2D(Double_t *x, Double_t * par)
64{
65 //Cosh function -needde by the generic function object
66 return ( 1/(TMath::CosH(3.14159*x[0]/(2*par[0]))*
67 TMath::CosH(3.14159*x[1]/(2*par[1]))));
68}
69
70static Double_t funGati2D(Double_t *x, Double_t * par)
71{
72 //Gati function -needde by the generic function object
73 Float_t k3=par[1];
74 Float_t k3R=TMath::Sqrt(k3);
75 Float_t k2=(TMath::Pi()/2)*(1-k3R/2.);
76 Float_t k1=k2*k3R/(4*TMath::ATan(k3R));
77 Float_t l=x[0]/par[0];
78 Float_t tan2=TMath::TanH(k2*l);
79 tan2*=tan2;
80 Float_t res = k1*(1-tan2)/(1+k3*tan2);
81 //par[4] = is equal to k3Y
82 k3=par[4];
83 k3R=TMath::Sqrt(k3);
84 k2=(TMath::Pi()/2)*(1-k3R/2.);
85 k1=k2*k3R/(4*TMath::ATan(k3R));
86 l=x[1]/par[0];
87 tan2=TMath::TanH(k2*l);
88 tan2*=tan2;
89 res = res*k1*(1-tan2)/(1+k3*tan2);
90 return res;
91}
92
93///////////////////////////////////////////////////////////////////////////
94///////////////////////////////////////////////////////////////////////////
95
96ClassImp(AliTPCPRF2D)
97
98AliTPCPRF2D::AliTPCPRF2D()
99 :TObject(),
100 fcharge(0),
101 fY1(0.),
102 fY2(0.),
103 fNYdiv(0),
104 fNChargeArray(0),
105 fChargeArray(0),
106 fHeightFull(0.),
107 fHeightS(0.),
108 fShiftY(0.),
109 fWidth(0.),
110 fK(0.),
111 fNPRF(0),
112 fNdiv(5),
113 fDStep(0.),
114 fKNorm(1.),
115 fInteg(0.),
116 fGRF(0),
117 fK3X(0.),
118 fK3Y(0.),
119 fPadDistance(0.),
120 fOrigSigmaX(0.),
121 fOrigSigmaY(0.),
122 fChargeAngle(0.),
123 fPadAngle(0.),
124 fSigmaX(0.),
125 fSigmaY(0.),
126 fMeanX(0.),
127 fMeanY(0.),
128 fInterX(0),
129 fInterY(0),
130 fCurrentY(0.),
131 fDYtoWire(0.),
132 fDStepM1(0.)
133{
134 //default constructor for response function object
135
136 fNPRF =fgkNPRF ;
137
138 //chewron default values
139 SetPad(0.8,0.8);
140 SetChevron(0.2,0.0,1.0);
141 SetY(-0.2,0.2,2);
142 SetInterpolationType(2,0);
143}
144
145AliTPCPRF2D::~AliTPCPRF2D()
146{
147 if (fChargeArray!=0) delete [] fChargeArray;
148 if (fGRF !=0 ) fGRF->Delete();
149}
150
151void AliTPCPRF2D::SetY(Float_t y1, Float_t y2, Int_t nYdiv)
152{
153 //
154 //set virtual line position
155 //first and last line and number of lines
156 fNYdiv = nYdiv;
157 fY1=y1;
158 fY2=y2;
159}
160
161void AliTPCPRF2D::SetPad(Float_t width, Float_t height)
162{
163 //set base chevron parameters
164 fHeightFull=height;
165 fWidth=width;
166}
167void AliTPCPRF2D::SetChevron(Float_t hstep,
168 Float_t shifty,
169 Float_t fac)
170{
171 //set shaping of chewron parameters
172 fHeightS=hstep;
173 fShiftY=shifty;
174 fK=fac;
175}
176
177void AliTPCPRF2D::SetChParam(Float_t width, Float_t height,
178 Float_t hstep, Float_t shifty, Float_t fac)
179{
180 SetPad(width,height);
181 SetChevron(hstep,shifty,fac);
182}
183
184
185Float_t AliTPCPRF2D::GetPRF(Float_t xin, Float_t yin)
186{
187 //function which return pad response
188 //for the charge in distance xin
189 //return cubic aproximation of PRF or PRF at nearest virtual wire
190 if (fChargeArray==0) return 0;
191 //transform position to "wire position"
192 Float_t y=fDYtoWire*(yin-fY1);
193 if (fNYdiv == 1) y=fY1;
194 //normaly it find nearest line charge
195 if (fInterY ==0){
196 Int_t i=Int_t(0.5+y);
197 if (y<0) i=Int_t(-0.5+y);
198 if ((i<0) || (i>=fNYdiv) ) return 0;
199 fcharge = &(fChargeArray[i*fNPRF]);
200 return GetPRFActiv(xin);
201 }
202 //make interpolation from more fore lines
203 Int_t i= Int_t(y);
204 Float_t res;
205 if ((i<0) || (i>=fNYdiv) ) return 0;
206 Float_t z0=0;
207 Float_t z1=0;
208 Float_t z2=0;
209 Float_t z3=0;
210 if (i>0) {
211 fcharge =&(fChargeArray[(i-1)*fNPRF]);
212 z0 = GetPRFActiv(xin);
213 }
214 fcharge =&(fChargeArray[i*fNPRF]);
215 z1=GetPRFActiv(xin);
216 if ((i+1)<fNYdiv){
217 fcharge =&(fChargeArray[(i+1)*fNPRF]);
218 z2 = GetPRFActiv(xin);
219 }
220 if ((i+2)<fNYdiv){
221 fcharge =&(fChargeArray[(i+2)*fNPRF]);
222 z3 = GetPRFActiv(xin);
223 }
224 Float_t a,b,c,d,k,l;
225 a=z1;
226 b=(z2-z0)/2.;
227 k=z2-a-b;
228 l=(z3-z1)/2.-b;
229 d=l-2*k;
230 c=k-d;
231 Float_t dy=y-Float_t(i);
232
233 res = a+b*dy+c*dy*dy+d*dy*dy*dy;
234 return res;
235}
236
237
238Float_t AliTPCPRF2D::GetPRFActiv(Float_t xin)
239{
240 //GEt response function on given charege line
241 //return spline aproximaton
242 Float_t x = (xin*fDStepM1)+fNPRF/2;
243 Int_t i = Int_t(x);
244
245 if ( (i>1) && ((i+2)<fNPRF)) {
246 Float_t a,b,c,d,k,l;
247 a = fcharge[i];
248 b = (fcharge[i+1]-fcharge[i-1])*0.5;
249 k = fcharge[i+1]-a-b;
250 l = (fcharge[i+2]-fcharge[i])*0.5-b;
251 d=l-2.*k;
252 c=k-d;
253 Float_t dx=x-Float_t(i);
254 Float_t res = a+b*dx+c*dx*dx+d*dx*dx*dx;
255 return res;
256 }
257 else return 0;
258}
259
260
261Float_t AliTPCPRF2D::GetGRF(Float_t xin, Float_t yin)
262{
263 //function which returnoriginal charge distribution
264 //this function is just normalised for fKnorm
265 if (GetGRF() != 0 )
266 return fKNorm*GetGRF()->Eval(xin,yin)/fInteg;
267 else
268 return 0.;
269}
270
271
272void AliTPCPRF2D::SetParam( TF2 * GRF, Float_t kNorm,
273 Float_t sigmaX, Float_t sigmaY)
274{
275 //adjust parameters of the original charge distribution
276 //and pad size parameters
277 if (fGRF !=0 ) fGRF->Delete();
278 fGRF = GRF;
279 fKNorm = kNorm;
280 sprintf(fType,"User");
281 if (sigmaX ==0) sigmaX=(fWidth*(1+TMath::Abs(fK)))/fgkSQRT12;
282 if (sigmaY ==0) sigmaY=(fWidth*(1+TMath::Abs(fK)))/fgkSQRT12;
283 fOrigSigmaX=sigmaX;
284 fOrigSigmaY=sigmaY;
285 Double_t estimsigma =
286 TMath::Sqrt(sigmaX*sigmaX+(fWidth*fWidth*(1+TMath::Abs(fK))/12)+
287 TMath::Tan(fPadAngle*fgkDegtoRad)*TMath::Tan(fPadAngle*fgkDegtoRad)*fHeightFull*fHeightFull/12);
288 if (estimsigma < 5*sigmaX) {
289 fDStep = estimsigma/10.;
290 fNPRF = Int_t(estimsigma*8./fDStep);
291 }
292 else{
293 fDStep = sigmaX;
294 Double_t width = fWidth*(1+TMath::Abs(fK))+TMath::Abs(TMath::Tan(fPadAngle*fgkDegtoRad))*fHeightFull;
295 fNPRF = Int_t((width+8.*sigmaX)/fDStep);
296 };
297
298}
299
300
301void AliTPCPRF2D::SetGauss(Float_t sigmaX, Float_t sigmaY,
302 Float_t kNorm)
303{
304 //
305 // set parameters for Gauss generic charge distribution
306 //
307 fKNorm = kNorm;
308 fOrigSigmaX=sigmaX;
309 fOrigSigmaY=sigmaY;
310 sprintf(fType,"Gauss");
311 if (fGRF !=0 ) fGRF->Delete();
312 fGRF = new TF2("funGauss2D",funGauss2D,-5.,5.,-5.,5.,4);
313
314 funParam[0]=sigmaX;
315 funParam[1]=sigmaY;
316 funParam[2]=fK;
317 funParam[3]=fHeightS;
318
319 fGRF->SetParameters(funParam);
320 Double_t estimsigma =
321 TMath::Sqrt(sigmaX*sigmaX+(fWidth*fWidth*(1+TMath::Abs(fK))/12)+
322 TMath::Tan(fPadAngle)*TMath::Tan(fPadAngle*fgkDegtoRad)*fHeightFull*fHeightFull/12);
323 if (estimsigma < 5*sigmaX) {
324 fDStep = estimsigma/10.;
325 fNPRF = Int_t(estimsigma*8./fDStep);
326 }
327 else{
328 fDStep = sigmaX;
329 Double_t width = fWidth*(1+TMath::Abs(fK))+TMath::Abs(TMath::Tan(fPadAngle*fgkDegtoRad))*fHeightFull;
330 fNPRF = Int_t((width+8.*sigmaX)/fDStep);
331 };
332
333
334}
335void AliTPCPRF2D::SetCosh(Float_t sigmaX, Float_t sigmaY,
336 Float_t kNorm)
337{
338 // set parameters for Cosh generic charge distribution
339 //
340 fKNorm = kNorm;
341 fOrigSigmaX=sigmaX;
342 fOrigSigmaY=sigmaY;
343 sprintf(fType,"Cosh");
344 if (fGRF !=0 ) fGRF->Delete();
345 fGRF = new TF2("funCosh2D", funCosh2D,-5.,5.,-5.,5.,4);
346 funParam[0]=sigmaX;
347 funParam[1]=sigmaY;
348 funParam[2]=fK;
349 funParam[3]=fHeightS;
350 fGRF->SetParameters(funParam);
351
352 Double_t estimsigma = TMath::Sqrt(sigmaX*sigmaX+fWidth*fWidth*(1+TMath::Abs(fK))/12);
353 if (estimsigma < 5*sigmaX) {
354 fDStep = estimsigma/10.;
355 fNPRF = Int_t(estimsigma*8./fDStep);
356 }
357 else{
358 fDStep = sigmaX;
359 fNPRF = Int_t((1.2*fWidth*(1+TMath::Abs(fK))+8.*sigmaX)/fDStep);
360 };
361
362}
363
364void AliTPCPRF2D::SetGati(Float_t K3X, Float_t K3Y,
365 Float_t padDistance,
366 Float_t kNorm)
367{
368 // set parameters for Gati generic charge distribution
369 //
370 fKNorm = kNorm;
371 fK3X=K3X;
372 fK3Y=K3Y;
373 fPadDistance=padDistance;
374 sprintf(fType,"Gati");
375 if (fGRF !=0 ) fGRF->Delete();
376 fGRF = new TF2("funGati2D", funGati2D,-5.,5.,-5.,5.,5);
377
378 funParam[0]=padDistance;
379 funParam[1]=K3X;
380 funParam[2]=fK;
381 funParam[3]=fHeightS;
382 funParam[4]=K3Y;
383 fGRF->SetParameters(funParam);
384 fOrigSigmaX=padDistance;
385 fOrigSigmaY=padDistance;
386 Float_t sigmaX = fOrigSigmaX;
387 Double_t estimsigma = TMath::Sqrt(sigmaX*sigmaX+fWidth*fWidth*(1+TMath::Abs(fK))/12);
388 if (estimsigma < 5*sigmaX) {
389 fDStep = estimsigma/10.;
390 fNPRF = Int_t(estimsigma*8./fDStep);
391 }
392 else{
393 fDStep = sigmaX;
394 fNPRF = Int_t((1.2*fWidth*(1+TMath::Abs(fK))+8.*sigmaX)/fDStep);
395 };
396}
397
398
399
400void AliTPCPRF2D::Update()
401{
402 //
403 //update fields with interpolated values for
404 //PRF calculation
405
406 if ( fGRF == 0 ) return;
407 //initialize interpolated values to 0
408 Int_t i;
409 if (fChargeArray!=0) delete [] fChargeArray;
410 fChargeArray = new Float_t[fNPRF*fNYdiv];
411 fNChargeArray = fNPRF*fNYdiv;
412 for (i =0; i<fNPRF*fNYdiv;i++) fChargeArray[i] = 0;
413 //firstly calculate total integral of charge
414
415 ////////////////////////////////////////////////////////
416 //I'm waiting for normal integral
417 //in this moment only sum
418 Float_t x2= 4*fOrigSigmaX;
419 Float_t y2= 4*fOrigSigmaY;
420 Float_t dx = fOrigSigmaX/Float_t(fNdiv*6);
421 Float_t dy = fOrigSigmaY/Float_t(fNdiv*6);
422 Int_t nx = Int_t(0.5+x2/dx);
423 Int_t ny = Int_t(0.5+y2/dy);
424 Int_t ix,iy;
425 fInteg = 0;
426 Double_t dInteg =0;
427 for (ix=-nx;ix<=nx;ix++)
428 for ( iy=-ny;iy<=ny;iy++)
429 dInteg+=fGRF->Eval(Float_t(ix)*dx,Float_t(iy)*dy)*dx*dy;
430 /////////////////////////////////////////////////////
431 fInteg =dInteg;
432 if ( fInteg == 0 ) fInteg = 1;
433
434 for (i=0; i<fNYdiv; i++){
435 if (fNYdiv == 1) fCurrentY = fY1;
436 else
437 fCurrentY = fY1+Double_t(i)*(fY2-fY1)/Double_t(fNYdiv-1);
438 fcharge = &(fChargeArray[i*fNPRF]);
439 Update1();
440 }
441 //calculate conversion coefitient to convert position to virtual wire
442 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
443 fDStepM1=1/fDStep;
444 UpdateSigma();
445}
446
447void AliTPCPRF2D::Update1()
448{
449 //
450 //update fields with interpolated values for
451 //PRF calculation for given charge line
452 Int_t i;
453 Double_t cos = TMath::Cos(fChargeAngle);
454 Double_t sin = TMath::Sin(fChargeAngle);
455 const Double_t kprec =0.00000001;
456 //integrate charge over pad for different distance of pad
457 for (i =0; i<fNPRF;i++){
458 //x in cm fWidth in cm
459 //calculate integral
460 Double_t xch = fDStep * (Double_t)(i-fNPRF/2);
461 fcharge[i]=0;
462 Double_t k=1;
463
464
465 for (Double_t ym=-fHeightFull/2.-fShiftY; ym<fHeightFull/2.-kprec;ym+=fHeightS){
466 Double_t y2chev=TMath::Min((ym+fHeightS),Double_t(fHeightFull/2.)); // end of chevron step
467 Double_t y1chev= ym; //beginning of chevron step
468 Double_t y2 = TMath::Min(y2chev,fCurrentY+3.5*fOrigSigmaY);
469 Double_t y1 = TMath::Max((y1chev),Double_t(-fHeightFull/2.));
470 y1 = TMath::Max(y1chev,fCurrentY-3.5*fOrigSigmaY);
471
472 Double_t x0 = fWidth*(-1.-(Double_t(k)*fK))*0.5+ym*TMath::Tan(fPadAngle*fgkDegtoRad);
473 Double_t kx = Double_t(k)*(fK*fWidth)/fHeightS;
474 kx = TMath::Tan(TMath::ATan(kx))+TMath::Tan(fPadAngle*fgkDegtoRad);
475
476 Int_t ny = TMath::Max(Int_t(fNdiv*TMath::Exp(-(y1-fCurrentY)*(y1-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),4);
477 Double_t dy = TMath::Min(fOrigSigmaY/Double_t(ny),y2-y1);
478 Double_t ndy = dy;
479
480 //loop over different y strips with variable step size dy
481 if (y2>(y1+kprec)) for (Double_t y = y1; y<y2+kprec;){
482 //new step SIZE
483
484 ny = TMath::Max(Int_t(fNdiv*TMath::Exp(-(y-fCurrentY)*(y-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),5);
485 ndy = fOrigSigmaY/Double_t(ny);
486 if (ndy>(y2-y-dy)) {
487 ndy =y2-y-dy;
488 if (ndy<kprec) ndy=2*kprec; //calculate new delta y
489 }
490 //
491 Double_t sumch=0;
492 //calculation of x borders and initial step
493 Double_t deltay = (y-y1chev);
494
495 Double_t xp1 = x0+deltay*kx;
496 //x begining of pad at position y
497 Double_t xp2 =xp1+fWidth; //x end of pad at position y
498 Double_t xp3 =xp1+kx*dy; //...at position y+dy
499 Double_t xp4 =xp2+kx*dy; //..
500
501 Double_t x1 = TMath::Min(xp1,xp3);
502 x1 = TMath::Max(xp1,xch-3.5*fOrigSigmaX); //beging of integration
503 Double_t x2 = TMath::Max(xp2,xp4);
504 x2 = TMath::Min(xp2+dy*kx,xch+3.5*fOrigSigmaX); //end of integration
505
506 Int_t nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x1-xch)*(x1-xch)/(2*fOrigSigmaX*fOrigSigmaX))*
507 TMath::Exp(-(y1-fCurrentY)*(y1-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),2);
508 Double_t dx = TMath::Min(fOrigSigmaX/Double_t(nx),x2-x1)/5.; //on the border more iteration
509 Double_t ndx=dx;
510
511 if (x2>(x1+kprec)) {
512 for (Double_t x = x1; x<x2+kprec ;){
513 //new step SIZE
514 nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x-xch)*(x-xch)/(2*fOrigSigmaX*fOrigSigmaX))),3);
515 ndx = fOrigSigmaX/Double_t(nx);
516 if (ndx>(x2-x-dx)) {
517 ndx =x2-x-dx;
518 }
519 if ( ( (x+dx+ndx)<TMath::Max(xp3,xp1)) || ( (x+dx+ndx)>TMath::Min(xp4,xp2))) {
520 ndx/=5.;
521 }
522 if (ndx<kprec) ndx=2*kprec;
523 //INTEGRAL APROXIMATION
524 Double_t ddx,ddy,dddx,dddy;
525 ddx = xch-(x+dx/2.);
526 ddy = fCurrentY-(y+dy/2.);
527 dddx = cos*ddx-sin*ddy;
528 dddy = sin*ddx+cos*ddy;
529 Double_t z0=fGRF->Eval(dddx,dddy); //middle point
530
531 ddx = xch-(x+dx/2.);
532 ddy = fCurrentY-(y);
533 dddx = cos*ddx-sin*ddy;
534 dddy = sin*ddx+cos*ddy;
535 Double_t z1=fGRF->Eval(dddx,dddy); //point down
536
537 ddx = xch-(x+dx/2.);
538 ddy = fCurrentY-(y+dy);
539 dddx = cos*ddx-sin*ddy;
540 dddy = sin*ddx+cos*ddy;
541 Double_t z3=fGRF->Eval(dddx,dddy); //point up
542
543 ddx = xch-(x);
544 ddy = fCurrentY-(y+dy/2.);
545 dddx = cos*ddx-sin*ddy;
546 dddy = sin*ddx+cos*ddy;
547 Double_t z2=fGRF->Eval(dddx,dddy); //point left
548
549 ddx = xch-(x+dx);
550 ddy = fCurrentY-(y+dy/2.);
551 dddx = cos*ddx-sin*ddy;
552 dddy = sin*ddx+cos*ddy;
553 Double_t z4=fGRF->Eval(dddx,dddy); //point right
554
555
556 if (z0<0) {z0=0;z1=0;z2=0;z3=0;z4=0;}
557
558 Double_t f2x= (z3+z1-2*z0)*4.;//second derivation in y
559 Double_t f2y= (z2+z4-2*z0)*4.;//second derivation in x
560 Double_t f1y= (z3-z1);
561 Double_t z ;
562 z = (z0+f2x/6.+f2y/6.);//second order aproxiation of integral
563 if (kx>kprec){ //positive derivation
564 if (x<(xp1+dy*kx)){ //calculate volume at left border
565 Double_t xx1 = x;
566 Double_t xx2 = TMath::Min(x+dx,xp1+dy*kx);
567 Double_t yy1 = y+(xx1-xp1)/kx;
568 Double_t yy2 = TMath::Min(y+(xx2-xp1)/kx,y+dy);
569 z=z0;
570 if (yy2<y+dy) {
571 z-= z0*(y+dy-yy2)/dy; //constant part rectangle
572 z-= f1y*(xx2-xx1)*(y+dy-yy2)*(y+dy-yy2)/(2.*dx*dy);
573 }
574 z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part rectangle
575
576 }
577 if (x>xp2){ //calculate volume at right border
578 Double_t xx1 = x;
579 Double_t xx2 = x+dx;
580 Double_t yy1 = y+(xx1-xp2)/kx;
581 Double_t yy2 = y+(xx2-xp2)/kx;
582 z=z0;
583 //rectangle part
584 z-=z0*(yy1-y)/dy; //constant part
585 z-=f1y*(xx2-xx1)*(yy1-y)*(yy1-y)/(2*dx*dy);
586 //triangle part
587 z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part
588 }
589 }
590 if (kx<-kprec){ //negative derivation
591 if (x<(xp1+dy*kx)){ //calculate volume at left border
592 Double_t xx1 = x;
593 Double_t xx2 = TMath::Min(x+dx,xp3-dy/kx);
594 Double_t yy1 = y+(xx1-xp1)/kx;
595 Double_t yy2 = TMath::Max(y,yy1+(xx2-xx1)/kx); //yy2<yy1
596 z = z0;
597 z-= z0*(yy2-y)/dy; // constant part rectangle
598 z-= f1y*(xx2-xx1)*(yy2-y)*(yy2-y)/(2.*dx*dy);
599 z-=z0*(xx2-xx1)*(yy1-yy2)/(2*dx*dy); //constant part triangle
600 }
601 if (x>xp2){ //calculate volume at right border
602 Double_t xx1 = TMath::Max(x,xp2+dy*kx);
603 Double_t xx2 = x+dx;
604 Double_t yy1 = TMath::Min(y+dy,y-(xp2-xx1)/kx);
605 Double_t yy2 = y-(xp2-xx2)/kx;
606 z=z0;
607 z-=z0*(yy2-y)/dy; //constant part rextangle
608 z-= f1y*(xx2-xx1)*(yy2-y)*(yy2-y)/(2.*dx*dy);
609 z-=z0*(xx2-xx1)*(yy1-yy2)/(2*dx*dy); //constant part triangle
610 }
611 }
612
613 if (z>0.) sumch+=fKNorm*z*dx*dy/fInteg;
614
615 x+=dx;
616 dx = ndx;
617 }; //loop over x
618 fcharge[i]+=sumch;
619 }//if x2>x1
620 y+=dy;
621 dy =ndy;
622 }//step over different y
623 k*=-1.;
624 }//step over chevron
625
626 }//step over different points on line NPRF
627}
628
629void AliTPCPRF2D::UpdateSigma()
630{
631 //
632 //calulate effective sigma X and sigma y of PRF
633 fMeanX = 0;
634 fMeanY = 0;
635 fSigmaX = 0;
636 fSigmaY = 0;
637
638 Float_t sum =0;
639 Int_t i;
640 Float_t x,y;
641
642 for (i=-1; i<=fNYdiv; i++){
643 if (fNYdiv == 1) y = fY1;
644 else
645 y = fY1+Float_t(i)*(fY2-fY1)/Float_t(fNYdiv-1);
646 for (x =-fNPRF*fDStep; x<fNPRF*fDStep;x+=fDStep)
647 {
648 //x in cm fWidth in cm
649 Float_t weight = GetPRF(x,y);
650 fSigmaX+=x*x*weight;
651 fSigmaY+=y*y*weight;
652 fMeanX+=x*weight;
653 fMeanY+=y*weight;
654 sum+=weight;
655 };
656 }
657 if (sum>0){
658 fMeanX/=sum;
659 fMeanY/=sum;
660 fSigmaX = TMath::Sqrt(fSigmaX/sum-fMeanX*fMeanX);
661 fSigmaY = TMath::Sqrt(fSigmaY/sum-fMeanY*fMeanY);
662 }
663 else fSigmaX=0;
664}
665
666
667void AliTPCPRF2D::Streamer(TBuffer &R__b)
668{
669 // Stream an object of class AliTPCPRF2D
670
671 if (R__b.IsReading()) {
672 UInt_t R__s, R__c;
673 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
674 AliTPCPRF2D::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
675 //read functions
676 if (strncmp(fType,"User",3)!=0){
677 delete fGRF;
678 if (strncmp(fType,"Gauss",3)==0)
679 fGRF = new TF2("funGauss2D",funGauss2D,-5.,5.,-5.,5.,4);
680 if (strncmp(fType,"Cosh",3)==0)
681 fGRF = new TF2("funCosh2D",funCosh2D,-5.,5.,-5.,5.,4);
682 if (strncmp(fType,"Gati",3)==0)
683 fGRF = new TF2("funGati2D",funGati2D,-5.,5.,-5.,5.,5);
684 if (fGRF!=0) fGRF->SetParameters(funParam);
685 }
686 //calculate conversion coefitient to convert position to virtual wire
687 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
688 fDStepM1=1/fDStep;
689 } else {
690 AliTPCPRF2D::Class()->WriteBuffer(R__b,this);
691 }
692}
693
694
695TH1F * AliTPCPRF2D::GenerDrawXHisto(Float_t x1, Float_t x2,Float_t y)
696{
697 //gener one dimensional hist of pad response function
698 // at position y
699 char s[100];
700 const Int_t kn=200;
701 sprintf(s,"Pad Response Function");
702 TH1F * hPRFc = new TH1F("hPRFc",s,kn+1,x1,x2);
703 Float_t x=x1;
704 Float_t y1;
705
706 for (Int_t i = 0;i<kn+1;i++)
707 {
708 x+=(x2-x1)/Float_t(kn);
709 y1 = GetPRF(x,y);
710 hPRFc->Fill(x,y1);
711 };
712 hPRFc->SetXTitle("pad (cm)");
713 return hPRFc;
714}
715
716AliH2F * AliTPCPRF2D::GenerDrawHisto(Float_t x1, Float_t x2, Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
717{
718 //
719 //gener two dimensional histogram with PRF
720 //
721 char s[100];
722 sprintf(s,"Pad Response Function");
723 AliH2F * hPRFc = new AliH2F("hPRFc",s,Nx,x1,x2,Ny,y1,y2);
724 Float_t dx=(x2-x1)/Float_t(Nx);
725 Float_t dy=(y2-y1)/Float_t(Ny) ;
726 Float_t x,y,z;
727 x = x1;
728 y = y1;
729 for ( Int_t i = 0;i<=Nx;i++,x+=dx){
730 y=y1;
731 for (Int_t j = 0;j<=Ny;j++,y+=dy){
732 z = GetPRF(x,y);
733 hPRFc->SetCellContent(i,j,z);
734 };
735 };
736 hPRFc->SetXTitle("pad direction (cm)");
737 hPRFc->SetYTitle("pad row direction (cm)");
738 hPRFc->SetTitleOffset(1.5,"X");
739 hPRFc->SetTitleOffset(1.5,"Y");
740 return hPRFc;
741}
742
743
744AliH2F * AliTPCPRF2D::GenerDrawDistHisto(Float_t x1, Float_t x2, Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
745{
746 //return histogram with distortion
747 const Float_t kminth=0.00001;
748 if (thr<kminth) thr=kminth;
749 char s[100];
750 sprintf(s,"COG distortion of PRF (threshold=%2.2f)",thr);
751 AliH2F * hPRFDist = new AliH2F("hDistortion",s,Nx,x1,x2,Ny,y1,y2);
752 Float_t dx=(x2-x1)/Float_t(Nx);
753 Float_t dy=(y2-y1)/Float_t(Ny) ;
754 Float_t x,y,z,ddx;
755 x=x1;
756 for ( Int_t i = 0;i<=Nx;i++,x+=dx){
757 y=y1;
758 for(Int_t j = 0;j<=Ny;j++,y+=dy)
759 {
760 Float_t sumx=0;
761 Float_t sum=0;
762 for (Int_t k=-3;k<=3;k++)
763 {
764 Float_t padx=Float_t(k)*fWidth;
765 z = GetPRF(x-padx,y);
766 if (z>thr){
767 sum+=z;
768 sumx+=z*padx;
769 }
770 };
771 if (sum>kminth)
772 {
773 ddx = (x-(sumx/sum));
774 }
775 else ddx=-1;
776 if (TMath::Abs(ddx)<10) hPRFDist->SetCellContent(i,j,ddx);
777 }
778 }
779
780 hPRFDist->SetXTitle("pad direction (cm)");
781 hPRFDist->SetYTitle("pad row direction (cm)");
782 hPRFDist->SetTitleOffset(1.5,"X");
783 hPRFDist->SetTitleOffset(1.5,"Y");
784 return hPRFDist;
785}
786
787
788
789
790
791void AliTPCPRF2D::DrawX(Float_t x1 ,Float_t x2,Float_t y1,Float_t y2, Int_t N)
792{
793 //
794 //draw pad response function at interval <x1,x2> at given y position
795 //
796 if (N<0) return;
797 TCanvas * c1 = new TCanvas("PRFX","Pad response function",700,900);
798 c1->cd();
799
800 TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
801 comment->SetTextAlign(12);
802 comment->SetFillColor(42);
803 DrawComment(comment);
804 comment->Draw();
805 c1->cd();
806
807 TPad * pad2 = new TPad("pPRF","",0.05,0.22,0.95,0.95);
808 pad2->Divide(2,(N+1)/2);
809 pad2->Draw();
810 gStyle->SetOptFit(1);
811 gStyle->SetOptStat(1);
812 for (Int_t i=0;i<N;i++){
813 char ch[200];
814 Float_t y;
815 if (N==1) y=y1;
816 else y = y1+i*(y2-y1)/Float_t(N-1);
817 pad2->cd(i+1);
818 TH1F * hPRFc =GenerDrawXHisto(x1, x2,y);
819 sprintf(ch,"PRF at wire position: %2.3f",y);
820 hPRFc->SetTitle(ch);
821 sprintf(ch,"PRF %d",i);
822 hPRFc->SetName(ch);
823 hPRFc->Fit("gaus");
824 }
825
826}
827
828
829
830void AliTPCPRF2D::DrawPRF(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
831{
832 //
833 //
834 TCanvas * c1 = new TCanvas("canPRF","Pad response function",700,900);
835 c1->cd();
836 TPad * pad2 = new TPad("pad2PRF","",0.05,0.22,0.95,0.95);
837 pad2->Draw();
838 gStyle->SetOptFit(1);
839 gStyle->SetOptStat(1);
840 TH2F * hPRFc = GenerDrawHisto(x1, x2, y1, y2, Nx,Ny);
841 pad2->cd();
842 hPRFc->Draw("surf");
843 c1->cd();
844 TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
845 comment->SetTextAlign(12);
846 comment->SetFillColor(42);
847 DrawComment(comment);
848 comment->Draw();
849}
850
851void AliTPCPRF2D::DrawDist(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
852{
853 //
854 //draw distortion of the COG method - for different threshold parameter
855 TCanvas * c1 = new TCanvas("padDistortion","COG distortion",700,900);
856 c1->cd();
857 TPad * pad1 = new TPad("dist","",0.05,0.55,0.95,0.95,21);
858 pad1->Draw();
859 TPad * pad2 = new TPad("dist","",0.05,0.22,0.95,0.53,21);
860 pad2->Draw();
861 gStyle->SetOptFit(1);
862 gStyle->SetOptStat(0);
863
864 AliH2F * hPRFDist = GenerDrawDistHisto(x1, x2, y1, y2, Nx,Ny,thr);
865
866 pad1->cd();
867 hPRFDist->Draw("surf");
868 Float_t distmax =hPRFDist->GetMaximum();
869 Float_t distmin =hPRFDist->GetMinimum();
870 gStyle->SetOptStat(1);
871
872 TH1F * dist = hPRFDist->GetAmplitudes(distmin,distmax,distmin-1);
873 pad2->cd();
874 dist->Draw();
875 c1->cd();
876 TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
877 comment->SetTextAlign(12);
878 comment->SetFillColor(42);
879 DrawComment(comment);
880 comment->Draw();
881}
882
883void AliTPCPRF2D::DrawComment(TPaveText *comment)
884{
885 //
886 //function to write comment to picture
887
888 char s[100];
889 //draw comments to picture
890 TText * title = comment->AddText("Pad Response Function parameters:");
891 title->SetTextSize(0.03);
892 sprintf(s,"Height of pad: %2.2f cm",fHeightFull);
893 comment->AddText(s);
894 sprintf(s,"Width pad: %2.2f cm",fWidth);
895 comment->AddText(s);
896 sprintf(s,"Pad Angle: %2.2f ",fPadAngle);
897 comment->AddText(s);
898
899 if (TMath::Abs(fK)>0.0001){
900 sprintf(s,"Height of one chevron unit h: %2.2f cm",2*fHeightS);
901 comment->AddText(s);
902 sprintf(s,"Overlap factor: %2.2f",fK);
903 comment->AddText(s);
904 }
905
906 if (strncmp(fType,"User",3)==0){
907 sprintf(s,"Charge distribution - user defined function %s ",fGRF->GetTitle());
908 comment->AddText(s);
909 sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
910 comment->AddText(s);
911 sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
912 comment->AddText(s);
913 }
914 if (strncmp(fType,"Gauss",3)==0){
915 sprintf(s,"Gauss charge distribution");
916 comment->AddText(s);
917 sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
918 comment->AddText(s);
919 sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
920 comment->AddText(s);
921 }
922 if (strncmp(fType,"Gati",3)==0){
923 sprintf(s,"Gati charge distribution");
924 comment->AddText(s);
925 sprintf(s,"K3X of Gati : %2.2f ",fK3X);
926 comment->AddText(s);
927 sprintf(s,"K3Y of Gati: %2.2f ",fK3Y);
928 comment->AddText(s);
929 sprintf(s,"Wire to Pad Distance: %2.2f ",fPadDistance);
930 comment->AddText(s);
931 }
932 if (strncmp(fType,"Cosh",3)==0){
933 sprintf(s,"Cosh charge distribution");
934 comment->AddText(s);
935 sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
936 comment->AddText(s);
937 sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
938 comment->AddText(s);
939 }
940 sprintf(s,"Normalisation: %2.2f ",fKNorm);
941 comment->AddText(s);
942}
943