]>
Commit | Line | Data |
---|---|---|
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 | |
f2231ce8 | 128 | // #sigma_{y}^{||} = #sigma_{y}^{0} -a*exp(1/(x-b)) |
9462866a | 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(); | |
162637e4 | 145 | // cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); |
9462866a | 146 | // cdb->SetRun(0); |
147 | // // initialize magnetic field. | |
6fc46cba | 148 | // AliMagFCheb *field=new AliMagFCheb("Maps","Maps", 2, 1., 10., AliMagFCheb::k5kG); |
9462866a | 149 | // AliTracker::SetFieldMap(field, kTRUE); |
150 | // | |
151 | // AliTRDclusterResolution *res = new AliTRDclusterResolution(); | |
152 | // res->SetMCdata(); | |
153 | // res->Load("TRD.TaskClErrParam.root"); | |
6fc46cba | 154 | // res->SetExB(); |
155 | // res->SetVisual(); | |
156 | // //res->SetSaveAs(); | |
157 | // res->SetProcessCharge(kFALSE); | |
158 | // res->SetProcessCenterPad(kFALSE); | |
159 | // //res->SetProcessMean(kFALSE); | |
160 | // res->SetProcessSigma(kFALSE); | |
9462866a | 161 | // if(!res->PostProcess()) return; |
162 | // new TCanvas; | |
163 | // res->GetRefFigure(fig); | |
164 | // } | |
165 | // | |
166 | // Authors: // | |
167 | // Alexandru Bercuci <A.Bercuci@gsi.de> // | |
168 | //////////////////////////////////////////////////////////////////////////// | |
169 | ||
b2dc316d | 170 | #include "AliTRDclusterResolution.h" |
171 | #include "AliTRDtrackInfo/AliTRDclusterInfo.h" | |
5198d8c6 | 172 | #include "AliTRDgeometry.h" |
6bc4a8f4 | 173 | #include "AliTRDcalibDB.h" |
a076fc2f | 174 | #include "AliTRDCommonParam.h" |
6bc4a8f4 | 175 | #include "Cal/AliTRDCalROC.h" |
176 | #include "Cal/AliTRDCalDet.h" | |
b2dc316d | 177 | |
178 | #include "AliLog.h" | |
6bc4a8f4 | 179 | #include "AliTracker.h" |
180 | #include "AliCDBManager.h" | |
b2dc316d | 181 | |
f2231ce8 | 182 | #include "TROOT.h" |
b2dc316d | 183 | #include "TObjArray.h" |
184 | #include "TAxis.h" | |
fc0946a7 | 185 | #include "TF1.h" |
186 | #include "TGraphErrors.h" | |
f2231ce8 | 187 | #include "TLine.h" |
b2dc316d | 188 | #include "TH2I.h" |
189 | #include "TMath.h" | |
9462866a | 190 | #include "TLinearFitter.h" |
191 | ||
192 | #include "TCanvas.h" | |
193 | #include "TSystem.h" | |
194 | ||
b2dc316d | 195 | |
196 | ClassImp(AliTRDclusterResolution) | |
197 | ||
ec3f0161 | 198 | const Float_t AliTRDclusterResolution::fgkTimeBinLength = 1./ AliTRDCommonParam::Instance()->GetSamplingFrequency(); |
b2dc316d | 199 | //_______________________________________________________ |
6fc46cba | 200 | AliTRDclusterResolution::AliTRDclusterResolution(const char *name) |
201 | : AliTRDrecoTask(name, "Cluster Error Parametrization") | |
f2231ce8 | 202 | ,fCanvas(0x0) |
b2dc316d | 203 | ,fInfo(0x0) |
fc0946a7 | 204 | ,fResults(0x0) |
205 | ,fAt(0x0) | |
206 | ,fAd(0x0) | |
f2231ce8 | 207 | ,fStatus(0) |
208 | ,fDet(-1) | |
fc0946a7 | 209 | ,fExB(0.) |
9462866a | 210 | ,fVdrift(0.) |
b2dc316d | 211 | { |
ec3f0161 | 212 | // time drift axis |
213 | fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength); | |
214 | // z axis spans the drift cell 2.5 mm | |
fc0946a7 | 215 | fAd = new TAxis(kND, 0., .25); |
f2231ce8 | 216 | |
217 | // By default register all analysis | |
218 | // The user can switch them off in his steering macro | |
ec3f0161 | 219 | SetProcess(kQRes); |
220 | SetProcess(kCenter); | |
221 | SetProcess(kMean); | |
222 | SetProcess(kSigm); | |
b2dc316d | 223 | } |
224 | ||
225 | //_______________________________________________________ | |
226 | AliTRDclusterResolution::~AliTRDclusterResolution() | |
227 | { | |
f2231ce8 | 228 | if(fCanvas) delete fCanvas; |
fc0946a7 | 229 | if(fAd) delete fAd; |
230 | if(fAt) delete fAt; | |
231 | if(fResults){ | |
232 | fResults->Delete(); | |
233 | delete fResults; | |
234 | } | |
b2dc316d | 235 | } |
236 | ||
237 | //_______________________________________________________ | |
238 | void AliTRDclusterResolution::ConnectInputData(Option_t *) | |
239 | { | |
240 | fInfo = dynamic_cast<TObjArray *>(GetInputData(0)); | |
241 | } | |
242 | ||
243 | //_______________________________________________________ | |
244 | void AliTRDclusterResolution::CreateOutputObjects() | |
245 | { | |
246 | OpenFile(0, "RECREATE"); | |
247 | fContainer = Histos(); | |
248 | } | |
249 | ||
250 | //_______________________________________________________ | |
e15179be | 251 | Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig) |
b2dc316d | 252 | { |
e15179be | 253 | if(!fResults) return kFALSE; |
fc0946a7 | 254 | |
f2231ce8 | 255 | TList *l = 0x0; |
fc0946a7 | 256 | TObjArray *arr = 0x0; |
d667707c | 257 | TH2 *h2 = 0x0;TH1 *h1 = 0x0; |
9462866a | 258 | TGraphErrors *gm(0x0), *gs(0x0), *gp(0x0); |
fc0946a7 | 259 | switch(ifig){ |
260 | case kQRes: | |
261 | if(!(arr = (TObjArray*)fResults->At(kQRes))) break; | |
262 | if(!(gm = (TGraphErrors*)arr->At(0))) break; | |
263 | if(!(gs = (TGraphErrors*)arr->At(1))) break; | |
9462866a | 264 | if(!(gp = (TGraphErrors*)arr->At(2))) break; |
fc0946a7 | 265 | gs->Draw("apl"); |
9462866a | 266 | gs->GetHistogram()->GetYaxis()->SetRangeUser(-.01, .6); |
267 | gs->GetHistogram()->SetXTitle("Q [a.u.]"); | |
268 | gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [mm] / freq"); | |
269 | gm->Draw("pl"); | |
270 | gp->Draw("pl"); | |
e15179be | 271 | return kTRUE; |
9462866a | 272 | case kCenter: |
273 | if(!(arr = (TObjArray*)fResults->At(kCenter))) break; | |
f2231ce8 | 274 | gPad->Divide(3, 1); l = gPad->GetListOfPrimitives(); |
275 | for(Int_t ipad = 3; ipad--;){ | |
276 | if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE; | |
277 | ((TVirtualPad*)l->At(ipad))->cd(); | |
278 | h2->Draw("lego2fb"); | |
279 | } | |
e15179be | 280 | return kTRUE; |
9462866a | 281 | case kSigm: |
282 | if(!(arr = (TObjArray*)fResults->At(kSigm))) break; | |
f2231ce8 | 283 | gPad->Divide(2, 1); l = gPad->GetListOfPrimitives(); |
d667707c | 284 | if(!(h2 = (TH2F*)arr->At(0))) return kFALSE; |
285 | ((TVirtualPad*)l->At(0))->cd(); | |
286 | h1 = h2->ProjectionY("hsx_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24); | |
287 | h1->SetYTitle("<#sigma_{x}> [#mum]"); | |
288 | h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc"); | |
289 | ||
290 | if(!(h2 = (TH2F*)arr->At(1))) return kFALSE; | |
291 | ((TVirtualPad*)l->At(1))->cd(); | |
292 | h1 = h2->ProjectionY("hsy_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24); | |
293 | h1->SetYTitle("<#sigma_{y}> [#mum]"); | |
294 | h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc"); | |
9462866a | 295 | return kTRUE; |
296 | case kMean: | |
297 | if(!(arr = (TObjArray*)fResults->At(kMean))) break; | |
f2231ce8 | 298 | gPad->Divide(2, 1); l = gPad->GetListOfPrimitives(); |
d667707c | 299 | ((TVirtualPad*)l->At(0))->cd(); |
300 | if(!(gm = (TGraphErrors*)arr->At(0))) return kFALSE; | |
301 | gm->Draw("apl"); | |
302 | gm->GetHistogram()->SetXTitle("t_{drift} [#mus]"); | |
303 | gm->GetHistogram()->SetYTitle("dx [#mum]"); | |
304 | ||
305 | ((TVirtualPad*)l->At(1))->cd(); | |
306 | if(!(gm = (TGraphErrors*)arr->At(1))) return kFALSE; | |
307 | gm->Draw("apl"); | |
308 | gm->GetHistogram()->SetXTitle("t_{drift} [#mus]"); | |
309 | gm->GetHistogram()->SetYTitle("dy [#mum]"); | |
310 | ||
e15179be | 311 | return kTRUE; |
fc0946a7 | 312 | default: |
313 | break; | |
314 | } | |
9462866a | 315 | AliWarning("No container/data found."); |
e15179be | 316 | return kFALSE; |
b2dc316d | 317 | } |
318 | ||
319 | //_______________________________________________________ | |
320 | TObjArray* AliTRDclusterResolution::Histos() | |
321 | { | |
322 | if(fContainer) return fContainer; | |
6fc46cba | 323 | fContainer = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainer)); |
b2dc316d | 324 | //fContainer->SetOwner(kTRUE); |
325 | ||
5198d8c6 | 326 | TH2I *h2 = 0x0; |
327 | TObjArray *arr = 0x0; | |
328 | ||
329 | fContainer->AddAt(h2 = new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), kQRes); | |
330 | h2->SetXTitle("log(q) [a.u.]"); | |
331 | h2->SetYTitle("#Delta y[cm]"); | |
332 | h2->SetZTitle("entries"); | |
333 | ||
9462866a | 334 | fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kCenter); |
6fc46cba | 335 | arr->SetName("y(PadWidth)"); |
5198d8c6 | 336 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ |
6bc4a8f4 | 337 | arr->AddAt(h2 = new TH2I(Form("h_y%d", ily), Form("Ly[%d]", ily), 51, -.51, .51, 100, -.5, .5), ily); |
5198d8c6 | 338 | h2->SetXTitle("y_{w} [w]"); |
339 | h2->SetYTitle("#Delta y[cm]"); | |
340 | h2->SetZTitle("entries"); | |
341 | } | |
342 | ||
9462866a | 343 | fContainer->AddAt(arr = new TObjArray(kN), kSigm); |
6fc46cba | 344 | arr->SetName("Resolution"); |
b2dc316d | 345 | Int_t ih = 0; |
fc0946a7 | 346 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ |
347 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ | |
ec3f0161 | 348 | arr->AddAt(h2 = new TH2I(Form("hr_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++); |
5198d8c6 | 349 | h2->SetXTitle("tg#phi"); |
350 | h2->SetYTitle("#Delta y[cm]"); | |
351 | h2->SetZTitle("entries"); | |
b2dc316d | 352 | } |
353 | } | |
9462866a | 354 | |
355 | fContainer->AddAt(arr = new TObjArray(kN), kMean); | |
6fc46cba | 356 | arr->SetName("Systematics"); |
9462866a | 357 | ih = 0; |
358 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ | |
359 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ | |
ec3f0161 | 360 | arr->AddAt(h2 = new TH2I(Form("hs_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++); |
9462866a | 361 | h2->SetXTitle("tg#phi-h*tg(#theta)"); |
362 | h2->SetYTitle("#Delta y[cm]"); | |
363 | h2->SetZTitle("entries"); | |
364 | } | |
365 | } | |
366 | ||
b2dc316d | 367 | return fContainer; |
368 | } | |
369 | ||
370 | //_______________________________________________________ | |
371 | void AliTRDclusterResolution::Exec(Option_t *) | |
372 | { | |
6bc4a8f4 | 373 | if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the task."); |
374 | ||
87b166d3 | 375 | Int_t det, t; |
fc0946a7 | 376 | Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3]; |
b2dc316d | 377 | TH2I *h2 = 0x0; |
5198d8c6 | 378 | |
9462866a | 379 | // define limits around ExB for which x contribution is negligible |
380 | const Float_t kDtgPhi = 3.5e-2; //(+- 2 deg) | |
381 | ||
382 | TObjArray *arr0 = (TObjArray*)fContainer->At(kCenter); | |
383 | TObjArray *arr1 = (TObjArray*)fContainer->At(kSigm); | |
384 | TObjArray *arr2 = (TObjArray*)fContainer->At(kMean); | |
5198d8c6 | 385 | |
b2dc316d | 386 | const AliTRDclusterInfo *cli = 0x0; |
387 | TIterator *iter=fInfo->MakeIterator(); | |
388 | while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){ | |
6bc4a8f4 | 389 | cli->GetCluster(det, x, y, z, q, t, covcl); |
390 | if(fDet>=0 && fDet!=det) continue; | |
9462866a | 391 | |
392 | dy = cli->GetResolution(); | |
393 | cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]); | |
394 | ||
395 | // resolution as a function of cluster charge | |
396 | // only for phi equal exB | |
397 | if(TMath::Abs(dydx-fExB) < kDtgPhi){ | |
398 | h2 = (TH2I*)fContainer->At(kQRes); | |
399 | h2->Fill(TMath::Log(q), dy); | |
400 | } | |
401 | ||
402 | // do not use problematic clusters in resolution analysis | |
403 | // TODO define limits as calibration aware (gain) !! | |
404 | if(q<20. || q>250.) continue; | |
405 | ||
406 | // resolution as a function of y displacement from pad center | |
407 | // only for phi equal exB | |
408 | if(TMath::Abs(dydx-fExB) < kDtgPhi){ | |
409 | h2 = (TH2I*)arr0->At(AliTRDgeometry::GetLayer(det)); | |
410 | h2->Fill(cli->GetYDisplacement(), dy); | |
411 | } | |
6bc4a8f4 | 412 | |
ec3f0161 | 413 | Int_t it = fAt->FindBin((t+.5)*fgkTimeBinLength); |
fc0946a7 | 414 | if(it==0 || it == fAt->GetNbins()+1){ |
ec3f0161 | 415 | AliWarning(Form("Drift time %f outside allowed range", t)); |
b2dc316d | 416 | continue; |
417 | } | |
fc0946a7 | 418 | Int_t id = fAd->FindBin(cli->GetAnisochronity()); |
419 | if(id==0 || id == fAd->GetNbins()+1){ | |
b2dc316d | 420 | AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity())); |
421 | continue; | |
422 | } | |
9462866a | 423 | // calculate index of cluster position in array |
424 | Int_t hid = (id-1)*kNTB+it-1; | |
b2dc316d | 425 | |
9462866a | 426 | // fill histo for resolution (sigma) |
427 | h2 = (TH2I*)arr1->At(hid); | |
6bc4a8f4 | 428 | h2->Fill(dydx, dy); |
b2dc316d | 429 | |
9462866a | 430 | // fill histo for systematic (mean) |
431 | h2 = (TH2I*)arr2->At(hid); | |
432 | h2->Fill(dydx-cli->GetTilt()*dzdx, dy); | |
b2dc316d | 433 | } |
434 | PostData(0, fContainer); | |
435 | } | |
436 | ||
437 | ||
438 | //_______________________________________________________ | |
439 | Bool_t AliTRDclusterResolution::PostProcess() | |
440 | { | |
fc0946a7 | 441 | if(!fContainer) return kFALSE; |
6bc4a8f4 | 442 | if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the post processing."); |
fc0946a7 | 443 | |
5198d8c6 | 444 | TObjArray *arr = 0x0; |
fc0946a7 | 445 | TH2 *h2 = 0x0; |
446 | if(!fResults){ | |
447 | TGraphErrors *g = 0x0; | |
6fc46cba | 448 | fResults = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainer)); |
fc0946a7 | 449 | fResults->SetOwner(); |
9462866a | 450 | fResults->AddAt(arr = new TObjArray(3), kQRes); |
fc0946a7 | 451 | arr->SetOwner(); |
452 | arr->AddAt(g = new TGraphErrors(), 0); | |
453 | g->SetLineColor(kBlue); g->SetMarkerColor(kBlue); | |
454 | g->SetMarkerStyle(7); | |
455 | arr->AddAt(g = new TGraphErrors(), 1); | |
456 | g->SetLineColor(kRed); g->SetMarkerColor(kRed); | |
457 | g->SetMarkerStyle(23); | |
9462866a | 458 | arr->AddAt(g = new TGraphErrors(), 2); |
459 | g->SetLineColor(kGreen); g->SetMarkerColor(kGreen); | |
fc0946a7 | 460 | g->SetMarkerStyle(7); |
fc0946a7 | 461 | |
fa7831fb | 462 | fResults->AddAt(arr = new TObjArray(3), kCenter); |
9462866a | 463 | arr->SetOwner(); |
f2231ce8 | 464 | |
465 | if(!(h2 = (TH2F*)gROOT->FindObject("hYM"))){ | |
466 | h2 = new TH2F("hYM", "", | |
467 | AliTRDgeometry::kNlayer, -.5, AliTRDgeometry::kNlayer-.5, 51, -.51, .51); | |
468 | } | |
469 | arr->AddAt(h2, 0); | |
9462866a | 470 | h2->SetXTitle("ly"); |
471 | h2->SetYTitle("y [w]"); | |
f2231ce8 | 472 | h2->SetZTitle("#mu_{x} [cm]"); |
9462866a | 473 | arr->AddAt(h2 = (TH2F*)h2->Clone("hYS"), 1); |
f2231ce8 | 474 | h2->SetZTitle("#sigma_{x} [cm]"); |
fa7831fb | 475 | arr->AddAt(h2 = (TH2F*)h2->Clone("hYP"), 2); |
476 | h2->SetZTitle("entries"); | |
9462866a | 477 | |
478 | fResults->AddAt(arr = new TObjArray(2), kSigm); | |
479 | arr->SetOwner(); | |
f2231ce8 | 480 | if(!(h2 = (TH2F*)gROOT->FindObject("hSX"))){ |
481 | h2 = new TH2F("hSX", "", | |
fc0946a7 | 482 | fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(), |
f2231ce8 | 483 | fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()); |
484 | } | |
485 | arr->AddAt(h2, 0); | |
81a6494d | 486 | h2->SetXTitle("d [cm]"); |
ec3f0161 | 487 | h2->SetYTitle("t_{drift} [#mus]"); |
f2231ce8 | 488 | h2->SetZTitle("#sigma_{x} [cm]"); |
9462866a | 489 | arr->AddAt(h2 = (TH2F*)h2->Clone("hSY"), 1); |
f2231ce8 | 490 | h2->SetZTitle("#sigma_{y} [cm]"); |
9462866a | 491 | |
d667707c | 492 | fResults->AddAt(arr = new TObjArray(4), kMean); |
9462866a | 493 | arr->SetOwner(); |
d667707c | 494 | arr->AddAt(g = new TGraphErrors(), 0); |
495 | g->SetLineColor(kBlue); g->SetMarkerColor(kBlue); | |
496 | g->SetMarkerStyle(24); | |
497 | arr->AddAt(g = new TGraphErrors(), 1); | |
498 | g->SetLineColor(kRed); g->SetMarkerColor(kRed); | |
499 | g->SetMarkerStyle(24); | |
500 | arr->AddAt(h2 = (TH2F*)h2->Clone("hDX"), 2); | |
f2231ce8 | 501 | h2->SetZTitle("dx [cm]"); |
d667707c | 502 | arr->AddAt(h2 = (TH2F*)h2->Clone("hDY"), 3); |
81a6494d | 503 | h2->SetZTitle("dy [cm]"); |
fc0946a7 | 504 | } else { |
505 | TObject *o = 0x0; | |
506 | TIterator *iter=fResults->MakeIterator(); | |
5198d8c6 | 507 | while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point |
fc0946a7 | 508 | } |
9462866a | 509 | |
510 | // precalculated value of tg^2(alpha_L) | |
511 | Double_t exb2 = fExB*fExB; | |
512 | // square of the mean value of sigma drift length. | |
513 | // has to come from previous calibration | |
514 | //Double_t sxd2 = 1.;// [mm^2] | |
fc0946a7 | 515 | |
9462866a | 516 | printf("ExB[%e] ExB2[%e]\n", fExB, exb2); |
fc0946a7 | 517 | |
518 | // process resolution dependency on charge | |
ec3f0161 | 519 | if(HasProcess(kQRes)) ProcessCharge(); |
fc0946a7 | 520 | |
fc0946a7 | 521 | // process resolution dependency on y displacement |
ec3f0161 | 522 | if(HasProcess(kCenter)) ProcessCenterPad(); |
fc0946a7 | 523 | |
5198d8c6 | 524 | // process resolution dependency on drift legth and drift cell width |
ec3f0161 | 525 | if(HasProcess(kSigm)) ProcessSigma(); |
9462866a | 526 | |
527 | // process systematic shift on drift legth and drift cell width | |
ec3f0161 | 528 | if(HasProcess(kMean)) ProcessMean(); |
fc0946a7 | 529 | |
530 | return kTRUE; | |
b2dc316d | 531 | } |
532 | ||
6bc4a8f4 | 533 | //_______________________________________________________ |
9462866a | 534 | Bool_t AliTRDclusterResolution::SetExB(Int_t det, Int_t col, Int_t row) |
6bc4a8f4 | 535 | { |
536 | // check OCDB | |
537 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
538 | if(cdb->GetRun() < 0){ | |
539 | AliError("OCDB manager not properly initialized"); | |
540 | return kFALSE; | |
541 | } | |
542 | ||
543 | // check magnetic field | |
544 | if(TMath::Abs(AliTracker::GetBz()) < 1.e-10){ | |
545 | AliWarning("B=0. Magnetic field may not be initialized. Continue if you know what you are doing !"); | |
546 | } | |
547 | ||
548 | // set reference detector if any | |
549 | if(det>=0 && det<AliTRDgeometry::kNdet) fDet = det; | |
550 | else det = 0; | |
551 | ||
552 | AliTRDcalibDB *fCalibration = AliTRDcalibDB::Instance(); | |
553 | AliTRDCalROC *fCalVdriftROC = fCalibration->GetVdriftROC(det); | |
554 | const AliTRDCalDet *fCalVdriftDet = fCalibration->GetVdriftDet(); | |
555 | ||
9462866a | 556 | fVdrift = fCalVdriftDet->GetValue(det) * fCalVdriftROC->GetValue(col, row); |
a076fc2f | 557 | fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVdrift); |
6bc4a8f4 | 558 | SetBit(kExB); |
559 | return kTRUE; | |
560 | } | |
b2dc316d | 561 | |
f2231ce8 | 562 | //_______________________________________________________ |
563 | void AliTRDclusterResolution::SetVisual() | |
564 | { | |
565 | if(fCanvas) return; | |
566 | fCanvas = new TCanvas("clResCanvas", "Cluster Resolution Visualization", 10, 10, 600, 600); | |
567 | } | |
568 | ||
9462866a | 569 | //_______________________________________________________ |
570 | void AliTRDclusterResolution::ProcessCharge() | |
571 | { | |
572 | TH2I *h2 = 0x0; | |
573 | if((h2 = (TH2I*)fContainer->At(kQRes))) { | |
574 | AliWarning("Missing dy=f(Q) histo"); | |
575 | return; | |
576 | } | |
577 | TF1 f("f", "gaus", -.5, .5); | |
578 | TAxis *ax = 0x0; | |
579 | TH1D *h1 = 0x0; | |
580 | ||
581 | TObjArray *arr = (TObjArray*)fResults->At(kQRes); | |
582 | TGraphErrors *gqm = (TGraphErrors*)arr->At(0); | |
583 | TGraphErrors *gqs = (TGraphErrors*)arr->At(1); | |
584 | TGraphErrors *gqp = (TGraphErrors*)arr->At(2); | |
585 | Double_t q, n = 0., entries; | |
586 | ax = h2->GetXaxis(); | |
587 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
2a12b21f | 588 | q = TMath::Exp(ax->GetBinCenter(ix)); |
9462866a | 589 | if(q<20. || q>250.) continue; // ?! |
590 | ||
591 | h1 = h2->ProjectionY("py", ix, ix); | |
592 | entries = h1->GetEntries(); | |
593 | if(entries < 50) continue; | |
594 | Adjust(&f, h1); | |
595 | h1->Fit(&f, "Q"); | |
596 | ||
597 | // Fill sy^2 = f(q) | |
598 | Int_t ip = gqm->GetN(); | |
599 | gqm->SetPoint(ip, q, 10.*f.GetParameter(1)); | |
600 | gqm->SetPointError(ip, 0., 10.*f.GetParError(1)); | |
601 | ||
602 | // correct sigma for ExB effect | |
603 | gqs->SetPoint(ip, q, 1.e1*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/); | |
604 | gqs->SetPointError(ip, 0., 1.e1*f.GetParError(2)/**f.GetParameter(2)*/); | |
605 | ||
606 | // save probability | |
607 | n += entries; | |
608 | gqp->SetPoint(ip, q, entries); | |
609 | gqp->SetPointError(ip, 0., 0./*TMath::Sqrt(entries)*/); | |
610 | } | |
611 | ||
612 | // normalize probability and get mean sy | |
613 | Double_t sm = 0., sy; | |
614 | for(Int_t ip=gqp->GetN(); ip--;){ | |
615 | gqp->GetPoint(ip, q, entries); | |
616 | entries/=n; | |
617 | gqp->SetPoint(ip, q, entries); | |
618 | gqs->GetPoint(ip, q, sy); | |
619 | sm += entries*sy; | |
620 | } | |
621 | ||
622 | // error parametrization s(q) = <sy> + b(1/q-1/q0) | |
623 | TF1 fq("fq", "[0] + [1]/x", 20., 250.); | |
624 | gqs->Fit(&fq); | |
625 | printf("sy(Q) :: sm[%f] b[%f] 1/q0[%f]\n", sm, fq.GetParameter(1), (sm-fq.GetParameter(0))/fq.GetParameter(1)); | |
626 | } | |
627 | ||
628 | //_______________________________________________________ | |
629 | void AliTRDclusterResolution::ProcessCenterPad() | |
630 | { | |
631 | TObjArray *arr = (TObjArray*)fContainer->At(kCenter); | |
632 | if(!arr) { | |
633 | AliWarning("Missing dy=f(y_d) container"); | |
634 | return; | |
635 | } | |
636 | TF1 f("f", "gaus", -.5, .5); | |
637 | TAxis *ax = 0x0; | |
f2231ce8 | 638 | TH1D *h1 = 0x0, *h11 = 0x0; |
9462866a | 639 | TH2I *h2 = 0x0; |
640 | ||
641 | TObjArray *arrg = (TObjArray*)fResults->At(kCenter); | |
642 | TH2F *hym = (TH2F*)arrg->At(0); | |
643 | TH2F *hys = (TH2F*)arrg->At(1); | |
fa7831fb | 644 | TH2F *hyp = (TH2F*)arrg->At(2); |
9462866a | 645 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ |
646 | if(!(h2 = (TH2I*)arr->At(ily))) continue; | |
647 | ax = h2->GetXaxis(); | |
648 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
649 | Float_t yd = ax->GetBinCenter(ix); | |
650 | h1 = h2->ProjectionY("py", ix, ix); | |
fa7831fb | 651 | Int_t entries = (Int_t)h1->GetEntries(); |
652 | if(entries < 50) continue; | |
653 | //Adjust(&f, h1); | |
f2231ce8 | 654 | h1->Fit(&f, "QN"); |
9462866a | 655 | |
656 | // Fill sy = f(y_w) | |
fa7831fb | 657 | hyp->Fill(ily, yd, entries); |
9462866a | 658 | hym->Fill(ily, yd, f.GetParameter(1)); |
fa7831fb | 659 | //hym->SetPointError(ip, 0., f.GetParError(1)); |
9462866a | 660 | hys->Fill(ily, yd, f.GetParameter(2)); |
fa7831fb | 661 | //hys->SetPointError(ip, 0., f.GetParError(2)); |
9462866a | 662 | } |
663 | } | |
fa7831fb | 664 | |
f2231ce8 | 665 | // POSTPROCESS SPECTRA |
666 | // Found correction for systematic deviation | |
667 | TF1 fprf("fprf", "[0]+[1]*sin([2]*x)", -.5, .5); | |
668 | fprf.SetParameter(0, 0.); | |
669 | fprf.SetParameter(1, 1.1E-2); | |
670 | fprf.SetParameter(2, -TMath::PiOver2()/0.25); | |
671 | printf(" const Float_t cy[AliTRDgeometry::kNlayer][3] = {\n"); | |
672 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ | |
673 | h1 = hym->ProjectionY("hym_pyy", ily+1, ily+1); | |
674 | // adjust errors | |
675 | for(Int_t ib=h1->GetNbinsX(); ib--;) h1->SetBinError(ib, 0.002); | |
676 | h1->Fit(&fprf, "Q"); | |
677 | printf(" {%5.3e, %5.3e, %5.3e},\n", fprf.GetParameter(0), fprf.GetParameter(1), fprf.GetParameter(2)); | |
678 | ||
679 | if(!fCanvas) continue; | |
680 | fCanvas->cd(); | |
681 | h1->SetMinimum(-0.02);h1->SetMaximum(0.02);h1->Draw("e1"); | |
682 | h11 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1); | |
683 | h11->Scale(.8/h11->Integral()); | |
684 | h11->SetLineColor(kBlue); h11->Draw("csame"); | |
685 | fCanvas->Modified(); fCanvas->Update(); | |
686 | if(IsSaveAs()) | |
687 | fCanvas->SaveAs(Form("Figures/ProcessCenterPad_M_Ly%d.gif", ily)); | |
688 | else gSystem->Sleep(500); | |
689 | } | |
690 | printf(" };\n"); | |
691 | ||
692 | // Parameterization for sigma PRF | |
fa7831fb | 693 | TF1 fgaus("fgaus", "gaus", -.5, .5); |
f2231ce8 | 694 | printf(" const Float_t scy[AliTRDgeometry::kNlayer][4] = {\n"); |
fa7831fb | 695 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ |
f2231ce8 | 696 | h1 = hys->ProjectionY("hys_pyy", ily+1, ily+1); |
697 | // adjust errors | |
698 | for(Int_t ib=h1->GetNbinsX(); ib--;) h1->SetBinError(ib, 0.002); | |
fa7831fb | 699 | |
fa7831fb | 700 | h1->Fit(&fgaus, "Q"); |
f2231ce8 | 701 | printf(" {%5.3e, %5.3e, %5.3e, ", fgaus.GetParameter(0), fgaus.GetParameter(1), fgaus.GetParameter(2)); |
fa7831fb | 702 | |
f2231ce8 | 703 | // calculate mean sigma on the pad center distribution |
fa7831fb | 704 | Float_t sy = 0.; |
f2231ce8 | 705 | h1 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1); |
fa7831fb | 706 | for(Int_t ix=1; ix<=h1->GetNbinsX(); ix++){ |
707 | sy += fgaus.Eval(h1->GetBinCenter(ix))*h1->GetBinContent(ix); | |
fa7831fb | 708 | } |
709 | sy /= h1->GetEntries(); | |
710 | printf("%5.3e},\n", sy); | |
fa7831fb | 711 | |
f2231ce8 | 712 | if(!fCanvas) continue; |
713 | fCanvas->cd(); | |
714 | h1->SetMinimum(0.01);h1->SetMaximum(0.04);h1->Draw("e1"); | |
715 | h11 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1); | |
716 | h11->Scale(1./h11->Integral()); | |
717 | h11->SetLineColor(kBlue); h11->Draw("csame"); | |
718 | fCanvas->Modified(); fCanvas->Update(); | |
719 | if(IsSaveAs()) | |
720 | fCanvas->SaveAs(Form("Figures/ProcessCenterPad_S_Ly%d.gif", ily)); | |
721 | else gSystem->Sleep(500); | |
722 | } | |
723 | printf(" };\n"); | |
9462866a | 724 | } |
725 | ||
726 | //_______________________________________________________ | |
727 | void AliTRDclusterResolution::ProcessSigma() | |
728 | { | |
729 | TObjArray *arr = (TObjArray*)fContainer->At(kSigm); | |
730 | if(!arr){ | |
731 | AliWarning("Missing dy=f(x_d, d_wire) container"); | |
732 | return; | |
733 | } | |
734 | TLinearFitter gs(1,"pol1"); | |
735 | TF1 f("f", "gaus", -.5, .5); | |
736 | TAxis *ax = 0x0; | |
737 | TH1D *h1 = 0x0; | |
738 | TH2I *h2 = 0x0; | |
d667707c | 739 | TH1 *hFrame=0x0; |
740 | ||
f2231ce8 | 741 | // init visualization |
742 | TGraphErrors *ggs = 0x0; | |
d667707c | 743 | TGraph *line = 0x0; |
f2231ce8 | 744 | if(fCanvas){ |
745 | ggs = new TGraphErrors(); | |
d667707c | 746 | line = new TGraph(); |
747 | line->SetLineColor(kRed);line->SetLineWidth(2); | |
f2231ce8 | 748 | } |
9462866a | 749 | |
d667707c | 750 | Double_t d(0.), x(0.), sx, sy, exb2=0.;//fExB*fExB; |
9462866a | 751 | TObjArray *arrr = (TObjArray*)fResults->At(kSigm); |
752 | TH2F *hsx = (TH2F*)arrr->At(0); | |
753 | TH2F *hsy = (TH2F*)arrr->At(1); | |
754 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ | |
755 | d = fAd->GetBinCenter(id); //[mm] | |
f2231ce8 | 756 | printf(" Doing d = %5.3f [mm]\n", d); |
9462866a | 757 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ |
758 | x = fAt->GetBinCenter(it); //[mm] | |
9462866a | 759 | Int_t idx = (id-1)*kNTB+it-1; |
760 | ||
761 | // retrieve data histogram | |
762 | if(!(h2 = (TH2I*)arr->At(idx))) { | |
763 | AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d)); | |
764 | continue; | |
765 | } | |
f2231ce8 | 766 | |
767 | if(fCanvas){ | |
768 | new(ggs) TGraphErrors(); | |
769 | ggs->SetMarkerStyle(7); | |
770 | } | |
9462866a | 771 | gs.ClearPoints(); |
772 | ax = h2->GetXaxis(); | |
773 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
774 | Float_t dydx = ax->GetBinCenter(ix); | |
d667707c | 775 | //if(TMath::Abs(dydx)>0.18) continue; |
9462866a | 776 | Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); |
777 | Double_t tgg2 = tgg*tgg; | |
778 | h1 = h2->ProjectionY("py", ix, ix); | |
779 | if(h1->GetEntries() < 100) continue; | |
780 | //Adjust(&f, h1); | |
f2231ce8 | 781 | //printf("\tFit ix[%d] on %s entries [%d]\n", ix, h2->GetName(), (Int_t)h1->GetEntries()); |
9462866a | 782 | h1->Fit(&f, "QN"); |
f2231ce8 | 783 | Double_t s2 = f.GetParameter(2)*f.GetParameter(2); |
784 | Double_t s2e = 2.*f.GetParameter(2)*f.GetParError(2); | |
9462866a | 785 | // Fill sy^2 = f(tg^2(phi-a_L)) |
f2231ce8 | 786 | gs.AddPoint(&tgg2, s2, s2e); |
787 | ||
788 | if(!ggs) continue; | |
789 | Int_t ip = ggs->GetN(); | |
790 | ggs->SetPoint(ip, tgg2, s2); | |
791 | ggs->SetPointError(ip, 0., s2e); | |
9462866a | 792 | } |
793 | if(gs.Eval()) continue; | |
d667707c | 794 | |
795 | // s^2_x = s0^2_x - x^2*tg^2(a_L)/12 | |
796 | sx = gs.GetParameter(1)/* - x*x*exb2/12.*/; | |
797 | if(sx<0.) continue; | |
f2231ce8 | 798 | hsx->SetBinContent(id, it, TMath::Sqrt(sx)); |
d667707c | 799 | //hsx->SetBinError(id, it, .5*gs.GetParError(1)/TMath::Sqrt(sx)); |
9462866a | 800 | |
d667707c | 801 | // s^2_y = s0^2_y + tg^2(a_L) * s^2_x |
802 | // s0^2_y = f(D_L)*x + s_PRF^2 | |
803 | sy= gs.GetParameter(0)/*-exb2*sx*/; | |
804 | if(sy <0.) continue; | |
805 | hsy->SetBinContent(id, it, TMath::Sqrt(sy)); | |
9462866a | 806 | //hsy->SetBinError(id, it, sig.GetParError(0)+exb2*exb2*sig.GetParError(1)); |
f2231ce8 | 807 | |
d667707c | 808 | if(!fCanvas) continue; |
809 | fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy(); | |
810 | if(!hFrame){ | |
811 | hFrame=new TH1I("hFrame", "", 100, 0., .3); | |
812 | hFrame->SetMinimum(0.);hFrame->SetMaximum(.005); | |
813 | hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})"); | |
814 | hFrame->SetYTitle("#sigma^{2}y[cm^{2}]"); | |
815 | hFrame->SetLineColor(1);hFrame->SetLineWidth(1); | |
816 | hFrame->Draw(); | |
817 | } else hFrame->Reset(); | |
818 | Double_t xx = 0., dxx=.2/50; | |
819 | for(Int_t ip=0;ip<50;ip++){ | |
820 | line->SetPoint(ip, xx, gs.GetParameter(0)+xx*gs.GetParameter(1)); | |
821 | xx+=dxx; | |
822 | } | |
823 | ggs->Draw("pl"); line->Draw("l"); | |
f2231ce8 | 824 | fCanvas->Modified(); fCanvas->Update(); |
825 | if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_D%d_T%02d.gif", id, it)); | |
826 | else gSystem->Sleep(100); | |
827 | ||
d667707c | 828 | printf(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", x, TMath::Sqrt(sx), TMath::Sqrt(sy)); |
9462866a | 829 | } |
830 | } | |
9462866a | 831 | |
d667707c | 832 | printf(" const Double_t sx[%d][%d]={\n", kNTB-1, kND); |
833 | for(Int_t iy=1; iy<kNTB; iy++){ | |
834 | printf(" {"); | |
835 | for(Int_t ix=1; ix<kND; ix++){ | |
836 | printf("%5.3e, ", hsx->GetBinContent(ix, iy)); | |
837 | } | |
838 | printf("%5.3e}", hsx->GetBinContent(kND, iy)); | |
839 | printf("%c\n", iy==(kNTB-1)?' ':','); | |
f2231ce8 | 840 | } |
d667707c | 841 | printf(" };\n"); |
9462866a | 842 | |
d667707c | 843 | printf(" const Double_t sy[%d][%d]={\n", kNTB-1, kND); |
844 | for(Int_t iy=1; iy<kNTB; iy++){ | |
845 | printf(" {"); | |
846 | for(Int_t ix=1; ix<kND; ix++){ | |
847 | printf("%5.3e, ", hsy->GetBinContent(ix, iy)); | |
848 | } | |
849 | printf("%5.3e}", hsy->GetBinContent(kND, iy)); | |
850 | printf("%c\n", iy==(kNTB-1)?' ':','); | |
f2231ce8 | 851 | } |
d667707c | 852 | printf(" };\n"); |
853 | ||
854 | return; | |
9462866a | 855 | } |
856 | ||
857 | //_______________________________________________________ | |
858 | void AliTRDclusterResolution::ProcessMean() | |
859 | { | |
860 | TObjArray *arr = (TObjArray*)fContainer->At(kMean); | |
861 | if(!arr){ | |
862 | AliWarning("Missing dy=f(x_d, d_wire) container"); | |
863 | return; | |
864 | } | |
865 | TGraphErrors *gm = new TGraphErrors(); | |
866 | TF1 f("f", "gaus", -.5, .5); | |
81a6494d | 867 | TF1 line("l", "[0]+[1]*x", -.15, .15); |
868 | Double_t d(0.), x(0.), dx, dy; | |
9462866a | 869 | TAxis *ax = 0x0; |
870 | TH1D *h1 = 0x0; | |
871 | TH2I *h2 = 0x0; | |
81a6494d | 872 | TH1 *hFrame=0x0; |
9462866a | 873 | |
9462866a | 874 | TObjArray *arrr = (TObjArray*)fResults->At(kMean); |
d667707c | 875 | TH2F *hdx = (TH2F*)arrr->At(2); |
876 | TH2F *hdy = (TH2F*)arrr->At(3); | |
9462866a | 877 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ |
878 | d = fAd->GetBinCenter(id); //[mm] | |
f2231ce8 | 879 | printf(" Doing d = %5.3f [mm]\n", d); |
9462866a | 880 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ |
881 | x = fAt->GetBinCenter(it); //[mm] | |
9462866a | 882 | Int_t idx = (id-1)*kNTB+it-1; |
883 | if(!(h2 = (TH2I*)arr->At(idx))) { | |
884 | AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d)); | |
885 | continue; | |
886 | } | |
f2231ce8 | 887 | |
888 | new(gm) TGraphErrors(); | |
889 | gm->SetMarkerStyle(7); | |
9462866a | 890 | ax = h2->GetXaxis(); |
891 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
892 | Double_t dydx = ax->GetBinCenter(ix); | |
893 | h1 = h2->ProjectionY("py", ix, ix); | |
81a6494d | 894 | if(h1->GetEntries() < 200) continue; |
9462866a | 895 | //Adjust(&f, h1); |
9462866a | 896 | h1->Fit(&f, "QN"); |
897 | ||
898 | // Fill <Dy> = f(dydx - h*dzdx) | |
f2231ce8 | 899 | Int_t ip = gm->GetN(); |
900 | gm->SetPoint(ip, dydx, f.GetParameter(1)); | |
901 | gm->SetPointError(ip, 0., f.GetParError(1)); | |
9462866a | 902 | ip++; |
903 | } | |
904 | if(gm->GetN()<4) continue; | |
905 | ||
906 | gm->Fit(&line, "QN"); | |
81a6494d | 907 | dx = line.GetParameter(1); |
908 | Double_t xs = line.GetParameter(0); | |
909 | dy = xs + fExB*dx; // xs = dy - tg(a_L)*dx | |
9462866a | 910 | hdx->SetBinContent(id, it, dx); |
81a6494d | 911 | hdy->SetBinContent(id, it, dy); |
f2231ce8 | 912 | |
913 | if(!fCanvas) continue; | |
914 | fCanvas->cd(); | |
81a6494d | 915 | if(!hFrame){ |
916 | hFrame=new TH1I("hFrame", "", 100, -.3, .3); | |
917 | hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1); | |
918 | hFrame->SetXTitle("tg#phi-htg#theta"); | |
919 | hFrame->SetYTitle("#Deltay[cm]"); | |
920 | hFrame->SetLineColor(1);hFrame->SetLineWidth(1); | |
921 | hFrame->Draw(); | |
922 | } else hFrame->Reset(); | |
923 | gm->Draw("pl"); line.Draw("same"); | |
f2231ce8 | 924 | fCanvas->Modified(); fCanvas->Update(); |
925 | if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_D%d_T%02d.gif", id, it)); | |
926 | else gSystem->Sleep(100); | |
81a6494d | 927 | printf(" xd=%4.1f[cm] dx=%5.3e[cm] dy=%5.3e[cm] xs=%5.3e[cm]\n", x, dx, dy, xs); |
9462866a | 928 | } |
929 | } | |
930 | ||
d667707c | 931 | // dump to stdout correction map |
932 | printf(" const Double_t dx[%d][%d]={\n", kNTB-1, kND); | |
933 | for(Int_t iy=1; iy<kNTB; iy++){ | |
934 | printf(" {"); | |
935 | for(Int_t ix=1; ix<kND; ix++){ | |
936 | printf("%+5.3e,", hdx->GetBinContent(ix, iy)); | |
937 | } | |
938 | printf("%+5.3e}", hdx->GetBinContent(kND, iy)); | |
939 | printf("%c\n", iy==(kNTB-1)?' ':','); | |
f2231ce8 | 940 | } |
941 | printf(" };\n"); | |
942 | ||
d667707c | 943 | // Collapse the z direction |
944 | TAxis *ay = hdx->GetYaxis(); | |
945 | // radial systematics | |
946 | TGraphErrors *g = (TGraphErrors*)arrr->At(0); | |
947 | for(Int_t iy = 1; iy<kNTB; iy++){ | |
948 | Double_t m=0., rms=0.; | |
949 | for(Int_t ix = 1; ix<=kND; ix++){ | |
950 | d = hdx->GetBinContent(ix, iy); | |
951 | m += d; rms += (d*d); | |
952 | } | |
953 | m /= Int_t(kND); rms = TMath::Sqrt(rms/Int_t(kND)-m*m); | |
954 | g->SetPoint(iy-1, ay->GetBinCenter(iy), 1.e4*m); | |
955 | g->SetPointError(iy-1, 0., 1.e4*rms); | |
956 | } | |
957 | ||
958 | // r-phi systematics | |
959 | g = (TGraphErrors*)arrr->At(1); | |
960 | for(Int_t iy = 1; iy<kNTB; iy++){ | |
961 | Double_t m=0., rms=0.; | |
962 | for(Int_t ix = 1; ix<=kND; ix++){ | |
963 | d = hdy->GetBinContent(ix, iy); | |
964 | m += d; rms += (d*d); | |
965 | } | |
966 | m /= Int_t(kND); rms = TMath::Sqrt(rms/Int_t(kND)-m*m); | |
967 | g->SetPoint(iy-1, ay->GetBinCenter(iy), 1.e4*m); | |
968 | g->SetPointError(iy-1, 0., 1.e4*rms); | |
969 | } | |
f2231ce8 | 970 | // delete gm; TODO memory leak ? |
9462866a | 971 | } |