AliPhysics  9fe175b (9fe175b)
 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 
14 #include "AliAnalysisManager.h"
15 #include "AliInputEventHandler.h"
16 #include "AliEventPoolManager.h"
17 #include "AliBasicParticle.h"
18 #include "AliVTrack.h"
19 #include "AliEmcalJet.h"
20 #include "AliTLorentzVector.h"
21 
22 #include "AliClusterContainer.h"
23 #include "AliTrackContainer.h"
24 
26 
27 //________________________________________________________________________
30  fTrackBias(5),
31  fClusterBias(5),
32  fDoEventMixing(kFALSE),
33  fNMixingTracks(50000), fMinNTracksMixedEvents(5000), fMinNEventsMixedEvents(5), fNCentBinsMixedEvent(10),
34  fPoolMgr(0),
35  fTriggerType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
36  fDoEffCorrection(0), fEffFunctionCorrection(0),
37  fEmbeddingCorrectionHist(0),
38  fDoLessSparseAxes(0), fDoWiderTrackBin(0),
39  fHistTrackPt(0),
40  fHistJetEtaPhi(0),
41  fHistJHPsi(0),
42  fHistJetHEtaPhi(0),
43  fhnMixedEvents(0),
44  fhnJH(0)
45 {
46  // Default Constructor
47  InitializeArraysToZero();
48 }
49 
50 //________________________________________________________________________
52  AliAnalysisTaskEmcalJet(name, kTRUE),
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  fDoEffCorrection(0), fEffFunctionCorrection(0),
60  fEmbeddingCorrectionHist(0),
61  fDoLessSparseAxes(0), fDoWiderTrackBin(0),
62  fHistTrackPt(0),
63  fHistJetEtaPhi(0),
64  fHistJHPsi(0),
65  fHistJetHEtaPhi(0),
66  fhnMixedEvents(0),
67  fhnJH(0)
68 {
69  // Constructor
71 }
72 
73 //________________________________________________________________________
75 {
76  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; trackPtBin++){
77  fHistTrackEtaPhi[trackPtBin]=0;
78  }
79  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
80  fHistJetPt[centralityBin]=0;
81  fHistJetPtBias[centralityBin]=0;
82  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
83  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
84  fHistJetH[centralityBin][jetPtBin][etaBin]=0;
85  fHistJetHBias[centralityBin][jetPtBin][etaBin]=0;
86  }
87  }
88  }
89 }
90 
91 //________________________________________________________________________
93  // Called once
95  OpenFile(1);
96 
97  // Create histograms
98  fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
99  fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
100  fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
101 
102  fHistJHPsi = new TH3F("fHistJHPsi","Jet-Hadron ntr-trpt-dpsi",20,0,100,200,0,20,120,0,180);
103 
104  fOutput->Add(fHistTrackPt);
105  fOutput->Add(fHistJetEtaPhi);
106  fOutput->Add(fHistJetHEtaPhi);
107  fOutput->Add(fHistJHPsi);
108 
109  TString name;
110 
111  for(Int_t trackPtBin = 0; trackPtBin < kMaxTrackPtBins; ++trackPtBin){
112  name = Form("fHistTrackEtaPhi_%i", trackPtBin);
113  fHistTrackEtaPhi[trackPtBin] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
114  fOutput->Add(fHistTrackEtaPhi[trackPtBin]);
115  }
116 
117  for(Int_t centralityBin = 0; centralityBin < kMaxCentralityBins; ++centralityBin){
118  name = Form("fHistJetPt_%i",centralityBin);
119  fHistJetPt[centralityBin] = new TH1F(name,name,200,0,200);
120  fOutput->Add(fHistJetPt[centralityBin]);
121 
122  name = Form("fHistJetPtBias_%i",centralityBin);
123  fHistJetPtBias[centralityBin] = new TH1F(name,name,200,0,200);
124  fOutput->Add(fHistJetPtBias[centralityBin]);
125 
126  for(Int_t jetPtBin = 0; jetPtBin < kMaxJetPtBins; ++jetPtBin){
127  for(Int_t etaBin = 0; etaBin < kMaxEtaBins; ++etaBin){
128  name = Form("fHistJetH_%i_%i_%i",centralityBin,jetPtBin,etaBin);
129  fHistJetH[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
130  fOutput->Add(fHistJetH[centralityBin][jetPtBin][etaBin]);
131 
132  name = Form("fHistJetHBias_%i_%i_%i",centralityBin,jetPtBin,etaBin);
133  fHistJetHBias[centralityBin][jetPtBin][etaBin]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
134  fOutput->Add(fHistJetHBias[centralityBin][jetPtBin][etaBin]);
135  }
136  }
137  }
138 
139  UInt_t cifras = 0; // bit coded, see GetDimParams() below
140  if(fDoLessSparseAxes) {
141  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
142  } else {
143  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
144  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
145  }
146  fhnJH = NewTHnSparseF("fhnJH", cifras);
147  fhnJH->Sumw2();
148  fOutput->Add(fhnJH);
149 
150  if(fDoEventMixing){
151  if(fDoLessSparseAxes) {
152  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
153  } else {
154  cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7;
155  //cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
156  }
157  fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
158  fhnMixedEvents->Sumw2();
159  fOutput->Add(fhnMixedEvents);
160  }
161 
162  PostData(1, fOutput);
163 
164  // Event Mixing
165  Int_t poolSize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
166  // ZVertex
167  Int_t nZVertexBins = 10;
168  Double_t* zVertexBins = GenerateFixedBinArray(nZVertexBins, -10, 10);
169  // Event activity (centrality of multiplicity)
170  Int_t nEventActivityBins = 8;
171  Double_t* eventActivityBins = 0;
172  // +1 to accomodate the fact that we define bins rather than array entries.
173  Double_t multiplicityBins[kMixedEventMulitplictyBins+1] = {0., 4., 9., 15., 25., 35., 55., 100., 500.};
174 
175  if (fForceBeamType != kpp ) { //all besides pp
176  // Event Activity is centrality in AA, pA
177  nEventActivityBins = fNCentBinsMixedEvent;
178  eventActivityBins = GenerateFixedBinArray(nEventActivityBins, 0, 100);
179  }
180  else if (fForceBeamType == kpp) { //for pp only
181  // Event Activity is multiplicity in pp
182  eventActivityBins = multiplicityBins;
183  }
184 
185  fPoolMgr = new AliEventPoolManager(poolSize, fNMixingTracks, nEventActivityBins, eventActivityBins, nZVertexBins, zVertexBins);
186 }
187 
188 //________________________________________________________________________
189 Int_t AliAnalysisTaskEmcalJetHMEC::GetEtaBin(Double_t eta) const {
190  // Get eta bin for histos.
191 
192  Int_t etabin = -1;
193  eta = TMath::Abs(eta);
194  if (eta <= 0.4) etabin = 0;
195  else if (eta > 0.4 && eta < 0.8) etabin = 1;
196  else if (eta >= 0.8) etabin = 2;
197  return etabin;
198 }
199 
200 //________________________________________________________________________
202 {
203  Int_t ptBin = -1;
204  if (pt < 0.5) ptBin = 0;
205  else if (pt < 1 ) ptBin = 1;
206  else if (pt < 2 ) ptBin = 2;
207  else if (pt < 3 ) ptBin = 3;
208  else if (pt < 5 ) ptBin = 4;
209  else if (pt < 8 ) ptBin = 5;
210  else if (pt < 20 ) ptBin = 6;
211 
212  return ptBin;
213 }
214 
215 //________________________________________________________________________
217 {
218  // Get jet pt bin for histos.
219 
220  Int_t ptBin = -1;
221  if (pt >= 15 && pt < 20) ptBin = 0;
222  else if (pt >= 20 && pt < 25) ptBin = 1;
223  else if (pt >= 25 && pt < 30) ptBin = 2;
224  else if (pt >= 30 && pt < 60) ptBin = 3;
225  else if (pt >= 60) ptBin = 4;
226 
227  return ptBin;
228 }
229 
230 //________________________________________________________________________
233 }
234 
235 //________________________________________________________________________
237  // Main loop called for each event
238 
239  // Retrieve clusters
241  if (!clusters) {
242  AliError(Form("%s: Unable to retrieve clusters!", GetName()));
243  return kTRUE;
244  }
245 
246  // Retrieve tracks
247  AliTrackContainer * tracks = static_cast<AliTrackContainer * >(GetParticleContainer("tracksForCorrelations"));
248  if (!tracks) {
249  AliError(Form("%s: Unable to retrieve tracks!", GetName()));
250  return kTRUE;
251  }
252 
253  // Retrieve jets
254  AliJetContainer * jets = GetJetContainer(0);
255  if (!jets) {
256  AliError(Form("%s: Unable to retrieve jets!", GetName()));
257  return kTRUE;
258  }
259 
260 
261  // Used to calculate the angle betwene the jet and the hadron
262  TVector3 jetVector;
263  // Get z vertex
264  Double_t zVertex=fVertex[2];
265  // Flags
266  Bool_t biasedJet = kFALSE;
267  Bool_t leadJet = kFALSE;
268  // Relative angles and distances
269  Double_t deltaPhi = 0;
270  Double_t deltaEta = 0;
271  Double_t deltaR = 0;
272  // Event activity (centrality or multipilicity)
273  Double_t eventActivity = 0;
274  // Efficiency correction
275  Double_t efficiency = -999;
276  // Determining bins for histogram indices
277  Int_t jetPtBin = -1;
278  Int_t etaBin = -1;
279  // For comparison to the current jet
280  AliEmcalJet * leadingJet = jets->GetLeadingJet();
281  // For getting the proper properties of tracks
282  AliTLorentzVector track;
283 
284  // Determine the trigger for the current event
285  UInt_t eventTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
286 
287  for (auto jet : jets->accepted()) {
288  // Selects only events that we are interested in (ie triggered)
289  if (!(eventTrigger & fTriggerType)) continue;
290 
291  // Jet properties
292  // Determine if we have the lead jet
293  leadJet = kFALSE;
294  if (jet == leadingJet) leadJet = kTRUE;
295 
296  // Determine if the jet is biased
297  biasedJet = BiasedJet(jet);
298 
299  // Calculate vector
300  jetVector.SetXYZ(jet->Px(), jet->Py(), jet->Pz());
301 
302  // Fill jet properties
303  FillHist(fHistJetPt[fCentBin], jet->Pt());
304  if (biasedJet == kTRUE) {
305  FillHist(fHistJetPtBias[fCentBin], jet->Pt());
306  }
307 
308  fHistJetEtaPhi->Fill(jet->Eta(), jet->Phi());
309 
310  if (jet->Pt() > 15) {
311 
312  for (auto trackIter : tracks->accepted_momentum()) {
313 
314  // Get proper track proeprties
315  track.Clear();
316  track = trackIter.first;
317 
318  // Determine relative angles and distances and set the respective variables
319  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
320 
321  // Determine bins for filling histograms
322  // jet Pt
323  jetPtBin = GetJetPtBin(jet->Pt());
324  if (jetPtBin < 0)
325  {
326  AliError(Form("Jet Pt Bin negative: %f", jet->Pt()));
327  continue;
328  }
329  // eta
330  etaBin = GetEtaBin(deltaEta);
331  if (etaBin < 0) {
332  AliError(Form("Eta Bin negative: %f", deltaEta));
333  continue;
334  }
335 
336  // Fill track properties
337  fHistTrackPt->Fill(track.Pt());
338 
339  if ( (jet->Pt() > 20.) && (jet->Pt() < 60.) ) {
340  fHistJHPsi->Fill(tracks->GetNTracks(), track.Pt(), track.Vect().Angle(jetVector) * TMath::RadToDeg() );
341  }
342 
343  fHistJetH[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
344  fHistJetHEtaPhi->Fill(deltaEta, deltaPhi);
345 
346  // Calculate single particle tracking efficiency for correlations
347  efficiency = EffCorrection(track.Eta(), track.Pt(), fDoEffCorrection);
348 
349  if (biasedJet == kTRUE) {
350  fHistJetHBias[fCentBin][jetPtBin][etaBin]->Fill(deltaPhi, track.Pt());
351 
352  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
353  eventActivity = fCent;
354  }
355  else if (fBeamType == kpp) {
356  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
357  }
358 
359  if(fDoLessSparseAxes) { // check if we want all dimensions
360  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
361  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
362  } else {
363  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
364  FillHist(fhnJH, triggerEntries, 1.0/efficiency);
365  }
366  }
367 
368  } //track loop
369  }//jet pt cut
370  }//jet loop
371 
372  //Prepare to do event mixing
373 
374  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
375  TObjArray* tracksClone = 0;
376 
377  if(fDoEventMixing == kTRUE){
378 
379  // event mixing
380 
381  // 1. First get an event pool corresponding in mult (cent) and
382  // zvertex to the current event. Once initialized, the pool
383  // should contain nMix (reduced) events. This routine does not
384  // pre-scan the chain. The first several events of every chain
385  // will be skipped until the needed pools are filled to the
386  // specified depth. If the pool categories are not too rare, this
387  // should not be a problem. If they are rare, you could lose
388  // statistics.
389 
390  // 2. Collect the whole pool's content of tracks into one TObjArray
391  // (bgTracks), which is effectively a single background super-event.
392 
393  // 3. The reduced and bgTracks arrays must both be passed into
394  // FillCorrelations(). Also nMix should be passed in, so a weight
395  // of 1./nMix can be applied.
396 
397  AliEventPool *pool = 0;
398  if (fBeamType == kAA || fBeamType == kpA) {//everything but pp
399  pool = fPoolMgr->GetEventPool(fCent, zVertex);
400  }
401  else if (fBeamType == kpp) {//pp only
402  pool = fPoolMgr->GetEventPool(static_cast<Double_t>(tracks->GetNTracks()), zVertex);
403  }
404 
405  if (!pool){
406  if (fBeamType == kAA || fBeamType == kpA) AliFatal(Form("No pool found for centrality = %f, zVertex = %f", fCent, zVertex));
407  else if (fBeamType == kpp) AliFatal(Form("No pool found for ntracks_pp = %d, zVertex = %f", tracks->GetNTracks(), zVertex));
408  return kTRUE;
409  }
410 
411  // The number of events in the pool
412  Int_t nMix = pool->GetCurrentNEvents();
413 
414  if(eventTrigger & fTriggerType) {
415  // check for a trigger jet
416  if (pool->IsReady() || pool->NTracksInPool() >= fMinNTracksMixedEvents || nMix >= fMinNEventsMixedEvents) {
417 
418  for (auto jet : jets->accepted()) {
419 
420  // Jet properties
421  // Determine if we have the lead jet
422  leadJet = kFALSE;
423  if (jet == leadingJet) { leadJet = kTRUE; }
424 
425  // Determine if the jet is biased
426  biasedJet = BiasedJet(jet);
427 
428  // Make sure event contains jet above our threshold (reduce stats of sparse)
429  if (jet->Pt() < 15) continue;
430 
431  // Fill for biased jet triggers only
432  if (biasedJet == kTRUE) {
433 
434  // Fill mixed-event histos here
435  for (Int_t jMix=0; jMix < nMix; jMix++) {
436  TObjArray* bgTracks = pool->GetEvent(jMix);
437 
438  for(Int_t ibg=0; ibg < bgTracks->GetEntries(); ibg++){
439 
440  AliBasicParticle *bgTrack = static_cast<AliBasicParticle*>(bgTracks->At(ibg));
441  if(!bgTrack)
442  {
443  AliError(Form("%s:Failed to retrieve tracks from mixed events", GetName()));
444  }
445 
446  // Fill into TLorentzVector for use with functions below
447  track.Clear();
448  track.SetPtEtaPhiE(bgTrack->Pt(), bgTrack->Eta(), bgTrack->Phi(), 0);
449 
450  // Calculate single particle tracking efficiency of mixed events for correlations
451  efficiency = EffCorrection(track.Eta(), track.Pt(), fDoEffCorrection);
452 
453  // Phi is [-0.5*TMath::Pi(), 3*TMath::Pi()/2.]
454  GetDeltaEtaDeltaPhiDeltaR(track, jet, deltaEta, deltaPhi, deltaR);
455 
456  if (fBeamType == kAA || fBeamType == kpA) { //pA and AA
457  eventActivity = fCent;
458  }
459  else if (fBeamType == kpp) {
460  eventActivity = static_cast<Double_t>(tracks->GetNTracks());
461  }
462 
463  if(fDoLessSparseAxes) { // check if we want all the axis filled
464  Double_t triggerEntries[6] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet)};
465  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), kTRUE);
466  } else {
467  Double_t triggerEntries[7] = {eventActivity, jet->Pt(), track.Pt(), deltaEta, deltaPhi, static_cast<Double_t>(leadJet), deltaR};
468  FillHist(fhnMixedEvents, triggerEntries, 1./(nMix*efficiency), kTRUE);
469  }
470  }
471  }
472  }
473  }
474  }
475  }
476 
477  if(eventTrigger & fMixingEventType) {
478  tracksClone = CloneAndReduceTrackList();
479 
480  //update pool if jet in event or not
481  pool->UpdatePool(tracksClone);
482  }
483 
484  } // end of event mixing
485 
486  return kTRUE;
487 }
488 
489 //________________________________________________________________________
491 {
492  //just terminate
493 }
494 
495 //________________________________________________________________________
497 {
498  if ((jet->MaxTrackPt() > fTrackBias) || (jet->MaxClusterPt() > fClusterBias))
499  {
500  return kTRUE;
501  }
502  return kFALSE;
503 }
504 
505 //________________________________________________________________________
506 void AliAnalysisTaskEmcalJetHMEC::GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector & particleOne, AliVParticle * particleTwo, Double_t & deltaEta, Double_t & deltaPhi, Double_t & deltaR)
507 {
508  // TODO: Understand order of arguments to DeltaPhi vs DeltaEta
509  // Returns deltaPhi in symmetric range so that we can calculate DeltaR.
510  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -1.0*TMath::Pi(), TMath::Pi());
511  deltaEta = particleOne.Eta() - particleTwo->Eta();
512  deltaR = TMath::Sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
513 
514  // Adjust to the normal range after the DeltaR caluclation
515  deltaPhi = DeltaPhi(particleTwo->Phi(), particleOne.Phi(), -0.5*TMath::Pi(), 3*TMath::Pi()/2.);
516 }
517 
518 //________________________________________________________________________
519 THnSparse* AliAnalysisTaskEmcalJetHMEC::NewTHnSparseF(const char* name, UInt_t entries){
520  // generate new THnSparseF, axes are defined in GetDimParams()
521 
522  Int_t count = 0;
523  UInt_t tmp = entries;
524  while(tmp!=0){
525  count++;
526  tmp = tmp &~ -tmp; // clear lowest bit
527  }
528 
529  TString hnTitle(name);
530  const Int_t dim = count;
531  Int_t nbins[dim];
532  Double_t xmin[dim];
533  Double_t xmax[dim];
534 
535  Int_t i=0;
536  Int_t c=0;
537  while(c<dim && i<32){
538  if(entries&(1<<i)){
539 
540  TString label("");
541  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
542  hnTitle += Form(";%s",label.Data());
543  c++;
544  }
545 
546  i++;
547  }
548  hnTitle += ";";
549 
550  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
551 }
552 
553 //________________________________________________________________________
554 void AliAnalysisTaskEmcalJetHMEC::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
555 {
556  // stores label and binning of axis for THnSparse
557 
558  const Double_t pi = TMath::Pi();
559 
560  switch(iEntry){
561 
562  case 0:
563  label = "V0 centrality (%)";
564  nbins = 10;
565  xmin = 0.;
566  xmax = 100.;
567  break;
568 
569  case 1:
570  label = "corrected jet pt";
571  nbins = 20;
572  xmin = 0.;
573  xmax = 200.;
574  break;
575 
576  case 2:
577  if(fDoWiderTrackBin) {
578  label = "track pT";
579  nbins = 40;
580  xmin = 0.;
581  xmax = 10.;
582  } else {
583  label = "track pT";
584  nbins = 100;
585  xmin = 0.;
586  xmax = 10;
587  }
588  break;
589 
590  case 3:
591  label = "deltaEta";
592  nbins = 24;
593  xmin = -1.2;
594  xmax = 1.2;
595  break;
596 
597  case 4:
598  label = "deltaPhi";
599  nbins = 72;
600  xmin = -0.5*pi;
601  xmax = 1.5*pi;
602  break;
603 
604  case 5:
605  label = "leading jet";
606  nbins = 3;
607  xmin = -0.5;
608  xmax = 2.5;
609  break;
610 
611  case 6:
612  label = "trigger track";
613  nbins =10;
614  xmin = 0;
615  xmax = 50;
616  break;
617 
618  case 7:
619  label = "deltaR";
620  nbins = 10;
621  xmin = 0.;
622  xmax = 5.0;
623  break;
624 
625  case 8:
626  label = "leading track";
627  nbins = 13;
628  xmin = 0;
629  xmax = 50;
630  break;
631  }
632 }
633 
634 //_________________________________________________
635 // From CF event mixing code PhiCorrelations
637 {
638  // clones a track list by using AliBasicTrack which uses much less memory (used for event mixing)
639  TObjArray* tracksClone = new TObjArray;
640  tracksClone->SetOwner(kTRUE);
641 
642  // Loop over all tracks
643  AliBasicParticle * clone = 0;
644  AliTrackContainer * tracks = GetTrackContainer("tracksForCorrelations");
645 
646  for (auto particle : tracks->accepted())
647  {
648  // Fill some QA information about the tracks
649  Int_t trackPtBin = GetTrackPtBin(particle->Pt());
650  if(trackPtBin > -1) fHistTrackEtaPhi[trackPtBin]->Fill(particle->Eta(),particle->Phi());
651 
652  // Create new particle
653  clone = new AliBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
654  // Set so that we can do comparisons using the IsEqual() function.
655  clone ->SetUniqueID(particle->GetUniqueID());
656 
657  tracksClone->Add(clone);
658  }
659 
660  return tracksClone;
661 }
662 
663 //________________________________________________________________________
664 Double_t AliAnalysisTaskEmcalJetHMEC::EffCorrection(Double_t trackETA, Double_t trackPT, Int_t effSwitch) const {
665  // default (current) parameters
666  // x-variable = track pt, y-variable = track eta
667  Double_t x = trackPT;
668  Double_t y = trackETA;
669  Double_t TRefficiency = -999;
670  Int_t runNUM = fCurrentRunNumber;
671  Int_t runSwitchGood = -999;
672  //Int_t centbin = -99;
673 
674  Double_t etaaxis = 0;
675  Double_t ptaxis = 0;
676 
677  if(effSwitch < 1) {
678  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;
679 
680  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;
681 
682  // Determine which efficiency to use.
683  // This is just a way to map all possible values of the cent bin and runSwitchGood to a unique flag.
684  // 4 is the number of cent bins, and we want to index the effSwitch starting at 2.
685  effSwitch = 2 + runSwitchGood*4 + fCentBin;
686  }
687 
688  // 0-10% centrality: Semi-Good Runs
689  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};
690  // 10-30% centrality: Semi-Good Runs
691  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};
692  // 30-50% centrality: Semi-Good Runs
693  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};
694  // 50-90% centrality: Semi-Good Runs
695  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};
696 
697  // 0-10% centrality: Good Runs
698  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};
699  // 10-30% centrality: Good Runs
700  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};
701  // 30-50% centrality: Good Runs
702  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};
703  // 50-90% centrality: Good Runs
704  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};
705 
706  // set up a switch for different parameter values...
707  switch(effSwitch) {
708  case 1 :
709  // first switch value - TRefficiency not used so = 1
710  TRefficiency = 1.0;
711  break;
712 
713  case 2 :
714  // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
715  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);
716  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])));
717  TRefficiency = ptaxis*etaaxis;
718  break;
719 
720  case 3 :
721  // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
722  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);
723  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])));
724  TRefficiency = ptaxis*etaaxis;
725  break;
726 
727  case 4 :
728  // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
729  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);
730  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])));
731  TRefficiency = ptaxis*etaaxis;
732  break;
733 
734  case 5 :
735  // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
736  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);
737  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])));
738  TRefficiency = ptaxis*etaaxis;
739  break;
740 
741  case 6 :
742  // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
743  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);
744  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])));
745  TRefficiency = ptaxis*etaaxis;
746  break;
747 
748  case 7 :
749  // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
750  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);
751  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])));
752  TRefficiency = ptaxis*etaaxis;
753  break;
754 
755  case 8 :
756  // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
757  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);
758  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])));
759  TRefficiency = ptaxis*etaaxis;
760  break;
761 
762  case 9 :
763  // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
764  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);
765  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])));
766  TRefficiency = ptaxis*etaaxis;
767  break;
768 
769  default :
770  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
771  TRefficiency = 1.0;
772 
773  }
774 
775  return TRefficiency;
776 }
777 
778 //________________________________________________________________________
779 void AliAnalysisTaskEmcalJetHMEC::FillHist(TH1 * hist, Double_t fillValue, Double_t weight, Bool_t noCorrection)
780 {
781  if (fEmbeddingCorrectionHist == 0 || noCorrection == kTRUE)
782  {
783  hist->Fill(fillValue, weight);
784  }
785  else
786  {
787  // Determine where to get the values in the correction hist
788  Int_t xBin = fEmbeddingCorrectionHist->GetXaxis()->FindBin(fillValue);
789 
790  std::vector <Double_t> yBinsContent;
791  accessSetOfYBinValues(fEmbeddingCorrectionHist, xBin, yBinsContent);
792 
793  // Loop over all possible bins to contribute.
794  // If content is 0 then calling Fill won't make a difference
795  for (Int_t index = 1; index <= fEmbeddingCorrectionHist->GetYaxis()->GetNbins(); index++)
796  {
797  // Determine the value to fill based on the center of the bins.
798  // This in principle allows the binning between the correction and hist to be different
799  Double_t fillLocation = fEmbeddingCorrectionHist->GetYaxis()->GetBinCenter(index);
800  Printf("fillLocation: %f, weight: %f", fillLocation, yBinsContent.at(index-1));
801  // minus 1 since loop starts at 1
802  hist->Fill(fillLocation, weight*yBinsContent.at(index-1));
803  }
804 
805  //TEMP
806  //hist->Draw();
807  //END TEMP
808  }
809 }
810 
811 //________________________________________________________________________
812 void AliAnalysisTaskEmcalJetHMEC::FillHist(THnSparse * hist, Double_t *fillValue, Double_t weight, Bool_t noCorrection)
813 {
814  if (fEmbeddingCorrectionHist == 0 || noCorrection == kTRUE)
815  {
816  hist->Fill(fillValue, weight);
817  }
818  else
819  {
820  // Jet pt is always located in the second position
821  Double_t jetPt = fillValue[1];
822 
823  // Determine where to get the values in the correction hist
824  Int_t xBin = fEmbeddingCorrectionHist->GetXaxis()->FindBin(jetPt);
825 
826  std::vector <Double_t> yBinsContent;
827  accessSetOfYBinValues(fEmbeddingCorrectionHist, xBin, yBinsContent);
828 
829  // Loop over all possible bins to contribute.
830  // If content is 0 then calling Fill won't make a difference
831  for (Int_t index = 1; index <= fEmbeddingCorrectionHist->GetYaxis()->GetNbins(); index++)
832  {
833  // Determine the value to fill based on the center of the bins.
834  // This in principle allows the binning between the correction and hist to be different
835  fillValue[1] = fEmbeddingCorrectionHist->GetYaxis()->GetBinCenter(index);
836  Printf("fillValue[1]: %f, weight: %f", fillValue[1], yBinsContent.at(index-1));
837  // minus 1 since loop starts at 1
838  hist->Fill(fillValue, weight*yBinsContent.at(index-1));
839  }
840  }
841 }
842 
843 //________________________________________________________________________
844 void AliAnalysisTaskEmcalJetHMEC::accessSetOfYBinValues(TH2F * hist, Int_t xBin, std::vector <Double_t> & yBinsContent, Double_t scaleFactor)
845 {
846  for (Int_t index = 1; index <= hist->GetYaxis()->GetNbins(); index++)
847  {
848  //yBinsContent[index-1] = hist->GetBinContent(hist->GetBin(xBin,index));
849  yBinsContent.push_back(hist->GetBinContent(hist->GetBin(xBin,index)));
850 
851  if (scaleFactor >= 0)
852  {
853  // -1 since index starts at 1
854  hist->SetBinContent(hist->GetBin(xBin,index), yBinsContent.at(index-1)/scaleFactor);
855  }
856  }
857 }
858 
void accessSetOfYBinValues(TH2F *hist, Int_t xBin, std::vector< Double_t > &yBinsContent, Double_t scaleFactor=-1.0)
void GetDeltaEtaDeltaPhiDeltaR(AliTLorentzVector &particleOne, AliVParticle *particleTwo, Double_t &deltaEta, Double_t &deltaPhi, Double_t &deltaR)
Bool_t fDoEventMixing
flag to do evt mixing
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.
Int_t GetNTracks() const
AliJetContainer * GetJetContainer(Int_t i=0) const
virtual Int_t GetTrackPtBin(Double_t pt) const
const AliTrackIterableContainer accepted() const
Container with name, TClonesArray and cuts for particles.
Declaration of class AliTLorentzVector.
AliEventPoolManager * fPoolMgr
! Event pool manager
Int_t fCentBin
!event centrality bin
virtual Double_t Pt() const
virtual Int_t GetJetPtBin(Double_t pt) const
virtual Double_t Eta() const
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
AliEmcalJet * GetLeadingJet(const char *opt="")
TH3 * fHistJHPsi
! Psi angle distribution
BeamType fForceBeamType
forced beam type
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:139
virtual Int_t GetEtaBin(Double_t eta) const
BeamType fBeamType
!event beam type
void FillHist(TH1 *hist, Double_t fillValue, Double_t weight=1.0, Bool_t noCorrection=kTRUE)
Double_t fCent
!event centrality
Double_t fClusterBias
Jet cluster bias.
ClassImp(AliAnalysisTaskEmcalJetHMEC) AliAnalysisTaskEmcalJetHMEC
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliEmcalList * fOutput
!output list
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:138
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const AliTrackIterableMomentumContainer accepted_momentum() const
const Double_t pi
const Int_t nbins
const AliJetIterableContainer accepted() const
Int_t fNMixingTracks
size of track buffer for event mixing
virtual Double_t EffCorrection(Double_t trkETA, Double_t trkPT, Int_t effswitch) const
virtual Double_t Phi() const
Container structure for EMCAL clusters.
Container for jet within the EMCAL jet framework.
THnSparse * fhnJH
! JetH THnSparse
virtual void GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)