MUON/AliMUONRecoCheck.cxx

   1 /**************************************************************************
   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 **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 /// \class AliMUONRecoCheck
  20 /// Utility class to check reconstruction
  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.
  24 //-----------------------------------------------------------------------------
  25
  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 // 25 Jan 2008:
  32 // Use the new ESDInterface to create MUON objects from ESD data
  33 // - Philippe Pillot
  34 //
  35
  36 #include "AliMUONRecoCheck.h"
  37 #include "AliMUONTrack.h"
  38 #include "AliMUONVTrackStore.h"
  39 #include "AliMUONVCluster.h"
  40 #include "AliMUONVClusterStore.h"
  41 #include "AliMUONConstants.h"
  42 #include "AliMUONESDInterface.h"
  43 #include "AliMUONTrackParam.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliTrackReference.h"
  48 #include "AliLog.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESDMuonTrack.h"
  51
  52 #include <TFile.h>
  53 #include <TTree.h>
  54 #include <TParticle.h>
  55 #include <TParticlePDG.h>
  56 #include <Riostream.h>
  57
  58 /// \cond CLASSIMP
  59 ClassImp(AliMUONRecoCheck)
  60 /// \endcond
  61
  62 //_____________________________________________________________________________
  63 AliMUONRecoCheck::AliMUONRecoCheck(const Char_t *esdFileName, const Char_t *pathSim)
  64 : TObject(),
  65 fMCEventHandler(new AliMCEventHandler()),
  66 fESDEvent(new AliESDEvent()),
  67 fESDTree (0x0),
  68 fESDFile (0x0),
  69 fCurrentEvent(0),
  70 fTrackRefStore(0x0),
  71 fRecoTrackRefStore(0x0),
  72 fRecoTrackStore(0x0),
  73 fESDEventOwner(kTRUE)
  74 {
  75   /// Normal ctor
  76
  77   // TrackRefs and Particules
  78   fMCEventHandler->SetInputPath(pathSim);
  79   fMCEventHandler->InitIO("");
  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
  96 }
  97
  98 //_____________________________________________________________________________
  99 AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent *esdEvent, AliMCEventHandler *mcEventHandler)
 100 : TObject(),
 101 fMCEventHandler(0),
 102 fESDEvent(0),
 103 fESDTree (0x0),
 104 fESDFile (0x0),
 105 fCurrentEvent(0),
 106 fTrackRefStore(0x0),
 107 fRecoTrackRefStore(0x0),
 108 fRecoTrackStore(0x0),
 109 fESDEventOwner(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 //_____________________________________________________________________________
 122 AliMUONRecoCheck::~AliMUONRecoCheck()
 123 {
 124   /// Destructor
 125   if (fESDEventOwner) {
 126     delete fMCEventHandler;
 127     delete fESDEvent;
 128     if (fESDFile) fESDFile->Close();
 129   }
 130   ResetStores();
 131 }
 132
 133 //_____________________________________________________________________________
 134 void 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 //_____________________________________________________________________________
 143 Int_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 //_____________________________________________________________________________
 158 Int_t AliMUONRecoCheck::NumberOfEvents() const
 159 {
 160   /// Return the number of events
 161   if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
 162   return 0;
 163 }
 164
 165 //_____________________________________________________________________________
 166 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event, Bool_t refit)
 167 {
 168   /// Return a track store containing the reconstructed tracks (converted into
 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.
 172
 173   if (!fESDEventOwner) {
 174     if (fRecoTrackStore == 0x0) MakeReconstructedTracks(refit);
 175     return fRecoTrackStore;
 176   }
 177
 178   if (event != fCurrentEvent) {
 179     ResetStores();
 180     fCurrentEvent = event;
 181   }
 182
 183   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
 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;
 189     }
 190     MakeReconstructedTracks(refit);
 191     return fRecoTrackStore;
 192   }
 193 }
 194
 195 //_____________________________________________________________________________
 196 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
 197 {
 198   /// Return a track store containing the track references (converted into
 199   /// MUONTrack objects) for a given event
 200
 201   if (!fESDEventOwner) {
 202     if (fTrackRefStore == 0x0) MakeTrackRefs();
 203     return fTrackRefStore;
 204   }
 205
 206   if (event != fCurrentEvent) {
 207     ResetStores();
 208     fCurrentEvent = event;
 209   }
 210
 211   if (fTrackRefStore != 0x0) return fTrackRefStore;
 212   else {
 213     if (!fMCEventHandler->GetEvent(event)) {
 214       AliError(Form("fails to read MC objects for event %d", event));
 215       return 0x0;
 216     }
 217     MakeTrackRefs();
 218     return fTrackRefStore;
 219   }
 220 }
 221
 222 //_____________________________________________________________________________
 223 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event, UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
 224 {
 225   /// Return a track store containing the reconstructible tracks for a given event,
 226   /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5.
 227
 228   if (!fESDEventOwner) {
 229     if (fRecoTrackRefStore == 0x0) {
 230       if (TrackRefs(event) == 0x0) return 0x0;
 231       MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
 232     }
 233     return fRecoTrackRefStore;
 234   }
 235
 236   if (event != fCurrentEvent) {
 237     ResetStores();
 238     fCurrentEvent = event;
 239   }
 240
 241   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
 242   else {
 243     if (TrackRefs(event) == 0x0) return 0x0;
 244     MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
 245     return fRecoTrackRefStore;
 246   }
 247 }
 248
 249 //_____________________________________________________________________________
 250 void AliMUONRecoCheck::MakeReconstructedTracks(Bool_t refit)
 251 {
 252   /// Make reconstructed tracks
 253   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
 254
 255   // loop over all reconstructed tracks and add them to the store (skip ghosts)
 256   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
 257   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 258     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
 259     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore, refit);
 260   }
 261
 262 }
 263
 264 //_____________________________________________________________________________
 265 void AliMUONRecoCheck::MakeTrackRefs()
 266 {
 267   /// Make reconstructible tracks
 268   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
 269   if (!tmpTrackRefStore) return;
 270
 271   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
 272   TParticle* particle;
 273   TClonesArray* trackRefs;
 274   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
 275   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
 276   if (!cStore) return;
 277   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
 278
 279   // loop over simulated tracks
 280   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
 281     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
 282
 283     // skip empty trackRefs
 284     if (nHits < 1) continue;
 285
 286     // get the particle charge for further calculation
 287     TParticlePDG* ppdg = particle->GetPDG();
 288     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
 289
 290     AliMUONTrack track;
 291
 292     // loop over simulated track hits
 293     for (Int_t iHit = 0; iHit < nHits; ++iHit) {
 294       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
 295
 296       // skip trackRefs not in MUON
 297       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
 298
 299       // Get track parameters of current hit
 300       x = trackReference->X();
 301       y = trackReference->Y();
 302       z = trackReference->Z();
 303       pX = trackReference->Px();
 304       pY = trackReference->Py();
 305       pZ = trackReference->Pz();
 306
 307       // check chamberId of current trackReference
 308       Int_t detElemId = trackReference->UserId();
 309       Int_t chamberId = detElemId / 100 - 1;
 310       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
 311
 312       // set hit parameters
 313       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
 314       hit->SetXYZ(x,y,z);
 315       hit->SetErrXY(0.,0.);
 316
 317       // compute track parameters at hit
 318       bendingSlope = 0;
 319       nonBendingSlope = 0;
 320       inverseBendingMomentum = 0;
 321       if (TMath::Abs(pZ) > 0) {
 322         bendingSlope = pY/pZ;
 323         nonBendingSlope = pX/pZ;
 324       }
 325       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 326       if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 327       inverseBendingMomentum *= charge;
 328
 329       // set track parameters at hit
 330       AliMUONTrackParam trackParam;
 331       trackParam.SetNonBendingCoor(x);
 332       trackParam.SetBendingCoor(y);
 333       trackParam.SetZ(z);
 334       trackParam.SetBendingSlope(bendingSlope);
 335       trackParam.SetNonBendingSlope(nonBendingSlope);
 336       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
 337
 338       // add track parameters at current hit to the track
 339       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
 340     }
 341
 342     // if none of the track hits was in MUON, goto the next track
 343     if (track.GetNClusters() < 1) continue;
 344
 345     // get track parameters at particle's vertex
 346     x = particle->Vx();
 347     y = particle->Vy();
 348     z = particle->Vz();
 349     pX = particle->Px();
 350     pY = particle->Py();
 351     pZ = particle->Pz();
 352
 353     // compute rest of track parameters at particle's vertex
 354     bendingSlope = 0;
 355     nonBendingSlope = 0;
 356     inverseBendingMomentum = 0;
 357     if (TMath::Abs(pZ) > 0) {
 358       bendingSlope = pY/pZ;
 359       nonBendingSlope = pX/pZ;
 360     }
 361     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 362     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 363     inverseBendingMomentum *= charge;
 364
 365     // set track parameters at particle's vertex
 366     AliMUONTrackParam trackParamAtVertex;
 367     trackParamAtVertex.SetNonBendingCoor(x);
 368     trackParamAtVertex.SetBendingCoor(y);
 369     trackParamAtVertex.SetZ(z);
 370     trackParamAtVertex.SetBendingSlope(bendingSlope);
 371     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
 372     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
 373
 374     // add track parameters at vertex
 375     track.SetTrackParamAtVertex(&trackParamAtVertex);
 376
 377     // store the track
 378     track.SetUniqueID(iTrackRef);
 379     tmpTrackRefStore->Add(track);
 380   }
 381
 382   CleanMuonTrackRef(tmpTrackRefStore);
 383
 384   delete hit;
 385   delete cStore;
 386   delete tmpTrackRefStore;
 387 }
 388
 389 //_____________________________________________________________________________
 390 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
 391 {
 392   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
 393   /// each chamber (one when particle is entering + one when particle is leaving
 394   /// the sensitive volume)
 395   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
 396
 397   Double_t maxGasGap = 1.; // cm
 398   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
 399   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
 400   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
 401   if (!cStore) return;
 402   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
 403
 404   // create iterator
 405   TIter next(tmpTrackRefStore->CreateIterator());
 406
 407   // loop over tmpTrackRef
 408   AliMUONTrack* track;
 409   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
 410
 411     AliMUONTrack newTrack;
 412
 413     // loop over tmpTrackRef's hits
 414     Int_t iHit1 = 0;
 415     Int_t nTrackHits = track->GetNClusters();
 416     while (iHit1 < nTrackHits) {
 417       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
 418
 419       // get track parameters at hit1
 420       x1  = trackParam1->GetNonBendingCoor();
 421       y1  = trackParam1->GetBendingCoor();
 422       z1  = trackParam1->GetZ();
 423       pX1 = trackParam1->Px();
 424       pY1 = trackParam1->Py();
 425       pZ1 = trackParam1->Pz();
 426
 427       // prepare new track parameters
 428       x  = x1;
 429       y  = y1;
 430       z  = z1;
 431       pX = pX1;
 432       pY = pY1;
 433       pZ = pZ1;
 434
 435       // loop over next tmpTrackRef's hits
 436       Int_t nCombinedHits = 1;
 437       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
 438         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
 439
 440         // get z position of hit2
 441         z2 = trackParam2->GetZ();
 442
 443         // complete new track parameters if hit2 is on the same detection element
 444         if ( TMath::Abs(z2-z1) < maxGasGap ) {
 445           x  += trackParam2->GetNonBendingCoor();
 446           y  += trackParam2->GetBendingCoor();
 447           z  += z2;
 448           pX += trackParam2->Px();
 449           pY += trackParam2->Py();
 450           pZ += trackParam2->Pz();
 451           nCombinedHits++;
 452           iHit1 = iHit2;
 453         }
 454
 455       }
 456
 457       // finalize new track parameters
 458       x  /= (Double_t)nCombinedHits;
 459       y  /= (Double_t)nCombinedHits;
 460       z  /= (Double_t)nCombinedHits;
 461       pX /= (Double_t)nCombinedHits;
 462       pY /= (Double_t)nCombinedHits;
 463       pZ /= (Double_t)nCombinedHits;
 464       bendingSlope = 0;
 465       nonBendingSlope = 0;
 466       inverseBendingMomentum = 0;
 467       if (TMath::Abs(pZ) > 0) {
 468         bendingSlope = pY/pZ;
 469         nonBendingSlope = pX/pZ;
 470       }
 471       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 472       if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 473       inverseBendingMomentum *= trackParam1->GetCharge();
 474
 475       // set hit parameters
 476       hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
 477       hit->SetXYZ(x,y,z);
 478       hit->SetErrXY(0.,0.);
 479
 480       // set new track parameters at new hit
 481       AliMUONTrackParam trackParam;
 482       trackParam.SetNonBendingCoor(x);
 483       trackParam.SetBendingCoor(y);
 484       trackParam.SetZ(z);
 485       trackParam.SetBendingSlope(bendingSlope);
 486       trackParam.SetNonBendingSlope(nonBendingSlope);
 487       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
 488
 489       // add track parameters at current hit to the track
 490       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
 491
 492       iHit1++;
 493     }
 494
 495     newTrack.SetUniqueID(track->GetUniqueID());
 496     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
 497     fTrackRefStore->Add(newTrack);
 498
 499   }
 500
 501   delete hit;
 502   delete cStore;
 503 }
 504
 505 //_____________________________________________________________________________
 506 void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
 507 {
 508   /// Isolate the reconstructible tracks
 509   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
 510
 511   // create iterator on trackRef
 512   TIter next(fTrackRefStore->CreateIterator());
 513
 514   // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore
 515   AliMUONTrack* track;
 516   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
 517     if (track->IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(*track);
 518   }
 519
 520 }
 521
 522 //_____________________________________________________________________________
 523 AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore,
 524                                                     Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut)
 525 {
 526   /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters.
 527   /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel".
 528   /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track
 529
 530   AliMUONTrack *matchedTrack = 0x0;
 531   nMatchClusters = 0;
 532
 533   if (useLabel) { // by using the MC label
 534
 535     // get the corresponding simulated track if any
 536     Int_t label = track.GetMCLabel();
 537     matchedTrack = (AliMUONTrack*) trackStore.FindObject(label);
 538
 539     // get the fraction of matched clusters
 540     if (matchedTrack) {
 541       Int_t nClusters = track.GetNClusters();
 542       for (Int_t iCl = 0; iCl < nClusters; iCl++)
 543         if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
 544           nMatchClusters++;
 545     }
 546
 547   } else { // by comparing cluster/TrackRef positions
 548
 549     // look for the corresponding simulated track if any
 550     TIter next(trackStore.CreateIterator());
 551     AliMUONTrack* track2;
 552     while ( ( track2 = static_cast<AliMUONTrack*>(next()) ) ) {
 553
 554       // check compatibility
 555       Int_t n = 0;
 556       if (track.Match(*track2, sigmaCut, n)) {
 557         matchedTrack = track2;
 558         nMatchClusters = n;
 559         break;
 560       }
 561
 562     }
 563
 564   }
 565
 566   return matchedTrack;
 567
 568 }
 569