EMCAL/AliEMCALTracker.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 //                       Class AliEMCALTracker
  17 //                      -----------------------
  18 // Implementation of the track matching method between barrel tracks and
  19 // EMCAL clusters.
  20 // Besides algorithm implementation, some cuts are required to be set
  21 // in order to define, for each track, an acceptance window where clusters
  22 // are searched to find best match (if any).
  23 // The class accepts as input an ESD container, and works directly on it,
  24 // simply setting, for each of its tracks, the fEMCALindex flag, for each
  25 // track which is matched to a cluster.
  26 // In order to use method, one must launch PropagateBack().
  27 //
  28 // ------------------------------------------------------------------------
  29 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  30 // Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch)
  31 //=========================================================================
  32
  33 #include <Riostream.h>
  34 #include <iomanip>
  35
  36 #include <TFile.h>
  37 #include <TTree.h>
  38 #include <TList.h>
  39 #include <TString.h>
  40 #include <TVector3.h>
  41 #include <TClonesArray.h>
  42 #include <TGeoMatrix.h>
  43
  44 #include "AliLog.h"
  45 #include "AliESDEvent.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDCaloCluster.h"
  48 #include "AliEMCALRecPoint.h"
  49 #include "AliRunLoader.h"
  50 #include "AliEMCALTrack.h"
  51 #include "AliEMCALLoader.h"
  52 #include "AliEMCALGeometry.h"
  53 #include "AliEMCALReconstructor.h"
  54 #include "AliEMCALRecParam.h"
  55 #include "AliCDBEntry.h"
  56 #include "AliCDBManager.h"
  57 #include "AliEMCALReconstructor.h"
  58 #include "AliEMCALRecoUtils.h"
  59
  60 #include "AliEMCALTracker.h"
  61
  62 using std::cerr;
  63 using std::endl;
  64 ClassImp(AliEMCALTracker)
  65
  66 //
  67 //------------------------------------------------------------------------------
  68 //
  69 AliEMCALTracker::AliEMCALTracker() :
  70   AliTracker(),
  71   fCutPt(0),
  72   fCutNITS(0),
  73   fCutNTPC(50),
  74   fStep(20),
  75   fTrackCorrMode(kTrackCorrMMB),
  76   fEMCalSurfaceDistance(440),
  77   fClusterWindow(50),
  78   fCutEta(0.025),
  79   fCutPhi(0.05),
  80   fITSTrackSA(kFALSE),
  81   fTracks(0),
  82   fClusters(0),
  83   fGeom(0)
  84 {
  85   //
  86   // Default constructor.
  87   // Initializes all simple data members to default values,
  88    // and all collections to NULL.
  89   // Output file name is set to a default value.
  90   //
  91   InitParameters();
  92 }
  93 //
  94 //------------------------------------------------------------------------------
  95 //
  96 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) :
  97   AliTracker(),
  98   fCutPt(copy.fCutPt),
  99   fCutNITS(copy.fCutNITS),
 100   fCutNTPC(copy.fCutNTPC),
 101   fStep(copy.fStep),
 102   fTrackCorrMode(copy.fTrackCorrMode),
 103   fEMCalSurfaceDistance(copy.fEMCalSurfaceDistance),
 104   fClusterWindow(copy.fClusterWindow),
 105   fCutEta(copy.fCutEta),
 106   fCutPhi(copy.fCutPhi),
 107   fITSTrackSA(copy.fITSTrackSA),
 108   fTracks((TObjArray*)copy.fTracks->Clone()),
 109   fClusters((TObjArray*)copy.fClusters->Clone()),
 110   fGeom(copy.fGeom)
 111 {
 112   //
 113   // Copy constructor
 114   // Besides copying all parameters, duplicates all collections.
 115   //
 116 }
 117 //
 118 //------------------------------------------------------------------------------
 119 //
 120 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
 121 { // assignment operator; use copy ctor
 122   if (&source == this) return *this;
 123
 124   new (this) AliEMCALTracker(source);
 125   return *this;
 126 }
 127 //
 128 //------------------------------------------------------------------------------
 129 //
 130 void AliEMCALTracker::InitParameters()
 131 {
 132   //
 133   // Retrieve initialization parameters
 134   //
 135
 136   // Check if the instance of AliEMCALRecParam exists,
 137   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 138
 139   if (!recParam) {
 140     AliFatal("Reconstruction parameters for EMCAL not set!");
 141   } else {
 142     fCutEta  =  recParam->GetMthCutEta();
 143     fCutPhi  =  recParam->GetMthCutPhi();
 144     fStep    =  recParam->GetExtrapolateStep();
 145     fCutPt   =  recParam->GetTrkCutPt();
 146     fCutNITS =  recParam->GetTrkCutNITS();
 147     fCutNTPC =  recParam->GetTrkCutNTPC();
 148   }
 149 }
 150 //
 151 //------------------------------------------------------------------------------
 152 //
 153 void AliEMCALTracker::Clear(Option_t* option)
 154 {
 155   //
 156   // Clearing method
 157   // Deletes all objects in arrays and the arrays themselves
 158   //
 159
 160   TString opt(option);
 161   Bool_t clearTracks = opt.Contains("TRACKS");
 162   Bool_t clearClusters = opt.Contains("CLUSTERS");
 163   if (opt.Contains("ALL")) {
 164     clearTracks = kTRUE;
 165     clearClusters = kTRUE;
 166   }
 167
 168   //fTracks is a collection of esdTrack
 169   //When clearing this array, the linked objects should not be deleted
 170   if (fTracks != 0x0 && clearTracks) {
 171     fTracks->Clear();
 172     delete fTracks;
 173     fTracks = 0;
 174   }
 175   if (fClusters != 0x0 && clearClusters) {
 176     fClusters->Delete();
 177     delete fClusters;
 178     fClusters = 0;
 179   }
 180 }
 181 //
 182 //------------------------------------------------------------------------------
 183 //
 184 Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
 185 {
 186   //
 187   // Load EMCAL clusters in the form of AliEMCALRecPoint,
 188   // from simulation temporary files.
 189   // (When included in reconstruction chain, this method is used automatically)
 190   //
 191
 192   Clear("CLUSTERS");
 193
 194   cTree->SetBranchStatus("*",0); //disable all branches
 195   cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 196
 197   TBranch *branch = cTree->GetBranch("EMCALECARP");
 198   if (!branch) {
 199     AliError("Can't get the branch with the EMCAL clusters");
 200     return 1;
 201   }
 202
 203   TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
 204   branch->SetAddress(&clusters);
 205
 206   //cTree->GetEvent(0);
 207   branch->GetEntry(0);
 208   Int_t nClusters = (Int_t)clusters->GetEntries();
 209   if (fClusters) fClusters->Delete();
 210   else fClusters = new TObjArray(0);
 211   for (Int_t i = 0; i < nClusters; i++) {
 212     AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
 213     if (!cluster) continue;
 214     AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 215     fClusters->AddLast(matchCluster);
 216   }
 217
 218   branch->SetAddress(0);
 219   clusters->Delete();
 220   delete clusters;
 221
 222   AliDebug(1,Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
 223
 224   return 0;
 225 }
 226 //
 227 //------------------------------------------------------------------------------
 228 //
 229 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
 230 {
 231   //
 232   // Load EMCAL clusters in the form of AliESDCaloClusters,
 233   // from an AliESD object.
 234   //
 235
 236   // make sure that tracks/clusters collections are empty
 237   Clear("CLUSTERS");
 238   fClusters = new TObjArray(0);
 239
 240   Int_t nClusters = esd->GetNumberOfCaloClusters();
 241   for (Int_t i=0; i<nClusters; i++) {
 242     AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
 243     if (!cluster || !cluster->IsEMCAL()) continue ;
 244     AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 245     fClusters->AddLast(matchCluster);
 246   }
 247
 248   AliDebug(1,Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 249   return 0;
 250 }
 251 //
 252 //------------------------------------------------------------------------------
 253 //
 254 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
 255 {
 256   //
 257   // Load ESD tracks.
 258   //
 259
 260   UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
 261   UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
 262   Bool_t desc1 = (mask1 >> 3) & 0x1;
 263   Bool_t desc2 = (mask2 >> 3) & 0x1;
 264   if (desc1==0 || desc2==0) {
 265 //     AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
 266 //                   mask1, esd->GetESDRun()->GetDetectorsInReco(),
 267 //                   mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
 268     fITSTrackSA = kTRUE;
 269   }
 270
 271   Clear("TRACKS");
 272   fTracks = new TObjArray(0);
 273
 274   Int_t nTracks = esd->GetNumberOfTracks();
 275   //Bool_t isKink=kFALSE;
 276   for (Int_t i = 0; i < nTracks; i++) {
 277     AliESDtrack *esdTrack = esd->GetTrack(i);
 278     // set by default the value corresponding to "no match"
 279     esdTrack->SetEMCALcluster(kUnmatched);
 280     esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
 281
 282     //Select good quaulity tracks
 283     if (esdTrack->Pt()<fCutPt) continue;
 284     if (!fITSTrackSA)
 285       if (esdTrack->GetNcls(1)<fCutNTPC) continue;
 286
 287     //Loose geometric cut
 288     Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
 289     if (TMath::Abs(esdTrack->Eta())>0.9 || phi <= 10 || phi >= 250) continue;
 290     fTracks->AddLast(esdTrack);
 291   }
 292
 293   AliDebug(1,Form("Collected %d tracks", fTracks->GetEntries()));
 294   return 0;
 295 }
 296 //
 297 //------------------------------------------------------------------------------
 298 //
 299 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
 300 {
 301   //
 302   // Set track correction mode
 303   // gest the choice in string format and converts into
 304   // internal enum
 305   //
 306
 307   TString opt(option);
 308   opt.ToUpper();
 309
 310   if (!opt.CompareTo("NONE")) {
 311     fTrackCorrMode = kTrackCorrNone;
 312   } else if (!opt.CompareTo("MMB")) {
 313     fTrackCorrMode = kTrackCorrMMB;
 314   } else  {
 315     cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
 316   }
 317 }
 318 //
 319 //------------------------------------------------------------------------------
 320 //
 321 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
 322 {
 323   //
 324   // Main operation method.
 325   // Gets external AliESD containing tracks to be matched.
 326   // After executing match finding, stores in the same ESD object all infos
 327   // and releases the object for further reconstruction steps.
 328   //
 329   //
 330   // Note: should always return 0=OK, because otherwise all tracking
 331   // is aborted for this event
 332
 333   if (!esd)
 334   {
 335     AliError("NULL ESD passed");
 336     return 1;
 337   }
 338
 339   // step 1: collect clusters
 340   Int_t okLoadClusters = 0;
 341   if (!fClusters || (fClusters && fClusters->IsEmpty()))
 342     okLoadClusters = LoadClusters(esd);
 343
 344   Int_t nClusters = fClusters->GetEntries();
 345
 346   // step 2: collect ESD tracks
 347   Int_t nTracks, okLoadTracks;
 348   okLoadTracks = LoadTracks(esd);
 349   nTracks = fTracks->GetEntries();
 350
 351   AliDebug(5,Form("Propagate back %d tracks ok %d, for %d clusters ok %d",
 352                   nTracks,okLoadTracks,nClusters,okLoadClusters));
 353
 354   // step 3: for each track, find the closest cluster as matched within residual cuts
 355   Int_t index=-1;
 356   for (Int_t it = 0; it < nTracks; it++)
 357   {
 358     AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
 359     index = FindMatchedCluster(track);
 360     if (index>-1)
 361     {
 362       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
 363       track->SetEMCALcluster(cluster->Index());
 364       track->SetStatus(AliESDtrack::kEMCALmatch);
 365     }
 366   }
 367
 368   return 0;
 369 }
 370
 371 //
 372 //------------------------------------------------------------------------------
 373 //
 374 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
 375 {
 376   //
 377   // For each track, extrapolate it to all the clusters
 378   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
 379   //
 380
 381   Float_t maxEta=fCutEta;
 382   Float_t maxPhi=fCutPhi;
 383   Int_t index = -1;
 384
 385   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
 386   // Otherwise use the TPCInner point
 387   AliExternalTrackParam *trkParam = 0;
 388
 389   if (!fITSTrackSA)
 390   {
 391     const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
 392
 393     if (friendTrack && friendTrack->GetTPCOut())
 394       trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
 395     else if (track->GetInnerParam())
 396       trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
 397   }
 398   else
 399     trkParam = new AliExternalTrackParam(*track);
 400
 401   if (!trkParam) return index;
 402
 403   AliExternalTrackParam trkParamTmp(*trkParam);
 404   Float_t eta, phi, pt;
 405   if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))
 406   {
 407     if (fITSTrackSA) delete trkParam;
 408     return index;
 409   }
 410
 411   track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
 412
 413   if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
 414   {
 415     if (fITSTrackSA) delete trkParam;
 416     return index;
 417   }
 418
 419   //Perform extrapolation
 420   Double_t trkPos[3];
 421   trkParamTmp.GetXYZ(trkPos);
 422   Int_t nclusters = fClusters->GetEntries();
 423   for (Int_t ic=0; ic<nclusters; ic++)
 424   {
 425     AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 426
 427     Float_t clsPos[3] = {static_cast<Float_t>(cluster->X()),
 428                          static_cast<Float_t>(cluster->Y()),
 429                          static_cast<Float_t>(cluster->Z())};
 430
 431     Double_t dR = TMath::Sqrt(TMath::Power(trkPos[0]-clsPos[0],2)+TMath::Power(trkPos[1]-clsPos[1],2)+TMath::Power(trkPos[2]-clsPos[2],2));
 432     //printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
 433
 434     if (dR > fClusterWindow) continue;
 435
 436     AliExternalTrackParam trkParTmp(trkParamTmp);
 437
 438     Float_t tmpEta, tmpPhi;
 439     if (!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
 440
 441     if (TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
 442     {
 443       maxPhi=tmpPhi;
 444       maxEta=tmpEta;
 445       index=ic;
 446     }
 447   }
 448
 449   if (fITSTrackSA) delete trkParam;
 450
 451   return index;
 452 }
 453
 454 //
 455 //------------------------------------------------------------------------------
 456 //
 457 void AliEMCALTracker::UnloadClusters()
 458 {
 459   //
 460   // Free memory from all arrays
 461   // This method is called after the local tracking step
 462   // so we can safely delete everything
 463   //
 464
 465   Clear();
 466 }
 467
 468 //
 469 //------------------------------------------------------------------------------
 470 //
 471 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint) :
 472   fIndex(index),
 473   fX(0.),
 474   fY(0.),
 475   fZ(0.)
 476 {
 477   //
 478   // Translates an AliEMCALRecPoint object into the internal format.
 479   // Index of passed cluster in its native array must be specified.
 480   //
 481   TVector3 clpos;
 482   recPoint->GetGlobalPosition(clpos);
 483
 484   fX = clpos.X();
 485   fY = clpos.Y();
 486   fZ = clpos.Z();
 487 }
 488 //
 489 //------------------------------------------------------------------------------
 490 //
 491 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster) :
 492   fIndex(index),
 493   fX(0.),
 494   fY(0.),
 495   fZ(0.)
 496 {
 497   //
 498   // Translates an AliESDCaloCluster object into the internal format.
 499   // Index of passed cluster in its native array must be specified.
 500   //
 501   Float_t clpos[3]= {0., 0., 0.};
 502   caloCluster->GetPosition(clpos);
 503
 504   fX = (Double_t)clpos[0];
 505   fY = (Double_t)clpos[1];
 506   fZ = (Double_t)clpos[2];
 507 }