PWG4/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$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.3  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  * Revision 1.2  2007/08/17 12:40:04  schutz
  24  * New analysis classes by Gustavo Conesa
  25  *
  26  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  27  * new analysis classes in the the new analysis framework
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Class that contains methods to select candidate pairs to neutral meson
  34 // 2 main selections, invariant mass around pi0 (also any other mass),
  35 // apperture angle to distinguish from combinatorial.
  36 // There is a 3rd cut based on the gamma correlation on phi or pt.
  37 //-- Author: Gustavo Conesa (INFN-LNF)
  38
  39 // --- ROOT system ---
  40 #include <TParticle.h>
  41 #include <TLorentzVector.h>
  42 #include <TH2.h>
  43 #include <TList.h>
  44 #include <TArrayD.h>
  45
  46 //---- AliRoot system ----
  47 #include "AliNeutralMesonSelection.h"
  48 #include "AliLog.h"
  49
  50 ClassImp(AliNeutralMesonSelection)
  51
  52
  53 //____________________________________________________________________________
  54   AliNeutralMesonSelection::AliNeutralMesonSelection() :
  55     TObject(), fSelect(0), fM(0),
  56     fInvMassMaxCut(0.), fInvMassMinCut(0.),
  57     fAngleMaxParam(),  fMinPt(0),
  58     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
  59     fRatioMaxCut(0), fRatioMinCut(0),  fKeepNeutralMesonHistos(0),
  60     fhAnglePairNoCut(0),  fhAnglePairCorrelationCut(0),
  61     fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0),
  62     fhInvMassPairNoCut(0),   fhInvMassPairCorrelationCut(0),
  63     fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0)
  64 {
  65   //Default Ctor
  66
  67   //Initialize parameters
  68
  69   // kGammaHadron and kGammaJet
  70   fAngleMaxParam.Set(4) ;
  71   fAngleMaxParam.Reset(0.);
  72
  73   //Initialize parameters
  74   InitParameters();
  75 }
  76
  77 //____________________________________________________________________________
  78 AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
  79   TObject(),
  80   fSelect(g.fSelect), fM(g.fM),
  81   fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
  82   fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
  83   fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
  84   fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut),
  85   fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
  86   fhAnglePairNoCut(g. fhAnglePairNoCut),
  87   fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut),
  88   fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
  89   fhAnglePairAllCut(g. fhAnglePairAllCut),
  90   fhInvMassPairNoCut(g.fhInvMassPairNoCut),
  91   fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut),
  92   fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
  93   fhInvMassPairAllCut(g.fhInvMassPairAllCut)
  94 {
  95   // cpy ctor
  96 }
  97
  98 //_________________________________________________________________________
  99 AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
 100 {
 101   // assignment operator
 102
 103   if(this == &source)return *this;
 104   ((TObject *)this)->operator=(source);
 105
 106   fSelect = source.fSelect ;
 107   fM = source.fM ;
 108   fInvMassMaxCut = source.fInvMassMaxCut ;
 109   fInvMassMinCut = source.fInvMassMinCut ;
 110   fAngleMaxParam = source.fAngleMaxParam ;
 111   fMinPt = source.fMinPt ;
 112   fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
 113   fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
 114   fRatioMaxCut = source.fRatioMaxCut ;
 115   fRatioMinCut = source.fRatioMinCut ;
 116   fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
 117
 118   fhAnglePairNoCut = source. fhAnglePairNoCut ;
 119   fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ;
 120   fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
 121   fhAnglePairAllCut = source. fhAnglePairAllCut ;
 122   fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
 123   fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ;
 124   fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
 125   fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
 126
 127   return *this;
 128
 129 }
 130
 131 //____________________________________________________________________________
 132 AliNeutralMesonSelection::~AliNeutralMesonSelection()
 133 {
 134  // Remove all pointers except analysis output pointers.
 135
 136 }
 137
 138
 139
 140 //________________________________________________________________________
 141 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
 142 {
 143
 144   // Create histograms to be saved in output file and
 145   // store them in outputContainer
 146   TList * outputContainer = new TList() ;
 147   outputContainer->SetName("MesonDecayHistos") ;
 148
 149   fhAnglePairNoCut  = new TH2F
 150     ("AnglePairNoCut",
 151      "Angle between all #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2);
 152   fhAnglePairNoCut->SetYTitle("Angle (rad)");
 153   fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 154
 155   fhAnglePairOpeningAngleCut  = new TH2F
 156     ("AnglePairOpeningAngleCut",
 157      "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
 158      ,200,0,50,200,0,0.2);
 159   fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
 160   fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 161
 162   fhAnglePairAllCut  = new TH2F
 163     ("AnglePairAllCut",
 164      "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
 165      ,200,0,50,200,0,0.2);
 166   fhAnglePairAllCut->SetYTitle("Angle (rad)");
 167   fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 168
 169   //
 170   fhInvMassPairNoCut  = new TH2F
 171     ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
 172      120,0,120,360,0,0.5);
 173   fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 174   fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 175
 176   fhInvMassPairOpeningAngleCut  = new TH2F
 177     ("InvMassPairOpeningAngleCut",
 178      "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
 179      120,0,120,360,0,0.5);
 180   fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 181   fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
 182
 183   fhInvMassPairAllCut  = new TH2F
 184     ("InvMassPairAllCut",
 185      "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
 186      120,0,120,360,0,0.5);
 187   fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 188   fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
 189
 190   fhAnglePairCorrelationCut  = new TH2F
 191     ("AnglePairCorrelationCut",
 192      "Angle between correlated #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2);
 193   fhAnglePairCorrelationCut->SetYTitle("Angle (rad)");
 194   fhAnglePairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 195
 196   fhInvMassPairCorrelationCut  = new TH2F
 197     ("InvMassPairCorrelationCut","Invariant Mass of correlated #gamma pair vs E_{#pi^{0}}",
 198      120,0,120,360,0,0.5);
 199   fhInvMassPairCorrelationCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 200   fhInvMassPairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
 201
 202   outputContainer->Add(fhAnglePairNoCut) ;
 203   outputContainer->Add(fhAnglePairOpeningAngleCut) ;
 204   outputContainer->Add(fhAnglePairAllCut) ;
 205
 206   outputContainer->Add(fhInvMassPairNoCut) ;
 207   outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
 208   outputContainer->Add(fhInvMassPairAllCut) ;
 209
 210   outputContainer->Add(fhAnglePairCorrelationCut) ;
 211   outputContainer->Add(fhInvMassPairCorrelationCut) ;
 212
 213   return outputContainer;
 214 }
 215
 216  //____________________________________________________________________________
 217 void AliNeutralMesonSelection::InitParameters()
 218 {
 219
 220   //Initialize the parameters of the analysis.
 221   fKeepNeutralMesonHistos = kTRUE ;
 222
 223   //-------------kHadron, kJetLeadCone-----------------
 224   fAngleMaxParam.Set(4) ;
 225   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
 226   fAngleMaxParam.AddAt(-0.25,1) ;
 227   fAngleMaxParam.AddAt(0.025,2) ;
 228   fAngleMaxParam.AddAt(-2e-4,3) ;
 229
 230   fInvMassMaxCut  = 0.16 ;
 231   fInvMassMinCut  = 0.11 ;
 232
 233   fM = 0.1349766;//neutralMeson mass
 234
 235   fMinPt = 0.   ;
 236   fDeltaPhiMaxCut      = 4.5;
 237   fDeltaPhiMinCut      = 1.5 ;
 238   fRatioMaxCut    = 1.0 ;
 239   fRatioMinCut    = 0.1 ;
 240 }
 241
 242 //__________________________________________________________________________-
 243 Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
 244   //Check if the opening angle of the candidate pairs is inside
 245   //our selection windowd
 246
 247   Bool_t result = kFALSE;
 248   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
 249     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
 250   Double_t arg = (e*e-2*fM*fM)/(e*e);
 251   Double_t min = 100. ;
 252   if(arg>0.)
 253     min = TMath::ACos(arg);
 254
 255   if((angle<max)&&(angle>=min))
 256     result = kTRUE;
 257
 258   return result;
 259 }
 260
 261 //____________________________________________________________________________
 262 Bool_t  AliNeutralMesonSelection::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi)  const
 263 {
 264   //Select pair if delta
 265   Bool_t cut = kFALSE ;
 266
 267   if(fSelect == kNoSelectPhiPt) cut = kTRUE ;
 268   else if((phig-phi) > fDeltaPhiMinCut && ((phig-phi) < fDeltaPhiMaxCut)){
 269     //Cut on pt
 270     if((fSelect == kSelectPhiPtRatio && ptg > 0. && pt/ptg  > fRatioMinCut &&  pt/ptg  < fRatioMaxCut) ||
 271         (fSelect == kSelectPhiMinPt && pt > fMinPt)  )  cut = kTRUE ;
 272   }
 273   else cut = kFALSE ;
 274
 275   return cut ;
 276
 277 }
 278
 279 //____________________________________________________________________________
 280 Bool_t  AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector gammai, TLorentzVector gammaj)
 281 {
 282
 283   //Search for the neutral pion within selection cuts
 284
 285   Double_t ptg = pGamma->Pt();
 286   Double_t phig = pGamma->Phi() ;
 287   Bool_t goodpair = kFALSE ;
 288
 289   Double_t pt  = (gammai+gammaj).Pt();
 290   Double_t phi = (gammai+gammaj).Phi();
 291   if(phi < 0)
 292     phi+=TMath::TwoPi();
 293   Double_t invmass = (gammai+gammaj).M();
 294   Double_t angle   = gammaj.Angle(gammai.Vect());
 295   Double_t e       = (gammai+gammaj).E();
 296
 297   //Fill histograms with no cuts applied.
 298   fhAnglePairNoCut->Fill(e,angle);
 299   fhInvMassPairNoCut->Fill(e,invmass);
 300
 301   //Cut on phig-phi meson
 302   if(CutPtPhi(ptg, phig, pt, phi)){
 303
 304     fhAnglePairCorrelationCut     ->Fill(e,angle);
 305     fhInvMassPairCorrelationCut->Fill(e,invmass);
 306
 307     //Cut on the aperture of the pair
 308     if(IsAngleInWindow(angle,e)){
 309       fhAnglePairOpeningAngleCut     ->Fill(e,angle);
 310       fhInvMassPairOpeningAngleCut->Fill(e,invmass);
 311       AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
 312
 313       //Cut on the invariant mass of the pair
 314       if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
 315         fhInvMassPairAllCut  ->Fill(e,invmass);
 316         fhAnglePairAllCut       ->Fill(e,angle);
 317         goodpair = kTRUE;
 318         AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
 319       }//(invmass>0.125) && (invmass<0.145)
 320     }//Opening angle cut
 321   } // cut on pt and phi
 322
 323
 324   return goodpair;
 325
 326 }
 327
 328 //__________________________________________________________________
 329 void AliNeutralMesonSelection::Print(const Option_t * opt) const
 330 {
 331
 332   //Print some relevant parameters set for the analysis
 333   if(! opt)
 334     return;
 335
 336   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 337
 338   printf("mass : %f  \n", fM );
 339   printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
 340   printf("Angle selection param: \n");
 341   printf("p0 :     %f", fAngleMaxParam.At(0));
 342   printf("p1 :     %f", fAngleMaxParam.At(1));
 343   printf("p2 :     %f", fAngleMaxParam.At(2));
 344   printf("p3 :     %f", fAngleMaxParam.At(3));
 345
 346   printf("pT meson       >    %f\n", fMinPt) ;
 347   printf("Phi gamma-meson      <     %f\n", fDeltaPhiMaxCut) ;
 348   printf("Phi gamma-meson      >     %f\n", fDeltaPhiMinCut) ;
 349   printf("pT meson / pT Gamma             <     %f\n", fRatioMaxCut) ;
 350   printf("pT meson / pT Gamma             >     %f\n", fRatioMinCut) ;
 351   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
 352
 353 }