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