BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx File Reference
#include "TCanvas.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include <sstream>
#include <iostream>
#include <string>
#include <cstring>
#include <vector>
#include <cmath>
#include "curve.h"

Go to the source code of this file.

Macros

#define mass_e   0.511e-3
 
#define mass_mu   105.658e-3
 
#define mass_pi   139.570e-3
 
#define mass_k   493.677e-3
 
#define mass_p   938.272e-3
 
#define P_electron_min   421.97*mass_e
 
#define P_electron_max   3030*mass_e
 
#define P_electron_bin   80
 
#define P_muon_min   10.98001*mass_mu
 
#define P_muon_max   115.6*mass_mu
 
#define P_muon_bin   20
 
#define P_pion_min   0.8
 
#define P_pion_max   0.805
 
#define P_pion_bin   1
 
#define P_kaon_min   0.63001*mass_k
 
#define P_kaon_max   0.74*mass_k
 
#define P_kaon_bin   5
 
#define P_proton_min   0.26*mass_p
 
#define P_proton_max   0.63*mass_p
 
#define P_proton_bin   320
 

Functions

void histgen (string, int, TObjArray &, vector< double > &)
 
void histps (string)
 
bool event_validate (double, double, double, int)
 
double get_hist_range (int, double, int)
 
void hist_update ()
 
const bool ps_flag (false)
 
int main ()
 

Variables

const int na = 10
 
const int nc = 2
 
const double electron_step = (double)(P_electron_max - P_electron_min) / (double)P_electron_bin
 
const double muon_step = (double)(P_muon_max - P_muon_min) / (double)P_muon_bin
 
const double pion_step = (double)(P_pion_max - P_pion_min) / (double)P_pion_bin
 
const double kaon_step = (double)(P_kaon_max - P_kaon_min) / (double)P_kaon_bin
 
const double proton_step = (double)(P_proton_max - P_proton_min) / (double)P_proton_bin
 
const int run_flag [5] = {0,0,1,0,0}
 

Macro Definition Documentation

◆ mass_e

#define mass_e   0.511e-3

◆ mass_k

#define mass_k   493.677e-3

◆ mass_mu

#define mass_mu   105.658e-3

◆ mass_p

#define mass_p   938.272e-3

◆ mass_pi

#define mass_pi   139.570e-3

◆ P_electron_bin

#define P_electron_bin   80

◆ P_electron_max

#define P_electron_max   3030*mass_e

◆ P_electron_min

#define P_electron_min   421.97*mass_e

◆ P_kaon_bin

#define P_kaon_bin   5

◆ P_kaon_max

#define P_kaon_max   0.74*mass_k

◆ P_kaon_min

#define P_kaon_min   0.63001*mass_k

◆ P_muon_bin

#define P_muon_bin   20

◆ P_muon_max

#define P_muon_max   115.6*mass_mu

◆ P_muon_min

#define P_muon_min   10.98001*mass_mu

◆ P_pion_bin

#define P_pion_bin   1

◆ P_pion_max

#define P_pion_max   0.805

◆ P_pion_min

#define P_pion_min   0.8

◆ P_proton_bin

#define P_proton_bin   320

◆ P_proton_max

#define P_proton_max   0.63*mass_p

◆ P_proton_min

#define P_proton_min   0.26*mass_p

Function Documentation

◆ event_validate()

bool event_validate ( double  bitrunc,
double  ptrk,
double  costheta,
int  type 
)

Definition at line 348 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

349{
350 bool validate = false;
351
352 switch (type) {
353 case 0:
354 //if (bitrunc > 420 / pow(ptrk, 1.2))
355 if( ( (ptrk<=0.3 && bitrunc>18842.7-102517*ptrk+157204*ptrk*ptrk) || (ptrk>0.3 && bitrunc>5186.13-14056 *ptrk +11221*ptrk*ptrk) ) && (bitrunc<32695.9-136724*ptrk+217506*ptrk*ptrk-120748*ptrk*ptrk*ptrk) )
356 validate = true;
357 break;
358 case 1:
359 /*if (bitrunc > 60 / pow(ptrk, 1.8) + 300 &&
360 bitrunc < 60 / pow(ptrk, 2.4) + 600)*/
361 validate = true;
362 break;
363 case 2:
364 //if (bitrunc < 15 / ptrk / ptrk + 540)
365 validate = true;
366 break;
367 case 3:
368 //if (bitrunc < 2000)
369 validate = true;
370 break;
371 case 4:
372 if (bitrunc > 0 && bitrunc < 1200)
373 validate = true;
374 break;
375 case 5:
376 if (bitrunc < 700 && bitrunc > 400)
377 validate = true;
378 break;
379 case 6:
380 validate = true;
381 break;
382 default:
383 std::cout << "Error: no such sample type!" << std::endl;
384 break;
385 }
386
387 if (ptrk <= 0) validate = false;
388 return validate;
389}
float ptrk

Referenced by histgen().

◆ get_hist_range()

double get_hist_range ( int  type,
double  momentum,
int  lowhigh 
)

Definition at line 391 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

391 {
392 double mean_exp[5];
393 double sigma_exp[5];
394
395 const double theta(1), Nohit(15), t0(1200);
396 dedx_pid_exp(momentum, theta, Nohit, t0, mean_exp, sigma_exp);
397
398 switch (type) {
399 case 0:
400 // cout << "type " << type << " mean " << mean_exp[4] << " sigma " << sigma_exp[4] << endl;
401 if(lowhigh) return mean_exp[4]+22*(sigma_exp[4]);
402 else return mean_exp[4]-10*(sigma_exp[4]) ? mean_exp[4]-10*(sigma_exp[4]) : 0;
403 break;
404 case 1:
405 // cout << "type " << type << " mean " << mean_exp[3] << " sigma " << sigma_exp[3] << endl;
406 if(lowhigh) return mean_exp[3]+26*(sigma_exp[3]);
407 else return mean_exp[3]-10*(sigma_exp[3]) ? mean_exp[3]-10*(sigma_exp[3]) : 0;
408 break;
409 case 2:
410 cout << "type " << type << " mean " << mean_exp[2] << " sigma " << sigma_exp[2] << endl;
411 cout << "mean - 10*sigma" << mean_exp[2]-10*(sigma_exp[2]) << endl;
412 cout << "mean + 15*sigma" << mean_exp[2]+15*(sigma_exp[2]) << endl;
413 cout << "mean + 25*sigma" << mean_exp[2]+25*(sigma_exp[2]) << endl;
414 cout << "mean + 28*sigma" << mean_exp[2]+28*(sigma_exp[2]) << endl;
415 if(lowhigh) return mean_exp[2]+28*(sigma_exp[2]);
416 else return mean_exp[2]-10*(sigma_exp[2]) ? mean_exp[2]-10*(sigma_exp[2]) : 0;
417 break;
418 case 3:
419 // cout << "type " << type << " mean " << mean_exp[1] << " sigma " << sigma_exp[1] << endl;
420 if(lowhigh) return mean_exp[1]+25*(sigma_exp[1]);
421 else return mean_exp[1]-10*(sigma_exp[1]) ? mean_exp[1]-10*(sigma_exp[1]) : 0;
422 break;
423 case 4:
424 // cout << "type " << type << " mean " << mean_exp[0] << " sigma " << sigma_exp[0] << endl;
425 if(lowhigh) return mean_exp[0]+25*(sigma_exp[0]);
426 else return mean_exp[0]-10*(sigma_exp[0]) ? mean_exp[0]-10*(sigma_exp[0]) : 0;
427 break;
428 default:
429 cout << "no type is determined" << endl;
430 return mean_exp[2] * 8;
431 }
432}
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
void dedx_pid_exp(int vflag[3], float dedx, int trkalg, int Nohit, float mom, float theta, float t0, float lsamp, double dedx_exp[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5], std::vector< double > &par, std::vector< double > &sig_par)

Referenced by histgen().

◆ hist_update()

void hist_update ( )

Definition at line 435 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

435 {
436 char hname[200];
437 TH1F *h1 = new TH1F();
438 TTree *tree = new TTree("TH1F_Col", "TH1F_Col");
439 tree->Branch("TH1F_Col", "TH1F", &h1, 3200, 0);
440 TFile* f1 = new TFile("DedxSampleHists.root", "UPDATE");
441 TTree* bin = (TTree*)f1->Get("bin");
442 Int_t totNum;
443 bin->SetBranchAddress("totalNum", &totNum);
444 bin->GetEntry(0);
445
446 for (Int_t i = 0; i < totNum; i++){
447 if(i >= 10000) sprintf(hname, "h%05d", i);
448 if(i >= 1000) sprintf(hname, "h%04d", i);
449 if(i >= 100) sprintf(hname, "h%03d", i);
450 if(i >= 10) sprintf(hname, "h%02d", i);
451 if(i < 10) sprintf(hname, "h%01d", i);
452 h1 = (TH1F*)f1->Get(hname);
453 tree->Fill();
454 }
455 tree->Write();
456 f1->Close();
457}
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
TFile * f1
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85

Referenced by main().

◆ histgen()

void histgen ( string  filename,
int  type,
TObjArray &  hist,
vector< double > &  bg 
)

Definition at line 112 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

112 {
113 if(!(run_flag[type])){cout << "skip this type " << type << " particle" << endl; return ;}
114 double exraw[100];
115 float ptrk, dEdx_meas, costheta, charge;
116 float chip, chipi, chik, chimu, chie;
117 double pmin, pmax, pstep, mass;
118 int nbin, ndedxhit;
119
120 TFile* f = new TFile(filename.c_str());
121 TTree* trk = (TTree*)f->Get("n103");
122 trk->SetBranchAddress("dEdx_meas", &dEdx_meas);
123 trk->SetBranchAddress("ptrk", &ptrk);
124 trk->SetBranchAddress("costheta", &costheta);
125 trk->SetBranchAddress("charge", &charge);
126 trk->SetBranchAddress("dEdx_hit", exraw);
127 trk->SetBranchAddress("ndedxhit", &ndedxhit);
128 trk->SetBranchAddress("chidedx_p", &chip);
129 trk->SetBranchAddress("chidedx_mu", &chimu);
130 trk->SetBranchAddress("chidedx_pi", &chipi);
131 trk->SetBranchAddress("chidedx_k", &chik);
132 trk->SetBranchAddress("chidedx_e", &chie);
133
134 //initialize variables
135 switch (type) {
136 case 0:
137 pmin = P_proton_min;
138 pmax = P_proton_max;
139 pstep = proton_step;
140 nbin = P_proton_bin;
141 mass = mass_p;
142 break;
143 case 1:
144 pmin = P_kaon_min;
145 pmax = P_kaon_max;
146 pstep = kaon_step;
147 nbin = P_kaon_bin;
148 mass = mass_k;
149 break;
150 case 2:
151 pmin = P_pion_min;
152 pmax = P_pion_max;
153 pstep = pion_step;
154 nbin = P_pion_bin;
155 mass = mass_pi;
156 break;
157 case 3:
158 pmin = P_muon_min;
159 pmax = P_muon_max;
160 pstep = muon_step;
161 nbin = P_muon_bin;
162 mass = mass_mu;
163 break;
164 case 4:
165 pmin = P_electron_min;
166 pmax = P_electron_max;
167 pstep = electron_step;
168 nbin = P_electron_bin;
169 mass = mass_e;
170 break;
171 default:
172 std::cout << "Error: no such sample type!" << std::endl;
173 }
174
175 //define histograms
176 TH1F** h = new TH1F*[nbin * 2 * na];
177 stringstream ss;
178 double momentum;
179 for (int i = 0; i < nbin * 2 * na; i++){ // nbin is determined by pmin, pmax and pstep
180 ss.str("");
181 ss << "h" << histNum + i;
182 momentum = pmin + pstep * ((int)(i / 2 / na) + 0.5);
183 h[i] = new TH1F(ss.str().c_str(), ss.str().c_str(), 200, get_hist_range(type, momentum, 0),
184 get_hist_range(type, momentum, 1));
185 }
186 histNum += nbin * 2 * na;
187
188 //loop every hit of n103
189 int binidx = -1;
190 int angleidx = -1;
191 int chargeidx = -1;
192 int histidx = -1;
193 // int histidx1 = -1;
194
195 for (int i = 0; i < trk->GetEntries(); i++){
196 trk->GetEntry(i);
197 if (type == 0 && chip > 20) continue;
198 //if (type == 1 && (chik < -4 || chik > 7)) continue;
199 //if (type == 2 && chipi > 7) continue;
200 // if (type == 3 && ptrk < 10 && chimu > 10 ) continue;
201 if (type == 3 && ptrk < 10 && dEdx_meas > 600 ) continue;
202 if (!event_validate(dEdx_meas, ptrk, costheta, type)) continue;
203
204 // betagamma bin index
205 binidx = (int)((ptrk - pmin) / pstep);
206
207 // costheta bin index
208 if (fabs(costheta) >= 0.93)
209 angleidx = na - 1;
210 else
211 angleidx = (int)(fabs(costheta) * na / 0.93);
212
213 // charge bin index
214 if (charge == 0)
215 std::cout << "Warning: Neutral track exist!" << std::endl;
216 else
217 chargeidx = (charge > 0) ? 0 : 1;
218
219 // fill histograms
220 histidx = -1;
221 for (int j = 0; j < ndedxhit; j++) {
222 if (binidx >= 0 && binidx < nbin && angleidx >= 0 && angleidx < na){
223 bool isLowMomentumProton = (type == 0 && binidx< 100);
224 bool isLowMomentumKaon = (type == 1 && binidx <= 40 );
225 bool isHighMomentumKaon = false; //(type == 1 && binidx >= 40);
226 bool isLowMomentumPion = (type == 2 && binidx < 90);
227 bool isHighMomentumPion = (type == 2 && binidx >= 956);
228 bool isMiddleMomentumGG4pi = false; //(type == 6 && ptrk > 0.75 && ptrk < 0.9);
229 bool isHighMomentumGG4pi = false; //(type == 6 && ptrk > 0.9 && ptrk < 1.2);
230 bool isHighMomentumElectron = (type == 4 && ptrk > 3.0/2);
231 if (isLowMomentumProton) {
232 for(int idx_theta =(int)(angleidx/10)*10; idx_theta<(int)(angleidx/10)*10+10; idx_theta++){
233 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge<(int)(chargeidx/2)*2+2; idx_charge++){
234 histidx = binidx*2*na + idx_theta*2 + idx_charge;
235 h[histidx]->Fill(exraw[j]);
236 }
237 }
238 }// end of low momentum proton
239 else if(isLowMomentumKaon){
240 if(angleidx<na-2){
241 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
242 histidx = binidx*2*na + angleidx*2 + idx_charge;
243 h[histidx]->Fill(exraw[j]);
244 }
245 }
246 else{
247 for(int idx_theta=(int)(angleidx/10)*10; idx_theta < (int)(angleidx/10)*10+10; idx_theta++){
248 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
249 histidx = binidx*2*na + idx_theta*2 + idx_charge;
250 h[histidx]->Fill(exraw[j]);
251 } // for
252 } // for
253 }// else
254 }// end of low mometnum kaon
255 else if(isHighMomentumKaon){
256 for(int idx_bg = (int)(binidx/2)*2; idx_bg < (int)(binidx/2)*2+2; idx_bg++){
257 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
258 histidx = idx_bg*2*na + angleidx*2 + idx_charge;
259 h[histidx]->Fill(exraw[j]);
260 } // for
261 } // for
262 } // end of low momentum pion
263 else if(isLowMomentumPion){
264 for(int idx_theta = (int)(angleidx/10)*10; idx_theta < (int)(angleidx/10)*10+10; idx_theta++){
265 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
266 histidx = binidx*2*na + idx_theta*2 + idx_charge;
267 h[histidx]->Fill(exraw[j]);
268 } // for
269 } // for
270 } // end of low momentum pion
271 else if(isHighMomentumPion){
272 for(int idx_bg = (int)(binidx/2)*2; idx_bg < (int)(binidx/2)*2+2; idx_bg++){
273 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
274 histidx = idx_bg*2*na + angleidx*2 + idx_charge;
275 h[histidx]->Fill(exraw[j]);
276 } // for
277 } // for
278 } // end of high momentum pion
279 else if(isMiddleMomentumGG4pi){
280 for(int idx_bg = (int)(binidx/3)*3; idx_bg < (int)(binidx/3)*3+3; idx_bg++){
281 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
282 histidx = idx_bg*2*na + angleidx*2 + idx_charge;
283 h[histidx]->Fill(exraw[j]);
284 } // for
285 } // for
286 } // end of middle momentum gg4pi
287 else if(isHighMomentumGG4pi){
288 for(int idx_bg = (int)(binidx/6)*6; idx_bg < (int)(binidx/6)*6+6; idx_bg++){
289 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
290 histidx = idx_bg*2*na + angleidx*2 + idx_charge;
291 h[histidx]->Fill(exraw[j]);
292 } // for
293 } // for
294 } // end of high momentum gg4pi
295 else if(isHighMomentumElectron){
296 for(int idx_theta = (int)(angleidx/10)*10; idx_theta < (int)(angleidx/10)*10+10; idx_theta++){
297 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
298 histidx = binidx*2*na + idx_theta*2 + idx_charge;
299 h[histidx]->Fill(exraw[j]);
300 } // for
301 } // for
302 } // end of high momentum electron
303 else{
304 for(int idx_charge = (int)(chargeidx/2)*2; idx_charge < (int)(chargeidx/2)*2+2; idx_charge++){
305 histidx = binidx*2*na + angleidx*2 + idx_charge;
306 h[histidx]->Fill(exraw[j]);
307 }
308 } // end of else
309 } // end of loop in momentum and angle
310 } //loop j, ndedxhit
311 } //loop i, entries in trk
312
313 //push the histograms and betagammas to the vector
314 for (int i = 0; i < nbin * na * 2; i++)
315 {
316 hist.Add(h[i]);
317 }
318 for (int i = 0; i < nbin; i++)
319 {
320 bg.push_back((pmin + pstep * i) / mass);
321 bg.push_back((pmin + pstep * (i + 1)) / mass);
322 }
323}
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
float costheta
float charge
float dEdx_meas
float bg

Referenced by main().

◆ histps()

void histps ( string  name)

Definition at line 325 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

325 {
326 TFile f(name.c_str());
327 TTree* t = (TTree*)f.Get("bin");
328 Int_t total, bgb;
329 t->SetBranchAddress("totalNum", &total);
330 t->SetBranchAddress("betagammaBounds", &bgb);
331 t->GetEntry(0);
332 TH1F* h;
333 stringstream ss;
334 TCanvas c1("c1", "canvas", 800, 600);
335 c1.Print((name + ".ps[").c_str());
336 for (Int_t i = 0; i < total; i++) {
337 if(i%100==0) cout << "i: " << i << " ps is written now" << endl;
338 ss.str("");
339 ss << "h" << i;
340 h = (TH1F*)f.Get(ss.str().c_str());
341 h->Draw();
342 c1.Print((name + ".ps").c_str());
343 }
344 c1.Print((name + ".ps]").c_str());
345}
TTree * t
Definition: binning.cxx:23
char * c_str(Index i)
Definition: EvtCyclic3.cc:252

Referenced by main().

◆ main()

int main ( )

Definition at line 62 of file Mdc/DedxCalibAlg/DedxCalibAlg-00-02-01/share/template/Simulation/histgen/test.cxx.

62 {
63 std::cout << "in main()" << std::endl;
64 //call fuction histgen to generate histograms
65 TObjArray histCol;
66 vector<double> bg;
68
69 std::cout << "before histgen" << std::endl;
70 histgen("../hadron_track/proton/proton.root", 0, histCol, bg);
71 std::cout << "proton hists generated" << std::endl;
72 histgen("../hadron_track/kaon/kaon.root", 1, histCol, bg);
73 std::cout << "kaon hists generated" << std::endl;
74 histgen("../hadron_track/pion/pion.root", 2, histCol, bg);
75 std::cout << "pion hists generated" << std::endl;
76 histgen("../hadron_track/muon/muon.root", 3, histCol, bg);
77 std::cout << "muon hists generated" << std::endl;
78 histgen("../hadron_track/electron/electron.root", 4, histCol, bg);
79 std::cout << "electron hists generated" << std::endl;
80
81 //define file and write to it
82 string filename = "pion_test.root";
83 TFile* f = new TFile(filename.c_str(), "recreate");
84
85 int total = histCol.GetLast() + 1;
86 int bg_num = bg.size();
87 cout << "total: " << total << " bg num = " << bg_num << endl;
88 double betagamma[3000];
89 for (unsigned int i = 0; i < bg.size(); i++) betagamma[i] = bg[i];
90
91 TTree* t = new TTree("bin", "bin information");
92 t->Branch("totalNum", &total, "total/I");
93 t->Branch("betagammaBounds", &bg_num, "bg_num/I");
94 t->Branch("betagamma", betagamma, "betagamma[bg_num]/D");
95 t->Fill();
96 t->Write();
97
98 std::cout << "begin write hists" << std::endl;
99 for (int i=0; i<total; i++){
100 if(i%100==0) cout << "i: " << i << " histogram is written now" << endl;
101 ((TH1F*)histCol[i])->Write();
102 }
103 std::cout << "write hists finish" << std::endl;
104
105 f->Close();
106 if(ps_flag) histps(filename);
107 hist_update();
108 return 0;
109}
void histgen(string, int, TObjArray &, vector< double > &)
const bool ps_flag(false)
void ReadCurPara(int index)

◆ ps_flag()

const bool ps_flag ( false  )

Referenced by main().

Variable Documentation

◆ electron_step

const double electron_step = (double)(P_electron_max - P_electron_min) / (double)P_electron_bin

◆ kaon_step

const double kaon_step = (double)(P_kaon_max - P_kaon_min) / (double)P_kaon_bin

◆ muon_step

const double muon_step = (double)(P_muon_max - P_muon_min) / (double)P_muon_bin

◆ na

const int na = 10

◆ nc

◆ pion_step

const double pion_step = (double)(P_pion_max - P_pion_min) / (double)P_pion_bin

◆ proton_step

const double proton_step = (double)(P_proton_max - P_proton_min) / (double)P_proton_bin

◆ run_flag

const int run_flag[5] = {0,0,1,0,0}