Add cluster task for kine train
authormvl <marco.van.leeuwen@cern.ch>
Mon, 7 Apr 2014 15:44:06 +0000 (17:44 +0200)
committermvl <marco.van.leeuwen@cern.ch>
Mon, 7 Apr 2014 15:44:58 +0000 (17:44 +0200)
JETAN/AliAnalysisTaskJetClusterKine.cxx [new file with mode: 0644]
JETAN/AliAnalysisTaskJetClusterKine.h [new file with mode: 0644]
JETAN/CMakelibFASTJETAN.pkg
JETAN/FASTJETANLinkDef.h
PWGJE/macros/AddTaskJetClusterKine.C [new file with mode: 0644]

diff --git a/JETAN/AliAnalysisTaskJetClusterKine.cxx b/JETAN/AliAnalysisTaskJetClusterKine.cxx
new file mode 100644 (file)
index 0000000..3f184c7
--- /dev/null
@@ -0,0 +1,739 @@
+// **************************************
+// 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
+}
+
+
diff --git a/JETAN/AliAnalysisTaskJetClusterKine.h b/JETAN/AliAnalysisTaskJetClusterKine.h
new file mode 100644 (file)
index 0000000..f4e86e2
--- /dev/null
@@ -0,0 +1,183 @@
+#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
index d62975b..457a574 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS AliFastJetFinder.cxx AliFastJetHeaderV1.cxx AliFastJetInput.cxx AliFastJetBkg.cxx AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx AliAnalysisTaskJetBackgroundSubtract.cxx)
+set ( SRCS AliFastJetFinder.cxx AliFastJetHeaderV1.cxx AliFastJetInput.cxx AliFastJetBkg.cxx AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx AliAnalysisTaskJetClusterKine.cxx AliAnalysisTaskJetBackgroundSubtract.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index 96a29fc..c1525e7 100644 (file)
@@ -11,5 +11,6 @@
 #pragma        link C++ class AliSISConeJetFinder+;
 #pragma        link C++ class AliSISConeJetHeader+;
 #pragma        link C++ class AliAnalysisTaskJetCluster+;
+#pragma        link C++ class AliAnalysisTaskJetClusterKine+;
 #pragma        link C++ class AliAnalysisTaskJetBackgroundSubtract+;
 #endif
diff --git a/PWGJE/macros/AddTaskJetClusterKine.C b/PWGJE/macros/AddTaskJetClusterKine.C
new file mode 100644 (file)
index 0000000..1848fc4
--- /dev/null
@@ -0,0 +1,116 @@
+AliAnalysisTaskJetClusterKine *AddTaskJetClusterKine(char* bGen = "KINECHARGED",Char_t *jf = "ANTIKT", Float_t radius = 0.4, Int_t kWriteAOD = 1, char* deltaFile = "", Float_t ptTrackCut = 0.15, Float_t etaTrackWindow = 0.9, Float_t vertexWindow = 10.);
+
+Float_t kPtTrackCutCl     = 0.15;
+Float_t kTrackEtaWindowCl = 0.8;
+Float_t kVertexWindowCl   = 10.0;
+
+
+AliAnalysisTaskJetClusterKine *AddTaskJetClusterKine(char* bGen, Char_t *jf, Float_t radius, Int_t kWriteAOD, char *deltaFile, Float_t ptTrackCut, Float_t etaTrackWindow, Float_t vertexWindow){
+
+ // Creates a jet fider task, configures it and adds it to the analysis manager.
+   kPtTrackCutCl     = ptTrackCut;
+   kTrackEtaWindowCl = etaTrackWindow;
+   kVertexWindowCl   = vertexWindow;
+
+   TString outputFile(deltaFile);
+    // Get the pointer to the existing analysis manager via the static access method.
+    //==============================================================================
+    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+    if(!mgr){
+       ::Error("AddTaskJetClusterKine", "No analysis manager to connect to.");
+       return NULL;
+    }  
+
+    // Check the analysis type using the event handlers connected to the analysis manager.
+    //==============================================================================
+    if(!mgr->GetMCtruthEventHandler()){
+      ::Error("AddTaskJetClusterKine", "This task requires an input MC event handler");
+       return NULL;
+    }
+
+    TString typeGen(bGen);
+    typeGen.ToUpper();
+
+    // Create the task and configure it.
+    //===========================================================================
+
+    TString cAdd = "";
+    cAdd += Form("%02d_",TMath::Nint(radius*10.));
+    cAdd += Form("Cut%05d",TMath::Nint(1000.*kPtTrackCutCl));
+    
+    Printf("%s %s%s", typeGen.Data(), jf, cAdd.Data());
+
+    AliAnalysisTaskJetClusterKine* clus = new  AliAnalysisTaskJetClusterKine(Form("JetCluster%s_%s%s",bGen,jf,cAdd.Data()));
+      
+    // or a config file
+    clus->SetVtxCuts(kVertexWindowCl);
+
+    if(typeGen.Contains("KINECHARGED")){
+       clus->SetTrackTypeGen(AliAnalysisTaskJetClusterKine::kTrackKineCharged);
+       clus->SetTrackPtCut(kPtTrackCutCl);
+       clus->SetTrackEtaWindow(kTrackEtaWindowCl);
+
+    }else if(typeGen.Contains("KINEFULL")){
+       clus->SetTrackTypeGen(AliAnalysisTaskJetClusterKine::kTrackKineAll);
+       clus->SetTrackPtCut(kPtTrackCutCl);
+       clus->SetTrackEtaWindow(kTrackEtaWindowCl);
+    }
+
+   clus->SetRparam(radius);
+   clus->SetGhostArea(0.005);
+   clus->SetGhostEtamax(kTrackEtaWindowCl);
+
+   clus->SetDebugLevel(0);
+
+   switch (jf) {
+   case "ANTIKT":
+     clus->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh
+     break;
+   case "CA":
+     clus->SetAlgorithm(1); // CA from fastjet/JetDefinition.hh
+     break;
+   case "KT":
+     clus->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh
+     break;
+   default:
+     ::Error("AddTaskJetClusterKine", "Wrong jet finder selected\n");
+     return 0;
+   }
+
+   TString nameOutArray =  Form("clusters%s_%s%s",bGen,jf,cAdd.Data()); //FF//
+   if(kWriteAOD){
+     if(outputFile.Length())clus->SetJetOutputFile(outputFile);
+     Printf("Output branch: %s",nameOutArray.Data());//FF//  
+     clus->SetJetOutputBranch(nameOutArray.Data());//FF//
+     clus->SetJetOutputMinPt(0); // store only jets / clusters above a certain threshold
+   }
+   clus->SetJetOutputContainer(kWriteAOD); //0=no output 1=AOD 2=Exchange
+
+   mgr->AddTask(clus);
+
+   // Create ONLY the output containers for the data produced by the task.
+   // Get and connect other common input/output containers via the manager as below
+   //==============================================================================
+   AliAnalysisDataContainer *coutput1_clus = mgr->CreateContainer(Form("cluster_%s_%s%s",bGen,jf,cAdd.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,Form("%s:PWGJE_cluster_%s_%s%s",AliAnalysisManager::GetCommonFileName(),bGen,jf,cAdd.Data()));
+
+   mgr->ConnectInput  (clus, 0, mgr->GetCommonInputContainer());
+
+   if(kWriteAOD==1){//FF//
+      mgr->ConnectOutput (clus, 0, mgr->GetCommonOutputContainer());
+
+   }
+
+   mgr->ConnectOutput (clus, 1, coutput1_clus );
+
+
+
+   if(kWriteAOD==2){//FF//
+      AliAnalysisDataContainer *coutput2_clus = mgr->CreateContainer( nameOutArray.Data(), //??
+                                           TClonesArray::Class(), 
+                                           AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectOutput (clus, 2, coutput2_clus); //FF//
+   }
+
+
+   return clus;
+}