[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackESDReader.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.                  *
 **************************************************************************/
/* $Id:  $ */

//_________________________________________________________________________
// Class for reading data (ESDs) in order to do prompt gamma 
// or other particle identification and correlations
//
// It is a filtering class, transforms ESD tracks or CaloClusters
// into AOD tracks and calocluters, which are the basic input of the analysis
// classes in this frame.
// It is recommended to use the official filter AliAnalysisTaskESDfilter, and 
// then the reader for AODs AliCaloTrackAODReader.
//
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
//#include "Riostream.h"

//---- ANALYSIS system ----
#include "AliCaloTrackESDReader.h" 
#include "AliESDEvent.h"
#include "AliESDCaloCluster.h"
#include "AliAODCaloCluster.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "AliAODEvent.h"
#include "AliFidutialCut.h"

ClassImp(AliCaloTrackESDReader)

//____________________________________________________________________________
AliCaloTrackESDReader::AliCaloTrackESDReader() : 
AliCaloTrackReader()
{
	//Default Ctor
	
	//Initialize parameters
	fDataType=kESD;
	fReadStack          = kTRUE;
	fReadAODMCParticles = kFALSE;
	//We want tracks fitted in the detectors:
	fTrackStatus=AliESDtrack::kTPCrefit;
    fTrackStatus|=AliESDtrack::kITSrefit;

}

//____________________________________________________________________________
AliCaloTrackESDReader::AliCaloTrackESDReader(const AliCaloTrackESDReader & g) :   
  AliCaloTrackReader(g)
{
  // cpy ctor
}

//_________________________________________________________________________
//AliCaloTrackESDReader & AliCaloTrackESDReader::operator = (const AliCaloTrackESDReader & source)
//{
//	// assignment operator
//	
//	if(&source == this) return *this;
//	
//	return *this;
//	
//}

//____________________________________________________________________________
void AliCaloTrackESDReader::FillInputCTS() {
  //Return array with CTS tracks
  //TObjArray * fAODCTS = new TObjArray();
  Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  Int_t naod = 0;
  Double_t pos[3];
  Double_t p[3];
  Double_t covTr[21];
  Double_t pid[10];
  Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();

  Double_t      timezero        = 0;   //TO BE FIXED

  //To be replaced by call to AliEMCALGeoUtils when the class becomes available
  Double_t radius = 441.0; //[cm] EMCAL radius +13cm
  if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - org entries %d\n", nTracks);
  for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
    AliESDtrack * track = (AliESDtrack*) ((AliESDEvent*)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
    
    //We want tracks fitted in the detectors: TPCrefit, ITSrefit ... check the set bits.
	if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
      
	track->GetPxPyPz(p) ;
	TLorentzVector momentum(p[0],p[1],p[2],0);
      
	if(fCTSPtMin < momentum.Pt() &&fFidutialCut->IsInFidutialCut(momentum,"CTS")){
	
		if(fDebug > 3 && momentum.Pt() > 0.2) printf("AliCaloTrackESDReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	
		track->GetXYZ(pos);
		track->GetCovarianceXYZPxPyPz(covTr);
		track->GetESDpid(pid);
	
		Float_t impactXY, impactZ;
	
		track->GetImpactParameters(impactXY,impactZ);
	
		if (impactXY<3) {
		  // track inside the beam pipe
		  //Put new aod object in file in AOD tracks array
		  AliAODTrack *aodTrack = new((*(fOutputEvent->GetTracks()))[naod++]) 
		    AliAODTrack(track->GetID(), track->GetLabel(), p, kTRUE, pos, kFALSE,covTr, (Short_t)track->GetSign(), track->GetITSClusterMap(), 
				pid,
				0x0,//primary,
				kTRUE, // check if this is right
				kTRUE, // check if this is right
				AliAODTrack::kPrimary, 
				0);
		  
		  aodTrack->SetFlags(track->GetStatus());
		  aodTrack->ConvertAliPIDtoAODPID();

		  //fill needed AliAODPid information, including
		  //Extrapolate track to EMCAL surface for AOD-level track-cluster matching
		  AliAODPid *aodpid = new AliAODPid;
		  aodpid->SetITSsignal(track->GetITSsignal());
		  aodpid->SetTPCsignal(track->GetTPCsignal());
		  //n TRD planes = 6
		  Int_t nslices = track->GetNumberOfTRDslices()*6;
		  Double_t *trdslices = new Double_t[nslices];
		  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
		    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
		  }
		  aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
		  Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
		  aodpid->SetIntegratedTimes(times);

		  aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed 
		  aodpid->SetHMPIDsignal(track->GetHMPIDsignal());

		  Double_t emcpos[3] = {0.,0.,0.};
		  Double_t emcmom[3] = {0.,0.,0.};
		  aodpid->SetEMCALPosition(emcpos);
		  aodpid->SetEMCALMomentum(emcmom);
		  
		  AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();
		  if(!outerparam) continue;
		  
		  Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
		  Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
		  if(!(okpos && okmom)) continue;
		  
		  aodpid->SetEMCALPosition(emcpos);
		  aodpid->SetEMCALMomentum(emcmom);
		  
		  aodTrack->SetDetPID(aodpid);
		}
		else continue;   // outside the beam pipe: orphan track	
	}//Pt and Fidutial cut passed. 
  }// track loop

  //Put references to selected tracks in array
  for(Int_t itrack = 0; itrack < (fOutputEvent->GetTracks())->GetEntriesFast(); itrack++){
    AliAODTrack * track =  (AliAODTrack*) (fOutputEvent->GetTracks())->At(itrack);	
    fAODCTS->Add(track);				
  }	
  
  if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - aod entries %d\n", fAODCTS->GetEntriesFast());
}

//____________________________________________________________________________
void AliCaloTrackESDReader::FillInputEMCAL() {
  //Return array with EMCAL clusters in aod format
  
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  GetVertex(v);
  
  Float_t pos[3] ;
  Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
  Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
  
  //Loop to select clusters in fidutial cut and fill container with aodClusters
  for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
    AliESDCaloCluster * clus = 0;
    if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
      if (clus->IsEMCAL()){
	TLorentzVector momentum ;
	clus->GetMomentum(momentum, v);      
	if(fDebug > 3 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						   momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); 
	if(fEMCALPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"EMCAL")){
	  
	  if(fDebug > 2 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						     momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	  
	  clus->GetPosition(pos) ;
	  Int_t id = clus->GetID();
	  Int_t nLabel = clus->GetNLabels();
	  Int_t *labels=0x0;
	  if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
	  
	  Float_t energy = clus->E();
	  Char_t ttype= AliAODCluster::kEMCALClusterv1;
	  
	  //Put new aod object in file in AOD calo clusters array
	  AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++]) 
	    AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL,ttype,0);
	  
	  //       printf("Reader PID ESD: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f \n",
	  // 	     pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	  // 	     pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
	  caloCluster->SetPIDFromESD(clus->GetPid());
	  caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(), 
				      clus->GetM20(), clus->GetM02(),   
				      clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF()) ;
	  
	  
	  if(fDebug > 3 && momentum.E() > 0.2)
	    printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, NexMax %d, TOF %e\n",
		   clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
		   clus->GetNExMax(), clus->GetTOF());
	  
	  caloCluster->SetNCells(clus->GetNCells());
	  caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
	  caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
	  
	  TArrayI* matchedT = 	clus->GetTracksMatched();
	  if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {	
	    for (Int_t im = 0; im < matchedT->GetSize(); im++) {
	      Int_t iESDtrack = matchedT->At(im);
	      if(iESDtrack < nTracks && iESDtrack > -1)
		caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
	    }
	  }
	  //Fill reference array
	}//Pt and Fidutial cut passed.
      }//EMCAL cluster
    }//cluster exists
  }//esd cluster loop
  
  //Put references to selected clusters in array
  for(Int_t iclus = 0; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
    AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);	
    fAODEMCAL->Add(clus);				
  }
  if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
  
}

//____________________________________________________________________________
void AliCaloTrackESDReader::FillInputPHOS() {
  //Return array with PHOS clusters in aod format
  Int_t nEMCAL = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  GetVertex(v);
  
  Float_t pos[3] ;
  Double_t * pid = new Double_t[AliPID::kSPECIESN];
  Int_t naod      = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
  Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;

  //Loop to select clusters in fidutial cut and fill container with aodClusters
  for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
    AliESDCaloCluster * clus = 0;
    if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
      if (clus->IsPHOS()){
	TLorentzVector momentum ;
	clus->GetMomentum(momentum, v);      
	if(fDebug > 3 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						   momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	if(fPHOSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"PHOS")){
	  
	  if(fDebug > 2 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
						     momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	  
	  clus->GetPosition(pos) ;
	  Int_t id = clus->GetID();
	  Int_t nLabel = clus->GetNLabels();
	  Int_t *labels=0x0;
	  if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
	  Float_t energy = clus->E();
	  
	  //Phos cluster type
	  pid = clus->GetPid();
	  // printf("Reader PID ESD: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f \n",
	  // 	 pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	  // 	 pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
	  Char_t ttype= AliAODCluster::kPHOSNeutral;
	  Float_t wNeutral = pid[AliPID::kNeutron]+ pid[AliPID::kKaon0]+pid[AliPID::kPhoton]+pid[AliPID::kPi0];
	  Float_t wCharged = pid[AliPID::kMuon] + pid[AliPID::kElectron] + pid[AliPID::kEleCon]+ 
	    pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion];
	  if( wCharged > wNeutral)  ttype= AliAODCluster::kPHOSCharged;
	  
	  //Put new aod object in file in AOD calo clusters array
	  AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++]) 
	    AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0);


	  caloCluster->SetPIDFromESD(pid);   
	  caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(), 
				      clus->GetM20(), clus->GetM02(),  
				      clus->GetEmcCpvDistance(),  clus->GetNExMax()) ;
	  
	  if(fDebug > 3 && momentum.E() > 0.2)
	    printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, EmcCpvDist %3.3f, NexMax %d, TOF %e\n",
		   clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
		   clus->GetEmcCpvDistance(),  clus->GetNExMax(), clus->GetTOF());
	  
	  caloCluster->SetNCells(clus->GetNCells());
	  caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
	  caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());

	  TArrayI* matchedT = 	clus->GetTracksMatched();
	  if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {	
	    for (Int_t im = 0; im < matchedT->GetSize(); im++) {
	      Int_t iESDtrack = matchedT->At(im);
	      if(iESDtrack < nTracks && iESDtrack > -1)
		caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
	    }
	  }
	}//Pt and Fidutial cut passed.
      }//PHOS cluster
    }//cluster exists
  }//esd cluster loop
  
  //Put references to selected clusters in array
  for(Int_t iclus = nEMCAL; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
    AliAODCaloCluster * clus =  (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);	
    fAODPHOS->Add(clus);				
  }	
  if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
  
}

//____________________________________________________________________________
void AliCaloTrackESDReader::FillInputEMCALCells() {
  //Return array with EMCAL cells in esd format
  
  fEMCALCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetEMCALCells(); 
  
}

//____________________________________________________________________________
void AliCaloTrackESDReader::FillInputPHOSCells() {
  //Return array with PHOS cells in esd format
  
  fPHOSCells = (TNamed*) ((AliESDEvent*)fInputEvent)->GetPHOSCells(); 
  
}

//____________________________________________________________________________
void AliCaloTrackESDReader::GetVertex(Double_t  v[3]) const {
  //Return vertex position
  
  //((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
  ((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
	
}


//____________________________________________________________________________
Double_t AliCaloTrackESDReader::GetBField() const {
  //Return magnetic field

  Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();

  return bfield;
}


//____________________________________________________________________________
void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
  // Connect the data pointers
  
  if(strcmp(esd->GetName(),"AliESDEvent")){
    printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());
    abort();
  }
  
  SetInputEvent(esd);
  SetOutputEvent(aod);
  SetMC(mc);
  
}
Commit	Line	Data
1c5acb87	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16	/* $Id: $ */
	17
	18	//_________________________________________________________________________
	19	// Class for reading data (ESDs) in order to do prompt gamma
	20	// or other particle identification and correlations
	21	//
8dacfd76	22	// It is a filtering class, transforms ESD tracks or CaloClusters
	23	// into AOD tracks and calocluters, which are the basic input of the analysis
	24	// classes in this frame.
	25	// It is recommended to use the official filter AliAnalysisTaskESDfilter, and
	26	// then the reader for AODs AliCaloTrackAODReader.
591cc579	27	//
591cc579	28	//
1c5acb87	29	//*-- Author: Gustavo Conesa (LNF-INFN)
	30	//////////////////////////////////////////////////////////////////////////////
	31
	32
	33	// --- ROOT system ---
	34	//#include "Riostream.h"
	35
	36	//---- ANALYSIS system ----
	37	#include "AliCaloTrackESDReader.h"
	38	#include "AliESDEvent.h"
1c5acb87	39	#include "AliESDCaloCluster.h"
	40	#include "AliAODCaloCluster.h"
	41	#include "AliAODTrack.h"
bf1d1697	42	#include "AliAODPid.h"
1c5acb87	43	#include "AliAODEvent.h"
477d6cee	44	#include "AliFidutialCut.h"
1c5acb87	45
	46	ClassImp(AliCaloTrackESDReader)
	47
	48	//____________________________________________________________________________
	49	AliCaloTrackESDReader::AliCaloTrackESDReader() :
	50	AliCaloTrackReader()
	51	{
	52	//Default Ctor
	53
	54	//Initialize parameters
	55	fDataType=kESD;
591cc579	56	fReadStack = kTRUE;
	57	fReadAODMCParticles = kFALSE;
	58	//We want tracks fitted in the detectors:
	59	fTrackStatus=AliESDtrack::kTPCrefit;
	60	fTrackStatus\|=AliESDtrack::kITSrefit;
	61
1c5acb87	62	}
	63
	64	//____________________________________________________________________________
	65	AliCaloTrackESDReader::AliCaloTrackESDReader(const AliCaloTrackESDReader & g) :
477d6cee	66	AliCaloTrackReader(g)
1c5acb87	67	{
477d6cee	68	// cpy ctor
1c5acb87	69	}
	70
	71	//_________________________________________________________________________
	72	//AliCaloTrackESDReader & AliCaloTrackESDReader::operator = (const AliCaloTrackESDReader & source)
	73	//{
	74	// // assignment operator
	75	//
	76	// if(&source == this) return *this;
	77	//
	78	// return *this;
	79	//
	80	//}
	81
	82	//____________________________________________________________________________
	83	void AliCaloTrackESDReader::FillInputCTS() {
477d6cee	84	//Return array with CTS tracks
591cc579	85	//TObjArray * fAODCTS = new TObjArray();
477d6cee	86	Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
	87	Int_t naod = 0;
	88	Double_t pos[3];
	89	Double_t p[3];
	90	Double_t covTr[21];
	91	Double_t pid[10];
591cc579	92	Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
8dacfd76	93
bf1d1697	94	Double_t timezero = 0; //TO BE FIXED
591cc579	95
	96	//To be replaced by call to AliEMCALGeoUtils when the class becomes available
	97	Double_t radius = 441.0; //[cm] EMCAL radius +13cm
b551439d	98	if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - org entries %d\n", nTracks);
477d6cee	99	for (Int_t itrack = 0; itrack < nTracks; itrack++) {////////////// track loop
	100	AliESDtrack * track = (AliESDtrack) ((AliESDEvent)fInputEvent)->GetTrack(itrack) ; // retrieve track from esd
	101
591cc579	102	//We want tracks fitted in the detectors: TPCrefit, ITSrefit ... check the set bits.
591cc579	103	if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
477d6cee	104
591cc579	105	track->GetPxPyPz(p) ;
591cc579	106	TLorentzVector momentum(p[0],p[1],p[2],0);
477d6cee	107
591cc579	108	if(fCTSPtMin < momentum.Pt() &&fFidutialCut->IsInFidutialCut(momentum,"CTS")){
1c5acb87	109
591cc579	110	if(fDebug > 3 && momentum.Pt() > 0.2) printf("AliCaloTrackESDReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	111	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
1c5acb87	112
591cc579	113	track->GetXYZ(pos);
	114	track->GetCovarianceXYZPxPyPz(covTr);
	115	track->GetESDpid(pid);
477d6cee	116
591cc579	117	Float_t impactXY, impactZ;
477d6cee	118
591cc579	119	track->GetImpactParameters(impactXY,impactZ);
477d6cee	120
591cc579	121	if (impactXY<3) {
bf1d1697	122	// track inside the beam pipe
	123	//Put new aod object in file in AOD tracks array
	124	AliAODTrack aodTrack = new(((fOutputEvent->GetTracks()))[naod++])
	125	AliAODTrack(track->GetID(), track->GetLabel(), p, kTRUE, pos, kFALSE,covTr, (Short_t)track->GetSign(), track->GetITSClusterMap(),
	126	pid,
	127	0x0,//primary,
	128	kTRUE, // check if this is right
	129	kTRUE, // check if this is right
	130	AliAODTrack::kPrimary,
	131	0);
	132
	133	aodTrack->SetFlags(track->GetStatus());
	134	aodTrack->ConvertAliPIDtoAODPID();
	135
	136	//fill needed AliAODPid information, including
	137	//Extrapolate track to EMCAL surface for AOD-level track-cluster matching
	138	AliAODPid *aodpid = new AliAODPid;
	139	aodpid->SetITSsignal(track->GetITSsignal());
	140	aodpid->SetTPCsignal(track->GetTPCsignal());
	141	//n TRD planes = 6
	142	Int_t nslices = track->GetNumberOfTRDslices()*6;
	143	Double_t *trdslices = new Double_t[nslices];
	144	for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
	145	for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
	146	}
	147	aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
	148	Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
	149	aodpid->SetIntegratedTimes(times);
	150
	151	aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed
	152	aodpid->SetHMPIDsignal(track->GetHMPIDsignal());
	153
	154	Double_t emcpos[3] = {0.,0.,0.};
	155	Double_t emcmom[3] = {0.,0.,0.};
	156	aodpid->SetEMCALPosition(emcpos);
	157	aodpid->SetEMCALMomentum(emcmom);
	158
	159	AliExternalTrackParam outerparam = (AliExternalTrackParam)track->GetOuterParam();
b551439d	160	if(!outerparam) continue;
bf1d1697	161
	162	Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
	163	Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
b551439d	164	if(!(okpos && okmom)) continue;
bf1d1697	165
	166	aodpid->SetEMCALPosition(emcpos);
	167	aodpid->SetEMCALMomentum(emcmom);
	168
	169	aodTrack->SetDetPID(aodpid);
591cc579	170	}
591cc579	171	else continue; // outside the beam pipe: orphan track
bf1d1697	172	}//Pt and Fidutial cut passed.
477d6cee	173	}// track loop
b551439d	174
477d6cee	175	//Put references to selected tracks in array
	176	for(Int_t itrack = 0; itrack < (fOutputEvent->GetTracks())->GetEntriesFast(); itrack++){
	177	AliAODTrack * track = (AliAODTrack*) (fOutputEvent->GetTracks())->At(itrack);
	178	fAODCTS->Add(track);
	179	}
	180
a3aebfff	181	if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputCTS() - aod entries %d\n", fAODCTS->GetEntriesFast());
1c5acb87	182	}
	183
	184	//____________________________________________________________________________
	185	void AliCaloTrackESDReader::FillInputEMCAL() {
477d6cee	186	//Return array with EMCAL clusters in aod format
	187
	188	//Get vertex for momentum calculation
	189	Double_t v[3] ; //vertex ;
	190	GetVertex(v);
	191
477d6cee	192	Float_t pos[3] ;
a3aebfff	193	Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
	194	Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
	195
	196	//Loop to select clusters in fidutial cut and fill container with aodClusters
477d6cee	197	for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
	198	AliESDCaloCluster * clus = 0;
	199	if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
	200	if (clus->IsEMCAL()){
	201	TLorentzVector momentum ;
	202	clus->GetMomentum(momentum, v);
a3aebfff	203	if(fDebug > 3 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	204	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	205	if(fEMCALPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"EMCAL")){
	206
a3aebfff	207	if(fDebug > 2 && momentum.E() > 0.2) printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	208	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
a3aebfff	209
477d6cee	210	clus->GetPosition(pos) ;
477d6cee	211	Int_t id = clus->GetID();
	212	Int_t nLabel = clus->GetNLabels();
	213	Int_t *labels=0x0;
	214	if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray();
	215
	216	Float_t energy = clus->E();
	217	Char_t ttype= AliAODCluster::kEMCALClusterv1;
	218
	219	//Put new aod object in file in AOD calo clusters array
	220	AliAODCaloCluster caloCluster = new(((fOutputEvent->GetCaloClusters()))[naod++])
	221	AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL,ttype,0);
	222
a3aebfff	223	// printf("Reader PID ESD: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f \n",
	224	// pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	225	// pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
477d6cee	226	caloCluster->SetPIDFromESD(clus->GetPid());
	227	caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
	228	clus->GetM20(), clus->GetM02(),
	229	clus->GetEmcCpvDistance(), clus->GetNExMax(), clus->GetTOF()) ;
	230
	231
	232	if(fDebug > 3 && momentum.E() > 0.2)
a3aebfff	233	printf("AliCaloTrackESDReader::FillInputEMCAL() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, NexMax %d, TOF %e\n",
477d6cee	234	clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
	235	clus->GetNExMax(), clus->GetTOF());
	236
	237	caloCluster->SetNCells(clus->GetNCells());
	238	caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
	239	caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
	240
	241	TArrayI* matchedT = clus->GetTracksMatched();
a3aebfff	242	if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
477d6cee	243	for (Int_t im = 0; im < matchedT->GetSize(); im++) {
a3aebfff	244	Int_t iESDtrack = matchedT->At(im);
	245	if(iESDtrack < nTracks && iESDtrack > -1)
	246	caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
477d6cee	247	}
	248	}
	249	//Fill reference array
	250	}//Pt and Fidutial cut passed.
	251	}//EMCAL cluster
	252	}//cluster exists
	253	}//esd cluster loop
	254
	255	//Put references to selected clusters in array
	256	for(Int_t iclus = 0; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
	257	AliAODCaloCluster * clus = (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
	258	fAODEMCAL->Add(clus);
	259	}
a3aebfff	260	if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());
477d6cee	261
1c5acb87	262	}
	263
	264	//____________________________________________________________________________
	265	void AliCaloTrackESDReader::FillInputPHOS() {
477d6cee	266	//Return array with PHOS clusters in aod format
8896b363	267	Int_t nEMCAL = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
477d6cee	268	//Get vertex for momentum calculation
	269	Double_t v[3] ; //vertex ;
	270	GetVertex(v);
	271
477d6cee	272	Float_t pos[3] ;
477d6cee	273	Double_t * pid = new Double_t[AliPID::kSPECIESN];
a3aebfff	274	Int_t naod = (fOutputEvent->GetCaloClusters())->GetEntriesFast();
	275	Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
	276
	277	//Loop to select clusters in fidutial cut and fill container with aodClusters
477d6cee	278	for (Int_t iclus = 0; iclus < ((AliESDEvent*)fInputEvent)->GetNumberOfCaloClusters(); iclus++) {
	279	AliESDCaloCluster * clus = 0;
	280	if ( (clus = ((AliESDEvent*)fInputEvent)->GetCaloCluster(iclus)) ) {
	281	if (clus->IsPHOS()){
	282	TLorentzVector momentum ;
	283	clus->GetMomentum(momentum, v);
a3aebfff	284	if(fDebug > 3 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - all clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	285	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
	286	if(fPHOSPtMin < momentum.Pt() && fFidutialCut->IsInFidutialCut(momentum,"PHOS")){
	287
a3aebfff	288	if(fDebug > 2 && momentum.E() > 0.1)printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
477d6cee	289	momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
477d6cee	290
477d6cee	291	clus->GetPosition(pos) ;
	292	Int_t id = clus->GetID();
	293	Int_t nLabel = clus->GetNLabels();
	294	Int_t *labels=0x0;
	295	if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray();
	296	Float_t energy = clus->E();
	297
	298	//Phos cluster type
a3aebfff	299	pid = clus->GetPid();
	300	// printf("Reader PID ESD: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f \n",
	301	// pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	302	// pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]);
477d6cee	303	Char_t ttype= AliAODCluster::kPHOSNeutral;
	304	Float_t wNeutral = pid[AliPID::kNeutron]+ pid[AliPID::kKaon0]+pid[AliPID::kPhoton]+pid[AliPID::kPi0];
	305	Float_t wCharged = pid[AliPID::kMuon] + pid[AliPID::kElectron] + pid[AliPID::kEleCon]+
	306	pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion];
	307	if( wCharged > wNeutral) ttype= AliAODCluster::kPHOSCharged;
	308
	309	//Put new aod object in file in AOD calo clusters array
	310	AliAODCaloCluster caloCluster = new(((fOutputEvent->GetCaloClusters()))[naod++])
	311	AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0);
a3aebfff	312
	313
	314	caloCluster->SetPIDFromESD(pid);
477d6cee	315	caloCluster->SetCaloCluster(clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),
	316	clus->GetM20(), clus->GetM02(),
	317	clus->GetEmcCpvDistance(), clus->GetNExMax()) ;
	318
	319	if(fDebug > 3 && momentum.E() > 0.2)
a3aebfff	320	printf("AliCaloTrackESDReader::FillInputPHOS() - Selected clusters Distance BC %2.2f, dispersion %2.2f, M20 %f, M02 %3.2f, EmcCpvDist %3.3f, NexMax %d, TOF %e\n",
477d6cee	321	clus->GetDistanceToBadChannel(), clus->GetClusterDisp(),clus->GetM20(), clus->GetM02(),
	322	clus->GetEmcCpvDistance(), clus->GetNExMax(), clus->GetTOF());
	323
	324	caloCluster->SetNCells(clus->GetNCells());
	325	caloCluster->SetCellsAbsId(clus->GetCellsAbsId());
	326	caloCluster->SetCellsAmplitudeFraction(clus->GetCellsAmplitudeFraction());
a3aebfff	327
477d6cee	328	TArrayI* matchedT = clus->GetTracksMatched();
a3aebfff	329	if (nTracks > 0 && matchedT && clus->GetTrackMatched() >= 0) {
477d6cee	330	for (Int_t im = 0; im < matchedT->GetSize(); im++) {
a3aebfff	331	Int_t iESDtrack = matchedT->At(im);
	332	if(iESDtrack < nTracks && iESDtrack > -1)
	333	caloCluster->AddTrackMatched((fInputEvent->GetTrack(iESDtrack)));
477d6cee	334	}
	335	}
	336	}//Pt and Fidutial cut passed.
	337	}//PHOS cluster
	338	}//cluster exists
	339	}//esd cluster loop
	340
	341	//Put references to selected clusters in array
8896b363	342	for(Int_t iclus = nEMCAL; iclus < (fOutputEvent->GetCaloClusters())->GetEntriesFast(); iclus++){
477d6cee	343	AliAODCaloCluster * clus = (AliAODCaloCluster*) (fOutputEvent->GetCaloClusters())->At(iclus);
	344	fAODPHOS->Add(clus);
	345	}
dde5a268	346	if(fDebug > 1) printf("AliCaloTrackESDReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());
477d6cee	347
1c5acb87	348	}
	349
	350	//____________________________________________________________________________
	351	void AliCaloTrackESDReader::FillInputEMCALCells() {
477d6cee	352	//Return array with EMCAL cells in esd format
	353
	354	fEMCALCells = (TNamed) ((AliESDEvent)fInputEvent)->GetEMCALCells();
	355
1c5acb87	356	}
	357
	358	//____________________________________________________________________________
	359	void AliCaloTrackESDReader::FillInputPHOSCells() {
477d6cee	360	//Return array with PHOS cells in esd format
	361
	362	fPHOSCells = (TNamed) ((AliESDEvent)fInputEvent)->GetPHOSCells();
	363
1c5acb87	364	}
	365
	366	//____________________________________________________________________________
	367	void AliCaloTrackESDReader::GetVertex(Double_t v[3]) const {
477d6cee	368	//Return vertex position
477d6cee	369
8a587055	370	//((AliESDEvent*)fInputEvent)->GetVertex()->GetXYZ(v) ;//SPD
	371	((AliESDEvent*)fInputEvent)->GetPrimaryVertex()->GetXYZ(v);
	372
1c5acb87	373	}
	374
	375
8dacfd76	376	//____________________________________________________________________________
	377	Double_t AliCaloTrackESDReader::GetBField() const {
	378	//Return magnetic field
	379
	380	Double_t bfield = ((AliESDEvent*)fInputEvent)->GetMagneticField();
	381
	382	return bfield;
	383	}
	384
	385
1c5acb87	386	//____________________________________________________________________________
477d6cee	387	void AliCaloTrackESDReader::SetInputOutputMCEvent(AliVEvent* esd, AliAODEvent* aod, AliMCEvent* mc) {
	388	// Connect the data pointers
	389
	390	if(strcmp(esd->GetName(),"AliESDEvent")){
591cc579	391	printf("AliCaloTrackESDReader::SetInputOutputMCEvent() - STOP ::Wrong reader, here only ESDs. Input name: %s != AliESDEvent \n",esd->GetName());
477d6cee	392	abort();
	393	}
	394
	395	SetInputEvent(esd);
	396	SetOutputEvent(aod);
	397	SetMC(mc);
	398
1c5acb87	399	}