AliRoot Core  edcc906 (edcc906)
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 
101 ClassImp(AliTPCClusterParam)
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
186  fWaveCorrectionMirroredPad( kFALSE ),
187  fWaveCorrectionMirroredZ( 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
211  fWaveCorrectionMirroredPad( kFALSE ),
212  fWaveCorrectionMirroredZ( 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  if (value<0.01) value = 0.01;
938  return value;
939 }
940 
941 
942 
943 Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
945 
946  Float_t value=0;
947  Float_t length=0.75;
948  if (type==1) length=1;
949  if (type==2) length=1.5;
950  value += fParamS1[dim][0];
951  value += fParamS1[dim][1]*z/length;
952  value += fParamS1[dim][2]*angle*angle*length;
953  value = TMath::Sqrt(TMath::Abs(value));
954  return value;
955 }
956 
957 Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
959 
960  Float_t value=0;
961  value += fParamSQ[dim][type][0];
962  value += fParamSQ[dim][type][1]*z;
963  value += fParamSQ[dim][type][2]*angle*angle;
964  value += fParamSQ[dim][type][3]*z/Qmean;
965  value += fParamSQ[dim][type][4]*angle*angle/Qmean;
966  value = TMath::Sqrt(TMath::Abs(value));
967  return value;
968 
969 
970 }
971 
972 Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
974 
975  Float_t value=0;
976  value += fParamSQPar[dim][type][0];
977  value += fParamSQPar[dim][type][1]*z;
978  value += fParamSQPar[dim][type][2]*angle*angle;
979  value += fParamSQPar[dim][type][3]*z*z;
980  value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
981  value += fParamSQPar[dim][type][5]*z*angle*angle;
982  value += fParamSQPar[dim][type][6]*z/Qmean;
983  value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
984  value = TMath::Sqrt(TMath::Abs(value));
985  return value;
986 
987 
988 }
989 
990 Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
992 
993  Float_t value=0;
994  value += fParamSQPar[dim][type][0];
995  value += fParamSQPar[dim][type][1]*z;
996  value += fParamSQPar[dim][type][2]*angle*angle;
997  value += fParamSQPar[dim][type][3]*z*z;
998  value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
999  value += fParamSQPar[dim][type][5]*z*angle*angle;
1000  value += fParamSQPar[dim][type][6]*z/Qmean;
1001  value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
1002  Float_t valueMean = GetError0Par(dim,type,z,angle);
1003  value -= 0.35*0.35*valueMean*valueMean;
1004  value = TMath::Sqrt(TMath::Abs(value));
1005  return value;
1006 
1007 
1008 }
1009 
1010 Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
1012 
1013  Float_t value=0;
1014  value += fParamRMS0[dim][type][0];
1015  value += fParamRMS0[dim][type][1]*z;
1016  value += fParamRMS0[dim][type][2]*angle*angle;
1017  value = TMath::Sqrt(TMath::Abs(value));
1018  return value;
1019 }
1020 
1021 Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const {
1023 
1024  Float_t value=0;
1025  Float_t length=0.75;
1026  if (type==1) length=1;
1027  if (type==2) length=1.5;
1028  if (type==0){
1029  value += fParamRMS1[dim][0];
1030  }else{
1031  value += fParamRMS1[dim][1];
1032  }
1033  value += fParamRMS1[dim][2]*z;
1034  value += fParamRMS1[dim][3]*angle*angle*length*length;
1035  value = TMath::Sqrt(TMath::Abs(value));
1036  return value;
1037 }
1038 
1039 Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
1041 
1042  Float_t value=0;
1043  value += fParamRMSQ[dim][type][0];
1044  value += fParamRMSQ[dim][type][1]*z;
1045  value += fParamRMSQ[dim][type][2]*angle*angle;
1046  value += fParamRMSQ[dim][type][3]*z/Qmean;
1047  value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
1048  value = TMath::Sqrt(TMath::Abs(value));
1049  return value;
1050 }
1051 
1052 Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const {
1054 
1055  Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
1056  Float_t value = fRMSSigmaFit[dim][type][0];
1057  value+= fRMSSigmaFit[dim][type][1]*mean;
1058  return value;
1059 }
1060 
1061 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 {
1063 
1064  Double_t value =0;
1065  //
1066  Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
1067  if (rmsL<rmsMeanQ) return value;
1068  //
1069  Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
1070  //
1071  if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
1072  //1.5 sigma cut on mean
1073  value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
1074  }else{
1075  if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
1076  //3 sigma cut on local
1077  value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
1078  }
1079  }
1080  return TMath::Sqrt(TMath::Abs(value));
1081 }
1082 
1083 
1084 
1087 
1088  FitResol(tree);
1089  FitRMS(tree);
1090 
1091 }
1092 
1095 
1096  SetInstance(this);
1097  for (Int_t idir=0;idir<2; idir++){
1098  for (Int_t itype=0; itype<3; itype++){
1099  Float_t param0[10];
1100  Float_t error0[10];
1101  // model error param
1102  FitResol0(tree, idir, itype,param0,error0);
1103  printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1104  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1105  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1106  for (Int_t ipar=0;ipar<4; ipar++){
1107  fParamS0[idir][itype][ipar] = param0[ipar];
1108  fErrorS0[idir][itype][ipar] = param0[ipar];
1109  }
1110  // error param with parabolic correction
1111  FitResol0Par(tree, idir, itype,param0,error0);
1112  printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
1113  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]);
1114  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]);
1115  for (Int_t ipar=0;ipar<7; ipar++){
1116  fParamS0Par[idir][itype][ipar] = param0[ipar];
1117  fErrorS0Par[idir][itype][ipar] = param0[ipar];
1118  }
1119  //
1120  FitResolQ(tree, idir, itype,param0,error0);
1121  printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1122  printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1123  printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1124  for (Int_t ipar=0;ipar<6; ipar++){
1125  fParamSQ[idir][itype][ipar] = param0[ipar];
1126  fErrorSQ[idir][itype][ipar] = param0[ipar];
1127  }
1128  //
1129  FitResolQPar(tree, idir, itype,param0,error0);
1130  printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
1131  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]);
1132  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]);
1133  for (Int_t ipar=0;ipar<9; ipar++){
1134  fParamSQPar[idir][itype][ipar] = param0[ipar];
1135  fErrorSQPar[idir][itype][ipar] = param0[ipar];
1136  }
1137  }
1138  }
1139  //
1140  printf("Resol z-scaled\n");
1141  for (Int_t idir=0;idir<2; idir++){
1142  Float_t param0[4];
1143  Float_t error0[4];
1144  FitResol1(tree, idir,param0,error0);
1145  printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
1146  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1147  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1148  for (Int_t ipar=0;ipar<4; ipar++){
1149  fParamS1[idir][ipar] = param0[ipar];
1150  fErrorS1[idir][ipar] = param0[ipar];
1151  }
1152  }
1153 
1154  for (Int_t idir=0;idir<2; idir++){
1155  printf("\nDirection %d\n",idir);
1156  printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1157  for (Int_t itype=0; itype<3; itype++){
1158  Float_t length=0.75;
1159  if (itype==1) length=1;
1160  if (itype==2) length=1.5;
1161  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));
1162  }
1163  }
1164 }
1165 
1166 
1167 
1170 
1171  SetInstance(this);
1172  for (Int_t idir=0;idir<2; idir++){
1173  for (Int_t itype=0; itype<3; itype++){
1174  Float_t param0[6];
1175  Float_t error0[6];
1176  FitRMS0(tree, idir, itype,param0,error0);
1177  printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1178  printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1179  printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1180  for (Int_t ipar=0;ipar<4; ipar++){
1181  fParamRMS0[idir][itype][ipar] = param0[ipar];
1182  fErrorRMS0[idir][itype][ipar] = param0[ipar];
1183  }
1184  FitRMSQ(tree, idir, itype,param0,error0);
1185  printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1186  printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1187  printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1188  for (Int_t ipar=0;ipar<6; ipar++){
1189  fParamRMSQ[idir][itype][ipar] = param0[ipar];
1190  fErrorRMSQ[idir][itype][ipar] = param0[ipar];
1191  }
1192  }
1193  }
1194  //
1195  printf("RMS z-scaled\n");
1196  for (Int_t idir=0;idir<2; idir++){
1197  Float_t param0[5];
1198  Float_t error0[5];
1199  FitRMS1(tree, idir,param0,error0);
1200  printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1201  printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1202  printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1203  for (Int_t ipar=0;ipar<5; ipar++){
1204  fParamRMS1[idir][ipar] = param0[ipar];
1205  fErrorRMS1[idir][ipar] = param0[ipar];
1206  }
1207  }
1208 
1209  for (Int_t idir=0;idir<2; idir++){
1210  printf("\nDirection %d\n",idir);
1211  printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1212  for (Int_t itype=0; itype<3; itype++){
1213  Float_t length=0.75;
1214  if (itype==1) length=1;
1215  if (itype==2) length=1.5;
1216  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);
1217  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);
1218  }
1219  }
1220  //
1221  // Fit RMS sigma
1222  //
1223  printf("RMS fluctuation parameterization \n");
1224  for (Int_t idir=0;idir<2; idir++){
1225  for (Int_t itype=0; itype<3; itype++){
1226  Float_t param0[5];
1227  Float_t error0[5];
1228  FitRMSSigma(tree, idir,itype,param0,error0);
1229  printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1230  for (Int_t ipar=0;ipar<2; ipar++){
1231  fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1232  }
1233  }
1234  }
1235  //
1236  // store systematic error end RMS fluctuation parameterization
1237  //
1238  TH1F hratio("hratio","hratio",100,-0.1,0.1);
1239  tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1240  fErrorRMSSys[0] = hratio.GetRMS();
1241  tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1242  fErrorRMSSys[1] = hratio.GetRMS();
1243  TH1F hratioR("hratioR","hratioR",100,0,0.2);
1244  tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1245  fRMSSigmaRatio[0][0]=hratioR.GetMean();
1246  fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1247  tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1248  fRMSSigmaRatio[1][0]=hratioR.GetMean();
1249  fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1250 }
1251 
1252 void AliTPCClusterParam::Test(TTree * tree, const char *output){
1254 
1255  SetInstance(this);
1256  TFile f(output,"recreate");
1257  f.cd();
1258  //
1259  // 1D histograms - resolution
1260  //
1261  for (Int_t idim=0; idim<2; idim++){
1262  for (Int_t ipad=0; ipad<3; ipad++){
1263  char hname1[300];
1264  char hcut1[300];
1265  char hexp1[300];
1266  //
1267  snprintf(hname1,300,"Delta0 Dir %d Pad %d",idim,ipad);
1268  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1269  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1270  TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1271  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1272  tree->Draw(hexp1,hcut1,"");
1273  his1DRel0.Write();
1274  //
1275  snprintf(hname1,300,"Delta0Par Dir %d Pad %d",idim,ipad);
1276  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1277  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1278  TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1279  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1280  tree->Draw(hexp1,hcut1,"");
1281  his1DRel0Par.Write();
1282  //
1283  }
1284  }
1285  //
1286  // 2D histograms - resolution
1287  //
1288  for (Int_t idim=0; idim<2; idim++){
1289  for (Int_t ipad=0; ipad<3; ipad++){
1290  char hname1[300];
1291  char hcut1[300];
1292  char hexp1[300];
1293  //
1294  snprintf(hname1,300,"2DDelta0 Dir %d Pad %d",idim,ipad);
1295  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1296  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1297  TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1298  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1299  tree->Draw(hexp1,hcut1,"");
1300  profDRel0.Write();
1301  //
1302  snprintf(hname1,300,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1303  snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1304  snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1305  TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1306  snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1307  tree->Draw(hexp1,hcut1,"");
1308  profDRel0Par.Write();
1309  //
1310  }
1311  }
1312 }
1313 
1314 
1315 
1316 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1318 
1319  //
1320  // Error parameterization
1321  //
1322  printf("\nResolution Scaled factors\n");
1323  printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1324  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])),
1325  TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
1326  for (Int_t ipad=0; ipad<3; ipad++){
1327  Float_t length=0.75;
1328  if (ipad==1) length=1;
1329  if (ipad==2) length=1.5;
1330  printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1331  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1332  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1333  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1334  TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
1335  }
1336  for (Int_t ipad=0; ipad<3; ipad++){
1337  Float_t length=0.75;
1338  if (ipad==1) length=1;
1339  if (ipad==2) length=1.5;
1340  printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1341  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1342  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1343  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1344  TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
1345  }
1346  printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1347  TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1348 
1349  for (Int_t ipad=0; ipad<3; ipad++){
1350  Float_t length=0.75;
1351  if (ipad==1) length=1;
1352  if (ipad==2) length=1.5;
1353  printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1354  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1355  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1356  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1357  TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
1358  }
1359  for (Int_t ipad=0; ipad<3; ipad++){
1360  Float_t length=0.75;
1361  if (ipad==1) length=1;
1362  if (ipad==2) length=1.5;
1363  printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1364  TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1365  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1366  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1367  TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
1368  }
1369 
1370  //
1371  // RMS scaling
1372  //
1373  printf("\n");
1374  printf("\nRMS Scaled factors\n");
1375  printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1376  printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1377  TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1378  TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1379  TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1380  TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1381  TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
1382  for (Int_t ipad=0; ipad<3; ipad++){
1383  Float_t length=0.75;
1384  if (ipad==1) length=1;
1385  if (ipad==2) length=1.5;
1386  if (ipad==0){
1387  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1388  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1389  0.,
1390  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1391  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1392  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1393  }else{
1394  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1395  0.,
1396  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1397  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1398  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1399  TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
1400  }
1401  }
1402  printf("\n");
1403  printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1404  TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1405  TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1406  TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1407  TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1408  TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
1409  for (Int_t ipad=0; ipad<3; ipad++){
1410  Float_t length=0.75;
1411  if (ipad==1) length=1;
1412  if (ipad==2) length=1.5;
1413  if (ipad==0){
1414  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1415  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1416  0.,
1417  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1418  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1419  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1420  }else{
1421  printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1422  0.,
1423  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1424  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1425  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1426  TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
1427  }
1428  }
1429  printf("\ndEdx correction matrix used in GetQnormCorr\n");
1430  fQNormCorr->Print();
1431 
1432 }
1433 
1434 
1435 
1436 
1437 
1438 Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1444 
1445  if (fQNorm==0) return 0;
1446  TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1447  if (!norm) return 0;
1448  TVectorD &no = *norm;
1449  Float_t res =
1450  no[0]+
1451  no[1]*dr+
1452  no[2]*ty+
1453  no[3]*tz+
1454  no[4]*dr*ty+
1455  no[5]*dr*tz+
1456  no[6]*ty*tz+
1457  no[7]*dr*dr+
1458  no[8]*ty*ty+
1459  no[9]*tz*tz;
1460  res/=no[0];
1461  return res;
1462 }
1463 
1464 
1465 
1466 Float_t AliTPCClusterParam::QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t p2, Float_t p3){
1470 
1471  if (fQNormHis==0) return 0;
1472  TH3F * norm = (TH3F*)fQNormHis->At(4*itype+ipad);
1473  if (!norm) return 1;
1474  p2=TMath::Abs(p2);
1475  dr=TMath::Min(dr,Float_t(norm->GetXaxis()->GetXmax()-norm->GetXaxis()->GetBinWidth(0)));
1476  dr=TMath::Max(dr,Float_t(norm->GetXaxis()->GetXmin()+norm->GetXaxis()->GetBinWidth(0)));
1477  //
1478  p2=TMath::Min(p2,Float_t(norm->GetYaxis()->GetXmax()-norm->GetYaxis()->GetBinWidth(0)));
1479  p2=TMath::Max(p2,Float_t(norm->GetYaxis()->GetXmin()+norm->GetYaxis()->GetBinWidth(0)));
1480  //
1481  p3=TMath::Min(p3,Float_t(norm->GetZaxis()->GetXmax()-norm->GetZaxis()->GetBinWidth(0)));
1482  p3=TMath::Max(p3,Float_t(norm->GetZaxis()->GetXmin()+norm->GetZaxis()->GetBinWidth(0)));
1483  //
1484  Double_t res = norm->GetBinContent(norm->FindBin(dr,p2,p3));
1485  if (res==0) res = norm->GetBinContent(norm->FindBin(0.5,0.5,0.5)); // This is just hack - to be fixed entries without
1486 
1487  return res;
1488 }
1489 
1490 
1491 
1492 void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm){
1497 
1498  if (fQNorm==0) fQNorm = new TObjArray(6);
1499  fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1500 }
1501 
1504 
1505  if (!fQNormCorr) fQNormCorr= new TMatrixD(12,6);
1506  for (Int_t irow=0;irow<12; irow++)
1507  for (Int_t icol=0;icol<6; icol++){
1508  (*fQNormCorr)(irow,icol)=1.; // default - no correction
1509  if (irow>5) (*fQNormCorr)(irow,icol)=0.; // default - no correction
1510  }
1511 }
1512 
1513 void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val, Int_t mode){
1522 
1523  if (!fQNormCorr) {
1524  ResetQnormCorr();
1525  }
1526  //
1527  // eff shap parameterization matrix
1528  //
1529  // rows
1530  // itype*3+ipad - itype=0 qtot itype=1 qmax ipad=0
1531  //
1532  if (mode==1) {
1533  // mode introduced in 20.07.2014 - git describe ~ vAN-20140703-48-g3449a97 - to keep back compatibility with o
1534  (*fQNormCorr)(itype*3+ipad, corrType) = val; // set value
1535  return;
1536  }
1537  if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val; // multiplicative correction
1538  if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val; // additive correction
1539 }
1540 
1541 Double_t AliTPCClusterParam::GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const{
1543 
1544  if (!fQNormCorr) return 0;
1545  return (*fQNormCorr)(itype*3+ipad, corrType);
1546 }
1547 
1548 
1549 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){
1559 
1560  //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
1561  // Following variable used - correspondance to the our variable conventions
1562  //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1563  Double_t dp = ((pad-int(pad)-0.5)*2.);
1564  //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1565  Double_t dt = ((time-int(time)-0.5)*2.);
1566  //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1567  Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1568  //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1569  Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1570  //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1571  Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1572  //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1573  Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1574  //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1575  Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1576  //
1577  //
1578  //
1579  TVectorD * pvec = 0;
1580  if (isMax){
1581  pvec = fPosQMnorm[ipad];
1582  }else{
1583  pvec = fPosQTnorm[ipad];
1584  }
1585  TVectorD &param = *pvec;
1586  //
1587  // Eval part - in correspondance with fit part from debug streamer
1588  //
1589  Double_t result=param[0];
1590  Int_t index =1;
1591  //
1592  result+=dp*param[index++]; //1
1593  result+=dt*param[index++]; //2
1594  result+=dp*dp*param[index++]; //3
1595  result+=dt*dt*param[index++]; //4
1596  result+=dt*dt*dt*param[index++]; //5
1597  result+=dp*dt*param[index++]; //6
1598  result+=dp*dt*dt*param[index++]; //7
1599  result+=(dq0)*param[index++]; //8
1600  result+=(dq1)*param[index++]; //9
1601  //
1602  //
1603  result+=dp*dp*(di)*param[index++]; //10
1604  result+=dt*dt*(di)*param[index++]; //11
1605  result+=dp*dp*sy*param[index++]; //12
1606  result+=dt*sz*param[index++]; //13
1607  result+=dt*dt*sz*param[index++]; //14
1608  result+=dt*dt*dt*sz*param[index++]; //15
1609  //
1610  result+=dp*dp*1*sy*sz*param[index++]; //16
1611  result+=dt*sy*sz*param[index++]; //17
1612  result+=dt*dt*sy*sz*param[index++]; //18
1613  result+=dt*dt*dt*sy*sz*param[index++]; //19
1614  //
1615  result+=dp*dp*(dq0)*param[index++]; //20
1616  result+=dt*1*(dq0)*param[index++]; //21
1617  result+=dt*dt*(dq0)*param[index++]; //22
1618  result+=dt*dt*dt*(dq0)*param[index++]; //23
1619  //
1620  result+=dp*dp*(dq1)*param[index++]; //24
1621  result+=dt*(dq1)*param[index++]; //25
1622  result+=dt*dt*(dq1)*param[index++]; //26
1623  result+=dt*dt*dt*(dq1)*param[index++]; //27
1624 
1625  if (result<0.75) result=0.75;
1626  if (result>1.25) result=1.25;
1627 
1628  return result;
1629 
1630 }
1631 
1632 
1633 
1634 
1635 
1636 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*/){
1637 
1645 
1646  //
1647  //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
1648  //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
1649  //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
1650  //chainres->SetAlias("st","(sin(dt)-dt)");
1651  //
1652  //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1653 
1654  //
1655  // Derived variables
1656  //
1657  Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
1658  Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
1659  Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
1660  Double_t st = (TMath::Sin(dt)-dt);
1661  //
1662  Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.)));
1663  //
1664  //
1665  //
1666  TVectorD * pvec = 0;
1667  if (type==0){
1668  pvec = fPosYcor[ipad];
1669  }else{
1670  pvec = fPosZcor[ipad];
1671  }
1672  TVectorD &param = *pvec;
1673  //
1674  Double_t result=0;
1675  Int_t index =1;
1676 
1677  if (type==0){
1678  // y corr
1679  result+=(dp)*param[index++]; //1
1680  result+=(dp)*di*param[index++]; //2
1681  //
1682  result+=(sp)*param[index++]; //3
1683  result+=(sp)*di*param[index++]; //4
1684  }
1685  if (type==1){
1686  result+=(dt)*param[index++]; //1
1687  result+=(dt)*di*param[index++]; //2
1688  //
1689  result+=(st)*param[index++]; //3
1690  result+=(st)*di*param[index++]; //4
1691  }
1692  if (TMath::Abs(result)>0.05) return 0;
1693  return result;
1694 }
1695 
1696 
1697 
1698 Double_t AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
1705 
1706  const Double_t kEpsilon = 0.0001;
1707  const Double_t twoPi = TMath::TwoPi();
1708  const Double_t hnorm = 0.5/TMath::Sqrt(twoPi);
1709  const Double_t sqtwo = TMath::Sqrt(2.);
1710 
1711  if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
1712  // small angular effect
1713  Double_t val = TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1)/(s0*s1*twoPi);
1714  return val;
1715  }
1716  Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
1717  Double_t sigma = TMath::Sqrt(sigma2);
1718  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2.*sigma2));
1719  //
1720  Double_t sigmaErf = 1./(2.*s0*s1*sqtwo*sigma);
1721  Double_t k0s1s1 = 2.*k0*s1*s1;
1722  Double_t k1s0s0 = 2.*k1*s0*s0;
1723  Double_t erf0 = AliMathBase::ErfFast((sigma2-k0s1s1*x0-k1s0s0*x1)*sigmaErf);
1724  Double_t erf1 = AliMathBase::ErfFast((sigma2+k0s1s1*x0+k1s0s0*x1)*sigmaErf);
1725  Double_t norm = hnorm/sigma;
1726  Double_t val = norm*exp0*(erf0+erf1);
1727  return val;
1728 }
1729 
1730 
1731 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){
1738 
1739  Double_t sum =1, mean=0;
1740  // the COG of exponent
1741  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1742  mean+=iexp*TMath::Exp(-iexp/tau);
1743  sum +=TMath::Exp(-iexp/tau);
1744  }
1745  mean/=sum;
1746  //
1747  sum = 1;
1748  Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
1749  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1750  val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
1751  sum+=TMath::Exp(-iexp/tau);
1752  }
1753  return val/sum;
1754 }
1755 
1756 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){
1764 
1765  Double_t sum =0, mean=0;
1766  // the COG of G4
1767  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1768  Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1769  mean+=iexp*g4;
1770  sum +=g4;
1771  }
1772  mean/=sum;
1773  //
1774  sum = 0;
1775  Double_t val = 0;
1776  for (Float_t iexp=0;iexp<5;iexp+=0.2){
1777  Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
1778  val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
1779  sum+=g4;
1780  }
1781  return val/sum;
1782 }
1783 
1784 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){
1793 
1794  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1795  //
1796  // Gaus width sy and sz is determined by RF width and diffusion
1797  // Integral of Q is equal 1
1798  // Q max is calculated at position cpad, ctime
1799  // Example function:
1800  // TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000)
1801  //
1803  Double_t padLength= param->GetPadPitchLength(sector,row);
1804  Double_t padWidth = param->GetPadPitchWidth(sector);
1805  Double_t zwidth = param->GetZWidth();
1806  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1807 
1808  // diffusion in pad, time bin units
1809  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1810  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1811  diffT*=effDiff; //
1812  diffL*=effDiff; //
1813  //
1814  // transform angular effect to pad units
1815  //
1816  Double_t pky = ky*effLength/padWidth;
1817  Double_t pkz = kz*effLength/zwidth;
1818  // position in pad unit
1819  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1820  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1821  //
1822  //
1823  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1824  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1825  //return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
1826  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1827  return GaussConvolution(py,pz, pky,pkz,sy,sz)*length;
1828 }
1829 
1830 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){
1841 
1842  // Response function aproximated by convolution of gaussian with angular effect (integral=1)
1843  //
1844  // Gaus width sy and sz is determined by RF width and diffusion
1845  // Integral of Q is equal 1
1846  // Q max is calculated at position cpad, ctime
1847  //
1848  //
1849  //
1851  Double_t padLength= param->GetPadPitchLength(sector,row);
1852  Double_t padWidth = param->GetPadPitchWidth(sector);
1853  Double_t zwidth = param->GetZWidth();
1854  Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
1855  //
1856  // diffusion in pad units
1857  Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
1858  Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
1859  diffT*=effDiff; //
1860  diffL*=effDiff; //
1861  //
1862  // transform angular effect to pad units
1863  Double_t pky = ky*effLength/padWidth;
1864  Double_t pkz = kz*effLength/zwidth;
1865  // position in pad unit
1866  //
1867  Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
1868  Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
1869  //
1870  Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
1871  Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
1872  //
1873  //
1874  //
1875  Double_t sumAll=0,sumThr=0;
1876  //
1877  Double_t corr =1;
1878  Double_t qnorm=qtot;
1879  for (Float_t iy=-3;iy<=3;iy+=1.)
1880  for (Float_t iz=-4;iz<=4;iz+=1.){
1881  // Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);
1882  Double_t val = GaussConvolution(py-iy,pz-iz, pky,pkz, sy,sz);
1883  Double_t qlocal =qnorm*val;
1884  if (TMath::Abs(iy)<1.5&&TMath::Abs(iz)<1.5){
1885  sumThr+=qlocal; // Virtual charge used in cluster finder
1886  }
1887  else{
1888  if (qlocal>thr && TMath::Abs(iz)<2.5&&TMath::Abs(iy)<2.5) sumThr+=qlocal;
1889  }
1890  sumAll+=qlocal;
1891  }
1892  if (sumAll>0&&sumThr>0) {
1893  corr=(sumThr)/sumAll;
1894  }
1895  //
1896  Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
1897  return corr*length;
1898 }
1899 
1900 
1901 
1903 {
1905 
1906  delete fWaveCorrectionMap;
1907  fWaveCorrectionMap = 0;
1908  fWaveCorrectionMirroredPad = kFALSE;
1909  fWaveCorrectionMirroredZ = kFALSE;
1911  if( Map ){
1912  fWaveCorrectionMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
1913  if( fWaveCorrectionMap ){
1914  fWaveCorrectionMirroredPad = ( fWaveCorrectionMap->GetAxis(3)->FindFixBin(0.5)<=1 ); // Pad axis is mirrored at 0.5
1915  fWaveCorrectionMirroredZ = ( fWaveCorrectionMap->GetAxis(1)->FindFixBin(0.0)<=1); // Z axis is mirrored at 0
1916  fWaveCorrectionMirroredAngle = ( fWaveCorrectionMap->GetAxis(4)->FindFixBin(0.0)<=1 ); // Angle axis is mirrored at 0
1917  }
1918  }
1919 }
1920 
1922 {
1924 
1925  delete fResolutionYMap;
1926  fResolutionYMap = 0;
1927  if( Map ){
1928  fResolutionYMap = dynamic_cast<THnBase*>( Map->Clone(Map->GetName()));
1929  }
1930 }
1931 
1932 Float_t AliTPCClusterParam::GetWaveCorrection(Int_t Type, Float_t Z, Int_t QMax, Float_t Pad, Float_t angleY ) const
1933 {
1935 
1936  if( !fWaveCorrectionMap ) return 0;
1937  Bool_t swapY = kFALSE;
1938  Pad = Pad-(Int_t)Pad;
1939 
1940  if( TMath::Abs(Pad-0.5)<1.e-8 ){// one pad clusters a stored in underflow bins
1941  Pad = -1.;
1942  } else {
1943  if( fWaveCorrectionMirroredPad && (Pad<0.5) ){ // cog axis is mirrored at 0.5
1944  swapY = !swapY;
1945  Pad = 1.0 - Pad;
1946  }
1947  }
1948 
1949  if( fWaveCorrectionMirroredZ && (Z<0) ){ // Z axis is mirrored at 0
1950  swapY = !swapY;
1951  Z = -Z;
1952  }
1953 
1954  if( fWaveCorrectionMirroredAngle && (angleY<0) ){ // Angle axis is mirrored at 0
1955  angleY = -angleY;
1956  }
1957  double var[5] = { static_cast<double>(Type), Z, static_cast<double>(QMax), Pad, angleY };
1958  Long64_t bin = fWaveCorrectionMap->GetBin(var, kFALSE );
1959  if( bin<0 ) return 0;
1960  Double_t dY = fWaveCorrectionMap->GetBinContent(bin);
1961  return (swapY ?-dY :dY);
1962 }
1963 
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:302
Float_t fErrorRMSQ[2][3][6]
shape parameterization coeficients
TDatime dt
Definition: pdc06_config.C:130
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:338
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)
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:300
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
void FitRMS(TTree *tree)
Float_t fErrorS1[2][4]
error parameterization coeficients
TVectorD * fQpadMnorm
q pad normalization - Max charge
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
TVectorD * fQpadTnorm
q pad normalization - Total charge
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:367
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)
static Double_t ErfFast(Double_t x)
Definition: AliMathBase.h:39
Float_t GetDiffT() const
Definition: AliTPCParam.h:337
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)
TF1 * f
Definition: interpolTest.C:21
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:283
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
class TVectorT< Double_t > TVectorD
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)
void res(Char_t i)
Definition: Resolution.C:2
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
class TMatrixT< Double_t > TMatrixD
void SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm)