Add Spectrum Task and helper class, some minor mods
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Aug 2008 09:18:30 +0000 (09:18 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Aug 2008 09:18:30 +0000 (09:18 +0000)
PWG4/AliAnalysisHelperJetTasks.cxx [new file with mode: 0644]
PWG4/AliAnalysisHelperJetTasks.h [new file with mode: 0644]
PWG4/AliAnalysisTaskJetSpectrum.cxx [new file with mode: 0644]
PWG4/AliAnalysisTaskJetSpectrum.h [new file with mode: 0644]
PWG4/AliAnalysisTaskUE.h
PWG4/Makefile
PWG4/PROOF-INF.PWG4JetTasks/SETUP.C
PWG4/PWG4JetTasksLinkDef.h
PWG4/libPWG4JetTasks.pkg
PWG4/macros/AnalysisTrainCAF.C

diff --git a/PWG4/AliAnalysisHelperJetTasks.cxx b/PWG4/AliAnalysisHelperJetTasks.cxx
new file mode 100644 (file)
index 0000000..d0278a0
--- /dev/null
@@ -0,0 +1,43 @@
+
+#include "TROOT.h"
+#include "TList.h"
+#include "AliMCEvent.h"
+#include "AliGenEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+#include "AliGenPythiaEventHeader.h"
+#include <fstream>
+#include <iostream>
+#include "AliAnalysisHelperJetTasks.h"
+
+
+ClassImp(AliAnalysisHelperJetTasks)
+
+
+
+AliGenPythiaEventHeader*  AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){
+  
+  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
+  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
+  if(!pythiaGenHeader){
+    // cocktail ??
+    AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
+    
+    if (!genCocktailHeader) {
+      Printf("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__);
+      return 0;
+    }
+    TList* headerList = genCocktailHeader->GetHeaders();
+    for (Int_t i=0; i<headerList->GetEntries(); i++) {
+      pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
+      if (pythiaGenHeader)
+        break;
+    }
+    if(!pythiaGenHeader){
+      Printf("%s %d: PythiaHeader not found!",(char*)__FILE__,__LINE__);
+      return 0;
+    }
+  }
+  return pythiaGenHeader;
+
+}
diff --git a/PWG4/AliAnalysisHelperJetTasks.h b/PWG4/AliAnalysisHelperJetTasks.h
new file mode 100644 (file)
index 0000000..7f437e7
--- /dev/null
@@ -0,0 +1,24 @@
+#ifndef  ALIANALYSISHELPERJETTASKS_H
+#define  ALIANALYSISHELPERJETTASKS_H
+
+
+#include "TObject.h"
+class AliMCEvent;
+class AliGenPythiaEventHeader;
+
+// Helper Class that contains a lot of usefull static functions (i.e. for Flavor selection.
+
+class AliAnalysisHelperJetTasks : public TObject {
+ public:
+  AliAnalysisHelperJetTasks() : TObject() {;}
+  virtual ~AliAnalysisHelperJetTasks(){;}
+  
+  static AliGenPythiaEventHeader*  GetPythiaEventHeader(AliMCEvent *mcEvent);
+  
+
+  private:
+  
+  ClassDef(AliAnalysisHelperJetTasks, 1) // 
+};
+
+#endif // ALIANALYSISHELPERJETTASKS_H
diff --git a/PWG4/AliAnalysisTaskJetSpectrum.cxx b/PWG4/AliAnalysisTaskJetSpectrum.cxx
new file mode 100644 (file)
index 0000000..4ae3dc0
--- /dev/null
@@ -0,0 +1,588 @@
+/**************************************************************************
+ * 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 <TSystem.h>
+#include <TInterpreter.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include  "TDatabasePDG.h"
+
+#include "AliAnalysisTaskJetSpectrum.h"
+#include "AliAnalysisManager.h"
+#include "AliJetFinder.h"
+#include "AliJetReader.h"
+#include "AliJetReaderHeader.h"
+#include "AliUA1JetHeaderV1.h"
+#include "AliJet.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODTrack.h"
+#include "AliAODJet.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliJetKineReaderHeader.h"
+#include "AliGenCocktailEventHeader.h"
+
+#include "AliAnalysisHelperJetTasks.h"
+
+ClassImp(AliAnalysisTaskJetSpectrum)
+
+AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
+  fJetFinderRec(0x0),
+  fJetFinderGen(0x0),
+  fAOD(0x0),
+  fBranchRec("jets"),
+  fConfigRec("ConfigJets.C"),
+  fBranchGen(""),
+  fConfigGen(""),
+  fUseAODInput(kFALSE),
+  fUseExternalWeightOnly(kFALSE),
+  fAnalysisType(0),
+  fExternalWeight(1),    
+  fh1PtHard(0x0),
+  fh1PtHard_NoW(0x0),  
+  fh1PtHard_Trials(0x0),
+  fh1PtHard_Trials_NoW(0x0),
+  fh1NGenJets(0x0),
+  fh1NRecJets(0x0),
+  fHistList(0x0)
+{
+  // Default constructor
+    for(int ij  = 0;ij<kMaxJets;++ij){
+    fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
+    fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
+    fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3RecEtaPhiPt_NoFound[ij] =  fh3MCEtaPhiPt[ij] = 0;
+  }
+  
+}
+
+AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
+  AliAnalysisTaskSE(name),
+  fJetFinderRec(0x0),
+  fJetFinderGen(0x0),
+  fAOD(0x0),
+  fBranchRec("jets"),
+  fConfigRec("ConfigJets.C"),
+  fBranchGen(""),
+  fConfigGen(""),
+  fUseAODInput(kFALSE),
+  fUseExternalWeightOnly(kFALSE),
+  fAnalysisType(0),
+  fExternalWeight(1),    
+  fh1PtHard(0x0),
+  fh1PtHard_NoW(0x0),  
+  fh1PtHard_Trials(0x0),
+  fh1PtHard_Trials_NoW(0x0),
+  fh1NGenJets(0x0),
+  fh1NRecJets(0x0),
+  fHistList(0x0)
+{
+  // Default constructor
+  for(int ij  = 0;ij<kMaxJets;++ij){
+    fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
+    fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
+
+    fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3RecEtaPhiPt_NoFound[ij] =  fh3MCEtaPhiPt[ij] = 0;
+  }
+
+  DefineOutput(1, TList::Class());  
+}
+
+
+
+
+void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
+{
+
+  //
+  // Create the output container
+  //
+
+  
+  // Connect the AOD
+
+  if(fUseAODInput){
+    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
+      return;
+    }    
+  }
+  else{
+    //  assume that the AOD is in the general output...
+    fAOD  = AODEvent();
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
+      return;
+    }    
+  }
+
+
+  if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
+
+  OpenFile(1);
+  if(!fHistList)fHistList = new TList();
+
+  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;
+    }
+  }
+  const Int_t nBinEta = 22;
+  Double_t binLimitsEta[nBinEta+1];
+  for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
+    if(iEta==0){
+      binLimitsEta[iEta] = -1.1;
+    }
+    else{
+      binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
+    }
+  }
+
+  const Int_t nBinPhi = 360;
+  Double_t binLimitsPhi[nBinPhi+1];
+  for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
+    if(iPhi==0){
+      binLimitsPhi[iPhi] = 0;
+    }
+    else{
+      binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
+    }
+  }
+
+  const Int_t nBinFrag = 25;
+
+
+  fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
+
+  fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
+
+  fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+
+  fh1PtHard_Trials_NoW = new TH1F("fh1PtHard_Trials_NoW","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+
+  fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
+
+  fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
+
+
+  for(int ij  = 0;ij<kMaxJets;++ij){
+    fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
+    fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
+    fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
+    fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
+    fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
+
+
+    fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
+                            nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+
+
+
+
+
+
+    fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+
+
+
+    fh3PtRecGenHard_NoW[ij] = new TH3F(Form("fh3PtRecGenHard_NoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+
+       
+    fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
+                          nBinFrag,0.,1.,nBinPt,binLimitsPt);
+
+    fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
+                            nBinFrag,0.,10.,nBinPt,binLimitsPt);
+
+    fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+                                 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+
+
+    fh3RecEtaPhiPt_NoGen[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+                                       nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+
+    fh3RecEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoFound_g%d",ij),"No found for generated jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+                                       nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+
+
+    fh3MCEtaPhiPt[ij] = new TH3F(Form("fh3MCEtaPhiPt_j%d",ij),"MC eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+                                nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+
+  }
+
+  const Int_t saveLevel = 1; // large save level more histos
+
+  if(saveLevel>0){
+    fHistList->Add(fh1PtHard);
+    fHistList->Add(fh1PtHard_NoW);
+    fHistList->Add(fh1PtHard_Trials);
+    fHistList->Add(fh1PtHard_Trials_NoW);
+    fHistList->Add(fh1NGenJets);
+    fHistList->Add(fh1NRecJets);
+    for(int ij  = 0;ij<kMaxJets;++ij){
+      fHistList->Add(fh1E[ij]);
+      fHistList->Add(fh1PtRecIn[ij]);
+      fHistList->Add(fh1PtRecOut[ij]);
+      fHistList->Add(fh1PtGenIn[ij]);
+      fHistList->Add(fh1PtGenOut[ij]);
+      fHistList->Add(fh2PtFGen[ij]);
+      if(saveLevel>2){
+       fHistList->Add(fh3RecEtaPhiPt[ij]);
+       fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
+       fHistList->Add(fh3RecEtaPhiPt_NoFound[ij]);
+       fHistList->Add(fh3MCEtaPhiPt[ij]);      
+      }
+    }
+  }
+
+  // =========== 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){
+      // Printf("%s ",h1->GetName());
+      h1->Sumw2();
+    }
+  }
+
+  TH1::AddDirectory(oldStatus);
+
+}
+
+void AliAnalysisTaskJetSpectrum::Init()
+{
+  //
+  // Initialization
+  //
+
+  Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
+  if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
+
+}
+
+void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
+{
+  //
+  // Execute analysis for current event
+  //
+
+  if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
+
+  
+  AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+
+  if(!aodH){
+    Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
+    return;
+  }
+
+
+  //  aodH->SelectEvent(kTRUE);
+
+  // ========= These pointers need to be valid in any case ======= 
+
+
+  /*
+  AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
+  if(!jhRec){
+    Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
+    return;
+  }
+  */
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
+  if(!aodRecJets){
+    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
+    return;
+  }
+
+  // ==== General variables needed
+
+
+  // We use statice array, not to fragment the memory
+  AliAODJet genJets[kMaxJets];
+  Int_t nGenJets = 0;
+  AliAODJet recJets[kMaxJets];
+  Int_t nRecJets = 0;
+
+  Double_t eventW = 1;
+  Double_t ptHard = 0; 
+  Double_t nTrials = 1; // Trials for MC trigger weigth for real data
+
+  if(fUseExternalWeightOnly){
+    eventW = fExternalWeight;
+  }
+
+
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  if((fAnalysisType&kAnaMC)==kAnaMC){
+    // this is the part we only use when we have MC information
+    AliMCEvent* mcEvent = MCEvent();
+    //    AliStack *pStack = 0; 
+    if(!mcEvent){
+      Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
+      return;
+    }
+    AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
+    if(!pythiaGenHeader){
+      return;
+    }
+
+    nTrials = pythiaGenHeader->Trials();
+    ptHard  = pythiaGenHeader->GetPtHard();
+
+    if(!fUseExternalWeightOnly){
+       // case were we combine more than one p_T hard bin...
+       // eventW = AnalysisHelperCKB::GetXSectionWeight(pythiaGenHeader->GetPtHard(),fEnergy)*fExternalWeight;      
+    }
+
+    // fetch the pythia generated jets only to be used here
+    Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
+    AliAODJet pythiaGenJets[kMaxJets];
+
+    if(fBranchGen.Length()==0)nGenJets = nPythiaGenJets;
+    for(int ip = 0;ip < nPythiaGenJets;++ip){
+      if(ip>=kMaxJets)continue;
+      Float_t p[4];
+      pythiaGenHeader->TriggerJet(ip,p);
+      pythiaGenJets[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
+      if(fBranchGen.Length()==0){
+       // if we have MC particles and we do not read from the aod branch
+       // use the pythia jets
+       genJets[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
+      }
+    }
+
+
+  }// (fAnalysisType&kMC)==kMC)
+
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  fh1PtHard->Fill(ptHard,eventW);
+  fh1PtHard_NoW->Fill(ptHard,1);
+  fh1PtHard_Trials->Fill(ptHard,nTrials);
+
+
+  // If we set a second branch for the input jets fetch this 
+  if(fBranchGen.Length()>0){
+    TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
+    if(aodGenJets){
+      nGenJets = aodGenJets->GetEntries();
+      for(int ig = 0;ig < nGenJets;++ig){
+       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
+       if(!tmp)continue;
+       genJets[ig] = *tmp;
+      }
+    }
+    else{
+      Printf("%s:%d Generated jet branch %s not found",fBranchGen.Data());
+    }
+  }
+
+  fh1NGenJets->Fill(nGenJets);
+  // We do not want to exceed the maximum jet number
+  nGenJets = TMath::Min(nGenJets,kMaxJets);
+
+  // Fetch the reconstructed jets...
+
+
+  nRecJets = aodRecJets->GetEntries();
+  for(int ir = 0;ir < nRecJets;++ir){
+    AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
+    if(!tmp)continue;
+    recJets[ir] = *tmp;
+  }
+
+  fh1NRecJets->Fill(nRecJets);
+  nRecJets = TMath::Min(nRecJets,kMaxJets);
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  // Relate the jets
+  Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
+  Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
+  
+  for(int i = 0;i<kMaxJets;++i){
+    iGenIndex[i] = iRecIndex[i] = -1;
+  }
+
+
+  GetClosestJets(genJets,nGenJets,recJets,nRecJets,
+                iGenIndex,iRecIndex,fDebug);
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+  // loop over reconstructed jets
+  for(int ir = 0;ir < nRecJets;++ir){
+    Double_t ptRec = recJets[ir].Pt();
+    Double_t phiRec = recJets[ir].Phi();
+    if(phiRec<0)phiRec+=TMath::Pi()*2.;    
+    Double_t etaRec = recJets[ir].Eta();
+
+    fh1E[ir]->Fill(recJets[ir].E(),eventW);
+    fh1PtRecIn[ir]->Fill(ptRec,eventW);
+    fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
+    // Fill Correlation
+    Int_t ig = iGenIndex[ir];
+    if(ig>=0&&ig<nGenJets){
+      fh1PtRecOut[ir]->Fill(ptRec,eventW);
+      Double_t ptGen = genJets[ig].Pt();
+      fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
+      fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
+      fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
+    }
+  }// loop over reconstructed jets
+
+
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  for(int ig = 0;ig < nGenJets;++ig){
+    Double_t ptGen = genJets[ig].Pt();
+    // Fill Correlation
+    fh1PtGenIn[ig]->Fill(ptGen,eventW);
+    Int_t ir = iRecIndex[ig];
+    if(ir>=0){
+      fh1PtGenOut[ig]->Fill(ptGen,eventW);
+    }
+  }// loop over reconstructed jets
+
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  PostData(1, fHistList);
+}
+
+void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
+{
+// Terminate analysis
+//
+    if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
+}
+
+
+void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
+                                               AliAODJet *recJets,Int_t &nRecJets,
+                                               Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
+
+  //
+  // Relate the two input jet Arrays
+  //
+
+  //
+  // The association has to be unique
+  // So check in two directions
+  // find the closest rec to a gen
+  // and check if there is no other rec which is closer
+  // Caveat: Close low energy/split jets may disturb this correlation
+
+  // Idea: search in two directions generated e.g (a--e) and rec (1--3)
+  // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
+  // in the end we have something like this
+  //    1   2   3
+  // ------------
+  // a| 3   2   0
+  // b| 0   1   0
+  // c| 0   0   3
+  // d| 0   0   1
+  // e| 0   0   1
+  // Topology
+  //   1     2
+  //     a         b        
+  //
+  //  d      c
+  //        3     e
+  // Only entries with "3" match from both sides
+
+  Int_t iFlag[kMaxJets][kMaxJets];
+
+  for(int i = 0;i < kMaxJets;++i){
+    iRecIndex[i] = -1;
+    iGenIndex[i] = -1;
+    for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
+  }
+
+  if(nRecJets==0)return;
+  if(nGenJets==0)return;
+
+  const Float_t maxDist = 1.4;
+  // find the closest distance
+  for(int ig = 0;ig<nGenJets;++ig){
+    Float_t dist = maxDist;
+    if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
+    for(int ir = 0;ir<nRecJets;++ir){
+      Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
+      if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
+      if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
+      if(dR<dist){
+       iRecIndex[ig] = ir;
+       dist = dR;
+      }        
+    }
+    if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
+  }
+  // other way around
+  for(int ir = 0;ir<nRecJets;++ir){
+    Float_t dist = maxDist;
+    for(int ig = 0;ig<nGenJets;++ig){
+      Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
+      if(dR<dist){
+       iGenIndex[ir] = ig;
+       dist = dR;
+      }        
+    }
+    if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
+  }
+
+  // check for "true" correlations
+
+  if(iDebug>1)Printf(">>>>>> Matrix");
+
+  for(int ig = 0;ig<nGenJets;++ig){
+    for(int ir = 0;ir<nRecJets;++ir){
+      // Print
+      if(iDebug>1)printf("%d ",iFlag[ig][ir]);
+      if(iFlag[ig][ir]==3){
+       iGenIndex[ir] = ig;
+       iRecIndex[ig] = ir;
+      }
+      else{
+       iGenIndex[ir] = iRecIndex[ig] = -1;
+      }
+    }
+    if(iDebug>1)printf("\n");
+  }
+}
+
+
diff --git a/PWG4/AliAnalysisTaskJetSpectrum.h b/PWG4/AliAnalysisTaskJetSpectrum.h
new file mode 100644 (file)
index 0000000..90ef36a
--- /dev/null
@@ -0,0 +1,109 @@
+#ifndef ALIANALYSISTASKJETSPECTRUM_H
+#define ALIANALYSISTASKJETSPECTRUM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+#include "AliAnalysisTaskSE.h"
+class AliJetFinder;
+class AliESDEvent;
+class AliAODEvent;
+class AliAODJet;
+class AliGenPythiaEventHeader;
+
+class TList;
+class TChain;
+class TH2F;
+class TH3F;
+
+class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
+{
+ public:
+    AliAnalysisTaskJetSpectrum();
+    AliAnalysisTaskJetSpectrum(const char* name);
+    virtual ~AliAnalysisTaskJetSpectrum() {;}
+    // Implementation of interface methods
+    virtual void UserCreateOutputObjects();
+    virtual void Init();
+    virtual void LocalInit() { Init(); }
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *option);
+
+    virtual void SetExternalWeight(Float_t f){fExternalWeight = f;}
+    virtual void SetUseExternalWeightOnly(Bool_t b){fUseExternalWeightOnly = b;}
+    virtual void SetAODInput(Bool_t b){fUseAODInput = b;}
+    virtual void SetAnalysisType(Int_t b){fAnalysisType = b;}
+    virtual void SetBranchGen(char* c){fBranchGen = c;}
+    virtual void SetBranchRec(char* c){fBranchRec = c;}
+
+    // Helper
+    static void GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
+                              AliAODJet *recJets,Int_t &nRecJets,
+                              Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug = 0);
+
+  //
+
+    enum {kAnaMC =  0x1};
+
+ private:
+
+    AliAnalysisTaskJetSpectrum(const AliAnalysisTaskJetSpectrum&);
+    AliAnalysisTaskJetSpectrum& operator=(const AliAnalysisTaskJetSpectrum&);
+    
+
+    enum {kMaxJets =  5};
+
+    AliJetFinder *fJetFinderRec;
+    AliJetFinder *fJetFinderGen;
+    AliAODEvent  *fAOD; // where we take the jets from can be input or output AOD
+
+    TString       fBranchRec;  // AOD branch name for reconstructed
+    TString       fConfigRec;  // Name of the Config file 
+    TString       fBranchGen;  // AOD brnach for genereated
+    TString       fConfigGen;  // Name of the Config file (if any)
+
+    Bool_t        fUseAODInput;
+    Bool_t        fUseExternalWeightOnly;
+    Int_t         fAnalysisType;
+    Float_t       fExternalWeight;    
+
+    TH1F*         fh1PtHard;  // Pt har of the event...       
+    TH1F*         fh1PtHard_NoW;  // Pt har of the event...       
+    TH1F*         fh1PtHard_Trials;  // Number of trials 
+    TH1F*         fh1PtHard_Trials_NoW;  // Number of trials 
+    TH1F*         fh1NGenJets;
+    TH1F*         fh1NRecJets;
+    TH1F*         fh1E[kMaxJets];       // Jet Energy       
+    TH1F*         fh1PtRecIn[kMaxJets];       // Jet pt for all      
+    TH1F*         fh1PtRecOut[kMaxJets];      // Jet pt with corellated generated jet    
+    TH1F*         fh1PtGenIn[kMaxJets];       // Detection efficiency for given p_T.gen
+    TH1F*         fh1PtGenOut[kMaxJets];      //
+
+
+    
+    TH2F*         fh2PtFGen[kMaxJets];  //
+    TH2F*         fh2Frag[kMaxJets];    // fragmentation function
+    TH2F*         fh2FragLn[kMaxJets];  //
+
+    TH3F*         fh3PtRecGenHard[kMaxJets];  //                              
+    TH3F*         fh3PtRecGenHard_NoW[kMaxJets];  //                  
+    TH3F*         fh3RecEtaPhiPt[kMaxJets]; // 
+    TH3F*         fh3RecEtaPhiPt_NoGen[kMaxJets]; // 
+    TH3F*         fh3RecEtaPhiPt_NoFound[kMaxJets]; //                    
+    TH3F*         fh3MCEtaPhiPt[kMaxJets]; //                                                                                                        
+    // ========= Multiplicity dependence ======
+
+    // ==========TODO , flavaor dependence ========                           
+    // ============================================                       
+                                     
+
+    // ============= TODO , phi dependence ========
+    // ============================================                            
+  
+    TList *fHistList; // Output list
+
+
+    ClassDef(AliAnalysisTaskJetSpectrum, 1) // Analysis task for standard jet analysis
+};
+#endif
index 974fab08acea1f13ea9d6cfbe3df7301420e183d..0e3ad6274c3de7742d8e3768f95046d6c472178e 100644 (file)
@@ -132,7 +132,7 @@ AliAnalysisTaskUE&   operator=(const  AliAnalysisTaskUE &det);
               
       //        TH2F*  fhValidRegion; //! test to be canceled
          
-    ClassDef( AliAnalysisTaskUE, 1); // Analysis task for Underlying Event analysis
+    ClassDef( AliAnalysisTaskUE, 1) // Analysis task for Underlying Event analysis
 };
  
 #endif
index ccd2932d0d4a94a8a67328db42202fb8c7fc78bb..68561d58c47c6c9988104294c14f709eb700e880 100644 (file)
@@ -30,6 +30,10 @@ ifneq ($(ANALYSISalice_INCLUDE),)
   ALICEINC += -I../$(ANALYSISalice_INCLUDE)
 endif
 
+ifneq ($(JETAN_INCLUDE),)
+  ALICEINC += -I../$(JETAN_INCLUDE)
+endif
+
 ifneq ($(PWG4PartCorr_INCLUDE),)
   ALICEINC += -I../$(PWG4PartCorr_INCLUDE)
 endif
@@ -44,6 +48,7 @@ endif
 ifeq ($(ALICEINC),-I.)
   ifneq ($(ALICE_ROOT),)
     ALICEINC += -I$(ALICE_ROOT)/include
+    ALICEINC += -I$(ALICE_ROOT)/JETAN   # some extra includes
   endif
 endif
 
index 1d79478be0137380978ca634502721faa110aa53..2c73d5e06c6848388a9a232b6679d073df9449f8 100644 (file)
@@ -4,8 +4,12 @@ void SETUP()
    // Load the JET-Tasks library
    gSystem->Load("libPWG4JetTasks");
 
+
+   gSystem->AddIncludePath("-I$ROOTSYS/include -IJETAN");
+   gROOT->ProcessLine(".include JETAN");
+
    // Set the Include paths
-   gSystem->AddIncludePath("-I$ROOTSYS/include -IPWG4");
+   gSystem->AddIncludePath("-IPWG4");
    gROOT->ProcessLine(".include PWG4");
 
    // Set our location, so that other packages can find us
index 80b9442efd3ec59adcf8d9c51e022c6dcea76fd9..7b698ee871bd12e01e7ba3cda65ae1b61bf1a473 100644 (file)
@@ -5,6 +5,8 @@
 #pragma link off all functions;
 
 #pragma link C++ class AliAnalysisTaskUE+;
+#pragma link C++ class AliAnalysisTaskJetSpectrum+;
+#pragma link C++ class AliAnalysisHelperJetTasks+;
 
 
 #endif
index b41ddad5d074472a4d6423d6602997cb68fdc246..cb216c3b4cc7eb4f3d05b0dd9cf43e3ca656734b 100644 (file)
@@ -1,6 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS = AliAnalysisTaskUE.cxx 
+SRCS = AliAnalysisTaskUE.cxx AliAnalysisTaskJetSpectrum.cxx AliAnalysisHelperJetTasks.cxx
 
 HDRS:= $(SRCS:.cxx=.h) 
 
@@ -8,8 +8,10 @@ DHDR:= PWG4JetTasksLinkDef.h
 
 EXPORT:=$(SRCS:.cxx=.h)
 
+EINCLUDE=JETAN
+
 ifeq (win32gcc,$(ALICE_TARGET))
 PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
-                         -lESD -lAOD -lANALYSIS -lANALYSISalice \
+                         -lESD -lAOD -lANALYSIS -lANALYSISalice -lJETAN \
                          -L$(shell root-config --libdir) -lEG
 endif
index 27f62bdf7c8588230ea6cd20256c25700ced2efc..5cfb1ad38d31d999b9afac34df8eb005c2ce27b7 100644 (file)
@@ -1,23 +1,29 @@
 //______________________________________________________________________________
-void AnalysisTrainCAF()
+void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0)
 {
 // Example of running analysis train in CAF. To run in debug mode:
 //  - export ROOTSYS=debug  on your local client
 //  - un-comment gProof->ClearPackages()
-//  - un-comme lines with debugging info
+//  - un-comment lines with debugging info
 
 // WHY AOD is not a exchange container when running from ESD->AOD
 
     Bool_t debug         = kTRUE;
     Bool_t useMC         = kTRUE;
     Bool_t readTR        = kFALSE;
+    Bool_t bPROOF         = kTRUE;
+
+
     
     Int_t iAODanalysis   = 0;
     Int_t iAODhandler    = 1;
     Int_t iESDfilter     = 1;  // Only active if iAODanalysis=0
     Int_t iJETAN         = 1;
-    Int_t iJETANMC       = 1;
-    Int_t iPWG4UE      = 1;
+    Int_t iJETANAOD      = 0;
+    Int_t iJETANMC       = 0;
+    Int_t iDIJETAN       = 0;
+    Int_t iPWG4SPECTRUM  = 1;
+    Int_t iPWG4UE        = 0;
 
     if (iAODanalysis) {
        useMC = kFALSE;
@@ -25,11 +31,16 @@ void AnalysisTrainCAF()
        iESDfilter = 0;
        iMUONfilter = 0;
     }    
-    if (iJETAN) iESDfilter=1;
+    if (iJETANAOD) iESDfilter=1;
     if (iESDfilter) iAODhandler=1;
 
     // Dataset from CAF
     TString dataset = "/PWG4/arian/jetjet15-50";
+    // CKB quick hack
+    gROOT->LoadMacro("CreateESDChain.C");
+    TChain *chain = CreateESDChain("jetjet15-50.txt",10);
+
+
     printf("==================================================================\n");
     printf("===========    RUNNING ANALYSIS TRAIN IN CAF MODE    =============\n");
     printf("==================================================================\n");
@@ -37,8 +48,11 @@ void AnalysisTrainCAF()
     else              printf("=  ESD analysis on dataset: %s\n", dataset.Data());
     if (iESDfilter)   printf("=  ESD filter                                                    =\n");
     if (iJETAN)       printf("=  Jet analysis                                                  =\n");
-    if (iJETANMC)       printf("=  Jet analysis from Kinematics                                  =\n");
-    if (iPWG4UE)  printf("=  PWG4 UE                                                   =\n");
+    if (iJETANAOD)    printf("=  Jet analysis from AOD                                        =\n");
+    if (iJETANMC)     printf("=  Jet analysis from Kinematics                                  =\n");
+    if (iDIJETAN)     printf("=  DiJet analysis                                                =\n");
+    if (iPWG4SPECTRUM)printf("=  PWG4 Jet spectrum analysis                                    =\n");
+    if (iPWG4UE)      printf("=  PWG4 UE                                                       =\n");
     printf("==================================================================\n");
     if (useMC) printf(":: use MC    TRUE\n");
     else       printf(":: use MC    FALSE\n");
@@ -61,50 +75,65 @@ void AnalysisTrainCAF()
     const char* proofNode = "lxb6046";
     //    const char* proofNode = "kleinb@localhost";
 
-    TProof::Mgr(proofNode)->ShowROOTVersions();
-    //    TProof::Mgr(proofNode)->SetROOTVersion("vHEAD-r24503_dbg");
+
+
 
     // Connect to proof
-    TProof::Open(proofNode); // may be username@lxb6046 if user not the same as on local
-
-    // Clear packages if changing ROOT version on CAF or local
-    //    gProof->ClearPackages();
-    // Enable proof debugging if needed
-    //    gProof->SetLogLevel(5);
-    // To debug the train in PROOF mode, type in a root session:
-    // root[0] TProof::Mgr("lxb6064")->GetSessionLogs()->Display("*",0,10000);
-    // Common packages
-    // --- Enable the STEERBase Package
-    gProof->UploadPackage("pars/STEERBase.par");
-    gProof->EnablePackage("STEERBase");
-    // --- Enable the ESD Package
-    gProof->UploadPackage("pars/ESD.par");
-    gProof->EnablePackage("ESD");
-    // --- Enable the AOD Package
-    gProof->UploadPackage("pars/AOD.par");
-    gProof->EnablePackage("AOD");
-    // --- Enable the ANALYSIS Package
-    gProof->UploadPackage("pars/ANALYSIS.par");
-    gProof->EnablePackage("ANALYSIS");
-    // --- Enable the ANALYSISalice Package
-    gProof->UploadPackage("pars/ANALYSISalice.par");
-    gProof->EnablePackage("ANALYSISalice");
-
-    AliPDG::AddParticlesToPdgDataBase();
-
-    // --- Enable the JETAN Package
-    if (iJETAN||iJETANMC) {
-       gProof->UploadPackage("pars/JETAN.par");
-       gProof->EnablePackage("JETAN");
-    }   
-    // --- Enable particle correlation analysis
-    if (iPWG4UE) {
-       gProof->UploadPackage("pars/PWG4JetTasks.par");
-       gProof->EnablePackage("PWG4JetTasks");
-    }   
+    if(bPROOF){
+      TProof::Mgr(proofNode)->ShowROOTVersions();
+      TProof::Mgr(proofNode)->SetROOTVersion("v5-21-01-alice_dbg");
+      TProof::Open(proofNode); 
+
+      // Clear packages if changing ROOT version on CAF or local
+      //    gProof->ClearPackages();
+      // Enable proof debugging if needed
+      //    gProof->SetLogLevel(5);
+      // To debug the train in PROOF mode, type in a root session:
+      // root[0] TProof::Mgr("lxb6064")->GetSessionLogs()->Display("*",0,10000);
+      // Common packages
+      // --- Enable the STEERBase Package
+      gProof->UploadPackage("${ALICE_ROOT}/STEERBase.par");
+      gProof->EnablePackage("STEERBase");
+      // --- Enable the ESD Package
+      gProof->UploadPackage("${ALICE_ROOT}/ESD.par");
+      gProof->EnablePackage("ESD");
+      // --- Enable the AOD Package
+      gProof->UploadPackage("${ALICE_ROOT}/AOD.par");
+      gProof->EnablePackage("AOD");
+      // --- Enable the ANALYSIS Package
+      gProof->UploadPackage("${ALICE_ROOT}/ANALYSIS.par");
+      gProof->EnablePackage("ANALYSIS");
+      // --- Enable the ANALYSISalice Package
+      gProof->UploadPackage("${ALICE_ROOT}/ANALYSISalice.par");
+      gProof->EnablePackage("ANALYSISalice");
+
+
+      // --- Enable the JETAN Package
+      if (iJETAN||iJETANMC) {
+       gProof->UploadPackage("${ALICE_ROOT}/JETAN.par");
+       gProof->EnablePackage("JETAN");
+      }   
+      // --- Enable particle correlation analysis
+      if (iPWG4UE||iPWG4SPECTRUM) {
+       gProof->UploadPackage("${ALICE_ROOT}/PWG4JetTasks.par");
+       gProof->EnablePackage("PWG4JetTasks");
+      }   
+
+    }
+    else{
+      gSystem->Load("libSTEERBase");
+      gSystem->Load("libESD");
+      gSystem->Load("libAOD");
+      gSystem->Load("libANALYSIS");
+      gSystem->Load("libANALYSISalice");
+
+
+      // --- Enable the JETAN Package
+      if (iJETAN||iJETANMC) gSystem->Load("libJETAN");
+      // --- Enable particle correlation analysis
+      if (iPWG4UE||iPWG4SPECTRUM)gSystem->Load("libPWG4JetTasks");  
+    }
 
-    // Create the chain
-    //
 
     // Make the analysis manager
     AliAnalysisManager *mgr  = new AliAnalysisManager("Analysis Train", "A test setup for the analysis train");
@@ -134,17 +163,18 @@ void AnalysisTrainCAF()
        // AOD output handler
        AliAODHandler* aodHandler   = new AliAODHandler();
        mgr->SetOutputEventHandler(aodHandler);
-       aodHandler->SetOutputFileName("AliAODs.root");
-//       aodHandler->SetCreateNonStandardAOD();
+       aodHandler->SetOutputFileName(Form("AliAODs_CKB_%07d-%07d.root",nOffset,nOffset+nEvents));
+       //       aodHandler->SetSelectAll(kFALSE);
+       //       aodHandler->SetCreateNonStandardAOD();
        cout_aod = mgr->CreateContainer("cAOD", TTree::Class(),
-                                                             AliAnalysisManager::kOutputContainer, "default");
+                                      AliAnalysisManager::kOutputContainer, "default");
        cout_aod->SetSpecialOutput();
     }   
 
     // Debugging if needed
-    if (debug) mgr->SetDebugLevel(3);
+    if (debug) mgr->SetDebugLevel(0);
 //    AliLog::EnableDebug(kTRUE);
-//    AliLog::SetGlobalLogLevel(2);
+    AliLog::SetGlobalLogLevel(0);
 
 
     if (iESDfilter && !iAODanalysis) {
@@ -191,6 +221,25 @@ void AnalysisTrainCAF()
        // mgr->ConnectOutput (jetana,     0, cout_aod );
        mgr->ConnectOutput (jetana,     1, cout_jet );
     }   
+    // JETANALYSIS from the AOD
+    if (iJETANAOD) {
+       AliAnalysisTaskJets *jetanaAOD = new AliAnalysisTaskJets("AODJetAnalysis");
+       jetanaAOD->SetDebugLevel(10);
+       jetanaAOD->SetNonStdBranch("jetsAOD");    
+       jetanaAOD->SetConfigFile("${ALICE_ROOT}/JETAN/ConfigJetAnalysisAOD.C");
+       mgr->AddTask(jetanaAOD);
+       // Output histograms list for jet analysis                       
+       AliAnalysisDataContainer *cout_jetAOD = mgr->CreateContainer("jethistAOD", TList::Class(),
+                                                             AliAnalysisManager::kOutputContainer, "jethistAOD.root");
+       // Dummy AOD output container for jet analysis (no client yet)
+       c_aodjet0 = mgr->CreateContainer("cAODjet0", TTree::Class(),
+                           AliAnalysisManager::kExchangeContainer);
+       // Connect to data containers
+       mgr->ConnectInput  (jetanaAOD,     0, cout_aod  );
+       mgr->ConnectOutput (jetanaAOD,     0, c_aodjet0 );
+       // mgr->ConnectOutput (jetana,     0, cout_aod );
+       mgr->ConnectOutput (jetanaAOD,     1, cout_jetAOD );
+    }   
     // Jet analysisMC
     AliAnalysisDataContainer *c_aodjetMC = 0;
     if (iJETANMC && useMC) {
@@ -211,10 +260,51 @@ void AnalysisTrainCAF()
        // mgr->ConnectOutput (jetanaMC,     0, cout_aod );
        mgr->ConnectOutput (jetanaMC,     1, cout_jetMC );
     }   
+    // Dijet analysis
+    if(iDIJETAN){
+      AliAnalysisTaskDiJets *dijetana = new AliAnalysisTaskDiJets("DiJetAnalysis");
+      dijetana->SetDebugLevel(2);
+      
+      mgr->AddTask(dijetana);
+
+      //
+      // Create containers for input/output
+      AliAnalysisDataContainer *c_aod_dijet = mgr->CreateContainer("cAODdijet", TTree::Class(),
+                                                                  AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectInput  (dijetana,  0, cinput  );
+      mgr->ConnectOutput (dijetana,  0, c_aod_dijet);
+    }
+
+
+    if (iPWG4SPECTRUM) {
+      AliAnalysisTaskJetSpectrum* pwg4spec = new  AliAnalysisTaskJetSpectrum("Jet Spectrum");
+      
+      // default parameters use a switch via iPWGSPECTRUM
+      // or a config file
+      pwg4spec->SetAnalysisType(AliAnalysisTaskJetSpectrum::kAnaMC);
+      //      if(iAODanalysis)pwg4spec->SetAODInput(kTRUE);
+      pwg4spec->SetDebugLevel(11); 
+      //      pwg4spec->SetBranchRec("jetsMC"); 
+      //      pwg4spec->SetBranchGen("jetsAOD"); 
+      mgr->AddTask(pwg4spec);
+
+      AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer, "histos_pwg4spec.root");
+      // coutput1_Spec->SetSpecialOutput();
+      // Dummy AOD output container for jet analysis (no client yet)
+      c_aodSpec = mgr->CreateContainer("cAODjetSpec", TTree::Class(),
+                                       AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectInput  (pwg4spec,  0, cinput);    
+      // mgr->ConnectInput  (pwg4spec,  0, c_aodjet);    
+      mgr->ConnectOutput (pwg4spec,  0, c_aodSpec );
+      //       mgr->ConnectOutput (pwg4spec,  1, coutput1_Spec );
+       
+      mgr->ConnectOutput (pwg4spec,  1, coutput1_Spec );
+    }   
+
     
     // Particle correlation analysis
     if (iPWG4UE) {
-      AliAnalysisTaskUE* ueana = new  AliAnalysisTaskUE("Undelying Event");
+      AliAnalysisTaskUE* ueana = new  AliAnalysisTaskUE("Underlying Event");
       
 
       // default parameters use a switch via iPWGUE
@@ -242,8 +332,8 @@ void AnalysisTrainCAF()
 
       AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, "histosUE.root");
 
-      //      mgr->ConnectInput  (ueana,  0, cinput);    
-      mgr->ConnectInput  (ueana,  0, c_aodjet);    
+      mgr->ConnectInput  (ueana,  0, cinput);    
+      //      mgr->ConnectInput  (ueana,  0, c_aodjet);    
       mgr->ConnectOutput (ueana,     0, coutput1_UE );
     }   
 
@@ -251,6 +341,7 @@ void AnalysisTrainCAF()
     //    
     if (mgr->InitAnalysis()) {
        mgr->PrintStatus();
-       mgr->StartAnalysis("proof",dataset.Data(), 20000);
+       if(bPROOF)mgr->StartAnalysis("proof",dataset.Data(), nEvents,nOffset);
+       else mgr->StartAnalysis("local",chain);
     }   
 }