1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // -----------------------
19 // Class AliMUONRecoCheck
20 // -----------------------
21 // Utility class to check the muon reconstruction. Reconstructed tracks are compared
22 // to reference tracks. The reference tracks are built from AliTrackReference for the
23 // hit in chamber (0..9) and from kinematics for the vertex parameters.
26 #include "AliMUONRecoCheck.h"
27 #include "AliMUONTrack.h"
28 #include "AliMUONData.h"
29 #include "AliMUONConstants.h"
31 #include "AliRun.h" // for gAlice
32 #include "AliLoader.h"
33 #include "AliRunLoader.h"
34 #include "AliHeader.h"
36 #include "AliTrackReference.h"
39 #include <TParticle.h>
41 ClassImp(AliMUONRecoCheck)
43 //_____________________________________________________________________________
44 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
47 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
49 // open the run loader
50 fRunLoader = AliRunLoader::Open(chLoader);
52 AliError(Form("no run loader found in file %s","galice.root" ));
55 // initialize loader's
56 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
58 // initialize container
59 fMUONData = new AliMUONData(loader,"MUON","MUON");
61 // Loading AliRun master
62 if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
64 fRunLoader->LoadKinematics("READ");
65 fRunLoader->LoadTrackRefs("READ");
66 loader->LoadTracks("READ");
68 fReconstructibleTracks = 0;
72 //____________________________________________________________________
73 AliMUONRecoCheck::AliMUONRecoCheck(const AliMUONRecoCheck& rhs)
76 /// Protected copy constructor
78 AliFatal("Not implemented.");
81 //_____________________________________________________________________________
82 AliMUONRecoCheck::~AliMUONRecoCheck()
86 fRunLoader->UnloadKinematics();
87 fRunLoader->UnloadTrackRefs();
88 fRunLoader->UnloadTracks();
89 fMuonTrackRef->Delete();
94 //________________________________________________________________________
95 AliMUONRecoCheck& AliMUONRecoCheck::operator = (const AliMUONRecoCheck& rhs)
97 /// Protected assignement operator
99 if (this == &rhs) return *this;
101 AliFatal("Not implemented.");
106 //_____________________________________________________________________________
107 void AliMUONRecoCheck::MakeTrackRef()
109 /// Make reconstructible tracks
111 AliTrackReference *trackReference;
112 AliMUONTrack *muonTrack;
113 AliMUONTrackParam *trackParam;
114 AliMUONHitForRec *hitForRec;
115 Float_t x, y, z, pX, pY, pZ, pYZ;
116 Int_t track, trackSave;
118 Float_t bendingSlope = 0;
119 Float_t nonBendingSlope = 0;
120 Float_t inverseBendingMomentum = 0;
122 TTree* treeTR = fRunLoader->TreeTR();
123 if (treeTR == NULL) return;
125 TBranch* branch = treeTR->GetBranch("MUON");
126 if (branch == NULL) return;
128 TClonesArray* trackRefs = 0;
129 branch->SetAddress(&trackRefs);
130 branch->SetAutoDelete(kTRUE);
132 Int_t nTrackRef = (Int_t)branch->GetEntries();
134 track = trackSave = -999;
136 Int_t iHitMin, iChamber;
138 trackParam = new AliMUONTrackParam();
139 hitForRec = new AliMUONHitForRec();
140 muonTrack = new AliMUONTrack();
142 Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
143 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
145 branch->GetEntry(iTrackRef);
150 if (!trackRefs->GetEntries()) continue;
153 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
155 trackReference = (AliTrackReference*)trackRefs->At(iHit);
156 x = trackReference->X();
157 y = trackReference->Y();
158 z = trackReference->Z();
159 pX = trackReference->Px();
160 pY = trackReference->Py();
161 pZ = trackReference->Pz();
163 track = trackReference->GetTrack();
167 << "Track ID " << track
168 << " larger than max number of particles " << max << endl;
172 if (track != trackSave && iHit != 0) {
178 // track parameters at hit
179 trackParam->SetBendingCoor(y);
180 trackParam->SetNonBendingCoor(x);
183 if (TMath::Abs(pZ) > 0) {
184 bendingSlope = pY/pZ;
185 nonBendingSlope = pX/pZ;
187 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
188 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
190 trackParam->SetBendingSlope(bendingSlope);
191 trackParam->SetNonBendingSlope(nonBendingSlope);
192 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
194 hitForRec->SetBendingCoor(y);
195 hitForRec->SetNonBendingCoor(x);
197 hitForRec->SetBendingReso2(0.0);
198 hitForRec->SetNonBendingReso2(0.0);
199 iChamber = AliMUONConstants::ChamberNumber(z);
200 hitForRec->SetChamberNumber(iChamber);
202 muonTrack->AddTrackParamAtHit(trackParam);
203 muonTrack->AddHitForRecAtHit(hitForRec);
204 muonTrack->SetTrackID(track);
207 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
210 // track parameters at vertex
211 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
222 trackParam->SetBendingCoor(y);
223 trackParam->SetNonBendingCoor(x);
225 if (TMath::Abs(pZ) > 0) {
226 bendingSlope = pY/pZ;
227 nonBendingSlope = pX/pZ;
229 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
230 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
231 trackParam->SetBendingSlope(bendingSlope);
232 trackParam->SetNonBendingSlope(nonBendingSlope);
233 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
235 muonTrack->SetTrackParamAtVertex(trackParam);
238 AddMuonTrackReference(muonTrack);
239 muonTrack->ResetTrackParamAtHit();
240 muonTrack->ResetHitForRecAtHit();
242 } // end while isNewTrack
247 ReconstructibleTracks();
255 //____________________________________________________________________________
256 TClonesArray* AliMUONRecoCheck::GetTrackReco()
258 /// Return TClonesArray of reconstructed tracks
260 GetMUONData()->ResetRecTracks();
261 GetMUONData()->SetTreeAddress("RT");
262 fTrackReco = GetMUONData()->RecTracks();
263 GetMUONData()->GetRecTracks();
264 fRecoTracks = fTrackReco->GetEntriesFast();
267 //_____________________________________________________________________________
268 void AliMUONRecoCheck::PrintEvent() const
273 AliMUONHitForRec *hitForRec;
274 TClonesArray * hitForRecAtHit = 0;
275 Float_t xRec,yRec,zRec;
277 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
279 printf(" ******************************* \n");
280 printf(" nb of tracks %d \n",nTrackRef);
282 for (Int_t index = 0; index < nTrackRef; index++) {
283 track = (AliMUONTrack*)fMuonTrackRef->At(index);
284 hitForRecAtHit = track->GetHitForRecAtHit();
285 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
286 printf(" track number %d \n",index);
287 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
288 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
289 xRec = hitForRec->GetNonBendingCoor();
290 yRec = hitForRec->GetBendingCoor();
291 zRec = hitForRec->GetZ();
292 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
297 //_____________________________________________________________________________
298 void AliMUONRecoCheck::ResetTracks() const
302 if (fMuonTrackRef) fMuonTrackRef->Clear();
304 //_____________________________________________________________________________
305 void AliMUONRecoCheck::CleanMuonTrackRef()
307 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
308 /// each chamber (one when particle is entering + one when particle is leaving
309 /// the sensitive volume)
311 Float_t maxGasGap = 1.; // cm
312 AliMUONTrack *track, *trackNew;
313 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
314 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
315 TClonesArray * hitForRecAtHit = 0;
316 TClonesArray * trackParamAtHit = 0;
317 Float_t xRec,yRec,zRec;
318 Float_t xRec1,yRec1,zRec1;
319 Float_t xRec2,yRec2,zRec2;
320 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
321 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
322 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
323 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
327 Int_t nTrackHits = 0;
329 hitForRec = new AliMUONHitForRec();
330 trackParam = new AliMUONTrackParam();
331 trackNew = new AliMUONTrack();
333 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
334 for (Int_t index = 0; index < nTrackRef; index++) {
335 track = (AliMUONTrack*)fMuonTrackRef->At(index);
336 hitForRecAtHit = track->GetHitForRecAtHit();
337 trackParamAtHit = track->GetTrackParamAtHit();
338 trackParamAtVertex = track->GetTrackParamAtVertex();
339 nTrackHits = hitForRecAtHit->GetEntriesFast();
341 while (iHit1 < nTrackHits) {
342 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
343 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
344 xRec1 = hitForRec1->GetNonBendingCoor();
345 yRec1 = hitForRec1->GetBendingCoor();
346 zRec1 = hitForRec1->GetZ();
350 bendingSlope1 = trackParam1->GetBendingSlope();
351 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
352 bendingMomentum1 = 0;
353 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
354 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
355 bendingSlope = bendingSlope1;
356 nonBendingSlope = nonBendingSlope1;
357 bendingMomentum = bendingMomentum1;
359 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
360 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
361 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
362 xRec2 = hitForRec2->GetNonBendingCoor();
363 yRec2 = hitForRec2->GetBendingCoor();
364 zRec2 = hitForRec2->GetZ();
365 bendingSlope2 = trackParam2->GetBendingSlope();
366 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
367 bendingMomentum2 = 0;
368 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
369 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
371 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
377 bendingSlope += bendingSlope2;
378 nonBendingSlope += nonBendingSlope2;
379 bendingMomentum += bendingMomentum2;
384 xRec /= (Float_t)nRec;
385 yRec /= (Float_t)nRec;
386 zRec /= (Float_t)nRec;
387 bendingSlope /= (Float_t)nRec;
388 nonBendingSlope /= (Float_t)nRec;
389 bendingMomentum /= (Float_t)nRec;
391 hitForRec->SetNonBendingCoor(xRec);
392 hitForRec->SetBendingCoor(yRec);
393 hitForRec->SetZ(zRec);
394 iChamber = AliMUONConstants::ChamberNumber(zRec);
395 hitForRec->SetChamberNumber(iChamber);
396 hitForRec->SetBendingReso2(0.0);
397 hitForRec->SetNonBendingReso2(0.0);
398 trackParam->SetNonBendingCoor(xRec);
399 trackParam->SetBendingCoor(yRec);
400 trackParam->SetZ(zRec);
401 trackParam->SetNonBendingSlope(nonBendingSlope);
402 trackParam->SetBendingSlope(bendingSlope);
403 if (TMath::Abs(bendingMomentum) > 0)
404 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
406 trackNew->AddHitForRecAtHit(hitForRec);
407 trackNew->AddTrackParamAtHit(trackParam);
412 trackNew->SetTrackID(track->GetTrackID());
413 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
414 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
415 trackNew->ResetHitForRecAtHit();
416 trackNew->ResetTrackParamAtHit();
420 fMuonTrackRef->Clear();
421 nTrackRef = newMuonTrackRef->GetEntriesFast();
422 for (Int_t index = 0; index < nTrackRef; index++) {
423 track = (AliMUONTrack*)newMuonTrackRef->At(index);
424 AddMuonTrackReference(track);
431 newMuonTrackRef->Delete();
432 delete newMuonTrackRef;
436 //_____________________________________________________________________________
437 void AliMUONRecoCheck::ReconstructibleTracks()
439 /// Calculate the number of reconstructible tracks
442 TClonesArray* hitForRecAtHit = NULL;
443 AliMUONHitForRec* hitForRec;
445 Int_t nTrackRef, nTrackHits;
446 Int_t isChamberInTrack[10];
448 Bool_t isTrackOK = kTRUE;
450 fReconstructibleTracks = 0;
452 nTrackRef = fMuonTrackRef->GetEntriesFast();
453 for (Int_t index = 0; index < nTrackRef; index++) {
454 track = (AliMUONTrack*)fMuonTrackRef->At(index);
455 hitForRecAtHit = track->GetHitForRecAtHit();
456 nTrackHits = hitForRecAtHit->GetEntriesFast();
457 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
459 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
460 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
461 zRec = hitForRec->GetZ();
462 iChamber = hitForRec->GetChamberNumber();
463 if (iChamber < 0 || iChamber > 10) continue;
464 isChamberInTrack[iChamber] = 1;
466 // track is reconstructible if the particle is depositing a hit
467 // in the following chamber combinations:
470 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
471 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
472 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
473 Int_t nHitsInLastStations=0;
474 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
475 if (isChamberInTrack[ch]) nHitsInLastStations++;
476 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
478 if (isTrackOK) fReconstructibleTracks++;
479 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
481 fMuonTrackRef->Compress();