]> 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
Adding new Task - Krypton calibration (Jacek)
[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
55static Double_t funGauss2D(Double_t *x, Double_t * par)
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
63static Double_t funCosh2D(Double_t *x, Double_t * par)
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
70static Double_t funGati2D(Double_t *x, Double_t * par)
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 ;
179c6296 137
8c555625 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);
6e7b5431 142 SetInterpolationType(2,0);
8c555625 143}
144
145AliTPCPRF2D::~AliTPCPRF2D()
146{
6e7b5431 147 if (fChargeArray!=0) delete [] fChargeArray;
148 if (fGRF !=0 ) fGRF->Delete();
8c555625 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;
8c555625 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;
6e7b5431 174 fK=fac;
8c555625 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
6e7b5431 185Float_t AliTPCPRF2D::GetPRF(Float_t xin, Float_t yin)
8c555625 186{
cc80f89e 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
6e7b5431 190 if (fChargeArray==0) return 0;
8c555625 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
6e7b5431 195 if (fInterY ==0){
8c555625 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;
6e7b5431 199 fcharge = &(fChargeArray[i*fNPRF]);
8c555625 200 return GetPRFActiv(xin);
201 }
0798b21e 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;
8c555625 235}
236
237
238Float_t AliTPCPRF2D::GetPRFActiv(Float_t xin)
239{
cc80f89e 240 //GEt response function on given charege line
241 //return spline aproximaton
8c555625 242 Float_t x = (xin*fDStepM1)+fNPRF/2;
243 Int_t i = Int_t(x);
244
6e7b5431 245 if ( (i>1) && ((i+2)<fNPRF)) {
73042f01 246 Float_t a,b,c,d,k,l;
8c555625 247 a = fcharge[i];
248 b = (fcharge[i+1]-fcharge[i-1])*0.5;
73042f01 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;
8c555625 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{
cc80f89e 263 //function which returnoriginal charge distribution
264 //this function is just normalised for fKnorm
6e7b5431 265 if (GetGRF() != 0 )
266 return fKNorm*GetGRF()->Eval(xin,yin)/fInteg;
8c555625 267 else
268 return 0.;
269}
270
271
272void AliTPCPRF2D::SetParam( TF2 * GRF, Float_t kNorm,
273 Float_t sigmaX, Float_t sigmaY)
274{
cc80f89e 275 //adjust parameters of the original charge distribution
276 //and pad size parameters
8c555625 277 if (fGRF !=0 ) fGRF->Delete();
278 fGRF = GRF;
6e7b5431 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;
cc80f89e 283 fOrigSigmaX=sigmaX;
284 fOrigSigmaY=sigmaY;
6e7b5431 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
8c555625 298}
299
300
301void AliTPCPRF2D::SetGauss(Float_t sigmaX, Float_t sigmaY,
302 Float_t kNorm)
303{
cc80f89e 304 //
305 // set parameters for Gauss generic charge distribution
306 //
6e7b5431 307 fKNorm = kNorm;
308 fOrigSigmaX=sigmaX;
309 fOrigSigmaY=sigmaY;
310 sprintf(fType,"Gauss");
8c555625 311 if (fGRF !=0 ) fGRF->Delete();
2a336e15 312 fGRF = new TF2("funGauss2D",funGauss2D,-5.,5.,-5.,5.,4);
6e7b5431 313
8c555625 314 funParam[0]=sigmaX;
315 funParam[1]=sigmaY;
316 funParam[2]=fK;
317 funParam[3]=fHeightS;
6e7b5431 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
8c555625 334}
8c555625 335void AliTPCPRF2D::SetCosh(Float_t sigmaX, Float_t sigmaY,
336 Float_t kNorm)
cc80f89e 337{
338 // set parameters for Cosh generic charge distribution
339 //
6e7b5431 340 fKNorm = kNorm;
341 fOrigSigmaX=sigmaX;
342 fOrigSigmaY=sigmaY;
343 sprintf(fType,"Cosh");
8c555625 344 if (fGRF !=0 ) fGRF->Delete();
2a336e15 345 fGRF = new TF2("funCosh2D", funCosh2D,-5.,5.,-5.,5.,4);
8c555625 346 funParam[0]=sigmaX;
347 funParam[1]=sigmaY;
348 funParam[2]=fK;
349 funParam[3]=fHeightS;
350 fGRF->SetParameters(funParam);
6e7b5431 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
8c555625 362}
363
364void AliTPCPRF2D::SetGati(Float_t K3X, Float_t K3Y,
365 Float_t padDistance,
366 Float_t kNorm)
367{
cc80f89e 368 // set parameters for Gati generic charge distribution
369 //
6e7b5431 370 fKNorm = kNorm;
8c555625 371 fK3X=K3X;
372 fK3Y=K3Y;
6e7b5431 373 fPadDistance=padDistance;
374 sprintf(fType,"Gati");
375 if (fGRF !=0 ) fGRF->Delete();
2a336e15 376 fGRF = new TF2("funGati2D", funGati2D,-5.,5.,-5.,5.,5);
6e7b5431 377
8c555625 378 funParam[0]=padDistance;
379 funParam[1]=K3X;
380 funParam[2]=fK;
381 funParam[3]=fHeightS;
382 funParam[4]=K3Y;
383 fGRF->SetParameters(funParam);
cc80f89e 384 fOrigSigmaX=padDistance;
385 fOrigSigmaY=padDistance;
6e7b5431 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 };
8c555625 396}
397
398
399
400void AliTPCPRF2D::Update()
401{
cc80f89e 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;
6e7b5431 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;
cc80f89e 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;
8c555625 436 else
cc80f89e 437 fCurrentY = fY1+Double_t(i)*(fY2-fY1)/Double_t(fNYdiv-1);
6e7b5431 438 fcharge = &(fChargeArray[i*fNPRF]);
8c555625 439 Update1();
440 }
cc80f89e 441 //calculate conversion coefitient to convert position to virtual wire
442 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
443 fDStepM1=1/fDStep;
444 UpdateSigma();
8c555625 445}
446
8c555625 447void AliTPCPRF2D::Update1()
448{
cc80f89e 449 //
450 //update fields with interpolated values for
451 //PRF calculation for given charge line
8c555625 452 Int_t i;
cc80f89e 453 Double_t cos = TMath::Cos(fChargeAngle);
454 Double_t sin = TMath::Sin(fChargeAngle);
6e7b5431 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;
cc80f89e 463
6e7b5431 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
cc80f89e 483
6e7b5431 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; //..
8c555625 500
6e7b5431 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;
8c555625 518 }
6e7b5431 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
cc80f89e 625
6e7b5431 626 }//step over different points on line NPRF
cc80f89e 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
8c555625 638 Float_t sum =0;
cc80f89e 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;
8c555625 655 };
cc80f89e 656 }
8c555625 657 if (sum>0){
cc80f89e 658 fMeanX/=sum;
659 fMeanY/=sum;
660 fSigmaX = TMath::Sqrt(fSigmaX/sum-fMeanX*fMeanX);
661 fSigmaY = TMath::Sqrt(fSigmaY/sum-fMeanY*fMeanY);
8c555625 662 }
663 else fSigmaX=0;
8c555625 664}
665
cc80f89e 666
8c555625 667void AliTPCPRF2D::Streamer(TBuffer &R__b)
668{
669 // Stream an object of class AliTPCPRF2D
670
671 if (R__b.IsReading()) {
a8a6107b 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);
8c555625 675 //read functions
a8a6107b 676 if (strncmp(fType,"User",3)!=0){
677 delete fGRF;
678 if (strncmp(fType,"Gauss",3)==0)
2a336e15 679 fGRF = new TF2("funGauss2D",funGauss2D,-5.,5.,-5.,5.,4);
a8a6107b 680 if (strncmp(fType,"Cosh",3)==0)
2a336e15 681 fGRF = new TF2("funCosh2D",funCosh2D,-5.,5.,-5.,5.,4);
a8a6107b 682 if (strncmp(fType,"Gati",3)==0)
2a336e15 683 fGRF = new TF2("funGati2D",funGati2D,-5.,5.,-5.,5.,5);
a8a6107b 684 if (fGRF!=0) fGRF->SetParameters(funParam);
6e7b5431 685 }
8c555625 686 //calculate conversion coefitient to convert position to virtual wire
687 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
688 fDStepM1=1/fDStep;
689 } else {
a8a6107b 690 AliTPCPRF2D::Class()->WriteBuffer(R__b,this);
8c555625 691 }
692}
693
694
6e7b5431 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);
8c555625 703 Float_t x=x1;
704 Float_t y1;
8c555625 705
6e7b5431 706 for (Int_t i = 0;i<kn+1;i++)
8c555625 707 {
6e7b5431 708 x+=(x2-x1)/Float_t(kn);
709 y1 = GetPRF(x,y);
8c555625 710 hPRFc->Fill(x,y1);
711 };
6e7b5431 712 hPRFc->SetXTitle("pad (cm)");
713 return hPRFc;
714}
8c555625 715
6e7b5431 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
8c555625 800 TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
801 comment->SetTextAlign(12);
802 comment->SetFillColor(42);
6e7b5431 803 DrawComment(comment);
8c555625 804 comment->Draw();
6e7b5431 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
8c555625 826}
827
828
829
6e7b5431 830void AliTPCPRF2D::DrawPRF(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
8c555625 831{
6e7b5431 832 //
833 //
8c555625 834 TCanvas * c1 = new TCanvas("canPRF","Pad response function",700,900);
835 c1->cd();
6e7b5431 836 TPad * pad2 = new TPad("pad2PRF","",0.05,0.22,0.95,0.95);
837 pad2->Draw();
8c555625 838 gStyle->SetOptFit(1);
8c555625 839 gStyle->SetOptStat(1);
6e7b5431 840 TH2F * hPRFc = GenerDrawHisto(x1, x2, y1, y2, Nx,Ny);
8c555625 841 pad2->cd();
6e7b5431 842 hPRFc->Draw("surf");
8c555625 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);
6e7b5431 847 DrawComment(comment);
8c555625 848 comment->Draw();
849}
850
6e7b5431 851void AliTPCPRF2D::DrawDist(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
8c555625 852{
6e7b5431 853 //
854 //draw distortion of the COG method - for different threshold parameter
8c555625 855 TCanvas * c1 = new TCanvas("padDistortion","COG distortion",700,900);
856 c1->cd();
6e7b5431 857 TPad * pad1 = new TPad("dist","",0.05,0.55,0.95,0.95,21);
8c555625 858 pad1->Draw();
6e7b5431 859 TPad * pad2 = new TPad("dist","",0.05,0.22,0.95,0.53,21);
8c555625 860 pad2->Draw();
8c555625 861 gStyle->SetOptFit(1);
862 gStyle->SetOptStat(0);
6e7b5431 863
864 AliH2F * hPRFDist = GenerDrawDistHisto(x1, x2, y1, y2, Nx,Ny,thr);
865
8c555625 866 pad1->cd();
6e7b5431 867 hPRFDist->Draw("surf");
868 Float_t distmax =hPRFDist->GetMaximum();
869 Float_t distmin =hPRFDist->GetMinimum();
870 gStyle->SetOptStat(1);
8c555625 871
6e7b5431 872 TH1F * dist = hPRFDist->GetAmplitudes(distmin,distmax,distmin-1);
873 pad2->cd();
874 dist->Draw();
8c555625 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);
6e7b5431 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:");
8c555625 891 title->SetTextSize(0.03);
6e7b5431 892 sprintf(s,"Height of pad: %2.2f cm",fHeightFull);
8c555625 893 comment->AddText(s);
6e7b5431 894 sprintf(s,"Width pad: %2.2f cm",fWidth);
8c555625 895 comment->AddText(s);
6e7b5431 896 sprintf(s,"Pad Angle: %2.2f ",fPadAngle);
8c555625 897 comment->AddText(s);
8c555625 898
6e7b5431 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);
8c555625 942}
943