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