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