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