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