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