[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.cxx

//========================================================================
// 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.                  
//======================================================================== 
//                       
//                       Class AliEMCALTracker 
//                      -----------------------
// Implementation of the track matching method between barrel tracks and
// EMCAL clusters.
// Besides algorithm implementation, some cuts are required to be set
// in order to define, for each track, an acceptance window where clusters
// are searched to find best match (if any).
// The class accepts as input an ESD container, and works directly on it,
// simply setting, for each of its tracks, the fEMCALindex flag, for each
// track which is matched to a cluster.
// In order to use method, one must launch PropagateBack().
//
// ------------------------------------------------------------------------
// author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
// Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch)
//=========================================================================

#include <Riostream.h>
#include <iomanip>

#include <TFile.h>
#include <TTree.h>
#include <TList.h>
#include <TString.h>
#include <TVector3.h>
#include <TClonesArray.h>
#include <TGeoMatrix.h>

#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDCaloCluster.h"
#include "AliEMCALRecPoint.h"
#include "AliRunLoader.h"
#include "AliEMCALTrack.h"
#include "AliEMCALLoader.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALReconstructor.h"
#include "AliEMCALRecParam.h"
#include "AliCDBEntry.h"
#include "AliCDBManager.h"
#include "AliEMCALReconstructor.h"
#include "AliEMCALRecoUtils.h"

#include "AliEMCALTracker.h"

using std::cerr;
using std::endl;
ClassImp(AliEMCALTracker)

//
//------------------------------------------------------------------------------
//
AliEMCALTracker::AliEMCALTracker() 
: AliTracker(),
  fCutPt(0),
  fCutNITS(0),
  fCutNTPC(50),
  fStep(20),
  fTrackCorrMode(kTrackCorrMMB),
  fClusterWindow(50),
  fCutEta(0.025),
  fCutPhi(0.05),
  fITSTrackSA(kFALSE),
  fTracks(0),
  fClusters(0),
  fGeom(0)
{
  //
  // Default constructor.
  // Initializes all simple data members to default values,
   // and all collections to NULL.
  // Output file name is set to a default value.
  //
  InitParameters();
}
//
//------------------------------------------------------------------------------
//
AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) 
  : AliTracker(),
    fCutPt(copy.fCutPt),
    fCutNITS(copy.fCutNITS),
    fCutNTPC(copy.fCutNTPC),
    fStep(copy.fStep),
    fTrackCorrMode(copy.fTrackCorrMode),
    fClusterWindow(copy.fClusterWindow),
    fCutEta(copy.fCutEta),
    fCutPhi(copy.fCutPhi),
    fITSTrackSA(copy.fITSTrackSA), 
    fTracks((TObjArray*)copy.fTracks->Clone()),
    fClusters((TObjArray*)copy.fClusters->Clone()),
    fGeom(copy.fGeom)
{
  //
  // Copy constructor
  // Besides copying all parameters, duplicates all collections.
  //
}
//
//------------------------------------------------------------------------------
//
AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
{ // assignment operator; use copy ctor
  if (&source == this) return *this;

  new (this) AliEMCALTracker(source);
  return *this;
}
//
//------------------------------------------------------------------------------
//
void AliEMCALTracker::InitParameters()
{
  //
  // Retrieve initialization parameters
  //
	
  // Check if the instance of AliEMCALRecParam exists, 
  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();

  if(!recParam){
    AliFatal("Reconstruction parameters for EMCAL not set!");
  }
  else{
 
    fCutEta  =  recParam->GetMthCutEta();
    fCutPhi  =  recParam->GetMthCutPhi();
    fStep    =  recParam->GetExtrapolateStep();
    fCutPt   =  recParam->GetTrkCutPt();
    fCutNITS =  recParam->GetTrkCutNITS();
    fCutNTPC =  recParam->GetTrkCutNTPC();
  }
	
}

//
//------------------------------------------------------------------------------
//
void AliEMCALTracker::Clear(Option_t* option)
{
	//
	// Clearing method
        // Deletes all objects in arrays and the arrays themselves
	//

	TString opt(option);
	Bool_t clearTracks = opt.Contains("TRACKS");
	Bool_t clearClusters = opt.Contains("CLUSTERS");
	if (opt.Contains("ALL")) {
		clearTracks = kTRUE;
		clearClusters = kTRUE;
	}
	
	//fTracks is a collection of esdTrack
	//When clearing this array, the linked objects should not be deleted
	if (fTracks != 0x0 && clearTracks) {
	   fTracks->Clear();
	   delete fTracks;
           fTracks = 0;
	}
	if (fClusters != 0x0 && clearClusters) {
	   fClusters->Delete();
	   delete fClusters;
	   fClusters = 0;
	}
}
//
//------------------------------------------------------------------------------
//
Int_t AliEMCALTracker::LoadClusters(TTree *cTree) 
{
	//
	// Load EMCAL clusters in the form of AliEMCALRecPoint,
	// from simulation temporary files.
	// (When included in reconstruction chain, this method is used automatically)
	//
	
	Clear("CLUSTERS");

	cTree->SetBranchStatus("*",0); //disable all branches
	cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need

	TBranch *branch = cTree->GetBranch("EMCALECARP");
	if (!branch) {
		AliError("Can't get the branch with the EMCAL clusters");
		return 1;
	}
	
	TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
	branch->SetAddress(&clusters);
	
	//cTree->GetEvent(0);
	branch->GetEntry(0);
	Int_t nClusters = (Int_t)clusters->GetEntries();
	if(fClusters) fClusters->Delete();
	else fClusters = new TObjArray(0);
	for (Int_t i = 0; i < nClusters; i++) {
		AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
		if (!cluster) continue;
		AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
		fClusters->AddLast(matchCluster);
	}

	branch->SetAddress(0);
        clusters->Delete();
        delete clusters;

	AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));

	return 0;
}
//
//------------------------------------------------------------------------------
//
Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd) 
{
  //
  // Load EMCAL clusters in the form of AliESDCaloClusters,
  // from an AliESD object.
  //
  
  // make sure that tracks/clusters collections are empty
  Clear("CLUSTERS");
  fClusters = new TObjArray(0);
  
  Int_t nClusters = esd->GetNumberOfCaloClusters();       		
  for (Int_t i=0; i<nClusters; i++) 
    {
      AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
      if (!cluster || !cluster->IsEMCAL()) continue ; 
      AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
      fClusters->AddLast(matchCluster);
    }
  
  AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
  return 0;
}
//
//------------------------------------------------------------------------------
//
Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
{
  //
  // Load ESD tracks.
  //

  UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
  UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
  Bool_t desc1 = (mask1 >> 3) & 0x1;
  Bool_t desc2 = (mask2 >> 3) & 0x1;
  if (desc1==0 || desc2==0) {
    AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
                  mask1, esd->GetESDRun()->GetDetectorsInReco(),
                  mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
    fITSTrackSA = kTRUE;
  }
  
  Clear("TRACKS");
  fTracks = new TObjArray(0);
	
  Int_t nTracks = esd->GetNumberOfTracks();
  //Bool_t isKink=kFALSE;
  for (Int_t i = 0; i < nTracks; i++) 
    {
      AliESDtrack *esdTrack = esd->GetTrack(i);
      // set by default the value corresponding to "no match"
      esdTrack->SetEMCALcluster(kUnmatched);
      esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);

      //Select good quaulity tracks
      if(esdTrack->Pt()<fCutPt) continue;
      if(!fITSTrackSA)
	if(esdTrack->GetNcls(1)<fCutNTPC)continue;

      //Loose geometric cut
      Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
      if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;

      fTracks->AddLast(esdTrack);
    }

      AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
      return 0;
}
//
//------------------------------------------------------------------------------
//
void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
{
  //
  // Set track correction mode
  // gest the choice in string format and converts into 
  // internal enum
  //
  
  TString opt(option);
  opt.ToUpper();
  
  if (!opt.CompareTo("NONE")) 
    {
      fTrackCorrMode = kTrackCorrNone;
    }
  else if (!opt.CompareTo("MMB")) 
    {
      fTrackCorrMode = kTrackCorrMMB;
    }
  else 
    {
      cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
    }
}
//
//------------------------------------------------------------------------------
//
Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
{
	//
	// Main operation method.
	// Gets external AliESD containing tracks to be matched.
	// After executing match finding, stores in the same ESD object all infos
	// and releases the object for further reconstruction steps.
	//
        //
        // Note: should always return 0=OK, because otherwise all tracking
        // is aborted for this event
  
	if (!esd) {
		AliError("NULL ESD passed");
		return 1;
	}
	
	// step 1: collect clusters
	Int_t okLoadClusters, nClusters;
	if (!fClusters || (fClusters && fClusters->IsEmpty())) {
		okLoadClusters = LoadClusters(esd);
	}
	nClusters = fClusters->GetEntries();
		
	// step 2: collect ESD tracks
	Int_t nTracks, okLoadTracks;
	okLoadTracks = LoadTracks(esd);
	nTracks = fTracks->GetEntries();
	
	// step 3: for each track, find the closest cluster as matched within residual cuts
	Int_t index=-1;
	for (Int_t it = 0; it < nTracks; it++) 
	  {
	    AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
	    index = FindMatchedCluster(track);
	    if (index>-1) 
	      {
		AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
		track->SetEMCALcluster(cluster->Index());
		track->SetStatus(AliESDtrack::kEMCALmatch);
	      }
	  }

	return 0;
}

//
//------------------------------------------------------------------------------
//
Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
{         
  //
  // For each track, extrapolate it to all the clusters
  // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
  //

  Float_t maxEta=fCutEta;
  Float_t maxPhi=fCutPhi;
  Int_t index = -1;
  
  // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
  // Otherwise use the TPCInner point
  AliExternalTrackParam *trkParam = 0;
  
  if(!fITSTrackSA){ 
    const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
    if(friendTrack && friendTrack->GetTPCOut())
      trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
    else if(track->GetInnerParam())
      trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
  }
  else
    trkParam = new AliExternalTrackParam(*track);
  
  if(!trkParam) return index;
  
  AliExternalTrackParam trkParamTmp(*trkParam);
  Float_t eta, phi;
  if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(kTRUE), fStep, eta, phi))  return index;
  track->SetTrackPhiEtaOnEMCal(phi,eta);
  if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;

  //Perform extrapolation
  Double_t trkPos[3];
  trkParamTmp.GetXYZ(trkPos);
  Int_t nclusters = fClusters->GetEntries();
  for(Int_t ic=0; ic<nclusters; ic++)
    {
      AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
      Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
      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));
//       printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
      if(dR > fClusterWindow) continue;
      
      AliExternalTrackParam trkParTmp(trkParamTmp);

      Float_t tmpEta, tmpPhi;
      if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
      if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
        {
          maxPhi=tmpPhi;
          maxEta=tmpEta;
          index=ic;
        }
      }
  return index;
}

//
//------------------------------------------------------------------------------
//
void AliEMCALTracker::UnloadClusters() 
{
	//
	// Free memory from all arrays
        // This method is called after the local tracking step
        // so we can safely delete everything 
	//
	
  	Clear();
}

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