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