AliPhysics  c6e65cb (c6e65cb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetHCorrelations.cxx
Go to the documentation of this file.
1 //Measure Jet-hadron correlations
3 //Does event Mixing using AliEventPoolManager
5 
7 
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <THnSparse.h>
12 #include <TVector3.h>
13 #include <TFile.h>
14 #include <TGrid.h>
15 
16 #include <AliAnalysisManager.h>
17 #include <AliInputEventHandler.h>
18 #include <AliEventPoolManager.h>
19 #include <AliLog.h>
20 #include <AliVAODHeader.h>
21 #include <AliVTrack.h>
22 
23 #include "AliEmcalJet.h"
24 #include "AliTLorentzVector.h"
25 #include "AliBasicParticle.h"
26 #include "AliEmcalContainerUtils.h"
27 #include "AliClusterContainer.h"
28 #include "AliTrackContainer.h"
29 #include "AliJetContainer.h"
31 
35 
36 // 0-10% centrality: Semi-Good Runs
37 Double_t AliAnalysisTaskEmcalJetHCorrelations::p0_10SG[17] = {0.906767, 0.0754127, 1.11638, -0.0233078, 0.795454, 0.00935385, -0.000327857, 1.08903, 0.0107272, 0.443252, -0.143411, 0.965822, 0.359156, -0.581221, 1.0739, 0.00632828, 0.706356};
38 // 10-30% centrality: Semi-Good Runs
39 Double_t AliAnalysisTaskEmcalJetHCorrelations::p10_30SG[17] = {0.908011, 0.0769254, 1.11912, -0.0249449, 0.741488, 0.0361252, -0.00367954, 1.10424, 0.011472, 0.452059, -0.133282, 0.980633, 0.358222, -0.620256, 1.06871, 0.00564449, 0.753168};
40 // 30-50% centrality: Semi-Good Runs
41 Double_t AliAnalysisTaskEmcalJetHCorrelations::p30_50SG[17] = {0.958708, 0.0799197, 1.10817, -0.0357678, 0.75051, 0.0607808, -0.00929713, 0.998801, 0.00692244, 0.615452, -0.0480328, 0.968431, 0.321634, -0.619066, 1.03412, 0.00656201, 0.798666};
42 // 50-90% centrality: Semi-Good Runs
43 Double_t AliAnalysisTaskEmcalJetHCorrelations::p50_90SG[17] = {0.944565, 0.0807258, 1.12709, -0.0324746, 0.666452, 0.0842476, -0.00963837, 1.02829, 0.00666852, 0.549625, -0.0603107, 0.981374, 0.309374, -0.619181, 1.05367, 0.005925, 0.744887};
44 
45 // 0-10% centrality: Good Runs
46 Double_t AliAnalysisTaskEmcalJetHCorrelations::p0_10G[17] = {0.971679, 0.0767571, 1.13355, -0.0274484, 0.856652, 0.00536795, 3.90795e-05, 1.06889, 0.011007, 0.447046, -0.146626, 0.919777, 0.192601, -0.268515, 1.00243, 0.00620849, 0.709477};
47 // 10-30% centrality: Good Runs
48 Double_t AliAnalysisTaskEmcalJetHCorrelations::p10_30G[17] = {0.97929, 0.0776039, 1.12213, -0.0300645, 0.844722, 0.0134788, -0.0012333, 1.07955, 0.0116835, 0.456608, -0.132743, 0.930964, 0.174175, -0.267154, 0.993118, 0.00574892, 0.765256};
49 // 30-50% centrality: Good Runs
50 Double_t AliAnalysisTaskEmcalJetHCorrelations::p30_50G[17] = {0.997696, 0.0816769, 1.14341, -0.0353734, 0.752151, 0.0744259, -0.0102926, 1.01561, 0.00713274, 0.57203, -0.0640248, 0.947747, 0.102007, -0.194698, 0.999164, 0.00568476, 0.7237};
51 // 50-90% centrality: Good Runs
52 Double_t AliAnalysisTaskEmcalJetHCorrelations::p50_90G[17] = {0.97041, 0.0813559, 1.12151, -0.0368797, 0.709327, 0.0701501, -0.00784043, 1.06276, 0.00676173, 0.53607, -0.0703117, 0.982534, 0.0947881, -0.18073, 1.03229, 0.00580109, 0.737801};
53 
58  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetHCorrelations", kFALSE),
59  fTrackBias(5),
60  fClusterBias(5),
61  fDoEventMixing(kFALSE),
62  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
63  fPoolMgr(nullptr),
64  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
65  fDisableFastPartition(kFALSE),
66  fDoEffCorrection(0),
67  fNoMixedEventJESCorrection(kFALSE),
68  fJESCorrectionHist(nullptr),
69  fDoLessSparseAxes(kFALSE), fDoWiderTrackBin(kFALSE),
70  fHistTrackPt(nullptr),
71  fHistJetEtaPhi(nullptr),
72  fHistJetHEtaPhi(nullptr),
73  fHistJHPsi(nullptr),
74  fhnMixedEvents(nullptr),
75  fhnJH(nullptr)
76 {
77  // Default Constructor
79 }
80 
85  AliAnalysisTaskEmcalJet(name, kTRUE),
86  fTrackBias(5),
87  fClusterBias(5),
88  fDoEventMixing(kFALSE),
89  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
90  fPoolMgr(nullptr),
91  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
92  fDisableFastPartition(kFALSE),
93  fDoEffCorrection(0),
94  fNoMixedEventJESCorrection(kFALSE),
95  fJESCorrectionHist(nullptr),
96  fDoLessSparseAxes(kFALSE), fDoWiderTrackBin(kFALSE),
97  fHistTrackPt(nullptr),
98  fHistJetEtaPhi(nullptr),
99  fHistJetHEtaPhi(nullptr),
100  fHistJHPsi(nullptr),
101  fhnMixedEvents(nullptr),
102  fhnJH(nullptr)
103 {
104  // Constructor
106  // Ensure that additional general histograms are created
108 }
109 
114 {
115  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; trackPtBin++){
116  fHistTrackEtaPhi[trackPtBin] = nullptr;
117  }
118  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
119  fHistJetPt[centralityBin] = nullptr;
120  fHistJetPtBias[centralityBin] = nullptr;
121  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
122  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
123  fHistJetH[centralityBin][jetPtBin][etaBin] = nullptr;
124  fHistJetHBias[centralityBin][jetPtBin][etaBin] = nullptr;
125  }
126  }
127  }
128 }
129 
134  // Called once
136 
137  // Create histograms
138  fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
139  fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
140  fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
141 
142  fHistJHPsi = new TH3F("fHistJHPsi","Jet-Hadron ntr-trpt-dpsi",20,0,100,200,0,20,120,0,180);
143 
144  fOutput->Add(fHistTrackPt);
145  fOutput->Add(fHistJetEtaPhi);
146  fOutput->Add(fHistJetHEtaPhi);
147  fOutput->Add(fHistJHPsi);
148 
149  TString name;
150 
151  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; ++trackPtBin){
152  name = Form("fHistTrackEtaPhi_%i", trackPtBin);
153  fHistTrackEtaPhi[trackPtBin] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
154  fOutput->Add(fHistTrackEtaPhi[trackPtBin]);
155  }
156 
157  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
158  name = Form("fHistJetPt_%i",centralityBin);
159  fHistJetPt[centralityBin] = new TH1F(name,name,200,0,200);
160  fOutput->Add(fHistJetPt[centralityBin]);
161 
162  name = Form("fHistJetPtBias_%i",centralityBin);
163  fHistJetPtBias[centralityBin] = new TH1F(name,name,200,0,200);
164  fOutput->Add(fHistJetPtBias[centralityBin]);
165 
166  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
167  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
168  name = Form("fHistJetH_%i_%i_%i",centralityBin,jetPtBin,etaBin);
169  fHistJetH[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
170  fOutput->Add(fHistJetH[centralityBin][jetPtBin][etaBin]);
171 
172  name = Form("fHistJetHBias_%i_%i_%i",centralityBin,jetPtBin,etaBin);
173  fHistJetHBias[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
174  fOutput->Add(fHistJetHBias[centralityBin][jetPtBin][etaBin]);
175  }
176  }
177  }
178 
179  UInt_t cifras = 0; // bit coded, see GetDimParams() below
180  if(fDoLessSparseAxes) {
181  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
182  } else {
183  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
184  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
185  }
186  fhnJH = NewTHnSparseF("fhnJH", cifras);
187  fhnJH->Sumw2();
188  fOutput->Add(fhnJH);
189 
190  if(fDoEventMixing){
191  if(fDoLessSparseAxes) {
192  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
193  } else {
194  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
195  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
196  }
197  fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
198  fhnMixedEvents->Sumw2();
199  fOutput->Add(fhnMixedEvents);
200  }
201 
202  PostData(1, fOutput);
203 
204  // Event Mixing
205  Int_t poolSize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
206  // ZVertex
207  Int_t nZVertexBins = 10;
208  Double_t* zVertexBins = GenerateFixedBinArray(nZVertexBins, -10, 10);
209  // Event activity (centrality of multiplicity)
210  Int_t nEventActivityBins = 8;
211  Double_t* eventActivityBins = 0;
212  // +1 to accomodate the fact that we define bins rather than array entries.
213  Double_t multiplicityBins[kMixedEventMulitplictyBins+1] = {0., 4., 9., 15., 25., 35., 55., 100., 500.};
214 
215  // Cannot use GetBeamType() since it is not available until UserExec()
216  if (fForceBeamType != AliAnalysisTaskEmcal::kpp ) { //all besides pp
217  // Event Activity is centrality in AA, pA
218  nEventActivityBins = fNCentBinsMixedEvent;
219  eventActivityBins = GenerateFixedBinArray(nEventActivityBins, 0, 100);
220  }
221  else if (fForceBeamType == AliAnalysisTaskEmcal::kpp) { //for pp only
222  // Event Activity is multiplicity in pp
223  eventActivityBins = multiplicityBins;
224  }
225 
226  fPoolMgr = new AliEventPoolManager(poolSize, fNMixingTracks, nEventActivityBins, eventActivityBins, nZVertexBins, zVertexBins);
227 }
228 
236 {
237  // Get eta bin for histos.
238 
239  Int_t etabin = -1;
240  eta = TMath::Abs(eta);
241  if (eta <= 0.4) etabin = 0;
242  else if (eta > 0.4 && eta < 0.8) etabin = 1;
243  else if (eta >= 0.8) etabin = 2;
244  return etabin;
245 }
246 
254 {
255  Int_t ptBin = -1;
256  if (pt < 0.5) ptBin = 0;
257  else if (pt < 1 ) ptBin = 1;
258  else if (pt < 2 ) ptBin = 2;
259  else if (pt < 3 ) ptBin = 3;
260  else if (pt < 5 ) ptBin = 4;
261  else if (pt < 8 ) ptBin = 5;
262  else if (pt < 20 ) ptBin = 6;
263 
264  return ptBin;
265 }
266 
274 {
275  // Get jet pt bin for histos.
276 
277  Int_t ptBin = -1;
278  if (pt >= 15 && pt < 20) ptBin = 0;
279  else if (pt >= 20 && pt < 25) ptBin = 1;
280  else if (pt >= 25 && pt < 30) ptBin = 2;
281  else if (pt >= 30 && pt < 60) ptBin = 3;
282  else if (pt >= 60) ptBin = 4;
283 
284  return ptBin;
285 }
286 
293 {
294  UInt_t eventTrigger = 0;
295  if (fIsEmbedded) {
296  eventTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
297  }
298  else {
299  auto embeddingHelper = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
300  if (embeddingHelper) {
301  auto aodHeader = dynamic_cast<AliVAODHeader *>(embeddingHelper->GetEventHeader());
302  if (aodHeader) {
303  eventTrigger = aodHeader->GetOfflineTrigger();
304  }
305  else {
306  AliErrorStream() << "Failed to retrieve requested AOD header from embedding helper\n";
307  }
308  }
309  else {
310  AliErrorStream() << "Failed to retrieve requested embedding helper\n";
311  }
312  }
313 
314  return eventTrigger;
315 }
316 
321 {
322  // NOTE: Clusters are never used directly in the task, so the container is neither created not retrieved
323  // Retrieve tracks
324  AliTrackContainer * tracks = static_cast<AliTrackContainer * >(GetParticleContainer("tracksForCorrelations"));
325  if (!tracks) {
326  AliError(Form("%s: Unable to retrieve tracks!", GetName()));
327  return kFALSE;
328  }
329 
330  // Retrieve jets
331  AliJetContainer * jets = GetJetContainer(0);
332  if (!jets) {
333  AliError(Form("%s: Unable to retrieve jets!", GetName()));
334  return kFALSE;
335  }
336 
337  // Used to calculate the angle betwene the jet and the hadron
338  TVector3 jetVector;
339  // Get z vertex
340  Double_t zVertex=fVertex[2];
341  // Flags
342  Bool_t biasedJet = kFALSE;
343  Bool_t leadJet = kFALSE;
344  // Relative angles and distances
345  Double_t deltaPhi = 0;
346  Double_t deltaEta = 0;
347  Double_t deltaR = 0;
348  // Event activity (centrality or multipilicity)
349  Double_t eventActivity = 0;
350  // Efficiency correction
351  Double_t efficiency = -999;
352  // Determining bins for histogram indices
353  Int_t jetPtBin = -1;
354  Int_t etaBin = -1;
355  // For comparison to the current jet
356  AliEmcalJet * leadingJet = jets->GetLeadingJet();
357  // For getting the proper properties of tracks
358  AliTLorentzVector track;
359 
360  // Determine the trigger for the current event
361  UInt_t eventTrigger = RetrieveTriggerMask();
362 
363  // Handle fast partition if selected
364  if ((eventTrigger & AliVEvent::kFastOnly) && fDisableFastPartition) {
365  AliDebugStream(4) << GetName() << ": Fast partition disabled\n";
366  if (fGeneralHistograms) {
367  fHistEventRejection->Fill("Fast Partition", 1);
368  }
369  return kFALSE;
370  }
371 
372  for (auto jet : jets->accepted()) {
373  // Selects only events that we are interested in (ie triggered)
374  if (!(eventTrigger & fTriggerType)) continue;
375  AliDebugStream(5) << "Jet accepted!\nJet: " << jet->toString().Data() << "\n";
376 
377  // Require the found jet to be matched
378  // This match should be between detector and particle level MC
379  if (fIsEmbedded) {
380  if (jet->MatchedJet()) {
381  AliDebugStream(4) << "Jet is matched!\nJet: " << jet->toString().Data() << "\n";
382  }
383  else {
384  AliDebugStream(5) << "Rejected jet because it was not matched to a external event jet.\n";
385  continue;
386  }
387  }
388 
389  // Jet properties
390  // Determine if we have the lead jet
391  leadJet = kFALSE;
392  if (jet == leadingJet) leadJet = kTRUE;
393 
394  // Determine if the jet is biased
395  biasedJet = BiasedJet(jet);
396 
397  // Calculate vector
398  jetVector.SetXYZ(jet->Px(), jet->Py(), jet->Pz());
399 
400  // Fill jet properties
401  FillHist(fHistJetPt[fCentBin], jet->Pt());
402  if (biasedJet == kTRUE) {
403  FillHist(fHistJetPtBias[fCentBin], jet->Pt());
404  }
405 
406  fHistJetEtaPhi->Fill(jet->Eta(), jet->Phi());
407 
408  if (jet->Pt() > 15) {
409 
410  AliDebugStream(4) << "Passed min jet pt cut of 15. \nJet: " << jet->toString().Data() << "\n";
411  for (auto trackIter : tracks->accepted_momentum()) {
412 
413  // Get proper track proeprties
414  track.Clear();
415  track = trackIter.first;
416 
417  // Determine relative angles and distances and set the respective variables
418  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
419 
420  // Determine bins for filling histograms
421  // jet Pt
422  jetPtBin = GetJetPtBin(jet->Pt());
423  if (jetPtBin < 0)
424  {
425  AliErrorStream() << "Jet Pt Bin negative: " << jet->Pt() << "\n";
426  continue;
427  }
428  // eta
429  etaBin = GetEtaBin(deltaEta);
430  if (etaBin < 0) {
431  AliErrorStream() << "Eta Bin negative: " << deltaEta << "\n";
432  continue;
433  }
434 
435  // Fill track properties
436  fHistTrackPt->Fill(track.Pt());
437 
438  if ( (jet->Pt() > 20.) && (jet->Pt() < 60.) ) {
439  fHistJHPsi->Fill(tracks->GetNTracks(), track.Pt(), track.Vect().Angle(jetVector) * TMath::RadToDeg() );
440  }
441 
442  fHistJetH[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
443  fHistJetHEtaPhi->Fill(deltaEta, deltaPhi);
444 
445  // Calculate single particle tracking efficiency for correlations
446  efficiency = EffCorrection(track.Eta(), track.Pt());
447  AliDebugStream(6) << GetName() << ": efficiency: " << efficiency << "\n";
448 
449  if (biasedJet == kTRUE) {
450  fHistJetHBias[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
451 
452  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
453  eventActivity = fCent;
454  }
455  else if (fBeamType == kpp) {
456  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
457  }
458 
459  if(fDoLessSparseAxes) { // check if we want all dimensions
460  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
461  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
462  } else {
463  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
464  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
465  }
466  }
467 
468  } //track loop
469  }//jet pt cut
470  }//jet loop
471 
472  //Prepare to do event mixing
473 
474  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
475  TObjArray* tracksClone = 0;
476 
477  if(fDoEventMixing == kTRUE){
478 
479  // event mixing
480 
481  // 1. First get an event pool corresponding in mult (cent) and
482  // zvertex to the current event. Once initialized, the pool
483  // should contain nMix (reduced) events. This routine does not
484  // pre-scan the chain. The first several events of every chain
485  // will be skipped until the needed pools are filled to the
486  // specified depth. If the pool categories are not too rare, this
487  // should not be a problem. If they are rare, you could lose
488  // statistics.
489 
490  // 2. Collect the whole pool's content of tracks into one TObjArray
491  // (bgTracks), which is effectively a single background super-event.
492 
493  // 3. The reduced and bgTracks arrays must both be passed into
494  // FillCorrelations(). Also nMix should be passed in, so a weight
495  // of 1./nMix can be applied.
496 
497  AliEventPool *pool = 0;
498  if (fBeamType == kAA || fBeamType == kpA) {//everything but pp
499  pool = fPoolMgr->GetEventPool(fCent, zVertex);
500  }
501  else if (fBeamType == kpp) {//pp only
502  pool = fPoolMgr->GetEventPool(static_cast<Double_t>(tracks->GetNTracks()), zVertex);
503  }
504 
505  if (!pool){
506  if (fBeamType == kAA || fBeamType == kpA) AliFatal(Form("No pool found for centrality = %f, zVertex = %f", fCent, zVertex));
507  else if (fBeamType == kpp) AliFatal(Form("No pool found for ntracks_pp = %d, zVertex = %f", tracks->GetNTracks(), zVertex));
508  return kTRUE;
509  }
510 
511  // The number of events in the pool
512  Int_t nMix = pool->GetCurrentNEvents();
513 
514  if(eventTrigger & fTriggerType) {
515  // check for a trigger jet
516  if (pool->IsReady() || pool->NTracksInPool() >= fMinNTracksMixedEvents || nMix >= fMinNEventsMixedEvents) {
517 
518  for (auto jet : jets->accepted()) {
519 
520  // Jet properties
521  // Determine if we have the lead jet
522  leadJet = kFALSE;
523  if (jet == leadingJet) { leadJet = kTRUE; }
524 
525  // Determine if the jet is biased
526  biasedJet = BiasedJet(jet);
527 
528  // Make sure event contains jet above our threshold (reduce stats of sparse)
529  if (jet->Pt() < 15) continue;
530 
531  // Fill for biased jet triggers only
532  if (biasedJet == kTRUE) {
533 
534  // Fill mixed-event histos here
535  for (Int_t jMix=0; jMix < nMix; jMix++) {
536  TObjArray* bgTracks = pool->GetEvent(jMix);
537 
538  for(Int_t ibg=0; ibg < bgTracks->GetEntries(); ibg++){
539 
540  AliBasicParticle *bgTrack = static_cast<AliBasicParticle*>(bgTracks->At(ibg));
541  if(!bgTrack)
542  {
543  AliError(Form("%s:Failed to retrieve tracks from mixed events", GetName()));
544  }
545 
546  // Fill into TLorentzVector for use with functions below
547  track.Clear();
548  track.SetPtEtaPhiE(bgTrack->Pt(), bgTrack->Eta(), bgTrack->Phi(), 0);
549 
550  // Calculate single particle tracking efficiency of mixed events for correlations
551  efficiency = EffCorrection(track.Eta(), track.Pt());
552 
553  // Phi is [-0.5*TMath::Pi(), 3*TMath::Pi()/2.]
554  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
555 
556  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
557  eventActivity = fCent;
558  }
559  else if (fBeamType == kpp) {
560  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
561  }
562 
563  if(fDoLessSparseAxes) { // check if we want all the axis filled
564  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
565  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
566  } else {
567  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
568  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
569  }
570  }
571  }
572  }
573  }
574  }
575  }
576 
577  if(eventTrigger & fMixingEventType) {
578  tracksClone = CloneAndReduceTrackList();
579 
580  //update pool if jet in event or not
581  pool->UpdatePool(tracksClone);
582  }
583 
584  } // end of event mixing
585 
586  return kTRUE;
587 }
588 
596 {
597  if ((jet->MaxTrackPt() > fTrackBias) || (jet->MaxClusterPt() > fClusterBias))
598  {
599  return kTRUE;
600  }
601  return kFALSE;
602 }
603 
613 void AliAnalysisTaskEmcalJetHCorrelations::GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector & particleOne, AliVParticle * particleTwo, Double_t & deltaEta, Double_t & deltaPhi, Double_t & deltaR)
614 {
615  // TODO: Understand order of arguments to DeltaPhi vs DeltaEta
616  // Returns deltaPhi in symmetric range so that we can calculate DeltaR.
617  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -1.0*TMath::Pi(), TMath::Pi());
618  deltaEta = particleOne.Eta() - particleTwo->Eta();
619  deltaR = TMath::Sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
620 
621  // Adjust to the normal range after the DeltaR caluclation
622  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -0.5*TMath::Pi(), 3*TMath::Pi()/2.);
623 }
624 
632 THnSparse* AliAnalysisTaskEmcalJetHCorrelations::NewTHnSparseF(const char* name, UInt_t entries)
633 {
634  Int_t count = 0;
635  UInt_t tmp = entries;
636  while(tmp!=0){
637  count++;
638  tmp = tmp &~ -tmp; // clear lowest bit
639  }
640 
641  TString hnTitle(name);
642  const Int_t dim = count;
643  Int_t nbins[dim];
644  Double_t xmin[dim];
645  Double_t xmax[dim];
646 
647  Int_t i=0;
648  Int_t c=0;
649  while(c<dim && i<32){
650  if(entries&(1<<i)){
651 
652  TString label("");
653  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
654  hnTitle += Form(";%s",label.Data());
655  c++;
656  }
657 
658  i++;
659  }
660  hnTitle += ";";
661 
662  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
663 }
664 
675 {
676  const Double_t pi = TMath::Pi();
677 
678  switch(iEntry){
679 
680  case 0:
681  label = "V0 centrality (%)";
682  nbins = 10;
683  xmin = 0.;
684  xmax = 100.;
685  break;
686 
687  case 1:
688  label = "corrected jet pt";
689  nbins = 20;
690  xmin = 0.;
691  xmax = 200.;
692  break;
693 
694  case 2:
695  if(fDoWiderTrackBin) {
696  label = "track pT";
697  nbins = 40;
698  xmin = 0.;
699  xmax = 10.;
700  } else {
701  label = "track pT";
702  nbins = 100;
703  xmin = 0.;
704  xmax = 10;
705  }
706  break;
707 
708  case 3:
709  label = "deltaEta";
710  nbins = 24;
711  xmin = -1.2;
712  xmax = 1.2;
713  break;
714 
715  case 4:
716  label = "deltaPhi";
717  nbins = 72;
718  xmin = -0.5*pi;
719  xmax = 1.5*pi;
720  break;
721 
722  case 5:
723  label = "leading jet";
724  nbins = 3;
725  xmin = -0.5;
726  xmax = 2.5;
727  break;
728 
729  case 6:
730  label = "trigger track";
731  nbins =10;
732  xmin = 0;
733  xmax = 50;
734  break;
735 
736  case 7:
737  label = "deltaR";
738  nbins = 10;
739  xmin = 0.;
740  xmax = 5.0;
741  break;
742 
743  case 8:
744  label = "leading track";
745  nbins = 13;
746  xmin = 0;
747  xmax = 50;
748  break;
749  }
750 }
751 
760 {
761  // clones a track list by using AliBasicTrack which uses much less memory (used for event mixing)
762  TObjArray* tracksClone = new TObjArray;
763  tracksClone->SetOwner(kTRUE);
764 
765  // Loop over all tracks
766  AliBasicParticle * clone = 0;
767  AliTrackContainer * tracks = GetTrackContainer("tracksForCorrelations");
768 
769  for (auto particle : tracks->accepted())
770  {
771  // Fill some QA information about the tracks
772  Int_t trackPtBin = GetTrackPtBin(particle->Pt());
773  if(trackPtBin > -1) fHistTrackEtaPhi[trackPtBin]->Fill(particle->Eta(),particle->Phi());
774 
775  // Create new particle
776  clone = new AliBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
777  // Set so that we can do comparisons using the IsEqual() function.
778  clone ->SetUniqueID(particle->GetUniqueID());
779 
780  tracksClone->Add(clone);
781  }
782 
783  return tracksClone;
784 }
785 
797  return EffCorrection(trackETA, trackPT, fBeamType);
798 }
799 
816 {
817  // default (current) parameters
818  // x-variable = track pt, y-variable = track eta
819  Double_t x = trackPT;
820  Double_t y = trackETA;
821  Double_t TRefficiency = -999;
822  Int_t runNUM = fCurrentRunNumber;
823  Int_t runSwitchGood = -999;
824  //Int_t centbin = -99;
825 
826  Double_t etaaxis = 0;
827  Double_t ptaxis = 0;
828 
829  Int_t effSwitch = fDoEffCorrection;
830 
831  if (beamType != AliAnalysisTaskEmcal::kpp) {
832  if(effSwitch == 1) {
833  // Semi-Good OROC C08 Runlists
834  if ((runNUM == 169975 || runNUM == 169981 || runNUM == 170038 || runNUM == 170040 || runNUM == 170083 || runNUM == 170084 || runNUM == 170085 || runNUM == 170088 || runNUM == 170089 || runNUM == 170091 || runNUM == 170152 || runNUM == 170155 || runNUM == 170159 || runNUM == 170163 || runNUM == 170193 || runNUM == 170195 || runNUM == 170203 || runNUM == 170204 || runNUM == 170228 || runNUM == 170230 || runNUM == 170268 || runNUM == 170269 || runNUM == 170270 || runNUM == 170306 || runNUM == 170308 || runNUM == 170309)) runSwitchGood = 0;
835 
836  // Good Runlists
837  if ((runNUM == 167902 || runNUM == 167903 || runNUM == 167915 || runNUM == 167920 || runNUM == 167987 || runNUM == 167988 || runNUM == 168066 || runNUM == 168068 || runNUM == 168069 || runNUM == 168076 || runNUM == 168104 || runNUM == 168107 || runNUM == 168108 || runNUM == 168115 || runNUM == 168212 || runNUM == 168310 || runNUM == 168311 || runNUM == 168322 || runNUM == 168325 || runNUM == 168341 || runNUM == 168342 || runNUM == 168361 || runNUM == 168362 || runNUM == 168458 || runNUM == 168460 || runNUM == 168461 || runNUM == 168464 || runNUM == 168467 || runNUM == 168511 || runNUM == 168512 || runNUM == 168777 || runNUM == 168826 || runNUM == 168984 || runNUM == 168988 || runNUM == 168992 || runNUM == 169035 || runNUM == 169091 || runNUM == 169094 || runNUM == 169138 || runNUM == 169143 || runNUM == 169144 || runNUM == 169145 || runNUM == 169148 || runNUM == 169156 || runNUM == 169160 || runNUM == 169167 || runNUM == 169238 || runNUM == 169411 || runNUM == 169415 || runNUM == 169417 || runNUM == 169835 || runNUM == 169837 || runNUM == 169838 || runNUM == 169846 || runNUM == 169855 || runNUM == 169858 || runNUM == 169859 || runNUM == 169923 || runNUM == 169956 || runNUM == 170027 || runNUM == 170036 || runNUM == 170081)) runSwitchGood = 1;
838 
839  // Determine which efficiency to use.
840  // This is just a way to map all possible values of the cent bin and runSwitchGood to a unique flag.
841  // 4 is the number of cent bins, and we want to index the effSwitch starting at 2.
842  if (runSwitchGood != -999) {
843  effSwitch = 2 + runSwitchGood*4 + fCentBin;
844  }
845  }
846 
847  // set up a switch for different parameter values...
848  switch(effSwitch) {
849  case 1 :
850  // first switch value - TRefficiency not used so = 1
851  // In this case, the run number isn't in any run list, so efficiency = 1
852  TRefficiency = 1.0;
853  break;
854 
855  case 2 :
856  // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
857  ptaxis = (x<2.9)*(p0_10SG[0]*exp(-pow(p0_10SG[1]/x,p0_10SG[2])) + p0_10SG[3]*x) + (x>=2.9)*(p0_10SG[4] + p0_10SG[5]*x + p0_10SG[6]*x*x);
858  etaaxis = (y<-0.07)*(p0_10SG[7]*exp(-pow(p0_10SG[8]/TMath::Abs(y+0.91),p0_10SG[9])) + p0_10SG[10]*y) + (y>=-0.07 && y<=0.4)*(p0_10SG[11] + p0_10SG[12]*y + p0_10SG[13]*y*y) + (y>0.4)*(p0_10SG[14]*exp(-pow(p0_10SG[15]/TMath::Abs(-y+0.91),p0_10SG[16])));
859  TRefficiency = ptaxis*etaaxis;
860  break;
861 
862  case 3 :
863  // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
864  ptaxis = (x<2.9)*(p10_30SG[0]*exp(-pow(p10_30SG[1]/x,p10_30SG[2])) + p10_30SG[3]*x) + (x>=2.9)*(p10_30SG[4] + p10_30SG[5]*x + p10_30SG[6]*x*x);
865  etaaxis = (y<-0.07)*(p10_30SG[7]*exp(-pow(p10_30SG[8]/TMath::Abs(y+0.91),p10_30SG[9])) + p10_30SG[10]*y) + (y>=-0.07 && y<=0.4)*(p10_30SG[11] + p10_30SG[12]*y + p10_30SG[13]*y*y) + (y>0.4)*(p10_30SG[14]*exp(-pow(p10_30SG[15]/TMath::Abs(-y+0.91),p10_30SG[16])));
866  TRefficiency = ptaxis*etaaxis;
867  break;
868 
869  case 4 :
870  // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
871  ptaxis = (x<2.9)*(p30_50SG[0]*exp(-pow(p30_50SG[1]/x,p30_50SG[2])) + p30_50SG[3]*x) + (x>=2.9)*(p30_50SG[4] + p30_50SG[5]*x + p30_50SG[6]*x*x);
872  etaaxis = (y<-0.07)*(p30_50SG[7]*exp(-pow(p30_50SG[8]/TMath::Abs(y+0.91),p30_50SG[9])) + p30_50SG[10]*y) + (y>=-0.07 && y<=0.4)*(p30_50SG[11] + p30_50SG[12]*y + p30_50SG[13]*y*y) + (y>0.4)*(p30_50SG[14]*exp(-pow(p30_50SG[15]/TMath::Abs(-y+0.91),p30_50SG[16])));
873  TRefficiency = ptaxis*etaaxis;
874  break;
875 
876  case 5 :
877  // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
878  ptaxis = (x<2.9)*(p50_90SG[0]*exp(-pow(p50_90SG[1]/x,p50_90SG[2])) + p50_90SG[3]*x) + (x>=2.9)*(p50_90SG[4] + p50_90SG[5]*x + p50_90SG[6]*x*x);
879  etaaxis = (y<-0.07)*(p50_90SG[7]*exp(-pow(p50_90SG[8]/TMath::Abs(y+0.91),p50_90SG[9])) + p50_90SG[10]*y) + (y>=-0.07 && y<=0.4)*(p50_90SG[11] + p50_90SG[12]*y + p50_90SG[13]*y*y) + (y>0.4)*(p50_90SG[14]*exp(-pow(p50_90SG[15]/TMath::Abs(-y+0.91),p50_90SG[16])));
880  TRefficiency = ptaxis*etaaxis;
881  break;
882 
883  case 6 :
884  // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
885  ptaxis = (x<2.9)*(p0_10G[0]*exp(-pow(p0_10G[1]/x,p0_10G[2])) + p0_10G[3]*x) + (x>=2.9)*(p0_10G[4] + p0_10G[5]*x + p0_10G[6]*x*x);
886  etaaxis = (y<0.0)*(p0_10G[7]*exp(-pow(p0_10G[8]/TMath::Abs(y+0.91),p0_10G[9])) + p0_10G[10]*y) + (y>=0.0 && y<=0.4)*(p0_10G[11] + p0_10G[12]*y + p0_10G[13]*y*y) + (y>0.4)*(p0_10G[14]*exp(-pow(p0_10G[15]/TMath::Abs(-y+0.91),p0_10G[16])));
887  TRefficiency = ptaxis*etaaxis;
888  break;
889 
890  case 7 :
891  // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
892  ptaxis = (x<2.9)*(p10_30G[0]*exp(-pow(p10_30G[1]/x,p10_30G[2])) + p10_30G[3]*x) + (x>=2.9)*(p10_30G[4] + p10_30G[5]*x + p10_30G[6]*x*x);
893  etaaxis = (y<0.0)*(p10_30G[7]*exp(-pow(p10_30G[8]/TMath::Abs(y+0.91),p10_30G[9])) + p10_30G[10]*y) + (y>=0.0 && y<=0.4)*(p10_30G[11] + p10_30G[12]*y + p10_30G[13]*y*y) + (y>0.4)*(p10_30G[14]*exp(-pow(p10_30G[15]/TMath::Abs(-y+0.91),p10_30G[16])));
894  TRefficiency = ptaxis*etaaxis;
895  break;
896 
897  case 8 :
898  // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
899  ptaxis = (x<2.9)*(p30_50G[0]*exp(-pow(p30_50G[1]/x,p30_50G[2])) + p30_50G[3]*x) + (x>=2.9)*(p30_50G[4] + p30_50G[5]*x + p30_50G[6]*x*x);
900  etaaxis = (y<0.0)*(p30_50G[7]*exp(-pow(p30_50G[8]/TMath::Abs(y+0.91),p30_50G[9])) + p30_50G[10]*y) + (y>=0.0 && y<=0.4)*(p30_50G[11] + p30_50G[12]*y + p30_50G[13]*y*y) + (y>0.4)*(p30_50G[14]*exp(-pow(p30_50G[15]/TMath::Abs(-y+0.91),p30_50G[16])));
901  TRefficiency = ptaxis*etaaxis;
902  break;
903 
904  case 9 :
905  // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
906  ptaxis = (x<2.9)*(p50_90G[0]*exp(-pow(p50_90G[1]/x,p50_90G[2])) + p50_90G[3]*x) + (x>=2.9)*(p50_90G[4] + p50_90G[5]*x + p50_90G[6]*x*x);
907  etaaxis = (y<0.0)*(p50_90G[7]*exp(-pow(p50_90G[8]/TMath::Abs(y+0.91),p50_90G[9])) + p50_90G[10]*y) + (y>=0.0 && y<=0.4)*(p50_90G[11] + p50_90G[12]*y + p50_90G[13]*y*y) + (y>0.4)*(p50_90G[14]*exp(-pow(p50_90G[15]/TMath::Abs(-y+0.91),p50_90G[16])));
908  TRefficiency = ptaxis*etaaxis;
909  break;
910 
911  default :
912  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
913  // ie. The efficiency correction is disabled.
914  TRefficiency = 1.0;
915  }
916  }
917  else {
918  // Track efficiency for pp
919  // Calculated using LHC12f1a. See analysis note for more details!
920 
921  if (fDoEffCorrection != 0) {
922  // If the trackPt > 6 GeV, then all we need is this coefficient
923  Double_t coefficient = 0.898052; // p6
924  if (trackPT < 6) {
925  coefficient = (1 + -0.442232 * trackPT // p0
926  + 0.501831 * std::pow(trackPT, 2) // p1
927  + -0.252024 * std::pow(trackPT, 3) // p2
928  + 0.062964 * std::pow(trackPT, 4) // p3
929  + -0.007681 * std::pow(trackPT, 5) // p4
930  + 0.000365 * std::pow(trackPT, 6)); // p5
931  }
932 
933  // Calculate track eff
934  TRefficiency = coefficient * (1 + 0.402825 * std::abs(trackETA) // p7
935  + -2.213152 * std::pow(trackETA, 2) // p8
936  + 4.311098 * std::abs(std::pow(trackETA, 3)) // p9
937  + -2.778200 * std::pow(trackETA, 4)); // p10
938  }
939  else {
940  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
941  TRefficiency = 1;
942  }
943  }
944 
945  return TRefficiency;
946 }
947 
957 void AliAnalysisTaskEmcalJetHCorrelations::FillHist(TH1 * hist, Double_t fillValue, Double_t weight, Bool_t noCorrection)
958 {
959  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
960  {
961  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
962  hist->Fill(fillValue, weight);
963  }
964  else
965  {
966  // Determine where to get the values in the correction hist
967  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(fillValue);
968 
969  std::vector <Double_t> yBinsContent;
970  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), fillValue));
971  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
972  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
973 
974  // Loop over all possible bins to contribute.
975  // If content is 0 then calling Fill won't make a difference
976  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
977  {
978  // Don't bother trying to fill in the weight is 0
979  if (yBinsContent.at(index-1) > 0) {
980  // Determine the value to fill based on the center of the bins.
981  // This in principle allows the binning between the correction and hist to be different
982  Double_t fillLocation = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
983  AliDebug(4, TString::Format("fillLocation: %f, weight: %f", fillLocation, yBinsContent.at(index-1)));
984  // minus 1 since loop starts at 1
985  hist->Fill(fillLocation, weight*yBinsContent.at(index-1));
986  }
987  }
988 
989  //TEMP
990  //hist->Draw();
991  //END TEMP
992  }
993 }
994 
1005 void AliAnalysisTaskEmcalJetHCorrelations::FillHist(THnSparse * hist, Double_t *fillValue, Double_t weight, Bool_t noCorrection)
1006 {
1007  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
1008  {
1009  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
1010  hist->Fill(fillValue, weight);
1011  }
1012  else
1013  {
1014  // Jet pt is always located in the second position
1015  Double_t jetPt = fillValue[1];
1016 
1017  // Determine where to get the values in the correction hist
1018  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(jetPt);
1019 
1020  std::vector <Double_t> yBinsContent;
1021  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), jetPt));
1022  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
1023  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
1024 
1025  // Loop over all possible bins to contribute.
1026  // If content is 0 then calling Fill won't make a difference
1027  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
1028  {
1029  // Don't bother trying to fill in the weight is 0
1030  if (yBinsContent.at(index-1) > 0) {
1031  // Determine the value to fill based on the center of the bins.
1032  // This in principle allows the binning between the correction and hist to be different
1033  fillValue[1] = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
1034  AliDebug(4,TString::Format("fillValue[1]: %f, weight: %f", fillValue[1], yBinsContent.at(index-1)));
1035  // minus 1 since loop starts at 1
1036  hist->Fill(fillValue, weight*yBinsContent.at(index-1));
1037  }
1038  }
1039  }
1040 }
1041 
1051 void AliAnalysisTaskEmcalJetHCorrelations::AccessSetOfYBinValues(TH2D * hist, Int_t xBin, std::vector <Double_t> & yBinsContent, Double_t scaleFactor)
1052 {
1053  for (Int_t index = 1; index <= hist->GetYaxis()->GetNbins(); index++)
1054  {
1055  //yBinsContent[index-1] = hist->GetBinContent(hist->GetBin(xBin,index));
1056  yBinsContent.push_back(hist->GetBinContent(hist->GetBin(xBin,index)));
1057 
1058  if (scaleFactor >= 0)
1059  {
1060  // -1 since index starts at 1
1061  hist->SetBinContent(hist->GetBin(xBin,index), yBinsContent.at(index-1)/scaleFactor);
1062  }
1063  }
1064 }
1065 
1083 {
1084  // Initialize grid connection if necessary
1085  if (filename.Contains("alien://") && !gGrid) {
1086  TGrid::Connect("alien://");
1087  }
1088 
1089  // Setup hist name if a track or cluster bias was defined.
1090  // NOTE: This can always be disabled by setting kDisableBias.
1091  // We arbitrarily add 0.1 to test since the values are doubles and cannot be
1092  // tested directly for equality. If we are still less than disable bins, then
1093  // it has been set and we should format it.
1094  // NOTE: To ensure we can disable, we don't just take the member values!
1095  // NOTE: The histBaseName will be attempted if the formatted name cannot be found.
1096  TString histBaseName = histName;
1098  histName = TString::Format("%s_Track%.2f", histName.Data(), trackBias);
1099  }
1100  if (clusterBias + 0.1 < AliAnalysisTaskEmcalJetHCorrelations::kDisableBias) {
1101  histName = TString::Format("%s_Clus%.2f", histName.Data(), clusterBias);
1102  }
1103 
1104  // Open file containing the correction
1105  TFile * jesCorrectionFile = TFile::Open(filename);
1106  if (!jesCorrectionFile || jesCorrectionFile->IsZombie()) {
1107  AliError(TString::Format("%s: Could not open JES correction file %s", GetName(), filename.Data()));
1108  return kFALSE;
1109  }
1110 
1111  // Retrieve the histogram containing the correction and safely add it to the task.
1112  TH2D * JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histName.Data()));
1113  if (JESCorrectionHist) {
1114  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histName.Data(), filename.Data()));
1115  }
1116  else {
1117  AliError(TString::Format("%s: JES correction hist name \"%s\" not found in file %s.", GetName(), histName.Data(), filename.Data()));
1118 
1119  // Attempt the base name instead of the formatted hist name
1120  JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histBaseName.Data()));
1121  if (JESCorrectionHist) {
1122  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histBaseName.Data(), filename.Data()));
1123  histName = histBaseName;
1124  }
1125  else
1126  {
1127  AliError(TString::Format("%s: JES correction with base hist name %s not found in file %s.", GetName(), histBaseName.Data(), filename.Data()));
1128  return kFALSE;
1129  }
1130  }
1131 
1132  // Clone to ensure that the hist is available
1133  TH2D * tempHist = static_cast<TH2D *>(JESCorrectionHist->Clone());
1134  tempHist->SetDirectory(0);
1135  SetJESCorrectionHist(tempHist);
1136 
1137  // Close file
1138  jesCorrectionFile->Close();
1139 
1140  // Append to task name for clarity
1141  // Unfortunately, this doesn't change the name of the output list (it would need to be
1142  // changed in the AnalysisManager output container), so the suffix is still important
1143  // if this correction is manually configured!
1144  TString tempName = GetName();
1145  TString tag = "_JESCorr";
1146  // Append the tag if it isn't already included
1147  if (tempName.Index(tag) == -1) {
1148  // Insert before the suffix
1149  Ssiz_t suffixLocation = tempName.Last('_');
1150  tempName.Insert(suffixLocation, tag.Data());
1151 
1152  // Set the new name
1153  AliDebug(3, TString::Format("%s: Setting task name to %s", GetName(), tempName.Data()));
1154  SetName(tempName.Data());
1155  }
1156 
1157  // Successful
1158  return kTRUE;
1159 }
1160 
1166  const char *nTracks,
1167  const char *nCaloClusters,
1168  // Jet options
1169  const Double_t trackBias,
1170  const Double_t clusterBias,
1171  // Mixed event options
1172  const Int_t nTracksMixedEvent, // Additionally acts as a switch for enabling mixed events
1173  const Int_t minNTracksMixedEvent,
1174  const Int_t minNEventsMixedEvent,
1175  const UInt_t nCentBinsMixedEvent,
1176  // Triggers
1177  UInt_t trigEvent,
1178  UInt_t mixEvent,
1179  // Options
1180  const Bool_t lessSparseAxes,
1181  const Bool_t widerTrackBin,
1182  // Corrections
1183  const Int_t doEffCorrSW,
1184  const Bool_t JESCorrection,
1185  const char * JESCorrectionFilename,
1186  const char * JESCorrectionHistName,
1187  const char *suffix
1188  )
1189 {
1190  // Get the pointer to the existing analysis manager via the static access method.
1191  //==============================================================================
1192  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1193  if (!mgr)
1194  {
1195  AliErrorClass("No analysis manager to connect to.");
1196  return nullptr;
1197  }
1198 
1199  //-------------------------------------------------------
1200  // Init the task and do settings
1201  //-------------------------------------------------------
1202 
1203  // Determine cluster and track names
1204  TString trackName(nTracks);
1205  TString clusName(nCaloClusters);
1206 
1207  if (trackName == "usedefault") {
1209  }
1210 
1211  if (clusName == "usedefault") {
1213  }
1214 
1215  TString name("AliAnalysisTaskJetH");
1216  if (!trackName.IsNull()) {
1217  name += TString::Format("_%s", trackName.Data());
1218  }
1219  if (!clusName.IsNull()) {
1220  name += TString::Format("_%s", clusName.Data());
1221  }
1222  if (strcmp(suffix, "") != 0) {
1223  name += TString::Format("_%s", suffix);
1224  }
1225 
1227  // Set jet bias
1228  correlationTask->SetTrackBias(trackBias);
1229  correlationTask->SetClusterBias(clusterBias);
1230  // Mixed events
1231  correlationTask->SetEventMixing(static_cast<Bool_t>(nTracksMixedEvent));
1232  correlationTask->SetNumberOfMixingTracks(nTracksMixedEvent);
1233  correlationTask->SetMinNTracksForMixedEvents(minNTracksMixedEvent);
1234  correlationTask->SetMinNEventsForMixedEvents(minNEventsMixedEvent);
1235  correlationTask->SetNCentBinsMixedEvent(nCentBinsMixedEvent);
1236  // Triggers
1237  correlationTask->SetTriggerType(trigEvent);
1238  correlationTask->SetMixedEventTriggerType(mixEvent);
1239  // Options
1240  correlationTask->SetNCentBins(5);
1241  correlationTask->SetVzRange(-10,10);
1242  correlationTask->SetDoLessSparseAxes(lessSparseAxes);
1243  correlationTask->SetDoWiderTrackBin(widerTrackBin);
1244  // Corrections
1245  correlationTask->SetDoEffCorr(doEffCorrSW);
1246  if (JESCorrection == kTRUE)
1247  {
1248  Bool_t result = correlationTask->RetrieveAndInitializeJESCorrectionHist(JESCorrectionFilename, JESCorrectionHistName, correlationTask->GetTrackBias(), correlationTask->GetClusterBias());
1249  if (!result) {
1250  AliErrorClass("Failed to successfully retrieve and initialize the JES correction! Task initialization continuing without JES correction (can be set manually later).");
1251  }
1252  }
1253 
1254  //-------------------------------------------------------
1255  // Final settings, pass to manager and set the containers
1256  //-------------------------------------------------------
1257 
1258  mgr->AddTask(correlationTask);
1259 
1260  // Create containers for input/output
1261  mgr->ConnectInput (correlationTask, 0, mgr->GetCommonInputContainer() );
1262  AliAnalysisDataContainer * cojeth = mgr->CreateContainer(correlationTask->GetName(),
1263  TList::Class(),
1264  AliAnalysisManager::kOutputContainer,
1265  Form("%s", AliAnalysisManager::GetCommonFileName()));
1266  mgr->ConnectOutput(correlationTask, 1, cojeth);
1267 
1268  return correlationTask;
1269 }
1270 
1277  std::string clusName,
1278  const double jetConstituentPtCut,
1279  const double trackEta,
1280  const double jetRadius)
1281 {
1282  bool returnValue = false;
1283  AliInfoStream() << "Configuring Jet-H Correlations task for a standard analysis.\n";
1284 
1285  // Add Containers
1286  // Clusters
1287  if (clusName == "usedefault") {
1289  }
1290  // For jet finding
1291  AliClusterContainer * clustersForJets = new AliClusterContainer(clusName.c_str());
1292  clustersForJets->SetName("clustersForJets");
1293  clustersForJets->SetMinE(jetConstituentPtCut);
1294 
1295  // Tracks
1296  // For jet finding
1297  if (trackName == "usedefault") {
1299  }
1300  AliParticleContainer * particlesForJets = CreateParticleOrTrackContainer(trackName.c_str());
1301  particlesForJets->SetName("particlesForJets");
1302  particlesForJets->SetMinPt(jetConstituentPtCut);
1303  particlesForJets->SetEtaLimits(-1.0*trackEta, trackEta);
1304  // Don't need to adopt the container - we'll just use it to find the right jet collection
1305  // For correlations
1306  AliParticleContainer * particlesForCorrelations = CreateParticleOrTrackContainer(trackName.c_str());
1307  if (particlesForCorrelations)
1308  {
1309  particlesForCorrelations->SetName("tracksForCorrelations");
1310  particlesForCorrelations->SetMinPt(0.15);
1311  particlesForCorrelations->SetEtaLimits(-1.0*trackEta, trackEta);
1312  // Adopt the container
1313  this->AdoptParticleContainer(particlesForCorrelations);
1314  }
1315  else {
1316  AliWarningStream() << "No particle container was successfully created!\n";
1317  }
1318 
1319  // Jets
1323  jetRadius,
1325  particlesForJets,
1326  clustersForJets);
1327  // 0.6 * jet area
1328  jetContainer->SetJetAreaCut(jetRadius * jetRadius * TMath::Pi() * 0.6);
1329  jetContainer->SetMaxTrackPt(100);
1330  jetContainer->SetJetPtCut(0.1);
1331 
1332  // Successfully configured
1333  returnValue = true;
1334 
1335  return returnValue;
1336 }
1337 
1344  std::string clusName,
1345  const double jetConstituentPtCut,
1346  const double trackEta,
1347  const double jetRadius,
1348  const std::string & jetTag)
1349 {
1350  bool returnValue = false;
1351  AliInfoStream() << "Configuring Jet-H Correlations task for an embedding analysis.\n";
1352 
1353  // Set the task to know it that is embedded
1354  this->SetIsEmbedded(true);
1355 
1356  // Add Containers
1357  // Clusters
1358  if (clusName == "usedefault") {
1360  }
1361  // For jet finding
1362  AliClusterContainer * clustersForJets = new AliClusterContainer(clusName.c_str());
1363  clustersForJets->SetName("clustersForJets");
1364  clustersForJets->SetMinE(jetConstituentPtCut);
1365  // We need the combined clusters, which should be available in the internal event.
1366  // However, we don't need to adopt the container - we'll just use it to find the right jet collection
1367  // For correlations
1368  /*AliClusterContainer * clustersforCorrelations = new AliClusterContainer("usedefault");
1369  clustersForCorrelations->SetName("clustersForCorrelations");
1370  clustersForCorrelations->SetMinE(0.30);
1371  clustersForCorrelations->SetIsEmbedding(true);
1372  this->AdoptClusterContainer(clustersForCorrelations);*/
1373 
1374  // Tracks
1375  // For jet finding
1376  if (trackName == "usedefault") {
1378  }
1379  AliParticleContainer * particlesForJets = CreateParticleOrTrackContainer(trackName.c_str());
1380  particlesForJets->SetName("particlesForJets");
1381  particlesForJets->SetMinPt(jetConstituentPtCut);
1382  particlesForJets->SetEtaLimits(-1.0*trackEta, trackEta);
1383  // Don't need to adopt the container - we'll just use it to find the right jet collection
1384  // For correlations
1385  AliParticleContainer * particlesForCorrelations = CreateParticleOrTrackContainer(trackName.c_str());
1386  // Ensure that we don't operate on a null pointer
1387  if (particlesForCorrelations)
1388  {
1389  particlesForCorrelations->SetName("tracksForCorrelations");
1390  particlesForCorrelations->SetMinPt(0.15);
1391  particlesForCorrelations->SetEtaLimits(-1.0*trackEta, trackEta);
1392  particlesForCorrelations->SetIsEmbedding(true);
1393  // Adopt the container
1394  this->AdoptParticleContainer(particlesForCorrelations);
1395  }
1396  else {
1397  AliWarningStream() << "No particle container was successfully created!\n";
1398  }
1399 
1400  // Jets
1401  // The tag "hybridJets" is defined in the jet finder
1405  jetRadius,
1407  particlesForJets,
1408  clustersForJets,
1409  jetTag);
1410  // 0.6 * jet area
1411  jetContainer->SetJetAreaCut(jetRadius * jetRadius * TMath::Pi() * 0.6);
1412  jetContainer->SetMaxTrackPt(100);
1413  jetContainer->SetJetPtCut(0.1);
1414 
1415  // Successfully configured
1416  returnValue = true;
1417 
1418  return returnValue;
1419 }
1420 
1429 {
1430  AliParticleContainer * partCont = 0;
1432  AliTrackContainer * trackCont = new AliTrackContainer(collectionName.c_str());
1433  partCont = trackCont;
1434  }
1435  else if (collectionName != "") {
1436  partCont = new AliParticleContainer(collectionName.c_str());
1437  }
1438 
1439  return partCont;
1440 }
static AliAnalysisTaskEmcalJetHCorrelations * AddTaskEmcalJetHCorrelations(const char *nTracks="usedefault", const char *nCaloClusters="usedefault", const Double_t trackBias=5, const Double_t clusterBias=5, const Int_t nTracksMixedEvent=0, const Int_t minNTracksMixedEvent=5000, const Int_t minNEventsMixedEvent=5, const UInt_t nCentBinsMixedEvent=10, UInt_t trigEvent=AliVEvent::kAny, UInt_t mixEvent=AliVEvent::kAny, const Bool_t lessSparseAxes=kFALSE, const Bool_t widerTrackBin=kFALSE, const Int_t doEffCorrSW=0, const Bool_t JESCorrection=kFALSE, const char *JESCorrectionFilename="alien:///alice/cern.ch/user/r/rehlersi/JESCorrection.root", const char *JESCorrectionHistName="JESCorrection", const char *suffix="biased")
virtual void GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
Bool_t RetrieveAndInitializeJESCorrectionHist(TString filename, TString histName, Double_t trackBias=AliAnalysisTaskEmcalJetHCorrelations::kDisableBias, Double_t clusterBias=AliAnalysisTaskEmcalJetHCorrelations::kDisableBias)
void FillHist(TH1 *hist, Double_t fillValue, Double_t weight=1.0, Bool_t noCorrection=kFALSE)
const char * filename
Definition: TestFCM.C:1
virtual void SetTrackBias(Double_t b)
Require a track with pt > b in jet.
TH2 * fHistJetHBias[6][5][3]
! Jet-hadron correlations of jets which meet the constituent bias criteria (the arrays correspond to ...
double Double_t
Definition: External.C:58
Definition: External.C:260
static Double_t p50_90SG[17]
50-90% centrality semi-good runs
AliParticleContainer * CreateParticleOrTrackContainer(std::string const &collectionName) const
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: External.C:236
Int_t GetNTracks() const
Bool_t fDisableFastPartition
True if task should be disabled for the fast partition, where the EMCal is not included.
void AccessSetOfYBinValues(TH2D *hist, Int_t xBin, std::vector< Double_t > &yBinsContent, Double_t scaleFactor=-1.0)
AliJetContainer * GetJetContainer(Int_t i=0) const
void AdoptParticleContainer(AliParticleContainer *cont)
const AliTrackIterableContainer accepted() const
static Double_t p0_10SG[17]
0-10% centrality semi-good runs
Container with name, TClonesArray and cuts for particles.
Declaration of class AliTLorentzVector.
static Double_t p10_30G[17]
10-30% centrality good runs
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
virtual void SetMixedEventTriggerType(UInt_t me)
Set the mixed event trigger selection.
TH2D * fJESCorrectionHist
Histogram containing the jet energy scale correction.
TCanvas * c
Definition: TestFitELoss.C:172
TH2 * fHistJetH[6][5][3]
! Jet-hadron correlations (the arrays correspond to centrality, jet pt bins, and eta bins) ...
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Int_t fCentBin
!event centrality bin
void SetVzRange(Double_t min, Double_t max)
virtual Double_t Pt() const
TH1 * fHistJetPt[6]
! Jet pt spectrum (the array corresponds to centrality bins)
virtual Double_t Eta() const
TH2 * fHistJetHEtaPhi
! Eta-phi distribution of jets which are in jet-hadron correlations
UInt_t fTriggerType
Event selection for jets (ie triggered events).
Container for particles within the EMCAL framework.
Bool_t fIsEmbedded
trigger, embedded signal
BeamType
Switch for the beam type.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH2 * fHistTrackEtaPhi[7]
! Track eta-phi distribution (the array corresponds to track pt)
AliEmcalJet * GetLeadingJet(const char *opt="")
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
unsigned int UInt_t
Definition: External.C:33
BeamType fForceBeamType
forced beam type
Definition: External.C:228
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:155
BeamType fBeamType
!event beam type
Double_t fCent
!event centrality
virtual void SetTriggerType(UInt_t te)
Set the trigger event trigger selection.
Bool_t fDoWiderTrackBin
True if the track pt bins in the THnSparse should be wider.
THnSparse * fhnMixedEvents
! Mixed events THnSparse
UInt_t fMixingEventType
Event selection for mixed events.
static Double_t p50_90G[17]
50-90% centrality good runs
! Number of elements in mixed event multiplicity binned arrays
Int_t fMinNTracksMixedEvents
threshold to use event pool # tracks
virtual void SetNCentBins(Int_t n)
! Arbitrarily large value which can be used to disable the constituent bias. Can be used for either t...
bool ConfigureForStandardAnalysis(std::string trackName="usedefault", std::string clusName="usedefault", const double jetConstituentPtCut=3, const double trackEta=0.8, const double jetRadius=0.2)
Bool_t fNoMixedEventJESCorrection
True if the jet energy scale correction should be applied to mixed event histograms.
AliEventPoolManager * fPoolMgr
! Event pool manager
Int_t fDoEffCorrection
Control the efficiency correction. See EffCorrection() for meaning of values.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
virtual THnSparse * NewTHnSparseF(const char *name, UInt_t entries)
AliEmcalList * fOutput
!output list
static Double_t p0_10G[17]
0-10% centrality good runs
bool ConfigureForEmbeddingAnalysis(std::string trackName="usedefault", std::string clusName="caloClustersCombined", const double jetConstituentPtCut=3, const double trackEta=0.8, const double jetRadius=0.2, const std::string &jetTag="hybridJets")
Double_t fVertex[3]
!event vertex
AliTrackContainer * GetTrackContainer(Int_t i=0) const
TH1 * fHistEventRejection
!book keep reasons for rejecting event
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Double_t EffCorrection(Double_t trkETA, Double_t trkPT, AliAnalysisTaskEmcal::BeamType beamType) const
Double_t MaxClusterPt() const
Definition: AliEmcalJet.h:154
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
static std::string DetermineUseDefaultName(InputObject_t objType)
const AliTrackIterableMomentumContainer accepted_momentum() const
static Double_t p30_50SG[17]
30-50% centrality semi-good runs
void GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector &particleOne, AliVParticle *particleTwo, Double_t &deltaEta, Double_t &deltaPhi, Double_t &deltaR)
void UserCreateOutputObjects()
Main initialization function on the worker.
const Double_t pi
const Int_t nbins
Bool_t fDoLessSparseAxes
True if there should be fewer THnSparse axes.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Jet-hadron correlations analysis task for central Pb-Pb and pp.
UInt_t fNCentBinsMixedEvent
N cent bins for the event mixing pool.
void SetMaxTrackPt(Float_t b)
TH1 * fHistJetPtBias[6]
! Jet pt spectrum of jets which meet the constituent bias criteria (the array corresponds to centrali...
virtual void SetClusterBias(Double_t b)
Require a cluster with pt > b in jet.
virtual Double_t Phi() const
Container structure for EMCAL clusters.
Int_t fMinNEventsMixedEvents
threshold to use event pool # events
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Int_t fNMixingTracks
size of track buffer for event mixing
static Double_t p30_50G[17]
30-50% centrality good runs
Container for jet within the EMCAL jet framework.
static Double_t p10_30SG[17]
10-30% centrality semi-good runs
Definition: External.C:196
void SetJetAreaCut(Float_t cut)
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()