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