X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FAliAnalysisTaskJetResponseV2.cxx;h=b7a016eef42a14db42c7b3c4e53f69e315dabdf9;hb=291f2e16161d0a059222a86dc3a09a10cb327dc4;hp=6b4891f078c6173f88427ab71b686f7868414d1b;hpb=8c483a6b3d4c70a2dbcc2ad00e635c36f9a5001f;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/AliAnalysisTaskJetResponseV2.cxx b/PWGJE/AliAnalysisTaskJetResponseV2.cxx index 6b4891f078c..b7a016eef42 100644 --- a/PWGJE/AliAnalysisTaskJetResponseV2.cxx +++ b/PWGJE/AliAnalysisTaskJetResponseV2.cxx @@ -14,7 +14,7 @@ **************************************************************************/ // -// task compares jets in two branches, +// task compares jets in two or three branches, // written for analysis of jet embedding in HI events // @@ -45,1077 +45,1176 @@ #include "AliAODEvent.h" #include "AliAODJet.h" +#include "AliAODHandler.h" #include "AliAnalysisTaskJetResponseV2.h" +using std::cout; +using std::endl; + ClassImp(AliAnalysisTaskJetResponseV2) AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() : AliAnalysisTaskSE(), -fESD(0x0), -fAOD(0x0), -fBackgroundBranch(""), -fIsPbPb(kTRUE), -fOfflineTrgMask(AliVEvent::kAny), -fMinContribVtx(1), -fVtxZMin(-8.), -fVtxZMax(8.), -fEvtClassMin(0), -fEvtClassMax(4), -fCentMin(0.), -fCentMax(100.), -fNInputTracksMin(0), -fNInputTracksMax(-1), -fJetEtaMin(-.5), -fJetEtaMax(.5), -fJetPtMin(20.), -fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), -fJetPtFractionMin(0.5), -fNMatchJets(4), -fMatchMaxDist(0.8), -fKeepJets(kFALSE), -fkNbranches(2), -fkEvtClasses(12), -fOutputList(0x0), -fbEvent(kTRUE), -fbJetsMismatch1(kTRUE), -fbJetsMismatch2(kTRUE), -fbJetsRp(kTRUE), -fbJetsDeltaPt(kTRUE), -fbJetsEta(kTRUE), -fbJetsPhi(kTRUE), -fbJetsArea(kTRUE), -fbJetsBeforeCut1(kTRUE), -fbJetsBeforeCut2(kTRUE), -fHistEvtSelection(0x0), -fHistJetSelection(0x0), -fh2JetSelection(0x0), -fhnEvent(0x0), -fhnJetsMismatch1(0x0), -fhnJetsMismatch2(0x0), -fhnJetsRp(0x0), -fhnJetsDeltaPt(0x0), -fhnJetsEta(0x0), -fhnJetsPhi(0x0), -fhnJetsArea(0x0), -fhnJetsBeforeCut1(0x0), -fhnJetsBeforeCut2(0x0) + fESD(0x0), + fAOD(0x0), + fAODOut(0x0), + fAODExtension(0x0), + fNonStdFile(""), + fBackgroundBranch(""), + fIsPbPb(kTRUE), + fOfflineTrgMask(AliVEvent::kAny), + fMinContribVtx(1), + fVtxZMin(-8.), + fVtxZMax(8.), + fEvtClassMin(0), + fEvtClassMax(4), + fCentMin(0.), + fCentMax(100.), + fNInputTracksMin(0), + fNInputTracksMax(-1), + fJetEtaMin(-.5), + fJetEtaMax(.5), + fJetPtMin(20.), + fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), + fJetPtFractionMin(0.5), + fNMatchJets(4), + fMatchMaxDist(0.8), + fKeepJets(kFALSE), + fkNbranches(2), + fkEvtClasses(12), + fOutputList(0x0), + fbEvent(kTRUE), + fbJetsMismatch1(kTRUE), + fbJetsMismatch2(kTRUE), + fbJetsRp(kTRUE), + fbJetsDeltaPt(kTRUE), + fbJetsEta(kTRUE), + fbJetsPhi(kTRUE), + fbJetsArea(kTRUE), + fbJets3Branches(kFALSE), + fbJetsBeforeCut1(kTRUE), + fbJetsBeforeCut2(kTRUE), + fHistEvtSelection(0x0), + fHistJetSelection(0x0), + fh2JetSelection(0x0), + fhnEvent(0x0), + fhnJetsMismatch1(0x0), + fhnJetsMismatch2(0x0), + fhnJetsRp(0x0), + fhnJetsDeltaPt(0x0), + fhnJetsEta(0x0), + fhnJetsPhi(0x0), + fhnJetsArea(0x0), + fhnJets3Branches(0x0), + fhnJetsBeforeCut1(0x0), + fhnJetsBeforeCut2(0x0) { - // default Constructor + // default Constructor - fJetBranchName[0] = ""; - fJetBranchName[1] = ""; + fJetBranchName[0] = ""; + fJetBranchName[1] = ""; + fJetBranchName[2] = ""; - fListJets[0] = new TList; - fListJets[1] = new TList; + fListJets[0] = new TList; + fListJets[1] = new TList; + fListJets[2] = new TList; } AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) : -AliAnalysisTaskSE(name), -fESD(0x0), -fAOD(0x0), -fBackgroundBranch(""), -fIsPbPb(kTRUE), -fOfflineTrgMask(AliVEvent::kAny), -fMinContribVtx(1), -fVtxZMin(-8.), -fVtxZMax(8.), -fEvtClassMin(0), -fEvtClassMax(4), -fCentMin(0.), -fCentMax(100.), -fNInputTracksMin(0), -fNInputTracksMax(-1), -fJetEtaMin(-.5), -fJetEtaMax(.5), -fJetPtMin(20.), -fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), -fJetPtFractionMin(0.5), -fNMatchJets(4), -fMatchMaxDist(0.8), -fKeepJets(kFALSE), -fkNbranches(2), -fkEvtClasses(12), -fOutputList(0x0), -fbEvent(kTRUE), -fbJetsMismatch1(kTRUE), -fbJetsMismatch2(kTRUE), -fbJetsRp(kTRUE), -fbJetsDeltaPt(kTRUE), -fbJetsEta(kTRUE), -fbJetsPhi(kTRUE), -fbJetsArea(kTRUE), -fbJetsBeforeCut1(kTRUE), -fbJetsBeforeCut2(kTRUE), -fHistEvtSelection(0x0), -fHistJetSelection(0x0), -fh2JetSelection(0x0), -fhnEvent(0x0), -fhnJetsMismatch1(0x0), -fhnJetsMismatch2(0x0), -fhnJetsRp(0x0), -fhnJetsDeltaPt(0x0), -fhnJetsEta(0x0), -fhnJetsPhi(0x0), -fhnJetsArea(0x0), -fhnJetsBeforeCut1(0x0), -fhnJetsBeforeCut2(0x0) + AliAnalysisTaskSE(name), + fESD(0x0), + fAOD(0x0), + fAODOut(0x0), + fAODExtension(0x0), + fNonStdFile(""), + fBackgroundBranch(""), + fIsPbPb(kTRUE), + fOfflineTrgMask(AliVEvent::kAny), + fMinContribVtx(1), + fVtxZMin(-8.), + fVtxZMax(8.), + fEvtClassMin(0), + fEvtClassMax(4), + fCentMin(0.), + fCentMax(100.), + fNInputTracksMin(0), + fNInputTracksMax(-1), + fJetEtaMin(-.5), + fJetEtaMax(.5), + fJetPtMin(20.), + fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), + fJetPtFractionMin(0.5), + fNMatchJets(4), + fMatchMaxDist(0.8), + fKeepJets(kFALSE), + fkNbranches(2), + fkEvtClasses(12), + fOutputList(0x0), + fbEvent(kTRUE), + fbJetsMismatch1(kTRUE), + fbJetsMismatch2(kTRUE), + fbJetsRp(kTRUE), + fbJetsDeltaPt(kTRUE), + fbJetsEta(kTRUE), + fbJetsPhi(kTRUE), + fbJetsArea(kTRUE), + fbJets3Branches(kFALSE), + fbJetsBeforeCut1(kTRUE), + fbJetsBeforeCut2(kTRUE), + fHistEvtSelection(0x0), + fHistJetSelection(0x0), + fh2JetSelection(0x0), + fhnEvent(0x0), + fhnJetsMismatch1(0x0), + fhnJetsMismatch2(0x0), + fhnJetsRp(0x0), + fhnJetsDeltaPt(0x0), + fhnJetsEta(0x0), + fhnJetsPhi(0x0), + fhnJetsArea(0x0), + fhnJets3Branches(0x0), + fhnJetsBeforeCut1(0x0), + fhnJetsBeforeCut2(0x0) { - // Constructor + // Constructor - fJetBranchName[0] = ""; - fJetBranchName[1] = ""; + fJetBranchName[0] = ""; + fJetBranchName[1] = ""; + fJetBranchName[2] = ""; - fListJets[0] = new TList; - fListJets[1] = new TList; + fListJets[0] = new TList; + fListJets[1] = new TList; + fListJets[2] = new TList; - DefineOutput(1, TList::Class()); + DefineOutput(1, TList::Class()); } AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2() { - delete fListJets[0]; - delete fListJets[1]; + delete fListJets[0]; + delete fListJets[1]; + delete fListJets[2]; } -void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2) +void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3) { - fJetBranchName[0] = branch1; - fJetBranchName[1] = branch2; + fJetBranchName[0] = branch1; + fJetBranchName[1] = branch2; + fJetBranchName[2] = branch3; + + if(strlen(fJetBranchName[2].Data()) ) { + fkNbranches = 3; + fbJets3Branches = kTRUE; + } } void AliAnalysisTaskJetResponseV2::Init() { - // check for jet branches - if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){ - AliError("Jet branch name not set."); - } + // check for jet branches + if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){ + AliError("Jet branch name not set."); + } } void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects() { - // Create histograms - // Called once - OpenFile(1); - if(!fOutputList) fOutputList = new TList; - fOutputList->SetOwner(kTRUE); - - Bool_t oldStatus = TH1::AddDirectoryStatus(); - TH1::AddDirectory(kFALSE); - - - fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5); - fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); - fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN"); - fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)"); - fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)"); - fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)"); - fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)"); - - fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5); - fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); - fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN"); - fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet"); - fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list"); - fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut"); - fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut"); - fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut"); - fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask"); - - fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.); - fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); - fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN"); - fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet"); - fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list"); - fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut"); - fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut"); - fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut"); - fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask"); - - - UInt_t entries = 0; // bit coded, see GetDimParams() below - UInt_t opt = 0; // bit coded, default (0) or high resolution (1) - - if(fbEvent){ - entries = 1<<0 | 1<<1 | 1<<2 | 1<<26; // cent : nInpTrks : rp psi : pT hard bin - opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution - fhnEvent = NewTHnSparseF("fhnEvent", entries, opt); - } + // Create histograms + // Called once + OpenFile(1); + if(!fOutputList) fOutputList = new TList; + fOutputList->SetOwner(kTRUE); + + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + + fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5); + fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); + fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN"); + fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)"); + fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)"); + fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)"); + fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)"); + + fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5); + fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); + fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN"); + fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet"); + fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list"); + fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut"); + fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut"); + fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut"); + fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask"); + + fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.); + fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); + fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN"); + fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet"); + fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list"); + fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut"); + fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut"); + fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut"); + fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask"); + + + UInt_t entries = 0; // bit coded, see GetDimParams() below + UInt_t opt = 0; // bit coded, default (0) or high resolution (1) + + if(fbEvent){ + entries = 1<<0 | 1<<1 | 1<<2 | 1<<26; // cent : nInpTrks : rp psi : pT hard bin + opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution + fhnEvent = NewTHnSparseF("fhnEvent", entries, opt); + } - if(fbJetsMismatch1){ // real mismatch (no related rec jet found) - // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin - entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26; - opt = 1<<6 | 1<<8 | 1<<10; - fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt); - } - - if(fbJetsMismatch2){ // acceptance + fraction cut - // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin - entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26; - opt = 1<<6 | 1<<7 | 1<<8 | 1<<9; - fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt); - - } + if(fbJetsMismatch1){ // real mismatch (no related rec jet found) + // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin + entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26; + opt = 1<<6 | 1<<8 | 1<<10; + fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt); + } + + if(fbJetsMismatch2){ // acceptance + fraction cut + // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin + entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26; + opt = 1<<6 | 1<<7 | 1<<8 | 1<<9; + fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt); + + } - if(fbJetsRp){ - // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho - /* - entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25; - opt = 1<<4 | 1<<5;*/ - // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin - entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26; - opt = 1<<4 | 1<<5; - fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt); - } - - // cent : nInpTrks : rp bins: deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt) - if(fbJetsDeltaPt){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26; - opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18; - fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt); - } - - // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta) - if(fbJetsEta){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26; - opt = 1<<15 | 1<<8 | 1<<9; - fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt); - } - - // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin - if(fbJetsPhi){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26; - opt = 1<<10 | 1<<11; - fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt); - } - - // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area) - if(fbJetsArea){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26; - opt = 1<<18 | 1<<12 | 1<<13; - fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt); - } - - - //before cut - - // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible) - if(fbJetsBeforeCut1){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11; - opt = 0; - fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt); - } - - // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution) - if(fbJetsBeforeCut2){ - entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9; - opt = 0; - fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt); - } - - fOutputList->Add(fHistEvtSelection); - fOutputList->Add(fHistJetSelection); - fOutputList->Add(fh2JetSelection); - if(fbEvent) fOutputList->Add(fhnEvent); - if(fbJetsMismatch1) fOutputList->Add(fhnJetsMismatch1); - if(fbJetsMismatch2) fOutputList->Add(fhnJetsMismatch2); - if(fbJetsRp) fOutputList->Add(fhnJetsRp); - if(fbJetsDeltaPt) fOutputList->Add(fhnJetsDeltaPt); - if(fbJetsEta) fOutputList->Add(fhnJetsEta); - if(fbJetsPhi) fOutputList->Add(fhnJetsPhi); - if(fbJetsArea) fOutputList->Add(fhnJetsArea); - if(fbJetsBeforeCut1) fOutputList->Add(fhnJetsBeforeCut1); - if(fbJetsBeforeCut2) fOutputList->Add(fhnJetsBeforeCut2); - - // =========== Switch on Sumw2 for all histos =========== - for (Int_t i=0; iGetEntries(); ++i) { - TH1 *h1 = dynamic_cast(fOutputList->At(i)); - if (h1){ - h1->Sumw2(); - continue; - } - THnSparse *hn = dynamic_cast(fOutputList->At(i)); - if (hn){ - hn->Sumw2(); - } - } - TH1::AddDirectory(oldStatus); - - PostData(1, fOutputList); + if(fbJetsRp){ + // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho + /* + entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25; + opt = 1<<4 | 1<<5;*/ + // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin + entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26; + opt = 1<<4 | 1<<5; + fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt); + } + + // cent : nInpTrks : rp bins: deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt) + if(fbJetsDeltaPt){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26; + opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18; + fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt); + } + + // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta) + if(fbJetsEta){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26; + opt = 1<<15 | 1<<8 | 1<<9; + fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt); + } + + // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin + if(fbJetsPhi){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26; + opt = 1<<10 | 1<<11; + fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt); + } + + // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area) + if(fbJetsArea){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26; + opt = 1<<18 | 1<<12 | 1<<13; + fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt); + } + + // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin + if(fbJets3Branches){ + entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<17 | 1<<26; + opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28; + fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt); + } + + + + + //before cut + + // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible) + if(fbJetsBeforeCut1){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11; + opt = 0; + fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt); + } + + // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution) + if(fbJetsBeforeCut2){ + entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9; + opt = 0; + fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt); + } + + fOutputList->Add(fHistEvtSelection); + fOutputList->Add(fHistJetSelection); + fOutputList->Add(fh2JetSelection); + if(fbEvent) fOutputList->Add(fhnEvent); + if(fbJetsMismatch1) fOutputList->Add(fhnJetsMismatch1); + if(fbJetsMismatch2) fOutputList->Add(fhnJetsMismatch2); + if(fbJetsRp) fOutputList->Add(fhnJetsRp); + if(fbJetsDeltaPt) fOutputList->Add(fhnJetsDeltaPt); + if(fbJetsEta) fOutputList->Add(fhnJetsEta); + if(fbJetsPhi) fOutputList->Add(fhnJetsPhi); + if(fbJetsArea) fOutputList->Add(fhnJetsArea); + if(fbJets3Branches) fOutputList->Add(fhnJets3Branches); + if(fbJetsBeforeCut1) fOutputList->Add(fhnJetsBeforeCut1); + if(fbJetsBeforeCut2) fOutputList->Add(fhnJetsBeforeCut2); + + // =========== Switch on Sumw2 for all histos =========== + for (Int_t i=0; iGetEntries(); ++i) { + TH1 *h1 = dynamic_cast(fOutputList->At(i)); + if (h1){ + h1->Sumw2(); + continue; + } + THnSparse *hn = dynamic_cast(fOutputList->At(i)); + if (hn){ + hn->Sumw2(); + } + } + TH1::AddDirectory(oldStatus); + + PostData(1, fOutputList); } void AliAnalysisTaskJetResponseV2::UserExec(Option_t *) { - // load events, apply event cuts, then compare leading jets - - if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){ - AliError("Jet branch name not set."); - return; - } - - fESD=dynamic_cast(InputEvent()); - if (!fESD) { - AliError("ESD not available"); - fAOD = dynamic_cast(InputEvent()); - } else { - fAOD = dynamic_cast(AODEvent()); - } - if (!fAOD) { - AliError("AOD not available"); - return; - } - - // -- event selection -- - fHistEvtSelection->Fill(1); // number of events before event selection - - // physics selection - AliInputEventHandler* inputHandler = (AliInputEventHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); - if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){ - if(fDebug) Printf(" Trigger Selection: event REJECTED ... "); - fHistEvtSelection->Fill(2); - PostData(1, fOutputList); - return; - } - - // vertex selection - AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); - Int_t nTracksPrim = primVtx->GetNContributors(); - if ((nTracksPrim < fMinContribVtx) || - (primVtx->GetZ() < fVtxZMin) || - (primVtx->GetZ() > fVtxZMax) ){ - if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); - fHistEvtSelection->Fill(3); - PostData(1, fOutputList); - return; - } - - // event class selection (from jet helper task) - Int_t eventClass = AliAnalysisHelperJetTasks::EventClass(); - if(fDebug) Printf("Event class %d", eventClass); - if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){ - fHistEvtSelection->Fill(4); - PostData(1, fOutputList); - return; - } - - // centrality selection - AliCentrality *cent = 0x0; - Float_t centValue = 0.; - if(fESD) cent = fESD->GetCentrality(); - if(cent) centValue = cent->GetCentralityPercentile("V0M"); - if(fDebug) printf("centrality: %f\n", centValue); - if (centValue < fCentMin || centValue > fCentMax){ - fHistEvtSelection->Fill(4); - PostData(1, fOutputList); - return; - } - - - // multiplicity due to input tracks - Int_t nInputTracks = GetNInputTracks(); - - if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){ - fHistEvtSelection->Fill(5); - PostData(1, fOutputList); - return; - } + // load events, apply event cuts, then compare leading jets + + if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){ + AliError("Jet branch name not set."); + return; + } + + fESD=dynamic_cast(InputEvent()); + if (!fESD) { + AliError("ESD not available"); + + fAOD = dynamic_cast(InputEvent()); + // assume that the AOD is in the general output... + fAODOut = AODEvent(); + } else { + fAOD = dynamic_cast(AODEvent()); + } + if (!fAOD) { + AliError("AOD not available"); + return; + } + + if(fNonStdFile.Length()!=0){ + // case that we have an AOD extension we need can fetch the jets from the extended output + AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); + fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0); + if(!fAODExtension){ + if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data()); + } + } + + // -- event selection -- + fHistEvtSelection->Fill(1); // number of events before event selection + + // physics selection + AliInputEventHandler* inputHandler = (AliInputEventHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); + if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){ + if(fDebug) Printf(" Trigger Selection: event REJECTED ... "); + fHistEvtSelection->Fill(2); + PostData(1, fOutputList); + return; + } + + // vertex selection + AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); + Int_t nTracksPrim = primVtx->GetNContributors(); + if ((nTracksPrim < fMinContribVtx) || + (primVtx->GetZ() < fVtxZMin) || + (primVtx->GetZ() > fVtxZMax) ){ + if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); + fHistEvtSelection->Fill(3); + PostData(1, fOutputList); + return; + } + + // event class selection (from jet helper task) + Int_t eventClass = AliAnalysisHelperJetTasks::EventClass(); + if(fDebug) Printf("Event class %d", eventClass); + if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){ + fHistEvtSelection->Fill(4); + PostData(1, fOutputList); + return; + } + + // centrality selection + AliCentrality *cent = 0x0; + Float_t centValue = 0.; + if(fESD) cent = fESD->GetCentrality(); + if(cent) centValue = cent->GetCentralityPercentile("V0M"); + if(fDebug) printf("centrality: %f\n", centValue); + if (centValue < fCentMin || centValue > fCentMax){ + fHistEvtSelection->Fill(4); + PostData(1, fOutputList); + return; + } + + + // multiplicity due to input tracks + Int_t nInputTracks = GetNInputTracks(); + + if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){ + fHistEvtSelection->Fill(5); + PostData(1, fOutputList); + return; + } - fHistEvtSelection->Fill(0); // accepted events - // -- end event selection -- - - // pt hard - Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard(); - Int_t pthardbin = GetPtHardBin(pthard); - - // reaction plane - Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE); - if(fbEvent){ - Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp }; - fhnEvent->Fill(eventEntries); - } - - - // get background - AliAODJetEventBackground* externalBackground = 0; - if(!externalBackground&&fBackgroundBranch.Length()){ - externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data())); - //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; - } - Float_t rho = 0; - if(externalBackground)rho = externalBackground->GetBackground(0); - - - // fetch jets - TClonesArray *aodJets[2]; - aodJets[0] = dynamic_cast(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet - aodJets[1] = dynamic_cast(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE - - - - for (Int_t iJetType = 0; iJetType < 2; iJetType++) { - fListJets[iJetType]->Clear(); - if (!aodJets[iJetType]) continue; - - if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast()); + fHistEvtSelection->Fill(0); // accepted events + // -- end event selection -- + + // pt hard + Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard(); + Int_t pthardbin = GetPtHardBin(pthard); + + // reaction plane + Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE); + if(fbEvent){ + Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp }; + fhnEvent->Fill(eventEntries); + } + + + // get background + AliAODJetEventBackground* externalBackground = 0; + if(!externalBackground&&fBackgroundBranch.Length()){ + externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data())); + if(!externalBackground && fAODOut) externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data())); + if(!externalBackground && fAODExtension) externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); + //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; + } + Float_t rho = 0; + if(externalBackground) rho = externalBackground->GetBackground(0); + + + // fetch jets + TClonesArray *aodJets[3]; + aodJets[0] = dynamic_cast(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet + if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet + if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); + aodJets[1] = dynamic_cast(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1 + if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); + if( strlen(fJetBranchName[2].Data()) ) { + aodJets[2] = dynamic_cast(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2 + if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data())); + if(aodJets[2]) fkNbranches=3; + if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data()); + } + + if(fDebug>10) printf("fkNbranches %d\n",fkNbranches); + for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) { + fListJets[iJetType]->Clear(); + if (!aodJets[iJetType]) { + if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data()); + continue; + } + if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast()); - for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) { - AliAODJet *jet = dynamic_cast((*aodJets[iJetType])[iJet]); - if (jet) fListJets[iJetType]->Add(jet); - } - fListJets[iJetType]->Sort(); - } + for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) { + AliAODJet *jet = dynamic_cast((*aodJets[iJetType])[iJet]); + if (jet) fListJets[iJetType]->Add(jet); + } + fListJets[iJetType]->Sort(); + } - // jet matching - static TArrayI aMatchIndex(fListJets[0]->GetEntries()); - static TArrayF aPtFraction(fListJets[0]->GetEntries()); - if(aMatchIndex.GetSize()GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries()); - if(aPtFraction.GetSize()GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries()); - - // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction' - AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), - fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()), - aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2); + // jet matching + static TArrayI aMatchIndex(fListJets[0]->GetEntries()); + static TArrayF aPtFraction(fListJets[0]->GetEntries()); + if(aMatchIndex.GetSize()GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries()); + if(aPtFraction.GetSize()GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries()); + + // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction' + AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), + fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()), + aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2); + static TArrayI aMatchIndexv2(fListJets[0]->GetEntries()); + static TArrayF aPtFractionv2(fListJets[0]->GetEntries()); + if(aMatchIndexv2.GetSize()GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries()); + if(aPtFractionv2.GetSize()GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries()); + if( strlen(fJetBranchName[2].Data()) ) { + AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), + fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()), + aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2); + } - // loop over matched jets - Int_t ir = -1; // index of matched reconstruced jet - Float_t fraction = -1.; - AliAODJet *jet[2] = { 0x0, 0x0 }; - Float_t jetEta[2] = { -990., -990. }; - Float_t jetPhi[2] = { -990., -990. }; - Float_t jetPt[2] = { -990., -990. }; - Float_t jetArea[2] = { -990., -990. }; - Float_t rpJet[2] = { -990., -990. }; - Int_t rpBin = -990; + // loop over matched jets + Int_t ir = -1; // index of matched reconstruced jet + Int_t ir2 = -1; // index of matched reconstruced jet + Float_t fraction = -1.; + Float_t fraction2 = -1.; + AliAODJet *jet[3] = { 0x0, 0x0, 0x0 }; + Float_t jetEta[3] = { -990., -990., -990. }; + Float_t jetPhi[3] = { -990., -990., -990. }; + Float_t jetPt[3] = { -990., -990., -990. }; + Float_t jetArea[3] = { -990., -990., -990. }; + Float_t rpJet[3] = { -990., -990., -990. }; + Int_t rpBin = -990; - for(Int_t ig=0; igGetEntries(); ++ig){ - ir = aMatchIndex[ig]; - - //fetch jets - jet[0] = (AliAODJet*)(fListJets[0]->At(ig)); - if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir)); - else jet[1] = 0x0; - - for(Int_t i=0; iEta(); - jetPhi[i] = jet[i]->Phi(); - jetPt[i] = GetPt(jet[i], i); - jetArea[i] = jet[i]->EffectiveAreaCharged(); - rpJet[i] = TVector2::Phi_mpi_pi(rp-jetPhi[i]); - if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3); - } + for(Int_t ig=0; igGetEntries(); ++ig){ + ir = aMatchIndex[ig]; + if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig]; + else ir2 =0; + + //fetch jets + jet[0] = (AliAODJet*)(fListJets[0]->At(ig)); + if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir)); + else jet[1] = 0x0; + if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2)); + else jet[2] = 0x0; + + for(Int_t i=0; iEta(); + jetPhi[i] = jet[i]->Phi(); + jetPt[i] = GetPt(jet[i], i); + jetArea[i] = jet[i]->EffectiveAreaCharged(); + rpJet[i] = TVector2::Phi_mpi_pi(rp-jetPhi[i]); + if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3); } - fraction = aPtFraction[ig]; - - // jet statistics - fHistJetSelection->Fill(1); // all probe jets - if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets - - if(ir<0){ - fHistJetSelection->Fill(2); - if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]); + } + fraction = aPtFraction[ig]; + if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig]; + else fraction2 = 0.; + + // jet statistics + fHistJetSelection->Fill(1); // all probe jets + if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets + + if(ir<0){ + fHistJetSelection->Fill(2); + if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]); - if(fbJetsMismatch1){ - if(!jet[0]) continue; - Double_t jetEntriesMismatch1[7] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0] - }; - fhnJetsMismatch1->Fill(jetEntriesMismatch1); - } - - continue; + if(fbJetsMismatch1){ + if(!jet[0]) continue; + Double_t jetEntriesMismatch1[7] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0] + }; + fhnJetsMismatch1->Fill(jetEntriesMismatch1); } + + continue; + } - if(!jet[0] || !jet[1]){ - fHistJetSelection->Fill(3); - if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]); - continue; - } + if(!jet[0] || !jet[1]){ + fHistJetSelection->Fill(3); + if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]); + continue; + } - // look for distance to next rec jet - Float_t distNextJet = -0.01; // no neighbor - Float_t ptNextJet = -1.; - for(Int_t r=0; rGetEntries(); ++r){ - if(r==ir) continue; - Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r)); - if(distNextJet<0. || distNextJet>tmpDeltaR){ - distNextJet = tmpDeltaR; - if(fKeepJets) ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0); - else ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->Pt(); - } + // look for distance to next rec jet + Float_t distNextJet = -0.01; // no neighbor + Float_t ptNextJet = -1.; + for(Int_t r=0; rGetEntries(); ++r){ + if(r==ir) continue; + Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r)); + if(distNextJet<0. || distNextJet>tmpDeltaR){ + distNextJet = tmpDeltaR; + if(fKeepJets) ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0); + else ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->Pt(); } + } - //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.; - //Float_t relRho = rho>0. ? localRho / rho : 0.; + //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.; + //Float_t relRho = rho>0. ? localRho / rho : 0.; - // calculate parameters of associated jets - /* from Leticia, kT clusterizer - Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 }; - Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 }; - Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 }; - */ - // own, from embedded tracks - Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 }; - Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 }; - Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 }; + // calculate parameters of associated jets + /* from Leticia, kT clusterizer + Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 }; + Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 }; + Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 }; + */ + // own, from embedded tracks + Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 }; + Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 }; + Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 }; - Float_t rpCorr = 0.; + Float_t rpCorr = 0.; - if(eventClass>0&&eventClass<4){ - rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]); - rpCorr *= rho * jetArea[1]; - } - - Float_t deltaPt = jetPt[1]-jetPt[0]; - Float_t deltaEta = jetEta[1]-jetEta[0]; - Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]); - Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi); - Float_t deltaArea = jetArea[1]-jetArea[0]; + if(eventClass>0&&eventClass<4){ + rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]); + rpCorr *= rho * jetArea[1]; + } + + Float_t delta = jetPt[2]-jetPt[1]; + Float_t deltaPt = jetPt[1]-jetPt[0]; + Float_t deltaEta = jetEta[1]-jetEta[0]; + Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]); + Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi); + Float_t deltaArea = jetArea[1]-jetArea[0]; - // fill thnsparse before acceptance cut - if(fbJetsBeforeCut1){ - Double_t jetBeforeCutEntries1[10] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], - (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction - }; - fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1); - } + // fill thnsparse before acceptance cut + if(fbJetsBeforeCut1){ + Double_t jetBeforeCutEntries1[10] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], + (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction + }; + fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1); + } - if(fbJetsBeforeCut2){ - Double_t jetBeforeCutEntries2[10] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], - (Double_t)jetEta[0], (Double_t)jetEta[1], - (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR - }; - fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2); - } + if(fbJetsBeforeCut2){ + Double_t jetBeforeCutEntries2[10] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], + (Double_t)jetEta[0], (Double_t)jetEta[1], + (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR + }; + fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2); + } - // jet selection - Bool_t jetAccepted = kTRUE; - // minimum fraction required - if(fractionFill(4); - fh2JetSelection->Fill(4.,jetPt[0]); - jetAccepted = kFALSE; - } + // jet selection + Bool_t jetAccepted = kTRUE; + // minimum fraction required + if(fractionFill(4); + fh2JetSelection->Fill(4.,jetPt[0]); + jetAccepted = kFALSE; + } - if(jetAccepted){ - // jet acceptance + minimum pT check - if(jetEta[0]>fJetEtaMax || jetEta[0]fJetEtaMax || jetEta[1]fJetEtaMax || jetEta[0]fJetEtaMax || jetEta[1]Fill(5); - fh2JetSelection->Fill(5.,jetPt[0]); - jetAccepted = kFALSE; - } + if(fDebug){ + Printf("Jet not in eta acceptance."); + Printf("[0]: jet %d eta %.2f", ig, jetEta[0]); + Printf("[1]: jet %d eta %.2f", ir, jetEta[1]); + } + fHistJetSelection->Fill(5); + fh2JetSelection->Fill(5.,jetPt[0]); + jetAccepted = kFALSE; } - if(jetAccepted){ - if(jetPt[1] < fJetPtMin){ - if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]); - fHistJetSelection->Fill(6); - fh2JetSelection->Fill(6.,jetPt[0]); - jetAccepted = kFALSE; - } + } + if(jetAccepted){ + if(jetPt[0] < fJetPtMin){ + if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]); + fHistJetSelection->Fill(6); + fh2JetSelection->Fill(6.,jetPt[0]); + jetAccepted = kFALSE; } + } - if(jetAccepted){ - if(jet[1]->Trigger()&fJetTriggerExcludeMask){ - fHistJetSelection->Fill(7); - fh2JetSelection->Fill(7.,jetPt[0]); - jetAccepted = kFALSE; - } - + if(jetAccepted){ + if(jet[1]->Trigger()&fJetTriggerExcludeMask){ + fHistJetSelection->Fill(7); + fh2JetSelection->Fill(7.,jetPt[0]); + jetAccepted = kFALSE; } + + } - if(!jetAccepted){ - if(fbJetsMismatch2){ - Double_t jetEntriesMismatch2[10] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], - (Double_t)jetEta[0], (Double_t)jetEta[1], - (Double_t)deltaEta, (Double_t)deltaR, - (Double_t)fraction - }; - fhnJetsMismatch2->Fill(jetEntriesMismatch2); - } - continue; + if(!jetAccepted){ + if(fbJetsMismatch2){ + Double_t jetEntriesMismatch2[10] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], + (Double_t)jetEta[0], (Double_t)jetEta[1], + (Double_t)deltaEta, (Double_t)deltaR, + (Double_t)fraction + }; + fhnJetsMismatch2->Fill(jetEntriesMismatch2); } + continue; + } - // all accepted jets - fHistJetSelection->Fill(0); - fh2JetSelection->Fill(0.,jetPt[0]); + // all accepted jets + fHistJetSelection->Fill(0); + fh2JetSelection->Fill(0.,jetPt[0]); - // fill thnsparse - if(fbJetsRp){ - Double_t jetEntriesRp[11] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp, - (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], - (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr, - (Double_t)pthardbin - }; - fhnJetsRp->Fill(jetEntriesRp); - } + // fill thnsparse + if(fbJetsRp){ + Double_t jetEntriesRp[11] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp, + (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], + (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr, + (Double_t)pthardbin + }; + fhnJetsRp->Fill(jetEntriesRp); + } - if(fbJetsDeltaPt){ - Double_t jetEntriesDeltaPt[8] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea, - (Double_t)pthardbin - }; - fhnJetsDeltaPt->Fill(jetEntriesDeltaPt); - } + if(fbJetsDeltaPt){ + Double_t jetEntriesDeltaPt[8] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea, + (Double_t)pthardbin + }; + fhnJetsDeltaPt->Fill(jetEntriesDeltaPt); + } - if(fbJetsEta){ - Double_t jetEntriesEta[11] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], - (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR, - (Double_t)pthardbin - }; - fhnJetsEta->Fill(jetEntriesEta); - } + if(fbJetsEta){ + Double_t jetEntriesEta[11] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], + (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR, + (Double_t)pthardbin + }; + fhnJetsEta->Fill(jetEntriesEta); + } - if(fbJetsPhi){ - Double_t jetEntriesPhi[10] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1], - (Double_t)deltaPt, (Double_t)deltaPhi, - (Double_t)pthardbin - }; - fhnJetsPhi->Fill(jetEntriesPhi); - } + if(fbJetsPhi){ + Double_t jetEntriesPhi[10] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1], + (Double_t)deltaPt, (Double_t)deltaPhi, + (Double_t)pthardbin + }; + fhnJetsPhi->Fill(jetEntriesPhi); + } - if(fbJetsArea){ - Double_t jetEntriesArea[14] = { - (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, - (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1], - (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea, - (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet, - (Double_t)pthardbin - }; - fhnJetsArea->Fill(jetEntriesArea); - } + if(fbJetsArea){ + Double_t jetEntriesArea[14] = { + (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, + (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1], + (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea, + (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet, + (Double_t)pthardbin + }; + fhnJetsArea->Fill(jetEntriesArea); + } + + if(fbJets3Branches){ + Double_t jetEntries3Branches[14] = { + (Double_t)centValue, (Double_t)nInputTracks, + (Double_t)jetPt[0], (Double_t)jetPt[1], + (Double_t)jetArea[0], (Double_t)jetArea[1], + (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin, + (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR + }; + fhnJets3Branches->Fill(jetEntries3Branches); + } - } + } - PostData(1, fOutputList); + PostData(1, fOutputList); } void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *) { - // Draw result to the screen - // Called once at the end of the query + // Draw result to the screen + // Called once at the end of the query - if (!GetOutputData(1)) - return; + if (!GetOutputData(1)) + return; } Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks() { - // number of tracks in the event, obtained via jet finder - - Int_t nInputTracks = 0; - - TString jbname(fJetBranchName[1]); - //needs complete event, use jets without background subtraction - for(Int_t i=1; i<=3; ++i){ - if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0"); - } - // use only HI event - if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD"); - if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD"); - - if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data()); - TClonesArray *tmpAODjets = dynamic_cast(fAOD->FindListObject(jbname.Data())); - if(!tmpAODjets){ - Printf("Jet branch %s not found", jbname.Data()); - Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED"); - return -1; - } + // number of tracks in the event, obtained via jet finder + + Int_t nInputTracks = 0; + + TString jbname(fJetBranchName[1]); + //needs complete event, use jets without background subtraction + for(Int_t i=1; i<=3; ++i){ + if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0"); + } + // use only HI event + if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD"); + else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD"); + if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD"); + else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD"); + + if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data()); + TClonesArray *tmpAODjets = dynamic_cast(fAOD->FindListObject(jbname.Data())); + if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast(fAODOut->FindListObject(jbname.Data())); + if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast(fAODExtension->GetAOD()->FindListObject(jbname.Data())); + if(!tmpAODjets){ + Printf("Jet branch %s not found", jbname.Data()); + Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED"); + return -1; + } - for (Int_t iJet=0; iJetGetEntriesFast(); iJet++){ - AliAODJet *jet = dynamic_cast((*tmpAODjets)[iJet]); - if(!jet) continue; - TRefArray *trackList = jet->GetRefTracks(); - Int_t nTracks = trackList->GetEntriesFast(); - nInputTracks += nTracks; - if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks); - } - if(fDebug) Printf("---> input tracks: %d", nInputTracks); - - return nInputTracks; + for (Int_t iJet=0; iJetGetEntriesFast(); iJet++){ + AliAODJet *jet = dynamic_cast((*tmpAODjets)[iJet]); + if(!jet) continue; + TRefArray *trackList = jet->GetRefTracks(); + Int_t nTracks = trackList->GetEntriesFast(); + nInputTracks += nTracks; + if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks); + } + if(fDebug) Printf("---> input tracks: %d", nInputTracks); + + return nInputTracks; } THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt) { - // generate new THnSparseF, axes are defined in GetDimParams() - - Int_t count = 0; - UInt_t tmp = entries; - while(tmp!=0){ - count++; - tmp = tmp &~ -tmp; // clear lowest bit - } - - TString hnTitle(name); - const Int_t dim = count; - Int_t nbins[dim]; - Double_t xmin[dim]; - Double_t xmax[dim]; - - Int_t i=0; - Int_t c=0; - while(cGetPtSubtracted(0); - } - else{ - pt = j->Pt(); // if negative, pt=0.01 - } + if(fKeepJets && mode>0){ // background subtracted pt, can be negative + pt = j->GetPtSubtracted(0); + } + else{ + pt = j->Pt(); // if negative, pt=0.01 + } - return pt; + return pt; }