AliPhysics  8d00e07 (8d00e07)
AliAnalysisTaskEmcalTriggerJets.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
27 #include <THistManager.h>
28 #include <TLinearBinning.h>
29 #include <TMath.h>
30 #include <TString.h>
31 
32 #include "AliAnalysisManager.h"
33 #include "AliAODTrack.h"
34 #include "AliClusterContainer.h"
35 #include "AliEmcalJet.h"
36 #include "AliInputEventHandler.h"
37 #include "AliJetContainer.h"
38 #include "AliLog.h"
39 #include "AliPIDResponse.h"
40 #include "AliTrackContainer.h"
42 #include "AliVCluster.h"
43 #include "AliVEvent.h"
44 #include "AliVTrack.h"
45 
46 #include <array>
47 #include <iostream>
48 
52 
53 namespace EmcalTriggerJets {
54 
55 AliAnalysisTaskEmcalTriggerJets::AliAnalysisTaskEmcalTriggerJets():
57  fPIDResponse(nullptr),
58  fHistos(nullptr)
59 {
60 
61 }
62 
64  AliAnalysisTaskEmcalJet(name, true),
67 {
68 
69 }
70 
72  delete fHistos;
73 }
74 
77 
78  const std::array<TString, 5> kEmcalTriggers = {"INT7", "EJ1", "EJ2", "DJ1", "DJ2"};
79  const int kNJetPtBins = 9;
80  const int kNJetRadiusBins = 7;
81 
82  TLinearBinning jetptbinning(9, 20, 200), pbinning(300, 0., 30.), dEdxbinning(600, 0., 600.), massbinning(400., 0., 4.), eopbinning(150, 0., 1.5), radiusBinning(10, 0., 1.);
83  const TBinning *binningPID[5] {&jetptbinning, &pbinning, &dEdxbinning, &massbinning, &eopbinning},
84  *binningAssociate[6] = {&jetptbinning, &pbinning, &radiusBinning, &dEdxbinning, &massbinning, &eopbinning};
85  const std::array<int, 2> kJetRadii = {2, 4};
86  const std::array<TString, 3> kJetTypes = {"Charged", "Full", "Neutral"};
87  const std::array<TString, 2> kDetectors = {"EMCAL", "DCAL"};
88  const std::array<TString, 2> kConstituentType = {"Charged", "Neutral"};
89  fHistos = new THistManager("EmcalJetHistos");
90  for(auto t : kEmcalTriggers){
91  fHistos->CreateTH1("hEventCount" + t, "Event counter for trigger " + t, 1., 0.5, 1.5);
92  for(auto det: kDetectors){
93  for(auto jt : kJetTypes) {
94  for(auto radius : kJetRadii){
95  fHistos->CreateTH1("hPtRaw" + jt + "JetR" + Form("%02d", radius) + det + t,
96  "Raw pt spectrum for " + jt + " jets with R=" + Form("%.1f", double(radius)/10.) + " in " + det +" for trigger " + t,
97  200, 0., 200);
98  if(jt == "Full") {
99  // Neutral energy fraction vs. jet pt
100  fHistos->CreateTH2("hNefFullJetR" + TString::Format("%02d", radius) + det + t,
101  "Neutral energy fraction vs. jet pt for R=" + TString::Format("%.1f", double(radius)/10.) + " full jets in " + det + " for trigger " + t,
102  200, 0., 200., 100, 0., 1.);
103  // pt, z and pt_rel of the leading particle (cluster and track) vs jet pt
104  for(auto constituent : kConstituentType){
105  fHistos->CreateTH2("hLeading" + constituent + "PtFullJetR" + TString::Format("%02d", radius) + det + t,
106  "p{t_const} vs. p_{t, jet} for leading " + constituent + " constituents in full jets with R=" + TString::Format("%.1f", double(radius)/10.) + " in " + det + "for trigger " + t,
107  200, 0., 200., 200, 0., 200.);
108  fHistos->CreateTH2("hLeading" + constituent + "ZFullJetR" + TString::Format("%02d", radius) + det + t,
109  "z vs. p_{t, jet} for leading " + constituent + " constituents in full jets with R=" + TString::Format("%.1f", double(radius)/10.) + " in " + det + " for trigger " + t,
110  200, 0., 200., 100, 0., 1.);
111  fHistos->CreateTH2("hLeading" + constituent + "PtrelFullJetR" + TString::Format("%02d", radius) + det + t,
112  "p_{t,rel} vs. p_{t, jet} for leading " + constituent + " constituents in full jets with R=" + TString::Format("%.1f", double(radius)/10.) + " in " + det + " for trigger " + t,
113  200., 0., 200., 100, 0., 1.);
114  }
115  fHistos->CreateTHnSparse("hPIDConstituentFullJet" + TString::Format("R%02d", radius) + det + t, TString::Format("PID for full jet constituents for R=%0.2f in %s for trigger %s", double(radius)/10., det.Data(), t.Data()), 5, binningPID);
116  fHistos->CreateTHnSparse("hPIDLeadingFullJet" + TString::Format("R%02d", radius) + det + t, TString::Format("PID for full jet leading constituents for R=%0.2f in %s for trigger %s", double(radius)/10., det.Data(), t.Data()), 5, binningPID);
117  }
118  }
119  }
120  }
121  }
122  for(auto h : *(fHistos->GetListOfHistograms())){
123  fOutput->Add(h);
124  }
125  PostData(1, fOutput);
126 }
127 
129  fPIDResponse = fInputHandler->GetPIDResponse();
130  if(!fPIDResponse){
131  AliErrorStream() << "PID Response not available - PID plots will not be filled" << std::endl;
132  }
133 }
134 
136  std::vector<TString> triggers, kEmcalTriggers = {"EJ1", "EJ2", "DJ1", "DJ2"};
137  if(fInputHandler->IsEventSelected() & AliVEvent::kINT7) triggers.push_back("INT7");
138  if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE){
139  TString fired = fInputEvent->GetFiredTriggerClasses();
140  for(auto e : kEmcalTriggers){
141  if(fired.Contains(e)) triggers.push_back(e);
142  }
143  }
144  if(!triggers.size()) return false;
145 
146  for(auto t : triggers) fHistos->FillTH1("hEventCount" + t, 1);
147  std::vector<TString> jettypes = {"Full", "Charged", "Neutral"}, detectors = {"EMCAL", "DCAL"}, radii = {"R02", "R04"};
148  for(auto jt : jettypes) {
149  for(auto det : detectors){
150  for(auto r : radii) {
151  double radius = double(TString(r).ReplaceAll("R0","").Atoi())/10.;
152  TString namejcont = jt + "Jets" + r + det,
153  rawhistnamebase = "hPtRaw" + jt + "Jet" + r + det,
154  tagForLeading = "FullJet" + r + det;
155  AliJetContainer *c = this->GetJetContainer(namejcont);
156  if(!c) AliErrorStream() << "Not found jet container " << namejcont << std::endl;
157  bool doPID = fPIDResponse && (jt == "Full") && (r == "R04");
158  for(auto j : c->accepted()){
159  TLorentzVector jetvec(j->Px(), j->Py(), j->Pt(), j->E());
160  // get the leading particle
161  AliVCluster *leadingClust = (jt == "Full") ? j->GetLeadingCluster(this->GetClusterContainer("caloClusters")->GetArray()) : nullptr;
162  AliVTrack *leadingTrack = (jt == "Full") ? static_cast<AliVTrack *>(j->GetLeadingTrack(this->GetTrackContainer("tracks")->GetArray())) : nullptr;
163  TLorentzVector clustervec;
164  TVector3 trackvec;
165  if(leadingClust){
166  leadingClust->GetMomentum(clustervec, fVertex);
167  }
168  if(leadingTrack) {
169  trackvec.SetXYZ(leadingTrack->Px(), leadingTrack->Py(), leadingTrack->Pz());
170  }
171  double absJetPt = TMath::Abs(j->Pt());
172  for(auto t : triggers) {
173  fHistos->FillTH1(rawhistnamebase + t, absJetPt);
174  if(jt == "Full") {
175  AliDebugStream(1) << "Filling full jet leading histograms, constituents c[" << j->GetNumberOfTracks() << "], n[" << j->GetNumberOfClusters() << "]" << std::endl;
176  FillJetPIDPlots(j, radius, t, det);
177  fHistos->FillTH2("hNefFullJet" + r + det + t, absJetPt, j->NEF());
178  // fill also histograms for the neutral energy and the leading particles
179  if(leadingClust){
180  fHistos->FillTH2("hLeadingNeutralPt" + tagForLeading + t, absJetPt, TMath::Abs(clustervec.Pt()));
181  fHistos->FillTH2("hLeadingNeutralZ" + tagForLeading + t, absJetPt, j->GetZ(clustervec.Px(), clustervec.Py(), clustervec.Pz()));
182  fHistos->FillTH2("hLeadingNeutralPtrel" + tagForLeading + t, absJetPt, clustervec.Pt(jetvec.Vect()));
183  } else {
184  AliDebugStream(2) << "No leading cluster found" << std::endl;
185  }
186  if(leadingTrack){
187  fHistos->FillTH2("hLeadingChargedPt" + tagForLeading + t, absJetPt, TMath::Abs(leadingTrack->Pt()));
188  fHistos->FillTH2("hLeadingChargedZ" + tagForLeading + t, absJetPt, j->GetZ(leadingTrack));
189  fHistos->FillTH2("hLeadingChargedPtrel" + tagForLeading + t, absJetPt, trackvec.Pt(jetvec.Vect()));
190  AliDebugStream(2) << "Filling full jet leading PID histograms" << std::endl;
191  this->FillJetPIDPlotsLeading(leadingTrack, j->Pt(), radius, t, det);
192  } else {
193  AliDebugStream(2) << "No leading track found" << std::endl;
194  }
195  AliDebugStream(1) << "Filling full jet leading constituent histograms done" << std::endl;
196  }
197  }
198  }
199  }
200  }
201  }
202  return true;
203 }
204 
205 void AliAnalysisTaskEmcalTriggerJets::FillJetPIDPlots(const AliEmcalJet *jet, double radius, const char *trigger, const char *detector){
206  if(TMath::Abs(jet->Pt()) < 20 || TMath::Abs(jet->Pt()) > 200) return;
207  TString histname = TString::Format("hPIDConstituentFullJetR%02d%s%s", int(radius * 10), detector, trigger);
208  AliTrackContainer *tc = GetTrackContainer("tracks");
209 
210  for(int icharged = 0; icharged < jet->GetNumberOfTracks(); icharged++){
211  AliVTrack *constituent = static_cast<AliVTrack *>(jet->TrackAt(icharged, tc->GetArray()));
212  // Select only constituents with sufficient PID information in both TPC and TOF
213  if(constituent->GetTPCsignalN() < 30) continue;
214  if(constituent->GetTPCSharedMapPtr()->CountBits(0) > 70) continue;
215  if(!((constituent->GetStatus() & AliVTrack::kTOFout) && (constituent->GetStatus() & AliVTrack::kTIME))) continue;
216  Double_t trtime = (constituent->GetTOFsignal() - fPIDResponse->GetTOFResponse().GetTimeZero()) * 1e-12;
217  Double_t v = constituent->GetIntegratedLength()/(100. * trtime);
218  Double_t beta = v / TMath::C(), gamma = 1 / TMath::Sqrt(1-beta*beta), mtof = constituent->P() / (beta*gamma);
219  Double_t eop = -1.;
220  if(constituent->GetEMCALcluster() >= 0) {
221  AliVCluster *matched = static_cast<AliVCluster *>((*GetClusterContainer("caloClusters"))[constituent->GetEMCALcluster()]);
222  if(matched) eop = TMath::Abs(matched->GetNonLinCorrEnergy()/constituent->P());
223  }
224  Double_t datapoint[5] = {TMath::Abs(jet->Pt()), TMath::Abs(constituent->P()), constituent->GetTPCsignal(), mtof, eop};
225  fHistos->FillTHnSparse(histname, datapoint);
226  }
227 }
228 
229 void AliAnalysisTaskEmcalTriggerJets::FillJetPIDPlotsLeading(const AliVTrack *leading, double ptjet, double radius, const char *trigger, const char *detector){
230  if(ptjet < 20 || ptjet > 200) return;
231 
232  TString histname = TString::Format("hPIDLeadingFullJetR%02d%s%s", int(radius*10.), detector, trigger);
233  if(leading->GetTPCsignalN() < 30) return;
234  if(!((leading->GetStatus() & AliVTrack::kTOFout) && (leading->GetStatus() & AliVTrack::kTIME))) return;
235  if(leading->GetTPCSharedMapPtr()->CountBits(0) > 70) return;
236  Double_t trtime = (leading->GetTOFsignal() - fPIDResponse->GetTOFResponse().GetTimeZero()) * 1e-12;
237  Double_t v = leading->GetIntegratedLength()/(100. * trtime);
238  Double_t beta = v / TMath::C(), gamma = 1 / TMath::Sqrt(1-beta*beta), mtof = leading->P() / (beta*gamma);
239  Double_t eop = -1.;
240  if(leading->GetEMCALcluster() >= 0){
241  AliVCluster *matched = static_cast<AliVCluster *>((*GetClusterContainer("caloClusters"))[leading->GetEMCALcluster()]);
242  if(matched) eop = TMath::Abs(matched->GetNonLinCorrEnergy()/leading->P());
243  }
244  Double_t datapoint[5] = {ptjet, TMath::Abs(leading->P()), leading->GetTPCsignal(), mtof, eop};
245  fHistos->FillTHnSparse(histname, datapoint);
246 }
247 
249  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
250 
252  mgr->AddTask(task);
253 
254  // Adding containers for clusters and tracks
255  AliClusterContainer *clustercont = task->AddClusterContainer("caloClusters");
256  clustercont->SetMinE(0.3);
257  AliTrackContainer *trackcont = task->AddTrackContainer("tracks");
258  trackcont->SetMinPt(0.15);
259 
260  // Adding Jet containers
261  // - Using jet radii 0.2 and 0.4
262  // - Splitting between EMCAL and DCAL side
263  // - processing full, charged and neutral jets
264 
265  // Full jets, R=0.2, EMCAL
266  AliJetContainer *cont = task->AddJetContainer(
270  0.2,
272  trackcont, clustercont
273  );
274  cont->SetName("FullJetsR02EMCAL");
275 
276  // Full jets, R=0.4, EMCAL
277  cont = task->AddJetContainer(
281  0.4,
283  trackcont, clustercont
284  );
285  cont->SetName("FullJetsR04EMCAL");
286 
287  // Full jets, R=0.2, DCAL
288  cont = task->AddJetContainer(
292  0.2,
294  trackcont, clustercont
295  );
296  cont->SetName("FullJetsR02DCAL");
297 
298  // Full jets, R=0.4, EMCAL
299  cont = task->AddJetContainer(
303  0.4,
305  trackcont, clustercont
306  );
307  cont->SetName("FullJetsR04DCAL");
308 
309 
310  // Charged jets, R=0.2, EMCAL
311  cont = task->AddJetContainer(
315  0.2,
317  trackcont, nullptr
318  );
319  cont->SetName("ChargedJetsR02EMCAL");
320 
321  // Charged jets, R=0.4, EMCAL
322  cont = task->AddJetContainer(
326  0.4,
328  trackcont, nullptr
329  );
330  cont->SetName("ChargedJetsR04EMCAL");
331 
332  // Charged jets, R=0.2, DCAL
333  cont = task->AddJetContainer(
337  0.2,
339  trackcont, nullptr
340  );
341  cont->SetName("ChargedJetsR02DCAL");
342 
343  // Charged jets, R=0.4, DCAL
344  cont = task->AddJetContainer(
348  0.4,
350  trackcont, nullptr
351  );
352  cont->SetName("ChargedJetsR04DCAL");
353 
354 
355  // Neutral jets, R=0.2, EMCAL
356  cont = task->AddJetContainer(
360  0.2,
362  nullptr, clustercont
363  );
364  cont->SetName("NeutralJetsR02EMCAL");
365 
366  // Neutral jets, R=0.4, EMCAL
367  cont = task->AddJetContainer(
371  0.4,
373  nullptr, clustercont
374  );
375  cont->SetName("NeutralJetsR04EMCAL");
376 
377  // Neutral jets, R=0.2, DCAL
378  cont = task->AddJetContainer(
382  0.2,
384  nullptr, clustercont
385  );
386  cont->SetName("NeutralJetsR02DCAL");
387 
388  // Neutral jets, R=0.4, DCAL
389  cont = task->AddJetContainer(
393  0.4,
395  nullptr, clustercont
396  );
397  cont->SetName("NeutralJetsR04DCAL");
398 
399  // Connect Input / Output containers
400  TString outfilename = mgr->GetCommonFileName();
401  outfilename += ":EmcalTriggerJets";
402  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
403  mgr->ConnectOutput(task, 1, mgr->CreateContainer("HistsEmcalTriggerJets", TList::Class(), AliAnalysisManager::kOutputContainer, outfilename));
404 
405  return task;
406 }
407 
408 } /* namespace EmcalTriggerJets */
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
AliJetContainer * GetJetContainer(Int_t i=0) const
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
TCanvas * c
Definition: TestFitELoss.C:172
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Interface for binnings used by the histogram handler.
Definition: TBinning.h:21
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
void FillJetPIDPlots(const AliEmcalJet *jet, double radius, const char *trigger, const char *detector)
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
virtual void UserExecOnce()
Task initializations handled in user tasks.
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalList * fOutput
!output list
Double_t fVertex[3]
!event vertex
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Container class for histograms.
Definition: THistManager.h:99
void UserCreateOutputObjects()
Main initialization function on the worker.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:72
const AliJetIterableContainer accepted() const
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
void FillJetPIDPlotsLeading(const AliVTrack *leading, double ptjet, double radius, const char *trigger, const char *detector)
Container structure for EMCAL clusters.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
virtual bool Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
static AliAnalysisTaskEmcalTriggerJets * AddTaskEmcalTriggerJets(const char *name)