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