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