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