BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Ppjrhopi.cxx
Go to the documentation of this file.
1#include "CLHEP/Vector/LorentzVector.h"
2#include "CLHEP/Vector/ThreeVector.h"
3#include "GaudiKernel/Bootstrap.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/MsgStream.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include <vector>
8
9#include "EventModel/EventHeader.h"
10#include "EventModel/EventModel.h"
11#include "EvtRecEvent/EvtRecEvent.h"
12#include "EvtRecEvent/EvtRecTrack.h"
13#include "VertexDbSvc/IVertexDbSvc.h"
14#include "VertexFit/KinematicFit.h"
15#include "VertexFit/VertexFit.h"
16
17using CLHEP::Hep3Vector;
18using CLHEP::HepLorentzVector;
19#include "CLHEP/Geometry/Point3D.h"
20#ifndef ENABLE_BACKWARDS_COMPATIBILITY
21typedef HepGeom::Point3D<double> HepPoint3D;
22#endif
23#include "Ppjrhopi.h"
24
25// const double twopi = 6.2831853;
26// const double pi = 3.1415927;
27const double mpi = 0.13957;
28const double mk = 0.493677;
29const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
30// const double velc = 29.9792458; tof_path unit in cm.
31const double velc = 299.792458; // tof path unit in mm
32typedef std::vector<int> Vint;
33typedef std::vector<HepLorentzVector> Vp4;
34
35int Ncut0, Ncut1, Ncut2, Ncut3, Ncut4, Ncut5, Ncut6, Ncut7, Ncut8, Ncut9, Ncut10;
36
37/////////////////////////////////////////////////////////////////////////////
38
39DECLARE_COMPONENT( Ppjrhopi )
40Ppjrhopi::Ppjrhopi( const std::string& name, ISvcLocator* pSvcLocator )
41 : Algorithm( name, pSvcLocator ) {
42
43 // Declare the properties
44 declareProperty( "Vr0cut", m_vr0cut = 5.0 );
45 declareProperty( "Vz0cut", m_vz0cut = 20.0 );
46 declareProperty( "Vr1cut", m_vr1cut = 1.0 );
47 declareProperty( "Vz1cut", m_vz1cut = 5.0 );
48 declareProperty( "Vctcut", m_cthcut = 0.93 );
49 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
50 declareProperty( "GammaAngCut", m_gammaAngCut = 20.0 );
51 declareProperty( "Test4C", m_test4C = 1 );
52 declareProperty( "Test5C", m_test5C = 1 );
53 declareProperty( "CheckDedx", m_checkDedx = 1 );
54 declareProperty( "CheckTof", m_checkTof = 1 );
55}
56
57// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
58StatusCode Ppjrhopi::initialize() {
59 MsgStream log( msgSvc(), name() );
60
61 log << MSG::INFO << "in initialize()" << endmsg;
62
63 StatusCode status;
64
65 if ( m_test4C == 1 )
66 {
67 /*
68 NTuplePtr nt4(ntupleSvc(), "FILE1/fit4c");
69 if ( nt4 ) m_tuple4 = nt4;
70 else {
71 m_tuple4 = ntupleSvc()->book ("FILE1/fit4c", CLID_ColumnWiseTuple, "ks N-Tuple
72 example"); if ( m_tuple4 ) { status = m_tuple4->addItem ("run", m_run); status =
73 m_tuple4->addItem ("rec", m_rec); status = m_tuple4->addItem ("nch", m_nch);
74 status = m_tuple4->addItem ("nneu", m_nneu);
75 status = m_tuple4->addItem ("mgdgam", m_gdgam);
76 status = m_tuple4->addItem ("recpp", m_recpp);
77 status = m_tuple4->addItem ("chi2", m_chi1);
78 status = m_tuple4->addItem ("mpi0", m_mpi0);
79 status = m_tuple4->addItem ("mprho0", m_mprho0);
80 status = m_tuple4->addItem ("mprhop", m_mprhop);
81 status = m_tuple4->addItem ("mprhom", m_mprhom);
82 status = m_tuple4->addItem ("mpjjj", m_mpjjj);
83 status = m_tuple4->addItem ("mbepi0", m_bepi0);
84 status = m_tuple4->addItem ("mbe4cjpsi",m_be4cjpsi);
85 status = m_tuple4->addItem ("mp2pi1", m_mp2pi1);
86 status = m_tuple4->addItem ("mf2pi1g1", m_mf2pi1g1);
87 status = m_tuple4->addItem ("mf2pi1g2", m_mf2pi1g2);
88 status = m_tuple4->addItem ("mf2pi1pi0",m_mf2pi1pi0);
89 status = m_tuple4->addItem ("mt2pi2g1", m_mt2pi2g1);
90 status = m_tuple4->addItem ("mt2pi2g2", m_mt2pi2g2);
91 status = m_tuple4->addItem ("mp2pi3", m_mp2pi3);
92 status = m_tuple4->addItem ("mf2pi3g1", m_mf2pi3g1);
93 status = m_tuple4->addItem ("mf2pi3g2", m_mf2pi3g2);
94 status = m_tuple4->addItem ("mf2pi3pi0",m_mf2pi3pi0);
95 status = m_tuple4->addItem ("mp2pi4", m_mp2pi4);
96 status = m_tuple4->addItem ("mf2pi4g1", m_mf2pi4g1);
97 status = m_tuple4->addItem ("mf2pi4g2", m_mf2pi4g2);
98 status = m_tuple4->addItem ("mf2pi4pi0",m_mf2pi4pi0);
99 status = m_tuple4->addItem ("mp4pi", m_mp4pi);
100 status = m_tuple4->addItem ("mppptot", m_mppptot);
101 status = m_tuple4->addItem ("mp4pig1", m_mp4pig1);
102 status = m_tuple4->addItem ("mp4pig2", m_mp4pig2);
103 status = m_tuple4->addItem ("mpx1", m_mpx1);
104 status = m_tuple4->addItem ("mpy1", m_mpy1);
105 status = m_tuple4->addItem ("mpz1", m_mpz1);
106 status = m_tuple4->addItem ("mpe1", m_mpe1);
107 status = m_tuple4->addItem ("mpx2", m_mpx2);
108 status = m_tuple4->addItem ("mpy2", m_mpy2);
109 status = m_tuple4->addItem ("mpz2", m_mpz2);
110 status = m_tuple4->addItem ("mpe2", m_mpe2);
111 status = m_tuple4->addItem ("mpx3", m_mpx3);
112 status = m_tuple4->addItem ("mpy3", m_mpy3);
113 status = m_tuple4->addItem ("mpz3", m_mpz3);
114 status = m_tuple4->addItem ("mpe3", m_mpe3);
115 status = m_tuple4->addItem ("mpx4", m_mpx4);
116 status = m_tuple4->addItem ("mpy4", m_mpy4);
117 status = m_tuple4->addItem ("mpz4", m_mpz4);
118 status = m_tuple4->addItem ("mpe4", m_mpe4);
119 status = m_tuple4->addItem ("mpxg1", m_mpxg1);
120 status = m_tuple4->addItem ("mpyg1", m_mpyg1);
121 status = m_tuple4->addItem ("mpzg1", m_mpzg1);
122 status = m_tuple4->addItem ("mpeg1", m_mpeg1);
123 status = m_tuple4->addItem ("mpxg2", m_mpxg2);
124 status = m_tuple4->addItem ("mpyg2", m_mpyg2);
125 status = m_tuple4->addItem ("mpzg2", m_mpzg2);
126 status = m_tuple4->addItem ("mpeg2", m_mpeg2);
127 status = m_tuple4->addItem ("chikk", m_chikk);
128 status = m_tuple4->addItem ("p1vx", m_p1vx);
129 status = m_tuple4->addItem ("p1vy", m_p1vy);
130 status = m_tuple4->addItem ("p1vz", m_p1vz);
131 status = m_tuple4->addItem ("p1vr", m_p1vr);
132 status = m_tuple4->addItem ("p1vct", m_p1vct);
133 status = m_tuple4->addItem ("m1vx", m_m1vx);
134 status = m_tuple4->addItem ("m1vy", m_m1vy);
135 status = m_tuple4->addItem ("m1vz", m_m1vz);
136 status = m_tuple4->addItem ("m1vr", m_m1vr);
137 status = m_tuple4->addItem ("m1vct", m_m1vct);
138 status = m_tuple4->addItem ("p2vx", m_p2vx);
139 status = m_tuple4->addItem ("p2vy", m_p2vy);
140 status = m_tuple4->addItem ("p2vz", m_p2vz);
141 status = m_tuple4->addItem ("p2vr", m_p2vr);
142 status = m_tuple4->addItem ("p2vct", m_p2vct);
143 status = m_tuple4->addItem ("m2vx", m_m2vx);
144 status = m_tuple4->addItem ("m2vy", m_m2vy);
145 status = m_tuple4->addItem ("m2vz", m_m2vz);
146 status = m_tuple4->addItem ("m2vr", m_m2vr);
147 status = m_tuple4->addItem ("m2vct", m_m2vct);
148 status = m_tuple4->addItem ("mgood", m_good);
149 status = m_tuple4->addItem ("mgam", m_gam);
150 status = m_tuple4->addItem ("mpip", m_pip);
151 status = m_tuple4->addItem ("mpim", m_pim);
152 status = m_tuple4->addItem ("mp1ptot", m_p1ptot);
153 status = m_tuple4->addItem ("memcTp1", m_emcTp1);
154 status = m_tuple4->addItem ("mm1ptot", m_m1ptot);
155 status = m_tuple4->addItem ("memcTm1", m_emcTm1);
156 status = m_tuple4->addItem ("mp2ptot", m_p2ptot);
157 status = m_tuple4->addItem ("memcTp2", m_emcTp2);
158 status = m_tuple4->addItem ("mm2ptot", m_m2ptot);
159 status = m_tuple4->addItem ("memcTm2", m_emcTm2);
160 status = m_tuple4->addItem ("p1pxy", m_p1pxy);
161 status = m_tuple4->addItem ("m1pxy", m_m1pxy);
162 status = m_tuple4->addItem ("p2pxy", m_p2pxy);
163 status = m_tuple4->addItem ("m2pxy", m_m2pxy);
164
165 status = m_tuple4->addItem ("mpidpip", m_pidpip, 0, 10);
166 status = m_tuple4->addIndexedItem ("mipipin" , m_pidpip, m_ipipin);
167 status = m_tuple4->addItem ("mpidpim", m_pidpim, 0, 10);
168 status = m_tuple4->addIndexedItem ("mipimin" , m_pidpim, m_ipimin);
169
170 status = m_tuple4->addItem ("laypip1", m_laypip1);
171 status = m_tuple4->addItem ("laypim1", m_laypim1);
172 status = m_tuple4->addItem ("laypip2", m_laypip2);
173 status = m_tuple4->addItem ("laypim2", m_laypim2);
174 status = m_tuple4->addItem ("mangle", m_angle);
175 status = m_tuple4->addItem ("cosuubr",m_cosuubr );
176 status = m_tuple4->addItem ("cosmupbr", m_cosmupbr);
177 status = m_tuple4->addItem ("cosmumbr", m_cosmumbr);
178 status = m_tuple4->addItem ("phimupbr", m_phimupbr);
179 status = m_tuple4->addItem ("phimumbr", m_phimumbr);
180 status = m_tuple4->addItem ("ngch", m_ngch, 0, 10);
181 status = m_tuple4->addItem ("nggneu", m_nggneu,0, 10);
182 // modifiey by Zhu
183 // status = m_tuple4->addItem ("indx0", indx0, 0, 10);
184 status = m_tuple4->addIndexedItem ("mptrk" , m_ngch, m_ptrk);
185 status = m_tuple4->addIndexedItem ("chie", m_ngch, m_chie);
186 status = m_tuple4->addIndexedItem ("chimu", m_ngch,m_chimu);
187 status = m_tuple4->addIndexedItem ("chipi", m_ngch,m_chipi);
188 status = m_tuple4->addIndexedItem ("chik", m_ngch,m_chik);
189 status = m_tuple4->addIndexedItem ("chip", m_ngch,m_chip);
190 status = m_tuple4->addIndexedItem ("probPH", m_ngch,m_probPH);
191 status = m_tuple4->addIndexedItem ("normPH", m_ngch,m_normPH);
192 status = m_tuple4->addIndexedItem ("ghit", m_ngch,m_ghit);
193 status = m_tuple4->addIndexedItem ("thit", m_ngch,m_thit);
194
195 status = m_tuple4->addIndexedItem ("ptot_etof", m_ngch,m_ptot_etof);
196 status = m_tuple4->addIndexedItem ("cntr_etof", m_ngch,m_cntr_etof);
197 status = m_tuple4->addIndexedItem ("te_etof", m_ngch,m_te_etof);
198 status = m_tuple4->addIndexedItem ("tmu_etof", m_ngch,m_tmu_etof);
199 status = m_tuple4->addIndexedItem ("tpi_etof", m_ngch,m_tpi_etof);
200 status = m_tuple4->addIndexedItem ("tk_etof", m_ngch,m_tk_etof);
201 status = m_tuple4->addIndexedItem ("tp_etof", m_ngch,m_tp_etof);
202 status = m_tuple4->addIndexedItem ("ph_etof", m_ngch,m_ph_etof);
203 status = m_tuple4->addIndexedItem ("rhit_etof", m_ngch,m_rhit_etof);
204 status = m_tuple4->addIndexedItem ("qual_etof", m_ngch,m_qual_etof);
205 status = m_tuple4->addIndexedItem ("ec_toff_e", m_ngch,m_ec_toff_e);
206 status = m_tuple4->addIndexedItem ("ec_toff_mu",m_ngch,m_ec_toff_mu);
207 status = m_tuple4->addIndexedItem ("ec_toff_pi",m_ngch,m_ec_toff_pi);
208 status = m_tuple4->addIndexedItem ("ec_toff_k", m_ngch,m_ec_toff_k);
209 status = m_tuple4->addIndexedItem ("ec_toff_p", m_ngch,m_ec_toff_p);
210 status = m_tuple4->addIndexedItem ("ec_tsig_e", m_ngch,m_ec_tsig_e);
211 status = m_tuple4->addIndexedItem ("ec_tsig_mu",m_ngch,m_ec_tsig_mu);
212 status = m_tuple4->addIndexedItem ("ec_tsig_pi",m_ngch,m_ec_tsig_pi);
213 status = m_tuple4->addIndexedItem ("ec_tsig_k", m_ngch,m_ec_tsig_k);
214 status = m_tuple4->addIndexedItem ("ec_tsig_p", m_ngch,m_ec_tsig_p);
215 status = m_tuple4->addIndexedItem ("ec_tof", m_ngch,m_ec_tof);
216 status = m_tuple4->addIndexedItem ("ptot_btof1",m_ngch,m_ptot_btof1);
217 status = m_tuple4->addIndexedItem ("cntr_btof1",m_ngch,m_cntr_btof1);
218 status = m_tuple4->addIndexedItem ("te_btof1", m_ngch,m_te_btof1);
219 status = m_tuple4->addIndexedItem ("tmu_btof1", m_ngch,m_tmu_btof1);
220 status = m_tuple4->addIndexedItem ("tpi_btof1", m_ngch,m_tpi_btof1);
221 status = m_tuple4->addIndexedItem ("tk_btof1", m_ngch,m_tk_btof1);
222 status = m_tuple4->addIndexedItem ("tp_btof1", m_ngch,m_tp_btof1);
223 status = m_tuple4->addIndexedItem ("ph_btof1", m_ngch,m_ph_btof1);
224 status = m_tuple4->addIndexedItem ("zhit_btof1",m_ngch,m_zhit_btof1);
225 status = m_tuple4->addIndexedItem ("qual_btof1",m_ngch,m_qual_btof1);
226 status = m_tuple4->addIndexedItem ("b1_toff_e", m_ngch,m_b1_toff_e);
227 status = m_tuple4->addIndexedItem ("b1_toff_mu",m_ngch,m_b1_toff_mu);
228 status = m_tuple4->addIndexedItem ("b1_toff_pi",m_ngch,m_b1_toff_pi);
229 status = m_tuple4->addIndexedItem ("b1_toff_k", m_ngch,m_b1_toff_k);
230 status = m_tuple4->addIndexedItem ("b1_toff_p", m_ngch,m_b1_toff_p);
231 status = m_tuple4->addIndexedItem ("b1_tsig_e", m_ngch,m_b1_tsig_e);
232 status = m_tuple4->addIndexedItem ("b1_tsig_mu",m_ngch,m_b1_tsig_mu);
233 status = m_tuple4->addIndexedItem ("b1_tsig_pi",m_ngch,m_b1_tsig_pi);
234 status = m_tuple4->addIndexedItem ("b1_tsig_k", m_ngch,m_b1_tsig_k);
235 status = m_tuple4->addIndexedItem ("b1_tsig_p", m_ngch,m_b1_tsig_p);
236 status = m_tuple4->addIndexedItem ("b1_tof", m_ngch,m_b1_tof);
237
238 status = m_tuple4->addIndexedItem ("mdedx_pid", m_ngch,m_dedx_pid);
239 status = m_tuple4->addIndexedItem ("mtof1_pid", m_ngch,m_tof1_pid);
240 status = m_tuple4->addIndexedItem ("mtof2_pid", m_ngch,m_tof2_pid);
241 status = m_tuple4->addIndexedItem ("mprob_pid", m_ngch,m_prob_pid);
242 status = m_tuple4->addIndexedItem ("mptrk_pid", m_ngch,m_ptrk_pid);
243 status = m_tuple4->addIndexedItem ("mcost_pid", m_ngch,m_cost_pid);
244
245 status = m_tuple4->addIndexedItem ("numHits", m_nggneu,m_numHits); // Total
246 number of hits status = m_tuple4->addIndexedItem ("secondmoment",
247 m_nggneu,m_secondmoment); status = m_tuple4->addIndexedItem ("mx", m_nggneu,m_x); //
248 Shower coordinates and errors status = m_tuple4->addIndexedItem ("my", m_nggneu,m_y);
249 status = m_tuple4->addIndexedItem ("mz", m_nggneu,m_z);
250 status = m_tuple4->addIndexedItem ("cosemc", m_nggneu,m_cosemc); // Shower
251 Counter angles and errors status = m_tuple4->addIndexedItem ("phiemc",
252 m_nggneu,m_phiemc); status = m_tuple4->addIndexedItem ("energy", m_nggneu,m_energy); //
253 Total energy observed in Emc status = m_tuple4->addIndexedItem ("eseed",
254 m_nggneu,m_eSeed); status = m_tuple4->addIndexedItem ("me9", m_nggneu,m_e3x3);
255 status = m_tuple4->addIndexedItem ("me25", m_nggneu,m_e5x5);
256 status = m_tuple4->addIndexedItem ("mlat", m_nggneu,m_lat);
257 status = m_tuple4->addIndexedItem ("ma20", m_nggneu,m_a20);
258 status = m_tuple4->addIndexedItem ("ma42", m_nggneu,m_a42);
259
260
261 }
262 else {
263 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
264 return StatusCode::FAILURE;
265 }
266 }
267 */
268 } // test 4C
269
270 //
271 //--------end of book--------
272 //
273
274 log << MSG::INFO << "successfully return from initialize()" << endmsg;
275 return StatusCode::SUCCESS;
276}
277
278// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
279StatusCode Ppjrhopi::execute() {
280
281 // std::cout << "execute()" << std::endl;
282
283 MsgStream log( msgSvc(), name() );
284 log << MSG::INFO << "in execute()" << endmsg;
285
286 setFilterPassed( false );
287
288 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
289 int runNo = eventHeader->runNumber();
290 int event = eventHeader->eventNumber();
291 log << MSG::DEBUG << "runNo, evtnum = " << runNo << " , " << event << endmsg;
292
293 Ncut0++;
294
295 // FOR 6.4.0 EventModel::Recon--->EventModel::EvtRec IN 6.4.1
296 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
297 log << MSG::INFO << "get event tag OK" << endmsg;
298 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
299 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
300
301 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
302 //
303 // check x0, y0, z0, r0
304 // suggest cut: |z0|<5 && r0<1
305 //
306 if ( evtRecEvent->totalNeutral() > 100 ) { return StatusCode::SUCCESS; }
307
308 Vint iGood, ipip, ipim;
309 iGood.clear();
310 ipip.clear();
311 ipim.clear();
312 Vp4 ppip, ppim;
313 ppip.clear();
314 ppim.clear();
315
316 Hep3Vector xorigin( 0, 0, 0 );
317
318 // if (m_reader.isRunNumberValid(runNo)) {
319 IVertexDbSvc* vtxsvc;
320 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
321 if ( vtxsvc->isVertexValid() )
322 {
323 double* dbv = vtxsvc->PrimaryVertex();
324 double* vv = vtxsvc->SigmaPrimaryVertex();
325 // HepVector dbv = m_reader.PrimaryVertex(runNo);
326 // HepVector vv = m_reader.SigmaPrimaryVertex(runNo);
327 xorigin.setX( dbv[0] );
328 xorigin.setY( dbv[1] );
329 xorigin.setZ( dbv[2] );
330 }
331
332 int nCharge = 0;
333 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
334 {
335 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
336 if ( !( *itTrk )->isMdcTrackValid() ) continue;
337 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
338 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
339
340 double pch = mdcTrk->p();
341 double x0 = mdcTrk->x();
342 double y0 = mdcTrk->y();
343 double z0 = mdcTrk->z();
344 double phi0 = mdcTrk->helix( 1 );
345 double xv = xorigin.x();
346 double yv = xorigin.y();
347 double Rxy = fabs( ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 ) );
348 // 2009//4
349 double m_vx0 = x0;
350 double m_vy0 = y0;
351 double m_vz0 = z0 - xorigin.z();
352 double m_vr0 = Rxy;
353 double m_Vctc = z0 / sqrt( Rxy * Rxy + z0 * z0 );
354 double m_Vct = cos( mdcTrk->theta() );
355
356 // m_tuple1->write();
357 // 2009//4
358 if ( fabs( m_vz0 ) >= m_vz0cut ) continue;
359 if ( m_vr0 >= m_vr0cut ) continue;
360 // if(fabs(m_Vct)>=m_cthcut) continue;
361 iGood.push_back( ( *itTrk )->trackId() );
362 nCharge += mdcTrk->charge();
363 }
364
365 //
366 // Finish Good Charged Track Selection
367 //
368 int nGood = iGood.size();
369
370 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
371 if ( ( nGood != 4 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
372
373 Ncut1++;
374
375 Vint iGam;
376 iGam.clear();
377 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
378 {
379 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
380 if ( !( *itTrk )->isEmcShowerValid() ) continue;
381 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
382 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
383 // find the nearest charged track
384 double dthe = 200.;
385 double dphi = 200.;
386 double dang = 200.;
387 log << MSG::DEBUG << "liuf neu= " << i << endmsg;
388 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
389 {
390 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
391 if ( !( *jtTrk )->isExtTrackValid() ) continue;
392 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
393 if ( extTrk->emcVolumeNumber() == -1 ) continue;
394 Hep3Vector extpos = extTrk->emcPosition();
395 log << MSG::DEBUG << "liuf charge= " << j << endmsg;
396 // double ctht = extpos.cosTheta(emcpos);
397 double angd = extpos.angle( emcpos );
398 double thed = extpos.theta() - emcpos.theta();
399 double phid = extpos.deltaPhi( emcpos );
400 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
401 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
402
403 // if(fabs(thed) < fabs(dthe)) dthe = thed;
404 // if(fabs(phid) < fabs(dphi)) dphi = phid;
405 if ( angd < dang )
406 {
407 dang = angd;
408 dthe = thed;
409 dphi = phid;
410 }
411 }
412 if ( dang >= 200 ) continue;
413 double eraw = emcTrk->energy();
414 dthe = dthe * 180 / ( CLHEP::pi );
415 dphi = dphi * 180 / ( CLHEP::pi );
416 dang = dang * 180 / ( CLHEP::pi );
417 // 2009//4
418 double m_dthe = dthe;
419 double m_dphi = dphi;
420 double m_dang = dang;
421 double m_eraw = eraw;
422 // m_tuple2->write();
423 // 2009//4
424 log << MSG::DEBUG << "eraw dang= " << eraw << " , " << dang << "," << i << endmsg;
425 if ( eraw < m_energyThreshold ) continue;
426 if ( dang < m_gammaAngCut ) continue;
427 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
428 //
429 // good photon cut will be set here
430 //
431 iGam.push_back( ( *itTrk )->trackId() );
432 }
433
434 //
435 // Finish Good Photon Selection
436 //
437 int nGam = iGam.size();
438
439 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
440 << endmsg;
441 if ( nGam < 2 ) { return StatusCode::SUCCESS; }
442
443 Ncut2++;
444
445 //
446 // Assign 4-momentum to each photon
447 //
448
449 Vp4 pGam;
450 pGam.clear();
451 for ( int i = 0; i < nGam; i++ )
452 {
453 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
454 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
455 double eraw = emcTrk->energy();
456 double phi = emcTrk->phi();
457 double the = emcTrk->theta();
458 HepLorentzVector ptrk;
459 ptrk.setPx( eraw * sin( the ) * cos( phi ) );
460 ptrk.setPy( eraw * sin( the ) * sin( phi ) );
461 ptrk.setPz( eraw * cos( the ) );
462 ptrk.setE( eraw );
463
464 // ptrk = ptrk.boost(-0.011,0,0);// boost to cms
465
466 pGam.push_back( ptrk );
467 }
468
469 for ( int i = 0; i < nGood; i++ )
470 { // for rhopi without PID
471 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
472 if ( !( *itTrk )->isMdcTrackValid() ) continue;
473 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
474 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
475 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
476 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
477 if ( mdcKalTrk->charge() > 0 )
478 {
479 ipip.push_back( iGood[i] );
480 HepLorentzVector ptrk;
481 ptrk.setPx( mdcKalTrk->px() );
482 ptrk.setPy( mdcKalTrk->py() );
483 ptrk.setPz( mdcKalTrk->pz() );
484 double p3 = ptrk.mag();
485 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
486 ppip.push_back( ptrk );
487 }
488 else
489 {
490 ipim.push_back( iGood[i] );
491 HepLorentzVector ptrk;
492 ptrk.setPx( mdcKalTrk->px() );
493 ptrk.setPy( mdcKalTrk->py() );
494 ptrk.setPz( mdcKalTrk->pz() );
495 double p3 = ptrk.mag();
496 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
497 ppim.push_back( ptrk );
498 }
499 } // without PID
500
501 int npip = ipip.size();
502 int npim = ipim.size();
503 if ( npip != 2 || npim != 2 ) return StatusCode::SUCCESS;
504 /*
505 log << MSG::DEBUG << "ngood track ID = " << ipip[0] << " , "
506 << ipim[0]<< " , " << ipip[1] << " , " << ipim[1] << endmsg;
507 */
508 Ncut3++;
509
510 //
511 // find the two pi from the primary vetex
512 // ipip[0] && ipim[0] from ppsi
513 // ipip[1] && ipim[1] from jpsi
514 // should change track ID
515 //
516 HepLorentzVector pTot0( 0.011 * 3.6862, 0, 0, 3.6862 );
517 HepLorentzVector pTrec1, pTrec2, pTrec3, pTrec4;
518 HepLorentzVector pTrecf;
519 double m_recjpsi1, m_recjpsi2, m_recjpsi3, m_recjpsi4, m_recppf;
520 double deljp1, deljp2, deljp3, deljp4;
521 pTrec1 = pTot0 - ppip[0] - ppim[0];
522 pTrec2 = pTot0 - ppip[0] - ppim[1];
523 pTrec3 = pTot0 - ppip[1] - ppim[0];
524 pTrec4 = pTot0 - ppip[1] - ppim[1];
525 m_recjpsi1 = pTrec1.m();
526 m_recjpsi2 = pTrec2.m();
527 m_recjpsi3 = pTrec3.m();
528 m_recjpsi4 = pTrec4.m();
529 deljp1 = fabs( m_recjpsi1 - 3.097 );
530 deljp2 = fabs( m_recjpsi2 - 3.097 );
531 deljp3 = fabs( m_recjpsi3 - 3.097 );
532 deljp4 = fabs( m_recjpsi4 - 3.097 );
533
534 int itmp, itmp1, itmp2;
535 HepLorentzVector ptmp, ptmp1, ptmp2;
536
537 pTrecf = pTrec1;
538 m_recppf = pTrec1.m();
539
540 if ( deljp2 < deljp1 && deljp2 < deljp3 && deljp2 < deljp4 )
541 {
542 itmp = ipim[1];
543 ipim[1] = ipim[0];
544 ipim[0] = itmp;
545
546 ptmp = ppim[1];
547 ppim[1] = ppim[0];
548 ppim[0] = ptmp;
549
550 pTrecf = pTrec2;
551 m_recppf = pTrec2.m();
552 }
553
554 if ( deljp3 < deljp1 && deljp3 < deljp2 && deljp3 < deljp4 )
555 {
556 itmp = ipip[1];
557 ipip[1] = ipip[0];
558 ipip[0] = itmp;
559
560 ptmp = ppip[1];
561 ppip[1] = ppip[0];
562 ppip[0] = ptmp;
563
564 pTrecf = pTrec3;
565 m_recppf = pTrec3.m();
566 }
567
568 if ( deljp4 < deljp1 && deljp4 < deljp2 && deljp4 < deljp3 )
569 {
570 itmp1 = ipip[1];
571 ipip[1] = ipip[0];
572 ipip[0] = itmp1;
573 itmp2 = ipim[1];
574 ipim[1] = ipim[0];
575 ipim[0] = itmp2;
576
577 ptmp1 = ppip[1];
578 ppip[1] = ppip[0];
579 ppip[0] = ptmp1;
580 ptmp2 = ppim[1];
581 ppim[1] = ppim[0];
582 ppim[0] = ptmp2;
583
584 pTrecf = pTrec4;
585 m_recppf = pTrec4.m();
586 }
587
588 if ( fabs( m_recppf - 3.097 ) > 0.2 ) return StatusCode::SUCCESS;
589
590 log << MSG::DEBUG << "ngood track ID after jpsi = " << ipip[0] << " , " << ipim[0] << " , "
591 << ipip[1] << " , " << ipim[1] << endmsg;
592 Ncut4++;
593
594 HepLorentzVector ppi2_no1 = ppip[0] + ppim[0];
595 HepLorentzVector ppi2_no2 = ppip[1] + ppim[1];
596 HepLorentzVector ppi2_no3 = ppip[0] + ppim[1];
597 HepLorentzVector ppi2_no4 = ppip[1] + ppim[0];
598 HepLorentzVector p4pi_no = ppi2_no1 + ppi2_no2;
599
600 double emcTg1 = 0.0;
601 double emcTg2 = 0.0;
602 double emcTg3 = 0.0;
603 double emcTg4 = 0.0;
604 double laypi1 = -1.0;
605 double laypi2 = -1.0;
606 double laypi3 = -1.0;
607 double laypi4 = -1.0;
608
609 EvtRecTrackIterator itTrkp1 = evtRecTrkCol->begin() + ipip[0];
610 RecMdcTrack* mdcTrkp1 = ( *itTrkp1 )->mdcTrack();
611 RecMdcKalTrack* mdcKalTrkp1 = ( *itTrkp1 )->mdcKalTrack();
612 RecEmcShower* emcTrkp1 = ( *itTrkp1 )->emcShower();
613 RecMucTrack* mucTrkp1 = ( *itTrkp1 )->mucTrack();
614
615 double phi01 = mdcTrkp1->helix( 1 );
616 double m_p1vx = mdcTrkp1->x();
617 double m_p1vy = mdcTrkp1->y();
618 double m_p1vz = mdcTrkp1->z() - xorigin.z();
619 double m_p1vr = fabs( ( mdcTrkp1->x() - xorigin.x() ) * cos( phi01 ) +
620 ( mdcTrkp1->y() - xorigin.y() ) * sin( phi01 ) );
621 double m_p1vct = cos( mdcTrkp1->theta() );
622 double m_p1ptot = mdcKalTrkp1->p();
623 double m_p1pxy =
624 sqrt( mdcKalTrkp1->px() * mdcKalTrkp1->px() + mdcKalTrkp1->py() * mdcKalTrkp1->py() );
625
626 if ( ( *itTrkp1 )->isEmcShowerValid() ) { emcTg1 = emcTrkp1->energy(); }
627 if ( ( *itTrkp1 )->isMucTrackValid() ) { laypi1 = mucTrkp1->numLayers(); }
628 double m_laypip1 = laypi1;
629
630 EvtRecTrackIterator itTrkm1 = evtRecTrkCol->begin() + ipim[0];
631 RecMdcTrack* mdcTrkm1 = ( *itTrkm1 )->mdcTrack();
632 RecMdcKalTrack* mdcKalTrkm1 = ( *itTrkm1 )->mdcKalTrack();
633 RecEmcShower* emcTrkm1 = ( *itTrkm1 )->emcShower();
634 RecMucTrack* mucTrkm1 = ( *itTrkm1 )->mucTrack();
635
636 double phi02 = mdcTrkm1->helix( 1 );
637 double m_m1vx = mdcTrkm1->x();
638 double m_m1vy = mdcTrkm1->y();
639 double m_m1vz = mdcTrkm1->z() - xorigin.z();
640 double m_m1vr = fabs( ( mdcTrkm1->x() - xorigin.x() ) * cos( phi02 ) +
641 ( mdcTrkm1->y() - xorigin.y() ) * sin( phi02 ) );
642 double m_m1vct = cos( mdcTrkm1->theta() );
643 double m_m1ptot = mdcKalTrkm1->p();
644 double m_m1pxy =
645 sqrt( mdcKalTrkm1->px() * mdcKalTrkm1->px() + mdcKalTrkm1->py() * mdcKalTrkm1->py() );
646
647 if ( ( *itTrkm1 )->isEmcShowerValid() ) { emcTg2 = emcTrkm1->energy(); }
648 if ( ( *itTrkm1 )->isMucTrackValid() ) { laypi2 = mucTrkm1->numLayers(); }
649 double m_laypim1 = laypi2;
650
651 EvtRecTrackIterator itTrkp2 = evtRecTrkCol->begin() + ipip[1];
652 RecMdcTrack* mdcTrkp2 = ( *itTrkp2 )->mdcTrack();
653 RecMdcKalTrack* mdcKalTrkp2 = ( *itTrkp2 )->mdcKalTrack();
654 RecEmcShower* emcTrkp2 = ( *itTrkp2 )->emcShower();
655 RecMucTrack* mucTrkp2 = ( *itTrkp2 )->mucTrack();
656
657 double phi03 = mdcTrkp2->helix( 1 );
658 double m_p2vx = mdcTrkp2->x();
659 double m_p2vy = mdcTrkp2->y();
660 double m_p2vz = mdcTrkp2->z() - xorigin.z();
661 double m_p2vr = fabs( ( mdcTrkp2->x() - xorigin.x() ) * cos( phi03 ) +
662 ( mdcTrkp2->y() - xorigin.y() ) * sin( phi03 ) );
663 double m_p2vct = cos( mdcTrkp2->theta() );
664 double m_p2ptot = mdcKalTrkp2->p();
665 double m_p2pxy =
666 sqrt( mdcKalTrkp2->px() * mdcKalTrkp2->px() + mdcKalTrkp2->py() * mdcKalTrkp2->py() );
667
668 if ( ( *itTrkp2 )->isEmcShowerValid() ) { emcTg3 = emcTrkp2->energy(); }
669 if ( ( *itTrkp2 )->isMucTrackValid() ) { laypi3 = mucTrkp2->numLayers(); }
670 double m_laypip2 = laypi3;
671
672 EvtRecTrackIterator itTrkm2 = evtRecTrkCol->begin() + ipim[1];
673 RecMdcTrack* mdcTrkm2 = ( *itTrkm2 )->mdcTrack();
674 RecMdcKalTrack* mdcKalTrkm2 = ( *itTrkm2 )->mdcKalTrack();
675 RecEmcShower* emcTrkm2 = ( *itTrkm2 )->emcShower();
676 RecMucTrack* mucTrkm2 = ( *itTrkm2 )->mucTrack();
677
678 double phi04 = mdcTrkm2->helix( 1 );
679 double m_m2vx = mdcTrkm2->x();
680 double m_m2vy = mdcTrkm2->y();
681 double m_m2vz = mdcTrkm2->z() - xorigin.z();
682 double m_m2vr = fabs( ( mdcTrkm2->x() - xorigin.x() ) * cos( phi04 ) +
683 ( mdcTrkm2->y() - xorigin.y() ) * sin( phi04 ) );
684 double m_m2vct = cos( mdcTrkm2->theta() );
685 double m_m2ptot = mdcKalTrkm2->p();
686 double m_m2pxy =
687 sqrt( mdcKalTrkm2->px() * mdcKalTrkm2->px() + mdcKalTrkm2->py() * mdcKalTrkm2->py() );
688
689 if ( ( *itTrkm2 )->isEmcShowerValid() ) { emcTg4 = emcTrkm2->energy(); }
690 if ( ( *itTrkm2 )->isMucTrackValid() ) { laypi4 = mucTrkm2->numLayers(); }
691 double m_laypim2 = laypi4;
692
693 double m_emcTp1 = emcTg1;
694 double m_emcTm1 = emcTg2;
695 double m_emcTp2 = emcTg3;
696 double m_emcTm2 = emcTg4;
697
698 if ( fabs( m_p1vz ) >= m_vz1cut ) return StatusCode::SUCCESS;
699 if ( m_p1vr >= m_vr1cut ) return StatusCode::SUCCESS;
700 if ( fabs( m_p1vct ) >= m_cthcut ) return StatusCode::SUCCESS;
701
702 if ( fabs( m_m1vz ) >= m_vz1cut ) return StatusCode::SUCCESS;
703 if ( m_m1vr >= m_vr1cut ) return StatusCode::SUCCESS;
704 if ( fabs( m_m1vct ) >= m_cthcut ) return StatusCode::SUCCESS;
705 Ncut5++;
706
707 HepLorentzVector p4muonp = ppip[1];
708 HepLorentzVector p4muonm = ppim[1];
709 HepLorentzVector p4uu = pTrecf;
710
711 // Lorentz transformation : boost and rotate
712 Hep3Vector p3jpsiUnit = ( p4uu.vect() ).unit();
713 double jBeta = p4uu.beta(); // just same as the P/E
714
715 // std::cout << jBeta << " " << p4uu.beta() << std::endl;
716
717 //
718 // Loop each gamma pair, check ppi0, pTot
719 // and other mass from MDC momentum
720 //
721
722 HepLorentzVector pTot;
723 double minpi0 = 999.0;
724 for ( int i = 0; i < nGam - 1; i++ )
725 {
726 for ( int j = i + 1; j < nGam; j++ )
727 {
728 HepLorentzVector p2g = pGam[i] + pGam[j];
729 pTot = ppip[0] + ppim[0] + ppip[1] + ppim[1];
730 pTot += p2g;
731 if ( fabs( p2g.m() - 0.135 ) < minpi0 )
732 {
733 minpi0 = fabs( p2g.m() - 0.135 );
734 // 2009//4
735 double m_m2gg = p2g.m();
736 // 2009//4
737 HepLorentzVector prho0_no = ppi2_no2;
738 HepLorentzVector prhop_no = ppip[1] + p2g;
739 HepLorentzVector prhom_no = ppim[1] + p2g;
740 HepLorentzVector prho0pi0 = ppi2_no2 + p2g;
741 HepLorentzVector frho1pi0 = ppi2_no1 + p2g;
742 HepLorentzVector frho2pi0 = ppi2_no3 + p2g;
743 HepLorentzVector frho3pi0 = ppi2_no4 + p2g;
744 HepLorentzVector prho0g1 = ppi2_no2 + pGam[i];
745 HepLorentzVector prho0g2 = ppi2_no2 + pGam[j];
746 HepLorentzVector frho1g1 = ppi2_no1 + pGam[i];
747 HepLorentzVector frho1g2 = ppi2_no1 + pGam[j];
748 HepLorentzVector frho2g1 = ppi2_no3 + pGam[i];
749 HepLorentzVector frho2g2 = ppi2_no3 + pGam[j];
750 HepLorentzVector frho3g1 = ppi2_no4 + pGam[i];
751 HepLorentzVector frho3g2 = ppi2_no4 + pGam[j];
752 HepLorentzVector p5pi_no = p4pi_no + p2g;
753
754 // 2009//4
755 double m_prho0_no = prho0_no.m();
756 double m_prhop_no = prhop_no.m();
757 double m_prhom_no = prhom_no.m();
758 double m_prho0pi0 = prho0pi0.m();
759 double m_frho1pi0 = frho1pi0.m();
760 double m_frho2pi0 = frho2pi0.m();
761 double m_frho3pi0 = frho3pi0.m();
762 double m_prho0g1 = prho0g1.m();
763 double m_prho0g2 = prho0g2.m();
764 double m_frho1g1 = frho1g1.m();
765 double m_frho1g2 = frho1g2.m();
766 double m_frho2g1 = frho2g1.m();
767 double m_frho2g2 = frho2g2.m();
768 double m_frho3g1 = frho3g1.m();
769 double m_frho3g2 = frho3g2.m();
770 double m_p4pi_no = p4pi_no.m();
771 double m_p5pi_no = p5pi_no.m();
772 double m_mdcpx1 = ppip[0].px();
773 double m_mdcpy1 = ppip[0].py();
774 double m_mdcpz1 = ppip[0].pz();
775 double m_mdcpe1 = ppip[0].e();
776 double m_mdcpx2 = ppim[0].px();
777 double m_mdcpy2 = ppim[0].py();
778 double m_mdcpz2 = ppim[0].pz();
779 double m_mdcpe2 = ppim[0].e();
780 double m_mdcpx3 = ppip[1].px();
781 double m_mdcpy3 = ppip[1].py();
782 double m_mdcpz3 = ppip[1].pz();
783 double m_mdcpe3 = ppip[1].e();
784 double m_mdcpx4 = ppim[1].px();
785 double m_mdcpy4 = ppim[1].py();
786 double m_mdcpz4 = ppim[1].pz();
787 double m_mdcpe4 = ppim[1].e();
788 double m_mdcpxg1 = pGam[i].px();
789 double m_mdcpyg1 = pGam[i].py();
790 double m_mdcpzg1 = pGam[i].pz();
791 double m_mdcpeg1 = pGam[i].e();
792 double m_mdcpxg2 = pGam[j].px();
793 double m_mdcpyg2 = pGam[j].py();
794 double m_mdcpzg2 = pGam[j].pz();
795 double m_mdcpeg2 = pGam[j].e();
796 double m_etot = pTot.e();
797 double m_mrecjp1 = m_recjpsi1;
798 double m_mrecjp2 = m_recjpsi2;
799 double m_mrecjp3 = m_recjpsi3;
800 double m_mrecjp4 = m_recjpsi4;
801 // m_tuple3 -> write();
802 // 2009//4
803 }
804 }
805 }
806 Ncut6++;
807
808 //
809 // Test vertex fit
810 //
811
812 HepPoint3D vx( 0., 0., 0. );
813 HepSymMatrix Evx( 3, 0 );
814 double bx = 1E+6;
815 double by = 1E+6;
816 double bz = 1E+6;
817 Evx[0][0] = bx * bx;
818 Evx[1][1] = by * by;
819 Evx[2][2] = bz * bz;
820
821 VertexParameter vxpar;
822 vxpar.setVx( vx );
823 vxpar.setEvx( Evx );
824
825 VertexFit* vtxfit = VertexFit::instance();
826 vtxfit->init();
827
828 // assume charged tracks to pi
829
830 RecMdcKalTrack* pipTrk1 = ( *( evtRecTrkCol->begin() + ipip[0] ) )->mdcKalTrack();
831 RecMdcKalTrack* pimTrk1 = ( *( evtRecTrkCol->begin() + ipim[0] ) )->mdcKalTrack();
832 RecMdcKalTrack* pipTrk2 = ( *( evtRecTrkCol->begin() + ipip[1] ) )->mdcKalTrack();
833 RecMdcKalTrack* pimTrk2 = ( *( evtRecTrkCol->begin() + ipim[1] ) )->mdcKalTrack();
834
835 WTrackParameter wvpipTrk1, wvpimTrk1, wvpipTrk2, wvpimTrk2;
836 wvpipTrk1 = WTrackParameter( mpi, pipTrk1->getZHelix(), pipTrk1->getZError() );
837 wvpimTrk1 = WTrackParameter( mpi, pimTrk1->getZHelix(), pimTrk1->getZError() );
838 wvpipTrk2 = WTrackParameter( mpi, pipTrk2->getZHelix(), pipTrk2->getZError() );
839 wvpimTrk2 = WTrackParameter( mpi, pimTrk2->getZHelix(), pimTrk2->getZError() );
840
841 vtxfit->AddTrack( 0, wvpipTrk1 );
842 vtxfit->AddTrack( 1, wvpimTrk1 );
843 vtxfit->AddVertex( 0, vxpar, 0, 1 );
844 if ( !vtxfit->Fit( 0 ) ) return StatusCode::SUCCESS;
845 vtxfit->Swim( 0 );
846
847 Ncut7++;
848
849 WTrackParameter wpip1 = vtxfit->wtrk( 0 );
850 WTrackParameter wpim1 = vtxfit->wtrk( 1 );
851
852 KinematicFit* kmfit = KinematicFit::instance();
853
854 //
855 // Apply Kinematic 4C fit
856 //
857 int igbf1 = -1;
858 int igbf2 = -1;
859 HepLorentzVector pTgam1( 0, 0, 0, 0 );
860 HepLorentzVector pTgam2( 0, 0, 0, 0 );
861
862 if ( m_test4C == 1 )
863 {
864 // double ecms = 3.097;
865 HepLorentzVector ecms( 0.011 * 3.6862, 0, 0, 3.6862 );
866
867 //
868 // kinematic fit to pi pi K K pi0
869 //
870 WTrackParameter wvkipTrk2, wvkimTrk2;
871 wvkipTrk2 = WTrackParameter( mk, pipTrk2->getZHelixK(), pipTrk2->getZErrorK() );
872 wvkimTrk2 = WTrackParameter( mk, pimTrk2->getZHelixK(), pimTrk2->getZErrorK() );
873 double chisq = 9999.;
874 int ig11 = -1;
875 int ig21 = -1;
876 double chikk = 9999.;
877 for ( int i = 0; i < nGam - 1; i++ )
878 {
879 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
880 for ( int j = i + 1; j < nGam; j++ )
881 {
882 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
883 kmfit->init();
884 kmfit->setDynamicerror( 1 );
885 kmfit->AddTrack( 0, wpip1 );
886 kmfit->AddTrack( 1, wpim1 );
887 kmfit->AddTrack( 2, wvkipTrk2 );
888 kmfit->AddTrack( 3, wvkimTrk2 );
889 kmfit->AddTrack( 4, 0.0, g1Trk );
890 kmfit->AddTrack( 5, 0.0, g2Trk );
891 kmfit->AddFourMomentum( 0, ecms );
892 bool oksq = kmfit->Fit();
893 if ( oksq && kmfit->chisq() < chikk ) { chikk = kmfit->chisq(); }
894 }
895 }
896 Ncut8++;
897
898 //
899 // kinematic fit to pi pi pi pi pi0
900 //
901
902 chisq = 9999.;
903 int ig1 = -1;
904 int ig2 = -1;
905 for ( int i = 0; i < nGam - 1; i++ )
906 {
907 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
908 for ( int j = i + 1; j < nGam; j++ )
909 {
910 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
911 kmfit->init();
912 kmfit->setDynamicerror( 1 );
913 kmfit->AddTrack( 0, wpip1 );
914 kmfit->AddTrack( 1, wpim1 );
915 kmfit->AddTrack( 2, wvpipTrk2 );
916 kmfit->AddTrack( 3, wvpimTrk2 );
917 kmfit->AddTrack( 4, 0.0, g1Trk );
918 kmfit->AddTrack( 5, 0.0, g2Trk );
919 kmfit->AddFourMomentum( 0, ecms );
920 bool oksq = kmfit->Fit();
921 if ( oksq )
922 {
923 double chi2 = kmfit->chisq();
924 if ( chi2 < chisq )
925 {
926 chisq = chi2;
927 ig1 = iGam[i];
928 ig2 = iGam[j];
929 igbf1 = iGam[i];
930 igbf2 = iGam[j];
931 pTgam1 = pGam[i];
932 pTgam2 = pGam[j];
933 }
934 }
935 }
936 }
937 // log << MSG::DEBUG << "photon ID from 4c fit to 4pi+pi0 " << ig1 << " , "
938 // << ig2 << endmsg;
939 if ( chisq > 200 ) return StatusCode::SUCCESS;
940 Ncut9++;
941
942 // select charge track and nneu track
943 Vint jGood;
944 jGood.clear();
945 jGood.push_back( ipip[0] );
946 jGood.push_back( ipim[0] );
947 jGood.push_back( ipip[1] );
948 jGood.push_back( ipim[1] );
949
950 Vint jGgam;
951 jGgam.clear();
952 jGgam.push_back( igbf1 );
953 jGgam.push_back( igbf2 );
954
955 double chi1_pp = 9999.0;
956
957 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + ig1 ) )->emcShower();
958 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + ig2 ) )->emcShower();
959 kmfit->init();
960 kmfit->AddTrack( 0, wpip1 );
961 kmfit->AddTrack( 1, wpim1 );
962 kmfit->AddTrack( 2, wvpipTrk2 );
963 kmfit->AddTrack( 3, wvpimTrk2 );
964 kmfit->AddTrack( 4, 0.0, g1Trk );
965 kmfit->AddTrack( 5, 0.0, g2Trk );
966 kmfit->AddFourMomentum( 0, ecms );
967 bool oksq = kmfit->Fit();
968 if ( oksq )
969 {
970 chi1_pp = kmfit->chisq();
971 HepLorentzVector ppi0 = kmfit->pfit( 4 ) + kmfit->pfit( 5 );
972 HepLorentzVector prho0 = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
973 HepLorentzVector prhop = kmfit->pfit( 2 ) + ppi0;
974 HepLorentzVector prhom = kmfit->pfit( 3 ) + ppi0;
975 HepLorentzVector pjjj = prho0 + ppi0;
976
977 HepLorentzVector p2pi1 = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
978 HepLorentzVector f2pi1g1 = p2pi1 + kmfit->pfit( 4 );
979 HepLorentzVector f2pi1g2 = p2pi1 + kmfit->pfit( 5 );
980 HepLorentzVector f2pi1pi0 = p2pi1 + ppi0;
981
982 HepLorentzVector t2pi2g1 = prho0 + kmfit->pfit( 4 );
983 HepLorentzVector t2pi2g2 = prho0 + kmfit->pfit( 5 );
984
985 HepLorentzVector p2pi3 = kmfit->pfit( 0 ) + kmfit->pfit( 3 );
986 HepLorentzVector f2pi3g1 = p2pi3 + kmfit->pfit( 4 );
987 HepLorentzVector f2pi3g2 = p2pi3 + kmfit->pfit( 5 );
988 HepLorentzVector f2pi3pi0 = p2pi3 + ppi0;
989
990 HepLorentzVector p2pi4 = kmfit->pfit( 1 ) + kmfit->pfit( 2 );
991 HepLorentzVector f2pi4g1 = p2pi4 + kmfit->pfit( 4 );
992 HepLorentzVector f2pi4g2 = p2pi4 + kmfit->pfit( 5 );
993 HepLorentzVector f2pi4pi0 = p2pi4 + ppi0;
994
995 HepLorentzVector p4pi = p2pi1 + prho0;
996 HepLorentzVector p4pig1 = p4pi + kmfit->pfit( 4 );
997 HepLorentzVector p4pig2 = p4pi + kmfit->pfit( 5 );
998 HepLorentzVector ppptot = p4pi + ppi0;
999
1000 // add
1001 HepLorentzVector be4cpi0 = pTgam1 + pTgam2;
1002 HepLorentzVector be4c_ppi1 = ppip[0] + ppim[0];
1003 HepLorentzVector be4c_ppi2 = ppip[1] + ppim[1];
1004 HepLorentzVector be4cjp = be4cpi0 + be4c_ppi2;
1005
1006 //**********************************//
1007 // final event selection //
1008 // for pion control sample //
1009 //**********************************//
1010 /*
1011 if(fabs(ppi0.m()-0.135)>0.02) return StatusCode::SUCCESS;
1012 if(fabs(m_recppf-3.097)>0.01) return StatusCode::SUCCESS;
1013 if(fabs(m_emcTp2+m_emcTm2)>2.6) return StatusCode::SUCCESS;
1014 if(chi1_pp>chikk) return StatusCode::SUCCESS;
1015
1016 */
1017 //**********************************//
1018
1019 /*
1020 m_run = eventHeader->runNumber();
1021 m_rec = eventHeader->eventNumber();
1022 m_nch = evtRecEvent->totalCharged();
1023 m_nneu = evtRecEvent->totalNeutral();
1024
1025 m_gdgam=nGam;
1026 m_recpp=m_recppf;
1027 m_ngch = jGood.size();
1028 m_nggneu=jGgam.size();
1029
1030 m_chi1=chi1_pp;
1031 m_bepi0=be4cpi0.m();
1032 m_be4cjpsi = be4cjp.m();
1033 //
1034 m_mpi0 = ppi0.m();
1035 m_mprho0=prho0.m();
1036 m_mprhop=prhop.m();
1037 m_mprhom=prhom.m();
1038 m_mpjjj =pjjj.m();
1039 m_mp2pi1 = p2pi1.m();
1040 m_mf2pi1g1=f2pi1g1.m();
1041 m_mf2pi1g2=f2pi1g2.m();
1042 m_mf2pi1pi0=f2pi1pi0.m();
1043 m_mt2pi2g1=t2pi2g1.m();
1044 m_mt2pi2g2=t2pi2g2.m();
1045 m_mp2pi3 = p2pi3.m();
1046 m_mf2pi3g1=f2pi3g1.m();
1047 m_mf2pi3g2=f2pi3g2.m();
1048 m_mf2pi3pi0=f2pi3pi0.m();
1049 m_mp2pi4 = p2pi4.m();
1050 m_mf2pi4g1=f2pi4g1.m();
1051 m_mf2pi4g2=f2pi4g2.m();
1052 m_mf2pi4pi0=f2pi4pi0.m();
1053 m_mp4pi = p4pi.m();
1054 m_mppptot=ppptot.m();
1055 m_mp4pig1=p4pig1.m();
1056 m_mp4pig2=p4pig2.m();
1057 m_mpx1=kmfit->pfit(0).px();
1058 m_mpy1=kmfit->pfit(0).py();
1059 m_mpz1=kmfit->pfit(0).pz();
1060 m_mpe1=kmfit->pfit(0).e();
1061 m_mpx2=kmfit->pfit(1).px();
1062 m_mpy2=kmfit->pfit(1).py();
1063 m_mpz2=kmfit->pfit(1).pz();
1064 m_mpe2=kmfit->pfit(1).e();
1065 m_mpx3=kmfit->pfit(2).px();
1066 m_mpy3=kmfit->pfit(2).py();
1067 m_mpz3=kmfit->pfit(2).pz();
1068 m_mpe3=kmfit->pfit(2).e();
1069 m_mpx4=kmfit->pfit(3).px();
1070 m_mpy4=kmfit->pfit(3).py();
1071 m_mpz4=kmfit->pfit(3).pz();
1072 m_mpe4=kmfit->pfit(3).e();
1073 m_mpxg1=kmfit->pfit(4).px();
1074 m_mpyg1=kmfit->pfit(4).py();
1075 m_mpzg1=kmfit->pfit(4).pz();
1076 m_mpeg1=kmfit->pfit(4).e();
1077 m_mpxg2=kmfit->pfit(5).px();
1078 m_mpyg2=kmfit->pfit(5).py();
1079 m_mpzg2=kmfit->pfit(5).pz();
1080 m_mpeg2=kmfit->pfit(5).e();
1081 m_good = nGood;
1082 m_chikk=chikk;
1083 m_gam = nGam;
1084 m_pip = npip;
1085 m_pim = npim;
1086
1087 //
1088 // fill information of dedx and tof
1089 //
1090
1091 for(int jk = 0; jk < 2; jk++) {
1092 m_ipipin[jk]=0;
1093 m_ipimin[jk]=0;
1094 }
1095 int ipidpip=0;
1096 int ipidpim=0;
1097
1098 ParticleID *pid = ParticleID::instance();
1099 for(int i = 0; i < npip; i++) {
1100 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + ipip[i];
1101 pid->init();
1102 pid->setMethod(pid->methodProbability());
1103 pid->setChiMinCut(4);
1104 pid->setRecTrack(*itTrk);
1105 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1106 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1107 pid->calculate();
1108 if(!(pid->IsPidInfoValid())) continue;
1109 if(pid->probPion() > pid->probKaon()) {
1110 m_ipipin[ipidpip]=1;
1111 ipidpip++;
1112 }
1113 }
1114
1115 // ParticleID *pid = ParticleID::instance();
1116 for(int j = 0; j < npim; j++) {
1117 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
1118 pid->init();
1119 pid->setMethod(pid->methodProbability());
1120 pid->setChiMinCut(4);
1121 pid->setRecTrack(*jtTrk);
1122 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1123 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1124 pid->calculate();
1125 if(!(pid->IsPidInfoValid())) continue;
1126 if(pid->probPion() > pid->probKaon()) {
1127 m_ipimin[ipidpim]=1;
1128 ipidpim++;
1129 }
1130 }
1131
1132 m_pidpip=ipidpip;
1133 m_pidpim=ipidpim;
1134
1135 //
1136 // check dedx infomation
1137 //
1138
1139 for(int ii = 0; ii < 4; ii++) {
1140 // dedx
1141 m_ptrk[ii] = 9999.0;
1142 m_chie[ii] = 9999.0;
1143 m_chimu[ii] = 9999.0;
1144 m_chipi[ii] = 9999.0;
1145 m_chik[ii] = 9999.0;
1146 m_chip[ii] = 9999.0;
1147 m_ghit[ii] = 9999.0;
1148 m_thit[ii] = 9999.0;
1149 m_probPH[ii] = 9999.0;
1150 m_normPH[ii] = 9999.0;
1151
1152 //endtof
1153 m_cntr_etof[ii] = 9999.0;
1154 m_ptot_etof[ii] = 9999.0;
1155 m_ph_etof[ii] = 9999.0;
1156 m_rhit_etof[ii] = 9999.0;
1157 m_qual_etof[ii] = 9999.0;
1158 m_te_etof[ii] = 9999.0;
1159 m_tmu_etof[ii] = 9999.0;
1160 m_tpi_etof[ii] = 9999.0;
1161 m_tk_etof[ii] = 9999.0;
1162 m_tp_etof[ii] = 9999.0;
1163 m_ec_tof[ii] = 9999.0;
1164 m_ec_toff_e[ii] = 9999.0;
1165 m_ec_toff_mu[ii] = 9999.0;
1166 m_ec_toff_pi[ii] = 9999.0;
1167 m_ec_toff_k[ii] = 9999.0;
1168 m_ec_toff_p[ii] = 9999.0;
1169 m_ec_tsig_e[ii] = 9999.0;
1170 m_ec_tsig_mu[ii] = 9999.0;
1171 m_ec_tsig_pi[ii] = 9999.0;
1172 m_ec_tsig_k[ii] = 9999.0;
1173 m_ec_tsig_p[ii] = 9999.0;
1174
1175 // barrel tof
1176 m_cntr_btof1[ii] = 9999.0;
1177 m_ptot_btof1[ii] = 9999.0;
1178 m_ph_btof1[ii] = 9999.0;
1179 m_zhit_btof1[ii] = 9999.0;
1180 m_qual_btof1[ii] = 9999.0;
1181 m_te_btof1[ii] = 9999.0;
1182 m_tmu_btof1[ii] = 9999.0;
1183 m_tpi_btof1[ii] = 9999.0;
1184 m_tk_btof1[ii] = 9999.0;
1185 m_tp_btof1[ii] = 9999.0;
1186 m_b1_tof[ii] = 9999.0;
1187 m_b1_toff_e[ii] = 9999.0;
1188 m_b1_toff_mu[ii] = 9999.0;
1189 m_b1_toff_pi[ii] = 9999.0;
1190 m_b1_toff_k[ii] = 9999.0;
1191 m_b1_toff_p[ii] = 9999.0;
1192 m_b1_tsig_e[ii] = 9999.0;
1193 m_b1_tsig_mu[ii] = 9999.0;
1194 m_b1_tsig_pi[ii] = 9999.0;
1195 m_b1_tsig_k[ii] = 9999.0;
1196 m_b1_tsig_p[ii] = 9999.0;
1197 //pid
1198 m_dedx_pid[ii] = 9999.0;
1199 m_tof1_pid[ii] = 9999.0;
1200 m_tof2_pid[ii] = 9999.0;
1201 m_prob_pid[ii] = 9999.0;
1202 m_ptrk_pid[ii] = 9999.0;
1203 m_cost_pid[ii] = 9999.0;
1204 }
1205
1206
1207 int indx0=0;
1208 for(int i = 0; i < m_ngch; i++) {
1209 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1210 if(!(*itTrk)->isMdcTrackValid()) continue;
1211 if(!(*itTrk)->isMdcDedxValid())continue;
1212 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1213 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
1214 m_ptrk[indx0] = mdcTrk->p();
1215
1216 m_chie[indx0] = dedxTrk->chiE();
1217 m_chimu[indx0] = dedxTrk->chiMu();
1218 m_chipi[indx0] = dedxTrk->chiPi();
1219 m_chik[indx0] = dedxTrk->chiK();
1220 m_chip[indx0] = dedxTrk->chiP();
1221 m_ghit[indx0] = dedxTrk->numGoodHits();
1222 m_thit[indx0] = dedxTrk->numTotalHits();
1223 m_probPH[indx0] = dedxTrk->probPH();
1224 m_normPH[indx0] = dedxTrk->normPH();
1225 indx0++;
1226 // m_tuple7->write();
1227 }
1228
1229
1230 //
1231 // check TOF infomation
1232 //
1233
1234
1235 int indx1=0;
1236 for(int i = 0; i < m_ngch; i++) {
1237 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1238 if(!(*itTrk)->isMdcTrackValid()) continue;
1239 if(!(*itTrk)->isTofTrackValid()) continue;
1240
1241 RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
1242 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1243
1244 double ptrk = mdcTrk->p();
1245 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1246 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
1247 TofHitStatus *status = new TofHitStatus;
1248 status->setStatus((*iter_tof)->status());
1249 if(!(status->is_barrel())){//endcap
1250 if( !(status->is_counter()) ) continue; // ?
1251 if( status->layer()!=1 ) continue;//layer1
1252 double path=(*iter_tof)->path(); // ?
1253 double tof = (*iter_tof)->tof();
1254 double ph = (*iter_tof)->ph();
1255 double rhit = (*iter_tof)->zrhit();
1256 double qual = 0.0 + (*iter_tof)->quality();
1257 double cntr = 0.0 + (*iter_tof)->tofID();
1258 double texp[5];
1259 double tsig[5];
1260 for(int j = 0; j < 5; j++) {//0 e, 1 mu, 2 pi, 3 K, 4 p
1261 texp[j] = (*iter_tof)->texp(j);
1262 // tsig[j] = (*iter_tof)->sigma(j);
1263 // toffset[j] = (*iter_tof)->offset(j);
1264 }
1265 m_cntr_etof[indx1] = cntr;
1266 m_ptot_etof[indx1] = ptrk;
1267 m_ph_etof[indx1] = ph;
1268 m_rhit_etof[indx1] = rhit;
1269 m_qual_etof[indx1] = qual;
1270 m_te_etof[indx1] = tof - texp[0];
1271 m_tmu_etof[indx1] = tof - texp[1];
1272 m_tpi_etof[indx1] = tof - texp[2];
1273 m_tk_etof[indx1] = tof - texp[3];
1274 m_tp_etof[indx1] = tof - texp[4];
1275
1276 m_ec_tof[indx1] = tof;
1277
1278 m_ec_toff_e[indx1] = (*iter_tof)->toffset(0);
1279 m_ec_toff_mu[indx1] = (*iter_tof)->toffset(1);
1280 m_ec_toff_pi[indx1] = (*iter_tof)->toffset(2);
1281 m_ec_toff_k[indx1] = (*iter_tof)->toffset(3);
1282 m_ec_toff_p[indx1] = (*iter_tof)->toffset(4);
1283
1284 m_ec_tsig_e[indx1] = (*iter_tof)->sigma(0);
1285 m_ec_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1286 m_ec_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1287 m_ec_tsig_k[indx1] = (*iter_tof)->sigma(3);
1288 m_ec_tsig_p[indx1] = (*iter_tof)->sigma(4);
1289
1290 // m_tuple8->write();
1291 }
1292 else {//barrel
1293 if( !(status->is_cluster()) ) continue; // ?
1294 double path=(*iter_tof)->path(); // ?
1295 double tof = (*iter_tof)->tof();
1296 double ph = (*iter_tof)->ph();
1297 double rhit = (*iter_tof)->zrhit();
1298 double qual = 0.0 + (*iter_tof)->quality();
1299 double cntr = 0.0 + (*iter_tof)->tofID();
1300 double texp[5];
1301 for(int j = 0; j < 5; j++) {
1302 texp[j] = (*iter_tof)->texp(j);
1303 }
1304 m_cntr_btof1[indx1] = cntr;
1305 m_ptot_btof1[indx1] = ptrk;
1306 m_ph_btof1[indx1] = ph;
1307 m_zhit_btof1[indx1] = rhit;
1308 m_qual_btof1[indx1] = qual;
1309 m_te_btof1[indx1] = tof - texp[0];
1310 m_tmu_btof1[indx1] = tof - texp[1];
1311 m_tpi_btof1[indx1] = tof - texp[2];
1312 m_tk_btof1[indx1] = tof - texp[3];
1313 m_tp_btof1[indx1] = tof - texp[4];
1314
1315 m_b1_tof[indx1] = tof;
1316
1317 m_b1_toff_e[indx1] = (*iter_tof)->toffset(0);
1318 m_b1_toff_mu[indx1] = (*iter_tof)->toffset(1);
1319 m_b1_toff_pi[indx1] = (*iter_tof)->toffset(2);
1320 m_b1_toff_k[indx1] = (*iter_tof)->toffset(3);
1321 m_b1_toff_p[indx1] = (*iter_tof)->toffset(4);
1322
1323 m_b1_tsig_e[indx1] = (*iter_tof)->sigma(0);
1324 m_b1_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1325 m_b1_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1326 m_b1_tsig_k[indx1] = (*iter_tof)->sigma(3);
1327 m_b1_tsig_p[indx1] = (*iter_tof)->sigma(4);
1328
1329 // m_tuple9->write();
1330 }
1331 delete status;
1332 }
1333 indx1++;
1334 } // loop all charged track
1335
1336 //
1337 // Assign 4-momentum to each charged track
1338 //
1339 int indx2=0;
1340 // ParticleID *pid = ParticleID::instance();
1341 for(int i = 0; i < m_ngch; i++) {
1342 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1343 // if(pid) delete pid;
1344 pid->init();
1345 pid->setMethod(pid->methodProbability());
1346 // pid->setMethod(pid->methodLikelihood()); //for Likelihood Method
1347
1348 pid->setChiMinCut(4);
1349 pid->setRecTrack(*itTrk);
1350 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID
1351 sub-system pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1352 // pid->identify(pid->onlyPion());
1353 // pid->identify(pid->onlyKaon());
1354 pid->calculate();
1355 if(!(pid->IsPidInfoValid())) continue;
1356 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1357
1358 m_dedx_pid[indx2] = pid->chiDedx(2);
1359 m_tof1_pid[indx2] = pid->chiTof1(2);
1360 m_tof2_pid[indx2] = pid->chiTof2(2);
1361 m_prob_pid[indx2] = pid->probPion();
1362
1363 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
1364 // if(pid->probPion() < 0.001) continue;
1365 // if(pid->pdf(2)<pid->pdf(3)) continue; // for Likelihood Method(0=electron 1=muon
1366 2=pion 3=kaon 4=proton)
1367
1368 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use
1369 RecMdcKalTrack substitute RecMdcTrack RecMdcKalTrack::setPidType
1370 (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default
1371 setting is pion
1372
1373 m_ptrk_pid[indx2] = mdcKalTrk->p();
1374 m_cost_pid[indx2] = cos(mdcTrk->theta());
1375 }
1376
1377
1378 int iphoton = 0;
1379 for (int i=0; i<m_nggneu; i++)
1380 {
1381 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGgam[i];
1382 if (!(*itTrk)->isEmcShowerValid()) continue;
1383 RecEmcShower *emcTrk = (*itTrk)->emcShower();
1384 m_numHits[iphoton] = emcTrk->numHits();
1385 m_secondmoment[iphoton] = emcTrk->secondMoment();
1386 m_x[iphoton] = emcTrk->x();
1387 m_y[iphoton]= emcTrk->y();
1388 m_z[iphoton]= emcTrk->z();
1389 m_cosemc[iphoton] = cos(emcTrk->theta());
1390 m_phiemc[iphoton] = emcTrk->phi();
1391 m_energy[iphoton] = emcTrk->energy();
1392 m_eSeed[iphoton] = emcTrk->eSeed();
1393 m_e3x3[iphoton] = emcTrk->e3x3();
1394 m_e5x5[iphoton] = emcTrk->e5x5();
1395 m_lat[iphoton] = emcTrk->latMoment();
1396 m_a20[iphoton] = emcTrk->a20Moment();
1397 m_a42[iphoton] = emcTrk->a42Moment();
1398 iphoton++;
1399 }
1400 // m_tuple4->write();
1401 */
1402 Ncut10++;
1403 // log << MSG::DEBUG << "chisquare from 4c fit to 4pi+pi0 " << m_chi1 << endmsg;
1404 }
1405 }
1406
1407 //
1408 // Apply Kinematic 5C Fit
1409 //
1410
1411 // find the best combination over all possible pi+ pi- gamma gamma pair
1412
1413 setFilterPassed( true );
1414
1415 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setPartId( 2 );
1416 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setPartId( 2 );
1417 ( *( evtRecTrkCol->begin() + ipip[1] ) )->setPartId( 2 );
1418 ( *( evtRecTrkCol->begin() + ipim[1] ) )->setPartId( 2 );
1419
1420 return StatusCode::SUCCESS;
1421}
1422
1423// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1424StatusCode Ppjrhopi::finalize() {
1425 cout << "total number: " << Ncut0 << endl;
1426 cout << "nGood==4, nCharge==0: " << Ncut1 << endl;
1427 cout << "nGam>=2: " << Ncut2 << endl;
1428 cout << "Pass no Pid: " << Ncut3 << endl;
1429 cout << "ChangeID recfrom psp: " << Ncut4 << endl;
1430 cout << "vetex position: " << Ncut5 << endl;
1431 cout << "Mass from MDC: " << Ncut6 << endl;
1432 cout << "primary vetex fit: " << Ncut7 << endl;
1433 cout << "Pass 4C for ppkkp0: " << Ncut8 << endl;
1434 cout << "Pass 4C for 4pi+pi0: " << Ncut9 << endl;
1435 cout << "Pass 4C <200: " << Ncut10 << endl;
1436 MsgStream log( msgSvc(), name() );
1437 log << MSG::INFO << "in finalize()" << endmsg;
1438 return StatusCode::SUCCESS;
1439}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
runNo
int runNo
Definition DQA_TO_DB.cxx:13
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi
double mpi
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
xmass
const double xmass[5]
Definition Gam4pikp.cxx:35
velc
const double velc
Definition Gam4pikp.cxx:36
mk
const double mk
Definition Gam4pikp.cxx:33
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:37
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
Ncut2
int Ncut2
Definition Ppjrhopi.cxx:35
Ncut1
int Ncut1
Definition Ppjrhopi.cxx:35
Ncut0
int Ncut0
Definition Ppjrhopi.cxx:35
Ncut5
int Ncut5
Definition Ppjrhopi.cxx:35
Ncut10
int Ncut10
Definition Ppjrhopi.cxx:35
Ncut9
int Ncut9
Definition Ppjrhopi.cxx:35
Ncut3
int Ncut3
Definition Ppjrhopi.cxx:35
Ncut8
int Ncut8
Definition Ppjrhopi.cxx:35
Ncut6
int Ncut6
Definition Ppjrhopi.cxx:35
Ncut7
int Ncut7
Definition Ppjrhopi.cxx:35
Ncut4
int Ncut4
Definition Ppjrhopi.cxx:35
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:18
KinematicFit::setDynamicerror
void setDynamicerror(const bool dynamicerror=1)
Definition Analysis/VertexFit/include/VertexFit/KinematicFit.h:130
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
Ppjrhopi::finalize
StatusCode finalize()
Definition Ppjrhopi.cxx:1424
Ppjrhopi::initialize
StatusCode initialize()
Definition Ppjrhopi.cxx:58
Ppjrhopi::Ppjrhopi
Ppjrhopi(const std::string &name, ISvcLocator *pSvcLocator)
Definition Ppjrhopi.cxx:40
Ppjrhopi::execute
StatusCode execute()
Definition Ppjrhopi.cxx:279
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:91
VertexFit::init
void init()
Definition VertexFit.cxx:27
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288
ecms
const double ecms
Definition inclkstar.cxx:26