BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
WidgetGenerator Class Reference

#include <WidgetGenerator.h>

Public Member Functions

 WidgetGenerator ()
 
 WidgetGenerator (int nevents, double upperbg, double lowerbg)
 
virtual ~WidgetGenerator ()
 
void generateEvents (TString filename, double genmass)
 
void simulateDedx (TString infilename, TString filename, double genmass)
 
void simulateReconstruction (TString infilename, TString outfilename)
 

Detailed Description

Definition at line 34 of file WidgetGenerator.h.

Constructor & Destructor Documentation

◆ WidgetGenerator() [1/2]

WidgetGenerator::WidgetGenerator ( )

Definition at line 3 of file WidgetGenerator.cc.

3 :
4 m_nevents(1000000),
5 m_upperbg(5.0),
6 m_lowerbg(0.1)
7{}

◆ WidgetGenerator() [2/2]

WidgetGenerator::WidgetGenerator ( int  nevents,
double  upperbg,
double  lowerbg 
)

Definition at line 9 of file WidgetGenerator.cc.

9 :
10 m_nevents(nevents),
11 m_upperbg(upperbg),
12 m_lowerbg(lowerbg)
13{}

◆ ~WidgetGenerator()

virtual WidgetGenerator::~WidgetGenerator ( )
inlinevirtual

Definition at line 40 of file WidgetGenerator.h.

40{};

Member Function Documentation

◆ generateEvents()

void WidgetGenerator::generateEvents ( TString  filename,
double  genmass 
)

Definition at line 16 of file WidgetGenerator.cc.

16 {
17
18 std::cout << "WidgetGenerator: generating " << m_nevents << " particles of mass " << genmass << std::endl;
19 std::cout << m_upperbg << "\t" << m_lowerbg << std::endl;
20
21
22 double dedxpub; // dE/dx without electron saturation correction
23 double dedx; // dE/dx with electron saturation correction
24 double p; // track momentum
25 double bg; // track beta-gamma
26 double costh; // cosine of track polar angle
27 double db; // the nearest distance to the IP in the xy plane
28 double dz; // the nearest distance to the IP in the z plane
29 double chiPi; // PID chi value
30 double nhit; // number of hits on this track
31 double eop; // energy over momentum in the calorimeter
32
33 double dedx_cur, dedx_res;
34
35 // do not overwrite a file if it already exists
36 TFile* genfile = new TFile(filename,"CREATE");
37 TTree* gentree = new TTree("track","Fake track sample");
38
39 gentree->Branch("dedx",&dedxpub,"dedx/D");
40 gentree->Branch("dedxcur",&dedx_cur,"dedxcur/D");
41 gentree->Branch("dedxres",&dedx_res,"dedxres/D");
42 gentree->Branch("dedxsat",&dedx,"dedxsat/D");
43 gentree->Branch("pF",&p,"pF/D");
44 gentree->Branch("bg",&bg,"bg/D");
45 gentree->Branch("costh",&costh,"costh/D");
46 gentree->Branch("db",&db,"db/D");
47 gentree->Branch("dz",&dz,"dz/D");
48 gentree->Branch("chiPi",&chiPi,"chiPi/D");
49 gentree->Branch("numGoodHits",&nhit,"numGoodHits/D");
50 gentree->Branch("eopst",&eop,"eopst/D");
51
52 TH2F* dedx_costh = new TH2F("dedx_costh","dE/dx vs. cos(#theta);cos(#theta);dE/dx",
53 100,-1.0,1.0,100,0.0,10.0);
54 TH2F* dedx_bg = new TH2F("dedx_bg","dE/dx vs. #beta#gamma;#beta#gamma;dE/dx",
55 100,m_lowerbg,m_upperbg,100,0.0,10.0);
56
57 WidgetParameterization gpar;
58 HadronSaturation hadsat;
59
60 TRandom *r0 = new TRandom();
61 // set the seed to the current machine clock
62 r0->SetSeed(0);
63
64 if( genmass == 0 ){
65 std::cout << "WidgetGenerator::ERROR - please input non-zero particle mass" << std::endl;
66 return;
67 }
68
69 double mass = genmass;
70 double massvec[5];
71 massvec[0] = Widget::mpion;
72 massvec[1] = Widget::mkaon;
73 massvec[2] = Widget::mproton;
74 massvec[3] = Widget::mmuon;
75 massvec[4] = Widget::melectron;
76
77 // generate the events
78 for( int i = 0; i < m_nevents; ++i ){
79
80 bg = r0->Rndm(i)*(m_upperbg+m_lowerbg)+m_lowerbg;
81 while( bg < m_lowerbg || bg > m_upperbg )
82 bg = r0->Rndm(i)*(m_upperbg+m_lowerbg)+m_lowerbg;
83 p = bg*mass;
84
85 costh = 2*(r0->Rndm(i))-1;
86 while( costh < -1.0 || costh > 1.0 )
87 costh = 2*(r0->Rndm(i))-1;
88
89 nhit = r0->Gaus(25);
90 eop = 0.5;
91 chiPi = r0->Gaus(1.0);
92
93 dedx_cur = gpar.dedxPrediction(bg);
94 dedx_res = gpar.sigmaPrediction(bg,nhit,sqrt(1-costh*costh));
95 if( mass == Widget::melectron || mass == Widget::mmuon ) dedx_res = 0.05;
96
97 // fake the hadron saturation
98 dedx = r0->Gaus(hadsat.I2D(costh,dedx_cur)/hadsat.I2D(costh,1.0),dedx_res/3.0);
99 // if( dedx > 2 ) continue;
100
101 /* while( dedx > 2 )
102 dedx = r0->Gaus(hadsat.I2D(costh,dedx_cur)/hadsat.I2D(costh,1.0),dedx_res);*/
103
104 dedxpub = dedx;
105
106 dedx_costh->Fill(costh,dedx);
107 dedx_bg->Fill(bg,dedx);
108
109 gentree->Fill();
110 }
111 delete r0;
112
113 genfile->cd();
114 gentree->Write();
115 dedx_costh->Write();
116 dedx_bg->Write();
117 genfile->Close();
118}
mass
double mass
Definition: CosmicGenerator.cxx:138
HadronSaturation::I2D
double I2D(double cosTheta, double I, double alpha, double gamma, double delta, double power, double ratio) const
Definition: HadronSaturation.cc:133
WidgetParameterization::sigmaPrediction
double sigmaPrediction(double dedx, double nhit, double sin)
Definition: WidgetParameterization.cc:130
bg
float bg
Definition: prepare_data.cxx:28

Referenced by HadronInterface::GenerateSample().

◆ simulateDedx()

void WidgetGenerator::simulateDedx ( TString  infilename,
TString  filename,
double  genmass 
)

Definition at line 121 of file WidgetGenerator.cc.

121 {
122
123 TFile* infile = new TFile(infilename);
124 TTree* intree = (TTree*)infile->Get("track");
125
126 double dedx, bg, costh, p;
127 intree->SetBranchAddress("p",&p);
128 intree->SetBranchAddress("costh",&costh);
129
130 // do not overwrite a file that already exists
131 TFile* outfile = new TFile(outfilename,"CREATE");
132 TTree* tracks = new TTree("track","Fake track sample");
133 tracks->Branch("dedxsat",&dedx,"dedxsat/D");
134 tracks->Branch("p",&p,"p/D");
135 tracks->Branch("costh",&costh,"costh/D");
136
137 TH2F* dedx_costh = new TH2F("dedx_costh","dE/dx vs. cos(#theta);cos(#theta);dE/dx",
138 100,-1.0,1.0,100,0.0,10.0);
139 TH2F* dedx_bg = new TH2F("dedx_bg","dE/dx vs. #beta#gamma;#beta#gamma;dE/dx",
140 100,m_lowerbg,m_upperbg,100,0.0,10.0);
141
142 HadronSaturation hadsat;
143
144 TRandom *r0 = new TRandom();
145 // set the seed to the current machine clock
146 r0->SetSeed(0);
147
148 double mass = genmass;
149 double massvec[5];
150 massvec[0] = Widget::mpion;
151 massvec[1] = Widget::mkaon;
152 massvec[2] = Widget::mproton;
153 massvec[3] = Widget::mmuon;
154 massvec[4] = Widget::melectron;
155
156 // use the real momentum and angle, but replace the dE/dx measurement
157 for( int i = 0; i < intree->GetEntries(); ++i ){
158 intree->GetEvent(i);
159
160 if( i < 10 )
161 std::cout << "Event " << i << ": p = " << p << ", cos(theta) = " << costh << std::endl;
162
163 if( genmass == 0 ) mass = massvec[i%5];
164 bg = p/mass;
165
166 // fake the hadron saturation
167 dedx = r0->Gaus(hadsat.I2D(costh,0.8)/hadsat.I2D(costh,1.0),0.05);
168 while( dedx > 2 )
169 dedx = r0->Gaus(hadsat.I2D(costh,0.8)/hadsat.I2D(costh,1.0),0.05);
170
171 dedx_costh->Fill(costh,dedx);
172 dedx_bg->Fill(bg,dedx);
173
174 tracks->Fill();
175 }
176 delete r0;
177
178 outfile->cd();
179 tracks->Write();
180 dedx_costh->Write();
181 dedx_bg->Write();
182 outfile->Close();
183}

Referenced by HadronInterface::SimulateSample().

◆ simulateReconstruction()

void WidgetGenerator::simulateReconstruction ( TString  infilename,
TString  outfilename 
)

Definition at line 186 of file WidgetGenerator.cc.

186 {
187
188 TFile* infile = new TFile(infilename);
189 TTree* intree = (TTree*)infile->Get("track");
190
191 int nhits;
192 double dedxhit[100], adcraw[100], path[100];
193 intree->SetBranchAddress("dedxhit",dedxhit);
194 intree->SetBranchAddress("adcraw",adcraw);
195 intree->SetBranchAddress("path",path);
196
197 double layer[100], layerdedx[100], layerdx[100];
198 intree->SetBranchAddress("layer",layer);
199 intree->SetBranchAddress("layerdx",layerdx);
200 intree->SetBranchAddress("layerdedx",layerdedx);
201
202 double mean, dedx, dedxerr, costh;
203 intree->SetBranchAddress("numLayerHits",&nhits);
204 intree->SetBranchAddress("mean",&mean);
205 intree->SetBranchAddress("dedx",&dedx);
206 intree->SetBranchAddress("dedxerr",&dedxerr);
207 intree->SetBranchAddress("costh",&costh);
208
209 // do not overwrite a file that already exists
210 TFile* outfile = new TFile(outfilename,"CREATE");
211 TTree* tracks = intree->CloneTree(0);
212
213 double newdedxhit[100];
214 double newmean, newdedx, newerr;
215 tracks->Branch("newmean",&newmean,"newmean/D");
216 tracks->Branch("newdedx",&newdedx,"newdedx/D");
217 tracks->Branch("newerr",&newerr,"newerr/D");
218
219 for( int i = 0; i < intree->GetEntries(); ++i ){
220 intree->GetEntry(i);
221
222 // clear the simulated measurements
223 newmean = 0; newdedx = 0; newerr = 0;
224
225 // sort the dedx hits for truncation
226 std::vector<double> sortedDedx;
227 for( int j = 0; j < nhits; ++j ){
228 sortedDedx.push_back(layerdedx[j]);
229 }
230 std::sort(sortedDedx.begin(), sortedDedx.end());
231
232 // recalculate the mean values
233 int lb = int(nhits * 0.05 + 0.5);
234 int ub = int(nhits * 0.75 + 0.5);
235 double ntrunc = 0;
236 for( int j = 0; j < nhits; ++j ){
237 newmean += sortedDedx[j];
238 if( j >= lb && j < ub ){
239 newdedx += sortedDedx[j];
240 newerr += sortedDedx[j]*sortedDedx[j];
241 ntrunc++;
242 }
243 }
244
245 if( nhits != 0 ) newmean /= nhits;
246 if( ntrunc != 0 ) newdedx /= ntrunc;
247 else newdedx = newmean;
248
249 if( ntrunc > 1 ) newerr = sqrt(newerr/ntrunc - newdedx*newdedx)/(ntrunc-1);
250 else newerr = 0;
251
252 tracks->Fill();
253 }
254
255 tracks->AutoSave();
256 outfile->Close();
257 infile->Close();
258}
nhits
int nhits
Definition: prepare_data.cxx:29

Referenced by main(), and HadronInterface::SimulateReconstruction().

The documentation for this class was generated from the following files: