AliPhysics  863cc41 (863cc41)
 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  fSDZCut(0.1),
63  fSDBetaCut(0),
64  fReclusterizer(kCAAlgo),
65  fTriggerSelectionBits(AliVEvent::kAny),
66  fTriggerSelectionString("")
67 {
68  memset(fJetTreeData, 0, sizeof(Double_t) * kTNVar);
69 }
70 
72  AliAnalysisTaskEmcalJet(name, kTRUE),
73  fJetSubstructureTree(nullptr),
74  fSDZCut(0.1),
75  fSDBetaCut(0),
76  fReclusterizer(kCAAlgo),
77  fTriggerSelectionBits(AliVEvent::kAny),
78  fTriggerSelectionString("")
79 {
80  memset(fJetTreeData, 0, sizeof(Double_t) * kTNVar);
81  DefineOutput(2, TTree::Class());
82 }
83 
85 
86 }
87 
90 
91  OpenFile(2);
92  fJetSubstructureTree = new TTree("jetSubstructure", "Tree with jet substructure information");
93  TString varnames[kTNVar];
94  varnames[0] = "Radius";
95  varnames[1] = "EventWeight";
96  varnames[2] = "PtJetRec";
97  varnames[3] = "PtJetSim";
98  varnames[4] = "AreaRec";
99  varnames[5] = "AreaSim";
100  varnames[6] = "NEFRec";
101  varnames[7] = "NEFSim";
102  varnames[8] = "MassRec";
103  varnames[9] = "MassSim";
104  varnames[10] = "ZgMeasured";
105  varnames[11] = "ZgTrue";
106  varnames[12] = "RgMeasured";
107  varnames[13] = "RgTrue";
108  varnames[14] = "MgMeasured";
109  varnames[15] = "MgTrue";
110  varnames[16] = "PtgMeasured";
111  varnames[17] = "PtgTrue";
112  varnames[18] = "OneSubjettinessMeasured";
113  varnames[19] = "OneSubjettinessTrue";
114  varnames[20] = "TwoSubjettinessMeasured";
115  varnames[21] = "TwoSubjettinessTrue";
116  varnames[22] = "NCharged";
117  varnames[23] = "NNeutral";
118  varnames[24] = "NConstTrue";
119  varnames[25] = "NDroppedMeasured";
120  varnames[26] = "NDroppedTrue";
121 
122  for(int ib = 0; ib < kTNVar; ib++){
123  fJetSubstructureTree->Branch(varnames[ib], fJetTreeData + ib, Form("%s/D", varnames[ib].Data()));
124  }
125  PostData(1, fOutput);
126  PostData(2, fJetSubstructureTree);
127 }
128 
130  AliClusterContainer *clusters = GetClusterContainer("caloClusters");
131  AliTrackContainer *tracks = GetTrackContainer("tracks");
132  AliParticleContainer *particles = GetParticleContainer("mcparticles");
133 
134  AliJetContainer *mcjets = GetJetContainer("mcjets");
135  AliJetContainer *datajets = GetJetContainer("datajets");
136 
137  // Run trigger selection (not on pure MCgen train)
138  if(datajets){
139  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
140  if(fTriggerSelectionString.Length()) {
141  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
142  }
143  }
144 
145  Double_t weight = 1;
146 
147  AliSoftdropDefinition softdropSettings;
148  softdropSettings.fBeta = fSDBetaCut;
149  softdropSettings.fZ = fSDZCut;
150  switch(fReclusterizer) {
151  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
152  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
153  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
154  };
155 
156  AliNSubjettinessDefinition nsubjettinessSettings;
157  nsubjettinessSettings.fBeta = 1.;
158  nsubjettinessSettings.fRadius = 0.4;
159 
160  if(mcjets && !datajets) {
161  // pure MC (gen) train - run over MC jets
162  for(auto jet : mcjets->accepted()) {
163  try {
164  AliJetSubstructureData structure = MakeJetSubstructure(*jet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
165  FillTree(mcjets->GetJetRadius(), weight, nullptr, jet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness));
166  } catch (ReclusterizerException &e) {
167  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
168  }
169  }
170  }
171 
172  if(datajets) {
173  for(auto jet : datajets->accepted()) {
174  AliEmcalJet *associatedJet = jet->MatchedJet();
175  if(mcjets) {
176  if(!associatedJet) continue;
177  try {
178  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, mcjets->GetJetRadius() * 2., particles, nullptr, {softdropSettings, nsubjettinessSettings}),
179  structureMC = MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings});
180  FillTree(datajets->GetJetRadius(), weight, jet, associatedJet, &(structureData.fSoftDrop), &(structureMC.fSoftDrop), &(structureData.fNsubjettiness), &(structureMC.fNsubjettiness));
181  } catch(ReclusterizerException &e) {
182  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
183  }
184  } else {
185  try {
186  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
187  FillTree(datajets->GetJetRadius(), weight, jet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness), nullptr);
188  } catch(ReclusterizerException &e) {
189  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
190  }
191  }
192  }
193  }
194 
195  return true;
196 }
197 
199  const AliEmcalJet *datajet, const AliEmcalJet *mcjet,
200  AliSoftDropParameters *dataSoftdrop, AliSoftDropParameters *mcSoftdrop,
201  AliNSubjettinessParameters *dataSubjettiness, AliNSubjettinessParameters *mcSubjettiness){
202  fJetTreeData[kTRadius] = r;
203  fJetTreeData[kTWeight] = weight;
204  if(datajet) {
205  fJetTreeData[kTPtJetRec] = TMath::Abs(datajet->Pt());
207  fJetTreeData[kTNNeutral] = datajet->GetNumberOfClusters();
208  fJetTreeData[kTAreaRec] = datajet->Area();
209  fJetTreeData[kTNEFRec] = datajet->NEF();
210  fJetTreeData[kTMassRec] = datajet->M();
211  } else {
212  fJetTreeData[kTPtJetRec] = 0.;
215  fJetTreeData[kTAreaRec] = 0.;
216  fJetTreeData[kTNEFRec] = 0.;
217  fJetTreeData[kTMassRec] = 0.;
218  }
219 
220  if(mcjet) {
221  fJetTreeData[kTPtJetSim] = TMath::Abs(mcjet->Pt());
223  fJetTreeData[kTAreaSim] = mcjet->Area();
224  fJetTreeData[kTNEFSim] = mcjet->NEF();
225  fJetTreeData[kTMassSim] = mcjet->M();
226  } else {
227  fJetTreeData[kTPtJetSim] = 0.;
229  fJetTreeData[kTAreaSim] = 0.;
230  fJetTreeData[kTNEFSim] = 0.;
231  fJetTreeData[kTMassSim] = 0;
232  }
233 
234  if(dataSoftdrop) {
235  fJetTreeData[kTZgMeasured] = dataSoftdrop->fZg;
236  fJetTreeData[kTRgMeasured] = dataSoftdrop->fRg;
237  fJetTreeData[kTMgMeasured] = dataSoftdrop->fMg;
238  fJetTreeData[kTPtgMeasured] = dataSoftdrop->fPtg;
239  fJetTreeData[kTNDroppedMeasured] = dataSoftdrop->fNDropped;
240  } else {
246  }
247 
248  if(mcSoftdrop) {
249  fJetTreeData[kTZgTrue] = mcSoftdrop->fZg;
250  fJetTreeData[kTRgTrue] = mcSoftdrop->fRg;
251  fJetTreeData[kTMgTrue] = mcSoftdrop->fMg;
252  fJetTreeData[kTPtgTrue] = mcSoftdrop->fPtg;
253  fJetTreeData[kTNDroppedTrue] = mcSoftdrop->fNDropped;
254  } else {
255  fJetTreeData[kTZgTrue] = 0.;
256  fJetTreeData[kTRgTrue] = 0.;
257  fJetTreeData[kTMgTrue] = 0.;
258  fJetTreeData[kTPtgTrue] = 0.;
260  }
261 
262  if(dataSubjettiness) {
265  } else {
268  }
269 
270  if(mcSubjettiness) {
273  } else {
276  }
277 
278  fJetSubstructureTree->Fill();
279 }
280 
281 
283  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
284  std::vector<fastjet::PseudoJet> constituents;
285  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
286  AliVTrack *track = static_cast<AliVTrack *>(jet.TrackAt(itrk, tracks->GetArray()));
287  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
288  constituentTrack.set_user_index(jet.TrackAt(itrk));
289  constituents.push_back(constituentTrack);
290  }
291 
292  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
293  AliVCluster *cluster = jet.ClusterAt(icl, clusters->GetArray());
294  TLorentzVector clustervec;
295  cluster->GetMomentum(clustervec, fVertex);
296  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
297  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
298  constituents.push_back(constituentCluster);
299  }
300 
301  // Redo jet finding on constituents with a
302  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
303  std::vector<fastjet::PseudoJet> outputjets;
304  try {
305  fastjet::ClusterSequence jetfinder(constituents, jetdef);
306  outputjets = jetfinder.inclusive_jets(0);
307  AliJetSubstructureData result({MakeSoftDropParameters(outputjets[0], settings.fSoftdropSettings), MakeNsubjettinessParameters(outputjets[0], settings.fSubjettinessSettings)});
308  return result;
309  } catch (fastjet::Error &e) {
310  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
311  throw ReclusterizerException();
312  }
313 }
314 
316  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
317  softdropAlgorithm.set_verbose_structure(kTRUE);
318  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
319  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
320  fastjet::PseudoJet groomed = softdropAlgorithm(jet);
321 
322  AliSoftDropParameters result({groomed.structure_of<fastjet::contrib::SoftDrop>().symmetry(),
323  groomed.structure_of<fastjet::contrib::SoftDrop>().delta_R(),
324  groomed.structure_of<fastjet::contrib::SoftDrop>().mu(),
325  groomed.perp(),
326  groomed.structure_of<fastjet::contrib::SoftDrop>().dropped_count()});
327  return result;
328 }
329 
332  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(cut.fBeta,cut.fRadius)).result(jet),
333  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(cut.fBeta,cut.fRadius)).result(jet)
334  });
335  return result;
336 }
337 
339  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
340 
341  Bool_t isAOD(kFALSE);
342  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
343  if(inputhandler) {
344  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = kTRUE;
345  }
346 
347  AliAnalysisTaskEmcalJetSubstructureTree *treemaker = new AliAnalysisTaskEmcalJetSubstructureTree("JetSubstructureTreemaker_" + TString::Format("R%02d_", int(jetradius * 10.)) + trigger);
348  mgr->AddTask(treemaker);
349  treemaker->SetMakeGeneralHistograms(kTRUE);
350  treemaker->SetVzRange(-10., 10);
351 
352  // Adding containers
353  if(isMC) {
354  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
355 
356  AliJetContainer *mcjets = treemaker->AddJetContainer(
360  jetradius,
362  particles, nullptr);
363  mcjets->SetName("mcjets");
364  mcjets->SetJetPtCut(20.);
365  }
366 
367  if(isData) {
369  tracks->SetMinPt(0.15);
371  clusters->SetMinE(0.3); // 300 MeV E-cut
372 
373  AliJetContainer *datajets = treemaker->AddJetContainer(
377  jetradius,
379  tracks, clusters);
380  datajets->SetName("datajets");
381  datajets->SetJetPtCut(20.);
382 
383  treemaker->SetUseAliAnaUtils(true, true);
384 
385  // configure trigger selection
386  TString triggerstring(trigger);
387  if(triggerstring.Contains("INT7")) {
388  treemaker->SetTriggerBits(AliVEvent::kINT7);
389  } else if(triggerstring.Contains("EJ1")) {
390  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
391  treemaker->SetTriggerString("EJ1");
392  } else if(triggerstring.Contains("EJ2")) {
393  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
394  treemaker->SetTriggerString("EJ2");
395  }
396  }
397 
398  // Connecting containers
399  TString outputfile = mgr->GetCommonFileName();
400  outputfile += TString::Format(":JetSubstructure_R%02d_%s", int(jetradius * 10.), trigger);
401  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
402  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer("JetSubstructureHistos_" + TString::Format("R%0d_", int(jetradius * 10.)) + trigger, AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile));
403  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer("JetSubstuctureTree_" + TString::Format("R%0d_", int(jetradius * 10.)) + trigger, TTree::Class(), AliAnalysisManager::kOutputContainer, Form("JetSubstructureTree_%s.root", trigger)));
404 
405  return treemaker;
406 }
407 
408 } /* 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.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
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
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.
Double_t M() const
Definition: AliEmcalJet.h:113
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:63
Container for jet within the EMCAL jet framework.
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
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.