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