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