KaliVeda
Toolkit for HIC analysis
KVChargeSignal.cpp
1 //Created by KVClassFactory on Tue Jan 13 15:10:39 2015
2 //Author: ,,,
3 
4 #include "KVChargeSignal.h"
5 #include "TMath.h"
6 #include "TH1F.h"
7 #include "TEnv.h"
8 #include "KVDataSet.h"
9 
10 ClassImp(KVChargeSignal)
11 
12 // BEGIN_HTML <!--
14 /* -->
15 <h2>KVChargeSignal</h2>
16 <h4>digitized charge signal</h4>
17 <!-- */
18 // --> END_HTML
20 
22 
23 void KVChargeSignal::init()
24 {
25 
26  fFunc1 = fFunc2 = 0;
27  bidim = 0;
28  SetDefaultValues();
29 
30 }
31 
32 
33 
36 
37 KVChargeSignal::KVChargeSignal()
38 {
39  // Default constructor
40  init();
41 }
42 
43 
44 
45 
48 
49 KVChargeSignal::KVChargeSignal(const char* name) : KVFAZIASignal(name, "Charge")
50 {
51  // Write your code here
52  SetType(name);
53  init();
54 }
55 
56 
57 
59 
60 void KVChargeSignal::SetDefaultValues()
61 {
62 
63  if (fType.Contains("QL")) {
64  SetChannelWidth(4.);
65  SetTauRC(250.);
66  }
67  else {
68  SetChannelWidth(10.);
69  SetTauRC(40.);
70  }
71  SetBaseLineLength(100);
72 
73  SetTrapShaperParameters(2., 1.);
74  SetPoleZeroCorrection();
75 }
76 
77 
78 
115 
116 void KVChargeSignal::LoadPSAParameters()
117 {
118 
119 // TString spar;
120 // Double_t lval;
121 // //--- BaseLineLength
122 // spar.Form("%s.%s.BaseLineLength",dettype,GetName());
123 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
124 // else lval = gEnv->GetValue(spar.Data(),0.0);
125 // SetBaseLineLength(lval);
126 // //--- ChannelWidth
127 // spar.Form("%s.%s.ChannelWidth",dettype,GetName());
128 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
129 // else lval = gEnv->GetValue(spar.Data(),0.0);
130 // SetChannelWidth(lval);
131 // //--- TauRC
132 // spar.Form("%s.%s.TauRC",dettype,GetName());
133 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
134 // else lval = gEnv->GetValue(spar.Data(),0.0);
135 // SetTauRC(lval);
136 //
137 // //--- ShaperRiseTime for trapezoidal filter
138 // spar.Form("%s.%s.ShaperRiseTime",dettype,GetName());
139 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
140 // else lval = gEnv->GetValue(spar.Data(),0.0);
141 // Double_t risetime = lval;
142 // //--- ShaperFlatTop for trapezoidal filter
143 // spar.Form("%s.%s.ShaperFlatTop",dettype,GetName());
144 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
145 // else lval = gEnv->GetValue(spar.Data(),0.0);
146 // Double_t flattop = lval;
147 // SetTrapShaperParameters(risetime,flattop);
148 // //--- Pole-Zero Correction
149 // spar.Form("%s.%s.PZCorrection",dettype,GetName());
150 // if (gDataSet) lval = gDataSet->GetDataSetEnv(spar.Data(),0.0);
151 // else lval = gEnv->GetValue(spar.Data(),0.0);
152 // if (lval==1) SetPoleZeroCorrection(kTRUE);
153 // else SetPoleZeroCorrection(kFALSE);
154 
155 
156 }
157 
158 
159 // //________________________________________________________________
160 // KVPSAResult* KVChargeSignal::TreateSignal()
161 // {
162 // //to be implemented in child class
163 // if (GetN()==0) {
164 // //Info("TreateSignal","Empty signal %s",GetName());
165 // return 0;
166 // }
167 // KVPSAResult* psa = new KVPSAResult(GetName());
168 // psa->SetValue(Form("%s.%s.RawAmplitude",fDetName.Data(),fType.Data()),GetRawAmplitude());
169 // if (fAdc.fN==0) SetADCData();
170 //
171 // ComputeBaseLine();
172 // Add(-1.*fBaseLine);
173 // ApplyModifications();
174 //
175 // if(fWithPoleZeroCorrection) PoleZeroSuppression(fTauRC);
176 // FIR_ApplyTrapezoidal(fTrapRiseTime,fTrapFlatTop);
177 //
178 // ComputeAmplitude();
179 //
180 // //Init(); remplacer par SetADCData()
181 // //
182 // SetADCData();
183 // ComputeRiseTime();
184 //
185 // // storing result
186 // psa->SetValue(Form("%s.%s.BaseLine",fDetName.Data(),fType.Data()),fBaseLine);
187 // psa->SetValue(Form("%s.%s.SigmaBaseLine",fDetName.Data(),fType.Data()),fSigmaBase);
188 // psa->SetValue(Form("%s.%s.Amplitude",fDetName.Data(),fType.Data()),fAmplitude);
189 // psa->SetValue(Form("%s.%s.RiseTime",fDetName.Data(),fType.Data()),fRiseTime);
190 //
191 // return psa;
192 //
193 //
194 // /*
195 // Int_t xmin_bl = 0;
196 // Int_t xmax_bl = 200;
197 // psa->SetValue("BaseLine.Xmin",xmin_bl);
198 // psa->SetValue("BaseLine.Xmax",xmax_bl);
199 // Double_t xx,yy;
200 // GetPoint(0,xx,yy);
201 // */
202 // Int_t nn = GetN();
203 // if (!fFunc1){
204 // //Info("TreateSignal","Creation de la fonction de fit");
205 // fFunc1 = new TF1("KVChargeSignal1","[0]+[1]/(1+TMath::Exp([2]-x))",0,nn-1); fFunc1->SetNpx(nn);
206 // fFunc2 = new TF1("KVChargeSignal2","[0]+([1]+[3]*x)/(1+TMath::Exp([2]-x))",0,nn-2); fFunc2->SetNpx(nn);
207 // }
208 //
209 // Double_t par0 = GetMean(2)-2*GetRMS(2); //BaseLine
210 // Double_t par1 = GetRMS(2); //Amplitude
211 // Double_t par2 = GetMean(1); //Front
212 //
213 // fFunc1->SetParameters(par0,par1,par2);
214 // fFunc1->SetParLimits(1,0,1e5);
215 // fFunc1->SetParLimits(2,50,nn);
216 //
217 // //Determination ligne de base par fit pol0
218 // //
219 // Int_t times=0;
220 // Int_t nres1 = Int_t(Fit(fFunc1,"0WRQ"));
221 // if (nres1!=0 && times<3){
222 // nres1 = Int_t(Fit(fFunc1,"0WRQ"));
223 // times+=1;
224 // }
225 // if (nres1!=0){
226 // Error("TreateSignal","%s:%s : Fit #1 of the signal failed",fDetName.Data(),fType.Data());
227 // return 0;
228 // }
229 // for (Int_t nn=0;nn<fFunc1->GetNpar();nn+=1)
230 // fFunc2->SetParameter(nn,fFunc1->GetParameter(nn));
231 // fFunc2->SetParameter(fFunc2->GetNpar()-1,0.);
232 //
233 // Int_t nres2 = Int_t(Fit(fFunc2,"0WRQ"));
234 // if (nres2!=0){
235 // Warning("TreateSignal","%s:%s : Fit #2 of the signal failed",fDetName.Data(),fType.Data());
236 // psa->SetValue("BaseLine.Fit",fFunc1->GetParameter(0));
237 // psa->SetValue("Amplitude.Fit",fFunc1->GetParameter(1));
238 // psa->SetValue("Front0.Fit",fFunc1->GetParameter(2));
239 // psa->SetValue("Front1.Fit",fFunc1->GetParameter(3));
240 // }
241 // else{
242 // psa->SetValue("BaseLine.Fit",fFunc2->GetParameter(0));
243 // psa->SetValue("Amplitude.Fit",fFunc2->GetParameter(1));
244 // psa->SetValue("Front0.Fit",fFunc2->GetParameter(2));
245 // psa->SetValue("Front1.Fit",fFunc2->GetParameter(3));
246 // }
247 // //Min et Max de la ligne de base
248 // /*
249 // Double_t blmin = 1e6;
250 // Double_t blmax = -1e6;
251 // for (Int_t nn=xmin_bl;nn<xmax_bl;nn+=1)
252 // {
253 // GetPoint(nn,xx,yy);
254 // yy-=ybl;
255 // if (yy>blmax) blmax=yy;
256 // if (yy<blmin) blmin=yy;
257 // }
258 //
259 // Double_t trigger = 0;
260 // Double_t threshold = blmax;
261 //
262 // psa->SetValue("Noise.Threshold",blmax);
263 //
264 // Double_t max_signal=-1e6;
265 // Double_t tmax_signal=-1;
266 //
267 // for (Int_t nn=0;nn<GetN();nn+=1)
268 // {
269 // GetPoint(nn,xx,yy);
270 // yy-=ybl;
271 // if (yy>threshold){
272 // trigger += 1;
273 // if (yy>max_signal){
274 // max_signal = yy;
275 // tmax_signal = xx;
276 // }
277 // }
278 // SetPoint(nn,xx,yy);
279 // }
280 // psa->SetValue("Noise.Trigger",trigger/GetN());
281 // psa->SetValue("Signal.Max",max_signal);
282 // psa->SetValue("Signal.tMax",tmax_signal);
283 // */
284 // delete fFunc1;
285 // delete fFunc2;
286 // return psa;
287 //
288 // }
289 // //________________________________________________________________
290 // KVPSAResult* KVChargeSignal::TreateSignal(TF1* filter)
291 // {
292 // //to be implemented in child class
293 // KVPSAResult* psa = new KVPSAResult(GetName());
294 // Int_t nn = GetN();
295 // if (!filter){
296 // return 0;
297 // }
298 //
299 // Double_t xx,yy,y0;
300 // Double_t ymin=1e6,ymax=-1e6;
301 // GetPoint(0,xx,y0);
302 // for (Int_t nn=0;nn<GetN();nn+=1){
303 // GetPoint(nn,xx,yy);
304 // if (yy>ymax) ymax=yy;
305 // if (yy<ymin) ymin=yy;
306 // }
307 //
308 // Double_t width = ymax-ymin;
309 //
310 // Double_t par0 = y0; //GetMean(2)-3*GetRMS(2); //BaseLine
311 // Double_t par1 = width; //2*GetRMS(2); //Amplitude
312 // Double_t par2 = GetMean(1); //Front
313 //
314 // filter->SetParameters(par0,par1,par2,40);
315 // //filter->SetParameters(-7350,par1,590,70);
316 // filter->SetParLimits(1,0,width*1.1);
317 // filter->SetParLimits(2,50,nn);
318 // filter->SetParLimits(3,0.01,nn);
319 // //filter->SetParLimits(Parameter(3,1);
320 // //filter->FixParameter(4,1);
321 //
322 // /*
323 // filter->SetParLimits(3,-1,0.0001);
324 // filter->SetParLimits(4,0.000,100);
325 // */
326 // //Determination ligne de base par fit pol0
327 // //
328 // Int_t times=0;
329 // Int_t nres = Int_t(Fit(filter,"WRQ"));
330 // while (nres!=0 && times<10){
331 // nres = Int_t(Fit(filter,"WRQ"));
332 // times+=1;
333 // //printf("on recommence Fit #1 %d/10\n",times);
334 // }
335 // if (nres!=0){
336 // Error("TreateSignal","%s:%s : Fit #1 of the signal failed",fDetName.Data(),fType.Data());
337 // return 0;
338 // }
339 //
340 //
341 // /*
342 // filter->ReleaseParameter(3);
343 // filter->SetParLimits(3,-1,0.000);
344 // */
345 // /*
346 // filter->ReleaseParameter(4);
347 // filter->SetParLimits(4, 0,10);
348 // */
349 // /*
350 // times=0;
351 // nres = Int_t(Fit(filter,"WRQ"));
352 // while (nres!=0 && times<10){
353 // nres = Int_t(Fit(filter,"WRQ"));
354 // times+=1;
355 // printf("on recommence Fit #2 %d/10\n",times);
356 // }
357 // if (nres!=0){
358 // Error("TreateSignal","%s : Fit #2 of the signal failed",fDetName.Data(),fType.Data());
359 // return 0;
360 // }
361 // */
362 // //Int_t np=0;
363 // psa->SetValue("width",width);
364 //
365 // for (Int_t nn=0;nn<filter->GetNpar();nn+=1)
366 // filter->ReleaseParameter(nn);
367 //
368 // return psa;
369 //
370 // }
371 
372 
374 
375 Double_t KVChargeSignal::GetMaxFluctuationsWindow(Double_t* window, Int_t width)
376 {
377  if (!window) return 0;
378  TH1F* h2 = new TH1F("windows", "windows", (TMath::Nint(GetAmplitude()) + 1), GetYmin() - 0.5, GetYmax() + 0.5);
379  TGraph* gmoy = new TGraph();
380  TGraph* grms = new TGraph();
381  if (bidim) delete bidim;
382  bidim = new TGraph();
383  Double_t xx, yy;
384  Int_t nsamples = GetN() / width;
385 
386  printf("npoints=%d - nsamples=%d - width=%d\n", GetN(), nsamples, width);
387 
388  Double_t rms_max = 0;
389  Int_t irms_max = -1;
390  for (Int_t ii = 0; ii < nsamples; ii += 1) {
391  h2->Reset();
392  Double_t xinf = 0;
393  for (Int_t nn = 0; nn < width; nn += 1) {
394  GetPoint(width * ii + nn, xx, yy);
395  if (nn == 0) xinf = xx;
396  h2->Fill(yy);
397  }
398  Double_t xmoy = xinf;
399  Double_t rms = h2->GetRMS();
400  Double_t mean = h2->GetMean();
401  printf("%d %lf %lf %lf\n", ii, xmoy, rms, mean);
402  gmoy->SetPoint(ii, xmoy, mean);
403  grms->SetPoint(ii, xmoy, rms);
404  bidim->SetPoint(ii, mean, rms);
405  if (rms > rms_max) {
406  //endroit ou les fluctuations sont maximums
407  rms_max = h2->GetRMS();
408  irms_max = ii;
409  }
410  }
411 
412  Int_t ideb = irms_max;
413  Double_t ybefore = rms_max;
414  Double_t yafter = -1;
415  Int_t iparcours = ideb - 1;
416 
417  grms->GetPoint(iparcours, xx, yafter);
418  while (yafter < ybefore) {
419  ybefore = yafter;
420  iparcours -= 1;
421  grms->GetPoint(iparcours, xx, yafter);
422  }
423  Double_t xgauche = 0;
424  grms->GetPoint(iparcours + 1, xgauche, yy);
425  //xgauche-=width/2.;
426  Double_t ybas = 0;
427  gmoy->GetPoint(iparcours + 1, xx, ybas);
428 
429  ideb = irms_max;
430  ybefore = rms_max;
431  iparcours = ideb + 1;
432 
433  grms->GetPoint(iparcours, xx, yafter);
434  while (yafter < ybefore) {
435  ybefore = yafter;
436  iparcours += 1;
437  grms->GetPoint(iparcours, xx, yafter);
438  }
439  Double_t xdroite = 0;
440  grms->GetPoint(iparcours - 1, xdroite, yy);
441  //xdroite+=width/2.;
442  Double_t yhaut = 0;
443  gmoy->GetPoint(iparcours - 1, xx, yhaut);
444 
445  window[0] = xgauche;
446  window[1] = xdroite;
447  window[2] = ybas;
448  window[3] = yhaut;
449 
450  delete h2;
451  delete gmoy;
452  delete grms;
453  return rms_max;
454 }
455 
456 
Int_t
int Int_t
Double_t
double Double_t
width
Option_t Option_t width
name
char name[80]
KVChargeSignal::KVChargeSignal
KVChargeSignal()
Default constructor.
Definition: KVChargeSignal.cpp:37
KVChargeSignal::GetMaxFluctuationsWindow
Double_t GetMaxFluctuationsWindow(Double_t *window, Int_t width=10)
Definition: KVChargeSignal.cpp:375
KVChargeSignal::LoadPSAParameters
void LoadPSAParameters() override
Definition: KVChargeSignal.cpp:116
KVChargeSignal::SetDefaultValues
void SetDefaultValues() override
Definition: KVChargeSignal.cpp:60
KVSignal::SetTauRC
void SetTauRC(Int_t taurc)
Definition: KVSignal.h:270
KVSignal::fType
TString fType
string to identify the signal type
Definition: KVSignal.h:41
KVSignal::GetAmplitude
Double_t GetAmplitude() const
Definition: KVSignal.h:220
KVSignal::SetPoleZeroCorrection
void SetPoleZeroCorrection(Bool_t with=kTRUE)
Definition: KVSignal.h:266
KVSignal::SetBaseLineLength
void SetBaseLineLength(Int_t length, Int_t first=0)
Definition: KVSignal.h:166
KVSignal::SetTrapShaperParameters
void SetTrapShaperParameters(Double_t rise, Double_t flat)
Definition: KVSignal.h:228
KVSignal::SetType
void SetType(const Char_t *type)
Definition: KVSignal.h:77
KVSignal::SetChannelWidth
void SetChannelWidth(double width)
Definition: KVSignal.h:137
KVSignal::GetYmax
Double_t GetYmax() const
Definition: KVSignal.h:307
KVSignal::GetYmin
Double_t GetYmin() const
Definition: KVSignal.h:303
TGraph::SetPoint
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
TGraph::GetN
Int_t GetN() const
TGraph::GetPoint
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
TH1F::Reset
void Reset(Option_t *option="") override
TH1::GetMean
virtual Double_t GetMean(Int_t axis=1) const
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::GetRMS
Double_t GetRMS(Int_t axis=1) const
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TMath::Nint
Int_t Nint(T x)
ClassImp
ClassImp(TPyArg)