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