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