[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(){
  //Process analysis for this event

  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	//____________________________________________________________________________
17a1db27	268	Bool_t AliAnaGamma::ProcessEvent(){
17a1db27	269	//Process analysis for this event
bdcfac30	270
bdcfac30	271	if(!fOutputContainer)
	272	AliFatal("Histograms not initialized");
	273
	274	//CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
	275
	276	TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
	277	TClonesArray * plEMCAL = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter (EMCAL)
	278	TClonesArray * plPHOS = new TClonesArray("TParticle",1000); // All particles measured Gamma calorimeter
4b707925	279	TClonesArray * plPrimCTS = new TClonesArray("TParticle",1000); // primary tracks
	280	TClonesArray * plPrimEMCAL = new TClonesArray("TParticle",1000); // primary emcal clusters
	281	TClonesArray * plPrimPHOS = new TClonesArray("TParticle",1000); // primary phos clusters
bdcfac30	282	TClonesArray * plParton = new TClonesArray("TParticle",1000); // All partons
	283	//Fill lists with photons, neutral particles and charged particles
	284	//look for the highest energy photon in the event inside fCalorimeter
	285
	286	//Fill particle lists
	287	if(fReader->GetDataType() == AliGammaReader::kData){
	288	AliDebug(1,"Data analysis");
4b707925	289	fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL);
bdcfac30	290	}
	291	else if( fReader->GetDataType()== AliGammaReader::kMC){
	292	AliDebug(1,"Kinematics analysis");
4b707925	293	fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL);
bdcfac30	294	}
	295	else if(fReader->GetDataType() == AliGammaReader::kMCData) {
	296	AliDebug(1,"Data + Kinematics analysis");
4b707925	297	fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS);
bdcfac30	298	}
	299	else
	300	AliError("Option not implemented");
4b707925	301
	302	//Fill AOD with calorimeter
	303	//Temporal solution, just for testing
	304	//FillAODs(plPHOS,plEMCAL);
	305
3bb2c538	306	//Select particles to do the final analysis.
	307	if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
	308	if(fReader->IsPHOSOn()) fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
	309	if(fReader->IsCTSOn() && fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
	310
	311
bdcfac30	312	//Search highest energy prompt gamma in calorimeter
bdcfac30	313	if(fCalorimeter == "PHOS")
4b707925	314	MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ;
bdcfac30	315	else if (fCalorimeter == "EMCAL")
4b707925	316	MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ;
bdcfac30	317	else
	318	AliFatal("Wrong calorimeter name");
	319
	320	plCTS->Clear() ;
	321	plEMCAL->Clear() ;
	322	plPHOS->Clear() ;
	323	plParton->Clear() ;
4b707925	324	plPrimCTS->Clear() ;
	325	plPrimEMCAL->Clear() ;
	326	plPrimPHOS->Clear() ;
bdcfac30	327
	328	delete plCTS ;
	329	delete plPHOS ;
	330	delete plEMCAL ;
	331	delete plParton ;
4b707925	332	delete plPrimCTS ;
	333	delete plPrimPHOS ;
	334	delete plPrimEMCAL ;
bdcfac30	335
	336	return kTRUE;
	337
	338	}
	339
	340	//____________________________________________________________________________
4b707925	341	void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo) {
bdcfac30	342
	343	TParticle * pGamma = new TParticle ;
	344	Bool_t isInCalo = kFALSE ;
	345
	346	switch(fAnaType)
	347	{
	348
	349	//Only Prompt photon analysis
	350	case kPrompt:
	351	{
	352	AliDebug(1,"kPrompt analysis");
	353	switch(fGammaDirect->GetICMethod())
	354	{
	355
	356	case AliAnaGammaDirect::kSeveralIC:
	357	{
	358	fGammaDirect->MakeSeveralICAnalysis(plCalo, plCTS);
	359	AliDebug(1,"kSeveralIC analysis");
	360	}
	361	break;
	362
	363	default :
	364	{
4b707925	365	fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
bdcfac30	366	if(!isInCalo)
	367	AliDebug(1,"Prompt gamma not found");
	368	}
	369	break;
	370	}//IC method
	371	}// case kPrompt:
	372	break;
	373
	374	//Correlate prompt photon with something: parton, hadron, jet.
	375	case kCorrelation:
	376	{
	377	AliDebug(1,"kCorrelation analysis");
	378	//Find prompt photon
4b707925	379	fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
bdcfac30	380
	381	if(isInCalo){//If prompt photon found, do correlation
	382
	383	switch(fGammaCorrelation->GetCorrelationType())
	384	{
	385	case AliAnaGammaCorrelation::kParton:
	386	{
	387	AliDebug(1,"kParton correlation");
	388	fGammaCorrelation->MakeGammaCorrelation(pGamma, plParton, NULL);
	389	}// case kParton
	390	break;
	391
	392	case AliAnaGammaCorrelation::kHadron:
	393	{
	394	AliDebug(1,"kHadron correlation");
	395	fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	396	}// case kHadron
	397	break;
	398
	399	case AliAnaGammaCorrelation::kJetLeadCone:
	400	{
	401	AliDebug(1,"kJetLeadCone correlation");
	402	fGammaCorrelation->MakeGammaCorrelation(pGamma, plCTS, plNe);
	403	}// case kJetLeadCone
	404	break;
	405
	406	case AliAnaGammaCorrelation::kJetFinder:
	407	{
	408	AliDebug(1,"kJetFinder correlation");
	409	printf("Analysis not implemented \n");
	410	}// case kJetFinder
	411	break;
	412	}// switch correlation
	413	}// is in calo
	414	else AliDebug(2,"Prompt gamma not found");
	415	}// case kCorrelation
	416	break;
	417
	418	} //switch(fAnaType)
	419
	420	delete pGamma ;
	421
	422	}
4b707925	423
	424	//____________________________________________________
	425	void AliAnaGamma::AddCluster(AliAODCaloCluster p)
	426	{
	427	// Add new jet to the list
	428	cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
	429	cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
	430	new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
	431	}
	432
	433	//___________________________________________________
	434	void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
	435	{
	436	// Connect to the AOD
	437	fAODclusters = aod->GetCaloClusters();
	438	}
	439
	440	//____________________________________________________________________________
	441	void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
	442	//Fill AOD caloClusters
	443	//Temporal method, just for testing AOD creation
	444
	445	//Fill AOD with PHOS clusters
	446	Int_t nphos = plPHOS->GetEntries() ;
	447	cout<<"PHOS entries "<<nphos<<endl;
	448	for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
	449	TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
	450	Float_t pos []= {0,0,0};
	451	AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
	452	AddCluster(phos);
	453	}
	454
	455	//Fill AOD with EMCAL clusters
	456	cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
	457	for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
	458	TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
	459	Float_t pos []= {0,0,0};
	460	AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
	461	AddCluster(emcal);
	462	}
	463	}