PWG4/PartCorrBase/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 /* $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>
  27 #include <TArrayD.h>
  28
  29 //---- AliRoot system ----
  30 #include "AliNeutralMesonSelection.h"
  31 #include "AliLog.h"
  32
  33 ClassImp(AliNeutralMesonSelection)
  34
  35
  36 //____________________________________________________________________________
  37   AliNeutralMesonSelection::AliNeutralMesonSelection() :
  38     TObject(), fM(0),
  39     fInvMassMaxCut(0.), fInvMassMinCut(0.),
  40     fAngleMaxParam(),  fKeepNeutralMesonHistos(0),
  41     fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
  42     fhAnglePairAllCut(0),
  43     fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
  44     fhInvMassPairAllCut(0),
  45         fHistoNEBins(0),   fHistoEMax(0.),   fHistoEMin(0.),
  46         fHistoNPtBins(0),  fHistoPtMax(0.),  fHistoPtMin(0.),
  47         fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
  48         fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.){
  49   //Default Ctor
  50
  51   //Initialize parameters
  52
  53   // kGammaHadron and kGammaJet
  54   fAngleMaxParam.Set(4) ;
  55   fAngleMaxParam.Reset(0.);
  56
  57   //Initialize parameters
  58   InitParameters();
  59 }
  60
  61 //____________________________________________________________________________
  62 AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
  63   TObject(), fM(g.fM),
  64   fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
  65   fAngleMaxParam(g.fAngleMaxParam),
  66   fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
  67   fhAnglePairNoCut(g. fhAnglePairNoCut),
  68   fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
  69   fhAnglePairAllCut(g. fhAnglePairAllCut),
  70   fhInvMassPairNoCut(g.fhInvMassPairNoCut),
  71   fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
  72   fhInvMassPairAllCut(g.fhInvMassPairAllCut),
  73   fHistoNEBins(g.fHistoNEBins),   fHistoEMax(g.fHistoEMax),   fHistoEMin(g.fHistoEMin),
  74   fHistoNPtBins(g.fHistoNPtBins),   fHistoPtMax(g.fHistoPtMax),   fHistoPtMin(g.fHistoPtMin),
  75   fHistoNAngleBins(g.fHistoNAngleBins), fHistoAngleMax(g.fHistoAngleMax), fHistoAngleMin(g.fHistoAngleMin),
  76   fHistoNIMBins(g.fHistoNIMBins), fHistoIMMax(g.fHistoIMMax), fHistoIMMin(g.fHistoIMMin)
  77 {
  78   // cpy ctor
  79 }
  80
  81 //_________________________________________________________________________
  82 AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
  83 {
  84   // assignment operator
  85
  86   if(this == &source)return *this;
  87   ((TObject *)this)->operator=(source);
  88
  89   fM = source.fM ;
  90   fInvMassMaxCut = source.fInvMassMaxCut ;
  91   fInvMassMinCut = source.fInvMassMinCut ;
  92   fAngleMaxParam = source.fAngleMaxParam ;
  93   fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
  94
  95   fhAnglePairNoCut = source. fhAnglePairNoCut ;
  96   fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
  97   fhAnglePairAllCut = source. fhAnglePairAllCut ;
  98   fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
  99   fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
 100   fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
 101
 102   fHistoNEBins = source.fHistoNEBins;   fHistoEMax = source.fHistoEMax;   fHistoEMin = source.fHistoEMin;
 103   fHistoNPtBins = source.fHistoNPtBins;   fHistoPtMax = source.fHistoPtMax;   fHistoPtMin = source.fHistoPtMin;
 104   fHistoNAngleBins = source.fHistoNAngleBins; fHistoAngleMax = source.fHistoAngleMax; fHistoAngleMin = source.fHistoAngleMin;
 105   fHistoNIMBins = source.fHistoNIMBins; fHistoIMMax = source.fHistoIMMax; fHistoIMMin = source.fHistoIMMin;
 106
 107   return *this;
 108
 109 }
 110
 111 //____________________________________________________________________________
 112 AliNeutralMesonSelection::~AliNeutralMesonSelection()
 113 {
 114   //dtor
 115
 116   if(!fKeepNeutralMesonHistos){
 117     //Histograms initialized and filled but not passed to output container
 118     //delete here, I am not sure this is correct
 119
 120     if(fhAnglePairNoCut) delete fhAnglePairNoCut;
 121     if(fhAnglePairOpeningAngleCut) delete fhAnglePairOpeningAngleCut;
 122     if(fhAnglePairAllCut) delete fhAnglePairAllCut;
 123     if(fhInvMassPairNoCut) delete fhInvMassPairNoCut;
 124     if(fhInvMassPairOpeningAngleCut) delete fhInvMassPairOpeningAngleCut;
 125     if(fhInvMassPairAllCut) delete fhInvMassPairAllCut;
 126
 127   }
 128
 129 }
 130 //________________________________________________________________________
 131 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
 132 {
 133
 134         // Create histograms to be saved in output file and
 135         // store them in outputContainer of the analysis class that calls this class.
 136
 137         TList * outputContainer = new TList() ;
 138         outputContainer->SetName("MesonDecayHistos") ;
 139
 140         fhAnglePairNoCut  = new TH2F
 141     ("AnglePairNoCut",
 142      "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
 143         fhAnglePairNoCut->SetYTitle("Angle (rad)");
 144         fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 145
 146         fhAnglePairOpeningAngleCut  = new TH2F
 147     ("AnglePairOpeningAngleCut",
 148      "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
 149      ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
 150         fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
 151         fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 152
 153         fhAnglePairAllCut  = new TH2F
 154     ("AnglePairAllCut",
 155      "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
 156      ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
 157         fhAnglePairAllCut->SetYTitle("Angle (rad)");
 158         fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 159
 160         //
 161         fhInvMassPairNoCut  = new TH2F
 162     ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
 163      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
 164         fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 165         fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 166
 167         fhInvMassPairOpeningAngleCut  = new TH2F
 168     ("InvMassPairOpeningAngleCut",
 169      "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
 170      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
 171         fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 172         fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
 173
 174         fhInvMassPairAllCut  = new TH2F
 175     ("InvMassPairAllCut",
 176      "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
 177      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
 178         fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 179         fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
 180
 181         outputContainer->Add(fhAnglePairNoCut) ;
 182         outputContainer->Add(fhAnglePairOpeningAngleCut) ;
 183         outputContainer->Add(fhAnglePairAllCut) ;
 184
 185         outputContainer->Add(fhInvMassPairNoCut) ;
 186         outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
 187         outputContainer->Add(fhInvMassPairAllCut) ;
 188
 189         return outputContainer;
 190 }
 191
 192  //____________________________________________________________________________
 193 void AliNeutralMesonSelection::InitParameters()
 194 {
 195
 196   //Initialize the parameters of the analysis.
 197   fKeepNeutralMesonHistos = kFALSE ;
 198
 199   fAngleMaxParam.Set(4) ;
 200   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
 201   fAngleMaxParam.AddAt(-0.25,1) ;
 202   fAngleMaxParam.AddAt(0.025,2) ;
 203   fAngleMaxParam.AddAt(-2e-4,3) ;
 204
 205   fInvMassMaxCut  = 0.16 ;
 206   fInvMassMinCut  = 0.11 ;
 207
 208   fM = 0.1349766;//neutralMeson mass, pi0
 209
 210  //Histogrammes settings
 211   fHistoNEBins = 100 ;
 212   fHistoEMax   = 50 ;
 213   fHistoEMin   = 0.  ;
 214
 215   fHistoNPtBins = 240 ;
 216   fHistoPtMax   = 120 ;
 217   fHistoPtMin   = 0.  ;
 218
 219   fHistoNAngleBins = 200 ;
 220   fHistoAngleMax   = 0.2;
 221   fHistoAngleMin   = 0.  ;
 222
 223   fHistoNIMBins = 300 ;
 224   fHistoIMMax   = 0.5   ;
 225   fHistoIMMin   = 0.  ;
 226 }
 227
 228 //__________________________________________________________________________-
 229 Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
 230   //Check if the opening angle of the candidate pairs is inside
 231   //our selection windowd
 232
 233   Bool_t result = kFALSE;
 234   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
 235     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
 236   Double_t arg = (e*e-2*fM*fM)/(e*e);
 237   Double_t min = 100. ;
 238   if(arg>0.)
 239     min = TMath::ACos(arg);
 240
 241   if((angle<max)&&(angle>=min))
 242     result = kTRUE;
 243
 244   return result;
 245 }
 246
 247 //____________________________________________________________________________
 248 Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)
 249 {
 250
 251   //Search for the neutral pion within selection cuts
 252   Bool_t goodpair = kFALSE ;
 253
 254   Double_t pt  = (gammai+gammaj).Pt();
 255   Double_t phi = (gammai+gammaj).Phi();
 256   if(phi < 0)
 257     phi+=TMath::TwoPi();
 258   Double_t invmass = (gammai+gammaj).M();
 259   Double_t angle   = gammaj.Angle(gammai.Vect());
 260   Double_t e       = (gammai+gammaj).E();
 261
 262   //Fill histograms with no cuts applied.
 263   fhAnglePairNoCut->Fill(e,angle);
 264   fhInvMassPairNoCut->Fill(e,invmass);
 265
 266   //Cut on the aperture of the pair
 267   if(IsAngleInWindow(angle,e)){
 268     fhAnglePairOpeningAngleCut     ->Fill(e,angle);
 269     fhInvMassPairOpeningAngleCut->Fill(e,invmass);
 270     AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
 271
 272     //Cut on the invariant mass of the pair
 273     if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
 274       fhInvMassPairAllCut  ->Fill(e,invmass);
 275       fhAnglePairAllCut       ->Fill(e,angle);
 276       goodpair = kTRUE;
 277       AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
 278     }//(invmass>0.125) && (invmass<0.145)
 279   }//Opening angle cut
 280
 281   return goodpair;
 282
 283 }
 284
 285 //__________________________________________________________________
 286 void AliNeutralMesonSelection::Print(const Option_t * opt) const
 287 {
 288
 289   //Print some relevant parameters set for the analysis
 290   if(! opt)
 291     return;
 292
 293   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 294
 295   printf("mass : %f  \n", fM );
 296   printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
 297   printf("Angle selection param: \n");
 298   printf("p0 :     %f", fAngleMaxParam.At(0));
 299   printf("p1 :     %f", fAngleMaxParam.At(1));
 300   printf("p2 :     %f", fAngleMaxParam.At(2));
 301   printf("p3 :     %f", fAngleMaxParam.At(3));
 302
 303   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
 304
 305   if(fKeepNeutralMesonHistos){
 306         printf("Histograms: %3.1f < E  < %3.1f,  Nbin = %d\n", fHistoEMin,  fHistoEMax,  fHistoNEBins);
 307         printf("Histograms: %3.1f < pT < %3.1f,  Nbin = %d\n", fHistoPtMin,  fHistoPtMax,  fHistoNPtBins);
 308         printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
 309         printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
 310
 311   }
 312
 313 }