--- /dev/null
+
+// ******************************************
+// This task computes several jet observables like
+// the fraction of energy in inner and outer coronnas,
+// jet-track correlations,triggered jet shapes and
+// correlation strength distribution of particles inside jets.
+// Author: lcunquei@cern.ch
+// *******************************************
+
+
+/**************************************************************************
+ * 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 "TChain.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH1D.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 "AliAODMCParticle.h"
+//#include "AliAnalysisTaskFastEmbedding.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODJet.h"
+
+#include "AliAnalysisTaskJetCorePP.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskJetCorePP)
+
+//Filip Krizek 1st March 2013
+
+//---------------------------------------------------------------------
+AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
+AliAnalysisTaskSE(),
+fESD(0x0),
+fAODIn(0x0),
+fAODOut(0x0),
+fAODExtension(0x0),
+fJetBranchName(""),
+//fListJets(0x0),
+fNonStdFile(""),
+fSystem(0), //pp=0 pPb=1
+fJetParamR(0.4),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-10.0),
+fVtxZMax(10.0),
+fFilterMask(0),
+fCentMin(0.0),
+fCentMax(100.0),
+fJetEtaMin(-0.5),
+fJetEtaMax(0.5),
+fTriggerEtaCut(0.9),
+fTrackEtaCut(0.9),
+fTrackLowPtCut(0.15),
+fOutputList(0x0),
+fHistEvtSelection(0x0),
+fh2Ntriggers(0x0),
+fHJetSpec(0x0),
+fHJetDensity(0x0),
+fHJetDensityA4(0x0),
+fhJetPhi(0x0),
+fhTriggerPhi(0x0),
+fhJetEta(0x0),
+fhTriggerEta(0x0),
+fhVertexZ(0x0),
+fhVertexZAccept(0x0),
+fhContribVtx(0x0),
+fhContribVtxAccept(0x0),
+fhDphiTriggerJet(0x0),
+fhDphiTriggerJetAccept(0x0),
+fhCentrality(0x0),
+fhCentralityAccept(0x0),
+fkAcceptance(2.0*TMath::Pi()*1.8)
+{
+ // default Constructor
+ fListJets = new TList();
+}
+
+//---------------------------------------------------------------------
+
+AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
+AliAnalysisTaskSE(name),
+fESD(0x0),
+fAODIn(0x0),
+fAODOut(0x0),
+fAODExtension(0x0),
+fJetBranchName(""),
+//fListJets(0x0),
+fNonStdFile(""),
+fSystem(0), //pp=0 pPb=1
+fJetParamR(0.4),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-10.0),
+fVtxZMax(10.0),
+fFilterMask(0),
+fCentMin(0.0),
+fCentMax(100.0),
+fJetEtaMin(-0.5),
+fJetEtaMax(0.5),
+fTriggerEtaCut(0.9),
+fTrackEtaCut(0.9),
+fTrackLowPtCut(0.15),
+fOutputList(0x0),
+fHistEvtSelection(0x0),
+fh2Ntriggers(0x0),
+fHJetSpec(0x0),
+fHJetDensity(0x0),
+fHJetDensityA4(0x0),
+fhJetPhi(0x0),
+fhTriggerPhi(0x0),
+fhJetEta(0x0),
+fhTriggerEta(0x0),
+fhVertexZ(0x0),
+fhVertexZAccept(0x0),
+fhContribVtx(0x0),
+fhContribVtxAccept(0x0),
+fhDphiTriggerJet(0x0),
+fhDphiTriggerJetAccept(0x0),
+fhCentrality(0x0),
+fhCentralityAccept(0x0),
+fkAcceptance(2.0*TMath::Pi()*1.8)
+{
+ // Constructor
+ fListJets = new TList();
+
+ DefineOutput(1, TList::Class());
+}
+
+//--------------------------------------------------------------
+AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
+AliAnalysisTaskSE(a.GetName()),
+fESD(a.fESD),
+fAODIn(a.fAODIn),
+fAODOut(a.fAODOut),
+fAODExtension(a.fAODExtension),
+fJetBranchName(a.fJetBranchName),
+fListJets(a.fListJets),
+fNonStdFile(a.fNonStdFile),
+fSystem(a.fSystem),
+fJetParamR(a.fJetParamR),
+fOfflineTrgMask(a.fOfflineTrgMask),
+fMinContribVtx(a.fMinContribVtx),
+fVtxZMin(a.fVtxZMin),
+fVtxZMax(a.fVtxZMax),
+fFilterMask(a.fFilterMask),
+fCentMin(a.fCentMin),
+fCentMax(a.fCentMax),
+fJetEtaMin(a.fJetEtaMin),
+fJetEtaMax(a.fJetEtaMax),
+fTriggerEtaCut(a.fTriggerEtaCut),
+fTrackEtaCut(a.fTrackEtaCut),
+fTrackLowPtCut(a.fTrackLowPtCut),
+fOutputList(a.fOutputList),
+fHistEvtSelection(a.fHistEvtSelection),
+fh2Ntriggers(a.fh2Ntriggers),
+fHJetSpec(a.fHJetSpec),
+fHJetDensity(a.fHJetDensity),
+fHJetDensityA4(a.fHJetDensityA4),
+fhJetPhi(a.fhJetPhi),
+fhTriggerPhi(a.fhTriggerPhi),
+fhJetEta(a.fhJetEta),
+fhTriggerEta(a.fhTriggerEta),
+fhVertexZ(a.fhVertexZ),
+fhVertexZAccept(a.fhVertexZAccept),
+fhContribVtx(a.fhContribVtx),
+fhContribVtxAccept(a.fhContribVtxAccept),
+fhDphiTriggerJet(a.fhDphiTriggerJet),
+fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
+fhCentrality(a.fhCentrality),
+fhCentralityAccept(a.fhCentralityAccept),
+fkAcceptance(a.fkAcceptance)
+{
+ //Copy Constructor
+}
+//--------------------------------------------------------------
+
+AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
+ // assignment operator
+ this->~AliAnalysisTaskJetCorePP();
+ new(this) AliAnalysisTaskJetCorePP(a);
+ return *this;
+}
+//--------------------------------------------------------------
+
+AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
+{
+ //Destructor
+ delete fListJets;
+ delete fOutputList; // ????
+}
+
+//--------------------------------------------------------------
+
+void AliAnalysisTaskJetCorePP::Init()
+{
+ // check for jet branches
+ if(!strlen(fJetBranchName.Data())){
+ AliError("Jet branch name not set.");
+ }
+
+}
+
+//--------------------------------------------------------------
+
+void AliAnalysisTaskJetCorePP::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)");
+
+ fOutputList->Add(fHistEvtSelection);
+
+ Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
+
+ fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
+ nBinsCentrality,0.0,100.0,50,0.0,50.0);
+
+ //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
+ const Int_t dimSpec = 6;
+ const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 140, 50, 100, 50};
+ const Double_t lowBinSpec[dimSpec] = {0.0, 0.0,-80.0, 0.0, 0.0, 0.0};
+ const Double_t hiBinSpec[dimSpec] = {100.0, 1.0,200.0,50.0, 200, 20.0};
+ fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum [cent,A,pT-rho*A,pTtrig]",
+ dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
+ fOutputList->Add(fHJetSpec);
+
+ //------------------- HISTOS FOR DIAGNOSTIC ----------------------
+ //Jet number density histos [Trk Mult, jet density, pT trigger]
+ const Int_t dimJetDens = 3;
+ const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
+ const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
+ const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0 };
+
+ fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
+ dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
+
+ fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
+ dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
+
+ fOutputList->Add(fh2Ntriggers);
+ fOutputList->Add(fHJetDensity);
+ fOutputList->Add(fHJetDensityA4);
+
+
+ //inclusive azimuthal and pseudorapidity histograms
+ fhJetPhi = new TH1D("fhJetPhi","Azim dist jets",50,-TMath::Pi(),TMath::Pi());
+ fhTriggerPhi= new TH1D("fhTriggerPhi","azim dist trig had",50,-TMath::Pi(),TMath::Pi());
+ fhJetEta = new TH1D("fhJetEta","Eta dist jets",40,-0.9,0.9);
+ fhTriggerEta = new TH1D("fhTriggerEta","Eta dist trig had",40,-0.9,0.9);
+ fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
+ fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
+ fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
+ fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
+ fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
+ fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
+ fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
+ fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
+
+ fOutputList->Add(fhJetPhi);
+ fOutputList->Add(fhTriggerPhi);
+ fOutputList->Add(fhJetEta);
+ fOutputList->Add(fhTriggerEta);
+ fOutputList->Add(fhVertexZ);
+ fOutputList->Add(fhVertexZAccept);
+ fOutputList->Add(fhContribVtx);
+ fOutputList->Add(fhContribVtxAccept);
+ fOutputList->Add(fhDphiTriggerJet);
+ fOutputList->Add(fhDphiTriggerJetAccept);
+ fOutputList->Add(fhCentrality);
+ fOutputList->Add(fhCentralityAccept);
+
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for(Int_t i=0; i<fOutputList->GetEntries(); i++){
+ TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
+ if(h1){
+ h1->Sumw2();
+ continue;
+ }
+ THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
+ if(hn){
+ hn->Sumw2();
+ }
+ }
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1, fOutputList);
+}
+
+//--------------------------------------------------------------------
+
+void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
+{
+
+ //Event loop
+
+ if(TMath::Abs((Float_t) fJetParamR)<0.00001){
+ AliError("Cone radius is set to zero.");
+ return;
+ }
+ if(!strlen(fJetBranchName.Data())){
+ AliError("Jet branch name not set.");
+ return;
+ }
+
+ fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+ if(!fESD){
+ AliError("ESD not available");
+ fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
+ }
+
+ fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
+ AliAODEvent* aod = NULL;
+ // take all other information from the aod we take the tracks from
+ if(!aod){
+ if(!fESD) aod = fAODIn;
+ else aod = fAODOut;
+ }
+
+ if(fNonStdFile.Length()!=0){
+ // case that we have an AOD extension we can fetch the jets from the extended output
+ AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(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;
+ }
+
+ //check AOD pointer
+ if(!aod){
+ if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }
+
+ // vertex selection
+ AliAODVertex* primVtx = aod->GetPrimaryVertex();
+
+ if(!primVtx){
+ if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }
+
+ Int_t nTracksPrim = primVtx->GetNContributors();
+ Float_t vtxz = primVtx->GetZ();
+ //Input events
+ fhContribVtx->Fill(nTracksPrim);
+ if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
+
+ if((nTracksPrim < fMinContribVtx) ||
+ (vtxz < fVtxZMin) ||
+ (vtxz > fVtxZMax)){
+ if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
+ (char*)__FILE__,__LINE__,vtxz);
+ fHistEvtSelection->Fill(3);
+ PostData(1, fOutputList);
+ return;
+ }else{
+ //Accepted events
+ fhContribVtxAccept->Fill(nTracksPrim);
+ fhVertexZAccept->Fill(vtxz);
+ }
+ //FK// No event class selection imposed
+ // 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;
+ Double_t centValue = 0.0;
+ if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
+ if(fESD){
+ cent = fESD->GetCentrality();
+ if(cent) centValue = cent->GetCentralityPercentile("V0M");
+ }else{
+ centValue = aod->GetHeader()->GetCentrality();
+ }
+ if(fDebug) printf("centrality: %f\n", centValue);
+ //Input events
+ fhCentrality->Fill(centValue);
+
+ if(centValue < fCentMin || centValue > fCentMax){
+ fHistEvtSelection->Fill(4);
+ PostData(1, fOutputList);
+ return;
+ }else{
+ //Accepted events
+ fhCentralityAccept->Fill( centValue );
+ }
+ }
+
+ if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
+
+ fHistEvtSelection->Fill(0); // accepted events
+
+ // ------------------- end event selection --------------------
+
+ // fetch jets
+ TClonesArray *aodJets = 0x0;
+
+ if(fAODOut && !aodJets){
+ aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
+ }
+ if(fAODExtension && !aodJets){
+ aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
+ }
+ if(fAODIn && !aodJets){
+ aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
+ }
+
+ // ------------- Hadron trigger --------------
+ TList particleList; //list of tracks
+ Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
+
+ if(indexTrigg<0) return; // no trigger track found
+
+ AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
+ if(!triggerHadron){
+ PostData(1, fOutputList);
+ return;
+ }
+
+
+ fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
+
+ // if(triggerHadron->Pt() > 10.0){}
+ if(triggerHadron->Pt() > 3.0){
+ //Trigger Diagnostics---------------------------------
+ fhTriggerPhi->Fill(RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
+ fhTriggerEta->Fill(triggerHadron->Eta());
+ }
+
+ //--------- Fill list of jets --------------
+ fListJets->Clear();
+ if(aodJets){
+ if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(),
+ aodJets->GetEntriesFast());
+ for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
+ AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
+ if (jet) fListJets->Add(jet);
+ }
+ }
+
+ Double_t etaJet = 0.0;
+ Double_t pTJet = 0.0;
+ Double_t areaJet = 0.0;
+ Double_t phiJet = 0.0;
+ Int_t injet4 = 0;
+ Int_t injet = 0;
+ //---------- jet loop ---------
+ for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
+ AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
+ if(!jet) continue;
+ etaJet = jet->Eta();
+ phiJet = jet->Phi();
+ pTJet = jet->Pt();
+ if(pTJet==0) continue;
+
+ if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
+ areaJet = jet->EffectiveAreaCharged();
+
+ //Jet Diagnostics---------------------------------
+ fhJetPhi->Fill(RelativePhi(phiJet,0.0)); //phi -pi,pi
+ fhJetEta->Fill(etaJet);
+ if(areaJet >= 0.07) injet++;
+ if(areaJet >= 0.4) injet4++;
+ //--------------------------------------------------
+
+ Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
+
+ fhDphiTriggerJet->Fill(dphi); //Input
+ if(TMath::Abs((Double_t) dphi) < TMath::Pi()-0.6) continue;
+ fhDphiTriggerJetAccept->Fill(dphi); //Accepted
+
+ //Background w.r.t. jet axis
+ Double_t pTBckWrtJet =
+ GetBackgroundInPerpCone(fJetParamR, phiJet, etaJet, &particleList);
+
+ Double_t ratioOfAreas = areaJet/(TMath::Pi()*fJetParamR*fJetParamR);
+ Double_t rhoA = ratioOfAreas*pTBckWrtJet; //bg activity in a cone of similar size
+ Double_t ptcorr = pTJet - rhoA; //Pt Jet UE subtr
+
+
+ //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
+ Double_t fillspec[] = { centValue,
+ areaJet,
+ ptcorr,
+ triggerHadron->Pt(),
+ pTJet,
+ rhoA
+ };
+ fHJetSpec->Fill(fillspec);
+
+ }//jet loop
+
+ //Fill Jet Density In the Event A>0.07
+ if(injet>0){
+ Double_t filldens[]={ (Double_t) particleList.GetEntries(),
+ injet/fkAcceptance,
+ triggerHadron->Pt()
+ };
+ fHJetDensity->Fill(filldens);
+ }
+
+ //Fill Jet Density In the Event A>0.4
+ if(injet4>0){
+ Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
+ injet4/fkAcceptance,
+ triggerHadron->Pt()
+ };
+ fHJetDensityA4->Fill(filldens4);
+ }
+
+ PostData(1, fOutputList);
+}
+
+//----------------------------------------------------------------------------
+void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
+{
+ // Draw result to the screen
+ // Called once at the end of the query
+
+ if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
+ if(!GetOutputData(1)) return;
+}
+
+
+//----------------------------------------------------------------------------
+Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
+ //Fill the list of accepted tracks (passed track cut)
+ //return consecutive index of the hardest ch hadron in the list
+ Int_t iCount = 0;
+ AliAODEvent *aodevt = NULL;
+
+ if(!fESD) aodevt = fAODIn;
+ else aodevt = fAODOut;
+
+ if(!aodevt) return -1;
+
+ Int_t index = -1; //index of the highest particle in the list
+ Double_t ptmax = -10;
+
+ for(int it = 0; it < aodevt->GetNumberOfTracks(); it++){
+ AliAODTrack *tr = aodevt->GetTrack(it);
+
+ //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
+ if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
+ if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
+ if(tr->Pt() < fTrackLowPtCut) continue;
+ list->Add(tr);
+ if(tr->Pt()>ptmax){
+ ptmax = tr->Pt();
+ index = iCount;
+ }
+ iCount++;
+ }
+
+ if(index>-1){ //check pseudorapidity cut on trigger
+ AliAODTrack *trigger = (AliAODTrack*) list->At(index);
+ if(trigger && TMath::Abs((Float_t) trigger->Eta())< fTriggerEtaCut){ return index;}
+ return -1;
+ }else{
+ return -1;
+ }
+}
+
+//----------------------------------------------------------------------------
+
+Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
+ //calculate sum of track pT in the cone perpendicular in phi to the jet
+ //jetR = cone radius
+ // jetPhi, jetEta = direction of the jet
+ Int_t numberOfTrks = trkList->GetEntries();
+ Double_t pTsum = 0.0;
+ Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
+ for(Int_t it=0; it<numberOfTrks; it++){
+ AliVParticle *trk = (AliVParticle*) trkList->At(it);
+ Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
+ Double_t deta = trk->Eta()-jetEta;
+
+ if( (dphi*dphi + deta*deta)< (jetR*jetR)){
+ pTsum += trk->Pt();
+ }
+ }
+
+ return pTsum;
+}
+
+//----------------------------------------------------------------------------
+
+Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
+ //Get relative azimuthal angle of two particles -pi to pi
+ if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
+ else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
+
+ if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
+ else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
+
+ Double_t dphi = mphi - vphi;
+ if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
+ else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
+
+ return dphi;//dphi in [-Pi, Pi]
+}
+
+
--- /dev/null
+#ifndef ALIANALYSISTASKJETCOREPP_H
+#define ALIANALYSISTASKJETCOREPP_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// **************************************
+// This task performs hadron-trigger recoil jet correlations
+// Output pT spectrum of jet given trigger pT
+// Author: filip krizek 1st March 2013
+// *******************************************
+
+class TH1F;
+class TH1D;
+class TH1I;
+class TH2F;
+class TH3F;
+class THnSparse;
+class AliESDEvent;
+class AliAODExtension;
+class AliAODEvent;
+
+#include "AliAnalysisTaskSE.h"
+#include "AliVEvent.h"
+
+class AliAnalysisTaskJetCorePP : public AliAnalysisTaskSE {
+public:
+ AliAnalysisTaskJetCorePP();
+ AliAnalysisTaskJetCorePP(const char *name);
+ AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a);
+ AliAnalysisTaskJetCorePP& operator=(const AliAnalysisTaskJetCorePP& a); // not implemented
+ virtual ~AliAnalysisTaskJetCorePP();
+ virtual void LocalInit() {Init();}
+ virtual void Init();
+ virtual void UserCreateOutputObjects();
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(const Option_t*);
+
+ virtual void SetBranchName(const TString &name){ fJetBranchName = name; }
+ virtual void SetNonStdFile(char* c){fNonStdFile = c;}
+ virtual void SetSystem(Int_t sys) { fSystem = sys; }
+ virtual void SetJetR(Float_t jR) { fJetParamR = jR; }
+ virtual void SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+ virtual void SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
+ virtual void SetVtxZMin(Float_t z) { fVtxZMin = z; }
+ virtual void SetVtxZMax(Float_t z) { fVtxZMax = z; }
+ virtual void SetFilterMask(UInt_t i){fFilterMask = i;}
+ virtual void SetCentMin(Float_t cent) { fCentMin = cent; }
+ virtual void SetCentMax(Float_t cent) { fCentMax = cent; }
+ virtual void SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; }
+ virtual void SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; }
+ virtual void SetTriggerEtaCut(Float_t eta) { fTriggerEtaCut = eta; }
+ virtual void SetTrackEtaCut(Float_t eta) { fTrackEtaCut = eta; }
+ virtual void SetTrackLowPtCut(Float_t pt) { fTrackLowPtCut=pt; }
+
+ Double_t RelativePhi(Double_t angle1, Double_t angle2);
+
+private:
+ //private member functions
+ Int_t GetListOfTracks(TList *list); //returns index of trig and track list
+ Double_t GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList); //sums pT in the cone perp in phi to jet
+
+ //private member objects
+ AliESDEvent *fESD; //! ESD object
+ AliAODEvent *fAODIn; //! AOD event for AOD input tracks
+ AliAODEvent *fAODOut; //! AOD event
+ AliAODExtension *fAODExtension; //! where we take the jets from can be input or output AOD
+
+ // jets to compare
+ TString fJetBranchName; // name of jet branch
+ TList *fListJets; //! jet lists
+
+ TString fNonStdFile; // name of delta aod file to catch the extension
+
+ // event selection
+ Int_t fSystem; // collision system pp=0, pPb=1
+ Float_t fJetParamR; // jet cone resolution (radius) R
+ AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline trigs
+ Int_t fMinContribVtx; // min numb of trk contrib for prim vertex
+ Float_t fVtxZMin; // lower bound on vertex z
+ Float_t fVtxZMax; // upper bound on vertex z
+ UInt_t fFilterMask; // filter bit for slected tracks
+ Float_t fCentMin; // lower bound on centrality
+ Float_t fCentMax; // upper bound on centrality
+ Float_t fJetEtaMin; // lower bound on eta for found jets
+ Float_t fJetEtaMax; // upper bound on eta for found jets
+ Float_t fTriggerEtaCut; // lower bound on eta for trigger track
+ Float_t fTrackEtaCut; // upper bound on eta for trigger track
+ Float_t fTrackLowPtCut; // upper bound on eta for trigger track
+
+
+ TList *fOutputList; //! output data container
+ TH1I *fHistEvtSelection; //! event selection statistic
+ TH2F *fh2Ntriggers; //trigger pT versus centrality
+ THnSparse *fHJetSpec; //Recoil jet spectrum
+
+ //Diagnostics
+ THnSparse *fHJetDensity; //density of jet with A>0.07 //fk
+ THnSparse *fHJetDensityA4; //density of jets with A>0.4 //fk
+ TH1D *fhJetPhi; //Azimuthal distribution of jets
+ TH1D *fhTriggerPhi; //Azimuthal distribution of trigger hadron
+ TH1D *fhJetEta; //Pseudorapidity distribution of jets
+ TH1D *fhTriggerEta; //Pseudorapidity distribution of trigger hadron
+ TH1D *fhVertexZ; //z vertex distribution
+ TH1D *fhVertexZAccept; //z vertex distribution after cut
+ TH1D *fhContribVtx; //contributors to vertex
+ TH1D *fhContribVtxAccept; //contributors to vertex after cut
+ TH1D *fhDphiTriggerJet; //Deltaphi between trigger and jet
+ TH1D *fhDphiTriggerJetAccept; //Deltaphi between trigger and jet after cut
+ TH1D *fhCentrality; //Deltaphi between trigger and jet
+ TH1D *fhCentralityAccept; //Deltaphi between trigger and jet after cut
+
+ const Double_t fkAcceptance; //eta times phi Alice coverage
+
+ ClassDef(AliAnalysisTaskJetCorePP, 1); //has to end with number larger than 0
+};
+
+#endif
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff