18#include "CgemSimCheck/CgemSimCheck.h"
19#include "GaudiKernel/MsgStream.h"
20#include "GaudiKernel/AlgFactory.h"
21#include "GaudiKernel/ISvcLocator.h"
22#include "GaudiKernel/DataSvc.h"
23#include "GaudiKernel/SmartDataPtr.h"
24#include "GaudiKernel/IDataProviderSvc.h"
25#include "GaudiKernel/IDataManagerSvc.h"
26#include "GaudiKernel/PropertyMgr.h"
27#include "GaudiKernel/Bootstrap.h"
28#include "GaudiKernel/INTupleSvc.h"
29#include "EventModel/EventHeader.h"
31#include "RawEvent/RawDataUtil.h"
32#include "RawDataProviderSvc/IRawDataProviderSvc.h"
34#include "Identifier/CgemID.h"
35#include "Identifier/Identifier.h"
36#include "CgemRawEvent/CgemDigi.h"
43#include "EvTimeEvent/RecEsTime.h"
45#include "McTruth/McParticle.h"
46#include "McTruth/CgemMcHit.h"
59 Algorithm(name, pSvcLocator)
65 declareProperty(
"PostFix", m_postfix =
"");
70 cout <<
"INFO : CgemSimCheck::CgemSimCheck(), McTruth Hit will be output to a txt file!"
74 string outputFile1 =
"output_McTruth.txt";
75 cout<<
"output_McTruth : "<<outputFile1<<endl;
78 m_output_McTruth.open(outputFile1.c_str(), ios::out );
80 if (!m_output_McTruth.good())
82 cout <<
"ERROR : CgemSimCheck::CgemSimCheck(), Fail to open McTruth output file: "
83 << outputFile1 << endl;
92 cout <<
"INFO : CgemSimCheck::CgemSimCheck(), Digi Hit will be output to a txt file!"
96 string outputFile2 =
"output_Digi.txt";
97 cout<<
"Digi file : "<<outputFile2<<endl;
99 m_output_Digi.open(outputFile2.c_str(), ios::out );
101 if (!m_output_Digi.good())
103 cout <<
"ERROR : CgemSimCheck::CgemSimCheck(), Fail to open Digi output file: "
104 << outputFile2 << endl;
110 if (m_F_McParticle == 1 )
112 cout <<
"INFO : CgemSimCheck::CgemSimCheck(), MCParticle will be output to a txt file!"
116 string outputFile3 =
"output_McParticle.txt";
117 cout<<
"Digi file : "<<outputFile3<<endl;
119 m_output_McParticle.open(outputFile3.c_str(), ios::out );
121 if (!m_output_McParticle.good())
123 cout <<
"ERROR : CgemSimCheck::CgemSimCheck(), Fail to open McParticle output file: "
124 << outputFile3 << endl;
134 m_output_McTruth.close();
135 m_output_Digi.close();
141 MsgStream log(
msgSvc(), name());
142 log << MSG::INFO <<
"in beginRun()" << endreq;
143 return StatusCode::SUCCESS;
149 MsgStream log(
msgSvc(), name());
150 log << MSG::INFO <<
"in initialize()" << endreq;
154 sc = service(
"CgemGeomSvc", ISvcCgem);
155 if(sc != StatusCode::SUCCESS)
157 log << MSG::ERROR <<
"can not use CgemGeomSvc" << endreq;
158 return StatusCode::FAILURE;
163 int postfixlen = m_postfix.length();
164 char *foldername = (
char *)malloc((postfixlen + 6) *
sizeof(char));
170 foldername = (
char *)malloc((postfixlen + 12) *
sizeof(char));
171 sprintf(foldername,
"FILE135/Digi%s", m_postfix.c_str());
173 if ( nt ) m_tuple = nt;
175 m_tuple =
ntupleSvc->book(foldername, CLID_ColumnWiseTuple,
"mdcTrack");
177 m_tuple->addItem(
"run", m_run);
178 m_tuple->addItem(
"evt", m_evt);
179 m_tuple->addItem(
"nhit", m_nHit, 0, 100000);
180 m_tuple->addItem(
"X_nstrip", m_Xnstrip,0,30000);
181 m_tuple->addItem(
"V_nstrip", m_Vnstrip,0,30000);
183 m_tuple->addIndexedItem(
"hittrackID",m_nHit, m_ID_track);
184 m_tuple->addIndexedItem(
"CgemlayerID",m_nHit, m_ID_layer);
185 m_tuple->addIndexedItem(
"Readout_sheetID", m_nHit,m_ID_sheet);
188 m_tuple->addIndexedItem(
"Read_stripID", m_nHit,m_ID_strip);
189 m_tuple->addIndexedItem(
"Deposit_Energy", m_nHit,m_E_deposit);
190 m_tuple->addIndexedItem(
"Ex", m_Xnstrip,m_E_X);
191 m_tuple->addIndexedItem(
"Ev", m_Vnstrip,m_E_V);
192 m_tuple->addIndexedItem(
"Globletime", m_nHit,m_globle_time);
198 foldername = (
char *)malloc((postfixlen + 12) *
sizeof(char));
199 sprintf(foldername,
"FILE135/McTruth%s", m_postfix.c_str());
201 if ( nt1 ) mc_tuple = nt1;
203 mc_tuple =
ntupleSvc->book(foldername, CLID_ColumnWiseTuple,
"mdcTrack");
207 mc_tuple->addItem(
"MC_nhit", m_mc_nHit, 0, 100);
209 mc_tuple->addIndexedItem(
"hittrackID",m_mc_nHit, m_mc_ID_track);
210 mc_tuple->addIndexedItem(
"CgemlayerID",m_mc_nHit, m_mc_ID_layer);
211 mc_tuple->addIndexedItem(
"Deposit_Energy", m_mc_nHit,m_mc_E_deposit);
212 mc_tuple->addIndexedItem(
"X_pre_point", m_mc_nHit, m_mc_XYZ_pre_X );
213 mc_tuple->addIndexedItem(
"Y_pre_point", m_mc_nHit, m_mc_XYZ_pre_Y );
214 mc_tuple->addIndexedItem(
"Z_pre_point", m_mc_nHit, m_mc_XYZ_pre_Z );
215 mc_tuple->addIndexedItem(
"X_post_point", m_mc_nHit, m_mc_XYZ_post_X );
216 mc_tuple->addIndexedItem(
"Y_post_point", m_mc_nHit, m_mc_XYZ_post_Y );
217 mc_tuple->addIndexedItem(
"Z_post_point", m_mc_nHit, m_mc_XYZ_post_Z );
218 mc_tuple->addIndexedItem(
"P_X_pre_point", m_mc_nHit, m_mc_P_pre_X);
219 mc_tuple->addIndexedItem(
"P_Y_pre_point", m_mc_nHit, m_mc_P_pre_Y);
220 mc_tuple->addIndexedItem(
"P_Z_pre_point", m_mc_nHit, m_mc_P_pre_Z);
221 mc_tuple->addIndexedItem(
"P_X_post_point",m_mc_nHit, m_mc_P_post_X);
222 mc_tuple->addIndexedItem(
"P_Y_post_point",m_mc_nHit, m_mc_P_post_Y);
223 mc_tuple->addIndexedItem(
"P_Z_post_point",m_mc_nHit, m_mc_P_post_Z);
228 foldername = (
char *)malloc((postfixlen + 12) *
sizeof(char));
229 sprintf(foldername,
"FILE135/McPar%s", m_postfix.c_str());
231 if ( nt2 ) mcP_tuple = nt2;
233 mcP_tuple =
ntupleSvc->book(foldername, CLID_ColumnWiseTuple,
"mdcTrack");
237 mcP_tuple->addItem(
"nParticle", Nparticle, 0, 1000);
239 mcP_tuple->addIndexedItem(
"trkindex", Nparticle, m_trkindex);
240 mcP_tuple->addIndexedItem(
"x", Nparticle, m_mcParticle_x);
241 mcP_tuple->addIndexedItem(
"y", Nparticle, m_mcParticle_y);
242 mcP_tuple->addIndexedItem(
"z", Nparticle, m_mcParticle_z);
243 mcP_tuple->addIndexedItem(
"px", Nparticle, m_mcParticle_px);
244 mcP_tuple->addIndexedItem(
"py", Nparticle, m_mcParticle_py);
245 mcP_tuple->addIndexedItem(
"pz", Nparticle, m_mcParticle_pz);
246 mcP_tuple->addIndexedItem(
"E", Nparticle, m_mcParticle_E);
247 mcP_tuple->addIndexedItem(
"phi", Nparticle, m_mcParticle_phi);
248 mcP_tuple->addIndexedItem(
"costheta", Nparticle, m_mcParticle_costheta);
249 mcP_tuple->addIndexedItem(
"theta", Nparticle, m_mcParticle_theta);
250 mcP_tuple->addIndexedItem(
"pt", Nparticle, m_mcParticle_pt);
253 return StatusCode::SUCCESS;
259 MsgStream log(
msgSvc(), name());
260 log << MSG::INFO <<
"in execute()" << endreq;
262 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
265 log << MSG::FATAL <<
"Could not find Event Header" << endreq;
266 return StatusCode::FAILURE;
270 cout <<
"Event: " << eventHeader->eventNumber() << endl;
271 int iEvt = eventHeader->eventNumber();
272 int iRun = eventHeader->runNumber();
276 cout<<
" Run : "<<m_run<<
"Event : "<<m_evt<<endl;
278 if (m_F_McTruth == 1)
280 SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
281 if (!lv_CgemMcHitCol)
283 log << MSG::WARNING <<
"Could not find event" << endreq;
284 return StatusCode::FAILURE;
287 Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
289 for(; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth)
291 m_output_McTruth << left << setw( 5) << (*iter_truth)->GetTrackID()
292 << left << setw( 3) << (*iter_truth)->GetLayerID()
293 << left << setw(12) << (*iter_truth)->GetTotalEnergyDeposit()
294 << left << setw(12) << (*iter_truth)->GetPositionXOfPrePoint()
295 << left << setw(12) << (*iter_truth)->GetPositionYOfPrePoint()
296 << left << setw(12) << (*iter_truth)->GetPositionZOfPrePoint()
297 << left << setw(12) << (*iter_truth)->GetPositionXOfPostPoint()
298 << left << setw(12) << (*iter_truth)->GetPositionYOfPostPoint()
299 << left << setw(12) << (*iter_truth)->GetPositionZOfPostPoint()
300 << left << setw(12) << (*iter_truth)->GetMomentumXOfPrePoint()
301 << left << setw(12) << (*iter_truth)->GetMomentumYOfPrePoint()
302 << left << setw(12) << (*iter_truth)->GetMomentumZOfPrePoint()
303 << left << setw(12) << (*iter_truth)->GetMomentumXOfPostPoint()
304 << left << setw(12) << (*iter_truth)->GetMomentumYOfPostPoint()
305 << left << setw(12) << (*iter_truth)->GetMomentumZOfPostPoint()
309 m_mc_ID_track[i]=(*iter_truth)->GetTrackID();
310 m_mc_ID_layer[i]=(*iter_truth)->GetLayerID();
311 m_mc_E_deposit[i]=(*iter_truth)->GetTotalEnergyDeposit();
312 m_mc_XYZ_pre_X[i] =(*iter_truth)->GetPositionXOfPrePoint() ;
313 m_mc_XYZ_pre_Y[i] =(*iter_truth)->GetPositionYOfPrePoint() ;
314 m_mc_XYZ_pre_Z[i] =(*iter_truth)->GetPositionZOfPrePoint() ;
315 m_mc_XYZ_post_X[i]=(*iter_truth)->GetPositionXOfPostPoint();
316 m_mc_XYZ_post_Y[i]=(*iter_truth)->GetPositionYOfPostPoint();
317 m_mc_XYZ_post_Z[i]=(*iter_truth)->GetPositionZOfPostPoint();
318 m_mc_P_pre_X[i]=(*iter_truth)->GetMomentumXOfPrePoint() ;
319 m_mc_P_pre_Y[i]=(*iter_truth)->GetMomentumYOfPrePoint() ;
320 m_mc_P_pre_Z[i]=(*iter_truth)->GetMomentumZOfPrePoint() ;
321 m_mc_P_post_X[i]=(*iter_truth)->GetMomentumXOfPostPoint();
322 m_mc_P_post_Y[i]=(*iter_truth)->GetMomentumYOfPostPoint();
323 m_mc_P_post_Z[i]=(*iter_truth)->GetMomentumZOfPostPoint();
328 cout<<
" mc_Hit "<<m_mc_nHit<<endl;
335 SmartDataPtr<CgemDigiCol> lv_CgemDigiCol(eventSvc(),
"/Event/Digi/CgemDigiCol");
337 if ( !lv_CgemDigiCol )
339 log << MSG::WARNING <<
"Could not find event" << endreq;
340 return StatusCode::FAILURE;
343 CgemDigiCol::iterator iter_digi = lv_CgemDigiCol->begin();
348 for( ; iter_digi != lv_CgemDigiCol->end(); ++iter_digi)
351 m_output_Digi << left << setw( 8) << (*iter_digi)->getTrackIndex()
352 << left << setw( 8) <<
CgemID::layer((*iter_digi)->identify())
353 << left << setw( 8) <<
CgemID::sheet((*iter_digi)->identify())
354 << left << setw( 8) <<
CgemID::strip((*iter_digi)->identify())
355 << left << setw(16) << (*iter_digi)->getTimeChannel()
356 << left << setw(16) << (*iter_digi)->getChargeChannel()
358 xflag=cgemid->
is_xstrip((*iter_digi)->identify());
367 m_ID_track[j]=(*iter_digi)->getTrackIndex();
371 m_globle_time[j]= (*iter_digi)->getTimeChannel();
372 m_E_deposit[j]= (*iter_digi)->getChargeChannel();
375 m_E_X[m_X]=(*iter_digi)->getChargeChannel();
379 m_E_V[m_V]=(*iter_digi)->getChargeChannel();
393 if (m_F_McParticle == 1 )
395 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),
"/Event/MC/McParticleCol");
398 cout<<
"Could not retrieve McParticelCol" << endl;
401 Event::McParticleCol::iterator iter_mcP = mcParticleCol->begin();
403 for (; iter_mcP != mcParticleCol->end(); iter_mcP++)
405 m_trkindex[ii]=(*iter_mcP)->trackIndex();
406 m_mcParticle_x[ii]=((*iter_mcP)->initialPosition()).
x();
407 m_mcParticle_y[ii]=((*iter_mcP)->initialPosition()).y();
408 m_mcParticle_z[ii]=((*iter_mcP)->initialPosition()).z();
409 m_mcParticle_px[ii]=((*iter_mcP)->initialFourMomentum()).
x();
410 m_mcParticle_py[ii]=((*iter_mcP)->initialFourMomentum()).y();
411 m_mcParticle_pz[ii]=((*iter_mcP)->initialFourMomentum()).z();
412 m_mcParticle_E[ii]=((*iter_mcP)->initialFourMomentum()).e();
413 m_mcParticle_phi[ii]=((*iter_mcP)->initialFourMomentum()).phi();
414 m_mcParticle_costheta[ii]=((*iter_mcP)->initialFourMomentum()).cosTheta();
415 m_mcParticle_theta[ii]=((*iter_mcP)->initialFourMomentum()).theta();
416 m_mcParticle_pt[ii]=(((*iter_mcP)->initialFourMomentum()).mag())*
sin(m_mcParticle_theta[ii]);
418 m_output_McParticle << left << setw( 8) << (*iter_mcP)->trackIndex()
419 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).
x()
420 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).y()
421 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).z()
422 <<left<<setw(10)<<((*iter_mcP)->initialFourMomentum()).px()
423 <<left<<setw(10)<<((*iter_mcP)->initialFourMomentum()).py()
424 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).pz()
425 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).e()
426 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).phi()
427 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).cosTheta()
428 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).theta()
429 <<left<<setw(12)<<m_mcParticle_pt[ii]
432 ofstream outoneevt(
"evtput.txt",ios::app);
434 outoneevt<<left<<setw(7)<<m_evt
437 <<left << setw( 5) << (*iter_mcP)->trackIndex()
438 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).
x()
439 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).y()
440 <<left<<setw(10)<<((*iter_mcP)->initialPosition()).z()
441 <<left<<setw(10)<<((*iter_mcP)->initialFourMomentum()).px()
442 <<left<<setw(10)<<((*iter_mcP)->initialFourMomentum()).py()
443 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).pz()
444 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).e()
445 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).phi()
446 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).cosTheta()
447 <<left<<setw(12)<<((*iter_mcP)->initialFourMomentum()).theta()
448 <<left<<setw(12)<<m_mcParticle_pt[ii]
457 return StatusCode::SUCCESS;
463 MsgStream log(
msgSvc(), name());
464 log << MSG::INFO <<
"in finalize()" << endreq;
466 return StatusCode::SUCCESS;
double sin(const BesAngle a)
static int strip(const Identifier &id)
static int sheet(const Identifier &id)
static int layer(const Identifier &id)
static bool is_xstrip(const Identifier &id)
CgemSimCheck(const std::string &name, ISvcLocator *pSvcLocator)