BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
DQARhopi.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#include "McTruth/McParticle.h"
19
20#include "TH1F.h"
21#include "TMath.h"
22#include "GaudiKernel/INTupleSvc.h"
23#include "GaudiKernel/NTuple.h"
24#include "GaudiKernel/ITHistSvc.h"
25
26#include "GaudiKernel/Bootstrap.h"
27#include "GaudiKernel/IHistogramSvc.h"
28#include "CLHEP/Vector/ThreeVector.h"
29#include "CLHEP/Vector/LorentzVector.h"
30#include "CLHEP/Vector/TwoVector.h"
31using CLHEP::Hep3Vector;
32using CLHEP::Hep2Vector;
33using CLHEP::HepLorentzVector;
34#include "CLHEP/Geometry/Point3D.h"
35#ifndef ENABLE_BACKWARDS_COMPATIBILITY
36 typedef HepGeom::Point3D<double> HepPoint3D;
37#endif
38#include "DQARhopiAlg/DQARhopi.h"
39
40#include "VertexFit/KinematicFit.h"
41#include "VertexFit/VertexFit.h"
42#include "ParticleID/ParticleID.h"
43
44#include <vector>
45using namespace Event;
46//const double twopi = 6.2831853;
47//const double pi = 3.1415927;
48const double mpi = 0.13957;
49const double mk = 0.493677;
50const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
51//const double velc = 29.9792458; tof_path unit in cm.
52const double velc = 299.792458; // tof path unit in mm
53typedef std::vector<int> Vint;
54typedef std::vector<HepLorentzVector> Vp4;
55
56const HepLorentzVector ecms(0.034,0,0,3.097);
57const Hep3Vector u_cms = -ecms.boostVector();
58
59int Ncut0,Ncut1,Ncut2,Ncut3,Ncut4,Ncut5,Ncut6,Ncut7,Ncut8,Ncut9,Ncut10;
60
61/////////////////////////////////////////////////////////////////////////////
62
63DQARhopi::DQARhopi(const std::string& name, ISvcLocator* pSvcLocator) :
64 Algorithm(name, pSvcLocator) {
65
66 //Declare the properties
67 declareProperty("Vr0cut", m_vr0cut=1.0);
68 declareProperty("Vz0cut", m_vz0cut=5.0);
69 declareProperty("Vctcut", m_cthcut=0.93);
70 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
71 declareProperty("GammaAngCut", m_gammaAngCut=25.0);
72 declareProperty("Test4C", m_test4C = 1);
73 declareProperty("Test5C", m_test5C = 1);
74 declareProperty("CheckDedx", m_checkDedx = 1);
75 declareProperty("CheckTof", m_checkTof = 1);
76}
77
78// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
79StatusCode DQARhopi::initialize(){
80 MsgStream log(msgSvc(), name());
81
82 log << MSG::INFO << "in initialize()" << endmsg;
83
84 StatusCode status;
85
86 NTuplePtr nt4(ntupleSvc(), "DQAFILE/Rhopi");
87 if ( nt4 ) m_tuple4 = nt4;
88 else {
89 m_tuple4 = ntupleSvc()->book ("DQAFILE/Rhopi", CLID_ColumnWiseTuple, "ks N-Tuple example");
90 if ( m_tuple4 ) {
91 status = m_tuple4->addItem ("run", m_run);
92 status = m_tuple4->addItem ("rec", m_rec);
93 status = m_tuple4->addItem ("nch", m_nch);
94 status = m_tuple4->addItem ("nneu", m_nneu);
95 status = m_tuple4->addItem ("chi1", m_chi1);
96 status = m_tuple4->addItem ("mpi0", m_mpi0);
97 status = m_tuple4->addItem ("mprho0", m_prho0);
98 status = m_tuple4->addItem ("mprhop", m_prhop);
99 status = m_tuple4->addItem ("mprhom", m_prhom);
100 status = m_tuple4->addItem ("mgood", m_good);
101 status = m_tuple4->addItem ("mgam", m_gam);
102 status = m_tuple4->addItem ("mpip", m_pip);
103 status = m_tuple4->addItem ("mpim", m_pim);
104 status = m_tuple4->addItem ("m2gam", m_2gam);
105 status = m_tuple4->addItem ("ngch", m_ngch, 0, 10);
106 status = m_tuple4->addItem ("nggneu", m_nggneu,0, 10);
107 status = m_tuple4->addItem ("moutpi0",m_outpi0);
108 status = m_tuple4->addItem ("cosang", m_cosang);
109 status = m_tuple4->addItem ("moutpip",m_outpip);
110 status = m_tuple4->addItem ("moutpim",m_outpim);
111 status = m_tuple4->addItem ("menpip", m_enpip);
112 status = m_tuple4->addItem ("mdcpip", m_dcpip );
113 status = m_tuple4->addItem ("menpim", m_enpim );
114 status = m_tuple4->addItem ("mdcpim", m_dcpim );
115 status = m_tuple4->addItem ("mpipf", m_pipf);
116 status = m_tuple4->addItem ("mpimf", m_pimf);
117 status = m_tuple4->addItem ("mpi0f", m_pi0f);
118
119 status = m_tuple4->addItem ("mpmax", m_pmax);
120 status = m_tuple4->addItem ("mppx", m_ppx);
121 status = m_tuple4->addItem ("mppy", m_ppy);
122 status = m_tuple4->addItem ("mppz", m_ppz);
123 status = m_tuple4->addItem ("mcosthep",m_costhep);
124 status = m_tuple4->addItem ("mppxkal", m_ppxkal);
125 status = m_tuple4->addItem ("mppykal", m_ppykal);
126 status = m_tuple4->addItem ("mppzkal", m_ppzkal);
127 status = m_tuple4->addItem ("mmpx", m_mpx);
128 status = m_tuple4->addItem ("mmpy", m_mpy);
129 status = m_tuple4->addItem ("mmpz", m_mpz);
130 status = m_tuple4->addItem ("mcosthem",m_costhem);
131 status = m_tuple4->addItem ("mmpxkal", m_mpxkal);
132 status = m_tuple4->addItem ("mmpykal", m_mpykal);
133 status = m_tuple4->addItem ("mmpzkal", m_mpzkal);
134
135 status = m_tuple4->addItem ("mvxpin", m_vxpin);
136 status = m_tuple4->addItem ("mvypin", m_vypin);
137 status = m_tuple4->addItem ("mvzpin", m_vzpin);
138 status = m_tuple4->addItem ("mvrpin", m_vrpin);
139 status = m_tuple4->addItem ("mcosthepin",m_costhepin);
140 status = m_tuple4->addItem ("mvxmin", m_vxmin);
141 status = m_tuple4->addItem ("mvymin", m_vymin);
142 status = m_tuple4->addItem ("mvzmin", m_vzmin);
143 status = m_tuple4->addItem ("mvrmin", m_vrmin);
144 status = m_tuple4->addItem ("mcosthemin",m_costhemin);
145 status = m_tuple4->addItem ("mvxp", m_vxp);
146 status = m_tuple4->addItem ("mvyp", m_vyp);
147 status = m_tuple4->addItem ("mvzp", m_vzp);
148 status = m_tuple4->addItem ("mvrp", m_vrp);
149 status = m_tuple4->addItem ("mvxm", m_vxm);
150 status = m_tuple4->addItem ("mvym", m_vym);
151 status = m_tuple4->addItem ("mvzm", m_vzm);
152 status = m_tuple4->addItem ("mvrm", m_vrm);
153 status = m_tuple4->addItem ("nangecc",m_nangecc,0,10);
154 status = m_tuple4->addIndexedItem ("mdthec", m_nangecc, m_dthec);
155 status = m_tuple4->addIndexedItem ("mdphic", m_nangecc, m_dphic);
156 status = m_tuple4->addIndexedItem ("mdangc", m_nangecc, m_dangc);
157 status = m_tuple4->addIndexedItem ("mspippim", m_nangecc, m_mspippim);
158
159 status = m_tuple4->addItem ("angsg", dangsg);
160 status = m_tuple4->addItem ("thesg", dthesg);
161 status = m_tuple4->addItem ("phisg", dphisg);
162 status = m_tuple4->addItem ("cosgth1", cosgth1);
163 status = m_tuple4->addItem ("cosgth2", cosgth2);
164
165 status = m_tuple4->addItem ("mchi5", m_chi5);
166 status = m_tuple4->addItem ("mkpi0", m_kpi0);
167 status = m_tuple4->addItem ("mkpkm", m_kpkm);
168 status = m_tuple4->addItem ("mkpp0", m_kpp0);
169 status = m_tuple4->addItem ("mkmp0", m_kmp0);
170 status = m_tuple4->addItem ("mpgam2pi1",m_pgam2pi1);
171 status = m_tuple4->addItem ("mpgam2pi2",m_pgam2pi2);
172 status = m_tuple4->addItem ("cosva1", cosva1);
173 status = m_tuple4->addItem ("cosva2", cosva2);
174 status = m_tuple4->addItem ("laypi1", m_laypi1);
175 status = m_tuple4->addItem ("hit1", m_hit1);
176 status = m_tuple4->addItem ("laypi2", m_laypi2);
177 status = m_tuple4->addItem ("hit2", m_hit2);
178 status = m_tuple4->addItem ("anglepm", m_anglepm);
179
180 status = m_tuple4->addIndexedItem ("mptrk", m_ngch, m_ptrk);
181 status = m_tuple4->addIndexedItem ("chie", m_ngch, m_chie);
182 status = m_tuple4->addIndexedItem ("chimu", m_ngch,m_chimu);
183 status = m_tuple4->addIndexedItem ("chipi", m_ngch,m_chipi);
184 status = m_tuple4->addIndexedItem ("chik", m_ngch,m_chik);
185 status = m_tuple4->addIndexedItem ("chip", m_ngch,m_chip);
186 status = m_tuple4->addIndexedItem ("probPH", m_ngch,m_probPH);
187 status = m_tuple4->addIndexedItem ("normPH", m_ngch,m_normPH);
188 status = m_tuple4->addIndexedItem ("ghit", m_ngch,m_ghit);
189 status = m_tuple4->addIndexedItem ("thit", m_ngch,m_thit);
190
191 status = m_tuple4->addIndexedItem ("ptot_etof", m_ngch,m_ptot_etof);
192 status = m_tuple4->addIndexedItem ("cntr_etof", m_ngch,m_cntr_etof);
193 status = m_tuple4->addIndexedItem ("te_etof", m_ngch,m_te_etof);
194 status = m_tuple4->addIndexedItem ("tmu_etof", m_ngch,m_tmu_etof);
195 status = m_tuple4->addIndexedItem ("tpi_etof", m_ngch,m_tpi_etof);
196 status = m_tuple4->addIndexedItem ("tk_etof", m_ngch,m_tk_etof);
197 status = m_tuple4->addIndexedItem ("tp_etof", m_ngch,m_tp_etof);
198 status = m_tuple4->addIndexedItem ("ph_etof", m_ngch,m_ph_etof);
199 status = m_tuple4->addIndexedItem ("rhit_etof", m_ngch,m_rhit_etof);
200 status = m_tuple4->addIndexedItem ("qual_etof", m_ngch,m_qual_etof);
201 status = m_tuple4->addIndexedItem ("ec_toff_e", m_ngch,m_ec_toff_e);
202 status = m_tuple4->addIndexedItem ("ec_toff_mu",m_ngch,m_ec_toff_mu);
203 status = m_tuple4->addIndexedItem ("ec_toff_pi",m_ngch,m_ec_toff_pi);
204 status = m_tuple4->addIndexedItem ("ec_toff_k", m_ngch,m_ec_toff_k);
205 status = m_tuple4->addIndexedItem ("ec_toff_p", m_ngch,m_ec_toff_p);
206 status = m_tuple4->addIndexedItem ("ec_tsig_e", m_ngch,m_ec_tsig_e);
207 status = m_tuple4->addIndexedItem ("ec_tsig_mu",m_ngch,m_ec_tsig_mu);
208 status = m_tuple4->addIndexedItem ("ec_tsig_pi",m_ngch,m_ec_tsig_pi);
209 status = m_tuple4->addIndexedItem ("ec_tsig_k", m_ngch,m_ec_tsig_k);
210 status = m_tuple4->addIndexedItem ("ec_tsig_p", m_ngch,m_ec_tsig_p);
211 status = m_tuple4->addIndexedItem ("ec_tof", m_ngch,m_ec_tof);
212 status = m_tuple4->addIndexedItem ("ptot_btof1",m_ngch,m_ptot_btof1);
213 status = m_tuple4->addIndexedItem ("cntr_btof1",m_ngch,m_cntr_btof1);
214 status = m_tuple4->addIndexedItem ("te_btof1", m_ngch,m_te_btof1);
215 status = m_tuple4->addIndexedItem ("tmu_btof1", m_ngch,m_tmu_btof1);
216 status = m_tuple4->addIndexedItem ("tpi_btof1", m_ngch,m_tpi_btof1);
217 status = m_tuple4->addIndexedItem ("tk_btof1", m_ngch,m_tk_btof1);
218 status = m_tuple4->addIndexedItem ("tp_btof1", m_ngch,m_tp_btof1);
219 status = m_tuple4->addIndexedItem ("ph_btof1", m_ngch,m_ph_btof1);
220 status = m_tuple4->addIndexedItem ("zhit_btof1",m_ngch,m_zhit_btof1);
221 status = m_tuple4->addIndexedItem ("qual_btof1",m_ngch,m_qual_btof1);
222 status = m_tuple4->addIndexedItem ("b1_toff_e", m_ngch,m_b1_toff_e);
223 status = m_tuple4->addIndexedItem ("b1_toff_mu",m_ngch,m_b1_toff_mu);
224 status = m_tuple4->addIndexedItem ("b1_toff_pi",m_ngch,m_b1_toff_pi);
225 status = m_tuple4->addIndexedItem ("b1_toff_k", m_ngch,m_b1_toff_k);
226 status = m_tuple4->addIndexedItem ("b1_toff_p", m_ngch,m_b1_toff_p);
227 status = m_tuple4->addIndexedItem ("b1_tsig_e", m_ngch,m_b1_tsig_e);
228 status = m_tuple4->addIndexedItem ("b1_tsig_mu",m_ngch,m_b1_tsig_mu);
229 status = m_tuple4->addIndexedItem ("b1_tsig_pi",m_ngch,m_b1_tsig_pi);
230 status = m_tuple4->addIndexedItem ("b1_tsig_k", m_ngch,m_b1_tsig_k);
231 status = m_tuple4->addIndexedItem ("b1_tsig_p", m_ngch,m_b1_tsig_p);
232 status = m_tuple4->addIndexedItem ("b1_tof", m_ngch,m_b1_tof);
233
234 status = m_tuple4->addIndexedItem ("mdedx_pid", m_ngch,m_dedx_pid);
235 status = m_tuple4->addIndexedItem ("mtof1_pid", m_ngch,m_tof1_pid);
236 status = m_tuple4->addIndexedItem ("mtof2_pid", m_ngch,m_tof2_pid);
237 status = m_tuple4->addIndexedItem ("mprob_pid", m_ngch,m_prob_pid);
238 status = m_tuple4->addIndexedItem ("mptrk_pid", m_ngch,m_ptrk_pid);
239 status = m_tuple4->addIndexedItem ("mcost_pid", m_ngch,m_cost_pid);
240 status = m_tuple4->addItem ("mpnp", m_pnp);
241 status = m_tuple4->addItem ("mpnm", m_pnm);
242
243 status = m_tuple4->addIndexedItem ("numHits", m_nggneu,m_numHits); // Total number of hits
244 status = m_tuple4->addIndexedItem ("secondmoment", m_nggneu,m_secondmoment);
245 status = m_tuple4->addIndexedItem ("mx", m_nggneu,m_x); // Shower coordinates and errors
246 status = m_tuple4->addIndexedItem ("my", m_nggneu,m_y);
247 status = m_tuple4->addIndexedItem ("mz", m_nggneu,m_z);
248 status = m_tuple4->addIndexedItem ("cosemc", m_nggneu,m_cosemc); // Shower Counter angles and errors
249 status = m_tuple4->addIndexedItem ("phiemc", m_nggneu,m_phiemc);
250 status = m_tuple4->addIndexedItem ("energy", m_nggneu,m_energy); // Total energy observed in Emc
251 status = m_tuple4->addIndexedItem ("eseed", m_nggneu,m_eSeed);
252 status = m_tuple4->addIndexedItem ("me9", m_nggneu,m_e3x3);
253 status = m_tuple4->addIndexedItem ("me25", m_nggneu,m_e5x5);
254 status = m_tuple4->addIndexedItem ("mlat", m_nggneu,m_lat);
255 status = m_tuple4->addIndexedItem ("ma20", m_nggneu,m_a20);
256 status = m_tuple4->addIndexedItem ("ma42", m_nggneu,m_a42);
257
258 }
259 else {
260 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
261 return StatusCode::FAILURE;
262 }
263 }
264
265 if(service("THistSvc", m_thsvc).isFailure()) {
266 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
267 return StatusCode::FAILURE;
268 }
269
270 // "DQAHist" is fixed
271 // "Rhopi" is control sample name, just as ntuple case.
272
273 TH1F* mrho0 = new TH1F("mrho0", "mass for #rho^{0}->#pi^{+}#pi^{-}", 160, 0.0, 3.2);
274 if(m_thsvc->regHist("/DQAHist/Rhopi/mrho0", mrho0).isFailure()) {
275 log << MSG::ERROR << "Couldn't register mrho0" << endreq;
276 }
277
278 TH1F* mrhop = new TH1F("mrhop", "mass for #rho^{+}->#pi^{+}#pi^{0}", 160, 0.0,3.2);
279 if(m_thsvc->regHist("/DQAHist/Rhopi/mrhop", mrhop).isFailure()) {
280 log << MSG::ERROR << "Couldn't register mrhop" << endreq;
281 }
282
283 TH1F* mrhom = new TH1F("mrhom", "mass for #rho^{-}->#pi^{-}#pi^{0}", 160, 0.0, 3.2);
284 if(m_thsvc->regHist("/DQAHist/Rhopi/mrhom", mrhom).isFailure()) {
285 log << MSG::ERROR << "Couldn't register mrhom" << endreq;
286 }
287
288
289 TH1F* mpi0 = new TH1F("mpi0", "mass for #pi^{0}->#gamma #gamma", 50,0.08, 0.18);
290 if(m_thsvc->regHist("/DQAHist/Rhopi/mpi0", mpi0).isFailure()) {
291 log << MSG::ERROR << "Couldn't register mpi0" << endreq;
292 }
293
294 //
295 //--------end of book--------
296 //
297
298 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
299 return StatusCode::SUCCESS;
300
301}
302
303// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
304StatusCode DQARhopi::execute() {
305
306// std::cout << "execute()" << std::endl;
307
308 MsgStream log(msgSvc(), name());
309 log << MSG::INFO << "in execute()" << endreq;
310
311 setFilterPassed(false);
312
313 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
314 int run=eventHeader->runNumber();
315 int event=eventHeader->eventNumber();
316 log << MSG::DEBUG <<"run, evtnum = "
317 << run << " , "
318 << event <<endreq;
319
320 m_run = eventHeader->runNumber();
321 m_rec = eventHeader->eventNumber();
322
323// cout<<"event "<<event<<endl;
324 Ncut0++;
325 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
326 log << MSG::INFO << "get event tag OK" << endreq;
327 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
328 << evtRecEvent->totalCharged() << " , "
329 << evtRecEvent->totalNeutral() << " , "
330 << evtRecEvent->totalTracks() <<endreq;
331
332 m_nch = evtRecEvent->totalCharged();
333 m_nneu = evtRecEvent->totalNeutral();
334
335 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
336
337 //
338 // check x0, y0, z0, r0
339 // suggest cut: |z0|<5 && r0<1
340 //
341 Vint iGood, ipip, ipim, ipnp,ipnm;
342 iGood.clear();
343 ipip.clear();
344 ipim.clear();
345 ipnp.clear();
346 ipnm.clear();
347 Vp4 ppip, ppim , ppnp, ppnm;
348 ppip.clear();
349 ppim.clear();
350 ppnp.clear();
351 ppnm.clear();
352
353 Hep3Vector xorigin(0,0,0);
354
355
356 //if (m_reader.isRunNumberValid(runNo)) {
357 IVertexDbSvc* vtxsvc;
358 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
359 if(vtxsvc->isVertexValid()){
360 double* dbv = vtxsvc->PrimaryVertex();
361 double* vv = vtxsvc->SigmaPrimaryVertex();
362 xorigin.setX(dbv[0]);
363 xorigin.setY(dbv[1]);
364 xorigin.setZ(dbv[2]);
365 }
366
367 int nCharge = 0;
368 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
369 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
370 if(!(*itTrk)->isMdcTrackValid()) continue;
371 if (!(*itTrk)->isMdcKalTrackValid()) continue;
372 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
373 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
374
375 double pch =mdcTrk->p();
376 double x0 =mdcTrk->x();
377 double y0 =mdcTrk->y();
378 double z0 =mdcTrk->z();
379 double phi0=mdcTrk->helix(1);
380 double xv=xorigin.x();
381 double yv=xorigin.y();
382 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
383 double m_vx0 = x0;
384 double m_vy0 = y0;
385 double m_vz0 = z0-xorigin.z();
386 double m_vr0 = Rxy;
387 double m_Vct=cos(mdcTrk->theta());
388// m_tuple1->write();
389 if(fabs(m_vz0) >= m_vz0cut) continue;
390 if(m_vr0 >= m_vr0cut) continue;
391 if(fabs(m_Vct)>=m_cthcut) continue;
392
393 iGood.push_back((*itTrk)->trackId());
394 nCharge += mdcTrk->charge();
395
396 }
397
398 //
399 // Finish Good Charged Track Selection
400 //
401 int nGood = iGood.size();
402 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
403 if((nGood != 2)||(nCharge!=0)){
404 return StatusCode::SUCCESS;
405 }
406 Ncut1++;
407
408 Vint iGam;
409 iGam.clear();
410 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
411 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
412 if(!(*itTrk)->isEmcShowerValid()) continue;
413 RecEmcShower *emcTrk = (*itTrk)->emcShower();
414 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
415 // find the nearest charged track
416 double dthe = 200.;
417 double dphi = 200.;
418 double dang = 200.;
419 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
420 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
421 if(!(*jtTrk)->isExtTrackValid()) continue;
422 RecExtTrack *extTrk = (*jtTrk)->extTrack();
423 if(extTrk->emcVolumeNumber() == -1) continue;
424 Hep3Vector extpos = extTrk->emcPosition();
425 // double ctht = extpos.cosTheta(emcpos);
426 double angd = extpos.angle(emcpos);
427 double thed = extpos.theta() - emcpos.theta();
428 double phid = extpos.deltaPhi(emcpos);
429 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
430 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
431
432 if(fabs(thed) < fabs(dthe)) dthe = thed;
433 if(fabs(phid) < fabs(dphi)) dphi = phid;
434 if(angd < dang) dang = angd;
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 double m_dthe = dthe;
442 double m_dphi = dphi;
443 double m_dang = dang;
444 double m_eraw = eraw;
445// m_tuple2->write();
446 if(eraw < m_energyThreshold) continue;
447 if(dang < m_gammaAngCut) continue;
448 //
449 // good photon cut will be set here
450 //
451 iGam.push_back((*itTrk)->trackId());
452 }
453
454 //
455 // Finish Good Photon Selection
456 //
457 int nGam = iGam.size();
458
459 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
460 if(nGam<2){
461 return StatusCode::SUCCESS;
462 }
463 Ncut2++;
464
465
466 //
467 // Assign 4-momentum to each photon
468 //
469
470 Vp4 pGam;
471 pGam.clear();
472 for(int i = 0; i < nGam; i++) {
473 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
474 RecEmcShower* emcTrk = (*itTrk)->emcShower();
475 double eraw = emcTrk->energy();
476 double phi = emcTrk->phi();
477 double the = emcTrk->theta();
478 HepLorentzVector ptrk;
479 ptrk.setPx(eraw*sin(the)*cos(phi));
480 ptrk.setPy(eraw*sin(the)*sin(phi));
481 ptrk.setPz(eraw*cos(the));
482 ptrk.setE(eraw);
483
484// ptrk = ptrk.boost(-0.011,0,0);// boost to cms
485
486 pGam.push_back(ptrk);
487 }
488
489 log << MSG::DEBUG << "liuf1 Good Photon " <<endreq;
490
491
492 for(int i = 0; i < nGood; i++) {//for rhopi without PID
493 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
494 if(!(*itTrk)->isMdcTrackValid()) continue;
495 if (!(*itTrk)->isMdcKalTrackValid()) continue;
496 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
497 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
498 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
499
500 if(mdcTrk->charge() >0) {
501 ipip.push_back(iGood[i]);
502 HepLorentzVector ptrk;
503 ptrk.setPx(mdcKalTrk->px());
504 ptrk.setPy(mdcKalTrk->py());
505 ptrk.setPz(mdcKalTrk->pz());
506 double p3 = ptrk.mag();
507 ptrk.setE(sqrt(p3*p3+mpi*mpi));
508 ppip.push_back(ptrk);
509 } else {
510 ipim.push_back(iGood[i]);
511 HepLorentzVector ptrk;
512 ptrk.setPx(mdcKalTrk->px());
513 ptrk.setPy(mdcKalTrk->py());
514 ptrk.setPz(mdcKalTrk->pz());
515 double p3 = ptrk.mag();
516 ptrk.setE(sqrt(p3*p3+mpi*mpi));
517 ppim.push_back(ptrk);
518 }
519 }// without PID
520
521 int npip = ipip.size();
522 int npim = ipim.size();
523 log << MSG::DEBUG << "num of pion "<< ipip.size()<<","<<ipim.size()<<endreq;
524 Ncut3++;
525
526
527
528 //***********************************************//
529 // the angle between the two charged tracks //
530 //***********************************************//
531
532 int langcc=0;
533 double dthec = 200.;
534 double dphic = 200.;
535 double dangc = 200.;
536 for(int i=0; i < npip; i++) {
537 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + ipip[i] ;
538 RecMdcTrack* mdcTrk1 = (*itTrk)->mdcTrack();
539 RecMdcKalTrack* mdcKalTrk1 = (*itTrk)->mdcKalTrack();
540 Hep3Vector emcpos(ppip[i].px(), ppip[i].py(), ppip[i].pz());
541
542 for(int j = 0; j < npim; j++) {
543 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
544 RecMdcTrack* mdcTrk2 = (*jtTrk)->mdcTrack();
545 RecMdcKalTrack* mdcKalTrk2 = (*jtTrk)->mdcKalTrack();
546
547 HepLorentzVector pippim=ppip[i]+ppim[j];
548
549 Hep3Vector extpos(ppim[j].px(), ppim[j].py(), ppim[j].pz());
550
551 double angd = extpos.angle(emcpos);
552 double thed = extpos.theta() - emcpos.theta();
553 double phid = extpos.deltaPhi(emcpos);
554 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
555 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
556
557 m_dthec[langcc] = thed * 180 / (CLHEP::pi);
558 m_dphic[langcc] = phid * 180 / (CLHEP::pi);
559 m_dangc[langcc] = angd * 180 / (CLHEP::pi);
560 m_mspippim[langcc] =pippim.m();
561 langcc++;
562 }
563 }
564 m_nangecc=langcc;
565
566 //
567 // Loop each gamma pair, check ppi0 and pTot
568 //
569
570 double m_m2gg,m_momentpi0;
571 HepLorentzVector pTot,p2g;
572
573 log << MSG::DEBUG << "liuf2 Good Photon " <<langcc<<endreq;
574 //******************************************************//
575 // asign the momentum of MDC and KALFIT //
576 // the deposite energy of EMC for charged tracks //
577 //******************************************************//
578 double m_momentpip,m_momentpim,extmot1,extmot2;
579 double emcTg1=0.0;
580 double emcTg2=0.0;
581 double nlaypi1=0;
582 double nhit1=0;
583 double nlaypi2=0;
584 double nhit2=0;
585
586 EvtRecTrackIterator itTrk11 = evtRecTrkCol->begin() + ipip[0];
587 RecMdcTrack* mdcTrk11 = (*itTrk11)->mdcTrack();
588 RecMdcKalTrack* mdcKalTrk11 = (*itTrk11)->mdcKalTrack();
589 RecEmcShower* emcTrk11 = (*itTrk11)->emcShower();
590 RecMucTrack *mucTrk11=(*itTrk11)->mucTrack();
591 double phi01=mdcTrk11->helix(1);
592
593 EvtRecTrackIterator itTrk12 = evtRecTrkCol->begin() + ipim[0];
594 RecMdcTrack* mdcTrk12 = (*itTrk12)->mdcTrack();
595 RecMdcKalTrack* mdcKalTrk12 = (*itTrk12)->mdcKalTrack();
596 RecEmcShower* emcTrk12 = (*itTrk12)->emcShower();
597 RecMucTrack *mucTrk12=(*itTrk12)->mucTrack();
598 double phi02=mdcTrk12->helix(1);
599
600 m_vxpin = mdcTrk11->x();
601 m_vypin = mdcTrk11->y();
602 m_vzpin = mdcTrk11->z()-xorigin.z();
603 m_vrpin = fabs((mdcTrk11->x()-xorigin.x())*cos(phi01)+(mdcTrk11->y()-xorigin.y())*sin(phi01));
604 m_costhepin =cos(mdcTrk11->theta());
605
606 m_momentpip=mdcTrk11->p();
607 m_ppx =mdcTrk11->px();
608 m_ppy =mdcTrk11->py();
609 m_ppz =mdcTrk11->pz();
610
611 m_vxp = mdcKalTrk11->x();
612 m_vyp = mdcKalTrk11->y();
613 m_vzp = mdcKalTrk11->z()-xorigin.z();
614 m_vrp = fabs((mdcKalTrk11->x()-xorigin.x())*cos(phi01)+(mdcKalTrk11->y()-xorigin.y())*sin(phi01));
615 m_costhep =cos(mdcKalTrk11->theta());
616
617 extmot1=mdcKalTrk11->p();
618 m_ppxkal =mdcKalTrk11->px();
619 m_ppykal =mdcKalTrk11->py();
620 m_ppzkal =mdcKalTrk11->pz();
621
622 m_vxmin = mdcTrk12->x();
623 m_vymin = mdcTrk12->y();
624 m_vzmin = mdcTrk12->z();
625 m_vrmin = fabs((mdcTrk12->x()-xorigin.x())*cos(phi02)+(mdcTrk12->y()-xorigin.y())*sin(phi02));
626 m_costhemin =cos(mdcTrk12->theta());
627
628 m_momentpim=mdcTrk12->p();
629 m_mpx =mdcTrk12->px();
630 m_mpy =mdcTrk12->py();
631 m_mpz =mdcTrk12->pz();
632
633 m_vxm = mdcKalTrk12->x();
634 m_vym = mdcKalTrk12->y();
635 m_vzm = mdcKalTrk12->z();
636 m_vrm = fabs((mdcKalTrk12->x()-xorigin.x())*cos(phi02)+(mdcKalTrk12->y()-xorigin.y())*sin(phi02));
637 m_costhem =cos(mdcKalTrk12->theta());
638
639 extmot2 =mdcKalTrk12->p();
640 m_mpxkal =mdcKalTrk12->px();
641 m_mpykal =mdcKalTrk12->py();
642 m_mpzkal =mdcKalTrk12->pz();
643
644 if((*itTrk11)->isEmcShowerValid()){
645 emcTg1 = emcTrk11->energy();
646 }
647 if((*itTrk12)->isEmcShowerValid()){
648 emcTg2 = emcTrk12->energy();
649 }
650 if((*itTrk11)->isMucTrackValid()){
651 nlaypi1=mucTrk11->numLayers();
652 nhit1 =mucTrk11->numHits();
653 }
654 if((*itTrk12)->isMucTrackValid()){
655 nlaypi2=mucTrk12->numLayers();
656 nhit2 =mucTrk12->numHits();
657 }
658
659 m_laypi1=nlaypi1;
660 m_hit1 =nhit1;
661 m_laypi2=nlaypi2;
662 m_hit2 =nhit2;
663
664 log << MSG::DEBUG << "liuf3 Good Photon " << ipim[0] <<endreq;
665
666 RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
667 RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
668
669 log << MSG::DEBUG << "liuf4 Good Photon " <<endreq;
670
671 WTrackParameter wvpipTrk, wvpimTrk,wvkpTrk, wvkmTrk;
672 wvpipTrk = WTrackParameter(mpi, pipTrk->getZHelix(), pipTrk->getZError());
673 wvpimTrk = WTrackParameter(mpi, pimTrk->getZHelix(), pimTrk->getZError());
674
675 wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());//kaon
676 wvkmTrk = WTrackParameter(mk, pimTrk->getZHelixK(), pimTrk->getZErrorK());//kaon
677
678/* Default is pion, for other particles:
679 wvppTrk = WTrackParameter(mp, pipTrk->getZHelixP(), pipTrk->getZErrorP());//proton
680 wvmupTrk = WTrackParameter(mmu, pipTrk->getZHelixMu(), pipTrk->getZErrorMu());//muon
681 wvepTrk = WTrackParameter(me, pipTrk->getZHelixE(), pipTrk->getZErrorE());//electron
682 wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());//kaon
683*/
684 //
685 // Test vertex fit
686 //
687
688 HepPoint3D vx(0., 0., 0.);
689 HepSymMatrix Evx(3, 0);
690 double bx = 1E+6;
691 double by = 1E+6;
692 double bz = 1E+6;
693 Evx[0][0] = bx*bx;
694 Evx[1][1] = by*by;
695 Evx[2][2] = bz*bz;
696
697 VertexParameter vxpar;
698 vxpar.setVx(vx);
699 vxpar.setEvx(Evx);
700
701 //****************************************************//
702 // Test vertex fit //
703 //****************************************************//
704
705 VertexFit* vtxfit = VertexFit::instance();
706
707 //****************************************************//
708 // if the charged particle is kaon //
709 //****************************************************//
710 double chi5=9999.0;
711 double jkpi0=-0.5;
712 double jkpkm=0.0;
713 double jkpp0=0.0;
714 double jkmp0=0.0;
715 vtxfit->init();
716 vtxfit->AddTrack(0, wvkpTrk);
717 vtxfit->AddTrack(1, wvkmTrk);
718 vtxfit->AddVertex(0, vxpar,0, 1);
719 if(vtxfit->Fit(0)) {
720 vtxfit->Swim(0);
721 WTrackParameter wkp = vtxfit->wtrk(0);
722 WTrackParameter wkm = vtxfit->wtrk(1);
723
724 KinematicFit * kmfit = KinematicFit::instance();
725
726 //
727 // Apply Kinematic 4C fit
728 //
729
730 double chisq = 9999.;
731 for(int i = 0; i < nGam-1; i++) {
732 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
733 for(int j = i+1; j < nGam; j++) {
734 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
735 kmfit->init();
736 kmfit->setDynamicerror(1);
737 kmfit->AddTrack(0, wkp);
738 kmfit->AddTrack(1, wkm);
739 kmfit->AddTrack(2, 0.0, g1Trk);
740 kmfit->AddTrack(3, 0.0, g2Trk);
741 kmfit->AddFourMomentum(0, ecms);
742 bool oksq = kmfit->Fit();
743 if(oksq) {
744 double chi2 = kmfit->chisq();
745 if(chi2 < chi5) {
746 HepLorentzVector kpi0 = kmfit->pfit(2) + kmfit->pfit(3);
747 HepLorentzVector kpkm = kmfit->pfit(0) + kmfit->pfit(1);
748 HepLorentzVector kpp0 = kmfit->pfit(0) + kpi0;
749 HepLorentzVector kmp0 = kmfit->pfit(1) + kpi0;
750 chi5 = kmfit->chisq();
751 jkpi0 = kpi0.m();
752 jkpkm= kpkm.m();
753 jkpp0= kpp0.m();
754 jkmp0= kmp0.m();
755 }
756 }
757 }
758 }
759 } // vetex is true
760
761 //****************************************************//
762 // if the charged particles are pions for real //
763 //****************************************************//
764
765 vtxfit->init();
766 vtxfit->AddTrack(0, wvpipTrk);
767 vtxfit->AddTrack(1, wvpimTrk);
768 vtxfit->AddVertex(0, vxpar,0, 1);
769 if(!vtxfit->Fit(0)) return SUCCESS;
770 vtxfit->Swim(0);
771
772 WTrackParameter wpip = vtxfit->wtrk(0);
773 WTrackParameter wpim = vtxfit->wtrk(1);
774
775 log << MSG::DEBUG << "liuf5 Good Photon " <<endreq;
776
777 KinematicFit * kmfit = KinematicFit::instance();
778
779 //
780 // Apply Kinematic 4C fit
781 //
782
783 double chisq = 9999.;
784 int ig1 = -1;
785 int ig2 = -1;
786 HepLorentzVector gg1,gg2;
787 for(int i = 0; i < nGam-1; i++) {
788 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
789 for(int j = i+1; j < nGam; j++) {
790 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
791 kmfit->init();
792 kmfit->setDynamicerror(1);
793 kmfit->AddTrack(0, wpip);
794 kmfit->AddTrack(1, wpim);
795 kmfit->AddTrack(2, 0.0, g1Trk);
796 kmfit->AddTrack(3, 0.0, g2Trk);
797 kmfit->AddFourMomentum(0, ecms);
798 bool oksq = kmfit->Fit();
799 if(oksq) {
800 double chi2 = kmfit->chisq();
801 if(chi2 < chisq) {
802 chisq = chi2;
803 ig1 = iGam[i];
804 ig2 = iGam[j];
805 gg1 = pGam[i];
806 gg2 = pGam[j];
807 }
808 }
809 }
810 }
811
812 p2g = gg1 + gg2;
813 m_pmax=gg1.e()>gg2.e()?gg1.e():gg2.e();
814 m_m2gg = p2g.m();
815 m_cosang=(gg1.e()-gg2.e())/p2g.rho();
816 m_momentpi0=sqrt(p2g.px()*p2g.px()+p2g.py()*p2g.py()+p2g.pz()*p2g.pz());
817 log << MSG::DEBUG << " chisq for 4c " << chisq <<endreq;
818 if(chisq > 200) {
819 return StatusCode::SUCCESS;
820 }
821
822// select charge track and nneu track
823 Vint jGood;
824 jGood.clear();
825 jGood.push_back(ipip[0]);
826 jGood.push_back(ipim[0]);
827 m_ngch = jGood.size();
828
829 Vint jGgam;
830 jGgam.clear();
831 jGgam.push_back(ig1);
832 jGgam.push_back(ig2);
833 m_nggneu=jGgam.size();
834
835 HepLorentzVector ppip1=ppip[0];
836 HepLorentzVector ppim1=ppim[0];
837
838 HepLorentzVector Ppipboost = ppip1.boost(u_cms);
839 HepLorentzVector Ppimboost = ppim1.boost(u_cms);
840 Hep3Vector p3pi1 = Ppipboost.vect(); //pi1
841 Hep3Vector p3pi2 = Ppimboost.vect(); //pi2
842 m_anglepm = (p3pi1.angle(p3pi2))* 180 / (CLHEP::pi);
843
844//*******************************************************//
845
846 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
847 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
848 kmfit->init();
849 kmfit->setDynamicerror(1);
850 kmfit->AddTrack(0, wpip);
851 kmfit->AddTrack(1, wpim);
852 kmfit->AddTrack(2, 0.0, g1Trk);
853 kmfit->AddTrack(3, 0.0, g2Trk);
854 kmfit->AddFourMomentum(0, ecms);
855 bool oksq = kmfit->Fit();
856 if(!oksq) return SUCCESS;
857
858 HepLorentzVector ppi0 = kmfit->pfit(2) + kmfit->pfit(3);
859 HepLorentzVector prho0 = kmfit->pfit(0) + kmfit->pfit(1);
860 HepLorentzVector prhop = kmfit->pfit(0) + ppi0;
861 HepLorentzVector prhom = kmfit->pfit(1) + ppi0;
862 HepLorentzVector pgam2pi1 = prho0 + kmfit->pfit(2);
863 HepLorentzVector pgam2pi2 = prho0 + kmfit->pfit(3);
864 HepLorentzVector pgam11 =kmfit->pfit(2);
865 HepLorentzVector pgam12 =kmfit->pfit(3);
866
867/*
868 HepLorentzVector ppi0_a=ppi0;
869 HepLorentzVector pgam11_a =pgam11;
870 HepLorentzVector pgam12_a =pgam12;
871
872 Hep3Vector pi0_cms = -ppi0_a.boostVector();
873 HepLorentzVector pgam1 = pgam11_a.boost(pi0_cms);
874 HepLorentzVector pgam2 = pgam12_a.boost(pi0_cms);
875*/
876 m_chi1 = kmfit->chisq();
877 m_mpi0 = ppi0.m();
878 m_prho0= prho0.m();
879 m_prhop= prhop.m();
880 m_prhom= prhom.m();
881 m_good = nGood;
882 m_gam = nGam;
883
884 m_pip = npip;
885 m_pim = npim;
886 m_2gam = m_m2gg;
887 m_outpi0=m_momentpi0;
888 m_outpip=m_momentpip;
889 m_outpim=m_momentpim;
890 m_enpip =emcTg1;
891 m_dcpip =extmot1;
892 m_enpim =emcTg2;
893 m_dcpim =extmot2;
894 m_pipf=kmfit->pfit(0).rho();
895 m_pimf=kmfit->pfit(1).rho();
896 m_pi0f=ppi0.rho();
897
898 m_chi5=chi5;
899 m_kpi0=jkpi0;
900 m_kpkm=jkpkm;
901 m_kpp0=jkpp0;
902 m_kmp0=jkmp0;
903 m_pgam2pi1=pgam2pi1.m();
904 m_pgam2pi2=pgam2pi2.m();
905 cosva1=kmfit->pfit(2).rho();
906 cosva2=kmfit->pfit(3).rho();
907// m_pe1 =pe1;
908// m_pe2 =pe2;
909
910//
911// fill information of dedx and tof
912//
913
914 //
915 // check dedx infomation
916 //
917
918 for(int ii = 0; ii < m_ngch; ii++) {
919// dedx
920 m_ptrk[ii] = 9999.0;
921 m_chie[ii] = 9999.0;
922 m_chimu[ii] = 9999.0;
923 m_chipi[ii] = 9999.0;
924 m_chik[ii] = 9999.0;
925 m_chip[ii] = 9999.0;
926 m_ghit[ii] = 9999.0;
927 m_thit[ii] = 9999.0;
928 m_probPH[ii] = 9999.0;
929 m_normPH[ii] = 9999.0;
930
931//endtof
932 m_cntr_etof[ii] = 9999.0;
933 m_ptot_etof[ii] = 9999.0;
934 m_ph_etof[ii] = 9999.0;
935 m_rhit_etof[ii] = 9999.0;
936 m_qual_etof[ii] = 9999.0;
937 m_te_etof[ii] = 9999.0;
938 m_tmu_etof[ii] = 9999.0;
939 m_tpi_etof[ii] = 9999.0;
940 m_tk_etof[ii] = 9999.0;
941 m_tp_etof[ii] = 9999.0;
942 m_ec_tof[ii] = 9999.0;
943 m_ec_toff_e[ii] = 9999.0;
944 m_ec_toff_mu[ii] = 9999.0;
945 m_ec_toff_pi[ii] = 9999.0;
946 m_ec_toff_k[ii] = 9999.0;
947 m_ec_toff_p[ii] = 9999.0;
948 m_ec_tsig_e[ii] = 9999.0;
949 m_ec_tsig_mu[ii] = 9999.0;
950 m_ec_tsig_pi[ii] = 9999.0;
951 m_ec_tsig_k[ii] = 9999.0;
952 m_ec_tsig_p[ii] = 9999.0;
953
954// barrel tof
955 m_cntr_btof1[ii] = 9999.0;
956 m_ptot_btof1[ii] = 9999.0;
957 m_ph_btof1[ii] = 9999.0;
958 m_zhit_btof1[ii] = 9999.0;
959 m_qual_btof1[ii] = 9999.0;
960 m_te_btof1[ii] = 9999.0;
961 m_tmu_btof1[ii] = 9999.0;
962 m_tpi_btof1[ii] = 9999.0;
963 m_tk_btof1[ii] = 9999.0;
964 m_tp_btof1[ii] = 9999.0;
965 m_b1_tof[ii] = 9999.0;
966 m_b1_toff_e[ii] = 9999.0;
967 m_b1_toff_mu[ii] = 9999.0;
968 m_b1_toff_pi[ii] = 9999.0;
969 m_b1_toff_k[ii] = 9999.0;
970 m_b1_toff_p[ii] = 9999.0;
971 m_b1_tsig_e[ii] = 9999.0;
972 m_b1_tsig_mu[ii] = 9999.0;
973 m_b1_tsig_pi[ii] = 9999.0;
974 m_b1_tsig_k[ii] = 9999.0;
975 m_b1_tsig_p[ii] = 9999.0;
976//pid
977 m_dedx_pid[ii] = 9999.0;
978 m_tof1_pid[ii] = 9999.0;
979 m_tof2_pid[ii] = 9999.0;
980 m_prob_pid[ii] = 9999.0;
981 m_ptrk_pid[ii] = 9999.0;
982 m_cost_pid[ii] = 9999.0;
983 }
984
985
986 int indx0=0;
987 for(int i = 0; i < m_ngch; i++) {
988 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
989 if(!(*itTrk)->isMdcTrackValid()) continue;
990 if(!(*itTrk)->isMdcDedxValid())continue;
991 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
992 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
993 m_ptrk[indx0] = mdcTrk->p();
994
995 m_chie[indx0] = dedxTrk->chiE();
996 m_chimu[indx0] = dedxTrk->chiMu();
997 m_chipi[indx0] = dedxTrk->chiPi();
998 m_chik[indx0] = dedxTrk->chiK();
999 m_chip[indx0] = dedxTrk->chiP();
1000 m_ghit[indx0] = dedxTrk->numGoodHits();
1001 m_thit[indx0] = dedxTrk->numTotalHits();
1002 m_probPH[indx0] = dedxTrk->probPH();
1003 m_normPH[indx0] = dedxTrk->normPH();
1004 indx0++;
1005 }
1006
1007
1008 //
1009 // check TOF infomation
1010 //
1011
1012
1013 int indx1=0;
1014 for(int i = 0; i < m_ngch; i++) {
1015 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1016 if(!(*itTrk)->isMdcTrackValid()) continue;
1017 if(!(*itTrk)->isTofTrackValid()) continue;
1018
1019 RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
1020 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
1021
1022 double ptrk = mdcTrk->p();
1023 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1024 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
1025 TofHitStatus *status = new TofHitStatus;
1026 status->setStatus((*iter_tof)->status());
1027 if(!(status->is_barrel())){//endcap
1028 if( !(status->is_counter()) ) continue; // ?
1029 if( status->layer()!=1 ) continue;//layer1
1030 double path=(*iter_tof)->path(); // ?
1031 double tof = (*iter_tof)->tof();
1032 double ph = (*iter_tof)->ph();
1033 double rhit = (*iter_tof)->zrhit();
1034 double qual = 0.0 + (*iter_tof)->quality();
1035 double cntr = 0.0 + (*iter_tof)->tofID();
1036 double texp[5];
1037 double tsig[5];
1038 for(int j = 0; j < 5; j++) {//0 e, 1 mu, 2 pi, 3 K, 4 p
1039 texp[j] = (*iter_tof)->texp(j);
1040// tsig[j] = (*iter_tof)->sigma(j);
1041// toffset[j] = (*iter_tof)->offset(j);
1042 }
1043 m_cntr_etof[indx1] = cntr;
1044 m_ptot_etof[indx1] = ptrk;
1045 m_ph_etof[indx1] = ph;
1046 m_rhit_etof[indx1] = rhit;
1047 m_qual_etof[indx1] = qual;
1048 m_te_etof[indx1] = tof - texp[0];
1049 m_tmu_etof[indx1] = tof - texp[1];
1050 m_tpi_etof[indx1] = tof - texp[2];
1051 m_tk_etof[indx1] = tof - texp[3];
1052 m_tp_etof[indx1] = tof - texp[4];
1053
1054 m_ec_tof[indx1] = tof;
1055
1056 m_ec_toff_e[indx1] = (*iter_tof)->toffset(0);
1057 m_ec_toff_mu[indx1] = (*iter_tof)->toffset(1);
1058 m_ec_toff_pi[indx1] = (*iter_tof)->toffset(2);
1059 m_ec_toff_k[indx1] = (*iter_tof)->toffset(3);
1060 m_ec_toff_p[indx1] = (*iter_tof)->toffset(4);
1061
1062 m_ec_tsig_e[indx1] = (*iter_tof)->sigma(0);
1063 m_ec_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1064 m_ec_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1065 m_ec_tsig_k[indx1] = (*iter_tof)->sigma(3);
1066 m_ec_tsig_p[indx1] = (*iter_tof)->sigma(4);
1067
1068 }
1069 else {//barrel
1070 if( !(status->is_cluster()) ) continue; // ?
1071 double path=(*iter_tof)->path(); // ?
1072 double tof = (*iter_tof)->tof();
1073 double ph = (*iter_tof)->ph();
1074 double rhit = (*iter_tof)->zrhit();
1075 double qual = 0.0 + (*iter_tof)->quality();
1076 double cntr = 0.0 + (*iter_tof)->tofID();
1077 double texp[5];
1078 for(int j = 0; j < 5; j++) {
1079 texp[j] = (*iter_tof)->texp(j);
1080 }
1081 m_cntr_btof1[indx1] = cntr;
1082 m_ptot_btof1[indx1] = ptrk;
1083 m_ph_btof1[indx1] = ph;
1084 m_zhit_btof1[indx1] = rhit;
1085 m_qual_btof1[indx1] = qual;
1086 m_te_btof1[indx1] = tof - texp[0];
1087 m_tmu_btof1[indx1] = tof - texp[1];
1088 m_tpi_btof1[indx1] = tof - texp[2];
1089 m_tk_btof1[indx1] = tof - texp[3];
1090 m_tp_btof1[indx1] = tof - texp[4];
1091
1092 m_b1_tof[indx1] = tof;
1093
1094 m_b1_toff_e[indx1] = (*iter_tof)->toffset(0);
1095 m_b1_toff_mu[indx1] = (*iter_tof)->toffset(1);
1096 m_b1_toff_pi[indx1] = (*iter_tof)->toffset(2);
1097 m_b1_toff_k[indx1] = (*iter_tof)->toffset(3);
1098 m_b1_toff_p[indx1] = (*iter_tof)->toffset(4);
1099
1100 m_b1_tsig_e[indx1] = (*iter_tof)->sigma(0);
1101 m_b1_tsig_mu[indx1] = (*iter_tof)->sigma(1);
1102 m_b1_tsig_pi[indx1] = (*iter_tof)->sigma(2);
1103 m_b1_tsig_k[indx1] = (*iter_tof)->sigma(3);
1104 m_b1_tsig_p[indx1] = (*iter_tof)->sigma(4);
1105
1106 }
1107 delete status;
1108 }
1109 indx1++;
1110 } // loop all charged track
1111
1112 //
1113 // Assign 4-momentum to each charged track
1114 //
1115 int indx2=0;
1116 ParticleID *pid = ParticleID::instance();
1117 for(int i = 0; i < m_ngch; i++) {
1118 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1119 // if(pid) delete pid;
1120 pid->init();
1121 pid->setMethod(pid->methodProbability());
1122// pid->setMethod(pid->methodLikelihood()); //for Likelihood Method
1123
1124 pid->setChiMinCut(4);
1125 pid->setRecTrack(*itTrk);
1126 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID sub-system
1127 pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
1128 // pid->identify(pid->onlyPion());
1129 // pid->identify(pid->onlyKaon());
1130 pid->calculate();
1131 if(!(pid->IsPidInfoValid())) continue;
1132 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
1133
1134 m_dedx_pid[indx2] = pid->chiDedx(2);
1135 m_tof1_pid[indx2] = pid->chiTof1(2);
1136 m_tof2_pid[indx2] = pid->chiTof2(2);
1137 m_prob_pid[indx2] = pid->probPion();
1138
1139// if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
1140// if(pid->probPion() < 0.001) continue;
1141// if(pid->pdf(2)<pid->pdf(3)) continue; // for Likelihood Method(0=electron 1=muon 2=pion 3=kaon 4=proton)
1142
1143 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use RecMdcKalTrack substitute RecMdcTrack
1144 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default setting is pion
1145
1146 m_ptrk_pid[indx2] = mdcKalTrk->p();
1147 m_cost_pid[indx2] = cos(mdcTrk->theta());
1148
1149 if(mdcTrk->charge() >0 && pid->probPion() > pid->probKaon()) {
1150 ipnp.push_back(jGood[i]);
1151 HepLorentzVector ptrk;
1152 ptrk.setPx(mdcKalTrk->px());
1153 ptrk.setPy(mdcKalTrk->py());
1154 ptrk.setPz(mdcKalTrk->pz());
1155 double p3 = ptrk.mag();
1156 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1157// ptrk = ptrk.boost(-0.011,0,0);//boost to cms
1158
1159 ppnp.push_back(ptrk);
1160 } //plus charge with with PID
1161 if(mdcTrk->charge() <0 && pid->probPion() > pid->probKaon()) {
1162 ipnm.push_back(jGood[i]);
1163 HepLorentzVector ptrk;
1164 ptrk.setPx(mdcKalTrk->px());
1165 ptrk.setPy(mdcKalTrk->py());
1166 ptrk.setPz(mdcKalTrk->pz());
1167 double p3 = ptrk.mag();
1168 ptrk.setE(sqrt(p3*p3+mpi*mpi));
1169
1170// ptrk = ptrk.boost(-0.011,0,0);//boost to cms
1171
1172 ppnm.push_back(ptrk);
1173 } //minus charge with with PID
1174 }
1175 int npnp = ipnp.size();
1176 int npnm = ipnm.size();
1177
1178 m_pnp = npnp;
1179 m_pnm = npnm;
1180
1181 int iphoton = 0;
1182 for (int i=0; i<m_nggneu; i++)
1183 {
1184 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGgam[i];
1185 if (!(*itTrk)->isEmcShowerValid()) continue;
1186 RecEmcShower *emcTrk = (*itTrk)->emcShower();
1187 m_numHits[iphoton] = emcTrk->numHits();
1188 m_secondmoment[iphoton] = emcTrk->secondMoment();
1189 m_x[iphoton] = emcTrk->x();
1190 m_y[iphoton]= emcTrk->y();
1191 m_z[iphoton]= emcTrk->z();
1192 m_cosemc[iphoton] = cos(emcTrk->theta());
1193 m_phiemc[iphoton] = emcTrk->phi();
1194 m_energy[iphoton] = emcTrk->energy();
1195 m_eSeed[iphoton] = emcTrk->eSeed();
1196 m_e3x3[iphoton] = emcTrk->e3x3();
1197 m_e5x5[iphoton] = emcTrk->e5x5();
1198 m_lat[iphoton] = emcTrk->latMoment();
1199 m_a20[iphoton] = emcTrk->a20Moment();
1200 m_a42[iphoton] = emcTrk->a42Moment();
1201 iphoton++;
1202 }
1203 m_tuple4->write();
1204 Ncut4++;
1205
1206 if(kmfit->chisq()>=chi5) return SUCCESS;
1207 if(pgam2pi2.m()<=1.0) return SUCCESS;
1208 if(pgam2pi1.m()<=1.0) return SUCCESS;
1209 if(nGam!=2) return SUCCESS;
1210 if(emcTg1/extmot1>=0.8) return SUCCESS;
1211 if(emcTg2/extmot2>=0.8) return SUCCESS;
1212 Ncut5++;
1213
1214 // DQA
1215 TH1 *h(0);
1216 if (m_thsvc->getHist("/DQAHist/Rhopi/mpi0", h).isSuccess()) {
1217 h->Fill(ppi0.m());
1218 } else {
1219 log << MSG::ERROR << "Couldn't retrieve mpi0" << endreq;
1220 }
1221
1222 if(fabs(ppi0.m()-0.135)>=0.015) return SUCCESS;
1223 Ncut6++;
1224
1225 if (m_thsvc->getHist("/DQAHist/Rhopi/mrho0", h).isSuccess()) {
1226 h->Fill(prho0.m());
1227 } else {
1228 log << MSG::ERROR << "Couldn't retrieve mrho0" << endreq;
1229 }
1230
1231
1232 if (m_thsvc->getHist("/DQAHist/Rhopi/mrhop", h).isSuccess()) {
1233 h->Fill(prhop.m());
1234 } else {
1235 log << MSG::ERROR << "Couldn't retrieve mrhop" << endreq;
1236 }
1237
1238 if (m_thsvc->getHist("/DQAHist/Rhopi/mrhom", h).isSuccess()) {
1239 h->Fill(prhom.m());
1240 } else {
1241 log << MSG::ERROR << "Couldn't retrieve mrhom" << endreq;
1242 }
1243
1244
1245 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
1246
1247 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(3);
1248 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(3);
1249
1250 // Quality: defined by whether dE/dx or TOF is used to identify particle
1251 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
1252 // 1 - only dE/dx (can be used for TOF calibration)
1253 // 2 - only TOF (can be used for dE/dx calibration)
1254 // 3 - Both dE/dx and TOF
1255
1256 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0);
1257 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0);
1258
1259 setFilterPassed(true);
1260
1261 return StatusCode::SUCCESS;
1262}
1263
1264
1265// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1266StatusCode DQARhopi::finalize() {
1267 cout<<"total number: "<<Ncut0<<endl;
1268 cout<<"nGood==2, nCharge==0: "<<Ncut1<<endl;
1269 cout<<"nGam>=2: "<<Ncut2<<endl;
1270 cout<<"Pass Pid: "<<Ncut3<<endl;
1271 cout<<"Pass 4C: "<<Ncut4<<endl;
1272 cout<<"final cut without pi0:"<<Ncut5<<endl;
1273 cout<<"final cut with pi0: "<<Ncut6<<endl;
1274 MsgStream log(msgSvc(), name());
1275 log << MSG::INFO << "in finalize()" << endmsg;
1276 return StatusCode::SUCCESS;
1277}
1278
u_cms
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
mpi0
const double mpi0
Definition: DQAJpsi2PPbarAlg.cxx:52
Ncut2
int Ncut2
Definition: DQARhopi.cxx:59
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: DQARhopi.cxx:36
Ncut1
int Ncut1
Definition: DQARhopi.cxx:59
u_cms
const Hep3Vector u_cms
Definition: DQARhopi.cxx:57
Vp4
std::vector< HepLorentzVector > Vp4
Definition: DQARhopi.cxx:54
Ncut0
int Ncut0
Definition: DQARhopi.cxx:59
xmass
const double xmass[5]
Definition: DQARhopi.cxx:50
Ncut5
int Ncut5
Definition: DQARhopi.cxx:59
Ncut10
int Ncut10
Definition: DQARhopi.cxx:59
velc
const double velc
Definition: DQARhopi.cxx:52
mk
const double mk
Definition: DQARhopi.cxx:49
Ncut9
int Ncut9
Definition: DQARhopi.cxx:59
Ncut3
int Ncut3
Definition: DQARhopi.cxx:59
mpi
const double mpi
Definition: DQARhopi.cxx:48
Vint
std::vector< int > Vint
Definition: DQARhopi.cxx:53
Ncut8
int Ncut8
Definition: DQARhopi.cxx:59
Ncut6
int Ncut6
Definition: DQARhopi.cxx:59
Ncut7
int Ncut7
Definition: DQARhopi.cxx:59
Ncut4
int Ncut4
Definition: DQARhopi.cxx:59
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
Ncut1
int Ncut1
Definition: Ppjrhopi.cxx:53
Ncut0
int Ncut0
Definition: Ppjrhopi.cxx:53
Ncut5
int Ncut5
Definition: Ppjrhopi.cxx:53
Ncut3
int Ncut3
Definition: Ppjrhopi.cxx:53
Ncut6
int Ncut6
Definition: Ppjrhopi.cxx:53
Ncut4
int Ncut4
Definition: Ppjrhopi.cxx:53
DQARhopi::initialize
StatusCode initialize()
Definition: DQARhopi.cxx:79
DQARhopi::execute
StatusCode execute()
Definition: DQARhopi.cxx:304
DQARhopi::DQARhopi
DQARhopi(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DQARhopi.cxx:63
DQARhopi::finalize
StatusCode finalize()
Definition: DQARhopi.cxx:1266
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
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::instance
static ParticleID * instance()
Definition: ParticleID.cxx:22
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition: ParticleID.cxx:97
ParticleID::init
void init()
Definition: ParticleID.cxx:27
ParticleID::chiDedx
double chiDedx(int n) const
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
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