]> 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
Adding new class AliTPCcalibTime
[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),
b17540e4 147 fQNorm(0),
148 fQpadTnorm(0), // q pad normalization - Total charge
149 fQpadMnorm(0) // q pad normalization - Max charge
150 //
f1c2a4a3 151{
152 //
153 // Default constructor
154 //
b17540e4 155 fPosQTnorm[0] = 0; fPosQTnorm[1] = 0; fPosQTnorm[2] = 0;
156 fPosQMnorm[0] = 0; fPosQMnorm[1] = 0; fPosQMnorm[2] = 0;
2e5bcb67 157 //
158 fPosYcor[0] = 0; fPosYcor[1] = 0; fPosYcor[2] = 0;
159 fPosZcor[0] = 0; fPosZcor[1] = 0; fPosZcor[2] = 0;
f1c2a4a3 160}
38caa778 161
162AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
163 TObject(param),
164 fRatio(0),
b17540e4 165 fQNorm(0),
166 fQpadTnorm(new TVectorD(*(param.fQpadTnorm))), // q pad normalization - Total charge
167 fQpadMnorm(new TVectorD(*(param.fQpadMnorm))) // q pad normalization - Max charge
168
38caa778 169{
170 //
171 // copy constructor
172 //
173 memcpy(this, &param,sizeof(AliTPCClusterParam));
174 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
b17540e4 175 //
176 if (param.fPosQTnorm[0]){
177 fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
178 fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
179 fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
180 //
181 fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
182 fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
183 fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
184 }
2e5bcb67 185 if (param.fPosYcor[0]){
186 fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
187 fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
188 fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
189 //
190 fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
191 fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
192 fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
193 }
194
38caa778 195}
196
b17540e4 197
38caa778 198AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
199 //
200 // Assignment operator
201 //
202 if (this != &param) {
203 memcpy(this, &param,sizeof(AliTPCClusterParam));
204 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
b17540e4 205 if (param.fPosQTnorm[0]){
206 fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
207 fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
208 fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
209 //
210 fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
211 fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
212 fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
213 }
2e5bcb67 214 if (param.fPosYcor[0]){
215 fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
216 fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
217 fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
218 //
219 fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
220 fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
221 fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
222 }
38caa778 223 }
224 return *this;
225}
226
227
f1c2a4a3 228AliTPCClusterParam::~AliTPCClusterParam(){
229 //
230 // destructor
231 //
232 if (fQNorm) fQNorm->Delete();
233 delete fQNorm;
b17540e4 234 if (fPosQTnorm[0]){
235 delete fPosQTnorm[0];
236 delete fPosQTnorm[1];
237 delete fPosQTnorm[2];
238 //
239 delete fPosQMnorm[0];
240 delete fPosQMnorm[1];
241 delete fPosQMnorm[2];
242 }
2e5bcb67 243 if (fPosYcor[0]){
244 delete fPosYcor[0];
245 delete fPosYcor[1];
246 delete fPosYcor[2];
247 //
248 delete fPosZcor[0];
249 delete fPosZcor[1];
250 delete fPosZcor[2];
251 }
f1c2a4a3 252}
12ca5da1 253
254
255void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
256 //
257 // Fit z - angular dependence of resolution
258 //
259 // Int_t dim=0, type=0;
260 char varVal[100];
261 sprintf(varVal,"Resol:AngleM:Zm");
262 char varErr[100];
263 sprintf(varErr,"Sigma:AngleS:Zs");
264 char varCut[100];
265 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
266 //
267 Int_t entries = tree->Draw(varVal,varCut);
268 Float_t px[10000], py[10000], pz[10000];
269 Float_t ex[10000], ey[10000], ez[10000];
270 //
271 tree->Draw(varErr,varCut);
272 for (Int_t ipoint=0; ipoint<entries; ipoint++){
273 ex[ipoint]= tree->GetV3()[ipoint];
274 ey[ipoint]= tree->GetV2()[ipoint];
275 ez[ipoint]= tree->GetV1()[ipoint];
276 }
277 tree->Draw(varVal,varCut);
278 for (Int_t ipoint=0; ipoint<entries; ipoint++){
279 px[ipoint]= tree->GetV3()[ipoint];
280 py[ipoint]= tree->GetV2()[ipoint];
281 pz[ipoint]= tree->GetV1()[ipoint];
282 }
283
284 //
285 TLinearFitter fitter(3,"hyp2");
286 for (Int_t ipoint=0; ipoint<entries; ipoint++){
287 Float_t val = pz[ipoint]*pz[ipoint];
288 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
289 Double_t x[2];
290 x[0] = px[ipoint];
291 x[1] = py[ipoint]*py[ipoint];
292 fitter.AddPoint(x,val,err);
293 }
294 fitter.Eval();
295 TVectorD param(3);
296 fitter.GetParameters(param);
297 param0[0] = param[0];
298 param0[1] = param[1];
299 param0[2] = param[2];
300 Float_t chi2 = fitter.GetChisquare()/entries;
301 param0[3] = chi2;
302 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
303 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
304 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
305}
306
307
308void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
309 //
310 // Fit z - angular dependence of resolution
311 //
312 // Int_t dim=0, type=0;
313 char varVal[100];
314 sprintf(varVal,"Resol:AngleM:Zm");
315 char varErr[100];
316 sprintf(varErr,"Sigma:AngleS:Zs");
317 char varCut[100];
318 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
319 //
320 Int_t entries = tree->Draw(varVal,varCut);
321 Float_t px[10000], py[10000], pz[10000];
322 Float_t ex[10000], ey[10000], ez[10000];
323 //
324 tree->Draw(varErr,varCut);
325 for (Int_t ipoint=0; ipoint<entries; ipoint++){
326 ex[ipoint]= tree->GetV3()[ipoint];
327 ey[ipoint]= tree->GetV2()[ipoint];
328 ez[ipoint]= tree->GetV1()[ipoint];
329 }
330 tree->Draw(varVal,varCut);
331 for (Int_t ipoint=0; ipoint<entries; ipoint++){
332 px[ipoint]= tree->GetV3()[ipoint];
333 py[ipoint]= tree->GetV2()[ipoint];
334 pz[ipoint]= tree->GetV1()[ipoint];
335 }
336
337 //
338 TLinearFitter fitter(6,"hyp5");
339 for (Int_t ipoint=0; ipoint<entries; ipoint++){
340 Float_t val = pz[ipoint]*pz[ipoint];
341 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
342 Double_t x[6];
343 x[0] = px[ipoint];
344 x[1] = py[ipoint]*py[ipoint];
345 x[2] = x[0]*x[0];
346 x[3] = x[1]*x[1];
347 x[4] = x[0]*x[1];
348 fitter.AddPoint(x,val,err);
349 }
350 fitter.Eval();
351 TVectorD param(6);
352 fitter.GetParameters(param);
353 param0[0] = param[0];
354 param0[1] = param[1];
355 param0[2] = param[2];
356 param0[3] = param[3];
357 param0[4] = param[4];
358 param0[5] = param[5];
359 Float_t chi2 = fitter.GetChisquare()/entries;
360 param0[6] = chi2;
361 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
362 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
363 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
364 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
365 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
366 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
367}
368
369
370
371
372
373void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
374 //
375 // Fit z - angular dependence of resolution - pad length scaling
376 //
377 // Int_t dim=0, type=0;
378 char varVal[100];
379 sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
380 char varErr[100];
381 sprintf(varErr,"Sigma:AngleS:Zs");
382 char varCut[100];
383 sprintf(varCut,"Dim==%d&&QMean<0",dim);
384 //
385 Int_t entries = tree->Draw(varVal,varCut);
386 Float_t px[10000], py[10000], pz[10000];
387 Float_t ex[10000], ey[10000], ez[10000];
388 //
389 tree->Draw(varErr,varCut);
390 for (Int_t ipoint=0; ipoint<entries; ipoint++){
391 ex[ipoint]= tree->GetV3()[ipoint];
392 ey[ipoint]= tree->GetV2()[ipoint];
393 ez[ipoint]= tree->GetV1()[ipoint];
394 }
395 tree->Draw(varVal,varCut);
396 for (Int_t ipoint=0; ipoint<entries; ipoint++){
397 px[ipoint]= tree->GetV3()[ipoint];
398 py[ipoint]= tree->GetV2()[ipoint];
399 pz[ipoint]= tree->GetV1()[ipoint];
400 }
401
402 //
403 TLinearFitter fitter(3,"hyp2");
404 for (Int_t ipoint=0; ipoint<entries; ipoint++){
405 Float_t val = pz[ipoint]*pz[ipoint];
406 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
407 Double_t x[2];
408 x[0] = px[ipoint];
409 x[1] = py[ipoint]*py[ipoint];
410 fitter.AddPoint(x,val,err);
411 }
412 fitter.Eval();
413 TVectorD param(3);
414 fitter.GetParameters(param);
415 param0[0] = param[0];
416 param0[1] = param[1];
417 param0[2] = param[2];
418 Float_t chi2 = fitter.GetChisquare()/entries;
419 param0[3] = chi2;
420 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
421 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
422 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
423}
424
425void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
426 //
427 // Fit z - angular dependence of resolution - Q scaling
428 //
429 // Int_t dim=0, type=0;
430 char varVal[100];
431 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
432 char varVal0[100];
433 sprintf(varVal0,"Resol:AngleM:Zm");
434 //
435 char varErr[100];
436 sprintf(varErr,"Sigma:AngleS:Zs");
437 char varCut[100];
438 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
439 //
440 Int_t entries = tree->Draw(varVal,varCut);
441 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
442 Float_t ex[20000], ey[20000], ez[20000];
443 //
444 tree->Draw(varErr,varCut);
445 for (Int_t ipoint=0; ipoint<entries; ipoint++){
446 ex[ipoint]= tree->GetV3()[ipoint];
447 ey[ipoint]= tree->GetV2()[ipoint];
448 ez[ipoint]= tree->GetV1()[ipoint];
449 }
450 tree->Draw(varVal,varCut);
451 for (Int_t ipoint=0; ipoint<entries; ipoint++){
452 px[ipoint]= tree->GetV3()[ipoint];
453 py[ipoint]= tree->GetV2()[ipoint];
454 pz[ipoint]= tree->GetV1()[ipoint];
455 }
456 tree->Draw(varVal0,varCut);
457 for (Int_t ipoint=0; ipoint<entries; ipoint++){
458 pu[ipoint]= tree->GetV3()[ipoint];
459 pt[ipoint]= tree->GetV2()[ipoint];
460 }
461
462 //
463 TLinearFitter fitter(5,"hyp4");
464 for (Int_t ipoint=0; ipoint<entries; ipoint++){
465 Float_t val = pz[ipoint]*pz[ipoint];
466 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
467 Double_t x[4];
468 x[0] = pu[ipoint];
469 x[1] = pt[ipoint]*pt[ipoint];
470 x[2] = px[ipoint];
471 x[3] = py[ipoint]*py[ipoint];
472 fitter.AddPoint(x,val,err);
473 }
474
475 fitter.Eval();
476 TVectorD param(5);
477 fitter.GetParameters(param);
478 param0[0] = param[0];
479 param0[1] = param[1];
480 param0[2] = param[2];
481 param0[3] = param[3];
482 param0[4] = param[4];
483 Float_t chi2 = fitter.GetChisquare()/entries;
484 param0[5] = chi2;
485 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
486 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
487 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
488 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
489 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
490}
491
492void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
493 //
494 // Fit z - angular dependence of resolution - Q scaling - parabolic correction
495 //
496 // Int_t dim=0, type=0;
497 char varVal[100];
498 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
499 char varVal0[100];
500 sprintf(varVal0,"Resol:AngleM:Zm");
501 //
502 char varErr[100];
503 sprintf(varErr,"Sigma:AngleS:Zs");
504 char varCut[100];
505 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
506 //
507 Int_t entries = tree->Draw(varVal,varCut);
508 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
509 Float_t ex[20000], ey[20000], ez[20000];
510 //
511 tree->Draw(varErr,varCut);
512 for (Int_t ipoint=0; ipoint<entries; ipoint++){
513 ex[ipoint]= tree->GetV3()[ipoint];
514 ey[ipoint]= tree->GetV2()[ipoint];
515 ez[ipoint]= tree->GetV1()[ipoint];
516 }
517 tree->Draw(varVal,varCut);
518 for (Int_t ipoint=0; ipoint<entries; ipoint++){
519 px[ipoint]= tree->GetV3()[ipoint];
520 py[ipoint]= tree->GetV2()[ipoint];
521 pz[ipoint]= tree->GetV1()[ipoint];
522 }
523 tree->Draw(varVal0,varCut);
524 for (Int_t ipoint=0; ipoint<entries; ipoint++){
525 pu[ipoint]= tree->GetV3()[ipoint];
526 pt[ipoint]= tree->GetV2()[ipoint];
527 }
528
529 //
530 TLinearFitter fitter(8,"hyp7");
531 for (Int_t ipoint=0; ipoint<entries; ipoint++){
532 Float_t val = pz[ipoint]*pz[ipoint];
533 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
534 Double_t x[7];
535 x[0] = pu[ipoint];
536 x[1] = pt[ipoint]*pt[ipoint];
537 x[2] = x[0]*x[0];
538 x[3] = x[1]*x[1];
539 x[4] = x[0]*x[1];
540 x[5] = px[ipoint];
541 x[6] = py[ipoint]*py[ipoint];
542 //
543 fitter.AddPoint(x,val,err);
544 }
545
546 fitter.Eval();
547 TVectorD param(8);
548 fitter.GetParameters(param);
549 param0[0] = param[0];
550 param0[1] = param[1];
551 param0[2] = param[2];
552 param0[3] = param[3];
553 param0[4] = param[4];
554 param0[5] = param[5];
555 param0[6] = param[6];
556 param0[7] = param[7];
557
558 Float_t chi2 = fitter.GetChisquare()/entries;
559 param0[8] = chi2;
560 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
561 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
562 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
563 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
564 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
565 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
566 error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
567 error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
568}
569
570
571
572void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
573 //
574 // Fit z - angular dependence of resolution
575 //
576 // Int_t dim=0, type=0;
577 char varVal[100];
578 sprintf(varVal,"RMSm:AngleM:Zm");
579 char varErr[100];
580 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
581 char varCut[100];
582 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
583 //
584 Int_t entries = tree->Draw(varVal,varCut);
585 Float_t px[10000], py[10000], pz[10000];
586 Float_t ex[10000], ey[10000], ez[10000];
587 //
588 tree->Draw(varErr,varCut);
589 for (Int_t ipoint=0; ipoint<entries; ipoint++){
590 ex[ipoint]= tree->GetV3()[ipoint];
591 ey[ipoint]= tree->GetV2()[ipoint];
592 ez[ipoint]= tree->GetV1()[ipoint];
593 }
594 tree->Draw(varVal,varCut);
595 for (Int_t ipoint=0; ipoint<entries; ipoint++){
596 px[ipoint]= tree->GetV3()[ipoint];
597 py[ipoint]= tree->GetV2()[ipoint];
598 pz[ipoint]= tree->GetV1()[ipoint];
599 }
600
601 //
602 TLinearFitter fitter(3,"hyp2");
603 for (Int_t ipoint=0; ipoint<entries; ipoint++){
604 Float_t val = pz[ipoint]*pz[ipoint];
605 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
606 Double_t x[2];
607 x[0] = px[ipoint];
608 x[1] = py[ipoint]*py[ipoint];
609 fitter.AddPoint(x,val,err);
610 }
611 fitter.Eval();
612 TVectorD param(3);
613 fitter.GetParameters(param);
614 param0[0] = param[0];
615 param0[1] = param[1];
616 param0[2] = param[2];
617 Float_t chi2 = fitter.GetChisquare()/entries;
618 param0[3] = chi2;
619 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
620 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
621 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
622}
623
624void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
625 //
626 // Fit z - angular dependence of resolution - pad length scaling
627 //
628 // Int_t dim=0, type=0;
629 char varVal[100];
630 sprintf(varVal,"RMSm:AngleM*Length:Zm");
631 char varErr[100];
632 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
633 char varCut[100];
634 sprintf(varCut,"Dim==%d&&QMean<0",dim);
635 //
636 Int_t entries = tree->Draw(varVal,varCut);
637 Float_t px[10000], py[10000], pz[10000];
638 Float_t type[10000], ey[10000], ez[10000];
639 //
640 tree->Draw(varErr,varCut);
641 for (Int_t ipoint=0; ipoint<entries; ipoint++){
642 type[ipoint] = tree->GetV3()[ipoint];
643 ey[ipoint] = tree->GetV2()[ipoint];
644 ez[ipoint] = tree->GetV1()[ipoint];
645 }
646 tree->Draw(varVal,varCut);
647 for (Int_t ipoint=0; ipoint<entries; ipoint++){
648 px[ipoint]= tree->GetV3()[ipoint];
649 py[ipoint]= tree->GetV2()[ipoint];
650 pz[ipoint]= tree->GetV1()[ipoint];
651 }
652
653 //
654 TLinearFitter fitter(4,"hyp3");
655 for (Int_t ipoint=0; ipoint<entries; ipoint++){
656 Float_t val = pz[ipoint]*pz[ipoint];
657 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
658 Double_t x[3];
659 x[0] = (type[ipoint]<0.5)? 0.:1.;
660 x[1] = px[ipoint];
661 x[2] = py[ipoint]*py[ipoint];
662 fitter.AddPoint(x,val,err);
663 }
664 fitter.Eval();
665 TVectorD param(4);
666 fitter.GetParameters(param);
667 param0[0] = param[0];
668 param0[1] = param[0]+param[1];
669 param0[2] = param[2];
670 param0[3] = param[3];
671 Float_t chi2 = fitter.GetChisquare()/entries;
672 param0[4] = chi2;
673 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
674 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
675 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
676 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
677}
678
679void AliTPCClusterParam::FitRMSQ(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,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
686 char varVal0[100];
687 sprintf(varVal0,"RMSm:AngleM:Zm");
688 //
689 char varErr[100];
690 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
691 char varCut[100];
692 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
693 //
694 Int_t entries = tree->Draw(varVal,varCut);
695 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
696 Float_t ex[20000], ey[20000], ez[20000];
697 //
698 tree->Draw(varErr,varCut);
699 for (Int_t ipoint=0; ipoint<entries; ipoint++){
700 ex[ipoint]= tree->GetV3()[ipoint];
701 ey[ipoint]= tree->GetV2()[ipoint];
702 ez[ipoint]= tree->GetV1()[ipoint];
703 }
704 tree->Draw(varVal,varCut);
705 for (Int_t ipoint=0; ipoint<entries; ipoint++){
706 px[ipoint]= tree->GetV3()[ipoint];
707 py[ipoint]= tree->GetV2()[ipoint];
708 pz[ipoint]= tree->GetV1()[ipoint];
709 }
710 tree->Draw(varVal0,varCut);
711 for (Int_t ipoint=0; ipoint<entries; ipoint++){
712 pu[ipoint]= tree->GetV3()[ipoint];
713 pt[ipoint]= tree->GetV2()[ipoint];
714 }
715
716 //
717 TLinearFitter fitter(5,"hyp4");
718 for (Int_t ipoint=0; ipoint<entries; ipoint++){
719 Float_t val = pz[ipoint]*pz[ipoint];
720 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
721 Double_t x[4];
722 x[0] = pu[ipoint];
723 x[1] = pt[ipoint]*pt[ipoint];
724 x[2] = px[ipoint];
725 x[3] = py[ipoint]*py[ipoint];
726 fitter.AddPoint(x,val,err);
727 }
728
729 fitter.Eval();
730 TVectorD param(5);
731 fitter.GetParameters(param);
732 param0[0] = param[0];
733 param0[1] = param[1];
734 param0[2] = param[2];
735 param0[3] = param[3];
736 param0[4] = param[4];
737 Float_t chi2 = fitter.GetChisquare()/entries;
738 param0[5] = chi2;
739 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
740 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
741 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
742 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
743 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
744}
745
746
747void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
748 //
749 // Fit z - angular dependence of resolution - Q scaling
750 //
751 // Int_t dim=0, type=0;
752 char varVal[100];
753 sprintf(varVal,"RMSs:RMSm");
754 //
755 char varCut[100];
756 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
757 //
758 Int_t entries = tree->Draw(varVal,varCut);
759 Float_t px[20000], py[20000];
760 //
761 tree->Draw(varVal,varCut);
762 for (Int_t ipoint=0; ipoint<entries; ipoint++){
763 px[ipoint]= tree->GetV2()[ipoint];
764 py[ipoint]= tree->GetV1()[ipoint];
765 }
766 TLinearFitter fitter(2,"pol1");
767 for (Int_t ipoint=0; ipoint<entries; ipoint++){
768 Float_t val = py[ipoint];
769 Float_t err = fRatio*px[ipoint];
770 Double_t x[4];
771 x[0] = px[ipoint];
236a0d03 772 if (err>0) fitter.AddPoint(x,val,err);
12ca5da1 773 }
774 fitter.Eval();
775 param0[0]= fitter.GetParameter(0);
776 param0[1]= fitter.GetParameter(1);
777}
778
779
780
781Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
782 //
783 //
784 //
785 Float_t value=0;
786 value += fParamS0[dim][type][0];
787 value += fParamS0[dim][type][1]*z;
788 value += fParamS0[dim][type][2]*angle*angle;
789 value = TMath::Sqrt(TMath::Abs(value));
790 return value;
791}
792
793
794Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
795 //
796 //
797 //
798 Float_t value=0;
799 value += fParamS0Par[dim][type][0];
800 value += fParamS0Par[dim][type][1]*z;
801 value += fParamS0Par[dim][type][2]*angle*angle;
802 value += fParamS0Par[dim][type][3]*z*z;
803 value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
804 value += fParamS0Par[dim][type][5]*z*angle*angle;
805 value = TMath::Sqrt(TMath::Abs(value));
806 return value;
807}
808
809
810
811Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
812 //
813 //
814 //
815 Float_t value=0;
816 Float_t length=0.75;
817 if (type==1) length=1;
818 if (type==2) length=1.5;
819 value += fParamS1[dim][0];
820 value += fParamS1[dim][1]*z/length;
821 value += fParamS1[dim][2]*angle*angle*length;
822 value = TMath::Sqrt(TMath::Abs(value));
823 return value;
824}
825
826Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
827 //
828 //
829 //
830 Float_t value=0;
831 value += fParamSQ[dim][type][0];
832 value += fParamSQ[dim][type][1]*z;
833 value += fParamSQ[dim][type][2]*angle*angle;
834 value += fParamSQ[dim][type][3]*z/Qmean;
835 value += fParamSQ[dim][type][4]*angle*angle/Qmean;
836 value = TMath::Sqrt(TMath::Abs(value));
837 return value;
838
839
840}
841
842Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
843 //
844 //
845 //
846 Float_t value=0;
847 value += fParamSQPar[dim][type][0];
848 value += fParamSQPar[dim][type][1]*z;
849 value += fParamSQPar[dim][type][2]*angle*angle;
850 value += fParamSQPar[dim][type][3]*z*z;
851 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
852 value += fParamSQPar[dim][type][5]*z*angle*angle;
853 value += fParamSQPar[dim][type][6]*z/Qmean;
854 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
855 value = TMath::Sqrt(TMath::Abs(value));
856 return value;
857
858
859}
860
861Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
862 //
863 //
864 //
865 Float_t value=0;
866 value += fParamSQPar[dim][type][0];
867 value += fParamSQPar[dim][type][1]*z;
868 value += fParamSQPar[dim][type][2]*angle*angle;
869 value += fParamSQPar[dim][type][3]*z*z;
870 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
871 value += fParamSQPar[dim][type][5]*z*angle*angle;
872 value += fParamSQPar[dim][type][6]*z/Qmean;
873 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
874 Float_t valueMean = GetError0Par(dim,type,z,angle);
875 value -= 0.35*0.35*valueMean*valueMean;
876 value = TMath::Sqrt(TMath::Abs(value));
877 return value;
878
879
880}
881
882Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
883 //
884 // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
885 //
886 Float_t value=0;
887 value += fParamRMS0[dim][type][0];
888 value += fParamRMS0[dim][type][1]*z;
889 value += fParamRMS0[dim][type][2]*angle*angle;
890 value = TMath::Sqrt(TMath::Abs(value));
891 return value;
892}
893
894Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
895 //
896 // calculate mean RMS of cluster - z,angle - pad length scalling
897 //
898 Float_t value=0;
899 Float_t length=0.75;
900 if (type==1) length=1;
901 if (type==2) length=1.5;
902 if (type==0){
903 value += fParamRMS1[dim][0];
904 }else{
905 value += fParamRMS1[dim][1];
906 }
907 value += fParamRMS1[dim][2]*z;
908 value += fParamRMS1[dim][3]*angle*angle*length*length;
909 value = TMath::Sqrt(TMath::Abs(value));
910 return value;
911}
912
913Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
914 //
915 // calculate mean RMS of cluster - z,angle, Q dependence
916 //
917 Float_t value=0;
918 value += fParamRMSQ[dim][type][0];
919 value += fParamRMSQ[dim][type][1]*z;
920 value += fParamRMSQ[dim][type][2]*angle*angle;
921 value += fParamRMSQ[dim][type][3]*z/Qmean;
922 value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
923 value = TMath::Sqrt(TMath::Abs(value));
924 return value;
925}
926
927Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
928 //
929 // calculates RMS of signal shape fluctuation
930 //
931 Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
932 Float_t value = fRMSSigmaFit[dim][type][0];
933 value+= fRMSSigmaFit[dim][type][1]*mean;
934 return value;
935}
936
937Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
938 //
939 // calculates vallue - sigma distortion contribution
940 //
941 Double_t value =0;
942 //
943 Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
944 if (rmsL<rmsMeanQ) return value;
945 //
946 Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
947 //
948 if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
949 //1.5 sigma cut on mean
950 value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
951 }else{
952 if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
953 //3 sigma cut on local
954 value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
955 }
956 }
8076baa0 957 return TMath::Sqrt(TMath::Abs(value));
12ca5da1 958}
959
960
961
962void AliTPCClusterParam::FitData(TTree * tree){
963 //
964 // make fits for error param and shape param
965 //
966 FitResol(tree);
967 FitRMS(tree);
968
969}
970
971void AliTPCClusterParam::FitResol(TTree * tree){
972 //
973 SetInstance(this);
974 for (Int_t idir=0;idir<2; idir++){
975 for (Int_t itype=0; itype<3; itype++){
976 Float_t param0[10];
977 Float_t error0[10];
978 // model error param
979 FitResol0(tree, idir, itype,param0,error0);
980 printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
981 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
982 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
983 for (Int_t ipar=0;ipar<4; ipar++){
984 fParamS0[idir][itype][ipar] = param0[ipar];
985 fErrorS0[idir][itype][ipar] = param0[ipar];
986 }
987 // error param with parabolic correction
988 FitResol0Par(tree, idir, itype,param0,error0);
989 printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
990 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]);
991 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]);
992 for (Int_t ipar=0;ipar<7; ipar++){
993 fParamS0Par[idir][itype][ipar] = param0[ipar];
994 fErrorS0Par[idir][itype][ipar] = param0[ipar];
995 }
996 //
997 FitResolQ(tree, idir, itype,param0,error0);
998 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
999 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1000 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1001 for (Int_t ipar=0;ipar<6; ipar++){
1002 fParamSQ[idir][itype][ipar] = param0[ipar];
1003 fErrorSQ[idir][itype][ipar] = param0[ipar];
1004 }
1005 //
1006 FitResolQPar(tree, idir, itype,param0,error0);
1007 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
1008 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]);
1009 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]);
1010 for (Int_t ipar=0;ipar<9; ipar++){
1011 fParamSQPar[idir][itype][ipar] = param0[ipar];
1012 fErrorSQPar[idir][itype][ipar] = param0[ipar];
1013 }
1014 }
1015 }
1016 //
1017 printf("Resol z-scaled\n");
1018 for (Int_t idir=0;idir<2; idir++){
1019 Float_t param0[4];
1020 Float_t error0[4];
1021 FitResol1(tree, idir,param0,error0);
1022 printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
1023 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1024 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1025 for (Int_t ipar=0;ipar<4; ipar++){
1026 fParamS1[idir][ipar] = param0[ipar];
1027 fErrorS1[idir][ipar] = param0[ipar];
1028 }
1029 }
1030
1031 for (Int_t idir=0;idir<2; idir++){
1032 printf("\nDirection %d\n",idir);
1033 printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1034 for (Int_t itype=0; itype<3; itype++){
1035 Float_t length=0.75;
1036 if (itype==1) length=1;
1037 if (itype==2) length=1.5;
1038 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));
1039 }
1040 }
1041}
1042
1043
1044
1045void AliTPCClusterParam::FitRMS(TTree * tree){
1046 //
1047 SetInstance(this);
1048 for (Int_t idir=0;idir<2; idir++){
1049 for (Int_t itype=0; itype<3; itype++){
1050 Float_t param0[6];
1051 Float_t error0[6];
1052 FitRMS0(tree, idir, itype,param0,error0);
1053 printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1054 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1055 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1056 for (Int_t ipar=0;ipar<4; ipar++){
1057 fParamRMS0[idir][itype][ipar] = param0[ipar];
1058 fErrorRMS0[idir][itype][ipar] = param0[ipar];
1059 }
1060 FitRMSQ(tree, idir, itype,param0,error0);
1061 printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1062 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1063 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1064 for (Int_t ipar=0;ipar<6; ipar++){
1065 fParamRMSQ[idir][itype][ipar] = param0[ipar];
1066 fErrorRMSQ[idir][itype][ipar] = param0[ipar];
1067 }
1068 }
1069 }
1070 //
1071 printf("RMS z-scaled\n");
1072 for (Int_t idir=0;idir<2; idir++){
1073 Float_t param0[5];
1074 Float_t error0[5];
1075 FitRMS1(tree, idir,param0,error0);
1076 printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1077 printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1078 printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1079 for (Int_t ipar=0;ipar<5; ipar++){
1080 fParamRMS1[idir][ipar] = param0[ipar];
1081 fErrorRMS1[idir][ipar] = param0[ipar];
1082 }
1083 }
1084
1085 for (Int_t idir=0;idir<2; idir++){
1086 printf("\nDirection %d\n",idir);
1087 printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1088 for (Int_t itype=0; itype<3; itype++){
1089 Float_t length=0.75;
1090 if (itype==1) length=1;
1091 if (itype==2) length=1.5;
1092 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);
1093 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);
1094 }
1095 }
1096 //
1097 // Fit RMS sigma
1098 //
1099 printf("RMS fluctuation parameterization \n");
1100 for (Int_t idir=0;idir<2; idir++){
1101 for (Int_t itype=0; itype<3; itype++){
1102 Float_t param0[5];
1103 Float_t error0[5];
1104 FitRMSSigma(tree, idir,itype,param0,error0);
1105 printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1106 for (Int_t ipar=0;ipar<2; ipar++){
1107 fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1108 }
1109 }
1110 }
1111 //
1112 // store systematic error end RMS fluctuation parameterization
1113 //
1114 TH1F hratio("hratio","hratio",100,-0.1,0.1);
1115 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1116 fErrorRMSSys[0] = hratio.GetRMS();
1117 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1118 fErrorRMSSys[1] = hratio.GetRMS();
1119 TH1F hratioR("hratioR","hratioR",100,0,0.2);
1120 tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1121 fRMSSigmaRatio[0][0]=hratioR.GetMean();
1122 fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1123 tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1124 fRMSSigmaRatio[1][0]=hratioR.GetMean();
1125 fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1126}
1127
1128void AliTPCClusterParam::Test(TTree * tree, const char *output){
1129 //
1130 // Draw standard quality histograms to output file
1131 //
1132 SetInstance(this);
1133 TFile f(output,"recreate");
1134 f.cd();
1135 //
1136 // 1D histograms - resolution
1137 //
1138 for (Int_t idim=0; idim<2; idim++){
1139 for (Int_t ipad=0; ipad<3; ipad++){
1140 char hname1[300];
1141 char hcut1[300];
1142 char hexp1[300];
1143 //
1144 sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1145 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1146 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1147 TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1148 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1149 tree->Draw(hexp1,hcut1,"");
1150 his1DRel0.Write();
1151 //
1152 sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1153 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1154 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1155 TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1156 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1157 tree->Draw(hexp1,hcut1,"");
1158 his1DRel0Par.Write();
1159 //
1160 }
1161 }
1162 //
1163 // 2D histograms - resolution
1164 //
1165 for (Int_t idim=0; idim<2; idim++){
1166 for (Int_t ipad=0; ipad<3; ipad++){
1167 char hname1[300];
1168 char hcut1[300];
1169 char hexp1[300];
1170 //
1171 sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1172 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1173 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1174 TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1175 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1176 tree->Draw(hexp1,hcut1,"");
1177 profDRel0.Write();
1178 //
1179 sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1180 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1181 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1182 TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1183 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1184 tree->Draw(hexp1,hcut1,"");
1185 profDRel0Par.Write();
1186 //
1187 }
1188 }
1189}
1190
1191
1192
1193void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1194 //
1195 // Print param Information
1196 //
1197
1198 //
1199 // Error parameterization
1200 //
1201 printf("\nResolution Scaled factors\n");
1202 printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1203 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])),
1204 TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
12ca5da1 1205 for (Int_t ipad=0; ipad<3; ipad++){
1206 Float_t length=0.75;
1207 if (ipad==1) length=1;
1208 if (ipad==2) length=1.5;
1209 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1210 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
8076baa0 1211 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1212 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1213 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
12ca5da1 1214 }
1215 for (Int_t ipad=0; ipad<3; ipad++){
1216 Float_t length=0.75;
1217 if (ipad==1) length=1;
1218 if (ipad==2) length=1.5;
1219 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1220 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
8076baa0 1221 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1222 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1223 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
12ca5da1 1224 }
1225 printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1226 TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1227
1228 for (Int_t ipad=0; ipad<3; ipad++){
1229 Float_t length=0.75;
1230 if (ipad==1) length=1;
1231 if (ipad==2) length=1.5;
1232 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1233 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
8076baa0 1234 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1235 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1236 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
12ca5da1 1237 }
1238 for (Int_t ipad=0; ipad<3; ipad++){
1239 Float_t length=0.75;
1240 if (ipad==1) length=1;
1241 if (ipad==2) length=1.5;
1242 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
8076baa0 1243 TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1244 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1245 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1246 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
12ca5da1 1247 }
1248
1249 //
1250 // RMS scaling
1251 //
1252 printf("\n");
1253 printf("\nRMS Scaled factors\n");
1254 printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1255 printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1256 TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1257 TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1258 TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1259 TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1260 TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
12ca5da1 1261 for (Int_t ipad=0; ipad<3; ipad++){
1262 Float_t length=0.75;
1263 if (ipad==1) length=1;
1264 if (ipad==2) length=1.5;
1265 if (ipad==0){
1266 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1267 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1268 0.,
8076baa0 1269 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1270 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1271 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1272 }else{
1273 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1274 0.,
1275 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
8076baa0 1276 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1277 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1278 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1279 }
1280 }
1281 printf("\n");
8076baa0 1282 printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1283 TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1284 TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1285 TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1286 TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1287 TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
12ca5da1 1288 for (Int_t ipad=0; ipad<3; ipad++){
1289 Float_t length=0.75;
1290 if (ipad==1) length=1;
1291 if (ipad==2) length=1.5;
1292 if (ipad==0){
1293 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1294 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1295 0.,
8076baa0 1296 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1297 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1298 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1299 }else{
1300 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1301 0.,
1302 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
8076baa0 1303 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1304 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1305 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1306 }
1307 }
1308}
1309
1310
1311
1312
1313
0a65832b 1314Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1315 // get Q normalization
1316 // type - 0 Qtot 1 Qmax
1317 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1318 //
f1afff3b 1319 //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);
1320
f1c2a4a3 1321 if (fQNorm==0) return 0;
0a65832b 1322 TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1323 if (!norm) return 0;
f1afff3b 1324 TVectorD &no = *norm;
684602c8 1325 Float_t res =
1326 no[0]+
f1afff3b 1327 no[1]*dr+
1328 no[2]*ty+
1329 no[3]*tz+
1330 no[4]*dr*ty+
1331 no[5]*dr*tz+
1332 no[6]*ty*tz+
1333 no[7]*dr*dr+
1334 no[8]*ty*ty+
1335 no[9]*tz*tz;
1336 res/=no[0];
0a65832b 1337 return res;
1338}
1339
1340
1341
1342void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
1343 //
1344 // set normalization
1345 //
1346 // type - 0 Qtot 1 Qmax
1347 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1348 //
1349
1350 if (fQNorm==0) fQNorm = new TObjArray(6);
1351 fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1352}
236a0d03 1353
b17540e4 1354Float_t AliTPCClusterParam::QnormPos(Int_t ipad,Bool_t isMax, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm, Float_t qt){
1355 //
1356 // Make Q normalization as function of following parameters
1357 // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion
1358 // 1 - dp - relative pad position
1359 // 2 - dt - relative time position
1360 // 3 - di - drift length (norm to 1);
1361 // 4 - dq0 - Tot/Max charge
1362 // 5 - dq1 - Max/Tot charge
1363 // 6 - sy - sigma y - shape
1364 // 7 - sz - sigma z - shape
1365 //
1366
1367 //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
1368 // Following variable used - correspondance to the our variable conventions
1369 //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1370 Double_t dp = ((pad-int(pad)-0.5)*2.);
1371 //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1372 Double_t dt = ((time-int(time)-0.5)*2.);
1373 //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1374 Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1375 //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1376 Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1377 //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1378 Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1379 //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1380 Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1381 //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1382 Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1383 //
1384 //
1385 //
1386 TVectorD * pvec = 0;
1387 if (isMax){
1388 pvec = fPosQMnorm[ipad];
1389 }else{
1390 pvec = fPosQTnorm[ipad];
1391 }
1392 TVectorD &param = *pvec;
1393 //
1394 // Eval part - in correspondance with fit part from debug streamer
1395 //
1396 Double_t result=param[0];
1397 Int_t index =1;
1398 //
1399 result+=dp*param[index++]; //1
1400 result+=dt*param[index++]; //2
1401 result+=dp*dp*param[index++]; //3
1402 result+=dt*dt*param[index++]; //4
1403 result+=dt*dt*dt*param[index++]; //5
1404 result+=dp*dt*param[index++]; //6
1405 result+=dp*dt*dt*param[index++]; //7
1406 result+=(dq0)*param[index++]; //8
1407 result+=(dq1)*param[index++]; //9
1408 //
1409 //
1410 result+=dp*dp*(di)*param[index++]; //10
1411 result+=dt*dt*(di)*param[index++]; //11
1412 result+=dp*dp*sy*param[index++]; //12
1413 result+=dt*sz*param[index++]; //13
1414 result+=dt*dt*sz*param[index++]; //14
1415 result+=dt*dt*dt*sz*param[index++]; //15
1416 //
1417 result+=dp*dp*1*sy*sz*param[index++]; //16
1418 result+=dt*sy*sz*param[index++]; //17
1419 result+=dt*dt*sy*sz*param[index++]; //18
1420 result+=dt*dt*dt*sy*sz*param[index++]; //19
1421 //
1422 result+=dp*dp*(dq0)*param[index++]; //20
1423 result+=dt*1*(dq0)*param[index++]; //21
1424 result+=dt*dt*(dq0)*param[index++]; //22
1425 result+=dt*dt*dt*(dq0)*param[index++]; //23
1426 //
1427 result+=dp*dp*(dq1)*param[index++]; //24
1428 result+=dt*(dq1)*param[index++]; //25
1429 result+=dt*dt*(dq1)*param[index++]; //26
1430 result+=dt*dt*dt*(dq1)*param[index++]; //27
1431
2e5bcb67 1432 if (result<0.75) result=0.75;
1433 if (result>1.25) result=1.25;
1434
b17540e4 1435 return result;
1436
1437}
236a0d03 1438
1439
1440
236a0d03 1441
236a0d03 1442
2e5bcb67 1443Float_t AliTPCClusterParam::PosCorrection(Int_t type, Int_t ipad, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm){
1444
1445 //
1446 // Make postion correction
1447 // type - 0 - y correction
1448 // 1 - z correction
1449 // ipad - 0, 1, 2 - short, medium long pads
1450 // pad - float pad number
1451 // time - float time bin number
1452 // z - z of the cluster
1453 // sy2 - shape of the cluster -pad direction
1454 // sz2 - s - time direction
1455 // qm - maximum charge
1456
1457 //
1458 //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
1459 //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
1460 //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
1461 //chainres->SetAlias("st","(sin(dt)-dt)");
1462 //
1463 //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1464 //chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
1465 //chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1466 //chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1467
1468 //
1469 // Derived variables
1470 //
1471 Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
1472 Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
1473 Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
1474 Double_t st = (TMath::Sin(dt)-dt);
1475 //
1476 Double_t di = TMath::Sqrt(1.-TMath::Abs(z/250.));
1477 Double_t dq = TMath::Sqrt(15./(5.+qm));
1478 Double_t sy = (0.32/TMath::Sqrt(0.01*0.01+sy2));
1479 Double_t sz = (0.32/TMath::Sqrt(0.01*0.01+sz2));
1480 //
1481 //
1482 //
1483 TVectorD * pvec = 0;
1484 if (type==0){
1485 pvec = fPosYcor[ipad];
1486 }else{
1487 pvec = fPosZcor[ipad];
1488 }
1489 TVectorD &param = *pvec;
1490 //
1491 Double_t result=param[0];
1492 Int_t index =1;
1493
1494 if (type==0){
1495 // y corr
1496 result+=(dp)*param[index++]; //1
1497 result+=(dp)*di*param[index++]; //2
1498 result+=(dp)*dq*param[index++]; //3
1499 result+=(dp)*sy*param[index++]; //4
1500 //
1501 result+=(sp)*param[index++]; //5
1502 result+=(sp)*di*param[index++]; //6
1503 result+=(sp)*dq*param[index++]; //7
1504 result+=(sp)*sy*param[index++]; //8
1505 }
1506 if (type==1){
1507 result+=(dt)*param[index++]; //1
1508 result+=(dt)*di*param[index++]; //2
1509 result+=(dt)*dq*param[index++]; //3
1510 result+=(dt)*sz*param[index++]; //4
1511 //
1512 result+=(st)*param[index++]; //5
1513 result+=(st)*di*param[index++]; //6
1514 result+=(st)*dq*param[index++]; //7
1515 result+=(st)*sz*param[index++]; //8
1516 }
1517 return result;
1518}
1519
1520
1521