Belle II Software  release-05-01-25
eclComputePulseTemplates_Step0.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Sample scripts to compute hadron and diode templates *
6  * *
7  * Author: The Belle II Collaboration *
8  * Contributors: Savino Longo (longos@uvuc.ca) *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
13 #include<iostream>
14 #include<fstream>
15 #include<stdlib.h>
16 #include"TTree.h"
17 #include"TString.h"
18 #include"TFile.h"
19 #include <assert.h>
20 
21 using namespace std;
22 
23 /*
24 Note photon shape calibrations are needed to compute hadron and diode templates.
25 Photon calibrations should be placed in file named params_gamma_shape.dat in same directiory as this script.
26 Format of params_gamma_shape.dat:
27 
28 cellID PhotonScale 0 p0 p1 p2 p3 p4 p5 p6 p7 p8 p9
29 
30 where p0-p9 are 10 ShaperDSP shaper parameters for basf2 C++ implementation of ShaperDSP function.
31 
32 eclComputePulseTemplates_Step0.cc - places photon parameters in ntuple
33 eclComputePulseTemplates_Step1.cc - produce photon array input for eclComputePulseTemplates_Step2.py
34 eclComputePulseTemplates_Step2.py - produce hadron and diode array input for eclComputePulseTemplates_Step3.py
35 eclComputePulseTemplates_Step3.cc - computed hadron and diode shape parameters and organizes into temp ntuple
36 eclComputePulseTemplates_Step4.cc - organizes photon, hadron and diode templates into final ntuple to produce dbobject.
37 
38 To run template calculation:
39 
40 First "OutputDirectory" variable in eclComputePulseTemplates_Stepx.cc files should be modified to a temp directory before running scripts.
41 recompile with "scons ecl"
42 
43 1) execute eclComputePulseTemplates_Step0 0
44 2) for i=0 i<873 i++
45 3) execute eclComputePulseTemplates_Step1 i*10 (i+1)*10
46 4) basf2 eclComputePulseTemplates_Step2.py i*10 (i+1)*10
47 5) execute eclComputePulseTemplates_Step3 i*10 (i+1)*10
48 6) execute eclComputePulseTemplates_Step0 1
49 7) execute eclWriteWaveformParametersLocalDB
50 Final step will create local db object.
51 */
52 
54 struct crystalInfo {
55  vector<double> PhotonWaveformPars;
56  vector<double> HadronWaveformPars;
57  vector<double> DiodeWaveformPars;
58  double MaxResHadron;
59  double MaxValHadron;
60  double MaxResDiode;
61  double MaxValDiode;
62 };
63 
64 int main(int argc, char* argv[])
65 {
66  //
67  TString OutputDirectory = "";
68  if (OutputDirectory == "") {
69  std::cout << "Error set ouput directory" << std::endl;
70  return -1;
71  }
72  //
73  assert(argc == 2);
74  int Flag = atoi(argv[1]);
75  TString InputDirectory = OutputDirectory;
76  //
77  TString ParameterTreeOutputName = OutputDirectory + "PhotonWaveformParameters.root";
78  if (Flag == 1) ParameterTreeOutputName = "DigitWaveformParameters.root";
79  TFile* ParameterTreeOutputFile = new TFile(ParameterTreeOutputName, "RECREATE");
80  //
81  TTree* ParameterTree = new TTree("ParTree", "");
82  double PhotonWaveformPar[11];
83  double HadronWaveformPar[11];
84  double DiodeWaveformPar[11];
85  for (int k = 0; k < 11; k++) {
86  PhotonWaveformPar[k] = 0;
87  HadronWaveformPar[k] = 0;
88  DiodeWaveformPar[k] = 0;
89  }
90  double mHadronRes = 0;
91  double mDiodeRes = 0;
92  double mHadronMax = 0;
93  double mDiodeMax = 0;
94  ParameterTree->Branch("PhotonPar", &PhotonWaveformPar, "PhotonWaveformPar[11]/D");
95  ParameterTree->Branch("HadronPar", &HadronWaveformPar, "HadronWaveformPar[11]/D");
96  ParameterTree->Branch("DiodePar", &DiodeWaveformPar, "DiodeWaveformPar[11]/D");
97  ParameterTree->Branch("mHadronRes", &mHadronRes, "mHadronRes/D");
98  ParameterTree->Branch("mDiodeRes", &mDiodeRes, "mDiodeRes/D");
99  ParameterTree->Branch("mHadronMax", &mHadronMax, "mHadronMax/D");
100  ParameterTree->Branch("mDiodeMax", &mDiodeMax, "mDiodeMax/D");
101  //
102  std::vector<crystalInfo> cellIDcheck(8736);
103  ifstream PhotonFile("/home/belle2/longos/WaveformFitting/ecl/tools/params_gamma_shape.dat");
104  if (PhotonFile.is_open()) {
105  std::vector<double> templine(12);
106  for (int k = 0; k < 8736; k++) {
107  for (unsigned int j = 0; j < templine.size(); j++) PhotonFile >> templine[j];
108  cellIDcheck[int(templine[0]) - 1].PhotonWaveformPars.resize(11);
109  cellIDcheck[int(templine[0]) - 1].PhotonWaveformPars[0] = templine[1];
110  for (int j = 0; j < 10; j++) cellIDcheck[int(templine[0]) - 1].PhotonWaveformPars[j + 1] = templine[j + 2];
111  std::cout << int(templine[0]) << " " << templine[1] << endl;
112  }
113  PhotonFile.close();
114  } else {
115  std::cout << "ERROR cannont open photon param file." << endl;
116  return -1;
117  }
118  if (Flag == 0) {
119  for (unsigned int f = 0; f < cellIDcheck.size(); f++) {
120  for (int k = 0; k < 11; k++) PhotonWaveformPar[k] = cellIDcheck[f].PhotonWaveformPars[k];
121  ParameterTree->Fill();
122  }
123  } else if (Flag == 1) {
124  int batch = 10;
125  for (int k = 0; k < 874; k++) {
126  int low = k * batch;
127  int high = (k + 1) * batch;
128  //std::cout<<k<<std::endl;
129  TFile* TempFile = new TFile(InputDirectory + Form("HadronPars_Low%d_High%d.root", low, high), "READ");
130  TTree* TempTree = (TTree*) TempFile->Get("HadronWaveformInfo");
131  double tHadronShapePars_A[11];
132  double tDiodeShapePars_A[11];
133  double tMaxResidualHadron_A;
134  double tMaxResidualDiode_A;
135  double tMaxValDiode_A;
136  double tMaxValHadron_A;
137  TempTree->SetBranchAddress("TempHadronPar11_A", &tHadronShapePars_A);
138  TempTree->SetBranchAddress("TempDiodePar11_A", &tDiodeShapePars_A);
139  TempTree->SetBranchAddress("MaxResHadron_A", &tMaxResidualHadron_A);
140  TempTree->SetBranchAddress("MaxResDiode_A", &tMaxResidualDiode_A);
141  TempTree->SetBranchAddress("MaxValDiode_A", &tMaxValDiode_A);
142  TempTree->SetBranchAddress("MaxValHadron_A", &tMaxValHadron_A);
143  //
144  for (int j = 0; j < batch; j++) {
145  int tCellID = (k * batch) + j;
146  if (tCellID >= 8736)continue;
147  TempTree->GetEntry(j);
148  cellIDcheck[tCellID].HadronWaveformPars.resize(11);
149  cellIDcheck[tCellID].DiodeWaveformPars.resize(11);
150  for (int p = 0; p < 11; p++) {
151  if (tHadronShapePars_A[p] > 100 || tHadronShapePars_A[p] < -100) {
152  std::cout << "Warning large parameter for: " << tCellID << " " << tHadronShapePars_A[p] << std::endl;
153  for (int h = 0; h < 11; h++)std::cout << tHadronShapePars_A[h] << " , ";
154  std::cout << std::endl;
155  }
156  cellIDcheck[tCellID].HadronWaveformPars[p] = tHadronShapePars_A[p];
157  }
158  //std::cout<<tCellID<<std::endl;
159  //for(int h=0;h<11;h++ )std::cout<<tHadronShapePars_A[h]<<" , ";
160  //std::cout<<std::endl;
161  cellIDcheck[tCellID].MaxResHadron = tMaxResidualHadron_A;
162  cellIDcheck[tCellID].MaxValHadron = tMaxValHadron_A;
163  for (int p = 0; p < 11; p++) cellIDcheck[tCellID].DiodeWaveformPars[p] = tDiodeShapePars_A[p];
164  cellIDcheck[tCellID].MaxResDiode = tMaxResidualDiode_A;
165  cellIDcheck[tCellID].MaxValDiode = tMaxValDiode_A;
166  }
167  TempFile->Close();
168  }
169  //
170  for (unsigned int f = 0; f < 8736; f++) {
171  for (int k = 0; k < 11; k++) {
172  PhotonWaveformPar[k] = cellIDcheck[f].PhotonWaveformPars[k];
173  HadronWaveformPar[k] = cellIDcheck[f].HadronWaveformPars[k];
174  DiodeWaveformPar[k] = cellIDcheck[f].DiodeWaveformPars[k];
175  }
176  mHadronRes = cellIDcheck[f].MaxResHadron;
177  mDiodeRes = cellIDcheck[f].MaxResDiode;
178  mHadronMax = cellIDcheck[f].MaxValHadron;
179  mDiodeMax = cellIDcheck[f].MaxValDiode;
180  ParameterTree->Fill();
181  }
182  }
183  //
184  ParameterTreeOutputFile->cd();
185  ParameterTree->Write();
186  ParameterTreeOutputFile->Write();
187  ParameterTreeOutputFile->Close();
188  //
189  //
190  return 0;
191 }
crystalInfo
a struct to hold relavent Crystal information
Definition: eclComputePulseTemplates_Step0.cc:54
crystalInfo::HadronWaveformPars
vector< double > HadronWaveformPars
hadron waveform parameters
Definition: eclComputePulseTemplates_Step0.cc:56
crystalInfo::MaxValHadron
double MaxValHadron
max.
Definition: eclComputePulseTemplates_Step0.cc:59
crystalInfo::MaxResHadron
double MaxResHadron
max.
Definition: eclComputePulseTemplates_Step0.cc:58
crystalInfo::DiodeWaveformPars
vector< double > DiodeWaveformPars
diode waveform parameters
Definition: eclComputePulseTemplates_Step0.cc:57
crystalInfo::PhotonWaveformPars
vector< double > PhotonWaveformPars
photon waveform parameters
Definition: eclComputePulseTemplates_Step0.cc:55
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
crystalInfo::MaxResDiode
double MaxResDiode
max.
Definition: eclComputePulseTemplates_Step0.cc:60
crystalInfo::MaxValDiode
double MaxValDiode
max.
Definition: eclComputePulseTemplates_Step0.cc:61