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.1.2.1  2007/07/26 10:32:09  schutz
  21  * new analysis classes in the the new analysis framework
  22  *
  23  *
  24  */
  25
  26 //_________________________________________________________________________
  27 // Class that contains methods to select candidate pairs to neutral meson
  28 //-- Author: Gustavo Conesa (INFN-LNF)
  29
  30 // --- ROOT system ---
  31 #include <TParticle.h>
  32 #include <TLorentzVector.h>
  33 #include <TH2.h>
  34 #include <TList.h>
  35
  36 //---- AliRoot system ----
  37 #include "AliNeutralMesonSelection.h"
  38 #include "Riostream.h"
  39 #include "AliLog.h"
  40
  41 ClassImp(AliNeutralMesonSelection)
  42
  43
  44 //____________________________________________________________________________
  45   AliNeutralMesonSelection::AliNeutralMesonSelection() :
  46     TObject(), fSelect(0), fM(0),
  47     fInvMassMaxCut(0.), fInvMassMinCut(0.),
  48     fAngleMaxParam(),  fMinPt(0),
  49     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
  50     fRatioMaxCut(0), fRatioMinCut(0),  fKeepNeutralMesonHistos(0),
  51     fhAnglePairNoCut(0),  fhAnglePairCorrelationCut(0),
  52     fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0),
  53     fhInvMassPairNoCut(0),   fhInvMassPairCorrelationCut(0),
  54     fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0)
  55 {
  56   //Default Ctor
  57
  58   //Initialize parameters
  59
  60   // kGammaHadron and kGammaJet
  61   fAngleMaxParam.Set(4) ;
  62   fAngleMaxParam.Reset(0.);
  63
  64   //Initialize parameters
  65   InitParameters();
  66 }
  67
  68 //____________________________________________________________________________
  69 AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
  70   TObject(),
  71   fSelect(g.fSelect), fM(g.fM),
  72   fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
  73   fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
  74   fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
  75   fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut),
  76   fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
  77   fhAnglePairNoCut(g. fhAnglePairNoCut),
  78   fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut),
  79   fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
  80   fhAnglePairAllCut(g. fhAnglePairAllCut),
  81   fhInvMassPairNoCut(g.fhInvMassPairNoCut),
  82   fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut),
  83   fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
  84   fhInvMassPairAllCut(g.fhInvMassPairAllCut)
  85 {
  86   // cpy ctor
  87 }
  88
  89 //_________________________________________________________________________
  90 AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
  91 {
  92   // assignment operator
  93
  94   if(this == &source)return *this;
  95   ((TObject *)this)->operator=(source);
  96
  97   fSelect = source.fSelect ;
  98   fM = source.fM ;
  99   fInvMassMaxCut = source.fInvMassMaxCut ;
 100   fInvMassMinCut = source.fInvMassMinCut ;
 101   fAngleMaxParam = source.fAngleMaxParam ;
 102   fMinPt = source.fMinPt ;
 103   fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
 104   fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
 105   fRatioMaxCut = source.fRatioMaxCut ;
 106   fRatioMinCut = source.fRatioMinCut ;
 107   fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
 108
 109   fhAnglePairNoCut = source. fhAnglePairNoCut ;
 110   fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ;
 111   fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
 112   fhAnglePairAllCut = source. fhAnglePairAllCut ;
 113   fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
 114   fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ;
 115   fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
 116   fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
 117
 118   return *this;
 119
 120 }
 121
 122 //____________________________________________________________________________
 123 AliNeutralMesonSelection::~AliNeutralMesonSelection()
 124 {
 125   // Remove all pointers
 126
 127   delete fhAnglePairNoCut  ;
 128   delete fhAnglePairCorrelationCut  ;
 129   delete fhAnglePairOpeningAngleCut   ;
 130   delete fhAnglePairAllCut   ;
 131   delete fhInvMassPairNoCut    ;
 132   delete fhInvMassPairCorrelationCut    ;
 133   delete fhInvMassPairOpeningAngleCut  ;
 134   delete fhInvMassPairAllCut   ;
 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) {
 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)
 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 }