Some of the coding violations corrected
[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    snprintf(fType,5,"User");
287    if (sigmaX ==0) sigmaX=(fWidth*(1+TMath::Abs(fK)))/fgkSQRT12;
288    if (sigmaY ==0) sigmaY=(fWidth*(1+TMath::Abs(fK)))/fgkSQRT12;
289    fOrigSigmaX=sigmaX; 
290    fOrigSigmaY=sigmaY; 
291    Double_t estimsigma = 
292      TMath::Sqrt(sigmaX*sigmaX+(fWidth*fWidth*(1+TMath::Abs(fK))/12)+
293                  TMath::Tan(fPadAngle*fgkDegtoRad)*TMath::Tan(fPadAngle*fgkDegtoRad)*fHeightFull*fHeightFull/12);   
294    if (estimsigma < 5*sigmaX) {
295      fDStep = estimsigma/10.;
296      fNPRF  = Int_t(estimsigma*8./fDStep); 
297    }
298    else{
299      fDStep = sigmaX; 
300      Double_t width = fWidth*(1+TMath::Abs(fK))+TMath::Abs(TMath::Tan(fPadAngle*fgkDegtoRad))*fHeightFull;
301      fNPRF = Int_t((width+8.*sigmaX)/fDStep);
302    };
303
304 }
305   
306
307 void AliTPCPRF2D::SetGauss(Float_t sigmaX, Float_t sigmaY,
308                       Float_t kNorm)
309 {
310   // 
311   // set parameters for Gauss generic charge distribution
312   //
313   fKNorm = kNorm;
314   fOrigSigmaX=sigmaX;
315   fOrigSigmaY=sigmaY;
316   //sprintf(fType,"Gauss");
317   snprintf(fType,5,"Gauss");
318   if (fGRF !=0 ) fGRF->Delete();
319   fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
320   
321   funParam[0]=sigmaX;
322   funParam[1]=sigmaY;  
323   funParam[2]=fK;
324   funParam[3]=fHeightS;    
325  
326   fGRF->SetParameters(funParam); 
327   Double_t estimsigma = 
328      TMath::Sqrt(sigmaX*sigmaX+(fWidth*fWidth*(1+TMath::Abs(fK))/12)+
329                  TMath::Tan(fPadAngle)*TMath::Tan(fPadAngle*fgkDegtoRad)*fHeightFull*fHeightFull/12);   
330    if (estimsigma < 5*sigmaX) {
331      fDStep = estimsigma/10.;
332      fNPRF  = Int_t(estimsigma*8./fDStep); 
333    }
334    else{
335      fDStep = sigmaX; 
336      Double_t width = fWidth*(1+TMath::Abs(fK))+TMath::Abs(TMath::Tan(fPadAngle*fgkDegtoRad))*fHeightFull;
337      fNPRF = Int_t((width+8.*sigmaX)/fDStep);
338    };
339  
340   
341 }
342 void AliTPCPRF2D::SetCosh(Float_t sigmaX, Float_t sigmaY,
343                      Float_t kNorm)
344
345   // set parameters for Cosh generic charge distribution
346   //
347   fKNorm = kNorm;
348   fOrigSigmaX=sigmaX;
349   fOrigSigmaY=sigmaY; 
350   //  sprintf(fType,"Cosh");
351   snprintf(fType,5,"Cosh");
352   if (fGRF !=0 ) fGRF->Delete();
353   fGRF = new TF2("FunCosh2D",   FunCosh2D,-5.,5.,-5.,5.,4);   
354   funParam[0]=sigmaX;
355   funParam[1]=sigmaY;
356   funParam[2]=fK;  
357   funParam[3]=fHeightS;
358   fGRF->SetParameters(funParam);
359
360   Double_t estimsigma = TMath::Sqrt(sigmaX*sigmaX+fWidth*fWidth*(1+TMath::Abs(fK))/12);   
361   if (estimsigma < 5*sigmaX) {
362     fDStep = estimsigma/10.;
363     fNPRF  = Int_t(estimsigma*8./fDStep); 
364   }
365   else{
366     fDStep = sigmaX; 
367     fNPRF = Int_t((1.2*fWidth*(1+TMath::Abs(fK))+8.*sigmaX)/fDStep);
368   };  
369  
370 }
371
372 void AliTPCPRF2D::SetGati(Float_t K3X, Float_t K3Y,
373                      Float_t padDistance,
374                      Float_t kNorm)
375 {
376   // set parameters for Gati generic charge distribution
377   //
378   fKNorm = kNorm;
379   fK3X=K3X;
380   fK3Y=K3Y;
381   fPadDistance=padDistance;  
382   //sprintf(fType,"Gati");
383   snprintf(fType,5,"Gati");
384   if (fGRF !=0 ) fGRF->Delete();
385   fGRF = new TF2("FunGati2D",   FunGati2D,-5.,5.,-5.,5.,5);  
386  
387   funParam[0]=padDistance;
388   funParam[1]=K3X;
389   funParam[2]=fK;  
390   funParam[3]=fHeightS;
391   funParam[4]=K3Y;
392   fGRF->SetParameters(funParam);
393   fOrigSigmaX=padDistance;
394   fOrigSigmaY=padDistance;
395   Float_t sigmaX = fOrigSigmaX;
396   Double_t estimsigma = TMath::Sqrt(sigmaX*sigmaX+fWidth*fWidth*(1+TMath::Abs(fK))/12);   
397   if (estimsigma < 5*sigmaX) {
398     fDStep = estimsigma/10.;
399     fNPRF  = Int_t(estimsigma*8./fDStep); 
400   }
401   else{
402     fDStep = sigmaX; 
403     fNPRF = Int_t((1.2*fWidth*(1+TMath::Abs(fK))+8.*sigmaX)/fDStep);
404   };
405 }
406
407
408
409 void AliTPCPRF2D::Update()
410 {
411   //
412   //update fields  with interpolated values for
413   //PRF calculation
414
415   if ( fGRF == 0 ) return;  
416   //initialize interpolated values to 0
417   Int_t i;
418   if (fChargeArray!=0) delete [] fChargeArray;
419   fChargeArray = new Float_t[fNPRF*fNYdiv];
420   fNChargeArray = fNPRF*fNYdiv;
421   for (i =0; i<fNPRF*fNYdiv;i++)  fChargeArray[i] = 0;
422   //firstly calculate total integral of charge
423
424   ////////////////////////////////////////////////////////
425   //I'm waiting for normal integral
426   //in this moment only sum
427   Float_t x2=  4*fOrigSigmaX;
428   Float_t y2=  4*fOrigSigmaY;
429   Float_t dx = fOrigSigmaX/Float_t(fNdiv*6);
430   Float_t dy = fOrigSigmaY/Float_t(fNdiv*6);  
431   Int_t nx  = Int_t(0.5+x2/dx);
432   Int_t ny  = Int_t(0.5+y2/dy);
433   Int_t ix,iy;
434   fInteg  = 0;
435   Double_t dInteg =0;
436   for (ix=-nx;ix<=nx;ix++)
437     for ( iy=-ny;iy<=ny;iy++) 
438       dInteg+=fGRF->Eval(Float_t(ix)*dx,Float_t(iy)*dy)*dx*dy;  
439   /////////////////////////////////////////////////////
440   fInteg =dInteg;
441   if ( fInteg == 0 ) fInteg = 1; 
442
443   for (i=0; i<fNYdiv; i++){
444     if (fNYdiv == 1) fCurrentY = fY1;
445     else
446       fCurrentY = fY1+Double_t(i)*(fY2-fY1)/Double_t(fNYdiv-1);
447     fcharge   = &(fChargeArray[i*fNPRF]);
448     Update1();
449   }
450   //calculate conversion coefitient to convert position to virtual wire
451   fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
452   fDStepM1=1/fDStep;
453   UpdateSigma();
454 }
455
456 void AliTPCPRF2D::Update1()
457 {
458   //
459   //update fields  with interpolated values for
460   //PRF calculation for given charge line
461   Int_t i;
462   Double_t cos = TMath::Cos(fChargeAngle);
463   Double_t sin = TMath::Sin(fChargeAngle);
464   const Double_t kprec =0.00000001;
465   //integrate charge over pad for different distance of pad
466   for (i =0; i<fNPRF;i++){      
467     //x in cm fWidth in cm
468     //calculate integral 
469     Double_t xch = fDStep * (Double_t)(i-fNPRF/2);
470     fcharge[i]=0;
471     Double_t k=1;  
472     
473     
474     for (Double_t ym=-fHeightFull/2.-fShiftY;  ym<fHeightFull/2.-kprec;ym+=fHeightS){   
475       Double_t y2chev=TMath::Min((ym+fHeightS),Double_t(fHeightFull/2.)); // end of chevron step
476       Double_t y1chev= ym;  //beginning of chevron step
477       Double_t y2 = TMath::Min(y2chev,fCurrentY+3.5*fOrigSigmaY);
478       Double_t y1 = TMath::Max((y1chev),Double_t(-fHeightFull/2.));
479       y1 = TMath::Max(y1chev,fCurrentY-3.5*fOrigSigmaY);
480
481       Double_t x0 = fWidth*(-1.-(Double_t(k)*fK))*0.5+ym*TMath::Tan(fPadAngle*fgkDegtoRad);
482       Double_t kx  = Double_t(k)*(fK*fWidth)/fHeightS;     
483       kx = TMath::Tan(TMath::ATan(kx))+TMath::Tan(fPadAngle*fgkDegtoRad);     
484
485       Int_t ny = TMath::Max(Int_t(fNdiv*TMath::Exp(-(y1-fCurrentY)*(y1-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),4);
486       Double_t dy = TMath::Min(fOrigSigmaY/Double_t(ny),y2-y1);
487       Double_t ndy = dy;
488       
489       //loop over different y strips with variable step size  dy
490       if (y2>(y1+kprec)) for (Double_t y = y1; y<y2+kprec;){      
491         //new step SIZE 
492         
493         ny = TMath::Max(Int_t(fNdiv*TMath::Exp(-(y-fCurrentY)*(y-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),5);
494         ndy = fOrigSigmaY/Double_t(ny); 
495         if (ndy>(y2-y-dy)) {
496           ndy =y2-y-dy;
497           if (ndy<kprec) ndy=2*kprec; //calculate new delta y
498         }
499         //              
500         Double_t sumch=0;
501         //calculation of x borders and initial step
502         Double_t deltay = (y-y1chev);           
503
504         Double_t xp1  = x0+deltay*kx;
505                 //x begining of pad at position y
506         Double_t xp2 =xp1+fWidth;        //x end of pad at position y
507         Double_t xp3 =xp1+kx*dy; //...at position y+dy
508         Double_t xp4 =xp2+kx*dy; //..  
509         
510         Double_t x1 = TMath::Min(xp1,xp3);
511         x1 = TMath::Max(xp1,xch-3.5*fOrigSigmaX); //beging of integration
512         Double_t x2 = TMath::Max(xp2,xp4);
513          x2 = TMath::Min(xp2+dy*kx,xch+3.5*fOrigSigmaX); //end of integration
514
515         Int_t nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x1-xch)*(x1-xch)/(2*fOrigSigmaX*fOrigSigmaX))*
516                                     TMath::Exp(-(y1-fCurrentY)*(y1-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),2);
517         Double_t dx = TMath::Min(fOrigSigmaX/Double_t(nx),x2-x1)/5.; //on the border more iteration
518         Double_t ndx=dx;
519         
520         if (x2>(x1+kprec)) {
521           for (Double_t x = x1; x<x2+kprec ;){
522           //new step SIZE         
523           nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x-xch)*(x-xch)/(2*fOrigSigmaX*fOrigSigmaX))),3);       
524           ndx = fOrigSigmaX/Double_t(nx);
525           if (ndx>(x2-x-dx)) {
526             ndx =x2-x-dx;                  
527           }
528           if ( ( (x+dx+ndx)<TMath::Max(xp3,xp1)) || ( (x+dx+ndx)>TMath::Min(xp4,xp2))) {
529             ndx/=5.;
530           }       
531           if (ndx<kprec) ndx=2*kprec;
532           //INTEGRAL APROXIMATION
533           Double_t ddx,ddy,dddx,dddy;
534           ddx = xch-(x+dx/2.);
535           ddy = fCurrentY-(y+dy/2.);
536           dddx = cos*ddx-sin*ddy;
537           dddy = sin*ddx+cos*ddy;
538           Double_t z0=fGRF->Eval(dddx,dddy);  //middle point
539           
540           ddx = xch-(x+dx/2.);
541           ddy = fCurrentY-(y);
542           dddx = cos*ddx-sin*ddy;
543           dddy = sin*ddx+cos*ddy;
544           Double_t z1=fGRF->Eval(dddx,dddy);  //point down
545           
546           ddx = xch-(x+dx/2.);
547           ddy = fCurrentY-(y+dy);
548           dddx = cos*ddx-sin*ddy;
549           dddy = sin*ddx+cos*ddy;
550           Double_t z3=fGRF->Eval(dddx,dddy);  //point up
551           
552           ddx = xch-(x);
553           ddy = fCurrentY-(y+dy/2.);
554           dddx = cos*ddx-sin*ddy;
555           dddy = sin*ddx+cos*ddy;
556           Double_t z2=fGRF->Eval(dddx,dddy);  //point left  
557           
558           ddx = xch-(x+dx);
559           ddy = fCurrentY-(y+dy/2.);
560           dddx = cos*ddx-sin*ddy;
561           dddy = sin*ddx+cos*ddy;
562           Double_t z4=fGRF->Eval(dddx,dddy);  //point right
563           
564           
565           if (z0<0) {z0=0;z1=0;z2=0;z3=0;z4=0;}
566           
567           Double_t f2x= (z3+z1-2*z0)*4.;//second derivation in y
568           Double_t f2y= (z2+z4-2*z0)*4.;//second derivation in x
569           Double_t f1y= (z3-z1);
570           Double_t z ;    
571           z = (z0+f2x/6.+f2y/6.);//second order aproxiation of integral     
572           if (kx>kprec){  //positive derivation
573             if (x<(xp1+dy*kx)){                //calculate volume at left border 
574               Double_t xx1  = x;
575               Double_t xx2  = TMath::Min(x+dx,xp1+dy*kx);
576               Double_t yy1  = y+(xx1-xp1)/kx;
577               Double_t yy2  = TMath::Min(y+(xx2-xp1)/kx,y+dy);        
578               z=z0;
579               if (yy2<y+dy) {           
580                 z-= z0*(y+dy-yy2)/dy; //constant part rectangle
581                 z-= f1y*(xx2-xx1)*(y+dy-yy2)*(y+dy-yy2)/(2.*dx*dy);
582               }
583               z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part rectangle
584               
585             }
586             if (x>xp2){          //calculate volume at right  border 
587               Double_t xx1  = x;
588               Double_t xx2  = x+dx;
589               Double_t yy1  = y+(xx1-xp2)/kx;
590               Double_t yy2  = y+(xx2-xp2)/kx;                
591               z=z0;
592               //rectangle part
593               z-=z0*(yy1-y)/dy; //constant part
594               z-=f1y*(xx2-xx1)*(yy1-y)*(yy1-y)/(2*dx*dy);
595               //triangle part         
596               z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part            
597             }
598           }       
599           if (kx<-kprec){ //negative  derivation            
600             if (x<(xp1+dy*kx)){       //calculate volume at left border          
601               Double_t xx1  = x;
602               Double_t xx2  = TMath::Min(x+dx,xp3-dy/kx);
603               Double_t yy1  = y+(xx1-xp1)/kx;
604               Double_t yy2  = TMath::Max(y,yy1+(xx2-xx1)/kx); //yy2<yy1 
605               z = z0;
606               z-= z0*(yy2-y)/dy; // constant part rectangle 
607               z-= f1y*(xx2-xx1)*(yy2-y)*(yy2-y)/(2.*dx*dy); 
608               z-=z0*(xx2-xx1)*(yy1-yy2)/(2*dx*dy); //constant part triangle
609             }
610             if (x>xp2){       //calculate volume at right  border 
611               Double_t xx1  = TMath::Max(x,xp2+dy*kx);
612               Double_t xx2  = x+dx;
613               Double_t yy1  = TMath::Min(y+dy,y-(xp2-xx1)/kx);
614               Double_t yy2  = y-(xp2-xx2)/kx;
615               z=z0;
616               z-=z0*(yy2-y)/dy;  //constant part rextangle
617               z-= f1y*(xx2-xx1)*(yy2-y)*(yy2-y)/(2.*dx*dy); 
618               z-=z0*(xx2-xx1)*(yy1-yy2)/(2*dx*dy); //constant part triangle
619             }               
620           }     
621                
622           if (z>0.)           sumch+=fKNorm*z*dx*dy/fInteg;
623           
624           x+=dx;
625           dx = ndx;
626         }; //loop over x          
627         fcharge[i]+=sumch;
628         }//if x2>x1
629         y+=dy;
630         dy =ndy;
631       }//step over different y
632       k*=-1.;
633     }//step over chevron 
634     
635    }//step over different points on line NPRF
636 }
637
638 void AliTPCPRF2D::UpdateSigma()
639 {
640   //
641   //calulate effective sigma X and sigma y of PRF
642   fMeanX = 0;
643   fMeanY = 0;
644   fSigmaX = 0;
645   fSigmaY = 0;
646  
647   Float_t sum =0;
648   Int_t i;
649   Float_t x,y;
650
651   for (i=-1; i<=fNYdiv; i++){
652     if (fNYdiv == 1) y = fY1;
653     else
654       y = fY1+Float_t(i)*(fY2-fY1)/Float_t(fNYdiv-1);
655     for (x =-fNPRF*fDStep; x<fNPRF*fDStep;x+=fDStep)
656       {      
657         //x in cm fWidth in cm
658         Float_t weight = GetPRF(x,y);
659         fSigmaX+=x*x*weight; 
660         fSigmaY+=y*y*weight;
661         fMeanX+=x*weight;
662         fMeanY+=y*weight;
663         sum+=weight;
664     };  
665   }
666   if (sum>0){
667     fMeanX/=sum;
668     fMeanY/=sum;    
669     fSigmaX = TMath::Sqrt(fSigmaX/sum-fMeanX*fMeanX);
670     fSigmaY = TMath::Sqrt(fSigmaY/sum-fMeanY*fMeanY);   
671   }
672   else fSigmaX=0; 
673 }
674
675
676 void AliTPCPRF2D::Streamer(TBuffer &xRuub)
677 {
678    // Stream an object of class AliTPCPRF2D
679
680    if (xRuub.IsReading()) {
681       UInt_t xRuus, xRuuc;
682       Version_t xRuuv = xRuub.ReadVersion(&xRuus, &xRuuc);
683       AliTPCPRF2D::Class()->ReadBuffer(xRuub, this, xRuuv, xRuus, xRuuc);
684       //read functions
685       if (strncmp(fType,"User",3)!=0){
686         delete fGRF;  
687         if (strncmp(fType,"Gauss",3)==0) 
688           fGRF = new TF2("FunGauss2D",FunGauss2D,-5.,5.,-5.,5.,4);
689         if (strncmp(fType,"Cosh",3)==0) 
690           fGRF = new TF2("FunCosh2D",FunCosh2D,-5.,5.,-5.,5.,4);
691         if (strncmp(fType,"Gati",3)==0) 
692           fGRF = new TF2("FunGati2D",FunGati2D,-5.,5.,-5.,5.,5);      
693         if (fGRF!=0) fGRF->SetParameters(funParam);
694       }
695       //calculate conversion coefitient to convert position to virtual wire
696       fDYtoWire=Float_t(fNYdiv-1)/(fY2-fY1);
697       fDStepM1=1/fDStep;
698    } else {
699       AliTPCPRF2D::Class()->WriteBuffer(xRuub,this);
700    }
701 }
702
703
704 TH1F *  AliTPCPRF2D::GenerDrawXHisto(Float_t x1, Float_t x2,Float_t y)
705 {
706   //gener one dimensional hist of pad response function
707   //  at position y 
708   char s[100]; 
709   const Int_t kn=200;
710   //sprintf(s,"Pad Response Function"); 
711   snprintf(s,100,"Pad Response Function");  
712   TH1F * hPRFc = new TH1F("hPRFc",s,kn+1,x1,x2);
713   Float_t x=x1;
714   Float_t y1;
715
716   for (Int_t i = 0;i<kn+1;i++)
717     {
718       x+=(x2-x1)/Float_t(kn);
719       y1 = GetPRF(x,y);
720       hPRFc->Fill(x,y1);
721     };
722   hPRFc->SetXTitle("pad  (cm)");
723   return hPRFc;
724 }  
725
726 AliH2F * AliTPCPRF2D::GenerDrawHisto(Float_t x1, Float_t x2, Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
727 {
728   //
729   //gener two dimensional histogram with PRF
730   //
731   char s[100];
732   //sprintf(s,"Pad Response Function"); 
733   snprintf(s,100,"Pad Response Function");  
734   AliH2F * hPRFc = new AliH2F("hPRFc",s,Nx,x1,x2,Ny,y1,y2);
735   Float_t dx=(x2-x1)/Float_t(Nx);
736   Float_t dy=(y2-y1)/Float_t(Ny) ;
737   Float_t x,y,z; 
738   x = x1;
739   y = y1;
740   for ( Int_t i  = 0;i<=Nx;i++,x+=dx){
741     y=y1;
742     for (Int_t j  = 0;j<=Ny;j++,y+=dy){
743       z = GetPRF(x,y);
744       hPRFc->SetCellContent(i,j,z);
745     };
746   }; 
747   hPRFc->SetXTitle("pad direction (cm)");
748   hPRFc->SetYTitle("pad row  direction (cm)");
749   hPRFc->SetTitleOffset(1.5,"X");
750   hPRFc->SetTitleOffset(1.5,"Y");
751   return hPRFc;
752 }
753
754
755 AliH2F * AliTPCPRF2D::GenerDrawDistHisto(Float_t x1, Float_t x2, Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t  thr)
756 {
757   //return histogram with distortion
758   const Float_t kminth=0.00001;
759   if (thr<kminth) thr=kminth;
760   char s[100]; 
761   //sprintf(s,"COG distortion of PRF (threshold=%2.2f)",thr); 
762   snprintf(s,100,"COG distortion of PRF (threshold=%2.2f)",thr); 
763   AliH2F * hPRFDist = new AliH2F("hDistortion",s,Nx,x1,x2,Ny,y1,y2);
764   Float_t dx=(x2-x1)/Float_t(Nx);
765   Float_t dy=(y2-y1)/Float_t(Ny) ;
766   Float_t x,y,z,ddx;
767   x=x1;
768   for ( Int_t i  = 0;i<=Nx;i++,x+=dx){
769     y=y1;
770     for(Int_t j  = 0;j<=Ny;j++,y+=dy)      
771       {
772         Float_t sumx=0;
773         Float_t sum=0;
774         for (Int_t k=-3;k<=3;k++)
775           {         
776             Float_t padx=Float_t(k)*fWidth;
777             z = GetPRF(x-padx,y); 
778             if (z>thr){
779               sum+=z;
780               sumx+=z*padx;
781             }   
782           };    
783         if (sum>kminth)  
784           {
785             ddx = (x-(sumx/sum));
786           }
787         else ddx=-1;
788         if (TMath::Abs(ddx)<10)         hPRFDist->SetCellContent(i,j,ddx);
789       }
790   }
791
792   hPRFDist->SetXTitle("pad direction (cm)");
793   hPRFDist->SetYTitle("pad row  direction (cm)");
794   hPRFDist->SetTitleOffset(1.5,"X");
795   hPRFDist->SetTitleOffset(1.5,"Y");
796   return hPRFDist;
797 }  
798   
799
800
801
802
803 void AliTPCPRF2D::DrawX(Float_t x1 ,Float_t x2,Float_t y1,Float_t y2, Int_t N)
804
805   //
806   //draw pad response function at interval <x1,x2> at  given y position
807   //
808   if (N<0) return;
809   TCanvas  * c1 = new TCanvas("PRFX","Pad response function",700,900);
810   c1->cd();  
811   
812   TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
813   comment->SetTextAlign(12);
814   comment->SetFillColor(42);
815   DrawComment(comment);  
816   comment->Draw();
817   c1->cd(); 
818
819   TPad * pad2 = new TPad("pPRF","",0.05,0.22,0.95,0.95);
820   pad2->Divide(2,(N+1)/2);
821   pad2->Draw();
822   gStyle->SetOptFit(1);
823   gStyle->SetOptStat(1); 
824   for (Int_t i=0;i<N;i++){
825     char ch[200];
826     Float_t y;
827     if (N==1) y=y1;
828     else y = y1+i*(y2-y1)/Float_t(N-1);
829     pad2->cd(i+1);
830     TH1F * hPRFc =GenerDrawXHisto(x1, x2,y);
831     //sprintf(ch,"PRF at wire position: %2.3f",y);
832     snprintf(ch,40,"PRF at wire position: %2.3f",y);
833     hPRFc->SetTitle(ch);  
834     //sprintf(ch,"PRF %d",i);
835     snprintf(ch,15,"PRF %d",i);
836     hPRFc->SetName(ch);  
837      hPRFc->Fit("gaus");
838   }
839  
840 }
841
842
843
844 void AliTPCPRF2D::DrawPRF(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny)
845
846   //
847   //
848   TCanvas  * c1 = new TCanvas("canPRF","Pad response function",700,900);
849   c1->cd();
850   TPad * pad2 = new TPad("pad2PRF","",0.05,0.22,0.95,0.95);
851   pad2->Draw(); 
852   gStyle->SetOptFit(1);
853   gStyle->SetOptStat(1); 
854   TH2F * hPRFc = GenerDrawHisto(x1, x2, y1, y2, Nx,Ny);   
855   pad2->cd();
856   hPRFc->Draw("surf");
857   c1->cd(); 
858   TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
859   comment->SetTextAlign(12);
860   comment->SetFillColor(42);
861   DrawComment(comment);  
862   comment->Draw();
863 }
864
865 void AliTPCPRF2D::DrawDist(Float_t x1 ,Float_t x2,Float_t y1, Float_t y2, Int_t Nx, Int_t Ny, Float_t thr)
866
867   //
868   //draw distortion of the COG method - for different threshold parameter
869   TCanvas  * c1 = new TCanvas("padDistortion","COG distortion",700,900);
870   c1->cd();
871   TPad * pad1 = new TPad("dist","",0.05,0.55,0.95,0.95,21);
872   pad1->Draw();
873   TPad * pad2 = new TPad("dist","",0.05,0.22,0.95,0.53,21);
874   pad2->Draw();
875   gStyle->SetOptFit(1);
876   gStyle->SetOptStat(0); 
877   
878   AliH2F * hPRFDist = GenerDrawDistHisto(x1, x2, y1, y2, Nx,Ny,thr); 
879   
880   pad1->cd();
881   hPRFDist->Draw("surf");
882   Float_t distmax =hPRFDist->GetMaximum();
883   Float_t distmin =hPRFDist->GetMinimum();
884   gStyle->SetOptStat(1); 
885   
886   TH1F * dist = hPRFDist->GetAmplitudes(distmin,distmax,distmin-1);
887   pad2->cd();
888   dist->Draw();
889   c1->cd(); 
890   TPaveText * comment = new TPaveText(0.05,0.02,0.95,0.20,"NDC");
891   comment->SetTextAlign(12);
892   comment->SetFillColor(42);
893   DrawComment(comment);  
894   comment->Draw();
895 }
896
897 void AliTPCPRF2D::DrawComment(TPaveText *comment)
898 {
899   //
900   //function to write comment to picture 
901   
902   char s[100];
903   //draw comments to picture
904   TText * title = comment->AddText("Pad Response Function  parameters:");
905   title->SetTextSize(0.03);
906   //sprintf(s,"Height of pad:  %2.2f cm",fHeightFull);
907   snprintf(s,100,"Height of pad:  %2.2f cm",fHeightFull);
908   comment->AddText(s);
909   //sprintf(s,"Width pad:  %2.2f cm",fWidth);
910   snprintf(s,100,"Width pad:  %2.2f cm",fWidth);
911   comment->AddText(s);
912   //sprintf(s,"Pad Angle:  %2.2f ",fPadAngle);
913   snprintf(s,100,"Pad Angle:  %2.2f ",fPadAngle);
914   comment->AddText(s);
915   
916   if (TMath::Abs(fK)>0.0001){
917     //sprintf(s,"Height of one chevron unit h:  %2.2f cm",2*fHeightS);
918     snprintf(s,100,"Height of one chevron unit h:  %2.2f cm",2*fHeightS);
919     comment->AddText(s);
920     //sprintf(s,"Overlap factor:  %2.2f",fK);
921     snprintf(s,100,"Overlap factor:  %2.2f",fK);
922     comment->AddText(s); 
923   }
924
925   if (strncmp(fType,"User",3)==0){
926     //sprintf(s,"Charge distribution - user defined function  %s ",fGRF->GetTitle());
927     snprintf(s,100,"Charge distribution - user defined function  %s ",fGRF->GetTitle());
928     comment->AddText(s);  
929     //sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
930     snprintf(s,100,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);  
931     comment->AddText(s);  
932     //sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
933     snprintf(s,100,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
934     comment->AddText(s); 
935   }
936   if (strncmp(fType,"Gauss",3)==0){
937     //sprintf(s,"Gauss charge distribution");
938     snprintf(s,100,"Gauss charge distribution");
939     comment->AddText(s);  
940     //sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
941     snprintf(s,100,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
942     comment->AddText(s);  
943     //sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
944     snprintf(s,100,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
945     comment->AddText(s); 
946   }
947   if (strncmp(fType,"Gati",3)==0){
948     //sprintf(s,"Gati charge distribution");
949     snprintf(s,100,"Gati charge distribution");
950     comment->AddText(s);  
951     //sprintf(s,"K3X of Gati : %2.2f ",fK3X);
952     snprintf(s,100,"K3X of Gati : %2.2f ",fK3X);
953     comment->AddText(s);  
954     //sprintf(s,"K3Y of Gati: %2.2f ",fK3Y);
955     snprintf(s,100,"K3Y of Gati: %2.2f ",fK3Y);
956     comment->AddText(s); 
957     //sprintf(s,"Wire to Pad Distance: %2.2f ",fPadDistance);
958     snprintf(s,100,"Wire to Pad Distance: %2.2f ",fPadDistance);
959     comment->AddText(s); 
960   }
961   if (strncmp(fType,"Cosh",3)==0){
962     //sprintf(s,"Cosh charge distribution");
963     snprintf(s,100,"Cosh charge distribution");
964     comment->AddText(s);  
965     //sprintf(s,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
966     snprintf(s,100,"Sigma x of charge distribution: %2.2f ",fOrigSigmaX);
967     comment->AddText(s);  
968     //sprintf(s,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
969     snprintf(s,100,"Sigma y of charge distribution: %2.2f ",fOrigSigmaY);
970     comment->AddText(s); 
971   }
972   //sprintf(s,"Normalisation: %2.2f ",fKNorm);
973   snprintf(s,100,"Normalisation: %2.2f ",fKNorm);
974   comment->AddText(s);    
975 }
976