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