[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEmEtReconstructed.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for MC analysis
//  - MC output
//  implementation file
//
//*-- Author: Marcelo G. Munhoz (USP)
	      //_________________________________________________________________________

#include "AliAnalysisEmEtReconstructed.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliVEvent.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "TH2F.h"
#include "TParticle.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "TList.h"
#include "AliESDCaloCluster.h"
#include "TGeoGlobalMagField.h"
#include "AliMagF.h"
#include "AliEMCALTrack.h"
#include "AliESDtrackCuts.h"
#include "AliEMCALGeometry.h"
#include "AliExternalTrackParam.h"
#include "AliTrackerBase.h"
#include "TGeoManager.h"

  using namespace std;

ClassImp(AliAnalysisEmEtReconstructed);


// ctor
AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
							    ,fAllRectotETDep(0)
							    ,fElectronMatchtotETDep(0)
							    ,fNeutralRectotET(0)
							    ,fTotEMRectotET(0)

							    ,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
							    ,fTotChargedMatchtotETDep(0)

							    ,fTotalRectotETDep(0)

							    ,fESD(0)
							    ,fGeoUt(0)

							    ,fHistAllRecETDep(0) 
							    ,fHistAllRec(0) 
							    ,fHistAllRectotETDep(0) 

							    ,fHistElectronRecETDep(0) 
							    ,fHistElectronRec(0) 
							    ,fHistElectronMatchtotETDep(0) 
							    ,fHistElectronRecdEdxP(0)

							    ,fHistNeutralRectotET(0)  

							    ,fHistTotEMRectotET(0)

							    ,fHistMuonRecETDep(0) 
							    ,fHistMuonRec(0) 
							    ,fHistMuonMatchtotETDep(0) 
							    ,fHistMuonRecdEdxP(0)

							    ,fHistPionRecETDep(0) 
							    ,fHistPionRec(0) 
							    ,fHistPionMatchtotETDep(0) 
							    ,fHistPionRecdEdxP(0)

							    ,fHistKaonRecETDep(0) 
							    ,fHistKaonRec(0) 
							    ,fHistKaonMatchtotETDep(0) 
							    ,fHistKaonRecdEdxP(0)

							    ,fHistProtonRecETDep(0) 
							    ,fHistProtonRec(0) 
							    ,fHistProtonMatchtotETDep(0) 
							    ,fHistProtonRecdEdxP(0)

							    ,fHistTotChargedMatchtotETDep(0)

							    ,fHistTotalRectotETDep(0)

							    ,fHistDeltaRZ(0)
{//constructor
  fHistogramNameSuffix = TString("EmcalRec");
	
  fResCut = 0.02;
  //fResCut = fEmcalTrackDistanceCut;
	
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
  //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  TGeoManager::Import("geometry.root");
}

// dtor
AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed() 
{//Destructor 
  delete fGeoUt;

  delete fHistAllRecETDep;
  delete fHistAllRec;
  delete fHistAllRectotETDep;
	
  delete fHistElectronRecETDep;
  delete fHistElectronRec;
  delete fHistElectronMatchtotETDep; 
	
  delete fHistElectronRecdEdxP;

  delete fHistNeutralRectotET;  

  delete fHistTotEMRectotET;

  delete fHistMuonRecETDep;
  delete fHistMuonRec;
  delete fHistMuonMatchtotETDep; 

  delete fHistMuonRecdEdxP;
	
  delete fHistPionRecETDep;
  delete fHistPionRec;
  delete fHistPionMatchtotETDep; 

  delete fHistPionRecdEdxP;

  delete fHistKaonRecETDep;
  delete fHistKaonRec;
  delete fHistKaonMatchtotETDep; 

  delete fHistKaonRecdEdxP;
	
  delete fHistProtonRecETDep;
  delete fHistProtonRec;
  delete fHistProtonMatchtotETDep; 

  delete fHistProtonRecdEdxP;
	
  delete fHistTotChargedMatchtotETDep;
	
  delete fHistTotalRectotETDep;
	
  //few checks
  delete fHistDeltaRZ;
	
}

Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
{ // analyse MC and real event info
  if(!ev){
    AliError("ERROR: Event does not exist");   
    return 0;
  }
	
  fESD = dynamic_cast<AliESDEvent*>(ev);

  if(!fGeoUt){
    fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
    //fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
    AliInfo("Creating new AliEMCALGeometry");
  }
  //fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
  if(!fGeoUt){
    AliInfo("No fGeoUt!");
  }
  else{
    if(!fESD->GetEMCALMatrix(0)){
      AliInfo("No matrix!");
    }
    else{
      fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
    }
  }
	
  ResetEventValues();
	
  // get all emcal clusters
  TRefArray* caloClusters = new TRefArray();
  fESD->GetEMCALClusters( caloClusters );
	
  Int_t nCluster = caloClusters->GetEntries();
	
  Float_t pos[3] = {0};
  TVector3 caloPos(0,0,0);
  TVector3 trackPos(0,0,0);
  Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
  Double_t xCluster[4]={0}, xCharged[7]={0};
	
  AliESDtrack *track = 0;
    
  // loop the clusters
  for (int iCluster = 0; iCluster < nCluster; iCluster++ ) 
    {		
      // Retrieve calo cluster information
      AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
      Float_t caloE = caloCluster->E();
      caloCluster->GetPosition(pos);		
      caloPos.SetXYZ(pos[0],pos[1],pos[2]);
		
      // look for track that matches calo cluster  
      //track = FindMatch(caloCluster, res); // Marcelo's matching	

      // *********************
      // tender's matching
      delta_eta = caloCluster->GetTrackDz(); 
      delta_phi = caloCluster->GetTrackDx(); 

      if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
          track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
        
      //if (track)
      //   if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
      //       track = 0;
      // *********************
        
      // Retrieve track PID
      if (track)
	maxPid = GetTrackPID(track);
      else
	maxPid = -99;
		
      // calculate ET
      Double_t etDep = CalculateTransverseEnergy(caloCluster);
		
      // All clusters
      //fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
      //fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
		
      xCluster[0] = caloE;
      xCluster[1] = caloPos.Eta();
      xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
      xCluster[3] = caloCluster->GetNCells();
      fAllRectotETDep += etDep;		
		
	
      if(fMakeSparse){
	fHistAllRecETDep->Fill(xCluster,etDep);
	fHistAllRec->Fill(xCluster);
      }

      if (track)
	{
	  xCharged[0] = track->Eta();
	  xCharged[1] = track->Pt();
	}
      else
	{
	  xCharged[0] = -99;
	  xCharged[1] = -99;
	}
      xCharged[2] = caloE;
      xCharged[3] = caloPos.Eta();
      xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
      xCharged[5] = caloCluster->GetNCells();
      xCharged[6] = res;
		
      Bool_t isCharged = kFALSE;
		
      if (maxPid == AliPID::kProton)
	{
			
	  if(fMakeSparse){
	    fHistProtonRecETDep->Fill(xCharged,etDep);
	    fHistProtonRec->Fill(xCharged);
	  }

	  fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fProtonMatchtotETDep += etDep;

	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kPion)
	{
			
	  if(fMakeSparse){
	    fHistPionRecETDep->Fill(xCharged,etDep);
	    fHistPionRec->Fill(xCharged);
	  }

	  fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fPionMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kKaon)
	{
			
	  if(fMakeSparse){
	    fHistKaonRecETDep->Fill(xCharged,etDep);
	    fHistKaonRec->Fill(xCharged);
	  }

	  fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {

	      fKaonMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kMuon)
	{
			
	  if(fMakeSparse){
	    fHistMuonRecETDep->Fill(xCharged,etDep);
	    fHistMuonRec->Fill(xCharged);
	  }
	  fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {

	      fMuonMatchtotETDep += etDep;						
	      isCharged = kTRUE;
	    }
	}
      else if (maxPid == AliPID::kElectron)
	{
			
	  if(fMakeSparse){
	    fHistElectronRecETDep->Fill(xCharged,etDep);
	    fHistElectronRec->Fill(xCharged);
	  }

	  fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
			
	  if ((res>0.) && (res<fResCut))
	    {
	      fElectronMatchtotETDep += etDep;
	      isCharged = kTRUE;
	    }
	}
		
      if (!isCharged)
	{
	  fNeutralRectotET += etDep;			
	}
		
    } // end of loop over clusters	
	
  fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
  fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
  fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
	
  fHistAllRectotETDep->Fill(fAllRectotETDep);
	
  fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep); 
  fHistNeutralRectotET->Fill(fNeutralRectotET);
	
  fHistTotEMRectotET->Fill(fTotEMRectotET);
	
  fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep); 
  fHistPionMatchtotETDep->Fill(fPionMatchtotETDep); 
  fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep); 
  fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep); 
  fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
	
  fHistTotalRectotETDep->Fill(fTotalRectotETDep);
	
  delete caloClusters;
	
  return 0;    
}

void AliAnalysisEmEtReconstructed::Init()
{ // init
  AliAnalysisEt::Init();
}


void AliAnalysisEmEtReconstructed::ResetEventValues()
{ // reset event values
  AliAnalysisEt::ResetEventValues();
	
  // collision geometry defaults for p+p:
  fAllRectotETDep = 0;
	
  fElectronMatchtotETDep = 0;
  fNeutralRectotET = 0;
	
  fTotEMRectotET = 0;
	
  fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
  fTotChargedMatchtotETDep = 0;
	
  fTotalRectotETDep = 0;
}


void AliAnalysisEmEtReconstructed::CreateHistograms()
{ // histogram related additions
  //AliAnalysisEt::CreateHistograms();
	
  if(fMakeSparse){
    fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
    fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
  }
  TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
  fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);


  if(fMakeSparse){
    fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
    fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
  }
  histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
  fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);

  histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
  fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);

  histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
  fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
	
  histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
  fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);

  if(fMakeSparse){
    fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
    fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
  }
  histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
  fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);

  histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
  fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
    fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
  }
  histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
  fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
  histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
  fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
    fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
  }
  histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
  fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);

  histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
  fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  if(fMakeSparse){
    fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
    fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
  }
  histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
  fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);

  histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
  fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	
  histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
  fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
	
  histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
  fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);		
	
  histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
  fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
}

void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
{//Function for filling the output list
  //AliAnalysisEt::FillOutputList(list);
	
	
  if(fMakeSparse){
    list->Add(fHistAllRecETDep); 
    list->Add(fHistAllRec); 
  }
  list->Add(fHistAllRectotETDep); 

  if(fMakeSparse){
    list->Add(fHistElectronRecETDep); 
    list->Add(fHistElectronRec); 
  }
  list->Add(fHistElectronMatchtotETDep); 
  list->Add(fHistElectronRecdEdxP);
	
	
  list->Add(fHistTotEMRectotET); 

  list->Add(fHistMuonRec); 
  list->Add(fHistMuonRecdEdxP);
  list->Add(fHistMuonMatchtotETDep); 

  if(fMakeSparse){
    list->Add(fHistPionRecETDep); 
    list->Add(fHistPionRec); 
  }
  list->Add(fHistPionMatchtotETDep); 
  list->Add(fHistPionRecdEdxP);
	
  if(fMakeSparse){
    list->Add(fHistKaonRecETDep); 
    list->Add(fHistKaonRec); 
  }
  list->Add(fHistKaonMatchtotETDep); 
  list->Add(fHistKaonRecdEdxP);

  if(fMakeSparse){
    list->Add(fHistProtonRecETDep); 
    list->Add(fHistProtonRec); 
  }
  list->Add(fHistProtonMatchtotETDep); 
  list->Add(fHistProtonRecdEdxP);
	
  list->Add(fHistTotChargedMatchtotETDep); 
  list->Add(fHistTotalRectotETDep); 
	
  list->Add(fHistDeltaRZ);
}

//________________________________________________________________________
//project to a EMCal radius
Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliExternalTrackParam *trackParam, TVector3 &trackPos)
{//Gets the projection of the track
  Bool_t proj = kFALSE;
  Double_t emcalR = fGeoUt->GetEMCGeometry()->GetIPDistance();
	
  if (trackParam) //it is constructed from TParticle
    {
      Double_t trkPos[3] = {0};
		
      //Assume the track is a pion with mass 0.139GeV/c^2
      //Extrapolation step is 1cm
      if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, emcalR, 0.139, 1, kTRUE, 0.8) ) return proj;
		
      trackParam->GetXYZ(trkPos);
		
      trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
		
      proj = kTRUE;               
    }
	
  return proj;
}

//________________________________________________________________________
//project to a cluster position
Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliEMCALTrack* emcTrack, TVector3 &trackPos, TVector3 clusPos)
{//project to a cluster position
  Bool_t proj = kFALSE;
	
  if (emcTrack)
    {	
      Double_t trkPos[3] = {0};
		
      emcTrack->PropagateToGlobal(clusPos.X(),clusPos.Y(),clusPos.Z(),0.,0.);
      emcTrack->GetXYZ(trkPos);
		
      trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
		
      proj = kTRUE;
    }
	
  return proj;
}

//________________________________________________________________________	
AliESDtrack* AliAnalysisEmEtReconstructed::FindMatch(const AliESDCaloCluster *caloCluster, Double_t& resMin)
{//find a matched track
  Double_t res=0;
  resMin=999;
	
  TVector3 caloPos(0,0,0);
  Float_t pos[3] = {0};
  caloCluster->GetPosition(pos);		
  caloPos.SetXYZ(pos[0],pos[1],pos[2]);
	
  // loop over tracks
  TVector3 trackPos(0,0,0);
  TVector3 trackMatchPos(0,0,0);
  AliEMCALTrack *emcTrack = 0;	
  AliESDtrack *trackMatch = 0;	
	
  TObjArray* list = fEsdtrackCutsITSTPC->GetAcceptedTracks(fESD);;
  Int_t nGoodTracks = list->GetEntries();
	
  for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) 
    {
      AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
      if (!track)
        {
	  AliError(Form("ERROR: Could not get track %d", iTrack));
	  continue;
        }
		
      emcTrack = new AliEMCALTrack(*track);
		
      if (GetTrackProjection(emcTrack,trackPos,caloPos))
	{
	  res = sqrt(pow(trackPos.Phi()-caloPos.Phi(),2)+pow(trackPos.Eta()-caloPos.Eta(),2));
			
	  if (res < resMin)
	    {
	      resMin = res;
	      trackMatch = track;
	      trackMatchPos.SetXYZ(trackPos.X(),trackPos.Y(),trackPos.Z());
	    }
	}
		
      delete emcTrack;
    }		
	
  fHistDeltaRZ->Fill(trackMatchPos.Phi()-caloPos.Phi(),trackMatchPos.Eta()-caloPos.Eta());

  return trackMatch;
}

//________________________________________________________________________	
Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
{//Get the default track ID
  const Double_t *pidWeights = track->PID();
  Int_t maxpid = -1;
  Double_t maxpidweight = 0;
	
  if (pidWeights)
    {
      for (Int_t p =0; p < AliPID::kSPECIES; p++)
	{
	  if (pidWeights[p] > maxpidweight)
	    {
	      maxpidweight = pidWeights[p];
	      maxpid = p;
	    }
	}
    }
	
  return maxpid;
}
Commit	Line	Data
2e12941b	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Base class for MC analysis
	4	// - MC output
	5	// implementation file
	6	//
	7	//*-- Author: Marcelo G. Munhoz (USP)
f2adae2b	8	//_________________________________________________________________________
2e12941b	9
	10	#include "AliAnalysisEmEtReconstructed.h"
	11	#include "AliAnalysisEtCuts.h"
	12	#include "AliESDtrack.h"
	13	#include "AliStack.h"
	14	#include "AliVEvent.h"
	15	#include "AliMCEvent.h"
	16	#include "AliESDEvent.h"
	17	#include "TH2F.h"
	18	#include "TParticle.h"
	19	#include "AliGenHijingEventHeader.h"
	20	#include "AliGenPythiaEventHeader.h"
	21	#include "TList.h"
	22	#include "AliESDCaloCluster.h"
	23	#include "TGeoGlobalMagField.h"
	24	#include "AliMagF.h"
	25	#include "AliEMCALTrack.h"
	26	#include "AliESDtrackCuts.h"
	27	#include "AliEMCALGeometry.h"
	28	#include "AliExternalTrackParam.h"
	29	#include "AliTrackerBase.h"
	30	#include "TGeoManager.h"
	31
f2adae2b	32	using namespace std;
2e12941b	33
	34	ClassImp(AliAnalysisEmEtReconstructed);
	35
	36
	37	// ctor
	38	AliAnalysisEmEtReconstructed::AliAnalysisEmEtReconstructed():AliAnalysisEtReconstructed()
f2adae2b	39	,fAllRectotETDep(0)
	40	,fElectronMatchtotETDep(0)
	41	,fNeutralRectotET(0)
	42	,fTotEMRectotET(0)
	43
	44	,fMuonMatchtotETDep(0), fPionMatchtotETDep(0), fKaonMatchtotETDep(0), fProtonMatchtotETDep(0)
	45	,fTotChargedMatchtotETDep(0)
	46
	47	,fTotalRectotETDep(0)
	48
	49	,fESD(0)
	50	,fGeoUt(0)
	51
	52	,fHistAllRecETDep(0)
	53	,fHistAllRec(0)
	54	,fHistAllRectotETDep(0)
	55
	56	,fHistElectronRecETDep(0)
	57	,fHistElectronRec(0)
	58	,fHistElectronMatchtotETDep(0)
	59	,fHistElectronRecdEdxP(0)
	60
	61	,fHistNeutralRectotET(0)
	62
	63	,fHistTotEMRectotET(0)
	64
	65	,fHistMuonRecETDep(0)
	66	,fHistMuonRec(0)
	67	,fHistMuonMatchtotETDep(0)
	68	,fHistMuonRecdEdxP(0)
	69
	70	,fHistPionRecETDep(0)
	71	,fHistPionRec(0)
	72	,fHistPionMatchtotETDep(0)
	73	,fHistPionRecdEdxP(0)
	74
	75	,fHistKaonRecETDep(0)
	76	,fHistKaonRec(0)
	77	,fHistKaonMatchtotETDep(0)
	78	,fHistKaonRecdEdxP(0)
	79
	80	,fHistProtonRecETDep(0)
	81	,fHistProtonRec(0)
	82	,fHistProtonMatchtotETDep(0)
	83	,fHistProtonRecdEdxP(0)
	84
	85	,fHistTotChargedMatchtotETDep(0)
	86
	87	,fHistTotalRectotETDep(0)
	88
	89	,fHistDeltaRZ(0)
2e12941b	90	{//constructor
f2adae2b	91	fHistogramNameSuffix = TString("EmcalRec");
2e12941b	92
f2adae2b	93	fResCut = 0.02;
f2adae2b	94	//fResCut = fEmcalTrackDistanceCut;
2e12941b	95
f2adae2b	96	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG));
	97	//TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
	98	TGeoManager::Import("geometry.root");
2e12941b	99	}
	100
	101	// dtor
	102	AliAnalysisEmEtReconstructed::~AliAnalysisEmEtReconstructed()
f2adae2b	103	{//Destructor
a4d91e7c	104	delete fGeoUt;
a4d91e7c	105
f2adae2b	106	delete fHistAllRecETDep;
	107	delete fHistAllRec;
	108	delete fHistAllRectotETDep;
a4d91e7c	109
f2adae2b	110	delete fHistElectronRecETDep;
	111	delete fHistElectronRec;
	112	delete fHistElectronMatchtotETDep;
a4d91e7c	113
f2adae2b	114	delete fHistElectronRecdEdxP;
a4d91e7c	115
f2adae2b	116	delete fHistNeutralRectotET;
a4d91e7c	117
f2adae2b	118	delete fHistTotEMRectotET;
a4d91e7c	119
f2adae2b	120	delete fHistMuonRecETDep;
	121	delete fHistMuonRec;
	122	delete fHistMuonMatchtotETDep;
a4d91e7c	123
f2adae2b	124	delete fHistMuonRecdEdxP;
a4d91e7c	125
f2adae2b	126	delete fHistPionRecETDep;
	127	delete fHistPionRec;
	128	delete fHistPionMatchtotETDep;
a4d91e7c	129
f2adae2b	130	delete fHistPionRecdEdxP;
a4d91e7c	131
f2adae2b	132	delete fHistKaonRecETDep;
	133	delete fHistKaonRec;
	134	delete fHistKaonMatchtotETDep;
a4d91e7c	135
f2adae2b	136	delete fHistKaonRecdEdxP;
a4d91e7c	137
f2adae2b	138	delete fHistProtonRecETDep;
	139	delete fHistProtonRec;
	140	delete fHistProtonMatchtotETDep;
a4d91e7c	141
f2adae2b	142	delete fHistProtonRecdEdxP;
a4d91e7c	143
f2adae2b	144	delete fHistTotChargedMatchtotETDep;
a4d91e7c	145
f2adae2b	146	delete fHistTotalRectotETDep;
a4d91e7c	147
f2adae2b	148	//few checks
f2adae2b	149	delete fHistDeltaRZ;
a4d91e7c	150
2e12941b	151	}
	152
	153	Int_t AliAnalysisEmEtReconstructed::AnalyseEvent(AliVEvent* ev)
	154	{ // analyse MC and real event info
f2adae2b	155	if(!ev){
	156	AliError("ERROR: Event does not exist");
	157	return 0;
	158	}
	159
	160	fESD = dynamic_cast<AliESDEvent*>(ev);
	161
	162	if(!fGeoUt){
	163	fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");//new AliEMCALGeometry("EMCAL_FIRSTYEAR","EMCAL");
7687e5b3	164	//fGeoUt = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
f2adae2b	165	AliInfo("Creating new AliEMCALGeometry");
	166	}
	167	//fGeoUt = new AliEMCALGeometry("EMCAL_COMPLETE1","EMCAL");
	168	if(!fGeoUt){
	169	AliInfo("No fGeoUt!");
	170	}
	171	else{
	172	if(!fESD->GetEMCALMatrix(0)){
	173	AliInfo("No matrix!");
	174	}
	175	else{
	176	fGeoUt->SetMisalMatrix(fESD->GetEMCALMatrix(0),0);
	177	}
	178	}
2e12941b	179
f2adae2b	180	ResetEventValues();
2e12941b	181
f2adae2b	182	// get all emcal clusters
	183	TRefArray* caloClusters = new TRefArray();
	184	fESD->GetEMCALClusters( caloClusters );
2e12941b	185
f2adae2b	186	Int_t nCluster = caloClusters->GetEntries();
2e12941b	187
f2adae2b	188	Float_t pos[3] = {0};
	189	TVector3 caloPos(0,0,0);
	190	TVector3 trackPos(0,0,0);
7687e5b3	191	Double_t res=0, delta_eta=0, delta_phi=0, maxPid=-99;
f2adae2b	192	Double_t xCluster[4]={0}, xCharged[7]={0};
7687e5b3	193
	194	AliESDtrack *track = 0;
	195
f2adae2b	196	// loop the clusters
	197	for (int iCluster = 0; iCluster < nCluster; iCluster++ )
	198	{
	199	// Retrieve calo cluster information
	200	AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
	201	Float_t caloE = caloCluster->E();
	202	caloCluster->GetPosition(pos);
	203	caloPos.SetXYZ(pos[0],pos[1],pos[2]);
2e12941b	204
f2adae2b	205	// look for track that matches calo cluster
7687e5b3	206	//track = FindMatch(caloCluster, res); // Marcelo's matching
	207
	208	// *********************
	209	// tender's matching
	210	delta_eta = caloCluster->GetTrackDz();
	211	delta_phi = caloCluster->GetTrackDx();
	212
	213	if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
	214	track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
	215
	216	//if (track)
	217	// if ( !fEsdtrackCutsITSTPC->IsSelected(track) )
	218	// track = 0;
	219	// *********************
	220
f2adae2b	221	// Retrieve track PID
	222	if (track)
	223	maxPid = GetTrackPID(track);
	224	else
	225	maxPid = -99;
2e12941b	226
f2adae2b	227	// calculate ET
f2adae2b	228	Double_t etDep = CalculateTransverseEnergy(caloCluster);
2e12941b	229
f2adae2b	230	// All clusters
	231	//fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
	232	//fHistAllRecEtaETDep->Fill(etDep,caloPos.Eta());
2e12941b	233
f2adae2b	234	xCluster[0] = caloE;
	235	xCluster[1] = caloPos.Eta();
	236	xCluster[2] = TMath::RadToDeg()*caloPos.Phi();
	237	xCluster[3] = caloCluster->GetNCells();
	238	fAllRectotETDep += etDep;
2e12941b	239
f2adae2b	240
	241	if(fMakeSparse){
	242	fHistAllRecETDep->Fill(xCluster,etDep);
	243	fHistAllRec->Fill(xCluster);
	244	}
	245
	246	if (track)
	247	{
	248	xCharged[0] = track->Eta();
	249	xCharged[1] = track->Pt();
	250	}
	251	else
	252	{
	253	xCharged[0] = -99;
	254	xCharged[1] = -99;
	255	}
	256	xCharged[2] = caloE;
	257	xCharged[3] = caloPos.Eta();
	258	xCharged[4] = TMath::RadToDeg()*caloPos.Phi();
	259	xCharged[5] = caloCluster->GetNCells();
	260	xCharged[6] = res;
2e12941b	261
f2adae2b	262	Bool_t isCharged = kFALSE;
2e12941b	263
f2adae2b	264	if (maxPid == AliPID::kProton)
f2adae2b	265	{
2e12941b	266
f2adae2b	267	if(fMakeSparse){
	268	fHistProtonRecETDep->Fill(xCharged,etDep);
	269	fHistProtonRec->Fill(xCharged);
	270	}
2e12941b	271
f2adae2b	272	fHistProtonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	273
f2adae2b	274	if ((res>0.) && (res<fResCut))
	275	{
	276	fProtonMatchtotETDep += etDep;
	277
	278	isCharged = kTRUE;
	279	}
	280	}
	281	else if (maxPid == AliPID::kPion)
	282	{
2e12941b	283
f2adae2b	284	if(fMakeSparse){
	285	fHistPionRecETDep->Fill(xCharged,etDep);
	286	fHistPionRec->Fill(xCharged);
	287	}
2e12941b	288
f2adae2b	289	fHistPionRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	290
f2adae2b	291	if ((res>0.) && (res<fResCut))
	292	{
	293	fPionMatchtotETDep += etDep;
	294	isCharged = kTRUE;
	295	}
	296	}
	297	else if (maxPid == AliPID::kKaon)
	298	{
2e12941b	299
f2adae2b	300	if(fMakeSparse){
	301	fHistKaonRecETDep->Fill(xCharged,etDep);
	302	fHistKaonRec->Fill(xCharged);
	303	}
2e12941b	304
f2adae2b	305	fHistKaonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	306
f2adae2b	307	if ((res>0.) && (res<fResCut))
f2adae2b	308	{
2e12941b	309
f2adae2b	310	fKaonMatchtotETDep += etDep;
	311	isCharged = kTRUE;
	312	}
	313	}
	314	else if (maxPid == AliPID::kMuon)
	315	{
2e12941b	316
f2adae2b	317	if(fMakeSparse){
	318	fHistMuonRecETDep->Fill(xCharged,etDep);
	319	fHistMuonRec->Fill(xCharged);
	320	}
	321	fHistMuonRecdEdxP->Fill(track->P(),track->GetTPCsignal());
2e12941b	322
f2adae2b	323	if ((res>0.) && (res<fResCut))
f2adae2b	324	{
2e12941b	325
f2adae2b	326	fMuonMatchtotETDep += etDep;
	327	isCharged = kTRUE;
	328	}
	329	}
	330	else if (maxPid == AliPID::kElectron)
	331	{
2e12941b	332
f2adae2b	333	if(fMakeSparse){
	334	fHistElectronRecETDep->Fill(xCharged,etDep);
	335	fHistElectronRec->Fill(xCharged);
	336	}
	337
	338	fHistElectronRecdEdxP->Fill(track->P(),track->GetTPCsignal());
	339
	340	if ((res>0.) && (res<fResCut))
	341	{
	342	fElectronMatchtotETDep += etDep;
	343	isCharged = kTRUE;
	344	}
	345	}
2e12941b	346
f2adae2b	347	if (!isCharged)
	348	{
	349	fNeutralRectotET += etDep;
	350	}
2e12941b	351
f2adae2b	352	} // end of loop over clusters
2e12941b	353
f2adae2b	354	fTotEMRectotET = fElectronMatchtotETDep + fNeutralRectotET;
	355	fTotChargedMatchtotETDep = fMuonMatchtotETDep + fPionMatchtotETDep + fKaonMatchtotETDep + fProtonMatchtotETDep;
	356	fTotalRectotETDep = fTotEMRectotET + fTotChargedMatchtotETDep;
2e12941b	357
f2adae2b	358	fHistAllRectotETDep->Fill(fAllRectotETDep);
2e12941b	359
f2adae2b	360	fHistElectronMatchtotETDep->Fill(fElectronMatchtotETDep);
f2adae2b	361	fHistNeutralRectotET->Fill(fNeutralRectotET);
2e12941b	362
f2adae2b	363	fHistTotEMRectotET->Fill(fTotEMRectotET);
2e12941b	364
f2adae2b	365	fHistMuonMatchtotETDep->Fill(fMuonMatchtotETDep);
	366	fHistPionMatchtotETDep->Fill(fPionMatchtotETDep);
	367	fHistKaonMatchtotETDep->Fill(fKaonMatchtotETDep);
	368	fHistProtonMatchtotETDep->Fill(fProtonMatchtotETDep);
	369	fHistTotChargedMatchtotETDep->Fill(fTotChargedMatchtotETDep);
2e12941b	370
f2adae2b	371	fHistTotalRectotETDep->Fill(fTotalRectotETDep);
2e12941b	372
f2adae2b	373	delete caloClusters;
2e12941b	374
f2adae2b	375	return 0;
2e12941b	376	}
	377
	378	void AliAnalysisEmEtReconstructed::Init()
	379	{ // init
f2adae2b	380	AliAnalysisEt::Init();
2e12941b	381	}
	382
	383
	384	void AliAnalysisEmEtReconstructed::ResetEventValues()
	385	{ // reset event values
f2adae2b	386	AliAnalysisEt::ResetEventValues();
2e12941b	387
f2adae2b	388	// collision geometry defaults for p+p:
f2adae2b	389	fAllRectotETDep = 0;
2e12941b	390
f2adae2b	391	fElectronMatchtotETDep = 0;
f2adae2b	392	fNeutralRectotET = 0;
2e12941b	393
f2adae2b	394	fTotEMRectotET = 0;
2e12941b	395
f2adae2b	396	fMuonMatchtotETDep = 0; fPionMatchtotETDep = 0; fKaonMatchtotETDep = 0; fProtonMatchtotETDep = 0;
f2adae2b	397	fTotChargedMatchtotETDep = 0;
2e12941b	398
f2adae2b	399	fTotalRectotETDep = 0;
2e12941b	400	}
	401
	402
	403	void AliAnalysisEmEtReconstructed::CreateHistograms()
	404	{ // histogram related additions
f2adae2b	405	//AliAnalysisEt::CreateHistograms();
	406
	407	if(fMakeSparse){
	408	fHistAllRecETDep = CreateClusterHistoSparse("fHistAllRecETDep_","E_{T}, all particles");
	409	fHistAllRec = CreateClusterHistoSparse("fHistAllRec_","counts, all particles");
	410	}
	411	TString histname = "fHistAllRectotETDep_" + fHistogramNameSuffix;
	412	fHistAllRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
	413
	414
	415	if(fMakeSparse){
	416	fHistElectronRecETDep = CreateChargedPartHistoSparse("fHistElectronRecETDep_","E_{T}, electrons");
	417	fHistElectronRec = CreateChargedPartHistoSparse("fHistElectronRec_","counts, electrons");
	418	}
	419	histname = "fHistElectronMatchtotETDep_" + fHistogramNameSuffix;
	420	fHistElectronMatchtotETDep = new TH1F(histname.Data(),"total ET, MC primary Electrons",fgNumOfEBins, fgEAxis);
	421
	422	histname = "fHistElectronRecdEdxP_" + fHistogramNameSuffix;
	423	fHistElectronRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	424
	425	histname = "fHistNeutralRectotET_" + fHistogramNameSuffix;
	426	fHistNeutralRectotET = new TH1F(histname.Data(),"total ET, neutral particles",fgNumOfEBins, fgEAxis);
	427
	428	histname = "fHistTotEMRectotET_" + fHistogramNameSuffix;
	429	fHistTotEMRectotET = new TH1F(histname.Data(),"total electromagnetic ET",fgNumOfEBins, fgEAxis);
	430
	431	if(fMakeSparse){
	432	fHistMuonRecETDep = CreateChargedPartHistoSparse("fHistMuonRecETDep_","E_{T}, muons");
	433	fHistMuonRec = CreateChargedPartHistoSparse("fHistMuonRec_","counts, muons");
	434	}
	435	histname = "fHistMuonMatchtotETDep_" + fHistogramNameSuffix;
	436	fHistMuonMatchtotETDep = new TH1F(histname.Data(),"total ET, Muons",fgNumOfEBins, fgEAxis);
	437
	438	histname = "fHistMuonRecdEdxP_" + fHistogramNameSuffix;
	439	fHistMuonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	440
	441	if(fMakeSparse){
	442	fHistPionRecETDep = CreateChargedPartHistoSparse("fHistPionRecETDep_","E_{T}, pions");
	443	fHistPionRec = CreateChargedPartHistoSparse("fHistPionRec_","counts, pions");
	444	}
	445	histname = "fHistPionMatchtotETDep_" + fHistogramNameSuffix;
	446	fHistPionMatchtotETDep = new TH1F(histname.Data(),"total ET, Pions",fgNumOfEBins, fgEAxis);
	447	histname = "fHistPionRecdEdxP_" + fHistogramNameSuffix;
	448	fHistPionRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	449
	450	if(fMakeSparse){
	451	fHistKaonRecETDep = CreateChargedPartHistoSparse("fHistKaonRecETDep_","E_{T}, kaons");
	452	fHistKaonRec = CreateChargedPartHistoSparse("fHistKaonRec_","counts, kaons");
	453	}
	454	histname = "fHistKaonMatchtotETDep_" + fHistogramNameSuffix;
	455	fHistKaonMatchtotETDep = new TH1F(histname.Data(),"total ET, Kaons",fgNumOfEBins, fgEAxis);
	456
	457	histname = "fHistKaonRecdEdxP_" + fHistogramNameSuffix;
	458	fHistKaonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
	459
	460	if(fMakeSparse){
	461	fHistProtonRecETDep = CreateChargedPartHistoSparse("fHistProtonRecETDep_","E_{T}, protons");
	462	fHistProtonRec = CreateChargedPartHistoSparse("fHistProtonRec_","counts, protons");
	463	}
	464	histname = "fHistProtonMatchtotETDep_" + fHistogramNameSuffix;
	465	fHistProtonMatchtotETDep = new TH1F(histname.Data(),"total ET, Protons",fgNumOfEBins, fgEAxis);
	466
	467	histname = "fHistProtonRecdEdxP_" + fHistogramNameSuffix;
	468	fHistProtonRecdEdxP = new TH2F(histname,"TPC dEdx vs P",100,0.,10.,100,0.,200.);
469
470	histname = "fHistTotChargedMatchtotETDep_" + fHistogramNameSuffix;
471	fHistTotChargedMatchtotETDep = new TH1F(histname.Data(),"total ET, charged particles",fgNumOfEBins, fgEAxis);
472
473	histname = "fHistTotalRectotETDep_" + fHistogramNameSuffix;
474	fHistTotalRectotETDep = new TH1F(histname.Data(),"total ET, all particles",fgNumOfEBins, fgEAxis);
475
476	histname = "fHistDeltaRZ_" + fHistogramNameSuffix;
477	fHistDeltaRZ = new TH2F(histname,"#Delta#phi vs #Delta#eta (track projection - cluster position)",200,-0.1,0.1,200,-0.1,0.1);
2e12941b	478	}
	479
	480	void AliAnalysisEmEtReconstructed::FillOutputList(TList *list)
	481	{//Function for filling the output list
f2adae2b	482	//AliAnalysisEt::FillOutputList(list);
	483
	484
	485	if(fMakeSparse){
	486	list->Add(fHistAllRecETDep);
	487	list->Add(fHistAllRec);
	488	}
	489	list->Add(fHistAllRectotETDep);
	490
	491	if(fMakeSparse){
	492	list->Add(fHistElectronRecETDep);
	493	list->Add(fHistElectronRec);
	494	}
	495	list->Add(fHistElectronMatchtotETDep);
	496	list->Add(fHistElectronRecdEdxP);
	497
	498
	499	list->Add(fHistTotEMRectotET);
	500
	501	list->Add(fHistMuonRec);
	502	list->Add(fHistMuonRecdEdxP);
	503	list->Add(fHistMuonMatchtotETDep);
	504
	505	if(fMakeSparse){
	506	list->Add(fHistPionRecETDep);
	507	list->Add(fHistPionRec);
	508	}
	509	list->Add(fHistPionMatchtotETDep);
	510	list->Add(fHistPionRecdEdxP);
	511
	512	if(fMakeSparse){
	513	list->Add(fHistKaonRecETDep);
	514	list->Add(fHistKaonRec);
	515	}
	516	list->Add(fHistKaonMatchtotETDep);
	517	list->Add(fHistKaonRecdEdxP);
	518
	519	if(fMakeSparse){
	520	list->Add(fHistProtonRecETDep);
	521	list->Add(fHistProtonRec);
	522	}
	523	list->Add(fHistProtonMatchtotETDep);
	524	list->Add(fHistProtonRecdEdxP);
	525
	526	list->Add(fHistTotChargedMatchtotETDep);
	527	list->Add(fHistTotalRectotETDep);
	528
	529	list->Add(fHistDeltaRZ);
2e12941b	530	}
	531
	532	//________________________________________________________________________
f2adae2b	533	//project to a EMCal radius
2e12941b	534	Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliExternalTrackParam *trackParam, TVector3 &trackPos)
2e12941b	535	{//Gets the projection of the track
f2adae2b	536	Bool_t proj = kFALSE;
f2adae2b	537	Double_t emcalR = fGeoUt->GetEMCGeometry()->GetIPDistance();
2e12941b	538
f2adae2b	539	if (trackParam) //it is constructed from TParticle
2e12941b	540	{
f2adae2b	541	Double_t trkPos[3] = {0};
2e12941b	542
f2adae2b	543	//Assume the track is a pion with mass 0.139GeV/c^2
	544	//Extrapolation step is 1cm
	545	if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, emcalR, 0.139, 1, kTRUE, 0.8) ) return proj;
2e12941b	546
f2adae2b	547	trackParam->GetXYZ(trkPos);
2e12941b	548
f2adae2b	549	trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
2e12941b	550
f2adae2b	551	proj = kTRUE;
2e12941b	552	}
2e12941b	553
f2adae2b	554	return proj;
2e12941b	555	}
	556
	557	//________________________________________________________________________
f2adae2b	558	//project to a cluster position
2e12941b	559	Bool_t AliAnalysisEmEtReconstructed::GetTrackProjection(AliEMCALTrack* emcTrack, TVector3 &trackPos, TVector3 clusPos)
2e12941b	560	{//project to a cluster position
f2adae2b	561	Bool_t proj = kFALSE;
2e12941b	562
f2adae2b	563	if (emcTrack)
	564	{
	565	Double_t trkPos[3] = {0};
2e12941b	566
f2adae2b	567	emcTrack->PropagateToGlobal(clusPos.X(),clusPos.Y(),clusPos.Z(),0.,0.);
f2adae2b	568	emcTrack->GetXYZ(trkPos);
2e12941b	569
f2adae2b	570	trackPos.SetXYZ(trkPos[0],trkPos[1],trkPos[2]);
2e12941b	571
f2adae2b	572	proj = kTRUE;
f2adae2b	573	}
2e12941b	574
f2adae2b	575	return proj;
2e12941b	576	}
	577
	578	//________________________________________________________________________
	579	AliESDtrack* AliAnalysisEmEtReconstructed::FindMatch(const AliESDCaloCluster *caloCluster, Double_t& resMin)
	580	{//find a matched track
f2adae2b	581	Double_t res=0;
f2adae2b	582	resMin=999;
2e12941b	583
f2adae2b	584	TVector3 caloPos(0,0,0);
	585	Float_t pos[3] = {0};
	586	caloCluster->GetPosition(pos);
	587	caloPos.SetXYZ(pos[0],pos[1],pos[2]);
2e12941b	588
f2adae2b	589	// loop over tracks
	590	TVector3 trackPos(0,0,0);
	591	TVector3 trackMatchPos(0,0,0);
	592	AliEMCALTrack *emcTrack = 0;
	593	AliESDtrack *trackMatch = 0;
2e12941b	594
f2adae2b	595	TObjArray* list = fEsdtrackCutsITSTPC->GetAcceptedTracks(fESD);;
f2adae2b	596	Int_t nGoodTracks = list->GetEntries();
2e12941b	597
f2adae2b	598	for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
	599	{
	600	AliESDtrack track = dynamic_cast<AliESDtrack> (list->At(iTrack));
	601	if (!track)
2e12941b	602	{
f2adae2b	603	AliError(Form("ERROR: Could not get track %d", iTrack));
f2adae2b	604	continue;
2e12941b	605	}
2e12941b	606
f2adae2b	607	emcTrack = new AliEMCALTrack(*track);
2e12941b	608
f2adae2b	609	if (GetTrackProjection(emcTrack,trackPos,caloPos))
	610	{
	611	res = sqrt(pow(trackPos.Phi()-caloPos.Phi(),2)+pow(trackPos.Eta()-caloPos.Eta(),2));
2e12941b	612
f2adae2b	613	if (res < resMin)
	614	{
	615	resMin = res;
	616	trackMatch = track;
	617	trackMatchPos.SetXYZ(trackPos.X(),trackPos.Y(),trackPos.Z());
	618	}
	619	}
2e12941b	620
f2adae2b	621	delete emcTrack;
f2adae2b	622	}
2e12941b	623
f2adae2b	624	fHistDeltaRZ->Fill(trackMatchPos.Phi()-caloPos.Phi(),trackMatchPos.Eta()-caloPos.Eta());
2e12941b	625
f2adae2b	626	return trackMatch;
2e12941b	627	}
	628
	629	//________________________________________________________________________
	630	Double_t AliAnalysisEmEtReconstructed::GetTrackPID(const AliESDtrack *track) const
	631	{//Get the default track ID
f2adae2b	632	const Double_t *pidWeights = track->PID();
	633	Int_t maxpid = -1;
	634	Double_t maxpidweight = 0;
2e12941b	635
f2adae2b	636	if (pidWeights)
	637	{
	638	for (Int_t p =0; p < AliPID::kSPECIES; p++)
2e12941b	639	{
f2adae2b	640	if (pidWeights[p] > maxpidweight)
	641	{
	642	maxpidweight = pidWeights[p];
	643	maxpid = p;
	644	}
2e12941b	645	}
f2adae2b	646	}
2e12941b	647
f2adae2b	648	return maxpid;
2e12941b	649	}