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