AliPhysics  6133e27 (6133e27)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetSubstructureTree.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 <iostream>
28 #include <string>
29 #include <sstream>
30 #include <vector>
31 
32 #include <fastjet/ClusterSequence.hh>
33 #include <fastjet/contrib/Nsubjettiness.hh>
34 #include <fastjet/contrib/SoftDrop.hh>
35 
36 #include <TMath.h>
37 #include <TString.h>
38 
39 #include "AliAODInputHandler.h"
40 #include "AliAnalysisManager.h"
42 #include "AliClusterContainer.h"
43 #include "AliJetContainer.h"
45 #include "AliEmcalJet.h"
46 #include "AliEmcalList.h"
47 #include "AliLog.h"
48 #include "AliParticleContainer.h"
49 #include "AliTrackContainer.h"
50 #include "AliVCluster.h"
51 #include "AliVParticle.h"
52 
56 
57 namespace EmcalTriggerJets {
58 
59 AliAnalysisTaskEmcalJetSubstructureTree::AliAnalysisTaskEmcalJetSubstructureTree() :
61  fJetSubstructureTree(nullptr),
62  fJetSubstructureInfo(),
63  fSDZCut(0.1),
64  fSDBetaCut(0),
65  fReclusterizer(kCAAlgo),
66  fTriggerSelectionBits(AliVEvent::kAny),
67  fTriggerSelectionString("")
68 {
69 }
70 
72  AliAnalysisTaskEmcalJet(name, kTRUE),
73  fJetSubstructureTree(nullptr),
74  fJetSubstructureInfo(),
75  fSDZCut(0.1),
76  fSDBetaCut(0),
77  fReclusterizer(kCAAlgo),
78  fTriggerSelectionBits(AliVEvent::kAny),
79  fTriggerSelectionString("")
80 {
81 
82 }
83 
85 
86 }
87 
90 
91  fJetSubstructureTree = new TTree("jetSubstructure", "Tree with jet substructure information");
92  std::stringstream leaflist;
93  leaflist << "fR/D:"
94  << "fEventWeight:"
95  << "fPtJetRec:"
96  << "fPtJetSim:"
97  << "fAreaRec:"
98  << "fAreaSim:"
99  << "fNEFRec:"
100  << "fNEFSim:"
101  << "fZgMeasured:"
102  << "fZgTrue:"
103  << "fRgMeasured:"
104  << "fRgTrue:"
105  << "fMgMeasured:"
106  << "fMgTrue:"
107  << "fPtgMeasured:"
108  << "fPtgTrue:"
109  << "fOneSubjettinessMeasured:"
110  << "fOneSubjettinessTrue:"
111  << "fTwoSubjettinessMeasured:"
112  << "fTwoSubjettinessTrue:"
113  << "fNCharged/I:"
114  << "fNNeutral:"
115  << "fNTrueConst:"
116  << "fNDroppedMeasured:"
117  << "fNDroppedTrue";
118  std::string leafstring = leaflist.str();
119  printf("branch string: %s\n", leafstring.c_str());
120  fJetSubstructureTree->Branch("JetInfo", &fJetSubstructureInfo, leafstring.c_str(), sizeof(fJetSubstructureInfo));
122  PostData(1, fOutput);
123 }
124 
126  AliClusterContainer *clusters = GetClusterContainer("caloClusters");
127  AliTrackContainer *tracks = GetTrackContainer("tracks");
128  AliParticleContainer *particles = GetParticleContainer("mcparticles");
129 
130  AliJetContainer *mcjets = GetJetContainer("mcjets");
131  AliJetContainer *datajets = GetJetContainer("datajets");
132 
133  // Run trigger selection (not on pure MCgen train)
134  if(datajets){
135  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
136  if(fTriggerSelectionString.Length()) {
137  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
138  }
139  }
140 
141  Double_t weight = 1;
142 
143  AliSoftdropDefinition softdropSettings;
144  softdropSettings.fBeta = fSDBetaCut;
145  softdropSettings.fZ = fSDZCut;
146  switch(fReclusterizer) {
147  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
148  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
149  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
150  };
151 
152  AliNSubjettinessDefinition nsubjettinessSettings;
153  nsubjettinessSettings.fBeta = 1.;
154  nsubjettinessSettings.fRadius = 0.4;
155 
156  if(mcjets && !datajets) {
157  // pure MC (gen) train - run over MC jets
158  for(auto jet : mcjets->accepted()) {
159  try {
160  AliJetSubstructureData structure = MakeJetSubstructure(*jet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
161  FillTree(mcjets->GetJetRadius(), weight, nullptr, jet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness));
162  } catch (ReclusterizerException &e) {
163  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
164  }
165  }
166  }
167 
168  if(datajets) {
169  for(auto jet : datajets->accepted()) {
170  AliEmcalJet *associatedJet = jet->MatchedJet();
171  if(mcjets) {
172  if(!associatedJet) continue;
173  try {
174  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, mcjets->GetJetRadius() * 2., particles, nullptr, {softdropSettings, nsubjettinessSettings}),
175  structureMC = MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings});
176  FillTree(datajets->GetJetRadius(), weight, jet, associatedJet, &(structureData.fSoftDrop), &(structureMC.fSoftDrop), &(structureData.fNsubjettiness), &(structureMC.fNsubjettiness));
177  } catch(ReclusterizerException &e) {
178  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
179  }
180  } else {
181  try {
182  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
183  FillTree(datajets->GetJetRadius(), weight, jet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness), nullptr);
184  } catch(ReclusterizerException &e) {
185  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
186  }
187  }
188  }
189  }
190 
191  return true;
192 }
193 
195  const AliEmcalJet *datajet, const AliEmcalJet *mcjet,
196  AliSoftDropParameters *dataSoftdrop, AliSoftDropParameters *mcSoftdrop,
197  AliNSubjettinessParameters *dataSubjettiness, AliNSubjettinessParameters *mcSubjettiness){
200  if(datajet) {
201  fJetSubstructureInfo.fPtJetRec = TMath::Abs(datajet->Pt());
204  fJetSubstructureInfo.fAreaRec = datajet->Area();
205  fJetSubstructureInfo.fNEFRec = datajet->NEF();
206  } else {
212  }
213 
214  if(mcjet) {
215  fJetSubstructureInfo.fPtJetSim = TMath::Abs(mcjet->Pt());
218  fJetSubstructureInfo.fNEFSim = mcjet->NEF();
219  } else {
224  }
225 
226  if(dataSoftdrop) {
227  fJetSubstructureInfo.fZgMeasured = dataSoftdrop->fZg;
228  fJetSubstructureInfo.fRgMeasured = dataSoftdrop->fRg;
229  fJetSubstructureInfo.fMgMeasured = dataSoftdrop->fMg;
230  fJetSubstructureInfo.fPtgMeasured = dataSoftdrop->fPtg;
232  } else {
238  }
239 
240  if(mcSoftdrop) {
241  fJetSubstructureInfo.fZgTrue = mcSoftdrop->fZg;
242  fJetSubstructureInfo.fRgTrue = mcSoftdrop->fRg;
243  fJetSubstructureInfo.fMgTrue = mcSoftdrop->fMg;
244  fJetSubstructureInfo.fPtgTrue = mcSoftdrop->fPtg;
246  } else {
252  }
253 
254  if(dataSubjettiness) {
257  } else {
260  }
261 
262  if(mcSubjettiness) {
265  } else {
268  }
269 
270  fJetSubstructureTree->Fill();
271 }
272 
273 
275  const int kClusterOffset = 30000; // In order to handle tracks and clusters in the same index space the cluster index needs and offset, large enough so that there is no overlap with track indices
276  std::vector<fastjet::PseudoJet> constituents;
277  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
278  AliVTrack *track = static_cast<AliVTrack *>(jet.TrackAt(itrk, tracks->GetArray()));
279  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
280  constituentTrack.set_user_index(jet.TrackAt(itrk));
281  constituents.push_back(constituentTrack);
282  }
283 
284  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
285  AliVCluster *cluster = jet.ClusterAt(icl, clusters->GetArray());
286  TLorentzVector clustervec;
287  cluster->GetMomentum(clustervec, fVertex);
288  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
289  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
290  constituents.push_back(constituentCluster);
291  }
292 
293  // Redo jet finding on constituents with a
294  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
295  std::vector<fastjet::PseudoJet> outputjets;
296  try {
297  fastjet::ClusterSequence jetfinder(constituents, jetdef);
298  outputjets = jetfinder.inclusive_jets(0);
299  AliJetSubstructureData result({MakeSoftDropParameters(outputjets[0], settings.fSoftdropSettings), MakeNsubjettinessParameters(outputjets[0], settings.fSubjettinessSettings)});
300  return result;
301  } catch (fastjet::Error &e) {
302  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
303  throw ReclusterizerException();
304  }
305 }
306 
308  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
309  softdropAlgorithm.set_verbose_structure(kTRUE);
310  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
311  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
312  fastjet::PseudoJet groomed = softdropAlgorithm(jet);
313 
314  AliSoftDropParameters result({groomed.structure_of<fastjet::contrib::SoftDrop>().symmetry(),
315  groomed.structure_of<fastjet::contrib::SoftDrop>().delta_R(),
316  groomed.structure_of<fastjet::contrib::SoftDrop>().mu(),
317  groomed.perp(),
318  groomed.structure_of<fastjet::contrib::SoftDrop>().dropped_count()});
319  return result;
320 }
321 
324  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(cut.fBeta,cut.fRadius)).result(jet),
325  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(cut.fBeta,cut.fRadius)).result(jet)
326  });
327  return result;
328 }
329 
331  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
332 
333  Bool_t isAOD(kFALSE);
334  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
335  if(inputhandler) {
336  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = kTRUE;
337  }
338 
339  AliAnalysisTaskEmcalJetSubstructureTree *treemaker = new AliAnalysisTaskEmcalJetSubstructureTree("JetSubstructureTreemaker_" + TString::Format("R%02d_", int(jetradius * 10.)) + trigger);
340  mgr->AddTask(treemaker);
341  treemaker->SetMakeGeneralHistograms(kTRUE);
342  treemaker->SetVzRange(-10., 10);
343 
344  // Adding containers
345  if(isMC) {
346  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
347 
348  AliJetContainer *mcjets = treemaker->AddJetContainer(
352  jetradius,
354  particles, nullptr);
355  mcjets->SetName("mcjets");
356  mcjets->SetJetPtCut(20.);
357  }
358 
359  if(isData) {
361  tracks->SetMinPt(0.15);
363  clusters->SetMinE(0.3); // 300 MeV E-cut
364 
365  AliJetContainer *datajets = treemaker->AddJetContainer(
369  jetradius,
371  tracks, clusters);
372  datajets->SetName("datajets");
373  datajets->SetJetPtCut(20.);
374 
375  treemaker->SetUseAliAnaUtils(true, true);
376 
377  // configure trigger selection
378  TString triggerstring(trigger);
379  if(triggerstring.Contains("INT7")) {
380  treemaker->SetTriggerBits(AliVEvent::kINT7);
381  } else if(triggerstring.Contains("EJ1")) {
382  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
383  treemaker->SetTriggerString("EJ1");
384  } else if(triggerstring.Contains("EJ2")) {
385  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
386  treemaker->SetTriggerString("EJ2");
387  }
388  }
389 
390  // Connecting containers
391  TString outputfile = mgr->GetCommonFileName();
392  outputfile += TString::Format(":JetSubstructure_R%02d_%s", int(jetradius * 10.), trigger);
393  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
394  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer("JetSubstructure_" + TString::Format("R%0d_", int(jetradius * 10.)) + trigger, AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile));
395 
396  return treemaker;
397 }
398 
399 } /* namespace EmcalTriggerJets */
Double_t Area() const
Definition: AliEmcalJet.h:123
double Double_t
Definition: External.C:58
AliJetContainer * GetJetContainer(Int_t i=0) const
Container with name, TClonesArray and cuts for particles.
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
AliJetSubstructureData MakeJetSubstructure(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters, const AliJetSubstructureSettings &settings) const
AliEmcalJet * MatchedJet() const
Definition: AliEmcalJet.h:226
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:130
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Structure for results from the soft drop algorithm.
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.
void SetVzRange(Double_t min, Double_t max)
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
fastjet::JetAlgorithm fRecluserAlgo
Reclusterization algorithm.
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:131
void SetJetPtCut(Float_t cut)
void FillTree(double r, double weight, const AliEmcalJet *datajet, const AliEmcalJet *mcjet, AliSoftDropParameters *dataSoftdrop, AliSoftDropParameters *mcsoftdrop, AliNSubjettinessParameters *dataSubjettiness, AliNSubjettinessParameters *mcSubjettiness)
AliNSubjettinessParameters MakeNsubjettinessParameters(const fastjet::PseudoJet &jet, const AliNSubjettinessDefinition &cut) const
Definition for the algorithm obtaining the softdrop parameters.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
Double_t Pt() const
Definition: AliEmcalJet.h:102
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
Bool_t isMC
Float_t GetJetRadius() const
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
AliJetSubstructureInfo fJetSubstructureInfo
! Jet Substructure information to be filled in the tree
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
static AliAnalysisTaskEmcalJetSubstructureTree * AddEmcalJetSubstructureTreeMaker(Bool_t isMC, Bool_t isData, Double_t jetradius, const char *name)
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
AliSoftDropParameters MakeSoftDropParameters(const fastjet::PseudoJet &jet, const AliSoftdropDefinition &cut) const
void UserCreateOutputObjects()
Main initialization function on the worker.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Double_t NEF() const
Definition: AliEmcalJet.h:141
Container structure for EMCAL clusters.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:63
Container for jet within the EMCAL jet framework.
static TString ClusterContainerNameFactory(Bool_t isAOD)
Get name of the default cluster container.
static TString TrackContainerNameFactory(Bool_t isAOD)
Get name of the default track container.