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