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 /// 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.
24 #include "AliMUONRecoCheck.h"
25 #include "AliMUONHitForRec.h"
26 #include "AliMUONTrack.h"
27 #include "AliMUONConstants.h"
28 #include "AliMUONMCDataInterface.h"
29 #include "AliMUONDataInterface.h"
31 #include "AliTrackReference.h"
33 #include "AliMUONTrackStoreV1.h"
34 #include <Riostream.h>
35 #include <TParticle.h>
36 #include <TParticlePDG.h>
39 ClassImp(AliMUONRecoCheck)
42 //_____________________________________________________________________________
43 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
45 fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
46 fDataInterface(new AliMUONDataInterface(chLoader))
51 //_____________________________________________________________________________
52 AliMUONRecoCheck::~AliMUONRecoCheck()
55 delete fMCDataInterface;
56 delete fDataInterface;
59 //_____________________________________________________________________________
61 AliMUONRecoCheck::ReconstructedTracks(Int_t event)
63 /// Return the reconstructed track store for a given event
64 return fDataInterface->TrackStore(event);
67 //_____________________________________________________________________________
69 AliMUONRecoCheck::TrackRefs(Int_t event)
71 /// Return a track store containing the track references (converted into
72 /// MUONTrack objects) for a given event
73 return MakeTrackRefs(event);
76 //_____________________________________________________________________________
78 AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
80 /// Return a track store containing the reconstructible tracks for a given event
81 AliMUONVTrackStore* tmp = MakeTrackRefs(event);
82 AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
84 return reconstructible;
87 //_____________________________________________________________________________
89 AliMUONRecoCheck::MakeTrackRefs(Int_t event)
91 /// Make reconstructible tracks
93 AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
95 Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
97 Int_t trackSave(-999);
102 AliStack* stack = fMCDataInterface->Stack(event);
103 Int_t max = stack->GetNtrack();
105 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef)
107 TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
112 if (!trackRefs->GetEntries()) continue;
116 AliMUONTrack muonTrack;
118 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit)
120 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
122 Float_t x = trackReference->X();
123 Float_t y = trackReference->Y();
124 Float_t z = trackReference->Z();
125 Float_t pX = trackReference->Px();
126 Float_t pY = trackReference->Py();
127 Float_t pZ = trackReference->Pz();
129 Int_t track = trackReference->GetTrack();
134 << "Track ID " << track
135 << " larger than max number of particles " << max << endl;
139 if (track != trackSave && iHit != 0) {
145 Float_t bendingSlope = 0;
146 Float_t nonBendingSlope = 0;
147 Float_t inverseBendingMomentum = 0;
149 AliMUONTrackParam trackParam;
151 // track parameters at hit
152 trackParam.SetBendingCoor(y);
153 trackParam.SetNonBendingCoor(x);
156 if (TMath::Abs(pZ) > 0)
158 bendingSlope = pY/pZ;
159 nonBendingSlope = pX/pZ;
161 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
162 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
164 trackParam.SetBendingSlope(bendingSlope);
165 trackParam.SetNonBendingSlope(nonBendingSlope);
166 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
168 AliMUONHitForRec hitForRec;
170 hitForRec.SetBendingCoor(y);
171 hitForRec.SetNonBendingCoor(x);
173 hitForRec.SetBendingReso2(0.0);
174 hitForRec.SetNonBendingReso2(0.0);
175 Int_t detElemId = hitForRec.GetDetElemId();
179 iChamber = detElemId / 100 - 1;
183 iChamber = AliMUONConstants::ChamberNumber(z);
185 hitForRec.SetChamberNumber(iChamber);
187 muonTrack.AddTrackParamAtHit(&trackParam,0);
188 muonTrack.AddHitForRecAtHit(&hitForRec);
189 muonTrack.SetTrackID(track);
192 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
195 // track parameters at vertex
196 TParticle* particle = stack->Particle(muonTrack.GetTrackID());
200 Float_t x = particle->Vx();
201 Float_t y = particle->Vy();
202 Float_t z = particle->Vz();
203 Float_t pX = particle->Px();
204 Float_t pY = particle->Py();
205 Float_t pZ = particle->Pz();
207 AliMUONTrackParam trackParam;
209 trackParam.SetBendingCoor(y);
210 trackParam.SetNonBendingCoor(x);
213 Float_t bendingSlope = 0;
214 Float_t nonBendingSlope = 0;
215 Float_t inverseBendingMomentum = 0;
217 if (TMath::Abs(pZ) > 0)
219 bendingSlope = pY/pZ;
220 nonBendingSlope = pX/pZ;
223 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
224 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
226 TParticlePDG* ppdg = particle->GetPDG(1);
227 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
228 inverseBendingMomentum *= charge;
230 trackParam.SetBendingSlope(bendingSlope);
231 trackParam.SetNonBendingSlope(nonBendingSlope);
232 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
234 muonTrack.SetTrackParamAtVertex(&trackParam);
237 trackRefStore->Add(muonTrack);
238 } // end while isNewTrack
242 AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
244 delete trackRefStore;
249 //_____________________________________________________________________________
251 AliMUONRecoCheck::NumberOfEvents() const
253 /// Return the number of events
254 if ( fDataInterface )
256 return fDataInterface->NumberOfEvents();
261 //_____________________________________________________________________________
263 AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
265 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
266 /// each chamber (one when particle is entering + one when particle is leaving
267 /// the sensitive volume)
269 Float_t maxGasGap = 1.; // cm
270 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
271 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
272 TClonesArray * hitForRecAtHit = 0;
273 TClonesArray * trackParamAtHit = 0;
274 Float_t xRec,yRec,zRec;
275 Float_t xRec1,yRec1,zRec1;
276 Float_t xRec2,yRec2,zRec2;
277 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
278 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
279 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
281 AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
283 Int_t iChamber = 0, detElemId = 0;
285 Int_t nTrackHits = 0;
287 hitForRec = new AliMUONHitForRec();
288 trackParam = new AliMUONTrackParam();
290 TIter next(trackRefs.CreateIterator());
293 while ( ( track = static_cast<AliMUONTrack*>(next())) )
295 hitForRecAtHit = track->GetHitForRecAtHit();
296 trackParamAtHit = track->GetTrackParamAtHit();
297 trackParamAtVertex = track->GetTrackParamAtVertex();
298 nTrackHits = hitForRecAtHit->GetEntriesFast();
299 AliMUONTrack trackNew;
301 while (iHit1 < nTrackHits)
303 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
304 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
305 xRec1 = hitForRec1->GetNonBendingCoor();
306 yRec1 = hitForRec1->GetBendingCoor();
307 zRec1 = hitForRec1->GetZ();
311 bendingSlope1 = trackParam1->GetBendingSlope();
312 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
313 bendingMomentum1 = 0;
314 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
315 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
316 bendingSlope = bendingSlope1;
317 nonBendingSlope = nonBendingSlope1;
318 bendingMomentum = bendingMomentum1;
320 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
322 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
323 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
324 xRec2 = hitForRec2->GetNonBendingCoor();
325 yRec2 = hitForRec2->GetBendingCoor();
326 zRec2 = hitForRec2->GetZ();
327 bendingSlope2 = trackParam2->GetBendingSlope();
328 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
329 bendingMomentum2 = 0;
330 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
331 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
333 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
339 bendingSlope += bendingSlope2;
340 nonBendingSlope += nonBendingSlope2;
341 bendingMomentum += bendingMomentum2;
346 xRec /= (Float_t)nRec;
347 yRec /= (Float_t)nRec;
348 zRec /= (Float_t)nRec;
349 bendingSlope /= (Float_t)nRec;
350 nonBendingSlope /= (Float_t)nRec;
351 bendingMomentum /= (Float_t)nRec;
353 hitForRec->SetNonBendingCoor(xRec);
354 hitForRec->SetBendingCoor(yRec);
355 hitForRec->SetZ(zRec);
356 detElemId = hitForRec->GetDetElemId();
357 if (detElemId) iChamber = detElemId / 100 - 1;
358 else iChamber = AliMUONConstants::ChamberNumber(zRec);
359 hitForRec->SetChamberNumber(iChamber);
360 hitForRec->SetBendingReso2(0.0);
361 hitForRec->SetNonBendingReso2(0.0);
362 trackParam->SetNonBendingCoor(xRec);
363 trackParam->SetBendingCoor(yRec);
364 trackParam->SetZ(zRec);
365 trackParam->SetNonBendingSlope(nonBendingSlope);
366 trackParam->SetBendingSlope(bendingSlope);
367 if (TMath::Abs(bendingMomentum) > 0)
368 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
370 trackNew.AddHitForRecAtHit(hitForRec);
371 trackNew.AddTrackParamAtHit(trackParam,0);
376 trackNew.SetTrackID(track->GetTrackID());
377 trackNew.SetTrackParamAtVertex(trackParamAtVertex);
378 newMuonTrackRef->Add(trackNew);
384 return newMuonTrackRef;
387 //_____________________________________________________________________________
389 AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
391 /// Calculate the number of reconstructible tracks
393 AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
395 TClonesArray* hitForRecAtHit = NULL;
396 AliMUONHitForRec* hitForRec;
399 Int_t isChamberInTrack[10];
401 Bool_t isTrackOK = kTRUE;
403 TIter next(trackRefs.CreateIterator());
406 while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
408 hitForRecAtHit = track->GetHitForRecAtHit();
409 nTrackHits = hitForRecAtHit->GetEntriesFast();
410 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
412 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
413 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
414 zRec = hitForRec->GetZ();
415 iChamber = hitForRec->GetChamberNumber();
416 if (iChamber < 0 || iChamber > 10) continue;
417 isChamberInTrack[iChamber] = 1;
419 // track is reconstructible if the particle is depositing a hit
420 // in the following chamber combinations:
423 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
424 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
425 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
426 Int_t nHitsInLastStations=0;
427 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
428 if (isChamberInTrack[ch]) nHitsInLastStations++;
429 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
433 reconstructibleStore->Add(*track);
437 return reconstructibleStore;