Script supersceeded by AliForwarddNdetaTask.C and
[u/mrichter/AliRoot.git] / 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/OCDB");
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 <TH3F.h>
95 #include <TProfile2D.h>
96 #include <TVectorD.h>
97 #include <TObjArray.h>
98 #include "AliTPCcalibDB.h"
99 #include "AliTPCParam.h"
100
101 #include "AliMathBase.h"
102
103 ClassImp(AliTPCClusterParam)
104
105
106 AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
107
108
109 /*
110   Example usage fitting parameterization:
111   TFile fres("resol.root");    //tree with resolution and shape 
112   TTree * treeRes =(TTree*)fres.Get("Resol");
113   
114   AliTPCClusterParam param;
115   param.SetInstance(&param);
116   param.FitResol(treeRes);
117   param.FitRMS(treeRes);
118   TFile fparam("TPCClusterParam.root","recreate");
119   param.Write("Param");
120   //
121   //
122   TFile fparam("TPCClusterParam.root");
123   AliTPCClusterParam *param2  =  (AliTPCClusterParam *) fparam.Get("Param"); 
124   param2->SetInstance(param2);
125   param2->Test(treeRes);
126   
127
128   treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
129
130 */
131
132
133
134
135 //_ singleton implementation __________________________________________________
136 AliTPCClusterParam* AliTPCClusterParam::Instance()
137 {
138   //
139   // Singleton implementation
140   // Returns an instance of this class, it is created if neccessary
141   //
142   if (fgInstance == 0){
143     fgInstance = new AliTPCClusterParam();
144   }
145   return fgInstance;
146 }
147
148
149 AliTPCClusterParam::AliTPCClusterParam():
150   TObject(),
151   fRatio(0),
152   fQNorm(0),
153   fQNormCorr(0),
154   fQNormHis(0),
155   fQpadTnorm(0),           // q pad normalization - Total charge
156   fQpadMnorm(0)            // q pad normalization - Max charge
157   //
158 {
159   //
160   // Default constructor
161   //
162   fPosQTnorm[0] = 0;   fPosQTnorm[1] = 0;   fPosQTnorm[2] = 0; 
163   fPosQMnorm[0] = 0;   fPosQMnorm[1] = 0;   fPosQMnorm[2] = 0; 
164   //
165   fPosYcor[0]   = 0;   fPosYcor[1]   = 0;   fPosYcor[2]   = 0; 
166   fPosZcor[0]   = 0;   fPosZcor[1]   = 0;   fPosZcor[2]   = 0; 
167 }
168
169 AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
170   TObject(param),
171   fRatio(0),
172   fQNorm(0),
173   fQNormCorr(0),
174   fQNormHis(0),
175   fQpadTnorm(new TVectorD(*(param.fQpadTnorm))),           // q pad normalization - Total charge
176   fQpadMnorm(new TVectorD(*(param.fQpadMnorm)))            // q pad normalization - Max charge
177
178 {
179   //
180   // copy constructor
181   //
182   memcpy(this, &param,sizeof(AliTPCClusterParam));
183   if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
184   if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
185   //
186   if (param.fPosQTnorm[0]){
187     fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
188     fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
189     fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
190     //
191     fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
192     fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
193     fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
194   }
195   if (param.fPosYcor[0]){
196     fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
197     fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
198     fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
199     //
200     fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
201     fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
202     fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
203   }
204   
205 }
206
207
208 AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
209   //
210   // Assignment operator
211   //
212   if (this != &param) {
213     memcpy(this, &param,sizeof(AliTPCClusterParam));
214     if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
215     if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
216     if (param.fPosQTnorm[0]){
217       fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
218       fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
219       fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
220       //
221       fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
222       fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
223       fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
224     }
225     if (param.fPosYcor[0]){
226       fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
227       fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
228       fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
229       //
230       fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
231       fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
232       fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
233     }
234   }
235   return *this;
236 }
237
238
239 AliTPCClusterParam::~AliTPCClusterParam(){
240   //
241   // destructor
242   //
243   if (fQNorm) fQNorm->Delete();
244   if (fQNormCorr) delete fQNormCorr;
245   if (fQNormHis) fQNormHis->Delete();
246   delete fQNorm;
247   delete fQNormHis;
248   if (fPosQTnorm[0]){
249     delete fPosQTnorm[0];
250     delete fPosQTnorm[1];
251     delete fPosQTnorm[2];
252     //
253     delete fPosQMnorm[0];
254     delete fPosQMnorm[1];
255     delete fPosQMnorm[2];
256   }
257   if (fPosYcor[0]){
258     delete fPosYcor[0];
259     delete fPosYcor[1];
260     delete fPosYcor[2];
261     //
262     delete fPosZcor[0];
263     delete fPosZcor[1];
264     delete fPosZcor[2];
265   }
266 }
267
268
269 void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
270   //
271   // Fit z - angular dependence of resolution 
272   //
273   // Int_t dim=0, type=0;
274   TString varVal;
275   varVal="Resol:AngleM:Zm";
276   TString varErr;
277   varErr="Sigma:AngleS:Zs";
278   TString varCut;
279   varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
280   //
281   Int_t entries = tree->Draw(varVal.Data(),varCut);
282   Float_t px[10000], py[10000], pz[10000];
283   Float_t ex[10000], ey[10000], ez[10000];
284   //
285   tree->Draw(varErr.Data(),varCut);  
286   for (Int_t ipoint=0; ipoint<entries; ipoint++){
287     ex[ipoint]= tree->GetV3()[ipoint];
288     ey[ipoint]= tree->GetV2()[ipoint];
289     ez[ipoint]= tree->GetV1()[ipoint];
290   } 
291   tree->Draw(varVal.Data(),varCut);
292   for (Int_t ipoint=0; ipoint<entries; ipoint++){
293     px[ipoint]= tree->GetV3()[ipoint];
294     py[ipoint]= tree->GetV2()[ipoint];
295     pz[ipoint]= tree->GetV1()[ipoint];
296   }
297   
298   //  
299   TLinearFitter fitter(3,"hyp2");
300   for (Int_t ipoint=0; ipoint<entries; ipoint++){
301     Float_t val = pz[ipoint]*pz[ipoint];
302     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
303     Double_t x[2];
304     x[0] = px[ipoint];
305     x[1] = py[ipoint]*py[ipoint];
306     fitter.AddPoint(x,val,err);
307   }
308   fitter.Eval();
309   TVectorD param(3);
310   fitter.GetParameters(param);
311   param0[0] = param[0];
312   param0[1] = param[1];
313   param0[2] = param[2];
314   Float_t chi2 =  fitter.GetChisquare()/entries;
315   param0[3] = chi2;
316   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
317   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
318   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
319 }
320
321
322 void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
323   //
324   // Fit z - angular dependence of resolution 
325   //
326   // Int_t dim=0, type=0;
327  TString varVal;
328   varVal="Resol:AngleM:Zm";
329  TString varErr;
330   varErr="Sigma:AngleS:Zs";
331  TString varCut;
332   varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
333   //
334   Int_t entries = tree->Draw(varVal.Data(),varCut);
335   Float_t px[10000], py[10000], pz[10000];
336   Float_t ex[10000], ey[10000], ez[10000];
337   //
338   tree->Draw(varErr.Data(),varCut);  
339   for (Int_t ipoint=0; ipoint<entries; ipoint++){
340     ex[ipoint]= tree->GetV3()[ipoint];
341     ey[ipoint]= tree->GetV2()[ipoint];
342     ez[ipoint]= tree->GetV1()[ipoint];
343   } 
344   tree->Draw(varVal.Data(),varCut);
345   for (Int_t ipoint=0; ipoint<entries; ipoint++){
346     px[ipoint]= tree->GetV3()[ipoint];
347     py[ipoint]= tree->GetV2()[ipoint];
348     pz[ipoint]= tree->GetV1()[ipoint];
349   }
350   
351   //  
352   TLinearFitter fitter(6,"hyp5");
353   for (Int_t ipoint=0; ipoint<entries; ipoint++){
354     Float_t val = pz[ipoint]*pz[ipoint];
355     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
356     Double_t x[6];
357     x[0] = px[ipoint];
358     x[1] = py[ipoint]*py[ipoint];
359     x[2] = x[0]*x[0];
360     x[3] = x[1]*x[1];
361     x[4] = x[0]*x[1];
362     fitter.AddPoint(x,val,err);
363   }
364   fitter.Eval();
365   TVectorD param(6);
366   fitter.GetParameters(param);
367   param0[0] = param[0];
368   param0[1] = param[1];
369   param0[2] = param[2];
370   param0[3] = param[3];
371   param0[4] = param[4];
372   param0[5] = param[5];
373   Float_t chi2 =  fitter.GetChisquare()/entries;
374   param0[6] = chi2;
375   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
376   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
377   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
378   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
379   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
380   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
381 }
382
383
384
385
386
387 void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
388   //
389   // Fit z - angular dependence of resolution - pad length scaling 
390   //
391   // Int_t dim=0, type=0;
392  TString varVal;
393   varVal="Resol:AngleM*sqrt(Length):Zm/Length";
394  TString varErr;
395   varErr="Sigma:AngleS:Zs";
396  TString varCut;
397   varCut=Form("Dim==%d&&QMean<0",dim);
398   //
399   Int_t entries = tree->Draw(varVal.Data(),varCut);
400   Float_t px[10000], py[10000], pz[10000];
401   Float_t ex[10000], ey[10000], ez[10000];
402   //
403   tree->Draw(varErr.Data(),varCut);  
404   for (Int_t ipoint=0; ipoint<entries; ipoint++){
405     ex[ipoint]= tree->GetV3()[ipoint];
406     ey[ipoint]= tree->GetV2()[ipoint];
407     ez[ipoint]= tree->GetV1()[ipoint];
408   } 
409   tree->Draw(varVal.Data(),varCut);
410   for (Int_t ipoint=0; ipoint<entries; ipoint++){
411     px[ipoint]= tree->GetV3()[ipoint];
412     py[ipoint]= tree->GetV2()[ipoint];
413     pz[ipoint]= tree->GetV1()[ipoint];
414   }
415   
416   //  
417   TLinearFitter fitter(3,"hyp2");
418   for (Int_t ipoint=0; ipoint<entries; ipoint++){
419     Float_t val = pz[ipoint]*pz[ipoint];
420     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
421     Double_t x[2];
422     x[0] = px[ipoint];
423     x[1] = py[ipoint]*py[ipoint];
424     fitter.AddPoint(x,val,err);
425   }
426   fitter.Eval();
427   TVectorD param(3);
428   fitter.GetParameters(param);
429   param0[0] = param[0];
430   param0[1] = param[1];
431   param0[2] = param[2];
432   Float_t chi2 =  fitter.GetChisquare()/entries;
433   param0[3] = chi2;
434   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
435   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
436   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
437 }
438
439 void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
440   //
441   // Fit z - angular dependence of resolution - Q scaling 
442   //
443   // Int_t dim=0, type=0;
444  TString varVal;
445   varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
446   char varVal0[100];
447   sprintf(varVal0,"Resol:AngleM:Zm");
448   //
449  TString varErr;
450   varErr="Sigma:AngleS:Zs";
451  TString varCut;
452   varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
453   //
454   Int_t entries = tree->Draw(varVal.Data(),varCut);
455   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
456   Float_t ex[20000], ey[20000], ez[20000];
457   //
458   tree->Draw(varErr.Data(),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.Data(),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   tree->Draw(varVal0,varCut);
471   for (Int_t ipoint=0; ipoint<entries; ipoint++){
472     pu[ipoint]= tree->GetV3()[ipoint];
473     pt[ipoint]= tree->GetV2()[ipoint];
474   }
475   
476   //  
477   TLinearFitter fitter(5,"hyp4");
478   for (Int_t ipoint=0; ipoint<entries; ipoint++){
479     Float_t val = pz[ipoint]*pz[ipoint];
480     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
481     Double_t x[4];
482     x[0] = pu[ipoint];
483     x[1] = pt[ipoint]*pt[ipoint];
484     x[2] = px[ipoint];
485     x[3] = py[ipoint]*py[ipoint];
486     fitter.AddPoint(x,val,err);
487   }
488
489   fitter.Eval();
490   TVectorD param(5);
491   fitter.GetParameters(param);
492   param0[0] = param[0];
493   param0[1] = param[1];
494   param0[2] = param[2];
495   param0[3] = param[3];
496   param0[4] = param[4];
497   Float_t chi2 =  fitter.GetChisquare()/entries;
498   param0[5] = chi2;
499   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
500   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
501   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
502   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
503   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
504 }
505
506 void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
507   //
508   // Fit z - angular dependence of resolution - Q scaling  - parabolic correction
509   //
510   // Int_t dim=0, type=0;
511  TString varVal;
512   varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
513   char varVal0[100];
514   sprintf(varVal0,"Resol:AngleM:Zm");
515   //
516  TString varErr;
517   varErr="Sigma:AngleS:Zs";
518  TString varCut;
519   varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
520   //
521   Int_t entries = tree->Draw(varVal.Data(),varCut);
522   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
523   Float_t ex[20000], ey[20000], ez[20000];
524   //
525   tree->Draw(varErr.Data(),varCut);  
526   for (Int_t ipoint=0; ipoint<entries; ipoint++){
527     ex[ipoint]= tree->GetV3()[ipoint];
528     ey[ipoint]= tree->GetV2()[ipoint];
529     ez[ipoint]= tree->GetV1()[ipoint];
530   } 
531   tree->Draw(varVal.Data(),varCut);
532   for (Int_t ipoint=0; ipoint<entries; ipoint++){
533     px[ipoint]= tree->GetV3()[ipoint];
534     py[ipoint]= tree->GetV2()[ipoint];
535     pz[ipoint]= tree->GetV1()[ipoint];
536   }
537   tree->Draw(varVal0,varCut);
538   for (Int_t ipoint=0; ipoint<entries; ipoint++){
539     pu[ipoint]= tree->GetV3()[ipoint];
540     pt[ipoint]= tree->GetV2()[ipoint];
541   }
542   
543   //  
544   TLinearFitter fitter(8,"hyp7");
545   for (Int_t ipoint=0; ipoint<entries; ipoint++){
546     Float_t val = pz[ipoint]*pz[ipoint];
547     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
548     Double_t x[7];
549     x[0] = pu[ipoint];
550     x[1] = pt[ipoint]*pt[ipoint];
551     x[2] = x[0]*x[0];
552     x[3] = x[1]*x[1];
553     x[4] = x[0]*x[1];
554     x[5] = px[ipoint];
555     x[6] = py[ipoint]*py[ipoint];
556     //
557     fitter.AddPoint(x,val,err);
558   }
559
560   fitter.Eval();
561   TVectorD param(8);
562   fitter.GetParameters(param);
563   param0[0] = param[0];
564   param0[1] = param[1];
565   param0[2] = param[2];
566   param0[3] = param[3];
567   param0[4] = param[4];
568   param0[5] = param[5];
569   param0[6] = param[6];
570   param0[7] = param[7];
571
572   Float_t chi2 =  fitter.GetChisquare()/entries;
573   param0[8] = chi2;
574   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
575   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
576   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
577   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
578   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
579   error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
580   error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
581   error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
582 }
583
584
585
586 void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
587   //
588   // Fit z - angular dependence of resolution 
589   //
590   // Int_t dim=0, type=0;
591  TString varVal;
592   varVal="RMSm:AngleM:Zm";
593  TString varErr;
594   varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
595  TString varCut;
596   varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
597   //
598   Int_t entries = tree->Draw(varVal.Data(),varCut);
599   Float_t px[10000], py[10000], pz[10000];
600   Float_t ex[10000], ey[10000], ez[10000];
601   //
602   tree->Draw(varErr.Data(),varCut);  
603   for (Int_t ipoint=0; ipoint<entries; ipoint++){
604     ex[ipoint]= tree->GetV3()[ipoint];
605     ey[ipoint]= tree->GetV2()[ipoint];
606     ez[ipoint]= tree->GetV1()[ipoint];
607   } 
608   tree->Draw(varVal.Data(),varCut);
609   for (Int_t ipoint=0; ipoint<entries; ipoint++){
610     px[ipoint]= tree->GetV3()[ipoint];
611     py[ipoint]= tree->GetV2()[ipoint];
612     pz[ipoint]= tree->GetV1()[ipoint];
613   }
614   
615   //  
616   TLinearFitter fitter(3,"hyp2");
617   for (Int_t ipoint=0; ipoint<entries; ipoint++){
618     Float_t val = pz[ipoint]*pz[ipoint];
619     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
620     Double_t x[2];
621     x[0] = px[ipoint];
622     x[1] = py[ipoint]*py[ipoint];
623     fitter.AddPoint(x,val,err);
624   }
625   fitter.Eval();
626   TVectorD param(3);
627   fitter.GetParameters(param);
628   param0[0] = param[0];
629   param0[1] = param[1];
630   param0[2] = param[2];
631   Float_t chi2 =  fitter.GetChisquare()/entries;
632   param0[3] = chi2;
633   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
634   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
635   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
636 }
637
638 void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
639   //
640   // Fit z - angular dependence of resolution - pad length scaling 
641   //
642   // Int_t dim=0, type=0;
643  TString varVal;
644   varVal="RMSm:AngleM*Length:Zm";
645  TString varErr;
646   varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad";
647  TString varCut;
648   varCut=Form("Dim==%d&&QMean<0",dim);
649   //
650   Int_t entries = tree->Draw(varVal.Data(),varCut);
651   Float_t px[10000], py[10000], pz[10000];
652   Float_t type[10000], ey[10000], ez[10000];
653   //
654   tree->Draw(varErr.Data(),varCut);  
655   for (Int_t ipoint=0; ipoint<entries; ipoint++){
656     type[ipoint] = tree->GetV3()[ipoint];
657     ey[ipoint]   = tree->GetV2()[ipoint];
658     ez[ipoint]   = tree->GetV1()[ipoint];
659   } 
660   tree->Draw(varVal.Data(),varCut);
661   for (Int_t ipoint=0; ipoint<entries; ipoint++){
662     px[ipoint]= tree->GetV3()[ipoint];
663     py[ipoint]= tree->GetV2()[ipoint];
664     pz[ipoint]= tree->GetV1()[ipoint];
665   }
666   
667   //  
668   TLinearFitter fitter(4,"hyp3");
669   for (Int_t ipoint=0; ipoint<entries; ipoint++){
670     Float_t val = pz[ipoint]*pz[ipoint];
671     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
672     Double_t x[3];
673     x[0] = (type[ipoint]<0.5)? 0.:1.;
674     x[1] = px[ipoint];
675     x[2] = py[ipoint]*py[ipoint];
676     fitter.AddPoint(x,val,err);
677   }
678   fitter.Eval();
679   TVectorD param(4);
680   fitter.GetParameters(param);
681   param0[0] = param[0];
682   param0[1] = param[0]+param[1];
683   param0[2] = param[2];
684   param0[3] = param[3];
685   Float_t chi2 =  fitter.GetChisquare()/entries;
686   param0[4] = chi2;
687   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
688   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
689   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
690   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
691 }
692
693 void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
694   //
695   // Fit z - angular dependence of resolution - Q scaling 
696   //
697   // Int_t dim=0, type=0;
698  TString varVal;
699   varVal="RMSm:AngleM/sqrt(QMean):Zm/QMean";
700   char varVal0[100];
701   sprintf(varVal0,"RMSm:AngleM:Zm");
702   //
703  TString varErr;
704   varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
705  TString varCut;
706   varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
707   //
708   Int_t entries = tree->Draw(varVal.Data(),varCut);
709   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
710   Float_t ex[20000], ey[20000], ez[20000];
711   //
712   tree->Draw(varErr.Data(),varCut);  
713   for (Int_t ipoint=0; ipoint<entries; ipoint++){
714     ex[ipoint]= tree->GetV3()[ipoint];
715     ey[ipoint]= tree->GetV2()[ipoint];
716     ez[ipoint]= tree->GetV1()[ipoint];
717   } 
718   tree->Draw(varVal.Data(),varCut);
719   for (Int_t ipoint=0; ipoint<entries; ipoint++){
720     px[ipoint]= tree->GetV3()[ipoint];
721     py[ipoint]= tree->GetV2()[ipoint];
722     pz[ipoint]= tree->GetV1()[ipoint];
723   }
724   tree->Draw(varVal0,varCut);
725   for (Int_t ipoint=0; ipoint<entries; ipoint++){
726     pu[ipoint]= tree->GetV3()[ipoint];
727     pt[ipoint]= tree->GetV2()[ipoint];
728   }
729   
730   //  
731   TLinearFitter fitter(5,"hyp4");
732   for (Int_t ipoint=0; ipoint<entries; ipoint++){
733     Float_t val = pz[ipoint]*pz[ipoint];
734     Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
735     Double_t x[4];
736     x[0] = pu[ipoint];
737     x[1] = pt[ipoint]*pt[ipoint];
738     x[2] = px[ipoint];
739     x[3] = py[ipoint]*py[ipoint];
740     fitter.AddPoint(x,val,err);
741   }
742
743   fitter.Eval();
744   TVectorD param(5);
745   fitter.GetParameters(param);
746   param0[0] = param[0];
747   param0[1] = param[1];
748   param0[2] = param[2];
749   param0[3] = param[3];
750   param0[4] = param[4];
751   Float_t chi2 =  fitter.GetChisquare()/entries;
752   param0[5] = chi2;
753   error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
754   error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
755   error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
756   error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
757   error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
758 }
759
760
761 void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
762   //
763   // Fit z - angular dependence of resolution - Q scaling 
764   //
765   // Int_t dim=0, type=0;
766   TString varVal;
767   varVal="RMSs:RMSm";
768   //
769  TString varCut;
770   varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
771   //
772   Int_t entries = tree->Draw(varVal.Data(),varCut);
773   Float_t px[20000], py[20000];
774   //
775   tree->Draw(varVal.Data(),varCut);
776   for (Int_t ipoint=0; ipoint<entries; ipoint++){
777     px[ipoint]= tree->GetV2()[ipoint];
778     py[ipoint]= tree->GetV1()[ipoint];
779   }
780   TLinearFitter fitter(2,"pol1");
781   for (Int_t ipoint=0; ipoint<entries; ipoint++){
782     Float_t val = py[ipoint];
783     Float_t err = fRatio*px[ipoint];
784     Double_t x[4];
785     x[0] = px[ipoint];
786     if (err>0) fitter.AddPoint(x,val,err);
787   }
788   fitter.Eval();
789   param0[0]= fitter.GetParameter(0);
790   param0[1]= fitter.GetParameter(1);
791 }
792
793
794
795 Float_t  AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
796   //
797   //
798   //
799   Float_t value=0;
800   value += fParamS0[dim][type][0];
801   value += fParamS0[dim][type][1]*z;
802   value += fParamS0[dim][type][2]*angle*angle;
803   value  = TMath::Sqrt(TMath::Abs(value)); 
804   return value;
805 }
806
807
808 Float_t  AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
809   //
810   //
811   //
812   Float_t value=0;
813   value += fParamS0Par[dim][type][0];
814   value += fParamS0Par[dim][type][1]*z;
815   value += fParamS0Par[dim][type][2]*angle*angle;
816   value += fParamS0Par[dim][type][3]*z*z;
817   value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
818   value += fParamS0Par[dim][type][5]*z*angle*angle;
819   value  = TMath::Sqrt(TMath::Abs(value)); 
820   return value;
821 }
822
823
824
825 Float_t  AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
826   //
827   //
828   //
829   Float_t value=0;
830   Float_t length=0.75;
831   if (type==1) length=1;
832   if (type==2) length=1.5;
833   value += fParamS1[dim][0];
834   value += fParamS1[dim][1]*z/length;
835   value += fParamS1[dim][2]*angle*angle*length;
836   value  = TMath::Sqrt(TMath::Abs(value)); 
837   return value;
838 }
839
840 Float_t  AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
841   //
842   //
843   //
844   Float_t value=0;
845   value += fParamSQ[dim][type][0];
846   value += fParamSQ[dim][type][1]*z;
847   value += fParamSQ[dim][type][2]*angle*angle;
848   value += fParamSQ[dim][type][3]*z/Qmean;
849   value += fParamSQ[dim][type][4]*angle*angle/Qmean;
850   value  = TMath::Sqrt(TMath::Abs(value)); 
851   return value;
852
853
854 }
855
856 Float_t  AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
857   //
858   //
859   //
860   Float_t value=0;
861   value += fParamSQPar[dim][type][0];
862   value += fParamSQPar[dim][type][1]*z;
863   value += fParamSQPar[dim][type][2]*angle*angle;
864   value += fParamSQPar[dim][type][3]*z*z;
865   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
866   value += fParamSQPar[dim][type][5]*z*angle*angle;
867   value += fParamSQPar[dim][type][6]*z/Qmean;
868   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
869   value  = TMath::Sqrt(TMath::Abs(value)); 
870   return value;
871
872
873 }
874
875 Float_t  AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
876   //
877   //
878   //
879   Float_t value=0;
880   value += fParamSQPar[dim][type][0];
881   value += fParamSQPar[dim][type][1]*z;
882   value += fParamSQPar[dim][type][2]*angle*angle;
883   value += fParamSQPar[dim][type][3]*z*z;
884   value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
885   value += fParamSQPar[dim][type][5]*z*angle*angle;
886   value += fParamSQPar[dim][type][6]*z/Qmean;
887   value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
888   Float_t valueMean = GetError0Par(dim,type,z,angle);
889   value -= 0.35*0.35*valueMean*valueMean; 
890   value  = TMath::Sqrt(TMath::Abs(value)); 
891   return value;
892
893
894 }
895
896 Float_t  AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
897   //
898   // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
899   //
900   Float_t value=0;
901   value += fParamRMS0[dim][type][0];
902   value += fParamRMS0[dim][type][1]*z;
903   value += fParamRMS0[dim][type][2]*angle*angle;
904   value  = TMath::Sqrt(TMath::Abs(value)); 
905   return value;
906 }
907
908 Float_t  AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
909   //
910   // calculate mean RMS of cluster - z,angle - pad length scalling
911   //
912   Float_t value=0;
913   Float_t length=0.75;
914   if (type==1) length=1;
915   if (type==2) length=1.5;
916   if (type==0){
917     value += fParamRMS1[dim][0];
918   }else{
919     value += fParamRMS1[dim][1];
920   }
921   value += fParamRMS1[dim][2]*z;
922   value += fParamRMS1[dim][3]*angle*angle*length*length;
923   value  = TMath::Sqrt(TMath::Abs(value)); 
924   return value;
925 }
926
927 Float_t  AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
928   //
929   // calculate mean RMS of cluster - z,angle, Q dependence
930   //
931   Float_t value=0;
932   value += fParamRMSQ[dim][type][0];
933   value += fParamRMSQ[dim][type][1]*z;
934   value += fParamRMSQ[dim][type][2]*angle*angle;
935   value += fParamRMSQ[dim][type][3]*z/Qmean;
936   value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
937   value  = TMath::Sqrt(TMath::Abs(value)); 
938   return value;
939 }
940
941 Float_t  AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
942   //
943   // calculates RMS of signal shape fluctuation
944   //
945   Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
946   Float_t value  = fRMSSigmaFit[dim][type][0];
947   value+=  fRMSSigmaFit[dim][type][1]*mean;
948   return value;
949 }
950
951 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) const {
952   //
953   // calculates vallue - sigma distortion contribution
954   //
955   Double_t value =0;
956   //
957   Float_t rmsMeanQ  = GetRMSQ(dim,type,z,angle,Qmean);
958   if (rmsL<rmsMeanQ) return value;
959   //
960   Float_t rmsSigma  = GetRMSSigma(dim,type,z,angle,Qmean);
961   //
962   if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
963     //1.5 sigma cut on mean
964     value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
965   }else{
966     if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
967       //3 sigma cut on local
968       value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
969     }
970   }
971   return TMath::Sqrt(TMath::Abs(value));
972 }
973
974
975
976 void AliTPCClusterParam::FitData(TTree * tree){
977   //
978   // make fits for error param and shape param
979   //
980   FitResol(tree);
981   FitRMS(tree);
982
983 }
984
985 void AliTPCClusterParam::FitResol(TTree * tree){
986   //
987   SetInstance(this);
988   for (Int_t idir=0;idir<2; idir++){    
989     for (Int_t itype=0; itype<3; itype++){
990       Float_t param0[10];
991       Float_t error0[10];
992       // model error param
993       FitResol0(tree, idir, itype,param0,error0);
994       printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
995       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
996       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
997       for (Int_t ipar=0;ipar<4; ipar++){
998         fParamS0[idir][itype][ipar] = param0[ipar];     
999         fErrorS0[idir][itype][ipar] = param0[ipar];     
1000       } 
1001       // error param with parabolic correction
1002       FitResol0Par(tree, idir, itype,param0,error0);
1003       printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
1004       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]);
1005       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]);
1006       for (Int_t ipar=0;ipar<7; ipar++){
1007         fParamS0Par[idir][itype][ipar] = param0[ipar];  
1008         fErrorS0Par[idir][itype][ipar] = param0[ipar];  
1009       }
1010       //
1011       FitResolQ(tree, idir, itype,param0,error0);
1012       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1013       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1014       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1015       for (Int_t ipar=0;ipar<6; ipar++){
1016         fParamSQ[idir][itype][ipar] = param0[ipar];     
1017         fErrorSQ[idir][itype][ipar] = param0[ipar];     
1018       }
1019       //
1020       FitResolQPar(tree, idir, itype,param0,error0);
1021       printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
1022       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]);
1023       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]);
1024       for (Int_t ipar=0;ipar<9; ipar++){
1025         fParamSQPar[idir][itype][ipar] = param0[ipar];  
1026         fErrorSQPar[idir][itype][ipar] = param0[ipar];  
1027       }
1028     }
1029   }
1030   //
1031   printf("Resol z-scaled\n");
1032   for (Int_t idir=0;idir<2; idir++){    
1033     Float_t param0[4];
1034     Float_t error0[4];
1035     FitResol1(tree, idir,param0,error0);
1036     printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
1037     printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1038     printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1039     for (Int_t ipar=0;ipar<4; ipar++){
1040       fParamS1[idir][ipar] = param0[ipar];      
1041       fErrorS1[idir][ipar] = param0[ipar];      
1042     }
1043   }
1044
1045   for (Int_t idir=0;idir<2; idir++){
1046     printf("\nDirection %d\n",idir);
1047     printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1048     for (Int_t itype=0; itype<3; itype++){
1049       Float_t length=0.75;
1050       if (itype==1) length=1;
1051       if (itype==2) length=1.5;
1052       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));
1053     }
1054   }  
1055 }
1056
1057
1058
1059 void AliTPCClusterParam::FitRMS(TTree * tree){
1060   //
1061   SetInstance(this);
1062   for (Int_t idir=0;idir<2; idir++){    
1063     for (Int_t itype=0; itype<3; itype++){
1064       Float_t param0[6];
1065       Float_t error0[6];
1066       FitRMS0(tree, idir, itype,param0,error0);
1067       printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1068       printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1069       printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1070       for (Int_t ipar=0;ipar<4; ipar++){
1071         fParamRMS0[idir][itype][ipar] = param0[ipar];   
1072         fErrorRMS0[idir][itype][ipar] = param0[ipar];   
1073       }
1074       FitRMSQ(tree, idir, itype,param0,error0);
1075       printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1076       printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1077       printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1078       for (Int_t ipar=0;ipar<6; ipar++){
1079         fParamRMSQ[idir][itype][ipar] = param0[ipar];   
1080         fErrorRMSQ[idir][itype][ipar] = param0[ipar];   
1081       }
1082     }
1083   }
1084   //
1085   printf("RMS z-scaled\n");
1086   for (Int_t idir=0;idir<2; idir++){    
1087     Float_t param0[5];
1088     Float_t error0[5];
1089     FitRMS1(tree, idir,param0,error0);
1090     printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1091     printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1092     printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1093     for (Int_t ipar=0;ipar<5; ipar++){
1094       fParamRMS1[idir][ipar] = param0[ipar];    
1095       fErrorRMS1[idir][ipar] = param0[ipar];    
1096     }
1097   }
1098
1099   for (Int_t idir=0;idir<2; idir++){
1100     printf("\nDirection %d\n",idir);
1101     printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1102     for (Int_t itype=0; itype<3; itype++){
1103       Float_t length=0.75;
1104       if (itype==1) length=1;
1105       if (itype==2) length=1.5;
1106       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);
1107       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);
1108     }
1109   }  
1110   //
1111   // Fit RMS sigma
1112   //
1113   printf("RMS fluctuation  parameterization \n");
1114   for (Int_t idir=0;idir<2; idir++){    
1115     for (Int_t itype=0; itype<3; itype++){ 
1116       Float_t param0[5];
1117       Float_t error0[5];
1118       FitRMSSigma(tree, idir,itype,param0,error0); 
1119       printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1120       for (Int_t ipar=0;ipar<2; ipar++){
1121         fRMSSigmaFit[idir][itype][ipar] = param0[ipar]; 
1122       }
1123     }
1124   } 
1125   //
1126   // store systematic error end RMS fluctuation parameterization
1127   //
1128   TH1F hratio("hratio","hratio",100,-0.1,0.1);
1129   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1130   fErrorRMSSys[0] = hratio.GetRMS();
1131   tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1132   fErrorRMSSys[1] = hratio.GetRMS();
1133   TH1F hratioR("hratioR","hratioR",100,0,0.2);
1134   tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1135   fRMSSigmaRatio[0][0]=hratioR.GetMean();
1136   fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1137   tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1138   fRMSSigmaRatio[1][0]=hratioR.GetMean();
1139   fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1140 }
1141
1142 void AliTPCClusterParam::Test(TTree * tree, const char *output){
1143   //
1144   // Draw standard quality histograms to output file
1145   //
1146   SetInstance(this);
1147   TFile f(output,"recreate");
1148   f.cd();
1149   //
1150   // 1D histograms - resolution
1151   //
1152   for (Int_t idim=0; idim<2; idim++){
1153     for (Int_t ipad=0; ipad<3; ipad++){
1154       char hname1[300];
1155       char hcut1[300];
1156       char hexp1[300];
1157       //
1158       sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1159       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1160       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1161       TH1F  his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1162       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1163       tree->Draw(hexp1,hcut1,"");
1164       his1DRel0.Write();
1165       //
1166       sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1167       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1168       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1169       TH1F  his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1170       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1171       tree->Draw(hexp1,hcut1,"");
1172       his1DRel0Par.Write();
1173       //
1174     }
1175   }
1176   //
1177   // 2D histograms - resolution
1178   //
1179   for (Int_t idim=0; idim<2; idim++){
1180     for (Int_t ipad=0; ipad<3; ipad++){
1181       char hname1[300];
1182       char hcut1[300];
1183       char hexp1[300];
1184       //
1185       sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1186       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1187       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1188       TProfile2D  profDRel0(hname1, hname1, 6,0,250,6,0,1);
1189       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1190       tree->Draw(hexp1,hcut1,"");
1191       profDRel0.Write();
1192       //
1193       sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1194       sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1195       sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1196       TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1197       sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1198       tree->Draw(hexp1,hcut1,"");
1199       profDRel0Par.Write();
1200       //
1201     }
1202   }
1203 }
1204
1205
1206
1207 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1208   //
1209   // Print param Information
1210   //
1211
1212   //
1213   // Error parameterization
1214   //
1215   printf("\nResolution Scaled factors\n");
1216   printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1217   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])),
1218          TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
1219   for (Int_t ipad=0; ipad<3; ipad++){
1220     Float_t length=0.75;
1221     if (ipad==1) length=1;
1222     if (ipad==2) length=1.5;    
1223     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad, 
1224            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1225            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1226            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1227            TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][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%dPar\t%f\t%f\t%f\t%f\n", ipad, 
1234            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1235            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1236            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1237            TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
1238   }
1239   printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1240          TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1241   
1242   for (Int_t ipad=0; ipad<3; ipad++){
1243     Float_t length=0.75;
1244     if (ipad==1) length=1;
1245     if (ipad==2) length=1.5;    
1246     printf("\t%d\t%f\t%f\t%f\t%f\n", ipad, 
1247            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1248            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1249            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1250            TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
1251   }
1252   for (Int_t ipad=0; ipad<3; ipad++){
1253     Float_t length=0.75;
1254     if (ipad==1) length=1;
1255     if (ipad==2) length=1.5;        
1256     printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad, 
1257            TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1258            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1259            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1260            TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
1261   }
1262   
1263   //
1264   // RMS scaling
1265   //
1266   printf("\n");
1267   printf("\nRMS Scaled factors\n");
1268   printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1269   printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n", 
1270          TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1271          TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1272          TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1273          TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1274          TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
1275   for (Int_t ipad=0; ipad<3; ipad++){
1276     Float_t length=0.75;
1277     if (ipad==1) length=1;
1278     if (ipad==2) length=1.5;    
1279     if (ipad==0){
1280       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, 
1281              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1282              0.,
1283              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1284              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1285              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1286     }else{
1287       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, 
1288              0.,
1289              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1290              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1291              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1292              TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));  
1293     }
1294   }
1295   printf("\n");
1296   printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n", 
1297          TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1298          TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1299          TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1300          TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1301          TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
1302   for (Int_t ipad=0; ipad<3; ipad++){
1303     Float_t length=0.75;
1304     if (ipad==1) length=1;
1305     if (ipad==2) length=1.5;    
1306     if (ipad==0){
1307       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, 
1308              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1309              0.,
1310              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1311              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1312              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1313     }else{
1314       printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, 
1315              0.,
1316              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1317              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1318              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1319              TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));  
1320     }
1321   }
1322 }
1323
1324
1325
1326
1327
1328 Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1329   // get Q normalization
1330   // type - 0 Qtot 1 Qmax
1331   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1332   //
1333   //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++++dr*dr++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
1334
1335   if (fQNorm==0) return 0;
1336   TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1337   if (!norm) return 0;
1338   TVectorD &no  = *norm;
1339   Float_t   res = 
1340     no[0]+
1341     no[1]*dr+
1342     no[2]*ty+
1343     no[3]*tz+
1344     no[4]*dr*ty+
1345     no[5]*dr*tz+
1346     no[6]*ty*tz+
1347     no[7]*dr*dr+
1348     no[8]*ty*ty+
1349     no[9]*tz*tz;
1350   res/=no[0];
1351   return res;
1352 }
1353
1354
1355
1356 Float_t AliTPCClusterParam::QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t p2, Float_t p3){
1357   // get Q normalization
1358   // type - 0 Qtot 1 Qmax
1359   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1360   //
1361
1362   if (fQNormHis==0) return 0;
1363   TH3F * norm = (TH3F*)fQNormHis->At(4*itype+ipad);
1364   if (!norm) return 1;
1365   p2=TMath::Abs(p2);
1366   dr=TMath::Min(dr,Float_t(norm->GetXaxis()->GetXmax()-norm->GetXaxis()->GetBinWidth(0)));
1367   dr=TMath::Max(dr,Float_t(norm->GetXaxis()->GetXmin()+norm->GetXaxis()->GetBinWidth(0)));
1368   //
1369   p2=TMath::Min(p2,Float_t(norm->GetYaxis()->GetXmax()-norm->GetYaxis()->GetBinWidth(0)));
1370   p2=TMath::Max(p2,Float_t(norm->GetYaxis()->GetXmin()+norm->GetYaxis()->GetBinWidth(0)));
1371   //
1372   p3=TMath::Min(p3,Float_t(norm->GetZaxis()->GetXmax()-norm->GetZaxis()->GetBinWidth(0)));
1373   p3=TMath::Max(p3,Float_t(norm->GetZaxis()->GetXmin()+norm->GetZaxis()->GetBinWidth(0)));
1374   //
1375   Double_t res = norm->GetBinContent(norm->FindBin(dr,p2,p3));
1376   if (res==0) res = norm->GetBinContent(norm->FindBin(0.5,0.5,0.5));  // This is just hack - to be fixed entries without 
1377
1378   return res;
1379 }
1380
1381
1382
1383 void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm){
1384   //
1385   // set normalization
1386   //
1387   // type - 0 Qtot 1 Qmax
1388   // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1389   //
1390
1391   if (fQNorm==0) fQNorm = new TObjArray(6);
1392   fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1393 }
1394
1395 void AliTPCClusterParam::ResetQnormCorr(){
1396   //
1397   //
1398   //
1399   if (!fQNormCorr) fQNormCorr= new TMatrixD(12,6);
1400   for (Int_t irow=0;irow<12; irow++)
1401     for (Int_t icol=0;icol<6; icol++){
1402       (*fQNormCorr)(irow,icol)=1.;             // default - no correction
1403       if (irow>5) (*fQNormCorr)(irow,icol)=0.; // default - no correction
1404     } 
1405 }
1406
1407 void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val){
1408   //
1409   // ipad        - pad type
1410   // itype       - 0- qtot 1-qmax
1411   // corrType    - 0 - s0y corr     - eff. PRF corr
1412   //             - 1 - s0z corr     - eff. TRF corr
1413   //             - 2 - d0y          - eff. diffusion correction y
1414   //             - 3 - d0z          - eff. diffusion correction
1415   //             - 4 - eff length   - eff.length - wire pitch + x diffsion
1416   //             - 5 - pad type normalization
1417   if (!fQNormCorr) {
1418     ResetQnormCorr();
1419   }
1420   //
1421   // eff shap parameterization matrix
1422   //
1423   // rows
1424   // itype*3+ipad  - itype=0 qtot itype=1 qmax ipad=0
1425   // 
1426   if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val;  // multiplicative correction
1427   if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val;  // additive       correction  
1428 }
1429
1430 Double_t  AliTPCClusterParam::GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const{
1431   //
1432   // see AliTPCClusterParam::SetQnormCorr
1433   //
1434   if (!fQNormCorr) return 0;
1435   return  (*fQNormCorr)(itype*3+ipad, corrType);
1436 }
1437
1438
1439 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){
1440   //
1441   // Make Q normalization as function of following parameters
1442   // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion
1443   // 1 - dp   - relative pad position 
1444   // 2 - dt   - relative time position
1445   // 3 - di   - drift length (norm to 1);
1446   // 4 - dq0  - Tot/Max charge
1447   // 5 - dq1  - Max/Tot charge
1448   // 6 - sy   - sigma y - shape
1449   // 7 - sz   - sigma z - shape
1450   //  
1451   
1452   //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain - 
1453   // Following variable used - correspondance to the our variable conventions  
1454   //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1455   Double_t dp = ((pad-int(pad)-0.5)*2.);
1456   //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1457   Double_t dt = ((time-int(time)-0.5)*2.);
1458   //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1459   Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1460   //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1461   Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1462   //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1463   Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1464   //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1465   Double_t sy  = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1466   //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1467   Double_t sz  = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1468   //
1469   //
1470   //
1471   TVectorD * pvec = 0;
1472   if (isMax){
1473     pvec = fPosQMnorm[ipad];
1474   }else{
1475     pvec = fPosQTnorm[ipad];    
1476   }
1477   TVectorD &param = *pvec;
1478   //
1479   // Eval part  - in correspondance with fit part from debug streamer
1480   // 
1481   Double_t result=param[0];
1482   Int_t index =1;
1483   //
1484   result+=dp*param[index++];                               //1
1485   result+=dt*param[index++];                               //2
1486   result+=dp*dp*param[index++];                             //3
1487   result+=dt*dt*param[index++];                             //4
1488   result+=dt*dt*dt*param[index++];                             //5
1489   result+=dp*dt*param[index++];                            //6
1490   result+=dp*dt*dt*param[index++];                          //7
1491   result+=(dq0)*param[index++];                            //8
1492   result+=(dq1)*param[index++];                            //9
1493   //
1494   //
1495   result+=dp*dp*(di)*param[index++];                        //10
1496   result+=dt*dt*(di)*param[index++];                        //11
1497   result+=dp*dp*sy*param[index++];                          //12
1498   result+=dt*sz*param[index++];                          //13
1499   result+=dt*dt*sz*param[index++];                          //14
1500   result+=dt*dt*dt*sz*param[index++];                          //15
1501   //
1502   result+=dp*dp*1*sy*sz*param[index++];                     //16
1503   result+=dt*sy*sz*param[index++];                       //17
1504   result+=dt*dt*sy*sz*param[index++];                       //18
1505   result+=dt*dt*dt*sy*sz*param[index++];                       //19
1506   //
1507   result+=dp*dp*(dq0)*param[index++];                       //20
1508   result+=dt*1*(dq0)*param[index++];                       //21
1509   result+=dt*dt*(dq0)*param[index++];                       //22
1510   result+=dt*dt*dt*(dq0)*param[index++];                       //23
1511   //
1512   result+=dp*dp*(dq1)*param[index++];                       //24
1513   result+=dt*(dq1)*param[index++];                       //25
1514   result+=dt*dt*(dq1)*param[index++];                       //26
1515   result+=dt*dt*dt*(dq1)*param[index++];                       //27
1516
1517   if (result<0.75) result=0.75;
1518   if (result>1.25) result=1.25;
1519
1520   return result;
1521   
1522 }
1523
1524
1525
1526
1527
1528 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*/){
1529
1530   //
1531   // Make postion correction
1532   // type - 0 - y correction
1533   //        1 - z correction
1534   // ipad - 0, 1, 2 - short, medium long pads 
1535   // pad  - float pad number          
1536   // time - float time bin number
1537   //    z - z of the cluster
1538   
1539   //
1540   //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
1541   //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
1542   //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
1543   //chainres->SetAlias("st","(sin(dt)-dt)");
1544   //
1545   //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1546
1547   //
1548   // Derived variables
1549   //
1550   Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
1551   Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
1552   Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
1553   Double_t st = (TMath::Sin(dt)-dt);
1554   //
1555   Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.)));
1556   //
1557   //
1558   //
1559   TVectorD * pvec = 0;
1560   if (type==0){
1561     pvec = fPosYcor[ipad];
1562   }else{
1563     pvec = fPosZcor[ipad];    
1564   }
1565   TVectorD &param = *pvec;
1566   //
1567   Double_t result=0;
1568   Int_t index =1;
1569
1570   if (type==0){
1571     // y corr
1572     result+=(dp)*param[index++];             //1
1573     result+=(dp)*di*param[index++];          //2
1574     //
1575     result+=(sp)*param[index++];             //3
1576     result+=(sp)*di*param[index++];          //4
1577   }
1578   if (type==1){
1579     result+=(dt)*param[index++];             //1
1580     result+=(dt)*di*param[index++];          //2
1581     //
1582     result+=(st)*param[index++];             //3
1583     result+=(st)*di*param[index++];          //4
1584   }
1585   if (TMath::Abs(result)>0.05) return 0;
1586   return result;
1587 }
1588
1589
1590
1591 Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
1592   //
1593   // 2 D gaus convoluted with angular effect
1594   // See in mathematica: 
1595   //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
1596   // 
1597   //TF1 f1("f1","AliTPCClusterParam::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
1598   //TF2 f2("f2","AliTPCClusterParam::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
1599   //
1600   const Float_t kEpsilon = 0.0001;
1601   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
1602     // small angular effect
1603     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
1604     return val;
1605   }
1606   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
1607   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
1608   //
1609   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);                        
1610   Double_t erf0 = AliMathBase::ErfFast( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
1611   Double_t erf1 = AliMathBase::ErfFast( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
1612   Double_t norm = 1./TMath::Sqrt(sigma2);
1613   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
1614   Double_t val  = norm*exp0*(erf0+erf1);
1615   return val;
1616 }
1617
1618
1619 Double_t  AliTPCClusterParam::GaussConvolutionTail(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
1620   //
1621   // 2 D gaus convoluted with angular effect and exponential tail in z-direction
1622   // tail integrated numerically 
1623   // Integral normalized to one
1624   // Mean at 0
1625   // 
1626   // TF1 f1t("f1t","AliTPCClusterParam::GaussConvolutionTail(0,x,0,0,0.5,0.5,0.9)",-5,5)
1627   Double_t sum =1, mean=0;
1628   // the COG of exponent
1629   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1630     mean+=iexp*TMath::Exp(-iexp/tau);
1631     sum +=TMath::Exp(-iexp/tau);
1632   }
1633   mean/=sum;
1634   //
1635   sum = 1;
1636   Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
1637   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1638     val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
1639     sum+=TMath::Exp(-iexp/tau);
1640   }
1641   return val/sum;
1642 }
1643
1644 Double_t  AliTPCClusterParam::GaussConvolutionGamma4(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
1645   //
1646   // 2 D gaus convoluted with angular effect and exponential tail in z-direction
1647   // tail integrated numerically 
1648   // Integral normalized to one
1649   // Mean at 0
1650   // 
1651   // TF1 f1g4("f1g4","AliTPCClusterParam::GaussConvolutionGamma4(0,x,0,0,0.5,0.2,1.6)",-5,5)
1652   // TF2 f2g4("f2g4","AliTPCClusterParam::GaussConvolutionGamma4(y,x,0,0,0.5,0.2,1.6)",-5,5,-5,5)
1653   Double_t sum =0, mean=0;
1654   // the COG of G4
1655   for (Float_t iexp=0;iexp<5;iexp+=0.2){
1656     Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1657     mean+=iexp*g4;
1658     sum +=g4;
1659   }
1660   mean/=sum;
1661   //
1662   sum = 0;
1663   Double_t val = 0;
1664   for (Float_t iexp=0;iexp<5;iexp+=0.2){ 
1665     Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1666     val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
1667     sum+=g4;
1668   }
1669   return val/sum;
1670 }
1671
1672 Double_t  AliTPCClusterParam::QmaxCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0, Float_t effPad, Float_t effDiff){
1673   //
1674   //
1675   // cpad      - pad (y) coordinate
1676   // ctime     - time(z) coordinate
1677   // ky        - dy/dx
1678   // kz        - dz/dx
1679   // rmsy0     - RF width in pad units
1680   // rmsz0     - RF width in time bin  units
1681   // effLength - contibution of PRF and diffusion
1682   // effDiff   - overwrite diffusion
1683
1684   // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1685   //  
1686   // Gaus width sy and sz is determined by RF width and diffusion 
1687   // Integral of Q is equal 1
1688   // Q max is calculated at position cpad, ctime
1689   // Example function:         
1690   //  TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000) 
1691   //
1692   AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters(); 
1693   Double_t padLength= param->GetPadPitchLength(sector,row);
1694   Double_t padWidth = param->GetPadPitchWidth(sector);
1695   Double_t zwidth   = param->GetZWidth();
1696   Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1697
1698   // diffusion in pad, time bin  units
1699   Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1700   Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1701   diffT*=effDiff;  //
1702   diffL*=effDiff;  //
1703   //
1704   // transform angular effect to pad units
1705   //
1706   Double_t pky   = ky*effLength/padWidth;
1707   Double_t pkz   = kz*effLength/zwidth;
1708   // position in pad unit
1709   Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1710   Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1711   //
1712   //
1713   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1714   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL); 
1715   //return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
1716   Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1717   return GaussConvolution(py,pz, pky,pkz,sy,sz)*length;
1718 }
1719
1720 Double_t  AliTPCClusterParam::QtotCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0,  Float_t qtot, Float_t thr, Float_t effPad, Float_t effDiff){
1721   //
1722   //
1723   // cpad      - pad (y) coordinate
1724   // ctime     - time(z) coordinate
1725   // ky        - dy/dx
1726   // kz        - dz/dx
1727   // rmsy0     - RF width in pad units
1728   // rmsz0     - RF width in time bin  units
1729   // qtot      - the sum of signal in cluster - without thr correction
1730   // thr       - threshold
1731   // effLength - contibution of PRF and diffusion
1732   // effDiff   - overwrite diffusion
1733
1734   // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1735   //  
1736   // Gaus width sy and sz is determined by RF width and diffusion 
1737   // Integral of Q is equal 1
1738   // Q max is calculated at position cpad, ctime
1739   //          
1740   //  
1741   //
1742   AliTPCParam * param   = AliTPCcalibDB::Instance()->GetParameters(); 
1743   Double_t padLength= param->GetPadPitchLength(sector,row);
1744   Double_t padWidth = param->GetPadPitchWidth(sector);
1745   Double_t zwidth   = param->GetZWidth();
1746   Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1747   //
1748   // diffusion in pad units
1749   Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1750   Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1751   diffT*=effDiff;  //
1752   diffL*=effDiff;  //
1753   //
1754   // transform angular effect to pad units 
1755   Double_t pky   = ky*effLength/padWidth;
1756   Double_t pkz   = kz*effLength/zwidth;
1757   // position in pad unit
1758   //  
1759   Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1760   Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1761   //
1762   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1763   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL); 
1764   //
1765   //
1766   //
1767   Double_t sumAll=0,sumThr=0;
1768   //
1769   Double_t corr =1;
1770   Double_t qnorm=qtot;
1771   for (Float_t iy=-3;iy<=3;iy+=1.)
1772     for (Float_t iz=-4;iz<=4;iz+=1.){
1773       //      Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);      
1774       Double_t val = GaussConvolution(py-iy,pz-iz, pky,pkz, sy,sz);      
1775       Double_t qlocal =qnorm*val;
1776       if (TMath::Abs(iy)<1.5&&TMath::Abs(iz)<1.5){
1777         sumThr+=qlocal;   // Virtual charge used in cluster finder
1778       }
1779       else{
1780         if (qlocal>thr && TMath::Abs(iz)<2.5&&TMath::Abs(iy)<2.5) sumThr+=qlocal;
1781       }
1782       sumAll+=qlocal;
1783     }
1784   if (sumAll>0&&sumThr>0) {
1785     corr=(sumThr)/sumAll;
1786   }
1787   //
1788   Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1789   return corr*length;
1790 }
1791
1792
1793
1794
1795
1796
1797