[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 "AliEMCALTracker.h"

ClassImp(AliEMCALTracker)

//
//------------------------------------------------------------------------------
//
AliEMCALTracker::AliEMCALTracker() 
  : AliTracker(),
    fCutPt(0),
    fCutNITS(0),
    fCutNTPC(50),
    fStep(50),
    fTrackCorrMode(kTrackCorrMMB),	
    fCutEta(0.025),
    fCutPhi(0.05),
    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),
    fCutEta(copy.fCutEta),
    fCutPhi(copy.fCutPhi),
    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& copy)
{
  //
  // Assignment operator.
  // Besides copying all parameters, duplicates all collections.	
  //

  fCutPt  = copy.fCutPt;
  fCutEta = copy.fCutEta;
  fCutPhi = copy.fCutPhi;	
  fStep = copy.fStep;
  fTrackCorrMode = copy.fTrackCorrMode;

  fCutNITS = copy.fCutNITS;
  fCutNTPC = copy.fCutNTPC;
  
  fTracks = (TObjArray*)copy.fTracks->Clone();
  fClusters = (TObjArray*)copy.fClusters->Clone();
  fGeom = copy.fGeom;
  
  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.
  //
	
  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);

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

      //Reject kink daughters
      isKink = kFALSE;
      for (Int_t j = 0; j < 3; j++)
	{
	  if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
	}
      if (isKink) continue;

      //Loose geometric cut
      Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
      if(TMath::Abs(esdTrack->Eta())>1 || phi <= 60 || phi >= 200 ) 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
  //

  Double_t maxEta=fCutEta;
  Double_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;
  const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
  if(friendTrack && friendTrack->GetTPCOut())
    trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
  else
    trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
  if(!trkParam) return index;

  //Perform extrapolation
  Double_t trkPos[3];
  Int_t nclusters = fClusters->GetEntries();
  for(Int_t ic=0; ic<nclusters; ic++)
    {
      //AliExternalTrackParam *trkParamTmp = new AliExternalTrackParam(*trkParam);
      AliExternalTrackParam trkParamTmp(*trkParam);
      AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
      TVector3 vec(cluster->X(),cluster->Y(),cluster->Z());
      Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
      //Rotate to proper local coordinate
      vec.RotateZ(-alpha); 
      trkParamTmp.Rotate(alpha); 
      //extrapolation is done here
      if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), fStep, kFALSE, 0.8, -1)) continue; 

      //Calculate the residuals
      trkParamTmp.GetXYZ(trkPos);        
      TVector3 clsPosVec(cluster->X(),cluster->Y(),cluster->Z());
      TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
      
      Double_t clsPhi = clsPosVec.Phi();
      if(clsPhi<0) clsPhi+=2*TMath::Pi();
      Double_t trkPhi = trkPosVec.Phi();
      if(trkPhi<0) trkPhi+=2*TMath::Pi();

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