AliPhysics  aaf9c62 (aaf9c62)
AliAnalysisTaskGenerateThermalBackgroundMC.cxx
Go to the documentation of this file.
1 /**********************************************************************************
2 * Copyright (C) 2016, Copyright Holders of the ALICE Collaboration *
3 * All rights reserved. *
4 * *
5 * Redistribution and use in source and binary forms, with or without *
6 * modification, are permitted provided that the following conditions are met: *
7 * * Redistributions of source code must retain the above copyright *
8 * notice, this list of conditions and the following disclaimer. *
9 * * Redistributions in binary form must reproduce the above copyright *
10 * notice, this list of conditions and the following disclaimer in the *
11 * documentation and/or other materials provided with the distribution. *
12 * * Neither the name of the <organization> nor the *
13 * names of its contributors may be used to endorse or promote products *
14 * derived from this software without specific prior written permission. *
15 * *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19 * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26 * *********************************************************************************/
27 
29 
30 #include <TClonesArray.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TF1.h>
34 
35 #include <AliAnalysisManager.h>
36 #include <AliVEventHandler.h>
37 #include <AliTLorentzVector.h>
38 #include <AliBasicParticle.h>
39 
43 
49  fEventInitialized(false),
50  fEvent(nullptr),
51  fOutputCollectionName(""),
52  fThermalParticlesArray(),
53  fChargedParticleFraction(0.67),
54  fRandom(0),
55  fPtDistribution(nullptr),
56  fAlpha(2),
57  fBeta(0.3),
58  fMinChargedPt(0.15),
59  fMinNeutralPt(0.3),
60  fMaxPt(100.),
61  fUseGaussianForN(true),
62  fNGaussianMean(3500),
63  fNGaussianSigma(500),
64  fMinN(3000),
65  fMaxN(4500),
66  fMinEta(-0.9),
67  fMaxEta(0.9),
68  fOutput(0),
69  fHistManager()
70 {}
71 
77  AliAnalysisTaskSE(name),
78  fEventInitialized(false),
79  fEvent(nullptr),
83  fRandom(0),
85  fAlpha(2),
86  fBeta(0.3),
87  fMinChargedPt(0.15),
88  fMinNeutralPt(0.3),
89  fMaxPt(100.),
90  fUseGaussianForN(true),
91  fNGaussianMean(3500),
92  fNGaussianSigma(500),
93  fMinN(3000),
94  fMaxN(4500),
95  fMinEta(-0.9),
96  fMaxEta(0.9),
97  fOutput(0),
98  fHistManager(name)
99 {
100  DefineOutput(1, TList::Class());
101 }
102 
107 
112 {
113  // Define gamma distribution f_gamma to sample single particle pT.
114  // We use the convention in TMath::GammaDist(x, alpha, mu, beta). We fix mu = 0.
115  // Then f_gamma(x;a,b) = x^(a-1) * (1/b)^a * e^(-x/b) ~ x^(a-1) * e^(-x/b).
116  fPtDistribution = new TF1("fGamma", "TMath::GammaDist(x, [0], 0, [1])", TMath::Min(fMinChargedPt, fMinNeutralPt), fMaxPt);
117  fPtDistribution->SetParameter(0, fAlpha);
118  fPtDistribution->SetParameter(1, fBeta);
119 
121  // Allocate histograms
122 
123  // Number of particles per event vs. mean pT per event
124  TString histname = "hNvsMeanPt";
125  TString title = histname + ";#it{N}_{particles};<#it{p}_{T}>";
126  if (fUseGaussianForN) {
127  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, fNGaussianMean-5*fNGaussianSigma, fNGaussianMean+5*fNGaussianSigma, 200, 0., 4*fBeta);
128  }
129  else {
130  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, fMinN*0.9, fMaxN*1.1, 200, 0., 4*fBeta);
131  }
132 
133  // Eta-phi of particles
134  histname = "hEtaPhi";
135  title = histname + ";#eta;#phi";
136  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, -1., 1., 100, 0, 2*TMath::Pi());
137 
138  // Particle pT spectrum
139  histname = "hPt";
140  title = histname + ";#it{p}_{T}";
141  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 20);
142 
143  // Plot highest-pT particle in the event
144  histname = "hLeadingPt";
145  title = histname + ";#it{p}_{T,leading}";
146  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 20);
147 
148  // Background density
149  histname = "hBackgroundDensity";
150  title = histname + ";#rho (GeV/#it{c});counts";
151  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 500);
152 
153  // DeltaPt
154  histname = "hDeltaPt02";
155  title = histname + ";#delta#it{p}_{T} (GeV/#it{c});counts";
156  fHistManager.CreateTH1(histname.Data(), title.Data(), 400, -100, 100);
157 
158  histname = "hDeltaPt04";
159  title = histname + ";#delta#it{p}_{T} (GeV/#it{c});counts";
160  fHistManager.CreateTH1(histname.Data(), title.Data(), 400, -100, 100);
161 
163 
164  // Allow for output files
165  OpenFile(1);
166  fOutput = new TList();
167  fOutput->SetOwner();
168 
169  TIter next(fHistManager.GetListOfHistograms());
170  TObject* obj = 0;
171  while ((obj = next())) {
172  fOutput->Add(obj);
173  }
174 
175  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
176 }
177 
182 {
183  fEvent = InputEvent();
184  if (!fEvent) {
185  AliError("Could not retrieve event! Returning!");
186  return;
187  }
188 
189  // Initialize random number generator to random seed
190  fRandom.SetSeed(0);
191  gRandom = &fRandom; // It seems this is necessary in order to use the TF1::GetRandom() function...
192 
193  // Create a new collection during the first event
195 
196  fEventInitialized = true;
197 }
198 
203 {
204  // Check to ensure that we are not trying to create a new branch on top of the old branch
205  TObject* existingObject = fEvent->FindListObject(fOutputCollectionName.c_str());
206  if (existingObject) {
207  AliFatal(TString::Format("Attempted to create a new branch, \"%s\", with the same name as an existing branch!", fOutputCollectionName.c_str()));
208  }
209 
210  // Create new branch and add it to the event. It will then be there for every event, and cleared for each event.
211  fThermalParticlesArray = new TClonesArray("AliBasicParticle", fMaxN+1);
213  fEvent->AddObject(fThermalParticlesArray);
214 }
215 
220 {
221  // Initialize the event if not initialized
222  if (!fEventInitialized) {
223  ExecOnce();
224  }
225 
226  // Only continue if we are initialized successfully
227  if (!fEventInitialized) {
228  return;
229  }
230 
231  // Generate the thermal particles, and write them to the event
232  Run();
233 
234  // Plot some properties of the generated thermal particles
235  FillHistograms();
236 
237 }
238 
243 {
244  // Get the TClonesArray
245  fThermalParticlesArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fOutputCollectionName.c_str()));
246  if (!fThermalParticlesArray) {
247  AliError(Form("%s: Could not retrieve array with name %s!", GetName(), fOutputCollectionName.c_str()));
248  return;
249  }
250 
251  // Set number of particles in event
252  Int_t N = 0;
253  if (fUseGaussianForN) {
255  }
256  else {
257  N = fRandom.Uniform(fMinN, fMaxN);
258  }
259 
260  // Loop through particles to:
261  // (1) Sample their pT from the thermal distribution, and assign a random eta/phi
262  // (2) Create an AliBasicParticle, and write it to the TClonesArray
263  for (Int_t i=0; i<N; i++) {
264 
265  Double_t eta = fRandom.Uniform(fMinEta, fMaxEta);
266  Double_t phi = fRandom.Uniform(0., 2*TMath::Pi());
267  Double_t pT = fPtDistribution->GetRandom();
268 
269  Int_t charge = 1;
270  if (fRandom.Rndm() < fChargedParticleFraction) {
271  if (pT < fMinChargedPt) {
272  continue;
273  }
274  }
275  else {
276  charge = 0;
277  if (pT < fMinNeutralPt) {
278  continue;
279  }
280  }
281 
282  new((*fThermalParticlesArray)[i]) AliBasicParticle(eta, phi, pT, charge);
283 
284  }
285  fThermalParticlesArray->Compress();
286 }
287 
292 {
293  // Initialize the various sums to 0
294  Double_t ptSum = 0;
295  Double_t ptSumRC02 = 0;
296  Double_t ptSumRC04 = 0;
297  Double_t maxPt = 0;
298 
299  // Generate random cone eta-phi
300  Double_t jetR02 = 0.2;
301  Double_t jetR04 = 0.4;
302  Double_t etaRC02 = fRandom.Uniform(fMinEta + jetR02, fMaxEta - jetR02);
303  Double_t phiRC02 = fRandom.Uniform(jetR02, 2*TMath::Pi() - jetR02);
304  Double_t etaRC04 = fRandom.Uniform(fMinEta + jetR04, fMaxEta - jetR04);
305  Double_t phiRC04 = fRandom.Uniform(jetR04, 2*TMath::Pi() - jetR04);
306 
307  // Loop over particles.
308  Int_t N = fThermalParticlesArray->GetEntries();
309  for (Int_t i=0; i<N; i++) {
310 
311  AliBasicParticle* thermalParticle = dynamic_cast<AliBasicParticle*>(fThermalParticlesArray->At(i));
312  Double_t pT = thermalParticle->Pt();
313  Double_t eta = thermalParticle->Eta();
314  Double_t phi = thermalParticle->Phi();
315 
316  TString histname = "hPt";
317  fHistManager.FillTH1(histname, pT);
318 
319  histname = "hEtaPhi";
320  fHistManager.FillTH2(histname, eta, phi);
321 
322  ptSum += pT;
323 
324  if (pT > maxPt) {
325  maxPt = pT;
326  }
327 
328  if (GetDeltaR(eta, phi, etaRC02, phiRC02) < jetR02) {
329  ptSumRC02 += pT;
330  }
331  if (GetDeltaR(eta, phi, etaRC04, phiRC04) < jetR04) {
332  ptSumRC04 += pT;
333  }
334  }
335 
336  // Compute mean pT
337  Double_t meanPt = ptSum / N;
338  TString histname = "hNvsMeanPt";
339  fHistManager.FillTH2(histname, N, meanPt);
340 
341  // Compute leading particle pT in the event
342  histname = "hLeadingPt";
343  fHistManager.FillTH1(histname, maxPt);
344 
345  // Compute the background density
346  Double_t area = (fMaxEta - fMinEta) * 2 * TMath::Pi();
347  Double_t backgroundDensity = ptSum / area;
348  histname = "hBackgroundDensity";
349  fHistManager.FillTH1(histname, backgroundDensity);
350 
351  // Compute delta pT
352  Double_t deltaPt02 = ptSumRC02 - backgroundDensity * TMath::Pi() * jetR02 * jetR02;
353  histname = "hDeltaPt02";
354  fHistManager.FillTH1(histname, deltaPt02);
355 
356  Double_t deltaPt04 = ptSumRC04 - backgroundDensity * TMath::Pi() * jetR04 * jetR04;
357  histname = "hDeltaPt04";
358  fHistManager.FillTH1(histname, deltaPt04);
359 
360 }
361 
366 {
367  Double_t deltaPhi = TMath::Abs(phi - phiRef);
368  Double_t deltaEta = TMath::Abs(eta - etaRef);
369  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
370  return deltaR;
371 }
372 
377  const char *outputCollectionName,
378  const Double_t beta,
379  const char *suffix)
380 {
381  // Get the pointer to the existing analysis manager via the static access method.
382  //==============================================================================
383  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
384  if (!mgr)
385  {
386  ::Error("AddTaskGenerateThermalBackgroundMC", "No analysis manager to connect to.");
387  return 0;
388  }
389 
390  // Check the analysis type using the event handlers connected to the analysis manager.
391  //==============================================================================
392  AliVEventHandler* handler = mgr->GetInputEventHandler();
393  if (!handler)
394  {
395  ::Error("AddTaskGenerateThermalBackgroundMC", "This task requires an input event handler");
396  return 0;
397  }
398 
399  enum EDataType_t {
400  kUnknown,
401  kESD,
402  kAOD
403  };
404 
405  EDataType_t dataType = kUnknown;
406 
407  if (handler->InheritsFrom("AliESDInputHandler")) {
408  dataType = kESD;
409  }
410  else if (handler->InheritsFrom("AliAODInputHandler")) {
411  dataType = kAOD;
412  }
413 
414  //-------------------------------------------------------
415  // Init the task and do settings
416  //-------------------------------------------------------
417  TString name("AliAnalysisTaskGenerateThermalBackgroundMC");
418  if (strcmp(outputCollectionName,"") != 0) {
419  name += "_";
420  name += outputCollectionName;
421  }
422  if (strcmp(suffix,"") != 0) {
423  name += "_";
424  name += suffix;
425  }
426 
428  // Configure the task
430 
431  task->SetOutputCollectionName(outputCollectionName);
432  task->SetBeta(beta);
433 
434  //-------------------------------------------------------
435  // Final settings, pass to manager and set the containers
436  //-------------------------------------------------------
437 
438  mgr->AddTask(task);
439 
440  // Create containers for input/output
441  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
442  TString contname(name);
443  contname += "_histos";
444  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
445  TList::Class(),AliAnalysisManager::kOutputContainer,
446  Form("%s", AliAnalysisManager::GetCommonFileName()));
447  mgr->ConnectInput (task, 0, cinput1 );
448  mgr->ConnectOutput (task, 1, coutput1 );
449 
450  return task;
451 }
Int_t charge
Double_t fBeta
Value of b in the gamma distribution f_gamma ~ x^(a-1) * e^(-x/b)
double Double_t
Definition: External.C:58
std::string fOutputCollectionName
Name of TClonesArray output the thermal particles to the event.
const char * title
Definition: MakeQAPdf.C:27
static AliAnalysisTaskGenerateThermalBackgroundMC * AddTaskGenerateThermalBackgroundMC(const char *outputName="thermalparticles", const Double_t beta=0.3, const char *suffix="")
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
virtual Double_t Pt() const
bool fUseGaussianForN
Use a Gaussian for number of particles per event. Otherwise, use flat distribution.
TRandom * gRandom
virtual Double_t Eta() const
bool fEventInitialized
If the event is initialized properly.
Double_t fMinChargedPt
Min pT of charged thermal particles.
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
int Int_t
Definition: External.C:63
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t GetDeltaR(Double_t eta, Double_t phi, Double_t etaRef, Double_t phiRef)
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Int_t fMinN
Min number of particles in thermal model (for flat distribution)
This task uses a thermal model of background to generate a collection of truth-level background parti...
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t fAlpha
Value of a in the gamma distribution f_gamma ~ x^(a-1) * e^(-x/b)
TClonesArray * fThermalParticlesArray
! Thermal particle collection
const char Option_t
Definition: External.C:48
Double_t fMinNeutralPt
Min pT of neutral thermal particles.
Int_t fMaxN
Max number of particles in thermal model (for flat distribution)
Double_t fChargedParticleFraction
Fraction of thermal particles set to have nonzero charge.
TF1 * fPtDistribution
! Distribution to sample particle pT from
virtual Double_t Phi() const
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65