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