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