WFAnalysis.cpp 21.2 KB
Newer Older
1
2
/** @ingroup ana
 *  @file WFAnalysis.cpp
3
4
 *  @brief Implementation of WFAnalysis.
 *
5
 *  Function definitions for WFAnalysis are provided.
6
7
 *  This class is the main  class for the waveform analysis.
 *  It initializes the histograms used for output of processed events.
8
 *  Also includes methods that accept all waveforms in an event.
9
10
 *  This is where the analysis should be done.
 *
11
 *  @author Sheng Yang, Chad Lantz, Yakov Kulinich
12
13
14
15
16
17
 *  @bug No known bugs.
 */

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
18
#include <TRandom.h>
19
#include <TTree.h>
20
#include <TCanvas.h>
21
22
23
24
25
26
#include <iostream>

#include "WFAnalysis.h"

/** @brief Default Constructor for WFAnalysis.
 */
Yakov Kulinich's avatar
Yakov Kulinich committed
27
WFAnalysis::WFAnalysis( ){
28
29
30

}

31

32
33
34
35
/** @brief Destructor for WFAnalysis.
 */
WFAnalysis::~WFAnalysis( ){

36
37
    delete f;
    delete nlFit;
38

39
40
}

41

42
43
/** @brief Initialization method for WFAnalysis
 *
44
 *  Can add other things here that you would
45
46
47
48
49
 *  perhaps not put into the constructor.
 *  I.e. a TTree, some tools. Etc.
 *
 */
void WFAnalysis::Initialize( ){
50

51
     nlFit = new TF1("nl1", "(x<660)*(-0.188891+1.03623*x) + \(x>660 && x<850)*( -51210.2 + 260.781*x - 0.4916*pow(x,2) + 0.00041227*pow(x,3) - (1.29681e-7)*pow(x,4) ) + \(x>850)*( -526857 + 1844.37*x - 2.1493*pow(x,2) + 0.000834971*pow(x,3) )", 0, 1000);
52
53
54
55
56
57
58
59
60
61
62

}


/** @brief Historgam Setup method for WFAnalysis
 *
 *  Should instantiate any histograms you wish to output here.
 *
 *  @return none
 */
void WFAnalysis::SetupHistograms( ){
63

clantz's avatar
clantz committed
64
  f = new TF1("f","gaus",-50,50);
65
  hDiffHisto = new TH1D("diffhist","First Derivative values histogrammed",75,-75,75);
66
  hPed = new TH1D("ped","ped", 1024, -512, 512);
67

68
69
70
71
}


/** @brief Analyze Events method for WFAnalysis
72
 *  @param vWF 1D vector of all waveforms for N ch
73
74
75
76
77
78
79
80
81
82
83
84
 *
 *  Here is the event-based analysis code.
 *  A const 1D vector of N TH1* is received.
 *  N is the number of channels
 *
 */
void WFAnalysis::AnalyzeEvent( const std::vector< TH1* >& vWFH ){

  // example... you can loop through all histos as follows
  for( unsigned int ch = 0; ch < vWFH.size(); ch++ ){
    TH1* h = vWFH[ch];
    for( unsigned int samp = 0; samp < h->GetNbinsX(); samp++ ){
85
      //Unused
86
87
88
89
90
91
    }
  }
}


/** @brief Analyze Events method for WFAnalysis
92
 *  @param vWF const 2D vector for all waveforms in event
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
 *
 *  Here is the event-based analysis code.
 *  A const 2D vector of NxM is received.
 *  N is the number of channels
 *  M is the number of samples per channel
 *
 *  @return none
 */
void WFAnalysis::AnalyzeEvent( const std::vector< std::vector< float > >& vWF ){

  // example... you can loop through each sample in each channel
  // Can do the same in previous AnalyzeEvent( .. ) function, just
  // looking at vWFH[ch]->GetBinContent( samp + 1 );
  for( unsigned int ch = 0; ch < vWF.size(); ch++ ){
    for( unsigned int samp = 0; samp < vWF[ch].size(); samp++ ){
108
      //Unused
109
110
111
112
    }
  }
}

113

114
115
/**
 * @brief Analyze Event method for WF analysis
116
117
 * @param vCh A vector of pointers to Channel objects
 *
118
 *  Receives a const vector of Channels. Loops over all Channels within
119
120
121
 *  the vector for analysis. Inverts channels with negative polarity before processing \n  \n
 *
 *  Processing occurs in the following order: \n
122
123
 *  Median filtering - If selected by the user, apply median filtering to the waveform \n
 *  FFT low pass filtering - If selected by the user, apply FFT low pass filter to the waveform \n
124
125
126
127
128
129
 *  Hit detection - Finds a hit within the waveform and gives the time if the hit was present \n
 *  Pedestal finding - Histograms the waveform outside of the hit window to determine pedestal and pedestal rms \n
 *  Pedestal subtraction - Removes the pedestal from the processed waveform \n
 *  Zero suppression - Sets all bins outside of the hit window and below a threshold to zero \n
 *  Integration - Integrates the processed waveform within the hit window accoutinging for bin timing calibration. Provides a value in pC \n
 *  Peak information - Uses the information gathered from the previous functions to save information about the waveform peak and derivative peak to the given channel
130
131
132
 *
 */
void WFAnalysis::AnalyzeEvent( const std::vector< Channel* > vCh ){
133
  for( Channel* Ch : vCh ){
134
    if( !Ch->is_on ) continue;
135
136
137
138
139
140
141
142
143
144
145
146
147
148

    //Clear old hits
    Ch->was_hit = false;
    Ch->nHits = 0;
    Ch->hit_window.first.clear();
    Ch->hit_window.second.clear();
    Ch->Peak_center.clear();
    Ch->Diff_Peak_center.clear();
    Ch->Peak_time.clear();
    Ch->Diff_Peak_time.clear();
    Ch->Peak_max.clear();
    Ch->Diff_max.clear();
    Ch->Charge.clear();

149
150
151
152
153
154
155
156
157
158
159
160
161
162
    //retrieving information for each channel as a histogram
    TH1D* h = Ch->WF_histo;
    TH1D* hProcessed = Ch->PWF_histo;
    TH1D* hDiff = Ch->FirstDerivative;

    //Invert the signal if it has negative polarity
    if( !Ch->pos_polarity ){
      for(int bin = 0; bin < h->GetNbinsX(); bin++){
        h->SetBinContent( bin+1, -1*Ch->pWF->at(bin) );
        Ch->pWF->at(bin) = -1*Ch->pWF->at(bin) ;
      }
      h->ResetStats();
      h->GetYaxis()->SetRange(0, 0);
    }
163

164
165
    //Apply filters as requested by the user
    if( Ch->doMedianFilter ) MedianFilter( Ch );
166
    if( Ch->doLPF ) LowPassFilter( Ch, h );
167

168
169
170
    //If the hit window has been fixed by the user, set was_hit to true and set the
    //hit window to the user provided values
    //Otherwise, find the hit window
171
172
    if( m_FixHitWindow ){
      Ch->was_hit = true;
173
174
175
      Ch->nHits = 1;
      Ch->hit_window.first.push_back( m_fixed_hit_window.first );
      Ch->hit_window.second.push_back( m_fixed_hit_window.second );
176
177
178
179
180
181
182
183
    }else{
      //Get the first derivative and determine the hit window from it
      //if a valid hit window is not found was_hit is set to false
      GetDifferential( Ch );

      //We set the offset of the DRS4 channels to -250 at first. This lead to clipping
      //of the baseline. If that's the case, estimate the RMS to be 5 based on good runs.
      Ch->FirstDerivativeRMS = ( Ch->offset != -250 ) ? GetRMS( Ch ) : 5.0;
184
      FindHitWindow( Ch );
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    }//end if else fixed hit window


    //Get and subtract the pedestal from the data outside the hit window
    GetPedestal( Ch );

    //If the channel was hit, proceed with processing
    if( Ch->was_hit ){
      SubtractPedestal( Ch );

      //The DRS4 saturates around 800-900 mV. If we have those values, flag it as saturated
      //std::cout << "MAX = " << Ch->PWF_histo->GetMaximum() << std::endl;
      //THIS IS ONLY FOR 2018 TEST BEAM - VERY RARE
      //if( Ch->PWF_histo->GetMaximum() > 890.0 ){
      //  Ch->saturated = true;
      //  std::cout << "====SATURATED=====" << std::endl;
      //}
202
      //Calibrate out the DRS4 non-linearity if requested by the user
203
204
205
206
      if( Ch->DRS4NLcomp ) DRS4Cal( Ch );
      //Zero Suppress the processed waveform vector and retrieve the energy related values from it
      ZeroSuppress( Ch );

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
      for(int hit = 0; hit < Ch->nHits; hit++){
        Ch->Charge.push_back( GetCharge( Ch, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );

        //Get peak information from the waveform using the hit window
        Ch->Peak_center.push_back( GetMaximumBin( hProcessed, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );
        Ch->Peak_max.push_back( hProcessed->GetBinContent( Ch->Peak_center.back() ) );
        Ch->Peak_time.push_back( Ch->pTimeVec->at( Ch->Peak_center.back() ) );

        //Set the range of the fit function to just around the peak and feed it some starting parameters
        //Then perform the fit
        // int range = 0.05*(Ch->hit_window.second[hit] - Ch->hit_window.first[hit]);
        // Ch->FitFunc->SetRange( Ch->Peak_center.back() - range, Ch->Peak_center.back() + range);
        // Ch->FitFunc->SetParameters(Ch->Peak_max, 3, 0.5, Ch->hit_window.first[hit]);
        // hProcessed->Fit(Ch->FitFunc,"QNR");//Fit (Q)uietly, do (N)ot draw and use the specified (R)ange
      }//end loop over hits
222
    }//end if channel->was_hit
223
224
225
226
227
228
    hProcessed->ResetStats();
    hProcessed->GetYaxis()->SetRangeUser(0, 0);
    hDiff->ResetStats();
    hDiff->GetYaxis()->SetRange(0, 0);
  }//end channel loop
}//end AnalyzeEvent
229

230

231
232
/** @brief MedianFilter method for WFAnalysis
 *  @param Ch channel to perform Median filtering on
233
 *
234
235
236
 *  Uses the median of the raw waveform in a sliding
 *  window to determine the value of the bin at the
 *  center of the window
237
238
 *
 */
239
240
241
242
243
244
245
246
247
248
void WFAnalysis::MedianFilter( Channel* Ch ){
  int nBins = Ch->WF.size();
  std::vector< float > output( nBins, 0 );
  std::vector< float > window;

  //Fill the front and back of the output vector since we can't filter
  //those bins
  for(int bin = 0; bin < Ch->medFiltHalfWindow; bin++){
    output[bin] = Ch->WF[bin];
    output[nBins -1 - bin] = Ch->WF[nBins -1 - bin];
249
  }
250
251
252
253
254
255
  //Go through each bin except for the first and last halfWindow bins
  //Set the filtered bin value to the median of values of the bins in
  //the window around the current bin
  for(int bin = Ch->medFiltHalfWindow; bin < nBins - Ch->medFiltHalfWindow; bin++){
    for(int i = 0; i < 2*Ch->medFiltHalfWindow; i++){
      window.push_back( Ch->WF[ bin - Ch->medFiltHalfWindow + i ] );
256
    }
257
258
259
260
261
262
    output[bin] = TMath::Median(2*Ch->medFiltHalfWindow, &window[0]);
    window.clear();
  }
  //Replace the raw waveform with the filtered one
  for(int bin = 0; bin < nBins; bin++){
    Ch->WF[bin] = output[bin];
263
264
265
  }
}

266

267
/** @brief GetDifferential method for WFAnalysis
268
 *  @param Ch channel who's waveform will be differentiated
269
 *
270
 * The derivative is calculated as the difference of the summation N points before and after the ith point
271
 *
272
273
274
 *  @f[
 *       \delta_i(N)= \sum_{k=1}^{N} (s_i + k)  -  \sum_{k=1}^{N} (s_i - k)
 *  @f]
275
 * N is given by Ch::Threshold and is set by Detector::SetWFAthresholds()
276
 */
277
void WFAnalysis::GetDifferential( Channel* Ch ){
278

279
280


281
  //Pad the front of the vector with 0s
282
  for(int bin = 0; bin < Ch->diffSmoothing; bin ++) Ch->vDiff.push_back( 0.0 );
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297


  //decalare temp variable to store the sum before and after the i data point
  double sum_before = 0;
  double sum_after = 0;

  // Generate 2 sliding window, sum_before and sum_after
  // Calculate the sum before and after the start data point
  for (int i = 1; i <= Ch->diffSmoothing; i++  ){
    sum_before += Ch->WF.at( Ch->diffSmoothing - i );
    sum_after  += Ch->WF.at( Ch->diffSmoothing + i );
  }

  // Loop over histogram
  for( unsigned int bin = Ch->diffSmoothing; bin < Ch->WF_histo->GetNbinsX() - Ch->diffSmoothing - 1 ; bin++ ){
298
        //set the bin to the calculated derivative value
299
        Ch->vDiff.push_back( sum_after - sum_before);
300
        Ch->FirstDerivative->SetBinContent(bin,(sum_after - sum_before));
301

302
303
304
305
306
307
308
309
        // Move the sum_before window forward
        sum_before += Ch->WF.at( bin );
        sum_before -= Ch->WF.at( bin - Ch->diffSmoothing );

        // Move the sum_after window forward
        sum_after -= Ch->WF.at( bin + 1);
        sum_after += Ch->WF.at( bin + Ch->diffSmoothing + 1);

310
    }//end derivative loop
311
312

    //Pad the back of the vector with 0s
313
    for(int bin = 0; bin <= Ch->diffSmoothing; bin ++) Ch->vDiff.push_back( 0.0 );
314
315
}

316
/** @brief GetRMS method for WFAnalysis
317
318
 *  @param Input Channel
 *  @return Width of gaussian fit
319
 *
320
321
 *  Histgrams the second derivative for a given Channel and gets the baseline RMS
 *  from a gaussian fit which ignores tails created by the peaks.
322
323
 *
 */
324
325
double WFAnalysis::GetRMS( Channel* Ch ){

326
    Double_t xmin,xmax;
327
    Ch->FirstDerivative->GetMinimumAndMaximum(xmin,xmax);
328
    if(xmax == 0) return 0;
329

330
331
    hDiffHisto->Reset();
    hDiffHisto->ResetStats();
332
    hDiffHisto->SetBins( (xmax - xmin)/12, xmin, xmax);
333

334
    //Loop over the histogram excluding the window used for differentiating to fill hRMS
335
    Int_t nbins = Ch->FirstDerivative->GetNbinsX();
336
    for(int bin = Ch->diffSmoothing; bin < nbins - Ch->diffSmoothing; bin++){
337
        hDiffHisto->Fill( Ch->vDiff.at( bin ) );
338
    }
339

340
    //Set the range of the gaussian fit to something reasonable, then perform the fit
341
    f->SetRange( hDiffHisto->GetMean() - 1.5*hDiffHisto->GetRMS(), hDiffHisto->GetMean() + 1.5*hDiffHisto->GetRMS() );
342
    f->SetParameters( hDiffHisto->GetMaximum() , hDiffHisto->GetMean(), hDiffHisto->GetRMS()/2.0);
343
    hDiffHisto->Fit("f","qR");
344

345
    //Return parameter 2. "gaus" is [0]*exp(-0.5*((x-[1])/[2])**2)
clantz's avatar
clantz committed
346
347
    //And reset the histogram
    return f->GetParameter(2);
348

349
350
}

351

352
/** @brief Defines the hit window for a given channel
353
354
355
356
357
358
359
360
361
362
  * @param ch Channel to be processed
  *
  * This function first checks if the FirstDerivative has any bins above threshold. If not, there will be no hit recorded.
  * In this case, we set datum to 0.0 in the case of float/doubles, and -1 in the case of ints.
  * If the derivative has bins above threshold the function steps through the waveform first derivative bin by bin looking for the first bin above threshold.
  * Once it finds the first bin above threshold it looks for the following landmarks in order:
  * FirstDerivative maximum, FirstDerivative zero crossing, FirstDerivative minimum, FirstDerivative zero crossing.
  * Which give us the corresponding datum:
  * Hit window start, Peak center, N/A, Hit window end.
  */
363
void WFAnalysis::FindHitWindow( Channel* ch ){
364

365
  double threshold = ch->Threshold*ch->FirstDerivativeRMS;
366
367

  //If there are no bins above threshold, we won't find a hit
368
  if( ch->FirstDerivative->GetMaximum() < threshold ||  ch->FirstDerivative->GetMinimum() > -1*threshold ) return;
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391

  hPed->Reset();
  int min = ch->WF_histo->GetMinimum();
  int max = ch->WF_histo->GetMaximum();
  hPed->SetBins( (max - min)/2, min, max); //Set the range of the histogram tailored to each waveform


  int nBins = ch->WF.size();
  for(int bin = 0; bin < nBins; bin++ ){
    hPed->Fill( ch->WF[bin] );
  }

  double mean = hPed->GetMean();
  double rms = hPed->GetRMS();

  hPed->Reset();


  float currentVal, prevVal;
  int startBin=0, maxBin=0, zeroBin=0, minBin=0, endBin=0;
  //Search for the first bin above threshold
  for(int bin = 1; bin < nBins; bin++){
    startLooking:
392
    startBin = maxBin = zeroBin = minBin = endBin = 0;
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
    if( ch->vDiff.at(bin) < threshold )continue;

    startBin = bin;

    //From the start of the hit window, look for the maximum of the derivative
    prevVal = ch->vDiff.at(bin);
    for( bin; bin < nBins; bin++){
      currentVal = ch->vDiff.at(bin);

      //If we have passed the maximum, the previous value was the max
      if( currentVal < prevVal ){
        maxBin = bin - 1;
        max = prevVal;
        break;
      }else{
        prevVal = currentVal;
      }//end if else currentVal < prevVal
    }//end looking for maxBin

    //From the maxBin, start looking for the zero crossing
    for(bin; bin < nBins; bin++){
      // if( ch->FirstDerivative->GetBinContent(bin) < 0.0 && ch->WF.at(bin) > 1.5*rms ){
      if( ch->vDiff.at(bin) < 0.0 ){
        zeroBin = bin - 1;
        break;
      }
    }//end looking for zero crossing

    //From the zero bin, start looking for the minimum bin
    //If the minimum bin is below threshold (technically above) go back to the beginning
    prevVal = ch->vDiff.at( zeroBin );
    for(bin; bin < nBins; bin++){
      currentVal = ch->vDiff.at( bin );
      if( currentVal > prevVal ){

428
        // if( prevVal > -1*threshold )continue;
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
        if( prevVal > -1*threshold )goto startLooking;

        minBin = bin - 1;
        min = prevVal;
        break;
      }else{
        prevVal = currentVal;
      }//end if currentVal > prevVal
    }//end looking for minimum

    if(minBin == 0){
      startBin = maxBin = zeroBin = 0;
      continue;
    }

444
    //Find the next time the derivative goes to zero,
445
446
447
448
449
450
451
452
    //this is the end of the hit
    for(bin; bin < nBins; bin++){
      if( ch->vDiff.at(bin) > 0.0 ){
        endBin = bin - 1;
        break;
      }
    }

453
454
455
456
    //If any of the values remain at default, start the loop over
    //This also exits the loop if we made it to the last bin
    if( startBin == 0 || maxBin == 0 || zeroBin == 0 || minBin == 0 || endBin == 0 ) continue;

457
458
459
460
    //If we've gotten here, we should have a hit
    //For now just assign values and break
    //TODO: Add multiple hit storage
    ch->was_hit = true;
461
462
463
464
465
466
    ch->nHits++;
    ch->hit_window.first.push_back( startBin );
    ch->hit_window.second.push_back( endBin );
    ch->Diff_Peak_center.push_back( maxBin );
    ch->Diff_Peak_time.push_back( ch->pTimeVec->at(maxBin) );
    ch->Diff_max.push_back( ch->vDiff.at(maxBin) );
467
468
469
470
471
472

  }//end loop over all bins

}//end FindHitWindowNew


473
474
475
476
477
478
/** @brief Finds the pedestal and pedestal rms of the raw waveform
 *  @param ch Channel to be processed
 *
 *  Histograms data from the raw waveform excluding the hit window.
 *  Saves the mean and rms of that histogram as PedMean and PedRMS.
 */
479
480
void WFAnalysis::GetPedestal( Channel* ch ){
    int nBins = ch->WF_histo->GetNbinsX();
481
482
483
    int min = ch->WF_histo->GetMinimum();
    int max = ch->WF_histo->GetMaximum();
    hPed->SetBins( (max - min)/2, min, max); //Set the range of the histogram tailored to each waveform
484

485
486
487
488
489
    //If the channel was hit, we only calculate the pedestal until the first hit
    //Otherwise use every bin
    int end = ( ch->was_hit ) ? ch->hit_window.first.front() : nBins;

    for(int bin = 0; bin < end; bin++){
490
          hPed->Fill( ch->pWF->at(bin) );
491
    }
492

493
494
    ch->PedMean = hPed->GetMean();
    ch->PedRMS  = hPed->GetRMS();
495

496
497
}

498
499
500
501
502


/** @brief Subtracts channel pedestal from the processed waveform histogram
*  @param ch Channel to be processed
*
503
*  Saves a pedestal subtracted version of the waveform to PWF_histo
504
505
506
507
508
509
510
511
512
513
514
*/
void WFAnalysis::SubtractPedestal( Channel* ch ){

  int nBins = ch->WF.size();
  for(int bin = 0; bin < nBins; bin++){
    ch->PWF_histo->SetBinContent( bin, ch->pWF->at(bin) - ch->PedMean );
  }

}


515
516
/** @brief Supress values below PedRMS
 *
517
 *  Changes all values outside of the hit window smaller than PedRMS after
518
519
520
521
522
523
 *  pedestal subtraction to zero.
 */
void WFAnalysis::ZeroSuppress( Channel* ch ){
    int nBins = ch->PWF_histo->GetNbinsX();
    for(int bin = 0; bin < nBins; bin++){
        double content = ch->PWF_histo->GetBinContent(bin);
524

525
526
527
528
529
530
        //Skip the hit windows
        for(int hit = 0; hit < ch->nHits; hit++){
          if( bin == ch->hit_window.first[hit] ){
            bin = ch->hit_window.second[hit] + 1;
          }
        }
531
532

        if( content <= 1.5*ch->PedRMS ){ ch->PWF_histo->SetBinContent( bin, 0 ); }
533
534
535
536
537
538
539
540
541
542
543
544
545
    }
}

/** @brief Applies a low pass filter to the processed waveform
 *  @param ch Channel to be processed
 *
 *  Applies FFT to the raw waveform, removes high frequencies
 *  and reconstructs the signal in PWF_histo.
 */
void WFAnalysis::LowPassFilter( Channel* ch, TH1D* hIn ){
   if(!hIn){ hIn = ch->WF_histo; }
   Int_t n = hIn->GetNbinsX();
   double re[n], im[n], reOut[n], imOut[n];
546

547
548
   TH1 *hm =0;
   hm = hIn->FFT(hm, "MAG");
549

550
   // Apply the filter
551
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
552
   fft->GetPointsComplex(re,im);
553
   for(int i = 0; i< ch->LPFfreq; i++){
554
555
556
       reOut[i] = re[i];
       imOut[i] = im[i];
   }
557
   for(int i = ch->LPFfreq; i < n; i++){
558
559
560
561
562
563
564
565
566
       reOut[i] = 0.0;
       imOut[i] = 0.0;
   }

   //Make a new TVirtualFFT with an inverse transform. Set the filtered points and apply the transform
   TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &n, "C2R M K");
   fft_back->SetPointsComplex(reOut,imOut);
   fft_back->Transform();
   ch->PWF_histo = (TH1D*)TH1::TransformHisto(fft_back,ch->PWF_histo,"MAG");
567

568
569
   //The transform scales the signal by sqrt(n) each time, so rescale by 1/n
   ch->PWF_histo->Scale(1.0/n);
570

571
572
573
   delete hm;
   delete fft;
   delete fft_back;
574

575
576
}

577
578
579
580
581
582

/** @brief Provides calibration for DRS4 voltage response
 *  @param ch Channel to be processed
 *
 *  Uses a TF1 with the DRS4 response to convert from the recorded value in mV
 *  to the actual value of the waveform in mV and reconstructs the signal in PWF_histo.
583
 *s
584
585
 */
void WFAnalysis::DRS4Cal( Channel* ch ){
586

587
588
589
590
591
592
593
594
    for(int bin = 0; bin < ch->PWF_histo->GetNbinsX(); bin++){
        ch->PWF_histo->SetBinContent(bin, nlFit->Eval( ch->PWF_histo->GetBinContent(bin) ) );
    }
}

/** @brief Uses time vector and input resistance to find the actual charge detected by the DRS4
 *  @param ch Channel to be processed
 *
595
596
 *  Integrates the processed waveform over the hit window. Uses calibrated timing information
 *  from the DRS4 to determine charge in pC
597
 */
598
599
600
601
602
603
604
605
606
double WFAnalysis::GetCharge( Channel* ch, int start, int end){
  double charge = 0.0;
  for(int bin = start; bin < end; bin++){
    //Charge for a given time bin = dt*I = (t(bin+1)-t(bin))*V/R
    //Units are Charge(pC), time(ns), Voltage(mV), resistance(ohm), current(Amp = C/s)
    charge += (ch->pTimeVec->at(bin+1) - ch->pTimeVec->at(bin)) * ch->PWF_histo->GetBinContent(bin)/Rin;
  }//end hit window
  return charge;
}//end GetCharge
607

608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631

/** @brief Get the maximum bin in a given range
 *  @return maxBin
 *  @param h - histogram under inspection
 *  @param _start - First bin of the desired range
 *  @param _end - Last bin in the desired range
 *
 *  Finds the maximum bin for the given histogram in the given range (in bins)
 *
 */
int WFAnalysis::GetMaximumBin( TH1* h, int _start, int _end ){
  float max = h->GetBinContent(_start);
  int maxBin = _start;

  for(int bin = _start+1; bin < _end; bin++){
    //If this bin isn't larger than the maximum identified, continue to the next
    if( h->GetBinContent(bin) < max) continue;

    max = h->GetBinContent(bin);
    maxBin = bin;
  }
  return maxBin;
}

632
/** @brief Finalize method for WFAnalysis
633
 *  @return none
634
 *
Yakov Kulinich's avatar
Yakov Kulinich committed
635
 *  Write histograms, TTree if it exists.
636
637
638
639
640
641
642
643
 *
 */
void WFAnalysis::Finalize( ){

  // If these exist...
  // m_tree->Write();
  // m_hist->Write();
}