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