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) : KVSignal(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 
162 
163 KVChargeSignal::~KVChargeSignal()
164 {
165  // Destructor
166 }
167 
168 
169 
170 
178 
179 void KVChargeSignal::Copy(TObject& obj) const
180 {
181  // This method copies the current state of 'this' object into 'obj'
182  // You should add here any member variables, for example:
183  // (supposing a member variable KVChargeSignal::fToto)
184  // CastedObj.fToto = fToto;
185  // or
186  // CastedObj.SetToto( GetToto() );
187 
188  KVSignal::Copy(obj);
189  //KVChargeSignal& CastedObj = (KVChargeSignal&)obj;
190 }
191 
192 
193 // //________________________________________________________________
194 // KVPSAResult* KVChargeSignal::TreateSignal()
195 // {
196 // //to be implemented in child class
197 // if (GetN()==0) {
198 // //Info("TreateSignal","Empty signal %s",GetName());
199 // return 0;
200 // }
201 // KVPSAResult* psa = new KVPSAResult(GetName());
202 // psa->SetValue(Form("%s.%s.RawAmplitude",fDetName.Data(),fType.Data()),GetRawAmplitude());
203 // if (fAdc.fN==0) SetADCData();
204 //
205 // ComputeBaseLine();
206 // Add(-1.*fBaseLine);
207 // ApplyModifications();
208 //
209 // if(fWithPoleZeroCorrection) PoleZeroSuppression(fTauRC);
210 // FIR_ApplyTrapezoidal(fTrapRiseTime,fTrapFlatTop);
211 //
212 // ComputeAmplitude();
213 //
214 // //Init(); remplacer par SetADCData()
215 // //
216 // SetADCData();
217 // ComputeRiseTime();
218 //
219 // // storing result
220 // psa->SetValue(Form("%s.%s.BaseLine",fDetName.Data(),fType.Data()),fBaseLine);
221 // psa->SetValue(Form("%s.%s.SigmaBaseLine",fDetName.Data(),fType.Data()),fSigmaBase);
222 // psa->SetValue(Form("%s.%s.Amplitude",fDetName.Data(),fType.Data()),fAmplitude);
223 // psa->SetValue(Form("%s.%s.RiseTime",fDetName.Data(),fType.Data()),fRiseTime);
224 //
225 // return psa;
226 //
227 //
228 // /*
229 // Int_t xmin_bl = 0;
230 // Int_t xmax_bl = 200;
231 // psa->SetValue("BaseLine.Xmin",xmin_bl);
232 // psa->SetValue("BaseLine.Xmax",xmax_bl);
233 // Double_t xx,yy;
234 // GetPoint(0,xx,yy);
235 // */
236 // Int_t nn = GetN();
237 // if (!fFunc1){
238 // //Info("TreateSignal","Creation de la fonction de fit");
239 // fFunc1 = new TF1("KVChargeSignal1","[0]+[1]/(1+TMath::Exp([2]-x))",0,nn-1); fFunc1->SetNpx(nn);
240 // fFunc2 = new TF1("KVChargeSignal2","[0]+([1]+[3]*x)/(1+TMath::Exp([2]-x))",0,nn-2); fFunc2->SetNpx(nn);
241 // }
242 //
243 // Double_t par0 = GetMean(2)-2*GetRMS(2); //BaseLine
244 // Double_t par1 = GetRMS(2); //Amplitude
245 // Double_t par2 = GetMean(1); //Front
246 //
247 // fFunc1->SetParameters(par0,par1,par2);
248 // fFunc1->SetParLimits(1,0,1e5);
249 // fFunc1->SetParLimits(2,50,nn);
250 //
251 // //Determination ligne de base par fit pol0
252 // //
253 // Int_t times=0;
254 // Int_t nres1 = Int_t(Fit(fFunc1,"0WRQ"));
255 // if (nres1!=0 && times<3){
256 // nres1 = Int_t(Fit(fFunc1,"0WRQ"));
257 // times+=1;
258 // }
259 // if (nres1!=0){
260 // Error("TreateSignal","%s:%s : Fit #1 of the signal failed",fDetName.Data(),fType.Data());
261 // return 0;
262 // }
263 // for (Int_t nn=0;nn<fFunc1->GetNpar();nn+=1)
264 // fFunc2->SetParameter(nn,fFunc1->GetParameter(nn));
265 // fFunc2->SetParameter(fFunc2->GetNpar()-1,0.);
266 //
267 // Int_t nres2 = Int_t(Fit(fFunc2,"0WRQ"));
268 // if (nres2!=0){
269 // Warning("TreateSignal","%s:%s : Fit #2 of the signal failed",fDetName.Data(),fType.Data());
270 // psa->SetValue("BaseLine.Fit",fFunc1->GetParameter(0));
271 // psa->SetValue("Amplitude.Fit",fFunc1->GetParameter(1));
272 // psa->SetValue("Front0.Fit",fFunc1->GetParameter(2));
273 // psa->SetValue("Front1.Fit",fFunc1->GetParameter(3));
274 // }
275 // else{
276 // psa->SetValue("BaseLine.Fit",fFunc2->GetParameter(0));
277 // psa->SetValue("Amplitude.Fit",fFunc2->GetParameter(1));
278 // psa->SetValue("Front0.Fit",fFunc2->GetParameter(2));
279 // psa->SetValue("Front1.Fit",fFunc2->GetParameter(3));
280 // }
281 // //Min et Max de la ligne de base
282 // /*
283 // Double_t blmin = 1e6;
284 // Double_t blmax = -1e6;
285 // for (Int_t nn=xmin_bl;nn<xmax_bl;nn+=1)
286 // {
287 // GetPoint(nn,xx,yy);
288 // yy-=ybl;
289 // if (yy>blmax) blmax=yy;
290 // if (yy<blmin) blmin=yy;
291 // }
292 //
293 // Double_t trigger = 0;
294 // Double_t threshold = blmax;
295 //
296 // psa->SetValue("Noise.Threshold",blmax);
297 //
298 // Double_t max_signal=-1e6;
299 // Double_t tmax_signal=-1;
300 //
301 // for (Int_t nn=0;nn<GetN();nn+=1)
302 // {
303 // GetPoint(nn,xx,yy);
304 // yy-=ybl;
305 // if (yy>threshold){
306 // trigger += 1;
307 // if (yy>max_signal){
308 // max_signal = yy;
309 // tmax_signal = xx;
310 // }
311 // }
312 // SetPoint(nn,xx,yy);
313 // }
314 // psa->SetValue("Noise.Trigger",trigger/GetN());
315 // psa->SetValue("Signal.Max",max_signal);
316 // psa->SetValue("Signal.tMax",tmax_signal);
317 // */
318 // delete fFunc1;
319 // delete fFunc2;
320 // return psa;
321 //
322 // }
323 // //________________________________________________________________
324 // KVPSAResult* KVChargeSignal::TreateSignal(TF1* filter)
325 // {
326 // //to be implemented in child class
327 // KVPSAResult* psa = new KVPSAResult(GetName());
328 // Int_t nn = GetN();
329 // if (!filter){
330 // return 0;
331 // }
332 //
333 // Double_t xx,yy,y0;
334 // Double_t ymin=1e6,ymax=-1e6;
335 // GetPoint(0,xx,y0);
336 // for (Int_t nn=0;nn<GetN();nn+=1){
337 // GetPoint(nn,xx,yy);
338 // if (yy>ymax) ymax=yy;
339 // if (yy<ymin) ymin=yy;
340 // }
341 //
342 // Double_t width = ymax-ymin;
343 //
344 // Double_t par0 = y0; //GetMean(2)-3*GetRMS(2); //BaseLine
345 // Double_t par1 = width; //2*GetRMS(2); //Amplitude
346 // Double_t par2 = GetMean(1); //Front
347 //
348 // filter->SetParameters(par0,par1,par2,40);
349 // //filter->SetParameters(-7350,par1,590,70);
350 // filter->SetParLimits(1,0,width*1.1);
351 // filter->SetParLimits(2,50,nn);
352 // filter->SetParLimits(3,0.01,nn);
353 // //filter->SetParLimits(Parameter(3,1);
354 // //filter->FixParameter(4,1);
355 //
356 // /*
357 // filter->SetParLimits(3,-1,0.0001);
358 // filter->SetParLimits(4,0.000,100);
359 // */
360 // //Determination ligne de base par fit pol0
361 // //
362 // Int_t times=0;
363 // Int_t nres = Int_t(Fit(filter,"WRQ"));
364 // while (nres!=0 && times<10){
365 // nres = Int_t(Fit(filter,"WRQ"));
366 // times+=1;
367 // //printf("on recommence Fit #1 %d/10\n",times);
368 // }
369 // if (nres!=0){
370 // Error("TreateSignal","%s:%s : Fit #1 of the signal failed",fDetName.Data(),fType.Data());
371 // return 0;
372 // }
373 //
374 //
375 // /*
376 // filter->ReleaseParameter(3);
377 // filter->SetParLimits(3,-1,0.000);
378 // */
379 // /*
380 // filter->ReleaseParameter(4);
381 // filter->SetParLimits(4, 0,10);
382 // */
383 // /*
384 // times=0;
385 // nres = Int_t(Fit(filter,"WRQ"));
386 // while (nres!=0 && times<10){
387 // nres = Int_t(Fit(filter,"WRQ"));
388 // times+=1;
389 // printf("on recommence Fit #2 %d/10\n",times);
390 // }
391 // if (nres!=0){
392 // Error("TreateSignal","%s : Fit #2 of the signal failed",fDetName.Data(),fType.Data());
393 // return 0;
394 // }
395 // */
396 // //Int_t np=0;
397 // psa->SetValue("width",width);
398 //
399 // for (Int_t nn=0;nn<filter->GetNpar();nn+=1)
400 // filter->ReleaseParameter(nn);
401 //
402 // return psa;
403 //
404 // }
405 
406 
408 
409 Double_t KVChargeSignal::GetMaxFluctuationsWindow(Double_t* window, Int_t width)
410 {
411  if (!window) return 0;
412  TH1F* h2 = new TH1F("windows", "windows", (TMath::Nint(GetAmplitude()) + 1), GetYmin() - 0.5, GetYmax() + 0.5);
413  TGraph* gmoy = new TGraph();
414  TGraph* grms = new TGraph();
415  if (bidim) delete bidim;
416  bidim = new TGraph();
417  Double_t xx, yy;
418  Int_t nsamples = GetN() / width;
419 
420  printf("npoints=%d - nsamples=%d - width=%d\n", GetN(), nsamples, width);
421 
422  Double_t rms_max = 0;
423  Int_t irms_max = -1;
424  for (Int_t ii = 0; ii < nsamples; ii += 1) {
425  h2->Reset();
426  Double_t xinf = 0;
427  for (Int_t nn = 0; nn < width; nn += 1) {
428  GetPoint(width * ii + nn, xx, yy);
429  if (nn == 0) xinf = xx;
430  h2->Fill(yy);
431  }
432  Double_t xmoy = xinf;
433  Double_t rms = h2->GetRMS();
434  Double_t mean = h2->GetMean();
435  printf("%d %lf %lf %lf\n", ii, xmoy, rms, mean);
436  gmoy->SetPoint(ii, xmoy, mean);
437  grms->SetPoint(ii, xmoy, rms);
438  bidim->SetPoint(ii, mean, rms);
439  if (rms > rms_max) {
440  //endroit ou les fluctuations sont maximums
441  rms_max = h2->GetRMS();
442  irms_max = ii;
443  }
444  }
445 
446  Int_t ideb = irms_max;
447  Double_t ybefore = rms_max;
448  Double_t yafter = -1;
449  Int_t iparcours = ideb - 1;
450 
451  grms->GetPoint(iparcours, xx, yafter);
452  while (yafter < ybefore) {
453  ybefore = yafter;
454  iparcours -= 1;
455  grms->GetPoint(iparcours, xx, yafter);
456  }
457  Double_t xgauche = 0;
458  grms->GetPoint(iparcours + 1, xgauche, yy);
459  //xgauche-=width/2.;
460  Double_t ybas = 0;
461  gmoy->GetPoint(iparcours + 1, xx, ybas);
462 
463  ideb = irms_max;
464  ybefore = rms_max;
465  iparcours = ideb + 1;
466 
467  grms->GetPoint(iparcours, xx, yafter);
468  while (yafter < ybefore) {
469  ybefore = yafter;
470  iparcours += 1;
471  grms->GetPoint(iparcours, xx, yafter);
472  }
473  Double_t xdroite = 0;
474  grms->GetPoint(iparcours - 1, xdroite, yy);
475  //xdroite+=width/2.;
476  Double_t yhaut = 0;
477  gmoy->GetPoint(iparcours - 1, xx, yhaut);
478 
479  window[0] = xgauche;
480  window[1] = xdroite;
481  window[2] = ybas;
482  window[3] = yhaut;
483 
484  delete h2;
485  delete gmoy;
486  delete grms;
487  return rms_max;
488 }
489 
490 
KVChargeSignal::KVChargeSignal
KVChargeSignal()
Default constructor.
Definition: KVChargeSignal.cpp:37
KVChargeSignal::~KVChargeSignal
virtual ~KVChargeSignal()
Destructor.
Definition: KVChargeSignal.cpp:163
KVChargeSignal::GetMaxFluctuationsWindow
Double_t GetMaxFluctuationsWindow(Double_t *window, Int_t width=10)
Definition: KVChargeSignal.cpp:409
KVChargeSignal::LoadPSAParameters
virtual void LoadPSAParameters()
Definition: KVChargeSignal.cpp:116
KVChargeSignal::SetDefaultValues
virtual void SetDefaultValues()
Definition: KVChargeSignal.cpp:60
KVChargeSignal::Copy
void Copy(TObject &obj) const
Definition: KVChargeSignal.cpp:179
KVSignal::SetTauRC
void SetTauRC(Int_t taurc)
Definition: KVSignal.h:300
KVSignal::fType
TString fType
string to identify the signal type : "QH1", "I2" etc ...
Definition: KVSignal.h:32
KVSignal::Copy
void Copy(TObject &obj) const
Definition: KVSignal.cpp:169
KVSignal::GetAmplitude
Double_t GetAmplitude() const
Definition: KVSignal.h:250
KVSignal::SetPoleZeroCorrection
void SetPoleZeroCorrection(Bool_t with=kTRUE)
Definition: KVSignal.h:296
KVSignal::SetBaseLineLength
void SetBaseLineLength(Int_t length, Int_t first=0)
Definition: KVSignal.h:196
KVSignal::SetTrapShaperParameters
void SetTrapShaperParameters(Double_t rise, Double_t flat)
Definition: KVSignal.h:258
KVSignal::SetType
void SetType(const Char_t *type)
Definition: KVSignal.h:94
KVSignal::SetChannelWidth
void SetChannelWidth(double width)
Definition: KVSignal.h:167
KVSignal::GetYmax
Double_t GetYmax() const
Definition: KVSignal.h:337
KVSignal::GetYmin
Double_t GetYmin() const
Definition: KVSignal.h:333