New class added to check the reconstruction of muon tracks using reference tracks
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
CommitLineData
b8dc484b 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-commercial 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//////////////////////////////////////////////////////////////////////////////////////
17//
18// Utility class to check the muon reconstruction. Reconstructed tracks are compared
19// to reference tracks. The reference tracks are built from AliTrackReference for the
20// hit in chamber (0..9) and from kinematics for the vertex parameters.
21//
22//////////////////////////////////////////////////////////////////////////////////////
23
24
25#include <TParticle.h>
26
27#include "AliRun.h" // for gAlice
28#include "AliLog.h"
29#include "AliLoader.h"
30#include "AliRunLoader.h"
31#include "AliTrackReference.h"
32#include "AliHeader.h"
33#include "AliMC.h"
34#include "AliStack.h"
35#include "AliMUON.h"
36#include "AliMUONRecoCheck.h"
37#include "AliMUONTrack.h"
38#include "AliMUONData.h"
39#include "AliMUONChamber.h"
40#include "AliMUONConstants.h"
41
42ClassImp(AliMUONRecoCheck)
43
44//_____________________________________________________________________________
45AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
46{
47// Constructor
48 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
49
50 // open the run loader
51 fRunLoader = AliRunLoader::Open(chLoader);
52 if (!fRunLoader) {
53 AliError(Form("no run loader found in file %s","galice.root" ));
54 return;
55 }
56 // initialize loader's
57 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
58
59 // initialize container
60 fMUONData = new AliMUONData(loader,"MUON","MUON");
61
62 // Loading AliRun master
63 if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
64
65 fRunLoader->LoadKinematics("READ");
66 fRunLoader->LoadTrackRefs("READ");
67 loader->LoadTracks("READ");
68
69 fReconstructibleTracks = 0;
70 fRecoTracks = 0;
71}
72
73//_____________________________________________________________________________
74AliMUONRecoCheck::~AliMUONRecoCheck()
75{
76 fRunLoader->UnloadKinematics();
77 fRunLoader->UnloadTrackRefs();
78 fRunLoader->UnloadTracks();
79 delete fMuonTrackRef;
80 delete fMUONData;
81}
82
83//_____________________________________________________________________________
84void AliMUONRecoCheck::MakeTrackRef()
85{
86 // Make reconstructible tracks
87
88 AliTrackReference *trackReference;
89 AliMUONTrack *muonTrack;
90 AliMUONTrackParam *trackParam;
91 AliMUONHitForRec *hitForRec;
92 Float_t x, y, z, pX, pY, pZ, pYZ;
93 Int_t track, trackSave;
94 TParticle *particle;
95 Float_t bendingSlope = 0;
96 Float_t nonBendingSlope = 0;
97 Float_t inverseBendingMomentum = 0;
98
99 TTree* treeTR = fRunLoader->TreeTR();
100 if (treeTR == NULL) return;
101
102 TBranch* branch = treeTR->GetBranch("MUON");
103 if (branch == NULL) return;
104
105 TClonesArray* trackRefs = new TClonesArray("AliTrackReference", 10);
106 branch->SetAddress(&trackRefs);
107
108 Int_t nTrackRef = (Int_t)treeTR->GetEntries();
109
110 track = trackSave = -999;
111 Bool_t isNewTrack;
112 Int_t iHitMin, iChamber;
113
114 trackParam = new AliMUONTrackParam();
115 hitForRec = new AliMUONHitForRec();
116 muonTrack = new AliMUONTrack();
117
118 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
119 treeTR->GetEntry(iTrackRef);
120
121 iHitMin = 0;
122 isNewTrack = kTRUE;
123
124 while (isNewTrack) {
125
126 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
127
128 trackReference = (AliTrackReference*)trackRefs->At(iHit);
129 x = trackReference->X();
130 y = trackReference->Y();
131 z = trackReference->Z();
132 pX = trackReference->Px();
133 pY = trackReference->Py();
134 pZ = trackReference->Pz();
135
136 track = trackReference->GetTrack();
137
138 if (track != trackSave && iHit != 0) {
139 iHitMin = iHit;
140 trackSave = track;
141 break;
142 }
143
144 // track parameters at hit
145 trackParam->SetBendingCoor(y);
146 trackParam->SetNonBendingCoor(x);
147 trackParam->SetZ(z);
148
149 if (TMath::Abs(pZ) > 0) {
150 bendingSlope = pY/pZ;
151 nonBendingSlope = pX/pZ;
152 }
153 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
154 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
155
156 trackParam->SetBendingSlope(bendingSlope);
157 trackParam->SetNonBendingSlope(nonBendingSlope);
158 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
159
160 hitForRec->SetBendingCoor(y);
161 hitForRec->SetNonBendingCoor(x);
162 hitForRec->SetZ(z);
163 hitForRec->SetBendingReso2(0.0);
164 hitForRec->SetNonBendingReso2(0.0);
165 iChamber = ChamberNumber(z);
166 hitForRec->SetChamberNumber(iChamber);
167
168 muonTrack->AddTrackParamAtHit(trackParam);
169 muonTrack->AddHitForRecAtHit(hitForRec);
170 muonTrack->SetTrackID(track);
171
172 trackSave = track;
173 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
174 }
175
176 // track parameters at vertex
177 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
178 if (particle) {
179 x = particle->Vx();
180 y = particle->Vy();
181 z = particle->Vz();
182 pX = particle->Px();
183 pY = particle->Py();
184 pZ = particle->Pz();
185
186 trackParam->SetBendingCoor(y);
187 trackParam->SetNonBendingCoor(x);
188 trackParam->SetZ(z);
189 if (TMath::Abs(pZ) > 0) {
190 bendingSlope = pY/pZ;
191 nonBendingSlope = pX/pZ;
192 }
193 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
194 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
195 trackParam->SetBendingSlope(bendingSlope);
196 trackParam->SetNonBendingSlope(nonBendingSlope);
197 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
198
199 muonTrack->SetTrackParamAtVertex(trackParam);
200 }
201
202 AddMuonTrackReference(muonTrack);
203 muonTrack->ResetTrackParamAtHit();
204 muonTrack->ResetHitForRecAtHit();
205
206 } // end while isNewTrack
207
208 }
209
210 CleanMuonTrackRef();
211
212 ReconstructibleTracks();
213
214 delete muonTrack;
215 delete trackParam;
216 delete hitForRec;
217 delete trackRefs;
218
219}
220
221//____________________________________________________________________________
222TClonesArray* AliMUONRecoCheck::GetTrackReco()
223{
224 // Return TClonesArray of reconstructed tracks
225
226 GetMUONData()->SetTreeAddress("RT");
227 fTrackReco = GetMUONData()->RecTracks();
228 GetMUONData()->GetRecTracks();
229 fRecoTracks = fTrackReco->GetEntriesFast();
230 return fTrackReco;
231}
232//_____________________________________________________________________________
233void AliMUONRecoCheck::PrintEvent() const
234{
235 // debug facility
236
237 AliMUONTrack *track;
238 AliMUONHitForRec *hitForRec;
239 TClonesArray * hitForRecAtHit = 0;
240 Float_t xRec,yRec,zRec;
241
242 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
243
244 printf(" ******************************* \n");
245 printf(" nb of tracks %d \n",nTrackRef);
246
247 for (Int_t index = 0; index < nTrackRef; index++) {
248 track = (AliMUONTrack*)fMuonTrackRef->At(index);
249 hitForRecAtHit = track->GetHitForRecAtHit();
250 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
251 printf(" track number %d \n",index);
252 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
253 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
254 xRec = hitForRec->GetNonBendingCoor();
255 yRec = hitForRec->GetBendingCoor();
256 zRec = hitForRec->GetZ();
257 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
258 }
259 }
260}
261
262//_____________________________________________________________________________
263void AliMUONRecoCheck::ResetTracks() const
264{
265 if (fMuonTrackRef) fMuonTrackRef->Clear();
266}
267//_____________________________________________________________________________
268void AliMUONRecoCheck::CleanMuonTrackRef()
269{
270 // Re-calculate hits parameters because two AliTrackReferences are recorded for
271 // each chamber (one when particle is entering + one when particle is leaving
272 // the sensitive volume)
273
274 Float_t maxGasGap = 1.; // cm
275 AliMUONTrack *track, *trackNew;
276 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
277 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
278 TClonesArray * hitForRecAtHit = 0;
279 TClonesArray * trackParamAtHit = 0;
280 Float_t xRec,yRec,zRec;
281 Float_t xRec1,yRec1,zRec1;
282 Float_t xRec2,yRec2,zRec2;
283 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
284 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
285 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
286 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
287 Int_t iHit1;
288 Int_t iChamber = 0;
289 Int_t nRec = 0;
290 Int_t nTrackHits = 0;
291
292 hitForRec = new AliMUONHitForRec();
293 trackParam = new AliMUONTrackParam();
294 trackNew = new AliMUONTrack();
295
296 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
297
298 for (Int_t index = 0; index < nTrackRef; index++) {
299 track = (AliMUONTrack*)fMuonTrackRef->At(index);
300 hitForRecAtHit = track->GetHitForRecAtHit();
301 trackParamAtHit = track->GetTrackParamAtHit();
302 trackParamAtVertex = track->GetTrackParamAtVertex();
303 nTrackHits = hitForRecAtHit->GetEntriesFast();
304 iHit1 = 0;
305 while (iHit1 < nTrackHits) {
306 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
307 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
308 xRec1 = hitForRec1->GetNonBendingCoor();
309 yRec1 = hitForRec1->GetBendingCoor();
310 zRec1 = hitForRec1->GetZ();
311 xRec = xRec1;
312 yRec = yRec1;
313 zRec = zRec1;
314 bendingSlope1 = trackParam1->GetBendingSlope();
315 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
316 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
317 bendingSlope = bendingSlope1;
318 nonBendingSlope = nonBendingSlope1;
319 bendingMomentum = bendingMomentum1;
320 nRec = 1;
321 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
322 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
323 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
324 xRec2 = hitForRec2->GetNonBendingCoor();
325 yRec2 = hitForRec2->GetBendingCoor();
326 zRec2 = hitForRec2->GetZ();
327 bendingSlope2 = trackParam2->GetBendingSlope();
328 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
329 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
330
331 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
332 nRec++;
333 xRec += xRec2;
334 yRec += yRec2;
335 zRec += zRec2;
336 bendingSlope += bendingSlope2;
337 nonBendingSlope += nonBendingSlope2;
338 bendingMomentum += bendingMomentum2;
339 iHit1 = iHit2;
340 }
341
342 } // end iHit2
343 xRec /= (Float_t)nRec;
344 yRec /= (Float_t)nRec;
345 zRec /= (Float_t)nRec;
346 bendingSlope /= (Float_t)nRec;
347 nonBendingSlope /= (Float_t)nRec;
348 bendingMomentum /= (Float_t)nRec;
349
350 hitForRec->SetNonBendingCoor(xRec);
351 hitForRec->SetBendingCoor(yRec);
352 hitForRec->SetZ(zRec);
353 iChamber = ChamberNumber(zRec);
354 hitForRec->SetChamberNumber(iChamber);
355 hitForRec->SetBendingReso2(0.0);
356 hitForRec->SetNonBendingReso2(0.0);
357 trackParam->SetNonBendingCoor(xRec);
358 trackParam->SetBendingCoor(yRec);
359 trackParam->SetZ(zRec);
360 trackParam->SetNonBendingSlope(nonBendingSlope);
361 trackParam->SetBendingSlope(bendingSlope);
362 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
363
364 trackNew->AddHitForRecAtHit(hitForRec);
365 trackNew->AddTrackParamAtHit(trackParam);
366
367 iHit1++;
368 } // end iHit1
369
370 trackNew->SetTrackID(track->GetTrackID());
371 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
372 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
373 trackNew->ResetHitForRecAtHit();
374 trackNew->ResetTrackParamAtHit();
375
376 } // end trackRef
377
378 fMuonTrackRef->Clear();
379 nTrackRef = newMuonTrackRef->GetEntriesFast();
380 for (Int_t index = 0; index < nTrackRef; index++) {
381 track = (AliMUONTrack*)newMuonTrackRef->At(index);
382 AddMuonTrackReference(track);
383 }
384
385
386 delete trackNew;
387 delete hitForRec;
388 delete trackParam;
389 delete newMuonTrackRef;
390
391}
392
393//_____________________________________________________________________________
394void AliMUONRecoCheck::ReconstructibleTracks()
395{
396 // calculate the number of reconstructible tracks
397
398 AliMUONTrack* track;
399 TClonesArray* hitForRecAtHit = NULL;
400 AliMUONHitForRec* hitForRec;
401 Float_t zRec;
402 Int_t nTrackRef, nTrackHits;
403 Int_t isChamberInTrack[10];
404 Int_t iChamber = 0;
405 Bool_t isTrackOK = kTRUE;
406
407 fReconstructibleTracks = 0;
408
409 nTrackRef = fMuonTrackRef->GetEntriesFast();
410 for (Int_t index = 0; index < nTrackRef; index++) {
411 track = (AliMUONTrack*)fMuonTrackRef->At(index);
412 hitForRecAtHit = track->GetHitForRecAtHit();
413 nTrackHits = hitForRecAtHit->GetEntriesFast();
414 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
415
416 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
417 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
418 zRec = hitForRec->GetZ();
419 iChamber = hitForRec->GetChamberNumber();
420 if (iChamber < 0 || iChamber > 10) continue;
421 isChamberInTrack[iChamber] = 1;
422
423 }
424 // track is reconstructible if the particle is crossing every tracking chambers
425 isTrackOK = kTRUE;
426 for (Int_t ch = 0; ch < 10; ch++) {
427 if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
428 }
429 if (isTrackOK) fReconstructibleTracks++;
430 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
431 }
432 fMuonTrackRef->Compress();
433}
434
435
436//_____________________________________________________________________________
437Int_t AliMUONRecoCheck::ChamberNumber(Float_t z) const
438{
439 // return chamber number according z position of hit. Should be taken from geometry ?
440
441 Float_t dMaxChamber = AliMUONConstants::DzSlat() + AliMUONConstants::DzCh() + 0.25; // cm st 3 &4 & 5
442 if ( z > (AliMUONConstants::DefaultChamberZ(4)+50.)) dMaxChamber = 7.; // cm stations 1 & 2
443 Int_t iChamber;
444
445 for (iChamber = 0; iChamber < 10; iChamber++) {
446
447 if (TMath::Abs(z-AliMUONConstants::DefaultChamberZ(iChamber)) < dMaxChamber) {
448 return iChamber;
449 }
450 }
451 return -1;
452}