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