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