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