AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTPCClusterParam.cxx
Go to the documentation of this file.
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 
77 
78 #include "AliTPCClusterParam.h"
79 #include "TMath.h"
80 #include "TFile.h"
81 #include "TTree.h"
82 #include <TVectorF.h>
83 #include <TLinearFitter.h>
84 #include <TH1F.h>
85 #include <TH3F.h>
86 #include <TProfile2D.h>
87 #include <TVectorD.h>
88 #include <TObjArray.h>
89 #include "AliTPCcalibDB.h"
90 #include "AliTPCParam.h"
91 #include "THnBase.h"
92 #include <RVersion.h>
93 #if ROOT_VERSION_CODE < ROOT_VERSION(6,3,3)
94 #include "TFormulaPrimitive.h"
95 #else
96 #include "v5/TFormulaPrimitive.h"
97 #endif
98 #include "AliMathBase.h"
99 
103 
104 
105 AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
106 
107 
108 /*
109  Example usage fitting parameterization:
110  TFile fres("resol.root"); //tree with resolution and shape
111  TTree * treeRes =(TTree*)fres.Get("Resol");
112 
113  AliTPCClusterParam param;
114  param.SetInstance(&param);
115  param.FitResol(treeRes);
116  param.FitRMS(treeRes);
117  TFile fparam("TPCClusterParam.root","recreate");
118  param.Write("Param");
119  //
120  //
121  TFile fparam("TPCClusterParam.root");
122  AliTPCClusterParam *param2 = (AliTPCClusterParam *) fparam.Get("Param");
123  param2->SetInstance(param2);
124  param2->Test(treeRes);
125 
126 
127  treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
128 
129 */
130 
131 
132 
133 Double_t AliTPCClusterParamClusterResolutionY(Double_t x, Double_t s0, Double_t s1){
134  //
135  // cluster resoultion
136  Double_t result = AliTPCClusterParam::SGetError0Par(0,s0,x,s1);
137  return result;
138 }
139 
140 Double_t AliTPCClusterParamClusterResolutionZ(Double_t x, Double_t s0, Double_t s1){
141  //
142  // cluster resoultion
143  Double_t result = AliTPCClusterParam::SGetError0Par(1,s0,x,s1);
144  return result;
145 }
146 
147 
148 
149 //_ singleton implementation __________________________________________________
151 {
154 
155  if (fgInstance == 0){
157  //
158  ::Info("AliTPCClusterParam::Instance()","registering of the cluster param function primitives");
159  ::Info("AliTPCClusterParam::Instance()","clusterResolutionYAliRoot");
160 #if ROOT_VERSION_CODE < ROOT_VERSION(6,3,3)
161  TFormulaPrimitive::AddFormula(new TFormulaPrimitive("clusterResolutionYAliRoot","clusterResolutionYAliRoot",AliTPCClusterParamClusterResolutionY));
162 #else
163  ROOT::v5::TFormulaPrimitive::AddFormula(new ROOT::v5::TFormulaPrimitive("clusterResolutionYAliRoot","clusterResolutionYAliRoot",AliTPCClusterParamClusterResolutionY));
164 #endif
165  ::Info("AliTPCClusterParam::Instance()","clusterResolutionZAliRoot");
166 #if ROOT_VERSION_CODE < ROOT_VERSION(6,3,3)
167  TFormulaPrimitive::AddFormula(new TFormulaPrimitive("clusterResolutionZAliRoot","clusterResolutionZAliRoot",AliTPCClusterParamClusterResolutionZ));
168 #else
169  ROOT::v5::TFormulaPrimitive::AddFormula(new ROOT::v5::TFormulaPrimitive("clusterResolutionZAliRoot","clusterResolutionZAliRoot",AliTPCClusterParamClusterResolutionZ));
170 #endif
171 
172  }
173  return fgInstance;
174 }
175 
176 
178  TObject(),
179  fRatio(0),
180  fQNorm(0),
181  fQNormCorr(0),
182  fQNormHis(0),
183  fQpadTnorm(0), // q pad normalization - Total charge
184  fQpadMnorm(0), // q pad normalization - Max charge
185  fWaveCorrectionMap(0),
186  fWaveCorrectionMirroredPad( kFALSE ),
187  fWaveCorrectionMirroredZ( kFALSE ),
188  fWaveCorrectionMirroredAngle( kFALSE ),
189  fResolutionYMap(0)
190  //
191 {
193 
194  fPosQTnorm[0] = 0; fPosQTnorm[1] = 0; fPosQTnorm[2] = 0;
195  fPosQMnorm[0] = 0; fPosQMnorm[1] = 0; fPosQMnorm[2] = 0;
196  //
197  fPosYcor[0] = 0; fPosYcor[1] = 0; fPosYcor[2] = 0;
198  fPosZcor[0] = 0; fPosZcor[1] = 0; fPosZcor[2] = 0;
199  fErrorRMSSys[0]=0; fErrorRMSSys[1]=0;
200 }
201 
203  TObject(param),
204  fRatio(0),
205  fQNorm(0),
206  fQNormCorr(0),
207  fQNormHis(0),
208  fQpadTnorm(new TVectorD(*(param.fQpadTnorm))), // q pad normalization - Total charge
209  fQpadMnorm(new TVectorD(*(param.fQpadMnorm))), // q pad normalization - Max charge
210  fWaveCorrectionMap(0),
211  fWaveCorrectionMirroredPad( kFALSE ),
212  fWaveCorrectionMirroredZ( kFALSE ),
213  fWaveCorrectionMirroredAngle( kFALSE ),
214  fResolutionYMap(0)
215 {
217 
218  if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
219  if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
220  //
221  if (param.fPosQTnorm[0]){
222  fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
223  fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
224  fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
225  //
226  fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
227  fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
228  fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
229  }
230  if (param.fPosYcor[0]){
231  fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
232  fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
233  fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
234  //
235  fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
236  fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
237  fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
238  }
239 
240  for (Int_t ii = 0; ii < 2; ++ii) {
241  for (Int_t jj = 0; jj < 3; ++jj) {
242  for (Int_t kk = 0; kk < 4; ++kk) {
243  fParamS0[ii][jj][kk] = param.fParamS0[ii][jj][kk];
244  fErrorS0[ii][jj][kk] = param.fErrorS0[ii][jj][kk];
245  fParamRMS0[ii][jj][kk] = param.fParamRMS0[ii][jj][kk];
246  fErrorRMS0[ii][jj][kk] = param.fErrorRMS0[ii][jj][kk];
247  }
248  for (Int_t kk = 0; kk < 7; ++kk) {
249  fParamS0Par[ii][jj][kk] = param.fParamS0Par[ii][jj][kk];
250  fErrorS0Par[ii][jj][kk] = param.fErrorS0Par[ii][jj][kk];
251  }
252  for (Int_t kk = 0; kk < 6; ++kk) {
253  fParamSQ[ii][jj][kk] = param.fParamSQ[ii][jj][kk];
254  fErrorSQ[ii][jj][kk] = param.fErrorSQ[ii][jj][kk];
255  fParamRMSQ[ii][jj][kk] = param.fParamRMSQ[ii][jj][kk];
256  fErrorRMSQ[ii][jj][kk] = param.fErrorRMSQ[ii][jj][kk];
257  }
258  for (Int_t kk = 0; kk < 9; ++kk) {
259  fParamSQPar[ii][jj][kk] = param.fParamSQPar[ii][jj][kk];
260  fErrorSQPar[ii][jj][kk] = param.fErrorSQPar[ii][jj][kk];
261  }
262  for (Int_t kk = 0; kk < 2; ++kk) {
263  fRMSSigmaFit[ii][jj][kk] = param.fRMSSigmaFit[ii][jj][kk];
264  }
265  }
266  for (Int_t jj = 0; jj < 4; ++jj) {
267  fParamS1[ii][jj] = param.fParamS1[ii][jj];
268  fErrorS1[ii][jj] = param.fErrorS1[ii][jj];
269  }
270  for (Int_t jj = 0; jj < 5; ++jj) {
271  fParamRMS1[ii][jj] = param.fParamRMS1[ii][jj];
272  fErrorRMS1[ii][jj] = param.fErrorRMS1[ii][jj];
273  }
274  fErrorRMSSys[ii] = param.fErrorRMSSys[ii];
275  for (Int_t jj = 0; jj < 2; ++jj){
276  fRMSSigmaRatio[ii][jj] = param.fRMSSigmaRatio[ii][jj];
277  }
278  }
279 
282 }
283 
284 
287 
288  if (this != &param) {
289  if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
290  if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
291  if (param.fPosQTnorm[0]){
292  fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
293  fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
294  fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
295  //
296  fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
297  fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
298  fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
299  }
300  if (param.fPosYcor[0]){
301  fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
302  fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
303  fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
304  //
305  fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
306  fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
307  fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
308  }
309 
310  for (Int_t ii = 0; ii < 2; ++ii) {
311  for (Int_t jj = 0; jj < 3; ++jj) {
312  for (Int_t kk = 0; kk < 4; ++kk) {
313  fParamS0[ii][jj][kk] = param.fParamS0[ii][jj][kk];
314  fErrorS0[ii][jj][kk] = param.fErrorS0[ii][jj][kk];
315  fParamRMS0[ii][jj][kk] = param.fParamRMS0[ii][jj][kk];
316  fErrorRMS0[ii][jj][kk] = param.fErrorRMS0[ii][jj][kk];
317  }
318  for (Int_t kk = 0; kk < 7; ++kk) {
319  fParamS0Par[ii][jj][kk] = param.fParamS0Par[ii][jj][kk];
320  fErrorS0Par[ii][jj][kk] = param.fErrorS0Par[ii][jj][kk];
321  }
322  for (Int_t kk = 0; kk < 6; ++kk) {
323  fParamSQ[ii][jj][kk] = param.fParamSQ[ii][jj][kk];
324  fErrorSQ[ii][jj][kk] = param.fErrorSQ[ii][jj][kk];
325  fParamRMSQ[ii][jj][kk] = param.fParamRMSQ[ii][jj][kk];
326  fErrorRMSQ[ii][jj][kk] = param.fErrorRMSQ[ii][jj][kk];
327  }
328  for (Int_t kk = 0; kk < 9; ++kk) {
329  fParamSQPar[ii][jj][kk] = param.fParamSQPar[ii][jj][kk];
330  fErrorSQPar[ii][jj][kk] = param.fErrorSQPar[ii][jj][kk];
331  }
332  for (Int_t kk = 0; kk < 2; ++kk) {
333  fRMSSigmaFit[ii][jj][kk] = param.fRMSSigmaFit[ii][jj][kk];
334  }
335  }
336  for (Int_t jj = 0; jj < 4; ++jj) {
337  fParamS1[ii][jj] = param.fParamS1[ii][jj];
338  fErrorS1[ii][jj] = param.fErrorS1[ii][jj];
339  }
340  for (Int_t jj = 0; jj < 5; ++jj) {
341  fParamRMS1[ii][jj] = param.fParamRMS1[ii][jj];
342  fErrorRMS1[ii][jj] = param.fErrorRMS1[ii][jj];
343  }
344  fErrorRMSSys[ii] = param.fErrorRMSSys[ii];
345  for (Int_t jj = 0; jj < 2; ++jj){
346  fRMSSigmaRatio[ii][jj] = param.fRMSSigmaRatio[ii][jj];
347  }
348  }
349 
352  }
353  return *this;
354 }
355 
356 
359 
360  if (fQNorm) fQNorm->Delete();
361  if (fQNormCorr) delete fQNormCorr;
362  if (fQNormHis) fQNormHis->Delete();
363  delete fQNorm;
364  delete fQNormHis;
365  if (fPosQTnorm[0]){
366  delete fPosQTnorm[0];
367  delete fPosQTnorm[1];
368  delete fPosQTnorm[2];
369  //
370  delete fPosQMnorm[0];
371  delete fPosQMnorm[1];
372  delete fPosQMnorm[2];
373  }
374  if (fPosYcor[0]){
375  delete fPosYcor[0];
376  delete fPosYcor[1];
377  delete fPosYcor[2];
378  //
379  delete fPosZcor[0];
380  delete fPosZcor[1];
381  delete fPosZcor[2];
382  }
383  delete fWaveCorrectionMap;
384  delete fResolutionYMap;
385 }
386 
387 
388 void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
392 
393  TString varVal;
394  varVal="Resol:AngleM:Zm";
395  TString varErr;
396  varErr="Sigma:AngleS:Zs";
397  TString varCut;
398  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
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 
441 void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
445 
446  TString varVal;
447  varVal="Resol:AngleM:Zm";
448  TString varErr;
449  varErr="Sigma:AngleS:Zs";
450  TString varCut;
451  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
452  //
453  Int_t entries = tree->Draw(varVal.Data(),varCut);
454  Float_t px[10000], py[10000], pz[10000];
455  Float_t ex[10000], ey[10000], ez[10000];
456  //
457  tree->Draw(varErr.Data(),varCut);
458  for (Int_t ipoint=0; ipoint<entries; ipoint++){
459  ex[ipoint]= tree->GetV3()[ipoint];
460  ey[ipoint]= tree->GetV2()[ipoint];
461  ez[ipoint]= tree->GetV1()[ipoint];
462  }
463  tree->Draw(varVal.Data(),varCut);
464  for (Int_t ipoint=0; ipoint<entries; ipoint++){
465  px[ipoint]= tree->GetV3()[ipoint];
466  py[ipoint]= tree->GetV2()[ipoint];
467  pz[ipoint]= tree->GetV1()[ipoint];
468  }
469 
470  //
471  TLinearFitter fitter(6,"hyp5");
472  for (Int_t ipoint=0; ipoint<entries; ipoint++){
473  Float_t val = pz[ipoint]*pz[ipoint];
474  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
475  Double_t x[6];
476  x[0] = px[ipoint];
477  x[1] = py[ipoint]*py[ipoint];
478  x[2] = x[0]*x[0];
479  x[3] = x[1]*x[1];
480  x[4] = x[0]*x[1];
481  fitter.AddPoint(x,val,err);
482  }
483  fitter.Eval();
484  TVectorD param(6);
485  fitter.GetParameters(param);
486  param0[0] = param[0];
487  param0[1] = param[1];
488  param0[2] = param[2];
489  param0[3] = param[3];
490  param0[4] = param[4];
491  param0[5] = param[5];
492  Float_t chi2 = fitter.GetChisquare()/entries;
493  param0[6] = chi2;
494  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
495  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
496  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
497  error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
498  error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
499  error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
500 }
501 
502 
503 
504 
505 
506 void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
510 
511  TString varVal;
512  varVal="Resol:AngleM*sqrt(Length):Zm/Length";
513  TString varErr;
514  varErr="Sigma:AngleS:Zs";
515  TString varCut;
516  varCut=Form("Dim==%d&&QMean<0",dim);
517  //
518  Int_t entries = tree->Draw(varVal.Data(),varCut);
519  Float_t px[10000], py[10000], pz[10000];
520  Float_t ex[10000], ey[10000], ez[10000];
521  //
522  tree->Draw(varErr.Data(),varCut);
523  for (Int_t ipoint=0; ipoint<entries; ipoint++){
524  ex[ipoint]= tree->GetV3()[ipoint];
525  ey[ipoint]= tree->GetV2()[ipoint];
526  ez[ipoint]= tree->GetV1()[ipoint];
527  }
528  tree->Draw(varVal.Data(),varCut);
529  for (Int_t ipoint=0; ipoint<entries; ipoint++){
530  px[ipoint]= tree->GetV3()[ipoint];
531  py[ipoint]= tree->GetV2()[ipoint];
532  pz[ipoint]= tree->GetV1()[ipoint];
533  }
534 
535  //
536  TLinearFitter fitter(3,"hyp2");
537  for (Int_t ipoint=0; ipoint<entries; ipoint++){
538  Float_t val = pz[ipoint]*pz[ipoint];
539  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
540  Double_t x[2];
541  x[0] = px[ipoint];
542  x[1] = py[ipoint]*py[ipoint];
543  fitter.AddPoint(x,val,err);
544  }
545  fitter.Eval();
546  TVectorD param(3);
547  fitter.GetParameters(param);
548  param0[0] = param[0];
549  param0[1] = param[1];
550  param0[2] = param[2];
551  Float_t chi2 = fitter.GetChisquare()/entries;
552  param0[3] = chi2;
553  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
554  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
555  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
556 }
557 
558 void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
562 
563  TString varVal;
564  varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
565  char varVal0[100];
566  snprintf(varVal0,100,"Resol:AngleM:Zm");
567  //
568  TString varErr;
569  varErr="Sigma:AngleS:Zs";
570  TString varCut;
571  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
572  //
573  Int_t entries = tree->Draw(varVal.Data(),varCut);
574  Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
575  Float_t ex[20000], ey[20000], ez[20000];
576  //
577  tree->Draw(varErr.Data(),varCut);
578  for (Int_t ipoint=0; ipoint<entries; ipoint++){
579  ex[ipoint]= tree->GetV3()[ipoint];
580  ey[ipoint]= tree->GetV2()[ipoint];
581  ez[ipoint]= tree->GetV1()[ipoint];
582  }
583  tree->Draw(varVal.Data(),varCut);
584  for (Int_t ipoint=0; ipoint<entries; ipoint++){
585  px[ipoint]= tree->GetV3()[ipoint];
586  py[ipoint]= tree->GetV2()[ipoint];
587  pz[ipoint]= tree->GetV1()[ipoint];
588  }
589  tree->Draw(varVal0,varCut);
590  for (Int_t ipoint=0; ipoint<entries; ipoint++){
591  pu[ipoint]= tree->GetV3()[ipoint];
592  pt[ipoint]= tree->GetV2()[ipoint];
593  }
594 
595  //
596  TLinearFitter fitter(5,"hyp4");
597  for (Int_t ipoint=0; ipoint<entries; ipoint++){
598  Float_t val = pz[ipoint]*pz[ipoint];
599  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
600  Double_t x[4];
601  x[0] = pu[ipoint];
602  x[1] = pt[ipoint]*pt[ipoint];
603  x[2] = px[ipoint];
604  x[3] = py[ipoint]*py[ipoint];
605  fitter.AddPoint(x,val,err);
606  }
607 
608  fitter.Eval();
609  TVectorD param(5);
610  fitter.GetParameters(param);
611  param0[0] = param[0];
612  param0[1] = param[1];
613  param0[2] = param[2];
614  param0[3] = param[3];
615  param0[4] = param[4];
616  Float_t chi2 = fitter.GetChisquare()/entries;
617  param0[5] = chi2;
618  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
619  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
620  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
621  error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
622  error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
623 }
624 
625 void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
629 
630  TString varVal;
631  varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
632  char varVal0[100];
633  snprintf(varVal0,100,"Resol:AngleM:Zm");
634  //
635  TString varErr;
636  varErr="Sigma:AngleS:Zs";
637  TString varCut;
638  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
639  //
640  Int_t entries = tree->Draw(varVal.Data(),varCut);
641  Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
642  Float_t ex[20000], ey[20000], ez[20000];
643  //
644  tree->Draw(varErr.Data(),varCut);
645  for (Int_t ipoint=0; ipoint<entries; ipoint++){
646  ex[ipoint]= tree->GetV3()[ipoint];
647  ey[ipoint]= tree->GetV2()[ipoint];
648  ez[ipoint]= tree->GetV1()[ipoint];
649  }
650  tree->Draw(varVal.Data(),varCut);
651  for (Int_t ipoint=0; ipoint<entries; ipoint++){
652  px[ipoint]= tree->GetV3()[ipoint];
653  py[ipoint]= tree->GetV2()[ipoint];
654  pz[ipoint]= tree->GetV1()[ipoint];
655  }
656  tree->Draw(varVal0,varCut);
657  for (Int_t ipoint=0; ipoint<entries; ipoint++){
658  pu[ipoint]= tree->GetV3()[ipoint];
659  pt[ipoint]= tree->GetV2()[ipoint];
660  }
661 
662  //
663  TLinearFitter fitter(8,"hyp7");
664  for (Int_t ipoint=0; ipoint<entries; ipoint++){
665  Float_t val = pz[ipoint]*pz[ipoint];
666  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
667  Double_t x[7];
668  x[0] = pu[ipoint];
669  x[1] = pt[ipoint]*pt[ipoint];
670  x[2] = x[0]*x[0];
671  x[3] = x[1]*x[1];
672  x[4] = x[0]*x[1];
673  x[5] = px[ipoint];
674  x[6] = py[ipoint]*py[ipoint];
675  //
676  fitter.AddPoint(x,val,err);
677  }
678 
679  fitter.Eval();
680  TVectorD param(8);
681  fitter.GetParameters(param);
682  param0[0] = param[0];
683  param0[1] = param[1];
684  param0[2] = param[2];
685  param0[3] = param[3];
686  param0[4] = param[4];
687  param0[5] = param[5];
688  param0[6] = param[6];
689  param0[7] = param[7];
690 
691  Float_t chi2 = fitter.GetChisquare()/entries;
692  param0[8] = chi2;
693  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
694  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
695  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
696  error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
697  error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
698  error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
699  error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
700  error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
701 }
702 
703 
704 
705 void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
709 
710  TString varVal;
711  varVal="RMSm:AngleM:Zm";
712  TString varErr;
713  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
714  TString varCut;
715  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
716  //
717  Int_t entries = tree->Draw(varVal.Data(),varCut);
718  Float_t px[10000], py[10000], pz[10000];
719  Float_t ex[10000], ey[10000], ez[10000];
720  //
721  tree->Draw(varErr.Data(),varCut);
722  for (Int_t ipoint=0; ipoint<entries; ipoint++){
723  ex[ipoint]= tree->GetV3()[ipoint];
724  ey[ipoint]= tree->GetV2()[ipoint];
725  ez[ipoint]= tree->GetV1()[ipoint];
726  }
727  tree->Draw(varVal.Data(),varCut);
728  for (Int_t ipoint=0; ipoint<entries; ipoint++){
729  px[ipoint]= tree->GetV3()[ipoint];
730  py[ipoint]= tree->GetV2()[ipoint];
731  pz[ipoint]= tree->GetV1()[ipoint];
732  }
733 
734  //
735  TLinearFitter fitter(3,"hyp2");
736  for (Int_t ipoint=0; ipoint<entries; ipoint++){
737  Float_t val = pz[ipoint]*pz[ipoint];
738  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
739  Double_t x[2];
740  x[0] = px[ipoint];
741  x[1] = py[ipoint]*py[ipoint];
742  fitter.AddPoint(x,val,err);
743  }
744  fitter.Eval();
745  TVectorD param(3);
746  fitter.GetParameters(param);
747  param0[0] = param[0];
748  param0[1] = param[1];
749  param0[2] = param[2];
750  Float_t chi2 = fitter.GetChisquare()/entries;
751  param0[3] = chi2;
752  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
753  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
754  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
755 }
756 
757 void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
761 
762  TString varVal;
763  varVal="RMSm:AngleM*Length:Zm";
764  TString varErr;
765  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad";
766  TString varCut;
767  varCut=Form("Dim==%d&&QMean<0",dim);
768  //
769  Int_t entries = tree->Draw(varVal.Data(),varCut);
770  Float_t px[10000], py[10000], pz[10000];
771  Float_t type[10000], ey[10000], ez[10000];
772  //
773  tree->Draw(varErr.Data(),varCut);
774  for (Int_t ipoint=0; ipoint<entries; ipoint++){
775  type[ipoint] = tree->GetV3()[ipoint];
776  ey[ipoint] = tree->GetV2()[ipoint];
777  ez[ipoint] = tree->GetV1()[ipoint];
778  }
779  tree->Draw(varVal.Data(),varCut);
780  for (Int_t ipoint=0; ipoint<entries; ipoint++){
781  px[ipoint]= tree->GetV3()[ipoint];
782  py[ipoint]= tree->GetV2()[ipoint];
783  pz[ipoint]= tree->GetV1()[ipoint];
784  }
785 
786  //
787  TLinearFitter fitter(4,"hyp3");
788  for (Int_t ipoint=0; ipoint<entries; ipoint++){
789  Float_t val = pz[ipoint]*pz[ipoint];
790  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
791  Double_t x[3];
792  x[0] = (type[ipoint]<0.5)? 0.:1.;
793  x[1] = px[ipoint];
794  x[2] = py[ipoint]*py[ipoint];
795  fitter.AddPoint(x,val,err);
796  }
797  fitter.Eval();
798  TVectorD param(4);
799  fitter.GetParameters(param);
800  param0[0] = param[0];
801  param0[1] = param[0]+param[1];
802  param0[2] = param[2];
803  param0[3] = param[3];
804  Float_t chi2 = fitter.GetChisquare()/entries;
805  param0[4] = chi2;
806  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
807  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
808  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
809  error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
810 }
811 
812 void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
816 
817  TString varVal;
818  varVal="RMSm:AngleM/sqrt(QMean):Zm/QMean";
819  char varVal0[100];
820  snprintf(varVal0,100,"RMSm:AngleM:Zm");
821  //
822  TString varErr;
823  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
824  TString varCut;
825  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
826  //
827  Int_t entries = tree->Draw(varVal.Data(),varCut);
828  Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
829  Float_t ex[20000], ey[20000], ez[20000];
830  //
831  tree->Draw(varErr.Data(),varCut);
832  for (Int_t ipoint=0; ipoint<entries; ipoint++){
833  ex[ipoint]= tree->GetV3()[ipoint];
834  ey[ipoint]= tree->GetV2()[ipoint];
835  ez[ipoint]= tree->GetV1()[ipoint];
836  }
837  tree->Draw(varVal.Data(),varCut);
838  for (Int_t ipoint=0; ipoint<entries; ipoint++){
839  px[ipoint]= tree->GetV3()[ipoint];
840  py[ipoint]= tree->GetV2()[ipoint];
841  pz[ipoint]= tree->GetV1()[ipoint];
842  }
843  tree->Draw(varVal0,varCut);
844  for (Int_t ipoint=0; ipoint<entries; ipoint++){
845  pu[ipoint]= tree->GetV3()[ipoint];
846  pt[ipoint]= tree->GetV2()[ipoint];
847  }
848 
849  //
850  TLinearFitter fitter(5,"hyp4");
851  for (Int_t ipoint=0; ipoint<entries; ipoint++){
852  Float_t val = pz[ipoint]*pz[ipoint];
853  Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
854  Double_t x[4];
855  x[0] = pu[ipoint];
856  x[1] = pt[ipoint]*pt[ipoint];
857  x[2] = px[ipoint];
858  x[3] = py[ipoint]*py[ipoint];
859  fitter.AddPoint(x,val,err);
860  }
861 
862  fitter.Eval();
863  TVectorD param(5);
864  fitter.GetParameters(param);
865  param0[0] = param[0];
866  param0[1] = param[1];
867  param0[2] = param[2];
868  param0[3] = param[3];
869  param0[4] = param[4];
870  Float_t chi2 = fitter.GetChisquare()/entries;
871  param0[5] = chi2;
872  error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
873  error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
874  error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
875  error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
876  error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
877 }
878 
879 
880 void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
884 
885  TString varVal;
886  varVal="RMSs:RMSm";
887  //
888  TString varCut;
889  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
890  //
891  Int_t entries = tree->Draw(varVal.Data(),varCut);
892  Float_t px[20000], py[20000];
893  //
894  tree->Draw(varVal.Data(),varCut);
895  for (Int_t ipoint=0; ipoint<entries; ipoint++){
896  px[ipoint]= tree->GetV2()[ipoint];
897  py[ipoint]= tree->GetV1()[ipoint];
898  }
899  TLinearFitter fitter(2,"pol1");
900  for (Int_t ipoint=0; ipoint<entries; ipoint++){
901  Float_t val = py[ipoint];
902  Float_t err = fRatio*px[ipoint];
903  Double_t x[4];
904  x[0] = px[ipoint];
905  if (err>0) fitter.AddPoint(x,val,err);
906  }
907  fitter.Eval();
908  param0[0]= fitter.GetParameter(0);
909  param0[1]= fitter.GetParameter(1);
910 }
911 
912 
913 
914 Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
916 
917  Float_t value=0;
918  value += fParamS0[dim][type][0];
919  value += fParamS0[dim][type][1]*z;
920  value += fParamS0[dim][type][2]*angle*angle;
921  value = TMath::Sqrt(TMath::Abs(value));
922  return value;
923 }
924 
925 
926 Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
928 
929  Float_t value=0;
930  value += fParamS0Par[dim][type][0];
931  value += fParamS0Par[dim][type][1]*z;
932  value += fParamS0Par[dim][type][2]*angle*angle;
933  value += fParamS0Par[dim][type][3]*z*z;
934  value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
935  value += fParamS0Par[dim][type][5]*z*angle*angle;
936  value = TMath::Sqrt(TMath::Abs(value));
937  return value;
938 }
939 
940 
941 
942 Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
944 
945  Float_t value=0;
946  Float_t length=0.75;
947  if (type==1) length=1;
948  if (type==2) length=1.5;
949  value += fParamS1[dim][0];
950  value += fParamS1[dim][1]*z/length;
951  value += fParamS1[dim][2]*angle*angle*length;
952  value = TMath::Sqrt(TMath::Abs(value));
953  return value;
954 }
955 
956 Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
958 
959  Float_t value=0;
960  value += fParamSQ[dim][type][0];
961  value += fParamSQ[dim][type][1]*z;
962  value += fParamSQ[dim][type][2]*angle*angle;
963  value += fParamSQ[dim][type][3]*z/Qmean;
964  value += fParamSQ[dim][type][4]*angle*angle/Qmean;
965  value = TMath::Sqrt(TMath::Abs(value));
966  return value;
967 
968 
969 }
970 
971 Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
973 
974  Float_t value=0;
975  value += fParamSQPar[dim][type][0];
976  value += fParamSQPar[dim][type][1]*z;
977  value += fParamSQPar[dim][type][2]*angle*angle;
978  value += fParamSQPar[dim][type][3]*z*z;
979  value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
980  value += fParamSQPar[dim][type][5]*z*angle*angle;
981  value += fParamSQPar[dim][type][6]*z/Qmean;
982  value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
983  value = TMath::Sqrt(TMath::Abs(value));
984  return value;
985 
986 
987 }
988 
989 Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
991 
992  Float_t value=0;
993  value += fParamSQPar[dim][type][0];
994  value += fParamSQPar[dim][type][1]*z;
995  value += fParamSQPar[dim][type][2]*angle*angle;
996  value += fParamSQPar[dim][type][3]*z*z;
997  value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
998  value += fParamSQPar[dim][type][5]*z*angle*angle;
999  value += fParamSQPar[dim][type][6]*z/Qmean;
1000  value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
1001  Float_t valueMean = GetError0Par(dim,type,z,angle);
1002  value -= 0.35*0.35*valueMean*valueMean;
1003  value = TMath::Sqrt(TMath::Abs(value));
1004  return value;
1005 
1006 
1007 }
1008 
1009 Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
1011 
1012  Float_t value=0;
1013  value += fParamRMS0[dim][type][0];
1014  value += fParamRMS0[dim][type][1]*z;
1015  value += fParamRMS0[dim][type][2]*angle*angle;
1016  value = TMath::Sqrt(TMath::Abs(value));
1017  return value;
1018 }
1019 
1020 Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
1022 
1023  Float_t value=0;
1024  Float_t length=0.75;
1025  if (type==1) length=1;
1026  if (type==2) length=1.5;
1027  if (type==0){
1028  value += fParamRMS1[dim][0];
1029  }else{
1030  value += fParamRMS1[dim][1];
1031  }
1032  value += fParamRMS1[dim][2]*z;
1033  value += fParamRMS1[dim][3]*angle*angle*length*length;
1034  value = TMath::Sqrt(TMath::Abs(value));
1035  return value;
1036 }
1037 
1038 Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
1040 
1041  Float_t value=0;
1042  value += fParamRMSQ[dim][type][0];
1043  value += fParamRMSQ[dim][type][1]*z;
1044  value += fParamRMSQ[dim][type][2]*angle*angle;
1045  value += fParamRMSQ[dim][type][3]*z/Qmean;
1046  value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
1047  value = TMath::Sqrt(TMath::Abs(value));
1048  return value;
1049 }
1050 
1051 Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
1053 
1054  Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
1055  Float_t value = fRMSSigmaFit[dim][type][0];
1056  value+= fRMSSigmaFit[dim][type][1]*mean;
1057  return value;
1058 }
1059 
1060 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 {
1062 
1063  Double_t value =0;
1064  //
1065  Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
1066  if (rmsL<rmsMeanQ) return value;
1067  //
1068  Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
1069  //
1070  if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
1071  //1.5 sigma cut on mean
1072  value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
1073  }else{
1074  if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
1075  //3 sigma cut on local
1076  value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
1077  }
1078  }
1079  return TMath::Sqrt(TMath::Abs(value));
1080 }
1081 
1082 
1083 
1086 
1087  FitResol(tree);
1088  FitRMS(tree);
1089 
1090 }
1091 
1094 
1095  SetInstance(this);
1096  for (Int_t idir=0;idir<2; idir++){
1097  for (Int_t itype=0; itype<3; itype++){
1098  Float_t param0[10];
1099  Float_t error0[10];
1100  // model error param
1101  FitResol0(tree, idir, itype,param0,error0);
1102  printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1103  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1104  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1105  for (Int_t ipar=0;ipar<4; ipar++){
1106  fParamS0[idir][itype][ipar] = param0[ipar];
1107  fErrorS0[idir][itype][ipar] = param0[ipar];
1108  }
1109  // error param with parabolic correction
1110  FitResol0Par(tree, idir, itype,param0,error0);
1111  printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
1112  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]);
1113  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]);
1114  for (Int_t ipar=0;ipar<7; ipar++){
1115  fParamS0Par[idir][itype][ipar] = param0[ipar];
1116  fErrorS0Par[idir][itype][ipar] = param0[ipar];
1117  }
1118  //
1119  FitResolQ(tree, idir, itype,param0,error0);
1120  printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1121  printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1122  printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1123  for (Int_t ipar=0;ipar<6; ipar++){
1124  fParamSQ[idir][itype][ipar] = param0[ipar];
1125  fErrorSQ[idir][itype][ipar] = param0[ipar];
1126  }
1127  //
1128  FitResolQPar(tree, idir, itype,param0,error0);
1129  printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
1130  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]);
1131  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]);
1132  for (Int_t ipar=0;ipar<9; ipar++){
1133  fParamSQPar[idir][itype][ipar] = param0[ipar];
1134  fErrorSQPar[idir][itype][ipar] = param0[ipar];
1135  }
1136  }
1137  }
1138  //
1139  printf("Resol z-scaled\n");
1140  for (Int_t idir=0;idir<2; idir++){
1141  Float_t param0[4];
1142  Float_t error0[4];
1143  FitResol1(tree, idir,param0,error0);
1144  printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
1145  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1146  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1147  for (Int_t ipar=0;ipar<4; ipar++){
1148  fParamS1[idir][ipar] = param0[ipar];
1149  fErrorS1[idir][ipar] = param0[ipar];
1150  }
1151  }
1152 
1153  for (Int_t idir=0;idir<2; idir++){
1154  printf("\nDirection %d\n",idir);
1155  printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1156  for (Int_t itype=0; itype<3; itype++){
1157  Float_t length=0.75;
1158  if (itype==1) length=1;
1159  if (itype==2) length=1.5;
1160  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));
1161  }
1162  }
1163 }
1164 
1165 
1166 
1169 
1170  SetInstance(this);
1171  for (Int_t idir=0;idir<2; idir++){
1172  for (Int_t itype=0; itype<3; itype++){
1173  Float_t param0[6];
1174  Float_t error0[6];
1175  FitRMS0(tree, idir, itype,param0,error0);
1176  printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1177  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1178  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1179  for (Int_t ipar=0;ipar<4; ipar++){
1180  fParamRMS0[idir][itype][ipar] = param0[ipar];
1181  fErrorRMS0[idir][itype][ipar] = param0[ipar];
1182  }
1183  FitRMSQ(tree, idir, itype,param0,error0);
1184  printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1185  printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1186  printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1187  for (Int_t ipar=0;ipar<6; ipar++){
1188  fParamRMSQ[idir][itype][ipar] = param0[ipar];
1189  fErrorRMSQ[idir][itype][ipar] = param0[ipar];
1190  }
1191  }
1192  }
1193  //
1194  printf("RMS z-scaled\n");
1195  for (Int_t idir=0;idir<2; idir++){
1196  Float_t param0[5];
1197  Float_t error0[5];
1198  FitRMS1(tree, idir,param0,error0);
1199  printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1200  printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1201  printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1202  for (Int_t ipar=0;ipar<5; ipar++){
1203  fParamRMS1[idir][ipar] = param0[ipar];
1204  fErrorRMS1[idir][ipar] = param0[ipar];
1205  }
1206  }
1207 
1208  for (Int_t idir=0;idir<2; idir++){
1209  printf("\nDirection %d\n",idir);
1210  printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1211  for (Int_t itype=0; itype<3; itype++){
1212  Float_t length=0.75;
1213  if (itype==1) length=1;
1214  if (itype==2) length=1.5;
1215  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);
1216  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);
1217  }
1218  }
1219  //
1220  // Fit RMS sigma
1221  //
1222  printf("RMS fluctuation parameterization \n");
1223  for (Int_t idir=0;idir<2; idir++){
1224  for (Int_t itype=0; itype<3; itype++){
1225  Float_t param0[5];
1226  Float_t error0[5];
1227  FitRMSSigma(tree, idir,itype,param0,error0);
1228  printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1229  for (Int_t ipar=0;ipar<2; ipar++){
1230  fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1231  }
1232  }
1233  }
1234  //
1235  // store systematic error end RMS fluctuation parameterization
1236  //
1237  TH1F hratio("hratio","hratio",100,-0.1,0.1);
1238  tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1239  fErrorRMSSys[0] = hratio.GetRMS();
1240  tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1241  fErrorRMSSys[1] = hratio.GetRMS();
1242  TH1F hratioR("hratioR","hratioR",100,0,0.2);
1243  tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1244  fRMSSigmaRatio[0][0]=hratioR.GetMean();
1245  fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1246  tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1247  fRMSSigmaRatio[1][0]=hratioR.GetMean();
1248  fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1249 }
1250 
1251 void AliTPCClusterParam::Test(TTree * tree, const char *output){
1253 
1254  SetInstance(this);
1255  TFile f(output,"recreate");
1256  f.cd();
1257  //
1258  // 1D histograms - resolution
1259  //
1260  for (Int_t idim=0; idim<2; idim++){
1261  for (Int_t ipad=0; ipad<3; ipad++){
1262  char hname1[300];
1263  char hcut1[300];
1264  char hexp1[300];
1265  //
1266  snprintf(hname1,300,"Delta0 Dir %d Pad %d",idim,ipad);
1267  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1268  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1269  TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1270  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1271  tree->Draw(hexp1,hcut1,"");
1272  his1DRel0.Write();
1273  //
1274  snprintf(hname1,300,"Delta0Par Dir %d Pad %d",idim,ipad);
1275  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1276  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1277  TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1278  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1279  tree->Draw(hexp1,hcut1,"");
1280  his1DRel0Par.Write();
1281  //
1282  }
1283  }
1284  //
1285  // 2D histograms - resolution
1286  //
1287  for (Int_t idim=0; idim<2; idim++){
1288  for (Int_t ipad=0; ipad<3; ipad++){
1289  char hname1[300];
1290  char hcut1[300];
1291  char hexp1[300];
1292  //
1293  snprintf(hname1,300,"2DDelta0 Dir %d Pad %d",idim,ipad);
1294  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1295  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1296  TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1297  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1298  tree->Draw(hexp1,hcut1,"");
1299  profDRel0.Write();
1300  //
1301  snprintf(hname1,300,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1302  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1303  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1304  TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1305  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1306  tree->Draw(hexp1,hcut1,"");
1307  profDRel0Par.Write();
1308  //
1309  }
1310  }
1311 }
1312 
1313 
1314 
1315 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1317 
1318  //
1319  // Error parameterization
1320  //
1321  printf("\nResolution Scaled factors\n");
1322  printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1323  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])),
1324  TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
1325  for (Int_t ipad=0; ipad<3; ipad++){
1326  Float_t length=0.75;
1327  if (ipad==1) length=1;
1328  if (ipad==2) length=1.5;
1329  printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1330  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1331  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1332  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1333  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
1334  }
1335  for (Int_t ipad=0; ipad<3; ipad++){
1336  Float_t length=0.75;
1337  if (ipad==1) length=1;
1338  if (ipad==2) length=1.5;
1339  printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1340  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1341  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1342  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1343  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
1344  }
1345  printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1346  TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1347 
1348  for (Int_t ipad=0; ipad<3; ipad++){
1349  Float_t length=0.75;
1350  if (ipad==1) length=1;
1351  if (ipad==2) length=1.5;
1352  printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1353  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1354  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1355  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1356  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
1357  }
1358  for (Int_t ipad=0; ipad<3; ipad++){
1359  Float_t length=0.75;
1360  if (ipad==1) length=1;
1361  if (ipad==2) length=1.5;
1362  printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1363  TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1364  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1365  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1366  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
1367  }
1368 
1369  //
1370  // RMS scaling
1371  //
1372  printf("\n");
1373  printf("\nRMS Scaled factors\n");
1374  printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1375  printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1376  TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1377  TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1378  TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1379  TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1380  TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
1381  for (Int_t ipad=0; ipad<3; ipad++){
1382  Float_t length=0.75;
1383  if (ipad==1) length=1;
1384  if (ipad==2) length=1.5;
1385  if (ipad==0){
1386  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1387  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1388  0.,
1389  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1390  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1391  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1392  }else{
1393  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1394  0.,
1395  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1396  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1397  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1398  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1399  }
1400  }
1401  printf("\n");
1402  printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1403  TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1404  TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1405  TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1406  TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1407  TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
1408  for (Int_t ipad=0; ipad<3; ipad++){
1409  Float_t length=0.75;
1410  if (ipad==1) length=1;
1411  if (ipad==2) length=1.5;
1412  if (ipad==0){
1413  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1414  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1415  0.,
1416  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1417  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1418  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1419  }else{
1420  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1421  0.,
1422  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1423  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1424  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1425  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1426  }
1427  }
1428  printf("\ndEdx correction matrix used in GetQnormCorr\n");
1429  fQNormCorr->Print();
1430 
1431 }
1432 
1433 
1434 
1435 
1436 
1437 Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1443 
1444  if (fQNorm==0) return 0;
1445  TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1446  if (!norm) return 0;
1447  TVectorD &no = *norm;
1448  Float_t res =
1449  no[0]+
1450  no[1]*dr+
1451  no[2]*ty+
1452  no[3]*tz+
1453  no[4]*dr*ty+
1454  no[5]*dr*tz+
1455  no[6]*ty*tz+
1456  no[7]*dr*dr+
1457  no[8]*ty*ty+
1458  no[9]*tz*tz;
1459  res/=no[0];
1460  return res;
1461 }
1462 
1463 
1464 
1465 Float_t AliTPCClusterParam::QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t p2, Float_t p3){
1469 
1470  if (fQNormHis==0) return 0;
1471  TH3F * norm = (TH3F*)fQNormHis->At(4*itype+ipad);
1472  if (!norm) return 1;
1473  p2=TMath::Abs(p2);
1474  dr=TMath::Min(dr,Float_t(norm->GetXaxis()->GetXmax()-norm->GetXaxis()->GetBinWidth(0)));
1475  dr=TMath::Max(dr,Float_t(norm->GetXaxis()->GetXmin()+norm->GetXaxis()->GetBinWidth(0)));
1476  //
1477  p2=TMath::Min(p2,Float_t(norm->GetYaxis()->GetXmax()-norm->GetYaxis()->GetBinWidth(0)));
1478  p2=TMath::Max(p2,Float_t(norm->GetYaxis()->GetXmin()+norm->GetYaxis()->GetBinWidth(0)));
1479  //
1480  p3=TMath::Min(p3,Float_t(norm->GetZaxis()->GetXmax()-norm->GetZaxis()->GetBinWidth(0)));
1481  p3=TMath::Max(p3,Float_t(norm->GetZaxis()->GetXmin()+norm->GetZaxis()->GetBinWidth(0)));
1482  //
1483  Double_t res = norm->GetBinContent(norm->FindBin(dr,p2,p3));
1484  if (res==0) res = norm->GetBinContent(norm->FindBin(0.5,0.5,0.5)); // This is just hack - to be fixed entries without
1485 
1486  return res;
1487 }
1488 
1489 
1490 
1491 void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm){
1496 
1497  if (fQNorm==0) fQNorm = new TObjArray(6);
1498  fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1499 }
1500 
1503 
1504  if (!fQNormCorr) fQNormCorr= new TMatrixD(12,6);
1505  for (Int_t irow=0;irow<12; irow++)
1506  for (Int_t icol=0;icol<6; icol++){
1507  (*fQNormCorr)(irow,icol)=1.; // default - no correction
1508  if (irow>5) (*fQNormCorr)(irow,icol)=0.; // default - no correction
1509  }
1510 }
1511 
1512 void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val, Int_t mode){
1521 
1522  if (!fQNormCorr) {
1523  ResetQnormCorr();
1524  }
1525  //
1526  // eff shap parameterization matrix
1527  //
1528  // rows
1529  // itype*3+ipad - itype=0 qtot itype=1 qmax ipad=0
1530  //
1531  if (mode==1) {
1532  // mode introduced in 20.07.2014 - git describe ~ vAN-20140703-48-g3449a97 - to keep back compatibility with o
1533  (*fQNormCorr)(itype*3+ipad, corrType) = val; // set value
1534  return;
1535  }
1536  if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val; // multiplicative correction
1537  if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val; // additive correction
1538 }
1539 
1540 Double_t AliTPCClusterParam::GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const{
1542 
1543  if (!fQNormCorr) return 0;
1544  return (*fQNormCorr)(itype*3+ipad, corrType);
1545 }
1546 
1547 
1548 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){
1558 
1559  //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
1560  // Following variable used - correspondance to the our variable conventions
1561  //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1562  Double_t dp = ((pad-int(pad)-0.5)*2.);
1563  //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1564  Double_t dt = ((time-int(time)-0.5)*2.);
1565  //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1566  Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1567  //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1568  Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1569  //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1570  Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1571  //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1572  Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1573  //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1574  Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1575  //
1576  //
1577  //
1578  TVectorD * pvec = 0;
1579  if (isMax){
1580  pvec = fPosQMnorm[ipad];
1581  }else{
1582  pvec = fPosQTnorm[ipad];
1583  }
1584  TVectorD &param = *pvec;
1585  //
1586  // Eval part - in correspondance with fit part from debug streamer
1587  //
1588  Double_t result=param[0];
1589  Int_t index =1;
1590  //
1591  result+=dp*param[index++]; //1
1592  result+=dt*param[index++]; //2
1593  result+=dp*dp*param[index++]; //3
1594  result+=dt*dt*param[index++]; //4
1595  result+=dt*dt*dt*param[index++]; //5
1596  result+=dp*dt*param[index++]; //6
1597  result+=dp*dt*dt*param[index++]; //7
1598  result+=(dq0)*param[index++]; //8
1599  result+=(dq1)*param[index++]; //9
1600  //
1601  //
1602  result+=dp*dp*(di)*param[index++]; //10
1603  result+=dt*dt*(di)*param[index++]; //11
1604  result+=dp*dp*sy*param[index++]; //12
1605  result+=dt*sz*param[index++]; //13
1606  result+=dt*dt*sz*param[index++]; //14
1607  result+=dt*dt*dt*sz*param[index++]; //15
1608  //
1609  result+=dp*dp*1*sy*sz*param[index++]; //16
1610  result+=dt*sy*sz*param[index++]; //17
1611  result+=dt*dt*sy*sz*param[index++]; //18
1612  result+=dt*dt*dt*sy*sz*param[index++]; //19
1613  //
1614  result+=dp*dp*(dq0)*param[index++]; //20
1615  result+=dt*1*(dq0)*param[index++]; //21
1616  result+=dt*dt*(dq0)*param[index++]; //22
1617  result+=dt*dt*dt*(dq0)*param[index++]; //23
1618  //
1619  result+=dp*dp*(dq1)*param[index++]; //24
1620  result+=dt*(dq1)*param[index++]; //25
1621  result+=dt*dt*(dq1)*param[index++]; //26
1622  result+=dt*dt*dt*(dq1)*param[index++]; //27
1623 
1624  if (result<0.75) result=0.75;
1625  if (result>1.25) result=1.25;
1626 
1627  return result;
1628 
1629 }
1630 
1631 
1632 
1633 
1634 
1635 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*/){
1636 
1644 
1645  //
1646  //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
1647  //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
1648  //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
1649  //chainres->SetAlias("st","(sin(dt)-dt)");
1650  //
1651  //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1652 
1653  //
1654  // Derived variables
1655  //
1656  Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
1657  Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
1658  Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
1659  Double_t st = (TMath::Sin(dt)-dt);
1660  //
1661  Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.)));
1662  //
1663  //
1664  //
1665  TVectorD * pvec = 0;
1666  if (type==0){
1667  pvec = fPosYcor[ipad];
1668  }else{
1669  pvec = fPosZcor[ipad];
1670  }
1671  TVectorD &param = *pvec;
1672  //
1673  Double_t result=0;
1674  Int_t index =1;
1675 
1676  if (type==0){
1677  // y corr
1678  result+=(dp)*param[index++]; //1
1679  result+=(dp)*di*param[index++]; //2
1680  //
1681  result+=(sp)*param[index++]; //3
1682  result+=(sp)*di*param[index++]; //4
1683  }
1684  if (type==1){
1685  result+=(dt)*param[index++]; //1
1686  result+=(dt)*di*param[index++]; //2
1687  //
1688  result+=(st)*param[index++]; //3
1689  result+=(st)*di*param[index++]; //4
1690  }
1691  if (TMath::Abs(result)>0.05) return 0;
1692  return result;
1693 }
1694 
1695 
1696 
1697 Double_t AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
1704 
1705  const Double_t kEpsilon = 0.0001;
1706  const Double_t twoPi = TMath::TwoPi();
1707  const Double_t hnorm = 0.5/TMath::Sqrt(twoPi);
1708  const Double_t sqtwo = TMath::Sqrt(2.);
1709 
1710  if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
1711  // small angular effect
1712  Double_t val = TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1)/(s0*s1*twoPi);
1713  return val;
1714  }
1715  Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
1716  Double_t sigma = TMath::Sqrt(sigma2);
1717  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2.*sigma2));
1718  //
1719  Double_t sigmaErf = 1./(2.*s0*s1*sqtwo*sigma);
1720  Double_t k0s1s1 = 2.*k0*s1*s1;
1721  Double_t k1s0s0 = 2.*k1*s0*s0;
1722  Double_t erf0 = AliMathBase::ErfFast((sigma2-k0s1s1*x0-k1s0s0*x1)*sigmaErf);
1723  Double_t erf1 = AliMathBase::ErfFast((sigma2+k0s1s1*x0+k1s0s0*x1)*sigmaErf);
1724  Double_t norm = hnorm/sigma;
1725  Double_t val = norm*exp0*(erf0+erf1);
1726  return val;
1727 }
1728 
1729 
1730 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){
1737 
1738  Double_t sum =1, mean=0;
1739  // the COG of exponent
1740  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1741  mean+=iexp*TMath::Exp(-iexp/tau);
1742  sum +=TMath::Exp(-iexp/tau);
1743  }
1744  mean/=sum;
1745  //
1746  sum = 1;
1747  Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
1748  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1749  val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
1750  sum+=TMath::Exp(-iexp/tau);
1751  }
1752  return val/sum;
1753 }
1754 
1755 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){
1763 
1764  Double_t sum =0, mean=0;
1765  // the COG of G4
1766  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1767  Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1768  mean+=iexp*g4;
1769  sum +=g4;
1770  }
1771  mean/=sum;
1772  //
1773  sum = 0;
1774  Double_t val = 0;
1775  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1776  Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1777  val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
1778  sum+=g4;
1779  }
1780  return val/sum;
1781 }
1782 
1783 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){
1792 
1793  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1794  //
1795  // Gaus width sy and sz is determined by RF width and diffusion
1796  // Integral of Q is equal 1
1797  // Q max is calculated at position cpad, ctime
1798  // Example function:
1799  // TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000)
1800  //
1802  Double_t padLength= param->GetPadPitchLength(sector,row);
1803  Double_t padWidth = param->GetPadPitchWidth(sector);
1804  Double_t zwidth = param->GetZWidth();
1805  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1806 
1807  // diffusion in pad, time bin units
1808  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1809  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1810  diffT*=effDiff; //
1811  diffL*=effDiff; //
1812  //
1813  // transform angular effect to pad units
1814  //
1815  Double_t pky = ky*effLength/padWidth;
1816  Double_t pkz = kz*effLength/zwidth;
1817  // position in pad unit
1818  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1819  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1820  //
1821  //
1822  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1823  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1824  //return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
1825  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1826  return GaussConvolution(py,pz, pky,pkz,sy,sz)*length;
1827 }
1828 
1829 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){
1840 
1841  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1842  //
1843  // Gaus width sy and sz is determined by RF width and diffusion
1844  // Integral of Q is equal 1
1845  // Q max is calculated at position cpad, ctime
1846  //
1847  //
1848  //
1850  Double_t padLength= param->GetPadPitchLength(sector,row);
1851  Double_t padWidth = param->GetPadPitchWidth(sector);
1852  Double_t zwidth = param->GetZWidth();
1853  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1854  //
1855  // diffusion in pad units
1856  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1857  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1858  diffT*=effDiff; //
1859  diffL*=effDiff; //
1860  //
1861  // transform angular effect to pad units
1862  Double_t pky = ky*effLength/padWidth;
1863  Double_t pkz = kz*effLength/zwidth;
1864  // position in pad unit
1865  //
1866  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1867  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1868  //
1869  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1870  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1871  //
1872  //
1873  //
1874  Double_t sumAll=0,sumThr=0;
1875  //
1876  Double_t corr =1;
1877  Double_t qnorm=qtot;
1878  for (Float_t iy=-3;iy<=3;iy+=1.)
1879  for (Float_t iz=-4;iz<=4;iz+=1.){
1880  // Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);
1881  Double_t val = GaussConvolution(py-iy,pz-iz, pky,pkz, sy,sz);
1882  Double_t qlocal =qnorm*val;
1883  if (TMath::Abs(iy)<1.5&&TMath::Abs(iz)<1.5){
1884  sumThr+=qlocal; // Virtual charge used in cluster finder
1885  }
1886  else{
1887  if (qlocal>thr && TMath::Abs(iz)<2.5&&TMath::Abs(iy)<2.5) sumThr+=qlocal;
1888  }
1889  sumAll+=qlocal;
1890  }
1891  if (sumAll>0&&sumThr>0) {
1892  corr=(sumThr)/sumAll;
1893  }
1894  //
1895  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1896  return corr*length;
1897 }
1898 
1899 
1900 
1902 {
1904 
1905  delete fWaveCorrectionMap;
1906  fWaveCorrectionMap = 0;
1907  fWaveCorrectionMirroredPad = kFALSE;
1908  fWaveCorrectionMirroredZ = kFALSE;
1910  if( Map ){
1911  fWaveCorrectionMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
1912  if( fWaveCorrectionMap ){
1913  fWaveCorrectionMirroredPad = ( fWaveCorrectionMap->GetAxis(3)->FindFixBin(0.5)<=1 ); // Pad axis is mirrored at 0.5
1914  fWaveCorrectionMirroredZ = ( fWaveCorrectionMap->GetAxis(1)->FindFixBin(0.0)<=1); // Z axis is mirrored at 0
1915  fWaveCorrectionMirroredAngle = ( fWaveCorrectionMap->GetAxis(4)->FindFixBin(0.0)<=1 ); // Angle axis is mirrored at 0
1916  }
1917  }
1918 }
1919 
1921 {
1923 
1924  delete fResolutionYMap;
1925  fResolutionYMap = 0;
1926  if( Map ){
1927  fResolutionYMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
1928  }
1929 }
1930 
1931 Float_t AliTPCClusterParam::GetWaveCorrection(Int_t Type, Float_t Z, Int_t QMax, Float_t Pad, Float_t angleY ) const
1932 {
1934 
1935  if( !fWaveCorrectionMap ) return 0;
1936  Bool_t swapY = kFALSE;
1937  Pad = Pad-(Int_t)Pad;
1938 
1939  if( TMath::Abs(Pad-0.5)<1.e-8 ){// one pad clusters a stored in underflow bins
1940  Pad = -1.;
1941  } else {
1942  if( fWaveCorrectionMirroredPad && (Pad<0.5) ){ // cog axis is mirrored at 0.5
1943  swapY = !swapY;
1944  Pad = 1.0 - Pad;
1945  }
1946  }
1947 
1948  if( fWaveCorrectionMirroredZ && (Z<0) ){ // Z axis is mirrored at 0
1949  swapY = !swapY;
1950  Z = -Z;
1951  }
1952 
1953  if( fWaveCorrectionMirroredAngle && (angleY<0) ){ // Angle axis is mirrored at 0
1954  angleY = -angleY;
1955  }
1956  double var[5] = { static_cast<double>(Type), Z, static_cast<double>(QMax), Pad, angleY };
1957  Long64_t bin = fWaveCorrectionMap->GetBin(var, kFALSE );
1958  if( bin<0 ) return 0;
1959  Double_t dY = fWaveCorrectionMap->GetBinContent(bin);
1960  return (swapY ?-dY :dY);
1961 }
1962 
static AliTPCcalibDB * Instance()
Float_t fErrorRMSSys[2]
systematic relative error of the parametererization
Float_t GetPadPitchLength(Int_t isector=0, Int_t padrow=0) const
Definition: AliTPCParam.h:301
Float_t fErrorRMSQ[2][3][6]
shape parameterization coeficients
Float_t fParamS0[2][3][4]
error parameterization coeficients
Float_t fErrorRMS1[2][5]
shape parameterization coeficients
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TMatrixD * fQNormCorr
q norm correction for analytica correction
const Double_t kEpsilon
void FitRMS0(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
Float_t GetDiffL() const
Definition: AliTPCParam.h:337
void FitData(TTree *tree)
Float_t fErrorSQ[2][3][6]
error parameterization coeficients
Float_t GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const
Float_t GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const
Double_t AliTPCClusterParamClusterResolutionY(Double_t x, Double_t s0, Double_t s1)
#define TObjArray
static Double_t 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 effLength=0, Float_t effDiff=1)
void FitRMS1(TTree *tree, Int_t dim, Float_t *param0, Float_t *error)
TVectorD * fPosQTnorm[3]
q position normalization
Manager and of geomety classes for set: TPC.
Definition: AliTPCParam.h:18
static Double_t 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 effLength=0, Float_t effDiff=1)
TFile f("CalibObjects.root")
THnBase * fResolutionYMap
Map of resolution in Y.
static AliTPCClusterParam * Instance()
Float_t fParamRMS1[2][5]
shape parameterization coeficients
void FitResol0(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
Float_t GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const
void FitRMSSigma(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
void SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val, Int_t mode=1)
AliTPCParam * GetParameters() const
TVectorD * fPosZcor[3]
position correction parameterization
Float_t GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const
Float_t GetPadPitchWidth(Int_t isector=0) const
Definition: AliTPCParam.h:299
static Double_t GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1)
TTree * tree
Float_t fErrorS0Par[2][3][7]
error parameterization coeficients
Float_t fRatio
ratio of values constibution to error
Float_t fRMSSigmaFit[2][3][2]
mean value of the varation of RMS to RMS
ClassImp(TPCGenInfo)
Definition: AliTPCCmpNG.C:254
void FitRMS(TTree *tree)
Float_t fErrorS1[2][4]
error parameterization coeficients
static Float_t SGetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle)
void sum()
Float_t GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) const
void dY()
Definition: CalibAlign.C:547
Bool_t fWaveCorrectionMirroredPad
flag is the cog axis mirrored at 0.5
Float_t GetWaveCorrection(Int_t Type, Float_t Z, Int_t QMax, Float_t Pad, Float_t angleY) const
Float_t QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz)
TObjArray * fQNormHis
q norm correction for analytical correction
Float_t fParamS0Par[2][3][7]
error parameterization coeficients
TPC cluster error, shape and charge parameterization as function of drift length and inclination angl...
Double_t chi2
Definition: AnalyzeLaser.C:7
void FitResolQPar(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
Float_t GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const
Float_t GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const
TVectorD * fPosQMnorm[3]
q position normalization
Double_t GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const
void SetResolutionYMap(THnBase *ResolutionYMap)
void SetInstance(AliTPCClusterParam *const param)
Float_t GetZWidth() const
Definition: AliTPCParam.h:366
Float_t fParamRMSQ[2][3][6]
shape parameterization coeficients
Float_t Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz)
Float_t fParamS1[2][4]
error parameterization coeficients
Float_t 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)
AliTPCClusterParam & operator=(const AliTPCClusterParam &param)
void FitRMSQ(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
Float_t GetDiffT() const
Definition: AliTPCParam.h:336
Bool_t fWaveCorrectionMirroredAngle
flag is the Angle axis mirrored at 0
Float_t GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const
Float_t fParamSQPar[2][3][9]
error parameterization coeficients
Float_t fErrorSQPar[2][3][9]
error parameterization coeficients
void FitResol1(TTree *tree, Int_t dim, Float_t *param0, Float_t *error)
void FitResol(TTree *tree)
Bool_t fWaveCorrectionMirroredZ
flag is the Z axis mirrored at 0
virtual void Print(Option_t *option="") const
static AliTPCClusterParam * fgInstance
Double_t AliTPCClusterParamClusterResolutionZ(Double_t x, Double_t s0, Double_t s1)
Float_t GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM) const
Float_t GetWWPitch(Int_t isector=0) const
Definition: AliTPCParam.h:282
static Double_t GaussConvolutionGamma4(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau)
Float_t fRMSSigmaRatio[2][2]
mean value of the varation of RMS to RMS
void FitResolQ(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
void Test(TTree *tree, const char *output="TestClusterParam.root")
void FitResol0Par(TTree *tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error)
static Double_t GaussConvolutionTail(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau)
TObjArray * fQNorm
q norm paramters
THnBase * fWaveCorrectionMap
dY with respect to the distance to the center of the pad
Float_t fParamSQ[2][3][6]
error parameterization coeficients
Float_t GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle) const
Float_t fErrorRMS0[2][3][4]
shape parameterization coeficients
void SetWaveCorrectionMap(THnBase *WaveCorrectionMap)
Float_t fErrorS0[2][3][4]
error parameterization coeficients
TVectorD * fPosYcor[3]
position correction parameterization
Float_t GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const
Float_t 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)
Float_t fParamRMS0[2][3][4]
shape parameterization coeficients
void SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm)