]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/GammaConv/AliConversionMesonCuts.cxx
Conversion Task able to run on AOD's, added different trigger selection, implemented...
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionMesonCuts.cxx
index 0a74603c85dc51deb7a5165fdb499e8a1418e3d1..b081698079322a53390fc1d10f374ca195f39c22 100644 (file)
 #include "AliESDEvent.h"
 #include "AliCentrality.h"
 #include "TList.h"
+#include "TPDGCode.h"
+#include "TDatabasePDG.h"
+#include "AliAODMCParticle.h"
+
 class iostream;
 
 using namespace std;
@@ -172,13 +176,15 @@ AliConversionMesonCuts::~AliConversionMesonCuts() {
 void AliConversionMesonCuts::InitCutHistograms(TString name){
 
    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
-
+   TH1::AddDirectory(kFALSE);
+   
    if(fHistograms != NULL){
       delete fHistograms;
       fHistograms=NULL;
    }
    if(fHistograms==NULL){
       fHistograms=new TList();
+      fHistograms->SetOwner(kTRUE);
       if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
    }
@@ -203,11 +209,12 @@ void AliConversionMesonCuts::InitCutHistograms(TString name){
    hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
    hMesonBGCuts->GetXaxis()->SetBinLabel(7,"out");
    fHistograms->Add(hMesonBGCuts);
+   
+   TH1::AddDirectory(kTRUE);
 }
 
-
 //________________________________________________________________________
-Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack){
+Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack, 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
    
@@ -218,9 +225,9 @@ Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *
       
       Double_t rapidity = 10.;
       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
-         rapidity=8.;
+         rapidity=8.-fRapidityShift;
       } else{
-         rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
+         rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
       }        
       
       // Rapidity Cut
@@ -242,61 +249,175 @@ Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *
    }
    return kFALSE;
 }
+//________________________________________________________________________
+Bool_t AliConversionMesonCuts::MesonIsSelectedAODMC(AliAODMCParticle *MCMother,TClonesArray *AODMCArray, 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(!AODMCArray)return kFALSE;
+
+   if(MCMother->GetPdgCode()==111 || MCMother->GetPdgCode()==221){
+      Double_t rMeson = sqrt( (MCMother->Xv()*MCMother->Xv()) + (MCMother->Yv()*MCMother->Yv()) ) ;
+      if(rMeson>fMaxR) return kFALSE; // cuts on distance from collision point
+
+      Double_t rapidity = 10.;
+      if(MCMother->E() - MCMother->Pz() == 0 || MCMother->E() + MCMother->Pz() == 0){
+         rapidity=8.-fRapidityShift;
+      } else{
+         rapidity = 0.5*(TMath::Log((MCMother->E()+MCMother->Pz()) / (MCMother->E()-MCMother->Pz())))-fRapidityShift;
+      }
+
+      // Rapidity Cut
+      if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
+
+      // Select only -> 2y decay channel
+      if(MCMother->GetNDaughters()!=2)return kFALSE;
 
+      for(Int_t i=0;i<2;i++){
+         AliAODMCParticle *MDaughter=static_cast<AliAODMCParticle*>(AODMCArray->At(MCMother->GetDaughter(i)));
+         // Is Daughter a Photon?
+         if(MDaughter->GetPdgCode()!=22)return kFALSE;
+         // Is Photon in Acceptance?
+         //   if(bMCDaughtersInAcceptance){
+         //    if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
+         //   }
+      }
+      return kTRUE;
+   }
+   return kFALSE;
+}
 //________________________________________________________________________
-Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack){
+Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelelectron, Int_t &labelpositron, 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( !fMCStack )return kFALSE;
        
-   if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
+   if( fMCMother->GetPdgCode() != 111 && fMCMother->GetPdgCode() != 221 ) return kFALSE;
                
-      if(fMCMother->R()>fMaxR) return kFALSE; // cuts on distance from collision point
+   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.;
-      }
-      else{
-         rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
-      }        
+   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;
+   if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 
       // Select only -> Dalitz decay channel
-      if(fMCMother->GetNDaughters()!=3)return kFALSE;
+   if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 
-      Int_t daughterPDGs[3] = {0,0,0};
-      Int_t index = 0;
+   TParticle *positron = 0x0;
+   TParticle *electron = 0x0;
+   TParticle    *gamma = 0x0;
+       
+   for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
+    
+     TParticle* temp = (TParticle*)fMCStack->Particle( index );
+               
+               switch( temp->GetPdgCode() ) {
+               case ::kPositron:
+                       positron      =  temp;
+                       labelpositron = index;
+                       break;
+               case ::kElectron:
+                       electron      =  temp;
+                       labelelectron = index;
+                       break;
+               case ::kGamma:
+                       gamma         =  temp;
+                       labelgamma    = index;
+                       break;
+               }
+  }
+  
+  if( positron && electron && 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
+   // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
+
+   if(!fMCStack)return kFALSE;
+        // if(fMCMother->GetPdgCode()==20443 ){
+        //      return kFALSE;
+        // }
+   if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
+                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;
 
-      //                iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter()
+                // Select only -> ChiC radiative (JPsi+gamma) decay channel
+      if(fMCMother->GetNDaughters()!=2)return kFALSE;
 
-      for(Int_t i=fMCMother->GetFirstDaughter(); i<= fMCMother->GetLastDaughter();i++){
-         TParticle *MDaughter=fMCStack->Particle(i);
-         // Is Daughter a Photon or an electron?
-         daughterPDGs[index]=MDaughter->GetPdgCode();
-         index++;
-      }
-      for (Int_t j=0;j<2;j++){
-
-         for (Int_t i=0;i<2;i++){
-            if (daughterPDGs[i] > daughterPDGs[i+1]){
-               Int_t interMed = daughterPDGs[i] ; 
-               daughterPDGs[i] = daughterPDGs[i+1];
-               daughterPDGs[i+1] = interMed;
-            }
-         }
-      }
-      if (daughterPDGs[0] == -11 && daughterPDGs[1] == 11 && daughterPDGs[2] == 22) return kTRUE;
+                       TParticle *jpsi         = 0x0;
+                       TParticle *gamma        = 0x0;
+                       TParticle *positron = 0x0;
+                       TParticle *electron = 0x0;
+
+                       Int_t labeljpsiChiC = -1;
+       
+                       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){                           
+                               
+                               TParticle* temp = (TParticle*)fMCStack->Particle( index );
+                               
+                               switch( temp->GetPdgCode() ) {
+                               case 443:
+                                       jpsi =  temp;
+                                       labeljpsiChiC = index;
+                                       break;
+                               case 22:
+                                       gamma    =  temp;
+                                       labelgammaChiC = index;
+                                       break;
+                               }
+                       }
+
+      if ( !jpsi || ! gamma) return kFALSE;
+                       if(jpsi->GetNDaughters()!=2)return kFALSE;
+                
+
+                       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){                             
+                               TParticle* temp = (TParticle*)fMCStack->Particle( index );
+                               switch( temp->GetPdgCode() ) {
+                               case -11:
+                                       electron =  temp;
+                                       labelelectronChiC = index;
+                                       break;
+                               case 11:
+                                       positron =  temp;
+                                       labelpositronChiC = index;
+                                       break;
+                               }
+                       }
+                       if( !electron || !positron) return kFALSE;
+                       if( positron && electron && gamma) return kTRUE;
    }
    return kFALSE;
 }
 
 
+
 ///________________________________________________________________________
-Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal)
+Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
 {
    // Selection of reconstructed Meson candidates
    // Use flag IsSignal in order to fill Fill different
@@ -320,7 +441,7 @@ Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_
    else{
       // PseudoRapidity Cut --> But we cut on Rapidity !!!
       cutIndex++;
-      if(abs(pi0->Rapidity())>fRapidityCutMeson){
+      if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
          if(hist)hist->Fill(cutIndex);
          return kFALSE;
       }
@@ -361,13 +482,10 @@ Bool_t AliConversionMesonCuts::UpdateCutString() {
    ///Update the cut string (if it has been created yet)
 
    if(fCutString && fCutString->GetString().Length() == kNCuts) {
-      //         cout << "Updating cut id in spot number " << cutID << " to " << value << endl;
       fCutString->SetString(GetCutNumber());
    } else {
-      //         cout << "fCutString not yet initialized, will not be updated" << endl;
       return kFALSE;
    }
-   //   cout << fCutString->GetString().Data() << endl;
    return kTRUE;
 }
 
@@ -651,7 +769,15 @@ Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
    case 6:  //
       fRapidityCutMeson   = 0.75;
       break;
-
+   case 7:  //
+      fRapidityCutMeson   = 0.3;
+      break;
+   case 8:  //
+      fRapidityCutMeson   = 0.35;
+      break;
+   case 9:  //
+      fRapidityCutMeson   = 0.4;
+      break;
    default:
       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
       return kFALSE;