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