[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.                  *
 **************************************************************************/
/* $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(),             fAsymmetryCut(1),                fUseAsymmetryCut(0),
    fM(0),                 fInvMassMaxCut(0.),              fInvMassMinCut(0.),           fInvMassMaxCutParam(),
    fAngleMaxParam(),      fUseAngleCut(0),                 
    fShiftMinAngle(),      fKeepNeutralMesonHistos(0), 
    fhAnglePairNoCut(0),   fhAnglePairOpeningAngleCut(0),   fhAnglePairAsymmetryCut(0),   fhAnglePairAllCut(0), 
    fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0), fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
    fhAsymmetryNoCut(0),   fhAsymmetryOpeningAngleCut(0),   fhAsymmetryAllCut(0),
    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(const Float_t angle,const 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, TString 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;
  if(calo=="EMCAL" && 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;
    //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
    //       e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
  }
  
  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;
  
}

//____________________________________________________________________________
void  AliNeutralMesonSelection::SetParticle(TString particleName){
  // Set some default parameters for selection of pi0 or eta
  
  if(particleName=="Pi0"){
    
    fM               = 0.135 ; // GeV
    fInvMassMaxCut   = 0.16  ; // GeV
    fInvMassMinCut   = 0.11  ; // 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) ;
    
    fHistoNIMBins    = 150 ;
    fHistoIMMax      = 0.3 ;
    fHistoIMMin      = 0.  ;  
    
  }else if(particleName=="Eta"){
  
    fM               = 0.547 ; // GeV
    fInvMassMaxCut   = 0.590 ; // GeV
    fInvMassMinCut   = 0.510 ; // GeV
    
    fInvMassMaxCutParam[0] = 0.0 ;
    fInvMassMaxCutParam[1] = 0.0 ;
    fInvMassMaxCutParam[2] = 0.0 ;
    
    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
  
    fHistoNIMBins    = 200  ; // GeV
    fHistoIMMax      = 0.75 ; // GeV
    fHistoIMMin      = 0.35 ; // GeV 
    
  }
  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("mass : %f  \n", 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	**************************************************************************/
	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() :
9542ac75	36	TObject(), fAsymmetryCut(1), fUseAsymmetryCut(0),
3bfcb597	37	fM(0), fInvMassMaxCut(0.), fInvMassMinCut(0.), fInvMassMaxCutParam(),
9542ac75	38	fAngleMaxParam(), fUseAngleCut(0),
7287751b	39	fShiftMinAngle(), fKeepNeutralMesonHistos(0),
9542ac75	40	fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0), fhAnglePairAsymmetryCut(0), fhAnglePairAllCut(0),
	41	fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0), fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
	42	fhAsymmetryNoCut(0), fhAsymmetryOpeningAngleCut(0), fhAsymmetryAllCut(0),
	43	fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
9542ac75	44	fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
9542ac75	45	fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
477d6cee	46	{
1c5acb87	47	//Default Ctor
1c5acb87	48
1c5acb87	49	//Initialize parameters
	50	InitParameters();
	51	}
1c5acb87	52
1c5acb87	53	//________________________________________________________________________
	54	TList * AliNeutralMesonSelection::GetCreateOutputObjects()
	55	{
477d6cee	56	// Create histograms to be saved in output file and
	57	// store them in outputContainer of the analysis class that calls this class.
	58
	59	TList * outputContainer = new TList() ;
	60	outputContainer->SetName("MesonDecayHistos") ;
9542ac75	61
fad7739b	62	if(fKeepNeutralMesonHistos){
	63
	64	outputContainer->SetOwner(kFALSE);
	65
	66	fhAnglePairNoCut = new TH2F
	67	("AnglePairNoCut",
	68	"Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	69	fhAnglePairNoCut->SetYTitle("Angle (rad)");
	70	fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75	71
	72	fhAsymmetryNoCut = new TH2F
	73	("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
3bfcb597	74	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	75	fhAsymmetryNoCut->SetYTitle("Asymmetry");
	76	fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	77
	78	fhInvMassPairNoCut = new TH2F
	79	("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
3bfcb597	80	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	81	fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	82	fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	83
	84	outputContainer->Add(fhAnglePairNoCut) ;
	85	outputContainer->Add(fhAsymmetryNoCut) ;
	86	outputContainer->Add(fhInvMassPairNoCut) ;
	87
	88	if(fUseAngleCut) {
	89	fhAnglePairOpeningAngleCut = new TH2F
	90	("AnglePairOpeningAngleCut",
	91	"Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
	92	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	93	fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
	94	fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	95
	96	fhAsymmetryOpeningAngleCut = new TH2F
	97	("AsymmetryOpeningAngleCut",
	98	"Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597	99	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	100	fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
	101	fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	102
	103	fhInvMassPairOpeningAngleCut = new TH2F
	104	("InvMassPairOpeningAngleCut",
	105	"Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597	106	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	107	fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	108	fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
	109
	110	outputContainer->Add(fhAnglePairOpeningAngleCut) ;
	111	outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
	112	outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
	113	}
	114
	115	if(fUseAsymmetryCut) {
	116	fhAnglePairAsymmetryCut = new TH2F
	117	("AnglePairAsymmetryCut",
	118	"Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
	119	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	120	fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
	121	fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
	122
	123	fhInvMassPairAsymmetryCut = new TH2F
	124	("InvMassPairAsymmetryCut",
	125	"Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
3bfcb597	126	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75	127	fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	128	fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");
	129
	130	outputContainer->Add(fhAnglePairAsymmetryCut) ;
	131	outputContainer->Add(fhInvMassPairAsymmetryCut) ;
	132
	133	}
	134
fad7739b	135	fhAnglePairAllCut = new TH2F
fad7739b	136	("AnglePairAllCut",
9542ac75	137	"Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
fad7739b	138	,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
	139	fhAnglePairAllCut->SetYTitle("Angle (rad)");
	140	fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75	141
fad7739b	142	fhInvMassPairAllCut = new TH2F
fad7739b	143	("InvMassPairAllCut",
9542ac75	144	"Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
3bfcb597	145	fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
fad7739b	146	fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
	147	fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	148
9542ac75	149	fhAsymmetryAllCut = new TH2F
	150	("AsymmetryAllCut",
	151	"Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
3bfcb597	152	fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75	153	fhAsymmetryAllCut->SetYTitle("Asymmetry");
	154	fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
	155
	156	outputContainer->Add(fhAnglePairAllCut) ;
	157	outputContainer->Add(fhAsymmetryAllCut) ;
	158	outputContainer->Add(fhInvMassPairAllCut) ;
	159
fad7739b	160	}
9542ac75	161
477d6cee	162	return outputContainer;
9542ac75	163
1c5acb87	164	}
1c5acb87	165
477d6cee	166	//____________________________________________________________________________
1c5acb87	167	void AliNeutralMesonSelection::InitParameters()
1c5acb87	168	{
477d6cee	169
1c5acb87	170	//Initialize the parameters of the analysis.
1c5acb87	171	fAngleMaxParam.Set(4) ;
9542ac75	172	fAngleMaxParam.Reset(0.);
1c5acb87	173
3bfcb597	174	SetParticle("Pi0");
1c5acb87	175
3bfcb597	176	//Histogrammes settings
	177	fHistoNEBins = 200 ;
	178	fHistoEMax = 50 ;
a5fb4114	179	fHistoEMin = 0. ;
1c5acb87	180
3bfcb597	181	fHistoNAngleBins = 200 ;
3bfcb597	182	fHistoAngleMax = 0.5 ;
1c5acb87	183	fHistoAngleMin = 0. ;
1c5acb87	184
1c5acb87	185	}
	186
	187	//__________________________________________________________________________-
a5fb4114	188	Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const
	189	{
	190
	191	// Check if the opening angle of the candidate pairs is inside
	192	// our selection window
3bfcb597	193	// Attention, only valid for Pi0, if needed for Eta need to revise max angle function or change parameters
9542ac75	194
1c5acb87	195	Double_t max = fAngleMaxParam.At(0)TMath::Exp(fAngleMaxParam.At(1)e)
	196	+fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
	197	Double_t arg = (ee-2fMfM)/(ee);
	198	Double_t min = 100. ;
	199	if(arg>0.)
7287751b	200	min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
3bfcb597	201
9542ac75	202	if((angle<max)&&(angle>=min)) return kTRUE ;
	203	else return kFALSE ;
	204
1c5acb87	205	}
	206
	207	//____________________________________________________________________________
3bfcb597	208	Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj, TString calo)
1c5acb87	209	{
	210
	211	//Search for the neutral pion within selection cuts
1c5acb87	212
9542ac75	213	// Double_t pt = (gammai+gammaj).Pt();
1c5acb87	214	Double_t phi = (gammai+gammaj).Phi();
	215	if(phi < 0)
	216	phi+=TMath::TwoPi();
	217	Double_t invmass = (gammai+gammaj).M();
	218	Double_t angle = gammaj.Angle(gammai.Vect());
	219	Double_t e = (gammai+gammaj).E();
9542ac75	220	Double_t asy = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
9542ac75	221
1c5acb87	222	//Fill histograms with no cuts applied.
fad7739b	223	if(fKeepNeutralMesonHistos){
	224	fhAnglePairNoCut ->Fill(e,angle);
	225	fhInvMassPairNoCut->Fill(e,invmass);
9542ac75	226	fhAsymmetryNoCut ->Fill(e,asy);
fad7739b	227	}
9542ac75	228
1c5acb87	229	//Cut on the aperture of the pair
9542ac75	230	if(fUseAngleCut){
	231	if(IsAngleInWindow(angle,e)){
	232	if(fKeepNeutralMesonHistos ){
	233	fhAnglePairOpeningAngleCut ->Fill(e,angle);
	234	fhInvMassPairOpeningAngleCut->Fill(e,invmass);
	235	fhAsymmetryOpeningAngleCut ->Fill(e,asy);
	236	}
	237	//AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
	238	} else return kFALSE;
	239	}
	240
	241	// Asymmetry cut
	242	if(fUseAsymmetryCut){
	243	if(fAsymmetryCut > asy){
	244	if(fKeepNeutralMesonHistos){
	245	fhInvMassPairAsymmetryCut->Fill(e,invmass);
	246	fhAnglePairAsymmetryCut ->Fill(e,angle);
	247	}
	248	} else return kFALSE;
	249	}
1c5acb87	250
9542ac75	251
9542ac75	252	//Cut on the invariant mass of the pair
3bfcb597	253
	254	Float_t invmassmaxcut = fInvMassMaxCut;
	255	if(calo=="EMCAL" && e > 6.){ // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
	256
	257	invmassmaxcut = (fInvMassMaxCutParam[0]+fInvMassMaxCut)+fInvMassMaxCutParam[1]e+fInvMassMaxCutParam[2]e*e;
	258	//printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
	259	// e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
	260	}
	261
	262	if((invmass > fInvMassMinCut) && (invmass < invmassmaxcut)){
9542ac75	263	if(fKeepNeutralMesonHistos){
	264	fhInvMassPairAllCut->Fill(e,invmass);
	265	fhAnglePairAllCut ->Fill(e,angle);
	266	fhAsymmetryAllCut ->Fill(e,asy);
	267	}
	268
	269	//AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
	270	return kTRUE;
	271
	272	}//(invmass>0.125) && (invmass<0.145)
	273	else return kFALSE;
1c5acb87	274
	275	}
	276
3bfcb597	277	//____________________________________________________________________________
	278	void AliNeutralMesonSelection::SetParticle(TString particleName){
	279	// Set some default parameters for selection of pi0 or eta
	280
	281	if(particleName=="Pi0"){
	282
	283	fM = 0.135 ; // GeV
	284	fInvMassMaxCut = 0.16 ; // GeV
	285	fInvMassMinCut = 0.11 ; // GeV
	286
	287	fInvMassMaxCutParam[0] = 0.0 ;
	288	fInvMassMaxCutParam[1] =-7.e-5 ;
	289	fInvMassMaxCutParam[2] = 8.e-5 ;
	290
7287751b	291	fShiftMinAngle[0] =-0.03 ;
	292	fShiftMinAngle[1] = 0.0025;
	293
	294	fAngleMaxParam.AddAt( 0.8, 0) ;
	295	fAngleMaxParam.AddAt(-1, 1) ;
	296	fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09
	297	fAngleMaxParam.AddAt(-2.e-3,3) ;
3bfcb597	298
	299	fHistoNIMBins = 150 ;
	300	fHistoIMMax = 0.3 ;
	301	fHistoIMMin = 0. ;
	302
	303	}else if(particleName=="Eta"){
	304
	305	fM = 0.547 ; // GeV
	306	fInvMassMaxCut = 0.590 ; // GeV
	307	fInvMassMinCut = 0.510 ; // GeV
	308
	309	fInvMassMaxCutParam[0] = 0.0 ;
	310	fInvMassMaxCutParam[1] = 0.0 ;
7287751b	311	fInvMassMaxCutParam[2] = 0.0 ;
	312
	313	fShiftMinAngle[0] =-0.03 ;
	314	fShiftMinAngle[0] = 0.00 ;
3bfcb597	315
7287751b	316	fAngleMaxParam.AddAt( 0.80, 0) ; // Same as pi0
	317	fAngleMaxParam.AddAt(-0.25, 1) ; // Same as pi0
	318	fAngleMaxParam.AddAt( 0.12, 2) ; // Shifted with respect to pi0
	319	fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
3bfcb597	320
	321	fHistoNIMBins = 200 ; // GeV
	322	fHistoIMMax = 0.75 ; // GeV
	323	fHistoIMMin = 0.35 ; // GeV
	324
	325	}
	326	else
	327	printf("AliAnaNeutralMesonSelection::SetParticle(%s) * Particle NOT defined (Pi0 or Eta), Pi0 settings by default * \n",particleName.Data());
	328
	329
	330	}
	331
1c5acb87	332	//__________________________________________________________________
	333	void AliNeutralMesonSelection::Print(const Option_t * opt) const
	334	{
	335
	336	//Print some relevant parameters set for the analysis
	337	if(! opt)
	338	return;
	339
a3aebfff	340	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
1c5acb87	341
	342	printf("mass : %f \n", fM );
	343	printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
9542ac75	344
	345	printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
	346
	347	printf("Use angle cut? : %d \n", fUseAngleCut );
	348	if(fUseAngleCut){
	349	printf("Angle selection param: \n");
	350	printf("p0 : %f\n", fAngleMaxParam.At(0));
	351	printf("p1 : %f\n", fAngleMaxParam.At(1));
	352	printf("p2 : %f\n", fAngleMaxParam.At(2));
	353	printf("p3 : %f\n", fAngleMaxParam.At(3));
7287751b	354	printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
9542ac75	355	}
9542ac75	356
1c5acb87	357	printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
	358
	359	if(fKeepNeutralMesonHistos){
a5fb4114	360	printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
477d6cee	361	printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
3bfcb597	362	printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
1c5acb87	363	}
477d6cee	364
1c5acb87	365	}