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