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