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 /// Utility class to check reconstruction
21 /// Reconstructed tracks are compared to reference tracks.
22 /// The reference tracks are built from AliTrackReference for the
23 /// hit in chamber (0..9) and from kinematics for the vertex parameters.
24 //-----------------------------------------------------------------------------
27 // Added a method to create a list of reconstructed AliMUONTrack objects from
28 // ESD data. This is necessary since the track objects that are actually created
29 // during offline reconstruction are no longer stored to disk.
30 // - Artur Szostak <artursz@iafrica.com>
32 // Use the new ESDInterface to create MUON objects from ESD data
36 #include "AliMUONRecoCheck.h"
37 #include "AliMUONTrack.h"
38 #include "AliMUONVTrackStore.h"
39 #include "AliMUONVCluster.h"
40 #include "AliMUONVClusterStore.h"
41 #include "AliMUONConstants.h"
42 #include "AliMUONESDInterface.h"
43 #include "AliMCEventHandler.h"
44 #include "AliMCEvent.h"
46 #include "AliTrackReference.h"
48 #include "AliESDEvent.h"
49 #include "AliESDMuonTrack.h"
53 #include <TParticle.h>
54 #include <TParticlePDG.h>
55 #include <Riostream.h>
58 ClassImp(AliMUONRecoCheck)
61 //_____________________________________________________________________________
62 AliMUONRecoCheck::AliMUONRecoCheck(const Char_t *esdFileName, const Char_t *pathSim)
64 fMCEventHandler(new AliMCEventHandler()),
65 fESDEvent(new AliESDEvent()),
70 fRecoTrackRefStore(0x0),
76 // TrackRefs and Particules
77 fMCEventHandler->SetInputPath(pathSim);
78 fMCEventHandler->InitIO("");
81 fESDFile = TFile::Open(esdFileName); // open the file
82 if (!fESDFile || !fESDFile->IsOpen()) {
83 AliError(Form("opening ESD file %s failed", esdFileName));
87 fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
89 AliError("no ESD tree found");
94 fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
97 //_____________________________________________________________________________
98 AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent *esdEvent, AliMCEventHandler *mcEventHandler)
106 fRecoTrackRefStore(0x0),
107 fRecoTrackStore(0x0),
108 fESDEventOwner(kFALSE)
112 // TrackRefs and Particules
113 fMCEventHandler = mcEventHandler;
116 fESDEvent = esdEvent;
120 //_____________________________________________________________________________
121 AliMUONRecoCheck::~AliMUONRecoCheck()
124 if (fESDEventOwner) {
125 delete fMCEventHandler;
127 if (fESDFile) fESDFile->Close();
132 //_____________________________________________________________________________
133 void AliMUONRecoCheck::ResetStores()
135 /// Deletes all the store objects that have been created and resets the pointers to 0x0
136 delete fTrackRefStore; fTrackRefStore = 0x0;
137 delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
138 delete fRecoTrackStore; fRecoTrackStore = 0x0;
141 //_____________________________________________________________________________
142 Int_t AliMUONRecoCheck::NumberOfEvents() const
144 /// Return the number of events
145 if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
149 //_____________________________________________________________________________
150 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
152 /// Return a track store containing the reconstructed tracks (converted into
153 /// MUONTrack objects) for a given event
155 if (!fESDEventOwner) {
156 MakeReconstructedTracks();
157 return fRecoTrackStore;
160 if (event != fCurrentEvent) {
162 fCurrentEvent = event;
165 if (fRecoTrackStore != 0x0) return fRecoTrackStore;
167 if (!fESDTree) return 0x0;
168 if (fESDTree->GetEvent(event) <= 0) {
169 AliError(Form("fails to read ESD object for event %d", event));
172 MakeReconstructedTracks();
173 return fRecoTrackStore;
177 //_____________________________________________________________________________
178 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
180 /// Return a track store containing the track references (converted into
181 /// MUONTrack objects) for a given event
183 if (!fESDEventOwner) {
185 return fTrackRefStore;
188 if (event != fCurrentEvent) {
190 fCurrentEvent = event;
193 if (fTrackRefStore != 0x0) return fTrackRefStore;
195 if (!fMCEventHandler->GetEvent(event)) {
196 AliError(Form("fails to read MC objects for event %d", event));
200 return fTrackRefStore;
204 //_____________________________________________________________________________
205 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
207 /// Return a track store containing the reconstructible tracks for a given event
209 if (!fESDEventOwner) {
210 if (TrackRefs(event) == 0x0) return 0x0;
211 MakeReconstructibleTracks();
212 return fRecoTrackRefStore;
215 if (event != fCurrentEvent) {
217 fCurrentEvent = event;
220 if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
222 if (TrackRefs(event) == 0x0) return 0x0;
223 MakeReconstructibleTracks();
224 return fRecoTrackRefStore;
228 //_____________________________________________________________________________
229 void AliMUONRecoCheck::MakeReconstructedTracks()
231 /// Make reconstructed tracks
232 if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
234 // loop over all reconstructed tracks and add them to the store (skip ghosts)
235 Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
236 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
237 AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
238 if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore);
243 //_____________________________________________________________________________
244 void AliMUONRecoCheck::MakeTrackRefs()
246 /// Make reconstructible tracks
247 AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
248 if (!tmpTrackRefStore) return;
250 Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
252 TClonesArray* trackRefs;
253 Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
254 AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
256 AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
258 // loop over simulated tracks
259 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
260 Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
262 // skip empty trackRefs
263 if (nHits < 1) continue;
265 // get the particle charge for further calculation
266 TParticlePDG* ppdg = particle->GetPDG();
267 Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
271 // loop over simulated track hits
272 for (Int_t iHit = 0; iHit < nHits; ++iHit) {
273 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
275 // skip trackRefs not in MUON
276 if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
278 // Get track parameters of current hit
279 x = trackReference->X();
280 y = trackReference->Y();
281 z = trackReference->Z();
282 pX = trackReference->Px();
283 pY = trackReference->Py();
284 pZ = trackReference->Pz();
286 // check chamberId of current trackReference
287 Int_t detElemId = trackReference->UserId();
288 Int_t chamberId = detElemId / 100 - 1;
289 if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
291 // set hit parameters
292 hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
294 hit->SetErrXY(0.,0.);
296 // compute track parameters at hit
299 inverseBendingMomentum = 0;
300 if (TMath::Abs(pZ) > 0) {
301 bendingSlope = pY/pZ;
302 nonBendingSlope = pX/pZ;
304 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
305 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
306 inverseBendingMomentum *= charge;
308 // set track parameters at hit
309 AliMUONTrackParam trackParam;
310 trackParam.SetNonBendingCoor(x);
311 trackParam.SetBendingCoor(y);
313 trackParam.SetBendingSlope(bendingSlope);
314 trackParam.SetNonBendingSlope(nonBendingSlope);
315 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
317 // add track parameters at current hit to the track
318 track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
321 // if none of the track hits was in MUON, goto the next track
322 if (track.GetNClusters() < 1) continue;
324 // get track parameters at particle's vertex
332 // compute rest of track parameters at particle's vertex
335 inverseBendingMomentum = 0;
336 if (TMath::Abs(pZ) > 0) {
337 bendingSlope = pY/pZ;
338 nonBendingSlope = pX/pZ;
340 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
341 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
342 inverseBendingMomentum *= charge;
344 // set track parameters at particle's vertex
345 AliMUONTrackParam trackParamAtVertex;
346 trackParamAtVertex.SetNonBendingCoor(x);
347 trackParamAtVertex.SetBendingCoor(y);
348 trackParamAtVertex.SetZ(z);
349 trackParamAtVertex.SetBendingSlope(bendingSlope);
350 trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
351 trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
353 // add track parameters at vertex
354 track.SetTrackParamAtVertex(&trackParamAtVertex);
357 track.SetTrackID(iTrackRef);
358 tmpTrackRefStore->Add(track);
361 CleanMuonTrackRef(tmpTrackRefStore);
365 delete tmpTrackRefStore;
368 //_____________________________________________________________________________
369 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
371 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
372 /// each chamber (one when particle is entering + one when particle is leaving
373 /// the sensitive volume)
374 if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
376 Double_t maxGasGap = 1.; // cm
377 Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
378 Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
379 AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
381 AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
384 TIter next(tmpTrackRefStore->CreateIterator());
386 // loop over tmpTrackRef
388 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
390 AliMUONTrack newTrack;
392 // loop over tmpTrackRef's hits
394 Int_t nTrackHits = track->GetNClusters();
395 while (iHit1 < nTrackHits) {
396 AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
398 // get track parameters at hit1
399 x1 = trackParam1->GetNonBendingCoor();
400 y1 = trackParam1->GetBendingCoor();
401 z1 = trackParam1->GetZ();
402 pX1 = trackParam1->Px();
403 pY1 = trackParam1->Py();
404 pZ1 = trackParam1->Pz();
406 // prepare new track parameters
414 // loop over next tmpTrackRef's hits
415 Int_t nCombinedHits = 1;
416 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
417 AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
419 // get z position of hit2
420 z2 = trackParam2->GetZ();
422 // complete new track parameters if hit2 is on the same detection element
423 if ( TMath::Abs(z2-z1) < maxGasGap ) {
424 x += trackParam2->GetNonBendingCoor();
425 y += trackParam2->GetBendingCoor();
427 pX += trackParam2->Px();
428 pY += trackParam2->Py();
429 pZ += trackParam2->Pz();
436 // finalize new track parameters
437 x /= (Double_t)nCombinedHits;
438 y /= (Double_t)nCombinedHits;
439 z /= (Double_t)nCombinedHits;
440 pX /= (Double_t)nCombinedHits;
441 pY /= (Double_t)nCombinedHits;
442 pZ /= (Double_t)nCombinedHits;
445 inverseBendingMomentum = 0;
446 if (TMath::Abs(pZ) > 0) {
447 bendingSlope = pY/pZ;
448 nonBendingSlope = pX/pZ;
450 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
451 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
452 inverseBendingMomentum *= trackParam1->GetCharge();
454 // set hit parameters
455 hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
457 hit->SetErrXY(0.,0.);
459 // set new track parameters at new hit
460 AliMUONTrackParam trackParam;
461 trackParam.SetNonBendingCoor(x);
462 trackParam.SetBendingCoor(y);
464 trackParam.SetBendingSlope(bendingSlope);
465 trackParam.SetNonBendingSlope(nonBendingSlope);
466 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
468 // add track parameters at current hit to the track
469 newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
474 newTrack.SetTrackID(track->GetTrackID());
475 newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
476 fTrackRefStore->Add(newTrack);
484 //_____________________________________________________________________________
485 void AliMUONRecoCheck::MakeReconstructibleTracks()
487 /// Isolate the reconstructible tracks
488 if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
490 // create iterator on trackRef
491 TIter next(fTrackRefStore->CreateIterator());
493 // loop over trackRef
495 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
497 Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
498 for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
500 // loop over trackRef's hits to get hit chambers
501 Int_t nTrackHits = track->GetNClusters();
502 for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
503 AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
504 chamberInTrack[hit->GetChamberId()] = kTRUE;
507 // track is reconstructible if the particle is depositing a hit
508 // in the following chamber combinations:
509 Bool_t trackOK = kTRUE;
510 if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
511 if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
512 if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
513 Int_t nHitsInLastStations = 0;
514 for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
515 if (chamberInTrack[iCh]) nHitsInLastStations++;
516 if(nHitsInLastStations < 3) trackOK = kFALSE;
518 // Add reconstructible tracks to fRecoTrackRefStore
519 if (trackOK) fRecoTrackRefStore->Add(*track);
521 delete [] chamberInTrack;