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