--- /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 "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; 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 AliAnalysisTaskJetAntenna::UserExec(Option_t *)
+{
+
+
+ if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].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());
+
+ 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<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
+
+
+ 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<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
+ if(!(inputHandler->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<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
+ aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));
+ }
+ if(fAODExtension && !aodJets[0]){
+ aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
+ aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
+ }
+ if(fAODIn&&!aodJets[0]){
+ aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
+ aodJets[1] = dynamic_cast<TClonesArray*>(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<AliAODJet*>((*aodJets[iJetType])[iJet]);
+ if (jet) fListJets[iJetType]->Add(jet);
+ }
+ }
+
+
+
+ for(Int_t i=0; i<fListJets[0]->GetEntries(); ++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((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
+
+ if(fSemigoodCorrect){
+ if((phibig>fHolePos-fHoleWidth) && (phibig<fHolePos+fHoleWidth)) continue;
+ }
+
+
+ //two vectors perpendicular to the jet axis
+ pxbig=jetbig->Px();
+ 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;it<nT;++it){
+ AliVParticle *track = (AliVParticle*)ParticleList.At(it);
+ 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 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;
+}
--- /dev/null
+#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
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff