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
35 #include "AliMUONRecoCheck.h"
36 #include "AliMUONTrack.h"
37 #include "AliMUONVTrackStore.h"
38 #include "AliMUONVCluster.h"
39 #include "AliMUONVClusterStore.h"
40 #include "AliMUONConstants.h"
41 #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(Char_t *esdFileName, Char_t *pathSim)
64 fMCEventHandler(new AliMCEventHandler()),
65 fESDEvent(new AliESDEvent()),
70 fRecoTrackRefStore(0x0),
75 // TrackRefs and Particules
76 fMCEventHandler->SetInputPath(pathSim);
77 fMCEventHandler->InitIO("");
80 fESDFile = TFile::Open(esdFileName); // open the file
81 if (!fESDFile || !fESDFile->IsOpen()) {
82 AliError(Form("opening ESD file %s failed", esdFileName));
86 fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
88 AliError("no ESD tree found");
93 fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
96 //_____________________________________________________________________________
97 AliMUONRecoCheck::~AliMUONRecoCheck()
100 delete fMCEventHandler;
102 if (fESDFile) fESDFile->Close();
106 //_____________________________________________________________________________
107 void AliMUONRecoCheck::ResetStores()
109 /// Deletes all the store objects that have been created and resets the pointers to 0x0
110 delete fTrackRefStore; fTrackRefStore = 0x0;
111 delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
112 delete fRecoTrackStore; fRecoTrackStore = 0x0;
115 //_____________________________________________________________________________
116 Int_t AliMUONRecoCheck::NumberOfEvents() const
118 /// Return the number of events
119 if (fESDTree) return fESDTree->GetEntries();
123 //_____________________________________________________________________________
124 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
126 /// Return a track store containing the reconstructed tracks (converted into
127 /// MUONTrack objects) for a given event
128 if (event != fCurrentEvent) {
130 fCurrentEvent = event;
133 if (fRecoTrackStore != 0x0) return fRecoTrackStore;
135 if (!fESDTree) return 0x0;
136 if (fESDTree->GetEvent(event) <= 0) {
137 AliError(Form("fails to read ESD object for event %d", event));
140 MakeReconstructedTracks();
141 return fRecoTrackStore;
145 //_____________________________________________________________________________
146 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
148 /// Return a track store containing the track references (converted into
149 /// MUONTrack objects) for a given event
150 if (event != fCurrentEvent) {
152 fCurrentEvent = event;
155 if (fTrackRefStore != 0x0) return fTrackRefStore;
157 if (!fMCEventHandler->GetEvent(event)) {
158 AliError(Form("fails to read MC objects for event %d", event));
162 return fTrackRefStore;
166 //_____________________________________________________________________________
167 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
169 /// Return a track store containing the reconstructible tracks for a given event
170 if (event != fCurrentEvent) {
172 fCurrentEvent = event;
175 if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
177 if (TrackRefs(event) == 0x0) return 0x0;
178 MakeReconstructibleTracks();
179 return fRecoTrackRefStore;
183 //_____________________________________________________________________________
184 void AliMUONRecoCheck::MakeReconstructedTracks()
186 /// Make reconstructed tracks
187 if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
189 // loop over all reconstructed tracks and add them to the store
190 Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
191 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
192 AliMUONESDInterface::Add(*(fESDEvent->GetMuonTrack(iTrack)), *fRecoTrackStore);
196 //_____________________________________________________________________________
197 void AliMUONRecoCheck::MakeTrackRefs()
199 /// Make reconstructible tracks
200 AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
201 if (!tmpTrackRefStore) return;
203 Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
205 TClonesArray* trackRefs;
206 Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
207 AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
209 AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
211 // loop over simulated tracks
212 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
213 Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
215 // skip empty trackRefs
216 if (nHits < 1) continue;
218 // get the particle charge for further calculation
219 TParticlePDG* ppdg = particle->GetPDG();
220 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
224 // loop over simulated track hits
225 for (Int_t iHit = 0; iHit < nHits; ++iHit) {
226 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
228 // skip trackRefs not in MUON
229 if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
231 // Get track parameters of current hit
232 x = trackReference->X();
233 y = trackReference->Y();
234 z = trackReference->Z();
235 pX = trackReference->Px();
236 pY = trackReference->Py();
237 pZ = trackReference->Pz();
239 // check chamberId of current trackReference
240 Int_t chamberId = AliMUONConstants::ChamberNumber(z);
241 if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
243 // set hit parameters
244 hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, 0, 0));
246 hit->SetErrXY(0.,0.);
248 // compute track parameters at hit
249 Double_t bendingSlope = 0;
250 Double_t nonBendingSlope = 0;
251 Double_t inverseBendingMomentum = 0;
252 if (TMath::Abs(pZ) > 0) {
253 bendingSlope = pY/pZ;
254 nonBendingSlope = pX/pZ;
256 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
257 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
258 inverseBendingMomentum *= charge;
260 // set track parameters at hit
261 AliMUONTrackParam trackParam;
262 trackParam.SetNonBendingCoor(x);
263 trackParam.SetBendingCoor(y);
265 trackParam.SetBendingSlope(bendingSlope);
266 trackParam.SetNonBendingSlope(nonBendingSlope);
267 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
269 // add track parameters at current hit to the track
270 track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
273 // if none of the track hits was in MUON, goto the next track
274 if (track.GetNClusters() < 1) continue;
276 // get track parameters at particle's vertex
284 // compute rest of track parameters at particle's vertex
287 inverseBendingMomentum = 0;
288 if (TMath::Abs(pZ) > 0) {
289 bendingSlope = pY/pZ;
290 nonBendingSlope = pX/pZ;
292 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
293 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
294 inverseBendingMomentum *= charge;
296 // set track parameters at particle's vertex
297 AliMUONTrackParam trackParamAtVertex;
298 trackParamAtVertex.SetNonBendingCoor(x);
299 trackParamAtVertex.SetBendingCoor(y);
300 trackParamAtVertex.SetZ(z);
301 trackParamAtVertex.SetBendingSlope(bendingSlope);
302 trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
303 trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
305 // add track parameters at vertex
306 track.SetTrackParamAtVertex(&trackParamAtVertex);
309 track.SetTrackID(iTrackRef);
310 tmpTrackRefStore->Add(track);
313 CleanMuonTrackRef(tmpTrackRefStore);
317 delete tmpTrackRefStore;
320 //_____________________________________________________________________________
321 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
323 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
324 /// each chamber (one when particle is entering + one when particle is leaving
325 /// the sensitive volume)
326 if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
328 Double_t maxGasGap = 1.; // cm
329 Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
330 Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
331 AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
333 AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
336 TIter next(tmpTrackRefStore->CreateIterator());
338 // loop over tmpTrackRef
340 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
342 AliMUONTrack newTrack;
344 // loop over tmpTrackRef's hits
346 Int_t nTrackHits = track->GetNClusters();
347 while (iHit1 < nTrackHits) {
348 AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
350 // get track parameters at hit1
351 x1 = trackParam1->GetNonBendingCoor();
352 y1 = trackParam1->GetBendingCoor();
353 z1 = trackParam1->GetZ();
354 pX1 = trackParam1->Px();
355 pY1 = trackParam1->Py();
356 pZ1 = trackParam1->Pz();
358 // prepare new track parameters
366 // loop over next tmpTrackRef's hits
367 Int_t nCombinedHits = 1;
368 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
369 AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
371 // get z position of hit2
372 z2 = trackParam2->GetZ();
374 // complete new track parameters if hit2 is on the same detection element
375 if ( TMath::Abs(z2-z1) < maxGasGap ) {
376 x += trackParam2->GetNonBendingCoor();
377 y += trackParam2->GetBendingCoor();
379 pX += trackParam2->Px();
380 pY += trackParam2->Py();
381 pZ += trackParam2->Pz();
388 // finalize new track parameters
389 x /= (Double_t)nCombinedHits;
390 y /= (Double_t)nCombinedHits;
391 z /= (Double_t)nCombinedHits;
392 pX /= (Double_t)nCombinedHits;
393 pY /= (Double_t)nCombinedHits;
394 pZ /= (Double_t)nCombinedHits;
397 inverseBendingMomentum = 0;
398 if (TMath::Abs(pZ) > 0) {
399 bendingSlope = pY/pZ;
400 nonBendingSlope = pX/pZ;
402 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
403 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
404 inverseBendingMomentum *= trackParam1->GetCharge();
406 // set hit parameters
407 hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0));
409 hit->SetErrXY(0.,0.);
411 // set new track parameters at new hit
412 AliMUONTrackParam trackParam;
413 trackParam.SetNonBendingCoor(x);
414 trackParam.SetBendingCoor(y);
416 trackParam.SetBendingSlope(bendingSlope);
417 trackParam.SetNonBendingSlope(nonBendingSlope);
418 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
420 // add track parameters at current hit to the track
421 newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
426 newTrack.SetTrackID(track->GetTrackID());
427 newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
428 fTrackRefStore->Add(newTrack);
436 //_____________________________________________________________________________
437 void AliMUONRecoCheck::MakeReconstructibleTracks()
439 /// Isolate the reconstructible tracks
440 if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
442 // create iterator on trackRef
443 TIter next(fTrackRefStore->CreateIterator());
445 // loop over trackRef
447 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
449 Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
450 for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
452 // loop over trackRef's hits to get hit chambers
453 Int_t nTrackHits = track->GetNClusters();
454 for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
455 AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
456 chamberInTrack[hit->GetChamberId()] = kTRUE;
459 // track is reconstructible if the particle is depositing a hit
460 // in the following chamber combinations:
461 Bool_t trackOK = kTRUE;
462 if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
463 if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
464 if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
465 Int_t nHitsInLastStations = 0;
466 for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
467 if (chamberInTrack[iCh]) nHitsInLastStations++;
468 if(nHitsInLastStations < 3) trackOK = kFALSE;
470 // Add reconstructible tracks to fRecoTrackRefStore
471 if (trackOK) fRecoTrackRefStore->Add(*track);
473 delete [] chamberInTrack;