CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
McTestAlg.cxx
Go to the documentation of this file.
1////---------------------------------------------------------------------------//
2////Description:
3////Author: Dengzy
4////Created: Mar, 2004
5////Modified:
6////Comment:
7//
8#include "McTestAlg/McTestAlg.h"
9
10#include "GaudiKernel/IDataProviderSvc.h"
11#include "GaudiKernel/ISvcLocator.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/SmartDataPtr.h"
14#include "GaudiKernel/RegistryEntry.h"
15#include "GaudiKernel/MsgStream.h"
16
17#include "CLHEP/Vector/LorentzVector.h"
18#include "CLHEP/Geometry/Point3D.h"
19
20#include "MdcRawEvent/MdcDigi.h"
21#include "TofRawEvent/TofDigi.h"
22#include "EmcRawEvent/EmcDigi.h"
23#include "MucRawEvent/MucDigi.h"
24
25#include "McTruth/McParticle.h"
26#include "McTruth/MdcMcHit.h"
27#include "McTruth/TofMcHit.h"
28#include "McTruth/EmcMcHit.h"
29#include "McTruth/MucMcHit.h"
30
31#include "Identifier/Identifier.h"
32#include "Identifier/MdcID.h"
33#include "Identifier/TofID.h"
34#include "Identifier/EmcID.h"
35#include "Identifier/MucID.h"
36
37#include "RawEvent/RawDataUtil.h"
38#include "RawEvent/DigiEvent.h"
39#include "McTruth/McEvent.h"
40#include "EventModel/EventModel.h"
41#include "EventModel/EventHeader.h"
42
43McTestAlg::McTestAlg(const std::string& name, ISvcLocator* pSvcLocator) :
44Algorithm(name, pSvcLocator)
45{
46 declareProperty("ParticleRootFlag",m_particleRootFlag=false);
47 declareProperty("MdcRootFlag",m_mdcRootFlag=false);
48 declareProperty("TofRootFlag",m_tofRootFlag=false);
49}
50
51StatusCode McTestAlg::initialize()
52{
53 MsgStream log(msgSvc(), name());
54 log << MSG::INFO << " McTestAlg initialize()" << endreq;
55
56 if(m_particleRootFlag)
57 {
58 StatusCode sc;
59 NTuplePtr ntp(ntupleSvc(), "FILE900/particle");
60 if(ntp) tupleParticle = ntp;
61 else {
62 tupleParticle = ntupleSvc()->book("FILE900/particle",CLID_ColumnWiseTuple,"McTestAlg");
63 if(tupleParticle)
64 sc = tupleParticle->addItem("me",me);
65 }
66 }
67
68 if(m_mdcRootFlag)
69 {
70 StatusCode sc;
71 NTuplePtr nt1(ntupleSvc(), "FILE901/n1"); //for Mdc McTruth
72 if(nt1) tupleMdc1 = nt1;
73 else {
74 tupleMdc1 = ntupleSvc()->book("FILE901/n1",CLID_ColumnWiseTuple,"McTestAlg");
75 if(tupleMdc1)
76 {
77 sc = tupleMdc1->addItem("truthIndex",truthMdcIndex);
78 sc = tupleMdc1->addItem("truthLayer",truthMdcLayer);
79 sc = tupleMdc1->addItem("truthWire",truthMdcWire);
80 sc = tupleMdc1->addItem("truthEdep",truthMdcEdep);
81 sc = tupleMdc1->addItem("truthDriftD",truthMdcDriftD);
82 sc = tupleMdc1->addItem("truthX",truthMdcX);
83 sc = tupleMdc1->addItem("truthY",truthMdcY);
84 sc = tupleMdc1->addItem("truthZ",truthMdcZ);
85 sc = tupleMdc1->addItem("ntruth",ntruthMdc);
86 }
87 else { // did not manage to book the N tuple....
88 log << MSG::ERROR <<"Cannot book MDC N-tuple:" << long(tupleMdc1) << endmsg;
89 return StatusCode::FAILURE;
90 }
91 }
92
93 NTuplePtr nt2(ntupleSvc(), "FILE901/n2"); //for Mdc digit
94 if(nt2) tupleMdc2 = nt2;
95 else {
96 tupleMdc2 = ntupleSvc()->book("FILE901/n2",CLID_ColumnWiseTuple,"McTestAlg");
97 if(tupleMdc2)
98 {
99 sc = tupleMdc2->addItem("layer",m_layer);
100 sc = tupleMdc2->addItem("cell",m_cell);
101 sc = tupleMdc2->addItem("ADC",m_charge);
102 sc = tupleMdc2->addItem("TDC",m_time);
103 }
104 }
105 }
106
107 if(m_tofRootFlag)
108 {
109 StatusCode sc;
110 NTuplePtr nt(ntupleSvc(), "FILE902/lr");
111 if(nt) tupleTof = nt;
112 else {
113 tupleTof=ntupleSvc()->book("FILE902/lr",CLID_ColumnWiseTuple,"McTestAlg");
114 if(tupleTof)
115 {
116 sc = tupleTof->addItem("truthIndex",truthIndex);
117 sc = tupleTof->addItem("truthPartId",truthPartId);
118 sc = tupleTof->addItem("truthLayer",truthLayer);
119 sc = tupleTof->addItem("truthScinNb",truthScinNb);
120 sc = tupleTof->addItem("truthX",truthX);
121 sc = tupleTof->addItem("truthY",truthY);
122 sc = tupleTof->addItem("truthZ",truthZ);
123 sc = tupleTof->addItem("ntruth",ntruth);
124 sc = tupleTof->addItem("tleft",tleft);
125 sc = tupleTof->addItem("tright",tright);
126 sc = tupleTof->addItem("qleft",qleft);
127 sc = tupleTof->addItem("qright",qright);
128 }
129 else { // did not manage to book the N tuple....
130 log << MSG::ERROR <<"Cannot book N-tuple:" << long(tupleTof) << endmsg;
131 return StatusCode::FAILURE;
132 }
133 }
134 }
135 return StatusCode::SUCCESS;
136}
137
138StatusCode McTestAlg::finalize()
139{
140 MsgStream log(msgSvc(), name());
141 log << MSG::INFO << "McTestAlg finalize()" << endreq;
142
143 return StatusCode::SUCCESS;
144
145}
146
147
148StatusCode McTestAlg::execute()
149{
150 //interface to event data service
151 ISvcLocator* svcLocator = Gaudi::svcLocator();
152 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
153 if (sc.isFailure())
154 std::cout<<"Could not accesss EventDataSvc!"<<std::endl;
155
156 SmartDataPtr<Event::EventHeader> eventHeader(m_evtSvc,"/Event/EventHeader");
157 if(!eventHeader)
158 std::cout<<"Could not retrieve EventHeader"<<std::endl;
159
160 int event=eventHeader->eventNumber();
161 std::cout<<"event: "<<event<<std::endl;
162
163 RetrieveMcParticle();
164 RetrieveMdc();
165 RetrieveTof();
166 RetrieveEmc();
167 RetrieveMuc();
168
169 return StatusCode::SUCCESS;
170
171}
172
173void McTestAlg::RetrieveMcParticle()
174{
175 SmartDataPtr<Event::McParticleCol> mcParticleCol(m_evtSvc,"/Event/MC/McParticleCol");
176 if(!mcParticleCol)
177 std::cout<<"Could not retrieve McParticelCol"<<std::endl;
178 else
179 {
180 int pdgcode;
181 double px,py,pz,E,mass;
182 int nflag=0;
183 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
184 for (;iter_mc != mcParticleCol->end(); iter_mc++)
185 {
186 //Event::McParticle mp = (*iter_mc)->mother();
187 pdgcode = (*iter_mc)->particleProperty();
188 if((*iter_mc)->trackIndex()<0)
189 std::cout<<"ERROR! trackIndex<0"<<std::endl;
190 px=(*iter_mc)->initialFourMomentum().x();
191 py=(*iter_mc)->initialFourMomentum().y();
192 pz=(*iter_mc)->initialFourMomentum().z();
193 E=(*iter_mc)->initialFourMomentum().t();
194 if(E*E-px*px-py*py-pz*pz>=0)
195 mass=sqrt(E*E-px*px-py*py-pz*pz);
196 else
197 mass=0;
198
199 if(m_particleRootFlag)
200 {
201 if(abs(pdgcode)==11)
202 me=mass;
203 tupleParticle->write();
204 }
205 if(abs(pdgcode)==2212||abs(pdgcode)==211)
206 nflag++;
207 }
208 if(nflag!=4)
209 std::cout<<"nflag!=4"<<std::endl;
210 }
211}
212
213void McTestAlg::MdcInit()
214{
215 truthMdcIndex = -9;
216 truthMdcLayer = -9;
217 truthMdcWire = -9;
218 truthMdcEdep = -9;
219 truthMdcDriftD = -9;
220 truthMdcX = -9;
221 truthMdcY = -9;
222 truthMdcZ = -9;
223 ntruthMdc = 0;
224
225 m_layer = -9;
226 m_cell = -9;
227 m_charge = -9;
228 m_time = -9;
229}
230
231
232void McTestAlg::RetrieveMdc()
233{
234 if(m_mdcRootFlag)
235 MdcInit();
236
237 //retrieve MDC McTruth from TDS
238 SmartDataPtr<Event::MdcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MdcMcHitCol");
239 if(!aMcHitCol)
240 std::cout<<"Could not retrieve MDC McTruth collection"<<std::endl;
241 else
242 {
243 Event::MdcMcHitCol::iterator iMcHitCol;
244 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
245 {
246 const Identifier ident = (*iMcHitCol)->identify();
247 //std::cout<<(*iMcHitCol)->getTrackIndex();
248 //std::cout<<" "<<MdcID::layer(ident);
249 //std::cout<<" "<<MdcID::wire(ident);
250 //std::cout<<" "<<(*iMcHitCol)->getDepositEnergy();
251 //std::cout<<" "<<(*iMcHitCol)->getDriftDistance();
252 //std::cout<<" "<<(*iMcHitCol)->getPositionX();
253 //std::cout<<" "<<(*iMcHitCol)->getPositionY();
254 //std::cout<<" "<<(*iMcHitCol)->getPositionZ();
255 //std::cout<<std::endl;
256
257 if(m_mdcRootFlag)
258 {
259 truthMdcIndex = (*iMcHitCol)->getTrackIndex();
260 truthMdcLayer = MdcID::layer(ident);
261 truthMdcWire = MdcID::wire(ident);
262 truthMdcEdep = (*iMcHitCol)->getDepositEnergy();
263 truthMdcDriftD = (*iMcHitCol)->getDriftDistance();
264 truthMdcX = (*iMcHitCol)->getPositionX();
265 truthMdcY = (*iMcHitCol)->getPositionY();
266 truthMdcZ = (*iMcHitCol)->getPositionZ();
267 ntruthMdc++;
268 tupleMdc1->write();
269 }
270 }
271 //std::cout<<"end of retrieve MDC McTruth collection"<<std::endl;
272 }
273
274 //retrieve MDC digits from TDS
275 SmartDataPtr<MdcDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/MdcDigiCol");
276 if(!aDigiCol)
277 std::cout<<"Could not retrieve MDC digi collection"<<std::endl;
278
279 else
280 {
281 MdcDigiCol::iterator iDigiCol;
282 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
283 {
284 const Identifier ident = (*iDigiCol)->identify();
285 //std::cout<<"layer: "<<MdcID::layer(ident);
286 //std::cout<<" cell: "<<MdcID::wire(ident);
287 //std::cout<<" charge: "<<(*iDigiCol)->getChargeChannel();
288 //std::cout<<" time: "<<(*iDigiCol)->getTimeChannel()<<std::endl;
289
290 if(m_mdcRootFlag){
291 m_layer = MdcID::layer(ident);
292 m_cell = MdcID::wire(ident);
293 m_charge = (*iDigiCol)->getChargeChannel()/1.0e6;
294 m_time = (*iDigiCol)->getTimeChannel()/1.0e5;
295 tupleMdc2->write();
296 }
297 }
298 //std::cout<<"end of retrieve MDC digi collection"<<std::endl;
299 }
300}
301
302void McTestAlg::TofInit()
303{
304 truthIndex = -9;
305 truthPartId = -9;
306 truthLayer = -9;
307 truthScinNb = -9;
308 truthX = -9;
309 truthY = -9;
310 truthZ = -9;
311 ntruth = 0;
312 tleft = -9;
313 tright = -9;
314 qleft = -9;
315 qright = -9;
316}
317
318
319void McTestAlg::RetrieveTof()
320{
321 int partId,layer,scinNb,end;
322 double charge,time;
323 partId = layer = scinNb = end = -9;
324 charge = time = -9;
325 if(m_tofRootFlag)
326 TofInit();
327
328 //retrieve TOF McTruth from TDS
329 SmartDataPtr<Event::TofMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/TofMcHitCol");
330 if(!aMcHitCol)
331 std::cout<<"Could not retrieve TOF McTruth collection"<<std::endl;
332 else
333 {
334 Event::TofMcHitCol::iterator iMcHitCol;
335 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
336 {
337 const Identifier ident = (*iMcHitCol)->identify();
338 /*std::cout<<(*iMcHitCol)->getTrackIndex();
339 std::cout<<" "<<TofID::barrel_ec(ident);
340 std::cout<<" "<<TofID::layer(ident);
341 std::cout<<" "<<TofID::phi_module(ident);
342 std::cout<<" "<<(*iMcHitCol)->getPositionX();
343 std::cout<<" "<<(*iMcHitCol)->getPositionY();
344 std::cout<<" "<<(*iMcHitCol)->getPositionZ();
345 std::cout<<" "<<(*iMcHitCol)->getPx();
346 std::cout<<" "<<(*iMcHitCol)->getPy();
347 std::cout<<" "<<(*iMcHitCol)->getPz();
348 std::cout<<" "<<(*iMcHitCol)->getTrackLength();
349 std::cout<<" "<<(*iMcHitCol)->getFlightTime();
350 std::cout<<std::endl;*/
351 if(m_tofRootFlag)
352 {
353 truthIndex = (*iMcHitCol)->getTrackIndex();
354 truthPartId = TofID::barrel_ec(ident);
355 truthLayer = TofID::layer(ident);
356 truthScinNb = TofID::phi_module(ident);
357 truthX = (*iMcHitCol)->getPositionX();
358 truthY = (*iMcHitCol)->getPositionY();
359 truthZ = (*iMcHitCol)->getPositionZ();
360 ntruth++;
361 }
362 }
363 }
364
365 //retrieve TOF digits from TDS
366 SmartDataPtr<TofDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/TofDigiCol");
367 if(!aDigiCol)
368 std::cout<<"Could not retrieve TOF digi collection"<<std::endl;
369 else
370 {
371 TofDigiCol::iterator iDigiCol;
372 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
373 {
374 const Identifier ident = (*iDigiCol)->identify();
375 //std::cout<<"partId: "<<TofID::barrel_ec(ident);
376 //std::cout<<" layer: "<<TofID::layer(ident);
377 //std::cout<<" scinNb: "<<TofID::phi_module(ident);
378 //std::cout<<" end: "<<TofID::end(ident);
379 //std::cout<<std::endl;
380 //std::cout<<" charge: "<<(*iDigiCol)->getChargeChannel();
381 //std::cout<<" time: "<<(*iDigiCol)->getTimeChannel()<<std::endl;
382 //if(TofID::barrel_ec(ident)==barrel_ec && layer == TofID::layer(ident) &&
383 // phi_module == TofID::phi_module(ident) )
384 partId=TofID::barrel_ec(ident);
385 layer=TofID::layer(ident);
386 scinNb=TofID::phi_module(ident);
387 end=TofID::end(ident);
388 charge = (*iDigiCol)->getChargeChannel()/1.0e6;
389 time = (*iDigiCol)->getTimeChannel()/1.0e6;
390 if(m_tofRootFlag)
391 {
392 if(truthPartId==partId && truthLayer==layer && truthScinNb==scinNb)
393 {
394 if(end==0) {qright = charge; tright=time;}
395 else {qleft = charge; tleft = time;}
396 //std::cout<<partId<<" "<<scinNb<<" "<<charge<<" "<<time<<std::endl;
397 }
398 else
399 std::cout<<"digi doesn't match"<<std::endl;
400 }
401 }
402 if(m_tofRootFlag)
403 {
404 if(tleft>0&&tright>0&&qleft>0&&qright>0)
405 tupleTof->write();
406 else
407 std::cout<<"no digi match MCtruth"<<std::endl;
408 }
409
410 //std::cout<<"end of retrieve TOF digits"<<std::endl;
411 }
412}
413
414void McTestAlg::RetrieveEmc()
415{
416 //retrieve EMC McTruth from TDS
417 SmartDataPtr<Event::EmcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/EmcMcHitCol");
418 if(!aMcHitCol)
419 std::cout<<"Could not retrieve EMC McTruth collection"<<std::endl;
420 else
421 {
422 Event::EmcMcHitCol::iterator iMcHitCol;
423 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
424 {
425 const Identifier ident = (*iMcHitCol)->identify();
426 //std::cout<<(*iMcHitCol)->getTrackIndex();
427 //std::cout<<" "<<EmcID::barrel_ec(ident);
428 //std::cout<<" "<<EmcID::theta_module(ident);
429 //std::cout<<" "<<EmcID::phi_module(ident);
430 //std::cout<<" "<<(*iMcHitCol)->getPositionX();
431 //std::cout<<" "<<(*iMcHitCol)->getPositionY();
432 //std::cout<<" "<<(*iMcHitCol)->getPositionZ();
433 //std::cout<<" "<<(*iMcHitCol)->getPx();
434 //std::cout<<" "<<(*iMcHitCol)->getPy();
435 //std::cout<<" "<<(*iMcHitCol)->getPz();
436 //std::cout<<" "<<(*iMcHitCol)->getDepositEnergy();
437 //std::cout<<std::endl;
438 }
439 //std::cout<<"end of retrieve EMC McTruth"<<std::endl;
440 }
441
442 //retrieve EMC digits from TDS
443 SmartDataPtr<EmcDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/EmcDigiCol");
444 if(!aDigiCol)
445 std::cout<<"Could not retrieve EMC digi collection"<<std::endl;
446
447 else
448 {
449 EmcDigiCol::iterator iDigiCol;
450 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
451 {
452 const Identifier ident = (*iDigiCol)->identify();
453 //std::cout<<"barrel_ec: "<<EmcID::barrel_ec(ident);
454 //std::cout<<" theta: "<<EmcID::theta_module(ident);
455 //std::cout<<" phi: "<<EmcID::phi_module(ident);
456 //std::cout<<" charge: "<<(*iDigiCol)->getChargeChannel();
457 //std::cout<<" time: "<<(*iDigiCol)->getTimeChannel()<<std::endl;
458 }
459 //std::cout<<"end of retrieve EMC digits"<<std::endl;
460 }
461}
462
463void McTestAlg::RetrieveMuc()
464{
465 //retrieve MUC McTruth from TDS
466 SmartDataPtr<Event::MucMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MucMcHitCol");
467 if(!aMcHitCol)
468 std::cout<<"Could not retrieve MUC McTruth collection"<<std::endl;
469 else
470 {
471 Event::MucMcHitCol::iterator iMcHitCol;
472 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
473 {
474 const Identifier ident = (*iMcHitCol)->identify();
475 //std::cout<<(*iMcHitCol)->getTrackIndex();
476 //std::cout<<" "<<MucID::part(ident);
477 //std::cout<<" "<<MucID::seg(ident);
478 //std::cout<<" "<<MucID::gap(ident);
479 //std::cout<<" "<<MucID::strip(ident);
480 //std::cout<<" "<<(*iMcHitCol)->getPositionX();
481 //std::cout<<" "<<(*iMcHitCol)->getPositionY();
482 //std::cout<<" "<<(*iMcHitCol)->getPositionZ();
483 //std::cout<<" "<<(*iMcHitCol)->getPx();
484 //std::cout<<" "<<(*iMcHitCol)->getPy();
485 //std::cout<<" "<<(*iMcHitCol)->getPz();
486 //std::cout<<std::endl;
487 }
488 //std::cout<<"end of retrieve MUC McTruth"<<std::endl;
489 }
490
491 //retrieve MUC digits from TDS
492 SmartDataPtr<MucDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/MucDigiCol");
493 if(!aDigiCol)
494 std::cout<<"Could not retrieve MUC digi collection"<<std::endl;
495
496 else
497 {
498 MucDigiCol::iterator iDigiCol;
499 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
500 {
501 const Identifier ident = (*iDigiCol)->identify();
502 //std::cout<<"Part: "<<MucID::part(ident);
503 //std::cout<<" Seg: "<<MucID::seg(ident);
504 //std::cout<<" Gap: "<<MucID::gap(ident);
505 //std::cout<<" Strip: "<<MucID::strip(ident)<<std::endl;
506 }
507 //std::cout<<"end of retrieve MUC digits"<<std::endl;
508 }
509}
510
511
512
mass
double mass
Definition CosmicGenerator.cxx:138
time
Double_t time
Definition DataBase/tau_mode.c:7
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
McTestAlg::McTestAlg
McTestAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition McTestAlg.cxx:43
McTestAlg::MdcInit
void MdcInit()
Definition McTestAlg.cxx:213
McTestAlg::RetrieveEmc
void RetrieveEmc()
Definition McTestAlg.cxx:414
McTestAlg::RetrieveMuc
void RetrieveMuc()
Definition McTestAlg.cxx:463
McTestAlg::execute
StatusCode execute()
Definition McTestAlg.cxx:148
McTestAlg::RetrieveMcParticle
void RetrieveMcParticle()
Definition McTestAlg.cxx:173
McTestAlg::RetrieveMdc
void RetrieveMdc()
Definition McTestAlg.cxx:232
McTestAlg::RetrieveTof
void RetrieveTof()
Definition McTestAlg.cxx:319
McTestAlg::TofInit
void TofInit()
Definition McTestAlg.cxx:302
McTestAlg::initialize
StatusCode initialize()
Definition McTestAlg.cxx:51
McTestAlg::finalize
StatusCode finalize()
Definition McTestAlg.cxx:138
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:54
TofID::end
static int end(const Identifier &id)
Definition TofID.cxx:129
TofID::phi_module
static int phi_module(const Identifier &id)
Definition TofID.cxx:117
TofID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition TofID.cxx:95
TofID::layer
static int layer(const Identifier &id)
Definition TofID.cxx:109