AliPhysics  bdbde52 (bdbde52)
 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 <set>
30 #include <sstream>
31 #include <vector>
32 
33 #include <fastjet/ClusterSequence.hh>
34 #include <fastjet/contrib/Nsubjettiness.hh>
35 #include <fastjet/contrib/SoftDrop.hh>
36 
37 #include <THistManager.h>
38 #include <TLinearBinning.h>
39 #include <TLorentzVector.h>
40 #include <TMath.h>
41 #include <TString.h>
42 #include <TVector3.h>
43 
44 #include "AliAODInputHandler.h"
45 #include "AliAnalysisManager.h"
47 #include "AliClusterContainer.h"
48 #include "AliJetContainer.h"
50 #include "AliEmcalJet.h"
51 #include "AliEmcalList.h"
53 #include "AliLog.h"
54 #include "AliParticleContainer.h"
55 #include "AliTrackContainer.h"
56 #include "AliRhoParameter.h"
57 #include "AliVCluster.h"
58 #include "AliVParticle.h"
59 
60 #ifdef EXPERIMENTAL_JETCONSTITUENTS
61 #include "AliEmcalClusterJetConstituent.h"
62 #include "AliEmcalParticleJetConstituent.h"
63 #endif
64 
65 
69 
70 namespace EmcalTriggerJets {
71 
74  fJetSubstructureTree(nullptr),
75  fQAHistos(nullptr),
76  fSDZCut(0.1),
77  fSDBetaCut(0),
78  fReclusterizer(kCAAlgo),
79  fTriggerSelectionBits(AliVEvent::kAny),
80  fTriggerSelectionString("")
81 {
82  memset(fJetTreeData, 0, sizeof(Double_t) * kTNVar);
83 }
84 
86  AliAnalysisTaskEmcalJet(name, kTRUE),
87  fJetSubstructureTree(nullptr),
88  fQAHistos(nullptr),
89  fSDZCut(0.1),
90  fSDBetaCut(0),
91  fReclusterizer(kCAAlgo),
92  fTriggerSelectionBits(AliVEvent::kAny),
93  fTriggerSelectionString("")
94 {
95  memset(fJetTreeData, 0, sizeof(Double_t) * kTNVar);
96  DefineOutput(2, TTree::Class());
97 }
98 
100 
101 }
102 
105 
106  // Make QA for constituent clusters
107  TLinearBinning jetptbinning(9, 20, 200),
108  clusterenergybinning(200, 0., 200),
109  timebinning(1000, -500., 500.),
110  m02binning(100, 0., 1.),
111  ncellbinning(101, -0.5, 100.5);
112  fQAHistos = new THistManager("QAhistos");
113  fQAHistos->CreateTH2("hClusterConstE", "EMCAL cluster energy vs jet pt; p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
114  fQAHistos->CreateTH2("hClusterConstTime", "EMCAL cluster time vs. jet pt; p_{t, jet} (GeV/c); t_{cl} (ns)", jetptbinning, timebinning);
115  fQAHistos->CreateTH2("hClusterConstM02", "EMCAL cluster M02 vs. jet pt; p{t, jet} (GeV/c); M02", jetptbinning, m02binning);
116  fQAHistos->CreateTH2("hClusterConstNcell", "EMCAL cluster ncell vs. jet pt; p{t, jet} (GeV/c); Number of cells", jetptbinning, ncellbinning);
117 
118  // Test of constituent QA
119 #ifdef EXPERIMENTAL_JETCONSTITUENTS
120  fQAHistos->CreateTH2("hChargedConstituentPt", "charged constituent pt vs jet pt (via constituent map); p_{t,jet} (GeV/c); p_{t,ch} (GeV/c)", jetptbinning, clusterenergybinning);
121  fQAHistos->CreateTH2("hChargedIndexPt", "charged constituent pt vs jet pt (via index map); p_{t, jet} (GeV/c); p_{t, ch} (GeV/c)", jetptbinning, clusterenergybinning);
122 
123  fQAHistos->CreateTH2("hClusterConstituentEDefault", "cluster constituent default energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
124  fQAHistos->CreateTH2("hClusterConstituentENLC", "cluster constituent non-linearity-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
125  fQAHistos->CreateTH2("hClusterConstituentEHC", "cluster constituent hadronic-corrected energy vs. jet pt (va constituent map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
126  fQAHistos->CreateTH2("hClusterIndexENLC", "cluster constituent non-linearity-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
127  fQAHistos->CreateTH2("hClusterIndexEHC", "cluster constituent hadronic-corrected energy vs. jet pt (via index map); p_{t, jet} (GeV/c); E_{cl} (GeV)", jetptbinning, clusterenergybinning);
128 #endif
129  for(auto h : *(fQAHistos->GetListOfHistograms())) fOutput->Add(h);
130 
131  OpenFile(2);
132  fJetSubstructureTree = new TTree("jetSubstructure", "Tree with jet substructure information");
133  TString varnames[kTNVar];
134  varnames[0] = "Radius";
135  varnames[1] = "EventWeight";
136  varnames[2] = "PtJetRec";
137  varnames[3] = "PtJetSim";
138  varnames[4] = "EJetRec";
139  varnames[5] = "EJetSim";
140  varnames[6] = "RhoPtRec";
141  varnames[7] = "RhoPtSim";
142  varnames[8] = "RhoMassRec";
143  varnames[9] = "RhoMassSim";
144  varnames[10] = "AreaRec";
145  varnames[11] = "AreaSim";
146  varnames[12] = "NEFRec";
147  varnames[13] = "NEFSim";
148  varnames[14] = "MassRec";
149  varnames[15] = "MassSim";
150  varnames[16] = "ZgMeasured";
151  varnames[17] = "ZgTrue";
152  varnames[18] = "RgMeasured";
153  varnames[19] = "RgTrue";
154  varnames[20] = "MgMeasured";
155  varnames[21] = "MgTrue";
156  varnames[22] = "PtgMeasured";
157  varnames[23] = "PtgTrue";
158  varnames[24] = "MugMeasured";
159  varnames[25] = "MugTrue";
160  varnames[26] = "OneSubjettinessMeasured";
161  varnames[27] = "OneSubjettinessTrue";
162  varnames[28] = "TwoSubjettinessMeasured";
163  varnames[29] = "TwoSubjettinessTrue";
164  varnames[30] = "AngularityMeasured";
165  varnames[31] = "AngularityTrue";
166  varnames[32] = "PtDMeasured";
167  varnames[33] = "PtDTrue";
168  varnames[34] = "NCharged";
169  varnames[35] = "NNeutral";
170  varnames[36] = "NConstTrue";
171  varnames[37] = "NDroppedMeasured";
172  varnames[38] = "NDroppedTrue";
173 
174  for(int ib = 0; ib < kTNVar; ib++){
175  fJetSubstructureTree->Branch(varnames[ib], fJetTreeData + ib, Form("%s/D", varnames[ib].Data()));
176  }
177  PostData(1, fOutput);
178  PostData(2, fJetSubstructureTree);
179 }
180 
182  AliClusterContainer *clusters = GetClusterContainer("caloClusters");
183  AliTrackContainer *tracks = GetTrackContainer("tracks");
184  AliParticleContainer *particles = GetParticleContainer("mcparticles");
185 
186  AliJetContainer *mcjets = GetJetContainer("mcjets");
187  AliJetContainer *datajets = GetJetContainer("datajets");
188 
189  //for(auto e : *(fInputEvent->GetList())) std::cout << e->GetName() << std::endl;
190 
191  TString rhoTagData = datajets ? TString::Format("R%02d", static_cast<Int_t>(datajets->GetJetRadius() * 10.)) : "",
192  rhoTagMC = mcjets ? TString::Format("R%02d", static_cast<Int_t>(mcjets->GetJetRadius() * 10.)) : "";
193 
194  AliRhoParameter *rhoPtRec = GetRhoFromEvent("RhoSparse_Full_" + rhoTagData),
195  *rhoMassRec = GetRhoFromEvent("RhoMassSparse_Full_" + rhoTagData),
196  *rhoPtSim = GetRhoFromEvent("RhoSparse_Full_" + rhoTagMC),
197  *rhoMassSim = GetRhoFromEvent("RhoMassSparse_Full_" + rhoTagMC);
198  AliDebugStream(2) << "Found rho parameter for reconstructed pt: " << (rhoPtRec ? "yes" : "no") << ", value: " << (rhoPtRec ? rhoPtRec->GetVal() : 0.) << std::endl;
199  AliDebugStream(2) << "Found rho parameter for sim pt: " << (rhoPtSim ? "yes" : "no") << ", value: " << (rhoPtSim ? rhoPtSim->GetVal() : 0.) << std::endl;
200  AliDebugStream(2) << "Found rho parameter for reconstructed Mass: " << (rhoMassRec ? "yes" : "no") << ", value: " << (rhoMassRec ? rhoMassRec->GetVal() : 0.) << std::endl;
201  AliDebugStream(2) << "Found rho parameter for sim Mass: " << (rhoMassSim ? "yes" : "no") << ", value: " << (rhoMassSim ? rhoMassSim->GetVal() : 0.) << std::endl;
202 
203  // Run trigger selection (not on pure MCgen train)
204  if(datajets){
205  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
206  if(!mcjets){
207  // Pure data - do EMCAL trigger selection from selection string
208  if(fTriggerSelectionString.Length()) {
209  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
210  }
211  } else {
212  // Simulation - do EMCAL trigger selection from trigger selection object
213  PWG::EMCAL::AliEmcalTriggerDecisionContainer *mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject("EmcalTriggerDecision"));
214  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
215  if(fTriggerSelectionString.Length()){
216  if(!mctrigger){
217  AliErrorStream() << "Trigger decision container not found in event - not possible to select EMCAL triggers" << std::endl;
218  return false;
219  }
220  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
221  }
222  }
223  }
224 
225  Double_t weight = 1;
226  Double_t rhoparameters[4]; memset(rhoparameters, 0, sizeof(Double_t) * 4);
227  if(rhoPtRec) rhoparameters[0] = rhoPtRec->GetVal();
228  if(rhoPtSim) rhoparameters[1] = rhoPtSim->GetVal();
229  if(rhoMassRec) rhoparameters[2] = rhoMassRec->GetVal();
230  if(rhoMassSim) rhoparameters[3] = rhoMassSim->GetVal();
231 
232  AliSoftdropDefinition softdropSettings;
233  softdropSettings.fBeta = fSDBetaCut;
234  softdropSettings.fZ = fSDZCut;
235  switch(fReclusterizer) {
236  case kCAAlgo: softdropSettings.fRecluserAlgo = fastjet::cambridge_aachen_algorithm; break;
237  case kKTAlgo: softdropSettings.fRecluserAlgo = fastjet::kt_algorithm; break;
238  case kAKTAlgo: softdropSettings.fRecluserAlgo = fastjet::antikt_algorithm; break;
239  };
240 
241  AliNSubjettinessDefinition nsubjettinessSettings;
242  nsubjettinessSettings.fBeta = 1.;
243  nsubjettinessSettings.fRadius = 0.4;
244 
245  std::set<AliEmcalJet *> taglist;
246  if(datajets) {
247  AliDebugStream(1) << "In data jets branch: found" << datajets->GetNJets() << " jets, " << datajets->GetNAcceptedJets() << " were accepted\n";
248  AliDebugStream(1) << "Having MC information: " << (mcjets ? TString::Format("yes, with %d jets", mcjets->GetNJets()) : "no") << std::endl;
249  for(auto jet : datajets->accepted()) {
250  double pt = jet->Pt(), pz = jet->Pz(), E = jet->E(), M = TMath::Sqrt(E*E - pt*pt - pz*pz);
251  AliDebugStream(2) << "Next jet: pt:" << jet->Pt() << ", E: " << E << ", pz: " << pz << ", M(self): " << M << "M(fj)" << jet->M() << std::endl;
252  AliEmcalJet *associatedJet = jet->ClosestJet();
253 
254  if(mcjets) {
255  if(!associatedJet) continue;
256  taglist.insert(associatedJet);
257  try {
258  DoConstituentQA(jet, tracks, clusters);
259  AliJetSubstructureData structureData = MakeJetSubstructure(*jet, datajets->GetJetRadius() * 2., tracks, clusters, {softdropSettings, nsubjettinessSettings}),
260  structureMC = MakeJetSubstructure(*associatedJet, mcjets->GetJetRadius() * 2, particles, nullptr, {softdropSettings, nsubjettinessSettings});
261  Double_t angularity[2] = {MakeAngularity(*jet, tracks, clusters), MakeAngularity(*associatedJet, particles, nullptr)},
262  ptd[2] = {MakePtD(*jet, tracks, clusters), MakePtD(*associatedJet, particles, nullptr)};
263  FillTree(datajets->GetJetRadius(), weight, jet, associatedJet, &(structureData.fSoftDrop), &(structureMC.fSoftDrop), &(structureData.fNsubjettiness), &(structureMC.fNsubjettiness), angularity, ptd, rhoparameters);
264  } catch(ReclusterizerException &e) {
265  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
266  }
267  } else {
268  try {
269  DoConstituentQA(jet, tracks, clusters);
270  AliJetSubstructureData structure = MakeJetSubstructure(*jet, 0.4, tracks, clusters, {softdropSettings, nsubjettinessSettings});
271  Double_t angularity[2] = {MakeAngularity(*jet, tracks, clusters), 0.},
272  ptd[2] = {MakePtD(*jet, tracks, clusters), 0.};
273  FillTree(datajets->GetJetRadius(), weight, jet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness), nullptr, angularity, ptd, rhoparameters);
274  } catch(ReclusterizerException &e) {
275  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
276  }
277  }
278  }
279  }
280 
281  if(mcjets) {
282  // process non-matched true jets
283  // - for efficiency studies
284  // - for MCgen train (taglist will be empty in this case)
285  AliDebugStream(1) << "In MC pure jet branch: found " << mcjets->GetNJets() << " jets, " << mcjets->GetNAcceptedJets() << " were accepted\n";
286  for(auto j : mcjets->accepted()){
287  AliEmcalJet *mcjet = static_cast<AliEmcalJet *>(j);
288  if(taglist.find(mcjet) != taglist.end()) continue;
289  try {
290  AliJetSubstructureData structure = MakeJetSubstructure(*mcjet, mcjets->GetJetRadius() * 2., particles, nullptr,{softdropSettings, nsubjettinessSettings});
291  Double_t angularity[2] = {0., MakeAngularity(*mcjet, particles, nullptr)},
292  ptd[2] = {0., MakePtD(*mcjet, particles, nullptr)};
293  FillTree(mcjets->GetJetRadius(), weight, nullptr, mcjet, nullptr, &(structure.fSoftDrop), nullptr, &(structure.fNsubjettiness), angularity, ptd, rhoparameters);
294  } catch (ReclusterizerException &e) {
295  AliErrorStream() << "Error in reclusterization - skipping jet" << std::endl;
296  }
297  }
298  }
299 
300  return true;
301 }
302 
304  const AliEmcalJet *datajet, const AliEmcalJet *mcjet,
305  AliSoftDropParameters *dataSoftdrop, AliSoftDropParameters *mcSoftdrop,
306  AliNSubjettinessParameters *dataSubjettiness, AliNSubjettinessParameters *mcSubjettiness,
307  Double_t *angularity, Double_t *ptd, Double_t *rhoparameters){
308  fJetTreeData[kTRadius] = r;
309  fJetTreeData[kTWeight] = weight;
310  fJetTreeData[kTRhoPtRec] = rhoparameters[0];
311  fJetTreeData[kTRhoPtSim] = rhoparameters[1];
312  fJetTreeData[kTRhoMassRec] = rhoparameters[2];
313  fJetTreeData[kTRhoMassSim] = rhoparameters[3];
314  if(datajet) {
315  fJetTreeData[kTPtJetRec] = TMath::Abs(datajet->Pt());
317  fJetTreeData[kTNNeutral] = datajet->GetNumberOfClusters();
318  fJetTreeData[kTAreaRec] = datajet->Area();
319  fJetTreeData[kTNEFRec] = datajet->NEF();
320  fJetTreeData[kTMassRec] = datajet->M();
321  fJetTreeData[kTEJetRec] = datajet->E();
322  } else {
323  fJetTreeData[kTPtJetRec] = 0.;
326  fJetTreeData[kTAreaRec] = 0.;
327  fJetTreeData[kTNEFRec] = 0.;
328  fJetTreeData[kTMassRec] = 0.;
329  fJetTreeData[kTEJetRec] = 0.;
330  }
331 
332  if(mcjet) {
333  fJetTreeData[kTPtJetSim] = TMath::Abs(mcjet->Pt());
335  fJetTreeData[kTAreaSim] = mcjet->Area();
336  fJetTreeData[kTNEFSim] = mcjet->NEF();
337  fJetTreeData[kTMassSim] = mcjet->M();
338  fJetTreeData[kTEJetSim] = mcjet->E();
339  } else {
340  fJetTreeData[kTPtJetSim] = 0.;
342  fJetTreeData[kTAreaSim] = 0.;
343  fJetTreeData[kTNEFSim] = 0.;
344  fJetTreeData[kTMassSim] = 0.;
345  fJetTreeData[kTEJetSim] = 0.;
346  }
347 
348  if(dataSoftdrop) {
349  fJetTreeData[kTZgMeasured] = dataSoftdrop->fZg;
350  fJetTreeData[kTRgMeasured] = dataSoftdrop->fRg;
351  fJetTreeData[kTMgMeasured] = dataSoftdrop->fMg;
352  fJetTreeData[kTPtgMeasured] = dataSoftdrop->fPtg;
353  fJetTreeData[kTMugMeasured] = dataSoftdrop->fMug;
354  fJetTreeData[kTNDroppedMeasured] = dataSoftdrop->fNDropped;
355  } else {
362  }
363 
364  if(mcSoftdrop) {
365  fJetTreeData[kTZgTrue] = mcSoftdrop->fZg;
366  fJetTreeData[kTRgTrue] = mcSoftdrop->fRg;
367  fJetTreeData[kTMgTrue] = mcSoftdrop->fMg;
368  fJetTreeData[kTPtgTrue] = mcSoftdrop->fPtg;
369  fJetTreeData[kTMugTrue] = mcSoftdrop->fMug;
370  fJetTreeData[kTNDroppedTrue] = mcSoftdrop->fNDropped;
371  } else {
372  fJetTreeData[kTZgTrue] = 0.;
373  fJetTreeData[kTRgTrue] = 0.;
374  fJetTreeData[kTMgTrue] = 0.;
375  fJetTreeData[kTPtgTrue] = 0.;
376  fJetTreeData[kTMugTrue] = 0.;
378  }
379 
380  if(dataSubjettiness) {
383  } else {
386  }
387 
388  if(mcSubjettiness) {
391  } else {
394  }
395 
396  fJetTreeData[kTAngularityMeasured] = angularity[0];
397  fJetTreeData[kTAngularityTrue] = angularity[1];
398  fJetTreeData[kTPtDMeasured] = ptd[0];
399  fJetTreeData[kTPtDTrue] = ptd[1];
400 
401  fJetSubstructureTree->Fill();
402 }
403 
404 
406  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
407  std::vector<fastjet::PseudoJet> constituents;
408  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
409  AliVTrack *track = static_cast<AliVTrack *>(jet.TrackAt(itrk, tracks->GetArray()));
410  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
411  constituentTrack.set_user_index(jet.TrackAt(itrk));
412  constituents.push_back(constituentTrack);
413  }
414 
415  if(clusters){
416  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
417  AliVCluster *cluster = jet.ClusterAt(icl, clusters->GetArray());
418  TLorentzVector clustervec;
419  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
420  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
421  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
422  constituents.push_back(constituentCluster);
423  }
424  }
425 
426  // Redo jet finding on constituents with a
427  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
428  std::vector<fastjet::PseudoJet> outputjets;
429  try {
430  fastjet::ClusterSequence jetfinder(constituents, jetdef);
431  outputjets = jetfinder.inclusive_jets(0);
432  AliJetSubstructureData result({MakeSoftDropParameters(outputjets[0], settings.fSoftdropSettings), MakeNsubjettinessParameters(outputjets[0], settings.fSubjettinessSettings)});
433  return result;
434  } catch (fastjet::Error &e) {
435  AliErrorStream() << " FJ Exception caught: " << e.message() << std::endl;
436  throw ReclusterizerException();
437  }
438 }
439 
441  fastjet::contrib::SoftDrop softdropAlgorithm(cutparameters.fBeta, cutparameters.fZ);
442  softdropAlgorithm.set_verbose_structure(kTRUE);
443  std::unique_ptr<fastjet::contrib::Recluster> reclusterizer(new fastjet::contrib::Recluster(cutparameters.fRecluserAlgo, 1, true));
444  softdropAlgorithm.set_reclustering(kTRUE, reclusterizer.get());
445  fastjet::PseudoJet groomed = softdropAlgorithm(jet);
446 
447  AliSoftDropParameters result({groomed.structure_of<fastjet::contrib::SoftDrop>().symmetry(),
448  groomed.m(),
449  groomed.structure_of<fastjet::contrib::SoftDrop>().delta_R(),
450  groomed.perp(),
451  groomed.structure_of<fastjet::contrib::SoftDrop>().mu(),
452  groomed.structure_of<fastjet::contrib::SoftDrop>().dropped_count()});
453  return result;
454 }
455 
458  fastjet::contrib::Nsubjettiness (1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet),
459  fastjet::contrib::Nsubjettiness (2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedCutoffMeasure(cut.fBeta, cut.fRadius, 1e100)).result(jet)
460  });
461  return result;
462 }
463 
465  if(!jet.GetNumberOfTracks()) return 0;
466  TVector3 jetvec(jet.Px(), jet.Py(), jet.Pz());
467  Double_t den(0.), num(0.);
468  if(tracks){
469  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
470  AliVParticle *track = jet.TrackAt(itrk, tracks->GetArray());
471  if(!track){
472  AliErrorStream() << "Associated constituent particle / track not found\n";
473  continue;
474  }
475  TVector3 trackvec(track->Px(), track->Py(), track->Pz());
476 
477  num += track->Pt() * trackvec.DrEtaPhi(jetvec);
478  den += +track->Pt();
479  }
480  }
481  if(clusters) {
482  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
483  AliVCluster *clust = jet.ClusterAt(icl, clusters->GetArray());
484  if(!clust) {
485  AliErrorStream() << "Associated constituent cluster not found\n";
486  continue;
487  }
488  TLorentzVector clusterp;
489  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
490 
491  num += clusterp.Pt() * clusterp.Vect().DrEtaPhi(jetvec);
492  den += clusterp.Pt();
493  }
494  }
495  return num/den;
496 }
497 
499  if (!jet.GetNumberOfTracks()) return 0;
500  Double_t den(0.), num(0.);
501  if(particles){
502  for(UInt_t itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++) {
503  AliVParticle *trk = jet.TrackAt(itrk, particles->GetArray());
504  if(!trk){
505  AliErrorStream() << "Associated constituent particle / track not found\n";
506  continue;
507  }
508  num += trk->Pt() * trk->Pt();
509  den += trk->Pt();
510  }
511  }
512  if(clusters){
513  for(UInt_t icl = 0; icl < jet.GetNumberOfClusters(); icl++){
514  AliVCluster *clust = jet.ClusterAt(icl, clusters->GetArray());
515  if(!clust) {
516  AliErrorStream() << "Associated constituent cluster not found\n";
517  continue;
518  }
519  TLorentzVector clusterp;
520  clust->GetMomentum(clusterp, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
521  num += clusterp.Pt() * clusterp.Pt();
522  den += clusterp.Pt();
523  }
524  }
525  return TMath::Sqrt(num)/den;
526 }
527 
529  for(int icl = 0; icl < jet->GetNumberOfClusters(); icl++){
530  AliVCluster *clust = jet->ClusterAt(icl, clusters->GetArray());
531  AliDebugStream(3) << "cluster time " << clust->GetTOF() << std::endl;
532  fQAHistos->FillTH2("hClusterConstE", jet->Pt(),clust->GetUserDefEnergy(clusters->GetDefaultClusterEnergy()));
533  fQAHistos->FillTH2("hClusterConstTime", jet->Pt(), clust->GetTOF());
534  fQAHistos->FillTH2("hClusterConstM02", jet->Pt(), clust->GetM02());
535  fQAHistos->FillTH2("hClusterConstNcell", jet->Pt(), clust->GetNCells());
536 
537 #ifdef EXPERIMENTAL_JETCONSTITUENTS
538  fQAHistos->FillTH2("hClusterIndexENLC", jet->Pt(), clust->GetNonLinCorrEnergy());
539  fQAHistos->FillTH2("hClusterIndexEHC", jet->Pt(), clust->GetHadCorrEnergy());
540 #endif
541  }
542 
543 #ifdef EXPERIMENTAL_JETCONSTITUENTS
544  // Loop over charged particles - fill test histogram
545  for(int itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
546  AliVParticle *part = jet->TrackAt(itrk, cont->GetArray());
547  fQAHistos->FillTH2("hChargedIndexPt", jet->Pt(), part->Pt());
548  }
549 
550  // Look over charged constituents
551  AliDebugStream(2) << "Jet: Number of particle constituents: " << jet->GetParticleConstituents().GetEntriesFast() << std::endl;
552  for(auto pconst : jet->GetParticleConstituents()) {
553  PWG::JETFW::AliEmcalParticleJetConstituent *part = static_cast<PWG::JETFW::AliEmcalParticleJetConstituent *>(pconst);
554  AliDebugStream(3) << "Found particle constituent with pt " << part->Pt() << ", from VParticle " << part->GetParticle()->Pt() << std::endl;
555  fQAHistos->FillTH2("hChargedConstituentPt", jet->Pt(), part->Pt());
556  }
557 
558  // Loop over neutral constituents
559  AliDebugStream(2) << "Jet: Number of cluster constituents: " << jet->GetClusterConstituents().GetEntriesFast() << std::endl;
560  for(auto cconst : jet->GetClusterConstituents()){
561  PWG::JETFW::AliEmcalClusterJetConstituent *clust = static_cast<PWG::JETFW::AliEmcalClusterJetConstituent *>(cconst);
562  AliDebugStream(3) << "Found cluster constituent with energy " << clust->E() << " using energy definition " << static_cast<int>(clust->GetDefaultEnergyType()) << std::endl;
563  fQAHistos->FillTH2("hClusterConstituentEDefault", jet->Pt(), clust->E());
564  fQAHistos->FillTH2("hClusterConstituentENLC", jet->Pt(), clust->GetCluster()->GetNonLinCorrEnergy());
565  fQAHistos->FillTH2("hClusterConstituentEHC", jet->Pt(), clust->GetCluster()->GetHadCorrEnergy());
566  }
567 #endif
568 }
569 
571  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
572 
573  Bool_t isAOD(kFALSE);
574  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
575  if(inputhandler) {
576  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = kTRUE;
577  }
578 
579  AliAnalysisTaskEmcalJetSubstructureTree *treemaker = new AliAnalysisTaskEmcalJetSubstructureTree("JetSubstructureTreemaker_" + TString::Format("R%02d_", int(jetradius * 10.)) + trigger);
580  mgr->AddTask(treemaker);
581  treemaker->SetMakeGeneralHistograms(kTRUE);
582 
583  // Adding containers
584  if(isMC) {
585  AliParticleContainer *particles = treemaker->AddMCParticleContainer("mcparticles");
586  particles->SetMinPt(0.);
587 
588  AliJetContainer *mcjets = treemaker->AddJetContainer(
591  recombinationScheme,
592  jetradius,
594  particles, nullptr);
595  mcjets->SetName("mcjets");
596  mcjets->SetJetPtCut(20.);
597  }
598 
599  if(isData) {
601  tracks->SetMinPt(0.15);
602  AliClusterContainer *clusters(nullptr);
603  if(jettype == kFull){
604  std::cout << "Using full jets ..." << std::endl;
606  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
607  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
608  } else {
609  std::cout << "Using charged jets ... " << std::endl;
610  }
611 
612  AliJetContainer *datajets = treemaker->AddJetContainer(
615  recombinationScheme,
616  jetradius,
618  tracks, clusters);
619  datajets->SetName("datajets");
620  datajets->SetJetPtCut(20.);
621 
622  treemaker->SetUseAliAnaUtils(true, true);
623  treemaker->SetVzRange(-10., 10);
624 
625  // configure trigger selection
626  TString triggerstring(trigger);
627  if(triggerstring.Contains("INT7")) {
628  treemaker->SetTriggerBits(AliVEvent::kINT7);
629  } else if(triggerstring.Contains("EJ1")) {
630  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
631  treemaker->SetTriggerString("EJ1");
632  } else if(triggerstring.Contains("EJ2")) {
633  treemaker->SetTriggerBits(AliVEvent::kEMCEJE);
634  treemaker->SetTriggerString("EJ2");
635  }
636  }
637 
638  // Connecting containers
639  TString outputfile = mgr->GetCommonFileName();
640  outputfile += TString::Format(":JetSubstructure_R%02d_%s", int(jetradius * 10.), trigger);
641  mgr->ConnectInput(treemaker, 0, mgr->GetCommonInputContainer());
642  mgr->ConnectOutput(treemaker, 1, mgr->CreateContainer("JetSubstructureHistos_" + TString::Format("R%0d_", int(jetradius * 10.)) + trigger, AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile));
643  mgr->ConnectOutput(treemaker, 2, mgr->CreateContainer("JetSubstuctureTree_" + TString::Format("R%0d_", int(jetradius * 10.)) + trigger, TTree::Class(), AliAnalysisManager::kOutputContainer, Form("JetSubstructureTree_R%02d_%s.root", static_cast<int>(jetradius*10.), trigger)));
644 
645  return treemaker;
646 }
647 
648 } /* namespace EmcalTriggerJets */
Double_t Area() const
Definition: AliEmcalJet.h:123
bool IsEventSelected(const char *name) const
Checks whether the events is selected for a given trigger type.
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
Double_t MakePtD(const AliEmcalJet &jet, const AliParticleContainer *const particles, const AliClusterContainer *const clusters) const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:222
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Py() const
Definition: AliEmcalJet.h:100
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
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
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.
Double_t E() const
Definition: AliEmcalJet.h:112
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 GetDefaultClusterEnergy() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
static AliAnalysisTaskEmcalJetSubstructureTree * AddEmcalJetSubstructureTreeMaker(Bool_t isMC, Bool_t isData, Double_t jetradius, JetType_t jettype, AliJetContainer::ERecoScheme_t recombinationScheme, const char *name)
AliRhoParameter * GetRhoFromEvent(const char *name)
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.
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:131
void SetJetPtCut(Float_t cut)
TPC acceptance.
Definition: AliEmcalJet.h:60
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Int_t GetNJets() const
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 MakeAngularity(const AliEmcalJet &jet, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const
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
Bool_t isMC
void DoConstituentQA(const AliEmcalJet *jet, const AliParticleContainer *tracks, const AliClusterContainer *clusters)
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.
void FillTree(double r, double weight, const AliEmcalJet *datajet, const AliEmcalJet *mcjet, AliSoftDropParameters *dataSoftdrop, AliSoftDropParameters *mcsoftdrop, AliNSubjettinessParameters *dataSubjettiness, AliNSubjettinessParameters *mcSubjettiness, Double_t *angularity, Double_t *ptd, Double_t *rhoparameters)
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Container class for histograms.
Definition: THistManager.h:99
Double_t Pz() const
Definition: AliEmcalJet.h:101
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
void SetDefaultClusterEnergy(Int_t d)
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
void SetClusHadCorrEnergyCut(Double_t cut)
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.