Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
Femto1DCFAnaMapMCRoco.cxx
1/*
2 * Femto1DCFAnaMapMCRoco.cxx
3 *
4 * Created on: 2 kwi 2018
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "Femto1DCFAnaMapMCRoco.h"
11
12#include "Const.h"
13#include "CorrFitMapKstarRstar.h"
14#include "Cout.h"
15#include "DividedHisto.h"
16#include "FastHist.h"
17#include "Femto1DCF.h"
18#include "FemtoPair.h"
19#include "FemtoSourceDensity.h"
20#include "FemtoSourceModel.h"
21#include "FemtoSourceModelNumerical1D.h"
22#include "StdHist.h"
23
24#include <TAxis.h>
25#include <TDatabasePDG.h>
26#include <TFile.h>
27#include <TH1.h>
28#include <TH2.h>
29#include <TLorentzVector.h>
30#include <TMath.h>
31#include <TParticlePDG.h>
32#include <TRandom.h>
33#include <TVector3.h>
34#include <iostream>
35
36
37namespace Hal {
38 Femto1DCFAnaMapMCRoco::Femto1DCFAnaMapMCRoco() : Femto1DMapGenerator() {}
39
40 Femto1DCFAnaMapMCRoco::~Femto1DCFAnaMapMCRoco() {
41 if (fGeneratorIntegrated) delete fGeneratorIntegrated;
42 if (fSourceParams) delete fSourceParams;
43 if (fSampleRandom) delete fSampleRandom;
44 }
45
46 void Femto1DCFAnaMapMCRoco::Exec(Int_t pairs_per_bin, Bool_t autoscale) {
47 if (autoscale) pairs_per_bin = (Double_t) pairs_per_bin * fIntegralScale;
48 const Int_t pointsQ = fMap->GetNum()->GetNbinsX() + 1;
49 Double_t* kstar = new Double_t[pointsQ];
50 Double_t* kfill = new Double_t[pointsQ];
51 Double_t** parametrizations = new Double_t*[fRBins];
52 FemtoSourceModel* sourceModel = fGenerator->GetSourceModel();
53 FemtoSourceModel* sourceModelntegrated = fGeneratorIntegrated->GetSourceModel();
54 FemtoSourceDensity* densityModel = sourceModel->GetDensityProb();
55 FemtoSourceDensity* densityIntegratedModel = sourceModelntegrated->GetDensityProb();
56 Int_t sourceParamsNo = sourceModel->GetNpar();
57 for (int i = 0; i < fRBins; i++) {
58 parametrizations[i] = new Double_t[sourceParamsNo];
59 }
60
61 for (int i = 1; i <= fMap->GetNum()->GetNbinsX(); i++) {
62 kstar[i] = fMap->GetNum()->GetXaxis()->GetBinCenter(i);
63 kfill[i] = kstar[i];
64 if (fKinematics == Femto::EKinematics::kLCMS) { kstar[i] = kstar[i] * 0.5; }
65 }
66 TVector3 boost(0.1, 0.1, 0.1);
67
68 const Double_t scale = 1.0 / TMath::Sqrt(3.0);
69
70
71 for (int i = 0; i < fRBins; i++) {
72 for (int j = 0; j < sourceParamsNo; j++) {
73 parametrizations[i][j] = fSourceParams[j];
74 }
75 }
76 switch (fModelType) {
77 case EModelType::k1dModel: {
78 for (int i = 0; i < fRBins; i++) {
79 parametrizations[i][0] = fRadiiBins[i];
80 }
81 } break;
82 case EModelType::k3dModel: {
83 for (int i = 0; i < fRBins; i++) {
84 parametrizations[i][0] = fRadiiBins[i] * scale;
85 parametrizations[i][1] = fRadiiBins[i] * scale;
86 parametrizations[i][2] = fRadiiBins[i] * scale;
87 }
88 } break;
89 case EModelType::kOther: {
90 // DO NOTHING TODO check
91 } break;
92 };
93 // Int_t mainBin = fRBins * 0.5;
94 // if (refRadius != 0) { mainBin = fMap->GetNum()->GetYaxis()->FindBin(refRadius); }
95 // mainBin = 0;
96 // fGenerator->GetSourceModel()->SetRadius(fRadiiBins[mainBin]); // generate pairs for only one parametrization
97
98 TH1D* monGaus1 = new TH1D("monG1", "monG", 500, 0, 50);
99 TH1D* monGaus2 = new TH1D("monG3", "monG", 500, 0, 50);
100 TH1D* monGaus3 = new TH1D("monG6", "monG", 500, 0, 50);
101 TH1D* monGaus4 = new TH1D("monG9", "monG", 500, 0, 50);
102 TH1D* monGaus5 = new TH1D("monG10", "monG", 500, 0, 50);
103 TH1D* monRaw = new TH1D("monR", "monR", 500, 0, 50);
104
105 TDatabasePDG* pdg = TDatabasePDG::Instance();
106 TParticlePDG* pid1 = pdg->GetParticle(fPid1);
107 TParticlePDG* pid2 = pdg->GetParticle(fPid2);
108 Double_t m1 = Const::PionPlusMass();
109 Double_t m2 = Const::PionPlusMass();
110 if (pid1) m1 = pid1->Mass();
111 if (pid2) m2 = pid2->Mass();
112 m1 = m1 * m1;
113 m2 = m2 * m2;
114 Int_t nbinsX, nbinsY;
115 Double_t minX, maxX, minY, maxY;
116 Hal::Std::GetAxisPar(*fMap->GetNum(), nbinsX, minX, maxX, "x");
117 Hal::Std::GetAxisPar(*fMap->GetNum(), nbinsY, minY, maxY, "y");
118
119 FastHist2D* num1 = new FastHist2D("2dnum", "2dnum", nbinsX, minX, maxX, nbinsY, minY, maxY);
120 FastHist2D* num2 = new FastHist2D("2dden", "2dden", nbinsX, minX, maxX, nbinsY, minY, maxY);
121 for (int ikst = 1; ikst <= fMap->GetNum()->GetNbinsX(); ikst++) {
122 Double_t E1 = TMath::Sqrt(m1 + kstar[ikst] * kstar[ikst]);
123 Double_t E2 = TMath::Sqrt(m2 + kstar[ikst] * kstar[ikst]);
124 Double_t px, py, pz;
125 gRandom->Sphere(px, py, pz, kstar[ikst]);
126 TLorentzVector p1(px, py, pz, E1);
127 TLorentzVector p2(-px, -py, -pz, E2);
128 p1.Boost(boost);
129 p2.Boost(boost);
130 fPair->SetTrueMomenta1(p1.X(), p1.Y(), p1.Z(), p1.T());
131 fPair->SetTrueMomenta2(p2.X(), p2.Y(), p2.Z(), p2.T());
132
133 for (int i = 0; i < pairs_per_bin; i++) {
134 fGeneratorIntegrated->GenerateFreezoutCooordinates(fPair);
135 Double_t weight = fWeight->GenerateWeight(fPair);
136 TVector3 Radius(sourceModelntegrated->GetROut(), sourceModelntegrated->GetRSide(), sourceModelntegrated->GetRLong());
137 // Double_t RadiusFor1D[3] = {0, 0, 0};
138 Double_t Rinv = Radius.Mag();
139 // RadiusFor1D[0] = Rinv * TMath::Sqrt(3);
140 Double_t refWeight = 1; // 1.0 / sourceModelntegrated->GetProbDensity3d(Radius, nullptr);
141 refWeight = 1.0 / densityIntegratedModel->GetProbDensity1d(Rinv, nullptr);
142 if (TMath::IsNaN(refWeight)) {
143 i--;
144 continue;
145 }
146 // std::cout << "\t" << std::endl;
147 for (int r_bin = 0; r_bin < fRBins; r_bin++) {
148 Double_t R = fRadiiBins[r_bin];
149 Double_t newWeight = 1; // sourceModel->GetProbDensity3d(Radius, parametrizations[r_bin]); // why?
150 newWeight = densityModel->GetProbDensity1d(Rinv, parametrizations[r_bin]); // why?
151 Double_t effWeight = newWeight * refWeight;
152 if (fDebugDistribution) {
153 if (r_bin == 1) {
154 monRaw->Fill(Rinv, 1);
155 monGaus1->Fill(Rinv, effWeight);
156 }
157 if (r_bin == 10) {
158 // monRaw->Fill(Rinv, 1);
159 monGaus2->Fill(Rinv, effWeight);
160 }
161 if (r_bin == 25) {
162 // monRaw->Fill(Rinv, 1);
163 monGaus3->Fill(Rinv, effWeight);
164 }
165 if (r_bin == 50) {
166 // monRaw->Fill(Rinv, 1);
167 monGaus4->Fill(Rinv, effWeight);
168 }
169 if (r_bin == 75) {
170 // monRaw->Fill(Rinv, 1);
171 monGaus5->Fill(Rinv, effWeight);
172 }
173 }
174 // std::cout << r_bin << " " << newWeight << " " << refWeight << std::endl;
175 Int_t bin = num1->FindBin(kfill[ikst], R);
176 num1->IncrementRawBinContent(bin, weight * effWeight);
177 num2->IncrementRawBinContent(bin, effWeight);
178 }
179 }
180 }
181 for (int i = 0; i <= nbinsX + 1; i++) {
182 for (int j = 0; j <= nbinsY + 1; j++) {
183 fMap->GetNum()->SetBinContent(i, j, num1->GetBinContent(i, j));
184 fMap->GetDen()->SetBinContent(i, j, num2->GetBinContent(i, j));
185 }
186 }
187 delete num1;
188 delete num2;
189 if (fDebugDistribution) {
190 TFile* fx = new TFile("ctrl.root", "recreate");
191 monRaw->Write();
192 monGaus1->Write();
193 monGaus2->Write();
194 monGaus3->Write();
195 monGaus4->Write();
196 monGaus5->Write();
197 fx->Close();
198 } else {
199 delete monRaw;
200 delete monGaus1;
201 delete monGaus2;
202 delete monGaus3;
203 delete monGaus4;
204 delete monGaus5;
205 }
206 /*
207 mainBin = fRbins * 0.5;
208 fGenerator->GetSourceModel()->SetRadius(fRadiiBins[mainBin]);
209 for (int ikst = 1; ikst <= fMap->GetNum()->GetNbinsX(); ikst++) {
210 Double_t E1 = TMath::Sqrt(fM1 + kstar[ikst] * kstar[ikst]);
211 Double_t E2 = TMath::Sqrt(fM2 + kstar[ikst] * kstar[ikst]);
212 Double_t px, py, pz;
213 gRandom->Sphere(px, py, pz, kstar[ikst]);
214 TLorentzVector p1(px, py, pz, E1);
215 TLorentzVector p2(-px, -py, -pz, E2);
216 p1.Boost(boost);
217 p2.Boost(boost);
218 fPair->SetTrueMomenta1(p1.X(), p1.Y(), p1.Z(), p1.T());
219 fPair->SetTrueMomenta2(p2.X(), p2.Y(), p2.Z(), p2.T());
220
221 for (int i = 0; i < pairs_per_bin; i++) {
222 fGenerator->GenerateFreezoutCooordinates(fPair);
223 Double_t weight = fWeight->GenerateWeight(fPair);
224 Double_t Radius[3] = {sourceModel->GetROut(), sourceModel->GetRSide(), sourceModel->GetRLong()};
225 Double_t refWeight = 1.0 / sourceModel->GetProbDensity(Radius, parametrizations[mainBin]);
226 for (int r_bin = 0; r_bin < fRbins; r_bin++) {
227 Double_t R = fRadiiBins[r_bin];
228 Double_t newWeight = sourceModel->GetProbDensity(Radius, parametrizations[r_bin]);
229 Double_t effWeight = newWeight * refWeight;
230 ((TH2*) fMap->GetNum())->Fill(kfill[ikst], R, weight * effWeight);
231 ((TH2*) fMap->GetDen())->Fill(kfill[ikst], R, effWeight);
232 }
233 }
234 }
235
236
237 mainBin = fRbins * 0.9;
238 fGenerator->GetSourceModel()->SetRadius(fRadiiBins[mainBin]);
239 for (int ikst = 1; ikst <= fMap->GetNum()->GetNbinsX(); ikst++) {
240 Double_t E1 = TMath::Sqrt(fM1 + kstar[ikst] * kstar[ikst]);
241 Double_t E2 = TMath::Sqrt(fM2 + kstar[ikst] * kstar[ikst]);
242 Double_t px, py, pz;
243 gRandom->Sphere(px, py, pz, kstar[ikst]);
244 TLorentzVector p1(px, py, pz, E1);
245 TLorentzVector p2(-px, -py, -pz, E2);
246 p1.Boost(boost);
247 p2.Boost(boost);
248 fPair->SetTrueMomenta1(p1.X(), p1.Y(), p1.Z(), p1.T());
249 fPair->SetTrueMomenta2(p2.X(), p2.Y(), p2.Z(), p2.T());
250
251 for (int i = 0; i < pairs_per_bin; i++) {
252 fGenerator->GenerateFreezoutCooordinates(fPair);
253 Double_t weight = fWeight->GenerateWeight(fPair);
254 Double_t Radius[3] = {sourceModel->GetROut(), sourceModel->GetRSide(), sourceModel->GetRLong()};
255 Double_t refWeight = 1.0 / sourceModel->GetProbDensity(Radius, parametrizations[mainBin]);
256 for (int r_bin = 0; r_bin < fRbins; r_bin++) {
257 Double_t R = fRadiiBins[r_bin];
258 Double_t newWeight = sourceModel->GetProbDensity(Radius, parametrizations[r_bin]);
259 Double_t effWeight = newWeight * refWeight;
260 ((TH2*) fMap->GetNum())->Fill(kfill[ikst], R, weight * effWeight);
261 ((TH2*) fMap->GetDen())->Fill(kfill[ikst], R, effWeight);
262 }
263 }
264 }
265 */
266
267 for (int i = 0; i < fRBins; i++) {
268 delete[] parametrizations[i];
269 }
270 delete[] parametrizations;
271 delete[] kfill;
272 delete[] kstar;
273 }
274
275 void Femto1DCFAnaMapMCRoco::SaveMap(TString filename) {
276 TFile* file = new TFile(filename, "recreate");
277 CorrFitMapKstarRstarDiv* map = new CorrFitMapKstarRstarDiv(*fMap, fKinematics);
278 map->Write("map");
279 file->Close();
280 }
281
282 Bool_t Femto1DCFAnaMapMCRoco::Init() {
283 if (fGenerator == nullptr) return kFALSE;
284 if (fGenerator->GetSourceModel()->GetModelNumProp() != FemtoSourceModel::ENumProperty::kFullyAnalytical) {
285 Cout::PrintInfo(" Femto1DCFAnaMapMCRoco::Init - cannot use nonanalytical source emission function", EInfo::kLowWarning);
286 return kFALSE;
287 }
288 FemtoSourceModel* sourceModel = fGenerator->GetSourceModel();
289 Int_t sourceParamsNo = sourceModel->GetNpar();
290 fSourceParams = new Double_t[sourceParamsNo];
291 for (int i = 0; i < sourceParamsNo; i++) {
292 fSourceParams[i] = sourceModel->GetParameter(i);
293 }
294
295 if (dynamic_cast<FemtoSourceModel1D*>(sourceModel)) {
296 fModelType = EModelType::k1dModel;
297 } else if (dynamic_cast<FemtoSourceModel3D*>(sourceModel)) {
298 fModelType = EModelType::k3dModel;
299 } else {
300 Cout::PrintInfo(" Femto1DCFAnaMapMCRoco::Init - cannot use analytical model that not base from 1d or 3d source model",
301 EInfo::kLowWarning);
302 return kFALSE;
303 }
304
305 fPair = Femto::MakePair(fKinematics, kFALSE);
306 fMap = new DividedHisto2D("map", fKStarBins, fKStarMin, fKStarMax, fRBins, fRMin, fRMax, 'D');
307 fMap->SetDirectory(nullptr);
308 fRStep = (fRMax - fRMin) / ((Double_t) fRBins);
309 fRMinEff = fRMin + 0.5 * fRStep;
310 fRadiiBins.MakeBigger(fRBins);
311 for (int r_bin = 0; r_bin < fRBins; r_bin++) {
312 fRadiiBins[r_bin] = fRMinEff + ((double) r_bin) * fRStep;
313 }
314
315 TDatabasePDG* pdg = TDatabasePDG::Instance();
316 Int_t pid1 = fWeight->GetPdg1();
317 Int_t pid2 = fWeight->GetPdg2();
318 TParticlePDG* part1 = pdg->GetParticle(pid1);
319 TParticlePDG* part2 = pdg->GetParticle(pid2);
320 if (part1 == nullptr) return kFALSE;
321 if (part2 == nullptr) return kFALSE;
322 fPair->SetPdg1(pid1);
323 fPair->SetPdg2(pid2);
324 /* fM1 = part1->Mass();
325 fM2 = part2->Mass();
326 fM1 = fM1 * fM1;
327 fM2 = fM2 * fM2;*/
328 Double_t rmin = 0;
329 Double_t rmax = fRadiiBins[fRBins - 1] * 5;
330 if (fSampleRandom) {
331 delete fSampleRandom;
332 fSampleRandom = nullptr;
333 }
334 fSampleRandom = new TH1D("samram", "samram", 50, rmin, rmax);
335 Double_t params[1] = {0};
336 Double_t minIntegral = 1E+9;
337 FemtoSourceDensity* sourceBase = fGenerator->GetSourceModel()->GetDensityProb();
338 for (int r_bin = 0; r_bin < fRBins; r_bin++) {
339 params[0] = fRadiiBins[r_bin];
340 Double_t localIntegral = 0;
341 for (int i = 1; i <= fSampleRandom->GetNbinsX(); i++) {
342 Double_t r = fSampleRandom->GetXaxis()->GetBinCenter(i);
343 Double_t val = sourceBase->GetProbDensity1d(r, params);
344 localIntegral += val;
345 if (val > fSampleRandom->GetBinContent(i)) { fSampleRandom->SetBinContent(i, val); }
346 }
347 if (localIntegral < minIntegral) { minIntegral = localIntegral; }
348 }
349 Double_t maxIntegral = 0;
350 for (int i = 1; i <= fSampleRandom->GetNbinsX(); i++) {
351 maxIntegral += fSampleRandom->GetBinContent(i);
352 }
353 fIntegralScale = maxIntegral / minIntegral;
354 Cout::Text(Form("Integral scale = %4.2f ", fIntegralScale), "L", kYellow);
355 fGeneratorIntegrated = fGenerator->MakeCopy();
356 FemtoSourceModelNumerical1D model;
357 model.SetRadiusDistribution(*fSampleRandom);
358 fGeneratorIntegrated->SetSourceModel(model);
359 delete fSampleRandom;
360 fSampleRandom = nullptr;
361 return kTRUE;
362 }
363} // namespace Hal
Hal::FemtoSourceDensity::GetProbDensity1d
virtual Double_t GetProbDensity1d(const Double_t, const Double_t *) const
Definition FemtoSourceDensity.h:60
Hal::FemtoSourceModel::GetDensityProb
FemtoSourceDensity * GetDensityProb() const
Definition FemtoSourceModel.h:144
Hal::FemtoSourceModel::GetParameter
Double_t GetParameter(Int_t n) const
Definition FemtoSourceModel.h:132
Hal::FemtoSourceModel::GetRLong
Double_t GetRLong() const
Definition FemtoSourceModel.h:88
Hal::FemtoSourceModel::GetRSide
Double_t GetRSide() const
Definition FemtoSourceModel.h:83
Hal::FemtoSourceModel::GetROut
Double_t GetROut() const
Definition FemtoSourceModel.h:78