From 95392a578f466991a62d6a23648e3ba5ab3ebada Mon Sep 17 00:00:00 2001 From: kleinb Date: Sun, 9 Aug 2009 16:36:43 +0000 Subject: [PATCH] Added tasks from Sona --- PWG4/CMake_libPWG4JetTasks.txt | 2 + .../AliAnalysisTaskJetCorrections.cxx | 1002 +++++++++++++ PWG4/JetTasks/AliAnalysisTaskJetCorrections.h | 121 ++ PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx | 1333 +++++++++++++++++ PWG4/JetTasks/AliAnalysisTaskThreeJets.h | 146 ++ PWG4/PWG4JetTasksLinkDef.h | 2 + PWG4/libPWG4JetTasks.pkg | 2 +- PWG4/macros/AddTaskJetCorrections.C | 63 + PWG4/macros/AddTaskThreeJets.C | 73 + PWG4/macros/AnalysisTrainCAF.C | 29 +- 10 files changed, 2765 insertions(+), 8 deletions(-) create mode 100644 PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx create mode 100644 PWG4/JetTasks/AliAnalysisTaskJetCorrections.h create mode 100644 PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx create mode 100644 PWG4/JetTasks/AliAnalysisTaskThreeJets.h create mode 100644 PWG4/macros/AddTaskJetCorrections.C create mode 100644 PWG4/macros/AddTaskThreeJets.C diff --git a/PWG4/CMake_libPWG4JetTasks.txt b/PWG4/CMake_libPWG4JetTasks.txt index c730b3c6391..060a636c5cd 100644 --- a/PWG4/CMake_libPWG4JetTasks.txt +++ b/PWG4/CMake_libPWG4JetTasks.txt @@ -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 index 00000000000..e2e1dd06bcf --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx @@ -0,0 +1,1002 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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(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(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(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(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(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(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(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(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(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=0&&ig 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 index 00000000000..2fd2f4dcb14 --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskJetCorrections.h @@ -0,0 +1,121 @@ +#ifndef ALIANALYSISTASKJETCORRECTIONS_H +#define ALIANALYSISTASKJETCORRECTIONS_H + +#include "AliAnalysisTaskSE.h" +#include +#include + +class AliJetFinder; +class AliESDEvent; +class AliAODEvent; +class AliAODJet; +class AliGenPythiaEventHeader; +class AliAODPid; + +class TList; +class TArrayD; +class TChain; +class TH1; +class TH2; +class TH1F; +class TH2F; +class TH2I; +class TH3D; +class TTree; +class TProfile; +class TLorentzVector; +class TVector3; +class TParticle; + +class AliAnalysisTaskJetCorrections : public AliAnalysisTaskSE +{ + public: + AliAnalysisTaskJetCorrections(); + AliAnalysisTaskJetCorrections(const char * name); + virtual ~AliAnalysisTaskJetCorrections() {;} + + //Implementation of interface methods + virtual Bool_t Notify(); + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() { Init(); }; + virtual void UserExec(Option_t * option); + virtual void Terminate(Option_t * option); + + virtual void SetAODInput(Bool_t b){fUseAODInput = b;} + + virtual void SetBranchGen(char* c){fBranchGen = c;} + virtual void SetBranchRec(char* c){fBranchRec = c;} + + virtual Double_t SetR(Double_t b){fR = b; return fR;} + + virtual void GetThrustAxis(TVector3 &n01, TVector3 * pTrack, Int_t &nTracks); + private: + AliAnalysisTaskJetCorrections(const AliAnalysisTaskJetCorrections&); + AliAnalysisTaskJetCorrections& operator = (const AliAnalysisTaskJetCorrections&); + + enum {kMaxJets = 6}; + enum {kMaxEvents = 10}; + enum {kJets = 3}; + enum {kTracks = 1000}; + + AliAODEvent *fAOD; // where we take the jets from can be input or output AOD + + TString fBranchRec; // AOD branch name for reconstructed + TString fBranchGen; // AOD brnach for genereated + + Bool_t fUseAODInput; + + Double_t fR; + TList * fList; + + Int_t fGlobVar; + Double_t fXsection; + + + TH1F * fhEGen; + TH1F * fhERec; + TH1F * fhEGenRest; + TH1F * fhERecRest; + TH1F * fhEsumGenRest; + TH1F * fhEsumRecRest; + + TH2F * fhE2vsE1Gen; + TH2F * fhE2vsE1Rec; + TH2F * fhE2E1vsEsumGen; + TH2F * fhE2E1vsEsumRec; + TH2F * fhE2E1vsE1Gen; + TH2F * fhE2E1vsE1Rec; + TH2F * fhE2E1vsdPhiGen; + TH2F * fhE2E1vsdPhiRec; + + TH2F * fhTrackBalance2; + TH2F * fhTrackBalance3; + + TH2F * fhEt1Et22; + TH2F * fhEt1Et23; + + TProfile * fhECorrJet10[3]; + TProfile * fhECorrJet05[3]; + TProfile * fhECorrJet01[3]; + TProfile * fhECorrJet001[3]; + + TProfile * fhdEvsErec10[3]; + TProfile * fhdEvsErec05[3]; + TProfile * fhdEvsErec01[3]; + TProfile * fhdEvsErec001[3]; + + TH2F * fhdPhidEta10[3]; + TH2F * fhdPhidEta05[3]; + TH2F * fhdPhidEta01[3]; + TH2F * fhdPhidEta001[3]; + + TH2F * fhdPhidEtaPt10[3]; + TH2F * fhdPhidEtaPt05[3]; + TH2F * fhdPhidEtaPt01[3]; + TH2F * fhdPhidEtaPt001[3]; + + ClassDef(AliAnalysisTaskJetCorrections, 1) +}; + +#endif diff --git a/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx new file mode 100644 index 00000000000..96e856164f6 --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx @@ -0,0 +1,1333 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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(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(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(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(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(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(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(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(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 index 00000000000..c1ac94ac8c2 --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskThreeJets.h @@ -0,0 +1,146 @@ +#ifndef ALIANALYSISTASKTHREEJETS_H +#define ALIANALYSISTASKTHREEJETS_H + +#include "AliAnalysisTaskSE.h" +#include +#include + +//class AliJetHistos; +//class AliJetCorrHistos; +class AliJetFinder; +class AliESDEvent; +class AliAODEvent; +class AliAODJet; +class AliGenPythiaEventHeader; +class AliAODPid; + +class TList; +class TArrayD; +class TChain; +class TH1; +class TH2; +class TH1F; +class TH2F; +class TH2I; +class TH3D; +class TTree; +class TProfile; +class TLorentzVector; +class TVector3; +class TParticle; + +class AliAnalysisTaskThreeJets : public AliAnalysisTaskSE +{ + public: + AliAnalysisTaskThreeJets(); + AliAnalysisTaskThreeJets(const char * name); + virtual ~AliAnalysisTaskThreeJets() {;} + + //Implementation of interface methods + virtual Bool_t Notify(); + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() { Init(); }; + virtual void UserExec(Option_t * option); + virtual void Terminate(Option_t * option); + + virtual void SetAODInput(Bool_t b){fUseAODInput = b;} + + virtual void SetBranchGen(char* c){fBranchGen = c;} + virtual void SetBranchRec(char* c){fBranchRec = c;} + +/* virtual void SetUseBkg(Bool_t b) {fUseBkg = b;} */ +/* virtual void SetUseSimTPC(Bool_t b) {fUseFastTPC = b;} */ + + virtual void FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection); + + virtual Double_t SetR(Double_t b){fR = b; return fR;} +// TArrayD * GetThrustParamMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom, TArrayD * evsh); + virtual void GetThrustAxis(TVector3 &n01, TVector3 * p_track, Int_t &nTracks); + virtual void GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes); + virtual Bool_t IsPrimChar(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug); + + Double_t Exponent(Double_t x, Double_t * par); + Double_t Exponent2(Double_t x, Double_t * par); + Double_t Gauss(Double_t x, Double_t * par); + Double_t Total(Double_t x, Double_t * par); + + private: + AliAnalysisTaskThreeJets(const AliAnalysisTaskThreeJets&); + AliAnalysisTaskThreeJets& operator = (const AliAnalysisTaskThreeJets&); + + enum {kMaxJets = 6}; + enum {kMaxEvents = 10}; + enum {kJets = 3}; + enum {kTracks = 1000}; + + AliAODEvent *fAOD; // where we take the jets from can be input or output AOD + + TString fBranchRec; // AOD branch name for reconstructed + TString fBranchGen; // AOD brnach for genereated + + Bool_t fUseAODInput; + + Double_t fR; + TList * fList; + + Int_t fGlobVar; + Double_t fXsection; + + TH2F * fhX3X4Rec; + TH2F * fhX3X4Gen; + + TH1F * fhMu35Rec; + TH1F * fhMu34Rec; + TH1F * fhMu45Rec; + + TH1F * fhMu35Gen; + TH1F * fhMu34Gen; + TH1F * fhMu45Gen; + + TH1I * fhInOut; + TH1F * fhThrustRec2; + TH1F * fhThrustRec3; + + TH1F * fhThrustGen2; + TH1F * fhThrustGen3; + + TH1F * fhCGen2; + TH1F * fhCGen3; + + TH1F * fhSGen2; + TH1F * fhSGen3; + + TH1F * fhAGen2; + TH1F * fhAGen3; + + TH1F * fhCRec2; + TH1F * fhCRec3; + + TH1F * fhSRec2; + TH1F * fhSRec3; + + TH1F * fhARec2; + TH1F * fhARec3; + + TH2F * fhX3; + TH2F * fhX4; + TH2F * fhX5; + + TProfile * fhXSec; + TH2F * fhX3X4Rec60; + TH2F * fhX3X4Rec60100; + TH2F * fhX3X4Rec100; + TH2F * fhX3X4Gen60; + TH2F * fhX3X4Gen60100; + TH2F * fhX3X4Gen100; + + TH2F * fhdPhiThrustGen; + TH2F * fhdPhiThrustGenALL; + TH2F * fhdPhiThrustRec; + TH2F * fhdPhiThrustRecALL; + + ClassDef(AliAnalysisTaskThreeJets, 1) +}; + +#endif diff --git a/PWG4/PWG4JetTasksLinkDef.h b/PWG4/PWG4JetTasksLinkDef.h index 272de46faa2..e52a4195d56 100644 --- a/PWG4/PWG4JetTasksLinkDef.h +++ b/PWG4/PWG4JetTasksLinkDef.h @@ -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+; diff --git a/PWG4/libPWG4JetTasks.pkg b/PWG4/libPWG4JetTasks.pkg index e09877d182f..5fb5da32e13 100644 --- a/PWG4/libPWG4JetTasks.pkg +++ b/PWG4/libPWG4JetTasks.pkg @@ -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 index 00000000000..bf0e4f14d34 --- /dev/null +++ b/PWG4/macros/AddTaskJetCorrections.C @@ -0,0 +1,63 @@ +AliAnalysisTaskJetCorrections * AddTaskJetCorrections() +{ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskJetSpectrum", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskJetSpectrum", "This task requires an input event handler"); + return NULL; + } + + AliAnalysisTaskJetCorrections * jetCorr = new AliAnalysisTaskJetCorrections("Jet Corrections"); + + jetCorr->SetBranchGen("jetsMC"); + jetCorr->SetBranchRec("jets"); + jetCorr->SetR(.5); + mgr->AddTask(jetCorr); + + AliAnalysisDataContainer *coutput1_Corr = mgr->CreateContainer("jetCorr", TList::Class(),AliAnalysisManager::kOutputContainer,"jetCorr.root"); + + mgr->ConnectInput (jetCorr, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput (jetCorr, 0, mgr->GetCommonOutputContainer()); + mgr->ConnectOutput (jetCorr, 1, coutput1_Corr ); + + return jetCorr; +} + + +AliAnalysisTaskJetCorrections * AddTaskJetCorrections(AliAnalysisManager* mgr ,AliAnalysisDataContainer *cinput) +{ + if (!mgr) { + ::Error("AddTaskJetSpectrum", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskJetSpectrum", "This task requires an input event handler"); + return NULL; + } + + AliAnalysisTaskJetCorrections * jetCorr = new AliAnalysisTaskJetCorrections("Jet Corrections"); + + jetCorr->SetBranchGen("jetsMC"); + jetCorr->SetBranchRec("jets"); + jetCorr->SetR(.5); + mgr->AddTask(jetCorr); + + AliAnalysisDataContainer * coutpu0 = mgr->CreateContainer("coutpu0", TTree::Class(), + AliAnalysisManager::kExchangeContainer); + AliAnalysisDataContainer *coutput1_jetCorr = mgr->CreateContainer("jetCorr", TList::Class(),AliAnalysisManager::kOutputContainer,"jetCorr.root"); + + mgr->ConnectInput (jetCorr, 0, cinput); + mgr->ConnectOutput (jetCorr, 0, coutpu0); + mgr->ConnectOutput (jetCorr, 1, coutput1_Corr ); + + return jetCorr; +} diff --git a/PWG4/macros/AddTaskThreeJets.C b/PWG4/macros/AddTaskThreeJets.C new file mode 100644 index 00000000000..94019a14fad --- /dev/null +++ b/PWG4/macros/AddTaskThreeJets.C @@ -0,0 +1,73 @@ +AliAnalysisTaskThreeJets * AddTaskThreeJets() +{ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskThreeJets", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskThreeJets", "This task requires an input event handler"); + return NULL; + } + + // Create the task and configure it. + //=========================================================================== + + AliAnalysisTaskThreeJets * threeJets = new AliAnalysisTaskThreeJets("Three Jet Analysis"); + + threeJets->SetBranchGen("jetsMC"); + threeJets->SetBranchRec("jets"); + threeJets->SetR(.5); + mgr->AddTask(threeJets); + + + + + // Create ONLY the output containers for the data produced by the task. + // Get and connect other common input/output containers via the manager as below + //============================================================================== + AliAnalysisDataContainer *coutput1_Corr = mgr->CreateContainer("threeJets", TList::Class(),AliAnalysisManager::kOutputContainer,"threeJets.root"); + + mgr->ConnectInput (threeJets, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput (threeJets, 0, mgr->GetCommonOutputContainer()); + mgr->ConnectOutput (threeJets, 1, coutput1_Corr ); + + return threeJets; +} + + +AliAnalysisTaskThreeJets * AddTaskJetCorrections(AliAnalysisManager* mgr,AliAnalysisDataContainer *cinput) +{ + if (!mgr) { + ::Error("AddTaskJetSpectrum", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskJetSpectrum", "This task requires an input event handler"); + return NULL; + } + + AliAnalysisTaskThreeJets * threeJets = new AliAnalysisTaskThreeJets("ThreeJetAnalysis"); + + threeJets->SetBranchGen("jetsMC"); + threeJets->SetBranchRec("jets"); + threeJets->SetR(.5); + mgr->AddTask(threeJets); + + AliAnalysisDataContainer * coutpu0 = mgr->CreateContainer("coutpu0", TTree::Class(), + AliAnalysisManager::kExchangeContainer); + AliAnalysisDataContainer *coutput1_threeJets = mgr->CreateContainer("threeJets", TList::Class(),AliAnalysisManager::kOutputContainer,"threeJets.root"); + + mgr->ConnectInput (threeJets, 0, cinput); + mgr->ConnectOutput (threeJets, 0, coutpu0); + mgr->ConnectOutput (threeJets, 1, coutput1_Corr ); + + return threeJets; +} + diff --git a/PWG4/macros/AnalysisTrainCAF.C b/PWG4/macros/AnalysisTrainCAF.C index 6c95390098c..534dcfcbc0a 100644 --- a/PWG4/macros/AnalysisTrainCAF.C +++ b/PWG4/macros/AnalysisTrainCAF.C @@ -10,7 +10,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG Bool_t debug = kTRUE; Bool_t useMC = kTRUE; Bool_t readTR = kFALSE; - Bool_t bPROOF = kFALSE; + Bool_t bPROOF = kTRUE; Bool_t bLOCALPAR = kFALSE; // flag that swtiches on loading of local par files insead of loading libs, needed for grid and local testing @@ -24,6 +24,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG Int_t iDIJETAN = 1; Int_t iPWG4SPECTRUM = 1; Int_t iPWG4JFSYSTEMATICS = 1; + Int_t iPWG4JETCORRECTION = 1; + Int_t iPWG4THREEJETS = 1; Int_t iPWG4UE = 0; Int_t iPWG4PID = 0; @@ -40,8 +42,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG TString dataset(ds); TChain *chain = 0; // CKB quick hack for local analysis - gROOT->LoadMacro("CreateESDChain.C"); - TChain *chain = CreateESDChain("tmp.txt",1000); + // gROOT->LoadMacro("CreateESDChain.C"); + // TChain *chain = CreateESDChain("tmp.txt",1000); printf("==================================================================\n"); @@ -57,6 +59,9 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG if (iDIJETAN) printf("= DiJet analysis =\n"); if (iPWG4SPECTRUM)printf("= PWG4 Jet spectrum analysis =\n"); if (iPWG4JFSYSTEMATICS)printf("= PWG4 Jet Finder systematics =\n"); + if (iPWG4JETCORRECTION)printf("= PWG4 Jet Correction =\n"); + if (iPWG4THREEJETS)printf("= PWG4 Three Jets =\n"); + if (iPWG4UE) printf("= PWG4 UE =\n"); printf("==================================================================\n"); if (useMC) printf(":: use MC TRUE\n"); @@ -77,8 +82,8 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG // TProof::Reset("alicecaf"); // One may enable a different ROOT version on CAF - const char* proofNode = "localhost"; - // const char* proofNode = "alicecaf"; + // const char* proofNode = "localhost"; + const char* proofNode = "alicecaf"; // Connect to proof if(bPROOF){ @@ -116,7 +121,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG gProof->EnablePackage("JETAN"); } // --- Enable particle correlation analysis - if (iPWG4UE||iPWG4SPECTRUM||iPWGJFSYSTEMATICS) { + if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS||iPWG4JETCORRECTION||iPWG4THREEJETS) { gProof->UploadPackage("PWG4JetTasks.par"); gProof->EnablePackage("PWG4JetTasks"); } @@ -154,7 +159,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG // --- Enable the JETAN Package if (iJETAN||iJETANESD||iJETANMC||iJETANMC2) gSystem->Load("libJETAN"); // --- Enable particle correlation analysis - if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS)gSystem->Load("libPWG4JetTasks"); + if (iPWG4UE||iPWG4SPECTRUM||iPWG4JFSYSTEMATICS||iPWG4THREEJETS)gSystem->Load("libPWG4JetTasks"); } } @@ -245,6 +250,16 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG AliAnalysisTaskJFSystematics* pwg4jfs = AddTaskJFSystematics(); pwg4jfs->SetDebugLevel(11); } + if (iPWG4JETCORRECTION) { + gROOT->LoadMacro("AddTaskJetCorrections.C"); + AliAnalysisTaskJetCorrections* pwg4jc = AddTaskJetCorrections(); + pwg4jc->SetDebugLevel(11); + } + if (iPWG4THREEJETS) { + gROOT->LoadMacro("AddTaskThreeJets.C"); + AliAnalysisTaskThreeJets* pwg4jjj = AddTaskThreeJets(); + pwg4jjj->SetDebugLevel(11); + } if (iPWG4UE) { gROOT->LoadMacro("AddTaskUE.C"); AliAnalysisTaskUE* ueana = AddTaskUE(); -- 2.43.0