AliEMCALRecoUtilsBase.cxx

Go to the documentation of this file.

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

// ROOT includes
#include <TGeoManager.h>
#include <TGeoBBox.h>

// STEER includes
#include "AliVCluster.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliExternalTrackParam.h"
#include "AliTrackerBase.h"

// EMCAL includes
#include "AliEMCALRecoUtilsBase.h"
#include "AliEMCALGeometry.h"

ClassImp(AliEMCALRecoUtilsBase) ;

//_____________________________________
AliEMCALRecoUtilsBase::AliEMCALRecoUtilsBase()
{
}

//
// Copy constructor.
//
//______________________________________________________________________
AliEMCALRecoUtilsBase::AliEMCALRecoUtilsBase(const AliEMCALRecoUtilsBase & reco) 
: TNamed(reco)
{  
}

//______________________________________________________________________
AliEMCALRecoUtilsBase & AliEMCALRecoUtilsBase::operator = (const AliEMCALRecoUtilsBase & reco) 
{  
  if (this == &reco)return *this;
  ((TNamed *)this)->operator=(reco);
  
  return *this;
}

//------------------------------------------------------------------------------------
Bool_t AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
                                                             Double_t emcalR, 
                                                             Double_t mass,
                                                             Double_t step, 
                                                             Double_t minpt,
                                                             Bool_t useMassForTracking, 
                                                             Bool_t useDCA)
{ 
  track->SetTrackPhiEtaPtOnEMCal(-999, -999, -999);
  
  if ( track->Pt() < minpt )
    return kFALSE;
  
  if ( TMath::Abs(track->Eta()) > 0.9 ) return kFALSE;
  
  // Save some time and memory in case of no DCal present
  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
  
  Int_t       nSupMod = AliEMCALGeoParams::fgkEMCALModules;
  if ( geom ) nSupMod = geom->GetNumberOfSuperModules();
  
  if ( nSupMod < 13 ) // Run1 10 (12, 2 not active but present)
  {
    Double_t phi = track->Phi()*TMath::RadToDeg();
    if ( phi <= 10 || phi >= 250 ) return kFALSE;
  }
  
  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(track);
  AliAODTrack *aodt = 0;
  if (!esdt) 
  {
    aodt = dynamic_cast<AliAODTrack*>(track);
    
    if (!aodt)
      return kFALSE;
  }
  
  // Select the mass hypothesis
  if ( mass < 0 )
  {
    Bool_t onlyTPC = kFALSE;
    if ( mass == -99 ) onlyTPC=kTRUE;
    
    if (esdt)
    {
      if ( useMassForTracking ) mass = esdt->GetMassForTracking();
      else                      mass = esdt->GetMass(onlyTPC);
    }
    else
    {
      if ( useMassForTracking ) mass = aodt->GetMassForTracking();
      else                      mass = aodt->M();
    }
  }
  
  AliExternalTrackParam *trackParam = 0;
  if (esdt) 
  {
    const AliExternalTrackParam *in = esdt->GetInnerParam();
    
    if (!in)
      return kFALSE;
    
    trackParam = new AliExternalTrackParam(*in);
  } 
  else 
  {
    Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0};
    aodt->PxPyPz(pxpypz);
    if (useDCA) {
      // Note: useDCA is by default false in this function only for backwards compatibility
      aodt->GetXYZ(xyz);
    }
    else {
      aodt->XvYvZv(xyz);
    }
    aodt->GetCovarianceXYZPxPyPz(cv);  
    trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge());
  }
  
  if (!trackParam)
    return kFALSE;
  
  Float_t etaout=-999, phiout=-999, ptout=-999;
  Bool_t ret = ExtrapolateTrackToEMCalSurface(trackParam, 
                                              emcalR,
                                              mass,
                                              step,
                                              etaout, 
                                              phiout,
                                              ptout);
  
  delete trackParam;
  
  if (!ret) return kFALSE;
  
  if ( TMath::Abs(etaout) > 0.75 ) return kFALSE;
  
  // Save some time and memory in case of no DCal present
  if ( nSupMod < 13 ) // Run1 10 (12, 2 not active but present)
  {
    if ( (phiout < 70*TMath::DegToRad()) || (phiout > 190*TMath::DegToRad()) )  return kFALSE;
  }
  
  track->SetTrackPhiEtaPtOnEMCal(phiout, etaout, ptout);
  
  return kTRUE;
}


//---------------------------------------------------------------------------------------
Bool_t AliEMCALRecoUtilsBase::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, 
                                                             Double_t emcalR,
                                                             Double_t mass, 
                                                             Double_t step, 
                                                             Float_t &eta, 
                                                             Float_t &phi,
                                                             Float_t &pt)
{
  eta = -999, phi = -999, pt = -999;
  
  if (!trkParam) return kFALSE;
  
  if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
  
  Double_t trkPos[3] = {0.,0.,0.};
  
  if (!trkParam->GetXYZ(trkPos)) return kFALSE;
  
  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
  
  eta = trkPosVec.Eta();
  phi = trkPosVec.Phi();
  pt  = trkParam->Pt();
  
  if ( phi < 0 )
    phi += TMath::TwoPi();
  
  return kTRUE;
}

//-----------------------------------------------------------------------------------
Bool_t AliEMCALRecoUtilsBase::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, 
                                                         const Float_t *clsPos, 
                                                         Double_t mass, 
                                                         Double_t step, 
                                                         Float_t &tmpEta, 
                                                         Float_t &tmpPhi)
{
  tmpEta = -999;
  tmpPhi = -999;
  
  if (!trkParam) return kFALSE;
  
  Double_t trkPos[3] = {0.,0.,0.};
  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
  
  Float_t phi = vec.Phi();
  
  if ( phi < 0 )
    phi += TMath::TwoPi();
  
  // Rotate the cluster to the local extrapolation coordinate system
  Double_t alpha = ((int)(phi*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
  
  vec.RotateZ(-alpha); 
  
  if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) 
    return kFALSE;
  
  if (!trkParam->GetXYZ(trkPos)) 
    return kFALSE; // Get the extrapolated global position
  
  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
  
  // Track cluster matching
  tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
  tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
  
  return kTRUE;
}

//----------------------------------------------------------------------------------
Bool_t AliEMCALRecoUtilsBase::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, 
                                                        const AliVCluster *cluster, 
                                                        Double_t mass, 
                                                        Double_t step, 
                                                        Float_t &tmpEta, 
                                                        Float_t &tmpPhi)
{
  tmpEta = -999;
  tmpPhi = -999;
  
  if (!cluster || !trkParam) 
    return kFALSE;
  
  Float_t clsPos[3] = {0.,0.,0.};
  cluster->GetPosition(clsPos);
  
  return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
}

//_________________________________________________________
Float_t  AliEMCALRecoUtilsBase::GetDepth(Float_t energy, 
                                         Int_t   iParticle, 
                                         Int_t   iSM ) const 
{
  // parameters 
  Float_t x0    = 1.31;
  Float_t ecr   = 8;
  Float_t depth = 0;
  Float_t arg   = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
  
  switch ( iParticle )
  {
    case kPhoton:
      if (arg < 1) 
        depth = 0;
      else
        depth = x0 * (TMath::Log(arg) + 0.5); 
      break;
      
    case kElectron:
      if (arg < 1) 
        depth = 0;
      else
        depth = x0 * (TMath::Log(arg) - 0.5); 
      break;
      
    case kHadron:
      // hadron 
      // boxes anc. here
      if (gGeoManager) 
      {
        gGeoManager->cd("ALIC_1/XEN1_1");
        TGeoNode        *geoXEn1    = gGeoManager->GetCurrentNode();
        TGeoNodeMatrix  *geoSM      = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
        if (geoSM) 
        {
          TGeoVolume      *geoSMVol   = geoSM->GetVolume(); 
          TGeoShape       *geoSMShape = geoSMVol->GetShape();
          TGeoBBox        *geoBox     = dynamic_cast<TGeoBBox *>(geoSMShape);
          if (geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
          else AliFatal("Null GEANT box");
        }
        else AliFatal("NULL  GEANT node matrix");
      }
      else
      {//electron
        if (arg < 1) 
          depth = 0;
        else
          depth = x0 * (TMath::Log(arg) - 0.5); 
      }
      
      break;
      
    default://photon
      if (arg < 1) 
        depth = 0;
      else
        depth = x0 * (TMath::Log(arg) + 0.5);
  }  
  
  return depth;
}