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