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),
45 fLeadingHadronType(0),
50 fParticleContainer(0),
58 // Default constructor.
60 fClassName = "AliEmcalJet";
63 //________________________________________________________________________
64 AliJetContainer::AliJetContainer(const char *name):
65 AliEmcalContainer(name),
66 fJetAcceptanceType(kUser),
87 fLeadingHadronType(0),
92 fParticleContainer(0),
100 // Standard constructor.
102 fClassName = "AliEmcalJet";
105 //________________________________________________________________________
106 void AliJetContainer::SetArray(AliVEvent *event)
110 AliEmcalContainer::SetArray(event);
112 if(fJetAcceptanceType==kTPC) {
113 AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
116 else if(fJetAcceptanceType==kEMCAL) {
117 AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
123 //________________________________________________________________________
124 void AliJetContainer::SetEMCALGeometry() {
125 fGeom = AliEMCALGeometry::GetInstance();
127 AliError(Form("%s: Can not create geometry", GetName()));
132 //________________________________________________________________________
133 void AliJetContainer::LoadRho(AliVEvent *event)
137 if (!fRhoName.IsNull() && !fRho) {
138 fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
140 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
146 //________________________________________________________________________
147 void AliJetContainer::LoadLocalRho(AliVEvent *event)
151 if (!fLocalRhoName.IsNull() && !fLocalRho) {
152 fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
154 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
160 //________________________________________________________________________
161 void AliJetContainer::LoadRhoMass(AliVEvent *event)
165 if (!fRhoMassName.IsNull() && !fRhoMass) {
166 fRhoMass = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoMassName));
168 AliError(Form("%s: Could not retrieve rho_mass %s!", GetName(), fRhoMassName.Data()));
175 //________________________________________________________________________
176 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
178 // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
183 Int_t tempID = fCurrentID;
185 AliEmcalJet *jetMax = GetNextAcceptJet(0);
186 AliEmcalJet *jet = 0;
188 if (option.Contains("rho")) {
189 while ((jet = GetNextAcceptJet())) {
190 if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
195 while ((jet = GetNextAcceptJet())) {
196 if (jet->Pt() > jetMax->Pt()) jetMax = jet;
205 //________________________________________________________________________
206 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
208 //Get i^th jet in array
210 if(i<0 || i>fClArray->GetEntriesFast()) return 0;
211 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
216 //________________________________________________________________________
217 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
219 //Only return jet if is accepted
221 AliEmcalJet *jet = GetJet(i);
222 if(!AcceptJet(jet)) return 0;
227 //________________________________________________________________________
228 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
230 //Get particle with label lab in array
232 Int_t i = GetIndexFromLabel(lab);
236 //________________________________________________________________________
237 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
239 //Get particle with label lab in array
241 Int_t i = GetIndexFromLabel(lab);
242 return GetAcceptJet(i);
245 //________________________________________________________________________
246 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
248 //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
250 if (i>=0) fCurrentID = i;
252 const Int_t njets = GetNEntries();
253 AliEmcalJet *jet = 0;
254 while (fCurrentID < njets && !jet) {
255 jet = GetAcceptJet(fCurrentID);
262 //________________________________________________________________________
263 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
265 //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
267 if (i>=0) fCurrentID = i;
269 const Int_t njets = GetNEntries();
270 AliEmcalJet *jet = 0;
271 while (fCurrentID < njets && !jet) {
272 jet = GetJet(fCurrentID);
279 //________________________________________________________________________
280 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
281 AliEmcalJet *jet = GetJet(i);
283 return jet->Pt() - fRho->GetVal()*jet->Area();
286 //________________________________________________________________________
287 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
288 AliEmcalJet *jet = GetJet(i);
290 return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
293 //________________________________________________________________________
294 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
296 //Get momentum of the i^th jet in array
298 AliEmcalJet *jet = GetJet(i);
299 if(jet) jet->GetMom(mom);
302 //________________________________________________________________________
303 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
305 // Accept jet with a bias.
307 if (fLeadingHadronType == 0) {
308 if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
310 else if (fLeadingHadronType == 1) {
311 if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
314 if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
322 //________________________________________________________________________
323 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
326 // Return true if jet is accepted.
329 AliDebug(11,"No jet found");
333 if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap) {
334 AliDebug(11,"Cut rejecting jet: Bit map");
338 if (jet->Pt() <= fJetPtCut) {
339 AliDebug(11,"Cut rejecting jet: JetPtCut");
343 if (jet->Area() <= fJetAreaCut) {
344 AliDebug(11,"Cut rejecting jet: Area");
348 if (jet->AreaEmc() < fAreaEmcCut) {
349 AliDebug(11,"Cut rejecting jet: AreaEmc");
353 if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut) {
354 AliDebug(11,"Cut rejecting jet: ZLeading");
358 if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut) {
359 AliDebug(11,"Cut rejecting jet: ZLeadEmc");
363 if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut) {
364 AliDebug(11,"Cut rejecting jet: NEF");
368 if (!AcceptBiasJet(jet)) {
369 AliDebug(11,"Cut rejecting jet: Bias");
373 if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt) {
374 AliDebug(11,"Cut rejecting jet: MaxTrClPt");
378 if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection)) {
379 AliDebug(11,"Cut rejecting jet: Flavour");
383 if(fTagStatus>-1 && jet->GetTagStatus()!=fTagStatus) {
384 AliDebug(11,"Cut rejecting jet: tag status");
388 Double_t jetPhi = jet->Phi();
389 Double_t jetEta = jet->Eta();
391 if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
392 jetPhi -= TMath::Pi() * 2;
394 return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
397 //________________________________________________________________________
398 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
400 if (fLeadingHadronType == 0) // charged leading hadron
401 return jet->MaxTrackPt();
402 else if (fLeadingHadronType == 1) // neutral leading hadron
403 return jet->MaxClusterPt();
404 else // charged or neutral
405 return jet->MaxPartPt();
408 //________________________________________________________________________
409 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
411 Double_t maxClusterPt = 0;
412 Double_t maxClusterEta = 0;
413 Double_t maxClusterPhi = 0;
415 Double_t maxTrackPt = 0;
416 Double_t maxTrackEta = 0;
417 Double_t maxTrackPhi = 0;
419 if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
420 AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
422 TLorentzVector nPart;
423 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
425 maxClusterEta = nPart.Eta();
426 maxClusterPhi = nPart.Phi();
427 maxClusterPt = nPart.Pt();
431 if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
432 AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
434 maxTrackEta = track->Eta();
435 maxTrackPhi = track->Phi();
436 maxTrackPt = track->Pt();
440 if (maxTrackPt > maxClusterPt)
441 mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
443 mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
446 //________________________________________________________________________
447 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
450 if (fClusterContainer && fClusterContainer->GetArray()) {
453 AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
455 cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
457 return GetZ(jet,mom);
466 //________________________________________________________________________
467 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
470 if (fParticleContainer && fParticleContainer->GetArray() ) {
473 AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
475 mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
477 return GetZ(jet,mom);
486 //________________________________________________________________________
487 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
490 Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
493 return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
495 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
501 //________________________________________________________________________
502 void AliJetContainer::SetJetEtaPhiEMCAL()
504 //Set default cuts for full jets
506 if(!fGeom) SetEMCALGeometry();
508 SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
510 if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
511 SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
513 SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
517 AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
518 SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
519 SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
523 //________________________________________________________________________
524 void AliJetContainer::SetJetEtaPhiTPC()
526 //Set default cuts for charged jets
528 SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
529 SetJetPhiLimits(-10, 10);
532 //________________________________________________________________________
533 void AliJetContainer::PrintCuts()
535 // Reset cuts to default values
537 Printf("PtBiasJetTrack: %f",fPtBiasJetTrack);
538 Printf("PtBiasJetClus: %f",fPtBiasJetClus);
539 Printf("JetPtCut: %f", fJetPtCut);
540 Printf("JetAreaCut: %f",fJetAreaCut);
541 Printf("AreaEmcCut: %f",fAreaEmcCut);
542 Printf("JetMinEta: %f", fJetMinEta);
543 Printf("JetMaxEta: %f", fJetMaxEta);
544 Printf("JetMinPhi: %f", fJetMinPhi);
545 Printf("JetMaxPhi: %f", fJetMaxPhi);
546 Printf("MaxClusterPt: %f",fMaxClusterPt);
547 Printf("MaxTrackPt: %f",fMaxTrackPt);
548 Printf("LeadingHadronType: %d",fLeadingHadronType);
549 Printf("ZLeadingEmcCut: %f",fZLeadingEmcCut);
550 Printf("ZLeadingChCut: %f",fZLeadingChCut);
554 //________________________________________________________________________
555 void AliJetContainer::ResetCuts()
557 // Reset cuts to default values
568 fMaxClusterPt = 1000;
570 fLeadingHadronType = 0;
571 fZLeadingEmcCut = 10.;
572 fZLeadingChCut = 10.;
575 //________________________________________________________________________
576 void AliJetContainer::SetClassName(const char *clname)
578 // Set the class name
581 if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
582 else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
585 //________________________________________________________________________
586 Double_t AliJetContainer::GetFractionSharedPt(AliEmcalJet *jet1) const
589 // Get fraction of shared pT between matched full and charged jet
590 // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
591 // Only works if tracks array of both jets is the same
594 AliEmcalJet *jet2 = jet1->ClosestJet();
597 Double_t fraction = 0.;
598 Double_t jetPt2 = jet2->Pt();
602 AliVParticle *vpf = 0x0;
604 for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
605 Int_t idx = (Int_t)jet2->TrackAt(icc);
607 for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
608 if(idx == jet1->TrackAt(icf) && iFound==0 ) {
610 vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
611 if(vpf) sumPt += vpf->Pt();
616 fraction = sumPt/jetPt2;