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