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