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