From e59c88bd4e4d9b2f3b5f11e0d26a0d228c3e3a79 Mon Sep 17 00:00:00 2001 From: lcunquei Date: Thu, 22 May 2014 17:28:31 +0200 Subject: [PATCH] New antenna radiation task --- PWGJE/CMakelibPWGJE.pkg | 3 +- PWGJE/PWGJELinkDef.h | 1 + PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx | 629 ++++++++++++++++++ PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h | 168 +++++ PWGJE/macros/AddTaskJetAntenna.C | 52 ++ 5 files changed, 852 insertions(+), 1 deletion(-) create mode 100644 PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx create mode 100644 PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h create mode 100644 PWGJE/macros/AddTaskJetAntenna.C diff --git a/PWGJE/CMakelibPWGJE.pkg b/PWGJE/CMakelibPWGJE.pkg index 40c1f035455..5446c08d694 100644 --- a/PWGJE/CMakelibPWGJE.pkg +++ b/PWGJE/CMakelibPWGJE.pkg @@ -53,7 +53,8 @@ set ( SRCS UserTasks/AliAnalysisTaskIDFFTCF.cxx UserTasks/AliIDFFUtils.cxx UserTasks/AliAnalysisTaskPPJetSpectra.cxx - ) + UserTasks/AliAnalysisTaskJetAntenna.cxx + ) # Add code that needs fastjet or FJWrapper here if (FASTJET_FOUND) diff --git a/PWGJE/PWGJELinkDef.h b/PWGJE/PWGJELinkDef.h index 2ce323adb99..320810645d0 100644 --- a/PWGJE/PWGJELinkDef.h +++ b/PWGJE/PWGJELinkDef.h @@ -44,6 +44,7 @@ #pragma link C++ class AliAnalysisTaskIDFFTCF::AliFragFuncQAJetHistos+; #pragma link C++ class AliIDFFUtils+; #pragma link C++ class AliAnalysisTaskPPJetSpectra+; +#pragma link C++ class AliAnalysisTaskJetAntenna+; #ifdef HAVE_FASTJET #pragma link C++ class AliAnalysisTaskCheckSingleTrackJetRejection+; diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx b/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx new file mode 100644 index 00000000000..bb1c42ec904 --- /dev/null +++ b/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx @@ -0,0 +1,629 @@ +// ****************************************** +// 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 "TH2F.h" +#include "TH3F.h" +#include "THnSparse.h" +#include "TMatrixD.h" +#include "TMatrixDSym.h" +#include "TMatrixDSymEigen.h" +#include "TCanvas.h" +#include "TRandom3.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 "AliAnalysisTaskJetAntenna.h" + +using std::cout; +using std::endl; + +ClassImp(AliAnalysisTaskJetAntenna) + +AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna() : +AliAnalysisTaskSE(), +fESD(0x0), +fAODIn(0x0), +fAODOut(0x0), +fAODExtension(0x0), +fBackgroundBranch(""), +fNonStdFile(""), +fIsPbPb(kTRUE), +fOfflineTrgMask(AliVEvent::kAny), +fMinContribVtx(1), +fVtxZMin(-10.), +fVtxZMax(10.), +fEvtClassMin(0), +fEvtClassMax(4), +fFilterMask(0), +fFilterMaskBestPt(0), +fFilterType(0), +fCentMin(0.), +fCentMax(100.), +fNInputTracksMin(0), +fNInputTracksMax(-1), +fRequireITSRefit(0), +fApplySharedClusterCut(0), +fTrackTypeRec(kTrackUndef), +fRPAngle(0), +fNRPBins(50), +fSemigoodCorrect(0), +fHolePos(4.71), +fHoleWidth(0.2), +fCutTM(0.15), +fJetEtaMin(-.5), +fJetEtaMax(.5), +fNevents(0), +fTindex(0), +fJetPtMin(20.), +fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), +fJetPtFractionMin(0.5), +fNMatchJets(4), +fMatchMaxDist(0.8), +fKeepJets(kFALSE), +fkNbranches(2), +fkEvtClasses(12), +fOutputList(0x0), +fHistEvtSelection(0x0), +fh1JetEntries(0x0), +fh2Circularity(0x0), +fhnJetTM(0x0) +{ + // default Constructor + + fJetBranchName[0] = ""; + fJetBranchName[1] = ""; + + fListJets[0] = new TList; + fListJets[1] = new TList; +} + +AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna(const char *name) : +AliAnalysisTaskSE(name), +fESD(0x0), +fAODIn(0x0), +fAODOut(0x0), +fAODExtension(0x0), +fBackgroundBranch(""), +fNonStdFile(""), +fIsPbPb(kTRUE), +fOfflineTrgMask(AliVEvent::kAny), +fMinContribVtx(1), +fVtxZMin(-10.), +fVtxZMax(10.), +fEvtClassMin(0), +fEvtClassMax(4), +fFilterMask(0), +fFilterMaskBestPt(0), +fFilterType(0), +fCentMin(0.), +fCentMax(100.), +fNInputTracksMin(0), +fNInputTracksMax(-1), +fRequireITSRefit(0), +fApplySharedClusterCut(0), +fTrackTypeRec(kTrackUndef), +fRPAngle(0), +fNRPBins(50), +fSemigoodCorrect(0), +fHolePos(4.71), +fHoleWidth(0.2), +fCutTM(0.15), +fJetEtaMin(-.5), +fJetEtaMax(.5), +fNevents(0), +fTindex(0), +fJetPtMin(20.), +fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), +fJetPtFractionMin(0.5), +fNMatchJets(4), +fMatchMaxDist(0.8), +fKeepJets(kFALSE), +fkNbranches(2), +fkEvtClasses(12), +fOutputList(0x0), +fHistEvtSelection(0x0), +fh1JetEntries(0x0), +fh2Circularity(0x0), +fhnJetTM(0x0) + { + // Constructor + + + fJetBranchName[0] = ""; + fJetBranchName[1] = ""; + + fListJets[0] = new TList; + fListJets[1] = new TList; + + DefineOutput(1, TList::Class()); +} + +AliAnalysisTaskJetAntenna::~AliAnalysisTaskJetAntenna() +{ + delete fListJets[0]; + delete fListJets[1]; +} + +void AliAnalysisTaskJetAntenna::SetBranchNames(const TString &branch1, const TString &branch2) +{ + fJetBranchName[0] = branch1; + fJetBranchName[1] = branch2; +} + +void AliAnalysisTaskJetAntenna::Init() +{ + + // check for jet branches + if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){ + AliError("Jet branch name not set."); + } + +} + +void AliAnalysisTaskJetAntenna::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); + fh1JetEntries=new TH1F("JetEntries","",150,0,150); + fOutputList->Add(fh1JetEntries); + fh2Circularity=new TH2F("Circcularity","",10,0,1,150,0,150); + fOutputList->Add(fh2Circularity); + Int_t nbinsJet[6]={15,30,9,360,10,50}; + Double_t binlowJet[6]= {0, 0, 0,-0.5*TMath::Pi(),0,0}; + Double_t binupJet[6]= {1.5, 150,150,1.5*TMath::Pi(),1,200}; + fhnJetTM = new THnSparseF("fhnJetTM", "fhnJetTM; dr;pt_jet;pt_track;phi;",6,nbinsJet,binlowJet,binupJet); + Double_t xPt3[10]; + xPt3[0] = 0.; + for(Int_t i = 1;i<=9;i++){ + if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5 + else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12 + else xPt3[i] = xPt3[i-1] + 150.; // 18 + } + fhnJetTM->SetBinEdges(2,xPt3); + fOutputList->Add(fhnJetTM); + + // =========== 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 AliAnalysisTaskJetAntenna::UserExec(Option_t *) +{ + + + 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"); + fAODIn = dynamic_cast(InputEvent()); + } + fAODOut = dynamic_cast(AODEvent()); + + static AliAODEvent* aod = 0; + // 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 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 + + + Bool_t selected=kTRUE; + selected = AliAnalysisHelperJetTasks::Selected(); + if(!selected){ + // no selection by the service task, we continue + PostData(1,fOutputList); + return;} + + + + // physics selection: this is now redundant, all should appear as accepted after service task selection + AliInputEventHandler* inputHandler = (AliInputEventHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); + std::cout<IsEventSelected()<<" "<IsEventSelected() & fOfflineTrgMask)){ + if(fDebug) Printf(" Trigger Selection: event REJECTED ... "); + fHistEvtSelection->Fill(2); + PostData(1, fOutputList); + return; + } + + + + // vertex selection + if(!aod){ + if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); + fHistEvtSelection->Fill(3); + PostData(1, fOutputList); + } + 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(); + 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; + } + + + + // centrality selection + AliCentrality *cent = 0x0; + Double_t centValue = 0.; + if(fIsPbPb){ + if(fESD) {cent = fESD->GetCentrality(); + if(cent) centValue = cent->GetCentralityPercentile("V0M");} + else centValue=aod->GetHeader()->GetCentrality(); + + if(fDebug) printf("centrality: %f\n", centValue); + if (centValue < fCentMin || centValue > fCentMax){ + fHistEvtSelection->Fill(4); + PostData(1, fOutputList); + return; + }} + + + fHistEvtSelection->Fill(0); + // accepted events + // -- end event selection -- + + // get background + AliAODJetEventBackground* externalBackground = 0; + if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){ + externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data())); + if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; + } + if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){ + externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); + if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; + } + + if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){ + externalBackground = (AliAODJetEventBackground*)(fAODIn->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]=0; + if(fAODOut&&!aodJets[0]){ + aodJets[0] = dynamic_cast(fAODOut->FindListObject(fJetBranchName[0].Data())); + aodJets[1] = dynamic_cast(fAODOut->FindListObject(fJetBranchName[1].Data())); + } + if(fAODExtension && !aodJets[0]){ + aodJets[0] = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); + aodJets[1] = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); + } + if(fAODIn&&!aodJets[0]){ + aodJets[0] = dynamic_cast(fAODIn->FindListObject(fJetBranchName[0].Data())); + aodJets[1] = dynamic_cast(fAODIn->FindListObject(fJetBranchName[1].Data())); + } + + + + Int_t nT=0; + TList ParticleList; + nT = GetListOfTracks(&ParticleList); + if(nT<0){ + PostData(1, fOutputList); + return; + } + + 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()); + for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) { + AliAODJet *jet = dynamic_cast((*aodJets[iJetType])[iJet]); + if (jet) fListJets[iJetType]->Add(jet); + } + } + + + + for(Int_t i=0; iGetEntries(); ++i){ + + Double_t etabig=0; + Double_t ptbig=0; + Double_t areabig=0; + Double_t phibig=0.; + Double_t pxbig,pybig,pzbig; + + + AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i)); + etabig = jetbig->Eta(); + phibig = jetbig->Phi(); + ptbig = jetbig->Pt(); + if(ptbig==0) continue; + areabig = jetbig->EffectiveAreaCharged(); + ptbig=ptbig-rho*areabig; + if((etabigfJetEtaMax)) continue; + + if(fSemigoodCorrect){ + if((phibig>fHolePos-fHoleWidth) && (phibigPx(); + pybig=jetbig->Py(); + pzbig=jetbig->Pz(); + TVector3 ppJ1(pxbig, pybig, pzbig); + TVector3 ppJ3(- pxbig * pzbig, - pybig * pzbig, pxbig * pxbig + pybig * pybig); + ppJ3.SetMag(1.); + TVector3 ppJ2(-pybig, pxbig, 0); + ppJ2.SetMag(1.); + + Float_t mxx = 0.; + Float_t myy = 0.; + Float_t mxy = 0.; + Int_t nc = 0; + Float_t sump2 = 0.; + + for(int it = 0;itPx(), track->Py(), track->Pz()); + Float_t phi = track->Phi(); + Float_t eta = track->Eta(); + Float_t pt = track->Pt(); + Float_t deta = eta - etabig; + Float_t dphi = RelativePhi(phi,phibig); + if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue; + Float_t r = TMath::Sqrt(dphi * dphi + deta * deta); + if (r < 0.4 && pt>fCutTM) { + //longitudinal and perpendicular component of the track pT in the + //local frame + TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1; + TVector3 pPerp = pp - pLong; + //projection onto the two perpendicular vectors defined above + Float_t ppjX = pPerp.Dot(ppJ2); + Float_t ppjY = pPerp.Dot(ppJ3); + Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY); + //components of the 2D symmetrical sphericity matrix + mxx += (ppjX * ppjX / ppjT); + myy += (ppjY * ppjY / ppjT); + mxy += (ppjX * ppjY / ppjT); + nc++; + sump2 += ppjT;} + // max pt + if(nc<2) continue; + + } // 1st Track Loop + + + // Sphericity Matrix + const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2}; + TMatrixDSym m0(2,ele); + // Find eigenvectors + TMatrixDSymEigen m(m0); + TVectorD eval(2); + TMatrixD evecm = m.GetEigenVectors(); + eval = m.GetEigenValues(); + // Largest eigenvector + Int_t jev = 0; + if (eval[0] < eval[1]) jev = 1; + TVectorD evec0(2); + // Principle axis + evec0 = TMatrixDColumn(evecm, jev); + TVector2 evec(evec0[0], evec0[1]); + Float_t circ=0; + if(jev==1) circ=2*eval[0]; + if(jev==0) circ=2*eval[1]; + fh2Circularity->Fill(circ,ptbig); + fh1JetEntries->Fill(ptbig); + for (Int_t ip = 0; ip < nT; ip++) { + AliVParticle *track = (AliVParticle*)ParticleList.At(ip); + TVector3 pp(track->Px(), track->Py(), track->Pz()); + Float_t phi = track->Phi(); + Float_t eta = track->Eta(); + Float_t pt = track->Pt(); + + Float_t deta = eta - etabig; + Float_t dphi = RelativePhi(phi,phibig); + if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue; + + Float_t dRR = TMath::Sqrt(dphi * dphi + deta * deta); + TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1; + TVector3 pPerp = pp - pLong; + Float_t ppjX = pPerp.Dot(ppJ2); + Float_t ppjY = pPerp.Dot(ppJ3); + TVector2 vr(ppjX, ppjY) ; + //and this is the angle between the particle and the TM axis. + float phistr=evec.Phi()-vr.Phi(); + if(phistr>2*TMath::Pi()) phistr -= 2*TMath::Pi(); + if(phistr<-2*TMath::Pi()) phistr += 2*TMath::Pi(); + if(phistr<-0.5*TMath::Pi()) phistr += 2*TMath::Pi(); + if(phistr>1.5*TMath::Pi()) phistr -= 2*TMath::Pi(); + + double jetEntries[6] = {dRR,ptbig,pt,phistr,circ,nc}; + fhnJetTM->Fill(jetEntries); + + } // 2nd Track loop + }//jet loop + + PostData(1, fOutputList); +} + +void AliAnalysisTaskJetAntenna::Terminate(const Option_t *) +{ + // Draw result to the screen + // Called once at the end of the query + + if (!GetOutputData(1)) + return; +} + + +Int_t AliAnalysisTaskJetAntenna::GetListOfTracks(TList *list){ + + Int_t iCount = 0; + AliAODEvent *aod = 0; + + if(!fESD)aod = fAODIn; + else aod = fAODOut; + + if(!aod)return 0; + + Int_t index=-1; + Double_t ptmax=-10; + + + + for(int it = 0;it < aod->GetNumberOfTracks();++it){ + AliAODTrack *tr = aod->GetTrack(it); + Bool_t bGood = false; + if(fFilterType == 0)bGood = true; + else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal(); + else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal(); + if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue; + if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;} + if(bGood==false) continue; + if (fApplySharedClusterCut) { + Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls()); + if (frac > 0.4) continue; + } + if(TMath::Abs(tr->Eta())>0.9)continue; + if(tr->Pt()<0.15)continue; + list->Add(tr); + iCount++; + if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks + if(tr->TestFilterBit(fFilterMaskBestPt)){ + if(tr->Pt()>ptmax){ + ptmax=tr->Pt(); + index=iCount-1; + } + } + } + else{ + if(tr->Pt()>ptmax){ + ptmax=tr->Pt(); + index=iCount-1; + } + } + } + + return index; +} + + +Double_t AliAnalysisTaskJetAntenna::RelativePhi(Double_t mphi,Double_t vphi){ + + if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi()); + else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi()); + if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi()); + else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi()); + double dphi = mphi-vphi; + if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi()); + else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi()); + return dphi;//dphi in [-Pi, Pi] +} + +Int_t AliAnalysisTaskJetAntenna::GetPhiBin(Double_t phi) +{ + Int_t phibin=-1; + if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;} + Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi))); + phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi())); + if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");} + return phibin; +} diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h b/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h new file mode 100644 index 00000000000..bd93fb8056b --- /dev/null +++ b/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h @@ -0,0 +1,168 @@ +#ifndef ALIANALYSISTASKJETANTENNA_H +#define ALIANALYSISTASKJETANTENNA_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +// ************************************** +// This task computes several jet observables like +// the fraction of energy in inner and outer coronnas, +// the distance from track to jet axis and a +// correlation strength distribution of particles inside jets. +// Author: lcunquei@cern.ch +// ******************************************* + +class TH1F; +class TH1I; +class TH2F; +class TH3F; +class THnSparse; +class TRandom3; +class AliESDEvent; +class AliAODExtension; +class AliAODEvent; + +#include "AliAnalysisTaskSE.h" +#include "AliVEvent.h" + +class AliAnalysisTaskJetAntenna : public AliAnalysisTaskSE { +public: + AliAnalysisTaskJetAntenna(); + AliAnalysisTaskJetAntenna(const char *name); + virtual ~AliAnalysisTaskJetAntenna(); + virtual void LocalInit() {Init();} + virtual void Init(); + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(const Option_t*); + + Double_t RelativePhi(Double_t angle1,Double_t angle2); + Int_t GetPhiBin(Double_t phi); + + + virtual AliVEvent::EOfflineTriggerTypes GetOfflineTrgMask() const { return fOfflineTrgMask; } + virtual void GetBranchNames(TString &branch1, TString &branch2) const { branch1 = fJetBranchName[0]; branch2 = fJetBranchName[1]; } + virtual Bool_t GetIsPbPb() const { return fIsPbPb; } + virtual Int_t GetMinContribVtx() const { return fMinContribVtx; }; + virtual Float_t GetVtxZMin() const { return fVtxZMin; } + virtual Float_t GetVtxZMax() const { return fVtxZMax; } + virtual Int_t GetEvtClassMin() const { return fEvtClassMin; } + virtual Int_t GetEvtClassMax() const { return fEvtClassMax; } + virtual Float_t GetCentMin() const { return fCentMin; } + virtual Float_t GetCentMax() const { return fCentMax; } + virtual Int_t GetNInputTracksMin() const { return fNInputTracksMin; } + virtual Int_t GetNInputTracksMax() const { return fNInputTracksMax; } + virtual Float_t GetJetEtaMin() const { return fJetEtaMin; } + virtual Float_t GetJetEtaMax() const { return fJetEtaMax; } + virtual Float_t GetJetPtMin() const { return fJetPtMin; } + virtual Float_t GetJetPtFractionMin() const { return fJetPtFractionMin; } + virtual Int_t GetNMatchJets() const { return fNMatchJets; } + virtual void SetBranchNames(const TString &branch1, const TString &branch2); + virtual void SetBackgroundBranch(TString &branch) { fBackgroundBranch = branch;} + virtual void SetIsPbPb(Bool_t b=kTRUE) { fIsPbPb = b; } + 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 SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; } + virtual void SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; } + virtual void SetFilterMask(UInt_t i){fFilterMask = i;} + virtual void SetFilterMaskBestPt(UInt_t i){fFilterMaskBestPt = i;} + virtual void SetFilterType(Int_t iType){fFilterType=iType;} + + virtual void SetCentMin(Float_t cent) { fCentMin = cent; } + virtual void SetCentMax(Float_t cent) { fCentMax = cent; } + virtual void SetNInputTracksMin(Int_t nTr) { fNInputTracksMin = nTr; } + virtual void SetNInputTracksMax(Int_t nTr) { fNInputTracksMax = nTr; } + virtual void SetRequireITSRefit(Int_t nref) {fRequireITSRefit=nref;} + virtual void SetSharedClusterCut(Int_t docut){fApplySharedClusterCut=docut;} + + virtual void SetNRPBins(Int_t bins){fNRPBins=bins;} + virtual void SetSemigoodCorrect(Int_t yesno){fSemigoodCorrect=yesno;} + virtual void SetHolePos(Float_t poshole) { fHolePos = poshole; } + virtual void SetHoleWidth(Float_t holewidth) { fHoleWidth = holewidth; } + virtual void SetTrackTypeRec(Int_t i){fTrackTypeRec = i;} + virtual void SetTMCut(Int_t i){fCutTM=i;} + virtual void SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; } + virtual void SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; } + virtual void SetJetPtMin(Float_t pt) { fJetPtMin = pt; } + virtual void SetJetTriggerExclude(UChar_t i) { fJetTriggerExcludeMask = i; } + virtual void SetJetPtFractionMin(Float_t frac) { fJetPtFractionMin = frac; } + virtual void SetNMatchJets(Int_t n) { fNMatchJets = n; } + virtual void SetKeepJets(Bool_t b = kTRUE) { fKeepJets = b; } + virtual void SetNonStdFile(char* c){fNonStdFile = c;} + enum {kTrackUndef = 0, kTrackAOD, kTrackKineAll,kTrackKineCharged, kTrackAODMCAll, kTrackAODMCCharged, kTrackAODMCChargedAcceptance}; + +private: + // ESD/AOD events + 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 + Int_t GetListOfTracks(TList *list); + Int_t SelectTrigger(TList *list,Double_t minT,Double_t maxT,Int_t &number); + Int_t GetHardestTrackBackToJet(AliAODJet *jet); + Int_t GetListOfTracksCloseToJet(TList *list,AliAODJet *jet); + // jets to compare + TString fJetBranchName[2]; // name of jet branches to compare + TList *fListJets[2]; //! jet lists + + TString fBackgroundBranch; + TString fNonStdFile; // name of delta aod file to catch the extension + // event selection + Bool_t fIsPbPb; // is Pb-Pb (fast embedding) or p-p (detector response) + AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept + Int_t fMinContribVtx; // minimum number of track contributors for primary vertex + Float_t fVtxZMin; // lower bound on vertex z + Float_t fVtxZMax; // upper bound on vertex z + Int_t fEvtClassMin; // lower bound on event class + Int_t fEvtClassMax; // upper bound on event class + UInt_t fFilterMask; // filter bit for slecected tracks + UInt_t fFilterMaskBestPt; // filter bit for selected hig pt tracks (best quality) + UInt_t fFilterType; // type of slected tracks parrallel to filtermask + Float_t fCentMin; // lower bound on centrality + Float_t fCentMax; // upper bound on centrality + Int_t fNInputTracksMin; // lower bound of nb. of input tracks + Int_t fNInputTracksMax; // upper bound of nb. of input tracks + Int_t fRequireITSRefit; + Int_t fApplySharedClusterCut; // flag to apply shared cluster cut (needed for some AODs where this cut was not applied in the filtering) + Int_t fTrackTypeRec; + Float_t fRPAngle; + Int_t fNRPBins; + Int_t fSemigoodCorrect; + Float_t fHolePos; + Float_t fHoleWidth; + Float_t fCutTM; //lower pt cut for particles entering the TM axis calculation + Float_t fJetEtaMin; // lower bound on eta for found jets + Float_t fJetEtaMax; // upper bound on eta for found jets + Int_t fNevents; // number of events + Int_t fTindex; // index reference + Float_t fJetPtMin; // minimum jet pT + UChar_t fJetTriggerExcludeMask; // mask for jet triggeres to exclude + Float_t fJetPtFractionMin; // minimum fraction for positiv match of jets + Int_t fNMatchJets; // maximal nb. of jets taken for matching + Double_t fMatchMaxDist; // maximal distance of matching jets + Bool_t fKeepJets; // keep jets with negative pt after background subtraction + + + // output objects + const Int_t fkNbranches; //! number of branches to be read + const Int_t fkEvtClasses; //! number of event classes + + TList *fOutputList; //! output data container + + TH1I *fHistEvtSelection; //! event selection statistic + TH1F* fh1JetEntries; //centrality bias of triggers + TH2F* fh2Circularity; //jet density + + THnSparse *fhnJetTM; //Recoil jet spectrum + + + AliAnalysisTaskJetAntenna(const AliAnalysisTaskJetAntenna&); // not implemented + AliAnalysisTaskJetAntenna& operator=(const AliAnalysisTaskJetAntenna&); // not implemented + + ClassDef(AliAnalysisTaskJetAntenna, 6); +}; + +#endif + diff --git a/PWGJE/macros/AddTaskJetAntenna.C b/PWGJE/macros/AddTaskJetAntenna.C new file mode 100644 index 00000000000..f76b862fec0 --- /dev/null +++ b/PWGJE/macros/AddTaskJetAntenna.C @@ -0,0 +1,52 @@ + + +AliAnalysisTaskJetAntenna* AddTaskJetAntenna(const char* bRec1,const char* bRec2, UInt_t filterMask = 272 , Float_t ptTrackMin = 0.15, Int_t kTriggerMask=0, Float_t fCutTM=0.15){ + + Printf("adding task jet antenna\n"); + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if(!mgr){ + ::Error("AddTaskJetAntenna", "No analysis manager to connect to."); + return NULL; + } + if(!mgr->GetInputEventHandler()){ + ::Error("AddTaskJetAntenna", "This task requires an input event handler."); + return NULL; + } + + + + + TString typeRec(bRec1); + TString typeGen(bRec2); + + AliAnalysisTaskJetAntenna *task = new AliAnalysisTaskJetAntenna(Form("JetAntenna_%s_%s_%d_%f",bRec1,bRec2,kTriggerMask,fCutTM)); + + + + task->SetBranchNames(bRec1,bRec2); + task->SetOfflineTrgMask(kTriggerMask); + task->SetCentMin(0.); + task->SetCentMax(100.); + task->SetFilterMask(filterMask); + task->SetTMCut(fCutTM); + + + + + + mgr->AddTask(task); + + + AliAnalysisDataContainer *coutputJetAntenna = mgr->CreateContainer(Form("pwgjejetantenna_%s_%s_%d_%f",bRec1,bRec2,kTriggerMask,fCutTM), TList::Class(),AliAnalysisManager::kOutputContainer,Form("%s:PWGJE_jetantenna_%s_%s_%d_%f",AliAnalysisManager::GetCommonFileName(),bRec1,bRec2,kTriggerMask,fCutTM)); + + + + + + mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer()); + mgr->ConnectOutput(task, 1, coutputJetAntenna); + + return task; +} -- 2.43.0