| 1 | #include <TROOT.h>
|
|---|
| 2 | #include <TCanvas.h>
|
|---|
| 3 | #include <TProfile.h>
|
|---|
| 4 | #include <TH1.h>
|
|---|
| 5 | #include <TH2.h>
|
|---|
| 6 | #include <TF1.h>
|
|---|
| 7 |
|
|---|
| 8 | #include <stdint.h>
|
|---|
| 9 | #include <cstdio>
|
|---|
| 10 |
|
|---|
| 11 | #define HAVE_ZLIB
|
|---|
| 12 | #include "fits.h"
|
|---|
| 13 | #include "FOpenDataFile.c"
|
|---|
| 14 | #include "FOpenCalibFile.c"
|
|---|
| 15 | #include "FPedestal.c"
|
|---|
| 16 | #include "FOscilloscope.c"
|
|---|
| 17 |
|
|---|
| 18 | #define n_peaks 4
|
|---|
| 19 |
|
|---|
| 20 | int peakfinderscope(const char *name, const char *drsname, size_t pixelnr)
|
|---|
| 21 | {
|
|---|
| 22 | //******************************************************************************
|
|---|
| 23 | //Read a datafile and plot the DRS-calibrated data
|
|---|
| 24 | //ATTENTION: only works for ROI=1024
|
|---|
| 25 | // (array indices of the calibration wrong otherwise)
|
|---|
| 26 | //******************************************************************************
|
|---|
| 27 |
|
|---|
| 28 | gROOT->SetStyle("Plain");
|
|---|
| 29 |
|
|---|
| 30 | //-------------------------------------------
|
|---|
| 31 | //Define the data variables
|
|---|
| 32 | //-------------------------------------------
|
|---|
| 33 | vector<int16_t> data;
|
|---|
| 34 | vector<int16_t> data_offset;
|
|---|
| 35 | unsigned int data_num;
|
|---|
| 36 | size_t data_n;
|
|---|
| 37 | UInt_t data_px;
|
|---|
| 38 | UInt_t data_roi;
|
|---|
| 39 |
|
|---|
| 40 | //-------------------------------------------
|
|---|
| 41 | //Get the DRS calibration
|
|---|
| 42 | //-------------------------------------------
|
|---|
| 43 | size_t drs_n;
|
|---|
| 44 | vector<float> drs_basemean;
|
|---|
| 45 | vector<float> drs_gainmean;
|
|---|
| 46 | vector<float> drs_triggeroffsetmean;
|
|---|
| 47 | FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
|
|---|
| 48 |
|
|---|
| 49 | //-------------------------------------------
|
|---|
| 50 | //Check the sizes of the data columns
|
|---|
| 51 | //-------------------------------------------
|
|---|
| 52 | // if(drs_n!=data_n)
|
|---|
| 53 | // {
|
|---|
| 54 | // cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
|
|---|
| 55 | // return 1;
|
|---|
| 56 | // }
|
|---|
| 57 |
|
|---|
| 58 | //-------------------------------------------
|
|---|
| 59 | //Create the canvas, plot and title
|
|---|
| 60 | //-------------------------------------------
|
|---|
| 61 | char title[500];
|
|---|
| 62 | std::sprintf(title,"Data: %s, DRS: %s, Px %i",name,drsname,pixelnr);
|
|---|
| 63 | TCanvas *canv = new TCanvas( "canv", "Waveform", 200, 10, 700, 500 );
|
|---|
| 64 | // TProfile *pulseshape = new TProfile("pulseshape", title, 70, -10.5, 59.5);
|
|---|
| 65 | TH2F *pulseshape = new TH2F("pulseshape",title,70,-10.5,59.5,350,-10.5,59.5);
|
|---|
| 66 | pulseshape->GetXaxis()->SetTitle("Time after trigger (bins @2 GHz)");
|
|---|
| 67 | pulseshape->GetYaxis()->SetTitle("Amplitude (mV)");
|
|---|
| 68 |
|
|---|
| 69 | char title_base[500];
|
|---|
| 70 | std::sprintf(title_base,"All samples of Px %i (data: %s, DRS: %s)",pixelnr,name,drsname);
|
|---|
| 71 |
|
|---|
| 72 | char title_spect[500];
|
|---|
| 73 | std::sprintf(title_spect,"Spectrum of all pixels (data: %s, DRS: %s)",name,drsname);
|
|---|
| 74 | TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
|
|---|
| 75 | TH1F *spectrum = new TH1F("spectrum", title_spect,700,-200,1200);
|
|---|
| 76 | spectrum->GetXaxis()->SetTitle("Pulse size (a.u.)");
|
|---|
| 77 | spectrum->GetYaxis()->SetTitle("Counts");
|
|---|
| 78 |
|
|---|
| 79 | //-------------------------------------------
|
|---|
| 80 | //Define the fit
|
|---|
| 81 | //-------------------------------------------
|
|---|
| 82 | int fit_start = 30;
|
|---|
| 83 | int fit_end = 380;
|
|---|
| 84 | TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp((-x-[11])/[12])",fit_start,fit_end);
|
|---|
| 85 | Double_t par[5+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,10000,30,30};
|
|---|
| 86 | // TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
|
|---|
| 87 | // Double_t par[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
|
|---|
| 88 | f_peaks->SetParameters(par);
|
|---|
| 89 | f_peaks->SetLineColor(2);
|
|---|
| 90 |
|
|---|
| 91 | TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
|
|---|
| 92 | f_pedestal->SetLineColor(2);
|
|---|
| 93 | f_pedestal->SetLineWidth(1);
|
|---|
| 94 |
|
|---|
| 95 | TF1* f_peak[n_peaks];
|
|---|
| 96 | for(int i=0; i<n_peaks; i++) {
|
|---|
| 97 | f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
|
|---|
| 98 | f_peak[i]->SetLineColor(2);
|
|---|
| 99 | f_peak[i]->SetLineWidth(1);
|
|---|
| 100 | }
|
|---|
| 101 |
|
|---|
| 102 | //-------------------------------------------
|
|---|
| 103 | //Open the file
|
|---|
| 104 | //-------------------------------------------
|
|---|
| 105 | fits datafile(name);
|
|---|
| 106 | if (!datafile)
|
|---|
| 107 | {
|
|---|
| 108 | cout << "Couldn't properly open the datafile." << endl;
|
|---|
| 109 | return 1;
|
|---|
| 110 | }
|
|---|
| 111 | FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
|
|---|
| 112 |
|
|---|
| 113 | //-------------------------------------------
|
|---|
| 114 | //Find the baseline
|
|---|
| 115 | //-------------------------------------------
|
|---|
| 116 | float pedestal_mean, pedestal_rms;
|
|---|
| 117 | // FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
|
|---|
| 118 | //Return values
|
|---|
| 119 | //Value of maximal probability: -2.225 +- 4.3144
|
|---|
| 120 | pedestal_mean = -2.225;
|
|---|
| 121 | pedestal_rms = 4.3144;
|
|---|
| 122 |
|
|---|
| 123 | //-------------------------------------------
|
|---|
| 124 | //Oscilloscope
|
|---|
| 125 | //-------------------------------------------
|
|---|
| 126 | // float threshold = pedestal_mean+pedestal_rms;
|
|---|
| 127 | float threshold = 2.5;
|
|---|
| 128 |
|
|---|
| 129 | //***********************
|
|---|
| 130 | FOscilloscope(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
|
|---|
| 131 | canv->cd();
|
|---|
| 132 | pulseshape->Draw();
|
|---|
| 133 | canv->Modified();
|
|---|
| 134 | canv->Update();
|
|---|
| 135 |
|
|---|
| 136 | canv_spect->cd();
|
|---|
| 137 | spectrum->Draw();
|
|---|
| 138 | canv_spect->Modified();
|
|---|
| 139 | canv_spect->Update();
|
|---|
| 140 | //***********************
|
|---|
| 141 | spectrum->Fit(f_peaks,"R+");
|
|---|
| 142 | f_peaks->GetParameters(&par[0]);
|
|---|
| 143 |
|
|---|
| 144 | f_peak[0]->SetParameters(par[2],par[0],par[3]);
|
|---|
| 145 | f_peak[1]->SetParameters(par[4],par[0]+par[1],par[5]);
|
|---|
| 146 | f_peak[2]->SetParameters(par[6],par[0]+2*par[1],par[7]);
|
|---|
| 147 | f_peak[3]->SetParameters(par[8],par[0]+3*par[1],par[9]);
|
|---|
| 148 | // f_peak[4]->SetParameters(par[10],par[12],par[11]);
|
|---|
| 149 | f_pedestal->SetParameters(par[10],par[11],par[12]);
|
|---|
| 150 |
|
|---|
| 151 | std::cout << "Crosstalk propability approx. " << ((par[4]*par[5])/(par[2]*par[3]+par[4]*par[5])) << std::endl;
|
|---|
| 152 | //-------------------------------------------
|
|---|
| 153 | //Draw the data
|
|---|
| 154 | //-------------------------------------------
|
|---|
| 155 | canv->cd();
|
|---|
| 156 | pulseshape->Draw();
|
|---|
| 157 | canv->Modified();
|
|---|
| 158 | canv->Update();
|
|---|
| 159 |
|
|---|
| 160 | canv_spect->cd();
|
|---|
| 161 | canv_spect->SetLogy();
|
|---|
| 162 | spectrum->Draw();
|
|---|
| 163 | f_peaks->Draw("same");
|
|---|
| 164 |
|
|---|
| 165 | for(int i=0; i<n_peaks; i++) {
|
|---|
| 166 | f_peak[i]->Draw("same");
|
|---|
| 167 | }
|
|---|
| 168 |
|
|---|
| 169 | f_pedestal->Draw("same");
|
|---|
| 170 |
|
|---|
| 171 | canv_spect->Modified();
|
|---|
| 172 | canv_spect->Update();
|
|---|
| 173 |
|
|---|
| 174 | return 0;
|
|---|
| 175 | }
|
|---|
