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