Patch for the tracker
[u/mrichter/AliRoot.git] / TPC / AliTPCPRF2D.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //  AliTPCPRF2D -                                                            //
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 ///////////////////////////////////////////////////////////////////////////////
31
32 #include <Riostream.h>
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>
42 #include <string.h>
43
44 #include "AliH2F.h"
45 #include "AliTPCPRF2D.h"
46
47
48 extern TStyle * gStyle;
49
50 const Double_t AliTPCPRF2D::fgkDegtoRad = 0.01745329251994;
51 const Double_t AliTPCPRF2D::fgkSQRT12=3.464101;
52 const Int_t   AliTPCPRF2D::fgkNPRF = 100;
53
54
55 static Double_t FunGauss2D(const Double_t *const x, const Double_t *const par)
56
57 //Gauss function  -needde by the generic function object 
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
63 static Double_t FunCosh2D(const Double_t *const x, const Double_t *const par)
64 {
65  //Cosh function  -needde by the generic function object 
66   return ( 1/(TMath::CosH(3.14159*x[0]/(2*par[0]))*
67            TMath::CosH(3.14159*x[1]/(2*par[1]))));
68 }    
69
70 static Double_t FunGati2D(const Double_t *const x, const Double_t *const par)
71 {
72   //Gati function  -needde by the generic function object 
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));
77   Float_t l=x[0]/par[0];
78   Float_t tan2=TMath::TanH(k2*l);
79   tan2*=tan2;
80   Float_t res = k1*(1-tan2)/(1+k3*tan2);
81  //par[4] = is equal to k3Y
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));
86   l=x[1]/par[0];
87   tan2=TMath::TanH(k2*l); 
88   tan2*=tan2;
89   res = res*k1*(1-tan2)/(1+k3*tan2);   
90   return res;  
91 }   
92
93 ///////////////////////////////////////////////////////////////////////////
94 ///////////////////////////////////////////////////////////////////////////
95
96 ClassImp(AliTPCPRF2D)
97
98 AliTPCPRF2D::AliTPCPRF2D()
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.)             
133 {
134   //default constructor for response function object
135
136   fNPRF =fgkNPRF ;
137   for(Int_t i=0;i<5;i++){
138         funParam[i]=0.;
139         fType[i]=0;
140   }
141   
142
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);
147   SetInterpolationType(2,0);
148 }
149
150 AliTPCPRF2D::~AliTPCPRF2D()
151 {
152   if (fChargeArray!=0) delete [] fChargeArray;
153   if (fGRF !=0 ) fGRF->Delete(); 
154 }
155
156 void 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;
162   fY1=y1;
163   fY2=y2;
164 }
165
166 void AliTPCPRF2D::SetPad(Float_t width, Float_t height)
167 {
168   //set base chevron parameters
169  fHeightFull=height;
170  fWidth=width;
171 }
172 void 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;
179   fK=fac;
180 }
181
182 void 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
190 Float_t AliTPCPRF2D::GetPRF(Float_t xin, Float_t yin)
191 {
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
195    if (fChargeArray==0) return 0;
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
200   if (fInterY ==0){   
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;
204     fcharge   = &(fChargeArray[i*fNPRF]);
205     return GetPRFActiv(xin);
206   }
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;            
240
241
242
243 Float_t AliTPCPRF2D::GetPRFActiv(Float_t xin)
244 {
245   //GEt response function on given charege line 
246   //return spline aproximaton 
247   Float_t x = (xin*fDStepM1)+fNPRF/2;
248   Int_t i = Int_t(x);
249   
250   if  ( (i>1) && ((i+2)<fNPRF)) {
251     Float_t a,b,c,d,k,l;
252     a = fcharge[i];
253     b = (fcharge[i+1]-fcharge[i-1])*0.5; 
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;
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
266 Float_t  AliTPCPRF2D::GetGRF(Float_t xin, Float_t yin)
267 {  
268   //function which returnoriginal charge distribution
269   //this function is just normalised for fKnorm
270   if (GetGRF() != 0 ) 
271     return fKNorm*GetGRF()->Eval(xin,yin)/fInteg;
272       else
273     return 0.;
274 }
275
276    
277 void AliTPCPRF2D::SetParam( TF2 *const GRF,  Float_t kNorm, 
278                        Float_t sigmaX, Float_t sigmaY)
279 {
280   //adjust parameters of the original charge distribution
281   //and pad size parameters
282    if (fGRF !=0 ) fGRF->Delete();
283    fGRF = GRF;
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;
288    fOrigSigmaX=sigmaX; 
289    fOrigSigmaY=sigmaY; 
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
303 }
304   
305
306 void AliTPCPRF2D::SetGauss(Float_t sigmaX, Float_t sigmaY,
307                       Float_t kNorm)
308 {
309   // 
310   // set parameters for Gauss generic charge distribution
311   //
312   fKNorm = kNorm;
313   fOrigSigmaX=sigmaX;
314   fOrigSigmaY=sigmaY;
315   sprintf(fType,"Gauss");
316   if (fGRF !=0 ) fGRF->Delete();
317   fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
318   
319   funParam[0]=sigmaX;
320   funParam[1]=sigmaY;  
321   funParam[2]=fK;
322   funParam[3]=fHeightS;    
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   
339 }
340 void AliTPCPRF2D::SetCosh(Float_t sigmaX, Float_t sigmaY,
341                      Float_t kNorm)
342
343   // set parameters for Cosh generic charge distribution
344   //
345   fKNorm = kNorm;
346   fOrigSigmaX=sigmaX;
347   fOrigSigmaY=sigmaY; 
348   sprintf(fType,"Cosh");
349   if (fGRF !=0 ) fGRF->Delete();
350   fGRF = new TF2("FunCosh2D",   FunCosh2D,-5.,5.,-5.,5.,4);   
351   funParam[0]=sigmaX;
352   funParam[1]=sigmaY;
353   funParam[2]=fK;  
354   funParam[3]=fHeightS;
355   fGRF->SetParameters(funParam);
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  
367 }
368
369 void AliTPCPRF2D::SetGati(Float_t K3X, Float_t K3Y,
370                      Float_t padDistance,
371                      Float_t kNorm)
372 {
373   // set parameters for Gati generic charge distribution
374   //
375   fKNorm = kNorm;
376   fK3X=K3X;
377   fK3Y=K3Y;
378   fPadDistance=padDistance;  
379   sprintf(fType,"Gati");
380   if (fGRF !=0 ) fGRF->Delete();
381   fGRF = new TF2("FunGati2D",   FunGati2D,-5.,5.,-5.,5.,5);  
382  
383   funParam[0]=padDistance;
384   funParam[1]=K3X;
385   funParam[2]=fK;  
386   funParam[3]=fHeightS;
387   funParam[4]=K3Y;
388   fGRF->SetParameters(funParam);
389   fOrigSigmaX=padDistance;
390   fOrigSigmaY=padDistance;
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   };
401 }
402
403
404
405 void AliTPCPRF2D::Update()
406 {
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;
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;
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;
441     else
442       fCurrentY = fY1+Double_t(i)*(fY2-fY1)/Double_t(fNYdiv-1);
443     fcharge   = &(fChargeArray[i*fNPRF]);
444     Update1();
445   }
446   //calculate conversion coefitient to convert position to virtual wire
447   fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
448   fDStepM1=1/fDStep;
449   UpdateSigma();
450 }
451
452 void AliTPCPRF2D::Update1()
453 {
454   //
455   //update fields  with interpolated values for
456   //PRF calculation for given charge line
457   Int_t i;
458   Double_t cos = TMath::Cos(fChargeAngle);
459   Double_t sin = TMath::Sin(fChargeAngle);
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;  
468     
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 
488         
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; //..  
505         
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;                  
523           }
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 
630     
631    }//step over different points on line NPRF
632 }
633
634 void 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  
643   Float_t sum =0;
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;
660     };  
661   }
662   if (sum>0){
663     fMeanX/=sum;
664     fMeanY/=sum;    
665     fSigmaX = TMath::Sqrt(fSigmaX/sum-fMeanX*fMeanX);
666     fSigmaY = TMath::Sqrt(fSigmaY/sum-fMeanY*fMeanY);   
667   }
668   else fSigmaX=0; 
669 }
670
671
672 void AliTPCPRF2D::Streamer(TBuffer &xRuub)
673 {
674    // Stream an object of class AliTPCPRF2D
675
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);
680       //read functions
681       if (strncmp(fType,"User",3)!=0){
682         delete fGRF;  
683         if (strncmp(fType,"Gauss",3)==0) 
684           fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
685         if (strncmp(fType,"Cosh",3)==0) 
686           fGRF = new TF2("FunCosh2D",FunCosh2D,-5.,5.,-5.,5.,4);
687         if (strncmp(fType,"Gati",3)==0) 
688           fGRF = new TF2("FunGati2D",FunGati2D,-5.,5.,-5.,5.,5);      
689         if (fGRF!=0) fGRF->SetParameters(funParam);
690       }
691       //calculate conversion coefitient to convert position to virtual wire
692       fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
693       fDStepM1=1/fDStep;
694    } else {
695       AliTPCPRF2D::Class()->WriteBuffer(xRuub,this);
696    }
697 }
698
699
700 TH1F *  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);
708   Float_t x=x1;
709   Float_t y1;
710
711   for (Int_t i = 0;i<kn+1;i++)
712     {
713       x+=(x2-x1)/Float_t(kn);
714       y1 = GetPRF(x,y);
715       hPRFc->Fill(x,y1);
716     };
717   hPRFc->SetXTitle("pad  (cm)");
718   return hPRFc;
719 }  
720
721 AliH2F * 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
749 AliH2F * 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
796 void 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   
805   TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
806   comment->SetTextAlign(12);
807   comment->SetFillColor(42);
808   DrawComment(comment);  
809   comment->Draw();
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  
831 }
832
833
834
835 void AliTPCPRF2D::DrawPRF(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
836
837   //
838   //
839   TCanvas  * c1 = new TCanvas("canPRF","Pad response function",700,900);
840   c1->cd();
841   TPad * pad2 = new TPad("pad2PRF","",0.05,0.22,0.95,0.95);
842   pad2->Draw(); 
843   gStyle->SetOptFit(1);
844   gStyle->SetOptStat(1); 
845   TH2F * hPRFc = GenerDrawHisto(x1, x2, y1, y2, Nx,Ny);   
846   pad2->cd();
847   hPRFc->Draw("surf");
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);
852   DrawComment(comment);  
853   comment->Draw();
854 }
855
856 void AliTPCPRF2D::DrawDist(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
857
858   //
859   //draw distortion of the COG method - for different threshold parameter
860   TCanvas  * c1 = new TCanvas("padDistortion","COG distortion",700,900);
861   c1->cd();
862   TPad * pad1 = new TPad("dist","",0.05,0.55,0.95,0.95,21);
863   pad1->Draw();
864   TPad * pad2 = new TPad("dist","",0.05,0.22,0.95,0.53,21);
865   pad2->Draw();
866   gStyle->SetOptFit(1);
867   gStyle->SetOptStat(0); 
868   
869   AliH2F * hPRFDist = GenerDrawDistHisto(x1, x2, y1, y2, Nx,Ny,thr); 
870   
871   pad1->cd();
872   hPRFDist->Draw("surf");
873   Float_t distmax =hPRFDist->GetMaximum();
874   Float_t distmin =hPRFDist->GetMinimum();
875   gStyle->SetOptStat(1); 
876   
877   TH1F * dist = hPRFDist->GetAmplitudes(distmin,distmax,distmin-1);
878   pad2->cd();
879   dist->Draw();
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);
884   DrawComment(comment);  
885   comment->Draw();
886 }
887
888 void 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:");
896   title->SetTextSize(0.03);
897   sprintf(s,"Height of pad:  %2.2f cm",fHeightFull);
898   comment->AddText(s);
899   sprintf(s,"Width pad:  %2.2f cm",fWidth);
900   comment->AddText(s);
901   sprintf(s,"Pad Angle:  %2.2f ",fPadAngle);
902   comment->AddText(s);
903   
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);    
947 }
948