Added tasks from Sona
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 9 Aug 2009 16:36:43 +0000 (16:36 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 9 Aug 2009 16:36:43 +0000 (16:36 +0000)
PWG4/CMake_libPWG4JetTasks.txt
PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskJetCorrections.h [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskThreeJets.h [new file with mode: 0644]
PWG4/PWG4JetTasksLinkDef.h
PWG4/libPWG4JetTasks.pkg
PWG4/macros/AddTaskJetCorrections.C [new file with mode: 0644]
PWG4/macros/AddTaskThreeJets.C [new file with mode: 0644]
PWG4/macros/AnalysisTrainCAF.C

index c730b3c..060a636 100644 (file)
@@ -4,6 +4,8 @@ set(SRCS
        JetTasks/AliAnalysisTaskUE.cxx
        JetTasks/AliAnalysisTaskJetSpectrum.cxx
        JetTasks/AliAnalysisTaskJFSystematics.cxx
+       JetTasks/AliAnalysisTaskJetCorrections.cxx;
+       JetTasks/AliAnalysisTaskThreeJets.cxx;
        JetTasks/AliAnalysisHelperJetTasks.cxx
        JetTasks/AliAnaESDSpectraQA.cxx
        JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx
diff --git a/PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx b/PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx
new file mode 100644 (file)
index 0000000..e2e1dd0
--- /dev/null
@@ -0,0 +1,1002 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TInterpreter.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TList.h>
+#include <TTree.h>
+#include <TBranch.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
+#include <TRefArray.h>
+#include <TArrayD.h>
+#include <fstream>
+#include <TVector3.h>
+#include <TVectorD.h>
+#include <TMatrixDSym.h>
+#include <TMatrixDSymEigen.h>
+#include <TStyle.h>
+#include <TProfile.h>
+#include  "TDatabasePDG.h"
+
+#include "AliAnalysisTaskJetCorrections.h"
+#include "AliAnalysisManager.h"
+#include "AliJetFinder.h"
+#include "AliJetReader.h"
+#include "AliJetHeader.h"
+#include "AliJetReaderHeader.h"
+#include "AliUA1JetHeaderV1.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.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 "AliAODPid.h"
+#include "AliExternalTrackParam.h"
+
+#include "AliAnalysisTaskJetSpectrum.h"
+#include "AliAnalysisTaskThreeJets.h"
+
+#include "AliAnalysisHelperJetTasks.h"
+
+
+ClassImp(AliAnalysisTaskJetCorrections)
+
+  AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections() : AliAnalysisTaskSE(),
+                                                                  
+                                                                  fAOD(0x0),
+                                                                  
+                                                                  fBranchRec(""),
+                                                                  fBranchGen(""),
+                                                                  
+                                                                  fUseAODInput(kFALSE),
+
+                                                                  fR(0x0),
+                                                                  fList(0x0),
+                                                                  
+                                                                  fGlobVar(1),
+                                                                  fXsection(1),
+                                                                  
+                                                                  fhEGen(0x0),
+                                                                  fhERec(0x0),
+                                                                  
+                                                                  fhEGenRest(0x0),
+                                                                  fhERecRest(0x0),
+
+                                                                  fhEsumGenRest(0x0),
+                                                                  fhEsumRecRest(0x0),
+                                                                  
+                                                                  fhE2vsE1Gen(0x0),
+                                                                  fhE2vsE1Rec(0x0),
+                                                                  fhE2E1vsEsumGen(0x0),
+                                                                  fhE2E1vsEsumRec(0x0),
+                                                                  fhE2E1vsE1Gen(0x0),
+                                                                  fhE2E1vsE1Rec(0x0),
+                                                                  fhE2E1vsdPhiGen(0x0),
+                                                                  fhE2E1vsdPhiRec(0x0),
+
+                                                                  fhTrackBalance2(0x0),
+                                                                  fhTrackBalance3(0x0),
+
+                                                                  fhEt1Et22(0x0),
+                                                                  fhEt1Et23(0x0)
+
+{
+  
+  for (Int_t i = 0; i < 3; i++)
+    {
+      fhECorrJet10[i] = 0;    
+      fhECorrJet05[i] = 0;    
+      fhECorrJet01[i] = 0;    
+      fhECorrJet001[i] = 0;
+      fhdEvsErec10[i] = 0;
+      fhdEvsErec05[i] = 0;
+      fhdEvsErec01[i] = 0;
+      fhdEvsErec001[i] = 0;
+      fhdPhidEta10[i] = 0;
+      fhdPhidEta05[i] = 0;
+      fhdPhidEta01[i] = 0;
+      fhdPhidEta001[i] = 0;
+      fhdPhidEtaPt10[i] = 0;
+      fhdPhidEtaPt05[i] = 0;
+      fhdPhidEtaPt01[i] = 0;
+      fhdPhidEtaPt001[i] = 0;
+    }
+}
+
+AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections(const char * name):
+  AliAnalysisTaskSE(name),
+  
+  fAOD(0x0),
+  
+  fBranchRec(""),
+  fBranchGen(""),
+  
+  fUseAODInput(kFALSE),
+
+  fR(0x0),
+  fList(0x0),
+
+  fGlobVar(1),
+  fXsection(1),
+
+  fhEGen(0x0),
+  fhERec(0x0),
+
+  fhEGenRest(0x0),
+  fhERecRest(0x0),
+  
+  fhEsumGenRest(0x0),
+  fhEsumRecRest(0x0),
+  
+  fhE2vsE1Gen(0x0),
+  fhE2vsE1Rec(0x0),
+  fhE2E1vsEsumGen(0x0),
+  fhE2E1vsEsumRec(0x0),
+  fhE2E1vsE1Gen(0x0),
+  fhE2E1vsE1Rec(0x0),
+  fhE2E1vsdPhiGen(0x0),
+  fhE2E1vsdPhiRec(0x0),
+
+  fhTrackBalance2(0x0),
+  fhTrackBalance3(0x0),
+
+  fhEt1Et22(0x0),
+  fhEt1Et23(0x0)
+{
+  for (Int_t i = 0; i < 3; i++)
+    {
+      fhECorrJet10[i] = 0;    
+      fhECorrJet05[i] = 0;    
+      fhECorrJet01[i] = 0;    
+      fhECorrJet001[i] = 0;
+      fhdEvsErec10[i] = 0;
+      fhdEvsErec05[i] = 0;
+      fhdEvsErec01[i] = 0;
+      fhdEvsErec001[i] = 0;
+      fhdPhidEta10[i] = 0;
+      fhdPhidEta05[i] = 0;
+      fhdPhidEta01[i] = 0;
+      fhdPhidEta001[i] = 0;
+      fhdPhidEtaPt10[i] = 0;
+      fhdPhidEtaPt05[i] = 0;
+      fhdPhidEtaPt01[i] = 0;
+      fhdPhidEtaPt001[i] = 0;
+    }
+  DefineOutput(1, TList::Class());
+}
+
+
+
+Bool_t AliAnalysisTaskJetCorrections::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  UInt_t   ntrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "pyxsec.root");
+    }
+    TFile *fxsec = TFile::Open(fileName.Data());
+    if(!fxsec){
+      Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
+      // no a severe condition
+      return kTRUE;
+    }
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree){
+      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+    }
+    xtree->SetBranchAddress("xsection",&fXsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+  }
+  
+  return kTRUE;
+}
+
+
+//___________________________________________________________________________________________________________________________________
+void AliAnalysisTaskJetCorrections::UserCreateOutputObjects()
+{
+  //
+  // Create the output container
+  //
+  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
+
+  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;
+    }    
+  }
+
+  printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
+
+  fList = new TList();
+
+  fhEGen = new TH1F("EGen", "", 100, 0, 200);
+  fhEGen->Sumw2();
+  fList->Add(fhEGen);
+
+  fhERec = new TH1F("ERec", "", 100, 0, 200);
+  fhERec->Sumw2();
+  fList->Add(fhERec);
+
+  fhEGenRest = new TH1F("EGenRest", "", 100, 0, 200);
+  fhEGenRest->Sumw2();
+  fList->Add(fhEGenRest);
+
+  fhERecRest = new TH1F("ERecRest", "", 100, 0, 200);
+  fhERecRest->Sumw2();
+  fList->Add(fhERecRest);
+
+  fhEsumGenRest = new TH1F("EsumGenRest", "", 100, 0, 200);
+  fhEsumGenRest->Sumw2();
+  fList->Add(fhEsumGenRest);
+
+  fhEsumRecRest = new TH1F("EsumRecRest", "", 100, 0, 200);
+  fhEsumRecRest->Sumw2();
+  fList->Add(fhEsumRecRest);
+
+  fhE2vsE1Gen = new TH2F("E2vsE1Gen", "", 100, 0, 200, 100, 0, 200);
+  fhE2vsE1Gen->Sumw2();
+  fList->Add(fhE2vsE1Gen);
+
+  fhE2vsE1Rec = new TH2F("E2vsE1Rec", "", 100, 0, 200, 100, 0, 200);
+  fhE2vsE1Rec->Sumw2();
+  fList->Add(fhE2vsE1Rec);
+
+  fhE2E1vsEsumGen = new TH2F("E2E1vsEsumGen", "", 100, 0, 200, 25, 0, 1);
+  fhE2E1vsEsumGen->Sumw2();
+  fList->Add(fhE2E1vsEsumGen);
+
+  fhE2E1vsEsumRec = new TH2F("E2E1vsEsumRec", "", 100, 0, 200, 25, 0, 1);
+  fhE2E1vsEsumRec->Sumw2();
+  fList->Add(fhE2E1vsEsumRec);
+
+  fhE2E1vsE1Gen = new TH2F("E2E1vsE1Gen", "", 100, 0, 200, 25, 0, 1);
+  fhE2E1vsE1Gen->Sumw2();
+  fList->Add(fhE2E1vsE1Gen);
+
+  fhE2E1vsE1Rec = new TH2F("E2E1vsE1Rec", "", 100, 0, 200, 25, 0, 1);
+  fhE2E1vsE1Rec->Sumw2();
+  fList->Add(fhE2E1vsE1Rec);
+
+  fhE2E1vsdPhiGen =  new TH2F("E2E1vsdPhiGen", "", 64, -3.20, 3.20, 25, 0, 1);
+  fList->Add(fhE2E1vsdPhiGen);
+
+  fhE2E1vsdPhiRec =  new TH2F("E2E1vsdPhiRec", "", 64, -3.20, 3.20, 25, 0, 1);
+  fList->Add(fhE2E1vsdPhiRec);
+  
+  fhTrackBalance2 = new TH2F("TrackBalance2", "", 60, 0, 30, 60, 0, 30);
+  fhTrackBalance2->Sumw2();
+  fList->Add(fhTrackBalance2);
+  
+  fhTrackBalance3 = new TH2F("TrackBalance3", "", 60, 0, 30, 60, 0, 30);
+  fhTrackBalance3->Sumw2();
+  fList->Add(fhTrackBalance3);
+
+  fhEt1Et22 = new TH2F("Et1Et22", "", 100, 0, 50, 100, 0, 50);
+  fhEt1Et22->Sumw2();
+  fList->Add(fhEt1Et22);
+
+  fhEt1Et23 = new TH2F("Et1Et23", "", 100, 0, 50, 100, 0, 50);
+  fhEt1Et23->Sumw2();
+  fList->Add(fhEt1Et23);
+
+  for(Int_t i = 0; i < 3; i++)
+    {
+      fhECorrJet10[i] = new TProfile(Form("ECorrJet10%d", i+1), "", 100, 0, 200, 0, 10);
+      fhECorrJet10[i]->SetXTitle("E_{rec} [GeV]");
+      fhECorrJet10[i]->SetYTitle("C=E_{gen}/E_{rec}");
+      fhECorrJet10[i]->Sumw2();
+
+      fhECorrJet05[i] = new TProfile(Form("ECorrJet05%d", i+1), "", 100, 0, 200, 0, 10);
+      fhECorrJet05[i]->SetXTitle("E_{rec} [GeV]");
+      fhECorrJet05[i]->SetYTitle("C=E_{gen}/E_{rec}");
+      fhECorrJet05[i]->Sumw2();
+
+      fhECorrJet01[i] = new TProfile(Form("ECorrJet01%d", i+1), "", 100, 0, 200, 0, 10);
+      fhECorrJet01[i]->SetXTitle("E_{rec} [GeV]");
+      fhECorrJet01[i]->SetYTitle("C=E_{gen}/E_{rec}");
+      fhECorrJet01[i]->Sumw2();
+
+      fhECorrJet001[i] = new TProfile(Form("ECorrJet001%d", i+1), "", 100, 0, 200, 0, 10);
+      fhECorrJet001[i]->SetXTitle("E_{rec} [GeV]");
+      fhECorrJet001[i]->SetYTitle("C=E_{gen}/E_{rec}");
+      fhECorrJet001[i]->Sumw2();
+
+      fhdEvsErec10[i] = new TProfile(Form("dEvsErec10_%d", i+1),"", 100, 0, 200, -1, 10);
+      fhdEvsErec10[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
+      fhdEvsErec10[i]->SetXTitle("E_{rec} [GeV]");
+      fhdEvsErec10[i]->Sumw2();
+
+      fhdEvsErec05[i] = new TProfile(Form("dEvsErec05_%d", i+1),"", 100, 0, 200, -1, 10);
+      fhdEvsErec05[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
+      fhdEvsErec05[i]->SetXTitle("E_{rec} [GeV]");
+      fhdEvsErec05[i]->Sumw2();
+
+      fhdEvsErec01[i] = new TProfile(Form("dEvsErec01_%d", i+1),"", 100, 0, 200, -1, 10);
+      fhdEvsErec01[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
+      fhdEvsErec01[i]->SetXTitle("E_{rec} [GeV]");
+      fhdEvsErec01[i]->Sumw2();
+
+      fhdEvsErec001[i] = new TProfile(Form("dEvsErec001_%d", i+1),"", 100, 0, 200, -1, 10);
+      fhdEvsErec001[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
+      fhdEvsErec001[i]->SetXTitle("E_{rec} [GeV]");
+      fhdEvsErec001[i]->Sumw2();
+
+      fhdPhidEta10[i] = new TH2F(Form("dPhidEta10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEta10[i]->SetXTitle("#phi [rad]");
+      fhdPhidEta10[i]->SetYTitle("#eta");
+      fhdPhidEta10[i]->Sumw2();
+
+      fhdPhidEta05[i] = new TH2F(Form("dPhidEta05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEta05[i]->SetXTitle("#phi [rad]");
+      fhdPhidEta05[i]->SetYTitle("#eta");
+      fhdPhidEta05[i]->Sumw2();
+
+      fhdPhidEta01[i] = new TH2F(Form("dPhidEta01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEta01[i]->SetXTitle("#phi [rad]");
+      fhdPhidEta01[i]->SetYTitle("#eta");
+      fhdPhidEta01[i]->Sumw2();
+
+      fhdPhidEta001[i] = new TH2F(Form("dPhidEta001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEta001[i]->SetXTitle("#phi [rad]");
+      fhdPhidEta001[i]->SetYTitle("#eta");
+      fhdPhidEta001[i]->Sumw2();
+
+      fhdPhidEtaPt10[i] = new TH2F(Form("dPhidEtaPt10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEtaPt10[i]->SetXTitle("#phi [rad]");
+      fhdPhidEtaPt10[i]->SetYTitle("#eta");
+      fhdPhidEtaPt10[i]->Sumw2();
+
+      fhdPhidEtaPt05[i] = new TH2F(Form("dPhidEtaPt05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEtaPt05[i]->SetXTitle("#phi [rad]");
+      fhdPhidEtaPt05[i]->SetYTitle("#eta");
+      fhdPhidEtaPt05[i]->Sumw2();
+
+      fhdPhidEtaPt01[i] = new TH2F(Form("dPhidEtaPt01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEtaPt01[i]->SetXTitle("#phi [rad]");
+      fhdPhidEtaPt01[i]->SetYTitle("#eta");
+      fhdPhidEtaPt01[i]->Sumw2();
+
+      fhdPhidEtaPt001[i] = new TH2F(Form("dPhidEtaPt001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
+      fhdPhidEtaPt001[i]->SetXTitle("#phi [rad]");
+      fhdPhidEtaPt001[i]->SetYTitle("#eta");
+      fhdPhidEtaPt001[i]->Sumw2();
+
+      fList->Add(fhECorrJet10[i]);
+      fList->Add(fhECorrJet05[i]);
+      fList->Add(fhECorrJet01[i]);
+      fList->Add(fhECorrJet001[i]);
+      fList->Add(fhdEvsErec10[i]);
+      fList->Add(fhdEvsErec05[i]);
+      fList->Add(fhdEvsErec01[i]);
+      fList->Add(fhdEvsErec001[i]);
+      fList->Add(fhdPhidEta10[i]);
+      fList->Add(fhdPhidEta05[i]);
+      fList->Add(fhdPhidEta01[i]);
+      fList->Add(fhdPhidEta001[i]);
+      fList->Add(fhdPhidEtaPt10[i]);
+      fList->Add(fhdPhidEtaPt05[i]);
+      fList->Add(fhdPhidEtaPt01[i]);
+      fList->Add(fhdPhidEtaPt001[i]);
+    }
+    
+  Printf("UserCreateOutputObjects finished\n");
+}
+
+//__________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskJetCorrections::Init()
+{
+  printf("AliAnalysisJetCut::Init() \n");
+}
+
+//____________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskJetCorrections::UserExec(Option_t * )
+{
+//  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
+
+  //create an AOD handler
+  AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+  
+  if(!aodH)
+    {
+      Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
+      return;
+    }
+  
+  AliMCEvent* mcEvent =MCEvent();
+  if(!mcEvent){
+    Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
+    return;
+  }
+  
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+ //primary vertex
+  AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
+  if(!pvtx) return;
+
+ AliAODJet genJets[kMaxJets];
+ Int_t nGenJets = 0;
+ AliAODJet recJets[kMaxJets];
+ Int_t nRecJets = 0;
+
+  //array of reconstructed jets from the AOD input
+ 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;
+ }
+ // reconstructed jets
+ nRecJets = aodRecJets->GetEntries(); 
+ nRecJets = TMath::Min(nRecJets, kMaxJets);
+ for(int ir = 0;ir < nRecJets;++ir)
+   {
+     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
+     if(!tmp)continue;
+     recJets[ir] = *tmp;
+    }
+ // If we set a second branch for the input jets fetch this
+ TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
+ if(!aodGenJets)
+   {
+     printf("NO MC jets branch with name %s  Found \n",fBranchGen.Data());
+     return;
+   }
+ //   //Generated jets
+ nGenJets = aodGenJets->GetEntries();
+ nGenJets = TMath::Min(nGenJets, kMaxJets);
+ for(Int_t ig =0 ; ig < nGenJets; ++ig)
+   {
+     AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
+     if(!tmp)continue;
+     genJets[ig] = * tmp;
+   }
+
+ Double_t eRec[kMaxJets];
+ Double_t eGen[kMaxJets];
+ Double_t eRecRest[kMaxJets];
+ Double_t eGenRest[kMaxJets];
+
+// AliAODJet jetRec[kMaxJets];
+ AliAODJet jetGen[kMaxJets];
+ Int_t idxRec[kMaxJets];
+ Int_t idxGen[kMaxJets];
+
+// Double_t EsumRec = 0;
+ // Double_t EsumGen =0;
+ TLorentzVector vRec[kMaxJets];
+ TLorentzVector vGen[kMaxJets];
+
+ TLorentzVector vsumRec;
+ TLorentzVector vsumGen;
+ TVector3 pRec[kMaxJets];
+ TVector3 pGen[kMaxJets];
+ // HISTOS FOR MC
+ Int_t nGenSel = 0;
+ Int_t counter = 0;
+ Int_t tag = 0;
+ AliAODJet selJets[kMaxJets];
+ // loop for applying the separation cut
+ for(Int_t i = 0; i < nGenJets; i++)
+   {
+     if(nGenJets == 1)
+       {
+        selJets[nGenSel] = genJets[i];
+        nGenSel++;
+       }
+     else
+       {
+        counter = 0;
+        tag = 0;
+        for(Int_t j = 0; j < nGenJets; j++)
+          {
+            if(i!=j)
+              {
+                Double_t dRij = genJets[i].DeltaR(&genJets[j]);
+                counter++;
+                if(dRij > 2*fR) tag++;
+              }
+          }
+        if(counter!=0)
+          {
+            if(tag/counter == 1)
+              {
+                selJets[nGenSel] = genJets[i];
+                nGenSel++;
+              }
+          }
+       }
+   }
+ for (Int_t gj = 0; gj < nGenSel; gj++)
+   {
+     eGen[gj] = selJets[gj].E();
+     fhEGen->Fill(eGen[gj], fXsection);
+   }
+ TMath::Sort(nGenSel, eGen, idxGen);
+ for (Int_t ig = 0; ig < nGenSel; ig++)
+   jetGen[ig] = selJets[idxGen[ig]];
+ //rest frame MC jets
+ for (Int_t i = 0; i < nGenSel; ++i)
+   {
+     vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
+     pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
+     vsumGen += vGen[i];
+   }
+ if(nGenSel > 1 && pGen[0].DeltaPhi(pGen[1]) > 2.8)
+   {
+     fhE2vsE1Gen->Fill(jetGen[0].E(), jetGen[1].E(), fXsection);
+     fhE2E1vsEsumGen->Fill(jetGen[0].E()+jetGen[1].E(), TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
+     fhE2E1vsE1Gen->Fill(jetGen[0].E(),  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
+     Double_t deltaPhi = (jetGen[0].Phi()-jetGen[1].Phi());
+     if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
+     if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
+     fhE2E1vsdPhiGen->Fill(deltaPhi,  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
+   }
+     
+ Double_t fPxGen = vsumGen.Px();
+ Double_t fPyGen = vsumGen.Py();
+ Double_t fPzGen = vsumGen.Pz();
+ Double_t fEGen = vsumGen.E();
+ Double_t eSumGenRest = 0; 
+ for (Int_t j = 0; j < nGenSel; j++)
+   {
+     vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
+     eGenRest[j] = vGen[j].E();
+     if(nGenSel > 1)
+       fhEGenRest->Fill(eGenRest[j], fXsection);
+     eSumGenRest += eGenRest[j];
+   }
+ if(nGenSel > 1)    
+   fhEsumGenRest->Fill(eSumGenRest, fXsection);
+ //END VARIABLES FOR MC JETS ---------------
+ //   }
+ //AOD JET VARIABLES
+ Int_t nRecSel = 0;
+ Int_t counter1 = 0;
+ Int_t tag1 = 0;
+ AliAODJet recSelJets[kMaxJets];
+ for(Int_t i = 0; i < nRecJets; i++)
+   {
+     if(nRecJets == 1)
+       {
+        recSelJets[nRecSel] = recJets[i];
+        nRecSel++;
+       }
+     else
+       {
+        counter1 = 0;
+        tag1 = 0;
+        for(Int_t j = 0; j < nRecJets; j++)
+          {
+            if(i!=j)
+              {
+                Double_t dRij = recJets[i].DeltaR(&recJets[j]);
+                counter1++;
+                if(dRij > 2*fR) tag1++;
+              }
+          }
+        if(counter1!=0)
+          {
+            if(tag1/counter1 == 1)
+              {
+                recSelJets[nRecSel] = recJets[i];
+                nRecSel++;
+              }
+          }
+       }
+   } 
+  
+ if(nRecSel == 0) return;
+ Printf("******NUMBER OF JETS AFTER DELTA R CUT : %d **********\n", nRecSel);
+ //sort rec/gen jets by energy in C.M.S
+ AliAODJet jetRecTmp[kMaxJets];
+ Int_t nAccJets = 0;
+ Double_t jetTrackPt[kTracks];
+ TLorentzVector jetTrackTmp[kTracks];
+ Int_t nTracks = 0;
+ for (Int_t rj = 0; rj < nRecSel; rj++)
+   {
+     TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(recSelJets[rj].GetRefTracks());
+     if(!jetTracksAOD) continue;
+     if(jetTracksAOD->GetEntries() < 3) continue;
+     Int_t nJetTracks = 0;
+     for(Int_t j = 0; j < jetTracksAOD->GetEntries(); j++)
+       {
+        AliAODTrack * track = dynamic_cast<AliAODTrack*>(jetTracksAOD->At(j));
+        if(!track) continue;
+        Double_t cv[21];
+        track->GetCovarianceXYZPxPyPz(cv);
+        if(cv[14] > 1000.) continue;
+        jetTrackPt[nTracks] = track->Pt();
+        jetTrackTmp[nTracks].SetPxPyPzE(track->Px(),track->Py(),track->Pz(),track->E());
+        nTracks++;
+        nJetTracks++;
+       }
+     if(nJetTracks < 4) continue;
+     jetRecTmp[nAccJets] = recSelJets[rj];
+     eRec[nAccJets] = recSelJets[rj].E();
+     fhERec->Fill(eRec[nAccJets], fXsection);
+     nAccJets++;
+   }
+
+ if(nAccJets == 0) return;
+ if(nTracks == 0) return;
+
+ Printf(" ************ Number of accepted jets : %d ************ \n", nAccJets);
+
+ AliAODJet jetRecAcc[kMaxJets];
+ TMath::Sort(nAccJets, eRec, idxRec);
+ for (Int_t rj = 0; rj < nAccJets; rj++)
+   jetRecAcc[rj] = jetRecTmp[idxRec[rj]];
+ //rest frame for reconstructed jets
+ for (Int_t i = 0; i < nAccJets; i++)
+   {
+     vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
+     pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
+     vsumRec += vRec[i];
+   }
+
+ //check balance of two leading hadrons, deltaPhi > 2.
+ Int_t idxTrack[kTracks];
+ TMath::Sort(nTracks, jetTrackPt, idxTrack);
+
+ TLorentzVector jetTrack[kTracks];
+ for(Int_t iTr = 0; iTr < nTracks; iTr++)
+   jetTrack[iTr] = jetTrackTmp[idxTrack[iTr]];
+ Int_t n = 1;
+ while(jetTrack[0].DeltaPhi(jetTrack[n]) < 2.8)
+   n++;
+
+ Double_t et1 = 0;
+ Double_t et2 = 0;
+ for(Int_t iTr = 0; iTr < nTracks; iTr++)
+   {
+     if(TMath::Abs(jetTrack[0].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != 0)
+       et1 += jetTrack[iTr].Et();
+
+     if(TMath::Abs(jetTrack[n].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != n)
+       et2 += jetTrack[iTr].Et();
+   }
+
+ if(nAccJets == 2)
+   {
+     fhTrackBalance2->Fill(jetTrack[0].Et(), jetTrack[n].Et());
+     fhEt1Et22->Fill(et1, et2);
+   }
+ if(nAccJets == 3)
+   {
+     fhTrackBalance3->Fill(jetTrack[0].Et(), jetTrack[n].Et());
+     fhEt1Et23->Fill(et1, et2);
+   }
+    
+ if(nAccJets > 1 && pRec[0].DeltaPhi(pRec[1]) > 2.8)
+   {
+     fhE2vsE1Rec->Fill(jetRecAcc[0].E(), jetRecAcc[1].E(), fXsection);
+     fhE2E1vsEsumRec->Fill(jetRecAcc[0].E()+jetRecAcc[1].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
+     fhE2E1vsE1Rec->Fill(jetRecAcc[0].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
+     Double_t deltaPhi = (jetRecAcc[0].Phi()-jetRecAcc[1].Phi());
+     if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
+     if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
+     fhE2E1vsdPhiRec->Fill(-1*deltaPhi, TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
+   }
+ Double_t fPx = vsumRec.Px();
+ Double_t fPy = vsumRec.Py();
+ Double_t fPz = vsumRec.Pz();
+ Double_t fE = vsumRec.E();
+
+ Double_t eSumRecRest = 0;
+ for (Int_t j = 0; j < nAccJets; j++)
+   {
+     vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
+     eRecRest[j] = vRec[j].E();
+     if(nAccJets > 1)
+       fhERecRest->Fill(eRecRest[j], fXsection);
+     eSumRecRest += eRecRest[j];
+   }
+ if(nAccJets > 1)
+   fhEsumRecRest->Fill(eSumRecRest, fXsection);
+ // 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;
+ }
+ Double_t dR = 1.4;
+ if(nAccJets == nGenSel)
+   {
+    AliAnalysisHelperJetTasks::GetClosestJets(jetGen,nGenSel,jetRecAcc,nAccJets,
+                                                iGenIndex,iRecIndex,0);
+      
+     for(int ir = 0;ir < nAccJets;ir++){
+       Int_t ig = iGenIndex[ir];
+       if(ig>=0&&ig<nGenSel){
+        Double_t dPhi = TMath::Abs(jetRecAcc[ir].Phi()-jetGen[ig].Phi());
+        if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
+        if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
+        Double_t sigma = TMath::Abs(jetGen[ig].E()-jetRecAcc[ir].E())/jetGen[ig].E();
+        dR = jetRecAcc[ir].DeltaR(&jetGen[ig]);
+        if(dR < 2*fR && dR >= fR) 
+          {
+            switch(nAccJets)
+              {
+              case 1:
+                {
+                  fhdEvsErec10[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet10[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 2:
+                {
+                  fhdEvsErec10[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet10[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 3:
+                {
+                  fhdEvsErec10[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet10[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              }
+          }
+        if(dR < fR && dR >= 0.1) 
+          {
+            switch(nAccJets)
+              {
+              case 1:
+                {
+                  fhdEvsErec05[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet05[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 2:
+                {
+                  fhdEvsErec05[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet05[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 3:
+                {
+                  fhdEvsErec05[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet05[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              }
+          }
+        if(dR < 0.1 && dR >= 0.01)
+          {
+            switch(nAccJets)
+              {
+              case 1:
+                {
+                  fhdEvsErec01[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet01[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 2:
+                {
+                  fhdEvsErec01[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet01[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              case 3:
+                {
+                  fhdEvsErec01[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+                  fhECorrJet01[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+                  for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                    {
+                      fhdPhidEtaPt01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                      fhdPhidEta01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                    }
+                }
+                break;
+              }
+          }
+        if(dR > 0.01) continue;
+        switch(nAccJets)
+          {
+          case 1:
+            {
+              fhECorrJet001[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+              fhdEvsErec001[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+              for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                {
+                  fhdPhidEtaPt001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                  fhdPhidEta001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                }
+            }
+            break;
+          case 2:
+            {
+              fhECorrJet001[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+              fhdEvsErec001[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+              for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                {
+                  fhdPhidEtaPt001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                  fhdPhidEta001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                }
+            }
+            break;
+          case 3:
+            {
+              fhECorrJet001[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
+              fhdEvsErec001[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
+              for(Int_t iTr = 0; iTr < nTracks; iTr++)
+                {
+                  fhdPhidEtaPt001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
+                  fhdPhidEta001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
+                }
+            }
+            break;
+          }
+       }
+     }
+     // loop over reconstructed jets
+   }
+
+ Printf("%s:%d",(char*)__FILE__,__LINE__);
+ PostData(1, fList);
+ Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
+  
+}
+
+
+//__________________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
+{
+  printf("AnalysisJetCorrelation::Terminate()");
+
+} 
+
+//_______________________________________User defined functions_____________________________________________________________________________________
+void AliAnalysisTaskJetCorrections::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, Int_t &nTracks)
+{
+  TVector3 psum;
+  Double_t psum1 = 0;
+  Double_t psum2 = 0;
+  Double_t thrust[kTracks];
+  Double_t th = -3;
+  
+  for(Int_t j = 0; j < nTracks; j++)
+    {
+      psum.SetXYZ(0., 0., 0.);
+      psum1 = 0;
+      psum2 = 0;
+      for(Int_t i = 0; i < nTracks; i++)
+       {
+         psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
+         psum2 += pTrack[i].Mag();
+         
+         if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
+         if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
+       }
+      
+      thrust[j] = psum1/psum2;
+      
+      if(th == thrust[j]) 
+       break;
+      
+      th = thrust[j];
+      
+      n01 = psum.Unit();
+    }
+}
+//______________________________________________________________________________________________________
+
+  
diff --git a/PWG4/JetTasks/AliAnalysisTaskJetCorrections.h b/PWG4/JetTasks/AliAnalysisTaskJetCorrections.h
new file mode 100644 (file)
index 0000000..2fd2f4d
--- /dev/null
@@ -0,0 +1,121 @@
+#ifndef ALIANALYSISTASKJETCORRECTIONS_H\r
+#define ALIANALYSISTASKJETCORRECTIONS_H\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include <iostream>\r
+#include <fstream>\r
+\r
+class AliJetFinder;\r
+class AliESDEvent;\r
+class AliAODEvent;\r
+class AliAODJet;\r
+class AliGenPythiaEventHeader;\r
+class AliAODPid;\r
+\r
+class TList;\r
+class TArrayD;\r
+class TChain;\r
+class TH1;\r
+class TH2;\r
+class TH1F;\r
+class TH2F;\r
+class TH2I;\r
+class TH3D;\r
+class TTree;\r
+class TProfile;\r
+class TLorentzVector;\r
+class TVector3;\r
+class TParticle;\r
+\r
+class AliAnalysisTaskJetCorrections : public AliAnalysisTaskSE\r
+{\r
+ public:\r
+  AliAnalysisTaskJetCorrections();\r
+  AliAnalysisTaskJetCorrections(const char * name);\r
+  virtual ~AliAnalysisTaskJetCorrections() {;}\r
+  \r
+  //Implementation of interface methods\r
+  virtual Bool_t Notify(); \r
+  virtual void UserCreateOutputObjects();\r
+  virtual void Init();\r
+  virtual void LocalInit() { Init(); };\r
+  virtual void UserExec(Option_t * option);\r
+  virtual void Terminate(Option_t * option);\r
+\r
+  virtual void SetAODInput(Bool_t b){fUseAODInput = b;}\r
+\r
+  virtual void SetBranchGen(char* c){fBranchGen = c;}\r
+  virtual void SetBranchRec(char* c){fBranchRec = c;}\r
+\r
+  virtual Double_t SetR(Double_t b){fR = b; return fR;} \r
+\r
+  virtual void GetThrustAxis(TVector3 &n01, TVector3 * pTrack, Int_t &nTracks);\r
+ private:\r
+  AliAnalysisTaskJetCorrections(const AliAnalysisTaskJetCorrections&);\r
+  AliAnalysisTaskJetCorrections& operator = (const AliAnalysisTaskJetCorrections&);\r
+\r
+  enum {kMaxJets = 6};\r
+  enum {kMaxEvents = 10};\r
+  enum {kJets = 3};\r
+  enum {kTracks = 1000};\r
+\r
+  AliAODEvent  *fAOD; // where we take the jets from can be input or output AOD\r
+  \r
+  TString       fBranchRec;  // AOD branch name for reconstructed\r
+  TString       fBranchGen;  // AOD brnach for genereated\r
+  \r
+  Bool_t        fUseAODInput;\r
+\r
+  Double_t fR;\r
+  TList * fList;\r
+\r
+  Int_t fGlobVar;\r
+  Double_t fXsection;\r
+\r
+\r
+  TH1F * fhEGen;\r
+  TH1F * fhERec;\r
+  TH1F * fhEGenRest;\r
+  TH1F * fhERecRest;\r
+  TH1F * fhEsumGenRest;\r
+  TH1F * fhEsumRecRest;\r
+\r
+  TH2F * fhE2vsE1Gen;\r
+  TH2F * fhE2vsE1Rec;\r
+  TH2F * fhE2E1vsEsumGen;\r
+  TH2F * fhE2E1vsEsumRec;\r
+  TH2F * fhE2E1vsE1Gen;\r
+  TH2F * fhE2E1vsE1Rec;\r
+  TH2F * fhE2E1vsdPhiGen;\r
+  TH2F * fhE2E1vsdPhiRec;\r
+\r
+  TH2F * fhTrackBalance2;\r
+  TH2F * fhTrackBalance3;\r
+\r
+  TH2F * fhEt1Et22;\r
+  TH2F * fhEt1Et23;\r
+\r
+  TProfile * fhECorrJet10[3];\r
+  TProfile * fhECorrJet05[3];\r
+  TProfile * fhECorrJet01[3];\r
+  TProfile * fhECorrJet001[3];\r
+  \r
+  TProfile * fhdEvsErec10[3];\r
+  TProfile * fhdEvsErec05[3];\r
+  TProfile * fhdEvsErec01[3];\r
+  TProfile * fhdEvsErec001[3];\r
+\r
+  TH2F * fhdPhidEta10[3];\r
+  TH2F * fhdPhidEta05[3];\r
+  TH2F * fhdPhidEta01[3];\r
+  TH2F * fhdPhidEta001[3];\r
+\r
+  TH2F * fhdPhidEtaPt10[3];\r
+  TH2F * fhdPhidEtaPt05[3];\r
+  TH2F * fhdPhidEtaPt01[3];\r
+  TH2F * fhdPhidEtaPt001[3];\r
+\r
+  ClassDef(AliAnalysisTaskJetCorrections, 1)\r
+};\r
+\r
+#endif\r
diff --git a/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx
new file mode 100644 (file)
index 0000000..96e8561
--- /dev/null
@@ -0,0 +1,1333 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TInterpreter.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TList.h>
+#include <TTree.h>
+#include <TBranch.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
+#include <TRefArray.h>
+#include <TArrayD.h>
+#include <fstream>
+#include <TVector3.h>
+#include <TVectorD.h>
+#include <TMatrixDSym.h>
+#include <TMatrixDSymEigen.h>
+#include <TStyle.h>
+#include <TProfile.h>
+#include  "TDatabasePDG.h"
+
+#include "AliAnalysisTaskThreeJets.h"
+#include "AliAnalysisManager.h"
+#include "AliJetFinder.h"
+#include "AliJetReader.h"
+#include "AliJetHeader.h"
+#include "AliJetReaderHeader.h"
+#include "AliUA1JetHeaderV1.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.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 "AliAODPid.h"
+#include "AliExternalTrackParam.h"
+
+#include "AliAnalysisTaskJetSpectrum.h"
+
+#include "AliAnalysisHelperJetTasks.h"
+
+
+ClassImp(AliAnalysisTaskThreeJets)
+
+  AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets() : AliAnalysisTaskSE(),
+                                                  
+                                                        fAOD(0x0),
+                                                
+                                                        fBranchRec(""),
+                                                        fBranchGen(""),
+                                               
+                                                        fUseAODInput(kFALSE),
+
+                                                        fR(0x0),
+                                                        fList(0x0),
+
+                                                        fGlobVar(1),
+                                                        fXsection(1),
+
+                                                        fhX3X4Rec(0x0),
+                                                        fhX3X4Gen(0x0),
+                                                        
+                                                        fhMu35Rec(0x0),
+                                                        fhMu34Rec(0x0),
+                                                        fhMu45Rec(0x0),
+    
+                                                        fhMu35Gen(0x0),
+                                                        fhMu34Gen(0x0),
+                                                        fhMu45Gen(0x0),
+       
+                                                        fhInOut(0x0),
+
+                                                        fhThrustRec2(0x0),
+                                                        fhThrustRec3(0x0),
+
+                                                        fhThrustGen2(0x0),
+                                                        fhThrustGen3(0x0),
+
+                                                        fhCGen2(0x0),
+                                                        fhCGen3(0x0),
+
+                                                        fhSGen2(0x0),
+                                                        fhSGen3(0x0),
+
+                                                        fhAGen2(0x0),
+                                                        fhAGen3(0x0),
+
+                                                        fhCRec2(0x0),
+                                                        fhCRec3(0x0),
+
+                                                        fhSRec2(0x0),
+                                                        fhSRec3(0x0),
+
+                                                        fhARec2(0x0),
+                                                        fhARec3(0x0),
+
+                                                        fhX3(0x0),
+                                                        fhX4(0x0),
+                                                        fhX5(0x0),
+
+                                                        fhXSec(0x0),
+
+                                                        fhX3X4Rec60(0x0),
+                                                        fhX3X4Rec60100(0x0),
+                                                        fhX3X4Rec100(0x0),
+
+                                                        fhX3X4Gen60(0x0),
+                                                        fhX3X4Gen60100(0x0),
+                                                        fhX3X4Gen100(0x0),
+
+                                                        fhdPhiThrustGen(0x0),
+                                                        fhdPhiThrustGenALL(0x0),
+
+                                                        fhdPhiThrustRec(0x0),
+                                                        fhdPhiThrustRecALL(0x0)
+{
+
+}
+
+AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets(const char * name):
+  AliAnalysisTaskSE(name),
+  
+  fAOD(0x0),
+  
+  fBranchRec(""),
+  fBranchGen(""),
+  
+  fUseAODInput(kFALSE),
+
+  fR(0x0),
+  fList(0x0),
+
+  fGlobVar(1),
+  fXsection(1),
+
+  fhX3X4Rec(0x0),
+  fhX3X4Gen(0x0),
+  
+  fhMu35Rec(0x0),
+  fhMu34Rec(0x0),
+  fhMu45Rec(0x0),
+  
+  fhMu35Gen(0x0),
+  fhMu34Gen(0x0),
+  fhMu45Gen(0x0),
+  
+  fhInOut(0x0),
+  
+  fhThrustRec2(0x0),
+  fhThrustRec3(0x0),
+  
+  fhThrustGen2(0x0),
+  fhThrustGen3(0x0),
+  
+  fhCGen2(0x0),
+  fhCGen3(0x0),
+  
+  fhSGen2(0x0),
+  fhSGen3(0x0),
+
+  fhAGen2(0x0),
+  fhAGen3(0x0),
+  
+  fhCRec2(0x0),
+  fhCRec3(0x0),
+  
+  fhSRec2(0x0),
+  fhSRec3(0x0),
+  
+  fhARec2(0x0),
+  fhARec3(0x0),
+  
+  fhX3(0x0),
+  fhX4(0x0),
+  fhX5(0x0),
+  
+  fhXSec(0x0),
+  
+  fhX3X4Rec60(0x0),
+  fhX3X4Rec60100(0x0),
+  fhX3X4Rec100(0x0),
+  
+  fhX3X4Gen60(0x0),
+  fhX3X4Gen60100(0x0),
+  fhX3X4Gen100(0x0),
+  
+  fhdPhiThrustGen(0x0),
+  fhdPhiThrustGenALL(0x0),
+  
+  fhdPhiThrustRec(0x0),
+  fhdPhiThrustRecALL(0x0)
+{
+  DefineOutput(1, TList::Class());
+}
+
+
+
+Bool_t AliAnalysisTaskThreeJets::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  UInt_t   ntrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "pyxsec.root");
+    }
+    TFile *fxsec = TFile::Open(fileName.Data());
+    if(!fxsec){
+      Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
+      // no a severe condition
+      return kTRUE;
+    }
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree){
+      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+    }
+    xtree->SetBranchAddress("xsection",&fXsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+  }
+  
+  return kTRUE;
+}
+
+
+//___________________________________________________________________________________________________________________________________
+void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
+{
+  //
+  // Create the output container
+  //
+  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
+
+  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;
+    }    
+  }
+
+  printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
+
+  fList = new TList();
+
+  fhX3X4Gen = new TH2F("X3vsX4Gen", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Gen->SetXTitle("X_{3}");
+  fhX3X4Gen->SetYTitle("X_{4}");
+  fhX3X4Gen->Sumw2();
+  fList->Add(fhX3X4Gen);
+     
+  fhX3X4Rec = new TH2F("X3vsX4Rec", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Rec->SetXTitle("X_{3}");
+  fhX3X4Rec->SetYTitle("X_{4}");
+  fhX3X4Rec->Sumw2();
+  fList->Add(fhX3X4Rec);
+     
+  fhMu35Rec = new TH1F("Mu35Rec", "", 20,0.1,0.8);
+  fhMu35Rec->Sumw2();
+  fhMu35Rec->SetXTitle("#mu_{35}");
+  fhMu35Rec->SetYTitle("#frac{dN}{d#mu_{35Rec}}");  
+  fList->Add(fhMu35Rec);
+         
+  fhMu34Rec = new TH1F("Mu34Rec", "", 20,0.5,1);
+  fhMu34Rec->Sumw2();
+  fhMu34Rec->SetXTitle("#mu_{34}");
+  fhMu34Rec->SetYTitle("#frac{dN}{d#mu_{34}}");
+  fList->Add(fhMu34Rec);
+  
+  fhMu45Rec = new TH1F("Mu45Rec", "", 20,0,0.65);
+  fhMu45Rec->Sumw2();
+  fhMu45Rec->SetXTitle("#mu_{45}");
+  fhMu45Rec->SetYTitle("#frac{dN}{d#mu_{45}}");
+  fList->Add(fhMu45Rec);
+     
+  fhMu35Gen = new TH1F("Mu35Gen", "", 20,0.1,0.8);
+  fhMu35Gen->Sumw2();
+  fhMu35Gen->SetXTitle("#mu_{35Gen}");
+  fhMu35Gen->SetYTitle("#frac{dN}{d#mu_{35Gen}}");  
+  fList->Add(fhMu35Gen);
+         
+  fhMu34Gen = new TH1F("Mu34Gen", "", 20,0.5,1);
+  fhMu34Gen->Sumw2();
+  fhMu34Gen->SetXTitle("#mu_{34Gen}");
+  fhMu34Gen->SetYTitle("#frac{dN}{d#mu_{34Gen}}");
+  fList->Add(fhMu34Gen);
+  
+  fhMu45Gen = new TH1F("Mu45Gen", "", 20,0,0.65);
+  fhMu45Gen->Sumw2();
+  fhMu45Gen->SetXTitle("#mu_{45Gen}");
+  fhMu45Gen->SetYTitle("#frac{dN}{d#mu_{45}}");
+  fList->Add(fhMu45Gen);
+
+  fhInOut = new TH1I("InOut", "", 6, 0, 6);
+  fhInOut->SetXTitle("#RecJets_{GenJets=3}");
+  fhInOut->SetYTitle("#Entries");
+  fList->Add(fhInOut);
+
+  fhThrustGen2 = new TH1F("ThrustGen2", "", 50, 0.5, 1);
+  fhThrustGen2->Sumw2();
+  fList->Add(fhThrustGen2);
+
+  fhThrustGen3 = new TH1F("ThrustGen3", "", 50, 0.5, 1);
+  fhThrustGen3->Sumw2();
+  fList->Add(fhThrustGen3);
+
+  fhThrustRec2 = new TH1F("ThrustRec2", "", 50, 0.5, 1);
+  fhThrustRec2->Sumw2();
+  fList->Add(fhThrustRec2);
+
+  fhThrustRec3 = new TH1F("ThrustRec3", "", 50, 0.5, 1);
+  fhThrustRec3->Sumw2();
+  fList->Add(fhThrustRec3);
+
+  fhCGen2 = new TH1F("CGen2", "", 100, 0, 1);
+  fhCGen2->Sumw2();
+  fList->Add(fhCGen2);
+
+  fhCGen3 = new TH1F("CGen3", "", 100, 0, 1);
+  fhCGen3->Sumw2();
+  fList->Add(fhCGen3);
+
+  fhCRec2 = new TH1F("CRec2", "", 100, 0, 1);
+  fhCRec2->Sumw2();
+  fList->Add(fhCRec2);
+
+  fhCRec3 = new TH1F("CRec3", "", 100, 0, 1);
+  fhCRec3->Sumw2();
+  fList->Add(fhCRec3);
+
+  fhSGen2 = new TH1F("SGen2", "", 100, 0, 1);
+  fList->Add(fhSGen2);
+
+  fhSGen3 = new TH1F("SGen3", "", 100, 0, 1);
+  fList->Add(fhSGen3);
+
+  fhSRec2 = new TH1F("SRec2", "", 100, 0, 1);
+  fList->Add(fhSRec2);
+
+  fhSRec3 = new TH1F("SRec3", "", 100, 0, 1);
+  fList->Add(fhSRec3);
+
+  fhAGen2 = new TH1F("AGen2", "", 50, 0, 0.5);
+  fList->Add(fhAGen2);
+
+  fhAGen3 = new TH1F("AGen3", "", 50, 0, 0.5);
+  fList->Add(fhAGen3);
+
+  fhARec2 = new TH1F("ARec2", "", 50, 0, 0.5);
+  fList->Add(fhARec2);
+
+  fhARec3 = new TH1F("ARec3", "", 50, 0, 0.5);
+  fList->Add(fhARec3);
+
+  fhX3 = new TH2F("X3", "", 22, 0.6, 1.02, 100, 0, 1);
+  fhX3->SetYTitle("|X_{3}^{MC} - X_{3}^{AOD}|/X_{3}^{MC}");
+  fhX3->SetXTitle("X_{3}");
+  fhX3->Sumw2();
+  fList->Add(fhX3);
+
+  fhX4 = new TH2F("X4", "",33, 0.4, 1.02, 100, 0, 1);
+  fhX4->SetYTitle("|X_{4}^{MC} - X_{4}^{AOD}|/X_{4}^{MC}");
+  fhX4->SetXTitle("X_{4}");
+  fhX4->Sumw2();
+  fList->Add(fhX4);
+
+  fhX5 = new TH2F("X5", "",100, 0., 1., 100, 0, 1);
+  fhX5->SetYTitle("|X_{5}^{MC} - X_{5}^{AOD}|/X_{5}^{MC}");
+  fhX5->SetXTitle("X_{5}");
+  fhX5->Sumw2();
+  fList->Add(fhX5);
+
+  fhXSec = new TProfile("XSec", "", 200, 0, 200, 0, 1);
+  fList->Add(fhXSec);
+
+  fhX3X4Rec60 = new TH2F("X3vsX4Rec60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Rec60->SetXTitle("X_{3}");
+  fhX3X4Rec60->SetYTitle("X_{4}");
+  fhX3X4Rec60->Sumw2();
+  fList->Add(fhX3X4Rec60);
+
+  fhX3X4Rec60100 = new TH2F("X3vsX4Rec60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Rec60100->SetXTitle("X_{3}");
+  fhX3X4Rec60100->SetYTitle("X_{4}");
+  fhX3X4Rec60100->Sumw2();
+  fList->Add(fhX3X4Rec60100);
+
+  fhX3X4Rec100 = new TH2F("X3vsX4Rec100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Rec100->SetXTitle("X_{3}");
+  fhX3X4Rec100->SetYTitle("X_{4}");
+  fhX3X4Rec100->Sumw2();
+  fList->Add(fhX3X4Rec100);
+
+  fhX3X4Gen60 = new TH2F("X3vsX4Gen60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Gen60->SetXTitle("X_{3}");
+  fhX3X4Gen60->SetYTitle("X_{4}");
+  fhX3X4Gen60->Sumw2();
+  fList->Add(fhX3X4Gen60);
+
+  fhX3X4Gen60100 = new TH2F("X3vsX4Gen60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Gen60100->SetXTitle("X_{3}");
+  fhX3X4Gen60100->SetYTitle("X_{4}");
+  fhX3X4Gen60100->Sumw2();
+  fList->Add(fhX3X4Gen60100);
+
+  fhX3X4Gen100 = new TH2F("X3vsX4Gen100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
+  fhX3X4Gen100->SetXTitle("X_{3}");
+  fhX3X4Gen100->SetYTitle("X_{4}");
+  fhX3X4Gen100->Sumw2();
+  fList->Add(fhX3X4Gen100);
+
+  fhdPhiThrustGen = new TH2F("dPhiThrustGen", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
+  fhdPhiThrustGen->Sumw2();
+  fList->Add(fhdPhiThrustGen);
+
+  fhdPhiThrustGenALL = new TH2F("dPhiThrustGenALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0,  150);
+  fhdPhiThrustGenALL->Sumw2();
+  fList->Add(fhdPhiThrustGenALL);
+
+  fhdPhiThrustRec = new TH2F("dPhiThrustRec", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
+  fhdPhiThrustRec->Sumw2();
+  fList->Add(fhdPhiThrustRec);
+
+  fhdPhiThrustRecALL = new TH2F("dPhiThrustRecALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
+  fhdPhiThrustRecALL->Sumw2();
+  fList->Add(fhdPhiThrustRecALL);
+    
+  Printf("UserCreateOutputObjects finished\n");
+}
+
+//__________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskThreeJets::Init()
+{
+  printf("AliAnalysisJetCut::Init() \n");
+}
+
+//____________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskThreeJets::UserExec(Option_t * )
+{
+  //  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
+
+  
+  //create an AOD handler
+  AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+  
+  if(!aodH)
+    {
+      Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
+      return;
+    }
+  
+  AliMCEvent* mcEvent =MCEvent();
+  if(!mcEvent){
+    Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
+    return;
+  }
+  
+  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+  
+  //primary vertex
+  AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
+  if(!pvtx) return;
+  
+  AliAODJet genJetsPythia[kMaxJets];
+  Int_t nPythiaGenJets = 0;
+  
+  AliAODJet genJets[kMaxJets];
+  Int_t nGenJets = 0;
+  
+  AliAODJet recJets[kMaxJets];
+  Int_t nRecJets = 0;
+  
+  //array of reconstructed jets from the AOD input
+  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;
+  }
+  
+  // reconstructed jets
+  nRecJets = aodRecJets->GetEntries(); 
+  nRecJets = TMath::Min(nRecJets, kMaxJets);
+  
+  for(int ir = 0;ir < nRecJets;++ir)
+    {
+     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
+     if(!tmp)continue;
+     recJets[ir] = *tmp;
+    }
+  // If we set a second branch for the input jets fetch this
+  TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
+  
+  if(!aodGenJets)
+    {
+      printf("NO MC jets Found\n");
+      return;
+    }
+  
+  //   //Generated jets
+  nGenJets = aodGenJets->GetEntries();
+  nGenJets = TMath::Min(nGenJets, kMaxJets);
+  
+  for(Int_t ig =0 ; ig < nGenJets; ++ig)
+    {
+      AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
+      if(!tmp)continue;
+      genJets[ig] = * tmp;
+    }
+  
+  AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
+  if(!pythiaGenHeader){
+    Printf("!!!NO GEN HEADER AVALABLE!!!");
+    return;
+  }
+  
+  // Int_t ProcessType = pythiaGenHeader->ProcessType();
+  // if(ProcessType != 28) return;
+  fhXSec->Fill(pythiaGenHeader->GetPtHard(), fXsection);
+  nPythiaGenJets = pythiaGenHeader->NTriggerJets();
+  nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
+   
+  Double_t eRec[kMaxJets];
+  Double_t eGen[kMaxJets];
+  Double_t eJetRec[kMaxJets];
+  //  Double_t EJetGen[kMaxJets];
+
+  AliAODJet jetRec[kMaxJets];
+  AliAODJet jetGen[kMaxJets];
+
+  Int_t idxRec[kMaxJets];
+  Int_t idxGen[kMaxJets];
+   
+  Double_t xRec[kMaxJets]; 
+  Double_t xGen[kMaxJets];
+
+  Double_t eSumRec = 0;
+  Double_t eSumGen = 0;
+  
+  TLorentzVector vRec[kMaxJets];
+  TLorentzVector vRestRec[kMaxJets];
+  
+  TLorentzVector vGen[kMaxJets];
+  TLorentzVector vRestGen[kMaxJets];
+
+  TLorentzVector vsumRec;
+  TLorentzVector vsumGen;
+
+  TVector3 pRec[kMaxJets];
+  TVector3 pGen[kMaxJets];
+
+  TVector3 pTrack[kTracks];
+
+  TVector3 pRestRec[kMaxJets];
+  TVector3 pRestGen[kMaxJets];
+
+  Double_t psumRestRec = 0;
+  //  Double_t psumRestGen = 0;
+  //Pythia_________________________________________________________________________________________________________________
+
+  for(int ip = 0;ip < nPythiaGenJets;++ip)
+    {
+      if(ip>=kMaxJets)continue;
+      Float_t p[4];
+      pythiaGenHeader->TriggerJet(ip,p);
+      genJetsPythia[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
+    }
+  
+//_________________________________________________________________________________________________________________________
+//________histos for MC___________________________________________________________________________________________________________
+
+  Int_t nGenSel = 0;
+  Int_t counter = 0;
+  Int_t tag = 0;
+
+  AliAODJet selJets[kMaxJets];
+
+  for(Int_t i = 0; i < nGenJets; i++)
+    {
+      if(nGenJets == 1)
+       {
+         selJets[nGenSel] = genJets[i];
+         nGenSel++;
+       }
+      else
+       {
+         counter = 0;
+         tag = 0;
+         for(Int_t j = 0; j < nGenJets; j++)
+           {
+             if(i!=j)
+               {
+                 Double_t dRij = genJets[i].DeltaR(&genJets[j]);
+                 counter++;
+                 if(dRij > 2*fR) tag++;
+               }
+           }
+         if(counter!=0)
+           {
+             if(tag/counter == 1)
+               {
+                 selJets[nGenSel] = genJets[i];
+                 nGenSel++;
+               }
+           }
+       }
+    }
+
+  if(nGenSel == 0) return;
+
+  for (Int_t gj = 0; gj < nGenSel; gj++)
+    {
+      eGen[gj] = selJets[gj].E();
+    }
+  
+  TMath::Sort(nGenSel, eGen, idxGen);
+  for (Int_t ig = 0; ig < nGenSel; ig++)
+    {
+      jetGen[ig] = selJets[idxGen[ig]];
+    }
+  AliStack * stack = mcEvent->Stack();
+  Int_t nMCtracks = stack->GetNprimary();
+  Double_t * eTracksMC = new Double_t[kTracks];
+  Double_t pTracksMC[kTracks];
+  Int_t * idxTracksMC = new Int_t[kTracks];
+  TLorentzVector jetTracksMC[kTracks];
+  TLorentzVector jetTracksSortMC[kTracks];
+  TVector3 pTrackMC[kTracks];
+  TLorentzVector vTrackMCAll[kTracks];
+  Double_t pTrackMCAll[kTracks];
+  TLorentzVector vTrackMC[kTracks];
+  TVector3 pTrackMCBoost[kTracks];
+  Double_t eventShapes[4];
+
+  Int_t nAccTr = 0;
+  Int_t nInJet[kMaxJets];
+  TLorentzVector inJetPartV[kMaxJets][kTracks];
+  Int_t nAllTracksMC = 0;
+
+  for(Int_t iTrack = 0; iTrack < nMCtracks; iTrack++)
+    {
+      TParticle * part = (TParticle*)stack->Particle(iTrack);
+      if (!part) continue;
+      Double_t fEta = part->Eta();
+      if(TMath::Abs(fEta) > .9) continue;
+      if(!IsPrimChar(part, nMCtracks, 0)) continue;
+      vTrackMCAll[nAllTracksMC].SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
+      pTrackMCAll[nAllTracksMC] = part->Pt();
+      nAllTracksMC++;
+    }
+  if(nAllTracksMC == 0) return;
+  for(Int_t iJet = 0; iJet < nGenSel; iJet++)
+    {
+      Int_t nJetTracks = 0;
+      for(Int_t i = 0; i < nAllTracksMC; i++)
+       {
+         Double_t dPhi = (jetGen[iJet].Phi()-vTrackMCAll[i].Phi());
+         if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
+         if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
+         Double_t dEta = (jetGen[iJet].Eta()-vTrackMCAll[i].Eta());
+         Double_t deltaR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
+         if(deltaR < fR && vTrackMCAll[i].Pt() > 1.5)
+           {
+             jetTracksMC[nAccTr] = vTrackMCAll[i];
+             eTracksMC[nAccTr] = vTrackMCAll[i].E();
+             pTracksMC[nAccTr] = vTrackMCAll[i].Pt();
+             inJetPartV[iJet][nJetTracks].SetPxPyPzE(vTrackMCAll[i].Px(), vTrackMCAll[i].Py(), vTrackMCAll[i].Pz(),vTrackMCAll[i].E());
+             nAccTr++;
+             nJetTracks++;
+           }
+       }
+      nInJet[iJet] = nJetTracks;
+    }
+
+  if(nAccTr == 0) return;
+  Printf("*********** Number of Jets : %d ***************\n", nGenSel);  
+  Double_t pTav[kMaxJets];
+  for(Int_t i = 0; i < nGenSel; i++)
+    {
+      Double_t pTsum = 0;
+      Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
+      for(Int_t iT = 0; iT < nInJet[i]; iT++)
+       {
+         Double_t pt = inJetPartV[i][iT].Pt();
+         pTsum += pt;
+       }
+      pTav[i] = pTsum/nInJet[i]; 
+    }
+  
+  TMath::Sort(nAccTr, pTracksMC, idxTracksMC);
+  for(Int_t i = 0; i < nAccTr; i++)
+    {
+      jetTracksSortMC[i] = jetTracksMC[idxTracksMC[i]];
+      pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
+      vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
+    }
+
+  TVector3 n01MC = pTrackMC[0].Unit();
+  //Thrust calculation, iterative method
+  if(nGenSel > 1)
+    { 
+//       if(fGlobVar == 1)
+//     {
+      GetEventShapes(n01MC, pTrackMC, nAccTr, eventShapes);
+//     }
+      
+      Double_t s = eventShapes[1];
+      Double_t a = eventShapes[2];
+      Double_t c = eventShapes[3];
+
+      switch(nGenSel)
+       {
+       case 2:
+         {
+           fhAGen2->Fill(a);
+           fhSGen2->Fill(s);
+           fhCGen2->Fill(c);
+         }
+         break;
+       case 3:
+         {
+           fhAGen3->Fill(a);
+           fhSGen3->Fill(s);
+           fhCGen3->Fill(c);
+         }
+         break;
+       }
+      Double_t thrust01MC = eventShapes[0];
+      
+      switch(nGenSel)
+       {
+       case 2:
+         fhThrustGen2->Fill(thrust01MC, fXsection);
+         break;
+       case 3:
+         fhThrustGen3->Fill(thrust01MC, fXsection);
+         break;
+       }
+    }
+  
+  
+  //rest frame MC jets
+  for (Int_t i = 0; i < nGenSel; ++i)
+    {
+      vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
+      pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
+      vsumGen += vGen[i];
+    }
+  if(eventShapes[0] >0.8 )
+    {
+      for(Int_t i = 0; i < nGenSel; i++)
+       fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
+    }
+  if(eventShapes[0] <= 0.8)   
+    {
+      for(Int_t i = 0; i < nGenSel; i++)
+       fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
+    }
+      
+  Double_t fPxGen = vsumGen.Px();
+  Double_t fPyGen = vsumGen.Py();
+  Double_t fPzGen = vsumGen.Pz();
+  Double_t fEGen = vsumGen.E();
+
+  Double_t eRestGen[kMaxJets];
+  for (Int_t j = 0; j < nGenSel; j++)
+    {
+      vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
+      eRestGen[j] = vGen[j].E();
+    }
+
+  for (Int_t j = 0; j < nAccTr; j++)
+    {
+      vTrackMC[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
+      pTrackMCBoost[j].SetXYZ(vTrackMC[j].Px(),vTrackMC[j].Py(),vTrackMC[j].Pz());
+    }
+      
+  Int_t idxRestGen[kMaxJets];
+  TMath::Sort(nGenSel, eRestGen, idxRestGen);
+  for(Int_t j = 0; j < nGenSel; j++)
+    {
+      vRestGen[j] = vGen[idxRestGen[j]];
+      eSumGen += vRestGen[j].E();
+    }
+
+  if (nGenSel == 3)
+    {
+      //      if(nInJet[0] < 3 || nInJet[1] < 3 || nInJet[2] < 3) return;
+      //       if(pRestGen[1].DeltaPhi(pRestGen[2]) > 0.95 && pRestGen[1].DeltaPhi(pRestGen[2]) < 1.15)
+      //       {
+      
+      for(Int_t i = 0; i < nGenSel; i++)
+       {
+         xGen[i] = 2*vRestGen[i].E()/eSumGen;
+       }
+      
+      Printf("***************** Values of Dalitz variables are : %f, %f, %f ****************\n", xGen[0], xGen[1], xGen[2]);
+      
+      Printf("***************** fXSection = %f ******************\n", fXsection);
+      if(eSumGen <= 60)
+       fhX3X4Gen60->Fill(xGen[0], xGen[1], fXsection);
+      
+      if(eSumGen > 60 && eSumGen <= 100)
+       fhX3X4Gen60100->Fill(xGen[0], xGen[1], fXsection);
+      
+      if(eSumGen > 100)
+       fhX3X4Gen100->Fill(xGen[0], xGen[1], fXsection);
+         
+      FillTopology(fhX3X4Gen, fhMu34Gen, fhMu45Gen, fhMu35Gen, xGen, pRestGen, fXsection);
+    }
+  
+
+
+//_______________________________________________histos for MC_____________________________________________________
+
+
+//_______________________________________________histos AOD________________________________________________________
+
+//   Printf("Event Number : %d, Number of gen jets : %d ", fEntry, nGenJets);
+
+  Int_t nRecSel = 0;
+  Int_t counter1 = 0;
+  Int_t tag1 = 0;
+
+  AliAODJet recSelJets[kMaxJets];
+
+  for(Int_t i = 0; i < nRecJets; i++)
+    {
+      if(nRecJets == 1)
+       {
+         recSelJets[nRecSel] = recJets[i];
+         nRecSel++;
+       }
+      else
+       {
+         counter1 = 0;
+         tag1 = 0;
+         for(Int_t j = 0; j < nRecJets; j++)
+           {
+             if(i!=j)
+               {
+                 Double_t dRij = recJets[i].DeltaR(&recJets[j]);
+                 counter1++;
+                 if(dRij > 2*fR) tag1++;
+               }
+           }
+         if(counter1!=0)
+           {
+             if(tag1/counter1 == 1)
+               {
+                 recSelJets[nRecSel] = recJets[i];
+                 nRecSel++;
+               }
+           }
+       }
+    } 
+  
+  if(nRecSel == 0) return;
+  
+  //sort rec/gen jets by energy in C.M.S
+  for (Int_t rj = 0; rj < nRecSel; rj++)
+    {
+      eRec[rj] = recSelJets[rj].E();
+    }
+  //  Int_t nAODtracks = fAOD->GetNumberOfTracks();
+  Int_t nTracks = 0; //tracks accepted in the whole event
+  //  Int_t nTracksALL = 0;
+  TLorentzVector jetTracks[kTracks];
+  TLorentzVector jetTracksSort[kTracks];
+  Double_t * eTracks = new Double_t[kTracks];
+  Double_t pTracks[kTracks];
+  Int_t * idxTracks = new Int_t[kTracks];
+  Double_t eventShapesRec[4];
+  Int_t jetMult[kMaxJets];
+  //  TLorentzVector vTracksAll[kTracks];
+  //  Double_t pTracksAll[kTracks];
+  Int_t nAccJets = 0;
+  AliAODJet jetRecAcc[kMaxJets];
+  Int_t nJetTracks = 0;
+
+  AliAODTrack jetTrack[kTracks];
+  Double_t * cv = new Double_t[21];
+  TMath::Sort(nRecSel, eRec, idxRec);
+  for (Int_t rj = 0; rj < nRecSel; rj++)
+    {
+      nJetTracks = 0;
+      eJetRec[rj] = eRec[idxRec[rj]];
+      jetRec[rj] = recSelJets[idxRec[rj]];
+      TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(jetRec[rj].GetRefTracks());
+      if(!jetTracksAOD) continue;
+      if(jetTracksAOD->GetEntries() < 3) continue;
+      for(Int_t i = 0; i < jetTracksAOD->GetEntries(); i++)
+       {
+         AliAODTrack * track = (AliAODTrack*)jetTracksAOD->At(i);
+         track->GetCovarianceXYZPxPyPz(cv);
+         if(cv[14] > 1000.) continue;
+         jetTrack[nTracks] = *track;
+         jetTracks[nTracks].SetPxPyPzE(jetTrack[nTracks].Px(), jetTrack[nTracks].Py(), jetTrack[nTracks].Pz(), jetTrack[nTracks].E());
+         eTracks[nTracks] = jetTracks[nTracks].E();
+         pTracks[nTracks] = jetTracks[nTracks].Pt();
+         nTracks++;
+         nJetTracks++;
+       }
+      if(nJetTracks < 3) continue;
+      jetRecAcc[nAccJets] = jetRec[rj];
+      jetMult[nAccJets] = jetTracksAOD->GetEntries();
+      nAccJets++;
+    }
+  
+  if (nAccJets == 0) return;
+
+  TLorentzVector vTrack[kTracks];
+  TMath::Sort(nTracks, pTracks, idxTracks);
+  for(Int_t i = 0; i < nTracks; i++)
+    {
+      jetTracksSort[i] = jetTracks[idxTracks[i]];
+      pTrack[i].SetXYZ(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz());
+      vTrack[i].SetPxPyPzE(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz(), jetTracksSort[i].E());
+    }
+
+  for (Int_t i = 0; i < nAccJets; ++i)
+    {
+      vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
+      pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
+      vsumRec += vRec[i];
+    }
+
+  //Thrust, iterative method, AODs
+  TVector3 n01 = pTrack[0].Unit();
+  if(nAccJets > 1)
+    {
+//       if(fGlobVar == 1)
+//     {
+      GetEventShapes(n01, pTrack, nTracks, eventShapesRec);
+//     }
+//       fGlobVar = 0;      
+//       Double_t Max3 = TMath::Max(eventShapesRec0[0],eventShapesRec1[0]);
+//       Double_t Max4 = TMath::Max(eventShapesRec3[0],eventShapesRec2[0]);
+      Double_t thrust = eventShapesRec[0];
+      
+      if(eventShapesRec[0] > 0.8)
+       {
+         for(Int_t i = 0; i < nAccJets; i++)
+           fhdPhiThrustRec->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
+         
+       }
+      
+      if(eventShapesRec[0] <= 0.8)
+       {
+         for(Int_t i = 0; i < nAccJets; i++)
+           fhdPhiThrustRecALL->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
+       }
+
+      switch(nAccJets)
+       {
+       case 2:
+         fhThrustRec2->Fill(thrust, fXsection);
+         break;
+       case 3:
+         fhThrustRec3->Fill(thrust, fXsection);
+         break;
+       }
+      
+      switch(nAccJets)
+       {
+       case 2:
+         {
+           fhARec2->Fill(eventShapesRec[2], fXsection);
+           fhSRec2->Fill(eventShapesRec[1], fXsection);
+           fhCRec2->Fill(eventShapesRec[3], fXsection);
+         }
+         break;
+       case 3:
+         {
+           fhARec3->Fill(eventShapesRec[2], fXsection);
+           fhSRec3->Fill(eventShapesRec[1], fXsection);
+           fhCRec3->Fill(eventShapesRec[3], fXsection);
+         }
+         break;
+       }
+      
+    }
+
+  //rest frame for reconstructed jets
+  Double_t fPx = vsumRec.Px();
+  Double_t fPy = vsumRec.Py();
+  Double_t fPz = vsumRec.Pz();
+  Double_t fE = vsumRec.E();
+                 
+  TVector3 pTrackRest[kTracks];
+  for(Int_t j = 0; j < nTracks; j++)
+    {
+      vTrack[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
+      pTrackRest[j].SetXYZ(vTrack[j].Px(), vTrack[j].Py(),vTrack[j].Pz());
+    }
+      
+  Double_t eRestRec[kMaxJets];
+  Int_t idxRestRec[kMaxJets];
+  for (Int_t j = 0; j < nAccJets; j++)
+    {
+      vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
+      eRestRec[j] = vRec[j].E();  
+      eSumRec += vRec[j].E();
+    }
+
+  TMath::Sort(nAccJets, eRestRec, idxRestRec);
+  for (Int_t i = 0; i < nAccJets; i++)
+    {
+      vRestRec[i] = vRec[idxRestRec[i]]; 
+      pRestRec[i].SetXYZ(vRestRec[i].Px(), vRestRec[i].Py(), vRestRec[i].Pz());
+      psumRestRec += pRestRec[i].Perp();
+    } 
+
+  if(nAccJets == 3)
+    {
+//       if(pRest[1].DeltaPhi(pRest[2]) > 0.95 && pRest[1].DeltaPhi(pRest[2]) < 1.15)
+//     {
+         fhInOut->Fill(nGenSel);
+//       for(Int_t j = 0; j < nTracksALL; j++)
+//         {
+//           vTracksAll[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
+//           pTracksAll[j].SetXYZ(vTracksAll[j].Px(), vTracksAll[j].Py(),vTracksAll[j].Pz());
+//           fhdPhiRec->Fill(pRest[0].DeltaPhi(pTracksAll[j]), pTracksAll[j].Perp(), fXsection);             
+//         }
+         //and the Dalitz variables and Energy distributions in the rest frame
+         for (Int_t i = 0; i < nAccJets; i++)
+           xRec[i] = 2*vRestRec[i].E()/eSumRec;
+         
+         if(eSumRec <= 60)
+           fhX3X4Rec60->Fill(xRec[0], xRec[1], fXsection);
+         
+         if(eSumRec > 60 && eSumRec <= 100)
+           fhX3X4Rec60100->Fill(xRec[0], xRec[1], fXsection);
+         
+         if(eSumRec > 100)
+           fhX3X4Rec100->Fill(xRec[0], xRec[1], fXsection);
+         
+         if(nAccJets == 3 && nAccJets == nGenJets)
+           {
+             fhX3->Fill(xGen[0], TMath::Abs(xGen[0]-xRec[0])/xGen[0], fXsection);
+             fhX4->Fill(xGen[1], TMath::Abs(xGen[1]-xRec[1])/xGen[1], fXsection);
+             fhX5->Fill(xGen[2], TMath::Abs(xGen[2]-xRec[2])/xGen[2], fXsection);
+           }
+         
+         FillTopology(fhX3X4Rec, fhMu34Rec, fhMu45Rec, fhMu35Rec, xRec, pRestRec, fXsection);
+    }
+  Printf("%s:%d",(char*)__FILE__,__LINE__);
+  
+  PostData(1, fList);
+
+  Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
+  
+}
+
+//__________________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskThreeJets::Terminate(Option_t *)
+{
+  printf("AnalysisJetCorrelation::Terminate()");
+
+} 
+
+//_______________________________________User defined functions_____________________________________________________________________________________
+void AliAnalysisTaskThreeJets::FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection)
+{
+  Dalitz->Fill(x[0], x[1], xsection);
+  fhMu35->Fill(TMath::Sqrt(x[0]*x[2]*(1-(pRest[0].Unit()).Dot(pRest[2].Unit()))/2), xsection);
+  fhMu34->Fill(TMath::Sqrt(x[0]*x[1]*(1-(pRest[0].Unit()).Dot(pRest[1].Unit()))/2), xsection);
+  fhMu45->Fill(TMath::Sqrt(x[1]*x[2]*(1-(pRest[1].Unit()).Dot(pRest[2].Unit()))/2), xsection);
+}
+
+//_____________________________________________________________________________________________________________________________
+
+Bool_t AliAnalysisTaskThreeJets::IsPrimChar(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
+{
+  //
+  // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
+  // shall be counted as a primary particle
+  //
+  // This function or a equivalent should be available in some common place of AliRoot
+  //
+  // WARNING: Call this function only for particles that are among the particles from the event generator!
+  // --> stack->Particle(id) with id < stack->GetNprimary()
+
+  // if the particle has a daughter primary, we do not want to count it
+  if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
+  {
+    if (adebug)
+      printf("Dropping particle because it has a daughter among the primaries.\n");
+    return kFALSE;
+  }
+
+  Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
+  
+
+  // skip quarks and gluon
+  if (pdgCode <= 10 || pdgCode == 21)
+  {
+    if (adebug)
+      printf("Dropping particle because it is a quark or gluon.\n");
+    return kFALSE;
+  }
+
+  Int_t status = aParticle->GetStatusCode();
+  // skip non final state particles..
+  if(status!=1){
+    if (adebug)
+      printf("Dropping particle because it is not a final state particle.\n");
+    return kFALSE;
+  }
+
+  if (strcmp(aParticle->GetName(),"XXX") == 0)
+  {
+    Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
+    return kFALSE;
+  }
+
+  TParticlePDG* pdgPart = aParticle->GetPDG();
+
+  if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
+  {
+    Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
+    return kFALSE;
+  }
+
+  if (pdgPart->Charge() == 0)
+  {
+    if (adebug)
+      printf("Dropping particle because it is not charged.\n");
+    return kFALSE;
+  }
+
+  return kTRUE;
+}
+
+//______________________________________________________________________________________________________
+
+void AliAnalysisTaskThreeJets::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, Int_t &nTracks)
+{
+  TVector3 psum;
+  Double_t psum1 = 0;
+  Double_t psum2 = 0;
+  Double_t thrust[kTracks];
+  Double_t th = -3;
+  
+  for(Int_t j = 0; j < nTracks; j++)
+    {
+      psum.SetXYZ(0., 0., 0.);
+      psum1 = 0;
+      psum2 = 0;
+      for(Int_t i = 0; i < nTracks; i++)
+       {
+         psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
+         psum2 += pTrack[i].Mag();
+         
+         if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
+         if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
+       }
+      
+      thrust[j] = psum1/psum2;
+      
+      if(th == thrust[j]) 
+       break;
+      
+      th = thrust[j];
+      
+      n01 = psum.Unit();
+    }
+}
+
+//___________________________________________________________________________________________________________
+
+void AliAnalysisTaskThreeJets::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
+{       
+  TVector3 psum;
+  Double_t psum1 = 0;
+  Double_t psum2 = 0;
+  Double_t thrust[kTracks];
+  Double_t th = -3;
+  
+  //Sphericity calculation
+  TMatrixDSym m(3);
+  Double_t s00 = 0;
+  Double_t s01 = 0;
+  Double_t s02 = 0;
+  
+  Double_t s10 = 0;
+  Double_t s11 = 0;
+  Double_t s12 = 0;
+  
+  Double_t s20 = 0;
+  Double_t s21 = 0;
+  Double_t s22 = 0;
+  
+  Double_t ptot = 0;
+  
+  Double_t c = -10;
+  
+  for(Int_t j = 0; j < nTracks; j++)
+    {  
+      psum.SetXYZ(0., 0., 0.);
+      psum1 = 0;
+      psum2 = 0;
+      for(Int_t i = 0; i < nTracks; i++)
+       {
+         psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
+         psum2 += pTrack[i].Mag();
+         
+         if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
+         if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
+       }
+      
+      thrust[j] = psum1/psum2;
+      if(th == thrust[j]) 
+       break;
+      
+      th = thrust[j];
+      
+      n01 = psum.Unit();
+    }
+
+  eventShapes[0] = th;
+
+  for(Int_t j = 0; j < nTracks; j++)
+    {  
+      s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
+      s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
+      s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
+      
+      s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
+      s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
+      s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
+      
+      s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
+      s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
+      s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
+      
+      ptot += pTrack[j].Mag();
+    }
+
+  if(ptot !=0)
+    {
+      m(0,0) = s00/ptot;
+      m(0,1) = s01/ptot;
+      m(0,2) = s02/ptot;
+
+      m(1,0) = s10/ptot;
+      m(1,1) = s11/ptot;
+      m(1,2) = s12/ptot;
+
+      m(2,0) = s20/ptot;
+      m(2,1) = s21/ptot;
+      m(2,2) = s22/ptot;
+
+      TMatrixDSymEigen eigen(m);
+      TVectorD eigenVal = eigen.GetEigenValues();
+
+      Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
+      eventShapes[1] = sphericity;
+
+      Double_t aplaarity = (3/2)*(eigenVal(2));
+      eventShapes[2] = aplaarity;
+
+      c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
+      eventShapes[3] = c;
+    }
+}
+  
+//__________________________________________________________________________________________________________________________
+
+Double_t AliAnalysisTaskThreeJets::Exponent(Double_t x, Double_t * par)
+{
+  return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2])+0.5*TMath::Power((x-par[3])/par[0], 2))+par[4]*x;
+}
+
+Double_t AliAnalysisTaskThreeJets::Exponent2(Double_t x, Double_t * par)
+{
+  return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2]))+par[3]*x;
+}
+
+Double_t AliAnalysisTaskThreeJets::Gauss(Double_t x, Double_t * par)
+{
+  return 1/(par[1])*TMath::Power(1/TMath::E(), 0.5*(x-par[0])*(x-par[0])/(par[1]*par[1]));
+}
+
+Double_t AliAnalysisTaskThreeJets::Total(Double_t x, Double_t * par)
+{
+  return Exponent(x, par)+Gauss(x, &par[4]);
+}
+
+  
diff --git a/PWG4/JetTasks/AliAnalysisTaskThreeJets.h b/PWG4/JetTasks/AliAnalysisTaskThreeJets.h
new file mode 100644 (file)
index 0000000..c1ac94a
--- /dev/null
@@ -0,0 +1,146 @@
+#ifndef ALIANALYSISTASKTHREEJETS_H\r
+#define ALIANALYSISTASKTHREEJETS_H\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include <iostream>\r
+#include <fstream>\r
+\r
+//class AliJetHistos;\r
+//class AliJetCorrHistos;\r
+class AliJetFinder;\r
+class AliESDEvent;\r
+class AliAODEvent;\r
+class AliAODJet;\r
+class AliGenPythiaEventHeader;\r
+class AliAODPid;\r
+\r
+class TList;\r
+class TArrayD;\r
+class TChain;\r
+class TH1;\r
+class TH2;\r
+class TH1F;\r
+class TH2F;\r
+class TH2I;\r
+class TH3D;\r
+class TTree;\r
+class TProfile;\r
+class TLorentzVector;\r
+class TVector3;\r
+class TParticle;\r
+\r
+class AliAnalysisTaskThreeJets : public AliAnalysisTaskSE\r
+{\r
+ public:\r
+  AliAnalysisTaskThreeJets();\r
+  AliAnalysisTaskThreeJets(const char * name);\r
+  virtual ~AliAnalysisTaskThreeJets() {;}\r
+  \r
+  //Implementation of interface methods\r
+  virtual Bool_t Notify(); \r
+  virtual void UserCreateOutputObjects();\r
+  virtual void Init();\r
+  virtual void LocalInit() { Init(); };\r
+  virtual void UserExec(Option_t * option);\r
+  virtual void Terminate(Option_t * option);\r
+\r
+  virtual void SetAODInput(Bool_t b){fUseAODInput = b;}\r
+\r
+  virtual void SetBranchGen(char* c){fBranchGen = c;}\r
+  virtual void SetBranchRec(char* c){fBranchRec = c;}\r
+\r
+/*   virtual void SetUseBkg(Bool_t b) {fUseBkg = b;} */\r
+/*   virtual void SetUseSimTPC(Bool_t b) {fUseFastTPC = b;} */\r
+\r
+  virtual void FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection);\r
+\r
+  virtual Double_t SetR(Double_t b){fR = b; return fR;} \r
+//  TArrayD * GetThrustParamMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom, TArrayD * evsh);\r
+  virtual void GetThrustAxis(TVector3 &n01, TVector3 * p_track, Int_t &nTracks);\r
+  virtual void GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes);\r
+  virtual Bool_t IsPrimChar(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug);\r
+\r
+  Double_t Exponent(Double_t x, Double_t * par);\r
+  Double_t Exponent2(Double_t x, Double_t * par);\r
+  Double_t Gauss(Double_t x, Double_t * par);\r
+  Double_t Total(Double_t x, Double_t * par);\r
+\r
+ private:\r
+  AliAnalysisTaskThreeJets(const AliAnalysisTaskThreeJets&);\r
+  AliAnalysisTaskThreeJets& operator = (const AliAnalysisTaskThreeJets&);\r
+\r
+  enum {kMaxJets = 6};\r
+  enum {kMaxEvents = 10};\r
+  enum {kJets = 3};\r
+  enum {kTracks = 1000};\r
+\r
+  AliAODEvent  *fAOD; // where we take the jets from can be input or output AOD\r
+  \r
+  TString       fBranchRec;  // AOD branch name for reconstructed\r
+  TString       fBranchGen;  // AOD brnach for genereated\r
+  \r
+  Bool_t        fUseAODInput;\r
+\r
+  Double_t fR;\r
+  TList * fList;\r
+\r
+  Int_t fGlobVar;\r
+  Double_t fXsection;\r
+\r
+  TH2F * fhX3X4Rec;\r
+  TH2F * fhX3X4Gen;\r
+  \r
+  TH1F * fhMu35Rec; \r
+  TH1F * fhMu34Rec; \r
+  TH1F * fhMu45Rec; \r
+  \r
+  TH1F * fhMu35Gen; \r
+  TH1F * fhMu34Gen; \r
+  TH1F * fhMu45Gen; \r
+\r
+  TH1I * fhInOut;\r
+  TH1F * fhThrustRec2;\r
+  TH1F * fhThrustRec3;\r
+\r
+  TH1F * fhThrustGen2;\r
+  TH1F * fhThrustGen3;\r
+\r
+  TH1F * fhCGen2;\r
+  TH1F * fhCGen3;\r
+\r
+  TH1F * fhSGen2;\r
+  TH1F * fhSGen3;\r
+\r
+  TH1F * fhAGen2;\r
+  TH1F * fhAGen3;\r
+\r
+  TH1F * fhCRec2;\r
+  TH1F * fhCRec3;\r
+\r
+  TH1F * fhSRec2;\r
+  TH1F * fhSRec3;\r
+\r
+  TH1F * fhARec2;\r
+  TH1F * fhARec3;\r
+\r
+  TH2F * fhX3;\r
+  TH2F * fhX4;\r
+  TH2F * fhX5;\r
+\r
+  TProfile * fhXSec;\r
+  TH2F * fhX3X4Rec60;\r
+  TH2F * fhX3X4Rec60100;\r
+  TH2F * fhX3X4Rec100;\r
+  TH2F * fhX3X4Gen60;\r
+  TH2F * fhX3X4Gen60100;\r
+  TH2F * fhX3X4Gen100;\r
+\r
+  TH2F * fhdPhiThrustGen;\r
+  TH2F * fhdPhiThrustGenALL;\r
+  TH2F * fhdPhiThrustRec;\r
+  TH2F * fhdPhiThrustRecALL;\r
+\r
+  ClassDef(AliAnalysisTaskThreeJets, 1)\r
+};\r
+\r
+#endif\r
index 272de46..e52a419 100644 (file)
@@ -8,6 +8,8 @@
 #pragma link C++ class AliAnalysisTaskJetSpectrum+;
 #pragma link C++ class AliAnalysisTaskJFSystematics+;
 #pragma link C++ class AliAnalysisHelperJetTasks+;
+#pragma link C++ class AliAnalysisTaskJetCorrections+;
+#pragma link C++ class AliAnalysisTaskThreeJets+;
 #pragma link C++ class AliAnaESDSpectraQA+;
 #pragma link C++ class AliAnalysisTaskPWG4PidDetEx+;
 #pragma link C++ class AliJetSpectrumUnfolding+;
index e09877d..5fb5da3 100644 (file)
@@ -1,6 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx
+SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx JetTasks/AliAnalysisTaskJetCorrections.cxx JetTasks/AliAnalysisTaskThreeJets.cxx
 
 HDRS:= $(SRCS:.cxx=.h) 
 
diff --git a/PWG4/macros/AddTaskJetCorrections.C b/PWG4/macros/AddTaskJetCorrections.C
new file mode 100644 (file)
index 0000000..bf0e4f1
--- /dev/null
@@ -0,0 +1,63 @@
+AliAnalysisTaskJetCorrections * AddTaskJetCorrections()\r
+{\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskJetSpectrum", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }  \r
+  \r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //==============================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskJetSpectrum", "This task requires an input event handler");\r
+    return NULL;\r
+   }\r
+  \r
+  AliAnalysisTaskJetCorrections * jetCorr = new  AliAnalysisTaskJetCorrections("Jet Corrections");\r
+  \r
+  jetCorr->SetBranchGen("jetsMC"); \r
+  jetCorr->SetBranchRec("jets");\r
+  jetCorr->SetR(.5); \r
+  mgr->AddTask(jetCorr);\r
+   \r
+  AliAnalysisDataContainer *coutput1_Corr = mgr->CreateContainer("jetCorr", TList::Class(),AliAnalysisManager::kOutputContainer,"jetCorr.root");\r
+\r
+   mgr->ConnectInput  (jetCorr, 0, mgr->GetCommonInputContainer());\r
+   mgr->ConnectOutput (jetCorr, 0, mgr->GetCommonOutputContainer());\r
+   mgr->ConnectOutput (jetCorr,  1, coutput1_Corr );\r
+   \r
+   return jetCorr;\r
+}\r
+\r
+\r
+AliAnalysisTaskJetCorrections * AddTaskJetCorrections(AliAnalysisManager* mgr ,AliAnalysisDataContainer *cinput)\r
+{\r
+  if (!mgr) {\r
+    ::Error("AddTaskJetSpectrum", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }  \r
+  \r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //==============================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskJetSpectrum", "This task requires an input event handler");\r
+    return NULL;\r
+   }\r
+  \r
+  AliAnalysisTaskJetCorrections * jetCorr = new  AliAnalysisTaskJetCorrections("Jet Corrections");\r
+  \r
+  jetCorr->SetBranchGen("jetsMC"); \r
+  jetCorr->SetBranchRec("jets");\r
+  jetCorr->SetR(.5); \r
+  mgr->AddTask(jetCorr);\r
+\r
+  AliAnalysisDataContainer * coutpu0 = mgr->CreateContainer("coutpu0", TTree::Class(),\r
+                                 AliAnalysisManager::kExchangeContainer);\r
+  AliAnalysisDataContainer *coutput1_jetCorr = mgr->CreateContainer("jetCorr", TList::Class(),AliAnalysisManager::kOutputContainer,"jetCorr.root");\r
+\r
+   mgr->ConnectInput  (jetCorr, 0, cinput);\r
+   mgr->ConnectOutput (jetCorr, 0, coutpu0);\r
+   mgr->ConnectOutput (jetCorr,  1, coutput1_Corr );\r
+   \r
+   return jetCorr;\r
+}\r
diff --git a/PWG4/macros/AddTaskThreeJets.C b/PWG4/macros/AddTaskThreeJets.C
new file mode 100644 (file)
index 0000000..94019a1
--- /dev/null
@@ -0,0 +1,73 @@
+AliAnalysisTaskThreeJets * AddTaskThreeJets()\r
+{\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskThreeJets", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }  \r
+  \r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //==============================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskThreeJets", "This task requires an input event handler");\r
+    return NULL;\r
+   }\r
+  \r
+  // Create the task and configure it.\r
+  //===========================================================================\r
+  \r
+  AliAnalysisTaskThreeJets * threeJets = new  AliAnalysisTaskThreeJets("Three Jet Analysis");\r
+  \r
+  threeJets->SetBranchGen("jetsMC"); \r
+  threeJets->SetBranchRec("jets");\r
+  threeJets->SetR(.5); \r
+  mgr->AddTask(threeJets);\r
+   \r
+\r
+\r
+      \r
+   // Create ONLY the output containers for the data produced by the task.\r
+   // Get and connect other common input/output containers via the manager as below\r
+   //==============================================================================\r
+   AliAnalysisDataContainer *coutput1_Corr = mgr->CreateContainer("threeJets", TList::Class(),AliAnalysisManager::kOutputContainer,"threeJets.root");\r
+\r
+   mgr->ConnectInput  (threeJets, 0, mgr->GetCommonInputContainer());\r
+   mgr->ConnectOutput (threeJets, 0, mgr->GetCommonOutputContainer());\r
+   mgr->ConnectOutput (threeJets,  1, coutput1_Corr );\r
+   \r
+   return threeJets;\r
+}\r
+\r
+\r
+AliAnalysisTaskThreeJets * AddTaskJetCorrections(AliAnalysisManager* mgr,AliAnalysisDataContainer *cinput)\r
+{\r
+  if (!mgr) {\r
+    ::Error("AddTaskJetSpectrum", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }  \r
+  \r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //==============================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskJetSpectrum", "This task requires an input event handler");\r
+    return NULL;\r
+   }\r
+  \r
+  AliAnalysisTaskThreeJets * threeJets = new  AliAnalysisTaskThreeJets("ThreeJetAnalysis");\r
+  \r
+  threeJets->SetBranchGen("jetsMC"); \r
+  threeJets->SetBranchRec("jets");\r
+  threeJets->SetR(.5); \r
+  mgr->AddTask(threeJets);\r
+\r
+  AliAnalysisDataContainer * coutpu0 = mgr->CreateContainer("coutpu0", TTree::Class(),\r
+                                 AliAnalysisManager::kExchangeContainer);\r
+  AliAnalysisDataContainer *coutput1_threeJets = mgr->CreateContainer("threeJets", TList::Class(),AliAnalysisManager::kOutputContainer,"threeJets.root");\r
+\r
+   mgr->ConnectInput  (threeJets, 0, cinput);\r
+   mgr->ConnectOutput (threeJets, 0, coutpu0);\r
+   mgr->ConnectOutput (threeJets,  1, coutput1_Corr );\r
+   \r
+   return threeJets;\r
+}\r
+\r
index 6c95390..534dcfc 100644 (file)
@@ -10,7 +10,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
   Bool_t debug         = kTRUE;\r
   Bool_t useMC         = kTRUE;\r
   Bool_t readTR        = kFALSE;\r
-  Bool_t bPROOF        = kFALSE;\r
+  Bool_t bPROOF        = kTRUE;\r
   Bool_t bLOCALPAR     = kFALSE;  // flag that swtiches on loading of local par files insead of loading libs, needed for grid and local testing\r
 \r
     \r
@@ -24,6 +24,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
   Int_t iDIJETAN       = 1;\r
   Int_t iPWG4SPECTRUM  = 1;\r
   Int_t iPWG4JFSYSTEMATICS  = 1;\r
+  Int_t iPWG4JETCORRECTION  = 1;\r
+  Int_t iPWG4THREEJETS  = 1;\r
   Int_t iPWG4UE        = 0;\r
   Int_t iPWG4PID        = 0;\r
 \r
@@ -40,8 +42,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
   TString dataset(ds);\r
   TChain *chain = 0;\r
   // CKB quick hack for local analysis\r
-  gROOT->LoadMacro("CreateESDChain.C");\r
-  TChain *chain = CreateESDChain("tmp.txt",1000);\r
+  // gROOT->LoadMacro("CreateESDChain.C");\r
+  // TChain *chain = CreateESDChain("tmp.txt",1000);\r
 \r
  \r
   printf("==================================================================\n");\r
@@ -57,6 +59,9 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
   if (iDIJETAN)     printf("=  DiJet analysis                                                 =\n");\r
   if (iPWG4SPECTRUM)printf("=  PWG4 Jet spectrum analysis                                    =\n");\r
   if (iPWG4JFSYSTEMATICS)printf("=  PWG4 Jet Finder systematics                                   =\n");\r
+  if (iPWG4JETCORRECTION)printf("=  PWG4 Jet Correction                                   =\n");\r
+  if (iPWG4THREEJETS)printf("=  PWG4 Three Jets                                   =\n");\r
+\r
   if (iPWG4UE)      printf("=  PWG4 UE                                                        =\n");\r
   printf("==================================================================\n");\r
   if (useMC) printf(":: use MC    TRUE\n");\r
@@ -77,8 +82,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
   // TProof::Reset("alicecaf"); \r
   // One may enable a different ROOT version on CAF\r
   \r
-  const char* proofNode = "localhost";\r
-  //  const char* proofNode = "alicecaf";\r
+  //  const char* proofNode = "localhost";\r
+  const char* proofNode = "alicecaf";\r
   \r
   // Connect to proof\r
   if(bPROOF){\r
@@ -116,7 +121,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
       gProof->EnablePackage("JETAN");\r
     }   \r
     // --- Enable particle correlation analysis\r
-    if (iPWG4UE||iPWG4SPECTRUM||iPWGJFSYSTEMATICS) {\r
+    if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS||iPWG4JETCORRECTION||iPWG4THREEJETS) {\r
       gProof->UploadPackage("PWG4JetTasks.par");\r
       gProof->EnablePackage("PWG4JetTasks");\r
     }   \r
@@ -154,7 +159,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
       // --- Enable the JETAN Package\r
       if (iJETAN||iJETANESD||iJETANMC||iJETANMC2) gSystem->Load("libJETAN");\r
       // --- Enable particle correlation analysis\r
-      if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS)gSystem->Load("libPWG4JetTasks"); \r
+      if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS||iPWG4THREEJETS)gSystem->Load("libPWG4JetTasks"); \r
     }\r
 \r
   }\r
@@ -245,6 +250,16 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG
       AliAnalysisTaskJFSystematics* pwg4jfs = AddTaskJFSystematics();\r
       pwg4jfs->SetDebugLevel(11);\r
     }   \r
+    if (iPWG4JETCORRECTION) {\r
+      gROOT->LoadMacro("AddTaskJetCorrections.C");\r
+      AliAnalysisTaskJetCorrections* pwg4jc = AddTaskJetCorrections();\r
+      pwg4jc->SetDebugLevel(11);\r
+    }   \r
+    if (iPWG4THREEJETS) {\r
+      gROOT->LoadMacro("AddTaskThreeJets.C");\r
+      AliAnalysisTaskThreeJets* pwg4jjj = AddTaskThreeJets();\r
+      pwg4jjj->SetDebugLevel(11);\r
+    }   \r
     if (iPWG4UE) {\r
       gROOT->LoadMacro("AddTaskUE.C");\r
       AliAnalysisTaskUE* ueana = AddTaskUE();\r