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