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"
50#include "GaudiKernel/INTupleSvc.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 &idx_digi = 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 idx_digi = 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)[idx_digi]->GetEnergyDeposit() + lvd_E_average;
559 unsigned int ident =
CgemID::getIntID(lvi_ID_layer,lvi_ID_sheet,jj,lvi_ID_strip);
563 (*m_digi_collection)[idx_digi]->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 &idx_digi=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 idx_digi = m_digi_collection->entries() - 1;
736 ident_tmp.push_back(ident);
742 lvd_E_strip += (*m_digi_collection)[idx_digi]->GetEnergyDeposit();
749 (*m_digi_collection)[idx_digi]->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;
788 m_cgemHitTimePnt[i][j][k][l] = NULL;
795 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
798 G4int lvi_ID_HC = -1;
799 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID(
"BesCgemHitsCollection");
809 log<< MSG::ERROR <<
"CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
820 G4int lvi_N_hit = lv_HC->entries();
821 for (G4int ii=0; ii<lvi_N_hit; ii++)
827 G4int lvi_RdtElectron = lv_hit->
GetRdtEl ();
828 if(lvi_RdtElectron ==1)
continue;
849 cout <<
"..............................................................................."
851 cout<<
"CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
852 <<
" of " << lvi_N_hit <<
", E="<<lvd_E_deposit<<
", layer "<<lvi_ID_layer<<
", nsheet="<<lvi_N_sheet<< endl;
853 cout<<
" lvi_pdg "<<lvi_pdg<<
", from "<<lv3_XYZ_pre<<
" to "<<lv3_XYZ_post<<endl;
867 int particle_temp = -1;
869 if(lvi_pdg == 11){particle_temp=0;charge_temp=-1;}
870 if(lvi_pdg == -11){particle_temp=0;charge_temp=1;}
871 if(lvi_pdg == -13){particle_temp=1;charge_temp=1;}
872 if(lvi_pdg == 13){particle_temp=1;charge_temp=-1;}
873 if(lvi_pdg == 211){particle_temp=2;charge_temp=1;}
874 if(lvi_pdg == -211){particle_temp=2;charge_temp=-1;}
875 if(lvi_pdg == 321){particle_temp=3;charge_temp=1;}
876 if(lvi_pdg == -321){particle_temp=3;charge_temp=-1;}
877 if(lvi_pdg == 2212){particle_temp=4;charge_temp=1;}
878 if(lvi_pdg == -2212){particle_temp=4;charge_temp=-1;}
879 if(lvi_pdg == 1011){particle_temp=5;charge_temp=1;}
880 if(lvi_pdg == 1022){particle_temp=6;charge_temp=2;}
882 if(particle_temp==-1) {
883 cout<<
"BesCgemDigitizer::Digitize_v3(): no digitization for pdg "<<lvi_pdg<<endl;
889 posIn[0] = lv3_XYZ_pre.x();
890 posIn[1] = lv3_XYZ_pre.y();
891 posIn[2] = lv3_XYZ_pre.z();
892 posOut[0] = lv3_XYZ_post.x();
893 posOut[1] = lv3_XYZ_post.y();
894 posOut[2] = lv3_XYZ_post.z();
899 if(printFlag) { std::cout <<
"just before m_cgemDigiSvc->setTrack "<<endl;}
900 m_cgemDigiSvc->
setTrack(lvi_ID_layer, particle_temp, charge_temp, lvi_momentum.mag()/1000., posIn, posOut);
902 std::cout <<
"layer "<<lvi_ID_layer <<
", particle " << particle_temp <<
", Q= " << charge_temp <<
", p= " << lvi_momentum.mag()/1000.
903 <<
" GeV/c, (" <<posIn[0] <<
", " <<posIn[1]<<
", " <<posIn[2]<<
") to (" <<posOut[0] <<
", " <<posOut[1]<<
", " <<posOut[2] <<
")"<< std::endl;
904 std::cout <<
"R_in, R_out = "<< sqrt(posIn[0]*posIn[0]+posIn[1]*posIn[1]) <<
" " << sqrt(posOut[0]*posOut[0]+posOut[1]*posOut[1]) << std::endl;
905 std::cout <<
"check--> NXStrips, NVStrips = " << m_cgemDigiSvc->
getNXstrips() << m_cgemDigiSvc->
getNVstrips()<< std::endl;
908 vector<int> vecID[2][2];
909 vector<double> vecQ[2][2];
910 vector<double> vecT[2][2];
929 vecID[sheetID][0].push_back(m_cgemDigiSvc->
getXstripID(i));
930 vecQ[sheetID][0].push_back(m_cgemDigiSvc->
getXstripQ(i));
931 vecT[sheetID][0].push_back(m_cgemDigiSvc->
getXstripT(i));
941 vecID[sheetID][1].push_back(m_cgemDigiSvc->
getVstripID(i));
942 vecQ[sheetID][1].push_back(m_cgemDigiSvc->
getVstripQ(i));
943 vecT[sheetID][1].push_back(m_cgemDigiSvc->
getVstripT(i));
946 for(
int i=0; i<lvi_N_sheet; i++)
948 Ntot[0]+=vecID[i][0].size();
949 Ntot[1]+=vecID[i][1].size();
953 cout <<
" fired "<<Ntot[0]<<
" X strips, "<<Ntot[1]<<
" V strips" << endl;
957 vector<int> ident_tmp;
960 double Q_saturation= 45.0;
961 for(
int i=0; i<lvi_N_sheet; i++)
967 for(
unsigned int iStrip=0; iStrip<vecID[i][j].size(); iStrip++)
969 int ID=vecID[i][j][iStrip];
970 G4int &idx_digi=m_F_hit[lvi_ID_layer][i][j][ID];
1020 double EnergyDeposit_temp=vecQ[i][j][iStrip];
1022 if(EnergyDeposit_temp > Q_saturation) EnergyDeposit_temp = Q_saturation;
1025 m_digi_collection->insert(lv_new_digi);
1026 idx_digi = m_digi_collection->entries() - 1;
1031 m_cgemHitTimePnt[lvi_ID_layer][i][j][ID]=lv_hit;
1035 ident_tmp.push_back(ident);
1091 double EnergyDeposit_temp = (*m_digi_collection)[idx_digi]->GetEnergyDeposit()+vecQ[i][j][iStrip];
1093 if(EnergyDeposit_temp > Q_saturation) EnergyDeposit_temp = Q_saturation;
1094 (*m_digi_collection)[idx_digi]->SetEnergyDeposit(EnergyDeposit_temp);
1097 double GlobalTime_temp = (*m_digi_collection)[idx_digi]->GetGlobalTime();
1098 if(vecT[i][j][iStrip]<GlobalTime_temp)
1100 (*m_digi_collection)[idx_digi]->SetGlobalTime(vecT[i][j][iStrip]);
1102 BesCgemHit* lastTCgemHit = m_cgemHitTimePnt[lvi_ID_layer][i][j][ID];
1104 m_cgemHitTimePnt[lvi_ID_layer][i][j][ID]=lv_hit;
1123 int array_tmp[2000];
1124 int size_tmp=ident_tmp.size();
1125 if(size_tmp>2000) size_tmp=2000;
1126 for(
unsigned int hh =0; hh < size_tmp; hh++){
1127 array_tmp[hh] = ident_tmp[hh];
1135 StoreDigiCollection(m_digi_collection);
1139 for(
int i=0; i<m_digi_collection->entries(); i++){
1140 std::cout <<
"check CGEM digitization in simulation --> layer:" << (*m_digi_collection)[i]->GetLayerID()
1141 <<
" sheet:" << (*m_digi_collection)[i]->GetSheetID()
1142 <<
" x/v:" << (*m_digi_collection)[i]->GetStripType()
1143 <<
" ID:" << (*m_digi_collection)[i]->GetStripID()
1144 <<
" Q:" << (*m_digi_collection)[i]->GetEnergyDeposit()
1145 <<
" T:" << (*m_digi_collection)[i]->GetGlobalTime() << std::endl;
1160 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
1161 MsgStream log(
msgSvc,
"BesCgemDigitizer::DoLorentzDiffusion()");
1166 G4double lvd_A_phi_pre = f_XYZ_pre.phi();
1168 G4double lvd_A_phi_post = f_XYZ_post.phi();
1171 G4double lvd_XYZ_x_pre = f_XYZ_pre.
x() ;
1172 G4double lvd_XYZ_x_post = f_XYZ_post.
x() ;
1173 G4double lvd_XYZ_z_pre = f_XYZ_pre. z() ;
1174 G4double lvd_XYZ_z_post = f_XYZ_post. z() ;
1179 log<< MSG::INFO <<
"Position before Lorentz Diffuse: " << endreq;
1180 log<< MSG::INFO <<
"PreX=" << f_XYZ_pre.x()
1181 <<
" PreY=" << f_XYZ_pre.y()
1182 <<
" PreZ=" << f_XYZ_pre.z()
1183 <<
" PostX=" << f_XYZ_post.x()
1184 <<
" PostY=" << f_XYZ_post.y()
1185 <<
" PostZ=" << f_XYZ_post.z()
1187 log<< MSG::INFO <<
"" << endreq;
1190 f_XYZ_pre.setX ( lvd_R_layer *
cos ( lvd_A_phi_pre ));
1191 f_XYZ_pre.setY ( lvd_R_layer *
sin ( lvd_A_phi_pre ));
1194 f_XYZ_post.setX ( lvd_R_layer *
cos ( lvd_A_phi_post ));
1195 f_XYZ_post.setY ( lvd_R_layer *
sin ( lvd_A_phi_post ));
1201 log<< MSG::INFO <<
"Position after Lorentz Diffuse: " << endreq;
1202 log<< MSG::INFO <<
"PreX=" << f_XYZ_pre.x()
1203 <<
" PreY=" << f_XYZ_pre.y()
1204 <<
" PreZ=" << f_XYZ_pre.z()
1205 <<
" PostX=" << f_XYZ_post.x()
1206 <<
" PostY=" << f_XYZ_post.y()
1207 <<
" PostZ=" << f_XYZ_post.z()
1209 log<< MSG::INFO <<
"" << endreq;
1211 log<< MSG::INFO <<
"Postition Delta before and after Lorentz Diffuse: " << endreq;
1212 log<< MSG::INFO <<
"DXold=" <<
abs(lvd_XYZ_x_pre - lvd_XYZ_x_post)
1213 <<
" DXnew=" <<
abs(f_XYZ_pre.x() - f_XYZ_post.x())
1214 <<
" DZold=" <<
abs(lvd_XYZ_z_pre - lvd_XYZ_z_post)
1215 <<
" DZnew=" <<
abs(f_XYZ_pre.z() - f_XYZ_post.z())
1217 log<< MSG::INFO <<
"" << endreq;
1234 G4double &f_x, G4double &f_v, G4int &f_ID_sheet, G4int &f_ID_x, G4int &f_ID_v)
1237 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
1238 MsgStream log(
msgSvc,
"BesCgemDigitizer::GetIDFromXYZ()");
1249 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1252 G4double lvd_A_phi = f_XYZ.phi();
1257 if (lvd_A_phi >= twopi)
1263 G4double lvd_x = f_XYZ.rho() * lvd_A_phi;
1264 G4double lvd_v = f_XYZ.z() + lvd_L_layer/2.;
1268 if (lvd_x < 0 or lvd_x > lvi_N_sheet*lvd_W_sheet)
1270 log<< MSG::ERROR <<
"ERROR : CgemSim::BesCgemDigitizer::GetIDFromXYZ(), X overlaps range!" << endreq;
1271 log<< MSG::ERROR <<
"The overlaps X=" << lvd_x
1272 <<
" Its range should be < " << lvi_N_sheet*lvd_W_sheet
1273 <<
" at layerID=" << f_ID_layer
1279 for (G4int i=0; i<lvi_N_sheet; i++)
1281 if (lvd_x >= i*lvd_W_sheet and lvd_x < (i+1)*lvd_W_sheet)
1289 f_x = lvd_x - f_ID_sheet * lvd_W_sheet;
1291 f_ID_x = floor( f_x / lvd_W_pitch_x );
1292 f_ID_v = floor((f_x*
cos(lvd_A_stero) + f_v*
sin(lvd_A_stero)) / lvd_W_pitch_v);
1298 const G4int f_ID_x_post,
const G4int f_ID_v_post,
1299 G4int &f_N_hit_strip_x , G4int &f_N_hit_strip_v ,
1300 G4int &f_ID_x_start , G4int &f_ID_v_start )
1303 f_N_hit_strip_x =
abs(f_ID_x_post - f_ID_x_pre) + 1;
1304 f_N_hit_strip_v =
abs(f_ID_v_post - f_ID_v_pre) + 1;
1307 f_ID_x_start = (f_ID_x_pre <= f_ID_x_post) ? f_ID_x_pre : f_ID_x_post;
1308 f_ID_v_start = (f_ID_v_pre <= f_ID_v_post) ? f_ID_v_pre : f_ID_v_post;
1314 const G4int f_ID_sheet_pre ,
const G4double f_x_pre ,
const G4double f_v_pre ,
1315 const G4int f_ID_sheet_post,
const G4double f_x_post,
const G4double f_v_post,
1316 G4int &f_ID_sheet_mid_pre , G4int &f_ID_x_mid_pre , G4int &f_ID_v_mid_pre,
1317 G4int &f_ID_sheet_mid_post , G4int &f_ID_x_mid_post , G4int &f_ID_v_mid_post)
1326 G4double lvd_v_mid = 0;
1328 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1330 G4int lvi_N_strip_x = ceil(lvd_W_sheet / lvd_W_pitch_x);
1333 f_ID_sheet_mid_pre = f_ID_sheet_pre;
1334 f_ID_sheet_mid_post = f_ID_sheet_post;
1337 if (f_x_pre > f_x_post)
1340 f_ID_x_mid_pre = lvi_N_strip_x - 1;
1341 f_ID_x_mid_post = 0;
1344 lvd_v_mid = f_x_post*(f_v_pre-f_v_post)/(lvd_W_sheet-f_x_pre+f_x_post) + f_v_post;
1346 f_ID_v_mid_pre = floor((lvd_W_sheet*
cos(lvd_A_stero)+lvd_v_mid*
sin(lvd_A_stero)) / lvd_W_pitch_v);
1347 f_ID_v_mid_post = floor(lvd_v_mid*
sin(lvd_A_stero) / lvd_W_pitch_v);
1353 f_ID_x_mid_post = lvi_N_strip_x - 1;
1356 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;
1358 f_ID_v_mid_pre = floor(lvd_v_mid*
sin(lvd_A_stero) / lvd_W_pitch_v);
1359 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 sin(const BesAngle a)
double cos(const BesAngle a)
G4TDigiCollection< BesCgemDigi > BesCgemDigisCollection
G4THitsCollection< BesCgemHit > BesCgemHitsCollection
double abs(const EvtComplex &c)
NTuple::Array< double > m_layer
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)
void AddDigiIdxQ(G4int id)
G4ThreeVector GetPositionOfPrePoint() const
void AddIdentifier(G4int f_ID_Identifier[2000], G4int N_dim)
G4double GetGlobalTime() const
void SetDigiIdxT(G4int id)
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
double getQsaturation(int layer, int sheet, int view, int id)
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)