Bug fix
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
CommitLineData
b8dc484b 1/**************************************************************************
48217459 2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
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**************************************************************************/
b8dc484b 15
d19b6003 16/* $Id$ */
17
3d1463c8 18//-----------------------------------------------------------------------------
5398f946 19/// \class AliMUONRecoCheck
78649106 20/// Utility class to check reconstruction
5398f946 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.
3d1463c8 24//-----------------------------------------------------------------------------
b8dc484b 25
cc056de1 26#include "AliMUONRecoCheck.h"
b16234ab 27#include "AliMUONHitForRec.h"
cc056de1 28#include "AliMUONTrack.h"
cc056de1 29#include "AliMUONConstants.h"
48217459 30#include "AliMUONMCDataInterface.h"
31#include "AliMUONDataInterface.h"
b8dc484b 32#include "AliStack.h"
cc056de1 33#include "AliTrackReference.h"
34#include "AliLog.h"
48217459 35#include "AliMUONTrackStoreV1.h"
36#include <Riostream.h>
cc056de1 37#include <TParticle.h>
5b1afea0 38#include <TParticlePDG.h>
b8dc484b 39
5398f946 40/// \cond CLASSIMP
b8dc484b 41ClassImp(AliMUONRecoCheck)
5398f946 42/// \endcond
b8dc484b 43
44//_____________________________________________________________________________
48217459 45AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
46: TObject(),
47fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
48fDataInterface(new AliMUONDataInterface(chLoader))
068e4c36 49{
48217459 50 /// Normal ctor
068e4c36 51}
52
53//_____________________________________________________________________________
48217459 54AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b 55{
48217459 56 /// Destructor
57 delete fMCDataInterface;
58 delete fDataInterface;
b8dc484b 59}
60
48217459 61//_____________________________________________________________________________
62AliMUONVTrackStore*
63AliMUONRecoCheck::ReconstructedTracks(Int_t event)
64{
65 /// Return the reconstructed track store for a given event
66 return fDataInterface->TrackStore(event);
67}
45a05f07 68
b8dc484b 69//_____________________________________________________________________________
48217459 70AliMUONVTrackStore*
71AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b 72{
48217459 73 /// Return a track store containing the track references (converted into
74 /// MUONTrack objects) for a given event
75 return MakeTrackRefs(event);
b8dc484b 76}
77
78//_____________________________________________________________________________
48217459 79AliMUONVTrackStore*
80AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b 81{
48217459 82 /// Return a track store containing the reconstructible tracks for a given event
83 AliMUONVTrackStore* tmp = MakeTrackRefs(event);
84 AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
85 delete tmp;
86 return reconstructible;
87}
b8dc484b 88
48217459 89//_____________________________________________________________________________
90AliMUONVTrackStore*
91AliMUONRecoCheck::MakeTrackRefs(Int_t event)
92{
93 /// Make reconstructible tracks
94
95 AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
96
97 Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
98
99 Int_t trackSave(-999);
100 Int_t iHitMin(0);
b8dc484b 101
b8dc484b 102 Bool_t isNewTrack;
48217459 103
104 AliStack* stack = fMCDataInterface->Stack(event);
105 Int_t max = stack->GetNtrack();
106
107 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef)
108 {
109 TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
b8dc484b 110
111 iHitMin = 0;
112 isNewTrack = kTRUE;
dba3c809 113
29fc2c86 114 if (!trackRefs->GetEntries()) continue;
48217459 115
29fc2c86 116 while (isNewTrack) {
068e4c36 117
48217459 118 AliMUONTrack muonTrack;
b8dc484b 119
48217459 120 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit)
121 {
122 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
123
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();
130
131 Int_t track = trackReference->GetTrack();
132
133 if ( track >= max )
134 {
135 AliWarningStream()
136 << "Track ID " << track
137 << " larger than max number of particles " << max << endl;
138 isNewTrack = kFALSE;
139 break;
140 }
141 if (track != trackSave && iHit != 0) {
142 iHitMin = iHit;
143 trackSave = track;
144 break;
145 }
146
147 Float_t bendingSlope = 0;
148 Float_t nonBendingSlope = 0;
149 Float_t inverseBendingMomentum = 0;
150
151 AliMUONTrackParam trackParam;
152
153 // track parameters at hit
154 trackParam.SetBendingCoor(y);
155 trackParam.SetNonBendingCoor(x);
156 trackParam.SetZ(z);
157
158 if (TMath::Abs(pZ) > 0)
159 {
160 bendingSlope = pY/pZ;
161 nonBendingSlope = pX/pZ;
162 }
163 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
164 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
165
166 trackParam.SetBendingSlope(bendingSlope);
167 trackParam.SetNonBendingSlope(nonBendingSlope);
168 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
169
170 AliMUONHitForRec hitForRec;
171
172 hitForRec.SetBendingCoor(y);
173 hitForRec.SetNonBendingCoor(x);
174 hitForRec.SetZ(z);
175 hitForRec.SetBendingReso2(0.0);
176 hitForRec.SetNonBendingReso2(0.0);
177 Int_t detElemId = hitForRec.GetDetElemId();
178 Int_t iChamber;
179 if (detElemId)
180 {
181 iChamber = detElemId / 100 - 1;
182 }
183 else
184 {
185 iChamber = AliMUONConstants::ChamberNumber(z);
186 }
187 hitForRec.SetChamberNumber(iChamber);
188
189 muonTrack.AddTrackParamAtHit(&trackParam,0);
190 muonTrack.AddHitForRecAtHit(&hitForRec);
191 muonTrack.SetTrackID(track);
192
193 trackSave = track;
194 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
b8dc484b 195 }
48217459 196
b8dc484b 197 // track parameters at vertex
48217459 198 TParticle* particle = stack->Particle(muonTrack.GetTrackID());
b8dc484b 199
48217459 200 if (particle)
201 {
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();
208
209 AliMUONTrackParam trackParam;
210
211 trackParam.SetBendingCoor(y);
212 trackParam.SetNonBendingCoor(x);
213 trackParam.SetZ(z);
214
215 Float_t bendingSlope = 0;
216 Float_t nonBendingSlope = 0;
217 Float_t inverseBendingMomentum = 0;
218
219 if (TMath::Abs(pZ) > 0)
220 {
221 bendingSlope = pY/pZ;
222 nonBendingSlope = pX/pZ;
223 }
224
225 Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
226 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
227
228 TParticlePDG* ppdg = particle->GetPDG(1);
229 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
230 inverseBendingMomentum *= charge;
231
232 trackParam.SetBendingSlope(bendingSlope);
233 trackParam.SetNonBendingSlope(nonBendingSlope);
234 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
235
236 muonTrack.SetTrackParamAtVertex(&trackParam);
b8dc484b 237 }
cc056de1 238
48217459 239 trackRefStore->Add(muonTrack);
b8dc484b 240 } // end while isNewTrack
b8dc484b 241 }
b8dc484b 242
48217459 243 AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
244
245 delete trackRefStore;
b8dc484b 246
48217459 247 return rv;
b8dc484b 248}
249
b8dc484b 250//_____________________________________________________________________________
48217459 251Int_t
252AliMUONRecoCheck::NumberOfEvents() const
b8dc484b 253{
48217459 254 /// Return the number of events
255 if ( fDataInterface )
256 {
257 return fDataInterface->NumberOfEvents();
b8dc484b 258 }
48217459 259 return 0;
b8dc484b 260}
261
262//_____________________________________________________________________________
48217459 263AliMUONVTrackStore*
264AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
b8dc484b 265{
48217459 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)
269
b8dc484b 270 Float_t maxGasGap = 1.; // cm
b8dc484b 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;
48217459 281
282 AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b 283 Int_t iHit1;
30d3ea8a 284 Int_t iChamber = 0, detElemId = 0;
b8dc484b 285 Int_t nRec = 0;
286 Int_t nTrackHits = 0;
48217459 287
b8dc484b 288 hitForRec = new AliMUONHitForRec();
289 trackParam = new AliMUONTrackParam();
48217459 290
291 TIter next(trackRefs.CreateIterator());
292 AliMUONTrack* track;
293
294 while ( ( track = static_cast<AliMUONTrack*>(next())) )
295 {
b8dc484b 296 hitForRecAtHit = track->GetHitForRecAtHit();
297 trackParamAtHit = track->GetTrackParamAtHit();
298 trackParamAtVertex = track->GetTrackParamAtVertex();
299 nTrackHits = hitForRecAtHit->GetEntriesFast();
48217459 300 AliMUONTrack trackNew;
b8dc484b 301 iHit1 = 0;
48217459 302 while (iHit1 < nTrackHits)
303 {
b8dc484b 304 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
305 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
306 xRec1 = hitForRec1->GetNonBendingCoor();
307 yRec1 = hitForRec1->GetBendingCoor();
308 zRec1 = hitForRec1->GetZ();
309 xRec = xRec1;
310 yRec = yRec1;
311 zRec = zRec1;
312 bendingSlope1 = trackParam1->GetBendingSlope();
313 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1 314 bendingMomentum1 = 0;
315 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
48217459 316 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b 317 bendingSlope = bendingSlope1;
318 nonBendingSlope = nonBendingSlope1;
319 bendingMomentum = bendingMomentum1;
320 nRec = 1;
48217459 321 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
322 {
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();
333
334 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
335
336 nRec++;
337 xRec += xRec2;
338 yRec += yRec2;
339 zRec += zRec2;
340 bendingSlope += bendingSlope2;
341 nonBendingSlope += nonBendingSlope2;
342 bendingMomentum += bendingMomentum2;
343 iHit1 = iHit2;
344 }
345
b8dc484b 346 } // end iHit2
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;
353
354 hitForRec->SetNonBendingCoor(xRec);
355 hitForRec->SetBendingCoor(yRec);
356 hitForRec->SetZ(zRec);
30d3ea8a 357 detElemId = hitForRec->GetDetElemId();
358 if (detElemId) iChamber = detElemId / 100 - 1;
359 else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b 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);
de876dd1 368 if (TMath::Abs(bendingMomentum) > 0)
48217459 369 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
370
371 trackNew.AddHitForRecAtHit(hitForRec);
372 trackNew.AddTrackParamAtHit(trackParam,0);
b8dc484b 373
374 iHit1++;
375 } // end iHit1
48217459 376
377 trackNew.SetTrackID(track->GetTrackID());
378 trackNew.SetTrackParamAtVertex(trackParamAtVertex);
379 newMuonTrackRef->Add(trackNew);
b8dc484b 380
381 } // end trackRef
b8dc484b 382
b8dc484b 383 delete hitForRec;
384 delete trackParam;
48217459 385 return newMuonTrackRef;
b8dc484b 386}
387
388//_____________________________________________________________________________
48217459 389AliMUONVTrackStore*
390AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
b8dc484b 391{
48217459 392 /// Calculate the number of reconstructible tracks
393
394 AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b 395
b8dc484b 396 TClonesArray* hitForRecAtHit = NULL;
397 AliMUONHitForRec* hitForRec;
398 Float_t zRec;
48217459 399 Int_t nTrackHits;
b8dc484b 400 Int_t isChamberInTrack[10];
401 Int_t iChamber = 0;
402 Bool_t isTrackOK = kTRUE;
48217459 403
404 TIter next(trackRefs.CreateIterator());
405 AliMUONTrack* track;
b8dc484b 406
48217459 407 while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
408 {
b8dc484b 409 hitForRecAtHit = track->GetHitForRecAtHit();
410 nTrackHits = hitForRecAtHit->GetEntriesFast();
411 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
48217459 412
b8dc484b 413 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
414 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
415 zRec = hitForRec->GetZ();
48217459 416 iChamber = hitForRec->GetChamberNumber();
b8dc484b 417 if (iChamber < 0 || iChamber > 10) continue;
418 isChamberInTrack[iChamber] = 1;
cc056de1 419 }
420 // track is reconstructible if the particle is depositing a hit
421 // in the following chamber combinations:
48217459 422
b8dc484b 423 isTrackOK = kTRUE;
cc056de1 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;
48217459 431
432 if (isTrackOK)
433 {
434 reconstructibleStore->Add(*track);
435 }
b8dc484b 436 }
48217459 437
438 return reconstructibleStore;
b8dc484b 439}
440