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