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
  32 #include "AliMUONRecoCheck.h"
  33 #include "AliMUONRawClusterV2.h"
  34 #include "AliMUONTrack.h"
  35 #include "AliMUONConstants.h"
  36 #include "AliMUONTrackStoreV1.h"
  37
  38 #include "AliMCEventHandler.h"
  39 #include "AliMCEvent.h"
  40 #include "AliStack.h"
  41 #include "AliTrackReference.h"
  42 #include "AliLog.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDMuonTrack.h"
  45
  46 #include <TFile.h>
  47 #include <TTree.h>
  48 #include <TParticle.h>
  49 #include <TParticlePDG.h>
  50 #include <Riostream.h>
  51
  52 /// \cond CLASSIMP
  53 ClassImp(AliMUONRecoCheck)
  54 /// \endcond
  55
  56 //_____________________________________________________________________________
  57 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *esdFileName, Char_t *pathSim)
  58 : TObject(),
  59 fMCEventHandler(new AliMCEventHandler()),
  60 fESDEvent(new AliESDEvent()),
  61 fESDTree (0x0),
  62 fESDFile (0x0),
  63 fCurrentEvent(0),
  64 fTrackRefStore(0x0),
  65 fRecoTrackRefStore(0x0),
  66 fRecoTrackStore(0x0)
  67 {
  68   /// Normal ctor
  69
  70   // TrackRefs and Particules
  71   fMCEventHandler->SetInputPath(pathSim);
  72   fMCEventHandler->InitIO("");
  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
  89 }
  90
  91 //_____________________________________________________________________________
  92 AliMUONRecoCheck::~AliMUONRecoCheck()
  93 {
  94   /// Destructor
  95   delete fMCEventHandler;
  96   delete fESDEvent;
  97   if (fESDFile) fESDFile->Close();
  98   ResetStores();
  99 }
 100
 101 //_____________________________________________________________________________
 102 void 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 //_____________________________________________________________________________
 111 Int_t AliMUONRecoCheck::NumberOfEvents() const
 112 {
 113   /// Return the number of events
 114   if (fESDTree) return fESDTree->GetEntries();
 115   return 0;
 116 }
 117
 118 //_____________________________________________________________________________
 119 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
 120 {
 121   /// Return a track store containing the reconstructed tracks (converted into
 122   /// MUONTrack objects) for a given event
 123   if (event != fCurrentEvent) {
 124     ResetStores();
 125     fCurrentEvent = event;
 126   }
 127
 128   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
 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;
 134     }
 135     MakeReconstructedTracks();
 136     return fRecoTrackStore;
 137   }
 138 }
 139
 140 //_____________________________________________________________________________
 141 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
 142 {
 143   /// Return a track store containing the track references (converted into
 144   /// MUONTrack objects) for a given event
 145   if (event != fCurrentEvent) {
 146     ResetStores();
 147     fCurrentEvent = event;
 148   }
 149
 150   if (fTrackRefStore != 0x0) return fTrackRefStore;
 151   else {
 152     if (!fMCEventHandler->GetEvent(event)) {
 153       AliError(Form("fails to read MC objects for event %d", event));
 154       return 0x0;
 155     }
 156     MakeTrackRefs();
 157     return fTrackRefStore;
 158   }
 159 }
 160
 161 //_____________________________________________________________________________
 162 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
 163 {
 164   /// Return a track store containing the reconstructible tracks for a given event
 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   }
 176 }
 177
 178 //_____________________________________________________________________________
 179 void 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
 199 //_____________________________________________________________________________
 200 void AliMUONRecoCheck::MakeTrackRefs()
 201 {
 202   /// Make reconstructible tracks
 203   AliMUONVTrackStore *tmpTrackRefStore = new AliMUONTrackStoreV1();
 204
 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();
 209
 210   // loop over simulated tracks
 211   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
 212     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
 213
 214     // skip empty trackRefs
 215     if (nHits < 1) continue;
 216
 217     // get the particle charge for further calculation
 218     TParticlePDG* ppdg = particle->GetPDG();
 219     Int_t charge = (Int_t)(ppdg->Charge()/3.0);
 220
 221     AliMUONTrack track;
 222
 223     // loop over simulated track hits
 224     for (Int_t iHit = 0; iHit < nHits; ++iHit) {
 225       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
 226
 227       // skip trackRefs not in MUON
 228       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
 229
 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();
 237
 238       // check chamberId of current trackReference
 239       Int_t chamberId = AliMUONConstants::ChamberNumber(z);
 240       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
 241
 242       // set hit parameters
 243       AliMUONRawClusterV2 hit(chamberId, 0, 0);
 244       hit.SetXYZ(x,y,z);
 245       hit.SetErrXY(0.,0.);
 246
 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;
 254       }
 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
 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
 308     track.GetTrackParamAtCluster()->Sort();
 309     track.SetTrackID(iTrackRef);
 310     tmpTrackRefStore->Add(track);
 311   }
 312
 313   CleanMuonTrackRef(tmpTrackRefStore);
 314
 315   delete tmpTrackRefStore;
 316 }
 317
 318 //_____________________________________________________________________________
 319 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
 320 {
 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)
 324   fTrackRefStore = new AliMUONTrackStoreV1();
 325
 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;
 329
 330   // create iterator
 331   TIter next(tmpTrackRefStore->CreateIterator());
 332
 333   // loop over tmpTrackRef
 334   AliMUONTrack* track;
 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);
 365
 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++;
 378           iHit1 = iHit2;
 379         }
 380
 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();
 400
 401       // set hit parameters
 402       AliMUONRawClusterV2 hit(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0);
 403       hit.SetXYZ(x,y,z);
 404       hit.SetErrXY(0.,0.);
 405
 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);
 417
 418       iHit1++;
 419     }
 420
 421     newTrack.GetTrackParamAtCluster()->Sort();
 422     newTrack.SetTrackID(track->GetTrackID());
 423     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
 424     fTrackRefStore->Add(newTrack);
 425
 426   }
 427
 428 }
 429
 430 //_____________________________________________________________________________
 431 void AliMUONRecoCheck::MakeReconstructibleTracks()
 432 {
 433   /// Isolate the reconstructible tracks
 434   fRecoTrackRefStore = new AliMUONTrackStoreV1();
 435
 436   // create iterator on trackRef
 437   TIter next(fTrackRefStore->CreateIterator());
 438
 439   // loop over trackRef
 440   AliMUONTrack* track;
 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;
 445
 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;
 451     }
 452
 453     // track is reconstructible if the particle is depositing a hit
 454     // in the following chamber combinations:
 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;
 463
 464     // Add reconstructible tracks to fRecoTrackRefStore
 465     if (trackOK) fRecoTrackRefStore->Add(*track);
 466
 467     delete [] chamberInTrack;
 468   }
 469
 470 }
 471