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