]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCClusterParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Memory leak problems (TLinearFitter)- way around
[u/mrichter/AliRoot.git] / TPC / AliTPCClusterParam.cxx
CommitLineData
12ca5da1 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////////
18// //
dcf3a564 19// TPC cluster error, shape and charge parameterization as function
20// of drift length, and inclination angle //
21//
22// Following notation is used in following
23// Int_t dim 0 - y direction
24// 1 - z direction
25//
26// Int_t type 0 - short pads
27// 1 - medium pads
28// 2 - long pads
29// Float_t z - drift length
30//
d028aade 31// Float_t angle - tangent of inclination angle at given dimension
dcf3a564 32//
33// Implemented parameterization
34//
35//
36// 1. Resolution as function of drift length and inclination angle
37// 1.a) GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle)
38// Simple error parameterization as derived from analytical formula
39// only linear term in drift length and angle^2
40// The formula is valid only with precission +-5%
41// Separate parameterization for differnt pad geometry
42// 1.b) GetError0Par
43// Parabolic term correction - better precision
44//
45// 1.c) GetError1 - JUST FOR Study
46// Similar to GetError1
47// The angular and diffusion effect is scaling with pad length
48// common parameterization for different pad length
49//
96305e49 50// 2. Error parameterization using charge
51// 2.a) GetErrorQ
52// GetError0+
53// adding 1/Q component to diffusion and angluar part
54// 2.b) GetErrorQPar
55// GetError0Par+
56// adding 1/Q component to diffusion and angluar part
57// 2.c) GetErrorQParScaled - Just for study
58// One parameterization for all pad shapes
59// Smaller precission as previous one
dcf3a564 60//
96305e49 61//
d028aade 62// Example how to retrieve the paramterization:
63/*
64 AliCDBManager::Instance()->SetRun(1)
65 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
66 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
67
68 //
69 //
70 AliTPCClusterParam::SetInstance(param);
71 TF1 f1("f1","AliTPCClusterParam::SGetError0Par(1,0,x,0)",0,250);
72
73*/
96305e49 74//
d028aade 75// //
12ca5da1 76///////////////////////////////////////////////////////////////////////////////
77#include "AliTPCClusterParam.h"
78#include "TMath.h"
79#include "TFile.h"
80#include "TTree.h"
81#include <TVectorF.h>
82#include <TLinearFitter.h>
83#include <TH1F.h>
84#include <TProfile2D.h>
0a65832b 85#include <TVectorD.h>
86#include <TObjArray.h>
12ca5da1 87
88ClassImp(AliTPCClusterParam)
89
90
91AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
92
93
94/*
95 Example usage fitting parameterization:
96 TFile fres("resol.root"); //tree with resolution and shape
97 TTree * treeRes =(TTree*)fres.Get("Resol");
98
99 AliTPCClusterParam param;
100 param.SetInstance(&param);
101 param.FitResol(treeRes);
102 param.FitRMS(treeRes);
103 TFile fparam("TPCClusterParam.root","recreate");
104 param.Write("Param");
105 //
106 //
107 TFile fparam("TPCClusterParam.root");
108 AliTPCClusterParam *param2 = (AliTPCClusterParam *) fparam.Get("Param");
109 param2->SetInstance(param2);
110 param2->Test(treeRes);
111
112
113 treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
114
115*/
116
117
118
119
120//_ singleton implementation __________________________________________________
121AliTPCClusterParam* AliTPCClusterParam::Instance()
122{
123 //
124 // Singleton implementation
125 // Returns an instance of this class, it is created if neccessary
126 //
127 if (fgInstance == 0){
128 fgInstance = new AliTPCClusterParam();
129 }
130 return fgInstance;
131}
132
133
f1c2a4a3 134AliTPCClusterParam::AliTPCClusterParam():
135 TObject(),
38caa778 136 fRatio(0),
f1c2a4a3 137 fQNorm(0)
138{
139 //
140 // Default constructor
141 //
142}
38caa778 143
144AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
145 TObject(param),
146 fRatio(0),
147 fQNorm(0)
148{
149 //
150 // copy constructor
151 //
152 memcpy(this, &param,sizeof(AliTPCClusterParam));
153 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
154}
155
156AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
157 //
158 // Assignment operator
159 //
160 if (this != &param) {
161 memcpy(this, &param,sizeof(AliTPCClusterParam));
162 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
163 }
164 return *this;
165}
166
167
f1c2a4a3 168AliTPCClusterParam::~AliTPCClusterParam(){
169 //
170 // destructor
171 //
172 if (fQNorm) fQNorm->Delete();
173 delete fQNorm;
174}
12ca5da1 175
176
177void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
178 //
179 // Fit z - angular dependence of resolution
180 //
181 // Int_t dim=0, type=0;
182 char varVal[100];
183 sprintf(varVal,"Resol:AngleM:Zm");
184 char varErr[100];
185 sprintf(varErr,"Sigma:AngleS:Zs");
186 char varCut[100];
187 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
188 //
189 Int_t entries = tree->Draw(varVal,varCut);
190 Float_t px[10000], py[10000], pz[10000];
191 Float_t ex[10000], ey[10000], ez[10000];
192 //
193 tree->Draw(varErr,varCut);
194 for (Int_t ipoint=0; ipoint<entries; ipoint++){
195 ex[ipoint]= tree->GetV3()[ipoint];
196 ey[ipoint]= tree->GetV2()[ipoint];
197 ez[ipoint]= tree->GetV1()[ipoint];
198 }
199 tree->Draw(varVal,varCut);
200 for (Int_t ipoint=0; ipoint<entries; ipoint++){
201 px[ipoint]= tree->GetV3()[ipoint];
202 py[ipoint]= tree->GetV2()[ipoint];
203 pz[ipoint]= tree->GetV1()[ipoint];
204 }
205
206 //
207 TLinearFitter fitter(3,"hyp2");
208 for (Int_t ipoint=0; ipoint<entries; ipoint++){
209 Float_t val = pz[ipoint]*pz[ipoint];
210 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
211 Double_t x[2];
212 x[0] = px[ipoint];
213 x[1] = py[ipoint]*py[ipoint];
214 fitter.AddPoint(x,val,err);
215 }
216 fitter.Eval();
217 TVectorD param(3);
218 fitter.GetParameters(param);
219 param0[0] = param[0];
220 param0[1] = param[1];
221 param0[2] = param[2];
222 Float_t chi2 = fitter.GetChisquare()/entries;
223 param0[3] = chi2;
224 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
225 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
226 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
227}
228
229
230void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
231 //
232 // Fit z - angular dependence of resolution
233 //
234 // Int_t dim=0, type=0;
235 char varVal[100];
236 sprintf(varVal,"Resol:AngleM:Zm");
237 char varErr[100];
238 sprintf(varErr,"Sigma:AngleS:Zs");
239 char varCut[100];
240 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
241 //
242 Int_t entries = tree->Draw(varVal,varCut);
243 Float_t px[10000], py[10000], pz[10000];
244 Float_t ex[10000], ey[10000], ez[10000];
245 //
246 tree->Draw(varErr,varCut);
247 for (Int_t ipoint=0; ipoint<entries; ipoint++){
248 ex[ipoint]= tree->GetV3()[ipoint];
249 ey[ipoint]= tree->GetV2()[ipoint];
250 ez[ipoint]= tree->GetV1()[ipoint];
251 }
252 tree->Draw(varVal,varCut);
253 for (Int_t ipoint=0; ipoint<entries; ipoint++){
254 px[ipoint]= tree->GetV3()[ipoint];
255 py[ipoint]= tree->GetV2()[ipoint];
256 pz[ipoint]= tree->GetV1()[ipoint];
257 }
258
259 //
260 TLinearFitter fitter(6,"hyp5");
261 for (Int_t ipoint=0; ipoint<entries; ipoint++){
262 Float_t val = pz[ipoint]*pz[ipoint];
263 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
264 Double_t x[6];
265 x[0] = px[ipoint];
266 x[1] = py[ipoint]*py[ipoint];
267 x[2] = x[0]*x[0];
268 x[3] = x[1]*x[1];
269 x[4] = x[0]*x[1];
270 fitter.AddPoint(x,val,err);
271 }
272 fitter.Eval();
273 TVectorD param(6);
274 fitter.GetParameters(param);
275 param0[0] = param[0];
276 param0[1] = param[1];
277 param0[2] = param[2];
278 param0[3] = param[3];
279 param0[4] = param[4];
280 param0[5] = param[5];
281 Float_t chi2 = fitter.GetChisquare()/entries;
282 param0[6] = chi2;
283 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
284 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
285 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
286 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
287 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
288 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
289}
290
291
292
293
294
295void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
296 //
297 // Fit z - angular dependence of resolution - pad length scaling
298 //
299 // Int_t dim=0, type=0;
300 char varVal[100];
301 sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
302 char varErr[100];
303 sprintf(varErr,"Sigma:AngleS:Zs");
304 char varCut[100];
305 sprintf(varCut,"Dim==%d&&QMean<0",dim);
306 //
307 Int_t entries = tree->Draw(varVal,varCut);
308 Float_t px[10000], py[10000], pz[10000];
309 Float_t ex[10000], ey[10000], ez[10000];
310 //
311 tree->Draw(varErr,varCut);
312 for (Int_t ipoint=0; ipoint<entries; ipoint++){
313 ex[ipoint]= tree->GetV3()[ipoint];
314 ey[ipoint]= tree->GetV2()[ipoint];
315 ez[ipoint]= tree->GetV1()[ipoint];
316 }
317 tree->Draw(varVal,varCut);
318 for (Int_t ipoint=0; ipoint<entries; ipoint++){
319 px[ipoint]= tree->GetV3()[ipoint];
320 py[ipoint]= tree->GetV2()[ipoint];
321 pz[ipoint]= tree->GetV1()[ipoint];
322 }
323
324 //
325 TLinearFitter fitter(3,"hyp2");
326 for (Int_t ipoint=0; ipoint<entries; ipoint++){
327 Float_t val = pz[ipoint]*pz[ipoint];
328 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
329 Double_t x[2];
330 x[0] = px[ipoint];
331 x[1] = py[ipoint]*py[ipoint];
332 fitter.AddPoint(x,val,err);
333 }
334 fitter.Eval();
335 TVectorD param(3);
336 fitter.GetParameters(param);
337 param0[0] = param[0];
338 param0[1] = param[1];
339 param0[2] = param[2];
340 Float_t chi2 = fitter.GetChisquare()/entries;
341 param0[3] = chi2;
342 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
343 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
344 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
345}
346
347void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
348 //
349 // Fit z - angular dependence of resolution - Q scaling
350 //
351 // Int_t dim=0, type=0;
352 char varVal[100];
353 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
354 char varVal0[100];
355 sprintf(varVal0,"Resol:AngleM:Zm");
356 //
357 char varErr[100];
358 sprintf(varErr,"Sigma:AngleS:Zs");
359 char varCut[100];
360 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
361 //
362 Int_t entries = tree->Draw(varVal,varCut);
363 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
364 Float_t ex[20000], ey[20000], ez[20000];
365 //
366 tree->Draw(varErr,varCut);
367 for (Int_t ipoint=0; ipoint<entries; ipoint++){
368 ex[ipoint]= tree->GetV3()[ipoint];
369 ey[ipoint]= tree->GetV2()[ipoint];
370 ez[ipoint]= tree->GetV1()[ipoint];
371 }
372 tree->Draw(varVal,varCut);
373 for (Int_t ipoint=0; ipoint<entries; ipoint++){
374 px[ipoint]= tree->GetV3()[ipoint];
375 py[ipoint]= tree->GetV2()[ipoint];
376 pz[ipoint]= tree->GetV1()[ipoint];
377 }
378 tree->Draw(varVal0,varCut);
379 for (Int_t ipoint=0; ipoint<entries; ipoint++){
380 pu[ipoint]= tree->GetV3()[ipoint];
381 pt[ipoint]= tree->GetV2()[ipoint];
382 }
383
384 //
385 TLinearFitter fitter(5,"hyp4");
386 for (Int_t ipoint=0; ipoint<entries; ipoint++){
387 Float_t val = pz[ipoint]*pz[ipoint];
388 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
389 Double_t x[4];
390 x[0] = pu[ipoint];
391 x[1] = pt[ipoint]*pt[ipoint];
392 x[2] = px[ipoint];
393 x[3] = py[ipoint]*py[ipoint];
394 fitter.AddPoint(x,val,err);
395 }
396
397 fitter.Eval();
398 TVectorD param(5);
399 fitter.GetParameters(param);
400 param0[0] = param[0];
401 param0[1] = param[1];
402 param0[2] = param[2];
403 param0[3] = param[3];
404 param0[4] = param[4];
405 Float_t chi2 = fitter.GetChisquare()/entries;
406 param0[5] = chi2;
407 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
408 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
409 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
410 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
411 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
412}
413
414void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
415 //
416 // Fit z - angular dependence of resolution - Q scaling - parabolic correction
417 //
418 // Int_t dim=0, type=0;
419 char varVal[100];
420 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
421 char varVal0[100];
422 sprintf(varVal0,"Resol:AngleM:Zm");
423 //
424 char varErr[100];
425 sprintf(varErr,"Sigma:AngleS:Zs");
426 char varCut[100];
427 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
428 //
429 Int_t entries = tree->Draw(varVal,varCut);
430 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
431 Float_t ex[20000], ey[20000], ez[20000];
432 //
433 tree->Draw(varErr,varCut);
434 for (Int_t ipoint=0; ipoint<entries; ipoint++){
435 ex[ipoint]= tree->GetV3()[ipoint];
436 ey[ipoint]= tree->GetV2()[ipoint];
437 ez[ipoint]= tree->GetV1()[ipoint];
438 }
439 tree->Draw(varVal,varCut);
440 for (Int_t ipoint=0; ipoint<entries; ipoint++){
441 px[ipoint]= tree->GetV3()[ipoint];
442 py[ipoint]= tree->GetV2()[ipoint];
443 pz[ipoint]= tree->GetV1()[ipoint];
444 }
445 tree->Draw(varVal0,varCut);
446 for (Int_t ipoint=0; ipoint<entries; ipoint++){
447 pu[ipoint]= tree->GetV3()[ipoint];
448 pt[ipoint]= tree->GetV2()[ipoint];
449 }
450
451 //
452 TLinearFitter fitter(8,"hyp7");
453 for (Int_t ipoint=0; ipoint<entries; ipoint++){
454 Float_t val = pz[ipoint]*pz[ipoint];
455 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
456 Double_t x[7];
457 x[0] = pu[ipoint];
458 x[1] = pt[ipoint]*pt[ipoint];
459 x[2] = x[0]*x[0];
460 x[3] = x[1]*x[1];
461 x[4] = x[0]*x[1];
462 x[5] = px[ipoint];
463 x[6] = py[ipoint]*py[ipoint];
464 //
465 fitter.AddPoint(x,val,err);
466 }
467
468 fitter.Eval();
469 TVectorD param(8);
470 fitter.GetParameters(param);
471 param0[0] = param[0];
472 param0[1] = param[1];
473 param0[2] = param[2];
474 param0[3] = param[3];
475 param0[4] = param[4];
476 param0[5] = param[5];
477 param0[6] = param[6];
478 param0[7] = param[7];
479
480 Float_t chi2 = fitter.GetChisquare()/entries;
481 param0[8] = chi2;
482 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
483 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
484 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
485 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
486 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
487 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
488 error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
489 error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
490}
491
492
493
494void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
495 //
496 // Fit z - angular dependence of resolution
497 //
498 // Int_t dim=0, type=0;
499 char varVal[100];
500 sprintf(varVal,"RMSm:AngleM:Zm");
501 char varErr[100];
502 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
503 char varCut[100];
504 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
505 //
506 Int_t entries = tree->Draw(varVal,varCut);
507 Float_t px[10000], py[10000], pz[10000];
508 Float_t ex[10000], ey[10000], ez[10000];
509 //
510 tree->Draw(varErr,varCut);
511 for (Int_t ipoint=0; ipoint<entries; ipoint++){
512 ex[ipoint]= tree->GetV3()[ipoint];
513 ey[ipoint]= tree->GetV2()[ipoint];
514 ez[ipoint]= tree->GetV1()[ipoint];
515 }
516 tree->Draw(varVal,varCut);
517 for (Int_t ipoint=0; ipoint<entries; ipoint++){
518 px[ipoint]= tree->GetV3()[ipoint];
519 py[ipoint]= tree->GetV2()[ipoint];
520 pz[ipoint]= tree->GetV1()[ipoint];
521 }
522
523 //
524 TLinearFitter fitter(3,"hyp2");
525 for (Int_t ipoint=0; ipoint<entries; ipoint++){
526 Float_t val = pz[ipoint]*pz[ipoint];
527 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
528 Double_t x[2];
529 x[0] = px[ipoint];
530 x[1] = py[ipoint]*py[ipoint];
531 fitter.AddPoint(x,val,err);
532 }
533 fitter.Eval();
534 TVectorD param(3);
535 fitter.GetParameters(param);
536 param0[0] = param[0];
537 param0[1] = param[1];
538 param0[2] = param[2];
539 Float_t chi2 = fitter.GetChisquare()/entries;
540 param0[3] = chi2;
541 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
542 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
543 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
544}
545
546void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
547 //
548 // Fit z - angular dependence of resolution - pad length scaling
549 //
550 // Int_t dim=0, type=0;
551 char varVal[100];
552 sprintf(varVal,"RMSm:AngleM*Length:Zm");
553 char varErr[100];
554 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
555 char varCut[100];
556 sprintf(varCut,"Dim==%d&&QMean<0",dim);
557 //
558 Int_t entries = tree->Draw(varVal,varCut);
559 Float_t px[10000], py[10000], pz[10000];
560 Float_t type[10000], ey[10000], ez[10000];
561 //
562 tree->Draw(varErr,varCut);
563 for (Int_t ipoint=0; ipoint<entries; ipoint++){
564 type[ipoint] = tree->GetV3()[ipoint];
565 ey[ipoint] = tree->GetV2()[ipoint];
566 ez[ipoint] = tree->GetV1()[ipoint];
567 }
568 tree->Draw(varVal,varCut);
569 for (Int_t ipoint=0; ipoint<entries; ipoint++){
570 px[ipoint]= tree->GetV3()[ipoint];
571 py[ipoint]= tree->GetV2()[ipoint];
572 pz[ipoint]= tree->GetV1()[ipoint];
573 }
574
575 //
576 TLinearFitter fitter(4,"hyp3");
577 for (Int_t ipoint=0; ipoint<entries; ipoint++){
578 Float_t val = pz[ipoint]*pz[ipoint];
579 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
580 Double_t x[3];
581 x[0] = (type[ipoint]<0.5)? 0.:1.;
582 x[1] = px[ipoint];
583 x[2] = py[ipoint]*py[ipoint];
584 fitter.AddPoint(x,val,err);
585 }
586 fitter.Eval();
587 TVectorD param(4);
588 fitter.GetParameters(param);
589 param0[0] = param[0];
590 param0[1] = param[0]+param[1];
591 param0[2] = param[2];
592 param0[3] = param[3];
593 Float_t chi2 = fitter.GetChisquare()/entries;
594 param0[4] = chi2;
595 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
596 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
597 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
598 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
599}
600
601void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
602 //
603 // Fit z - angular dependence of resolution - Q scaling
604 //
605 // Int_t dim=0, type=0;
606 char varVal[100];
607 sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
608 char varVal0[100];
609 sprintf(varVal0,"RMSm:AngleM:Zm");
610 //
611 char varErr[100];
612 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
613 char varCut[100];
614 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
615 //
616 Int_t entries = tree->Draw(varVal,varCut);
617 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
618 Float_t ex[20000], ey[20000], ez[20000];
619 //
620 tree->Draw(varErr,varCut);
621 for (Int_t ipoint=0; ipoint<entries; ipoint++){
622 ex[ipoint]= tree->GetV3()[ipoint];
623 ey[ipoint]= tree->GetV2()[ipoint];
624 ez[ipoint]= tree->GetV1()[ipoint];
625 }
626 tree->Draw(varVal,varCut);
627 for (Int_t ipoint=0; ipoint<entries; ipoint++){
628 px[ipoint]= tree->GetV3()[ipoint];
629 py[ipoint]= tree->GetV2()[ipoint];
630 pz[ipoint]= tree->GetV1()[ipoint];
631 }
632 tree->Draw(varVal0,varCut);
633 for (Int_t ipoint=0; ipoint<entries; ipoint++){
634 pu[ipoint]= tree->GetV3()[ipoint];
635 pt[ipoint]= tree->GetV2()[ipoint];
636 }
637
638 //
639 TLinearFitter fitter(5,"hyp4");
640 for (Int_t ipoint=0; ipoint<entries; ipoint++){
641 Float_t val = pz[ipoint]*pz[ipoint];
642 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
643 Double_t x[4];
644 x[0] = pu[ipoint];
645 x[1] = pt[ipoint]*pt[ipoint];
646 x[2] = px[ipoint];
647 x[3] = py[ipoint]*py[ipoint];
648 fitter.AddPoint(x,val,err);
649 }
650
651 fitter.Eval();
652 TVectorD param(5);
653 fitter.GetParameters(param);
654 param0[0] = param[0];
655 param0[1] = param[1];
656 param0[2] = param[2];
657 param0[3] = param[3];
658 param0[4] = param[4];
659 Float_t chi2 = fitter.GetChisquare()/entries;
660 param0[5] = chi2;
661 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
662 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
663 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
664 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
665 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
666}
667
668
669void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
670 //
671 // Fit z - angular dependence of resolution - Q scaling
672 //
673 // Int_t dim=0, type=0;
674 char varVal[100];
675 sprintf(varVal,"RMSs:RMSm");
676 //
677 char varCut[100];
678 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
679 //
680 Int_t entries = tree->Draw(varVal,varCut);
681 Float_t px[20000], py[20000];
682 //
683 tree->Draw(varVal,varCut);
684 for (Int_t ipoint=0; ipoint<entries; ipoint++){
685 px[ipoint]= tree->GetV2()[ipoint];
686 py[ipoint]= tree->GetV1()[ipoint];
687 }
688 TLinearFitter fitter(2,"pol1");
689 for (Int_t ipoint=0; ipoint<entries; ipoint++){
690 Float_t val = py[ipoint];
691 Float_t err = fRatio*px[ipoint];
692 Double_t x[4];
693 x[0] = px[ipoint];
694 fitter.AddPoint(x,val,err);
695 }
696 fitter.Eval();
697 param0[0]= fitter.GetParameter(0);
698 param0[1]= fitter.GetParameter(1);
699}
700
701
702
703Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
704 //
705 //
706 //
707 Float_t value=0;
708 value += fParamS0[dim][type][0];
709 value += fParamS0[dim][type][1]*z;
710 value += fParamS0[dim][type][2]*angle*angle;
711 value = TMath::Sqrt(TMath::Abs(value));
712 return value;
713}
714
715
716Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
717 //
718 //
719 //
720 Float_t value=0;
721 value += fParamS0Par[dim][type][0];
722 value += fParamS0Par[dim][type][1]*z;
723 value += fParamS0Par[dim][type][2]*angle*angle;
724 value += fParamS0Par[dim][type][3]*z*z;
725 value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
726 value += fParamS0Par[dim][type][5]*z*angle*angle;
727 value = TMath::Sqrt(TMath::Abs(value));
728 return value;
729}
730
731
732
733Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
734 //
735 //
736 //
737 Float_t value=0;
738 Float_t length=0.75;
739 if (type==1) length=1;
740 if (type==2) length=1.5;
741 value += fParamS1[dim][0];
742 value += fParamS1[dim][1]*z/length;
743 value += fParamS1[dim][2]*angle*angle*length;
744 value = TMath::Sqrt(TMath::Abs(value));
745 return value;
746}
747
748Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
749 //
750 //
751 //
752 Float_t value=0;
753 value += fParamSQ[dim][type][0];
754 value += fParamSQ[dim][type][1]*z;
755 value += fParamSQ[dim][type][2]*angle*angle;
756 value += fParamSQ[dim][type][3]*z/Qmean;
757 value += fParamSQ[dim][type][4]*angle*angle/Qmean;
758 value = TMath::Sqrt(TMath::Abs(value));
759 return value;
760
761
762}
763
764Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
765 //
766 //
767 //
768 Float_t value=0;
769 value += fParamSQPar[dim][type][0];
770 value += fParamSQPar[dim][type][1]*z;
771 value += fParamSQPar[dim][type][2]*angle*angle;
772 value += fParamSQPar[dim][type][3]*z*z;
773 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
774 value += fParamSQPar[dim][type][5]*z*angle*angle;
775 value += fParamSQPar[dim][type][6]*z/Qmean;
776 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
777 value = TMath::Sqrt(TMath::Abs(value));
778 return value;
779
780
781}
782
783Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
784 //
785 //
786 //
787 Float_t value=0;
788 value += fParamSQPar[dim][type][0];
789 value += fParamSQPar[dim][type][1]*z;
790 value += fParamSQPar[dim][type][2]*angle*angle;
791 value += fParamSQPar[dim][type][3]*z*z;
792 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
793 value += fParamSQPar[dim][type][5]*z*angle*angle;
794 value += fParamSQPar[dim][type][6]*z/Qmean;
795 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
796 Float_t valueMean = GetError0Par(dim,type,z,angle);
797 value -= 0.35*0.35*valueMean*valueMean;
798 value = TMath::Sqrt(TMath::Abs(value));
799 return value;
800
801
802}
803
804Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
805 //
806 // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
807 //
808 Float_t value=0;
809 value += fParamRMS0[dim][type][0];
810 value += fParamRMS0[dim][type][1]*z;
811 value += fParamRMS0[dim][type][2]*angle*angle;
812 value = TMath::Sqrt(TMath::Abs(value));
813 return value;
814}
815
816Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
817 //
818 // calculate mean RMS of cluster - z,angle - pad length scalling
819 //
820 Float_t value=0;
821 Float_t length=0.75;
822 if (type==1) length=1;
823 if (type==2) length=1.5;
824 if (type==0){
825 value += fParamRMS1[dim][0];
826 }else{
827 value += fParamRMS1[dim][1];
828 }
829 value += fParamRMS1[dim][2]*z;
830 value += fParamRMS1[dim][3]*angle*angle*length*length;
831 value = TMath::Sqrt(TMath::Abs(value));
832 return value;
833}
834
835Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
836 //
837 // calculate mean RMS of cluster - z,angle, Q dependence
838 //
839 Float_t value=0;
840 value += fParamRMSQ[dim][type][0];
841 value += fParamRMSQ[dim][type][1]*z;
842 value += fParamRMSQ[dim][type][2]*angle*angle;
843 value += fParamRMSQ[dim][type][3]*z/Qmean;
844 value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
845 value = TMath::Sqrt(TMath::Abs(value));
846 return value;
847}
848
849Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
850 //
851 // calculates RMS of signal shape fluctuation
852 //
853 Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
854 Float_t value = fRMSSigmaFit[dim][type][0];
855 value+= fRMSSigmaFit[dim][type][1]*mean;
856 return value;
857}
858
859Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
860 //
861 // calculates vallue - sigma distortion contribution
862 //
863 Double_t value =0;
864 //
865 Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
866 if (rmsL<rmsMeanQ) return value;
867 //
868 Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
869 //
870 if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
871 //1.5 sigma cut on mean
872 value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
873 }else{
874 if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
875 //3 sigma cut on local
876 value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
877 }
878 }
8076baa0 879 return TMath::Sqrt(TMath::Abs(value));
12ca5da1 880}
881
882
883
884void AliTPCClusterParam::FitData(TTree * tree){
885 //
886 // make fits for error param and shape param
887 //
888 FitResol(tree);
889 FitRMS(tree);
890
891}
892
893void AliTPCClusterParam::FitResol(TTree * tree){
894 //
895 SetInstance(this);
896 for (Int_t idir=0;idir<2; idir++){
897 for (Int_t itype=0; itype<3; itype++){
898 Float_t param0[10];
899 Float_t error0[10];
900 // model error param
901 FitResol0(tree, idir, itype,param0,error0);
902 printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
903 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
904 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
905 for (Int_t ipar=0;ipar<4; ipar++){
906 fParamS0[idir][itype][ipar] = param0[ipar];
907 fErrorS0[idir][itype][ipar] = param0[ipar];
908 }
909 // error param with parabolic correction
910 FitResol0Par(tree, idir, itype,param0,error0);
911 printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
912 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]);
913 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]);
914 for (Int_t ipar=0;ipar<7; ipar++){
915 fParamS0Par[idir][itype][ipar] = param0[ipar];
916 fErrorS0Par[idir][itype][ipar] = param0[ipar];
917 }
918 //
919 FitResolQ(tree, idir, itype,param0,error0);
920 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
921 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
922 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
923 for (Int_t ipar=0;ipar<6; ipar++){
924 fParamSQ[idir][itype][ipar] = param0[ipar];
925 fErrorSQ[idir][itype][ipar] = param0[ipar];
926 }
927 //
928 FitResolQPar(tree, idir, itype,param0,error0);
929 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
930 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]);
931 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]);
932 for (Int_t ipar=0;ipar<9; ipar++){
933 fParamSQPar[idir][itype][ipar] = param0[ipar];
934 fErrorSQPar[idir][itype][ipar] = param0[ipar];
935 }
936 }
937 }
938 //
939 printf("Resol z-scaled\n");
940 for (Int_t idir=0;idir<2; idir++){
941 Float_t param0[4];
942 Float_t error0[4];
943 FitResol1(tree, idir,param0,error0);
944 printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
945 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
946 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
947 for (Int_t ipar=0;ipar<4; ipar++){
948 fParamS1[idir][ipar] = param0[ipar];
949 fErrorS1[idir][ipar] = param0[ipar];
950 }
951 }
952
953 for (Int_t idir=0;idir<2; idir++){
954 printf("\nDirection %d\n",idir);
955 printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
956 for (Int_t itype=0; itype<3; itype++){
957 Float_t length=0.75;
958 if (itype==1) length=1;
959 if (itype==2) length=1.5;
960 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));
961 }
962 }
963}
964
965
966
967void AliTPCClusterParam::FitRMS(TTree * tree){
968 //
969 SetInstance(this);
970 for (Int_t idir=0;idir<2; idir++){
971 for (Int_t itype=0; itype<3; itype++){
972 Float_t param0[6];
973 Float_t error0[6];
974 FitRMS0(tree, idir, itype,param0,error0);
975 printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
976 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
977 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
978 for (Int_t ipar=0;ipar<4; ipar++){
979 fParamRMS0[idir][itype][ipar] = param0[ipar];
980 fErrorRMS0[idir][itype][ipar] = param0[ipar];
981 }
982 FitRMSQ(tree, idir, itype,param0,error0);
983 printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
984 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
985 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
986 for (Int_t ipar=0;ipar<6; ipar++){
987 fParamRMSQ[idir][itype][ipar] = param0[ipar];
988 fErrorRMSQ[idir][itype][ipar] = param0[ipar];
989 }
990 }
991 }
992 //
993 printf("RMS z-scaled\n");
994 for (Int_t idir=0;idir<2; idir++){
995 Float_t param0[5];
996 Float_t error0[5];
997 FitRMS1(tree, idir,param0,error0);
998 printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
999 printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1000 printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1001 for (Int_t ipar=0;ipar<5; ipar++){
1002 fParamRMS1[idir][ipar] = param0[ipar];
1003 fErrorRMS1[idir][ipar] = param0[ipar];
1004 }
1005 }
1006
1007 for (Int_t idir=0;idir<2; idir++){
1008 printf("\nDirection %d\n",idir);
1009 printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1010 for (Int_t itype=0; itype<3; itype++){
1011 Float_t length=0.75;
1012 if (itype==1) length=1;
1013 if (itype==2) length=1.5;
1014 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);
1015 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);
1016 }
1017 }
1018 //
1019 // Fit RMS sigma
1020 //
1021 printf("RMS fluctuation parameterization \n");
1022 for (Int_t idir=0;idir<2; idir++){
1023 for (Int_t itype=0; itype<3; itype++){
1024 Float_t param0[5];
1025 Float_t error0[5];
1026 FitRMSSigma(tree, idir,itype,param0,error0);
1027 printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1028 for (Int_t ipar=0;ipar<2; ipar++){
1029 fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1030 }
1031 }
1032 }
1033 //
1034 // store systematic error end RMS fluctuation parameterization
1035 //
1036 TH1F hratio("hratio","hratio",100,-0.1,0.1);
1037 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1038 fErrorRMSSys[0] = hratio.GetRMS();
1039 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1040 fErrorRMSSys[1] = hratio.GetRMS();
1041 TH1F hratioR("hratioR","hratioR",100,0,0.2);
1042 tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1043 fRMSSigmaRatio[0][0]=hratioR.GetMean();
1044 fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1045 tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1046 fRMSSigmaRatio[1][0]=hratioR.GetMean();
1047 fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1048}
1049
1050void AliTPCClusterParam::Test(TTree * tree, const char *output){
1051 //
1052 // Draw standard quality histograms to output file
1053 //
1054 SetInstance(this);
1055 TFile f(output,"recreate");
1056 f.cd();
1057 //
1058 // 1D histograms - resolution
1059 //
1060 for (Int_t idim=0; idim<2; idim++){
1061 for (Int_t ipad=0; ipad<3; ipad++){
1062 char hname1[300];
1063 char hcut1[300];
1064 char hexp1[300];
1065 //
1066 sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1067 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1068 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1069 TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1070 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1071 tree->Draw(hexp1,hcut1,"");
1072 his1DRel0.Write();
1073 //
1074 sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1075 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1076 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1077 TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1078 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1079 tree->Draw(hexp1,hcut1,"");
1080 his1DRel0Par.Write();
1081 //
1082 }
1083 }
1084 //
1085 // 2D histograms - resolution
1086 //
1087 for (Int_t idim=0; idim<2; idim++){
1088 for (Int_t ipad=0; ipad<3; ipad++){
1089 char hname1[300];
1090 char hcut1[300];
1091 char hexp1[300];
1092 //
1093 sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1094 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1095 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1096 TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1097 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1098 tree->Draw(hexp1,hcut1,"");
1099 profDRel0.Write();
1100 //
1101 sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1102 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1103 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1104 TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1105 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1106 tree->Draw(hexp1,hcut1,"");
1107 profDRel0Par.Write();
1108 //
1109 }
1110 }
1111}
1112
1113
1114
1115void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1116 //
1117 // Print param Information
1118 //
1119
1120 //
1121 // Error parameterization
1122 //
1123 printf("\nResolution Scaled factors\n");
1124 printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1125 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])),
1126 TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
12ca5da1 1127 for (Int_t ipad=0; ipad<3; ipad++){
1128 Float_t length=0.75;
1129 if (ipad==1) length=1;
1130 if (ipad==2) length=1.5;
1131 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1132 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
8076baa0 1133 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1134 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1135 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
12ca5da1 1136 }
1137 for (Int_t ipad=0; ipad<3; ipad++){
1138 Float_t length=0.75;
1139 if (ipad==1) length=1;
1140 if (ipad==2) length=1.5;
1141 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1142 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
8076baa0 1143 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1144 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1145 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
12ca5da1 1146 }
1147 printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1148 TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1149
1150 for (Int_t ipad=0; ipad<3; ipad++){
1151 Float_t length=0.75;
1152 if (ipad==1) length=1;
1153 if (ipad==2) length=1.5;
1154 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1155 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
8076baa0 1156 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1157 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1158 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
12ca5da1 1159 }
1160 for (Int_t ipad=0; ipad<3; ipad++){
1161 Float_t length=0.75;
1162 if (ipad==1) length=1;
1163 if (ipad==2) length=1.5;
1164 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
8076baa0 1165 TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1166 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1167 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1168 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
12ca5da1 1169 }
1170
1171 //
1172 // RMS scaling
1173 //
1174 printf("\n");
1175 printf("\nRMS Scaled factors\n");
1176 printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1177 printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1178 TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1179 TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1180 TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1181 TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1182 TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
12ca5da1 1183 for (Int_t ipad=0; ipad<3; ipad++){
1184 Float_t length=0.75;
1185 if (ipad==1) length=1;
1186 if (ipad==2) length=1.5;
1187 if (ipad==0){
1188 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1189 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1190 0.,
8076baa0 1191 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1192 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1193 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1194 }else{
1195 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1196 0.,
1197 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
8076baa0 1198 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1199 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1200 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1201 }
1202 }
1203 printf("\n");
8076baa0 1204 printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1205 TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1206 TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1207 TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1208 TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1209 TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
12ca5da1 1210 for (Int_t ipad=0; ipad<3; ipad++){
1211 Float_t length=0.75;
1212 if (ipad==1) length=1;
1213 if (ipad==2) length=1.5;
1214 if (ipad==0){
1215 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1216 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1217 0.,
8076baa0 1218 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1219 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1220 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1221 }else{
1222 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1223 0.,
1224 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
8076baa0 1225 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1226 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1227 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1228 }
1229 }
1230}
1231
1232
1233
1234
1235
0a65832b 1236Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1237 // get Q normalization
1238 // type - 0 Qtot 1 Qmax
1239 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1240 //
f1afff3b 1241 //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
1242
f1c2a4a3 1243 if (fQNorm==0) return 0;
0a65832b 1244 TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1245 if (!norm) return 0;
f1afff3b 1246 TVectorD &no = *norm;
1247 Float_t res = no[0]+
1248 no[1]*dr+
1249 no[2]*ty+
1250 no[3]*tz+
1251 no[4]*dr*ty+
1252 no[5]*dr*tz+
1253 no[6]*ty*tz+
1254 no[7]*dr*dr+
1255 no[8]*ty*ty+
1256 no[9]*tz*tz;
1257 res/=no[0];
0a65832b 1258 return res;
1259}
1260
1261
1262
1263void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
1264 //
1265 // set normalization
1266 //
1267 // type - 0 Qtot 1 Qmax
1268 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1269 //
1270
1271 if (fQNorm==0) fQNorm = new TObjArray(6);
1272 fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1273}