]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/qaRec/AliTRDclusterResolution.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add processing of position/error for y center pad dependence
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDclusterResolution.cxx
CommitLineData
9462866a 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-commercialf 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/* $Id: AliTRDclusterResolution.cxx */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// TRD cluster error parameterization //
21// //
22// This class is designed to produce the reference plots for a detailed study//
23// and parameterization of TRD cluster errors. The following effects are taken//
24// into account : //
25// - dependence with the total charge of the cluster //
26// - dependence with the distance from the center pad. This is monitored
27// for each layer individually since the pad size varies with layer
28// - dependence with the drift length - here the influence of anisochronity
29// and diffusion are searched
30// - dependence with the distance to the anode wire - anisochronity effects
31// - dependence with track angle (for y resolution)
32// The correlation between effects is taken into account.
33//
34// Since magnetic field plays a very important role in the TRD measurement
35// the ExB correction is forced by the setter function SetExB(Int_t). The
36// argument is the detector index, if none is specified all will be
37// considered.
38//
39// Two cases are of big importance.
40// - comparison with MC
41// - comparison with Kalman fit. In this case the covariance matrix of the
42// Kalman fit are needed.
43//
44// The functionalities implemented in this class are based on the storage
45// class AliTRDclusterInfo.
46//
47// The Method
48// ----------
49//
50// The method to disentangle s_y and s_x is based on the relation (see also fig.)
51// BEGIN_LATEX
52// #sigma^{2} = #sigma^{2}_{y} + tg^{2}(#alpha_{L})*#sigma^{2}_{x_{d}} + tg^{2}(#phi-#alpha_{L})*(#sigma^{2}_{x_{d}}+#sigma^{2}_{x_{c}})
53// END_LATEX
54// with
55// BEGIN_LATEX
56// #sigma^{2}_{x_{c}} #approx 0
57// END_LATEX
58// we suppose the chamber is well calibrated for t_{0} and aligned in
59// radial direction.
60//
fa7831fb 61// Clusters can be radially shifted due to three causes:
9462866a 62// - globally shifted - due to residual misalignment/miscalibration(t0)
63// - locally shifted - due to different local drift velocity from the mean
64// - randomly shifted - due to neighboring (radial direction) clusters
65// charge induced by asymmetry of the TRF.
66//
67// We estimate this effects by the relations:
68// BEGIN_LATEX
69// #mu_{y} = tg(#alpha_{L})*#Delta x_{d}(...) + tg(#phi-#alpha_{L})*(#Delta x_{c}(...) + #Delta x_{d}(...))
70// END_LATEX
71// where
72// BEGIN_LATEX
73// #Delta x_{d}(...) = (<v_{d}> + #delta v_{d}(x_{d}, d)) * (t + t^{*}(Q))
74// END_LATEX
75// and we specified explicitely the variation of drift velocity parallel
76// with the track (x_{d}) and perpendicular to it due to anisochronity (d).
77//
78// For estimating the contribution from asymmetry of TRF the following
79// parameterization is being used
80// BEGIN_LATEX
81// t^{*}(Q) = #delta_{0} * #frac{Q_{t+1} - Q_{t-1}}{Q_{t-1} + Q_{t} + Q_{t+1}}
82// END_LATEX
83//
84//
fa7831fb 85// Clusters can also be r-phi shifted due to:
86// - wrong PRF or wrong cuts at digits level
87//The following correction is applied :
88// BEGIN_LATEX
89// <#Delta y> = a + b * sin(c*y_{pw})
90// END_LATEX
91
9462866a 92// The Models
93//
94// Parameterization against total charge
95//
96// Obtained for B=0T at phi=0. All other effects integrated out.
97// BEGIN_LATEX
98// #sigma^{2}_{y}(Q) = #sigma^{2}_{y}(...) + b(#frac{1}{Q} - #frac{1}{Q_{0}})
99// END_LATEX
100// For B diff 0T the error of the average ExB correction error has to be subtracted !!
101//
102// Parameterization Sx
103//
104// The parameterization of the error in the x direction can be written as
105// BEGIN_LATEX
106// #sigma_{x} = #sigma_{x}^{||} + #sigma_{x}^{#perp}
107// END_LATEX
108//
109// where the parallel component is given mainly by the TRF width while
110// the perpendicular component by the anisochronity. The model employed for
111// the parallel is gaus(0)+expo(3) with the following parameters
112// 1 C 5.49018e-01 1.23854e+00 3.84540e-04 -8.21084e-06
113// 2 M 7.82999e-01 6.22531e-01 2.71272e-04 -6.88485e-05
114// 3 S 2.74451e-01 1.13815e+00 2.90667e-04 1.13493e-05
115// 4 E1 2.53596e-01 1.08646e+00 9.95591e-05 -2.11625e-05
116// 5 E2 -2.40078e-02 4.26520e-01 4.67153e-05 -2.35392e-04
117//
118// and perpendicular to the track is pol2 with the parameters
119//
120// Par_0 = 0.190676 +/- 0.41785
121// Par_1 = -3.9269 +/- 7.49862
122// Par_2 = 14.7851 +/- 27.8012
123//
124// Parameterization Sy
125//
126// The parameterization of the error in the y direction along track uses
127// BEGIN_LATEX
128// #sigma_{y}^{||} = #sigma_{y}^{0} + exp(-a*(x-b))
129// END_LATEX
130//
131// with following values for the parameters:
132// 1 sy0 2.60967e-01 2.99652e-03 7.82902e-06 -1.89636e-04
133// 2 a -7.68941e+00 1.87883e+00 3.84539e-04 9.38268e-07
134// 3 b -3.41160e-01 7.72850e-02 1.63231e-05 2.51602e-05
135//
136//==========================================================================
137// Example how to retrive reference plots from the task
138// void steerClErrParam(Int_t fig=0)
139// {
140// gSystem->Load("libANALYSIS.so");
141// gSystem->Load("libTRDqaRec.so");
142//
143// // initialize DB manager
144// AliCDBManager *cdb = AliCDBManager::Instance();
145// cdb->SetDefaultStorage("local://$ALICE_ROOT");
146// cdb->SetRun(0);
147// // initialize magnetic field.
148// AliMagFMaps *field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
149// AliTracker::SetFieldMap(field, kTRUE);
150//
151// AliTRDclusterResolution *res = new AliTRDclusterResolution();
152// res->SetMCdata();
153// res->Load("TRD.TaskClErrParam.root");
154// res->SetExB();
155// if(!res->PostProcess()) return;
156// new TCanvas;
157// res->GetRefFigure(fig);
158// }
159//
160// Authors: //
161// Alexandru Bercuci <A.Bercuci@gsi.de> //
162////////////////////////////////////////////////////////////////////////////
163
b2dc316d 164#include "AliTRDclusterResolution.h"
165#include "AliTRDtrackInfo/AliTRDclusterInfo.h"
5198d8c6 166#include "AliTRDgeometry.h"
6bc4a8f4 167#include "AliTRDcalibDB.h"
168#include "Cal/AliTRDCalROC.h"
169#include "Cal/AliTRDCalDet.h"
b2dc316d 170
171#include "AliLog.h"
6bc4a8f4 172#include "AliTracker.h"
173#include "AliCDBManager.h"
b2dc316d 174
175#include "TObjArray.h"
176#include "TAxis.h"
fc0946a7 177#include "TF1.h"
178#include "TGraphErrors.h"
b2dc316d 179#include "TH2I.h"
180#include "TMath.h"
9462866a 181#include "TLinearFitter.h"
182
183#include "TCanvas.h"
184#include "TSystem.h"
185
b2dc316d 186
187ClassImp(AliTRDclusterResolution)
188
189
190//_______________________________________________________
191AliTRDclusterResolution::AliTRDclusterResolution()
5198d8c6 192 : AliTRDrecoTask("ClErrParam", "Cluster Error Parametrization")
b2dc316d 193 ,fInfo(0x0)
fc0946a7 194 ,fResults(0x0)
195 ,fAt(0x0)
196 ,fAd(0x0)
197 ,fExB(0.)
9462866a 198 ,fVdrift(0.)
6bc4a8f4 199 ,fDet(-1)
b2dc316d 200{
9462866a 201 // ideal equidistant binning
202 //fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15);
203
204 // equidistant binning for scaled for X0-MC (track ref)
205 fAt = new TAxis(kNTB, -0.075+.088, (kNTB-.5)*.15+.088);
206
207 // non-equidistant binning after vdrift correction
208// const Double_t x[kNTB+1] = {
209// 0.000000, 0.185084, 0.400490, 0.519701, 0.554653, 0.653150, 0.805063, 0.990261,
210// 1.179610, 1.356406, 1.524094, 1.685499, 1.843083, 1.997338, 2.148077, 2.298274,
211// 2.448656, 2.598639, 2.747809, 2.896596, 3.045380, 3.195000, 3.340303, 3.461688,
212// 3.540094, 3.702677};
213// fAt = new TAxis(kNTB, x);
214
fc0946a7 215 fAd = new TAxis(kND, 0., .25);
b2dc316d 216}
217
218//_______________________________________________________
219AliTRDclusterResolution::~AliTRDclusterResolution()
220{
fc0946a7 221 if(fAd) delete fAd;
222 if(fAt) delete fAt;
223 if(fResults){
224 fResults->Delete();
225 delete fResults;
226 }
b2dc316d 227}
228
229//_______________________________________________________
230void AliTRDclusterResolution::ConnectInputData(Option_t *)
231{
232 fInfo = dynamic_cast<TObjArray *>(GetInputData(0));
233}
234
235//_______________________________________________________
236void AliTRDclusterResolution::CreateOutputObjects()
237{
238 OpenFile(0, "RECREATE");
239 fContainer = Histos();
240}
241
242//_______________________________________________________
e15179be 243Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
b2dc316d 244{
e15179be 245 if(!fResults) return kFALSE;
fc0946a7 246
247 TObjArray *arr = 0x0;
248 TH2 *h2 = 0x0;
9462866a 249 TGraphErrors *gm(0x0), *gs(0x0), *gp(0x0);
fc0946a7 250 switch(ifig){
251 case kQRes:
252 if(!(arr = (TObjArray*)fResults->At(kQRes))) break;
253 if(!(gm = (TGraphErrors*)arr->At(0))) break;
254 if(!(gs = (TGraphErrors*)arr->At(1))) break;
9462866a 255 if(!(gp = (TGraphErrors*)arr->At(2))) break;
fc0946a7 256 gs->Draw("apl");
9462866a 257 gs->GetHistogram()->GetYaxis()->SetRangeUser(-.01, .6);
258 gs->GetHistogram()->SetXTitle("Q [a.u.]");
259 gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [mm] / freq");
260 gm->Draw("pl");
261 gp->Draw("pl");
e15179be 262 return kTRUE;
9462866a 263 case kCenter:
264 if(!(arr = (TObjArray*)fResults->At(kCenter))) break;
265 if(!(h2 = (TH2F*)arr->At(0))) break;
fc0946a7 266 h2->Draw("lego2");
e15179be 267 return kTRUE;
9462866a 268 case kSigm:
269 if(!(arr = (TObjArray*)fResults->At(kSigm))) break;
270 gPad->Divide(2, 1);
271 gPad->cd(1);
272 if(!(h2 = (TH2F*)arr->At(0))) break;
273 h2->Draw("lego2fb");
274 gPad->cd(2);
275 if(!(h2 = (TH2F*)arr->At(1))) break;
276 h2->Draw("lego2fb");
277 return kTRUE;
278 case kMean:
279 if(!(arr = (TObjArray*)fResults->At(kMean))) break;
280 arr->ls();
281/* gPad->Divide(2, 1);
282 gPad->cd(1);*/
283 if(!(h2 = (TH2F*)arr->At(0))) break;
284 h2->Draw("lego2 fb");
285/* gPad->cd(2);
286 if(!(h2 = (TH2F*)arr->At(1))) break;
287 h2->Draw("lego2fb");*/
e15179be 288 return kTRUE;
fc0946a7 289 default:
290 break;
291 }
9462866a 292 AliWarning("No container/data found.");
e15179be 293 return kFALSE;
b2dc316d 294}
295
296//_______________________________________________________
297TObjArray* AliTRDclusterResolution::Histos()
298{
299 if(fContainer) return fContainer;
9462866a 300 fContainer = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainers));
b2dc316d 301 //fContainer->SetOwner(kTRUE);
302
5198d8c6 303 TH2I *h2 = 0x0;
304 TObjArray *arr = 0x0;
305
306 fContainer->AddAt(h2 = new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), kQRes);
307 h2->SetXTitle("log(q) [a.u.]");
308 h2->SetYTitle("#Delta y[cm]");
309 h2->SetZTitle("entries");
310
9462866a 311 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kCenter);
5198d8c6 312 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
6bc4a8f4 313 arr->AddAt(h2 = new TH2I(Form("h_y%d", ily), Form("Ly[%d]", ily), 51, -.51, .51, 100, -.5, .5), ily);
5198d8c6 314 h2->SetXTitle("y_{w} [w]");
315 h2->SetYTitle("#Delta y[cm]");
316 h2->SetZTitle("entries");
317 }
318
9462866a 319 fContainer->AddAt(arr = new TObjArray(kN), kSigm);
b2dc316d 320 Int_t ih = 0;
fc0946a7 321 for(Int_t id=1; id<=fAd->GetNbins(); id++){
322 for(Int_t it=1; it<=fAt->GetNbins(); it++){
9462866a 323 arr->AddAt(h2 = new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
5198d8c6 324 h2->SetXTitle("tg#phi");
325 h2->SetYTitle("#Delta y[cm]");
326 h2->SetZTitle("entries");
b2dc316d 327 }
328 }
9462866a 329
330 fContainer->AddAt(arr = new TObjArray(kN), kMean);
331 ih = 0;
332 for(Int_t id=1; id<=fAd->GetNbins(); id++){
333 for(Int_t it=1; it<=fAt->GetNbins(); it++){
334 arr->AddAt(h2 = new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
335 h2->SetXTitle("tg#phi-h*tg(#theta)");
336 h2->SetYTitle("#Delta y[cm]");
337 h2->SetZTitle("entries");
338 }
339 }
340
b2dc316d 341 return fContainer;
342}
343
344//_______________________________________________________
345void AliTRDclusterResolution::Exec(Option_t *)
346{
6bc4a8f4 347 if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the task.");
348
87b166d3 349 Int_t det, t;
fc0946a7 350 Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
b2dc316d 351 TH2I *h2 = 0x0;
5198d8c6 352
9462866a 353 // define limits around ExB for which x contribution is negligible
354 const Float_t kDtgPhi = 3.5e-2; //(+- 2 deg)
355
356 TObjArray *arr0 = (TObjArray*)fContainer->At(kCenter);
357 TObjArray *arr1 = (TObjArray*)fContainer->At(kSigm);
358 TObjArray *arr2 = (TObjArray*)fContainer->At(kMean);
5198d8c6 359
b2dc316d 360 const AliTRDclusterInfo *cli = 0x0;
361 TIterator *iter=fInfo->MakeIterator();
362 while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
6bc4a8f4 363 cli->GetCluster(det, x, y, z, q, t, covcl);
364 if(fDet>=0 && fDet!=det) continue;
9462866a 365
366 dy = cli->GetResolution();
367 cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
368
369 // resolution as a function of cluster charge
370 // only for phi equal exB
371 if(TMath::Abs(dydx-fExB) < kDtgPhi){
372 h2 = (TH2I*)fContainer->At(kQRes);
373 h2->Fill(TMath::Log(q), dy);
374 }
375
376 // do not use problematic clusters in resolution analysis
377 // TODO define limits as calibration aware (gain) !!
378 if(q<20. || q>250.) continue;
379
380 // resolution as a function of y displacement from pad center
381 // only for phi equal exB
382 if(TMath::Abs(dydx-fExB) < kDtgPhi){
383 h2 = (TH2I*)arr0->At(AliTRDgeometry::GetLayer(det));
384 h2->Fill(cli->GetYDisplacement(), dy);
385 }
6bc4a8f4 386
fc0946a7 387 Int_t it = fAt->FindBin(cli->GetDriftLength());
388 if(it==0 || it == fAt->GetNbins()+1){
b2dc316d 389 AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength()));
390 continue;
391 }
fc0946a7 392 Int_t id = fAd->FindBin(cli->GetAnisochronity());
393 if(id==0 || id == fAd->GetNbins()+1){
b2dc316d 394 AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity()));
395 continue;
396 }
9462866a 397 // calculate index of cluster position in array
398 Int_t hid = (id-1)*kNTB+it-1;
b2dc316d 399
9462866a 400 // fill histo for resolution (sigma)
401 h2 = (TH2I*)arr1->At(hid);
6bc4a8f4 402 h2->Fill(dydx, dy);
b2dc316d 403
9462866a 404 // fill histo for systematic (mean)
405 h2 = (TH2I*)arr2->At(hid);
406 h2->Fill(dydx-cli->GetTilt()*dzdx, dy);
b2dc316d 407 }
408 PostData(0, fContainer);
409}
410
411
412//_______________________________________________________
413Bool_t AliTRDclusterResolution::PostProcess()
414{
fc0946a7 415 if(!fContainer) return kFALSE;
6bc4a8f4 416 if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the post processing.");
fc0946a7 417
5198d8c6 418 TObjArray *arr = 0x0;
fc0946a7 419 TH2 *h2 = 0x0;
420 if(!fResults){
421 TGraphErrors *g = 0x0;
9462866a 422 fResults = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainers));
fc0946a7 423 fResults->SetOwner();
9462866a 424 fResults->AddAt(arr = new TObjArray(3), kQRes);
fc0946a7 425 arr->SetOwner();
426 arr->AddAt(g = new TGraphErrors(), 0);
427 g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
428 g->SetMarkerStyle(7);
429 arr->AddAt(g = new TGraphErrors(), 1);
430 g->SetLineColor(kRed); g->SetMarkerColor(kRed);
431 g->SetMarkerStyle(23);
9462866a 432 arr->AddAt(g = new TGraphErrors(), 2);
433 g->SetLineColor(kGreen); g->SetMarkerColor(kGreen);
fc0946a7 434 g->SetMarkerStyle(7);
fc0946a7 435
fa7831fb 436 fResults->AddAt(arr = new TObjArray(3), kCenter);
9462866a 437 arr->SetOwner();
438 arr->AddAt(h2 = new TH2F("hYM", "",
439 AliTRDgeometry::kNlayer, -.5, AliTRDgeometry::kNlayer-.5, 51, -.51, .51), 0);
440 h2->SetXTitle("ly");
441 h2->SetYTitle("y [w]");
442 h2->SetZTitle("#mu_{x} [mm]");
443 arr->AddAt(h2 = (TH2F*)h2->Clone("hYS"), 1);
444 h2->SetZTitle("#sigma_{x} [mm]");
fa7831fb 445 arr->AddAt(h2 = (TH2F*)h2->Clone("hYP"), 2);
446 h2->SetZTitle("entries");
9462866a 447
448 fResults->AddAt(arr = new TObjArray(2), kSigm);
449 arr->SetOwner();
450 arr->AddAt(h2 = new TH2F("hSX", "",
fc0946a7 451 fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(),
9462866a 452 fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()), 0);
fc0946a7 453 h2->SetXTitle("d [mm]");
454 h2->SetYTitle("x [mm]");
9462866a 455 h2->SetZTitle("#sigma_{x} [mm]");
456 arr->AddAt(h2 = (TH2F*)h2->Clone("hSY"), 1);
457 h2->SetZTitle("#sigma_{y} [mm]");
458
459 fResults->AddAt(arr = new TObjArray(2), kMean);
460 arr->SetOwner();
461 arr->AddAt(h2 = (TH2F*)h2->Clone("hDX"), 0);
462 h2->SetZTitle("dx [mm]");
463 arr->AddAt(h2 = (TH2F*)h2->Clone("hX0"), 1);
464 h2->SetZTitle("x0 [mm]");
fc0946a7 465 } else {
466 TObject *o = 0x0;
467 TIterator *iter=fResults->MakeIterator();
5198d8c6 468 while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point
fc0946a7 469 }
9462866a 470
471 // precalculated value of tg^2(alpha_L)
472 Double_t exb2 = fExB*fExB;
473 // square of the mean value of sigma drift length.
474 // has to come from previous calibration
475 //Double_t sxd2 = 1.;// [mm^2]
fc0946a7 476
9462866a 477 printf("ExB[%e] ExB2[%e]\n", fExB, exb2);
fc0946a7 478
479 // process resolution dependency on charge
fa7831fb 480 ProcessCharge();
fc0946a7 481
fc0946a7 482 // process resolution dependency on y displacement
fa7831fb 483 ProcessCenterPad();
fc0946a7 484
5198d8c6 485 // process resolution dependency on drift legth and drift cell width
9462866a 486 ProcessSigma();
487
488 // process systematic shift on drift legth and drift cell width
fa7831fb 489 ProcessMean();
fc0946a7 490
491 return kTRUE;
b2dc316d 492}
493
6bc4a8f4 494//_______________________________________________________
9462866a 495Bool_t AliTRDclusterResolution::SetExB(Int_t det, Int_t col, Int_t row)
6bc4a8f4 496{
497 // check OCDB
498 AliCDBManager *cdb = AliCDBManager::Instance();
499 if(cdb->GetRun() < 0){
500 AliError("OCDB manager not properly initialized");
501 return kFALSE;
502 }
503
504 // check magnetic field
505 if(TMath::Abs(AliTracker::GetBz()) < 1.e-10){
506 AliWarning("B=0. Magnetic field may not be initialized. Continue if you know what you are doing !");
507 }
508
509 // set reference detector if any
510 if(det>=0 && det<AliTRDgeometry::kNdet) fDet = det;
511 else det = 0;
512
513 AliTRDcalibDB *fCalibration = AliTRDcalibDB::Instance();
514 AliTRDCalROC *fCalVdriftROC = fCalibration->GetVdriftROC(det);
515 const AliTRDCalDet *fCalVdriftDet = fCalibration->GetVdriftDet();
516
9462866a 517 fVdrift = fCalVdriftDet->GetValue(det) * fCalVdriftROC->GetValue(col, row);
518 fExB = fCalibration->GetOmegaTau(fVdrift, -0.1*AliTracker::GetBz());
6bc4a8f4 519 SetBit(kExB);
520 return kTRUE;
521}
b2dc316d 522
9462866a 523//_______________________________________________________
524void AliTRDclusterResolution::ProcessCharge()
525{
526 TH2I *h2 = 0x0;
527 if((h2 = (TH2I*)fContainer->At(kQRes))) {
528 AliWarning("Missing dy=f(Q) histo");
529 return;
530 }
531 TF1 f("f", "gaus", -.5, .5);
532 TAxis *ax = 0x0;
533 TH1D *h1 = 0x0;
534
535 TObjArray *arr = (TObjArray*)fResults->At(kQRes);
536 TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
537 TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
538 TGraphErrors *gqp = (TGraphErrors*)arr->At(2);
539 Double_t q, n = 0., entries;
540 ax = h2->GetXaxis();
541 for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
542 Double_t q = TMath::Exp(ax->GetBinCenter(ix));
543 if(q<20. || q>250.) continue; // ?!
544
545 h1 = h2->ProjectionY("py", ix, ix);
546 entries = h1->GetEntries();
547 if(entries < 50) continue;
548 Adjust(&f, h1);
549 h1->Fit(&f, "Q");
550
551 // Fill sy^2 = f(q)
552 Int_t ip = gqm->GetN();
553 gqm->SetPoint(ip, q, 10.*f.GetParameter(1));
554 gqm->SetPointError(ip, 0., 10.*f.GetParError(1));
555
556 // correct sigma for ExB effect
557 gqs->SetPoint(ip, q, 1.e1*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
558 gqs->SetPointError(ip, 0., 1.e1*f.GetParError(2)/**f.GetParameter(2)*/);
559
560 // save probability
561 n += entries;
562 gqp->SetPoint(ip, q, entries);
563 gqp->SetPointError(ip, 0., 0./*TMath::Sqrt(entries)*/);
564 }
565
566 // normalize probability and get mean sy
567 Double_t sm = 0., sy;
568 for(Int_t ip=gqp->GetN(); ip--;){
569 gqp->GetPoint(ip, q, entries);
570 entries/=n;
571 gqp->SetPoint(ip, q, entries);
572 gqs->GetPoint(ip, q, sy);
573 sm += entries*sy;
574 }
575
576 // error parametrization s(q) = <sy> + b(1/q-1/q0)
577 TF1 fq("fq", "[0] + [1]/x", 20., 250.);
578 gqs->Fit(&fq);
579 printf("sy(Q) :: sm[%f] b[%f] 1/q0[%f]\n", sm, fq.GetParameter(1), (sm-fq.GetParameter(0))/fq.GetParameter(1));
580}
581
582//_______________________________________________________
583void AliTRDclusterResolution::ProcessCenterPad()
584{
585 TObjArray *arr = (TObjArray*)fContainer->At(kCenter);
586 if(!arr) {
587 AliWarning("Missing dy=f(y_d) container");
588 return;
589 }
590 TF1 f("f", "gaus", -.5, .5);
591 TAxis *ax = 0x0;
592 TH1D *h1 = 0x0;
593 TH2I *h2 = 0x0;
594
595 TObjArray *arrg = (TObjArray*)fResults->At(kCenter);
596 TH2F *hym = (TH2F*)arrg->At(0);
597 TH2F *hys = (TH2F*)arrg->At(1);
fa7831fb 598 TH2F *hyp = (TH2F*)arrg->At(2);
9462866a 599 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
600 if(!(h2 = (TH2I*)arr->At(ily))) continue;
601 ax = h2->GetXaxis();
602 for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
603 Float_t yd = ax->GetBinCenter(ix);
604 h1 = h2->ProjectionY("py", ix, ix);
fa7831fb 605 Int_t entries = (Int_t)h1->GetEntries();
606 if(entries < 50) continue;
607 //Adjust(&f, h1);
9462866a 608 h1->Fit(&f, "Q");
609
610 // Fill sy = f(y_w)
fa7831fb 611 hyp->Fill(ily, yd, entries);
9462866a 612 hym->Fill(ily, yd, f.GetParameter(1));
fa7831fb 613 //hym->SetPointError(ip, 0., f.GetParError(1));
9462866a 614 hys->Fill(ily, yd, f.GetParameter(2));
fa7831fb 615 //hys->SetPointError(ip, 0., f.GetParError(2));
9462866a 616 }
617 }
fa7831fb 618
619 // postprocess spectra
620 TCanvas *c = new TCanvas;
621 TF1 fgaus("fgaus", "gaus", -.5, .5);
622 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
623 h1 = hys->ProjectionY("pyy", ily+1, ily+1);
624
625 c->cd();
626 h1->Fit(&fgaus, "Q");
627 c->Modified(); c->Update(); gSystem->Sleep(500);
628 printf(" {%5.3e, %5.3e, %5.3e, ", fgaus.GetParameter(0), fgaus.GetParameter(1), fgaus.GetParameter(2));
629
630 Float_t sy = 0.;
631 h1 = hyp->ProjectionY("pyy", ily+1, ily+1);
632 for(Int_t ix=1; ix<=h1->GetNbinsX(); ix++){
633 sy += fgaus.Eval(h1->GetBinCenter(ix))*h1->GetBinContent(ix);
634// printf("ix[%d] gaus[%e] p[%d]\n", ix, fgaus.Eval(h1->GetBinCenter(ix)), (Int_t)h1->GetBinContent(ix));
635 }
636 sy /= h1->GetEntries();
637 printf("%5.3e},\n", sy);
638 }
639
9462866a 640}
641
642//_______________________________________________________
643void AliTRDclusterResolution::ProcessSigma()
644{
645 TObjArray *arr = (TObjArray*)fContainer->At(kSigm);
646 if(!arr){
647 AliWarning("Missing dy=f(x_d, d_wire) container");
648 return;
649 }
650 TLinearFitter gs(1,"pol1");
651 TF1 f("f", "gaus", -.5, .5);
652 TAxis *ax = 0x0;
653 TH1D *h1 = 0x0;
654 TH2I *h2 = 0x0;
655
656 Double_t d(0.), x(0.), sx, sy;
657 TObjArray *arrr = (TObjArray*)fResults->At(kSigm);
658 TH2F *hsx = (TH2F*)arrr->At(0);
659 TH2F *hsy = (TH2F*)arrr->At(1);
660 for(Int_t id=1; id<=fAd->GetNbins(); id++){
661 d = fAd->GetBinCenter(id); //[mm]
662 printf(" Doing d = %f[mm]\n", d);
663 for(Int_t it=1; it<=fAt->GetNbins(); it++){
664 x = fAt->GetBinCenter(it); //[mm]
665 printf(" Doing xd = %f[mm]\n", x);
666 Int_t idx = (id-1)*kNTB+it-1;
667
668 // retrieve data histogram
669 if(!(h2 = (TH2I*)arr->At(idx))) {
670 AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
671 continue;
672 }
673 gs.ClearPoints();
674 ax = h2->GetXaxis();
675 for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
676 Float_t dydx = ax->GetBinCenter(ix);
677 Double_t tgg = (dydx-fExB)/(1.+dydx*fExB);
678 Double_t tgg2 = tgg*tgg;
679 h1 = h2->ProjectionY("py", ix, ix);
680 if(h1->GetEntries() < 100) continue;
681 //Adjust(&f, h1);
682 printf("\tFit ix[%d] on %s entries [%d]\n", ix, h2->GetName(), (Int_t)h1->GetEntries());
683 h1->Fit(&f, "QN");
684
685 // Fill sy^2 = f(tg^2(phi-a_L))
686 gs.AddPoint(&tgg2, 1.e2*f.GetParameter(2)*f.GetParameter(2), 2.e2*f.GetParameter(2)*f.GetParError(2));
687 }
688 if(gs.Eval()) continue;
689 sx = TMath::Sqrt(gs.GetParameter(1));
690 sy = TMath::Sqrt(gs.GetParameter(0));
691 printf("sx[%f] sy[%f] [mm]\n", sx, sy);
692
693 hsx->SetBinContent(id, it, sx);
694 //hsx->SetBinError(id, it, sig.GetParError(1));
695
696 // s^2_y = s0^2_y + tg^2(a_L) * s^2_x
697 hsy->SetBinContent(id, it, sy/* - exb2*sig.GetParameter(1)*/);
698 //hsy->SetBinError(id, it, sig.GetParError(0)+exb2*exb2*sig.GetParError(1));
699 }
700 }
701
702 // fit sy parallel to the drift
703 h1 = hsy->ProjectionY("pyy"); h1->Scale(1./hsy->GetNbinsX());
704 TF1 fsyD("fsy", "[0]+exp([1]*(x+[2]))", 0.5, 3.5);
705 h1->Fit(&fsyD);
706
707
708 // fit sx parallel to the drift
709 h1 = hsx->ProjectionY("pyy"); h1->Scale(1./hsx->GetNbinsX());
710 TF1 fsxD("fsx", "gaus(0)+expo(3)", 0., 3.5);
711 h1->Fit(&fsxD);
712
713 // fit sx perpendicular to the drift
714 Int_t nb = 0;
715 TAxis *ay = hsx->GetYaxis();
716 TH1D *hx = hsx->ProjectionX("px", 1,1);
717 TH1D *hAn = (TH1D*)hx->Clone("hAn"); hAn->Reset();
718 for(Int_t ib = 11; ib<24; ib++, nb++){
719 hx = hsx->ProjectionX("px", ib, ib);
720 for(Int_t id=1; id<=hx->GetNbinsX(); id++)
721 hx->SetBinContent(id, hx->GetBinContent(id)-fsxD.Eval(ay->GetBinCenter(ib)));
722 hAn->Add(hx);
723 }
724 hAn->Scale(1./nb);
725 hAn->Fit("pol2");
726}
727
728//_______________________________________________________
729void AliTRDclusterResolution::ProcessMean()
730{
731 TObjArray *arr = (TObjArray*)fContainer->At(kMean);
732 if(!arr){
733 AliWarning("Missing dy=f(x_d, d_wire) container");
734 return;
735 }
736 TGraphErrors *gm = new TGraphErrors();
737 TF1 f("f", "gaus", -.5, .5);
738 TF1 line("l", Form("%6.4f*[0]+[1]", fExB), -.12, .12);
739 Double_t d(0.), x(0.), dx, x0;
740 TAxis *ax = 0x0;
741 TH1D *h1 = 0x0;
742 TH2I *h2 = 0x0;
743
744 TCanvas *c=new TCanvas;
745
746 TObjArray *arrr = (TObjArray*)fResults->At(kMean);
747 TH2F *hdx = (TH2F*)arrr->At(0);
748 TH2F *hx0 = (TH2F*)arrr->At(1);
749 for(Int_t id=1; id<=fAd->GetNbins(); id++){
750 d = fAd->GetBinCenter(id); //[mm]
751 printf(" Doing d = %f[mm]\n", d);
752 for(Int_t it=1; it<=fAt->GetNbins(); it++){
753 x = fAt->GetBinCenter(it); //[mm]
754 printf(" Doing xd = %f[mm]\n", x);
755 Int_t idx = (id-1)*kNTB+it-1;
756 if(!(h2 = (TH2I*)arr->At(idx))) {
757 AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
758 continue;
759 }
760 Int_t ip = 0;
761 ax = h2->GetXaxis();
762 for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
763 Double_t dydx = ax->GetBinCenter(ix);
764 h1 = h2->ProjectionY("py", ix, ix);
765 if(h1->GetEntries() < 50) continue;
766 //Adjust(&f, h1);
767 //printf("\tFit ix[%d] on %s entries [%d]\n", ix, h2->GetName(), (Int_t)h1->GetEntries());
768 h1->Fit(&f, "QN");
769
770 // Fill <Dy> = f(dydx - h*dzdx)
771 gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
772 gm->SetPointError(ip, 0., 10.*f.GetParError(1));
773 ip++;
774 }
775 if(gm->GetN()<4) continue;
776
777 gm->Fit(&line, "QN");
778 c->cd();
779 gm->Draw("apl");
780 c->Modified(); c->Update(); gSystem->Sleep(500);
781 dx = line.GetParameter(1); // = (fVdrift + dVd)*(fT + tCorr);
782 x0 = line.GetParameter(0); // = tg(a_L)*dx
783 printf(" dx[%f] x0[%f] [mm]\n", dx, x0);
784
785 hdx->SetBinContent(id, it, dx);
786 hx0->SetBinContent(id, it, x0);
787 }
788 }
789
790 TH1D *h=hdx->ProjectionY("py"); h->Scale(1./hdx->GetNbinsX());
791}