AliAnalysisTaskFastEmbedding:
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Oct 2011 18:42:39 +0000 (18:42 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 25 Oct 2011 18:42:39 +0000 (18:42 +0000)
- add event selection
- new random selection of extra AODs and events, taking into account the number of events in AODs
- additional mc information (pt hard, x-section, etc.)

AliAnalysisTaskJetResponseV2:
- add trigger exclude mask (e.g. jets with track pT > 100 GeV)
- add pt hard bin
- keep negative jet pT (after background subtraction)

JETAN/AliAnalysisTaskFastEmbedding.cxx
JETAN/AliAnalysisTaskFastEmbedding.h
PWG4/JetTasks/AliAnalysisTaskJetResponseV2.cxx
PWG4/JetTasks/AliAnalysisTaskJetResponseV2.h
PWG4/macros/AddTaskFastEmbedding.C
PWG4/macros/AddTaskJetResponseV2.C

index e131bc4..c092b06 100644 (file)
 #include <TRandom3.h>
 #include <TH1F.h>
 #include <TH2F.h>
+#include <TProfile.h>
+#include <TKey.h>
 
 
 #include "AliAnalysisTaskFastEmbedding.h"
 #include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliAODEvent.h"
 #include "AliAODTrack.h"
 #include "AliAODJet.h"
 #include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliInputEventHandler.h"
 
 #include "AliLog.h"
 
@@ -49,330 +54,509 @@ ClassImp(AliAnalysisTaskFastEmbedding)
 
 //__________________________________________________________________________
 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
-    : AliAnalysisTaskSE()
-      ,fAODout(0)
-      ,fAODevent(0)
-      ,fAODtree(0)
-      ,fAODfile(0)
-      ,rndm(0)
-      ,fAODPathArray(0)
-      ,fAODPath("AliAOD.root")
-      ,fTrackBranch("aodExtraTracks")
-      ,fMCparticlesBranch("aodExtraMCparticles")
-      ,fJetBranch("")
-      ,fEntry(0)
-      ,fEmbedMode(0)
-      ,fEvtSelecMode(0)
-      ,fEvtSelMinJetPt(-1)
-      ,fEvtSelMaxJetPt(-1)
-      ,fEvtSelMinJetEta(-999.)
-      ,fEvtSelMaxJetEta( 999.)
-      ,fEvtSelMinJetPhi(0.)
-      ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
-      ,fToyMinNbOfTracks(1)
-      ,fToyMaxNbOfTracks(1)
-      ,fToyMinTrackPt(50.)
-      ,fToyMaxTrackPt(50.)
-      ,fToyDistributionTrackPt(0.)
-      ,fToyMinTrackEta(-.5)
-      ,fToyMaxTrackEta(.5)
-      ,fToyMinTrackPhi(0.)
-      ,fToyMaxTrackPhi(2*TMath::Pi())
-      ,fToyFilterMap(0)
-      ,fHistList(0)
-      ,fh1TrackPt(0)
-      ,fh2TrackEtaPhi(0)
-      ,fh1TrackN(0)
-      ,fh1JetPt(0)
-      ,fh2JetEtaPhi(0)
-      ,fh1JetN(0)
-      ,fh1MCTrackPt(0)
-      ,fh2MCTrackEtaPhi(0)
-      ,fh1MCTrackN(0)
-      ,fh1AODfile(0)
-
-
+: AliAnalysisTaskSE()
+,fESD(0)
+,fAODout(0)
+,fAODevent(0)
+,fAODtree(0)
+,fAODfile(0)
+,mcHeader(0)
+,rndm(0)
+,fInputEntries(0)
+,fAODPathArray(0)
+,fAODEntriesArray(0)
+,fAODPath("AliAOD.root")
+,fAODEntries(-1)
+,fAODEntriesSum(0)
+,fAODEntriesMax(0)
+,fOfflineTrgMask(AliVEvent::kAny)
+,fMinContribVtx(1)
+,fVtxZMin(-8.)
+,fVtxZMax(8.)
+,fEvtClassMin(0)
+,fEvtClassMax(4)
+,fCentMin(0.)
+,fCentMax(100.)
+,fNInputTracksMin(0)
+,fNInputTracksMax(-1)
+,fTrackBranch("aodExtraTracks")
+,fMCparticlesBranch("aodExtraMCparticles")
+,fJetBranch("")
+,fFileId(-1)
+,fAODEntry(0)
+,fCountEvents(-1)
+,fEmbedMode(0)
+,fEvtSelecMode(0)
+,fEvtSelMinJetPt(-1)
+,fEvtSelMaxJetPt(-1)
+,fEvtSelMinJetEta(-999.)
+,fEvtSelMaxJetEta( 999.)
+,fEvtSelMinJetPhi(0.)
+,fEvtSelMaxJetPhi(TMath::Pi()*2.)
+,fToyMinNbOfTracks(1)
+,fToyMaxNbOfTracks(1)
+,fToyMinTrackPt(50.)
+,fToyMaxTrackPt(50.)
+,fToyDistributionTrackPt(0.)
+,fToyMinTrackEta(-.5)
+,fToyMaxTrackEta(.5)
+,fToyMinTrackPhi(0.)
+,fToyMaxTrackPhi(2*TMath::Pi())
+,fToyFilterMap(0)
+,fTrackFilterMap(0)
+,fNPtHard(10)
+,fPtHard(0)
+,fPtHardBin(-1)
+,fAODJets(0x0)
+,fNevents(0)
+,fXsection(0)
+,fAvgTrials(0)
+,fHistList(0)
+,fHistEvtSelection(0)
+,fh1Xsec(0)
+,fh1Trials(0)
+,fh1TrialsEvtSel(0)
+,fh2PtHard(0)
+,fh2PtHardEvtSel(0)
+,fh2PtHardTrials(0)
+,fh1TrackPt(0)
+,fh2TrackEtaPhi(0)
+,fh1TrackN(0)
+,fh1JetPt(0)
+,fh2JetEtaPhi(0)
+,fh1JetN(0)
+,fh1MCTrackPt(0)
+,fh2MCTrackEtaPhi(0)
+,fh1MCTrackN(0)
+,fh1AODfile(0)
+,fh2AODevent(0)
 {
-    // default constructor
+   // default constructor
 
 }
 
 
 //__________________________________________________________________________
 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
-    : AliAnalysisTaskSE(name)
-      ,fAODout(0)
-      ,fAODevent(0)
-      ,fAODtree(0)
-      ,fAODfile(0)
-      ,rndm(0)
-      ,fAODPathArray(0)
-      ,fAODPath("AliAOD.root")
-      ,fTrackBranch("aodExtraTracks")
-      ,fMCparticlesBranch("aodExtraMCparticles")
-      ,fJetBranch("")
-      ,fEntry(0)
-      ,fEmbedMode(0)
-      ,fEvtSelecMode(0)
-      ,fEvtSelMinJetPt(-1)
-      ,fEvtSelMaxJetPt(-1)
-      ,fEvtSelMinJetEta(-999.)
-      ,fEvtSelMaxJetEta( 999.)
-      ,fEvtSelMinJetPhi(0.)
-      ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
-      ,fToyMinNbOfTracks(1)
-      ,fToyMaxNbOfTracks(1)
-      ,fToyMinTrackPt(50.)
-      ,fToyMaxTrackPt(50.)
-      ,fToyDistributionTrackPt(0.)
-      ,fToyMinTrackEta(-.5)
-      ,fToyMaxTrackEta(.5)
-      ,fToyMinTrackPhi(0.)
-      ,fToyMaxTrackPhi(2*TMath::Pi())
-      ,fToyFilterMap(0)
-      ,fHistList(0)
-      ,fh1TrackPt(0)
-      ,fh2TrackEtaPhi(0)
-      ,fh1TrackN(0)
-      ,fh1JetPt(0)
-      ,fh2JetEtaPhi(0)
-      ,fh1JetN(0)
-      ,fh1MCTrackPt(0)
-      ,fh2MCTrackEtaPhi(0)
-      ,fh1MCTrackN(0)
-      ,fh1AODfile(0)
+: AliAnalysisTaskSE(name)
+,fESD(0)
+,fAODout(0)
+,fAODevent(0)
+,fAODtree(0)
+,fAODfile(0)
+,mcHeader(0)
+,rndm(0)
+,fInputEntries(0)
+,fAODPathArray(0)
+,fAODEntriesArray(0)
+,fAODPath("AliAOD.root")
+,fAODEntries(-1)
+,fAODEntriesSum(0)
+,fAODEntriesMax(0)
+,fOfflineTrgMask(AliVEvent::kAny)
+,fMinContribVtx(1)
+,fVtxZMin(-8.)
+,fVtxZMax(8.)
+,fEvtClassMin(0)
+,fEvtClassMax(4)
+,fCentMin(0.)
+,fCentMax(100.)
+,fNInputTracksMin(0)
+,fNInputTracksMax(-1)
+,fTrackBranch("aodExtraTracks")
+,fMCparticlesBranch("aodExtraMCparticles")
+,fJetBranch("")
+,fFileId(-1)
+,fAODEntry(0)
+,fCountEvents(-1)
+,fEmbedMode(0)
+,fEvtSelecMode(0)
+,fEvtSelMinJetPt(-1)
+,fEvtSelMaxJetPt(-1)
+,fEvtSelMinJetEta(-999.)
+,fEvtSelMaxJetEta( 999.)
+,fEvtSelMinJetPhi(0.)
+,fEvtSelMaxJetPhi(TMath::Pi()*2.)
+,fToyMinNbOfTracks(1)
+,fToyMaxNbOfTracks(1)
+,fToyMinTrackPt(50.)
+,fToyMaxTrackPt(50.)
+,fToyDistributionTrackPt(0.)
+,fToyMinTrackEta(-.5)
+,fToyMaxTrackEta(.5)
+,fToyMinTrackPhi(0.)
+,fToyMaxTrackPhi(2*TMath::Pi())
+,fToyFilterMap(0)
+,fTrackFilterMap(0)
+,fNPtHard(10)
+,fPtHard(0)
+,fPtHardBin(-1)
+,fAODJets(0x0)
+,fNevents(0)
+,fXsection(0)
+,fAvgTrials(0)
+,fHistList(0)
+,fHistEvtSelection(0)
+,fh1Xsec(0)
+,fh1Trials(0)
+,fh1TrialsEvtSel(0)
+,fh2PtHard(0)
+,fh2PtHardEvtSel(0)
+,fh2PtHardTrials(0)
+,fh1TrackPt(0)
+,fh2TrackEtaPhi(0)
+,fh1TrackN(0)
+,fh1JetPt(0)
+,fh2JetEtaPhi(0)
+,fh1JetN(0)
+,fh1MCTrackPt(0)
+,fh2MCTrackEtaPhi(0)
+,fh1MCTrackN(0)
+,fh1AODfile(0)
+,fh2AODevent(0)
 {
-    // constructor
-    DefineOutput(1, TList::Class());
+   // constructor
+   DefineOutput(1, TList::Class());
 }
 
 
 //__________________________________________________________________________
 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy)
-    : AliAnalysisTaskSE()
-      ,fAODout(copy.fAODout)
-      ,fAODevent(copy.fAODevent)
-      ,fAODtree(copy.fAODtree)
-      ,fAODfile(copy.fAODfile)
-      ,rndm(copy.rndm)
-      ,fAODPathArray(copy.fAODPathArray)
-      ,fAODPath(copy.fAODPath)
-      ,fTrackBranch(copy.fTrackBranch)
-      ,fMCparticlesBranch(copy.fMCparticlesBranch)
-      ,fJetBranch(copy.fJetBranch)
-      ,fEntry(copy.fEntry)
-      ,fEmbedMode(copy.fEmbedMode)
-      ,fEvtSelecMode(copy.fEvtSelecMode)
-      ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
-      ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
-      ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
-      ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
-      ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
-      ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
-      ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
-      ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
-      ,fToyMinTrackPt(copy.fToyMinTrackPt)
-      ,fToyMaxTrackPt(copy.fToyMaxTrackPt)
-      ,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
-      ,fToyMinTrackEta(copy.fToyMinTrackEta)
-      ,fToyMaxTrackEta(copy.fToyMaxTrackEta)
-      ,fToyMinTrackPhi(copy.fToyMinTrackPhi)
-      ,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
-      ,fToyFilterMap(copy.fToyFilterMap)
-      ,fHistList(copy.fHistList)
-      ,fh1TrackPt(copy.fh1TrackPt)
-      ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
-      ,fh1TrackN(copy.fh1TrackN)
-      ,fh1JetPt(copy.fh1JetPt)
-      ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
-      ,fh1JetN(copy.fh1JetN)
-      ,fh1MCTrackPt(copy.fh1MCTrackPt)
-      ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
-      ,fh1MCTrackN(copy.fh1MCTrackN)
-      ,fh1AODfile(copy.fh1AODfile)
+: AliAnalysisTaskSE()
+,fESD(copy.fESD)
+,fAODout(copy.fAODout)
+,fAODevent(copy.fAODevent)
+,fAODtree(copy.fAODtree)
+,fAODfile(copy.fAODfile)
+,mcHeader(copy.mcHeader)
+,rndm(copy.rndm)
+,fInputEntries(copy.fInputEntries)
+,fAODPathArray(copy.fAODPathArray)
+,fAODEntriesArray(copy.fAODEntriesArray)
+,fAODPath(copy.fAODPath)
+,fAODEntries(copy.fAODEntries)
+,fAODEntriesSum(copy.fAODEntriesSum)
+,fAODEntriesMax(copy.fAODEntriesMax)
+,fOfflineTrgMask(copy.fOfflineTrgMask)
+,fMinContribVtx(copy.fMinContribVtx)
+,fVtxZMin(copy.fVtxZMin)
+,fVtxZMax(copy.fVtxZMax)
+,fEvtClassMin(copy.fEvtClassMin)
+,fEvtClassMax(copy.fEvtClassMax)
+,fCentMin(copy.fCentMin)
+,fCentMax(copy.fCentMax)
+,fNInputTracksMin(copy.fNInputTracksMin)
+,fNInputTracksMax(copy.fNInputTracksMax)
+,fTrackBranch(copy.fTrackBranch)
+,fMCparticlesBranch(copy.fMCparticlesBranch)
+,fJetBranch(copy.fJetBranch)
+,fFileId(copy.fFileId)
+,fAODEntry(copy.fAODEntry)
+,fCountEvents(copy.fCountEvents)
+,fEmbedMode(copy.fEmbedMode)
+,fEvtSelecMode(copy.fEvtSelecMode)
+,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
+,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
+,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
+,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
+,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
+,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
+,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
+,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
+,fToyMinTrackPt(copy.fToyMinTrackPt)
+,fToyMaxTrackPt(copy.fToyMaxTrackPt)
+,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
+,fToyMinTrackEta(copy.fToyMinTrackEta)
+,fToyMaxTrackEta(copy.fToyMaxTrackEta)
+,fToyMinTrackPhi(copy.fToyMinTrackPhi)
+,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
+,fToyFilterMap(copy.fToyFilterMap)
+,fTrackFilterMap(copy.fTrackFilterMap)
+,fNPtHard(copy.fNPtHard)
+,fPtHard(copy.fPtHard)
+,fPtHardBin(copy.fPtHardBin)
+,fAODJets(copy.fAODJets)
+,fNevents(copy.fNevents)
+,fXsection(copy.fXsection)
+,fAvgTrials(copy.fAvgTrials)
+,fHistList(copy.fHistList)
+,fHistEvtSelection(copy.fHistEvtSelection)
+,fh1Xsec(copy.fh1Xsec)
+,fh1Trials(copy.fh1Trials)
+,fh1TrialsEvtSel(copy.fh1TrialsEvtSel)
+,fh2PtHard(copy.fh2PtHard)
+,fh2PtHardEvtSel(copy.fh2PtHardEvtSel)
+,fh2PtHardTrials(copy.fh2PtHardTrials)
+,fh1TrackPt(copy.fh1TrackPt)
+,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
+,fh1TrackN(copy.fh1TrackN)
+,fh1JetPt(copy.fh1JetPt)
+,fh2JetEtaPhi(copy.fh2JetEtaPhi)
+,fh1JetN(copy.fh1JetN)
+,fh1MCTrackPt(copy.fh1MCTrackPt)
+,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
+,fh1MCTrackN(copy.fh1MCTrackN)
+,fh1AODfile(copy.fh1AODfile)
+,fh2AODevent(copy.fh2AODevent)
 {
-    // copy constructor
+   // copy constructor
 }
 
 
 //__________________________________________________________________________
 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
 {
-    // assignment
-
-    if(this!=&o){
-       AliAnalysisTaskSE::operator=(o);
-       fAODout            = o.fAODout;
-       fAODevent          = o.fAODevent;
-       fAODtree           = o.fAODtree;
-       fAODfile           = o.fAODfile;
-        rndm               = o.rndm;
-        fAODPathArray       = o.fAODPathArray;
-        fAODPath           = o.fAODPath;
-        fTrackBranch       = o.fTrackBranch;
-        fMCparticlesBranch = o.fMCparticlesBranch;
-        fJetBranch         = o.fJetBranch;
-        fEntry             = o.fEntry;
-        fEmbedMode         = o.fEmbedMode;
-        fEvtSelecMode      = o.fEvtSelecMode;
-        fEvtSelMinJetPt    = o.fEvtSelMinJetPt;
-        fEvtSelMaxJetPt    = o.fEvtSelMaxJetPt;
-        fEvtSelMinJetEta   = o.fEvtSelMinJetEta;
-        fEvtSelMaxJetEta   = o.fEvtSelMaxJetEta;
-        fEvtSelMinJetPhi   = o.fEvtSelMinJetPhi;
-        fEvtSelMaxJetPhi   = o.fEvtSelMaxJetPhi;
-        fToyMinNbOfTracks  = o.fToyMinNbOfTracks;
-        fToyMaxNbOfTracks  = o.fToyMaxNbOfTracks;
-        fToyMinTrackPt     = o.fToyMinTrackPt;
-        fToyMaxTrackPt     = o.fToyMaxTrackPt;
-        fToyDistributionTrackPt = o.fToyDistributionTrackPt;
-        fToyMinTrackEta    = o.fToyMinTrackEta;
-        fToyMaxTrackEta    = o.fToyMaxTrackEta;
-        fToyMinTrackPhi    = o.fToyMinTrackPhi;
-        fToyMaxTrackPhi    = o.fToyMaxTrackPhi;
-        fToyFilterMap      = o.fToyFilterMap;
-        fHistList          = o.fHistList;
-        fh1TrackPt         = o.fh1TrackPt;
-        fh2TrackEtaPhi     = o.fh2TrackEtaPhi;
-        fh1TrackN          = o.fh1TrackN;
-       fh1JetPt           = o.fh1JetPt;
-        fh2JetEtaPhi       = o.fh2JetEtaPhi;
-        fh1JetN            = o.fh1JetN;
-        fh1MCTrackPt       = o.fh1MCTrackPt;
-        fh2MCTrackEtaPhi   = o.fh2MCTrackEtaPhi;
-        fh1MCTrackN        = o.fh1MCTrackN;
-        fh1AODfile         = o.fh1AODfile;
-    }
-
-    return *this;
+   // assignment
+
+   if(this!=&o){
+      AliAnalysisTaskSE::operator=(o);
+      fESD               = o.fESD;
+      fAODout            = o.fAODout;
+      fAODevent          = o.fAODevent;
+      fAODtree           = o.fAODtree;
+      fAODfile           = o.fAODfile;
+      mcHeader           = o.mcHeader;
+      rndm               = o.rndm;
+      fInputEntries      = o.fInputEntries;
+      fAODPathArray      = o.fAODPathArray;
+      fAODEntriesArray   = o.fAODEntriesArray;
+      fAODPath           = o.fAODPath;
+      fAODEntries        = o.fAODEntries;
+      fAODEntriesSum     = o.fAODEntriesSum;
+      fAODEntriesMax     = o.fAODEntriesMax;
+      fOfflineTrgMask    = o.fOfflineTrgMask;
+      fMinContribVtx     = o.fMinContribVtx;
+      fVtxZMin           = o.fVtxZMin;
+      fVtxZMax           = o.fVtxZMax;
+      fEvtClassMin       = o.fEvtClassMin;
+      fEvtClassMax       = o.fEvtClassMax;
+      fCentMin           = o.fCentMin;
+      fCentMax           = o.fCentMax;
+      fNInputTracksMin   = o.fNInputTracksMin;
+      fNInputTracksMax   = o.fNInputTracksMax;
+      fTrackBranch       = o.fTrackBranch;
+      fMCparticlesBranch = o.fMCparticlesBranch;
+      fJetBranch         = o.fJetBranch;
+      fFileId            = o.fFileId;
+      fAODEntry          = o.fAODEntry;
+      fCountEvents       = o.fCountEvents;
+      fEmbedMode         = o.fEmbedMode;
+      fEvtSelecMode      = o.fEvtSelecMode;
+      fEvtSelMinJetPt    = o.fEvtSelMinJetPt;
+      fEvtSelMaxJetPt    = o.fEvtSelMaxJetPt;
+      fEvtSelMinJetEta   = o.fEvtSelMinJetEta;
+      fEvtSelMaxJetEta   = o.fEvtSelMaxJetEta;
+      fEvtSelMinJetPhi   = o.fEvtSelMinJetPhi;
+      fEvtSelMaxJetPhi   = o.fEvtSelMaxJetPhi;
+      fToyMinNbOfTracks  = o.fToyMinNbOfTracks;
+      fToyMaxNbOfTracks  = o.fToyMaxNbOfTracks;
+      fToyMinTrackPt     = o.fToyMinTrackPt;
+      fToyMaxTrackPt     = o.fToyMaxTrackPt;
+      fToyDistributionTrackPt = o.fToyDistributionTrackPt;
+      fToyMinTrackEta    = o.fToyMinTrackEta;
+      fToyMaxTrackEta    = o.fToyMaxTrackEta;
+      fToyMinTrackPhi    = o.fToyMinTrackPhi;
+      fToyMaxTrackPhi    = o.fToyMaxTrackPhi;
+      fToyFilterMap      = o.fToyFilterMap;
+      fTrackFilterMap    = o.fTrackFilterMap;
+      fNPtHard           = o.fNPtHard;
+      fPtHard            = o.fPtHard;
+      fPtHardBin         = o.fPtHardBin;
+      fAODJets           = o.fAODJets;
+      fNevents           = o.fNevents;
+      fXsection          = o.fXsection;
+      fAvgTrials         = o.fAvgTrials;
+      fHistList          = o.fHistList;
+      fHistEvtSelection  = o.fHistEvtSelection;
+      fh1Xsec            = o.fh1Xsec;
+      fh1Trials          = o.fh1Trials;
+      fh1TrialsEvtSel    = o.fh1TrialsEvtSel;
+      fh2PtHard          = o.fh2PtHard;
+      fh2PtHardEvtSel    = o.fh2PtHardEvtSel;
+      fh2PtHardTrials    = o.fh2PtHardTrials;
+      fh1TrackPt         = o.fh1TrackPt;
+      fh2TrackEtaPhi     = o.fh2TrackEtaPhi;
+      fh1TrackN          = o.fh1TrackN;
+      fh1JetPt           = o.fh1JetPt;
+      fh2JetEtaPhi       = o.fh2JetEtaPhi;
+      fh1JetN            = o.fh1JetN;
+      fh1MCTrackPt       = o.fh1MCTrackPt;
+      fh2MCTrackEtaPhi   = o.fh2MCTrackEtaPhi;
+      fh1MCTrackN        = o.fh1MCTrackN;
+      fh1AODfile         = o.fh1AODfile;
+      fh2AODevent        = o.fh2AODevent;
+   }
+
+   return *this;
 }
 
 
 //__________________________________________________________________________
 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
 {
-    // destructor
-    delete rndm;
+   // destructor
+   delete rndm;
 }
 
 
 //__________________________________________________________________________
 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
 {
-    // create output objects
-    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
-
-       OpenFile(1);
-    if(!fHistList) fHistList = new TList();
-    fHistList->SetOwner(kTRUE);
-       
-
-    // set seed
-    rndm = new TRandom3();
-    Int_t id = GetJobID();
-    if(id>-1) rndm->SetSeed(id);
-    else      rndm->SetSeed();   // a TTUID is generated and used for seed
-    AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
-
-    // embed mode with AOD
-    Int_t rc = -1;
-    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-
-       // open input AOD
-       rc = OpenAODfile();
-       if(rc<0){
-           PostData(1, fHistList);
-          return;
-          }
-    } //end: embed mode with AOD
-
-
-    // connect output aod
-    // create a new branch for extra tracks
-    fAODout = AODEvent();
-    if(!fAODout){
-        AliError("Output AOD not found.");
-               PostData(1, fHistList);
-        return;
-    }
-    if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
-        AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
-        TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
-        tracks->SetName(fTrackBranch.Data());
-        AddAODBranch("TClonesArray", &tracks);
-    }
-    // create new branch for extra mcparticle if available as input
-    if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
-       AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
-       TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
-       mcparticles->SetName(fMCparticlesBranch.Data());
-       AddAODBranch("TClonesArray", &mcparticles);
-    }
-
-
-
-    //qa histograms
-    Bool_t oldStatus = TH1::AddDirectoryStatus();
-    TH1::AddDirectory(kFALSE);
-
-    fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
-    fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
-    fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
-
-    fHistList->Add(fh1TrackPt);
-    fHistList->Add(fh2TrackEtaPhi);
-    fHistList->Add(fh1TrackN);
-       
-       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-       
-           fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 250, 0., 250.);
-           fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
-           fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
-       
-           fHistList->Add(fh1JetPt);
-        fHistList->Add(fh2JetEtaPhi);
-        fHistList->Add(fh1JetN);
-    }
-
-    if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){ 
-
-       fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
-       fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
-       fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
-          
-       fHistList->Add(fh1MCTrackPt);
-       fHistList->Add(fh2MCTrackEtaPhi);
-       fHistList->Add(fh1MCTrackN);
-         
-    }
-       
-    fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 100, 0, 100);
-    fHistList->Add(fh1AODfile);
-
-    // =========== Switch on Sumw2 for all histos ===========
-    for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
-        TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
-        if (h1){
-          h1->Sumw2();
-          continue;
-        }
-    }
-
-    TH1::AddDirectory(oldStatus);
-
-
-    fh1AODfile->Fill(rc);
-    
-    PostData(1, fHistList);
+   // create output objects
+   if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
+   AliLog::SetClassDebugLevel("AliAnalysisTaskFastEmbedding", AliLog::kInfo);
+
+   OpenFile(1);
+   if(!fHistList) fHistList = new TList();
+   fHistList->SetOwner(kTRUE);
+   
+
+   // set seed
+   rndm = new TRandom3();
+   Int_t id = GetJobID();
+   if(id>-1) rndm->SetSeed(id);
+   else      rndm->SetSeed();   // a TTUID is generated and used for seed
+   AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
+
+
+
+   // embed mode with AOD
+   if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+
+      // open input AOD
+      fFileId = OpenAODfile();
+      fNevents = 0; // force to open another aod in UserExec()
+      if(fFileId<0){
+         AliError("");
+         PostData(1, fHistList);
+         return;
+      }
+   } //end: embed mode with AOD
+
+
+   // connect output aod
+   // create a new branch for extra tracks
+   fAODout = AODEvent();
+   if(!fAODout){
+      AliError("Output AOD not found.");
+      PostData(1, fHistList);
+      return;
+   }
+   if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
+      AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
+      TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
+      tracks->SetName(fTrackBranch.Data());
+      AddAODBranch("TClonesArray", &tracks);
+   }
+   // create new branch for extra mcparticle if available as input
+   if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
+      AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
+      TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
+      mcparticles->SetName(fMCparticlesBranch.Data());
+      AddAODBranch("TClonesArray", &mcparticles);
+   }
+
+
+
+   //qa histograms
+   Bool_t oldStatus = TH1::AddDirectoryStatus();
+   TH1::AddDirectory(kFALSE);
+   
+   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
+   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
+   
+   fh1Xsec         = new TProfile("fh1Xsec","xsec from pyxsec.root;p_{T,hard} bin;<#sigma>",fNPtHard+1,-1.5,fNPtHard-0.5);
+   fh1Trials       = new TH1F("fh1Trials","trials (simulation) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
+   fh1TrialsEvtSel = new TH1F("fh1TrialsEvtSel","trials (event selection) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
+   fh2PtHard       = new TH2F("fh2PtHard","PYTHIA Pt hard;p_{T,hard} bin;p_{T,hard}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
+   fh2PtHardEvtSel = new TH2F("fh2PtHardEvtSel","PYTHIA Pt hard (event selection);p_{T,hard} bin;p_{T,hard}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
+   fh2PtHardTrials = new TH2F("fh2PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard} bin;#sum{p_{T,hard}}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
+   
+   fHistList->Add(fHistEvtSelection);
+   fHistList->Add(fh1Xsec);
+   fHistList->Add(fh1Trials);
+   fHistList->Add(fh1TrialsEvtSel);
+   fHistList->Add(fh2PtHard);
+   fHistList->Add(fh2PtHardEvtSel);
+   fHistList->Add(fh2PtHardTrials);
+
+   fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
+   fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
+   fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
+
+   fHistList->Add(fh1TrackPt);
+   fHistList->Add(fh2TrackEtaPhi);
+   fHistList->Add(fh1TrackN);
+   
+   if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+      
+      fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 120, 0., 120.);
+      fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
+      fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
+      
+      fHistList->Add(fh1JetPt);
+      fHistList->Add(fh2JetEtaPhi);
+      fHistList->Add(fh1JetN);
+   }
+
+
+   if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){ 
+
+      fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
+      fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
+      fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
+      
+      fHistList->Add(fh1MCTrackPt);
+      fHistList->Add(fh2MCTrackEtaPhi);
+      fHistList->Add(fh1MCTrackN);
+      
+   }
+   
+   fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 2300, -0.5, 2299.5);
+   fh2AODevent = new TH2I("fh1AODevent","selected events;file;event", 2500,-.5,2499.5,5000,-.5,4999.5);
+   fHistList->Add(fh1AODfile);
+   fHistList->Add(fh2AODevent);
+
+   // =========== Switch on Sumw2 for all histos ===========
+   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
+      TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
+      if (h1){
+         h1->Sumw2();
+         continue;
+      }
+   }
+
+   TH1::AddDirectory(oldStatus);
+
+   PostData(1, fHistList);
 }
 
 
 //__________________________________________________________________________
 void AliAnalysisTaskFastEmbedding::Init()
 {
-    // Initialization
-    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
+   // Initialization
+   if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
+
+}
+
+//__________________________________________________________________________
+Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
+{
+
+   if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserNotify()");
+
+   // get total nb of events in tree (of this subjob)
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+   fInputEntries = inputHandler->GetTree()->GetEntriesFast();
+   AliInfo(Form("Total nb. of events: %d", fInputEntries));
+
+   return kTRUE;
 
 }
 
@@ -380,322 +564,408 @@ void AliAnalysisTaskFastEmbedding::Init()
 //__________________________________________________________________________
 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
 {
-    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
-
-    if(!fAODout){
-       AliError("Need output AOD, but is not connected."); 
-       PostData(1, fHistList);
-       return;
-    }
-
-    // connect aod out
-    TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
-    if(!tracks){
-        AliError("Extra track branch not found in output.");
-        PostData(1, fHistList);
-        return;
-    }
-    tracks->Delete();
-    Int_t nAODtracks=0;
-
-
-    TRef dummy;
-
-    // === embed mode with AOD ===
-    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-       if(!fAODevent){
-          AliError("Need input AOD, but is not connected."); 
-          PostData(1, fHistList);
-          return;
-       }
-
-       // fetch jets
-       TClonesArray *aodJets = 0;
-       if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
-       else                    aodJets = fAODevent->GetJets();
-       if(!aodJets){
-          AliError("Could not find jets in AOD. Check jet branch when indicated.");
-          PostData(1, fHistList);
-          return;
-       }
-        
-
-       Int_t nEvents = fAODtree->GetEntries();
-
-       Bool_t useEntry = kFALSE;
-       while(!useEntry){  // protection need, if no event fulfills requierment
-          if(fEntry>=nEvents){
-              fEntry=0;
-              if(!fAODPathArray){
-                 AliDebug(AliLog::kDebug, "Last event in AOD reached, start from entry 0 again.");
-              } 
-              else {
-                 AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
-
-                 Int_t rc = OpenAODfile();
-                 if(rc<0) {
-                    PostData(1, fHistList);
-                    return;
-                 }
-                 fh1AODfile->Fill(rc);
-
-                // new file => we must use the new jet array
-                if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
-                else                    aodJets = fAODevent->GetJets();
-                if(!aodJets){
-                  AliError("Could not find jets in AOD. Check jet branch when indicated.");
-                   PostData(1, fHistList);
-                  return;
-                }
-              }
-          }
-    
-          fAODtree->GetEvent(fEntry);
-
-         // jet pt selection
-         if(fEvtSelecMode==kEventsJetPt){
-             Int_t nJets = aodJets->GetEntries();
-             for(Int_t ij=0; ij<nJets; ++ij){
-                 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
-                 if(!jet) continue;
-
-                 if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
-                   && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
-                    && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
-                    && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
-                    useEntry = kTRUE;
-                    break;
-                 } 
-             }
-         }
-
-          // no selection
-         if(fEvtSelecMode==kEventsAll){
-             useEntry = kTRUE;
-         }
-
-         fEntry++;
-       }
-       AliInfo(Form("Use entry %d from extra AOD.", fEntry-1));
-
-
-       TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
-       TClonesArray *mcpartOUT = 0x0;
-       if(mcpartIN){
-          mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
-       if(!mcpartOUT){
-         AliError("Extra MC particles branch not found in output.");
+   if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
+
+   if(!fAODout){
+      AliError("Need output AOD, but is not connected."); 
+      PostData(1, fHistList);
+      return;
+   }
+
+   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
+   if (!fESD) {
+      AliError("ESD not available");
+      PostData(1, fHistList);
+      return;
+   }
+
+   // -- event selection --
+   fHistEvtSelection->Fill(1); // number of events before event selection
+
+   // physics selection
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+   ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
+      if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
+      fHistEvtSelection->Fill(2);
+      PostData(1, fHistList);
+      return;
+   } 
+
+   // vertex selection
+   AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
+   Int_t nTracksPrim = primVtx->GetNContributors();
+   if ((nTracksPrim < fMinContribVtx) ||
+         (primVtx->GetZ() < fVtxZMin) ||
+         (primVtx->GetZ() > fVtxZMax) ){
+      if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+      fHistEvtSelection->Fill(3);
+      PostData(1, fHistList);
+      return;
+   }
+
+   /** takes wrong value, since jet service tasks runs later
+   // event class selection (from jet helper task)
+   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
+   if(fDebug) Printf("Event class %d", eventClass);
+   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
+   fHistEvtSelection->Fill(4);
+   PostData(1, fHistList);
+   return;
+   }*/
+
+   // centrality selection
+   AliCentrality *cent = 0x0;
+   Float_t centValue = 0.; 
+   if(fESD) cent = fESD->GetCentrality();
+   if(cent) centValue = cent->GetCentralityPercentile("V0M");
+   if(fDebug) printf("centrality: %f\n", centValue);
+   if (centValue < fCentMin || centValue > fCentMax){
+      fHistEvtSelection->Fill(4);
+      PostData(1, fHistList);
+      return;
+   }
+
+
+   /*  ** not implemented **
+   // multiplicity due to input tracks
+   Int_t nInputTracks = GetNInputTracks();
+
+   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
+   fHistEvtSelection->Fill(5);
+   PostData(1, fHistList);
+   return;
+   }
+   */
+
+   fHistEvtSelection->Fill(0); // accepted events  
+   // -- end event selection --
+
+   // connect aod out
+   TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
+   if(!tracks){
+      AliError("Extra track branch not found in output.");
+      PostData(1, fHistList);
+      return;
+   }
+   tracks->Delete();
+   Int_t nAODtracks=0;
+
+   TRef dummy;
+
+   // === embed mode with AOD ===
+   if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+      if(!fAODevent){
+         AliError("Need input AOD, but is not connected."); 
          PostData(1, fHistList);
          return;
-       }
-          mcpartOUT->Delete();
-       } else {
-          AliInfo("No extra MC particles found.");
-       }
-    
-
-       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
-          // loop over input tracks
-          // add to output aod
-          Int_t nTracks = 0;
-          Int_t nJets = aodJets->GetEntries();
-          Int_t nSelectedJets = 0;
-       AliAODJet *leadJet = 0x0; // used for jet tracks only
-           
-           if(fEmbedMode==kAODFull){
-                          nTracks = fAODevent->GetNumberOfTracks();
-                                  
-                              for(Int_t ij=0; ij<nJets; ++ij){
-                       AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
-                       if(!jet) continue;
-                       if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
-                          && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
-                          && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
-                          && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
-                                                 
-                                                     fh1JetPt->Fill(jet->Pt());
-                                                         fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
-                                                         nSelectedJets++;
-                                                         
-                                          }
-                   }                              
-                                  fh1JetN->Fill(nSelectedJets);
-                  }
-
-           if(fEmbedMode==kAODJetTracks){
-              // look for leading jet within selection
-              // get reference tracks for this jet
-              Float_t leadJetPt = 0.;
-              for(Int_t ij=0; ij<nJets; ++ij){
-                  AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
-                  if(!jet) continue;
-                  if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+      }
+
+
+
+      Bool_t useEntry = kFALSE;
+      while(!useEntry){  // protection need, if no event fulfills requierment
+
+         fAODEntry++; // go to next event 
+         fCountEvents++;
+         if(fAODEntry>=fNevents) fAODEntry=0; 
+
+         // open new aod file
+         if(fCountEvents>=fNevents){ 
+            if(!fAODPathArray){
+               AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
+            } 
+            else {
+               AliDebug(AliLog::kDebug, "Select new AOD file ...");
+
+               fFileId = OpenAODfile();
+               if(fFileId<0) {
+                  PostData(1, fHistList);
+                  return;
+               }
+               fh1AODfile->Fill(fFileId);
+               if(fAODEntries!=fNevents){
+                  AliError("File with list of AODs and file with nb. of entries inconsistent");
+                  PostData(1, fHistList);
+                  return;
+               }
+            }
+         }
+
+
+
+         // get next event
+         fAODtree->GetEvent(fAODEntry);
+
+         // get pt hard
+         if(mcHeader){
+            fPtHard = mcHeader->GetPtHard();
+            GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
+            AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
+         }
+         else{
+            AliDebug(AliLog::kDebug,"No mcHeader");
+         }
+         fPtHardBin = GetPtHardBin(fPtHard);
+
+         //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials); 
+
+         // fill event variables for each event
+         fh1Xsec->Fill(fPtHardBin,fXsection);
+         fh2PtHard->Fill(fPtHardBin,fPtHard);
+
+         fh1Trials->Fill(fPtHardBin, fAvgTrials);
+         fh1TrialsEvtSel->Fill(fPtHardBin);
+         fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
+
+         // jet pt selection
+         if(fEvtSelecMode==kEventsJetPt){
+            Int_t nJets = fAODJets->GetEntries();
+            //Printf("nb. jets: %d",nJets);
+            for(Int_t ij=0; ij<nJets; ++ij){
+               AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
+               if(!jet) continue;
+
+               if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+                     && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+                     && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+                     && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
+                  useEntry = kTRUE;
+                  break;
+               } 
+            }
+         }
+
+         // no selection
+         if(fEvtSelecMode==kEventsAll){
+            useEntry = kTRUE;
+         }
+
+      }
+      AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
+
+      fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
+      fh2AODevent->Fill(fFileId,fAODEntry);
+
+      TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
+      TClonesArray *mcpartOUT = 0x0;
+      if(mcpartIN){
+         mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
+         if(!mcpartOUT){
+            AliError("Extra MC particles branch not found in output.");
+            PostData(1, fHistList);
+            return;
+         }
+         mcpartOUT->Delete();
+      } else {
+         AliInfo("No extra MC particles found.");
+      }
+
+
+      if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
+         // loop over input tracks
+         // add to output aod
+         Int_t nTracks = 0;
+         Int_t nJets = fAODJets->GetEntries();
+         Int_t nSelectedJets = 0;
+         AliAODJet *leadJet = 0x0; // used for jet tracks only
+
+         if(fEmbedMode==kAODFull){
+            nTracks = fAODevent->GetNumberOfTracks();
+
+            for(Int_t ij=0; ij<nJets; ++ij){
+               AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
+               if(!jet) continue;
+               if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
-                     if(jet->Pt()>leadJetPt){
-                         leadJetPt = jet->Pt();
-                         leadJet = jet;
-                     } 
-                  }
-               }
-               if(leadJet){
-                          nTracks = leadJet->GetRefTracks()->GetEntriesFast();
-                                  fh1JetPt->Fill(leadJet->Pt());
-                   fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
-                   fh1JetN->Fill(1);                              
-                          }
-           }
-           fh1TrackN->Fill((Float_t)nTracks);
-
-          for(Int_t it=0; it<nTracks; ++it){
-              AliAODTrack *tmpTr = 0x0;
-               if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
-               if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
-               if(!tmpTr) continue; 
-
-              tmpTr->SetStatus(AliESDtrack::kEmbedded);
-
-              new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
-              dummy = (*tracks)[nAODtracks-1];
 
-               fh1TrackPt->Fill(tmpTr->Pt());
-               fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
-          }
-
-          if(mcpartIN){
-              Int_t nMCpart = mcpartIN->GetEntriesFast();
-
-               Int_t nPhysicalPrimary=0;
-              Int_t nAODmcpart=0;
-              for(Int_t ip=0; ip<nMCpart; ++ip){
-                  AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
-
-                   if(fEmbedMode==kAODJetTracks){
-                      // jet track? else continue
-                      // not implemented yet
-                      continue;
-                   } 
-
-                  new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
-                  dummy = (*mcpartOUT)[nAODmcpart-1];
-
-                   if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
-                      fh1MCTrackPt->Fill(tmpPart->Pt());
-                      fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
-                      nPhysicalPrimary++;
-                   }
-              }
-               fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
-               
-          }
-       } // end: embed all tracks or jet tracks
-
-
-       if(fEmbedMode==kAODJet4Mom){
-
-          // loop over jets
-          Int_t nJets = aodJets->GetEntries();
-           fh1TrackN->Fill((Float_t)nJets);
-          for(Int_t ij=0; ij<nJets; ++ij){
-               AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+                  fh1JetPt->Fill(jet->Pt());
+                  fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
+                  nSelectedJets++;
+
+               }
+            }                             
+            fh1JetN->Fill(nSelectedJets);
+         }
+
+         if(fEmbedMode==kAODJetTracks){
+            // look for leading jet within selection
+            // get reference tracks for this jet
+            Float_t leadJetPt = 0.;
+            for(Int_t ij=0; ij<nJets; ++ij){
+               AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
                if(!jet) continue;
-              AliAODTrack *tmpTr = (AliAODTrack*)jet;
-              tmpTr->SetFlags(AliESDtrack::kEmbedded);
+               if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+                     && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+                     && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+                     && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
+                  if(jet->Pt()>leadJetPt){
+                     leadJetPt = jet->Pt();
+                     leadJet = jet;
+                  } 
+               }
+            }
+            if(leadJet){
+               nTracks = leadJet->GetRefTracks()->GetEntriesFast();
+               fh1JetPt->Fill(leadJet->Pt());
+               fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
+               fh1JetN->Fill(1);                                  
+            }
+         }
 
-              new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
-              dummy = (*tracks)[nAODtracks-1]; 
+         fh1TrackN->Fill((Float_t)nTracks);
 
+         for(Int_t it=0; it<nTracks; ++it){
+            AliAODTrack *tmpTr = 0x0;
+            if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
+            if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
+            if(!tmpTr) continue;
+
+            tmpTr->SetStatus(AliESDtrack::kEmbedded);
+
+            new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
+            dummy = (*tracks)[nAODtracks-1];
+
+            if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
                fh1TrackPt->Fill(tmpTr->Pt());
                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
-          }
-
-       } // end: embed jets as 4-momenta
-
-
-    } //end: embed mode with AOD
-
-
-    // === embed mode with toy ===
-    if(fEmbedMode==kToyTracks){
-        Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
-
-        fh1TrackN->Fill((Float_t)nT);
-
-        for(Int_t i=0; i<nT; ++i){
-
-          Double_t pt = -1.;
-           if(fToyMinTrackPt!=fToyMaxTrackPt){
-              if(fToyDistributionTrackPt==0){
-                 pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
-              } else {
-                 while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
-                    pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
-                    pt += fToyMinTrackPt;
-                 }
-              }
-           } else {
-              pt = fToyMinTrackPt;
-           }
-           Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
-          Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
-          phi = TVector2::Phi_0_2pi(phi);
-
-          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
-
-          Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
-          Float_t mom[3];
-          mom[0] = pt;
-          mom[1] = phi;
-          mom[2] = theta;
-
-          Float_t xyz[3];
-          xyz[0] = -999.;
-          xyz[1] = -999.;
-          xyz[2] = -999.;
-       
-          AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
-                                             -999,   // label
-                                             mom,    // momentum[3]
-                                             kFALSE, // cartesian
-                                             xyz,    // position
-                                             kFALSE, // DCA
-                                             NULL,   // covMatrix,
-                                             -99,    // charge
-                                             0,      // itsClusMap
-                                             NULL,   // pid 
-                                             NULL,   // prodVertex
-                                             kFALSE, // used for vertex fit
-                                             kFALSE, // used for prim vtx fit
-                                             AliAODTrack::kUndef, // type
-                                             fToyFilterMap,  // select info
-                                             -999.    // chi2 per NDF
-                                           );
-          tmpTr->SetFlags(AliESDtrack::kEmbedded);
-
-           new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
-           dummy = (*tracks)[nAODtracks-1];
-
-           fh1TrackPt->Fill(pt);
-           fh2TrackEtaPhi->Fill(eta,phi);
-
-          delete tmpTr;
-       }
-    } //end: embed mode with toy
-
-
-    PostData(1, fHistList);
+            }
+         }
+
+         if(mcpartIN){
+            Int_t nMCpart = mcpartIN->GetEntriesFast();
+
+            Int_t nPhysicalPrimary=0;
+            Int_t nAODmcpart=0;
+            for(Int_t ip=0; ip<nMCpart; ++ip){
+               AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
+
+               if(fEmbedMode==kAODJetTracks){
+                  // jet track? else continue
+                  // not implemented yet
+                  continue;
+               } 
+
+               new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
+               dummy = (*mcpartOUT)[nAODmcpart-1];
+
+               if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
+                  fh1MCTrackPt->Fill(tmpPart->Pt());
+                  fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
+                  nPhysicalPrimary++;
+               }
+            }
+            fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
+
+         }
+      } // end: embed all tracks or jet tracks
+
+
+      if(fEmbedMode==kAODJet4Mom){
+
+         // loop over jets
+         Int_t nJets = fAODJets->GetEntries();
+         fh1TrackN->Fill((Float_t)nJets);
+         for(Int_t ij=0; ij<nJets; ++ij){
+            AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
+            if(!jet) continue;
+            AliAODTrack *tmpTr = (AliAODTrack*)jet;
+            tmpTr->SetFlags(AliESDtrack::kEmbedded);
+
+            new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
+            dummy = (*tracks)[nAODtracks-1]; 
+
+            fh1TrackPt->Fill(tmpTr->Pt());
+            fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
+         }
+
+      } // end: embed jets as 4-momenta
+
+
+   } //end: embed mode with AOD
+
+
+   // === embed mode with toy ===
+   if(fEmbedMode==kToyTracks){
+      Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
+
+      fh1TrackN->Fill((Float_t)nT);
+
+      for(Int_t i=0; i<nT; ++i){
+
+         Double_t pt = -1.;
+         if(fToyMinTrackPt!=fToyMaxTrackPt){
+            if(fToyDistributionTrackPt==0){
+               pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
+            } else {
+               while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
+                  pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
+                  pt += fToyMinTrackPt;
+               }
+            }
+         } else {
+            pt = fToyMinTrackPt;
+         }
+         Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
+         Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
+         phi = TVector2::Phi_0_2pi(phi);
+
+         if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
+
+         Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
+         Float_t mom[3];
+         mom[0] = pt;
+         mom[1] = phi;
+         mom[2] = theta;
+
+         Float_t xyz[3];
+         xyz[0] = -999.;
+         xyz[1] = -999.;
+         xyz[2] = -999.;
+
+         AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
+         -999,   // label
+         mom,    // momentum[3]
+         kFALSE, // cartesian
+         xyz,    // position
+         kFALSE, // DCA
+         NULL,   // covMatrix,
+         -99,    // charge
+         0,      // itsClusMap
+         NULL,   // pid 
+         NULL,   // prodVertex
+         kFALSE, // used for vertex fit
+         kFALSE, // used for prim vtx fit
+         AliAODTrack::kUndef, // type
+         fToyFilterMap,  // select info
+         -999.    // chi2 per NDF
+         );
+         tmpTr->SetFlags(AliESDtrack::kEmbedded);
+
+         new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
+         dummy = (*tracks)[nAODtracks-1];
+
+         fh1TrackPt->Fill(pt);
+         fh2TrackEtaPhi->Fill(eta,phi);
+
+         delete tmpTr;
+      }
+   } //end: embed mode with toy
+
+
+   PostData(1, fHistList);
 }
 
 
 //__________________________________________________________________________
 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
 {
-    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
+   if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
 
-    if(fAODfile && fAODfile->IsOpen())
-    fAODfile->Close();  
+   if(fAODfile && fAODfile->IsOpen())
+   fAODfile->Close();  
 
 }
 
@@ -708,11 +978,11 @@ Int_t AliAnalysisTaskFastEmbedding::GetJobID()
    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
 
    if(env && strlen(env)){
-       id= atoi(env);
-       AliInfo(Form("Job index %d", id));
+      id= atoi(env);
+      AliInfo(Form("Job index %d", id));
    }
    else{
-       AliInfo("Job index not found. Okay if running locally.");
+      AliInfo("Job index not found. Okay if running locally.");
    }
 
    return id;
@@ -722,60 +992,205 @@ Int_t AliAnalysisTaskFastEmbedding::GetJobID()
 
 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
 
-     Int_t nFiles = fAODPathArray->GetEntries();
-     Int_t n = rndm->Integer(nFiles);
-     AliInfo(Form("Select AOD file %d", n));
-     TObjString *objStr = (TObjString*) fAODPathArray->At(n);
-     if(!objStr){
-          AliError("Could not get path of aod file.");
-          return -1;
-     }
-     fAODPath = objStr->GetString();
-
-     return n;
+   Int_t nFiles = fAODPathArray->GetEntries();
+
+   // choose randomly file and event
+   Int_t n = -1;
+   Float_t tmpProp = -1.;
+   Int_t pendingEvents = fInputEntries-Entry();
+   //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
+   while(rndm->Rndm()>tmpProp){
+      n = rndm->Integer(nFiles);
+      fAODEntries = fAODEntriesArray->At(n);
+      Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
+      tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
+   }
+   fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
+
+   AliInfo(Form("Select AOD file %d", n));
+   TObjString *objStr = (TObjString*) fAODPathArray->At(n);
+   if(!objStr){
+      AliError("Could not get path of aod file.");
+      return -1;
+   }
+   fAODPath = objStr->GetString();
+
+   return n;
 }
 
 //__________________________________________________________________________
 
 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
 
-    if(trial>3){
-        AliError("Could not open AOD files, give up ...");
-        return -1;
-    }
-       
-       Int_t rc = 0;
-       if(fAODPathArray) rc = SelectAODfile();
-       if(rc<0) return -1;
-
-    TDirectory *owd = gDirectory;
-    if (fAODfile)
-      fAODfile->Close();
-    fAODfile = TFile::Open(fAODPath.Data());
-    owd->cd();
-    if(!fAODfile){
-       
-       rc = -1;
-       if(fAODPathArray){
-           AliError(Form("Could not open AOD file %s, try another from the list ...", fAODPath.Data()));
-           rc = OpenAODfile(trial+1);
-       } else {
-              AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
-          }
-        
-       return rc;
-    }
-
-    fAODtree = (TTree*)fAODfile->Get("aodTree");
-
-    if(!fAODtree){
-       AliError("AOD tree not found.");
-       return -1;
-    }
-
-    delete fAODevent;
-    fAODevent = new AliAODEvent();
-    fAODevent->ReadFromTree(fAODtree);
-    return rc;  // file position in AOD path array, if array available
+   if(fAODPathArray) fFileId = SelectAODfile();
+   if(fFileId<0) return -1;
+
+   TDirectory *owd = gDirectory;
+   if (fAODfile) fAODfile->Close();
+   fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
+   owd->cd();
+   if(!fAODfile){
+
+      fFileId = -1;
+      if(fAODPathArray){
+         if(trial<=3){ 
+            AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
+            fFileId = OpenAODfile(trial+1);
+         } else {
+            AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
+            return -1;
+         }
+      } else {
+         AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
+      }
+
+      return fFileId;
+   }
+
+   fAODtree = (TTree*)fAODfile->Get("aodTree");
+
+   if(!fAODtree){
+      AliError("AOD tree not found.");
+      return -1;
+   }
+
+   /*
+      fAODtree->SetBranchStatus("*",0);
+      fAODtree->SetBranchStatus("header",1);
+      fAODtree->SetBranchStatus("tracks*",1);
+      fAODtree->SetBranchStatus("jets*",1);
+      fAODtree->SetBranchStatus("clusters*",1);
+      fAODtree->SetBranchStatus("mcparticles*",1);
+      fAODtree->SetBranchStatus("mcHeader*",1);
+   */
+
+   delete fAODevent;
+   fAODevent = new AliAODEvent();
+   fAODevent->ReadFromTree(fAODtree);
+
+
+   // fetch header, jets, etc. from new file
+   fNevents = fAODtree->GetEntries();
+   mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
+
+   if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
+   else                    fAODJets = fAODevent->GetJets();
+   if(!fAODJets){
+      AliError("Could not find jets in AOD. Check jet branch when indicated.");
+      return -1;
+   }
+
+   TFile *curfile = fAODtree->GetCurrentFile();
+   if (!curfile) {
+      AliError("No current file.");
+      return -1;
+   }
+
+   Float_t trials = 1.;
+   Float_t xsec = 0.;
+   PythiaInfoFromFile(curfile->GetName(),xsec,trials);
+   fXsection = xsec;
+
+   // construct a poor man average trials 
+   Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
+   if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
+
+   if(fFileId>=0){
+      AliInfo(Form("Read successfully AOD event from file %d",fFileId));
+   }
+
+   fCountEvents=0; // new file, reset counter
+
+   return fFileId;  // file position in AOD path array, if array available
+}
+
+
+//____________________________________________________________________________
+Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
+
+   static Float_t ptHard = -1.;
+   if(bSet) ptHard = newValue;
+
+   return ptHard;
 }
 
+//____________________________________________________________________________
+Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
+
+   const Int_t nBins = 10;
+   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
+
+   Int_t bin = -1;
+   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
+      bin++;
+   }
+
+   return bin;
+}
+
+//____________________________________________________________________________
+Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
+   //
+   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
+   // This should provide the path to the AOD/ESD file
+   // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks) 
+
+   TString file(currFile);  
+   fXsec = 0;
+   fTrials = 1;
+
+   if(file.Contains("root_archive.zip#")){
+      Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
+      Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
+      file.Replace(pos+1,20,"");
+   }
+   else {
+      // not an archive take the basename....
+      file.ReplaceAll(gSystem->BaseName(file.Data()),"");
+   }
+
+   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
+   if(!fxsec){
+      // next trial fetch the histgram file
+      fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
+      if(!fxsec){
+         // not a severe condition but inciate that we have no information
+         AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
+         return kFALSE;
+      }
+      else{
+         // find the tlist we want to be independtent of the name so use the Tkey
+         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
+         if(!key){
+            fxsec->Close();
+            return kFALSE;
+         }
+         TList *list = dynamic_cast<TList*>(key->ReadObj());
+         if(!list){
+            fxsec->Close();
+            return kFALSE;
+         }
+         fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+         fxsec->Close();
+      }
+   } // no tree pyxsec.root
+   else {
+      TTree *xtree = (TTree*)fxsec->Get("Xsection");
+      if(!xtree){
+         fxsec->Close();
+         return kFALSE;
+      }
+      UInt_t   ntrials  = 0;
+      Double_t  xsection  = 0;
+      xtree->SetBranchAddress("xsection",&xsection);
+      xtree->SetBranchAddress("ntrials",&ntrials);
+      xtree->GetEntry(0);
+      fTrials = ntrials;
+      fXsec = xsection;
+      fxsec->Close();
+   }
+   return kTRUE;
+}
+
+
index db40f3b..d17a43c 100644 (file)
@@ -2,12 +2,13 @@
 #define ALIANALYSISTASKFASTEMBEDDING_H
 
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
+* See cxx source for full Copyright notice                               */
 
 /* $Id$ */
 
 #include "AliAnalysisTaskSE.h"
 
+class AliESDEvent;
 class AliAODEvent;
 class TTree;
 class TFile;
@@ -17,117 +18,179 @@ class TObjString;
 class TRandom3;
 class TH1F;
 class TH2F;
+class TProfile;
+class AliAODMCHeader;
 
 class AliAnalysisTaskFastEmbedding : public AliAnalysisTaskSE {
 
-    public:
-       
-       AliAnalysisTaskFastEmbedding();
-       AliAnalysisTaskFastEmbedding(const char *name);
-       AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy);
-       AliAnalysisTaskFastEmbedding& operator=(const AliAnalysisTaskFastEmbedding &o);
-       virtual ~AliAnalysisTaskFastEmbedding();
-
-       virtual void UserCreateOutputObjects();
-       virtual void LocalInit() { Init(); }
-       virtual void Init();
-        virtual void UserExec(Option_t*);
-       virtual void Terminate(Option_t */*option*/);
-
-       void SetAODPath(TString path) {fAODPath = path;}
-       void SetArrayOfAODPaths(TObjArray* arr) {fAODPathArray = arr;}
-        void SetTrackBranch(TString name) {fTrackBranch = name;}
-        void SetMCparticlesBranch(TString name) {fMCparticlesBranch = name;}
-        void SetJetBranch(TString name) {fJetBranch = name;}
-
-       void SetEmbedMode(Int_t m) {fEmbedMode = m;}
-       Int_t GetEmbedMode() {return fEmbedMode;} 
-       void SetEvtSelecMode(Int_t s) {fEvtSelecMode = s;}
-       Int_t GetEvtSelecMode() {return fEvtSelecMode;}
-
-       void SetEvtSelJetPtRange(Float_t minPt, Float_t maxPt) {fEvtSelMinJetPt = minPt; fEvtSelMaxJetPt = maxPt;}
-       void SetEvtSelJetEtaRange(Float_t minEta, Float_t maxEta) {fEvtSelMinJetEta = minEta; fEvtSelMaxJetEta = maxEta;}
-       void SetEvtSelJetPhiRange(Float_t minPhi, Float_t maxPhi) {fEvtSelMinJetPhi = minPhi; fEvtSelMaxJetPhi = maxPhi;}
-       
-        void SetToyNumberOfTrackRange(Int_t minN = 1, Int_t maxN = 1){ fToyMinNbOfTracks = minN, fToyMaxNbOfTracks = maxN; }
-       void SetToyTrackRanges(Double_t minPt = 50., Double_t maxPt = 50., Double_t ptDistr=0,
-               Double_t minEta = -.5, Double_t maxEta = .5,
-               Double_t minPhi = 0., Double_t maxPhi = 2*TMath::Pi())
-                {
-               fToyMinTrackPt = minPt; fToyMaxTrackPt = maxPt; fToyDistributionTrackPt = ptDistr;
-               fToyMinTrackEta = minEta; fToyMaxTrackEta = maxEta;
-               fToyMinTrackPhi = minPhi; fToyMaxTrackPhi = maxPhi;}
-        void SetToyFilterMap(UInt_t f) {fToyFilterMap = f;}
-
-
-        // embedding modes
-       enum {kAODFull=0, kAODJetTracks, kAODJet4Mom, kToyTracks};
-       // event selection from AOD
-       enum {kEventsAll=0, kEventsJetPt};
-
-
-    private:
-
-        AliAODEvent* fAODout;    //! AOD out
-       AliAODEvent* fAODevent;  //! AOD in
-       TTree* fAODtree;         //! AODin tree
-        TFile* fAODfile;         //! AODin file
-       TRandom3* rndm;           //! random nummer generator
-
-       TObjArray* fAODPathArray;  // array of paths of AOD in file
-       TString fAODPath;  // path of AOD in file
-
-        TString fTrackBranch; // name of branch for extra tracks in AOD out
-        TString fMCparticlesBranch; // name of branch for extra mcparticles in AOD out
-        TString fJetBranch; // name of branch for extra jets AOD in
-
-        Int_t fEntry; // entry of extra AOD
-
-       Int_t fEmbedMode;
-       Int_t fEvtSelecMode;
-
-       // event selection from AOD
-       Float_t fEvtSelMinJetPt;       // minimum pt of the leading jet
-       Float_t fEvtSelMaxJetPt;       // maximum pt of the leading jet
-        Float_t fEvtSelMinJetEta;      // minimum eta of the leading jet
-        Float_t fEvtSelMaxJetEta;      // maximum eta of the leading jet
-        Float_t fEvtSelMinJetPhi;      // minimum phi of the leading jet
-        Float_t fEvtSelMaxJetPhi;      // maximum phi of the leading jet
-        
-         
-        // settings for toy "track generation"
-        Int_t    fToyMinNbOfTracks;             // minimum nb. of tracks per event
-        Int_t    fToyMaxNbOfTracks;             // maximum nb. of tracks per event
-       Double_t fToyMinTrackPt;                // minimum track pT
-       Double_t fToyMaxTrackPt;                // maximum track pT
-        Double_t fToyDistributionTrackPt;       // distribution of track pt
-       Double_t fToyMinTrackEta;               // minimum eta of tracks
-       Double_t fToyMaxTrackEta;               // maximum eta of tracks
-       Double_t fToyMinTrackPhi;               // minimum phi of tracks
-       Double_t fToyMaxTrackPhi;               // maximum phi of tracks
-       UInt_t fToyFilterMap;                   // filter map of tracks
-
-
-        // qa histos
-        TList *fHistList;          //  list of histograms
-        TH1F  *fh1TrackPt;         //! track pt
-        TH2F  *fh2TrackEtaPhi;     //! track eta-phi
-        TH1F  *fh1TrackN;          //! nb. of tracks
-        TH1F  *fh1JetPt;           //! jet pt
-        TH2F  *fh2JetEtaPhi;       //! jet eta-phi
-        TH1F  *fh1JetN;            //! nb. of jets
-        TH1F  *fh1MCTrackPt;       //! MC track pt
-        TH2F  *fh2MCTrackEtaPhi;   //! MC track eta-phi
-        TH1F  *fh1MCTrackN;        //! nb. of MC tracks
-        TH1I  *fh1AODfile;         //! used AOD files from AODPathArray
-               
-
-       Int_t GetJobID();    // get job id (sub-job id on the GRID)
-        Int_t SelectAODfile();
-        Int_t OpenAODfile(Int_t trial = 0);
-
-
-       ClassDef(AliAnalysisTaskFastEmbedding, 4);
+public:
+   
+   AliAnalysisTaskFastEmbedding();
+   AliAnalysisTaskFastEmbedding(const char *name);
+   AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy);
+   AliAnalysisTaskFastEmbedding& operator=(const AliAnalysisTaskFastEmbedding &o);
+   virtual ~AliAnalysisTaskFastEmbedding();
+
+   virtual void UserCreateOutputObjects();
+   virtual void LocalInit() { Init(); }
+   virtual void Init();
+   virtual Bool_t UserNotify();
+   virtual void UserExec(Option_t*);
+   virtual void Terminate(Option_t */*option*/);
+
+   void SetAODPath(TString path) {fAODPath = path;}
+   void SetArrayOfAODPaths(TObjArray* arr) {fAODPathArray = arr;}
+   void SetArrayOfAODEntries(TArrayI* arr) {fAODEntriesArray = arr;}
+   void SetAODEntriesSum(Int_t i){ fAODEntriesSum = i;}
+   void SetAODEntriesMax(Int_t i){ fAODEntriesMax = i;}
+   
+   virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+   virtual void     SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
+   virtual void     SetVtxZMin(Float_t z) { fVtxZMin = z; }
+   virtual void     SetVtxZMax(Float_t z) { fVtxZMax = z; }
+   virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
+   virtual void     SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; }
+   virtual void     SetCentMin(Float_t cent) { fCentMin = cent; }
+   virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
+   virtual void     SetNInputTracksMin(Int_t nTr) { fNInputTracksMin = nTr; }
+   virtual void     SetNInputTracksMax(Int_t nTr) { fNInputTracksMax = nTr; }
+   
+   void SetTrackBranch(TString name) {fTrackBranch = name;}
+   void SetMCparticlesBranch(TString name) {fMCparticlesBranch = name;}
+   void SetJetBranch(TString name) {fJetBranch = name;}
+
+   void SetEmbedMode(Int_t m) {fEmbedMode = m;}
+   Int_t GetEmbedMode() {return fEmbedMode;} 
+   void SetEvtSelecMode(Int_t s) {fEvtSelecMode = s;}
+   Int_t GetEvtSelecMode() {return fEvtSelecMode;}
+
+   void SetEvtSelJetPtRange(Float_t minPt, Float_t maxPt) {fEvtSelMinJetPt = minPt; fEvtSelMaxJetPt = maxPt;}
+   void SetEvtSelJetEtaRange(Float_t minEta, Float_t maxEta) {fEvtSelMinJetEta = minEta; fEvtSelMaxJetEta = maxEta;}
+   void SetEvtSelJetPhiRange(Float_t minPhi, Float_t maxPhi) {fEvtSelMinJetPhi = minPhi; fEvtSelMaxJetPhi = maxPhi;}
+   
+   void SetToyNumberOfTrackRange(Int_t minN = 1, Int_t maxN = 1){ fToyMinNbOfTracks = minN, fToyMaxNbOfTracks = maxN; }
+   void SetToyTrackRanges(Double_t minPt = 50., Double_t maxPt = 50., Double_t ptDistr=0,
+   Double_t minEta = -.5, Double_t maxEta = .5,
+   Double_t minPhi = 0., Double_t maxPhi = 2*TMath::Pi())
+   {
+      fToyMinTrackPt = minPt; fToyMaxTrackPt = maxPt; fToyDistributionTrackPt = ptDistr;
+      fToyMinTrackEta = minEta; fToyMaxTrackEta = maxEta;
+      fToyMinTrackPhi = minPhi; fToyMaxTrackPhi = maxPhi;}
+   void SetToyFilterMap(UInt_t f) {fToyFilterMap = f;}
+   void SetTrackFilterMap(UInt_t f) {fTrackFilterMap = f;}
+
+   static Float_t GetPtHard(Bool_t bSet=kFALSE, Float_t newValue = 0.);
+   
+   virtual Int_t      GetPtHardBin(Double_t ptHard);
+   virtual Bool_t     PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials);
+
+   // embedding modes
+   enum {kAODFull=0, kAODJetTracks, kAODJet4Mom, kToyTracks};
+   // event selection from AOD
+   enum {kEventsAll=0, kEventsJetPt};
+
+
+private:
+
+   AliESDEvent    *fESD;        //! ESD object
+   AliAODEvent    *fAODout;     //! AOD out
+   AliAODEvent    *fAODevent;   //! AOD in
+   TTree          *fAODtree;    //! AODin tree
+   TFile          *fAODfile;    //! AODin file
+   AliAODMCHeader *mcHeader;    //! mc header
+   TRandom3       *rndm;        //! random nummer generator
+   Int_t          fInputEntries; // total nb. of events (for this subjob)
+
+   TObjArray *fAODPathArray;    // array of paths of AOD in file
+   TArrayI   *fAODEntriesArray; // array of entries of AODs 
+   TString    fAODPath;         // path of AOD in file
+   Int_t      fAODEntries;      // entries of AOD
+   Int_t      fAODEntriesSum;   // sum of all entries of AODs
+   Int_t      fAODEntriesMax;   // maximum entries of AODs
+
+   AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept
+   Int_t   fMinContribVtx; // minimum number of track contributors for primary vertex
+   Float_t fVtxZMin;         // lower bound on vertex z
+   Float_t fVtxZMax;         // upper bound on vertex z
+   Int_t   fEvtClassMin;       // lower bound on event class
+   Int_t   fEvtClassMax;       // upper bound on event class
+   Float_t fCentMin;         // lower bound on centrality
+   Float_t fCentMax;         // upper bound on centrality
+   Int_t   fNInputTracksMin;  // lower bound of nb. of input tracks
+   Int_t   fNInputTracksMax;  // upper bound of nb. of input tracks
+   
+   TString fTrackBranch;       // name of branch for extra tracks in AOD out
+   TString fMCparticlesBranch; // name of branch for extra mcparticles in AOD out
+   TString fJetBranch;         // name of branch for extra jets AOD in
+
+   Int_t fFileId;   // nb. of file from the list
+   Int_t fAODEntry; // entry of extra AOD
+   Int_t fCountEvents; // count processed events in this file
+
+   Int_t fEmbedMode;
+   Int_t fEvtSelecMode;
+
+   // event selection from AOD
+   Float_t fEvtSelMinJetPt;       // minimum pt of the leading jet
+   Float_t fEvtSelMaxJetPt;       // maximum pt of the leading jet
+   Float_t fEvtSelMinJetEta;      // minimum eta of the leading jet
+   Float_t fEvtSelMaxJetEta;      // maximum eta of the leading jet
+   Float_t fEvtSelMinJetPhi;      // minimum phi of the leading jet
+   Float_t fEvtSelMaxJetPhi;      // maximum phi of the leading jet
+   
+   
+   // settings for toy "track generation"
+   Int_t    fToyMinNbOfTracks;             // minimum nb. of tracks per event
+   Int_t    fToyMaxNbOfTracks;             // maximum nb. of tracks per event
+   Float_t  fToyMinTrackPt;                // minimum track pT
+   Float_t  fToyMaxTrackPt;                // maximum track pT
+   Float_t  fToyDistributionTrackPt;       // distribution of track pt
+   Float_t  fToyMinTrackEta;               // minimum eta of tracks
+   Float_t  fToyMaxTrackEta;               // maximum eta of tracks
+   Float_t  fToyMinTrackPhi;               // minimum phi of tracks
+   Float_t  fToyMaxTrackPhi;               // maximum phi of tracks
+   UInt_t   fToyFilterMap;                 // filter map of tracks
+   UInt_t   fTrackFilterMap;               // filter map of tracks for QA plots
+   
+   Int_t         fNPtHard;      // nb. of pT hard bins
+   Double_t      fPtHard;       // pT hard
+   Int_t         fPtHardBin;    // pT hard bin
+   TClonesArray* fAODJets;      //! array of jets from aod
+   Int_t         fNevents;      // number of events in aod
+   Float_t       fXsection;     // average xsection of the event
+   Float_t       fAvgTrials;    // average number of trials per event
+
+
+   // histos
+   TList *fHistList;          //  list of histograms
+   TH1I  *fHistEvtSelection;  //! stastic of event selection
+   TProfile *fh1Xsec;         //! cross-section
+   TH1F  *fh1Trials;          //! nb. of trials (simulation)
+   TH1F  *fh1TrialsEvtSel;    //! nb. of trials (event selection, e.g. jet pT)
+   TH2F  *fh2PtHard;          //! pT hard bin
+   TH2F  *fh2PtHardEvtSel;    //! pT hard bin (event selection)
+   TH2F  *fh2PtHardTrials;     //! pT hard bin, weighted by nb. of trials
+   
+   // qa histos
+   TH1F  *fh1TrackPt;         //! track pt
+   TH2F  *fh2TrackEtaPhi;     //! track eta-phi
+   TH1F  *fh1TrackN;          //! nb. of tracks
+   TH1F  *fh1JetPt;           //! jet pt
+   TH2F  *fh2JetEtaPhi;       //! jet eta-phi
+   TH1F  *fh1JetN;            //! nb. of jets
+   TH1F  *fh1MCTrackPt;       //! MC track pt
+   TH2F  *fh2MCTrackEtaPhi;   //! MC track eta-phi
+   TH1F  *fh1MCTrackN;        //! nb. of MC tracks
+   TH1I  *fh1AODfile;         //! used AOD files from AODPathArray
+   TH2I  *fh2AODevent;        //! selected events in AODs
+   
+
+   Int_t GetJobID();    // get job id (sub-job id on the GRID)
+   Int_t SelectAODfile();
+   Int_t OpenAODfile(Int_t trial = 0);
+
+
+   ClassDef(AliAnalysisTaskFastEmbedding, 5);
 };
 
 #endif
index 16d3d3d..d9f5dbf 100644 (file)
@@ -19,6 +19,7 @@
 #include "AliAnalysisHelperJetTasks.h"
 #include "AliInputEventHandler.h"
 #include "AliAODJetEventBackground.h"
+#include "AliAnalysisTaskFastEmbedding.h"
 
 #include "AliAODEvent.h"
 #include "AliAODJet.h"
 ClassImp(AliAnalysisTaskJetResponseV2)
 
 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
-  AliAnalysisTaskSE(),
-  fESD(0x0),
-  fAOD(0x0),
-  fBackgroundBranch(""),
-  fIsPbPb(kTRUE),
-  fOfflineTrgMask(AliVEvent::kAny),
-  fMinContribVtx(1),
-  fVtxZMin(-8.),
-  fVtxZMax(8.),
-  fEvtClassMin(0),
-  fEvtClassMax(4),
-  fCentMin(0.),
-  fCentMax(100.),
-  fNInputTracksMin(0),
-  fNInputTracksMax(-1),
-  fJetEtaMin(-.5),
-  fJetEtaMax(.5),
-  fJetPtMin(20.),
-  fJetPtFractionMin(0.5),
-  fNMatchJets(4),
-  fMatchMaxDist(0.8),
-  fkNbranches(2),
-  fkEvtClasses(12),
-  fOutputList(0x0),
-  fbEvent(kTRUE),
-  fbJetsMismatch1(kTRUE),
-  fbJetsMismatch2(kTRUE),
-  fbJetsRp(kTRUE),
-  fbJetsDeltaPt(kTRUE),
-  fbJetsEta(kTRUE),
-  fbJetsPhi(kTRUE),
-  fbJetsArea(kTRUE),
-  fbJetsBeforeCut1(kTRUE),
-  fbJetsBeforeCut2(kTRUE),
-  fHistEvtSelection(0x0),
-  fHistJetSelection(0x0),
-  fhnEvent(0x0),
-  fhnJetsMismatch1(0x0),
-  fhnJetsMismatch2(0x0),
-  fhnJetsRp(0x0),
-  fhnJetsDeltaPt(0x0),
-  fhnJetsEta(0x0),
-  fhnJetsPhi(0x0),
-  fhnJetsArea(0x0),
-  fhnJetsBeforeCut1(0x0),
-  fhnJetsBeforeCut2(0x0)
+AliAnalysisTaskSE(),
+fESD(0x0),
+fAOD(0x0),
+fBackgroundBranch(""),
+fIsPbPb(kTRUE),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-8.),
+fVtxZMax(8.),
+fEvtClassMin(0),
+fEvtClassMax(4),
+fCentMin(0.),
+fCentMax(100.),
+fNInputTracksMin(0),
+fNInputTracksMax(-1),
+fJetEtaMin(-.5),
+fJetEtaMax(.5),
+fJetPtMin(20.),
+fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
+fJetPtFractionMin(0.5),
+fNMatchJets(4),
+fMatchMaxDist(0.8),
+fKeepJets(kFALSE),
+fkNbranches(2),
+fkEvtClasses(12),
+fOutputList(0x0),
+fbEvent(kTRUE),
+fbJetsMismatch1(kTRUE),
+fbJetsMismatch2(kTRUE),
+fbJetsRp(kTRUE),
+fbJetsDeltaPt(kTRUE),
+fbJetsEta(kTRUE),
+fbJetsPhi(kTRUE),
+fbJetsArea(kTRUE),
+fbJetsBeforeCut1(kTRUE),
+fbJetsBeforeCut2(kTRUE),
+fHistEvtSelection(0x0),
+fHistJetSelection(0x0),
+fh2JetSelection(0x0),
+fhnEvent(0x0),
+fhnJetsMismatch1(0x0),
+fhnJetsMismatch2(0x0),
+fhnJetsRp(0x0),
+fhnJetsDeltaPt(0x0),
+fhnJetsEta(0x0),
+fhnJetsPhi(0x0),
+fhnJetsArea(0x0),
+fhnJetsBeforeCut1(0x0),
+fhnJetsBeforeCut2(0x0)
 {
-  // default Constructor
+   // default Constructor
 
-  fJetBranchName[0] = "";
-  fJetBranchName[1] = "";
+   fJetBranchName[0] = "";
+   fJetBranchName[1] = "";
 
-  fListJets[0] = new TList;
-  fListJets[1] = new TList;
+   fListJets[0] = new TList;
+   fListJets[1] = new TList;
 }
 
 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
-  AliAnalysisTaskSE(name),
-  fESD(0x0),
-  fAOD(0x0),
-  fBackgroundBranch(""),
-  fIsPbPb(kTRUE),
-  fOfflineTrgMask(AliVEvent::kAny),
-  fMinContribVtx(1),
-  fVtxZMin(-8.),
-  fVtxZMax(8.),
-  fEvtClassMin(0),
-  fEvtClassMax(4),
-  fCentMin(0.),
-  fCentMax(100.),
-  fNInputTracksMin(0),
-  fNInputTracksMax(-1),
-  fJetEtaMin(-.5),
-  fJetEtaMax(.5),
-  fJetPtMin(20.),
-  fJetPtFractionMin(0.5),
-  fNMatchJets(4),
-  fMatchMaxDist(0.8),
-  fkNbranches(2),
-  fkEvtClasses(12),
-  fOutputList(0x0),
-  fbEvent(kTRUE),
-  fbJetsMismatch1(kTRUE),
-  fbJetsMismatch2(kTRUE),
-  fbJetsRp(kTRUE),
-  fbJetsDeltaPt(kTRUE),
-  fbJetsEta(kTRUE),
-  fbJetsPhi(kTRUE),
-  fbJetsArea(kTRUE),
-  fbJetsBeforeCut1(kTRUE),
-  fbJetsBeforeCut2(kTRUE),
-  fHistEvtSelection(0x0),
-  fHistJetSelection(0x0),
-  fhnEvent(0x0),
-  fhnJetsMismatch1(0x0),
-  fhnJetsMismatch2(0x0),
-  fhnJetsRp(0x0),
-  fhnJetsDeltaPt(0x0),
-  fhnJetsEta(0x0),
-  fhnJetsPhi(0x0),
-  fhnJetsArea(0x0),
-  fhnJetsBeforeCut1(0x0),
-  fhnJetsBeforeCut2(0x0)
+AliAnalysisTaskSE(name),
+fESD(0x0),
+fAOD(0x0),
+fBackgroundBranch(""),
+fIsPbPb(kTRUE),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-8.),
+fVtxZMax(8.),
+fEvtClassMin(0),
+fEvtClassMax(4),
+fCentMin(0.),
+fCentMax(100.),
+fNInputTracksMin(0),
+fNInputTracksMax(-1),
+fJetEtaMin(-.5),
+fJetEtaMax(.5),
+fJetPtMin(20.),
+fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
+fJetPtFractionMin(0.5),
+fNMatchJets(4),
+fMatchMaxDist(0.8),
+fKeepJets(kFALSE),
+fkNbranches(2),
+fkEvtClasses(12),
+fOutputList(0x0),
+fbEvent(kTRUE),
+fbJetsMismatch1(kTRUE),
+fbJetsMismatch2(kTRUE),
+fbJetsRp(kTRUE),
+fbJetsDeltaPt(kTRUE),
+fbJetsEta(kTRUE),
+fbJetsPhi(kTRUE),
+fbJetsArea(kTRUE),
+fbJetsBeforeCut1(kTRUE),
+fbJetsBeforeCut2(kTRUE),
+fHistEvtSelection(0x0),
+fHistJetSelection(0x0),
+fh2JetSelection(0x0),
+fhnEvent(0x0),
+fhnJetsMismatch1(0x0),
+fhnJetsMismatch2(0x0),
+fhnJetsRp(0x0),
+fhnJetsDeltaPt(0x0),
+fhnJetsEta(0x0),
+fhnJetsPhi(0x0),
+fhnJetsArea(0x0),
+fhnJetsBeforeCut1(0x0),
+fhnJetsBeforeCut2(0x0)
 {
-  // Constructor
+   // Constructor
 
-  fJetBranchName[0] = "";
-  fJetBranchName[1] = "";
+   fJetBranchName[0] = "";
+   fJetBranchName[1] = "";
 
-  fListJets[0] = new TList;
-  fListJets[1] = new TList;
+   fListJets[0] = new TList;
+   fListJets[1] = new TList;
 
-  DefineOutput(1, TList::Class());
+   DefineOutput(1, TList::Class());
 }
 
 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
 {
-  delete fListJets[0];
-  delete fListJets[1];
+   delete fListJets[0];
+   delete fListJets[1];
 }
 
 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2)
 {
-  fJetBranchName[0] = branch1;
-  fJetBranchName[1] = branch2;
+   fJetBranchName[0] = branch1;
+   fJetBranchName[1] = branch2;
 }
 
 void AliAnalysisTaskJetResponseV2::Init()
@@ -167,550 +174,599 @@ void AliAnalysisTaskJetResponseV2::Init()
 
 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
 {
-  // Create histograms
-  // Called once
-  OpenFile(1);
-  if(!fOutputList) fOutputList = new TList;
-  fOutputList->SetOwner(kTRUE);
-
-  Bool_t oldStatus = TH1::AddDirectoryStatus();
-  TH1::AddDirectory(kFALSE);
-
-
-  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
-  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
-  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
-  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
-  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
-  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
-  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
-  
-  fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 7, -0.5, 6.5);
-  fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
-  fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
-  fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
-  fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
-  fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
-  fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
-  fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
-
-  
-  UInt_t entries = 0; // bit coded, see GetDimParams() below
-  UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
-  
-  if(fbEvent){
-   entries = 1<<0 | 1<<1 | 1<<2;  // cent : nInpTrks : rp psi
-   opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
-   fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
-  }
+   // Create histograms
+   // Called once
+   OpenFile(1);
+   if(!fOutputList) fOutputList = new TList;
+   fOutputList->SetOwner(kTRUE);
+
+   Bool_t oldStatus = TH1::AddDirectoryStatus();
+   TH1::AddDirectory(kFALSE);
+
+
+   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
+   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
+
+   fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
+   fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
+   fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
+   fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
+   fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
+
+   fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
+   fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
+   fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
+   fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
+   fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
+
+
+   UInt_t entries = 0; // bit coded, see GetDimParams() below
+   UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
+
+   if(fbEvent){
+      entries = 1<<0 | 1<<1 | 1<<2 | 1<<26;  // cent : nInpTrks : rp psi : pT hard bin
+      opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
+      fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
+   }
+   
+   if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
+      // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
+      opt =  1<<6 | 1<<8 | 1<<10;
+      fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
+   }
+
+   if(fbJetsMismatch2){  // acceptance + fraction cut
+      // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
+      opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
+      fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
+
+   }
    
-  if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
-   // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12;
-   opt =  1<<6 | 1<<8 | 1<<10;
-   fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
-  }
-  
-  if(fbJetsMismatch2){  // acceptance + fraction cut
-  // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19;
-   opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
-   fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
-  
-  }
    
-     
-  if(fbJetsRp){
-   // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
-   /*
+   if(fbJetsRp){
+      // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
+      /*
    entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25;
    opt = 1<<4 | 1<<5;*/
-   // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23;
-   opt = 1<<4 | 1<<5;
-   fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
-  }
-  
-  // cent : nInpTrks : rp bins:  deltaPt : jetPt(2x) : deltaArea (hr delta pt)
-  if(fbJetsDeltaPt){
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
-   opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
-   fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
-  }
-  
-  // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (hr for eta)
-  if(fbJetsEta){
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
-   opt = 1<<15 | 1<<8 | 1<<9;
-   fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
-  }
-  
-  // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi
-  if(fbJetsPhi){
-    entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16;
-    opt = 1<<10 | 1<<11;
-    fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
-  }
-  
-  // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) (hr for area)
-  if(fbJetsArea){
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7;
-   opt = 1<<18 | 1<<12 | 1<<13;
-   fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
-  }
-  
-  
-  //before cut
-  
-  // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
-  if(fbJetsBeforeCut1){
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
-   opt = 0;
-   fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
-  }
-  
-  // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
-  if(fbJetsBeforeCut2){
-   entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
-   opt = 0;
-   fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
-  }
-  
-  fOutputList->Add(fHistEvtSelection);
-  fOutputList->Add(fHistJetSelection);
-  if(fbEvent)           fOutputList->Add(fhnEvent);
-  if(fbJetsMismatch1)   fOutputList->Add(fhnJetsMismatch1);
-  if(fbJetsMismatch2)   fOutputList->Add(fhnJetsMismatch2);
-  if(fbJetsRp)          fOutputList->Add(fhnJetsRp);
-  if(fbJetsDeltaPt)     fOutputList->Add(fhnJetsDeltaPt);
-  if(fbJetsEta)         fOutputList->Add(fhnJetsEta);
-  if(fbJetsPhi)         fOutputList->Add(fhnJetsPhi);
-  if(fbJetsArea)        fOutputList->Add(fhnJetsArea);
-  if(fbJetsBeforeCut1)  fOutputList->Add(fhnJetsBeforeCut1);
-  if(fbJetsBeforeCut2)  fOutputList->Add(fhnJetsBeforeCut2);
-
-  // =========== Switch on Sumw2 for all histos ===========
-  for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
-    TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
-    if (h1){
-      h1->Sumw2();
-      continue;
-    }
-       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
-    if (hn){
-      hn->Sumw2();
-    }    
-  }
-  TH1::AddDirectory(oldStatus);
-
-  PostData(1, fOutputList);
+      // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
+      opt = 1<<4 | 1<<5;
+      fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
+   }
+
+   // cent : nInpTrks : rp bins:  deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
+   if(fbJetsDeltaPt){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
+      opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
+      fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
+   }
+
+   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
+   if(fbJetsEta){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
+      opt = 1<<15 | 1<<8 | 1<<9;
+      fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
+   }
+
+   // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
+   if(fbJetsPhi){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
+      opt = 1<<10 | 1<<11;
+      fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
+   }
+
+   // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area)
+   if(fbJetsArea){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26;
+      opt = 1<<18 | 1<<12 | 1<<13;
+      fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
+   }
+
+
+   //before cut
+
+   // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
+   if(fbJetsBeforeCut1){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
+      opt = 0;
+      fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
+   }
+
+   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
+   if(fbJetsBeforeCut2){
+      entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
+      opt = 0;
+      fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
+   }
+
+   fOutputList->Add(fHistEvtSelection);
+   fOutputList->Add(fHistJetSelection);
+   fOutputList->Add(fh2JetSelection);
+   if(fbEvent)           fOutputList->Add(fhnEvent);
+   if(fbJetsMismatch1)   fOutputList->Add(fhnJetsMismatch1);
+   if(fbJetsMismatch2)   fOutputList->Add(fhnJetsMismatch2);
+   if(fbJetsRp)          fOutputList->Add(fhnJetsRp);
+   if(fbJetsDeltaPt)     fOutputList->Add(fhnJetsDeltaPt);
+   if(fbJetsEta)         fOutputList->Add(fhnJetsEta);
+   if(fbJetsPhi)         fOutputList->Add(fhnJetsPhi);
+   if(fbJetsArea)        fOutputList->Add(fhnJetsArea);
+   if(fbJetsBeforeCut1)  fOutputList->Add(fhnJetsBeforeCut1);
+   if(fbJetsBeforeCut2)  fOutputList->Add(fhnJetsBeforeCut2);
+
+   // =========== Switch on Sumw2 for all histos ===========
+   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
+      TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
+      if (h1){
+         h1->Sumw2();
+         continue;
+      }
+      THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
+      if (hn){
+         hn->Sumw2();
+      }          
+   }
+   TH1::AddDirectory(oldStatus);
+
+   PostData(1, fOutputList);
 }
 
 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
 {
-  // load events, apply event cuts, then compare leading jets
+   // load events, apply event cuts, then compare leading jets
 
-  if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
+   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
       AliError("Jet branch name not set.");
       return;
    }
 
-  fESD=dynamic_cast<AliESDEvent*>(InputEvent());
-  if (!fESD) {
-    AliError("ESD not available");
-    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
-  } else {
-    fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
-  }
-  if (!fAOD) {
-    AliError("AOD not available");
-    return;
-  }
-
-  // -- event selection --
-  fHistEvtSelection->Fill(1); // number of events before event selection
-
-  // physics selection
-  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
-    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
-  if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
-    if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
-    fHistEvtSelection->Fill(2);
-    PostData(1, fOutputList);
-    return;
-  }
-
-  // vertex selection
-  AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
-  Int_t nTracksPrim = primVtx->GetNContributors();
-  if ((nTracksPrim < fMinContribVtx) ||
-      (primVtx->GetZ() < fVtxZMin) ||
-      (primVtx->GetZ() > fVtxZMax) ){
-    if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
-    fHistEvtSelection->Fill(3);
-    PostData(1, fOutputList);
-    return;
-  }
-
-  // event class selection (from jet helper task)
-  Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
-  if(fDebug) Printf("Event class %d", eventClass);
-  if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
-    fHistEvtSelection->Fill(4);
-    PostData(1, fOutputList);
-    return;
-  }
-
-  // centrality selection
-  AliCentrality *cent = 0x0;
-  Float_t centValue = 0.; 
-  if(fESD) cent = fESD->GetCentrality();
-  if(cent) centValue = cent->GetCentralityPercentile("V0M");
-  if(fDebug) printf("centrality: %f\n", centValue);
-  if (centValue < fCentMin || centValue > fCentMax){
-    fHistEvtSelection->Fill(4);
-    PostData(1, fOutputList);
-    return;
-  }
-
-
-  // multiplicity due to input tracks
-  Int_t nInputTracks = GetNInputTracks();
-  
-  if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
-       fHistEvtSelection->Fill(5);
-       PostData(1, fOutputList);
-       return;
-  }
-  
+   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
+   if (!fESD) {
+      AliError("ESD not available");
+      fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+   } else {
+      fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
+   }
+   if (!fAOD) {
+      AliError("AOD not available");
+      return;
+   }
+
+   // -- event selection --
+   fHistEvtSelection->Fill(1); // number of events before event selection
+
+   // physics selection
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+   ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
+      if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
+      fHistEvtSelection->Fill(2);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // vertex selection
+   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
+   Int_t nTracksPrim = primVtx->GetNContributors();
+   if ((nTracksPrim < fMinContribVtx) ||
+         (primVtx->GetZ() < fVtxZMin) ||
+         (primVtx->GetZ() > fVtxZMax) ){
+      if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+      fHistEvtSelection->Fill(3);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // event class selection (from jet helper task)
+   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
+   if(fDebug) Printf("Event class %d", eventClass);
+   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
+      fHistEvtSelection->Fill(4);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // centrality selection
+   AliCentrality *cent = 0x0;
+   Float_t centValue = 0.; 
+   if(fESD) cent = fESD->GetCentrality();
+   if(cent) centValue = cent->GetCentralityPercentile("V0M");
+   if(fDebug) printf("centrality: %f\n", centValue);
+   if (centValue < fCentMin || centValue > fCentMax){
+      fHistEvtSelection->Fill(4);
+      PostData(1, fOutputList);
+      return;
+   }
+
+
+   // multiplicity due to input tracks
+   Int_t nInputTracks = GetNInputTracks();
+
+   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
+      fHistEvtSelection->Fill(5);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   
+   fHistEvtSelection->Fill(0); // accepted events  
+   // -- end event selection --
+
+   // pt hard
+   Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
+   Int_t pthardbin = GetPtHardBin(pthard);
+
+   // reaction plane
+   Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
+   if(fbEvent){
+      Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
+      fhnEvent->Fill(eventEntries);
+   }
+
+
+   // get background
+   AliAODJetEventBackground* externalBackground = 0;
+   if(!externalBackground&&fBackgroundBranch.Length()){
+      externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
+      //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
+   }
+   Float_t rho = 0;
+   if(externalBackground)rho = externalBackground->GetBackground(0);
+
+
+   // fetch jets
+   TClonesArray *aodJets[2];
+   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
+   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
+
+
+
+   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
+      fListJets[iJetType]->Clear();
+      if (!aodJets[iJetType]) continue;
+
+      if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
+      
+      for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
+         AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
+         if (jet) fListJets[iJetType]->Add(jet);
+      }
+      fListJets[iJetType]->Sort();
+   }
+   
+   
+   // jet matching 
+   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
+   static TArrayF aPtFraction(fListJets[0]->GetEntries());
+   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
+   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
+
+   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
+   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
+   fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
+   aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
    
-  fHistEvtSelection->Fill(0); // accepted events  
-  // -- end event selection --
-
-  Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
-  if(fbEvent){
-   Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
-   fhnEvent->Fill(eventEntries);
-  }
-  
-  
-  // get background
-  AliAODJetEventBackground* externalBackground = 0;
-  if(!externalBackground&&fBackgroundBranch.Length()){
-    externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
-    //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
-  }
-  Float_t rho = 0;
-  if(externalBackground)rho = externalBackground->GetBackground(0);
-  
-  
-  // fetch jets
-  TClonesArray *aodJets[2];
-  aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
-  aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
-  
-  
-
-  for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
-    fListJets[iJetType]->Clear();
-    if (!aodJets[iJetType]) continue;
-
-    if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
-       
-    for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
-      AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
-      if (jet) fListJets[iJetType]->Add(jet);
-    }
-    fListJets[iJetType]->Sort();
-  }
-    
-    
-  // jet matching 
-  static TArrayI aMatchIndex(fListJets[0]->GetEntries());
-  static TArrayF aPtFraction(fListJets[0]->GetEntries());
-  if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
-  if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
-  
-  // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
-  AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
-                                            fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
-                                            aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
-                                                                                       
-  // loop over matched jets
-  Int_t ir = -1; // index of matched reconstruced jet
-  Float_t fraction = -1.;
-  AliAODJet *jet[2]  = { 0x0, 0x0 };
-  Float_t jetEta[2]  = { -990., -990. };
-  Float_t jetPhi[2]  = { -990., -990. };
-  Float_t jetPt[2]   = { -990., -990. };
-  Float_t jetArea[2] = { -990., -990. };
-  Float_t rpJet[2]   = { -990., -990. };
+   // loop over matched jets
+   Int_t ir = -1; // index of matched reconstruced jet
+   Float_t fraction = -1.;
+   AliAODJet *jet[2]  = { 0x0, 0x0 };
+   Float_t jetEta[2]  = { -990., -990. };
+   Float_t jetPhi[2]  = { -990., -990. };
+   Float_t jetPt[2]   = { -990., -990. };
+   Float_t jetArea[2] = { -990., -990. };
+   Float_t rpJet[2]   = { -990., -990. };
+   Int_t rpBin = -990;
    
-  for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
-    fHistJetSelection->Fill(1); // all probe jets
-    ir = aMatchIndex[ig];
-        if(ir<0){
-      fHistJetSelection->Fill(2);
-      
-      if(fbJetsMismatch1){
-         jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
-         if(!jet[0]) continue;
-         jetEta[0]  = jet[0]->Eta();
-         jetPhi[0]  = jet[0]->Phi();
-         jetPt[0]   = jet[0]->Pt();
-         jetArea[0] = jet[0]->EffectiveAreaCharged();
-         rpJet[0]   = TVector2::Phi_mpi_pi(rp-jetPhi[0]);
-      
-         Int_t rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[0]), 3);
-      
-         Double_t jetEntriesMismatch1[7] = {
-           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
-           (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
-         };             
-         fhnJetsMismatch1->Fill(jetEntriesMismatch1);
-      }
-      
-      continue;
-        }
-        fraction = aPtFraction[ig];
-        
-        // fetch jets
-        jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
-        jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
-        if(!jet[0] || !jet[1]){
-      fHistJetSelection->Fill(3);
-      continue;
-        }
-        
-        // look for distance to next rec jet
-        Float_t distNextJet = -0.01; // no neighbor
-        Float_t ptNextJet = -1.;
-        for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
-           if(r==ir) continue;
-               Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
-               if(distNextJet<0. || distNextJet>tmpDeltaR){
-               distNextJet = tmpDeltaR;
-                   ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
-               }
-        }
-        
-        // get parameters
-        for(Int_t i=0; i<fkNbranches; ++i){
-      jetEta[i]  = jet[i]->Eta();
-      jetPhi[i]  = jet[i]->Phi();
-      jetPt[i]   = jet[i]->Pt();
-      jetArea[i] = jet[i]->EffectiveAreaCharged();
-      rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
-        }
-        Int_t rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[1]), 3);
-    //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
-    //Float_t relRho = rho>0. ? localRho / rho : 0.;
-
-        
-        // calculate parameters of associated jets
-    /* from Leticia, kT clusterizer
-    Float_t par0[4] = { 0.00409,       0.01229,      0.05,         0.26 };
-    Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
-    Float_t par2[4] = { 1.02865e-01,   1.49039e-01,  1.53910e-01,  1.51109e-01 };
-    */
-    // own, from embedded tracks
-    Float_t par0[4] = { 0.02841,       0.05039,      0.09092,      0.24089     };
-    Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
-    Float_t par2[4] = { 4.95415e-02,   9.79538e-02,  1.32814e-01,  1.71743e-01 };
-    
-    Float_t rpCorr = 0.;         
-    
-    if(eventClass>0&&eventClass<4){
-       rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
-       rpCorr *= rho * jetArea[1];
-    }
-
-        Float_t deltaPt    = jetPt[1]-jetPt[0];
-        Float_t deltaEta   = jetEta[1]-jetEta[0];
-        Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
-        Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
-        Float_t deltaArea  = jetArea[1]-jetArea[0];
-        
-       
-        // fill thnsparse before acceptance cut
-    if(fbJetsBeforeCut1){
-      Double_t jetBeforeCutEntries1[10] = { 
-         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
-         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
-         (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
-         };
-      fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
-    }
-        
-    if(fbJetsBeforeCut2){
-      Double_t jetBeforeCutEntries2[10] = {
-         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, 
-         (Double_t)jetPt[0], (Double_t)jetPt[1],
-         (Double_t)jetEta[0], (Double_t)jetEta[1], 
-         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
-         };
-          fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
-    }
-        
-        Bool_t jetAccepted = kTRUE;
-        // minimum fraction required
-        if(fraction<fJetPtFractionMin){
-       fHistJetSelection->Fill(4);
-       jetAccepted = kFALSE;
-        }
-        
-    if(jetAccepted){
-      // jet acceptance + minimum pT check
-      if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
-         jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
-               
-         if(fDebug){
-            Printf("Jet not in eta acceptance.");
-            Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
-            Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
+   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
+      ir = aMatchIndex[ig];
+
+      //fetch jets
+      jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
+      if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
+      else      jet[1] = 0x0;
+
+      for(Int_t i=0; i<fkNbranches; ++i){
+         if(!jet[i]){
+            jetEta[i]  = -990;
+            jetPhi[i]  = -990.;
+            jetPt[i]   = -990.;
+            jetArea[i] = -990.;
+            rpJet[i]   = -990.;
+            if(i==1) rpBin = -990;
+         } 
+         else {
+            jetEta[i]  = jet[i]->Eta();
+            jetPhi[i]  = jet[i]->Phi();
+            jetPt[i]   = GetPt(jet[i], i);
+            jetArea[i] = jet[i]->EffectiveAreaCharged();
+            rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
+            if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
          }
-         fHistJetSelection->Fill(5);
-         jetAccepted = kFALSE;
       }
-    }
-    if(jetAccepted){
-      if(jetPt[1] < fJetPtMin){
-         if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
-         fHistJetSelection->Fill(6);
-         jetAccepted = kFALSE;
+      fraction = aPtFraction[ig];
+
+      // jet statistics
+      fHistJetSelection->Fill(1); // all probe jets
+      if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
+
+      if(ir<0){
+         fHistJetSelection->Fill(2);
+         if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
+         
+         if(fbJetsMismatch1){
+            if(!jet[0]) continue;
+            Double_t jetEntriesMismatch1[7] = {
+               (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+               (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
+            };          
+            fhnJetsMismatch1->Fill(jetEntriesMismatch1);
+         }
+         
+         continue;
+      }
+      
+      if(!jet[0] || !jet[1]){
+         fHistJetSelection->Fill(3);
+         if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
+         continue;
+      }
+      
+      // look for distance to next rec jet
+      Float_t distNextJet = -0.01; // no neighbor
+      Float_t ptNextJet = -1.;
+      for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
+         if(r==ir) continue;
+         Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
+         if(distNextJet<0. || distNextJet>tmpDeltaR){
+            distNextJet = tmpDeltaR;
+            if(fKeepJets) ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
+            else          ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
+         }
       }
-    }
-    
-    if(!jetAccepted){
-      if(fbJetsMismatch2){
-         Double_t jetEntriesMismatch2[10] = {
+      
+
+
+      //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
+      //Float_t relRho = rho>0. ? localRho / rho : 0.;
+
+      
+      // calculate parameters of associated jets
+      /* from Leticia, kT clusterizer
+   Float_t par0[4] = { 0.00409,       0.01229,      0.05,         0.26 };
+   Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
+   Float_t par2[4] = { 1.02865e-01,   1.49039e-01,  1.53910e-01,  1.51109e-01 };
+   */
+      // own, from embedded tracks
+      Float_t par0[4] = { 0.02841,       0.05039,      0.09092,      0.24089     };
+      Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
+      Float_t par2[4] = { 4.95415e-02,   9.79538e-02,  1.32814e-01,  1.71743e-01 };
+      
+      Float_t rpCorr = 0.;         
+      
+      if(eventClass>0&&eventClass<4){
+         rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
+         rpCorr *= rho * jetArea[1];
+      }
+
+      Float_t deltaPt    = jetPt[1]-jetPt[0];
+      Float_t deltaEta   = jetEta[1]-jetEta[0];
+      Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
+      Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
+      Float_t deltaArea  = jetArea[1]-jetArea[0];
+      
+      
+      // fill thnsparse before acceptance cut
+      if(fbJetsBeforeCut1){
+         Double_t jetBeforeCutEntries1[10] = { 
             (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+            (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
+            (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
+         };
+         fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
+      }
+      
+      if(fbJetsBeforeCut2){
+         Double_t jetBeforeCutEntries2[10] = {
+            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, 
             (Double_t)jetPt[0], (Double_t)jetPt[1],
-            (Double_t)jetEta[0], (Double_t)jetEta[1],
-            (Double_t)deltaEta, (Double_t)deltaR,
-            (Double_t)fraction
+            (Double_t)jetEta[0], (Double_t)jetEta[1], 
+            (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
          };
-         fhnJetsMismatch2->Fill(jetEntriesMismatch2);     
-      }
-      continue;
-    }
-    
-    // all accepted jets
-    fHistJetSelection->Fill(0);
-        
-        // fill thnsparse
-    if(fbJetsRp){
-      Double_t jetEntriesRp[10] = {
-          (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
-          (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], 
-          (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr
-          };
-      fhnJetsRp->Fill(jetEntriesRp);
-    }
-        
-    if(fbJetsDeltaPt){
-      Double_t jetEntriesDeltaPt[7] = {
-           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
-           (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea
-           };           
-      fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
-    }
-        
-    if(fbJetsEta){
-      Double_t jetEntriesEta[10] = {
-           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
-           (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
-           (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
-           };                           
-      fhnJetsEta->Fill(jetEntriesEta);
-    }
-    
-    if(fbJetsPhi){
-      Double_t jetEntriesPhi[9] = {
+         fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
+      }
+      
+      // jet selection
+      Bool_t jetAccepted = kTRUE;
+      // minimum fraction required
+      if(fraction<fJetPtFractionMin){
+         fHistJetSelection->Fill(4);
+         fh2JetSelection->Fill(4.,jetPt[0]);
+         jetAccepted = kFALSE;
+      }
+      
+      if(jetAccepted){
+         // jet acceptance + minimum pT check
+         if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
+               jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
+            
+            if(fDebug){
+               Printf("Jet not in eta acceptance.");
+               Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
+               Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
+            }
+            fHistJetSelection->Fill(5);
+            fh2JetSelection->Fill(5.,jetPt[0]);
+            jetAccepted = kFALSE;
+         }
+      }
+      if(jetAccepted){
+         if(jetPt[1] < fJetPtMin){
+            if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
+            fHistJetSelection->Fill(6);
+            fh2JetSelection->Fill(6.,jetPt[0]);
+            jetAccepted = kFALSE;
+         }
+      }
+
+      if(jetAccepted){
+         if(jet[1]->Trigger()&fJetTriggerExcludeMask){
+            fHistJetSelection->Fill(7);
+            fh2JetSelection->Fill(7.,jetPt[0]);
+            jetAccepted = kFALSE;
+         }
+         
+      }
+      
+      if(!jetAccepted){
+         if(fbJetsMismatch2){
+            Double_t jetEntriesMismatch2[10] = {
+               (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+               (Double_t)jetPt[0], (Double_t)jetPt[1],
+               (Double_t)jetEta[0], (Double_t)jetEta[1],
+               (Double_t)deltaEta, (Double_t)deltaR,
+               (Double_t)fraction
+            };
+            fhnJetsMismatch2->Fill(jetEntriesMismatch2);     
+         }
+         continue;
+      }
+      
+      // all accepted jets
+      fHistJetSelection->Fill(0);
+      fh2JetSelection->Fill(0.,jetPt[0]);
+      
+      // fill thnsparse
+      if(fbJetsRp){
+         Double_t jetEntriesRp[11] = {
+            (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
+            (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], 
+            (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
+            (Double_t)pthardbin
+         };
+         fhnJetsRp->Fill(jetEntriesRp);
+      }
+      
+      if(fbJetsDeltaPt){
+         Double_t jetEntriesDeltaPt[8] = {
+            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+            (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
+            (Double_t)pthardbin
+         };             
+         fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
+      }
+      
+      if(fbJetsEta){
+         Double_t jetEntriesEta[11] = {
+            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+            (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
+            (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
+            (Double_t)pthardbin
+         };                             
+         fhnJetsEta->Fill(jetEntriesEta);
+      }
+      
+      if(fbJetsPhi){
+         Double_t jetEntriesPhi[10] = {
             (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
             (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
-            (Double_t)deltaPt, (Double_t)deltaPhi
-      };
-      fhnJetsPhi->Fill(jetEntriesPhi);
-    }
-        
-    if(fbJetsArea){
-      Double_t jetEntriesArea[13] = {
-           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
-           (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1], 
-           (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea, 
-           (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet
-           };                           
-      fhnJetsArea->Fill(jetEntriesArea);
-    }
-        
-  }
-
-  PostData(1, fOutputList);
+            (Double_t)deltaPt, (Double_t)deltaPhi,
+            (Double_t)pthardbin
+         };
+         fhnJetsPhi->Fill(jetEntriesPhi);
+      }
+      
+      if(fbJetsArea){
+         Double_t jetEntriesArea[14] = {
+            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
+            (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1], 
+            (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea, 
+            (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
+            (Double_t)pthardbin
+         };                             
+         fhnJetsArea->Fill(jetEntriesArea);
+      }
+      
+   }
+
+   PostData(1, fOutputList);
 }
 
 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
 {
-  // Draw result to the screen
-  // Called once at the end of the query
+   // Draw result to the screen
+   // Called once at the end of the query
 
-  if (!GetOutputData(1))
-    return;
+   if (!GetOutputData(1))
+   return;
 }
 
 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
 {
 
-  Int_t nInputTracks = 0;
-  
-  TString jbname(fJetBranchName[1]);
-  //needs complete event, use jets without background subtraction
-  for(Int_t i=1; i<=3; ++i){
-    if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
-  }
-  // use only HI event
-  if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
-  if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
-  
-  if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
-  TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
-  if(!tmpAODjets){
-    Printf("Jet branch %s not found", jbname.Data());
-       Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
-       return -1;
-  }
-      
-  for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
+   Int_t nInputTracks = 0;
+
+   TString jbname(fJetBranchName[1]);
+   //needs complete event, use jets without background subtraction
+   for(Int_t i=1; i<=3; ++i){
+      if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
+   }
+   // use only HI event
+   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
+   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
+
+   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
+   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
+   if(!tmpAODjets){
+      Printf("Jet branch %s not found", jbname.Data());
+      Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
+      return -1;
+   }
+   
+   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
-         if(!jet) continue;
-         TRefArray *trackList = jet->GetRefTracks();
-         Int_t nTracks = trackList->GetEntriesFast();
-         nInputTracks += nTracks;
-         if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
-  }
-  if(fDebug) Printf("---> input tracks: %d", nInputTracks);
-  
-  return nInputTracks;  
+      if(!jet) continue;
+      TRefArray *trackList = jet->GetRefTracks();
+      Int_t nTracks = trackList->GetEntriesFast();
+      nInputTracks += nTracks;
+      if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
+   }
+   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
+
+   return nInputTracks;  
 }
 
 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
 {
-  Int_t count = 0;
-  UInt_t tmp = entries;
-  while(tmp!=0){
-     count++;
-     tmp = tmp &~ -tmp;  // clear lowest bit
-  }
-
-  TString hnTitle(name);
-  const Int_t dim = count;
-  Int_t nbins[dim];
-  Double_t xmin[dim];
-  Double_t xmax[dim];
-  
-  Int_t i=0;
-  Int_t c=0;
-  while(c<dim && i<32){
+   Int_t count = 0;
+   UInt_t tmp = entries;
+   while(tmp!=0){
+      count++;
+      tmp = tmp &~ -tmp;  // clear lowest bit
+   }
+
+   TString hnTitle(name);
+   const Int_t dim = count;
+   Int_t nbins[dim];
+   Double_t xmin[dim];
+   Double_t xmax[dim];
+
+   Int_t i=0;
+   Int_t c=0;
+   while(c<dim && i<32){
       if(entries&(1<<i)){
          Bool_t highres = opt&(1<<i);
          TString label("");
@@ -718,12 +774,12 @@ THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t
          hnTitle += Form(";%s",label.Data());
          c++;
       }
-     
-     i++;
-  }
-  hnTitle += ";";
-  
-  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
+      
+      i++;
+   }
+   hnTitle += ";";
+
+   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 }
 
 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
@@ -732,267 +788,304 @@ void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString
    const Double_t pi = TMath::Pi();
    
    switch(iEntry){
-   
-      case 0:
-         label = "V0 centrality (%)";
-         if(hr){
-            nbins = 100;
-            xmin = 0.;
-            xmax = 100.;
-         } else {
-            nbins = 8;
-            xmin = 0.;
-            xmax = 80.;
-         }
+      
+   case 0:
+      label = "V0 centrality (%)";
+      if(hr){
+         nbins = 100;
+         xmin = 0.;
+         xmax = 100.;
+      } else {
+         nbins = 8;
+         xmin = 0.;
+         xmax = 80.;
+      }
       break;
       
       
-      case 1:
-         label = "nb. of input tracks";
-         if(fIsPbPb){
-            if(hr){
-               nbins = 400;
-               xmin = 0.;
-               xmax = 4000.;
-            } else {
-               nbins = 40;
-               xmin = 0.;
-               xmax = 4000.;
-            }
+   case 1:
+      label = "nb. of input tracks";
+      if(fIsPbPb){
+         if(hr){
+            nbins = 400;
+            xmin = 0.;
+            xmax = 4000.;
          } else {
             nbins = 40;
             xmin = 0.;
-            xmax = 400.;
+            xmax = 4000.;
          }
+      } else {
+         nbins = 40;
+         xmin = 0.;
+         xmax = 400.;
+      }
       break;
       
       
-      case 2:
-         label = "event plane #psi";
-         if(hr){
-            nbins = 30;
-            xmin = 0.;
-            xmax = pi;
-         } else {
-            nbins = 30;
-            xmin = 0.;
-            xmax = pi;
-         }
+   case 2:
+      label = "event plane #psi";
+      if(hr){
+         nbins = 30;
+         xmin = 0.;
+         xmax = pi;
+      } else {
+         nbins = 30;
+         xmin = 0.;
+         xmax = pi;
+      }
       break;
       
       
-      case 3:
-         label = "event plane bin";
-        nbins = 3;
-         xmin = -.5;
-         xmax = 2.5;
+   case 3:
+      label = "event plane bin";
+      nbins = 3;
+      xmin = -.5;
+      xmax = 2.5;
       break;
       
       
-      case 4:
-      case 5:
-         if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
-         if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
-         nbins = 48;
-         xmin = -pi;
-         xmax =  pi;
+   case 4:
+   case 5:
+      if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
+      if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
+      nbins = 48;
+      xmin = -pi;
+      xmax =  pi;
       break;
       
       
-      case 6:
-      case 7:
-         if(iEntry==6)label = "probe p_{T} (GeV/c)";
-         if(iEntry==7)label = "rec p_{T} (GeV/c)";
-         if(hr){
-            nbins = 250;
-            xmin = 0.;
-            xmax = 250.;
-         } else {
-            nbins = 50;
-            xmin = 0.;
-            xmax = 250.;
-         }
+   case 6:
+   case 7:
+      if(iEntry==6)label = "probe p_{T} (GeV/c)";
+      if(iEntry==7)label = "rec p_{T} (GeV/c)";
+      if(hr){
+         nbins = 300;
+         xmin = -50.;
+         xmax = 250.;
+      } else {
+         nbins = 50;
+         xmin = 0.;
+         xmax = 250.;
+      }
       break;
       
       
-      case 8:
-      case 9:
-         if(iEntry==8)label = "probe #eta";
-         if(iEntry==9)label = "rec #eta";
-         if(hr){
-            nbins = 56;
-            xmin = -.7;
-            xmax =  .7;
-         } else {
-            nbins = 28;
-            xmin = -.7;
-            xmax =  .7;
-         }
+   case 8:
+   case 9:
+      if(iEntry==8)label = "probe #eta";
+      if(iEntry==9)label = "rec #eta";
+      if(hr){
+         nbins = 56;
+         xmin = -.7;
+         xmax =  .7;
+      } else {
+         nbins = 28;
+         xmin = -.7;
+         xmax =  .7;
+      }
       break;
       
       
-      case 10:
-      case 11:
-         if(iEntry==10)label = "probe #phi";
-         if(iEntry==11)label = "rec #phi";
-         if(hr){
-            nbins = 90;  // modulo 18 (sectors)
-            xmin = 0.;
-            xmax = 2*pi;
-         } else {
-            nbins = 25;
-            xmin = 0.;
-            xmax = 2*pi;
-         }
+   case 10:
+   case 11:
+      if(iEntry==10)label = "probe #phi";
+      if(iEntry==11)label = "rec #phi";
+      if(hr){
+         nbins = 90;  // modulo 18 (sectors)
+         xmin = 0.;
+         xmax = 2*pi;
+      } else {
+         nbins = 25;
+         xmin = 0.;
+         xmax = 2*pi;
+      }
       break;
       
       
-      case 12:
-      case 13:
-         if(iEntry==12)label = "probe area";
-         if(iEntry==13)label = "rec area";
-         if(hr){
-            nbins = 100;
-            xmin = 0.;
-            xmax = 1.;
-         } else {
-            nbins = 25;
-            xmin = 0.;
-            xmax = 1.;
-         }
+   case 12:
+   case 13:
+      if(iEntry==12)label = "probe area";
+      if(iEntry==13)label = "rec area";
+      if(hr){
+         nbins = 100;
+         xmin = 0.;
+         xmax = 1.;
+      } else {
+         nbins = 25;
+         xmin = 0.;
+         xmax = 1.;
+      }
       break;
       
-      case 14:
-         label = "#Delta p_{T}";
-         if(hr){
-            nbins = 241;
-            xmin = -120.5;
-            xmax =  120.5;
-         } else {
-            nbins = 101;
-            xmin = -101.;
-            xmax =  101.;
-         }
+   case 14:
+      label = "#Delta p_{T}";
+      if(hr){
+         nbins = 241;
+         xmin = -120.5;
+         xmax =  120.5;
+      } else {
+         nbins = 101;
+         xmin = -101.;
+         xmax =  101.;
+      }
       break;
       
-      case 15:
-         label = "#Delta#eta";
-         if(hr){
-            nbins = 51;
-            xmin = -1.02;
-            xmax =  1.02;
-         } else {
-            nbins = 51;
-            xmin = -1.02;
-            xmax =  1.02;
-         }
+   case 15:
+      label = "#Delta#eta";
+      if(hr){
+         nbins = 51;
+         xmin = -1.02;
+         xmax =  1.02;
+      } else {
+         nbins = 51;
+         xmin = -1.02;
+         xmax =  1.02;
+      }
       break;
       
       
-      case 16:
-         label = "#Delta#phi";
-         if(hr){
-            nbins = 45;
-            xmin = -pi;
-            xmax =  pi;
-         } else {
-            nbins = 45;
-            xmin = -pi;
-            xmax =  pi;
-         }
+   case 16:
+      label = "#Delta#phi";
+      if(hr){
+         nbins = 45;
+         xmin = -pi;
+         xmax =  pi;
+      } else {
+         nbins = 45;
+         xmin = -pi;
+         xmax =  pi;
+      }
       break;
       
       
-      case 17:
-         label = "#DeltaR";
-         if(hr){
-            nbins = 50;
-            xmin = 0.;
-            xmax = 1.;
-         } else {
-            nbins = 50;
-            xmin = 0.;
-            xmax = 1.;
-         }
+   case 17:
+      label = "#DeltaR";
+      if(hr){
+         nbins = 50;
+         xmin = 0.;
+         xmax = 1.;
+      } else {
+         nbins = 50;
+         xmin = 0.;
+         xmax = 1.;
+      }
       break;
       
       
-      case 18:
-         label = "#Deltaarea";
-         if(hr){
-            nbins = 81;
-            xmin = -.81;
-            xmax =  .81;
-         } else {
-            nbins = 33;
-            xmin = -.825;
-            xmax =  .825;
-         }
+   case 18:
+      label = "#Deltaarea";
+      if(hr){
+         nbins = 81;
+         xmin = -.81;
+         xmax =  .81;
+      } else {
+         nbins = 33;
+         xmin = -.825;
+         xmax =  .825;
+      }
       break;
       
       
-      case 19:
-         label = "fraction";
-         if(hr){
-            nbins = 52;
-            xmin = 0.;
-            xmax = 1.04;
-         } else {
-            nbins = 52;
-            xmin = 0.;
-            xmax = 1.04;
-         }
+   case 19:
+      label = "fraction";
+      if(hr){
+         nbins = 52;
+         xmin = 0.;
+         xmax = 1.04;
+      } else {
+         nbins = 52;
+         xmin = 0.;
+         xmax = 1.04;
+      }
       break;
       
       
-      case 20:
-         label = "distance to closest rec jet";
-         if(hr){
-            nbins = 51;
-            xmin = -0.02;
-            xmax =  1.;
-         } else {
-            nbins = 51;
-            xmin = -0.02;
-            xmax = 1.;
-         }
+   case 20:
+      label = "distance to closest rec jet";
+      if(hr){
+         nbins = 51;
+         xmin = -0.02;
+         xmax =  1.;
+      } else {
+         nbins = 51;
+         xmin = -0.02;
+         xmax = 1.;
+      }
       break;
       
       
-      case 21:
-         label = "p_{T} of closest rec jet";
-         nbins = 100;
-         xmin =  0.;
-         xmax =  200.;
+   case 21:
+      label = "p_{T} of closest rec jet";
+      nbins = 100;
+      xmin =  0.;
+      xmax =  200.;
       break;
       
-      case 22:
-         label = "#rho";
-         nbins = 125;
-         xmin  = 0.;
-         xmax  = 250.;
+   case 22:
+      label = "#rho";
+      nbins = 125;
+      xmin  = 0.;
+      xmax  = 250.;
       break;
       
-      case 23:
-         label = "abs. correction of #rho for RP";
-         nbins =  51;
-         xmin  = -51.;
-         xmax  =  51.;
+   case 23:
+      label = "abs. correction of #rho for RP";
+      nbins =  51;
+      xmin  = -51.;
+      xmax  =  51.;
       break;
       
-      case 24:
-         label = "local #rho";
-         nbins = 125;
-         xmin  = 0.;
-         xmax  = 250.;
+   case 24:
+      label = "local #rho";
+      nbins = 125;
+      xmin  = 0.;
+      xmax  = 250.;
       break;
       
-      case 25:
-         label = "local #rho / #rho";
-         nbins = 500;
-         xmin  = 0.;
-         xmax  = 5.;
+   case 25:
+      label = "local #rho / #rho";
+      nbins = 500;
+      xmin  = 0.;
+      xmax  = 5.;
+      break;
+      
+   case 26:
+      label = "p_{T,hard} bin";
+      nbins = 10;
+      xmin  = -.5;
+      xmax  = 9.5;
+      break;
+      
+   }
+
+}
+
+//____________________________________________________________________________
+Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
+
+   const Int_t nBins = 10;
+   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
    
+   Int_t bin = -1;
+   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
+      bin++;
+   }
+   
+   return bin;
+}
+
+//____________________________________________________________________________
+Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
+
+   Double_t pt = 0.;
+
+   if(fKeepJets && mode==1){  // background subtracted pt, can be negative
+      pt = j->GetPtSubtracted(0);
+   }
+   else{
+      pt = j->Pt();  // if negative, pt=0.01
    }
 
+   return pt;
 }
index 1383c6a..b847223 100644 (file)
@@ -12,130 +12,138 @@ class AliAODEvent;
 #include "AliVEvent.h"
 
 class AliAnalysisTaskJetResponseV2 : public AliAnalysisTaskSE {
- public:
-  AliAnalysisTaskJetResponseV2();
-  AliAnalysisTaskJetResponseV2(const char *name);
-  virtual ~AliAnalysisTaskJetResponseV2();
-
-  virtual void     LocalInit() {Init();}
-  virtual void     Init();
-  virtual void     UserCreateOutputObjects();
-  virtual void     UserExec(Option_t *option);
-  virtual void     Terminate(const Option_t*);
-
-  virtual Int_t      GetNInputTracks();
-  virtual THnSparse* NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt);
-  virtual void       GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax);
-
-  virtual AliVEvent::EOfflineTriggerTypes GetOfflineTrgMask() const { return fOfflineTrgMask; }
-  virtual void     GetBranchNames(TString &branch1, TString &branch2) const { branch1 = fJetBranchName[0]; branch2 = fJetBranchName[1]; }
-  virtual Bool_t   GetIsPbPb() const { return fIsPbPb; }
-  virtual Int_t    GetMinContribVtx() const { return fMinContribVtx; };
-  virtual Float_t  GetVtxZMin() const { return fVtxZMin; }
-  virtual Float_t  GetVtxZMax() const { return fVtxZMax; }
-  virtual Int_t    GetEvtClassMin() const { return fEvtClassMin; }
-  virtual Int_t    GetEvtClassMax() const { return fEvtClassMax; }
-  virtual Float_t  GetCentMin() const { return fCentMin; }
-  virtual Float_t  GetCentMax() const { return fCentMax; }
-  virtual Int_t    GetNInputTracksMin() const { return fNInputTracksMin; }
-  virtual Int_t    GetNInputTracksMax() const { return fNInputTracksMax; } 
-  virtual Float_t  GetJetEtaMin() const { return fJetEtaMin; }
-  virtual Float_t  GetJetEtaMax() const { return fJetEtaMax; }
-  virtual Float_t  GetJetPtMin() const { return fJetPtMin; }
-  virtual Float_t  GetJetPtFractionMin() const { return fJetPtFractionMin; }
-  virtual Int_t    GetNMatchJets() const { return fNMatchJets; }
-
-  virtual void     SetBranchNames(const TString &branch1, const TString &branch2);
-  virtual void     SetBackgroundBranch(TString &branch) { fBackgroundBranch = branch;}
-  virtual void     SetIsPbPb(Bool_t b=kTRUE) { fIsPbPb = b; }
-  virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
-  virtual void     SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
-  virtual void     SetVtxZMin(Float_t z) { fVtxZMin = z; }
-  virtual void     SetVtxZMax(Float_t z) { fVtxZMax = z; }
-  virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
-  virtual void     SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; }
-  virtual void     SetCentMin(Float_t cent) { fCentMin = cent; }
-  virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
-  virtual void     SetNInputTracksMin(Int_t nTr) { fNInputTracksMin = nTr; }
-  virtual void     SetNInputTracksMax(Int_t nTr) { fNInputTracksMax = nTr; }
-  virtual void     SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; }
-  virtual void     SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; }
-  virtual void     SetJetPtMin(Float_t pt) { fJetPtMin = pt; }
-  virtual void     SetJetPtFractionMin(Float_t frac) { fJetPtFractionMin = frac; }
-  virtual void     SetNMatchJets(Int_t n) { fNMatchJets = n; }
-  virtual void     SetFillEvent(Bool_t b) { fbEvent = b; }
-  virtual void     SetFillJetsMismatch1(Bool_t b) { fbJetsMismatch1 = b; }
-  virtual void     SetFillJetsMismatch2(Bool_t b) { fbJetsMismatch2 = b; }
-  virtual void     SetFillJetsRp(Bool_t b) { fbJetsRp = b; }
-  virtual void     SetFillJetsDeltaPt(Bool_t b) { fbJetsDeltaPt = b; }
-  virtual void     SetFillJetsEta(Bool_t b) { fbJetsEta = b; }
-  virtual void     SetFillJetsPhi(Bool_t b) { fbJetsPhi = b; }
-  virtual void     SetFillJetsArea(Bool_t b) { fbJetsArea = b; }
-  virtual void     SetFillJetsBeforeCut1(Bool_t b) { fbJetsBeforeCut1 = b; }
-  virtual void     SetFillJetsBeforeCut2(Bool_t b) { fbJetsBeforeCut2 = b; }
-  
- private:
-  // ESD/AOD events
-  AliESDEvent *fESD;    //! ESD object
-  AliAODEvent *fAOD;    //! AOD event
-
-  // jets to compare
-  TString fJetBranchName[2]; //  name of jet branches to compare
-  TList *fListJets[2];       //! jet lists
-  
-  TString fBackgroundBranch;
-
-  // event selection
-  Bool_t fIsPbPb;         // is Pb-Pb (fast embedding) or p-p (detector response)
-  AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept
-  Int_t   fMinContribVtx; // minimum number of track contributors for primary vertex
-  Float_t fVtxZMin;      // lower bound on vertex z
-  Float_t fVtxZMax;      // upper bound on vertex z
-  Int_t   fEvtClassMin;          // lower bound on event class
-  Int_t   fEvtClassMax;          // upper bound on event class
-  Float_t fCentMin;      // lower bound on centrality
-  Float_t fCentMax;      // upper bound on centrality
-  Int_t   fNInputTracksMin;  // lower bound of nb. of input tracks
-  Int_t   fNInputTracksMax;  // upper bound of nb. of input tracks
-  Float_t fJetEtaMin;     // lower bound on eta for found jets
-  Float_t fJetEtaMax;     // upper bound on eta for found jets
-  Float_t fJetPtMin;      // minimum jet pT
-  Float_t fJetPtFractionMin; // minimum fraction for positiv match of jets
-  Int_t   fNMatchJets;       // maximal nb. of jets taken for matching
-  Double_t fMatchMaxDist;     // maximal distance of matching jets
-  
-
-  // output objects
-  const Int_t fkNbranches;                   //! number of branches to be read
-  const Int_t fkEvtClasses;                  //! number of event classes
-  TList *fOutputList;                        //! output data container
-  Bool_t fbEvent;                            // fill fhnEvent
-  Bool_t fbJetsMismatch1;                    // fill fhnJetsMismatch1
-  Bool_t fbJetsMismatch2;                    // fill fhnJetsMismatch2
-  Bool_t fbJetsRp;                           // fill fhnJetsRp
-  Bool_t fbJetsDeltaPt;                      // fill fhnJetsDeltaPt
-  Bool_t fbJetsEta;                          // fill fhnJetsEta
-  Bool_t fbJetsPhi;                          // fill fhnJetsEta
-  Bool_t fbJetsArea;                         // fill fhnJetsArea
-  Bool_t fbJetsBeforeCut1;                   // fill fhnJetsBeforeCut1
-  Bool_t fbJetsBeforeCut2;                   // fill fhnJetsBeforeCut2
-  TH1I  *fHistEvtSelection;                  //! event selection statistic
-  TH1I  *fHistJetSelection;                  //! jet selection statistic
-  THnSparse *fhnEvent;                       //! variables per event
-  THnSparse *fhnJetsMismatch1;               //! variables per jet
-  THnSparse *fhnJetsMismatch2;               //! variables per jet
-  THnSparse *fhnJetsRp;                      //! variables per jet
-  THnSparse *fhnJetsDeltaPt;                 //! variables per jet
-  THnSparse *fhnJetsEta;                     //! variables per jet
-  THnSparse *fhnJetsPhi;                     //! variables per jet
-  THnSparse *fhnJetsArea;                    //! variables per jet
-  THnSparse *fhnJetsBeforeCut1;               //! variables per jet before acceptance cut
-  THnSparse *fhnJetsBeforeCut2;               //! variables per jet before acceptance cut
-  
-  AliAnalysisTaskJetResponseV2(const AliAnalysisTaskJetResponseV2&); // not implemented
-  AliAnalysisTaskJetResponseV2& operator=(const AliAnalysisTaskJetResponseV2&); // not implemented
-
-  ClassDef(AliAnalysisTaskJetResponseV2, 3);
+public:
+   AliAnalysisTaskJetResponseV2();
+   AliAnalysisTaskJetResponseV2(const char *name);
+   virtual ~AliAnalysisTaskJetResponseV2();
+
+   virtual void     LocalInit() {Init();}
+   virtual void     Init();
+   virtual void     UserCreateOutputObjects();
+   virtual void     UserExec(Option_t *option);
+   virtual void     Terminate(const Option_t*);
+
+   virtual Int_t      GetNInputTracks();
+   virtual THnSparse* NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt);
+   virtual void       GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax);
+   virtual Int_t      GetPtHardBin(Double_t ptHard);
+   virtual Double_t   GetPt(AliAODJet* j, Int_t mode);
+
+   virtual AliVEvent::EOfflineTriggerTypes GetOfflineTrgMask() const { return fOfflineTrgMask; }
+   virtual void     GetBranchNames(TString &branch1, TString &branch2) const { branch1 = fJetBranchName[0]; branch2 = fJetBranchName[1]; }
+   virtual Bool_t   GetIsPbPb() const { return fIsPbPb; }
+   virtual Int_t    GetMinContribVtx() const { return fMinContribVtx; };
+   virtual Float_t  GetVtxZMin() const { return fVtxZMin; }
+   virtual Float_t  GetVtxZMax() const { return fVtxZMax; }
+   virtual Int_t    GetEvtClassMin() const { return fEvtClassMin; }
+   virtual Int_t    GetEvtClassMax() const { return fEvtClassMax; }
+   virtual Float_t  GetCentMin() const { return fCentMin; }
+   virtual Float_t  GetCentMax() const { return fCentMax; }
+   virtual Int_t    GetNInputTracksMin() const { return fNInputTracksMin; }
+   virtual Int_t    GetNInputTracksMax() const { return fNInputTracksMax; } 
+   virtual Float_t  GetJetEtaMin() const { return fJetEtaMin; }
+   virtual Float_t  GetJetEtaMax() const { return fJetEtaMax; }
+   virtual Float_t  GetJetPtMin() const { return fJetPtMin; }
+   virtual Float_t  GetJetPtFractionMin() const { return fJetPtFractionMin; }
+   virtual Int_t    GetNMatchJets() const { return fNMatchJets; }
+
+   virtual void     SetBranchNames(const TString &branch1, const TString &branch2);
+   virtual void     SetBackgroundBranch(TString &branch) { fBackgroundBranch = branch;}
+   virtual void     SetIsPbPb(Bool_t b=kTRUE) { fIsPbPb = b; }
+   virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+   virtual void     SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
+   virtual void     SetVtxZMin(Float_t z) { fVtxZMin = z; }
+   virtual void     SetVtxZMax(Float_t z) { fVtxZMax = z; }
+   virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
+   virtual void     SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; }
+   virtual void     SetCentMin(Float_t cent) { fCentMin = cent; }
+   virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
+   virtual void     SetNInputTracksMin(Int_t nTr) { fNInputTracksMin = nTr; }
+   virtual void     SetNInputTracksMax(Int_t nTr) { fNInputTracksMax = nTr; }
+   virtual void     SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; }
+   virtual void     SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; }
+   virtual void     SetJetPtMin(Float_t pt) { fJetPtMin = pt; }
+   virtual void     SetJetTriggerExclude(UChar_t i) { fJetTriggerExcludeMask = i; }
+   virtual void     SetJetPtFractionMin(Float_t frac) { fJetPtFractionMin = frac; }
+   virtual void     SetNMatchJets(Int_t n) { fNMatchJets = n; }
+   virtual void     SetFillEvent(Bool_t b) { fbEvent = b; }
+   virtual void     SetFillJetsMismatch1(Bool_t b) { fbJetsMismatch1 = b; }
+   virtual void     SetFillJetsMismatch2(Bool_t b) { fbJetsMismatch2 = b; }
+   virtual void     SetFillJetsRp(Bool_t b) { fbJetsRp = b; }
+   virtual void     SetFillJetsDeltaPt(Bool_t b) { fbJetsDeltaPt = b; }
+   virtual void     SetFillJetsEta(Bool_t b) { fbJetsEta = b; }
+   virtual void     SetFillJetsPhi(Bool_t b) { fbJetsPhi = b; }
+   virtual void     SetFillJetsArea(Bool_t b) { fbJetsArea = b; }
+   virtual void     SetFillJetsBeforeCut1(Bool_t b) { fbJetsBeforeCut1 = b; }
+   virtual void     SetFillJetsBeforeCut2(Bool_t b) { fbJetsBeforeCut2 = b; }
+   virtual void     SetKeepJets(Bool_t b = kTRUE) { fKeepJets = b; }
+
+private:
+   // ESD/AOD events
+   AliESDEvent *fESD;    //! ESD object
+   AliAODEvent *fAOD;    //! AOD event
+
+   // jets to compare
+   TString fJetBranchName[2]; //  name of jet branches to compare
+   TList *fListJets[2];       //! jet lists
+
+   TString fBackgroundBranch;
+
+   // event selection
+   Bool_t fIsPbPb;         // is Pb-Pb (fast embedding) or p-p (detector response)
+   AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept
+   Int_t   fMinContribVtx; // minimum number of track contributors for primary vertex
+   Float_t fVtxZMin;     // lower bound on vertex z
+   Float_t fVtxZMax;     // upper bound on vertex z
+   Int_t   fEvtClassMin;         // lower bound on event class
+   Int_t   fEvtClassMax;         // upper bound on event class
+   Float_t fCentMin;     // lower bound on centrality
+   Float_t fCentMax;     // upper bound on centrality
+   Int_t   fNInputTracksMin;  // lower bound of nb. of input tracks
+   Int_t   fNInputTracksMax;  // upper bound of nb. of input tracks
+   Float_t fJetEtaMin;     // lower bound on eta for found jets
+   Float_t fJetEtaMax;     // upper bound on eta for found jets
+   Float_t fJetPtMin;      // minimum jet pT
+   UChar_t fJetTriggerExcludeMask; // mask for jet triggeres to exclude
+   Float_t fJetPtFractionMin; // minimum fraction for positiv match of jets
+   Int_t   fNMatchJets;       // maximal nb. of jets taken for matching
+   Double_t fMatchMaxDist;     // maximal distance of matching jets
+   Bool_t  fKeepJets;          // keep jets with negative pt after background subtraction
+
+
+   // output objects
+   const Int_t fkNbranches;                   //! number of branches to be read
+   const Int_t fkEvtClasses;                  //! number of event classes
+   TList *fOutputList;                        //! output data container
+   Bool_t fbEvent;                            // fill fhnEvent
+   Bool_t fbJetsMismatch1;                    // fill fhnJetsMismatch1
+   Bool_t fbJetsMismatch2;                    // fill fhnJetsMismatch2
+   Bool_t fbJetsRp;                           // fill fhnJetsRp
+   Bool_t fbJetsDeltaPt;                      // fill fhnJetsDeltaPt
+   Bool_t fbJetsEta;                          // fill fhnJetsEta
+   Bool_t fbJetsPhi;                          // fill fhnJetsEta
+   Bool_t fbJetsArea;                         // fill fhnJetsArea
+   Bool_t fbJetsBeforeCut1;                   // fill fhnJetsBeforeCut1
+   Bool_t fbJetsBeforeCut2;                   // fill fhnJetsBeforeCut2
+   TH1I  *fHistEvtSelection;                  //! event selection statistic
+   TH1I  *fHistJetSelection;                  //! jet selection statistic
+   TH2F  *fh2JetSelection;                    //! jet selection statistic, with probe jet pt
+   THnSparse *fhnEvent;                       //! variables per event
+   THnSparse *fhnJetsMismatch1;               //! variables per jet
+   THnSparse *fhnJetsMismatch2;               //! variables per jet
+   THnSparse *fhnJetsRp;                      //! variables per jet
+   THnSparse *fhnJetsDeltaPt;                 //! variables per jet
+   THnSparse *fhnJetsEta;                     //! variables per jet
+   THnSparse *fhnJetsPhi;                     //! variables per jet
+   THnSparse *fhnJetsArea;                    //! variables per jet
+   THnSparse *fhnJetsBeforeCut1;               //! variables per jet before acceptance cut
+   THnSparse *fhnJetsBeforeCut2;               //! variables per jet before acceptance cut
+
+   AliAnalysisTaskJetResponseV2(const AliAnalysisTaskJetResponseV2&); // not implemented
+   AliAnalysisTaskJetResponseV2& operator=(const AliAnalysisTaskJetResponseV2&); // not implemented
+
+   ClassDef(AliAnalysisTaskJetResponseV2, 4);
 };
 
 #endif
+
index 557162a..fa3c4ed 100644 (file)
 
 AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(){
 
-    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-    if(!mgr){
-       ::Error("AddTaskCentralitySelection", "No analysis manager to connect ot.");
-       return NULL;
-    }
-    if(!mgr->GetInputEventHandler()){
-        ::Error("AddTaskCentralitySelection", "This task requires an input event handler.");
-       return NULL;
-    }
-
-
-    AliAnalysisTaskFastEmbedding *task = new AliAnalysisTaskFastEmbedding("FastEmbedding");
-    // ## set embedding mode ##
-    // kAODFull=0, kAODJetTracks, kAODJet4Mom, kToySingle4Mom
-    task->SetEmbedMode(AliAnalysisTaskFastEmbedding::kToyTracks);
-
-    // ## set ranges for toy ##
-    //SetToyTrackRanges(
-    Double_t minPt = 20.;   Double_t maxPt = 200.;
-    Double_t minEta = -0.5; Double_t maxEta = 0.5;
-    //Double_t minEta = -0.4; Double_t maxEta = 0.4;  // for LHC10h pass1
-    Double_t minPhi = 0.;   Double_t maxPhi = 2*TMath::Pi();
-    //fToyDistributionTrackPt: 0 = uniform distribution
-    //                         else = exponential / power law (not implemented yet)
-    //task->SetToyNumberOfTrackRange(4,4);
-    //task->SetToyTrackRanges(0.15, 300., 5,-.9, .9, 0., 2*TMath::Pi());
-    task->SetToyTrackRanges(minPt,maxPt,0.,minEta,maxEta,minPhi,maxPhi);
-    task->SetToyFilterMap((1<<32)-1);
-
-    // ## set event selection for events of the addition AOD ##
-    // kEventsAll=0; kEventsJetPt
-    task->SetEvtSelecMode(AliAnalysisTaskFastEmbedding::kEventsJetPt);
-
-    // ## set jet pT range for event selection ##
-    // SetEvtSelJetPtRange(Float_t minPt, Float_t maxPt)
-    task->SetEvtSelJetPtRange(40.,-1.);
-    //task->SetEvtSelJetEtaRange(-0.4, 0.4); // smaller eta window for LHC10h pass1
-
-    mgr->AddTask(task);
-
-    // ## create the output containers ##
-    AliAnalysisDataContainer *coutputFastEmbedding = mgr->CreateContainer(
-         "fastembedding", TList::Class(), AliAnalysisManager::kOutputContainer,
-         Form("%s:PWG4_FastEmbedding", AliAnalysisManager::GetCommonFileName()));
-
-    mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer());
-    mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
-    mgr->ConnectOutput(task, 1, coutputFastEmbedding);
-
-
-    return task;
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if(!mgr){
+      ::Error("AddTaskCentralitySelection", "No analysis manager to connect ot.");
+      return NULL;
+   }
+   if(!mgr->GetInputEventHandler()){
+      ::Error("AddTaskCentralitySelection", "This task requires an input event handler.");
+      return NULL;
+   }
+
+
+   AliAnalysisTaskFastEmbedding *task = new AliAnalysisTaskFastEmbedding("FastEmbedding");
+   // ## set embedding mode ##
+   // kAODFull=0, kAODJetTracks, kAODJet4Mom, kToySingle4Mom
+   task->SetEmbedMode(AliAnalysisTaskFastEmbedding::kToyTracks);
+
+   // ## set ranges for toy ##
+   //SetToyTrackRanges(
+   Double_t minPt = 50.;   Double_t maxPt = 100.;
+   Double_t minEta = -0.5; Double_t maxEta = 0.5;
+   //Double_t minEta = -0.4; Double_t maxEta = 0.4;  // for LHC10h pass1
+   Double_t minPhi = 0.;   Double_t maxPhi = 2*TMath::Pi();
+   //fToyDistributionTrackPt: 0 = uniform distribution
+   //                         else = exponential / power law (not implemented yet)
+   //task->SetToyNumberOfTrackRange(4,4);
+   //task->SetToyTrackRanges(0.15, 300., 5,-.9, .9, 0., 2*TMath::Pi());
+   task->SetToyTrackRanges(minPt,maxPt,0.,minEta,maxEta,minPhi,maxPhi);
+   task->SetToyFilterMap((1<<32)-1);
+
+   // ## set event selection for events of the addition AOD ##
+   // kEventsAll=0; kEventsJetPt
+   task->SetEvtSelecMode(AliAnalysisTaskFastEmbedding::kEventsJetPt);
+
+   // ## set jet pT range for event selection ##
+   // SetEvtSelJetPtRange(Float_t minPt, Float_t maxPt)
+   task->SetEvtSelJetPtRange(50.,110.);
+   //task->SetEvtSelJetEtaRange(-0.4, 0.4); // smaller eta window for LHC10h pass1
+   task->SetEvtSelJetEtaRange(-0.5, 0.5);
+
+   task->SetTrackFilterMap(272);
+
+   // event selection
+   task->SetOfflineTrgMask(AliVEvent::kMB);
+   task->SetCentMin(0.);
+   task->SetCentMax(80.);
+
+   //task->SetVtxMin(-10.);
+   //task->SetVtxMax(10.);
+
+   mgr->AddTask(task);
+
+   // ## create the output containers ##
+   AliAnalysisDataContainer *coutputFastEmbedding = mgr->CreateContainer(
+   "fastembedding", TList::Class(), AliAnalysisManager::kOutputContainer,
+   Form("%s:PWG4_FastEmbedding", AliAnalysisManager::GetCommonFileName()));
+
+   mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer());
+   mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
+   mgr->ConnectOutput(task, 1, coutputFastEmbedding);
+
+
+   return task;
 
 }
 
 AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(TObjArray* aodarray){
 
-    AliAnalysisTaskFastEmbedding *task = AddTaskFastEmbedding();
-    if(aodarray){
+   AliAnalysisTaskFastEmbedding *task = AddTaskFastEmbedding();
+   if(aodarray){
       task->SetArrayOfAODPaths(aodarray);
       task->SetEmbedMode(AliAnalysisTaskFastEmbedding::kAODFull);
-    }
+   }
 
-    return task;
+   return task;
 }
 
 
 AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(const char* filepath, Int_t mode = 0){
 
-    AliAnalysisTaskFastEmbedding *task = AddTaskFastEmbedding();
-    if(strlen(filepath)){
-
-       if(mode==0){ // path to single AOD
-          task->SetAODPath(filepath);
-       }
-       if(mode==1){ // path to text file with list of paths of multiple AODs
-           Printf("Read aod paths from file %s", filepath);
-           TObjArray* array = new TObjArray();
-           TObjString* ostr = 0;
-           TString line;
-           ifstream in;
-           in.open(filepath);
-           while(in.good()){
-              in >> line;
-              if(line.Length() == 0) continue;
-              Printf("found aod path %s", line.Data());
-              ostr = new TObjString(line.Data());
-              array->Add(ostr);
-           }
-           Printf("-> %d aod paths found", array->GetEntries());
-           
-           task->SetArrayOfAODPaths(array);
-       }
-
-       task->SetEmbedMode(AliAnalysisTaskFastEmbedding::kAODFull);
-    }
-
-    return task;
+   AliAnalysisTaskFastEmbedding *task = AddTaskFastEmbedding();
+   if(strlen(filepath)){
+
+      if(mode==0){ // path to single AOD
+         task->SetAODPath(filepath);
+      }
+      if(mode==1){ // path to text file with list of paths of multiple AODs
+         Printf("Read aod paths from file %s", filepath);
+         TObjArray* array = new TObjArray();
+         TObjString* ostr = 0;
+         TString line;
+         ifstream in;
+         in.open(filepath);
+         while(in.good()){
+            in >> line;
+            if(line.Length() == 0) continue;
+            Printf("found aod path %s", line.Data());
+            ostr = new TObjString(line.Data());
+            array->Add(ostr);
+         }
+         Printf("-> %d aod paths found", array->GetEntries());
+         
+         task->SetArrayOfAODPaths(array);
+      }
+
+      task->SetEmbedMode(AliAnalysisTaskFastEmbedding::kAODFull);
+   }
+
+   return task;
+}
+
+AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(const char* aodpath, const char* entriespath){
+
+   AliAnalysisTaskFastEmbedding *task = AddTaskFastEmbedding(aodpath, 1);
+
+   Printf("Read entries of aod files from %s", entriespath);
+   TArrayI* array = new TArrayI();
+   Int_t count = 0;
+   Int_t iEntry = -1;
+   Int_t iEntrySum = 0;
+   Int_t iEntryMax = 0;
+   TString line;
+   ifstream in;
+   in.open(entriespath);
+   while(in.good()){
+      in >> line;
+      iEntry = line.Atoi();
+
+      array->Set(count+1);
+      array->AddAt(iEntry,count);
+      count++;
+      iEntrySum += iEntry;
+      if(iEntry>iEntryMax) iEntryMax = iEntry;
+   }
+
+   task->SetArrayOfAODEntries(array);
+   task->SetAODEntriesSum(iEntrySum);
+   task->SetAODEntriesMax(iEntryMax);
+
+   return task;
 }
 
index e1605cc..5f04a22 100644 (file)
@@ -9,23 +9,23 @@ AliAnalysisTaskJetResponseV2* AddTaskJetResponseV2(Char_t* type = "clusters", Ch
 
 AliAnalysisTaskJetResponseV2* AddTaskJetResponseV2(Bool_t emb = kTRUE, Char_t* type = "clusters", Char_t* jf = "FASTKT", Float_t radius = 0.4, UInt_t filterMask = 256 , Float_t ptTrackMin = 0.15, Int_t iBack = 1, Int_t eventClassMin = 0, Int_t eventClassMax = 4){
 
-  Printf("adding task jet response\n");
-
-    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-    if(!mgr){
-       ::Error("AddTaskJetResponseV2", "No analysis manager to connect to.");
-       return NULL;
-    }
-    if(!mgr->GetInputEventHandler()){
-        ::Error("AddTaskJetResponseV2", "This task requires an input event handler.");
-       return NULL;
-    }
+   Printf("adding task jet response\n");
+
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if(!mgr){
+      ::Error("AddTaskJetResponseV2", "No analysis manager to connect to.");
+      return NULL;
+   }
+   if(!mgr->GetInputEventHandler()){
+      ::Error("AddTaskJetResponseV2", "This task requires an input event handler.");
+      return NULL;
+   }
 
    TString branch1 = "";
    TString branch2 = "";
    TString suffix  = "";
    TString suffix2 = "";
-    
+   
    if(emb){
 
       // embedding in HI event
@@ -43,12 +43,12 @@ AliAnalysisTaskJetResponseV2* AddTaskJetResponseV2(Bool_t emb = kTRUE, Char_t* t
       suffix2 += Form("_Filter%05d", filterMask);
       suffix2 += Form("_Cut%05d", (int)((1000.*ptTrackMin)));
       if(type=="clusters") suffix2 += Form("_Skip00");
-    
+      
       branch1 = Form("%sAODextraonly%s",type, suffix.Data());
       branch2 = Form("%sAODextra%s",type, suffix2.Data());
       
    } else {
-   
+      
       // p-p detector response
       suffix += Form("_%s", jf);
       suffix += Form("%02d", (int)((radius+0.01)*10.));
@@ -66,43 +66,48 @@ AliAnalysisTaskJetResponseV2* AddTaskJetResponseV2(Bool_t emb = kTRUE, Char_t* t
       
       branch1 = Form("%sAODMC2%s",type, suffix.Data()); // MC truth
       branch2 = Form("%sAOD%s",type, suffix2.Data());    // MC reconstucted
-         
+      
    }
    
    AliAnalysisTaskJetResponseV2 *task = new AliAnalysisTaskJetResponseV2(Form("JetResponseV2%s", suffix2.Data()));
-    
-    Printf("Branch1: %s",branch1.Data());
-    Printf("Branch2: %s",branch2.Data());
-    task->SetBranchNames(branch1,branch2);
-    task->SetOfflineTrgMask(AliVEvent::kMB);
-
-    task->SetEvtClassMin(eventClassMin);
-    task->SetEvtClassMax(eventClassMax);
-    task->SetCentMin(0.);
-    task->SetCentMax(100.);
    
-    task->SetJetPtMin(0.);   // min jet pt is implicit a cut on delta pT!!
+   Printf("Branch1: %s",branch1.Data());
+   Printf("Branch2: %s",branch2.Data());
+
+   task->SetBranchNames(branch1,branch2);
+   task->SetOfflineTrgMask(AliVEvent::kMB);
 
-    //task->SetNMatchJets(1); // leading jets only
+   task->SetEvtClassMin(eventClassMin);
+   task->SetEvtClassMax(eventClassMax);
+   task->SetCentMin(0.);
+   task->SetCentMax(100.);
+
+   //task->SetVtxMin(-10.);
+   //task->SetVtxMax(10.);
+   
+   task->SetJetPtMin(0.);   // min jet pt is implicit a cut on delta pT!!
 
-    if(!emb){
-       task->SetIsPbPb(kFALSE);
-       task->SetJetPtFractionMin(0.01);
-       task->SetNMatchJets(999);
-    }
+   task->SetKeepJets(kTRUE);
+
+   //task->SetNMatchJets(1); // leading jets only
+
+   if(!emb){
+      task->SetIsPbPb(kFALSE);
+      task->SetJetPtFractionMin(0.01);
+      task->SetNMatchJets(999);
+   }
 
 
-    mgr->AddTask(task);
+   mgr->AddTask(task);
 
 
-    AliAnalysisDataContainer *coutputJetResponseV2 = mgr->CreateContainer(
-         Form("jetresponseV2_%s%s", type,suffix2.Data()), TList::Class(), AliAnalysisManager::kOutputContainer,
-         Form("%s:PWG4_JetResponseV2_%s%s", AliAnalysisManager::GetCommonFileName(), type, suffix2.Data()));
+   AliAnalysisDataContainer *coutputJetResponseV2 = mgr->CreateContainer(
+   Form("jetresponseV2_%s%s", type,suffix2.Data()), TList::Class(), AliAnalysisManager::kOutputContainer,
+   Form("%s:PWG4_JetResponseV2_%s%s", AliAnalysisManager::GetCommonFileName(), type, suffix2.Data()));
 
-    mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer());
-    mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
-    mgr->ConnectOutput(task, 1, coutputJetResponseV2);
+   mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer());
+   mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
+   mgr->ConnectOutput(task, 1, coutputJetResponseV2);
 
-    return task;
+   return task;
 }