[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),
  fEMCalSurfaceDistance(440),
  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),
    fEMCalSurfaceDistance(copy.fEMCalSurfaceDistance),
    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, pt;
  if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))  {
	if(fITSTrackSA) delete trkParam;
	return index;
  }
  track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
  if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
	 if(fITSTrackSA) delete trkParam;
	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;
        }
      }

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