[u/mrichter/AliRoot.git] / PWG4 / AliAnaGamma.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.4  2007/10/29 13:48:42  gustavo
 * Corrected coding violations
 *
 * Revision 1.2  2007/08/17 12:40:04  schutz
 * New analysis classes by Gustavo Conesa
 *
 * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
 * new analysis classes in the the new analysis framework
 *
 *
 */

//_________________________________________________________________________
// Base class for gamma and correlation analysis
// It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
// got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes 
// AliAnaGammaDirect, AliAnaGammaCorrelation ....
//
//*-- Author: Gustavo Conesa (INFN-LNF)

// --- ROOT system ---

#include <TParticle.h>
#include <TH2.h>

//---- AliRoot system ---- 
#include "AliAnaGamma.h" 
#include "AliGammaReader.h" 
#include "AliAnaGammaDirect.h" 
#include "AliAnaGammaCorrelation.h" 
#include "AliAnaGammaSelection.h" 
#include "AliNeutralMesonSelection.h"
#include "AliAODCaloCluster.h"
#include "AliAODEvent.h"
#include "Riostream.h"
#include "AliLog.h"

ClassImp(AliAnaGamma)


//____________________________________________________________________________
  AliAnaGamma::AliAnaGamma() : 
    TObject(),
    fOutputContainer(0x0), 
    fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0), 
    fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
    fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
{
  //Default Ctor

  //Initialize parameters, pointers and histograms
  if(!fReader)
    fReader = new AliGammaReader();
  if(!fGammaDirect)
    fGammaDirect = new AliAnaGammaDirect();
  if(!fGammaCorrelation)
    fGammaCorrelation = new AliAnaGammaCorrelation();
  if(!fGammaSelection)
    fGammaSelection = new AliAnaGammaSelection();
  if(!fNeutralMesonSelection)
    fNeutralMesonSelection = new AliNeutralMesonSelection();

  fAODclusters = 0;

  InitParameters();
  
}

//____________________________________________________________________________
AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :   
  TObject(),
  fOutputContainer(g. fOutputContainer), 
  fAnaType(g.fAnaType),  fCalorimeter(g.fCalorimeter), 
  fData(g.fData), fKine(g.fKine),fReader(g.fReader),
  fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
  fGammaSelection(g.fGammaSelection),
  fNeutralMesonSelection(g.fNeutralMesonSelection),  
  fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
{
  // cpy ctor
  
}

//_________________________________________________________________________
AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
{
  // assignment operator

  if(this == &source)return *this;
  ((TObject *)this)->operator=(source);

  fOutputContainer = source.fOutputContainer ;
  fAnaType = source.fAnaType;
  fCalorimeter = source.fCalorimeter ;
  fData = source.fData ; 
  fKine = source.fKine ;
  fReader = source.fReader ;
  fGammaDirect = source.fGammaDirect ;
  fGammaCorrelation = source.fGammaCorrelation ;
  fGammaSelection = source.fGammaSelection ;
  fNeutralMesonSelection = source.fNeutralMesonSelection ;

  return *this;

}

//____________________________________________________________________________
AliAnaGamma::~AliAnaGamma() 
{
  // Remove all pointers.

  // Protection added in case of NULL pointers (MG)
  if (fOutputContainer) {
     fOutputContainer->Clear();
     delete fOutputContainer ;
  }   

  if (fData) delete fData ; 
  if (fKine) delete fKine ;
  if (fReader) delete fReader ;
  if (fGammaDirect) delete fGammaDirect ;
  if (fGammaCorrelation) delete fGammaCorrelation ;
  if (fGammaSelection) delete fGammaSelection ;
  if (fNeutralMesonSelection) delete fNeutralMesonSelection ;

}

//________________________________________________________________________
void AliAnaGamma::Init()
{  

  //Init container histograms and other common variables

  //Histograms container
  fOutputContainer = new TList ;
  
  //Fill container with appropriate histograms
  
  //Set selection  analysis histograms
  TList * selectcontainer =  fGammaSelection->GetCreateOutputObjects(); 
  for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
    Bool_t  add = kTRUE ;
    TString name = (selectcontainer->At(i))->GetName();   
    if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
    if(!fReader->IsPHOSOn() && name.Contains("PHOS"))   add = kFALSE;
    if(!fReader->IsCTSOn() &&  !fGammaSelection->FillCTS() && name.Contains("CTS"))   add = kFALSE;
    if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
  }  //Set selection  analysis histograms
 
  
  //Set prompt photon analysis histograms
  TList * promptcontainer =  fGammaDirect->GetCreateOutputObjects(); 
  for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
    fOutputContainer->Add(promptcontainer->At(i)) ;
  
  //Check if selected options are correct or set them when necessary
  if(fReader->GetDataType() == AliGammaReader::kMCData){
    fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
    fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
  }
  
  if(fAnaType == kCorrelation){
    
    //Check if selected options are correct
    if (fGammaDirect->GetICMethod()==AliAnaGammaDirect::kSeveralIC)
      AliFatal("Correlation not allowed with multiple isolation cuts, kCorrelation and kSeveralIC do not go together");
    
    if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kParton && 
       fReader->GetDataType() != AliGammaReader::kMC)
      AliFatal("kParton must be analyzed with data kMC");
    
    //Set the parameters for the neutral pair selection depending on the analysis, 
    fNeutralMesonSelection->SetDeltaPhiCutRange(fGammaCorrelation->GetDeltaPhiMinCut(), 
						fGammaCorrelation->GetDeltaPhiMaxCut());  
    

    if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kHadron){
      fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiMinPt);
      fNeutralMesonSelection->SetMinPt(fGammaCorrelation->GetMinPtHadron());
      
    }
    
    if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kJetLeadCone){
      fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiPtRatio);
      fNeutralMesonSelection->SetRatioCutRange(fGammaCorrelation->GetRatioMinCut(), 
					       fGammaCorrelation->GetRatioMaxCut());
    }

    //Set the neutral mesosn selection histograms
    TList * neutralmesoncontainer =  fNeutralMesonSelection->GetCreateOutputObjects();
    if(fNeutralMesonSelection->AreNeutralMesonSelectionHistosKept()){
      for(Int_t i = 0; i < neutralmesoncontainer->GetEntries(); i++)
	fOutputContainer->Add(neutralmesoncontainer->At(i)) ;
    }
    
    //Set correlation histograms
    TList * correlationcontainer =  fGammaCorrelation->GetCreateOutputObjects();
    for(Int_t i = 0; i < correlationcontainer->GetEntries(); i++)
      fOutputContainer->Add(correlationcontainer->At(i)) ;
    fGammaCorrelation->SetOutputContainer(fOutputContainer);
    fGammaCorrelation->SetNeutralMesonSelection(fNeutralMesonSelection);
 
  }//kCorrelation  
  
}

//____________________________________________________________________________
void AliAnaGamma::InitParameters()
{

  //Init data members
 
  fAnaType = kPrompt;
  fCalorimeter = "EMCAL";

}

//__________________________________________________________________
void AliAnaGamma::Print(const Option_t * opt) const
{

  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  Info("Print", "%s %s", GetName(), GetTitle() ) ;
  printf("Analysis type           =     %d\n", fAnaType) ;
  printf("Calorimeter           =     %s\n", fCalorimeter.Data()) ;

  switch(fAnaType)
    {
    case kPrompt:
      {
	fGammaDirect->Print("");
      }// case kIsolationCut
      break;
   
    case kCorrelation:
      {
	fGammaCorrelation->Print("");
      }//  case kCorrelation
      break;
      
    }//switch
} 


//____________________________________________________________________________
Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){

  AliDebug(1,Form("Entry %d",entry));
  //cout<<"Event >>>>>>>>>>>>> "<<entry<<endl;
  if(!fOutputContainer)
    AliFatal("Histograms not initialized");

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

  TClonesArray * plCTS      = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
  TClonesArray * plEMCAL    = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
  TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured  Gamma calorimeter
  TClonesArray * plPrimCTS      = new TClonesArray("TParticle",1000); // primary tracks
  TClonesArray * plPrimEMCAL    = new TClonesArray("TParticle",1000);   // primary emcal clusters
  TClonesArray * plPrimPHOS     = new TClonesArray("TParticle",1000);  // primary phos clusters
  TClonesArray * plParton   = new TClonesArray("TParticle",1000);  // All partons
  //Fill lists with photons, neutral particles and charged particles
  //look for the highest energy photon in the event inside fCalorimeter
    
  //Fill particle lists 
  if(fReader->GetDataType() == AliGammaReader::kData){
    AliDebug(1,"Data analysis");
    fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL); 
  }
  else if( fReader->GetDataType()== AliGammaReader::kMC){
    AliDebug(1,"Kinematics analysis");
    fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL); 
  }
  else if(fReader->GetDataType() == AliGammaReader::kMCData) {
   AliDebug(1,"Data + Kinematics analysis");
   fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS); 
  }
  else
    AliError("Option not implemented");

  //Fill AOD with calorimeter
  //Temporal solution, just for testing
  //FillAODs(plPHOS,plEMCAL);  

  //Select particles to do the final analysis.
  if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
  if(fReader->IsPHOSOn())  fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
  if(fReader->IsCTSOn() &&  fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);


  //Search highest energy prompt gamma in calorimeter
  if(fCalorimeter == "PHOS")
    MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ; 
  else if (fCalorimeter == "EMCAL")
    MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ; 
  else
    AliFatal("Wrong calorimeter name");

  plCTS->Clear() ;
  plEMCAL->Clear() ;
  plPHOS->Clear() ;
  plParton->Clear() ;
  plPrimCTS->Clear() ;
  plPrimEMCAL->Clear() ;
  plPrimPHOS->Clear() ;

  delete plCTS ;
  delete plPHOS ;
  delete plEMCAL ;
  delete plParton ;
  delete plPrimCTS ;
  delete plPrimPHOS ;
  delete plPrimEMCAL ;

  return kTRUE;

}

//____________________________________________________________________________
void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo)  {
  
  TParticle * pGamma = new TParticle ;
  Bool_t isInCalo = kFALSE ;
  
  switch(fAnaType)
    {
      
      //Only Prompt photon analysis
    case kPrompt:
      {	  
	AliDebug(1,"kPrompt analysis");
	switch(fGammaDirect->GetICMethod())
	  {
	    
	  case AliAnaGammaDirect::kSeveralIC:
	    {
	      fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
	      AliDebug(1,"kSeveralIC analysis");
	    }
	    break;
	    
	  default :
	    {
	      fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
	      if(!isInCalo)
		AliDebug(1,"Prompt gamma not found");
	    }
	    break;
	  }//IC method
      }// case kPrompt:
      break;
      
      //Correlate prompt photon with something: parton, hadron, jet.
    case kCorrelation:
      {
	AliDebug(1,"kCorrelation analysis");
	//Find prompt photon	
	fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);

	if(isInCalo){//If prompt photon found, do correlation
	  
	  switch(fGammaCorrelation->GetCorrelationType())
	    {
	    case AliAnaGammaCorrelation::kParton:
	      {
		AliDebug(1,"kParton correlation");
		fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
	      }//  case kParton
	      break;
      
	    case AliAnaGammaCorrelation::kHadron:
	      {
		AliDebug(1,"kHadron correlation");
		fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	      }//  case kHadron
	      break;

	    case AliAnaGammaCorrelation::kJetLeadCone:
	      {		
		AliDebug(1,"kJetLeadCone correlation");
		fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	      }//  case kJetLeadCone
	      break;
	      
	    case AliAnaGammaCorrelation::kJetFinder:
	      {	
		AliDebug(1,"kJetFinder correlation");
		printf("Analysis not implemented \n");
	      }//  case kJetFinder
	      break;
	    }// switch correlation
	}// is in calo
	else  AliDebug(2,"Prompt gamma not found");
      }// case kCorrelation
      break;

    } //switch(fAnaType)

  delete pGamma ; 
  
}

//____________________________________________________
void AliAnaGamma::AddCluster(AliAODCaloCluster p)
{
// Add new jet to the list
  cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
  cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
  new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
}

//___________________________________________________
void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
{
// Connect to the AOD
  fAODclusters = aod->GetCaloClusters();
}

//____________________________________________________________________________
void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
  //Fill AOD caloClusters
  //Temporal method, just for testing AOD creation
  
  //Fill AOD with PHOS clusters
  Int_t nphos =  plPHOS->GetEntries() ;
  cout<<"PHOS entries "<<nphos<<endl;
  for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
    TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
    Float_t pos []= {0,0,0};
    AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
    AddCluster(phos);
  }
  
  //Fill AOD with EMCAL clusters
  cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
  for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
    TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
    Float_t pos []= {0,0,0};
    AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
    AddCluster(emcal);
  }
}
Commit	Line	Data
bdcfac30	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$
60fecfb6	20	* Revision 1.4 2007/10/29 13:48:42 gustavo
	21	* Corrected coding violations
	22	*
4b707925	23	* Revision 1.2 2007/08/17 12:40:04 schutz
	24	* New analysis classes by Gustavo Conesa
	25	*
bdcfac30	26	* Revision 1.1.2.1 2007/07/26 10:32:09 schutz
	27	* new analysis classes in the the new analysis framework
	28	*
	29	*
	30	*/
	31
	32	//_________________________________________________________________________
3bb2c538	33	// Base class for gamma and correlation analysis
	34	// It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
	35	// got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes
	36	// AliAnaGammaDirect, AliAnaGammaCorrelation ....
	37	//
bdcfac30	38	//*-- Author: Gustavo Conesa (INFN-LNF)
	39
	40	// --- ROOT system ---
	41
	42	#include <TParticle.h>
	43	#include <TH2.h>
	44
	45	//---- AliRoot system ----
	46	#include "AliAnaGamma.h"
	47	#include "AliGammaReader.h"
	48	#include "AliAnaGammaDirect.h"
	49	#include "AliAnaGammaCorrelation.h"
3bb2c538	50	#include "AliAnaGammaSelection.h"
bdcfac30	51	#include "AliNeutralMesonSelection.h"
4b707925	52	#include "AliAODCaloCluster.h"
4b707925	53	#include "AliAODEvent.h"
bdcfac30	54	#include "Riostream.h"
	55	#include "AliLog.h"
	56
	57	ClassImp(AliAnaGamma)
	58
	59
	60	//____________________________________________________________________________
	61	AliAnaGamma::AliAnaGamma() :
	62	TObject(),
	63	fOutputContainer(0x0),
	64	fAnaType(0), fCalorimeter(0), fData(0x0), fKine(0x0),
3bb2c538	65	fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
4b707925	66	fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
bdcfac30	67	{
	68	//Default Ctor
	69
	70	//Initialize parameters, pointers and histograms
	71	if(!fReader)
	72	fReader = new AliGammaReader();
	73	if(!fGammaDirect)
	74	fGammaDirect = new AliAnaGammaDirect();
	75	if(!fGammaCorrelation)
	76	fGammaCorrelation = new AliAnaGammaCorrelation();
3bb2c538	77	if(!fGammaSelection)
3bb2c538	78	fGammaSelection = new AliAnaGammaSelection();
bdcfac30	79	if(!fNeutralMesonSelection)
	80	fNeutralMesonSelection = new AliNeutralMesonSelection();
	81
4b707925	82	fAODclusters = 0;
4b707925	83
bdcfac30	84	InitParameters();
	85
	86	}
	87
	88	//____________________________________________________________________________
	89	AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :
	90	TObject(),
	91	fOutputContainer(g. fOutputContainer),
	92	fAnaType(g.fAnaType), fCalorimeter(g.fCalorimeter),
	93	fData(g.fData), fKine(g.fKine),fReader(g.fReader),
	94	fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
3bb2c538	95	fGammaSelection(g.fGammaSelection),
4b707925	96	fNeutralMesonSelection(g.fNeutralMesonSelection),
4b707925	97	fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
bdcfac30	98	{
	99	// cpy ctor
	100
	101	}
	102
	103	//_________________________________________________________________________
	104	AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
	105	{
	106	// assignment operator
	107
	108	if(this == &source)return *this;
	109	((TObject *)this)->operator=(source);
	110
	111	fOutputContainer = source.fOutputContainer ;
	112	fAnaType = source.fAnaType;
	113	fCalorimeter = source.fCalorimeter ;
	114	fData = source.fData ;
	115	fKine = source.fKine ;
	116	fReader = source.fReader ;
	117	fGammaDirect = source.fGammaDirect ;
	118	fGammaCorrelation = source.fGammaCorrelation ;
3bb2c538	119	fGammaSelection = source.fGammaSelection ;
bdcfac30	120	fNeutralMesonSelection = source.fNeutralMesonSelection ;
	121
	122	return *this;
	123
	124	}
	125
	126	//____________________________________________________________________________
	127	AliAnaGamma::~AliAnaGamma()
	128	{
60fecfb6	129	// Remove all pointers.
	130
	131	// Protection added in case of NULL pointers (MG)
	132	if (fOutputContainer) {
	133	fOutputContainer->Clear();
	134	delete fOutputContainer ;
	135	}
	136
	137	if (fData) delete fData ;
	138	if (fKine) delete fKine ;
	139	if (fReader) delete fReader ;
	140	if (fGammaDirect) delete fGammaDirect ;
	141	if (fGammaCorrelation) delete fGammaCorrelation ;
	142	if (fGammaSelection) delete fGammaSelection ;
	143	if (fNeutralMesonSelection) delete fNeutralMesonSelection ;
bdcfac30	144
	145	}
	146
	147	//________________________________________________________________________
	148	void AliAnaGamma::Init()
	149	{
	150
	151	//Init container histograms and other common variables
	152
	153	//Histograms container
	154	fOutputContainer = new TList ;
	155
	156	//Fill container with appropriate histograms
	157
3bb2c538	158	//Set selection analysis histograms
	159	TList * selectcontainer = fGammaSelection->GetCreateOutputObjects();
	160	for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
	161	Bool_t add = kTRUE ;
	162	TString name = (selectcontainer->At(i))->GetName();
	163	if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
	164	if(!fReader->IsPHOSOn() && name.Contains("PHOS")) add = kFALSE;
	165	if(!fReader->IsCTSOn() && !fGammaSelection->FillCTS() && name.Contains("CTS")) add = kFALSE;
	166	if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
	167	} //Set selection analysis histograms
	168
	169
bdcfac30	170	//Set prompt photon analysis histograms
	171	TList * promptcontainer = fGammaDirect->GetCreateOutputObjects();
	172	for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
	173	fOutputContainer->Add(promptcontainer->At(i)) ;
	174
4b707925	175	//Check if selected options are correct or set them when necessary
3bb2c538	176	if(fReader->GetDataType() == AliGammaReader::kMCData){
4b707925	177	fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
3bb2c538	178	fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
3bb2c538	179	}
4b707925	180
bdcfac30	181	if(fAnaType == kCorrelation){
	182
	183	//Check if selected options are correct
	184	if (fGammaDirect->GetICMethod()==AliAnaGammaDirect::kSeveralIC)
	185	AliFatal("Correlation not allowed with multiple isolation cuts, kCorrelation and kSeveralIC do not go together");
	186
	187	if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kParton &&
	188	fReader->GetDataType() != AliGammaReader::kMC)
	189	AliFatal("kParton must be analyzed with data kMC");
	190
	191	//Set the parameters for the neutral pair selection depending on the analysis,
	192	fNeutralMesonSelection->SetDeltaPhiCutRange(fGammaCorrelation->GetDeltaPhiMinCut(),
	193	fGammaCorrelation->GetDeltaPhiMaxCut());
	194
	195
	196	if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kHadron){
	197	fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiMinPt);
	198	fNeutralMesonSelection->SetMinPt(fGammaCorrelation->GetMinPtHadron());
	199
	200	}
	201
	202	if(fGammaCorrelation->GetCorrelationType() == AliAnaGammaCorrelation::kJetLeadCone){
	203	fNeutralMesonSelection->SetPhiPtSelection(AliNeutralMesonSelection::kSelectPhiPtRatio);
	204	fNeutralMesonSelection->SetRatioCutRange(fGammaCorrelation->GetRatioMinCut(),
	205	fGammaCorrelation->GetRatioMaxCut());
	206	}
	207
	208	//Set the neutral mesosn selection histograms
	209	TList * neutralmesoncontainer = fNeutralMesonSelection->GetCreateOutputObjects();
	210	if(fNeutralMesonSelection->AreNeutralMesonSelectionHistosKept()){
	211	for(Int_t i = 0; i < neutralmesoncontainer->GetEntries(); i++)
	212	fOutputContainer->Add(neutralmesoncontainer->At(i)) ;
	213	}
	214
	215	//Set correlation histograms
	216	TList * correlationcontainer = fGammaCorrelation->GetCreateOutputObjects();
	217	for(Int_t i = 0; i < correlationcontainer->GetEntries(); i++)
	218	fOutputContainer->Add(correlationcontainer->At(i)) ;
	219	fGammaCorrelation->SetOutputContainer(fOutputContainer);
	220	fGammaCorrelation->SetNeutralMesonSelection(fNeutralMesonSelection);
	221
	222	}//kCorrelation
	223
	224	}
	225
	226	//____________________________________________________________________________
	227	void AliAnaGamma::InitParameters()
	228	{
	229
	230	//Init data members
	231
	232	fAnaType = kPrompt;
	233	fCalorimeter = "EMCAL";
	234
	235	}
	236
	237	//__________________________________________________________________
	238	void AliAnaGamma::Print(const Option_t * opt) const
	239	{
	240
	241	//Print some relevant parameters set for the analysis
	242	if(! opt)
	243	return;
	244
245	Info("Print", "%s %s", GetName(), GetTitle() ) ;
246	printf("Analysis type = %d\n", fAnaType) ;
247	printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
248
249	switch(fAnaType)
250	{
251	case kPrompt:
252	{
253	fGammaDirect->Print("");
254	}// case kIsolationCut
255	break;
256
257	case kCorrelation:
258	{
259	fGammaCorrelation->Print("");
260	}// case kCorrelation
261	break;
262
263	}//switch
264	}
265
266
267	//____________________________________________________________________________
268	Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
269
270	AliDebug(1,Form("Entry %d",entry));
3bb2c538	271	//cout<<"Event >>>>>>>>>>>>> "<<entry<<endl;
bdcfac30	272	if(!fOutputContainer)
	273	AliFatal("Histograms not initialized");
	274
	275	//CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
	276
	277	TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
	278	TClonesArray * plEMCAL = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter (EMCAL)
	279	TClonesArray * plPHOS = new TClonesArray("TParticle",1000); // All particles measured Gamma calorimeter
4b707925	280	TClonesArray * plPrimCTS = new TClonesArray("TParticle",1000); // primary tracks
	281	TClonesArray * plPrimEMCAL = new TClonesArray("TParticle",1000); // primary emcal clusters
	282	TClonesArray * plPrimPHOS = new TClonesArray("TParticle",1000); // primary phos clusters
bdcfac30	283	TClonesArray * plParton = new TClonesArray("TParticle",1000); // All partons
	284	//Fill lists with photons, neutral particles and charged particles
	285	//look for the highest energy photon in the event inside fCalorimeter
	286
	287	//Fill particle lists
	288	if(fReader->GetDataType() == AliGammaReader::kData){
	289	AliDebug(1,"Data analysis");
4b707925	290	fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL);
bdcfac30	291	}
	292	else if( fReader->GetDataType()== AliGammaReader::kMC){
	293	AliDebug(1,"Kinematics analysis");
4b707925	294	fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL);
bdcfac30	295	}
	296	else if(fReader->GetDataType() == AliGammaReader::kMCData) {
	297	AliDebug(1,"Data + Kinematics analysis");
4b707925	298	fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS);
bdcfac30	299	}
	300	else
	301	AliError("Option not implemented");
4b707925	302
	303	//Fill AOD with calorimeter
	304	//Temporal solution, just for testing
	305	//FillAODs(plPHOS,plEMCAL);
	306
3bb2c538	307	//Select particles to do the final analysis.
	308	if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
	309	if(fReader->IsPHOSOn()) fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
	310	if(fReader->IsCTSOn() && fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
	311
	312
bdcfac30	313	//Search highest energy prompt gamma in calorimeter
bdcfac30	314	if(fCalorimeter == "PHOS")
4b707925	315	MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ;
bdcfac30	316	else if (fCalorimeter == "EMCAL")
4b707925	317	MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ;
bdcfac30	318	else
	319	AliFatal("Wrong calorimeter name");
	320
	321	plCTS->Clear() ;
	322	plEMCAL->Clear() ;
	323	plPHOS->Clear() ;
	324	plParton->Clear() ;
4b707925	325	plPrimCTS->Clear() ;
	326	plPrimEMCAL->Clear() ;
	327	plPrimPHOS->Clear() ;
bdcfac30	328
	329	delete plCTS ;
	330	delete plPHOS ;
	331	delete plEMCAL ;
	332	delete plParton ;
4b707925	333	delete plPrimCTS ;
	334	delete plPrimPHOS ;
	335	delete plPrimEMCAL ;
bdcfac30	336
	337	return kTRUE;
	338
	339	}
	340
	341	//____________________________________________________________________________
4b707925	342	void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo) {
bdcfac30	343
	344	TParticle * pGamma = new TParticle ;
	345	Bool_t isInCalo = kFALSE ;
	346
	347	switch(fAnaType)
	348	{
	349
	350	//Only Prompt photon analysis
	351	case kPrompt:
	352	{
	353	AliDebug(1,"kPrompt analysis");
	354	switch(fGammaDirect->GetICMethod())
	355	{
	356
	357	case AliAnaGammaDirect::kSeveralIC:
	358	{
	359	fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
	360	AliDebug(1,"kSeveralIC analysis");
	361	}
	362	break;
	363
	364	default :
	365	{
4b707925	366	fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
bdcfac30	367	if(!isInCalo)
	368	AliDebug(1,"Prompt gamma not found");
	369	}
	370	break;
	371	}//IC method
	372	}// case kPrompt:
	373	break;
	374
	375	//Correlate prompt photon with something: parton, hadron, jet.
	376	case kCorrelation:
	377	{
	378	AliDebug(1,"kCorrelation analysis");
	379	//Find prompt photon
4b707925	380	fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
bdcfac30	381
	382	if(isInCalo){//If prompt photon found, do correlation
	383
	384	switch(fGammaCorrelation->GetCorrelationType())
	385	{
	386	case AliAnaGammaCorrelation::kParton:
	387	{
	388	AliDebug(1,"kParton correlation");
	389	fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
	390	}// case kParton
	391	break;
	392
	393	case AliAnaGammaCorrelation::kHadron:
	394	{
	395	AliDebug(1,"kHadron correlation");
	396	fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	397	}// case kHadron
	398	break;
	399
	400	case AliAnaGammaCorrelation::kJetLeadCone:
	401	{
	402	AliDebug(1,"kJetLeadCone correlation");
	403	fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	404	}// case kJetLeadCone
	405	break;
	406
	407	case AliAnaGammaCorrelation::kJetFinder:
	408	{
	409	AliDebug(1,"kJetFinder correlation");
	410	printf("Analysis not implemented \n");
	411	}// case kJetFinder
	412	break;
	413	}// switch correlation
	414	}// is in calo
	415	else AliDebug(2,"Prompt gamma not found");
	416	}// case kCorrelation
	417	break;
	418
	419	} //switch(fAnaType)
	420
	421	delete pGamma ;
	422
	423	}
4b707925	424
	425	//____________________________________________________
	426	void AliAnaGamma::AddCluster(AliAODCaloCluster p)
	427	{
	428	// Add new jet to the list
	429	cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
	430	cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
	431	new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
	432	}
	433
	434	//___________________________________________________
	435	void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
	436	{
	437	// Connect to the AOD
	438	fAODclusters = aod->GetCaloClusters();
	439	}
	440
	441	//____________________________________________________________________________
	442	void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
	443	//Fill AOD caloClusters
	444	//Temporal method, just for testing AOD creation
	445
	446	//Fill AOD with PHOS clusters
	447	Int_t nphos = plPHOS->GetEntries() ;
	448	cout<<"PHOS entries "<<nphos<<endl;
	449	for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
	450	TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
	451	Float_t pos []= {0,0,0};
	452	AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
	453	AddCluster(phos);
	454	}
	455
	456	//Fill AOD with EMCAL clusters
	457	cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
	458	for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
	459	TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
	460	Float_t pos []= {0,0,0};
	461	AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
	462	AddCluster(emcal);
	463	}
	464	}