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 /// \class AliMUONRecoCheck
19 /// \brief Utility class to check reconstruction
20 /// Reconstructed tracks are compared to reference tracks.
21 /// The reference tracks are built from AliTrackReference for the
22 /// hit in chamber (0..9) and from kinematics for the vertex parameters.
25 #include "AliMUONRecoCheck.h"
26 #include "AliMUONTrack.h"
27 #include "AliMUONData.h"
28 #include "AliMUONConstants.h"
30 #include "AliLoader.h"
31 #include "AliRunLoader.h"
32 #include "AliHeader.h"
34 #include "AliTrackReference.h"
37 #include <TParticle.h>
40 ClassImp(AliMUONRecoCheck)
43 //_____________________________________________________________________________
44 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
50 fReconstructibleTracks(0),
52 fIsLoadConstructor(kTRUE)
54 /// Constructor using "galice.root",
55 /// takes care of loading/unloading Kinematics, TrackRefs and Tracks
57 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
60 fRunLoader = AliRunLoader::Open(chLoader);
62 AliError(Form("no run loader found " ));
67 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
70 fMUONData = new AliMUONData(loader,"MUON","MUON");
72 AliError(Form("no MUONData found " ));
76 fRunLoader->LoadKinematics("READ");
77 fRunLoader->LoadTrackRefs("READ");
78 loader->LoadTracks("READ");
82 //_____________________________________________________________________________
83 AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONData *muondata)
89 fReconstructibleTracks(0),
91 fIsLoadConstructor(kFALSE)
93 /// Constructor from AliRunLoader and AliMUONData,
94 /// does not load/unload Kinematics, TrackRefs and Tracks internally
95 /// it should be done in the execution program (e.g. see MUONRecoCheck.C)
97 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
100 fRunLoader = runloader;
102 AliError(Form("no run loader found " ));
107 fMUONData = muondata;
109 AliError(Form("no MUONData found " ));
116 //_____________________________________________________________________________
117 AliMUONRecoCheck::~AliMUONRecoCheck()
121 delete fMuonTrackRef;
123 if(fIsLoadConstructor){
124 fRunLoader->UnloadKinematics();
125 fRunLoader->UnloadTrackRefs();
126 fRunLoader->UnloadTracks();
132 //_____________________________________________________________________________
133 void AliMUONRecoCheck::MakeTrackRef()
135 /// Make reconstructible tracks
137 AliTrackReference *trackReference;
138 AliMUONTrack *muonTrack;
139 AliMUONTrackParam *trackParam;
140 AliMUONHitForRec *hitForRec;
141 Float_t x, y, z, pX, pY, pZ, pYZ;
142 Int_t track, trackSave;
144 Float_t bendingSlope = 0;
145 Float_t nonBendingSlope = 0;
146 Float_t inverseBendingMomentum = 0;
148 TTree* treeTR = fRunLoader->TreeTR();
149 if (treeTR == NULL) return;
151 TBranch* branch = treeTR->GetBranch("MUON");
152 if (branch == NULL) return;
154 TClonesArray* trackRefs = 0;
155 branch->SetAddress(&trackRefs);
156 branch->SetAutoDelete(kTRUE);
158 Int_t nTrackRef = (Int_t)branch->GetEntries();
160 track = trackSave = -999;
162 Int_t iHitMin, iChamber, detElemId;
164 trackParam = new AliMUONTrackParam();
165 hitForRec = new AliMUONHitForRec();
167 Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
168 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
170 branch->GetEntry(iTrackRef);
175 if (!trackRefs->GetEntries()) continue;
179 muonTrack = new AliMUONTrack();
181 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
183 trackReference = (AliTrackReference*)trackRefs->At(iHit);
184 x = trackReference->X();
185 y = trackReference->Y();
186 z = trackReference->Z();
187 pX = trackReference->Px();
188 pY = trackReference->Py();
189 pZ = trackReference->Pz();
191 track = trackReference->GetTrack();
195 << "Track ID " << track
196 << " larger than max number of particles " << max << endl;
200 if (track != trackSave && iHit != 0) {
206 // track parameters at hit
207 trackParam->SetBendingCoor(y);
208 trackParam->SetNonBendingCoor(x);
211 if (TMath::Abs(pZ) > 0) {
212 bendingSlope = pY/pZ;
213 nonBendingSlope = pX/pZ;
215 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
216 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
218 trackParam->SetBendingSlope(bendingSlope);
219 trackParam->SetNonBendingSlope(nonBendingSlope);
220 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
222 hitForRec->SetBendingCoor(y);
223 hitForRec->SetNonBendingCoor(x);
225 hitForRec->SetBendingReso2(0.0);
226 hitForRec->SetNonBendingReso2(0.0);
227 detElemId = hitForRec->GetDetElemId();
228 if (detElemId) iChamber = detElemId / 100 - 1;
229 else iChamber = AliMUONConstants::ChamberNumber(z);
230 hitForRec->SetChamberNumber(iChamber);
232 muonTrack->AddTrackParamAtHit(trackParam,0);
233 muonTrack->AddHitForRecAtHit(hitForRec);
234 muonTrack->SetTrackID(track);
237 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
240 // track parameters at vertex
241 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
252 trackParam->SetBendingCoor(y);
253 trackParam->SetNonBendingCoor(x);
255 if (TMath::Abs(pZ) > 0) {
256 bendingSlope = pY/pZ;
257 nonBendingSlope = pX/pZ;
259 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
260 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
261 trackParam->SetBendingSlope(bendingSlope);
262 trackParam->SetNonBendingSlope(nonBendingSlope);
263 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
265 muonTrack->SetTrackParamAtVertex(trackParam);
268 AddMuonTrackReference(muonTrack);
272 } // end while isNewTrack
277 ReconstructibleTracks();
284 //____________________________________________________________________________
285 TClonesArray* AliMUONRecoCheck::GetTrackReco()
287 /// Return TClonesArray of reconstructed tracks
289 fMUONData->ResetRecTracks();
290 fMUONData->SetTreeAddress("RT");
291 fTrackReco = fMUONData->RecTracks();
292 fMUONData->GetRecTracks();
293 fRecoTracks = fTrackReco->GetEntriesFast();
296 //_____________________________________________________________________________
297 void AliMUONRecoCheck::PrintEvent() const
302 AliMUONHitForRec *hitForRec;
303 TClonesArray * hitForRecAtHit = 0;
304 Float_t xRec,yRec,zRec;
306 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
308 printf(" ******************************* \n");
309 printf(" nb of tracks %d \n",nTrackRef);
311 for (Int_t index = 0; index < nTrackRef; index++) {
312 track = (AliMUONTrack*)fMuonTrackRef->At(index);
313 hitForRecAtHit = track->GetHitForRecAtHit();
314 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
315 printf(" track number %d \n",index);
316 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
317 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
318 xRec = hitForRec->GetNonBendingCoor();
319 yRec = hitForRec->GetBendingCoor();
320 zRec = hitForRec->GetZ();
321 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
326 //_____________________________________________________________________________
327 void AliMUONRecoCheck::ResetTracks() const
331 if (fMuonTrackRef) fMuonTrackRef->Clear();
333 //_____________________________________________________________________________
334 void AliMUONRecoCheck::CleanMuonTrackRef()
336 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
337 /// each chamber (one when particle is entering + one when particle is leaving
338 /// the sensitive volume)
340 Float_t maxGasGap = 1.; // cm
341 AliMUONTrack *track, *trackNew;
342 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
343 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
344 TClonesArray * hitForRecAtHit = 0;
345 TClonesArray * trackParamAtHit = 0;
346 Float_t xRec,yRec,zRec;
347 Float_t xRec1,yRec1,zRec1;
348 Float_t xRec2,yRec2,zRec2;
349 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
350 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
351 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
352 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
354 Int_t iChamber = 0, detElemId = 0;
356 Int_t nTrackHits = 0;
358 hitForRec = new AliMUONHitForRec();
359 trackParam = new AliMUONTrackParam();
361 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
362 for (Int_t index = 0; index < nTrackRef; index++) {
363 track = (AliMUONTrack*)fMuonTrackRef->At(index);
364 hitForRecAtHit = track->GetHitForRecAtHit();
365 trackParamAtHit = track->GetTrackParamAtHit();
366 trackParamAtVertex = track->GetTrackParamAtVertex();
367 nTrackHits = hitForRecAtHit->GetEntriesFast();
368 trackNew = new AliMUONTrack();
370 while (iHit1 < nTrackHits) {
371 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
372 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
373 xRec1 = hitForRec1->GetNonBendingCoor();
374 yRec1 = hitForRec1->GetBendingCoor();
375 zRec1 = hitForRec1->GetZ();
379 bendingSlope1 = trackParam1->GetBendingSlope();
380 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
381 bendingMomentum1 = 0;
382 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
383 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
384 bendingSlope = bendingSlope1;
385 nonBendingSlope = nonBendingSlope1;
386 bendingMomentum = bendingMomentum1;
388 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
389 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
390 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
391 xRec2 = hitForRec2->GetNonBendingCoor();
392 yRec2 = hitForRec2->GetBendingCoor();
393 zRec2 = hitForRec2->GetZ();
394 bendingSlope2 = trackParam2->GetBendingSlope();
395 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
396 bendingMomentum2 = 0;
397 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
398 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
400 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
406 bendingSlope += bendingSlope2;
407 nonBendingSlope += nonBendingSlope2;
408 bendingMomentum += bendingMomentum2;
413 xRec /= (Float_t)nRec;
414 yRec /= (Float_t)nRec;
415 zRec /= (Float_t)nRec;
416 bendingSlope /= (Float_t)nRec;
417 nonBendingSlope /= (Float_t)nRec;
418 bendingMomentum /= (Float_t)nRec;
420 hitForRec->SetNonBendingCoor(xRec);
421 hitForRec->SetBendingCoor(yRec);
422 hitForRec->SetZ(zRec);
423 detElemId = hitForRec->GetDetElemId();
424 if (detElemId) iChamber = detElemId / 100 - 1;
425 else iChamber = AliMUONConstants::ChamberNumber(zRec);
426 hitForRec->SetChamberNumber(iChamber);
427 hitForRec->SetBendingReso2(0.0);
428 hitForRec->SetNonBendingReso2(0.0);
429 trackParam->SetNonBendingCoor(xRec);
430 trackParam->SetBendingCoor(yRec);
431 trackParam->SetZ(zRec);
432 trackParam->SetNonBendingSlope(nonBendingSlope);
433 trackParam->SetBendingSlope(bendingSlope);
434 if (TMath::Abs(bendingMomentum) > 0)
435 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
437 trackNew->AddHitForRecAtHit(hitForRec);
438 trackNew->AddTrackParamAtHit(trackParam,0);
443 trackNew->SetTrackID(track->GetTrackID());
444 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
445 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
451 fMuonTrackRef->Clear();
452 nTrackRef = newMuonTrackRef->GetEntriesFast();
453 for (Int_t index = 0; index < nTrackRef; index++) {
454 track = (AliMUONTrack*)newMuonTrackRef->At(index);
455 AddMuonTrackReference(track);
460 newMuonTrackRef->Delete();
461 delete newMuonTrackRef;
465 //_____________________________________________________________________________
466 void AliMUONRecoCheck::ReconstructibleTracks()
468 /// Calculate the number of reconstructible tracks
471 TClonesArray* hitForRecAtHit = NULL;
472 AliMUONHitForRec* hitForRec;
474 Int_t nTrackRef, nTrackHits;
475 Int_t isChamberInTrack[10];
477 Bool_t isTrackOK = kTRUE;
479 fReconstructibleTracks = 0;
481 nTrackRef = fMuonTrackRef->GetEntriesFast();
482 for (Int_t index = 0; index < nTrackRef; index++) {
483 track = (AliMUONTrack*)fMuonTrackRef->At(index);
484 hitForRecAtHit = track->GetHitForRecAtHit();
485 nTrackHits = hitForRecAtHit->GetEntriesFast();
486 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
488 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
489 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
490 zRec = hitForRec->GetZ();
491 iChamber = hitForRec->GetChamberNumber();
492 if (iChamber < 0 || iChamber > 10) continue;
493 isChamberInTrack[iChamber] = 1;
495 // track is reconstructible if the particle is depositing a hit
496 // in the following chamber combinations:
499 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
500 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
501 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
502 Int_t nHitsInLastStations=0;
503 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
504 if (isChamberInTrack[ch]) nHitsInLastStations++;
505 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
507 if (isTrackOK) fReconstructibleTracks++;
508 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
510 fMuonTrackRef->Compress();