adding task used for cross section and two-body analysis as well as the code of the...
authorfreidt <freidt@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Mar 2013 15:31:38 +0000 (15:31 +0000)
committerfreidt <freidt@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Mar 2013 15:31:38 +0000 (15:31 +0000)
PWGUD/CMakelibPWGUDdiffractive.pkg
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.cxx [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.h [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.cxx [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.h [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.cxx [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.h [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.cxx [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.h [new file with mode: 0644]
PWGUD/DIFFRACTIVE/xsAndTwoProng/runMesonGrid.C [new file with mode: 0644]
PWGUD/PWGUDdiffractiveLinkDef.h

index 27f4a19..266ed7a 100644 (file)
@@ -30,6 +30,10 @@ set ( SRCS
   DIFFRACTIVE/example/AliAnalysisTaskCDex.cxx
   DIFFRACTIVE/example/AliCDMesonUtilsStripped.cxx
   DIFFRACTIVE/example/AliCDMesonBaseStripped.cxx
+  DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.cxx
+  DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.cxx
+  DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.cxx
+  DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.cxx
   DIFFRACTIVE/xsAndTwoProng/AliCDMesonTracks.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.cxx b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.cxx
new file mode 100644 (file)
index 0000000..7e622b8
--- /dev/null
@@ -0,0 +1,2134 @@
+/*************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Select events according to gap conditions, analyze two track events in pp
+// collisions
+//
+// Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TList.h>
+#include <THnSparse.h>
+
+#include <TRandom3.h>
+
+#include "AliAODInputHandler.h"
+#include "AliAODHandler.h"
+#include "AliESDInputHandler.h"
+#include "AliPIDResponse.h"
+#include "AliPhysicsSelection.h"
+#include "AliGenEventHeader.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliMCEvent.h"
+#include "AliESDtrack.h"
+#include "AliStack.h"
+
+#include "AliCDMesonBase.h"
+#include "AliCDMesonUtils.h"
+#include "AliCDMesonTracks.h"
+#include "AliAnalysisTaskCDMeson.h"
+
+
+//------------------------------------------------------------------------------
+AliAnalysisTaskCDMeson::AliAnalysisTaskCDMeson(const char* name, Long_t state):
+       AliAnalysisTaskSE(name)
+       , fDoAOD(kFALSE)
+       , fAnalysisStatus(state)
+       , fAnalysisStatusString(new TObjString(Form("0x%lx", fAnalysisStatus)))
+       , fNoGapFraction(0.01)
+       , fReducedGapFraction(0.02)
+       , fMaxVtxDst(0.5) // value to be checked with the vertex study histograms
+
+       , fESDEvent(0x0)
+       , fAODEvent(0x0)
+       , fPIDResponse(0x0)
+       , fPhysicsSelection(0x0)
+       , fTracks(new AliCDMesonTracks())
+       , fVtxDst(-1.)
+       , fVtxZ(-20)
+       , fResidualTracks(0)
+       , fResidualTracklets(0)
+       , fMCprocessType(0)
+       , fMCprocess(-1)
+
+       , fRun(-999)
+       , fPIDmode(0)
+
+       , fTheta(0.)
+       , fPhi(0.)
+       , fMass(0.)
+       , fCurrentGapCondition(0)
+
+       , fGapRun(0x0)
+       , fHist(0x0)
+       , fThnMother(0x0)
+       , fThnMotherSoft(0x0)
+       , fThnMultiplicity(0x0)
+       , fThnMC(0x0)
+       , fThnEmptyEvents(0x0)
+       , fPWAtree(0x0)
+
+       , fv0ntrk(0x0)
+       , fv0fmdntrk(0x0)
+       , fv0fmdspdntrk(0x0)
+       , fv0fmdspdtpcntrk(0x0)
+
+       , fv0Rmntrk(0x0)
+       , fv0fmdRmntrk(0x0)
+       , fv0fmdspdRmntrk(0x0)
+       , fv0fmdspdtpcRmntrk(0x0)
+
+       , fMultStudy(0x0)
+       , fMultStudyV0dg(0x0)
+       , fMultStudyV0FMDdg(0x0)
+       , fMultStudyV0FMDSPDdg(0x0)
+       , fMultStudyV0FMDSPDTPCdg(0x0)
+
+       , fMultResponseMC(0x0)
+       , fMultResponseV0dgMC(0x0)
+       , fMultResponseV0FMDdgMC(0x0)
+       , fMultResponseV0FMDSPDdgMC(0x0)
+       , fMultResponseV0FMDSPDTPCdgMC(0x0)
+
+       , fMultRegionsMC(0x0)
+
+       , fhspdV0dg(0x0)
+       , fhspdV0FMDdg(0x0)
+       , fhspdV0FMDSPDdg(0x0)
+       , fhspdV0FMDSPDTPCdg(0x0)
+       , fhspdAfterCuts(0x0)
+
+       , fGapResponseMCv0Dg(0x0)
+       , fGapResponseMCv0fmdDg(0x0)
+       , fGapResponseMCv0fmdspdDg(0x0)
+       , fGapResponseMCv0fmdspdtpcDg(0x0)
+
+       , fhspd(0x0)
+       , fhfo(0x0)
+       , fhpriVtxDist(0x0)
+       , fhpriVtxPos(0x0)
+       , fhpriv(0x0)
+       , fhntrk(0x0)
+       , fhStatsFlow(0x0)
+       , fhFOchans(0x0)
+
+       , fVZEROhists(0x0)
+
+       , fv0Map(0x0)
+       , fv0fmdMap(0x0)
+       , fv0fmdspdMap(0x0)
+       , fv0fmdspdtpcMap(0x0)
+
+       , fv0MapCutted(0x0)
+       , fv0fmdMapCutted(0x0)
+       , fv0fmdspdMapCutted(0x0)
+       , fv0fmdspdtpcMapCutted(0x0)
+
+       , fHitMapSPDinner(0x0)
+       , fHitMapSPDouter(0x0)
+       , fHitMapSPDtrklt(0x0)
+       , fHitMapFMDa(0x0)
+       , fHitMapFMDc(0x0)
+
+       , fFMDsum1I(0x0)
+       , fFMDsum2I(0x0)
+       , fFMDsum2O(0x0)
+       , fFMDsum3I(0x0)
+       , fFMDsum3O(0x0)
+
+       , fPriVtxX(0x0)
+       , fPriVtxY(0x0)
+       , fPriVtxZ(0x0)
+       , fPriVtxDst(0x0)
+       , fPriVtxDstV0dg(0x0)
+       , fPriVtxDstV0FMDdg(0x0)
+       , fPriVtxDstV0FMDSPDdg(0x0)
+       , fPriVtxDstV0FMDSPDTPCdg(0x0)
+
+       , fTPCGapDCAaSide(0x0)
+       , fTPCGapDCAcSide(0x0)
+
+       , fComb2trkPIDuls(0x0)
+       , fComb2trkPIDls(0x0)
+       , fComb2trkPIDulsDG(0x0)
+       , fComb2trkPIDlsDG(0x0)
+       , fComb2trkPIDulsMC(0x0)
+       , fComb2trkPIDlsMC(0x0)
+       , fComb2trkPIDulsDGmc(0x0)
+       , fComb2trkPIDlsDGmc(0x0)
+       , fMCProcessUls(0x0)
+       , fMCProcessLs(0x0)
+
+       , fMAllTrackMass(0x0)
+{
+       //
+       // standard constructor (the one which should be used)
+       //
+       // slot in TaskSE must start from 1
+       Int_t iOutputSlot = 1;
+       DefineOutput(iOutputSlot++, TList::Class());
+       if (!(fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) {
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMother)) {
+                       DefineOutput(iOutputSlot++, THnSparse::Class());
+               }
+               if (!fAnalysisStatus
+                   || ((fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)
+                       && (fAnalysisStatus & AliCDMesonBase::kBitTHnMother))) {
+                       DefineOutput(iOutputSlot++, THnSparse::Class());
+               }
+               if (fAnalysisStatus & AliCDMesonBase::kBitMultStudy) {
+                       DefineOutput(iOutputSlot++, THnSparse::Class());
+               }
+               if (fAnalysisStatus & AliCDMesonBase::kBitTHnMC) {
+                       DefineOutput(iOutputSlot++, THnSparse::Class());
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitPWAtree)) {
+                       DefineOutput(iOutputSlot++, TTree::Class());
+               }
+       }
+       else { // create empty event THnSparse
+               DefineOutput(iOutputSlot++, THnSparse::Class());
+       }
+
+       if (fAnalysisStatus & AliCDMesonBase::kBitEEStudy) { // empty event study
+               fPhysicsSelection = new AliPhysicsSelection;
+               fPhysicsSelection->SetSkipTriggerClassSelection();
+               // accept all (A,C,E and I) triggers of a certain class
+       }
+
+       for (Int_t iGap = 0; iGap < kMax; ++iGap) {
+               fGapInformation[iGap] = 0;
+       }
+
+       // reset track pair pointers
+       fTrkPair[0] = 0x0;
+       fTrkPair[1] = 0x0;
+}
+
+
+//------------------------------------------------------------------------------
+AliAnalysisTaskCDMeson::AliAnalysisTaskCDMeson():
+       AliAnalysisTaskSE()
+       , fDoAOD(kFALSE)
+       , fAnalysisStatus(0x0)
+       , fAnalysisStatusString(0x0)
+       , fNoGapFraction(0.01)
+       , fReducedGapFraction(0.02)
+       , fMaxVtxDst(0.5)
+
+       , fESDEvent(0x0)
+       , fAODEvent(0x0)
+       , fPIDResponse(0x0)
+       , fPhysicsSelection(0x0)
+       , fTracks(0x0)
+       , fVtxDst(-1.)
+       , fVtxZ(-20)
+       , fResidualTracks(0)
+       , fResidualTracklets(0)
+       , fMCprocessType(0)
+       , fMCprocess(-1)
+
+       , fRun(-999)
+       , fPIDmode(0)
+
+       , fTheta(0.)
+       , fPhi(0.)
+       , fMass(0.)
+       , fCurrentGapCondition(0)
+
+       , fGapRun(0x0)
+       , fHist(0x0)
+       , fThnMother(0x0)
+       , fThnMotherSoft(0x0)
+       , fThnMultiplicity(0x0)
+       , fThnMC(0x0)
+       , fThnEmptyEvents(0x0)
+       , fPWAtree(0x0)
+
+       , fv0ntrk(0x0)
+       , fv0fmdntrk(0x0)
+       , fv0fmdspdntrk(0x0)
+       , fv0fmdspdtpcntrk(0x0)
+
+       , fv0Rmntrk(0x0)
+       , fv0fmdRmntrk(0x0)
+       , fv0fmdspdRmntrk(0x0)
+       , fv0fmdspdtpcRmntrk(0x0)
+
+       , fMultStudy(0x0)
+       , fMultStudyV0dg(0x0)
+       , fMultStudyV0FMDdg(0x0)
+       , fMultStudyV0FMDSPDdg(0x0)
+       , fMultStudyV0FMDSPDTPCdg(0x0)
+
+       , fMultResponseMC(0x0)
+       , fMultResponseV0dgMC(0x0)
+       , fMultResponseV0FMDdgMC(0x0)
+       , fMultResponseV0FMDSPDdgMC(0x0)
+       , fMultResponseV0FMDSPDTPCdgMC(0x0)
+
+       , fMultRegionsMC(0x0)
+
+       , fhspdV0dg(0x0)
+       , fhspdV0FMDdg(0x0)
+       , fhspdV0FMDSPDdg(0x0)
+       , fhspdV0FMDSPDTPCdg(0x0)
+       , fhspdAfterCuts(0x0)
+
+       , fGapResponseMCv0Dg(0x0)
+       , fGapResponseMCv0fmdDg(0x0)
+       , fGapResponseMCv0fmdspdDg(0x0)
+       , fGapResponseMCv0fmdspdtpcDg(0x0)
+
+       , fhspd(0x0)
+       , fhfo(0x0)
+       , fhpriVtxDist(0x0)
+       , fhpriVtxPos(0x0)
+       , fhpriv(0x0)
+       , fhntrk(0x0)
+       , fhStatsFlow(0x0)
+       , fhFOchans(0x0)
+
+       , fVZEROhists(0x0)
+
+       , fv0Map(0x0)
+       , fv0fmdMap(0x0)
+       , fv0fmdspdMap(0x0)
+       , fv0fmdspdtpcMap(0x0)
+
+       , fv0MapCutted(0x0)
+       , fv0fmdMapCutted(0x0)
+       , fv0fmdspdMapCutted(0x0)
+       , fv0fmdspdtpcMapCutted(0x0)
+
+       , fHitMapSPDinner(0x0)
+       , fHitMapSPDouter(0x0)
+       , fHitMapSPDtrklt(0x0)
+       , fHitMapFMDa(0x0)
+       , fHitMapFMDc(0x0)
+
+       , fFMDsum1I(0x0)
+       , fFMDsum2I(0x0)
+       , fFMDsum2O(0x0)
+       , fFMDsum3I(0x0)
+       , fFMDsum3O(0x0)
+
+       , fPriVtxX(0x0)
+       , fPriVtxY(0x0)
+       , fPriVtxZ(0x0)
+       , fPriVtxDst(0x0)
+       , fPriVtxDstV0dg(0x0)
+       , fPriVtxDstV0FMDdg(0x0)
+       , fPriVtxDstV0FMDSPDdg(0x0)
+       , fPriVtxDstV0FMDSPDTPCdg(0x0)
+
+       , fTPCGapDCAaSide(0x0)
+       , fTPCGapDCAcSide(0x0)
+
+       , fComb2trkPIDuls(0x0)
+       , fComb2trkPIDls(0x0)
+       , fComb2trkPIDulsDG(0x0)
+       , fComb2trkPIDlsDG(0x0)
+       , fComb2trkPIDulsMC(0x0)
+       , fComb2trkPIDlsMC(0x0)
+       , fComb2trkPIDulsDGmc(0x0)
+       , fComb2trkPIDlsDGmc(0x0)
+       , fMCProcessUls(0x0)
+       , fMCProcessLs(0x0)
+
+       , fMAllTrackMass(0x0)
+{
+       //
+       // default constructor (should not be used)
+       //
+
+       fTrkPair[0] = 0x0;
+       fTrkPair[1] = 0x0;
+}
+
+
+//------------------------------------------------------------------------------
+AliAnalysisTaskCDMeson::~AliAnalysisTaskCDMeson()
+{
+       //
+       //Destructor
+       //
+       //delete fESDEvent;
+       //delete fESDpid;
+
+       if (!(fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) { // normal operation
+               if ((!fAnalysisStatus || (fAnalysisStatus && AliCDMesonBase::kBitTHnMother))
+                   && fThnMother
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fThnMother;
+                       fThnMother = 0x0;
+               }
+               if ((!fAnalysisStatus
+                    || ((fAnalysisStatus && AliCDMesonBase::kBitSoftTracks)
+                        && (fAnalysisStatus && AliCDMesonBase::kBitTHnMother)))
+                   && fThnMotherSoft
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fThnMotherSoft;
+                       fThnMotherSoft = 0x0;
+               }
+               if ((!fAnalysisStatus
+                    || (fAnalysisStatus && AliCDMesonBase::kBitMultStudy))
+                   && fThnMultiplicity
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fThnMultiplicity;
+                       fThnMultiplicity = 0x0;
+               }
+               if ((!fAnalysisStatus || (fAnalysisStatus && AliCDMesonBase::kBitTHnMC))
+                   && fThnMC
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fThnMC;
+                       fThnMC = 0x0;
+               }
+       }
+       else { // empty event study
+               if (fPhysicsSelection
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fPhysicsSelection;
+                       fPhysicsSelection = 0x0;
+               }
+               if (fThnEmptyEvents
+                   && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+                       AliAnalysisManager::kProofAnalysis)) {
+                       delete fThnEmptyEvents;
+                       fThnEmptyEvents = 0x0;
+               }
+       }
+
+       /*
+       if (fHist
+           && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+               AliAnalysisManager::kProofAnalysis)) {
+               fHist->Clear();
+               delete fHist;
+               fHist = 0x0;
+       }
+       */
+
+       if (fTracks
+           && (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
+               AliAnalysisManager::kProofAnalysis)) {
+               delete fTracks;
+               fTracks = 0x0;
+       }
+
+       if (fVZEROhists) {
+               delete fVZEROhists;
+               fVZEROhists = 0x0;
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::UserCreateOutputObjects()
+{
+       //
+       //CreateOutputObjects
+       //
+
+       //= THnSparse ================================================================
+       if (!(fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) { // normal operation
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMother)) {
+                       fThnMother = AliCDMesonBase::GetThnMother(/*"CDMeson_Mother"*/);
+               }
+               if (!fAnalysisStatus
+                   || ((fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)
+                       && (fAnalysisStatus & AliCDMesonBase::kBitTHnMother))) {
+                       fThnMotherSoft =  AliCDMesonBase::GetThnMother("CDMeson_MotherSoft");
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitMultStudy)){
+                       fThnMultiplicity =  AliCDMesonBase::GetThnMultiplicity();
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMC)) {
+                       fThnMC = AliCDMesonBase::GetThnMother("CDMeson_MotherMC");
+               }
+       }
+       else { // empty event studies
+               fThnEmptyEvents = AliCDMesonBase::GetThnEmptyEvents();
+       }
+
+       //= TList for Histograms =====================================================
+       fHist = new TList;
+       fHist->SetOwner(); // ensures that the histograms are all deleted on exit!
+
+       //fHist->Add(fAnalysisStatusString); // WARNING: CRASHES MERGING
+
+       //= GAP RUN =
+       if (!(fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) { // normal operation
+               // LHC10 data
+               // runa1=114649, runb1=117630, runc0=117631, runc1=121526, rund0=121527,
+               // rund1=126460, rune0=126461, rune1 = 130930; runf0=130931
+               // all checked on elog
+               // LHC11 data
+               //const Int_t rune0 = 160676;
+               //const Int_t runf1 = 165462;
+
+               // LHC10 + LHC11 data
+               const Int_t runb0 = 114650; // first run of LHC10b
+               const Int_t runf1 = 165462; // last run of LHC11f
+
+               const Int_t run0 = runb0-1, run1 = runf1 + 1;
+
+               const Int_t nrun = (Int_t)(run1-run0);
+
+               const Int_t bins[2] = {nrun, (Double_t)AliCDMesonBase::kBitGapMax - 1.};
+               const Double_t xmin[2] = {(Double_t)run0, 1.};
+               const Double_t xmax[2] = {
+                       (Double_t)run1, (Double_t)AliCDMesonBase::kBitGapMax
+               };
+
+               fGapRun = new THnSparseI("GapRun", ";run number;gap condition", 2, bins,
+                                        xmin, xmax);
+               fHist->Add(fGapRun);
+       }
+
+       //= MULTIPLICITY PER GAP CONDITION =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitMultPerGapHists)) {
+               fv0ntrk = new TH2D("b00_v0ntrk", ";number of tracks;gap condition",
+                                  80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0
+               fHist->Add(fv0ntrk);
+
+               fv0fmdntrk = new TH2D("b01_v0fmdntrk", ";number of tracks;gap condition",
+                                     80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMD
+               fHist->Add(fv0fmdntrk);
+
+               fv0fmdspdntrk =
+                       new TH2D("b02_v0fmdspdntrk", ";number of tracks;gap condition",
+                        80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMDSPD
+               fHist->Add(fv0fmdspdntrk);
+
+               fv0fmdspdtpcntrk =
+                       new TH2D("b03_v0fmdspdtpcntrk", ";number of tracks;gap condition",
+                                80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMDSPDTPC
+               fHist->Add(fv0fmdspdtpcntrk);
+       }
+
+       //= MULTIPLICITY REMOVED PER GAP CONDITION =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitRmMultPerGapHists)) {
+               fv0Rmntrk = new TH2F("b05_v0Rmntrk",
+                                    ";number of removed tracks due to cuts;gap condition",
+                                    80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0
+               fHist->Add(fv0Rmntrk);
+
+               fv0fmdRmntrk =
+                       new TH2F("b06_v0fmdRmntrk",
+                                ";number of removed tracks due to cuts;gap condition",
+                                80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMD
+               fHist->Add(fv0fmdRmntrk);
+
+               fv0fmdspdRmntrk =
+                       new TH2F("b07_v0fmdspdRmntrk",
+                                ";number of removed tracks due to cuts;gap condition",
+                                80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMDSPD
+               fHist->Add(fv0fmdspdRmntrk);
+
+               fv0fmdspdtpcRmntrk =
+                       new TH2F("b08_v0fmdspdtpcRmntrk",
+                                ";number of removed tracks due to cuts;gap condition",
+                                80, 0., 80., 4, 1., 5.);
+               //x: ntrk; y: V0FMDSPDTPC
+               fHist->Add(fv0fmdspdtpcRmntrk);
+       }
+
+       //= SOFT TRACK INFLUENCE =
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)) {
+               fMultStudy = new TH2F("b10_MultStudySoftInfluence",
+                                     ";multiplicity;multiplicity (including ITSsa)",
+                                     20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultStudy);
+
+               fMultStudyV0dg =
+                       new TH2F("b11_MultStudySoftInfluence_V0_DoubleGap",
+                                ";multiplicity;multiplicity (including ITSsa)",
+                                20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultStudyV0dg);
+
+               fMultStudyV0FMDdg =
+                       new TH2F("b12_MultStudySoftInfluence_V0FMD_DoubleGap",
+                                ";multiplicity;multiplicity (including ITSsa)",
+                                20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultStudyV0FMDdg);
+
+               fMultStudyV0FMDSPDdg
+                       = new TH2F("b13_MultStudySoftInfluence_V0FMDSPD_DoubleGap",
+                                  ";multiplicity;multiplicity (including ITSsa)",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultStudyV0FMDSPDdg);
+
+               fMultStudyV0FMDSPDTPCdg
+                       = new TH2F("b14_MultStudySoftInfluence_V0FMDSPDTPC_DoubleGap",
+                                  ";multiplicity;multiplicity (including ITSsa)",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultStudyV0FMDSPDTPCdg);
+       }
+
+       //= MULTIPLICITY RESPONSE DATA VS MC =
+       if (fAnalysisStatus & AliCDMesonBase::kBitMultResponseMC) {
+               fMultResponseMC
+                       = new TH2F("b15_MultResponseMC",
+                                  ";multiplicity from MC truth; reconstructed multiplicity",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultResponseMC);
+
+               fMultResponseV0dgMC
+                       = new TH2F("b16_MultResponseV0dgMC",
+                                  ";multiplicity from MC truth; reconstructed multiplicity",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultResponseV0dgMC);
+
+               fMultResponseV0FMDdgMC
+                       = new TH2F("b17_MultResponseV0FMDdgMC",
+                                  ";multiplicity from MC truth; reconstructed multiplicity",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultResponseV0FMDdgMC);
+
+               fMultResponseV0FMDSPDdgMC
+                       = new TH2F("b18_MultResponseV0FMDSPDdgMC",
+                                  ";multiplicity from MC truth; reconstructed multiplicity",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultResponseV0FMDSPDdgMC);
+
+               fMultResponseV0FMDSPDTPCdgMC
+                       = new TH2F("b19_MultResponseV0FMDSPDTPCdgMC",
+                                  ";multiplicity from MC truth; reconstructed multiplicity",
+                                  20, 0., 20., 20, 0., 20.);
+               fHist->Add(fMultResponseV0FMDSPDTPCdgMC);
+
+               fMultRegionsMC
+                       = new TH2F("b20_fMultRegionsMC", ";multiplicity;",
+                                  50, 0., 50., 4, 0., 4.);
+               TAxis* a = fMultRegionsMC->GetYaxis();
+               a->SetBinLabel(1, "A-side");
+               a->SetBinLabel(2, "Center");
+               a->SetBinLabel(3, "C-side");
+               a->SetBinLabel(4, "A and C side");
+               fHist->Add(fMultRegionsMC);
+
+               fGapResponseMCv0Dg =
+                       new TH1I("b26_fGapResponseMCv0Dg", ";;counts", 3, 0., 3.);
+               a = fGapResponseMCv0Dg->GetXaxis();
+               a->SetBinLabel(1, "DG in MC, no DG in data");
+               a->SetBinLabel(2, "DG in both MC and data");
+               a->SetBinLabel(3, "DG in data, no DG in MC");
+               fHist->Add(fGapResponseMCv0Dg);
+               fGapResponseMCv0fmdDg =
+                       new TH1I("b27_fGapResponseMCv0fmdDg", ";;counts", 3, 0., 3.);
+               a = fGapResponseMCv0fmdDg->GetXaxis();
+               a->SetBinLabel(1, "DG in MC, no DG in data");
+               a->SetBinLabel(2, "DG in both MC and data");
+               a->SetBinLabel(3, "DG in data, no DG in MC");
+               fHist->Add(fGapResponseMCv0fmdDg);
+               fGapResponseMCv0fmdspdDg =
+                       new TH1I("b28_fGapResponseMCv0fmdspdDg", ";;counts", 3, 0., 3.);
+               a = fGapResponseMCv0fmdspdDg->GetXaxis();
+               a->SetBinLabel(1, "DG in MC, no DG in data");
+               a->SetBinLabel(2, "DG in both MC and data");
+               a->SetBinLabel(3, "DG in data, no DG in MC");
+               fHist->Add(fGapResponseMCv0fmdspdDg);
+               fGapResponseMCv0fmdspdtpcDg
+                       = new TH1I("b29_fGapResponseMCv0fmdspdtpcDg", ";;counts", 3, 0., 3.);
+               a = fGapResponseMCv0fmdspdtpcDg->GetXaxis();
+               a->SetBinLabel(1, "DG in MC, no DG in data");
+               a->SetBinLabel(2, "DG in both MC and data");
+               a->SetBinLabel(3, "DG in data, no DG in MC");
+               fHist->Add(fGapResponseMCv0fmdspdtpcDg);
+       }
+
+       //= FAST-OR MULTIPLICITY STUDY =
+       if (fAnalysisStatus & AliCDMesonBase::kBitFastORmultStudy) {
+               fhspdV0dg =
+                       new TH1I("b21_spd_V0dg",
+                                "V0 double gap;fired SPD FastOR chips; counts", 50, 0., 50.);
+               fHist->Add(fhspdV0dg);
+               fhspdV0FMDdg =
+                       new TH1I("b22_spd_V0FMDdg",
+                                "V0-FMD double gap;fired SPD FastOR chips; counts", 50, 0., 50.);
+               fHist->Add(fhspdV0FMDdg);
+               fhspdV0FMDSPDdg =
+                       new TH1I("b23_spd_V0FMDSPDdg",
+                                "V0-FMD-SPD double gap;fired SPD FastOR chips; counts",
+                                50, 0., 50.);
+               fHist->Add(fhspdV0FMDSPDdg);
+               fhspdV0FMDSPDTPCdg =
+                       new TH1I("b24_spd_V0FMDSPDTPCdg",
+                                "V0-FMD-SPD-TPC double gap;fired SPD FastOR chips; counts",
+                                50, 0., 50.);
+               fHist->Add(fhspdV0FMDSPDTPCdg);
+               fhspdAfterCuts =
+                       new TH1I("b25_spd_afterCuts",
+                                "after cuts;fired SPD FastOR chips; counts",
+                                50, 0., 50.);
+               fHist->Add(fhspdAfterCuts);
+       }
+
+       //= STATISTICS FLOW =
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhspd = new TH1I("a00_spd", ";fired SPD FastOR chips; counts", 50, 0., 50.);
+               fHist->Add(fhspd);
+
+               if (fAnalysisStatus & AliCDMesonBase::kBitFastORStudy) {
+                       fhfo =
+                               new TH2I("a05_fo", ";fired chips inner layer;fired chips outer layer",
+                                        50, 0., 50., 50, 0., 50.);
+                       fHist->Add(fhfo);
+
+                       fhFOchans =
+                               new TH1F("a06_foChans", ";channel key; counts", 1200, 0., 1200.);
+                       fHist->Add(fhFOchans);
+               }
+
+               fhpriVtxDist = new TH1I("a08_priVtxDist", ";vertex position (cm); counts",
+                                       601, -30.05, 30.05);
+               fHist->Add(fhpriVtxDist);
+
+               fhpriVtxPos = new TH1I("a09_priVtxPos", ";vertex position (boolean);counts",
+                                 2, 0., 2.);
+               fhpriVtxPos->GetXaxis()->SetBinLabel(1,"inside");
+               fhpriVtxPos->GetXaxis()->SetBinLabel(2,"outside");
+               fHist->Add(fhpriVtxPos);
+
+
+               fhpriv = new TH1I("a10_priv", ";vertex quality (boolean);counts",
+                                 2, 0., 2.);
+               fhpriv->GetXaxis()->SetBinLabel(1,"bad vertex");
+               fhpriv->GetXaxis()->SetBinLabel(2,"good vertex");
+               fHist->Add(fhpriv);
+
+               fhntrk = new TH1I("a20_ntrk", ";track multiplicity after cuts;counts",
+                                 10, 0., 10.);
+               fHist->Add(fhntrk);
+
+               fhStatsFlow = AliCDMesonBase::GetHistStatsFlow();
+               fHist->Add(fhStatsFlow);
+       }
+
+       //= ETA-PHI MAPS FOR TRACKS =
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitEtaPhiMaps)) {
+               fv0Map = new TH2F("d0_v0Map", ";#eta;#phi", 40, -2., 2.,
+                                 60, 0., TMath::TwoPi());
+               fHist->Add(fv0Map);
+
+               fv0fmdMap = new TH2F("d0_v0fmdMap", ";#eta;#phi", 40, -2., 2.,
+                                    60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdMap);
+
+               fv0fmdspdMap = new TH2F("d0_v0fmdspdMap", ";#eta;#phi", 40, -2., 2.,
+                                       60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdspdMap);
+
+               fv0fmdspdtpcMap = new TH2F("d0_v0fmdspdtpcMap", ";#eta;#phi", 30, -1.5, 1.5,
+                                          60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdspdtpcMap);
+       }
+
+       //= ETA-PHI MAPS FOR CUTTED TRACKS =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitEtaPhiMapsWithCuts)) {
+               fv0MapCutted = new TH2F("d0_v0MapCutted", ";#eta;#phi", 24, -1.2, 1.2,
+                                       60, 0., TMath::TwoPi());
+               fHist->Add(fv0MapCutted);
+
+               fv0fmdMapCutted = new TH2F("d0_v0fmdMapCutted", ";#eta;#phi", 24, -1.2, 1.2,
+                                          60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdMapCutted);
+
+               fv0fmdspdMapCutted = new TH2F("d0_v0fmdspdMapCutted", ";#eta;#phi",
+                                             24, -1.2, 1.2, 60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdspdMapCutted);
+
+               fv0fmdspdtpcMapCutted = new TH2F("d0_v0fmdspdtpcMapCutted", ";#eta;#phi",
+                                                24, -1.2, 1.2, 60, 0., TMath::TwoPi());
+               fHist->Add(fv0fmdspdtpcMapCutted);
+       }
+
+       //= SPD HIT MAPS =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitHitMapSPD) ||
+           (fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) {
+               fHitMapSPDinner = new TH2F("d1_HitMapSPDinner", ";#eta;#phi", 72, -3.6, 3.6,
+                                          20, 0., TMath::TwoPi());
+               fHist->Add(fHitMapSPDinner);
+
+               fHitMapSPDouter = new TH2F("d2_HitMapSPDouter", ";#eta;#phi", 48, -2.4, 2.4,
+                                          40, 0., TMath::TwoPi());
+               fHist->Add(fHitMapSPDouter);
+
+               fHitMapSPDtrklt = new TH2F("d3_HitMapSPDtrklt", ";#eta;#phi", 60, -3., 3.,
+                                          72, 0., TMath::TwoPi());
+               fHist->Add(fHitMapSPDtrklt);
+       }
+
+       //= FMD HIT MAPS =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitHitMapFMD) ||
+           (fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) {
+               fHitMapFMDa = new TH2F("d4_HitMapFMDa", ";#eta;#phi",
+                                      36, 1.5, 5.1, 40, 0., TMath::TwoPi());
+               fHist->Add(fHitMapFMDa);
+               fHitMapFMDc = new TH2F("d4_HitMapFMDc", ";#eta;#phi",
+                                      20, -3.5, -1.5, 40, 0., TMath::TwoPi());
+               fHist->Add(fHitMapFMDc);
+       }
+
+       //= FMD SUMMATION =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitFMDsum) ||
+           (fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) {
+               fFMDsum1I = new TH1F("d5_FMDsum1I", ";sum;counts", 50, 0., 50.);
+               fHist->Add(fFMDsum1I);
+               fFMDsum2I = new TH1F("d5_FMDsum2I", ";sum;counts", 50, 0., 50.);
+               fHist->Add(fFMDsum2I);
+               fFMDsum2O = new TH1F("d5_FMDsum2O", ";sum;counts", 50, 0., 50.);
+               fHist->Add(fFMDsum2O);
+               fFMDsum3I = new TH1F("d5_FMDsum3I", ";sum;counts", 50, 0., 50.);
+               fHist->Add(fFMDsum3I);
+               fFMDsum3O = new TH1F("d5_FMDsum3O", ";sum;counts", 50, 0., 50.);
+               fHist->Add(fFMDsum3O);
+       }
+
+       //= VERTEX STUDIES =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitVtxStudies)) {
+               fPriVtxX = new TH1I("e0_PriVtxX", ";x (cm);counts", 100, -10., 10.);
+               fHist->Add(fPriVtxX);
+               fPriVtxY = new TH1I("e0_PriVtxY", ";y (cm);counts", 100, -10., 10.);
+               fHist->Add(fPriVtxY);
+               fPriVtxZ = new TH1I("e0_PriVtxZ", ";z (cm);counts", 100, -10., 10.);
+               fHist->Add(fPriVtxZ);
+
+               fPriVtxDst = new TH1I("e1_PriVtxDst", ";dst (cm);counts", 101, 0., 2.525);
+               fHist->Add(fPriVtxDst);
+
+               fPriVtxDstV0dg =
+                       new TH1I("e2_PriVtxDstV0dg", ";dst (cm);counts", 101, 0., 2.525);
+               fHist->Add(fPriVtxDstV0dg);
+
+               fPriVtxDstV0FMDdg =
+                       new TH1I("e3_PriVtxDstV0FMDdg", ";dst (cm);counts", 101, 0., 2.525);
+               fHist->Add(fPriVtxDstV0FMDdg);
+
+               fPriVtxDstV0FMDSPDdg =
+                       new TH1I("e4_PriVtxDstV0FMDSPDdg", ";dst (cm);counts", 101, 0., 2.525);
+               fHist->Add(fPriVtxDstV0FMDSPDdg);
+
+               fPriVtxDstV0FMDSPDTPCdg =
+                       new TH1I("e5_PriVtxDstV0FMDSPDTPCdg", ";dst (cm);counts", 101, 0., 2.525);
+               fHist->Add(fPriVtxDstV0FMDSPDTPCdg);
+       }
+
+       //= TPC GAP STUDY =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitTPCGapStudy)) {
+               if (!fDoAOD) {
+                       fTPCGapDCAaSide = new TH2F("i00_TPCGapDCAaSide", ";dca_{xy};dca_{z}",
+                                                  2000, 0., 1., 1000, 0., 5.);
+                       fHist->Add(fTPCGapDCAaSide);
+                       fTPCGapDCAcSide = new TH2F("i01_TPCGapDCAcSide", ";dca_{xy};dca_{z}",
+                                                  2000, 0., 1., 1000, 0., 5.);
+                       fHist->Add(fTPCGapDCAcSide);
+               }
+       }
+
+       //= VZERO TRIGGER STUDY =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitVZEROStudy)) {
+               if (!fDoAOD) { // not yet possible for AODs
+                       fVZEROhists = AliCDMesonBase::GetHistVZEROStudies(fHist);
+               }
+       }
+
+       //= PID MATRIX FOR DATA AND MC =
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitPIDStudy)) {
+               fComb2trkPIDuls = AliCDMesonBase::GetHistPIDStudies("f0_Comb2trkPIDuls");
+               fHist->Add(fComb2trkPIDuls);
+
+               fComb2trkPIDls = AliCDMesonBase::GetHistPIDStudies("f1_Comb2trkPIDls");
+               fHist->Add(fComb2trkPIDls);
+
+               fComb2trkPIDulsDG =
+                       AliCDMesonBase::GetHistPIDStudies("f2_Comb2trkPIDulsDG");
+               fHist->Add(fComb2trkPIDulsDG);
+
+               fComb2trkPIDlsDG = AliCDMesonBase::GetHistPIDStudies("f3_Comb2trkPIDlsDG");
+               fHist->Add(fComb2trkPIDlsDG);
+
+               if (fAnalysisStatus & AliCDMesonBase::kBitTHnMC) {
+                       fComb2trkPIDulsMC =
+                               AliCDMesonBase::GetHistPIDStudies("f4_Comb2trkPIDulsMC");
+                       fHist->Add(fComb2trkPIDulsMC);
+
+                       fComb2trkPIDlsMC = AliCDMesonBase::GetHistPIDStudies("f5_Comb2trkPIDlsMC");
+                       fHist->Add(fComb2trkPIDlsMC);
+
+                       fComb2trkPIDulsDGmc =
+                               AliCDMesonBase::GetHistPIDStudies("f6_Comb2trkPIDulsDGmc");
+                       fHist->Add(fComb2trkPIDulsDGmc);
+
+                       fComb2trkPIDlsDGmc =
+                               AliCDMesonBase::GetHistPIDStudies("f7_Comb2trkPIDlsDGmc");
+                       fHist->Add(fComb2trkPIDlsDGmc);
+               }
+       }
+
+       //= MC PROCESS IDs =
+       if (!fAnalysisStatus || fAnalysisStatus & AliCDMesonBase::kBitMCProcess) {
+               fMCProcessUls = new TH2F("g0_MCProcessUls",
+                                        ";MC Process ID;double-gap topology",
+                                        105, 1., 106., 12, 0., 12.);
+               TAxis* axis = fMCProcessUls->GetYaxis();
+               axis->SetBinLabel(1, "V0");
+               axis->SetBinLabel(2, "FMD");
+               axis->SetBinLabel(3, "SPD");
+               axis->SetBinLabel(4, "TPC");
+               axis->SetBinLabel(5, "V0-FMD");
+               axis->SetBinLabel(6, "V0-FMD-SPD");
+               axis->SetBinLabel(7, "V0-FMD-SPD-TPC");
+               axis->SetBinLabel(8, "TPC-SPD");
+               axis->SetBinLabel(9, "TPC-SPD-FMD");
+               axis->SetBinLabel(10, "TPC-SPD-FMD-V0");
+               axis->SetBinLabel(11, "SPD-FMD");
+               axis->SetBinLabel(12, "SPD-FMD-V0");
+               fHist->Add(fMCProcessUls);
+
+               fMCProcessLs = new TH2F("g1_MCProcessLs",
+                                       ";MC Process ID;double-gap topology",
+                                       105, 1., 106., 12, 0., 12.);
+               axis = fMCProcessLs->GetYaxis();
+               axis->SetBinLabel(1, "V0");
+               axis->SetBinLabel(2, "FMD");
+               axis->SetBinLabel(3, "SPD");
+               axis->SetBinLabel(4, "TPC");
+               axis->SetBinLabel(5, "V0-FMD");
+               axis->SetBinLabel(6, "V0-FMD-SPD");
+               axis->SetBinLabel(7, "V0-FMD-SPD-TPC");
+               axis->SetBinLabel(8, "TPC-SPD");
+               axis->SetBinLabel(9, "TPC-SPD-FMD");
+               axis->SetBinLabel(10, "TPC-SPD-FMD-V0");
+               axis->SetBinLabel(11, "SPD-FMD");
+               axis->SetBinLabel(12, "SPD-FMD-V0");
+               fHist->Add(fMCProcessLs);
+       }
+
+       //= ALL TRACK INVARIANT MASS =================================================
+       if (fAnalysisStatus & AliCDMesonBase::kBitAllTrackMass) {
+               fMAllTrackMass = new TH1F("j00_AllTrackInvMass",
+                                         ";all track invariant mass (GeV);counts",
+                                         95, 0.25, 5.);
+               fHist->Add(fMAllTrackMass);
+       }
+       //= PWA TREE =================================================================
+       if (!fAnalysisStatus || fAnalysisStatus & AliCDMesonBase::kBitPWAtree) {
+               fPWAtree = new TTree("freidt_pwa", "pwa");
+               fPWAtree->Branch("theta", &fTheta);
+               fPWAtree->Branch("phi", &fPhi);
+               fPWAtree->Branch("m", &fMass);
+               fPWAtree->Branch("theta", &fTheta);
+               fPWAtree->Branch("pX", &fMomentum[0]);
+               fPWAtree->Branch("pY", &fMomentum[1]);
+               fPWAtree->Branch("pZ", &fMomentum[2]);
+               fPWAtree->Branch("gapCondition", &fCurrentGapCondition);
+       }
+
+       PostOutputs();
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::UserExec(Option_t *)
+{
+       //
+       // Executed for every event which passed the physics selection
+       // 
+
+       //printf("Entry: %ld\n", (Long_t)Entry()); // print current event number
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinTotalInput); // stats flow
+       }
+
+       //= INPUT DATA SANITY TESTS ==================================================
+       if (!CheckInput()) {
+               PostOutputs();
+               return;
+       }
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinGoodInput); // stats flow
+       }
+
+       //= ANALYZE ONLY PHOJET CD EVENTS ============================================
+       if (fAnalysisStatus & AliCDMesonBase::kBitCDonly) {
+               if (fMCEvent) {
+                       AliGenEventHeader* header = fMCEvent->GenEventHeader();
+                       if (header) {
+                               // Phojet
+                               if (TString(header->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
+                                       fMCprocess = ((AliGenDPMjetEventHeader*)header)->ProcessType();
+                                       if (fMCprocess != 4) { // no central diffraction
+                                               PostOutputs();
+                                               return;
+                                       }
+                                       else{ // CD found
+                                               fhStatsFlow->Fill(AliCDMesonBase::kBinCDonlyEvents); // statsflow
+                                       }
+                               }
+                       }
+               }
+       }
+
+
+       //= V0-AND/OR ================================================================
+       if (!fDoAOD) {
+               if (AliCDMesonUtils::V0AND(fESDEvent))
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinV0OR);
+               if (AliCDMesonUtils::V0OR(fESDEvent))
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinV0AND);
+       }
+
+       //= EMPTY EVENT STUDY ========================================================
+       // only implemented on ESDs!
+       if (fAnalysisStatus & AliCDMesonBase::kBitEEStudy) {
+               DoEmptyEventStudy();
+               PostOutputs();
+               return;
+       }
+
+       //= MC TRUTH =================================================================
+       // for now only implemented on ESDs
+       Int_t nMCprimaries = 0;
+       DetermineMCprocessType();
+       if (fAnalysisStatus & AliCDMesonBase::kBitTHnMC) {
+               nMCprimaries = DoMCTruth();
+       }
+
+       //= EVENT SELECTION ==========================================================
+       Int_t kfo =
+               (fAnalysisStatus & AliCDMesonBase::kBitFastORStudy & !fDoAOD) ? 1 : 0;
+       Int_t ninnerp=-999, nouterp=-999;
+       Bool_t eventIsValid = (fDoAOD) ?
+               AliCDMesonUtils::CutEvent(fAODEvent, fhpriv, fPriVtxX, fPriVtxY, fPriVtxZ,
+                                         fhpriVtxPos, fhpriVtxDist) :
+               AliCDMesonUtils::CutEvent(fESDEvent, fhspd, fhpriv, fhpriVtxPos,
+                                         fhpriVtxDist, fhfo, fhFOchans, kfo, ninnerp,
+                                         nouterp, fPriVtxX, fPriVtxY, fPriVtxZ);
+       if (!eventIsValid) {
+               PostOutputs();
+               return;
+       }
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinEventsAfterCuts); // stats flow
+       }
+
+       //= PILE UP ==================================================================
+       const Bool_t isPileup = (fDoAOD) ?
+               fAODEvent->IsPileupFromSPD(2, 0.8, 3., 2., 5.) :
+               fESDEvent->IsPileupFromSPD(2, 0.8, 3., 2., 5.);
+       // using only 2 instead of three contributors
+
+       if (isPileup) {
+               PostOutputs();
+               return;
+       }
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinEventsWithOutPileUp); // stats flow
+       }
+
+       //= VZERO TRIGGER STUDY ======================================================
+       // dtermine the VZERO signal distributions
+       if (fVZEROhists) {
+               AliCDMesonUtils::DoVZEROStudy(fESDEvent, fVZEROhists, fRun);
+       }
+
+       //= GAP ======================================================================
+       // determine the complete gap configuration (for all gap-tagging detectors)
+       if (!DetermineGap()) {
+               PostOutputs();
+               return;
+       }
+
+       // fill gaprun (determines the gap fraction per run)
+       const Double_t x[2] = {(Double_t)fRun, (Double_t)fCurrentGapCondition};
+       if (fGapRun) fGapRun->Fill(x);
+
+       //= VERTEX COINCIDENCE AND POSITION ==========================================
+       AnalyzeVtx();
+
+       //= TRACK CUTS ===============================================================
+       Bool_t doSoft = !fAnalysisStatus ||
+               (fAnalysisStatus & AliCDMesonBase::kBitSoftTracks);
+
+       fTracks->ProcessEvent(fAODEvent, fESDEvent, doSoft); // apply cuts
+       DoMultiplicityStudy(nMCprimaries); // fill corresponding histograms
+       if (fMAllTrackMass &&
+           (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG)) {
+               // calculate the invariant mass of all tracks in the event including soft
+               // ones if the event is a full double gap event
+               fMAllTrackMass->Fill(fTracks->GetInvariantMass());
+       }
+
+       // is multiplicity within the desired range of  2 to 3?
+       Int_t nch = fTracks->GetTracks();
+       Int_t ncombined = fTracks->GetCombinedTracks();
+       Bool_t wSoft = ((ncombined >= 2) && (ncombined <=3)); // including soft tracks
+       Bool_t woSoft = ((nch >= 2) && (nch <= 3)); // excluding soft tracks
+
+       if (!wSoft && !woSoft) {
+               // multiplicity out of range (both with soft and without soft track)
+               PostOutputs();
+               return;
+       }
+
+       //= TRACK PAIRS ==============================================================
+       // loop over all track combinations
+       for(Int_t ii=0; ii<ncombined; ii++){
+               for(Int_t jj=ii+1; jj<ncombined; jj++){
+                       // assign current tracks
+                       fTrkPair[0] = fTracks->GetTrack(ii);
+                       fTrkPair[1] = fTracks->GetTrack(jj);
+
+                       // analyze track pairs
+                       if (wSoft) DoTrackPair(kTRUE);
+                       if (woSoft && (ii < nch) && (jj < nch)) DoTrackPair(kFALSE);
+               }
+       }
+
+       //============================================================================
+       PostOutputs();
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::PostOutputs()
+{
+       //
+       // PostData
+       //
+
+       Int_t iOutputSlot = 1; // dynamic output slot number handling
+       PostData(iOutputSlot++, fHist);
+       if (!(fAnalysisStatus & AliCDMesonBase::kBitEEStudy)) { // normal operation
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMother)) {
+                       // THnSparse for invariant mass distributions
+                       PostData(iOutputSlot++, fThnMother);
+               }
+               if (!fAnalysisStatus
+                   || ((fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)
+                       && (fAnalysisStatus & AliCDMesonBase::kBitTHnMother))) {
+                       // same including soft tracks
+                       PostData(iOutputSlot++, fThnMotherSoft);
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitMultStudy)) {
+                       // multiplicity study
+                       PostData(iOutputSlot++, fThnMultiplicity);
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMC)) {
+                       // MC truth study is active
+                       PostData(iOutputSlot++, fThnMC);
+               }
+               if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitPWAtree)) {
+                       // PWA tree active
+                       PostData(iOutputSlot++, fPWAtree);
+               }
+       }
+       else if (fAnalysisStatus & AliCDMesonBase::kBitEEStudy) {
+               // empty event study
+               PostData(iOutputSlot++, fThnEmptyEvents);
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::FillEtaPhiMaps() const
+{
+       //
+       // Controls which event contributes to which map
+       //
+
+       AliVEvent* event = (fDoAOD) ? (AliVEvent*)fAODEvent : (AliVEvent*)fESDEvent;
+
+       if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdspdtpcMap,
+                                              fv0fmdspdtpcMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdspdMap,
+                                              fv0fmdspdMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdMap, fv0fmdMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0Map, fv0MapCutted);
+       }
+       else if (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdspdMap,
+                                              fv0fmdspdMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdMap, fv0fmdMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0Map, fv0MapCutted);
+       }
+       else if (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG) {
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0fmdMap, fv0fmdMapCutted);
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0Map, fv0MapCutted);
+       }
+       else if (fGapInformation[kV0] == AliCDMesonBase::kBinDG) {
+               AliCDMesonUtils::FillEtaPhiMap(event, fTracks, fv0Map, fv0MapCutted);
+       }
+}
+
+
+//------------------------------------------------------------------------------
+Bool_t AliAnalysisTaskCDMeson::CheckInput()
+{
+       //
+       //general protection
+       //
+       if (const AliESDInputHandler *esdH =
+           dynamic_cast<AliESDInputHandler*>(fInputHandler)){
+               fESDEvent = esdH->GetEvent();
+       }
+       else if (const AliAODInputHandler *aodH =
+                dynamic_cast<AliAODInputHandler*>(fInputHandler)){
+               fAODEvent = aodH->GetEvent();
+               fDoAOD = kTRUE;
+       }
+       fNoGapFraction = (fDoAOD) ? 1. : fNoGapFraction; // process all running on AOD
+       fPIDResponse = (AliPIDResponse*)fInputHandler->GetPIDResponse();
+
+       if(!fESDEvent && !fAODEvent){
+               printf("freidtlog No valid event\n");
+               return kFALSE;
+       }
+       
+       if(!fPIDResponse){
+               printf("freidtlog No PIDd\n");
+               // PID is fixed to unknown
+               //return kFALSE;
+       }
+
+       if(fDoAOD && fabs(fAODEvent->GetMagneticField())<1){
+               printf("freidtlog strange Bfield! %f\n", fAODEvent->GetMagneticField());
+               return kFALSE;
+       }
+       else if((!fDoAOD) && fabs(fESDEvent->GetMagneticField())<1){
+               printf("freidtlog strange Bfield! %f\n", fESDEvent->GetMagneticField());
+               return kFALSE;
+       }
+
+       Int_t tmprun = 0;
+       if (fDoAOD) {
+               tmprun = fAODEvent->GetRunNumber();
+       }
+       else {
+               tmprun = fESDEvent->GetRunNumber();
+       }
+
+       if(fRun!=tmprun){
+               fRun = tmprun;
+               AliCDMesonUtils::SPDLoadGeom(fRun);
+
+               // change PID cuts if necessary
+               if ((fRun >= 162933) && (fRun <=  165462)) {
+                       fPIDmode = 1;
+               }
+               else {
+                       fPIDmode = 0;
+               }
+       }
+
+       // get MC event
+       fMCEvent = MCEvent();
+
+       return kTRUE;
+}
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::DoEmptyEventStudy() {
+       // Analyses the observables needed for the gap determination in empty events
+       // for cross checks same observables are stored for A/C/I(B)-triggers as well
+       //
+
+
+       if (!fAnalysisStatus & AliCDMesonBase::kBitEEStudy) {
+               // check whether empty event analysis is activated
+               return;
+       }
+
+       if(!fESDEvent || fDoAOD) { // ensure that we are running on ESDs
+               return;
+       }
+
+       Int_t eventType = AliCDMesonUtils::GetEventType(fESDEvent);
+       TRandom3 rnd(0);
+       if ((eventType != AliCDMesonBase::kBinEventUnknown) &&
+           ((rnd.Rndm() > 0.95) ||
+            (AliCDMesonBase::kBinEventE == eventType))) {
+               Int_t fmdA = 0, fmdC = 0, spdIA = 0, spdIC = 0, spdOA = 0, spdOC = 0,
+                       spdTrklSum = 0, spdTrklCentral = 0, spdTrklForwardA = 0,
+                       spdTrklForwardC = 0;
+
+               Float_t fmdSums[] = { 0., 0., 0., 0., 0. };
+               AliCDMesonUtils::GetMultFMD(fESDEvent, fmdA, fmdC, fmdSums);
+               // Obtain FMD multiplicity using FMDHitCombinations
+               AliCDMesonUtils::GetMultSPD(fESDEvent, spdIA, spdIC, spdOA, spdOC);
+               // Obtain SPD multiplicity using FastOR information
+               AliCDMesonUtils::GetSPDTrackletMult(fESDEvent, spdTrklSum,
+                                                   spdTrklForwardA, spdTrklForwardC,
+                                                   spdTrklCentral);
+               // obtain SPD multiplicity using tracklets via AliMultiplicity
+
+               AliCDMesonBase::FillThnEmptyEvents(fThnEmptyEvents, eventType, fmdA, fmdC,
+                                                  spdIA, spdIC, spdOA, spdOC,
+                                                  spdTrklForwardA, spdTrklForwardC,
+                                                  (Int_t)fmdSums[0], (Int_t)fmdSums[1],
+                                                  (Int_t)fmdSums[2], (Int_t)fmdSums[3],
+                                                  (Int_t)fmdSums[4]);
+       }
+       // the following code is only executed for empty events, this ensures, that
+       // all the histograms contain only information on empty events, while the
+       // THnSparse contains all information on A/C/I events
+       if (AliCDMesonBase::kBinEventE == eventType) {
+               TH1F* fmdSums[] = { fFMDsum1I, fFMDsum2I, fFMDsum2O, fFMDsum3I,
+                                   fFMDsum3O };
+               AliCDMesonUtils::GetGapConfig(fESDEvent, fHitMapSPDinner, fHitMapSPDouter,
+                                             fHitMapSPDtrklt, fHitMapFMDa, fHitMapFMDc,
+                                             (TH1**)fmdSums, fTPCGapDCAaSide,
+                                             fTPCGapDCAcSide);
+               // fill hit maps
+       }
+       return;
+}
+
+
+//------------------------------------------------------------------------------
+Bool_t AliAnalysisTaskCDMeson::DetermineGap()
+{
+       // determines the gap configuration for all gap tagging detectors based on the
+       // data set which is available
+       //
+
+       fCurrentGapCondition = 0x0; // initialize gap condition
+
+       if (fDoAOD) {
+               AliAODHandler* aodHandler =
+                       (AliAODHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+               TTree *aodTree = aodHandler->GetTree();
+               if (aodTree) {
+                       aodTree->SetBranchAddress("gapCondition", &fCurrentGapCondition);
+                       aodTree->GetEvent(Entry()); // seems to be needed! (loads current event)
+               }
+               if (!fCurrentGapCondition) {
+                       fCurrentGapCondition = 0xfffe;
+                       puts("freidtlog - error while gap condition determination using AODs\n");
+                       return kFALSE;
+               }
+               // DEBUGGING STUFF
+               //aodTree->ls();
+               //aodTree->MakeCode();
+               //aodTree->MakeClass();
+               //aodTree->GetEvent(Entry());
+               //aodTree->Show();
+               //TBranch* branch = aodTree->GetBranch("gapCondition");
+               //printf("freidtlog - branch=x%x\n",branch);
+               //if (!fCurrentGapCondition) {
+               //branch->SetAddress(&fCurrentGapCondition);
+                       //}
+               //Int_t entry = Entry();
+               //printf("Entry()=%d\n", entry);
+               //branch->GetEvent(entry);
+               //printf("freidtlog - gapcondition=0x%x\n", fCurrentGapCondition);
+       }
+       else {
+               if (fAnalysisStatus & AliCDMesonBase::kBitReadPreprocessedGap) {
+                       // retrieve preprocessed gap information
+                       AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+                       TTree* esdTree = am->GetInputEventHandler()->GetTree(); // get ESD tree
+                       if (esdTree) {
+                               esdTree->SetBranchAddress("gapCondition", &fCurrentGapCondition);
+                               esdTree->GetEvent(Entry()); // seems to be needed! (loads current event)
+                       }
+               }
+               else {
+                       // determine gaps from ESDs
+                       TH1F* fmdSums[] = {fFMDsum1I, fFMDsum2I, fFMDsum2O, fFMDsum3I, fFMDsum3O};
+                       fCurrentGapCondition = AliCDMesonUtils::GetGapConfig(fESDEvent,
+                                                                            fHitMapSPDinner,
+                                                                            fHitMapSPDouter,
+                                                                            fHitMapSPDtrklt,
+                                                                            fHitMapFMDa,
+                                                                            fHitMapFMDc,
+                                                                            (TH1**)fmdSums,
+                                                                            fTPCGapDCAaSide,
+                                                                            fTPCGapDCAcSide);
+               }
+               if (!fCurrentGapCondition) {
+                       fCurrentGapCondition = 0xfffe;
+                       puts("freidtlog - error while gap condition determination using ESDs\n");
+                       return kFALSE;
+               }
+       }
+       //printf("freidtlog - Event: %ld, gapCondition: 0x%x\n", (Long_t)Entry(),
+       //       fCurrentGapCondition);
+
+
+       // disentagle the contributions to the gap conditions of different "tightness"
+       fGapInformation[kV0] =
+               AliCDMesonBase::GetGapBin("V0", fCurrentGapCondition, kFALSE);
+       fGapInformation[kV0FMD] =
+               AliCDMesonBase::GetGapBin("V0FMD", fCurrentGapCondition, kFALSE);
+       fGapInformation[kV0FMDSPD] =
+               AliCDMesonBase::GetGapBin("V0FMDSPD", fCurrentGapCondition, kFALSE);
+       fGapInformation[kV0FMDSPDTPC] =
+               AliCDMesonBase::GetGapBin("V0FMDSPDTPC", fCurrentGapCondition, kFALSE);
+       fGapInformation[kFMD] =
+               AliCDMesonBase::GetGapBin("FMD",fCurrentGapCondition, kFALSE);
+       fGapInformation[kSPD] =
+               AliCDMesonBase::GetGapBin("SPD",fCurrentGapCondition, kFALSE);
+       fGapInformation[kTPC] =
+               AliCDMesonBase::GetGapBin("TPC",fCurrentGapCondition, kFALSE);
+       fGapInformation[kTPCSPD] =
+               AliCDMesonBase::GetGapBin("TPCSPD",fCurrentGapCondition, kFALSE);
+       fGapInformation[kTPCSPDFMD] =
+               AliCDMesonBase::GetGapBin("TPCSPDFMD",fCurrentGapCondition, kFALSE);
+       fGapInformation[kTPCSPDFMDV0] =
+               AliCDMesonBase::GetGapBin("TPCSPDFMDV0",fCurrentGapCondition, kFALSE);
+       fGapInformation[kSPDFMD] =
+               AliCDMesonBase::GetGapBin("SPDFMD",fCurrentGapCondition, kFALSE);
+       fGapInformation[kSPDFMDV0] =
+               AliCDMesonBase::GetGapBin("SPDFMDV0",fCurrentGapCondition, kFALSE);
+
+       fGapInformationWCent[kV0] = AliCDMesonBase::GetGapBin("V0", fCurrentGapCondition);
+       fGapInformationWCent[kV0FMD] =
+               AliCDMesonBase::GetGapBin("V0FMD", fCurrentGapCondition);
+       fGapInformationWCent[kV0FMDSPD] =
+               AliCDMesonBase::GetGapBin("V0FMDSPD", fCurrentGapCondition);
+       fGapInformationWCent[kV0FMDSPDTPC] =
+               AliCDMesonBase::GetGapBin("V0FMDSPDTPC", fCurrentGapCondition);
+       fGapInformationWCent[kFMD] =
+               AliCDMesonBase::GetGapBin("FMD",fCurrentGapCondition);
+       fGapInformationWCent[kSPD] =
+               AliCDMesonBase::GetGapBin("SPD",fCurrentGapCondition);
+       fGapInformationWCent[kTPC] =
+               AliCDMesonBase::GetGapBin("TPC",fCurrentGapCondition);
+       fGapInformationWCent[kTPCSPD] =
+               AliCDMesonBase::GetGapBin("TPCSPD",fCurrentGapCondition);
+       fGapInformationWCent[kTPCSPDFMD] =
+               AliCDMesonBase::GetGapBin("TPCSPDFMD",fCurrentGapCondition);
+       fGapInformationWCent[kTPCSPDFMDV0] =
+               AliCDMesonBase::GetGapBin("TPCSPDFMDV0",fCurrentGapCondition);
+       fGapInformationWCent[kSPDFMD] =
+               AliCDMesonBase::GetGapBin("SPDFMD",fCurrentGapCondition);
+       fGapInformationWCent[kSPDFMDV0] =
+               AliCDMesonBase::GetGapBin("SPDFMDV0",fCurrentGapCondition);
+
+       return kTRUE;
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::DoMultiplicityStudy(Int_t nMCprimaries)
+{
+       // stores the multiplicity distributions for different gap conditions and
+       // compares the multiplicity distribution for "normal primary tracks (ITSTPC)"
+       // and the multiplicity distribution furthermore taking also soft tracks into
+       // account
+       //
+
+       // retrieve values from the track object
+       Int_t ntrk0 = fTracks->GetTracksBeforeCuts();
+       Int_t nch = fTracks->GetTracks();
+       Int_t ncombined = fTracks->GetCombinedTracks();
+       Int_t nITSpureSA = fTracks->GetITSpureSACount();
+
+       // determine the residual tracks / tracklets
+       fResidualTracks = ntrk0 - ncombined - nITSpureSA;
+       fResidualTracklets = fTracks->GetRemainingTrackletsCentralBarrel();
+
+       // soft track specific stuff
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)) {
+               fMultStudy->Fill((Double_t)nch, (Double_t)ncombined);
+               if (fGapInformation[kV0] == AliCDMesonBase::kBinDG) {
+                       fMultStudyV0dg->Fill((Double_t)nch, (Double_t)ncombined);
+               }
+               if (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG) {
+                       fMultStudyV0FMDdg->Fill((Double_t)nch, (Double_t)ncombined);
+               }
+               if (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+                       fMultStudyV0FMDSPDdg->Fill((Double_t)nch, (Double_t)ncombined);
+               }
+               if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+                       fMultStudyV0FMDSPDTPCdg->Fill((Double_t)nch, (Double_t)ncombined);
+               }
+       }
+
+
+       // multiplicity distributions for different gaps
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitMultPerGapHists)) {
+               fv0ntrk->Fill(ncombined, fGapInformation[kV0]);
+               fv0fmdntrk->Fill(ncombined, fGapInformation[kV0FMD]);
+               fv0fmdspdntrk->Fill(ncombined, fGapInformation[kV0FMDSPD]);
+               fv0fmdspdtpcntrk->Fill(ncombined, fGapInformation[kV0FMDSPDTPC]);
+       }
+
+       // multiplicity removed by cuts for different gaps
+       if (!fAnalysisStatus ||
+           (fAnalysisStatus & AliCDMesonBase::kBitRmMultPerGapHists)) {
+               fv0Rmntrk->Fill(ntrk0-ncombined-nITSpureSA, fGapInformation[kV0]);
+               fv0fmdRmntrk->Fill(ntrk0-ncombined-nITSpureSA, fGapInformation[kV0FMD]);
+               fv0fmdspdRmntrk->Fill(ntrk0-ncombined-nITSpureSA,
+                                     fGapInformation[kV0FMDSPD]);
+               fv0fmdspdtpcRmntrk->Fill(ntrk0-ncombined-nITSpureSA,
+                                        fGapInformation[kV0FMDSPDTPC]);
+       }
+
+
+       // fill maps with eta phi information for QA purposes
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitEtaPhiMaps) ||
+           (fAnalysisStatus & AliCDMesonBase::kBitEtaPhiMapsWithCuts)) {
+               FillEtaPhiMaps();
+       }
+
+       // fill stats and general multiplicity histograms
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               fhntrk->Fill(ncombined); // multiplicity distribution
+
+               if (fGapInformationWCent[kV0] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinv0Gap);
+               }
+               if (fGapInformationWCent[kV0FMD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinv0fmdGap);
+               }
+               if (fGapInformationWCent[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinv0fmdspdGap);
+               }
+               if (fGapInformationWCent[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinv0fmdspdtpcGap);
+               }
+               if (fGapInformationWCent[kTPC] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBintpcGap);
+               }
+               if (fGapInformationWCent[kTPCSPD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBintpcspdGap);
+               }
+               if (fGapInformationWCent[kTPCSPDFMD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBintpcspdfmdGap);
+               }
+               if (fGapInformationWCent[kTPCSPDFMDV0] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBintpcspdfmdv0Gap);
+               }
+               if (fGapInformationWCent[kFMD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinfmdGap);
+               }
+               if (fGapInformationWCent[kSPD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinspdGap);
+               }
+               if (fGapInformationWCent[kSPDFMD] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinspdfmdGap);
+               }
+               if (fGapInformationWCent[kSPDFMDV0] == AliCDMesonBase::kBinDG) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinspdfmdv0Gap);
+               }
+               if (nch == 2) fhStatsFlow->Fill(AliCDMesonBase::kBinTwoTrackEvents);
+               else if (nch == 3) fhStatsFlow->Fill(AliCDMesonBase::kBinThreeTrackEvents);
+       }
+
+       // fill multiplicity response histograms (data vs. MC)
+       if (fAnalysisStatus & AliCDMesonBase::kBitMultResponseMC) {
+               fMultResponseMC->Fill(nMCprimaries, ncombined);
+               if (fGapInformation[kV0] == AliCDMesonBase::kBinDG) {
+                       fMultResponseV0dgMC->Fill(nMCprimaries, ncombined);
+               }
+               if (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG) {
+                       fMultResponseV0FMDdgMC->Fill(nMCprimaries, ncombined);
+               }
+               if (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+                       fMultResponseV0FMDSPDdgMC->Fill(nMCprimaries, ncombined);
+               }
+               if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+                       fMultResponseV0FMDSPDTPCdgMC->Fill(nMCprimaries, ncombined);
+               }
+       }
+
+       // event cleanliness
+       if (fResidualTracks == 0 && fhStatsFlow) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinResidualTracks);
+       }
+       if (fResidualTracklets == 0 && fhStatsFlow) {
+               fhStatsFlow->Fill(AliCDMesonBase::kBinResidualTracklets);
+       }
+
+
+       // multiplicity THnSparse
+       if (!fAnalysisStatus
+           || (fAnalysisStatus & AliCDMesonBase::kBitMultStudy)) {
+               AliCDMesonBase::FillThnMultiplicity(fThnMultiplicity, nch,
+                                                   ncombined-nch, ncombined,
+                                                   fGapInformation[kV0],
+                                                   fGapInformation[kFMD],
+                                                   fGapInformation[kSPD],
+                                                   fGapInformation[kTPC],
+                                                   fResidualTracks, fResidualTracklets,
+                                                   fVtxZ, fVtxDst, fMCprocessType);
+       }
+
+       if (fAnalysisStatus & AliCDMesonBase::kBitFastORmultStudy) {
+               if (fhspdV0dg && fhspdV0FMDdg && fhspdV0FMDSPDdg && fhspdV0FMDSPDTPCdg &&
+                   !fDoAOD) {
+                       Int_t mult = AliCDMesonUtils::GetFastORmultiplicity(fESDEvent);
+
+                       if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+                               fhspdV0FMDSPDTPCdg->Fill(mult);
+                               fhspdV0FMDSPDdg->Fill(mult);
+                               fhspdV0FMDdg->Fill(mult);
+                               fhspdV0dg->Fill(mult);
+                               fhspdAfterCuts->Fill(mult);
+                       }
+                       else if (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+                               fhspdV0FMDSPDdg->Fill(mult);
+                               fhspdV0FMDdg->Fill(mult);
+                               fhspdV0dg->Fill(mult);
+                               fhspdAfterCuts->Fill(mult);
+                       }
+                       else if (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG) {
+                               fhspdV0FMDdg->Fill(mult);
+                               fhspdV0dg->Fill(mult);
+                               fhspdAfterCuts->Fill(mult);
+                       }
+                       else if (fGapInformation[kV0] == AliCDMesonBase::kBinDG) {
+                               fhspdV0dg->Fill(mult);
+                               fhspdAfterCuts->Fill(mult);
+                       }
+                       else {
+                               fhspdAfterCuts->Fill(mult);
+                       }
+               }
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::DoTrackPair(Bool_t soft /* = kFALSE */)
+{
+       //
+       // process a two track pair
+       //
+
+       AliCDMesonUtils::SwapTrack((const AliVTrack**)fTrkPair); // random order
+
+       const Int_t combCh =
+               AliCDMesonUtils::GetCombCh((const AliVTrack**)fTrkPair); // uls or ls?
+
+       const Int_t combPID = DoPID(combCh); // retrieve PID information
+
+       // analyze vertex quality
+       Int_t vtxInRng = 0;
+       Int_t vtxCoincidence = 0;
+
+       vtxInRng = (fabs(fVtxZ) < 4.) ? 1 : 0;
+       // 4cm range, determined from detector geometry
+
+       vtxCoincidence = (fVtxDst < fMaxVtxDst) ? 1 : 0;
+       // vertex coincidence of the track and SPD vertex
+
+       // initialize kinematic values
+       Double_t mass = -999., pt = -999., cts = -999., oa = -999.;
+       const TVector3 p1(fTrkPair[0]->Px(), fTrkPair[0]->Py(), fTrkPair[0]->Pz());
+       const TVector3 p2(fTrkPair[1]->Px(), fTrkPair[1]->Py(), fTrkPair[1]->Pz());
+       const TVector3* momenta[2] = { &p1, &p2 };
+       AliCDMesonUtils::GetMassPtCtsOA(combPID, momenta, mass, pt, cts, oa);
+
+       //= THnMother ================================================================
+       if (!fAnalysisStatus ||
+           ((fAnalysisStatus & AliCDMesonBase::kBitTHnMother)
+            && (!soft || fAnalysisStatus & AliCDMesonBase::kBitSoftTracks))) {
+               TRandom3 rnd(0);
+               if (((fGapInformation[kV0] == AliCDMesonBase::kBinDG)
+                    && (!(fAnalysisStatus & AliCDMesonBase::kBitReduceGapEvents)
+                        || (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG)
+                        || (rnd.Rndm() > (1. - fReducedGapFraction))))
+                   || (rnd.Rndm() > (1. - fNoGapFraction)) || fDoAOD
+                   || (vtxInRng && vtxCoincidence && (fResidualTracks == 0)
+                       && (fResidualTracklets == 0))) {
+                       // reduce amount of no-gap events by a factor of 10 in order to spare
+                       // memory and shrink the final data volume
+
+                       // do switching needed to process either events including soft or without
+                       THnSparseD* thn = (soft) ? fThnMotherSoft : fThnMother; // output variable
+                       Int_t mult = (soft) ? fTracks->GetCombinedTracks() : fTracks->GetTracks();
+                       // correct multiplicity value
+                       Int_t residualTracks = (soft) ?
+                               fResidualTracks : fResidualTracks + fTracks->GetSoftTracks();
+                       // if the analysis is done without soft tracks, they will be added to the
+                       // residualTracks
+
+                       // fill MC histograms
+                       if (fMCEvent && (!fAnalysisStatus ||
+                                        fAnalysisStatus & AliCDMesonBase::kBitMCProcess)) {
+                               if (mult == 2) {
+                                       FillMChists(combCh);
+                               }
+                       }
+
+                       AliCDMesonBase::FillThnMother(thn, (Double_t)mult, (Double_t)combCh,
+                                                     (Double_t)combPID,
+                                                     (Double_t)fGapInformation[kV0],
+                                                     (Double_t)fGapInformation[kFMD],
+                                                     (Double_t)fGapInformation[kSPD],
+                                                     (Double_t)fGapInformation[kTPC],
+                                                     mass, pt, cts, oa, fTrkPair[0]->Pt(),
+                                                     (residualTracks == 0) ? 0. : 1.,
+                                                     (Double_t)vtxInRng,
+                                                     (Double_t)fMCprocessType,
+                                                     (Double_t)vtxCoincidence,
+                                                     (fResidualTracklets == 0) ? 0. : 1.);
+               }
+       }
+
+       //= PWA ======================================================================
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitPWAtree)) {
+               if ((soft && (fTracks->GetCombinedTracks() == 2))
+                   || ((fTracks->GetTracks() == 2)
+                       && !(!fAnalysisStatus
+                            || (fAnalysisStatus & AliCDMesonBase::kBitSoftTracks)))) {
+                       // TODO refine this condition
+                       // TODO add flag for 3 track events and soft tracks and the other quality
+                       // flags
+                       // PWA is only done for two track events
+                       AliCDMesonUtils::GetPWAinfo(combPID, (const AliVTrack**)fTrkPair, fTheta,
+                                                   fPhi, fMass, fMomentum);
+                       TRandom3 rnd(0);
+                       if (((fGapInformation[kV0] == AliCDMesonBase::kBinDG)
+                            && (!(fAnalysisStatus & AliCDMesonBase::kBitReduceGapEvents)
+                                || (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG)
+                                || (rnd.Rndm() > (1. - fReducedGapFraction)))) ||
+                           (rnd.Rndm() > (1. - fNoGapFraction))) {
+                               fPWAtree->Fill();
+                       }
+               }
+       }
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliAnalysisTaskCDMeson::DoPID(Int_t combCh)
+{
+       // determine the PID of the two tracks for full double gap events, store
+       // the results in another histogram than for the remaining events
+       //
+
+       Int_t combPID = 0;
+       if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+               // double-gap events
+               combPID = (combCh == AliCDMesonBase::kBinPM) ? // assignment to uls/ls histo
+                       AliCDMesonUtils::GetCombPID(fPIDResponse, (const AliVTrack**)fTrkPair,
+                                                   fPIDmode, fComb2trkPIDulsDG) :
+                       AliCDMesonUtils::GetCombPID(fPIDResponse, (const AliVTrack**)fTrkPair,
+                                                   fPIDmode, fComb2trkPIDlsDG);
+       }
+       else { // non full double-gap events
+               combPID = (combCh == AliCDMesonBase::kBinPM) ? // assignment to uls/ls histo
+                       AliCDMesonUtils::GetCombPID(fPIDResponse, (const AliVTrack**)fTrkPair,
+                                                   fPIDmode, fComb2trkPIDuls) :
+                       AliCDMesonUtils::GetCombPID(fPIDResponse, (const AliVTrack**)fTrkPair,
+                                                   fPIDmode, fComb2trkPIDls);
+       }
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitStatsFlow)) {
+               if (combPID == AliCDMesonBase::kBinPion ||
+                   combPID == AliCDMesonBase::kBinPionE) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinPionEvents);
+               }
+               else if (combPID == AliCDMesonBase::kBinKaon ||
+                                        combPID == AliCDMesonBase::kBinKaonE) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinKaonEvents);
+               }
+               else if (combPID == AliCDMesonBase::kBinProton ||
+                        combPID == AliCDMesonBase::kBinProtonE) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinProtonEvents);
+               }
+               else if (combPID == AliCDMesonBase::kBinElectron ||
+                        combPID == AliCDMesonBase::kBinElectronE) {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinElectronEvents);
+               }
+               else {
+                       fhStatsFlow->Fill(AliCDMesonBase::kBinUnknownPIDEvents);
+               }
+       }
+
+       return combPID;
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::DetermineMCprocessType()
+{
+       //
+       // retrieves the MC process type from the AliGenEventHeader and classifies
+       // them
+       //
+
+       // get MC information
+       fMCprocess = -1; //detailed MC sub process information
+       fMCprocessType = AliCDMesonBase::kBinND; // ND is default, also for data
+
+       if (fMCEvent) {
+               AliGenEventHeader* header = fMCEvent->GenEventHeader();
+               if (header) {
+                       // Pythia6
+                       if (TString(header->IsA()->GetName()) == "AliGenPythiaEventHeader") {
+                               fMCprocess = ((AliGenPythiaEventHeader*)header)->ProcessType();
+                               switch(fMCprocess) {
+                               case 92:
+                               case 93:
+                               case 103:
+                               case 104: fMCprocessType = AliCDMesonBase::kBinSD; break;
+                               case 94:
+                               case 105: fMCprocessType = AliCDMesonBase::kBinDD; break;
+                               default: fMCprocessType = AliCDMesonBase::kBinND; break;
+                               }
+                       }
+                       // Phojet
+                       else if (TString(header->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
+                               fMCprocess = ((AliGenDPMjetEventHeader*)header)->ProcessType();
+                               switch(fMCprocess) {
+                               case 5:
+                               case 6: fMCprocessType = AliCDMesonBase::kBinSD; break;
+                               case 7: fMCprocessType = AliCDMesonBase::kBinDD; break;
+                               case 4: fMCprocessType = AliCDMesonBase::kBinCD; break;
+                               default: fMCprocessType = AliCDMesonBase::kBinND; break;
+                               }
+                       }
+               }
+       }
+}
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::FillMChists(Int_t combCh)
+{
+       // fill the MC histograms for different gap conditions. Only real two track
+       // events are taken into account
+       //
+
+       if (fMCprocess > 0) { // MC process
+               for (Int_t iGap = kV0; iGap < kMax; ++iGap) {
+                       if (fGapInformation[iGap] == AliCDMesonBase::kBinDG) {
+                               if (combCh == AliCDMesonBase::kBinPM) {
+                                       fMCProcessUls->Fill(fMCprocess, iGap);
+                               }
+                               else if (combCh == AliCDMesonBase::kBinPPMM) {
+                                       fMCProcessLs->Fill(fMCprocess, iGap);
+                               }
+                       }
+               }
+       }
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliAnalysisTaskCDMeson::DoMCTruth()
+{
+       //
+       // analyses the MC truth and does a gap selection based on the eta values of
+       // the primaries
+       //
+
+       if (!fMCEvent) return -1;
+       AliStack* stack = fMCEvent->Stack();
+       if (!stack) return -1;
+
+       //= Multiplicity =============================================================
+       // determine number of charged physical primaries on the stack
+       Int_t nPhysicalPrimaries = 0;
+       Int_t nPrimaries = stack->GetNprimary();
+       for (Int_t iTracks = 0; iTracks < nPrimaries; ++iTracks) {
+               TParticle* part = stack->Particle(iTracks);
+               if (stack->IsPhysicalPrimary(iTracks) && (part->GetPDG()->Charge() != 0.) &&
+                   (part->Eta() < 0.9) && (part->Eta() > -0.9)) { // TODO add parameter
+                       ++nPhysicalPrimaries;
+               }
+       }
+
+
+       //= Track Gap && Multiplicity in Region Bins =================================
+       Int_t gapA = 0;
+       Int_t gapAv0 = 0;
+       Int_t gapAv0fmd = 0;
+       Int_t gapC = 0;
+       Int_t gapCv0 = 0;
+       Int_t gapCv0fmd = 0;
+       Int_t central = 0;
+       for (Int_t iTracks = 0; iTracks <  nPrimaries; ++iTracks) {
+               TParticle* part = (TParticle*)stack->Particle(iTracks);
+               if (part && stack->IsPhysicalPrimary(iTracks) &&
+                   (part->GetPDG()->Charge() != 0.)) {
+                       if (part->Eta() > -0.9 && part->Eta() < 0.9) central++;
+                       if (part->Eta() > 0.9 && part->Eta() < 5.1) gapA++;
+                       if (part->Eta() > 2.8 && part->Eta() < 5.1) gapAv0++;
+                       if (part->Eta() > 1.7 && part->Eta() < 5.1) gapAv0fmd++;
+                       if (part->Eta() < -0.9 && part->Eta() > -3.7) gapC++;
+                       if (part->Eta() < -1.9 && part->Eta() > -3.7) gapCv0++;
+                       if (part->Eta() < -1.9 && part->Eta() > -3.7) gapCv0fmd++;
+               }
+       }
+
+       if (fMultRegionsMC) {
+               // multiplicity distribution separated in A-side, central barrel and C-side
+               fMultRegionsMC->Fill(gapA, 0);
+               fMultRegionsMC->Fill(central, 1);
+               fMultRegionsMC->Fill(gapC, 2);
+               fMultRegionsMC->Fill(gapA+gapC, 3);
+       }
+
+       // WARNING: gap response determination based on primaries only, method too
+       // simple to obtain reliable results; at least decays should be taken into
+       // account properly
+       if (fGapResponseMCv0Dg) {
+               Bool_t gapFromMC = (gapAv0 == 0) && (gapCv0 == 0);
+               Bool_t gapFromData = (fGapInformation[kV0] == AliCDMesonBase::kBinDG);
+
+               if (gapFromMC && !gapFromData) {
+                       fGapResponseMCv0Dg->Fill(0);
+               }
+               else if (gapFromMC && gapFromData) {
+                       fGapResponseMCv0Dg->Fill(1);
+               }
+               else if (!gapFromMC && gapFromData) {
+                       fGapResponseMCv0Dg->Fill(2);
+               }
+       }
+       if (fGapResponseMCv0fmdDg) {
+               Bool_t gapFromMC = (gapAv0fmd == 0) && (gapCv0fmd == 0);
+               Bool_t gapFromData = (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG);
+
+               if (gapFromMC && !gapFromData) {
+                       fGapResponseMCv0fmdDg->Fill(0);
+               }
+               else if (gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdDg->Fill(1);
+               }
+               else if (!gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdDg->Fill(2);
+               }
+       }
+       if (fGapResponseMCv0fmdspdDg) {
+               Bool_t gapFromMC = (gapA == 0) && (gapC == 0);
+               Bool_t gapFromData = (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG);
+
+               if (gapFromMC && !gapFromData) {
+                       fGapResponseMCv0fmdspdDg->Fill(0);
+               }
+               else if (gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdspdDg->Fill(1);
+               }
+               else if (!gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdspdDg->Fill(2);
+               }
+       }
+       if (fGapResponseMCv0fmdspdtpcDg) {
+               Bool_t gapFromMC = (gapA == 0) && (gapC == 0);
+               Bool_t gapFromData =
+                       (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG);
+
+               if (gapFromMC && !gapFromData) {
+                       fGapResponseMCv0fmdspdtpcDg->Fill(0);
+               }
+               else if (gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdspdtpcDg->Fill(1);
+               }
+               else if (!gapFromMC && gapFromData) {
+                       fGapResponseMCv0fmdspdtpcDg->Fill(2);
+               }
+       }
+
+       if ((nPhysicalPrimaries < 2 ) || (nPhysicalPrimaries >= 4 ))
+               return nPhysicalPrimaries;
+
+       Int_t gapCond = AliCDMesonBase::kBinNG;
+       if ((gapA == 0) && (gapC == 0) && (central > 0))
+               gapCond = AliCDMesonBase::kBinDG;
+       else if ((gapA == 0) && (gapC > 0)) gapCond = AliCDMesonBase::kBinGA;
+       else if ((gapA > 0) && (gapC == 0)) gapCond = AliCDMesonBase::kBinGC;
+       else if ((gapA > 0) && (gapC > 0)) gapCond = AliCDMesonBase::kBinNG;
+
+       //= Two Track ================================================================
+       for (Int_t iTracks = 0; iTracks < stack->GetNprimary(); ++iTracks) {
+               TParticle* part1 =0x0;
+               part1 = (TParticle*)stack->Particle(iTracks);
+               if (part1) {
+                       if (stack->IsPhysicalPrimary(iTracks) &&
+                           (part1->GetPDG()->Charge() != 0.) && (part1->Eta() < 0.9) &&
+                           (part1->Eta() > -0.9)) {
+                               for (Int_t jTracks = iTracks; jTracks < stack->GetNprimary(); ++jTracks) {
+                                       TParticle* part2 = 0x0;
+                                       part2 = (TParticle*)stack->Particle(jTracks);
+                                       if (part2) {
+                                               if (stack->IsPhysicalPrimary(jTracks) &&
+                                                   (part2->GetPDG()->Charge() != 0.) && (part2->Eta() < 0.9) &&
+                                                   (part2->Eta() > -0.9)) {
+                                                       TParticle* particles[2] = { part1, part2 };
+                                                       DoMCTrackPair(particles, gapCond, nPhysicalPrimaries);
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+       return nPhysicalPrimaries;
+}
+
+
+//------------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::DoMCTrackPair(TParticle* particles[],
+                                           Int_t gapCond, Int_t multiplicity)
+{
+       //
+       // analyses a two track pair from MC truth
+       //
+
+       // swap tracks randomly
+       TRandom3 rnd(0);
+       if (rnd.Rndm() >= 0.5) {
+               TParticle* t = particles[1];
+               particles[1] = particles[0];
+               particles[0] = t;
+       }
+
+       // determine whether the pair is like sign or unlike sign
+       const Int_t combCh =
+               (particles[0]->GetPDG()->Charge()*particles[1]->GetPDG()->Charge() > 0.) ?
+               AliCDMesonBase::kBinPPMM : AliCDMesonBase::kBinPM;
+
+       // determine PID
+       Int_t combPID = 0;
+       if (gapCond == AliCDMesonBase::kBinDG) {
+               combPID = (combCh == AliCDMesonBase::kBinPM) ? // assignment to uls/ls histo
+                       AliCDMesonUtils::GetCombPID((const TParticle**)particles,
+                                                   fComb2trkPIDulsDGmc) :
+                       AliCDMesonUtils::GetCombPID((const TParticle**)particles,
+                                                   fComb2trkPIDlsDGmc);
+       }
+       else {
+               combPID = (combCh == AliCDMesonBase::kBinPM) ? // assignment to uls/ls histo
+                       AliCDMesonUtils::GetCombPID((const TParticle**)particles,
+                                                   fComb2trkPIDulsMC) :
+                       AliCDMesonUtils::GetCombPID((const TParticle**)particles,
+                                                   fComb2trkPIDlsMC);
+       }
+
+       // vertex quality is always nice for MC truth ;)
+       Int_t vtxInRng = 1;
+       Int_t vtxCoincidence = 1;
+
+       // initialize kinematic values
+       Double_t mass = -999., pt = -999., cts = -999., oa = -999.;
+       const TVector3 p1(particles[0]->Px(), particles[0]->Py(), particles[0]->Pz());
+       const TVector3 p2(particles[1]->Px(), particles[1]->Py(), particles[1]->Pz());
+       const TVector3* momenta[2] = { &p1, &p2 };
+       AliCDMesonUtils::GetMassPtCtsOA(combPID, momenta, mass, pt, cts, oa);
+
+       //= THnMother ================================================================
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitTHnMother)) {
+               if ((gapCond != AliCDMesonBase::kBinNG) ||
+                   (rnd.Rndm() > (1. - fNoGapFraction)) || fDoAOD) {
+
+                       AliCDMesonBase::FillThnMother(fThnMC, multiplicity, (Double_t)combCh,
+                                                     (Double_t)combPID, (Double_t)gapCond,
+                                                     (Double_t)gapCond, (Double_t)gapCond,
+                                                     (Double_t)gapCond, mass, pt, cts, oa,
+                                                     particles[0]->Pt(), 0., (Double_t)vtxInRng,
+                                                     (Double_t)fMCprocessType,
+                                                     (Double_t)vtxCoincidence, 1.);
+               }
+       }
+
+       //= PWA ======================================================================
+       // TODO ? should this be implemented???
+}
+
+
+//--------------------------------------------------------------------------
+void AliAnalysisTaskCDMeson::AnalyzeVtx()
+{
+       // calculates the distance between the vertex obtain from tracks and the
+       // vertex obtain from spd tracklets
+       // stores the z position of the primary vertex from tracks
+
+       fVtxDst = 0.; // reset distance
+
+       // retrieve the pointers of the current primary vertices
+       AliVVertex* trackVtx = (fDoAOD) ?
+               (AliVVertex*)fAODEvent->GetPrimaryVertex() :
+               (AliVVertex*)fESDEvent->GetPrimaryVertexTracks();
+       AliVVertex* spdVtx = (fDoAOD) ?
+               (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() :
+               (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
+
+       fVtxZ = trackVtx->GetZ(); // store the vertex z position
+
+       if (fDoAOD && (trackVtx == spdVtx)) { // no primary track vertex in the AOD
+               fVtxDst = -5; // set arbitrary distance (counted in underflow bin!)
+       }
+       else { // do proper calculation of the geometrical distance
+               fVtxDst += (trackVtx->GetX()-spdVtx->GetX())
+                       * (trackVtx->GetX()-spdVtx->GetX());
+               fVtxDst += (trackVtx->GetY()-spdVtx->GetY())
+                       * (trackVtx->GetY()-spdVtx->GetY());
+               fVtxDst += (trackVtx->GetZ()-spdVtx->GetZ())
+                       * (trackVtx->GetZ()-spdVtx->GetZ());
+               fVtxDst = TMath::Sqrt(fVtxDst);
+       }
+
+       if (!fAnalysisStatus || (fAnalysisStatus & AliCDMesonBase::kBitVtxStudies)) {
+               fPriVtxDst->Fill(fVtxDst);
+               if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBase::kBinDG) {
+                       fPriVtxDstV0FMDSPDTPCdg->Fill(fVtxDst);
+                       fPriVtxDstV0FMDSPDdg->Fill(fVtxDst);
+                       fPriVtxDstV0FMDdg->Fill(fVtxDst);
+                       fPriVtxDstV0dg->Fill(fVtxDst);
+               }
+               else if (fGapInformation[kV0FMDSPD] == AliCDMesonBase::kBinDG) {
+                       fPriVtxDstV0FMDSPDdg->Fill(fVtxDst);
+                       fPriVtxDstV0FMDdg->Fill(fVtxDst);
+                       fPriVtxDstV0dg->Fill(fVtxDst);
+               }
+               else if (fGapInformation[kV0FMD] == AliCDMesonBase::kBinDG) {
+                       fPriVtxDstV0FMDdg->Fill(fVtxDst);
+                       fPriVtxDstV0dg->Fill(fVtxDst);
+               }
+               else if (fGapInformation[kV0] == AliCDMesonBase::kBinDG) {
+                       fPriVtxDstV0dg->Fill(fVtxDst);
+               }
+       }
+}
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.h b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDMeson.h
new file mode 100644 (file)
index 0000000..bd5bc2f
--- /dev/null
@@ -0,0 +1,281 @@
+/*************************************************************************
+* Copyright(c) 1998-1999, 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.                  *
+**************************************************************************/
+//
+// Select events according to gap conditions, analyze two track events in pp
+// collisions
+//
+// Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#ifndef ALIANALYSISTASKCDMESON_H
+#define ALIANALYSISTASKCDMESON_H
+
+#ifndef ALIANALYSISTASK_H
+#include "AliAnalysisTaskSE.h"
+#endif
+
+#include "THnSparse.h" // forward declaration is not possible
+
+class AliESDEvent;
+class AliVTrack;
+class AliPIDResponse;
+class AliPhysicsSelection;
+
+class TH1I;
+class TH1F;
+class TH2I;
+class TH2F;
+class TH2D;
+class TList;
+class TParticle;
+class THnSparse;
+class TObjString;
+class TObjArray;
+
+class AliCDMesonTracks;
+
+class AliAnalysisTaskCDMeson : public AliAnalysisTaskSE
+{
+public:
+       AliAnalysisTaskCDMeson(const char* name, Long_t state = 0x0);
+       AliAnalysisTaskCDMeson();
+       virtual ~AliAnalysisTaskCDMeson();
+
+       virtual void UserCreateOutputObjects();
+       virtual void UserExec(Option_t *);
+
+private:
+       enum { //
+               kV0 = 0,
+               kFMD,
+               kSPD,
+               kTPC,
+               kV0FMD,
+               kV0FMDSPD,
+               kV0FMDSPDTPC,
+               kTPCSPD,
+               kTPCSPDFMD,
+               kTPCSPDFMDV0,
+               kSPDFMD,
+               kSPDFMDV0,
+               kMax
+       };
+
+       AliAnalysisTaskCDMeson(const AliAnalysisTaskCDMeson  &p);
+       AliAnalysisTaskCDMeson& operator=(const AliAnalysisTaskCDMeson  &p);
+
+       void FillEtaPhiMaps() const; // controls which event contributes to which map
+
+
+       // functions called by the UserExec(...), not to be called elsewhere!
+       //-------------------------------------------------------------------
+       Bool_t CheckInput();
+       void PostOutputs(); // cares about posting the output before exiting UserExec,
+       // WARNING: PostOutputs should only be used directly in the UserExec!!
+       void DoEmptyEventStudy(); // does the empty event study
+       Bool_t DetermineGap(); // determines the gap of all available detectors
+       void AnalyzeVtx(); // calcs the distance of the pri. vertex from tracks and SPD
+       void DoMultiplicityStudy(Int_t nMCprimaries);
+       // compares ITSTPC and ITSTPC+soft distribution and the multiplicity response
+       void DoTrackPair(Bool_t soft = kFALSE);
+       // analyses track pairs and prepares THnMother(Soft) input
+       Int_t DoMCTruth(); // analyses the MCtruth for corrections, returns #Primaries
+       void DoMCTrackPair(TParticle* particles[], Int_t gapCond, Int_t multiplicity);
+       // analyzes track pairs in MC truth
+       void DetermineMCprocessType(); // determines the MC process ID
+
+       // functions called by the DoTrackPair(...), not to be called elsewhere!
+       // ---------------------------------------------------------------------
+       Int_t DoPID(Int_t combCh); //retrieves the PID information for the current two
+       void FillMChists(Int_t combCh); //fills the MC histograms for two track events
+
+       // analysis task status
+       //---------------------
+       Bool_t fDoAOD; // true for running on AODs
+       Long_t fAnalysisStatus; // stores the status bits used to distinguish
+       // which histograms should be generated
+       TObjString *fAnalysisStatusString; // to be stored in the output list
+       Double_t fNoGapFraction; // fraction of no-gap events to be processed only
+                                // used in case of ESD input
+       Double_t fReducedGapFraction; // fraction of events stored altough they
+                                     // show only the requested and not the refined
+                                     // gap
+       Double_t fMaxVtxDst; // maximum distance of the track and SPD vertex
+
+       // event information
+       //------------------
+       AliESDEvent *fESDEvent; // esd event object
+       AliAODEvent *fAODEvent; // esd event object
+       AliPIDResponse *fPIDResponse; // esd pid object
+       AliPhysicsSelection *fPhysicsSelection; // physics selection object
+       AliCDMesonTracks *fTracks; // object taking care about the track cuts
+       Double_t fVtxDst; // distance of the primary vertex from tracks and from SPD
+       Double_t fVtxZ; // z-position of the primary vertex from tracks
+       Int_t fResidualTracks; // tracks rejected by cuts within the event
+       Int_t fResidualTracklets; // SPD tracklets not assigned to tracks
+       Int_t fMCprocessType; // MC process type, 0 for data
+       Int_t fMCprocess; // detailed MC sub process information
+
+       // information about the trackpair which is currently processed
+       AliVTrack* fTrkPair[2]; // track objects
+
+       Int_t fRun; // number of the run which is about to be processed
+       Int_t fPIDmode; // selects set of PID cuts, 0 for 3sigma standard cuts,
+       //1 for LHC11f
+
+       // information stored for the PWA (addresses needed for the tree)
+       Float_t fTheta; // theta angle in the helicity frame
+       Float_t fPhi; // phi angle in the helicity frame
+       Float_t fMass; // mass of the two track system
+       Float_t fMomentum[3]; // momentum of the two track system
+       Int_t fCurrentGapCondition; // gap condition of the current event
+       Int_t fGapInformationWCent[kMax]; // gap condition for different detectors
+                                         // individually and their combinations
+                                         // requiring central act. from SPD FastOR
+       Int_t fGapInformation[kMax]; // same as above, without the requirement
+                                    // on central activity
+
+       // output objects
+       //---------------
+       THnSparseI *fGapRun; //contains the gap information of the different detectors
+       // per run (x-axis: run number, y-axis gap-state (bit encoded)
+
+       TList *fHist; // output list (contains all histograms)
+       THnSparseD *fThnMother; // THnSparse for mother pt and mass
+       THnSparseD *fThnMotherSoft; // same as above, using soft tracks as well
+       THnSparseD *fThnMultiplicity; // ThnSparse for multiplicity studies
+       THnSparseD *fThnMC; // THnSparse for MC mother pt and mass
+
+       THnSparseD *fThnEmptyEvents; // THnSparse used to analyse empty events
+       TTree *fPWAtree; // Tree containing information needed for PWA
+
+       // Multiplicity distributions for the different gap conditions
+       TH2D *fv0ntrk; //v0bit vs. nch
+       TH2D *fv0fmdntrk; //v0fmdbit vs. nch
+       TH2D *fv0fmdspdntrk; //v0fmdspdbit vs. nch
+       TH2D *fv0fmdspdtpcntrk; //v0fmdspdtpcbit vs. nch
+
+       // How many tracks are removed by the current cuts
+       TH2F *fv0Rmntrk; // V0 gap events vs. of removed tracks
+       TH2F *fv0fmdRmntrk; // V0-FMD gap events vs. of removed tracks
+       TH2F *fv0fmdspdRmntrk; // V0-FMD-SPD gap events vs. of removed tracks
+       TH2F *fv0fmdspdtpcRmntrk; // V0-FMD-SPD-TPC gap events vs. of removed tracks
+
+       TH2F *fMultStudy; // multiplicity after standard cuts and soft cuts vs. mult
+       TH2F *fMultStudyV0dg; // same for with events V0 Double-Gap
+       TH2F *fMultStudyV0FMDdg; // same for events with V0-FMD DG
+       TH2F *fMultStudyV0FMDSPDdg; // same for events with V0-FMD-SPD DG
+       TH2F *fMultStudyV0FMDSPDTPCdg; // same for events with V0-FMD-SPD-TPC DG
+
+       TH2F *fMultResponseMC; // multiplicity response MC to reconstruction
+       TH2F *fMultResponseV0dgMC; // same for events with V0 double gap
+       TH2F *fMultResponseV0FMDdgMC; // with V0-FMD DG
+       TH2F *fMultResponseV0FMDSPDdgMC; // with V0-FMD-SPD DG
+       TH2F *fMultResponseV0FMDSPDTPCdgMC; // with V0-FMD-SPD-TPC DG
+
+       TH2F *fMultRegionsMC; // multiplicity in the 3 regions of interest
+
+       TH1I *fhspdV0dg; // FastOR multiplicity in SPD FastOR with V0 double gap
+       TH1I *fhspdV0FMDdg; // with V0-FMD double gap
+       TH1I *fhspdV0FMDSPDdg; // with V0-FMD-SPD double gap
+       TH1I *fhspdV0FMDSPDTPCdg; // with V0-FMD-SPD double gap
+       TH1I *fhspdAfterCuts; // fast OR multiplicity after the event (vertex) cuts
+
+       TH1I *fGapResponseMCv0Dg; // gap status from MC and reconstruction for V0
+       TH1I *fGapResponseMCv0fmdDg; // for V0-FMD
+       TH1I *fGapResponseMCv0fmdspdDg; // for V0-FMD-SPD
+       TH1I *fGapResponseMCv0fmdspdtpcDg; // for V0-FMD-SPD-TPC
+
+       // Statistics flow diagrams
+       TH1I *fhspd; // distribution of fastOR-multiplicity
+       //SPD fastOR offline esdEvent->GetMultiplicity->TestFastOrFiredChips()
+       TH2I *fhfo; // histogram containing the SPD FastOR information ??
+       TH1I *fhpriVtxDist; // primary vertex distribution
+       TH1I *fhpriVtxPos; // position of the primary vertex (0 within range, 1 else)
+       TH1I *fhpriv; //primary vertex cut effect
+       TH1I *fhntrk; //n-trk-after-cut effect w/o requiring a gap-condition
+       TH1F *fhStatsFlow; // stepwise statistics flow
+       TH1F *fhFOchans; // hits in the different fastOR channels
+
+       TObjArray* fVZEROhists; // histograms for the VZERO trigger study
+
+       // eta-phi maps containing the ESDtracks
+       TH2F *fv0Map; // V0 gap events
+       TH2F *fv0fmdMap; // V0-FMD gap events
+       TH2F *fv0fmdspdMap; // V0-FMD-SPD gap events
+       TH2F *fv0fmdspdtpcMap; // V0-FMD-SPD-TPC gap events
+
+       // eta-phi maps containing the ESDtracks after the track cut
+       TH2F *fv0MapCutted; // V0 gap events
+       TH2F *fv0fmdMapCutted; // V0-FMD gap events
+       TH2F *fv0fmdspdMapCutted; // V0-FMD-SPD gap events
+       TH2F *fv0fmdspdtpcMapCutted; // V0-FMD-SPD-TPC gap events
+
+       TH2F *fHitMapSPDinner; // eta-phi map containing SPD FastOR hits (inner layer)
+       TH2F *fHitMapSPDouter; // eta-phi map containing SPD FastOR hits (outer layer)
+       TH2F *fHitMapSPDtrklt; // eta-phi map containing SPD hits based on tracklets
+       TH2F *fHitMapFMDa; // eta-phi-multiplicity map containing FMD hits at A side
+       TH2F *fHitMapFMDc; // eta-phi-multiplicity map containing FMD hits at C side
+
+       // Summed FMD
+       TH1F *fFMDsum1I; // FMD 1
+       TH1F *fFMDsum2I; // FMD 2 inner ring
+       TH1F *fFMDsum2O; // FMD 2 outer ring
+       TH1F *fFMDsum3I; // FMD 3 inner ring
+       TH1F *fFMDsum3O; // FMD 3 outer ring
+
+       // Vertex
+       TH1I *fPriVtxX; // x distribution of the primary vertex
+       TH1I *fPriVtxY; // y distribution of the primary vertex
+       TH1I *fPriVtxZ; // z distribution of the primary vertex
+       TH1I *fPriVtxDst; // distance distribution of the pri. vtx from SPD and tracks
+       TH1I *fPriVtxDstV0dg; // with V0 DG
+       TH1I *fPriVtxDstV0FMDdg; // with V0 FMD DG
+       TH1I *fPriVtxDstV0FMDSPDdg; // with V0 FMD SPD DG
+       TH1I *fPriVtxDstV0FMDSPDTPCdg; //  with V0 FMD SPD TPC DG
+
+       // TPC Gap
+       TH2F *fTPCGapDCAaSide; // dca distribution for tracks destroying the TPC gap
+       TH2F *fTPCGapDCAcSide; // for V0-FMD-SPD gap events on A and C side
+
+       // PID
+       // All events but double gap events
+       TH2F *fComb2trkPIDuls; // PID matrix on unlike-sign two particles events
+       TH2F *fComb2trkPIDls; // PID matrix on like-sign about two particles
+       // Full-gap events (V0-FMD-SPD-TPC)
+       TH2F *fComb2trkPIDulsDG; // PID matrix on unlike-sign two particles events
+       TH2F *fComb2trkPIDlsDG; // PID matrix on like-sign about two particles
+       // 'PID' in MC
+       // All events but double gap events
+       TH2F *fComb2trkPIDulsMC; // PID matrix on unlike-sign two particles events
+       TH2F *fComb2trkPIDlsMC; // PID matrix on like-sign about two particles
+       // Full-gap events (V0-FMD-SPD-TPC)
+       TH2F *fComb2trkPIDulsDGmc; // PID matrix on unlike-sign two particles events
+       TH2F *fComb2trkPIDlsDGmc; // PID matrix on like-sign about two particles
+
+       // MC Process
+       TH2F *fMCProcessUls; // Process fulfilling a gap condition with two
+       // unlike-sign particles
+       TH2F *fMCProcessLs; // same for like-sign particles
+
+       TH1F *fMAllTrackMass; // histogram for the invariant mass dist. of all tracks
+
+       ClassDef(AliAnalysisTaskCDMeson, 1);
+};
+
+
+#endif
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.cxx b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.cxx
new file mode 100644 (file)
index 0000000..06d4c10
--- /dev/null
@@ -0,0 +1,528 @@
+/*************************************************************************
+* Copyright(c) 1998-1999, 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$
+//
+// Task to skim ESD files, basic idea found at PWGGA/EMCALTasks as
+// AliEsdSkimTask. This version is modified to produce a exact copy of the ESD
+// stream containing only events which are of interest for the central
+// diffractive analysis (PWGUD) and a control sample. In order to be able to
+// process MC data as well, the usual concept using an output slot is for
+// output data is violated and an output file is directly created without
+// a corresponding output slot.
+//
+// Author:
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+// ROOT headers
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TH1.h>
+#include <TAxis.h>
+#include <TRandom3.h>
+//#include <TKey.h>
+//#include <TROOT.h>
+
+// AliRoot headers
+#include "AliAnalysisManager.h"
+#include "AliCentrality.h"
+#include "AliESDEvent.h"
+#include "AliESDFMD.h"
+#include "AliEventplane.h"
+#include "AliInputEventHandler.h"
+#include "AliMultiplicity.h"
+#include "AliVEventHandler.h"
+#include "AliMCEventHandler.h"
+//#include "AliMCEvent.h"
+//#include "AliStack.h"
+//#include "AliGenEventHeader.h"
+//#include "AliRunLoader.h"
+
+// own classes
+#include "AliCDMesonBase.h"
+#include "AliCDMesonUtils.h"
+
+// own header
+#include "AliAnalysisTaskCDskimESD.h"
+
+ClassImp(AliAnalysisTaskCDskimESD)
+
+//______________________________________________________________________________
+AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD(const char *opt,
+                                                   Bool_t reduceGapEvents)
+       : AliAnalysisTaskSE(opt)
+       , fInputTree(0x0)
+       , fInputTreeCopy(0x0)
+       , fOutputTree(0x0)
+       , fDoMC(kFALSE)
+       , fWrkDir(0x0)
+       , fKinematicsFile(0x0)
+       , fTrackRefsFile(0x0)
+       , fTEinput(0x0)
+       , fTE(0x0)
+       , fgaliceFile(0x0)
+       , fDoGapCond(kTRUE)
+       , fCurrentGapCondition(0)
+       , fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
+                               // only events containing a V0 double gap or
+       , fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
+       , fReduceGapEvents(reduceGapEvents)
+       , fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
+                              AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
+       , fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
+       , fStatsFlow(0x0)
+       , fSkimmingList(0x0)
+       , fFileName()
+       , fEventNumberInFile(-1)
+       , fRunNumber(-1)
+       , fEventTime(0)
+{
+       //
+       // Constructor.
+       //
+
+       if (!opt)
+               return;
+
+       //fBranchNames = "ESD:AliESDHeader.,AliESDRun."; // TODO don't we need all?
+
+       // check whether we are running on MC data
+       AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+       AliMCEventHandler* mcH =
+               dynamic_cast<AliMCEventHandler*>(am->GetMCtruthEventHandler());
+       if(mcH) {
+               fDoMC = kTRUE;
+
+               fWrkDir = gDirectory; // save the old directory
+               fKinematicsFile = new TFile("Kinematics.root", "recreate");
+               fKinematicsFile->Close();
+               delete fKinematicsFile;
+               fKinematicsFile = 0x0;
+               gDirectory = fWrkDir; // restore the old directory
+
+               fTrackRefsFile = new TFile("TrackRefs.root", "recreate");
+               fTrackRefsFile->Close();
+               delete fTrackRefsFile;
+               fTrackRefsFile = 0x0;
+               gDirectory = fWrkDir; // restore the old directory
+       }
+
+       DefineOutput(1, TTree::Class());
+       DefineOutput(2, TH1::Class());
+       DefineOutput(3, TTree::Class());
+       if (fDoMC) DefineOutput(4, TTree::Class());
+}
+
+
+//______________________________________________________________________________
+AliAnalysisTaskCDskimESD::AliAnalysisTaskCDskimESD()
+       : AliAnalysisTaskSE()
+       , fInputTree(0x0)
+       , fInputTreeCopy(0x0)
+       , fOutputTree(0x0)
+       , fDoMC(kFALSE)
+       , fWrkDir(0x0)
+       , fKinematicsFile(0x0)
+       , fTrackRefsFile(0x0)
+       , fTEinput(0x0)
+       , fTE(0x0)
+       , fgaliceFile(0x0)
+       , fDoGapCond(kTRUE)
+       , fCurrentGapCondition(0)
+       , fRequestedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C)
+                               // only events containing a V0 double gap or
+       , fNoGapFraction(0.006) // 0.6% of the minimum bias events are stored
+       , fReduceGapEvents(kFALSE)
+       , fRefinedGapCondition(AliCDMesonBase::kBitV0A + AliCDMesonBase::kBitV0C +
+                              AliCDMesonBase::kBitFMDA + AliCDMesonBase::kBitFMDC)
+       , fReducedGapFraction(0.02) // 2% of the events with a loose gap are stored
+       , fStatsFlow(0x0)
+       , fSkimmingList(0x0)
+       , fFileName()
+       , fEventNumberInFile(-1)
+       , fRunNumber(-1)
+       , fEventTime(0)
+{
+       //
+       // Default Constructor
+       //
+}
+
+
+//______________________________________________________________________________
+AliAnalysisTaskCDskimESD::~AliAnalysisTaskCDskimESD()
+{
+       //
+       // deconstructor
+       //
+
+       /* // this delete lead to chrases caused by the garbage collection ...
+       if (fOutputTree) {
+               delete fOutputTree;
+               fOutputTree = 0x0;
+               }
+       if (fInputTreeCopy) {
+               delete fInputTreeCopy;
+               fInputTreeCopy = 0x0;
+       }
+       if (fStatsFlow) {
+               delete fStatsFlow;
+               fStatsFlow = 0x0;
+       } */
+       if (fDoMC && fKinematicsFile) {
+               fKinematicsFile->Close();
+               delete fKinematicsFile;
+               fKinematicsFile = 0x0;
+       }
+       if (fDoMC && fTrackRefsFile) {
+               fTrackRefsFile->Close();
+               delete fTrackRefsFile;
+               fTrackRefsFile = 0x0;
+       }
+}
+
+
+//______________________________________________________________________________
+void AliAnalysisTaskCDskimESD::UserExec(Option_t */*opt*/)
+{
+       //
+  // Process event.
+       //
+
+       fStatsFlow->Fill(0.); // UserExec(...) <= bin name
+
+       if (fDoMC) {
+               // replace dummy tree by the correct tree
+               AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+               if(!fTEinput) {
+                       fTEinput =
+                               ((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
+                       if (fTEinput) {
+                               delete fTE;
+                               fTE = fTEinput->CloneTree(0);
+                               fTE->SetDirectory(fgaliceFile);
+                               fTE->SetAutoFlush(-10*1024*1024);
+                               //CopygALICE(am);
+                       }
+               }
+               else if (fTEinput != am->GetMCtruthEventHandler()->GetTree()) {
+                       fTEinput =
+                               ((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
+                       fTEinput->CopyAddresses(fTE);
+               }
+       }
+
+       // check whether the ESD input tree is still the same - NOT NECESSARY
+       //TTree* currentInput =
+       //      AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree();
+       //if (fInputTree != currentInput) { // ESD input tree has changed
+       //      puts("CHANGED INPUT TREE");
+       //      fInputTree = currentInput;
+       //      fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");
+       //      fInputTreeCopy->CopyAddresses(fOutputTree);
+       //}
+
+       AliESDEvent *esdin = dynamic_cast<AliESDEvent*>(InputEvent());
+       if (!esdin || !fInputTree || !fInputTreeCopy || !fOutputTree || !fTEinput
+           || !fTE) {
+               // check whether everything is ready
+               PostData(1, fOutputTree);
+               PostData(2, fStatsFlow);
+               PostData(3, fSkimmingList);
+               if (fDoMC) PostData(4, fTE);
+               return;
+       }
+
+       fStatsFlow->Fill(1.); // Ready <= bin name
+
+       //============================================================================
+       // event selection, every event passing the next lines is written to the new
+       // ESD stream
+       //
+       // modify this lines in order to adjust the ESD selection to your needs
+
+       // determine current gap condition
+       fCurrentGapCondition = AliCDMesonUtils::GetGapConfig(esdin, 0x0, 0x0, 0x0,
+                                                            0x0, 0x0, 0x0, 0x0, 0x0);
+
+       TRandom3 rnd(0);
+       if (!(AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) ||
+            (fRequestedGapCondition & fCurrentGapCondition) || fReduceGapEvents) {
+               // check whether the event satifies the required gap condition and whether
+               // all double-gap events will be stored
+
+               if (fReduceGapEvents &&
+                   (AliCDMesonBase::kBitBaseLine & fCurrentGapCondition) &&
+                   !(fRequestedGapCondition & fCurrentGapCondition)) {
+                       // not all double gap events can be stored, but the gap determination
+                       // was successful and at least a according to the less stringent condition
+                       // a double gap was found
+                       if (!(fRefinedGapCondition & fCurrentGapCondition)) {
+                               // event fulfilled the more stringent gap and hence is will be stored
+                               fStatsFlow->Fill(5.);
+                       }
+                       else if (rnd.Rndm() > (1.-fReducedGapFraction)) {
+                               // randomly select a fraction of
+                               fStatsFlow->Fill(4.); // Double Gap Event <= bin name
+                       }
+                       else {
+                               // reject event (no refined gap and not randomly selected)
+                               fStatsFlow->Fill(3.); // Event Rejected <= bin name
+                               PostData(1, fOutputTree);
+                               PostData(2, fStatsFlow);
+                               PostData(3, fSkimmingList);
+                               if (fDoMC) PostData(4, fTE);
+                               return;
+                       }
+               }
+               else if (rnd.Rndm() > (1.-fNoGapFraction)) { // randomly selected
+                       fStatsFlow->Fill(2.); // Control Sample Event <= bin name
+               }
+               else {
+                       fStatsFlow->Fill(3.); // Event Rejected <= bin name
+                       PostData(1, fOutputTree);
+                       PostData(2, fStatsFlow);
+                       PostData(3, fSkimmingList);
+                       if (fDoMC) PostData(4, fTE);
+                       return;
+               }
+       }
+       else {
+               fStatsFlow->Fill(4.); // Double Gap Event <= bin name
+       }
+       // end of event selection
+       //============================================================================
+
+       // load the current event in the cloned input tree
+       //
+       // unfortunately the Entry() gives the entry number with in the current input
+       // tree, not within the chain for MC data
+       // => use GetReadEntry() of the ESD tree (value return by Entry() for real
+       //data)
+       fInputTreeCopy->GetEvent(fInputTree->GetReadEntry());
+       fOutputTree->Fill(); // fill the current event into the new ESD stream
+       //printf("Entry()=%lld\nEntries=%lld\nChainOffset=%lld\nReadEntry=%lld\n",
+       //       Entry(), fInputTreeCopy->GetEntries(),
+       //       fInputTreeCopy->GetChainOffset(), fInputTree->GetReadEntry());
+
+       // MC specific stuff ---------------------------------------------------------
+       AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+
+       if (fDoMC) { // do the MC related part of the skimming
+               Long64_t iEntry = fOutputTree->GetEntries() - 1;
+               AliMCEventHandler* mcHandler =
+                       (AliMCEventHandler*)am->GetMCtruthEventHandler();
+               TTree* treeK = mcHandler->TreeK();
+               TTree* treeTR = mcHandler->TreeTR();
+               if (!treeK) AliFatal("TreeK not found!");
+               if (!treeTR) AliFatal("TreeTR not found!");
+
+               fWrkDir = gDirectory;
+               gDirectory = fKinematicsFile;
+               TString dirName = TString(Form("Event%ld", (Long_t)iEntry));
+               gDirectory->mkdir(dirName);
+               gDirectory->cd(dirName);
+               TTree* outputTreeK = treeK->CloneTree(); // copy the whole tree
+               outputTreeK->Write();
+               treeK->CopyAddresses(outputTreeK, kTRUE); // separate clone again
+               treeK->GetListOfClones()->Remove((TObject*)outputTreeK);
+               outputTreeK->Delete();
+               outputTreeK = 0x0;
+
+               gDirectory = fTrackRefsFile;
+               gDirectory->mkdir(dirName);
+               gDirectory->cd(dirName);
+               TTree* outputTreeTR = treeTR->CloneTree(); // copy the whole tree
+               outputTreeTR->Write();
+               treeTR->CopyAddresses(outputTreeTR, kTRUE); // separate clone again
+               treeTR->GetListOfClones()->Remove((TObject*)outputTreeTR);
+               outputTreeTR->Delete();
+               outputTreeTR = 0x0;
+               gDirectory = fWrkDir;
+
+               TTree* currentTEinput =
+                       ((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
+               if (fTEinput != currentTEinput) {
+                       // when the input file is changed, the kinematics tree address changes
+                       // hence our output tree fTE has to be linked to the new one!
+                       fTEinput = currentTEinput;
+                       fTEinput->CopyAddresses(fTE);
+               }
+               fTE->Fill(); // fill MC event headers
+       }
+
+       fStatsFlow->Fill(6.); // Tree Filled <= bin name
+
+       // remember which event we stored
+       fEventNumberInFile = esdin->GetEventNumberInFile();
+       fRunNumber = esdin->GetRunNumber();
+       fFileName = fInputTree->GetCurrentFile()->GetName();
+       fEventTime = esdin->GetTimeStamp();
+       fSkimmingList->Fill();
+
+       PostData(1, fOutputTree);
+       PostData(2, fStatsFlow);
+       PostData(3, fSkimmingList);
+       if (fDoMC) PostData(4, fTE);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskCDskimESD::UserCreateOutputObjects()
+{
+       // Create output objects.
+       fStatsFlow = (TH1*)(new TH1F("skimmingStatsFlow", "", 7, 0., 7.));
+       TAxis* x = fStatsFlow->GetXaxis();
+       x->SetBinLabel(1, "UserExec(...)");
+       x->SetBinLabel(2, "Ready");
+       x->SetBinLabel(3, "Control Sample Event");
+       x->SetBinLabel(4, "Rejected Event");
+       x->SetBinLabel(5, "Double Gap Event");
+       x->SetBinLabel(6, "Refined Double Gap Event");
+       x->SetBinLabel(7, "Tree Filled");
+
+       PostData(2, fStatsFlow);
+
+       // TODO implement fSkimmingList
+       fSkimmingList = new TTree("SkimmingList", "SkimmingList");
+       fSkimmingList->Branch("FileName", &fFileName);
+       fSkimmingList->Branch("EventNumberInFile", &fEventNumberInFile);
+       fSkimmingList->Branch("RunNumber", &fRunNumber);
+       fSkimmingList->Branch("TimeStamp", &fEventTime);
+       PostData(3, fSkimmingList);
+
+       // Get input information
+       AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+       fInputTree = am->GetInputEventHandler()->GetTree();
+       if (!fInputTree) {
+               puts("AliAnalysisTaskCDskimESD: input tree not found\n");
+               return;
+       }
+       fInputTreeCopy = (TTree*)fInputTree->Clone("esdTree");
+
+
+       // prevent the task from being run on proof
+       if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
+           AliAnalysisManager::kProofAnalysis) {
+               AliFatal("AliAnalysisTaskCDskimESD: cannot be run on PROOF!");
+       }
+
+       TFile *file = 0x0;
+       file = OpenFile(1); // open file for the first output slot
+       if (!file) {
+               puts("AliAnalysisTaskCDskimESD: file not opened\n");
+                       return;
+       }
+
+       fOutputTree = fInputTreeCopy->CloneTree(0);
+       if (!fOutputTree) {
+               puts("AliAnalysisTaskCDskimESD: cloning tree not successful\n");
+               return;
+       }
+
+       file->SetCompressionLevel(5); // caused a crash
+       fOutputTree->SetDirectory(file);
+       fOutputTree->SetAutoFlush(-10*1024*1024);
+
+       if (fDoGapCond) {
+               if(!fOutputTree->GetBranch("gapCondition")) {
+                       fOutputTree->Branch("gapCondition", &fCurrentGapCondition);
+               }
+       }
+
+       PostData(1, fOutputTree);
+
+       // process the mc information
+       if(fDoMC) {
+               fWrkDir = gDirectory; // save the old directory
+               fKinematicsFile = new TFile("Kinematics.root", "update");
+               gDirectory = fWrkDir; // restore the old directory
+
+               fTrackRefsFile = new TFile("TrackRefs.root", "update");
+               gDirectory = fWrkDir; // restore the old directory
+
+               fTEinput =
+                       ((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetTree();
+               if (!fTEinput) {
+                       // this trick is done in order to post an output tree, although the tree
+                       // which is cloned is not available during the UserCreateOutputObjects()
+                       fTE = new TTree("TE", "TE");
+               }
+               else {
+                       fTE = fTEinput->CloneTree(0);
+                       //CopygALICE(am); // not needed, not properly working/tested
+               }
+               fgaliceFile = OpenFile(4);
+               if (!file) {
+                       puts("AliAnalysisTaskCDskimESD: galice.root file not opened\n");
+                       return;
+               }
+               fgaliceFile->SetCompressionLevel(5);
+               fTE->SetDirectory(fgaliceFile);
+               fTE->SetAutoFlush(-10*1024*1024);
+
+
+               PostData(4, fTE);
+       }
+}
+
+
+//______________________________________________________________________________
+/* // not needed => not yet working / properly tested
+void AliAnalysisTaskCDskimESD::CopygALICE(AliAnalysisManager* am)
+{
+       //
+       // copy all objects contained in the galice root excluding the TE tree
+       //
+
+       TString galiceInput =
+               *(((AliMCEventHandler*)am->GetMCtruthEventHandler())->GetInputPath());
+       galiceInput += "galice.root";
+       printf("galiceInput=%s\n", galiceInput.Data());
+       fWrkDir = gDirectory;
+       TFile* input = new TFile(galiceInput.Data(), "READ");
+
+       input->cd();
+       gDirectory->ls();
+
+       //loop on all entries of of the input directory
+       TKey *key;
+       TIter nextkey(input->GetListOfKeys());
+       input->GetListOfKeys()->ls();
+       while ((key = (TKey*)nextkey())) {
+               TString name = key->GetName();
+               printf("name=%s\n", name.Data());
+               const TString classname = key->GetClassName();
+               printf("classname=%s\n", classname.Data());
+               TClass *cl = gROOT->GetClass(classname);
+               if (!cl) continue;
+               if (!name.Contains("TE") && !classname.Contains("AliRun")) {
+                       TObject *obj = key->ReadObj();
+                       fgaliceFile->cd();
+                       obj->Write();
+                       delete obj;
+                       input->cd();
+               }
+               else if (cl->InheritsFrom("AliRunLoader")) {
+                       // TODO write the runloader stuff
+                       //AliRunLoader* rl = AliRunLoader::Instance();
+                       //printf("FILENAME=%s\n", rl->GetFileName().Data());
+               }
+               fgaliceFile->ls();
+       }
+       fgaliceFile->SaveSelf(kTRUE);
+       fWrkDir->cd();
+}
+*/
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.h b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliAnalysisTaskCDskimESD.h
new file mode 100644 (file)
index 0000000..b2c7ef3
--- /dev/null
@@ -0,0 +1,105 @@
+/*************************************************************************
+* Copyright(c) 1998-1999, 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$
+//
+// Task to skim ESD files, basic idea found at PWGGA/EMCALTasks as
+// AliEsdSkimTask. This version is modified to produce a exact copy of the ESD
+// stream containing only events which are of interest for the central
+// diffractive analysis (PWGUD) and a control sample. In order to be able to
+// process MC data as well, the usual concept using an output slot is for
+// output data is violated and an output file is directly created without
+// a corresponding output slot.
+//
+// Author:
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#ifndef ALIANALYSISTASKCDSKIMESD_H
+#define  ALIANALYSISTASKCDSKIMESD_H
+
+#include "AliAnalysisTaskSE.h"
+
+// $Id$
+
+class TTree;
+class TH1;
+class TFile;
+class TDirectory;
+
+class AliAnalysisTaskCDskimESD : public AliAnalysisTaskSE {
+public:
+       AliAnalysisTaskCDskimESD(const char *opt, Bool_t reduceGapEvents);
+       AliAnalysisTaskCDskimESD();
+       virtual ~AliAnalysisTaskCDskimESD();
+
+       void UserExec(Option_t *opt);
+       void UserCreateOutputObjects();
+
+       void SetGapCondition(Int_t cond) { fRequestedGapCondition=cond; }
+       void SetRefinedGapCondition(Int_t cond) {fRefinedGapCondition=cond; }
+
+protected:
+       // variables concerning the reconstructed information in the ESDs
+       TTree*      fInputTree;             // input tree
+       TTree*      fInputTreeCopy;         // input tree copy
+       TTree*      fOutputTree;            // output tree
+
+       // variables needed for the reproduction of the MC information
+       Bool_t      fDoMC;                  // do we process MC information?
+       TDirectory* fWrkDir;                // directory the analysis task used before
+                                           // it was manipulated to store the
+                                           // kinemactics truth
+       TFile*      fKinematicsFile;        // output file containing the skimmed
+                                           // kinematics information
+       TFile*      fTrackRefsFile;         // output file containing the skimmed
+                                           // TrackRefs information
+       TTree*      fTEinput;               // input tree containing MC eventt headers
+       TTree*      fTE;                    // tree containing the MC event headers
+       TFile*      fgaliceFile;            // pointer to the galice.root output file
+
+       // variables concerning the gaps used for the selection
+       Bool_t      fDoGapCond;             // if true the gap condition for central
+                                           // diffractive analysis are stored
+       Int_t       fCurrentGapCondition;   // gap condition
+       Int_t       fRequestedGapCondition; // gap condition required for selection
+       Double_t    fNoGapFraction;         // fraction of no-gap events to be stored
+       Bool_t      fReduceGapEvents;       // do not store all V0 double-gap events
+                                           // necessary due to high double gap rates
+                                           // in MC productions (LHC10d{1,2,4,4a})
+       Int_t       fRefinedGapCondition;   // gap condition with more stringent
+                                           // selection applied in case
+                                           // fReduceDgEvents is active
+       Double_t    fReducedGapFraction;    // fraction of events stored altough they
+                                           // show only the requested and not the
+                                           // refined gap
+
+       // selection monitoring
+       TH1*        fStatsFlow;             // histogram indicating the number of
+                                           // input events
+       TTree*      fSkimmingList;          // which events of the intial ones are
+                                           // skimmed?
+       TString     fFileName;              // name of the input file
+       Int_t       fEventNumberInFile;     // event index in the ESD file
+       Int_t       fRunNumber;             // run number of the event
+       UInt_t      fEventTime;             // event timestamp
+
+private:
+       //void CopygALICE(AliAnalysisManager* am); // copy objects in galice.root
+       // not implemented:
+       AliAnalysisTaskCDskimESD(const AliAnalysisTaskCDskimESD&);
+       AliAnalysisTaskCDskimESD &operator=(const AliAnalysisTaskCDskimESD&);
+       ClassDef(AliAnalysisTaskCDskimESD, 1); // Esd skimming task
+};
+#endif // ALIANALYSISTASKCDSKIMESD_H
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.cxx b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.cxx
new file mode 100644 (file)
index 0000000..308dbcc
--- /dev/null
@@ -0,0 +1,610 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+//
+// AliDDMesonBase
+// for
+// AliAnalysisTaskDDMeson
+//
+//  Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#include "TH2.h"
+#include "TH2D.h"
+#include "TH1I.h"
+#include "TAxis.h"
+#include "THnSparse.h"
+#include "TString.h"
+#include "TList.h"
+
+#include "AliCDMesonBase.h"
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonBase::GetGapBin(TString tag, const Int_t gapcg,
+                                Bool_t checkCentralActivity /* = kTRUE */)
+{
+       //
+       // retrieve gap topology for a given string
+       //
+
+       tag.ToUpper();
+
+       Bool_t ka = kFALSE, kc = kFALSE;
+
+       if(tag.Contains("V0")){
+               ka = ka || (gapcg & kBitV0A);
+               kc = kc || (gapcg & kBitV0C);
+       }
+       if(tag.Contains("FMD")){
+               ka = ka || (gapcg & kBitFMDA);
+               kc = kc || (gapcg & kBitFMDC);
+       }
+       if(tag.Contains("SPD")){
+               ka = ka || (gapcg & kBitSPDA);
+               kc = kc || (gapcg & kBitSPDC);
+       }
+       if(tag.Contains("TPC")){
+               ka = ka || (gapcg & kBitTPCA);
+               kc = kc || (gapcg & kBitTPCC);
+       }
+       if(tag.Contains("ZDC")){
+               ka = ka || (gapcg & kBitZDCA);
+               kc = kc || (gapcg & kBitZDCC);
+       }
+
+       if(ka && kc)
+               return kBinNG;
+       else{
+               if(!ka && !kc)
+                       if (((gapcg & kBitCentAct) && checkCentralActivity) ||
+                           !checkCentralActivity) {
+                               return kBinDG; // central activity seen (or not required)
+                       }
+                       else {
+                               return kBinNG; // no central activity
+                       }
+               else if(!kc)
+                       return kBinGC;
+               else
+                       return kBinGA;
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::CheckRange(Double_t &var, const Double_t min,
+                                const Double_t max)
+{
+       //
+       // check whether the value is with in the range specified with min and max
+       //
+
+       const Double_t eps = 1e-3;
+       if( var >= max ) var = max - eps;
+       if( var <= min ) var = min + eps;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonBase::GetAxis(TString thntit, TString name)
+{
+       //
+       // get the number of an axis, derived from the ThnSparse title (list of axes)
+       //
+
+       thntit.ToUpper();
+       thntit.ReplaceAll(" ","");
+       const Int_t nmax = 20;
+       TString tits[nmax];
+       Int_t counter = 0;
+       while(thntit.Contains(",")){
+               const Int_t pos = thntit.First(",");
+               tits[counter] = thntit(0, pos);
+               thntit = thntit(pos+1, thntit.Length()-pos);
+               counter++;
+               if(counter>=nmax-1){
+                       printf("freidtlog AliCDMesonBase::GetAxis too small nmax! %d %d\n",
+                              counter, nmax);
+                       return -1; //exit(1); // TODO
+               }
+       }
+       tits[counter++] = thntit;
+       
+       //----------------
+
+       name.ToUpper();
+
+       for(Int_t ii=0; ii<counter; ii++){
+         if( tits[ii] == name )
+                 return ii;
+       }
+       printf("freidtlog AliCDMesonBase::GetAxis !%s! %s not found!\n", name.Data(),
+              thntit.Data());
+       for(Int_t ii=0; ii<counter; ii++){
+               printf("*************** AliCDMesonBase::GetAxis *****************\n");
+               printf("freidtlog AliCDMesonBase::GetAxis %d !%s!\n", ii, tits[ii].Data());
+               printf("\n");
+       }
+       return -1; //exit(1); // TODO
+}
+
+
+//------------------------------------------------------------------------------
+TString AliCDMesonBase::GetTitleMother()
+{
+       //
+       // get the title of ThnMother (= list of axes)
+       //
+
+       TString title = "Ncombined, CombCh, CombPID, V0, FMD, SPD, TPC,";
+       title += " Mass, Pt, CTS, OA, DaughterPt, TrackResiduals, VertexZinRng,";
+       title += " ProcessType, VertexCoincidence, TrackletResiduals";
+       return title;
+}
+
+
+//------------------------------------------------------------------------------
+THnSparseD * AliCDMesonBase::GetThnMother(TString name /* = "CDMeson_Mother" */)
+{
+       // creates the THnSparse called Mother (used for the generation of invariant
+       // mass distributions)
+       //
+
+       const Int_t nbin[] = {
+               fgkNNcombined, fgkNCombCh, fgkNCombPID, fgkNGapConfig, fgkNGapConfig,
+               fgkNGapConfig, fgkNGapConfig, fgkNMass, fgkNMotherPt, fgkNCTS, fgkNOA,
+               fgkNDaughterPt, fgkNTrackResiduals, fgkNVertexZinRng, fgkNProcessType,
+               fgkNVertexCoincidence, fgkNTrackletResiduals
+       };
+       const Double_t binmin[] = {
+               fgkMinNcombined, fgkMinCombCh, fgkMinCombPID, fgkMinGapConfig,
+               fgkMinGapConfig, fgkMinGapConfig, fgkMinGapConfig, fgkMinMass,
+               fgkMinMotherPt, fgkMinCTS, fgkMinOA, fgkMinDaughterPt,
+               fgkMinTrackResiduals, fgkMinVertexZinRng, fgkMinProcessType,
+               fgkMinVertexCoincidence, fgkMinTrackletResiduals
+       };
+       const Double_t binmax[] = {
+               fgkMaxNcombined, fgkMaxCombCh, fgkMaxCombPID, fgkMaxGapConfig,
+               fgkMaxGapConfig, fgkMaxGapConfig, fgkMaxGapConfig, fgkMaxMass,
+               fgkMaxMotherPt, fgkMaxCTS, fgkMaxOA, fgkMaxDaughterPt,
+               fgkMaxTrackResiduals, fgkMaxVertexZinRng, fgkMaxProcessType,
+               fgkMaxVertexCoincidence, fgkMaxTrackletResiduals
+       };
+       
+       const Int_t npar = sizeof(nbin)/sizeof(Int_t);
+
+       return new THnSparseD(name.Data(), GetTitleMother(), npar, nbin, binmin,
+                             binmax);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::FillThnMother(THnSparseD *thn, const Double_t vNch,
+                                   const Double_t vCombCh,
+                                   const Double_t vCombPID, const Double_t vV0,
+                                   const Double_t vFMD, const Double_t vSPD,
+                                   const Double_t vTPC, const Double_t vMass,
+                                   const Double_t vPt,  const Double_t vOA,
+                                   const Double_t vCTS,
+                                   const Double_t vDaughterPt,
+                                   const Double_t vTrackResiduals,
+                                   const Double_t vVertexZ,
+                                   const Double_t vProcessType,
+                                   const Double_t vVertexCoincidence,
+                                   const Double_t vTrkltResiduals)
+{
+       //
+       // fill ThnMother
+       //
+
+       Double_t var[]={
+               vNch, vCombCh, vCombPID, vV0, vFMD, vSPD, vTPC, vMass, vPt, vOA, vCTS,
+               vDaughterPt, vTrackResiduals, vVertexZ, vProcessType, vVertexCoincidence,
+               vTrkltResiduals
+       };
+       const Int_t nv = sizeof(var)/sizeof(Double_t);
+       if(nv!=thn->GetNdimensions()){
+               printf("AliCDMesonBase::FillThnMother nv error!! %d %d\n", nv,
+                      thn->GetNdimensions());
+               return; //exit(1);
+       }
+
+       CheckRange(var[7], fgkMinMass, fgkMaxMass);
+       CheckRange(var[8], fgkMinMotherPt, fgkMaxMotherPt);
+       CheckRange(var[11], fgkMinDaughterPt, fgkMaxDaughterPt);
+
+       thn->Fill(var);
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonBase::GetAxisMother(const TString name)
+{
+       //
+       // return axis number corresponding to the name
+       //
+
+       return GetAxis(GetTitleMother(), name);
+}
+
+
+//------------------------------------------------------------------------------
+TString AliCDMesonBase::GetTitleEmptyEvents()
+{
+       //
+       // title / list of axes of the empty event study THnSparse
+       //
+
+       TString title = "EventType, FMD-A, FMD-C, SPD-I-A, SPD-I-C, SPD-O-A, SPD-O-C";
+       title += ", SPDtrkltA, SPDtrkltC";
+       title += ", fmdSum1I, fmdSum2I, fmdSum2O, fmdSum3I, fmdSum3O";
+       //title += ", TPC_all, TPC_diffVertex";
+       return title;
+}
+
+
+//------------------------------------------------------------------------------
+THnSparseD* AliCDMesonBase::GetThnEmptyEvents()
+{
+       // creates the THnSparse called Empty Events
+       // EventType, FMD-A, FMD-C, SPD-I-A, SPD-I-C, SPD-O-A, SPD-O-C, FMD1I..FMD3O
+       // TPC_all, TPC_diffVertex (not used so far)
+
+       const Int_t nbin[] = {
+               fgkNEventType, fgkNMultW, fgkNMultW, fgkNMult, fgkNMult, fgkNMult, fgkNMult,
+               fgkNMult, fgkNMult, fgkNMult, fgkNMult, fgkNMult, fgkNMult, fgkNMult
+               //, fgkNMult, fgkNMult
+       };
+       const Double_t binmin[] = {
+               fgkMinEventType, fgkMinMultW, fgkMinMultW, fgkMinMult, fgkMinMult,
+               fgkMinMult, fgkMinMult, fgkMinMult, fgkMinMult, fgkMinMult, fgkMinMult,
+               fgkMinMult, fgkMinMult, fgkMinMult //, fgkMinMult, fgkMinMult
+       };
+       const Double_t binmax[] = {
+               fgkMaxEventType, fgkMaxMultW, fgkMaxMultW, fgkMaxMult, fgkMaxMult,
+               fgkMaxMult, fgkMaxMult, fgkMaxMult, fgkMaxMult, fgkMaxMult, fgkMaxMult,
+               fgkMaxMult, fgkMaxMult, fgkMaxMult //, fgkMaxMult, fgkMaxMult
+       };
+       
+       const Int_t npar = sizeof(nbin)/sizeof(Int_t);
+
+       return new THnSparseD("CDMeson_EmptyEvents", GetTitleEmptyEvents(), npar,
+                             nbin, binmin, binmax);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::FillThnEmptyEvents(THnSparseD * thn, const Int_t eventType,
+                                        const Int_t multFMDA,
+                                        const Int_t multFMDC,
+                                        const Int_t multSPDIA,
+                                        const Int_t multSPDIC,
+                                        const Int_t multSPDOA,
+                                        const Int_t multSPDOC,
+                                        const Int_t multSPDtrkltA,
+                                        const Int_t multSPDtrkltC,
+                                        const Int_t fmdSum1I,
+                                        const Int_t fmdSum2I,
+                                        const Int_t fmdSum2O,
+                                        const Int_t fmdSum3I,
+                                        const Int_t fmdSum3O/*,
+                                        const Int_t multTPC,
+                                        const Int_t multTPCdiffVertex */)
+{
+       //
+       // Fill ThnEmptyEvents
+       //
+
+       Double_t var[]={
+               eventType, multFMDA, multFMDC, multSPDIA, multSPDIC, multSPDOA,
+               multSPDOC, multSPDtrkltA, multSPDtrkltC, fmdSum1I, fmdSum2I, fmdSum2O,
+               fmdSum3I, fmdSum3O//, multTPC, multTPCdiffVertex
+       };
+       const Int_t nv = sizeof(var)/sizeof(Double_t);
+       if(nv!=thn->GetNdimensions()){
+               printf("AliCDMesonBase::FillThnEmptyEvents nv error!! %d %d\n", nv,
+                      thn->GetNdimensions());
+               return; //exit(1); // TODO
+       }
+
+       thn->Fill(var);
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonBase::GetAxisEmptyEvents(const TString name)
+{
+       //
+       // return axis number corresponding to the name
+       //
+
+       return GetAxis(GetTitleEmptyEvents(), name);
+}
+
+
+//------------------------------------------------------------------------------
+TString AliCDMesonBase::GetTitleMultiplicity()
+{
+       //
+       // get title of the multiplicity study ThnSparse
+       //
+
+       TString title = "Nch, Nsoft, Ncombined, V0, FMD, SPD, TPC, NresidualTracks";
+       title += ", NresidualTracklets, VertexZ, VerticesDistance, ProcessType";
+       return title;
+}
+
+
+//------------------------------------------------------------------------------
+THnSparseD* AliCDMesonBase::GetThnMultiplicity()
+{
+       //
+       // creates the THnSparse called Multiplicity
+       //
+
+       const Int_t nbin[] = {
+               fgkNNch, fgkNNsoft, fgkNNcomb,fgkNGapConfig, fgkNGapConfig,
+               fgkNGapConfig, fgkNGapConfig, fgkNNresidualTracks, fgkNNresidualTracklets,
+               fgkNVertexZ, fgkNVerticesDistance, fgkNProcessType
+       };
+       const Double_t binmin[] = {
+               fgkMinNch, fgkMinNsoft, fgkMinNcomb, fgkMinGapConfig, fgkMinGapConfig,
+               fgkMinGapConfig, fgkMinGapConfig, fgkMinNresidualTracks,
+               fgkMinNresidualTracklets, fgkMinVertexZ, fgkMinVerticesDistance,
+               fgkMinProcessType
+       };
+       const Double_t binmax[] = {
+               fgkMaxNch, fgkMaxNsoft, fgkMaxNcomb, fgkMaxGapConfig, fgkMaxGapConfig,
+               fgkMaxGapConfig, fgkMaxGapConfig, fgkMaxNresidualTracks,
+               fgkMaxNresidualTracklets, fgkMaxVertexZ, fgkMaxVerticesDistance,
+               fgkMaxProcessType
+       };
+
+       const Int_t npar = sizeof(nbin)/sizeof(Int_t);
+
+       return new THnSparseD("CDMeson_Multiplicity", GetTitleMultiplicity(), npar,
+                             nbin, binmin, binmax);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::FillThnMultiplicity(THnSparseD *thn, const Double_t vNch,
+                                         const Double_t vNsoft,
+                                         const Double_t vNcombined,
+                                         const Double_t vV0,
+                                         const Double_t vFMD,
+                                         const Double_t vSPD,
+                                         const Double_t vTPC,
+                                         const Double_t vNresidualTracks,
+                                         const Double_t vNresidualTracklets,
+                                         const Double_t vVertexZ,
+                                         const Double_t vVerticesDistance,
+                                         const Double_t vProcessType)
+{
+       // fill ThnMultiplicity
+       // input list copied from GetTitle
+       // var[] copied from GetTitle
+
+       Double_t var[]={
+               vNch, vNsoft, vNcombined, vV0, vFMD, vSPD, vTPC, vNresidualTracks,
+               vNresidualTracklets, vVertexZ, vVerticesDistance, vProcessType
+       };
+       const Int_t nv = sizeof(var)/sizeof(Double_t);
+       if(nv!=thn->GetNdimensions()){
+               printf("AliCDMesonBase::FillThnMultiplicity nv error!! %d %d\n", nv,
+                      thn->GetNdimensions());
+               return; //exit(1);
+       }
+
+       thn->Fill(var);
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonBase::GetAxisMultiplicity(const TString name)
+{
+       //
+       // return axis number corresponding to the name
+       //
+
+       return GetAxis(GetTitleMultiplicity(), name);
+}
+
+
+//------------------------------------------------------------------------------
+TH1F* AliCDMesonBase::GetHistStatsFlow()
+{
+       //
+       // setup the stats flow histogram
+       //
+
+       TH1F *hist = new TH1F("c00_statsFlow", "", AliCDMesonBase::kBinLastValue,
+                       0, AliCDMesonBase::kBinLastValue);
+       TAxis* axis = hist->GetXaxis();
+       axis->SetBinLabel(AliCDMesonBase::kBinTotalInput+1, "total Input");
+       axis->SetBinLabel(AliCDMesonBase::kBinGoodInput+1, "good ESDs");
+       axis->SetBinLabel(AliCDMesonBase::kBinV0OR+1, "V0-OR");
+       axis->SetBinLabel(AliCDMesonBase::kBinV0AND+1, "V0-AND");
+       axis->SetBinLabel(AliCDMesonBase::kBinEventsAfterCuts+1, "after cuts");
+       axis->SetBinLabel(AliCDMesonBase::kBinEventsWithOutPileUp+1, "w/o pile up");
+       axis->SetBinLabel(AliCDMesonBase::kBinv0Gap+1, "with V0 DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinv0fmdGap+1, "with V0-FMD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinv0fmdspdGap+1,
+                         "with V0-FMD-SPD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinv0fmdspdtpcGap+1,
+                         "with V0-FMD-SPD-TPC DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinv0fmdspdtpczdcGap+1,
+                         "with V0-FMD-SPD-TPC-ZDC DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinfmdGap+1, "with FMD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinspdGap+1, "with SPD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBintpcGap+1, "with TPC DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBintpcspdGap+1, "width TPC-SPD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBintpcspdfmdGap+1,
+                         "width TPC-SPD-FMD DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBintpcspdfmdv0Gap+1,
+                         "width TPC-SPD-FMD-V0 DG gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinspdfmdGap+1, "with SPD FMD gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinspdfmdv0Gap+1, "with SPD FMD V0 gap");
+       axis->SetBinLabel(AliCDMesonBase::kBinTwoTrackEvents+1, "with two tracks");
+       axis->SetBinLabel(AliCDMesonBase::kBinThreeTrackEvents+1,
+                         "with three tracks");
+       axis->SetBinLabel(AliCDMesonBase::kBinPionEvents+1, "with two pions");
+       axis->SetBinLabel(AliCDMesonBase::kBinKaonEvents+1, "with two kaons");
+       axis->SetBinLabel(AliCDMesonBase::kBinProtonEvents+1, "with two proton");
+       axis->SetBinLabel(AliCDMesonBase::kBinElectronEvents+1, "with two electron");
+       axis->SetBinLabel(AliCDMesonBase::kBinUnknownPIDEvents+1, "with unknown PID");
+       axis->SetBinLabel(AliCDMesonBase::kBinResidualTracks+1,
+                         "without residual tracks");
+       axis->SetBinLabel(AliCDMesonBase::kBinResidualTracklets+1,
+                         "without residual tracklets");
+       axis->SetBinLabel(AliCDMesonBase::kBinCDonlyEvents+1, "CD only events");
+       return hist;
+}
+
+
+//------------------------------------------------------------------------------
+TH2F* AliCDMesonBase::GetHistPIDStudies(TString name)
+{
+       //
+       // setup the PID studies histogram
+       //
+
+       TH2F *hist = new TH2F(name.Data(), ";particle 1;particle 2",
+                             AliCDMesonBase::kBinPIDUnknown,
+                             AliCDMesonBase::kBinPionE,
+                             AliCDMesonBase::kBinPIDUnknown+1,
+                             AliCDMesonBase::kBinPIDUnknown,
+                             AliCDMesonBase::kBinPionE,
+                             AliCDMesonBase::kBinPIDUnknown+1);
+       TAxis* x = hist->GetXaxis();
+       TAxis* y = hist->GetYaxis();
+       x->SetBinLabel(AliCDMesonBase::kBinPionE, "#pi (ex)");
+       x->SetBinLabel(AliCDMesonBase::kBinPion, "#pi");
+       x->SetBinLabel(AliCDMesonBase::kBinSinglePion, "-");
+       x->SetBinLabel(AliCDMesonBase::kBinKaonE, "K (ex)");
+       x->SetBinLabel(AliCDMesonBase::kBinKaon, "K");
+       x->SetBinLabel(AliCDMesonBase::kBinSingleKaon, ",");
+       x->SetBinLabel(AliCDMesonBase::kBinProtonE, "p (ex)");
+       x->SetBinLabel(AliCDMesonBase::kBinProton, "p");
+       x->SetBinLabel(AliCDMesonBase::kBinSingleProton, "_");
+       x->SetBinLabel(AliCDMesonBase::kBinElectronE, "e (ex)");
+       x->SetBinLabel(AliCDMesonBase::kBinElectron, "e");
+       x->SetBinLabel(AliCDMesonBase::kBinSingleElectron, ".");
+       x->SetBinLabel(AliCDMesonBase::kBinPIDUnknown, "X");
+       y->SetBinLabel(AliCDMesonBase::kBinPionE, "#pi (ex)");
+       y->SetBinLabel(AliCDMesonBase::kBinPion, "#pi");
+       y->SetBinLabel(AliCDMesonBase::kBinSinglePion, "-");
+       y->SetBinLabel(AliCDMesonBase::kBinKaonE, "K (ex)");
+       y->SetBinLabel(AliCDMesonBase::kBinKaon, "K");
+       y->SetBinLabel(AliCDMesonBase::kBinSingleKaon, ",");
+       y->SetBinLabel(AliCDMesonBase::kBinProtonE, "p (ex)");
+       y->SetBinLabel(AliCDMesonBase::kBinProton, "p");
+       y->SetBinLabel(AliCDMesonBase::kBinSingleProton, "_");
+       y->SetBinLabel(AliCDMesonBase::kBinElectronE, "e (ex)");
+       y->SetBinLabel(AliCDMesonBase::kBinElectron, "e");
+       y->SetBinLabel(AliCDMesonBase::kBinSingleElectron, ".");
+       y->SetBinLabel(AliCDMesonBase::kBinPIDUnknown, "X");
+       return hist;
+}
+
+
+//------------------------------------------------------------------------------
+TObjArray* AliCDMesonBase::GetHistVZEROStudies(TList* l)
+{
+       //
+       // Create Histograms for the VZERO trigger studies
+       //
+
+       TObjArray* arr = new TObjArray(130);
+
+       TObject* o = 0x0;
+       for (Int_t iPMT = 0; iPMT < 64; ++iPMT) {
+               o = (TObject*)new TH2F(Form("h00_%02d_ADC_TriggerThr", iPMT),
+                                      ";ADC Counts;Trigger Threshold (ADC Counts)",
+                                      400, 0., 4000., 48, 2., 50.);
+               arr->Add(o);
+               l->Add(o);
+       }
+       for (Int_t iPMT = 0; iPMT < 64; ++iPMT) {
+               o = (TObject*)new TH2F(Form("h01_%02d_ADC_Multiplicity", iPMT),
+                                      ";ADC Counts;Multiplicity",
+                                      400, 0., 4000., 250, 0., 250.);
+               arr->Add(o);
+               l->Add(o);
+       }
+       /*
+         // not of use for pp data in 2010
+       o = (TObject*)new TH2F("h02_TriggerChargeA_Trigger",
+                              ";Trigger Charge A;Trigger Decision A",
+                              250, 0., 5000., 2, 0., 2.);
+       arr->Add(o);
+       l->Add(o);
+       o = (TObject*)new TH2F("h03__TriggerChargeC_Trigger",
+                              ";Trigger Charge C;Trigger Decision C",
+                              250, 0., 5000., 2, 0., 2.);
+       arr->Add(o);
+       l->Add(o);
+       */
+
+       return arr;
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::GetGapTriggers(THnSparseI* gaprun, Int_t gapCondition,
+                                    const Int_t run, Double_t& triggers,
+                                    Double_t& total)
+{
+       // determine the number of certain triggers with a gap in the detectors
+       // specified by gapCondition in run and the total number of events
+       // surviving all cuts (including pile-up rejection)
+
+       triggers = 0;
+       total = 0;
+       Int_t nTuple[] = { gaprun->GetAxis(0)->FindBin(run), 0 };
+       for (Int_t i = 0; i< kBitGapMax; ++i) {
+               nTuple[1] = i;
+               Double_t temp = gaprun->GetBinContent(nTuple);
+               if (!(i & gapCondition)) {
+                       triggers += temp;
+               }
+               total += temp;
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonBase::GetNoGapTriggers(THnSparseI* gaprun, Int_t gapCondition,
+                                      const Int_t run, Double_t& triggers,
+                                      Double_t& total)
+{
+       // determine the number of certain triggers with a NO gap in the detectors
+       // specified by the gapCondition in run and the total number of events
+       // surviving all cuts (including pile-up rejection)
+
+       triggers = 0;
+       total = 0;
+       Int_t nTuple[] = { gaprun->GetAxis(0)->FindBin(run), 0 };
+       for (Int_t i = 0; i< kBitGapMax; ++i) {
+               nTuple[1] = i;
+               Double_t temp = gaprun->GetBinContent(nTuple);
+               if (i & gapCondition) {
+                       triggers += temp;
+               }
+               total += temp;
+       }
+}
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.h b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonBase.h
new file mode 100644 (file)
index 0000000..3ad9a4c
--- /dev/null
@@ -0,0 +1,376 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+//
+// AliCDMesonBase
+// for
+// AliAnalysisTaskCDMeson
+//
+//
+//  Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#ifndef ALICDMESONBASE_H
+#define ALICDMESONBASE_H
+
+#include "THnSparse.h"
+
+class TH2;
+class TH2F;
+class TH1F;
+class TList;
+class TObjArray;
+
+class AliCDMesonBase
+{
+public:
+       enum{
+               // combined charged
+               kBinPM = 1, // unlike-sign
+               kBinPPMM, // like-sign
+
+               // pid results (used for single particles and pairs)
+               kBinPionE = 1, // pions with exclusion cut
+               kBinPion, // pions without exclusion cut but TPC and TOF
+               kBinSinglePion, // one pion identified, the other particle unidentified
+               kBinKaonE, // same for kaons ...
+               kBinKaon,
+               kBinSingleKaon,
+               kBinProtonE, // protons
+               kBinProton,
+               kBinSingleProton,
+               kBinElectronE, // and electrons
+               kBinElectron,
+               kBinSingleElectron,
+               kBinPIDUnknown, // could not be identified
+
+               // gap Conditions
+               kBinDG = 1, // double gap
+               kBinGC, // single gap c side
+               kBinGA, // single gap a side
+               kBinNG, // no gap
+
+
+               // event types used for the empty event studies
+               kBinEventI = 1, // Beam-Beam Interaction
+               kBinEventA, // Beam from A-side
+               kBinEventC, // Beam from C-side
+               kBinEventAC, // Beam from one side (there is no separation in 2011)
+               kBinEventE, // no beam (empty event)
+               kBinEventUnknown,
+
+               // StatsFlow histogram entries
+               // names for the bins are specified in the .cxx-file
+               kBinTotalInput = 0,
+               kBinGoodInput,
+               kBinV0OR,
+               kBinV0AND,
+               kBinEventsAfterCuts,
+               kBinEventsWithOutPileUp,
+               kBinv0Gap,
+               kBinv0fmdGap,
+               kBinv0fmdspdGap,
+               kBinv0fmdspdtpcGap,
+               kBinv0fmdspdtpczdcGap,
+               kBinfmdGap,
+               kBinspdGap,
+               kBintpcGap,
+               kBintpcspdGap,
+               kBintpcspdfmdGap,
+               kBintpcspdfmdv0Gap,
+               kBinspdfmdGap,
+               kBinspdfmdv0Gap,
+               kBinTwoTrackEvents,
+               kBinThreeTrackEvents,
+               kBinPionEvents,
+               kBinKaonEvents,
+               kBinProtonEvents,
+               kBinElectronEvents,
+               kBinUnknownPIDEvents,
+               kBinResidualTracks,
+               kBinResidualTracklets,
+               kBinCDonlyEvents,
+               kBinLastValue, // used to specify the correct histogram width
+
+
+               // MC event/process types
+               kBinND = 0, // used for data as well, "Non-Diffractive"
+               kBinCD, // central diffractive
+               kBinSD, // single diffractive
+               kBinDD, // double diffractive
+
+
+               // gap-condition bits used in AliAnalysisTaskCDMeson::fGapRun
+               kBitBaseLine = (1<<0),
+
+               kBitV0A  = (1<<1),
+               kBitV0C  = (1<<2),
+               kBitFMDA = (1<<3),
+               kBitFMDC = (1<<4),
+
+               kBitSPDA = (1<<5),
+               kBitSPDC = (1<<6),
+               kBitTPCA = (1<<7),
+               kBitTPCC = (1<<8),
+
+               kBitZDCA   = (1<<9),
+               kBitZDCC   = (1<<10),
+               kBitCentAct = (1<<11),
+               kBitGapMax = (1<<12),
+
+
+               // analysis task status bits
+               // do not change the order in order to be backward compatible!
+               kBitConfigurationSet = (1<<0), // if not set everything is active
+               kBitImposeZDCGap = (1<<1), // excluding the ZDC gap
+               kBitEtaPhiMaps = (1<<2),
+               kBitEtaPhiMapsWithCuts = (1<<3),
+               kBitStatsFlow = (1<<4),
+               kBitMultPerGapHists = (1<<5),
+               kBitRmMultPerGapHists = (1<<6),
+               kBitTHnMother = (1<<7),
+               kBitFastORStudy = (1<<8), // also not active by default
+               kBitHitMapSPD = (1<<9),
+               kBitHitMapFMD = (1<<10),
+               kBitVtxStudies = (1<<11),
+               kBitEEStudy = (1<<12),
+               kBitPIDStudy = (1<<13),
+               kBitMCProcess = (1<<14),
+               kBitFMDsum = (1<<15),
+               kBitSoftTracks = (1<<16),
+               kBitPWAtree = (1<<17),
+               kBitTHnMC = (1<<18),
+               kBitMultResponseMC = (1<<19),
+               kBitMultStudy = (1<<20),
+               kBitReadPreprocessedGap = (1<<21), // not active by default as well
+               kBitReduceGapEvents = (1<<22),
+               kBitVZEROStudy = (1<<23),
+               kBitTPCGapStudy = (1<<24),
+               kBitAllTrackMass = (1<<25),
+               kBitCDonly = (1<<26),
+               kBitFastORmultStudy = (1<<27),
+               kBitConfigurationVersion = (1<<28) // always set, last bit
+       };
+
+       static Int_t GetGapBin(const TString tag, const Int_t gapcg,
+                              Bool_t checkCentralActivity = kTRUE);
+
+       static THnSparseD* GetThnMother(TString name = "CDMeson_Mother");
+       static void FillThnMother(THnSparseD *thn, const Double_t vNch,
+                                 const Double_t vCombCh, const Double_t vCombPID,
+                                 const Double_t vV0, const Double_t vFMD,
+                                 const Double_t vSPD, const Double_t vTPC,
+                                 const Double_t vMass, const Double_t vPt,
+                                 const Double_t vOA, const Double_t vCTS,
+                                 const Double_t vDaughterPt,
+                                 const Double_t vTrackResiduals,
+                                 const Double_t vVertexZ,
+                                 const Double_t vProcessType,
+                                 const Double_t vVertexCoincidence,
+                                 const Double_t vTrkltResiduals);
+       static Int_t GetAxisMother(const TString name);
+
+       static THnSparseD* GetThnEmptyEvents();
+       static void FillThnEmptyEvents(THnSparseD *thn, const Int_t eventType,
+                                      const Int_t multFMDA, const Int_t multFMDC,
+                                      const Int_t multSPDIA, const Int_t multSPDIC,
+                                      const Int_t multSPDOA, const Int_t multSPDOC,
+                                      const Int_t multSPDtrklA,
+                                      const Int_t multSPDtrklC, const Int_t fmdSum1I,
+                                      const Int_t fmdSum2I, const Int_t fmdSum20,
+                                      const Int_t fmdSum3I, const Int_t fmdSum3O/*,
+                                      const Int_t multTPC,
+                                      const Int_t multTPCdiffVertex */);
+       static Int_t GetAxisEmptyEvents(const TString name);
+
+       static THnSparseD* GetThnMultiplicity();
+       static void FillThnMultiplicity(THnSparseD *thn, const Double_t vNch,
+                                       const Double_t vNsoft,
+                                       const Double_t vNcombined,
+                                       const Double_t vV0, const Double_t vFMD,
+                                       const Double_t vSPD, const Double_t vTPC,
+                                       const Double_t vNresidualTracks,
+                                       const Double_t vNresidualTracklets,
+                                       const Double_t vVertexZ,
+                                       const Double_t vVerticesDistance,
+                                       const Double_t vProcessType);
+       static Int_t GetAxisMultiplicity(const TString name);
+
+       static TH1F* GetHistStatsFlow();
+       static TH2F* GetHistPIDStudies(TString name);
+       static TObjArray* GetHistVZEROStudies(TList* l);
+
+       static void GetGapTriggers(THnSparseI* gaprun, Int_t gapCondition,
+                                  const Int_t run, Double_t& triggers,
+                                  Double_t& total);
+
+       static void GetNoGapTriggers(THnSparseI* gaprun, Int_t gapCondition,
+                                    const Int_t run, Double_t& triggers,
+                                    Double_t& total);
+private:
+       static void CheckRange(Double_t &var, const Double_t min, const Double_t max);
+       static Int_t GetAxis(TString thntit, TString name);
+       static TString GetTitleMother();
+       static TString GetTitleEmptyEvents();
+       static TString GetTitleMultiplicity();
+
+       //== INVARIANT MASS DISTRIBUTIONS (ThnMother) ================================
+       //-- event characteristics
+       // number of charged primary particles - combined => can be w/ soft tracks
+       static const Int_t fgkNNcombined      = 2; // number of bins
+       static const Double_t fgkMinNcombined = 2; // lower border
+       static const Double_t fgkMaxNcombined = 4; // upper border (#bins + lower border)
+
+       // unlike sign or like sign charged?
+       static const Int_t fgkNCombCh      = 2; // kBinPPMM = number of bins
+       static const Double_t fgkMinCombCh = 1.; // kBinPM = lower border
+       static const Double_t fgkMaxCombCh = 3.; // kBinPPMM + kBinPM = upper border
+
+       // two track PID, take care the enum changed from time to time
+       static const Int_t fgkNCombPID      = 13; // kBinPIDUnknown = number of bins
+       static const Double_t fgkMinCombPID = 1.; // kBinPionE = lower border
+       static const Double_t fgkMaxCombPID = 14.;// ...  = upper border
+
+       // gap configuration, used for different detectors
+       static const Int_t fgkNGapConfig      = 4; // kBinNG = number of bins
+       static const Double_t fgkMinGapConfig = 1.; // kBinDG = lower border
+       static const Double_t fgkMaxGapConfig = 5.; // kBinNG + kBinDG = upper border
+
+       //-- mother kinematics
+       // invariant mass of the two-track system
+       static const Int_t fgkNMass      = 1024; // number of bins
+       static const Double_t fgkMinMass = 0.; // lower border
+       static const Double_t fgkMaxMass = 5.12; // upper border
+
+       // transverse momentum of the two-track system
+       static const Int_t fgkNMotherPt      = 128;//10; // number of bins
+       static const Double_t fgkMinMotherPt = 0.; // lower border
+       static const Double_t fgkMaxMotherPt = .64; //6.4 //2; // upper border
+
+       // cosine theta* (opening angle of the two daugther tracks in the
+       // centre-of-mass system of the two-track/mother system)
+       // **no meaning full variable**
+       static const Int_t fgkNCTS      = 2; // number of bins
+       static const Double_t fgkMinCTS = -1.; // lower border
+       static const Double_t fgkMaxCTS = -0.9; // upper border
+
+       // opening angle in the lab frame
+       static const Int_t fgkNOA      = 20; // number of bins
+       static const Double_t fgkMinOA = -1.; // lower border
+       static const Double_t fgkMaxOA = 1.; // upper border
+
+       //-- daughter kinematics
+       // transverse momentum of one of the two daughter particles
+       // (randomly selected)
+       static const Int_t fgkNDaughterPt      = 128; // number of bins
+       static const Double_t fgkMinDaughterPt = 0.; // lower border
+       static const Double_t fgkMaxDaughterPt = 6.4; // upper border
+
+       // pseudo rapidity of one of the two daughter particles
+       // (randomly selected)
+       //static const Int_t fgkNDaughterEta      = 64; // number of bins
+       //static const Double_t fgkMinDaughterEta = -1.28; // lower border
+       //static const Double_t fgkMaxDaughterEta =  1.28; // upper border
+
+       //-- Event quality information
+       // boolean values to reduce output size
+
+       // are there tracks in addition to the ones selected using AliCDMesonTracks
+       // (ITSTPC, ITSsa, ITSpureSA) (0 = no, 1 = yes)
+       static const Int_t fgkNTrackResiduals      = 2; // number of bins
+       static const Double_t fgkMinTrackResiduals = 0.; // lower border
+       static const Double_t fgkMaxTrackResiduals = 2.; // upper border
+
+       // vertex with in +/-4cm (0 = no, 1 = yes)
+       static const Int_t fgkNVertexZinRng      = 2; // number of bins
+       static const Double_t fgkMinVertexZinRng = 0.; // lower border
+       static const Double_t fgkMaxVertexZinRng = 2.; // upper border
+
+       // are the vertices from SPD and tracks within 0.5cm? (0 = no, 1 = yes)
+       static const Int_t fgkNVertexCoincidence      = 2; // number of bins
+       static const Double_t fgkMinVertexCoincidence = 0.; // lower border
+       static const Double_t fgkMaxVertexCoincidence = 2.; // upper border
+
+       // are there SPD tracklets which are not assigned to tracks? (0 = no, 1 = yes)
+       static const Int_t fgkNTrackletResiduals      = 2; // number of bins
+       static const Double_t fgkMinTrackletResiduals = 0.; // lower border
+       static const Double_t fgkMaxTrackletResiduals = 2.; // upper border
+
+       //-- MC event information
+       static const Int_t fgkNProcessType      = 4; // kBinDD = number of bins
+       static const Double_t fgkMinProcessType = 0.; // kBinND = lower border
+       static const Double_t fgkMaxProcessType = 4.; // kBinDD = upper border
+
+
+       //== EMPTY EVENT STUDY (ThnEmptyEvents) ======================================
+       // event type
+       static const Int_t fgkNEventType      = 5; // kBinEventE = number of bins (5)
+       static const Double_t fgkMinEventType = 1.; // kBinEventI = lower border (1)
+       static const Double_t fgkMaxEventType = 6.; // kBinEventE+kBinEventI = u. b.
+
+       // multiplicities (reused for different detectors and ways of counting)
+       static const Int_t fgkNMult      = 32; // number of bins
+       static const Double_t fgkMinMult = 0.; // lower border
+       static const Double_t fgkMaxMult = 31.; // upper border
+
+       // multplicities - extended range
+       // (reused for different detectors and ways of counting)
+       static const Int_t fgkNMultW      = 64; // number of bins
+       static const Double_t fgkMinMultW = 0; // lower border
+       static const Double_t fgkMaxMultW = 63; // upper border
+
+       //== MULTIPLICITY STUDY (TnnMultiplicity) ====================================
+       // number of ITSTPC tracks in event
+       static const Int_t fgkNNch      = 51; // number of bins
+       static const Double_t fgkMinNch = 0.; // lower border
+       static const Double_t fgkMaxNch = 51.; // upper border
+
+       // number of ITS standalone tracks in event
+       static const Int_t fgkNNsoft      = 11; // number of bins
+       static const Double_t fgkMinNsoft = 0.; // lower border
+       static const Double_t fgkMaxNsoft = 11.; // upper border
+
+       // combined multiplicity
+       static const Int_t fgkNNcomb      = 61; // number of bins
+       static const Double_t fgkMinNcomb = 0.; // lower border
+       static const Double_t fgkMaxNcomb = 61.; // upper border
+
+       // gap configuration is reused from THnMother
+
+       // number of residual tracks
+       static const Int_t fgkNNresidualTracks      = 11; // number of bins
+       static const Double_t fgkMinNresidualTracks = 0.; // lower border
+       static const Double_t fgkMaxNresidualTracks = 11.; // upper border
+
+       // number of residual tracklets
+       static const Int_t fgkNNresidualTracklets      = 21; // number of bins
+       static const Double_t fgkMinNresidualTracklets = 0.; // lower border
+       static const Double_t fgkMaxNresidualTracklets = 21.; // upper border
+
+       // vertex z-position
+       static const Int_t fgkNVertexZ      = 20; // number of bins
+       static const Double_t fgkMinVertexZ = -10.; // lower border
+       static const Double_t fgkMaxVertexZ = 10.; // upper border
+
+       // SPD and track vertex distance
+       static const Int_t fgkNVerticesDistance      = 10; // number of bins
+       static const Double_t fgkMinVerticesDistance = 0.; // lower border
+       static const Double_t fgkMaxVerticesDistance = 5.; // upper border
+
+       // MC process type is resued from THnMother
+};
+
+#endif
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.cxx b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.cxx
new file mode 100644 (file)
index 0000000..4aa4ad2
--- /dev/null
@@ -0,0 +1,1425 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+//
+// AliCDMesonUtils
+// for
+// AliAnalysisTaskCDMeson
+//
+//  Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TGeoMatrix.h>
+#include <THnSparse.h>
+#include <TLorentzVector.h>
+#include <TRandom3.h>
+#include <TParticle.h>
+#include <TObjArray.h>
+
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliCDBStorage.h"
+#include "AliESDEvent.h"
+#include "AliPIDResponse.h"
+#include "AliVTrack.h"
+#include "AliVParticle.h"
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDFMD.h"
+#include "AliESDVZERO.h"
+#include "AliGeomManager.h"
+#include "AliITSAlignMille2Module.h"
+#include "AliITSsegmentationSPD.h"
+#include "AliMultiplicity.h"
+#include "AliPIDResponse.h"
+#include "AliSPDUtils.h"
+#include "AliTriggerAnalysis.h"
+#include "AliVZEROTriggerData.h"
+#include "AliVZEROCalibData.h"
+
+#include "AliAODTracklets.h"
+#include "AliAODEvent.h"
+
+#include "AliCDMesonBase.h"
+#include "AliCDMesonTracks.h"
+
+#include "AliCDMesonUtils.h"
+
+
+//==============================================================================
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetMassPtCtsOA(const Int_t pid,
+                                     const TVector3* momenta[], Double_t & mass,
+                                     Double_t &pt, Double_t &cts, Double_t &oa)
+{
+       //
+       // Get Mass Pt Cts OA
+       //
+
+       Double_t tmpp1[3], tmpp2[3];
+       momenta[0]->GetXYZ(tmpp1);
+       momenta[1]->GetXYZ(tmpp2);
+
+       Double_t masstrk = -999;
+       if(pid==AliCDMesonBase::kBinPionE || pid==AliCDMesonBase::kBinPion ||
+          pid==AliCDMesonBase::kBinSinglePion)
+               masstrk = AliPID::ParticleMass(AliPID::kPion);
+       else if(pid==AliCDMesonBase::kBinKaonE || pid==AliCDMesonBase::kBinKaon ||
+               pid==AliCDMesonBase::kBinSingleKaon)
+               masstrk = AliPID::ParticleMass(AliPID::kKaon);
+       else if(pid==AliCDMesonBase::kBinProtonE || pid==AliCDMesonBase::kBinProton ||
+               pid==AliCDMesonBase::kBinSingleProton)
+               masstrk = AliPID::ParticleMass(AliPID::kProton);
+       else if(pid==AliCDMesonBase::kBinElectronE ||
+               pid==AliCDMesonBase::kBinElectron ||
+               pid==AliCDMesonBase::kBinSingleElectron)
+               masstrk = AliPID::ParticleMass(AliPID::kElectron);
+       else
+               masstrk = AliPID::ParticleMass(AliPID::kPion);
+
+       const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masstrk, masstrk,
+                                                 cts);
+       mass = sumv.M();
+       pt = sumv.Pt();
+
+       oa = GetOA(tmpp1, tmpp2);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetPWAinfo(const Int_t pid, const AliVTrack *trks[],
+                                 Float_t& theta, Float_t& phi, Float_t& mass,
+                                 Float_t momentum[])
+{
+       // transforms the coordiante system to the helicity frame and determines
+       // invariant mass and 3-vector of the two track system
+       //
+
+       Double_t tmpp1[3], tmpp2[3]; // 3-vectors of the daughter tracks
+       trks[0]->GetPxPyPz(tmpp1);
+       trks[1]->GetPxPyPz(tmpp2);
+
+       // determine mass of the daugther tracks
+       Double_t masstrk = -999;
+       if(pid==AliCDMesonBase::kBinPion || pid==AliCDMesonBase::kBinSinglePion)
+               masstrk = AliPID::ParticleMass(AliPID::kPion);
+       else if(pid==AliCDMesonBase::kBinKaon || pid==AliCDMesonBase::kBinSingleKaon)
+               masstrk = AliPID::ParticleMass(AliPID::kKaon);
+       else if(pid==AliCDMesonBase::kBinProton ||
+               pid==AliCDMesonBase::kBinSingleProton)
+               masstrk = AliPID::ParticleMass(AliPID::kProton);
+       else if(pid==AliCDMesonBase::kBinElectron ||
+               pid==AliCDMesonBase::kBinSingleElectron)
+               masstrk = AliPID::ParticleMass(AliPID::kElectron);
+       else
+               masstrk = AliPID::ParticleMass(AliPID::kPion);
+
+       TLorentzVector va, vb; // 4-vectors of the two tracks
+       va.SetXYZM(tmpp1[0], tmpp1[1], tmpp1[2], masstrk);
+       vb.SetXYZM(tmpp2[0], tmpp2[1], tmpp2[2], masstrk);
+
+       TLorentzVector sumv = va+vb; // 4-vector of the pair
+       TVector3 sumXZ(sumv.X(), 0, sumv.Z()); // its projection to the xz-plane
+
+       // mass and momentum
+       mass = sumv.M();
+       momentum[0] = sumv.X();
+       momentum[1] = sumv.Y();
+       momentum[2] = sumv.Z();
+
+       // coordinate axis vectors
+       const TVector3 x(1.,0,0);
+       const TVector3 y(0,1.,0);
+       const TVector3 z(0,0,1.);
+
+       const Double_t alpha = -z.Angle(sumXZ);
+
+       sumXZ.Rotate(alpha, y);
+       sumv.Rotate(alpha, y);
+       va.Rotate(alpha, y);
+       vb.Rotate(alpha, y);
+
+       const Double_t beta = z.Angle(sumv.Vect());
+       sumv.Rotate(beta, x);
+       va.Rotate(beta, x);
+       vb.Rotate(beta, x);
+
+
+       const TVector3 bv = -sumv.BoostVector();
+       va.Boost(bv);
+       vb.Boost(bv);
+
+       // angles
+       theta = va.Theta();
+       phi = va.Phi();
+}
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetEventType(const AliESDEvent *ESDEvent)
+{
+       // checks of which type a event is:
+       // beam-beam interaction (I), beam-gas (A/C), empty (E)
+       //
+
+       TString firedTriggerClasses = ESDEvent->GetFiredTriggerClasses();
+
+       if (firedTriggerClasses.Contains("CINT1A-ABCE-NOPF-ALL")) { // A
+               return AliCDMesonBase::kBinEventA;
+       }
+       if (firedTriggerClasses.Contains("CINT1C-ABCE-NOPF-ALL")) { // C
+               return AliCDMesonBase::kBinEventC;
+       }
+       if (firedTriggerClasses.Contains("CINT1B-ABCE-NOPF-ALL")) { // I
+               return AliCDMesonBase::kBinEventI;
+       }
+       if (firedTriggerClasses.Contains("CINT1-E-NOPF-ALL")) { // E
+               return AliCDMesonBase::kBinEventE;
+       }
+       if (firedTriggerClasses.Contains("CDG5-E")) { // E
+               return AliCDMesonBase::kBinEventE;
+       }
+       if (firedTriggerClasses.Contains("CDG5-I")) { // I
+               return AliCDMesonBase::kBinEventI;
+       }
+       if (firedTriggerClasses.Contains("CDG5-AC")) { // AC
+               return AliCDMesonBase::kBinEventAC;
+       }
+       return AliCDMesonBase::kBinEventUnknown;
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::SwapTrack(const AliVTrack* trks[])
+{
+       //
+       // swap two esdtracks, needed since only the properties of one a stored and
+       // AliRoot sorts them
+       //
+
+       TRandom3 tmprd(0);
+       if(tmprd.Rndm()>0.5)
+               return;
+
+       const AliVTrack* tmpt = trks[0];
+       trks[0]=trks[1];
+       trks[1]=tmpt;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetCombCh(const AliVTrack * trks[])
+{
+       //
+       // get combination of charges
+       //
+
+       const Double_t ch1 = trks[0]->Charge();
+       const Double_t ch2 = trks[1]->Charge();
+
+       if(ch1*ch2<0){
+               return AliCDMesonBase::kBinPM;
+       }
+       else 
+               return AliCDMesonBase::kBinPPMM;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetCombPID(AliPIDResponse* pid, const AliVTrack* trks[],
+                                  const Int_t mode, TH2* comb2trkPID /*= 0x0 */)
+{
+       //
+       // Get PID for a two track system (data)
+       //
+
+       if (!pid){
+               return AliCDMesonBase::kBinPIDUnknown;
+       }
+
+       // get PID for the single tracks
+       Int_t kpid[2];
+       for(Int_t ii=0; ii<2; ii++){
+               kpid[ii] = GetPID(pid, trks[ii], mode);
+       }
+
+       // PID QA histogram
+       if (comb2trkPID) {
+               comb2trkPID->Fill(kpid[0], kpid[1]);
+       }
+
+       return CombinePID(kpid);
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetCombPID(const TParticle* particles[],
+                                  TH2 *comb2trkPID /* = 0x0 */)
+{
+       //
+       // Get PID for a two track system (MC)
+       //
+
+       // get PID for the single tracks
+       Int_t kpid[2];
+       for(Int_t ii=0; ii<2; ii++){
+               kpid[ii] = GetPID(particles[ii]->GetPdgCode());
+       }
+
+       // PID QA histogram
+       if (comb2trkPID) {
+               comb2trkPID->Fill(kpid[0], kpid[1]);
+       }
+
+       return CombinePID(kpid);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetSPDTrackletMult(const AliESDEvent *ESDEvent,
+                                         Int_t& sum, Int_t& forwardA,
+                                         Int_t& forwardC, Int_t& central)
+{
+       // obtain the multiplicity for eta < -.9 (forwardC), eta > 0.9 (fowardA), in
+       // the central barrel (central) and its sum from SPD tracklets with the
+       // AliMultiplicity class
+       // code simular to PWGPP/ITS/AliAnalysisTaskSPD.cxx
+
+       sum = forwardA = forwardC = central = 0; // initialize values
+
+       const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
+       for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
+            iTracklet++) {
+               float_t eta = mult->GetEta(iTracklet);
+               if (eta < -0.9) {
+                       forwardC++;
+               }
+               else if (eta < 0.9) {
+                       central++;
+               }
+               else {
+                       forwardA++;
+               }
+               sum++;
+       }
+}
+
+
+//------------------------------------------------------------------------------
+Bool_t AliCDMesonUtils::CutEvent(const AliESDEvent *ESDEvent, TH1 *hspd,
+                                 TH1 *hpriv, TH1 *hpriVtxPos, TH1 *hpriVtxDist,
+                                 TH2 *hfo, TH1* hfochans, Int_t &kfo,
+                                 Int_t &nip, Int_t &nop, TH1 *hpriVtxX,
+                                 TH1 *hpriVtxY, TH1 *hpriVtxZ)
+{
+       //
+       // CutEvent
+       //
+
+       AliTriggerAnalysis triggerAnalysis;
+       /*
+       //printf("freidtlog active triggers: %s\n",
+       //       ESDEvent->GetHeader()->GetActiveTriggerInputs().Data());
+       //printf("freidtlog fired triggers: %s\n",
+       //       ESDEvent->GetHeader()->GetFiredTriggerInputs().Data());
+       //*/
+       //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
+       //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
+       // returns the number of fired chips in the SPD
+       // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
+       // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
+
+
+       // SPD FastOR cut
+       const Int_t fastORhw = triggerAnalysis.SPDFiredChips(ESDEvent, 1);
+       if (hspd) hspd->Fill(fastORhw);
+       /*
+       if(fastORhw<1)
+               return kFALSE;
+       */
+
+       if (hfochans) {
+               const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
+
+               for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
+                       if (mult->TestFastOrFiredChips(iChipKey)) {
+                               hfochans->Fill((Double_t)iChipKey);
+                       }
+               }
+       }
+
+       if(kfo){ // spd trigger study
+               Int_t nfoctr[10];
+               GetNFO(ESDEvent, "[1.0]", nfoctr, 0x0, 0x0);
+               nip = nfoctr[kInnerPixel];
+               nop = nfoctr[kOuterPixel];
+
+               if(AliCDMesonBase::GetGapBin("V0", GetV0(ESDEvent)) ==
+                   AliCDMesonBase::kBinDG) {
+                       hfo->Fill(nip, nop);
+               }
+
+               if(nop<2)
+                       kfo = 1;
+               else
+                       kfo = nop;
+
+               if(kfo>=10)
+                       kfo=9;
+       }
+
+       // collision vertex cut
+       // A cut in XY is implicitly done during the reconstruction by constraining
+       // the vertex to the beam diamond.
+
+       // Primary vertex
+       Bool_t kpr0 = kTRUE;
+       const AliESDVertex *vertex = ESDEvent->GetPrimaryVertexTracks();
+       if(vertex->GetNContributors()<1) {
+               // SPD vertex
+               vertex = ESDEvent->GetPrimaryVertexSPD();
+               if(vertex->GetNContributors()<1) {
+                       // NO VERTEX, SKIP EVENT
+                       kpr0 = kFALSE;
+               }
+       }
+       const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ()) < 10.);
+       // 10 is the common value, unit: cm
+       if (hpriv) hpriv->Fill(kpriv);
+       if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
+       if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
+       if(!kpriv)
+               return kFALSE;
+
+       if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
+       if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
+       if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
+       return kTRUE;
+}
+
+//------------------------------------------------------------------------------
+Bool_t AliCDMesonUtils::CutEvent(const AliAODEvent *AODEvent, TH1 *hpriv,
+                                 TH1 *hpriVtxX, TH1 *hpriVtxY, TH1 *hpriVtxZ,
+                                 TH1 *hpriVtxPos, TH1 *hpriVtxDist)
+{
+       //
+       // Cut Event for AOD Events, to be combined with the ESD Track Cut
+       //
+
+       // TODO: no idea about fast or yet, to be thought of
+
+       // Primary vertex
+       Bool_t kpr0 = kTRUE;
+       const AliAODVertex *vertex = AODEvent->GetPrimaryVertex();
+       if(vertex->GetNContributors()<1) {
+               // SPD vertex
+               vertex = AODEvent->GetPrimaryVertexSPD();
+               if(vertex->GetNContributors()<1) {
+                       // NO GOOD VERTEX, SKIP EVENT
+                       kpr0 = kFALSE;
+               }
+       }
+       const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ())<10.);
+       // 10 is the common value, unit: cm
+       hpriv->Fill(kpriv);
+       if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
+       if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
+       if(!kpriv)
+               return kFALSE;
+
+       if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
+       if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
+       if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
+       return kTRUE;
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::DoVZEROStudy(const AliESDEvent *ESDEvent,
+                                   TObjArray* hists, const Int_t run)
+{
+       //
+       //
+       // IMPORTANT: the order of the histograms here and in
+       // AliCDMesonBase::GetHistVZEROStudies(..) has to match
+
+       const AliESDVZERO* esdV0 = ESDEvent->GetVZEROData();
+       if (!esdV0) {
+               Printf("ERROR: esd V0  not available");
+               return;
+       }
+
+       // determine trigger decision
+       AliTriggerAnalysis triggerAnalysis;
+       //const Bool_t khw = kFALSE;
+
+       //Float_t v0a = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
+       //             AliTriggerAnalysis::kV0BB) ? 1. : 0.;
+       //Float_t v0c = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
+       //             AliTriggerAnalysis::kV0BB) ? 1. : 0.;
+
+       // obtain OCDB objects
+       AliCDBManager *man = AliCDBManager::Instance();
+       TString cdbpath;
+       if (man->IsDefaultStorageSet()) {
+               const AliCDBStorage *dsto = man->GetDefaultStorage();
+               cdbpath = TString(dsto->GetBaseFolder());
+       }
+       else { //should not be used!
+               man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
+               cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
+       }
+       man->SetSpecificStorage("VZERO/Trigger/Data",cdbpath);
+       man->SetSpecificStorage("VZERO/Calib/Data",cdbpath);
+       man->SetRun(run);
+
+       AliCDBEntry *ent1 = man->Get("VZERO/Trigger/Data");
+       AliVZEROTriggerData *trigData = (AliVZEROTriggerData*)ent1->GetObject();
+       if (!trigData) {
+               printf("freidtlog failed loading VZERO trigger data from OCDB\n");
+               return;
+       }
+
+       AliCDBEntry *ent2 = man->Get("VZERO/Calib/Data");
+       AliVZEROCalibData *calData = (AliVZEROCalibData*)ent2->GetObject();
+       if (!calData) {
+               printf("freidtlog failed loading VZERO calibration data from OCDB\n");
+               return;
+       }
+
+       // fill histograms
+       Int_t pmtHist = 0;
+       Int_t iPMT = 0;
+       for (iPMT = 0; iPMT < 64; ++iPMT) {
+               Int_t board   = AliVZEROCalibData::GetBoardNumber(iPMT);
+               Int_t channel = AliVZEROCalibData::GetFEEChannelNumber(iPMT);
+               ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
+                                                     trigData->GetPedestalCut(0, board, channel));
+       //                                      calData->GetDiscriThr(iPMT));
+       }
+       pmtHist = iPMT;
+       for (iPMT = 0; iPMT < 64; ++iPMT) {
+               ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
+                                                     esdV0->GetMultiplicity(iPMT));
+       }
+       /*
+       // not used in 2010 for pp
+       ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeA(), v0a);
+       ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeC(), v0c);
+       */
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetGapConfig(const AliESDEvent *ESDEvent,
+                                    TH2 *hitMapSPDinner, TH2 *hitMapSPDouter,
+                                    TH2 *hitMapSPDtrklt, TH2 *hitMapFMDa,
+                                    TH2 *hitMapFMDc, TH1 **fmdSums,
+                                    TH2 *TPCGapDCAaSide, TH2 *TPCGapDCAcSide)
+{
+       //
+       // GetGapConfigAndTracks
+       //
+       // retrieves the gap configuration of a track and returns it as
+       // an bit vector
+       // kBaseLine ensures, that this event is valid
+       // + is equivalent to | in this case
+       Int_t gapConfig = AliCDMesonBase::kBitBaseLine + GetV0(ESDEvent)
+               + GetFMD(ESDEvent, hitMapFMDa, hitMapFMDc, fmdSums)
+               + GetSPD(ESDEvent, hitMapSPDinner, hitMapSPDouter, hitMapSPDtrklt);
+       if (gapConfig == AliCDMesonBase::kBitBaseLine) {
+               gapConfig += GetTPC(ESDEvent, TPCGapDCAaSide, TPCGapDCAcSide);
+       }
+       else {
+               gapConfig += GetTPC(ESDEvent, 0x0, 0x0);
+       }
+       if (GetFastORmultiplicity(ESDEvent) > 0) {
+               gapConfig += AliCDMesonBase::kBitCentAct;
+       }
+       return gapConfig; // + GetZDC(ESDEvent);
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::FillEtaPhiMap(const AliVEvent *event,
+                                    const AliCDMesonTracks* tracks, TH2 *map,
+                                    TH2 *map_c)
+{
+       //
+       // Fills the eta phi information about all tracks and about the tracks surving
+       // the track cuts (CutTrack) into two separate histograms
+       //
+
+       // all tracks
+       for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); ++itrk) {
+               if (AliVParticle *trk = event->GetTrack(itrk)) {
+                       if (map) {
+                               map->Fill(trk->Eta(), trk->Phi());
+                       }
+               }
+       }
+
+       // tracks that survived the cuts
+       for (Int_t itrk = 0; itrk < tracks->GetCombinedTracks(); ++itrk) {
+               if (map_c) {
+                       map_c->Fill(tracks->GetTrack(itrk)->Eta(), tracks->GetTrack(itrk)->Phi());
+               }
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetMultFMD(const AliESDEvent *ESDEvent, Int_t& fmdA,
+                                 Int_t& fmdC, Float_t  *fmdSums)
+{
+       //
+       // Multiplicity seen by FMD
+       //
+       // WARNING: this function is only working with a modified AliRoot so far
+
+#ifdef STD_ALIROOT
+       fmdA = FMDHitCombinations(ESDEvent, 0);
+       fmdC = FMDHitCombinations(ESDEvent, 1);
+#else
+       AliTriggerAnalysis triggerAnalysis;
+       triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
+       triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide, &fmdA);
+       triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide, &fmdC);
+#endif
+
+       if (fmdSums) {
+               const AliESDFMD* fmdData =
+                       (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
+               for (UShort_t det=1; det<=3; det++) {
+                       Int_t nRings = (det == 1 ? 1 : 2);
+                       for (UShort_t ir = 0; ir < nRings; ir++) {
+                               Char_t   ring = (ir == 0 ? 'I' : 'O');
+                               UShort_t nsec = (ir == 0 ? 20  : 40);
+                               UShort_t nstr = (ir == 0 ? 512 : 256);
+                               for (UShort_t sec =0; sec < nsec;  sec++) {
+                                       for (UShort_t strip = 0; strip < nstr; strip++) {
+                                               Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
+                                               if (mult == AliESDFMD::kInvalidMult) continue;
+
+                                               if (det == 1 && ring == 'I') fmdSums[0] += mult;
+                                               else if (det == 2 && ring == 'I') fmdSums[1] += mult;
+                                               else if (det == 2 && ring == 'O') fmdSums[2] += mult;
+                                               else if (det == 3 && ring == 'I') fmdSums[3] += mult;
+                                               else if (det == 3 && ring == 'O') fmdSums[4] += mult;
+                                       }
+                               }
+                       }
+               }
+       }
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetMultSPD(const AliESDEvent * ESDEvent, Int_t& spdIA,
+                                 Int_t& spdIC, Int_t& spdOA, Int_t& spdOC)
+{
+       //
+       // Retrieves the multiplicity seen by the SPD FastOR
+       //
+
+       Int_t nfoctr[10];
+       GetNFO(ESDEvent, "]0.9[", nfoctr, 0x0, 0x0);
+
+       spdIA = nfoctr[kIPA];
+       spdIC = nfoctr[kIPC];
+       spdOA = nfoctr[kOPA];
+       spdOC = nfoctr[kOPC];
+}
+
+
+//==============================================================================
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetV0(const AliESDEvent * ESDEvent)
+{
+       //
+       //GetV0
+       //
+
+       AliTriggerAnalysis triggerAnalysis;
+       const Bool_t khw = kFALSE;
+       const Bool_t v0A =
+               (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
+                AliTriggerAnalysis::kV0BB);
+       const Bool_t v0C =
+               (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
+                AliTriggerAnalysis::kV0BB);
+
+       return v0A * AliCDMesonBase::kBitV0A + v0C * AliCDMesonBase::kBitV0C;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetFMD(const AliESDEvent *ESDEvent, TH2 *hitMapFMDa,
+                              TH2 *hitMapFMDc, TH1 **fmdSums)
+{
+       //
+       // GetFMD
+       //
+
+       AliTriggerAnalysis triggerAnalysis;
+       triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
+       const Bool_t fmdA =
+               triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide);
+       const Bool_t fmdC =
+               triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide);
+
+       //printf("FR - GetFMD\n");
+
+       // prepartions for a charge summation algorithm
+       Bool_t hitMaps = (Bool_t)(hitMapFMDa && hitMapFMDc);
+       Bool_t calcSum = fmdSums ?
+               (fmdSums[0] && fmdSums[1] && fmdSums[2] && fmdSums[3] && fmdSums[4]) :
+               kFALSE; // are the histograms defined?
+       Float_t sum[] = { 0., 0., 0., 0., 0. };
+       // summed multiplicity in the FMD detectors
+
+       // hit map generation
+       if (hitMaps || calcSum) {
+               const AliESDFMD* fmdData =
+                       (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
+
+               for (UShort_t det=1; det<=3;det++) {
+                       Int_t nRings = (det == 1 ? 1 : 2);
+                       for (UShort_t ir = 0; ir < nRings; ir++) {
+                               Char_t   ring = (ir == 0 ? 'I' : 'O');
+                               UShort_t nsec = (ir == 0 ? 20  : 40);
+                               UShort_t nstr = (ir == 0 ? 512 : 256);
+                               for (UShort_t sec =0; sec < nsec;  sec++) {
+                                       for (UShort_t strip = 0; strip < nstr; strip++) {
+                                               Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
+                                               if (mult == AliESDFMD::kInvalidMult) continue;
+
+                                               if (calcSum) {
+                                                       if (det == 1 && ring == 'I') sum[0] += mult;
+                                                       else if (det == 2 && ring == 'I') sum[1] += mult;
+                                                       else if (det == 2 && ring == 'O') sum[2] += mult;
+                                                       else if (det == 3 && ring == 'I') sum[3] += mult;
+                                                       else if (det == 3 && ring == 'O') sum[4] += mult;
+                                               }
+
+                                               if (hitMaps) { // care about hit map specific information
+                                                       const Float_t eta = fmdData->Eta(det,ring,sec,strip);
+                                                       const Float_t phi =
+                                                               fmdData->Phi(det,ring,sec,strip) / 180. * TMath::Pi();
+                                                       //printf("FR - GetFMD: %f %f %f\n", eta, phi, mult);
+                                                       if (eta != AliESDFMD::kInvalidEta) {
+                                                               if ((-3.5 < eta) && (eta < -1.5)) {
+                                                                       hitMapFMDc->Fill(eta, phi, mult); // use mult as weight
+                                                               }
+                                                               else if ((1.5 < eta) && (eta < 5.5)) {
+                                                                       hitMapFMDa->Fill(eta, phi, mult); // use mult as weight
+                                                               }
+                                                       }
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+
+       if (calcSum) {
+               //printf("DEBUG -- SUM(%f,%f,%f,%f,%f)\n", sum[0], sum[1], sum[2], sum[3],
+               //       sum[4]);
+               for (UInt_t i = 0; i < 5; i++) { // 
+                       fmdSums[i]->Fill(sum[i]);
+               }
+       }
+
+       return fmdA * AliCDMesonBase::kBitFMDA + fmdC * AliCDMesonBase::kBitFMDC;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetSPD(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
+                              TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt)
+{
+       //
+       // GetSPD
+       //
+
+       Int_t nfoctr[10];
+       GetNFO(ESDEvent, "]0.9[", nfoctr, hitMapSPDinner, hitMapSPDouter);
+       // get multiplicity from fastOR and fill corresponding hit maps
+
+       if (hitMapSPDtrklt) FillSPDtrkltMap(ESDEvent, hitMapSPDtrklt);
+       // fill tracklet hit map
+
+       const Int_t ipA = nfoctr[kIPA]; // inner layer A side
+       const Int_t ipC = nfoctr[kIPC]; // inner layer C side
+       const Int_t opA = nfoctr[kOPA]; // outer layer A side
+       const Int_t opC = nfoctr[kOPC]; // outer layer C side
+
+       const Bool_t spdA = ipA + opA; // A side hit?
+       const Bool_t spdC = ipC + opC; // C side hit?
+
+       return spdA * AliCDMesonBase::kBitSPDA + spdC * AliCDMesonBase::kBitSPDC;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetTPC(const AliESDEvent * ESDEvent, TH2 *TPCGapDCAaSide,
+                              TH2 *TPCGapDCAcSide)
+{
+       //
+       //GetTPC
+       //
+
+       const Double_t etacut = 0.9;
+       Int_t nA = 0;
+       Int_t nC = 0;
+
+       AliESDtrackCuts cuts;
+       cuts.SetMaxDCAToVertexXY(0.1);
+       cuts.SetMaxDCAToVertexZ(2.);
+
+       for(Int_t itrack = 0; itrack < ESDEvent->GetNumberOfTracks(); itrack++){
+               const AliESDtrack* esdtrack = ESDEvent->GetTrack(itrack);
+               Float_t b[2];
+               Float_t bCov[3];
+               esdtrack->GetImpactParameters(b,bCov);
+
+               if (bCov[0]<=0 || bCov[2]<=0) {
+                       printf("freidtlog - Estimated b resolution lower or equal zero!\n");
+                       bCov[0]=0;
+                       bCov[2]=0;
+               }
+
+               Float_t dcaToVertexXY = b[0];
+               Float_t dcaToVertexZ = b[1];
+               if (esdtrack->Eta() > etacut) {
+                       if (cuts.AcceptTrack(esdtrack)) nA++;
+                       if (TPCGapDCAaSide) {
+                               TPCGapDCAaSide->Fill(dcaToVertexXY, dcaToVertexZ);
+                       }
+               }
+               else if (esdtrack->Eta() < -etacut) {
+                       if (cuts.AcceptTrack(esdtrack)) nC++;
+                       if (TPCGapDCAcSide) {
+                               TPCGapDCAcSide->Fill(dcaToVertexXY, dcaToVertexZ);
+                       }
+               }
+       }
+
+       const Bool_t tpcA = nA;
+       const Bool_t tpcC = nC;
+
+       return tpcA * AliCDMesonBase::kBitTPCA + tpcC * AliCDMesonBase::kBitTPCC;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetZDC(const AliESDEvent * ESDEvent)
+{
+       //
+       //GetZDC
+       //
+
+       const Int_t qa = ESDEvent->GetESDZDC()->GetESDQuality();
+       Bool_t zdcA = kFALSE, zdcC = kFALSE;
+       for(Int_t ii=0; ii<6; ii++){
+               if(qa & (1<<ii)){
+                       if(ii<4) zdcA = kTRUE;
+                       else zdcC = kTRUE;
+               }
+       }
+
+       return zdcA * AliCDMesonBase::kBitZDCA + zdcC * AliCDMesonBase::kBitZDCC;
+}
+
+
+//==============================================================================
+//------------------------------------------------------------------------------
+#ifdef STD_ALIROOT
+Int_t AliCDMesonUtils::FMDHitCombinations(const AliESDEvent* ESDEvent,
+                                          Int_t side)
+{
+       //
+       // copy of the FMDHitCombinations function originating from AliTriggerAnalysis
+       //
+       // side == 0 -> A side, side > 0 -> C side
+
+       // workaround for AliESDEvent::GetFMDData is not const!
+       const AliESDFMD* fmdData =
+               (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
+       if (!fmdData)
+               {
+                       puts("AliESDFMD not available");
+                       return -1;
+               }
+
+       Int_t detFrom = (side == 0) ? 1 : 3;
+       Int_t detTo   = (side == 0) ? 2 : 3;
+
+       Int_t triggers = 0;
+       Float_t totalMult = 0;
+       for (UShort_t det=detFrom;det<=detTo;det++) {
+               Int_t nRings = (det == 1 ? 1 : 2);
+               for (UShort_t ir = 0; ir < nRings; ir++) {
+                       Char_t   ring = (ir == 0 ? 'I' : 'O');
+                       UShort_t nsec = (ir == 0 ? 20  : 40);
+                       UShort_t nstr = (ir == 0 ? 512 : 256);
+                       for (UShort_t sec =0; sec < nsec;  sec++) {
+                               for (UShort_t strip = 0; strip < nstr; strip++) {
+                                       Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
+                                       if (mult == AliESDFMD::kInvalidMult) continue;
+
+                                       if (mult > 0.3) // fFMDLowCut
+                                               totalMult = totalMult + mult;
+                                       else {
+                                               if (totalMult > 0.5) // fFMDHitCut
+                                                       triggers++;
+                                       }
+                               }
+                       }
+               }
+       }
+       return triggers;
+}
+#endif
+
+
+//==============================================================================
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::SPDLoadGeom(const Int_t run)
+{
+       // method to get the gGeomanager
+       // it is called at the CreatedOutputObject stage
+       // to comply with the CAF environment
+
+       AliCDBManager *man = AliCDBManager::Instance();
+
+       TString cdbpath;
+       if (man->IsDefaultStorageSet()) {
+               const AliCDBStorage *dsto = man->GetDefaultStorage();
+               cdbpath = TString(dsto->GetBaseFolder());
+               //printf("default was set to: %s\n", cdbpath.Data());
+       }
+       else { //should not be used!
+               // man->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
+               // would be needed on grid
+               man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
+               cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
+       }
+
+       man->SetSpecificStorage("ITS/Align/Data",cdbpath);
+       man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
+       man->SetRun(run);
+
+       AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
+       if (!obj) {
+               printf("freidtlog failed loading geometry object\n");
+               return;
+       }
+       AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
+       AliGeomManager::ApplyAlignObjsFromCDB("ITS"); 
+}
+
+//------------------------------------------------------------------------------
+Bool_t AliCDMesonUtils::SPDLoc2Glo(const Int_t id, const Double_t *loc,
+                                   Double_t *glo) 
+{
+       //
+       //SPDLoc2Glo, do not touch
+       //
+
+       static TGeoHMatrix mat;
+       Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
+       if (vid<0) {
+               printf("freidtlog Did not find module with such ID %d\n",id);
+               return kFALSE;
+       }
+       AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
+       mat.LocalToMaster(loc,glo);
+       return kTRUE;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::CheckChipEta(const Int_t chipKey, const TString scut,
+                                    const Double_t vtxPos[],
+                                    TH2 *hitMapSPDinner, TH2 *hitMapSPDouter)
+{
+       //
+       //CheckChipEta
+       //
+
+       // retrieves the position in eta for a given chip and applies the cut
+       // results:
+       // 0 <= out of range
+       // -1 <= negative pseudo-rapidity position, in range (C-Side)
+       // 1 <= positive pseudo-rapidity position, in range (A-Side)
+       //
+       // scut: "[0.9" or "]0.9", only 3 digits for the value!!
+
+
+       const Bool_t kincl = (scut[0] == '[');
+       const TString cutval = scut(1,3);
+       const Double_t etacut = fabs(cutval.Atof());
+
+       //no eta cut, save time
+       if(kincl && etacut>=2)
+               return kTRUE;
+
+       Int_t etaside = 1;
+       //------------------------------- NOT TO TOUCH ------------------------>>
+       UInt_t module=999, offchip=999;
+       AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
+       UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
+       if(hs<2) offchip = 4 - offchip; // inversion  in the inner layer...
+
+       const Int_t col[]={
+               hs<2? 0 : 31, 
+               hs<2? 31 : 0, 
+               hs<2? 31 : 0, 
+               hs<2? 0 : 31};
+       const Int_t aa[]={0, 0, 255, 255};
+       const AliITSsegmentationSPD seg;
+
+       for(Int_t ic=0; ic<4; ic++){
+               Float_t localchip[3]={0.,0.,0.};
+               seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]);
+               // local coordinate of the chip center
+               //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
+               const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
+               Double_t glochip[3]={0.,0.,0.};
+               if(!SPDLoc2Glo(module,local,glochip)){
+                       return kFALSE;
+               }
+
+               //-------------------------------------------------------------------<<
+
+               const TVector3 pos(glochip[0]-vtxPos[0], glochip[1]-vtxPos[1],
+                                  glochip[2]-vtxPos[2]);
+               //pos.Print();
+
+               if (chipKey < 400) { // inner SPD layer
+                       if (hitMapSPDinner) {
+                               Double_t phi = pos.Phi(); // output in the range -Pi +Pi
+                               if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
+                               const Double_t eta = pos.Eta();
+                               hitMapSPDinner->Fill(eta, phi);
+                       }
+               }
+               else {
+                       if (hitMapSPDouter) { // outer SPD layer
+                               Double_t phi = pos.Phi(); // output in the range -Pi +Pi
+                               if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
+                               const Double_t eta = pos.Eta();
+                               hitMapSPDouter->Fill(eta, phi);
+                       }
+               }
+
+               if( kincl && fabs(pos.Eta()) > etacut)
+                       return kFALSE;
+
+               if(!kincl){
+                       if(fabs(pos.Eta()) < etacut)
+                               return kFALSE;
+                       else if(pos.Eta()<0)
+                               etaside = -1;
+                       else
+                               etaside = 1;
+               }
+       }
+
+       return etaside;
+}
+
+
+//------------------------------------------------------------------------------
+void AliCDMesonUtils::GetNFO(const AliESDEvent *ESDEvent, const TString etacut,
+                             Int_t ctr[], TH2 *hitMapSPDinner,
+                             TH2 *hitMapSPDouter)
+{
+       //
+       // GetNFO
+       //
+       // analyzes the SPD fastOR for a given eta range and returns
+       // an array with the number of hits in:
+
+       Int_t ninner=0; // inner layer
+       Int_t nouter=0; // outer layer
+       Int_t ipA = 0; // inner layer A side
+       Int_t ipC = 0; // inner layer C side
+       Int_t opA = 0; // outer layer A side
+       Int_t opC = 0; // outer layer C side
+
+       const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
+
+       // position of the primary vertex
+       Double_t tmp[3] = { 0., 0., 0. };
+       ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
+       Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
+
+
+       for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
+               if(mult->TestFastOrFiredChips(iChipKey)){
+                       // here you check if the FastOr bit is 1 or 0
+                       const Int_t iseta = CheckChipEta(iChipKey, etacut, vtxPos, hitMapSPDinner,
+                                                        hitMapSPDouter);
+                       if(iseta==0)
+                               continue;
+
+                       if(iChipKey<400) {
+                               ninner++;  // here you count the FastOr bits in the inner layer
+                               if(iseta>0)
+                                       ipA ++;
+                               else
+                                       ipC ++;
+                       }
+                       else {
+                               nouter++;  // here you count the FastOr bits in the outer layer
+                               if(iseta>0)
+                                       opA ++;
+                               else
+                                       opC ++;
+                       }
+               }
+       }
+
+       ctr[kInnerPixel]= ninner;
+       ctr[kOuterPixel]= nouter;
+       ctr[kIPA]= ipA;
+       ctr[kIPC]= ipC;
+       ctr[kOPA]= opA;
+       ctr[kOPC]= opC;
+
+       return;
+}
+
+
+//--------------------------------------------------------------------------
+Int_t AliCDMesonUtils::GetFastORmultiplicity(const AliESDEvent* ESDEvent)
+{
+       // determine the number of fired fastOR chips in both layers within
+       // -0.9 < eta < 0.9
+       //
+
+       const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
+
+       // position of the primary vertex
+       Double_t tmp[3] = { 0., 0., 0. };
+       ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
+       Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
+
+       Int_t multiplicity = 0;
+
+       for (Int_t iChipKey=0; iChipKey < 1200; iChipKey++) {
+               if(mult->TestFastOrFiredChips(iChipKey)){
+                       // here you check if the FastOr bit is 1 or 0
+                       const Int_t iseta = CheckChipEta(iChipKey, "[0.9]", vtxPos, 0x0, 0x0);
+                       if(iseta==0)
+                               continue;
+                       else
+                               ++multiplicity;
+               }
+       }
+
+       return multiplicity;
+}
+
+
+//--------------------------------------------------------------------------
+void AliCDMesonUtils::FillSPDtrkltMap(const AliVEvent* event,
+                                      TH2 *hitMapSPDtrklt)
+{
+       //
+       // fill eta phi map of SPD tracklets
+       //
+
+       if (hitMapSPDtrklt) {
+               if (TString(event->ClassName()) == "AliESDEvent") {
+                       const AliMultiplicity *mult = ((AliESDEvent*)event)->GetMultiplicity();
+                       for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
+                            iTracklet++) {
+                               Double_t eta = mult->GetEta(iTracklet);
+                               Double_t phi = mult->GetPhi(iTracklet);
+                               hitMapSPDtrklt->Fill(eta, phi);
+                       }
+               }
+               else if (TString(event->ClassName()) == "AliAODEvent") {
+                       const AliAODTracklets* mult = ((AliAODEvent*)event)->GetTracklets();
+                       for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
+                            iTracklet++) {
+                               Double_t eta = -TMath::Log(TMath::Tan(mult->GetTheta(iTracklet)/2.));
+                               Double_t phi = mult->GetPhi(iTracklet);
+                               hitMapSPDtrklt->Fill(eta, phi);
+                       }
+               }
+       }
+}
+
+
+//==========================================================================
+Int_t AliCDMesonUtils::GetPID(AliPIDResponse *pid, const AliVTrack *trk,
+                              const Int_t mode /* = 0 */)
+{
+       // determines PID for ESDs and AODs
+       // 
+       //
+
+       if (!pid) return AliCDMesonBase::kBinPIDUnknown; // no pid available
+
+       Double_t tpcpion = -999., tpckaon = -999., tpcproton = -999.,
+               tpcelectron = -999.;
+       Double_t tofpion = -999., tofkaon = -999., tofproton = -999.,
+               tofelectron = -999.;
+       Double_t itspion = -999., itskaon = -999., itsproton = -999.,
+               itselectron = -999.;
+
+
+       // check whether the track was pure ITS standalone track (has not left TPC)
+       const Bool_t kits = !(trk->GetStatus() & AliESDtrack::kTPCout);
+
+       if (kits) { // do ITS pid
+               const Double_t nin = 3.;
+
+               itspion = pid->NumberOfSigmasITS( trk, AliPID::kPion );
+               itskaon = pid->NumberOfSigmasITS( trk, AliPID::kKaon );
+               itsproton = pid->NumberOfSigmasITS( trk, AliPID::kProton );
+               itselectron = pid->NumberOfSigmasITS( trk, AliPID::kElectron );
+
+               if(fabs(itspion) < nin) // inclusion only
+                       return AliCDMesonBase::kBinPion;
+               else if(fabs(itskaon) < nin)
+                       return AliCDMesonBase::kBinKaon;
+               else if(fabs(itsproton) < nin)
+                       return AliCDMesonBase::kBinProton;
+               else if(fabs(itselectron) < nin)
+                       return AliCDMesonBase::kBinElectron;
+               else{
+                       return AliCDMesonBase::kBinPIDUnknown;
+               }
+       }
+
+       tpcpion     = pid->NumberOfSigmasTPC( trk, AliPID::kPion );
+       tpckaon     = pid->NumberOfSigmasTPC( trk, AliPID::kKaon );
+       tpcproton   = pid->NumberOfSigmasTPC( trk, AliPID::kProton );
+       tpcelectron = pid->NumberOfSigmasTPC( trk, AliPID::kElectron );
+
+       // derive information, whether tof pid is available
+       const Bool_t ka = !(trk->GetStatus() & AliESDtrack::kTOFmismatch); 
+       const Bool_t kb =  (trk->GetStatus() & AliESDtrack::kTOFpid);
+       const Bool_t ktof = ka && kb;
+
+       // retrieve TOF information if available
+       if(ktof){
+               tofpion = pid->NumberOfSigmasTOF(trk, AliPID::kPion);
+               tofkaon = pid->NumberOfSigmasTOF(trk, AliPID::kKaon);
+               tofproton = pid->NumberOfSigmasTOF(trk, AliPID::kProton);
+               tofelectron = pid->NumberOfSigmasTOF(trk, AliPID::kElectron);
+       }
+
+       if (mode == 0) {
+               // 2010 cuts for PID
+               const Double_t nin = 3.;
+
+               // inclusion + exclusion (either TPC or TOF)
+               if((fabs(tpcpion) < nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
+                   && fabs(tpcelectron) > nin) ||
+                  (fabs(tofpion) < nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
+                   && fabs(tofelectron) > nin))
+                       return AliCDMesonBase::kBinPionE;
+               else if((fabs(tpcpion) > nin && fabs(tpckaon) < nin && fabs(tpcproton) > nin
+                        && fabs(tpcelectron) > nin) ||
+                       (fabs(tofpion) > nin && fabs(tofkaon) < nin && fabs(tofproton) > nin
+                        && fabs(tofelectron) > nin))
+                       return AliCDMesonBase::kBinKaonE;
+               else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) < nin
+                        && fabs(tpcelectron) > nin) ||
+                       (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) < nin
+                        && fabs(tofelectron) > nin))
+                       return AliCDMesonBase::kBinProtonE;
+               else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
+                        && fabs(tpcelectron) < nin) ||
+                       (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
+                        && fabs(tofelectron) < nin))
+                       return AliCDMesonBase::kBinElectronE;
+               else if(fabs(tpcpion) < nin && fabs(tofpion) < nin) // inclusion (TPC + TOF)
+                       return AliCDMesonBase::kBinPion;
+               else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
+                       return AliCDMesonBase::kBinKaon;
+               else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
+                       return AliCDMesonBase::kBinProton;
+               else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
+                       return AliCDMesonBase::kBinElectron;
+               else{
+                       return AliCDMesonBase::kBinPIDUnknown;
+               }
+       }
+       else if (mode == 1) {
+               // 2011 cuts for PID in LHC11f
+               // TPC: [-3,5] sigma (pion)
+               // TOF: 3 sigma for all,
+               // only Pion is tuned!
+               // ONLY INCLUSION CUTS NO EXCLUSION
+               const Double_t nin = 3.;
+
+               if(tpcpion < 4. && tpcpion > -2. && fabs(tofpion) < -3. )
+                       return AliCDMesonBase::kBinPion;
+               else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
+                       return AliCDMesonBase::kBinKaon;
+               else if(fabs(tpcproton) < nin &&  fabs(tofproton) < nin)
+                       return AliCDMesonBase::kBinProton;
+               else if(fabs(tpcelectron) < nin &&  fabs(tofelectron) < nin)
+                       return AliCDMesonBase::kBinElectron;
+               else{
+                       return AliCDMesonBase::kBinPIDUnknown;
+               }
+       }
+       return AliCDMesonBase::kBinPIDUnknown;
+}
+
+
+//==========================================================================
+Int_t AliCDMesonUtils::GetPID(const Int_t pdgCode)
+{
+       //
+       // determine particle type based on PDG code
+       //
+
+       if (TMath::Abs(pdgCode) == 211) return AliCDMesonBase::kBinPionE;
+       else if (TMath::Abs(pdgCode) == 321) return AliCDMesonBase::kBinKaonE;
+       else if (TMath::Abs(pdgCode) == 2212) return AliCDMesonBase::kBinProtonE;
+       else if (TMath::Abs(pdgCode) == 11) return AliCDMesonBase::kBinElectronE;
+       else return AliCDMesonBase::kBinPIDUnknown;
+}
+
+
+//------------------------------------------------------------------------------
+Int_t AliCDMesonUtils::CombinePID(const Int_t pid[])
+{
+       //
+       // combine the PID result
+       //
+
+       // determine return value
+       if (pid[0] == pid[1]) { // same result for both tracks
+               return pid[0];
+       }
+       // one track identified with exclusion the other only without
+       else if ((pid[0] == AliCDMesonBase::kBinPionE &&
+                 pid[1] == AliCDMesonBase::kBinPion) ||
+                (pid[1] == AliCDMesonBase::kBinPionE &&
+                 pid[0] == AliCDMesonBase::kBinPion)) {
+               return AliCDMesonBase::kBinPion;
+       }
+       else if ((pid[0] == AliCDMesonBase::kBinKaonE &&
+                 pid[1] == AliCDMesonBase::kBinKaon) ||
+                (pid[1] == AliCDMesonBase::kBinKaonE &&
+                 pid[0] == AliCDMesonBase::kBinKaon)) {
+               return AliCDMesonBase::kBinKaon;
+       }
+       else if ((pid[0] == AliCDMesonBase::kBinProtonE &&
+                 pid[1] == AliCDMesonBase::kBinProton) ||
+                (pid[1] == AliCDMesonBase::kBinProtonE &&
+                 pid[0] == AliCDMesonBase::kBinProton)) {
+               return AliCDMesonBase::kBinProton;
+       }
+       else if ((pid[0] == AliCDMesonBase::kBinElectronE &&
+                 pid[1] == AliCDMesonBase::kBinElectron) ||
+                (pid[1] == AliCDMesonBase::kBinElectronE &&
+                 pid[0] == AliCDMesonBase::kBinElectron)) {
+               return AliCDMesonBase::kBinElectron;
+       }
+       // one track identified and one not
+       else if (((pid[0] == AliCDMesonBase::kBinPionE ||
+                  pid[0] == AliCDMesonBase::kBinPion) &&
+                 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
+                (pid[1] == AliCDMesonBase::kBinPion &&
+                 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
+               return AliCDMesonBase::kBinSinglePion;
+       }
+       else if (((pid[0] == AliCDMesonBase::kBinKaonE ||
+                  pid[0] == AliCDMesonBase::kBinKaon)&&
+                 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
+                (pid[1] == AliCDMesonBase::kBinKaon &&
+                 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
+               return AliCDMesonBase::kBinSingleKaon;
+       }
+       else if (((pid[0] == AliCDMesonBase::kBinProtonE ||
+                  pid[0] == AliCDMesonBase::kBinProton) &&
+                 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
+                (pid[1] == AliCDMesonBase::kBinProton &&
+                 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
+               return AliCDMesonBase::kBinSingleProton;
+       }
+       else if (((pid[0] == AliCDMesonBase::kBinElectronE ||
+                  pid[0] == AliCDMesonBase::kBinElectron) &&
+                 pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
+                (pid[1] == AliCDMesonBase::kBinElectron &&
+                 pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
+               return AliCDMesonBase::kBinSingleElectron;
+       }
+       else
+               return AliCDMesonBase::kBinPIDUnknown;
+}
+
+
+//------------------------------------------------------------------------------
+TLorentzVector AliCDMesonUtils::GetKinematics(const Double_t *pa,
+                                              const Double_t *pb,
+                                              const Double_t ma,
+                                              const Double_t mb, Double_t& cts)
+{
+       //
+       //get kinematics, cts = cos(theta_{#star})
+       //
+
+       TLorentzVector va, vb;
+       va.SetXYZM(pa[0], pa[1], pa[2], ma);
+       vb.SetXYZM(pb[0], pb[1], pb[2], mb);
+       const TLorentzVector sumv = va+vb;
+
+       const TVector3 bv = -sumv.BoostVector();
+
+       va.Boost(bv);
+       vb.Boost(bv);
+
+       // 3-vectors in the restframe of the mother particle
+       const TVector3 pra = va.Vect();
+       const TVector3 prb = vb.Vect();
+       const TVector3 diff = pra - prb; // their difference
+
+       cts = (diff.Mag2()-pra.Mag2()-prb.Mag2()) / (-2.*pra.Mag()*prb.Mag());
+       // cosine theta star, calculated according to the law-of-cosines
+
+       return sumv;
+}
+
+
+//------------------------------------------------------------------------------
+Double_t AliCDMesonUtils::GetOA(const Double_t *pa, const Double_t *pb)
+{
+       //
+       //cosOpeningAngle
+       //
+
+       TVector3 va, vb;
+       va.SetXYZ(pa[0], pa[1], pa[2]);
+       vb.SetXYZ(pb[0], pb[1], pb[2]);
+
+       const TVector3 ua = va.Unit();
+       const TVector3 ub = vb.Unit();
+
+       return ua.Dot(ub);
+}
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.h b/PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.h
new file mode 100644 (file)
index 0000000..53e2e75
--- /dev/null
@@ -0,0 +1,183 @@
+/*************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+//
+// AliCDMesonUtils
+// for
+// AliAnalysisTaskCDMeson
+//
+//
+//  Author:
+//  Xianguo Lu <lu@physi.uni-heidelberg.de>
+//  continued by
+//  Felix Reidt <Felix.Reidt@cern.ch>
+
+#ifndef ALICDMESONUTILS_H
+#define ALICDMESONUTILS_H
+
+#include <THnSparse.h> // forward declaration not possible
+
+
+#define STD_ALIROOT
+// if this is defined the code should run with a standard aliroot
+
+
+class TH1;
+class TH2;
+class TLorentzVector;
+class TParticle;
+class TVector3;
+class TObjArray;
+
+class AliVEvent;
+class AliESDEvent;
+class AliAODEvent;
+class AliESDtrack;
+class AliESDtrackCuts;
+class AliVTrack;
+class AliPIDResponse;
+
+class AliCDMesonTracks;
+
+class AliCDMesonUtils
+{
+public:
+       enum{
+               kInnerPixel = 0,
+               kOuterPixel,
+               kIPA,
+               kIPC,
+               kOPA,
+               kOPC
+       };
+
+       // ESD only
+       //---------
+
+       // cuts for ESD analysis
+       static Bool_t CutEvent(const AliESDEvent *ESDEvent, TH1 *hspd, TH1 *hpriv,
+                              TH1* hpriVtxPos, TH1* hpriVtxDist, TH2 *hfo,
+                              TH1* hfochans, Int_t &kfo, Int_t &nip, Int_t &nop,
+                              TH1* hpriVtxX, TH1* hpriVtxY, TH1* hpriVtxZ);
+
+       // V0-AND/OR
+       static Bool_t V0AND(const AliESDEvent *ESDEvent) {
+               // were both V0 detectors A/C firing?
+               // (4 = AliCDMesonBase::kBinNG)
+               return (GetV0(ESDEvent) == 4);
+       }
+
+       static Bool_t V0OR(const AliESDEvent *ESDEvent) {
+               // did one of the V0 detectors fire?
+               // (1-4 => single or both detectors active)
+               return (GetV0(ESDEvent) > 0);
+       }
+
+       // VZERO study
+       static void DoVZEROStudy(const AliESDEvent *ESDEvent, TObjArray* hists,
+                                const Int_t run);
+
+       // gap determination
+       static Int_t GetGapConfig(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
+                                 TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt,
+                                 TH2 *hitMapFMDa, TH2 *hitMapFMDc,
+                                 TH1 **fmdSums, TH2 *TPCGapDCAaSide,
+                                 TH2 *TPCGapDCAcSide);
+
+       static Int_t GetFastORmultiplicity(const AliESDEvent *ESDEvent);
+
+       static void GetSPDTrackletMult(const AliESDEvent *ESDEvent, Int_t& sum,
+                                      Int_t& forwardA, Int_t& forwardC,
+                                      Int_t& central);
+
+       static void SPDLoadGeom(const Int_t run); // only needed for ESDs, not in AODs
+
+       // functions needed for the empty event study (ESD only)
+       static Int_t GetEventType(const AliESDEvent *ESDEvent);
+       static void GetMultFMD(const AliESDEvent *ESDEvent, Int_t& fmdA,
+                              Int_t& fmdC, Float_t *fmdSums);
+       static void GetMultSPD(const AliESDEvent *ESDEvent, Int_t& spdIA,
+                              Int_t& spdIC, Int_t& spdOA, Int_t& spdOC);
+
+
+       // AOD only
+       //---------
+       static Bool_t CutEvent(const AliAODEvent *AODEvent, TH1 *hpriv, TH1* hpriVtxX,
+                              TH1* hpriVtxY, TH1* hpriVtxZ, TH1* hpriVtxPos,
+                              TH1* hpriVtxDist);
+
+
+       // independent
+       //------------
+       static void FillEtaPhiMap(const AliVEvent *event,
+                                 const AliCDMesonTracks* tracks, TH2 *map,
+                                 TH2 *map_c);
+       static void SwapTrack(const AliVTrack *trks[]);
+       static Int_t GetCombCh(const AliVTrack *trks[]);
+       static Int_t GetCombPID(AliPIDResponse *pid, const AliVTrack *trks[],
+                               const Int_t mode, TH2 *comb2trkPID = 0x0);
+       static Int_t GetCombPID(const TParticle* particles[], TH2 *comb2trkPID = 0x0);
+       static void GetMassPtCtsOA(const Int_t pid, const TVector3* momenta[],
+                                  Double_t& mass, Double_t& pt, Double_t& cts,
+                                  Double_t& oa);
+       static void GetPWAinfo(const Int_t pid, const AliVTrack *trks[],
+                              Float_t& theta, Float_t& phi, Float_t& mass,
+                              Float_t momentum[]);
+       static void FillSPDtrkltMap(const AliVEvent *event, TH2 *hitMapSPDtrklt);
+
+private:
+       // ESD only
+       //---------
+
+       // Gap determination functions
+       static Int_t GetV0(const AliESDEvent *ESDEvent);
+       static Int_t GetFMD(const AliESDEvent *ESDEvent, TH2 *hitMapFMDa,
+                           TH2 *hitMapFMDc, TH1** fmdSums);
+       static Int_t GetSPD(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
+                           TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt);
+       static Int_t GetTPC(const AliESDEvent *ESDEvent, TH2 *TPCGapDCAaSide,
+                           TH2 *TPCGapDCAcSide);
+       static Int_t GetZDC(const AliESDEvent *ESDEvent); // not used so far
+
+#ifdef STD_ALIROOT
+       // helpers for the FMD gap determination
+       static Int_t FMDHitCombinations(const AliESDEvent* ESDEvent, Int_t side);
+#endif
+
+       // helpers for the SPD gap determination
+       static Bool_t SPDLoc2Glo(const Int_t id, const Double_t *loc, Double_t *glo);
+       static Int_t CheckChipEta(const Int_t chipKey, const TString scut,
+                                 const Double_t vtxPos[], TH2 *hitMapSPDinner,
+                                 TH2 *hitMapSPDouter);
+       static void GetNFO(const AliESDEvent *ESDEvent, const TString etacut,
+                          Int_t ctr[], TH2 *hitMapSPDinner, TH2 *hitMapSPDouter);
+
+       // AOD only
+       //---------
+
+
+       // independent
+       //----------
+       static Int_t GetPID(AliPIDResponse *pid, const AliVTrack *trk,
+                           const Int_t mode = 0);
+       static Int_t GetPID(const Int_t pdgCode);
+       static Int_t CombinePID(const Int_t pid[]);
+
+       static TLorentzVector GetKinematics(const Double_t *pa, const Double_t *pb,
+                                           const Double_t ma, const Double_t mb,
+                                           Double_t & cts);
+       static Double_t GetOA(const Double_t *pa, const Double_t *pb);
+};
+
+#endif
diff --git a/PWGUD/DIFFRACTIVE/xsAndTwoProng/runMesonGrid.C b/PWGUD/DIFFRACTIVE/xsAndTwoProng/runMesonGrid.C
new file mode 100644 (file)
index 0000000..e07472f
--- /dev/null
@@ -0,0 +1,399 @@
+// runMeson.C
+//
+// run macro for central diffractive meson analysis
+//
+// Author: Felix Reidt <felix.reidt@cern.ch>
+//
+// based on:
+// ---
+// Template run macro for AliBasicTask.cxx/.h with example layout of
+// physics selections and options, in macro and task.
+//
+// Author: Arvinder Palaha
+//
+class AliAnalysisGrid;
+
+//______________________________________________________________________________
+void runMeson(
+             const char* runtype = "grid", // local, proof or grid
+             const char *gridmode = "full", // Set the run mode (can be "full", "test", "offline", "submit" or "terminate"). Full & Test work for proof
+             const bool bMCphyssel = 0, // 1 = looking at MC truth or reconstructed, 0 = looking at real data
+             const Long64_t nentries = 2000, // for local and proof mode, ignored in grid mode. Set to 1234567890 for all events.
+             const Long64_t firstentry = 0, // for local and proof mode, ignored in grid mode
+             const char *proofdataset = "/alice/data/LHC10c_000120821_p1", // path to dataset on proof cluster, for proof analysis
+             const char *proofcluster = "alice-caf.cern.ch", // which proof cluster to use in proof mode
+             const char *taskname = "CDMeson" // sets name of grid generated macros
+             )
+{
+       // check run type
+       if(runtype != "local" && runtype != "proof" && runtype != "grid"){
+               Printf("\n\tIncorrect run option, check first argument of run macro");
+               Printf("\tint runtype = local, proof or grid\n");
+               return;
+       }
+       Printf("%s analysis chosen",runtype);
+
+       // load libraries (to be optimized)
+       gSystem->Load("libCore.so");
+       gSystem->Load("libTree.so");
+       gSystem->Load("libPhysics");
+       gSystem->Load("libMinuit");
+       gSystem->Load("libProof");
+       gSystem->Load("libmicrocern");
+       gSystem->Load("liblhapdf");
+       gSystem->Load("libpythia6");
+       gSystem->Load("libEG");
+       gSystem->Load("libGeom");
+       gSystem->Load("libVMC");
+       gSystem->Load("libEGPythia6");
+       gSystem->Load("libSTEERBase");
+       gSystem->Load("libESD");
+       gSystem->Load("libCDB");
+       gSystem->Load("libRAWDatabase");
+       gSystem->Load("libRAWDatarec");
+       gSystem->Load("libAOD");
+       gSystem->Load("libANALYSIS");
+       gSystem->Load("libANALYSISalice");
+       gSystem->Load("libSTEER");
+       gSystem->Load("libTENDER");
+       gSystem->Load("libRAWDatasim");
+       gSystem->Load("libFASTSIM");
+       gSystem->Load("libEVGEN");
+       gSystem->Load("libAliPythia6");
+       gSystem->Load("libSTAT");
+       gSystem->Load("libhijing");
+       gSystem->Load("libTHijing");
+       gSystem->Load("libSTRUCT");
+       gSystem->Load("libPHOSUtils");
+       gSystem->Load("libPHOSbase");
+       gSystem->Load("libPHOSsim");
+       gSystem->Load("libPHOSrec");
+       gSystem->Load("libMUONcore");
+       gSystem->Load("libMUONmapping");
+       gSystem->Load("libMUONgeometry");
+       gSystem->Load("libMUONcalib");
+       gSystem->Load("libMUONraw");
+       gSystem->Load("libMUONtrigger");
+       gSystem->Load("libMUONbase");
+       gSystem->Load("libMUONsim");
+       gSystem->Load("libMUONrec");
+       gSystem->Load("libMUONevaluation");
+       gSystem->Load("libFMDbase");
+       gSystem->Load("libFMDsim");
+       gSystem->Load("libFMDrec");
+       gSystem->Load("libPMDbase");
+       gSystem->Load("libPMDsim");
+       gSystem->Load("libPMDrec");
+       gSystem->Load("libHMPIDbase");
+       gSystem->Load("libHMPIDsim");
+       gSystem->Load("libHMPIDrec");
+       gSystem->Load("libT0base");
+       gSystem->Load("libT0sim");
+       gSystem->Load("libT0rec");
+       gSystem->Load("libZDCbase");
+       gSystem->Load("libZDCsim");
+       gSystem->Load("libZDCrec");
+       gSystem->Load("libACORDEbase");
+       gSystem->Load("libACORDErec");
+       gSystem->Load("libACORDEsim");
+       gSystem->Load("libVZERObase");
+       gSystem->Load("libVZEROrec");
+       gSystem->Load("libVZEROsim");
+       gSystem->Load("libEMCALraw");
+       gSystem->Load("libEMCALUtils");
+       gSystem->Load("libEMCALbase");
+       gSystem->Load("libEMCALsim");
+       gSystem->Load("libEMCALrec");
+       gSystem->Load("libTPCbase");
+       gSystem->Load("libTPCrec");
+       gSystem->Load("libTPCsim");
+       gSystem->Load("libTPCfast");
+       gSystem->Load("libITSbase");
+       gSystem->Load("libITSsim");
+       gSystem->Load("libITSrec");
+       gSystem->Load("libTRDbase");
+       gSystem->Load("libTRDsim");
+       gSystem->Load("libTRDrec");
+       gSystem->Load("libTOFbase");
+       gSystem->Load("libTOFrec");
+       gSystem->Load("libTOFsim");
+       gSystem->Load("libHLTbase");
+       gSystem->Load("libHLTinterface");
+       gSystem->Load("libHLTsim");
+       gSystem->Load("libHLTrec");
+       gSystem->Load("libPWGPP");
+
+       // add aliroot include path
+       gROOT->ProcessLine(Form(".include %s/include",
+                               gSystem->ExpandPathName("$ALICE_ROOT")));
+       gROOT->ProcessLine(Form(".include $ALICE_ROOT/include",
+                               gSystem->ExpandPathName("$ALICE_ROOT")));
+       gROOT->ProcessLine(Form(".include $ALICE_ROOT/ITS",
+                               gSystem->ExpandPathName("$ALICE_ROOT")));
+       gROOT->ProcessLine(Form(".include $ALICE_ROOT/PWGPP/ITS",
+                               gSystem->ExpandPathName("$ALICE_ROOT")));
+       gROOT->ProcessLine(Form(".include $ALICE_ROOT/VZERO",
+                               gSystem->ExpandPathName("$ALICE_ROOT")));
+
+       gROOT->SetStyle("Plain");
+
+       // create the alien handler and attach it to the manager
+       AliAnalysisGrid *plugin =
+               CreateAlienHandler(taskname, gridmode, proofcluster, proofdataset);
+
+       // analysis manager
+       AliAnalysisManager* mgr = new AliAnalysisManager("CDMeson-Manager");
+       mgr->SetGridHandler(plugin);
+
+       AliESDInputHandler* esdH = new AliESDInputHandler();
+       mgr->SetInputEventHandler(esdH);
+
+       // === Physics Selection Task ===
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+       AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(bMCphyssel);
+       if(!physSelTask) { Printf("no physSelTask"); return; }
+
+       // === Add PID Response Task ===
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+       AliAnalysisTaskPIDResponse *pidResponseTask = AddTaskPIDResponse(bMCphyssel);
+       if(!pidResponseTask) { Printf("no pidResponseTask"); return; }
+
+
+       // === create user task ===
+       gROOT->LoadMacro("AliCDMesonBase.cxx+g");
+       gROOT->LoadMacro("AliCDMesonTracks.cxx+g");
+       gROOT->LoadMacro("AliCDMesonUtils.cxx+g");
+       gROOT->LoadMacro("AliAnalysisTaskCDMeson.cxx+g");
+
+       Long_t taskConfig = AliCDMesonBase::kBitConfigurationSet; // has to be set
+       taskConfig |= AliCDMesonBase::kBitEtaPhiMaps;
+       taskConfig |= AliCDMesonBase::kBitEtaPhiMapsWithCuts;
+       taskConfig |= AliCDMesonBase::kBitStatsFlow;
+       taskConfig |= AliCDMesonBase::kBitMultPerGapHists;
+       taskConfig |= AliCDMesonBase::kBitRmMultPerGapHists;
+       taskConfig |= AliCDMesonBase::kBitTHnMother;
+       taskConfig |= AliCDMesonBase::kBitFastORStudy;
+       taskConfig |= AliCDMesonBase::kBitHitMapSPD;
+       taskConfig |= AliCDMesonBase::kBitHitMapFMD;
+       taskConfig |= AliCDMesonBase::kBitVtxStudies;
+       taskConfig |= AliCDMesonBase::kBitPIDStudy;
+       taskConfig |= AliCDMesonBase::kBitFMDsum;
+       taskConfig |= AliCDMesonBase::kBitSoftTracks;
+       taskConfig |= AliCDMesonBase::kBitPWAtree;
+       taskConfig |= AliCDMesonBase::kBitMultStudy;
+       //taskConfig |= AliCDMesonBase::kBitReadPreprocessedGap;
+       taskConfig |= AliCDMesonBase::kBitVZEROStudy;
+       taskConfig |= AliCDMesonBase::kBitTPCGapStudy;
+       if(bMCphyssel) taskConfig |= AliCDMesonBase::kBitMCProcess; // only for MC data
+       if(bMCphyssel) taskConfig |= AliCDMesonBase::kBitTHnMC;
+       if(bMCphyssel) taskConfig |= AliCDMesonBase::kBitMultResponseMC;
+       //taskConfig |= AliCDMesonBase::kBitReduceGapEvents;
+       taskConfig |= AliCDMesonBase::kBitConfigurationVersion;
+
+       printf("taskConfig=0x%x\n", taskConfig);
+
+       AliAnalysisTaskSE* task = new AliAnalysisTaskCDMeson(taskname, taskConfig);
+       task->SelectCollisionCandidates(AliVEvent::kMB);
+       mgr->AddTask(task);
+
+       // INPUT ---------------------------------------------------------------------
+       AliAnalysisDataContainer *cinput  = mgr->GetCommonInputContainer();
+       mgr->ConnectInput(task, 0, cinput );
+
+       // OUTPUT --------------------------------------------------------------------
+       // output filename
+       Char_t foutname[100];
+       Char_t foutnamePWA[100];
+       sprintf(foutname,"freidt_%s.root",taskname);
+       sprintf(foutnamePWA,"freidt_%sPWA.root",taskname);
+
+       // output containers
+       // in AnalysisTaskSE, slot 0 reserved, must start from 1
+       // since the number of output containers is flexible, the slot is assigned
+       // dynamically
+       AliAnalysisDataContainer* outlist[6] = { 0x0, 0x0, 0x0, 0x0, 0x0, 0x0 };
+       Int_t nOutputs = 0;
+
+       outlist[nOutputs] =
+               mgr->CreateContainer("freidt_Hist", TList::Class(),
+                                    AliAnalysisManager::kOutputContainer,foutname);
+       nOutputs++;
+       mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+
+       if (!taskConfig
+           || ((taskConfig & AliCDMesonBase::kBitSoftTracks)
+               && (taskConfig & AliCDMesonBase::kBitTHnMother))) {
+               outlist[nOutputs] =
+                       mgr->CreateContainer("freidt_ThnMother", THnSparse::Class(),
+                                            AliAnalysisManager::kOutputContainer,foutname);
+               nOutputs++;
+               mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+       }
+
+       if (!taskConfig || (taskConfig & AliCDMesonBase::kBitSoftTracks)) {
+               outlist[nOutputs] =
+                       mgr->CreateContainer("freidt_ThnMotherSoft", THnSparse::Class(),
+                                            AliAnalysisManager::kOutputContainer,foutname);
+               nOutputs++;
+               mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+       }
+
+       if (!taskConfig || (taskConfig & AliCDMesonBase::kBitMultStudy)) {
+               outlist[nOutputs] =
+                       mgr->CreateContainer("freidt_ThnMultiplicity", THnSparse::Class(),
+                                            AliAnalysisManager::kOutputContainer,foutname);
+               nOutputs++;
+               mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+       }
+
+       if (!taskConfig || (taskConfig & AliCDMesonBase::kBitTHnMC)) {
+               outlist[nOutputs] =
+                       mgr->CreateContainer("freidt_ThnMotherMC", THnSparse::Class(),
+                                            AliAnalysisManager::kOutputContainer,foutname);
+               nOutputs++;
+               mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+       }
+       if (!taskConfig || (taskConfig & AliCDMesonBase::kBitPWAtree)) {
+               outlist[nOutputs] =
+                       mgr->CreateContainer("freidt_PWA", TTree::Class(),
+                                            AliAnalysisManager::kOutputContainer,foutnamePWA);
+               nOutputs++;
+               mgr->ConnectOutput(task, nOutputs, outlist[nOutputs-1]);
+       }
+
+       // enable debug printouts
+       mgr->SetDebugLevel(2);
+       //mgr->SetNSysInfo(100);
+       if (!mgr->InitAnalysis()) return;
+       mgr->PrintStatus();
+
+       // start analysis
+       Printf("Starting Analysis....");
+       mgr->StartAnalysis(runtype); //,nentries,firstentry);
+}
+
+//______________________________________________________________________________
+AliAnalysisGrid* CreateAlienHandler(const char *taskname,
+                                    const char *gridmode,
+                                    const char *proofcluster,
+                                    const char *proofdataset)
+{
+       AliAnalysisAlien *plugin = new AliAnalysisAlien();
+       // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")
+       plugin->SetOverwriteMode();
+       plugin->SetRunMode(gridmode);
+
+       plugin->SetMergeViaJDL(kTRUE);
+
+       // Set versions of used packages
+       plugin->SetAPIVersion("V1.1x");
+       plugin->SetROOTVersion("v5-34-01-1");
+       plugin->SetAliROOTVersion("v5-03-50-AN");
+
+       // Declare input data to be processed.
+       //plugin->SetCheckCopy(kFALSE);
+
+       // Method 1: Create automatically XML collections using alien 'find' command.
+       // Define production directory LFN
+       plugin->SetGridDataDir("/alice/sim/2012/LHC12d4a");
+       // On real reconstructed data:
+       // plugin->SetGridDataDir("/alice/data/2009/LHC09d");
+       // Set data search pattern
+       plugin->SetDataPattern("*AliESDs.root"); // THIS CHOOSES ALL PASSES
+       // Data pattern for reconstructed data
+       //plugin->SetDataPattern("*ESDs/pass2/*ESDs.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH
+       //    plugin->SetDataPattern("ESDs/pass2/AOD038/*AliAOD.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH
+       plugin->SetRunPrefix(""); //000");   // real data
+       // ...then add run numbers to be considered
+
+       //Int_t runlist[15]={117039, 146859, 146858, 146856, 146824, 146817, 146806, 146805, 146804, 146803, 146802, 146801, 146748, 146747, 146746};  
+       //for (Int_t ind=0; ind<1; ind++) {
+       //      plugin->AddRunNumber(runlist[ind]);
+       //}
+       //plugin->SetRunRange(114917,115322);
+       //plugin->AddRunNumber(117050);
+       plugin->SetRunRange(115393,126408);
+
+       plugin->SetNrunsPerMaster(1);
+
+       // Define alien work directory where all files will be copied. Relative to alien $HOME.
+       plugin->SetGridWorkingDir(taskname);
+
+       // Declare alien output directory. Relative to working directory.
+       plugin->SetGridOutputDir("out"); // In this case will be $HOME/taskname/out
+
+       // Declare the analysis source files names separated by blancs. To be compiled runtime
+       // using ACLiC on the worker nodes.
+       plugin->SetAnalysisSource("AliCDMesonBase.cxx AliCDMesonTracks.cxx AliCDMesonUtils.cxx AliAnalysisTaskCDMeson.cxx");
+
+       // Declare all libraries (other than the default ones for the framework. These will be
+       // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.
+       plugin->SetAdditionalLibs("libGui.so libCore.so libTree.so libPhysics.so libMinuit.so libProof.so libmicrocern.so liblhapdf.so libpythia6.so libEG.so libGeom.so libVMC.so libEGPythia6.so libSTEERBase.so libESD.so libRAWDatabase.so libRAWDatarec.so libAOD.so libANALYSIS.so libANALYSISalice.so libCDB.so libSTEER.so libRAWDatasim.so libFASTSIM.so libEVGEN.so libAliPythia6.so libSTAT.so libhijing.so libTHijing.so libSTRUCT.so libPHOSUtils.so libPHOSbase.so libPHOSsim.so libPHOSrec.so libMUONcore.so libMUONmapping.so libMUONgeometry.so libMUONcalib.so libMUONraw.so libMUONtrigger.so libMUONbase.so libMUONsim.so libMUONrec.so libMUONevaluation.so libFMDbase.so libFMDsim.so libFMDrec.so libPMDbase.so libPMDsim.so libPMDrec.so libHMPIDbase.so libHMPIDsim.so libHMPIDrec.so libT0base.so libT0sim.so libT0rec.so libZDCbase.so libZDCsim.so libZDCrec.so libACORDEbase.so libACORDErec.so libACORDEsim.so libVZERObase.so libVZEROrec.so libVZEROsim.so libEMCALraw.so libEMCALUtils.so libEMCALbase.so libEMCALsim.so libEMCALrec.so libTPCbase.so libTPCrec.so libTPCsim.so libTPCfast.so libITSbase.so libITSsim.so libITSrec.so libTRDbase.so libTRDsim.so libTRDrec.so libTOFbase.so libTOFrec.so libTOFsim.so libHLTbase.so libHLTinterface.so libHLTsim.so libHLTrec.so libTENDER.so libTENDERSupplies.so libPWGPP.so AliCDMesonBase.h AliCDMesonBase.cxx AliCDMesonTracks.h AliCDMesonTracks.cxx AliCDMesonUtils.h AliCDMesonUtils.cxx AliAnalysisTaskCDMeson.h AliAnalysisTaskCDMeson.cxx");
+
+       plugin->AddIncludePath("-I$ALICE_ROOT/ITS -I$ALICE_ROOT/PWGPP/ITS");
+
+       // Declare the output file names separated by blancs.
+       // (can be like: file.root or file.root@ALICE::Niham::File)
+       // To only save certain files, use SetDefaultOutputs(kFALSE), and then
+       // SetOutputFiles("list.root other.filename") to choose which files to save
+       plugin->SetDefaultOutputs();
+
+       // Optionally set a name for the generated analysis macro (default MyAnalysis.C)
+       plugin->SetAnalysisMacro("CDMeson.C");
+
+       // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)
+       plugin->SetSplitMaxInputFileNumber(100);
+
+       // Optionally modify the executable name (default analysis.sh)
+       plugin->SetExecutable("CDMeson.sh");
+
+       // set number of test files to use in "test" mode
+       plugin->SetNtestFiles(1);
+
+       // Optionally resubmit threshold.
+       //plugin->SetMasterResubmitThreshold(90);
+
+       // Optionally set time to live (default 30000 sec)
+       //plugin->SetTTL(30000);
+
+       // Optionally set input format (default xml-single)
+       //plugin->SetInputFormat("xml-single");
+
+       // Optionally modify the name of the generated JDL (default analysis.jdl)
+       plugin->SetJDLName("CDMeson.jdl");
+
+       // Optionally modify job price (default 1)
+       plugin->SetPrice(1);
+
+       // Optionally modify split mode (default 'se')
+       plugin->SetSplitMode("se");
+
+       //----------------------------------------------------------
+       //---      PROOF MODE SPECIFIC SETTINGS         ------------
+       //----------------------------------------------------------
+       // Proof cluster
+       plugin->SetProofCluster(proofcluster);
+       // Dataset to be used
+       plugin->SetProofDataSet(proofdataset);
+       // May need to reset proof. Supported modes: 0-no reset, 1-soft, 2-hard
+       plugin->SetProofReset(0);
+       // May limit number of workers
+       plugin->SetNproofWorkers(0);
+       // May limit the number of workers per slave
+       plugin->SetNproofWorkersPerSlave(1);
+       // May use a specific version of root installed in proof
+       plugin->SetRootVersionForProof("current");
+       // May set the aliroot mode. Check http://aaf.cern.ch/node/83
+       plugin->SetAliRootMode("default"); // Loads AF libs by default
+       // May request ClearPackages (individual ClearPackage not supported)
+       plugin->SetClearPackages(kFALSE);
+       // Plugin test mode works only providing a file containing test file locations, used in "local" mode also
+       plugin->SetFileForTestMode("files.txt"); // file should contain path name to a local directory containg *ESDs.root etc
+       // Request connection to alien upon connection to grid
+       plugin->SetProofConnectGrid(kFALSE);
+       // Other PROOF specific parameters
+       plugin->SetProofParameter("PROOF_UseMergers","-1");
+       printf("Using: PROOF_UseMergers   : %s\n", plugin->GetProofParameter("PROOF_UseMergers"));
+
+       return plugin;
+}
index af4eb1f..7cbdc8a 100644 (file)
@@ -4,9 +4,14 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
+#pragma link C++ class AliAnalysisTaskCDMeson+;
+#pragma link C++ class AliAnalysisTaskCDskimESD+;
 #pragma link C++ class AliAnalysisTaskCDex+;
 #pragma link C++ class AliCDMesonBaseStripped+;
+#pragma link C++ class AliCDMesonBase+;
 #pragma link C++ class AliCDMesonTracks+;
 #pragma link C++ class AliCDMesonUtilsStripped+;
+#pragma link C++ class AliCDMesonUtils+;
+
 
 #endif