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