Preprocessor with new DA
[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
9bdda5f6 26// 13 Nov 2007:
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>
31
cc056de1 32#include "AliMUONRecoCheck.h"
96ebe67e 33#include "AliMUONRawClusterV2.h"
cc056de1 34#include "AliMUONTrack.h"
cc056de1 35#include "AliMUONConstants.h"
48217459 36#include "AliMUONDataInterface.h"
96ebe67e 37#include "AliMUONTrackStoreV1.h"
38
39#include "AliMCEventHandler.h"
40#include "AliMCEvent.h"
b8dc484b 41#include "AliStack.h"
cc056de1 42#include "AliTrackReference.h"
9bdda5f6 43#include "AliLog.h"
44#include "AliESDEvent.h"
45#include "AliESDMuonTrack.h"
96ebe67e 46
cc056de1 47#include <TParticle.h>
5b1afea0 48#include <TParticlePDG.h>
96ebe67e 49#include <Riostream.h>
50
9bdda5f6 51#include <cassert>
52
5398f946 53/// \cond CLASSIMP
b8dc484b 54ClassImp(AliMUONRecoCheck)
5398f946 55/// \endcond
b8dc484b 56
57//_____________________________________________________________________________
96ebe67e 58AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *pathSim)
48217459 59: TObject(),
96ebe67e 60fMCEventHandler(new AliMCEventHandler()),
61fDataInterface(new AliMUONDataInterface(chLoader)),
62fCurrentEvent(0),
63fTrackRefStore(0x0),
64fRecoTrackRefStore(0x0),
65fRecoTrackStore(0x0)
068e4c36 66{
48217459 67 /// Normal ctor
96ebe67e 68 fMCEventHandler->SetInputPath(pathSim);
69 fMCEventHandler->InitIO("");
068e4c36 70}
71
72//_____________________________________________________________________________
48217459 73AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b 74{
48217459 75 /// Destructor
96ebe67e 76 delete fMCEventHandler;
48217459 77 delete fDataInterface;
96ebe67e 78 ResetStores();
79}
80
81//_____________________________________________________________________________
82void AliMUONRecoCheck::ResetStores()
83{
84 /// Deletes all the store objects that have been created and resets the pointers to 0x0
85 delete fTrackRefStore; fTrackRefStore = 0x0;
86 delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
87 delete fRecoTrackStore; fRecoTrackStore = 0x0;
88}
89
90//_____________________________________________________________________________
91Int_t AliMUONRecoCheck::NumberOfEvents() const
92{
93 /// Return the number of events
94 if (fDataInterface->IsValid()) return fDataInterface->NumberOfEvents();
95 return 0;
b8dc484b 96}
97
48217459 98//_____________________________________________________________________________
96ebe67e 99AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
48217459 100{
101 /// Return the reconstructed track store for a given event
96ebe67e 102 if (!fDataInterface->IsValid()) return new AliMUONTrackStoreV1();
103
104 if (event != fCurrentEvent) {
105 ResetStores();
106 fCurrentEvent = event;
107 }
108
109 if (fRecoTrackStore != 0x0) return fRecoTrackStore;
110
111 fRecoTrackStore = new AliMUONTrackStoreV1();
112
113 return fRecoTrackStore;
48217459 114}
45a05f07 115
b8dc484b 116//_____________________________________________________________________________
9bdda5f6 117AliMUONVTrackStore*
118AliMUONRecoCheck::ReconstructedTracks(AliESDEvent* esd, Bool_t padMissing)
119{
120 /// Create and return a list of reconstructed tracks from ESD data.
121 /// Note: the returned object must be deleted by the caller with the delete operator!
122 /// For example:
123 ///
124 /// AliMUONVTrackStore* tracks = AliMUONRecoCheck::ReconstructedTracks(esd);
125 /// delete tracks;
126 ///
127 /// If the padMissing flag is set to kTRUE then the missing hits on chambers [2..14]
128 /// are added with the coordinate (0, 0, 0). This should be fixed in the MUON ESD
129 /// structure, hopefully soon.
130
131 assert( esd != NULL );
132
133 AliMUONVTrackStore* store = new AliMUONTrackStoreV1;
134
135 Int_t nTracks = Int_t(esd->GetNumberOfMuonTracks());
136 for (Int_t i = 0; i < nTracks; i++)
137 {
138 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(i);
139 AliMUONTrack track;
140
141 // Fill the track parameters at the vertex.
142 AliMUONTrackParam vertexParams, ch1Params;
143 vertexParams.GetParamFrom(*esdTrack);
144 ch1Params.GetParamFromUncorrected(*esdTrack);
145 ch1Params.GetCovFrom(*esdTrack);
146
147 // Create the hit on chamber 1 (assume the same values as uncorrected track params)
148 AliMUONRawClusterV2 hit(0, 0, 0); // ChamberID = 0, but unknown DE or cluster index so set those to zero.
149 hit.SetXYZ(ch1Params.GetNonBendingCoor(), ch1Params.GetBendingCoor(), ch1Params.GetZ());
150 hit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());
151
152 track.SetTrackParamAtVertex(&vertexParams);
153 track.AddTrackParamAtCluster(ch1Params, hit, kTRUE);
154
155 if (padMissing)
156 {
157 // We need to add fake hits so that AliMUONTrack::CompatibleTrack will be able to
158 // match these tracks to the reference tracks returned by the TrackRefs method.
159 for (Int_t chamber = 1; chamber < AliMUONConstants::NTrackingCh(); chamber++)
160 {
161 if (esdTrack->IsInMuonClusterMap(chamber))
162 {
163 AliMUONTrackParam fakeParams;
164 fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
165 fakeParams.SetBendingCoor(0);
166 fakeParams.SetNonBendingCoor(0);
167 AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
168 fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
169 fakeHit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());
170 track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
171 }
172 }
173 for (Int_t chamber = AliMUONConstants::NTrackingCh(); chamber < AliMUONConstants::NCh(); chamber++)
174 {
175 Int_t bit = 3 - (chamber - AliMUONConstants::NTrackingCh());
176 if ( ((esdTrack->GetHitsPatternInTrigCh() >> bit) & 0x11) != 0 )
177 {
178 AliMUONTrackParam fakeParams;
179 fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
180 fakeParams.SetBendingCoor(0);
181 fakeParams.SetNonBendingCoor(0);
182 AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
183 fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
184 fakeHit.SetErrXY(1, 1);
185 track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
186 }
187 }
188 }
189 store->Add(track);
190 }
191
192 return store;
193}
194
195//_____________________________________________________________________________
96ebe67e 196AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b 197{
48217459 198 /// Return a track store containing the track references (converted into
199 /// MUONTrack objects) for a given event
96ebe67e 200 if (event != fCurrentEvent) {
201 ResetStores();
202 fCurrentEvent = event;
203 }
204
205 if (fTrackRefStore != 0x0) return fTrackRefStore;
206 else {
207 if (!fMCEventHandler->GetEvent(event)) return 0x0;
208 MakeTrackRefs();
209 return fTrackRefStore;
210 }
b8dc484b 211}
212
213//_____________________________________________________________________________
96ebe67e 214AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b 215{
48217459 216 /// Return a track store containing the reconstructible tracks for a given event
96ebe67e 217 if (event != fCurrentEvent) {
218 ResetStores();
219 fCurrentEvent = event;
220 }
221
222 if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
223 else {
224 if (TrackRefs(event) == 0x0) return 0x0;
225 MakeReconstructibleTracks();
226 return fRecoTrackRefStore;
227 }
48217459 228}
b8dc484b 229
48217459 230//_____________________________________________________________________________
96ebe67e 231void AliMUONRecoCheck::MakeTrackRefs()
48217459 232{
233 /// Make reconstructible tracks
96ebe67e 234 AliMUONVTrackStore *tmpTrackRefStore = new AliMUONTrackStoreV1();
48217459 235
96ebe67e 236 Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
237 TParticle* particle;
238 TClonesArray* trackRefs;
239 Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
b8dc484b 240
96ebe67e 241 // loop over simulated tracks
242 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
243 Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
48217459 244
96ebe67e 245 // skip empty trackRefs
246 if (nHits < 1) continue;
247
248 // skip trackRefs not in MUON
249 AliTrackReference* trackReference0 = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(0));
250 if (trackReference0->DetectorId() != AliTrackReference::kMUON) continue;
251
252 AliMUONTrack track;
b8dc484b 253
96ebe67e 254 // get track parameters at particle's vertex
255 x = particle->Vx();
256 y = particle->Vy();
257 z = particle->Vz();
258 pX = particle->Px();
259 pY = particle->Py();
260 pZ = particle->Pz();
dba3c809 261
96ebe67e 262 // compute rest of track parameters at particle's vertex
263 bendingSlope = 0;
264 nonBendingSlope = 0;
265 inverseBendingMomentum = 0;
266 if (TMath::Abs(pZ) > 0) {
267 bendingSlope = pY/pZ;
268 nonBendingSlope = pX/pZ;
269 }
270 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
271 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
272 TParticlePDG* ppdg = particle->GetPDG(1);
273 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
274 inverseBendingMomentum *= charge;
275
276 // set track parameters at particle's vertex
277 AliMUONTrackParam trackParamAtVertex;
278 trackParamAtVertex.SetNonBendingCoor(x);
279 trackParamAtVertex.SetBendingCoor(y);
280 trackParamAtVertex.SetZ(z);
281 trackParamAtVertex.SetBendingSlope(bendingSlope);
282 trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
283 trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
284
285 // add track parameters at vertex
286 track.SetTrackParamAtVertex(&trackParamAtVertex);
48217459 287
96ebe67e 288 // loop over simulated track hits
289 for (Int_t iHit = 0; iHit < nHits; ++iHit) {
290 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
068e4c36 291
96ebe67e 292 // Get track parameters of current hit
293 x = trackReference->X();
294 y = trackReference->Y();
295 z = trackReference->Z();
296 pX = trackReference->Px();
297 pY = trackReference->Py();
298 pZ = trackReference->Pz();
b8dc484b 299
96ebe67e 300 // check chamberId of current trackReference
301 Int_t chamberId = AliMUONConstants::ChamberNumber(z);
302 if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
48217459 303
96ebe67e 304 // set hit parameters
305 AliMUONRawClusterV2 hit(chamberId, 0, 0);
306 hit.SetXYZ(x,y,z);
307 hit.SetErrXY(0.,0.);
b8dc484b 308
96ebe67e 309 // compute track parameters at hit
310 Double_t bendingSlope = 0;
311 Double_t nonBendingSlope = 0;
312 Double_t inverseBendingMomentum = 0;
313 if (TMath::Abs(pZ) > 0) {
314 bendingSlope = pY/pZ;
315 nonBendingSlope = pX/pZ;
b8dc484b 316 }
96ebe67e 317 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
318 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
319 inverseBendingMomentum *= charge;
320
321 // set track parameters at hit
322 AliMUONTrackParam trackParam;
323 trackParam.SetNonBendingCoor(x);
324 trackParam.SetBendingCoor(y);
325 trackParam.SetZ(z);
326 trackParam.SetBendingSlope(bendingSlope);
327 trackParam.SetNonBendingSlope(nonBendingSlope);
328 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
329
330 // add track parameters at current hit to the track
331 track.AddTrackParamAtCluster(trackParam,hit,kTRUE);
332 }
333
334 track.GetTrackParamAtCluster()->Sort();
335 track.SetTrackID(iTrackRef);
336 tmpTrackRefStore->Add(track);
b8dc484b 337 }
b8dc484b 338
96ebe67e 339 CleanMuonTrackRef(tmpTrackRefStore);
48217459 340
96ebe67e 341 delete tmpTrackRefStore;
b8dc484b 342}
343
344//_____________________________________________________________________________
96ebe67e 345void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
b8dc484b 346{
48217459 347 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
348 /// each chamber (one when particle is entering + one when particle is leaving
349 /// the sensitive volume)
96ebe67e 350 fTrackRefStore = new AliMUONTrackStoreV1();
48217459 351
96ebe67e 352 Double_t maxGasGap = 1.; // cm
353 Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
354 Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
48217459 355
96ebe67e 356 // create iterator
357 TIter next(tmpTrackRefStore->CreateIterator());
48217459 358
96ebe67e 359 // loop over tmpTrackRef
48217459 360 AliMUONTrack* track;
96ebe67e 361 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
362
363 AliMUONTrack newTrack;
364
365 // loop over tmpTrackRef's hits
366 Int_t iHit1 = 0;
367 Int_t nTrackHits = track->GetNClusters();
368 while (iHit1 < nTrackHits) {
369 AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
370
371 // get track parameters at hit1
372 x1 = trackParam1->GetNonBendingCoor();
373 y1 = trackParam1->GetBendingCoor();
374 z1 = trackParam1->GetZ();
375 pX1 = trackParam1->Px();
376 pY1 = trackParam1->Py();
377 pZ1 = trackParam1->Pz();
378
379 // prepare new track parameters
380 x = x1;
381 y = y1;
382 z = z1;
383 pX = pX1;
384 pY = pY1;
385 pZ = pZ1;
386
387 // loop over next tmpTrackRef's hits
388 Int_t nCombinedHits = 1;
389 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
390 AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
48217459 391
96ebe67e 392 // get z position of hit2
393 z2 = trackParam2->GetZ();
394
395 // complete new track parameters if hit2 is on the same detection element
396 if ( TMath::Abs(z2-z1) < maxGasGap ) {
397 x += trackParam2->GetNonBendingCoor();
398 y += trackParam2->GetBendingCoor();
399 z += z2;
400 pX += trackParam2->Px();
401 pY += trackParam2->Py();
402 pZ += trackParam2->Pz();
403 nCombinedHits++;
48217459 404 iHit1 = iHit2;
405 }
406
96ebe67e 407 }
408
409 // finalize new track parameters
410 x /= (Double_t)nCombinedHits;
411 y /= (Double_t)nCombinedHits;
412 z /= (Double_t)nCombinedHits;
413 pX /= (Double_t)nCombinedHits;
414 pY /= (Double_t)nCombinedHits;
415 pZ /= (Double_t)nCombinedHits;
416 bendingSlope = 0;
417 nonBendingSlope = 0;
418 inverseBendingMomentum = 0;
419 if (TMath::Abs(pZ) > 0) {
420 bendingSlope = pY/pZ;
421 nonBendingSlope = pX/pZ;
422 }
423 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
424 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
425 inverseBendingMomentum *= trackParam1->GetCharge();
b8dc484b 426
96ebe67e 427 // set hit parameters
428 AliMUONRawClusterV2 hit(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0);
429 hit.SetXYZ(x,y,z);
430 hit.SetErrXY(0.,0.);
48217459 431
96ebe67e 432 // set new track parameters at new hit
433 AliMUONTrackParam trackParam;
434 trackParam.SetNonBendingCoor(x);
435 trackParam.SetBendingCoor(y);
436 trackParam.SetZ(z);
437 trackParam.SetBendingSlope(bendingSlope);
438 trackParam.SetNonBendingSlope(nonBendingSlope);
439 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
440
441 // add track parameters at current hit to the track
442 newTrack.AddTrackParamAtCluster(trackParam,hit,kTRUE);
b8dc484b 443
444 iHit1++;
96ebe67e 445 }
48217459 446
96ebe67e 447 newTrack.GetTrackParamAtCluster()->Sort();
448 newTrack.SetTrackID(track->GetTrackID());
449 newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
450 fTrackRefStore->Add(newTrack);
b8dc484b 451
96ebe67e 452 }
b8dc484b 453
b8dc484b 454}
455
456//_____________________________________________________________________________
96ebe67e 457void AliMUONRecoCheck::MakeReconstructibleTracks()
b8dc484b 458{
48217459 459 /// Calculate the number of reconstructible tracks
96ebe67e 460 fRecoTrackRefStore = new AliMUONTrackStoreV1();
48217459 461
96ebe67e 462 // create iterator on trackRef
463 TIter next(fTrackRefStore->CreateIterator());
b8dc484b 464
96ebe67e 465 // loop over trackRef
48217459 466 AliMUONTrack* track;
96ebe67e 467 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
468
469 Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
470 for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
48217459 471
96ebe67e 472 // loop over trackRef's hits to get hit chambers
473 Int_t nTrackHits = track->GetNClusters();
474 for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
475 AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
476 chamberInTrack[hit->GetChamberId()] = kTRUE;
cc056de1 477 }
96ebe67e 478
cc056de1 479 // track is reconstructible if the particle is depositing a hit
480 // in the following chamber combinations:
96ebe67e 481 Bool_t trackOK = kTRUE;
482 if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
483 if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
484 if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
485 Int_t nHitsInLastStations = 0;
486 for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
487 if (chamberInTrack[iCh]) nHitsInLastStations++;
488 if(nHitsInLastStations < 3) trackOK = kFALSE;
48217459 489
96ebe67e 490 // Add reconstructible tracks to fRecoTrackRefStore
491 if (trackOK) fRecoTrackRefStore->Add(*track);
48217459 492
96ebe67e 493 delete [] chamberInTrack;
b8dc484b 494 }
48217459 495
b8dc484b 496}
497