PWG/CaloTrackCorrBase/AliNeutralMesonSelection.cxx

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