7 #include <TProfile2D.h>
12 #include <TClonesArray.h>
16 #include <TLorentzVector.h>
17 #include <TParameter.h>
19 #include "AliAODEvent.h"
20 #include "AliAODInputHandler.h"
21 #include "AliAnalysisManager.h"
22 #include "AliAnalysisTask.h"
23 #include "AliCentrality.h"
25 #include "AliESDEvent.h"
26 #include "AliESDInputHandler.h"
27 #include "AliVParticle.h"
28 #include "AliVTrack.h"
29 #include "AliInputEventHandler.h"
30 #include "AliMCEvent.h"
32 #include "AliGenEventHeader.h"
33 #include "AliGenPythiaEventHeader.h"
54 fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
55 fMCParticleArrName(""), fMCParticleArray(0x0),
56 fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1),
57 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*
pi),
59 fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
60 fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
61 fRhoName(""), fRho(0x0), fRhoValue(0),
62 fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
63 fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
64 fRunPL(kFALSE), fRunDL(kTRUE),
65 fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
75 for(
Int_t i=0; i<kNTrig; i++)
78 fhCentrality[i] = 0x0;
81 fhDLJetPtVsCent[i] = 0x0;
82 fhPLJetPtVsCent[i] = 0x0;
92 fhJetPtGeoMatch[i] = 0x0;
93 fhJetPtEnMatch[i] = 0x0;
94 fhJetPhiGeoMatch[i] = 0x0;
95 fhJetPhiEnMatch[i] = 0x0;
97 for(
Int_t i=0; i<kNTT; i++)
100 { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
102 { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
110 fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
111 fMCParticleArrName(
""), fMCParticleArray(0x0),
112 fTrackArrName(
""), fTrackArray(0x0), fTriggerTrkIndex(-1),
113 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*
pi),
115 fRadius(0.4), fJetArrName(
""), fPLJetArrName(
""), fDLJetArrName(
""),
116 fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
117 fRhoName(
""), fRho(0x0), fRhoValue(0),
118 fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
119 fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
120 fRunPL(kFALSE), fRunDL(kTRUE),
121 fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
125 DefineInput(0, TChain::Class());
126 DefineOutput(1, TList::Class());
171 const Int_t nJetPtBins = 40;
172 const Float_t lowJetPtBin=-50, upJetPtBin=150;
174 const Int_t nTrkPtBins = 100;
175 const Float_t lowTrkPtBin=0, upTrkPtBin=100;
179 const Int_t dimJetqa = 3;
180 const Int_t nBinsJetqa[dimJetqa] = {nJetPtBins, 30, 11};
181 const Double_t lowBinJetqa[dimJetqa] = {lowJetPtBin, 0, 0};
182 const Double_t hiBinJetqa[dimJetqa] = {upJetPtBin, 30, 11};
185 const Int_t dimTT = 3;
186 const Int_t nBinsTT[dimTT] = {nTrkPtBins, 30, 11};
187 const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0, 0};
188 const Double_t hiBinTT[dimTT] = {upTrkPtBin, 30, 11};
190 const Int_t dimHJet = 6;
191 const Int_t nBinsHJet[dimHJet] = {nTrkPtBins, nJetPtBins, 140, 8, 30, 11};
192 const Double_t lowBinHJet[dimHJet] = {lowTrkPtBin, lowJetPtBin,
pi-4.95, 0, 0, 0};
193 const Double_t hiBinHJet[dimHJet] = {upTrkPtBin, upJetPtBin,
pi+2.05, 0.8, 30, 11};
196 const Int_t dimMthPt = 6;
197 const Int_t nBinsMthPt[dimMthPt] = {nJetPtBins, nJetPtBins, 20, 20, 30, 11};
198 const Double_t lowBinMthPt[dimMthPt] = {lowJetPtBin, lowJetPtBin, -0.95, 0, 0, 0};
199 const Double_t hiBinMthPt[dimMthPt] = {upJetPtBin, upJetPtBin, 1.05, 1, 30, 11};
201 const Int_t dimMthPhi = 7;
202 const Int_t nBinsMthPhi[dimMthPhi] = {nTrkPtBins, nJetPtBins/2, 70, 70, 10, 1, 11};
203 const Double_t lowBinMthPhi[dimMthPhi] = {lowTrkPtBin, lowJetPtBin,
pi-4.95,
pi-4.95, 0, 0, 0};
204 const Double_t hiBinMthPhi[dimMthPhi] = {upTrkPtBin, upJetPtBin,
pi+2.05,
pi+2.05, 1, 10, 11};
210 fhEventStat =
new TH1F(
"fhEventStat",
"Event statistics for jet analysis",9,0,9);
214 fhEventStat->GetXaxis()->SetBinLabel(4,
"Vtx+10cm");
216 fhEventStat->GetXaxis()->SetBinLabel(6,
"kCentral");
217 fhEventStat->GetXaxis()->SetBinLabel(7,
"kSemiCentral");
222 fhPtHardBins =
new TH1F(
"fhPtHardBins",
"Number of events in each pT hard bin",11,0,11);
225 const char *triggerName[
kNTrig] = {
"kMB",
"kEGA",
"kEJE"};
231 fhVtxZ[i] =
new TH1F(Form(
"%s_fhVtxZ",triggerName[i]),Form(
"%s: z distribution of event vertexz;z(cm)",triggerName[i]),400,-20,20);
234 fhCentrality[i] =
new TH1F(Form(
"%s_fhCentrality",triggerName[i]),Form(
"%s: Event centrality;centrality",triggerName[i]),100,0,100);
237 fhRhoVsCent[i] =
new TH2F(Form(
"%s_fhRhoVsCent",triggerName[i]),Form(
"%s: Rho vs centrality (R=%1.1f);centrality;Rho",triggerName[i],
fRadius),100,0,100,300,0,300);
242 fhPLJetPtVsCent[i] =
new THnSparseF(Form(
"%s_fhPLJetPtVsCent",triggerName[i]),Form(
"PYTHIA: jet p_{T} vs centrality vs pT hard bin (particle-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",
fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
245 fhDLJetPtVsCent[i] =
new THnSparseF(Form(
"%s_fhDLJetPtVsCent",triggerName[i]),Form(
"PYTHIA: jet p_{T} vs centrality vs pT hard bin (detector-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",
fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
251 fhPLTT[i] =
new THnSparseF(Form(
"%s_fhPLTT",triggerName[i]),Form(
"PYTHIA: TT p_{T} vs centrality vs pT hard bin (particle-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
254 fhDLTT[i] =
new THnSparseF(Form(
"%s_fhDLTT",triggerName[i]),Form(
"PYTHIA: TT p_{T} vs centrality vs pT hard bin (detector-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
257 fhTTPt[i] =
new THnSparseF(Form(
"%s_fhTTPt",triggerName[i]),Form(
"Embedded: TT p_{T} vs centrality vs pT hard bin;p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
260 fhPLHJet[i] =
new THnSparseF(Form(
"%s_fhPLHJet",triggerName[i]),Form(
"PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (particle-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",
fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
263 fhDLHJet[i] =
new THnSparseF(Form(
"%s_fhDLHJet",triggerName[i]),Form(
"PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (detector-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",
fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
266 fhHJet[i] =
new THnSparseF(Form(
"%s_fhHJet",triggerName[i]),Form(
"Embedded: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",
fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
272 fhJetPtGeoMatch[i] =
new THnSparseF(Form(
"%s_fhJetPtGeoMatch",triggerName[i]),Form(
"Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",
fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
275 fhJetPtEnMatch[i] =
new THnSparseF(Form(
"%s_fhJetPtEnMatch",triggerName[i]),Form(
"Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",
fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
278 fhJetPhiGeoMatch[i] =
new THnSparseF(Form(
"%s_fhJetPhiGeoMatch",triggerName[i]),Form(
"Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",
fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
281 fhJetPhiEnMatch[i] =
new THnSparseF(Form(
"%s_fhJetPhiEnMatch",triggerName[i]),Form(
"Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",
fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
288 for(
Int_t i=0; i<nObj; i++)
291 if (obj->IsA()->InheritsFrom(
"THnSparse" ))
293 THnSparseF *hn = (THnSparseF*)obj;
309 AliDebug(5,
"Entering UserExec");
314 AliError(
"Input event not available");
317 AliDebug(5,
"Got the input event");
324 AliError(
"The object for pt hard bin information is not available!");
329 fileName.Remove(0,50);
330 fileName.Remove(fileName.Index(
"/"));
333 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
335 if(
fPeriod.Contains(
"lhc11h",TString::kIgnoreCase))
338 else if (trigger & AliVEvent::kCentral) {
fTriggerType=0; }
339 else if (trigger & AliVEvent::kSemiCentral) {
fTriggerType=0; }
340 else if (trigger & AliVEvent::kEMCEGA) {
fTriggerType=1; }
341 else if (trigger & AliVEvent::kEMCEJE) {
fTriggerType=2; }
343 else if(
fPeriod.Contains(
"lhc10h",TString::kIgnoreCase))
347 else if(
fPeriod.Contains(
"lhc12a15a",TString::kIgnoreCase))
356 const AliVVertex* vtx =
fEvent->GetPrimaryVertex();
357 if (!vtx || vtx->GetNContributors()<1)
return;
360 if (TMath::Abs(vtx->GetZ())>
fMaxVtxZ)
return;
365 if(trigger & AliVEvent::kCentral)
fhEventStat->Fill(5.5);
366 else if (trigger & AliVEvent::kCentral)
fhEventStat->Fill(6.5);
377 fCentrality = centrality->GetCentralityPercentile(
"V0M");
392 AliError(Form(
"Could not retrieve tracks %s!",
fTrackArrName.Data()));
395 if (!
fTrackArray->GetClass()->GetBaseClass(
"AliVParticle"))
397 AliError(Form(
"%s: Collection %s does not contain AliVParticle objects!", GetName(),
fTrackArrName.Data()));
422 AliError(Form(
"Could not retrieve rho %s!",
fRhoName.Data()));
439 AliError(Form(
"%s: Could not retrieve jets %s!", GetName(),
fJetArrName.Data()));
442 if (!
fJetArray->GetClass()->GetBaseClass(
"AliEmcalJet"))
444 AliError(Form(
"%s: Collection %s does not contain AliEmcalJet objects!", GetName(),
fJetArrName.Data()));
455 AliError(Form(
"%s: Could not retrieve jets %s!", GetName(),
fPLJetArrName.Data()));
458 if (!
fPLJetArray->GetClass()->GetBaseClass(
"AliEmcalJet"))
460 AliError(Form(
"%s: Collection %s does not contain AliEmcalJet objects!", GetName(),
fPLJetArrName.Data()));
471 AliError(Form(
"%s: Could not retrieve jets %s!", GetName(),
fDLJetArrName.Data()));
474 if (!
fDLJetArray->GetClass()->GetBaseClass(
"AliEmcalJet"))
476 AliError(Form(
"%s: Collection %s does not contain AliEmcalJet objects!", GetName(),
fDLJetArrName.Data()));
512 for(
Int_t iPart=0; iPart<nParticles; iPart++)
520 if (pt<maxPt && pt>=minPt)
522 arr.AddAt(iPart,counter);
538 if(counter==0) indexPL = -1;
539 else if(counter==1) indexPL = arr.At(0);
543 indexPL = arr.At(TMath::FloorNint(pro));
556 for (
Int_t iTracks = 0; iTracks < Ntracks; ++iTracks)
558 AliVParticle *t =
static_cast<AliVParticle*
>(
fTrackArray->At(iTracks));
559 if(!t || t->Charge()==0)
continue;
567 if (pt<maxPt && pt>=minPt)
569 arr.AddAt(iTracks,counter);
586 if(counter==0) indexDL = -1;
587 else if(counter==1) indexDL = arr.At(0);
591 indexDL = arr.At(TMath::FloorNint(pro));
596 AliDebug(2,Form(
"TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
609 if(leadingIndex<0)
return;
611 AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
615 AliDebug(2,Form(
"Found a trigger with pt = %2.2f",triggerPt));
618 if(triggerPhi<0) triggerPhi += 2*
pi;
619 Int_t nJets = jetArray->GetEntries();
620 for(
Int_t ij=0; ij<nJets; ij++)
630 if(!isBkg) fill[1] = jetPt;
631 AliDebug(10,
"Fill the histograms");
639 if(leadingIndex<0)
return;
643 AliWarning(
"Jet array is not available.");
647 AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
650 for(
Int_t ij=0; ij<nJets; ij++)
656 if(jetPt<10)
continue;
660 if(mthJetIndexEn>-1 && fraction>0.5)
677 if(!jetArray)
return -1;
680 Int_t nJets = jetArray->GetEntries();
682 for(
Int_t ij=0; ij<nJets; ij++)
685 if(!jetTmp)
continue;
686 if(TMath::Abs(jetTmp->
Eta())>1)
continue;
689 Double_t dRTmp = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
705 if(!jetArray || !jet)
return -1;
708 Int_t nJets = jetArray->GetEntries();
712 for(
Int_t ij=0; ij<nJets; ij++)
715 if(!jetTmp)
continue;
716 if(TMath::Abs(jetTmp->
Eta())>1)
continue;
721 for(
Int_t ic=0; ic<nc; ic++)
723 for(
Int_t ijc=0; ijc<nJetC; ijc++)
754 if(dPhi>2*
pi) dPhi -= 2*
pi;
755 if(dPhi<-2*
pi) dPhi += 2*
pi;
756 if(dPhi<-0.5*
pi) dPhi += 2*
pi;
757 if(dPhi>1.5*
pi) dPhi -= 2*
pi;
764 Double_t dPhi = TMath::Abs(phi1-phi2);
765 if(dPhi>2*
pi) dPhi -= 2*
pi;
766 if(dPhi>
pi) dPhi = 2*
pi - dPhi;
775 return TMath::Sqrt(dPhi*dPhi+dEta*dEta);
783 AliWarning(Form(
"Particle-level jet array is not available: %s\n",
fPLJetArrName.Data()));
788 for(
Int_t ij=0; ij<nPLJets; ij++)
796 AliDebug(5, Form(
"PL jet %d has (pt,eta,phi) = (%2.2f,%2.2f,%2.2f)",ij,jetPt,jet->
Eta(),jet->
Phi()));
802 AliWarning(Form(
"Detector-level jet array is not available: %s\n",
fDLJetArrName.Data()));
807 for(
Int_t ij=0; ij<nDLJets; ij++)
815 AliDebug(5, Form(
"DL jet %d has pt = %2.2f",ij,jetPt));
834 if(TMath::Abs(jet->
Eta())>etaCut)
return kFALSE;
843 const char *decision[2] = {
"no",
"yes"};
844 const char *TTtype[2] = {
"Single inclusive",
"Leading"};
845 printf(
"\n\n===== h-jet analysis configuration =====\n");
849 printf(
"Track phi range: %2.0f < phi < %2.0f\n",
fMinTrkPhi*TMath::RadToDeg(),
fMaxTrkPhi*TMath::RadToDeg());
850 printf(
"TT type: %s\n", TTtype[
fTTtype]);
853 printf(
"Run QA: %s\n",decision[
fRunQA]);
854 printf(
"Run particle level: %s\n",decision[
fRunPL]);
855 printf(
"Run detector level: %s\n",decision[
fRunDL]);
856 printf(
"Run h+jet: %s\n",decision[
fRunHJet]);
857 printf(
"Run matching: %s\n",decision[
fRunMatch]);
858 printf(
"=======================================\n\n");
864 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
virtual ~AliAnalysisTaskHJetEmbed()
Bool_t AcceptTrack(const AliVParticle *track)
THnSparse * fhJetPtEnMatch[kNTrig]
Double_t CalculateDPhi(const Double_t phi1, const Double_t phi2)
Double_t fMinTrkPt
Index of the trigger track in the event.
THnSparse * fhPLJetPtVsCent[kNTrig]
TH2F * fhRhoVsCent[kNTrig]
TClonesArray * fDLJetArray
Array of the embedded PYTHIA jet array on particle level.
Double_t GetDPhi(const Double_t phi1, const Double_t phi2)
THnSparse * fhDLTT[kNTrig]
Double_t GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet *jet2)
TClonesArray * fPLJetArray
Array of the found jets.
TString fRhoName
Array of the embedded PYTHIA jet array on detector level.
void UserExec(Option_t *option)
THnSparse * fhJetPhiEnMatch[kNTrig]
THnSparse * fhPLHJet[kNTrig]
TH1F * fhCentrality[kNTrig]
void UserCreateOutputObjects()
Int_t fTriggerType
Input event.
THnSparse * fhTTPtQA[kNTrig]
UShort_t GetNumberOfConstituents() const
TH1F * fhEventStat
Output list.
Int_t TrackAt(Int_t idx) const
AliNamedString * fPtHardBinName
Value of the rho parameter.
THnSparse * fhPLTT[kNTrig]
Double_t fMaxVtxZ
V0M for current event.
void RunHJet(const Double_t minPt, const Double_t maxPt)
Int_t fPtHardBin
Pt hard bin param.
THnSparse * fhHJet[kNTrig]
THnSparse * fhDLJetPtVsCent[kNTrig]
Int_t FindGeoMatchedJet(const AliEmcalJet *jet, const TClonesArray *jetArray, Double_t &dR)
Int_t FindEnergyMatchedJet(const AliEmcalJet *jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
THnSparse * fhTTPt[kNTrig]
Double_t fRhoValue
Rho parameter.
THnSparse * fhJetPhiGeoMatch[kNTrig]
AliAnalysisTaskHJetEmbed()
TClonesArray * fMCParticleArray
Double_t GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz)
Double_t fCentrality
Trigger type of the event.
THnSparse * fhJetPtGeoMatch[kNTrig]
Represent a jet reconstructed using the EMCal jet framework.
void FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg=kFALSE)
TClonesArray * fTrackArray
void RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
Bool_t IsGoodJet(const AliEmcalJet *jet)
void Terminate(Option_t *)
TString fMCParticleArrName
THnSparse * fhDLHJet[kNTrig]
TList * OpenFile(const char *fname)