Belle II Software development
eclComputePulseTemplates_Step3.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <TF1.h>
9#include <math.h>
10#include <TTree.h>
11#include <TAxis.h>
12#include <TFile.h>
13#include <TCanvas.h>
14#include <vector>
15#include <TGraph.h>
16#include <iostream>
17#include <ecl/digitization/OfflineFitFunction.h>
18#include <assert.h>
19
20//
21// See eclComputePulseTemplates_Step0.cc for README instructions.
22//
23//Check Maximum residual for fit to template shape.
24double GetMaxRes(TGraph* gin, TF1* fit)
25{
26 double MaxResidual = 0;
27 double MaxData = 0;
28 for (int j = 0; j < 1000; j++) {
29 if (gin->GetX()[j] < 10) {
30 double tempVal = gin->GetY()[j];
31 double tempResidual = tempVal - fit->Eval(gin->GetX()[j]);
32 if (gin->GetY()[j] > MaxData) MaxData = gin->GetY()[j];
33 if (abs(tempResidual) > abs(MaxResidual)) MaxResidual = tempResidual;
34 }
35 }
36 return abs(MaxResidual) / MaxData;
37}
38//
39TF1* FitPulse(TGraph* gin, int ShapeFlag, double* pulseInputPara)
40{
41 //ShapeFlag = 0 for hadron
42 //ShapeFlag = 1 for diode
43 //
44 double highRange = 10;
45 TF1* ShpFloat = new TF1(Form("Shp_%d", ShapeFlag), Belle2::ECL::WaveFuncTwoComponent, 0, highRange, 26);
46 ShpFloat->FixParameter(0, 0);
47 ShpFloat->FixParameter(1, 0);
48 ShpFloat->FixParameter(2, 1);
49 ShpFloat->FixParameter(3, 0);
50 for (int k = 0; k < 10; k++) {
51 ShpFloat->SetParameter(4 + k, pulseInputPara[k + 1]);
52 ShpFloat->FixParameter(10 + 4 + k, pulseInputPara[k + 1]);
53 }
54 ShpFloat->SetParameter(24, pulseInputPara[0]);
55 ShpFloat->FixParameter(25, 1);
56 ShpFloat->SetNpx(1000);
57 gin->Fit(Form("Shp_%d", ShapeFlag), "Q N 0 R", "", 0, 30);
58 //
59 //Try Fit and Check residual
60 double Check = GetMaxRes(gin, ShpFloat);
61 double CheckMin = Check;
62 double ParMin11[11];
63 ParMin11[0] = ShpFloat->GetParameter(24);
64 for (int k = 0; k < 10; k++)ParMin11[k + 1] = ShpFloat->GetParameter(k + 4);
65 int nFits = 0;
66 int Attempt = 0;
67 //keep fitting until residual is less than 1%
68 std::cout << nFits << " " << Check << std::endl;
69 //while (false){
70 while ((Check > 0.01) && nFits < 20) {
71 //
72 ShpFloat->SetParameter(24, ParMin11[0] * (1.0 + (0.01 * nFits)));
73 for (int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1] * (1.0 + (0.01 * nFits)));
74 //
75 gin->Fit(Form("Shp_%d", ShapeFlag), "Q N 0 R W", "", 0.01, highRange);
76 nFits++;
77 Check = GetMaxRes(gin, ShpFloat);
78 //
79 if (Check < CheckMin) {
80 //save neew best parameters:w
81 CheckMin = Check;
82 ParMin11[0] = ShpFloat->GetParameter(24);
83 for (int k = 0; k < 10; k++)ParMin11[k + 1] = ShpFloat->GetParameter(k + 4);
84 }
85 std::cout << "nFit=" << nFits << " " << Check << std::endl;;
86 //
87 if (nFits == 20 && Attempt == 0 && (Check > 0.01)) {
88 //Try new initial conditions
89 if (ShapeFlag == 1) {
90 const double ParMin11t[11] = {36.1232, -0.284876, 0.350343, 0.432839, 0.445749, 0.27693, 0.00899611, 6.11111, 0.788569, 0.570159, -0.411252};
91 for (int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
92 } else {
93 const double ParMin11t[11] = {10, 0.031, 4.2e-5, 0.74, 0.43, 0.61, 0.03, 3.8, 0.81, 0.77, 0.59};
94 for (int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
95 }
96 Attempt++;
97 nFits = 0;
98 }
99 }
100 //set to best parameters
101 ShpFloat->SetParameter(24, ParMin11[0]);
102 for (int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1]);
103 //
104 return ShpFloat;
105}
106//
107// cppcheck-suppress constParameter ; main() keeps its standard signature
108int main(int argc, char* argv[])
109{
110 //
111 TString OutputDirectory = "./";
112 if (OutputDirectory == "") {
113 std::cout << "Error set output directory" << std::endl;
114 return -1;
115 }
116 //
117 assert(argc == 3);
118 int LowCellID = atoi(argv[1]);
119 int HighCellID = atoi(argv[2]);
120 //
121 double TimeAll_A[1000];
122 double PhotonValues_A[1000];
123 double HadronValues_A[1000];
124 double DiodeValues_A[1000];
125 TFile* HadronShapeFile = new TFile(OutputDirectory + Form("HadronShapes_Low%d_High%d.root", LowCellID, HighCellID));
126 TTree* chain = static_cast<TTree*>(HadronShapeFile->Get("HadronTree"));
127 chain->SetBranchAddress("TimeAll_A", &TimeAll_A);
128 chain->SetBranchAddress("ValuePhoton_A", &PhotonValues_A);
129 chain->SetBranchAddress("ValueHadron_A", &HadronValues_A);
130 chain->SetBranchAddress("ValueDiode_A", &DiodeValues_A);
131 //
132 double TempHadronPar11_A[11];
133 double TempDiodePar11_A[11];
134 double MaxResDiode_A = 0;
135 double MaxValDiode_A = 0;
136 double MaxResHadron_A = 0;
137 double MaxValHadron_A = 0;
138 //
139 TFile* f = new TFile(OutputDirectory + Form("HadronPars_Low%d_High%d.root", LowCellID, HighCellID), "RECREATE");
140 TTree* WaveformParametersTree = new TTree("HadronWaveformInfo", "");
141 WaveformParametersTree->Branch("TempHadronPar11_A", &TempHadronPar11_A, "TempHadronPar11_A[11]/D");
142 WaveformParametersTree->Branch("TempDiodePar11_A", &TempDiodePar11_A, "TempDiodePar11_A[11]/D");
143 WaveformParametersTree->Branch("MaxResDiode_A", &MaxResDiode_A, "MaxResDiode_A/D");
144 WaveformParametersTree->Branch("MaxResHadron_A", &MaxResHadron_A, "MaxResHadron_A/D");
145 WaveformParametersTree->Branch("MaxValHadron_A", &MaxValHadron_A, "MaxValHadron_A/D");
146 WaveformParametersTree->Branch("MaxValDiode_A", &MaxValDiode_A, "MaxValDiode_A/D");
147 std::vector<TCanvas*> AllFits;
148 //
149 Long64_t nentries = chain->GetEntriesFast();
150 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
151 chain->GetEntry(jentry);
152 //
153 for (int l = 0; l < 11; l++) {
154 TempDiodePar11_A[l] = 0;
155 TempHadronPar11_A[l] = 0;
156 }
157 int CellID = ((int)jentry) + LowCellID;
158 //
159 if (TimeAll_A[0] > -100) {
160 //
161 TGraph* G1_A = new TGraph(1000, TimeAll_A, PhotonValues_A);
162 TGraph* G2_A = new TGraph(1000, TimeAll_A, HadronValues_A);
163 TGraph* G3_A = new TGraph(1000, TimeAll_A, DiodeValues_A);
164 //
165 double InputPara_A_H[11] = {1.501407e+01, 0.000000e+00, 4.693325e-03, 7.538656e-01, 4.371189e-01, 1.182754e+00, 3.242071e-03, 6.669907e+00, 8.454105e-01, 7.240517e-01, 4.097126e-01};
166 double InputPara_A_D[11] = {1.501407e+01, 0.000000e+00, 4.693325e-03, 7.538656e-01, 4.371189e-01, 1.182754e+00, 3.242071e-03, 6.669907e+00, 8.454105e-01, 7.240517e-01, 4.097126e-01};
167 //
168 TF1* F2_A = FitPulse(G2_A, 0, InputPara_A_H);
169 TF1* F3_A = FitPulse(G3_A, 1, InputPara_A_D);
170 //
171 TempHadronPar11_A[0] = F2_A->GetParameter(24);
172 TempDiodePar11_A[0] = F3_A->GetParameter(24);
173 for (int k = 0; k < 10; k++) {
174 TempHadronPar11_A[k + 1] = F2_A->GetParameter(k + 4);
175 TempDiodePar11_A[k + 1] = F3_A->GetParameter(k + 4);
176 }
177 //
178 MaxResHadron_A = GetMaxRes(G2_A, F2_A);
179 MaxResDiode_A = GetMaxRes(G3_A, F3_A);
180 //
181 MaxValHadron_A = F2_A->GetMaximum(1, 10);
182 MaxValDiode_A = F3_A->GetMaximum(1, 10);
183 //
184 std::cout << CellID << " " << MaxResHadron_A << " " << MaxResDiode_A << std::endl;
185 WaveformParametersTree->Fill();
186 //
187 //Set to true to draw fit results in output file
188 //if (false) {
189 if (true) {
190 G2_A->SetMarkerColor(kBlue);
191 G3_A->SetMarkerColor(kRed);
192 G1_A->GetYaxis()->SetRangeUser(-0.5, 2);
193 //
194 AllFits.push_back(new TCanvas(Form("A_all_%lld", jentry), ""));
195 G1_A->Draw("AP");
196 G2_A->Draw("same P");
197 G3_A->Draw("sameP");
198 AllFits.push_back(new TCanvas(Form("A_fitHad_%lld", jentry), ""));
199 AllFits[AllFits.size() - 1]->cd();
200 G2_A->Draw("AP");
201 F2_A->Draw("same l");
202 AllFits.push_back(new TCanvas(Form("A_fitDio_%lld", jentry), ""));
203 AllFits[AllFits.size() - 1]->cd();
204 G3_A->Draw("AP");
205 F3_A->Draw("same l");
206 AllFits.push_back(new TCanvas(Form("MC_all_%lld", jentry), ""));
207 G1_A->Draw("AP");
208 //TF1* MCgamam = new TF1(Form("MC_%lld", jentry), Wa, 0, 32, 11);
209 //MCgamam->SetParameter(0, 27.7221);
210 //MCgamam->SetParameter(1, 0.);
211 //MCgamam->SetParameter(2, 0.648324);
212 //MCgamam->SetParameter(3, 0.401711);
213 //MCgamam->SetParameter(4, 0.374167);
214 //MCgamam->SetParameter(5, 0.849417);
215 //MCgamam->SetParameter(6, 0.00144548);
216 //MCgamam->SetParameter(7, 4.70722);
217 //MCgamam->SetParameter(8, 0.815639);
218 //MCgamam->SetParameter(9, 0.555605);
219 //MCgamam->SetParameter(10, 0.2752);
220 //MCgamam->Draw("same l");
221 } else {
222 G1_A->Delete();
223 G2_A->Delete();
224 G3_A->Delete();
225 F2_A->Delete();
226 F3_A->Delete();
227 }
228 }
229 }
230 //
231 std::cout << "eclComputePulseTemplates_Step3: Writing File " << f->GetName() << std::endl;
232 f->cd();
233 WaveformParametersTree->Write();
234 for (unsigned int k = 0; k < AllFits.size(); k++) AllFits[k]->Write();
235 f->Write();
236 //
237 return 1;
238 //
239}