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