TPC/AliTPCClusterParam.cxx

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