--- /dev/null
+// **************************************
+// Task used for jet finding in the Kine Train (generation and analysis on the fly, no detector effects)
+// Output is stored in an exchance container
+// *******************************************
+
+
+/**************************************************************************
+ * 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 <TROOT.h>
+#include <TRandom3.h>
+#include <TSystem.h>
+#include <TInterpreter.h>
+#include <TChain.h>
+#include <TRefArray.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TProfile.h>
+#include <TF1.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include "TDatabasePDG.h"
+#include <TGrid.h>
+
+#include "AliAnalysisTaskJetClusterKine.h"
+#include "AliAnalysisManager.h"
+#include "AliJetFinder.h"
+#include "AliJetHeader.h"
+#include "AliJetReader.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODExtension.h"
+#include "AliAODTrack.h"
+#include "AliAODJet.h"
+#include "AliVParticle.h" //FK//
+#include "AliAODMCParticle.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliGenEventHeader.h" //FK//
+#include "AliGenPythiaEventHeader.h"
+#include "AliJetKineReaderHeader.h"
+#include "AliGenCocktailEventHeader.h"
+#include "AliInputEventHandler.h"
+#include "AliAODJetEventBackground.h"
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+// get info on how fastjet was configured
+#include "fastjet/config.h"
+
+using std::vector;
+
+ClassImp(AliAnalysisTaskJetClusterKine)
+
+AliAnalysisTaskJetClusterKine::~AliAnalysisTaskJetClusterKine(){
+ //
+ // Destructor
+ //
+
+ delete fRef;
+
+ if(fTCAJetsOut)fTCAJetsOut->Delete();
+ delete fTCAJetsOut;
+
+}
+
+//_____________________________________________________________________
+
+AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine():
+ AliAnalysisTaskSE(),
+ fMcEvent(0x0),
+ fMcHandler(0x0),
+ fRef(new TRefArray),
+ fTrackTypeGen(kTrackKineCharged), //Kine charged?
+ fAvgTrials(1),
+ fTrackEtaWindow(0.9),
+ fTrackPtCut(0.),
+ fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
+ fVtxZCut(10.0),
+ fNonStdBranch(""),
+ fOutContainer(kNoOutput), //FF//
+ fNonStdFile(""),
+ fRparam(1.0),
+ fAlgorithm(fastjet::kt_algorithm),
+ fStrategy(fastjet::Best),
+ fRecombScheme(fastjet::BIpt_scheme),
+ fAreaType(fastjet::active_area),
+ fGhostArea(0.01),
+ fActiveAreaRepeats(1),
+ fGhostEtamax(1.5),
+ fTCAJetsOut(0x0),
+ fh1Xsec(0x0),
+ fh1Trials(0x0),
+ fh1PtHard(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
+ fh1NJetsGen(0x0),
+ fh1NConstGen(0x0),
+ fh1NConstLeadingGen(0x0),
+ fh1PtJetsGenIn(0x0),
+ fh1PtJetsLeadingGenIn(0x0),
+ fh1PtJetConstGen(0x0),
+ fh1PtJetConstLeadingGen(0x0),
+ fh1PtTracksGenIn(0x0),
+ fh1Nch(0x0),
+ fh1Z(0x0),
+ fh2NConstPt(0x0),
+ fh2NConstLeadingPt(0x0),
+ fh2JetPhiEta(0x0),
+ fh2LeadingJetPhiEta(0x0),
+ fh2JetEtaPt(0x0),
+ fh2LeadingJetEtaPt(0x0),
+ fh2TrackEtaPt(0x0),
+ fh2JetsLeadingPhiEta(0x0),
+ fh2JetsLeadingPhiPt(0x0),
+ fh2JetsLeadingPhiPtW(0x0),
+ fHistList(0x0)
+{
+ //
+ // Constructor
+ //
+}
+
+//_____________________________________________________________________
+
+AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine(const char* name):
+ AliAnalysisTaskSE(name),
+ fMcEvent(0x0),
+ fMcHandler(0x0),
+ fRef(new TRefArray),
+ fTrackTypeGen(kTrackKineCharged), //kine charged?
+ fAvgTrials(1),
+ fTrackEtaWindow(0.9),
+ fTrackPtCut(0.),
+ fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
+ fVtxZCut(10.0),
+ fNonStdBranch(""),
+ fOutContainer(kNoOutput),//FF//
+ fNonStdFile(""),
+ fRparam(1.0),
+ fAlgorithm(fastjet::kt_algorithm),
+ fStrategy(fastjet::Best),
+ fRecombScheme(fastjet::BIpt_scheme),
+ fAreaType(fastjet::active_area),
+ fGhostArea(0.01),
+ fActiveAreaRepeats(1),
+ fGhostEtamax(1.5),
+ fTCAJetsOut(0x0),
+ fh1Xsec(0x0),
+ fh1Trials(0x0),
+ fh1PtHard(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
+ fh1NJetsGen(0x0),
+ fh1NConstGen(0x0),
+ fh1NConstLeadingGen(0x0),
+ fh1PtJetsGenIn(0x0),
+ fh1PtJetsLeadingGenIn(0x0),
+ fh1PtJetConstGen(0x0),
+ fh1PtJetConstLeadingGen(0x0),
+ fh1PtTracksGenIn(0x0),
+ fh1Nch(0x0),
+ fh1Z(0x0),
+ fh2NConstPt(0x0),
+ fh2NConstLeadingPt(0x0),
+ fh2JetPhiEta(0x0),
+ fh2LeadingJetPhiEta(0x0),
+ fh2JetEtaPt(0x0),
+ fh2LeadingJetEtaPt(0x0),
+ fh2TrackEtaPt(0x0),
+ fh2JetsLeadingPhiEta(0x0),
+ fh2JetsLeadingPhiPt(0x0),
+ fh2JetsLeadingPhiPtW(0x0),
+ fHistList(0x0)
+{
+ //
+ // named ctor
+ //
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TClonesArray::Class());
+}
+
+
+//_____________________________________________________________________
+
+Bool_t AliAnalysisTaskJetClusterKine::Notify(){
+ //
+ // Implemented Notify() to read the cross sections
+ // and number of trials from pyxsec.root
+ //
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+
+void AliAnalysisTaskJetClusterKine::UserCreateOutputObjects(){
+
+ //
+ // Create the output container
+ //
+
+ // Connect the AOD
+
+
+ if(fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
+
+
+ if(fNonStdBranch.Length()!=0){
+ // only create the output branch if we have a name
+ // Create a new branch for jets...
+ // -> cleared in the UserExec....
+ // here we can also have the case that the brnaches are written to a separate file
+ fTCAJetsOut = new TClonesArray("AliAODJet", 0);
+ fTCAJetsOut->SetName(fNonStdBranch.Data());
+ if(fOutContainer==kAODBranch){ //FF//
+ AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
+ }
+
+
+ //if(fNonStdFile.Length()!=0){
+ //
+ // case that we have an AOD extension we need to fetch the jets from the extended output
+ // we identify the extension aod event by looking for the branchname
+ //AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+ // case that we have an AOD extension we need can fetch the background maybe from the extended output
+ //fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
+ //}
+ }
+
+
+ if(!fHistList) fHistList = new TList();
+ fHistList->SetOwner(kTRUE);
+ PostData(1, fHistList); // post data in any case once
+
+
+ if(fOutContainer==kExchCont){
+ fTCAJetsOut->SetOwner();
+ PostData(2, fTCAJetsOut); //FF// post data in any case once
+ }
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+
+ //
+ // Histogram
+
+ const Int_t nBinPt = 100;
+ Double_t binLimitsPt[nBinPt+1];
+ for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
+ if(iPt == 0){
+ binLimitsPt[iPt] = 0.0;
+ }else {// 1.0
+ binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
+ }
+ }
+
+ const Int_t nBinPhi = 90;
+ Double_t binLimitsPhi[nBinPhi+1];
+ for(Int_t iPhi = 0; iPhi<=nBinPhi; iPhi++){
+ if(iPhi==0){
+ binLimitsPhi[iPhi] = -1.*TMath::Pi();
+ }else{
+ binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
+ }
+ }
+
+ const Int_t nBinEta = 40;
+ Double_t binLimitsEta[nBinEta+1];
+ for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
+ if(iEta==0){
+ binLimitsEta[iEta] = -2.0;
+ }else{
+ binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
+ }
+ }
+
+ const int nChMax = 5000;
+
+ fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+ fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+
+ fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
+ fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+
+
+ fh1NJetsGen = new TH1F("fh1NJetsGen","N reconstructed jets",120,-0.5,119.5);
+
+ fh1NConstGen = new TH1F("fh1NConstGen","# jet constituents",120,-0.5,119.5);
+ fh1NConstLeadingGen = new TH1F("fh1NConstLeadingGen","jet constituents",120,-0.5,119.5);
+
+
+ fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
+ fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
+ fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+
+ fh1PtJetsGenIn = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
+ fh1PtJetsLeadingGenIn = new TH1F("fh1PtJetsLeadingGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
+ fh1PtJetConstGen = new TH1F("fh1PtJetsConstGen","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
+ fh1PtJetConstLeadingGen = new TH1F("fh1PtJetsConstLeadingGen","Gen jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
+ fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("Gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
+ fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
+
+
+ fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
+
+
+
+ fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
+ fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
+
+
+ fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
+ nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
+ fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
+ nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
+
+ fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
+ nBinEta,binLimitsEta,nBinPt,binLimitsPt);
+ fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
+ nBinEta,binLimitsEta,nBinPt,binLimitsPt);
+
+ fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
+ nBinEta,binLimitsEta,nBinPt,binLimitsPt);
+
+
+
+ fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
+ nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
+ fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
+ nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+ fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
+ nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+
+
+
+ const Int_t saveLevel = 3; // large save level more histos
+ if(saveLevel>0){
+ fHistList->Add(fh1Xsec);
+ fHistList->Add(fh1Trials);
+
+ fHistList->Add(fh1NJetsGen);
+ fHistList->Add(fh1NConstGen);
+ fHistList->Add(fh1NConstLeadingGen);
+ fHistList->Add(fh1PtJetsGenIn);
+ fHistList->Add(fh1PtJetsLeadingGenIn);
+ fHistList->Add(fh1PtTracksGenIn);
+ fHistList->Add(fh1PtJetConstGen);
+ fHistList->Add(fh1PtJetConstLeadingGen);
+ fHistList->Add(fh1Nch);
+ fHistList->Add(fh1Z);
+
+ fHistList->Add(fh2NConstPt);
+ fHistList->Add(fh2NConstLeadingPt);
+ fHistList->Add(fh2JetPhiEta);
+ fHistList->Add(fh2LeadingJetPhiEta);
+ fHistList->Add(fh2JetEtaPt);
+ fHistList->Add(fh2LeadingJetEtaPt);
+ fHistList->Add(fh2TrackEtaPt);
+ fHistList->Add(fh2JetsLeadingPhiEta);
+ fHistList->Add(fh2JetsLeadingPhiPt);
+ fHistList->Add(fh2JetsLeadingPhiPtW);
+ }
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for(Int_t i=0; i<fHistList->GetEntries(); ++i){
+ TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
+ if(h1){
+ h1->Sumw2();
+ continue;
+ }
+ THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
+ if(hn) hn->Sumw2();
+ }
+ TH1::AddDirectory(oldStatus);
+}
+//_________________________________________________________________________
+
+void AliAnalysisTaskJetClusterKine::Init(){
+ // MC handler
+ if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Init() \n");
+ fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); //FK//
+}
+
+//_________________________________________________________________________
+void AliAnalysisTaskJetClusterKine::LocalInit(){
+ // MC handler
+ //
+ // Initialization
+ //
+
+ if(fDebug > 1) printf("AnalysisTaskJetClusterKine::LocalInit() \n");
+ Init();
+}
+
+//_________________________________________________________________________
+void AliAnalysisTaskJetClusterKine::UserExec(Option_t* /*option*/){
+ // handle and reset the output jet branch
+
+ if(fDebug > 1) printf("AliAnalysisTaskJetClusterKine::UserExec \n");
+ if(fTCAJetsOut) fTCAJetsOut->Delete(); //clean TClonesArray
+
+ //
+ // Execute analysis for current event
+ //
+ Init();
+ if(fMcHandler){
+ fMcEvent = fMcHandler->MCEvent();
+ }else{
+ if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Handler() fMcHandler=NULL\n");
+ PostData(1, fHistList);
+
+ if(fOutContainer==kExchCont){
+ PostData(2, fTCAJetsOut); //FF//
+ }
+
+ return;
+ }
+ if(!fMcEvent){
+ if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Exec() fMcEvent=NULL \n");
+ PostData(1, fHistList);
+
+ if(fOutContainer==kExchCont){
+ PostData(2, fTCAJetsOut); //FF//
+ }
+
+ return;
+ }
+
+ const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
+ Float_t zVtx = vtxMC->GetZ();
+ if(TMath::Abs(zVtx) > fVtxZCut){ //vertex cut
+ PostData(1, fHistList);
+
+ if(fOutContainer==kExchCont){
+ PostData(2, fTCAJetsOut); //FF//
+ }
+
+ return;
+ }
+
+ fh1Z->Fill(zVtx);
+
+ fh1Trials->Fill("#sum{ntrials}",1);
+ if(fDebug > 10) Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+ // we simply fetch the mc particles as a list of AliVParticles
+
+ TList genParticles; //list of particles with above min pT and good eta range
+
+ Int_t nParticles = GetListOfTracks(&genParticles, fTrackTypeGen); //fill the list
+ fh1Nch->Fill(nParticles);
+
+ if(fDebug>2) Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__, nParticles, genParticles.GetEntries());
+
+ // find the jets....
+
+ vector<fastjet::PseudoJet> inputParticlesKine;
+
+ for(int i = 0; i < genParticles.GetEntries(); i++){ //loop over generated particles
+ AliVParticle *vp = (AliVParticle*) genParticles.At(i);
+ fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
+ jInp.set_user_index(i);
+ inputParticlesKine.push_back(jInp);
+
+ if(fTCAJetsOut){
+ if(i == 0){
+ fRef->Delete(); // make sure to delete before placement new...
+ new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
+ }
+ fRef->Add(vp); //TRefArray does not work with toy model ...
+ }
+ }//generated particles
+
+ if(inputParticlesKine.size()==0){ //FK//
+ if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
+ PostData(1, fHistList);
+
+ if(fOutContainer==kExchCont){
+ PostData(2, fTCAJetsOut); //FF//
+ }
+
+ return;
+ } //FK//
+
+ // run fast jet
+ // employ setters for these...
+ // now create the object that holds info about ghosts
+
+ fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
+ fastjet::AreaType areaType = fastjet::active_area;
+ fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
+ fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
+ fastjet::ClusterSequenceArea clustSeq(inputParticlesKine, jetDef, areaDef); //FK//
+
+ //range where to compute background
+ Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+ phiMin = 0;
+ phiMax = 2*TMath::Pi();
+ rapMax = fGhostEtamax - fRparam;
+ rapMin = - fGhostEtamax + fRparam;
+ fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
+
+
+ const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
+ const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
+
+
+ fh1NJetsGen->Fill(sortedJets.size()); //the number of found jets
+
+ // loop over all jets an fill information, first on is the leading jet
+
+ Int_t nGen = inclusiveJets.size();
+ if(inclusiveJets.size()>0){
+
+ //leading Jet
+ AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
+ Double_t area = clustSeq.area(sortedJets[0]);
+ leadingJet.SetEffArea(area,0);
+ Float_t pt = leadingJet.Pt();
+ Float_t phi = leadingJet.Phi();
+ if(phi<0) phi += TMath::Pi()*2.;
+ Float_t eta = leadingJet.Eta();
+
+ //inclusive jet
+ Int_t nAodOutJets = 0;
+ AliAODJet *aodOutJet = NULL;
+
+
+ // correlation of leading jet with tracks
+ //FK// TIterator *particleIter = genParticles.MakeIterator();//FK//
+ //FK// particleIter->Reset();
+ //FK// AliVParticle *tmpParticle = NULL;
+ //FK// while((tmpParticle = (AliVParticle*)(particleIter->Next()))){
+ //FK// Float_t tmpPt = tmpParticle->Pt();
+ //FK// // correlation
+ //FK// Float_t tmpPhi = tmpParticle->Phi();
+ //FK// if(tmpPhi<0) tmpPhi+= TMath::Pi()*2.;
+ //FK// Float_t dPhi = phi - tmpPhi;
+ //FK// if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
+ //FK// if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); //-pi,pi
+ //FK// fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
+ //FK// fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
+ //FK// }
+
+ TLorentzVector vecareab;
+ for(int j = 0; j < nGen;j++){ //loop over inclusive jets
+ AliAODJet tmpGenJet (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
+ aodOutJet = NULL;
+ Float_t tmpPt = tmpGenJet.Pt(); //incl jet pt
+
+ if((tmpPt > fJetOutputMinPt) && fTCAJetsOut){// cut on the non-background subtracted...
+ aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpGenJet);
+ aodOutJet->GetRefTracks()->Clear();
+ Double_t area1 = clustSeq.area(sortedJets[j]);
+ aodOutJet->SetEffArea(area1,0);
+ fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
+ vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
+ aodOutJet->SetVectorAreaCharged(&vecareab);
+ }
+
+ fh1PtJetsGenIn->Fill(tmpPt); //incl jet Pt
+ // Fill Spectra with constituentsemacs
+ const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
+
+ fh1NConstGen->Fill(constituents.size()); //number of constituents
+ fh2NConstPt->Fill(tmpPt,constituents.size()); //number of constituents vs jet pt
+
+ // loop over constiutents and fill spectrum
+ AliVParticle *partLead = 0;
+ Float_t ptLead = -1;
+
+ for(unsigned int ic = 0; ic < constituents.size(); ic++){
+ AliVParticle *part = (AliVParticle*) genParticles.At(constituents[ic].user_index());
+ if(!part) continue;
+ fh1PtJetConstGen->Fill(part->Pt()); //pt of constituents
+
+ if(aodOutJet){
+ aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));//FK//
+
+ if(part->Pt()>fMaxTrackPtInJet){
+ aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+ }
+ }
+
+ if(part->Pt()>ptLead){
+ ptLead = part->Pt();
+ partLead = part;
+ }
+
+ if(j==0) fh1PtJetConstLeadingGen->Fill(part->Pt()); //pt of leading jet constituents
+ }
+
+ if(partLead){
+ if(aodOutJet){
+ //set pT of leading constituent of jet
+ aodOutJet->SetPtLeading(partLead->Pt());
+ }
+ }
+
+ // correlation
+ Float_t tmpPhi = tmpGenJet.Phi(); //incl jet phi
+ Float_t tmpEta = tmpGenJet.Eta(); //incl jet eta
+
+ if(tmpPhi<0) tmpPhi += TMath::Pi()*2.;
+ fh2JetPhiEta->Fill(tmpGenJet.Phi(),tmpEta);
+ fh2JetEtaPt->Fill(tmpEta,tmpPt);
+
+ if(j==0){ //leading jet
+ fh1PtJetsLeadingGenIn->Fill(tmpPt); //leading jet pt
+ fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
+ fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
+ fh1NConstLeadingGen->Fill(constituents.size()); //number of constituents in leading jet
+ fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
+ continue;
+ }
+
+ Float_t dPhi = phi - tmpPhi;
+ if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
+ if(dPhi<(-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi(); //-pi.pi
+ Float_t dEta = eta - tmpGenJet.Eta();
+ fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
+ fh2JetsLeadingPhiPt->Fill(dPhi,pt);
+ fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
+ }// loop over reconstructed jets
+ //delete particleIter;
+ }//number of jets>0
+
+
+ // fill track information
+ Int_t nTrackOver = genParticles.GetSize();
+ // do the same for tracks and jets
+
+ if(nTrackOver>0){
+ TIterator *particleIter = genParticles.MakeIterator();
+ AliVParticle *tmpGen = 0;
+ particleIter->Reset();
+
+ while((tmpGen = (AliVParticle*)(particleIter->Next()))){
+ Float_t tmpPt = tmpGen->Pt();
+ Float_t tmpEta = tmpGen->Eta();
+ fh1PtTracksGenIn->Fill(tmpPt);
+ fh2TrackEtaPt->Fill(tmpEta,tmpPt);
+ }
+ delete particleIter;
+ }
+
+
+
+
+
+ if(fDebug > 2){
+ if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
+ }
+ PostData(1, fHistList);
+
+ if(fOutContainer==kExchCont){
+ PostData(2, fTCAJetsOut); //FF//
+ }
+
+
+}
+//____________________________________________________________________________
+
+void AliAnalysisTaskJetClusterKine::Terminate(Option_t* /*option*/){
+ //
+ // Terminate analysis
+ //
+ if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
+
+
+}
+
+//_____________________________________________________________________________________
+
+Int_t AliAnalysisTaskJetClusterKine::GetListOfTracks(TList *list, Int_t type){
+
+ //
+ // get list of tracks/particles for different types
+ //
+
+ if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
+
+ Int_t iCount = 0; //number of particles
+
+ if(type == kTrackKineAll || type == kTrackKineCharged){
+ if(! fMcEvent) return iCount;
+
+ //we want to have alivpartilces so use get track
+ for(int it = 0; it < fMcEvent->GetNumberOfTracks(); ++it){
+ if(!fMcEvent->IsPhysicalPrimary(it)) continue;
+
+ AliMCParticle* part = (AliMCParticle*)fMcEvent->GetTrack(it);
+ if(TMath::Abs(part->Eta()) > fTrackEtaWindow) continue; //apply eta cut
+ if(part->Pt() < fTrackPtCut) continue; //apply pT cut
+
+ if(type == kTrackKineAll){ // full jets
+ list->Add(part);
+ iCount++;
+
+ }else if(type == kTrackKineCharged){ //charged jets
+ if(part->Particle()->GetPDG()->Charge()==0) continue;
+ list->Add(part);
+ iCount++;
+ }
+ }
+ }
+
+ list->Sort(); //?//
+
+ return iCount; //number of particles in the list
+}
+
+
--- /dev/null
+#ifndef ALIANALYSISTASKJETCLUSTERKINE_H
+#define ALIANALYSISTASKJETCLUSTERKINE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// *******************************************
+// Jet Finding in the Kine train
+// *******************************************
+
+#include "AliAnalysisTaskSE.h"
+#include "THnSparse.h" // cannot forward declare ThnSparseF
+#include "AliVParticle.h" //FK//
+#ifndef __CINT__
+# include "fastjet/ClusterSequenceArea.hh"
+# include "fastjet/AreaDefinition.hh"
+# include "fastjet/JetDefinition.hh"
+#else
+namespace fastjet {
+ enum JetAlgorithm;
+ enum Strategy;
+ enum RecombinationScheme;
+ enum AreaType;
+}
+#endif
+
+////////////////
+class AliJetHeader;
+class AliESDEvent;
+class AliAODEvent;
+class AliAODExtension;
+class AliAODJet;
+class AliGenPythiaEventHeader;
+class AliJetFinder;
+class AliAODMCParticle;
+class AliMCEvent; //FK//
+class AliMCEventHandler; //FK//
+class AliVParticle; //FK//
+class AliGenEventHeader; //FK//
+class TList;
+class TChain;
+class TH2F;
+class TH1F;
+class TH3F;
+class TRefArray;
+class TClonesArray;
+class TF1;
+class TProfile;
+
+class AliAnalysisTaskJetClusterKine : public AliAnalysisTaskSE
+{
+ public:
+ AliAnalysisTaskJetClusterKine();
+ AliAnalysisTaskJetClusterKine(const char* name);
+ virtual ~AliAnalysisTaskJetClusterKine();
+ // Implementation of interface methods
+ virtual void UserCreateOutputObjects();
+ virtual void Init(); //FK//
+ virtual void LocalInit();
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *option);
+ virtual Bool_t Notify();
+
+
+ virtual void SetTrackEtaWindow(Float_t f){fTrackEtaWindow = f;}
+ virtual void SetTrackTypeGen(Int_t i){fTrackTypeGen = i;}
+ virtual void SetTrackPtCut(Float_t x){fTrackPtCut = x;}
+ virtual void SetVtxCuts(Float_t z){fVtxZCut = z;}
+
+ virtual void SetJetOutputBranch(const char *c){fNonStdBranch = c;}
+ virtual void SetJetOutputContainer(Int_t c){fOutContainer = c;} //FF//
+ virtual const char* GetJetOutputBranch(){return fNonStdBranch.Data();}
+ virtual void SetJetOutputFile(const char *c){fNonStdFile = c;}
+ virtual const char* GetJetOutputFile(){return fNonStdFile.Data();}
+ virtual void SetMaxTrackPtInJet(Float_t x){fMaxTrackPtInJet = x;}
+ virtual void SetJetOutputMinPt(Float_t x){fJetOutputMinPt = x;}
+
+
+
+ // for Fast Jet
+ fastjet::JetAlgorithm GetAlgorithm() const {return fAlgorithm;}
+ fastjet::Strategy GetStrategy() const {return fStrategy;}
+ fastjet::RecombinationScheme GetRecombScheme() const {return fRecombScheme;}
+ fastjet::AreaType GetAreaType() const {return fAreaType;}
+ // Setters
+ void SetRparam(Double_t f) {fRparam = f;}
+ // Temporary change to integer; problem with dictionary generation?
+ //void SetAlgorithm(fastjet::JetAlgorithm f) {fAlgorithm = f;}
+ void SetAlgorithm(Int_t f) {fAlgorithm = (fastjet::JetAlgorithm) f;}
+ void SetStrategy(fastjet::Strategy f) {fStrategy = f;}
+ void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;}
+ void SetAreaType(fastjet::AreaType f) {fAreaType = f;}
+ void SetGhostArea(Double_t f) {fGhostArea = f;}
+ void SetActiveAreaRepeats(Int_t f) {fActiveAreaRepeats = f;}
+ void SetGhostEtamax(Double_t f) {fGhostEtamax = f;}
+
+
+
+ // Helper
+ //
+
+ // we have different cases
+ // AOD reading -> MC from AOD
+ // ESD reading -> MC from Kinematics
+ // this has to match with our selection of input events
+ enum {kTrackKineAll=0, kTrackKineCharged};
+ enum { kNoOutput=0, kAODBranch=1, kExchCont=2 }; //FF//
+
+
+ private:
+
+ AliAnalysisTaskJetClusterKine(const AliAnalysisTaskJetClusterKine&);
+ AliAnalysisTaskJetClusterKine& operator=(const AliAnalysisTaskJetClusterKine&);
+
+ Int_t GetListOfTracks(TList *list,Int_t type);
+
+ AliMCEvent* fMcEvent; //! MC event //FK// FFF
+ AliInputEventHandler* fMcHandler; //! MCEventHandler //FK// FFF
+ TRefArray *fRef; // ! trefarray for track references within the jet FFF
+ Int_t fTrackTypeGen; // type of tracks used for 0=full, 1=charged FFF
+ Float_t fAvgTrials; // Average nimber of trials
+ Float_t fTrackEtaWindow; // eta window used for corraltion plots between rec and gen FFF
+ Float_t fTrackPtCut; // minimum track pt to be accepted FFF
+ Float_t fJetOutputMinPt; // minimum p_t for jets to be written out FFF
+ Float_t fMaxTrackPtInJet; // maximum track pt within a jet for flagging... FFF
+ Float_t fVtxZCut; // zvtx cut
+
+ // output configurartion
+ TString fNonStdBranch; // the name of the non-std branch name, if empty no branch is filled FFF
+ Int_t fOutContainer; //FF//output container 1=AOD Branch 2=Exchange container
+ TString fNonStdFile; // The optional name of the output file the non-std branch is written to FFF
+
+
+
+ // Fast jet
+ Double_t fRparam; // fastjet distance parameter FFF
+ fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm FFF
+ fastjet::Strategy fStrategy; //= fastjet::Best; FFF
+ fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme; FFF
+ fastjet::AreaType fAreaType; // fastjet area type
+ Double_t fGhostArea; // fasjet ghost area FFF
+ Int_t fActiveAreaRepeats; // fast jet active area repeats FFF
+ Double_t fGhostEtamax; // fast jet ghost area FFF
+
+ TClonesArray *fTCAJetsOut; //! TCA of output jets FFF
+
+ TProfile* fh1Xsec; //! pythia cross section and trials
+ TH1F* fh1Trials; //! trials are added FFF
+ TH1F* fh1PtHard; //! Pt har of the event...
+ TH1F* fh1PtHardNoW; //! Pt har of the event without weigt
+ TH1F* fh1PtHardTrials; //! Number of trials
+
+ TH1F* fh1NJetsGen; //! number of generator level jets FFF
+ TH1F* fh1NConstGen;//! number of constiutens in leading jet FFF
+ TH1F* fh1NConstLeadingGen;//! number of constiutens in leading jet FFF
+ TH1F* fh1PtJetsGenIn; //! Jet pt for all jets FFF
+ TH1F* fh1PtJetsLeadingGenIn; //! Jet pt for the leading jets FFF
+ TH1F* fh1PtJetConstGen;//! pt of constituents FFF
+ TH1F* fh1PtJetConstLeadingGen;// pt of constituents FFF
+ TH1F* fh1PtTracksGenIn; //! track pt for all tracks FFF
+
+ TH1F* fh1Nch; //! (charged) particle mult FFF
+ TH1F* fh1Z; // ! centrality of anaylsed events FFF
+
+
+ TH2F* fh2NConstPt; //! number of constituents vs. pt FFF
+ TH2F* fh2NConstLeadingPt; //! number of constituents vs. pt FFF
+ TH2F* fh2JetPhiEta; //! jet phi eta FFF
+ TH2F* fh2LeadingJetPhiEta; //! leading jet phi eta FFF
+ TH2F* fh2JetEtaPt; //! leading jet eta FFF
+ TH2F* fh2LeadingJetEtaPt; //! leading jet eta FFF
+ TH2F* fh2TrackEtaPt; //! track eta for all tracks FFF
+ TH2F* fh2JetsLeadingPhiEta; //! jet delta phi delta eta w.r.t. leading jet FFF
+ TH2F* fh2JetsLeadingPhiPt; //! jet correlation with leading jet FFF
+ TH2F* fh2JetsLeadingPhiPtW; //! jet correlation with leading jet FFF
+
+ TList *fHistList; //!leading tracks to be skipped in the randomized event Output list
+
+
+ ClassDef(AliAnalysisTaskJetClusterKine, 1)
+};
+
+#endif