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