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