Belle II Software  release-06-02-00
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.
24 double 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 //
39 TF1* 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 ((Check > 0.01) && nFits < 20) {
70  //
71  ShpFloat->SetParameter(24, ParMin11[0] * (1.0 + (0.01 * nFits)));
72  for (int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1] * (1.0 + (0.01 * nFits)));
73  //
74  gin->Fit(Form("Shp_%d", ShapeFlag), "Q N 0 R W", "", 0.01, highRange);
75  nFits++;
76  Check = GetMaxRes(gin, ShpFloat);
77  //
78  if (Check < CheckMin) {
79  //save neew best parameters:w
80  CheckMin = Check;
81  ParMin11[0] = ShpFloat->GetParameter(24);
82  for (int k = 0; k < 10; k++)ParMin11[k + 1] = ShpFloat->GetParameter(k + 4);
83  }
84  std::cout << "nFit=" << nFits << " " << Check << std::endl;;
85  //
86  if (nFits == 20 && Attempt == 0 && (Check > 0.01)) {
87  //Try new initial conditions
88  if (ShapeFlag == 1) {
89  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};
90  for (int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
91  } else {
92  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};
93  for (int k = 0; k < 11; k++)ParMin11[k] = ParMin11t[k];
94  }
95  Attempt++;
96  nFits = 0;
97  }
98  }
99  //set to best parameters
100  ShpFloat->SetParameter(24, ParMin11[0]);
101  for (int k = 0; k < 10; k++) ShpFloat->SetParameter(k + 4, ParMin11[k + 1]);
102  //
103  return ShpFloat;
104 }
105 //
106 int main(int argc, char* argv[])
107 {
108  //
109  TString OutputDirectory = "";
110  if (OutputDirectory == "") {
111  std::cout << "Error set ouput directory" << std::endl;
112  return -1;
113  }
114  //
115  assert(argc == 3);
116  int LowCellID = atoi(argv[1]);
117  int HighCellID = atoi(argv[2]);
118  //
119  double TimeAll_A[1000];
120  double PhotonValues_A[1000];
121  double HadronValues_A[1000];
122  double DiodeValues_A[1000];
123  TFile* HadronShapeFile = new TFile(OutputDirectory + Form("HadronShapes_Low%d_High%d.root", LowCellID, HighCellID));
124  TTree* chain = (TTree*) HadronShapeFile->Get("HadronTree");
125  chain->SetBranchAddress("TimeAll_A", &TimeAll_A);
126  chain->SetBranchAddress("ValuePhoton_A", &PhotonValues_A);
127  chain->SetBranchAddress("ValueHadron_A", &HadronValues_A);
128  chain->SetBranchAddress("ValueDiode_A", &DiodeValues_A);
129  //
130  double TempHadronPar11_A[11];
131  double TempDiodePar11_A[11];
132  double MaxResDiode_A = 0;
133  double MaxValDiode_A = 0;
134  double MaxResHadron_A = 0;
135  double MaxValHadron_A = 0;
136  //
137  TFile* f = new TFile(OutputDirectory + Form("HadronPars_Low%d_High%d.root", LowCellID, HighCellID), "RECREATE");
138  TTree* WaveformParametersTree = new TTree("HadronWaveformInfo", "");
139  WaveformParametersTree->Branch("TempHadronPar11_A", &TempHadronPar11_A, "TempHadronPar11_A[11]/D");
140  WaveformParametersTree->Branch("TempDiodePar11_A", &TempDiodePar11_A, "TempDiodePar11_A[11]/D");
141  WaveformParametersTree->Branch("MaxResDiode_A", &MaxResDiode_A, "MaxResDiode_A/D");
142  WaveformParametersTree->Branch("MaxResHadron_A", &MaxResHadron_A, "MaxResHadron_A/D");
143  WaveformParametersTree->Branch("MaxValHadron_A", &MaxValHadron_A, "MaxValHadron_A/D");
144  WaveformParametersTree->Branch("MaxValDiode_A", &MaxValDiode_A, "MaxValDiode_A/D");
145  std::vector<TCanvas*> AllFits;
146  //
147  Long64_t nentries = chain->GetEntriesFast();
148  for (Long64_t jentry = 0; jentry < nentries; jentry++) {
149  chain->GetEntry(jentry);
150  //
151  for (int l = 0; l < 11; l++) {
152  TempDiodePar11_A[l] = 0;
153  TempHadronPar11_A[l] = 0;
154  }
155  int CellID = ((int)jentry) + LowCellID + 1;
156  //
157  if (TimeAll_A[0] > -100) {
158  //
159  TGraph* G1_A = new TGraph(1000, TimeAll_A, PhotonValues_A);
160  TGraph* G2_A = new TGraph(1000, TimeAll_A, HadronValues_A);
161  TGraph* G3_A = new TGraph(1000, TimeAll_A, DiodeValues_A);
162  //
163  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};
164  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};
165  //
166  TF1* F2_A = FitPulse(G2_A, 0, InputPara_A_H);
167  TF1* F3_A = FitPulse(G3_A, 1, InputPara_A_D);
168  //
169  TempHadronPar11_A[0] = F2_A->GetParameter(24);
170  TempDiodePar11_A[0] = F3_A->GetParameter(24);
171  for (int k = 0; k < 10; k++) {
172  TempHadronPar11_A[k + 1] = F2_A->GetParameter(k + 4);
173  TempDiodePar11_A[k + 1] = F3_A->GetParameter(k + 4);
174  }
175  //
176  MaxResHadron_A = GetMaxRes(G2_A, F2_A);
177  MaxResDiode_A = GetMaxRes(G3_A, F3_A);
178  //
179  MaxValHadron_A = F2_A->GetMaximum(1, 10);
180  MaxValDiode_A = F3_A->GetMaximum(1, 10);
181  //
182  std::cout << CellID << " " << MaxResHadron_A << " " << MaxResDiode_A << std::endl;
183  WaveformParametersTree->Fill();
184  //
185  //Set to true to draw fit results in output file
186  if (false) {
187  G2_A->SetMarkerColor(kBlue);
188  G3_A->SetMarkerColor(kRed);
189  G1_A->GetYaxis()->SetRangeUser(-0.5, 2);
190  //
191  AllFits.push_back(new TCanvas(Form("A_all_%lld", jentry), ""));
192  G1_A->Draw("AP");
193  G2_A->Draw("same P");
194  G3_A->Draw("sameP");
195  AllFits.push_back(new TCanvas(Form("A_fitHad_%lld", jentry), ""));
196  AllFits[AllFits.size() - 1]->cd();
197  G2_A->Draw("AP");
198  F2_A->Draw("same l");
199  AllFits.push_back(new TCanvas(Form("A_fitDio_%lld", jentry), ""));
200  AllFits[AllFits.size() - 1]->cd();
201  G3_A->Draw("AP");
202  F3_A->Draw("same l");
203  AllFits.push_back(new TCanvas(Form("MC_all_%lld", jentry), ""));
204  G1_A->Draw("AP");
205  //TF1* MCgamam = new TF1(Form("MC_%lld", jentry), Wa, 0, 32, 11);
206  //MCgamam->SetParameter(0, 27.7221);
207  //MCgamam->SetParameter(1, 0.);
208  //MCgamam->SetParameter(2, 0.648324);
209  //MCgamam->SetParameter(3, 0.401711);
210  //MCgamam->SetParameter(4, 0.374167);
211  //MCgamam->SetParameter(5, 0.849417);
212  //MCgamam->SetParameter(6, 0.00144548);
213  //MCgamam->SetParameter(7, 4.70722);
214  //MCgamam->SetParameter(8, 0.815639);
215  //MCgamam->SetParameter(9, 0.555605);
216  //MCgamam->SetParameter(10, 0.2752);
217  //MCgamam->Draw("same l");
218  } else {
219  G1_A->Delete();
220  G2_A->Delete();
221  G3_A->Delete();
222  F2_A->Delete();
223  F3_A->Delete();
224  }
225  }
226  }
227  //
228  f->cd();
229  WaveformParametersTree->Write();
230  for (unsigned int k = 0; k < AllFits.size(); k++) AllFits[k]->Write();
231  f->Write();
232  //
233  return 1;
234  //
235 }
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:75