CGEM BOSS 6.6.5.h
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#include "GaudiKernel/PropertyMgr.h"
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#include "TTree.h"
22
23#include "TRandom.h"
24#include "TSystem.h"
25//#define INDUCTIONGAR_DEBUG_RECTANGLE2_FAST_TEST
26//uncomment the previous line to check correctness of results of rectangle2_fast.
27//it will be slow.
28
29//#define INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
30//uncomment the previous line to check correctness of results of Convolution_Tbranch_fft and Convolution_Ebranch_fft
31//it will be very slow.
32
33//#define INDUCTIONGAR_DEBUG_CONV1PERGRID_FFT_TEST
34//uncomment the previous line to check equality of 2 conv1PerGrid.
35
36const static int electrons_select_method_threhold=2500;
37
38typedef std::vector<int> Vint;
39using namespace std;
40static void checkMemorySize();
41
42//-----------------------Convolution of response function---------------------------------
43
44double rectangle2(double t, double *x, double *y, int Npx)
45{
46
47 double I = y[0];
48 int id = -1;
49 for(int i = 0; i < Npx-1; i++){
50 if(t>=x[i]&&t<x[i+1]){id = i; break;}
51 }
52 if(id != -1){I=y[id]+(y[id+1]-y[id])*(t-x[id])/(x[id+1]-x[id]);}
53
54 return I;
55}
56
57/**
58 * @brief calc estimate using linear interpolation of curve {x,y} , at point t; mimic rectangle2() but faster version
59 * assuming x: from 0 to 1000 and 2000 points.
60 * @param t
61 * @param x
62 * @param y
63 * @param Npx
64 * @return double
65 */
66double rectangle2_fast(double t, double *x,const double *y, int Npx) {
67
68 const int nbinsx=2000;
69 const double xmin=0;
70 const double xmax=1000;
71 const double binWidth=(xmax-xmin)/nbinsx;
72 double I = y[0];
73 int id = -1;
74
75 if (t > xmin+binWidth/2 && t< xmax-binWidth/2){
76 id=(int)floor((t-binWidth/2-xmin)/binWidth);
77 }
78#ifdef INDUCTIONGAR_DEBUG_RECTANGLE2_FAST_TEST // checking results of new and old version. for speed it should not run
79 int id2 = -1;
80 for (int i = 0; i < Npx - 1; i++) {
81 if (t >= x[i] && t < x[i + 1]) {
82 id2 = i;
83 break;
84 } // find the closist id, to make t=x[id]; find root.
85 }
86 if (id2 != id) {
87 cout<<
88 "CgemDigitizerSvc::InductionGar::rectangle2_fast: 2 method calced id doesnot correspond, "<<
89 "id new : "<< id2 <<" id old : "<<id<<
90 " t : "<<t<<" x at id old "<< x[id2]<< endl;
91 }
92#endif
93 if (id != -1) {
94 I = y[id] + (y[id + 1] - y[id]) * (t - x[id]) / (x[id + 1] - x[id]);
95 } // the tangent of y(i),x(i) at i=id, I is on tangent and at x=t;
96
97 return I;
98}
99
100double rectangle2(double t, std::vector<double>& x,const double *y, int Npx)
101{
102
103 double I = y[0];
104 int id = -1;
105 for(int i = 0; i < Npx-1; i++){
106 if(t>=x[i]&&t<x[i+1]){id = i; break;}
107 }
108 if(id != -1){I=y[id]+(y[id+1]-y[id])*(t-x[id])/(x[id+1]-x[id]);}
109
110 return I;
111}
112
113double T_branch2(double t)
114{
115 double result = 0.;
116 t=t-10.;
117 if(t>0) {
118 result = 2000.*(0.00181928*exp(-t/3.)-0.0147059*exp(-t/20.)+0.0128866*exp(-t/100.));
119 }
120 return result;
121}
122
123double E_branch2(double t)
124{
125 double result = 0.;
126 t=t-5;
127 if(t>0) {
128 result = 0.000627357*(1358.7*exp(-t*0.0385647)*t+1358.7*exp(-t*0.0114353)*t+100164.*exp(-t*0.0385647)-100164.*exp(-0.0114353*t));
129 }
130 return result;
131}
132
133TH1F Convolution_Tbranch(const double *Input_Curr_plot_001,double xmin,double xmax,int Npx=2000) // return TIGER output
134{
135
136 TH1F h_signalT("h_signalT","",Npx,xmin,xmax);
137 std::vector<double> Input_Time_plot_001(Npx);
138 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signalT.GetBinCenter(i+1);
139
140 double xmin_f1(0.), xmax_f1(1000);
141 int nStep_f1(200);
142 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
143 double x_f1(0.0);
144
145 std::vector<double> x(Npx);
146 std::vector<double> y_001(Npx);
147
148 for(int i=0; i<Npx; i++)
149 {
150 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
151 double fun_001=0.0;
152
153 for(int j=0; j<nStep_f1; j++)
154 {
155 x_f1=xmin_f1+step_f1*(j+0.5);
156 fun_001+=rectangle2(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*T_branch2(x[i]-x_f1)*step_f1;
157
158 }
159 y_001[i]=fun_001;
160 }
161
162 for(int init=0; init<Npx; init++){
163 h_signalT.SetBinContent(init+1,-y_001[init]);
164 }
165
166 return h_signalT;
167}
168
169TH1F Convolution_Tbranch(const double *Input_Curr_plot_001) // return TIGER output
170{
171 double xmin(0), xmax(1000);
172 const int Npx(2000);
173
174 TH1F h_signalT("h_signalT","",Npx,xmin,xmax);
175 double Input_Time_plot_001[Npx];
176 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signalT.GetBinCenter(i+1);
177
178 double xmin_f1(0.), xmax_f1(1000);
179 int nStep_f1(200);
180 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
181 double x_f1(0.0);
182
183 double x[Npx];
184 double y_001[Npx];
185
186 for(int i=0; i<Npx; i++)
187 {
188 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
189 double fun_001=0.0;
190
191 for(int j=0; j<nStep_f1; j++)
192 {
193 x_f1=xmin_f1+step_f1*(j+0.5);
194 fun_001+=rectangle2_fast(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*T_branch2(x[i]-x_f1)*step_f1;
195
196 }
197 y_001[i]=fun_001;
198 }
199
200 for(int init=0; init<Npx; init++){
201 h_signalT.SetBinContent(init+1,-y_001[init]);
202 }
203
204 return h_signalT;
205}
206
207TH1F Convolution_Ebranch(const double *Input_Curr_plot_001,double xmin,double xmax,int Npx=2000) // return TIGER output
208{
209 TH1F h_signalE("h_signalE","",Npx,xmin,xmax);
210 std::vector<double> Input_Time_plot_001(Npx);
211 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signalE.GetBinCenter(i+1);
212
213 double xmin_f1(0.), xmax_f1(1000);
214 int nStep_f1(200);
215 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
216 double x_f1(0.0);
217
218 std::vector<double> x(Npx);
219 std::vector<double> y_001(Npx);
220
221 for(int i=0; i<Npx; i++)
222 {
223 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
224 double fun_001=0.0;
225
226 for(int j=0; j<nStep_f1; j++)
227 {
228 x_f1=xmin_f1+step_f1*(j+0.5);
229 fun_001+=rectangle2(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*E_branch2(x[i]-x_f1)*step_f1;
230
231 }
232 y_001[i]=fun_001;
233 }
234
235 for(int init=0; init<Npx; init++){
236 h_signalE.SetBinContent(init+1,-y_001[init]);
237 }
238
239 return h_signalE;
240}
241
242TH1F Convolution_Ebranch(const double *Input_Curr_plot_001) // return TIGER output
243{
244 double xmin(0), xmax(1000);
245 const int Npx(2000);
246
247 TH1F h_signalE("h_signalE","",Npx,xmin,xmax);
248 double Input_Time_plot_001[Npx];
249 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signalE.GetBinCenter(i+1);
250
251 double xmin_f1(0.), xmax_f1(1000);
252 int nStep_f1(200);
253 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
254 double x_f1(0.0);
255
256 double x[Npx];
257 double y_001[Npx];
258
259 for(int i=0; i<Npx; i++)
260 {
261 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
262 double fun_001=0.0;
263
264 for(int j=0; j<nStep_f1; j++)
265 {
266 x_f1=xmin_f1+step_f1*(j+0.5);
267 fun_001+=rectangle2_fast(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*E_branch2(x[i]-x_f1)*step_f1;
268
269 }
270 y_001[i]=fun_001;
271 }
272
273 for(int init=0; init<Npx; init++){
274 h_signalE.SetBinContent(init+1,-y_001[init]);
275 }
276
277 return h_signalE;
278}
279
280#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
281/*
282 * @brief works exactly like Convolution_Tbranch, but works on every point rather than 1 in 10.
283 *
284 * @param Input_Curr_plot_001
285 * @return TH1F
286 */
287TH1F Convolution_Tbranch_noskip(const double *Input_Curr_plot_001) // return TIGER output
288{
289 double xmin(0), xmax(1000);
290 const int Npx(2000);
291
292 TH1F h_signalT("h_signalT", "", Npx, xmin, xmax);
293 double Input_Time_plot_001[Npx];
294 for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalT.GetBinCenter(i + 1);
295
296 double xmin_f1(0.), xmax_f1(1000);
297 int nStep_f1(2000);
298 double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
299 double x_f1(0.0);
300
301 double x[Npx];
302 double y_001[Npx];
303
304 for (int i = 0; i < Npx; i++) {
305 x[i] = xmin + (xmax - xmin) / Npx * (i + 0.5);
306 double fun_001 = 0.0;
307
308 for (int j = 0; j < nStep_f1; j++) {
309 x_f1 = xmin_f1 + step_f1 * (j + 0.5);
310 fun_001 += rectangle2_fast(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * T_branch2(x[i] - x_f1) * step_f1;
311 }
312 y_001[i] = fun_001;
313 }
314
315 for (int init = 0; init < Npx; init++) {
316 h_signalT.SetBinContent(init + 1, -y_001[init]);
317 }
318
319 return h_signalT;
320}
321
322TH1F Convolution_Ebranch_noskip(const double *Input_Curr_plot_001) // return TIGER output
323{
324 double xmin(0), xmax(1000);
325 const int Npx(2000);
326
327 TH1F h_signalE("h_signalE","",Npx,xmin,xmax);
328 double Input_Time_plot_001[Npx];
329 for(int i=0; i<Npx; i++) Input_Time_plot_001[i]=h_signalE.GetBinCenter(i+1);
330
331 double xmin_f1(0.), xmax_f1(1000);
332 int nStep_f1(2000);
333 double step_f1=(xmax_f1-xmin_f1)/nStep_f1;
334 double x_f1(0.0);
335
336 double x[Npx];
337 double y_001[Npx];
338
339 for(int i=0; i<Npx; i++)
340 {
341 x[i]=xmin+(xmax-xmin)/Npx*(i+0.5);
342 double fun_001=0.0;
343
344 for(int j=0; j<nStep_f1; j++)
345 {
346 x_f1=xmin_f1+step_f1*(j+0.5);
347 fun_001+=rectangle2_fast(x_f1,Input_Time_plot_001,Input_Curr_plot_001,Npx)*E_branch2(x[i]-x_f1)*step_f1;
348
349 }
350 y_001[i]=fun_001;
351 }
352
353 for(int init=0; init<Npx; init++){
354 h_signalE.SetBinContent(init+1,-y_001[init]);
355 }
356
357 return h_signalE;
358}
359#endif
360CTF::CTF(){
361 double xmin(0), xmax(1000);
362 const int Npx(2000);
363 TBranch.resize(Npx);// TBranch at x<0 is 0
364 double *p_TBranch=&(TBranch.front());
365 for (int i = 0; i < Npx; i++) {
366 double x=xmin+(xmax-xmin)*i/Npx;
367 TBranch[i]=T_branch2(x);
368 }// set initial TBranch values.
369 // basically, a N and a (2N-1) needs (3N-2) bins to avoid overlap.
370 conv=new ConvolveWithConst( p_TBranch,Npx,Npx);
371 //TBranch_fft.resize(Npx*3-2);
372
373 //conv.convolve
374}
375CTF::~CTF(){
376 delete conv;
377}
378/**
379 * @brief return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Tbranch
380 * has a option to test correctness.
381 *
382 * @param Input_Curr_plot_001
383 * @return TH1D
384 */
385TH1D CTF::Convolution_Tbranch_fft(const double *Input_Curr_plot_001) // return TIGER output; note it returns DOUBLE typed hist, rather than original float.
386{
387 double xmin(0), xmax(1000);
388 const int Npx(2000);
389
390 TH1D h_signal("h_signal", "", Npx, xmin, xmax);
391 double *p_signal=h_signal.GetArray()+1;//considering underflow.
392 this->conv->convolve(p_signal,Input_Curr_plot_001,0,Npx,-0.5);// the '-' before 0.5: from original function.
393
394#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
395 const double dx_abs_threhold=1e-6;// think about the accurancy of float; warn if |x| > threhold
396 const double dx_ref_threhold=1e-6;//think about the accurancy of float; warn if |dx/x| > threhold
397 TH1F h_signal_ref=Convolution_Tbranch_noskip(Input_Curr_plot_001);
398 bool found_mismatch=false;
399 for (int i=1; i<=Npx; i++){
400 if (abs((h_signal[i]-h_signal_ref[i])/h_signal_ref[i]) > dx_ref_threhold
401 and abs(h_signal[i]-h_signal_ref[i]) > dx_abs_threhold ) {
402 if (not found_mismatch){
403 std::cout<<"CgemDigitizerSvc::InductionGar::Convolution_Tbranch_fft found results not match: ";
404 found_mismatch=true;
405 }
406 std::cout<<"ref: "<<h_signal_ref[i]<<"; this: "<<h_signal[i]<<"; at i="<<i<<";threhold abs="<<dx_abs_threhold<<"; ref="<<dx_ref_threhold<<std::endl;
407 }
408 }
409#endif
410
411 return h_signal;
412}
413
414
415CEF::CEF(){
416 double xmin(0), xmax(1000);
417 const int Npx(2000);
418 EBranch.resize(Npx);// EBranch at x<0 is 0
419 double *p_EBranch=&(EBranch.front());
420 for (int i = 0; i < Npx; i++) {
421 double x=xmin+(xmax-xmin)*i/Npx;
422 EBranch[i]=E_branch2(x);
423 }// set initial EBranch values.
424 // basically, a N and a (2N-1) needs (3N-2) bins to avoid overlap.
425 conv=new ConvolveWithConst( p_EBranch,Npx,Npx);
426 //EBranch_fft.resize(Npx*3-2);
427
428 //conv.convolve
429}
430CEF::~CEF(){
431 delete conv;
432}
433/**
434 * @brief return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Ebranch
435 * has a option to test correctness.
436 *
437 * @param Input_Curr_plot_001
438 * @return TH1D
439 */
440TH1D CEF::Convolution_Ebranch_fft(const double *Input_Curr_plot_001) // return TIGER output; note it returns DOUBLE typed hist, rather than original float.
441{
442 double xmin(0), xmax(1000);
443 const int Npx(2000);
444
445 TH1D h_signal("h_signalE", "", Npx, xmin, xmax);
446 double *p_signal=h_signal.GetArray()+1;//considering underflow.
447 this->conv->convolve(p_signal,Input_Curr_plot_001,0,Npx,-0.5);
448
449#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
450 const double dx_abs_threhold=1e-6;// think about the accurancy of float; warn if |x| > threhold
451 const double dx_ref_threhold=1e-6;//think about the accurancy of float; warn if |dx/x| > threhold
452 TH1F h_signal_ref=Convolution_Ebranch_noskip(Input_Curr_plot_001);
453 bool found_mismatch=false;
454 for (int i=1; i<=Npx; i++){
455 if (abs((h_signal[i]-h_signal_ref[i])/h_signal_ref[i]) > dx_ref_threhold
456 and abs(h_signal[i]-h_signal_ref[i]) > dx_abs_threhold ) {
457 if (not found_mismatch){
458 std::cout<<"CgemDigitizerSvc::InductionGar::Convolution_Ebranch_fft found results not match: ";
459 found_mismatch=true;
460 }
461 std::cout<<"ref: "<<h_signal_ref[i]<<"; this: "<<h_signal[i]<<"; at i="<<i<<";threhold abs="<<dx_abs_threhold<<"; ref="<<dx_ref_threhold<<std::endl;
462 }
463 }
464#endif
465
466 return h_signal;
467}
468
469
470
471InductionGar::InductionGar(){
472 string testPath = getenv("CGEMDIGITIZERSVCROOT");
473 m_LUTFilePath = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_10_to_17.root";
474 m_V2EfineFile = testPath +"testFIle/V2Efine.root";
475 sample_delay = 162.5;
476 storeFlag = false;
477 m_saturation=true;
478
479 m_Ngaps_microSector=40;
480 m_gapShift_microSector[0][0]=0.;
481 m_gapShift_microSector[0][1]=0.;
482 m_gapShift_microSector[1][0]=0.;
483 m_gapShift_microSector[1][1]=0.;
484 m_gapShift_microSector[2][0]=0.;
485 m_gapShift_microSector[2][1]=0.;
486 m_microSector_width[0][0]=2*M_PI/m_Ngaps_microSector;
487 m_microSector_width[0][1]=2*M_PI/m_Ngaps_microSector;
488 m_microSector_width[1][0]=2*M_PI/m_Ngaps_microSector/2;
489 m_microSector_width[1][1]=2*M_PI/m_Ngaps_microSector/2;
490 m_microSector_width[2][0]=2*M_PI/m_Ngaps_microSector/2;
491 m_microSector_width[2][1]=2*M_PI/m_Ngaps_microSector/2;
492 m_QinGausSigma[0]=2;
493 m_QinGausSigma[1]=2;
494 m_gap_microSector=1.1;//in mm
495}
496
499
500void InductionGar::init(ICgemGeomSvc* geomSvc, double magConfig){
501 m_geomSvc = geomSvc;
502 m_magConfig = magConfig;
503 std::cout<<"InductionGar sampling time : "<<sample_delay<<std::endl;
504
505 string filePath_ratio = getenv("CGEMDIGITIZERSVCROOT");
506 string fileName;
507 if(0==m_magConfig) fileName = filePath_ratio + "/dat/par_InductionGar_0T.txt";
508 else fileName = filePath_ratio + "/dat/par_InductionGar_1T.txt";
509 ifstream fin(fileName.c_str(), ios::in);
510 cout << "InductionGar fileName: " << fileName << endl;
511
512 string strline;
513 string strpar;
514 if( ! fin.is_open() ){
515 cout << "ERROR: can not open file " << fileName << endl;
516 }
517
518 for(int m=0; m<3; m++){
519 for(int l=0; l<5; l++){
520 for(int k=0; k<5; k++){
521 for(int i=0; i<4; i++){
522 fin>>RatioX[m][l][k][i];
523 }
524 for(int j=0; j<3; j++){
525 fin>>RatioV[m][l][k][j];
526 }
527 }
528 }
529 }
530 fin.close();
531
532 string fileName1 = filePath_ratio + "/dat/VQ_relation.txt";
533 ifstream VQin(fileName1.c_str(), ios::in);
534 VQin>>VQ_slope>>VQ_const;
535 VQin.close();
536 //cout<<VQ_slope<<" "<<VQ_const<<endl;
537
538 string filePath = getenv("CGEMDIGITIZERSVCROOT");
539 for(int i=0; i<3; i++){ // layer
540 for(int j=0; j<5; j++){ // Grid-x
541 for(int k=0; k<5; k++){ // Grid-v
542 char temp1[1]; sprintf(temp1, "%i", (i+1));
543 char temp2[2]; sprintf(temp2, "%i", ((j+1)*10+k+1));
544 if(0==m_magConfig) fileName = filePath + "/dat/InducedCurrent_Root_0T/layer" + temp1 + "_" + temp2 + ".root";
545 else fileName = filePath + "/dat/InducedCurrent_Root_1T/layer" + temp1 + "_" + temp2 + ".root";
546 TFile* file_00 = TFile::Open(fileName.c_str(),"read");
547 TH1* readThis_x3 = 0;
548 TH1* readThis_x4 = 0;
549 TH1* readThis_x5 = 0;
550 TH1* readThis_x6 = 0;
551 TH1* readThis_v5 = 0;
552 TH1* readThis_v6 = 0;
553 TH1* readThis_v7 = 0;
554 file_00->GetObject("readThis_x3", readThis_x3);
555 file_00->GetObject("readThis_x4", readThis_x4); file_00->GetObject("readThis_v5", readThis_v5);
556 file_00->GetObject("readThis_x5", readThis_x5); file_00->GetObject("readThis_v6", readThis_v6);
557 file_00->GetObject("readThis_x6", readThis_x6); file_00->GetObject("readThis_v7", readThis_v7);
558
559 double scaleX=m_ScaleSignalX;
560 double scaleV=(1-scaleX*(RatioX[i][j][k][0]+RatioX[i][j][k][1]+RatioX[i][j][k][2]+RatioX[i][j][k][3]))/(RatioV[i][j][k][0]+RatioV[i][j][k][1]+RatioV[i][j][k][2]);
561 for(int init = 0; init<100; init++){
562 SignalX[i][j][k][0][init] = scaleX*readThis_x3->GetBinContent(init+1);
563 SignalX[i][j][k][1][init] = scaleX*readThis_x4->GetBinContent(init+1);
564 SignalX[i][j][k][2][init] = scaleX*readThis_x5->GetBinContent(init+1);
565 SignalX[i][j][k][3][init] = scaleX*readThis_x6->GetBinContent(init+1);
566 SignalV[i][j][k][0][init] = scaleV*readThis_v5->GetBinContent(init+1);
567 SignalV[i][j][k][1][init] = scaleV*readThis_v6->GetBinContent(init+1);
568 SignalV[i][j][k][2][init] = scaleV*readThis_v7->GetBinContent(init+1);
569 }
570 double temp[7][200];
571 for (int init = 0; init < 100; init++) { // since it was used like 2 bins eventually.
572 temp[0][init*2]=SignalX[i][j][k][0][init];
573 temp[1][init*2]=SignalX[i][j][k][1][init];
574 temp[2][init*2]=SignalX[i][j][k][2][init];
575 temp[3][init*2]=SignalX[i][j][k][3][init];
576 temp[4][init*2]=SignalV[i][j][k][0][init];
577 temp[5][init*2]=SignalV[i][j][k][1][init];
578 temp[6][init*2]=SignalV[i][j][k][2][init];
579 temp[0][init*2+1]=SignalX[i][j][k][0][init];
580 temp[1][init*2+1]=SignalX[i][j][k][1][init];
581 temp[2][init*2+1]=SignalX[i][j][k][2][init];
582 temp[3][init*2+1]=SignalX[i][j][k][3][init];
583 temp[4][init*2+1]=SignalV[i][j][k][0][init];
584 temp[5][init*2+1]=SignalV[i][j][k][1][init];
585 temp[6][init*2+1]=SignalV[i][j][k][2][init];
586 }
587 // since we actually used only index 1~80.
588 Signal_ConvolverX[i][j][k][0].init(temp[0]+2,160,2000);
589 Signal_ConvolverX[i][j][k][1].init(temp[1]+2,160,2000);
590 Signal_ConvolverX[i][j][k][2].init(temp[2]+2,160,2000);
591 Signal_ConvolverX[i][j][k][3].init(temp[3]+2,160,2000);
592 Signal_ConvolverV[i][j][k][0].init(temp[4]+2,160,2000);
593 Signal_ConvolverV[i][j][k][1].init(temp[5]+2,160,2000);
594 Signal_ConvolverV[i][j][k][2].init(temp[6]+2,160,2000);
595
596 }
597 }
598 }
599
600
601 TFile *LUTFile = TFile::Open(m_LUTFilePath.c_str(),"read");
602 TTree *tree = (TTree*)LUTFile->Get("tree");
603 Float_t V_th_T,V_th_E,QDC_a,QDC_b,QDC_saturation;
604 Int_t layer,sheet,strip_v,strip_x;
605 tree->SetBranchAddress("strip_x_boss", &strip_x);
606 tree->SetBranchAddress("strip_v_boss", &strip_v);
607 tree->SetBranchAddress("layer", &layer);
608 tree->SetBranchAddress("sheet", &sheet);
609 tree->SetBranchAddress("calib_QDC_slope", &QDC_a);
610 tree->SetBranchAddress("calib_QDC_const", &QDC_b);
611 tree->SetBranchAddress("calib_QDC_saturation", &QDC_saturation);
612 tree->SetBranchAddress("v_thr_T_mV", &V_th_T);
613 tree->SetBranchAddress("v_thr_E_mV", &V_th_E);
614
615 double thresholdMin_X_T=999;
616 double thresholdMin_V_T=999;
617 double thresholdMin_X_Q=999;
618 double thresholdMin_V_Q=999;
619 for(int i=0;i<tree->GetEntries();i++){
620 tree->GetEntry(i);
621 if(strip_x!=-1){
622 T_thr_V[layer][sheet][0][strip_x] = V_th_T;
623 if(V_th_T>0&&V_th_T<thresholdMin_X_T) thresholdMin_X_T=V_th_T;
624 //if(V_th_T<0) T_thr_V[layer][sheet][0][strip_x] =12.;
625 E_thr_V[layer][sheet][0][strip_x] = V_th_E;
626 if(V_th_E>0&&V_th_E<thresholdMin_X_Q) thresholdMin_X_Q=V_th_E;
627 //if(V_th_E<0) E_thr_V[layer][sheet][0][strip_x] =12.;
628 QDC_slope[layer][sheet][0][strip_x] = QDC_a;
629 QDC_const[layer][sheet][0][strip_x] = QDC_b;
630 Qsaturation[layer][sheet][0][strip_x] = QDC_saturation;
631 }
632 if(strip_v!=-1){
633 T_thr_V[layer][sheet][1][strip_v] = V_th_T;
634 if(V_th_T>0&&V_th_T<thresholdMin_V_T) thresholdMin_V_T=V_th_T;
635 //if(V_th_T<0) T_thr_V[layer][sheet][1][strip_v] =12.;
636 E_thr_V[layer][sheet][1][strip_v] = V_th_E;
637 if(V_th_E>0&&V_th_E<thresholdMin_V_Q) thresholdMin_V_Q=V_th_E;
638 //if(V_th_E<0) E_thr_V[layer][sheet][1][strip_v] =12.;
639 QDC_slope[layer][sheet][1][strip_v] = QDC_a;
640 QDC_const[layer][sheet][1][strip_v] = QDC_b;
641 Qsaturation[layer][sheet][1][strip_v] = QDC_saturation;
642 }
643 }
644 for(int i=0;i<tree->GetEntries();i++){
645 tree->GetEntry(i);
646 if(strip_x!=-1){
647 if(V_th_T<=0) T_thr_V[layer][sheet][0][strip_x] = thresholdMin_X_T;
648 if(V_th_E<=0) E_thr_V[layer][sheet][0][strip_x] = thresholdMin_X_Q;
649 }
650 if(strip_v!=-1){
651 if(V_th_T<=0) T_thr_V[layer][sheet][1][strip_v] = thresholdMin_V_T;
652 if(V_th_E<=0) E_thr_V[layer][sheet][1][strip_v] = thresholdMin_V_Q;
653 }
654 }
655
656 //3rd layer
657 for(int i=0;i<832;i++){
658 T_thr_V[2][0][0][i] = 12.;
659 E_thr_V[2][0][0][i] = 12.;
660 QDC_slope[2][0][0][i] = -10.47;
661 QDC_const[2][0][0][i] = 482.3;
662 Qsaturation[2][0][0][i] = 45.;
663 T_thr_V[2][1][0][i] = 10.;
664 E_thr_V[2][1][0][i] = 10.;
665 QDC_slope[2][1][0][i] = -10.47;
666 QDC_const[2][1][0][i] = 482.3;
667 Qsaturation[2][1][0][i] = 45.;
668 }
669
670 for(int i=0;i<1395;i++){
671 T_thr_V[2][0][1][i] = 12.;
672 E_thr_V[2][0][1][i] = 12.;
673 QDC_slope[2][0][1][i] = -10.47;
674 QDC_const[2][0][1][i] = 482.3;
675 Qsaturation[2][0][1][i] = 45.;
676 T_thr_V[2][1][1][i] = 10.;
677 E_thr_V[2][1][1][i] = 10.;
678 QDC_slope[2][1][1][i] = -10.47;
679 QDC_const[2][1][1][i] = 482.3;
680 Qsaturation[2][1][1][i] = 45.;
681 }
682
683 LUTFile->Close();
684 /*
685 TFile *V_Efine_File = new TFile(m_V2EfineFile.c_str());
686 TTree *tree_V2E = (TTree*)V_Efine_File->Get("tree");
687 Double_t V_ref,V2E_slope,V2E_const;
688 Int_t layer_,sheet_,stripv,stripx;
689 tree_V2E->SetBranchAddress("strip_x", &stripx);
690 tree_V2E->SetBranchAddress("strip_v", &stripv);
691 tree_V2E->SetBranchAddress("layer", &layer_);
692 tree_V2E->SetBranchAddress("sheet", &sheet_);
693 tree_V2E->SetBranchAddress("V2E_slope", &V2E_slope);
694 tree_V2E->SetBranchAddress("V2E_const", &V2E_const);
695 tree_V2E->SetBranchAddress("V_ref", &V_ref);
696
697
698 for(int i=0;i<tree_V2E->GetEntries();i++){
699 tree_V2E->GetEntry(i);
700 if(stripx!=-1){
701 Vref[layer_][sheet_][0][stripx] = V_ref;
702 toE_slope[layer_][sheet_][0][stripx] = V2E_slope;
703 toE_const[layer_][sheet_][0][stripx] = V2E_const;
704 }
705 if(stripv!=-1){
706 Vref[layer_][sheet_][1][stripv] = V_ref;
707 toE_slope[layer_][sheet_][1][stripv] = V2E_slope;
708 toE_const[layer_][sheet_][1][stripv] = V2E_const;
709 }
710 }
711
712 V_Efine_File->Close();
713 */
714
715 cout<<"InductionGar: "<<endl;
716 cout<<"m_Ngaps_microSector="<<m_Ngaps_microSector<<endl;
717 cout<<"m_gap_microSector="<<m_gap_microSector<<endl;
718 cout<<"m_gapShift_microSector: "<<m_gapShift_microSector[0][0]
719 <<", "<<m_gapShift_microSector[0][1]
720 <<", "<<m_gapShift_microSector[1][0]
721 <<", "<<m_gapShift_microSector[1][1]
722 <<", "<<m_gapShift_microSector[2][0]
723 <<", "<<m_gapShift_microSector[2][1]
724 <<endl;
725 cout<<"microSector_width="<<m_microSector_width[0]
726 <<", "<<m_microSector_width[1]
727 <<", "<<m_microSector_width[2]
728 <<", "<<m_microSector_width[3]
729 <<", "<<m_microSector_width[4]
730 <<", "<<m_microSector_width[5]
731 <<endl;
732 cout<<"QinGausSigma="<<m_QinGausSigma[0]
733 <<", "<<m_QinGausSigma[1]
734 <<endl;
735 //m_deadStrip[3][2]
736 m_deadStrip[0][0][0].insert(40);//x_0_down
737 m_deadStrip[0][0][0].insert(189);
738 m_deadStrip[0][0][0].insert(322);
739 m_deadStrip[0][0][0].insert(350);
740 m_deadStrip[0][0][0].insert(457);//x_0_up
741 m_deadStrip[0][0][1].insert(282);//v_0_down
742 m_deadStrip[0][0][1].insert(467);
743 m_deadStrip[0][0][1].insert(715);
744 m_deadStrip[0][0][1].insert(773);
745 m_deadStrip[0][0][1].insert(795);
746 m_deadStrip[0][0][1].insert(803);
747 m_deadStrip[1][0][1].insert(305);//v_1_down
748 m_deadStrip[1][0][1].insert(307);
749 m_deadStrip[1][0][1].insert(381);
750 m_deadStrip[1][0][1].insert(443);
751 m_deadStrip[1][0][1].insert(486);
752 m_deadStrip[1][0][1].insert(550);
753 m_deadStrip[1][0][1].insert(620);
754 m_deadStrip[1][0][1].insert(631);
755 m_deadStrip[1][0][1].insert(660);
756 m_deadStrip[1][0][1].insert(681);
757 m_deadStrip[1][1][1].insert(425);//v_1_up
758 m_deadStrip[1][1][1].insert(455);
759 m_deadStrip[1][1][1].insert(459);
760 m_deadStrip[1][1][1].insert(461);
761 m_deadStrip[1][1][1].insert(535);
762 m_deadStrip[1][1][1].insert(536);
763 m_deadStrip[1][1][1].insert(614);
764 m_deadStrip[1][1][1].insert(672);
765}
766
767void InductionGar::setMultiElectrons(int layer, int nElectrons, const std::vector<Float_t>& x, const std::vector<Float_t>& y, const std::vector<Float_t> &z, const std::vector<Float_t> &t){
768
769 m_xstripSheet.clear();
770 m_xstripID.clear();
771 m_vstripSheet.clear();
772 m_vstripID.clear();
773 m_xstripQ.clear();
774 m_vstripQ.clear();
775 m_xstripT_Branch.clear();
776 m_vstripT_Branch.clear();
777 m_xstripQ_Branch.clear();
778 m_vstripQ_Branch.clear();
779 m_xTfirst.clear();
780 m_vTfirst.clear();
781
782 int m000_nXstrips;
783 int m000_nVstrips;
784 std::vector<int> m000_xstripSheet;
785 std::vector<int> m000_xstripID;
786 std::vector<int> m000_vstripSheet;
787 std::vector<int> m000_vstripID;
788 std::vector<double> m000_xstripQ;
789 std::vector<double> m000_vstripQ;
790 std::vector<double> m000_xstripT_Branch;
791 std::vector<double> m000_vstripT_Branch;
792 std::vector<double> m000_xstripQ_Branch;
793 std::vector<double> m000_vstripQ_Branch;
794
795 m000_xstripSheet.clear();
796 m000_xstripID.clear();
797 m000_vstripSheet.clear();
798 m000_vstripID.clear();
799 m000_xstripQ.clear();
800 m000_vstripQ.clear();
801 m000_xstripT_Branch.clear();
802 m000_vstripT_Branch.clear();
803 m000_xstripQ_Branch.clear();
804 m000_vstripQ_Branch.clear();
805
806 CgemGeoLayer* CgemLayer;
807 CgemGeoReadoutPlane* CgemReadoutPlane;
808 CgemLayer = m_geomSvc->getCgemLayer(layer);
809 int NSheet = CgemLayer->getNumberOfSheet();
810
811 Vint Save_Sheetx, Save_Sheetv;
812 Vint Save_FinalstripIDx, Save_FinalstripIDv;
813 Vint Save_Gridx, Save_Gridv;
814 Vint Save_IDx, Save_IDv;
815 Vint Save_Tx, Save_Tv;
816 Save_Tx.clear();
817 Save_Tv.clear();
818 Save_IDx.clear();
819 Save_IDv.clear();
820 Save_Gridx.clear();
821 Save_Gridv.clear();
822 Save_Sheetx.clear();
823 Save_Sheetv.clear();
824 Save_FinalstripIDx.clear();
825 Save_FinalstripIDv.clear();
826 for(int i=0; i<2; i++){
827 for(int k=0; k<2; k++){
828 m_mapQ[i][k].clear();
829 //m_mapT[i][k].clear();
830 }
831 }
832
833 //std::cout << "nElectrons = " << nElectrons << std::endl;
834 //if(nElectrons>10000000) std::cout << "big nElectrons = " << nElectrons << std::endl;
835 for(int i=0; i<nElectrons; i++){
836 //cout<<i<<endl;
837 double phi_electron = atan2(y[i], x[i]);
838 double z_electron = z[i];
839 if(inMicroSectorGap(phi_electron,layer)) {
840 //cout<<"skip electron of layer "<<layer<<" at phi = "<<phi_electron<<endl;
841 continue;// skip electrons in the micro-sector gaps
842 }
843 G4ThreeVector pos(x[i], y[i], z[i]);
844 for(int j=0; j<NSheet; j++)
845 {
846 CgemReadoutPlane = m_geomSvc->getReadoutPlane(layer, j);
847 if(CgemReadoutPlane->OnThePlane(phi_electron,z_electron))
848 {
849 int xID;
850 double distX = CgemReadoutPlane->getDist2ClosestXStripCenter(phi_electron, xID);
851 int vID;
852 double distV = CgemReadoutPlane->getDist2ClosestVStripCenter(pos, vID);
853
854 //std::cout << "-------------" << std::endl;
855 //std::cout << "xID, distX = " << xID << "," << distX << std::endl;
856 //std::cout << "vID, distV = " << vID << "," << distV << std::endl;
857
858 int grid_x = 1;
859 int grid_v = 1;
860
861 //
862 //--- Calculation of grid index-------------
863 //
864 double L_XPitch = CgemReadoutPlane->getXPitch();
865 double L_VPitch = CgemReadoutPlane->getVPitch();
866 grid_v = floor(distX/L_XPitch*5.+2.5);
867 grid_x = floor(distV/L_VPitch*5.+2.5);
868 //std::cout << i << " " << j << " " << grid_x << " " << grid_v << std::endl;
869 //std::cout << "L_XPitch,L_VPitch = " << L_XPitch << "," << L_VPitch << std::endl;
870 //std::cout << "grid_x, grid_v = " << grid_x << "," << grid_v << std::endl;
871
872 if(xID>=0 && vID>=0) {
873 map<int, double>::iterator itx = m_mapQ[j][0].find(xID);
874 map<int, double>::iterator itv = m_mapQ[j][1].find(vID);
875 if(RatioX[layer][grid_x][grid_v][2]!=0){
876 if(itx==m_mapQ[j][0].end())
877 {
878 m_mapQ[j][0][xID]=RatioX[layer][grid_x][grid_v][2];
879 Save_Gridx.push_back(grid_x*100+grid_v*10+2);
880 Save_IDx.push_back(j*10000+xID);
881 Save_Tx.push_back(t[i]);
882 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][2] << std::endl;
883 }
884 else {
885 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][2];
886 m_mapQ[j][0][xID]=Q;
887 Save_Gridx.push_back(grid_x*100+grid_v*10+2);
888 Save_IDx.push_back(j*10000+xID);
889 Save_Tx.push_back(t[i]);
890 //std::cout << i << " " << j << "x-strip : Add " << Q << std::endl;
891 }
892 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID);
893 if(result_Sheetx==Save_Sheetx.end()){
894 Save_Sheetx.push_back(j*10000+xID);
895 }
896 }
897 if(RatioV[layer][grid_x][grid_v][1]!=0){
898 if(itv==m_mapQ[j][1].end())
899 {
900 m_mapQ[j][1][vID]=RatioV[layer][grid_x][grid_v][1];
901 Save_Gridv.push_back(grid_x*100+grid_v*10+1);
902 Save_IDv.push_back(j*10000+vID);
903 Save_Tv.push_back(t[i]);
904 }
905 else {
906 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][1];
907 m_mapQ[j][1][vID]=Q;
908 Save_Gridv.push_back(grid_x*100+grid_v*10+1);
909 Save_IDv.push_back(j*10000+vID);
910 Save_Tv.push_back(t[i]);
911 }
912 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID);
913 if(result_Sheetv==Save_Sheetv.end()){
914 Save_Sheetv.push_back(j*10000+vID);
915 }
916 }
917 }
918
919 if((xID>=0 && vID>=0) && (xID+1)<CgemReadoutPlane->getNXstrips()) {
920 map<int, double>::iterator itx = m_mapQ[j][0].find(xID+1);
921 if(RatioX[layer][grid_x][grid_v][3]!=0){
922 if(itx==m_mapQ[j][0].end())
923 {
924 m_mapQ[j][0][xID+1]=RatioX[layer][grid_x][grid_v][3];
925 Save_Gridx.push_back(grid_x*100+grid_v*10+3);
926 Save_IDx.push_back(j*10000+xID+1);
927 Save_Tx.push_back(t[i]);
928 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][3] << std::endl;
929 }
930 else {
931 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][3];
932 m_mapQ[j][0][xID+1]=Q;
933 Save_Gridx.push_back(grid_x*100+grid_v*10+3);
934 Save_IDx.push_back(j*10000+xID+1);
935 Save_Tx.push_back(t[i]);
936 }
937
938 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID+1);
939 if(result_Sheetx==Save_Sheetx.end()){
940 Save_Sheetx.push_back(j*10000+xID+1);
941 }
942 }
943 }
944
945 if((xID>=0 && vID>=0) && (vID+1)<CgemReadoutPlane->getNVstrips()) {
946 map<int, double>::iterator itv = m_mapQ[j][1].find(vID+1);
947 if(RatioV[layer][grid_x][grid_v][2]!=0){
948 if(itv==m_mapQ[j][1].end())
949 {
950 m_mapQ[j][1][vID+1]=RatioV[layer][grid_x][grid_v][2];
951 Save_Gridv.push_back(grid_x*100+grid_v*10+2);
952 Save_IDv.push_back(j*10000+vID+1);
953 Save_Tv.push_back(t[i]);
954 }
955 else {
956 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][2];
957 m_mapQ[j][1][vID+1]=Q;
958 Save_Gridv.push_back(grid_x*100+grid_v*10+2);
959 Save_IDv.push_back(j*10000+vID+1);
960 Save_Tv.push_back(t[i]);
961 }
962
963 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID+1);
964 if(result_Sheetv==Save_Sheetv.end()){
965 Save_Sheetv.push_back(j*10000+vID+1);
966 }
967 }
968 }
969
970 if((xID>=0 && vID>=0) && (xID-2)>=0) {
971 map<int, double>::iterator itx = m_mapQ[j][0].find(xID-2);
972 if(RatioX[layer][grid_x][grid_v][0]!=0){
973 if(itx==m_mapQ[j][0].end())
974 {
975 m_mapQ[j][0][xID-2]=RatioX[layer][grid_x][grid_v][0];
976 Save_Gridx.push_back(grid_x*100+grid_v*10+0);
977 Save_IDx.push_back(j*10000+xID-2);
978 Save_Tx.push_back(t[i]);
979 //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][0] << std::endl;
980 }
981 else {
982 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][0];
983 m_mapQ[j][0][xID-2]=Q;
984 Save_Gridx.push_back(grid_x*100+grid_v*10+0);
985 Save_IDx.push_back(j*10000+xID-2);
986 Save_Tx.push_back(t[i]);
987 }
988
989 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID-2);
990 if(result_Sheetx==Save_Sheetx.end()){
991 Save_Sheetx.push_back(j*10000+xID-2);
992 }
993 }
994 }
995
996 if((xID>=0 && vID>=0) && (xID-1)>=0) {
997 map<int, double>::iterator itx = m_mapQ[j][0].find(xID-1);
998 if(RatioX[layer][grid_x][grid_v][1]!=0){
999 if(itx==m_mapQ[j][0].end())
1000 {
1001 m_mapQ[j][0][xID-1]=RatioX[layer][grid_x][grid_v][1];
1002 Save_Gridx.push_back(grid_x*100+grid_v*10+1);
1003 Save_IDx.push_back(j*10000+xID-1);
1004 Save_Tx.push_back(t[i]);
1005 }
1006 else {
1007 double Q = (itx->second) + RatioX[layer][grid_x][grid_v][1];
1008 m_mapQ[j][0][xID-1]=Q;
1009 Save_Gridx.push_back(grid_x*100+grid_v*10+1);
1010 Save_IDx.push_back(j*10000+xID-1);
1011 Save_Tx.push_back(t[i]);
1012 }
1013
1014 vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j*10000+xID-1);
1015 if(result_Sheetx==Save_Sheetx.end()){
1016 Save_Sheetx.push_back(j*10000+xID-1);
1017 }
1018 }
1019 }
1020
1021 if((xID>=0 && vID>=0) && (vID-1)>=0) {
1022 map<int, double>::iterator itv = m_mapQ[j][1].find(vID-1);
1023 if(RatioV[layer][grid_x][grid_v][0]!=0){
1024 if(itv==m_mapQ[j][1].end())
1025 {
1026 m_mapQ[j][1][vID-1]=RatioV[layer][grid_x][grid_v][0];
1027 Save_Gridv.push_back(grid_x*100+grid_v*10+0);
1028 Save_IDv.push_back(j*10000+vID-1);
1029 Save_Tv.push_back(t[i]);
1030 }
1031 else {
1032 double Q = (itv->second) + RatioV[layer][grid_x][grid_v][0];
1033 m_mapQ[j][1][vID-1]=Q;
1034 Save_Gridv.push_back(grid_x*100+grid_v*10+0);
1035 Save_IDv.push_back(j*10000+vID-1);
1036 Save_Tv.push_back(t[i]);
1037 }
1038
1039 vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j*10000+vID-1);
1040 if(result_Sheetv==Save_Sheetv.end()){
1041 Save_Sheetv.push_back(j*10000+vID-1);
1042 }
1043 }
1044 }
1045 }
1046 }
1047 }
1048 //std::cout << "loop electron finished " << std::endl;
1049
1050 //-----Q Threshold = 45ADC = 1.5fC = 1.5*1e-15 C --------
1051 //-----Q Satruation = 45fC
1052 //-----Q Resolution = 0.256fC
1053 //-----T jitter = 2ns
1054
1055 double Q_threshold = 1.5e-15/1.602176565e-19;
1056 double Q_saturation= 45.e-15/1.602176565e-19;
1057 double Q_resolution= 0.256e-15/1.602176565e-19;
1058 double T_threshold = 12; // 12 mV
1059 //double T_threshold = 24; // 24 mV
1060 //double T_threshold = 18; // 18 mV
1061
1062 // cout << "---------------------" << endl;
1063 /*
1064 // consider charge resolution and saturation
1065 for(int i=0; i<Save_Sheetx.size(); i++){
1066 int sheetx_id = Save_Sheetx[i]/10000;
1067 int stripx_id = Save_Sheetx[i]%10000;
1068 //cout << "before " << m_mapQ[sheetx_id][0][stripx_id] << endl;
1069 m_mapQ[sheetx_id][0][stripx_id]=gRandom->Gaus(m_mapQ[sheetx_id][0][stripx_id],Q_resolution);
1070 if(m_mapQ[sheetx_id][0][stripx_id]>Q_saturation) m_mapQ[sheetx_id][0][stripx_id] = Q_saturation;
1071 //cout << "after " << m_mapQ[sheetx_id][0][stripx_id] << endl;
1072 }
1073
1074 for(int i=0; i<Save_Sheetv.size(); i++){
1075 int sheetv_id = Save_Sheetv[i]/10000;
1076 int stripv_id = Save_Sheetv[i]%10000;
1077 //cout << "before " << m_mapQ[sheetv_id][1][stripv_id] << endl;
1078 m_mapQ[sheetv_id][1][stripv_id]=gRandom->Gaus(m_mapQ[sheetv_id][1][stripv_id],Q_resolution);
1079 if(m_mapQ[sheetv_id][1][stripv_id]>Q_saturation) m_mapQ[sheetv_id][1][stripv_id] = Q_saturation;
1080 //cout << "after " << m_mapQ[sheetv_id][1][stripv_id] << endl;
1081 }
1082
1083
1084 int Nx_below_threshold = 0;
1085 for(int i=0; i<Save_Sheetx.size(); i++){
1086 int sheetx_id = Save_Sheetx[i]/10000;
1087 int stripx_id = Save_Sheetx[i]%10000;
1088 //std::cout << m_mapQ[sheetx_id][0][stripx_id] << std::endl;
1089 if (m_mapQ[sheetx_id][0][stripx_id]<Q_threshold) Nx_below_threshold++;
1090 }
1091
1092 int Nv_below_threshold = 0;
1093 for(int i=0; i<Save_Sheetv.size(); i++){
1094 int sheetv_id = Save_Sheetv[i]/10000;
1095 int stripv_id = Save_Sheetv[i]%10000;
1096 //std::cout << m_mapQ[sheetv_id][1][stripv_id] << std::endl;
1097 if (m_mapQ[sheetv_id][1][stripv_id]<Q_threshold) Nv_below_threshold++;
1098 }
1099 */
1100 //-------------------------------------------------------------
1101
1102 //m000_nXstrips = m_mapQ[0][0].size() + m_mapQ[1][0].size() - Nx_below_threshold;
1103 //m000_nVstrips = m_mapQ[0][1].size() + m_mapQ[1][1].size() - Nv_below_threshold;
1104 m000_nXstrips = m_mapQ[0][0].size() + m_mapQ[1][0].size();
1105 m000_nVstrips = m_mapQ[0][1].size() + m_mapQ[1][1].size();
1106
1107 //cout << m000_nXstrips << " " << Nx_below_threshold << " " << m000_nVstrips << " " << Nv_below_threshold << endl;
1108
1109
1110 for(int i=0; i<Save_Sheetx.size(); i++){
1111 int sheetx_id = Save_Sheetx[i]/10000;
1112 int stripx_id = Save_Sheetx[i]%10000;
1113 //if(m_mapQ[sheetx_id][0][stripx_id]>=Q_threshold){
1114 Save_FinalstripIDx.push_back(Save_Sheetx[i]);
1115 //std::cout << "zhaojy check InductionGar sheet -> " << m000_xstripSheet.size() << " " << sheetx_id << " " << stripx_id << std::endl;
1116 m000_xstripSheet.push_back(sheetx_id);
1117 m000_xstripID.push_back(stripx_id);
1118 m000_xstripQ.push_back(m_mapQ[sheetx_id][0][stripx_id]);
1119 //}
1120 }
1121
1122 for(int i=0; i<Save_Sheetv.size(); i++){
1123 int sheetv_id = Save_Sheetv[i]/10000;
1124 int stripv_id = Save_Sheetv[i]%10000;
1125 //if(m_mapQ[sheetv_id][1][stripv_id]>=Q_threshold){
1126 Save_FinalstripIDv.push_back(Save_Sheetv[i]);
1127 m000_vstripSheet.push_back(sheetv_id);
1128 m000_vstripID.push_back(stripv_id);
1129 m000_vstripQ.push_back(m_mapQ[sheetv_id][1][stripv_id]);
1130 //}
1131 }
1132
1133 // string FilePath = getenv("TESTCGEMDIGISVCALGROOT");
1134 //
1135 // for(int h=0; h<Save_FinalstripIDx.size(); h++){
1136 // const int _low = 0;
1137 // const int _up = 1000;
1138 // const int _nbins = 2000;
1139 // double binsize = (double)(_up-_low)/_nbins;
1140 // int sheetx_id = Save_FinalstripIDx[h]/10000;
1141 // int stripx_id = Save_FinalstripIDx[h]%10000;
1142 // double histX_bin_Content[_nbins];
1143 // for(int i=0; i<_nbins; i++){
1144 // histX_bin_Content[i] = 0;
1145 // }
1146 // for(int i=0; i<Save_Gridx.size(); i++){
1147 // int z_grid_x = Save_Gridx[i]/100;
1148 // int z_grid_v = Save_Gridx[i]%100/10;
1149 // int z_index = Save_Gridx[i]%10;
1150 // int z_IDx = Save_IDx[i];
1151 // int start_bin = Save_Tx[i]/binsize;
1152 // std::cout << "start_bin x" << start_bin << std::endl;
1153 // if(z_IDx == Save_FinalstripIDx[h]){
1154 // for(int addi = start_bin; addi < start_bin+160; addi+=2){
1155 // histX_bin_Content[addi]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1156 // histX_bin_Content[addi+1]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1157 // }
1158 // }
1159 // }
1160 //
1161 // TH1F h_signal = Convolution_Tbranch(histX_bin_Content);
1162 // T_threshold = T_thr_V[layer][sheetx_id][0][stripx_id];
1163 // TH1D h_signal = ctf.Convolution_Tbranch_fft(histX_bin_Content);
1164 // TH1F h_signal = Convolution_Tbranch(histX_bin_Content);
1165 //
1166 // int flg_xT = 0;
1167 // for(int init=1; init<_nbins; init++){
1168 // if(flg_xT==0 && h_signal.GetBinContent(init+1)<=-T_threshold){
1169 // m000_xstripT_Branch.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
1170 // flg_xT=1;
1171 // }
1172 // if(flg_xT==0 && h_signal.GetBinContent(init+1)>=T_threshold){
1173 // m000_xstripT_Branch.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
1174 // flg_xT=1;
1175 // }
1176 // }
1177 // if(flg_xT==0) m000_xstripT_Branch.push_back(-99);
1178 //
1179 // char temp1[1]; sprintf(temp1, "%i", sheetx_id);
1180 // char temp2[3]; sprintf(temp2, "%i", stripx_id);
1181 // string FILEname = FilePath + "/share/output/X-" + temp1 + '-' + temp2 + ".root";
1182 // TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); h_signal.Write(); f_x1->Close(); delete f_x1;
1183 // }
1184 //
1185 // for(int h=0; h<Save_FinalstripIDv.size(); h++){
1186 // const int _low = 0;
1187 // const int _up = 1000;
1188 // const int _nbins = 2000;
1189 // double binsize = (double)(_up-_low)/_nbins;
1190 // int sheetv_id = Save_FinalstripIDv[h]/10000;
1191 // int stripv_id = Save_FinalstripIDv[h]%10000;
1192 // double histV_bin_Content[_nbins];
1193 // for(int i=0; i<_nbins; i++){
1194 // histV_bin_Content[i] = 0;
1195 // }
1196 // for(int i=0; i<Save_Gridv.size(); i++){
1197 // int z_grid_x = Save_Gridv[i]/100;
1198 // int z_grid_v = Save_Gridv[i]%100/10;
1199 // int z_index = Save_Gridv[i]%10;
1200 // int z_IDv = Save_IDv[i];
1201 // int start_bin = Save_Tv[i]/binsize;
1202 // std::cout << "start_bin v" << start_bin << std::endl;
1203 //if(z_IDv == Save_FinalstripIDv[h]){
1204 // for(int addi = start_bin; addi < start_bin+160; addi+=2){
1205 // histV_bin_Content[addi]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1206 // histV_bin_Content[addi+1]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1207 // }
1208 //}
1209 //}
1210 //
1211 //TH1F h_signal = Convolution_Tbranch(histV_bin_Content);
1212 //T_threshold = T_thr_V[layer][sheetv_id][1][stripv_id];
1213 //TH1F h_signal = Convolution_Tbranch(histV_bin_Content);
1214 //TH1D h_signal = ctf.Convolution_Tbranch_fft(histV_bin_Content);
1215 //
1216 //int flg_vT = 0;
1217 //for(int init=1; init<_nbins; init++){
1218 // if(flg_vT==0 && h_signal.GetBinContent(init+1)<=-T_threshold){
1219 // m000_vstripT_Branch.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
1220 // flg_vT=1;
1221 // }
1222 // if(flg_vT==0 && h_signal.GetBinContent(init+1)>=T_threshold){
1223 // m000_vstripT_Branch.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
1224 // flg_vT=1;
1225 // }
1226 //}
1227 //if(flg_vT==0) m000_vstripT_Branch.push_back(-99);
1228 //
1229 //char temp1[1]; sprintf(temp1, "%i", sheetv_id);
1230 //char temp2[3]; sprintf(temp2, "%i", stripv_id);
1231 //string FILEname = FilePath + "/share/output/V-" + temp1 + '-' + temp2 + ".root";
1232 //TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); h_signal.Write(); f_x1->Close(); delete f_x1;
1233 //}
1234 //std::cout << "loop Q finished " << std::endl;
1235 //string FilePath = getenv("TESTCGEMDIGISVCALGROOT");
1236 //std::cout<<"start signal "<<std::endl;
1237 double Tmin;
1238 std::vector<double> xTfirst;//time when first electron leave 3rd gem
1239 std::vector<double> vTfirst;
1240 //T_Branch & E_Branch
1241 for(int h=0; h<Save_FinalstripIDx.size(); h++){
1242 std::map<int,std::vector<int> > electron_hit_tmp;
1243 const int _low = 0;
1244 const int _up = 1000;
1245 const int _nbins = 2000;
1246 double binsize = (double)(_up-_low)/_nbins;
1247 int sheetx_id = Save_FinalstripIDx[h]/10000;
1248 int stripx_id = Save_FinalstripIDx[h]%10000;
1249 double histX_bin_Content[_nbins];
1250 for(int i=0; i<_nbins; i++){
1251 histX_bin_Content[i] = 0;
1252 }
1253
1254 Tmin = 999999;
1255 if (isDeadStrip(layer,sheetx_id,0,stripx_id)) {
1256 xTfirst.push_back(Tmin);
1257 double T_th = -99.;
1258 m000_xstripT_Branch.push_back(T_th);
1259 double Qin = -99.;
1260 m000_xstripQ_Branch.push_back(Qin);
1261 continue;
1262 };
1263
1264 for(int i=0; i<Save_Gridx.size(); i++){
1265 //int z_grid_x = Save_Gridx[i]/100;
1266 //int z_grid_v = Save_Gridx[i]%100/10;
1267 //int z_index = Save_Gridx[i]%10;
1268 int z_IDx = Save_IDx[i];
1269 int start_bin = Save_Tx[i]/binsize;
1270 // std::cout << "start_bin x" << start_bin << std::endl;
1271
1272 if(z_IDx == Save_FinalstripIDx[h]){
1273 if(Save_Tx[i]<Tmin) Tmin = Save_Tx[i];
1274 if (start_bin<_nbins and start_bin>=0)electron_hit_tmp[Save_Gridx[i]].push_back(start_bin); //prepare the hit list. should have boundary check
1275 //discards evt >= _nbins or evt < 0.
1276
1277 //for(int addi = start_bin; addi < start_bin+160; addi+=2){
1278 // histX_bin_Content[addi]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1279 // histX_bin_Content[addi+1]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1280 //}
1281 }
1282
1283 }
1284
1285 for (std::map<int, std::vector<int> >::iterator it = electron_hit_tmp.begin();
1286 it != electron_hit_tmp.end(); it++) {
1287 std::vector<int> &hitlist = it->second;
1288 const int &Save_Grid = it->first;
1289 // legacy
1290 if (hitlist.size() < electrons_select_method_threhold) {
1291 conv1PerGrid_legacy(Save_Grid, histX_bin_Content, hitlist, layer, 0);
1292 } else {
1293 conv1PerGrid_fft(Save_Grid, histX_bin_Content, hitlist, layer, 0);
1294 }
1295 }
1296
1297 xTfirst.push_back(Tmin);
1298 TH1F sig_current("current","",_nbins,_low,_up);
1299 for(int i=0;i<_nbins;i++){
1300 sig_current.SetBinContent(i+1,histX_bin_Content[i]);
1301 }
1302
1303 //TH1F h_signalT = Convolution_Tbranch(histX_bin_Content);
1304 TH1D h_signalT = ctf.Convolution_Tbranch_fft(histX_bin_Content);
1305
1306
1307 int flg_xT = 0;
1308 double T_th = 0.;
1309 T_threshold = T_thr_V[layer][sheetx_id][0][stripx_id];
1310 //if(T_threshold<=0) T_threshold=thresholdMin_X_T;
1311 //T_threshold = 46.7;
1312 for(int init=1; init<_nbins; init++){
1313 if(flg_xT==0 && fabs(h_signalT.GetBinContent(init+1))>=fabs(T_threshold)){
1314 T_th = _low+binsize*(init-1)+fabs(T_threshold-h_signalT.GetBinContent(init))*binsize/fabs(h_signalT.GetBinContent(init+1)-h_signalT.GetBinContent(init))+gRandom->Gaus(0,2);
1315 m000_xstripT_Branch.push_back(T_th); // jitter 2ns
1316 flg_xT=1;
1317 }
1318 }
1319 if(flg_xT==0) {
1320 T_th = -99.;
1321 m000_xstripT_Branch.push_back(T_th);
1322 }
1323 double V_sampling = 0.;
1324 double Qin = 0.;
1325 double Efine = 0.;
1326
1327
1328 if(T_th>0){
1329 //TH1D h_signalE = Convolution_Ebranch(histX_bin_Content);//debug; do not need to do if V_T(max) < T_threshold(mV)
1330 TH1D h_signalE = cef.Convolution_Ebranch_fft(histX_bin_Content);//debug; do not need to do if V_T(max) < T_threshold(mV)
1331 double T_sample = T_th+sample_delay;//ns
1332 int sample_bin = T_sample/binsize;
1333 double sample_bin_ = T_sample/binsize;
1334 V_sampling = h_signalE.GetBinContent(sample_bin)+(h_signalE.GetBinContent(sample_bin+1)-h_signalE.GetBinContent(sample_bin))*double(sample_bin_-sample_bin);
1335 int over_th = 0;
1336
1337 //for(int iBin = 1;iBin<=_nbins;iBin++){
1338 // if(h_signalE.GetBinContent(iBin)>E_thr_V[layer][sheetx_id][0][stripx_id]){
1339 // over_th = 1;
1340 // break;
1341 // }
1342 //}
1343 if(h_signalE.GetMaximum()>E_thr_V[layer][sheetx_id][0][stripx_id]) over_th=1;
1344 /*
1345 if(over_th==0){
1346 Qin = -999.;
1347 }
1348 else{
1349 Efine = toE_const[layer][sheetx_id][0][stripx_id]+V_sampling*toE_slope[layer][sheetx_id][0][stripx_id];//????
1350 if(Efine<-16){
1351 Qin = Qsaturation[layer][sheetx_id][0][stripx_id];
1352 }
1353 else if(Efine>=1008){
1354 Qin = 0.;
1355 }
1356 else {
1357 Qin = (Efine-QDC_const[layer][sheetx_id][0][stripx_id])/QDC_slope[layer][sheetx_id][0][stripx_id];
1358 }
1359 }
1360 */
1361 // m000_xstripQ_Branch.push_back(Qin);
1362 if(over_th==1) {
1363 Qin = V_sampling*VQ_slope + VQ_const+gRandom->Gaus(0,m_QinGausSigma[0]);
1364 if(m_saturation&&Qin>Qsaturation[layer][sheetx_id][0][stripx_id])
1365 {
1366 Qin = Qsaturation[layer][sheetx_id][0][stripx_id];
1367 }
1368 }
1369 else Qin = -99;
1370 }
1371 else{
1372 V_sampling = -99.;
1373 Qin = -99.;
1374 //m000_xstripQ_Branch.push_back(Qin);
1375 }
1376
1377 m000_xstripQ_Branch.push_back(Qin);
1378 //cout<<"X layer:"<<layer<<" sheet:"<<sheetx_id<<" strip:"<<stripx_id<<" T_Branch:"<<T_th<<" V_sample:"<<V_sampling<<endl;
1379 /*
1380 char temp0[1]; sprintf(temp0, "%i", layer);
1381 char temp1[1]; sprintf(temp1, "%i", sheetx_id);
1382 char temp2[3]; sprintf(temp2, "%i", stripx_id);
1383 string FILEname = FilePath + "/share/output/X-" + temp0 + '-' + temp1 + '-' + temp2 + ".root";
1384 TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); sig_current.Write(); h_signalT.Write(); f_x1->Close(); delete f_x1;
1385 */
1386 }
1387 //std::cout << "loop X_Q finished " << std::endl;
1388
1389 for(int h=0; h<Save_FinalstripIDv.size(); h++){
1390 std::map<int,std::vector<int> > electron_hit_tmp;
1391 const int _low = 0;
1392 const int _up = 1000;
1393 const int _nbins = 2000;
1394 double binsize = (double)(_up-_low)/_nbins;
1395 int sheetv_id = Save_FinalstripIDv[h]/10000;
1396 int stripv_id = Save_FinalstripIDv[h]%10000;
1397 double histV_bin_Content[_nbins];
1398 for(int i=0; i<_nbins; i++){
1399 histV_bin_Content[i] = 0;
1400 }
1401
1402 Tmin = 999999;
1403 if (isDeadStrip(layer,sheetv_id,1,stripv_id)) {
1404 vTfirst.push_back(Tmin);
1405 double T_th = -99.;
1406 m000_vstripT_Branch.push_back(T_th);
1407 double Qin = -99.;
1408 m000_vstripQ_Branch.push_back(Qin);
1409 continue;
1410 };
1411
1412 for(int i=0; i<Save_Gridv.size(); i++){
1413 //int z_grid_x = Save_Gridv[i]/100;
1414 //int z_grid_v = Save_Gridv[i]%100/10;
1415 //int z_index = Save_Gridv[i]%10;
1416 int z_IDv = Save_IDv[i];
1417 int start_bin = Save_Tv[i]/binsize;
1418
1419 //std::cout << "start_bin v" << start_bin << std::endl;
1420 if(z_IDv == Save_FinalstripIDv[h]){
1421 if(Save_Tv[i]<Tmin) Tmin = Save_Tv[i];
1422 if (start_bin<_nbins and start_bin>=0)
1423 electron_hit_tmp[Save_Gridv[i]].push_back(start_bin); //prepare the hit list.
1424
1425 //for(int addi = start_bin; addi < start_bin+160; addi+=2){
1426 // histV_bin_Content[addi]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1427 // histV_bin_Content[addi+1]+=SignalV[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
1428 //}
1429 }
1430 }
1431
1432
1433 for (std::map<int,std::vector<int> >::iterator it=electron_hit_tmp.begin();
1434 it !=electron_hit_tmp.end(); it++){
1435 std::vector<int> & hitlist=it->second;
1436 const int & Save_Grid=it->first;
1437 // legacy
1438 if ( hitlist.size()< electrons_select_method_threhold){
1439 conv1PerGrid_legacy(Save_Grid,histV_bin_Content,hitlist,layer,1);
1440 }
1441 else{// fft
1442 conv1PerGrid_fft(Save_Grid,histV_bin_Content,hitlist,layer,1);
1443 }
1444 }
1445
1446
1447
1448 vTfirst.push_back(Tmin);
1449 TH1F sig_current("current","",_nbins,_low,_up);
1450 for(int i=0;i<_nbins;i++){
1451 sig_current.SetBinContent(i+1,histV_bin_Content[i]);
1452 }
1453
1454 //TH1F h_signalT = Convolution_Tbranch(histV_bin_Content);
1455 TH1D h_signalT = ctf.Convolution_Tbranch_fft(histV_bin_Content);
1456
1457 int flg_vT = 0;
1458 double T_th = 0.;
1459 T_threshold = T_thr_V[layer][sheetv_id][1][stripv_id];
1460 //if(T_threshold<=0) T_threshold=thresholdMin_V_T;
1461 //T_threshold = 46.7;
1462 for(int init=1; init<_nbins; init++){
1463 if(flg_vT==0 && fabs(h_signalT.GetBinContent(init+1))>=fabs(T_threshold)){
1464 T_th = _low+binsize*(init-1)+fabs(T_threshold-h_signalT.GetBinContent(init))*binsize/fabs(h_signalT.GetBinContent(init+1)-h_signalT.GetBinContent(init))+gRandom->Gaus(0,2);
1465 m000_vstripT_Branch.push_back(T_th); // jitter 2ns
1466 flg_vT=1;
1467 }
1468 }
1469 if(flg_vT==0) {
1470 T_th = -99.;
1471 m000_vstripT_Branch.push_back(T_th);
1472 }
1473 double V_sampling = 0.;
1474 double Qin = 0.;
1475 double Efine = 0.;
1476 if(T_th>0){
1477 //TH1F h_signalE = Convolution_Ebranch(histV_bin_Content);
1478 TH1D h_signalE = cef.Convolution_Ebranch_fft(histV_bin_Content);
1479 double T_sample = T_th+sample_delay;//ns
1480 int sample_bin = T_sample/binsize;
1481 double sample_bin_ = T_sample/binsize;
1482 V_sampling = h_signalE.GetBinContent(sample_bin)+(h_signalE.GetBinContent(sample_bin+1)-h_signalE.GetBinContent(sample_bin))*double(sample_bin_-sample_bin);
1483 int over_th = 0;
1484 //for(int iBin = 1;iBin<=_nbins;iBin++){
1485 // if(h_signalE.GetBinContent(iBin)>E_thr_V[layer][sheetv_id][1][stripv_id]){
1486 // over_th = 1;
1487 // break;
1488 // }
1489 //}
1490 if(h_signalE.GetMaximum()>E_thr_V[layer][sheetv_id][1][stripv_id]) over_th=1;
1491 /*
1492 if(over_th==0){
1493 Qin = -999.;
1494 }
1495 else{
1496 Efine = toE_const[layer][sheetv_id][1][stripv_id]+V_sampling*toE_slope[layer][sheetv_id][1][stripv_id];//????
1497 if(Efine<-16){
1498 Qin = Qsaturation[layer][sheetv_id][1][stripv_id];
1499 }
1500 else if(Efine>=1008){
1501 Qin = 0.;
1502 }
1503 else {
1504 Qin = (Efine-QDC_const[layer][sheetv_id][1][stripv_id])/QDC_slope[layer][sheetv_id][1][stripv_id];
1505 }
1506 }
1507 */
1508 //m000_vstripQ_Branch.push_back(Qin);
1509 if(over_th==1) {
1510 Qin = V_sampling*VQ_slope + VQ_const+gRandom->Gaus(0,m_QinGausSigma[1]);
1511 if(m_saturation&&Qin>Qsaturation[layer][sheetv_id][1][stripv_id])
1512 Qin = Qsaturation[layer][sheetv_id][1][stripv_id];
1513 }
1514 else Qin = -99;
1515 }
1516 else{
1517 V_sampling = -99.;
1518 Qin = -99.;
1519 //m000_vstripQ_Branch.push_back(Qin);
1520 }
1521
1522 m000_vstripQ_Branch.push_back(Qin);
1523 // cout<<"V layer:"<<layer<<" sheet:"<<sheetv_id<<" strip:"<<stripv_id<<" T_Branch:"<<T_th<<" V_sample:"<<V_sampling<<endl;
1524 /*
1525 int layer_tem = layer;
1526 char temp0[1]; sprintf(temp0, "%i", layer_tem);
1527 char temp1[1]; sprintf(temp1, "%i", sheetv_id);
1528 char temp2[3]; sprintf(temp2, "%i", stripv_id);
1529 string FILEname = FilePath + "/share/output/V-" + temp0 + '-' + temp1 + '-' + temp2 + ".root";
1530 TFile* f_x1 = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); sig_current.Write(); h_signalT.Write(); f_x1->Close(); delete f_x1;
1531 */
1532 }
1533 //std::cout << "loop V_Q finished " << std::endl;
1534
1535 //-------------------------------------------------------------------------------
1536 //cout<<"push back"<<endl;
1537 double q_to_fC_factor = 1.602176565e-4;
1538 //std::cout << "loop V_Q finished " << std::endl;
1539
1540 //cout<<"X"<<endl;
1541 if(storeFlag){
1542 m_nXstrips = 0;
1543 for(int i=0; i<m000_nXstrips; i++){
1544 m_xstripSheet.push_back(m000_xstripSheet[i]);
1545 m_xstripID.push_back(m000_xstripID[i]);
1546 m_xstripQ.push_back(m000_xstripQ[i]*q_to_fC_factor);
1547 m_xstripT_Branch.push_back(m000_xstripT_Branch[i]);
1548 m_xstripQ_Branch.push_back(m000_xstripQ_Branch[i]);
1549 m_xTfirst.push_back(xTfirst[i]);
1550 m_nXstrips++;
1551 }
1552 m_nVstrips = 0;
1553 //cout<<"V"<<endl;
1554 for(int i=0; i<m000_nVstrips; i++){
1555 m_vstripSheet.push_back(m000_vstripSheet[i]);
1556 m_vstripID.push_back(m000_vstripID[i]);
1557 m_vstripQ.push_back(m000_vstripQ[i]*q_to_fC_factor);
1558 m_vstripT_Branch.push_back(m000_vstripT_Branch[i]);
1559 m_vstripQ_Branch.push_back(m000_vstripQ_Branch[i]);
1560 m_vTfirst.push_back(vTfirst[i]);
1561 m_nVstrips++;
1562 }
1563 }
1564 else{
1565 m_nXstrips = 0;
1566 for(int i=0; i<m000_nXstrips; i++){
1567 if(m000_xstripQ_Branch[i]>=0){
1568 m_xstripSheet.push_back(m000_xstripSheet[i]);
1569 m_xstripID.push_back(m000_xstripID[i]);
1570 m_xstripQ.push_back(m000_xstripQ[i]*q_to_fC_factor);
1571 m_xstripT_Branch.push_back(m000_xstripT_Branch[i]);
1572 m_xstripQ_Branch.push_back(m000_xstripQ_Branch[i]);
1573 m_xTfirst.push_back(xTfirst[i]);
1574 m_nXstrips++;
1575 // cout<<m000_xstripQ[i]*q_to_fC_factor<<" "<<m000_xstripQ_Branch[i]<<endl;
1576 }
1577 }
1578 m_nVstrips = 0;
1579 //cout<<"V"<<endl;
1580 for(int i=0; i<m000_nVstrips; i++){
1581 if(m000_vstripQ_Branch[i]>=0){
1582 m_vstripSheet.push_back(m000_vstripSheet[i]);
1583 m_vstripID.push_back(m000_vstripID[i]);
1584 m_vstripQ.push_back(m000_vstripQ[i]*q_to_fC_factor);
1585 m_vstripT_Branch.push_back(m000_vstripT_Branch[i]);
1586 m_vstripQ_Branch.push_back(m000_vstripQ_Branch[i]);
1587 m_vTfirst.push_back(vTfirst[i]);
1588 m_nVstrips++;
1589 // cout<<m000_vstripQ[i]*q_to_fC_factor<<" "<<m000_vstripQ_Branch[i]<<endl;
1590 }
1591 }
1592 }
1593
1594 //cout<<"InductionGar::setMultiElectrons() ends"<<endl;
1595
1596}
1597
1598
1599
1600//i may want to warp this into a function
1601/**
1602 * @brief legacy method of doing conv1 at a single grid for a strip. additive on hist_bin_Content
1603 * discards t > _nbins or t <0 events
1604 * @param Save_Grid [grid_x]*100+[grid_v]*10+[z_index]*1
1605 * @param hist_bin_content to be summed on
1606 * @param hitlist list of hits
1607 * @param layer layer.
1608 * @param axis 0:X, 1:V
1609 */
1610void InductionGar::conv1PerGrid_legacy(int Save_Grid, double * hist_bin_Content,
1611 const std::vector<int> & hitlist ,int layer, int axis
1612 ){
1613 //double ***** Signal=0;
1614 // this should not work, since you get no dimention info
1615 // and this following one might work.
1616
1617 double(&Signal)[3][5][5][4][100] = axis ? SignalV : SignalX;
1618 // our old C++ compiler cling(from old ROOT) complains about this.
1619 // thankfully, our g++ does not.
1620
1621 const int _nbins = 2000;
1622 int z_grid_x = Save_Grid / 100;
1623 int z_grid_v = Save_Grid % 100 / 10;
1624 int z_index = Save_Grid % 10;
1625 for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
1626 if (*it_start_bin <0 ) continue;
1627 for (int addi = *it_start_bin; (addi < *it_start_bin + 160) and (addi < _nbins - 1); addi += 2) { // 160 cell-layer move. i believe more hits in a semi-cell make fft-convolve relatively faster.
1628 hist_bin_Content[addi] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1]; // caution, its indice starts from...1?
1629 hist_bin_Content[addi + 1] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1];
1630 } // here its content is like stage with step size 2;
1631 }
1632}
1633
1634/**
1635 * @brief fft method of doing conv1 at a single grid for a strip. same usage as conv1PerGrid_legacy
1636 * discards t > _nbins or t <0 events
1637 * @param Save_Grid [grid_x]*100+[grid_v]*10+[z_index]*1
1638 * @param hist_bin_content to be summed on
1639 * @param hitlist list of hits
1640 * @param layer layer.
1641 * @param axis 0:X, 1:V
1642 */
1643void InductionGar::conv1PerGrid_fft(int Save_Grid, double * hist_bin_Content,
1644 const std::vector<int> & hitlist ,int layer, int axis
1645 ){
1646 ConvolveWithConst(&Signal_Convolver)[3][5][5][4] = axis ? Signal_ConvolverV : Signal_ConvolverX;
1647
1648 const int _nbins = 2000;
1649 double hist_tmp[_nbins] = {0};
1650 int z_grid_x = Save_Grid / 100;
1651 int z_grid_v = Save_Grid % 100 / 10;
1652 int z_index = Save_Grid % 10;
1653
1654 //bool found_mismatch = false;
1655 double electron_hit_tmp_hist[_nbins] = {0};
1656
1657 for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
1658 if(*it_start_bin<_nbins and *it_start_bin>=0) electron_hit_tmp_hist[*it_start_bin] += 1;
1659 }
1660 Signal_Convolver[layer][z_grid_x][z_grid_v][z_index].convolve(hist_tmp, electron_hit_tmp_hist, 0, _nbins);
1661
1662#ifdef INDUCTIONGAR_DEBUG_CONV1PERGRID_FFT_TEST
1663 {
1664 double hist_tmp2[_nbins] = {0};
1665 bool found_mismatch = false;
1666 const double Accept_Threshold = 1e-5;
1667 const double Accept_Threshold_ref = 1e-5;
1668 conv1PerGrid_legacy(Save_Grid, hist_tmp2, hitlist, layer, axis);
1669 for (int i = 0; i < _nbins; i++) {
1670 if (abs(hist_tmp[i] - hist_tmp2[i]) >= Accept_Threshold
1671 and abs((hist_tmp[i] - hist_tmp2[i]) / hist_tmp2[i]) >= Accept_Threshold_ref
1672 ) {
1673 if (not found_mismatch) {
1674 std::cout << "CgemDigitizerSvc::InductionGar::conv1PerGrid_fft found results not match: " << endl;
1675 found_mismatch = true;
1676 }
1677 std::cout << " previous: " << hist_tmp[i] << "; new " << hist_tmp2[i] << "; at i=" << i << std::endl;
1678 }
1679 }
1680 }
1681#endif
1682
1683 //add
1684 for (int i = 0; i < _nbins; i++) {
1685 hist_bin_Content[i] += hist_tmp[i];
1686 }
1687}
1688static void checkMemorySize()
1689{
1690 ProcInfo_t info;
1691 gSystem->GetProcInfo(&info);
1692 double currentMemory = info.fMemResident/1024;// MB
1693 cout<<"CgemDigitizerSvc::InductionGar: current memory size "<<currentMemory<<" MB"<<endl;
1694}
1695
1696bool InductionGar::inMicroSectorGap(double phi, int layer)
1697{
1698 double dphi=phi+M_PI;
1699 while(dphi<0) dphi+=2*M_PI;
1700 while(dphi>2*M_PI) dphi-=2*M_PI;
1701
1702 int i_sheet=0;
1703 if(layer>0) {
1704 //nGaps*=2;
1705 if(dphi>M_PI) i_sheet=1;
1706 }
1707 double width=m_microSector_width[layer][i_sheet];
1708 int iSector=floor((phi-m_gapShift_microSector[layer][i_sheet])/width);
1709 CgemGeoLayer* CgemLayer = m_geomSvc->getCgemLayer(layer);
1710 double R = CgemLayer->getInnerROfGapI();
1711 double dist = ((iSector+1)*width-(phi-m_gapShift_microSector[layer][i_sheet]))*R;
1712 bool in=false;
1713 if(dist<m_gap_microSector) in=true;
1714 return in;
1715}
1716
1717bool InductionGar::isDeadStrip(int layer, int sheet, int XVview, int strip)
1718{
1719 std::set<int> &deadStrips= m_deadStrip[layer][sheet][XVview];
1720 //cout<<"layer "<<layer<<" sheet "<<sheet<<" XVview "<<XVview<<" strip "<<strip;
1721 if (deadStrips.count(strip)==1){
1722 //cout<<"is dead strip"<<endl;
1723 return true;
1724 }
1725 //cout<<"is not dead"<<endl;
1726 return false;
1727}
Float_t
float Float_t
Definition CgemDigitizerSvc.h:42
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
I
const DifComplex I
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:52
Convolution_Tbranch
TH1F Convolution_Tbranch(const double *Input_Curr_plot_001, double xmin, double xmax, int Npx=2000)
Definition InductionGar.cxx:133
rectangle2
double rectangle2(double t, double *x, double *y, int Npx)
Definition InductionGar.cxx:44
Convolution_Ebranch
TH1F Convolution_Ebranch(const double *Input_Curr_plot_001, double xmin, double xmax, int Npx=2000)
Definition InductionGar.cxx:207
rectangle2_fast
double rectangle2_fast(double t, double *x, const double *y, int Npx)
calc estimate using linear interpolation of curve {x,y} , at point t; mimic rectangle2() but faster v...
Definition InductionGar.cxx:66
E_branch2
double E_branch2(double t)
Definition InductionGar.cxx:123
T_branch2
double T_branch2(double t)
Definition InductionGar.cxx:113
Vint
std::vector< int > Vint
Definition InductionGar.cxx:38
M_PI
#define M_PI
Definition TConstant.h:4
CEF::Convolution_Ebranch_fft
TH1D Convolution_Ebranch_fft(const double *Input_Curr_plot_001)
return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Ebranch ...
Definition InductionGar.cxx:440
CEF::EBranch
std::vector< double > EBranch
Definition InductionGar.h:31
CEF::conv
ConvolveWithConst * conv
Definition InductionGar.h:30
CTF::conv
ConvolveWithConst * conv
Definition InductionGar.h:18
CTF::TBranch
std::vector< double > TBranch
Definition InductionGar.h:19
CTF::Convolution_Tbranch_fft
TH1D Convolution_Tbranch_fft(const double *Input_Curr_plot_001)
return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Tbranch ...
Definition InductionGar.cxx:385
CgemGeoLayer::getInnerROfGapI
double getInnerROfGapI() const
Definition CgemGeoLayer.h:110
CgemGeoLayer::getNumberOfSheet
int getNumberOfSheet() const
Definition CgemGeoLayer.h:215
CgemGeoReadoutPlane::getDist2ClosestXStripCenter
double getDist2ClosestXStripCenter(double phi, int &id)
Definition CgemGeoReadoutPlane.cxx:143
CgemGeoReadoutPlane::OnThePlane
bool OnThePlane(double phi, double z) const
Definition CgemGeoReadoutPlane.h:154
CgemGeoReadoutPlane::getDist2ClosestVStripCenter
double getDist2ClosestVStripCenter(G4ThreeVector pos, int &id)
Definition CgemGeoReadoutPlane.cxx:217
ConvolveWithConst
basically, convolve with fft. convolves 1-d. if size of 2 inputs are fixed, it may yield optimal spee...
Definition ConvolveWithConst.h:18
ConvolveWithConst::init
void init(const double *ConstArray, const int ConstArrayLength, const int Blength, TransformOptimOpt opt=opt_EX, TransformLengthOpt optL=opt_AsShortAsPsb)
Definition ConvolveWithConst.cxx:60
ConvolveWithConst::convolve
void convolve(double *output, const double *B, const int leftIndex, const int sizeofOut, double factor=1) const
do a convolve of stored const A and B, and put results to output.
Definition ConvolveWithConst.cxx:166
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:500
InductionGar::setMultiElectrons
void setMultiElectrons(int layer, int nElectrons, const std::vector< Float_t > &x, const std::vector< Float_t > &y, const std::vector< Float_t > &z, const std::vector< Float_t > &t)
Definition InductionGar.cxx:767
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:27
t
int t()
Definition t.c:1