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