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