Merge remote-tracking branch 'origin/master' into flatdev

[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionMesonCuts.cxx

diff --git a/PWGGA/GammaConv/AliConversionMesonCuts.cxx b/PWGGA/GammaConv/AliConversionMesonCuts.cxx
index a1c12f4748143d6b1d3fdf42cafc499ff746424c..ab13e8749ad7d7a3123996500e8af950c242fe59 100644 (file)
--- a/PWGGA/GammaConv/AliConversionMesonCuts.cxx
+++ b/PWGGA/GammaConv/AliConversionMesonCuts.cxx
@@ -94,6 +94,8 @@ AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *titl
    fPSigSmearingCte(0),
    fBrem(NULL),
    fRandom(0),
+   fFAlphaCut(0),
+   fAlphaPtDepCut(kFALSE),
    fElectronLabelArraySize(500),
    fElectronLabelArray(NULL),
    fDCAGammaGammaCut(1000),
@@ -149,6 +151,8 @@ AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref
    fPSigSmearingCte(ref.fPSigSmearingCte),
    fBrem(NULL),
    fRandom(ref.fRandom),
+   fFAlphaCut(NULL),
+   fAlphaPtDepCut(ref.fAlphaPtDepCut),
    fElectronLabelArraySize(ref.fElectronLabelArraySize),
    fElectronLabelArray(NULL),
    fDCAGammaGammaCut(ref.fDCAGammaGammaCut),
@@ -334,6 +338,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelectedAODMC(AliAODMCParticle *MCMother,T
    }
    return kFALSE;
 }
+
 //________________________________________________________________________
 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelelectron, Int_t &labelpositron, Int_t &labelgamma, Double_t fRapidityShift){
 
@@ -387,9 +392,63 @@ Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliS
 
    if( positron && electron && gamma) return kTRUE;
    return kFALSE;
+}
+
+//________________________________________________________________________
+Bool_t AliConversionMesonCuts::MesonIsSelectedMCEtaPiPlPiMiGamma(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelGamma, Double_t fRapidityShift){
+
+       // Returns true for all pions within acceptance cuts for decay into 2 photons
+       // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
+
+       if( !fMCStack )return kFALSE;
+
+       if( fMCMother->GetPdgCode() != 221 ) return kFALSE;
+
+       if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
+
+       Double_t rapidity = 10.;
+
+       if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
+               rapidity=8.-fRapidityShift;
+       }
+       else{
+               rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
+       }
+
+       // Rapidity Cut
+       if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
+
+       // Select only -> Dalitz decay channel
+       if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 
+       TParticle *posPion = 0x0;
+       TParticle *negPion = 0x0;
+       TParticle    *gamma = 0x0;
 
+       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
+
+               TParticle* temp = (TParticle*)fMCStack->Particle( index );
+
+               switch( temp->GetPdgCode() ) {
+               case 211:
+                       posPion      =  temp;
+                       labelPosPion = index;
+                       break;
+               case -211:
+                       negPion      =  temp;
+                       labelNegPion = index;
+                       break;
+               case ::kGamma:
+                       gamma         =  temp;
+                       labelGamma    = index;
+                       break;
+               }
+       }
+
+       if( posPion && negPion && gamma) return kTRUE;
+       return kFALSE;
 }
+
 //________________________________________________________________________
 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
@@ -485,6 +544,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
       if(hist)hist->Fill(cutIndex);
       cutIndex++;
+//       cout << "undefined rapidity" << endl;
       return kFALSE;
    }
    else{
@@ -492,6 +552,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
       cutIndex++;
       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
          if(hist)hist->Fill(cutIndex);
+//              cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
          return kFALSE;
       }
    }
@@ -500,13 +561,21 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
    // Opening Angle Cut
    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
    if( pi0->GetOpeningAngle() < fOpeningAngle){
+//       cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
    cutIndex++;
 
+   if ( fAlphaPtDepCut == kTRUE ) {
+ 
+       fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
+   }
+   
+   
    // Alpha Max Cut
    if(pi0->GetAlpha()>fAlphaCutMeson){
+//        cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
@@ -514,6 +583,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
 
    // Alpha Min Cut
    if(pi0->GetAlpha()<fAlphaMinCutMeson){
+//       cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
@@ -523,12 +593,14 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
 
    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
+//       cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
    cutIndex++;
 
    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
+//        cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
@@ -538,6 +610,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
 
    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
+//       cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl; 
       if(hist)hist->Fill(cutIndex);
       return kFALSE;
    }
@@ -796,42 +869,64 @@ Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
    case 0:     // 0- 0.7
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.7;
-      break;
-   case 1:     // 0-0.5
-      fAlphaMinCutMeson         = 0.0;
-      fAlphaCutMeson    = 0.5;
-      break;
-   case 2:     // 0.5-1
-      fAlphaMinCutMeson         = 0.5;
-      fAlphaCutMeson    = 1.;
+      fAlphaPtDepCut = kFALSE;
+      break;
+   case 1:     // Updated 31 October 2013 before 0.0 - 0.5
+      if( fFAlphaCut ) delete fFAlphaCut;
+      fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
+      fFAlphaCut->SetParameter(0,0.7);
+      fFAlphaCut->SetParameter(1,1.2);
+      fAlphaMinCutMeson         =  0.0;
+      fAlphaCutMeson    = -1.0;
+      fAlphaPtDepCut = kTRUE;
+      break;
+   case 2:     // Updated 31 October 2013 before 0.5-1  
+      if( fFAlphaCut ) delete fFAlphaCut;
+      fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
+      fFAlphaCut->SetParameter(0,0.8);
+      fFAlphaCut->SetParameter(1,1.2);
+      fAlphaMinCutMeson         =  0.0;
+      fAlphaCutMeson    = -1.0;
+      fAlphaPtDepCut = kTRUE;
       break;
    case 3:     // 0.0-1
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 1.;
+      fAlphaPtDepCut = kFALSE;
       break;
    case 4:     // 0-0.65
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.65;
+      fAlphaPtDepCut = kFALSE;
       break;
    case 5:     // 0-0.75
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.75;
+      fAlphaPtDepCut = kFALSE;
       break;
    case 6:     // 0-0.8
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.8;
+      fAlphaPtDepCut = kFALSE;
       break;
    case 7:     // 0.0-0.85
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.85;
+      fAlphaPtDepCut = kFALSE;
       break;
    case 8:     // 0.0-0.6
       fAlphaMinCutMeson         = 0.0;
       fAlphaCutMeson    = 0.6;
+      fAlphaPtDepCut = kFALSE;
       break;
-   case 9:     // 0.0-0.3
-      fAlphaMinCutMeson         = 0.0;
-      fAlphaCutMeson    = 0.3;
+   case 9: // Updated 11 November 2013 before 0.0 - 0.3
+      if( fFAlphaCut ) delete fFAlphaCut;
+      fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
+      fFAlphaCut->SetParameter(0,0.65);
+      fFAlphaCut->SetParameter(1,1.2);
+      fAlphaMinCutMeson  =  0.0;
+      fAlphaCutMeson     = -1.0;
+      fAlphaPtDepCut = kTRUE;
       break;
    default:
       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
@@ -844,8 +939,8 @@ Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
    // Set Cut
    switch(RapidityMesonCut){
-   case 0:  //
-      fRapidityCutMeson   = 0.9;
+   case 0:  // changed from 0.9 to 1.35
+      fRapidityCutMeson   = 1.35;
       break;
    case 1:  //
       fRapidityCutMeson   = 0.8;