TPC/AliTPCClusterParam.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  TPC cluster error, shape and charge parameterization as function
  20 //  of drift length, and inclination angle                                   //
  21 //
  22 //  Following notation is used in following
  23 //  Int_t dim 0 - y direction
  24 //            1 - z direction
  25 //
  26 //  Int_t type 0 - short pads
  27 //             1 - medium pads
  28 //             2 - long pads
  29 //  Float_t z    - drift length
  30 //
  31 //  Float_t angle - tangent of inclination angle at given dimension
  32 //
  33 //  Implemented parameterization
  34 //
  35 //
  36 //  1. Resolution as function of drift length and inclination angle
  37 //     1.a) GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle)
  38 //          Simple error parameterization as derived from analytical formula
  39 //          only linear term in drift length and angle^2
  40 //          The formula is valid only with precission +-5%
  41 //          Separate parameterization for differnt pad geometry
  42 //     1.b) GetError0Par
  43 //          Parabolic term correction - better precision
  44 //
  45 //     1.c) GetError1 - JUST FOR Study
  46 //          Similar to GetError1
  47 //          The angular and diffusion effect is scaling with pad length
  48 //          common parameterization for different pad length
  49 //
  50 //  2. Error parameterization using charge
  51 //     2.a) GetErrorQ
  52 //          GetError0+
  53 //          adding 1/Q component to diffusion and angluar part
  54 //     2.b) GetErrorQPar
  55 //          GetError0Par+
  56 //          adding 1/Q component to diffusion and angluar part
  57 //     2.c) GetErrorQParScaled - Just for study
  58 //          One parameterization for all pad shapes
  59 //          Smaller precission as previous one
  60 //
  61 //
  62 //  Example how to retrieve the paramterization:
  63 /*
  64       AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  65       AliCDBManager::Instance()->SetRun(0)
  66       AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
  67
  68       //
  69       //
  70       AliTPCClusterParam::SetInstance(param);
  71       TF1 f1("f1","AliTPCClusterParam::SGetError0Par(1,0,x,0)",0,250);
  72 */
  73
  74 // EXAMPLE hot to create parameterization
  75 /*
  76 // Note resol is the resolution tree created by AliTPCcalibTracks
  77 //
  78 AliTPCClusterParam  *param = new AliTPCClusterParam;
  79 param->FitData(Resol);
  80 AliTPCClusterParam::SetInstance(param);
  81
  82 */
  83
  84 //
  85 //                                                                     //
  86 ///////////////////////////////////////////////////////////////////////////////
  87 #include "AliTPCClusterParam.h"
  88 #include "TMath.h"
  89 #include "TFile.h"
  90 #include "TTree.h"
  91 #include <TVectorF.h>
  92 #include <TLinearFitter.h>
  93 #include <TH1F.h>
  94 #include <TProfile2D.h>
  95 #include <TVectorD.h>
  96 #include <TObjArray.h>
  97 #include "AliTPCcalibDB.h"
  98 #include "AliTPCParam.h"
  99
 100 ClassImp(AliTPCClusterParam)
 101
 102
 103 AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
 104
 105
 106 /*
 107   Example usage fitting parameterization:
 108   TFile fres("resol.root");    //tree with resolution and shape
 109   TTree * treeRes =(TTree*)fres.Get("Resol");
 110
 111   AliTPCClusterParam param;
 112   param.SetInstance(&param);
 113   param.FitResol(treeRes);
 114   param.FitRMS(treeRes);
 115   TFile fparam("TPCClusterParam.root","recreate");
 116   param.Write("Param");
 117   //
 118   //
 119   TFile fparam("TPCClusterParam.root");
 120   AliTPCClusterParam *param2  =  (AliTPCClusterParam *) fparam.Get("Param");
 121   param2->SetInstance(param2);
 122   param2->Test(treeRes);
 123
 124
 125   treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
 126
 127 */
 128
 129
 130
 131
 132 //_ singleton implementation __________________________________________________
 133 AliTPCClusterParam* AliTPCClusterParam::Instance()
 134 {
 135   //
 136   // Singleton implementation
 137   // Returns an instance of this class, it is created if neccessary
 138   //
 139   if (fgInstance == 0){
 140     fgInstance = new AliTPCClusterParam();
 141   }
 142   return fgInstance;
 143 }
 144
 145
 146 AliTPCClusterParam::AliTPCClusterParam():
 147   TObject(),
 148   fRatio(0),
 149   fQNorm(0),
 150   fQpadTnorm(0),           // q pad normalization - Total charge
 151   fQpadMnorm(0)            // q pad normalization - Max charge
 152   //
 153 {
 154   //
 155   // Default constructor
 156   //
 157   fPosQTnorm[0] = 0;   fPosQTnorm[1] = 0;   fPosQTnorm[2] = 0;
 158   fPosQMnorm[0] = 0;   fPosQMnorm[1] = 0;   fPosQMnorm[2] = 0;
 159   //
 160   fPosYcor[0]   = 0;   fPosYcor[1]   = 0;   fPosYcor[2]   = 0;
 161   fPosZcor[0]   = 0;   fPosZcor[1]   = 0;   fPosZcor[2]   = 0;
 162 }
 163
 164 AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
 165   TObject(param),
 166   fRatio(0),
 167   fQNorm(0),
 168   fQpadTnorm(new TVectorD(*(param.fQpadTnorm))),           // q pad normalization - Total charge
 169   fQpadMnorm(new TVectorD(*(param.fQpadMnorm)))            // q pad normalization - Max charge
 170
 171 {
 172   //
 173   // copy constructor
 174   //
 175   memcpy(this, &param,sizeof(AliTPCClusterParam));
 176   if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
 177   //
 178   if (param.fPosQTnorm[0]){
 179     fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
 180     fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
 181     fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
 182     //
 183     fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
 184     fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
 185     fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
 186   }
 187   if (param.fPosYcor[0]){
 188     fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
 189     fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
 190     fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
 191     //
 192     fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
 193     fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
 194     fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
 195   }
 196
 197 }
 198
 199
 200 AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
 201   //
 202   // Assignment operator
 203   //
 204   if (this != &param) {
 205     memcpy(this, &param,sizeof(AliTPCClusterParam));
 206     if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
 207     if (param.fPosQTnorm[0]){
 208       fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
 209       fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
 210       fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
 211       //
 212       fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
 213       fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
 214       fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
 215     }
 216     if (param.fPosYcor[0]){
 217       fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
 218       fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
 219       fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
 220       //
 221       fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
 222       fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
 223       fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
 224     }
 225   }
 226   return *this;
 227 }
 228
 229
 230 AliTPCClusterParam::~AliTPCClusterParam(){
 231   //
 232   // destructor
 233   //
 234   if (fQNorm) fQNorm->Delete();
 235   delete fQNorm;
 236   if (fPosQTnorm[0]){
 237     delete fPosQTnorm[0];
 238     delete fPosQTnorm[1];
 239     delete fPosQTnorm[2];
 240     //
 241     delete fPosQMnorm[0];
 242     delete fPosQMnorm[1];
 243     delete fPosQMnorm[2];
 244   }
 245   if (fPosYcor[0]){
 246     delete fPosYcor[0];
 247     delete fPosYcor[1];
 248     delete fPosYcor[2];
 249     //
 250     delete fPosZcor[0];
 251     delete fPosZcor[1];
 252     delete fPosZcor[2];
 253   }
 254 }
 255
 256
 257 void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 258   //
 259   // Fit z - angular dependence of resolution
 260   //
 261   // Int_t dim=0, type=0;
 262   char varVal[100];
 263   sprintf(varVal,"Resol:AngleM:Zm");
 264   char varErr[100];
 265   sprintf(varErr,"Sigma:AngleS:Zs");
 266   char varCut[100];
 267   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 268   //
 269   Int_t entries = tree->Draw(varVal,varCut);
 270   Float_t px[10000], py[10000], pz[10000];
 271   Float_t ex[10000], ey[10000], ez[10000];
 272   //
 273   tree->Draw(varErr,varCut);
 274   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 275     ex[ipoint]= tree->GetV3()[ipoint];
 276     ey[ipoint]= tree->GetV2()[ipoint];
 277     ez[ipoint]= tree->GetV1()[ipoint];
 278   }
 279   tree->Draw(varVal,varCut);
 280   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 281     px[ipoint]= tree->GetV3()[ipoint];
 282     py[ipoint]= tree->GetV2()[ipoint];
 283     pz[ipoint]= tree->GetV1()[ipoint];
 284   }
 285
 286   //
 287   TLinearFitter fitter(3,"hyp2");
 288   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 289     Float_t val = pz[ipoint]*pz[ipoint];
 290     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 291     Double_t x[2];
 292     x[0] = px[ipoint];
 293     x[1] = py[ipoint]*py[ipoint];
 294     fitter.AddPoint(x,val,err);
 295   }
 296   fitter.Eval();
 297   TVectorD param(3);
 298   fitter.GetParameters(param);
 299   param0[0] = param[0];
 300   param0[1] = param[1];
 301   param0[2] = param[2];
 302   Float_t chi2 =  fitter.GetChisquare()/entries;
 303   param0[3] = chi2;
 304   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 305   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 306   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 307 }
 308
 309
 310 void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 311   //
 312   // Fit z - angular dependence of resolution
 313   //
 314   // Int_t dim=0, type=0;
 315   char varVal[100];
 316   sprintf(varVal,"Resol:AngleM:Zm");
 317   char varErr[100];
 318   sprintf(varErr,"Sigma:AngleS:Zs");
 319   char varCut[100];
 320   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 321   //
 322   Int_t entries = tree->Draw(varVal,varCut);
 323   Float_t px[10000], py[10000], pz[10000];
 324   Float_t ex[10000], ey[10000], ez[10000];
 325   //
 326   tree->Draw(varErr,varCut);
 327   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 328     ex[ipoint]= tree->GetV3()[ipoint];
 329     ey[ipoint]= tree->GetV2()[ipoint];
 330     ez[ipoint]= tree->GetV1()[ipoint];
 331   }
 332   tree->Draw(varVal,varCut);
 333   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 334     px[ipoint]= tree->GetV3()[ipoint];
 335     py[ipoint]= tree->GetV2()[ipoint];
 336     pz[ipoint]= tree->GetV1()[ipoint];
 337   }
 338
 339   //
 340   TLinearFitter fitter(6,"hyp5");
 341   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 342     Float_t val = pz[ipoint]*pz[ipoint];
 343     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 344     Double_t x[6];
 345     x[0] = px[ipoint];
 346     x[1] = py[ipoint]*py[ipoint];
 347     x[2] = x[0]*x[0];
 348     x[3] = x[1]*x[1];
 349     x[4] = x[0]*x[1];
 350     fitter.AddPoint(x,val,err);
 351   }
 352   fitter.Eval();
 353   TVectorD param(6);
 354   fitter.GetParameters(param);
 355   param0[0] = param[0];
 356   param0[1] = param[1];
 357   param0[2] = param[2];
 358   param0[3] = param[3];
 359   param0[4] = param[4];
 360   param0[5] = param[5];
 361   Float_t chi2 =  fitter.GetChisquare()/entries;
 362   param0[6] = chi2;
 363   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 364   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 365   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 366   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 367   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 368   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
 369 }
 370
 371
 372
 373
 374
 375 void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
 376   //
 377   // Fit z - angular dependence of resolution - pad length scaling
 378   //
 379   // Int_t dim=0, type=0;
 380   char varVal[100];
 381   sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
 382   char varErr[100];
 383   sprintf(varErr,"Sigma:AngleS:Zs");
 384   char varCut[100];
 385   sprintf(varCut,"Dim==%d&&QMean<0",dim);
 386   //
 387   Int_t entries = tree->Draw(varVal,varCut);
 388   Float_t px[10000], py[10000], pz[10000];
 389   Float_t ex[10000], ey[10000], ez[10000];
 390   //
 391   tree->Draw(varErr,varCut);
 392   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 393     ex[ipoint]= tree->GetV3()[ipoint];
 394     ey[ipoint]= tree->GetV2()[ipoint];
 395     ez[ipoint]= tree->GetV1()[ipoint];
 396   }
 397   tree->Draw(varVal,varCut);
 398   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 399     px[ipoint]= tree->GetV3()[ipoint];
 400     py[ipoint]= tree->GetV2()[ipoint];
 401     pz[ipoint]= tree->GetV1()[ipoint];
 402   }
 403
 404   //
 405   TLinearFitter fitter(3,"hyp2");
 406   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 407     Float_t val = pz[ipoint]*pz[ipoint];
 408     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 409     Double_t x[2];
 410     x[0] = px[ipoint];
 411     x[1] = py[ipoint]*py[ipoint];
 412     fitter.AddPoint(x,val,err);
 413   }
 414   fitter.Eval();
 415   TVectorD param(3);
 416   fitter.GetParameters(param);
 417   param0[0] = param[0];
 418   param0[1] = param[1];
 419   param0[2] = param[2];
 420   Float_t chi2 =  fitter.GetChisquare()/entries;
 421   param0[3] = chi2;
 422   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 423   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 424   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 425 }
 426
 427 void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 428   //
 429   // Fit z - angular dependence of resolution - Q scaling
 430   //
 431   // Int_t dim=0, type=0;
 432   char varVal[100];
 433   sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
 434   char varVal0[100];
 435   sprintf(varVal0,"Resol:AngleM:Zm");
 436   //
 437   char varErr[100];
 438   sprintf(varErr,"Sigma:AngleS:Zs");
 439   char varCut[100];
 440   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 441   //
 442   Int_t entries = tree->Draw(varVal,varCut);
 443   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 444   Float_t ex[20000], ey[20000], ez[20000];
 445   //
 446   tree->Draw(varErr,varCut);
 447   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 448     ex[ipoint]= tree->GetV3()[ipoint];
 449     ey[ipoint]= tree->GetV2()[ipoint];
 450     ez[ipoint]= tree->GetV1()[ipoint];
 451   }
 452   tree->Draw(varVal,varCut);
 453   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 454     px[ipoint]= tree->GetV3()[ipoint];
 455     py[ipoint]= tree->GetV2()[ipoint];
 456     pz[ipoint]= tree->GetV1()[ipoint];
 457   }
 458   tree->Draw(varVal0,varCut);
 459   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 460     pu[ipoint]= tree->GetV3()[ipoint];
 461     pt[ipoint]= tree->GetV2()[ipoint];
 462   }
 463
 464   //
 465   TLinearFitter fitter(5,"hyp4");
 466   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 467     Float_t val = pz[ipoint]*pz[ipoint];
 468     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 469     Double_t x[4];
 470     x[0] = pu[ipoint];
 471     x[1] = pt[ipoint]*pt[ipoint];
 472     x[2] = px[ipoint];
 473     x[3] = py[ipoint]*py[ipoint];
 474     fitter.AddPoint(x,val,err);
 475   }
 476
 477   fitter.Eval();
 478   TVectorD param(5);
 479   fitter.GetParameters(param);
 480   param0[0] = param[0];
 481   param0[1] = param[1];
 482   param0[2] = param[2];
 483   param0[3] = param[3];
 484   param0[4] = param[4];
 485   Float_t chi2 =  fitter.GetChisquare()/entries;
 486   param0[5] = chi2;
 487   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 488   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 489   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 490   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 491   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 492 }
 493
 494 void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 495   //
 496   // Fit z - angular dependence of resolution - Q scaling  - parabolic correction
 497   //
 498   // Int_t dim=0, type=0;
 499   char varVal[100];
 500   sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
 501   char varVal0[100];
 502   sprintf(varVal0,"Resol:AngleM:Zm");
 503   //
 504   char varErr[100];
 505   sprintf(varErr,"Sigma:AngleS:Zs");
 506   char varCut[100];
 507   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 508   //
 509   Int_t entries = tree->Draw(varVal,varCut);
 510   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 511   Float_t ex[20000], ey[20000], ez[20000];
 512   //
 513   tree->Draw(varErr,varCut);
 514   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 515     ex[ipoint]= tree->GetV3()[ipoint];
 516     ey[ipoint]= tree->GetV2()[ipoint];
 517     ez[ipoint]= tree->GetV1()[ipoint];
 518   }
 519   tree->Draw(varVal,varCut);
 520   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 521     px[ipoint]= tree->GetV3()[ipoint];
 522     py[ipoint]= tree->GetV2()[ipoint];
 523     pz[ipoint]= tree->GetV1()[ipoint];
 524   }
 525   tree->Draw(varVal0,varCut);
 526   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 527     pu[ipoint]= tree->GetV3()[ipoint];
 528     pt[ipoint]= tree->GetV2()[ipoint];
 529   }
 530
 531   //
 532   TLinearFitter fitter(8,"hyp7");
 533   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 534     Float_t val = pz[ipoint]*pz[ipoint];
 535     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 536     Double_t x[7];
 537     x[0] = pu[ipoint];
 538     x[1] = pt[ipoint]*pt[ipoint];
 539     x[2] = x[0]*x[0];
 540     x[3] = x[1]*x[1];
 541     x[4] = x[0]*x[1];
 542     x[5] = px[ipoint];
 543     x[6] = py[ipoint]*py[ipoint];
 544     //
 545     fitter.AddPoint(x,val,err);
 546   }
 547
 548   fitter.Eval();
 549   TVectorD param(8);
 550   fitter.GetParameters(param);
 551   param0[0] = param[0];
 552   param0[1] = param[1];
 553   param0[2] = param[2];
 554   param0[3] = param[3];
 555   param0[4] = param[4];
 556   param0[5] = param[5];
 557   param0[6] = param[6];
 558   param0[7] = param[7];
 559
 560   Float_t chi2 =  fitter.GetChisquare()/entries;
 561   param0[8] = chi2;
 562   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 563   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 564   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 565   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 566   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 567   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
 568   error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
 569   error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
 570 }
 571
 572
 573
 574 void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 575   //
 576   // Fit z - angular dependence of resolution
 577   //
 578   // Int_t dim=0, type=0;
 579   char varVal[100];
 580   sprintf(varVal,"RMSm:AngleM:Zm");
 581   char varErr[100];
 582   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
 583   char varCut[100];
 584   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 585   //
 586   Int_t entries = tree->Draw(varVal,varCut);
 587   Float_t px[10000], py[10000], pz[10000];
 588   Float_t ex[10000], ey[10000], ez[10000];
 589   //
 590   tree->Draw(varErr,varCut);
 591   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 592     ex[ipoint]= tree->GetV3()[ipoint];
 593     ey[ipoint]= tree->GetV2()[ipoint];
 594     ez[ipoint]= tree->GetV1()[ipoint];
 595   }
 596   tree->Draw(varVal,varCut);
 597   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 598     px[ipoint]= tree->GetV3()[ipoint];
 599     py[ipoint]= tree->GetV2()[ipoint];
 600     pz[ipoint]= tree->GetV1()[ipoint];
 601   }
 602
 603   //
 604   TLinearFitter fitter(3,"hyp2");
 605   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 606     Float_t val = pz[ipoint]*pz[ipoint];
 607     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 608     Double_t x[2];
 609     x[0] = px[ipoint];
 610     x[1] = py[ipoint]*py[ipoint];
 611     fitter.AddPoint(x,val,err);
 612   }
 613   fitter.Eval();
 614   TVectorD param(3);
 615   fitter.GetParameters(param);
 616   param0[0] = param[0];
 617   param0[1] = param[1];
 618   param0[2] = param[2];
 619   Float_t chi2 =  fitter.GetChisquare()/entries;
 620   param0[3] = chi2;
 621   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 622   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 623   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 624 }
 625
 626 void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
 627   //
 628   // Fit z - angular dependence of resolution - pad length scaling
 629   //
 630   // Int_t dim=0, type=0;
 631   char varVal[100];
 632   sprintf(varVal,"RMSm:AngleM*Length:Zm");
 633   char varErr[100];
 634   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
 635   char varCut[100];
 636   sprintf(varCut,"Dim==%d&&QMean<0",dim);
 637   //
 638   Int_t entries = tree->Draw(varVal,varCut);
 639   Float_t px[10000], py[10000], pz[10000];
 640   Float_t type[10000], ey[10000], ez[10000];
 641   //
 642   tree->Draw(varErr,varCut);
 643   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 644     type[ipoint] = tree->GetV3()[ipoint];
 645     ey[ipoint]   = tree->GetV2()[ipoint];
 646     ez[ipoint]   = tree->GetV1()[ipoint];
 647   }
 648   tree->Draw(varVal,varCut);
 649   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 650     px[ipoint]= tree->GetV3()[ipoint];
 651     py[ipoint]= tree->GetV2()[ipoint];
 652     pz[ipoint]= tree->GetV1()[ipoint];
 653   }
 654
 655   //
 656   TLinearFitter fitter(4,"hyp3");
 657   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 658     Float_t val = pz[ipoint]*pz[ipoint];
 659     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 660     Double_t x[3];
 661     x[0] = (type[ipoint]<0.5)? 0.:1.;
 662     x[1] = px[ipoint];
 663     x[2] = py[ipoint]*py[ipoint];
 664     fitter.AddPoint(x,val,err);
 665   }
 666   fitter.Eval();
 667   TVectorD param(4);
 668   fitter.GetParameters(param);
 669   param0[0] = param[0];
 670   param0[1] = param[0]+param[1];
 671   param0[2] = param[2];
 672   param0[3] = param[3];
 673   Float_t chi2 =  fitter.GetChisquare()/entries;
 674   param0[4] = chi2;
 675   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 676   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 677   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 678   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 679 }
 680
 681 void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 682   //
 683   // Fit z - angular dependence of resolution - Q scaling
 684   //
 685   // Int_t dim=0, type=0;
 686   char varVal[100];
 687   sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
 688   char varVal0[100];
 689   sprintf(varVal0,"RMSm:AngleM:Zm");
 690   //
 691   char varErr[100];
 692   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
 693   char varCut[100];
 694   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 695   //
 696   Int_t entries = tree->Draw(varVal,varCut);
 697   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 698   Float_t ex[20000], ey[20000], ez[20000];
 699   //
 700   tree->Draw(varErr,varCut);
 701   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 702     ex[ipoint]= tree->GetV3()[ipoint];
 703     ey[ipoint]= tree->GetV2()[ipoint];
 704     ez[ipoint]= tree->GetV1()[ipoint];
 705   }
 706   tree->Draw(varVal,varCut);
 707   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 708     px[ipoint]= tree->GetV3()[ipoint];
 709     py[ipoint]= tree->GetV2()[ipoint];
 710     pz[ipoint]= tree->GetV1()[ipoint];
 711   }
 712   tree->Draw(varVal0,varCut);
 713   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 714     pu[ipoint]= tree->GetV3()[ipoint];
 715     pt[ipoint]= tree->GetV2()[ipoint];
 716   }
 717
 718   //
 719   TLinearFitter fitter(5,"hyp4");
 720   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 721     Float_t val = pz[ipoint]*pz[ipoint];
 722     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 723     Double_t x[4];
 724     x[0] = pu[ipoint];
 725     x[1] = pt[ipoint]*pt[ipoint];
 726     x[2] = px[ipoint];
 727     x[3] = py[ipoint]*py[ipoint];
 728     fitter.AddPoint(x,val,err);
 729   }
 730
 731   fitter.Eval();
 732   TVectorD param(5);
 733   fitter.GetParameters(param);
 734   param0[0] = param[0];
 735   param0[1] = param[1];
 736   param0[2] = param[2];
 737   param0[3] = param[3];
 738   param0[4] = param[4];
 739   Float_t chi2 =  fitter.GetChisquare()/entries;
 740   param0[5] = chi2;
 741   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 742   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 743   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 744   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 745   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 746 }
 747
 748
 749 void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
 750   //
 751   // Fit z - angular dependence of resolution - Q scaling
 752   //
 753   // Int_t dim=0, type=0;
 754   char varVal[100];
 755   sprintf(varVal,"RMSs:RMSm");
 756   //
 757   char varCut[100];
 758   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 759   //
 760   Int_t entries = tree->Draw(varVal,varCut);
 761   Float_t px[20000], py[20000];
 762   //
 763   tree->Draw(varVal,varCut);
 764   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 765     px[ipoint]= tree->GetV2()[ipoint];
 766     py[ipoint]= tree->GetV1()[ipoint];
 767   }
 768   TLinearFitter fitter(2,"pol1");
 769   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 770     Float_t val = py[ipoint];
 771     Float_t err = fRatio*px[ipoint];
 772     Double_t x[4];
 773     x[0] = px[ipoint];
 774     if (err>0) fitter.AddPoint(x,val,err);
 775   }
 776   fitter.Eval();
 777   param0[0]= fitter.GetParameter(0);
 778   param0[1]= fitter.GetParameter(1);
 779 }
 780
 781
 782
 783 Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
 784   //
 785   //
 786   //
 787   Float_t value=0;
 788   value += fParamS0[dim][type][0];
 789   value += fParamS0[dim][type][1]*z;
 790   value += fParamS0[dim][type][2]*angle*angle;
 791   value  = TMath::Sqrt(TMath::Abs(value));
 792   return value;
 793 }
 794
 795
 796 Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
 797   //
 798   //
 799   //
 800   Float_t value=0;
 801   value += fParamS0Par[dim][type][0];
 802   value += fParamS0Par[dim][type][1]*z;
 803   value += fParamS0Par[dim][type][2]*angle*angle;
 804   value += fParamS0Par[dim][type][3]*z*z;
 805   value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
 806   value += fParamS0Par[dim][type][5]*z*angle*angle;
 807   value  = TMath::Sqrt(TMath::Abs(value));
 808   return value;
 809 }
 810
 811
 812
 813 Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
 814   //
 815   //
 816   //
 817   Float_t value=0;
 818   Float_t length=0.75;
 819   if (type==1) length=1;
 820   if (type==2) length=1.5;
 821   value += fParamS1[dim][0];
 822   value += fParamS1[dim][1]*z/length;
 823   value += fParamS1[dim][2]*angle*angle*length;
 824   value  = TMath::Sqrt(TMath::Abs(value));
 825   return value;
 826 }
 827
 828 Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 829   //
 830   //
 831   //
 832   Float_t value=0;
 833   value += fParamSQ[dim][type][0];
 834   value += fParamSQ[dim][type][1]*z;
 835   value += fParamSQ[dim][type][2]*angle*angle;
 836   value += fParamSQ[dim][type][3]*z/Qmean;
 837   value += fParamSQ[dim][type][4]*angle*angle/Qmean;
 838   value  = TMath::Sqrt(TMath::Abs(value));
 839   return value;
 840
 841
 842 }
 843
 844 Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 845   //
 846   //
 847   //
 848   Float_t value=0;
 849   value += fParamSQPar[dim][type][0];
 850   value += fParamSQPar[dim][type][1]*z;
 851   value += fParamSQPar[dim][type][2]*angle*angle;
 852   value += fParamSQPar[dim][type][3]*z*z;
 853   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
 854   value += fParamSQPar[dim][type][5]*z*angle*angle;
 855   value += fParamSQPar[dim][type][6]*z/Qmean;
 856   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
 857   value  = TMath::Sqrt(TMath::Abs(value));
 858   return value;
 859
 860
 861 }
 862
 863 Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 864   //
 865   //
 866   //
 867   Float_t value=0;
 868   value += fParamSQPar[dim][type][0];
 869   value += fParamSQPar[dim][type][1]*z;
 870   value += fParamSQPar[dim][type][2]*angle*angle;
 871   value += fParamSQPar[dim][type][3]*z*z;
 872   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
 873   value += fParamSQPar[dim][type][5]*z*angle*angle;
 874   value += fParamSQPar[dim][type][6]*z/Qmean;
 875   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
 876   Float_t valueMean = GetError0Par(dim,type,z,angle);
 877   value -= 0.35*0.35*valueMean*valueMean;
 878   value  = TMath::Sqrt(TMath::Abs(value));
 879   return value;
 880
 881
 882 }
 883
 884 Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
 885   //
 886   // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
 887   //
 888   Float_t value=0;
 889   value += fParamRMS0[dim][type][0];
 890   value += fParamRMS0[dim][type][1]*z;
 891   value += fParamRMS0[dim][type][2]*angle*angle;
 892   value  = TMath::Sqrt(TMath::Abs(value));
 893   return value;
 894 }
 895
 896 Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
 897   //
 898   // calculate mean RMS of cluster - z,angle - pad length scalling
 899   //
 900   Float_t value=0;
 901   Float_t length=0.75;
 902   if (type==1) length=1;
 903   if (type==2) length=1.5;
 904   if (type==0){
 905     value += fParamRMS1[dim][0];
 906   }else{
 907     value += fParamRMS1[dim][1];
 908   }
 909   value += fParamRMS1[dim][2]*z;
 910   value += fParamRMS1[dim][3]*angle*angle*length*length;
 911   value  = TMath::Sqrt(TMath::Abs(value));
 912   return value;
 913 }
 914
 915 Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 916   //
 917   // calculate mean RMS of cluster - z,angle, Q dependence
 918   //
 919   Float_t value=0;
 920   value += fParamRMSQ[dim][type][0];
 921   value += fParamRMSQ[dim][type][1]*z;
 922   value += fParamRMSQ[dim][type][2]*angle*angle;
 923   value += fParamRMSQ[dim][type][3]*z/Qmean;
 924   value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
 925   value  = TMath::Sqrt(TMath::Abs(value));
 926   return value;
 927 }
 928
 929 Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 930   //
 931   // calculates RMS of signal shape fluctuation
 932   //
 933   Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
 934   Float_t value  = fRMSSigmaFit[dim][type][0];
 935   value+=  fRMSSigmaFit[dim][type][1]*mean;
 936   return value;
 937 }
 938
 939 Float_t  AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
 940   //
 941   // calculates vallue - sigma distortion contribution
 942   //
 943   Double_t value =0;
 944   //
 945   Float_t rmsMeanQ  = GetRMSQ(dim,type,z,angle,Qmean);
 946   if (rmsL<rmsMeanQ) return value;
 947   //
 948   Float_t rmsSigma  = GetRMSSigma(dim,type,z,angle,Qmean);
 949   //
 950   if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
 951     //1.5 sigma cut on mean
 952     value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
 953   }else{
 954     if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
 955       //3 sigma cut on local
 956       value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
 957     }
 958   }
 959   return TMath::Sqrt(TMath::Abs(value));
 960 }
 961
 962
 963
 964 void AliTPCClusterParam::FitData(TTree * tree){
 965   //
 966   // make fits for error param and shape param
 967   //
 968   FitResol(tree);
 969   FitRMS(tree);
 970
 971 }
 972
 973 void AliTPCClusterParam::FitResol(TTree * tree){
 974   //
 975   SetInstance(this);
 976   for (Int_t idir=0;idir<2; idir++){
 977     for (Int_t itype=0; itype<3; itype++){
 978       Float_t param0[10];
 979       Float_t error0[10];
 980       // model error param
 981       FitResol0(tree, idir, itype,param0,error0);
 982       printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
 983       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
 984       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
 985       for (Int_t ipar=0;ipar<4; ipar++){
 986         fParamS0[idir][itype][ipar] = param0[ipar];
 987         fErrorS0[idir][itype][ipar] = param0[ipar];
 988       }
 989       // error param with parabolic correction
 990       FitResol0Par(tree, idir, itype,param0,error0);
 991       printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
 992       printf("%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5]);
 993       printf("%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5]);
 994       for (Int_t ipar=0;ipar<7; ipar++){
 995         fParamS0Par[idir][itype][ipar] = param0[ipar];
 996         fErrorS0Par[idir][itype][ipar] = param0[ipar];
 997       }
 998       //
 999       FitResolQ(tree, idir, itype,param0,error0);
1000       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1001       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1002       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1003       for (Int_t ipar=0;ipar<6; ipar++){
1004         fParamSQ[idir][itype][ipar] = param0[ipar];
1005         fErrorSQ[idir][itype][ipar] = param0[ipar];
1006       }
1007       //
1008       FitResolQPar(tree, idir, itype,param0,error0);
1009       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
1010       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5],param0[6],param0[7]);
1011       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5],error0[6],error0[7]);
1012       for (Int_t ipar=0;ipar<9; ipar++){
1013         fParamSQPar[idir][itype][ipar] = param0[ipar];
1014         fErrorSQPar[idir][itype][ipar] = param0[ipar];
1015       }
1016     }
1017   }
1018   //
1019   printf("Resol z-scaled\n");
1020   for (Int_t idir=0;idir<2; idir++){
1021     Float_t param0[4];
1022     Float_t error0[4];
1023     FitResol1(tree, idir,param0,error0);
1024     printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
1025     printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1026     printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1027     for (Int_t ipar=0;ipar<4; ipar++){
1028       fParamS1[idir][ipar] = param0[ipar];
1029       fErrorS1[idir][ipar] = param0[ipar];
1030     }
1031   }
1032
1033   for (Int_t idir=0;idir<2; idir++){
1034     printf("\nDirection %d\n",idir);
1035     printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1036     for (Int_t itype=0; itype<3; itype++){
1037       Float_t length=0.75;
1038       if (itype==1) length=1;
1039       if (itype==2) length=1.5;
1040       printf("%d\t%f\t%f\t%f\n", itype,fParamS0[idir][itype][0],fParamS0[idir][itype][1]*TMath::Sqrt(length),fParamS0[idir][itype][2]/TMath::Sqrt(length));
1041     }
1042   }
1043 }
1044
1045
1046
1047 void AliTPCClusterParam::FitRMS(TTree * tree){
1048   //
1049   SetInstance(this);
1050   for (Int_t idir=0;idir<2; idir++){
1051     for (Int_t itype=0; itype<3; itype++){
1052       Float_t param0[6];
1053       Float_t error0[6];
1054       FitRMS0(tree, idir, itype,param0,error0);
1055       printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1056       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1057       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1058       for (Int_t ipar=0;ipar<4; ipar++){
1059         fParamRMS0[idir][itype][ipar] = param0[ipar];
1060         fErrorRMS0[idir][itype][ipar] = param0[ipar];
1061       }
1062       FitRMSQ(tree, idir, itype,param0,error0);
1063       printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1064       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1065       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1066       for (Int_t ipar=0;ipar<6; ipar++){
1067         fParamRMSQ[idir][itype][ipar] = param0[ipar];
1068         fErrorRMSQ[idir][itype][ipar] = param0[ipar];
1069       }
1070     }
1071   }
1072   //
1073   printf("RMS z-scaled\n");
1074   for (Int_t idir=0;idir<2; idir++){
1075     Float_t param0[5];
1076     Float_t error0[5];
1077     FitRMS1(tree, idir,param0,error0);
1078     printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1079     printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1080     printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1081     for (Int_t ipar=0;ipar<5; ipar++){
1082       fParamRMS1[idir][ipar] = param0[ipar];
1083       fErrorRMS1[idir][ipar] = param0[ipar];
1084     }
1085   }
1086
1087   for (Int_t idir=0;idir<2; idir++){
1088     printf("\nDirection %d\n",idir);
1089     printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1090     for (Int_t itype=0; itype<3; itype++){
1091       Float_t length=0.75;
1092       if (itype==1) length=1;
1093       if (itype==2) length=1.5;
1094       if (itype==0) printf("%d\t%f\t\t\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length);
1095       if (itype>0) printf("%d\t\t\t%f\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length);
1096     }
1097   }
1098   //
1099   // Fit RMS sigma
1100   //
1101   printf("RMS fluctuation  parameterization \n");
1102   for (Int_t idir=0;idir<2; idir++){
1103     for (Int_t itype=0; itype<3; itype++){
1104       Float_t param0[5];
1105       Float_t error0[5];
1106       FitRMSSigma(tree, idir,itype,param0,error0);
1107       printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1108       for (Int_t ipar=0;ipar<2; ipar++){
1109         fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1110       }
1111     }
1112   }
1113   //
1114   // store systematic error end RMS fluctuation parameterization
1115   //
1116   TH1F hratio("hratio","hratio",100,-0.1,0.1);
1117   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1118   fErrorRMSSys[0] = hratio.GetRMS();
1119   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1120   fErrorRMSSys[1] = hratio.GetRMS();
1121   TH1F hratioR("hratioR","hratioR",100,0,0.2);
1122   tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1123   fRMSSigmaRatio[0][0]=hratioR.GetMean();
1124   fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1125   tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1126   fRMSSigmaRatio[1][0]=hratioR.GetMean();
1127   fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1128 }
1129
1130 void AliTPCClusterParam::Test(TTree * tree, const char *output){
1131   //
1132   // Draw standard quality histograms to output file
1133   //
1134   SetInstance(this);
1135   TFile f(output,"recreate");
1136   f.cd();
1137   //
1138   // 1D histograms - resolution
1139   //
1140   for (Int_t idim=0; idim<2; idim++){
1141     for (Int_t ipad=0; ipad<3; ipad++){
1142       char hname1[300];
1143       char hcut1[300];
1144       char hexp1[300];
1145       //
1146       sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1147       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1148       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1149       TH1F  his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1150       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1151       tree->Draw(hexp1,hcut1,"");
1152       his1DRel0.Write();
1153       //
1154       sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1155       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1156       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1157       TH1F  his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1158       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1159       tree->Draw(hexp1,hcut1,"");
1160       his1DRel0Par.Write();
1161       //
1162     }
1163   }
1164   //
1165   // 2D histograms - resolution
1166   //
1167   for (Int_t idim=0; idim<2; idim++){
1168     for (Int_t ipad=0; ipad<3; ipad++){
1169       char hname1[300];
1170       char hcut1[300];
1171       char hexp1[300];
1172       //
1173       sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1174       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1175       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1176       TProfile2D  profDRel0(hname1, hname1, 6,0,250,6,0,1);
1177       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1178       tree->Draw(hexp1,hcut1,"");
1179       profDRel0.Write();
1180       //
1181       sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1182       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1183       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1184       TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1185       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1186       tree->Draw(hexp1,hcut1,"");
1187       profDRel0Par.Write();
1188       //
1189     }
1190   }
1191 }
1192
1193
1194
1195 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1196   //
1197   // Print param Information
1198   //
1199
1200   //
1201   // Error parameterization
1202   //
1203   printf("\nResolution Scaled factors\n");
1204   printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1205   printf("Y\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[0][0])),TMath::Sqrt(TMath::Abs(fParamS1[0][1])),
1206          TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
1207   for (Int_t ipad=0; ipad<3; ipad++){
1208     Float_t length=0.75;
1209     if (ipad==1) length=1;
1210     if (ipad==2) length=1.5;
1211     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1212            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1213            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1214            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1215            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
1216   }
1217   for (Int_t ipad=0; ipad<3; ipad++){
1218     Float_t length=0.75;
1219     if (ipad==1) length=1;
1220     if (ipad==2) length=1.5;
1221     printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1222            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1223            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1224            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1225            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
1226   }
1227   printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1228          TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1229
1230   for (Int_t ipad=0; ipad<3; ipad++){
1231     Float_t length=0.75;
1232     if (ipad==1) length=1;
1233     if (ipad==2) length=1.5;
1234     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1235            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1236            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1237            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1238            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
1239   }
1240   for (Int_t ipad=0; ipad<3; ipad++){
1241     Float_t length=0.75;
1242     if (ipad==1) length=1;
1243     if (ipad==2) length=1.5;
1244     printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1245            TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1246            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1247            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1248            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
1249   }
1250
1251   //
1252   // RMS scaling
1253   //
1254   printf("\n");
1255   printf("\nRMS Scaled factors\n");
1256   printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1257   printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1258          TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1259          TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1260          TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1261          TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1262          TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
1263   for (Int_t ipad=0; ipad<3; ipad++){
1264     Float_t length=0.75;
1265     if (ipad==1) length=1;
1266     if (ipad==2) length=1.5;
1267     if (ipad==0){
1268       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1269              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1270              0.,
1271              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1272              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1273              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1274     }else{
1275       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1276              0.,
1277              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1278              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1279              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1280              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1281     }
1282   }
1283   printf("\n");
1284   printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1285          TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1286          TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1287          TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1288          TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1289          TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
1290   for (Int_t ipad=0; ipad<3; ipad++){
1291     Float_t length=0.75;
1292     if (ipad==1) length=1;
1293     if (ipad==2) length=1.5;
1294     if (ipad==0){
1295       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1296              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1297              0.,
1298              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1299              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1300              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1301     }else{
1302       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1303              0.,
1304              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1305              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1306              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1307              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1308     }
1309   }
1310 }
1311
1312
1313
1314
1315
1316 Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1317   // get Q normalization
1318   // type - 0 Qtot 1 Qmax
1319   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1320   //
1321   //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
1322
1323   if (fQNorm==0) return 0;
1324   TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1325   if (!norm) return 0;
1326   TVectorD &no  = *norm;
1327   Float_t   res =
1328     no[0]+
1329     no[1]*dr+
1330     no[2]*ty+
1331     no[3]*tz+
1332     no[4]*dr*ty+
1333     no[5]*dr*tz+
1334     no[6]*ty*tz+
1335     no[7]*dr*dr+
1336     no[8]*ty*ty+
1337     no[9]*tz*tz;
1338   res/=no[0];
1339   return res;
1340 }
1341
1342
1343
1344 void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
1345   //
1346   // set normalization
1347   //
1348   // type - 0 Qtot 1 Qmax
1349   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1350   //
1351
1352   if (fQNorm==0) fQNorm = new TObjArray(6);
1353   fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1354 }
1355
1356 Float_t AliTPCClusterParam::QnormPos(Int_t ipad,Bool_t isMax, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm, Float_t qt){
1357   //
1358   // Make Q normalization as function of following parameters
1359   // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion
1360   // 1 - dp   - relative pad position
1361   // 2 - dt   - relative time position
1362   // 3 - di   - drift length (norm to 1);
1363   // 4 - dq0  - Tot/Max charge
1364   // 5 - dq1  - Max/Tot charge
1365   // 6 - sy   - sigma y - shape
1366   // 7 - sz   - sigma z - shape
1367   //
1368
1369   //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
1370   // Following variable used - correspondance to the our variable conventions
1371   //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1372   Double_t dp = ((pad-int(pad)-0.5)*2.);
1373   //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1374   Double_t dt = ((time-int(time)-0.5)*2.);
1375   //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1376   Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1377   //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1378   Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1379   //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1380   Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1381   //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1382   Double_t sy  = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1383   //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1384   Double_t sz  = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1385   //
1386   //
1387   //
1388   TVectorD * pvec = 0;
1389   if (isMax){
1390     pvec = fPosQMnorm[ipad];
1391   }else{
1392     pvec = fPosQTnorm[ipad];
1393   }
1394   TVectorD &param = *pvec;
1395   //
1396   // Eval part  - in correspondance with fit part from debug streamer
1397   //
1398   Double_t result=param[0];
1399   Int_t index =1;
1400   //
1401   result+=dp*param[index++];                               //1
1402   result+=dt*param[index++];                               //2
1403   result+=dp*dp*param[index++];                             //3
1404   result+=dt*dt*param[index++];                             //4
1405   result+=dt*dt*dt*param[index++];                             //5
1406   result+=dp*dt*param[index++];                            //6
1407   result+=dp*dt*dt*param[index++];                          //7
1408   result+=(dq0)*param[index++];                            //8
1409   result+=(dq1)*param[index++];                            //9
1410   //
1411   //
1412   result+=dp*dp*(di)*param[index++];                        //10
1413   result+=dt*dt*(di)*param[index++];                        //11
1414   result+=dp*dp*sy*param[index++];                          //12
1415   result+=dt*sz*param[index++];                          //13
1416   result+=dt*dt*sz*param[index++];                          //14
1417   result+=dt*dt*dt*sz*param[index++];                          //15
1418   //
1419   result+=dp*dp*1*sy*sz*param[index++];                     //16
1420   result+=dt*sy*sz*param[index++];                       //17
1421   result+=dt*dt*sy*sz*param[index++];                       //18
1422   result+=dt*dt*dt*sy*sz*param[index++];                       //19
1423   //
1424   result+=dp*dp*(dq0)*param[index++];                       //20
1425   result+=dt*1*(dq0)*param[index++];                       //21
1426   result+=dt*dt*(dq0)*param[index++];                       //22
1427   result+=dt*dt*dt*(dq0)*param[index++];                       //23
1428   //
1429   result+=dp*dp*(dq1)*param[index++];                       //24
1430   result+=dt*(dq1)*param[index++];                       //25
1431   result+=dt*dt*(dq1)*param[index++];                       //26
1432   result+=dt*dt*dt*(dq1)*param[index++];                       //27
1433
1434   if (result<0.75) result=0.75;
1435   if (result>1.25) result=1.25;
1436
1437   return result;
1438
1439 }
1440
1441
1442
1443
1444
1445 Float_t AliTPCClusterParam::PosCorrection(Int_t type, Int_t ipad,  Float_t pad, Float_t time, Float_t z, Float_t /*sy2*/, Float_t /*sz2*/, Float_t /*qm*/){
1446
1447   //
1448   // Make postion correction
1449   // type - 0 - y correction
1450   //        1 - z correction
1451   // ipad - 0, 1, 2 - short, medium long pads
1452   // pad  - float pad number
1453   // time - float time bin number
1454   //    z - z of the cluster
1455
1456   //
1457   //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
1458   //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
1459   //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
1460   //chainres->SetAlias("st","(sin(dt)-dt)");
1461   //
1462   //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1463
1464   //
1465   // Derived variables
1466   //
1467   Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
1468   Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
1469   Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
1470   Double_t st = (TMath::Sin(dt)-dt);
1471   //
1472   Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.)));
1473   //
1474   //
1475   //
1476   TVectorD * pvec = 0;
1477   if (type==0){
1478     pvec = fPosYcor[ipad];
1479   }else{
1480     pvec = fPosZcor[ipad];
1481   }
1482   TVectorD &param = *pvec;
1483   //
1484   Double_t result=0;
1485   Int_t index =1;
1486
1487   if (type==0){
1488     // y corr
1489     result+=(dp)*param[index++];             //1
1490     result+=(dp)*di*param[index++];          //2
1491     //
1492     result+=(sp)*param[index++];             //3
1493     result+=(sp)*di*param[index++];          //4
1494   }
1495   if (type==1){
1496     result+=(dt)*param[index++];             //1
1497     result+=(dt)*di*param[index++];          //2
1498     //
1499     result+=(st)*param[index++];             //3
1500     result+=(st)*di*param[index++];          //4
1501   }
1502   if (TMath::Abs(result)>0.05) return 0;
1503   return result;
1504 }
1505
1506
1507
1508 Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
1509   //
1510   // 2 D gaus convoluted with angular effect
1511   // See in mathematica:
1512   //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
1513   //
1514   //TF1 f1("f1","AliTPCClusterParam::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
1515   //TF2 f2("f2","AliTPCClusterParam::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
1516   //
1517   const Float_t kEpsilon = 0.0001;
1518   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
1519     // small angular effect
1520     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
1521     return val;
1522   }
1523   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
1524   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
1525   //
1526   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);
1527   Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
1528   Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
1529   Double_t norm = 1./TMath::Sqrt(sigma2);
1530   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
1531   Double_t val  = norm*exp0*(erf0+erf1);
1532   return val;
1533 }
1534
1535
1536 Double_t  AliTPCClusterParam::GaussConvolutionTail(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
1537   //
1538   // 2 D gaus convoluted with angular effect and exponential tail in z-direction
1539   // tail integrated numerically
1540   // Integral normalized to one
1541   // Mean at 0
1542   //
1543   // TF1 f1t("f1t","AliTPCClusterParam::GaussConvolutionTail(0,x,0,0,0.5,0.5,0.9)",-5,5)
1544   Double_t sum =1, mean=0;
1545   // the COG of exponent
1546   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1547     mean+=iexp*TMath::Exp(-iexp/tau);
1548     sum +=TMath::Exp(-iexp/tau);
1549   }
1550   mean/=sum;
1551   //
1552   sum = 1;
1553   Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
1554   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1555     val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
1556     sum+=TMath::Exp(-iexp/tau);
1557   }
1558   return val/sum;
1559 }
1560
1561 Double_t  AliTPCClusterParam::GaussConvolutionGamma4(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
1562   //
1563   // 2 D gaus convoluted with angular effect and exponential tail in z-direction
1564   // tail integrated numerically
1565   // Integral normalized to one
1566   // Mean at 0
1567   //
1568   // TF1 f1g4("f1g4","AliTPCClusterParam::GaussConvolutionGamma4(0,x,0,0,0.5,0.2,1.6)",-5,5)
1569   // TF2 f2g4("f2g4","AliTPCClusterParam::GaussConvolutionGamma4(y,x,0,0,0.5,0.2,1.6)",-5,5,-5,5)
1570   Double_t sum =0, mean=0;
1571   // the COG of G4
1572   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1573     Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1574     mean+=iexp*g4;
1575     sum +=g4;
1576   }
1577   mean/=sum;
1578   //
1579   sum = 0;
1580   Double_t val = 0;
1581   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1582     Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1583     val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
1584     sum+=g4;
1585   }
1586   return val/sum;
1587 }
1588
1589 Double_t  AliTPCClusterParam::QmaxCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0, Float_t tau){
1590   //
1591   //
1592   // cpad      - pad (y) coordinate
1593   // ctime     - time(z) coordinate
1594   // ky        - dy/dx
1595   // kz        - dz/dx
1596   // rmsy0     - RF width in pad units
1597   // rmsz0     - RF width in pad units
1598   // tau       - correction for the assymetric TRF shape - decay length
1599
1600   // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1601   //
1602   // Gaus width sy and sz is determined by RF width and diffusion
1603   // Integral of Q is equal 1
1604   // Q max is calculated at position cpad, ctime
1605   // Example function:
1606   //  TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000)
1607   //
1608   AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters();
1609   Double_t padLength= param->GetPadPitchLength(sector,row);
1610   Double_t padWidth = param->GetPadPitchWidth(sector);
1611   // diffusion in pad units
1612   Double_t diffT=TMath::Sqrt(ctime*param->GetZWidth())*param->GetDiffT()/padWidth;
1613   Double_t diffL=TMath::Sqrt(ctime*param->GetZWidth())*param->GetDiffL()/padWidth;
1614   //
1615   // transform angular effect to pad units
1616   Double_t pky   = ky*padLength/padWidth;
1617   Double_t pkz   = kz*padLength/padWidth;
1618   // position in pad unit
1619   Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1620   Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1621   //
1622   //
1623   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1624   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1625   return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
1626 }
1627
1628 Double_t  AliTPCClusterParam::QtotCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0,  Float_t qtot, Float_t thr, Float_t tau){
1629   //
1630   //
1631   // cpad      - pad (y) coordinate
1632   // ctime     - time(z) coordinate
1633   // ky        - dy/dx
1634   // kz        - dz/dx
1635   // rmsy0     - RF width in pad units
1636   // rmsz0     - RF width in pad units
1637   // qtot      - the sum of signal in cluster - without thr correction
1638   // thr       - threshold
1639   // tau       - correction for the assymetric TRF shape - decay length
1640
1641   // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1642   //
1643   // Gaus width sy and sz is determined by RF width and diffusion
1644   // Integral of Q is equal 1
1645   // Q max is calculated at position cpad, ctime
1646   //
1647   //
1648   //
1649   AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters();
1650   Double_t padLength= param->GetPadPitchLength(sector,row);
1651   Double_t padWidth = param->GetPadPitchWidth(sector);
1652   //
1653   // diffusion in pad units
1654   Double_t diffT=TMath::Sqrt(ctime*param->GetZWidth())*param->GetDiffT()/padWidth;
1655   Double_t diffL=TMath::Sqrt(ctime*param->GetZWidth())*param->GetDiffL()/padWidth;
1656   //
1657   // transform angular effect to pad units
1658   Double_t pky   = ky*padLength/padWidth;
1659   Double_t pkz   = kz*padLength/padWidth;
1660   // position in pad unit
1661   //
1662   Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1663   Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1664   //
1665   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1666   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1667   //
1668   //
1669   //
1670   Double_t sumAll=0,sumThr=0;
1671   //
1672   Double_t corr =1;
1673   Double_t qnorm=qtot;
1674   for (Float_t iy=-2;iy<=2;iy++)
1675     for (Float_t iz=-2;iz<=2;iz++){
1676       Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);
1677       Double_t qlocal =qnorm*val;
1678       if (TMath::Abs(iy)<2&&TMath::Abs(iz)<2){
1679         sumThr+=qlocal;   // using Virtual charge in cluster finder
1680       }
1681       else{
1682         if (qlocal>thr) sumThr+=qlocal;
1683       }
1684       sumAll+=qlocal;
1685     }
1686   if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
1687   //
1688   return corr;
1689 }
1690
1691
1692
1693
1694
1695
1696