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