adding HBOM task for back extrapolation of deltapt distributions (M. Zimmermann)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Jun 2012 18:26:30 +0000 (18:26 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 12 Jun 2012 18:26:30 +0000 (18:26 +0000)
PWGJE/CMakelibPWGJE.pkg
PWGJE/PWGJELinkDef.h
PWGJE/UserTasks/AliAnalysisTaskJetHBOM.cxx [new file with mode: 0644]
PWGJE/UserTasks/AliAnalysisTaskJetHBOM.h [new file with mode: 0644]
PWGJE/macros/AddTaskJetHBOM.C [new file with mode: 0644]

index d9ecf30..74ce1c1 100644 (file)
@@ -44,7 +44,8 @@ set ( SRCS
     AliPWG4HighPtSpectra.cxx 
     AliPWG4HighPtTrackQA.cxx 
     UserTasks/AliAnalysisTaskCheckSingleTrackJetRejection.cxx 
-    UserTasks/AliAnalysisTaskJetHadronCorrelation.cxx                        
+    UserTasks/AliAnalysisTaskJetHadronCorrelation.cxx
+    UserTasks/AliAnalysisTaskJetHBOM.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index 8e5873d..bf27f78 100644 (file)
@@ -31,4 +31,5 @@
 #pragma link C++ class AliAnalysisTaskPartonDisc+;
 #pragma link C++ class AliAnalysisTaskCheckSingleTrackJetRejection+;
 #pragma link C++ class  AliAnalysisTaskJetHadronCorrelation+;
+#pragma link C++ class  AliAnalysisTaskJetHBOM+;
 #endif
diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetHBOM.cxx b/PWGJE/UserTasks/AliAnalysisTaskJetHBOM.cxx
new file mode 100644 (file)
index 0000000..9e04069
--- /dev/null
@@ -0,0 +1,958 @@
+// **************************************
+// Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
+// *******************************************
+
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+#include <TROOT.h>
+#include <TRandom3.h>
+#include <TSystem.h>
+#include <TInterpreter.h>
+#include <TChain.h>
+#include <TRefArray.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TProfile.h>
+#include <TF1.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include  "TDatabasePDG.h"
+
+#include "AliAnalysisTaskJetHBOM.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODExtension.h"
+#include "AliAODTrack.h"
+#include "AliAODJet.h"
+#include "AliAODMCParticle.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+#include "AliInputEventHandler.h"
+#include "AliAODJetEventBackground.h"
+
+
+ClassImp(AliAnalysisTaskJetHBOM)
+
+AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
+  //
+  // Destructor
+  //
+
+  delete fRef;
+  delete fRandom;
+
+  if(fTCARandomConesOut)fTCARandomConesOut->Delete();
+  delete fTCARandomConesOut;
+  
+  if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
+  delete fTCARandomConesOutRan;
+
+}
+
+AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(): 
+  AliAnalysisTaskSE(),
+  fAOD(0x0),
+  fAODExtension(0x0),
+  fRef(new TRefArray),
+  fUseAODTrackInput(kFALSE),
+  fUseAODMCInput(kFALSE),
+  fEventSelection(kFALSE),     
+  fFilterMask(0),
+  fFilterMaskBestPt(0),
+  fFilterType(0),
+  fJetTypes(kJet),
+  fTrackTypeRec(kTrackUndef),
+  fTrackTypeGen(kTrackUndef),  
+  fNSkipLeadingRan(0),
+  fNSkipLeadingCone(0),
+  fNRandomCones(0),
+  fNHBOM(0),
+  fTrackEtaWindow(0.9),    
+  //fRecEtaWindow(0.5),
+  fTrackPtCut(0.),                                                     
+  fJetOutputMinPt(0.150),
+  fMaxTrackPtInJet(100.),
+  fVtxZCut(8),
+  fVtxR2Cut(1),
+  fCentCutUp(0),
+  fCentCutLo(0),
+  fNonStdBranch(""),
+  fBackgroundBranch(""),
+  fNonStdFile(""),
+  fMomResH1(0x0),
+  fMomResH2(0x0),
+  fMomResH3(0x0),
+  fMomResH1Fit(0x0),
+  fMomResH2Fit(0x0),
+  fMomResH3Fit(0x0),
+  fhEffH1(0x0),
+  fhEffH2(0x0),
+  fhEffH3(0x0),
+  fUseTrMomentumSmearing(kFALSE),
+  fUseDiceEfficiency(kFALSE),
+  fRparam(1.0), 
+  fAlgorithm(fastjet::kt_algorithm),
+  fStrategy(fastjet::Best),
+  fRecombScheme(fastjet::BIpt_scheme),
+  fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
+  background(0),
+  fTCARandomConesOut(0x0),
+  fTCARandomConesOutRan(0x0),
+  fRandom(0),
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
+  fh1PtHard(0x0),
+  fh1PtHardNoW(0x0),  
+  fh1PtHardTrials(0x0),
+  fh1Nch(0x0),
+  fh1CentralityPhySel(0x0), 
+  fh1Centrality(0x0), 
+  fh1DeltapT(0x0),
+  fh1Rho(0x0),
+  fh1PtRandCone(0x0),
+  fh1Area(0x0),
+  fh1efficiencyPt(0x0),
+  fh2efficiencyPhi(0x0),
+  fh1ZPhySel(0x0), 
+  fh1Z(0x0), 
+  fHistList(0x0)  
+{
+  
+}
+
+AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
+  AliAnalysisTaskSE(name),
+  fAOD(0x0),
+  fAODExtension(0x0),
+  fRef(new TRefArray),
+  fUseAODTrackInput(kFALSE),
+  fUseAODMCInput(kFALSE),
+  fEventSelection(kFALSE),                                                       
+  fFilterMask(0),
+  fFilterMaskBestPt(0),
+  fFilterType(0),
+  fJetTypes(kJet),
+  fTrackTypeRec(kTrackUndef),
+  fTrackTypeGen(kTrackUndef),
+  fNSkipLeadingRan(0),
+  fNSkipLeadingCone(0),
+  fNRandomCones(0),
+  fNHBOM(0),
+  fTrackEtaWindow(0.9),    
+  //fRecEtaWindow(0.5),
+  fTrackPtCut(0.),                                                     
+  fJetOutputMinPt(0.150),
+  fMaxTrackPtInJet(100.),
+  fVtxZCut(8),
+  fVtxR2Cut(1),
+  fCentCutUp(0),
+  fCentCutLo(0),
+  fNonStdBranch(""),
+  fBackgroundBranch(""),
+  fNonStdFile(""),
+  fMomResH1(0x0),
+  fMomResH2(0x0),
+  fMomResH3(0x0),
+  fMomResH1Fit(0x0),
+  fMomResH2Fit(0x0),
+  fMomResH3Fit(0x0),
+  fhEffH1(0x0),
+  fhEffH2(0x0),
+  fhEffH3(0x0),
+  fUseTrMomentumSmearing(kFALSE),
+  fUseDiceEfficiency(kFALSE),
+  fRparam(1.0), 
+  fAlgorithm(fastjet::kt_algorithm),
+  fStrategy(fastjet::Best),
+  fRecombScheme(fastjet::BIpt_scheme),
+  fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
+  background(0),
+  fTCARandomConesOut(0x0),
+  fTCARandomConesOutRan(0x0),
+  fRandom(0),
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
+  fh1PtHard(0x0),
+  fh1PtHardNoW(0x0),  
+  fh1PtHardTrials(0x0),
+  fh1Nch(0x0),
+  fh1CentralityPhySel(0x0), 
+  fh1Centrality(0x0), 
+  fh1DeltapT(0x0),
+  fh1Rho(0x0),
+  fh1PtRandCone(0x0),
+  fh1Area(0x0),
+  fh1efficiencyPt(0x0),
+  fh2efficiencyPhi(0x0),
+  fh1ZPhySel(0x0), 
+  fh1Z(0x0), 
+  fHistList(0x0)
+{
+ DefineOutput(1, TList::Class());  
+}
+
+
+
+Bool_t AliAnalysisTaskJetHBOM::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // 
+  return kTRUE;
+}
+
+void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
+{
+
+  //
+  // Create the output container
+  //
+
+  fRandom = new TRandom3(0);
+
+
+  // Connect the AOD
+
+
+  if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
+
+  
+
+  if(fNonStdBranch.Length()!=0)
+    {
+      // only create the output branch if we have a name
+      // Create a new branch for jets...
+      //  -> cleared in the UserExec....
+      // here we can also have the case that the brnaches are written to a separate file
+      
+      // create the branch for the random cones with the same R 
+      TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
+      if(fUseDiceEfficiency || fUseTrMomentumSmearing)
+       cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
+    
+      //create array for the random cones; Until now only one cone per event is used
+      if(!AODEvent()->FindListObject(cName.Data())){
+       fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
+       fTCARandomConesOut->SetName(cName.Data());
+       AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
+      }
+      // create the branch with the random for the random cones on the random event
+      cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+      if(!AODEvent()->FindListObject(cName.Data())){
+       fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
+       fTCARandomConesOutRan->SetName(cName.Data());
+       AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
+      }
+      
+      if(fNonStdFile.Length()!=0){
+       // 
+       // case that we have an AOD extension we need to fetch the jets from the extended output
+       // we identify the extension aod event by looking for the branchname
+       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+       // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
+       fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
+      }
+    }
+
+  //  FitMomentumResolution();
+
+
+  if(!fHistList)fHistList = new TList();
+  fHistList->SetOwner();
+  PostData(1, fHistList); // post data in any case once
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  //
+  //  Histogram
+    
+  const Int_t nBinPt = 100;
+  Double_t binLimitsPt[nBinPt+1];
+  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
+    if(iPt == 0){
+      binLimitsPt[iPt] = 0.0;
+    }
+    else {// 1.0
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
+    }
+  }
+  
+  const Int_t nBinPhi = 90;
+  Double_t binLimitsPhi[nBinPhi+1];
+  for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
+    if(iPhi==0){
+      binLimitsPhi[iPhi] = -1.*TMath::Pi();
+    }
+    else{
+      binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
+    }
+  }
+
+
+
+  const Int_t nBinEta = 40;
+  Double_t binLimitsEta[nBinEta+1];
+  for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
+    if(iEta==0){
+      binLimitsEta[iEta] = -2.0;
+    }
+    else{
+      binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
+    }
+  }
+
+  const int nChMax = 5000;
+
+  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+
+  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+
+
+  fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
+  fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
+  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+
+  fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
+
+  fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
+  fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
+
+  fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
+  fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
+  
+  fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
+  fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
+  fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
+  fh1Area = new TH1F("fh1Area","area",100,0,1);
+
+  const Int_t saveLevel = 3; // large save level more histos
+  if(saveLevel>0){
+    fHistList->Add(fh1Xsec);
+    fHistList->Add(fh1Trials);
+
+    fHistList->Add(fh1Nch);
+    fHistList->Add(fh1Centrality);
+    fHistList->Add(fh1CentralityPhySel);
+    fHistList->Add(fh1Z);
+    fHistList->Add(fh1ZPhySel);
+    fHistList->Add(fh1DeltapT);
+    fHistList->Add(fh1Rho);
+    fHistList->Add(fh1PtRandCone);
+    fHistList->Add(fh1Area);
+  }
+
+  // =========== 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;
+    }
+    THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
+    if(hn)hn->Sumw2();
+  }
+  TH1::AddDirectory(oldStatus);
+}
+
+void AliAnalysisTaskJetHBOM::Init()
+{
+  //
+  // Initialization
+  //
+
+  if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
+
+  FitMomentumResolution();
+
+}
+
+void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
+{
+
+  // handle and reset the output jet branch 
+
+  if(fTCARandomConesOut)fTCARandomConesOut->Delete();
+  if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
+
+  //
+  // Execute analysis for current event
+  //
+  AliESDEvent *fESD = 0;
+  if(fUseAODTrackInput){    
+    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
+      return;
+    }
+    // fetch the header
+  }
+  else{
+    //  assume that the AOD is in the general output...
+    fAOD  = AODEvent();
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
+      return;
+    }
+    if(fDebug>0){
+      fESD = dynamic_cast<AliESDEvent*> (InputEvent());
+    }
+  }
+
+  //Check if information is provided detector level effects
+  if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
+  if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
+  
+  Bool_t selectEvent =  false;
+  Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
+
+  Float_t cent = 0;
+  Float_t zVtx  = 0;
+
+  if(fAOD){
+    const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
+    TString vtxTitle(vtxAOD->GetTitle());
+    zVtx = vtxAOD->GetZ();
+
+    cent = fAOD->GetHeader()->GetCentrality();
+    if(physicsSelection){
+      fh1CentralityPhySel->Fill(cent);
+      fh1ZPhySel->Fill(zVtx);
+    }
+    // zVertex and centrality selection
+    if(fEventSelection){
+      if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
+       Float_t yvtx = vtxAOD->GetY();
+       Float_t xvtx = vtxAOD->GetX();
+       Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
+       if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){//usual fVtxZCut=10 and fVtxR2Cut=1  // apply vertex cut later on 
+         if(physicsSelection){
+           selectEvent = true;
+         }
+       }
+      }
+      //centrality selection
+      if(fCentCutUp>0){
+       if(cent<fCentCutLo||cent>fCentCutUp){
+         selectEvent = false;
+       }
+      }
+    }else{
+      selectEvent = true;
+    }
+  }
+  
+
+  if(!selectEvent){
+    PostData(1, fHistList);
+    return;
+  }
+  fh1Centrality->Fill(cent);  
+  fh1Z->Fill(zVtx);
+  fh1Trials->Fill("#sum{ntrials}",1);
+  
+
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+  // ==== General variables needed
+
+
+
+  // we simply fetch the tracks/mc particles as a list of AliVParticles
+
+  //reconstructed particles
+  TList recParticles;
+  Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
+  Float_t nCh = recParticles.GetEntries(); 
+  fh1Nch->Fill(nCh);
+  if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
+  //nT = GetListOfTracks(&genParticles,fTrackTypeGen);
+  //if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
+
+
+  //apply efficiency fNHBOM times
+  if(fNHBOM>0){
+    for(int particle=0;particle<recParticles.GetEntries();particle++){
+      // hier Effizienzen laden und ├╝berpr├╝fen ob das Teilchen nachgewiesen wird.
+      AliVParticle *vp = (AliVParticle*)recParticles.At(particle);
+      
+      Double_t pT = vp->Pt();
+      Double_t phi = vp->Phi();
+      
+      //load efficiency
+      Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
+      Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
+      if(pT>10){
+       efficiencyPt = 0.857; //this is the result for the fit with pT>10; statistic is low for pT>10
+       //if the efficiency is not from fastMCInput_LHC10h_110719a.root this should be calculated new
+       if(fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))>0.9||fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))<0.8){
+         efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
+         Printf("%s:%d Wrong efficiency input. Check efficiency for pt>10GeV",(char*)__FILE__,__LINE__);
+       }
+      }
+      Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
+      
+      // if ran<eff -> particle is detected eff^fNHBOM = efficiency to detect particle fNHBOM times
+      Double_t ran = fRandom->Rndm();
+      if(ran>TMath::Power(eff,fNHBOM)){
+       recParticles.Remove(vp);
+      }
+    }
+  }
+
+
+  // find the jets....
+
+  vector<fastjet::PseudoJet> inputParticlesRec;
+  vector<fastjet::PseudoJet> inputParticlesRecRan;
+  
+  //randomize particles
+  AliAODJet vTmpRan(1,0,0,1);
+  for(int i = 0; i < recParticles.GetEntries(); i++){
+    AliVParticle *vp = (AliVParticle*)recParticles.At(i);
+
+    // Carefull energy is not well determined in real data, should not matter for p_T scheme?
+    // we take total momentum here
+
+      //Add particles to fastjet in case we are not running toy model
+      fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
+      jInp.set_user_index(i);
+      inputParticlesRec.push_back(jInp);
+
+    // the randomized input changes eta and phi, but keeps the p_T
+      Double_t pT = vp->Pt();
+      Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
+      Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
+      
+      Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
+      Double_t pZ = pT/TMath::Tan(theta);
+
+      Double_t pX = pT * TMath::Cos(phi);
+      Double_t pY = pT * TMath::Sin(phi);
+      Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
+      fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
+
+      jInpRan.set_user_index(i);
+      inputParticlesRecRan.push_back(jInpRan);
+      vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
+
+  }// recparticles
+
+  if(inputParticlesRec.size()==0){
+    if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
+    PostData(1, fHistList);
+    return;
+  }
+  
+  // run fast jet
+  // employ setters for these...
+
+  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
+  fastjet::AreaType areaType =   fastjet::active_area;
+  fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
+  fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
+  fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
+  
+  //range where to compute background
+  Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+  phiMin = 0;
+  phiMax = 2*TMath::Pi();
+  rapMax = fGhostEtamax - fRparam;
+  rapMin = - fGhostEtamax + fRparam;
+  fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
+
+  const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
+  const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
+
+ // loop over all jets and fill information, first one is the leading jet
+
+  if(inclusiveJets.size()>0){
+    
+    //background estimates:all bckg jets wo the 2 hardest
+    vector<fastjet::PseudoJet> jets2=sortedJets;
+    if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); //removes the two jets with the highest pT; +2 is correct ro remove 2 jets
+    Double_t bkg1=0;
+    Double_t sigma1=0.;
+    Double_t meanarea1=0.;
+    clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
+    background = bkg1;//sets background variable of the task to the correct value
+
+
+    // generate random cones
+    if(fTCARandomConesOut){       
+      // create a random jet within the acceptance
+      Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
+      Int_t nCone = 0;
+      Int_t nConeRan = 0;
+      Double_t pTC = 1; // small number
+      Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
+      Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
+      // massless jet
+      Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+      Double_t pZC = pTC/TMath::Tan(thetaC);
+      Double_t pXC = pTC * TMath::Cos(phiC);
+      Double_t pYC = pTC * TMath::Sin(phiC);
+      Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+      AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
+      
+      tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
+      if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
+      if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
+  
+      // loop over the reconstructed particles and add up the pT in the random cones
+      // maybe better to loop over randomized particles not in the real jets...
+      // but this by definition brings dow average energy in the whole  event
+      AliAODJet vTmpRanR(1,0,0,1);
+      for(int i = 0; i < recParticles.GetEntries(); i++){
+       AliVParticle *vp = (AliVParticle*)recParticles.At(i);
+       //add up energy in cone
+       if(fTCARandomConesOut){
+         AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(0);  
+         if(jC&&jC->DeltaR(vp)<fRparam){
+           if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+           jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
+         }
+       }// add up energy in cone
+       
+       // the randomized input changes eta and phi, but keeps the p_T
+       if(fTCARandomConesOutRan){
+         Double_t pTR = vp->Pt();
+         Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
+         Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
+         
+         Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
+         Double_t pZR = pTR/TMath::Tan(thetaR);
+       
+         Double_t pXR = pTR * TMath::Cos(phiR);
+         Double_t pYR = pTR * TMath::Sin(phiR);
+         Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
+         vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
+         
+         AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(0);  
+         if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
+           if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+           jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
+         }
+       }
+      }// loop over recparticles
+    }  //fTCARandomConesOut
+    Float_t jetArea = fRparam*fRparam*TMath::Pi();
+    if(fTCARandomConesOut){
+      // rescale the momentum vectors for the random cones
+      AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(0);
+      if(rC){
+       Double_t etaC = rC->Eta();
+       Double_t phiC = rC->Phi();
+       // massless jet, unit vector
+       Double_t pTC = rC->ChargedBgEnergy();
+       if(pTC<=0)pTC = 0.001; // for almost empty events
+       Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+       Double_t pZC = pTC/TMath::Tan(thetaC);
+       Double_t pXC = pTC * TMath::Cos(phiC);
+       Double_t pYC = pTC * TMath::Sin(phiC);
+       Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+       rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
+       rC->SetBgEnergy(0,0);
+       rC->SetEffArea(jetArea,0);
+      }
+    }
+    if(fTCARandomConesOutRan){
+      AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(0);
+      // same with random
+      if(rC){
+       Double_t etaC = rC->Eta();
+       Double_t phiC = rC->Phi();
+       // massless jet, unit vector
+       Double_t pTC = rC->ChargedBgEnergy();
+       if(pTC<=0)pTC = 0.001;// for almost empty events
+       Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+       Double_t pZC = pTC/TMath::Tan(thetaC);
+       Double_t pXC = pTC * TMath::Cos(phiC);
+       Double_t pYC = pTC * TMath::Sin(phiC);
+       Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+       rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
+       rC->SetBgEnergy(0,0);
+       rC->SetEffArea(jetArea,0);
+      }
+    }//find the random jets
+  }//inclusive Jets > 0
+
+ //Calculate delta pT
+ AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
+ if(randCone){
+   //background is the backbround density per area and area=pi*0.4^2 -> backgroundCone is the background energy under the cone
+   Float_t backgroundCone = background * randCone->EffectiveAreaCharged();
+   //calculates difference between expected and measured energy density
+   Float_t ptSub = randCone->Pt() - backgroundCone;
+   fh1DeltapT->Fill(ptSub);// delta pT
+   fh1Rho->Fill(background);// background rho
+   fh1PtRandCone->Fill(randCone->Pt());// pT of random cone
+   fh1Area->Fill(randCone->EffectiveAreaCharged()); // area of random cone; should always be pi*0.4^2 = 0.5
+ }else{
+   if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
+ }
+
+ if (fDebug > 2){
+   if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
+   if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
+ }
+ PostData(1, fHistList);
+}
+
+void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
+{
+  //
+  // Terminate analysis
+  //
+    if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
+
+    if(fMomResH1Fit) delete fMomResH1Fit;
+    if(fMomResH2Fit) delete fMomResH2Fit;
+    if(fMomResH3Fit) delete fMomResH3Fit;
+    
+}
+
+
+Int_t  AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
+
+  //
+  // get list of tracks/particles for different types
+  //
+
+  if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
+
+  Int_t iCount = 0;
+   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
+    if(type!=kTrackAODextraonly) {
+      AliAODEvent *aod = 0;
+      if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+      else aod = AODEvent();
+      if(!aod){
+       if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
+       return iCount;
+      }
+
+      for(int it = 0;it < aod->GetNumberOfTracks();++it){
+       AliAODTrack *tr = aod->GetTrack(it);
+       Bool_t bGood = false;
+       if(fFilterType == 0)bGood = true;
+       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
+       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
+       if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
+         if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
+         continue;
+       }
+       if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
+         if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
+         continue;
+       }
+       if(tr->Pt()<fTrackPtCut){
+         if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
+         continue;
+       }
+       if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
+       list->Add(tr);
+       iCount++;
+      }
+    }
+    if(type==kTrackAODextra || type==kTrackAODextraonly) {
+      AliAODEvent *aod = 0;
+      if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+      else aod = AODEvent();
+      
+      if(!aod){
+       return iCount;
+      }
+      TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
+      if(!aodExtraTracks)return iCount;
+      for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
+       AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
+       if (!track) continue;
+
+       AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
+       if(!trackAOD)continue;
+       Bool_t bGood = false;
+       if(fFilterType == 0)bGood = true;
+       else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
+       else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
+       if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
+        if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
+       if(trackAOD->Pt()<fTrackPtCut) continue;
+       list->Add(trackAOD);
+       iCount++;
+      }
+    }
+  }
+  else if (type ==  kTrackKineAll||type == kTrackKineCharged){
+    AliMCEvent* mcEvent = MCEvent();
+    if(!mcEvent)return iCount;
+    // we want to have alivpartilces so use get track
+    for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
+      if(!mcEvent->IsPhysicalPrimary(it))continue;
+      AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
+      if(type == kTrackKineAll){
+       if(part->Pt()<fTrackPtCut)continue;
+       list->Add(part);
+       iCount++;
+      }
+      else if(type == kTrackKineCharged){
+       if(part->Particle()->GetPDG()->Charge()==0)continue;
+       if(part->Pt()<fTrackPtCut)continue;
+       list->Add(part);
+       iCount++;
+      }
+    }
+  }
+  else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
+    AliAODEvent *aod = 0;
+    if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+    else aod = AODEvent();
+    if(!aod)return iCount;
+    TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
+    if(!tca)return iCount;
+    for(int it = 0;it < tca->GetEntriesFast();++it){
+      AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
+      if(!part->IsPhysicalPrimary())continue;
+      if(type == kTrackAODMCAll){
+       if(part->Pt()<fTrackPtCut)continue;
+       list->Add(part);
+       iCount++;
+      }
+      else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
+       if(part->Charge()==0)continue;
+       if(part->Pt()<fTrackPtCut)continue;
+       if(kTrackAODMCCharged){
+         list->Add(part);
+       }
+       else {
+         if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
+         list->Add(part);
+       }
+       iCount++;
+      }
+    }
+  }// AODMCparticle
+  list->Sort();
+  return iCount;
+}
+
+void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
+
+  //
+  // set mom res profiles
+  //
+
+  fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
+  fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
+  fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
+}
+
+void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
+  //
+  // set tracking efficiency histos
+  //
+
+  fhEffH1 = (TH1*)h1->Clone("fhEffH1");
+  fhEffH2 = (TH1*)h2->Clone("fhEffH2");
+  fhEffH3 = (TH1*)h3->Clone("fhEffH3");
+}
+
+Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
+
+  //
+  // Get smearing on generated momentum
+  //
+
+  //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
+
+  TProfile *fMomRes = 0x0;
+  if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
+  if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
+  if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
+
+  if(!fMomRes) {
+    return 0.;
+  }
+
+
+  Double_t smear = 0.;
+
+  if(pt>20.) {
+    if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
+    if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
+    if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
+  }
+  else {
+
+    Int_t bin = fMomRes->FindBin(pt);
+
+    smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
+
+  }
+
+  if(fMomRes) delete fMomRes;
+  
+  return smear;
+}
+
+void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
+  //
+  // Fit linear function on momentum resolution at high pT
+  //
+
+  if(!fMomResH1Fit && fMomResH1) {
+    fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
+    fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
+    fMomResH1Fit ->SetRange(5.,100.);
+  }
+
+  if(!fMomResH2Fit && fMomResH2) {
+    fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
+    fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
+    fMomResH2Fit ->SetRange(5.,100.);
+  }
+
+  if(!fMomResH3Fit && fMomResH3) {
+    fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
+    fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
+    fMomResH3Fit ->SetRange(5.,100.);
+  }
+
+}
+
diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetHBOM.h b/PWGJE/UserTasks/AliAnalysisTaskJetHBOM.h
new file mode 100644 (file)
index 0000000..f2eddb5
--- /dev/null
@@ -0,0 +1,228 @@
+#ifndef ALIANALYSISTASKJETHBOM_H
+#define ALIANALYSISTASKJETHBOM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// **************************************
+// Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
+// *******************************************
+
+#include  "AliAnalysisTaskSE.h"
+#include  "THnSparse.h" // cannot forward declare ThnSparseF
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.hh"
+
+
+////////////////
+class AliJetHeader;
+class AliESDEvent;
+class AliAODEvent;
+class AliAODExtension;
+class AliAODJet;
+class AliGenPythiaEventHeader;
+class AliCFManager;
+class AliAODJetEventBackground;
+class AliJetFinder;
+class TList;
+class TChain;
+class TH2F;
+class TH1F;
+class TH3F;
+class TProfile;
+class TRandom3;
+class TRefArray;
+class TClonesArray;
+class TF1;
+
+class AliAnalysisTaskJetHBOM : public AliAnalysisTaskSE
+{
+ public:
+    AliAnalysisTaskJetHBOM();
+    AliAnalysisTaskJetHBOM(const char* name);
+    virtual ~AliAnalysisTaskJetHBOM();
+    // Implementation of interface methods
+    virtual void UserCreateOutputObjects();
+    virtual void Init();
+    virtual void LocalInit() { Init(); }
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *option);
+    virtual Bool_t Notify();
+
+    
+
+    virtual void SetAODTrackInput(Bool_t b){fUseAODTrackInput = b;}
+    virtual void SetAODMCInput(Bool_t b){fUseAODMCInput = b;}
+    virtual void SetEventSelection(Bool_t b){fEventSelection = b;}
+    virtual void SetTrackEtaWindow(Float_t f){fTrackEtaWindow = f;}
+    virtual void SetTrackTypeGen(Int_t i){fTrackTypeGen = i;}
+    virtual void SetTrackTypeRec(Int_t i){fTrackTypeRec = i;}
+    virtual void SetTrackPtCut(Float_t x){fTrackPtCut = x;}
+    virtual void SetCentralityCut(Float_t xLo,Float_t xUp){fCentCutLo = xLo; fCentCutUp = xUp;}
+    virtual void SetFilterMask(UInt_t i,Int_t iType = 0){fFilterMask = i;
+      fFilterType = iType;}
+    virtual void SetJetTypes(UInt_t i){fJetTypes = i;}
+    virtual void SetVtxCuts(Float_t z,Float_t r = 1){fVtxZCut = z; fVtxR2Cut = r *r;}    
+    virtual void SetBackgroundBranch(const char* c){fBackgroundBranch = c;}
+    virtual const char* GetBackgroundBranch(){return fBackgroundBranch.Data();}    
+    virtual void SetNSkipLeadingRan(Int_t x){fNSkipLeadingRan = x;}
+    virtual void SetNSkipLeadingCone(Int_t x){fNSkipLeadingCone = x;}
+    virtual void SetNRandomCones(Int_t x){fNRandomCones = x;}
+
+    virtual void SetfNHBOM(Int_t x){fNHBOM = x;};
+    virtual void SetEfficiencyPt(TH1F *h){fh1efficiencyPt = h;}
+    virtual void SetEfficiencyPhi(TH2D *h){fh2efficiencyPhi = h;}
+
+    virtual void SetJetOutputBranch(const char *c){fNonStdBranch = c;}
+    virtual const char* GetJetOutputBranch(){return fNonStdBranch.Data();}
+    virtual void SetJetOutputFile(const char *c){fNonStdFile = c;}
+    virtual const char* GetJetOutputFile(){return fNonStdFile.Data();}
+    virtual void SetMaxTrackPtInJet(Float_t x){fMaxTrackPtInJet = x;}
+    virtual void SetJetOutputMinPt(Float_t x){fJetOutputMinPt = x;}
+
+    //Setters for detector level effects
+    virtual void SetSmearResolution(Bool_t b){fUseTrMomentumSmearing = b;} 
+    virtual void SetDiceEfficiency(Bool_t b){fUseDiceEfficiency = b;} 
+    virtual void SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3);
+    virtual void SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3);
+
+    Double_t GetMomentumSmearing(Int_t cat, Double_t pt);
+    void FitMomentumResolution();
+
+
+    // for Fast Jet
+    fastjet::JetAlgorithm        GetAlgorithm()         const {return fAlgorithm;}
+    fastjet::Strategy            GetStrategy()          const {return fStrategy;}
+    fastjet::RecombinationScheme GetRecombScheme()      const {return fRecombScheme;}
+    fastjet::AreaType            GetAreaType()          const {return fAreaType;}
+    // Setters
+    void SetRparam(Double_t f)                           {fRparam = f;}
+    void SetAlgorithm(fastjet::JetAlgorithm f)           {fAlgorithm = f;}
+    void SetStrategy(fastjet::Strategy f)                {fStrategy = f;}
+    void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;}
+    void SetAreaType(fastjet::AreaType f)                {fAreaType = f;}
+    void SetGhostArea(Double_t f) {fGhostArea = f;}
+    void SetActiveAreaRepeats(Int_t f) {fActiveAreaRepeats = f;}
+    void SetGhostEtamax(Double_t f) {fGhostEtamax = f;}
+
+
+
+    // Helper
+    //
+
+    // we have different cases
+    // AOD reading -> MC from AOD
+    // ESD reading -> MC from Kinematics
+    // this has to match with our selection of input events
+    enum {kTrackUndef = 0, kTrackAOD, kTrackKineAll,kTrackKineCharged, kTrackAODMCAll, kTrackAODMCCharged, kTrackAODMCChargedAcceptance, kTrackAODextra, kTrackAODextraonly};
+    enum {kMaxJets = 4};
+    enum {kMaxCorrelation =  3};
+    enum {kMaxRadius =       5};
+    enum {kMaxCent =         4};
+    enum {kJet = 1<<0,
+         kJetRan = 1<<1,         
+         kRC = 1<<2,
+         kRCRan = 1<<3
+    };
+    
+
+ private:
+
+    AliAnalysisTaskJetHBOM(const AliAnalysisTaskJetHBOM&);
+    AliAnalysisTaskJetHBOM& operator=(const AliAnalysisTaskJetHBOM&);
+
+    Int_t GetListOfTracks(TList *list,Int_t type);
+
+    AliAODEvent     *fAOD;                // ! where we take the jets from can be input or output AOD
+    AliAODExtension *fAODExtension;       // ! AOD extension in case we write a non-sdt branch to a separate file and the aod is standard
+    TRefArray       *fRef;               // ! trefarray for track references within the jet
+    Bool_t        fUseAODTrackInput;      // take track from input AOD not from ouptu AOD
+    Bool_t        fUseAODMCInput;         // take MC from input AOD not from ouptu AOD
+    Bool_t        fEventSelection;        // use the event selection of this task, otherwise analyse all
+    UInt_t        fFilterMask;            // filter bit for slecected tracks
+    UInt_t        fFilterMaskBestPt;      // filter bit to mark jets with high quality leading tracks
+
+    UInt_t        fFilterType;            // filter type 0 = all, 1 = ITSTPC, 2 = TPC
+    UInt_t        fJetTypes;              // 1<<0 regular jets, 1<<1 << randomized event 1<<2 random cones 1<<3 random cones randomiuzed event
+    Int_t         fTrackTypeRec;          // type of tracks used for FF 
+    Int_t         fTrackTypeGen;          // type of tracks used for FF 
+    Int_t         fNSkipLeadingRan;       // number of leading tracks to be skipped in the randomized event
+    Int_t         fNSkipLeadingCone;      // number of leading jets to be for the random cones
+    Int_t         fNRandomCones;          // number of generated random cones
+    Int_t         fNHBOM;                 // number of detector runs
+    Float_t       fTrackEtaWindow;        // eta window used for corraltion plots between rec and gen 
+    //    Float_t       fRecEtaWindow;          // eta window used for corraltion plots between rec and gen 
+    Float_t       fTrackPtCut;            // minimum track pt to be accepted
+    Float_t       fJetOutputMinPt;        // minimum p_t for jets to be written out
+    Float_t       fMaxTrackPtInJet;       // maximum track pt within a jet for flagging...
+    //    Float_t       fJetTriggerPtCut;       // minimum jwt pT for AOD to be written
+    Float_t       fVtxZCut;               // zvtx cut
+    Float_t       fVtxR2Cut;              // R vtx cut (squared) 
+    Float_t       fCentCutUp;             // upper limit on centrality
+    Float_t       fCentCutLo;             // lower limit on centrality
+    // output configurartion
+    TString       fNonStdBranch;      // the name of the non-std branch name, if empty no branch is filled
+    TString       fBackgroundBranch;  // name of the branch used for background subtraction
+    TString       fNonStdFile;        // The optional name of the output file the non-std branch is written to
+
+    //Detector level effects
+    TProfile *fMomResH1; // Momentum resolution from TrackQA Hybrid Category 1
+    TProfile *fMomResH2; // Momentum resolution from TrackQA Hybrid Category 2
+    TProfile *fMomResH3; // Momentum resolution from TrackQA Hybrid Category 3
+    TF1 *fMomResH1Fit; //fit
+    TF1 *fMomResH2Fit; //fit
+    TF1 *fMomResH3Fit; //fit
+
+    TH1      *fhEffH1;        // Efficiency for Spectra Hybrid Category 1
+    TH1      *fhEffH2;        // Efficiency for Spectra Hybrid Category 2
+    TH1      *fhEffH3;        // Efficiency for Spectra Hybrid Category 3
+    Bool_t    fUseTrMomentumSmearing;     // Apply momentum smearing on track level
+    Bool_t    fUseDiceEfficiency;         // Apply efficiency on track level by dicing
+
+    // Fast jet
+    Double_t fRparam;                  // fastjet distance parameter
+    fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm
+    fastjet::Strategy fStrategy;  //= fastjet::Best;
+    fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme;
+    fastjet::AreaType fAreaType;  // fastjet area type
+    Double_t fGhostArea;          // fasjet ghost area
+    Int_t fActiveAreaRepeats;     // fast jet active area repeats
+    Double_t fGhostEtamax;        // fast jet ghost area
+
+    Double_t background; //background rho in the event
+
+    TClonesArray  *fTCARandomConesOut;    //! TCA of output jets in randomized event
+    TClonesArray  *fTCARandomConesOutRan; //! TCA of output jets in randomized event
+    AliAODJetEventBackground *fAODJetBackgroundOut; //! jet background to be written out
+
+    TRandom3*     fRandom;   //! random number generator
+    TProfile*     fh1Xsec;   //! pythia cross section and trials
+    TH1F*         fh1Trials; //! trials are added
+    TH1F*         fh1PtHard;  //! Pt har of the event...       
+    TH1F*         fh1PtHardNoW;  //! Pt har of the event without weigt       
+    TH1F*         fh1PtHardTrials;  //! Number of trials 
+
+    TH1F*         fh1Nch;            //! charged particle mult
+    TH1F*         fh1CentralityPhySel;          // ! centrality of anaylsed events 
+    TH1F*         fh1Centrality;                // ! centrality of selected events 
+    TH1F*         fh1DeltapT;        // pT of random Cone - background energy
+    TH1F*         fh1Rho;            //background rho
+    TH1F*         fh1PtRandCone;     //pT of random Cone
+    TH1F*         fh1Area;           //area of random jet
+
+    TH1F*         fh1efficiencyPt;          //here efficiency is stored
+    TH2D*         fh2efficiencyPhi;         //here efficiency is stored
+
+    TH1F*         fh1ZPhySel;          // ! zvtx of anaylsed events 
+    TH1F*         fh1Z;                // ! zvtx of selected events 
+
+
+    TList *fHistList; //!leading tracks to be skipped in the randomized event Output list
+   
+
+    ClassDef(AliAnalysisTaskJetHBOM, 1) 
+};
+#endif
diff --git a/PWGJE/macros/AddTaskJetHBOM.C b/PWGJE/macros/AddTaskJetHBOM.C
new file mode 100644 (file)
index 0000000..de36e40
--- /dev/null
@@ -0,0 +1,172 @@
+AliAnalysisTaskJetHBOM *AddTaskJetHBOM(char* bRec = "AOD",char* bGen = "",UInt_t filterMask = 272, UInt_t iPhysicsSelectionFlag = AliVEvent::kAny,Char_t *jf = "KT", Float_t radius = 0.4,Int_t kWriteAOD = kFALSE,char* deltaFile = "",Float_t ptTrackCut = 0.15, Float_t etaTrackWindow = 0.9,Float_t vertexWindow = 10.,Int_t fNHBOM = 0);
+
+Int_t kBackgroundModeCl = 0;
+Float_t kPtTrackCutCl = 0.15;
+Float_t kTrackEtaWindowCl = 0.8;
+Float_t kVertexWindowCl = 10;
+
+
+AliAnalysisTaskJetHBOM *AddTaskJetHBOM(char* bRec,char* bGen ,UInt_t filterMask,UInt_t iPhysicsSelectionFlag,Char_t *jf,Float_t radius,Int_t kWriteAOD,char *deltaFile,Float_t ptTrackCut,Float_t etaTrackWindow,Float_t vertexWindow,Int_t fNHBOM)
+ {
+ // Creates a jet fider task, configures it and adds it to the analysis manager.
+   kPtTrackCutCl = ptTrackCut;
+   kTrackEtaWindowCl = etaTrackWindow;
+   kVertexWindowCl = vertexWindow;
+
+   TString outputFile(deltaFile);
+    // Get the pointer to the existing analysis manager via the static access method.
+    //==============================================================================
+    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+    if (!mgr) {
+       ::Error("AddTaskJetHBOM", "No analysis manager to connect to.");
+       return NULL;
+    }  
+
+    // Check the analysis type using the event handlers connected to the analysis manager.
+    //==============================================================================
+    if (!mgr->GetInputEventHandler()) {
+      ::Error("AddTaskJetHBOM", "This task requires an input event handler");
+       return NULL;
+    }
+
+    TString type = mgr->GetInputEventHandler()->GetDataType();
+    TString typeRec(bRec);
+    TString typeGen(bGen);
+    if(!typeRec.Contains("AODextra")) {
+      typeGen.ToUpper();
+      typeRec.ToUpper();
+    }
+    cout << "typeRec: " << typeRec << endl;
+    // Create the task and configure it.
+    //===========================================================================
+
+
+    TString cAdd = "";
+    cAdd += Form("%02d_",(int)((radius+0.01)*10.));
+    cAdd += Form("B%d",(int)kBackgroundModeCl);
+    cAdd += Form("_Filter%05d",filterMask);
+    cAdd += Form("_Cut%05d",(int)(1000.*kPtTrackCutCl));
+    cAdd += Form("_hbom%02d",fNHBOM);
+    Printf("%s %E",cAdd.Data(),kPtTrackCutCl);
+    AliAnalysisTaskJetHBOM* hbom = new  AliAnalysisTaskJetHBOM(Form("JetHBOM%s_%s%s",bRec,jf,cAdd.Data()));
+      
+   hbom->SetFilterMask(filterMask); 
+   //   hbom->SetUseGlobalSelection(kTRUE); 
+   hbom->SetVtxCuts(kVertexWindowCl,1);//sets fVtxZCut and fVtxR2Cut
+   if(type == "AOD"){
+     // Assume all jet are produced already
+     hbom->SetAODTrackInput(kTRUE);
+     hbom->SetAODMCInput(kTRUE);
+   }
+
+   if(typeRec.Contains("AODMC2b")){// work down from the top AODMC2b -> AODMC2 -> AODMC -> AOD
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAODMCChargedAcceptance);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(kTrackEtaWindowCl);
+   }
+   else if (typeRec.Contains("AODMC2")){
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAODMCCharged);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(5);
+   }
+   else if (typeRec.Contains("AODMC")){
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAODMCAll);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(5);
+   }
+   else if (typeRec.Contains("AODextraonly")) {
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAODextraonly);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(kTrackEtaWindowCl);
+   }
+   else if (typeRec.Contains("AODextra")) {
+     cout << "AliAnalysisTaskJetHBOM::kTrackAODextra: " << AliAnalysisTaskJetHBOM::kTrackAODextra << endl;
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAODextra);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(kTrackEtaWindowCl);
+   }
+   else if (typeRec.Contains("AOD")) {
+     hbom->SetTrackTypeRec(AliAnalysisTaskJetHBOM::kTrackAOD);
+     hbom->SetTrackPtCut(kPtTrackCutCl);
+     hbom->SetTrackEtaWindow(kTrackEtaWindowCl);
+   }
+
+   hbom->SetRparam(radius);
+   hbom->SetGhostArea(0.005);
+   hbom->SetGhostEtamax(kTrackEtaWindowCl);
+
+   switch (jf) {
+   case "ANTIKT":
+     hbom->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh
+     break;
+   case "CA":
+     hbom->SetAlgorithm(1); // CA from fastjet/JetDefinition.hh
+     break;
+   case "KT":
+     hbom->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh
+     break;
+   default:
+     ::Error("AddTaskJetHBOM", "Wrong jet finder selected\n");
+     return 0;
+   }
+
+   
+   if(kWriteAOD){
+     if(outputFile.Length())hbom->SetJetOutputFile(outputFile);
+     hbom->SetJetOutputBranch(Form("hbom%s_%s%s",bRec,jf,cAdd.Data()));
+     hbom->SetJetOutputMinPt(0); // store only jets above a certain threshold
+   }
+
+   //sets number of detector hits
+   hbom->SetfNHBOM(fNHBOM);
+
+   //physics Selection
+   //if(iPhysicsSelectionFlag)hbom->SelectCollisionCandidates(iPhysicsSelectionFlag);
+
+   mgr->AddTask(hbom);
+
+   // Create ONLY the output containers for the data produced by the task.
+   // Get and connect other common input/output containers via the manager as below
+   //==============================================================================
+   AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer(Form("hbom_%s_%s_%s%s",bRec,bGen,jf,cAdd.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,Form("%s:PWGJE_hbom_%s_%s_%s%s",AliAnalysisManager::GetCommonFileName(),bRec,bGen,jf,cAdd.Data()));
+
+   mgr->ConnectInput  (hbom, 0, mgr->GetCommonInputContainer());
+   mgr->ConnectOutput (hbom, 0, mgr->GetCommonOutputContainer());
+   mgr->ConnectOutput (hbom,  1, coutput1_Spec );
+   
+   //loads efficiencies
+   hbom->SetEfficiencyPt(GetEfficiencyPt());
+   hbom->SetEfficiencyPhi(GetEfficiencyPhi());
+
+   return hbom;
+}
+
+//loads single track pT efficiency from root file
+TH1F *GetEfficiencyPt(){
+  static TFile *fIn = 0;
+  static TH1F *hEffPt = 0; 
+
+  TString effLoc = "$ALICE_ROOT/OADB/PWGJE/HBOM/fastMCInput_LHC10h_110719a.root";
+  if(!fIn)fIn = TFile::Open(effLoc.Data());
+  if(!hEffPt)hEffPt = (TH1F*)fIn->Get("hSingleTrackEffPt");
+  if(!hEffPt) cout<<"Could not load hSingleTrackEffPt"<<endl; 
+
+  if(hEffPt){return hEffPt;}
+  return 0;
+
+}
+
+//loads phi-pT efficiency from root file
+TH2D *GetEfficiencyPhi(){
+  static TFile *fIn = 0;
+  static TH2D *hPhiPt = 0; 
+
+  TString effLoc = "$ALICE_ROOT/OADB/PWGJE/HBOM/fastMCInput_LHC10h_110719a.root";
+  if(!fIn)fIn = TFile::Open(effLoc.Data());
+  if(!hPhiPt)hPhiPt = (TH2D*)fIn->Get("h2TrackPtPhiNorm");
+  if(!hPhiPt) cout<<"Could not load h2TrackPtPhiNorm"<<endl; 
+
+  if(hPhiPt){return hPhiPt;}
+  return 0;
+
+}