Fixes for #93216: Request to port to Release: AliESDEvent.h. Replacing UncheckedAt...
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliEsdJetTask.cxx
CommitLineData
48d874ff 1// $Id$
2
3#include "AliEsdJetTask.h"
4#include <TClonesArray.h>
5#include <TParticle.h>
7efbea04 6#include "AliESDJet.h"
48d874ff 7#include "AliAnalysisManager.h"
8#include "AliESDtrack.h"
48d874ff 9#include "AliFJWrapper.h"
10#include "AliESDCaloCluster.h"
11
12ClassImp(AliEsdJetTask)
13
48d874ff 14//________________________________________________________________________
15AliEsdJetTask::AliEsdJetTask(const char *name) :
16 AliAnalysisTaskSE("AliEsdJetTask"),
7efbea04 17 fTracksName("Tracks"),
48d874ff 18 fCaloName("CaloClusters"),
7efbea04 19 fJetsName("Jets"),
48d874ff 20 fAlgo(1),
21 fRadius(0.4),
22 fType(0),
23 fHadCorr(0),
7efbea04 24 fJets(0)
48d874ff 25{
26 // Standard constructor.
7efbea04 27
48d874ff 28 if (!name)
29 return;
7efbea04 30
48d874ff 31 SetName(name);
7efbea04 32 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff 33}
34
35//________________________________________________________________________
36AliEsdJetTask::~AliEsdJetTask()
37{
38 // Destructor
39}
40
41//________________________________________________________________________
42void AliEsdJetTask::UserCreateOutputObjects()
43{
44 // Create user objects.
45
7efbea04 46 fJets = new TClonesArray("AliESDJet");
48d874ff 47 fJets->SetName(fJetsName);
48d874ff 48}
49
50//________________________________________________________________________
51void AliEsdJetTask::UserExec(Option_t *)
52{
53 // Main loop, called for each event.
7efbea04 54 // Add jets to event if not yet there
48d874ff 55
48d874ff 56 if (!(InputEvent()->FindListObject(fJetsName)))
57 InputEvent()->AddObject(fJets);
58
59 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
60 TClonesArray *tracks = 0;
61 TClonesArray *clus = 0;
62 TList *l = InputEvent()->GetList();
63 if ((fType==0)||(fType==1)) {
7efbea04 64 if (fTracksName == "Tracks")
48d874ff 65 am->LoadBranch("Tracks");
7efbea04 66 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
48d874ff 67 if (!tracks) {
7efbea04 68 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
48d874ff 69 return;
70 }
71 }
72 if ((fType==0)||(fType==2)) {
73 if (fCaloName == "CaloClusters")
74 am->LoadBranch("CaloClusters");
75 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
76 if (!clus) {
77 AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
78 return;
79 }
48d874ff 80 }
81
82 FindJets(tracks, clus, fAlgo, fRadius);
83}
84
85//________________________________________________________________________
86void AliEsdJetTask::Terminate(Option_t *)
87{
88 // Called once at the end of the analysis.
89
90}
91
92//________________________________________________________________________
93void AliEsdJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius)
94{
95 // Find jets.
96
97 TString name("kt");
98 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
99 if (algo>=1) {
100 name = "antikt";
101 jalgo = fastjet::antikt_algorithm;
102 }
103
104 AliFJWrapper fjw(name, name);
105 fjw.SetR(radius);
106 fjw.SetAlgorithm(jalgo);
107 fjw.SetMaxRap(0.9);
108 fjw.Clear();
109
110 if (tracks) {
111 const Int_t Ntracks = tracks->GetEntries();
112 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
7efbea04 113 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
48d874ff 114 if (!t)
115 continue;
48d874ff 116 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
117 }
118 }
119 if (clus) {
120 Double_t vertex[3] = {0, 0, 0};
121 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
122 const Int_t Nclus = clus->GetEntries();
123 for (Int_t iClus = 0, iN = 0; iClus < Nclus; ++iClus) {
7efbea04 124 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
48d874ff 125 if (!c->IsEMCAL())
126 continue;
127 TLorentzVector nPart;
128 c->GetMomentum(nPart, vertex);
129 Double_t energy = nPart.P();
130 if (fHadCorr>0) {
131 Int_t imin = static_cast<Int_t>(c->GetEmcCpvDistance());
132 if (imin>=0) {
133 Double_t dPhiMin = c->GetTrackDx();
134 Double_t dEtaMin = c->GetTrackDz();
7efbea04 135 if (dPhiMin<0.05 && dEtaMin<0.025) { // pp cuts!!!
136 AliVTrack *t = dynamic_cast<AliESDtrack*>(tracks->At(imin));
48d874ff 137 if (t) {
7efbea04 138 energy -= fHadCorr*t->P();
139 if (energy<0)
140 continue;
48d874ff 141 }
142 }
143 }
144 }
145 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), energy, -iN-1);
48d874ff 146 }
147 }
148
7efbea04 149 // run jet finder
48d874ff 150 fjw.Run();
48d874ff 151
152 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
153 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
154 if (jets_incl[ij].perp()<1)
155 continue;
900f5135 156 //cout << jets_incl[ij].perp() << " " << jets_incl[ij].eta() << " " << jets_incl[ij].phi() << " " << jets_incl[ij].m() << endl;
157 AliESDJet *jet = new ((*fJets)[jetCount])
158 AliESDJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
159 jet->SetArea(fjw.GetJetArea(ij));
160#if 0 /*tobedone*/
48d874ff 161 jet->GetRefTracks()->Clear();
162 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
163 Double_t neutralE = 0;
164 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
165 Int_t uid = constituents[ic].user_index();
166 if (uid>=0)
167 jet->AddTrack(tracks->At(uid));
168 else {
169 TLorentzVector *nP = static_cast<TLorentzVector*>(fNeutrals->At(-(uid+1)));
170 neutralE += nP->E();
171 jet->AddTrack(nP);
172 }
173 }
174 jet->SetNEF(neutralE/jet->E());
7efbea04 175#endif
48d874ff 176 //jet->Print("");
177 jetCount++;
178 }
179}
48d874ff 180