10 #include <ecl/digitization/OfflineFitFunction.h>
17 double GetMaxRes(TGraph* gin, TF1* fit)
19 double MaxResidual = 0;
21 for (
int j = 0; j < 1000; j++) {
22 if (gin->GetX()[j] < 10) {
23 double tempVal = gin->GetY()[j];
24 double tempResidual = tempVal - fit->Eval(gin->GetX()[j]);
25 if (gin->GetY()[j] > MaxData) MaxData = gin->GetY()[j];
26 if (abs(tempResidual) > abs(MaxResidual)) MaxResidual = tempResidual;
29 return abs(MaxResidual) / MaxData;
32 TF1* FitPulse(TGraph* gin,
int ShapeFlag,
double* pulseInputPara)
37 double highRange = 10;
38 TF1* ShpFloat =
new TF1(Form(
"Shp_%d", ShapeFlag), Belle2::ECL::WaveFuncTwoComponent, 0, highRange, 26);
39 ShpFloat->FixParameter(0, 0);
40 ShpFloat->FixParameter(1, 0);
41 ShpFloat->FixParameter(2, 1);
42 ShpFloat->FixParameter(3, 0);
43 for (
int k = 0; k < 10; k++) {
44 ShpFloat->SetParameter(4 + k, pulseInputPara[k + 1]);
45 ShpFloat->FixParameter(10 + 4 + k, pulseInputPara[k + 1]);
47 ShpFloat->SetParameter(24, pulseInputPara[0]);
48 ShpFloat->FixParameter(25, 1);
49 ShpFloat->SetNpx(1000);
50 gin->Fit(Form(
"Shp_%d", ShapeFlag),
"Q N 0 R",
"", 0, 30);
53 double Check = GetMaxRes(gin, ShpFloat);
54 double CheckMin = Check;
56 ParMin11[0] = ShpFloat->GetParameter(24);
57 for (
int k = 0; k < 10; k++)ParMin11[k + 1] = ShpFloat->GetParameter(k + 4);
61 std::cout << nFits <<
" " << Check << std::endl;
62 while ((Check > 0.01) && nFits < 20) {
64 ShpFloat->SetParameter(24, ParMin11[0] * (1.0 + (0.01 * nFits)));
65 for (
int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1] * (1.0 + (0.01 * nFits)));
67 gin->Fit(Form(
"Shp_%d", ShapeFlag),
"Q N 0 R W",
"", 0.01, highRange);
69 Check = GetMaxRes(gin, ShpFloat);
71 if (Check < CheckMin) {
74 ParMin11[0] = ShpFloat->GetParameter(24);
75 for (
int k = 0; k < 10; k++)ParMin11[k + 1] = ShpFloat->GetParameter(k + 4);
77 std::cout <<
"nFit=" << nFits <<
" " << Check << std::endl;;
79 if (nFits == 20 && Attempt == 0 && (Check > 0.01)) {
82 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};
83 for (
int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
85 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};
86 for (
int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
93 ShpFloat->SetParameter(24, ParMin11[0]);
94 for (
int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1]);
99 int main(
int argc,
char* argv[])
102 TString OutputDirectory =
"";
103 if (OutputDirectory ==
"") {
104 std::cout <<
"Error set ouput directory" << std::endl;
109 int LowCellID = atoi(argv[1]);
110 int HighCellID = atoi(argv[2]);
112 double TimeAll_A[1000];
113 double PhotonValues_A[1000];
114 double HadronValues_A[1000];
115 double DiodeValues_A[1000];
116 TFile* HadronShapeFile =
new TFile(OutputDirectory + Form(
"HadronShapes_Low%d_High%d.root", LowCellID, HighCellID));
117 TTree* chain = (TTree*) HadronShapeFile->Get(
"HadronTree");
118 chain->SetBranchAddress(
"TimeAll_A", &TimeAll_A);
119 chain->SetBranchAddress(
"ValuePhoton_A", &PhotonValues_A);
120 chain->SetBranchAddress(
"ValueHadron_A", &HadronValues_A);
121 chain->SetBranchAddress(
"ValueDiode_A", &DiodeValues_A);
123 double TempHadronPar11_A[11];
124 double TempDiodePar11_A[11];
125 double MaxResDiode_A = 0;
126 double MaxValDiode_A = 0;
127 double MaxResHadron_A = 0;
128 double MaxValHadron_A = 0;
130 TFile* f =
new TFile(OutputDirectory + Form(
"HadronPars_Low%d_High%d.root", LowCellID, HighCellID),
"RECREATE");
131 TTree* WaveformParametersTree =
new TTree(
"HadronWaveformInfo",
"");
132 WaveformParametersTree->Branch(
"TempHadronPar11_A", &TempHadronPar11_A,
"TempHadronPar11_A[11]/D");
133 WaveformParametersTree->Branch(
"TempDiodePar11_A", &TempDiodePar11_A,
"TempDiodePar11_A[11]/D");
134 WaveformParametersTree->Branch(
"MaxResDiode_A", &MaxResDiode_A,
"MaxResDiode_A/D");
135 WaveformParametersTree->Branch(
"MaxResHadron_A", &MaxResHadron_A,
"MaxResHadron_A/D");
136 WaveformParametersTree->Branch(
"MaxValHadron_A", &MaxValHadron_A,
"MaxValHadron_A/D");
137 WaveformParametersTree->Branch(
"MaxValDiode_A", &MaxValDiode_A,
"MaxValDiode_A/D");
138 std::vector<TCanvas*> AllFits;
140 Long64_t nentries = chain->GetEntriesFast();
141 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
142 chain->GetEntry(jentry);
144 for (
int l = 0; l < 11; l++) {
145 TempDiodePar11_A[l] = 0;
146 TempHadronPar11_A[l] = 0;
148 int CellID = ((int)jentry) + LowCellID + 1;
150 if (TimeAll_A[0] > -100) {
152 TGraph* G1_A =
new TGraph(1000, TimeAll_A, PhotonValues_A);
153 TGraph* G2_A =
new TGraph(1000, TimeAll_A, HadronValues_A);
154 TGraph* G3_A =
new TGraph(1000, TimeAll_A, DiodeValues_A);
156 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};
157 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};
159 TF1* F2_A = FitPulse(G2_A, 0, InputPara_A_H);
160 TF1* F3_A = FitPulse(G3_A, 1, InputPara_A_D);
162 TempHadronPar11_A[0] = F2_A->GetParameter(24);
163 TempDiodePar11_A[0] = F3_A->GetParameter(24);
164 for (
int k = 0; k < 10; k++) {
165 TempHadronPar11_A[k + 1] = F2_A->GetParameter(k + 4);
166 TempDiodePar11_A[k + 1] = F3_A->GetParameter(k + 4);
169 MaxResHadron_A = GetMaxRes(G2_A, F2_A);
170 MaxResDiode_A = GetMaxRes(G3_A, F3_A);
172 MaxValHadron_A = F2_A->GetMaximum(1, 10);
173 MaxValDiode_A = F3_A->GetMaximum(1, 10);
175 std::cout << CellID <<
" " << MaxResHadron_A <<
" " << MaxResDiode_A << std::endl;
176 WaveformParametersTree->Fill();
180 G2_A->SetMarkerColor(kBlue);
181 G3_A->SetMarkerColor(kRed);
182 G1_A->GetYaxis()->SetRangeUser(-0.5, 2);
184 AllFits.push_back(
new TCanvas(Form(
"A_all_%lld", jentry),
""));
186 G2_A->Draw(
"same P");
188 AllFits.push_back(
new TCanvas(Form(
"A_fitHad_%lld", jentry),
""));
189 AllFits[AllFits.size() - 1]->cd();
191 F2_A->Draw(
"same l");
192 AllFits.push_back(
new TCanvas(Form(
"A_fitDio_%lld", jentry),
""));
193 AllFits[AllFits.size() - 1]->cd();
195 F3_A->Draw(
"same l");
196 AllFits.push_back(
new TCanvas(Form(
"MC_all_%lld", jentry),
""));
222 WaveformParametersTree->Write();
223 for (
unsigned int k = 0; k < AllFits.size(); k++) AllFits[k]->Write();