--- /dev/null
+// **************************************
+// 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.);
+ }
+
+}
+
--- /dev/null
+#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
--- /dev/null
+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;
+
+}