Skip to content

Commit abcc778

Browse files
joachimckhjoachimckh
authored
Common: ESE splines for FT0-A and FV0-A (#8149)
1 parent f01452c commit abcc778

4 files changed

Lines changed: 240 additions & 555 deletions

File tree

‎Common/Core/FFitWeights.cxx‎

Lines changed: 85 additions & 209 deletions
Original file line numberDiff line numberDiff line change
@@ -10,42 +10,37 @@
1010
// or submit itself to any jurisdiction.
1111

1212
/// \file FFitWeights.cxx
13-
/// \brief Implementation file for FFitWeights.h, see the header fore more information
13+
/// \brief Implementation file for FFitWeights.h, see the header for more information
1414
///
15-
/// \author Joachim C. K. B. Hansen, Lund University
15+
/// \author Joachim C. K. B. Hansen
1616

1717
#include "FFitWeights.h"
1818

19-
// using std::complex;
20-
// using std::pair;
21-
// using std::string;
22-
// using std::vector;
19+
#include <TSpline.h>
2320

2421
ClassImp(FFitWeights)
2522

2623
FFitWeights::FFitWeights() : TNamed("", ""),
2724
fW_data{nullptr},
28-
vGain{0},
2925
CentBin{100},
30-
ChIDBin{220},
31-
sAmpl{nullptr},
32-
sqVec{nullptr},
33-
sqCorVec{nullptr}
26+
qAxis{nullptr},
27+
nResolution{3000},
28+
qnTYPE{0}
3429
{
3530
}
3631

3732
FFitWeights::FFitWeights(const char* name) : TNamed(name, name),
3833
fW_data{nullptr},
39-
vGain{0},
4034
CentBin{100},
41-
ChIDBin{220},
42-
sAmpl{nullptr},
43-
sqVec{nullptr},
44-
sqCorVec{nullptr} {}
35+
qAxis{nullptr},
36+
nResolution{3000},
37+
qnTYPE{0} {}
4538

4639
FFitWeights::~FFitWeights()
4740
{
4841
delete fW_data;
42+
if (qAxis)
43+
delete qAxis;
4944
};
5045

5146
void FFitWeights::Init()
@@ -54,241 +49,122 @@ void FFitWeights::Init()
5449
fW_data->SetName("FFitWeights_Data");
5550
fW_data->SetOwner(kTRUE);
5651

57-
this->SetBinAxis(1000, 0, 5000, 0);
58-
this->SetBinAxis(250, -3500, 3500, 1);
59-
this->SetBinAxis(250, -250, 250, 2);
60-
61-
const char* tnd = "FT0Ampl";
62-
fW_data->Add(new TH2F(tnd, ";channel;amplitude", ChIDBin, 0, ChIDBin, sAmpl->GetNbins(), sAmpl->GetXmin(), sAmpl->GetXmax()));
63-
fW_data->Add(new TH2F(Form("%sCorr", tnd), ";channel;amplitude", ChIDBin, 0, ChIDBin, sAmpl->GetNbins(), sAmpl->GetXmin(), sAmpl->GetXmax()));
64-
65-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C", "x", 2), ";Centrality;Qx_{2}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
66-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C", "y", 2), ";Centrality;Qy_{2}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
67-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C", "x", 3), ";Centrality;Qx_{3}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
68-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C", "y", 3), ";Centrality;Qy_{3}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
69-
70-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_Rec", "x", 2), ";Centrality;Qx_{2}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
71-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_Rec", "y", 2), ";Centrality;Qy_{2}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
72-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_Rec", "x", 3), ";Centrality;Qx_{3}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
73-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_Rec", "y", 3), ";Centrality;Qy_{3}", CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
74-
75-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_RecTot", "x", 2), ";Centrality;Qx_{2}", CentBin, 0, CentBin, sqCorVec->GetNbins(), sqCorVec->GetXmin(), sqCorVec->GetXmax()));
76-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_RecTot", "y", 2), ";Centrality;Qy_{2}", CentBin, 0, CentBin, sqCorVec->GetNbins(), sqCorVec->GetXmin(), sqCorVec->GetXmax()));
77-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_RecTot", "x", 3), ";Centrality;Qx_{3}", CentBin, 0, CentBin, sqCorVec->GetNbins(), sqCorVec->GetXmin(), sqCorVec->GetXmax()));
78-
fW_data->Add(new TH2F(Form("hQ%s%i_FT0C_RecTot", "y", 3), ";Centrality;Qy_{3}", CentBin, 0, CentBin, sqCorVec->GetNbins(), sqCorVec->GetXmin(), sqCorVec->GetXmax()));
52+
if (!qAxis)
53+
this->SetBinAxis(500, 0, 25);
54+
for (const auto& qn : qnTYPE) {
55+
fW_data->Add(new TH2D(this->GetQName(qn.first, qn.second.c_str()), this->GetAxisName(qn.first, qn.second.c_str()), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax()));
56+
}
7957
};
8058

81-
void FFitWeights::FillFT0(std::size_t iCh, float amplitude, float GainCst)
59+
void FFitWeights::Fill(float centrality, float qn, int nh, const char* pf)
8260
{
8361
TObjArray* tar{nullptr};
8462

8563
tar = fW_data;
8664
if (!tar)
8765
return;
8866

89-
TH2F* th2 = reinterpret_cast<TH2F*>(tar->FindObject("FT0Ampl"));
67+
TH2D* th2 = reinterpret_cast<TH2D*>(tar->FindObject(this->GetQName(nh, pf)));
9068
if (!th2) {
91-
tar->Add(new TH2F("FT0Ampl", ";channel;amplitude", ChIDBin, 0, ChIDBin, sAmpl->GetNbins(), sAmpl->GetXmin(), sAmpl->GetXmax()));
92-
th2 = reinterpret_cast<TH2F*>(tar->At(tar->GetEntries() - 1));
69+
tar->Add(new TH2D(this->GetQName(nh, pf), this->GetAxisName(nh, pf), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax()));
70+
th2 = reinterpret_cast<TH2D*>(tar->At(tar->GetEntries() - 1));
9371
}
94-
th2->Fill(iCh, amplitude);
95-
96-
TH2F* th2Cor = reinterpret_cast<TH2F*>(tar->FindObject("FT0AmplCorr"));
97-
if (!th2Cor) {
98-
tar->Add(new TH2F("FT0AmplCorr", ";channel;amplitude", ChIDBin, 0, ChIDBin, sAmpl->GetNbins(), sAmpl->GetXmin(), sAmpl->GetXmax()));
99-
th2Cor = reinterpret_cast<TH2F*>(tar->At(tar->GetEntries() - 1));
100-
}
101-
th2Cor->Fill(iCh, amplitude / GainCst);
72+
th2->Fill(centrality, qn);
10273
};
10374

104-
void FFitWeights::FillQ(float mult, float vec, int nHarm, const char* coord, const char* qType)
75+
Long64_t FFitWeights::Merge(TCollection* collist)
10576
{
106-
TObjArray* tar{nullptr};
107-
108-
tar = fW_data;
109-
if (!tar)
110-
return;
111-
112-
TH2F* th2 = reinterpret_cast<TH2F*>(tar->FindObject(Form("hQ%s%i_FT0C%s", coord, nHarm, qType)));
113-
if (!th2) {
114-
tar->Add(new TH2F(Form("hQ%s%i_FT0C%s", coord, nHarm, qType), Form(";Centrality;Q%s_{%i}", coord, nHarm), CentBin, 0, CentBin, sqVec->GetNbins(), sqVec->GetXmin(), sqVec->GetXmax()));
115-
th2 = reinterpret_cast<TH2F*>(tar->At(tar->GetEntries() - 1));
77+
Long64_t nmerged = 0;
78+
if (!fW_data) {
79+
fW_data = new TObjArray();
80+
fW_data->SetName("FFitWeights_Data");
81+
fW_data->SetOwner(kTRUE);
82+
}
83+
FFitWeights* l_w = 0;
84+
TIter all_w(collist);
85+
while ((l_w = (reinterpret_cast<FFitWeights*>(all_w())))) {
86+
AddArray(fW_data, l_w->GetDataArray());
87+
nmerged++;
11688
}
117-
th2->Fill(mult, vec);
89+
return nmerged;
11890
};
119-
120-
void FFitWeights::CreateGain()
91+
void FFitWeights::AddArray(TObjArray* targ, TObjArray* sour)
12192
{
122-
vGain.clear();
123-
124-
TH1D* h1;
125-
if (fW_data->GetEntries() < 1)
93+
if (!sour) {
94+
printf("Source array does not exist!\n");
12695
return;
127-
128-
TH2F* hGain = reinterpret_cast<TH2F*>(fW_data->At(0)->Clone("FT0Ampl"));
129-
double vMean{0}; // = hGain->GetMean(2);
130-
// c-side first bin ich: 97 (last 208)'
131-
// gain split for c-side 144
132-
hGain->GetXaxis()->SetRangeUser(96, 144);
133-
float vMeanC_inner = hGain->GetMean(2);
134-
135-
hGain->GetXaxis()->SetRangeUser(144, 208);
136-
float vMeanC_outer = hGain->GetMean(2);
137-
138-
hGain->GetXaxis()->SetRangeUser(0, -1);
139-
140-
for (int iCh{0}; iCh < hGain->GetNbinsX(); iCh++) {
141-
h1 = static_cast<TH1D*>(hGain->ProjectionY(Form("proj%i", iCh), iCh, iCh));
142-
double mean = h1->GetMean();
143-
double meanErr = h1->GetMeanError();
144-
145-
if (iCh > 95 && iCh < 144)
146-
vMean = vMeanC_inner;
147-
else if (iCh > 144)
148-
vMean = vMeanC_outer;
149-
double fWeight = mean / vMean;
150-
151-
if (fWeight > 0) {
152-
vGain.push_back(fWeight);
96+
}
97+
for (int i = 0; i < sour->GetEntries(); i++) {
98+
TH2D* sourh = reinterpret_cast<TH2D*>(sour->At(i));
99+
TH2D* targh = reinterpret_cast<TH2D*>(targ->FindObject(sourh->GetName()));
100+
if (!targh) {
101+
targh = reinterpret_cast<TH2D*>(sourh->Clone(sourh->GetName()));
102+
targh->SetDirectory(0);
103+
targ->Add(targh);
153104
} else {
154-
vGain.push_back(1.0);
105+
targh->Add(sourh);
155106
}
156-
157-
TObjArray* tar = fW_data;
158-
if (!tar)
159-
return;
160-
161-
fW_data->Add(new TH1F("FT0MultCorr", ";channel", ChIDBin, 0, ChIDBin));
162-
TH1F* htmp = reinterpret_cast<TH1F*>(tar->FindObject("FT0MultCorr"));
163-
if (!htmp)
164-
return;
165-
166-
htmp->SetBinContent(iCh, mean);
167-
htmp->SetBinError(iCh, meanErr);
168107
}
169-
170-
// note to self if FT0A is implemented this has to be done differently
171-
};
172-
173-
std::vector<float> FFitWeights::GetGain()
174-
{
175-
return vGain;
176108
};
177109

178-
void FFitWeights::CreateRecenter(const char* xy)
110+
void FFitWeights::qSelectionSpline(std::vector<int> nhv, std::vector<std::string> stv) /* only execute OFFLINE */
179111
{
180-
if (fW_data->GetEntries() < 1)
181-
return;
182-
183112
TObjArray* tar{nullptr};
184113

185114
tar = fW_data;
186115
if (!tar)
187116
return;
188117

189-
for (int nHarm{2}; nHarm <= 3; nHarm++) {
190-
fW_data->Add(new TH1F(Form("havgQ%s%i_FT0C", xy, nHarm), "", CentBin, 0, CentBin));
191-
TH1F* hRec = reinterpret_cast<TH1F*>(tar->FindObject(Form("havgQ%s%i_FT0C", xy, nHarm)));
192-
if (!hRec)
193-
return;
194-
195-
TH2F* hQ = reinterpret_cast<TH2F*>(fW_data->At(0)->Clone(Form("hQ%s%i_FT0C", xy, nHarm)));
196-
TH1D* h1;
197-
198-
for (int i{1}; i < hQ->GetXaxis()->GetNbins() + 1; i++) {
199-
h1 = static_cast<TH1D*>(hQ->ProjectionY(Form("proj%i_Q%s%is", i, xy, nHarm), i, i));
200-
201-
double mean = h1->GetMean();
202-
double meanErr = h1->GetMeanError();
203-
int binc = hRec->GetXaxis()->GetBinCenter(i);
204-
205-
hRec->SetBinContent(binc, mean);
206-
hRec->SetBinError(binc, meanErr);
118+
for (const auto& pf : stv) {
119+
for (const auto& nh : nhv) {
120+
TH2D* th2 = reinterpret_cast<TH2D*>(tar->FindObject(this->GetQName(nh, pf.c_str())));
121+
if (!th2) {
122+
printf("qh not found!\n");
123+
return;
124+
}
125+
126+
TH1D* tmp = nullptr;
127+
TGraph* tmpgr = nullptr;
128+
TSpline3* spline = nullptr;
129+
for (int iSP{0}; iSP < 90; iSP++) {
130+
tmp = th2->ProjectionY(Form("q%i_%i_%i", nh, iSP, iSP + 1), iSP + 1, iSP + 1);
131+
std::vector<double> xq(nResolution);
132+
std::vector<double> yq(nResolution);
133+
for (int i{0}; i < nResolution; i++)
134+
xq[i] = static_cast<double>(i + 1) / static_cast<double>(nResolution);
135+
tmp->GetQuantiles(nResolution, yq.data(), xq.data());
136+
tmpgr = new TGraph(nResolution, yq.data(), xq.data());
137+
spline = new TSpline3(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP), tmpgr);
138+
spline->SetName(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP));
139+
fW_data->Add(spline);
140+
}
207141
}
208142
}
209143
};
210144

211-
void FFitWeights::CreateRMS(const char* xy)
145+
float FFitWeights::EvalSplines(float centr, const float& dqn, const int nh, const char* pf)
212146
{
213-
if (fW_data->GetEntries() < 1)
214-
return;
215-
216147
TObjArray* tar{nullptr};
217148

218149
tar = fW_data;
219150
if (!tar)
220-
return;
151+
return -1;
221152

222-
for (int nHarm{2}; nHarm <= 3; nHarm++) {
223-
fW_data->Add(new TH1F(Form("hrmsQ%s%i_FT0C", xy, nHarm), "", CentBin, 0, CentBin));
224-
TH1F* hRec = reinterpret_cast<TH1F*>(tar->FindObject(Form("hrmsQ%s%i_FT0C", xy, nHarm)));
225-
if (!hRec)
226-
return;
227-
228-
TH2F* hQ = reinterpret_cast<TH2F*>(fW_data->At(0)->Clone(Form("hQ%s%i_FT0C", xy, nHarm)));
229-
TH1D* h1;
230-
231-
for (int i{1}; i < hQ->GetXaxis()->GetNbins() + 1; i++) {
232-
h1 = static_cast<TH1D*>(hQ->ProjectionY(Form("proj%i_Q%s%is", i, xy, nHarm), i, i));
233-
234-
double mean = h1->GetRMS();
235-
double meanErr = h1->GetRMSError();
236-
int binc = hRec->GetXaxis()->GetBinCenter(i);
237-
238-
hRec->SetBinContent(binc, mean);
239-
hRec->SetBinError(binc, meanErr);
240-
}
153+
int isp = static_cast<int>(centr);
154+
if (isp < 0 || isp > 90) {
155+
return -1;
241156
}
242-
};
243-
244-
float FFitWeights::GetRecVal(int cent, const char* xy, const int nHarm)
245-
{
246-
TObjArray* tar{nullptr};
247-
248-
tar = fW_data;
249-
if (!tar)
250-
return -999;
251-
252-
TH1F* htmp = reinterpret_cast<TH1F*>(tar->FindObject(Form("havgQ%s%i_FT0C", xy, nHarm)));
253-
254-
return htmp->GetBinContent(cent);
255-
};
256157

257-
float FFitWeights::GetRMSVal(int cent, const char* xy, const int nHarm)
258-
{
259-
TObjArray* tar{nullptr};
260-
261-
tar = fW_data;
262-
if (!tar)
263-
return -999;
264-
265-
TH1F* htmp = reinterpret_cast<TH1F*>(tar->FindObject(Form("hrmsQ%s%i_FT0C", xy, nHarm)));
266-
267-
return htmp->GetBinContent(cent);
268-
};
269-
270-
TH1F* FFitWeights::GetRecHist(const char* xy, const int nHarm)
271-
{
272-
TObjArray* tar{nullptr};
273-
274-
tar = fW_data;
275-
if (!tar)
276-
return nullptr;
277-
278-
return reinterpret_cast<TH1F*>(tar->FindObject(Form("havgQ%s%i_FT0C", xy, nHarm)));
279-
};
280-
281-
TH1F* FFitWeights::GetRmsHist(const char* xy, const int nHarm)
282-
{
283-
TObjArray* tar{nullptr};
158+
TSpline3* spline = nullptr;
159+
spline = reinterpret_cast<TSpline3*>(tar->FindObject(Form("sp_q%i%s_%i", nh, pf, isp)));
160+
if (!spline) {
161+
return -1;
162+
}
284163

285-
tar = fW_data;
286-
if (!tar)
287-
return nullptr;
164+
float qn_val = 100. * spline->Eval(dqn);
165+
if (qn_val < 0) {
166+
return -1;
167+
}
288168

289-
return reinterpret_cast<TH1F*>(tar->FindObject(Form("hrmsQ%s%i_FT0C", xy, nHarm)));
169+
return qn_val;
290170
};
291-
292-
// float FFitWeights::EventPlane(const float& x, const float& y, const float& nHarm) {
293-
// return 1/nHarm * TMath::ATan2(y, x);
294-
// };

0 commit comments

Comments
 (0)