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