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