Add task for single muon flow estimation with the EP method (Diego)
authorpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 22 Jul 2012 13:53:23 +0000 (13:53 +0000)
committerpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 22 Jul 2012 13:53:23 +0000 (13:53 +0000)
PWG/CMakelibPWGmuon.pkg
PWG/PWGmuonLinkDef.h
PWG/muon/AddTaskFlowSingleMu.C [new file with mode: 0644]
PWG/muon/AliAnalysisTaskFlowSingleMu.cxx [new file with mode: 0644]
PWG/muon/AliAnalysisTaskFlowSingleMu.h [new file with mode: 0644]

index c9f50b7..73e2c79 100644 (file)
@@ -61,6 +61,7 @@ set ( SRCS
     muon/AliMuonPairCuts.cxx 
     muon/AliMergeableCollection.cxx 
     muon/AliVAnalysisMuon.cxx 
+    muon/AliAnalysisTaskFlowSingleMu.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index 628c3d4..66e18c1 100644 (file)
@@ -42,5 +42,6 @@
 #pragma link C++ class AliMergeableCollection+;
 #pragma link C++ class AliMergeableCollectionIterator+;
 #pragma link C++ class AliVAnalysisMuon+;
+#pragma link C++ class AliAnalysisTaskFlowSingleMu+;
 #endif
 
diff --git a/PWG/muon/AddTaskFlowSingleMu.C b/PWG/muon/AddTaskFlowSingleMu.C
new file mode 100644 (file)
index 0000000..ab95a5e
--- /dev/null
@@ -0,0 +1,51 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TString.h"
+#include "TObjArray.h"
+
+#include "AliLog.h"
+#include "AliVEventHandler.h"
+
+#include "AliAnalysisManager.h"
+#include "AliAnalysisDataContainer.h"
+
+#include "AliMuonTrackCuts.h"
+#include "AliAnalysisTaskFlowSingleMu.h"
+#endif
+
+AliAnalysisTaskFlowSingleMu* AddTaskFlowSingleMu(Bool_t isMC = kFALSE)
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddtaskFlowSingleMu", "No analysis manager to connect to.");
+    return NULL;
+  }
+
+  TString type = mgr->GetInputEventHandler()->GetDataType();
+  if (!type.Contains("ESD") && !type.Contains("AOD")) {
+    ::Error("AddtaskFlowSingleMu", "FlowSingleMu task needs the manager to have an ESD or AOD input handler.");
+    return NULL;
+  }
+
+  // Create container
+  TString currName = "";
+  TString outputfile = mgr->GetCommonFileName();
+  if ( ! outputfile.IsNull() ) outputfile += ":PWGHF_FlowSingleMu";
+  else outputfile = "FlowSingleMuAnalysis.root";
+
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("FlowSingleMuOut",TObjArray::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+
+  // Create cuts
+  AliMuonTrackCuts* muonTrackCuts = new AliMuonTrackCuts("TestStandardMuonTrackCuts", "TestStandardMuonTrackCuts");
+  muonTrackCuts->SetIsMC(isMC);
+
+  // Create task
+  AliAnalysisTaskFlowSingleMu *flowSingleMuTask = new AliAnalysisTaskFlowSingleMu("FlowSingleMuTask", *muonTrackCuts);
+  if ( isMC ) flowSingleMuTask->SetTrigClassPatterns("ANY");
+  mgr->AddTask(flowSingleMuTask);
+
+  // Connect containers
+  mgr->ConnectInput  (flowSingleMuTask,  0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput (flowSingleMuTask,  1, coutput1);
+
+  return flowSingleMuTask;
+}
diff --git a/PWG/muon/AliAnalysisTaskFlowSingleMu.cxx b/PWG/muon/AliAnalysisTaskFlowSingleMu.cxx
new file mode 100644 (file)
index 0000000..7a9bced
--- /dev/null
@@ -0,0 +1,796 @@
+/**************************************************************************
+ * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliAnalysisTaskFlowSingleMu.cxx 55545 2012-04-04 07:16:39Z pcrochet $ */
+
+//-----------------------------------------------------------------------------
+/// \class AliAnalysisTaskFlowSingleMu
+/// Analysis task for flow of single muons in the spectrometer.
+/// The output is a list of histograms and CF containers.
+/// The macro class can run on AODs or ESDs.
+/// If Monte Carlo information is present, some basics checks are performed.
+///
+/// \author Diego Stocco
+//-----------------------------------------------------------------------------
+
+#define AliAnalysisTaskFlowSingleMu_cxx
+
+#include "AliAnalysisTaskFlowSingleMu.h"
+
+// ROOT includes
+#include "TROOT.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TGraphErrors.h"
+#include "TAxis.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#include "TObjString.h"
+#include "TObjArray.h"
+#include "TF1.h"
+#include "TStyle.h"
+//#include "TMCProcess.h"
+#include "TFitResultPtr.h"
+#include "TFile.h"
+#include "TArrayI.h"
+#include "TArrayD.h"
+#include "TSystem.h"
+#include "TGraphErrors.h"
+#include "TLatex.h"
+
+// STEER includes
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliVHeader.h"
+#include "AliAODMCHeader.h"
+#include "AliStack.h"
+#include "AliEventplane.h"
+
+// ANALYSIS includes
+#include "AliAnalysisManager.h"
+
+// CORRFW includes
+#include "AliCFContainer.h"
+#include "AliCFGridSparse.h"
+#include "AliCFEffGrid.h"
+
+// PWG includes
+#include "AliVAnalysisMuon.h"
+#include "AliMergeableCollection.h"
+#include "AliCounterCollection.h"
+#include "AliMuonTrackCuts.h"
+
+
+/// \cond CLASSIMP
+ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
+/// \endcond
+
+
+//________________________________________________________________________
+AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
+  AliVAnalysisMuon(),
+  fEPKeys(0x0),
+  fRandom(0x0)
+{
+  /// Default ctor.
+}
+
+//________________________________________________________________________
+AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
+  AliVAnalysisMuon(name, cuts),
+  fEPKeys(0x0),
+  fRandom(0x0)
+{
+  //
+  /// Constructor.
+  //
+  TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
+  fEPKeys = epMethods.Tokenize(" ");
+}
+
+
+//________________________________________________________________________
+AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
+{
+  //
+  /// Destructor
+  //
+  delete fEPKeys;
+  delete fRandom;
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
+{
+  //
+  /// Create outputs
+  //
+
+  if ( ! fRandom ) fRandom = new TRandom3();
+  // Set seed here, to avoid having the same seed for all jobs in grid
+  fRandom->SetSeed(0);
+
+  Int_t nPtBins = 160;
+  Double_t ptMin = 0., ptMax = 80.;
+  TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
+  
+  Int_t nEtaBins = 25;
+  Double_t etaMin = -4.5, etaMax = -2.;
+  TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
+  
+  Int_t nPhiBins = 36;
+  Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
+  TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
+
+  Int_t nDeltaPhiBins = 18;
+  Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
+  TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
+    
+  Int_t nChargeBins = 2;
+  Double_t chargeMin = -2, chargeMax = 2.;
+  TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
+  
+  Int_t nMotherTypeBins = kNtrackSources;
+  Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
+  TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
+
+  Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
+  Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
+  Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
+  TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
+  TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
+
+  Int_t iobj=0;
+  TH1* histo = 0x0;
+  TString currTitle = "";
+  for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
+    TString currEP = fEPKeys->At(iep)->GetName();
+    AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
+  
+    TString currName = "";
+    for ( Int_t idim = 0; idim<kNvars; idim++){
+      currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
+      currName.ReplaceAll("()","");
+
+      cfContainer->SetVarTitle(idim, currName.Data());
+      cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
+    }
+    
+    TString stepTitle[kNsteps] = {"reconstructed", "generated"};
+
+    TAxis* currAxis = 0x0;
+    for (Int_t istep=0; istep<kNsteps; istep++){
+      cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
+      AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
+
+      currAxis = gridSparse->GetAxis(kHvarMotherType);
+      for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
+        currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
+      }
+    }
+
+    AddObjectToCollection(cfContainer, iobj++);
+
+    Bool_t isTPC = currEP.Contains("TPC");
+    Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
+    Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
+    
+    histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
+    histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
+    AddObjectToCollection(histo, iobj++);
+
+    for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
+      TString auxEP = fEPKeys->At(jep)->GetName();
+      currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
+      histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
+      histo->GetXaxis()->SetTitle(currTitle.Data());
+      AddObjectToCollection(histo, iobj++);
+    }
+
+    currTitle = Form("cos(#psi_{1} - #psi_{2})");
+    histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
+    histo->GetXaxis()->SetTitle(currTitle.Data());
+
+    AddObjectToCollection(histo, iobj++);
+  } // loop on EPs
+
+  fMuonTrackCuts->Print("mask");
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
+{
+  //
+  /// Fill output objects
+  //
+
+  //
+  // Base event cuts
+  //
+  if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
+  //Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes
+  //if ( isPileupFromSPD ) return;
+  if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) return;
+
+  TArrayD psiPlane(fEPKeys->GetEntriesFast());
+  psiPlane.Reset(-999.);
+  const Double_t kLowestVal = -10.;
+  Int_t harmonic = 2;
+
+  //
+  // Fill EP info
+  //
+  AliEventplane* eventPlane = InputEvent()->GetEventplane();
+  psiPlane[0] = eventPlane->GetEventplane("V0A",InputEvent(),harmonic);
+  psiPlane[1] = eventPlane->GetEventplane("V0C",InputEvent(),harmonic);
+  psiPlane[2] = eventPlane->GetEventplane("Q");
+  if ( psiPlane[2] < 0.) psiPlane[2] = -999.;
+  Double_t qx = 0., qy = 0.;
+  psiPlane[3] = eventPlane->CalculateVZEROEventPlane(InputEvent(),0,harmonic,qx,qy); // V0C ring 1
+  psiPlane[4] = eventPlane->CalculateVZEROEventPlane(InputEvent(),3,harmonic,qx,qy); // V0C ring 3
+  //psiPlane[3] = fRandom->Uniform(-TMath::Pi()/2., TMath::Pi()/2.);
+
+  Bool_t noValidEP = kTRUE;
+  TArrayD psiPlaneCalc(fEPKeys->GetEntriesFast());
+  for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
+    psiPlaneCalc[iep] = psiPlane[iep];
+    if ( psiPlane[iep] < kLowestVal ) continue;
+    if ( psiPlane[iep] < 0. ) psiPlaneCalc[iep] += TMath::Pi();
+    noValidEP = kFALSE;
+  }
+  if ( noValidEP ) return;
+
+  for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
+    if ( psiPlane[iep] < kLowestVal ) continue; 
+    for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
+      TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
+      TString currEP = fEPKeys->At(iep)->GetName();
+      ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hEventPlane%s",currEP.Data())))->Fill(psiPlane[iep]);
+      for ( Int_t jep=iep+1; jep<fEPKeys->GetEntriesFast(); jep++ ) {
+        if ( psiPlane[jep] < kLowestVal ) continue; 
+        ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub3%s%s",currEP.Data(), fEPKeys->At(jep)->GetName())))->Fill(TMath::Cos(2.*(psiPlaneCalc[iep]-psiPlaneCalc[jep])));
+      } // loop on auxialiary EP
+      if ( iep == 2 ) {
+        TVector2* qsub1 = eventPlane->GetQsub1();
+        TVector2* qsub2 = eventPlane->GetQsub2();
+        if ( qsub1 && qsub2 )
+          ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.));
+      }
+    } // loop on trigger
+  } // loop on EP
+
+  //
+  // Fill track info
+  //
+  Double_t containerInput[kNvars];
+  AliVParticle* track = 0x0;
+  for ( Int_t istep = 0; istep<2; ++istep ) {
+    Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
+    for (Int_t itrack = 0; itrack < nTracks; itrack++) {
+      track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
+      
+      Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
+      if ( ! isSelected ) continue;
+      
+      Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
+
+      containerInput[kHvarPt]         = track->Pt();
+      containerInput[kHvarEta]        = track->Eta();
+      containerInput[kHvarPhi]        = track->Phi();
+      containerInput[kHvarCharge]     = track->Charge()/3.;
+      containerInput[kHvarMotherType] = (Double_t)trackSrc;
+
+      for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
+        if ( psiPlane[iep] < kLowestVal ) continue; 
+        TString currEP = fEPKeys->At(iep)->GetName();
+        Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
+
+        // The difference between phi - psiPlane must be between (0, pi)
+        // Since the reaction plain is symmetric in pi,
+        // we can do in such a way that:
+        // -pi < delta_phi + N*pi < pi
+        // We choose N accrodingly
+        Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
+        deltaPhi += nPi * TMath::Pi();
+        containerInput[kHvarDeltaPhi] = deltaPhi;
+
+        for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
+          TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
+          if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
+          ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
+        } // loop on selected trigger classes
+      } // loop on event plane
+    } // loop on tracks
+  } // loop on container steps
+}
+
+
+//________________________________________________________________________
+TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
+{
+  //
+  /// Get the range from string
+  //
+  TArrayD centralityRange(2);
+  centralityRange.Reset(-999.);
+  TObjArray* array = sRange.Tokenize("_");
+  if ( array->GetEntries() == 2 ) {
+    for ( Int_t ientry=0; ientry<2; ientry++ ) {
+      centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
+    }
+  }
+  delete array;
+  return centralityRange;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
+  //
+  /// Draw some histograms at the end.
+  //
+
+  AliVAnalysisMuon::Terminate("");
+
+  if ( ! fMergeableCollection ) return;
+  
+  TString physSel = fTerminateOptions->At(0)->GetName();
+  TString trigClassName = fTerminateOptions->At(1)->GetName();
+  TString centralityRange = fTerminateOptions->At(2)->GetName();
+  TString furtherOpt = fTerminateOptions->At(3)->GetName();
+
+  TObjArray* centralityClasses = centralityRange.Tokenize(",");
+  TArrayD centralityArray(centralityClasses->GetEntries()+1);
+  Int_t ib = 0;
+  for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+    TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
+    if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
+    centralityArray.AddAt(range[1],ib++);
+  }
+  
+  TString selectEP = "";
+  TObjArray resoEParray;
+  resoEParray.SetOwner();
+  
+  TString outFilename = "";
+  TObjArray* optArr = furtherOpt.Tokenize(" ");
+  TString currName = "";
+  TString resoFilename = "";
+  for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
+    currName = optArr->At(iopt)->GetName();
+    if ( currName.Contains(".root") ) outFilename = currName;
+    else if ( currName.Contains(".txt") ) resoFilename = currName;
+    else {
+      for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
+        if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
+      }
+    }
+  }
+  delete optArr;
+
+  if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
+
+  furtherOpt.ToUpper();
+
+  Bool_t v2only = furtherOpt.Contains("V2ONLY");
+
+  TAxis* customBins = 0x0;
+  
+  //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
+  Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
+  //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
+  //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
+  //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
+  Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
+  customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
+
+  TCanvas *can = 0x0, *showCan = 0x0;
+  Int_t xshift = 100;
+  Int_t yshift = 20;
+  Int_t igroup1 = 0;
+  Int_t igroup2 = 0;
+
+  //
+  /// Read plane resolution file (if present)
+  //
+  TObjArray resoCentrality;
+  resoCentrality.SetOwner();
+  TArrayD resolution[3];
+  for ( Int_t ires=0; ires<3; ires++ ) {
+    resolution[ires].Set(50);
+  }
+
+  Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
+  Bool_t internalReso = kTRUE;
+  if ( externalReso ) {
+    // If plane resolution file exists, fill resolution
+    TString currLine = "";
+    ifstream inFile(resoFilename.Data());
+    if (inFile.is_open()) {
+      ib = 0;
+      while (! inFile.eof() ) {
+        currLine.ReadLine(inFile,kFALSE);
+        if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
+        TObjArray* array = currLine.Tokenize(" ");
+        resoCentrality.Add(new TObjString(array->At(0)->GetName()));
+        for ( Int_t ires=0; ires<3; ires++ ) {
+          resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
+        }
+        delete array;
+        ib++;
+      } // ! EOF
+    } // file is open
+    for ( Int_t ires=0; ires<3; ires++ ) {
+      resolution[ires].Set(ib);
+    }
+  } // file exists
+  else {
+    // If no plane resolution file exist,
+    // use the centralities in terminateOption
+    // and set the resolution to 1 and errors to 0.
+
+    for ( Int_t ires=0; ires<3; ires++ ) {
+      resolution[ires].Set(centralityClasses->GetEntries());
+      Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
+      resolution[ires].Reset(resetVal);
+    }
+    TArrayD currCos(3);
+    for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+      resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
+      Int_t icos = 0;
+      Double_t sumRelErrSquare = 0.;
+      for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
+        for ( Int_t jsub=0; jsub<resoEParray.GetEntries(); jsub++ ) {
+          TH1* histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName())));
+          if ( ! histo || histo->Integral() == 0 ) continue;
+          currCos[icos] = histo->GetMean();
+          if ( currCos[icos] <= 0. ) continue;
+          Double_t relErr = histo->GetMeanError() / currCos[icos];
+          sumRelErrSquare += relErr * relErr;
+          icos++;
+        } // loop on sub events
+      } // loop on sub events
+      if ( icos != 3  ) {
+        printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
+        internalReso = kFALSE;
+        continue;
+      }
+      resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
+      resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
+      printf("Resolution %s  %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
+    }
+  }
+
+  Bool_t hasResolution = ( externalReso || internalReso );
+
+
+  TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
+  ib = 0;
+  for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
+    TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
+    if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
+    resoCentrArray.AddAt(range[1], ib++);
+  }
+
+  if ( hasResolution ) {
+    currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
+    can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+    TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries());
+    for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
+
+      resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
+      resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
+    }
+    resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
+    resoGraph->GetXaxis()->SetTitle("Centrality");
+    resoGraph->GetYaxis()->SetTitle("Resolution");
+    resoGraph->Draw("apz");
+  }
+
+  const Int_t kNresoCentr = resoCentrality.GetEntries();
+
+  //
+  // Create gridSparse array
+  //
+  AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
+  TString stepName[2];
+  for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
+    for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+      gridSparseArray[istep*kNresoCentr+icent] = 0x0;
+    }
+
+    AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
+    if ( ! cfContainer ) continue;
+    for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+      stepName[istep] = cfContainer->GetStepTitle(istep);
+      gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
+    } // loop on step
+  } // loop on centrality
+
+  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
+  
+  Bool_t isMC = furtherOpt.Contains("MC");
+  Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
+  Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
+  if ( ! isMC ) srcColors[kUnidentified] = 1;
+
+  TList outList;
+
+  TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
+  Int_t drawOrder[3] = {0, 2, 1};
+  TString drawOpt[3] = {"apz","pz","a2"};
+  Int_t nTypes = ( hasResolution ) ? 3 : 1;
+
+  TString histoName = "", histoPattern = "";
+  ///////////
+  // Flow  //
+  ///////////
+  gStyle->SetOptFit(1111);
+  TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))";
+  //
+  TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
+  if  ( v2only ) func->SetParNames("scale", "v2");
+  else func->SetParNames("scale","v1", "v2", "v3");
+  Int_t v2par = ( v2only ) ? 1 : 2;
+
+  for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+    for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
+      TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
+      TGraphErrors* flowVsCentrality[3];
+      for ( Int_t itype=0; itype<nTypes; itype++ ) {
+        histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
+        flowVsCentrality[itype] = new TGraphErrors();
+        flowVsCentrality[itype]->SetName(histoName.Data());
+        flowVsCentrality[itype]->SetTitle(histoName.Data());
+      }
+      for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+        TH2* hDeltaPhi2D = 0x0;
+        TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
+        TArrayD weightedReso(3);
+        weightedReso.Reset(0.);
+        Double_t nMuons = 0.;
+        TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
+        for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
+          if ( resoCentrArray[rcent] < centralityArray[icent] ||
+               resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
+          AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
+          if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
+
+          SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
+          SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
+          TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
+
+          Double_t currYield = auxHisto->Integral();
+          if ( currYield == 0. ) continue;
+          nMuons += currYield;
+          debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
+          for ( Int_t ires=0; ires<3; ires++ ) {
+            weightedReso[ires] += resolution[ires][rcent] * currYield;
+            debugString += Form("  %g", resolution[ires][rcent]);
+          }
+          histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data());
+          if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast<TH2*>(auxHisto->Clone(histoName.Data()));
+          else hDeltaPhi2D->Add(auxHisto);
+          delete auxHisto;
+        }
+
+        if ( ! hDeltaPhi2D ) continue;
+
+        debugString += Form("\nWeighted reso ");
+        for ( Int_t ires=0; ires<3; ires++ ) {
+          weightedReso[ires] = weightedReso[ires]/nMuons;
+          debugString += Form("  %g", weightedReso[ires]);
+        }
+        AliDebug(1,debugString.Data());
+
+        TAxis* ptAxis = hDeltaPhi2D->GetYaxis();
+
+        Int_t ipad = 0;
+        TGraphErrors* flowVsPt[3];
+        for ( Int_t itype=0; itype<nTypes; itype++ ) {
+          histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
+          flowVsPt[itype] = new TGraphErrors();
+          flowVsPt[itype]->SetName(histoName.Data());
+          flowVsPt[itype]->SetTitle(histoName.Data());
+        }
+
+        if ( ! customBins ) customBins = static_cast<TAxis*>(ptAxis->Clone());
+
+
+        for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
+          Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
+          Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
+          Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
+          Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
+          Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
+          Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
+          TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin);
+          projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
+          if ( projHisto->Integral() < 50 ) break;
+          if ( ipad%4 == 0 ) {
+            currName = projHisto->GetName();
+            currName.Append("_can");
+            showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+            showCan->Divide(2,2);
+            ipad = 0;
+          }
+          showCan->cd(++ipad);
+          if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
+          else {
+            func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
+            func->FixParameter(1,0.);
+          }
+          projHisto->Fit(func,"R");
+          for ( Int_t itype=0; itype<nTypes; itype++ ) {
+            TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
+            Double_t resoVal = ( itype == 0 ) ? 1. : weightedReso[0];
+            Double_t resoErr = ( itype == 0 ) ? 0. : weightedReso[itype];
+            Double_t rawVal = func->GetParameter(v2par);
+            Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
+            Double_t yVal = rawVal/resoVal;
+            Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
+            Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
+            Int_t ipoint = currGraph->GetN();
+            currGraph->SetPoint(ipoint,xVal,yVal);
+            if ( itype > 0 && ipt != 0 ) {
+              // For the plot vs. pt the resolution error is fully correlated.
+              // Hence we do not use add it bin by bin, but rather write it
+              // on the title
+              resoErr = 0.;
+              TString graphTitle = currGraph->GetTitle();
+              if ( ! graphTitle.Contains("|") ) {
+                graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]);
+                currGraph->SetTitle(graphTitle.Data());
+              }
+            }
+            Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
+            Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
+            Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
+            currGraph->SetPointError(ipoint,xErr,yErr);
+          } // loop on type
+        } // loop on pt bins
+        for ( Int_t itype=0; itype<nTypes; itype++ ) {
+          Int_t currType = drawOrder[itype];
+          if ( itype < 2 ) {
+            currName = flowVsPt[currType]->GetName();
+            currName.Append("Can");
+            can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+            igroup2++;
+            can->cd();
+          }
+          flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
+          flowVsPt[currType]->GetYaxis()->SetTitle("v2");
+          flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
+          flowVsPt[currType]->SetFillStyle(0);
+          flowVsPt[currType]->Draw(drawOpt[currType].Data());
+          outList.Add(flowVsPt[currType]);
+        } // loop on types
+      } // loop on centrality
+      for ( Int_t itype=0; itype<nTypes; itype++ ) {
+        Int_t currType = drawOrder[itype];
+        if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
+        if ( itype < 2 ) {
+          currName = flowVsCentrality[currType]->GetName();
+          currName.Append("Can");
+          can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+          can->cd();
+        }
+        flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
+        flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
+        flowVsCentrality[currType]->SetFillStyle(0);
+        flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
+        flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
+        outList.Add(flowVsCentrality[currType]);
+        igroup2++;
+      } // loop on type
+      igroup1++;
+      igroup2 = 0;
+    } // loop on track sources      
+  } // loop on steps
+
+  delete customBins;
+
+
+  ///////////////////////
+  // Event plane check //
+  ///////////////////////
+  igroup1++;
+  igroup2 = 0;
+  Int_t icolor=0;
+  currName ="eventPlaneDistrib";
+  can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+  can->SetGridy();
+  igroup2++;
+
+  Bool_t addOffsetToCentrality = kFALSE;
+
+  TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
+  leg->SetBorderSize(1);
+  for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+    TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
+    if ( ! histo ) continue;
+    histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
+    if ( histo->Integral() < 50. ) continue;
+    if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
+    histo->Fit("pol0","R0");
+    TF1* pol0 = histo->GetFunction("pol0");
+    histo->Scale(1./pol0->GetParameter(0));
+    Int_t offset = ( addOffsetToCentrality ) ? icolor : 0;
+    TGraphErrors* graph = new TGraphErrors();
+    graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
+    for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
+      Int_t ibin = ipoint+1;
+      graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
+      graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
+    }
+    graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
+    graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
+    TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
+    if ( addOffsetToCentrality ) {
+      graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries());
+      graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
+      legTitle += Form(" ( + %i)", offset);
+    }
+    Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
+    graph->SetLineColor(currColor);
+    graph->SetMarkerColor(currColor);
+    graph->SetMarkerStyle(20+icolor);
+    TString currDrawOpt = "pz";
+    if ( icolor == 0 ) currDrawOpt += "a";
+    graph->Draw(currDrawOpt.Data());
+    legTitle.ReplaceAll("_"," - ");
+    leg->AddEntry(graph,legTitle.Data(),"lp");
+    if ( addOffsetToCentrality ) {
+      TLatex* latex = new TLatex(histo->GetXaxis()->GetXmax()+0.5*histo->GetXaxis()->GetBinWidth(1), 1.+(Double_t)offset, Form("#chi^{2}/NDF = %.3g", pol0->GetChisquare()/(Double_t)pol0->GetNDF()));
+      latex->SetTextColor(currColor);
+      latex->SetTextSize(0.025);
+      latex->Draw("same");
+    }
+    icolor++;
+  }
+  leg->Draw("same");
+
+
+  //////////////////////
+  // Event statistics //
+  //////////////////////
+  printf("\nTotal analyzed events:\n");
+  TString evtSel = Form("trigger:%s", trigClassName.Data());
+  fEventCounters->PrintSum(evtSel.Data());
+  printf("Physics selected analyzed events:\n");
+  evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
+  fEventCounters->PrintSum(evtSel.Data());
+  
+  TString countPhysSel = "any";
+  if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
+  else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
+  countPhysSel.Prepend("selected:");
+  printf("Analyzed events vs. centrality:\n");
+  evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
+  fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
+
+  if ( ! outFilename.IsNull() ) {
+    printf("\nWriting output file %s\n", outFilename.Data());
+    TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
+    outList.Write();
+    file->Close();
+  }
+}
diff --git a/PWG/muon/AliAnalysisTaskFlowSingleMu.h b/PWG/muon/AliAnalysisTaskFlowSingleMu.h
new file mode 100644 (file)
index 0000000..8a2dfaf
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALIANALYSISTASKFLOWSINGLEMU_H
+#define ALIANALYSISTASKFLOWSINGLEMU_H
+
+/* $Id: AliAnalysisTaskFlowSingleMu.h 55545 2012-04-04 07:16:39Z pcrochet $ */ 
+
+//
+// AliAnalysisTaskFlowSingleMu
+// Analysis task for flow of single muons in the spectrometer
+//
+//  Author: Diego Stocco
+//
+
+#include "AliVAnalysisMuon.h"
+#include "TRandom3.h"
+
+class TObjArray;
+class TString;
+class TArrayD;
+class TAxis;
+class AliMuonTrackCuts;
+
+class AliAnalysisTaskFlowSingleMu : public AliVAnalysisMuon {
+ public:
+  AliAnalysisTaskFlowSingleMu();
+  AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts);
+  virtual ~AliAnalysisTaskFlowSingleMu();
+
+  virtual void   Terminate(Option_t *option);
+
+  void MyUserCreateOutputObjects();
+  void ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality);
+
+  /*
+  enum {
+    kEPV0A, ///< EP form V0A
+    kEPTPC, ///< EP form TPC
+    kEPrandom ///< Random EP
+  };
+
+  void SetEPtype ( Int_t epType = kEPV0A ) { fEPtype = epType; }
+  */
+
+ private:
+
+  AliAnalysisTaskFlowSingleMu(const AliAnalysisTaskFlowSingleMu&);
+  AliAnalysisTaskFlowSingleMu& operator=(const AliAnalysisTaskFlowSingleMu&);
+
+  TArrayD GetCentralityRange(TString sRange);
+
+  /*
+  enum {
+    kTrackContainer, ///< CF container for tracks
+    kHistoEP,      ///< Event plane distribution
+    kNobjectTypes    ///< Number of objects
+  };
+  */
+  
+  enum {
+    kStepReconstructed,  ///< Reconstructed tracks
+    kStepGeneratedMC,    ///< Generated tracks (MC)
+    kNsteps              ///< Number of steps
+  };  
+  
+  enum {
+    kHvarPt,         ///< Pt at vertex
+    kHvarEta,        ///< Pseudo-Rapidity
+    kHvarPhi,        ///< Phi
+    kHvarDeltaPhi,   ///< Phi_mu - Psi_plane
+    kHvarCharge,     ///< Particle charge
+    kHvarMotherType, ///< Mother type (MC only)
+    kNvars           ///< THnSparse dimensions
+  };
+
+  TObjArray* fEPKeys; ///< EP keys
+  TRandom3* fRandom; //!< Random number generator
+
+  ClassDef(AliAnalysisTaskFlowSingleMu, 1); // Single muon analysis
+};
+
+#endif