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