00001 #include <fstream>
00002 #include <cmath>
00003 #include <cstdlib>
00004
00005 #include "TFile.h"
00006 #include "TTree.h"
00007 #include "TMultiLayerPerceptron.h"
00008 #include "TMLPAnalyzer.h"
00009
00010 #include "ParticleID/MucPID.h"
00011
00012 #ifndef BEAN
00013 #include "EvtRecEvent/EvtRecTrack.h"
00014 #include "MdcRecEvent/RecMdcTrack.h"
00015 #include "MucRecEvent/RecMucTrack.h"
00016 #endif
00017
00018 void CALG(double Px, double& e2);
00019
00020 MucPID * MucPID::m_pointer = 0;
00021
00022 MucPID * MucPID::instance() {
00023
00024 if(!m_pointer) m_pointer = new MucPID();
00025 return m_pointer;
00026 }
00027
00028 MucPID::MucPID():ParticleIDBase() {
00029 m_trainFile_muc = 0;
00030 m_trainTree_muc = 0;
00031 m_mlp_muc = 0;
00032 m_mlpa_muc = 0;
00033
00034
00035 std::string pi_muc_file = path + "/share/pion_muc_hist.txt";
00036
00037 ifstream input1(pi_muc_file.c_str(),std::ios_base::in);
00038 if ( !input1 ) {
00039 cout << " can not open: " << pi_muc_file << endl;
00040 exit(1);
00041 }
00042 for(int i=0; i<13; i++) {
00043 for(int j=0; j<300; j++) {
00044 input1>>m_p_h[i][j];
00045 }
00046 }
00047
00048 std::string mu_muc_file = path + "/share/muon_muc_hist.txt";
00049
00050
00051 ifstream input2(mu_muc_file.c_str(),std::ios_base::in);
00052 if ( !input2 ) {
00053 cout << " can not open: " << mu_muc_file << endl;
00054 exit(1);
00055 }
00056 for(int i=0; i<13; i++) {
00057 for(int j=0; j<300; j++) {
00058 input2>>m_m_h[i][j];
00059 }
00060 }
00061
00062
00063
00064
00065 }
00066
00067 void MucPID::init() {
00068 for(int i = 0; i < 5; i++) {
00069 m_chi[i] = 99.0;
00070 m_prob[i] = -1.0;
00071 }
00072 m_chimin = 99.;
00073 m_ndof = 0;
00074 m_hits = -99;
00075 m_depth = -999;
00076 m_val_muc=-99;
00077
00078 std::string neural = path + "/share/neural.root";
00079 std::string muc = path + "/share/muc.txt";
00080 double ptrk,phi,theta,depth,hits,chi2,distance,muc_delta_phi;
00081 int type;
00082 if(!m_trainTree_muc) {
00083 m_trainTree_muc = new TTree("m_trainTree_muc","m_trainTree_muc");
00084 m_trainTree_muc->Branch("ptrk",&ptrk,"ptrk/D");
00085 m_trainTree_muc->Branch("phi",&phi,"phi/D");
00086 m_trainTree_muc->Branch("theta",&theta,"theta/D");
00087 m_trainTree_muc->Branch("depth",&depth,"depth/D");
00088 m_trainTree_muc->Branch("hits",&hits,"hits/D");
00089 m_trainTree_muc->Branch("chi2",&chi2,"chi2/D");
00090 m_trainTree_muc->Branch("distance",&distance,"distance/D");
00091 m_trainTree_muc->Branch("muc_delta_phi",&muc_delta_phi,"muc_delta_phi/D");
00092 m_trainTree_muc->Branch("type",&type,"type/I");
00093 }
00094 if(!m_mlp_muc) {
00095 m_mlp_muc = new TMultiLayerPerceptron("@ptrk,@phi,@theta,@depth,@hits,@chi2,@distance,@muc_delta_phi:16:type",m_trainTree_muc);
00096 m_mlp_muc->LoadWeights(muc.c_str());
00097 }
00098 }
00099
00100
00101 void MucPID::calculate() {
00102 if(particleIDCalculation() == 0) m_ndof = 1;
00103 }
00104
00105 int MucPID::particleIDCalculation() {
00106 int irc = -1;
00107 EvtRecTrack* recTrk = PidTrk();
00108 if(!(recTrk->isMdcTrackValid())) return irc;
00109 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
00110
00111 m_depth = -99;
00112 m_hits = -99;
00113 m_chi2 =-99;
00114 m_distance =-99;
00115 m_muc_delta_phi =-99;
00116
00117 double ptrk = mdcTrk->p();
00118 double m_ptrk = ptrk;
00119 double m_pt = mdcTrk->pxy();
00120 double phi = mdcTrk->phi();
00121 double theta = mdcTrk->theta();
00122 double cost = cos(mdcTrk->theta());
00123 if(ptrk<0.5) return irc;
00124 if(fabs(cost)>0.83) return irc;
00125 if(!(recTrk->isExtTrackValid())) return irc;
00126 RecExtTrack* extTrk = recTrk->extTrack();
00127 if(extTrk->mucVolumeNumber() == -1) return irc;
00128 if (!(recTrk->isMucTrackValid())) return irc;
00129 RecMucTrack* mucTrk = recTrk->mucTrack();
00130
00131
00132 if(mucTrk->depth()>100000) return irc;
00133
00134 m_hits = mucTrk->maxHitsInLayer();
00135 m_depth = mucTrk->depth();
00136 m_distance = mucTrk->distance();
00137 m_chi2 = mucTrk->chi2();
00138
00139
00140
00141
00142
00143
00144 m_muc_delta_phi= mucTrk->deltaPhi();
00145 m_muc_delta_phi=3.14159-m_muc_delta_phi;
00146 theta = cos(theta);
00147 params_muc1[0] = m_ptrk;
00148 params_muc1[1] = phi;
00149 params_muc1[2] = theta;
00150 params_muc1[3] = m_depth;
00151 params_muc1[4] = m_hits;
00152 params_muc1[5] = m_chi2;
00153 params_muc1[6] = m_distance;
00154 params_muc1[7] = m_muc_delta_phi;
00155
00156 m_val_muc = m_mlp_muc->Evaluate(0,params_muc1);
00157 if(m_pt<0) m_pt = -m_pt;
00158 int pindex = int((m_ptrk-0.5)/0.1);
00159 int bindex = int((m_val_muc-0.5)/0.01);
00160 if(bindex>300||bindex<0) return irc;
00161 if(pindex>11) pindex=11;
00162 if(pindex<0) pindex=0;
00163 m_prob[0] = m_p_h[pindex][bindex];;
00164 m_prob[1] = m_m_h[pindex][bindex];
00165 m_prob[2] = m_p_h[pindex][bindex];
00166 m_prob[3] = m_p_h[pindex][bindex];
00167 m_prob[4] = m_p_h[pindex][bindex];
00168 for(int i =0; i<5; i++) {
00169 if(m_prob[i]==0) m_prob[i] = 0.001;
00170 }
00171 float ppp[5];
00172 for(int i =0; i<5; i++) {
00173 ppp[i] = 0.0;
00174 }
00175
00176 for(int j=0; j<bindex; j++) {
00177 ppp[1]+= m_m_h[pindex][j]*0.01;
00178 ppp[2]+= m_p_h[pindex][j]*0.01;
00179 }
00180 if(ppp[1]>0&&ppp[1]<1)
00181 CALG(ppp[1],m_chi[1]);
00182 if(ppp[2]>0&&ppp[2]<1)
00183 CALG(ppp[2],m_chi[2]);
00184 if(ppp[1]<=0||ppp[1]>=1)
00185 m_chi[1]=-99;
00186 if(ppp[2]<=0||ppp[2]>=1)
00187 m_chi[2]=-99;
00188
00189 m_chi[3]=m_chi[2];
00190 m_chi[4]=m_chi[2];
00191 m_chi[0] =m_chi[2];
00192 m_ndof = 1;
00193 irc = 0;
00194 return irc;
00195 }