[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliNeutralMesonSelection.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: AliNeutralMesonSelection.cxx 27413 2008-07-18 13:28:12Z gconesab $ */

//_________________________________________________________________________
// Class that contains methods to select candidate pairs to neutral meson 
// 2 main selections, invariant mass around pi0 (also any other mass),
// apperture angle to distinguish from combinatorial.
//-- Author: Gustavo Conesa (INFN-LNF)

// --- ROOT system ---
#include <TLorentzVector.h>
#include <TH2.h>
#include <TList.h>
#include <TArrayD.h>
 
//---- AliRoot system ----
#include "AliNeutralMesonSelection.h" 
#include "AliLog.h"

ClassImp(AliNeutralMesonSelection)
  
  
//____________________________________________________________________________
  AliNeutralMesonSelection::AliNeutralMesonSelection() : 
    TObject(), fM(0),
    fInvMassMaxCut(0.), fInvMassMinCut(0.),
    fAngleMaxParam(),  fKeepNeutralMesonHistos(0), 
    fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0), 
    fhAnglePairAllCut(0), 
    fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0), 
    fhInvMassPairAllCut(0),
	fHistoNEBins(0),   fHistoEMax(0.),   fHistoEMin(0.),
	fHistoNPtBins(0),  fHistoPtMax(0.),  fHistoPtMin(0.),
	fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
	fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.){
  //Default Ctor
  
  //Initialize parameters
  
  // kGammaHadron and kGammaJet 
  fAngleMaxParam.Set(4) ;
  fAngleMaxParam.Reset(0.);
  
  //Initialize parameters
  InitParameters();
}

//____________________________________________________________________________
AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :   
  TObject(), fM(g.fM),
  fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
  fAngleMaxParam(g.fAngleMaxParam),
  fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
  fhAnglePairNoCut(g. fhAnglePairNoCut), 
  fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut), 
  fhAnglePairAllCut(g. fhAnglePairAllCut), 
  fhInvMassPairNoCut(g.fhInvMassPairNoCut),  
  fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut), 
  fhInvMassPairAllCut(g.fhInvMassPairAllCut),
  fHistoNEBins(g.fHistoNEBins),   fHistoEMax(g.fHistoEMax),   fHistoEMin(g.fHistoEMin),
  fHistoNPtBins(g.fHistoNPtBins),   fHistoPtMax(g.fHistoPtMax),   fHistoPtMin(g.fHistoPtMin),
  fHistoNAngleBins(g.fHistoNAngleBins), fHistoAngleMax(g.fHistoAngleMax), fHistoAngleMin(g.fHistoAngleMin),
  fHistoNIMBins(g.fHistoNIMBins), fHistoIMMax(g.fHistoIMMax), fHistoIMMin(g.fHistoIMMin)
{
  // cpy ctor
}

//_________________________________________________________________________
AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
{
  // assignment operator
  
  if(this == &source)return *this;
  ((TObject *)this)->operator=(source);

  fM = source.fM ;
  fInvMassMaxCut = source.fInvMassMaxCut ; 
  fInvMassMinCut = source.fInvMassMinCut ;
  fAngleMaxParam = source.fAngleMaxParam ;
  fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
 
  fhAnglePairNoCut = source. fhAnglePairNoCut ; 
  fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ; 
  fhAnglePairAllCut = source. fhAnglePairAllCut ; 
  fhInvMassPairNoCut = source.fhInvMassPairNoCut ; 
  fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ; 
  fhInvMassPairAllCut = source.fhInvMassPairAllCut ; 
  
  fHistoNEBins = source.fHistoNEBins;   fHistoEMax = source.fHistoEMax;   fHistoEMin = source.fHistoEMin;
  fHistoNPtBins = source.fHistoNPtBins;   fHistoPtMax = source.fHistoPtMax;   fHistoPtMin = source.fHistoPtMin;
  fHistoNAngleBins = source.fHistoNAngleBins; fHistoAngleMax = source.fHistoAngleMax; fHistoAngleMin = source.fHistoAngleMin;
  fHistoNIMBins = source.fHistoNIMBins; fHistoIMMax = source.fHistoIMMax; fHistoIMMin = source.fHistoIMMin;
  
  return *this;
  
}

//____________________________________________________________________________
AliNeutralMesonSelection::~AliNeutralMesonSelection() 
{
  //dtor

  if(!fKeepNeutralMesonHistos){
    //Histograms initialized and filled but not passed to output container
    //delete here, I am not sure this is correct
    
    if(fhAnglePairNoCut) delete fhAnglePairNoCut;
    if(fhAnglePairOpeningAngleCut) delete fhAnglePairOpeningAngleCut; 
    if(fhAnglePairAllCut) delete fhAnglePairAllCut;
    if(fhInvMassPairNoCut) delete fhInvMassPairNoCut;
    if(fhInvMassPairOpeningAngleCut) delete fhInvMassPairOpeningAngleCut;
    if(fhInvMassPairAllCut) delete fhInvMassPairAllCut; 

  }
  
}
//________________________________________________________________________
TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
{  
	
	// Create histograms to be saved in output file and 
	// store them in outputContainer of the analysis class that calls this class.
	
	TList * outputContainer = new TList() ; 
	outputContainer->SetName("MesonDecayHistos") ; 
			  
	fhAnglePairNoCut  = new TH2F
    ("AnglePairNoCut",
     "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
	fhAnglePairNoCut->SetYTitle("Angle (rad)");
	fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	
	fhAnglePairOpeningAngleCut  = new TH2F
    ("AnglePairOpeningAngleCut",
     "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
     ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
	fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
	fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	
	fhAnglePairAllCut  = new TH2F
    ("AnglePairAllCut",
     "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
     ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
	fhAnglePairAllCut->SetYTitle("Angle (rad)");
	fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
	
	//
	fhInvMassPairNoCut  = new TH2F
    ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
     fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
	fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	
	fhInvMassPairOpeningAngleCut  = new TH2F
    ("InvMassPairOpeningAngleCut",
     "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
     fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
	fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	
	fhInvMassPairAllCut  = new TH2F
    ("InvMassPairAllCut",
     "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
     fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
	fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	
	outputContainer->Add(fhAnglePairNoCut) ; 
	outputContainer->Add(fhAnglePairOpeningAngleCut) ;
	outputContainer->Add(fhAnglePairAllCut) ; 
	
	outputContainer->Add(fhInvMassPairNoCut) ; 
	outputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
	outputContainer->Add(fhInvMassPairAllCut) ; 
	
	return outputContainer;
}

 //____________________________________________________________________________
void AliNeutralMesonSelection::InitParameters()
{
 
  //Initialize the parameters of the analysis.
  fKeepNeutralMesonHistos = kFALSE ;

  fAngleMaxParam.Set(4) ;
  fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
  fAngleMaxParam.AddAt(-0.25,1) ;
  fAngleMaxParam.AddAt(0.025,2) ;
  fAngleMaxParam.AddAt(-2e-4,3) ;

  fInvMassMaxCut  = 0.16 ;
  fInvMassMinCut  = 0.11 ;

  fM = 0.1349766;//neutralMeson mass, pi0
  
 //Histogrammes settings
  fHistoNEBins = 100 ;
  fHistoEMax   = 50 ;
  fHistoEMin   = 0.  ;  
  
  fHistoNPtBins = 240 ;
  fHistoPtMax   = 120 ;
  fHistoPtMin   = 0.  ;

  fHistoNAngleBins = 200 ;
  fHistoAngleMax   = 0.2;
  fHistoAngleMin   = 0.  ;

  fHistoNIMBins = 300 ;
  fHistoIMMax   = 0.5   ;
  fHistoIMMin   = 0.  ;  
}

//__________________________________________________________________________-
Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
  //Check if the opening angle of the candidate pairs is inside 
  //our selection windowd

  Bool_t result = kFALSE;
  Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
    +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
  Double_t arg = (e*e-2*fM*fM)/(e*e);
  Double_t min = 100. ;
  if(arg>0.)
    min = TMath::ACos(arg);

  if((angle<max)&&(angle>=min))
    result = kTRUE;
 
  return result;
}

//____________________________________________________________________________
Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)  
{  
  
  //Search for the neutral pion within selection cuts
  Bool_t goodpair = kFALSE ;
  
  Double_t pt  = (gammai+gammaj).Pt();
  Double_t phi = (gammai+gammaj).Phi();
  if(phi < 0)
    phi+=TMath::TwoPi();
  Double_t invmass = (gammai+gammaj).M();
  Double_t angle   = gammaj.Angle(gammai.Vect());
  Double_t e       = (gammai+gammaj).E();
  
  //Fill histograms with no cuts applied.
  fhAnglePairNoCut->Fill(e,angle);
  fhInvMassPairNoCut->Fill(e,invmass);
    
  //Cut on the aperture of the pair
  if(IsAngleInWindow(angle,e)){
    fhAnglePairOpeningAngleCut     ->Fill(e,angle);
    fhInvMassPairOpeningAngleCut->Fill(e,invmass);
    AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
    
    //Cut on the invariant mass of the pair
    if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
      fhInvMassPairAllCut  ->Fill(e,invmass);
      fhAnglePairAllCut       ->Fill(e,angle);
      goodpair = kTRUE;
      AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
    }//(invmass>0.125) && (invmass<0.145)
  }//Opening angle cut
  
  return goodpair; 
  
}

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

  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  Info("Print", "%s %s", GetName(), GetTitle() ) ;

  printf("mass : %f  \n", fM );
  printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
  printf("Angle selection param: \n");
  printf("p0 :     %f", fAngleMaxParam.At(0));
  printf("p1 :     %f", fAngleMaxParam.At(1));
  printf("p2 :     %f", fAngleMaxParam.At(2));
  printf("p3 :     %f", fAngleMaxParam.At(3));

  printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
  
  if(fKeepNeutralMesonHistos){
	printf("Histograms: %3.1f < E  < %3.1f,  Nbin = %d\n", fHistoEMin,  fHistoEMax,  fHistoNEBins);
  	printf("Histograms: %3.1f < pT < %3.1f,  Nbin = %d\n", fHistoPtMin,  fHistoPtMax,  fHistoNPtBins);
	printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
	printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
  
  }

}
Commit	Line	Data
1c5acb87	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: AliNeutralMesonSelection.cxx 27413 2008-07-18 13:28:12Z gconesab $ */
	16
	17	//_________________________________________________________________________
	18	// Class that contains methods to select candidate pairs to neutral meson
	19	// 2 main selections, invariant mass around pi0 (also any other mass),
	20	// apperture angle to distinguish from combinatorial.
	21	//-- Author: Gustavo Conesa (INFN-LNF)
	22
	23	// --- ROOT system ---
	24	#include <TLorentzVector.h>
	25	#include <TH2.h>
	26	#include <TList.h>
	27	#include <TArrayD.h>
	28
	29	//---- AliRoot system ----
	30	#include "AliNeutralMesonSelection.h"
	31	#include "AliLog.h"
	32
	33	ClassImp(AliNeutralMesonSelection)
	34
	35
	36	//____________________________________________________________________________
	37	AliNeutralMesonSelection::AliNeutralMesonSelection() :
	38	TObject(), fM(0),
	39	fInvMassMaxCut(0.), fInvMassMinCut(0.),
	40	fAngleMaxParam(), fKeepNeutralMesonHistos(0),
	41	fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
	42	fhAnglePairAllCut(0),
	43	fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
	44	fhInvMassPairAllCut(0),
	45	fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
	46	fHistoNPtBins(0), fHistoPtMax(0.), fHistoPtMin(0.),
	47	fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
	48	fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.){
	49	//Default Ctor
	50
	51	//Initialize parameters
	52
	53	// kGammaHadron and kGammaJet
	54	fAngleMaxParam.Set(4) ;
	55	fAngleMaxParam.Reset(0.);
	56
	57	//Initialize parameters
	58	InitParameters();
	59	}
	60
	61	//____________________________________________________________________________
	62	AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
	63	TObject(), fM(g.fM),
	64	fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
65	fAngleMaxParam(g.fAngleMaxParam),
66	fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
67	fhAnglePairNoCut(g. fhAnglePairNoCut),
68	fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
69	fhAnglePairAllCut(g. fhAnglePairAllCut),
70	fhInvMassPairNoCut(g.fhInvMassPairNoCut),
71	fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
72	fhInvMassPairAllCut(g.fhInvMassPairAllCut),
73	fHistoNEBins(g.fHistoNEBins), fHistoEMax(g.fHistoEMax), fHistoEMin(g.fHistoEMin),
74	fHistoNPtBins(g.fHistoNPtBins), fHistoPtMax(g.fHistoPtMax), fHistoPtMin(g.fHistoPtMin),
75	fHistoNAngleBins(g.fHistoNAngleBins), fHistoAngleMax(g.fHistoAngleMax), fHistoAngleMin(g.fHistoAngleMin),
76	fHistoNIMBins(g.fHistoNIMBins), fHistoIMMax(g.fHistoIMMax), fHistoIMMin(g.fHistoIMMin)
77	{
78	// cpy ctor
79	}
80
81	//_________________________________________________________________________
82	AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
83	{
84	// assignment operator
85
86	if(this == &source)return *this;
87	((TObject *)this)->operator=(source);
88
89	fM = source.fM ;
90	fInvMassMaxCut = source.fInvMassMaxCut ;
91	fInvMassMinCut = source.fInvMassMinCut ;
92	fAngleMaxParam = source.fAngleMaxParam ;
93	fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
94
95	fhAnglePairNoCut = source. fhAnglePairNoCut ;
96	fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
97	fhAnglePairAllCut = source. fhAnglePairAllCut ;
98	fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
99	fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
100	fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
101
102	fHistoNEBins = source.fHistoNEBins; fHistoEMax = source.fHistoEMax; fHistoEMin = source.fHistoEMin;
103	fHistoNPtBins = source.fHistoNPtBins; fHistoPtMax = source.fHistoPtMax; fHistoPtMin = source.fHistoPtMin;
104	fHistoNAngleBins = source.fHistoNAngleBins; fHistoAngleMax = source.fHistoAngleMax; fHistoAngleMin = source.fHistoAngleMin;
105	fHistoNIMBins = source.fHistoNIMBins; fHistoIMMax = source.fHistoIMMax; fHistoIMMin = source.fHistoIMMin;
106
107	return *this;
108
109	}
110
111	//____________________________________________________________________________
112	AliNeutralMesonSelection::~AliNeutralMesonSelection()
113	{
114	//dtor
115
116	if(!fKeepNeutralMesonHistos){
117	//Histograms initialized and filled but not passed to output container
118	//delete here, I am not sure this is correct
119
120	if(fhAnglePairNoCut) delete fhAnglePairNoCut;
121	if(fhAnglePairOpeningAngleCut) delete fhAnglePairOpeningAngleCut;
122	if(fhAnglePairAllCut) delete fhAnglePairAllCut;
123	if(fhInvMassPairNoCut) delete fhInvMassPairNoCut;
124	if(fhInvMassPairOpeningAngleCut) delete fhInvMassPairOpeningAngleCut;
125	if(fhInvMassPairAllCut) delete fhInvMassPairAllCut;
126
127	}
128
129	}
130	//________________________________________________________________________
131	TList * AliNeutralMesonSelection::GetCreateOutputObjects()
132	{
133
134	// Create histograms to be saved in output file and
135	// store them in outputContainer of the analysis class that calls this class.
136
137	TList * outputContainer = new TList() ;
138	outputContainer->SetName("MesonDecayHistos") ;
139
140	fhAnglePairNoCut = new TH2F
141	("AnglePairNoCut",
142	"Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
143	fhAnglePairNoCut->SetYTitle("Angle (rad)");
144	fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
145
146	fhAnglePairOpeningAngleCut = new TH2F
147	("AnglePairOpeningAngleCut",
148	"Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
149	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
150	fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
151	fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
152
153	fhAnglePairAllCut = new TH2F
154	("AnglePairAllCut",
155	"Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
156	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
157	fhAnglePairAllCut->SetYTitle("Angle (rad)");
158	fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
159
160	//
161	fhInvMassPairNoCut = new TH2F
162	("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
163	fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
164	fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
165	fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
166
167	fhInvMassPairOpeningAngleCut = new TH2F
168	("InvMassPairOpeningAngleCut",
169	"Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
170	fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
171	fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
172	fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
173
174	fhInvMassPairAllCut = new TH2F
175	("InvMassPairAllCut",
176	"Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
177	fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
178	fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
179	fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
180
181	outputContainer->Add(fhAnglePairNoCut) ;
182	outputContainer->Add(fhAnglePairOpeningAngleCut) ;
183	outputContainer->Add(fhAnglePairAllCut) ;
184
185	outputContainer->Add(fhInvMassPairNoCut) ;
186	outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
187	outputContainer->Add(fhInvMassPairAllCut) ;
188
189	return outputContainer;
190	}
191
192	//____________________________________________________________________________
193	void AliNeutralMesonSelection::InitParameters()
194	{
195
196	//Initialize the parameters of the analysis.
197	fKeepNeutralMesonHistos = kFALSE ;
198
199	fAngleMaxParam.Set(4) ;
200	fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
201	fAngleMaxParam.AddAt(-0.25,1) ;
202	fAngleMaxParam.AddAt(0.025,2) ;
203	fAngleMaxParam.AddAt(-2e-4,3) ;
204
205	fInvMassMaxCut = 0.16 ;
206	fInvMassMinCut = 0.11 ;
207
208	fM = 0.1349766;//neutralMeson mass, pi0
209
210	//Histogrammes settings
211	fHistoNEBins = 100 ;
212	fHistoEMax = 50 ;
213	fHistoEMin = 0. ;
214
215	fHistoNPtBins = 240 ;
216	fHistoPtMax = 120 ;
217	fHistoPtMin = 0. ;
218
219	fHistoNAngleBins = 200 ;
220	fHistoAngleMax = 0.2;
221	fHistoAngleMin = 0. ;
222
223	fHistoNIMBins = 300 ;
224	fHistoIMMax = 0.5 ;
225	fHistoIMMin = 0. ;
226	}
227
228	//__________________________________________________________________________-
229	Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
230	//Check if the opening angle of the candidate pairs is inside
231	//our selection windowd
232
233	Bool_t result = kFALSE;
234	Double_t max = fAngleMaxParam.At(0)TMath::Exp(fAngleMaxParam.At(1)e)
235	+fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
236	Double_t arg = (ee-2fMfM)/(ee);
237	Double_t min = 100. ;
238	if(arg>0.)
239	min = TMath::ACos(arg);
240
241	if((angle<max)&&(angle>=min))
242	result = kTRUE;
243
244	return result;
245	}
246
247	//____________________________________________________________________________
248	Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)
249	{
250
251	//Search for the neutral pion within selection cuts
252	Bool_t goodpair = kFALSE ;
253
254	Double_t pt = (gammai+gammaj).Pt();
255	Double_t phi = (gammai+gammaj).Phi();
256	if(phi < 0)
257	phi+=TMath::TwoPi();
258	Double_t invmass = (gammai+gammaj).M();
259	Double_t angle = gammaj.Angle(gammai.Vect());
260	Double_t e = (gammai+gammaj).E();
261
262	//Fill histograms with no cuts applied.
263	fhAnglePairNoCut->Fill(e,angle);
264	fhInvMassPairNoCut->Fill(e,invmass);
265
266	//Cut on the aperture of the pair
267	if(IsAngleInWindow(angle,e)){
268	fhAnglePairOpeningAngleCut ->Fill(e,angle);
269	fhInvMassPairOpeningAngleCut->Fill(e,invmass);
270	AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
271
272	//Cut on the invariant mass of the pair
273	if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
274	fhInvMassPairAllCut ->Fill(e,invmass);
275	fhAnglePairAllCut ->Fill(e,angle);
276	goodpair = kTRUE;
277	AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
278	}//(invmass>0.125) && (invmass<0.145)
279	}//Opening angle cut
280
281	return goodpair;
282
283	}
284
285	//__________________________________________________________________
286	void AliNeutralMesonSelection::Print(const Option_t * opt) const
287	{
288
289	//Print some relevant parameters set for the analysis
290	if(! opt)
291	return;
292
293	Info("Print", "%s %s", GetName(), GetTitle() ) ;
294
295	printf("mass : %f \n", fM );
296	printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
297	printf("Angle selection param: \n");
298	printf("p0 : %f", fAngleMaxParam.At(0));
299	printf("p1 : %f", fAngleMaxParam.At(1));
300	printf("p2 : %f", fAngleMaxParam.At(2));
301	printf("p3 : %f", fAngleMaxParam.At(3));
302
303	printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
304
305	if(fKeepNeutralMesonHistos){
306	printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
307	printf("Histograms: %3.1f < pT < %3.1f, Nbin = %d\n", fHistoPtMin, fHistoPtMax, fHistoNPtBins);
308	printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
309	printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
310
311	}
312
313	}