From Salvatore and 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 }
49
50 //________________________________________________________________________
51 AliJetContainer::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),
70   fJetTrigger(0),
71   fParticleContainer(0),
72   fClusterContainer(0),
73   fRho(0),
74   fGeom(0),
75   fRunNumber(0)
76 {
77   // Standard constructor.
78
79 }
80
81 //________________________________________________________________________
82 void AliJetContainer::SetJetArray(AliVEvent *event) 
83 {
84   // Set jet array
85
86   SetArray(event, "AliEmcalJet");
87
88   if(fJetAcceptanceType==kTPC) {
89     AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
90     SetJetEtaPhiTPC();
91   }
92   else if(fJetAcceptanceType==kEMCAL) {
93     AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
94     SetJetEtaPhiEMCAL();
95  }
96 }
97
98
99 //________________________________________________________________________
100 void AliJetContainer::SetEMCALGeometry() {
101   fGeom = AliEMCALGeometry::GetInstance();
102   if (!fGeom) {
103     AliError(Form("%s: Can not create geometry", GetName()));
104     return;
105   }
106 }
107
108 //________________________________________________________________________
109 void AliJetContainer::LoadRho(AliVEvent *event)
110 {
111   // Load rho
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
122 //________________________________________________________________________
123 AliEmcalJet* 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
147 //________________________________________________________________________
148 AliEmcalJet* 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 //________________________________________________________________________
160 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
161   AliEmcalJet *jet = GetJet(i);
162
163   return jet->Pt() - fRho->GetVal()*jet->Area();
164 }
165
166 //________________________________________________________________________
167 AliEmcalJet* 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
177 //________________________________________________________________________
178 AliEmcalJet* 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 //________________________________________________________________________
196 void 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
204 //________________________________________________________________________
205 Bool_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 //________________________________________________________________________
225 Bool_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 //________________________________________________________________________
259 Double_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
269 //________________________________________________________________________
270 void 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 }
306
307 //________________________________________________________________________
308 void 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);
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
321   }
322   else {
323     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
324     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
325     SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
326   }
327 }
328
329 //________________________________________________________________________
330 void AliJetContainer::SetJetEtaPhiTPC()
331 {
332   //Set default cuts for charged jets
333
334   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
335   SetJetPhiLimits(-10, 10);
336 }
337
338 //________________________________________________________________________
339 void 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 }