BOSS 6.6.4.p03
BESIII Offline Software System
Loading...
Searching...
No Matches
Gam4pikp.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/AlgFactory.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/PropertyMgr.h"
7
8#include "EventModel/EventModel.h"
9#include "EventModel/Event.h"
10
11#include "EvtRecEvent/EvtRecEvent.h"
12#include "EvtRecEvent/EvtRecTrack.h"
13#include "DstEvent/TofHitStatus.h"
14#include "EventModel/EventHeader.h"
15#include "McTruth/McParticle.h"
16#include "HltEvent/HltInf.h"
17#include "HltEvent/DstHltInf.h"
18
19#include "TMath.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "GaudiKernel/NTuple.h"
22#include "GaudiKernel/Bootstrap.h"
23#include "GaudiKernel/IHistogramSvc.h"
24#include "CLHEP/Vector/ThreeVector.h"
25#include "CLHEP/Vector/LorentzVector.h"
26#include "CLHEP/Vector/TwoVector.h"
27#include "VertexFit/IVertexDbSvc.h"
28#include "GaudiKernel/Bootstrap.h"
29#include "GaudiKernel/ISvcLocator.h"
30
31
32using CLHEP::Hep3Vector;
33using CLHEP::Hep2Vector;
34using CLHEP::HepLorentzVector;
35#include "CLHEP/Geometry/Point3D.h"
36#ifndef ENABLE_BACKWARDS_COMPATIBILITY
37typedef HepGeom::Point3D<double> HepPoint3D;
38#endif
39#include "Gam4pikpAlg/Gam4pikp.h"
40
41#include "VertexFit/KalmanKinematicFit.h"
42#include "VertexFit/VertexFit.h"
43#include "VertexFit/SecondVertexFit.h"
44#include "VertexFit/Helix.h"
45#include "ParticleID/ParticleID.h"
46#include <vector>
47const double mpi = 0.13957;
48const double mk = 0.493677;
49const double mpro = 0.938272;
50const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
51const double velc = 299.792458; // tof path unit in mm
52typedef std::vector<int> Vint;
53typedef std::vector<HepLorentzVector> Vp4;
54typedef std::vector<double> Vdouble;
55typedef std::vector<WTrackParameter> Vw;
56typedef std::vector<VertexParameter> Vv;
57static int Ncut[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
58static int Mcut[10]={0,0,0,0,0,0,0,0,0,0};
59/////////////////////////////////////////////////////////////////////////////
60
61Gam4pikp::Gam4pikp(const std::string& name, ISvcLocator* pSvcLocator) :
62 Algorithm(name, pSvcLocator) {
63
64 //Declare the properties
65 declareProperty("skim4pi", m_skim4pi=false);
66 declareProperty("skim4k", m_skim4k=false);
67 declareProperty("skim2pi2pr", m_skim2pi2pr=false);
68 declareProperty("rootput", m_rootput=false);
69 declareProperty("Ecms", m_ecms=3.6862);
70 declareProperty("Vr0cut", m_vr0cut=1.0);
71 declareProperty("Vz0cut", m_vz0cut=5.0);
72 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
73 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
74 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
75 declareProperty("GammaDangCut", m_gammadangcut=20.0);
76// declareProperty("Test4C", m_test4C = 1);
77// declareProperty("Test5C", m_test5C = 1);
78// declareProperty("CheckDedx", m_checkDedx = 1);
79// declareProperty("CheckTof", m_checkTof = 1);
80 }
81
82// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
83StatusCode Gam4pikp::initialize(){
84 MsgStream log(msgSvc(), name());
85
86 log << MSG::INFO << "in initialize()" << endmsg;
87// setFilterPassed(false);
88 StatusCode status;
89
90if(m_rootput)
91{
92 NTuplePtr nt1(ntupleSvc(), "FILE1/total4c");
93 if ( nt1 ) m_tuple1 = nt1;
94 else {
95 m_tuple1 = ntupleSvc()->book ("FILE1/total4c", CLID_ColumnWiseTuple, "ks N-Tuple example");
96 if ( m_tuple1 ) {
97
98 status = m_tuple1->addItem ("run", m_run );
99 status = m_tuple1->addItem ("rec", m_rec );
100 status = m_tuple1->addItem ("mpprecall", m_mpprecall );
101 status = m_tuple1->addItem ("meeall", m_meeall );
102 status = m_tuple1->addItem ("ncgjs", m_ncgjs );
103 status = m_tuple1->addItem ("cla2kpi", m_cla2kpi );
104 status = m_tuple1->addItem("indexmc", m_idxmc, 0, 100);
105 status = m_tuple1->addIndexedItem("pdgid", m_idxmc, m_pdgid);
106
107 status = m_tuple1->addIndexedItem("motheridx", m_idxmc, m_motheridx);
108 status = m_tuple1->addItem("indexmdc", m_idxmdc, 0, 5000);
109 status = m_tuple1->addIndexedItem ("x0js", m_idxmdc, m_x0js);
110 status = m_tuple1->addIndexedItem ("y0js", m_idxmdc, m_y0js);
111 status = m_tuple1->addIndexedItem ("z0js",m_idxmdc, m_z0js);
112 status = m_tuple1->addIndexedItem ("r0js",m_idxmdc, m_r0js);
113 status = m_tuple1->addIndexedItem ("Rxyjs",m_idxmdc, m_Rxyjs);
114 status = m_tuple1->addIndexedItem ("Rzjs",m_idxmdc, m_Rzjs);
115 status = m_tuple1->addIndexedItem ("Rnxyjs",m_idxmdc, m_Rnxyjs);
116 status = m_tuple1->addIndexedItem ("phinjs",m_idxmdc, m_phinjs);
117 status = m_tuple1->addIndexedItem ("Rnzjs",m_idxmdc, m_Rnzjs);
118 status = m_tuple1->addItem ("ncy20", m_ncy20);
119 status = m_tuple1->addItem ("ncy30", m_ncy30);
120 status = m_tuple1->addIndexedItem("angjs5", m_idxmdc, m_angjs5);
121 status = m_tuple1->addIndexedItem("nearjs5", m_idxmdc, m_nearjs5);
122 status = m_tuple1->addIndexedItem("angjs6", m_idxmdc, m_angjs6);
123 status = m_tuple1->addIndexedItem("nearjs6", m_idxmdc, m_nearjs6);
124 status = m_tuple1->addIndexedItem("ang4pi5", m_idxmdc, m_ang4pi5);
125 status = m_tuple1->addIndexedItem("near4pi5", m_idxmdc, m_near4pi5);
126 status = m_tuple1->addIndexedItem("ang4pi6", m_idxmdc, m_ang4pi6);
127 status = m_tuple1->addIndexedItem("near4pi6", m_idxmdc, m_near4pi6);
128 status = m_tuple1->addIndexedItem("ppmdcjs", m_idxmdc, m_ppmdcjs);
129 status = m_tuple1->addIndexedItem("pxmdcjs", m_idxmdc, m_pxmdcjs);
130 status = m_tuple1->addIndexedItem("pymdcjs", m_idxmdc, m_pymdcjs);
131 status = m_tuple1->addIndexedItem("pzmdcjs", m_idxmdc, m_pzmdcjs);
132 status = m_tuple1->addIndexedItem("ppkaljs", m_idxmdc, m_ppkaljs);
133 status = m_tuple1->addIndexedItem("ptmdcjs", m_idxmdc, m_ptmdcjs);
134 status = m_tuple1->addIndexedItem("ptkaljs", m_idxmdc, m_ptkaljs);
135 status = m_tuple1->addIndexedItem("ppmdc2kpi", m_idxmdc, m_ppmdc2kpi);
136 status = m_tuple1->addIndexedItem("pxmdc2kpi", m_idxmdc, m_pxmdc2kpi);
137 status = m_tuple1->addIndexedItem("pymdc2kpi", m_idxmdc, m_pymdc2kpi);
138 status = m_tuple1->addIndexedItem("pzmdc2kpi", m_idxmdc, m_pzmdc2kpi);
139 status = m_tuple1->addIndexedItem("ppkal2kpi", m_idxmdc, m_ppkal2kpi);
140 status = m_tuple1->addIndexedItem("ptmdc2kpi", m_idxmdc, m_ptmdc2kpi);
141 status = m_tuple1->addIndexedItem("charge2kpi", m_idxmdc, m_charge2kpi);
142 status = m_tuple1->addIndexedItem("ptkal2kpi", m_idxmdc, m_ptkal2kpi);
143 status = m_tuple1->addItem ("cy2pi", m_cy2kpi, 0, 100 );
144 status = m_tuple1->addIndexedItem("comcs2kpi", m_cy2kpi, m_comcs2kpi);
145 status = m_tuple1->addItem ("chiejs", m_idxmdc, m_chiejs);
146 status = m_tuple1->addItem ("chimujs", m_idxmdc, m_chimujs);
147 status = m_tuple1->addItem ("chipijs", m_idxmdc, m_chipijs);
148 status = m_tuple1->addItem ("chikjs", m_idxmdc, m_chikjs);
149 status = m_tuple1->addItem ("chipjs", m_idxmdc, m_chipjs);
150 status = m_tuple1->addItem ("ghitjs", m_idxmdc, m_ghitjs);
151 status = m_tuple1->addItem ("thitjs", m_idxmdc, m_thitjs);
152 status = m_tuple1->addIndexedItem("probphjs", m_idxmdc, m_probphjs);
153 status = m_tuple1->addIndexedItem("normphjs", m_idxmdc, m_normphjs);
154 status = m_tuple1->addItem ("pdg", m_idxmdc, m_pdg);
155 status = m_tuple1->addItem ("cbmc", m_idxmdc, m_cbmc);
156 status = m_tuple1->addIndexedItem("sigmaetof2kpi", m_idxmdc, m_sigmaetof2kpi);
157 status = m_tuple1->addIndexedItem("sigmamutof2kpi", m_idxmdc, m_sigmamutof2kpi);
158 status = m_tuple1->addIndexedItem("sigmapitof2kpi", m_idxmdc, m_sigmapitof2kpi);
159 status = m_tuple1->addIndexedItem("sigmaktof2kpi", m_idxmdc, m_sigmaktof2kpi);
160 status = m_tuple1->addIndexedItem("sigmaprtof2kpi", m_idxmdc, m_sigmaprtof2kpi);
161 status = m_tuple1->addIndexedItem("t0tof2kpi", m_idxmdc, m_t0tof2kpi);
162 status = m_tuple1->addIndexedItem("errt0tof2kpi", m_idxmdc, m_errt0tof2kpi);
163
164 status = m_tuple1->addItem ("chie2kpi", m_idxmdc, m_chie2kpi);
165 status = m_tuple1->addItem ("chimu2kpi", m_idxmdc, m_chimu2kpi);
166 status = m_tuple1->addItem ("chipi2kpi", m_idxmdc, m_chipi2kpi);
167 status = m_tuple1->addItem ("chik2kpi", m_idxmdc, m_chik2kpi);
168 status = m_tuple1->addItem ("chip2kpi", m_idxmdc, m_chip2kpi);
169 status = m_tuple1->addItem ("ghit2kpi", m_idxmdc, m_ghit2kpi);
170 status = m_tuple1->addItem ("thit2kpi", m_idxmdc, m_thit2kpi);
171 status = m_tuple1->addIndexedItem("probph2kpi", m_idxmdc, m_probph2kpi);
172 status = m_tuple1->addIndexedItem("normph2kpi", m_idxmdc, m_normph2kpi);
173 status = m_tuple1->addIndexedItem("pidnum2kpi", m_idxmdc, m_pidnum2kpi);
174 status = m_tuple1->addIndexedItem("bjmucjs", m_idxmdc, m_bjmucjs);
175 status = m_tuple1->addIndexedItem("bjmuc2kpi", m_idxmdc, m_bjmuc2kpi);
176 status = m_tuple1->addIndexedItem("bjemcjs", m_idxmdc, m_bjemcjs);
177 status = m_tuple1->addIndexedItem("bjemc2kpi", m_idxmdc, m_bjemc2kpi);
178 status = m_tuple1->addIndexedItem("bjtofjs", m_idxmdc, m_bjtofjs);
179 status = m_tuple1->addIndexedItem("bjtof2kpi", m_idxmdc, m_bjtof2kpi);
180 status = m_tuple1->addIndexedItem("bjtofvaljs", m_idxmdc, m_bjtofvaljs);
181 status = m_tuple1->addIndexedItem("bjtofval2kpi", m_idxmdc, m_bjtofval2kpi);
182
183 status = m_tuple1->addIndexedItem("emcjs", m_idxmdc, m_emcjs);
184 status = m_tuple1->addIndexedItem("evpjs", m_idxmdc, m_evpjs);
185 status = m_tuple1->addIndexedItem("timecgjs", m_idxmdc, m_timecgjs);
186 status = m_tuple1->addIndexedItem("depthjs", m_idxmdc, m_depthmucjs);
187 status = m_tuple1->addIndexedItem("layermucjs", m_idxmdc, m_layermucjs);
188
189 status = m_tuple1->addIndexedItem("emc2kpi", m_idxmdc, m_emc2kpi);
190 status = m_tuple1->addIndexedItem("evp2kpi", m_idxmdc, m_evp2kpi);
191 status = m_tuple1->addIndexedItem("timecg2kpi", m_idxmdc, m_timecg2kpi);
192 status = m_tuple1->addIndexedItem("depth2kpi", m_idxmdc, m_depthmuc2kpi);
193 status = m_tuple1->addIndexedItem("layermuc2kpi", m_idxmdc, m_layermuc2kpi);
194
195 status = m_tuple1->addIndexedItem("cotof1js", m_idxmdc, m_cotof1js);
196 status = m_tuple1->addIndexedItem("cotof2js", m_idxmdc, m_cotof2js);
197 status = m_tuple1->addIndexedItem("counterjs", m_idxmdc, m_counterjs);
198 status = m_tuple1->addIndexedItem("barreljs", m_idxmdc, m_barreljs);
199 status = m_tuple1->addIndexedItem("layertofjs", m_idxmdc, m_layertofjs);
200 status = m_tuple1->addIndexedItem("readoutjs", m_idxmdc, m_readoutjs);
201 status = m_tuple1->addIndexedItem("clusterjs", m_idxmdc, m_clusterjs);
202 status = m_tuple1->addIndexedItem("betajs", m_idxmdc, m_betajs);
203 status = m_tuple1->addIndexedItem("tofjs", m_idxmdc, m_tofjs);
204 status = m_tuple1->addIndexedItem("tofpathjs", m_idxmdc, m_tofpathjs);
205 status = m_tuple1->addIndexedItem("zhitjs", m_idxmdc, m_zhitjs);
206 status = m_tuple1->addIndexedItem("tofIDjs", m_idxmdc, m_tofIDjs);
207 status = m_tuple1->addIndexedItem("clusterIDjs", m_idxmdc, m_clusterIDjs);
208 status = m_tuple1->addIndexedItem("texejs", m_idxmdc, m_texejs);
209 status = m_tuple1->addIndexedItem("texmujs", m_idxmdc, m_texmujs);
210 status = m_tuple1->addIndexedItem("texpijs", m_idxmdc, m_texpijs);
211 status = m_tuple1->addIndexedItem("texkjs", m_idxmdc, m_texkjs);
212 status = m_tuple1->addIndexedItem("texprjs", m_idxmdc, m_texprjs);
213 status = m_tuple1->addIndexedItem("dtejs", m_idxmdc, m_dtejs);
214 status = m_tuple1->addIndexedItem("dtmujs", m_idxmdc, m_dtmujs);
215 status = m_tuple1->addIndexedItem("dtpijs", m_idxmdc, m_dtpijs);
216 status = m_tuple1->addIndexedItem("dtkjs", m_idxmdc, m_dtkjs);
217 status = m_tuple1->addIndexedItem("dtprjs", m_idxmdc, m_dtprjs);
218 status = m_tuple1->addIndexedItem("sigmaetofjs", m_idxmdc, m_sigmaetofjs);
219 status = m_tuple1->addIndexedItem("sigmamutofjs", m_idxmdc, m_sigmamutofjs);
220 status = m_tuple1->addIndexedItem("sigmapitofjs", m_idxmdc, m_sigmapitofjs);
221 status = m_tuple1->addIndexedItem("sigmaktofjs", m_idxmdc, m_sigmaktofjs);
222 status = m_tuple1->addIndexedItem("sigmaprtofjs", m_idxmdc, m_sigmaprtofjs);
223 status = m_tuple1->addIndexedItem("t0tofjs", m_idxmdc,m_t0tofjs);
224 status = m_tuple1->addIndexedItem("errt0tofjs", m_idxmdc,m_errt0tofjs);
225 status = m_tuple1->addIndexedItem("cotof12kpi", m_idxmdc, m_cotof12kpi);
226 status = m_tuple1->addIndexedItem("cotof22kpi", m_idxmdc, m_cotof22kpi);
227 status = m_tuple1->addIndexedItem("counter2kpi", m_idxmdc, m_counter2kpi);
228 status = m_tuple1->addIndexedItem("barrel2kpi", m_idxmdc, m_barrel2kpi);
229 status = m_tuple1->addIndexedItem("layertof2kpi", m_idxmdc, m_layertof2kpi);
230 status = m_tuple1->addIndexedItem("readout2kpi", m_idxmdc, m_readout2kpi);
231 status = m_tuple1->addIndexedItem("cluster2kpi", m_idxmdc, m_cluster2kpi);
232 status = m_tuple1->addIndexedItem("beta2kpi", m_idxmdc, m_beta2kpi);
233 status = m_tuple1->addIndexedItem("tof2kpi", m_idxmdc, m_tof2kpi);
234 status = m_tuple1->addIndexedItem("tofpath2kpi", m_idxmdc, m_tofpath2kpi);
235 status = m_tuple1->addIndexedItem("zhit2kpi", m_idxmdc, m_zhit2kpi);
236 status = m_tuple1->addIndexedItem("tofID2kpi", m_idxmdc, m_tofID2kpi);
237 status = m_tuple1->addIndexedItem("clusterID2kpi", m_idxmdc, m_clusterID2kpi);
238 status = m_tuple1->addIndexedItem("texe2kpi", m_idxmdc, m_texe2kpi);
239 status = m_tuple1->addIndexedItem("texmu2kpi", m_idxmdc, m_texmu2kpi);
240 status = m_tuple1->addIndexedItem("texpi2kpi", m_idxmdc, m_texpi2kpi);
241 status = m_tuple1->addIndexedItem("texk2kpi", m_idxmdc, m_texk2kpi);
242 status = m_tuple1->addIndexedItem("texpr2kpi", m_idxmdc, m_texpr2kpi);
243 status = m_tuple1->addIndexedItem("dte2kpi", m_idxmdc, m_dte2kpi);
244 status = m_tuple1->addIndexedItem("dtmu2kpi", m_idxmdc, m_dtmu2kpi);
245 status = m_tuple1->addIndexedItem("dtpi2kpi", m_idxmdc, m_dtpi2kpi);
246 status = m_tuple1->addIndexedItem("dtk2kpi", m_idxmdc, m_dtk2kpi);
247 status = m_tuple1->addIndexedItem("dtpr2kpi", m_idxmdc, m_dtpr2kpi);
248 status = m_tuple1->addIndexedItem("costpid2kpi", m_idxmdc, m_costpid2kpi);
249 status = m_tuple1->addIndexedItem("dedxpid2kpi", m_idxmdc, m_dedxpid2kpi);
250 status = m_tuple1->addIndexedItem("tof1pid2kpi", m_idxmdc, m_tof1pid2kpi);
251 status = m_tuple1->addIndexedItem("tof2pid2kpi", m_idxmdc, m_tof2pid2kpi);
252 status = m_tuple1->addIndexedItem("probe2kpi", m_idxmdc, m_probe2kpi);
253 status = m_tuple1->addIndexedItem("probmu2kpi", m_idxmdc, m_probmu2kpi);
254 status = m_tuple1->addIndexedItem("probpi2kpi", m_idxmdc, m_probpi2kpi);
255 status = m_tuple1->addIndexedItem("probk2kpi", m_idxmdc, m_probk2kpi);
256 status = m_tuple1->addIndexedItem("probpr2kpi", m_idxmdc, m_probpr2kpi);
257
258 status = m_tuple1->addIndexedItem("chipidxpid2kpi", m_idxmdc, m_chipidxpid2kpi);
259 status = m_tuple1->addIndexedItem("chipitof1pid2kpi", m_idxmdc, m_chipitof1pid2kpi);
260 status = m_tuple1->addIndexedItem("chipitof2pid2kpi", m_idxmdc, m_chipitof2pid2kpi);
261 status = m_tuple1->addIndexedItem("chipitofpid2kpi", m_idxmdc, m_chipitofpid2kpi);
262 status = m_tuple1->addIndexedItem("chipitofepid2kpi", m_idxmdc, m_chipitofepid2kpi);
263 status = m_tuple1->addIndexedItem("chipitofqpid2kpi", m_idxmdc, m_chipitofqpid2kpi);
264 status = m_tuple1->addIndexedItem("probpidxpid2kpi", m_idxmdc, m_probpidxpid2kpi);
265 status = m_tuple1->addIndexedItem("probpitofpid2kpi", m_idxmdc, m_probpitofpid2kpi);
266 status = m_tuple1->addIndexedItem("chikdxpid2kpi", m_idxmdc, m_chikdxpid2kpi);
267 status = m_tuple1->addIndexedItem("chiktof1pid2kpi", m_idxmdc, m_chiktof1pid2kpi);
268 status = m_tuple1->addIndexedItem("chiktof2pid2kpi", m_idxmdc, m_chiktof2pid2kpi);
269 status = m_tuple1->addIndexedItem("chiktofpid2kpi", m_idxmdc, m_chiktofpid2kpi);
270 status = m_tuple1->addIndexedItem("chiktofepid2kpi", m_idxmdc, m_chiktofepid2kpi);
271 status = m_tuple1->addIndexedItem("chiktofqpid2kpi", m_idxmdc, m_chiktofqpid2kpi);
272 status = m_tuple1->addIndexedItem("probkdxpid2kpi", m_idxmdc, m_probkdxpid2kpi);
273 status = m_tuple1->addIndexedItem("probktofpid2kpi", m_idxmdc, m_probktofpid2kpi);
274
275 status = m_tuple1->addIndexedItem("chiprdxpid2kpi", m_idxmdc, m_chiprdxpid2kpi);
276 status = m_tuple1->addIndexedItem("chiprtof1pid2kpi", m_idxmdc, m_chiprtof1pid2kpi);
277 status = m_tuple1->addIndexedItem("chiprtof2pid2kpi", m_idxmdc, m_chiprtof2pid2kpi);
278 status = m_tuple1->addIndexedItem("chiprtofpid2kpi", m_idxmdc, m_chiprtofpid2kpi);
279 status = m_tuple1->addIndexedItem("chiprtofepid2kpi", m_idxmdc, m_chiprtofepid2kpi);
280 status = m_tuple1->addIndexedItem("chiprtofqpid2kpi", m_idxmdc, m_chiprtofqpid2kpi);
281 status = m_tuple1->addIndexedItem("probprdxpid2kpi", m_idxmdc, m_probprdxpid2kpi);
282 status = m_tuple1->addIndexedItem("probprtofpid2kpi", m_idxmdc, m_probprtofpid2kpi);
283
284 status = m_tuple1->addIndexedItem("cosmdcjs", m_idxmdc, m_cosmdcjs);
285 status = m_tuple1->addIndexedItem("phimdcjs", m_idxmdc, m_phimdcjs);
286 status = m_tuple1->addIndexedItem("cosmdc2kpi", m_idxmdc, m_cosmdc2kpi);
287 status = m_tuple1->addIndexedItem("phimdc2kpi", m_idxmdc, m_phimdc2kpi);
288
289 status = m_tuple1->addIndexedItem("dedxpidjs", m_idxmdc, m_dedxpidjs);
290 status = m_tuple1->addIndexedItem("tof1pidjs", m_idxmdc, m_tof1pidjs);
291 status = m_tuple1->addIndexedItem("tof2pidjs", m_idxmdc, m_tof2pidjs);
292 status = m_tuple1->addIndexedItem("probejs", m_idxmdc, m_probejs);
293 status = m_tuple1->addIndexedItem("probmujs", m_idxmdc, m_probmujs);
294 status = m_tuple1->addIndexedItem("probpijs", m_idxmdc, m_probpijs);
295 status = m_tuple1->addIndexedItem("probkjs", m_idxmdc, m_probkjs);
296 status = m_tuple1->addIndexedItem("probprjs", m_idxmdc, m_probprjs);
297 status = m_tuple1->addItem ("mchic2kpi", m_mchic2kpi);
298 status = m_tuple1->addItem ("mpsip2kpi", m_mpsip2kpi);
299 status = m_tuple1->addItem ("chis2kpi", m_chis2kpi);
300 status = m_tuple1->addItem ("mchic4c2kpi", m_mchic4c2kpi);
301 status = m_tuple1->addItem ("mpsip4c2kpi", m_mpsip4c2kpi);
302 status = m_tuple1->addItem ("chis4c2kpi", m_chis4c2kpi);
303
304 status = m_tuple1->addItem("indexemc", m_idxemc, 0, 5000);
305 status = m_tuple1->addIndexedItem("numHits", m_idxemc, m_numHits);
306 status = m_tuple1->addIndexedItem("secmom", m_idxemc, m_secondmoment);
307 status = m_tuple1->addIndexedItem("latmom", m_idxemc, m_latmoment);
308 status = m_tuple1->addIndexedItem("timegm", m_idxemc, m_timegm);
309 status = m_tuple1->addIndexedItem("cellId", m_idxemc, m_cellId);
310 status = m_tuple1->addIndexedItem("module", m_idxemc, m_module);
311 status = m_tuple1->addIndexedItem("a20Moment", m_idxemc, m_a20Moment);
312 status = m_tuple1->addIndexedItem("a42Moment", m_idxemc, m_a42Moment);
313 status = m_tuple1->addIndexedItem("getEAll", m_idxemc, m_getEAll);
314 status = m_tuple1->addIndexedItem("getShowerId", m_idxemc, m_getShowerId);
315 status = m_tuple1->addIndexedItem("getClusterId", m_idxemc, m_getClusterId);
316 status = m_tuple1->addIndexedItem("x", m_idxemc, m_x);
317 status = m_tuple1->addIndexedItem("y", m_idxemc, m_y);
318 status = m_tuple1->addIndexedItem("z", m_idxemc, m_z);
319 status = m_tuple1->addIndexedItem("cosemc", m_idxemc, m_cosemc);
320 status = m_tuple1->addIndexedItem("phiemc", m_idxemc, m_phiemc);
321 status = m_tuple1->addIndexedItem("energy", m_idxemc, m_energy);
322 status = m_tuple1->addIndexedItem("e1", m_idxemc, m_eSeed);
323 status = m_tuple1->addIndexedItem("e9", m_idxemc, m_e3x3);
324 status = m_tuple1->addIndexedItem("e25", m_idxemc, m_e5x5);
325 status = m_tuple1->addIndexedItem("dang4c", m_idxemc, m_dang4c);
326 status = m_tuple1->addIndexedItem("dthe4c", m_idxemc, m_dthe4c);
327 status = m_tuple1->addIndexedItem("dphi4c", m_idxemc, m_dphi4c);
328 status = m_tuple1->addIndexedItem("dang4crt", m_idxemc, m_dang4crt);
329 status = m_tuple1->addIndexedItem("dthe4crt", m_idxemc, m_dthe4crt);
330 status = m_tuple1->addIndexedItem("dphi4crt", m_idxemc, m_dphi4crt);
331 status = m_tuple1->addIndexedItem("phtof", m_idxemc, 3, m_phgmtof,0.0,10000.0);
332 status = m_tuple1->addIndexedItem("phgmtof0", m_idxemc, m_phgmtof0);
333 status = m_tuple1->addIndexedItem("phgmtof1", m_idxemc, m_phgmtof1);
334 status = m_tuple1->addIndexedItem("phgmtof2", m_idxemc, m_phgmtof2);
335
336 }
337 else {
338 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
339 return StatusCode::FAILURE;
340 }
341 }
342}
343
344StatusCode sc;
345
346 if(m_skim4pi)
347 {
348 sc = createSubAlgorithm( "EventWriter", "Selectgam4pi", m_subalg1);
349 if( sc.isFailure() ) {
350 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4pi" <<endreq;
351 return sc;
352 } else {
353 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4pi" <<endreq;
354 }
355 }
356
357
358 if(m_skim4k)
359 {
360 sc = createSubAlgorithm( "EventWriter", "Selectgam4k", m_subalg2);
361 if( sc.isFailure() ) {
362 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam4k" <<endreq;
363 return sc;
364 } else {
365 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam4k" <<endreq;
366 }
367 }
368
369 if(m_skim2pi2pr)
370 {
371 sc = createSubAlgorithm( "EventWriter", "Selectgam2pi2pr", m_subalg3);
372 if( sc.isFailure() ) {
373 log << MSG::ERROR << "Error creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
374 return sc;
375 } else {
376 log << MSG::INFO << "Success creating Sub-Algorithm Selectgam2pi2pr" <<endreq;
377 }
378 }
379
380
381
382
383 //
384 //--------end of book--------
385 //
386
387 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
388 return StatusCode::SUCCESS;
389
390}
391
392// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
393StatusCode Gam4pikp::execute() {
394
395// std::cout << "execute()" << std::endl;
396 setFilterPassed(false);
397
398
399 MsgStream log(msgSvc(), name());
400 log << MSG::INFO << "in execute()" << endreq;
401
402 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
403 int run=eventHeader->runNumber();
404 int event=eventHeader->eventNumber();
405 log << MSG::DEBUG <<"run, evtnum = "
406 << run << " , "
407 << event <<endreq;
408// cout<<"run "<<run<<endl;
409// cout<<"event "<<event<<endl;
410if(m_rootput)
411{
412 m_run=run;
413 m_rec=event;
414}
415 Ncut[0]++;
416 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
417 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
418 << evtRecEvent->totalCharged() << " , "
419 << evtRecEvent->totalNeutral() << " , "
420 << evtRecEvent->totalTracks() <<endreq;
421
422 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
423
424 if(m_rootput)
425{
426 Gam4pikp::InitVar();
427}
428
429 // for MC topo analysis
430if(m_rootput)
431{
432 if (eventHeader->runNumber()<0)
433 {
434 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
435 int m_numParticle = 0;
436 if (!mcParticleCol)
437 {
438// std::cout << "Could not retrieve McParticelCol" << std::endl;
439 return StatusCode::FAILURE;
440 }
441 else
442 {
443 bool psipDecay = false;
444 int rootIndex = -1;
445 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
446 for (; iter_mc != mcParticleCol->end(); iter_mc++)
447 {
448 if ((*iter_mc)->primaryParticle()) continue;
449 if (!(*iter_mc)->decayFromGenerator()) continue;
450 if ((*iter_mc)->particleProperty()==100443)
451 {
452 psipDecay = true;
453 rootIndex = (*iter_mc)->trackIndex();
454 }
455 if (!psipDecay) continue;
456 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
457 int pdgid = (*iter_mc)->particleProperty();
458 m_pdgid[m_numParticle] = pdgid;
459 m_motheridx[m_numParticle] = mcidx;
460 m_numParticle += 1;
461 }
462 }
463 m_idxmc = m_numParticle;
464 }
465}
466
467 Vint iGood, ipip, ipim;
468 iGood.clear();
469 ipip.clear();
470 ipim.clear();
471 Vp4 ppip, ppim;
472 ppip.clear();
473 ppim.clear();
474 int nCharge = 0;
475 int ngdcgx=0;
476 int ncgx=0;
477 Vdouble pTrkCh;
478 pTrkCh.clear();
479 Hep3Vector v(0,0,0);
480 Hep3Vector vv(0,0,0);
481
482 IVertexDbSvc* vtxsvc;
483 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
484 if(vtxsvc->isVertexValid())
485 {
486 double* vertex = vtxsvc->PrimaryVertex();
487 double* vertexsigma = vtxsvc->SigmaPrimaryVertex();
488 v.setX(vertex[0]);
489 v.setY(vertex[1]);
490 v.setZ(vertex[2]);
491 vv.setX(vertexsigma[0]);
492 vv.setY(vertexsigma[1]);
493 vv.setZ(vertexsigma[2]);
494 }
495
496
497 if(run<0)
498 {
499 v[0]=0.;
500 v[1]=0.;
501 v[2]=0.;
502 vv[0]=0.;
503 vv[1]=0.;
504 vv[2]=0.;
505
506 }
507
508
509 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
510 if(i>=evtRecTrkCol->size()) break;
511 ncgx++;
512 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
513 if(!(*itTrk)->isMdcTrackValid()) continue;
514 if (!(*itTrk)->isMdcKalTrackValid()) continue;
515 ngdcgx++;
516 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
517
518 HepVector a = mdcTrk->helix();
519 HepSymMatrix Ea = mdcTrk->err();
520 HepPoint3D pivot(0.,0.,0.);
521 HepPoint3D IP(v[0],v[1],v[2]);
522 VFHelix helixp(pivot,a,Ea);
523 helixp.pivot(IP);
524 HepVector vec = helixp.a();
525 pTrkCh.push_back(mdcTrk->p()*mdcTrk->charge());
526 iGood.push_back((*itTrk)->trackId());
527 nCharge += mdcTrk->charge();
528 }
529
530
531 int nGood = iGood.size();
532 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
533 if((nGood < 4))
534 {
535 return StatusCode::SUCCESS;
536 }
537
538 Ncut[1]++;
539 Vint iGam;
540 int ngamx=0;
541 int ngdgamx=0;
542 iGam.clear();
543 Vint iGamnofit;
544 iGamnofit.clear();
545 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
546 ngamx++;
547 if(i>=evtRecTrkCol->size()) break;
548 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
549 if(!(*itTrk)->isEmcShowerValid()) continue;
550 RecEmcShower *emcTrk = (*itTrk)->emcShower();
551 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
552 // find the nearest charged track
553// double dthe = 200.;
554// double dphi = 200.;
555// double dang = 200.;
556 double dthe = 200.;
557 double dphi = 200.;
558 double dang = 200.;
559
560 ngdgamx++;
561 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
562 if(j>=evtRecTrkCol->size()) break;
563 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
564 if(!(*jtTrk)->isExtTrackValid()) continue;
565 RecExtTrack *extTrk = (*jtTrk)->extTrack();
566 if(extTrk->emcVolumeNumber() == -1) continue;
567 Hep3Vector extpos = extTrk->emcPosition();
568 // double ctht = extpos.cosTheta(emcpos);
569 double angd = extpos.angle(emcpos);
570 double thed = extpos.theta() - emcpos.theta();
571 double phid = extpos.deltaPhi(emcpos);
572 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
573 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
574
575// if(fabs(thed) < fabs(dthe)) dthe = thed;
576// if(fabs(phid) < fabs(dphi)) dphi = phid;
577// if(angd < dang) dang = angd;
578
579 if(angd < dang){ dang=angd;
580 dthe=thed;
581 dphi=phid;
582 }
583
584 }
585 // if(dang>=200) continue;
586 double eraw = emcTrk->energy();
587 dthe = dthe * 180 / (CLHEP::pi);
588 dphi = dphi * 180 / (CLHEP::pi);
589 dang = dang * 180 / (CLHEP::pi);
590// dthert = dthert * 180 / (CLHEP::pi);
591// dphirt = dphirt * 180 / (CLHEP::pi);
592// dangrt = dangrt * 180 / (CLHEP::pi);
593 // if(eraw < m_energyThreshold) continue;
594 iGam.push_back((*itTrk)->trackId());
595 if(eraw < m_energyThreshold) continue;
596 if(dang < 20.0) continue;
597 iGamnofit.push_back((*itTrk)->trackId());
598
599 }
600 // Finish Good Photon Selection
601 //
602 int nGam = iGam.size();
603 int nGamnofit = iGamnofit.size();
604
605 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
606 if(nGam<1){
607 return StatusCode::SUCCESS;
608 }
609 Ncut[2]++;
610
611 if(nGood>20||nGam>30) return StatusCode::SUCCESS;
612
613if(m_rootput)
614{
615m_idxmdc = nGood;
616m_idxemc=nGam;
617}
618//
619 Vp4 pGam;
620 pGam.clear();
621 for(int i = 0; i < nGam; i++) {
622 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
623 RecEmcShower* emcTrk = (*itTrk)->emcShower();
624 double eraw = emcTrk->energy();
625 double phi = emcTrk->phi();
626 double the = emcTrk->theta();
627 HepLorentzVector ptrk;
628 ptrk.setPx(eraw*sin(the)*cos(phi));
629 ptrk.setPy(eraw*sin(the)*sin(phi));
630 ptrk.setPz(eraw*cos(the));
631 ptrk.setE(eraw);
632
633 // ptrk = ptrk.boost(-0.011,0,0);// boost to cms
634
635 pGam.push_back(ptrk);
636 }
637 //
638 // Assign 4-momentum to each charged track
639 //
640 ParticleID *pid = ParticleID::instance();
641
642
643 for(int i = 0; i < nGood; i++) {
644 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
645 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
646
647 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use RecMdcKalTrack substitute RecMdcTrack
648 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default setting is pion
649
650 if(mdcKalTrk->charge() >0) {
651 ipip.push_back(iGood[i]);
652 HepLorentzVector ptrk;
653 ptrk.setPx(mdcKalTrk->px());
654 ptrk.setPy(mdcKalTrk->py());
655 ptrk.setPz(mdcKalTrk->pz());
656 double p3 = ptrk.mag();
657 ptrk.setE(sqrt(p3*p3+mpi*mpi));
658
659
660 ppip.push_back(ptrk);
661 } else {
662 ipim.push_back(iGood[i]);
663 HepLorentzVector ptrk;
664 ptrk.setPx(mdcKalTrk->px());
665 ptrk.setPy(mdcKalTrk->py());
666 ptrk.setPz(mdcKalTrk->pz());
667 double p3 = ptrk.mag();
668 ptrk.setE(sqrt(p3*p3+mpi*mpi));
669
670 // ptrk = ptrk.boost(-0.011,0,0);//boost to cms
671
672 ppim.push_back(ptrk);
673 }
674 }
675 int npip = ipip.size();
676 int npim = ipim.size();
677 // if(npip*npim != 4) return SUCCESS;
678 if((npip < 2)||(npim < 2)) return SUCCESS;
679
680 Ncut[3]++;
681
682 //
683 // Loop each gamma pair, check ppi0 and pTot
684 //
685
686 int icgp1=-1;
687 int icgp2=-1;
688 int icgm1=-1;
689 int icgm2=-1;
690 int igam=-1;
691 double chisqtrk = 9999.;
692 WTrackParameter wcgp1;
693 WTrackParameter wcgp2;
694 WTrackParameter wcgm1;
695 WTrackParameter wcgm2;
696 int ipip1js=-1;
697 int ipip2js=-1;
698 int ipim1js=-1;
699 int ipim2js=-1;
700 double mcompall=9999;
701 double mppreclst=9999;
702 double meelst=9999;;
703 double mchic2kpilst=9999;
704 double chis4c2kpilst=9999;
705 int type2kpilst=0;
706 double dtpr2kpilst[4]={9999,9999,9999,9999};
707
708 double mc1=mpi,mc2=mpi,mc3=mpi,mc4=mpi;
709 double mchic01=0.0;
710 HepLorentzVector pgam=(0,0,0,0);
711 HepPoint3D xorigin(0.,0.,0.);
712 xorigin[0]=9999.;
713 xorigin[1]=9999.;
714 xorigin[2]=9999.;
715 HepSymMatrix xem(3,0);
716
717 int cl4pi=0;
718 int clajs=0;
719 HepLorentzVector p4psipjs(0.011*m_ecms, 0.0, 0.0, m_ecms);
720 double psipBetajs = (p4psipjs.vect()).mag()/(p4psipjs.e());
721 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.0, m_ecms);
722 double psipBeta = (p4psip.vect()).mag()/(p4psip.e()); // beta = P/E
723
724 for(int ii = 0; ii < npip-1; ii++) {
725 RecMdcKalTrack *pip1Trk = (*(evtRecTrkCol->begin()+ipip[ii]))->mdcKalTrack();
726 for(int ij = ii+1; ij < npip; ij++) {
727 RecMdcKalTrack *pip2Trk = (*(evtRecTrkCol->begin()+ipip[ij]))->mdcKalTrack();
728 for(int ik = 0; ik < npim-1; ik++) {
729 RecMdcKalTrack *pim1Trk = (*(evtRecTrkCol->begin()+ipim[ik]))->mdcKalTrack();
730 for(int il = ik+1; il < npim; il++) {
731 RecMdcKalTrack *pim2Trk = (*(evtRecTrkCol->begin()+ipim[il]))->mdcKalTrack();
732 double squar[3]={9999.,9999.,9999.};
733 double squarkpi[6]={9999.,9999.,9999.,9999.,9999.,9999.};
734 WTrackParameter wvpip1Trk, wvpim1Trk, wvpip2Trk, wvpim2Trk;
735
736 Vint iGoodjs;
737 iGoodjs.clear();
738 Vdouble pTrkjs;
739 pTrkjs.clear();
740
741
742
743 pTrkjs.push_back(pip1Trk->p()*pip1Trk->charge());
744 pTrkjs.push_back(pip2Trk->p()*pip2Trk->charge());
745 pTrkjs.push_back(pim1Trk->p()*pim1Trk->charge());
746 pTrkjs.push_back(pim2Trk->p()*pim2Trk->charge());
747 iGoodjs.push_back(ipip[ii]);
748 iGoodjs.push_back(ipip[ij]);
749 iGoodjs.push_back(ipim[ik]);
750 iGoodjs.push_back(ipim[il]);
751
752 Gam4pikp::BubbleSort(pTrkjs, iGoodjs);
753 Vint jGoodjs;
754 jGoodjs.clear();
755 Vp4 p4chTrkjs;
756 p4chTrkjs.clear();
757 Vint i4cpip1js, i4cpip2js, i4cpim1js, i4cpim2js;
758 i4cpip1js.clear();
759 i4cpip2js.clear();
760 i4cpim1js.clear();
761 i4cpim2js.clear();
762 i4cpip1js.push_back(iGoodjs[2]);
763 i4cpip2js.push_back(iGoodjs[3]);
764 i4cpim1js.push_back(iGoodjs[1]);
765 i4cpim2js.push_back(iGoodjs[0]);
766 jGoodjs.push_back(i4cpip1js[0]);
767 jGoodjs.push_back(i4cpim1js[0]);
768 jGoodjs.push_back(i4cpip2js[0]);
769 jGoodjs.push_back(i4cpim2js[0]);
770
771 for (int i = 0; i < 4; i++)
772 {
773 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGoodjs[i];
774 if (!(*itTrk)->isMdcTrackValid()) continue;
775 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
776 if (!(*itTrk)->isMdcKalTrackValid()) continue;
777 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
778 double ptrk;
779 HepLorentzVector p4trk;
780 if (i<2)
781 {
782 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
783 ptrk = mdcKalTrk->p();
784 p4trk.setPx(mdcKalTrk->px());
785 p4trk.setPy(mdcKalTrk->py());
786 p4trk.setPz(mdcKalTrk->pz());
787 p4trk.setE(sqrt(ptrk*ptrk+mpi*mpi));
788 }
789 else{
790 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
791 ptrk = mdcKalTrk->p();
792 p4trk.setPx(mdcKalTrk->px());
793 p4trk.setPy(mdcKalTrk->py());
794 p4trk.setPz(mdcKalTrk->pz());
795 p4trk.setE(sqrt(ptrk*ptrk+xmass[0]*xmass[0]));
796 }
797 p4trk.boost(-1.0*psipBetajs, 0.0, 0.0);
798 p4chTrkjs.push_back(p4trk);
799 }
800 p4psipjs.boost(-1.0*psipBetajs, 0.0, 0.0);
801
802 HepLorentzVector p4pipijs = p4chTrkjs[0] + p4chTrkjs[1];
803 HepLorentzVector p4eejs = p4chTrkjs[2] + p4chTrkjs[3];
804 HepLorentzVector p4pipirecjs = p4psipjs - p4pipijs;
805 double mpprecjs=p4pipirecjs.m();
806 double mpipijs=p4pipijs.m();
807 double meejs=p4eejs.m();
808 double mcomp=sqrt((mpprecjs-3.097)*(mpprecjs-3.097)+(meejs-3.097)*(meejs-3.097));
809 if(mcomp<mcompall)
810 {
811 mcompall=mcomp;
812 ipip1js=i4cpip1js[0];
813 ipim1js=i4cpim1js[0];
814 ipip2js=i4cpip2js[0];
815 ipim2js=i4cpim2js[0];
816 mppreclst=mpprecjs;
817 meelst=meejs;
818 }
819
820if(m_rootput)
821{
822 m_mpprecall=mppreclst;
823 m_meeall=meelst;
824}
825 }
826
827 }
828 }
829 }
830
831 //{
832// HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.0, m_ecms);
833// double psipBeta = (p4psip.vect()).mag()/(p4psip.e()); // beta = P/E
834 Vint iGood4c;
835 iGood4c.clear();
836 Vdouble pTrk4c;
837 pTrk4c.clear();
838 Vint jGood;
839 jGood.clear();
840 Vint jsGood;
841 jsGood.clear();
842
843 if(mcompall>9997)
844 return StatusCode::SUCCESS;
845
846 jsGood.push_back(ipip1js);
847 jsGood.push_back(ipim1js);
848 jsGood.push_back(ipip2js);
849 jsGood.push_back(ipim2js);
850
851 for(int i = 0; i < evtRecEvent->totalCharged(); i++)
852 {
853 if(i>=evtRecTrkCol->size()) break;
854 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
855 if(!(*itTrk)->isMdcTrackValid()) continue;
856 if (!(*itTrk)->isMdcKalTrackValid()) continue;
857 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
858 if((i!=ipip1js)&&(i!=ipim1js)&&(i!=ipip2js)&&(i!=ipim2js))
859 {
860 jsGood.push_back(i);
861 }
862 }
863
864 int njsGood=jsGood.size();
865 //ParticleID *pid = ParticleID::instance();
866if(m_rootput)
867{
868 for (int i = 0; i < njsGood; i++)
869 {
870 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jsGood[i];
871 if (!(*itTrk)->isMdcTrackValid()) continue;
872 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
873 if (!(*itTrk)->isMdcKalTrackValid()) continue;
874 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
875 double ptrk;
876 ptrk = mdcKalTrk->p();
877 m_x0js[i] = mdcTrk->x();
878 m_y0js[i] = mdcTrk->y();
879 m_z0js[i] = mdcTrk->z();
880 m_r0js[i] = mdcTrk->r();
881 m_ppmdcjs[i] = mdcTrk->p();
882 m_pxmdcjs[i] = mdcTrk->px();
883 m_pymdcjs[i] = mdcTrk->py();
884 m_pzmdcjs[i] = mdcTrk->pz();
885 m_ppkaljs[i] = mdcKalTrk->p();
886 Hep3Vector p3jsi(mdcTrk->px(),mdcTrk->py(),mdcTrk->pz());
887 if(njsGood>4){
888 EvtRecTrackIterator itTrk5 = evtRecTrkCol->begin() + jsGood[4];
889 RecMdcTrack *mdcTrk5 = (*itTrk5)->mdcTrack();
890 Hep3Vector p3js5(mdcTrk5->px(),mdcTrk5->py(),mdcTrk5->pz());
891 m_angjs5[i]=p3jsi.angle(p3js5);
892 m_nearjs5[i]=p3jsi.howNear(p3js5);
893 }
894 if(njsGood>5){
895 EvtRecTrackIterator itTrk6 = evtRecTrkCol->begin() + jsGood[5];
896 RecMdcTrack *mdcTrk6 = (*itTrk6)->mdcTrack();
897 Hep3Vector p3js6(mdcTrk6->px(),mdcTrk6->py(),mdcTrk6->pz());
898 m_angjs6[i]=p3jsi.angle(p3js6);
899 m_nearjs6[i]=p3jsi.howNear(p3js6);
900 }
901
902
903 m_ptmdcjs[i] = mdcTrk->pxy();
904 m_ptkaljs[i] = mdcKalTrk->pxy();
905 double x0=mdcTrk->x();
906 double y0=mdcTrk->y();
907 double z0=mdcTrk->z();
908 double phi0=mdcTrk->helix(1);
909 double xv=v[0];
910 double yv=v[1];
911 double zv=v[2];
912 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
913 double Rz=z0-zv;
914 m_Rxyjs[i]=Rxy;
915 m_Rzjs[i]=Rz;
916 HepVector a = mdcTrk->helix();
917 HepSymMatrix Ea = mdcTrk->err();
918 HepPoint3D pivot(0.,0.,0.);
919 HepPoint3D IP(v[0],v[1],v[2]);
920 VFHelix helixp(pivot,a,Ea);
921 helixp.pivot(IP);
922 HepVector vec = helixp.a();
923 m_Rnxyjs[i]=vec[0];
924 m_Rnzjs[i]=vec[3];
925 m_phinjs[i]=vec[1];
926 //tof
927 if ((*itTrk)->isTofTrackValid())
928 { m_bjtofvaljs[i]=1;
929 m_cotof1js[i]=1;
930 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
931 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
932 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
933 TofHitStatus *status = new TofHitStatus;
934 status->setStatus((*iter_tof)->status());
935 if(status->is_cluster()){
936 m_bjtofjs[i]=1;
937 m_counterjs[i] = status->is_counter();
938 m_barreljs[i] = status->is_barrel();
939 m_layertofjs[i] = status->layer();
940 m_readoutjs[i] = status->is_readout();
941 m_clusterjs[i] = status->is_cluster();
942 m_cotof2js[i]=2;
943 m_betajs[i]= (*iter_tof)->beta();
944 m_tofjs[i] = (*iter_tof)->tof();
945 m_tofpathjs[i] = (*iter_tof)->path();
946 m_zhitjs[i]= (*iter_tof)->zrhit();
947 m_texejs[i] = (*iter_tof)->texpElectron();
948 m_texmujs[i] = (*iter_tof)->texpMuon();
949 m_texpijs[i] = (*iter_tof)->texpPion();
950 m_texkjs[i] = (*iter_tof)->texpKaon();
951 m_texprjs[i] = (*iter_tof)->texpProton();
952 m_dtejs[i] = m_tofjs[i] - m_texejs[i];
953 m_dtmujs[i] = m_tofjs[i] - m_texmujs[i];
954 m_dtpijs[i] = m_tofjs[i] - m_texpijs[i];
955 m_dtkjs[i] = m_tofjs[i] - m_texkjs[i];
956 m_dtprjs[i] = m_tofjs[i] - m_texprjs[i];
957
958 m_sigmaetofjs[i]=(*iter_tof)->sigma(0);
959 m_sigmamutofjs[i]=(*iter_tof)->sigma(1);
960 m_sigmapitofjs[i]=(*iter_tof)->sigma(2);
961 m_sigmaktofjs[i]=(*iter_tof)->sigma(3);
962 m_sigmaprtofjs[i]=(*iter_tof)->sigma(4);
963 m_t0tofjs[i]=(*iter_tof)->t0();
964 m_errt0tofjs[i]=(*iter_tof)->errt0();
965
966 m_tofIDjs[i]=(*iter_tof)->tofID();
967 m_clusterIDjs[i]=(*iter_tof)->tofTrackID();
968 }
969 delete status;
970 }
971 }
972 //dedx
973 if ((*itTrk)->isMdcDedxValid())
974 {
975
976 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
977 m_chiejs[i] = dedxTrk->chiE();
978 m_chimujs[i] = dedxTrk->chiMu();
979 m_chipijs[i] = dedxTrk->chiPi();
980 m_chikjs[i] = dedxTrk->chiK();
981 m_chipjs[i] = dedxTrk->chiP();
982 m_ghitjs[i] = dedxTrk->numGoodHits();
983 m_thitjs[i] = dedxTrk->numTotalHits();
984 m_probphjs[i] = dedxTrk->probPH();
985 m_normphjs[i] = dedxTrk->normPH();
986 }
987
988
989
990 //emc
991 if ( (*itTrk)->isEmcShowerValid() )
992 {
993 RecEmcShower* emcTrk = (*itTrk)->emcShower();
994 m_bjemcjs[i]=1;
995 m_emcjs[i] = emcTrk->energy();
996 m_evpjs[i] = emcTrk->energy()/ptrk;
997 m_timecgjs[i] = emcTrk->time();
998 // totalEnergy += emcTrk->energy();
999 }
1000
1001 //muc
1002 if ( (*itTrk)->isMucTrackValid() )
1003 {
1004 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
1005 double dpthp = mucTrk->depth()/25.0;//why?
1006 // m_depthmucjs[i] = dpthp<10.0 ? dpthp : 10.0;
1007 m_bjmucjs[i]=1;
1008 m_depthmucjs[i]=mucTrk->depth();
1009 m_layermucjs[i] = mucTrk->numLayers();
1010 }
1011
1012 m_cosmdcjs[i] = cos(mdcTrk->theta());
1013 m_phimdcjs[i] = mdcTrk->phi();
1014 //pid
1015 pid->init();
1016 pid->setMethod(pid->methodProbability());
1017 pid->setChiMinCut(4);
1018 pid->setRecTrack(*itTrk);
1019 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
1020 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
1021 pid->identify(pid->onlyElectron() | pid->onlyMuon());
1022 pid->calculate();
1023 if(!(pid->IsPidInfoValid())) continue;
1024 // m_costpidjs[i] = cos(mdcTrk->theta());
1025 m_dedxpidjs[i] = pid->chiDedx(2);
1026 m_tof1pidjs[i] = pid->chiTof1(2);
1027 m_tof2pidjs[i] = pid->chiTof2(2);
1028 m_probejs[i] = pid->probElectron();
1029 m_probmujs[i] = pid->probMuon();
1030 m_probpijs[i] = pid->probPion();
1031 m_probkjs[i] = pid->probKaon();
1032 m_probprjs[i] = pid->probProton();
1033
1034
1035
1036 }
1037
1038}
1039
1040 Vint jGam2kpi;
1041 jGam2kpi.clear();
1042
1043 Vint iGood2kpi, ipip2kpi, ipim2kpi;
1044 iGood2kpi.clear();
1045 ipip2kpi.clear();
1046 ipim2kpi.clear();
1047 Vp4 ppip2kpi, ppim2kpi;
1048 ppip2kpi.clear();
1049 ppim2kpi.clear();
1050
1051 Vint ipipnofit, ipimnofit, ikpnofit, ikmnofit, ipropnofit, ipromnofit;
1052 ipipnofit.clear();
1053 ipimnofit.clear();
1054 ikpnofit.clear();
1055 ikmnofit.clear();
1056 ipropnofit.clear();
1057 ipromnofit.clear();
1058 Vp4 ppipnofit, ppimnofit, pkpnofit, pkmnofit, ppropnofit, ppromnofit;
1059 ppipnofit.clear();
1060 ppimnofit.clear();
1061 pkpnofit.clear();
1062 pkmnofit.clear();
1063 ppropnofit.clear();
1064 ppromnofit.clear();
1065
1066 Vdouble p3pip2kpi;
1067 p3pip2kpi.clear();
1068 Vdouble p3pim2kpi;
1069 p3pim2kpi.clear();
1070
1071 Vint itrak2kpi;
1072 itrak2kpi.clear();
1073 int cls2kpi;
1074 double chis2kpi=9999.;
1075 double m4piall=9999.;
1076 double m4kall=9999.;
1077 double mgam4piall=9999.;
1078 double mgam4kall=9999.;
1079 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
1080 if(i>=evtRecTrkCol->size()) break;
1081 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
1082 if(!(*itTrk)->isMdcTrackValid()) continue;
1083 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1084 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1085 double z02kpi = mdcTrk->z();
1086 double r02kpi = mdcTrk->r();
1087 HepVector a = mdcTrk->helix();
1088 HepSymMatrix Ea = mdcTrk->err();
1089 HepPoint3D pivot(0.,0.,0.);
1090 HepPoint3D IP(v[0],v[1],v[2]);
1091 VFHelix helixp(pivot,a,Ea);
1092 helixp.pivot(IP);
1093 HepVector vec = helixp.a();
1094 double Rnxy02kpi=vec[0];
1095 double Rnz02kpi=vec[3];
1096 // if(fabs(Rnxy02kpi>1.0)) continue;
1097 // if(fabs(Rnz02kpi>10.0)) continue;
1098 iGood2kpi.push_back((*itTrk)->trackId());
1099 }
1100 int nGood2kpi = iGood2kpi.size();
1101 if(nGood2kpi==4)
1102 {
1103 for(int i = 0; i < nGood2kpi; i++) {
1104 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood2kpi[i];
1105 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1106 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
1107 if(mdcKalTrk->charge() >0) {
1108 ipip2kpi.push_back(iGood2kpi[i]);
1109 p3pip2kpi.push_back(mdcKalTrk->p());
1110 HepLorentzVector ptrk;
1111 ptrk.setPx(mdcKalTrk->px());
1112 ptrk.setPy(mdcKalTrk->py());
1113 ptrk.setPz(mdcKalTrk->pz());
1114 double p3 = ptrk.mag();
1115 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1116 ppip2kpi.push_back(ptrk);
1117 } else {
1118 ipim2kpi.push_back(iGood2kpi[i]);
1119 p3pim2kpi.push_back(mdcKalTrk->p());
1120 HepLorentzVector ptrk;
1121 ptrk.setPx(mdcKalTrk->px());
1122 ptrk.setPy(mdcKalTrk->py());
1123 ptrk.setPz(mdcKalTrk->pz());
1124 double p3 = ptrk.mag();
1125 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1126 ppim2kpi.push_back(ptrk);
1127 }
1128 }
1129 int npip2kpi = ipip2kpi.size();
1130 int npim2kpi = ipim2kpi.size();
1131
1132
1133
1134 ParticleID *pid = ParticleID::instance();
1135
1136if(m_rootput)
1137{
1138 m_cy2kpi=6;
1139}
1140
1141 if((npip2kpi == 2)&&(npim2kpi == 2))
1142 {
1143 //if(nGamnofit>=1)
1144
1145
1146Ncut[4]++;
1147
1148 HepPoint3D xorigin2kpi(0.,0.,0.);
1149 xorigin2kpi[0]=9999.;
1150 xorigin2kpi[1]=9999.;
1151 xorigin2kpi[2]=9999.;
1152 HepSymMatrix xem2kpi(3,0);
1153
1154 Gam4pikp::BubbleSort(p3pip2kpi, ipip2kpi);
1155 Gam4pikp::BubbleSort(p3pim2kpi, ipim2kpi);
1156
1157 Mcut[1]++;
1158 RecMdcKalTrack *pip12kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[0]))->mdcKalTrack();
1159 RecMdcKalTrack *pip22kpiTrk = (*(evtRecTrkCol->begin()+ipip2kpi[1]))->mdcKalTrack();
1160 RecMdcKalTrack *pim12kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[0]))->mdcKalTrack();
1161 RecMdcKalTrack *pim22kpiTrk = (*(evtRecTrkCol->begin()+ipim2kpi[1]))->mdcKalTrack();
1162 double squar2kpi[10]={9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
1163 double mc12kpi,mc22kpi,mc32kpi,mc42kpi;
1164 // double chis2kpi=9999.;
1165 WTrackParameter wcgp12kpi;
1166 WTrackParameter wcgp22kpi;
1167 WTrackParameter wcgm12kpi;
1168 WTrackParameter wcgm22kpi;
1169 int icgp12kpi=-1;
1170 int icgp22kpi=-1;
1171 int icgm12kpi=-1;
1172 int icgm22kpi=-1;
1173 int igam2kpi=-1;
1174 double m2kpi=9999;
1175
1176 int n20=0;
1177 int n30=0;
1178 WTrackParameter wvpip12kpiTrk, wvpim12kpiTrk, wvpip22kpiTrk, wvpim22kpiTrk;
1179 for(int k=0;k<6;k++)
1180 {
1181 if(k==0){mc12kpi=mpi;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpi;}
1182 if(k==1){mc12kpi=mk;mc22kpi=mk;mc32kpi=mk;mc42kpi=mk;}
1183 if(k==2){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpi;mc42kpi=mpro;}
1184 if(k==3){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpro;mc42kpi=mpi;}
1185 if(k==4){mc12kpi=mpi;mc22kpi=mpro;mc32kpi=mpro;mc42kpi=mpi; }
1186 if(k==5){mc12kpi=mpro;mc22kpi=mpi;mc32kpi=mpi;mc42kpi=mpro;}
1187
1188
1189
1190 wvpip12kpiTrk = WTrackParameter(mc12kpi, pip12kpiTrk->getZHelix(), pip12kpiTrk->getZError());
1191 wvpip22kpiTrk = WTrackParameter(mc22kpi, pip22kpiTrk->getZHelix(), pip22kpiTrk->getZError());
1192 wvpim12kpiTrk = WTrackParameter(mc32kpi, pim12kpiTrk->getZHelix(), pim12kpiTrk->getZError());
1193 wvpim22kpiTrk = WTrackParameter(mc42kpi, pim22kpiTrk->getZHelix(), pim22kpiTrk->getZError());
1194 HepPoint3D vx(0., 0., 0.);
1195 HepSymMatrix Evx(3, 0);
1196 double bx = 1E+6;
1197 double by = 1E+6;
1198 double bz = 1E+6;
1199 Evx[0][0] = bx*bx;
1200 Evx[1][1] = by*by;
1201 Evx[2][2] = bz*bz;
1202 VertexParameter vxpar;
1203 vxpar.setVx(vx);
1204 vxpar.setEvx(Evx);
1205 VertexFit* vtxfit = VertexFit::instance();
1206 vtxfit->init();
1207 vtxfit->AddTrack(0, wvpip12kpiTrk);
1208 vtxfit->AddTrack(1, wvpim12kpiTrk);
1209 vtxfit->AddTrack(2, wvpip22kpiTrk);
1210 vtxfit->AddTrack(3, wvpim22kpiTrk);
1211 vtxfit->AddVertex(0, vxpar,0, 1, 2, 3);
1212 if(!vtxfit->Fit(0)) continue;
1213 vtxfit->Swim(0);
1214 WTrackParameter wpip12kpi = vtxfit->wtrk(0);
1215 WTrackParameter wpim12kpi = vtxfit->wtrk(1);
1216 WTrackParameter wpip22kpi = vtxfit->wtrk(2);
1217 WTrackParameter wpim22kpi = vtxfit->wtrk(3);
1218 WTrackParameter wpip12kpiys = vtxfit->wtrk(0);
1219 WTrackParameter wpim12kpiys = vtxfit->wtrk(1);
1220 WTrackParameter wpip22kpiys = vtxfit->wtrk(2);
1221 WTrackParameter wpim22kpiys = vtxfit->wtrk(3);
1222 xorigin2kpi = vtxfit->vx(0);
1223 xem2kpi = vtxfit->Evx(0);
1224 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
1225
1226 HepLorentzVector ecms(0.040547,0,0,3.68632);
1227
1228 double chisq = 9999.;
1229 int ig1 = -1;
1230 for(int i = 0; i < nGam; i++) {
1231 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
1232 kmfit->init();
1233 kmfit->setBeamPosition(xorigin2kpi);
1234 kmfit->setVBeamPosition(xem2kpi);
1235 kmfit->AddTrack(0, wpip12kpi);
1236 kmfit->AddTrack(1, wpim12kpi);
1237 kmfit->AddTrack(2, wpip22kpi);
1238 kmfit->AddTrack(3, wpim22kpi);
1239 kmfit->AddTrack(4, 0.0, g1Trk);
1240 kmfit->AddFourMomentum(0, p4psip);
1241 // if(!kmfit->Fit(0)) continue;
1242 bool oksq = kmfit->Fit();
1243 if(oksq) {
1244 double chi2 = kmfit->chisq();
1245 if(chi2 < chisq) {
1246 chisq = chi2;
1247 squar2kpi[k]=chi2;
1248 ig1 = iGam[i];
1249
1250 }
1251 }
1252 }
1253 if(squar2kpi[k]<200&&squar2kpi[k]<chis2kpi)
1254 {
1255// m_comcs2kpi[k]=squar2kpi[k];
1256 chis2kpi=squar2kpi[k];
1257 if(squar2kpi[k]<20) n20=n20+1;
1258 if(squar2kpi[k]<30) n30=n30+1;
1259
1260 icgp12kpi=ipip2kpi[0];
1261 icgp22kpi=ipip2kpi[1];
1262 icgm12kpi=ipim2kpi[0];
1263 icgm22kpi=ipim2kpi[1];
1264 igam2kpi=ig1;
1265 wcgp12kpi=wpip12kpiys;
1266 wcgp22kpi=wpip22kpiys;
1267 wcgm12kpi=wpim12kpiys;
1268 wcgm22kpi=wpim22kpiys;
1269 cls2kpi=k;
1270
1271 if(k==0){
1272 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1273 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1274 }
1275
1276 if(k==1){
1277 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1278 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1279 }
1280
1281
1282 if(k==2){
1283 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1284 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1285 }
1286
1287 if(k==3){
1288 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[1]);
1289 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[0]);
1290 }
1291
1292 if(k==4){
1293 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
1294 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
1295 }
1296
1297 if(k==5){
1298 itrak2kpi.push_back(ipip2kpi[1]); itrak2kpi.push_back(ipim2kpi[0]);
1299 itrak2kpi.push_back(ipip2kpi[0]); itrak2kpi.push_back(ipim2kpi[1]);
1300 }
1301
1302
1303
1304
1305
1306 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+igam2kpi))->emcShower();
1307 kmfit->init();
1308 kmfit->setBeamPosition(xorigin2kpi);
1309 kmfit->setVBeamPosition(xem2kpi);
1310 kmfit->AddTrack(0, wpip12kpi);
1311 kmfit->AddTrack(1, wpim12kpi);
1312 kmfit->AddTrack(2, wpip22kpi);
1313 kmfit->AddTrack(3, wpim22kpi);
1314 kmfit->AddTrack(4, 0.0, g1Trk);
1315 kmfit->AddFourMomentum(0, p4psip);
1316 bool oksq = kmfit->Fit();
1317 if(oksq){
1318 HepLorentzVector pchic2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3);
1319 HepLorentzVector ppsip2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3) + kmfit->pfit(4);
1320 mchic2kpilst=pchic2kpi.m();
1321 chis4c2kpilst=kmfit->chisq();
1322if(m_rootput)
1323{
1324 m_mchic2kpi = pchic2kpi.m();
1325 m_chis2kpi = kmfit->chisq();
1326 m_mpsip2kpi=ppsip2kpi.m();
1327}
1328
1329
1330
1331 }
1332
1333
1334 }
1335 }
1336
1337
1338 if(chis2kpi<200)
1339 {
1340Ncut[5]++;
1341 if(m_rootput)
1342 {
1343 m_ncy20=n20;
1344 m_ncy30=n30;
1345 m_cla2kpi=cls2kpi;
1346 }
1347 type2kpilst=cls2kpi;
1348
1349 Vp4 p2kpi;
1350 p2kpi.clear();
1351
1352 for (int i = 0; i < 4; i++)
1353 {
1354 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1355 if (!(*itTrk)->isMdcTrackValid()) continue;
1356 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1357 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1358 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1359 double ptrk2kpi;
1360 HepLorentzVector p4trk2kpi;
1361 if(cls2kpi==1)
1362 {
1363 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
1364 ptrk2kpi = mdcKalTrk->p();
1365 p4trk2kpi.setPx(mdcKalTrk->px());
1366 p4trk2kpi.setPy(mdcKalTrk->py());
1367 p4trk2kpi.setPz(mdcKalTrk->pz());
1368 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mk*mk));
1369 p2kpi.push_back(p4trk2kpi);
1370 }
1371
1372 if(cls2kpi==2)
1373 {
1374 if (i<2)
1375 {
1376 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
1377 ptrk2kpi = mdcKalTrk->p();
1378 p4trk2kpi.setPx(mdcKalTrk->px());
1379 p4trk2kpi.setPy(mdcKalTrk->py());
1380 p4trk2kpi.setPz(mdcKalTrk->pz());
1381 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
1382 p2kpi.push_back(p4trk2kpi);
1383
1384 }
1385 else{
1386 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
1387 ptrk2kpi = mdcKalTrk->p();
1388 p4trk2kpi.setPx(mdcKalTrk->px());
1389 p4trk2kpi.setPy(mdcKalTrk->py());
1390 p4trk2kpi.setPz(mdcKalTrk->pz());
1391 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpro*mpro));
1392 p2kpi.push_back(p4trk2kpi);
1393
1394 }
1395
1396 }
1397 if(cls2kpi!=1&&cls2kpi!=2)
1398 {
1399 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
1400 ptrk2kpi = mdcKalTrk->p();
1401 p4trk2kpi.setPx(mdcKalTrk->px());
1402 p4trk2kpi.setPy(mdcKalTrk->py());
1403 p4trk2kpi.setPz(mdcKalTrk->pz());
1404 p4trk2kpi.setE(sqrt(ptrk2kpi*ptrk2kpi+mpi*mpi));
1405 p2kpi.push_back(p4trk2kpi);
1406 }
1407
1408
1409
1410
1411 }
1412
1413 for (int i = 0; i < 4; i++)
1414 {
1415 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1416 if (!(*itTrk)->isMdcTrackValid()) continue;
1417 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1418 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1419 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1420 if ((*itTrk)->isTofTrackValid())
1421 {
1422 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1423 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1424 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
1425 TofHitStatus *status = new TofHitStatus;
1426 status->setStatus((*iter_tof)->status());
1427 if(status->is_cluster()){ dtpr2kpilst[i]=((*iter_tof)->tof()-(*iter_tof)->texpProton());
1428 }
1429 delete status;
1430 }
1431 }
1432 }
1433
1434
1435 /*
1436 HepLorentzVector ecmsb(0.040547,0,0,3.68632);
1437 double chisq = 9999.;
1438 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+igam2kpi))->emcShower();
1439 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
1440 kmfit->init();
1441 kmfit->setBeamPosition(xorigin2kpi);
1442 kmfit->setVBeamPosition(xem2kpi);
1443 kmfit->AddTrack(0, wcgp12kpi);
1444 kmfit->AddTrack(1, wcgp22kpi);
1445 kmfit->AddTrack(2, wcgm12kpi);
1446 kmfit->AddTrack(3, wcgm22kpi);
1447 kmfit->AddTrack(4, 0.0, g1Trk);
1448 kmfit->AddFourMomentum(0, p4psip);
1449 bool oksq = kmfit->Fit();
1450 if(oksq) {
1451 Mcut[3]++;
1452 HepLorentzVector pchic4c2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3);
1453 HepLorentzVector ppsip4c2kpi = kmfit->pfit(0) + kmfit->pfit(1)+kmfit->pfit(2) + kmfit->pfit(3) + kmfit->pfit(4);
1454
1455 m_mchic4c2kpi = pchic4c2kpi.m();
1456// m2kpi=m_mchic4c2kpi;
1457 m_chis4c2kpi = kmfit->chisq();
1458 m_mpsip4c2kpi=ppsip4c2kpi.m();
1459
1460
1461
1462
1463 }
1464
1465*/
1466 Mcut[2]++;
1467if(m_rootput)
1468{
1469 for (int i = 0; i < 4; i++)
1470 {
1471 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + itrak2kpi[i];
1472 if (!(*itTrk)->isMdcTrackValid()) continue;
1473 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
1474 if (!(*itTrk)->isMdcKalTrackValid()) continue;
1475 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
1476 m_ppmdc2kpi[i]=mdcTrk->p();
1477 m_pxmdc2kpi[i]=mdcTrk->px();
1478 m_pymdc2kpi[i]=mdcTrk->py();
1479 m_pzmdc2kpi[i]=mdcTrk->pz();
1480 m_ppkal2kpi[i]=mdcKalTrk->p();
1481 m_charge2kpi[i]=mdcKalTrk->charge();
1482 double ptrk;
1483 ptrk=mdcKalTrk->p();
1484
1485 if (eventHeader->runNumber()<0)
1486 {double mcall=9999;
1487 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
1488 int m_numParticle = 0;
1489 if (!mcParticleCol)
1490 {
1491// std::cout << "Could not retrieve McParticelCol" << std::endl;
1492 return StatusCode::FAILURE;
1493 }
1494 else
1495 {
1496 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
1497 for (; iter_mc != mcParticleCol->end(); iter_mc++)
1498 { double commc;
1499 int pdgcode = (*iter_mc)->particleProperty();
1500 float px,py,pz;
1501 px=(*iter_mc)->initialFourMomentum().x();
1502 py=(*iter_mc)->initialFourMomentum().y();
1503 pz=(*iter_mc)->initialFourMomentum().z();
1504
1505 commc=ptrk*ptrk-px*px-py*py-pz*pz;
1506 if(fabs(commc)<fabs(mcall))
1507 {
1508 mcall=commc;
1509 m_pdg[i]=pdgcode;
1510 m_cbmc[i]=commc;
1511 }
1512
1513 }
1514
1515 }
1516
1517 }
1518
1519 if ((*itTrk)->isTofTrackValid())
1520 { m_bjtofval2kpi[i]=1;
1521 m_cotof12kpi[i]=1;
1522 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1523 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1524 for (;iter_tof != tofTrkCol.end(); iter_tof++ ){
1525 TofHitStatus *status = new TofHitStatus;
1526 status->setStatus((*iter_tof)->status());
1527
1528 if(status->is_cluster()){
1529
1530
1531 m_bjtof2kpi[i]=1;
1532 m_counter2kpi[i] = status->is_counter();
1533 m_barrel2kpi[i] = status->is_barrel();
1534 m_layertof2kpi[i] = status->layer();
1535 m_readout2kpi[i] = status->is_readout();
1536 m_cluster2kpi[i] = status->is_cluster();
1537 m_cotof22kpi[i]=2;
1538 m_beta2kpi[i]= (*iter_tof)->beta();
1539 m_tof2kpi[i] = (*iter_tof)->tof();
1540 m_tofpath2kpi[i] = (*iter_tof)->path();
1541 m_zhit2kpi[i]= (*iter_tof)->zrhit();
1542 m_texe2kpi[i] = (*iter_tof)->texpElectron();
1543 m_texmu2kpi[i] = (*iter_tof)->texpMuon();
1544 m_texpi2kpi[i] = (*iter_tof)->texpPion();
1545 m_texk2kpi[i] = (*iter_tof)->texpKaon();
1546 m_texpr2kpi[i] = (*iter_tof)->texpProton();
1547 m_dte2kpi[i] = m_tof2kpi[i] - m_texe2kpi[i];
1548 m_dtmu2kpi[i] = m_tof2kpi[i] - m_texmu2kpi[i];
1549 m_dtpi2kpi[i] = m_tof2kpi[i] - m_texpi2kpi[i];
1550 m_dtk2kpi[i] = m_tof2kpi[i] - m_texk2kpi[i];
1551 m_dtpr2kpi[i] = m_tof2kpi[i] - m_texpr2kpi[i];
1552 m_tofID2kpi[i]=(*iter_tof)->tofID();
1553 m_clusterID2kpi[i]=(*iter_tof)->tofTrackID();
1554 m_sigmaetof2kpi[i]=(*iter_tof)->sigma(0);
1555 m_sigmamutof2kpi[i]=(*iter_tof)->sigma(1);
1556 m_sigmapitof2kpi[i]=(*iter_tof)->sigma(2);
1557 m_sigmaktof2kpi[i]=(*iter_tof)->sigma(3);
1558 m_sigmaprtof2kpi[i]=(*iter_tof)->sigma(4);
1559 m_t0tof2kpi[i]=(*iter_tof)->t0();
1560 m_errt0tof2kpi[i]=(*iter_tof)->errt0();
1561
1562 }
1563 delete status;
1564 }
1565 }
1566
1567
1568
1569
1570
1571
1572 if ((*itTrk)->isMdcDedxValid())
1573 {
1574 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
1575 m_chie2kpi[i] = dedxTrk->chiE();
1576 m_chimu2kpi[i] = dedxTrk->chiMu();
1577 m_chipi2kpi[i] = dedxTrk->chiPi();
1578 m_chik2kpi[i] = dedxTrk->chiK();
1579 m_chip2kpi[i] = dedxTrk->chiP();
1580 m_ghit2kpi[i] = dedxTrk->numGoodHits();
1581 m_thit2kpi[i] = dedxTrk->numTotalHits();
1582 m_probph2kpi[i] = dedxTrk->probPH();
1583 m_normph2kpi[i] = dedxTrk->normPH();
1584 }
1585
1586 //emc
1587 if ( (*itTrk)->isEmcShowerValid() )
1588 {
1589 RecEmcShower* emcTrk = (*itTrk)->emcShower();
1590 m_bjemc2kpi[i]=1;
1591 m_emc2kpi[i] = emcTrk->energy();
1592 m_evp2kpi[i] = emcTrk->energy()/ptrk;
1593 m_timecg2kpi[i] = emcTrk->time();
1594 }
1595
1596
1597 //muc
1598 if ( (*itTrk)->isMucTrackValid() )
1599 {
1600 RecMucTrack* mucTrk = (*itTrk)->mucTrack();
1601 double dpthp = mucTrk->depth()/25.0;//why?
1602 m_bjmuc2kpi[i]=1;
1603 m_depthmuc2kpi[i]=mucTrk->depth();
1604 m_layermuc2kpi[i] = mucTrk->numLayers();
1605 }
1606
1607 m_cosmdc2kpi[i] = cos(mdcTrk->theta());
1608 m_phimdc2kpi[i] = mdcTrk->phi();
1609
1610 m_pidnum2kpi[i]=(*itTrk)->partId();
1611
1612 if(m_skim4pi)
1613 {
1614 if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1615 {
1616 if(i<4) (*itTrk)->setPartId(3);
1617 }
1618
1619if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1&&dtpr2kpilst[0]<-0.4&&dtpr2kpilst[1]<-0.4&&dtpr2kpilst[2]<-0.4&&dtpr2kpilst[3]<-0.4)
1620 {
1621 if(i<4) (*itTrk)->setPartId(4);
1622 }
1623
1624
1625if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
1626 {
1627 if(i==0||i==1) (*itTrk)->setPartId(3);
1628 if(i==2||i==3) (*itTrk)->setPartId(5);
1629 }
1630
1631 }
1632 //pid
1633 ParticleID *pid = ParticleID::instance();
1634 pid->init();
1635 pid->setMethod(pid->methodProbability());
1636 pid->setChiMinCut(4);
1637 pid->setRecTrack(*itTrk);
1638 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
1639 pid->identify(pid->onlyPion() | pid->onlyKaon() | pid->onlyProton());
1640 pid->identify(pid->onlyElectron() | pid->onlyMuon());
1641 pid->calculate();
1642 if(!(pid->IsPidInfoValid())) continue;
1643 m_costpid2kpi[i] = cos(mdcTrk->theta());
1644
1645
1646 m_probe2kpi[i] = pid->probElectron();
1647 m_probmu2kpi[i] = pid->probMuon();
1648 m_probpi2kpi[i] = pid->probPion();
1649 m_probk2kpi[i] = pid->probKaon();
1650 m_probpr2kpi[i] = pid->probProton();
1651
1652
1653
1654
1655
1656
1657 }
1658}
1659
1660 }
1661// Ncut[10]++;
1662
1663if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1664{
1665 jGam2kpi.push_back(igam2kpi);
1666}
1667
1668 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
1669 if(i>=evtRecTrkCol->size()) break;
1670 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
1671 if(!(*itTrk)->isEmcShowerValid()) continue;
1672if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1673 {
1674 if(i!=igam2kpi) jGam2kpi.push_back((*itTrk)->trackId());
1675 }
1676 else{
1677 jGam2kpi.push_back((*itTrk)->trackId());
1678 }
1679
1680 }
1681
1682int ngam2kpi=jGam2kpi.size();
1683
1684if(m_rootput)
1685{
1686for(int i = 0; i< ngam2kpi; i++) {
1687 if(i>=evtRecTrkCol->size()) break;
1688 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGam2kpi[i];
1689 if(!(*itTrk)->isEmcShowerValid()) continue;
1690 RecEmcShower *emcTrk = (*(evtRecTrkCol->begin()+jGam2kpi[i]))->emcShower();
1691 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
1692 double dthe = 200.;
1693 double dphi = 200.;
1694 double dang = 200.;
1695// double dthert = 200.;
1696// double dphirt = 200.;
1697// double dangrt = 200.;
1698 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
1699 if(j>=evtRecTrkCol->size()) break;
1700 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
1701 if(!(*jtTrk)->isExtTrackValid()) continue;
1702 RecExtTrack *extTrk = (*jtTrk)->extTrack();
1703 if(extTrk->emcVolumeNumber() == -1) continue;
1704 Hep3Vector extpos = extTrk->emcPosition();
1705 double angd = extpos.angle(emcpos);
1706 double thed = extpos.theta() - emcpos.theta();
1707 double phid = extpos.deltaPhi(emcpos);
1708 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
1709 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
1710// if(fabs(thed) < fabs(dthe)) dthe = thed;
1711// if(fabs(phid) < fabs(dphi)) dphi = phid;
1712// if(angd < dang) dang = angd;
1713 if(angd < dang) { dang = angd;
1714 dthe = thed;
1715 dphi = phid;
1716 }
1717
1718 }
1719 if(dang>=200) continue;
1720 double eraw = emcTrk->energy();
1721 dthe = dthe * 180 / (CLHEP::pi);
1722 dphi = dphi * 180 / (CLHEP::pi);
1723 dang = dang * 180 / (CLHEP::pi);
1724// dthert = dthert * 180 / (CLHEP::pi);
1725// dphirt = dphirt * 180 / (CLHEP::pi);
1726// dangrt = dangrt * 180 / (CLHEP::pi);
1727 m_numHits[i]= emcTrk->numHits();
1728 m_secondmoment[i] = emcTrk->secondMoment();
1729 m_latmoment[i] = emcTrk->latMoment();
1730 m_timegm[i] = emcTrk->time();
1731 m_cellId[i]=emcTrk->cellId();
1732 m_module[i]=emcTrk->module();
1733 m_a20Moment[i]=emcTrk->a20Moment();
1734 m_a42Moment[i]=emcTrk->a42Moment();
1735 m_getShowerId[i]=emcTrk->getShowerId();
1736 m_getClusterId[i]=emcTrk->getClusterId();
1737 m_getEAll[i]=emcTrk->getEAll();
1738 m_x[i]= emcTrk->x();
1739 m_y[i]= emcTrk->y();
1740 m_z[i]= emcTrk->z();
1741 m_cosemc[i] = cos(emcTrk->theta());
1742 m_phiemc[i] = emcTrk->phi();
1743 m_energy[i] = emcTrk->energy();
1744 m_eSeed[i] = emcTrk->eSeed();
1745 m_e3x3[i] = emcTrk->e3x3();
1746 m_e5x5[i] = emcTrk->e5x5();
1747 m_dang4c[i] = dang;
1748 m_dthe4c[i] = dthe;
1749 m_dphi4c[i] = dphi;
1750// m_dang4crt[i] = dangrt;
1751// m_dthe4crt[i] = dthert;
1752// m_dphi4crt[i] = dphirt;
1753
1754}
1755}
1756}
1757
1758}
1759
1760
1761
1762
1763
1764 //4pi
1765 if(m_skim4pi)
1766 {
1767
1768if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==0)
1769 {m_subalg1->execute(); // save gam4pi data
1770 Ncut[6]++;
1771 }
1772 }
1773
1774 //4k
1775 if(m_skim4k)
1776 {
1777
1778if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==1&&dtpr2kpilst[0]<-0.4&&dtpr2kpilst[1]<-0.4&&dtpr2kpilst[2]<-0.4&&dtpr2kpilst[3]<-0.4)
1779 {m_subalg2->execute(); // save gam4k data
1780 Ncut[7]++;
1781 }
1782 }
1783
1784 if(m_skim2pi2pr)
1785 {
1786if(mppreclst<3.06&&chis4c2kpilst<20&&((mchic2kpilst>3.39&&mchic2kpilst<3.44)||(mchic2kpilst>3.5&&mchic2kpilst<3.57))&&type2kpilst==2)
1787 {m_subalg3->execute(); // save gam 2(pi p) data
1788 // pi+ pi- with low momentum and p pbar with high momentum.
1789 Ncut[8]++;
1790 }
1791 }
1792
1793
1794
1795//cout<<"chis4c2kpilst="<<chis4c2kpilst<<endl;
1796if(m_rootput)
1797{
1798
1799if((mppreclst<3.06&&chis4c2kpilst<40)||((meelst>3.06&&meelst<3.12)&&(fabs(mppreclst-3.097)<0.01)))
1800 { Ncut[9]++;
1801 m_tuple1->write();
1802
1803 }
1804}
1805
1806 return StatusCode::SUCCESS;
1807 }
1808
1809
1810 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1811 StatusCode Gam4pikp::finalize() {
1812// cout<<"total number: "<<Ncut[0]<<endl;
1813// cout<<"nGood>=4, nCharge==0: "<<Ncut[1]<<endl;
1814// cout<<"nGam>=1: "<<Ncut[2]<<endl;
1815// cout<<"cgp>2 cgm>2: "<<Ncut[3]<<endl;
1816// cout<<"2+ 2-: "<<Ncut[4]<<endl;
1817// cout<<"all cg cycle, chisq<200: "<<Ncut[5]<<endl;
1818// cout<<"4pi ok: "<<Ncut[6]<<endl;
1819// cout<<"4k ok: "<<Ncut[7]<<endl;
1820// cout<<"2pi 2p ok: "<<Ncut[8]<<endl;
1821// cout<<"ntuple write: "<<Ncut[9]<<endl;
1822 MsgStream log(msgSvc(), name());
1823 log << MSG::INFO << "in finalize()" << endmsg;
1824 return StatusCode::SUCCESS;
1825 }
1826 //*************************************************************************
1827 void Gam4pikp::InitVar()
1828 {
1829
1830
1831 m_chis4c2kpi=9999.;
1832 m_chis2kpi=9999.;
1833 m_cla2kpi=9999.;
1834 m_ncy20=9999;
1835 m_ncy30=9999;
1836 m_meeall=9999;
1837 m_mpprecall=9999;
1838 for (int i=0; i<100; i++)
1839 {
1840
1841 m_angjs5[i]=9999;
1842 m_nearjs5[i]=9999;
1843 m_angjs6[i]=9999;
1844 m_nearjs6[i]=9999;
1845 m_ang4pi5[i]=9999;
1846 m_near4pi5[i]=9999;
1847 m_ang4pi6[i]=9999;
1848 m_near4pi6[i]=9999;
1849 m_probe2kpi[i]=-1;
1850 m_probmu2kpi[i]=-1;
1851 m_probpi2kpi[i]=-1;
1852 m_probk2kpi[i]=-1;
1853 m_probpr2kpi[i]=-1;
1854
1855 m_probejs[i]=-1;
1856 m_probmujs[i]=-1;
1857 m_probpijs[i]=-1;
1858 m_probkjs[i]=-1;
1859 m_probprjs[i]=-1;
1860
1861 m_cbmc[i]=9999;
1862 m_bjemcjs[i]=0;
1863 m_bjemc2kpi[i]=0;
1864 m_bjtofjs[i]=0;
1865 m_bjtof2kpi[i]=0;
1866 m_bjmucjs[i]=0;
1867 m_bjmuc2kpi[i]=0;
1868 m_charge2kpi[i]=9999;
1869 m_ghitjs[i] = 9999;
1870 m_thitjs[i] = 9999;
1871 m_probphjs[i] = 99999;
1872 m_normphjs[i] = 99999;
1873
1874 m_ghit2kpi[i] = 9999;
1875 m_thit2kpi[i] = 9999;
1876 m_probph2kpi[i] = 99999;
1877 m_normph2kpi[i] = 99999;
1878
1879
1880
1881 m_counterjs[i] = 9999;
1882 m_barreljs[i] = 9999;
1883 m_layertofjs[i] = 9999;
1884 m_readoutjs[i] = 9999;
1885 m_clusterjs[i] = 9999;
1886
1887 m_counter2kpi[i] = 9999;
1888 m_barrel2kpi[i] = 9999;
1889 m_layertof2kpi[i] = 9999;
1890 m_readout2kpi[i] = 9999;
1891 m_cluster2kpi[i] = 9999;
1892
1893 m_comcs2kpi[i]=9999;
1894 m_dte2kpi[i]=9999;
1895 m_dtmu2kpi[i]=9999;
1896 m_dtpi2kpi[i]=9999;
1897 m_dtk2kpi[i]=9999;
1898 m_dtpr2kpi[i]=9999;
1899 m_sigmaetof2kpi[i]=9999;
1900 m_sigmamutof2kpi[i]=9999;
1901 m_sigmapitof2kpi[i]=9999;
1902 m_sigmaktof2kpi[i]=9999;
1903 m_sigmaprtof2kpi[i]=9999;
1904 m_t0tof2kpi[i]=9999;
1905 m_errt0tof2kpi[i]=9999;
1906
1907 m_sigmaetofjs[i]=9999;
1908 m_sigmamutofjs[i]=9999;
1909 m_sigmapitofjs[i]=9999;
1910 m_sigmaktofjs[i]=9999;
1911 m_sigmaprtofjs[i]=9999;
1912 m_t0tofjs[i]=9999;
1913 m_errt0tofjs[i]=9999;
1914
1915 m_dtejs[i]=9999;
1916 m_dtmujs[i]=9999;
1917 m_dtpijs[i]=9999;
1918 m_dtkjs[i]=9999;
1919 m_dtprjs[i]=9999;
1920
1921 m_chie2kpi[i]=9999;
1922 m_chimu2kpi[i]=9999;
1923 m_chipi2kpi[i]=9999;
1924 m_chik2kpi[i]=9999;
1925 m_chip2kpi[i]=9999;
1926 m_pidnum2kpi[i]=9999;
1927 m_chiejs[i]=9999;
1928 m_chimujs[i]=9999;
1929 m_chipijs[i]=9999;
1930 m_chikjs[i]=9999;
1931 m_chipjs[i]=9999;
1932
1933
1934
1935
1936 }
1937
1938
1939 }
1940
1941 void Gam4pikp::BubbleSort(std::vector<double> &p, std::vector<int> &trkid)
1942 {
1943 if ( (int) p.size() != (int) trkid.size() )
1944 {
1945 std::cout << "the two vector have different length!" << std::endl;
1946 std::cout << "the size of p is " << (int) p.size() << std::endl;
1947 std::cout << "the size of trkid is " << (int) trkid.size() << std::endl;
1948 std::cout << "Please check your code!" << std::endl;
1949 exit(1);
1950 }
1951 unsigned int vsize = p.size();
1952 double ptemp;
1953 int idtemp;
1954 for ( unsigned int i=0; i<vsize-1; i++ )
1955 {
1956 for ( unsigned int j=i+1; j<vsize; j++ )
1957 {
1958 if ( p[i] > p[j] )
1959 {
1960 ptemp = p[i]; idtemp = trkid[i];
1961 p[i] = p[j]; trkid[i] = trkid[j];
1962 p[j] = ptemp; trkid[j] = idtemp;
1963 }
1964 } // for j
1965 } // for i
1966
1967 }
1968
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
Vv
std::vector< VertexParameter > Vv
Definition: Gam4pikp.cxx:56
Vp4
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
xmass
const double xmass[5]
Definition: Gam4pikp.cxx:50
mpro
const double mpro
Definition: Gam4pikp.cxx:49
velc
const double velc
Definition: Gam4pikp.cxx:51
mk
const double mk
Definition: Gam4pikp.cxx:48
Vw
std::vector< WTrackParameter > Vw
Definition: Gam4pikp.cxx:55
Vdouble
std::vector< double > Vdouble
Definition: Gam4pikp.cxx:54
mpi
const double mpi
Definition: Gam4pikp.cxx:47
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
DstEmcShower::cellId
int cellId() const
Definition: DstEmcShower.h:32
DstEmcShower::latMoment
double latMoment() const
Definition: DstEmcShower.h:52
DstEmcShower::a42Moment
double a42Moment() const
Definition: DstEmcShower.h:54
DstEmcShower::eSeed
double eSeed() const
Definition: DstEmcShower.h:47
DstEmcShower::theta
double theta() const
Definition: DstEmcShower.h:38
DstEmcShower::module
int module() const
Definition: DstEmcShower.h:33
DstEmcShower::e3x3
double e3x3() const
Definition: DstEmcShower.h:48
DstEmcShower::phi
double phi() const
Definition: DstEmcShower.h:39
DstEmcShower::secondMoment
double secondMoment() const
Definition: DstEmcShower.h:51
DstEmcShower::x
double x() const
Definition: DstEmcShower.h:35
DstEmcShower::e5x5
double e5x5() const
Definition: DstEmcShower.h:49
DstEmcShower::time
double time() const
Definition: DstEmcShower.h:50
DstEmcShower::z
double z() const
Definition: DstEmcShower.h:37
DstEmcShower::numHits
int numHits() const
Definition: DstEmcShower.h:30
DstEmcShower::a20Moment
double a20Moment() const
Definition: DstEmcShower.h:53
DstEmcShower::energy
double energy() const
Definition: DstEmcShower.h:45
DstEmcShower::y
double y() const
Definition: DstEmcShower.h:36
DstExtTrack::emcPosition
const Hep3Vector emcPosition() const
Definition: DstExtTrack.h:126
DstExtTrack::emcVolumeNumber
const int emcVolumeNumber() const
Definition: DstExtTrack.h:132
DstMdcDedx::probPH
double probPH() const
Definition: DstMdcDedx.h:66
DstMdcDedx::chiE
double chiE() const
Definition: DstMdcDedx.h:59
DstMdcDedx::numTotalHits
int numTotalHits() const
Definition: DstMdcDedx.h:65
DstMdcDedx::numGoodHits
int numGoodHits() const
Definition: DstMdcDedx.h:64
DstMdcDedx::normPH
double normPH() const
Definition: DstMdcDedx.h:67
DstMdcDedx::chiPi
double chiPi() const
Definition: DstMdcDedx.h:61
DstMdcDedx::chiK
double chiK() const
Definition: DstMdcDedx.h:62
DstMdcDedx::chiMu
double chiMu() const
Definition: DstMdcDedx.h:60
DstMdcDedx::chiP
double chiP() const
Definition: DstMdcDedx.h:63
DstMdcKalTrack::setPx
void setPx(const double px, const int pid)
Definition: DstMdcKalTrack.h:112
DstMdcKalTrack::px
const double px() const
Definition: DstMdcKalTrack.h:67
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition: DstMdcKalTrack.h:59
DstMdcKalTrack::pz
const double pz() const
Definition: DstMdcKalTrack.h:69
DstMdcKalTrack::py
const double py() const
Definition: DstMdcKalTrack.h:68
DstMdcKalTrack::p
const double p() const
Definition: DstMdcKalTrack.h:82
DstMdcKalTrack::charge
const int charge() const
Definition: DstMdcKalTrack.h:65
DstMdcKalTrack::pxy
const double pxy() const
Definition: DstMdcKalTrack.h:66
DstMdcTrack::theta
const double theta() const
Definition: DstMdcTrack.h:59
DstMdcTrack::r
const double r() const
Definition: DstMdcTrack.h:64
DstMdcTrack::py
const double py() const
Definition: DstMdcTrack.h:56
DstMdcTrack::err
const HepSymMatrix err() const
Definition: DstMdcTrack.cxx:155
DstMdcTrack::charge
const int charge() const
Definition: DstMdcTrack.h:53
DstMdcTrack::px
const double px() const
Definition: DstMdcTrack.h:55
DstMdcTrack::phi
const double phi() const
Definition: DstMdcTrack.h:60
DstMdcTrack::pz
const double pz() const
Definition: DstMdcTrack.h:57
DstMdcTrack::pxy
const double pxy() const
Definition: DstMdcTrack.h:54
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
DstMdcTrack::z
const double z() const
Definition: DstMdcTrack.h:63
DstMdcTrack::p
const double p() const
Definition: DstMdcTrack.h:58
DstMdcTrack::y
const double y() const
Definition: DstMdcTrack.h:62
DstMdcTrack::x
const double x() const
Definition: DstMdcTrack.h:61
DstMucTrack::depth
double depth() const
Definition: DstMucTrack.h:45
DstMucTrack::numLayers
int numLayers() const
Definition: DstMucTrack.h:42
Gam4pikp::initialize
StatusCode initialize()
Definition: Gam4pikp.cxx:83
Gam4pikp::Gam4pikp
Gam4pikp(const std::string &name, ISvcLocator *pSvcLocator)
Definition: Gam4pikp.cxx:61
Gam4pikp::finalize
StatusCode finalize()
Definition: Gam4pikp.cxx:1811
Gam4pikp::execute
StatusCode execute()
Definition: Gam4pikp.cxx:393
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
KalmanKinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition: KalmanKinematicFit.cxx:353
KalmanKinematicFit::pfit
HepLorentzVector pfit(int n) const
Definition: KalmanKinematicFit.h:154
KalmanKinematicFit::instance
static KalmanKinematicFit * instance()
Definition: KalmanKinematicFit.cxx:20
ParticleIDBase::useTof2
int useTof2() const
Definition: ParticleIDBase.h:66
ParticleIDBase::onlyProton
int onlyProton() const
Definition: ParticleIDBase.h:83
ParticleIDBase::methodProbability
int methodProbability() const
Definition: ParticleIDBase.h:87
ParticleIDBase::useDedx
int useDedx() const
Definition: ParticleIDBase.h:63
ParticleIDBase::onlyMuon
int onlyMuon() const
Definition: ParticleIDBase.h:80
ParticleIDBase::onlyKaon
int onlyKaon() const
Definition: ParticleIDBase.h:82
ParticleIDBase::onlyElectron
int onlyElectron() const
Definition: ParticleIDBase.h:79
ParticleIDBase::onlyPion
int onlyPion() const
Definition: ParticleIDBase.h:81
ParticleIDBase::useTof1
int useTof1() const
Definition: ParticleIDBase.h:65
ParticleIDBase::setChiMinCut
void setChiMinCut(const double chi=4)
Definition: ParticleIDBase.h:39
ParticleIDBase::setRecTrack
void setRecTrack(EvtRecTrack *trk)
Definition: ParticleIDBase.h:37
ParticleID::probKaon
double probKaon() const
Definition: ParticleID.h:124
ParticleID::setMethod
void setMethod(const int method)
Definition: ParticleID.h:94
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::identify
void identify(const int pidcase)
Definition: ParticleID.h:103
ParticleID::probMuon
double probMuon() const
Definition: ParticleID.h:122
ParticleID::probElectron
double probElectron() const
Definition: ParticleID.h:121
ParticleID::usePidSys
void usePidSys(const int pidsys)
Definition: ParticleID.h:97
ParticleID::instance
static ParticleID * instance()
Definition: ParticleID.cxx:22
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::probPion
double probPion() const
Definition: ParticleID.h:123
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition: ParticleID.cxx:94
ParticleID::init
void init()
Definition: ParticleID.cxx:27
ParticleID::probProton
double probProton() const
Definition: ParticleID.h:125
ParticleID::chiDedx
double chiDedx(int n) const
RecEmcShower::getClusterId
RecEmcID getClusterId() const
Definition: RecEmcShower.h:59
RecEmcShower::getShowerId
RecEmcID getShowerId() const
Definition: RecEmcShower.h:55
RecEmcShower::getEAll
RecEmcEnergy getEAll() const
Definition: RecEmcShower.h:92
RecMdcKalTrack::getZHelix
const HepVector & getZHelix() const
Definition: RecMdcKalTrack.h:66
RecMdcKalTrack::getZError
const HepSymMatrix & getZError() const
Definition: RecMdcKalTrack.h:67
TofHitStatus::is_barrel
bool is_barrel() const
Definition: TofHitStatus.h:26
TofHitStatus::layer
unsigned int layer() const
Definition: TofHitStatus.h:28
TofHitStatus::is_cluster
bool is_cluster() const
Definition: TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:126
TofHitStatus::is_counter
bool is_counter() const
Definition: TofHitStatus.h:24
TofHitStatus::is_readout
bool is_readout() const
Definition: TofHitStatus.h:23
TrackPool::setVBeamPosition
void setVBeamPosition(const HepSymMatrix VBeamPosition)
Definition: TrackPool.h:144
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
TrackPool::setBeamPosition
void setBeamPosition(const HepPoint3D BeamPosition)
Definition: TrackPool.h:143
VFHelix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Analysis/VertexFit/VertexFit-00-02-71/VertexFit/Helix.h:186
VFHelix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Analysis/VertexFit/VertexFit-00-02-71/VertexFit/Helix.h:240
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition: VertexFit.h:78
VertexFit::init
void init()
Definition: VertexFit.cxx:29
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition: VertexFit.cxx:89
VertexFit::Evx
HepSymMatrix Evx(int n) const
Definition: VertexFit.h:86
VertexFit::vx
HepPoint3D vx(int n) const
Definition: VertexFit.h:84
VertexFit::instance
static VertexFit * instance()
Definition: VertexFit.cxx:15
VertexFit::Swim
void Swim(int n)
Definition: VertexFit.h:58
VertexFit::Fit
bool Fit()
Definition: VertexFit.cxx:301
VertexParameter::setEvx
void setEvx(const HepSymMatrix &eVx)
Definition: VertexParameter.h:29
VertexParameter::setVx
void setVx(const HepPoint3D &vx)
Definition: VertexParameter.h:54
ecms
const double ecms
Definition: inclkstar.cxx:42
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:111
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:112
vec
dble_vec_t vec[12]
Definition: ranlxd.c:372
pi
const float pi
Definition: vector3.h:133