]>
Commit | Line | Data |
---|---|---|
77203477 | 1 | /************************************************************************** |
d2381af5 | 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 | **************************************************************************/ | |
77203477 | 15 | |
16 | /* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */ | |
17 | ||
18 | //////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // Reconstruction QA // | |
21 | // // | |
22 | // Authors: // | |
23 | // Markus Fasel <M.Fasel@gsi.de> // | |
24 | // // | |
25 | //////////////////////////////////////////////////////////////////////////// | |
26 | ||
aaf47b30 | 27 | #include <cstring> |
28 | ||
29 | ||
77203477 | 30 | #include <TObjArray.h> |
7102d1b1 | 31 | #include <TH2.h> |
32 | #include <TH1.h> | |
33 | #include <TF1.h> | |
77203477 | 34 | #include <TProfile.h> |
7102d1b1 | 35 | #include <TGraphErrors.h> |
77203477 | 36 | #include <TMath.h> |
aaf47b30 | 37 | #include "TTreeStream.h" |
38 | #include "TGeoManager.h" | |
77203477 | 39 | |
40 | #include "AliAnalysisManager.h" | |
77203477 | 41 | #include "AliTrackReference.h" |
aaf47b30 | 42 | #include "AliTrackPointArray.h" |
43 | #include "AliCDBManager.h" | |
44 | ||
9605ce80 | 45 | #include "AliTRDSimParam.h" |
46 | #include "AliTRDgeometry.h" | |
47 | #include "AliTRDpadPlane.h" | |
aaf47b30 | 48 | #include "AliTRDcluster.h" |
49 | #include "AliTRDseedV1.h" | |
50 | #include "AliTRDtrackV1.h" | |
51 | #include "AliTRDtrackerV1.h" | |
52 | #include "AliTRDReconstructor.h" | |
53 | #include "AliTRDrecoParam.h" | |
77203477 | 54 | |
55 | #include "AliTRDtrackInfo/AliTRDtrackInfo.h" | |
56 | #include "AliTRDtrackingResolution.h" | |
57 | ||
77203477 | 58 | ClassImp(AliTRDtrackingResolution) |
59 | ||
60 | //________________________________________________________ | |
3d86166d | 61 | AliTRDtrackingResolution::AliTRDtrackingResolution() |
62 | :AliTRDrecoTask("Resolution", "Tracking Resolution") | |
aaf47b30 | 63 | ,fReconstructor(0x0) |
9605ce80 | 64 | ,fGeo(0x0) |
77203477 | 65 | { |
aaf47b30 | 66 | fReconstructor = new AliTRDReconstructor(); |
67 | fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam()); | |
9605ce80 | 68 | fGeo = new AliTRDgeometry(); |
77203477 | 69 | } |
70 | ||
ed383798 | 71 | //________________________________________________________ |
72 | AliTRDtrackingResolution::~AliTRDtrackingResolution() | |
73 | { | |
9605ce80 | 74 | delete fGeo; |
ed383798 | 75 | delete fReconstructor; |
2c0cf367 | 76 | if(gGeoManager) delete gGeoManager; |
ed383798 | 77 | } |
78 | ||
77203477 | 79 | |
80 | //________________________________________________________ | |
81 | void AliTRDtrackingResolution::CreateOutputObjects() | |
82 | { | |
83 | // spatial resolution | |
77203477 | 84 | OpenFile(0, "RECREATE"); |
39779ce6 | 85 | |
3d86166d | 86 | fContainer = Histos(); |
874acced | 87 | |
88 | // cluster to tracklet residuals [2] | |
765bd0ab | 89 | fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -21., 21., 100, -.5, .5), kClusterYResidual); |
874acced | 90 | // tracklet to Riemann fit residuals [2] |
765bd0ab | 91 | fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual); |
92 | fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual); | |
93 | fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual); | |
94 | fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual); | |
39779ce6 | 95 | |
96 | // Resolution histos | |
4b8f8a35 | 97 | if(HasMCdata()){ |
d2381af5 | 98 | // tracklet resolution [0] |
765bd0ab | 99 | fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 21, -21., 21., 100, -.5, .5), kTrackletYResolution); |
d2381af5 | 100 | // tracklet angular resolution [1] |
765bd0ab | 101 | fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackletAngleResolution); |
874acced | 102 | |
103 | // Riemann track resolution [y, z, angular] | |
765bd0ab | 104 | fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution); |
105 | fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution); | |
106 | fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution); | |
874acced | 107 | |
108 | // Kalman track resolution [y, z, angular] | |
3d86166d | 109 | fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution); |
110 | fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution); | |
111 | fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution); | |
d2381af5 | 112 | } |
d85cd79c | 113 | |
114 | // CREATE GRAPHS for DISPLAY | |
115 | ||
116 | // define iterator over graphs | |
117 | Int_t jgraph = (Int_t)kGraphStart; | |
118 | TH2I *h2 = (TH2I *)(fContainer->At(kClusterYResidual)); | |
119 | // clusters tracklet residuals (mean-phi) | |
765bd0ab | 120 | TH1 *h = new TH1I("h", "", 100, -40., 40.); |
121 | h->GetXaxis()->SetTitle("#Phi [deg]"); | |
122 | h->GetYaxis()->SetTitle("Clusters Residuals : #sigma/#mu [mm]"); | |
123 | h->GetYaxis()->SetRangeUser(-.05, 1.); | |
124 | fContainer->AddAt(h, jgraph++); | |
125 | ||
d85cd79c | 126 | TGraphErrors *g = new TGraphErrors(h2->GetNbinsX()); |
127 | g->SetLineColor(kGreen); | |
128 | g->SetMarkerStyle(22); | |
129 | g->SetMarkerColor(kGreen); | |
130 | g->SetNameTitle("clm", "Residuals Clusters-Tracklet Mean"); | |
131 | fContainer->AddAt(g, jgraph++); | |
132 | ||
133 | // clusters tracklet residuals (sigma-phi) | |
134 | g = new TGraphErrors(h2->GetNbinsX()); | |
135 | g->SetLineColor(kRed); | |
136 | g->SetMarkerStyle(23); | |
137 | g->SetMarkerColor(kRed); | |
138 | g->SetNameTitle("cls", "Residuals Clusters-Tracklet Sigma"); | |
139 | fContainer->AddAt(g, jgraph++); | |
140 | ||
141 | if(HasMCdata()){ | |
142 | // tracklet y resolution | |
143 | h2 = (TH2I*)fContainer->At(kTrackletYResolution); | |
765bd0ab | 144 | h = new TH1I("h", "", 100, -40., 40.); |
145 | h->GetXaxis()->SetTitle("#Phi [deg]"); | |
146 | h->GetYaxis()->SetTitle("Tracklet Resolution : #sigma/#mu [mm]"); | |
147 | h->GetYaxis()->SetRangeUser(-.05, 1.); | |
148 | fContainer->AddAt(h, jgraph++); | |
149 | ||
d85cd79c | 150 | g = new TGraphErrors(h2->GetNbinsX()); |
151 | g->SetLineColor(kGreen); | |
152 | g->SetMarkerStyle(22); | |
153 | g->SetMarkerColor(kGreen); | |
154 | g->SetNameTitle("trkltym", "Resolution Tracklet Y Mean"); | |
155 | fContainer->AddAt(g, jgraph++); | |
156 | g = new TGraphErrors(h2->GetNbinsX()); | |
157 | g->SetLineColor(kRed); | |
158 | g->SetMarkerStyle(22); | |
159 | g->SetMarkerColor(kRed); | |
160 | g->SetNameTitle("trkltys", "Resolution Tracklet Y Sigma"); | |
161 | fContainer->AddAt(g, jgraph++); | |
162 | ||
163 | // tracklet phi resolution | |
164 | h2 = (TH2I*)fContainer->At(kTrackletAngleResolution); | |
765bd0ab | 165 | h = new TH1I("h", "", 100, -40., 40.); |
166 | h->GetXaxis()->SetTitle("#Phi [deg]"); | |
167 | h->GetYaxis()->SetTitle("Tracklet Angular Resolution : #sigma/#mu [deg]"); | |
28efdace | 168 | h->GetYaxis()->SetRangeUser(-.05, .2); |
765bd0ab | 169 | fContainer->AddAt(h, jgraph++); |
170 | ||
d85cd79c | 171 | g = new TGraphErrors(h2->GetNbinsX()); |
172 | g->SetLineColor(kGreen); | |
173 | g->SetMarkerStyle(22); | |
174 | g->SetMarkerColor(kGreen); | |
175 | g->SetNameTitle("trkltam", "Resolution Tracklet Y Mean"); | |
176 | fContainer->AddAt(g, jgraph++); | |
177 | g = new TGraphErrors(h2->GetNbinsX()); | |
178 | g->SetLineColor(kRed); | |
179 | g->SetMarkerStyle(22); | |
180 | g->SetMarkerColor(kRed); | |
181 | g->SetNameTitle("trkltas", "Angle Resolution Tracklet Sigma"); | |
182 | fContainer->AddAt(g, jgraph++); | |
183 | } | |
77203477 | 184 | } |
185 | ||
186 | //________________________________________________________ | |
7102d1b1 | 187 | void AliTRDtrackingResolution::Exec(Option_t *) |
188 | { | |
77203477 | 189 | // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire}) |
190 | // angular Resolution: res = Tracklet angle - TrackRef Angle | |
7102d1b1 | 191 | |
39779ce6 | 192 | Int_t nTrackInfos = fTracks->GetEntriesFast(); |
874acced | 193 | if(fDebugLevel>=2) printf("Number of Histograms: %d\n", Histos()->GetEntries()); |
aaf47b30 | 194 | |
9605ce80 | 195 | Int_t pdg; |
765bd0ab | 196 | Double_t p, dy, dphi, dymc, dzmc, dphimc; |
9605ce80 | 197 | Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer |
765bd0ab | 198 | Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map |
199 | ||
9605ce80 | 200 | AliTRDpadPlane *pp = 0x0; |
39779ce6 | 201 | AliTrackPoint tr[kNLayers], tk[kNLayers]; |
765bd0ab | 202 | AliExternalTrackParam *fOp = 0x0; |
aaf47b30 | 203 | AliTRDtrackV1 *fTrack = 0x0; |
77203477 | 204 | AliTRDtrackInfo *fInfo = 0x0; |
205 | if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos); | |
206 | for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){ | |
77203477 | 207 | // check if ESD and MC-Information are available |
39779ce6 | 208 | if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue; |
aaf47b30 | 209 | if(!(fTrack = fInfo->GetTRDtrack())) continue; |
765bd0ab | 210 | if(!(fOp = fInfo->GetOuterParam())) continue; |
9605ce80 | 211 | pdg = fInfo->GetPDG(); |
212 | ||
7102d1b1 | 213 | if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs()); |
214 | ||
765bd0ab | 215 | p = fOp->P(); |
39779ce6 | 216 | Int_t npts = 0; |
765bd0ab | 217 | memset(fP, 0, kNLayers*sizeof(Float_t)); |
9605ce80 | 218 | memset(fX, 0, kNLayers*sizeof(Float_t)); |
765bd0ab | 219 | memset(fY, 0, kNLayers*sizeof(Float_t)); |
220 | memset(fZ, 0, kNLayers*sizeof(Float_t)); | |
221 | memset(fPhi, 0, kNLayers*sizeof(Float_t)); | |
222 | memset(fTheta, 0, kNLayers*sizeof(Float_t)); | |
223 | memset(fLayerMap, 0, kNLayers*sizeof(Bool_t)); | |
224 | memset(fMCMap, 0, kNLayers*sizeof(Bool_t)); | |
aaf47b30 | 225 | AliTRDseedV1 *fTracklet = 0x0; |
77203477 | 226 | for(Int_t iplane = 0; iplane < kNLayers; iplane++){ |
aaf47b30 | 227 | if(!(fTracklet = fTrack->GetTracklet(iplane))) continue; |
228 | if(!fTracklet->IsOK()) continue; | |
229 | ||
39779ce6 | 230 | // Book point arrays |
765bd0ab | 231 | fLayerMap[iplane] = kTRUE; |
39779ce6 | 232 | tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.); |
233 | tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0)); | |
234 | npts++; | |
235 | ||
7102d1b1 | 236 | if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN()); |
237 | ||
765bd0ab | 238 | // define reference values |
239 | fP[iplane] = p; | |
9605ce80 | 240 | fX[iplane] = fTracklet->GetX0(); |
765bd0ab | 241 | fY[iplane] = fTracklet->GetYref(0); |
242 | fZ[iplane] = fTracklet->GetZref(0); | |
243 | fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1)); | |
244 | fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1)); | |
245 | ||
7102d1b1 | 246 | |
39779ce6 | 247 | // RESOLUTION (compare to MC) |
4b8f8a35 | 248 | if(HasMCdata()){ |
d2381af5 | 249 | if(fInfo->GetNTrackRefs() >= 2){ |
765bd0ab | 250 | Double_t pmc, ymc, zmc, phiMC, thetaMC; |
251 | if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){ | |
252 | fMCMap[iplane] = kTRUE; | |
253 | fP[iplane] = pmc; | |
254 | fY[iplane] = ymc; | |
255 | fZ[iplane] = zmc; | |
256 | fPhi[iplane] = phiMC; | |
257 | fTheta[iplane] = thetaMC; | |
258 | } | |
d2381af5 | 259 | } |
260 | } | |
765bd0ab | 261 | Float_t phi = fPhi[iplane]*TMath::RadToDeg(); |
9605ce80 | 262 | //Float_t theta = fTheta[iplane]*TMath::RadToDeg(); |
aaf47b30 | 263 | |
39779ce6 | 264 | // Do clusters residuals |
d2381af5 | 265 | if(!fTracklet->Fit(kFALSE)) continue; |
aaf47b30 | 266 | AliTRDcluster *c = 0x0; |
267 | for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){ | |
268 | if(!(c = fTracklet->GetClusters(ic))) continue; | |
9605ce80 | 269 | |
39779ce6 | 270 | dy = fTracklet->GetYat(c->GetX()) - c->GetY(); |
765bd0ab | 271 | ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy); |
9605ce80 | 272 | |
3d86166d | 273 | if(fDebugLevel>=2){ |
39779ce6 | 274 | Float_t q = c->GetQ(); |
9605ce80 | 275 | // Get z-position with respect to anode wire |
276 | AliTRDSimParam *simParam = AliTRDSimParam::Instance(); | |
277 | Int_t det = c->GetDetector(); | |
278 | Float_t x = c->GetX(); | |
279 | Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]); | |
280 | Int_t stack = fGeo->GetStack(det); | |
281 | pp = fGeo->GetPadPlane(iplane, stack); | |
282 | Float_t row0 = pp->GetRow0(); | |
283 | Float_t d = row0 - z + simParam->GetAnodeWireOffset(); | |
284 | d -= ((Int_t)(2 * d)) / 2.0; | |
285 | //if (d > 0.25) d = 0.5 - d; | |
286 | ||
765bd0ab | 287 | (*fDebugStream) << "ResidualClusters" |
288 | << "ly=" << iplane | |
9605ce80 | 289 | << "stk=" << stack |
3c3d9ff1 | 290 | << "pdg=" << pdg |
9605ce80 | 291 | << "phi=" << fPhi[iplane] |
292 | << "tht=" << fTheta[iplane] | |
765bd0ab | 293 | << "q=" << q |
9605ce80 | 294 | << "x=" << x |
295 | << "z=" << z | |
296 | << "d=" << d | |
765bd0ab | 297 | << "dy=" << dy |
39779ce6 | 298 | << "\n"; |
299 | } | |
aaf47b30 | 300 | } |
9605ce80 | 301 | pp = 0x0; |
aaf47b30 | 302 | } |
39779ce6 | 303 | |
d2381af5 | 304 | |
39779ce6 | 305 | // this protection we might drop TODO |
306 | if(fTrack->GetNumberOfTracklets() < 6) continue; | |
307 | ||
308 | AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr); | |
765bd0ab | 309 | Int_t iref = 0; |
310 | for(Int_t ip=0; ip<kNLayers; ip++){ | |
311 | if(!fLayerMap[ip]) continue; | |
312 | fTracklet = fTrack->GetTracklet(ip); | |
313 | // recalculate fit based on the new tilt correction | |
314 | fTracklet->Fit(); | |
874acced | 315 | |
765bd0ab | 316 | dy = fTracklet->GetYfit(0) - tr[iref].GetY(); |
317 | ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy); | |
318 | ||
319 | dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1); | |
320 | ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi); | |
321 | ||
322 | if(HasMCdata()){ | |
323 | dymc = fY[ip] - tr[iref].GetY(); | |
324 | ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc); | |
325 | ||
326 | dzmc = fZ[ip] - tr[iref].GetZ(); | |
327 | ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc); | |
328 | ||
329 | dphimc = fPhi[ip] - fTracklet->GetYfit(1); | |
330 | ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc); | |
331 | } | |
874acced | 332 | |
765bd0ab | 333 | iref++; |
874acced | 334 | |
3d86166d | 335 | if(fDebugLevel>=2){ |
765bd0ab | 336 | (*fDebugStream) << "RiemannTrack" |
337 | << "ly=" << ip | |
338 | << "mc=" << fMCMap[ip] | |
339 | << "p=" << fP[ip] | |
340 | << "phi=" << fPhi[ip] | |
341 | << "tht=" << fTheta[ip] | |
342 | << "dy=" << dy | |
343 | << "dphi=" << dphi | |
344 | << "dymc=" << dymc | |
345 | << "dzmc=" << dzmc | |
346 | << "dphimc="<< dphimc | |
aaf47b30 | 347 | << "\n"; |
348 | } | |
349 | } | |
350 | ||
2c0cf367 | 351 | // if(!gGeoManager) TGeoManager::Import("geometry.root"); |
aaf47b30 | 352 | // AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr); |
353 | // for(Int_t ip=0; ip<nc; ip++){ | |
354 | // dy = cl[ip].GetY() - tr[ip].GetY(); | |
3d86166d | 355 | // ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy); |
aaf47b30 | 356 | // dz = cl[ip].GetZ() - tr[ip].GetZ(); |
357 | // if(fDebugLevel>=1){ | |
765bd0ab | 358 | // (*fDebugStream) << "KalmanTrack" |
aaf47b30 | 359 | // << "dy=" << dy |
360 | // << "dz=" << dz | |
361 | // /* << "phi=" << phi | |
362 | // << "theta=" << theta | |
363 | // << "dphi=" << dphi*/ | |
364 | // << "\n"; | |
365 | // } | |
366 | // } | |
39779ce6 | 367 | |
368 | ||
77203477 | 369 | } |
3d86166d | 370 | PostData(0, fContainer); |
77203477 | 371 | } |
372 | ||
d85cd79c | 373 | //________________________________________________________ |
28efdace | 374 | void AliTRDtrackingResolution::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t */*opt*/) |
d85cd79c | 375 | { |
376 | switch(ifig){ | |
377 | case 0: | |
765bd0ab | 378 | first = (Int_t)kGraphStart; last = first+3; |
d85cd79c | 379 | break; |
380 | case 1: | |
765bd0ab | 381 | first = (Int_t)kGraphStart+3; last = first+3; |
d85cd79c | 382 | break; |
383 | case 2: | |
765bd0ab | 384 | first = (Int_t)kGraphStart+6; last = first+3; |
d85cd79c | 385 | break; |
386 | default: | |
387 | first = (Int_t)kGraphStart; last = first; | |
388 | break; | |
389 | } | |
390 | } | |
391 | ||
39779ce6 | 392 | |
393 | //________________________________________________________ | |
765bd0ab | 394 | Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &p, Double_t &ymc, Double_t &zmc, Double_t &phi, Double_t &theta) |
39779ce6 | 395 | { |
396 | ||
397 | AliTrackReference *fTrackRefs[2] = {0x0, 0x0}, *tempTrackRef = 0x0; | |
398 | ||
399 | // check for 2 track ref where the radial position has a distance less than 3.7mm | |
400 | Int_t nFound = 0; | |
401 | for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){ | |
402 | if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue; | |
403 | if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX()); | |
404 | if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue; | |
405 | fTrackRefs[nFound++] = tempTrackRef; | |
406 | if(nFound == 2) break; | |
407 | } | |
408 | if(nFound < 2){ | |
409 | if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.); | |
410 | return kFALSE; | |
411 | } | |
412 | // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation | |
413 | ||
414 | ||
415 | // RESOLUTION | |
416 | Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX(); | |
417 | if(dx <= 0.){ | |
418 | if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX()); | |
419 | return kFALSE; | |
420 | } | |
421 | Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx; | |
422 | Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx; | |
423 | Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0(); | |
765bd0ab | 424 | ymc = fTrackRefs[1]->LocalY() - dydx*dx0; |
425 | zmc = fTrackRefs[1]->Z() - dzdx*dx0; | |
39779ce6 | 426 | |
427 | // recalculate tracklet based on the MC info | |
84b3ecda | 428 | tracklet->SetZref(0, zmc); |
765bd0ab | 429 | tracklet->SetZref(1, -dzdx); // TODO |
39779ce6 | 430 | tracklet->Fit(); |
431 | Double_t dy = tracklet->GetYfit(0) - ymc; | |
432 | Double_t dz = tracklet->GetZfit(0) - zmc; | |
433 | ||
434 | //res_y *= 100; // in mm | |
765bd0ab | 435 | p = fTrackRefs[0]->P(); |
39779ce6 | 436 | |
437 | phi = TMath::ATan(dydx); | |
765bd0ab | 438 | theta = TMath::ATan(dzdx); |
39779ce6 | 439 | Double_t dphi = TMath::ATan(tracklet->GetYfit(1)) - phi; |
440 | if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi); | |
441 | ||
442 | // Fill Histograms | |
443 | if(TMath::Abs(dx-3.7)<1.E-3){ | |
3d86166d | 444 | ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(phi*TMath::RadToDeg(), dy); |
445 | ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg()); | |
39779ce6 | 446 | } |
447 | // Fill Debug Tree | |
3d86166d | 448 | if(fDebugLevel>=2){ |
39779ce6 | 449 | Int_t iplane = tracklet->GetPlane(); |
9605ce80 | 450 | Int_t pdg = fInfo->GetPDG(); |
765bd0ab | 451 | (*fDebugStream) << "ResolutionTrklt" |
452 | << "ly=" << iplane | |
9605ce80 | 453 | << "pdg=" << pdg |
765bd0ab | 454 | << "p=" << p |
455 | << "phi=" << phi | |
456 | << "tht=" << theta | |
457 | << "ymc=" << ymc | |
458 | << "zmc=" << zmc | |
39779ce6 | 459 | << "dx=" << dx |
460 | << "dy=" << dy | |
461 | << "dz=" << dz | |
39779ce6 | 462 | << "dphi=" << dphi |
463 | << "\n"; | |
464 | } | |
465 | ||
466 | return kTRUE; | |
467 | } | |
468 | ||
77203477 | 469 | //________________________________________________________ |
d85cd79c | 470 | Bool_t AliTRDtrackingResolution::PostProcess() |
7102d1b1 | 471 | { |
d85cd79c | 472 | //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0)); |
765bd0ab | 473 | fNRefFigures = 0; |
d85cd79c | 474 | if (!fContainer) { |
475 | Printf("ERROR: list not available"); | |
476 | return kFALSE; | |
3d86166d | 477 | } |
7102d1b1 | 478 | |
d2381af5 | 479 | TH2I *h2 = 0x0; |
480 | TH1D *h = 0x0; | |
d85cd79c | 481 | TGraphErrors *gm = 0x0, *gs = 0x0; |
d2381af5 | 482 | TF1 f("f1", "gaus", -.5, .5); |
874acced | 483 | // define iterator over graphs |
484 | Int_t jgraph = (Int_t)kGraphStart; | |
485 | ||
486 | //PROCESS RESIDUAL DISTRIBUTIONS | |
487 | ||
488 | // Clusters residuals | |
3c3d9ff1 | 489 | // define model |
490 | TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5); | |
3d86166d | 491 | h2 = (TH2I *)(fContainer->At(kClusterYResidual)); |
765bd0ab | 492 | jgraph++; //skip the frame histo |
d85cd79c | 493 | gm = (TGraphErrors*)fContainer->At(jgraph++); |
494 | gs = (TGraphErrors*)fContainer->At(jgraph++); | |
874acced | 495 | for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){ |
496 | Double_t phi = h2->GetXaxis()->GetBinCenter(ibin); | |
497 | Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2; | |
498 | h = h2->ProjectionY("py", ibin, ibin); | |
3c3d9ff1 | 499 | Fit(h, &fc); |
500 | gm->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(1)); | |
501 | gm->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(1)); | |
502 | gs->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(2)); | |
503 | gs->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(2)); | |
874acced | 504 | } |
765bd0ab | 505 | fNRefFigures++; |
d2381af5 | 506 | |
874acced | 507 | |
508 | //PROCESS RESOLUTION DISTRIBUTIONS | |
509 | if(HasMCdata()){ | |
510 | // tracklet y resolution | |
3d86166d | 511 | h2 = (TH2I*)fContainer->At(kTrackletYResolution); |
765bd0ab | 512 | jgraph++; //skip the frame histo |
d85cd79c | 513 | gm = (TGraphErrors*)fContainer->At(jgraph++); |
514 | gs = (TGraphErrors*)fContainer->At(jgraph++); | |
d2381af5 | 515 | for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ |
516 | Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); | |
517 | f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2); | |
518 | h = h2->ProjectionY("py", iphi, iphi); | |
3c3d9ff1 | 519 | Fit(h, &fc); |
d2381af5 | 520 | Int_t jphi = iphi -1; |
765bd0ab | 521 | gm->SetPoint(jphi, phi, 10.*f.GetParameter(1)); |
522 | gm->SetPointError(jphi, 0., 10.*f.GetParError(1)); | |
523 | gs->SetPoint(jphi, phi, 10.*f.GetParameter(2)); | |
524 | gs->SetPointError(jphi, 0., 10.*f.GetParError(2)); | |
d2381af5 | 525 | } |
765bd0ab | 526 | fNRefFigures++; |
d2381af5 | 527 | |
874acced | 528 | // tracklet phi resolution |
3d86166d | 529 | h2 = (TH2I*)fContainer->At(kTrackletAngleResolution); |
765bd0ab | 530 | jgraph++; //skip the frame histo |
d85cd79c | 531 | gm = (TGraphErrors*)fContainer->At(jgraph++); |
532 | gs = (TGraphErrors*)fContainer->At(jgraph++); | |
d2381af5 | 533 | for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ |
534 | Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); | |
535 | f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2); | |
536 | h = h2->ProjectionY("py", iphi, iphi); | |
537 | h->Fit(&f, "QN", "", -.5, .5); | |
538 | Int_t jphi = iphi -1; | |
539 | gm->SetPoint(jphi, phi, f.GetParameter(1)); | |
540 | gm->SetPointError(jphi, 0., f.GetParError(1)); | |
541 | gs->SetPoint(jphi, phi, f.GetParameter(2)); | |
542 | gs->SetPointError(jphi, 0., f.GetParError(2)); | |
543 | } | |
765bd0ab | 544 | fNRefFigures++; |
d2381af5 | 545 | } |
765bd0ab | 546 | |
d85cd79c | 547 | return kTRUE; |
548 | } | |
549 | ||
550 | ||
551 | //________________________________________________________ | |
552 | void AliTRDtrackingResolution::Terminate(Option_t *) | |
553 | { | |
554 | if(fDebugStream){ | |
555 | delete fDebugStream; | |
556 | fDebugStream = 0x0; | |
557 | fDebugLevel = 0; | |
558 | } | |
559 | if(HasPostProcess()) PostProcess(); | |
874acced | 560 | } |
d2381af5 | 561 | |
3c3d9ff1 | 562 | //________________________________________________________ |
563 | void AliTRDtrackingResolution::Fit(TH1 *h, TF1 *f) | |
564 | { | |
565 | // Helper function to avoid duplication of code | |
566 | // Make first guesses on the fit parameters | |
567 | ||
568 | // find the intial parameters of the fit !! (thanks George) | |
569 | Int_t nbinsy = Int_t(.5*h->GetNbinsX()); | |
570 | Double_t sum = 0.; | |
571 | for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.; | |
572 | f->SetParLimits(0, 0., 3.*sum); | |
573 | f->SetParameter(0, .9*sum); | |
574 | ||
575 | f->SetParLimits(1, -.1, .1); | |
576 | f->SetParameter(1, 0.); | |
577 | ||
578 | f->SetParLimits(2, 0., 1.e-1); | |
579 | f->SetParameter(2, 2.e-2); | |
580 | ||
581 | f->SetParLimits(3, 0., sum); | |
582 | f->SetParameter(3, .1*sum); | |
583 | ||
584 | f->SetParLimits(4, -.3, .3); | |
585 | f->SetParameter(4, 0.); | |
586 | ||
587 | f->SetParLimits(5, 0., 1.e2); | |
588 | f->SetParameter(5, 2.e-1); | |
589 | ||
590 | h->Fit(f, "QN", "", -0.5, 0.5); | |
591 | } | |
592 | ||
874acced | 593 | //________________________________________________________ |
3d86166d | 594 | TObjArray* AliTRDtrackingResolution::Histos() |
874acced | 595 | { |
765bd0ab | 596 | if(!fContainer) fContainer = new TObjArray(25); |
3d86166d | 597 | return fContainer; |
77203477 | 598 | } |
599 | ||
aaf47b30 | 600 | |
601 | //________________________________________________________ | |
602 | void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r) | |
603 | { | |
3d86166d | 604 | |
aaf47b30 | 605 | fReconstructor->SetRecoParam(r); |
606 | } |