]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliJetContainer.cxx
From Salvatore and Marta:
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx
CommitLineData
421b12a3 1//
2// Jet Container
3//
4// Author: M. Verweij
5
421b12a3 6#include <TClonesArray.h>
6421eeb0 7
421b12a3 8#include "AliEmcalJet.h"
9#include "AliVEvent.h"
10#include "AliLog.h"
421b12a3 11#include "AliEMCALGeometry.h"
6421eeb0 12#include "AliParticleContainer.h"
13#include "AliClusterContainer.h"
421b12a3 14
15#include "AliJetContainer.h"
16
17ClassImp(AliJetContainer)
18
19//________________________________________________________________________
20AliJetContainer::AliJetContainer():
21 AliEmcalContainer("AliJetContainer"),
22 fJetAcceptanceType(kUser),
23 fJetRadius(0),
24 fRhoName(),
25 fPtBiasJetTrack(0),
26 fPtBiasJetClus(0),
27 fJetPtCut(1),
28 fJetAreaCut(-1),
29 fAreaEmcCut(0),
30 fJetMinEta(-0.9),
31 fJetMaxEta(0.9),
32 fJetMinPhi(-10),
33 fJetMaxPhi(10),
34 fMaxClusterPt(1000),
35 fMaxTrackPt(100),
36 fLeadingHadronType(0),
37 fNLeadingJets(1),
38 fJetBitMap(0),
6421eeb0 39 fJetTrigger(0),
40 fParticleContainer(0),
41 fClusterContainer(0),
421b12a3 42 fRho(0),
56f370b0 43 fGeom(0),
44 fRunNumber(0)
421b12a3 45{
46 // Default constructor.
47
48}
49
50//________________________________________________________________________
51AliJetContainer::AliJetContainer(const char *name):
52 AliEmcalContainer(name),
53 fJetAcceptanceType(kUser),
54 fJetRadius(0),
55 fRhoName(),
56 fPtBiasJetTrack(0),
57 fPtBiasJetClus(0),
58 fJetPtCut(1),
59 fJetAreaCut(-1),
60 fAreaEmcCut(0),
61 fJetMinEta(-0.9),
62 fJetMaxEta(0.9),
63 fJetMinPhi(-10),
64 fJetMaxPhi(10),
65 fMaxClusterPt(1000),
66 fMaxTrackPt(100),
67 fLeadingHadronType(0),
68 fNLeadingJets(1),
69 fJetBitMap(0),
6421eeb0 70 fJetTrigger(0),
71 fParticleContainer(0),
72 fClusterContainer(0),
421b12a3 73 fRho(0),
56f370b0 74 fGeom(0),
75 fRunNumber(0)
421b12a3 76{
77 // Standard constructor.
78
79}
80
421b12a3 81//________________________________________________________________________
82void AliJetContainer::SetJetArray(AliVEvent *event)
83{
84 // Set jet array
85
421b12a3 86 SetArray(event, "AliEmcalJet");
87
56f370b0 88 if(fJetAcceptanceType==kTPC) {
89 AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
421b12a3 90 SetJetEtaPhiTPC();
56f370b0 91 }
92 else if(fJetAcceptanceType==kEMCAL) {
93 AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
421b12a3 94 SetJetEtaPhiEMCAL();
56f370b0 95 }
421b12a3 96}
97
56f370b0 98
421b12a3 99//________________________________________________________________________
100void AliJetContainer::SetEMCALGeometry() {
101 fGeom = AliEMCALGeometry::GetInstance();
102 if (!fGeom) {
103 AliError(Form("%s: Can not create geometry", GetName()));
104 return;
105 }
106}
107
108//________________________________________________________________________
109void AliJetContainer::LoadRho(AliVEvent *event)
110{
6421eeb0 111 // Load rho
421b12a3 112
113 if (!fRhoName.IsNull() && !fRho) {
114 fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
115 if (!fRho) {
116 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
117 return;
118 }
119 }
120}
121
6421eeb0 122//________________________________________________________________________
123AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt) const
124{
125 // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
126
127 TString option(opt);
128 option.ToLower();
129
130 AliEmcalJet *jetMax = GetNextAcceptJet(0);
131 AliEmcalJet *jet = 0;
132
133 if (option.Contains("rho")) {
134 while ((jet = GetNextAcceptJet())) {
135 if (jet->Pt()-jet->Area()*fRho->GetVal() > jetMax->Pt()-jetMax->Area()*fRho->GetVal()) jetMax = jet;
136 }
137 }
138 else {
139 while ((jet = GetNextAcceptJet())) {
140 if (jet->Pt() > jetMax->Pt()) jetMax = jet;
141 }
142 }
143
144 return jetMax;
145}
146
421b12a3 147//________________________________________________________________________
148AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
149
150 //Get i^th jet in array
151
152 if(i<0 || i>fClArray->GetEntriesFast()) return 0;
153 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
154 return jet;
155
156}
157
158
159//________________________________________________________________________
160Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
161 AliEmcalJet *jet = GetJet(i);
162
163 return jet->Pt() - fRho->GetVal()*jet->Area();
164}
165
166//________________________________________________________________________
167AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
168
169 //Only return jet if is accepted
170
171 AliEmcalJet *jet = GetJet(i);
172 if(!AcceptJet(jet)) return 0;
173
174 return jet;
175}
176
6421eeb0 177//________________________________________________________________________
178AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) const {
179
180 //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
181
182 static Int_t counter = -1;
183 if (i>=0) counter = i;
184
185 const Int_t njets = GetNEntries();
186 AliEmcalJet *jet = 0;
187 while (counter < njets && !jet) {
188 jet = GetAcceptJet(counter);
189 counter++;
190 }
191
192 return jet;
193}
194
195//________________________________________________________________________
196void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
197{
198 //Get momentum of the i^th jet in array
199
200 AliEmcalJet *jet = GetJet(i);
201 if(jet) jet->GetMom(mom);
202}
203
421b12a3 204//________________________________________________________________________
205Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
206{
207 // Accept jet with a bias.
208
209 if (fLeadingHadronType == 0) {
210 if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
211 }
212 else if (fLeadingHadronType == 1) {
213 if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
214 }
215 else {
216 if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
217 }
218
219 return kTRUE;
220
221
222}
223
224//________________________________________________________________________
225Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
226{
227 // Return true if jet is accepted.
228 if (!jet)
229 return kFALSE;
230
231 if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
232 return kFALSE;
233
234 if (jet->Pt() <= fJetPtCut)
235 return kFALSE;
236
237 if (jet->Area() <= fJetAreaCut)
238 return kFALSE;
239
240 if (jet->AreaEmc()<fAreaEmcCut)
241 return kFALSE;
242
243 if (!AcceptBiasJet(jet))
244 return kFALSE;
245
246 if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
247 return kFALSE;
248
249 Double_t jetPhi = jet->Phi();
250 Double_t jetEta = jet->Eta();
251
252 if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
253 jetPhi -= TMath::Pi() * 2;
254
255 return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
256}
257
258//________________________________________________________________________
259Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
260{
261 if (fLeadingHadronType == 0) // charged leading hadron
262 return jet->MaxTrackPt();
263 else if (fLeadingHadronType == 1) // neutral leading hadron
264 return jet->MaxClusterPt();
265 else // charged or neutral
266 return jet->MaxPartPt();
267}
268
6421eeb0 269//________________________________________________________________________
270void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
271{
272 Double_t maxClusterPt = 0;
273 Double_t maxClusterEta = 0;
274 Double_t maxClusterPhi = 0;
275
276 Double_t maxTrackPt = 0;
277 Double_t maxTrackEta = 0;
278 Double_t maxTrackPhi = 0;
279
280 if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
281 AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
282 if (cluster) {
283 TLorentzVector nPart;
284 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
285
286 maxClusterEta = nPart.Eta();
287 maxClusterPhi = nPart.Phi();
288 maxClusterPt = nPart.Pt();
289 }
290 }
291
292 if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
293 AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
294 if (track) {
295 maxTrackEta = track->Eta();
296 maxTrackPhi = track->Phi();
297 maxTrackPt = track->Pt();
298 }
299 }
300
301 if (maxTrackPt > maxClusterPt)
302 mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
303 else
304 mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
305}
421b12a3 306
307//________________________________________________________________________
308void AliJetContainer::SetJetEtaPhiEMCAL()
309{
310 //Set default cuts for full jets
311
312 if(!fGeom) SetEMCALGeometry();
313 if(fGeom) {
314 SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
56f370b0 315
316 if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
317 SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
318 else
319 SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
320
421b12a3 321 }
322 else {
56f370b0 323 AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
421b12a3 324 SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
325 SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
326 }
421b12a3 327}
328
329//________________________________________________________________________
330void AliJetContainer::SetJetEtaPhiTPC()
331{
56f370b0 332 //Set default cuts for charged jets
421b12a3 333
56f370b0 334 SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
421b12a3 335 SetJetPhiLimits(-10, 10);
421b12a3 336}
337
338//________________________________________________________________________
339void AliJetContainer::ResetCuts()
340{
341 // Reset cuts to default values
342
343 fPtBiasJetTrack = 0;
344 fPtBiasJetClus = 0;
345 fJetPtCut = 1;
346 fJetAreaCut = -1;
347 fAreaEmcCut = 0;
348 fJetMinEta = -0.9;
349 fJetMaxEta = 0.9;
350 fJetMinPhi = -10;
351 fJetMaxPhi = 10;
352 fMaxClusterPt = 1000;
353 fMaxTrackPt = 100;
354 fLeadingHadronType = 0;
355
356}