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

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.3  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
 *
 *
 */

//_________________________________________________________________________
// 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.
// There is a 3rd cut based on the gamma correlation on phi or pt.
//-- Author: Gustavo Conesa (INFN-LNF)

// --- ROOT system ---
#include <TParticle.h>
#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(), fSelect(0), fM(0),
    fInvMassMaxCut(0.), fInvMassMinCut(0.),
    fAngleMaxParam(),  fMinPt(0),
    fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
    fRatioMaxCut(0), fRatioMinCut(0),  fKeepNeutralMesonHistos(0),
    fhAnglePairNoCut(0),  fhAnglePairCorrelationCut(0),
    fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0), 
    fhInvMassPairNoCut(0),   fhInvMassPairCorrelationCut(0), 
    fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0) 
{
  //Default Ctor
  
  //Initialize parameters
  
  // kGammaHadron and kGammaJet 
  fAngleMaxParam.Set(4) ;
  fAngleMaxParam.Reset(0.);
  
  //Initialize parameters
  InitParameters();
}

//____________________________________________________________________________
AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :   
  TObject(),  
  fSelect(g.fSelect), fM(g.fM),
  fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
  fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
  fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut), 
  fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut), 
  fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
  fhAnglePairNoCut(g. fhAnglePairNoCut), 
  fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut), 
  fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut), 
  fhAnglePairAllCut(g. fhAnglePairAllCut), 
  fhInvMassPairNoCut(g.fhInvMassPairNoCut),  
  fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut), 
  fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut), 
  fhInvMassPairAllCut(g.fhInvMassPairAllCut)
{
  // cpy ctor
}

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

  fSelect = source.fSelect ;
  fM = source.fM ;
  fInvMassMaxCut = source.fInvMassMaxCut ; 
  fInvMassMinCut = source.fInvMassMinCut ;
  fAngleMaxParam = source.fAngleMaxParam ;
  fMinPt = source.fMinPt ;
  fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; 
  fDeltaPhiMinCut = source.fDeltaPhiMinCut ; 
  fRatioMaxCut = source.fRatioMaxCut ; 
  fRatioMinCut = source.fRatioMinCut ;
  fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
 
  fhAnglePairNoCut = source. fhAnglePairNoCut ; 
  fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ; 
  fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ; 
  fhAnglePairAllCut = source. fhAnglePairAllCut ; 
  fhInvMassPairNoCut = source.fhInvMassPairNoCut ; 
  fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ; 
  fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ; 
  fhInvMassPairAllCut = source.fhInvMassPairAllCut ; 
  
  return *this;
  
}

//____________________________________________________________________________
AliNeutralMesonSelection::~AliNeutralMesonSelection() 
{
 // Remove all pointers except analysis output pointers.

}


//________________________________________________________________________
TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
{  

  // Create histograms to be saved in output file and 
  // store them in outputContainer
  TList * outputContainer = new TList() ; 
  outputContainer->SetName("MesonDecayHistos") ; 
  
  fhAnglePairNoCut  = new TH2F
    ("AnglePairNoCut",
     "Angle between all #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2); 
  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}}"
     ,200,0,50,200,0,0.2); 
  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}}"
     ,200,0,50,200,0,0.2); 
  fhAnglePairAllCut->SetYTitle("Angle (rad)");
  fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
  
  //
  fhInvMassPairNoCut  = new TH2F
    ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
     120,0,120,360,0,0.5); 
  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}}",
     120,0,120,360,0,0.5); 
  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}}",
     120,0,120,360,0,0.5); 
  fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
  fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");

  fhAnglePairCorrelationCut  = new TH2F
    ("AnglePairCorrelationCut",
     "Angle between correlated #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2); 
  fhAnglePairCorrelationCut->SetYTitle("Angle (rad)");
  fhAnglePairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");

  fhInvMassPairCorrelationCut  = new TH2F
    ("InvMassPairCorrelationCut","Invariant Mass of correlated #gamma pair vs E_{#pi^{0}}",
     120,0,120,360,0,0.5); 
  fhInvMassPairCorrelationCut->SetYTitle("Invariant Mass (GeV/c^{2})");
  fhInvMassPairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
  
  outputContainer->Add(fhAnglePairNoCut) ; 
  outputContainer->Add(fhAnglePairOpeningAngleCut) ;
  outputContainer->Add(fhAnglePairAllCut) ; 
  
  outputContainer->Add(fhInvMassPairNoCut) ; 
  outputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
  outputContainer->Add(fhInvMassPairAllCut) ; 

  outputContainer->Add(fhAnglePairCorrelationCut) ; 
  outputContainer->Add(fhInvMassPairCorrelationCut) ; 
  
  return outputContainer;
}

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

  //-------------kHadron, kJetLeadCone-----------------
  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

  fMinPt = 0.   ;
  fDeltaPhiMaxCut      = 4.5;
  fDeltaPhiMinCut      = 1.5 ;
  fRatioMaxCut    = 1.0 ;
  fRatioMinCut    = 0.1 ; 
}

//__________________________________________________________________________-
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::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi)  const
{ 
  //Select pair if delta
  Bool_t cut = kFALSE ;
 
  if(fSelect == kNoSelectPhiPt) cut = kTRUE ;
  else if((phig-phi) > fDeltaPhiMinCut && ((phig-phi) < fDeltaPhiMaxCut)){
    //Cut on pt
    if((fSelect == kSelectPhiPtRatio && ptg > 0. && pt/ptg  > fRatioMinCut &&  pt/ptg  < fRatioMaxCut) ||
	(fSelect == kSelectPhiMinPt && pt > fMinPt)  )  cut = kTRUE ;
  }
  else cut = kFALSE ;
  
  return cut ;
  
}

//____________________________________________________________________________
Bool_t  AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector gammai, TLorentzVector gammaj)  
{  
  
  //Search for the neutral pion within selection cuts
  
  Double_t ptg = pGamma->Pt();
  Double_t phig = pGamma->Phi() ;
  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 phig-phi meson
  if(CutPtPhi(ptg, phig, pt, phi)){
    
    fhAnglePairCorrelationCut     ->Fill(e,angle);
    fhInvMassPairCorrelationCut->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
  } // cut on pt and phi
  
  
  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("pT meson       >    %f\n", fMinPt) ; 
  printf("Phi gamma-meson      <     %f\n", fDeltaPhiMaxCut) ; 
  printf("Phi gamma-meson      >     %f\n", fDeltaPhiMinCut) ;
  printf("pT meson / pT Gamma             <     %f\n", fRatioMaxCut) ; 
  printf("pT meson / pT Gamma             >     %f\n", fRatioMinCut) ;
  printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);

}
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$
3788def4	20	* Revision 1.3 2007/10/29 13:48:42 gustavo
	21	* Corrected coding violations
	22	*
3bb2c538	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	//_________________________________________________________________________
	33	// Class that contains methods to select candidate pairs to neutral meson
3bb2c538	34	// 2 main selections, invariant mass around pi0 (also any other mass),
	35	// apperture angle to distinguish from combinatorial.
	36	// There is a 3rd cut based on the gamma correlation on phi or pt.
bdcfac30	37	//-- Author: Gustavo Conesa (INFN-LNF)
	38
	39	// --- ROOT system ---
	40	#include <TParticle.h>
	41	#include <TLorentzVector.h>
	42	#include <TH2.h>
	43	#include <TList.h>
3bb2c538	44	#include <TArrayD.h>
3bb2c538	45
bdcfac30	46	//---- AliRoot system ----
bdcfac30	47	#include "AliNeutralMesonSelection.h"
bdcfac30	48	#include "AliLog.h"
	49
	50	ClassImp(AliNeutralMesonSelection)
	51
	52
	53	//____________________________________________________________________________
	54	AliNeutralMesonSelection::AliNeutralMesonSelection() :
	55	TObject(), fSelect(0), fM(0),
	56	fInvMassMaxCut(0.), fInvMassMinCut(0.),
	57	fAngleMaxParam(), fMinPt(0),
	58	fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
	59	fRatioMaxCut(0), fRatioMinCut(0), fKeepNeutralMesonHistos(0),
	60	fhAnglePairNoCut(0), fhAnglePairCorrelationCut(0),
	61	fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0),
	62	fhInvMassPairNoCut(0), fhInvMassPairCorrelationCut(0),
	63	fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0)
	64	{
	65	//Default Ctor
	66
	67	//Initialize parameters
	68
	69	// kGammaHadron and kGammaJet
	70	fAngleMaxParam.Set(4) ;
	71	fAngleMaxParam.Reset(0.);
	72
	73	//Initialize parameters
	74	InitParameters();
	75	}
	76
	77	//____________________________________________________________________________
	78	AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
	79	TObject(),
	80	fSelect(g.fSelect), fM(g.fM),
	81	fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
	82	fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
	83	fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
	84	fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut),
	85	fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
	86	fhAnglePairNoCut(g. fhAnglePairNoCut),
	87	fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut),
	88	fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
	89	fhAnglePairAllCut(g. fhAnglePairAllCut),
	90	fhInvMassPairNoCut(g.fhInvMassPairNoCut),
	91	fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut),
	92	fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
	93	fhInvMassPairAllCut(g.fhInvMassPairAllCut)
	94	{
	95	// cpy ctor
	96	}
	97
	98	//_________________________________________________________________________
	99	AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
	100	{
	101	// assignment operator
	102
	103	if(this == &source)return *this;
	104	((TObject *)this)->operator=(source);
	105
	106	fSelect = source.fSelect ;
	107	fM = source.fM ;
	108	fInvMassMaxCut = source.fInvMassMaxCut ;
	109	fInvMassMinCut = source.fInvMassMinCut ;
	110	fAngleMaxParam = source.fAngleMaxParam ;
	111	fMinPt = source.fMinPt ;
112	fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
113	fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
114	fRatioMaxCut = source.fRatioMaxCut ;
115	fRatioMinCut = source.fRatioMinCut ;
116	fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
117
118	fhAnglePairNoCut = source. fhAnglePairNoCut ;
119	fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ;
120	fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
121	fhAnglePairAllCut = source. fhAnglePairAllCut ;
122	fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
123	fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ;
124	fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
125	fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
126
127	return *this;
128
129	}
130
131	//____________________________________________________________________________
132	AliNeutralMesonSelection::~AliNeutralMesonSelection()
133	{
3788def4	134	// Remove all pointers except analysis output pointers.
3788def4	135
bdcfac30	136	}
	137
	138
	139
	140	//________________________________________________________________________
	141	TList * AliNeutralMesonSelection::GetCreateOutputObjects()
	142	{
	143
	144	// Create histograms to be saved in output file and
	145	// store them in outputContainer
	146	TList * outputContainer = new TList() ;
	147	outputContainer->SetName("MesonDecayHistos") ;
	148
	149	fhAnglePairNoCut = new TH2F
	150	("AnglePairNoCut",
	151	"Angle between all #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2);
	152	fhAnglePairNoCut->SetYTitle("Angle (rad)");
	153	fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	154
	155	fhAnglePairOpeningAngleCut = new TH2F
	156	("AnglePairOpeningAngleCut",
	157	"Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
	158	,200,0,50,200,0,0.2);
	159	fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
	160	fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	161
	162	fhAnglePairAllCut = new TH2F
	163	("AnglePairAllCut",
	164	"Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
	165	,200,0,50,200,0,0.2);
	166	fhAnglePairAllCut->SetYTitle("Angle (rad)");
	167	fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	168
	169	//
	170	fhInvMassPairNoCut = new TH2F
	171	("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
	172	120,0,120,360,0,0.5);
	173	fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	174	fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	175
	176	fhInvMassPairOpeningAngleCut = new TH2F
	177	("InvMassPairOpeningAngleCut",
	178	"Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
	179	120,0,120,360,0,0.5);
	180	fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	181	fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	182
	183	fhInvMassPairAllCut = new TH2F
	184	("InvMassPairAllCut",
	185	"Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
	186	120,0,120,360,0,0.5);
	187	fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	188	fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	189
	190	fhAnglePairCorrelationCut = new TH2F
	191	("AnglePairCorrelationCut",
	192	"Angle between correlated #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2);
	193	fhAnglePairCorrelationCut->SetYTitle("Angle (rad)");
	194	fhAnglePairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	195
	196	fhInvMassPairCorrelationCut = new TH2F
	197	("InvMassPairCorrelationCut","Invariant Mass of correlated #gamma pair vs E_{#pi^{0}}",
	198	120,0,120,360,0,0.5);
	199	fhInvMassPairCorrelationCut->SetYTitle("Invariant Mass (GeV/c^{2})");
200	fhInvMassPairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
201
202	outputContainer->Add(fhAnglePairNoCut) ;
203	outputContainer->Add(fhAnglePairOpeningAngleCut) ;
204	outputContainer->Add(fhAnglePairAllCut) ;
205
206	outputContainer->Add(fhInvMassPairNoCut) ;
207	outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
208	outputContainer->Add(fhInvMassPairAllCut) ;
209
210	outputContainer->Add(fhAnglePairCorrelationCut) ;
211	outputContainer->Add(fhInvMassPairCorrelationCut) ;
212
213	return outputContainer;
214	}
215
216	//____________________________________________________________________________
217	void AliNeutralMesonSelection::InitParameters()
218	{
219
220	//Initialize the parameters of the analysis.
221	fKeepNeutralMesonHistos = kTRUE ;
222
223	//-------------kHadron, kJetLeadCone-----------------
224	fAngleMaxParam.Set(4) ;
225	fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
226	fAngleMaxParam.AddAt(-0.25,1) ;
227	fAngleMaxParam.AddAt(0.025,2) ;
228	fAngleMaxParam.AddAt(-2e-4,3) ;
229
230	fInvMassMaxCut = 0.16 ;
231	fInvMassMinCut = 0.11 ;
232
233	fM = 0.1349766;//neutralMeson mass
234
235	fMinPt = 0. ;
236	fDeltaPhiMaxCut = 4.5;
237	fDeltaPhiMinCut = 1.5 ;
238	fRatioMaxCut = 1.0 ;
239	fRatioMinCut = 0.1 ;
240	}
241
242	//__________________________________________________________________________-
3bb2c538	243	Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
bdcfac30	244	//Check if the opening angle of the candidate pairs is inside
	245	//our selection windowd
	246
	247	Bool_t result = kFALSE;
	248	Double_t max = fAngleMaxParam.At(0)TMath::Exp(fAngleMaxParam.At(1)e)
	249	+fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
	250	Double_t arg = (ee-2fMfM)/(ee);
	251	Double_t min = 100. ;
	252	if(arg>0.)
	253	min = TMath::ACos(arg);
	254
	255	if((angle<max)&&(angle>=min))
	256	result = kTRUE;
	257
	258	return result;
	259	}
	260
	261	//____________________________________________________________________________
3bb2c538	262	Bool_t AliNeutralMesonSelection::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi) const
bdcfac30	263	{
	264	//Select pair if delta
	265	Bool_t cut = kFALSE ;
	266
	267	if(fSelect == kNoSelectPhiPt) cut = kTRUE ;
	268	else if((phig-phi) > fDeltaPhiMinCut && ((phig-phi) < fDeltaPhiMaxCut)){
	269	//Cut on pt
	270	if((fSelect == kSelectPhiPtRatio && ptg > 0. && pt/ptg > fRatioMinCut && pt/ptg < fRatioMaxCut) \|\|
	271	(fSelect == kSelectPhiMinPt && pt > fMinPt) ) cut = kTRUE ;
	272	}
	273	else cut = kFALSE ;
	274
	275	return cut ;
	276
	277	}
	278
	279	//____________________________________________________________________________
	280	Bool_t AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector gammai, TLorentzVector gammaj)
	281	{
	282
	283	//Search for the neutral pion within selection cuts
	284
	285	Double_t ptg = pGamma->Pt();
	286	Double_t phig = pGamma->Phi() ;
	287	Bool_t goodpair = kFALSE ;
	288
	289	Double_t pt = (gammai+gammaj).Pt();
	290	Double_t phi = (gammai+gammaj).Phi();
	291	if(phi < 0)
	292	phi+=TMath::TwoPi();
	293	Double_t invmass = (gammai+gammaj).M();
	294	Double_t angle = gammaj.Angle(gammai.Vect());
	295	Double_t e = (gammai+gammaj).E();
	296
	297	//Fill histograms with no cuts applied.
	298	fhAnglePairNoCut->Fill(e,angle);
	299	fhInvMassPairNoCut->Fill(e,invmass);
	300
	301	//Cut on phig-phi meson
	302	if(CutPtPhi(ptg, phig, pt, phi)){
	303
	304	fhAnglePairCorrelationCut ->Fill(e,angle);
	305	fhInvMassPairCorrelationCut->Fill(e,invmass);
	306
	307	//Cut on the aperture of the pair
	308	if(IsAngleInWindow(angle,e)){
	309	fhAnglePairOpeningAngleCut ->Fill(e,angle);
	310	fhInvMassPairOpeningAngleCut->Fill(e,invmass);
	311	AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
	312
	313	//Cut on the invariant mass of the pair
	314	if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
	315	fhInvMassPairAllCut ->Fill(e,invmass);
	316	fhAnglePairAllCut ->Fill(e,angle);
	317	goodpair = kTRUE;
	318	AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
	319	}//(invmass>0.125) && (invmass<0.145)
	320	}//Opening angle cut
	321	} // cut on pt and phi
	322
	323
	324	return goodpair;
	325
	326	}
327
328	//__________________________________________________________________
329	void AliNeutralMesonSelection::Print(const Option_t * opt) const
330	{
331
332	//Print some relevant parameters set for the analysis
333	if(! opt)
334	return;
335
336	Info("Print", "%s %s", GetName(), GetTitle() ) ;
337
338	printf("mass : %f \n", fM );
339	printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
340	printf("Angle selection param: \n");
341	printf("p0 : %f", fAngleMaxParam.At(0));
342	printf("p1 : %f", fAngleMaxParam.At(1));
343	printf("p2 : %f", fAngleMaxParam.At(2));
344	printf("p3 : %f", fAngleMaxParam.At(3));
345
346	printf("pT meson > %f\n", fMinPt) ;
347	printf("Phi gamma-meson < %f\n", fDeltaPhiMaxCut) ;
348	printf("Phi gamma-meson > %f\n", fDeltaPhiMinCut) ;
349	printf("pT meson / pT Gamma < %f\n", fRatioMaxCut) ;
350	printf("pT meson / pT Gamma > %f\n", fRatioMinCut) ;
351	printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
352
353	}