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