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