]>
Commit | Line | Data |
---|---|---|
b2dc316d | 1 | #include "AliTRDclusterResolution.h" |
2 | #include "AliTRDtrackInfo/AliTRDclusterInfo.h" | |
5198d8c6 | 3 | #include "AliTRDgeometry.h" |
4 | #include "AliTRDpadPlane.h" | |
b2dc316d | 5 | |
6 | #include "AliLog.h" | |
7 | ||
8 | #include "TObjArray.h" | |
9 | #include "TAxis.h" | |
fc0946a7 | 10 | #include "TF1.h" |
11 | #include "TGraphErrors.h" | |
b2dc316d | 12 | #include "TH2I.h" |
13 | #include "TMath.h" | |
14 | ||
15 | ClassImp(AliTRDclusterResolution) | |
16 | ||
17 | ||
18 | //_______________________________________________________ | |
19 | AliTRDclusterResolution::AliTRDclusterResolution() | |
5198d8c6 | 20 | : AliTRDrecoTask("ClErrParam", "Cluster Error Parametrization") |
b2dc316d | 21 | ,fInfo(0x0) |
fc0946a7 | 22 | ,fResults(0x0) |
23 | ,fAt(0x0) | |
24 | ,fAd(0x0) | |
25 | ,fExB(0.) | |
b2dc316d | 26 | { |
fc0946a7 | 27 | fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15); |
28 | fAd = new TAxis(kND, 0., .25); | |
b2dc316d | 29 | } |
30 | ||
31 | //_______________________________________________________ | |
32 | AliTRDclusterResolution::~AliTRDclusterResolution() | |
33 | { | |
fc0946a7 | 34 | if(fAd) delete fAd; |
35 | if(fAt) delete fAt; | |
36 | if(fResults){ | |
37 | fResults->Delete(); | |
38 | delete fResults; | |
39 | } | |
b2dc316d | 40 | } |
41 | ||
42 | //_______________________________________________________ | |
43 | void AliTRDclusterResolution::ConnectInputData(Option_t *) | |
44 | { | |
45 | fInfo = dynamic_cast<TObjArray *>(GetInputData(0)); | |
46 | } | |
47 | ||
48 | //_______________________________________________________ | |
49 | void AliTRDclusterResolution::CreateOutputObjects() | |
50 | { | |
51 | OpenFile(0, "RECREATE"); | |
52 | fContainer = Histos(); | |
53 | } | |
54 | ||
55 | //_______________________________________________________ | |
fc0946a7 | 56 | void AliTRDclusterResolution::GetRefFigure(Int_t ifig) |
b2dc316d | 57 | { |
fc0946a7 | 58 | if(!fResults) return /*kFALSE*/; |
59 | ||
60 | TObjArray *arr = 0x0; | |
61 | TH2 *h2 = 0x0; | |
62 | TGraphErrors *gm(0x0), *gs(0x0); | |
63 | switch(ifig){ | |
64 | case kQRes: | |
65 | if(!(arr = (TObjArray*)fResults->At(kQRes))) break; | |
66 | if(!(gm = (TGraphErrors*)arr->At(0))) break; | |
67 | if(!(gs = (TGraphErrors*)arr->At(1))) break; | |
68 | gs->Draw("apl"); | |
69 | gs->GetHistogram()->SetXTitle("Log(Q) [a.u.]"); | |
70 | gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]"); | |
71 | return /*kTRUE*/; | |
72 | case kYRes: | |
73 | if(!(arr = (TObjArray*)fResults->At(kYRes))) break; | |
74 | if(!(gm = (TGraphErrors*)arr->At(0))) break; | |
75 | if(!(gs = (TGraphErrors*)arr->At(1))) break; | |
76 | gs->Draw("apl"); | |
77 | gs->GetHistogram()->SetXTitle("y [w]"); | |
78 | gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]"); | |
79 | return /*kTRUE*/; | |
80 | case kSXRes: | |
81 | if(!(h2 = (TH2F*)fResults->At(kSXRes))) break; | |
82 | h2->Draw("lego2"); | |
83 | return /*kTRUE*/; | |
84 | case kSYRes: | |
85 | if(!(h2 = (TH2F*)fResults->At(kSYRes))) break; | |
86 | h2->Draw("lego2"); | |
87 | return /*kTRUE*/; | |
88 | default: | |
89 | break; | |
90 | } | |
91 | ||
92 | return /*kFALSE*/; | |
b2dc316d | 93 | } |
94 | ||
95 | //_______________________________________________________ | |
96 | TObjArray* AliTRDclusterResolution::Histos() | |
97 | { | |
98 | if(fContainer) return fContainer; | |
5198d8c6 | 99 | fContainer = new TObjArray(3); |
b2dc316d | 100 | //fContainer->SetOwner(kTRUE); |
101 | ||
5198d8c6 | 102 | TH2I *h2 = 0x0; |
103 | TObjArray *arr = 0x0; | |
104 | ||
105 | fContainer->AddAt(h2 = new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), kQRes); | |
106 | h2->SetXTitle("log(q) [a.u.]"); | |
107 | h2->SetYTitle("#Delta y[cm]"); | |
108 | h2->SetZTitle("entries"); | |
109 | ||
110 | Double_t w = 0.; | |
111 | AliTRDgeometry geo; | |
112 | fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kYRes); | |
113 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ | |
114 | w = .5*geo.GetPadPlane(ily, 2)->GetWidthIPad(); | |
115 | arr->AddAt(h2 = new TH2I(Form("h_y%d", ily), Form("Ly[%d]", ily), 50, -w, w, 100, -.5, .5), ily); | |
116 | h2->SetXTitle("y_{w} [w]"); | |
117 | h2->SetYTitle("#Delta y[cm]"); | |
118 | h2->SetZTitle("entries"); | |
119 | } | |
120 | ||
121 | fContainer->AddAt(arr = new TObjArray(kN), kSXRes); | |
b2dc316d | 122 | Int_t ih = 0; |
fc0946a7 | 123 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ |
124 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ | |
5198d8c6 | 125 | 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)), 30, -.15, .15, 100, -.5, .5), ih++); |
126 | h2->SetXTitle("tg#phi"); | |
127 | h2->SetYTitle("#Delta y[cm]"); | |
128 | h2->SetZTitle("entries"); | |
b2dc316d | 129 | } |
130 | } | |
b2dc316d | 131 | return fContainer; |
132 | } | |
133 | ||
134 | //_______________________________________________________ | |
135 | void AliTRDclusterResolution::Exec(Option_t *) | |
136 | { | |
87b166d3 | 137 | Int_t det, t; |
fc0946a7 | 138 | Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3]; |
b2dc316d | 139 | TH2I *h2 = 0x0; |
5198d8c6 | 140 | |
141 | TObjArray *arr0 = (TObjArray*)fContainer->At(kYRes); | |
142 | TObjArray *arr1 = (TObjArray*)fContainer->At(kSXRes); | |
143 | ||
b2dc316d | 144 | const AliTRDclusterInfo *cli = 0x0; |
145 | TIterator *iter=fInfo->MakeIterator(); | |
146 | while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){ | |
147 | dy = cli->GetResolution(); | |
fc0946a7 | 148 | Int_t it = fAt->FindBin(cli->GetDriftLength()); |
149 | if(it==0 || it == fAt->GetNbins()+1){ | |
b2dc316d | 150 | AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength())); |
151 | continue; | |
152 | } | |
fc0946a7 | 153 | Int_t id = fAd->FindBin(cli->GetAnisochronity()); |
154 | if(id==0 || id == fAd->GetNbins()+1){ | |
b2dc316d | 155 | AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity())); |
156 | continue; | |
157 | } | |
5198d8c6 | 158 | if(!(h2 = (TH2I*)arr1->At((id-1)*kNTB+it-1))){ |
b2dc316d | 159 | AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", (id-1)*kNTB+it-1, id, it, cli->GetDriftLength(), cli->GetAnisochronity())); |
160 | continue; | |
161 | } | |
162 | ||
163 | cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]); | |
fc0946a7 | 164 | h2->Fill(dydx, cov[0]!=0. ? dy/TMath::Sqrt(cov[0]) : dy); |
b2dc316d | 165 | |
fc0946a7 | 166 | // resolution as a function of: |
167 | // - cluster charge and | |
168 | // - y displacement | |
b2dc316d | 169 | // only for phi equal exB |
fc0946a7 | 170 | if(TMath::Abs(dydx-fExB)<.01){ |
171 | cli->GetCluster(det, x, y, z, q, t, covcl); | |
5198d8c6 | 172 | h2 = (TH2I*)fContainer->At(kQRes); |
b2dc316d | 173 | h2->Fill(TMath::Log(q), dy); |
fc0946a7 | 174 | |
5198d8c6 | 175 | h2 = (TH2I*)arr0->At(AliTRDgeometry::GetLayer(det)); |
fc0946a7 | 176 | h2->Fill(cli->GetYDisplacement(), dy); |
b2dc316d | 177 | } |
178 | } | |
179 | PostData(0, fContainer); | |
180 | } | |
181 | ||
182 | ||
183 | //_______________________________________________________ | |
184 | Bool_t AliTRDclusterResolution::PostProcess() | |
185 | { | |
fc0946a7 | 186 | if(!fContainer) return kFALSE; |
187 | ||
5198d8c6 | 188 | TObjArray *arr = 0x0; |
fc0946a7 | 189 | TH2 *h2 = 0x0; |
190 | if(!fResults){ | |
191 | TGraphErrors *g = 0x0; | |
192 | fResults = new TObjArray(4); | |
193 | fResults->SetOwner(); | |
fc0946a7 | 194 | fResults->AddAt(arr = new TObjArray(2), kQRes); |
195 | arr->SetOwner(); | |
196 | arr->AddAt(g = new TGraphErrors(), 0); | |
197 | g->SetLineColor(kBlue); g->SetMarkerColor(kBlue); | |
198 | g->SetMarkerStyle(7); | |
199 | arr->AddAt(g = new TGraphErrors(), 1); | |
200 | g->SetLineColor(kRed); g->SetMarkerColor(kRed); | |
201 | g->SetMarkerStyle(23); | |
202 | ||
203 | fResults->AddAt(arr = new TObjArray(2), kYRes); | |
204 | arr->SetOwner(); | |
205 | arr->AddAt(g = new TGraphErrors(), 0); | |
206 | g->SetLineColor(kBlue); g->SetMarkerColor(kBlue); | |
207 | g->SetMarkerStyle(7); | |
208 | arr->AddAt(g = new TGraphErrors(), 1); | |
209 | g->SetLineColor(kRed); g->SetMarkerColor(kRed); | |
210 | g->SetMarkerStyle(23); | |
211 | ||
212 | fResults->AddAt(h2 = new TH2F("hSX", "", | |
213 | fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(), | |
214 | fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()), kSXRes); | |
215 | h2->SetXTitle("d [mm]"); | |
216 | h2->SetYTitle("x [mm]"); | |
217 | h2->SetZTitle("#sigma_{x}^{2} [mm^{2}]"); | |
218 | fResults->AddAt(h2 = (TH2F*)h2->Clone("hSY"), kSYRes); | |
219 | h2->SetZTitle("#sigma_{y}^{2} [mm^{2}]"); | |
220 | } else { | |
221 | TObject *o = 0x0; | |
222 | TIterator *iter=fResults->MakeIterator(); | |
5198d8c6 | 223 | while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point |
fc0946a7 | 224 | } |
225 | ||
226 | TF1 f("f", "gaus", -.5, .5); | |
227 | TF1 sig("sig", "pol1", 0., .0225); | |
228 | ||
229 | TAxis *ax = 0x0; | |
230 | TH1D *h1 = 0x0; | |
231 | ||
232 | // process resolution dependency on charge | |
5198d8c6 | 233 | if((h2 = (TH2I*)fContainer->At(kQRes))) { |
fc0946a7 | 234 | TObjArray *arr = (TObjArray*)fResults->At(kQRes); |
235 | TGraphErrors *gqm = (TGraphErrors*)arr->At(0); | |
236 | TGraphErrors *gqs = (TGraphErrors*)arr->At(1); | |
237 | ax = h2->GetXaxis(); | |
238 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
239 | Float_t logq = ax->GetBinCenter(ix); | |
240 | h1 = h2->ProjectionY("py", ix, ix); | |
241 | if(h1->GetEntries() < 50) continue; | |
242 | Adjust(&f, h1); | |
243 | h1->Fit(&f, "Q"); | |
244 | ||
245 | // Fill sy^2 = f(log(q)) | |
246 | Int_t ip = gqm->GetN(); | |
247 | gqm->SetPoint(ip, logq, 10.*f.GetParameter(1)); | |
248 | gqm->SetPointError(ip, 0., 10.*f.GetParError(1)); | |
249 | gqs->SetPoint(ip, logq, 1.e2*f.GetParameter(2)*f.GetParameter(2)); | |
250 | gqs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2)); | |
251 | } | |
252 | } else AliWarning("Missing dy=f(Q) histo"); | |
253 | ||
254 | // process resolution dependency on y displacement | |
5198d8c6 | 255 | if((arr = (TObjArray*)fContainer->At(kYRes))) { |
256 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ | |
257 | if(!(h2 = (TH2I*)arr->At(ily))) continue; | |
258 | TObjArray *arrg = (TObjArray*)fResults->At(kYRes); | |
259 | TGraphErrors *gym = (TGraphErrors*)arrg->At(0); | |
260 | TGraphErrors *gys = (TGraphErrors*)arrg->At(1); | |
fc0946a7 | 261 | ax = h2->GetXaxis(); |
262 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
5198d8c6 | 263 | Float_t yd = ax->GetBinCenter(ix); |
fc0946a7 | 264 | h1 = h2->ProjectionY("py", ix, ix); |
265 | if(h1->GetEntries() < 50) continue; | |
266 | Adjust(&f, h1); | |
5198d8c6 | 267 | h1->Fit(&f, "Q"); |
268 | ||
269 | // Fill sy^2 = f(log(q)) | |
270 | Int_t ip = gym->GetN(); | |
271 | gym->SetPoint(ip, yd, 10.*f.GetParameter(1)); | |
272 | gym->SetPointError(ip, 0., 10.*f.GetParError(1)); | |
273 | gys->SetPoint(ip, yd, 1.e2*f.GetParameter(2)*f.GetParameter(2)); | |
274 | gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2)); | |
275 | } | |
276 | } | |
277 | } else AliWarning("Missing dy=f(y_d) container"); | |
fc0946a7 | 278 | |
5198d8c6 | 279 | // process resolution dependency on drift legth and drift cell width |
280 | TGraphErrors *gm = new TGraphErrors(), *gs = new TGraphErrors(); | |
281 | Float_t d(0.), x(0.); | |
282 | TH2F *hsx = (TH2F*)fResults->At(kSXRes); | |
283 | TH2F *hsy = (TH2F*)fResults->At(kSYRes); | |
284 | if((arr = (TObjArray*)fContainer->At(kSXRes))) { | |
285 | for(Int_t id=1; id<=fAd->GetNbins(); id++){ | |
286 | d = fAd->GetBinCenter(id); //[mm] | |
287 | printf("Doing d=%f[mm]\n", d); | |
288 | for(Int_t it=1; it<=fAt->GetNbins(); it++){ | |
289 | x = fAt->GetBinCenter(it); //[mm] | |
290 | printf("Doing x=%f[mm]\n", x); | |
291 | Int_t idx = (id-1)*kNTB+it-1; | |
292 | if(!(h2 = (TH2I*)arr->At(idx))) { | |
293 | AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d)); | |
294 | continue; | |
295 | } | |
296 | ||
297 | Int_t ip = 0; | |
298 | ax = h2->GetXaxis(); | |
299 | for(Int_t ix=1; ix<=ax->GetNbins(); ix++){ | |
300 | Float_t dydx = ax->GetBinCenter(ix) - fExB; | |
301 | if(dydx<0.) continue; | |
302 | h1 = h2->ProjectionY("py", ix, ix); | |
303 | if(h1->GetEntries() < 50) continue; | |
304 | Adjust(&f, h1); | |
305 | h1->Fit(&f, "Q", "", -.2, .2); | |
306 | ||
307 | // Fill sy^2 = f(tg_phi^2) | |
308 | gm->SetPoint(ip, dydx, 10.*f.GetParameter(1)); | |
309 | gm->SetPointError(ip, 0., 10.*f.GetParError(1)); | |
310 | gs->SetPoint(ip, dydx*dydx, 1.e2*f.GetParameter(2)*f.GetParameter(2)); | |
311 | gs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2)); | |
312 | ip++; | |
313 | } | |
314 | for(;ip<gm->GetN(); ip++){ | |
315 | gm->RemovePoint(ip);gs->RemovePoint(ip); | |
316 | } | |
317 | gs->Fit(&sig, "Q"); | |
318 | hsx->SetBinContent(id, it, sig.GetParameter(1)); | |
319 | hsx->SetBinError(id, it, sig.GetParError(1)); | |
320 | hsy->SetBinContent(id, it, sig.GetParameter(0)); | |
321 | hsy->SetBinError(id, it, sig.GetParError(0)); | |
fc0946a7 | 322 | } |
fc0946a7 | 323 | } |
5198d8c6 | 324 | } else AliWarning("Missing dy=f(x_d, d_wire) container"); |
fc0946a7 | 325 | delete gm; delete gs; |
326 | ||
327 | return kTRUE; | |
b2dc316d | 328 | } |
329 | ||
330 |