AliPhysics  aaf9c62 (aaf9c62)
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  this->SetMakeGeneralHistograms(true);
85 }
86 
88  if(fHistos) delete fHistos;
89 }
90 
93 
94  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.),
95  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());
96 
97  const TBinning *jetbinning[4] = {&binningptjet, &binningnef, &multbinning, &multbinning},
98  *chargedbinning[7] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningz, &binningR},
99  *neutralbinning[9] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningptconst, &binningz, &binningR, &binningptconst},
100  *binningHighZClusters[7] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningNCell, &binningFracCellLeading, &binningM02},
101  *leadingchargedbinning[5] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningR},
102  *leadingneutralbinning[6] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningR, &binningptconst},
103  *leadingjetvecbinning[4] = {&binningptjet, &etabinning, &phibinning, &binningptjet};
104 
105  fHistos = new THistManager(Form("histos_%s", GetName()));
106  for(auto c : fNamesJetContainers){
107  auto contname = dynamic_cast<TObjString *>(c);
108  if(!contname) continue;
109  fHistos->CreateTHnSparse(Form("hJetCounter%s", contname->String().Data()), Form("jet counter for jets %s", contname->String().Data()), 4, jetbinning);
110  fHistos->CreateTHnSparse(Form("hPtEtaPhiELeadingJet%s", contname->String().Data()), Form("Momemtum vector of leading jets %s", contname->String().Data()), 4, leadingjetvecbinning);
112  fHistos->CreateTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), Form("charged constituents in jets %s", contname->String().Data()), 7, chargedbinning);
113  fHistos->CreateTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), Form("leading charged constituent in jets %s", contname->String().Data()), 5, leadingchargedbinning);
114  fHistos->CreateTHnSparse(Form("hLeadingJetLeadingTrack%s", contname->String().Data()), Form("leading charged constituent in jets %s", contname->String().Data()), 5, leadingchargedbinning);
115  }
117  fHistos->CreateTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), Form("neutral constituents in jets %s", contname->String().Data()), 9, neutralbinning);
118  fHistos->CreateTHnSparse(Form("hHighZClusters%s", contname->String().Data()), "Properties of high-z clusters", 7, binningHighZClusters);
119  fHistos->CreateTHnSparse(Form("hLeadingCluster%s", contname->String().Data()), Form("leading neutral constituent in jets %s", contname->String().Data()), 6, leadingneutralbinning);
120  fHistos->CreateTHnSparse(Form("hLeadingJetLeadingCluster%s", contname->String().Data()), Form("leading neutral constituent in jets %s", contname->String().Data()), 6, leadingneutralbinning);
121  }
122  }
123 
124  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
125  PostData(1, fOutput);
126 }
127 
130  if(!tracks) tracks = GetParticleContainer(fNameTrackContainer);
131  const auto clusters = GetClusterContainer(fNameClusterContainer);
132  if(fNameTrackContainer.Length() && !tracks){
133  AliErrorStream() << "Track container " << fNameTrackContainer << " required but missing ..." << std::endl;
134  return kFALSE;
135  }
136  if(fNameClusterContainer.Length() && !clusters){
137  AliErrorStream() << "Cluster container " << fNameClusterContainer << " required but missing ..." << std::endl;
138  return kFALSE;
139  }
140 
141 
142  for(auto jc : fNamesJetContainers){
143  auto contname = dynamic_cast<TObjString *>(jc);
144  if(!contname) {
145  AliErrorStream() << "Non-string object in the list of jet container names" << std::endl;
146  continue;
147  }
148  const auto jetcont = GetJetContainer(contname->String().Data());
149  if(!jetcont){
150  AliErrorStream() << "Jet container with name " << contname->String() << " not found in the list of jet containers" << std::endl;
151  continue;
152  }
153  AliDebugStream(2) << "Reading " << jetcont->GetArray()->GetName() << std::endl;
154 
155  AliEmcalJet *leadingjet(nullptr);
156  for(auto jet : jetcont->accepted()){
157  if(!leadingjet || jet->Pt() > leadingjet->Pt()) leadingjet = jet;
158  AliDebugStream(3) << "Next accepted jet, found " << jet->GetNumberOfTracks() << " tracks and " << jet->GetNumberOfClusters() << " clusters." << std::endl;
159  Double_t pointjet[4] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters())},
160  pointcharged[7] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1.},
161  pointneutral[9] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1., -1.},
162  pointHighZCluster[7] = {std::abs(jet->Pt()), jet->NEF(), -1., -1., -1., -1., -1.};
163  fHistos->FillTHnSparse(Form("hJetCounter%s", contname->String().Data()), pointjet);
164  TVector3 jetvec{jet->Px(), jet->Py(), jet->Pz()};
166  for(decltype(jet->GetNumberOfTracks()) itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
167  const auto trk = jet->TrackAt(itrk, tracks->GetArray());
168  if(!trk) continue;
169  if(trk->Charge()){
170  pointcharged[4] = std::abs(trk->Pt());
171  pointcharged[5] = std::abs(jet->GetZ(trk));
172  pointcharged[6] = jet->DeltaR(trk);
173  fHistos->FillTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), pointcharged);
174  } else {
175  // particle level jets
176  pointneutral[4] = pointneutral[5] = std::abs(trk->E());
177  pointneutral[6] = std::abs(jet->GetZ(trk));
178  pointneutral[7] = jet->DeltaR(trk);
179  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
180  }
181  }
182  // Leading track
183  auto leadingtrack = jet->GetLeadingTrack(tracks->GetArray());
184  if(leadingtrack) {
185  double ltrackpoint[5] = {std::abs(jet->Pt()), jet->NEF(), std::abs(leadingtrack->Pt()), jet->GetZ(leadingtrack), jet->DeltaR(leadingtrack)};
186  fHistos->FillTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), ltrackpoint);
187  }
188  }
190  for(decltype(jet->GetNumberOfClusters()) icl = 0; icl < jet->GetNumberOfClusters(); icl++){
191  const auto clust = jet->ClusterAt(icl, clusters->GetArray());
192  std::vector<double> fracamp(clust->GetNCells());
193  memcpy(fracamp.data(), clust->GetCellsAmplitudeFraction(), sizeof(double) * clust->GetNCells());
194  if(!clust) continue;
195  TLorentzVector ptvec;
196  clust->GetMomentum(ptvec, this->fVertex, AliVCluster::kHadCorr);
197  pointneutral[4] = std::abs(clust->GetHadCorrEnergy());
198  pointneutral[5] = std::abs(clust->GetNonLinCorrEnergy());
199  pointneutral[6] = jet->GetZ(ptvec.Px(), ptvec.Py(), ptvec.Pz());
200  pointneutral[7] = jetvec.DeltaR(ptvec.Vect());
201  pointneutral[8] = std::abs(clust->E() * (*std::max_element(fracamp.begin(), fracamp.end())));
202  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
203 
204  if(pointneutral[6] > 0.95) {
205  pointHighZCluster[2] = pointneutral[4];
206  pointHighZCluster[3] = pointneutral[6];
207  pointHighZCluster[4] = clust->GetNCells();
208  pointHighZCluster[5] = *std::max_element(clust->GetCellsAmplitudeFraction(), clust->GetCellsAmplitudeFraction()+clust->GetNCells());
209  pointHighZCluster[6] = clust->GetM02();
210  fHistos->FillTHnSparse(Form("hHighZClusters%s", contname->String().Data()), pointHighZCluster);
211  }
212  }
213  // Leading cluster
214  auto leadingcluster = jet->GetLeadingCluster(clusters->GetArray());
215  if(leadingcluster){
216  TLorentzVector pvect;
217  leadingcluster->GetMomentum(pvect, fVertex);
218  std::vector<double> fracamp(leadingcluster->GetNCells());
219  memcpy(fracamp.data(), leadingcluster->GetCellsAmplitudeFraction(), sizeof(double) * leadingcluster->GetNCells());
220  double lclusterpoint[6] = {std::abs(jet->Pt()), jet->NEF(), std::abs(pvect.Pt()), jet->GetZ(pvect.Px(), pvect.Py(), pvect.Pz()), jetvec.DeltaR(pvect.Vect()), std::abs(leadingcluster->E() * (*std::max_element(fracamp.begin(), fracamp.end())))};
221  fHistos->FillTHnSparse(Form("hLeadingCluster%s", contname->String().Data()), lclusterpoint);
222  }
223  }
224  }
225 
226  // Look at leading particles in the leading jet
227  double leadingvec[4] = {0., 0., 0., 0.};
228  if(leadingjet) {
229  TVector3 leadingjetvec{leadingjet->Px(), leadingjet->Py(), leadingjet->Pz()};
230  leadingvec[0] = std::abs(leadingjet->Pt());
231  leadingvec[1] = leadingjet->Eta();
232  leadingvec[2] = leadingjet->Phi();
233  if(leadingvec[2] < 0) leadingvec[2] += TMath::TwoPi();
234  leadingvec[3] = leadingjet->E();
236  auto leadingtrack = leadingjet->GetLeadingTrack(tracks->GetArray());
237  if(leadingtrack) {
238  double ltrackpoint[5] = {std::abs(leadingjet->Pt()), leadingjet->NEF(), std::abs(leadingtrack->Pt()), leadingjet->GetZ(leadingtrack), leadingjet->DeltaR(leadingtrack)};
239  fHistos->FillTHnSparse(Form("hLeadingTrack%s", contname->String().Data()), ltrackpoint);
240  }
241  }
243  auto leadingcluster = leadingjet->GetLeadingCluster(clusters->GetArray());
244  if(leadingcluster){
245  TLorentzVector pvect;
246  leadingcluster->GetMomentum(pvect, fVertex);
247  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())};
248  fHistos->FillTHnSparse(Form("hLeadingJetLeadingCluster%s", contname->String().Data()), lclusterpoint);
249  }
250  }
251  }
252  fHistos->FillTHnSparse(Form("hPtEtaPhiELeadingJet%s", contname->String().Data()), leadingvec);
253  }
254 
255  return kTRUE;
256 }
257 
259  // Event selection
260  AliDebugStream(1) << "Trigger selection string: " << fTriggerSelectionString << ", fired trigger classes: " << fInputEvent->GetFiredTriggerClasses() << std::endl;
261  if(fTriggerSelectionString.Contains("INT7")){
262  // INT7 trigger
263  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
264  } else if(fTriggerSelectionString.Contains("EJ")){
265  auto triggerclass = fTriggerSelectionString(fTriggerSelectionString.Index("EJ"),3); // cleanup trigger string from possible tags
266  AliDebugStream(1) << "Inspecting trigger class " << triggerclass << std::endl;
267  // EMCAL JET trigger
268  if(!fMCEvent){ // Running on Data
269  if(!(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)) return false;
270  if(!fInputEvent->GetFiredTriggerClasses().Contains(triggerclass)) return false;
271  }
272 
274  auto triggerdecisions = dynamic_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer.Data()));
275  if(!triggerdecisions) {
276  AliErrorStream() << "No offline trigger selection available" << std::endl;
277  return false;
278  }
279  else if(!triggerdecisions->IsEventSelected(triggerclass.Data())) return false;
280  }
281  } else return false;
282 
283  AliDebugStream(1) << "Event is selected" << std::endl;
284  return true;
285 }
286 
288  using AnalysisHelpers = EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory;
289  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
290  if(!mgr) {
291  std::cerr << "[AliAnalysisTaskJetConstituentQA::AddTaskEmcalJetConstituentQA(EE)] No analysis manager provided ..." << std::endl;
292  return nullptr;
293  }
294 
295  std::string jettypestring;
296  switch(jettype){
297  case AliJetContainer::kFullJet: jettypestring = "fulljets"; break;
298  case AliJetContainer::kChargedJet: jettypestring = "chargedjets"; break;
299  case AliJetContainer::kNeutralJet: jettypestring = "neutraljets"; break;
300  };
301 
302  std::stringstream taskname;
303  taskname << "constituentQA_" << jettypestring << "_" << trigger;
304  auto task = new AliAnalysisTaskEmcalJetConstituentQA(taskname.str().data());
305  task->SetTriggerSelection(trigger);
306  task->SetJetType(jettype);
307  mgr->AddTask(task);
308 
309  auto inputhandler = mgr->GetInputEventHandler();
310  auto isAOD = false;
311  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = true;
312 
313  TString tracksname, clustername;
314  AliParticleContainer *tracks(nullptr);
315  AliClusterContainer *clusters(nullptr);
316  if(!partmode) {
317  if(jettype == AliJetContainer::kChargedJet || jettype == AliJetContainer::kFullJet){
318  tracksname = AnalysisHelpers::TrackContainerNameFactory(isAOD);
319  tracks = task->AddTrackContainer(tracksname);
320  task->SetNameTrackContainer(tracksname);
321  tracks->SetMinPt(0.15);
322  }
323 
324  if(jettype == AliJetContainer::kNeutralJet || jettype == AliJetContainer::kFullJet){
325  clustername = AnalysisHelpers::ClusterContainerNameFactory(isAOD);
326  clusters = task->AddClusterContainer(clustername);
327  task->SetNameClusterContainer(clustername);
328  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
329  clusters->SetClusHadCorrEnergyCut(0.3);
330  }
331  } else {
332  tracksname = "mcparticles";
333  tracks = task->AddParticleContainer(tracksname);
334  task->SetNameTrackContainer(tracksname);
335  tracks->SetMinPt(0.);
336  }
337 
338  // create jet containers for R02 and R04 jets
339  std::array<double, 4> jetradii = {{0.2, 0.3, 0.4, 0.5}};
340  for(auto r : jetradii) {
341  std::stringstream contname;
342  contname << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(r*10.);
343  auto jcont = task->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, r, AliJetContainer::kEMCALfid, tracks, clusters, "Jet");
344  jcont->SetName(contname.str().data());
345  task->AddNameJetContainer(contname.str().data());
346  jcont->SetMinPt(20.);
347  }
348 
349  std::stringstream contname, outfilename;
350  contname << "JetConstituentQA_" << jettypestring << "_" << trigger;
351  outfilename << mgr->GetCommonFileName() << ":JetConstituentQA_" << trigger;
352  if(partmode) {
353  contname << "_part";
354  outfilename << "_part";
355  }
356  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
357  mgr->ConnectOutput(task, 1, mgr->CreateContainer(contname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
358 
359  return task;
360 }
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
void SetMakeGeneralHistograms(Bool_t g)
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.
bool Bool_t
Definition: External.C:53
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)