3 // Container with name, TClonesArray and cuts for jets
7 #include <TClonesArray.h>
9 #include "AliEmcalJet.h"
10 #include "AliVEvent.h"
12 #include "AliEMCALGeometry.h"
13 #include "AliParticleContainer.h"
14 #include "AliClusterContainer.h"
15 #include "AliLocalRhoParameter.h"
17 #include "AliJetContainer.h"
19 ClassImp(AliJetContainer)
21 //________________________________________________________________________
22 AliJetContainer::AliJetContainer():
23 AliEmcalContainer("AliJetContainer"),
24 fJetAcceptanceType(kUser),
44 fLeadingHadronType(0),
48 fParticleContainer(0),
55 // Default constructor.
57 fClassName = "AliEmcalJet";
60 //________________________________________________________________________
61 AliJetContainer::AliJetContainer(const char *name):
62 AliEmcalContainer(name),
63 fJetAcceptanceType(kUser),
83 fLeadingHadronType(0),
87 fParticleContainer(0),
94 // Standard constructor.
96 fClassName = "AliEmcalJet";
99 //________________________________________________________________________
100 void AliJetContainer::SetArray(AliVEvent *event)
104 AliEmcalContainer::SetArray(event);
106 if(fJetAcceptanceType==kTPC) {
107 AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
110 else if(fJetAcceptanceType==kEMCAL) {
111 AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
117 //________________________________________________________________________
118 void AliJetContainer::SetEMCALGeometry() {
119 fGeom = AliEMCALGeometry::GetInstance();
121 AliError(Form("%s: Can not create geometry", GetName()));
126 //________________________________________________________________________
127 void AliJetContainer::LoadRho(AliVEvent *event)
131 if (!fRhoName.IsNull() && !fRho) {
132 fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
134 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
140 //________________________________________________________________________
141 void AliJetContainer::LoadLocalRho(AliVEvent *event)
145 if (!fLocalRhoName.IsNull() && !fLocalRho) {
146 fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
148 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
154 //________________________________________________________________________
155 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
157 // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
162 Int_t tempID = fCurrentID;
164 AliEmcalJet *jetMax = GetNextAcceptJet(0);
165 AliEmcalJet *jet = 0;
167 if (option.Contains("rho")) {
168 while ((jet = GetNextAcceptJet())) {
169 if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) ) jetMax = jet;
173 while ((jet = GetNextAcceptJet())) {
174 if (jet->Pt() > jetMax->Pt()) jetMax = jet;
183 //________________________________________________________________________
184 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
186 //Get i^th jet in array
188 if(i<0 || i>fClArray->GetEntriesFast()) return 0;
189 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
194 //________________________________________________________________________
195 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
197 //Only return jet if is accepted
199 AliEmcalJet *jet = GetJet(i);
200 if(!AcceptJet(jet)) return 0;
205 //________________________________________________________________________
206 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
208 //Get particle with label lab in array
210 Int_t i = GetIndexFromLabel(lab);
214 //________________________________________________________________________
215 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
217 //Get particle with label lab in array
219 Int_t i = GetIndexFromLabel(lab);
220 return GetAcceptJet(i);
223 //________________________________________________________________________
224 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
226 //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
228 if (i>=0) fCurrentID = i;
230 const Int_t njets = GetNEntries();
231 AliEmcalJet *jet = 0;
232 while (fCurrentID < njets && !jet) {
233 jet = GetAcceptJet(fCurrentID);
240 //________________________________________________________________________
241 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
243 //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
245 if (i>=0) fCurrentID = i;
247 const Int_t njets = GetNEntries();
248 AliEmcalJet *jet = 0;
249 while (fCurrentID < njets && !jet) {
250 jet = GetJet(fCurrentID);
257 //________________________________________________________________________
258 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
259 AliEmcalJet *jet = GetJet(i);
261 return jet->Pt() - fRho->GetVal()*jet->Area();
264 //________________________________________________________________________
265 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
266 AliEmcalJet *jet = GetJet(i);
268 return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
271 //________________________________________________________________________
272 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
274 //Get momentum of the i^th jet in array
276 AliEmcalJet *jet = GetJet(i);
277 if(jet) jet->GetMom(mom);
280 //________________________________________________________________________
281 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
283 // Accept jet with a bias.
285 if (fLeadingHadronType == 0) {
286 if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
288 else if (fLeadingHadronType == 1) {
289 if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
292 if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
300 //________________________________________________________________________
301 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
304 // Return true if jet is accepted.
309 if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
312 if (jet->Pt() <= fJetPtCut)
315 if (jet->Area() <= fJetAreaCut)
318 if (jet->AreaEmc() < fAreaEmcCut)
321 if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
324 if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
327 if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
330 if (!AcceptBiasJet(jet))
333 if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
336 if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
339 Double_t jetPhi = jet->Phi();
340 Double_t jetEta = jet->Eta();
342 if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
343 jetPhi -= TMath::Pi() * 2;
345 return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
348 //________________________________________________________________________
349 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
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();
359 //________________________________________________________________________
360 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
362 Double_t maxClusterPt = 0;
363 Double_t maxClusterEta = 0;
364 Double_t maxClusterPhi = 0;
366 Double_t maxTrackPt = 0;
367 Double_t maxTrackEta = 0;
368 Double_t maxTrackPhi = 0;
370 if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
371 AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
373 TLorentzVector nPart;
374 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
376 maxClusterEta = nPart.Eta();
377 maxClusterPhi = nPart.Phi();
378 maxClusterPt = nPart.Pt();
382 if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
383 AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
385 maxTrackEta = track->Eta();
386 maxTrackPhi = track->Phi();
387 maxTrackPt = track->Pt();
391 if (maxTrackPt > maxClusterPt)
392 mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
394 mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
397 //________________________________________________________________________
398 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
401 if (fClusterContainer && fClusterContainer->GetArray()) {
404 AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
406 cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
408 return GetZ(jet,mom);
417 //________________________________________________________________________
418 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
421 if (fParticleContainer && fParticleContainer->GetArray() ) {
424 AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
426 mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
428 return GetZ(jet,mom);
437 //________________________________________________________________________
438 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
441 Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
444 return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
446 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
452 //________________________________________________________________________
453 void AliJetContainer::SetJetEtaPhiEMCAL()
455 //Set default cuts for full jets
457 if(!fGeom) SetEMCALGeometry();
459 SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
461 if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
462 SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
464 SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
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);
474 //________________________________________________________________________
475 void AliJetContainer::SetJetEtaPhiTPC()
477 //Set default cuts for charged jets
479 SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
480 SetJetPhiLimits(-10, 10);
483 //________________________________________________________________________
484 void AliJetContainer::ResetCuts()
486 // Reset cuts to default values
497 fMaxClusterPt = 1000;
499 fLeadingHadronType = 0;
500 fZLeadingEmcCut = 10.;
501 fZLeadingChCut = 10.;
504 //________________________________________________________________________
505 void AliJetContainer::SetClassName(const char *clname)
507 // Set the class name
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));