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