]> 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
Additional protection
[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;
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;
cc80f89e 288 fOrigSigmaX=sigmaX;
289 fOrigSigmaY=sigmaY;
6e7b5431 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
8c555625 303}
304
305
306void AliTPCPRF2D::SetGauss(Float_t sigmaX, Float_t sigmaY,
307 Float_t kNorm)
308{
cc80f89e 309 //
310 // set parameters for Gauss generic charge distribution
311 //
6e7b5431 312 fKNorm = kNorm;
313 fOrigSigmaX=sigmaX;
314 fOrigSigmaY=sigmaY;
315 sprintf(fType,"Gauss");
8c555625 316 if (fGRF !=0 ) fGRF->Delete();
798017c7 317 fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
6e7b5431 318
8c555625 319 funParam[0]=sigmaX;
320 funParam[1]=sigmaY;
321 funParam[2]=fK;
322 funParam[3]=fHeightS;
6e7b5431 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
8c555625 339}
8c555625 340void AliTPCPRF2D::SetCosh(Float_t sigmaX, Float_t sigmaY,
341 Float_t kNorm)
cc80f89e 342{
343 // set parameters for Cosh generic charge distribution
344 //
6e7b5431 345 fKNorm = kNorm;
346 fOrigSigmaX=sigmaX;
347 fOrigSigmaY=sigmaY;
348 sprintf(fType,"Cosh");
8c555625 349 if (fGRF !=0 ) fGRF->Delete();
798017c7 350 fGRF = new TF2("FunCosh2D", FunCosh2D,-5.,5.,-5.,5.,4);
8c555625 351 funParam[0]=sigmaX;
352 funParam[1]=sigmaY;
353 funParam[2]=fK;
354 funParam[3]=fHeightS;
355 fGRF->SetParameters(funParam);
6e7b5431 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
8c555625 367}
368
369void AliTPCPRF2D::SetGati(Float_t K3X, Float_t K3Y,
370 Float_t padDistance,
371 Float_t kNorm)
372{
cc80f89e 373 // set parameters for Gati generic charge distribution
374 //
6e7b5431 375 fKNorm = kNorm;
8c555625 376 fK3X=K3X;
377 fK3Y=K3Y;
6e7b5431 378 fPadDistance=padDistance;
379 sprintf(fType,"Gati");
380 if (fGRF !=0 ) fGRF->Delete();
798017c7 381 fGRF = new TF2("FunGati2D", FunGati2D,-5.,5.,-5.,5.,5);
6e7b5431 382
8c555625 383 funParam[0]=padDistance;
384 funParam[1]=K3X;
385 funParam[2]=fK;
386 funParam[3]=fHeightS;
387 funParam[4]=K3Y;
388 fGRF->SetParameters(funParam);
cc80f89e 389 fOrigSigmaX=padDistance;
390 fOrigSigmaY=padDistance;
6e7b5431 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 };
8c555625 401}
402
403
404
405void AliTPCPRF2D::Update()
406{
cc80f89e 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;
6e7b5431 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;
cc80f89e 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;
8c555625 441 else
cc80f89e 442 fCurrentY = fY1+Double_t(i)*(fY2-fY1)/Double_t(fNYdiv-1);
6e7b5431 443 fcharge = &(fChargeArray[i*fNPRF]);
8c555625 444 Update1();
445 }
cc80f89e 446 //calculate conversion coefitient to convert position to virtual wire
447 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
448 fDStepM1=1/fDStep;
449 UpdateSigma();
8c555625 450}
451
8c555625 452void AliTPCPRF2D::Update1()
453{
cc80f89e 454 //
455 //update fields with interpolated values for
456 //PRF calculation for given charge line
8c555625 457 Int_t i;
cc80f89e 458 Double_t cos = TMath::Cos(fChargeAngle);
459 Double_t sin = TMath::Sin(fChargeAngle);
6e7b5431 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;
cc80f89e 468
6e7b5431 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
cc80f89e 488
6e7b5431 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; //..
8c555625 505
6e7b5431 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;
8c555625 523 }
6e7b5431 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
cc80f89e 630
6e7b5431 631 }//step over different points on line NPRF
cc80f89e 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
8c555625 643 Float_t sum =0;
cc80f89e 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;
8c555625 660 };
cc80f89e 661 }
8c555625 662 if (sum>0){
cc80f89e 663 fMeanX/=sum;
664 fMeanY/=sum;
665 fSigmaX = TMath::Sqrt(fSigmaX/sum-fMeanX*fMeanX);
666 fSigmaY = TMath::Sqrt(fSigmaY/sum-fMeanY*fMeanY);
8c555625 667 }
668 else fSigmaX=0;
8c555625 669}
670
cc80f89e 671
798017c7 672void AliTPCPRF2D::Streamer(TBuffer &xRuub)
8c555625 673{
674 // Stream an object of class AliTPCPRF2D
675
798017c7 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);
8c555625 680 //read functions
a8a6107b 681 if (strncmp(fType,"User",3)!=0){
682 delete fGRF;
683 if (strncmp(fType,"Gauss",3)==0)
798017c7 684 fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
a8a6107b 685 if (strncmp(fType,"Cosh",3)==0)
798017c7 686 fGRF = new TF2("FunCosh2D",FunCosh2D,-5.,5.,-5.,5.,4);
a8a6107b 687 if (strncmp(fType,"Gati",3)==0)
798017c7 688 fGRF = new TF2("FunGati2D",FunGati2D,-5.,5.,-5.,5.,5);
a8a6107b 689 if (fGRF!=0) fGRF->SetParameters(funParam);
6e7b5431 690 }
8c555625 691 //calculate conversion coefitient to convert position to virtual wire
692 fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
693 fDStepM1=1/fDStep;
694 } else {
798017c7 695 AliTPCPRF2D::Class()->WriteBuffer(xRuub,this);
8c555625 696 }
697}
698
699
6e7b5431 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);
8c555625 708 Float_t x=x1;
709 Float_t y1;
8c555625 710
6e7b5431 711 for (Int_t i = 0;i<kn+1;i++)
8c555625 712 {
6e7b5431 713 x+=(x2-x1)/Float_t(kn);
714 y1 = GetPRF(x,y);
8c555625 715 hPRFc->Fill(x,y1);
716 };
6e7b5431 717 hPRFc->SetXTitle("pad (cm)");
718 return hPRFc;
719}
8c555625 720
6e7b5431 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
8c555625 805 TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
806 comment->SetTextAlign(12);
807 comment->SetFillColor(42);
6e7b5431 808 DrawComment(comment);
8c555625 809 comment->Draw();
6e7b5431 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
8c555625 831}
832
833
834
6e7b5431 835void AliTPCPRF2D::DrawPRF(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
8c555625 836{
6e7b5431 837 //
838 //
8c555625 839 TCanvas * c1 = new TCanvas("canPRF","Pad response function",700,900);
840 c1->cd();
6e7b5431 841 TPad * pad2 = new TPad("pad2PRF","",0.05,0.22,0.95,0.95);
842 pad2->Draw();
8c555625 843 gStyle->SetOptFit(1);
8c555625 844 gStyle->SetOptStat(1);
6e7b5431 845 TH2F * hPRFc = GenerDrawHisto(x1, x2, y1, y2, Nx,Ny);
8c555625 846 pad2->cd();
6e7b5431 847 hPRFc->Draw("surf");
8c555625 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);
6e7b5431 852 DrawComment(comment);
8c555625 853 comment->Draw();
854}
855
6e7b5431 856void AliTPCPRF2D::DrawDist(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
8c555625 857{
6e7b5431 858 //
859 //draw distortion of the COG method - for different threshold parameter
8c555625 860 TCanvas * c1 = new TCanvas("padDistortion","COG distortion",700,900);
861 c1->cd();
6e7b5431 862 TPad * pad1 = new TPad("dist","",0.05,0.55,0.95,0.95,21);
8c555625 863 pad1->Draw();
6e7b5431 864 TPad * pad2 = new TPad("dist","",0.05,0.22,0.95,0.53,21);
8c555625 865 pad2->Draw();
8c555625 866 gStyle->SetOptFit(1);
867 gStyle->SetOptStat(0);
6e7b5431 868
869 AliH2F * hPRFDist = GenerDrawDistHisto(x1, x2, y1, y2, Nx,Ny,thr);
870
8c555625 871 pad1->cd();
6e7b5431 872 hPRFDist->Draw("surf");
873 Float_t distmax =hPRFDist->GetMaximum();
874 Float_t distmin =hPRFDist->GetMinimum();
875 gStyle->SetOptStat(1);
8c555625 876
6e7b5431 877 TH1F * dist = hPRFDist->GetAmplitudes(distmin,distmax,distmin-1);
878 pad2->cd();
879 dist->Draw();
8c555625 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);
6e7b5431 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:");
8c555625 896 title->SetTextSize(0.03);
6e7b5431 897 sprintf(s,"Height of pad: %2.2f cm",fHeightFull);
8c555625 898 comment->AddText(s);
6e7b5431 899 sprintf(s,"Width pad: %2.2f cm",fWidth);
8c555625 900 comment->AddText(s);
6e7b5431 901 sprintf(s,"Pad Angle: %2.2f ",fPadAngle);
8c555625 902 comment->AddText(s);
8c555625 903
6e7b5431 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);
8c555625 947}
948