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   }
 242
 243   AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
 244
 245   delete trackRefStore;
 246
 247   return rv;
 248 }
 249
 250 //_____________________________________________________________________________
 251 Int_t
 252 AliMUONRecoCheck::NumberOfEvents() const
 253 {
 254   /// Return the number of events
 255   if ( fDataInterface )
 256   {
 257     return fDataInterface->NumberOfEvents();
 258   }
 259   return 0;
 260 }
 261
 262 //_____________________________________________________________________________
 263 AliMUONVTrackStore*
 264 AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
 265 {
 266   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
 267   /// each chamber (one when particle is entering + one when particle is leaving
 268   /// the sensitive volume)
 269
 270   Float_t maxGasGap = 1.; // cm
 271   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
 272   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
 273   TClonesArray *  hitForRecAtHit = 0;
 274   TClonesArray *  trackParamAtHit = 0;
 275   Float_t xRec,yRec,zRec;
 276   Float_t xRec1,yRec1,zRec1;
 277   Float_t xRec2,yRec2,zRec2;
 278   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
 279   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
 280   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
 281
 282   AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
 283   Int_t iHit1;
 284   Int_t iChamber = 0, detElemId = 0;
 285   Int_t nRec = 0;
 286   Int_t nTrackHits = 0;
 287
 288   hitForRec = new AliMUONHitForRec();
 289   trackParam = new AliMUONTrackParam();
 290
 291   TIter next(trackRefs.CreateIterator());
 292   AliMUONTrack* track;
 293
 294   while ( ( track = static_cast<AliMUONTrack*>(next())) )
 295   {
 296     hitForRecAtHit = track->GetHitForRecAtHit();
 297     trackParamAtHit = track->GetTrackParamAtHit();
 298     trackParamAtVertex = track->GetTrackParamAtVertex();
 299     nTrackHits = hitForRecAtHit->GetEntriesFast();
 300     AliMUONTrack trackNew;
 301     iHit1 = 0;
 302     while (iHit1 < nTrackHits)
 303     {
 304       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
 305       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
 306       xRec1  = hitForRec1->GetNonBendingCoor();
 307       yRec1  = hitForRec1->GetBendingCoor();
 308       zRec1  = hitForRec1->GetZ();
 309       xRec   = xRec1;
 310       yRec   = yRec1;
 311       zRec   = zRec1;
 312       bendingSlope1 = trackParam1->GetBendingSlope();
 313       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
 314       bendingMomentum1 = 0;
 315       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
 316         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
 317       bendingSlope = bendingSlope1;
 318       nonBendingSlope = nonBendingSlope1;
 319       bendingMomentum = bendingMomentum1;
 320       nRec = 1;
 321       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
 322       {
 323         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
 324         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
 325         xRec2  = hitForRec2->GetNonBendingCoor();
 326         yRec2  = hitForRec2->GetBendingCoor();
 327         zRec2  = hitForRec2->GetZ();
 328         bendingSlope2 = trackParam2->GetBendingSlope();
 329         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
 330         bendingMomentum2 = 0;
 331         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
 332           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
 333
 334         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
 335
 336           nRec++;
 337           xRec += xRec2;
 338           yRec += yRec2;
 339           zRec += zRec2;
 340           bendingSlope += bendingSlope2;
 341           nonBendingSlope += nonBendingSlope2;
 342           bendingMomentum += bendingMomentum2;
 343           iHit1 = iHit2;
 344         }
 345
 346       } // end iHit2
 347       xRec /= (Float_t)nRec;
 348       yRec /= (Float_t)nRec;
 349       zRec /= (Float_t)nRec;
 350       bendingSlope /= (Float_t)nRec;
 351       nonBendingSlope /= (Float_t)nRec;
 352       bendingMomentum /= (Float_t)nRec;
 353
 354       hitForRec->SetNonBendingCoor(xRec);
 355       hitForRec->SetBendingCoor(yRec);
 356       hitForRec->SetZ(zRec);
 357       detElemId = hitForRec->GetDetElemId();
 358       if (detElemId) iChamber = detElemId / 100 - 1;
 359       else iChamber = AliMUONConstants::ChamberNumber(zRec);
 360       hitForRec->SetChamberNumber(iChamber);
 361       hitForRec->SetBendingReso2(0.0);
 362       hitForRec->SetNonBendingReso2(0.0);
 363       trackParam->SetNonBendingCoor(xRec);
 364       trackParam->SetBendingCoor(yRec);
 365       trackParam->SetZ(zRec);
 366       trackParam->SetNonBendingSlope(nonBendingSlope);
 367       trackParam->SetBendingSlope(bendingSlope);
 368       if (TMath::Abs(bendingMomentum) > 0)
 369         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
 370
 371       trackNew.AddHitForRecAtHit(hitForRec);
 372       trackNew.AddTrackParamAtHit(trackParam,0);
 373
 374       iHit1++;
 375     } // end iHit1
 376
 377     trackNew.SetTrackID(track->GetTrackID());
 378     trackNew.SetTrackParamAtVertex(trackParamAtVertex);
 379     newMuonTrackRef->Add(trackNew);
 380
 381   } // end trackRef
 382
 383   delete hitForRec;
 384   delete trackParam;
 385   return newMuonTrackRef;
 386 }
 387
 388 //_____________________________________________________________________________
 389 AliMUONVTrackStore*
 390 AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
 391 {
 392   /// Calculate the number of reconstructible tracks
 393
 394   AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
 395
 396   TClonesArray* hitForRecAtHit = NULL;
 397   AliMUONHitForRec* hitForRec;
 398   Float_t zRec;
 399   Int_t nTrackHits;
 400   Int_t isChamberInTrack[10];
 401   Int_t iChamber = 0;
 402   Bool_t isTrackOK = kTRUE;
 403
 404   TIter next(trackRefs.CreateIterator());
 405   AliMUONTrack* track;
 406
 407   while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
 408   {
 409     hitForRecAtHit = track->GetHitForRecAtHit();
 410     nTrackHits = hitForRecAtHit->GetEntriesFast();
 411     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
 412
 413     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
 414       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
 415       zRec  = hitForRec->GetZ();
 416       iChamber = hitForRec->GetChamberNumber();
 417       if (iChamber < 0 || iChamber > 10) continue;
 418       isChamberInTrack[iChamber] = 1;
 419     }
 420     // track is reconstructible if the particle is depositing a hit
 421     // in the following chamber combinations:
 422
 423     isTrackOK = kTRUE;
 424     if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
 425     if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
 426     if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
 427     Int_t nHitsInLastStations=0;
 428     for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
 429       if (isChamberInTrack[ch]) nHitsInLastStations++;
 430     if(nHitsInLastStations < 3) isTrackOK = kFALSE;
 431
 432     if (isTrackOK)
 433     {
 434       reconstructibleStore->Add(*track);
 435     }
 436   }
 437
 438   return reconstructibleStore;
 439 }
 440