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