2 // Detector response jet mass analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.h"
24 #include "AliEmcalParticle.h"
25 #include "AliMCEvent.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliAODMCHeader.h"
28 #include "AliMCEvent.h"
29 #include "AliAnalysisManager.h"
30 #include "AliJetContainer.h"
32 #include "AliAODEvent.h"
34 #include "AliAnalysisTaskJetMassResponseDet.h"
36 ClassImp(AliAnalysisTaskJetMassResponseDet)
38 //________________________________________________________________________
39 AliAnalysisTaskJetMassResponseDet::AliAnalysisTaskJetMassResponseDet() :
40 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMassResponseDet", kTRUE),
44 fh2PtVsMassJetPartAll(0),
45 fh2PtVsMassJetPartMatch(0),
46 fh2PtVsMassJetPartTagged(0),
47 fh2PtVsMassJetPartTaggedMatch(0),
48 fh2PtVsMassJetDetAll(0),
49 fh2PtVsMassJetDetTagged(0),
50 fh2EtaPhiMatchedDet(0),
51 fh2EtaPhiMatchedPart(0),
56 // Default constructor.
58 SetMakeGeneralHistograms(kTRUE);
61 //________________________________________________________________________
62 AliAnalysisTaskJetMassResponseDet::AliAnalysisTaskJetMassResponseDet(const char *name) :
63 AliAnalysisTaskEmcalJet(name, kTRUE),
67 fh2PtVsMassJetPartAll(0),
68 fh2PtVsMassJetPartMatch(0),
69 fh2PtVsMassJetPartTagged(0),
70 fh2PtVsMassJetPartTaggedMatch(0),
71 fh2PtVsMassJetDetAll(0),
72 fh2PtVsMassJetDetTagged(0),
73 fh2EtaPhiMatchedDet(0),
74 fh2EtaPhiMatchedPart(0),
79 // Standard constructor.
81 SetMakeGeneralHistograms(kTRUE);
84 //________________________________________________________________________
85 AliAnalysisTaskJetMassResponseDet::~AliAnalysisTaskJetMassResponseDet()
90 //________________________________________________________________________
91 void AliAnalysisTaskJetMassResponseDet::UserCreateOutputObjects()
93 // Create user output.
95 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
97 Bool_t oldStatus = TH1::AddDirectoryStatus();
98 TH1::AddDirectory(kFALSE);
100 const Int_t nBinsPt = 200;
101 const Double_t minPt = 0.;
102 const Double_t maxPt = 200.;
104 const Int_t nBinsM = 100;
105 const Double_t minM = 0.;
106 const Double_t maxM = 50.;
108 const Int_t nBinsConstEff = 40;
109 const Double_t minConstEff = 0.;
110 const Double_t maxConstEff = 2.;
112 // const Int_t nBinsConst = 26;
113 // const Double_t minConst = -5.5;
114 // const Double_t maxConst = 20.5;
116 //Binning for THnSparse
117 const Int_t nBinsSparse0 = 5;
118 const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsConstEff};
119 const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minConstEff};
120 const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxConstEff};
123 TString histName = "";
124 TString histTitle = "";
126 histName = "fh2PtVsMassJetPartAll";
127 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
128 fh2PtVsMassJetPartAll = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
129 fOutput->Add(fh2PtVsMassJetPartAll);
131 histName = "fh2PtVsMassJetPartMatch";
132 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
133 fh2PtVsMassJetPartMatch = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
134 fOutput->Add(fh2PtVsMassJetPartMatch);
136 histName = "fh2PtVsMassJetPartTagged";
137 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
138 fh2PtVsMassJetPartTagged = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
139 fOutput->Add(fh2PtVsMassJetPartTagged);
141 histName = "fh2PtVsMassJetPartTaggedMatch";
142 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
143 fh2PtVsMassJetPartTaggedMatch = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
144 fOutput->Add(fh2PtVsMassJetPartTaggedMatch);
146 histName = "fh2PtVsMassJetDetAll";
147 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
148 fh2PtVsMassJetDetAll = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
149 fOutput->Add(fh2PtVsMassJetDetAll);
151 histName = "fh2PtVsMassJetDetTagged";
152 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
153 fh2PtVsMassJetDetTagged = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
154 fOutput->Add(fh2PtVsMassJetDetTagged);
156 histName = "fh2EtaPhiMatchedDet";
157 histTitle = TString::Format("%s;#eta;#varphi",histName.Data());
158 fh2EtaPhiMatchedDet = new TH2F(histName.Data(),histTitle.Data(),100,-1.,1.,72,0.,TMath::TwoPi());
159 fOutput->Add(fh2EtaPhiMatchedDet);
161 histName = "fh2EtaPhiMatchedPart";
162 histTitle = TString::Format("%s;#eta;#varphi",histName.Data());
163 fh2EtaPhiMatchedPart = new TH2F(histName.Data(),histTitle.Data(),100,-1.,1.,72,0.,TMath::TwoPi());
164 fOutput->Add(fh2EtaPhiMatchedPart);
166 histName = "fhnMassResponse";
167 histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part};#it{N}_{const}^{det}/#it{N}_{const}^{part}",histName.Data());
168 fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
169 fOutput->Add(fhnMassResponse);
171 fh1AreaPartAll = new TH1D("fh1AreaPartAll","fh1AreaPartAll",100.,0.,1.);
172 fOutput->Add(fh1AreaPartAll);
174 fh1AreaDetAll = new TH1D("fh1AreaDetAll","fh1AreaDetAll",100.,0.,1.);
175 fOutput->Add(fh1AreaDetAll);
178 // =========== Switch on Sumw2 for all histos ===========
179 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
180 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
185 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
189 TH1::AddDirectory(oldStatus);
191 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
194 //________________________________________________________________________
195 Bool_t AliAnalysisTaskJetMassResponseDet::Run()
197 // Run analysis code here, if needed. It will be executed before FillHistograms().
201 //________________________________________________________________________
202 Bool_t AliAnalysisTaskJetMassResponseDet::FillHistograms()
206 AliJetContainer *cPart = GetJetContainer(fContainerPart);
207 AliJetContainer *cDet = GetJetContainer(fContainerDet);
208 AliEmcalJet* jPart = NULL;
209 AliEmcalJet* jDet = NULL;
211 //loop on particle level jets
214 cPart->ResetCurrentID();
215 while((jPart = cPart->GetNextAcceptJet())) {
216 fh2PtVsMassJetPartAll->Fill(jPart->Pt(),jPart->M());
217 fh1AreaPartAll->Fill(jPart->Area());
219 jDet = jPart->ClosestJet();
220 if(jDet) fh2PtVsMassJetPartMatch->Fill(jPart->Pt(),jPart->M());
221 if(jPart->GetTagStatus()<1 || !jPart->GetTaggedJet())
223 fh2PtVsMassJetPartTagged->Fill(jPart->Pt(),jPart->M());
224 if(jDet) fh2PtVsMassJetPartTaggedMatch->Fill(jPart->Pt(),jPart->M());
228 //loop on detector level jets
231 cDet->ResetCurrentID();
232 while((jDet = cDet->GetNextAcceptJet())) {
233 Double_t mjet = GetJetMass(jDet);
234 fh2PtVsMassJetDetAll->Fill(jDet->Pt(),mjet);
235 fh1AreaDetAll->Fill(jDet->Area());
237 if(jDet->GetTagStatus()>=1 && jDet->GetTaggedJet())
238 fh2PtVsMassJetDetTagged->Fill(jDet->Pt(),mjet);
240 //fill detector response
241 jPart = jDet->ClosestJet();
243 fh2EtaPhiMatchedDet->Fill(jDet->Eta(),jDet->Phi());
244 fh2EtaPhiMatchedPart->Fill(jPart->Eta(),jPart->Phi());
246 Int_t nConstPart = jPart->GetNumberOfConstituents();
247 Int_t nConstDet = jDet->GetNumberOfConstituents();
249 if(nConstPart>0) eff = (Double_t)nConstDet/((Double_t)nConstPart);
250 Double_t var[5] = {GetJetMass(jDet),jPart->M(),jDet->Pt(),jPart->Pt(),eff};
251 fhnMassResponse->Fill(var);
259 //________________________________________________________________________
260 Double_t AliAnalysisTaskJetMassResponseDet::GetJetMass(AliEmcalJet *jet) {
261 //calc subtracted jet mass
262 if(fJetMassType==kRaw)
264 else if(fJetMassType==kDeriv)
265 return jet->GetSecondOrderSubtracted();
270 //________________________________________________________________________
271 Bool_t AliAnalysisTaskJetMassResponseDet::RetrieveEventObjects() {
273 // retrieve event objects
274 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
280 //_______________________________________________________________________
281 void AliAnalysisTaskJetMassResponseDet::Terminate(Option_t *)
283 // Called once at the end of the analysis.