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