new Muon-Hadron correlation task in DHC (Constantin Loizides <cloizides@lbl.gov>...
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 24 Jan 2013 07:22:42 +0000 (07:22 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 24 Jan 2013 07:22:42 +0000 (07:22 +0000)
PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg
PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskDhc.C [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskMuonEffMC.C [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx [new file with mode: 0644]
PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.h [new file with mode: 0644]
PWGCF/PWGCFCorrelationsDPhiLinkDef.h

index 4df0efed0d12625bc416e855b952718bb1458a6b..5987ba237b8e1f75d4eb6523fe942e44ac7e92ad 100644 (file)
@@ -37,6 +37,7 @@ set ( SRCS
     Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx 
     Correlations/DPhi/FourierDecomposition/AliPool.cxx
     Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
+    Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx
     Correlations/DPhi/AliLeadingV0Correlation.cxx 
     Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.cxx 
     )
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskDhc.C b/PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskDhc.C
new file mode 100644 (file)
index 0000000..a6e77f7
--- /dev/null
@@ -0,0 +1,87 @@
+AliDhcTask *AddTaskDhc(Int_t iAna = 1)
+{
+  const char *nTracks     = "PicoTracks";
+  const char *inputTracks = "HybridTracks";
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskDhc", "No analysis manager found.");
+    return;
+  }
+
+  // Track Cuts
+  gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalEsdTpcTrack.C");
+  AliEmcalEsdTpcTrackTask *hybTask = AddTaskEmcalEsdTpcTrack(inputTracks,"Hybrid_LHC11h");
+
+  // Pico Tracks
+  gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPicoTrackMaker.C");
+  AliEmcalPicoTrackMaker *pTrackTask = AddTaskEmcalPicoTrackMaker(nTracks, inputTracks, "LHC11h");
+
+  // Binning
+  Double_t arPt[5] = {0.5, 1.0, 2.0, 4.0};
+  TAxis *axPt = new TAxis(3,arPt);
+  Double_t arCent[5] = {0.0, 20.0, 40.0, 60.0, 100.0};
+  TAxis *axCent = new TAxis(4,arCent);
+  TAxis *axZvtx = new TAxis(1,-10.0,10.0);
+  Double_t arCentMix[9] = {0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0};
+  TAxis *axCentMix = new TAxis(8,arCentMix);
+  TAxis *axZvtxMix = new TAxis(8,-10.0,10.0);
+
+  // Efficiency correction files
+  TFile *fiHEff  = 0;
+  TFile *fiMuEff = 0;
+  fiHEff = TFile::Open("alien:///alice/cern.ch/user/t/tschuste/correction_hybrid_nulled.root","OLD");
+  if (!fiHEff){
+    cout << "Requested file:" << fiHEff << " was not opened. ABORT." << endl;
+    return;
+  }
+  THnF* hHEff = (THnF*) fiHEff->Get("correction");
+  
+  fiMuEff = TFile::Open("alien:///alice/cern.ch/user/t/tschuste/correction_muon.root","OLD");
+  if (!fiMuEff){
+    cout << "Requested file:" << fiMuEff << " was not opened. ABORT." << endl;
+    return;
+  }
+  THnF* hMuEff = (THnF*) fiMuEff->Get("correction");
+  
+
+  AliDhcTask *dhcTask_NW = new AliDhcTask("Task_tschuste_Dhc_NW");
+  if (iAna==1) { // h-h
+    Int_t nDetaBins = 40;
+    Int_t nDPhiBins = 72;
+    dhcTask_NW->SetAnaMode(AliDhcTask::kHH);
+    dhcTask_NW->SetHEffT(hHEff);
+    dhcTask_NW->SetHEffA(hHEff);
+    dhcTask_NW->SetEtaMax(1.2);
+  } else if (iAna==2) { // mu-h
+    Int_t nDetaBins = 100;
+    Int_t nDPhiBins = 36;
+    dhcTask_NW->SetAnaMode(AliDhcTask::kMuH);
+    dhcTask_NW->SetHEffT(hMuEff);
+    dhcTask_NW->SetHEffA(hHEff);
+    dhcTask_NW->SetEtaMax(5.0);
+  }
+  dhcTask_NW->SetTracksName(nTracks);
+  dhcTask_NW->SetDoWeights(kFALSE);
+  dhcTask_NW->SetCentMethod("V0M");
+  dhcTask_NW->SetDEtaDPhiBins(nDetaBins,nDPhiBins);
+  dhcTask_NW->SetPtTBins(axPt);
+  dhcTask_NW->SetPtABins(axPt);
+  dhcTask_NW->SetCentBins(axCent);
+  dhcTask_NW->SetZVtxBins(axZvtx);
+  dhcTask_NW->SetCentMixBins(axCentMix);
+  dhcTask_NW->SetZVtxMixBins(axZvtxMix);
+  dhcTask_NW->SelectCollisionCandidates(AliVEvent::kINT7);
+  dhcTask_NW->SetVerbosity(10);
+  mgr->AddTask(dhcTask_NW);
+  
+  AliAnalysisDataContainer *co_Dhc_NW = mgr->CreateContainer("Cont_tschuste_DhcAna_NW", 
+                                                             TList::Class(),
+                                                             AliAnalysisManager::kOutputContainer,
+                                                             Form("%s:PWGCF.outDhc_NW_%d.root", AliAnalysisManager::GetCommonFileName(), iAna));
+  mgr->ConnectInput(dhcTask_NW,0,mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(dhcTask_NW,1,co_Dhc_NW);
+  
+  return dhcTask_NW;
+
+}
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskMuonEffMC.C b/PWGCF/Correlations/DPhi/FourierDecomposition/AddTaskMuonEffMC.C
new file mode 100644 (file)
index 0000000..8e36286
--- /dev/null
@@ -0,0 +1,68 @@
+
+AliMuonEffMC* AddTaskMuonEFFMC(Bool_t MDProcess = kTRUE,
+                              TString centralityEstimator = "V0M",
+                              const Int_t NEtaBins = 15,
+                              const Int_t NpTBins = 100,
+                              const Int_t NCentBins = 1,
+                              const Int_t NZvtxBins = 10,
+                              const Int_t NPhiBins = 12,
+                              const char* outputFileName = 0,
+                              const char* folderName = "Muon_TrkEff")
+{
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskSOH", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskSOH", "This task requires an input event handler");
+    return NULL;
+  }
+
+  TString analysisType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  if(dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())) analysisType = "MC";
+
+  Bool_t IsMc = kTRUE;
+  if (analysisType !="MC") IsMc = kFALSE;
+
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  AliMuonEffMC *MuonEff = new AliMuonEffMC("MuonEffMC");
+  MuonEff->SetMcAna(IsMc);
+  MuonEff->SetMDProcess(MDProcess);
+  MuonEff->SetCentEstimator(centralityEstimator);
+  MuonEff->SetNEtaBins(NEtaBins);
+  MuonEff->SetNpTBins(NpTBins);
+  MuonEff->SetNCentBins(NCentBins);
+  MuonEff->SetNZvtxBins(NZvtxBins);
+  MuonEff->SetNPhiBins(NPhiBins);
+  MuonEff->SelectCollisionCandidates(AliVEvent::kAnyINT);
+
+  // Add task(s)
+  mgr->AddTask(MuonEff); 
+
+  // Create containers for input/output
+  if (!outputFileName) 
+    outputFileName = AliAnalysisManager::GetCommonFileName();
+
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+  AliAnalysisDataContainer *coutputpt = mgr->CreateContainer(Form("MuonEff_%s",centralityEstimator.Data()), 
+                                                             TList::Class(), 
+                                                             AliAnalysisManager::kOutputContainer, 
+                                                             Form("%s:%s", outputFileName, folderName)););
+
+  // Connect input/output
+  mgr->ConnectInput(MuonEff, 0, cinput);
+  mgr->ConnectOutput(MuonEff, 1, coutputpt);
+
+  return MuonEff;
+}
index b6a1423cd1815da09bec009aed81aa11914911ce..9572d6c437def088156df64b0f5f5404a9637551 100644 (file)
@@ -716,6 +716,7 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
 
     Float_t pta  = a.Pt();
     Float_t etaa = a.Eta();
+    Float_t phia = a.Phi();
     Int_t abin = fHPtTrg->FindBin(pta);
     if (fHPtTrg->IsBinOverflow(abin) || fHPtTrg->IsBinUnderflow(abin))
       continue;
@@ -724,19 +725,23 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
       continue;
 
     // efficiency weighting
-    Int_t effBin[4];
     Double_t effWtT = 1.0;
     if (fHEffT) {
       // trigger particle
-      effBin[0] = fHEffT->GetAxis(0)->FindBin(etaa);
-      effBin[1] = fHEffT->GetAxis(1)->FindBin(pta);
-      effBin[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
-      effBin[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
-      effWtT = fHEffT->GetBinContent(effBin);
+      const Int_t nEffDimT = fHEffT->GetNdimensions();
+      Int_t effBinT[nEffDimT];
+      effBinT[0] = fHEffT->GetAxis(0)->FindBin(etaa);
+      effBinT[1] = fHEffT->GetAxis(1)->FindBin(pta);
+      effBinT[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
+      effBinT[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
+      if (nEffDimT>4) {
+        effBinT[4] = fHEffT->GetAxis(4)->FindBin(phia);
+      }
+      effWtT = fHEffT->GetBinContent(effBinT);
     }
     
     if (pairing == kSameEvt) {
-      fHTrk->Fill(a.Phi(),etaa);
+      fHTrk->Fill(phia,etaa);
       fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
     } else {
       fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
@@ -752,6 +757,7 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
       
       Float_t ptb  = b.Pt();
       Float_t etab = b.Eta();
+      Float_t phib = b.Phi();
       if (pta < ptb)
            continue;
 
@@ -762,7 +768,7 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
       if (etab<fEtaALo || etab>fEtaAHi)
         continue;
 
-      Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
+      Float_t dphi = DeltaPhi(phia, phib);
       Float_t deta = etaa - etab;
 
       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
@@ -784,13 +790,17 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
       }
       if (fHEffA) {
         // associated particle
-        effBin[0] = fHEffA->GetAxis(0)->FindBin(etab);
-        effBin[1] = fHEffA->GetAxis(1)->FindBin(ptb);
-        effBin[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
-        effBin[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
-        weight *= fHEffA->GetBinContent(effBin);
+        const Int_t nEffDimA = fHEffA->GetNdimensions();
+        Int_t effBinA[nEffDimA];
+        effBinA[0] = fHEffA->GetAxis(0)->FindBin(etab);
+        effBinA[1] = fHEffA->GetAxis(1)->FindBin(ptb);
+        effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
+        effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
+        if (nEffDimA>4) {
+          effBinA[4] = fHEffA->GetAxis(4)->FindBin(phib);
+        }
+        weight *= fHEffA->GetBinContent(effBinA);
       }
-
       hist[index]->Fill(dphi,deta,weight);
       bCountTrg = kTRUE;
 
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx b/PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx
new file mode 100644 (file)
index 0000000..f5b04a6
--- /dev/null
@@ -0,0 +1,435 @@
+// MUON track QA referring AliMuonEffMC.cxx
+// Author : Saehanseul Oh
+
+#include "AliMuonEffMC.h"
+
+#include <TList.h>
+#include <TH1D.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THn.h>
+#include <TChain.h>
+#include <TFile.h>
+
+#include "AliAnalysisManager.h"
+#include "AliAnalysisTask.h"
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDVertex.h"
+#include "AliAODVertex.h"
+#include "AliCentrality.h"
+#include "AliVParticle.h"
+#include "AliMCParticle.h"
+#include "AliAnalysisHelperJetTasks.h"
+#include "AliMCEvent.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliMuonEffMC)
+
+//________________________________________________________________________
+AliMuonEffMC::AliMuonEffMC() :
+  AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
+  fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
+  fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
+  fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
+{
+  // Constructor
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+
+  for(Int_t i=0; i<4; i++)
+  {
+    fHMuMotherGenPt[i] = 0x0;
+    fHMuMotherRecPt[i] = 0x0;
+    fHMuMotherGenPhi[i] = 0x0;
+    fHMuMotherRecPhi[i] = 0x0;
+    fHMuMotherGenEta[i] = 0x0;
+    fHMuMotherRecEta[i] = 0x0;
+    fHMuDCA[i] = 0x0;
+  }
+}
+
+//________________________________________________________________________
+AliMuonEffMC::AliMuonEffMC(const char *name) :
+  AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
+  fHEventStat(0), fHEvt(0x0),  fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
+  fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
+  fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
+{
+  // Constructor
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+
+  for(Int_t i=0; i<4; i++)
+  {
+    fHMuMotherGenPt[i] = 0x0;
+    fHMuMotherRecPt[i] = 0x0;
+    fHMuMotherGenPhi[i] = 0x0;
+    fHMuMotherRecPhi[i] = 0x0;
+    fHMuMotherGenEta[i] = 0x0;
+    fHMuMotherRecEta[i] = 0x0;
+    fHMuDCA[i] = 0x0;
+  }
+}
+
+//________________________________________________________________________
+AliMuonEffMC::~AliMuonEffMC()
+{
+  //Destructor
+}
+
+//________________________________________________________________________
+void AliMuonEffMC::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once (per slave on PROOF!)
+
+  fOutputList = new TList();
+  fOutputList->SetOwner(1);
+  fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
+  fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
+  fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
+  fHEventStat->GetXaxis()->SetBinLabel(3,"File");
+  fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
+  fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
+  fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
+  fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
+  fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
+  fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
+  fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
+  fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
+  fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
+  fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
+  fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
+  fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
+  fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
+  fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
+  fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
+  fOutputList->Add(fHEventStat);
+
+  fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
+  fOutputList->Add(fHEvt);
+
+
+  Int_t iTrackBin[5];
+  Double_t* trackBins[5];
+  const char* trackAxisTitle[5];
+
+  // eta
+  Double_t etaBins[fNEtaBins+1];
+  for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
+  iTrackBin[0] = fNEtaBins;
+  trackBins[0] = etaBins;
+  trackAxisTitle[0] = "#eta";
+
+  // p_T
+  Double_t pTBins[fNpTBins+1];
+  for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
+  iTrackBin[1] = fNpTBins;
+  trackBins[1] = pTBins;
+  trackAxisTitle[1] = "p_{T} (GeV/c)";
+
+  // centrality
+  Double_t CentBins[fNCentBins+1];
+  for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
+  iTrackBin[2] = fNCentBins;
+  trackBins[2] = CentBins;
+  trackAxisTitle[2] = "Cent";
+
+  // Z-vertex
+  Double_t ZvtxBins[fNZvtxBins+1];
+  for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
+  iTrackBin[3] = fNZvtxBins;
+  trackBins[3] = ZvtxBins;
+  trackAxisTitle[3] = "Zvtx";
+
+  // phi
+  Double_t phiBins[fNPhiBins+1];
+  for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
+  iTrackBin[4] = fNPhiBins;
+  trackBins[4] = phiBins;
+  trackAxisTitle[4] = "#phi";
+
+  // THn for tracking efficiency
+  fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
+  for (Int_t i=0; i<5; i++)
+  {
+    fHMuonParGen->SetBinEdges(i, trackBins[i]);
+    fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
+  }
+  fHMuonParGen->Sumw2();
+  fOutputList->Add(fHMuonParGen);
+
+  fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
+  fHMuonDetGen->Sumw2();
+  fOutputList->Add(fHMuonDetGen);
+
+  fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
+  fHMuonDetRec->Sumw2();
+  fOutputList->Add(fHMuonDetRec);
+
+  fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
+  fHEtcDetRec->Sumw2();
+  fOutputList->Add(fHEtcDetRec);
+
+  if(fMDProcess)
+  {
+    const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
+    
+    for(Int_t i=0; i<4; i++)
+    {
+      fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
+      fOutputList->Add(fHMuMotherGenPt[i]);
+      
+      fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
+      fOutputList->Add(fHMuMotherRecPt[i]);
+      
+      fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
+      fOutputList->Add(fHMuMotherGenPhi[i]);
+      
+      fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
+      fOutputList->Add(fHMuMotherRecPhi[i]);
+      
+      fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
+      fOutputList->Add(fHMuMotherGenEta[i]);
+      
+      fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
+      fOutputList->Add(fHMuMotherRecEta[i]);
+      
+      fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
+      fOutputList->Add(fHMuDCA[i]);
+    }
+  }
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliMuonEffMC::UserExec(Option_t *) 
+{
+  // Main loop, Called for each event
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!fESD) 
+  {
+    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+  }
+
+  if (!fESD && !fAOD) 
+  {
+    AliError("Neither fESD nor fAOD available");
+    return;
+  }
+       
+  fHEventStat->Fill(0.5);
+  if(fIsMc) 
+  {
+    fMC = MCEvent();
+    if (!fMC) {
+      printf("ERROR: fMC not available\n");
+      return;
+    }
+  }
+
+  const AliESDVertex* vertex = fESD->GetPrimaryVertex();
+  fZVertex = vertex->GetZ();
+  if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile("fCentralityEstimator");
+
+  // Fill Event histogram
+  fHEvt->Fill(fZVertex, fCentrality);
+
+   // Centrality, vertex, other event variables...
+  if (!VertexOk(fESD)) { 
+    AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));  
+    return; 
+  }
+
+  if (fCentrality > 100. || fCentrality < -1.5) { 
+    AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); 
+    return; 
+  }
+  if(fCentrality < 0) fCentrality = 1.0;
+  fHEventStat->Fill(1.5);
+  ULong64_t trigword=fESD->GetTriggerMask();
+  if (trigword & 0x01) fHEventStat->Fill(17.5);
+  if (trigword & 0x02) fHEventStat->Fill(3.5);
+  if (trigword & 0x04) fHEventStat->Fill(4.5);
+  if (trigword & 0x08) fHEventStat->Fill(5.5);      
+  if (trigword & 0x010) fHEventStat->Fill(6.5);
+  if (trigword & 0x020) fHEventStat->Fill(7.5);
+  if (trigword & 0x040) fHEventStat->Fill(8.5);
+  if (trigword & 0x080) fHEventStat->Fill(9.5);
+  if (trigword & 0x100) fHEventStat->Fill(10.5);
+  if (trigword & 0x200) fHEventStat->Fill(11.5);
+  if (trigword & 0x400) fHEventStat->Fill(12.5);
+  if (trigword & 0x800) fHEventStat->Fill(13.5);
+  if (trigword & 0x1000) fHEventStat->Fill(14.5);
+  if (trigword & 0x2000) fHEventStat->Fill(15.5);
+  if (trigword & 0x4000) fHEventStat->Fill(16.5);
+
+  if(fIsMc)
+  {
+    for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++) 
+    { 
+      // loop on muon tracks
+      AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
+      if(muonTrack)
+      {
+       if(!IsGoodMUONtrack(*muonTrack)) continue;
+       Double_t fillArrayDetRec[5] = { muonTrack->Eta(), muonTrack->Pt(), fCentrality, fZVertex, muonTrack->Phi() };
+       fHMuonDetRec->Fill(fillArrayDetRec); 
+       
+       Int_t label = TMath::Abs(muonTrack->GetLabel());
+       AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
+       if(TMath::Abs(McParticle->PdgCode()) != 13) 
+       {
+         fHEtcDetRec->Fill(fillArrayDetRec); 
+         continue;
+       }
+       Double_t fillArrayDetGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
+       fHMuonDetGen->Fill(fillArrayDetGen);
+       
+       if(fMDProcess)
+       {
+         if(McParticle->GetMother() > 0)
+         {
+           AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
+           Int_t motherlabel = TMath::Abs(MotherParticle->PdgCode());
+           
+           if(motherlabel==411 || motherlabel==413 || motherlabel==421 || motherlabel==423 || motherlabel==431 || motherlabel==433 || motherlabel==10413 || motherlabel==10411 || motherlabel==10423 || motherlabel==10421 || motherlabel==10433 || motherlabel==10431 || motherlabel==20413 || motherlabel==415 || motherlabel==20423 || motherlabel==425 || motherlabel==20433 || motherlabel==435)
+           {  
+             fHMuMotherGenPt[2]->Fill(McParticle->Pt(), MotherParticle->Pt());
+             fHMuMotherRecPt[2]->Fill(muonTrack->Pt(), MotherParticle->Pt());
+             fHMuMotherGenPhi[2]->Fill(McParticle->Phi(), MotherParticle->Phi());
+             fHMuMotherRecPhi[2]->Fill(muonTrack->Phi(), MotherParticle->Phi());
+             fHMuMotherGenEta[2]->Fill(McParticle->Eta(), MotherParticle->Eta());
+             fHMuMotherRecEta[2]->Fill(muonTrack->Eta(), MotherParticle->Eta());
+             fHMuDCA[2]->Fill(muonTrack->GetDCA());
+           }
+           else if(motherlabel==211) 
+           {
+             fHMuMotherGenPt[0]->Fill(McParticle->Pt(), MotherParticle->Pt());
+             fHMuMotherRecPt[0]->Fill(muonTrack->Pt(), MotherParticle->Pt());
+             fHMuMotherGenPhi[0]->Fill(McParticle->Phi(), MotherParticle->Phi());
+             fHMuMotherRecPhi[0]->Fill(muonTrack->Phi(), MotherParticle->Phi());
+             fHMuMotherGenEta[0]->Fill(McParticle->Eta(), MotherParticle->Eta());
+             fHMuMotherRecEta[0]->Fill(muonTrack->Eta(), MotherParticle->Eta());
+             fHMuDCA[0]->Fill(muonTrack->GetDCA());
+           }
+           else if(motherlabel==321)
+           {
+             fHMuMotherGenPt[1]->Fill(McParticle->Pt(), MotherParticle->Pt());
+             fHMuMotherRecPt[1]->Fill(muonTrack->Pt(), MotherParticle->Pt());
+             fHMuMotherGenPhi[1]->Fill(McParticle->Phi(), MotherParticle->Phi());
+             fHMuMotherRecPhi[1]->Fill(muonTrack->Phi(), MotherParticle->Phi());
+             fHMuMotherGenEta[1]->Fill(McParticle->Eta(), MotherParticle->Eta());
+             fHMuMotherRecEta[1]->Fill(muonTrack->Eta(), MotherParticle->Eta());
+             fHMuDCA[1]->Fill(muonTrack->GetDCA());
+           }   
+           else 
+           {
+             fHMuMotherGenPt[3]->Fill(McParticle->Pt(), MotherParticle->Pt());
+             fHMuMotherRecPt[3]->Fill(muonTrack->Pt(), MotherParticle->Pt());
+             fHMuMotherGenPhi[3]->Fill(McParticle->Phi(), MotherParticle->Phi());
+             fHMuMotherRecPhi[3]->Fill(muonTrack->Phi(), MotherParticle->Phi());
+             fHMuMotherGenEta[3]->Fill(McParticle->Eta(), MotherParticle->Eta());
+             fHMuMotherRecEta[3]->Fill(muonTrack->Eta(), MotherParticle->Eta());
+             fHMuDCA[3]->Fill(muonTrack->GetDCA());
+           }
+         }    
+       } // (mother hadron) : (daughter muon) QA 
+      } 
+    }
+    
+    for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
+    {
+      AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
+      if(TMath::Abs(McParticle->PdgCode())==13) 
+      {
+       Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
+       fHMuonParGen->Fill(fillArrayParGen);
+      }
+    }  
+  }
+  PostData(1, fOutputList);
+  return;
+}
+
+//________________________________________________________________________
+void AliMuonEffMC::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    AliError("Output list not available");
+    return;
+  }
+}
+
+//________________________________________________________________________
+Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
+{
+  // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
+  
+  Int_t nContributors  = 0;
+  Double_t zVertex     = 999;
+  TString name("");
+  
+  if (obj->InheritsFrom("AliESDEvent")) {
+    AliESDEvent* esdevt = (AliESDEvent*) obj;
+    const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
+    if (!vtx)
+      return 0;
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  else if (obj->InheritsFrom("AliAODEvent")) { 
+    AliAODEvent* aodevt = (AliAODEvent*) obj;
+    if (aodevt->GetNumberOfVertices() < 1)
+      return 0;
+    const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  
+  // Reject if TPC-only vertex
+  if (name.CompareTo("TPCVertex")==0)
+    return kFALSE;
+  
+  // Check # contributors and range...
+  if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
+    return kFALSE;
+  }
+  
+  return kTRUE;
+}
+
+
+//________________________________________________________________________
+Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
+{
+  // Applying track cuts for MUON tracks
+  if(!track.ContainTrackerData()) return kFALSE;
+  if(!track.ContainTriggerData()) return kFALSE;
+
+  Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
+  Double_t eta = track.Eta();
+
+  // Theta cut at absorber end
+  if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
+  // Eta cut
+  if(eta <= -4. || eta >= -2.5) return kFALSE;
+  if(track.GetMatchTrigger() <= 0) return kFALSE;
+
+  return kTRUE;
+}
diff --git a/PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.h b/PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.h
new file mode 100644 (file)
index 0000000..ee52235
--- /dev/null
@@ -0,0 +1,88 @@
+// MUON tracking efficiency + (mother hadron) : (daughter muon) kinematic relation 
+// Author : Saehanseul Oh
+
+
+#ifndef AliMuonEffMC_h
+#define AliMuonEffMC_h
+
+class TH1F;
+class TH1D;
+class TH2F;
+class TH3F;
+class THn;
+class TList;
+class TObjArray;
+class TObject;
+
+class AliAODEvent;
+class AliESDEvent;
+class AliESDtrackCuts;
+class AliEvtPoolManager;
+class AliMCEvent;
+class AliESDMuonTrack;
+#include "AliAnalysisTaskSE.h"
+
+
+class AliMuonEffMC : public AliAnalysisTaskSE {
+ public:
+  AliMuonEffMC();
+  AliMuonEffMC(const char *name);
+  virtual ~AliMuonEffMC();
+  void         UserCreateOutputObjects();
+  void         UserExec(Option_t *option);
+  void         Terminate(Option_t *);
+
+  void         SetMcAna(Bool_t IsMc)               { fIsMc = IsMc;               }
+  void         SetMDProcess(Bool_t MDProcess)      { fMDProcess = MDProcess;     }
+  void         SetCentEstimator(TString Cent)      { fCentralityEstimator = Cent;}
+  void         SetNEtaBins(Int_t NEtaBins)         { fNEtaBins = NEtaBins;       }
+  void         SetNpTBins(Int_t NpTBins)           { fNpTBins = NpTBins;         }
+  void         SetNCentBins(Int_t NCentBins)       { fNCentBins = NCentBins;     }
+  void         SetNZvtxBins(Int_t NZvtxBins)       { fNZvtxBins = NZvtxBins;     }
+  void         SetNPhiBins(Int_t NPhiBins)         { fNPhiBins = NPhiBins;       }
+
+ protected:
+  Bool_t       VertexOk(TObject* obj) const;
+  Bool_t       IsGoodMUONtrack(AliESDMuonTrack &track);
+
+ private:
+  AliESDEvent *fESD;               //! ESD object
+  AliAODEvent *fAOD;               //! AOD object 
+  AliMCEvent  *fMC;                //! MC object
+  Double_t     fCentrality;        //! Of current event
+  Double_t     fZVertex;           //! Of current event
+  TList       *fOutputList;        //! Output list
+  TH1D        *fHEventStat;        //! statistics histo
+  TH2F        *fHEvt;              //! Cent, vtx
+
+  Bool_t       fIsMc;              //! 
+  Bool_t       fMDProcess;         //! (mother hadron) : (daughter muon) QA 
+
+  TString      fCentralityEstimator;//! 
+  Int_t        fNEtaBins;          //! number of eta bins
+  Int_t        fNpTBins;           //! number of p_T bins
+  Int_t        fNCentBins;         //! number of centrality bins
+  Int_t        fNZvtxBins;         //! number of Z-vertex bins
+  Int_t        fNPhiBins;          //! number of phi bins
+
+  THn         *fHMuonParGen;       //! truth muon track eta, p_T, Centrality, Z-vertex, phi 
+  THn         *fHMuonDetGen;       //! detector level muon track generated eta, p_T, Centrality, Z-vertex, phi 
+  THn         *fHMuonDetRec;       //! reconstructed muon track eta, p_T, Centrality, Z-vertex, phi 
+  THn         *fHEtcDetRec;        //! particle reconstructed at MUON detector, but not muon
+
+  TH2F        *fHMuMotherGenPt[4]; //! particle-level muon p_T vs. mother p_T
+  TH2F        *fHMuMotherRecPt[4]; //! detector-level muon p_T vs. mother p_T
+  TH2F        *fHMuMotherGenPhi[4];//! particle-level muon phi vs. mother phi
+  TH2F        *fHMuMotherRecPhi[4];//! detector-level muon phi vs. mother phi
+  TH2F        *fHMuMotherGenEta[4];//! particle-level muon eta vs. mother eta
+  TH2F        *fHMuMotherRecEta[4];//! detector-level muon eta vs. mother eta
+  TH1F        *fHMuDCA[4];         //! muon DCA
+
+  AliMuonEffMC(const AliMuonEffMC&);            // not implemented
+  AliMuonEffMC &operator=(const AliMuonEffMC&); // not implemented
+
+  ClassDef(AliMuonEffMC, 1);
+};
+
+#endif
index c977d7f1971d4a365512d417ee4805f806819b97..aa72fbd2f1a2771c2c8d5ad55bad437691ce2ddf 100644 (file)
@@ -17,6 +17,7 @@
 #pragma link C++ class AliEvtPool+;
 #pragma link C++ class AliMiniTrack+;
 #pragma link C++ class AliDhcTask+;
+#pragma link C++ class AliMuonEffMC+;
 #pragma link C++ class AliLeadingV0Correlation+;
 #pragma link C++ class AliLeadingBasicParticle+;
 #pragma link C++ class AliAnalysisTaskLongRangeCorrelations+;