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 //-----------------------------------------------------------------------------
26 #include "AliMUONRecoCheck.h"
27 #include "AliMUONHitForRec.h"
28 #include "AliMUONTrack.h"
29 #include "AliMUONConstants.h"
30 #include "AliMUONMCDataInterface.h"
31 #include "AliMUONDataInterface.h"
33 #include "AliTrackReference.h"
35 #include "AliMUONTrackStoreV1.h"
36 #include <Riostream.h>
37 #include <TParticle.h>
38 #include <TParticlePDG.h>
41 ClassImp(AliMUONRecoCheck)
44 //_____________________________________________________________________________
45 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
47 fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
48 fDataInterface(new AliMUONDataInterface(chLoader))
53 //_____________________________________________________________________________
54 AliMUONRecoCheck::~AliMUONRecoCheck()
57 delete fMCDataInterface;
58 delete fDataInterface;
61 //_____________________________________________________________________________
63 AliMUONRecoCheck::ReconstructedTracks(Int_t event)
65 /// Return the reconstructed track store for a given event
66 return fDataInterface->TrackStore(event);
69 //_____________________________________________________________________________
71 AliMUONRecoCheck::TrackRefs(Int_t event)
73 /// Return a track store containing the track references (converted into
74 /// MUONTrack objects) for a given event
75 return MakeTrackRefs(event);
78 //_____________________________________________________________________________
80 AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
82 /// Return a track store containing the reconstructible tracks for a given event
83 AliMUONVTrackStore* tmp = MakeTrackRefs(event);
84 AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
86 return reconstructible;
89 //_____________________________________________________________________________
91 AliMUONRecoCheck::MakeTrackRefs(Int_t event)
93 /// Make reconstructible tracks
95 AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
97 Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
99 Int_t trackSave(-999);
104 AliStack* stack = fMCDataInterface->Stack(event);
105 Int_t max = stack->GetNtrack();
107 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef)
109 TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
114 if (!trackRefs->GetEntries()) continue;
118 AliMUONTrack muonTrack;
120 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit)
122 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
124 Float_t x = trackReference->X();
125 Float_t y = trackReference->Y();
126 Float_t z = trackReference->Z();
127 Float_t pX = trackReference->Px();
128 Float_t pY = trackReference->Py();
129 Float_t pZ = trackReference->Pz();
131 Int_t track = trackReference->GetTrack();
136 << "Track ID " << track
137 << " larger than max number of particles " << max << endl;
141 if (track != trackSave && iHit != 0) {
147 Float_t bendingSlope = 0;
148 Float_t nonBendingSlope = 0;
149 Float_t inverseBendingMomentum = 0;
151 AliMUONTrackParam trackParam;
153 // track parameters at hit
154 trackParam.SetBendingCoor(y);
155 trackParam.SetNonBendingCoor(x);
158 if (TMath::Abs(pZ) > 0)
160 bendingSlope = pY/pZ;
161 nonBendingSlope = pX/pZ;
163 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
164 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
166 trackParam.SetBendingSlope(bendingSlope);
167 trackParam.SetNonBendingSlope(nonBendingSlope);
168 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
170 AliMUONHitForRec hitForRec;
172 hitForRec.SetBendingCoor(y);
173 hitForRec.SetNonBendingCoor(x);
175 hitForRec.SetBendingReso2(0.0);
176 hitForRec.SetNonBendingReso2(0.0);
177 Int_t detElemId = hitForRec.GetDetElemId();
181 iChamber = detElemId / 100 - 1;
185 iChamber = AliMUONConstants::ChamberNumber(z);
187 hitForRec.SetChamberNumber(iChamber);
189 muonTrack.AddTrackParamAtHit(&trackParam,0);
190 muonTrack.AddHitForRecAtHit(&hitForRec);
191 muonTrack.SetTrackID(track);
194 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
197 // track parameters at vertex
198 TParticle* particle = stack->Particle(muonTrack.GetTrackID());
202 Float_t x = particle->Vx();
203 Float_t y = particle->Vy();
204 Float_t z = particle->Vz();
205 Float_t pX = particle->Px();
206 Float_t pY = particle->Py();
207 Float_t pZ = particle->Pz();
209 AliMUONTrackParam trackParam;
211 trackParam.SetBendingCoor(y);
212 trackParam.SetNonBendingCoor(x);
215 Float_t bendingSlope = 0;
216 Float_t nonBendingSlope = 0;
217 Float_t inverseBendingMomentum = 0;
219 if (TMath::Abs(pZ) > 0)
221 bendingSlope = pY/pZ;
222 nonBendingSlope = pX/pZ;
225 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
226 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
228 TParticlePDG* ppdg = particle->GetPDG(1);
229 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
230 inverseBendingMomentum *= charge;
232 trackParam.SetBendingSlope(bendingSlope);
233 trackParam.SetNonBendingSlope(nonBendingSlope);
234 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
236 muonTrack.SetTrackParamAtVertex(&trackParam);
239 trackRefStore->Add(muonTrack);
240 } // end while isNewTrack
243 AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
245 delete trackRefStore;
250 //_____________________________________________________________________________
252 AliMUONRecoCheck::NumberOfEvents() const
254 /// Return the number of events
255 if ( fDataInterface )
257 return fDataInterface->NumberOfEvents();
262 //_____________________________________________________________________________
264 AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
266 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
267 /// each chamber (one when particle is entering + one when particle is leaving
268 /// the sensitive volume)
270 Float_t maxGasGap = 1.; // cm
271 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
272 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
273 TClonesArray * hitForRecAtHit = 0;
274 TClonesArray * trackParamAtHit = 0;
275 Float_t xRec,yRec,zRec;
276 Float_t xRec1,yRec1,zRec1;
277 Float_t xRec2,yRec2,zRec2;
278 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
279 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
280 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
282 AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
284 Int_t iChamber = 0, detElemId = 0;
286 Int_t nTrackHits = 0;
288 hitForRec = new AliMUONHitForRec();
289 trackParam = new AliMUONTrackParam();
291 TIter next(trackRefs.CreateIterator());
294 while ( ( track = static_cast<AliMUONTrack*>(next())) )
296 hitForRecAtHit = track->GetHitForRecAtHit();
297 trackParamAtHit = track->GetTrackParamAtHit();
298 trackParamAtVertex = track->GetTrackParamAtVertex();
299 nTrackHits = hitForRecAtHit->GetEntriesFast();
300 AliMUONTrack trackNew;
302 while (iHit1 < nTrackHits)
304 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
305 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
306 xRec1 = hitForRec1->GetNonBendingCoor();
307 yRec1 = hitForRec1->GetBendingCoor();
308 zRec1 = hitForRec1->GetZ();
312 bendingSlope1 = trackParam1->GetBendingSlope();
313 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
314 bendingMomentum1 = 0;
315 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
316 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
317 bendingSlope = bendingSlope1;
318 nonBendingSlope = nonBendingSlope1;
319 bendingMomentum = bendingMomentum1;
321 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
323 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
324 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
325 xRec2 = hitForRec2->GetNonBendingCoor();
326 yRec2 = hitForRec2->GetBendingCoor();
327 zRec2 = hitForRec2->GetZ();
328 bendingSlope2 = trackParam2->GetBendingSlope();
329 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
330 bendingMomentum2 = 0;
331 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
332 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
334 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
340 bendingSlope += bendingSlope2;
341 nonBendingSlope += nonBendingSlope2;
342 bendingMomentum += bendingMomentum2;
347 xRec /= (Float_t)nRec;
348 yRec /= (Float_t)nRec;
349 zRec /= (Float_t)nRec;
350 bendingSlope /= (Float_t)nRec;
351 nonBendingSlope /= (Float_t)nRec;
352 bendingMomentum /= (Float_t)nRec;
354 hitForRec->SetNonBendingCoor(xRec);
355 hitForRec->SetBendingCoor(yRec);
356 hitForRec->SetZ(zRec);
357 detElemId = hitForRec->GetDetElemId();
358 if (detElemId) iChamber = detElemId / 100 - 1;
359 else iChamber = AliMUONConstants::ChamberNumber(zRec);
360 hitForRec->SetChamberNumber(iChamber);
361 hitForRec->SetBendingReso2(0.0);
362 hitForRec->SetNonBendingReso2(0.0);
363 trackParam->SetNonBendingCoor(xRec);
364 trackParam->SetBendingCoor(yRec);
365 trackParam->SetZ(zRec);
366 trackParam->SetNonBendingSlope(nonBendingSlope);
367 trackParam->SetBendingSlope(bendingSlope);
368 if (TMath::Abs(bendingMomentum) > 0)
369 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
371 trackNew.AddHitForRecAtHit(hitForRec);
372 trackNew.AddTrackParamAtHit(trackParam,0);
377 trackNew.SetTrackID(track->GetTrackID());
378 trackNew.SetTrackParamAtVertex(trackParamAtVertex);
379 newMuonTrackRef->Add(trackNew);
385 return newMuonTrackRef;
388 //_____________________________________________________________________________
390 AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
392 /// Calculate the number of reconstructible tracks
394 AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
396 TClonesArray* hitForRecAtHit = NULL;
397 AliMUONHitForRec* hitForRec;
400 Int_t isChamberInTrack[10];
402 Bool_t isTrackOK = kTRUE;
404 TIter next(trackRefs.CreateIterator());
407 while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
409 hitForRecAtHit = track->GetHitForRecAtHit();
410 nTrackHits = hitForRecAtHit->GetEntriesFast();
411 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
413 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
414 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
415 zRec = hitForRec->GetZ();
416 iChamber = hitForRec->GetChamberNumber();
417 if (iChamber < 0 || iChamber > 10) continue;
418 isChamberInTrack[iChamber] = 1;
420 // track is reconstructible if the particle is depositing a hit
421 // in the following chamber combinations:
424 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
425 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
426 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
427 Int_t nHitsInLastStations=0;
428 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
429 if (isChamberInTrack[ch]) nHitsInLastStations++;
430 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
434 reconstructibleStore->Add(*track);
438 return reconstructibleStore;