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 "AliMUONHitForRec.h"
  28 #include "AliMUONTrack.h"
  29 #include "AliMUONConstants.h"
  30 #include "AliMUONMCDataInterface.h"
  31 #include "AliMUONDataInterface.h"
  32 #include "AliStack.h"
  33 #include "AliTrackReference.h"
  34 #include "AliLog.h"
  35 #include "AliMUONTrackStoreV1.h"
  36 #include <Riostream.h>
  37 #include <TParticle.h>
  38 #include <TParticlePDG.h>
  39
  40 /// \cond CLASSIMP
  41 ClassImp(AliMUONRecoCheck)
  42 /// \endcond
  43
  44 //_____________________________________________________________________________
  45 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
  46 : TObject(),
  47 fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
  48 fDataInterface(new AliMUONDataInterface(chLoader))
  49 {
  50   /// Normal ctor
  51 }
  52
  53 //_____________________________________________________________________________
  54 AliMUONRecoCheck::~AliMUONRecoCheck()
  55 {
  56   /// Destructor
  57   delete fMCDataInterface;
  58   delete fDataInterface;
  59 }
  60
  61 //_____________________________________________________________________________
  62 AliMUONVTrackStore*
  63 AliMUONRecoCheck::ReconstructedTracks(Int_t event)
  64 {
  65   /// Return the reconstructed track store for a given event
  66   return fDataInterface->TrackStore(event);
  67 }
  68
  69 //_____________________________________________________________________________
  70 AliMUONVTrackStore*
  71 AliMUONRecoCheck::TrackRefs(Int_t event)
  72 {
  73   /// Return a track store containing the track references (converted into
  74   /// MUONTrack objects) for a given event
  75   return MakeTrackRefs(event);
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliMUONVTrackStore*
  80 AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
  81 {
  82   /// Return a track store containing the reconstructible tracks for a given event
  83   AliMUONVTrackStore* tmp = MakeTrackRefs(event);
  84   AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
  85   delete tmp;
  86   return reconstructible;
  87 }
  88
  89 //_____________________________________________________________________________
  90 AliMUONVTrackStore*
  91 AliMUONRecoCheck::MakeTrackRefs(Int_t event)
  92 {
  93   /// Make reconstructible tracks
  94
  95   AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
  96
  97   Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
  98
  99   Int_t trackSave(-999);
 100   Int_t iHitMin(0);
 101
 102   Bool_t isNewTrack;
 103
 104   AliStack* stack = fMCDataInterface->Stack(event);
 105   Int_t max = stack->GetNtrack();
 106
 107   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef)
 108   {
 109     TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
 110
 111     iHitMin = 0;
 112     isNewTrack = kTRUE;
 113
 114     if (!trackRefs->GetEntries()) continue;
 115
 116     while (isNewTrack) {
 117
 118       AliMUONTrack muonTrack;
 119
 120       for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit)
 121       {
 122         AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
 123
 124         Float_t x = trackReference->X();
 125         Float_t y = trackReference->Y();
 126         Float_t z = trackReference->Z();
 127         Float_t pX = trackReference->Px();
 128         Float_t pY = trackReference->Py();
 129         Float_t pZ = trackReference->Pz();
 130
 131         Int_t track = trackReference->GetTrack();
 132
 133         if ( track >= max )
 134         {
 135           AliWarningStream()
 136           << "Track ID " << track
 137           << " larger than max number of particles " << max << endl;
 138           isNewTrack = kFALSE;
 139           break;
 140         }
 141         if (track != trackSave && iHit != 0) {
 142           iHitMin = iHit;
 143           trackSave = track;
 144           break;
 145         }
 146
 147         Float_t bendingSlope = 0;
 148         Float_t nonBendingSlope = 0;
 149         Float_t inverseBendingMomentum = 0;
 150
 151         AliMUONTrackParam trackParam;
 152
 153         // track parameters at hit
 154         trackParam.SetBendingCoor(y);
 155         trackParam.SetNonBendingCoor(x);
 156         trackParam.SetZ(z);
 157
 158         if (TMath::Abs(pZ) > 0)
 159         {
 160           bendingSlope = pY/pZ;
 161           nonBendingSlope = pX/pZ;
 162         }
 163         Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 164         if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 165
 166         trackParam.SetBendingSlope(bendingSlope);
 167         trackParam.SetNonBendingSlope(nonBendingSlope);
 168         trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
 169
 170         AliMUONHitForRec hitForRec;
 171
 172         hitForRec.SetBendingCoor(y);
 173         hitForRec.SetNonBendingCoor(x);
 174         hitForRec.SetZ(z);
 175         hitForRec.SetBendingReso2(0.0);
 176         hitForRec.SetNonBendingReso2(0.0);
 177         Int_t detElemId = hitForRec.GetDetElemId();
 178         Int_t iChamber;
 179         if (detElemId)
 180         {
 181           iChamber = detElemId / 100 - 1;
 182         }
 183         else
 184         {
 185           iChamber = AliMUONConstants::ChamberNumber(z);
 186         }
 187         hitForRec.SetChamberNumber(iChamber);
 188
 189         muonTrack.AddTrackParamAtHit(&trackParam,0);
 190         muonTrack.AddHitForRecAtHit(&hitForRec);
 191         muonTrack.SetTrackID(track);
 192
 193         trackSave = track;
 194         if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
 195       }
 196
 197       // track parameters at vertex
 198       TParticle* particle = stack->Particle(muonTrack.GetTrackID());
 199
 200       if (particle)
 201       {
 202         Float_t x = particle->Vx();
 203         Float_t y = particle->Vy();
 204         Float_t z = particle->Vz();
 205         Float_t pX = particle->Px();
 206         Float_t pY = particle->Py();
 207         Float_t pZ = particle->Pz();
 208
 209         AliMUONTrackParam trackParam;
 210
 211         trackParam.SetBendingCoor(y);
 212         trackParam.SetNonBendingCoor(x);
 213         trackParam.SetZ(z);
 214
 215         Float_t bendingSlope = 0;
 216         Float_t nonBendingSlope = 0;
 217         Float_t inverseBendingMomentum = 0;
 218
 219         if (TMath::Abs(pZ) > 0)
 220         {
 221           bendingSlope = pY/pZ;
 222           nonBendingSlope = pX/pZ;
 223         }
 224
 225         Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 226         if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 227
 228         TParticlePDG* ppdg = particle->GetPDG(1);
 229         Int_t charge = (Int_t)(ppdg->Charge()/3.0);
 230         inverseBendingMomentum *= charge;
 231
 232         trackParam.SetBendingSlope(bendingSlope);
 233         trackParam.SetNonBendingSlope(nonBendingSlope);
 234         trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
 235
 236         muonTrack.SetTrackParamAtVertex(&trackParam);
 237       }
 238
 239       trackRefStore->Add(muonTrack);
 240     } // end while isNewTrack
 241     delete trackRefs;
 242   }
 243
 244   AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
 245
 246   delete trackRefStore;
 247
 248   return rv;
 249 }
 250
 251 //_____________________________________________________________________________
 252 Int_t
 253 AliMUONRecoCheck::NumberOfEvents() const
 254 {
 255   /// Return the number of events
 256   if ( fDataInterface )
 257   {
 258     return fDataInterface->NumberOfEvents();
 259   }
 260   return 0;
 261 }
 262
 263 //_____________________________________________________________________________
 264 AliMUONVTrackStore*
 265 AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
 266 {
 267   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
 268   /// each chamber (one when particle is entering + one when particle is leaving
 269   /// the sensitive volume)
 270
 271   Float_t maxGasGap = 1.; // cm
 272   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
 273   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
 274   TClonesArray *  hitForRecAtHit = 0;
 275   TClonesArray *  trackParamAtHit = 0;
 276   Float_t xRec,yRec,zRec;
 277   Float_t xRec1,yRec1,zRec1;
 278   Float_t xRec2,yRec2,zRec2;
 279   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
 280   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
 281   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
 282
 283   AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
 284   Int_t iHit1;
 285   Int_t iChamber = 0, detElemId = 0;
 286   Int_t nRec = 0;
 287   Int_t nTrackHits = 0;
 288
 289   hitForRec = new AliMUONHitForRec();
 290   trackParam = new AliMUONTrackParam();
 291
 292   TIter next(trackRefs.CreateIterator());
 293   AliMUONTrack* track;
 294
 295   while ( ( track = static_cast<AliMUONTrack*>(next())) )
 296   {
 297     hitForRecAtHit = track->GetHitForRecAtHit();
 298     trackParamAtHit = track->GetTrackParamAtHit();
 299     trackParamAtVertex = track->GetTrackParamAtVertex();
 300     nTrackHits = hitForRecAtHit->GetEntriesFast();
 301     AliMUONTrack trackNew;
 302     iHit1 = 0;
 303     while (iHit1 < nTrackHits)
 304     {
 305       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
 306       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
 307       xRec1  = hitForRec1->GetNonBendingCoor();
 308       yRec1  = hitForRec1->GetBendingCoor();
 309       zRec1  = hitForRec1->GetZ();
 310       xRec   = xRec1;
 311       yRec   = yRec1;
 312       zRec   = zRec1;
 313       bendingSlope1 = trackParam1->GetBendingSlope();
 314       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
 315       bendingMomentum1 = 0;
 316       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
 317         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
 318       bendingSlope = bendingSlope1;
 319       nonBendingSlope = nonBendingSlope1;
 320       bendingMomentum = bendingMomentum1;
 321       nRec = 1;
 322       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
 323       {
 324         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
 325         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
 326         xRec2  = hitForRec2->GetNonBendingCoor();
 327         yRec2  = hitForRec2->GetBendingCoor();
 328         zRec2  = hitForRec2->GetZ();
 329         bendingSlope2 = trackParam2->GetBendingSlope();
 330         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
 331         bendingMomentum2 = 0;
 332         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
 333           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
 334
 335         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
 336
 337           nRec++;
 338           xRec += xRec2;
 339           yRec += yRec2;
 340           zRec += zRec2;
 341           bendingSlope += bendingSlope2;
 342           nonBendingSlope += nonBendingSlope2;
 343           bendingMomentum += bendingMomentum2;
 344           iHit1 = iHit2;
 345         }
 346
 347       } // end iHit2
 348       xRec /= (Float_t)nRec;
 349       yRec /= (Float_t)nRec;
 350       zRec /= (Float_t)nRec;
 351       bendingSlope /= (Float_t)nRec;
 352       nonBendingSlope /= (Float_t)nRec;
 353       bendingMomentum /= (Float_t)nRec;
 354
 355       hitForRec->SetNonBendingCoor(xRec);
 356       hitForRec->SetBendingCoor(yRec);
 357       hitForRec->SetZ(zRec);
 358       detElemId = hitForRec->GetDetElemId();
 359       if (detElemId) iChamber = detElemId / 100 - 1;
 360       else iChamber = AliMUONConstants::ChamberNumber(zRec);
 361       hitForRec->SetChamberNumber(iChamber);
 362       hitForRec->SetBendingReso2(0.0);
 363       hitForRec->SetNonBendingReso2(0.0);
 364       trackParam->SetNonBendingCoor(xRec);
 365       trackParam->SetBendingCoor(yRec);
 366       trackParam->SetZ(zRec);
 367       trackParam->SetNonBendingSlope(nonBendingSlope);
 368       trackParam->SetBendingSlope(bendingSlope);
 369       if (TMath::Abs(bendingMomentum) > 0)
 370         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
 371
 372       trackNew.AddHitForRecAtHit(hitForRec);
 373       trackNew.AddTrackParamAtHit(trackParam,0);
 374
 375       iHit1++;
 376     } // end iHit1
 377
 378     trackNew.SetTrackID(track->GetTrackID());
 379     trackNew.SetTrackParamAtVertex(trackParamAtVertex);
 380     newMuonTrackRef->Add(trackNew);
 381
 382   } // end trackRef
 383
 384   delete hitForRec;
 385   delete trackParam;
 386   return newMuonTrackRef;
 387 }
 388
 389 //_____________________________________________________________________________
 390 AliMUONVTrackStore*
 391 AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
 392 {
 393   /// Calculate the number of reconstructible tracks
 394
 395   AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
 396
 397   TClonesArray* hitForRecAtHit = NULL;
 398   AliMUONHitForRec* hitForRec;
 399   Float_t zRec;
 400   Int_t nTrackHits;
 401   Int_t isChamberInTrack[10];
 402   Int_t iChamber = 0;
 403   Bool_t isTrackOK = kTRUE;
 404
 405   TIter next(trackRefs.CreateIterator());
 406   AliMUONTrack* track;
 407
 408   while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
 409   {
 410     hitForRecAtHit = track->GetHitForRecAtHit();
 411     nTrackHits = hitForRecAtHit->GetEntriesFast();
 412     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
 413
 414     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
 415       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
 416       zRec  = hitForRec->GetZ();
 417       iChamber = hitForRec->GetChamberNumber();
 418       if (iChamber < 0 || iChamber > 10) continue;
 419       isChamberInTrack[iChamber] = 1;
 420     }
 421     // track is reconstructible if the particle is depositing a hit
 422     // in the following chamber combinations:
 423
 424     isTrackOK = kTRUE;
 425     if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
 426     if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
 427     if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
 428     Int_t nHitsInLastStations=0;
 429     for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
 430       if (isChamberInTrack[ch]) nHitsInLastStations++;
 431     if(nHitsInLastStations < 3) isTrackOK = kFALSE;
 432
 433     if (isTrackOK)
 434     {
 435       reconstructibleStore->Add(*track);
 436     }
 437   }
 438
 439   return reconstructibleStore;
 440 }
 441