CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
BesCgemDigitizer.cc
Go to the documentation of this file.
1//-------------------------------------------------------------------------------------//
2// BOOST --- BESIII Object_Oriented Simulation Tool //
3//-------------------------------------------------------------------------------------//
4/*
5 * =====================================================================================
6 *
7 * Filename: BesCgemDigitizer.cc
8 * Description:
9 * Conventions:
10 * f_ : variable used in function parameter list
11 * gv_ : global variable
12 * lv_ : local variable used in function
13 * lvd_ : local variable double
14 * m_* : normal member of class
15 * sm_* : static member of class
16 * m_M_* : class data member, material of each layer
17 * m_N_* : class data member, number of layers (CgemLayer,GemFoil)
18 * m_R_* : class data member, radius of each (material) layer, and so on
19 * m_L_* : class data member, length of each layer or hole pitch
20 * m_T_* : class data member, thickness of each (material) layer
21 * m_A_* : class data member, angle of anode VStrip
22 * Version: CgemSim-01-00-00
23 * Created: 01/05/2014 09:49:21 PM
24 * Revision: CgemSim-00-00-01(in CMT version CgemBoss-0.0.1 written by xiuql)
25 * Compiler: gcc
26 * Author: [email protected]
27 * Organization: DG1,EPC,IHEP
28 * History:
29 * <Num> <Author> <Time> <Version> <remark>
30 * 0 juxd 20140105 01-00-00 created,
31 *
32 * =====================================================================================
33 */
34
35//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36/* Header file: BOSS */
37#include "BesCgemDigitizer.hh"
38#include "BesCgemHit.hh"
39#include "BesCgemDigi.hh"
40
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"
49
50#include "GaudiKernel/INTupleSvc.h"
51#include "Identifier/CgemID.h"
52#include "Identifier/Identifier.h"
53
54/* Header file: Geant4 */
55#include "G4DigiManager.hh"
56#include "Randomize.hh"
57#include "G4ios.hh"
58#include "globals.hh"
59
60/* Header file: C++ */
61#include <iostream>
62#include <iomanip>
63#include <vector>
64#include <string>
65#include <cmath>
66
67#include "TArrayI.h"
68using namespace std;
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72:G4VDigitizerModule(modName)
73{
74 /* Digi collection name handled by BesCgemDigitizer */
75 collectionName.push_back("BesCgemDigisCollection");
76
77 ISvcLocator* svcLocator = Gaudi::svcLocator();
78 ICgemGeomSvc* ISvc;
79 StatusCode sc=svcLocator->service("CgemGeomSvc", ISvc);
80
81 m_cgem_geomsvc=dynamic_cast<CgemGeomSvc *>(ISvc);
82
83 ISvcLocator* svcLocator2 = Gaudi::svcLocator();
84 ICgemDigitizerSvc* ISvc2;
85 StatusCode sc2=svcLocator2->service("CgemDigitizerSvc", ISvc2);
86
87 m_cgemDigiSvc=dynamic_cast<CgemDigitizerSvc *>(ISvc2);
88 m_E_threshold = 0.*eV; /* Energy threshold */
89
90 m_F_lorentz = 1; /* Flag of LorentzDiffusion(): 0-OFF; 1-ON */
91 m_F_smear = 0; /* Flag of Smear(): 0-OFF; 1-ON */
92 m_F_noise = 0; /* Flag of AddNoise(): 0-OFF; 1-ON */
93
94 m_F_printStrip = 0; /* Flag of Print Strip information: 0-OFF; 1-ON */
95 m_F_printHitStrip = 0; /* Flag of Print Hit Strip information: 0-OFF; 1-ON */
96 m_F_printDigi = 0; /* Flag of Print Digi information: 0-OFF; 1-ON */
97 m_F_ntuple = 0; /* Flag of fill ntuple: 0-OFF; 1-ON */
98
99 m_DigitizerVer = 2; /* Flag of Digitizer version: 2: Geo projection, 3: full digitization by CgemDigitizerSvc */
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;
104 //return sc;
105 }
106 else
107 {
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;
115 break;
116 }
117 }
118 }
119 }
120
121 ISvcLocator* iface = Gaudi::svcLocator();
122
123 if(m_F_ntuple){
124 SmartIF<INTupleSvc> ntupleSvc(iface->service("NTupleSvc"));
125 SmartIF<IDataProviderSvc> eds(iface->service("EventDataSvc"));
126
127 if (!ntupleSvc.isValid()) {
128 std::cout << "###################" << __LINE__ << " ntupleSvc is not valid"<< std::endl;
129 }
130
131 if (ntupleSvc){
132 std::cout << "###################" << __LINE__ << std::endl;
133 m_nt1= ntupleSvc->book("FILE101/hit",CLID_ColumnWiseTuple,"hit");
134 if( m_nt1){
135
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);
142
143 }
144 std::cout << "###################" << __LINE__ << std::endl;
145 }
146 }
147 std::cout << "###################" << __LINE__ << std::endl;
148}
149
150//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152{
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157{
158 if(m_DigitizerVer==3) Digitize_v3();
159 else if(m_DigitizerVer==2) Digitize_v2();
160}
162{
163 IMessageSvc* msgSvc;
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();
171 if(printFlag) {
172 log<< MSG::INFO << "|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" << endreq;
173 log<< MSG::INFO << "INFO : CgemSim::BesCgemDigitizer::Digitize(), Begin !!!" << endreq;
174 }
175
176 /* Flag to check over whether the strip hit before */
177 for (G4int i=0; i<3; i++) /* Layer: 3 */
178 {
179 for (G4int j=0; j<2; j++) /* Sheet: 2 */
180 {
181 for (G4int k=0; k<2; k++) /* Strip Type: 2 */
182 {
183 for (G4int l=0; l<1500; l++) /* StripID: 1419 */
184 {
185 m_F_hit[i][j][k][l] = -1;
186 }
187 }
188 }
189 } /* End of 'for (int i=0; i<3; i++)' */
190
191 /* Digi manager, singleton */
192 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
193
194 /* Retrieve Hit Collection ID */
195 G4int lvi_ID_HC = -1;
196 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID("BesCgemHitsCollection");
197
198 /* Retrieve Hit Collection by ID */
199 BesCgemHitsCollection *lv_HC = 0;
200 if (lvi_ID_HC >= 0)
201 {
202 lv_HC = (BesCgemHitsCollection*) (gv_digi_manager->GetHitsCollection(lvi_ID_HC));
203 }
204 else
205 {
206 log<< MSG::ERROR << "ERROR : CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
207 << endreq;
208 return ;
209 }
210
211 /* Local Variables used by hit and digi in the following */
212 /* Local Variables by hit */
213 G4int lvi_N_hit = 0; /* Number of hits in hit collection */
214 G4int lvi_ID_track = 0; /* Track ID of hit (track, event) */
215 G4int lvi_ID_layer = 0; /* CgemLayer ID of hit */
216 G4double lvd_global_time = 0.; /* Time since the event is created */
217 G4double lvd_E_deposit = 0.; /* Energy deposited */
218 G4double lvd_E_average = 0.; /* Energy average */
219 G4ThreeVector lv3_XYZ_pre (0., 0., 0.); /* Position of pre step point */
220 G4ThreeVector lv3_XYZ_post(0., 0., 0.); /* Position of post step point */
221 //Add by sunxh
222 //G4int lvi_ID_hit = 0;
223 G4int lvi_RdtElectron = 0;
224 //end
225
226 /* Local Variables by CgemLayer */
227 CgemGeoLayer *lv_cgem_layer; /* CgemGeoLayer of the SD belongs to */
228 G4int lvi_N_sheet = 0; /* Number of readout sheet */
229 G4int lvi_N_strip[2]; /* 0-X; 1-V */
230 G4double lvd_R_layer = 0.;
231 G4double lvd_L_layer = 0.;
232 G4double lvd_A_stero = 0.; /* Angle of readout V-X */
233 G4double lvd_W_sheet = 0.; /* Width of readout sheet */
234 G4double lvd_W_pitch_x = 0.; /* Width of readout pitch for x*/
235 G4double lvd_W_pitch_v = 0.; /* Width of readout pitch for v*/
236
237 /* Local Variables by Digi */
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];
250
251 G4int lvi_ID_sheet = 0;
252 G4int lvi_ID_strip = 0;
253 G4double lvd_E_strip = 0.;
254
255 G4int lvi_case = 0;
256 /* Create digi collection by hit collection */
257 if(lv_HC)
258 {
259 m_digi_collection = new BesCgemDigisCollection(moduleName, collectionName[0]);
260
261 /* Process hits and get digi */
262 lvi_N_hit = lv_HC->entries();
263 for (G4int ii=0; ii<lvi_N_hit; ii++)
264 {
265 /* Get information from the hit */
266 BesCgemHit *lv_hit = (*lv_HC)[ii];
267
268 if(printFlag)
269 {
270 log<< MSG::INFO << "..............................................................................."
271 << endreq;
272 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
273 << " of " << lvi_N_hit << endreq;
274 }
275
276 //Add by sunxh
277 lvi_RdtElectron = lv_hit->GetRdtEl ();
278 //lvi_ID_hit = lv_hit->GetHitID ();
279 //end
280
281 lvi_ID_track = lv_hit->GetTrackID ();
282 lvi_ID_layer = lv_hit->GetLayerID ();
283 lvd_global_time = lv_hit->GetGlobalTime ();
284 lvd_E_deposit = lv_hit->GetTotalEnergyDeposit ();
285 lv3_XYZ_pre = lv_hit->GetPositionOfPrePoint ();
286 lv3_XYZ_post = lv_hit->GetPositionOfPostPoint ();
287
288 lv_cgem_layer = m_cgem_geomsvc->getCgemLayer(lvi_ID_layer);
289 lvd_R_layer = lv_cgem_layer->getInnerROfAnode()*mm;
290 lvd_L_layer = lv_cgem_layer->getLengthOfCgemLayer()*mm;
291 lvi_N_sheet = lv_cgem_layer->getNumberOfSheet();
292 lvd_A_stero = lv_cgem_layer->getAngleOfStereo()*deg;
293 lvd_W_sheet = lv_cgem_layer->getWidthOfSheet()*mm;
294 lvd_W_pitch_x = lv_cgem_layer->getWidthOfPitchX()*mm;
295 lvd_W_pitch_v = lv_cgem_layer->getWidthOfPitchV()*mm;
296
297 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
298
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);
301
302 //Add by sunxh
303 if(lvi_RdtElectron ==1) continue;// 1: close the effect of delta electron
304 //End Add
305
306 /* Print Strip Information! */
307 if (m_F_printStrip ==1)
308 {
309 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Strip information:" << endreq;
310 log<< MSG::INFO << "ID_layer "
311 << "R "
312 << "N_sheet "
313 << "W_sheet "
314 << "N_X_strip "
315 << "N_V_strip "
316 << "A_stero "
317 << endreq;
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
325 << endreq;
326 log<< MSG::INFO <<"lv3_XYZ_pre, lv3_XYZ_post"<<lv3_XYZ_pre<<", "<<lv3_XYZ_post<<endreq;
327 log<< MSG::INFO << " " << endreq;
328 }
329
330 /* Lorentz diffusion */
331 if (m_F_lorentz == 1)
332 {
333 if(printFlag)
334 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Lorentz Diffuse!" << endreq;
335 DoLorentzDiffusion(lvi_ID_layer, lv3_XYZ_pre, lv3_XYZ_post);
336 }
337
338 /* Smear */
339 if (m_F_smear == 1)
340 {
341 if(printFlag)
342 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Smear!" << endreq;
343 Smear();
344 }
345
346 /* Add noise */
347 if (m_F_noise == 1)
348 {
349 if(printFlag)
350 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Add Noise!" << endreq;
351 AddNoise();
352 }
353
354 /* Get ID from position XYZ of readout layer */
355 /* Pre Point of the step */
356 GetIDFromXYZ(lvi_ID_layer, lv3_XYZ_pre,
357 lvd_x_pre, lvd_v_pre,
358 lvi_ID_sheet_pre, lvi_ID_x_pre, lvi_ID_v_pre);
359
360 /* Post Point of the step */
361 GetIDFromXYZ(lvi_ID_layer, lv3_XYZ_post,
362 lvd_x_post, lvd_v_post,
363 lvi_ID_sheet_post, lvi_ID_x_post, lvi_ID_v_post);
364
365 /* Suppose one step cannot cross over one total sheet */
366 /* Segment in the same sheet */
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))))
370 {
371 lvi_case = 0;
372 if(printFlag)
373 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Step in same sheet!" << endreq;
374
375 GetIDFromSegmentInSameSheet(lvi_ID_x_pre , lvi_ID_v_pre , /* Input */
376 lvi_ID_x_post , lvi_ID_v_post , /* Input */
377 lvi_N_hit_strip_x , lvi_N_hit_strip_v , /* Output */
378 lvi_ID_x_start , lvi_ID_v_start ); /* Output */
379
380 lvi_N_hit_strip[0] = lvi_N_hit_strip_x;
381 lvi_N_hit_strip[1] = lvi_N_hit_strip_v;
382
383 lvi_ID_strip_start[0] = lvi_ID_x_start;
384 lvi_ID_strip_start[1] = lvi_ID_v_start;
385
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;
390 }
391 /* Segment in different sheets, split it into two segemnts which are in same sheet */
392 else
393 {
394 if(printFlag)
395 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Step in different sheet!" << endreq;
396
397 lvi_case = 1;
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);
403
404 GetIDFromSegmentInSameSheet(lvi_ID_x_pre , lvi_ID_v_pre , /* Input */
405 lvi_ID_x_mid_pre , lvi_ID_v_mid_pre , /* Input */
406 lvi_N_hit_strip_x_pre , lvi_N_hit_strip_v_pre , /* Output */
407 lvi_ID_x_start_pre , lvi_ID_v_start_pre ); /* Output */
408
409 GetIDFromSegmentInSameSheet(lvi_ID_x_post , lvi_ID_v_post , /* Input */
410 lvi_ID_x_mid_post , lvi_ID_v_mid_post , /* Input */
411 lvi_N_hit_strip_x_post , lvi_N_hit_strip_v_post , /* Output */
412 lvi_ID_x_start_post , lvi_ID_v_start_post ); /* Output */
413
414 if(m_F_ntuple){
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;
419 }
420
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;
429 }
430
431
432 /* Print hit strip information */
433 if (m_F_printHitStrip == 1)
434 {
435 log<< MSG::INFO << " " << endreq;
436 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Hit strip information:"
437 << endreq;
438 log<< MSG::INFO << "PreSheet "
439 << "PostSheet "
440 << "PreX "
441 << "PostX "
442 << "PreV "
443 << "PostV "
444 << "N_X "
445 << "N_V "
446 << "StartX "
447 << "StartV "
448 << endreq;
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]
459 << endreq;
460 log<< MSG::INFO << " " << endreq;
461 }
462
463 /* uniform energy through sgement pre-post */
464 lvd_E_average = lvd_E_deposit / (lvi_N_hit_strip[0] + lvi_N_hit_strip[1]);
465
466 /* Print Information of Digi */
467 if (m_F_printDigi == 1)
468 {
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 "
477 << endreq;
478 }
479
480 /* Add by sunxh */
481 vector<int> ident_tmp;
482 //End Add
483
484 /* Create Digi */
485 for (G4int jj=0; jj<2; jj++) /* 0-X; 1-V */
486 {
487 for (G4int kk=0; kk<lvi_N_hit_strip[jj]; kk++)
488 {
489 if (lvi_case==0)
490 {
491 //G4cout << "### step in same sheet" << G4endl;
492 lvi_ID_sheet = lvi_ID_sheet_pre;
493 lvi_ID_strip = lvi_ID_strip_start[jj] + kk;
494 }
495 else
496 {
497 if (kk < lvi_N_hit_strip_pre[jj])
498 {
499 //G4cout << "*** step in different sheet" << G4endl;
500 lvi_ID_sheet = lvi_ID_sheet_pre;
501 lvi_ID_strip = lvi_ID_strip_start_pre[jj] + kk;
502 }
503 else
504 {
505 //G4cout << "*** step in different sheet" << G4endl;
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];
508 }
509 }
510
511 /* Compute energy of the hit strip */
512 G4int &lvi_F_hit = m_F_hit[lvi_ID_layer][lvi_ID_sheet][jj][lvi_ID_strip];
513 if (lvi_F_hit == -1)
514 {
515 lvd_E_strip = lvd_E_average ;
516
517 /* Signal should be larger than threshold */
518 if (lvd_E_strip >= m_E_threshold)
519 {
520 BesCgemDigi *lv_new_digi = new BesCgemDigi(); /* Create Digi */
521
522 lv_new_digi->SetTrackID ( lvi_ID_track );
523 lv_new_digi->SetLayerID ( lvi_ID_layer );
524 lv_new_digi->SetSheetID ( lvi_ID_sheet );
525 lv_new_digi->SetStripType ( jj );
526 lv_new_digi->SetStripID ( lvi_ID_strip );
527 lv_new_digi->SetGlobalTime ( lvd_global_time );
528 lv_new_digi->SetEnergyDeposit ( lvd_E_strip );
529
530 m_digi_collection->insert(lv_new_digi);
531
532 lvi_F_hit = m_digi_collection->entries() - 1;
533
534 //Add by sunxh
535 unsigned int ident = CgemID::getIntID(lvi_ID_layer,lvi_ID_sheet,jj,lvi_ID_strip);
536 ident_tmp.push_back(ident);
537 //End Add
538
539 /* Print Digi information! */
540 if (m_F_printDigi == 1)
541 {
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
549 << G4endl;
550 }
551 }
552 }
553 else
554 {
555 lvd_E_strip = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit() + lvd_E_average;
556
557
558 //Add by sunxh
559 unsigned int ident = CgemID::getIntID(lvi_ID_layer,lvi_ID_sheet,jj,lvi_ID_strip);
560 ident_tmp.push_back(Identifier(ident));
561 //End Add
562
563 (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(lvd_E_strip);
564 }
565
566 } /* End of 'for (G4int kk=0; kk<lvi_N_hit_strip[jj]; kk++)' */
567 } /* End of 'for (G4int jj=0; jj<2; jj++)' */
568
569 /* Add by sunxh */
570 int array_tmp[2000];
571 for(unsigned int hh =0; hh < ident_tmp.size(); hh++){
572 array_tmp[hh] = ident_tmp[hh];
573 }
574 lv_hit->AddIdentifier(array_tmp,ident_tmp.size());
575 /*End Add*/
576
577 } /* End of 'for (G4int ii=0; ii<lvi_N_hit; ii++)' */
578
579 if(m_F_ntuple) m_nhit = lvi_N_hit;
580 StoreDigiCollection(m_digi_collection); /* */
581
582 } /* End of 'if(lv_HC)' */
583 if(m_F_ntuple) m_nt1->write();
584}
585
586
588{
589 IMessageSvc* msgSvc;
590 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
591 MsgStream log(msgSvc, "BesCgemDigitizer::Digitize_v2()");
592 bool printFlag=true;
593
594 /* Flag to check over whether the strip hit before */
595 for (G4int i=0; i<3; i++) /* Layer: 3 */
596 {
597 for (G4int j=0; j<2; j++) /* Sheet: 2 */
598 {
599 for (G4int k=0; k<2; k++) /* Strip Type: 2 */
600 {
601 for (G4int l=0; l<1500; l++) /* StripID: 1419 */
602 {
603 m_F_hit[i][j][k][l] = -1;
604 }
605 }
606 }
607 } /* End of 'for (int i=0; i<3; i++)' */
608
609 /* Digi manager, singleton */
610 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
611
612 /* Retrieve Hit Collection ID */
613 G4int lvi_ID_HC = -1;
614 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID("BesCgemHitsCollection");
615
616 /* Retrieve Hit Collection by ID */
617 BesCgemHitsCollection *lv_HC = 0;
618 if (lvi_ID_HC >= 0)
619 {
620 lv_HC = (BesCgemHitsCollection*) (gv_digi_manager->GetHitsCollection(lvi_ID_HC));
621 }
622 else
623 {
624 log<< MSG::ERROR << "CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
625 << endreq;
626 return ;
627 }
628
629 /* Create digi collection by hit collection */
630 if(lv_HC)
631 {
632 m_digi_collection = new BesCgemDigisCollection(moduleName, collectionName[0]);
633
634 /* Process hits and get digi */
635 G4int lvi_N_hit = lv_HC->entries();
636 for (G4int ii=0; ii<lvi_N_hit; ii++)
637 {
638 /* Get information from the hit */
639 BesCgemHit *lv_hit = (*lv_HC)[ii];
640
641 //Add by sunxh
642 G4int lvi_RdtElectron = lv_hit->GetRdtEl ();
643 if(lvi_RdtElectron ==1) continue;// 1: close the effect of delta electron
644 //End Add
645
646 G4int lvi_ID_track = lv_hit->GetTrackID ();
647 G4int lvi_ID_layer = lv_hit->GetLayerID ();
648 G4int lvi_pdg = lv_hit->GetPDGCode ();
649 //G4int lvi_ID_parent = lv_hit->GetParentID ();
650 //if(lvi_pdg!=13) continue;
651 G4double lvd_global_time = lv_hit->GetGlobalTime ();
652 G4double lvd_E_deposit = lv_hit->GetTotalEnergyDeposit ();
653 G4ThreeVector lv3_XYZ_pre = lv_hit->GetPositionOfPrePoint ();
654 G4ThreeVector lv3_XYZ_post = lv_hit->GetPositionOfPostPoint ();
655 CgemGeoLayer *lv_cgem_layer = m_cgem_geomsvc->getCgemLayer(lvi_ID_layer);
656 G4int lvi_N_sheet = lv_cgem_layer->getNumberOfSheet();
657 if(printFlag) {
658 log<< MSG::INFO << "..............................................................................."
659 << endreq;
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;
663 }
664
665 /* Lorentz diffusion */
666 if (m_F_lorentz == 1)
667 {
668 if(printFlag)
669 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Lorentz Diffuse!" << endreq;
670 DoLorentzDiffusion(lvi_ID_layer, lv3_XYZ_pre, lv3_XYZ_post);
671 if(printFlag)
672 log<< MSG::INFO << " after Lorentz Diffuse: from "<<lv3_XYZ_pre<<" to "<<lv3_XYZ_post << endreq;
673 }// lv3_XYZ_pre, lv3_XYZ_post already on readout plane
674
675 /* Smear */
676 if (m_F_smear == 1)
677 {
678 if(printFlag)
679 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Smear!" << endreq;
680 Smear();
681 }
682
683 /* Add noise */
684 if (m_F_noise == 1)
685 {
686 if(printFlag)
687 log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Begin to Add Noise!" << endreq;
688 AddNoise();
689 }
690
691 vector<int> vecID[2][2];// [iSheet][isXV]
692 int Ntot[2]={0,0};// [isXV]
693 for(int i=0; i<lvi_N_sheet; i++)
694 {
695 CgemGeoReadoutPlane* readoutPlane=m_cgem_geomsvc->getReadoutPlane(lvi_ID_layer,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();
699 }
700 if(printFlag)
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]);
703
704 /* Add by sunxh */
705 vector<int> ident_tmp;
706 //End Add
707
708 for(int i=0; i<lvi_N_sheet; i++)
709 {
710 //cout<<"* sheet "<<i<<endl;
711 for(int j=0;j<2;j++)// isXV
712 {
713 //cout<<" * XV "<<j<<": "<<endl;
714 for(unsigned int iStrip=0; iStrip<vecID[i][j].size(); iStrip++)
715 {
716 int ID=vecID[i][j][iStrip];
717 G4int &lvi_F_hit=m_F_hit[lvi_ID_layer][i][j][ID];
718 BesCgemDigi *lv_new_digi;
719 double lvd_E_strip=lvd_E_average;
720 //cout<<" ID, E="<<ID<<", "<<lvd_E_strip<<endl;
721 if(lvi_F_hit == -1)
722 {
723 lv_new_digi = new BesCgemDigi();
724 lv_new_digi->SetTrackID ( lvi_ID_track );
725 lv_new_digi->SetLayerID ( lvi_ID_layer );
726 lv_new_digi->SetSheetID ( i);
727 lv_new_digi->SetStripType ( j );
728 lv_new_digi->SetStripID ( ID );
729 lv_new_digi->SetGlobalTime ( lvd_global_time );
730 lv_new_digi->SetEnergyDeposit ( lvd_E_strip );
731 m_digi_collection->insert(lv_new_digi);
732 lvi_F_hit = m_digi_collection->entries() - 1;
733
734 //Add by sunxh
735 unsigned int ident = CgemID::getIntID(lvi_ID_layer,i,j,ID);
736 ident_tmp.push_back(ident);
737 //End Add
738
739 }
740 else
741 {
742 lvd_E_strip += (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit();
743
744 //Add by sunxh
745 unsigned int ident = CgemID::getIntID(lvi_ID_layer,i,j,ID);
746 ident_tmp.push_back(Identifier(ident));
747 //End Add
748
749 (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(lvd_E_strip);
750 }
751 }// loop strip
752 }// loop XV
753 }// loop sheet
754
755 /* Add by sunxh */
756 int array_tmp[2000];
757 for(unsigned int hh =0; hh < ident_tmp.size(); hh++){
758 array_tmp[hh] = ident_tmp[hh];
759 }
760 lv_hit->AddIdentifier(array_tmp,ident_tmp.size());
761
762 // TArrayI tmp = lv_hit->GetIdentifier();
763
764 /*End Add*/
765
766 }// loop hits ii
767 StoreDigiCollection(m_digi_collection); /* */
768 }//if(lv_HC)
769}
770
772{
773 IMessageSvc* msgSvc;
774 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
775 MsgStream log(msgSvc, "BesCgemDigitizer::Digitize_v3()");
776 bool printFlag=false;
777
778 /* Flag to check over whether the strip hit before */
779 for (G4int i=0; i<3; i++) /* Layer: 3 */
780 {
781 for (G4int j=0; j<2; j++) /* Sheet: 2 */
782 {
783 for (G4int k=0; k<2; k++) /* Strip Type: 2 */
784 {
785 for (G4int l=0; l<1500; l++) /* max StripID: 1395=2790/2 */
786 {
787 m_F_hit[i][j][k][l] = -1;
788 }
789 }
790 }
791 } /* End of 'for (int i=0; i<3; i++)' */
792
793 /* Digi manager, singleton */
794 G4DigiManager *gv_digi_manager = G4DigiManager::GetDMpointer();
795
796 /* Retrieve Hit Collection ID */
797 G4int lvi_ID_HC = -1;
798 lvi_ID_HC = gv_digi_manager->GetHitsCollectionID("BesCgemHitsCollection");
799
800 /* Retrieve Hit Collection by ID */
801 BesCgemHitsCollection *lv_HC = 0;
802 if (lvi_ID_HC >= 0)
803 {
804 lv_HC = (BesCgemHitsCollection*) (gv_digi_manager->GetHitsCollection(lvi_ID_HC));
805 }
806 else
807 {
808 log<< MSG::ERROR << "CgemSim::BesCgemDigitizer::Digitize(), Fail to retrieve HitsCollection!"
809 << endreq;
810 return ;
811 }
812
813 /* Create digi collection by hit collection */
814 if(lv_HC)
815 {
816 m_digi_collection = new BesCgemDigisCollection(moduleName, collectionName[0]);
817
818 /* Process hits and get digi */
819 G4int lvi_N_hit = lv_HC->entries();
820 for (G4int ii=0; ii<lvi_N_hit; ii++)
821 {
822 /* Get information from the hit */
823 BesCgemHit *lv_hit = (*lv_HC)[ii];
824
825 //Add by sunxh
826 G4int lvi_RdtElectron = lv_hit->GetRdtEl ();
827 if(lvi_RdtElectron ==1) continue;// 1: close the effect of delta electron
828 //End Add
829
830 G4int lvi_ID_track = lv_hit->GetTrackID ();
831 G4int lvi_ID_layer = lv_hit->GetLayerID ();
832 G4int lvi_pdg = lv_hit->GetPDGCode ();
833 G4ThreeVector lvi_momentum = lv_hit->GetMomentumOfPrePoint();
834 //G4int lvi_ID_parent = lv_hit->GetParentID ();
835 //if(lvi_pdg!=13) continue;
836 //G4double lvd_global_time = lv_hit->GetGlobalTime ();
837 G4double lvd_E_deposit = lv_hit->GetTotalEnergyDeposit ();
838 G4ThreeVector lv3_XYZ_pre = lv_hit->GetPositionOfPrePoint ();
839 G4ThreeVector lv3_XYZ_post = lv_hit->GetPositionOfPostPoint ();
840 CgemGeoLayer *lv_cgem_layer = m_cgem_geomsvc->getCgemLayer(lvi_ID_layer);
841 G4int lvi_N_sheet = lv_cgem_layer->getNumberOfSheet();
842 if(printFlag) {
843 //log<< MSG::INFO << "..............................................................................."
844 // << endreq;
845 //log<< MSG::INFO << "CgemSim::BesCgemDigitizer::Digitize(), Digitize hit " << ii
846 // << " of " << lvi_N_hit <<", E="<<lvd_E_deposit<<", layer "<<lvi_ID_layer<<", nsheet="<<lvi_N_sheet<< endreq;
847 //log<< MSG::INFO << " lvi_pdg "<<lvi_pdg<<", from "<<lv3_XYZ_pre<<" to "<<lv3_XYZ_post<<endreq;
848 cout << "..............................................................................."
849 << endl;
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;
853 }
854
855 // if((0==particle)&&(-1==charge)) {particleName = "electron";} add p Lepton e- 11
856 // else if((0==particle)&&(1==charge)) {particleName = "e+";} add p Lepton e+ -11
857 // else if((1==particle)&&(1==charge)) {particleName = "mu+";} add p Lepton mu+ -13
858 // else if((1==particle)&&(-1==charge)) {particleName = "mu-";} add p Lepton mu- 13
859 // else if((2==particle)&&(1==charge)) {particleName = "pi+";} add p Meson pi+ 211
860 // else if((2==particle)&&(-1==charge)) {particleName = "pi-";} add p Meson pi- -211
861 // else if((3==particle)&&(1==charge)) {particleName = "K+";} add p Meson K+ 321
862 // else if((3==particle)&&(-1==charge)) {particleName = "K-";} add p Meson K- -321
863 // else if((4==particle)&&(1==charge)) {particleName = "p";} add p Baryon p+ 2212
864 // else if((4==particle)&&(-1==charge)) {particleName = "p-bar";} add p Baryon anti-p- -2212
865
866 int particle_temp = -1;
867 int charge_temp = 0;
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;}
878
879 double posIn[3];
880 double posOut[3];
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();
887
888 //if(sqrt(posIn[0]*posIn[0]+posIn[1]*posIn[1]) < sqrt(posOut[0]*posOut[0]+posOut[1]*posOut[1])){
889
890 /* Get information from digitization service */
891 m_cgemDigiSvc->setTrack(lvi_ID_layer, particle_temp, charge_temp, lvi_momentum.mag()/1000., posIn, posOut);
892 if(printFlag) {
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;
896 std::cout << "check-->" << m_cgemDigiSvc->getNXstrips() << m_cgemDigiSvc->getNVstrips()<< std::endl;
897 }
898
899 vector<int> vecID[2][2];// [iSheet][isXV]
900 vector<double> vecQ[2][2]; // [iSheet][isXV]
901 vector<double> vecT[2][2]; // [iSheet][isXV]
902 int Ntot[2]={0,0};// [isXV]
903 //for(int i=0; i<2; i++){
904 // for(int j=0; j<2; j++){
905 // vecID[i][j].clear();
906 // vecQ[i][j].clear();
907 // vecT[i][j].clear();
908 // }
909 //}
910
911 for(int i=0; i<m_cgemDigiSvc->getNXstrips(); i++){
912 //if(m_cgemDigiSvc->getXstripSheet(i)==0) vecID[0][0].push_back(m_cgemDigiSvc->getXstripID(i));
913 //if(m_cgemDigiSvc->getXstripSheet(i)==1) vecID[1][0].push_back(m_cgemDigiSvc->getXstripID(i));
914 //if(m_cgemDigiSvc->getXstripSheet(i)==0) vecQ[0][0].push_back(m_cgemDigiSvc->getXstripQ(i));
915 //if(m_cgemDigiSvc->getXstripSheet(i)==1) vecQ[1][0].push_back(m_cgemDigiSvc->getXstripQ(i));
916 //if(m_cgemDigiSvc->getXstripSheet(i)==0) vecT[0][0].push_back(m_cgemDigiSvc->getXstripT(i));
917 //if(m_cgemDigiSvc->getXstripSheet(i)==1) vecT[1][0].push_back(m_cgemDigiSvc->getXstripT(i));
918 //std::cout << "zhaojy check sheet" << m_cgemDigiSvc->getXstripSheet(i) << std::endl;
919 int sheetID=m_cgemDigiSvc->getXstripSheet(i);
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));
923 }
924 for(int i=0; i<m_cgemDigiSvc->getNVstrips(); i++){
925 //if(m_cgemDigiSvc->getVstripSheet(i)==0) vecID[0][1].push_back(m_cgemDigiSvc->getVstripID(i));
926 //if(m_cgemDigiSvc->getVstripSheet(i)==1) vecID[1][1].push_back(m_cgemDigiSvc->getVstripID(i));
927 //if(m_cgemDigiSvc->getVstripSheet(i)==0) vecQ[0][1].push_back(m_cgemDigiSvc->getVstripQ(i));
928 //if(m_cgemDigiSvc->getVstripSheet(i)==1) vecQ[1][1].push_back(m_cgemDigiSvc->getVstripQ(i));
929 //if(m_cgemDigiSvc->getVstripSheet(i)==0) vecT[0][1].push_back(m_cgemDigiSvc->getVstripT(i));
930 //if(m_cgemDigiSvc->getVstripSheet(i)==1) vecT[1][1].push_back(m_cgemDigiSvc->getVstripT(i));
931 int sheetID=m_cgemDigiSvc->getVstripSheet(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));
935 }
936
937 for(int i=0; i<lvi_N_sheet; i++)
938 {
939 Ntot[0]+=vecID[i][0].size();
940 Ntot[1]+=vecID[i][1].size();
941 }
942
943 if(printFlag)
944 cout << " fired "<<Ntot[0]<<" X strips, "<<Ntot[1]<<" V strips" << endl;
945 //double lvd_E_average=lvd_E_deposit/(Ntot[0]+Ntot[1]);
946
947 /* Add by sunxh */
948 vector<int> ident_tmp;
949 //End Add
950
951 double Q_saturation= 45.0; //-----Q Satruation
952 for(int i=0; i<lvi_N_sheet; i++)
953 {
954 //cout<<"* sheet "<<i<<endl;
955 for(int j=0;j<2;j++)// isXV
956 {
957 //cout<<" * XV "<<j<<": "<<endl;
958 for(unsigned int iStrip=0; iStrip<vecID[i][j].size(); iStrip++)
959 {
960 int ID=vecID[i][j][iStrip];
961 G4int &lvi_F_hit=m_F_hit[lvi_ID_layer][i][j][ID];
962 BesCgemDigi *lv_new_digi;
963 //double lvd_E_strip=lvd_E_average;
964 //cout<<" ID, E="<<ID<<", "<<lvd_E_strip<<endl;
965 if(lvi_F_hit == -1)
966 {
967 lv_new_digi = new BesCgemDigi();
968 lv_new_digi->SetTrackID ( lvi_ID_track );
969 lv_new_digi->SetLayerID ( lvi_ID_layer );
970 lv_new_digi->SetSheetID ( i);
971 lv_new_digi->SetStripType ( j );
972 lv_new_digi->SetStripID ( ID );
973
974 //if(j==0){ // X strip
975 // if(i==0){ // Sheet 0
976 // for(int k=0; k<vecID[0][0].size(); k++){
977 // if(vecID[0][0][k]==ID){
978 // lv_new_digi->SetGlobalTime ( vecT[0][0][k] );
979 // lv_new_digi->SetEnergyDeposit ( vecQ[0][0][k] );
980 // }
981 // }
982 // }
983 // else if(i==1){ // Sheet 1
984 // for(int k=0; k<vecID[1][0].size(); k++){
985 // if(vecID[1][0][k]==ID){
986 // lv_new_digi->SetGlobalTime ( vecT[1][0][k] );
987 // lv_new_digi->SetEnergyDeposit ( vecQ[1][0][k] );
988 // }
989 // }
990 // }
991 //}
992 //else if(j==1){ // V strip
993 // if(i==0){ // Sheet 0
994 // for(int k=0; k<vecID[0][1].size(); k++){
995 // if(vecID[0][1][k]==ID){
996 // lv_new_digi->SetGlobalTime ( vecT[0][1][k] );
997 // lv_new_digi->SetEnergyDeposit ( vecQ[0][1][k] );
998 // }
999 // }
1000 // }
1001 // else if(i==1){ // Sheet 1
1002 // for(int k=0; k<vecID[1][1].size(); k++){
1003 // if(vecID[1][1][k]==ID){
1004 // lv_new_digi->SetGlobalTime ( vecT[1][1][k] );
1005 // lv_new_digi->SetEnergyDeposit ( vecQ[1][1][k] );
1006 // }
1007 // }
1008 // }
1009 //}
1010 lv_new_digi->SetGlobalTime ( vecT[i][j][iStrip] );
1011 double EnergyDeposit_temp=vecQ[i][j][iStrip];
1012 if(EnergyDeposit_temp > Q_saturation) EnergyDeposit_temp = Q_saturation;
1013 lv_new_digi->SetEnergyDeposit ( EnergyDeposit_temp );
1014
1015 m_digi_collection->insert(lv_new_digi);
1016 lvi_F_hit = m_digi_collection->entries() - 1;
1017
1018 //Add by sunxh
1019 unsigned int ident = CgemID::getIntID(lvi_ID_layer,i,j,ID);
1020 ident_tmp.push_back(ident);
1021 //End Add
1022
1023 } // create a new digi
1024 else
1025 {
1026 //if(j==0){ // X strip
1027 // if(i==0){ // Sheet 0
1028 // for(int k=0; k<vecID[0][0].size(); k++){
1029 // if(vecID[0][0][k]==ID){
1030 // double EnergyDeposit_temp = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit()+vecQ[0][0][k];
1031 // double GlobalTime_temp = (*m_digi_collection)[lvi_F_hit]->GetGlobalTime();
1032 // if(GlobalTime_temp > Q_saturation) GlobalTime_temp = Q_saturation;
1033 // (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(EnergyDeposit_temp);
1034 // if(vecT[0][0][k]<GlobalTime_temp) (*m_digi_collection)[lvi_F_hit]->SetGlobalTime(vecT[0][0][k]);
1035 // }
1036 // }
1037 // }
1038 // else if(i==1){ // Sheet 1
1039 // for(int k=0; k<vecID[1][0].size(); k++){
1040 // if(vecID[1][0][k]==ID){
1041 // double EnergyDeposit_temp = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit()+vecQ[1][0][k];
1042 // double GlobalTime_temp = (*m_digi_collection)[lvi_F_hit]->GetGlobalTime();
1043 // if(GlobalTime_temp > Q_saturation) GlobalTime_temp = Q_saturation;
1044 // (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(EnergyDeposit_temp);
1045 // if(vecT[1][0][k]<GlobalTime_temp) (*m_digi_collection)[lvi_F_hit]->SetGlobalTime(vecT[1][0][k]);
1046 // }
1047 // }
1048 // }
1049 //}
1050 //else if(j==1){ // V strip
1051 // if(i==0){ // Sheet 0
1052 // for(int k=0; k<vecID[0][1].size(); k++){
1053 // if(vecID[0][1][k]==ID){
1054 // double EnergyDeposit_temp = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit()+vecQ[0][1][k];
1055 // double GlobalTime_temp = (*m_digi_collection)[lvi_F_hit]->GetGlobalTime();
1056 // if(GlobalTime_temp > Q_saturation) GlobalTime_temp = Q_saturation;
1057 // (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(EnergyDeposit_temp);
1058 // if(vecT[0][1][k]<GlobalTime_temp) (*m_digi_collection)[lvi_F_hit]->SetGlobalTime(vecT[0][1][k]);
1059 // }
1060 // }
1061 // }
1062 // else if(i==1){ // Sheet 1
1063 // for(int k=0; k<vecID[1][1].size(); k++){
1064 // if(vecID[1][1][k]==ID){
1065 // double EnergyDeposit_temp = (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit()+vecQ[1][1][k];
1066 // double GlobalTime_temp = (*m_digi_collection)[lvi_F_hit]->GetGlobalTime();
1067 // if(GlobalTime_temp > Q_saturation) GlobalTime_temp = Q_saturation;
1068 // (*m_digi_collection)[lvi_F_hit]->SetEnergyDeposit(EnergyDeposit_temp);
1069 // if(vecT[1][1][k]<GlobalTime_temp) (*m_digi_collection)[lvi_F_hit]->SetGlobalTime(vecT[1][1][k]);
1070 // }
1071 // }
1072 // }
1073 //}
1074
1075 // update energy
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);
1079
1080 // update time
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]);
1083 //-----check zhaojy---------------------
1084 //std::cout << "zhaojy check simulation --> layer:" << lvi_ID_layer << " sheet:" << i << " x/v:" << j << " ID:" << ID << " Q:" << (*m_digi_collection)[lvi_F_hit]->GetEnergyDeposit() << " T:" << (*m_digi_collection)[lvi_F_hit]->GetGlobalTime() << std::endl;
1085 //--------------------------------------
1086 //Add by sunxh
1087 unsigned int ident = CgemID::getIntID(lvi_ID_layer,i,j,ID);
1088 ident_tmp.push_back(Identifier(ident));
1089 //End Add
1090 }// else update an existing digi
1091 }// loop strip
1092 }// loop XV
1093 }// loop sheet
1094
1095 /* Add by sunxh */
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];
1101 }
1102 lv_hit->AddIdentifier(array_tmp,size_tmp);
1103 // TArrayI tmp = lv_hit->GetIdentifier();
1104 /*End Add*/
1105
1106 }// loop hits ii
1107
1108 StoreDigiCollection(m_digi_collection); /* */
1109
1110 //-----check zhaojy---------------------
1111 if(printFlag) {
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;
1119 }
1120 }
1121 //--------------------------------------
1122
1123 }//if(lv_HC)
1124}
1125
1126
1127
1128
1129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1130void BesCgemDigitizer::DoLorentzDiffusion(const G4int f_ID_layer, G4ThreeVector &f_XYZ_pre, G4ThreeVector &f_XYZ_post, const int f_F_print)
1131{
1132 IMessageSvc* msgSvc;
1133 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
1134 MsgStream log(msgSvc, "BesCgemDigitizer::DoLorentzDiffusion()");
1135 /* !!! Attention! JUST a SIMPLE-TEST-VERSION, it is not a real one!!! */
1136 CgemGeoLayer *f_cgem_layer = m_cgem_geomsvc->getCgemLayer(f_ID_layer);
1137 G4double lvd_R_layer = f_cgem_layer->getInnerROfAnode()*mm;
1138
1139 G4double lvd_A_phi_pre = f_XYZ_pre.phi();
1140 //G4double lvd_A_theta_pre = f_XYZ_pre.theta();
1141 G4double lvd_A_phi_post = f_XYZ_post.phi();
1142 //G4double lvd_A_theta_post = f_XYZ_post.theta();
1143
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() ;
1148
1149 /* Print information */
1150 if (f_F_print == 1)
1151 {
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()
1159 << endreq;
1160 log<< MSG::INFO << "" << endreq;
1161 }
1162
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 ));
1165 // f_XYZ_pre.setZ ( lvd_R_layer / tan ( lvd_A_theta_pre ));
1166
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 ));
1169 // f_XYZ_post.setZ ( lvd_R_layer / tan ( lvd_A_theta_post ));
1170
1171 /* Print information */
1172 if (f_F_print == 1)
1173 {
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()
1181 << endreq;
1182 log<< MSG::INFO << "" << endreq;
1183
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())
1189 << endreq;
1190 log<< MSG::INFO << "" << endreq;
1191 }
1192}
1193
1194//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1196{
1197}
1198
1199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1201{
1202}
1203
1204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1205/* Return ID from XYZ of readout layer */
1206void BesCgemDigitizer::GetIDFromXYZ(const G4int f_ID_layer, const G4ThreeVector f_XYZ,
1207 G4double &f_x, G4double &f_v, G4int &f_ID_sheet, G4int &f_ID_x, G4int &f_ID_v)
1208{
1209 IMessageSvc* msgSvc;
1210 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
1211 MsgStream log(msgSvc, "BesCgemDigitizer::GetIDFromXYZ()");
1212 /* CgemLayer information about readout strip */
1213 CgemGeoLayer *lv_cgem_layer = m_cgem_geomsvc->getCgemLayer(f_ID_layer);
1214 G4double lvd_R_layer = lv_cgem_layer->getInnerROfAnode()*mm;
1215 G4double lvd_L_layer = lv_cgem_layer->getLengthOfCgemLayer()*mm;
1216 G4double lvd_A_stero = lv_cgem_layer->getAngleOfStereo()*deg;
1217 G4double lvd_W_sheet = lv_cgem_layer->getWidthOfSheet()*mm;
1218 G4double lvd_W_pitch_x = lv_cgem_layer->getWidthOfPitchX()*mm;
1219 G4double lvd_W_pitch_v = lv_cgem_layer->getWidthOfPitchV()*mm;
1220 G4int lvi_N_sheet = lv_cgem_layer->getNumberOfSheet();
1221
1222 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1223
1224 /* lvd_A_phi should be in (0,twopi) */
1225 G4double lvd_A_phi = f_XYZ.phi();
1226 if (lvd_A_phi < 0)
1227 {
1228 lvd_A_phi += twopi;
1229 }
1230 if (lvd_A_phi >= twopi)
1231 {
1232 lvd_A_phi -= twopi;
1233 }
1234
1235 /* Transfer Cylinder to Rectangle */
1236 G4double lvd_x = f_XYZ.rho() * lvd_A_phi;
1237 G4double lvd_v = f_XYZ.z() + lvd_L_layer/2.;
1238
1239
1240 /* ID of Sheet */
1241 if (lvd_x < 0 or lvd_x > lvi_N_sheet*lvd_W_sheet)
1242 {
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
1247 << endreq;
1248 return ;
1249 }
1250 else
1251 {
1252 for (G4int i=0; i<lvi_N_sheet; i++)
1253 {
1254 if (lvd_x >= i*lvd_W_sheet and lvd_x < (i+1)*lvd_W_sheet)
1255 {
1256 f_ID_sheet = i;
1257 }
1258 }
1259 }
1260
1261 /* ID of Strip X and V */
1262 f_x = lvd_x - f_ID_sheet * lvd_W_sheet;
1263 f_v = lvd_v;
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);
1266}
1267
1268//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1269/* Return XYZ of readout layer from assigned ID */
1270void BesCgemDigitizer::GetIDFromSegmentInSameSheet(const G4int f_ID_x_pre , const G4int f_ID_v_pre ,
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 )
1274{
1275 /* Number of hit strips */
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;
1278
1279 /* Starting strip of hit strip */
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;
1282}
1283
1284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1285/* Return X and V strip ID of the middle point by the segment which crossed the sheet */
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)
1291{
1292 CgemGeoLayer *lv_cgem_layer = m_cgem_geomsvc->getCgemLayer(f_ID_layer);
1293 G4double lvd_R_layer = lv_cgem_layer->getInnerROfAnode()*mm;
1294 G4double lvd_W_sheet = lv_cgem_layer->getWidthOfSheet()*mm;
1295 G4int lvi_N_sheet = lv_cgem_layer->getNumberOfSheet();
1296 G4double lvd_W_pitch_x = lv_cgem_layer->getWidthOfPitchX()*mm;
1297 G4double lvd_W_pitch_v = lv_cgem_layer->getWidthOfPitchV()*mm;
1298 G4double lvd_A_stero = lv_cgem_layer->getAngleOfStereo()*deg;
1299 G4double lvd_v_mid = 0;
1300
1301 lvd_W_sheet = lvd_R_layer * twopi / lvi_N_sheet;
1302
1303 G4int lvi_N_strip_x = ceil(lvd_W_sheet / lvd_W_pitch_x);
1304
1305 /* Sheet ID of mid point */
1306 f_ID_sheet_mid_pre = f_ID_sheet_pre;
1307 f_ID_sheet_mid_post = f_ID_sheet_post;
1308
1309 /* Strip ID of mid point */
1310 if (f_x_pre > f_x_post)
1311 {
1312 /* X Strip ID of mid point */
1313 f_ID_x_mid_pre = lvi_N_strip_x - 1;
1314 f_ID_x_mid_post = 0;
1315
1316 /* V Strip ID of mid point */
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;
1318
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);
1321 }
1322 else
1323 {
1324 /* X Strip ID of mid point */
1325 f_ID_x_mid_pre = 0;
1326 f_ID_x_mid_post = lvi_N_strip_x - 1;
1327
1328 /* V Strip ID of mid point */
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;
1330
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);
1333 }
1334}
1335
1336//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1337/* Return XYZ of readout layer from assigned ID */
1339{
1340}
1341
1342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Double_t x[10]
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
G4TDigiCollection< BesCgemDigi > BesCgemDigisCollection
G4THitsCollection< BesCgemHit > BesCgemHitsCollection
double sin(const BesAngle a)
double cos(const BesAngle a)
void SetGlobalTime(G4double f_global_time)
void SetEnergyDeposit(G4double f_E_deposit)
BesCgemDigitizer(G4String modName)
virtual void Digitize()
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 AddIdentifier(G4int f_ID_Identifier[2000], G4int N_dim)
StatusCode setTrack(int layer, int particle, int charge, double p, double trkPosIn[], double trkPosOut[])
void getFiredStripID(G4ThreeVector pos1, G4ThreeVector pos2, vector< int > &vecXID, vector< int > &vecVID) 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)