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