AliPhysics  b0c77bb (b0c77bb)
AliAnalysisTaskEmcalSoftDropData.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2019, 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 <memory>
28 #include <sstream>
29 
30 #include <fastjet/ClusterSequence.hh>
31 #include <fastjet/JetDefinition.hh>
32 #include <fastjet/PseudoJet.hh>
33 #include <fastjet/contrib/SoftDrop.hh>
34 
35 #include <TArray.h>
36 #include <TCustomBinning.h>
37 #include <THistManager.h>
38 
40 #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 "AliVTrack.h"
52 
54 
55 using namespace PWGJE::EMCALJetTasks;
56 
59  fBinningMode(kSDModeINT7),
60  fTriggerBits(AliVEvent::kAny),
61  fTriggerString(""),
62  fUseDownscaleWeight(false),
63  fBeta(0.),
64  fZcut(0.1),
65  fReclusterizer(kCAAlgo),
66  fUseChargedConstituents(kTRUE),
67  fUseNeutralConstituents(kTRUE),
68  fHistos(nullptr),
69  fPtBinning(nullptr)
70 {
71 
72 }
73 
74 AliAnalysisTaskEmcalSoftDropData::AliAnalysisTaskEmcalSoftDropData(const char *name) :
75  AliAnalysisTaskEmcalJet(name, kTRUE),
76  fBinningMode(kSDModeINT7),
77  fTriggerBits(AliVEvent::kAny),
78  fTriggerString(""),
79  fUseDownscaleWeight(false),
80  fBeta(0.),
81  fZcut(0.1),
82  fReclusterizer(kCAAlgo),
83  fUseChargedConstituents(kTRUE),
84  fUseNeutralConstituents(kTRUE),
85  fHistos(nullptr),
86  fPtBinning(nullptr)
87 {
88 
89 }
90 
92  if(fPtBinning) delete fPtBinning;
93  if(fHistos) delete fHistos;
94 }
95 
98 
100  std::unique_ptr<TBinning> zgBinning(GetZgBinning());
101  TArrayD edgesPt;
102  fPtBinning->CreateBinEdges(edgesPt);
103  fJetPtMin = edgesPt[0];
104  fJetPtMax = edgesPt[edgesPt.GetSize()];
105 
106  fHistos = new THistManager("histosSoftdrop");
107  fHistos->CreateTH1("hEventCounter", "EventCounter", 1, 0., 1);
108  fHistos->CreateTH2("hZgVsPt", "zg vs pt", *zgBinning, *fPtBinning);
110  fHistos->CreateTH2("hZgVsPtWeighted", "zg vs pt", *zgBinning, *fPtBinning);
111  fHistos->CreateTH1("hEventCounterWeighted", "Event counter, weighted", 1., 0., 1.);
112  }
113 
114  for(auto h : *fHistos->GetListOfHistograms()) fOutput->Add(h);
115  PostData(1, fOutput);
116 }
117 
119  if(!(fInputHandler->IsEventSelected() & fTriggerBits)) return false;
120  if(fTriggerString.Length()) {
121  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerString)) return false;
122  }
123  return true;
124 }
125 
129 }
130 
132  auto jets = GetJetContainer("datajets");
133  if(!jets) {
134  AliErrorStream() << "Jet container not found" << std::endl;
135  return false;
136  }
137  auto clusters = GetClusterContainer(0);
138  if(fUseNeutralConstituents && !clusters) {
139  AliErrorStream() << "Cluster container not found, but neutral constituents requested" << std::endl;
140  }
141  auto tracks = GetTrackContainer(0);
142  if(fUseChargedConstituents && !tracks) {
143  AliErrorStream() << "Track container not found, but charged constituent requested." << std::endl;
144  return false;
145  }
146 
148  fHistos->FillTH1("hEventCounter", 1., weight);
149  if(fUseDownscaleWeight) fHistos->FillTH1("hEventCounterWeighted", 1., weight);
150 
151  for(auto jet : jets->accepted()){
152  if(jet->Pt() < fJetPtMin || jet->Pt() > fJetPtMax) continue;
153  auto zgparams = MakeSoftdrop(*jet, jets->GetJetRadius(), tracks, clusters);
154  fHistos->FillTH2("hZgVsPt", zgparams[0], jet->Pt());
155  if(fUseDownscaleWeight) fHistos->FillTH2("hZgVsPtWeighted", zgparams[0], jet->Pt(), weight);
156  }
157 }
158 
160  Double_t weight = 1.;
161  TString triggerclass;
162  if(fTriggerString == "INT7") triggerclass = "CINT7-B-NOPF-CENT";
163  else if(fTriggerString == "EJ1") triggerclass = "CEMC7EJ1-B-NOPF-CENTNOTRD";
164  else if(fTriggerString == "EJ2") triggerclass = "CEMC7EJ1-B-NOPF-CENT";
165  if(triggerclass.Length()) weight = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance()->GetDownscaleFactorForTriggerClass(triggerclass);
166  return weight;
167 }
168 
170  auto binning = new TCustomBinning;
171  switch(fBinningMode) {
172  case kSDModeINT7: {
173  binning->SetMinimum(20);
174  binning->AddStep(40., 5.);
175  binning->AddStep(60., 10.);
176  binning->AddStep(80., 20.);
177  binning->AddStep(120., 40.);
178  break;
179  }
180  case kSDModeEJ1: {
181  binning->SetMinimum(80.);
182  binning->AddStep(120., 5.);
183  binning->AddStep(160., 10.);
184  binning->AddStep(200., 20.);
185  break;
186  }
187  case kSDModeEJ2: {
188  binning->SetMinimum(70.);
189  binning->AddStep(100., 5.);
190  binning->AddStep(120., 10.);
191  binning->AddStep(140., 20.);
192  break;
193  }
194  };
195  return binning;
196 }
197 
199  auto binning = new TCustomBinning;
200  binning->SetMinimum(0.);
201  binning->AddStep(fZcut, fZcut);
202  binning->AddStep(0.5, 0.05);
203  return binning;
204 }
205 
206 std::vector<double> AliAnalysisTaskEmcalSoftDropData::MakeSoftdrop(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const {
207  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
208  std::vector<fastjet::PseudoJet> constituents;
209  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
210  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
211  if(tracks && (fUseChargedConstituents || isMC)){ // Neutral particles part of particle container in case of MC
212  AliDebugStream(1) << "Jet substructure: Using charged constituents" << std::endl;
213  for(int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++){
214  auto track = jet.TrackAt(itrk, tracks->GetArray());
215  if(!track->Charge() && !fUseNeutralConstituents) continue; // Reject neutral constituents in case of using only charged consituents
216  if(track->Charge() && !fUseChargedConstituents) continue; // Reject charged constituents in case of using only neutral consituents
217  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
218  constituentTrack.set_user_index(jet.TrackAt(itrk));
219  constituents.push_back(constituentTrack);
220  }
221  }
222 
223  if(clusters && fUseNeutralConstituents){
224  AliDebugStream(1) << "Jet substructure: Using neutral constituents" << std::endl;
225  for(int icl = 0; icl < jet.GetNumberOfClusters(); icl++) {
226  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
227  TLorentzVector clustervec;
228  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
229  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
230  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
231  constituents.push_back(constituentCluster);
232  }
233  }
234 
235  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
236  if(!constituents.size()){
237  AliErrorStream() << "Jet has 0 constituents." << std::endl;
238  throw 1;
239  }
240  // Redo jet finding on constituents with a
241  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
242  fastjet::ClusterSequence jetfinder(constituents, jetdef);
243  std::vector<fastjet::PseudoJet> outputjets = jetfinder.inclusive_jets(0);
244  auto sdjet = outputjets[0];
245  fastjet::contrib::SoftDrop softdropAlgorithm(fBeta, fZcut);
246  softdropAlgorithm.set_verbose_structure(kTRUE);
247  fastjet::JetAlgorithm reclusterizingAlgorithm;
248  switch(fReclusterizer) {
249  case kCAAlgo: reclusterizingAlgorithm = fastjet::cambridge_aachen_algorithm; break;
250  case kKTAlgo: reclusterizingAlgorithm = fastjet::kt_algorithm; break;
251  case kAKTAlgo: reclusterizingAlgorithm = fastjet::antikt_algorithm; break;
252  };
253 #if FASTJET_VERSION_NUMBER >= 30302
254  fastjet::Recluster reclusterizer(reclusterizingAlgorithm, 1, fastjet::Recluster::keep_only_hardest);
255 #else
256  fastjet::contrib::Recluster reclusterizer(reclusterizingAlgorithm, 1, true);
257 #endif
258  softdropAlgorithm.set_reclustering(kTRUE, &reclusterizer);
259  AliDebugStream(4) << "Jet has " << sdjet.constituents().size() << " constituents" << std::endl;
260  auto groomed = softdropAlgorithm(sdjet);
261  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
262 
263  std::vector<double> result ={softdropstruct.symmetry(),
264  groomed.m(),
265  softdropstruct.delta_R(),
266  groomed.perp(),
267  softdropstruct.mu(),
268  static_cast<double>(softdropstruct.dropped_count())};
269  return result;
270 }
271 
273  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
274 
275  Bool_t isAOD(kFALSE);
276  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
277  if(inputhandler) {
278  if(inputhandler->IsA() == AliAODInputHandler::Class()){
279  std::cout << "Analysing AOD events\n";
280  isAOD = kTRUE;
281  } else {
282  std::cout << "Analysing ESD events\n";
283  }
284  }
285 
286  std::stringstream taskname;
287  taskname << "SoftdropDataMaker_R" << std::setw(2) << std::setfill('0') << int(jetradius*10) << trigger;
288  AliAnalysisTaskEmcalSoftDropData *datamaker = new AliAnalysisTaskEmcalSoftDropData(taskname.str().data());
289  datamaker->SelectCollisionCandidates(AliVEvent::kINT7);
290  mgr->AddTask(datamaker);
291 
292  AliTrackContainer *tracks(nullptr);
293  if((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet)){
295  std::cout << "Track container name: " << tracks->GetName() << std::endl;
296  tracks->SetMinPt(0.15);
297  }
298  AliClusterContainer *clusters(nullptr);
299  if((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet)){
300  std::cout << "Using full or neutral jets ..." << std::endl;
302  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
303  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
304  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
305  } else {
306  std::cout << "Using charged jets ... " << std::endl;
307  }
308 
309  AliJetContainer *datajets = datamaker->AddJetContainer(
310  jettype,
312  recombinationScheme,
313  jetradius,
315  tracks, clusters);
316  datajets->SetName("detLevel");
317  datajets->SetJetPtCut(0.);
318  datajets->SetMaxTrackPt(1000.);
319 
320  std::string jettypestring;
321  switch(jettype) {
322  case AliJetContainer::kFullJet: jettypestring = "FullJets"; break;
323  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; break;
324  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; break;
325  default: jettypestring = "Undef";
326  };
327 
328  EBinningMode_t binmode(kSDModeINT7);
329  UInt_t triggerbits(AliVEvent::kINT7);
330  std::string triggerstring(trigger);
331  if(triggerstring == "EJ1") {
332  binmode = kSDModeEJ1;
333  triggerbits = AliVEvent::kEMCEJE;
334  } else if(triggerstring == "EJ2") {
335  binmode = kSDModeEJ2;
336  triggerbits = AliVEvent::kEMCEJE;
337  }
338  datamaker->SetBinningMode(binmode);
339 
340  // Connecting containers
341  std::stringstream outputfile, histname;
342  outputfile << mgr->GetCommonFileName() << ":SoftDropResponse_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
343  histname << "SoftDropResponseHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
344  mgr->ConnectInput(datamaker, 0, mgr->GetCommonInputContainer());
345  mgr->ConnectOutput(datamaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
346 
347  return datamaker;
348 }
double Double_t
Definition: External.C:58
AliJetContainer * GetJetContainer(Int_t i=0) const
virtual Bool_t IsTriggerSelected()
Selection of a hardware trigger.
Container with name, TClonesArray and cuts for particles.
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:137
static AliEmcalDownscaleFactorsOCDB * Instance()
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
virtual void CreateBinEdges(TArrayD &binedges) const =0
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
Interface for binnings used by the histogram handler.
Definition: TBinning.h:23
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
Int_t GetDefaultClusterEnergy() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:68
static AliAnalysisTaskEmcalSoftDropData * AddTaskEmcalSoftDropData(Double_t jetradius, AliJetContainer::EJetType_t jettype, AliJetContainer::ERecoScheme_t recombinationScheme, const char *trigger)
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.
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
void SetJetPtCut(Float_t cut)
virtual Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
virtual void RunChanged(Int_t newrun)
Process tasks relevant when a file with a different run number is processed.
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Helper class creating user defined custom binning.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t Pt() const
Definition: AliEmcalJet.h:109
Bool_t isMC
AliEmcalList * fOutput
!output list
std::vector< double > MakeSoftdrop(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const
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
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
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
void SetDefaultClusterEnergy(Int_t d)
void SetMaxTrackPt(Float_t b)
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
void SetClusHadCorrEnergyCut(Double_t cut)
void SetMinimum(Double_t min)
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.