TPC/AliTPCClusterParam.cxx

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