AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetHMEC.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 "AliBasicParticle.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliTLorentzVector.h"
23 #include "AliLog.h"
24 
25 #include "AliClusterContainer.h"
26 #include "AliTrackContainer.h"
27 
31 
32 // 0-10% centrality: Semi-Good Runs
33 Double_t AliAnalysisTaskEmcalJetHMEC::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};
34 // 10-30% centrality: Semi-Good Runs
35 Double_t AliAnalysisTaskEmcalJetHMEC::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};
36 // 30-50% centrality: Semi-Good Runs
37 Double_t AliAnalysisTaskEmcalJetHMEC::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};
38 // 50-90% centrality: Semi-Good Runs
39 Double_t AliAnalysisTaskEmcalJetHMEC::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};
40 
41 // 0-10% centrality: Good Runs
42 Double_t AliAnalysisTaskEmcalJetHMEC::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};
43 // 10-30% centrality: Good Runs
44 Double_t AliAnalysisTaskEmcalJetHMEC::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};
45 // 30-50% centrality: Good Runs
46 Double_t AliAnalysisTaskEmcalJetHMEC::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};
47 // 50-90% centrality: Good Runs
48 Double_t AliAnalysisTaskEmcalJetHMEC::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};
49 
50 //________________________________________________________________________
52  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetHMEC", kFALSE),
53  fTrackBias(5),
54  fClusterBias(5),
55  fDoEventMixing(kFALSE),
56  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
57  fPoolMgr(0),
58  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
59  fDisableFastPartition(kFALSE),
60  fDoEffCorrection(0),
61  fJESCorrectionHist(0),
62  fNoMixedEventJESCorrection(kFALSE),
63  fDoLessSparseAxes(0), fDoWiderTrackBin(0),
64  fHistTrackPt(0),
65  fHistJetEtaPhi(0),
66  fHistJHPsi(0),
67  fHistJetHEtaPhi(0),
68  fhnMixedEvents(0),
69  fhnJH(0)
70 {
71  // Default Constructor
73 }
74 
75 //________________________________________________________________________
77  AliAnalysisTaskEmcalJet(name, kTRUE),
78  fTrackBias(5),
79  fClusterBias(5),
80  fDoEventMixing(kFALSE),
81  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
82  fPoolMgr(0),
83  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
84  fDisableFastPartition(kFALSE),
85  fDoEffCorrection(0),
86  fJESCorrectionHist(0),
87  fNoMixedEventJESCorrection(kFALSE),
88  fDoLessSparseAxes(0), fDoWiderTrackBin(0),
89  fHistTrackPt(0),
90  fHistJetEtaPhi(0),
91  fHistJHPsi(0),
92  fHistJetHEtaPhi(0),
93  fhnMixedEvents(0),
94  fhnJH(0)
95 {
96  // Constructor
98  // Ensure that additional general histograms are created
100 }
101 
102 //________________________________________________________________________
104 {
105  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; trackPtBin++){
106  fHistTrackEtaPhi[trackPtBin]=0;
107  }
108  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
109  fHistJetPt[centralityBin]=0;
110  fHistJetPtBias[centralityBin]=0;
111  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
112  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
113  fHistJetH[centralityBin][jetPtBin][etaBin]=0;
114  fHistJetHBias[centralityBin][jetPtBin][etaBin]=0;
115  }
116  }
117  }
118 }
119 
120 //________________________________________________________________________
122  // Called once
124 
125  // Create histograms
126  fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
127  fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
128  fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
129 
130  fHistJHPsi = new TH3F("fHistJHPsi","Jet-Hadron ntr-trpt-dpsi",20,0,100,200,0,20,120,0,180);
131 
132  fOutput->Add(fHistTrackPt);
133  fOutput->Add(fHistJetEtaPhi);
134  fOutput->Add(fHistJetHEtaPhi);
135  fOutput->Add(fHistJHPsi);
136 
137  TString name;
138 
139  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; ++trackPtBin){
140  name = Form("fHistTrackEtaPhi_%i", trackPtBin);
141  fHistTrackEtaPhi[trackPtBin] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
142  fOutput->Add(fHistTrackEtaPhi[trackPtBin]);
143  }
144 
145  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
146  name = Form("fHistJetPt_%i",centralityBin);
147  fHistJetPt[centralityBin] = new TH1F(name,name,200,0,200);
148  fOutput->Add(fHistJetPt[centralityBin]);
149 
150  name = Form("fHistJetPtBias_%i",centralityBin);
151  fHistJetPtBias[centralityBin] = new TH1F(name,name,200,0,200);
152  fOutput->Add(fHistJetPtBias[centralityBin]);
153 
154  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
155  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
156  name = Form("fHistJetH_%i_%i_%i",centralityBin,jetPtBin,etaBin);
157  fHistJetH[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
158  fOutput->Add(fHistJetH[centralityBin][jetPtBin][etaBin]);
159 
160  name = Form("fHistJetHBias_%i_%i_%i",centralityBin,jetPtBin,etaBin);
161  fHistJetHBias[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
162  fOutput->Add(fHistJetHBias[centralityBin][jetPtBin][etaBin]);
163  }
164  }
165  }
166 
167  UInt_t cifras = 0; // bit coded, see GetDimParams() below
168  if(fDoLessSparseAxes) {
169  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
170  } else {
171  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
172  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
173  }
174  fhnJH = NewTHnSparseF("fhnJH", cifras);
175  fhnJH->Sumw2();
176  fOutput->Add(fhnJH);
177 
178  if(fDoEventMixing){
179  if(fDoLessSparseAxes) {
180  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
181  } else {
182  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
183  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
184  }
185  fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
186  fhnMixedEvents->Sumw2();
187  fOutput->Add(fhnMixedEvents);
188  }
189 
190  PostData(1, fOutput);
191 
192  // Event Mixing
193  Int_t poolSize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
194  // ZVertex
195  Int_t nZVertexBins = 10;
196  Double_t* zVertexBins = GenerateFixedBinArray(nZVertexBins, -10, 10);
197  // Event activity (centrality of multiplicity)
198  Int_t nEventActivityBins = 8;
199  Double_t* eventActivityBins = 0;
200  // +1 to accomodate the fact that we define bins rather than array entries.
201  Double_t multiplicityBins[kMixedEventMulitplictyBins+1] = {0., 4., 9., 15., 25., 35., 55., 100., 500.};
202 
203  // Cannot use GetBeamType() since it is not available until UserExec()
204  if (fForceBeamType != AliAnalysisTaskEmcal::kpp ) { //all besides pp
205  // Event Activity is centrality in AA, pA
206  nEventActivityBins = fNCentBinsMixedEvent;
207  eventActivityBins = GenerateFixedBinArray(nEventActivityBins, 0, 100);
208  }
209  else if (fForceBeamType == AliAnalysisTaskEmcal::kpp) { //for pp only
210  // Event Activity is multiplicity in pp
211  eventActivityBins = multiplicityBins;
212  }
213 
214  fPoolMgr = new AliEventPoolManager(poolSize, fNMixingTracks, nEventActivityBins, eventActivityBins, nZVertexBins, zVertexBins);
215 }
216 
217 //________________________________________________________________________
219  // Get eta bin for histos.
220 
221  Int_t etabin = -1;
222  eta = TMath::Abs(eta);
223  if (eta <= 0.4) etabin = 0;
224  else if (eta > 0.4 && eta < 0.8) etabin = 1;
225  else if (eta >= 0.8) etabin = 2;
226  return etabin;
227 }
228 
229 //________________________________________________________________________
231 {
232  Int_t ptBin = -1;
233  if (pt < 0.5) ptBin = 0;
234  else if (pt < 1 ) ptBin = 1;
235  else if (pt < 2 ) ptBin = 2;
236  else if (pt < 3 ) ptBin = 3;
237  else if (pt < 5 ) ptBin = 4;
238  else if (pt < 8 ) ptBin = 5;
239  else if (pt < 20 ) ptBin = 6;
240 
241  return ptBin;
242 }
243 
244 //________________________________________________________________________
246 {
247  // Get jet pt bin for histos.
248 
249  Int_t ptBin = -1;
250  if (pt >= 15 && pt < 20) ptBin = 0;
251  else if (pt >= 20 && pt < 25) ptBin = 1;
252  else if (pt >= 25 && pt < 30) ptBin = 2;
253  else if (pt >= 30 && pt < 60) ptBin = 3;
254  else if (pt >= 60) ptBin = 4;
255 
256  return ptBin;
257 }
258 
259 //________________________________________________________________________
262 }
263 
264 //________________________________________________________________________
266  // Main loop called for each event
267 
268  // Retrieve clusters
270  if (!clusters) {
271  AliError(Form("%s: Unable to retrieve clusters!", GetName()));
272  return kFALSE;
273  }
274 
275  // Retrieve tracks
276  AliTrackContainer * tracks = static_cast<AliTrackContainer * >(GetParticleContainer("tracksForCorrelations"));
277  if (!tracks) {
278  AliError(Form("%s: Unable to retrieve tracks!", GetName()));
279  return kFALSE;
280  }
281 
282  // Retrieve jets
283  AliJetContainer * jets = GetJetContainer(0);
284  if (!jets) {
285  AliError(Form("%s: Unable to retrieve jets!", GetName()));
286  return kFALSE;
287  }
288 
289  // Used to calculate the angle betwene the jet and the hadron
290  TVector3 jetVector;
291  // Get z vertex
292  Double_t zVertex=fVertex[2];
293  // Flags
294  Bool_t biasedJet = kFALSE;
295  Bool_t leadJet = kFALSE;
296  // Relative angles and distances
297  Double_t deltaPhi = 0;
298  Double_t deltaEta = 0;
299  Double_t deltaR = 0;
300  // Event activity (centrality or multipilicity)
301  Double_t eventActivity = 0;
302  // Efficiency correction
303  Double_t efficiency = -999;
304  // Determining bins for histogram indices
305  Int_t jetPtBin = -1;
306  Int_t etaBin = -1;
307  // For comparison to the current jet
308  AliEmcalJet * leadingJet = jets->GetLeadingJet();
309  // For getting the proper properties of tracks
310  AliTLorentzVector track;
311 
312  // Determine the trigger for the current event
313  UInt_t eventTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
314 
315  // Handle fast partition if selected
316  if ((eventTrigger & AliVEvent::kFastOnly) && fDisableFastPartition) {
317  AliDebug(4, TString::Format("%s: Fast partition disabled", GetName()));
318  if (fGeneralHistograms) {
319  fHistEventRejection->Fill("Fast Partition", 1);
320  }
321  return kFALSE;
322  }
323 
324  for (auto jet : jets->accepted()) {
325  // Selects only events that we are interested in (ie triggered)
326  if (!(eventTrigger & fTriggerType)) continue;
327  AliDebug(5, TString::Format("%s: Jet accepted!\nJet: %s", GetName(), jet->toString().Data()));
328 
329  // Jet properties
330  // Determine if we have the lead jet
331  leadJet = kFALSE;
332  if (jet == leadingJet) leadJet = kTRUE;
333 
334  // Determine if the jet is biased
335  biasedJet = BiasedJet(jet);
336 
337  // Calculate vector
338  jetVector.SetXYZ(jet->Px(), jet->Py(), jet->Pz());
339 
340  // Fill jet properties
341  FillHist(fHistJetPt[fCentBin], jet->Pt());
342  if (biasedJet == kTRUE) {
343  FillHist(fHistJetPtBias[fCentBin], jet->Pt());
344  }
345 
346  fHistJetEtaPhi->Fill(jet->Eta(), jet->Phi());
347 
348  if (jet->Pt() > 15) {
349 
350  for (auto trackIter : tracks->accepted_momentum()) {
351 
352  // Get proper track proeprties
353  track.Clear();
354  track = trackIter.first;
355 
356  // Determine relative angles and distances and set the respective variables
357  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
358 
359  // Determine bins for filling histograms
360  // jet Pt
361  jetPtBin = GetJetPtBin(jet->Pt());
362  if (jetPtBin < 0)
363  {
364  AliError(Form("Jet Pt Bin negative: %f", jet->Pt()));
365  continue;
366  }
367  // eta
368  etaBin = GetEtaBin(deltaEta);
369  if (etaBin < 0) {
370  AliError(Form("Eta Bin negative: %f", deltaEta));
371  continue;
372  }
373 
374  // Fill track properties
375  fHistTrackPt->Fill(track.Pt());
376 
377  if ( (jet->Pt() > 20.) && (jet->Pt() < 60.) ) {
378  fHistJHPsi->Fill(tracks->GetNTracks(), track.Pt(), track.Vect().Angle(jetVector) * TMath::RadToDeg() );
379  }
380 
381  fHistJetH[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
382  fHistJetHEtaPhi->Fill(deltaEta, deltaPhi);
383 
384  // Calculate single particle tracking efficiency for correlations
385  efficiency = EffCorrection(track.Eta(), track.Pt());
386  AliDebug(6, TString::Format("%s: efficiency: %f", GetName(), efficiency));
387 
388  if (biasedJet == kTRUE) {
389  fHistJetHBias[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
390 
391  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
392  eventActivity = fCent;
393  }
394  else if (fBeamType == kpp) {
395  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
396  }
397 
398  if(fDoLessSparseAxes) { // check if we want all dimensions
399  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
400  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
401  } else {
402  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
403  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
404  }
405  }
406 
407  } //track loop
408  }//jet pt cut
409  }//jet loop
410 
411  //Prepare to do event mixing
412 
413  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
414  TObjArray* tracksClone = 0;
415 
416  if(fDoEventMixing == kTRUE){
417 
418  // event mixing
419 
420  // 1. First get an event pool corresponding in mult (cent) and
421  // zvertex to the current event. Once initialized, the pool
422  // should contain nMix (reduced) events. This routine does not
423  // pre-scan the chain. The first several events of every chain
424  // will be skipped until the needed pools are filled to the
425  // specified depth. If the pool categories are not too rare, this
426  // should not be a problem. If they are rare, you could lose
427  // statistics.
428 
429  // 2. Collect the whole pool's content of tracks into one TObjArray
430  // (bgTracks), which is effectively a single background super-event.
431 
432  // 3. The reduced and bgTracks arrays must both be passed into
433  // FillCorrelations(). Also nMix should be passed in, so a weight
434  // of 1./nMix can be applied.
435 
436  AliEventPool *pool = 0;
437  if (fBeamType == kAA || fBeamType == kpA) {//everything but pp
438  pool = fPoolMgr->GetEventPool(fCent, zVertex);
439  }
440  else if (fBeamType == kpp) {//pp only
441  pool = fPoolMgr->GetEventPool(static_cast<Double_t>(tracks->GetNTracks()), zVertex);
442  }
443 
444  if (!pool){
445  if (fBeamType == kAA || fBeamType == kpA) AliFatal(Form("No pool found for centrality = %f, zVertex = %f", fCent, zVertex));
446  else if (fBeamType == kpp) AliFatal(Form("No pool found for ntracks_pp = %d, zVertex = %f", tracks->GetNTracks(), zVertex));
447  return kTRUE;
448  }
449 
450  // The number of events in the pool
451  Int_t nMix = pool->GetCurrentNEvents();
452 
453  if(eventTrigger & fTriggerType) {
454  // check for a trigger jet
455  if (pool->IsReady() || pool->NTracksInPool() >= fMinNTracksMixedEvents || nMix >= fMinNEventsMixedEvents) {
456 
457  for (auto jet : jets->accepted()) {
458 
459  // Jet properties
460  // Determine if we have the lead jet
461  leadJet = kFALSE;
462  if (jet == leadingJet) { leadJet = kTRUE; }
463 
464  // Determine if the jet is biased
465  biasedJet = BiasedJet(jet);
466 
467  // Make sure event contains jet above our threshold (reduce stats of sparse)
468  if (jet->Pt() < 15) continue;
469 
470  // Fill for biased jet triggers only
471  if (biasedJet == kTRUE) {
472 
473  // Fill mixed-event histos here
474  for (Int_t jMix=0; jMix < nMix; jMix++) {
475  TObjArray* bgTracks = pool->GetEvent(jMix);
476 
477  for(Int_t ibg=0; ibg < bgTracks->GetEntries(); ibg++){
478 
479  AliBasicParticle *bgTrack = static_cast<AliBasicParticle*>(bgTracks->At(ibg));
480  if(!bgTrack)
481  {
482  AliError(Form("%s:Failed to retrieve tracks from mixed events", GetName()));
483  }
484 
485  // Fill into TLorentzVector for use with functions below
486  track.Clear();
487  track.SetPtEtaPhiE(bgTrack->Pt(), bgTrack->Eta(), bgTrack->Phi(), 0);
488 
489  // Calculate single particle tracking efficiency of mixed events for correlations
490  efficiency = EffCorrection(track.Eta(), track.Pt());
491 
492  // Phi is [-0.5*TMath::Pi(), 3*TMath::Pi()/2.]
493  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
494 
495  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
496  eventActivity = fCent;
497  }
498  else if (fBeamType == kpp) {
499  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
500  }
501 
502  if(fDoLessSparseAxes) { // check if we want all the axis filled
503  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
504  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
505  } else {
506  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
507  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), fNoMixedEventJESCorrection);
508  }
509  }
510  }
511  }
512  }
513  }
514  }
515 
516  if(eventTrigger & fMixingEventType) {
517  tracksClone = CloneAndReduceTrackList();
518 
519  //update pool if jet in event or not
520  pool->UpdatePool(tracksClone);
521  }
522 
523  } // end of event mixing
524 
525  return kTRUE;
526 }
527 
528 //________________________________________________________________________
530 {
531  //just terminate
532 }
533 
534 //________________________________________________________________________
536 {
537  if ((jet->MaxTrackPt() > fTrackBias) || (jet->MaxClusterPt() > fClusterBias))
538  {
539  return kTRUE;
540  }
541  return kFALSE;
542 }
543 
544 //________________________________________________________________________
545 void AliAnalysisTaskEmcalJetHMEC::GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector & particleOne, AliVParticle * particleTwo, Double_t & deltaEta, Double_t & deltaPhi, Double_t & deltaR)
546 {
547  // TODO: Understand order of arguments to DeltaPhi vs DeltaEta
548  // Returns deltaPhi in symmetric range so that we can calculate DeltaR.
549  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -1.0*TMath::Pi(), TMath::Pi());
550  deltaEta = particleOne.Eta() - particleTwo->Eta();
551  deltaR = TMath::Sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
552 
553  // Adjust to the normal range after the DeltaR caluclation
554  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -0.5*TMath::Pi(), 3*TMath::Pi()/2.);
555 }
556 
557 //________________________________________________________________________
558 THnSparse* AliAnalysisTaskEmcalJetHMEC::NewTHnSparseF(const char* name, UInt_t entries){
559  // generate new THnSparseF, axes are defined in GetDimParams()
560 
561  Int_t count = 0;
562  UInt_t tmp = entries;
563  while(tmp!=0){
564  count++;
565  tmp = tmp &~ -tmp; // clear lowest bit
566  }
567 
568  TString hnTitle(name);
569  const Int_t dim = count;
570  Int_t nbins[dim];
571  Double_t xmin[dim];
572  Double_t xmax[dim];
573 
574  Int_t i=0;
575  Int_t c=0;
576  while(c<dim && i<32){
577  if(entries&(1<<i)){
578 
579  TString label("");
580  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
581  hnTitle += Form(";%s",label.Data());
582  c++;
583  }
584 
585  i++;
586  }
587  hnTitle += ";";
588 
589  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
590 }
591 
592 //________________________________________________________________________
594 {
595  // stores label and binning of axis for THnSparse
596 
597  const Double_t pi = TMath::Pi();
598 
599  switch(iEntry){
600 
601  case 0:
602  label = "V0 centrality (%)";
603  nbins = 10;
604  xmin = 0.;
605  xmax = 100.;
606  break;
607 
608  case 1:
609  label = "corrected jet pt";
610  nbins = 20;
611  xmin = 0.;
612  xmax = 200.;
613  break;
614 
615  case 2:
616  if(fDoWiderTrackBin) {
617  label = "track pT";
618  nbins = 40;
619  xmin = 0.;
620  xmax = 10.;
621  } else {
622  label = "track pT";
623  nbins = 100;
624  xmin = 0.;
625  xmax = 10;
626  }
627  break;
628 
629  case 3:
630  label = "deltaEta";
631  nbins = 24;
632  xmin = -1.2;
633  xmax = 1.2;
634  break;
635 
636  case 4:
637  label = "deltaPhi";
638  nbins = 72;
639  xmin = -0.5*pi;
640  xmax = 1.5*pi;
641  break;
642 
643  case 5:
644  label = "leading jet";
645  nbins = 3;
646  xmin = -0.5;
647  xmax = 2.5;
648  break;
649 
650  case 6:
651  label = "trigger track";
652  nbins =10;
653  xmin = 0;
654  xmax = 50;
655  break;
656 
657  case 7:
658  label = "deltaR";
659  nbins = 10;
660  xmin = 0.;
661  xmax = 5.0;
662  break;
663 
664  case 8:
665  label = "leading track";
666  nbins = 13;
667  xmin = 0;
668  xmax = 50;
669  break;
670  }
671 }
672 
673 //_________________________________________________
674 // From CF event mixing code PhiCorrelations
676 {
677  // clones a track list by using AliBasicTrack which uses much less memory (used for event mixing)
678  TObjArray* tracksClone = new TObjArray;
679  tracksClone->SetOwner(kTRUE);
680 
681  // Loop over all tracks
682  AliBasicParticle * clone = 0;
683  AliTrackContainer * tracks = GetTrackContainer("tracksForCorrelations");
684 
685  for (auto particle : tracks->accepted())
686  {
687  // Fill some QA information about the tracks
688  Int_t trackPtBin = GetTrackPtBin(particle->Pt());
689  if(trackPtBin > -1) fHistTrackEtaPhi[trackPtBin]->Fill(particle->Eta(),particle->Phi());
690 
691  // Create new particle
692  clone = new AliBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
693  // Set so that we can do comparisons using the IsEqual() function.
694  clone ->SetUniqueID(particle->GetUniqueID());
695 
696  tracksClone->Add(clone);
697  }
698 
699  return tracksClone;
700 }
701 
713  return EffCorrection(trackETA, trackPT, fBeamType);
714 }
715 
731  // default (current) parameters
732  // x-variable = track pt, y-variable = track eta
733  Double_t x = trackPT;
734  Double_t y = trackETA;
735  Double_t TRefficiency = -999;
736  Int_t runNUM = fCurrentRunNumber;
737  Int_t runSwitchGood = -999;
738  //Int_t centbin = -99;
739 
740  Double_t etaaxis = 0;
741  Double_t ptaxis = 0;
742 
743  Int_t effSwitch = fDoEffCorrection;
744 
745  if (beamType != AliAnalysisTaskEmcal::kpp) {
746  if(effSwitch == 1) {
747  // Semi-Good OROC C08 Runlists
748  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;
749 
750  // Good Runlists
751  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;
752 
753  // Determine which efficiency to use.
754  // This is just a way to map all possible values of the cent bin and runSwitchGood to a unique flag.
755  // 4 is the number of cent bins, and we want to index the effSwitch starting at 2.
756  if (runSwitchGood != -999) {
757  effSwitch = 2 + runSwitchGood*4 + fCentBin;
758  }
759  }
760 
761  // set up a switch for different parameter values...
762  switch(effSwitch) {
763  case 1 :
764  // first switch value - TRefficiency not used so = 1
765  // In this case, the run number isn't in any run list, so efficiency = 1
766  TRefficiency = 1.0;
767  break;
768 
769  case 2 :
770  // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
771  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);
772  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])));
773  TRefficiency = ptaxis*etaaxis;
774  break;
775 
776  case 3 :
777  // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
778  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);
779  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])));
780  TRefficiency = ptaxis*etaaxis;
781  break;
782 
783  case 4 :
784  // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
785  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);
786  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])));
787  TRefficiency = ptaxis*etaaxis;
788  break;
789 
790  case 5 :
791  // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
792  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);
793  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])));
794  TRefficiency = ptaxis*etaaxis;
795  break;
796 
797  case 6 :
798  // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
799  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);
800  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])));
801  TRefficiency = ptaxis*etaaxis;
802  break;
803 
804  case 7 :
805  // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
806  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);
807  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])));
808  TRefficiency = ptaxis*etaaxis;
809  break;
810 
811  case 8 :
812  // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
813  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);
814  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])));
815  TRefficiency = ptaxis*etaaxis;
816  break;
817 
818  case 9 :
819  // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
820  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);
821  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])));
822  TRefficiency = ptaxis*etaaxis;
823  break;
824 
825  default :
826  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
827  // ie. The efficiency correction is disabled.
828  TRefficiency = 1.0;
829  }
830  }
831  else {
832  // Track efficiency for pp
833  // Calculated using LHC12f1a. See analysis note for more details!
834 
835  if (fDoEffCorrection != 0) {
836  // If the trackPt > 6 GeV, then all we need is this coefficient
837  Double_t coefficient = 0.898052; // p6
838  if (trackPT < 6) {
839  coefficient = (1 + -0.442232 * trackPT // p0
840  + 0.501831 * std::pow(trackPT, 2) // p1
841  + -0.252024 * std::pow(trackPT, 3) // p2
842  + 0.062964 * std::pow(trackPT, 4) // p3
843  + -0.007681 * std::pow(trackPT, 5) // p4
844  + 0.000365 * std::pow(trackPT, 6)); // p5
845  }
846 
847  // Calculate track eff
848  TRefficiency = coefficient * (1 + 0.402825 * std::abs(trackETA) // p7
849  + -2.213152 * std::pow(trackETA, 2) // p8
850  + 4.311098 * std::abs(std::pow(trackETA, 3)) // p9
851  + -2.778200 * std::pow(trackETA, 4)); // p10
852  }
853  else {
854  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
855  TRefficiency = 1;
856  }
857  }
858 
859  return TRefficiency;
860 }
861 
862 //________________________________________________________________________
863 void AliAnalysisTaskEmcalJetHMEC::FillHist(TH1 * hist, Double_t fillValue, Double_t weight, Bool_t noCorrection)
864 {
865  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
866  {
867  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
868  hist->Fill(fillValue, weight);
869  }
870  else
871  {
872  // Determine where to get the values in the correction hist
873  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(fillValue);
874 
875  std::vector <Double_t> yBinsContent;
876  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), fillValue));
877  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
878  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
879 
880  // Loop over all possible bins to contribute.
881  // If content is 0 then calling Fill won't make a difference
882  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
883  {
884  // Don't bother trying to fill in the weight is 0
885  if (yBinsContent.at(index-1) > 0) {
886  // Determine the value to fill based on the center of the bins.
887  // This in principle allows the binning between the correction and hist to be different
888  Double_t fillLocation = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
889  AliDebug(4, TString::Format("fillLocation: %f, weight: %f", fillLocation, yBinsContent.at(index-1)));
890  // minus 1 since loop starts at 1
891  hist->Fill(fillLocation, weight*yBinsContent.at(index-1));
892  }
893  }
894 
895  //TEMP
896  //hist->Draw();
897  //END TEMP
898  }
899 }
900 
901 //________________________________________________________________________
902 void AliAnalysisTaskEmcalJetHMEC::FillHist(THnSparse * hist, Double_t *fillValue, Double_t weight, Bool_t noCorrection)
903 {
904  if (fJESCorrectionHist == 0 || noCorrection == kTRUE)
905  {
906  AliDebugStream(3) << GetName() << ":" << hist->GetName() << ": " << std::boolalpha << "Using normal weights: JESHist: " << (fJESCorrectionHist ? fJESCorrectionHist->GetName() : "Null") << ", noCorrection: " << noCorrection << std::endl;
907  hist->Fill(fillValue, weight);
908  }
909  else
910  {
911  // Jet pt is always located in the second position
912  Double_t jetPt = fillValue[1];
913 
914  // Determine where to get the values in the correction hist
915  Int_t xBin = fJESCorrectionHist->GetXaxis()->FindBin(jetPt);
916 
917  std::vector <Double_t> yBinsContent;
918  AliDebug(3, TString::Format("%s: Attempt to access weights from JES correction hist %s with jet pt %f!", GetName(), hist->GetName(), jetPt));
919  AccessSetOfYBinValues(fJESCorrectionHist, xBin, yBinsContent);
920  AliDebug(3, TString::Format("weights size: %zd", yBinsContent.size()));
921 
922  // Loop over all possible bins to contribute.
923  // If content is 0 then calling Fill won't make a difference
924  for (Int_t index = 1; index <= fJESCorrectionHist->GetYaxis()->GetNbins(); index++)
925  {
926  // Don't bother trying to fill in the weight is 0
927  if (yBinsContent.at(index-1) > 0) {
928  // Determine the value to fill based on the center of the bins.
929  // This in principle allows the binning between the correction and hist to be different
930  fillValue[1] = fJESCorrectionHist->GetYaxis()->GetBinCenter(index);
931  AliDebug(4,TString::Format("fillValue[1]: %f, weight: %f", fillValue[1], yBinsContent.at(index-1)));
932  // minus 1 since loop starts at 1
933  hist->Fill(fillValue, weight*yBinsContent.at(index-1));
934  }
935  }
936  }
937 }
938 
939 //________________________________________________________________________
940 void AliAnalysisTaskEmcalJetHMEC::AccessSetOfYBinValues(TH2D * hist, Int_t xBin, std::vector <Double_t> & yBinsContent, Double_t scaleFactor)
941 {
942  for (Int_t index = 1; index <= hist->GetYaxis()->GetNbins(); index++)
943  {
944  //yBinsContent[index-1] = hist->GetBinContent(hist->GetBin(xBin,index));
945  yBinsContent.push_back(hist->GetBinContent(hist->GetBin(xBin,index)));
946 
947  if (scaleFactor >= 0)
948  {
949  // -1 since index starts at 1
950  hist->SetBinContent(hist->GetBin(xBin,index), yBinsContent.at(index-1)/scaleFactor);
951  }
952  }
953 }
954 
972 {
973  // Initialize grid connection if necessary
974  if (filename.Contains("alien://") && !gGrid) {
975  TGrid::Connect("alien://");
976  }
977 
978  // Setup hist name if a track or cluster bias was defined.
979  // NOTE: This can always be disabled by setting kDisableBias.
980  // We arbitrarily add 0.1 to test since the values are doubles and cannot be
981  // tested directly for equality. If we are still less than disable bins, then
982  // it has been set and we should format it.
983  // NOTE: To ensure we can disable, we don't just take the member values!
984  // NOTE: The histBaseName will be attempted if the formatted name cannot be found.
985  TString histBaseName = histName;
986  if (trackBias + 0.1 < AliAnalysisTaskEmcalJetHMEC::kDisableBias) {
987  histName = TString::Format("%s_Track%.2f", histName.Data(), trackBias);
988  }
989  if (clusterBias + 0.1 < AliAnalysisTaskEmcalJetHMEC::kDisableBias) {
990  histName = TString::Format("%s_Clus%.2f", histName.Data(), clusterBias);
991  }
992 
993  // Open file containing the correction
994  TFile * jesCorrectionFile = TFile::Open(filename);
995  if (!jesCorrectionFile || jesCorrectionFile->IsZombie()) {
996  AliError(TString::Format("%s: Could not open JES correction file %s", GetName(), filename.Data()));
997  return kFALSE;
998  }
999 
1000  // Retrieve the histogram containing the correction and safely add it to the task.
1001  TH2D * JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histName.Data()));
1002  if (JESCorrectionHist) {
1003  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histName.Data(), filename.Data()));
1004  }
1005  else {
1006  AliError(TString::Format("%s: JES correction hist name \"%s\" not found in file %s.", GetName(), histName.Data(), filename.Data()));
1007 
1008  // Attempt the base name instead of the formatted hist name
1009  JESCorrectionHist = dynamic_cast<TH2D*>(jesCorrectionFile->Get(histBaseName.Data()));
1010  if (JESCorrectionHist) {
1011  AliInfo(TString::Format("%s: JES correction hist name \"%s\" loaded from file %s.", GetName(), histBaseName.Data(), filename.Data()));
1012  histName = histBaseName;
1013  }
1014  else
1015  {
1016  AliError(TString::Format("%s: JES correction with base hist name %s not found in file %s.", GetName(), histBaseName.Data(), filename.Data()));
1017  return kFALSE;
1018  }
1019  }
1020 
1021  // Clone to ensure that the hist is available
1022  TH2D * tempHist = static_cast<TH2D *>(JESCorrectionHist->Clone());
1023  tempHist->SetDirectory(0);
1024  SetJESCorrectionHist(tempHist);
1025 
1026  // Close file
1027  jesCorrectionFile->Close();
1028 
1029  // Append to task name for clarity
1030  // Unfortunately, this doesn't change the name of the output list (it would need to be
1031  // changed in the AnalysisManager output container), so the suffix is still important
1032  // if this correction is manually configured!
1033  TString tempName = GetName();
1034  TString tag = "_JESCorr";
1035  // Append the tag if it isn't already included
1036  if (tempName.Index(tag) == -1) {
1037  // Insert before the suffix
1038  Ssiz_t suffixLocation = tempName.Last('_');
1039  tempName.Insert(suffixLocation, tag.Data());
1040 
1041  // Set the new name
1042  AliDebug(3, TString::Format("%s: Setting task name to %s", GetName(), tempName.Data()));
1043  SetName(tempName.Data());
1044  }
1045 
1046  // Successful
1047  return kTRUE;
1048 }
1049 
1055  const char *nTracks,
1056  const char *nCaloClusters,
1057  // Jet options
1058  const Double_t trackBias,
1059  const Double_t clusterBias,
1060  const Double_t minJetArea,
1061  // Mixed event options
1062  const Int_t nTracksMixedEvent, // Additionally acts as a switch for enabling mixed events
1063  const Int_t minNTracksMixedEvent,
1064  const Int_t minNEventsMixedEvent,
1065  const UInt_t nCentBinsMixedEvent,
1066  // Triggers
1067  UInt_t trigEvent,
1068  UInt_t mixEvent,
1069  // Options
1070  const Double_t trackEta,
1071  const Bool_t lessSparseAxes,
1072  const Bool_t widerTrackBin,
1073  // Corrections
1074  const Int_t doEffCorrSW,
1075  const Bool_t JESCorrection,
1076  const char * JESCorrectionFilename,
1077  const char * JESCorrectionHistName,
1078  const char *suffix
1079  )
1080 {
1081  // Get the pointer to the existing analysis manager via the static access method.
1082  //==============================================================================
1083  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1084  if (!mgr)
1085  {
1086  AliErrorClass("No analysis manager to connect to.");
1087  return NULL;
1088  }
1089 
1090  // Check the analysis type using the event handlers connected to the analysis manager.
1091  //==============================================================================
1092  AliVEventHandler* handler = mgr->GetInputEventHandler();
1093  if (!handler)
1094  {
1095  AliErrorClass("This task requires an input event handler");
1096  return NULL;
1097  }
1098 
1099  //-------------------------------------------------------
1100  // Init the task and do settings
1101  //-------------------------------------------------------
1102 
1103  // Determine data type
1104  enum EDataType_t {
1105  kUnknown,
1106  kESD,
1107  kAOD
1108  };
1109 
1110  EDataType_t dataType = kUnknown;
1111 
1112  if (handler->InheritsFrom("AliESDInputHandler")) {
1113  dataType = kESD;
1114  }
1115  else if (handler->InheritsFrom("AliAODInputHandler")) {
1116  dataType = kAOD;
1117  }
1118 
1119  // Determine cluster and track names
1120  TString trackName(nTracks);
1121  TString clusName(nCaloClusters);
1122 
1123  if (trackName == "usedefault") {
1124  if (dataType == kESD) {
1125  trackName = "Tracks";
1126  }
1127  else if (dataType == kAOD) {
1128  trackName = "tracks";
1129  }
1130  else {
1131  trackName = "";
1132  }
1133  }
1134 
1135  if (clusName == "usedefault") {
1136  if (dataType == kESD) {
1137  clusName = "CaloClusters";
1138  }
1139  else if (dataType == kAOD) {
1140  clusName = "caloClusters";
1141  }
1142  else {
1143  clusName = "";
1144  }
1145  }
1146 
1147  TString name("AliAnalysisTaskJetH");
1148  if (!trackName.IsNull()) {
1149  name += TString::Format("_%s", trackName.Data());
1150  }
1151  if (!clusName.IsNull()) {
1152  name += TString::Format("_%s", clusName.Data());
1153  }
1154  if (strcmp(suffix, "") != 0) {
1155  name += TString::Format("_%s", suffix);
1156  }
1157 
1158  AliAnalysisTaskEmcalJetHMEC *correlationTask = new AliAnalysisTaskEmcalJetHMEC(name);
1159  // Set jet bias
1160  correlationTask->SetTrackBias(trackBias);
1161  correlationTask->SetClusterBias(clusterBias);
1162  // Mixed events
1163  correlationTask->SetEventMixing(static_cast<Bool_t>(nTracksMixedEvent));
1164  correlationTask->SetNumberOfMixingTracks(nTracksMixedEvent);
1165  correlationTask->SetMinNTracksForMixedEvents(minNTracksMixedEvent);
1166  correlationTask->SetMinNEventsForMixedEvents(minNEventsMixedEvent);
1167  correlationTask->SetNCentBinsMixedEvent(nCentBinsMixedEvent);
1168  // Triggers
1169  correlationTask->SetTriggerType(trigEvent);
1170  correlationTask->SetMixedEventTriggerType(mixEvent);
1171  // Options
1172  correlationTask->SetNCentBins(5);
1173  correlationTask->SetVzRange(-10,10);
1174  correlationTask->SetDoLessSparseAxes(lessSparseAxes);
1175  correlationTask->SetDoWiderTrackBin(widerTrackBin);
1176  // Corrections
1177  correlationTask->SetDoEffCorr(doEffCorrSW);
1178  if (JESCorrection == kTRUE)
1179  {
1180  Bool_t result = correlationTask->RetrieveAndInitializeJESCorrectionHist(JESCorrectionFilename, JESCorrectionHistName, correlationTask->GetTrackBias(), correlationTask->GetClusterBias());
1181  if (!result) {
1182  AliErrorClass("Failed to successfully retrieve and initialize the JES correction! Task initialization continuing without JES correction (can be set manually later).");
1183  }
1184  }
1185 
1186  // Jet parameters determined by how we ran the jet finder
1187  Double_t jetRadius = 0.2;
1188  Double_t minClusterPt = 3;
1189  Double_t minTrackPt = 3;
1190 
1191  // Add Containers
1192  // Clusters
1193  AliClusterContainer * clusterContainer = correlationTask->AddClusterContainer(clusName);
1194  clusterContainer->SetMinE(minClusterPt);
1195 
1196  // Tracks
1197  // For jet finding
1198  AliTrackContainer * tracksForJets = new AliTrackContainer(trackName);
1199  tracksForJets->SetName("tracksForJets");
1200  tracksForJets->SetMinPt(minTrackPt);
1201  tracksForJets->SetEtaLimits(-1.0*trackEta, trackEta);
1202  // Adopt the container
1203  correlationTask->AdoptParticleContainer(tracksForJets);
1204  // For correlations
1205  AliTrackContainer * tracksForCorrelations = new AliTrackContainer(trackName);
1206  tracksForCorrelations->SetName("tracksForCorrelations");
1207  tracksForCorrelations->SetMinPt(0.15);
1208  tracksForCorrelations->SetEtaLimits(-1.0*trackEta, trackEta);
1209  // Adopt the container
1210  correlationTask->AdoptParticleContainer(tracksForCorrelations);
1211 
1212  // Jets
1213  AliJetContainer * jetContainer = correlationTask->AddJetContainer(AliJetContainer::kFullJet,
1216  jetRadius,
1218  tracksForJets,
1219  clusterContainer);
1220  jetContainer->SetJetAreaCut(minJetArea);
1221  jetContainer->SetMaxTrackPt(100);
1222  jetContainer->SetJetPtCut(0.1);
1223 
1224  //-------------------------------------------------------
1225  // Final settings, pass to manager and set the containers
1226  //-------------------------------------------------------
1227 
1228  mgr->AddTask(correlationTask);
1229 
1230  // Create containers for input/output
1231  mgr->ConnectInput (correlationTask, 0, mgr->GetCommonInputContainer() );
1232  AliAnalysisDataContainer * cojeth = mgr->CreateContainer(correlationTask->GetName(),
1233  TList::Class(),
1234  AliAnalysisManager::kOutputContainer,
1235  Form("%s", AliAnalysisManager::GetCommonFileName()));
1236  mgr->ConnectOutput(correlationTask, 1, cojeth);
1237 
1238  return correlationTask;
1239 }
1240 
Jet-hadron correlations analysis task for central Pb-Pb and pp.
const char * filename
Definition: TestFCM.C:1
virtual void SetTriggerType(UInt_t te)
Set the trigger event trigger selection.
void GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector &particleOne, AliVParticle *particleTwo, Double_t &deltaEta, Double_t &deltaPhi, Double_t &deltaR)
Bool_t fDoEventMixing
flag to do evt mixing
double Double_t
Definition: External.C:58
Definition: External.C:260
Int_t fMinNEventsMixedEvents
threshold to use event pool # events
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
UInt_t fNCentBinsMixedEvent
N cent bins for the event mixing pool.
Definition: External.C:236
Int_t GetNTracks() const
TH2D * fJESCorrectionHist
Histogram containing the jet energy scale correction.
AliJetContainer * GetJetContainer(Int_t i=0) const
! Number of elements in track pt binned arrays
void AdoptParticleContainer(AliParticleContainer *cont)
! Number of elements in mixed event multiplicity binned arrays
static AliAnalysisTaskEmcalJetHMEC * AddTaskEmcalJetHMEC(const char *nTracks="usedefault", const char *nCaloClusters="usedefault", const Double_t trackBias=5, const Double_t clusterBias=5, const Double_t minJetArea=0.4, 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 Double_t trackEta=0.9, const Bool_t lessSparseAxes=0, const Bool_t widerTrackBin=0, 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")
TH1 * fHistTrackPt
! Track pt spectrum
const AliTrackIterableContainer accepted() const
Container with name, TClonesArray and cuts for particles.
! Number of elements in jet pt binned arrays
Declaration of class AliTLorentzVector.
Bool_t RetrieveAndInitializeJESCorrectionHist(TString filename, TString histName, Double_t trackBias=AliAnalysisTaskEmcalJetHMEC::kDisableBias, Double_t clusterBias=AliAnalysisTaskEmcalJetHMEC::kDisableBias)
static Double_t p50_90SG[17]
50-90% centrality semi-good runs
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
static Double_t p0_10SG[17]
0-10% centrality semi-good runs
TH2 * fHistJetHEtaPhi
! Eta-phi distribution of jets which are in jet-hadron correlations
TCanvas * c
Definition: TestFitELoss.C:172
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
AliEventPoolManager * fPoolMgr
! Event pool manager
Int_t fCentBin
!event centrality bin
TH2 * fHistJetHBias[6][5][3]
! Jet-hadron correlations of jets which meet the constituent bias criteria (the arrays correspond to ...
static Double_t p50_90G[17]
50-90% centrality good runs
void SetVzRange(Double_t min, Double_t max)
Int_t fDoEffCorrection
Control the efficiency correction. See EffCorrection() for meaning of values.
virtual Double_t Pt() const
Bool_t fDoWiderTrackBin
True if the track pt bins in the THnSparse should be wider.
virtual Double_t Eta() const
AliClusterContainer * AddClusterContainer(const char *n)
Int_t fMinNTracksMixedEvents
threshold to use event pool # tracks
virtual THnSparse * NewTHnSparseF(const char *name, UInt_t entries)
! Number of elements in eta binned arrays
! Arbitrarily large value which can be used to disable the constituent bias. Can be used for either t...
TH1 * fHistJetPtBias[6]
! Jet pt spectrum of jets which meet the constituent bias criteria (the array corresponds to centrali...
BeamType
Switch for the beam type.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
static Double_t p30_50SG[17]
30-50% centrality semi-good runs
TH1 * fHistJetPt[6]
! Jet pt spectrum (the array corresponds to centrality bins)
AliEmcalJet * GetLeadingJet(const char *opt="")
static Double_t p10_30G[17]
10-30% centrality good runs
static Double_t p0_10G[17]
0-10% centrality good runs
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
unsigned int UInt_t
Definition: External.C:33
TH3 * fHistJHPsi
! Psi angle distribution
Bool_t fNoMixedEventJESCorrection
True if the jet energy scale correction should be applied to mixed event histograms.
virtual void SetMinNTracksForMixedEvents(Int_t nmt)
void FillHist(TH1 *hist, Double_t fillValue, Double_t weight=1.0, Bool_t noCorrection=kFALSE)
UInt_t fTriggerType
Event selection for jets (ie triggered events).
BeamType fForceBeamType
forced beam type
Definition: External.C:228
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:148
BeamType fBeamType
!event beam type
static Double_t p30_50G[17]
30-50% centrality good runs
Double_t fCent
!event centrality
Double_t fClusterBias
Jet cluster bias.
virtual void SetClusterBias(Double_t b)
Require a cluster with pt > b in jet.
virtual void SetNCentBins(Int_t n)
virtual void SetEventMixing(Bool_t enable)
static Double_t p10_30SG[17]
10-30% centrality semi-good runs
Bool_t fDisableFastPartition
True if task should be disabled for the fast partition, where the EMCal is not included.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliEmcalList * fOutput
!output list
TH2 * fHistJetEtaPhi
! Jet eta-phi distribution
virtual void SetNumberOfMixingTracks(Int_t tracks)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t fVertex[3]
!event vertex
Double_t EffCorrection(Double_t trkETA, Double_t trkPT, AliAnalysisTaskEmcal::BeamType beamType) const
AliTrackContainer * GetTrackContainer(Int_t i=0) const
TH1 * fHistEventRejection
!book keep reasons for rejecting event
void SetMakeGeneralHistograms(Bool_t g)
! Number of elements in centrality binned arrays
THnSparse * fhnMixedEvents
! Mixed events THnSparse
void AccessSetOfYBinValues(TH2D *hist, Int_t xBin, std::vector< Double_t > &yBinsContent, Double_t scaleFactor=-1.0)
Base task in the EMCAL jet framework.
Double_t MaxClusterPt() const
Definition: AliEmcalJet.h:147
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
virtual void SetMixedEventTriggerType(UInt_t me)
Set the mixed event trigger selection.
Bool_t fDoLessSparseAxes
True if there should be fewer THnSparse axes.
const AliTrackIterableMomentumContainer accepted_momentum() const
const char Option_t
Definition: External.C:48
virtual void SetMinNEventsForMixedEvents(Int_t nme)
const Double_t pi
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Int_t fNMixingTracks
size of track buffer for event mixing
EDataType_t
Switch for the data type.
UInt_t fMixingEventType
Event selection for mixed events.
void SetMaxTrackPt(Float_t b)
virtual Double_t Phi() const
Container structure for EMCAL clusters.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:63
Container for jet within the EMCAL jet framework.
Definition: External.C:196
virtual void SetNCentBinsMixedEvent(Bool_t centbins)
virtual void SetTrackBias(Double_t b)
Require a track with pt > b in jet.
void SetJetAreaCut(Float_t cut)
TH2 * fHistJetH[6][5][3]
! Jet-hadron correlations (the arrays correspond to centrality, jet pt bins, and eta bins) ...
THnSparse * fhnJH
! JetH THnSparse
TH2 * fHistTrackEtaPhi[7]
! Track eta-phi distribution (the array corresponds to track pt)
virtual void GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)