AliPhysics  cda3415 (cda3415)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetMassResponseDet.cxx
Go to the documentation of this file.
1 //
2 // Detector response jet mass analysis task.
3 //
4 // Author: M.Verweij
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TList.h>
12 #include <TLorentzVector.h>
13 #include <TProfile.h>
14 #include <TChain.h>
15 #include <TSystem.h>
16 #include <TFile.h>
17 #include <TKey.h>
18 
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.h"
23 #include "AliLog.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"
31 
32 #include "AliAODEvent.h"
33 
35 
37 
38 //________________________________________________________________________
41  fContainerPart(0),
42  fContainerDet(0),
43  fJetMassType(kRaw),
44  fh2PtVsMassJetPartAll(0),
45  fh2PtVsMassJetPartMatch(0),
46  fh2PtVsMassJetPartTagged(0),
47  fh2PtVsMassJetPartTaggedMatch(0),
48  fh2PtVsMassJetDetAll(0),
49  fh2PtVsMassJetDetTagged(0),
50  fh2EtaPhiMatchedDet(0),
51  fh2EtaPhiMatchedPart(0),
52  fhnMassResponse(0),
53  fh1AreaPartAll(0),
54  fh1AreaDetAll(0)
55 {
56  // Default constructor.
57 
58  SetMakeGeneralHistograms(kTRUE);
59 }
60 
61 //________________________________________________________________________
63  AliAnalysisTaskEmcalJet(name, kTRUE),
64  fContainerPart(0),
65  fContainerDet(0),
66  fJetMassType(kRaw),
67  fh2PtVsMassJetPartAll(0),
68  fh2PtVsMassJetPartMatch(0),
69  fh2PtVsMassJetPartTagged(0),
70  fh2PtVsMassJetPartTaggedMatch(0),
71  fh2PtVsMassJetDetAll(0),
72  fh2PtVsMassJetDetTagged(0),
73  fh2EtaPhiMatchedDet(0),
74  fh2EtaPhiMatchedPart(0),
75  fhnMassResponse(0),
76  fh1AreaPartAll(0),
77  fh1AreaDetAll(0)
78 {
79  // Standard constructor.
80 
82 }
83 
84 //________________________________________________________________________
86 {
87  // Destructor.
88 }
89 
90 //________________________________________________________________________
92 {
93  // Create user output.
94 
96 
97  Bool_t oldStatus = TH1::AddDirectoryStatus();
98  TH1::AddDirectory(kFALSE);
99 
100  const Int_t nBinsPt = 200;
101  const Double_t minPt = 0.;
102  const Double_t maxPt = 200.;
103 
104  const Int_t nBinsM = 200;
105  const Double_t minM = 0.;
106  const Double_t maxM = 50.;
107 
108  const Int_t nBinsConstEff = 40;
109  const Double_t minConstEff = 0.;
110  const Double_t maxConstEff = 2.;
111 
112  // const Int_t nBinsConst = 26;
113  // const Double_t minConst = -5.5;
114  // const Double_t maxConst = 20.5;
115 
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};
121 
122  //Create histograms
123  TString histName = "";
124  TString histTitle = "";
125 
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);
130 
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);
135 
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);
140 
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);
145 
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);
150 
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);
155 
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);
160 
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);
165 
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);
170 
171  fh1AreaPartAll = new TH1D("fh1AreaPartAll","fh1AreaPartAll",100.,0.,1.);
172  fOutput->Add(fh1AreaPartAll);
173 
174  fh1AreaDetAll = new TH1D("fh1AreaDetAll","fh1AreaDetAll",100.,0.,1.);
175  fOutput->Add(fh1AreaDetAll);
176 
177 
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));
181  if (h1){
182  h1->Sumw2();
183  continue;
184  }
185  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
186  if(hn)hn->Sumw2();
187  }
188 
189  TH1::AddDirectory(oldStatus);
190 
191  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
192 }
193 
194 //________________________________________________________________________
196 {
197  // Run analysis code here, if needed. It will be executed before FillHistograms().
198  return kTRUE;
199 }
200 
201 //________________________________________________________________________
203 {
204  // Fill histograms.
205 
208  AliEmcalJet* jPart = NULL;
209  AliEmcalJet* jDet = NULL;
210 
211  //loop on particle level jets
212  Int_t nAccPart = 0;
213  if(cPart) {
214  cPart->ResetCurrentID();
215  while((jPart = cPart->GetNextAcceptJet())) {
216  fh2PtVsMassJetPartAll->Fill(jPart->Pt(),jPart->M());
217  fh1AreaPartAll->Fill(jPart->Area());
218  nAccPart++;
219  jDet = jPart->ClosestJet();
220  if(jDet) fh2PtVsMassJetPartMatch->Fill(jPart->Pt(),jPart->M());
221  if(jPart->GetTagStatus()<1 || !jPart->GetTaggedJet())
222  continue;
223  fh2PtVsMassJetPartTagged->Fill(jPart->Pt(),jPart->M());
224  if(jDet) fh2PtVsMassJetPartTaggedMatch->Fill(jPart->Pt(),jPart->M());
225  }
226  }
227 
228  //loop on detector level jets
229  Int_t nAccDet = 0;
230  if(cDet) {
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());
236  nAccDet++;
237  if(jDet->GetTagStatus()>=1 && jDet->GetTaggedJet())
238  fh2PtVsMassJetDetTagged->Fill(jDet->Pt(),mjet);
239 
240  //fill detector response
241  jPart = jDet->ClosestJet();
242  if(jPart) {
243  fh2EtaPhiMatchedDet->Fill(jDet->Eta(),jDet->Phi());
244  fh2EtaPhiMatchedPart->Fill(jPart->Eta(),jPart->Phi());
245 
246  Int_t nConstPart = jPart->GetNumberOfConstituents();
247  Int_t nConstDet = jDet->GetNumberOfConstituents();
248  Double_t eff = -1.;
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);
252  }
253  }
254  }
255 
256  return kTRUE;
257 }
258 
259 //________________________________________________________________________
261  //calc subtracted jet mass
262  if(fJetMassType==kRaw)
263  return jet->M();
264  else if(fJetMassType==kDeriv)
266 
267  return 0;
268 }
269 
270 //________________________________________________________________________
272  //
273  // retrieve event objects
275  return kFALSE;
276 
277  return kTRUE;
278 }
279 
280 //_______________________________________________________________________
282 {
283  // Called once at the end of the analysis.
284 }
285 
Double_t Area() const
Definition: AliEmcalJet.h:117
double Double_t
Definition: External.C:58
AliEmcalJet * GetTaggedJet() const
Definition: AliEmcalJet.h:223
Definition: External.C:236
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:213
AliJetContainer * GetJetContainer(Int_t i=0) const
Int_t GetTagStatus() const
Definition: AliEmcalJet.h:224
Double_t Eta() const
Definition: AliEmcalJet.h:108
Double_t Phi() const
Definition: AliEmcalJet.h:104
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:127
TH2F * fh2EtaPhiMatchedDet
pT vs mass of tagged detector level jets
THnSparse * fhnMassResponse
eta,phi of matched particle level jets
TH2F * fh2EtaPhiMatchedPart
eta,phi of matched detector level jets
TH2F * fh2PtVsMassJetPartTaggedMatch
pT vs mass of tagged particle level jets
int Int_t
Definition: External.C:63
Definition: External.C:212
ClassImp(AliAnalysisTaskJetMassResponseDet) AliAnalysisTaskJetMassResponseDet
TH2F * fh2PtVsMassJetDetTagged
pT vs mass of all detector level jets
AliEmcalJet * GetNextAcceptJet()
Double_t Pt() const
Definition: AliEmcalJet.h:96
AliEmcalList * fOutput
!output list
TH2F * fh2PtVsMassJetPartMatch
pT vs mass of all particle level jets
TH2F * fh2PtVsMassJetPartTagged
pT vs mass of all particle level jets matched to a detector level jet
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:247
TH2F * fh2PtVsMassJetDetAll
pT vs mass of tagged particle level jets matched to a detector level jet
TH1D * fh1AreaDetAll
area of all particle level jets
Double_t M() const
Definition: AliEmcalJet.h:107
Container for jet within the EMCAL jet framework.
Definition: External.C:196