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   }
 142   else{
 143
 144     fCutEta  =  recParam->GetMthCutEta();
 145     fCutPhi  =  recParam->GetMthCutPhi();
 146     fStep    =  recParam->GetExtrapolateStep();
 147     fCutPt   =  recParam->GetTrkCutPt();
 148     fCutNITS =  recParam->GetTrkCutNITS();
 149     fCutNTPC =  recParam->GetTrkCutNTPC();
 150   }
 151
 152 }
 153
 154 //
 155 //------------------------------------------------------------------------------
 156 //
 157 void AliEMCALTracker::Clear(Option_t* option)
 158 {
 159         //
 160         // Clearing method
 161         // Deletes all objects in arrays and the arrays themselves
 162         //
 163
 164         TString opt(option);
 165         Bool_t clearTracks = opt.Contains("TRACKS");
 166         Bool_t clearClusters = opt.Contains("CLUSTERS");
 167         if (opt.Contains("ALL")) {
 168                 clearTracks = kTRUE;
 169                 clearClusters = kTRUE;
 170         }
 171
 172         //fTracks is a collection of esdTrack
 173         //When clearing this array, the linked objects should not be deleted
 174         if (fTracks != 0x0 && clearTracks) {
 175            fTracks->Clear();
 176            delete fTracks;
 177            fTracks = 0;
 178         }
 179         if (fClusters != 0x0 && clearClusters) {
 180            fClusters->Delete();
 181            delete fClusters;
 182            fClusters = 0;
 183         }
 184 }
 185 //
 186 //------------------------------------------------------------------------------
 187 //
 188 Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
 189 {
 190         //
 191         // Load EMCAL clusters in the form of AliEMCALRecPoint,
 192         // from simulation temporary files.
 193         // (When included in reconstruction chain, this method is used automatically)
 194         //
 195
 196         Clear("CLUSTERS");
 197
 198         cTree->SetBranchStatus("*",0); //disable all branches
 199         cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 200
 201         TBranch *branch = cTree->GetBranch("EMCALECARP");
 202         if (!branch) {
 203                 AliError("Can't get the branch with the EMCAL clusters");
 204                 return 1;
 205         }
 206
 207         TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
 208         branch->SetAddress(&clusters);
 209
 210         //cTree->GetEvent(0);
 211         branch->GetEntry(0);
 212         Int_t nClusters = (Int_t)clusters->GetEntries();
 213         if(fClusters) fClusters->Delete();
 214         else fClusters = new TObjArray(0);
 215         for (Int_t i = 0; i < nClusters; i++) {
 216                 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
 217                 if (!cluster) continue;
 218                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 219                 fClusters->AddLast(matchCluster);
 220         }
 221
 222         branch->SetAddress(0);
 223         clusters->Delete();
 224         delete clusters;
 225
 226         AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
 227
 228         return 0;
 229 }
 230 //
 231 //------------------------------------------------------------------------------
 232 //
 233 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
 234 {
 235   //
 236   // Load EMCAL clusters in the form of AliESDCaloClusters,
 237   // from an AliESD object.
 238   //
 239
 240   // make sure that tracks/clusters collections are empty
 241   Clear("CLUSTERS");
 242   fClusters = new TObjArray(0);
 243
 244   Int_t nClusters = esd->GetNumberOfCaloClusters();
 245   for (Int_t i=0; i<nClusters; i++)
 246     {
 247       AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
 248       if (!cluster || !cluster->IsEMCAL()) continue ;
 249       AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 250       fClusters->AddLast(matchCluster);
 251     }
 252
 253   AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 254   return 0;
 255 }
 256 //
 257 //------------------------------------------------------------------------------
 258 //
 259 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
 260 {
 261   //
 262   // Load ESD tracks.
 263   //
 264
 265   UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
 266   UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
 267   Bool_t desc1 = (mask1 >> 3) & 0x1;
 268   Bool_t desc2 = (mask2 >> 3) & 0x1;
 269   if (desc1==0 || desc2==0) {
 270 //     AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
 271 //                   mask1, esd->GetESDRun()->GetDetectorsInReco(),
 272 //                   mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
 273     fITSTrackSA = kTRUE;
 274   }
 275
 276   Clear("TRACKS");
 277   fTracks = new TObjArray(0);
 278
 279   Int_t nTracks = esd->GetNumberOfTracks();
 280   //Bool_t isKink=kFALSE;
 281   for (Int_t i = 0; i < nTracks; i++)
 282     {
 283       AliESDtrack *esdTrack = esd->GetTrack(i);
 284       // set by default the value corresponding to "no match"
 285       esdTrack->SetEMCALcluster(kUnmatched);
 286       esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
 287
 288       //Select good quaulity tracks
 289       if(esdTrack->Pt()<fCutPt) continue;
 290       if(!fITSTrackSA)
 291         if(esdTrack->GetNcls(1)<fCutNTPC)continue;
 292
 293       //Loose geometric cut
 294       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
 295       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
 296
 297       fTracks->AddLast(esdTrack);
 298     }
 299
 300       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
 301       return 0;
 302 }
 303 //
 304 //------------------------------------------------------------------------------
 305 //
 306 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
 307 {
 308   //
 309   // Set track correction mode
 310   // gest the choice in string format and converts into
 311   // internal enum
 312   //
 313
 314   TString opt(option);
 315   opt.ToUpper();
 316
 317   if (!opt.CompareTo("NONE"))
 318     {
 319       fTrackCorrMode = kTrackCorrNone;
 320     }
 321   else if (!opt.CompareTo("MMB"))
 322     {
 323       fTrackCorrMode = kTrackCorrMMB;
 324     }
 325   else
 326     {
 327       cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
 328     }
 329 }
 330 //
 331 //------------------------------------------------------------------------------
 332 //
 333 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
 334 {
 335         //
 336         // Main operation method.
 337         // Gets external AliESD containing tracks to be matched.
 338         // After executing match finding, stores in the same ESD object all infos
 339         // and releases the object for further reconstruction steps.
 340         //
 341         //
 342         // Note: should always return 0=OK, because otherwise all tracking
 343         // is aborted for this event
 344
 345         if (!esd) {
 346                 AliError("NULL ESD passed");
 347                 return 1;
 348         }
 349
 350         // step 1: collect clusters
 351         Int_t okLoadClusters, nClusters;
 352         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
 353                 okLoadClusters = LoadClusters(esd);
 354         }
 355         nClusters = fClusters->GetEntries();
 356
 357         // step 2: collect ESD tracks
 358         Int_t nTracks, okLoadTracks;
 359         okLoadTracks = LoadTracks(esd);
 360         nTracks = fTracks->GetEntries();
 361
 362         // step 3: for each track, find the closest cluster as matched within residual cuts
 363         Int_t index=-1;
 364         for (Int_t it = 0; it < nTracks; it++)
 365           {
 366             AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
 367             index = FindMatchedCluster(track);
 368             if (index>-1)
 369               {
 370                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
 371                 track->SetEMCALcluster(cluster->Index());
 372                 track->SetStatus(AliESDtrack::kEMCALmatch);
 373               }
 374           }
 375
 376         return 0;
 377 }
 378
 379 //
 380 //------------------------------------------------------------------------------
 381 //
 382 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
 383 {
 384   //
 385   // For each track, extrapolate it to all the clusters
 386   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
 387   //
 388
 389   Float_t maxEta=fCutEta;
 390   Float_t maxPhi=fCutPhi;
 391   Int_t index = -1;
 392
 393   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
 394   // Otherwise use the TPCInner point
 395   AliExternalTrackParam *trkParam = 0;
 396
 397   if(!fITSTrackSA){
 398     const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
 399     if(friendTrack && friendTrack->GetTPCOut())
 400       trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
 401     else if(track->GetInnerParam())
 402       trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
 403   }
 404   else
 405     trkParam = new AliExternalTrackParam(*track);
 406
 407   if(!trkParam) return index;
 408
 409   AliExternalTrackParam trkParamTmp(*trkParam);
 410   Float_t eta, phi, pt;
 411   if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))  {
 412         if(fITSTrackSA) delete trkParam;
 413         return index;
 414   }
 415   track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
 416   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
 417          if(fITSTrackSA) delete trkParam;
 418         return index;
 419   }
 420
 421   //Perform extrapolation
 422   Double_t trkPos[3];
 423   trkParamTmp.GetXYZ(trkPos);
 424   Int_t nclusters = fClusters->GetEntries();
 425   for(Int_t ic=0; ic<nclusters; ic++)
 426     {
 427       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 428       Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
 429       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));
 430 //       printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
 431       if(dR > fClusterWindow) continue;
 432
 433       AliExternalTrackParam trkParTmp(trkParamTmp);
 434
 435       Float_t tmpEta, tmpPhi;
 436       if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
 437       if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
 438         {
 439           maxPhi=tmpPhi;
 440           maxEta=tmpEta;
 441           index=ic;
 442         }
 443       }
 444
 445   if(fITSTrackSA) delete trkParam;
 446   return index;
 447 }
 448
 449 //
 450 //------------------------------------------------------------------------------
 451 //
 452 void AliEMCALTracker::UnloadClusters()
 453 {
 454         //
 455         // Free memory from all arrays
 456         // This method is called after the local tracking step
 457         // so we can safely delete everything
 458         //
 459
 460         Clear();
 461 }
 462
 463 //
 464 //------------------------------------------------------------------------------
 465 //
 466 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
 467   : fIndex(index),
 468     fX(0.),
 469     fY(0.),
 470     fZ(0.)
 471 {
 472         //
 473         // Translates an AliEMCALRecPoint object into the internal format.
 474         // Index of passed cluster in its native array must be specified.
 475         //
 476         TVector3 clpos;
 477         recPoint->GetGlobalPosition(clpos);
 478
 479         fX = clpos.X();
 480         fY = clpos.Y();
 481         fZ = clpos.Z();
 482 }
 483 //
 484 //------------------------------------------------------------------------------
 485 //
 486 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
 487   : fIndex(index),
 488     fX(0.),
 489     fY(0.),
 490     fZ(0.)
 491 {
 492         //
 493         // Translates an AliESDCaloCluster object into the internal format.
 494         // Index of passed cluster in its native array must be specified.
 495         //
 496         Float_t clpos[3]= {0., 0., 0.};
 497         caloCluster->GetPosition(clpos);
 498
 499         fX = (Double_t)clpos[0];
 500         fY = (Double_t)clpos[1];
 501         fZ = (Double_t)clpos[2];
 502 }