AliPhysics  e6c8d43 (e6c8d43)
AliAnalysisTaskEmcalJetConstituentQA.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 <algorithm>
28 #include <array>
29 #include <cmath>
30 #include <iomanip>
31 #include <iostream>
32 #include <sstream>
33 #include <string>
34 #include <THistManager.h>
35 #include <TLinearBinning.h>
36 #include <TLorentzVector.h>
37 #include <TVector3.h>
38 
39 #include "AliAnalysisManager.h"
41 #include "AliAODInputHandler.h"
42 #include "AliClusterContainer.h"
45 #include "AliEmcalJet.h"
46 #include "AliInputEventHandler.h"
47 #include "AliJetContainer.h"
48 #include "AliLog.h"
49 #include "AliTrackContainer.h"
50 #include "AliVCluster.h"
51 #include "AliVEvent.h"
52 #include "AliVTrack.h"
53 
57 
58 using namespace EmcalTriggerJets;
59 
62  fHistos(nullptr),
63  fJetType(AliJetContainer::kFullJet),
64  fNameTrackContainer(""),
65  fNameClusterContainer(""),
66  fTriggerSelectionString(""),
67  fUseTriggerSelection(kFALSE),
68  fNameTriggerDecisionContainer("EmcalTriggerDecision")
69 {
70  this->SetUseAliAnaUtils(true);
71 }
72 
74  AliAnalysisTaskEmcalJet(name, true),
76  fJetType(AliJetContainer::kFullJet),
80  fUseTriggerSelection(kFALSE),
81  fNameTriggerDecisionContainer("EmcalTriggerDecision")
82 {
83  this->SetUseAliAnaUtils(true);
84 }
85 
87  if(fHistos) delete fHistos;
88 }
89 
92 
93  TLinearBinning binningz(50, 0., 1), multbinning(51, -0.5, 50.5), binningnef(50, 0., 1.), binningR(50, 0., 0.5), binningptconst(200, 0., 200.), binningptjet(20, 0., 200.),
94  binningNCell(101, -0.5, 100.5), binningFracCellLeading(100, 0., 1.), binningM02(100, 0., 1.), etabinning(100, -0.8, 0.8), phibinning(100, 0., TMath::TwoPi());
95 
96  const TBinning *jetbinning[4] = {&binningptjet, &binningnef, &multbinning, &multbinning},
97  *chargedbinning[7] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningz, &binningR},
98  *neutralbinning[8] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningptconst, &binningz, &binningR},
99  *binningHighZClusters[7] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningNCell, &binningFracCellLeading, &binningM02},
100  *leadingbinning[5] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningR},
101  *leadingjetvecbinning[4] = {&binningptjet, &etabinning, &phibinning, &binningptjet};
102 
103  fHistos = new THistManager(Form("histos_%s", GetName()));
104  for(auto c : fNamesJetContainers){
105  auto contname = dynamic_cast<TObjString *>(c);
106  if(!contname) continue;
107  fHistos->CreateTHnSparse(Form("hJetCounter%s", contname->String().Data()), Form("jet counter for jets %s", contname->String().Data()), 4, jetbinning);
108  fHistos->CreateTHnSparse(Form("hPtEtaPhiELeadingJet%s", contname->String().Data()), Form("Momemtum vector of leading jets %s", contname->String().Data()), 4, leadingjetvecbinning);
110  fHistos->CreateTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), Form("charged constituents in jets %s", contname->String().Data()), 7, chargedbinning);
111  fHistos->CreateTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), Form("leading charged constituent in jets %s", contname->String().Data()), 5, leadingbinning);
112  fHistos->CreateTHnSparse(Form("hLeadingJetLeadingTrack%s", contname->String().Data()), Form("leading charged constituent in jets %s", contname->String().Data()), 5, leadingbinning);
113  }
115  fHistos->CreateTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), Form("neutral constituents in jets %s", contname->String().Data()), 8, neutralbinning);
116  fHistos->CreateTHnSparse(Form("hHighZClusters%s", contname->String().Data()), "Properties of high-z clusters", 7, binningHighZClusters);
117  fHistos->CreateTHnSparse(Form("hLeadingCluster%s", contname->String().Data()), Form("leading neutral constituent in jets %s", contname->String().Data()), 5, leadingbinning);
118  fHistos->CreateTHnSparse(Form("hLeadingJetLeadingCluster%s", contname->String().Data()), Form("leading neutral constituent in jets %s", contname->String().Data()), 5, leadingbinning);
119  }
120  }
121 
122  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
123  PostData(1, fOutput);
124 }
125 
128  if(!tracks) tracks = GetParticleContainer(fNameTrackContainer);
129  const auto clusters = GetClusterContainer(fNameClusterContainer);
130  if(fNameTrackContainer.Length() && !tracks){
131  AliErrorStream() << "Track container " << fNameTrackContainer << " required but missing ..." << std::endl;
132  return kFALSE;
133  }
134  if(fNameClusterContainer.Length() && !clusters){
135  AliErrorStream() << "Cluster container " << fNameClusterContainer << " required but missing ..." << std::endl;
136  return kFALSE;
137  }
138 
139  // Event selection
140  AliDebugStream(1) << "Trigger selection string: " << fTriggerSelectionString << ", fired trigger classes: " << fInputEvent->GetFiredTriggerClasses() << std::endl;
141  if(fTriggerSelectionString.Contains("INT7")){
142  // INT7 trigger
143  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
144  } else if(fTriggerSelectionString.Contains("EJ")){
145  auto triggerclass = fTriggerSelectionString(fTriggerSelectionString.Index("EJ"),3); // cleanup trigger string from possible tags
146  AliDebugStream(1) << "Inspecting trigger class " << triggerclass << std::endl;
147  // EMCAL JET trigger
148  if(!fMCEvent){ // Running on Data
149  if(!(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)) return false;
150  if(!fInputEvent->GetFiredTriggerClasses().Contains(triggerclass)) return false;
151  }
152 
154  auto triggerdecisions = dynamic_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer.Data()));
155  if(!triggerdecisions) {
156  AliErrorStream() << "No offline trigger selection available" << std::endl;
157  return false;
158  }
159  else if(!triggerdecisions->IsEventSelected(triggerclass.Data())) return false;
160  }
161  } else return false;
162 
163  AliDebugStream(1) << "Event is selected" << std::endl;
164 
165  for(auto jc : fNamesJetContainers){
166  auto contname = dynamic_cast<TObjString *>(jc);
167  if(!contname) {
168  AliErrorStream() << "Non-string object in the list of jet container names" << std::endl;
169  continue;
170  }
171  const auto jetcont = GetJetContainer(contname->String().Data());
172  if(!jetcont){
173  AliErrorStream() << "Jet container with name " << contname->String() << " not found in the list of jet containers" << std::endl;
174  continue;
175  }
176  AliDebugStream(2) << "Reading " << jetcont->GetArray()->GetName() << std::endl;
177 
178  AliEmcalJet *leadingjet(nullptr);
179  for(auto jet : jetcont->accepted()){
180  if(!leadingjet || jet->Pt() > leadingjet->Pt()) leadingjet = jet;
181  AliDebugStream(3) << "Next accepted jet, found " << jet->GetNumberOfTracks() << " tracks and " << jet->GetNumberOfClusters() << " clusters." << std::endl;
182  Double_t pointjet[4] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters())},
183  pointcharged[7] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1.},
184  pointneutral[8] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1., -1.},
185  pointHighZCluster[7] = {std::abs(jet->Pt()), jet->NEF(), -1., -1., -1., -1., -1.};
186  fHistos->FillTHnSparse(Form("hJetCounter%s", contname->String().Data()), pointjet);
187  TVector3 jetvec{jet->Px(), jet->Py(), jet->Pz()};
189  for(decltype(jet->GetNumberOfTracks()) itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
190  const auto trk = jet->TrackAt(itrk, tracks->GetArray());
191  if(!trk) continue;
192  if(trk->Charge()){
193  pointcharged[4] = std::abs(trk->Pt());
194  pointcharged[5] = std::abs(jet->GetZ(trk));
195  pointcharged[6] = jet->DeltaR(trk);
196  fHistos->FillTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), pointcharged);
197  } else {
198  // particle level jets
199  pointneutral[4] = pointneutral[5] = std::abs(trk->E());
200  pointneutral[6] = std::abs(jet->GetZ(trk));
201  pointneutral[7] = jet->DeltaR(trk);
202  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
203  }
204  }
205  // Leading track
206  auto leadingtrack = jet->GetLeadingTrack(tracks->GetArray());
207  if(leadingtrack) {
208  double ltrackpoint[5] = {std::abs(jet->Pt()), jet->NEF(), std::abs(leadingtrack->Pt()), jet->GetZ(leadingtrack), jet->DeltaR(leadingtrack)};
209  fHistos->FillTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), ltrackpoint);
210  }
211  }
213  for(decltype(jet->GetNumberOfClusters()) icl = 0; icl < jet->GetNumberOfClusters(); icl++){
214  const auto clust = jet->ClusterAt(icl, clusters->GetArray());
215  if(!clust) continue;
216  TLorentzVector ptvec;
217  clust->GetMomentum(ptvec, this->fVertex, AliVCluster::kHadCorr);
218  pointneutral[4] = std::abs(clust->GetHadCorrEnergy());
219  pointneutral[5] = std::abs(clust->GetNonLinCorrEnergy());
220  pointneutral[6] = jet->GetZ(ptvec.Px(), ptvec.Py(), ptvec.Pz());
221  pointneutral[7] = jetvec.DeltaR(ptvec.Vect());
222  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
223 
224  if(pointneutral[6] > 0.95) {
225  pointHighZCluster[2] = pointneutral[4];
226  pointHighZCluster[3] = pointneutral[6];
227  pointHighZCluster[4] = clust->GetNCells();
228  pointHighZCluster[5] = *std::max_element(clust->GetCellsAmplitudeFraction(), clust->GetCellsAmplitudeFraction()+clust->GetNCells());
229  pointHighZCluster[6] = clust->GetM02();
230  fHistos->FillTHnSparse(Form("hHighZClusters%s", contname->String().Data()), pointHighZCluster);
231  }
232  }
233  // Leading cluster
234  auto leadingcluster = jet->GetLeadingCluster(clusters->GetArray());
235  if(leadingcluster){
236  TLorentzVector pvect;
237  leadingcluster->GetMomentum(pvect, fVertex);
238  double lclusterpoint[5] = {std::abs(jet->Pt()), jet->NEF(), std::abs(pvect.Pt()), jet->GetZ(pvect.Px(), pvect.Py(), pvect.Pz()), jetvec.DeltaR(pvect.Vect())};
239  fHistos->FillTHnSparse(Form("hLeadingCluster%s", contname->String().Data()), lclusterpoint);
240  }
241  }
242  }
243 
244  // Look at leading particles in the leading jet
245  double leadingvec[4] = {0., 0., 0., 0.};
246  if(leadingjet) {
247  TVector3 leadingjetvec{leadingjet->Px(), leadingjet->Py(), leadingjet->Pz()};
248  leadingvec[0] = std::abs(leadingjet->Pt());
249  leadingvec[1] = leadingjet->Eta();
250  leadingvec[2] = leadingjet->Phi();
251  if(leadingvec[2] < 0) leadingvec[2] += TMath::TwoPi();
252  leadingvec[3] = leadingjet->E();
254  auto leadingtrack = leadingjet->GetLeadingTrack(tracks->GetArray());
255  if(leadingtrack) {
256  double ltrackpoint[5] = {std::abs(leadingjet->Pt()), leadingjet->NEF(), std::abs(leadingtrack->Pt()), leadingjet->GetZ(leadingtrack), leadingjet->DeltaR(leadingtrack)};
257  fHistos->FillTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), ltrackpoint);
258  }
259  }
261  auto leadingcluster = leadingjet->GetLeadingCluster(clusters->GetArray());
262  if(leadingcluster){
263  TLorentzVector pvect;
264  leadingcluster->GetMomentum(pvect, fVertex);
265  double lclusterpoint[5] = {std::abs(leadingjet->Pt()), leadingjet->NEF(), std::abs(pvect.Pt()), leadingjet->GetZ(pvect.Px(), pvect.Py(), pvect.Pz()), leadingjetvec.DeltaR(pvect.Vect())};
266  fHistos->FillTHnSparse(Form("hLeadingJetLeadingCluster%s", contname->String().Data()), lclusterpoint);
267  }
268  }
269  }
270  fHistos->FillTHnSparse(Form("hPtEtaPhiELeadingJet%s", contname->String().Data()), leadingvec);
271  }
272 
273  return kTRUE;
274 }
275 
277  using AnalysisHelpers = EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory;
278  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
279  if(!mgr) {
280  std::cerr << "[AliAnalysisTaskJetConstituentQA::AddTaskEmcalJetConstituentQA(EE)] No analysis manager provided ..." << std::endl;
281  return nullptr;
282  }
283 
284  std::string jettypestring;
285  switch(jettype){
286  case AliJetContainer::kFullJet: jettypestring = "fulljets"; break;
287  case AliJetContainer::kChargedJet: jettypestring = "chargedjets"; break;
288  case AliJetContainer::kNeutralJet: jettypestring = "neutraljets"; break;
289  };
290 
291  std::stringstream taskname;
292  taskname << "constituentQA_" << jettypestring << "_" << trigger;
293  auto task = new AliAnalysisTaskEmcalJetConstituentQA(taskname.str().data());
294  task->SetTriggerSelection(trigger);
295  task->SetJetType(jettype);
296  mgr->AddTask(task);
297 
298  auto inputhandler = mgr->GetInputEventHandler();
299  auto isAOD = false;
300  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = true;
301 
302  TString tracksname, clustername;
303  AliParticleContainer *tracks(nullptr);
304  AliClusterContainer *clusters(nullptr);
305  if(!partmode) {
306  if(jettype == AliJetContainer::kChargedJet || jettype == AliJetContainer::kFullJet){
307  tracksname = AnalysisHelpers::TrackContainerNameFactory(isAOD);
308  tracks = task->AddTrackContainer(tracksname);
309  task->SetNameTrackContainer(tracksname);
310  tracks->SetMinPt(0.15);
311  }
312 
313  if(jettype == AliJetContainer::kNeutralJet || jettype == AliJetContainer::kFullJet){
314  clustername = AnalysisHelpers::ClusterContainerNameFactory(isAOD);
315  clusters = task->AddClusterContainer(clustername);
316  task->SetNameClusterContainer(clustername);
317  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
318  clusters->SetClusHadCorrEnergyCut(0.3);
319  }
320  } else {
321  tracksname = "mcparticles";
322  tracks = task->AddParticleContainer(tracksname);
323  task->SetNameTrackContainer(tracksname);
324  tracks->SetMinPt(0.);
325  }
326 
327  // create jet containers for R02 and R04 jets
328  std::array<double, 4> jetradii = {{0.2, 0.3, 0.4, 0.5}};
329  for(auto r : jetradii) {
330  std::stringstream contname;
331  contname << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(r*10.);
332  auto jcont = task->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, r, AliJetContainer::kEMCALfid, tracks, clusters, "Jet");
333  jcont->SetName(contname.str().data());
334  task->AddNameJetContainer(contname.str().data());
335  jcont->SetMinPt(20.);
336  }
337 
338  std::stringstream contname, outfilename;
339  contname << "JetConstituentQA_" << jettypestring << "_" << trigger;
340  outfilename << mgr->GetCommonFileName() << ":JetConstituentQA_" << trigger;
341  if(partmode) {
342  contname << "_part";
343  outfilename << "_part";
344  }
345  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
346  mgr->ConnectOutput(task, 1, mgr->CreateContainer(contname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
347 
348  return task;
349 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
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.
double Double_t
Definition: External.C:58
AliVParticle * GetLeadingTrack(TClonesArray *tracks=0) const
Class creating a linear binning, used in the histogram manager.
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
TObjArray fNamesJetContainers
Names of the connected jet container.
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
TCanvas * c
Definition: TestFitELoss.C:172
static AliAnalysisTaskEmcalJetConstituentQA * AddTaskEmcalJetConstituentQA(const char *trigger, AliJetContainer::EJetType_t jettype, bool parmode=kFALSE)
Double_t E() const
Definition: AliEmcalJet.h:119
Interface for binnings used by the histogram handler.
Definition: TBinning.h:21
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Collection of helper functions used to configure the analysis.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
AliVCluster * GetLeadingCluster(TClonesArray *clusters=0) const
Double_t GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz) const
Double_t DeltaR(const AliVParticle *part) const
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalList * fOutput
!output list
Double_t fVertex[3]
!event vertex
AliTrackContainer * GetTrackContainer(Int_t i=0) const
TString fNameTriggerDecisionContainer
Name of the trigger decision container.
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
Double_t Pz() const
Definition: AliEmcalJet.h:108
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t NEF() const
Definition: AliEmcalJet.h:148
void SetDefaultClusterEnergy(Int_t d)
Bool_t fUseTriggerSelection
Use trigger selection in addition to trigger string.
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.
Container structure for EMCAL clusters.
Container for jet within the EMCAL jet framework.
void SetClusHadCorrEnergyCut(Double_t cut)