prevent floating point
[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()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
170         jetMax = jet;
171     }
172   }
173   else {
174     while ((jet = GetNextAcceptJet())) {
175       if (jet->Pt() > jetMax->Pt()) jetMax = jet;
176     }
177   }
178
179   fCurrentID = tempID;
180
181   return jetMax;
182 }
183
184 //________________________________________________________________________
185 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
186
187   //Get i^th jet in array
188
189   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
190   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
191   return jet;
192
193 }
194
195 //________________________________________________________________________
196 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
197
198   //Only return jet if is accepted
199
200   AliEmcalJet *jet = GetJet(i);
201   if(!AcceptJet(jet)) return 0;
202
203   return jet;
204 }
205
206 //________________________________________________________________________
207 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
208
209   //Get particle with label lab in array
210   
211   Int_t i = GetIndexFromLabel(lab);
212   return GetJet(i);
213 }
214
215 //________________________________________________________________________
216 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
217
218   //Get particle with label lab in array
219   
220   Int_t i = GetIndexFromLabel(lab);
221   return GetAcceptJet(i);
222 }
223
224 //________________________________________________________________________
225 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
226
227   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
228
229   if (i>=0) fCurrentID = i;
230
231   const Int_t njets = GetNEntries();
232   AliEmcalJet *jet = 0;
233   while (fCurrentID < njets && !jet) { 
234     jet = GetAcceptJet(fCurrentID);
235     fCurrentID++;
236   }
237
238   return jet;
239 }
240
241 //________________________________________________________________________
242 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
243
244   //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
245
246   if (i>=0) fCurrentID = i;
247
248   const Int_t njets = GetNEntries();
249   AliEmcalJet *jet = 0;
250   while (fCurrentID < njets && !jet) { 
251     jet = GetJet(fCurrentID);
252     fCurrentID++;
253   }
254
255   return jet;
256 }
257
258 //________________________________________________________________________
259 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
260   AliEmcalJet *jet = GetJet(i);
261
262   return jet->Pt() - fRho->GetVal()*jet->Area();
263 }
264
265 //________________________________________________________________________
266 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
267   AliEmcalJet *jet = GetJet(i);
268
269   return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
270 }
271
272 //________________________________________________________________________
273 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
274 {
275   //Get momentum of the i^th jet in array
276
277   AliEmcalJet *jet = GetJet(i);
278   if(jet) jet->GetMom(mom);
279 }
280
281 //________________________________________________________________________
282 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
283
284   // Accept jet with a bias.
285
286   if (fLeadingHadronType == 0) {
287     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
288   }
289   else if (fLeadingHadronType == 1) {
290     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
291   }
292   else {
293     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
294   }
295
296   return kTRUE;
297
298
299 }
300
301 //________________________________________________________________________
302 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
303 {   
304
305    // Return true if jet is accepted.
306
307    if (!jet)
308       return kFALSE;
309
310    if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
311       return kFALSE;
312
313    if (jet->Pt() <= fJetPtCut) 
314       return kFALSE;
315
316    if (jet->Area() <= fJetAreaCut) 
317       return kFALSE;
318
319    if (jet->AreaEmc() < fAreaEmcCut)
320       return kFALSE;
321    
322    if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
323       return kFALSE;
324    
325    if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
326       return kFALSE;
327
328    if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
329       return kFALSE;
330    
331    if (!AcceptBiasJet(jet))
332       return kFALSE;
333    
334    if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
335       return kFALSE;
336    
337    if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
338       return kFALSE;
339    
340    Double_t jetPhi = jet->Phi();
341    Double_t jetEta = jet->Eta();
342    
343    if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
344       jetPhi -= TMath::Pi() * 2;
345    
346    return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
347 }
348
349 //________________________________________________________________________
350 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
351 {
352   if (fLeadingHadronType == 0)       // charged leading hadron
353     return jet->MaxTrackPt();
354   else if (fLeadingHadronType == 1)  // neutral leading hadron
355     return jet->MaxClusterPt();
356   else                               // charged or neutral
357     return jet->MaxPartPt();
358 }
359
360 //________________________________________________________________________
361 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
362 {
363   Double_t maxClusterPt = 0;
364   Double_t maxClusterEta = 0;
365   Double_t maxClusterPhi = 0;
366   
367   Double_t maxTrackPt = 0;
368   Double_t maxTrackEta = 0;
369   Double_t maxTrackPhi = 0;
370       
371   if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
372     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
373     if (cluster) {
374       TLorentzVector nPart;
375       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
376       
377       maxClusterEta = nPart.Eta();
378       maxClusterPhi = nPart.Phi();
379       maxClusterPt = nPart.Pt();
380     }
381   }
382       
383   if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
384     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
385     if (track) {
386       maxTrackEta = track->Eta();
387       maxTrackPhi = track->Phi();
388       maxTrackPt = track->Pt();
389     }
390   }
391       
392   if (maxTrackPt > maxClusterPt) 
393     mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
394   else 
395     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
396 }
397
398 //________________________________________________________________________
399 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
400 {
401
402   if (fClusterContainer && fClusterContainer->GetArray()) {
403     TLorentzVector mom;
404     
405     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
406     if (cluster) {
407       cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
408       
409       return GetZ(jet,mom);
410     }
411     else
412       return -1;
413   }
414   else
415     return -1;
416 }
417
418 //________________________________________________________________________
419 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
420 {
421
422   if (fParticleContainer && fParticleContainer->GetArray() ) {
423     TLorentzVector mom;
424     
425     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
426     if (track) {
427       mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
428       
429       return GetZ(jet,mom);
430     }
431     else
432       return -1;
433   }
434   else
435     return -1;
436 }
437
438 //________________________________________________________________________
439 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
440 {
441
442   Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
443
444   if(pJetSq>1e-6)
445     return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
446   else {
447     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
448     return -1;
449   }
450
451 }
452
453 //________________________________________________________________________
454 void AliJetContainer::SetJetEtaPhiEMCAL()
455 {
456   //Set default cuts for full jets
457
458   if(!fGeom) SetEMCALGeometry();
459   if(fGeom) {
460     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
461
462     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
463       SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
464     else
465       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
466
467   }
468   else {
469     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
470     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
471     SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
472   }
473 }
474
475 //________________________________________________________________________
476 void AliJetContainer::SetJetEtaPhiTPC()
477 {
478   //Set default cuts for charged jets
479
480   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
481   SetJetPhiLimits(-10, 10);
482 }
483
484 //________________________________________________________________________
485 void AliJetContainer::ResetCuts() 
486 {
487   // Reset cuts to default values
488
489   fPtBiasJetTrack = 0;
490   fPtBiasJetClus = 0;
491   fJetPtCut = 1;
492   fJetAreaCut = -1;
493   fAreaEmcCut = 0;
494   fJetMinEta = -0.9;
495   fJetMaxEta = 0.9;
496   fJetMinPhi = -10;
497   fJetMaxPhi = 10;
498   fMaxClusterPt = 1000;
499   fMaxTrackPt = 100;
500   fLeadingHadronType = 0;
501   fZLeadingEmcCut = 10.;
502   fZLeadingChCut = 10.;
503 }
504
505 //________________________________________________________________________
506 void AliJetContainer::SetClassName(const char *clname)
507 {
508   // Set the class name
509
510   TClass cls(clname);
511   if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
512   else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
513 }