AliPhysics  5364b50 (5364b50)
AliAnalysisTaskEmcalJetHCorrelations.cxx
Go to the documentation of this file.
1 //Measure Jet-hadron correlations
3 //Does event Mixing using AliEventPoolManager
5 
7 
8 #include <bitset>
9 
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TH3F.h>
13 #include <THnSparse.h>
14 #include <TVector3.h>
15 #include <TFile.h>
16 #include <TGrid.h>
17 
18 #include <AliAnalysisManager.h>
19 #include <AliInputEventHandler.h>
20 #include <AliEventPoolManager.h>
21 #include <AliLog.h>
22 #include <AliVAODHeader.h>
23 #include <AliVTrack.h>
24 
25 #include "AliEmcalJet.h"
26 #include "AliTLorentzVector.h"
27 #include "AliBasicParticle.h"
28 #include "AliEmcalContainerUtils.h"
29 #include "AliClusterContainer.h"
30 #include "AliTrackContainer.h"
31 #include "AliJetContainer.h"
33 
37 
38 // 0-10% centrality: Semi-Good Runs
39 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};
40 // 10-30% centrality: Semi-Good Runs
41 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};
42 // 30-50% centrality: Semi-Good Runs
43 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};
44 // 50-90% centrality: Semi-Good Runs
45 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};
46 
47 // 0-10% centrality: Good Runs
48 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};
49 // 10-30% centrality: Good Runs
50 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};
51 // 30-50% centrality: Good Runs
52 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};
53 // 50-90% centrality: Good Runs
54 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};
55 
60  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetHCorrelations", kFALSE),
61  fTrackBias(5),
62  fClusterBias(5),
63  fDoEventMixing(kFALSE),
64  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
65  fPoolMgr(nullptr),
66  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
67  fDisableFastPartition(kFALSE),
68  fDoEffCorrection(0),
69  fNoMixedEventJESCorrection(kFALSE),
70  fJESCorrectionHist(nullptr),
71  fDoLessSparseAxes(kFALSE), fDoWiderTrackBin(kFALSE),
72  fRequireMatchedJetWhenEmbedding(kTRUE),
73  fHistTrackPt(nullptr),
74  fHistJetEtaPhi(nullptr),
75  fHistJetHEtaPhi(nullptr),
76  fhnMixedEvents(nullptr),
77  fhnJH(nullptr),
78  fhnTrigger(nullptr)
79 {
80  // Default Constructor
82 }
83 
88  AliAnalysisTaskEmcalJet(name, kTRUE),
89  fTrackBias(5),
90  fClusterBias(5),
91  fDoEventMixing(kFALSE),
94  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
95  fDisableFastPartition(kFALSE),
99  fDoLessSparseAxes(kFALSE), fDoWiderTrackBin(kFALSE),
105  fhnJH(nullptr),
107 {
108  // Constructor
110  // Ensure that additional general histograms are created
112 }
113 
118 {
119  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; trackPtBin++){
120  fHistTrackEtaPhi[trackPtBin] = nullptr;
121  }
122  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
123  fHistJetPt[centralityBin] = nullptr;
124  fHistJetPtBias[centralityBin] = nullptr;
125  }
126 }
127 
132  // Called once
134 
135  // Create histograms
136  fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
137  fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
138  fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
139 
140  fOutput->Add(fHistTrackPt);
141  fOutput->Add(fHistJetEtaPhi);
142  fOutput->Add(fHistJetHEtaPhi);
143 
144  TString name;
145  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; ++trackPtBin){
146  name = Form("fHistTrackEtaPhi_%i", trackPtBin);
147  fHistTrackEtaPhi[trackPtBin] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
148  fOutput->Add(fHistTrackEtaPhi[trackPtBin]);
149  }
150 
151  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
152  name = Form("fHistJetPt_%i",centralityBin);
153  fHistJetPt[centralityBin] = new TH1F(name,name,200,0,200);
154  fOutput->Add(fHistJetPt[centralityBin]);
155 
156  name = Form("fHistJetPtBias_%i",centralityBin);
157  fHistJetPtBias[centralityBin] = new TH1F(name,name,200,0,200);
158  fOutput->Add(fHistJetPtBias[centralityBin]);
159  }
160 
161  UInt_t cifras = 0; // bit coded, see GetDimParams() below
162  if(fDoLessSparseAxes) {
163  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<9;
164  } else {
165  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7 | 1<<9;
166  }
167  fhnJH = NewTHnSparseF("fhnJH", cifras);
168  fhnJH->Sumw2();
169  fOutput->Add(fhnJH);
170 
171  if(fDoEventMixing){
172  // The event plane angle does not need to be included because the semi-central determined that the EP angle didn't change
173  // significantly for any of the EP orientations. However, it will be included so this can be demonstrated for the central
174  // analysis if so desired.
175  if(fDoLessSparseAxes) {
176  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<9;
177  } else {
178  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7 | 1<<9;
179  }
180  fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
181  fhnMixedEvents->Sumw2();
182  fOutput->Add(fhnMixedEvents);
183  }
184 
185  // Trigger THnSparse
186  cifras = 1<<0 | 1<<1 | 1<<9;
187  fhnTrigger = NewTHnSparseF("fhnTrigger", cifras);
188  fhnTrigger->Sumw2();
189  fOutput->Add(fhnTrigger);
190 
191  PostData(1, fOutput);
192 
193  // Event Mixing
194  Int_t poolSize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
195  // ZVertex
196  Int_t nZVertexBins = 10;
197  Double_t* zVertexBins = GenerateFixedBinArray(nZVertexBins, -10, 10);
198  // Event activity (centrality of multiplicity)
199  Int_t nEventActivityBins = 8;
200  Double_t* eventActivityBins = 0;
201  // +1 to accomodate the fact that we define bins rather than array entries.
202  Double_t multiplicityBins[kMixedEventMulitplictyBins+1] = {0., 4., 9., 15., 25., 35., 55., 100., 500.};
203 
204  // Cannot use GetBeamType() since it is not available until UserExec()
205  if (fForceBeamType != AliAnalysisTaskEmcal::kpp ) { //all besides pp
206  // Event Activity is centrality in AA, pA
207  nEventActivityBins = fNCentBinsMixedEvent;
208  eventActivityBins = GenerateFixedBinArray(nEventActivityBins, 0, 100);
209  }
210  else if (fForceBeamType == AliAnalysisTaskEmcal::kpp) { //for pp only
211  // Event Activity is multiplicity in pp
212  eventActivityBins = multiplicityBins;
213  }
214 
215  fPoolMgr = new AliEventPoolManager(poolSize, fNMixingTracks, nEventActivityBins, eventActivityBins, nZVertexBins, zVertexBins);
216 }
217 
225 {
226  Int_t ptBin = -1;
227  if (pt < 0.5) ptBin = 0;
228  else if (pt < 1 ) ptBin = 1;
229  else if (pt < 2 ) ptBin = 2;
230  else if (pt < 3 ) ptBin = 3;
231  else if (pt < 5 ) ptBin = 4;
232  else if (pt < 8 ) ptBin = 5;
233  else if (pt < 20 ) ptBin = 6;
234 
235  return ptBin;
236 }
237 
244 {
245  UInt_t eventTrigger = 0;
246  if (fIsEmbedded) {
247  auto embeddingHelper = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
248  if (embeddingHelper) {
249  auto aodHeader = dynamic_cast<AliVAODHeader *>(embeddingHelper->GetEventHeader());
250  if (aodHeader) {
251  eventTrigger = aodHeader->GetOfflineTrigger();
252  }
253  else {
254  AliErrorStream() << "Failed to retrieve requested AOD header from embedding helper\n";
255  }
256  }
257  else {
258  AliErrorStream() << "Failed to retrieve requested embedding helper\n";
259  }
260  }
261  else {
262  eventTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
263  }
264 
265  return eventTrigger;
266 }
267 
272 {
273  // NOTE: Clusters are never used directly in the task, so the container is neither created not retrieved
274  // Retrieve tracks
275  AliTrackContainer * tracks = static_cast<AliTrackContainer * >(GetParticleContainer("tracksForCorrelations"));
276  if (!tracks) {
277  AliError(Form("%s: Unable to retrieve tracks!", GetName()));
278  return kFALSE;
279  }
280 
281  // Retrieve jets
282  AliJetContainer * jets = GetJetContainer(0);
283  if (!jets) {
284  AliError(Form("%s: Unable to retrieve jets!", GetName()));
285  return kFALSE;
286  }
287 
288  // Get z vertex
289  Double_t zVertex=fVertex[2];
290  // Flags
291  Bool_t biasedJet = kFALSE;
292  Bool_t leadJet = kFALSE;
293  // Relative angles and distances
294  Double_t deltaPhi = 0;
295  Double_t deltaEta = 0;
296  Double_t deltaR = 0;
297  Double_t epAngle = 0;
298  // Event activity (centrality or multipilicity)
299  Double_t eventActivity = 0;
300  // Efficiency correction
301  Double_t efficiency = -999;
302  // For comparison to the current jet
303  AliEmcalJet * leadingJet = jets->GetLeadingJet();
304  // For getting the proper properties of tracks
305  AliTLorentzVector track;
306 
307  // Determine the trigger for the current event
308  UInt_t eventTrigger = RetrieveTriggerMask();
309 
310  AliDebugStream(5) << "Beginning main processing. Number of jets: " << jets->GetNJets() << ", accepted jets: " << jets->GetNAcceptedJets() << "\n";
311 
312  // Handle fast partition if selected
313  if ((eventTrigger & AliVEvent::kFastOnly) && fDisableFastPartition) {
314  AliDebugStream(4) << GetName() << ": Fast partition disabled\n";
315  if (fGeneralHistograms) {
316  fHistEventRejection->Fill("Fast Partition", 1);
317  }
318  return kFALSE;
319  }
320 
321  for (auto jet : jets->accepted()) {
322  // Selects only events that we are interested in (ie triggered)
323  if (!(eventTrigger & fTriggerType)) {
324  AliDebugStream(5) << "Rejected jets due to physics selection. Phys sel: " << std::bitset<32>(eventTrigger) << ", requested triggers: " << std::bitset<32>(fTriggerType) << " \n";
325  // We can break here - the physics selection is not going to change within an event.
326  break;
327  }
328 
329  AliDebugStream(5) << "Jet passed event selection!\nJet: " << jet->toString().Data() << "\n";
330 
331  // Require the found jet to be matched
332  // This match should be between detector and particle level MC
334  if (jet->MatchedJet()) {
335  AliDebugStream(4) << "Jet is matched!\nJet: " << jet->toString().Data() << "\n";
336  }
337  else {
338  AliDebugStream(5) << "Rejected jet because it was not matched to a external event jet.\n";
339  continue;
340  }
341  }
342 
343  // Determine event activity
344  if (fBeamType == kAA || fBeamType == kpA) {
345  eventActivity = fCent;
346  }
347  else if (fBeamType == kpp) {
348  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
349  }
350 
351  // Jet properties
352  // Determine if we have the lead jet
353  leadJet = kFALSE;
354  if (jet == leadingJet) leadJet = kTRUE;
355  biasedJet = BiasedJet(jet);
356  epAngle = GetRelativeEPAngle(jet->Phi(), fEPV0);
357 
358  // Fill jet properties
359  fHistJetEtaPhi->Fill(jet->Eta(), jet->Phi());
360  FillHist(fHistJetPt[fCentBin], jet->Pt());
361  if (biasedJet == kTRUE) {
362  FillHist(fHistJetPtBias[fCentBin], jet->Pt());
363 
364  const double triggerInfo[] = {eventActivity, jet->Pt(), epAngle};
365  fhnTrigger->Fill(triggerInfo);
366  }
367 
368  // Cut on jet pt of 15 to reduce the size of the sparses
369  if (jet->Pt() > 15) {
370 
371  AliDebugStream(4) << "Passed min jet pt cut of 15. Jet: " << jet->toString().Data() << "\n";
372  for (auto trackIter : tracks->accepted_momentum()) {
373 
374  // Get proper track proeprties
375  track.Clear();
376  track = trackIter.first;
377 
378  // Determine relative angles and distances and set the respective variables
379  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
380 
381  // Fill track properties
382  fHistTrackPt->Fill(track.Pt());
383  fHistJetHEtaPhi->Fill(deltaEta, deltaPhi);
384 
385  // Calculate single particle tracking efficiency for correlations
386  efficiency = EffCorrection(track.Eta(), track.Pt());
387  AliDebugStream(6) << GetName() << ": efficiency: " << efficiency << "\n";
388 
389  if (biasedJet == kTRUE) {
390  if(fDoLessSparseAxes) { // check if we want all dimensions
391  double triggerEntries[] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), epAngle};
392  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
393  } else {
394  double triggerEntries[] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR, epAngle};
395  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
396  }
397  }
398 
399  } //track loop
400  }//jet pt cut
401  }//jet loop
402 
403  //Prepare to do event mixing
404 
405  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
406  TObjArray* tracksClone = 0;
407 
408  if(fDoEventMixing == kTRUE){
409 
410  // event mixing
411 
412  // 1. First get an event pool corresponding in mult (cent) and
413  // zvertex to the current event. Once initialized, the pool
414  // should contain nMix (reduced) events. This routine does not
415  // pre-scan the chain. The first several events of every chain
416  // will be skipped until the needed pools are filled to the
417  // specified depth. If the pool categories are not too rare, this
418  // should not be a problem. If they are rare, you could lose
419  // statistics.
420 
421  // 2. Collect the whole pool's content of tracks into one TObjArray
422  // (bgTracks), which is effectively a single background super-event.
423 
424  // 3. The reduced and bgTracks arrays must both be passed into
425  // FillCorrelations(). Also nMix should be passed in, so a weight
426  // of 1./nMix can be applied.
427 
428  AliEventPool *pool = 0;
429  if (fBeamType == kAA || fBeamType == kpA) {//everything but pp
430  pool = fPoolMgr->GetEventPool(fCent, zVertex);
431  }
432  else if (fBeamType == kpp) {//pp only
433  pool = fPoolMgr->GetEventPool(static_cast<Double_t>(tracks->GetNTracks()), zVertex);
434  }
435 
436  if (!pool){
437  if (fBeamType == kAA || fBeamType == kpA) AliFatal(Form("No pool found for centrality = %f, zVertex = %f", fCent, zVertex));
438  else if (fBeamType == kpp) AliFatal(Form("No pool found for ntracks_pp = %d, zVertex = %f", tracks->GetNTracks(), zVertex));
439  return kTRUE;
440  }
441 
442  // The number of events in the pool
443  Int_t nMix = pool->GetCurrentNEvents();
444 
445  if(eventTrigger & fTriggerType) {
446  // check for a trigger jet
447  if (pool->IsReady() || pool->NTracksInPool() >= fMinNTracksMixedEvents || nMix >= fMinNEventsMixedEvents) {
448 
449  for (auto jet : jets->accepted()) {
450  // Require the found jet to be matched
451  // This match should be between detector and particle level MC
453  if (jet->MatchedJet()) {
454  AliDebugStream(4) << "Jet is matched!\nJet: " << jet->toString().Data() << "\n";
455  }
456  else {
457  AliDebugStream(5) << "Rejected jet because it was not matched to a external event jet.\n";
458  continue;
459  }
460  }
461 
462  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
463  eventActivity = fCent;
464  }
465  else if (fBeamType == kpp) {
466  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
467  }
468 
469  // Jet properties
470  // Determine if we have the lead jet
471  leadJet = kFALSE;
472  if (jet == leadingJet) { leadJet = kTRUE; }
473  biasedJet = BiasedJet(jet);
474  epAngle = GetRelativeEPAngle(jet->Phi(), fEPV0);
475 
476  // Make sure event contains a biased jet above our threshold (reduce stats of sparse)
477  if (jet->Pt() < 15 || biasedJet == kFALSE) continue;
478 
479  // Fill mixed-event histos here
480  for (Int_t jMix=0; jMix < nMix; jMix++) {
481  TObjArray* bgTracks = pool->GetEvent(jMix);
482 
483  for (Int_t ibg=0; ibg < bgTracks->GetEntries(); ibg++){
484  AliBasicParticle *bgTrack = static_cast<AliBasicParticle*>(bgTracks->At(ibg));
485  if(!bgTrack) {
486  AliError(Form("%s:Failed to retrieve tracks from mixed events", GetName()));
487  }
488 
489  // Fill into TLorentzVector for use with functions below
490  track.Clear();
491  track.SetPtEtaPhiE(bgTrack->Pt(), bgTrack->Eta(), bgTrack->Phi(), 0);
492 
493  // Calculate single particle tracking efficiency of mixed events for correlations
494  efficiency = EffCorrection(track.Eta(), track.Pt());
495 
496  // Phi is [-0.5*TMath::Pi(), 3*TMath::Pi()/2.]
497  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
498 
499  if (fDoLessSparseAxes) { // check if we want all the axis filled
500  double triggerEntries[] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), epAngle};
501  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
502  } else {
503  double triggerEntries[] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR, epAngle};
504  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
505  }
506  }
507  }
508  }
509  }
510  }
511 
512  if (eventTrigger & fMixingEventType) {
513  tracksClone = CloneAndReduceTrackList();
514 
515  //update pool if jet in event or not
516  pool->UpdatePool(tracksClone);
517  }
518 
519  } // end of event mixing
520 
521  return kTRUE;
522 }
523 
531 {
532  if ((jet->MaxTrackPt() > fTrackBias) || (jet->MaxClusterPt() > fClusterBias))
533  {
534  return kTRUE;
535  }
536  return kFALSE;
537 }
538 
548 void AliAnalysisTaskEmcalJetHCorrelations::GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector & particleOne, AliVParticle * particleTwo, Double_t & deltaEta, Double_t & deltaPhi, Double_t & deltaR)
549 {
550  // Define dPhi = jet.Phi() - particle.Phi() and similarly for dEta
551  // Returns deltaPhi in symmetric range so that we can calculate DeltaR.
552  deltaPhi = DeltaPhi(particleOne.Phi(), particleTwo->Phi(), -1.0*TMath::Pi(), TMath::Pi());
553  deltaEta = particleTwo->Eta() - particleOne.Eta();
554  deltaR = TMath::Sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
555 
556  // Adjust to the normal range after the DeltaR caluclation
557  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -0.5*TMath::Pi(), 3*TMath::Pi()/2.);
558 }
559 
570 {
571  Double_t dphi = (epAngle - jetAngle);
572 
573  // ran into trouble with a few dEP<-Pi so trying this...
574  if( dphi<-1*TMath::Pi() ) {
575  dphi = dphi + 1*TMath::Pi();
576  } // this assumes we are doing full jets currently
577 
578  if( (dphi>0) && (dphi<1*TMath::Pi()/2) ) {
579  // Do nothing! we are in quadrant 1
580  } else if ( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ) {
581  dphi = 1*TMath::Pi() - dphi;
582  } else if ( (dphi<0) && (dphi>-1*TMath::Pi()/2) ) {
583  dphi = std::abs(dphi);
584  } else if ( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ) {
585  dphi = dphi + 1*TMath::Pi();
586  }
587 
588  // Warn if we are not in the proper range
589  if ( dphi < 0 || dphi > TMath::Pi()/2 ) {
590  AliWarningStream() << GetName() << ": dPHI not in range [0, 0.5*Pi]!\n";
591  }
592 
593  return dphi; // dphi in [0, Pi/2]
594 }
595 
603 THnSparse* AliAnalysisTaskEmcalJetHCorrelations::NewTHnSparseF(const char* name, UInt_t entries)
604 {
605  Int_t count = 0;
606  UInt_t tmp = entries;
607  while(tmp!=0){
608  count++;
609  tmp = tmp &~ -tmp; // clear lowest bit
610  }
611 
612  TString hnTitle(name);
613  const Int_t dim = count;
614  Int_t nbins[dim];
615  Double_t xmin[dim];
616  Double_t xmax[dim];
617 
618  Int_t i=0;
619  Int_t c=0;
620  while(c<dim && i<32){
621  if(entries&(1<<i)){
622 
623  TString label("");
624  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
625  hnTitle += Form(";%s",label.Data());
626  c++;
627  }
628 
629  i++;
630  }
631  hnTitle += ";";
632 
633  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
634 }
635 
646 {
647  const Double_t pi = TMath::Pi();
648 
649  switch(iEntry){
650 
651  case 0:
652  label = "V0 centrality (%)";
653  nbins = 10;
654  xmin = 0.;
655  xmax = 100.;
656  break;
657 
658  case 1:
659  label = "Jet p_{T}";
660  nbins = 20;
661  xmin = 0.;
662  xmax = 200.;
663  break;
664 
665  case 2:
666  if(fDoWiderTrackBin) {
667  label = "Track p_{T}";
668  nbins = 40;
669  xmin = 0.;
670  xmax = 10.;
671  } else {
672  label = "Track p_{T}";
673  nbins = 100;
674  xmin = 0.;
675  xmax = 10;
676  }
677  break;
678 
679  case 3:
680  label = "#delta#eta";
681  nbins = 24;
682  xmin = -1.2;
683  xmax = 1.2;
684  break;
685 
686  case 4:
687  label = "#delta#phi";
688  nbins = 72;
689  xmin = -0.5*pi;
690  xmax = 1.5*pi;
691  break;
692 
693  case 5:
694  label = "Leading Jet";
695  nbins = 3;
696  xmin = -0.5;
697  xmax = 2.5;
698  break;
699 
700  case 6:
701  label = "Trigger track";
702  nbins = 10;
703  xmin = 0;
704  xmax = 50;
705  break;
706 
707  case 7:
708  label = "deltaR";
709  nbins = 10;
710  xmin = 0.;
711  xmax = 5.0;
712  break;
713 
714  case 8:
715  label = "Leading track";
716  nbins = 20;
717  xmin = 0;
718  xmax = 50;
719  break;
720 
721  case 9:
722  label = "Event plane angle";
723  nbins = 3;
724  xmin = 0;
725  xmax = TMath::Pi()/2.;
726  break;
727  }
728 }
729 
738 {
739  // clones a track list by using AliBasicTrack which uses much less memory (used for event mixing)
740  TObjArray* tracksClone = new TObjArray;
741  tracksClone->SetOwner(kTRUE);
742 
743  // Loop over all tracks
744  AliBasicParticle * clone = 0;
745  AliTrackContainer * tracks = GetTrackContainer("tracksForCorrelations");
746 
747  for (auto particle : tracks->accepted())
748  {
749  // Fill some QA information about the tracks
750  Int_t trackPtBin = GetTrackPtBin(particle->Pt());
751  if(trackPtBin > -1) fHistTrackEtaPhi[trackPtBin]->Fill(particle->Eta(),particle->Phi());
752 
753  // Create new particle
754  clone = new AliBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
755  // Set so that we can do comparisons using the IsEqual() function.
756  clone->SetUniqueID(particle->GetUniqueID());
757 
758  tracksClone->Add(clone);
759  }
760 
761  return tracksClone;
762 }
763 
775  return EffCorrection(trackETA, trackPT, fBeamType);
776 }
777 
794 {
795  // default (current) parameters
796  // x-variable = track pt, y-variable = track eta
797  Double_t x = trackPT;
798  Double_t y = trackETA;
799  Double_t TRefficiency = -999;
800  Int_t runNUM = fCurrentRunNumber;
801  Int_t runSwitchGood = -999;
802  //Int_t centbin = -99;
803 
804  Double_t etaaxis = 0;
805  Double_t ptaxis = 0;
806 
807  Int_t effSwitch = fDoEffCorrection;
808 
809  if (beamType != AliAnalysisTaskEmcal::kpp) {
810  if(effSwitch == 1) {
811  // Semi-Good OROC C08 Runlists
812  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;
813 
814  // Good Runlists
815  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;
816 
817  // Determine which efficiency to use.
818  // This is just a way to map all possible values of the cent bin and runSwitchGood to a unique flag.
819  // 4 is the number of cent bins, and we want to index the effSwitch starting at 2.
820  if (runSwitchGood != -999) {
821  effSwitch = 2 + runSwitchGood*4 + fCentBin;
822  }
823  }
824 
825  // set up a switch for different parameter values...
826  switch(effSwitch) {
827  case 1 :
828  // first switch value - TRefficiency not used so = 1
829  // In this case, the run number isn't in any run list, so efficiency = 1
830  TRefficiency = 1.0;
831  break;
832 
833  case 2 :
834  // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
835  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);
836  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])));
837  TRefficiency = ptaxis*etaaxis;
838  break;
839 
840  case 3 :
841  // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
842  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);
843  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])));
844  TRefficiency = ptaxis*etaaxis;
845  break;
846 
847  case 4 :
848  // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
849  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);
850  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])));
851  TRefficiency = ptaxis*etaaxis;
852  break;
853 
854  case 5 :
855  // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
856  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);
857  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])));
858  TRefficiency = ptaxis*etaaxis;
859  break;
860 
861  case 6 :
862  // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
863  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);
864  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])));
865  TRefficiency = ptaxis*etaaxis;
866  break;
867 
868  case 7 :
869  // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
870  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);
871  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])));
872  TRefficiency = ptaxis*etaaxis;
873  break;
874 
875  case 8 :
876  // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
877  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);
878  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])));
879  TRefficiency = ptaxis*etaaxis;
880  break;
881 
882  case 9 :
883  // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
884  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);
885  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])));
886  TRefficiency = ptaxis*etaaxis;
887  break;
888 
889  default :
890  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
891  // ie. The efficiency correction is disabled.
892  TRefficiency = 1.0;
893  }
894  }
895  else {
896  // Track efficiency for pp
897  // Calculated using LHC12f1a. See analysis note for more details!
898 
899  if (fDoEffCorrection != 0) {
900  // If the trackPt > 6 GeV, then all we need is this coefficient
901  Double_t coefficient = 0.898052; // p6
902  if (trackPT < 6) {
903  coefficient = (1 + -0.442232 * trackPT // p0
904  + 0.501831 * std::pow(trackPT, 2) // p1
905  + -0.252024 * std::pow(trackPT, 3) // p2
906  + 0.062964 * std::pow(trackPT, 4) // p3
907  + -0.007681 * std::pow(trackPT, 5) // p4
908  + 0.000365 * std::pow(trackPT, 6)); // p5
909  }
910 
911  // Calculate track eff
912  TRefficiency = coefficient * (1 + 0.402825 * std::abs(trackETA) // p7
913  + -2.213152 * std::pow(trackETA, 2) // p8
914  + 4.311098 * std::abs(std::pow(trackETA, 3)) // p9
915  + -2.778200 * std::pow(trackETA, 4)); // p10
916  }
917  else {
918  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
919  TRefficiency = 1;
920  }
921  }
922 
923  return TRefficiency;
924 }
925 
935 void AliAnalysisTaskEmcalJetHCorrelations::FillHist(TH1 * hist, Double_t fillValue, Double_t weight, Bool_t noCorrection)
936 {
937  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
938  {
939  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
940  hist->Fill(fillValue, weight);
941  }
942  else
943  {
944  // Determine where to get the values in the correction hist
945  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(fillValue);
946 
947  std::vector <Double_t> yBinsContent;
948  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), fillValue));
949  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
950  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
951 
952  // Loop over all possible bins to contribute.
953  // If content is 0 then calling Fill won't make a difference
954  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
955  {
956  // Don't bother trying to fill in the weight is 0
957  if (yBinsContent.at(index-1) > 0) {
958  // Determine the value to fill based on the center of the bins.
959  // This in principle allows the binning between the correction and hist to be different
960  Double_t fillLocation = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
961  AliDebug(4, TString::Format("fillLocation: %f, weight: %f", fillLocation, yBinsContent.at(index-1)));
962  // minus 1 since loop starts at 1
963  hist->Fill(fillLocation, weight*yBinsContent.at(index-1));
964  }
965  }
966 
967  //TEMP
968  //hist->Draw();
969  //END TEMP
970  }
971 }
972 
983 void AliAnalysisTaskEmcalJetHCorrelations::FillHist(THnSparse * hist, Double_t *fillValue, Double_t weight, Bool_t noCorrection)
984 {
985  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
986  {
987  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
988  hist->Fill(fillValue, weight);
989  }
990  else
991  {
992  // Jet pt is always located in the second position
993  Double_t jetPt = fillValue[1];
994 
995  // Determine where to get the values in the correction hist
996  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(jetPt);
997 
998  std::vector <Double_t> yBinsContent;
999  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), jetPt));
1000  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
1001  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
1002 
1003  // Loop over all possible bins to contribute.
1004  // If content is 0 then calling Fill won't make a difference
1005  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
1006  {
1007  // Don't bother trying to fill in the weight is 0
1008  if (yBinsContent.at(index-1) > 0) {
1009  // Determine the value to fill based on the center of the bins.
1010  // This in principle allows the binning between the correction and hist to be different
1011  fillValue[1] = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
1012  AliDebug(4,TString::Format("fillValue[1]: %f, weight: %f", fillValue[1], yBinsContent.at(index-1)));
1013  // minus 1 since loop starts at 1
1014  hist->Fill(fillValue, weight*yBinsContent.at(index-1));
1015  }
1016  }
1017  }
1018 }
1019 
1029 void AliAnalysisTaskEmcalJetHCorrelations::AccessSetOfYBinValues(TH2D * hist, Int_t xBin, std::vector <Double_t> & yBinsContent, Double_t scaleFactor)
1030 {
1031  for (Int_t index = 1; index <= hist->GetYaxis()->GetNbins(); index++)
1032  {
1033  //yBinsContent[index-1] = hist->GetBinContent(hist->GetBin(xBin,index));
1034  yBinsContent.push_back(hist->GetBinContent(hist->GetBin(xBin,index)));
1035 
1036  if (scaleFactor >= 0)
1037  {
1038  // -1 since index starts at 1
1039  hist->SetBinContent(hist->GetBin(xBin,index), yBinsContent.at(index-1)/scaleFactor);
1040  }
1041  }
1042 }
1043 
1061 {
1062  // Initialize grid connection if necessary
1063  if (filename.Contains("alien://") && !gGrid) {
1064  TGrid::Connect("alien://");
1065  }
1066 
1067  // Setup hist name if a track or cluster bias was defined.
1068  // NOTE: This can always be disabled by setting kDisableBias.
1069  // We arbitrarily add 0.1 to test since the values are doubles and cannot be
1070  // tested directly for equality. If we are still less than disable bins, then
1071  // it has been set and we should format it.
1072  // NOTE: To ensure we can disable, we don't just take the member values!
1073  // NOTE: The histBaseName will be attempted if the formatted name cannot be found.
1074  TString histBaseName = histName;
1076  histName = TString::Format("%s_Track%.2f", histName.Data(), trackBias);
1077  }
1078  if (clusterBias + 0.1 < AliAnalysisTaskEmcalJetHCorrelations::kDisableBias) {
1079  histName = TString::Format("%s_Clus%.2f", histName.Data(), clusterBias);
1080  }
1081 
1082  // Open file containing the correction
1083  TFile * jesCorrectionFile = TFile::Open(filename);
1084  if (!jesCorrectionFile || jesCorrectionFile->IsZombie()) {
1085  AliError(TString::Format("%s: Could not open JES correction file %s", GetName(), filename.Data()));
1086  return kFALSE;
1087  }
1088 
1089  // Retrieve the histogram containing the correction and safely add it to the task.
1090  TH2D * JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histName.Data()));
1091  if (JESCorrectionHist) {
1092  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histName.Data(), filename.Data()));
1093  }
1094  else {
1095  AliError(TString::Format("%s: JES correction hist name \"%s\" not found in file %s.", GetName(), histName.Data(), filename.Data()));
1096 
1097  // Attempt the base name instead of the formatted hist name
1098  JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histBaseName.Data()));
1099  if (JESCorrectionHist) {
1100  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histBaseName.Data(), filename.Data()));
1101  histName = histBaseName;
1102  }
1103  else
1104  {
1105  AliError(TString::Format("%s: JES correction with base hist name %s not found in file %s.", GetName(), histBaseName.Data(), filename.Data()));
1106  return kFALSE;
1107  }
1108  }
1109 
1110  // Clone to ensure that the hist is available
1111  TH2D * tempHist = static_cast<TH2D *>(JESCorrectionHist->Clone());
1112  tempHist->SetDirectory(0);
1113  SetJESCorrectionHist(tempHist);
1114 
1115  // Close file
1116  jesCorrectionFile->Close();
1117 
1118  // Append to task name for clarity
1119  // Unfortunately, this doesn't change the name of the output list (it would need to be
1120  // changed in the AnalysisManager output container), so the suffix is still important
1121  // if this correction is manually configured!
1122  TString tempName = GetName();
1123  TString tag = "_JESCorr";
1124  // Append the tag if it isn't already included
1125  if (tempName.Index(tag) == -1) {
1126  // Insert before the suffix
1127  Ssiz_t suffixLocation = tempName.Last('_');
1128  tempName.Insert(suffixLocation, tag.Data());
1129 
1130  // Set the new name
1131  AliDebug(3, TString::Format("%s: Setting task name to %s", GetName(), tempName.Data()));
1132  SetName(tempName.Data());
1133  }
1134 
1135  // Successful
1136  return kTRUE;
1137 }
1138 
1144  const char *nTracks,
1145  const char *nCaloClusters,
1146  // Jet options
1147  const Double_t trackBias,
1148  const Double_t clusterBias,
1149  // Mixed event options
1150  const Int_t nTracksMixedEvent, // Additionally acts as a switch for enabling mixed events
1151  const Int_t minNTracksMixedEvent,
1152  const Int_t minNEventsMixedEvent,
1153  const UInt_t nCentBinsMixedEvent,
1154  // Triggers
1155  UInt_t trigEvent,
1156  UInt_t mixEvent,
1157  // Options
1158  const Bool_t lessSparseAxes,
1159  const Bool_t widerTrackBin,
1160  // Corrections
1161  const Int_t doEffCorrSW,
1162  const Bool_t JESCorrection,
1163  const char * JESCorrectionFilename,
1164  const char * JESCorrectionHistName,
1165  const char *suffix
1166  )
1167 {
1168  // Get the pointer to the existing analysis manager via the static access method.
1169  //==============================================================================
1170  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1171  if (!mgr)
1172  {
1173  AliErrorClass("No analysis manager to connect to.");
1174  return nullptr;
1175  }
1176 
1177  //-------------------------------------------------------
1178  // Init the task and do settings
1179  //-------------------------------------------------------
1180 
1181  // Determine cluster and track names
1182  TString trackName(nTracks);
1183  TString clusName(nCaloClusters);
1184 
1185  if (trackName == "usedefault") {
1187  }
1188 
1189  if (clusName == "usedefault") {
1191  }
1192 
1193  TString name("AliAnalysisTaskJetH");
1194  if (!trackName.IsNull()) {
1195  name += TString::Format("_%s", trackName.Data());
1196  }
1197  if (!clusName.IsNull()) {
1198  name += TString::Format("_%s", clusName.Data());
1199  }
1200  if (strcmp(suffix, "") != 0) {
1201  name += TString::Format("_%s", suffix);
1202  }
1203 
1205  // Set jet bias
1206  correlationTask->SetTrackBias(trackBias);
1207  correlationTask->SetClusterBias(clusterBias);
1208  // Mixed events
1209  correlationTask->SetEventMixing(static_cast<Bool_t>(nTracksMixedEvent));
1210  correlationTask->SetNumberOfMixingTracks(nTracksMixedEvent);
1211  correlationTask->SetMinNTracksForMixedEvents(minNTracksMixedEvent);
1212  correlationTask->SetMinNEventsForMixedEvents(minNEventsMixedEvent);
1213  correlationTask->SetNCentBinsMixedEvent(nCentBinsMixedEvent);
1214  // Triggers
1215  correlationTask->SetTriggerType(trigEvent);
1216  correlationTask->SetMixedEventTriggerType(mixEvent);
1217  // Options
1218  correlationTask->SetNCentBins(5);
1219  correlationTask->SetVzRange(-10,10);
1220  correlationTask->SetDoLessSparseAxes(lessSparseAxes);
1221  correlationTask->SetDoWiderTrackBin(widerTrackBin);
1222  // Corrections
1223  correlationTask->SetDoEffCorr(doEffCorrSW);
1224  if (JESCorrection == kTRUE)
1225  {
1226  Bool_t result = correlationTask->RetrieveAndInitializeJESCorrectionHist(JESCorrectionFilename, JESCorrectionHistName, correlationTask->GetTrackBias(), correlationTask->GetClusterBias());
1227  if (!result) {
1228  AliErrorClass("Failed to successfully retrieve and initialize the JES correction! Task initialization continuing without JES correction (can be set manually later).");
1229  }
1230  }
1231 
1232  //-------------------------------------------------------
1233  // Final settings, pass to manager and set the containers
1234  //-------------------------------------------------------
1235 
1236  mgr->AddTask(correlationTask);
1237 
1238  // Create containers for input/output
1239  mgr->ConnectInput (correlationTask, 0, mgr->GetCommonInputContainer() );
1240  AliAnalysisDataContainer * cojeth = mgr->CreateContainer(correlationTask->GetName(),
1241  TList::Class(),
1242  AliAnalysisManager::kOutputContainer,
1243  Form("%s", AliAnalysisManager::GetCommonFileName()));
1244  mgr->ConnectOutput(correlationTask, 1, cojeth);
1245 
1246  return correlationTask;
1247 }
1248 
1255  std::string clusName,
1256  const double jetConstituentPtCut,
1257  const double trackEta,
1258  const double jetRadius)
1259 {
1260  bool returnValue = false;
1261  AliInfoStream() << "Configuring Jet-H Correlations task for a standard analysis.\n";
1262 
1263  // Add Containers
1264  // Clusters
1265  if (clusName == "usedefault") {
1267  }
1268  // For jet finding
1269  AliClusterContainer * clustersForJets = new AliClusterContainer(clusName.c_str());
1270  clustersForJets->SetName("clustersForJets");
1271  clustersForJets->SetMinE(jetConstituentPtCut);
1272 
1273  // Tracks
1274  // For jet finding
1275  if (trackName == "usedefault") {
1277  }
1278  AliParticleContainer * particlesForJets = CreateParticleOrTrackContainer(trackName.c_str());
1279  particlesForJets->SetName("particlesForJets");
1280  particlesForJets->SetMinPt(jetConstituentPtCut);
1281  particlesForJets->SetEtaLimits(-1.0*trackEta, trackEta);
1282  // Don't need to adopt the container - we'll just use it to find the right jet collection
1283  // For correlations
1284  AliParticleContainer * particlesForCorrelations = CreateParticleOrTrackContainer(trackName.c_str());
1285  if (particlesForCorrelations)
1286  {
1287  particlesForCorrelations->SetName("tracksForCorrelations");
1288  particlesForCorrelations->SetMinPt(0.15);
1289  particlesForCorrelations->SetEtaLimits(-1.0*trackEta, trackEta);
1290  // Adopt the container
1291  this->AdoptParticleContainer(particlesForCorrelations);
1292  }
1293  else {
1294  AliWarningStream() << "No particle container was successfully created!\n";
1295  }
1296 
1297  // Jets
1301  jetRadius,
1303  particlesForJets,
1304  clustersForJets);
1305  // 0.6 * jet area
1306  jetContainer->SetJetAreaCut(jetRadius * jetRadius * TMath::Pi() * 0.6);
1307  jetContainer->SetMaxTrackPt(100);
1308  jetContainer->SetJetPtCut(0.1);
1309 
1310  // Successfully configured
1311  returnValue = true;
1312 
1313  return returnValue;
1314 }
1315 
1322  std::string clusName,
1323  const double jetConstituentPtCut,
1324  const double trackEta,
1325  const double jetRadius,
1326  const std::string & jetTag,
1327  const std::string & correlationsTracksCutsPeriod)
1328 {
1329  bool returnValue = false;
1330  AliInfoStream() << "Configuring Jet-H Correlations task for an embedding analysis.\n";
1331 
1332  // Set the task to know it that is embedded
1333  this->SetIsEmbedded(true);
1334 
1335  // Add Containers
1336  // Clusters
1337  if (clusName == "usedefault") {
1339  }
1340  // For jet finding
1341  AliClusterContainer * clustersForJets = new AliClusterContainer(clusName.c_str());
1342  clustersForJets->SetName("clustersForJets");
1343  clustersForJets->SetMinE(jetConstituentPtCut);
1344  // We need the combined clusters, which should be available in the internal event.
1345  // However, we don't need to adopt the container - we'll just use it to find the right jet collection
1346  // For correlations
1347  /*AliClusterContainer * clustersforCorrelations = new AliClusterContainer("usedefault");
1348  clustersForCorrelations->SetName("clustersForCorrelations");
1349  clustersForCorrelations->SetMinE(0.30);
1350  clustersForCorrelations->SetIsEmbedding(true);
1351  this->AdoptClusterContainer(clustersForCorrelations);*/
1352 
1353  // Tracks
1354  // For jet finding
1355  if (trackName == "usedefault") {
1357  }
1358  AliParticleContainer * particlesForJets = CreateParticleOrTrackContainer(trackName.c_str());
1359  particlesForJets->SetName("particlesForJets");
1360  particlesForJets->SetMinPt(jetConstituentPtCut);
1361  particlesForJets->SetEtaLimits(-1.0*trackEta, trackEta);
1362  // Don't need to adopt the container - we'll just use it to find the right jet collection
1363  // For correlations
1364  AliParticleContainer * particlesForCorrelations = CreateParticleOrTrackContainer(trackName.c_str());
1365  // Ensure that we don't operate on a null pointer
1366  if (particlesForCorrelations)
1367  {
1368  particlesForCorrelations->SetName("tracksForCorrelations");
1369  particlesForCorrelations->SetMinPt(0.15);
1370  particlesForCorrelations->SetEtaLimits(-1.0*trackEta, trackEta);
1371  particlesForCorrelations->SetIsEmbedding(true);
1372  AliTrackContainer * trackCont = dynamic_cast<AliTrackContainer *>(particlesForCorrelations);
1373  if (trackCont) {
1374  // This option only exists for track containers
1375  trackCont->SetTrackCutsPeriod(correlationsTracksCutsPeriod.c_str());
1376  }
1377  // Adopt the container
1378  this->AdoptParticleContainer(particlesForCorrelations);
1379  }
1380  else {
1381  AliWarningStream() << "No particle container was successfully created!\n";
1382  }
1383 
1384  // Jets
1385  // The tag "hybridJets" is defined in the jet finder
1389  jetRadius,
1391  particlesForJets,
1392  clustersForJets,
1393  jetTag);
1394  // 0.6 * jet area
1395  jetContainer->SetJetAreaCut(jetRadius * jetRadius * TMath::Pi() * 0.6);
1396  jetContainer->SetMaxTrackPt(100);
1397  jetContainer->SetJetPtCut(0.1);
1398 
1399  // Successfully configured
1400  returnValue = true;
1401 
1402  return returnValue;
1403 }
1404 
1413 {
1414  AliParticleContainer * partCont = 0;
1416  AliTrackContainer * trackCont = new AliTrackContainer(collectionName.c_str());
1417  partCont = trackCont;
1418  }
1419  else if (collectionName != "") {
1420  partCont = new AliParticleContainer(collectionName.c_str());
1421  }
1422 
1423  return partCont;
1424 }
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)
Double_t GetRelativeEPAngle(Double_t jetAngle, Double_t epAngle) const
void FillHist(TH1 *hist, Double_t fillValue, Double_t weight=1.0, Bool_t noCorrection=kFALSE)
const char * filename
Definition: TestFCM.C:1
void SetTrackCutsPeriod(const char *period)
virtual void SetTrackBias(Double_t b)
Require a track with pt > b in jet.
double Double_t
Definition: External.C:58
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.
Bool_t fRequireMatchedJetWhenEmbedding
True if jets are required to be matched (ie. jet->MatchedJet() != nullptr)
Double_t fEPV0
!event plane V0
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
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
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", const std::string &correlationsTracksCutsPeriod="lhc11a")
Int_t GetNJets() const
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
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()