Bug correction
[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();
29fc2c86 79 fMuonTrackRef->Delete();
b8dc484b 80 delete fMuonTrackRef;
81 delete fMUONData;
82}
83
84//_____________________________________________________________________________
85void AliMUONRecoCheck::MakeTrackRef()
86{
87 // Make reconstructible tracks
88
89 AliTrackReference *trackReference;
90 AliMUONTrack *muonTrack;
91 AliMUONTrackParam *trackParam;
92 AliMUONHitForRec *hitForRec;
93 Float_t x, y, z, pX, pY, pZ, pYZ;
94 Int_t track, trackSave;
95 TParticle *particle;
96 Float_t bendingSlope = 0;
97 Float_t nonBendingSlope = 0;
98 Float_t inverseBendingMomentum = 0;
99
100 TTree* treeTR = fRunLoader->TreeTR();
101 if (treeTR == NULL) return;
102
103 TBranch* branch = treeTR->GetBranch("MUON");
104 if (branch == NULL) return;
105
106 TClonesArray* trackRefs = new TClonesArray("AliTrackReference", 10);
107 branch->SetAddress(&trackRefs);
108
109 Int_t nTrackRef = (Int_t)treeTR->GetEntries();
110
111 track = trackSave = -999;
112 Bool_t isNewTrack;
113 Int_t iHitMin, iChamber;
114
115 trackParam = new AliMUONTrackParam();
116 hitForRec = new AliMUONHitForRec();
117 muonTrack = new AliMUONTrack();
118
119 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
120 treeTR->GetEntry(iTrackRef);
121
122 iHitMin = 0;
123 isNewTrack = kTRUE;
dba3c809 124
29fc2c86 125 if (!trackRefs->GetEntries()) continue;
b8dc484b 126
29fc2c86 127 while (isNewTrack) {
128
b8dc484b 129 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
130
131 trackReference = (AliTrackReference*)trackRefs->At(iHit);
132 x = trackReference->X();
133 y = trackReference->Y();
134 z = trackReference->Z();
135 pX = trackReference->Px();
136 pY = trackReference->Py();
137 pZ = trackReference->Pz();
138
139 track = trackReference->GetTrack();
140
141 if (track != trackSave && iHit != 0) {
142 iHitMin = iHit;
143 trackSave = track;
144 break;
145 }
146
147 // track parameters at hit
148 trackParam->SetBendingCoor(y);
149 trackParam->SetNonBendingCoor(x);
150 trackParam->SetZ(z);
151
152 if (TMath::Abs(pZ) > 0) {
153 bendingSlope = pY/pZ;
154 nonBendingSlope = pX/pZ;
155 }
156 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
157 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
158
159 trackParam->SetBendingSlope(bendingSlope);
160 trackParam->SetNonBendingSlope(nonBendingSlope);
161 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
162
163 hitForRec->SetBendingCoor(y);
164 hitForRec->SetNonBendingCoor(x);
165 hitForRec->SetZ(z);
166 hitForRec->SetBendingReso2(0.0);
167 hitForRec->SetNonBendingReso2(0.0);
29fc2c86 168 iChamber = AliMUONConstants::ChamberNumber(z);
b8dc484b 169 hitForRec->SetChamberNumber(iChamber);
170
171 muonTrack->AddTrackParamAtHit(trackParam);
172 muonTrack->AddHitForRecAtHit(hitForRec);
173 muonTrack->SetTrackID(track);
174
175 trackSave = track;
176 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
177 }
178
179 // track parameters at vertex
180 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
181 if (particle) {
182 x = particle->Vx();
183 y = particle->Vy();
184 z = particle->Vz();
185 pX = particle->Px();
186 pY = particle->Py();
187 pZ = particle->Pz();
188
189 trackParam->SetBendingCoor(y);
190 trackParam->SetNonBendingCoor(x);
191 trackParam->SetZ(z);
192 if (TMath::Abs(pZ) > 0) {
193 bendingSlope = pY/pZ;
194 nonBendingSlope = pX/pZ;
195 }
196 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
197 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
198 trackParam->SetBendingSlope(bendingSlope);
199 trackParam->SetNonBendingSlope(nonBendingSlope);
200 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
201
202 muonTrack->SetTrackParamAtVertex(trackParam);
203 }
204
205 AddMuonTrackReference(muonTrack);
206 muonTrack->ResetTrackParamAtHit();
207 muonTrack->ResetHitForRecAtHit();
208
209 } // end while isNewTrack
210
211 }
212
213 CleanMuonTrackRef();
214
215 ReconstructibleTracks();
216
217 delete muonTrack;
218 delete trackParam;
219 delete hitForRec;
29fc2c86 220 trackRefs->Delete();
b8dc484b 221 delete trackRefs;
222
223}
224
225//____________________________________________________________________________
226TClonesArray* AliMUONRecoCheck::GetTrackReco()
227{
228 // Return TClonesArray of reconstructed tracks
229
29fc2c86 230 GetMUONData()->ResetRecTracks();
b8dc484b 231 GetMUONData()->SetTreeAddress("RT");
232 fTrackReco = GetMUONData()->RecTracks();
233 GetMUONData()->GetRecTracks();
234 fRecoTracks = fTrackReco->GetEntriesFast();
235 return fTrackReco;
236}
237//_____________________________________________________________________________
238void AliMUONRecoCheck::PrintEvent() const
239{
240 // debug facility
241
242 AliMUONTrack *track;
243 AliMUONHitForRec *hitForRec;
244 TClonesArray * hitForRecAtHit = 0;
245 Float_t xRec,yRec,zRec;
246
247 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
248
249 printf(" ******************************* \n");
250 printf(" nb of tracks %d \n",nTrackRef);
251
252 for (Int_t index = 0; index < nTrackRef; index++) {
253 track = (AliMUONTrack*)fMuonTrackRef->At(index);
254 hitForRecAtHit = track->GetHitForRecAtHit();
255 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
256 printf(" track number %d \n",index);
257 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
258 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
259 xRec = hitForRec->GetNonBendingCoor();
260 yRec = hitForRec->GetBendingCoor();
261 zRec = hitForRec->GetZ();
262 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
263 }
264 }
265}
266
267//_____________________________________________________________________________
268void AliMUONRecoCheck::ResetTracks() const
269{
270 if (fMuonTrackRef) fMuonTrackRef->Clear();
271}
272//_____________________________________________________________________________
273void AliMUONRecoCheck::CleanMuonTrackRef()
274{
275 // Re-calculate hits parameters because two AliTrackReferences are recorded for
276 // each chamber (one when particle is entering + one when particle is leaving
277 // the sensitive volume)
278
279 Float_t maxGasGap = 1.; // cm
280 AliMUONTrack *track, *trackNew;
281 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
282 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
283 TClonesArray * hitForRecAtHit = 0;
284 TClonesArray * trackParamAtHit = 0;
285 Float_t xRec,yRec,zRec;
286 Float_t xRec1,yRec1,zRec1;
287 Float_t xRec2,yRec2,zRec2;
288 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
289 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
290 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
291 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
292 Int_t iHit1;
293 Int_t iChamber = 0;
294 Int_t nRec = 0;
295 Int_t nTrackHits = 0;
296
297 hitForRec = new AliMUONHitForRec();
298 trackParam = new AliMUONTrackParam();
299 trackNew = new AliMUONTrack();
300
301 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
302
303 for (Int_t index = 0; index < nTrackRef; index++) {
304 track = (AliMUONTrack*)fMuonTrackRef->At(index);
305 hitForRecAtHit = track->GetHitForRecAtHit();
306 trackParamAtHit = track->GetTrackParamAtHit();
307 trackParamAtVertex = track->GetTrackParamAtVertex();
308 nTrackHits = hitForRecAtHit->GetEntriesFast();
309 iHit1 = 0;
310 while (iHit1 < nTrackHits) {
311 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
312 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
313 xRec1 = hitForRec1->GetNonBendingCoor();
314 yRec1 = hitForRec1->GetBendingCoor();
315 zRec1 = hitForRec1->GetZ();
316 xRec = xRec1;
317 yRec = yRec1;
318 zRec = zRec1;
319 bendingSlope1 = trackParam1->GetBendingSlope();
320 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
321 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
322 bendingSlope = bendingSlope1;
323 nonBendingSlope = nonBendingSlope1;
324 bendingMomentum = bendingMomentum1;
325 nRec = 1;
326 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
327 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
328 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
329 xRec2 = hitForRec2->GetNonBendingCoor();
330 yRec2 = hitForRec2->GetBendingCoor();
331 zRec2 = hitForRec2->GetZ();
332 bendingSlope2 = trackParam2->GetBendingSlope();
333 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
334 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
335
336 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
337 nRec++;
338 xRec += xRec2;
339 yRec += yRec2;
340 zRec += zRec2;
341 bendingSlope += bendingSlope2;
342 nonBendingSlope += nonBendingSlope2;
343 bendingMomentum += bendingMomentum2;
344 iHit1 = iHit2;
345 }
346
347 } // end iHit2
348 xRec /= (Float_t)nRec;
349 yRec /= (Float_t)nRec;
350 zRec /= (Float_t)nRec;
351 bendingSlope /= (Float_t)nRec;
352 nonBendingSlope /= (Float_t)nRec;
353 bendingMomentum /= (Float_t)nRec;
354
355 hitForRec->SetNonBendingCoor(xRec);
356 hitForRec->SetBendingCoor(yRec);
357 hitForRec->SetZ(zRec);
29fc2c86 358 iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b 359 hitForRec->SetChamberNumber(iChamber);
360 hitForRec->SetBendingReso2(0.0);
361 hitForRec->SetNonBendingReso2(0.0);
362 trackParam->SetNonBendingCoor(xRec);
363 trackParam->SetBendingCoor(yRec);
364 trackParam->SetZ(zRec);
365 trackParam->SetNonBendingSlope(nonBendingSlope);
366 trackParam->SetBendingSlope(bendingSlope);
367 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
368
369 trackNew->AddHitForRecAtHit(hitForRec);
370 trackNew->AddTrackParamAtHit(trackParam);
371
372 iHit1++;
373 } // end iHit1
374
375 trackNew->SetTrackID(track->GetTrackID());
376 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
377 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
378 trackNew->ResetHitForRecAtHit();
379 trackNew->ResetTrackParamAtHit();
380
381 } // end trackRef
382
383 fMuonTrackRef->Clear();
384 nTrackRef = newMuonTrackRef->GetEntriesFast();
385 for (Int_t index = 0; index < nTrackRef; index++) {
386 track = (AliMUONTrack*)newMuonTrackRef->At(index);
387 AddMuonTrackReference(track);
388 }
389
390
391 delete trackNew;
392 delete hitForRec;
393 delete trackParam;
29fc2c86 394 newMuonTrackRef->Delete();
b8dc484b 395 delete newMuonTrackRef;
396
397}
398
399//_____________________________________________________________________________
400void AliMUONRecoCheck::ReconstructibleTracks()
401{
402 // calculate the number of reconstructible tracks
403
404 AliMUONTrack* track;
405 TClonesArray* hitForRecAtHit = NULL;
406 AliMUONHitForRec* hitForRec;
407 Float_t zRec;
408 Int_t nTrackRef, nTrackHits;
409 Int_t isChamberInTrack[10];
410 Int_t iChamber = 0;
411 Bool_t isTrackOK = kTRUE;
412
413 fReconstructibleTracks = 0;
414
415 nTrackRef = fMuonTrackRef->GetEntriesFast();
416 for (Int_t index = 0; index < nTrackRef; index++) {
417 track = (AliMUONTrack*)fMuonTrackRef->At(index);
418 hitForRecAtHit = track->GetHitForRecAtHit();
419 nTrackHits = hitForRecAtHit->GetEntriesFast();
420 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
421
422 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
423 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
424 zRec = hitForRec->GetZ();
425 iChamber = hitForRec->GetChamberNumber();
426 if (iChamber < 0 || iChamber > 10) continue;
427 isChamberInTrack[iChamber] = 1;
428
429 }
430 // track is reconstructible if the particle is crossing every tracking chambers
431 isTrackOK = kTRUE;
432 for (Int_t ch = 0; ch < 10; ch++) {
433 if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
29fc2c86 434 } if (isTrackOK) fReconstructibleTracks++;
b8dc484b 435 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
436 }
437 fMuonTrackRef->Compress();
438}
439
440