CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcPID.cxx
Go to the documentation of this file.
1#include <fstream>
2#include <cmath>
3#include <cstdlib>
4
5#include "TTree.h"
6#include "TMultiLayerPerceptron.h"
7
8#include "ParticleID/EmcPID.h"
9
10#ifndef BEAN
11#include "EvtRecEvent/EvtRecTrack.h"
12#include "MdcRecEvent/RecMdcTrack.h"
13#include "ExtEvent/RecExtTrack.h"
14#endif
15
16void CALG(double Px, double& e2); // see "calg.cxx"
17
18EmcPID * EmcPID::m_pointer = 0;
19TMultiLayerPerceptron * EmcPID::m_mlp_emc=0;
20TTree * EmcPID::m_trainTree_emconly=0;
21
22
23EmcPID * EmcPID::instance() {
24 if(!m_pointer) m_pointer = new EmcPID();
25 return m_pointer;
26}
27
28EmcPID::EmcPID():ParticleIDBase() {
29 m_mlp_emc = 0;
30 std::string e_emc_file = path + "/share/elec_emc_hist.txt";
31 //std::string e_emc_file = "$PARTICLEIDROOT/share/elechist.txt";
32 ifstream input(e_emc_file.c_str(),std::ios_base::in);
33 if ( !input ) {
34 cout << " can not open: " << e_emc_file << endl;
35 exit(1);
36 }
37 for(int i=0; i<18; i++) {
38 for(int j=0; j<300; j++) {
39 input>>m_e_h[i][j];
40 }
41 }
42 // std::string pi_emc_file = "/ihepbatch/bes/huangb/boss610/Analysis/ParticleID/ParticleID-00-02-03/share/pionhist.txt";
43 std::string pi_emc_file = path + "/share/pion_emc_hist.txt";
44 //std::string pi_emc_file = "$PARTICLEIDROOT/share/pionhist.txt";
45 ifstream input1(pi_emc_file.c_str(),std::ios_base::in);
46 if ( !input1 ) {
47 cout << " can not open: " << pi_emc_file << endl;
48 exit(1);
49 }
50 for(int i=0; i<18; i++) {
51 for(int j=0; j<300; j++) {
52 input1>>m_p_h[i][j];
53 }
54 }
55 std::string mu_emc_file = path + "/share/muon_emc_hist.txt";
56 // std::string mu_emc_file = "/ihepbatch/bes/huangb/boss610/Analysis/ParticleID/ParticleID-00-02-03/share/muonhist.txt";
57 //std::string mu_emc_file = "$PARTICLEIDROOT/share/muonhist.txt";
58 ifstream input2(mu_emc_file.c_str(),std::ios_base::in);
59 if ( !input2 ) {
60 cout << " can not open: " << mu_emc_file << endl;
61 exit(1);
62 }
63 for(int i=0; i<18; i++) {
64 for(int j=0; j<300; j++) {
65 input2>>m_m_h[i][j];
66
67 }
68 }
69
70 if(!m_trainTree_emconly) {
71 m_trainTree_emconly = new TTree("m_trainTree_emconly","m_trainTree_emconly");
72 m_trainTree_emconly->Branch("ptrk",&m_ptrk,"ptrk/D");
73 m_trainTree_emconly->Branch("pt",&m_pt,"pt/D");
74 m_trainTree_emconly->Branch("type",&m_type,"type/D");
75 m_trainTree_emconly->Branch("energy",&m_energy,"energy/D");
76 m_trainTree_emconly->Branch("eseed",&m_eseed,"eseed/D");
77 m_trainTree_emconly->Branch("e3x3",&m_e3x3,"e3x3/D");
78 m_trainTree_emconly->Branch("e5x5",&m_e5x5,"e5x5/D");
79 m_trainTree_emconly->Branch("latmoment",&m_latmoment,"latmoment/D");
80 m_trainTree_emconly->Branch("a20moment",&m_a20moment,"a20moment/D");
81 m_trainTree_emconly->Branch("a42moment",&m_a42moment,"a42moment/D");
82 m_trainTree_emconly->Branch("secondmoment",&m_secondmoment,"secondmoment/D");
83 m_trainTree_emconly->Branch("delta_phi",&m_delta_phi,"delta_phi/D");
84 m_trainTree_emconly->Branch("delta_theta",&m_delta_theta,"delta_theta/D");
85 }
86 std::string emc = path + "/share/emc.txt";
87 if(!m_mlp_emc) {
88 // m_mlp_emc = new TMultiLayerPerceptron("ptrk,pt,energy,eseed,e3x3,e5x5,secondmoment,delta_theta,delta_phi:18:type",m_trainTree_emconly);
89 m_mlp_emc = new TMultiLayerPerceptron("ptrk,pt,energy,eseed,e3x3,e5x5,secondmoment,latmoment,a20moment,a42moment,delta_theta,delta_phi:24:type",m_trainTree_emconly);
90
91 m_mlp_emc->LoadWeights(emc.c_str());
92
93 }
94}
95void EmcPID::init() {
96 for(int i = 0; i < 5; i++) {
97 m_chi[i] = 99.0;
98 m_prob[i] = -1.0;
99 }
100 m_chimin = 99.;
101 m_ndof = 0;
102 m_energy = -99;
103 m_eseed = -99;
104 m_e3x3 = -99;
105 m_e5x5 = -99;
106 m_delta_theta = -99;
107 m_delta_phi = -99;
108 m_secondmoment = -99;
109 m_val_emc = -99;
110 // std::string emc = path + "/share/emc_epimu.txt";
111 /* if(!m_trainTree_emconly){
112 m_trainTree_emconly = new TTree("m_trainTree_emconly","m_trainTree_emconly");
113 m_trainTree_emconly->Branch("ptrk",&m_ptrk,"ptrk/D");
114 m_trainTree_emconly->Branch("pt",&m_pt,"pt/D");
115 m_trainTree_emconly->Branch("type",&m_type,"type/D");
116 m_trainTree_emconly->Branch("energy",&m_energy,"energy/D");
117 m_trainTree_emconly->Branch("eseed",&m_eseed,"eseed/D");
118 m_trainTree_emconly->Branch("e3x3",&m_e3x3,"e3x3/D");
119 m_trainTree_emconly->Branch("e5x5",&m_e5x5,"e5x5/D");
120 m_trainTree_emconly->Branch("secondmoment",&m_secondmoment,"secondmoment/D");
121 m_trainTree_emconly->Branch("delta_phi",&m_delta_phi,"delta_phi/D");
122 m_trainTree_emconly->Branch("delta_theta",&m_delta_theta,"delta_theta/D");
123 }
124 if(!m_mlp_emc){
125 m_mlp_emc = new TMultiLayerPerceptron("ptrk,pt,energy,eseed,e3x3,e5x5,secondmoment,delta_theta,delta_phi:18:type",m_trainTree_emconly);
126 m_mlp_emc->LoadWeights(emc.c_str());
127 }
128 */
129
130 /*if(!m_trainTree_emconly){
131 m_trainTree_emconly = new TTree("m_trainTree_emconly","m_trainTree_emconly");
132 m_trainTree_emconly->Branch("ptrk",&m_ptrk,"ptrk/D");
133 m_trainTree_emconly->Branch("pt",&m_pt,"pt/D");
134 m_trainTree_emconly->Branch("type",&m_type,"type/D");
135 m_trainTree_emconly->Branch("energy",&m_energy,"energy/D");
136 m_trainTree_emconly->Branch("eseed",&m_eseed,"eseed/D");
137 m_trainTree_emconly->Branch("e3x3",&m_e3x3,"e3x3/D");
138 m_trainTree_emconly->Branch("e5x5",&m_e5x5,"e5x5/D");
139 m_trainTree_emconly->Branch("secondmoment",&m_secondmoment,"secondmoment/D");
140 m_trainTree_emconly->Branch("delta_phi",&m_delta_phi,"delta_phi/D");
141 m_trainTree_emconly->Branch("delta_theta",&m_delta_theta,"delta_theta/D");
142 }
143 if(!m_mlp_emc){
144 m_mlp_emc = new TMultiLayerPerceptron("ptrk,pt,energy,eseed,e3x3,e5x5,secondmoment,delta_theta,delta_phi:18:type",m_trainTree_emconly);
145 m_mlp_emc->LoadWeights(emc.c_str());
146 }
147 */
148
149
150
151}
152
153void EmcPID::calculate() {
154 if(particleIDCalculation() == 0) m_ndof = 1;
155}
156
157int EmcPID::particleIDCalculation() {
158 int irc = -1;
159 EvtRecTrack* recTrk = PidTrk();
160 if(!(recTrk->isMdcTrackValid())) return irc;
161 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
162 if (!(recTrk)->isMdcKalTrackValid()) return irc;
163 RecMdcKalTrack* mdcKalTrk = (recTrk)->mdcKalTrack();
164 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
165 double ptrk=mdcKalTrk->p();
166 // double ptrk = mdcTrk->p();
167 m_ptrk = ptrk;
168
169 m_pt=mdcKalTrk->pxy();
170 m_pt = m_pt*mdcTrk->charge();
171 // double cost = cos(mdcTrk->theta());
172
173 if(!(recTrk->isExtTrackValid())) return irc;
174 RecExtTrack* extTrk = recTrk->extTrack();
175 if(extTrk->emcVolumeNumber() == -1) return irc;
176 if (!(recTrk->isEmcShowerValid())) return irc;
177 RecEmcShower* emcTrk = recTrk->emcShower();
178
179 m_energy = emcTrk->energy();
180 m_eseed = emcTrk->eSeed();
181 m_e3x3 = emcTrk->e3x3();
182 m_e5x5 = emcTrk->e5x5();
183
184 double m_emc_theta = emcTrk->theta();
185 double m_emc_phi = emcTrk->phi();
186
187 Hep3Vector mc = extTrk->emcPosition();
188 double m_ext_theta = mc.theta();
189 double m_ext_phi = mc.phi();
190
191
192 m_delta_theta = m_emc_theta - m_ext_theta;
193 m_delta_phi = m_emc_phi - m_ext_phi;
194 if(m_delta_phi>1) m_delta_phi=m_delta_phi-6.283;
195 if(m_delta_phi<-1) m_delta_phi=m_delta_phi+6.283;
196
197
198 m_secondmoment = emcTrk->secondMoment()/1000.;
199 m_a20moment=emcTrk->a20Moment();
200 m_latmoment=emcTrk->latMoment();
201 m_a42moment=emcTrk->a42Moment();
202
203 if(emcTrk->energy() <= 0) return irc;
204 //if(emcTrk->energy() > ptrk) return irc;
205
206 /* params_emc1[0] = m_ptrk;
207 params_emc1[1] = m_pt;
208 params_emc1[2] = m_energy;
209 params_emc1[3] = m_eseed;
210 params_emc1[4] = m_e3x3;
211 params_emc1[5] = m_e5x5;
212 params_emc1[6] = m_secondmoment;
213 params_emc1[7] = m_delta_theta;
214 params_emc1[8] = m_delta_phi;*/
215 params_emc1[0] =m_ptrk;
216 params_emc1[1] =m_pt;
217 params_emc1[2] =m_energy;
218 params_emc1[3] =m_eseed;
219 params_emc1[4] =m_e3x3;
220 params_emc1[5] =m_e5x5;
221 params_emc1[6] =m_secondmoment;
222 params_emc1[7] =m_latmoment;
223 params_emc1[8] =m_a20moment;
224 params_emc1[9] =m_a42moment;
225 params_emc1[10] =m_delta_theta;
226 params_emc1[11] =m_delta_phi;
227
228 m_val_emc = m_mlp_emc->Evaluate(0,params_emc1);
229 int pindex = int((m_ptrk-0.2)/0.1);
230 int bindex = int((m_val_emc-0.5)/0.01);
231 if(bindex>300||bindex<0) return irc;
232 if(pindex>17) pindex=17;
233 if(pindex<0) pindex=0;
234 // double bin_pos[3];
235 m_prob[0] = m_e_h[pindex][bindex];
236 m_prob[1] = m_m_h[pindex][bindex];
237 m_prob[2] = m_p_h[pindex][bindex];
238 m_prob[3] = m_p_h[pindex][bindex];
239 m_prob[4] = m_p_h[pindex][bindex];
240 for(int i =0; i<5; i++) {
241 if(m_prob[i]==0) m_prob[i] = 0.001;
242 }
243 //calculate the chisq value using GAUSIN
244 float ppp[5];
245 for(int i=0; i<5; i++) {
246 ppp[i]=0;
247 }
248 for(int j=0; j<=bindex; j++) {
249 ppp[0]+= m_e_h[pindex][j];
250 ppp[1]+= m_m_h[pindex][j];
251 ppp[2]+= m_p_h[pindex][j];
252 }
253 for(int i=0; i<3; i++) {
254 ppp[i]=ppp[i]*0.01;
255 if(ppp[i]>0&&ppp[i]<1) {
256 CALG(ppp[i],m_chi[i]);
257 }
258 if(ppp[i]<=0||ppp[i]>=1) m_chi[i]=-99;
259 }
260 // if(fabs(m_chi[2])==-99)
261 m_chi[3]=m_chi[2];
262 m_chi[4]=m_chi[2];
263
264 m_ndof = 1;
265 irc = 0;
266 return irc;
267}
Px
double Px(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:58
e2
Double_t e2
Definition: DataBase/tau_mode.c:7
CALG
void CALG(double Px, double &e2)
Definition: calg.cxx:31
CALG
void CALG(double Px, double &e2)
Definition: calg.cxx:31
EmcPID::calculate
void calculate()
Definition: EmcPID.cxx:153
EmcPID::init
void init()
Definition: EmcPID.cxx:95
EmcPID::particleIDCalculation
int particleIDCalculation()
Definition: EmcPID.cxx:157
EmcPID::instance
static EmcPID * instance()
Definition: EmcPID.cxx:23