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