AliPhysics  master (3d17d9d)
AliMCSpectraWeights.cxx
Go to the documentation of this file.
1 
7 #include "AliMCEvent.h"
8 #include "AliMCSpectraWeights.h"
9 #include "AliStack.h"
10 #include "TF1.h"
11 #include "TFile.h"
12 #include "TH1D.h"
13 #include "TH3F.h"
14 #include "TObjArray.h"
15 #include "TParticle.h"
16 #include "TParticlePDG.h"
17 #include <string>
18 
19 #if defined(__CLING__)
20 #include <algorithm>
21 #include <array>
22 #include <memory>
23 #endif
24 
25 int GetBinFromTH3(TH3F* h, std::array<float, 3> const & _values) {
26  return h->FindBin(_values[0], _values[1], _values[2]);
27 }
28 
29 void FillTH3WithValue(TH3F* h, std::array<float, 3> _xyzValue, float _value) {
30  auto const _iBin = GetBinFromTH3(h, _xyzValue);
31  h->SetBinContent(_iBin, _value);
32  h->SetBinError(_iBin, 1e-30);
33 }
34 
39  : fstCollisionSystem("pp"), fstFileMCSpectra(""), fstFilePublished(""),
40  fstSavedObjName("fMCSpectraWeights"), fstSavedListName("dNdPt_ParCor"),
41  fstPartTypes(0), fstCentralities(0), fBinsMultCent{0}, fBinsPt{0},
43  fHistMCFractions(nullptr), fHistMCWeights(nullptr), fMCEvent(nullptr),
45  fbTaskStatus(AliMCSpectraWeights::TaskState::kAllEmpty),
46  fFlag(AliMCSpectraWeights::SysFlag::kNominal), fUseMultiplicity(kTRUE) {}
47 
56 AliMCSpectraWeights::AliMCSpectraWeights(std::string collisionSystem,
57  std::string stName,
60  fstSavedObjName("fMCSpectraWeights"), fstSavedListName("dNdPt_ParCor"),
63  fHistMCFractions(nullptr), fHistMCWeights(nullptr), fMCEvent(nullptr),
65  fbTaskStatus(AliMCSpectraWeights::TaskState::kAllEmpty), fFlag(flag),
66  fUseMultiplicity(kTRUE) {
67  fstCollisionSystem = collisionSystem;
68  std::for_each(
70  [](char& c) { c = std::tolower(static_cast<unsigned char>(c)); });
71  // setting uniform name
72  if (fstCollisionSystem == "pp" || fstCollisionSystem == "p-p")
73  fstCollisionSystem = "pp";
74  else if (fstCollisionSystem == "ppb" || fstCollisionSystem == "p-pb")
75  fstCollisionSystem = "ppb";
76  else if (fstCollisionSystem == "pbpb" || fstCollisionSystem == "pb-pb")
77  fstCollisionSystem = "pbpb";
78  else if (fstCollisionSystem == "xexe" || fstCollisionSystem == "xe-xe")
79  fstCollisionSystem = "xexe";
80  else
81  fstCollisionSystem = "pp";
82 
83  // set default Binning
84  // pT binning
85  fBinsPt = {0.0, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5,
86  0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
87  1.0, 1.1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4,
88  2.6, 2.8, 3.0, 3.2, 3.6, 4.0, 5.0, 6.0, 8.0,
89  10.0, 13.0, 20.0, 30.0, 50.0, 80.0, 100.0, 200.0};
90  // multiplicity binning
91  if (fstCollisionSystem.find("pp") != std::string::npos) {
92  fNCentralities = 10;
93  std::vector<std::string> tmpCent{};
94  tmpCent.reserve(fNCentralities);
95  for (int icent = 0; icent < fNCentralities; icent++) {
96  tmpCent.push_back(Form("%d", icent));
97  }
98  fstCentralities = tmpCent;
99  fBinsMultCent.reserve(51);
100  for (int i = 0; i < 51; ++i) {
101  fBinsMultCent.push_back(i);
102  }
103  } else if (fstCollisionSystem.find("ppb") != std::string::npos) {
104  fBinsMultCent.reserve(301);
105  for (int i = 0; i < 301; ++i) {
106  fBinsMultCent.push_back(i);
107  }
108  } else if (fstCollisionSystem.find("pbpb") != std::string::npos) {
109  fBinsMultCent.reserve(201);
110  for (int i = 0; i < 201; ++i) {
111  fBinsMultCent.push_back(i * 25);
112  }
113  } else if (fstCollisionSystem.find("xexe") != std::string::npos) {
114  fBinsMultCent.reserve(201);
115  for (int i = 0; i < 201; ++i) {
116  fBinsMultCent.push_back(i * 25);
117  }
118  } else {
119  fBinsMultCent.reserve(101);
120  for (int i = 0; i < 101; ++i) {
121  fBinsMultCent.push_back(i);
122  }
123  }
124  if (fstCollisionSystem.find("pp") == std::string::npos) {
125  std::vector<std::string> cent{"MB", "c0005", "c0510", "c0010",
126  "c1020", "c2040", "c4060", "c6080"};
127  fNCentralities = 8;
128  fstCentralities.clear();
129  fstCentralities = cent;
130  }
131  const std::vector<std::string> tmpPart{"Pion", "Proton", "Kaon",
132  "SigmaMinus", "SigmaPlus", "Rest"};
133  fstPartTypes = tmpPart;
134  fbTaskStatus = AliMCSpectraWeights::TaskState::kAllEmpty;
136  "alien:///alice/cern.ch/user/p/phuhn/AllPublishedFractions.root";
137 }
138 
140  // if(fHistMCGenPrimTrackParticle) delete fHistMCGenPrimTrackParticle; //
141  // is put into train output in AnalysisTask if(fHistDataFractions) delete
142  // fHistDataFractions; if(fHistMCFractions) delete fHistMCFractions;
143  // if(fHistMCWeights) delete fHistMCWeights;
144  if (fMCEvent)
145  fMCEvent = 0;
146 }
147 
148 void AliMCSpectraWeights::SetBinsPt(std::vector<double> bins) {
149  fBinsPt.clear();
150  fBinsPt.reserve(bins.size());
151  for (auto& x : bins) {
152  fBinsPt.push_back(x);
153  }
154 }
155 
156 void AliMCSpectraWeights::SetBinsMultCent(std::vector<double> bins) {
157  fBinsMultCent.clear();
158  fBinsMultCent.reserve(bins.size());
159  for (auto& x : bins) {
160  fBinsMultCent.push_back(x);
161  }
162 }
163 
173  // Histograms
175 
176  // Loading output from previous train
177  if (fstFileMCSpectra.length() > 5) // *.root
178  {
179  // printf("AliMCSpectraWeights:: Opening %s \n",
180  // fstFileMCSpectra.Data());
181  TFile* fInput = TFile::Open(fstFileMCSpectra.c_str());
182  if (fInput) {
183  if (fInput->GetNkeys() != 1) {
184  if (!fInput->Get("dNdPt_ParCor"))
185  std::cerr << "AliMCSpectraWeights::WARNING: more than 1 "
186  "list in the streamed file; please specify; "
187  "using 1st list;\n\n";
188  } else
189  fstSavedListName = fInput->GetListOfKeys()->At(0)->GetName();
190  // printf("AliMCSpectraWeights:: Loading %s from list %s\n",
191  // fstSavedObjName.Data(), fstSavedListName.Data());
192  TList* listMC = (TList*)fInput->Get(fstSavedListName.c_str());
193  if (!listMC) {
194  std::cerr << "AliMCSpectraWeights::ERROR: could not load list "
195  "in streamed "
196  "file\n";
197  } else {
198  auto tmpHist =
199  (TH3F*)listMC->FindObject("fHistMCGenPrimTrackParticle");
200  if (!tmpHist) {
201  std::cerr << "AliMCSpectraWeights::WARNING: Couln't get "
202  "fHistMCGenPrimTrackParticle\n";
203  // AliMCSpectraWeights* inWeights =
204  // (AliMCSpectraWeights*)listMC->FindObject(
205  // fstSavedObjName.c_str());
206  // if (!inWeights) {
207  // std::cerr <<
208  // "AliMCSpectraWeights::ERROR:
209  // Couln't load "
210  // "object from "
211  // "previous output";
212  // return;
213  // }
214  // if
215  // (AliMCSpectraWeights::LoadFromAliMCSpectraWeight(
216  // inWeights))
217  // fbTaskStatus =
218  // AliMCSpectraWeights::TaskState::kMCSpectraObtained;
219  // else
220  // return;
221  }
223  (TH3F*)tmpHist->Clone("fHistMCGenPrimTrackParticle_prev");
224  fHistMCGenPrimTrackParticle->SetDirectory(0);
225  if (fHistMCGenPrimTrackParticle->GetEntries() > 0) {
226  std::cerr << "Previous train has mc tracks\n";
227  fbTaskStatus =
228  AliMCSpectraWeights::TaskState::kMCSpectraObtained;
229  }
230  // delete tmpHist; // gets deleted with
231  // fInput->Close() delete listMC;
232  }
233  fInput->Close();
234  delete fInput;
235  }
236  }
237 
238  // Loading measured fractions
239  if (fbTaskStatus == AliMCSpectraWeights::TaskState::kMCSpectraObtained) {
241  fbTaskStatus = AliMCSpectraWeights::TaskState::kDataFractionLoaded;
242  }
243 
244  // Calculating weight factors
245  if (fbTaskStatus == AliMCSpectraWeights::TaskState::kDataFractionLoaded) {
247  fbTaskStatus = AliMCSpectraWeights::TaskState::kMCWeightCalculated;
248  }
249  }
250 
251  std::cerr << "AliMCSpectraWeights::INFO: Init finished with status "
252  << fbTaskStatus << std::endl;
253 }
254 
259  // Initalizing histograms
260  // histogram charged patricles pt:multcent:type
261  std::array<float, 7> partArray{};
262  std::array<float, 6> partArrayDATA{};
263  for (int i = 0; i < 7; ++i) {
264  partArray[i] = -0.5 + i;
265  if (i < 6)
266  partArrayDATA[i] = -0.5 + i;
267  }
268 
270  new TH3F("fHistMCGenPrimTrackParticle",
271  "histogram for charged particle composition;#it{p}_{T} "
272  "(GeV/#it{c});multiplicity or centrality;Particle type",
273  static_cast<int>(fBinsPt.size()) - 1,
274  static_cast<float*>(fBinsPt.data()),
275  static_cast<int>(fBinsMultCent.size()) - 1,
276  static_cast<float*>(fBinsMultCent.data()),
277  static_cast<int>(partArray.size()) - 1,
278  static_cast<float*>(partArray.data()));
280 
282  new TH3F("fHistDataFractions",
283  "DATA fractions histogram;#it{p}_{T} "
284  "(GeV/#it{c});multiplicity or centrality;Particle type",
285  static_cast<int>(fBinsPt.size()) - 1,
286  static_cast<float*>(fBinsPt.data()),
287  static_cast<int>(fBinsMultCent.size()) - 1,
288  static_cast<float*>(fBinsMultCent.data()),
289  static_cast<int>(partArrayDATA.size()) - 1,
290  static_cast<float*>(partArrayDATA.data()));
291  fHistDataFractions->Sumw2();
292 
294  new TH3F("fHistMCFractions",
295  "MC fractions histogram;#it{p}_{T} (GeV/#it{c});multiplicity "
296  "or centrality;Particle type",
297  static_cast<int>(fBinsPt.size()) - 1,
298  static_cast<float*>(fBinsPt.data()),
299  static_cast<int>(fBinsMultCent.size()) - 1,
300  static_cast<float*>(fBinsMultCent.data()),
301  static_cast<int>(partArray.size()) - 1,
302  static_cast<float*>(partArray.data()));
303  fHistMCFractions->Sumw2();
304 
305  fHistMCWeights = new TH3F(
306  "fHistMCWeights",
307  "MC weight histogram for charged particle composition;#it{p}_{T} "
308  "(GeV/#it{c});multiplicity or centrality;Particle type",
309  static_cast<int>(fBinsPt.size()) - 1,
310  static_cast<float*>(fBinsPt.data()),
311  static_cast<int>(fBinsMultCent.size()) - 1,
312  static_cast<float*>(fBinsMultCent.data()),
313  static_cast<int>(partArrayDATA.size()) - 1,
314  static_cast<float*>(partArrayDATA.data()));
315  fHistMCWeights->Sumw2();
316  // printf("AliMCSpectraWeights: init histos successful\n"); // works
317 }
318 
323  // TFile *fMeasuredFile = AliDataFile::OpenOADB(fstFilePublished.Data());
324  // TFile *fMeasuredFile = TFile::Open(fstFilePublished.c_str(), "OPEN");
325  auto fMeasuredFile = TFile::Open(fstFilePublished.c_str());
326  if (!fMeasuredFile) {
327  std::cerr << "AliMCSpectraWeights::Error: Could not load measured "
328  "fractions in "
329  << fstFilePublished << "\n";
330  return;
331  }
332 
333  for (auto& part : fstPartTypes) {
334  if (part.find("Rest") != std::string::npos ||
335  part.find("rest") != std::string::npos)
336  continue; // there is no rest particles in measurement
337  int _iPart = GetPartTypeNumber(part);
338  for (auto& cent : fstCentralities) {
339  // CollisionSystem:ParticleType:CentNumber:Stat/Sys:Function:FunctionVar
340  std::string stHistName{fstCollisionSystem};
341  stHistName += part;
342  stHistName += cent;
343  stHistName += "Stat";
346  TH1D* hist = (TH1D*)fMeasuredFile->Get(stHistName.c_str());
347  if (!hist) {
348  std::cerr << "AliMCSpectraWeights::Error: could not find "
349  << stHistName << "\n";
350  continue;
351  }
352 
353  // hist-> pt:multcent:PartType
354  std::array<float, 3> binEntry{0.};
355  binEntry[2] = static_cast<float>(_iPart); // particle type
356  binEntry[1] = AliMCSpectraWeights::GetMultFromCent(cent);
357  if (!fUseMultiplicity)
358  binEntry[1] = static_cast<float>(GetCentFromString(cent));
359 
360  for (int ipt = 0; ipt < fHistDataFractions->GetNbinsX(); ++ipt) {
361  binEntry[0] = fHistDataFractions->GetXaxis()->GetBinCenter(ipt);
362  if (binEntry[0] <
363  0.1) // pT cut; measurements start at 0.15 at best
364  continue;
365 
366  auto const _FractionValue =
367  hist->GetBinContent(hist->FindBin(binEntry[0]));
368  FillTH3WithValue(fHistDataFractions, binEntry, _FractionValue);
369  }
370  }
371  }
372  fMeasuredFile->Close();
373 }
374 
380  // printf("AliMCSpectraWeights::DEBUG: Loading MC histogram from input
381  // object\n");
382  if (!obj)
383  return false;
386  std::cerr
387  << "AliMCSpectraWeights::ERROR: problem with loading from object\n";
388  return false;
389  }
390  fHistMCGenPrimTrackParticle->SetDirectory(0);
391  return true;
392 }
393 
402  return false;
403  std::array<std::array<TH1D*, 10>, 20> _histMCFractions{
404  nullptr};//FIXME: WARNING HARD CODED RANGES
405  std::array<TH1D*, 20> _h1pTMCAll{nullptr};
406  for (int icent = 0; icent < fNCentralities; ++icent) {
407  auto const multTuple = AliMCSpectraWeights::GetMultTupleFromCent(icent);
408  auto const _multFront = multTuple.front();
409  auto const _multBack = multTuple.back();
410 
411  auto const _multBin1 =
412  fHistMCGenPrimTrackParticle->GetYaxis()->FindBin(_multFront);
413  auto const _multBbin2 =
414  fHistMCGenPrimTrackParticle->GetYaxis()->FindBin(_multBack);
415 
416  for (int ipart = 0; ipart < fNPartTypes;
417  ++ipart) { // calculate pt spectra of pi+k+p+Sigma+Rest
418  int const _iPart =
420  int const _iPartBin =
421  fHistMCGenPrimTrackParticle->GetZaxis()->FindBin(_iPart);
422 
423  _histMCFractions[icent][ipart] =
424  static_cast<TH1D*>(fHistMCGenPrimTrackParticle->ProjectionX(
425  Form("h1MCFraction_%s_%d", fstPartTypes[ipart].c_str(),
426  icent),
427  _multBin1, _multBbin2, _iPartBin, _iPartBin, "e"));
428  if (!_histMCFractions[icent][ipart]) {
429  std::cerr << "AliMCSpectraWeights::ERROR could not create "
430  "h1MCFraction\n";
431  return false;
432  }
433  _histMCFractions[icent][ipart]->Scale(1, "width");
434  if (!_h1pTMCAll[icent]) {
435  _h1pTMCAll[icent] = (TH1D*)_histMCFractions[icent][ipart]->Clone(Form("h1pTMCAll_%d", icent));
436  } else {
437  _h1pTMCAll[icent]->Add(_histMCFractions[icent][ipart]);
438  }
439  }
440  // all hist calculated
441  // ------------------
442  // calculate fractions
443  for (auto _hist : _histMCFractions[icent]) {
444  if (nullptr == _hist) {
445  std::cerr << "AliMCSpectraWeights::ERROR could not calculate "
446  "fraction\n";
447  continue;
448  }
449 
450  int const ipart =
451  AliMCSpectraWeights::GetPartTypeNumber(_hist->GetName());
452  TH1D* h1MCFraction =
453  (TH1D*)_hist->Clone(Form("%s_fraction", _hist->GetName()));
454  h1MCFraction->Divide(_h1pTMCAll[icent]);
455  // Set content of fractions to fHistMCFractions
456  for (int ipt = 0; ipt < fHistMCFractions->GetNbinsX(); ++ipt) {
457  float const pt =
458  fHistMCFractions->GetXaxis()->GetBinCenter(ipt);
459  if (pt < 0)
460  continue;
461  // fHistMCFractions : pt-mult-ipart
462  std::array<float, 3> binEntry{
463  pt, static_cast<float>(AliMCSpectraWeights::GetMultFromCent(icent)),
464  static_cast<float>(AliMCSpectraWeights::GetPartTypeNumber(
465  fstPartTypes[ipart]))};
466  auto const _FractionValue =
467  h1MCFraction->GetBinContent(h1MCFraction->FindBin(pt));
468  FillTH3WithValue(fHistMCFractions, binEntry, _FractionValue);
469  }
470  delete h1MCFraction;
471  }
472  }
473  for (auto& histArray : _histMCFractions) {
474  for (auto& hist : histArray) {
475  if (hist){
476  delete hist;
477  hist = nullptr;}
478  }
479  }
480  for (auto& hist : _h1pTMCAll) {
481  if(hist) {delete hist;
482  hist = nullptr;}
483  }
484  return true;
485 }
486 
498  return false;
499 
500  for (int icent = 0; icent < fNCentralities; ++icent) {
501  auto multTuple = AliMCSpectraWeights::GetMultTupleFromCent(icent);
502  auto _multFront = multTuple.front();
503  auto _multBack = multTuple.back();
504 
505  auto _bin1 =
506  fHistMCGenPrimTrackParticle->GetYaxis()->FindBin(_multFront);
507  auto _bin2 =
508  fHistMCGenPrimTrackParticle->GetYaxis()->FindBin(_multBack);
509 
510  auto h1pTMCAll = (TH1D*)fHistMCGenPrimTrackParticle->ProjectionX(
511  "h1pTMCAll", _bin1, _bin2, 1,
512  fHistMCGenPrimTrackParticle->GetNbinsZ(), "e");
513  if (!h1pTMCAll) {
514  std::cerr
515  << "AliMCSpectraWeights::ERROR could not create h1pTMCAll\n";
516  return false;
517  }
518  auto const _iRestPos = AliMCSpectraWeights::GetPartTypeNumber("Rest");
519  auto const _RestBin = fHistMCGenPrimTrackParticle->GetZaxis()->FindBin(_iRestPos);
520  auto h1RestCorrFactor = fHistMCGenPrimTrackParticle->ProjectionX(
521  Form("h1RestCorrFactor_%d", icent), _bin1, _bin2, 1, _RestBin - 1,
522  "e");
523  h1RestCorrFactor->Divide(h1pTMCAll);
524  for (int ipart = 0; ipart < fNPartTypes; ++ipart) {
525  if ("Rest" == fstPartTypes[ipart])
526  continue;
527  for (int ipt = 0; ipt < fHistDataFractions->GetNbinsX(); ++ipt) {
528  float pt = fHistDataFractions->GetXaxis()->GetBinCenter(ipt);
529  // fHistMCFractions : pt-mult-ipart
530  std::array<float, 3> binEntry{
531  pt,
532  static_cast<float>(
534  static_cast<float>(AliMCSpectraWeights::GetPartTypeNumber(
535  fstPartTypes[ipart]))};
536 
537  int const _iBinFind = fHistDataFractions->FindBin(
538  binEntry[0], binEntry[1], binEntry[2]);
539 
540  float value = fHistDataFractions->GetBinContent(_iBinFind);
541  auto _RestCorrVal = h1RestCorrFactor->GetBinContent(
542  h1RestCorrFactor->FindBin(pt));
543  if (0 == _RestCorrVal)
544  _RestCorrVal = 1;
545 
546  fHistDataFractions->SetBinContent(_iBinFind,
547  value * _RestCorrVal);
548  }
549  }
550  delete h1pTMCAll;
551  delete h1RestCorrFactor;
552  }
553  return true;
554 }
555 
566  return false;
568  return false;
569  // correction of rest particles not measured in data fractions (see
570  // AnalysisNote)
571 
572  for (int icent = 0; icent < fNCentralities; icent++) {
573  for (int ipart = 0; ipart < fNPartTypes; ipart++) {
574  // if (ipart == GetPartTypeNumber("Rest"))
575  // continue;
576  for (int ipt = 0; ipt < static_cast<int>(fBinsPt.size()); ++ipt) {
577  float pt = fHistMCWeights->GetXaxis()->GetBinCenter(ipt);
578  if (pt < 0)
579  continue;
580  std::array<float, 3> binEntry{
581  pt, static_cast<float>(GetMultFromCent(icent)),
582  static_cast<float>(AliMCSpectraWeights::GetPartTypeNumber(
583  fstPartTypes[ipart]))};
584  // Double_t binEntryData[3] = {pt, static_cast<Double_t>(icent),
585  // static_cast<Double_t>(ipart)};
586 
587  auto const _iBinMC = fHistMCFractions->FindBin(
588  binEntry[0], binEntry[1], binEntry[2]);
589  auto const _iBinData = fHistDataFractions->FindBin(
590  binEntry[0], binEntry[1], binEntry[2]);
591 
592  float const dFractionMC =
593  fHistMCFractions->GetBinContent(_iBinMC);
594 
595  float const dFractionData =
596  fHistDataFractions->GetBinContent(_iBinData);
597 
598  auto const _iBinWeight = fHistMCWeights->FindBin(
599  binEntry[0], binEntry[1], binEntry[2]);
600 
601  if (ipart == GetPartTypeNumber("Rest")) {
602  fHistMCWeights->SetBinContent(_iBinWeight, 1);
603  } else if (dFractionMC != 0)
604  fHistMCWeights->SetBinContent(_iBinWeight,
605  dFractionData / dFractionMC);
606  else
607  fHistMCWeights->SetBinContent(_iBinWeight, 1);
608 
609  fHistMCWeights->SetBinError(_iBinWeight, 1e-30);
610  }
611  }
612  }
613  return true;
614 }
615 
624  fMultOrCent = 0;
625  float lowPtCut = 0.05;
626  float eta = 0.5;
627  // if (fstCollisionSystem.find("pp") != std::string::npos)
628  // eta = 0.5;
629  AliStack* fMCStack = fMCEvent->Stack();
630  for (int ipart = 0; ipart < fMCStack->GetNtrack(); ipart++) {
631  TParticle* mcGenParticle = fMCStack->Particle(ipart);
632  if (!fMCStack->IsPhysicalPrimary(ipart))
633  continue; // secondary rejection
634  if (TMath::Abs(mcGenParticle->GetPDG()->Charge()) < 0.01)
635  continue; // neutral rejection
636  float pEta = mcGenParticle->Eta();
637  if (TMath::Abs(pEta) > eta)
638  continue; // acceptance cut
639  if (mcGenParticle->Pt() < lowPtCut)
640  continue; // TODO: hard coded low pT cut
641  ++fMultOrCent;
642  }
643 }
644 
652  TParticle* mcGenParticle, float eventMultiplicityOrCentrality) {
653  float weight = 1;
654  if (!mcGenParticle->GetPDG()) {
655  return 1;
656  }
657  if (TMath::Abs(mcGenParticle->GetPDG()->Charge()) < 0.01) {
658  return 1;
659  }
660  int particleType = AliMCSpectraWeights::IdentifyMCParticle(mcGenParticle);
661  if (particleType == GetPartTypeNumber("Rest")) {
662  return 1;
663  }
664  if (fbTaskStatus == AliMCSpectraWeights::TaskState::kMCWeightCalculated) {
665  // rest particles can not be tuned
666  float icent =
667  AliMCSpectraWeights::GetCentFromMult(eventMultiplicityOrCentrality);
668  float pt = mcGenParticle->Pt();
669  if (pt < 0.15)
670  return 1;
671  if (pt >= 20)
672  pt = 19.9;
673  std::array<float, 3> binEntry{
674  pt, static_cast<float>(AliMCSpectraWeights::GetMultFromCent(icent)),
675  static_cast<float>(particleType)};
676 
677  auto const _iBin =
678  fHistMCWeights->FindBin(binEntry[0], binEntry[1], binEntry[2]);
679 
680  weight = fHistMCWeights->GetBinContent(_iBin);
681  if (weight <= 0)
682  weight = 1;
683  }
684  return weight;
685 }
686 
693 float AliMCSpectraWeights::GetMCSpectraWeight(TParticle* mcGenParticle,
694  AliMCEvent* mcEvent) {
695  if (mcEvent != fMCEvent) {
696  fMCEvent = mcEvent;
698  }
700 }
701 
707 void AliMCSpectraWeights::FillMCSpectra(AliMCEvent* mcEvent) {
708  if (fbTaskStatus >= AliMCSpectraWeights::TaskState::kMCSpectraObtained)
709  return;
710  if (mcEvent != fMCEvent) {
711  fMCEvent = mcEvent;
713  }
714 
715  AliStack* MCStack = fMCEvent->Stack();
716  if (!MCStack) {
717  printf("AliMCSpectraWeights::ERROR: fMCStack not available\n");
718  return;
719  }
720 
721  for (int iParticle = 0; iParticle < MCStack->GetNtrack(); ++iParticle) {
722  TParticle* mcGenParticle = MCStack->Particle(iParticle);
723  if (!mcGenParticle) {
724  printf(
725  "AliMCSpectraWeights::WARNING: mcGenParticle not available\n");
726  continue;
727  }
728  if (!mcGenParticle->GetPDG())
729  continue;
730  if (!MCStack->IsPhysicalPrimary(iParticle))
731  continue;
732  if (TMath::Abs(mcGenParticle->GetPDG()->Charge()) < 0.01)
733  continue;
734 
735  float partEta = mcGenParticle->Eta();
736  float _maxEta = 0.5; // hard coded max eta; in all papers 0.5
737 
738  if (TMath::Abs(partEta) > _maxEta)
739  continue; // apply same acceptance as in published spectra
740  int particleType =
742  if (particleType < 0)
743  continue;
744  // std::array<float, 3>
745  // binEntry{static_cast<float>(mcGenParticle->Pt()),
746  // fMultOrCent,
747  // static_cast<float>(particleType)};
749  static_cast<float>(mcGenParticle->Pt()), fMultOrCent,
750  static_cast<float>(particleType));
751  }
752 }
753 
759 int AliMCSpectraWeights::IdentifyMCParticle(TParticle* mcParticle) {
760  int ipdg = TMath::Abs(
761  mcParticle
762  ->GetPdgCode()); // Abs() because antiparticles are negaitve...
763  if (ipdg == 211)
765  if (ipdg == 321)
767  if (ipdg == 2212)
768  return GetPartTypeNumber(AliMCSpectraWeights::ParticleType::kProtons);
769  if (ipdg == 3222)
770  return GetPartTypeNumber(AliMCSpectraWeights::ParticleType::kSigmaPlus);
771  if (ipdg == 3112)
772  return GetPartTypeNumber(
773  AliMCSpectraWeights::ParticleType::kSigmaMinus);
774  // if(ipdg==3334) return AliMCSpectraWeights::ParticleType::kOmegaMinus;
775  // if(ipdg==3312) return AliMCSpectraWeights::ParticleType::kXiMinus;
776  // if(ipdg==11) return AliMCSpectraWeights::ParticleType::kElectron;
777  // if(ipdg==13) return AliMCSpectraWeights::ParticleType::kMuon;
778  // printf("AliMCSpectraWeights:: pdf code of rest particle %d\n", ipdg);
779  return GetPartTypeNumber(AliMCSpectraWeights::ParticleType::kRest);
780 }
781 
787 double AliMCSpectraWeights::GetMultFromCent(int CentBin) const {
788  if (fstCollisionSystem.find("pp") != std::string::npos &&
789  fstCollisionSystem.find("ppb") == std::string::npos) {
790  // for | eta | < 0.5
791  switch (CentBin) {
792  case 0:
793  return 21.3;
794  case 1:
795  return 16.5;
796  case 2:
797  return 13.5;
798  case 3:
799  return 11.5;
800  case 4:
801  return 10.1;
802  case 5:
803  return 8.45;
804  case 6:
805  return 6.72;
806  case 7:
807  return 5.4;
808  case 8:
809  return 3.9;
810  case 9:
811  return 2.26;
812  default:
813  return -2.0;
814  }
815  } else if (fstCollisionSystem.find("ppb") !=
816  std::string::npos)
817  {
818  // for | eta | < 0.5
819  switch (CentBin) {
820  case 0:
821  return 45.0;
822  case 1:
823  return 36.2;
824  case 2:
825  return 30.5;
826  case 3:
827  return 23.2;
828  case 4:
829  return 16.1;
830  case 5:
831  return 9.8;
832  case 6:
833  return 4.4;
834  default:
835  return -2.0;
836  }
837  } else if (fstCollisionSystem.find("pbpb") != std::string::npos) {
838  // for | eta | < 0.5
839  // Centrality 0–5% 5–10% 10–20% 20–30% 30–40% 40–50% 50–60%
840  // 60–70% 70–80% dNch /dη 1601 ± 60 1294 ± 49 966±37 649±23
841  // 426±15 261±9 149±6 76±4 35±2
842  switch (CentBin) {
843  case 0:
844  return 1601;
845  case 1:
846  return 1294;
847  case 2:
848  return 966;
849  case 3:
850  return 649;
851  case 4:
852  return 426;
853  case 5:
854  return 261;
855  case 6:
856  return 149;
857  case 7:
858  return 76;
859  case 8:
860  return 35;
861  default:
862  return -2.0;
863  }
864  } else if (fstCollisionSystem.find("xexe") != std::string::npos) {
865  // for | eta | < 0.5
866  switch (CentBin) {
867  case 0:
868  return 1167;
869  case 1:
870  return 939;
871  case 2:
872  return 706;
873  case 3:
874  return 478;
875  case 4:
876  return 315;
877  case 5:
878  return 198;
879  case 6:
880  return 118;
881  case 7:
882  return 64.7;
883  case 8:
884  return 32;
885  default:
886  return -2.0;
887  }
888  }
889 
890  return -1;
891 }
892 
898 double AliMCSpectraWeights::GetMultFromCent(std::string cent) {
899  return GetMultFromCent(GetCentFromString(cent));
900 }
901 
907 std::vector<float>
909  float const currentMult = AliMCSpectraWeights::GetMultFromCent(CentBin);
910  float nextMult = 0;
911  float previousMult = 0;
912  float dMultHigh = 0;
913  float dMultLow = 0;
914  if (CentBin < fNCentralities - 1)
915  nextMult = AliMCSpectraWeights::GetMultFromCent(CentBin + 1);
916  if (CentBin > 0)
917  previousMult = AliMCSpectraWeights::GetMultFromCent(CentBin - 1);
918 
919  if (CentBin == 0) {
920  if (fstCollisionSystem.find("pp") != std::string::npos && fstCollisionSystem.find("ppb") == std::string::npos)
921  dMultHigh = 49;
922  else if (fstCollisionSystem.find("ppb") != std::string::npos)
923  dMultHigh = 299;
924  else if (fstCollisionSystem.find("pbpb") != std::string::npos)
925  dMultHigh = 3499;
926  } else
927  dMultHigh = (previousMult + currentMult) / 2.;
928  if (CentBin == fNCentralities - 1)
929  dMultLow = 0.;
930  else
931  dMultLow = (currentMult + nextMult) / 2.;
932 
933  return {dMultLow, dMultHigh};
934 }
935 
941 // TODO: implement
943  for (int i = 0; i < static_cast<int>(fstCentralities.size()); ++i) {
944  if (fstCentralities[i].find(cent) != std::string::npos)
945  return i;
946  }
947  return -1;
948 }
949 
955 // TODO: implement
957  if (fstCollisionSystem.find("pp") != std::string::npos) {
958  if (dMult > 18)
959  return 0;
960  if (dMult > 14.5)
961  return 1;
962  if (dMult > 12)
963  return 2;
964  if (dMult > 10.9)
965  return 3;
966  if (dMult > 9)
967  return 4;
968  if (dMult > 7)
969  return 5;
970  if (dMult > 5.7)
971  return 6;
972  if (dMult > 4.3)
973  return 7;
974  if (dMult > 3)
975  return 8;
976  if (dMult > 0)
977  return 9;
978  } else if (fstCollisionSystem.find("ppb") !=
979  std::string::npos) // TODO: include other systems
980  {
981  // switch (dMult) {
982  // case 0:
983  // return 45.0;
984  // case 1:
985  // return 36.2;
986  // case 2:
987  // return 30.5;
988  // case 3:
989  // return 23.2;
990  // case 4:
991  // return 16.1;
992  // case 5:
993  // return 9.8;
994  // case 6:
995  // return 4.4;
996  // default:
997  // return -2.0;
998  // }
999 
1000  } else if (fstCollisionSystem.find("pbpb") != std::string::npos) {
1001 
1002  } else if (fstCollisionSystem.find("xexe") != std::string::npos) {
1003  }
1004 
1005  return -1;
1006 }
1007 
1015  switch (type) {
1018  case AliMCSpectraWeights::ParticleType::kProtons:
1019  return AliMCSpectraWeights::ParticleType::kProtons;
1022  case AliMCSpectraWeights::ParticleType::kSigmaMinus:
1023  return AliMCSpectraWeights::ParticleType::kSigmaMinus;
1024  case AliMCSpectraWeights::ParticleType::kSigmaPlus:
1025  return AliMCSpectraWeights::ParticleType::kSigmaPlus;
1026  case AliMCSpectraWeights::ParticleType::kRest:
1027  return AliMCSpectraWeights::ParticleType::kRest;
1028  default:
1029  return -1;
1030  break;
1031  }
1032 }
1033 
1039 int AliMCSpectraWeights::GetPartTypeNumber(std::string Particle) {
1040  if (Particle.find("Pion") != std::string::npos) {
1043  } else if (Particle.find("Proton") != std::string::npos) {
1045  AliMCSpectraWeights::ParticleType::kProtons);
1046  } else if (Particle.find("Kaon") != std::string::npos) {
1049  } else if (Particle.find("SigmaMinus") != std::string::npos) {
1051  AliMCSpectraWeights::ParticleType::kSigmaMinus);
1052  } else if (Particle.find("SigmaPlus") != std::string::npos) {
1054  AliMCSpectraWeights::ParticleType::kSigmaPlus);
1055  } else if (Particle.find("Rest") != std::string::npos) {
1057  AliMCSpectraWeights::ParticleType::kRest);
1058  } else
1059  return -1;
1060 }
1061 
1062 // TODO: implement flags
1069  switch (flag) {
1070  case SysFlag::kNominal:
1071  return "Bylinkin";
1072  case SysFlag::kPionUp:
1073  return "";
1074  case SysFlag::kPionDown:
1075  return "";
1076  case SysFlag::kProtonUp:
1077  return "";
1078  case SysFlag::kProtonDown:
1079  return "";
1080  case SysFlag::kKaonUp:
1081  return "";
1082  case SysFlag::kKaonDown:
1083  return "";
1084  case SysFlag::kSigmaPlusUp:
1085  return "";
1086  case SysFlag::kSigmaPlusDown:
1087  return "";
1088  case SysFlag::kSigmaMinusUp:
1089  return "";
1090  case SysFlag::kSigmaMinusDown:
1091  return "";
1092  case SysFlag::kBylinkinLower:
1093  return "";
1094  case SysFlag::kBylinkinUpper:
1095  return "";
1096  case SysFlag::kHagedorn:
1097  return "";
1098  case SysFlag::kHagedornUpper:
1099  return "";
1100  case SysFlag::kHagedornLower:
1101  return "";
1102  case SysFlag::kExponential:
1103  return "";
1104  case SysFlag::kExponentialUpper:
1105  return "";
1106  case SysFlag::kExponentialLower:
1107  return "";
1108  case SysFlag::kBlastwave:
1109  return "";
1110  case SysFlag::kBlastwaveUpper:
1111  return "";
1112  case SysFlag::kBlastwaveLower:
1113  return "";
1114  default:
1115  return "Bylinkin";
1116  }
1117 
1118  return "Bylinkin";
1119 }
1120 // TODO: implement
1122  switch (flag) {
1123  case SysFlag::kNominal:
1124  return "";
1125  case SysFlag::kPionUp:
1126  return "";
1127  case SysFlag::kPionDown:
1128  return "";
1129  case SysFlag::kProtonUp:
1130  return "";
1131  case SysFlag::kProtonDown:
1132  return "";
1133  case SysFlag::kKaonUp:
1134  return "";
1135  case SysFlag::kKaonDown:
1136  return "";
1137  case SysFlag::kSigmaPlusUp:
1138  return "";
1139  case SysFlag::kSigmaPlusDown:
1140  return "";
1141  case SysFlag::kSigmaMinusUp:
1142  return "";
1143  case SysFlag::kSigmaMinusDown:
1144  return "";
1145  case SysFlag::kBylinkinLower:
1146  return "";
1147  case SysFlag::kBylinkinUpper:
1148  return "";
1149  case SysFlag::kHagedorn:
1150  return "";
1151  case SysFlag::kHagedornUpper:
1152  return "";
1153  case SysFlag::kHagedornLower:
1154  return "";
1155  case SysFlag::kExponential:
1156  return "";
1157  case SysFlag::kExponentialUpper:
1158  return "";
1159  case SysFlag::kExponentialLower:
1160  return "";
1161  case SysFlag::kBlastwave:
1162  return "";
1163  case SysFlag::kBlastwaveUpper:
1164  return "";
1165  case SysFlag::kBlastwaveLower:
1166  return "";
1167  default:
1168  return "";
1169  }
1170 
1171  return "";
1172 }
void FillMCSpectra(AliMCEvent *mcEvent)
void Init()
Initialisation of object.
Definition: External.C:260
std::string fstCollisionSystem
void SetBinsPt(std::vector< double > bins)
function to set pt binning of all internal histograms
void FillTH3WithValue(TH3F *h, std::array< float, 3 > _xyzValue, float _value)
bool CalculateMCWeights()
calculate weight factors using internal information
TCanvas * c
Definition: TestFitELoss.C:172
void SetBinsMultCent(std::vector< double > bins)
function to set multiplicity binning of all internal histograms
int GetCentFromString(std::string cent)
bool LoadFromAliMCSpectraWeight(AliMCSpectraWeights *obj)
Load the information from a previous train output.
AliMCSpectraWeights()
default constructor
std::string GetFunctionFromSysFlag(SysFlag flag)
void LoadMeasuredFractions()
Load measured fractions (expert input) from alien.
bool CalcMCFractions()
calculate the relative fractions of all particle species in MC
std::vector< std::string > fstPartTypes
double GetCentFromMult(double dMult)
Definition: External.C:212
bool CorrectFractionsforRest()
correct the expert input for remaining particle species
int GetBinFromTH3(TH3F *h, std::array< float, 3 > const &_values)
std::vector< float > GetMultTupleFromCent(int CentBin) const
TH3F * GetHistMCGenPrimTrackParticles() const
double GetMultFromCent(int CentBin) const
TList * listMC
std::vector< std::string > fstCentralities
int IdentifyMCParticle(TParticle *mcParticle)
int GetPartTypeNumber(ParticleType type)
std::vector< float > fBinsPt
std::vector< float > fBinsMultCent
Class for re-weighting particle abundances in MC simulations.
std::string GetSysVarFromSysFlag(SysFlag flag)
float GetMCSpectraWeight(TParticle *mcGenParticle, float eventMultiplicityOrCentrality)
void InitHistos()
Create all internal histograms.
void CountEventMult()
count the number of charged particles in the current event