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