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