Index: /fact/tools/rootmacros/FIntFixedPosAllPx.c =================================================================== --- /fact/tools/rootmacros/FIntFixedPosAllPx.c (revision 12166) +++ /fact/tools/rootmacros/FIntFixedPosAllPx.c (revision 12166) @@ -0,0 +1,54 @@ +#include +int FIntFixedPosAllPx(fits &datafile, vector &data, vector &data_offset, unsigned int &data_num, UInt_t data_px, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_roi, UInt_t integration_size, UInt_t integration_delay, TH1F* spectrum[]) +//Data, calibration, data_roi, data_num, data_px, threshold, two histograms +{ +// UInt_t integration_size = 5; +// UInt_t integration_delay = 238; + +// float sample, last_sample, integral; + float integral; + UInt_t drs_calib_offset; + +//********************************************************************************** + for (size_t i=0; (iFill(integral); +// for(Int_t l=-10; l<60; l++) +// { +// drs_calib_offset = (l+data_offset[j])%data_roi; +// sample = (data[j*data_roi+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35; +// } +// k+=min_dist; +// +// drs_calib_offset = (start_sample-1+data_offset[j])%data_roi; +// last_sample = (data[j*data_roi+start_sample-1]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+start_sample-1])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35; +// for (UInt_t k=start_sample; kthreshold)) +// { +// } +// //TBD: mistake that after the deadtime the last_sample must be set new... +// //Solution: first process full pipeline +// last_sample = sample; +// } + } + } + return 0; +} Index: /fact/tools/rootmacros/FOpenCalibFile.c =================================================================== --- /fact/tools/rootmacros/FOpenCalibFile.c (revision 12166) +++ /fact/tools/rootmacros/FOpenCalibFile.c (revision 12166) @@ -0,0 +1,45 @@ +#include +int FOpenCalibFile(const char *drsname, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, size_t &drs_n) +{ +//------------------------------------------- +//Open the file +//------------------------------------------- + fits drsfile(drsname); + + if (!drsfile) + { + cout << "Couldn't properly open the drsfile." << endl; + return 1; + } + +//------------------------------------------- +//Print the header +//------------------------------------------- + cout << "----------------- DRS Calib Header -----------------" << endl; + drsfile.PrintKeys(); + cout << "---------------- DRS Calib Columns -----------------" << endl; + drsfile.PrintColumns(); + +//------------------------------------------- +//Get the sizes of the data column +//------------------------------------------- + drs_n = drsfile.GetN("BaselineMean"); + +//------------------------------------------- +//Set the sizes of the DRS data vectors +//------------------------------------------- + drs_basemean.resize(drs_n,0); + drs_gainmean.resize(drs_n,0); + drs_triggeroffsetmean.resize(drs_n,0); + +//------------------------------------------- +//Link the data to variables +//------------------------------------------- + drsfile.SetVecAddress("BaselineMean", drs_basemean); + drsfile.SetVecAddress("GainMean", drs_gainmean); + drsfile.SetVecAddress("TriggerOffsetMean", drs_triggeroffsetmean); + drsfile.GetRow(0); //Read the calibration data + + cout << "Reading calibration file successful..." << endl; + return 0; +} Index: /fact/tools/rootmacros/FOpenDataFile.c =================================================================== --- /fact/tools/rootmacros/FOpenDataFile.c (revision 12166) +++ /fact/tools/rootmacros/FOpenDataFile.c (revision 12166) @@ -0,0 +1,35 @@ +#include +int FOpenDataFile(fits &datafile, vector &data, vector &data_offset, unsigned int &data_num, size_t &data_n, UInt_t &data_roi, UInt_t &data_px) +{ +//------------------------------------------- +//Print the header +//------------------------------------------- +/* cout << "-------------------- Data Header -------------------" << endl; + datafile.PrintKeys(); + cout << "------------------- Data Columns -------------------" << endl; + datafile.PrintColumns(); + */ +//------------------------------------------- +//Get the size of the data column +//------------------------------------------- + data_roi = datafile.GetUInt("NROI"); // Value from header + data_px = datafile.GetUInt("NPIX"); + data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI + +//------------------------------------------- +//Set the sizes of the data vectors +//------------------------------------------- + data.resize(data_n,0); + data_offset.resize(data_px,0); + +//------------------------------------------- +//Link the data to variables +//------------------------------------------- + datafile.SetRefAddress("EventNum", data_num); + datafile.SetVecAddress("Data", data); + datafile.SetVecAddress("StartCellData", data_offset); + datafile.GetRow(0); + + cout << "Opening data file successful..." << endl; + return 0; +} Index: /fact/tools/rootmacros/FOscilloscope.c =================================================================== --- /fact/tools/rootmacros/FOscilloscope.c (revision 12166) +++ /fact/tools/rootmacros/FOscilloscope.c (revision 12166) @@ -0,0 +1,59 @@ +#include +int FOscilloscope(fits &datafile, vector &data, vector &data_offset, unsigned int &data_num, UInt_t data_px, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_roi, float threshold, TH2* pulseshape, TH1* spectrum) +//Data, calibration, data_roi, data_num, data_px, threshold, two histograms +{ + UInt_t start_sample = 10; + UInt_t end_sample = data_roi-60; + UInt_t min_dist = 150; + UInt_t integration_size = 10; + UInt_t integration_delay = 5; + + float sample, last_sample, integral; + UInt_t drs_calib_offset; + +//********************************************************************************** + for (size_t i=0; ithreshold)) + { + integral = 0; + for(UInt_t l=integration_delay; l60)&&(integral<=300)) { +// if((integral>60)&&(integral<=300)) { + for(Int_t l=-10; l<60; l++) + { + drs_calib_offset = (k+l+data_offset[j])%data_roi; + sample = (data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35; + pulseshape->Fill(l,sample); + } +// } + spectrum->Fill(integral); + k+=min_dist; + } + //TBD: mistake that after the deadtime the last_sample must be set new... + //Solution: first process full pipeline + last_sample = sample; + } + } + } + return 0; +} Index: /fact/tools/rootmacros/FOscilloscopeAllPx.c =================================================================== --- /fact/tools/rootmacros/FOscilloscopeAllPx.c (revision 12166) +++ /fact/tools/rootmacros/FOscilloscopeAllPx.c (revision 12166) @@ -0,0 +1,59 @@ +#include +int FOscilloscopeAllPx(fits &datafile, vector &data, vector &data_offset, unsigned int &data_num, UInt_t data_px, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_roi, float threshold, TH1F* spectrum[]) +//Data, calibration, data_roi, data_num, data_px, threshold, two histograms +{ + UInt_t start_sample = 10; + UInt_t end_sample = data_roi-60; + UInt_t min_dist = 150; + UInt_t integration_size = 10; + UInt_t integration_delay = 5; + + float sample, last_sample, integral; + UInt_t drs_calib_offset; + +//********************************************************************************** + for (size_t i=0; ithreshold)) + { + integral = 0; + for(UInt_t l=integration_delay; l60)&&(integral<=300)) { +// if((integral>60)&&(integral<=300)) { + //Iterate for the pulseshape +// for(Int_t l=-10; l<60; l++) +// { +// drs_calib_offset = (k+l+data_offset[j])%data_roi; +// sample = (data[j*data_roi+k+l]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k+l])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35; +// } +// } + spectrum[j]->Fill(integral); + k+=min_dist; + } + //TBD: mistake that after the deadtime the last_sample must be set new... + //Solution: first process full pipeline + last_sample = sample; + } + } + } + return 0; +} Index: /fact/tools/rootmacros/FPedestal.c =================================================================== --- /fact/tools/rootmacros/FPedestal.c (revision 12166) +++ /fact/tools/rootmacros/FPedestal.c (revision 12166) @@ -0,0 +1,35 @@ +#include +int FPedestal(char *title_base, fits &datafile, vector &data, vector &data_offset, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_roi, size_t pixelnr, float &pedestal_mean, float &pedestal_rms) +{ +// TCanvas *canv_base = new TCanvas( "canv_base", "Determine the baseline", 100, 520, 700, 500 ); + TH1F *baseline = new TH1F("baseline", title_base, 1000,-10,40); +// canv_base->SetLogy(); + baseline->GetXaxis()->SetTitle("Sample value (mV)"); + baseline->GetYaxis()->SetTitle("Entries"); + +//------------------------------------------- +//Find the baseline +//------------------------------------------- + for (size_t i=0; iFill(sample); + } + } + pedestal_mean = baseline->GetXaxis()->GetBinCenter(baseline->GetMaximumBin()); + pedestal_rms = baseline->GetRMS(); + std::cout << "Value of maximal probability: " << pedestal_mean << " +- " << pedestal_rms << std::endl; + +// canv_base->cd(); +// baseline->Draw(); +// canv_base->Modified(); +// canv_base->Update(); + + return 0; +} Index: /fact/tools/rootmacros/FPedestalAllPx.c =================================================================== --- /fact/tools/rootmacros/FPedestalAllPx.c (revision 12166) +++ /fact/tools/rootmacros/FPedestalAllPx.c (revision 12166) @@ -0,0 +1,62 @@ +#include +int FPedestalAllPx(fits &datafile, vector &data, vector &data_offset, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_px, UInt_t data_roi, float pedestal_mean[], float pedestal_rms[]) +{ +// TCanvas *canv_base = new TCanvas( "canv_base", "Determine the baseline", 100, 520, 700, 500 ); +// canv_base->SetLogy(); + + char title_base[500]; + char name_base[50]; + TH1F* baseline[data_px]; + for(int i=0; iGetXaxis()->SetTitle("Sample value (mV)"); + baseline[i]->GetYaxis()->SetTitle("Entries"); + } + +// TH1F *baseline = new TH1F("baseline", title_base, 1000,-10,40); +// baseline[i]->GetXaxis()->SetTitle("Sample value (mV)"); +// baseline[i]->GetYaxis()->SetTitle("Entries"); + +//------------------------------------------- +//Find the baseline +//------------------------------------------- + for (size_t i=0; iFill(sample); + + //Original code + //UInt_t drs_calib_offset = (k+data_offset[j])%data_roi; + //float sample = (data[j*data_roi+k]*2000/4096.-drs_basemean[j*data_roi+drs_calib_offset]-drs_triggeroffsetmean[j*data_roi+k])/drs_gainmean[j*data_roi+drs_calib_offset]*1907.35; + //baseline->Fill(sample); + } + } + } + + for (int j=0; jGetXaxis()->GetBinCenter(baseline[j]->GetMaximumBin()); + pedestal_rms[j] = baseline[j]->GetRMS(); + } + std::cout << "Value of maximal probability: " << pedestal_mean[1348] << " +- " << pedestal_rms[1348] << std::endl; + +// canv_base->cd(); +// baseline->Draw(); +// canv_base->Modified(); +// canv_base->Update(); + + return 0; +} Index: /fact/tools/rootmacros/TPKplotevent.c =================================================================== --- /fact/tools/rootmacros/TPKplotevent.c (revision 12166) +++ /fact/tools/rootmacros/TPKplotevent.c (revision 12166) @@ -0,0 +1,20 @@ +#include + +int TPKplotevent(char *title, vector &data, vector &data_offset, vector &drs_basemean, vector &drs_gainmean, vector &drs_triggeroffsetmean, UInt_t data_roi, size_t pixelnr) +{ + TCanvas *canv = new TCanvas( "canv", "Mean values of the first event", 100, 10, 700, 500 ); + TProfile *pix = new TProfile("pix", title, 1024, -0.5, 1023.5); + + for (UInt_t k=0; kFill(k,sample); + } + pix->Draw(); + canv->Modified(); + canv->Update(); + + cout << "Plotting successful..." << endl; + return 0; +} Index: /fact/tools/rootmacros/calscope.C =================================================================== --- /fact/tools/rootmacros/calscope.C (revision 12166) +++ /fact/tools/rootmacros/calscope.C (revision 12166) @@ -0,0 +1,84 @@ +#include +#include +#include + +#include +#include + +#include "fits.h" +#include "TPKplotevent.c" +#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" + +int calscope(const char *name = "../raw/20110916_025.fits", const char *drsname = "../raw/20110916_024.drs.fits", size_t eventnr = 0, size_t pixelnr = 0) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//Example call in ROOT: +//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348) +//T. Kr�henb�hl, August 2011, tpk@phys.ethz.ch +//****************************************************************************** + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +//Open the file +//------------------------------------------- + fits datafile(name); + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + +//------------------------------------------- +//Get the data +//------------------------------------------- + vector data; + vector data_offset; + unsigned int data_num; + size_t data_n; + UInt_t data_px; + UInt_t data_roi; + FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px); + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + size_t drs_n; + vector drs_basemean; + vector drs_gainmean; + vector drs_triggeroffsetmean; + FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- + if(drs_n!=data_n) + { + cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; + return 1; + } + + +//------------------------------------------- +//Create the title +//------------------------------------------- + char title[500]; + std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr); + +//------------------------------------------- +//Get the event +//------------------------------------------- + cout << "--------------------- Data --------------------" << endl; + datafile.GetRow(eventnr); + cout << "Event number: " << data_num << endl; + +//------------------------------------------- +//Draw the data +//------------------------------------------- + TPKplotevent(title, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr); + return 0; +} Index: /fact/tools/rootmacros/calscope_persistent.C =================================================================== --- /fact/tools/rootmacros/calscope_persistent.C (revision 12166) +++ /fact/tools/rootmacros/calscope_persistent.C (revision 12166) @@ -0,0 +1,153 @@ +#include +#include +#include +#include +#include +#include +#include "fits.h" + +int calscope_persistent(const char *name, const char *drsname, size_t pixelnr, UInt_t maxevents) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//Example call in ROOT: +//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348) +//T. Krähenbühl, August 2011, tpk@phys.ethz.ch +//****************************************************************************** + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +//Open the files +//------------------------------------------- + fits datafile(name); + fits drsfile(drsname); + + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + + if (!drsfile) + { + cout << "Couldn't properly open the drsfile." << endl; + return 1; + } + +//------------------------------------------- +//Print the headers +//------------------------------------------- + cout << "-------------------- Data Header -------------------" << endl; + datafile.PrintKeys(); + cout << "----------------- DRS Calib Header -----------------" << endl; + drsfile.PrintKeys(); + cout << "------------------- Data Columns -------------------" << endl; + datafile.PrintColumns(); + cout << "---------------- DRS Calib Columns -----------------" << endl; + drsfile.PrintColumns(); + cout << "--------------------- Data --------------------" << endl; + +//------------------------------------------- +//Get the sizes of the data columns +//------------------------------------------- + const size_t data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI + const size_t drs_n = drsfile.GetN("BaselineMean"); + if(drs_n!=data_n) + { + cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; + return 1; + } + + const UInt_t data_roi = datafile.GetUInt("NROI"); // Value from header + const UInt_t data_px = datafile.GetUInt("NPIX"); + +//------------------------------------------- +//Link the data to variables +//------------------------------------------- + unsigned int data_num; + datafile.SetRefAddress("EventNum", data_num); + vector data(data_n); + datafile.SetVecAddress("Data", data); + vector dataOffset(data_px); + datafile.SetVecAddress("StartCellData", dataOffset); + + vector drs_basemean(drs_n); + drsfile.SetVecAddress("BaselineMean", drs_basemean); + vector drs_gainmean(drs_n); + drsfile.SetVecAddress("GainMean", drs_gainmean); + vector drs_triggeroffsetmean(drs_n); + drsfile.SetVecAddress("TriggerOffsetMean", drs_triggeroffsetmean); + drsfile.GetRow(0); //Read the calibration data + +//------------------------------------------- +//Create the canvas, plot and title +//------------------------------------------- + char title[500]; + std::sprintf(title,"Data: %s, DRS: %s, Px %i",name,drsname,pixelnr); + TCanvas *canv = new TCanvas( "canv", "Mean values of the first event", 100, 10, 700, 500 ); + TProfile *pix = new TProfile("pix", title, 1024, -0.5, 1023.5); + pix->GetXaxis()->SetTitle("DRS bin (@2 GHz)"); + pix->GetYaxis()->SetTitle("Amplitude (mV)"); + + char ptitle[500]; + std::sprintf(ptitle,"Data: %s, DRS: %s, Px %i",name,drsname); + TCanvas *pcanv = new TCanvas( "pcanv", title, 800, 10, 700, 500 ); + TH2F *pers = new TH2F("pers",ptitle,1024,-0.5,1023.5,200,-49.5,150.5); + pers->GetXaxis()->SetTitle("DRS bin (@2 GHz)"); + pers->GetYaxis()->SetTitle("Amplitude (mV)"); + +//------------------------------------------- +//Iterate over the events +//------------------------------------------- +if(maxevents==0) maxevents = datafile.GetNumRows(); +// for (size_t i=0; iFill(k,sample); + pers->Fill(k,sample); + } +// } + } + +//------------------------------------------- +//Draw the data +//------------------------------------------- + canv->cd(); + pix->Draw(); + canv->Modified(); + canv->Update(); + pcanv->cd(); + pers->Draw("col"); + pcanv->Modified(); + pcanv->Update(); + return 0; +} Index: /fact/tools/rootmacros/factfir.C =================================================================== --- /fact/tools/rootmacros/factfir.C (revision 12166) +++ /fact/tools/rootmacros/factfir.C (revision 12166) @@ -0,0 +1,21 @@ +#include + +#ifndef FAD_MAX_SAMPLES + #warning FAD_MAX_SAMPLES not defined. defining: FAD_MAX_SAMPLES to 1024 + #define FAD_MAX_SAMPLES 1024 +#endif + +// source vector is +// float Ameas[FAD_MAX_SAMPLES]; +void factfir(double b, vector &a, int k, vector &source, vector &dest){ + //dest.clear(); + + for (int slice =0; slice < FAD_MAX_SAMPLES; slice++) { + float currentval = 0; + + for (int i=0; i < k ; i++){ + currentval += a[i] * source[ (slice - i + FAD_MAX_SAMPLES) % FAD_MAX_SAMPLES ]; + } + dest[slice] = currentval / b; + } +} Index: /fact/tools/rootmacros/fana.C =================================================================== --- /fact/tools/rootmacros/fana.C (revision 12166) +++ /fact/tools/rootmacros/fana.C (revision 12166) @@ -0,0 +1,444 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "fits.h" +//#include "TPKplotevent.c" +//#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" + +#include "zerosearch.C" + +#define NPIX 1440 +#define NCELL 1024 + +// data access and treatment +#define FAD_MAX_SAMPLES 1024 + +vector data; +vector data_offset; + +unsigned int data_num; + +size_t data_n; +UInt_t data_px; +UInt_t data_roi; +int NEvents; + +size_t drs_n; +vector drs_basemean; +vector drs_gainmean; +vector drs_triggeroffsetmean; + +int FOpenDataFile( fits & ); + + +vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors +vector N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors +vector Vcorr(FAD_MAX_SAMPLES); // corrected Values +vector Vdiff(FAD_MAX_SAMPLES); // numerical derivative + +vector Vslide(FAD_MAX_SAMPLES); // sliding average result +vector Vcfd(FAD_MAX_SAMPLES); // CDF result +vector Vcfd_slide(FAD_MAX_SAMPLES); // CDF result sliding averaged again + +#include "factfir.C" + +float getValue( int, int ); +void computeN1mean( int ); +void removeSpikes( int ); + +// histograms +const int Ntypes = 6; +const unsigned int tAmeas = 0, + tN1mean = 1, + tVcorr = 2, + tVslide = 3, + tVcfd = 4, + tVcfd_slide = 5; + + +TH1F* h; +TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts + // x = pixel id, y = DRS cell id +TH2F hPixelCellData("PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.); + +void BookHistos( int ); + +// Create a canvas +TCanvas* CW; + +int spikeDebug = 0; + + +int fana( const char *datafilename = "/media/daten_platte/FACT/data/20111009_013.fits", + const char *drsfilename = "/media/daten_platte/FACT/data/20111009_009.drs.fits", + int pixelnr = 0, + int firstevent = 0, + int nevents = -1 ) +{ +// sliding window filter settings + int k_slide = 16; + vector a_slide(k_slide, 1); + double b_slide = k_slide; + +// CFD filter settings + int k_cfd = 10; + vector a_cfd(k_cfd, 0); + double b_cfd = 1.; + a_cfd[0]=0.75; + a_cfd[k_cfd-1]=-1.; + + +// 2nd slinding window filter + int ks2 = 16; + vector as2(ks2, 1); + double bs2 = ks2; + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +// Open the file +//------------------------------------------- + fits datafile( datafilename ); + if (!datafile) { + printf( "Could not open the file: %s\n", datafilename ); + return 1; + } + + // access data + NEvents = FOpenDataFile( datafile ); + + printf("number of events in file: %d\n", NEvents); + + // compare the number of events in the data file with the nevents the + // the user would like to read. nevents = -1 means: read all + if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; + + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + + FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- + if(drs_n!=data_n) + { + cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; + return 1; + } +// Book the histograms + BookHistos( data_roi ); + +// Loop over events + cout << "--------------------- Data --------------------" << endl; + + float value; + float integral =0; + + TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5); + + for ( int ev = firstevent; ev < firstevent + nevents; ev++) { + + datafile.GetRow( ev ); + if (ev % 50 ==0) + cout << "Event number: " << data_num << endl; + + for ( int pix = 0; pix < 1440; pix++ ){ + hStartCell->Fill( pix, data_offset[pix] ); + } + + for ( unsigned int sl = 0; sl < data_roi; sl++){ + value = getValue(sl, pixelnr); + //printf("value = %f\n", value); + Ameas[ sl ] = value; + h[tAmeas].SetBinContent(sl, value); + + } + + computeN1mean( data_roi ); + removeSpikes( data_roi ); + + for ( unsigned int sl = 0; sl < data_roi; sl++){ + hPixelCellData.Fill(sl, Vcorr[ sl ]); + } + + // calculate sliding mean using FIR filter function + factfir(b_slide , a_slide, k_slide, Vcorr, Vslide); + // calculate CDF + factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); + // calculate i2nd sliding mean + factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd_slide); + + for ( unsigned int sl = 0; sl < data_roi; sl++){ + h[tVcfd].SetBinContent(sl, Vcfd[sl]); + h[tVslide].SetBinContent(sl, Vslide[sl]); + h[tVcfd_slide].SetBinContent(sl, Vcfd_slide[sl]); + } + + + vector * zeros = zerosearch(Vcfd_slide,-1,10,20); + if (zeros->size() ==0) + continue; + + // check value of Vside at zero position + for ( int i=0; isize(); i++){ + cout << zeros->at(i) << ":\t" << Vslide[zeros->at(i)]<Fill(Vslide[zeros->at(i)]); + } + +/* + + vector goodzeros; + if (zeros->size() > 1){ + for ( int i=0; isize(); i++){ + if (i==0) { + if ( abs(zeros->at(i) - zeros->at(i+1)) > 90 ) + goodzeros.push_back(zeros->at(i)); + } else if ( i==zeros->size()-1){ + if (abs(zeros->at(i) - zeros->at(i-1)) > 90) + goodzeros.push_back(zeros->at(i)); + } else { + if (abs(zeros->at(i) - zeros->at(i-1)) > 90 && + abs(zeros->at(i) - zeros->at(i+1)) > 90 ) + goodzeros.push_back(zeros->at(i)); + } + } + } else { + goodzeros.push_back(zeros->at(0)); + } + + // compute integrals around the good zeros. + // from -30 to +80 around it. + for (int i=0; i integral) + integral = Ameas[sl]; + } +// cout << "\t-->" << integral << endl; + sp->Fill(integral); + } + +*/ + + + if ( spikeDebug ){ + + CW->cd( tAmeas + 1); + gPad->SetGrid(); + h[tAmeas].Draw(); + + CW->cd( tN1mean + 1); + gPad->SetGrid(); + h[tN1mean].Draw(); + + CW->cd( tVcorr + 1); + gPad->SetGrid(); + h[tVcorr].Draw(); + + CW->cd( tVslide + 1); + gPad->SetGrid(); + h[tVslide].Draw(); + + CW->cd( tVcfd + 1); + gPad->SetGrid(); + h[tVcfd].Draw(); + TLine zeroline(0,0,1024,0); + zeroline.SetLineColor(kBlue); + zeroline.Draw(); + + CW->cd( tVcfd_slide + 1); + gPad->SetGrid(); + h[tVcfd_slide].Draw(); + zeroline.Draw(); + + CW->Update(); + + //Process gui events asynchronously during input + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + if (input=="q\n") break; + } + + + } + + TCanvas * cSpektrum = new TCanvas(); + cSpektrum->cd(); + sp->Draw(); + + + TCanvas * cStartCell = new TCanvas(); + cStartCell->cd(); + hStartCell->Draw(); + hPixelCellData.Draw(); + delete cStartCell; + return( 0 ); +} + +void removeSpikes(int Samples){ + + const float fract = 0.8; + float x, xp, xpp, x3p; + + // assume that there are no spikes + for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i]; + +// find the spike and replace it by mean value of neighbors + for ( int i = 0; i < Samples; i++) { + + // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] ); + + x = Ameas[i] - N1mean[i]; + + if ( x < -5. ){ // a spike candidate + // check consistency with a single channel spike + xp = Ameas[i+1] - N1mean[i+1]; + xpp = Ameas[i+2] - N1mean[i+2]; + x3p = Ameas[i+3] - N1mean[i+3]; + + // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p); + + if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){ + // printf("double spike candidate\n"); + Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.; + Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.; + // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]); + // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]); + i = i + 3; + } + else{ + + if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){ + Vcorr[i+1] = N1mean[i+1]; + // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]); + // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.); + N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.); + i = i + 2;//do not care about the next sample it was the spike + } + // treatment for the end of the pipeline must be added !!! + } + } + else{ + // do nothing + } + } // end of spike search and correction + for ( int i = 0; i < Samples; i++ ) h[ tVcorr ].SetBinContent( i, Vcorr[i] ); +} + +void computeN1mean( int Samples ){ + +// compute the mean of the left and right neighbors of a channel + + for( int i = 0; i < Samples; i++){ + if (i == 0){ // use right sample es mean + N1mean[i] = Ameas[i+1]; + } + else if ( i == Samples-1 ){ //use left sample as mean + N1mean[i] = Ameas[i-1]; + } + else{ + N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.; + } + h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]); + } +} // end of computeN1mean computation + +float getValue( int slice, int pixel ){ + + const float dconv = 2000/4096.0; + + float vraw, vcal; + + unsigned int pixel_pt; + unsigned int slice_pt; + unsigned int cal_pt; + unsigned int drs_cal_offset; + + // printf("pixel = %d, slice = %d\n", slice, pixel); + + pixel_pt = pixel * data_roi; + slice_pt = pixel_pt + slice; + drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi; + cal_pt = pixel_pt + drs_cal_offset; + + vraw = data[ slice_pt ] * dconv; + vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35; + + return( vcal ); +} + +void BookHistos( int Samples ){ +// booking and parameter settings for all histos + + h = new TH1F[ Ntypes ]; + + for ( int type = 0; type < Ntypes; type++){ + + h[ type ].SetBins(Samples, 0, Samples); + h[ type ].SetLineColor(1); + h[ type ].SetLineWidth(2); + + // set X axis paras + h[ type ].GetXaxis()->SetLabelSize(0.1); + h[ type ].GetXaxis()->SetTitleSize(0.1); + h[ type ].GetXaxis()->SetTitleOffset(1.2); + h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.)); + + // set Y axis paras + h[ type ].GetYaxis()->SetLabelSize(0.1); + h[ type ].GetYaxis()->SetTitleSize(0.1); + h[ type ].GetYaxis()->SetTitleOffset(0.3); + h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)"); + } + CW = new TCanvas("CW","DRS Waveform",10,10,800,600); + CW->Divide(1, Ntypes); + + hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024); + +} +int FOpenDataFile(fits &datafile){ + +/* cout << "-------------------- Data Header -------------------" << endl; + datafile.PrintKeys(); + cout << "------------------- Data Columns -------------------" << endl; + datafile.PrintColumns(); + */ + + //Get the size of the data column + data_roi = datafile.GetUInt("NROI"); // Value from header + data_px = datafile.GetUInt("NPIX"); + data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI + + //Set the sizes of the data vectors + data.resize(data_n,0); + data_offset.resize(data_px,0); + + //Link the data to variables + datafile.SetRefAddress("EventNum", data_num); + datafile.SetVecAddress("Data", data); + datafile.SetVecAddress("StartCellData", data_offset); + datafile.GetRow(0); + + cout << "Opening data file successful..." << endl; + + return datafile.GetNumRows() ; +} + Index: /fact/tools/rootmacros/fir.h =================================================================== --- /fact/tools/rootmacros/fir.h (revision 12166) +++ /fact/tools/rootmacros/fir.h (revision 12166) @@ -0,0 +1,8 @@ +#ifndef __FIR_H +#define __FIR_H + + #define FAD_MAX_SAMPLES 1024 +// extern float Vmeas[FAD_MAX_SAMPLES]; +// extern float Vfir[FAD_MAX_SAMPLES]; + +#endif Index: /fact/tools/rootmacros/fits.h =================================================================== --- /fact/tools/rootmacros/fits.h (revision 12166) +++ /fact/tools/rootmacros/fits.h (revision 12166) @@ -0,0 +1,657 @@ +#ifndef MARS_fits +#define MARS_fits + +#include + +#include +#include +#include +#include + +#ifdef __EXCEPTIONS +#include +#endif + +#ifdef __CINT__ +#define off_t size_t +#endif + +#ifndef __MARS__ +#include +#include +#include +#define gLog cerr +#define err "" +#define all "" +#else +#include "MLog.h" +#include "MLogManip.h" +#endif + +#ifdef HAVE_ZLIB +#include "izstream.h" +#else +#include +#define izstream ifstream +#warning Support for zipped FITS files disabled. +#endif + +namespace std +{ + +class fits : public izstream +{ +public: + struct Entry + { + char type; + string value; + string comment; + + template + T Get() const + { + T t; + + istringstream str(value); + str >> t; + + return t; + } + }; + + struct Table + { + off_t offset; + + string name; + size_t bytes_per_row; + size_t num_rows; + size_t num_cols; + + struct Column + { + size_t offset; + size_t num; + size_t size; + char type; + }; + + typedef map Keys; + typedef map Columns; + + Columns cols; + Keys keys; + + string Trim(const string &str, char c=' ') const + { + // Trim Both leading and trailing spaces + const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces + const size_t pend = str.find_last_not_of(c); // Find the first character position from reverse af + + // if all spaces or empty return an empty string + if (string::npos==pstart || string::npos==pend) + return string(); + + return str.substr(pstart, pend-pstart+1); + } + + bool Check(const string &key, char type, const string &value="") const + { + const Keys::const_iterator it = keys.find(key); + if (it==keys.end()) + { + ostringstream str; + str << "Key '" << key << "' not found."; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + if (it->second.type!=type) + { + ostringstream str; + str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << "."; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + if (!value.empty() && it->second.value!=value) + { + ostringstream str; + str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << "."; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + return true; + } + + Keys ParseBlock(const vector &vec) const + { + map rc; + + for (unsigned int i=0; i &vec, off_t off) : + offset(off), keys(ParseBlock(vec)) + { + if (!Check("XTENSION", 'T', "BINTABLE") || + !Check("NAXIS", 'I', "2") || + !Check("BITPIX", 'I', "8") || + !Check("PCOUNT", 'I', "0") || + !Check("GCOUNT", 'I', "1") || + !Check("EXTNAME", 'T') || + !Check("NAXIS1", 'I') || + !Check("NAXIS2", 'I') || + !Check("TFIELDS", 'I')) + return; + + bytes_per_row = Get("NAXIS1"); + num_rows = Get("NAXIS2"); + num_cols = Get("TFIELDS"); + + size_t bytes = 0; + for (size_t i=1; i<=num_cols; i++) + { + ostringstream num; + num << i; + + if (!Check("TTYPE"+num.str(), 'T') || + !Check("TFORM"+num.str(), 'T')) + return; + + const string id = Get("TTYPE"+num.str()); + const string fmt = Get("TFORM"+num.str()); + + istringstream sin(fmt); + int n = 0; + sin >> n; + if (!sin) + n = 1; + + const char type = fmt[fmt.length()-1]; + + size_t size = 0; + switch (type) + { + // We could use negative values to mark floats + // otheriwse we could just cast them to int64_t? + case 'L': size = 1; break; // logical + // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits) + case 'B': size = 1; break; // byte + case 'I': size = 2; break; // short + case 'J': size = 4; break; // int + case 'K': size = 8; break; // long long + case 'E': size = 4; break; // float + case 'D': size = 8; break; // double + // case 'C': size = 8; break; // complex float + // case 'M': size = 16; break; // complex double + // case 'P': size = 8; break; // array descriptor (32bit) + // case 'Q': size = 16; break; // array descriptor (64bit) + default: + { + ostringstream str; + str << "FITS format TFORM='" << fmt << "' not yet supported."; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return; +#endif + } + } + + const Table::Column col = { bytes, n, size, type }; + + cols[id] = col; + bytes += n*size; + } + + if (bytes!=bytes_per_row) + { +#ifdef __EXCEPTIONS + throw runtime_error("Column size mismatch"); +#else + gLog << err << "ERROR - Column size mismatch" << endl; + return; +#endif + } + + name = Get("EXTNAME"); + } + + void PrintKeys() const + { + for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++) + gLog << all << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl; + } + + void PrintColumns() const + { + typedef map, Column> Sorted; + + Sorted sorted; + + for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++) + sorted[make_pair(it->second.offset, it->first)] = it->second; + + for (Sorted::const_iterator it=sorted.begin(); it!=sorted.end(); it++) + gLog << all << setw(6) << it->second.offset << "| " << it->first.second << " [" << it->second.type << '=' << it->second.num << 'x' << it->second.size << "byte]" << endl; + } + + operator bool() const { return !name.empty(); } + + bool HasKey(const string &key) const + { + return keys.find(key)!=keys.end(); + } + + // Values of keys are always signed + template + T Get(const string &key) const + { + const map::const_iterator it = keys.find(key); + if (it==keys.end()) + { + ostringstream str; + str << "Key '" << key << "' not found." << endl; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return 0; +#endif + } + return it->second.Get(); + } + + size_t GetN(const string &key) const + { + const Columns::const_iterator it = cols.find(key); + return it==cols.end() ? 0 : it->second.num; + } + }; + +private: + Table fTable; + + typedef pair Address; + typedef vector
Addresses; + //map fAddresses; + Addresses fAddresses; + + vector fBufferRow; + vector fBufferDat; + + size_t fRow; + + vector ReadBlock(vector &vec) + { + bool endtag = false; + for (int i=0; i<36; i++) + { + char c[81]; + c[80] = 0; + read(c, 80); + if (!good()) + break; + + if (c[0]==0) + return vector(); + + string str(c); + +// if (!str.empty()) +// cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl; + + if (str=="END ") + { + endtag = true; + + // Make sure that no empty vector is returned + if (vec.size()%36==0) + vec.push_back(string("END = '' / ")); + } + + if (endtag) + continue; + + vec.push_back(str); + } + + return vec; + } + + string Compile(const string &key, int16_t i=-1) + { + if (i<0) + return key; + + ostringstream str; + str << key << i; + return str.str(); + } + +public: + fits(const string &fname) : izstream(fname.c_str()) + { + char simple[10]; + read(simple, 10); + if (!good()) + return; + + if (memcmp(simple, "SIMPLE = ", 10)) + { + clear(rdstate()|ios::badbit); +#ifdef __EXCEPTIONS + throw runtime_error("File is not a FITS file."); +#else + gLog << err << "ERROR - File is not a FITS file." << endl; + return; +#endif + } + + seekg(0); + + while (good()) + { + vector block; + while (1) + { + // FIXME: Set limit on memory consumption + ReadBlock(block); + if (!good()) + { + clear(rdstate()|ios::badbit); +#ifdef __EXCEPTIONS + throw runtime_error("FITS file corrupted."); +#else + gLog << err << "ERROR - FITS file corrupted." << endl; + return; +#endif + } + + if (block.size()%36) + break; + } + + if (block.size()==0) + break; + + if (block[0].substr(0, 9)=="SIMPLE =") + continue; + + if (block[0].substr(0, 9)=="XTENSION=") + { + // FIXME: Check for table name + + fTable = Table(block, tellg()); + fRow = (size_t)-1; + + //fTable.PrintKeys(); + + if (!fTable) + { + clear(rdstate()|ios::badbit); + return; + } + + fBufferRow.resize(fTable.bytes_per_row); + fBufferDat.resize(fTable.bytes_per_row); + + /* + // Next table should start at: + const size_t size = fTable.bytes_per_row*fTable.num_rows; + const size_t blks = size/(36*80); + const size_t rest = size%(36*80); + + seekg((blks+(rest>0?1:0))*(36*80), ios::cur); + if (!good()) + gLog << err << "File seems to be incomplete (less data than expected from header)." << endl; + + fRow = fTable.num_rows; + */ + + break; + } + } + } + + void ReadRow(size_t row) + { + // if (row!=fRow+1) // Fast seeking is ensured by izstream + seekg(fTable.offset+row*fTable.bytes_per_row); + + fRow = row; + + read(fBufferRow.data(), fBufferRow.size()); + //fin.clear(fin.rdstate()&~ios::eofbit); + } + + template + void revcpy(char *dest, const char *src, int num) + { + const char *pend = src + num*N; + for (const char *ptr = src; ptrsecond; + + const char *src = fBufferRow.data() + c.offset; + char *dest = reinterpret_cast(it->first); + + // Let the compiler do some optimization by + // knowing the we only have 1, 2, 4 and 8 + switch (c.size) + { + case 1: memcpy (dest, src, c.num*c.size); break; + case 2: revcpy<2>(dest, src, c.num); break; + case 4: revcpy<4>(dest, src, c.num); break; + case 8: revcpy<8>(dest, src, c.num); break; + } + } + + return good(); + } + + bool GetNextRow() + { + return GetRow(fRow+1); + } + + bool SkipNextRow() + { + seekg(fTable.offset+(++fRow)*fTable.bytes_per_row); + return good(); + } + + static bool Compare(const Address &p1, const Address &p2) + { + return p1.first>p2.first; + } + + template + bool SetPtrAddress(const string &name, T *ptr, size_t cnt) + { + if (fTable.cols.count(name)==0) + { + ostringstream str; + str <<"SetPtrAddress('" << name << "') - Column not found." << endl; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + if (sizeof(T)!=fTable.cols[name].size) + { + ostringstream str; + str << "SetPtrAddress('" << name << "') - Element size mismatch: expected " + << fTable.cols[name].size << " from header, got " << sizeof(T) << endl; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + if (cnt!=fTable.cols[name].num) + { + ostringstream str; + str << "SetPtrAddress('" << name << "') - Element count mismatch: expected " + << fTable.cols[name].num << " from header, got " << cnt << endl; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + // if (fAddresses.count(ptr)>0) + // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; + + //fAddresses[ptr] = fTable.cols[name]; + fAddresses.push_back(make_pair(ptr, fTable.cols[name])); + sort(fAddresses.begin(), fAddresses.end(), Compare); + return true; + } + + template + bool SetRefAddress(const string &name, T &ptr) + { + return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T)); + } + + template + bool SetVecAddress(const string &name, vector &vec) + { + return SetPtrAddress(name, vec.data(), vec.size()); + } + + template + T Get(const string &key) const + { + return fTable.Get(key); + } + + bool SetPtrAddress(const string &name, void *ptr) + { + if (fTable.cols.count(name)==0) + { + ostringstream str; + str <<"SetPtrAddress('" << name << "') - Column not found." << endl; +#ifdef __EXCEPTIONS + throw runtime_error(str.str()); +#else + gLog << err << "ERROR - " << str.str() << endl; + return false; +#endif + } + + // if (fAddresses.count(ptr)>0) + // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; + + //fAddresses[ptr] = fTable.cols[name]; + fAddresses.push_back(make_pair(ptr, fTable.cols[name])); + sort(fAddresses.begin(), fAddresses.end(), Compare); + return true; + } + + bool HasKey(const string &key) const { return fTable.HasKey(key); } + int64_t GetInt(const string &key) const { return fTable.Get(key); } + uint64_t GetUInt(const string &key) const { return fTable.Get(key); } + double GetFloat(const string &key) const { return fTable.Get(key); } + string GetStr(const string &key) const { return fTable.Get(key); } + + size_t GetN(const string &key) const + { + return fTable.GetN(key); + } + + size_t GetNumRows() const { return fTable.num_rows; } + size_t GetRow() const { return fRow; } + + operator bool() const { return fTable && fTable.offset!=0; } + + void PrintKeys() const { fTable.PrintKeys(); } + void PrintColumns() const { fTable.PrintColumns(); } +}; + +}; +#endif Index: /fact/tools/rootmacros/gainanalysis.C =================================================================== --- /fact/tools/rootmacros/gainanalysis.C (revision 12166) +++ /fact/tools/rootmacros/gainanalysis.C (revision 12166) @@ -0,0 +1,273 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "fits.h" +#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" +#include "FPedestal.c" +#include "FOscilloscopeAllPx.c" + +#define n_peaks 2 + +int gainanalysis(const char *name, const char *drsname, size_t pixelnr) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//****************************************************************************** + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +//Define the data variables +//------------------------------------------- + vector data; + vector data_offset; + unsigned int data_num; + size_t data_n; + UInt_t data_px; + UInt_t data_roi; + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + size_t drs_n; + vector drs_basemean; + vector drs_gainmean; + vector drs_triggeroffsetmean; + FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- +// if(drs_n!=data_n) +// { +// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; +// return 1; +// } + +//------------------------------------------- +//Open the file +//------------------------------------------- + fits datafile(name); + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px); + +//------------------------------------------- +//Create the canvas, plot and title +//------------------------------------------- + char title_spect[500]; + char name_spect[50]; + std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname); + TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 ); + TH1F* spectrum[data_px]; + for(int i=0; iGetXaxis()->SetTitle("Pulse size (a.u.)"); + spectrum[i]->GetYaxis()->SetTitle("Counts"); + } + + TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 ); +//------------------------------------------- +//Define the fit +//------------------------------------------- + int fit_start = 20; + int fit_end = 200; + 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-[7])/[8])",fit_start,fit_end); + Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30}; + Double_t parread_tmp[5+2*n_peaks]; + Double_t parread[data_px][5+2*n_peaks]; +// 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); +// Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0}; + f_peaks->SetParameters(par_init); + f_peaks->SetLineColor(2); + + TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end); + f_pedestal->SetLineColor(2); + f_pedestal->SetLineWidth(1); + + TF1* f_peak[n_peaks]; + for(int i=0; iSetLineColor(2); + f_peak[i]->SetLineWidth(1); + } + +//------------------------------------------- +//Find the baseline +//------------------------------------------- + float pedestal_mean, pedestal_rms; +// FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms); + //Return values + //Value of maximal probability: -2.225 +- 4.3144 + pedestal_mean = -2.225; + pedestal_rms = 4.3144; + +//------------------------------------------- +//Oscilloscope +//------------------------------------------- +// float threshold = pedestal_mean+pedestal_rms; + float threshold = 2.5; + +//*********************** + FOscilloscopeAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum); + + canv_spect->cd(); + spectrum[0]->Draw(); + canv_spect->Modified(); + canv_spect->Update(); +//*********************** +// fits datafile2("20110804_025.fits"); +// if (!datafile2) +// { +// cout << "Couldn't properly open the datafile2." << endl; +// return 1; +// } +// FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px); +// FOscilloscopeAllPx(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);// +// canv_spect->cd(); +// spectrum[0]->Draw(); +// canv_spect->Modified(); +// canv_spect->Update(); +////*********************** +// fits datafile3("20110804_026.fits"); +// if (!datafile3) +// { +// cout << "Couldn't properly open the datafile3." << endl; +// return 1; +// } +// FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px); +// FOscilloscopeAllPx(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum); +// +// canv_spect->cd(); +// spectrum[0]->Draw(); +// canv_spect->Modified(); +// canv_spect->Update(); +////*********************** +// fits datafile4("20110804_027.fits"); +// if (!datafile4) +// { +// cout << "Couldn't properly open the datafile4." << endl; +// return 1; +// } +// FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px); +// FOscilloscopeAllPx(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum); +// +// canv_spect->cd(); +// spectrum[0]->Draw(); +// canv_spect->Modified(); +// canv_spect->Update(); +//*********************** + Double_t parread_dist[data_px]; + Double_t parread_x[data_px]; + for(int i=0; iSetParameters(par_init); + if(!(spectrum[i]->Fit(f_peaks,"NR"))) { +// std::cout << "Fit successful." << std::endl; + f_peaks->GetParameters(&parread_tmp[0]); + for(int j=0; j<5+2*n_peaks; j++) { + parread[i][j] = parread_tmp[j]; + } + parread_dist[i] = parread_tmp[1]; + parread_x[i] = i; + } + else { + for(int j=0; j<5+2*n_peaks; j++) { + parread[i][j] = 0; + } + parread_dist[i] = 0; + parread_x[i] = i; + } + } + + canv_gain->cd(); + TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist); + plot_gain->GetYaxis()->SetRangeUser(0,150); + plot_gain->Draw("A*"); + canv_gain->Modified(); + canv_gain->Update(); + + canv_spect->cd(); + canv_spect->SetLogy(); + +//------------------------------------------- +//Draw the data where the gain is out of limits +//------------------------------------------- + char temp = 'x'; + std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl; + for(int i=0; i105)) { + std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl; + spectrum[i]->Draw(); + + f_peaks->SetParameters(&parread[i][0]); + f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]); + f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]); +// f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]); +// f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]); +// f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]); + f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]); + + f_peaks->Draw("same"); + for(int i=0; iDraw("same"); + } + f_pedestal->Draw("same"); + + canv_spect->Modified(); + canv_spect->Update(); + temp='x'; + while ((temp!='\n') && (temp!='a')) + cin.get(temp); + if(temp=='a') break; + } + } +//------------------------------------------- +//Draw the data +//------------------------------------------- + canv_spect->cd(); + Int_t plotspect_nr; + std::cout << "Enter the number of the spectrum to plot (-1 to quit): "; + std::cin >> plotspect_nr; + while(plotspect_nr>=0) { + plotspect_nr = plotspect_nr % data_px; + spectrum[plotspect_nr]->Draw(); + + f_peaks->SetParameters(&parread[plotspect_nr][0]); + f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]); + f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]); +// f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]); +// f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]); +// f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]); + f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]); + + f_peaks->Draw("same"); + for(int i=0; iDraw("same"); + } + f_pedestal->Draw("same"); + + canv_spect->Modified(); + canv_spect->Update(); + + std::cout << "Enter the number of the spectrum to plot (-1 to quit): "; + std::cin >> plotspect_nr; + } + + return 0; +} Index: /fact/tools/rootmacros/peakfinderscope.C =================================================================== --- /fact/tools/rootmacros/peakfinderscope.C (revision 12166) +++ /fact/tools/rootmacros/peakfinderscope.C (revision 12166) @@ -0,0 +1,231 @@ +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "fits.h" +#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" +#include "FPedestal.c" +#include "FOscilloscope.c" + +#define n_peaks 4 + +int peakfinderscope(const char *name, const char *drsname, size_t eventnr, size_t pixelnr) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//****************************************************************************** + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +//Define the data variables +//------------------------------------------- + vector data; + vector data_offset; + unsigned int data_num; + size_t data_n; + UInt_t data_px; + UInt_t data_roi; + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + size_t drs_n; + vector drs_basemean; + vector drs_gainmean; + vector drs_triggeroffsetmean; + FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- +// if(drs_n!=data_n) +// { +// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; +// return 1; +// } + +//------------------------------------------- +//Create the canvas, plot and title +//------------------------------------------- + char title[500]; + std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr); + TCanvas *canv = new TCanvas( "canv", "Waveform", 200, 10, 700, 500 ); +// TProfile *pulseshape = new TProfile("pulseshape", title, 70, -10.5, 59.5); + TH2F *pulseshape = new TH2F("pulseshape",title,70,-10.5,59.5,350,-10.5,59.5); + pulseshape->GetXaxis()->SetTitle("Time after trigger (bins @2 GHz)"); + pulseshape->GetYaxis()->SetTitle("Amplitude (mV)"); + + char title_base[500]; + std::sprintf(title_base,"All samples of Px %i (data: %s, DRS: %s)",pixelnr,name,drsname); + + char title_spect[500]; + std::sprintf(title_spect,"Spectrum of all pixels (data: %s, DRS: %s)",name,drsname); + TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 ); + TH1F *spectrum = new TH1F("spectrum", title_spect,300,0,600); + spectrum->GetXaxis()->SetTitle("Pulse size (a.u.)"); + spectrum->GetYaxis()->SetTitle("Counts"); + +//------------------------------------------- +//Define the fit +//------------------------------------------- + int fit_start = 30; + int fit_end = 380; + 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); + Double_t par[5+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,10000,30,30}; +// 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); +// Double_t par[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0}; + f_peaks->SetParameters(par); + f_peaks->SetLineColor(2); + + TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end); + f_pedestal->SetLineColor(2); + f_pedestal->SetLineWidth(1); + + TF1* f_peak[n_peaks]; + for(int i=0; iSetLineColor(2); + f_peak[i]->SetLineWidth(1); + } + +//------------------------------------------- +//Open the file +//------------------------------------------- + fits datafile(name); + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px); + +//------------------------------------------- +//Find the baseline +//------------------------------------------- + float pedestal_mean, pedestal_rms; +// FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms); + //Return values + //Value of maximal probability: -2.225 +- 4.3144 + pedestal_mean = -2.225; + pedestal_rms = 4.3144; + +//------------------------------------------- +//Oscilloscope +//------------------------------------------- +// float threshold = pedestal_mean+pedestal_rms; + float threshold = 2.5; + +//*********************** + FOscilloscope(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum); + canv->cd(); + pulseshape->Draw(); + canv->Modified(); + canv->Update(); + + canv_spect->cd(); + spectrum->Draw(); + canv_spect->Modified(); + canv_spect->Update(); +//*********************** + fits datafile2("20110804_025.fits"); + if (!datafile2) + { + cout << "Couldn't properly open the datafile2." << endl; + return 1; + } + FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px); + FOscilloscope(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum); + + canv->cd(); + pulseshape->Draw(); + canv->Modified(); + canv->Update(); + + canv_spect->cd(); + spectrum->Draw(); + canv_spect->Modified(); + canv_spect->Update(); +//*********************** + fits datafile3("20110804_026.fits"); + if (!datafile3) + { + cout << "Couldn't properly open the datafile3." << endl; + return 1; + } + FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px); + FOscilloscope(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum); + + canv->cd(); + pulseshape->Draw(); + canv->Modified(); + canv->Update(); + + canv_spect->cd(); + spectrum->Draw(); + canv_spect->Modified(); + canv_spect->Update(); +//*********************** + fits datafile4("20110804_027.fits"); + if (!datafile4) + { + cout << "Couldn't properly open the datafile4." << endl; + return 1; + } + FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px); + FOscilloscope(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum); + + canv->cd(); + pulseshape->Draw(); + canv->Modified(); + canv->Update(); + + canv_spect->cd(); + spectrum->Draw(); + canv_spect->Modified(); + canv_spect->Update(); +//*********************** + spectrum->Fit(f_peaks,"R+"); + f_peaks->GetParameters(&par[0]); + + f_peak[0]->SetParameters(par[2],par[0],par[3]); + f_peak[1]->SetParameters(par[4],par[0]+par[1],par[5]); + f_peak[2]->SetParameters(par[6],par[0]+2*par[1],par[7]); + f_peak[3]->SetParameters(par[8],par[0]+3*par[1],par[9]); +// f_peak[4]->SetParameters(par[10],par[12],par[11]); + f_pedestal->SetParameters(par[10],par[11],par[12]); + + std::cout << "Crosstalk propability approx. " << ((par[4]*par[5])/(par[2]*par[3]+par[4]*par[5])) << std::endl; +//------------------------------------------- +//Draw the data +//------------------------------------------- + canv->cd(); + pulseshape->Draw(); + canv->Modified(); + canv->Update(); + + canv_spect->cd(); + canv_spect->SetLogy(); + spectrum->Draw(); + f_peaks->Draw("same"); + + for(int i=0; iDraw("same"); + } + + f_pedestal->Draw("same"); + + canv_spect->Modified(); + canv_spect->Update(); + + return 0; +} Index: /fact/tools/rootmacros/poissonanalysis.C =================================================================== --- /fact/tools/rootmacros/poissonanalysis.C (revision 12166) +++ /fact/tools/rootmacros/poissonanalysis.C (revision 12166) @@ -0,0 +1,237 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "fits.h" +#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" +#include "FPedestal.c" +#include "FIntFixedPosAllPx.c" + +#define n_peaks 2 + +int poissonanalysis(const char *name, const char *drsname, UInt_t integration_delay, UInt_t integration_size, UInt_t spect_max ) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//****************************************************************************** + + gROOT->SetStyle("Plain"); +// TFile f("spectra.root","RECREATE"); +//------------------------------------------- +//Define the data variables +//------------------------------------------- + vector data; + vector data_offset; + unsigned int data_num; + size_t data_n; + UInt_t data_px; + UInt_t data_roi; + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + size_t drs_n; + vector drs_basemean; + vector drs_gainmean; + vector drs_triggeroffsetmean; + FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- +// if(drs_n!=data_n) +// { +// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; +// return 1; +// } + +//------------------------------------------- +//Open the file +//------------------------------------------- + fits datafile(name); + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px); + +//------------------------------------------- +//Create the canvas, plot and title +//------------------------------------------- + char title_spect[500]; + char name_spect[50]; +// std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname); + TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 ); + TH1F* spectrum[data_px]; + for(int i=0; iGetXaxis()->SetTitle("Pulse size (a.u.)"); + spectrum[i]->GetYaxis()->SetTitle("Counts"); + } + +// TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 ); +//------------------------------------------- +//Define the fit +//------------------------------------------- +// int fit_start = 20; +// int fit_end = 200; +// 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-[7])/[8])",fit_start,fit_end); +// Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30}; +// Double_t parread_tmp[5+2*n_peaks]; +// Double_t parread[data_px][5+2*n_peaks]; +//// 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); +//// Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0}; +// f_peaks->SetParameters(par_init); +// f_peaks->SetLineColor(2); +// +// TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end); +// f_pedestal->SetLineColor(2); +// f_pedestal->SetLineWidth(1); +// +// TF1* f_peak[n_peaks]; +// for(int i=0; iSetLineColor(2); +// f_peak[i]->SetLineWidth(1); +// } + +//------------------------------------------- +//Find the baseline +//------------------------------------------- +// float pedestal_mean, pedestal_rms; +//// FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms); +// //Return values +// //Value of maximal probability: -2.225 +- 4.3144 +// pedestal_mean = -2.225; +// pedestal_rms = 4.3144; + +//------------------------------------------- +//Oscilloscope +//------------------------------------------- +// float threshold = pedestal_mean+pedestal_rms; +// float threshold = 2.5; + +//*********************** +// UInt_t integration_size = 5; +// UInt_t integration_delay = 238; + + FIntFixedPosAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, integration_size, integration_delay, spectrum); + + canv_spect->cd(); + spectrum[0]->Draw(); + canv_spect->Modified(); + canv_spect->Update(); + +// Double_t parread_dist[data_px]; +// Double_t parread_x[data_px]; +// for(int i=0; iSetParameters(par_init); +// if(!(spectrum[i]->Fit(f_peaks,"NR"))) { +//// std::cout << "Fit successful." << std::endl; +// f_peaks->GetParameters(&parread_tmp[0]); +// for(int j=0; j<5+2*n_peaks; j++) { +// parread[i][j] = parread_tmp[j]; +// } +// parread_dist[i] = parread_tmp[1]; +// parread_x[i] = i; +// } +// else { +// for(int j=0; j<5+2*n_peaks; j++) { +// parread[i][j] = 0; +// } +// parread_dist[i] = 0; +// parread_x[i] = i; +// } +// } +// +// canv_gain->cd(); +// TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist); +// plot_gain->GetYaxis()->SetRangeUser(0,150); +// plot_gain->Draw("A*"); +// canv_gain->Modified(); +// canv_gain->Update(); +// +// canv_spect->cd(); +// canv_spect->SetLogy(); + +//------------------------------------------- +//Draw the data where the gain is out of limits +//------------------------------------------- +// char temp = 'x'; +// std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl; +// for(int i=0; i105)) { +// std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl; +// spectrum[i]->Draw(); +// +// f_peaks->SetParameters(&parread[i][0]); +// f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]); +// f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]); +//// f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]); +//// f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]); +//// f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]); +// f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]); +// +// f_peaks->Draw("same"); +// for(int i=0; iDraw("same"); +// } +// f_pedestal->Draw("same"); +// +// canv_spect->Modified(); +// canv_spect->Update(); +// temp='x'; +// while ((temp!='\n') && (temp!='a')) +// cin.get(temp); +// if(temp=='a') break; +// } +// } +//------------------------------------------- +//Draw the data +//------------------------------------------- + canv_spect->cd(); + Int_t plotspect_nr; + std::cout << "Enter the number of the spectrum to plot (-1 to quit): "; + std::cin >> plotspect_nr; + while(plotspect_nr>=0) { + plotspect_nr = plotspect_nr % data_px; + spectrum[plotspect_nr]->Draw(); + +// f_peaks->SetParameters(&parread[plotspect_nr][0]); +// f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]); +// f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]); +//// f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]); +//// f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]); +//// f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]); +// f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]); +// +// f_peaks->Draw("same"); +// for(int i=0; iDraw("same"); +// } +// f_pedestal->Draw("same"); + + canv_spect->Modified(); + canv_spect->Update(); + + std::cout << "Enter the number of the spectrum to plot (-1 to quit): "; + std::cin >> plotspect_nr; + } + +// f->Write(); + return 0; +} Index: /fact/tools/rootmacros/test.C =================================================================== --- /fact/tools/rootmacros/test.C (revision 12166) +++ /fact/tools/rootmacros/test.C (revision 12166) @@ -0,0 +1,47 @@ +#include "fits.h" + +#include + +int test(const char *name) +{ + fits file(name); + if (!file) + { + cout << "Couldn't properly open the file." << endl; + return 1; + } + + cout << "-------------------- Header -------------------" << endl; + file.PrintKeys(); + cout << "------------------- Columns -------------------" << endl; + file.PrintColumns(); + + cout << "--------------------- Data --------------------" << endl; + + unsigned int num; + short type; + + const size_t n = file.GetN("Data"); // Size of column "Data" + const UInt_t roi = file.GetUInt("NROI"); // Value from header + + vector data(n); + + file.SetRefAddress("EventNum", num); + file.SetRefAddress("TriggerType", type); + file.SetVecAddress("Data", data); + + for (size_t i=0; i +#include +#include + +#include +#include + +#include "fits.h" +#include "FOpenDataFile.c" +#include "FOpenCalibFile.c" +#include "FPedestalAllPx.c" + +int testFPedestalAllPx(const char *name, const char *drsname) +{ +//****************************************************************************** +//Read a datafile and plot the DRS-calibrated data +//ATTENTION: only works for ROI=1024 +// (array indices of the calibration wrong otherwise) +//Example call in ROOT: +//root [74] .x calscope.C++("20110804_024.fits","20110804_023.drs.fits",10,1348) +//T. Krähenbühl, August 2011, tpk@phys.ethz.ch +//****************************************************************************** + + gROOT->SetStyle("Plain"); + +//------------------------------------------- +//Open the file +//------------------------------------------- + fits datafile(name); + if (!datafile) + { + cout << "Couldn't properly open the datafile." << endl; + return 1; + } + +//------------------------------------------- +//Get the data +//------------------------------------------- + vector data; + vector data_offset; + unsigned int data_num; + size_t data_n; + UInt_t data_px; + UInt_t data_roi; + FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px); + +//------------------------------------------- +//Get the DRS calibration +//------------------------------------------- + size_t drs_n; + vector drs_basemean; + vector drs_gainmean; + vector drs_triggeroffsetmean; + FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n); + +//------------------------------------------- +//Check the sizes of the data columns +//------------------------------------------- + if(drs_n!=data_n) + { + cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl; + return 1; + } + + +//------------------------------------------- +//Create the title +//------------------------------------------- +// char title[500]; +// std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr); + +//------------------------------------------- +//Get the pedestals & RMS +//------------------------------------------- + float pedestal_mean[data_px]; + float pedestal_rms[data_px]; + FPedestalAllPx(datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_px, data_roi, pedestal_mean, pedestal_rms); +// vector calevent(data_px*data_roi); //Vector for the calibrated event +// FGetCalEvent(data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, calevent, data_px, data_roi); + +//------------------------------------------- +//Draw the data +//------------------------------------------- + TCanvas *canv_mean= new TCanvas( "canv_mean", "mean in mV", 30, 300, 900, 500 ); + TH1F* histo_mean = new TH1F("histo_mean","Value of maximal probability",2000,-99.5,100.5); + TCanvas *canv_rms= new TCanvas( "canv_rms", "RMS in mV", 960, 300, 900, 500 ); + TH1F* histo_rms = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5); + + for (int j=0; jFill(pedestal_mean[j]); + histo_rms->Fill(pedestal_rms[j]); + } +// histo_mean->GetYaxis()->SetRangeUser(0,150); + + canv_mean->cd(); + histo_mean->Draw(); + canv_mean->Modified(); + canv_mean->Update(); + canv_rms->cd(); + histo_rms->Draw(); + canv_rms->Modified(); + canv_rms->Update(); + + return 0; +} Index: /fact/tools/rootmacros/testfir.C =================================================================== --- /fact/tools/rootmacros/testfir.C (revision 12166) +++ /fact/tools/rootmacros/testfir.C (revision 12166) @@ -0,0 +1,62 @@ +#include "TRandom2.h" +#include "TH1F.h" +#include "TCanvas.h" + +#include +#include + + +#define FAD_MAX_SAMPLES 1024 + +vector Vmeas(FAD_MAX_SAMPLES); +vector Vfir(FAD_MAX_SAMPLES); +vector Vpulse(FAD_MAX_SAMPLES); +#include "factfir.C" + + +#define k 16 +double a[k]; +double b=k; + +void testfir(){ + + for (int j=0; jSetBinContent(i,Vmeas[i]); + fir->SetBinContent(i,Vfir[i]); + pulse->SetBinContent(i,Vpulse[i]); + } + + TCanvas *c = new TCanvas(); + c->Divide(1,3); + c->cd(1); + meas->Draw("HIST"); + c->cd(2); + fir->Draw("HIST"); + c->cd(3); + pulse->Draw("HIST"); + + +} Index: /fact/tools/rootmacros/zerosearch.C =================================================================== --- /fact/tools/rootmacros/zerosearch.C (revision 12166) +++ /fact/tools/rootmacros/zerosearch.C (revision 12166) @@ -0,0 +1,44 @@ +#include + +// searches for zero crossings in a given vector of floats +// for zero crossings in falling edge +// give edge = -1 +// +// +// returns pointer to vetctor of ints + +vector *zerosearch(vector &input, int edge =-1 , int pre=10, int post=30 ){ + + + vector * zeroPositions = new vector; + + for (int sl =pre ; sl < input.size()-post; sl++){ +//cout << "sl:" << sl << endl; + if (input[sl] * input[sl-1] < 0 || input[sl]==0.0 ){ // sign change --> zero crossing OR really zero + + if ( (input[sl-1] - input[sl]) * edge < 0){ // this is the crossing the user wanted + + // check if we go lower than a certain limit in the next few slices + for ( int lala=0; lalapush_back(sl); + sl += lala; + break; + } + } + + } else if ( (input[sl-1] - input[sl]) * edge > 0){ // this is the zero x-ing the user did not want + + // do nothing + + } else { // sl and sl-1 are equal .. don't know waht do to... + + } + + } + + } // end of loop over slices - between pre and post + + + return zeroPositions; +}