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