/************************************************************************** * 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. * **************************************************************************/ // // task compares jets in two or three branches, // written for analysis of jet embedding in HI events // #include "TChain.h" #include "TTree.h" #include "TMath.h" #include "TH1F.h" #include "TH1I.h" #include "TH2F.h" #include "TH3F.h" #include "THnSparse.h" #include "TCanvas.h" #include "AliLog.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliVEvent.h" #include "AliESDEvent.h" #include "AliESDInputHandler.h" #include "AliCentrality.h" #include "AliAnalysisHelperJetTasks.h" #include "AliInputEventHandler.h" #include "AliAODJetEventBackground.h" #include "AliAnalysisTaskFastEmbedding.h" #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), 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 fJetBranchName[0] = ""; fJetBranchName[1] = ""; fJetBranchName[2] = ""; fListJets[0] = new TList; fListJets[1] = new TList; fListJets[2] = new TList; } AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) : 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 fJetBranchName[0] = ""; fJetBranchName[1] = ""; fJetBranchName[2] = ""; fListJets[0] = new TList; fListJets[1] = new TList; fListJets[2] = new TList; DefineOutput(1, TList::Class()); } AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2() { delete fListJets[0]; delete fListJets[1]; delete fListJets[2]; } void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3) { 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."); } } 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); } 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); } // 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()); // 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 && 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(); } // 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 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]; 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]; 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(!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(); } } //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.; //Float_t relRho = rho>0. ? localRho / rho : 0.; // calculate parameters of associated jets /* from Leticia, kT clusterizer Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 }; Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 }; Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 }; */ // own, from embedded tracks Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 }; Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 }; Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 }; Float_t rpCorr = 0.; if(eventClass>0&&eventClass<4){ rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]); rpCorr *= rho * jetArea[1]; } Float_t 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); } 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; } if(jetAccepted){ // jet acceptance + minimum pT check if(jetEta[0]>fJetEtaMax || jetEta[0]fJetEtaMax || jetEta[1]Fill(5); fh2JetSelection->Fill(5.,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(fbJetsMismatch2){ Double_t jetEntriesMismatch2[10] = { (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], (Double_t)deltaEta, (Double_t)deltaR, (Double_t)fraction }; fhnJetsMismatch2->Fill(jetEntriesMismatch2); } continue; } // all accepted jets fHistJetSelection->Fill(0); fh2JetSelection->Fill(0.,jetPt[0]); // fill thnsparse if(fbJetsRp){ Double_t jetEntriesRp[11] = { (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp, (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr, (Double_t)pthardbin }; fhnJetsRp->Fill(jetEntriesRp); } if(fbJetsDeltaPt){ Double_t jetEntriesDeltaPt[8] = { (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea, (Double_t)pthardbin }; fhnJetsDeltaPt->Fill(jetEntriesDeltaPt); } if(fbJetsEta){ Double_t jetEntriesEta[11] = { (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR, (Double_t)pthardbin }; fhnJetsEta->Fill(jetEntriesEta); } if(fbJetsPhi){ Double_t jetEntriesPhi[10] = { (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)deltaPt, (Double_t)deltaPhi, (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(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); } void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *) { // Draw result to the screen // Called once at the end of the query 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"); 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; } 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(c0){ // background subtracted pt, can be negative pt = j->GetPtSubtracted(0); } else{ pt = j->Pt(); // if negative, pt=0.01 } return pt; }