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 {
 149   //
 150   // Default constructor
 151   //
 152 }
 153
 154 AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
 155   TObject(param),
 156   fRatio(0),
 157   fQNorm(0)
 158 {
 159   //
 160   // copy constructor
 161   //
 162   memcpy(this, &param,sizeof(AliTPCClusterParam));
 163   if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
 164 }
 165
 166 AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
 167   //
 168   // Assignment operator
 169   //
 170   if (this != &param) {
 171     memcpy(this, &param,sizeof(AliTPCClusterParam));
 172     if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
 173   }
 174   return *this;
 175 }
 176
 177
 178 AliTPCClusterParam::~AliTPCClusterParam(){
 179   //
 180   // destructor
 181   //
 182   if (fQNorm) fQNorm->Delete();
 183   delete fQNorm;
 184 }
 185
 186
 187 void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 188   //
 189   // Fit z - angular dependence of resolution
 190   //
 191   // Int_t dim=0, type=0;
 192   char varVal[100];
 193   sprintf(varVal,"Resol:AngleM:Zm");
 194   char varErr[100];
 195   sprintf(varErr,"Sigma:AngleS:Zs");
 196   char varCut[100];
 197   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 198   //
 199   Int_t entries = tree->Draw(varVal,varCut);
 200   Float_t px[10000], py[10000], pz[10000];
 201   Float_t ex[10000], ey[10000], ez[10000];
 202   //
 203   tree->Draw(varErr,varCut);
 204   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 205     ex[ipoint]= tree->GetV3()[ipoint];
 206     ey[ipoint]= tree->GetV2()[ipoint];
 207     ez[ipoint]= tree->GetV1()[ipoint];
 208   }
 209   tree->Draw(varVal,varCut);
 210   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 211     px[ipoint]= tree->GetV3()[ipoint];
 212     py[ipoint]= tree->GetV2()[ipoint];
 213     pz[ipoint]= tree->GetV1()[ipoint];
 214   }
 215
 216   //
 217   TLinearFitter fitter(3,"hyp2");
 218   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 219     Float_t val = pz[ipoint]*pz[ipoint];
 220     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 221     Double_t x[2];
 222     x[0] = px[ipoint];
 223     x[1] = py[ipoint]*py[ipoint];
 224     fitter.AddPoint(x,val,err);
 225   }
 226   fitter.Eval();
 227   TVectorD param(3);
 228   fitter.GetParameters(param);
 229   param0[0] = param[0];
 230   param0[1] = param[1];
 231   param0[2] = param[2];
 232   Float_t chi2 =  fitter.GetChisquare()/entries;
 233   param0[3] = chi2;
 234   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 235   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 236   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 237 }
 238
 239
 240 void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 241   //
 242   // Fit z - angular dependence of resolution
 243   //
 244   // Int_t dim=0, type=0;
 245   char varVal[100];
 246   sprintf(varVal,"Resol:AngleM:Zm");
 247   char varErr[100];
 248   sprintf(varErr,"Sigma:AngleS:Zs");
 249   char varCut[100];
 250   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 251   //
 252   Int_t entries = tree->Draw(varVal,varCut);
 253   Float_t px[10000], py[10000], pz[10000];
 254   Float_t ex[10000], ey[10000], ez[10000];
 255   //
 256   tree->Draw(varErr,varCut);
 257   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 258     ex[ipoint]= tree->GetV3()[ipoint];
 259     ey[ipoint]= tree->GetV2()[ipoint];
 260     ez[ipoint]= tree->GetV1()[ipoint];
 261   }
 262   tree->Draw(varVal,varCut);
 263   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 264     px[ipoint]= tree->GetV3()[ipoint];
 265     py[ipoint]= tree->GetV2()[ipoint];
 266     pz[ipoint]= tree->GetV1()[ipoint];
 267   }
 268
 269   //
 270   TLinearFitter fitter(6,"hyp5");
 271   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 272     Float_t val = pz[ipoint]*pz[ipoint];
 273     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 274     Double_t x[6];
 275     x[0] = px[ipoint];
 276     x[1] = py[ipoint]*py[ipoint];
 277     x[2] = x[0]*x[0];
 278     x[3] = x[1]*x[1];
 279     x[4] = x[0]*x[1];
 280     fitter.AddPoint(x,val,err);
 281   }
 282   fitter.Eval();
 283   TVectorD param(6);
 284   fitter.GetParameters(param);
 285   param0[0] = param[0];
 286   param0[1] = param[1];
 287   param0[2] = param[2];
 288   param0[3] = param[3];
 289   param0[4] = param[4];
 290   param0[5] = param[5];
 291   Float_t chi2 =  fitter.GetChisquare()/entries;
 292   param0[6] = chi2;
 293   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 294   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 295   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 296   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 297   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 298   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
 299 }
 300
 301
 302
 303
 304
 305 void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
 306   //
 307   // Fit z - angular dependence of resolution - pad length scaling
 308   //
 309   // Int_t dim=0, type=0;
 310   char varVal[100];
 311   sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
 312   char varErr[100];
 313   sprintf(varErr,"Sigma:AngleS:Zs");
 314   char varCut[100];
 315   sprintf(varCut,"Dim==%d&&QMean<0",dim);
 316   //
 317   Int_t entries = tree->Draw(varVal,varCut);
 318   Float_t px[10000], py[10000], pz[10000];
 319   Float_t ex[10000], ey[10000], ez[10000];
 320   //
 321   tree->Draw(varErr,varCut);
 322   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 323     ex[ipoint]= tree->GetV3()[ipoint];
 324     ey[ipoint]= tree->GetV2()[ipoint];
 325     ez[ipoint]= tree->GetV1()[ipoint];
 326   }
 327   tree->Draw(varVal,varCut);
 328   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 329     px[ipoint]= tree->GetV3()[ipoint];
 330     py[ipoint]= tree->GetV2()[ipoint];
 331     pz[ipoint]= tree->GetV1()[ipoint];
 332   }
 333
 334   //
 335   TLinearFitter fitter(3,"hyp2");
 336   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 337     Float_t val = pz[ipoint]*pz[ipoint];
 338     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 339     Double_t x[2];
 340     x[0] = px[ipoint];
 341     x[1] = py[ipoint]*py[ipoint];
 342     fitter.AddPoint(x,val,err);
 343   }
 344   fitter.Eval();
 345   TVectorD param(3);
 346   fitter.GetParameters(param);
 347   param0[0] = param[0];
 348   param0[1] = param[1];
 349   param0[2] = param[2];
 350   Float_t chi2 =  fitter.GetChisquare()/entries;
 351   param0[3] = chi2;
 352   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 353   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 354   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 355 }
 356
 357 void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 358   //
 359   // Fit z - angular dependence of resolution - Q scaling
 360   //
 361   // Int_t dim=0, type=0;
 362   char varVal[100];
 363   sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
 364   char varVal0[100];
 365   sprintf(varVal0,"Resol:AngleM:Zm");
 366   //
 367   char varErr[100];
 368   sprintf(varErr,"Sigma:AngleS:Zs");
 369   char varCut[100];
 370   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 371   //
 372   Int_t entries = tree->Draw(varVal,varCut);
 373   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 374   Float_t ex[20000], ey[20000], ez[20000];
 375   //
 376   tree->Draw(varErr,varCut);
 377   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 378     ex[ipoint]= tree->GetV3()[ipoint];
 379     ey[ipoint]= tree->GetV2()[ipoint];
 380     ez[ipoint]= tree->GetV1()[ipoint];
 381   }
 382   tree->Draw(varVal,varCut);
 383   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 384     px[ipoint]= tree->GetV3()[ipoint];
 385     py[ipoint]= tree->GetV2()[ipoint];
 386     pz[ipoint]= tree->GetV1()[ipoint];
 387   }
 388   tree->Draw(varVal0,varCut);
 389   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 390     pu[ipoint]= tree->GetV3()[ipoint];
 391     pt[ipoint]= tree->GetV2()[ipoint];
 392   }
 393
 394   //
 395   TLinearFitter fitter(5,"hyp4");
 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[4];
 400     x[0] = pu[ipoint];
 401     x[1] = pt[ipoint]*pt[ipoint];
 402     x[2] = px[ipoint];
 403     x[3] = py[ipoint]*py[ipoint];
 404     fitter.AddPoint(x,val,err);
 405   }
 406
 407   fitter.Eval();
 408   TVectorD param(5);
 409   fitter.GetParameters(param);
 410   param0[0] = param[0];
 411   param0[1] = param[1];
 412   param0[2] = param[2];
 413   param0[3] = param[3];
 414   param0[4] = param[4];
 415   Float_t chi2 =  fitter.GetChisquare()/entries;
 416   param0[5] = chi2;
 417   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 418   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 419   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 420   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 421   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 422 }
 423
 424 void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 425   //
 426   // Fit z - angular dependence of resolution - Q scaling  - parabolic correction
 427   //
 428   // Int_t dim=0, type=0;
 429   char varVal[100];
 430   sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
 431   char varVal0[100];
 432   sprintf(varVal0,"Resol:AngleM:Zm");
 433   //
 434   char varErr[100];
 435   sprintf(varErr,"Sigma:AngleS:Zs");
 436   char varCut[100];
 437   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 438   //
 439   Int_t entries = tree->Draw(varVal,varCut);
 440   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 441   Float_t ex[20000], ey[20000], ez[20000];
 442   //
 443   tree->Draw(varErr,varCut);
 444   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 445     ex[ipoint]= tree->GetV3()[ipoint];
 446     ey[ipoint]= tree->GetV2()[ipoint];
 447     ez[ipoint]= tree->GetV1()[ipoint];
 448   }
 449   tree->Draw(varVal,varCut);
 450   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 451     px[ipoint]= tree->GetV3()[ipoint];
 452     py[ipoint]= tree->GetV2()[ipoint];
 453     pz[ipoint]= tree->GetV1()[ipoint];
 454   }
 455   tree->Draw(varVal0,varCut);
 456   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 457     pu[ipoint]= tree->GetV3()[ipoint];
 458     pt[ipoint]= tree->GetV2()[ipoint];
 459   }
 460
 461   //
 462   TLinearFitter fitter(8,"hyp7");
 463   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 464     Float_t val = pz[ipoint]*pz[ipoint];
 465     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 466     Double_t x[7];
 467     x[0] = pu[ipoint];
 468     x[1] = pt[ipoint]*pt[ipoint];
 469     x[2] = x[0]*x[0];
 470     x[3] = x[1]*x[1];
 471     x[4] = x[0]*x[1];
 472     x[5] = px[ipoint];
 473     x[6] = py[ipoint]*py[ipoint];
 474     //
 475     fitter.AddPoint(x,val,err);
 476   }
 477
 478   fitter.Eval();
 479   TVectorD param(8);
 480   fitter.GetParameters(param);
 481   param0[0] = param[0];
 482   param0[1] = param[1];
 483   param0[2] = param[2];
 484   param0[3] = param[3];
 485   param0[4] = param[4];
 486   param0[5] = param[5];
 487   param0[6] = param[6];
 488   param0[7] = param[7];
 489
 490   Float_t chi2 =  fitter.GetChisquare()/entries;
 491   param0[8] = chi2;
 492   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 493   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 494   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 495   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 496   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 497   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
 498   error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
 499   error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
 500 }
 501
 502
 503
 504 void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 505   //
 506   // Fit z - angular dependence of resolution
 507   //
 508   // Int_t dim=0, type=0;
 509   char varVal[100];
 510   sprintf(varVal,"RMSm:AngleM:Zm");
 511   char varErr[100];
 512   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
 513   char varCut[100];
 514   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 515   //
 516   Int_t entries = tree->Draw(varVal,varCut);
 517   Float_t px[10000], py[10000], pz[10000];
 518   Float_t ex[10000], ey[10000], ez[10000];
 519   //
 520   tree->Draw(varErr,varCut);
 521   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 522     ex[ipoint]= tree->GetV3()[ipoint];
 523     ey[ipoint]= tree->GetV2()[ipoint];
 524     ez[ipoint]= tree->GetV1()[ipoint];
 525   }
 526   tree->Draw(varVal,varCut);
 527   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 528     px[ipoint]= tree->GetV3()[ipoint];
 529     py[ipoint]= tree->GetV2()[ipoint];
 530     pz[ipoint]= tree->GetV1()[ipoint];
 531   }
 532
 533   //
 534   TLinearFitter fitter(3,"hyp2");
 535   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 536     Float_t val = pz[ipoint]*pz[ipoint];
 537     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 538     Double_t x[2];
 539     x[0] = px[ipoint];
 540     x[1] = py[ipoint]*py[ipoint];
 541     fitter.AddPoint(x,val,err);
 542   }
 543   fitter.Eval();
 544   TVectorD param(3);
 545   fitter.GetParameters(param);
 546   param0[0] = param[0];
 547   param0[1] = param[1];
 548   param0[2] = param[2];
 549   Float_t chi2 =  fitter.GetChisquare()/entries;
 550   param0[3] = chi2;
 551   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 552   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 553   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 554 }
 555
 556 void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
 557   //
 558   // Fit z - angular dependence of resolution - pad length scaling
 559   //
 560   // Int_t dim=0, type=0;
 561   char varVal[100];
 562   sprintf(varVal,"RMSm:AngleM*Length:Zm");
 563   char varErr[100];
 564   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
 565   char varCut[100];
 566   sprintf(varCut,"Dim==%d&&QMean<0",dim);
 567   //
 568   Int_t entries = tree->Draw(varVal,varCut);
 569   Float_t px[10000], py[10000], pz[10000];
 570   Float_t type[10000], ey[10000], ez[10000];
 571   //
 572   tree->Draw(varErr,varCut);
 573   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 574     type[ipoint] = tree->GetV3()[ipoint];
 575     ey[ipoint]   = tree->GetV2()[ipoint];
 576     ez[ipoint]   = tree->GetV1()[ipoint];
 577   }
 578   tree->Draw(varVal,varCut);
 579   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 580     px[ipoint]= tree->GetV3()[ipoint];
 581     py[ipoint]= tree->GetV2()[ipoint];
 582     pz[ipoint]= tree->GetV1()[ipoint];
 583   }
 584
 585   //
 586   TLinearFitter fitter(4,"hyp3");
 587   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 588     Float_t val = pz[ipoint]*pz[ipoint];
 589     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 590     Double_t x[3];
 591     x[0] = (type[ipoint]<0.5)? 0.:1.;
 592     x[1] = px[ipoint];
 593     x[2] = py[ipoint]*py[ipoint];
 594     fitter.AddPoint(x,val,err);
 595   }
 596   fitter.Eval();
 597   TVectorD param(4);
 598   fitter.GetParameters(param);
 599   param0[0] = param[0];
 600   param0[1] = param[0]+param[1];
 601   param0[2] = param[2];
 602   param0[3] = param[3];
 603   Float_t chi2 =  fitter.GetChisquare()/entries;
 604   param0[4] = chi2;
 605   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 606   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 607   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 608   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 609 }
 610
 611 void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
 612   //
 613   // Fit z - angular dependence of resolution - Q scaling
 614   //
 615   // Int_t dim=0, type=0;
 616   char varVal[100];
 617   sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
 618   char varVal0[100];
 619   sprintf(varVal0,"RMSm:AngleM:Zm");
 620   //
 621   char varErr[100];
 622   sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
 623   char varCut[100];
 624   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
 625   //
 626   Int_t entries = tree->Draw(varVal,varCut);
 627   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
 628   Float_t ex[20000], ey[20000], ez[20000];
 629   //
 630   tree->Draw(varErr,varCut);
 631   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 632     ex[ipoint]= tree->GetV3()[ipoint];
 633     ey[ipoint]= tree->GetV2()[ipoint];
 634     ez[ipoint]= tree->GetV1()[ipoint];
 635   }
 636   tree->Draw(varVal,varCut);
 637   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 638     px[ipoint]= tree->GetV3()[ipoint];
 639     py[ipoint]= tree->GetV2()[ipoint];
 640     pz[ipoint]= tree->GetV1()[ipoint];
 641   }
 642   tree->Draw(varVal0,varCut);
 643   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 644     pu[ipoint]= tree->GetV3()[ipoint];
 645     pt[ipoint]= tree->GetV2()[ipoint];
 646   }
 647
 648   //
 649   TLinearFitter fitter(5,"hyp4");
 650   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 651     Float_t val = pz[ipoint]*pz[ipoint];
 652     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
 653     Double_t x[4];
 654     x[0] = pu[ipoint];
 655     x[1] = pt[ipoint]*pt[ipoint];
 656     x[2] = px[ipoint];
 657     x[3] = py[ipoint]*py[ipoint];
 658     fitter.AddPoint(x,val,err);
 659   }
 660
 661   fitter.Eval();
 662   TVectorD param(5);
 663   fitter.GetParameters(param);
 664   param0[0] = param[0];
 665   param0[1] = param[1];
 666   param0[2] = param[2];
 667   param0[3] = param[3];
 668   param0[4] = param[4];
 669   Float_t chi2 =  fitter.GetChisquare()/entries;
 670   param0[5] = chi2;
 671   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
 672   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
 673   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
 674   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
 675   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
 676 }
 677
 678
 679 void AliTPCClusterParam::FitRMSSigma(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,"RMSs:RMSm");
 686   //
 687   char varCut[100];
 688   sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
 689   //
 690   Int_t entries = tree->Draw(varVal,varCut);
 691   Float_t px[20000], py[20000];
 692   //
 693   tree->Draw(varVal,varCut);
 694   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 695     px[ipoint]= tree->GetV2()[ipoint];
 696     py[ipoint]= tree->GetV1()[ipoint];
 697   }
 698   TLinearFitter fitter(2,"pol1");
 699   for (Int_t ipoint=0; ipoint<entries; ipoint++){
 700     Float_t val = py[ipoint];
 701     Float_t err = fRatio*px[ipoint];
 702     Double_t x[4];
 703     x[0] = px[ipoint];
 704     if (err>0) fitter.AddPoint(x,val,err);
 705   }
 706   fitter.Eval();
 707   param0[0]= fitter.GetParameter(0);
 708   param0[1]= fitter.GetParameter(1);
 709 }
 710
 711
 712
 713 Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
 714   //
 715   //
 716   //
 717   Float_t value=0;
 718   value += fParamS0[dim][type][0];
 719   value += fParamS0[dim][type][1]*z;
 720   value += fParamS0[dim][type][2]*angle*angle;
 721   value  = TMath::Sqrt(TMath::Abs(value));
 722   return value;
 723 }
 724
 725
 726 Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
 727   //
 728   //
 729   //
 730   Float_t value=0;
 731   value += fParamS0Par[dim][type][0];
 732   value += fParamS0Par[dim][type][1]*z;
 733   value += fParamS0Par[dim][type][2]*angle*angle;
 734   value += fParamS0Par[dim][type][3]*z*z;
 735   value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
 736   value += fParamS0Par[dim][type][5]*z*angle*angle;
 737   value  = TMath::Sqrt(TMath::Abs(value));
 738   return value;
 739 }
 740
 741
 742
 743 Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
 744   //
 745   //
 746   //
 747   Float_t value=0;
 748   Float_t length=0.75;
 749   if (type==1) length=1;
 750   if (type==2) length=1.5;
 751   value += fParamS1[dim][0];
 752   value += fParamS1[dim][1]*z/length;
 753   value += fParamS1[dim][2]*angle*angle*length;
 754   value  = TMath::Sqrt(TMath::Abs(value));
 755   return value;
 756 }
 757
 758 Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 759   //
 760   //
 761   //
 762   Float_t value=0;
 763   value += fParamSQ[dim][type][0];
 764   value += fParamSQ[dim][type][1]*z;
 765   value += fParamSQ[dim][type][2]*angle*angle;
 766   value += fParamSQ[dim][type][3]*z/Qmean;
 767   value += fParamSQ[dim][type][4]*angle*angle/Qmean;
 768   value  = TMath::Sqrt(TMath::Abs(value));
 769   return value;
 770
 771
 772 }
 773
 774 Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 775   //
 776   //
 777   //
 778   Float_t value=0;
 779   value += fParamSQPar[dim][type][0];
 780   value += fParamSQPar[dim][type][1]*z;
 781   value += fParamSQPar[dim][type][2]*angle*angle;
 782   value += fParamSQPar[dim][type][3]*z*z;
 783   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
 784   value += fParamSQPar[dim][type][5]*z*angle*angle;
 785   value += fParamSQPar[dim][type][6]*z/Qmean;
 786   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
 787   value  = TMath::Sqrt(TMath::Abs(value));
 788   return value;
 789
 790
 791 }
 792
 793 Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 794   //
 795   //
 796   //
 797   Float_t value=0;
 798   value += fParamSQPar[dim][type][0];
 799   value += fParamSQPar[dim][type][1]*z;
 800   value += fParamSQPar[dim][type][2]*angle*angle;
 801   value += fParamSQPar[dim][type][3]*z*z;
 802   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
 803   value += fParamSQPar[dim][type][5]*z*angle*angle;
 804   value += fParamSQPar[dim][type][6]*z/Qmean;
 805   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
 806   Float_t valueMean = GetError0Par(dim,type,z,angle);
 807   value -= 0.35*0.35*valueMean*valueMean;
 808   value  = TMath::Sqrt(TMath::Abs(value));
 809   return value;
 810
 811
 812 }
 813
 814 Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
 815   //
 816   // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
 817   //
 818   Float_t value=0;
 819   value += fParamRMS0[dim][type][0];
 820   value += fParamRMS0[dim][type][1]*z;
 821   value += fParamRMS0[dim][type][2]*angle*angle;
 822   value  = TMath::Sqrt(TMath::Abs(value));
 823   return value;
 824 }
 825
 826 Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
 827   //
 828   // calculate mean RMS of cluster - z,angle - pad length scalling
 829   //
 830   Float_t value=0;
 831   Float_t length=0.75;
 832   if (type==1) length=1;
 833   if (type==2) length=1.5;
 834   if (type==0){
 835     value += fParamRMS1[dim][0];
 836   }else{
 837     value += fParamRMS1[dim][1];
 838   }
 839   value += fParamRMS1[dim][2]*z;
 840   value += fParamRMS1[dim][3]*angle*angle*length*length;
 841   value  = TMath::Sqrt(TMath::Abs(value));
 842   return value;
 843 }
 844
 845 Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 846   //
 847   // calculate mean RMS of cluster - z,angle, Q dependence
 848   //
 849   Float_t value=0;
 850   value += fParamRMSQ[dim][type][0];
 851   value += fParamRMSQ[dim][type][1]*z;
 852   value += fParamRMSQ[dim][type][2]*angle*angle;
 853   value += fParamRMSQ[dim][type][3]*z/Qmean;
 854   value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
 855   value  = TMath::Sqrt(TMath::Abs(value));
 856   return value;
 857 }
 858
 859 Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
 860   //
 861   // calculates RMS of signal shape fluctuation
 862   //
 863   Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
 864   Float_t value  = fRMSSigmaFit[dim][type][0];
 865   value+=  fRMSSigmaFit[dim][type][1]*mean;
 866   return value;
 867 }
 868
 869 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){
 870   //
 871   // calculates vallue - sigma distortion contribution
 872   //
 873   Double_t value =0;
 874   //
 875   Float_t rmsMeanQ  = GetRMSQ(dim,type,z,angle,Qmean);
 876   if (rmsL<rmsMeanQ) return value;
 877   //
 878   Float_t rmsSigma  = GetRMSSigma(dim,type,z,angle,Qmean);
 879   //
 880   if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
 881     //1.5 sigma cut on mean
 882     value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
 883   }else{
 884     if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
 885       //3 sigma cut on local
 886       value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
 887     }
 888   }
 889   return TMath::Sqrt(TMath::Abs(value));
 890 }
 891
 892
 893
 894 void AliTPCClusterParam::FitData(TTree * tree){
 895   //
 896   // make fits for error param and shape param
 897   //
 898   FitResol(tree);
 899   FitRMS(tree);
 900
 901 }
 902
 903 void AliTPCClusterParam::FitResol(TTree * tree){
 904   //
 905   SetInstance(this);
 906   for (Int_t idir=0;idir<2; idir++){
 907     for (Int_t itype=0; itype<3; itype++){
 908       Float_t param0[10];
 909       Float_t error0[10];
 910       // model error param
 911       FitResol0(tree, idir, itype,param0,error0);
 912       printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
 913       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
 914       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
 915       for (Int_t ipar=0;ipar<4; ipar++){
 916         fParamS0[idir][itype][ipar] = param0[ipar];
 917         fErrorS0[idir][itype][ipar] = param0[ipar];
 918       }
 919       // error param with parabolic correction
 920       FitResol0Par(tree, idir, itype,param0,error0);
 921       printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
 922       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]);
 923       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]);
 924       for (Int_t ipar=0;ipar<7; ipar++){
 925         fParamS0Par[idir][itype][ipar] = param0[ipar];
 926         fErrorS0Par[idir][itype][ipar] = param0[ipar];
 927       }
 928       //
 929       FitResolQ(tree, idir, itype,param0,error0);
 930       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
 931       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
 932       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
 933       for (Int_t ipar=0;ipar<6; ipar++){
 934         fParamSQ[idir][itype][ipar] = param0[ipar];
 935         fErrorSQ[idir][itype][ipar] = param0[ipar];
 936       }
 937       //
 938       FitResolQPar(tree, idir, itype,param0,error0);
 939       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
 940       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]);
 941       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]);
 942       for (Int_t ipar=0;ipar<9; ipar++){
 943         fParamSQPar[idir][itype][ipar] = param0[ipar];
 944         fErrorSQPar[idir][itype][ipar] = param0[ipar];
 945       }
 946     }
 947   }
 948   //
 949   printf("Resol z-scaled\n");
 950   for (Int_t idir=0;idir<2; idir++){
 951     Float_t param0[4];
 952     Float_t error0[4];
 953     FitResol1(tree, idir,param0,error0);
 954     printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
 955     printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
 956     printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
 957     for (Int_t ipar=0;ipar<4; ipar++){
 958       fParamS1[idir][ipar] = param0[ipar];
 959       fErrorS1[idir][ipar] = param0[ipar];
 960     }
 961   }
 962
 963   for (Int_t idir=0;idir<2; idir++){
 964     printf("\nDirection %d\n",idir);
 965     printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
 966     for (Int_t itype=0; itype<3; itype++){
 967       Float_t length=0.75;
 968       if (itype==1) length=1;
 969       if (itype==2) length=1.5;
 970       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));
 971     }
 972   }
 973 }
 974
 975
 976
 977 void AliTPCClusterParam::FitRMS(TTree * tree){
 978   //
 979   SetInstance(this);
 980   for (Int_t idir=0;idir<2; idir++){
 981     for (Int_t itype=0; itype<3; itype++){
 982       Float_t param0[6];
 983       Float_t error0[6];
 984       FitRMS0(tree, idir, itype,param0,error0);
 985       printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
 986       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
 987       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
 988       for (Int_t ipar=0;ipar<4; ipar++){
 989         fParamRMS0[idir][itype][ipar] = param0[ipar];
 990         fErrorRMS0[idir][itype][ipar] = param0[ipar];
 991       }
 992       FitRMSQ(tree, idir, itype,param0,error0);
 993       printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
 994       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
 995       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
 996       for (Int_t ipar=0;ipar<6; ipar++){
 997         fParamRMSQ[idir][itype][ipar] = param0[ipar];
 998         fErrorRMSQ[idir][itype][ipar] = param0[ipar];
 999       }
1000     }
1001   }
1002   //
1003   printf("RMS z-scaled\n");
1004   for (Int_t idir=0;idir<2; idir++){
1005     Float_t param0[5];
1006     Float_t error0[5];
1007     FitRMS1(tree, idir,param0,error0);
1008     printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1009     printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1010     printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1011     for (Int_t ipar=0;ipar<5; ipar++){
1012       fParamRMS1[idir][ipar] = param0[ipar];
1013       fErrorRMS1[idir][ipar] = param0[ipar];
1014     }
1015   }
1016
1017   for (Int_t idir=0;idir<2; idir++){
1018     printf("\nDirection %d\n",idir);
1019     printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1020     for (Int_t itype=0; itype<3; itype++){
1021       Float_t length=0.75;
1022       if (itype==1) length=1;
1023       if (itype==2) length=1.5;
1024       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);
1025       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);
1026     }
1027   }
1028   //
1029   // Fit RMS sigma
1030   //
1031   printf("RMS fluctuation  parameterization \n");
1032   for (Int_t idir=0;idir<2; idir++){
1033     for (Int_t itype=0; itype<3; itype++){
1034       Float_t param0[5];
1035       Float_t error0[5];
1036       FitRMSSigma(tree, idir,itype,param0,error0);
1037       printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1038       for (Int_t ipar=0;ipar<2; ipar++){
1039         fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1040       }
1041     }
1042   }
1043   //
1044   // store systematic error end RMS fluctuation parameterization
1045   //
1046   TH1F hratio("hratio","hratio",100,-0.1,0.1);
1047   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1048   fErrorRMSSys[0] = hratio.GetRMS();
1049   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1050   fErrorRMSSys[1] = hratio.GetRMS();
1051   TH1F hratioR("hratioR","hratioR",100,0,0.2);
1052   tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1053   fRMSSigmaRatio[0][0]=hratioR.GetMean();
1054   fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1055   tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1056   fRMSSigmaRatio[1][0]=hratioR.GetMean();
1057   fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1058 }
1059
1060 void AliTPCClusterParam::Test(TTree * tree, const char *output){
1061   //
1062   // Draw standard quality histograms to output file
1063   //
1064   SetInstance(this);
1065   TFile f(output,"recreate");
1066   f.cd();
1067   //
1068   // 1D histograms - resolution
1069   //
1070   for (Int_t idim=0; idim<2; idim++){
1071     for (Int_t ipad=0; ipad<3; ipad++){
1072       char hname1[300];
1073       char hcut1[300];
1074       char hexp1[300];
1075       //
1076       sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1077       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1078       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1079       TH1F  his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1080       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1081       tree->Draw(hexp1,hcut1,"");
1082       his1DRel0.Write();
1083       //
1084       sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1085       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1086       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1087       TH1F  his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1088       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1089       tree->Draw(hexp1,hcut1,"");
1090       his1DRel0Par.Write();
1091       //
1092     }
1093   }
1094   //
1095   // 2D histograms - resolution
1096   //
1097   for (Int_t idim=0; idim<2; idim++){
1098     for (Int_t ipad=0; ipad<3; ipad++){
1099       char hname1[300];
1100       char hcut1[300];
1101       char hexp1[300];
1102       //
1103       sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1104       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1105       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1106       TProfile2D  profDRel0(hname1, hname1, 6,0,250,6,0,1);
1107       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1108       tree->Draw(hexp1,hcut1,"");
1109       profDRel0.Write();
1110       //
1111       sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1112       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1113       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1114       TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1115       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1116       tree->Draw(hexp1,hcut1,"");
1117       profDRel0Par.Write();
1118       //
1119     }
1120   }
1121 }
1122
1123
1124
1125 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1126   //
1127   // Print param Information
1128   //
1129
1130   //
1131   // Error parameterization
1132   //
1133   printf("\nResolution Scaled factors\n");
1134   printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1135   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])),
1136          TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
1137   for (Int_t ipad=0; ipad<3; ipad++){
1138     Float_t length=0.75;
1139     if (ipad==1) length=1;
1140     if (ipad==2) length=1.5;
1141     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1142            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1143            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1144            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1145            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
1146   }
1147   for (Int_t ipad=0; ipad<3; ipad++){
1148     Float_t length=0.75;
1149     if (ipad==1) length=1;
1150     if (ipad==2) length=1.5;
1151     printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1152            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1153            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1154            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1155            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
1156   }
1157   printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1158          TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1159
1160   for (Int_t ipad=0; ipad<3; ipad++){
1161     Float_t length=0.75;
1162     if (ipad==1) length=1;
1163     if (ipad==2) length=1.5;
1164     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1165            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1166            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1167            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1168            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
1169   }
1170   for (Int_t ipad=0; ipad<3; ipad++){
1171     Float_t length=0.75;
1172     if (ipad==1) length=1;
1173     if (ipad==2) length=1.5;
1174     printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1175            TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1176            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1177            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1178            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
1179   }
1180
1181   //
1182   // RMS scaling
1183   //
1184   printf("\n");
1185   printf("\nRMS Scaled factors\n");
1186   printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1187   printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1188          TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1189          TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1190          TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1191          TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1192          TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
1193   for (Int_t ipad=0; ipad<3; ipad++){
1194     Float_t length=0.75;
1195     if (ipad==1) length=1;
1196     if (ipad==2) length=1.5;
1197     if (ipad==0){
1198       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1199              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1200              0.,
1201              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1202              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1203              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1204     }else{
1205       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1206              0.,
1207              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1208              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1209              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1210              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1211     }
1212   }
1213   printf("\n");
1214   printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1215          TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1216          TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1217          TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1218          TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1219          TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
1220   for (Int_t ipad=0; ipad<3; ipad++){
1221     Float_t length=0.75;
1222     if (ipad==1) length=1;
1223     if (ipad==2) length=1.5;
1224     if (ipad==0){
1225       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1226              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1227              0.,
1228              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1229              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1230              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1231     }else{
1232       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1233              0.,
1234              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1235              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1236              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1237              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1238     }
1239   }
1240 }
1241
1242
1243
1244
1245
1246 Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1247   // get Q normalization
1248   // type - 0 Qtot 1 Qmax
1249   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1250   //
1251   //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);
1252
1253   if (fQNorm==0) return 0;
1254   TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1255   if (!norm) return 0;
1256   TVectorD &no  = *norm;
1257   Float_t   res = no[0]+
1258     no[1]*dr+
1259     no[2]*ty+
1260     no[3]*tz+
1261     no[4]*dr*ty+
1262     no[5]*dr*tz+
1263     no[6]*ty*tz+
1264     no[7]*dr*dr+
1265     no[8]*ty*ty+
1266     no[9]*tz*tz;
1267   res/=no[0];
1268   return res;
1269 }
1270
1271
1272
1273 void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
1274   //
1275   // set normalization
1276   //
1277   // type - 0 Qtot 1 Qmax
1278   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1279   //
1280
1281   if (fQNorm==0) fQNorm = new TObjArray(6);
1282   fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1283 }
1284
1285
1286
1287
1288 /*
1289
1290  gSystem->Load("libSTAT.so")
1291  TStatToolkit toolkit;
1292  Double_t chi2;
1293  TVectorD fitParam;
1294  TMatrixD covMatrix;
1295  Int_t npoints;
1296
1297  TString fstring="";
1298  //
1299  fstring+="(Zm/250.)++"
1300  fstring+="(Zm/250.)^2++"
1301  fstring+="(Zm/250.)^3++"
1302  fstring+="AngleM++"
1303  fstring+="(AngleM)^2++"
1304  fstring+="(AngleM)^3++"
1305  fstring+="(Zm/250.)*AngleM++"
1306  fstring+="(Zm/250.)^2*AngleM++"
1307  fstring+="(Zm/250.)*AngleM^2++"
1308
1309  TString *strRMSY0 = toolkit.FitPlane(treeRes,"RMSm^2",fstring->Data(), "Dim==0&&Pad==0&&QMean<0", chi2,npoints,fitParam,covMatrix);
1310
1311  treeRes->SetAlias("RMSY0",strRMSY0->Data());
1312
1313 */