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

//_________________________________________________________________________
// 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(),             fAsymmetryCut(1),                
    fUseAsymmetryCut(0),   fM(0),                 
    fInvMassMaxCut(0.),    fInvMassMinCut(0.),   
    fLeftBandMinCut(0.),   fLeftBandMaxCut(0.),            
    fRightBandMinCut(0.),  fRightBandMaxCut(0.),              
    fAngleMaxParam(),      fUseAngleCut(0),       
    fKeepNeutralMesonHistos(0),
    fParticle(""),         fDecayBit(0),
    // histograms
    fhAnglePairNoCut(0),          fhAnglePairOpeningAngleCut(0),   
    fhAnglePairAsymmetryCut(0),   fhAnglePairAllCut(0), 
    fhInvMassPairNoCut(0),        fhInvMassPairOpeningAngleCut(0), 
    fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
    fhAsymmetryNoCut(0),          fhAsymmetryOpeningAngleCut(0),   
    fhAsymmetryAllCut(0),
    // histogram ranges and bins
    fHistoNEBins(0),       fHistoEMax(0.),                  fHistoEMin(0.),
    fHistoNAngleBins(0),   fHistoAngleMax(0.),              fHistoAngleMin(0.),
    fHistoNIMBins(0),      fHistoIMMax(0.),                 fHistoIMMin(0.)
{
  //Default Ctor
  
  //Initialize parameters
  InitParameters();
}

//_________________________________________________________
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)");
	      
    fhAsymmetryNoCut  = new TH2F
	  ("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
	   fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
	  fhAsymmetryNoCut->SetYTitle("Asymmetry");
	  fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
    
    fhInvMassPairNoCut  = new TH2F
	  ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
	   fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
	  fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	  fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
    
    outputContainer->Add(fhAnglePairNoCut) ; 
	  outputContainer->Add(fhAsymmetryNoCut) ; 
    outputContainer->Add(fhInvMassPairNoCut) ; 

    if(fUseAngleCut) {
      fhAnglePairOpeningAngleCut  = new TH2F
      ("AnglePairOpeningAngleCut",
       "Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
       ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
      fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
      fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
      
      fhAsymmetryOpeningAngleCut  = new TH2F
      ("AsymmetryOpeningAngleCut",
       "Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
       fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
      fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
      fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");   
      
      fhInvMassPairOpeningAngleCut  = new TH2F
      ("InvMassPairOpeningAngleCut",
       "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
       fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
      fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
      fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
      
      outputContainer->Add(fhAnglePairOpeningAngleCut) ;
      outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
      outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
    }
    
	  if(fUseAsymmetryCut) {
      fhAnglePairAsymmetryCut  = new TH2F
      ("AnglePairAsymmetryCut",
       "Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
       ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
      fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
      fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
      
      fhInvMassPairAsymmetryCut  = new TH2F
      ("InvMassPairAsymmetryCut",
       "Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
       fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
      fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
      fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");      
      
      outputContainer->Add(fhAnglePairAsymmetryCut) ;
      outputContainer->Add(fhInvMassPairAsymmetryCut) ;
      
    }
    
	  fhAnglePairAllCut  = new TH2F
	  ("AnglePairAllCut",
	   "Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
	   ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
	  fhAnglePairAllCut->SetYTitle("Angle (rad)");
	  fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
        
	  fhInvMassPairAllCut  = new TH2F
	  ("InvMassPairAllCut",
	   "Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
	   fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
	  fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	  fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	  
	  fhAsymmetryAllCut  = new TH2F
	  ("AsymmetryAllCut",
	   "Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
	   fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
	  fhAsymmetryAllCut->SetYTitle("Asymmetry");
	  fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
    
    outputContainer->Add(fhAnglePairAllCut) ; 
	  outputContainer->Add(fhAsymmetryAllCut) ; 
    outputContainer->Add(fhInvMassPairAllCut) ;     

  }
  
  return outputContainer;

}

//_____________________________________________
void AliNeutralMesonSelection::InitParameters()
{
  
  //Initialize the parameters of the analysis.
  fAngleMaxParam.Set(4) ;
  fAngleMaxParam.Reset(0.);

  SetParticle("Pi0");
  
  //Histogrammes settings
  fHistoNEBins     = 200 ;
  fHistoEMax       = 50  ;
  fHistoEMin       = 0.  ;  

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

}

//______________________________________________________________________________
Bool_t AliNeutralMesonSelection::IsAngleInWindow(Float_t angle, Float_t e) const
{
 
  // Check if the opening angle of the candidate pairs is inside 
  // our selection window
  // Attention, only valid for Pi0, if needed for Eta need to revise max angle function or change parameters
  	
  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)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
  
  if((angle<max)&&(angle>=min)) return kTRUE  ;
  else                          return kFALSE ;

}

//_________________________________________________________________
Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, 
                                             TLorentzVector gammaj, 
                                             Int_t calo)
{  
  
  //Search for the neutral pion within selection cuts
  
  //  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();
  Double_t asy     = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();

  //Fill histograms with no cuts applied.
  if(fKeepNeutralMesonHistos)
  {
	  fhAnglePairNoCut  ->Fill(e,angle);
	  fhInvMassPairNoCut->Fill(e,invmass);
	  fhAsymmetryNoCut  ->Fill(e,asy);
  }
  
  //Cut on the aperture of the pair
  if(fUseAngleCut)
  {
    if(IsAngleInWindow(angle,e))
    {
      if(fKeepNeutralMesonHistos )
      {
        fhAnglePairOpeningAngleCut  ->Fill(e,angle);
        fhInvMassPairOpeningAngleCut->Fill(e,invmass);
        fhAsymmetryOpeningAngleCut  ->Fill(e,asy);
      }
      //AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
    } else return kFALSE;
  }
  
  // Asymmetry cut
  if(fUseAsymmetryCut)
  {
    if(fAsymmetryCut > asy)
    {
      if(fKeepNeutralMesonHistos)
      {
        fhInvMassPairAsymmetryCut->Fill(e,invmass);
        fhAnglePairAsymmetryCut  ->Fill(e,angle);
      }
    } else return kFALSE;
  }
  
  
  //Cut on the invariant mass of the pair
  
  Float_t invmassmaxcut = fInvMassMaxCut;
  Float_t invmassRightBandMinCut = fRightBandMinCut;
  Float_t invmassRightBandMixCut = fRightBandMaxCut;

  // kEMCAL=0, kPHOS=1
  if(calo==0 && e > 6.)
  { // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
  
    invmassmaxcut = (fInvMassMaxCutParam[0]+fInvMassMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
    invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
    invmassRightBandMixCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;

    //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
    //       e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
  }
  
  // normal case, invariant mass selection around pi0/eta peak
  if( !fParticle.Contains("SideBand") )
  {
    if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
    { 
      if(fKeepNeutralMesonHistos)
      {
        fhInvMassPairAllCut->Fill(e,invmass);
        fhAnglePairAllCut  ->Fill(e,angle);
        fhAsymmetryAllCut  ->Fill(e,asy);
      }      
      
      //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
      return kTRUE;
      
    }//(invmass>0.125) && (invmass<0.145)
    else 
    {
      return kFALSE;
    }
  }// normal selection
 
  else // select a band around pi0/eta
  {
    if((invmass > fLeftBandMinCut  && invmass < fLeftBandMaxCut ) ||  
       (invmass > invmassRightBandMinCut && invmass < invmassRightBandMixCut))
    { 
      if(fKeepNeutralMesonHistos)
      {
        fhInvMassPairAllCut->Fill(e,invmass);
        fhAnglePairAllCut  ->Fill(e,angle);
        fhAsymmetryAllCut  ->Fill(e,asy);
      }      
      
      //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
      return kTRUE;
      
    }//(invmass>0.125) && (invmass<0.145)
    else 
    {
      return kFALSE;
    }
  }
  
}

//_______________________________________________________________
void  AliNeutralMesonSelection::SetParticle(TString particleName)
{
  // Set some default parameters for selection of pi0 or eta

  fParticle = particleName ;
  fDecayBit = kPi0;

  if(particleName.Contains("Pi0"))
  {
    fHistoNIMBins          = 150 ;
    fHistoIMMax            = 0.3 ;
    fHistoIMMin            = 0.  ;  
    
    fM                     = 0.135 ; // GeV
    fInvMassMaxCut         = 0.16  ; // GeV
    fInvMassMinCut         = 0.11  ; // GeV

    fLeftBandMinCut        = 0.05  ; // GeV
    fLeftBandMaxCut        = 0.09  ; // GeV
    fRightBandMinCut       = 0.160 ; // GeV
    fRightBandMaxCut       = 0.165 ; // GeV

    fInvMassMaxCutParam[0] = 0.0   ;
    fInvMassMaxCutParam[1] =-7.e-5 ;
    fInvMassMaxCutParam[2] = 8.e-5 ;    
         
    fShiftMinAngle[0]      =-0.03  ;
    fShiftMinAngle[1]      = 0.0025;
    
    fAngleMaxParam.AddAt( 0.8,  0) ;
    fAngleMaxParam.AddAt(-1,    1) ;
    fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09 
    fAngleMaxParam.AddAt(-2.e-3,3) ;
    
    fDecayBit = kPi0;
    if(particleName.Contains("Side")) fDecayBit = kPi0Side;
  
  }  
  else if(particleName.Contains("Eta"))
  {
    fHistoNIMBins          = 200  ; // GeV
    fHistoIMMax            = 0.75 ; // GeV
    fHistoIMMin            = 0.35 ; // GeV 
  
    fM                     = 0.547 ; // GeV
    fInvMassMaxCut         = 0.590 ; // GeV
    fInvMassMinCut         = 0.510 ; // GeV
    
    fLeftBandMinCut        = 0.450 ; // GeV
    fLeftBandMaxCut        = 0.500 ; // GeV
    fRightBandMinCut       = 0.600 ; // GeV
    fRightBandMaxCut       = 0.650 ; // GeV
    
    fInvMassMaxCutParam[0] = 0.00 ;
    fInvMassMaxCutParam[1] = 0.00 ;
    fInvMassMaxCutParam[2] = 0.00 ;
    
    fShiftMinAngle[0]      =-0.03 ;
    fShiftMinAngle[0]      = 0.00 ;
    
    fAngleMaxParam.AddAt( 0.80,  0) ; // Same as pi0
    fAngleMaxParam.AddAt(-0.25,  1) ; // Same as pi0
    fAngleMaxParam.AddAt( 0.12,  2) ; // Shifted with respect to pi0
    fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
    
    fDecayBit = kEta;
    if(particleName.Contains("Side")) fDecayBit = kEtaSide;
  }
  else 
    printf("AliAnaNeutralMesonSelection::SetParticle(%s) *** Particle NOT defined (Pi0 or Eta), Pi0 settings by default *** \n",particleName.Data());
  
}

//______________________________________________________________
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("Particle %s, bit %d, mass : %f  \n", fParticle.Data(), fDecayBit, fM );
  printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
  
  printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
  
  printf("Use angle cut? : %d  \n", fUseAngleCut );
  if(fUseAngleCut)
  {
    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("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
  }
  
  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 < 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	**************************************************************************/
1c5acb87	15
	16	//_________________________________________________________________________
	17	// Class that contains methods to select candidate pairs to neutral meson
	18	// 2 main selections, invariant mass around pi0 (also any other mass),
	19	// apperture angle to distinguish from combinatorial.
	20	//-- Author: Gustavo Conesa (INFN-LNF)
	21
	22	// --- ROOT system ---
	23	#include <TLorentzVector.h>
	24	#include <TH2.h>
	25	#include <TList.h>
3f150b4b	26
1c5acb87	27	//---- AliRoot system ----
1c5acb87	28	#include "AliNeutralMesonSelection.h"
1c5acb87	29
	30	ClassImp(AliNeutralMesonSelection)
	31
	32
b960c7eb	33	//______________________________________________________
1c5acb87	34	AliNeutralMesonSelection::AliNeutralMesonSelection() :
3f150b4b	35	TObject(), fAsymmetryCut(1),
	36	fUseAsymmetryCut(0), fM(0),
	37	fInvMassMaxCut(0.), fInvMassMinCut(0.),
	38	fLeftBandMinCut(0.), fLeftBandMaxCut(0.),
	39	fRightBandMinCut(0.), fRightBandMaxCut(0.),
	40	fAngleMaxParam(), fUseAngleCut(0),
a39c9b1b	41	fKeepNeutralMesonHistos(0),
	42	fParticle(""), fDecayBit(0),
	43	// histograms
3f150b4b	44	fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
	45	fhAnglePairAsymmetryCut(0), fhAnglePairAllCut(0),
	46	fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
	47	fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
	48	fhAsymmetryNoCut(0), fhAsymmetryOpeningAngleCut(0),
	49	fhAsymmetryAllCut(0),
	50	// histogram ranges and bins
9542ac75	51	fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
9542ac75	52	fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
9542ac75	53	fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
477d6cee	54	{
1c5acb87	55	//Default Ctor
1c5acb87	56
1c5acb87	57	//Initialize parameters
	58	InitParameters();
	59	}
1c5acb87	60
b960c7eb	61	//_________________________________________________________
1c5acb87	62	TList * AliNeutralMesonSelection::GetCreateOutputObjects()
1c5acb87	63	{
477d6cee	64	// Create histograms to be saved in output file and
	65	// store them in outputContainer of the analysis class that calls this class.
	66
	67	TList * outputContainer = new TList() ;
	68	outputContainer->SetName("MesonDecayHistos") ;
9542ac75	69
fad7739b	70	if(fKeepNeutralMesonHistos){
	71
	72	outputContainer->SetOwner(kFALSE);
	73
	74	fhAnglePairNoCut = new TH2F
	75	("AnglePairNoCut",
	76	"Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	77	fhAnglePairNoCut->SetYTitle("Angle (rad)");
	78	fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75	79
	80	fhAsymmetryNoCut = new TH2F
	81	("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
3bfcb597	82	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	83	fhAsymmetryNoCut->SetYTitle("Asymmetry");
	84	fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	85
	86	fhInvMassPairNoCut = new TH2F
	87	("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
3bfcb597	88	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	89	fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	90	fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	91
	92	outputContainer->Add(fhAnglePairNoCut) ;
	93	outputContainer->Add(fhAsymmetryNoCut) ;
	94	outputContainer->Add(fhInvMassPairNoCut) ;
	95
	96	if(fUseAngleCut) {
	97	fhAnglePairOpeningAngleCut = new TH2F
	98	("AnglePairOpeningAngleCut",
	99	"Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
	100	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	101	fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
	102	fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	103
	104	fhAsymmetryOpeningAngleCut = new TH2F
	105	("AsymmetryOpeningAngleCut",
	106	"Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597	107	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	108	fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
	109	fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	110
	111	fhInvMassPairOpeningAngleCut = new TH2F
	112	("InvMassPairOpeningAngleCut",
	113	"Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597	114	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	115	fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	116	fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	117
	118	outputContainer->Add(fhAnglePairOpeningAngleCut) ;
	119	outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
	120	outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
	121	}
	122
	123	if(fUseAsymmetryCut) {
	124	fhAnglePairAsymmetryCut = new TH2F
	125	("AnglePairAsymmetryCut",
	126	"Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
	127	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	128	fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
	129	fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	130
	131	fhInvMassPairAsymmetryCut = new TH2F
	132	("InvMassPairAsymmetryCut",
	133	"Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
3bfcb597	134	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	135	fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	136	fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");
	137
	138	outputContainer->Add(fhAnglePairAsymmetryCut) ;
	139	outputContainer->Add(fhInvMassPairAsymmetryCut) ;
	140
	141	}
	142
fad7739b	143	fhAnglePairAllCut = new TH2F
fad7739b	144	("AnglePairAllCut",
9542ac75	145	"Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
fad7739b	146	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	147	fhAnglePairAllCut->SetYTitle("Angle (rad)");
	148	fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75	149
fad7739b	150	fhInvMassPairAllCut = new TH2F
fad7739b	151	("InvMassPairAllCut",
9542ac75	152	"Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
3bfcb597	153	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
fad7739b	154	fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	155	fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	156
9542ac75	157	fhAsymmetryAllCut = new TH2F
	158	("AsymmetryAllCut",
	159	"Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
3bfcb597	160	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	161	fhAsymmetryAllCut->SetYTitle("Asymmetry");
	162	fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	163
	164	outputContainer->Add(fhAnglePairAllCut) ;
	165	outputContainer->Add(fhAsymmetryAllCut) ;
	166	outputContainer->Add(fhInvMassPairAllCut) ;
	167
fad7739b	168	}
9542ac75	169
477d6cee	170	return outputContainer;
9542ac75	171
1c5acb87	172	}
1c5acb87	173
b960c7eb	174	//_____________________________________________
1c5acb87	175	void AliNeutralMesonSelection::InitParameters()
1c5acb87	176	{
477d6cee	177
1c5acb87	178	//Initialize the parameters of the analysis.
1c5acb87	179	fAngleMaxParam.Set(4) ;
9542ac75	180	fAngleMaxParam.Reset(0.);
1c5acb87	181
3bfcb597	182	SetParticle("Pi0");
1c5acb87	183
3bfcb597	184	//Histogrammes settings
	185	fHistoNEBins = 200 ;
	186	fHistoEMax = 50 ;
a5fb4114	187	fHistoEMin = 0. ;
1c5acb87	188
3bfcb597	189	fHistoNAngleBins = 200 ;
3bfcb597	190	fHistoAngleMax = 0.5 ;
1c5acb87	191	fHistoAngleMin = 0. ;
1c5acb87	192
1c5acb87	193	}
1c5acb87	194
8a2dbbff	195	//______________________________________________________________________________
8a2dbbff	196	Bool_t AliNeutralMesonSelection::IsAngleInWindow(Float_t angle, Float_t e) const
a5fb4114	197	{
	198
	199	// Check if the opening angle of the candidate pairs is inside
	200	// our selection window
3bfcb597	201	// Attention, only valid for Pi0, if needed for Eta need to revise max angle function or change parameters
9542ac75	202
1c5acb87	203	Double_t max = fAngleMaxParam.At(0)TMath::Exp(fAngleMaxParam.At(1)e)
	204	+fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
	205	Double_t arg = (ee-2fMfM)/(ee);
	206	Double_t min = 100. ;
	207	if(arg>0.)
7287751b	208	min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
3bfcb597	209
9542ac75	210	if((angle<max)&&(angle>=min)) return kTRUE ;
	211	else return kFALSE ;
	212
1c5acb87	213	}
1c5acb87	214
b960c7eb	215	//_________________________________________________________________
	216	Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai,
	217	TLorentzVector gammaj,
bb787aa8	218	Int_t calo)
1c5acb87	219	{
	220
	221	//Search for the neutral pion within selection cuts
1c5acb87	222
9542ac75	223	// Double_t pt = (gammai+gammaj).Pt();
1c5acb87	224	Double_t phi = (gammai+gammaj).Phi();
	225	if(phi < 0)
	226	phi+=TMath::TwoPi();
b960c7eb	227
1c5acb87	228	Double_t invmass = (gammai+gammaj).M();
	229	Double_t angle = gammaj.Angle(gammai.Vect());
	230	Double_t e = (gammai+gammaj).E();
9542ac75	231	Double_t asy = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
9542ac75	232
1c5acb87	233	//Fill histograms with no cuts applied.
b960c7eb	234	if(fKeepNeutralMesonHistos)
b960c7eb	235	{
fad7739b	236	fhAnglePairNoCut ->Fill(e,angle);
fad7739b	237	fhInvMassPairNoCut->Fill(e,invmass);
9542ac75	238	fhAsymmetryNoCut ->Fill(e,asy);
fad7739b	239	}
9542ac75	240
1c5acb87	241	//Cut on the aperture of the pair
b960c7eb	242	if(fUseAngleCut)
	243	{
	244	if(IsAngleInWindow(angle,e))
	245	{
	246	if(fKeepNeutralMesonHistos )
	247	{
9542ac75	248	fhAnglePairOpeningAngleCut ->Fill(e,angle);
	249	fhInvMassPairOpeningAngleCut->Fill(e,invmass);
	250	fhAsymmetryOpeningAngleCut ->Fill(e,asy);
	251	}
	252	//AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
	253	} else return kFALSE;
	254	}
	255
	256	// Asymmetry cut
b960c7eb	257	if(fUseAsymmetryCut)
	258	{
	259	if(fAsymmetryCut > asy)
	260	{
	261	if(fKeepNeutralMesonHistos)
	262	{
9542ac75	263	fhInvMassPairAsymmetryCut->Fill(e,invmass);
	264	fhAnglePairAsymmetryCut ->Fill(e,angle);
	265	}
	266	} else return kFALSE;
	267	}
1c5acb87	268
9542ac75	269
9542ac75	270	//Cut on the invariant mass of the pair
3bfcb597	271
3bfcb597	272	Float_t invmassmaxcut = fInvMassMaxCut;
80515912	273	Float_t invmassRightBandMinCut = fRightBandMinCut;
	274	Float_t invmassRightBandMixCut = fRightBandMaxCut;
	275
bb787aa8	276	// kEMCAL=0, kPHOS=1
bb787aa8	277	if(calo==0 && e > 6.)
b960c7eb	278	{ // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
3bfcb597	279
3bfcb597	280	invmassmaxcut = (fInvMassMaxCutParam[0]+fInvMassMaxCut)+fInvMassMaxCutParam[1]e+fInvMassMaxCutParam[2]e*e;
80515912	281	invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]e+fInvMassMaxCutParam[2]e*e;
	282	invmassRightBandMixCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]e+fInvMassMaxCutParam[2]e*e;
	283
3bfcb597	284	//printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
	285	// e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
	286	}
	287
3f150b4b	288	// normal case, invariant mass selection around pi0/eta peak
	289	if( !fParticle.Contains("SideBand") )
	290	{
	291	if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
	292	{
	293	if(fKeepNeutralMesonHistos)
	294	{
	295	fhInvMassPairAllCut->Fill(e,invmass);
	296	fhAnglePairAllCut ->Fill(e,angle);
	297	fhAsymmetryAllCut ->Fill(e,asy);
	298	}
	299
	300	//AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
	301	return kTRUE;
	302
	303	}//(invmass>0.125) && (invmass<0.145)
	304	else
b960c7eb	305	{
3f150b4b	306	return kFALSE;
	307	}
	308	}// normal selection
	309
	310	else // select a band around pi0/eta
	311	{
	312	if((invmass > fLeftBandMinCut && invmass < fLeftBandMaxCut ) \|\|
80515912	313	(invmass > invmassRightBandMinCut && invmass < invmassRightBandMixCut))
3f150b4b	314	{
	315	if(fKeepNeutralMesonHistos)
	316	{
	317	fhInvMassPairAllCut->Fill(e,invmass);
	318	fhAnglePairAllCut ->Fill(e,angle);
	319	fhAsymmetryAllCut ->Fill(e,asy);
	320	}
	321
	322	//AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
	323	return kTRUE;
	324
	325	}//(invmass>0.125) && (invmass<0.145)
	326	else
	327	{
	328	return kFALSE;
	329	}
	330	}
1c5acb87	331
	332	}
	333
b960c7eb	334	//_______________________________________________________________
	335	void AliNeutralMesonSelection::SetParticle(TString particleName)
	336	{
3bfcb597	337	// Set some default parameters for selection of pi0 or eta
3f150b4b	338
3f150b4b	339	fParticle = particleName ;
a39c9b1b	340	fDecayBit = kPi0;
a39c9b1b	341
3f150b4b	342	if(particleName.Contains("Pi0"))
b960c7eb	343	{
	344	fHistoNIMBins = 150 ;
	345	fHistoIMMax = 0.3 ;
	346	fHistoIMMin = 0. ;
3bfcb597	347
b960c7eb	348	fM = 0.135 ; // GeV
	349	fInvMassMaxCut = 0.16 ; // GeV
	350	fInvMassMinCut = 0.11 ; // GeV
3bfcb597	351
3f150b4b	352	fLeftBandMinCut = 0.05 ; // GeV
	353	fLeftBandMaxCut = 0.09 ; // GeV
	354	fRightBandMinCut = 0.160 ; // GeV
	355	fRightBandMaxCut = 0.165 ; // GeV
	356
3bfcb597	357	fInvMassMaxCutParam[0] = 0.0 ;
	358	fInvMassMaxCutParam[1] =-7.e-5 ;
	359	fInvMassMaxCutParam[2] = 8.e-5 ;
	360
b960c7eb	361	fShiftMinAngle[0] =-0.03 ;
b960c7eb	362	fShiftMinAngle[1] = 0.0025;
7287751b	363
	364	fAngleMaxParam.AddAt( 0.8, 0) ;
	365	fAngleMaxParam.AddAt(-1, 1) ;
	366	fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09
	367	fAngleMaxParam.AddAt(-2.e-3,3) ;
a39c9b1b	368
	369	fDecayBit = kPi0;
	370	if(particleName.Contains("Side")) fDecayBit = kPi0Side;
3bfcb597	371
3f150b4b	372	}
3f150b4b	373	else if(particleName.Contains("Eta"))
b960c7eb	374	{
	375	fHistoNIMBins = 200 ; // GeV
	376	fHistoIMMax = 0.75 ; // GeV
	377	fHistoIMMin = 0.35 ; // GeV
	378
	379	fM = 0.547 ; // GeV
	380	fInvMassMaxCut = 0.590 ; // GeV
	381	fInvMassMinCut = 0.510 ; // GeV
3bfcb597	382
3f150b4b	383	fLeftBandMinCut = 0.450 ; // GeV
	384	fLeftBandMaxCut = 0.500 ; // GeV
	385	fRightBandMinCut = 0.600 ; // GeV
	386	fRightBandMaxCut = 0.650 ; // GeV
7287751b	387
3f150b4b	388	fInvMassMaxCutParam[0] = 0.00 ;
	389	fInvMassMaxCutParam[1] = 0.00 ;
	390	fInvMassMaxCutParam[2] = 0.00 ;
	391
	392	fShiftMinAngle[0] =-0.03 ;
	393	fShiftMinAngle[0] = 0.00 ;
3bfcb597	394
7287751b	395	fAngleMaxParam.AddAt( 0.80, 0) ; // Same as pi0
7287751b	396	fAngleMaxParam.AddAt(-0.25, 1) ; // Same as pi0
a39c9b1b	397	fAngleMaxParam.AddAt( 0.12, 2) ; // Shifted with respect to pi0
7287751b	398	fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
3bfcb597	399
a39c9b1b	400	fDecayBit = kEta;
a39c9b1b	401	if(particleName.Contains("Side")) fDecayBit = kEtaSide;
3bfcb597	402	}
	403	else
	404	printf("AliAnaNeutralMesonSelection::SetParticle(%s) * Particle NOT defined (Pi0 or Eta), Pi0 settings by default * \n",particleName.Data());
	405
3bfcb597	406	}
3bfcb597	407
b960c7eb	408	//______________________________________________________________
1c5acb87	409	void AliNeutralMesonSelection::Print(const Option_t * opt) const
	410	{
	411
	412	//Print some relevant parameters set for the analysis
	413	if(! opt)
	414	return;
	415
a3aebfff	416	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
1c5acb87	417
a39c9b1b	418	printf("Particle %s, bit %d, mass : %f \n", fParticle.Data(), fDecayBit, fM );
1c5acb87	419	printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
9542ac75	420
	421	printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
	422
	423	printf("Use angle cut? : %d \n", fUseAngleCut );
3f150b4b	424	if(fUseAngleCut)
3f150b4b	425	{
9542ac75	426	printf("Angle selection param: \n");
	427	printf("p0 : %f\n", fAngleMaxParam.At(0));
	428	printf("p1 : %f\n", fAngleMaxParam.At(1));
	429	printf("p2 : %f\n", fAngleMaxParam.At(2));
	430	printf("p3 : %f\n", fAngleMaxParam.At(3));
7287751b	431	printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
9542ac75	432	}
9542ac75	433
1c5acb87	434	printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
1c5acb87	435
3f150b4b	436	if(fKeepNeutralMesonHistos)
3f150b4b	437	{
a5fb4114	438	printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
477d6cee	439	printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
3bfcb597	440	printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
1c5acb87	441	}
477d6cee	442
1c5acb87	443	}