TPC/AliTPCPRF2D.cxx

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