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 = 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;
 115
 116   trackParam = new AliMUONTrackParam();
 117   hitForRec = new AliMUONHitForRec();
 118   muonTrack = new AliMUONTrack();
 119
 120   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
 121     branch->GetEntry(iTrackRef);
 122
 123     iHitMin = 0;
 124     isNewTrack = kTRUE;
 125
 126     if (!trackRefs->GetEntries()) continue;
 127
 128     while (isNewTrack) {
 129
 130       for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
 131
 132         trackReference = (AliTrackReference*)trackRefs->At(iHit);
 133         x = trackReference->X();
 134         y = trackReference->Y();
 135         z = trackReference->Z();
 136         pX = trackReference->Px();
 137         pY = trackReference->Py();
 138         pZ = trackReference->Pz();
 139
 140         track = trackReference->GetTrack();
 141
 142         if (track != trackSave && iHit != 0) {
 143           iHitMin = iHit;
 144           trackSave = track;
 145           break;
 146         }
 147
 148         // track parameters at hit
 149         trackParam->SetBendingCoor(y);
 150         trackParam->SetNonBendingCoor(x);
 151         trackParam->SetZ(z);
 152
 153         if (TMath::Abs(pZ) > 0) {
 154           bendingSlope = pY/pZ;
 155           nonBendingSlope = pX/pZ;
 156         }
 157         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 158         if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 159
 160         trackParam->SetBendingSlope(bendingSlope);
 161         trackParam->SetNonBendingSlope(nonBendingSlope);
 162         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
 163
 164         hitForRec->SetBendingCoor(y);
 165         hitForRec->SetNonBendingCoor(x);
 166         hitForRec->SetZ(z);
 167         hitForRec->SetBendingReso2(0.0);
 168         hitForRec->SetNonBendingReso2(0.0);
 169         iChamber = AliMUONConstants::ChamberNumber(z);
 170         hitForRec->SetChamberNumber(iChamber);
 171
 172         muonTrack->AddTrackParamAtHit(trackParam);
 173         muonTrack->AddHitForRecAtHit(hitForRec);
 174         muonTrack->SetTrackID(track);
 175
 176         trackSave = track;
 177         if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
 178       }
 179
 180       // track parameters at vertex
 181       particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
 182       if (particle) {
 183         x = particle->Vx();
 184         y = particle->Vy();
 185         z = particle->Vz();
 186         pX = particle->Px();
 187         pY = particle->Py();
 188         pZ = particle->Pz();
 189
 190         trackParam->SetBendingCoor(y);
 191         trackParam->SetNonBendingCoor(x);
 192         trackParam->SetZ(z);
 193         if (TMath::Abs(pZ) > 0) {
 194           bendingSlope = pY/pZ;
 195           nonBendingSlope = pX/pZ;
 196         }
 197         pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
 198         if (pYZ >0) inverseBendingMomentum = 1/pYZ;
 199         trackParam->SetBendingSlope(bendingSlope);
 200         trackParam->SetNonBendingSlope(nonBendingSlope);
 201         trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
 202
 203         muonTrack->SetTrackParamAtVertex(trackParam);
 204       }
 205
 206       AddMuonTrackReference(muonTrack);
 207       muonTrack->ResetTrackParamAtHit();
 208       muonTrack->ResetHitForRecAtHit();
 209
 210     } // end while isNewTrack
 211
 212   }
 213
 214   CleanMuonTrackRef();
 215
 216   ReconstructibleTracks();
 217
 218   delete muonTrack;
 219   delete trackParam;
 220   delete hitForRec;
 221
 222 }
 223
 224 //____________________________________________________________________________
 225 TClonesArray* AliMUONRecoCheck::GetTrackReco()
 226 {
 227   // Return TClonesArray of reconstructed tracks
 228
 229   GetMUONData()->ResetRecTracks();
 230   GetMUONData()->SetTreeAddress("RT");
 231   fTrackReco = GetMUONData()->RecTracks();
 232   GetMUONData()->GetRecTracks();
 233   fRecoTracks = fTrackReco->GetEntriesFast();
 234   return fTrackReco;
 235 }
 236 //_____________________________________________________________________________
 237 void AliMUONRecoCheck::PrintEvent() const
 238 {
 239   // debug facility
 240
 241   AliMUONTrack *track;
 242   AliMUONHitForRec *hitForRec;
 243   TClonesArray *  hitForRecAtHit = 0;
 244   Float_t xRec,yRec,zRec;
 245
 246   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
 247
 248   printf(" ******************************* \n");
 249   printf(" nb of tracks %d \n",nTrackRef);
 250
 251   for (Int_t index = 0; index < nTrackRef; index++) {
 252     track = (AliMUONTrack*)fMuonTrackRef->At(index);
 253     hitForRecAtHit = track->GetHitForRecAtHit();
 254     Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
 255     printf(" track number %d \n",index);
 256     for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
 257       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
 258       xRec  = hitForRec->GetNonBendingCoor();
 259       yRec  = hitForRec->GetBendingCoor();
 260       zRec  = hitForRec->GetZ();
 261       printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
 262     }
 263   }
 264 }
 265
 266 //_____________________________________________________________________________
 267 void AliMUONRecoCheck::ResetTracks() const
 268 {
 269   if (fMuonTrackRef) fMuonTrackRef->Clear();
 270 }
 271 //_____________________________________________________________________________
 272 void AliMUONRecoCheck::CleanMuonTrackRef()
 273 {
 274   // Re-calculate hits parameters because two AliTrackReferences are recorded for
 275   // each chamber (one when particle is entering + one when particle is leaving
 276   // the sensitive volume)
 277
 278   Float_t maxGasGap = 1.; // cm
 279   AliMUONTrack *track, *trackNew;
 280   AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
 281   AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
 282   TClonesArray *  hitForRecAtHit = 0;
 283   TClonesArray *  trackParamAtHit = 0;
 284   Float_t xRec,yRec,zRec;
 285   Float_t xRec1,yRec1,zRec1;
 286   Float_t xRec2,yRec2,zRec2;
 287   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
 288   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
 289   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
 290   TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
 291   Int_t iHit1;
 292   Int_t iChamber = 0;
 293   Int_t nRec = 0;
 294   Int_t nTrackHits = 0;
 295
 296   hitForRec = new AliMUONHitForRec();
 297   trackParam = new AliMUONTrackParam();
 298   trackNew = new AliMUONTrack();
 299
 300   Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
 301
 302   for (Int_t index = 0; index < nTrackRef; index++) {
 303     track = (AliMUONTrack*)fMuonTrackRef->At(index);
 304     hitForRecAtHit = track->GetHitForRecAtHit();
 305     trackParamAtHit = track->GetTrackParamAtHit();
 306     trackParamAtVertex = track->GetTrackParamAtVertex();
 307     nTrackHits = hitForRecAtHit->GetEntriesFast();
 308     iHit1 = 0;
 309     while (iHit1 < nTrackHits) {
 310       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
 311       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
 312       xRec1  = hitForRec1->GetNonBendingCoor();
 313       yRec1  = hitForRec1->GetBendingCoor();
 314       zRec1  = hitForRec1->GetZ();
 315       xRec   = xRec1;
 316       yRec   = yRec1;
 317       zRec   = zRec1;
 318       bendingSlope1 = trackParam1->GetBendingSlope();
 319       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
 320       bendingMomentum1 = 0;
 321       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
 322         bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
 323       bendingSlope = bendingSlope1;
 324       nonBendingSlope = nonBendingSlope1;
 325       bendingMomentum = bendingMomentum1;
 326       nRec = 1;
 327       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
 328         hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
 329         trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
 330         xRec2  = hitForRec2->GetNonBendingCoor();
 331         yRec2  = hitForRec2->GetBendingCoor();
 332         zRec2  = hitForRec2->GetZ();
 333         bendingSlope2 = trackParam2->GetBendingSlope();
 334         nonBendingSlope2 = trackParam2->GetNonBendingSlope();
 335         bendingMomentum2 = 0;
 336         if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
 337           bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
 338
 339         if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
 340           nRec++;
 341           xRec += xRec2;
 342           yRec += yRec2;
 343           zRec += zRec2;
 344           bendingSlope += bendingSlope2;
 345           nonBendingSlope += nonBendingSlope2;
 346           bendingMomentum += bendingMomentum2;
 347           iHit1 = iHit2;
 348         }
 349
 350       } // end iHit2
 351       xRec /= (Float_t)nRec;
 352       yRec /= (Float_t)nRec;
 353       zRec /= (Float_t)nRec;
 354       bendingSlope /= (Float_t)nRec;
 355       nonBendingSlope /= (Float_t)nRec;
 356       bendingMomentum /= (Float_t)nRec;
 357
 358       hitForRec->SetNonBendingCoor(xRec);
 359       hitForRec->SetBendingCoor(yRec);
 360       hitForRec->SetZ(zRec);
 361       iChamber = AliMUONConstants::ChamberNumber(zRec);
 362       hitForRec->SetChamberNumber(iChamber);
 363       hitForRec->SetBendingReso2(0.0);
 364       hitForRec->SetNonBendingReso2(0.0);
 365       trackParam->SetNonBendingCoor(xRec);
 366       trackParam->SetBendingCoor(yRec);
 367       trackParam->SetZ(zRec);
 368       trackParam->SetNonBendingSlope(nonBendingSlope);
 369       trackParam->SetBendingSlope(bendingSlope);
 370       if (TMath::Abs(bendingMomentum) > 0)
 371         trackParam->SetInverseBendingMomentum(1./bendingMomentum);
 372
 373       trackNew->AddHitForRecAtHit(hitForRec);
 374       trackNew->AddTrackParamAtHit(trackParam);
 375
 376       iHit1++;
 377     } // end iHit1
 378
 379     trackNew->SetTrackID(track->GetTrackID());
 380     trackNew->SetTrackParamAtVertex(trackParamAtVertex);
 381     {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
 382     trackNew->ResetHitForRecAtHit();
 383     trackNew->ResetTrackParamAtHit();
 384
 385   } // end trackRef
 386
 387   fMuonTrackRef->Clear();
 388   nTrackRef = newMuonTrackRef->GetEntriesFast();
 389   for (Int_t index = 0; index < nTrackRef; index++) {
 390     track = (AliMUONTrack*)newMuonTrackRef->At(index);
 391     AddMuonTrackReference(track);
 392   }
 393
 394
 395   delete trackNew;
 396   delete hitForRec;
 397   delete trackParam;
 398   newMuonTrackRef->Delete();
 399   delete newMuonTrackRef;
 400
 401 }
 402
 403 //_____________________________________________________________________________
 404 void AliMUONRecoCheck::ReconstructibleTracks()
 405 {
 406   // calculate the number of reconstructible tracks
 407
 408   AliMUONTrack* track;
 409   TClonesArray* hitForRecAtHit = NULL;
 410   AliMUONHitForRec* hitForRec;
 411   Float_t zRec;
 412   Int_t nTrackRef, nTrackHits;
 413   Int_t isChamberInTrack[10];
 414   Int_t iChamber = 0;
 415   Bool_t isTrackOK = kTRUE;
 416
 417   fReconstructibleTracks = 0;
 418
 419   nTrackRef = fMuonTrackRef->GetEntriesFast();
 420   for (Int_t index = 0; index < nTrackRef; index++) {
 421     track = (AliMUONTrack*)fMuonTrackRef->At(index);
 422     hitForRecAtHit = track->GetHitForRecAtHit();
 423     nTrackHits = hitForRecAtHit->GetEntriesFast();
 424     for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
 425
 426     for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
 427       hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
 428       zRec  = hitForRec->GetZ();
 429       iChamber = hitForRec->GetChamberNumber();
 430       if (iChamber < 0 || iChamber > 10) continue;
 431       isChamberInTrack[iChamber] = 1;
 432
 433     }
 434     // track is reconstructible if the particle is crossing every tracking chambers
 435     isTrackOK = kTRUE;
 436     for (Int_t ch = 0; ch < 10; ch++) {
 437       if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
 438     }    if (isTrackOK) fReconstructibleTracks++;
 439     if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
 440   }
 441   fMuonTrackRef->Compress();
 442 }
 443
 444