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