JetTasks/AliAnalysisTaskUE.cxx
JetTasks/AliAnalysisTaskJetSpectrum.cxx
JetTasks/AliAnalysisTaskJFSystematics.cxx
+ JetTasks/AliAnalysisTaskJetCorrections.cxx;
+ JetTasks/AliAnalysisTaskThreeJets.cxx;
JetTasks/AliAnalysisHelperJetTasks.cxx
JetTasks/AliAnaESDSpectraQA.cxx
JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx
--- /dev/null
+#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();
+ }
+}
+//______________________________________________________________________________________________________
+
+
--- /dev/null
+#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
--- /dev/null
+#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]);
+}
+
+
--- /dev/null
+#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
#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+;
#-*- 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)
--- /dev/null
+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
--- /dev/null
+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
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
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
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
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
// 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
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
// --- 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
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
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff