Adding simple example to load default debug streamer
[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/*
d028aade 64 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
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>
94#include <TProfile2D.h>
0a65832b 95#include <TVectorD.h>
96#include <TObjArray.h>
12ca5da1 97
98ClassImp(AliTPCClusterParam)
99
100
101AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0;
102
103
104/*
105 Example usage fitting parameterization:
106 TFile fres("resol.root"); //tree with resolution and shape
107 TTree * treeRes =(TTree*)fres.Get("Resol");
108
109 AliTPCClusterParam param;
110 param.SetInstance(&param);
111 param.FitResol(treeRes);
112 param.FitRMS(treeRes);
113 TFile fparam("TPCClusterParam.root","recreate");
114 param.Write("Param");
115 //
116 //
117 TFile fparam("TPCClusterParam.root");
118 AliTPCClusterParam *param2 = (AliTPCClusterParam *) fparam.Get("Param");
119 param2->SetInstance(param2);
120 param2->Test(treeRes);
121
122
123 treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
124
125*/
126
127
128
129
130//_ singleton implementation __________________________________________________
131AliTPCClusterParam* AliTPCClusterParam::Instance()
132{
133 //
134 // Singleton implementation
135 // Returns an instance of this class, it is created if neccessary
136 //
137 if (fgInstance == 0){
138 fgInstance = new AliTPCClusterParam();
139 }
140 return fgInstance;
141}
142
143
f1c2a4a3 144AliTPCClusterParam::AliTPCClusterParam():
145 TObject(),
38caa778 146 fRatio(0),
b17540e4 147 fQNorm(0),
148 fQpadTnorm(0), // q pad normalization - Total charge
149 fQpadMnorm(0) // q pad normalization - Max charge
150 //
f1c2a4a3 151{
152 //
153 // Default constructor
154 //
b17540e4 155 fPosQTnorm[0] = 0; fPosQTnorm[1] = 0; fPosQTnorm[2] = 0;
156 fPosQMnorm[0] = 0; fPosQMnorm[1] = 0; fPosQMnorm[2] = 0;
f1c2a4a3 157}
38caa778 158
159AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
160 TObject(param),
161 fRatio(0),
b17540e4 162 fQNorm(0),
163 fQpadTnorm(new TVectorD(*(param.fQpadTnorm))), // q pad normalization - Total charge
164 fQpadMnorm(new TVectorD(*(param.fQpadMnorm))) // q pad normalization - Max charge
165
38caa778 166{
167 //
168 // copy constructor
169 //
170 memcpy(this, &param,sizeof(AliTPCClusterParam));
171 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
b17540e4 172 //
173 if (param.fPosQTnorm[0]){
174 fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
175 fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
176 fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
177 //
178 fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
179 fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
180 fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
181 }
38caa778 182}
183
b17540e4 184
38caa778 185AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
186 //
187 // Assignment operator
188 //
189 if (this != &param) {
190 memcpy(this, &param,sizeof(AliTPCClusterParam));
191 if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
b17540e4 192 if (param.fPosQTnorm[0]){
193 fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
194 fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
195 fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
196 //
197 fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
198 fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
199 fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
200 }
38caa778 201 }
202 return *this;
203}
204
205
f1c2a4a3 206AliTPCClusterParam::~AliTPCClusterParam(){
207 //
208 // destructor
209 //
210 if (fQNorm) fQNorm->Delete();
211 delete fQNorm;
b17540e4 212 if (fPosQTnorm[0]){
213 delete fPosQTnorm[0];
214 delete fPosQTnorm[1];
215 delete fPosQTnorm[2];
216 //
217 delete fPosQMnorm[0];
218 delete fPosQMnorm[1];
219 delete fPosQMnorm[2];
220 }
f1c2a4a3 221}
12ca5da1 222
223
224void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
225 //
226 // Fit z - angular dependence of resolution
227 //
228 // Int_t dim=0, type=0;
229 char varVal[100];
230 sprintf(varVal,"Resol:AngleM:Zm");
231 char varErr[100];
232 sprintf(varErr,"Sigma:AngleS:Zs");
233 char varCut[100];
234 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
235 //
236 Int_t entries = tree->Draw(varVal,varCut);
237 Float_t px[10000], py[10000], pz[10000];
238 Float_t ex[10000], ey[10000], ez[10000];
239 //
240 tree->Draw(varErr,varCut);
241 for (Int_t ipoint=0; ipoint<entries; ipoint++){
242 ex[ipoint]= tree->GetV3()[ipoint];
243 ey[ipoint]= tree->GetV2()[ipoint];
244 ez[ipoint]= tree->GetV1()[ipoint];
245 }
246 tree->Draw(varVal,varCut);
247 for (Int_t ipoint=0; ipoint<entries; ipoint++){
248 px[ipoint]= tree->GetV3()[ipoint];
249 py[ipoint]= tree->GetV2()[ipoint];
250 pz[ipoint]= tree->GetV1()[ipoint];
251 }
252
253 //
254 TLinearFitter fitter(3,"hyp2");
255 for (Int_t ipoint=0; ipoint<entries; ipoint++){
256 Float_t val = pz[ipoint]*pz[ipoint];
257 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
258 Double_t x[2];
259 x[0] = px[ipoint];
260 x[1] = py[ipoint]*py[ipoint];
261 fitter.AddPoint(x,val,err);
262 }
263 fitter.Eval();
264 TVectorD param(3);
265 fitter.GetParameters(param);
266 param0[0] = param[0];
267 param0[1] = param[1];
268 param0[2] = param[2];
269 Float_t chi2 = fitter.GetChisquare()/entries;
270 param0[3] = chi2;
271 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
272 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
273 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
274}
275
276
277void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
278 //
279 // Fit z - angular dependence of resolution
280 //
281 // Int_t dim=0, type=0;
282 char varVal[100];
283 sprintf(varVal,"Resol:AngleM:Zm");
284 char varErr[100];
285 sprintf(varErr,"Sigma:AngleS:Zs");
286 char varCut[100];
287 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
288 //
289 Int_t entries = tree->Draw(varVal,varCut);
290 Float_t px[10000], py[10000], pz[10000];
291 Float_t ex[10000], ey[10000], ez[10000];
292 //
293 tree->Draw(varErr,varCut);
294 for (Int_t ipoint=0; ipoint<entries; ipoint++){
295 ex[ipoint]= tree->GetV3()[ipoint];
296 ey[ipoint]= tree->GetV2()[ipoint];
297 ez[ipoint]= tree->GetV1()[ipoint];
298 }
299 tree->Draw(varVal,varCut);
300 for (Int_t ipoint=0; ipoint<entries; ipoint++){
301 px[ipoint]= tree->GetV3()[ipoint];
302 py[ipoint]= tree->GetV2()[ipoint];
303 pz[ipoint]= tree->GetV1()[ipoint];
304 }
305
306 //
307 TLinearFitter fitter(6,"hyp5");
308 for (Int_t ipoint=0; ipoint<entries; ipoint++){
309 Float_t val = pz[ipoint]*pz[ipoint];
310 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
311 Double_t x[6];
312 x[0] = px[ipoint];
313 x[1] = py[ipoint]*py[ipoint];
314 x[2] = x[0]*x[0];
315 x[3] = x[1]*x[1];
316 x[4] = x[0]*x[1];
317 fitter.AddPoint(x,val,err);
318 }
319 fitter.Eval();
320 TVectorD param(6);
321 fitter.GetParameters(param);
322 param0[0] = param[0];
323 param0[1] = param[1];
324 param0[2] = param[2];
325 param0[3] = param[3];
326 param0[4] = param[4];
327 param0[5] = param[5];
328 Float_t chi2 = fitter.GetChisquare()/entries;
329 param0[6] = chi2;
330 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
331 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
332 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
333 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
334 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
335 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
336}
337
338
339
340
341
342void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
343 //
344 // Fit z - angular dependence of resolution - pad length scaling
345 //
346 // Int_t dim=0, type=0;
347 char varVal[100];
348 sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
349 char varErr[100];
350 sprintf(varErr,"Sigma:AngleS:Zs");
351 char varCut[100];
352 sprintf(varCut,"Dim==%d&&QMean<0",dim);
353 //
354 Int_t entries = tree->Draw(varVal,varCut);
355 Float_t px[10000], py[10000], pz[10000];
356 Float_t ex[10000], ey[10000], ez[10000];
357 //
358 tree->Draw(varErr,varCut);
359 for (Int_t ipoint=0; ipoint<entries; ipoint++){
360 ex[ipoint]= tree->GetV3()[ipoint];
361 ey[ipoint]= tree->GetV2()[ipoint];
362 ez[ipoint]= tree->GetV1()[ipoint];
363 }
364 tree->Draw(varVal,varCut);
365 for (Int_t ipoint=0; ipoint<entries; ipoint++){
366 px[ipoint]= tree->GetV3()[ipoint];
367 py[ipoint]= tree->GetV2()[ipoint];
368 pz[ipoint]= tree->GetV1()[ipoint];
369 }
370
371 //
372 TLinearFitter fitter(3,"hyp2");
373 for (Int_t ipoint=0; ipoint<entries; ipoint++){
374 Float_t val = pz[ipoint]*pz[ipoint];
375 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
376 Double_t x[2];
377 x[0] = px[ipoint];
378 x[1] = py[ipoint]*py[ipoint];
379 fitter.AddPoint(x,val,err);
380 }
381 fitter.Eval();
382 TVectorD param(3);
383 fitter.GetParameters(param);
384 param0[0] = param[0];
385 param0[1] = param[1];
386 param0[2] = param[2];
387 Float_t chi2 = fitter.GetChisquare()/entries;
388 param0[3] = chi2;
389 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
390 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
391 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
392}
393
394void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
395 //
396 // Fit z - angular dependence of resolution - Q scaling
397 //
398 // Int_t dim=0, type=0;
399 char varVal[100];
400 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
401 char varVal0[100];
402 sprintf(varVal0,"Resol:AngleM:Zm");
403 //
404 char varErr[100];
405 sprintf(varErr,"Sigma:AngleS:Zs");
406 char varCut[100];
407 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
408 //
409 Int_t entries = tree->Draw(varVal,varCut);
410 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
411 Float_t ex[20000], ey[20000], ez[20000];
412 //
413 tree->Draw(varErr,varCut);
414 for (Int_t ipoint=0; ipoint<entries; ipoint++){
415 ex[ipoint]= tree->GetV3()[ipoint];
416 ey[ipoint]= tree->GetV2()[ipoint];
417 ez[ipoint]= tree->GetV1()[ipoint];
418 }
419 tree->Draw(varVal,varCut);
420 for (Int_t ipoint=0; ipoint<entries; ipoint++){
421 px[ipoint]= tree->GetV3()[ipoint];
422 py[ipoint]= tree->GetV2()[ipoint];
423 pz[ipoint]= tree->GetV1()[ipoint];
424 }
425 tree->Draw(varVal0,varCut);
426 for (Int_t ipoint=0; ipoint<entries; ipoint++){
427 pu[ipoint]= tree->GetV3()[ipoint];
428 pt[ipoint]= tree->GetV2()[ipoint];
429 }
430
431 //
432 TLinearFitter fitter(5,"hyp4");
433 for (Int_t ipoint=0; ipoint<entries; ipoint++){
434 Float_t val = pz[ipoint]*pz[ipoint];
435 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
436 Double_t x[4];
437 x[0] = pu[ipoint];
438 x[1] = pt[ipoint]*pt[ipoint];
439 x[2] = px[ipoint];
440 x[3] = py[ipoint]*py[ipoint];
441 fitter.AddPoint(x,val,err);
442 }
443
444 fitter.Eval();
445 TVectorD param(5);
446 fitter.GetParameters(param);
447 param0[0] = param[0];
448 param0[1] = param[1];
449 param0[2] = param[2];
450 param0[3] = param[3];
451 param0[4] = param[4];
452 Float_t chi2 = fitter.GetChisquare()/entries;
453 param0[5] = chi2;
454 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
455 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
456 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
457 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
458 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
459}
460
461void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
462 //
463 // Fit z - angular dependence of resolution - Q scaling - parabolic correction
464 //
465 // Int_t dim=0, type=0;
466 char varVal[100];
467 sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
468 char varVal0[100];
469 sprintf(varVal0,"Resol:AngleM:Zm");
470 //
471 char varErr[100];
472 sprintf(varErr,"Sigma:AngleS:Zs");
473 char varCut[100];
474 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
475 //
476 Int_t entries = tree->Draw(varVal,varCut);
477 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
478 Float_t ex[20000], ey[20000], ez[20000];
479 //
480 tree->Draw(varErr,varCut);
481 for (Int_t ipoint=0; ipoint<entries; ipoint++){
482 ex[ipoint]= tree->GetV3()[ipoint];
483 ey[ipoint]= tree->GetV2()[ipoint];
484 ez[ipoint]= tree->GetV1()[ipoint];
485 }
486 tree->Draw(varVal,varCut);
487 for (Int_t ipoint=0; ipoint<entries; ipoint++){
488 px[ipoint]= tree->GetV3()[ipoint];
489 py[ipoint]= tree->GetV2()[ipoint];
490 pz[ipoint]= tree->GetV1()[ipoint];
491 }
492 tree->Draw(varVal0,varCut);
493 for (Int_t ipoint=0; ipoint<entries; ipoint++){
494 pu[ipoint]= tree->GetV3()[ipoint];
495 pt[ipoint]= tree->GetV2()[ipoint];
496 }
497
498 //
499 TLinearFitter fitter(8,"hyp7");
500 for (Int_t ipoint=0; ipoint<entries; ipoint++){
501 Float_t val = pz[ipoint]*pz[ipoint];
502 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
503 Double_t x[7];
504 x[0] = pu[ipoint];
505 x[1] = pt[ipoint]*pt[ipoint];
506 x[2] = x[0]*x[0];
507 x[3] = x[1]*x[1];
508 x[4] = x[0]*x[1];
509 x[5] = px[ipoint];
510 x[6] = py[ipoint]*py[ipoint];
511 //
512 fitter.AddPoint(x,val,err);
513 }
514
515 fitter.Eval();
516 TVectorD param(8);
517 fitter.GetParameters(param);
518 param0[0] = param[0];
519 param0[1] = param[1];
520 param0[2] = param[2];
521 param0[3] = param[3];
522 param0[4] = param[4];
523 param0[5] = param[5];
524 param0[6] = param[6];
525 param0[7] = param[7];
526
527 Float_t chi2 = fitter.GetChisquare()/entries;
528 param0[8] = chi2;
529 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
530 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
531 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
532 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
533 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
534 error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2));
535 error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2));
536 error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2));
537}
538
539
540
541void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
542 //
543 // Fit z - angular dependence of resolution
544 //
545 // Int_t dim=0, type=0;
546 char varVal[100];
547 sprintf(varVal,"RMSm:AngleM:Zm");
548 char varErr[100];
549 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
550 char varCut[100];
551 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
552 //
553 Int_t entries = tree->Draw(varVal,varCut);
554 Float_t px[10000], py[10000], pz[10000];
555 Float_t ex[10000], ey[10000], ez[10000];
556 //
557 tree->Draw(varErr,varCut);
558 for (Int_t ipoint=0; ipoint<entries; ipoint++){
559 ex[ipoint]= tree->GetV3()[ipoint];
560 ey[ipoint]= tree->GetV2()[ipoint];
561 ez[ipoint]= tree->GetV1()[ipoint];
562 }
563 tree->Draw(varVal,varCut);
564 for (Int_t ipoint=0; ipoint<entries; ipoint++){
565 px[ipoint]= tree->GetV3()[ipoint];
566 py[ipoint]= tree->GetV2()[ipoint];
567 pz[ipoint]= tree->GetV1()[ipoint];
568 }
569
570 //
571 TLinearFitter fitter(3,"hyp2");
572 for (Int_t ipoint=0; ipoint<entries; ipoint++){
573 Float_t val = pz[ipoint]*pz[ipoint];
574 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
575 Double_t x[2];
576 x[0] = px[ipoint];
577 x[1] = py[ipoint]*py[ipoint];
578 fitter.AddPoint(x,val,err);
579 }
580 fitter.Eval();
581 TVectorD param(3);
582 fitter.GetParameters(param);
583 param0[0] = param[0];
584 param0[1] = param[1];
585 param0[2] = param[2];
586 Float_t chi2 = fitter.GetChisquare()/entries;
587 param0[3] = chi2;
588 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
589 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
590 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
591}
592
593void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){
594 //
595 // Fit z - angular dependence of resolution - pad length scaling
596 //
597 // Int_t dim=0, type=0;
598 char varVal[100];
599 sprintf(varVal,"RMSm:AngleM*Length:Zm");
600 char varErr[100];
601 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
602 char varCut[100];
603 sprintf(varCut,"Dim==%d&&QMean<0",dim);
604 //
605 Int_t entries = tree->Draw(varVal,varCut);
606 Float_t px[10000], py[10000], pz[10000];
607 Float_t type[10000], ey[10000], ez[10000];
608 //
609 tree->Draw(varErr,varCut);
610 for (Int_t ipoint=0; ipoint<entries; ipoint++){
611 type[ipoint] = tree->GetV3()[ipoint];
612 ey[ipoint] = tree->GetV2()[ipoint];
613 ez[ipoint] = tree->GetV1()[ipoint];
614 }
615 tree->Draw(varVal,varCut);
616 for (Int_t ipoint=0; ipoint<entries; ipoint++){
617 px[ipoint]= tree->GetV3()[ipoint];
618 py[ipoint]= tree->GetV2()[ipoint];
619 pz[ipoint]= tree->GetV1()[ipoint];
620 }
621
622 //
623 TLinearFitter fitter(4,"hyp3");
624 for (Int_t ipoint=0; ipoint<entries; ipoint++){
625 Float_t val = pz[ipoint]*pz[ipoint];
626 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
627 Double_t x[3];
628 x[0] = (type[ipoint]<0.5)? 0.:1.;
629 x[1] = px[ipoint];
630 x[2] = py[ipoint]*py[ipoint];
631 fitter.AddPoint(x,val,err);
632 }
633 fitter.Eval();
634 TVectorD param(4);
635 fitter.GetParameters(param);
636 param0[0] = param[0];
637 param0[1] = param[0]+param[1];
638 param0[2] = param[2];
639 param0[3] = param[3];
640 Float_t chi2 = fitter.GetChisquare()/entries;
641 param0[4] = chi2;
642 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
643 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
644 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
645 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
646}
647
648void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
649 //
650 // Fit z - angular dependence of resolution - Q scaling
651 //
652 // Int_t dim=0, type=0;
653 char varVal[100];
654 sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
655 char varVal0[100];
656 sprintf(varVal0,"RMSm:AngleM:Zm");
657 //
658 char varErr[100];
659 sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
660 char varCut[100];
661 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
662 //
663 Int_t entries = tree->Draw(varVal,varCut);
664 Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
665 Float_t ex[20000], ey[20000], ez[20000];
666 //
667 tree->Draw(varErr,varCut);
668 for (Int_t ipoint=0; ipoint<entries; ipoint++){
669 ex[ipoint]= tree->GetV3()[ipoint];
670 ey[ipoint]= tree->GetV2()[ipoint];
671 ez[ipoint]= tree->GetV1()[ipoint];
672 }
673 tree->Draw(varVal,varCut);
674 for (Int_t ipoint=0; ipoint<entries; ipoint++){
675 px[ipoint]= tree->GetV3()[ipoint];
676 py[ipoint]= tree->GetV2()[ipoint];
677 pz[ipoint]= tree->GetV1()[ipoint];
678 }
679 tree->Draw(varVal0,varCut);
680 for (Int_t ipoint=0; ipoint<entries; ipoint++){
681 pu[ipoint]= tree->GetV3()[ipoint];
682 pt[ipoint]= tree->GetV2()[ipoint];
683 }
684
685 //
686 TLinearFitter fitter(5,"hyp4");
687 for (Int_t ipoint=0; ipoint<entries; ipoint++){
688 Float_t val = pz[ipoint]*pz[ipoint];
689 Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]);
690 Double_t x[4];
691 x[0] = pu[ipoint];
692 x[1] = pt[ipoint]*pt[ipoint];
693 x[2] = px[ipoint];
694 x[3] = py[ipoint]*py[ipoint];
695 fitter.AddPoint(x,val,err);
696 }
697
698 fitter.Eval();
699 TVectorD param(5);
700 fitter.GetParameters(param);
701 param0[0] = param[0];
702 param0[1] = param[1];
703 param0[2] = param[2];
704 param0[3] = param[3];
705 param0[4] = param[4];
706 Float_t chi2 = fitter.GetChisquare()/entries;
707 param0[5] = chi2;
708 error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2));
709 error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2));
710 error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2));
711 error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2));
712 error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2));
713}
714
715
716void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){
717 //
718 // Fit z - angular dependence of resolution - Q scaling
719 //
720 // Int_t dim=0, type=0;
721 char varVal[100];
722 sprintf(varVal,"RMSs:RMSm");
723 //
724 char varCut[100];
725 sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
726 //
727 Int_t entries = tree->Draw(varVal,varCut);
728 Float_t px[20000], py[20000];
729 //
730 tree->Draw(varVal,varCut);
731 for (Int_t ipoint=0; ipoint<entries; ipoint++){
732 px[ipoint]= tree->GetV2()[ipoint];
733 py[ipoint]= tree->GetV1()[ipoint];
734 }
735 TLinearFitter fitter(2,"pol1");
736 for (Int_t ipoint=0; ipoint<entries; ipoint++){
737 Float_t val = py[ipoint];
738 Float_t err = fRatio*px[ipoint];
739 Double_t x[4];
740 x[0] = px[ipoint];
236a0d03 741 if (err>0) fitter.AddPoint(x,val,err);
12ca5da1 742 }
743 fitter.Eval();
744 param0[0]= fitter.GetParameter(0);
745 param0[1]= fitter.GetParameter(1);
746}
747
748
749
750Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){
751 //
752 //
753 //
754 Float_t value=0;
755 value += fParamS0[dim][type][0];
756 value += fParamS0[dim][type][1]*z;
757 value += fParamS0[dim][type][2]*angle*angle;
758 value = TMath::Sqrt(TMath::Abs(value));
759 return value;
760}
761
762
763Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){
764 //
765 //
766 //
767 Float_t value=0;
768 value += fParamS0Par[dim][type][0];
769 value += fParamS0Par[dim][type][1]*z;
770 value += fParamS0Par[dim][type][2]*angle*angle;
771 value += fParamS0Par[dim][type][3]*z*z;
772 value += fParamS0Par[dim][type][4]*angle*angle*angle*angle;
773 value += fParamS0Par[dim][type][5]*z*angle*angle;
774 value = TMath::Sqrt(TMath::Abs(value));
775 return value;
776}
777
778
779
780Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){
781 //
782 //
783 //
784 Float_t value=0;
785 Float_t length=0.75;
786 if (type==1) length=1;
787 if (type==2) length=1.5;
788 value += fParamS1[dim][0];
789 value += fParamS1[dim][1]*z/length;
790 value += fParamS1[dim][2]*angle*angle*length;
791 value = TMath::Sqrt(TMath::Abs(value));
792 return value;
793}
794
795Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
796 //
797 //
798 //
799 Float_t value=0;
800 value += fParamSQ[dim][type][0];
801 value += fParamSQ[dim][type][1]*z;
802 value += fParamSQ[dim][type][2]*angle*angle;
803 value += fParamSQ[dim][type][3]*z/Qmean;
804 value += fParamSQ[dim][type][4]*angle*angle/Qmean;
805 value = TMath::Sqrt(TMath::Abs(value));
806 return value;
807
808
809}
810
811Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
812 //
813 //
814 //
815 Float_t value=0;
816 value += fParamSQPar[dim][type][0];
817 value += fParamSQPar[dim][type][1]*z;
818 value += fParamSQPar[dim][type][2]*angle*angle;
819 value += fParamSQPar[dim][type][3]*z*z;
820 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
821 value += fParamSQPar[dim][type][5]*z*angle*angle;
822 value += fParamSQPar[dim][type][6]*z/Qmean;
823 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
824 value = TMath::Sqrt(TMath::Abs(value));
825 return value;
826
827
828}
829
830Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
831 //
832 //
833 //
834 Float_t value=0;
835 value += fParamSQPar[dim][type][0];
836 value += fParamSQPar[dim][type][1]*z;
837 value += fParamSQPar[dim][type][2]*angle*angle;
838 value += fParamSQPar[dim][type][3]*z*z;
839 value += fParamSQPar[dim][type][4]*angle*angle*angle*angle;
840 value += fParamSQPar[dim][type][5]*z*angle*angle;
841 value += fParamSQPar[dim][type][6]*z/Qmean;
842 value += fParamSQPar[dim][type][7]*angle*angle/Qmean;
843 Float_t valueMean = GetError0Par(dim,type,z,angle);
844 value -= 0.35*0.35*valueMean*valueMean;
845 value = TMath::Sqrt(TMath::Abs(value));
846 return value;
847
848
849}
850
851Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){
852 //
853 // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly
854 //
855 Float_t value=0;
856 value += fParamRMS0[dim][type][0];
857 value += fParamRMS0[dim][type][1]*z;
858 value += fParamRMS0[dim][type][2]*angle*angle;
859 value = TMath::Sqrt(TMath::Abs(value));
860 return value;
861}
862
863Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){
864 //
865 // calculate mean RMS of cluster - z,angle - pad length scalling
866 //
867 Float_t value=0;
868 Float_t length=0.75;
869 if (type==1) length=1;
870 if (type==2) length=1.5;
871 if (type==0){
872 value += fParamRMS1[dim][0];
873 }else{
874 value += fParamRMS1[dim][1];
875 }
876 value += fParamRMS1[dim][2]*z;
877 value += fParamRMS1[dim][3]*angle*angle*length*length;
878 value = TMath::Sqrt(TMath::Abs(value));
879 return value;
880}
881
882Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
883 //
884 // calculate mean RMS of cluster - z,angle, Q dependence
885 //
886 Float_t value=0;
887 value += fParamRMSQ[dim][type][0];
888 value += fParamRMSQ[dim][type][1]*z;
889 value += fParamRMSQ[dim][type][2]*angle*angle;
890 value += fParamRMSQ[dim][type][3]*z/Qmean;
891 value += fParamRMSQ[dim][type][4]*angle*angle/Qmean;
892 value = TMath::Sqrt(TMath::Abs(value));
893 return value;
894}
895
896Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){
897 //
898 // calculates RMS of signal shape fluctuation
899 //
900 Float_t mean = GetRMSQ(dim,type,z,angle,Qmean);
901 Float_t value = fRMSSigmaFit[dim][type][0];
902 value+= fRMSSigmaFit[dim][type][1]*mean;
903 return value;
904}
905
906Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){
907 //
908 // calculates vallue - sigma distortion contribution
909 //
910 Double_t value =0;
911 //
912 Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean);
913 if (rmsL<rmsMeanQ) return value;
914 //
915 Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean);
916 //
917 if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){
918 //1.5 sigma cut on mean
919 value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ;
920 }else{
921 if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){
922 //3 sigma cut on local
923 value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
924 }
925 }
8076baa0 926 return TMath::Sqrt(TMath::Abs(value));
12ca5da1 927}
928
929
930
931void AliTPCClusterParam::FitData(TTree * tree){
932 //
933 // make fits for error param and shape param
934 //
935 FitResol(tree);
936 FitRMS(tree);
937
938}
939
940void AliTPCClusterParam::FitResol(TTree * tree){
941 //
942 SetInstance(this);
943 for (Int_t idir=0;idir<2; idir++){
944 for (Int_t itype=0; itype<3; itype++){
945 Float_t param0[10];
946 Float_t error0[10];
947 // model error param
948 FitResol0(tree, idir, itype,param0,error0);
949 printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
950 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
951 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
952 for (Int_t ipar=0;ipar<4; ipar++){
953 fParamS0[idir][itype][ipar] = param0[ipar];
954 fErrorS0[idir][itype][ipar] = param0[ipar];
955 }
956 // error param with parabolic correction
957 FitResol0Par(tree, idir, itype,param0,error0);
958 printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]);
959 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]);
960 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]);
961 for (Int_t ipar=0;ipar<7; ipar++){
962 fParamS0Par[idir][itype][ipar] = param0[ipar];
963 fErrorS0Par[idir][itype][ipar] = param0[ipar];
964 }
965 //
966 FitResolQ(tree, idir, itype,param0,error0);
967 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
968 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
969 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
970 for (Int_t ipar=0;ipar<6; ipar++){
971 fParamSQ[idir][itype][ipar] = param0[ipar];
972 fErrorSQ[idir][itype][ipar] = param0[ipar];
973 }
974 //
975 FitResolQPar(tree, idir, itype,param0,error0);
976 printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]);
977 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]);
978 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]);
979 for (Int_t ipar=0;ipar<9; ipar++){
980 fParamSQPar[idir][itype][ipar] = param0[ipar];
981 fErrorSQPar[idir][itype][ipar] = param0[ipar];
982 }
983 }
984 }
985 //
986 printf("Resol z-scaled\n");
987 for (Int_t idir=0;idir<2; idir++){
988 Float_t param0[4];
989 Float_t error0[4];
990 FitResol1(tree, idir,param0,error0);
991 printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]);
992 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
993 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
994 for (Int_t ipar=0;ipar<4; ipar++){
995 fParamS1[idir][ipar] = param0[ipar];
996 fErrorS1[idir][ipar] = param0[ipar];
997 }
998 }
999
1000 for (Int_t idir=0;idir<2; idir++){
1001 printf("\nDirection %d\n",idir);
1002 printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]);
1003 for (Int_t itype=0; itype<3; itype++){
1004 Float_t length=0.75;
1005 if (itype==1) length=1;
1006 if (itype==2) length=1.5;
1007 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));
1008 }
1009 }
1010}
1011
1012
1013
1014void AliTPCClusterParam::FitRMS(TTree * tree){
1015 //
1016 SetInstance(this);
1017 for (Int_t idir=0;idir<2; idir++){
1018 for (Int_t itype=0; itype<3; itype++){
1019 Float_t param0[6];
1020 Float_t error0[6];
1021 FitRMS0(tree, idir, itype,param0,error0);
1022 printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]);
1023 printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]);
1024 printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]);
1025 for (Int_t ipar=0;ipar<4; ipar++){
1026 fParamRMS0[idir][itype][ipar] = param0[ipar];
1027 fErrorRMS0[idir][itype][ipar] = param0[ipar];
1028 }
1029 FitRMSQ(tree, idir, itype,param0,error0);
1030 printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]);
1031 printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]);
1032 printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]);
1033 for (Int_t ipar=0;ipar<6; ipar++){
1034 fParamRMSQ[idir][itype][ipar] = param0[ipar];
1035 fErrorRMSQ[idir][itype][ipar] = param0[ipar];
1036 }
1037 }
1038 }
1039 //
1040 printf("RMS z-scaled\n");
1041 for (Int_t idir=0;idir<2; idir++){
1042 Float_t param0[5];
1043 Float_t error0[5];
1044 FitRMS1(tree, idir,param0,error0);
1045 printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]);
1046 printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]);
1047 printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]);
1048 for (Int_t ipar=0;ipar<5; ipar++){
1049 fParamRMS1[idir][ipar] = param0[ipar];
1050 fErrorRMS1[idir][ipar] = param0[ipar];
1051 }
1052 }
1053
1054 for (Int_t idir=0;idir<2; idir++){
1055 printf("\nDirection %d\n",idir);
1056 printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]);
1057 for (Int_t itype=0; itype<3; itype++){
1058 Float_t length=0.75;
1059 if (itype==1) length=1;
1060 if (itype==2) length=1.5;
1061 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);
1062 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);
1063 }
1064 }
1065 //
1066 // Fit RMS sigma
1067 //
1068 printf("RMS fluctuation parameterization \n");
1069 for (Int_t idir=0;idir<2; idir++){
1070 for (Int_t itype=0; itype<3; itype++){
1071 Float_t param0[5];
1072 Float_t error0[5];
1073 FitRMSSigma(tree, idir,itype,param0,error0);
1074 printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]);
1075 for (Int_t ipar=0;ipar<2; ipar++){
1076 fRMSSigmaFit[idir][itype][ipar] = param0[ipar];
1077 }
1078 }
1079 }
1080 //
1081 // store systematic error end RMS fluctuation parameterization
1082 //
1083 TH1F hratio("hratio","hratio",100,-0.1,0.1);
1084 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0");
1085 fErrorRMSSys[0] = hratio.GetRMS();
1086 tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0");
1087 fErrorRMSSys[1] = hratio.GetRMS();
1088 TH1F hratioR("hratioR","hratioR",100,0,0.2);
1089 tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0");
1090 fRMSSigmaRatio[0][0]=hratioR.GetMean();
1091 fRMSSigmaRatio[0][1]=hratioR.GetRMS();
1092 tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0");
1093 fRMSSigmaRatio[1][0]=hratioR.GetMean();
1094 fRMSSigmaRatio[1][1]=hratioR.GetRMS();
1095}
1096
1097void AliTPCClusterParam::Test(TTree * tree, const char *output){
1098 //
1099 // Draw standard quality histograms to output file
1100 //
1101 SetInstance(this);
1102 TFile f(output,"recreate");
1103 f.cd();
1104 //
1105 // 1D histograms - resolution
1106 //
1107 for (Int_t idim=0; idim<2; idim++){
1108 for (Int_t ipad=0; ipad<3; ipad++){
1109 char hname1[300];
1110 char hcut1[300];
1111 char hexp1[300];
1112 //
1113 sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
1114 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1115 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1116 TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2);
1117 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1118 tree->Draw(hexp1,hcut1,"");
1119 his1DRel0.Write();
1120 //
1121 sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
1122 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1123 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
1124 TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
1125 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1126 tree->Draw(hexp1,hcut1,"");
1127 his1DRel0Par.Write();
1128 //
1129 }
1130 }
1131 //
1132 // 2D histograms - resolution
1133 //
1134 for (Int_t idim=0; idim<2; idim++){
1135 for (Int_t ipad=0; ipad<3; ipad++){
1136 char hname1[300];
1137 char hcut1[300];
1138 char hexp1[300];
1139 //
1140 sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
1141 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1142 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1143 TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1);
1144 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1145 tree->Draw(hexp1,hcut1,"");
1146 profDRel0.Write();
1147 //
1148 sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
1149 sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
1150 sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
1151 TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
1152 sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
1153 tree->Draw(hexp1,hcut1,"");
1154 profDRel0Par.Write();
1155 //
1156 }
1157 }
1158}
1159
1160
1161
1162void AliTPCClusterParam::Print(Option_t* /*option*/) const{
1163 //
1164 // Print param Information
1165 //
1166
1167 //
1168 // Error parameterization
1169 //
1170 printf("\nResolution Scaled factors\n");
1171 printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1172 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])),
1173 TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
12ca5da1 1174 for (Int_t ipad=0; ipad<3; ipad++){
1175 Float_t length=0.75;
1176 if (ipad==1) length=1;
1177 if (ipad==2) length=1.5;
1178 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1179 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
8076baa0 1180 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
1181 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
1182 TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
12ca5da1 1183 }
1184 for (Int_t ipad=0; ipad<3; ipad++){
1185 Float_t length=0.75;
1186 if (ipad==1) length=1;
1187 if (ipad==2) length=1.5;
1188 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
1189 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
8076baa0 1190 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
1191 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
1192 TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
12ca5da1 1193 }
1194 printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
1195 TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
1196
1197 for (Int_t ipad=0; ipad<3; ipad++){
1198 Float_t length=0.75;
1199 if (ipad==1) length=1;
1200 if (ipad==2) length=1.5;
1201 printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
1202 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
8076baa0 1203 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
1204 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
1205 TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
12ca5da1 1206 }
1207 for (Int_t ipad=0; ipad<3; ipad++){
1208 Float_t length=0.75;
1209 if (ipad==1) length=1;
1210 if (ipad==2) length=1.5;
1211 printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
8076baa0 1212 TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
1213 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
1214 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
1215 TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
12ca5da1 1216 }
1217
1218 //
1219 // RMS scaling
1220 //
1221 printf("\n");
1222 printf("\nRMS Scaled factors\n");
1223 printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
8076baa0 1224 printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
1225 TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
1226 TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
1227 TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
1228 TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
1229 TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
12ca5da1 1230 for (Int_t ipad=0; ipad<3; ipad++){
1231 Float_t length=0.75;
1232 if (ipad==1) length=1;
1233 if (ipad==2) length=1.5;
1234 if (ipad==0){
1235 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1236 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
1237 0.,
8076baa0 1238 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1239 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1240 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1241 }else{
1242 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1243 0.,
1244 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
8076baa0 1245 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
1246 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
1247 TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
12ca5da1 1248 }
1249 }
1250 printf("\n");
8076baa0 1251 printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
1252 TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
1253 TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
1254 TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
1255 TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
1256 TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
12ca5da1 1257 for (Int_t ipad=0; ipad<3; ipad++){
1258 Float_t length=0.75;
1259 if (ipad==1) length=1;
1260 if (ipad==2) length=1.5;
1261 if (ipad==0){
1262 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1263 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
1264 0.,
8076baa0 1265 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1266 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1267 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1268 }else{
1269 printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
1270 0.,
1271 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
8076baa0 1272 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
1273 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
1274 TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
12ca5da1 1275 }
1276 }
1277}
1278
1279
1280
1281
1282
0a65832b 1283Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){
1284 // get Q normalization
1285 // type - 0 Qtot 1 Qmax
1286 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1287 //
f1afff3b 1288 //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
1289
f1c2a4a3 1290 if (fQNorm==0) return 0;
0a65832b 1291 TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
1292 if (!norm) return 0;
f1afff3b 1293 TVectorD &no = *norm;
684602c8 1294 Float_t res =
1295 no[0]+
f1afff3b 1296 no[1]*dr+
1297 no[2]*ty+
1298 no[3]*tz+
1299 no[4]*dr*ty+
1300 no[5]*dr*tz+
1301 no[6]*ty*tz+
1302 no[7]*dr*dr+
1303 no[8]*ty*ty+
1304 no[9]*tz*tz;
1305 res/=no[0];
0a65832b 1306 return res;
1307}
1308
1309
1310
1311void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){
1312 //
1313 // set normalization
1314 //
1315 // type - 0 Qtot 1 Qmax
1316 // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
1317 //
1318
1319 if (fQNorm==0) fQNorm = new TObjArray(6);
1320 fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
1321}
236a0d03 1322
b17540e4 1323Float_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){
1324 //
1325 // Make Q normalization as function of following parameters
1326 // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion
1327 // 1 - dp - relative pad position
1328 // 2 - dt - relative time position
1329 // 3 - di - drift length (norm to 1);
1330 // 4 - dq0 - Tot/Max charge
1331 // 5 - dq1 - Max/Tot charge
1332 // 6 - sy - sigma y - shape
1333 // 7 - sz - sigma z - shape
1334 //
1335
1336 //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
1337 // Following variable used - correspondance to the our variable conventions
1338 //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
1339 Double_t dp = ((pad-int(pad)-0.5)*2.);
1340 //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
1341 Double_t dt = ((time-int(time)-0.5)*2.);
1342 //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
1343 Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
1344 //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
1345 Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
1346 //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
1347 Double_t dq1 = 5.*(qm+2.)/(qt+2.);
1348 //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
1349 Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2);
1350 //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
1351 Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2);
1352 //
1353 //
1354 //
1355 TVectorD * pvec = 0;
1356 if (isMax){
1357 pvec = fPosQMnorm[ipad];
1358 }else{
1359 pvec = fPosQTnorm[ipad];
1360 }
1361 TVectorD &param = *pvec;
1362 //
1363 // Eval part - in correspondance with fit part from debug streamer
1364 //
1365 Double_t result=param[0];
1366 Int_t index =1;
1367 //
1368 result+=dp*param[index++]; //1
1369 result+=dt*param[index++]; //2
1370 result+=dp*dp*param[index++]; //3
1371 result+=dt*dt*param[index++]; //4
1372 result+=dt*dt*dt*param[index++]; //5
1373 result+=dp*dt*param[index++]; //6
1374 result+=dp*dt*dt*param[index++]; //7
1375 result+=(dq0)*param[index++]; //8
1376 result+=(dq1)*param[index++]; //9
1377 //
1378 //
1379 result+=dp*dp*(di)*param[index++]; //10
1380 result+=dt*dt*(di)*param[index++]; //11
1381 result+=dp*dp*sy*param[index++]; //12
1382 result+=dt*sz*param[index++]; //13
1383 result+=dt*dt*sz*param[index++]; //14
1384 result+=dt*dt*dt*sz*param[index++]; //15
1385 //
1386 result+=dp*dp*1*sy*sz*param[index++]; //16
1387 result+=dt*sy*sz*param[index++]; //17
1388 result+=dt*dt*sy*sz*param[index++]; //18
1389 result+=dt*dt*dt*sy*sz*param[index++]; //19
1390 //
1391 result+=dp*dp*(dq0)*param[index++]; //20
1392 result+=dt*1*(dq0)*param[index++]; //21
1393 result+=dt*dt*(dq0)*param[index++]; //22
1394 result+=dt*dt*dt*(dq0)*param[index++]; //23
1395 //
1396 result+=dp*dp*(dq1)*param[index++]; //24
1397 result+=dt*(dq1)*param[index++]; //25
1398 result+=dt*dt*(dq1)*param[index++]; //26
1399 result+=dt*dt*dt*(dq1)*param[index++]; //27
1400
1401 return result;
1402
1403}
236a0d03 1404
1405
1406
236a0d03 1407
236a0d03 1408