Index: fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C (revision 13465) @@ -0,0 +1,1890 @@ +/********************** FPulseTemplate *********************** +* read FACT raw data +* remove spikes +* calculate baseline +* find peaks (CFD and normal discriminator) +* compute pulse height and pulse integral spektrum of the peaks ++ search for Peaks in data ++ put peaks into different histos depending und Amplitude ++ draw histos for single, double, tripple, ....above 100mV Photonpulses ++ draw Parameterdevelopment depending on photon quantity ++ form a template (shape depending on number of photons) + so it can be used for detection of other peaks +*************************************************************/ + +//---------------------------------------------------------------------------- +// root libraries +//---------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +using namespace std; + +#define NPIX 1440 +#define NCELL 1024 +#define FAD_MAX_SAMPLES 1024 +#define MAX_PULS_ORDER 1 +#define HAVE_ZLIB + +//---------------------------------------------------------------------------- +// rootmacros +//---------------------------------------------------------------------------- + +#include "../fits.h" + +#include "../openFits.h" +#include "../openFits.c" + +#include "../discriminator.h" +#include "../discriminator.C" +#include "../zerosearch.h" +#include "../zerosearch.C" +#include "../factfir.C" + +#include "../DrsCalibration.C" +#include "../DrsCalibration.h" + +#include "../SpikeRemoval.h" +#include "../SpikeRemoval.C" + +//---------------------------------------------------------------------------- +// Decleration of global variables +//---------------------------------------------------------------------------- + +bool breakout=false; + +int gNEvents; +vector AllPixelDataVector; +vector StartCellVector; +unsigned int CurrentEventID; +size_t PXLxROI; +UInt_t gRegionOfInterest; +UInt_t NumberOfPixels; +TString histotitle; + +size_t TriggerOffsetROI, RC; +size_t drs_n; +vector drs_basemean; +vector drs_gainmean; +vector drs_triggeroffsetmean; + +vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector Vcorr(FAD_MAX_SAMPLES); // corrected Values +vector Vslide(FAD_MAX_SAMPLES); // sliding average result +vector Vcfd(FAD_MAX_SAMPLES); // CDF result +vector Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average + +float gGainMean = 9; +float gBSLMean = 0; + +typedef struct{ + double maxAmpl; + int countOfMax; + } OverlayMaximum; + +//typedef struct{ +// TString name; +// TString title; +// TString xTitle; +// TString yTitle; +// float yMax; +// float yMin; +//} histAttrib_t; + +//histAttrib_t* gHistoAttrib[MAX_PULS_ORDER]; +//histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER]; +//histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER]; +//histAttrib_t* gProfileAttrib[MAX_PULS_ORDER]; + + +// histograms +const int Number_Of_Debug_Histo_Types = 7; + +const unsigned int + Ameas_ = 0, + N1mean_ = 1, + Vcorr_ = 2, + Vtest_ = 3, + Vslide_ = 4, + Vcfd_ = 5, + Vcfd2_ = 6; + +//const char* gHistoTypes[8] = { +// "hPixelOverlay", +// "hPixelMax", +// "hPixelEdgeOverlay", +// "hPixelProfile", +// "hAllPixelOverlay", +// "hAllPixelMax", +// "hAllPixelMaxGaus", +// "hAllPixelProfile" +//}; + +//---------------------------------------------------------------------------- +// Initialisation of histograms +//---------------------------------------------------------------------------- + +// Temporary Objects +TH1F* debugHistos = NULL; +TH2F* hOverlayTemp = NULL; +TH1F* hMaximumTemp = NULL; +TH1D* hProjPeak = NULL; +TH1F* hTesthisto = NULL; +TH2F* hTesthisto2 = NULL; + +//Pixelwise Histograms +TH2F* hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks +TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL}; +TProfile* hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks +TProfile* hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks +TH1F* hPixelMax[MAX_PULS_ORDER] = {NULL}; +TH1F* hPixelMedian[MAX_PULS_ORDER] = {NULL}; +TH1F* hPixelMean[MAX_PULS_ORDER] = {NULL}; + +//All Pixel Histograms +TH2F* hAllPixelOverlay[MAX_PULS_ORDER] = {NULL}; +TProfile* hAllPixelProfile[MAX_PULS_ORDER] = {NULL}; +TProfile* hAllPixelProfile2[MAX_PULS_ORDER] = {NULL}; +TH1F* hAllPixelMax[MAX_PULS_ORDER] = {NULL}; +TH1F* hAllPixelMedian[MAX_PULS_ORDER] = {NULL}; +TH1F* hAllPixelMean[MAX_PULS_ORDER] = {NULL}; +TH2F* hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL}; + +//Histogram Parameters +Int_t gPixelOverlayXaxisLeft = 0; +Int_t gPixelOverlayXaxisRight = 0; + +//Root-File Objects +TObjArray* hList = NULL; +TObjArray* hAllPixelList = NULL; +TObjArray* hRootList = NULL; + +//Canvases +TCanvas* cgpPixelPulses[MAX_PULS_ORDER] = {NULL}; +TCanvas* cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL}; +TCanvas* cgpTestHistos = NULL; +//TCanvas* gpcDevelopment = NULL; + +//---------------------------------------------------------------------------- +// Functions +//---------------------------------------------------------------------------- + +void BookHistos(int ); +void BookPixelHistos(int ); +void BookTestHistos( int ); + +void DeletePixelHistos(int ); +void SaveHistograms( const char*, const char*, TObjArray*, int ); + +void FillHistograms(vector*, int, int, int); +void DrawPulseHistograms(int, int); +void DrawMaximumHistograms( int, int ); +void DrawTestHistograms( int); + +void PlotPulseShapeOfMaxPropability(unsigned int, int, int); +void FitMaxPropabilityPuls( TH1F* , int ); + +Double_t MedianOfH1( TH1 *h1 ); +void PlotMedianEachSliceOfPulse(unsigned int, int, int); + +void UpdateCanvases( int, int); +void DeletePixelCanvases( int ); + +void CreateRootFile(const char*, int ); +TFile *OpenRootFile(const char*, int ); +void CloseRootFile(TFile *tf); + +TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file +TString CreateSubDirName(const char* ); + +void WritePixelTemplateToCsv( TString, const char*, const char*, int, int); +void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int); +//void SetHistogrammSettings( const char*, int, int , int); +//void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int); +//void DeletePulseCanvas( TCanvas* , unsigned int); + + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +// MAIN - Funtion +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +int FCalcPulseTemplate( + char* datafilename = "../data/2011/11/09/20111109_006.fits.gz", + const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz", + const char* OutRootFileName = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root", + const char* OutPutPath = "../analysis/FPulseTemplate/20111109_006/Templates/", + bool ProduceGraphic = false, + int refresh_rate = 500, //refresh rate for canvases + bool spikeDebug = false, + bool debugPixel = false, + bool fitdata = false, + int verbosityLevel = 0, // different verbosity levels can be implemented here + int firstevent = 0, + int nevents = -1, + int firstpixel = 0, + int npixel = -1, + int sampleSetSize = 36, + int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms + float GainMean = 8, + float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother + int avg1 = 8, + int avg2 = 8, + int OverlayWindowLeft = 70, + int OverlayWindowRight = 230 + ) +{ + gGainMean = GainMean; + gBSLMean = BSLMean; +//---------------------------------------------------------------------------- +// Save-Root-File Settings +//---------------------------------------------------------------------------- + CreateRootFile( OutRootFileName, verbosityLevel ); + + +//---------------------------------------------------------------------------- +// root global Settings +//---------------------------------------------------------------------------- + + gPixelOverlayXaxisLeft = OverlayWindowLeft; + gPixelOverlayXaxisRight = OverlayWindowRight; + + gStyle->SetPalette(1,0); + gROOT->SetStyle("Plain"); +// gPad->SetGrid(); + + for ( + int pulse_order = MAX_PULS_ORDER - 1; + pulse_order >= 0 ; + pulse_order-- + ) + { + TString cName ="cgpPixelPulses"; + cName += pulse_order; + TString cTitle ="Pulses of Order: "; + cTitle += pulse_order; + cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800); + cgpPixelPulses[pulse_order]->Divide(2, 2); + cName = "cgpAllPixelPulses"; + cName += pulse_order; + cTitle ="AllPixel average of puls shapes of Order: "; + cTitle += pulse_order; + cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 ); + cgpAllPixelPulses[pulse_order]->Divide(2, 2); + } + + // Create (pointer to) Canvases, which are used in every run, + // also in 'non-debug' runs + // Canvases only need if spike Debug, but I want to deklare + // the pointers anyway ... + + TCanvas *cFiltered = NULL; + if (spikeDebug) + { + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300); + cFiltered->Divide(1, 3); + } + + cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 ); + cgpTestHistos->Divide(2,0); +//----------------------------------------------------------------------------- +// Filter-Settings +//----------------------------------------------------------------------------- +// 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.; + +//----------------------------------------------------------------------------- +// prepare datafile +//----------------------------------------------------------------------------- + +// Open the data file + + fits * datafile; + // Opens the raw data file and 'binds' the variables given as + // Parameters to the data file. So they are filled with + // raw data as soon as datafile->GetRow(int) is called. + gNEvents = openDataFits( + datafilename, + &datafile, + AllPixelDataVector, + StartCellVector, + CurrentEventID, + gRegionOfInterest, + NumberOfPixels, + PXLxROI, + verbosityLevel + ); + + if (gNEvents == 0) + { + cout << "return code of OpenDataFile:" << datafilename<< endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + if (verbosityLevel > 0) + { + cout << endl <<"number of events in file: "<< gNEvents << "\t"; + } + + if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< nevents << " will be processed"<< endl; + cout <<"Total # of Pixel: "<< NumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + RC = openCalibFits( + drsfilename, + drs_basemean, + drs_gainmean, + drs_triggeroffsetmean, + TriggerOffsetROI + ); + + if (RC == 0) + { + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + // Book the histograms + BookHistos(verbosityLevel ); + BookTestHistos( verbosityLevel ); +//------------------------------------- +// Loop over Pixel Sets +//------------------------------------- + for ( int firstPixelOfSample = firstpixel; + firstPixelOfSample < firstpixel + npixel; + firstPixelOfSample += sampleSetSize ) + { + + if (verbosityLevel >= 0) + { + cout << "------------------------------------------------" << endl + << "...processing Sample from Pixel: " + << firstPixelOfSample; + << "to Pixel: " + << firstPixelOfSample+sampleSetSize-1 << endl; + } + BookPixelHistos(verbosityLevel ); + +//-------------------------------------------------------------------- +// Loops over Every Event of Pixel +//-------------------------------------------------------------------- + for ( int ev = firstevent; ev < firstevent + nevents; ev++) + { + // Get an Event --> consists of 1440 Pixel ...erm....data + datafile->GetRow( ev ); + + if (verbosityLevel > 0) + { + cout << "-------------------------------------" << endl + << "...processing Event: " << CurrentEventID + << "/" << nevents << " of Pixel: " << pixel << endl; + } + +//-------------------------------------------------------------------- +// Loops over every Pixel of a Set of Pixels +//-------------------------------------------------------------------- + for ( int pixel = firstPixelOfSample; + pixel < sampleSetSize || pixel < firstpixel + npixel; + pixel++ ) + { + if (verbosityLevel > 0) + { + cout << "-------------------------------------" << endl + << "...processing Pixel: " << pixel + << "/" << nPixelSample << endl; + } + +//------------------------------------- +// Apply Calibration +//------------------------------------- + if (verbosityLevel > 2) cout << "...applying DrsCalibration"; + applyDrsCalibration( + Ameas, + pixel, + 12, + 12, + drs_basemean, + drs_gainmean, + drs_triggeroffsetmean, + gRegionOfInterest, + AllPixelDataVector, + StartCellVector + ); + if (verbosityLevel > 2) cout << "...done " << endl; + +//------------------------------------- +// Apply Filters +//------------------------------------- + // finds spikes in the raw data, and interpolates the value + // spikes are: 1 or 2 slice wide, positive non physical artifacts + if (verbosityLevel > 2) cout << "...removeing Spikes"; + removeSpikes (Ameas, Vcorr); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vcorr with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying sliding average filter"; + sliding_avg(Vcorr, Vslide, avg1); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vslide with CFD using FIR filter function + if (verbosityLevel > 2) cout << "...apllying factfir filter"; + factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vcfd with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; + sliding_avg(Vcfd, Vcfd2, avg2); + if (verbosityLevel > 2) cout << "...done " << endl; + +//------------------------------------- +// Search vor Zero crossings +//------------------------------------- + if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ; + // peaks in Ameas[] are found by searching for zero crossings + // in Vcfd2 + // first Argument 1 means ... *rising* edge + // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well + vector* pZXings = zerosearch( Vcfd2 , 1 , 8); + // pZXings means "zero cross ings" + EnlargeRegion(*pZXings, 10, 10); + findAbsMaxInRegions(*pZXings, Vslide); + removeMaximaBelow( *pZXings, 3.0); + removeRegionWithMaxOnEdge( *pZXings, 2); + removeRegionOnFallingEdge( *pZXings, 100); + findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel ); + if (verbosityLevel > 2) cout << "...done" << endl; + +//----------------------------------------------------------------------------- +// Fill Overlay Histos +//----------------------------------------------------------------------------- + FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel ); + +//----------------------------------------------------------------------------- +// Spike Debug +//----------------------------------------------------------------------------- + if ( spikeDebug ) + { + // TODO do this correct. The vectors should be the rigt ones... this is just luck + debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + + for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) + { + debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]); + debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]); + debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); + debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); + debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); + } + + + cFiltered->cd(1); + gPad->SetGrid(); + debugHistos[Ameas_].Draw(); + + cFiltered->cd(2); + gPad->SetGrid(); + debugHistos[Vcorr_].Draw(); + + cFiltered->cd(3); + gPad->SetGrid(); + debugHistos[Vslide_].Draw(); + + TBox *OneBox; + vector MyBoxes; + for (unsigned int i=0; isize(); i++){ + OneBox = new TBox( + pZXings->at(i).maxPos -10 , + pZXings->at(i).maxVal -0.5, + pZXings->at(i).maxPos +10 , + pZXings->at(i).maxVal +0.5); + OneBox->SetLineColor(kBlue); + OneBox->SetLineWidth(1); + OneBox->SetFillStyle(0); + OneBox->SetFillColor(kRed); + MyBoxes.push_back(OneBox); + OneBox->Draw(); + } + +// cFiltered->cd(3); +// gPad->SetGrid(); +// debugHistos[Vcfd2_].Draw(); +// TLine *zeroline = new TLine(0, 0, 1024, 0); +// zeroline->SetLineColor(kBlue); +// zeroline->Draw(); + + cFiltered->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") { + breakout=true; + } + + //TODO!!!!!!!!! + // do some Garbage collection ... + // all the Objects on the heap should be deleted here. + + }// end of if(spikeDebug) + + delete pZXings; + if (breakout) break; +//------------------------------------- +// Draw 1. Set of Pixel Histograms +//------------------------------------- + if ((ev % refresh_rate) == 0) + { + DrawPulseHistograms( + verbosityLevel, + MAX_PULS_ORDER + ); + DrawTestHistograms( + verbosityLevel + ); + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + MAX_PULS_ORDER + ); + } + } + + if (breakout) break; + + if (verbosityLevel > 2) + { + cout << endl << "-------------------------------------"<< endl; + } + + + + + }//End of Loop over Sample +//------------------------------------- +// end of Loops over Sample +//------------------------------------- + + }//End of Loop over Events +//------------------------------------- +// end of Loops over Events +//------------------------------------- + + + + +//------------------------------------- +// Histogramms of Maximas in Overlay Spectra +//------------------------------------- + + PlotPulseShapeOfMaxPropability( + MAX_PULS_ORDER, + fitdata, + verbosityLevel + ); + + for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...plotting most probable pulse of" + << "pulse order " << pulse_order; + PlotMedianEachSliceOfPulse( + hPixelOverlay[pulse_order], + hPixelMedian[pulse_order], + fitdata, + verbosityLevel + ); + + PlotMeanEachSliceOfPulse( + hPixelOverlay[pulse_order], + hPixelMean[pulse_order], + fitdata, + verbosityLevel + ); + } + WritePixelTemplateToCsv( + OutPutPath, + "PulseTemplate_PointSet", + "Maximum", + pixel, + verbosityLevel + ); + + if (verbosityLevel > 2) cout << "...done" << endl; + //here is what happends at the end of each loop over all pixels + //Save Histograms of Pixel into Output rootfile + DrawPulseHistograms( + verbosityLevel, + MAX_PULS_ORDER + ); + + DrawMaximumHistograms( + verbosityLevel, + MAX_PULS_ORDER + ); + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + MAX_PULS_ORDER + ); + } + +// SaveHistograms( +// OutRootFileName, +// CreateSubDirName(pixel), +// hList, +// verbosityLevel +// ); + + if (debugPixel) + { + //Process gui events asynchronously during input + cout << endl; + 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; + } + } + + DeletePixelHistos(verbosityLevel); + + if (verbosityLevel > 2) + { + cout << endl << "...End of Pixel" + << endl << "------------------------------------------------" + << endl; + } + } + // End of Loop over all Pixels + +//------------------------------------- +// Draw Histograms +//------------------------------------- + + SaveHistograms( //save histograms of all pixel into output root file + OutRootFileName, + CreateSubDirName("All"), + hAllPixelList, + verbosityLevel + ); + +// SaveHistograms( //save histograms of generell results into output root file +// OutRootFileName, +// "root", +// hRootList, +// verbosityLevel +// ); + + + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + MAX_PULS_ORDER + ); + } + + WriteAllPixelTemplateToCsv( + OutPutPath, + "PulseTemplate_PointSet", + "Maximum", + verbosityLevel + ); + + DeletePixelCanvases( verbosityLevel ); + return( 0 ); +} +//---------------------------------------------------------------------------- +// end of main function +//----------------------------------------------------------------------------- + + + + +//----------------------------------------------------------------------------- +// Funktions +//----------------------------------------------------------------------------- + +/* +ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel +wie viel speicher braucht das? +Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren +am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute. +ne klasse wäre super + +*/ +void +BookPixelHistos( + int sampleSize, + + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; + hList = new TObjArray; + + TString histoname; + for (int order = 0; order < MAX_PULS_ORDER; order ++) + { +// SetHistogramAttributes( +// "PeakOverlay", +// order, +// gOverlayAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hPixelOverlay"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelOverlay[order] = new TH2F( + histoname, + "Overlay of detected pulses of one pulse order for one Pixel", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 200.5 + ); + hPixelOverlay[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hPixelProfile->SetBit(TH2F::kCanRebin); + hList->Add( hPixelOverlay[order] ); + + //------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakMaximum", +// order, +// gMaximumAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hPixelMax"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelMax[order] = new TH1F( + histoname, + "Maximum value of each slice in overlay plot", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hPixelMax[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hList->Add( hPixelMax[order] ); + +//------------------------------------------------------------------------ + // SetHistogramAttributes( + // "PeakMaximum", + // order, + // gMaximumAttrib, + // MaxAmplOfFirstPulse + // ); + histoname = "hPixelMedian"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelMedian[order] = new TH1F( + histoname, + "Median of each slice in overlay plot", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hPixelMedian[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelMedian[order]->SetLineColor(kRed); + hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hList->Add( hPixelMedian[order] ); + +//------------------------------------------------------------------------ + // SetHistogramAttributes( + // "PeakMaximum", + // order, + // gMaximumAttrib, + // MaxAmplOfFirstPulse + // ); + histoname = "hPixelMean"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelMean[order] = new TH1F( + histoname, + "Mean of each slice in overlay plot", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hPixelMean[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelMean[order]->SetLineColor(kBlue); + hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hList->Add( hPixelMean[order] ); + +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakMaxGaus", +// order, +// gMaxGausAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hPixelEdgeOverlay"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelEdgeOverlay[order] = new TH2F( + histoname, + "Overlay at rising edge of detected pulses of one pulse order for one Pixel ", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 200.5 + ); + + hPixelEdgeOverlay[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hList->Add( hPixelEdgeOverlay[order] ); + + //------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakProfile", +// order, +// gProfileAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hPixelProfile"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelProfile[order] = new TProfile( + histoname, + "Mean value of each slice in overlay plot (Tprofile)", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx + (-1*gPixelOverlayXaxisLeft)-0.5, //xlow + gPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hPixelProfile[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hPixelProfile->SetBit(TH2F::kCanRebin); + hList->Add( hPixelProfile[order] ); + + + //------------------------------------------------------------------------ + // SetHistogramAttributes( + // "PeakProfile", + // order, + // gProfileAttrib, + // MaxAmplOfFirstPulse + // ); + histoname = "hPixelProfile2"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hPixelProfile2[order] = new TProfile( + histoname, + "Mean value of each slice in overlay plot (Tprofile)", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx + (-1*gPixelOverlayXaxisLeft)-0.5, //xlow + gPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hPixelProfile2[order]->SetLineColor(kRed); + hPixelProfile2[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hPixelProfile->SetBit(TH2F::kCanRebin); + + + + + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +//end of BookPixelHistos +//---------------------------------------------------------------------------- + + +void DeletePixelHistos( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ; + for (int order = 0; order < MAX_PULS_ORDER; order ++) + { + if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order; + delete hPixelOverlay[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order; + delete hPixelMax[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order; + delete hPixelMedian[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order; + delete hPixelMean[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order; + delete hPixelEdgeOverlay[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order; + delete hPixelProfile[order]; + if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order; + delete hPixelProfile2[order]; + } + if (verbosityLevel > 3) cout << endl << "...deleting hList"; + delete hList; + if (verbosityLevel > 2) cout << endl << "...done" << endl; +} +//end of DeletePixelHistos +//---------------------------------------------------------------------------- + +void BookTestHistos( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; + + hTesthisto = new TH1F ( + "hTesthisto", + "Deviation of rising edge and maximum", + 600, + -10.1, + 10.1 + ); + hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hTesthisto->GetYaxis()->SetTitle( "counts" ); + + hTesthisto2 = new TH2F ( + "hTesthisto2", + "Deviation of rising edge and maximum by event #", + gNEvents, + 250, + 800, + 600, + -10.1, + 10.1 + ); +// hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); +// hTesthisto2->GetYaxis()->SetTitle( "counts" ); +// hTesthisto2->SetMarkerStyle( 2 ); + + if (verbosityLevel > 2) cout << "...done" << endl; +} +//end of BookTestHistos +//---------------------------------------------------------------------------- + +void BookHistos( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...book histograms" << endl; + hAllPixelList = new TObjArray; + hRootList = new TObjArray; + debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ]; + for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){ + debugHistos[ type ].SetBins(1024, 0, 1024); + debugHistos[ type ].SetLineColor(1); + debugHistos[ type ].SetLineWidth(2); + debugHistos[ type ].SetStats(false); + + // set X axis paras + debugHistos[ type ].GetXaxis()->SetLabelSize(0.08); + debugHistos[ type ].GetXaxis()->SetTitleSize(0.08); + debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0); + debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.)); + + // set Y axis paras + debugHistos[ type ].GetYaxis()->SetLabelSize(0.08); + debugHistos[ type ].GetYaxis()->SetTitleSize(0.08); + debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3); + debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]"); + } + + // All Pixel Histograms + //------------------------------------------------------------------------ + TString histoname; + for (int order = 0; order < MAX_PULS_ORDER; order ++) + { +// SetHistogramAttributes( +// "AllPixelPeakOverlay", +// order, +// gOverlayAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelOverlay"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelOverlay[order] = new TH2F( + histoname, + "Overlay of detected Pulses of one Order for All Pixels", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 300.5 + ); + hAllPixelOverlay[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hAllPixelOverlay->SetBit(TH2F::kCanRebin); + hAllPixelList->Add( hAllPixelOverlay[order] ); + + //------------------------------------------------------------------------ + +// SetHistogramAttributes( +// "AllPixelPeakMax", +// order, +// gMaximumAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelMax"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelMax[order] = new TH1F( + histoname, + "Maximum value of each slice in overlay plot for All Pixels", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hAllPixelMax[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hAllPixelList->Add( hAllPixelMax[order] ); + +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakMaximum", +// order, +// gMaximumAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelMedian"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelMedian[order] = new TH1F( + histoname, + "Median of each slice in overlay plot for All Pixels", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hAllPixelMedian[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelMedian[order]->SetLineColor(kRed); + hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hAllPixelList->Add( hAllPixelMedian[order] ); + if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl; +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakMaximum", +// order, +// gMaximumAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelMean"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelMean[order] = new TH1F( + histoname, + "Mean of each slice in overlay plot for All Pixels", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 + ); + hAllPixelMean[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelMean[order]->SetLineColor(kBlue); + hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hAllPixelList->Add( hAllPixelMean[order] ); + +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "PeakMaxGaus", +// order, +// gMaxGausAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelEdgeOverlay"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelEdgeOverlay[order] = new TH2F( + histoname, + "Overlay at rising edge of detected pulses of one pulse order for All Pixels ", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , + (-1*gPixelOverlayXaxisLeft)-0.5, + gPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 200.5 + ); + + hAllPixelEdgeOverlay[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hAllPixelList->Add( hAllPixelEdgeOverlay[order] ); + +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "AllPixelPeakProfile", +// order, +// gProfileAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelProfile"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelProfile[order] = new TProfile( + histoname, + "Mean value of each slice in overlay plot for All Pixels (Tprofile)", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx + (-1*gPixelOverlayXaxisLeft)-0.5, //xlow + gPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hAllPixelProfile[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hPixelProfile->SetBit(TH2F::kCanRebin); + hAllPixelList->Add( hAllPixelProfile[order] ); + +//------------------------------------------------------------------------ +// SetHistogramAttributes( +// "AllPixelPeakProfile", +// order, +// gProfileAttrib, +// MaxAmplOfFirstPulse +// ); + histoname = "hAllPixelProfile2"; + histoname += order; + if (verbosityLevel > 3) cout << "...booking " << histoname << endl; + hAllPixelProfile2[order] = new TProfile( + histoname, + "Mean value of each slice in overlay plot for All Pixels (Tprofile)", + gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx + (-1*gPixelOverlayXaxisLeft)-0.5, //xlow + gPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hAllPixelProfile2[order]->SetAxisRange( + gBSLMean - 5, + (gGainMean*(order+1)) + 10, + "Y"); + hAllPixelProfile2[order]->SetLineColor(kRed); + hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hPixelProfile->SetBit(TH2F::kCanRebin); + hAllPixelList->Add( hAllPixelProfile2[order] ); + + + + }//end of for over order + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of BookHistos +//---------------------------------------------------------------------------- + + +void CreateRootFile(const char * loc_fname, int verbosityLevel) +{ + TFile* tf = new TFile(loc_fname,"UPDATE"); + if (tf->IsZombie()) + { + cout << "Error opening file" << endl; + exit(-1); + } + else { + if (verbosityLevel > 1) cout << "creating root-file successfull" << endl; + tf->Close(); + } +} +//end of CreateRootFile +//---------------------------------------------------------------------------- + + +TFile* OpenRootFile(const char * loc_fname, int verbosityLevel) +{ + TFile* tf = new TFile(loc_fname,"UPDATE"); + if (tf->IsZombie()) + { + cout << "Error opening file" << endl; + exit(-1); + } + else { + if (verbosityLevel > 1) cout << "...opening root-file:" + << loc_fname + << " successfull" << endl; + } + return tf; +} +//end of OpenRootFile +//---------------------------------------------------------------------------- + + +void CloseRootFile(TFile* tf) +{ + if (tf->IsZombie()) + { + cout << "Error closing file" << endl; + exit(-1); + } else { + tf->Close(); + } +} +//end of CloseRootFile +//---------------------------------------------------------------------------- + + +void SaveHistograms( + const char* loc_fname, + const char* subdirectory, + TObjArray* histList, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...saving pixel histograms to subdirectory: " + << subdirectory << endl ; + + TFile* tf = OpenRootFile(loc_fname, verbosityLevel); + if ( subdirectory != "root") + { + tf->cd(); + tf->mkdir(subdirectory,subdirectory); + tf->cd(subdirectory); + } + else + { + tf->cd(); + } + histList->Write(); // write the major histograms into the top level directory + if (verbosityLevel > 3) tf->ls(); + CloseRootFile( tf ); // close the file + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of SaveHistograms +//---------------------------------------------------------------------------- + + +TString CreateSubDirName(int pixel) +{ + TString path = "Pixel_"; + path += pixel; +// cout << "path " << path << endl ; + return path; +} +// end of CreateSubDirName +//---------------------------------------------------------------------------- + +TString CreateSubDirName(const char* title) +{ + TString path = title; + path += "_Pixel"; +// cout << "path " << path << endl ; + return path; +} +// end of CreateSubDirName +//---------------------------------------------------------------------------- + +void FillHistograms( + vector* pZXings, + int AmplWindowWidth, + int eventNumber, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; + bool use_this_peak=false; + int order_of_pulse=0; + vector::iterator reg; + + //Loop over all found zerocrossings in Region + for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++) + { + //skip those who are at the Rim of the ROI + if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12) + { + if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl; + continue; + } + //domstest->Fill(reg->maxPos); + + // Define axis range of Histogram + int Left = reg->maxPos - gPixelOverlayXaxisLeft; + int Right = reg->maxPos + gPixelOverlayXaxisRight; + int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft; + int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight; + + if (Left < 0) Left =0; + if (Right > (int)Ameas.size() ) Right=Ameas.size() ; + if (EdgeLeft < 0) EdgeLeft =0; + if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ; + + + hTesthisto->Fill( reg->slopeOfRisingEdge ) ; + hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ; + + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + + //determine order of pulse and which histogram shall be filled + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + for(int order = 0; order < MAX_PULS_ORDER; order++) + { + if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2)) + && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2)) + { + //hOverlayTemp = hPixelOverlay[order]; + if (verbosityLevel > 2) cout << "...#" << order ; + order_of_pulse = order; + use_this_peak = true; + break; + } + else if (order >= (MAX_PULS_ORDER - 1)){ + use_this_peak = false; + if (verbosityLevel > 2) cout << "...NONE" << endl ; + } + } + + //Fill overlay und profile histograms + if (use_this_peak){ + if (verbosityLevel > 2) cout << "...filling" ; + for ( int pos = Left; pos < Right; pos++){ +// if (); + hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + } + for ( int pos = EdgeLeft; pos < EdgeRight; pos++){ +// cout << endl << "###filling edge histo###" << endl; + hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; + hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; + hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; + hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; +// cout << endl << "######" << endl; + } + if (verbosityLevel > 2) cout << "...done" << endl; + } + //Breakout option + if (breakout) + break; + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of FillHistograms +//---------------------------------------------------------------------------- + +void DrawTestHistograms( + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; + + cgpTestHistos->cd(1); + hTesthisto->Draw(); + cgpTestHistos->cd(2); + hTesthisto2->Draw("BOX"); +} +// end of DrawTestHistograms +//---------------------------------------------------------------------------- + +void DrawPulseHistograms( + int verbosityLevel, + int max_pulse_order + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; + for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) + { + cgpPixelPulses[pulse_order]->cd(1); + hPixelOverlay[pulse_order]->Draw("COLZ"); + cgpPixelPulses[pulse_order]->cd(3); + hPixelProfile[pulse_order]->Draw("COLZ"); + hPixelProfile2[pulse_order]->Draw("SAME"); + + cgpPixelPulses[pulse_order]->cd(2); + hPixelEdgeOverlay[pulse_order]->Draw("COLZ"); +// cgpPixelPulses[pulse_order]->cd(4); +// hTesthisto->Draw(); + + cgpAllPixelPulses[pulse_order]->cd(1); + hAllPixelOverlay[pulse_order]->Draw("COLZ"); + cgpAllPixelPulses[pulse_order]->cd(2); + hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ"); + + cgpAllPixelPulses[pulse_order]->cd(3); + hAllPixelProfile[pulse_order]->Draw("COLZ"); + hAllPixelProfile2[pulse_order]->Draw("SAME"); + + } +} +// end of DrawPulseHistograms +//---------------------------------------------------------------------------- + + +void DrawMaximumHistograms( + int verbosityLevel, + int max_pulse_order + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ; + for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) + { + cgpPixelPulses[pulse_order]->cd(4); + hPixelMax[pulse_order]->Draw(); + hPixelMedian[pulse_order]->Draw("SAME"); + hPixelMean[pulse_order]->Draw("SAME"); + + cgpAllPixelPulses[pulse_order]->cd(4); + hAllPixelMax[pulse_order]->Draw(); + hAllPixelMedian[pulse_order]->Draw("SAME"); + hAllPixelMean[pulse_order]->Draw("SAME"); +// cgpAllPixelPulses[pulse_order]->cd(3); + // hAllPixelMaxGaus[pulse_order]->Draw("COLZ"); + } +} +// end of DrawMaximumHistograms +//---------------------------------------------------------------------------- + + +void UpdateCanvases( + int verbosityLevel, + int max_pulse_order + ) +{ + if (verbosityLevel > 3) cout << endl << "...updating canvases" ; + for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) + { + cgpPixelPulses[pulse_order]->Modified(); + cgpPixelPulses[pulse_order]->Update(); + cgpAllPixelPulses[pulse_order]->Modified(); + cgpAllPixelPulses[pulse_order]->Update(); + cgpTestHistos->Modified(); + cgpTestHistos->Update(); + } +} +// end of UpdateCanvases +//---------------------------------------------------------------------------- + + + +void +PlotMedianEachSliceOfPulse( + unsigned int max_pulse_order, + int fitdata, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of slice's Median" + << endl; + for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...calculation of " + << "pulse order " << pulse_order; + + Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins(); + + for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) { + TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice); + float median = MedianOfH1(h1); + float mean = h1->GetMean(); + if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median); + delete h1; + hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median ); + hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean ); + hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median ); + hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean ); + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + + if (fitdata) + { + FitMaxPropabilityPuls( + hPixelMedian[pulse_order], + verbosityLevel + ); + } + } +} +// end of PlotMedianEachSliceOfPulse +//---------------------------------------------------------------------------- + +Double_t MedianOfH1 ( + TH1* h1 + ) +{ + //compute the median for 1-d histogram h1 + Int_t nbins = h1->GetXaxis()->GetNbins(); + Double_t *x = new Double_t[nbins]; + Double_t *y = new Double_t[nbins]; + for (Int_t i=0;iGetXaxis()->GetBinCenter(i+1); + y[i] = h1->GetBinContent(i+1); + } + Double_t median = TMath::Median(nbins,x,y); + delete [] x; + delete [] y; + return median; +} +// end of PlotMedianEachSliceOfPulse +//---------------------------------------------------------------------------- + +void +PlotPulseShapeOfMaxPropability( + unsigned int max_pulse_order, + int fitdata, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of max. propability" + << endl; + for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...calculation of " + << "pulse order " << pulse_order; + // vector max_value_of to number of timeslices in Overlay Spectra + float max_value_of_slice; + + Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins(); + if (verbosityLevel > 2) cout << "...generating projections" << endl; + for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice; + //put maximumvalue of every Y-projection of every timeslice into vector + TH1D *h1 = hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice); + + max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() ); + + if (verbosityLevel > 2) cout << " with max value " + << max_value_of_slice << endl; + hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice ); + hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice ); + delete h1; + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + + if (fitdata) + { + FitMaxPropabilityPuls( + hPixelMax[pulse_order], + verbosityLevel + ); + } + } +} + + +void +FitMaxPropabilityPuls( + TH1F* hMaximumTemp, + int verbosityLevel + ) + { + if (verbosityLevel > 2) cout << "...fit Landau in histograms" ; + if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ; + hMaximumTemp->Fit("landau", "", "", -50, 250); + if (verbosityLevel > 2) cout << "...done" << endl; + } + + +//void SetHistogrammSettings( +// const char* histoType, +// int max_puls_order, +// int MaxAmplOfFirstPulse, +// int verbosityLevel +// ) +//{ +// if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl; +// for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++) +// { +// qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order); +// cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl; +// gHistSetting[pulse_order].name = "h"; +// gHistSetting[pulse_order].name += histoType; +// gHistSetting[pulse_order].name += "PeakOverlay"; +// gHistSetting[pulse_order].name += pulse_order; +// gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around "; +// gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax; +// gHistSetting[pulse_order].title += " mV"; +// gHistSetting[pulse_order].Mname = "h"; +// gHistSetting[pulse_order].Mname += histoType; +// gHistSetting[pulse_order].Mname += "PeakMaximum"; +// gHistSetting[pulse_order].Mname += pulse_order; +// gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of "; +// gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax; +// gHistSetting[pulse_order].Mtitle += " mV"; +// } +// if (verbosityLevel > 2) cout << "...done" << endl; +//} + + +//void CreatePulseCanvas( +// TCanvas* array_of_canvases, +// unsigned int max_pulse_order, +// const char* canvas_name, +// const char* canvas_title, +// const char* pixel, +// int verbosityLevel +// ) +//{ +// if (verbosityLevel > 2) cout << "...building Canvases" ; + +// for ( +// unsigned int pulse_order = 1; +// pulse_order <= max_pulse_order; +// pulse_order++ +// ) +// { +// TString name = canvas_name; +// name += " "; +// name += pulse_order; +// name += " "; + +// TString title = "Order "; +// title += pulse_order; +// title += " "; +// title += canvas_title; +// title += " of "; +// title += pixel; +// title += ". Pixel"; + +// array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700); +// array_of_canvases[pulse_order].Divide(2,2); +// } +// if (verbosityLevel > 2) cout << "...done" << endl; +//} + + +//void +//DeletePulseCanvas( +// TCanvas* array_of_canvases, +// unsigned int max_pulse_order +// ) +//{ +// for ( +// unsigned int pulse_order = 1; +// pulse_order <= max_pulse_order; +// pulse_order++ +// ) +// { +// delete array_of_canvases[pulse_order]; +// } +//} + + + + +void +DeletePixelCanvases( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ; + for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++) + { + delete cgpPixelPulses[pulse_order]; + } +} + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- + +//void SetHistogramAttributes( +// const char* histoType, +// int pulse_order, +// histAttrib_t* histoAttrib, +// int max_ampl_of_first_pulse +// ) +//{ +// histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order); +// histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2); +// histoAttrib[pulse_order].name = "h"; +// histoAttrib[pulse_order].name += histoType; +// histoAttrib[pulse_order].name = "_"; +// histoAttrib[pulse_order].name += pulse_order; + +// histoAttrib[pulse_order].title += histoType; +// histoAttrib[pulse_order].title += " with Amplitude around "; +// histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax; +// histoAttrib[pulse_order].title += " mV"; + +// histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]"; +// histoAttrib[pulse_order].yTitle += "Amplitude [mV]"; +//} + +void WritePixelTemplateToCsv( + TString path, + const char* csv_file_name, + const char* overlay_method, + int pixel, + int verbosityLevel + ) +{ + path += csv_file_name; + path += "_"; + path += pixel; + path += ".csv"; + + Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins(); + if (verbosityLevel > 0) + { + cout << "writing point-set to csv file: " ; + cout << path << endl; + cout << "...opening file" << endl; + } + if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl; + ofstream out; + out.open( path ); + out << "### point-set of a single photon pulse template" << endl; + out << "### template determined with pulse overlay at: " + << overlay_method << endl; + out << "### Slice's Amplitude determined by calculating the " << endl + << "### value of maximum propability of slice -> Amplitude1 " << endl + << "### mean of slice -> Amplitude2 " << endl + << "### median of slice -> Amplitude3 " << endl + << "### for each slice" << endl; + out << "### Pixel number (CHid): " << pixel << endl + << endl; + + out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl; + + for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + out << TimeSlice << "," ; + out << hPixelMax[0]->GetBinContent(TimeSlice) << ","; + out << hPixelMean[0]->GetBinContent(TimeSlice) << ","; + out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl; + } + out.close(); + if (verbosityLevel > 0) cout << "...file closed" << endl; +} + +void WriteAllPixelTemplateToCsv( + TString path, + const char* csv_file_name, + const char* overlay_method, + int verbosityLevel + ) +{ + path += csv_file_name; + path += "_AllPixel.csv"; + + Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins(); + if (verbosityLevel > 0) + { + cout << "writing point-set to csv file: " ; + cout << path << endl; + cout << "...opening file" << endl; + } + if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl; + ofstream out; + out.open( path ); + out << "### point-set of a single photon pulse template" << endl; + out << "### template determined with pulse overlay at: " + << overlay_method << "of all Pixels"; + out << "### Slice's Amplitude determined by calculating the " << endl + << "### value of maximum propability of slice -> Amplitude1 " << endl + << "### mean of slice -> Amplitude2 " << endl + << "### median of slice -> Amplitude3 " << endl + << "### for each slice" << endl + << endl; + + out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl; + + for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + out << TimeSlice << "," ; + out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ","; + out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ","; + out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl; + } + out.close(); + if (verbosityLevel > 0) cout << "...file closed" << endl; +} Index: fact/tools/rootmacros/PulseTemplates/FPulseOverlay.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/FPulseOverlay.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/FPulseOverlay.C (revision 13465) @@ -0,0 +1,904 @@ +/********************** FPulseOverlay *********************** +* read FACT raw data +* remove spikes +* calculate baseline +* find peaks (CFD and normal discriminator) ++ search for Peaks in data ++ put peaks into different histos depending on Amplitude ++ draw histos for single, double, tripple, ....Photonpulses ++ draw Parameterdevelopment + so it can be used for detection of other peaks +*************************************************************/ +//---------------------------------------------------------------------------- +// root libraries +//---------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +using namespace std; + +#define NPIX 1440 +#define NCELL 1024 +#define FAD_MAX_SAMPLES 1024 +#define HAVE_ZLIB + +//---------------------------------------------------------------------------- +// rootmacros +//---------------------------------------------------------------------------- +#include "rootfilehandler.h" +#include "rootfilehandler.C" + +#include "../fits.h" + +#include "../openFits.h" +#include "../openFits.c" + +#include "../discriminator.h" +#include "../discriminator.C" +#include "../zerosearch.h" +#include "../zerosearch.C" +#include "../factfir.C" + +#include "../DrsCalibration.C" +#include "../DrsCalibration.h" + +#include "../SpikeRemoval.h" +#include "../SpikeRemoval.C" + +#include "pixel.h" +#include "pixel.C" + +//---------------------------------------------------------------------------- +// Decleration of global variables +//---------------------------------------------------------------------------- + +bool breakout=false; + +vector AllPixelDataVector; +vector StartCellVector; +unsigned int CurrentEventID; +size_t PXLxROI; +UInt_t gRegionOfInterest; +UInt_t NumberOfPixels; +TString histotitle; + +size_t TriggerOffsetROI, RC; +size_t drs_n; +vector drs_basemean; +vector drs_gainmean; +vector drs_triggeroffsetmean; + +vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector Vcorr(FAD_MAX_SAMPLES); // corrected Values +vector Vslide(FAD_MAX_SAMPLES); // sliding average result +vector Vcfd(FAD_MAX_SAMPLES); // CDF result +vector Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average + +int gNEvents=0; +float gGainMean = 9; +float gBSLMean = 0; + +//typedef struct{ +// double maxAmpl; +// int countOfMax; +// } OverlayMaximum; + +//typedef struct{ +// TString name; +// TString title; +// TString xTitle; +// TString yTitle; +// float yMax; +// float yMin; +//} histAttrib_t; + +//histAttrib_t* gHistoAttrib[maxPulseOrder]; +//histAttrib_t* gProfileAttrib[maxPulseOrder]; + + +// histograms +const int Number_Of_Debug_Histo_Types = 7; + +const unsigned int + Ameas_ = 0, + N1mean_ = 1, + Vcorr_ = 2, + Vtest_ = 3, + Vslide_ = 4, + Vcfd_ = 5, + Vcfd2_ = 6; + +//const char* gHistoTypes[8] = { +// "hMaxOverlay", +// "hPixelMax", +// "hEdgeOverlay", +// "hMaxProfile", +// "hAllPixelOverlay", +// "hAllPixelMax", +// "hAllPixelMaxGaus", +// "hAllPixelProfile" +//}; + +//---------------------------------------------------------------------------- +// Initialisation of histograms +//---------------------------------------------------------------------------- + + // Temporary Objects + TH1F* debugHistos = NULL; + //TH2F* hOverlayTemp = NULL; + //TH1D* hProjPeak = NULL; + TH1F* hTesthisto = NULL; + TH2F* hTesthisto2 = NULL; + + //Histogram Parameters + Int_t gPixelOverlayXaxisLeft = 0; + Int_t gPixelOverlayXaxisRight = 0; + + //Root-File Objects + TObjArray* hList[sampleSetSize] = {NULL}; + TObjArray* hRootList = NULL; + + //Canvases + TCanvas* cgpPixelPulses[maxPulseOrder] = {NULL}; + TCanvas* cgpTestHistos = NULL; + //TCanvas* gpcDevelopment = NULL; + +//---------------------------------------------------------------------------- +// Functions +//---------------------------------------------------------------------------- + +void BookDebugHistos(int ); +void BookPixelHistos(int, int, int, int ); +void BookTestHistos( int ); + +void DeletePixelHistos(int ); +void SaveHistograms( const char*, const char*, TObjArray*, int ); + +void FillHistograms(vector*, int, int, int); +void DrawPulseHistograms(int, int); +void DrawTestHistograms( int); + +void UpdateCanvases( int, int); +void DeletePixelCanvases( int ); + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +// MAIN - Funtion +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +int FPulseTemplate( + char* datafilename = "../data/2011/11/09/20111109_006.fits.gz", + const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz", + const char* OutRootFileName = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root", + bool ProduceGraphic = false, + bool spikeDebug = false, + bool debugPixel = false, + int verbosityLevel = 0, // different verbosity levels can be implemented here + int firstevent = 0, + int nevents = -1, + int firstpixel = 0, + int npixel = -1, + int sampleSetSize = 36, + int maxPulseOrder = 1, + int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms + const char* histoOptions = "S", + float GainMean = 8, + float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother + int avg1 = 8, + int avg2 = 8, + int OverlayWindowLeft = 70, + int OverlayWindowRight = 230, + int refresh_rate = 500 //refresh rate for canvases + ) +{ +//---------------------------------------------------------------------------- +// set global variables +//---------------------------------------------------------------------------- + + gGainMean = GainMean; + gBSLMean = BSLMean; + +//---------------------------------------------------------------------------- +// Initialize Pixel +//---------------------------------------------------------------------------- + Pixel* pPixelSet = (Pixel*) calloc( sampleSetSize, sizeof(Pixel) ); + +//---------------------------------------------------------------------------- +// Save-Root-File Settings +//---------------------------------------------------------------------------- + CreateRootFile( OutRootFileName, verbosityLevel ); + + +//---------------------------------------------------------------------------- +// root global Settings +//---------------------------------------------------------------------------- + + gPixelOverlayXaxisLeft = OverlayWindowLeft; + gPixelOverlayXaxisRight = OverlayWindowRight; + + gStyle->SetPalette(1,0); + gROOT->SetStyle("Plain"); +// gPad->SetGrid(); + + int pulsesCanvasFrameNrs[4] = { + 1, // Top left + 2, // Top right + 3, // bottom left + 4 // bootom right + }; + + if (ProduceGraphic) + { + for ( + int pulse_order = maxPulseOrder - 1; + pulse_order >= 0 ; + pulse_order-- + ) + { + TString cName ="cgpPixelPulses"; + cName += pulse_order; + TString cTitle ="Pulses of Order: "; + cTitle += pulse_order; + cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800); + cgpPixelPulses[pulse_order]->Divide(2, 2); + cName = "cgpAllPixelPulses"; + cName += pulse_order; + cTitle ="AllPixel average of puls shapes of Order: "; + cTitle += pulse_order; + cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 ); + cgpAllPixelPulses[pulse_order]->Divide(2, 2); + } + + // Create (pointer to) Canvases, which are used in every run, + // also in 'non-debug' runs + // Canvases only need if spike Debug, but I want to deklare + // the pointers anyway ... + + TCanvas *cFiltered = NULL; + if (spikeDebug) + { + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300); + cFiltered->Divide(1, 3); + } + + cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 ); + cgpTestHistos->Divide(2,0); + } +//----------------------------------------------------------------------------- +// Filter-Settings +//----------------------------------------------------------------------------- +// 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.; + +//----------------------------------------------------------------------------- +// prepare datafile +//----------------------------------------------------------------------------- + +// Open the data file + + fits * datafile; + // Opens the raw data file and 'binds' the variables given as + // Parameters to the data file. So they are filled with + // raw data as soon as datafile->GetRow(int) is called. + gNEvents = openDataFits( + datafilename, + &datafile, + AllPixelDataVector, + StartCellVector, + CurrentEventID, + gRegionOfInterest, + NumberOfPixels, + PXLxROI, + verbosityLevel + ); + + if (gNEvents == 0) + { + cout << "return code of OpenDataFile:" << datafilename<< endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + if (verbosityLevel > 0) + { + cout << endl <<"number of events in file: "<< gNEvents << "\t"; + } + + if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< nevents << " will be processed"<< endl; + cout <<"Total # of Pixel: "<< NumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + RC = openCalibFits( + drsfilename, + drs_basemean, + drs_gainmean, + drs_triggeroffsetmean, + TriggerOffsetROI + ); + + if (RC == 0) + { + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + // Book the histograms + BookDebugHistos(verbosityLevel ); + BookTestHistos( verbosityLevel ); +//------------------------------------- +// Loop over Pixel Sets +//------------------------------------- + for ( int firstPixelOfSample = firstpixel; + firstPixelOfSample < firstpixel + npixel; + firstPixelOfSample += sampleSetSize ) + { + + if (verbosityLevel >= 0) + { + cout << "------------------------------------------------" << endl + << "...processing Sample from Pixel: " + << firstPixelOfSample; + << "to Pixel: " + << firstPixelOfSample+sampleSetSize-1 << endl; + } +// BookPixelHistos( +// sampleSetSize, +// firstPixelOfSample, +// maxPulseOrder, +// verbosityLevel +// ); + +//-------------------------------------------------------------------- +// Loops over Every Event of Pixel +//-------------------------------------------------------------------- + for ( int ev = firstevent; ev < firstevent + nevents; ev++) + { + // Get an Event --> consists of 1440 Pixel ...erm....data + datafile->GetRow( ev ); + + if (verbosityLevel > 0) + { + cout << "-------------------------------------" << endl + << "...processing Event: " << CurrentEventID + << "/" << nevents << " of Pixel: " << pixel << endl; + } + +//-------------------------------------------------------------------- +// Loops over every Pixel of a Set of Pixels +//-------------------------------------------------------------------- + for ( int pixel = firstPixelOfSample; + pixel < sampleSetSize || pixel < firstpixel + npixel; + pixel++ ) + { + if (verbosityLevel > 0) + { + cout << "-------------------------------------" << endl + << "...processing Pixel: " << pixel + << "/" << nPixelSample << endl; + } + + pPixelSet[pixel] = Pixel( + pixel, + maxPulseOrder, + verbosityLevel, + gPixelOverlayXaxisLeft, + gPixelOverlayXaxisRight, + gBSLMean, + gGainMean, + histoOptions + ); + +//------------------------------------- +// Apply Calibration +//------------------------------------- + if (verbosityLevel > 2) cout << "...applying DrsCalibration"; + applyDrsCalibration( + Ameas, + pixel, + 12, + 12, + drs_basemean, + drs_gainmean, + drs_triggeroffsetmean, + gRegionOfInterest, + AllPixelDataVector, + StartCellVector + ); + if (verbosityLevel > 2) cout << "...done " << endl; + +//------------------------------------- +// Apply Filters +//------------------------------------- + // finds spikes in the raw data, and interpolates the value + // spikes are: 1 or 2 slice wide, positive non physical artifacts + if (verbosityLevel > 2) cout << "...removeing Spikes"; + removeSpikes (Ameas, Vcorr); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vcorr with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying sliding average filter"; + sliding_avg(Vcorr, Vslide, avg1); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vslide with CFD using FIR filter function + if (verbosityLevel > 2) cout << "...apllying factfir filter"; + factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter Vcfd with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; + sliding_avg(Vcfd, Vcfd2, avg2); + if (verbosityLevel > 2) cout << "...done " << endl; + +//------------------------------------- +// Search vor Zero crossings +//------------------------------------- + if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ; + // peaks in Ameas[] are found by searching for zero crossings + // in Vcfd2 + // first Argument 1 means ... *rising* edge + // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well + vector* pZXings = zerosearch( Vcfd2 , 1 , 8); + // pZXings means "zero cross ings" + EnlargeRegion(*pZXings, 10, 10); + findAbsMaxInRegions(*pZXings, Vslide); + removeMaximaBelow( *pZXings, 3.0); + removeRegionWithMaxOnEdge( *pZXings, 2); + removeRegionOnFallingEdge( *pZXings, 100); + findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel ); + if (verbosityLevel > 2) cout << "...done" << endl; + +//----------------------------------------------------------------------------- +// Fill Overlay Histos +//----------------------------------------------------------------------------- + FillHistograms(pPixelSet[pixel], pZXings, AmplWindowWidth, ev, verbosityLevel ); + +//----------------------------------------------------------------------------- +// Spike Debug +//----------------------------------------------------------------------------- + if ( spikeDebug ) + { + // TODO do this correct. The vectors should be the rigt ones... this is just luck + debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); + + for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) + { + debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]); + debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]); + debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); + debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); + debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); + } + + + cFiltered->cd(1); + gPad->SetGrid(); + debugHistos[Ameas_].Draw(); + + cFiltered->cd(2); + gPad->SetGrid(); + debugHistos[Vcorr_].Draw(); + + cFiltered->cd(3); + gPad->SetGrid(); + debugHistos[Vslide_].Draw(); + + TBox *OneBox; + vector MyBoxes; + for (unsigned int i=0; isize(); i++){ + OneBox = new TBox( + pZXings->at(i).maxPos -10 , + pZXings->at(i).maxVal -0.5, + pZXings->at(i).maxPos +10 , + pZXings->at(i).maxVal +0.5); + OneBox->SetLineColor(kBlue); + OneBox->SetLineWidth(1); + OneBox->SetFillStyle(0); + OneBox->SetFillColor(kRed); + MyBoxes.push_back(OneBox); + OneBox->Draw(); + } + + // cFiltered->cd(3); + // gPad->SetGrid(); + // debugHistos[Vcfd2_].Draw(); + // TLine *zeroline = new TLine(0, 0, 1024, 0); + // zeroline->SetLineColor(kBlue); + // zeroline->Draw(); + + cFiltered->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") { + breakout=true; + } + + //TODO!!!!!!!!! + // do some Garbage collection ... + // all the Objects on the heap should be deleted here. + + }// end of if(spikeDebug) + + delete pZXings; + if (breakout) break; +//------------------------------------- +// Draw 1. Set of Pixel Histograms +//------------------------------------- + if ((ev % refresh_rate) == 0) + { + if (ProduceGraphic) + { + if (debugPixel) + { + pPixelSet[pixel] ->DrawHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + DrawTestHistograms( + verbosityLevel + ); + + UpdateCanvases( + verbosityLevel, + maxPulseOrder + ); + } + } + } + + if (breakout) break; + + if (verbosityLevel > 2) + { + cout << endl << "...End of Pixel" + << endl << "-------------------------------------"<< endl; + } + }//End of Loop over Sample +//------------------------------------- +// end of Loops over Sample +//------------------------------------- + + }//End of Loop over Events +//------------------------------------- +// end of Loops over Events +//------------------------------------- + + + +//------------------------------------- +// Draw Pixel Histogramms of Overlay Spectra +//------------------------------------- + + for ( int pixel = firstPixelOfSample; + pixel < sampleSetSize || pixel < firstpixel + npixel; + pixel++ ) + { + //here is what happends at the end of each loop over all Events + pPixelSet[pixel]->DrawHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + maxPulseOrder + ); + } + + //Save Histograms of Pixels into Output rootfile + SaveHistograms( + OutRootFileName, + CreateSubDirName( pPixelSet[pixel]->mChid ), + pPixelSet[pixel]->hList, + verbosityLevel + ); + + + if (debugPixel) + { + //Process gui events asynchronously during input + cout << endl; + 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; + } + } + + delete pPixelSet[pixel]; + + if (verbosityLevel > 2) + { + cout << endl << "...End of Pixel-Set" + << endl << "------------------------------------------------" + << endl; + } + } + // End of Loop over all Pixels + +//------------------------------------- +// Draw Histograms +//------------------------------------- + + SaveHistograms( //save histograms of generell results into output root file + OutRootFileName, + "root", + hRootList, + verbosityLevel + ); + + + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + maxPulseOrder + ); + } + + DeletePixelCanvases( verbosityLevel ); + return( 0 ); +} +//---------------------------------------------------------------------------- +// end of main function +//----------------------------------------------------------------------------- + + + + +//----------------------------------------------------------------------------- +// Funktions +//----------------------------------------------------------------------------- + +/* +ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel +wie viel speicher braucht das? +Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren +am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute. +ne klasse wäre super + +*/ + +void BookTestHistos( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; + + hTesthisto = new TH1F ( + "hTesthisto", + "Deviation of rising edge and maximum", + 600, + -10.1, + 10.1 + ); + hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hTesthisto->GetYaxis()->SetTitle( "counts" ); + + hTesthisto2 = new TH2F ( + "hTesthisto2", + "Deviation of rising edge and maximum by event #", + gNEvents, + 250, + 800, + 600, + -10.1, + 10.1 + ); +// hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); +// hTesthisto2->GetYaxis()->SetTitle( "counts" ); +// hTesthisto2->SetMarkerStyle( 2 ); + + if (verbosityLevel > 2) cout << "...done" << endl; +} +//end of BookTestHistos +//---------------------------------------------------------------------------- + + +void BookDebugHistos( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...book histograms" << endl; + hAllPixelList = new TObjArray; + hRootList = new TObjArray; + debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ]; + for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){ + debugHistos[ type ].SetBins(1024, 0, 1024); + debugHistos[ type ].SetLineColor(1); + debugHistos[ type ].SetLineWidth(2); + debugHistos[ type ].SetStats(false); + + // set X axis paras + debugHistos[ type ].GetXaxis()->SetLabelSize(0.08); + debugHistos[ type ].GetXaxis()->SetTitleSize(0.08); + debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0); + debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.)); + + // set Y axis paras + debugHistos[ type ].GetYaxis()->SetLabelSize(0.08); + debugHistos[ type ].GetYaxis()->SetTitleSize(0.08); + debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3); + debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]"); + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of BookDebugHistos +//---------------------------------------------------------------------------- + + +void FillHistograms( + Pixel* CurrentPixel, + vector* pZXings, + int AmplWindowWidth, + int eventNumber, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; + bool use_this_peak=false; + int order_of_pulse=0; + vector::iterator reg; + + //Loop over all found zerocrossings in Region + for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++) + { + //skip those who are at the Rim of the ROI + if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12) + { + if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl; + continue; + } + //domstest->Fill(reg->maxPos); + + // Define axis range of Histogram + int Left = reg->maxPos - gPixelOverlayXaxisLeft; + int Right = reg->maxPos + gPixelOverlayXaxisRight; + int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft; + int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight; + + if (Left < 0) Left =0; + if (Right > (int)Ameas.size() ) Right=Ameas.size() ; + if (EdgeLeft < 0) EdgeLeft =0; + if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ; + + + CurrentPixel->hSlopeRisingEdge->Fill( reg->slopeOfRisingEdge ) ; + hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ; + + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + + //determine order of pulse and which histogram shall be filled + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + for(int order = 0; order < maxPulseOrder; order++) + { + if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2)) + && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2)) + { + //hOverlayTemp = hMaxOverlay[order]; + if (verbosityLevel > 2) cout << "...#" << order ; + order_of_pulse = order; + use_this_peak = true; + break; + } + else if (order >= (maxPulseOrder - 1)){ + use_this_peak = false; + if (verbosityLevel > 2) cout << "...NONE" << endl ; + } + } + + //Fill overlay und profile histograms + if (use_this_peak){ + if (verbosityLevel > 2) cout << "...filling" ; + for ( int pos = Left; pos < Right; pos++){ +// if (); + CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; + } + for ( int pos = EdgeLeft; pos < EdgeRight; pos++){ +// cout << endl << "###filling edge histo###" << endl; + CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; + CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; +// cout << endl << "######" << endl; + } + if (verbosityLevel > 2) cout << "...done" << endl; + } + //Breakout option + if (breakout) + break; + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of FillHistograms +//---------------------------------------------------------------------------- + +void DrawTestHistograms( + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; + + cgpTestHistos->cd(1); + hTesthisto->Draw(); + cgpTestHistos->cd(2); + hTesthisto2->Draw("BOX"); +} +// end of DrawTestHistograms +//---------------------------------------------------------------------------- + +void UpdateCanvases( + int verbosityLevel, + int max_pulse_order + ) +{ + if (verbosityLevel > 3) cout << endl << "...updating canvases" ; + for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) + { + cgpPixelPulses[pulse_order]->Modified(); + cgpPixelPulses[pulse_order]->Update(); + cgpTestHistos->Modified(); + cgpTestHistos->Update(); + } +} +// end of UpdateCanvases +//---------------------------------------------------------------------------- + +void +DeletePixelCanvases( int verbosityLevel ) +{ + if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ; + for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order ++) + { + delete cgpPixelPulses[pulse_order]; + } +} + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- + + Index: fact/tools/rootmacros/PulseTemplates/SaveToCsv.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/SaveToCsv.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/SaveToCsv.C (revision 13465) @@ -0,0 +1,82 @@ +////////////////////////////////////////////////////////////// +//Functions for FPulsetemplate to write template Pulse to csv +////////////////////////////////////////////////////////////// + +#include "SaveToCsv.h" +#include +#include +#include + +TString +BuildPath( + TString path, + const char* csv_file_name, + int pixel + ) +{ + path += csv_file_name; + path += "_"; + + if (pixel = -1) + path += "AllPixel"; + else + path += pixel; + path += ".csv"; + + return path; +} + + +void +WritePixelTemplateToCsv( + TH1* phInputHistogram, + TString path, + const char* csv_file_name, + const char* overlay_method, + int pixel, + int verbosityLevel + ) +{ + path = BuildPath(path, csv_file_name, pixel); //if want to write ALL pixel= -1 + + Int_t nbins = phInputHistogram->GetXaxis()->GetNbins(); + + if (verbosityLevel > 0) + { + cout << "writing point-set to csv file: " ; + cout << path << endl; + cout << "...opening file" << endl; + } + if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl; + + ofstream out; + out.open( path ); + + out << "### point-set of a single photon pulse template" << endl + << "### template determined with pulse overlay at: " + << overlay_method; + if ( pixel = -1 ) + out << "of all Pixels"; + else + out << endl; + out << "### Slice's Amplitude determined by calculating the " << endl + << "### value of maximum propability of slice -> Amplitude1 " << endl + << "### mean of slice -> Amplitude2 " << endl + << "### median of slice -> Amplitude3 " << endl + << "### for each slice" << endl + << "### Pixel number (CHid): " << pixel << endl + << endl + << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl; + + for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + out << TimeSlice + << "," << phInputHistogram->GetBinContent(TimeSlice); +// << "," << hPixelMean[0]->GetBinContent(TimeSlice); +// << "," << hPixelMedian[0]->GetBinContent(TimeSlice) << endl; + } + + out.close(); + if (verbosityLevel > 0) cout << "...file closed" << endl; +} + Index: fact/tools/rootmacros/PulseTemplates/SaveToCsv.h =================================================================== --- fact/tools/rootmacros/PulseTemplates/SaveToCsv.h (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/SaveToCsv.h (revision 13465) @@ -0,0 +1,23 @@ +#ifndef SAVETOCSV_H +#define SAVETOCSV_H + +////////////////////////////////////////////////////////////// +//Functions for FPulsetemplate to write template Pulse to csv +////////////////////////////////////////////////////////////// + +void +WritePixelTemplateToCsv( + TString path, + const char* csv_file_name, + const char* overlay_method, + int pixel, + int verbosityLevel + ); + +TString +BuildPath( + TString path, + const char* csv_file_name, + int pixel + ); +#endif // SAVETOCSV_H Index: fact/tools/rootmacros/PulseTemplates/pixel.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/pixel.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/pixel.C (revision 13465) @@ -0,0 +1,321 @@ +#include "pixel.h" // class implemented +using namespace std; + +/////////////////////////////// PUBLIC /////////////////////////////////////// + +//============================= LIFECYCLE ==================================== + +//Pixel::Pixel(int pixelID) +//{ +// Chid = pixelID; +// max_puslse_order = 1; +// hMaxOverlay = NULL; +// hEdgeOverlay = NULL; +// hMaxProfile = NULL; +// hEdgeProfile = NULL; +// hList = NULL; + +//} + +Pixel::Pixel( + int pixelID = 0, + int maxPulsorder = 1, + int verbosityLevel = 0, + int pixelOverlayXaxisLeft, + int pixelOverlayXaxisRight, + int bSLMean, + int gainMean, + const char* options = "", + ) +{ + mChid = pixelID; + mMaxPulseOrder = maxPulsorder; + mVerbosityLevel = verbosityLevel; + mOptions = options; + mPixelOverlayXaxisLeft = pixelOverlayXaxisLeft; + mPixelOverlayXaxisRight = pixelOverlayXaxisRight; + mBSLMean = bSLMean; + mGainMean = gainMean; + + hMaxOverlay = (TH2F*) calloc(mMaxPulseOrder, sizeof(TH2F)); + hEdgeOverlay = (TH2F*) calloc(mMaxPulseOrder, sizeof(TH2F)); + hMaxProfile = (TProfile*) calloc(mMaxPulseOrder, sizeof(TProfile)); + hEdgeProfile = (TProfile*) calloc(mMaxPulseOrder, sizeof(TProfile)); + + if (mOptions == "S") + hSlopeRisingEdge= (TH1F*) calloc(mMaxPulseOrder, sizeof(TH1F)); + if (mOptions == "R") + hRisingToMaxEdge= (TH1F*) calloc(mMaxPulseOrder, sizeof(TH1F)); + if (mOptions == "P") + hPosOfMax = (TH1F*) calloc(mMaxPulseOrder, sizeof(TH1F)); + + hList = new TObjArray; + + BookPixelHistos(); + BookDistributionHistos() +} + +Pixel::~Pixel() +{ + if (mVerbosityLevel > 2) cout << endl << "...delete current pixel histograms" ; + for (int order = 0; + order < mMaxPulseOrder; + order ++) + { + if (mVerbosityLevel > 3) cout << endl << "...deleting hMaxOverlay" << mHistogramId; + delete hMaxOverlay[order]; + hMaxOverlay = NULL; + if (mVerbosityLevel > 3) cout << endl << "...deleting hEdgeOverlay" << mHistogramId; + delete hEdgeOverlay[order]; + hEdgeOverlay = NULL; + if (mVerbosityLevel > 3) cout << endl << "...deleting hMaxProfile" << mHistogramId; + delete hMaxProfile[order]; + hMaxProfile = NULL; + if (mVerbosityLevel > 3) cout << endl << "...deleting hMaxProfile2" << mHistogramId; + delete hEdgeProfile[order]; + hEdgeProfile = NULL; + } + if (mVerbosityLevel > 3) cout << endl << "...deleting hList"; + delete hList; + if (mVerbosityLevel > 2) cout << endl << "...done" << endl; +}// ~Pixel + + +//============================= OPERATORS ==================================== + +//XX& +//XX::operator=(const XX&); +//{ +// return *this; + +//}// = + +//============================= OPERATIONS =================================== +void +Pixel::BookPixelHistos() +{ + if (mVerbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; + TString histo_name; + for (int order = 0; order < mMaxPulseOrder; order ++) + { +/* + SetHistogramAttributes( + "PeakMaxGaus", + order, + gMaxGausAttrib, + MaxAmplOfFirstPulse + ); + */ + histo_name = "hMaxOverlay"; + histo_name += order; + if (mVerbosityLevel > 3) cout << "...booking " << histo_name << endl; + hMaxOverlay[order] = TH2F( + histo_name, + "Overlay of detected pulses of one pulse order for one Pixel", + mPixelOverlayXaxisLeft + mPixelOverlayXaxisRight , + (-1*mPixelOverlayXaxisLeft)-0.5, + mPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 200.5 + ); + hMaxOverlay[order]->SetAxisRange( + mBSLMean - 5, + (mGainMean*(order+1)) + 10, + "Y"); + hMaxOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hMaxOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hMaxProfile->SetBit(TH2F::kCanRebin); + hList->Add( hMaxOverlay[order] ); + +//------------------------------------------------------------------------ +/* + SetHistogramAttributes( + "PeakMaxGaus", + order, + gMaxGausAttrib, + MaxAmplOfFirstPulse + ); + */ + histo_name = "hEdgeOverlay"; + histo_name += order; + if (mVerbosityLevel > 3) cout << "...booking " << histo_name << endl; + hEdgeOverlay[order] = TH2F( + histo_name, + "Overlay at rising edge of detected pulses of one pulse order for one Pixel ", + mPixelOverlayXaxisLeft + mPixelOverlayXaxisRight , + (-1*mPixelOverlayXaxisLeft)-0.5, + mPixelOverlayXaxisRight-0.5 , + 512, + -55.5, + 200.5 + ); + + hEdgeOverlay[order]->SetAxisRange( + mBSLMean - 5, + (mGainMean*(order+1)) + 10, + "Y"); + hEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + hList->Add( hEdgeOverlay[order] ); + + //------------------------------------------------------------------------ +/* + SetHistogramAttributes( + "PeakMaxGaus", + order, + gMaxGausAttrib, + MaxAmplOfFirstPulse + ); + */ + histo_name = "hMaxProfile"; + histo_name += order; + if (mVerbosityLevel > 3) cout << "...booking " << histo_name << endl; + hMaxProfile[order] = TProfile( + histo_name, + "Mean value of each slice in overlay plot (Tprofile)", + mPixelOverlayXaxisLeft + mPixelOverlayXaxisRight ,//nbinsx + (-1*mPixelOverlayXaxisLeft)-0.5, //xlow + mPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hMaxProfile[order]->SetAxisRange( + mBSLMean - 5, + (mGainMean*(order+1)) + 10, + "Y"); + hMaxProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hMaxProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hMaxProfile->SetBit(TH2F::kCanRebin); + hList->Add( hMaxProfile[order] ); + + + //------------------------------------------------------------------------ +/* + SetHistogramAttributes( + "PeakMaxGaus", + order, + gMaxGausAttrib, + MaxAmplOfFirstPulse + ); + */ + histo_name = "hEdgeProfile"; + histo_name += order; + if (mVerbosityLevel > 3) cout << "...booking " << histo_name << endl; + hEdgeProfile[order] = TProfile( + histo_name, + "Mean value of each slice in overlay plot (Tprofile)", + mPixelOverlayXaxisLeft + mPixelOverlayXaxisRight ,//nbinsx + (-1*mPixelOverlayXaxisLeft)-0.5, //xlow + mPixelOverlayXaxisRight-0.5 , //xup + // 512, //nbinsy + -55.5, //ylow + 300.5, //yup + "s"); //option + hEdgeProfile[order]->SetLineColor(kRed); + hEdgeProfile[order]->SetAxisRange( + mBSLMean - 5, + (mGainMean*(order+1)) + 10, + "Y"); + hEdgeProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hEdgeProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); + //hMaxProfile->SetBit(TH2F::kCanRebin); + } + if (mVerbosityLevel > 2) cout << "...done" << endl; +} +//end of BookPixelHistos +//---------------------------------------------------------------------------- + + +void +Pixel::DrawHistograms( + TCanvas* pixelCanvas = NULL, + int* histoFrameNr = NULL + ) +{ + if (mVerbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; + for (int pulse_order = 0; pulse_order < mMaxPulseOrder; pulse_order++) + { + pixelCanvas[pulse_order]->cd( histoFrameNr[0] ); + hMaxOverlay[pulse_order]->Draw("COLZ"); + pixelCanvas[pulse_order]->cd( histoFrameNr[2] ); + hMaxProfile[pulse_order]->Draw("COLZ"); + hMaxProfile2[pulse_order]->Draw("SAME"); + + pixelCanvas[pulse_order]->cd( histoFrameNr[1] ); + hEdgeOverlay[pulse_order]->Draw("COLZ"); + } +} +// end of DrawPulseHistograms +//---------------------------------------------------------------------------- + + +void +Pixel::BookDistributionHistos( ) +{ + if (mOptions != "") + { + if (mVerbosityLevel > 2) cout << endl + << "...book distribution histograms" + << endl; + + if (mOptions != "S") + { + hSlopeRisingEdge = TH1F ( + "hSlopeRisingEdge", + "Distribution of rising edge", + 600, + -10.1, + 10.1 + ); + hSlopeRisingEdge->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hSlopeRisingEdge->GetYaxis()->SetTitle( "counts" ); + } + + if (mOptions != "R") + { + hRisingToMaxEdge = TH1F ( + "hRisingToMaxEdge", + "Deviation of rising edge and maximum", + 600, + -10.1, + 10.1 + ); + hRisingToMaxEdge->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hRisingToMaxEdge->GetYaxis()->SetTitle( "counts" ); + } + + if (mOptions != "P") + { + hPosOfMax = TH1F ( + "hPosOfMax", + "Deviation of rising edge and maximum", + 600, + -10.1, + 10.1 + ); + hPosOfMax->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); + hPosOfMax->GetYaxis()->SetTitle( "counts" ); + } + + if (mVerbosityLevel > 2) cout << "...done" << endl; + } +} + + + + +//============================= ACESS =================================== +//============================= INQUIRY =================================== +/////////////////////////////// PROTECTED /////////////////////////////////// + +/////////////////////////////// PRIVATE /////////////////////////////////// + + + + + + + + Index: fact/tools/rootmacros/PulseTemplates/pixel.h =================================================================== --- fact/tools/rootmacros/PulseTemplates/pixel.h (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/pixel.h (revision 13465) @@ -0,0 +1,115 @@ +/** A one line description of the class. + * + * #include "XX.h"
+ * -llib + * + * A longer description. + * + * @see something + */ + +#ifndef PIXEL_H +#define PIXEL_H + +// SYSTEM INCLUDES + +#include +#include +#include +#include + +// + +// PROJECT INCLUDES +// + +// LOCAL INCLUDES +// + +// FORWARD REFERENCES +// + + +class Pixel +{ +public: +// LIFECYCLE + + /** Default constructor. + */ +// Pixel(int pixelID); + Pixel( + int pixelID = 0, + int maxPulsorder = 1, + int verbosityLevel = 0, + int pixelOverlayXaxisLeft, + int pixelOverlayXaxisRight, + int bSLMean, + int gainMean, + const char *options = "" + ); + + +// /** Copy constructor. +// * +// * @param from The value to copy to this object. +// */ +// XX(const XX& from); + + + /** Destructor. + */ + ~Pixel(); + + +// OPERATORS + +// /** Assignment operator. +// * +// * @param from THe value to assign to this object. +// * +// * @return A reference to this object. +// */ +// XX& operator=(const XX& from); + +// OPERATIONS + + void BookPixelHistos(); + void BookDistributionHistos(); + + void DrawHistograms( + TCanvas *pixelCanvas, //Array of Canvases + int *histoFrameNR //Array with Canvas Frame numbers + ); + +// ACCESS +// INQUIRY + int mChid; + int mMaxPulseOrder; + int mVerbosityLevel; + int mPixelOverlayXaxisLeft; + int mPixelOverlayXaxisRight; + int mBSLMean; + int mGainMean; + const char* mOptions; + + TH2F* hMaxOverlay; //histogrammm for overlay of detected Peaks + TH2F* hEdgeOverlay; + TProfile* hMaxProfile; //histogrammm for Profile of detected Peaks + TProfile* hEdgeProfile; //histogrammm for Profile of detected Peaks + TH1F* hSlopeRisingEdge; + TH1F* hRisingToMaxEdge; + TH1F* hPosOfMax; + TObjArray* hList; + +protected: +private: +}; + +// INLINE METHODS +// + +// EXTERNAL REFERENCES +// + +#endif // _PIXEL_H_ Index: fact/tools/rootmacros/PulseTemplates/rootfilehandler.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/rootfilehandler.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/rootfilehandler.C (revision 13465) @@ -0,0 +1,165 @@ +#include "rootfilehandler.h" +#include +#include + +void +CreateRootFile( + const char * loc_fname, + int verbosityLevel + ) +{ + TFile* tf = new TFile(loc_fname,"UPDATE"); + if (tf->IsZombie()) + { + cout << "Error opening file" << endl; + exit(-1); + } + else { + if (verbosityLevel > 1) cout << "creating root-file successfull" << endl; + tf->Close(); + } +} +//end of CreateRootFile +//---------------------------------------------------------------------------- + + +TFile* +OpenRootFile( + const char * loc_fname, + int verbosityLevel + ) +{ + TFile* tf = new TFile(loc_fname,"UPDATE"); + if (tf->IsZombie()) + { + cout << "Error opening file" << endl; + exit(-1); + } + else { + if (verbosityLevel > 1) cout << "...opening root-file:" + << loc_fname + << " successfull" << endl; + } + return tf; +} +//end of OpenRootFile +//---------------------------------------------------------------------------- + + +void +CloseRootFile( + TFile* tf + ) +{ + if (tf->IsZombie()) + { + cout << "Error closing file" << endl; + exit(-1); + } else { + tf->Close(); + } +} +//end of CloseRootFile +//---------------------------------------------------------------------------- + + +void +SaveHistograms( + const char* loc_fname, + const char* subdirectory, + TObjArray* histList, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...saving pixel histograms to subdirectory: " + << subdirectory << endl ; + + TFile* tf = OpenRootFile(loc_fname, verbosityLevel); + if ( subdirectory != "root") + { + tf->cd(); + tf->mkdir(subdirectory,subdirectory); + tf->cd(subdirectory); + } + else + { + tf->cd(); + } + histList->Write(); // write the major histograms into the top level directory + if (verbosityLevel > 3) tf->ls(); + CloseRootFile( tf ); // close the file + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of SaveHistograms +//---------------------------------------------------------------------------- + + +TString +CreateSubDirName( + int pixel + ) +{ + TString path = "Pixel_"; + path += pixel; +// cout << "path " << path << endl ; + return path; +} +// end of CreateSubDirName +//---------------------------------------------------------------------------- + + +TString +CreateSubDirName( + const char* title + ) +{ + TString path = title; + path += "_Pixel"; +// cout << "path " << path << endl ; + return path; +} +// end of CreateSubDirName +//---------------------------------------------------------------------------- + + +int +CalculateHistogramId( + int pixelID, + int pulseOrder, + int maxPulseOrder + ) +{ + return (pixelID*maxPulseOrder + pulseOrder); +} +// end of CalculateHistogramId +//---------------------------------------------------------------------------- + +int +CalculatePixelId( + int mHistogramId, + int pulseOrder, + int maxPulseOrder + ) +{ + return (mHistogramId - pulseOrder)/(maxPulseOrder - 1); +} +// end of CalculatePixelIdFromHistogramId +//---------------------------------------------------------------------------- + +CalculatePulseOrder( + int mHistogramId, + int maxPulseOrder + ) +{ + if (mHistogramId < maxPulseOrder) + { + return mHistogramId; + } + else + { + return mHistogramId % (maxPulseOrder - 1); + } +} +// end of CalculatePulseOrderFromHistogramId +//---------------------------------------------------------------------------- Index: fact/tools/rootmacros/PulseTemplates/rootfilehandler.h =================================================================== --- fact/tools/rootmacros/PulseTemplates/rootfilehandler.h (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/rootfilehandler.h (revision 13465) @@ -0,0 +1,60 @@ +#ifndef ROOTFILEHANDLER_H +#define ROOTFILEHANDLER_H + +void +CreateRootFile( + const char* loc_fname, + int verbosityLevel + ); + +TFile* +OpenRootFile( + const char* loc_fname, + int verbosityLevel + ); + +void +CloseRootFile( + TFile* tf + ); + +void +SaveHistograms( + const char* loc_fname, + const char* subdirectory, + TObjArray* histList, + int verbosityLevel + ); + +TString +CreateSubDirName( + int pixel + ); + +TString +CreateSubDirName( + const char* title + ); + +int +CalculateHistogramId( + int pixelID, + int pulseOrder, + int maxPulseOrder + ); + +int +CalculatePixelId( + int histrogramID, + int pulseOrder, + int maxPulseOrder + ); + +int +CalculatePulseOrder( + int mHistogramId, + int maxPulseOrder + ); + + +#endif // ROOTFILEHANDLER_H Index: fact/tools/rootmacros/PulseTemplates/slicecalculation.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/slicecalculation.C (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/slicecalculation.C (revision 13465) @@ -0,0 +1,133 @@ +////////////////////////////////////////////////////////////// +// Functions for FPulsetemplate to calculate +// values for slices of TH2 histograms +////////////////////////////////////////////////////////////// + +#include "slicecalculation.h" +#include +#include +#include +#include + +//compute the median for 1-d histogram h1 +Double_t +MedianOfH1( + TH1* h1, + int verbosityLevel + ) +{ + if (verbosityLevel > 4) cout << endl + << "...calculating median of given TH1" + << endl; + Int_t nbins = h1->GetXaxis()->GetNbins(); + Double_t *x = new Double_t[nbins]; + Double_t *y = new Double_t[nbins]; + + for (Int_t i=0;iGetXaxis()->GetBinCenter(i+1); + y[i] = h1->GetBinContent(i+1); + } + Double_t median = TMath::Median(nbins,x,y); + + delete [] x; + delete [] y; + + if (verbosityLevel > 4) cout << "...done" << endl; + return median; +} + +//compute the median for each slice of a given TH2 histogram +void +PlotMedianEachSliceOfPulse( + TH2* phInputHistogram, + TH1* phOutputHistogram, + int fitdata, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of slice's Median" + << endl; + + Int_t nbins = phInputHistogram->GetXaxis()->GetNbins(); + + for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + TH1 *hProjection = phInputHistogram->ProjectionY("",TimeSlice,TimeSlice); + double median = MedianOfH1(hProjection); + + if (verbosityLevel > 4) + printf("Median of Slice %d, Median=%g\n",TimeSlice,median); + + delete hProjection; + + phOutputHistogram->SetBinContent(TimeSlice, median ); +// hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median ); + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + + if (fitdata) + { + FitMaxPropabilityPuls( + phOutputHistogram, + verbosityLevel + ); + } +} + +//compute the mean for each slice of a given TH2 histogram +void +PlotMeanEachSliceOfPulse( + TH2* phInputHistogram, + TH1* phOutputHistogram, + int fitdata, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of slice's mean" + << endl; + + Int_t nbins = phInputHistogram->GetXaxis()->GetNbins(); + + for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + TH1 *hProjection = phInputHistogram->ProjectionY("",TimeSlice,TimeSlice); + double mean = hProjection->GetMean(); + + if (verbosityLevel > 4) + printf("Mean of Slice %d, mean=%g\n",TimeSlice,median); + + delete hProjection; + + phOutputHistogram->SetBinContent(TimeSlice, mean ); +// hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean ); + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + + if (fitdata) + { + FitMaxPropabilityPuls( + phOutputHistogram, + verbosityLevel + ); + } +} + +//fit a function to the given histogram +void +FitMaxPropabilityPuls( + TH1* phInputHistogram, + float left_border, + float right_border, + int verbosityLevel + ) + { + if (verbosityLevel > 2) cout << "...fit Landau in histograms" ; + if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ; + phInputHistogram->Fit("landau", "", "", left_border, right_border); + if (verbosityLevel > 2) cout << "...done" << endl; + } Index: fact/tools/rootmacros/PulseTemplates/slicecalculation.h =================================================================== --- fact/tools/rootmacros/PulseTemplates/slicecalculation.h (revision 13465) +++ fact/tools/rootmacros/PulseTemplates/slicecalculation.h (revision 13465) @@ -0,0 +1,40 @@ +#ifndef __SLICECALCULATION_H +#define __SLICECALCULATION_H + +////////////////////////////////////////////////////////////// +// Functions for FPulsetemplate to calculate +// values for slices of TH2 histograms +////////////////////////////////////////////////////////////// + +Double_t +MedianOfH1( + TH1* h1, + int verbosityLevel + ); + +void +PlotMedianEachSliceOfPulse( + TH2* phInputHistogram, + TH1* phOutputHistogram, + int fitdata, + int verbosityLevel + ); + +//compute the mean for each slice of a given TH2 histogram +void +PlotMeanEachSliceOfPulse( + TH2* phInputHistogram, + TH1* phOutputHistogram , + int fitdata, + int verbosityLevel + ); + +void +FitMaxPropabilityPuls( + TH1* phInputHistogram, + float left_border=-50, + float right_border=250, + int verbosityLevel + ); + +#endif // SLICECALCULATION_H