[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.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 
//-- Author: Gustavo Conesa (INFN-LNF)

// --- ROOT system ---
#include <TParticle.h>
#include <TLorentzVector.h>
#include <TH2.h>
#include <TList.h>

//---- AliRoot system ----
#include "AliNeutralMesonSelection.h" 
#include "Riostream.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
    
  delete fhAnglePairNoCut  ; 
  delete fhAnglePairCorrelationCut  ; 
  delete fhAnglePairOpeningAngleCut   ; 
  delete fhAnglePairAllCut   ;  
  delete fhInvMassPairNoCut    ; 
  delete fhInvMassPairCorrelationCut    ; 
  delete fhInvMassPairOpeningAngleCut  ; 
  delete fhInvMassPairAllCut   ;    
  
}


//________________________________________________________________________
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) {
  //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)  
{ 
  //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$
	20	* Revision 1.1.2.1 2007/07/26 10:32:09 schutz
	21	* new analysis classes in the the new analysis framework
	22	*
	23	*
	24	*/
	25
	26	//_________________________________________________________________________
	27	// Class that contains methods to select candidate pairs to neutral meson
	28	//-- Author: Gustavo Conesa (INFN-LNF)
	29
	30	// --- ROOT system ---
	31	#include <TParticle.h>
	32	#include <TLorentzVector.h>
	33	#include <TH2.h>
	34	#include <TList.h>
	35
	36	//---- AliRoot system ----
	37	#include "AliNeutralMesonSelection.h"
	38	#include "Riostream.h"
	39	#include "AliLog.h"
	40
	41	ClassImp(AliNeutralMesonSelection)
	42
	43
	44	//____________________________________________________________________________
	45	AliNeutralMesonSelection::AliNeutralMesonSelection() :
	46	TObject(), fSelect(0), fM(0),
	47	fInvMassMaxCut(0.), fInvMassMinCut(0.),
	48	fAngleMaxParam(), fMinPt(0),
	49	fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
	50	fRatioMaxCut(0), fRatioMinCut(0), fKeepNeutralMesonHistos(0),
	51	fhAnglePairNoCut(0), fhAnglePairCorrelationCut(0),
	52	fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0),
	53	fhInvMassPairNoCut(0), fhInvMassPairCorrelationCut(0),
	54	fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0)
	55	{
	56	//Default Ctor
	57
	58	//Initialize parameters
	59
	60	// kGammaHadron and kGammaJet
	61	fAngleMaxParam.Set(4) ;
	62	fAngleMaxParam.Reset(0.);
	63
	64	//Initialize parameters
65	InitParameters();
66	}
67
68	//____________________________________________________________________________
69	AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
70	TObject(),
71	fSelect(g.fSelect), fM(g.fM),
72	fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
73	fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
74	fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
75	fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut),
76	fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
77	fhAnglePairNoCut(g. fhAnglePairNoCut),
78	fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut),
79	fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
80	fhAnglePairAllCut(g. fhAnglePairAllCut),
81	fhInvMassPairNoCut(g.fhInvMassPairNoCut),
82	fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut),
83	fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
84	fhInvMassPairAllCut(g.fhInvMassPairAllCut)
85	{
86	// cpy ctor
87	}
88
89	//_________________________________________________________________________
90	AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
91	{
92	// assignment operator
93
94	if(this == &source)return *this;
95	((TObject *)this)->operator=(source);
96
97	fSelect = source.fSelect ;
98	fM = source.fM ;
99	fInvMassMaxCut = source.fInvMassMaxCut ;
100	fInvMassMinCut = source.fInvMassMinCut ;
101	fAngleMaxParam = source.fAngleMaxParam ;
102	fMinPt = source.fMinPt ;
103	fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
104	fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
105	fRatioMaxCut = source.fRatioMaxCut ;
106	fRatioMinCut = source.fRatioMinCut ;
107	fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
108
109	fhAnglePairNoCut = source. fhAnglePairNoCut ;
110	fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ;
111	fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
112	fhAnglePairAllCut = source. fhAnglePairAllCut ;
113	fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
114	fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ;
115	fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
116	fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
117
118	return *this;
119
120	}
121
122	//____________________________________________________________________________
123	AliNeutralMesonSelection::~AliNeutralMesonSelection()
124	{
125	// Remove all pointers
126
127	delete fhAnglePairNoCut ;
128	delete fhAnglePairCorrelationCut ;
129	delete fhAnglePairOpeningAngleCut ;
130	delete fhAnglePairAllCut ;
131	delete fhInvMassPairNoCut ;
132	delete fhInvMassPairCorrelationCut ;
133	delete fhInvMassPairOpeningAngleCut ;
134	delete fhInvMassPairAllCut ;
135
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	//__________________________________________________________________________-
243	Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) {
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	//____________________________________________________________________________
262	Bool_t AliNeutralMesonSelection::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi)
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	}