1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // TPC cluster error and shape parameterization //
22 ///////////////////////////////////////////////////////////////////////////////
23 #include "AliTPCClusterParam.h"
28 #include <TLinearFitter.h>
30 #include <TProfile2D.h>
32 ClassImp(AliTPCClusterParam)
35 AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
39 Example usage fitting parameterization:
40 TFile fres("resol.root"); //tree with resolution and shape
41 TTree * treeRes =(TTree*)fres.Get("Resol");
43 AliTPCClusterParam param;
44 param.SetInstance(¶m);
45 param.FitResol(treeRes);
46 param.FitRMS(treeRes);
47 TFile fparam("TPCClusterParam.root","recreate");
51 TFile fparam("TPCClusterParam.root");
52 AliTPCClusterParam *param2 = (AliTPCClusterParam *) fparam.Get("Param");
53 param2->SetInstance(param2);
54 param2->Test(treeRes);
57 treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
64 //_ singleton implementation __________________________________________________
65 AliTPCClusterParam* AliTPCClusterParam::Instance()
68 // Singleton implementation
69 // Returns an instance of this class, it is created if neccessary
72 fgInstance = new AliTPCClusterParam();
80 void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
82 // Fit z - angular dependence of resolution
84 // Int_t dim=0, type=0;
86 sprintf(varVal,"Resol:AngleM:Zm");
88 sprintf(varErr,"Sigma:AngleS:Zs");
90 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
92 Int_t entries = tree->Draw(varVal,varCut);
93 Float_t px[10000], py[10000], pz[10000];
94 Float_t ex[10000], ey[10000], ez[10000];
96 tree->Draw(varErr,varCut);
97 for (Int_t ipoint=0; ipoint<entries; ipoint++){
98 ex[ipoint]= tree->GetV3()[ipoint];
99 ey[ipoint]= tree->GetV2()[ipoint];
100 ez[ipoint]= tree->GetV1()[ipoint];
102 tree->Draw(varVal,varCut);
103 for (Int_t ipoint=0; ipoint<entries; ipoint++){
104 px[ipoint]= tree->GetV3()[ipoint];
105 py[ipoint]= tree->GetV2()[ipoint];
106 pz[ipoint]= tree->GetV1()[ipoint];
110 TLinearFitter fitter(3,"hyp2");
111 for (Int_t ipoint=0; ipoint<entries; ipoint++){
112 Float_t val = pz[ipoint]*pz[ipoint];
113 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
116 x[1] = py[ipoint]*py[ipoint];
117 fitter.AddPoint(x,val,err);
121 fitter.GetParameters(param);
122 param0[0] = param[0];
123 param0[1] = param[1];
124 param0[2] = param[2];
125 Float_t chi2 = fitter.GetChisquare()/entries;
127 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
128 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
129 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
133 void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
135 // Fit z - angular dependence of resolution
137 // Int_t dim=0, type=0;
139 sprintf(varVal,"Resol:AngleM:Zm");
141 sprintf(varErr,"Sigma:AngleS:Zs");
143 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
145 Int_t entries = tree->Draw(varVal,varCut);
146 Float_t px[10000], py[10000], pz[10000];
147 Float_t ex[10000], ey[10000], ez[10000];
149 tree->Draw(varErr,varCut);
150 for (Int_t ipoint=0; ipoint<entries; ipoint++){
151 ex[ipoint]= tree->GetV3()[ipoint];
152 ey[ipoint]= tree->GetV2()[ipoint];
153 ez[ipoint]= tree->GetV1()[ipoint];
155 tree->Draw(varVal,varCut);
156 for (Int_t ipoint=0; ipoint<entries; ipoint++){
157 px[ipoint]= tree->GetV3()[ipoint];
158 py[ipoint]= tree->GetV2()[ipoint];
159 pz[ipoint]= tree->GetV1()[ipoint];
163 TLinearFitter fitter(6,"hyp5");
164 for (Int_t ipoint=0; ipoint<entries; ipoint++){
165 Float_t val = pz[ipoint]*pz[ipoint];
166 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
169 x[1] = py[ipoint]*py[ipoint];
173 fitter.AddPoint(x,val,err);
177 fitter.GetParameters(param);
178 param0[0] = param[0];
179 param0[1] = param[1];
180 param0[2] = param[2];
181 param0[3] = param[3];
182 param0[4] = param[4];
183 param0[5] = param[5];
184 Float_t chi2 = fitter.GetChisquare()/entries;
186 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
187 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
188 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
189 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
190 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
191 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
198 void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
200 // Fit z - angular dependence of resolution - pad length scaling
202 // Int_t dim=0, type=0;
204 sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
206 sprintf(varErr,"Sigma:AngleS:Zs");
208 sprintf(varCut,"Dim==%d&&QMean<0",dim);
210 Int_t entries = tree->Draw(varVal,varCut);
211 Float_t px[10000], py[10000], pz[10000];
212 Float_t ex[10000], ey[10000], ez[10000];
214 tree->Draw(varErr,varCut);
215 for (Int_t ipoint=0; ipoint<entries; ipoint++){
216 ex[ipoint]= tree->GetV3()[ipoint];
217 ey[ipoint]= tree->GetV2()[ipoint];
218 ez[ipoint]= tree->GetV1()[ipoint];
220 tree->Draw(varVal,varCut);
221 for (Int_t ipoint=0; ipoint<entries; ipoint++){
222 px[ipoint]= tree->GetV3()[ipoint];
223 py[ipoint]= tree->GetV2()[ipoint];
224 pz[ipoint]= tree->GetV1()[ipoint];
228 TLinearFitter fitter(3,"hyp2");
229 for (Int_t ipoint=0; ipoint<entries; ipoint++){
230 Float_t val = pz[ipoint]*pz[ipoint];
231 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
234 x[1] = py[ipoint]*py[ipoint];
235 fitter.AddPoint(x,val,err);
239 fitter.GetParameters(param);
240 param0[0] = param[0];
241 param0[1] = param[1];
242 param0[2] = param[2];
243 Float_t chi2 = fitter.GetChisquare()/entries;
245 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
246 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
247 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
250 void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
252 // Fit z - angular dependence of resolution - Q scaling
254 // Int_t dim=0, type=0;
256 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
258 sprintf(varVal0,"Resol:AngleM:Zm");
261 sprintf(varErr,"Sigma:AngleS:Zs");
263 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
265 Int_t entries = tree->Draw(varVal,varCut);
266 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
267 Float_t ex[20000], ey[20000], ez[20000];
269 tree->Draw(varErr,varCut);
270 for (Int_t ipoint=0; ipoint<entries; ipoint++){
271 ex[ipoint]= tree->GetV3()[ipoint];
272 ey[ipoint]= tree->GetV2()[ipoint];
273 ez[ipoint]= tree->GetV1()[ipoint];
275 tree->Draw(varVal,varCut);
276 for (Int_t ipoint=0; ipoint<entries; ipoint++){
277 px[ipoint]= tree->GetV3()[ipoint];
278 py[ipoint]= tree->GetV2()[ipoint];
279 pz[ipoint]= tree->GetV1()[ipoint];
281 tree->Draw(varVal0,varCut);
282 for (Int_t ipoint=0; ipoint<entries; ipoint++){
283 pu[ipoint]= tree->GetV3()[ipoint];
284 pt[ipoint]= tree->GetV2()[ipoint];
288 TLinearFitter fitter(5,"hyp4");
289 for (Int_t ipoint=0; ipoint<entries; ipoint++){
290 Float_t val = pz[ipoint]*pz[ipoint];
291 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
294 x[1] = pt[ipoint]*pt[ipoint];
296 x[3] = py[ipoint]*py[ipoint];
297 fitter.AddPoint(x,val,err);
302 fitter.GetParameters(param);
303 param0[0] = param[0];
304 param0[1] = param[1];
305 param0[2] = param[2];
306 param0[3] = param[3];
307 param0[4] = param[4];
308 Float_t chi2 = fitter.GetChisquare()/entries;
310 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
311 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
312 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
313 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
314 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
317 void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
319 // Fit z - angular dependence of resolution - Q scaling - parabolic correction
321 // Int_t dim=0, type=0;
323 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
325 sprintf(varVal0,"Resol:AngleM:Zm");
328 sprintf(varErr,"Sigma:AngleS:Zs");
330 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
332 Int_t entries = tree->Draw(varVal,varCut);
333 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
334 Float_t ex[20000], ey[20000], ez[20000];
336 tree->Draw(varErr,varCut);
337 for (Int_t ipoint=0; ipoint<entries; ipoint++){
338 ex[ipoint]= tree->GetV3()[ipoint];
339 ey[ipoint]= tree->GetV2()[ipoint];
340 ez[ipoint]= tree->GetV1()[ipoint];
342 tree->Draw(varVal,varCut);
343 for (Int_t ipoint=0; ipoint<entries; ipoint++){
344 px[ipoint]= tree->GetV3()[ipoint];
345 py[ipoint]= tree->GetV2()[ipoint];
346 pz[ipoint]= tree->GetV1()[ipoint];
348 tree->Draw(varVal0,varCut);
349 for (Int_t ipoint=0; ipoint<entries; ipoint++){
350 pu[ipoint]= tree->GetV3()[ipoint];
351 pt[ipoint]= tree->GetV2()[ipoint];
355 TLinearFitter fitter(8,"hyp7");
356 for (Int_t ipoint=0; ipoint<entries; ipoint++){
357 Float_t val = pz[ipoint]*pz[ipoint];
358 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
361 x[1] = pt[ipoint]*pt[ipoint];
366 x[6] = py[ipoint]*py[ipoint];
368 fitter.AddPoint(x,val,err);
373 fitter.GetParameters(param);
374 param0[0] = param[0];
375 param0[1] = param[1];
376 param0[2] = param[2];
377 param0[3] = param[3];
378 param0[4] = param[4];
379 param0[5] = param[5];
380 param0[6] = param[6];
381 param0[7] = param[7];
383 Float_t chi2 = fitter.GetChisquare()/entries;
385 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
386 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
387 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
388 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
389 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
390 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
391 error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
392 error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
397 void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
399 // Fit z - angular dependence of resolution
401 // Int_t dim=0, type=0;
403 sprintf(varVal,"RMSm:AngleM:Zm");
405 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
407 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
409 Int_t entries = tree->Draw(varVal,varCut);
410 Float_t px[10000], py[10000], pz[10000];
411 Float_t ex[10000], ey[10000], ez[10000];
413 tree->Draw(varErr,varCut);
414 for (Int_t ipoint=0; ipoint<entries; ipoint++){
415 ex[ipoint]= tree->GetV3()[ipoint];
416 ey[ipoint]= tree->GetV2()[ipoint];
417 ez[ipoint]= tree->GetV1()[ipoint];
419 tree->Draw(varVal,varCut);
420 for (Int_t ipoint=0; ipoint<entries; ipoint++){
421 px[ipoint]= tree->GetV3()[ipoint];
422 py[ipoint]= tree->GetV2()[ipoint];
423 pz[ipoint]= tree->GetV1()[ipoint];
427 TLinearFitter fitter(3,"hyp2");
428 for (Int_t ipoint=0; ipoint<entries; ipoint++){
429 Float_t val = pz[ipoint]*pz[ipoint];
430 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
433 x[1] = py[ipoint]*py[ipoint];
434 fitter.AddPoint(x,val,err);
438 fitter.GetParameters(param);
439 param0[0] = param[0];
440 param0[1] = param[1];
441 param0[2] = param[2];
442 Float_t chi2 = fitter.GetChisquare()/entries;
444 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
445 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
446 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
449 void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
451 // Fit z - angular dependence of resolution - pad length scaling
453 // Int_t dim=0, type=0;
455 sprintf(varVal,"RMSm:AngleM*Length:Zm");
457 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
459 sprintf(varCut,"Dim==%d&&QMean<0",dim);
461 Int_t entries = tree->Draw(varVal,varCut);
462 Float_t px[10000], py[10000], pz[10000];
463 Float_t type[10000], ey[10000], ez[10000];
465 tree->Draw(varErr,varCut);
466 for (Int_t ipoint=0; ipoint<entries; ipoint++){
467 type[ipoint] = tree->GetV3()[ipoint];
468 ey[ipoint] = tree->GetV2()[ipoint];
469 ez[ipoint] = tree->GetV1()[ipoint];
471 tree->Draw(varVal,varCut);
472 for (Int_t ipoint=0; ipoint<entries; ipoint++){
473 px[ipoint]= tree->GetV3()[ipoint];
474 py[ipoint]= tree->GetV2()[ipoint];
475 pz[ipoint]= tree->GetV1()[ipoint];
479 TLinearFitter fitter(4,"hyp3");
480 for (Int_t ipoint=0; ipoint<entries; ipoint++){
481 Float_t val = pz[ipoint]*pz[ipoint];
482 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
484 x[0] = (type[ipoint]<0.5)? 0.:1.;
486 x[2] = py[ipoint]*py[ipoint];
487 fitter.AddPoint(x,val,err);
491 fitter.GetParameters(param);
492 param0[0] = param[0];
493 param0[1] = param[0]+param[1];
494 param0[2] = param[2];
495 param0[3] = param[3];
496 Float_t chi2 = fitter.GetChisquare()/entries;
498 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
499 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
500 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
501 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
504 void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
506 // Fit z - angular dependence of resolution - Q scaling
508 // Int_t dim=0, type=0;
510 sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
512 sprintf(varVal0,"RMSm:AngleM:Zm");
515 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
517 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
519 Int_t entries = tree->Draw(varVal,varCut);
520 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
521 Float_t ex[20000], ey[20000], ez[20000];
523 tree->Draw(varErr,varCut);
524 for (Int_t ipoint=0; ipoint<entries; ipoint++){
525 ex[ipoint]= tree->GetV3()[ipoint];
526 ey[ipoint]= tree->GetV2()[ipoint];
527 ez[ipoint]= tree->GetV1()[ipoint];
529 tree->Draw(varVal,varCut);
530 for (Int_t ipoint=0; ipoint<entries; ipoint++){
531 px[ipoint]= tree->GetV3()[ipoint];
532 py[ipoint]= tree->GetV2()[ipoint];
533 pz[ipoint]= tree->GetV1()[ipoint];
535 tree->Draw(varVal0,varCut);
536 for (Int_t ipoint=0; ipoint<entries; ipoint++){
537 pu[ipoint]= tree->GetV3()[ipoint];
538 pt[ipoint]= tree->GetV2()[ipoint];
542 TLinearFitter fitter(5,"hyp4");
543 for (Int_t ipoint=0; ipoint<entries; ipoint++){
544 Float_t val = pz[ipoint]*pz[ipoint];
545 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
548 x[1] = pt[ipoint]*pt[ipoint];
550 x[3] = py[ipoint]*py[ipoint];
551 fitter.AddPoint(x,val,err);
556 fitter.GetParameters(param);
557 param0[0] = param[0];
558 param0[1] = param[1];
559 param0[2] = param[2];
560 param0[3] = param[3];
561 param0[4] = param[4];
562 Float_t chi2 = fitter.GetChisquare()/entries;
564 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
565 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
566 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
567 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
568 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
572 void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
574 // Fit z - angular dependence of resolution - Q scaling
576 // Int_t dim=0, type=0;
578 sprintf(varVal,"RMSs:RMSm");
581 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
583 Int_t entries = tree->Draw(varVal,varCut);
584 Float_t px[20000], py[20000];
586 tree->Draw(varVal,varCut);
587 for (Int_t ipoint=0; ipoint<entries; ipoint++){
588 px[ipoint]= tree->GetV2()[ipoint];
589 py[ipoint]= tree->GetV1()[ipoint];
591 TLinearFitter fitter(2,"pol1");
592 for (Int_t ipoint=0; ipoint<entries; ipoint++){
593 Float_t val = py[ipoint];
594 Float_t err = fRatio*px[ipoint];
597 fitter.AddPoint(x,val,err);
600 param0[0]= fitter.GetParameter(0);
601 param0[1]= fitter.GetParameter(1);
606 Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
611 value += fParamS0[dim][type][0];
612 value += fParamS0[dim][type][1]*z;
613 value += fParamS0[dim][type][2]*angle*angle;
614 value = TMath::Sqrt(TMath::Abs(value));
619 Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
624 value += fParamS0Par[dim][type][0];
625 value += fParamS0Par[dim][type][1]*z;
626 value += fParamS0Par[dim][type][2]*angle*angle;
627 value += fParamS0Par[dim][type][3]*z*z;
628 value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
629 value += fParamS0Par[dim][type][5]*z*angle*angle;
630 value = TMath::Sqrt(TMath::Abs(value));
636 Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
642 if (type==1) length=1;
643 if (type==2) length=1.5;
644 value += fParamS1[dim][0];
645 value += fParamS1[dim][1]*z/length;
646 value += fParamS1[dim][2]*angle*angle*length;
647 value = TMath::Sqrt(TMath::Abs(value));
651 Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
656 value += fParamSQ[dim][type][0];
657 value += fParamSQ[dim][type][1]*z;
658 value += fParamSQ[dim][type][2]*angle*angle;
659 value += fParamSQ[dim][type][3]*z/Qmean;
660 value += fParamSQ[dim][type][4]*angle*angle/Qmean;
661 value = TMath::Sqrt(TMath::Abs(value));
667 Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
672 value += fParamSQPar[dim][type][0];
673 value += fParamSQPar[dim][type][1]*z;
674 value += fParamSQPar[dim][type][2]*angle*angle;
675 value += fParamSQPar[dim][type][3]*z*z;
676 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
677 value += fParamSQPar[dim][type][5]*z*angle*angle;
678 value += fParamSQPar[dim][type][6]*z/Qmean;
679 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
680 value = TMath::Sqrt(TMath::Abs(value));
686 Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
691 value += fParamSQPar[dim][type][0];
692 value += fParamSQPar[dim][type][1]*z;
693 value += fParamSQPar[dim][type][2]*angle*angle;
694 value += fParamSQPar[dim][type][3]*z*z;
695 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
696 value += fParamSQPar[dim][type][5]*z*angle*angle;
697 value += fParamSQPar[dim][type][6]*z/Qmean;
698 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
699 Float_t valueMean = GetError0Par(dim,type,z,angle);
700 value -= 0.35*0.35*valueMean*valueMean;
701 value = TMath::Sqrt(TMath::Abs(value));
707 Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
709 // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
712 value += fParamRMS0[dim][type][0];
713 value += fParamRMS0[dim][type][1]*z;
714 value += fParamRMS0[dim][type][2]*angle*angle;
715 value = TMath::Sqrt(TMath::Abs(value));
719 Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
721 // calculate mean RMS of cluster - z,angle - pad length scalling
725 if (type==1) length=1;
726 if (type==2) length=1.5;
728 value += fParamRMS1[dim][0];
730 value += fParamRMS1[dim][1];
732 value += fParamRMS1[dim][2]*z;
733 value += fParamRMS1[dim][3]*angle*angle*length*length;
734 value = TMath::Sqrt(TMath::Abs(value));
738 Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
740 // calculate mean RMS of cluster - z,angle, Q dependence
743 value += fParamRMSQ[dim][type][0];
744 value += fParamRMSQ[dim][type][1]*z;
745 value += fParamRMSQ[dim][type][2]*angle*angle;
746 value += fParamRMSQ[dim][type][3]*z/Qmean;
747 value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
748 value = TMath::Sqrt(TMath::Abs(value));
752 Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
754 // calculates RMS of signal shape fluctuation
756 Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
757 Float_t value = fRMSSigmaFit[dim][type][0];
758 value+= fRMSSigmaFit[dim][type][1]*mean;
762 Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
764 // calculates vallue - sigma distortion contribution
768 Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
769 if (rmsL<rmsMeanQ) return value;
771 Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
773 if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
774 //1.5 sigma cut on mean
775 value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
777 if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
778 //3 sigma cut on local
779 value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
782 return TMath::Sqrt(value);
787 void AliTPCClusterParam::FitData(TTree * tree){
789 // make fits for error param and shape param
796 void AliTPCClusterParam::FitResol(TTree * tree){
799 for (Int_t idir=0;idir<2; idir++){
800 for (Int_t itype=0; itype<3; itype++){
804 FitResol0(tree, idir, itype,param0,error0);
805 printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
806 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
807 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
808 for (Int_t ipar=0;ipar<4; ipar++){
809 fParamS0[idir][itype][ipar] = param0[ipar];
810 fErrorS0[idir][itype][ipar] = param0[ipar];
812 // error param with parabolic correction
813 FitResol0Par(tree, idir, itype,param0,error0);
814 printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
815 printf("%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5]);
816 printf("%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5]);
817 for (Int_t ipar=0;ipar<7; ipar++){
818 fParamS0Par[idir][itype][ipar] = param0[ipar];
819 fErrorS0Par[idir][itype][ipar] = param0[ipar];
822 FitResolQ(tree, idir, itype,param0,error0);
823 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
824 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
825 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
826 for (Int_t ipar=0;ipar<6; ipar++){
827 fParamSQ[idir][itype][ipar] = param0[ipar];
828 fErrorSQ[idir][itype][ipar] = param0[ipar];
831 FitResolQPar(tree, idir, itype,param0,error0);
832 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
833 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5],param0[6],param0[7]);
834 printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5],error0[6],error0[7]);
835 for (Int_t ipar=0;ipar<9; ipar++){
836 fParamSQPar[idir][itype][ipar] = param0[ipar];
837 fErrorSQPar[idir][itype][ipar] = param0[ipar];
842 printf("Resol z-scaled\n");
843 for (Int_t idir=0;idir<2; idir++){
846 FitResol1(tree, idir,param0,error0);
847 printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
848 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
849 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
850 for (Int_t ipar=0;ipar<4; ipar++){
851 fParamS1[idir][ipar] = param0[ipar];
852 fErrorS1[idir][ipar] = param0[ipar];
856 for (Int_t idir=0;idir<2; idir++){
857 printf("\nDirection %d\n",idir);
858 printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
859 for (Int_t itype=0; itype<3; itype++){
861 if (itype==1) length=1;
862 if (itype==2) length=1.5;
863 printf("%d\t%f\t%f\t%f\n", itype,fParamS0[idir][itype][0],fParamS0[idir][itype][1]*TMath::Sqrt(length),fParamS0[idir][itype][2]/TMath::Sqrt(length));
870 void AliTPCClusterParam::FitRMS(TTree * tree){
873 for (Int_t idir=0;idir<2; idir++){
874 for (Int_t itype=0; itype<3; itype++){
877 FitRMS0(tree, idir, itype,param0,error0);
878 printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
879 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
880 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
881 for (Int_t ipar=0;ipar<4; ipar++){
882 fParamRMS0[idir][itype][ipar] = param0[ipar];
883 fErrorRMS0[idir][itype][ipar] = param0[ipar];
885 FitRMSQ(tree, idir, itype,param0,error0);
886 printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
887 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
888 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
889 for (Int_t ipar=0;ipar<6; ipar++){
890 fParamRMSQ[idir][itype][ipar] = param0[ipar];
891 fErrorRMSQ[idir][itype][ipar] = param0[ipar];
896 printf("RMS z-scaled\n");
897 for (Int_t idir=0;idir<2; idir++){
900 FitRMS1(tree, idir,param0,error0);
901 printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
902 printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
903 printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
904 for (Int_t ipar=0;ipar<5; ipar++){
905 fParamRMS1[idir][ipar] = param0[ipar];
906 fErrorRMS1[idir][ipar] = param0[ipar];
910 for (Int_t idir=0;idir<2; idir++){
911 printf("\nDirection %d\n",idir);
912 printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
913 for (Int_t itype=0; itype<3; itype++){
915 if (itype==1) length=1;
916 if (itype==2) length=1.5;
917 if (itype==0) printf("%d\t%f\t\t\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length);
918 if (itype>0) printf("%d\t\t\t%f\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length);
924 printf("RMS fluctuation parameterization \n");
925 for (Int_t idir=0;idir<2; idir++){
926 for (Int_t itype=0; itype<3; itype++){
929 FitRMSSigma(tree, idir,itype,param0,error0);
930 printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
931 for (Int_t ipar=0;ipar<2; ipar++){
932 fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
937 // store systematic error end RMS fluctuation parameterization
939 TH1F hratio("hratio","hratio",100,-0.1,0.1);
940 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
941 fErrorRMSSys[0] = hratio.GetRMS();
942 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
943 fErrorRMSSys[1] = hratio.GetRMS();
944 TH1F hratioR("hratioR","hratioR",100,0,0.2);
945 tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
946 fRMSSigmaRatio[0][0]=hratioR.GetMean();
947 fRMSSigmaRatio[0][1]=hratioR.GetRMS();
948 tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
949 fRMSSigmaRatio[1][0]=hratioR.GetMean();
950 fRMSSigmaRatio[1][1]=hratioR.GetRMS();
953 void AliTPCClusterParam::Test(TTree * tree, const char *output){
955 // Draw standard quality histograms to output file
958 TFile f(output,"recreate");
961 // 1D histograms - resolution
963 for (Int_t idim=0; idim<2; idim++){
964 for (Int_t ipad=0; ipad<3; ipad++){
969 sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
970 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
971 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
972 TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
973 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
974 tree->Draw(hexp1,hcut1,"");
977 sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
978 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
979 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
980 TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
981 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
982 tree->Draw(hexp1,hcut1,"");
983 his1DRel0Par.Write();
988 // 2D histograms - resolution
990 for (Int_t idim=0; idim<2; idim++){
991 for (Int_t ipad=0; ipad<3; ipad++){
996 sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
997 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
998 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
999 TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1000 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1001 tree->Draw(hexp1,hcut1,"");
1004 sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1005 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1006 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1007 TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1008 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1009 tree->Draw(hexp1,hcut1,"");
1010 profDRel0Par.Write();
1018 void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1020 // Print param Information
1024 // Error parameterization
1026 printf("\nResolution Scaled factors\n");
1027 printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
1028 printf("Y\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[0][0])),TMath::Sqrt(fParamS1[0][1]),
1029 TMath::Sqrt(fParamS1[0][2]),TMath::Sqrt(fParamS1[0][3]));
1030 for (Int_t ipad=0; ipad<3; ipad++){
1031 Float_t length=0.75;
1032 if (ipad==1) length=1;
1033 if (ipad==2) length=1.5;
1034 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1035 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
1036 TMath::Sqrt(fParamS0[0][ipad][1]*length),
1037 TMath::Sqrt(fParamS0[0][ipad][2]/length),
1038 TMath::Sqrt(fParamS0[0][ipad][3]));
1040 for (Int_t ipad=0; ipad<3; ipad++){
1041 Float_t length=0.75;
1042 if (ipad==1) length=1;
1043 if (ipad==2) length=1.5;
1044 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1045 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
1046 TMath::Sqrt(fParamS0Par[0][ipad][1]*length),
1047 TMath::Sqrt(fParamS0Par[0][ipad][2]/length),
1048 TMath::Sqrt(fParamS0Par[0][ipad][6]));
1050 printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1051 TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1053 for (Int_t ipad=0; ipad<3; ipad++){
1054 Float_t length=0.75;
1055 if (ipad==1) length=1;
1056 if (ipad==2) length=1.5;
1057 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1058 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
1059 TMath::Sqrt(fParamS0[1][ipad][1]*length),
1060 TMath::Sqrt(fParamS0[1][ipad][2]/length),
1061 TMath::Sqrt(fParamS0[1][ipad][3]));
1063 for (Int_t ipad=0; ipad<3; ipad++){
1064 Float_t length=0.75;
1065 if (ipad==1) length=1;
1066 if (ipad==2) length=1.5;
1067 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1068 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][0])),
1069 TMath::Sqrt(fParamS0Par[1][ipad][1]*length),
1070 TMath::Sqrt(fParamS0Par[1][ipad][2]/length),
1071 TMath::Sqrt(fParamS0Par[1][ipad][6]));
1078 printf("\nRMS Scaled factors\n");
1079 printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
1080 printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),TMath::Sqrt(fParamRMS1[0][1]),
1081 TMath::Sqrt(fParamRMS1[0][2]),TMath::Sqrt(fParamRMS1[0][3]),TMath::Sqrt(fParamRMS1[0][4]));
1082 for (Int_t ipad=0; ipad<3; ipad++){
1083 Float_t length=0.75;
1084 if (ipad==1) length=1;
1085 if (ipad==2) length=1.5;
1087 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1088 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1090 TMath::Sqrt(fParamRMS0[0][ipad][1]),
1091 TMath::Sqrt(fParamRMS0[0][ipad][2]/(length*length)),
1092 TMath::Sqrt(fParamRMS0[0][ipad][3]));
1094 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1096 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1097 TMath::Sqrt(fParamRMS0[0][ipad][1]),
1098 TMath::Sqrt(fParamRMS0[0][ipad][2]/(length*length)),
1099 TMath::Sqrt(fParamRMS0[0][ipad][3]));
1103 printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),TMath::Sqrt(fParamRMS1[1][1]),
1104 TMath::Sqrt(fParamRMS1[1][2]),TMath::Sqrt(fParamRMS1[1][3]),TMath::Sqrt(fParamRMS1[1][4]));
1105 for (Int_t ipad=0; ipad<3; ipad++){
1106 Float_t length=0.75;
1107 if (ipad==1) length=1;
1108 if (ipad==2) length=1.5;
1110 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1111 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1113 TMath::Sqrt(fParamRMS0[1][ipad][1]),
1114 TMath::Sqrt(fParamRMS0[1][ipad][2]/(length*length)),
1115 TMath::Sqrt(fParamRMS0[1][ipad][3]));
1117 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1119 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1120 TMath::Sqrt(fParamRMS0[1][ipad][1]),
1121 TMath::Sqrt(fParamRMS0[1][ipad][2]/(length*length)),
1122 TMath::Sqrt(fParamRMS0[1][ipad][3]));