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