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