New antenna radiation task
authorlcunquei <leticia.cunqueiro.mendez@cern.ch>
Thu, 22 May 2014 15:28:31 +0000 (17:28 +0200)
committermvl <marco.van.leeuwen@cern.ch>
Thu, 22 May 2014 15:29:23 +0000 (17:29 +0200)
PWGJE/CMakelibPWGJE.pkg
PWGJE/PWGJELinkDef.h
PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx [new file with mode: 0644]
PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h [new file with mode: 0644]
PWGJE/macros/AddTaskJetAntenna.C [new file with mode: 0644]

index 40c1f03..5446c08 100644 (file)
@@ -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)
index 2ce323a..3208106 100644 (file)
@@ -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 (file)
index 0000000..bb1c42e
--- /dev/null
@@ -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; 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;
+}
diff --git a/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h b/PWGJE/UserTasks/AliAnalysisTaskJetAntenna.h
new file mode 100644 (file)
index 0000000..bd93fb8
--- /dev/null
@@ -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 (file)
index 0000000..f76b862
--- /dev/null
@@ -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;
+}