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