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