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