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