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