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 "AliRun.h" // for gAlice
31 #include "AliLoader.h"
32 #include "AliRunLoader.h"
33 #include "AliHeader.h"
35 #include "AliTrackReference.h"
38 #include <TParticle.h>
41 ClassImp(AliMUONRecoCheck)
44 //_____________________________________________________________________________
45 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
51 fReconstructibleTracks(0),
53 fIsLoadConstructor(kTRUE)
55 /// Constructor using "galice.root",
56 /// takes care of loading/unloading Kinematics, TrackRefs and Tracks
58 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
61 fRunLoader = AliRunLoader::Open(chLoader);
63 AliError(Form("no run loader found " ));
68 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
71 fMUONData = new AliMUONData(loader,"MUON","MUON");
73 AliError(Form("no MUONData found " ));
77 fRunLoader->LoadKinematics("READ");
78 fRunLoader->LoadTrackRefs("READ");
79 loader->LoadTracks("READ");
83 //_____________________________________________________________________________
84 AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONData *muondata)
90 fReconstructibleTracks(0),
92 fIsLoadConstructor(kFALSE)
94 /// Constructor from AliRunLoader and AliMUONData,
95 /// does not load/unload Kinematics, TrackRefs and Tracks internally
96 /// it should be done in the execution program (e.g. see MUONRecoCheck.C)
98 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
101 fRunLoader = runloader;
103 AliError(Form("no run loader found " ));
108 fMUONData = muondata;
110 AliError(Form("no MUONData found " ));
117 //_____________________________________________________________________________
118 AliMUONRecoCheck::~AliMUONRecoCheck()
122 delete fMuonTrackRef;
124 if(fIsLoadConstructor){
125 fRunLoader->UnloadKinematics();
126 fRunLoader->UnloadTrackRefs();
127 fRunLoader->UnloadTracks();
133 //_____________________________________________________________________________
134 void AliMUONRecoCheck::MakeTrackRef()
136 /// Make reconstructible tracks
138 AliTrackReference *trackReference;
139 AliMUONTrack *muonTrack;
140 AliMUONTrackParam *trackParam;
141 AliMUONHitForRec *hitForRec;
142 Float_t x, y, z, pX, pY, pZ, pYZ;
143 Int_t track, trackSave;
145 Float_t bendingSlope = 0;
146 Float_t nonBendingSlope = 0;
147 Float_t inverseBendingMomentum = 0;
149 TTree* treeTR = fRunLoader->TreeTR();
150 if (treeTR == NULL) return;
152 TBranch* branch = treeTR->GetBranch("MUON");
153 if (branch == NULL) return;
155 TClonesArray* trackRefs = 0;
156 branch->SetAddress(&trackRefs);
157 branch->SetAutoDelete(kTRUE);
159 Int_t nTrackRef = (Int_t)branch->GetEntries();
161 track = trackSave = -999;
163 Int_t iHitMin, iChamber, detElemId;
165 trackParam = new AliMUONTrackParam();
166 hitForRec = new AliMUONHitForRec();
168 Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
169 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
171 branch->GetEntry(iTrackRef);
176 if (!trackRefs->GetEntries()) continue;
180 muonTrack = new AliMUONTrack();
182 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
184 trackReference = (AliTrackReference*)trackRefs->At(iHit);
185 x = trackReference->X();
186 y = trackReference->Y();
187 z = trackReference->Z();
188 pX = trackReference->Px();
189 pY = trackReference->Py();
190 pZ = trackReference->Pz();
192 track = trackReference->GetTrack();
196 << "Track ID " << track
197 << " larger than max number of particles " << max << endl;
201 if (track != trackSave && iHit != 0) {
207 // track parameters at hit
208 trackParam->SetBendingCoor(y);
209 trackParam->SetNonBendingCoor(x);
212 if (TMath::Abs(pZ) > 0) {
213 bendingSlope = pY/pZ;
214 nonBendingSlope = pX/pZ;
216 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
217 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
219 trackParam->SetBendingSlope(bendingSlope);
220 trackParam->SetNonBendingSlope(nonBendingSlope);
221 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
223 hitForRec->SetBendingCoor(y);
224 hitForRec->SetNonBendingCoor(x);
226 hitForRec->SetBendingReso2(0.0);
227 hitForRec->SetNonBendingReso2(0.0);
228 detElemId = hitForRec->GetDetElemId();
229 if (detElemId) iChamber = detElemId / 100 - 1;
230 else iChamber = AliMUONConstants::ChamberNumber(z);
231 hitForRec->SetChamberNumber(iChamber);
233 muonTrack->AddTrackParamAtHit(trackParam,0);
234 muonTrack->AddHitForRecAtHit(hitForRec);
235 muonTrack->SetTrackID(track);
238 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
241 // track parameters at vertex
242 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
253 trackParam->SetBendingCoor(y);
254 trackParam->SetNonBendingCoor(x);
256 if (TMath::Abs(pZ) > 0) {
257 bendingSlope = pY/pZ;
258 nonBendingSlope = pX/pZ;
260 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
261 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
262 trackParam->SetBendingSlope(bendingSlope);
263 trackParam->SetNonBendingSlope(nonBendingSlope);
264 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
266 muonTrack->SetTrackParamAtVertex(trackParam);
269 AddMuonTrackReference(muonTrack);
273 } // end while isNewTrack
278 ReconstructibleTracks();
285 //____________________________________________________________________________
286 TClonesArray* AliMUONRecoCheck::GetTrackReco()
288 /// Return TClonesArray of reconstructed tracks
290 fMUONData->ResetRecTracks();
291 fMUONData->SetTreeAddress("RT");
292 fTrackReco = fMUONData->RecTracks();
293 fMUONData->GetRecTracks();
294 fRecoTracks = fTrackReco->GetEntriesFast();
297 //_____________________________________________________________________________
298 void AliMUONRecoCheck::PrintEvent() const
303 AliMUONHitForRec *hitForRec;
304 TClonesArray * hitForRecAtHit = 0;
305 Float_t xRec,yRec,zRec;
307 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
309 printf(" ******************************* \n");
310 printf(" nb of tracks %d \n",nTrackRef);
312 for (Int_t index = 0; index < nTrackRef; index++) {
313 track = (AliMUONTrack*)fMuonTrackRef->At(index);
314 hitForRecAtHit = track->GetHitForRecAtHit();
315 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
316 printf(" track number %d \n",index);
317 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
318 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
319 xRec = hitForRec->GetNonBendingCoor();
320 yRec = hitForRec->GetBendingCoor();
321 zRec = hitForRec->GetZ();
322 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
327 //_____________________________________________________________________________
328 void AliMUONRecoCheck::ResetTracks() const
332 if (fMuonTrackRef) fMuonTrackRef->Clear();
334 //_____________________________________________________________________________
335 void AliMUONRecoCheck::CleanMuonTrackRef()
337 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
338 /// each chamber (one when particle is entering + one when particle is leaving
339 /// the sensitive volume)
341 Float_t maxGasGap = 1.; // cm
342 AliMUONTrack *track, *trackNew;
343 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
344 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
345 TClonesArray * hitForRecAtHit = 0;
346 TClonesArray * trackParamAtHit = 0;
347 Float_t xRec,yRec,zRec;
348 Float_t xRec1,yRec1,zRec1;
349 Float_t xRec2,yRec2,zRec2;
350 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
351 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
352 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
353 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
355 Int_t iChamber = 0, detElemId = 0;
357 Int_t nTrackHits = 0;
359 hitForRec = new AliMUONHitForRec();
360 trackParam = new AliMUONTrackParam();
362 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
363 for (Int_t index = 0; index < nTrackRef; index++) {
364 track = (AliMUONTrack*)fMuonTrackRef->At(index);
365 hitForRecAtHit = track->GetHitForRecAtHit();
366 trackParamAtHit = track->GetTrackParamAtHit();
367 trackParamAtVertex = track->GetTrackParamAtVertex();
368 nTrackHits = hitForRecAtHit->GetEntriesFast();
369 trackNew = new AliMUONTrack();
371 while (iHit1 < nTrackHits) {
372 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
373 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
374 xRec1 = hitForRec1->GetNonBendingCoor();
375 yRec1 = hitForRec1->GetBendingCoor();
376 zRec1 = hitForRec1->GetZ();
380 bendingSlope1 = trackParam1->GetBendingSlope();
381 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
382 bendingMomentum1 = 0;
383 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
384 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
385 bendingSlope = bendingSlope1;
386 nonBendingSlope = nonBendingSlope1;
387 bendingMomentum = bendingMomentum1;
389 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
390 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
391 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
392 xRec2 = hitForRec2->GetNonBendingCoor();
393 yRec2 = hitForRec2->GetBendingCoor();
394 zRec2 = hitForRec2->GetZ();
395 bendingSlope2 = trackParam2->GetBendingSlope();
396 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
397 bendingMomentum2 = 0;
398 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
399 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
401 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
407 bendingSlope += bendingSlope2;
408 nonBendingSlope += nonBendingSlope2;
409 bendingMomentum += bendingMomentum2;
414 xRec /= (Float_t)nRec;
415 yRec /= (Float_t)nRec;
416 zRec /= (Float_t)nRec;
417 bendingSlope /= (Float_t)nRec;
418 nonBendingSlope /= (Float_t)nRec;
419 bendingMomentum /= (Float_t)nRec;
421 hitForRec->SetNonBendingCoor(xRec);
422 hitForRec->SetBendingCoor(yRec);
423 hitForRec->SetZ(zRec);
424 detElemId = hitForRec->GetDetElemId();
425 if (detElemId) iChamber = detElemId / 100 - 1;
426 else iChamber = AliMUONConstants::ChamberNumber(zRec);
427 hitForRec->SetChamberNumber(iChamber);
428 hitForRec->SetBendingReso2(0.0);
429 hitForRec->SetNonBendingReso2(0.0);
430 trackParam->SetNonBendingCoor(xRec);
431 trackParam->SetBendingCoor(yRec);
432 trackParam->SetZ(zRec);
433 trackParam->SetNonBendingSlope(nonBendingSlope);
434 trackParam->SetBendingSlope(bendingSlope);
435 if (TMath::Abs(bendingMomentum) > 0)
436 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
438 trackNew->AddHitForRecAtHit(hitForRec);
439 trackNew->AddTrackParamAtHit(trackParam,0);
444 trackNew->SetTrackID(track->GetTrackID());
445 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
446 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
452 fMuonTrackRef->Clear();
453 nTrackRef = newMuonTrackRef->GetEntriesFast();
454 for (Int_t index = 0; index < nTrackRef; index++) {
455 track = (AliMUONTrack*)newMuonTrackRef->At(index);
456 AddMuonTrackReference(track);
461 newMuonTrackRef->Delete();
462 delete newMuonTrackRef;
466 //_____________________________________________________________________________
467 void AliMUONRecoCheck::ReconstructibleTracks()
469 /// Calculate the number of reconstructible tracks
472 TClonesArray* hitForRecAtHit = NULL;
473 AliMUONHitForRec* hitForRec;
475 Int_t nTrackRef, nTrackHits;
476 Int_t isChamberInTrack[10];
478 Bool_t isTrackOK = kTRUE;
480 fReconstructibleTracks = 0;
482 nTrackRef = fMuonTrackRef->GetEntriesFast();
483 for (Int_t index = 0; index < nTrackRef; index++) {
484 track = (AliMUONTrack*)fMuonTrackRef->At(index);
485 hitForRecAtHit = track->GetHitForRecAtHit();
486 nTrackHits = hitForRecAtHit->GetEntriesFast();
487 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
489 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
490 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
491 zRec = hitForRec->GetZ();
492 iChamber = hitForRec->GetChamberNumber();
493 if (iChamber < 0 || iChamber > 10) continue;
494 isChamberInTrack[iChamber] = 1;
496 // track is reconstructible if the particle is depositing a hit
497 // in the following chamber combinations:
500 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
501 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
502 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
503 Int_t nHitsInLastStations=0;
504 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
505 if (isChamberInTrack[ch]) nHitsInLastStations++;
506 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
508 if (isTrackOK) fReconstructibleTracks++;
509 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
511 fMuonTrackRef->Compress();