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