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