[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.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$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.1  2007/01/23 17:17:29  schutz
 * New Gamma package
 *
 *
 */

//_________________________________________________________________________
// Class for the analysis of gamma correlations (gamma-jet, 
// gamma-hadron and isolation cut.
// This class contains 3 main methods: one to fill lists of particles (ESDs) comming 
//  from the CTS (ITS+TPC) and the calorimeters;  the other one tags a candidate 
//  cluster as isolated;  the last one search in the 
//  corresponing calorimeter for the highest energy cluster, identified it as 
//  prompt photon;
//
//  Class created from old AliPHOSGammaJet 
//  (see AliRoot versions previous Release 4-09)
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---

#include <TFile.h>
#include <TParticle.h>
#include <TH2.h>
#include <TChain.h>
#include "AliAnaGammaDirect.h" 
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliESDCaloCluster.h"
#include "Riostream.h"
#include "AliLog.h"

ClassImp(AliAnaGammaDirect)

//____________________________________________________________________________
  AliAnaGammaDirect::AliAnaGammaDirect(const char *name) : 
    AliAnalysisTask(name,""), fChain(0), fESD(0),
    fOutputContainer(new TObjArray(100)), 
    fPrintInfo(0), fMinGammaPt(0.),
    fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
    fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
    fMakeICMethod(0)
{
  //Ctor        
  TList * list = gDirectory->GetListOfKeys() ; 
  TIter next(list) ; 
  TH2F * h = 0 ;
  Int_t index ; 
  for (index = 0 ; index < list->GetSize()-1 ; index++) { 
    //-1 to avoid GammaJet Task
    h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
    fOutputContainer->Add(h) ; 
  }
  
  // Input slot #0 works with an Ntuple
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TH1 container
  DefineOutput(0,  TObjArray::Class()) ; 
  
}


//____________________________________________________________________________
AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
  AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
  fOutputContainer(g. fOutputContainer),  fPrintInfo(g.fPrintInfo),
  fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
  fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
  fConeSize(g.fConeSize),
  fPtThreshold(g.fPtThreshold),
  fPtSumThreshold(g.fPtSumThreshold), 
  fMakeICMethod(g.fMakeICMethod)
{
  // cpy ctor
  SetName (g.GetName()) ; 
  SetTitle(g.GetTitle()) ; 
}

//____________________________________________________________________________
AliAnaGammaDirect::~AliAnaGammaDirect() 
{
  // Remove all pointers
  fOutputContainer->Clear() ; 
  delete fOutputContainer ;
  delete fhNGamma    ;  
  delete fhPhiGamma  ; 
  delete fhEtaGamma   ;  

}

//______________________________________________________________________________
void AliAnaGammaDirect::ConnectInputData(const Option_t*)
{
  // Initialisation of branch container and histograms 
    
  AliInfo(Form("*** Initialization of %s", GetName())) ; 
  
  // Get input data
  fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
  if (!fChain) {
    AliError(Form("Input 0 for %s not found\n", GetName()));
    return ;
  }
  
  // One should first check if the branch address was taken by some other task
  char ** address = (char **)GetBranchAddress(0, "ESD");
  if (address) {
    fESD = (AliESD*)(*address);
  } else {
    fESD = new AliESD();
    SetBranchAddress(0, "ESD", &fESD);
  }

}

//________________________________________________________________________
void AliAnaGammaDirect::CreateOutputObjects()
{  

  // Init parameteres and create histograms to be saved in output file and 
  // stores them in fOutputContainer
  InitParameters();

  fOutputContainer = new TObjArray(3) ;

  //Histograms of highest gamma identified in Event
  fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); 
  fhNGamma->SetYTitle("N");
  fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
  fOutputContainer->Add(fhNGamma) ; 
  
  fhPhiGamma  = new TH2F
    ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
  fhPhiGamma->SetYTitle("#phi");
  fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
  fOutputContainer->Add(fhPhiGamma) ; 
  
  fhEtaGamma  = new TH2F
    ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
  fhEtaGamma->SetYTitle("#eta");
  fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
  fOutputContainer->Add(fhEtaGamma) ;

}

//____________________________________________________________________________
void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl, 
					TClonesArray * plCTS, 
					TClonesArray * plEMCAL,  
					TClonesArray * plPHOS){
  
  //Create a list of particles from the ESD. These particles have been measured 
  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 
  
  Int_t index = pl->GetEntries() ; 
  Int_t npar  = 0 ;
  Float_t *pid = new Float_t[AliPID::kSPECIESN];  
  AliDebug(3,"Fill particle lists");
  if(fPrintInfo)
    AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));

  Double_t v[3] ; //vertex ;
  fESD->GetVertex()->GetXYZ(v) ; 

  //########### PHOS ##############
  if(fCalorimeter == "PHOS"){

    Int_t begphos = fESD->GetFirstPHOSCluster();  
    Int_t endphos = fESD->GetFirstPHOSCluster() + 
      fESD->GetNumberOfPHOSClusters() ;  
    Int_t indexNePHOS = plPHOS->GetEntries() ;
    AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));

    if(fCalorimeter == "PHOS"){
      for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
	AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd

	//Create a TParticle to fill the particle list
	
	Float_t en = clus->GetClusterEnergy() ;
	Float_t *p = new Float_t();
	clus->GetGlobalPosition(p) ;
	TVector3 pos(p[0],p[1],p[2]) ; 
	Double_t phi  = pos.Phi();
	Double_t theta= pos.Theta();
	Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
	Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
	Double_t pz = en*TMath::Cos(theta);

	TParticle * particle = new TParticle() ;
	particle->SetMomentum(px,py,pz,en) ;
	AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	
	//Select only photons
	
	pid=clus->GetPid();
	//cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
	if( !fPHOSPID)
	  new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
	else if( pid[AliPID::kPhoton] > 0.75)
	  new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
      }
    }
  }

  //########### CTS (TPC+ITS) #####################
  Int_t begtpc   = 0 ;  
  Int_t endtpc   = fESD->GetNumberOfTracks() ;
  Int_t indexCh  = plCTS->GetEntries() ;
  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
  
  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
    AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
    //We want tracks fitted in the detectors:
    ULong_t status=AliESDtrack::kTPCrefit;
    status|=AliESDtrack::kITSrefit;
   
    //We want tracks whose PID bit is set:
    //     ULong_t status =AliESDtrack::kITSpid;
    //     status|=AliESDtrack::kTPCpid;
  
    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
      // Do something with the tracks which were successfully
      // re-fitted 
      Double_t en = 0; //track ->GetTPCsignal() ;
      Double_t mom[3];
      track->GetPxPyPz(mom) ;
      Double_t px = mom[0];
      Double_t py = mom[1];
      Double_t pz = mom[2]; //Check with TPC people if this is correct.
      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
      //I just want to tag the particle as charged
       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
						 px, py, pz, en, v[0], v[1], v[2], 0);
  
      //TParticle * particle = new TParticle() ;
      //particle->SetMomentum(px,py,pz,en) ;
 
      new((*plCTS)[indexCh++])       TParticle(*particle) ;    
      new((*pl)[index++])           TParticle(*particle) ;
    }
  }
  
  //################ EMCAL ##############
  
  Int_t begem = fESD->GetFirstEMCALCluster();  
  Int_t endem = fESD->GetFirstEMCALCluster() + 
    fESD->GetNumberOfEMCALClusters() ;  
  Int_t indexNe  = plEMCAL->GetEntries() ; 
  
  AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
    
    for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
      AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
      Int_t clustertype= clus->GetClusterType();
      if(clustertype == AliESDCaloCluster::kClusterv1){
	Float_t en = clus->GetClusterEnergy() ;
	Float_t *p = new Float_t();
	clus->GetGlobalPosition(p) ;
	TVector3 pos(p[0],p[1],p[2]) ;
	Double_t phi  = pos.Phi();
	Double_t theta= pos.Theta();
	Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
	Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
	Double_t pz = en*TMath::Cos(theta);

	pid=clus->GetPid();
	if(fCalorimeter == "EMCAL")
	  {
	    TParticle * particle = new TParticle() ;
	    particle->SetMomentum(px,py,pz,en) ;
	    AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	    if(!fEMCALPID) //Only identified particles
	      new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
	    else if(pid[AliPID::kPhoton] > 0.75)
	      new((*plEMCAL)[indexNe++])       TParticle(*particle) ;	    
	  }
	else
	  {
	    Int_t pdg = 0;
	    if(fEMCALPID) 
	      {
		if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
		  pdg = 22;
		else if( pid[AliPID::kPi0] > 0.75)
		  pdg = 111;
	      }
	    else
	      pdg = 22; //No PID, assume all photons
	    
	    TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
						 px, py, pz, en, v[0], v[1], v[2], 0);
	    AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	    
	    new((*plEMCAL)[indexNe++])       TParticle(*particle) ; 
	    new((*pl)[index++])           TParticle(*particle) ;
	  }
      }
    }

    AliDebug(3,"Particle lists filled");
    
}


//____________________________________________________________________________
void AliAnaGammaDirect::Exec(Option_t *) 
{
  
  // Processing of one event

  //Get ESDs
  Long64_t entry = fChain->GetReadEntry() ;
  
  if (!fESD) {
    AliError("fESD is not connected to the input!") ; 
    return ; 
  }
  
  if (fPrintInfo) 
     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ; 

  //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.

  TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
  TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
  TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
  TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured in PHOS as Gamma calorimeter
  TClonesArray * plEMCAL   = new TClonesArray("TParticle",1000);  // All particles measured in EMCAL as Gamma calorimeter

  TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
 
  //Fill lists with photons, neutral particles and charged particles
  //look for the highest energy photon in the event inside fCalorimeter
  
  AliDebug(2, "Fill particle lists, get prompt gamma");
  
  //Fill particle lists 
  CreateParticleList(particleList, plCTS,plEMCAL,plPHOS); 
  
  if(fCalorimeter == "PHOS")
    plNe = plPHOS;
  if(fCalorimeter == "EMCAL")
    plNe = plEMCAL;
  
  
  Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
  //Search highest energy prompt gamma in calorimeter
  GetPromptGamma(plNe,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
  
  AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
  
  //If there is any photon candidate in fCalorimeter
  if(iIsInPHOSorEMCAL){
    if (fPrintInfo)
      AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
    
  }//Gamma in Calo
  
  AliDebug(2, "End of analysis, delete pointers");
  
  particleList->Delete() ; 
  plCTS->Delete() ;
  plNe->Delete() ;
  plEMCAL->Delete() ;
  plPHOS->Delete() ;
  pGamma->Delete();
 
  PostData(0, fOutputContainer);

 delete plNe ;
  delete plCTS ;
  //delete plPHOS ;
  //delete plEMCAL ;
  delete particleList ;

  //  delete pGamma;

}    


//____________________________________________________________________________
void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const 
{
  //Search for the prompt photon in Calorimeter with pt > fMinGammaPt

  Double_t pt = 0;
  Int_t index = -1; 
  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;

    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
      index = ipr ;
      pt = particle->Pt();
      pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
      AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
      Is  = kTRUE;
    }
  }

  //Do Isolation?
  if(fMakeICMethod && Is)
    {
      Float_t coneptsum = 0 ;
      Bool_t  icPtThres   = kFALSE;
      Bool_t  icPtSum     = kFALSE;
      MakeIsolationCut(plCTS,pl, pGamma, index, 
		       icPtThres, icPtSum,coneptsum);
      if(fMakeICMethod == 1) //Pt thres method
	Is = icPtThres ;
      if(fMakeICMethod == 2) //Pt cone sum method
	Is = icPtSum ;
    }
  
  if(Is){
    AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
    AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
    //Fill prompt gamma histograms
    fhNGamma  ->Fill(pGamma->Pt());
    fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
    fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
  }
  else
    AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
}

  //____________________________________________________________________________
void AliAnaGammaDirect::InitParameters()
{
 
  //Initialize the parameters of the analysis.
  fCalorimeter="PHOS";
  fPrintInfo           = kTRUE;
  fMinGammaPt  = 5. ;

  //Fill particle lists when PID is ok
  fEMCALPID = kFALSE;
  fPHOSPID = kFALSE;

  fConeSize             = 0.2 ; 
  fPtThreshold         = 2.0; 
  fPtSumThreshold  = 1.; 

  fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
}

//__________________________________________________________________
void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
					   TClonesArray * plNe, 
					   TParticle * pCandidate, 
					   Int_t index, 
					   Bool_t  &icmpt,  Bool_t  &icms, 
					   Float_t &coneptsum) const 
{  
  //Search in cone around a candidate particle if it is isolated 
  Float_t phiC  = pCandidate->Phi() ;
  Float_t etaC = pCandidate->Eta() ;
  Float_t pt     = -100. ;
  Float_t eta   = -100.  ;
  Float_t phi    = -100.  ;
  Float_t rad   = -100 ;
  Int_t    n        = 0 ;
  TParticle * particle  = new TParticle();

  coneptsum = 0; 
  icmpt = kFALSE;
  icms  = kFALSE;

  //Check charged particles in cone.
  for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
    particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
    pt    = particle->Pt();
    eta  = particle->Eta();
    phi  = particle->Phi() ;
    
    //Check if there is any particle inside cone with pt larger than  fPtThreshold
    rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
		      TMath::Power((phi-phiC),2));
    if(rad<fConeSize){
      AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
      coneptsum+=pt;
      if(pt > fPtThreshold ) n++;
    }
  }// charged particle loop
  
  //Check neutral particles in cone.
  for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
    if(ipr != index){//Do not count the candidate
      particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
      pt    = particle->Pt();
      eta  = particle->Eta();
      phi  = particle->Phi() ;
      
      //Check if there is any particle inside cone with pt larger than  fPtThreshold
      rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
			TMath::Power((phi-phiC),2));
      if(rad<fConeSize){
	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	coneptsum+=pt;
	if(pt > fPtThreshold ) n++;
      }
    }
  }// neutral particle loop
  
  if(n == 0) 
    icmpt =  kTRUE ;
  if(coneptsum < fPtSumThreshold)
    icms  =  kTRUE ;

}

void AliAnaGammaDirect::Print(const Option_t * opt) const
{

  //Print some relevant parameters set for the analysis
  if(! opt)
    return;

  Info("Print", "%s %s", GetName(), GetTitle() ) ;
  printf("IC method               =     %d\n", fMakeICMethod) ; 
  printf("Cone Size               =     %f\n", fConeSize) ; 
  printf("pT threshold           =     %f\n", fPtThreshold) ;
  printf("pT sum threshold   =     %f\n", fPtSumThreshold) ; 
  printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ; 
  printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ; 
} 

void AliAnaGammaDirect::Terminate(Option_t *)
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
    

}
Commit	Line	Data
f9cea31c	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	/* $Id$ */
	16
	17	/* History of cvs commits:
	18	*
	19	* $Log$
2a1d8a29	20	* Revision 1.1 2007/01/23 17:17:29 schutz
	21	* New Gamma package
	22	*
f9cea31c	23	*
	24	*/
	25
	26	//_________________________________________________________________________
	27	// Class for the analysis of gamma correlations (gamma-jet,
	28	// gamma-hadron and isolation cut.
	29	// This class contains 3 main methods: one to fill lists of particles (ESDs) comming
	30	// from the CTS (ITS+TPC) and the calorimeters; the other one tags a candidate
	31	// cluster as isolated; the last one search in the
	32	// corresponing calorimeter for the highest energy cluster, identified it as
	33	// prompt photon;
	34	//
	35	// Class created from old AliPHOSGammaJet
	36	// (see AliRoot versions previous Release 4-09)
	37	//
	38	//*-- Author: Gustavo Conesa (LNF-INFN)
	39	//////////////////////////////////////////////////////////////////////////////
	40
	41
	42	// --- ROOT system ---
	43
	44	#include <TFile.h>
	45	#include <TParticle.h>
	46	#include <TH2.h>
2a1d8a29	47	#include <TChain.h>
f9cea31c	48	#include "AliAnaGammaDirect.h"
	49	#include "AliESD.h"
	50	#include "AliESDtrack.h"
	51	#include "AliESDCaloCluster.h"
	52	#include "Riostream.h"
	53	#include "AliLog.h"
	54
	55	ClassImp(AliAnaGammaDirect)
	56
	57	//____________________________________________________________________________
	58	AliAnaGammaDirect::AliAnaGammaDirect(const char *name) :
	59	AliAnalysisTask(name,""), fChain(0), fESD(0),
	60	fOutputContainer(new TObjArray(100)),
	61	fPrintInfo(0), fMinGammaPt(0.),
	62	fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
	63	fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
2a1d8a29	64	fMakeICMethod(0)
f9cea31c	65	{
	66	//Ctor
	67	TList * list = gDirectory->GetListOfKeys() ;
	68	TIter next(list) ;
	69	TH2F * h = 0 ;
	70	Int_t index ;
	71	for (index = 0 ; index < list->GetSize()-1 ; index++) {
	72	//-1 to avoid GammaJet Task
	73	h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
	74	fOutputContainer->Add(h) ;
	75	}
	76
f9cea31c	77	// Input slot #0 works with an Ntuple
	78	DefineInput(0, TChain::Class());
	79	// Output slot #0 writes into a TH1 container
	80	DefineOutput(0, TObjArray::Class()) ;
	81
	82	}
	83
	84
	85	//____________________________________________________________________________
	86	AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
	87	AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
	88	fOutputContainer(g. fOutputContainer), fPrintInfo(g.fPrintInfo),
	89	fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
	90	fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
	91	fConeSize(g.fConeSize),
	92	fPtThreshold(g.fPtThreshold),
	93	fPtSumThreshold(g.fPtSumThreshold),
f9cea31c	94	fMakeICMethod(g.fMakeICMethod)
	95	{
	96	// cpy ctor
	97	SetName (g.GetName()) ;
	98	SetTitle(g.GetTitle()) ;
f9cea31c	99	}
	100
	101	//____________________________________________________________________________
	102	AliAnaGammaDirect::~AliAnaGammaDirect()
	103	{
	104	// Remove all pointers
	105	fOutputContainer->Clear() ;
	106	delete fOutputContainer ;
	107	delete fhNGamma ;
	108	delete fhPhiGamma ;
	109	delete fhEtaGamma ;
2a1d8a29	110
	111	}
	112
	113	//______________________________________________________________________________
	114	void AliAnaGammaDirect::ConnectInputData(const Option_t*)
	115	{
	116	// Initialisation of branch container and histograms
	117
	118	AliInfo(Form("*** Initialization of %s", GetName())) ;
	119
	120	// Get input data
	121	fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
	122	if (!fChain) {
	123	AliError(Form("Input 0 for %s not found\n", GetName()));
	124	return ;
	125	}
	126
	127	// One should first check if the branch address was taken by some other task
	128	char address = (char )GetBranchAddress(0, "ESD");
	129	if (address) {
	130	fESD = (AliESD)(address);
	131	} else {
	132	fESD = new AliESD();
	133	SetBranchAddress(0, "ESD", &fESD);
	134	}
	135
	136	}
	137
	138	//________________________________________________________________________
	139	void AliAnaGammaDirect::CreateOutputObjects()
	140	{
	141
	142	// Init parameteres and create histograms to be saved in output file and
	143	// stores them in fOutputContainer
	144	InitParameters();
	145
	146	fOutputContainer = new TObjArray(3) ;
	147
	148	//Histograms of highest gamma identified in Event
	149	fhNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
	150	fhNGamma->SetYTitle("N");
	151	fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
	152	fOutputContainer->Add(fhNGamma) ;
	153
	154	fhPhiGamma = new TH2F
	155	("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
	156	fhPhiGamma->SetYTitle("#phi");
	157	fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
	158	fOutputContainer->Add(fhPhiGamma) ;
	159
	160	fhEtaGamma = new TH2F
	161	("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
	162	fhEtaGamma->SetYTitle("#eta");
	163	fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
	164	fOutputContainer->Add(fhEtaGamma) ;
f9cea31c	165
	166	}
	167
	168	//____________________________________________________________________________
	169	void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl,
	170	TClonesArray * plCTS,
	171	TClonesArray * plEMCAL,
	172	TClonesArray * plPHOS){
	173
	174	//Create a list of particles from the ESD. These particles have been measured
	175	//by the Central Tracking system (TPC+ITS), PHOS and EMCAL
	176
	177	Int_t index = pl->GetEntries() ;
	178	Int_t npar = 0 ;
	179	Float_t *pid = new Float_t[AliPID::kSPECIESN];
	180	AliDebug(3,"Fill particle lists");
	181	if(fPrintInfo)
	182	AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
	183
	184	Double_t v[3] ; //vertex ;
	185	fESD->GetVertex()->GetXYZ(v) ;
	186
	187	//########### PHOS ##############
	188	if(fCalorimeter == "PHOS"){
	189
	190	Int_t begphos = fESD->GetFirstPHOSCluster();
	191	Int_t endphos = fESD->GetFirstPHOSCluster() +
	192	fESD->GetNumberOfPHOSClusters() ;
	193	Int_t indexNePHOS = plPHOS->GetEntries() ;
	194	AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
	195
	196	if(fCalorimeter == "PHOS"){
	197	for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
	198	AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
	199
	200	//Create a TParticle to fill the particle list
	201
	202	Float_t en = clus->GetClusterEnergy() ;
	203	Float_t *p = new Float_t();
	204	clus->GetGlobalPosition(p) ;
	205	TVector3 pos(p[0],p[1],p[2]) ;
	206	Double_t phi = pos.Phi();
	207	Double_t theta= pos.Theta();
	208	Double_t px = enTMath::Cos(phi)TMath::Sin(theta);;
	209	Double_t py = enTMath::Sin(phi)TMath::Sin(theta);
	210	Double_t pz = en*TMath::Cos(theta);
	211
	212	TParticle * particle = new TParticle() ;
	213	particle->SetMomentum(px,py,pz,en) ;
	214	AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
	215
	216	//Select only photons
	217
	218	pid=clus->GetPid();
	219	//cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
	220	if( !fPHOSPID)
	221	new((plPHOS)[indexNePHOS++]) TParticle(particle) ;
	222	else if( pid[AliPID::kPhoton] > 0.75)
	223	new((plPHOS)[indexNePHOS++]) TParticle(particle) ;
	224	}
	225	}
	226	}
	227
	228	//########### CTS (TPC+ITS) #####################
229	Int_t begtpc = 0 ;
230	Int_t endtpc = fESD->GetNumberOfTracks() ;
231	Int_t indexCh = plCTS->GetEntries() ;
232	AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
233
234	for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
235	AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
236	//We want tracks fitted in the detectors:
237	ULong_t status=AliESDtrack::kTPCrefit;
238	status\|=AliESDtrack::kITSrefit;
239
240	//We want tracks whose PID bit is set:
241	// ULong_t status =AliESDtrack::kITSpid;
242	// status\|=AliESDtrack::kTPCpid;
243
244	if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
245	// Do something with the tracks which were successfully
246	// re-fitted
247	Double_t en = 0; //track ->GetTPCsignal() ;
248	Double_t mom[3];
249	track->GetPxPyPz(mom) ;
250	Double_t px = mom[0];
251	Double_t py = mom[1];
252	Double_t pz = mom[2]; //Check with TPC people if this is correct.
253	Int_t pdg = 11; //Give any charged PDG code, in this case electron.
254	//I just want to tag the particle as charged
255	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
256	px, py, pz, en, v[0], v[1], v[2], 0);
257
258	//TParticle * particle = new TParticle() ;
259	//particle->SetMomentum(px,py,pz,en) ;
260
261	new((plCTS)[indexCh++]) TParticle(particle) ;
262	new((pl)[index++]) TParticle(particle) ;
263	}
264	}
265
266	//################ EMCAL ##############
267
268	Int_t begem = fESD->GetFirstEMCALCluster();
269	Int_t endem = fESD->GetFirstEMCALCluster() +
270	fESD->GetNumberOfEMCALClusters() ;
271	Int_t indexNe = plEMCAL->GetEntries() ;
272
273	AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
274
275	for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
276	AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
277	Int_t clustertype= clus->GetClusterType();
278	if(clustertype == AliESDCaloCluster::kClusterv1){
279	Float_t en = clus->GetClusterEnergy() ;
280	Float_t *p = new Float_t();
281	clus->GetGlobalPosition(p) ;
282	TVector3 pos(p[0],p[1],p[2]) ;
283	Double_t phi = pos.Phi();
284	Double_t theta= pos.Theta();
285	Double_t px = enTMath::Cos(phi)TMath::Sin(theta);;
286	Double_t py = enTMath::Sin(phi)TMath::Sin(theta);
287	Double_t pz = en*TMath::Cos(theta);
288
289	pid=clus->GetPid();
290	if(fCalorimeter == "EMCAL")
291	{
292	TParticle * particle = new TParticle() ;
293	particle->SetMomentum(px,py,pz,en) ;
294	AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
295	if(!fEMCALPID) //Only identified particles
296	new((plEMCAL)[indexNe++]) TParticle(particle) ;
297	else if(pid[AliPID::kPhoton] > 0.75)
298	new((plEMCAL)[indexNe++]) TParticle(particle) ;
299	}
300	else
301	{
302	Int_t pdg = 0;
303	if(fEMCALPID)
304	{
305	if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
306	pdg = 22;
307	else if( pid[AliPID::kPi0] > 0.75)
308	pdg = 111;
309	}
310	else
311	pdg = 22; //No PID, assume all photons
312
313	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
314	px, py, pz, en, v[0], v[1], v[2], 0);
315	AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
316
317	new((plEMCAL)[indexNe++]) TParticle(particle) ;
318	new((pl)[index++]) TParticle(particle) ;
319	}
320	}
321	}
322
323	AliDebug(3,"Particle lists filled");
324
325	}
326
327
328
329	//____________________________________________________________________________
330	void AliAnaGammaDirect::Exec(Option_t *)
331	{
332
333	// Processing of one event
2a1d8a29	334
f9cea31c	335	//Get ESDs
	336	Long64_t entry = fChain->GetReadEntry() ;
	337
	338	if (!fESD) {
	339	AliError("fESD is not connected to the input!") ;
	340	return ;
	341	}
	342
	343	if (fPrintInfo)
	344	AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
	345
	346	//CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
	347
	348	TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
	349	TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
	350	TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter (EMCAL)
	351	TClonesArray * plPHOS = new TClonesArray("TParticle",1000); // All particles measured in PHOS as Gamma calorimeter
	352	TClonesArray * plEMCAL = new TClonesArray("TParticle",1000); // All particles measured in EMCAL as Gamma calorimeter
	353
	354	TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
2a1d8a29	355
f9cea31c	356	//Fill lists with photons, neutral particles and charged particles
f9cea31c	357	//look for the highest energy photon in the event inside fCalorimeter
2a1d8a29	358
f9cea31c	359	AliDebug(2, "Fill particle lists, get prompt gamma");
2a1d8a29	360
f9cea31c	361	//Fill particle lists
f9cea31c	362	CreateParticleList(particleList, plCTS,plEMCAL,plPHOS);
2a1d8a29	363
f9cea31c	364	if(fCalorimeter == "PHOS")
	365	plNe = plPHOS;
	366	if(fCalorimeter == "EMCAL")
	367	plNe = plEMCAL;
2a1d8a29	368
	369
	370	Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
	371	//Search highest energy prompt gamma in calorimeter
	372	GetPromptGamma(plNe, plCTS, pGamma, iIsInPHOSorEMCAL) ;
	373
	374	AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
	375
	376	//If there is any photon candidate in fCalorimeter
	377	if(iIsInPHOSorEMCAL){
	378	if (fPrintInfo)
	379	AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
f9cea31c	380
2a1d8a29	381	}//Gamma in Calo
f9cea31c	382
	383	AliDebug(2, "End of analysis, delete pointers");
	384
	385	particleList->Delete() ;
	386	plCTS->Delete() ;
	387	plNe->Delete() ;
2a1d8a29	388	plEMCAL->Delete() ;
f9cea31c	389	plPHOS->Delete() ;
f9cea31c	390	pGamma->Delete();
2a1d8a29	391
2a1d8a29	392	PostData(0, fOutputContainer);
f9cea31c	393
2a1d8a29	394	delete plNe ;
f9cea31c	395	delete plCTS ;
2a1d8a29	396	//delete plPHOS ;
2a1d8a29	397	//delete plEMCAL ;
f9cea31c	398	delete particleList ;
2a1d8a29	399
f9cea31c	400	// delete pGamma;
f9cea31c	401
f9cea31c	402	}
	403
	404
	405	//____________________________________________________________________________
	406	void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const
	407	{
	408	//Search for the prompt photon in Calorimeter with pt > fMinGammaPt
	409
	410	Double_t pt = 0;
	411	Int_t index = -1;
	412	for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
	413	TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
	414
	415	if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
	416	index = ipr ;
	417	pt = particle->Pt();
	418	pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
	419	AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
	420	Is = kTRUE;
	421	}
	422	}
	423
	424	//Do Isolation?
	425	if(fMakeICMethod && Is)
	426	{
	427	Float_t coneptsum = 0 ;
	428	Bool_t icPtThres = kFALSE;
	429	Bool_t icPtSum = kFALSE;
	430	MakeIsolationCut(plCTS,pl, pGamma, index,
	431	icPtThres, icPtSum,coneptsum);
	432	if(fMakeICMethod == 1) //Pt thres method
	433	Is = icPtThres ;
	434	if(fMakeICMethod == 2) //Pt cone sum method
	435	Is = icPtSum ;
	436	}
	437
	438	if(Is){
	439	AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
	440	AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
	441	//Fill prompt gamma histograms
	442	fhNGamma ->Fill(pGamma->Pt());
	443	fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
	444	fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
	445	}
	446	else
	447	AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
	448	}
	449
	450	//____________________________________________________________________________
2a1d8a29	451	void AliAnaGammaDirect::InitParameters()
f9cea31c	452	{
f9cea31c	453
2a1d8a29	454	//Initialize the parameters of the analysis.
f9cea31c	455	fCalorimeter="PHOS";
f9cea31c	456	fPrintInfo = kTRUE;
2a1d8a29	457	fMinGammaPt = 5. ;
f9cea31c	458
	459	//Fill particle lists when PID is ok
	460	fEMCALPID = kFALSE;
	461	fPHOSPID = kFALSE;
	462
	463	fConeSize = 0.2 ;
	464	fPtThreshold = 2.0;
	465	fPtSumThreshold = 1.;
	466
2a1d8a29	467	fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
f9cea31c	468	}
	469
	470	//__________________________________________________________________
	471	void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
	472	TClonesArray * plNe,
	473	TParticle * pCandidate,
	474	Int_t index,
	475	Bool_t &icmpt, Bool_t &icms,
	476	Float_t &coneptsum) const
	477	{
	478	//Search in cone around a candidate particle if it is isolated
	479	Float_t phiC = pCandidate->Phi() ;
	480	Float_t etaC = pCandidate->Eta() ;
	481	Float_t pt = -100. ;
	482	Float_t eta = -100. ;
	483	Float_t phi = -100. ;
	484	Float_t rad = -100 ;
	485	Int_t n = 0 ;
	486	TParticle * particle = new TParticle();
	487
	488	coneptsum = 0;
	489	icmpt = kFALSE;
	490	icms = kFALSE;
	491
	492	//Check charged particles in cone.
	493	for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
	494	particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
	495	pt = particle->Pt();
	496	eta = particle->Eta();
	497	phi = particle->Phi() ;
	498
	499	//Check if there is any particle inside cone with pt larger than fPtThreshold
	500	rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
	501	TMath::Power((phi-phiC),2));
	502	if(rad<fConeSize){
	503	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	504	coneptsum+=pt;
	505	if(pt > fPtThreshold ) n++;
	506	}
	507	}// charged particle loop
	508
	509	//Check neutral particles in cone.
	510	for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
	511	if(ipr != index){//Do not count the candidate
	512	particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
	513	pt = particle->Pt();
	514	eta = particle->Eta();
	515	phi = particle->Phi() ;
	516
	517	//Check if there is any particle inside cone with pt larger than fPtThreshold
	518	rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
	519	TMath::Power((phi-phiC),2));
	520	if(rad<fConeSize){
	521	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	522	coneptsum+=pt;
	523	if(pt > fPtThreshold ) n++;
	524	}
	525	}
	526	}// neutral particle loop
	527
	528	if(n == 0)
	529	icmpt = kTRUE ;
	530	if(coneptsum < fPtSumThreshold)
	531	icms = kTRUE ;
532
533	}
534
f9cea31c	535	void AliAnaGammaDirect::Print(const Option_t * opt) const
	536	{
	537
	538	//Print some relevant parameters set for the analysis
	539	if(! opt)
	540	return;
	541
	542	Info("Print", "%s %s", GetName(), GetTitle() ) ;
2a1d8a29	543	printf("IC method = %d\n", fMakeICMethod) ;
f9cea31c	544	printf("Cone Size = %f\n", fConeSize) ;
	545	printf("pT threshold = %f\n", fPtThreshold) ;
	546	printf("pT sum threshold = %f\n", fPtSumThreshold) ;
2a1d8a29	547	printf("Min Gamma pT = %f\n", fMinGammaPt) ;
f9cea31c	548	printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
	549	}
	550
	551	void AliAnaGammaDirect::Terminate(Option_t *)
	552	{
	553	// The Terminate() function is the last function to be called during
	554	// a query. It always runs on the client, it can be used to present
	555	// the results graphically or save the results to file.
	556
	557
	558	}