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