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