CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
InductionGTS Class Reference

#include <InductionGTS.h>

+ Inheritance diagram for InductionGTS:

Public Member Functions

 InductionGTS ()
 
 ~InductionGTS ()
 
void init (ICgemGeomSvc *geomSvc, double magConfig)
 
void setDebugOutput (bool debugOutput)
 
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)
 
int getNXstrips () const
 
int getNVstrips () const
 
int getXstripSheet (int n) const
 
int getXstripID (int n) const
 
int getVstripSheet (int n) const
 
int getVstripID (int n) const
 
double getXstripQ (int n) const
 
double getVstripQ (int n) const
 
double getXstripT (int n) const
 
double getVstripT (int n) const
 
void clear ()
 
bool drive_to_anode (int, double, double, double, double, double &, double &, double &, double &)
 
bool useAPV (int stripid, int view, std::vector< double > stripidvec, std::vector< double > indchargevec, std::vector< double > indtimevec, double &charge, double &time, double &dtime)
 
void setVsampleDelay (double delay)
 
void setStoreFlag (bool flag)
 
void setLUTFilePath (std::string path)
 
void setSaturation (bool flag)
 
double getXstripT_Branch (int n) const
 
double getVstripT_Branch (int n) const
 
double getXstripQ_Branch (int n) const
 
double getVstripQ_Branch (int n) const
 
double getXfirstT (int n) const
 
double getVfirstT (int n) const
 
- Public Member Functions inherited from Induction
 Induction ()
 
virtual ~Induction ()
 

Detailed Description

Definition at line 23 of file InductionGTS.h.

Constructor & Destructor Documentation

◆ InductionGTS()

InductionGTS::InductionGTS ( )

Definition at line 31 of file InductionGTS.cxx.

31 {
32}

◆ ~InductionGTS()

InductionGTS::~InductionGTS ( )

Definition at line 34 of file InductionGTS.cxx.

35{
36
37 if(m_testing_ind == true) {
38 output->Write();
39 output->Close();
40 }
41
42 double tminbin = hcollected_charge->FindBin(apv_tstart);
43 // for(int it = tminbin; it < hnbin; it++) delete f[it];
44 delete hintegratore;
45 delete hcollected_charge;
46 delete hcharge;
47 delete f_FD;
48
49
50}
#define m_testing_ind
#define apv_tstart

Member Function Documentation

◆ clear()

void InductionGTS::clear ( )

Definition at line 125 of file InductionGTS.cxx.

125 {
126
127 // temp
128 stripid_x.clear();
129 sheetid_x.clear();
130 charge_x.clear();
131 time_x.clear();
132 stripid_v.clear();
133 sheetid_v.clear();
134 charge_v.clear();
135 time_v.clear();
136
137 m_NXstrips = 0;
138 m_NVstrips = 0;
139 m_XstripSheet.clear();
140 m_XstripID.clear();
141 m_VstripSheet.clear();
142 m_VstripID.clear();
143 m_XstripQ.clear();
144 m_VstripQ.clear();
145 m_XstripT.clear();
146 m_VstripT.clear();
147}

Referenced by setMultiElectrons().

◆ drive_to_anode()

bool InductionGTS::drive_to_anode ( int layer,
double xi,
double yi,
double zi,
double ti,
double & xf,
double & yf,
double & zf,
double & tf )

Definition at line 469 of file InductionGTS.cxx.

469 {
470
471 double R = TMath::Sqrt(xi*xi + yi*yi);
472 double phi = TMath::ATan2(yi, xi);
473
474 double shift_x_induct = diffusion->shift_x_induct();
475 double sigma_x_induct = diffusion->sigma_x_induct();
476
477 double gem3_rad_out = m_geomSvc->getCgemLayer(0)->getCgemFoil(2)->getOuterROfCgemFoil();
478 double shift_phi = shift_x_induct/gem3_rad_out;
479 double sigma_phi = sigma_x_induct/gem3_rad_out;
480 double phif = phi + gRandom->Gaus(shift_phi, sigma_phi);
481
482 double anode_rad_in = m_geomSvc->getCgemLayer(layer)->getInnerROfAnode();
483 xf = anode_rad_in * TMath::Cos(phif);
484 yf = anode_rad_in * TMath::Sin(phif);
485
486 // y
487 double shift_y_induct = diffusion->shift_y_induct();
488 double sigma_y_induct = diffusion->sigma_y_induct();
489 zf = zi + gRandom->Gaus(shift_y_induct, sigma_y_induct);
490
491 // cout << "xi " << xi << " yi " << yi << " zi " << zi << endl;
492 // cout << "R " << R << " phi " << phi << " phif " << phif << endl;
493 // cout << "xf " << xf << " yf " << yf << " zf " << zf << endl;
494
495 // t
496 double shift_t_induct = diffusion->shift_t_induct();
497 double sigma_t_induct = diffusion->sigma_t_induct();
498 tf = ti + gRandom->Gaus(shift_t_induct, sigma_t_induct);
499
500 double dphi = phif - phi;
501 double darc = anode_rad_in * dphi;
502
503
504 if(m_testing_ind) tree->Fill(dphi, darc);
505
506 return true;
507}
double getOuterROfCgemFoil() const
Definition CgemGeoFoil.h:34
CgemGeoFoil * getCgemFoil(int i) const
double getInnerROfAnode() const
double sigma_x_induct()
double shift_y_induct()
double shift_x_induct()
double shift_t_induct()
double sigma_y_induct()
double sigma_t_induct()
virtual CgemGeoLayer * getCgemLayer(int i) const =0
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:27

Referenced by setMultiElectrons().

◆ getNVstrips()

int InductionGTS::getNVstrips ( ) const
virtual

Implements Induction.

Definition at line 512 of file InductionGTS.cxx.

512{ return m_NVstrips; }

◆ getNXstrips()

int InductionGTS::getNXstrips ( ) const
virtual

Implements Induction.

Definition at line 511 of file InductionGTS.cxx.

511{ return m_NXstrips; }

◆ getVfirstT()

double InductionGTS::getVfirstT ( int n) const
inlinevirtual

Implements Induction.

Definition at line 59 of file InductionGTS.h.

59{return m_VstripT.at(n);}
const Int_t n

◆ getVstripID()

int InductionGTS::getVstripID ( int n) const
virtual

Implements Induction.

Definition at line 516 of file InductionGTS.cxx.

516{ return m_VstripID.at(n); }

◆ getVstripQ()

double InductionGTS::getVstripQ ( int n) const
virtual

Implements Induction.

Definition at line 518 of file InductionGTS.cxx.

518{ return m_VstripQ.at(n); }

◆ getVstripQ_Branch()

double InductionGTS::getVstripQ_Branch ( int n) const
inlinevirtual

Implements Induction.

Definition at line 57 of file InductionGTS.h.

57{return m_VstripQ.at(n);}

◆ getVstripSheet()

int InductionGTS::getVstripSheet ( int n) const
virtual

Implements Induction.

Definition at line 515 of file InductionGTS.cxx.

515{ return m_VstripSheet.at(n); }

◆ getVstripT()

double InductionGTS::getVstripT ( int n) const
virtual

Implements Induction.

Definition at line 520 of file InductionGTS.cxx.

520{ return m_VstripT.at(n); }

◆ getVstripT_Branch()

double InductionGTS::getVstripT_Branch ( int n) const
inlinevirtual

Implements Induction.

Definition at line 55 of file InductionGTS.h.

55{return m_VstripT.at(n);}

◆ getXfirstT()

double InductionGTS::getXfirstT ( int n) const
inlinevirtual

Implements Induction.

Definition at line 58 of file InductionGTS.h.

58{return m_XstripT.at(n);}

◆ getXstripID()

int InductionGTS::getXstripID ( int n) const
virtual

Implements Induction.

Definition at line 514 of file InductionGTS.cxx.

514{ return m_XstripID.at(n); }

◆ getXstripQ()

double InductionGTS::getXstripQ ( int n) const
virtual

Implements Induction.

Definition at line 517 of file InductionGTS.cxx.

517{ return m_XstripQ.at(n); }

◆ getXstripQ_Branch()

double InductionGTS::getXstripQ_Branch ( int n) const
inlinevirtual

Implements Induction.

Definition at line 56 of file InductionGTS.h.

56{return m_XstripQ.at(n);}

◆ getXstripSheet()

int InductionGTS::getXstripSheet ( int n) const
virtual

Implements Induction.

Definition at line 513 of file InductionGTS.cxx.

513{ return m_XstripSheet.at(n); }

◆ getXstripT()

double InductionGTS::getXstripT ( int n) const
virtual

Implements Induction.

Definition at line 519 of file InductionGTS.cxx.

519{ return m_XstripT.at(n); }

◆ getXstripT_Branch()

double InductionGTS::getXstripT_Branch ( int n) const
inlinevirtual

Implements Induction.

Definition at line 54 of file InductionGTS.h.

54{return m_XstripT.at(n);}

◆ init()

void InductionGTS::init ( ICgemGeomSvc * geomSvc,
double magConfig )
virtual

Implements Induction.

Definition at line 52 of file InductionGTS.cxx.

52 {
53 m_geomSvc = geomSvc;
54 m_magConfig = magConfig;
55 evt = 0;
56 // SETTINGS -------------------------------------------------
57 string filePath = getenv("CGEMDIGITIZERSVCROOT");
58 string fileName;
59
60 fileName = filePath + "/dat/par_settings_GTS.txt";
61 ifstream fin(fileName.c_str(), ios::in);
62
63 string strline;
64 string strpar;
65 if( ! fin.is_open() ){
66 std::cout << "InductionGTS::init ERROR: can not open file " << filePath << endl;
67 }
68 else std::cout << "InductionGTS::init: open file " << fileName << endl;
69
70
71 while(fin >> strpar){
72
73 if(strpar[0] == '#'){
74 getline(fin, strline);
75 }
76 else if(strpar == "magnetic_field") {
77 fin >> m_field;
78 }
79 else if(strpar == "tuning_factor_diff_perp") {
80 fin >> m_tuning_factor_diff_perp;
81 }
82 else if(strpar == "tuning_factor_diff_paral") {
83 fin >> m_tuning_factor_diff_paral;
84 }
85 }
86
87 cout << "InductionGTS::init: drift and avalanche parameters" << endl;
88 cout << "field " << m_field <<endl;
89 cout << "tuning factors for diffusion " << m_tuning_factor_diff_perp << " (for orthogonal) " << m_tuning_factor_diff_paral << " (for parallel)" << endl;
90
91 // diffusion orthogonal to mag field
92 if(m_field) fileName = filePath + "/dat/par_ArIso_1T_GTS.txt";
93 else fileName = filePath + "/dat/par_ArIso_0T_GTS.txt";
94
95
96 if(m_testing_ind == true) {
97 output = new TFile("induction.root", "RECREATE");
98 tree = new TNtuple("tree", "tree", "dphi:darc");
99 tree_strip = new TNtuple("tree_strip", "tree for strip", "evt:id:view:qind:q:t");
100 }
101
102 // diffusion parallel to mag field
103 string fileName2 = filePath + "/dat/par_ArIso_0T_GTS.txt";
104 hcollected_charge = new TH1F("hcollected_charge", "hcollected_charge", hnbin, hxmin, hxmax);
105 hintegratore = new TH1F("hintegratore", "", hnbin, hxmin, hxmax);
106 hcharge = new TH1F("hcharge", "ASIC measured charge", 28, hxmin, hxmax);
107 double tminbin = hcollected_charge->FindBin(apv_tstart);
108 double dt = hcollected_charge->GetXaxis()->GetBinWidth(1);
109
110 for(int it = tminbin; it < hnbin; it++) {
111 double t = apv_tstart + it * dt;
112 TString name = "f"; name += it;
113 f[it] = new TF1(name, "[0] * ((x - [1]) / [2]) * TMath::Exp(-(x - [1]) / [2])", t, hxmax);
114 }
115
116 f_FD = new TF1("f_FD", "[0] + [1]/(1+TMath::Exp(-(x - [2])/[3]))", hxmin, hxmax);
117
118
119
120 diffusion = new DiffusionGTS(m_geomSvc);
121 diffusion->readGasPerpParameters(fileName);
122 diffusion->readGasParalParameters(fileName2);
123}
TGraphErrors * dt
Definition AbsCor.cxx:72
#define hnbin
#define hxmax
#define hxmin
void readGasPerpParameters(std::string fileName)
void readGasParalParameters(std::string fileName)
int t()
Definition t.c:1

◆ setDebugOutput()

void InductionGTS::setDebugOutput ( bool debugOutput)
inlinevirtual

Implements Induction.

Definition at line 29 of file InductionGTS.h.

29{m_debugOutput = debugOutput;}

◆ setLUTFilePath()

void InductionGTS::setLUTFilePath ( std::string path)
inlinevirtual

Implements Induction.

Definition at line 52 of file InductionGTS.h.

52{return;}

◆ setMultiElectrons()

void InductionGTS::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 )
virtual

Implements Induction.

Definition at line 150 of file InductionGTS.cxx.

150 {
151 double time_spent = 0.;
152 clock_t begin = clock();
153
154 clear();
155
156 std::vector < double > stripid_x;
157 std::vector < double > ind_charge_x;
158 std::vector < double > ind_time_x;
159
160 std::vector < double > stripid_v;
161 std::vector < double > ind_charge_v;
162 std::vector < double > ind_time_v;
163
164 // cout << "NELETTRONI " << nElectrons << endl;
165 // ---------------------------------------------
166 // loop on ions in avalanche
167 for(int iele = 0; iele < nElectrons; iele++){
168
169 double x_start = x.at(iele);
170 double y_start = y.at(iele);
171 double z_start = z.at(iele);
172 double t_start = t.at(iele); // CHECK time here already contains the trigger time or not?? must be track_t0 + t_drift + t_diffusion
173
174 // diffusion is local
175 double x_end, y_end, z_end, t_end;
176 bool isanode = drive_to_anode(layer, x_start, y_start, z_start, t_start, x_end, y_end, z_end, t_end);
177
178 int sheet = 0; // CHECK THIS!!
179 double phi = TMath::ATan2(y_start, x_start);
180 if(layer > 0 && phi > TMath::Pi()) sheet = 1;
181
182 CgemGeoReadoutPlane* anode = m_geomSvc->getReadoutPlane(layer, sheet);
183 int X_ID;
184 int V_ID;
185 G4ThreeVector pos(x_end, y_end, z_end);
186
187 anode->getStripID(pos, X_ID, V_ID); // get the closest XID/VID
188
189 // cout << endl;
190 // cout << "start "<< x_start << " " << y_start << " " << z_start << " " << t_start << endl;
191 // cout << "end "<< x_end << " " << y_end << " " << z_end << " " << t_end << endl;
192 // cout << "X_ID " << X_ID << " V_ID " << V_ID << endl;
193
194 // cout << "->radius " << TMath::Sqrt(x_end*x_end + y_end*y_end)
195 // << " arco " << TMath::ATan2(y_end,x_end) * TMath::Sqrt(x_end*x_end + y_end*y_end)
196 // << " X_ID " << X_ID << endl;
197
198
199 // charge sharing // CHECK TO BE CHANGED
200 double percX = 0.5;
201 double percV = 0.5;
202 double qele = 1.6e-4;
203
204 // save all the vectors in parallel
205 stripid_x.push_back(X_ID);
206 ind_charge_x.push_back(percX * qele);
207 ind_time_x.push_back(t_end);
208 sheetid_x.push_back(sheet);
209 stripid_v.push_back(V_ID);
210 ind_charge_v.push_back(percV * qele);
211 ind_time_v.push_back(t_end);
212 sheetid_v.push_back(sheet);
213
214 }
215
216 // final output
217
218 // count and fill X stripIDs -----------------------------------------------------------
219 for (std::vector< double >::iterator it = stripid_x.begin() ; it != stripid_x.end(); it++) {
220 double stripid = *it;
221 int ipos = it - stripid_x.begin();
222 double sheetid = sheetid_x.at(ipos);
223
224 std::vector< double >::iterator it2 = std::find(m_XstripID.begin(), m_XstripID.end(), stripid);
225 std::vector< double >::iterator it3 = std::find(m_XstripSheet.begin(), m_XstripSheet.end(), sheetid);
226
227 if(it2 == m_XstripID.end() || it3 == m_XstripSheet.end()) { // if the strip is not present, then add stripid and sheetid
228 m_XstripID.push_back(stripid);
229 m_XstripSheet.push_back(sheetid);
230 }
231 else { // if the strip is present, then check if it belongs to the same sheetid
232 double sheetid3 = *it3;
233 if(sheetid != sheetid3) { // if the sheetid is not present, then add stripid and sheetid
234 m_XstripID.push_back(stripid);
235 m_XstripSheet.push_back(sheetid);
236 }
237 }
238 }
239 m_NXstrips = m_XstripID.size();
240 // cout << "nstrip X " << m_NXstrips << endl;
241
242 // measure charge, after electronics, now APV-25 ----------------------------------------
243 for(int istrip = 0; istrip < m_NXstrips; istrip++) {
244 double xID = m_XstripID.at(istrip);
245 int view = 0;
246 double strip_charge;
247 double strip_time;
248 double strip_dtime; // CHECK this!
249 bool isASIC = useAPV(xID, view, stripid_x, ind_charge_x, ind_time_x, strip_charge, strip_time, strip_dtime);
250 m_XstripQ.push_back(strip_charge);
251 m_XstripT.push_back(strip_time);
252 }
253
254
255 // count and fill V stripIDs -----------------------------------------------------------
256 for (std::vector< double >::iterator it = stripid_v.begin() ; it != stripid_v.end(); it++) {
257 double stripid = *it;
258 int ipos = it - stripid_v.begin();
259 double sheetid = sheetid_v.at(ipos);
260
261 std::vector< double >::iterator it2 = std::find(m_VstripID.begin(), m_VstripID.end(), stripid);
262 std::vector< double >::iterator it3 = std::find(m_VstripSheet.begin(), m_VstripSheet.end(), sheetid);
263
264 if(it2 == m_VstripID.end() || it3 == m_VstripSheet.end()) { // if the strip is not present, then add stripid and sheetid
265 m_VstripID.push_back(stripid);
266 m_VstripSheet.push_back(sheetid);
267 }
268 else { // if the strip is present, then check if it belongs to the same sheetid
269 double sheetid3 = *it3;
270 if(sheetid != sheetid3) { // if the sheetid is not present, then add stripid and sheetid
271 m_VstripID.push_back(stripid);
272 m_VstripSheet.push_back(sheetid);
273 }
274 }
275 }
276 m_NVstrips = m_VstripID.size();
277 // cout << "nstripv " << m_NVstrips << endl;
278
279
280 // count and fill V stripIDs -----------------------------------------------------------
281 for(int istrip = 0; istrip < m_NVstrips; istrip++) {
282 double vID = m_VstripID.at(istrip);
283 int view = 1;
284 // cout << "V stripid " << vID << " sheetid " << m_VstripSheet.at(istrip) << endl;
285 double strip_charge;
286 double strip_time;
287 double strip_dtime; // CHECK this!
288 bool isASIC = useAPV(vID, view, stripid_v, ind_charge_v, ind_time_v, strip_charge, strip_time, strip_dtime);
289 m_VstripQ.push_back(strip_charge);
290 m_VstripT.push_back(strip_time);
291 }
292
293 clock_t end = clock();
294 time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
295 cout << "InductionGTS::induzione " << time_spent << " seconds" << endl;
296
297 evt++;
298
299}
Double_t x[10]
void getStripID(G4ThreeVector pos, int &X_ID, int &V_ID) const
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
bool useAPV(int stripid, int view, std::vector< double > stripidvec, std::vector< double > indchargevec, std::vector< double > indtimevec, double &charge, double &time, double &dtime)
bool drive_to_anode(int, double, double, double, double, double &, double &, double &, double &)

◆ setSaturation()

void InductionGTS::setSaturation ( bool flag)
inlinevirtual

Implements Induction.

Definition at line 53 of file InductionGTS.h.

53{return;}

◆ setStoreFlag()

void InductionGTS::setStoreFlag ( bool flag)
inlinevirtual

Implements Induction.

Definition at line 51 of file InductionGTS.h.

51{return;}

◆ setVsampleDelay()

void InductionGTS::setVsampleDelay ( double delay)
inlinevirtual

Implements Induction.

Definition at line 50 of file InductionGTS.h.

50{return;}

◆ useAPV()

bool InductionGTS::useAPV ( int stripid,
int view,
std::vector< double > stripidvec,
std::vector< double > indchargevec,
std::vector< double > indtimevec,
double & charge,
double & time,
double & dtime )

Definition at line 303 of file InductionGTS.cxx.

304{
305
306 // cout << "--> x stripid " << stripid << " sheetid " << m_XstripSheet.at(istrip) << endl;
307
308 // cout << "APV for strip " << stripid << " no. " << istrip << "/" << m_NXstrips << endl;
309 // ----------- APV25 FAST SIMULATION ------------------
310
311 double tau_APV = 50; // ns
312 int ntime = 27; // time bin
313 double tstep = 25; // ns
314
315 hcollected_charge->Reset();
316 TString name;
317 if(m_testing_ind) {
318 name = "hcollected_charge";
319 name += stripid;
320 hcollected_charge->SetName(name);
321 }
322
323 std::vector< double >::iterator iter = stripidvec.begin();
324 while(1) {
325 if(iter != stripidvec.end()+1) iter = std::find(iter, stripidvec.end(), stripid);
326 else break;
327 int ipos = iter - stripidvec.begin();
328 std::vector< double >::iterator iter2 = indtimevec.begin() + ipos;
329 std::vector< double >::iterator iter3 = indchargevec.begin() + ipos;
330 hcollected_charge->Fill(*iter2, *iter3);
331 // cout << ipos << " time " << *iter2 << " charge " << *iter3 << endl;
332 iter++;
333 }
334 double ind_charge = hcollected_charge->Integral();
335
336 double dt = hcollected_charge->GetXaxis()->GetBinWidth(1);
337 // integrator
338 double tminbin = hcollected_charge->FindBin(apv_tstart);
339 for(int it = tminbin; it < hnbin; it++) {
340 double t = apv_tstart + it * dt;
341 double timebin = hcollected_charge->FindBin(t);
342 double integral = hcollected_charge->GetBinContent(timebin);
343
344 name = "f"; name += it; name += "_strip"; name += stripid;
345 f[it]->SetName(name);
346 // f[it] = new TF1(name, "[0] * ((x - [1]) / [2]) * TMath::Exp(-(x - [1]) / [2])", t, xmax);
347 f[it]->SetParameter(0, TMath::Exp(1.)*integral);
348 f[it]->SetParameter(1, t);
349 f[it]->SetParameter(2, tau_APV);
350 }
351
352 hintegratore->Reset();
353 if(m_testing_ind) {
354 name = "hintegratore";
355 name += stripid;
356 hintegratore->SetName(name);
357 }
358
359 for(int it = tminbin; it < hnbin; it++) {
360 double t = apv_tstart + it * dt;
361 double qmeas = 0;
362 for(int jt = tminbin; jt < hnbin; jt++) {
363 double t2 = apv_tstart + jt * dt;
364 if(t2 < t) qmeas += f[jt]->Eval(t);
365 }
366 hintegratore->SetBinContent(it, qmeas);
367 }
368
369 hcharge->Reset();
370 if(m_testing_ind) {
371 name = "hcharge";
372 name += stripid;
373 hcharge->SetName(name);
374 }
375
376 double jitter = 0.; // CHECK APV
377 for(int iapv = 0; iapv < ntime; iapv++) {
378 double t = jitter + apv_tstart + iapv * tstep;
379 double timebin = hintegratore->FindBin(t);
380 double qmeas = hintegratore->GetBinContent(timebin);
381 hcharge->AddBinContent(iapv+1, qmeas);
382 }
383
384 // output->Write();
385 // hcollected_charge->Write();
386 // hintegratore->Write();
387 // hcharge->Write();
388 // -----------------
389
390 // ------------FERMI DIRAC-------
391 // CHARGE
392 charge = hcharge->GetMaximum();
393 // cout << "charge " << charge << endl;
394 // TIME from FERMI DIRAC %%%%%%%%%%%%%%%%%%%%%%%%%%
395 // start parameters
396 double maxbin = hcharge->GetMaximumBin();
397 double QMaxHisto = hcharge->GetMaximum();
398
399 //find the lower edge of the rising @ 10% of the maximum
400 double minbin = 0;
401 for(int ibin = maxbin; ibin > 1; ibin--){
402 if(hcharge->GetBinContent(ibin) < 0.1 * QMaxHisto) {
403 minbin = ibin;
404 break;
405 }
406 }
407 double startvalues[5] = {0};
408 double fitlimup[5] = {0};
409 double fitlimlow[5] = {0};
410
411 // average of the three bins before minbin
412 double meanfirstbin = (hcharge->GetBinContent(minbin-1) + hcharge->GetBinContent(minbin-2) + hcharge->GetBinContent(minbin-3))/3.;
413 double sigma_MFB = TMath::Sqrt( ( pow(hcharge->GetBinContent(minbin-1) - meanfirstbin, 2) + pow(hcharge->GetBinContent(minbin-2) - meanfirstbin,2) + pow(hcharge->GetBinContent(minbin-3)- meanfirstbin,2)) / 3);
414 if(sigma_MFB < TMath::Abs(meanfirstbin)) sigma_MFB = TMath::Abs(meanfirstbin);
415
416 startvalues[0] = meanfirstbin;
417 startvalues[1] = QMaxHisto;
418 startvalues[2] = 0.5* (minbin + maxbin) * tstep;
419 startvalues[3] = 0.5*tstep;
420
421 // put limits
422 fitlimlow[0] = startvalues[0] - 2 * sigma_MFB;
423 fitlimlow[1] = startvalues[1] - 0.3 * startvalues[1];
424 fitlimlow[2] = minbin * tstep;
425 fitlimlow[3] = 0.1*tstep;
426
427 fitlimup[0] = startvalues[0] + 2 * sigma_MFB;
428 fitlimup[1] = startvalues[1] + 0.3 * startvalues[1];
429 fitlimup[2] = maxbin * tstep;
430 fitlimup[3] = 0.9*tstep;
431
432 // ..................
433 // the real Fermi Dirac
434 //
435 // FD(t) = K ----------------------------- + B
436 // 1 + exp( -(t - tFD)/sigFD )
437 // [0] = B
438 // [1] = K
439 // [2] = tFD
440 // [3] = sigFD
441 double minFD = minbin - 4;
442 double maxFD = maxbin + 0;
443 // cout<<"FERMI DIRAC FIT: \n";
444
445 f_FD->SetRange(minFD*tstep, maxFD*tstep);
446 name = "f_FD";
447 name += stripid;
448 f_FD->SetName(name);
449
450 f_FD->SetParameters(startvalues[0], startvalues[1], startvalues[2], startvalues[3]);
451 f_FD->SetParLimits(0, fitlimlow[0], fitlimup[0]);
452 f_FD->SetParLimits(1, fitlimlow[1], fitlimup[1]);
453 f_FD->SetParLimits(2, fitlimlow[2], fitlimup[2]);
454 f_FD->SetParLimits(3, fitlimlow[3], fitlimup[3]);
455
456 hcharge->Fit(f_FD,"WQRB0"); // quiet/range/params
457 time = 0;
458 time = gRandom->Gaus(f_FD->GetParameter(2), 8); //CHECK ??
459 // time = f_FD->GetParameter(2); // CHECK??
460 dtime = 0.;
461 dtime = f_FD->GetParError(2);
462
463 if(m_testing_ind) tree_strip->Fill(evt, stripid, view, ind_charge, charge, time);
464
465
466 //-------------
467}
#define tstep
Definition APV25GTS.cxx:12
Double_t time
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)

Referenced by setMultiElectrons().


The documentation for this class was generated from the following files: