37#include "BesCgemDigitizer.hh"
38#include "BesCgemHit.hh"
39#include "BesCgemDigi.hh"
41#include "GaudiKernel/SmartIF.h"
42#include "GaudiKernel/ISvcLocator.h"
43#include "GaudiKernel/Bootstrap.h"
44#include "GaudiKernel/IDataProviderSvc.h"
45#include "GaudiKernel/SmartDataPtr.h"
46#include "GaudiKernel/IJobOptionsSvc.h"
47#include "GaudiKernel/PropertyMgr.h"
48#include "EventModel/EventHeader.h"
50#include "GaudiKernel/INTupleSvc.h"
51#include "Identifier/CgemID.h"
52#include "Identifier/Identifier.h"
55#include "G4DigiManager.hh"
56#include "Randomize.hh"
72:G4VDigitizerModule(modName)
75 collectionName.push_back(
"BesCgemDigisCollection");
77 ISvcLocator* svcLocator = Gaudi::svcLocator();
79 StatusCode sc=svcLocator->service(
"CgemGeomSvc", ISvc);
83 ISvcLocator* svcLocator2 = Gaudi::svcLocator();
85 StatusCode sc2=svcLocator2->service(
"CgemDigitizerSvc", ISvc2);
88 m_E_threshold = 0.*eV;
95 m_F_printHitStrip = 0;
100 static IJobOptionsSvc* jobSvc = 0;
101 sc = svcLocator->service(
"JobOptionsSvc", jobSvc);
102 if ( sc.isFailure() ) {
103 std::cout <<
"BesCgemDigitizer: Can't get the JobOptionsSvc " << std::endl;
108 const vector<const Property*>* properties = jobSvc->getProperties(
"BesSim");
109 if ( properties != NULL ) {
110 for (
unsigned int i = 0; i < properties->size(); ++i ) {
111 if ( properties->at(i)->name() ==
"CgemDigitizer" ) {
112 string strRnd = properties->at(i)->toString();
113 sscanf(strRnd.c_str(),
"%i", &m_DigitizerVer);
114 cout<<
"BesCgemDigitizer: DigitizerVer is set to "<<m_DigitizerVer<<
" by jobOptions"<<endl;
121 ISvcLocator* iface = Gaudi::svcLocator();
124 SmartIF<INTupleSvc>
ntupleSvc(iface->service(
"NTupleSvc"));
125 SmartIF<IDataProviderSvc> eds(iface->service(
"EventDataSvc"));
128 std::cout <<
"###################" << __LINE__ <<
" ntupleSvc is not valid"<< std::endl;
132 std::cout <<
"###################" << __LINE__ << std::endl;
133 m_nt1=
ntupleSvc->book(
"FILE101/hit",CLID_ColumnWiseTuple,
"hit");
136 m_nt1->addItem(
"nevt",m_evt);
137 m_nt1->addItem(
"nhit",m_nhit,0,1000000);
138 m_nt1->addIndexedItem(
"layer",m_nhit,m_layer);
139 m_nt1->addIndexedItem(
"sheet",m_nhit,m_sheet);
140 m_nt1->addIndexedItem(
"phi",m_nhit,m_phi);
141 m_nt1->addIndexedItem(
"v",m_nhit,m_v);
144 std::cout <<
"###################" << __LINE__ << std::endl;
147 std::cout <<
"###################" << __LINE__ << std::endl;
164 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
165 MsgStream log(
msgSvc,
"BesCgemDigitizer::Digitize()");
166 bool printFlag=
false;
167 ISvcLocator* iface = Gaudi::svcLocator();
168 SmartIF<IDataProviderSvc> eds(iface->service(
"EventDataSvc"));
169 SmartDataPtr<Event::EventHeader> eventHeader(eds,
"/Event/EventHeader");
170 if(m_F_ntuple) m_evt = eventHeader->eventNumber();
172 log<< MSG::INFO <<
"|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" << endreq;
173 log<< MSG::INFO <<
"INFO : CgemSim::BesCgemDigitizer::Digitize(), Begin !!!" << endreq;
177 for (G4int i=0; i<3; i++)
179 for (G4int j=0; j<2; j++)
181 for (G4int k=0; k<2; k++)
183 for (G4int l=0; l<1500; l++)
185 m_F_hit[i][j][k][l] = -1;
192 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
195 G4int lvi_ID_HC = -1;
196 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID(
"BesCgemHitsCollection");
206 log<< MSG::ERROR <<
"ERROR : CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
214 G4int lvi_ID_track = 0;
215 G4int lvi_ID_layer = 0;
216 G4double lvd_global_time = 0.;
217 G4double lvd_E_deposit = 0.;
218 G4double lvd_E_average = 0.;
219 G4ThreeVector lv3_XYZ_pre (0., 0., 0.);
220 G4ThreeVector lv3_XYZ_post(0., 0., 0.);
223 G4int lvi_RdtElectron = 0;
228 G4int lvi_N_sheet = 0;
229 G4int lvi_N_strip[2];
230 G4double lvd_R_layer = 0.;
231 G4double lvd_L_layer = 0.;
232 G4double lvd_A_stero = 0.;
233 G4double lvd_W_sheet = 0.;
234 G4double lvd_W_pitch_x = 0.;
235 G4double lvd_W_pitch_v = 0.;
238 G4double lvd_x_pre, lvd_x_post;
239 G4double lvd_v_pre, lvd_v_post;
240 G4int lvi_ID_sheet_pre , lvi_ID_x_pre, lvi_ID_x_post;
241 G4int lvi_ID_sheet_post, lvi_ID_v_pre, lvi_ID_v_post;
242 G4int lvi_N_hit_strip_x, lvi_N_hit_strip_v;
243 G4int lvi_ID_x_start , lvi_ID_v_start ;
244 G4int lvi_ID_sheet_mid_pre , lvi_ID_x_mid_pre , lvi_ID_v_mid_pre ;
245 G4int lvi_ID_sheet_mid_post, lvi_ID_x_mid_post, lvi_ID_v_mid_post;
246 G4int lvi_N_hit_strip_x_pre , lvi_N_hit_strip_v_pre , lvi_ID_x_start_pre , lvi_ID_v_start_pre ;
247 G4int lvi_N_hit_strip_x_post, lvi_N_hit_strip_v_post, lvi_ID_x_start_post, lvi_ID_v_start_post;
248 G4int lvi_N_hit_strip[2], lvi_N_hit_strip_pre[2];
249 G4int lvi_ID_strip_start[2], lvi_ID_strip_start_pre[2], lvi_ID_strip_start_post[2];
251 G4int lvi_ID_sheet = 0;
252 G4int lvi_ID_strip = 0;
253 G4double lvd_E_strip = 0.;
262 lvi_N_hit = lv_HC->entries();
263 for (G4int ii=0; ii<lvi_N_hit; ii++)
270 log<< MSG::INFO <<
"..............................................................................."
272 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
273 <<
" of " << lvi_N_hit << endreq;
277 lvi_RdtElectron = lv_hit->
GetRdtEl ();
288 lv_cgem_layer = m_cgem_geomsvc->
getCgemLayer(lvi_ID_layer);
297 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
299 lvi_N_strip[0] = ceil(lvd_W_sheet / lvd_W_pitch_x);
300 lvi_N_strip[1] = ceil((lvd_W_sheet*
cos(lvd_A_stero)+lvd_L_layer*fabs(
sin(lvd_A_stero))) / lvd_W_pitch_v);
303 if(lvi_RdtElectron ==1)
continue;
307 if (m_F_printStrip ==1)
309 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Strip information:" << endreq;
310 log<< MSG::INFO <<
"ID_layer "
318 log<< MSG::INFO << left << setw( 9) << lvi_ID_layer
319 << left << setw( 8) << lvd_R_layer
320 << left << setw( 8) << lvi_N_sheet
321 << left << setw( 8) << lvd_W_sheet
322 << left << setw(10) << lvi_N_strip[0]
323 << left << setw(10) << lvi_N_strip[1]
324 << left << setw( 8) << lvd_A_stero
326 log<< MSG::INFO <<
"lv3_XYZ_pre, lv3_XYZ_post"<<lv3_XYZ_pre<<
", "<<lv3_XYZ_post<<endreq;
327 log<< MSG::INFO <<
" " << endreq;
331 if (m_F_lorentz == 1)
334 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Lorentz Diffuse!" << endreq;
342 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Smear!" << endreq;
350 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Add Noise!" << endreq;
357 lvd_x_pre, lvd_v_pre,
358 lvi_ID_sheet_pre, lvi_ID_x_pre, lvi_ID_v_pre);
362 lvd_x_post, lvd_v_post,
363 lvi_ID_sheet_post, lvi_ID_x_post, lvi_ID_v_post);
367 if (lvi_ID_sheet_pre == lvi_ID_sheet_post and
368 (
abs(lvi_ID_x_post - lvi_ID_x_pre) < (lvi_N_strip[0]-
abs(lvi_ID_x_post - lvi_ID_x_pre))) and
369 (
abs(lvi_ID_v_post - lvi_ID_v_pre) < (lvi_N_strip[1]-
abs(lvi_ID_v_post - lvi_ID_v_pre))))
373 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Step in same sheet!" << endreq;
376 lvi_ID_x_post , lvi_ID_v_post ,
377 lvi_N_hit_strip_x , lvi_N_hit_strip_v ,
378 lvi_ID_x_start , lvi_ID_v_start );
380 lvi_N_hit_strip[0] = lvi_N_hit_strip_x;
381 lvi_N_hit_strip[1] = lvi_N_hit_strip_v;
383 lvi_ID_strip_start[0] = lvi_ID_x_start;
384 lvi_ID_strip_start[1] = lvi_ID_v_start;
386 m_phi[ii]=(lvi_ID_x_pre+lvi_ID_x_post)/2*lvd_W_pitch_x;
387 m_v[ii] =(lvi_ID_v_pre+lvi_ID_v_post)/2*lvd_W_pitch_v;
388 m_layer[ii]=lvi_ID_layer;
389 m_sheet[ii]=lvi_ID_sheet_pre;
395 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Step in different sheet!" << endreq;
399 lvi_ID_sheet_pre , lvd_x_pre , lvd_v_pre ,
400 lvi_ID_sheet_post, lvd_x_post, lvd_v_post,
401 lvi_ID_sheet_mid_pre , lvi_ID_x_mid_pre , lvi_ID_v_mid_pre,
402 lvi_ID_sheet_mid_post , lvi_ID_x_mid_post , lvi_ID_v_mid_post);
405 lvi_ID_x_mid_pre , lvi_ID_v_mid_pre ,
406 lvi_N_hit_strip_x_pre , lvi_N_hit_strip_v_pre ,
407 lvi_ID_x_start_pre , lvi_ID_v_start_pre );
410 lvi_ID_x_mid_post , lvi_ID_v_mid_post ,
411 lvi_N_hit_strip_x_post , lvi_N_hit_strip_v_post ,
412 lvi_ID_x_start_post , lvi_ID_v_start_post );
415 m_phi[ii]=(lvi_ID_x_pre+lvi_ID_x_mid_pre)/2*0.65;
416 m_v[ii] =(lvi_ID_v_pre+lvi_ID_v_mid_pre)/2*0.65;
417 m_layer[ii]=lvi_ID_layer;
418 m_sheet[ii]=lvi_ID_sheet_mid_pre;
421 lvi_N_hit_strip[0] = lvi_N_hit_strip_x_pre + lvi_N_hit_strip_x_post;
422 lvi_N_hit_strip[1] = lvi_N_hit_strip_v_pre + lvi_N_hit_strip_v_post;
423 lvi_N_hit_strip_pre[0] = lvi_N_hit_strip_x_pre;
424 lvi_N_hit_strip_pre[1] = lvi_N_hit_strip_v_pre;
425 lvi_ID_strip_start_pre[0] = lvi_ID_x_start_pre;
426 lvi_ID_strip_start_pre[1] = lvi_ID_v_start_pre;
427 lvi_ID_strip_start_post[0] = lvi_ID_x_start_post;
428 lvi_ID_strip_start_post[1] = lvi_ID_v_start_post;
433 if (m_F_printHitStrip == 1)
435 log<< MSG::INFO <<
" " << endreq;
436 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Hit strip information:"
438 log<< MSG::INFO <<
"PreSheet "
449 log<< MSG::INFO << left << setw( 9) << lvi_ID_sheet_pre
450 << left << setw(10) << lvi_ID_sheet_post
451 << left << setw( 5) << lvi_ID_x_pre
452 << left << setw( 6) << lvi_ID_x_post
453 << left << setw( 5) << lvi_ID_v_pre
454 << left << setw( 6) << lvi_ID_v_post
455 << left << setw( 4) << lvi_N_hit_strip[0]
456 << left << setw( 4) << lvi_N_hit_strip[1]
457 << left << setw( 7) << lvi_ID_strip_start[0]
458 << left << setw( 7) << lvi_ID_strip_start[1]
460 log<< MSG::INFO <<
" " << endreq;
464 lvd_E_average = lvd_E_deposit / (lvi_N_hit_strip[0] + lvi_N_hit_strip[1]);
467 if (m_F_printDigi == 1)
469 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Digi information:" << endreq;
470 log<< MSG::INFO << left << setw( 8) <<
"TrackID "
471 << left << setw( 8) <<
"LayerID "
472 << left << setw( 8) <<
"SheetID "
473 << left << setw( 8) <<
"X-0,V-1 "
474 << left << setw( 8) <<
"StripID "
475 << left << setw(16) <<
"GlobalTime "
476 << left << setw(16) <<
"EnergyDeposit "
481 vector<int> ident_tmp;
485 for (G4int jj=0; jj<2; jj++)
487 for (G4int kk=0; kk<lvi_N_hit_strip[jj]; kk++)
492 lvi_ID_sheet = lvi_ID_sheet_pre;
493 lvi_ID_strip = lvi_ID_strip_start[jj] + kk;
497 if (kk < lvi_N_hit_strip_pre[jj])
500 lvi_ID_sheet = lvi_ID_sheet_pre;
501 lvi_ID_strip = lvi_ID_strip_start_pre[jj] + kk;
506 lvi_ID_sheet = lvi_ID_sheet_post;
507 lvi_ID_strip = lvi_ID_strip_start_post[jj] + kk - lvi_N_hit_strip_pre[jj];
512 G4int &lvi_F_hit = m_F_hit[lvi_ID_layer][lvi_ID_sheet][jj][lvi_ID_strip];
515 lvd_E_strip = lvd_E_average ;
518 if (lvd_E_strip >= m_E_threshold)
530 m_digi_collection->insert(lv_new_digi);
532 lvi_F_hit = m_digi_collection->entries() - 1;
535 unsigned int ident =
CgemID::getIntID(lvi_ID_layer,lvi_ID_sheet,jj,lvi_ID_strip);
536 ident_tmp.push_back(ident);
540 if (m_F_printDigi == 1)
542 G4cout << left << setw( 8) << lvi_ID_track
543 << left << setw( 8) << lvi_ID_layer
544 << left << setw( 8) << lvi_ID_sheet
545 << left << setw( 8) << jj
546 << left << setw( 8) << lvi_ID_strip
547 << left << setw(16) << lvd_global_time
548 << left << setw(16) << lvd_E_strip
555 lvd_E_strip = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit() + lvd_E_average;
559 unsigned int ident =
CgemID::getIntID(lvi_ID_layer,lvi_ID_sheet,jj,lvi_ID_strip);
563 (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(lvd_E_strip);
571 for(
unsigned int hh =0; hh < ident_tmp.size(); hh++){
572 array_tmp[hh] = ident_tmp[hh];
579 if(m_F_ntuple) m_nhit = lvi_N_hit;
580 StoreDigiCollection(m_digi_collection);
583 if(m_F_ntuple) m_nt1->write();
590 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
591 MsgStream log(
msgSvc,
"BesCgemDigitizer::Digitize_v2()");
595 for (G4int i=0; i<3; i++)
597 for (G4int j=0; j<2; j++)
599 for (G4int k=0; k<2; k++)
601 for (G4int l=0; l<1500; l++)
603 m_F_hit[i][j][k][l] = -1;
610 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
613 G4int lvi_ID_HC = -1;
614 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID(
"BesCgemHitsCollection");
624 log<< MSG::ERROR <<
"CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
635 G4int lvi_N_hit = lv_HC->entries();
636 for (G4int ii=0; ii<lvi_N_hit; ii++)
642 G4int lvi_RdtElectron = lv_hit->
GetRdtEl ();
643 if(lvi_RdtElectron ==1)
continue;
658 log<< MSG::INFO <<
"..............................................................................."
660 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
661 <<
" of " << lvi_N_hit <<
", E="<<lvd_E_deposit<<
", layer "<<lvi_ID_layer<<
", nsheet="<<lvi_N_sheet<< endreq;
662 log<< MSG::INFO <<
" lvi_pdg "<<lvi_pdg<<
", from "<<lv3_XYZ_pre<<
" to "<<lv3_XYZ_post<<endreq;
666 if (m_F_lorentz == 1)
669 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Lorentz Diffuse!" << endreq;
672 log<< MSG::INFO <<
" after Lorentz Diffuse: from "<<lv3_XYZ_pre<<
" to "<<lv3_XYZ_post << endreq;
679 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Smear!" << endreq;
687 log<< MSG::INFO <<
"CgemSim::BesCgemDigitizer::Digitize(), Begin to Add Noise!" << endreq;
691 vector<int> vecID[2][2];
693 for(
int i=0; i<lvi_N_sheet; i++)
696 readoutPlane->
getFiredStripID(lv3_XYZ_pre,lv3_XYZ_post,vecID[i][0],vecID[i][1]);
697 Ntot[0]+=vecID[i][0].size();
698 Ntot[1]+=vecID[i][1].size();
701 log<< MSG::INFO <<
" fired "<<Ntot[0]<<
" X strips, "<<Ntot[1]<<
" V strips" << endreq;
702 double lvd_E_average=lvd_E_deposit/(Ntot[0]+Ntot[1]);
705 vector<int> ident_tmp;
708 for(
int i=0; i<lvi_N_sheet; i++)
714 for(
unsigned int iStrip=0; iStrip<vecID[i][j].size(); iStrip++)
716 int ID=vecID[i][j][iStrip];
717 G4int &lvi_F_hit=m_F_hit[lvi_ID_layer][i][j][ID];
719 double lvd_E_strip=lvd_E_average;
731 m_digi_collection->insert(lv_new_digi);
732 lvi_F_hit = m_digi_collection->entries() - 1;
736 ident_tmp.push_back(ident);
742 lvd_E_strip += (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit();
749 (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(lvd_E_strip);
757 for(
unsigned int hh =0; hh < ident_tmp.size(); hh++){
758 array_tmp[hh] = ident_tmp[hh];
767 StoreDigiCollection(m_digi_collection);
774 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
775 MsgStream log(
msgSvc,
"BesCgemDigitizer::Digitize_v3()");
776 bool printFlag=
false;
779 for (G4int i=0; i<3; i++)
781 for (G4int j=0; j<2; j++)
783 for (G4int k=0; k<2; k++)
785 for (G4int l=0; l<1500; l++)
787 m_F_hit[i][j][k][l] = -1;
794 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
797 G4int lvi_ID_HC = -1;
798 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID(
"BesCgemHitsCollection");
808 log<< MSG::ERROR <<
"CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
819 G4int lvi_N_hit = lv_HC->entries();
820 for (G4int ii=0; ii<lvi_N_hit; ii++)
826 G4int lvi_RdtElectron = lv_hit->
GetRdtEl ();
827 if(lvi_RdtElectron ==1)
continue;
848 cout <<
"..............................................................................."
850 cout<<
"CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
851 <<
" of " << lvi_N_hit <<
", E="<<lvd_E_deposit<<
", layer "<<lvi_ID_layer<<
", nsheet="<<lvi_N_sheet<< endl;
852 cout<<
" lvi_pdg "<<lvi_pdg<<
", from "<<lv3_XYZ_pre<<
" to "<<lv3_XYZ_post<<endl;
866 int particle_temp = -1;
868 if(lvi_pdg == 11){particle_temp=0;charge_temp=-1;}
869 if(lvi_pdg == -11){particle_temp=0;charge_temp=1;}
870 if(lvi_pdg == -13){particle_temp=1;charge_temp=1;}
871 if(lvi_pdg == 13){particle_temp=1;charge_temp=-1;}
872 if(lvi_pdg == 211){particle_temp=2;charge_temp=1;}
873 if(lvi_pdg == -211){particle_temp=2;charge_temp=-1;}
874 if(lvi_pdg == 321){particle_temp=3;charge_temp=1;}
875 if(lvi_pdg == -321){particle_temp=3;charge_temp=-1;}
876 if(lvi_pdg == 2212){particle_temp=4;charge_temp=1;}
877 if(lvi_pdg == -2212){particle_temp=4;charge_temp=-1;}
881 posIn[0] = lv3_XYZ_pre.x();
882 posIn[1] = lv3_XYZ_pre.y();
883 posIn[2] = lv3_XYZ_pre.z();
884 posOut[0] = lv3_XYZ_post.x();
885 posOut[1] = lv3_XYZ_post.y();
886 posOut[2] = lv3_XYZ_post.z();
891 m_cgemDigiSvc->
setTrack(lvi_ID_layer, particle_temp, charge_temp, lvi_momentum.mag()/1000., posIn, posOut);
893 std::cout << lvi_ID_layer <<
" " << particle_temp <<
" " << charge_temp <<
" " << lvi_momentum.mag()/1000. <<
" " <<posIn[0]
894 <<
" " <<posIn[1]<<
" " <<posIn[2]<<
" " <<posOut[0] <<
" " <<posOut[1]<<
" " <<posOut[2] << std::endl;
895 std::cout << sqrt(posIn[0]*posIn[0]+posIn[1]*posIn[1]) <<
" " << sqrt(posOut[0]*posOut[0]+posOut[1]*posOut[1]) << std::endl;
899 vector<int> vecID[2][2];
900 vector<double> vecQ[2][2];
901 vector<double> vecT[2][2];
920 vecID[sheetID][0].push_back(m_cgemDigiSvc->
getXstripID(i));
921 vecQ[sheetID][0].push_back(m_cgemDigiSvc->
getXstripQ(i));
922 vecT[sheetID][0].push_back(m_cgemDigiSvc->
getXstripT(i));
932 vecID[sheetID][1].push_back(m_cgemDigiSvc->
getVstripID(i));
933 vecQ[sheetID][1].push_back(m_cgemDigiSvc->
getVstripQ(i));
934 vecT[sheetID][1].push_back(m_cgemDigiSvc->
getVstripT(i));
937 for(
int i=0; i<lvi_N_sheet; i++)
939 Ntot[0]+=vecID[i][0].size();
940 Ntot[1]+=vecID[i][1].size();
944 cout <<
" fired "<<Ntot[0]<<
" X strips, "<<Ntot[1]<<
" V strips" << endl;
948 vector<int> ident_tmp;
951 double Q_saturation= 45.0;
952 for(
int i=0; i<lvi_N_sheet; i++)
958 for(
unsigned int iStrip=0; iStrip<vecID[i][j].size(); iStrip++)
960 int ID=vecID[i][j][iStrip];
961 G4int &lvi_F_hit=m_F_hit[lvi_ID_layer][i][j][ID];
1011 double EnergyDeposit_temp=vecQ[i][j][iStrip];
1012 if(EnergyDeposit_temp > Q_saturation) EnergyDeposit_temp = Q_saturation;
1015 m_digi_collection->insert(lv_new_digi);
1016 lvi_F_hit = m_digi_collection->entries() - 1;
1020 ident_tmp.push_back(ident);
1076 double EnergyDeposit_temp = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit()+vecQ[i][j][iStrip];
1077 if(EnergyDeposit_temp > Q_saturation) EnergyDeposit_temp = Q_saturation;
1078 (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(EnergyDeposit_temp);
1081 double GlobalTime_temp = (*m_digi_collection)[lvi_F_hit]->GetGlobalTime();
1082 if(vecT[i][j][iStrip]<GlobalTime_temp) (*m_digi_collection)[lvi_F_hit]->SetGlobalTime(vecT[i][j][iStrip]);
1096 int array_tmp[2000];
1097 int size_tmp=ident_tmp.size();
1098 if(size_tmp>2000) size_tmp=2000;
1099 for(
unsigned int hh =0; hh < size_tmp; hh++){
1100 array_tmp[hh] = ident_tmp[hh];
1108 StoreDigiCollection(m_digi_collection);
1112 for(
int i=0; i<m_digi_collection->entries(); i++){
1113 std::cout <<
"check CGEM digitization in simulation --> layer:" << (*m_digi_collection)[i]->GetLayerID()
1114 <<
" sheet:" << (*m_digi_collection)[i]->GetSheetID()
1115 <<
" x/v:" << (*m_digi_collection)[i]->GetStripType()
1116 <<
" ID:" << (*m_digi_collection)[i]->GetStripID()
1117 <<
" Q:" << (*m_digi_collection)[i]->GetEnergyDeposit()
1118 <<
" T:" << (*m_digi_collection)[i]->GetGlobalTime() << std::endl;
1133 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
1134 MsgStream log(
msgSvc,
"BesCgemDigitizer::DoLorentzDiffusion()");
1139 G4double lvd_A_phi_pre = f_XYZ_pre.phi();
1141 G4double lvd_A_phi_post = f_XYZ_post.phi();
1144 G4double lvd_XYZ_x_pre = f_XYZ_pre.
x() ;
1145 G4double lvd_XYZ_x_post = f_XYZ_post.
x() ;
1146 G4double lvd_XYZ_z_pre = f_XYZ_pre. z() ;
1147 G4double lvd_XYZ_z_post = f_XYZ_post. z() ;
1152 log<< MSG::INFO <<
"Position before Lorentz Diffuse: " << endreq;
1153 log<< MSG::INFO <<
"PreX=" << f_XYZ_pre.x()
1154 <<
" PreY=" << f_XYZ_pre.y()
1155 <<
" PreZ=" << f_XYZ_pre.z()
1156 <<
" PostX=" << f_XYZ_post.x()
1157 <<
" PostY=" << f_XYZ_post.y()
1158 <<
" PostZ=" << f_XYZ_post.z()
1160 log<< MSG::INFO <<
"" << endreq;
1163 f_XYZ_pre.setX ( lvd_R_layer *
cos ( lvd_A_phi_pre ));
1164 f_XYZ_pre.setY ( lvd_R_layer *
sin ( lvd_A_phi_pre ));
1167 f_XYZ_post.setX ( lvd_R_layer *
cos ( lvd_A_phi_post ));
1168 f_XYZ_post.setY ( lvd_R_layer *
sin ( lvd_A_phi_post ));
1174 log<< MSG::INFO <<
"Position after Lorentz Diffuse: " << endreq;
1175 log<< MSG::INFO <<
"PreX=" << f_XYZ_pre.x()
1176 <<
" PreY=" << f_XYZ_pre.y()
1177 <<
" PreZ=" << f_XYZ_pre.z()
1178 <<
" PostX=" << f_XYZ_post.x()
1179 <<
" PostY=" << f_XYZ_post.y()
1180 <<
" PostZ=" << f_XYZ_post.z()
1182 log<< MSG::INFO <<
"" << endreq;
1184 log<< MSG::INFO <<
"Postition Delta before and after Lorentz Diffuse: " << endreq;
1185 log<< MSG::INFO <<
"DXold=" <<
abs(lvd_XYZ_x_pre - lvd_XYZ_x_post)
1186 <<
" DXnew=" <<
abs(f_XYZ_pre.x() - f_XYZ_post.x())
1187 <<
" DZold=" <<
abs(lvd_XYZ_z_pre - lvd_XYZ_z_post)
1188 <<
" DZnew=" <<
abs(f_XYZ_pre.z() - f_XYZ_post.z())
1190 log<< MSG::INFO <<
"" << endreq;
1207 G4double &f_x, G4double &f_v, G4int &f_ID_sheet, G4int &f_ID_x, G4int &f_ID_v)
1210 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
1211 MsgStream log(
msgSvc,
"BesCgemDigitizer::GetIDFromXYZ()");
1222 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1225 G4double lvd_A_phi = f_XYZ.phi();
1230 if (lvd_A_phi >= twopi)
1236 G4double lvd_x = f_XYZ.rho() * lvd_A_phi;
1237 G4double lvd_v = f_XYZ.z() + lvd_L_layer/2.;
1241 if (lvd_x < 0 or lvd_x > lvi_N_sheet*lvd_W_sheet)
1243 log<< MSG::ERROR <<
"ERROR : CgemSim::BesCgemDigitizer::GetIDFromXYZ(), X overlaps range!" << endreq;
1244 log<< MSG::ERROR <<
"The overlaps X=" << lvd_x
1245 <<
" Its range should be < " << lvi_N_sheet*lvd_W_sheet
1246 <<
" at layerID=" << f_ID_layer
1252 for (G4int i=0; i<lvi_N_sheet; i++)
1254 if (lvd_x >= i*lvd_W_sheet and lvd_x < (i+1)*lvd_W_sheet)
1262 f_x = lvd_x - f_ID_sheet * lvd_W_sheet;
1264 f_ID_x = floor( f_x / lvd_W_pitch_x );
1265 f_ID_v = floor((f_x*
cos(lvd_A_stero) + f_v*
sin(lvd_A_stero)) / lvd_W_pitch_v);
1271 const G4int f_ID_x_post,
const G4int f_ID_v_post,
1272 G4int &f_N_hit_strip_x , G4int &f_N_hit_strip_v ,
1273 G4int &f_ID_x_start , G4int &f_ID_v_start )
1276 f_N_hit_strip_x =
abs(f_ID_x_post - f_ID_x_pre) + 1;
1277 f_N_hit_strip_v =
abs(f_ID_v_post - f_ID_v_pre) + 1;
1280 f_ID_x_start = (f_ID_x_pre <= f_ID_x_post) ? f_ID_x_pre : f_ID_x_post;
1281 f_ID_v_start = (f_ID_v_pre <= f_ID_v_post) ? f_ID_v_pre : f_ID_v_post;
1287 const G4int f_ID_sheet_pre ,
const G4double f_x_pre ,
const G4double f_v_pre ,
1288 const G4int f_ID_sheet_post,
const G4double f_x_post,
const G4double f_v_post,
1289 G4int &f_ID_sheet_mid_pre , G4int &f_ID_x_mid_pre , G4int &f_ID_v_mid_pre,
1290 G4int &f_ID_sheet_mid_post , G4int &f_ID_x_mid_post , G4int &f_ID_v_mid_post)
1299 G4double lvd_v_mid = 0;
1301 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1303 G4int lvi_N_strip_x = ceil(lvd_W_sheet / lvd_W_pitch_x);
1306 f_ID_sheet_mid_pre = f_ID_sheet_pre;
1307 f_ID_sheet_mid_post = f_ID_sheet_post;
1310 if (f_x_pre > f_x_post)
1313 f_ID_x_mid_pre = lvi_N_strip_x - 1;
1314 f_ID_x_mid_post = 0;
1317 lvd_v_mid = f_x_post*(f_v_pre-f_v_post)/(lvd_W_sheet-f_x_pre+f_x_post) + f_v_post;
1319 f_ID_v_mid_pre = floor((lvd_W_sheet*
cos(lvd_A_stero)+lvd_v_mid*
sin(lvd_A_stero)) / lvd_W_pitch_v);
1320 f_ID_v_mid_post = floor(lvd_v_mid*
sin(lvd_A_stero) / lvd_W_pitch_v);
1326 f_ID_x_mid_post = lvi_N_strip_x - 1;
1329 lvd_v_mid = (lvd_W_sheet-f_x_post)*(f_v_pre-f_v_post)/(lvd_W_sheet-f_x_post+f_x_pre) + f_v_post;
1331 f_ID_v_mid_pre = floor(lvd_v_mid*
sin(lvd_A_stero) / lvd_W_pitch_v);
1332 f_ID_v_mid_post = floor((lvd_W_sheet*
cos(lvd_A_stero)+lvd_v_mid*
sin(lvd_A_stero)) / lvd_W_pitch_v);
double abs(const EvtComplex &c)
G4TDigiCollection< BesCgemDigi > BesCgemDigisCollection
G4THitsCollection< BesCgemHit > BesCgemHitsCollection
double sin(const BesAngle a)
double cos(const BesAngle a)
void SetLayerID(G4int f_ID_layer)
void SetGlobalTime(G4double f_global_time)
void SetStripType(G4int f_F_XV)
void SetEnergyDeposit(G4double f_E_deposit)
void SetStripID(G4int f_ID_strip)
void SetSheetID(G4int f_ID_sheet)
void SetTrackID(G4int f_ID_track)
BesCgemDigitizer(G4String modName)
void GetIDFromXYZ(const G4int f_ID_layer, const G4ThreeVector f_XYZ, G4double &f_x, G4double &f_v, G4int &f_ID_sheet, G4int &f_ID_x, G4int &f_ID_v)
void DoLorentzDiffusion(const G4int f_ID_layer, G4ThreeVector &f_XYZ_pre, G4ThreeVector &f_XYZ_post, const int f_F_print=0)
void GetIDFromSegmentInSameSheet(const G4int f_ID_x_pre, const G4int f_ID_v_pre, const G4int f_ID_x_post, const G4int f_ID_v_post, G4int &f_N_hit_strip_x, G4int &f_N_hit_strip_v, G4int &f_ID_x_start, G4int &f_ID_v_start)
void GetMiddleIDFromSegmentCrossedSheet(const G4int f_ID_layer, const G4int f_ID_sheet_pre, const G4double f_x_pre, const G4double f_v_pre, const G4int f_ID_sheet_post, const G4double f_x_post, const G4double f_v_post, G4int &f_ID_sheet_mid_pre, G4int &f_ID_x_mid_pre, G4int &f_ID_v_mid_pre, G4int &f_ID_sheet_mid_post, G4int &f_ID_x_mid_post, G4int &f_ID_v_mid_post)
G4ThreeVector GetPositionOfPrePoint() const
void AddIdentifier(G4int f_ID_Identifier[2000], G4int N_dim)
G4double GetGlobalTime() const
G4ThreeVector GetPositionOfPostPoint() const
G4double GetTotalEnergyDeposit() const
G4ThreeVector GetMomentumOfPrePoint() const
double getXstripT(int n) const
double getVstripQ(int n) const
int getVstripSheet(int n) const
int getVstripID(int n) const
int getXstripSheet(int n) const
StatusCode setTrack(int layer, int particle, int charge, double p, double trkPosIn[], double trkPosOut[])
double getVstripT(int n) const
double getXstripQ(int n) const
int getXstripID(int n) const
double getWidthOfPitchV() const
double getLengthOfCgemLayer() const
double getInnerROfAnode() const
int getNumberOfSheet() const
double getWidthOfSheet() const
double getAngleOfStereo() const
double getWidthOfPitchX() const
void getFiredStripID(G4ThreeVector pos1, G4ThreeVector pos2, vector< int > &vecXID, vector< int > &vecVID) const
CgemGeoLayer * getCgemLayer(int i) const
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
static value_type getIntID(unsigned int f_layer, unsigned int f_sheet, unsigned int f_strip_type, unsigned int f_strip)