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   fClusterWindow(50),
  77   fCutEta(0.025),
  78   fCutPhi(0.05),
  79   fTracks(0),
  80   fClusters(0),
  81   fGeom(0)
  82 {
  83   //
  84   // Default constructor.
  85   // Initializes all simple data members to default values,
  86    // and all collections to NULL.
  87   // Output file name is set to a default value.
  88   //
  89   InitParameters();
  90 }
  91 //
  92 //------------------------------------------------------------------------------
  93 //
  94 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
  95   : AliTracker(),
  96     fCutPt(copy.fCutPt),
  97     fCutNITS(copy.fCutNITS),
  98     fCutNTPC(copy.fCutNTPC),
  99     fStep(copy.fStep),
 100     fTrackCorrMode(copy.fTrackCorrMode),
 101     fClusterWindow(copy.fClusterWindow),
 102     fCutEta(copy.fCutEta),
 103     fCutPhi(copy.fCutPhi),
 104     fTracks((TObjArray*)copy.fTracks->Clone()),
 105     fClusters((TObjArray*)copy.fClusters->Clone()),
 106     fGeom(copy.fGeom)
 107 {
 108   //
 109   // Copy constructor
 110   // Besides copying all parameters, duplicates all collections.
 111   //
 112 }
 113 //
 114 //------------------------------------------------------------------------------
 115 //
 116 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
 117 { // assignment operator; use copy ctor
 118   if (&source == this) return *this;
 119
 120   new (this) AliEMCALTracker(source);
 121   return *this;
 122 }
 123 //
 124 //------------------------------------------------------------------------------
 125 //
 126 void AliEMCALTracker::InitParameters()
 127 {
 128   //
 129   // Retrieve initialization parameters
 130   //
 131
 132   // Check if the instance of AliEMCALRecParam exists,
 133   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 134
 135   if(!recParam){
 136     AliFatal("Reconstruction parameters for EMCAL not set!");
 137   }
 138   else{
 139
 140     fCutEta  =  recParam->GetMthCutEta();
 141     fCutPhi  =  recParam->GetMthCutPhi();
 142     fStep    =  recParam->GetExtrapolateStep();
 143     fCutPt   =  recParam->GetTrkCutPt();
 144     fCutNITS =  recParam->GetTrkCutNITS();
 145     fCutNTPC =  recParam->GetTrkCutNTPC();
 146   }
 147
 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         AliInfo(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     {
 243       AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
 244       if (!cluster || !cluster->IsEMCAL()) continue ;
 245       AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 246       fClusters->AddLast(matchCluster);
 247     }
 248
 249   AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 250   return 0;
 251 }
 252 //
 253 //------------------------------------------------------------------------------
 254 //
 255 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
 256 {
 257   //
 258   // Load ESD tracks.
 259   //
 260
 261   Clear("TRACKS");
 262   fTracks = new TObjArray(0);
 263
 264   Int_t nTracks = esd->GetNumberOfTracks();
 265   //Bool_t isKink=kFALSE;
 266   for (Int_t i = 0; i < nTracks; i++)
 267     {
 268       AliESDtrack *esdTrack = esd->GetTrack(i);
 269       // set by default the value corresponding to "no match"
 270       esdTrack->SetEMCALcluster(kUnmatched);
 271       esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
 272
 273       //Select good quaulity tracks
 274       if(esdTrack->Pt()<fCutPt) continue;
 275       if(esdTrack->GetNcls(1)<fCutNTPC)continue;
 276
 277       //Loose geometric cut
 278       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
 279       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
 280
 281       fTracks->AddLast(esdTrack);
 282     }
 283
 284       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
 285       return 0;
 286 }
 287 //
 288 //------------------------------------------------------------------------------
 289 //
 290 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
 291 {
 292   //
 293   // Set track correction mode
 294   // gest the choice in string format and converts into
 295   // internal enum
 296   //
 297
 298   TString opt(option);
 299   opt.ToUpper();
 300
 301   if (!opt.CompareTo("NONE"))
 302     {
 303       fTrackCorrMode = kTrackCorrNone;
 304     }
 305   else if (!opt.CompareTo("MMB"))
 306     {
 307       fTrackCorrMode = kTrackCorrMMB;
 308     }
 309   else
 310     {
 311       cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
 312     }
 313 }
 314 //
 315 //------------------------------------------------------------------------------
 316 //
 317 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
 318 {
 319         //
 320         // Main operation method.
 321         // Gets external AliESD containing tracks to be matched.
 322         // After executing match finding, stores in the same ESD object all infos
 323         // and releases the object for further reconstruction steps.
 324         //
 325         //
 326         // Note: should always return 0=OK, because otherwise all tracking
 327         // is aborted for this event
 328
 329         if (!esd) {
 330                 AliError("NULL ESD passed");
 331                 return 1;
 332         }
 333
 334         // step 1: collect clusters
 335         Int_t okLoadClusters, nClusters;
 336         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
 337                 okLoadClusters = LoadClusters(esd);
 338         }
 339         nClusters = fClusters->GetEntries();
 340
 341         // step 2: collect ESD tracks
 342         Int_t nTracks, okLoadTracks;
 343         okLoadTracks = LoadTracks(esd);
 344         nTracks = fTracks->GetEntries();
 345
 346         // step 3: for each track, find the closest cluster as matched within residual cuts
 347         Int_t index=-1;
 348         for (Int_t it = 0; it < nTracks; it++)
 349           {
 350             AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
 351             index = FindMatchedCluster(track);
 352             if (index>-1)
 353               {
 354                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
 355                 track->SetEMCALcluster(cluster->Index());
 356                 track->SetStatus(AliESDtrack::kEMCALmatch);
 357               }
 358           }
 359
 360         return 0;
 361 }
 362
 363 //
 364 //------------------------------------------------------------------------------
 365 //
 366 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
 367 {
 368   //
 369   // For each track, extrapolate it to all the clusters
 370   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
 371   //
 372
 373   Float_t maxEta=fCutEta;
 374   Float_t maxPhi=fCutPhi;
 375   Int_t index = -1;
 376
 377   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
 378   // Otherwise use the TPCInner point
 379   AliExternalTrackParam *trkParam = 0;
 380   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
 381   if(friendTrack && friendTrack->GetTPCOut())
 382     trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
 383   else
 384     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
 385   if(!trkParam) return index;
 386
 387
 388   AliExternalTrackParam trkParamTmp(*trkParam);
 389   Float_t eta, phi;
 390   if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(kTRUE), fStep, eta, phi))  return index;
 391   track->SetTrackPhiEtaOnEMCal(phi,eta);
 392   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
 393
 394   //Perform extrapolation
 395   Double_t trkPos[3];
 396   trkParamTmp.GetXYZ(trkPos);
 397   Int_t nclusters = fClusters->GetEntries();
 398   for(Int_t ic=0; ic<nclusters; ic++)
 399     {
 400       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 401       Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
 402       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));
 403       if(dR > fClusterWindow) continue;
 404
 405       AliExternalTrackParam trkParTmp(trkParamTmp);
 406
 407       Float_t tmpEta, tmpPhi;
 408       if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
 409       if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
 410         {
 411           maxPhi=tmpPhi;
 412           maxEta=tmpEta;
 413           index=ic;
 414         }
 415       }
 416   return index;
 417 }
 418
 419 //
 420 //------------------------------------------------------------------------------
 421 //
 422 void AliEMCALTracker::UnloadClusters()
 423 {
 424         //
 425         // Free memory from all arrays
 426         // This method is called after the local tracking step
 427         // so we can safely delete everything
 428         //
 429
 430         Clear();
 431 }
 432
 433 //
 434 //------------------------------------------------------------------------------
 435 //
 436 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
 437   : fIndex(index),
 438     fX(0.),
 439     fY(0.),
 440     fZ(0.)
 441 {
 442         //
 443         // Translates an AliEMCALRecPoint object into the internal format.
 444         // Index of passed cluster in its native array must be specified.
 445         //
 446         TVector3 clpos;
 447         recPoint->GetGlobalPosition(clpos);
 448
 449         fX = clpos.X();
 450         fY = clpos.Y();
 451         fZ = clpos.Z();
 452 }
 453 //
 454 //------------------------------------------------------------------------------
 455 //
 456 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
 457   : fIndex(index),
 458     fX(0.),
 459     fY(0.),
 460     fZ(0.)
 461 {
 462         //
 463         // Translates an AliESDCaloCluster object into the internal format.
 464         // Index of passed cluster in its native array must be specified.
 465         //
 466         Float_t clpos[3]= {0., 0., 0.};
 467         caloCluster->GetPosition(clpos);
 468
 469         fX = (Double_t)clpos[0];
 470         fY = (Double_t)clpos[1];
 471         fZ = (Double_t)clpos[2];
 472 }