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