CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
InductionGar.cxx
Go to the documentation of this file.
1#include "CgemDigitizerSvc/InductionGar.h"
2
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/Bootstrap.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include "GaudiKernel/DataSvc.h"
8
9#include "CLHEP/Units/PhysicalConstants.h"
10
11#include "G4DigiManager.hh"
12#include "Randomize.hh"
13#include "G4ios.hh"
14
15#include "CLHEP/Units/PhysicalConstants.h"
16
17#include <iomanip>
18#include <iostream>
19#include <fstream>
20#include <cmath>
21
22#include "TRandom.h"
23
24typedef std::vector<int> Vint;
25using namespace std;
26
27//-----------------------Convolution of response function---------------------------------
28
29double rectangle2(double t, double *x, double *y, int Npx)
30{
31
32 double I = y[0];
33 int id = -1;
34 for(int i = 0; i < Npx-1; i++){
35 if(t>=x[i]&&t<x[i+1]){id = i; break;}
36 }
37 if(id != -1){I=y[id]+(y[id+1]-y[id])*(t-x[id])/(x[id+1]-x[id]);}
38
39 return I;
40}
41
42double T_branch2(double t)
43{
44 double result = 0.;
45 t=t-10.;
46 if(t>0) {
47 result = 2000.*(0.00181928*exp(-t/3.)-0.0147059*exp(-t/20.)+0.0128866*exp(-t/100.));
48 }
49 return result;
50}
51
52TH1F Convolution_Tbranch(double *Input_Curr_plot_001) // return TIGER output
53{
54 double xmin(0), xmax(1000);
55 const int Npx(2000);
56
57 TH1F h_signal("h_signal","",Npx,xmin,xmax);
58 double Input_Time_plot_001[Npx];
59 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signal.GetBinCenter(i+1);
60
61 double xmin_f1(0.), xmax_f1(1000);
62 int nStep_f1(200);
63 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
64 double x_f1(0.0);
65
66 double x[Npx];
67 double y_001[Npx];
68
69 for(int i=0; i<Npx; i++)
70 {
71 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
72 double fun_001=0.0;
73
74 for(int j=0; j<nStep_f1; j++)
75 {
76 x_f1=xmin_f1+step_f1*(j+0.5);
77 fun_001+=rectangle2(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*T_branch2(x[i]-x_f1)*step_f1;
78
79 }
80 y_001[i]=fun_001;
81 }
82
83 for(int init=0; init<Npx; init++){
84 h_signal.SetBinContent(init+1,y_001[init]);
85 }
86
87 return h_signal;
88}
89
90
91InductionGar::InductionGar(){
92}
93
94InductionGar::~InductionGar(){
95}
96
97void InductionGar::init(ICgemGeomSvc* geomSvc, double magConfig){
98 m_geomSvc = geomSvc;
99 m_magConfig = magConfig;
100
101 string filePath_ratio = getenv("CGEMDIGITIZERSVCROOT");
102 string fileName;
103 if(0==m_magConfig) fileName = filePath_ratio + "/dat/par_InductionGar_0T.txt";
104 else fileName = filePath_ratio + "/dat/par_InductionGar_1T.txt";
105 ifstream fin(fileName.c_str(), ios::in);
106 cout << "fileName: " << fileName << endl;
107
108 string strline;
109 string strpar;
110 if( ! fin.is_open() ){
111 cout << "ERROR: can not open file " << fileName << endl;
112 }
113
114 for(int m=0; m<3; m++){
115 for(int l=0; l<5; l++){
116 for(int k=0; k<5; k++){
117 for(int i=0; i<4; i++){
118 fin>>RatioX[m][l][k][i];
119 }
120 for(int j=0; j<3; j++){
121 fin>>RatioV[m][l][k][j];
122 }
123 }
124 }
125 }
126 fin.close();
127
128 string filePath = getenv("CGEMDIGITIZERSVCROOT");
129 for(int i=0; i<3; i++){ // layer
130 for(int j=0; j<5; j++){ // Grid-x
131 for(int k=0; k<5; k++){ // Grid-v
132 char temp1[1]; sprintf(temp1, "%i", (i+1));
133 char temp2[2]; sprintf(temp2, "%i", ((j+1)*10+k+1));
134 if(0==m_magConfig) fileName = filePath + "/dat/InducedCurrent_Root_0T/layer" + temp1 + "_" + temp2 + ".root";
135 else fileName = filePath + "/dat/InducedCurrent_Root_1T/layer" + temp1 + "_" + temp2 + ".root";
136 TFile* file_00 = TFile::Open(fileName.c_str(),"read");
137 TH1* readThis_x3 = 0;
138 TH1* readThis_x4 = 0;
139 TH1* readThis_x5 = 0;
140 TH1* readThis_x6 = 0;
141 TH1* readThis_v5 = 0;
142 TH1* readThis_v6 = 0;
143 TH1* readThis_v7 = 0;
144 file_00->GetObject("readThis_x3", readThis_x3);
145 file_00->GetObject("readThis_x4", readThis_x4); file_00->GetObject("readThis_v5", readThis_v5);
146 file_00->GetObject("readThis_x5", readThis_x5); file_00->GetObject("readThis_v6", readThis_v6);
147 file_00->GetObject("readThis_x6", readThis_x6); file_00->GetObject("readThis_v7", readThis_v7);
148 for(int init = 0; init<100; init++){
149 SignalX[i][j][k][0][init] = readThis_x3->GetBinContent(init+1);
150 SignalX[i][j][k][1][init] = readThis_x4->GetBinContent(init+1);
151 SignalX[i][j][k][2][init] = readThis_x5->GetBinContent(init+1);
152 SignalX[i][j][k][3][init] = readThis_x6->GetBinContent(init+1);
153 SignalV[i][j][k][0][init] = readThis_v5->GetBinContent(init+1);
154 SignalV[i][j][k][1][init] = readThis_v6->GetBinContent(init+1);
155 SignalV[i][j][k][2][init] = readThis_v7->GetBinContent(init+1);
156 }
157 }
158 }
159 }
160}
161
162void InductionGar::setMultiElectrons(int layer, int nElectrons, std::vector<double> x, std::vector<double> y, std::vector<double> z, std::vector<double> t){
163
164 m_xstripSheet.clear();
165 m_xstripID.clear();
166 m_vstripSheet.clear();
167 m_vstripID.clear();
168 m_xstripQ.clear();
169 m_vstripQ.clear();
170 m_xstripT.clear();
171 m_vstripT.clear();
172
173 int m000_nXstrips;
174 int m000_nVstrips;
175 std::vector<int> m000_xstripSheet;
176 std::vector<int> m000_xstripID;
177 std::vector<int> m000_vstripSheet;
178 std::vector<int> m000_vstripID;
179 std::vector<double> m000_xstripQ;
180 std::vector<double> m000_vstripQ;
181 std::vector<double> m000_xstripT;
182 std::vector<double> m000_vstripT;
183
184 m000_xstripSheet.clear();
185 m000_xstripID.clear();
186 m000_vstripSheet.clear();
187 m000_vstripID.clear();
188 m000_xstripQ.clear();
189 m000_vstripQ.clear();
190 m000_xstripT.clear();
191 m000_vstripT.clear();
192
193 CgemGeoLayer* CgemLayer;
194 CgemGeoReadoutPlane* CgemReadoutPlane;
195 CgemLayer = m_geomSvc->getCgemLayer(layer);
196 int NSheet = CgemLayer->getNumberOfSheet();
197
198 Vint Save_Sheetx, Save_Sheetv;
199 Vint Save_FinalstripIDx, Save_FinalstripIDv;
200 Vint Save_Gridx, Save_Gridv;
201 Vint Save_IDx, Save_IDv;
202 Vint Save_Tx, Save_Tv;
203 Save_Tx.clear();
204 Save_Tv.clear();
205 Save_IDx.clear();
206 Save_IDv.clear();
207 Save_Gridx.clear();
208 Save_Gridv.clear();
209 Save_Sheetx.clear();
210 Save_Sheetv.clear();
211 Save_FinalstripIDx.clear();
212 Save_FinalstripIDv.clear();
213 for(int i=0; i<2; i++){
214 for(int k=0; k<2; k++){
215 m_mapQ[i][k].clear();
216 //m_mapT[i][k].clear();
217 }
218 }
219
220 // std::cout << "nElectrons = " << nElectrons << std::endl;
221 for(int i=0; i<nElectrons; i++){
222 double phi_electron = atan2(y[i], x[i]);
223 double z_electron = z[i];
224 G4ThreeVector pos(x[i], y[i], z[i]);
225 for(int j=0; j<NSheet; j++)
226 {
227 CgemReadoutPlane = m_geomSvc->getReadoutPlane(layer, j);
228 if(CgemReadoutPlane->OnThePlane(phi_electron,z_electron))
229 {
230 int xID;
231 double distX = CgemReadoutPlane->getDist2ClosestXStripCenter(phi_electron, xID);
232 int vID;
233 double distV = CgemReadoutPlane->getDist2ClosestVStripCenter(pos, vID);
234
235 //std::cout << "-------------" << std::endl;
236 //std::cout << "xID, distX = " << xID << "," << distX << std::endl;
237 //std::cout << "vID, distV = " << vID << "," << distV << std::endl;
238
239 int grid_x = 1;
240 int grid_v = 1;
241
242 //
243 //--- Calculation of grid index-------------
244 //
245 double L_XPitch = CgemReadoutPlane->getXPitch();
246 double L_VPitch = CgemReadoutPlane->getVPitch();
247 grid_v = floor(distX/L_XPitch*5.+2.5);
248 grid_x = floor(distV/L_VPitch*5.+2.5);
249 //std::cout << i << " " << j << " " << grid_x << " " << grid_v << std::endl;
250 //std::cout << "L_XPitch,L_VPitch = " << L_XPitch << "," << L_VPitch << std::endl;
251 //std::cout << "grid_x, grid_v = " << grid_x << "," << grid_v << std::endl;
252
253 if(xID>=0 && vID>=0) {
254 map<int, double>::iterator itx = m_mapQ[j][0].find(xID);
255 map<int, double>::iterator itv = m_mapQ[j][1].find(vID);
256 if(RatioX[layer][grid_x][grid_v][2]!=0){
257 if(itx==m_mapQ[j][0].end())
258 {
259 m_mapQ[j][0][xID]=RatioX[layer][grid_x][grid_v][2];
260 Save_Gridx.push_back(grid_x*100+grid_v*10+2);
261 Save_IDx.push_back(j*10000+xID);
262 Save_Tx.push_back(t[i]);
263 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][2] << std::endl;
264 }
265 else {
266 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][2];
267 m_mapQ[j][0][xID]=Q;
268 Save_Gridx.push_back(grid_x*100+grid_v*10+2);
269 Save_IDx.push_back(j*10000+xID);
270 Save_Tx.push_back(t[i]);
271 //std::cout << i << " " << j << "x-strip : Add " << Q << std::endl;
272 }
273 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID);
274 if(result_Sheetx==Save_Sheetx.end()){
275 Save_Sheetx.push_back(j*10000+xID);
276 }
277 }
278 if(RatioV[layer][grid_x][grid_v][1]!=0){
279 if(itv==m_mapQ[j][1].end())
280 {
281 m_mapQ[j][1][vID]=RatioV[layer][grid_x][grid_v][1];
282 Save_Gridv.push_back(grid_x*100+grid_v*10+1);
283 Save_IDv.push_back(j*10000+vID);
284 Save_Tv.push_back(t[i]);
285 }
286 else {
287 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][1];
288 m_mapQ[j][1][vID]=Q;
289 Save_Gridv.push_back(grid_x*100+grid_v*10+1);
290 Save_IDv.push_back(j*10000+vID);
291 Save_Tv.push_back(t[i]);
292 }
293 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID);
294 if(result_Sheetv==Save_Sheetv.end()){
295 Save_Sheetv.push_back(j*10000+vID);
296 }
297 }
298 }
299
300 if((xID>=0 && vID>=0) && (xID+1)<CgemReadoutPlane->getNXstrips()) {
301 map<int, double>::iterator itx = m_mapQ[j][0].find(xID+1);
302 if(RatioX[layer][grid_x][grid_v][3]!=0){
303 if(itx==m_mapQ[j][0].end())
304 {
305 m_mapQ[j][0][xID+1]=RatioX[layer][grid_x][grid_v][3];
306 Save_Gridx.push_back(grid_x*100+grid_v*10+3);
307 Save_IDx.push_back(j*10000+xID+1);
308 Save_Tx.push_back(t[i]);
309 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][3] << std::endl;
310 }
311 else {
312 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][3];
313 m_mapQ[j][0][xID+1]=Q;
314 Save_Gridx.push_back(grid_x*100+grid_v*10+3);
315 Save_IDx.push_back(j*10000+xID+1);
316 Save_Tx.push_back(t[i]);
317 }
318
319 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID+1);
320 if(result_Sheetx==Save_Sheetx.end()){
321 Save_Sheetx.push_back(j*10000+xID+1);
322 }
323 }
324 }
325
326 if((xID>=0 && vID>=0) && (vID+1)<CgemReadoutPlane->getNVstrips()) {
327 map<int, double>::iterator itv = m_mapQ[j][1].find(vID+1);
328 if(RatioV[layer][grid_x][grid_v][2]!=0){
329 if(itv==m_mapQ[j][1].end())
330 {
331 m_mapQ[j][1][vID+1]=RatioV[layer][grid_x][grid_v][2];
332 Save_Gridv.push_back(grid_x*100+grid_v*10+2);
333 Save_IDv.push_back(j*10000+vID+1);
334 Save_Tv.push_back(t[i]);
335 }
336 else {
337 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][2];
338 m_mapQ[j][1][vID+1]=Q;
339 Save_Gridv.push_back(grid_x*100+grid_v*10+2);
340 Save_IDv.push_back(j*10000+vID+1);
341 Save_Tv.push_back(t[i]);
342 }
343
344 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID+1);
345 if(result_Sheetv==Save_Sheetv.end()){
346 Save_Sheetv.push_back(j*10000+vID+1);
347 }
348 }
349 }
350
351 if((xID>=0 && vID>=0) && (xID-2)>=0) {
352 map<int, double>::iterator itx = m_mapQ[j][0].find(xID-2);
353 if(RatioX[layer][grid_x][grid_v][0]!=0){
354 if(itx==m_mapQ[j][0].end())
355 {
356 m_mapQ[j][0][xID-2]=RatioX[layer][grid_x][grid_v][0];
357 Save_Gridx.push_back(grid_x*100+grid_v*10+0);
358 Save_IDx.push_back(j*10000+xID-2);
359 Save_Tx.push_back(t[i]);
360 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][0] << std::endl;
361 }
362 else {
363 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][0];
364 m_mapQ[j][0][xID-2]=Q;
365 Save_Gridx.push_back(grid_x*100+grid_v*10+0);
366 Save_IDx.push_back(j*10000+xID-2);
367 Save_Tx.push_back(t[i]);
368 }
369
370 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID-2);
371 if(result_Sheetx==Save_Sheetx.end()){
372 Save_Sheetx.push_back(j*10000+xID-2);
373 }
374 }
375 }
376
377 if((xID>=0 && vID>=0) && (xID-1)>=0) {
378 map<int, double>::iterator itx = m_mapQ[j][0].find(xID-1);
379 if(RatioX[layer][grid_x][grid_v][1]!=0){
380 if(itx==m_mapQ[j][0].end())
381 {
382 m_mapQ[j][0][xID-1]=RatioX[layer][grid_x][grid_v][1];
383 Save_Gridx.push_back(grid_x*100+grid_v*10+1);
384 Save_IDx.push_back(j*10000+xID-1);
385 Save_Tx.push_back(t[i]);
386 }
387 else {
388 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][1];
389 m_mapQ[j][0][xID-1]=Q;
390 Save_Gridx.push_back(grid_x*100+grid_v*10+1);
391 Save_IDx.push_back(j*10000+xID-1);
392 Save_Tx.push_back(t[i]);
393 }
394
395 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID-1);
396 if(result_Sheetx==Save_Sheetx.end()){
397 Save_Sheetx.push_back(j*10000+xID-1);
398 }
399 }
400 }
401
402 if((xID>=0 && vID>=0) && (vID-1)>=0) {
403 map<int, double>::iterator itv = m_mapQ[j][1].find(vID-1);
404 if(RatioV[layer][grid_x][grid_v][0]!=0){
405 if(itv==m_mapQ[j][1].end())
406 {
407 m_mapQ[j][1][vID-1]=RatioV[layer][grid_x][grid_v][0];
408 Save_Gridv.push_back(grid_x*100+grid_v*10+0);
409 Save_IDv.push_back(j*10000+vID-1);
410 Save_Tv.push_back(t[i]);
411 }
412 else {
413 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][0];
414 m_mapQ[j][1][vID-1]=Q;
415 Save_Gridv.push_back(grid_x*100+grid_v*10+0);
416 Save_IDv.push_back(j*10000+vID-1);
417 Save_Tv.push_back(t[i]);
418 }
419
420 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID-1);
421 if(result_Sheetv==Save_Sheetv.end()){
422 Save_Sheetv.push_back(j*10000+vID-1);
423 }
424 }
425 }
426 }
427 }
428 }
429 // std::cout << "loop electron finished " << std::endl;
430
431 //-----Q Threshold = 45ADC = 1.5fC = 1.5*1e-15 C --------
432 //-----Q Satruation = 45fC
433 //-----Q Resolution = 0.256fC
434 //-----T jitter = 2ns
435
436 double Q_threshold = 1.5e-15/1.602176565e-19;
437 double Q_saturation= 45.e-15/1.602176565e-19;
438 double Q_resolution= 0.256e-15/1.602176565e-19;
439 double T_threshold = 12; // 12 mV
440 //double T_threshold = 24; // 24 mV
441 //double T_threshold = 18; // 18 mV
442
443 // cout << "---------------------" << endl;
444 // consider charge resolution and saturation
445 for(int i=0; i<Save_Sheetx.size(); i++){
446 int sheetx_id = Save_Sheetx[i]/10000;
447 int stripx_id = Save_Sheetx[i]%10000;
448 //cout << "before " << m_mapQ[sheetx_id][0][stripx_id] << endl;
449 m_mapQ[sheetx_id][0][stripx_id]=gRandom->Gaus(m_mapQ[sheetx_id][0][stripx_id],Q_resolution);
450 if(m_mapQ[sheetx_id][0][stripx_id]>Q_saturation) m_mapQ[sheetx_id][0][stripx_id] = Q_saturation;
451 //cout << "after " << m_mapQ[sheetx_id][0][stripx_id] << endl;
452 }
453
454 for(int i=0; i<Save_Sheetv.size(); i++){
455 int sheetv_id = Save_Sheetv[i]/10000;
456 int stripv_id = Save_Sheetv[i]%10000;
457 //cout << "before " << m_mapQ[sheetv_id][1][stripv_id] << endl;
458 m_mapQ[sheetv_id][1][stripv_id]=gRandom->Gaus(m_mapQ[sheetv_id][1][stripv_id],Q_resolution);
459 if(m_mapQ[sheetv_id][1][stripv_id]>Q_saturation) m_mapQ[sheetv_id][1][stripv_id] = Q_saturation;
460 //cout << "after " << m_mapQ[sheetv_id][1][stripv_id] << endl;
461 }
462
463
464 int Nx_below_threshold = 0;
465 for(int i=0; i<Save_Sheetx.size(); i++){
466 int sheetx_id = Save_Sheetx[i]/10000;
467 int stripx_id = Save_Sheetx[i]%10000;
468 //std::cout << m_mapQ[sheetx_id][0][stripx_id] << std::endl;
469 if (m_mapQ[sheetx_id][0][stripx_id]<Q_threshold) Nx_below_threshold++;
470 }
471
472 int Nv_below_threshold = 0;
473 for(int i=0; i<Save_Sheetv.size(); i++){
474 int sheetv_id = Save_Sheetv[i]/10000;
475 int stripv_id = Save_Sheetv[i]%10000;
476 //std::cout << m_mapQ[sheetv_id][1][stripv_id] << std::endl;
477 if (m_mapQ[sheetv_id][1][stripv_id]<Q_threshold) Nv_below_threshold++;
478 }
479
480 //-------------------------------------------------------------
481
482 m000_nXstrips = m_mapQ[0][0].size() + m_mapQ[1][0].size() - Nx_below_threshold;
483 m000_nVstrips = m_mapQ[0][1].size() + m_mapQ[1][1].size() - Nv_below_threshold;
484
485 //cout << m000_nXstrips << " " << Nx_below_threshold << " " << m000_nVstrips << " " << Nv_below_threshold << endl;
486
487
488 for(int i=0; i<Save_Sheetx.size(); i++){
489 int sheetx_id = Save_Sheetx[i]/10000;
490 int stripx_id = Save_Sheetx[i]%10000;
491 if(m_mapQ[sheetx_id][0][stripx_id]>=Q_threshold){
492 Save_FinalstripIDx.push_back(Save_Sheetx[i]);
493 //std::cout << "zhaojy check InductionGar sheet -> " << m000_xstripSheet.size() << " " << sheetx_id << " " << stripx_id << std::endl;
494 m000_xstripSheet.push_back(sheetx_id);
495 m000_xstripID.push_back(stripx_id);
496 m000_xstripQ.push_back(m_mapQ[sheetx_id][0][stripx_id]);
497 }
498 }
499
500 for(int i=0; i<Save_Sheetv.size(); i++){
501 int sheetv_id = Save_Sheetv[i]/10000;
502 int stripv_id = Save_Sheetv[i]%10000;
503 if(m_mapQ[sheetv_id][1][stripv_id]>=Q_threshold){
504 Save_FinalstripIDv.push_back(Save_Sheetv[i]);
505 m000_vstripSheet.push_back(sheetv_id);
506 m000_vstripID.push_back(stripv_id);
507 m000_vstripQ.push_back(m_mapQ[sheetv_id][1][stripv_id]);
508 }
509 }
510
511
512 //string FilePath = getenv("TESTCGEMDIGISVCALGROOT");
513 for(int h=0; h<Save_FinalstripIDx.size(); h++){
514 const int _low = 0;
515 const int _up = 1000;
516 const int _nbins = 2000;
517 double binsize = (double)(_up-_low)/_nbins;
518 int sheetx_id = Save_FinalstripIDx[h]/10000;
519 int stripx_id = Save_FinalstripIDx[h]%10000;
520 double histX_bin_Content[_nbins];
521 for(int i=0; i<_nbins; i++){
522 histX_bin_Content[i] = 0;
523 }
524 for(int i=0; i<Save_Gridx.size(); i++){
525 int z_grid_x = Save_Gridx[i]/100;
526 int z_grid_v = Save_Gridx[i]%100/10;
527 int z_index = Save_Gridx[i]%10;
528 int z_IDx = Save_IDx[i];
529 int start_bin = Save_Tx[i]/binsize;
530 //std::cout << "start_bin x" << start_bin << std::endl;
531 if(z_IDx == Save_FinalstripIDx[h]){
532 for(int addi = start_bin; addi < start_bin+160; addi+=2){
533 histX_bin_Content[addi]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
534 histX_bin_Content[addi+1]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
535 }
536 }
537 }
538
539 TH1F h_signal = Convolution_Tbranch(histX_bin_Content);
540
541 int flg_xT = 0;
542 for(int init=1; init<_nbins; init++){
543 if(flg_xT==0 && h_signal.GetBinContent(init+1)<=-T_threshold){
544 m000_xstripT.push_back(_low+binsize*(init-1)+fabs(-T_threshold-h_signal.GetBinContent(init))*binsize/fabs(h_signal.GetBinContent(init+1)-h_signal.GetBinContent(init))+gRandom->Gaus(0,2)); // jitter 2ns
545 flg_xT=1;
546 }
547 if(flg_xT==0 && h_signal.GetBinContent(init+1)>=T_threshold){
548 m000_xstripT.push_back(_low+binsize*(init-1)+fabs(T_threshold-h_signal.GetBinContent(init))*binsize/fabs(h_signal.GetBinContent(init+1)-h_signal.GetBinContent(init))+gRandom->Gaus(0,2)); // jitter 2ns
549 flg_xT=1;
550 }
551 }
552 if(flg_xT==0) m000_xstripT.push_back(-99);
553
554 //char temp1[1]; sprintf(temp1, "%i", sheetx_id);
555 //char temp2[3]; sprintf(temp2, "%i", stripx_id);
556 //string FILEname = FilePath + "/share/output/X-" + temp1 + '-' + temp2 + ".root";
557 //TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); h_signal.Write(); f_x1->Close(); delete f_x1;
558 }
559
560 for(int h=0; h<Save_FinalstripIDv.size(); h++){
561 const int _low = 0;
562 const int _up = 1000;
563 const int _nbins = 2000;
564 double binsize = (double)(_up-_low)/_nbins;
565 int sheetv_id = Save_FinalstripIDv[h]/10000;
566 int stripv_id = Save_FinalstripIDv[h]%10000;
567 double histV_bin_Content[_nbins];
568 for(int i=0; i<_nbins; i++){
569 histV_bin_Content[i] = 0;
570 }
571 for(int i=0; i<Save_Gridv.size(); i++){
572 int z_grid_x = Save_Gridv[i]/100;
573 int z_grid_v = Save_Gridv[i]%100/10;
574 int z_index = Save_Gridv[i]%10;
575 int z_IDv = Save_IDv[i];
576 int start_bin = Save_Tv[i]/binsize;
577 //std::cout << "start_bin v" << start_bin << std::endl;
578 if(z_IDv == Save_FinalstripIDv[h]){
579 for(int addi = start_bin; addi < start_bin+160; addi+=2){
580 histV_bin_Content[addi]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
581 histV_bin_Content[addi+1]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
582 }
583 }
584 }
585
586 TH1F h_signal = Convolution_Tbranch(histV_bin_Content);
587
588 int flg_vT = 0;
589 for(int init=1; init<_nbins; init++){
590 if(flg_vT==0 && h_signal.GetBinContent(init+1)<=-T_threshold){
591 m000_vstripT.push_back(_low+binsize*(init-1)+fabs(-T_threshold-h_signal.GetBinContent(init))*binsize/fabs(h_signal.GetBinContent(init+1)-h_signal.GetBinContent(init))+gRandom->Gaus(0,2)); // jitter 2ns
592 flg_vT=1;
593 }
594 if(flg_vT==0 && h_signal.GetBinContent(init+1)>=T_threshold){
595 m000_vstripT.push_back(_low+binsize*(init-1)+fabs(T_threshold-h_signal.GetBinContent(init))*binsize/fabs(h_signal.GetBinContent(init+1)-h_signal.GetBinContent(init))+gRandom->Gaus(0,2)); // jitter 2ns
596 flg_vT=1;
597 }
598 }
599 if(flg_vT==0) m000_vstripT.push_back(-99);
600
601 //char temp1[1]; sprintf(temp1, "%i", sheetv_id);
602 //char temp2[3]; sprintf(temp2, "%i", stripv_id);
603 //string FILEname = FilePath + "/share/output/V-" + temp1 + '-' + temp2 + ".root";
604 //TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); h_signal.Write(); f_x1->Close(); delete f_x1;
605 }
606
607 //-------------------------------------------------------------------------------
608
609 double q_to_fC_factor = 1.602176565e-4;
610
611 m_nXstrips = 0;
612 for(int i=0; i<m000_nXstrips; i++){
613 if(m000_xstripT[i]>=0){
614 m_xstripSheet.push_back(m000_xstripSheet[i]);
615 m_xstripID.push_back(m000_xstripID[i]);
616 m_xstripQ.push_back(m000_xstripQ[i]*q_to_fC_factor);
617 m_xstripT.push_back(m000_xstripT[i]);
618 m_nXstrips++;
619 }
620 }
621 m_nVstrips = 0;
622 for(int i=0; i<m000_nVstrips; i++){
623 if(m000_vstripT[i]>=0){
624 m_vstripSheet.push_back(m000_vstripSheet[i]);
625 m_vstripID.push_back(m000_vstripID[i]);
626 m_vstripQ.push_back(m000_vstripQ[i]*q_to_fC_factor);
627 m_vstripT.push_back(m000_vstripT[i]);
628 m_nVstrips++;
629 }
630 }
631}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
rectangle2
double rectangle2(double t, double *x, double *y, int Npx)
Definition: InductionGar.cxx:29
Convolution_Tbranch
TH1F Convolution_Tbranch(double *Input_Curr_plot_001)
Definition: InductionGar.cxx:52
T_branch2
double T_branch2(double t)
Definition: InductionGar.cxx:42
Vint
std::vector< int > Vint
Definition: InductionGar.cxx:24
I
const DifComplex I
CgemGeoReadoutPlane::getDist2ClosestXStripCenter
double getDist2ClosestXStripCenter(double phi, int &id)
Definition: CgemGeoReadoutPlane.cxx:143
CgemGeoReadoutPlane::getDist2ClosestVStripCenter
double getDist2ClosestVStripCenter(G4ThreeVector pos, int &id)
Definition: CgemGeoReadoutPlane.cxx:217
ICgemGeomSvc::getReadoutPlane
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
ICgemGeomSvc::getCgemLayer
virtual CgemGeoLayer * getCgemLayer(int i) const =0
InductionGar::init
void init(ICgemGeomSvc *geomSvc, double magConfig)
Definition: InductionGar.cxx:97
InductionGar::setMultiElectrons
void setMultiElectrons(int layer, int nElectrons, std::vector< double > x, std::vector< double > y, std::vector< double > z, std::vector< double > t)
Definition: InductionGar.cxx:162
t
int t()
Definition: t.c:1