Index: fact/tools/rootmacros/PulseTemplates/FPulseOverlay_filelist.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/FPulseOverlay_filelist.C (revision 14782) +++ fact/tools/rootmacros/PulseTemplates/FPulseOverlay_filelist.C (revision 14782) @@ -0,0 +1,1576 @@ +/********************** FPulseOverlay *********************** +* read FACT raw data +* remove spikes +* calculate baseline +* find peaks (CFD and normal discriminator) ++ seagRCh 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 +#include +#include +#include +using namespace std; + +#define NPIX 1440 +#define FAD_MAX_SAMPLES 1024 +#define HAVE_ZLIB + +//---------------------------------------------------------------------------- +// rootmacros +//---------------------------------------------------------------------------- + +//#include "../fits.h" + +#include "../openFits.h" +//#include "../openFits.c" + +#include "../discriminator.h" +//#include "../discriminator.C" +#include "../zerosearch.h" +//#include "../zerosearch.h" +# include "../factfir.h" +//#include "../factfir.C" + +//#include "../DrsCalibration.C" +#include "../DrsCalibration.h" + +#include "../SpikeRemoval.h" +//#include "../SpikeRemoval.C" + +#include "rootfilehandler.h" +//#include "rootfilehandler.C" + +#include "pixel.h" +//#include "pixel.C" + +#include "configfile.h" + + +//---------------------------------------------------------------------------- +// Decleration of global variables +//---------------------------------------------------------------------------- + +bool gBreakout=false; + +vector gAllPixelDataVector; +vector gStartCellVector; +unsigned int gCurrentEventID; +size_t gPXLxROI; +UInt_t gRegionOfInterest; +UInt_t gNumberOfPixels; +TString ghistotitle; + +size_t gTriggerOffsetROI; +size_t gRC; +size_t gDrs_n; +vector gDrsBaseMean; +vector gDrsGainMean; +vector gDrsTriggeroffsetMean; + +vector gAmeas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector gVcorr(FAD_MAX_SAMPLES); // corrected Values +vector gVslide(FAD_MAX_SAMPLES); // sliding average result +vector gVcfd(FAD_MAX_SAMPLES); // CDF result +vector gVcfd2(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 + gAmeas_ = 0, + N1mean_ = 1, + gVcorr_ = 2, + Vtest_ = 3, + gVslide_ = 4, + gVcfd_ = 5, + gVcfd2_ = 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[pixelSetSize] = {NULL}; + TList* hRootList = NULL; + + TCanvas** cgpPixelPulses = NULL; + TCanvas** cgpDistributions = NULL; + TCanvas* cFiltered = NULL; + TCanvas* cgpTestHistos = NULL; + +//---------------------------------------------------------------------------- +// Functions +//---------------------------------------------------------------------------- + +void BookDebugHistos(int ); +void BookTestHistos( int ); + +//void DeletePixelHistos(int ); +//void SaveHistograms( const char*, const char*, TObjArray*, int ); + +void FillHistograms(Pixel*, vector*, int, int, int , int ); +void DrawTestHistograms( int); +bool ProduceDebugHistos( vector *pZXings); +bool UseThisPulse( int, int, int, float, int, int); +void UpdateCanvases( int, int, bool); +void DeletePixelCanvases( int, int ); + +int HandleFitsFile( + fits * datafile, + TString datafilename, + TString drsfilename, + int nevents, + int npixel, + int verbosityLevel + ); + +void ReadSequenzFile( + TString fileName, + vector *fitsFiles, + vector *drsFiles + ); + +TString gDataDirectory; +TString gHomeDirectory; +TString ghostname; + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +// MAIN - Funtion +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +int FPulseOverlay_filelist( + TString inputPath = "raw/2012/03/09/", + TString sequenzFileName = "/home/jbuss/macros/isdc_rootmacros/PulseTemplates/seq", +// TString dataFileSuffix = ".fits.gz", +// int firstRunId = 18, +// int nRuns = 50, +// TString drsfilename = "20120309_012.drs.fits.gz", + TString OutputPath = "analysis/FPulseTemplate/20120309_018/Overlay/", + TString OutRootFileName = "20120309_018.root", + int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!! + int nevents = 100, + int firstpixel = 0, + int npixel = 2, + int pixelSetSize = 150, + int maxPulseOrder = 1, + int verbosityLevel = 6, // different verbosity levels can be implemented here + TString histoOptions = "SRM", + int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms + float GainMean = 9, + 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, + bool ProduceGraphic = true, + bool spikeDebug = false, + bool debugPixel = false, + bool stats = false, + bool saveResults = false, + bool testmode = false, + int refresh_rate = 500 //refresh rate for canvases + ) +{ +//----------------------------------------------------------------------------- +// prepare datafiles +//----------------------------------------------------------------------------- + + inputPath = SetHostsPaths(false, inputPath ); + + // declare the data file + fits * datafile = NULL; + + // 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. + + //vector of fits files to read + vector fitsFileNames; + vector drsFileNames; + + ReadSequenzFile( + sequenzFileName, + &fitsFileNames, + &drsFileNames + ); + + cout << "filelist\n" ; + + //fill vector of fits files to read + cout << "generating file list:" << endl; + for (unsigned int fileNr = 0; fileNr < fitsFileNames.size(); fileNr++) + { + if ( fileNr >= drsFileNames.size() ) + { + cout << "too less drs-files in vector, skipping rest" << endl; + break; + + drsFileNames.resize(fileNr); + fitsFileNames.resize(fileNr); + } + + cout << "...checking raw files"; +// cout << "\t fits:\t" << fitsFileNames.at(fileNr); +// cout << "\n\t\t drs:\t" << drsFileNames.at(fileNr) << endl; + + //filename/path construction + TString fits_filename = inputPath; + fits_filename.Append(fitsFileNames.at(fileNr)); + TString drs_filename = inputPath; + drs_filename.Append(drsFileNames.at(fileNr)); + cout << "\n\t fits:\t" << fits_filename; + cout << "\n\t drs:\t" << drs_filename << endl; + + //check if files are readable + ifstream fitsFile(fits_filename); + ifstream drsFile(drs_filename); + if (fitsFile.good() && drsFile.good()) + { + //fill vector + cout << "...added\n" << endl; + fitsFileNames.at(fileNr) = fits_filename; + drsFileNames.at(fileNr) = drs_filename; + } + else { + cout << "...FAILED:"; + cout << "\t deleting line" +// << fitsFileNames.at(fileNr) +// << " and " +// << drsFileNames.at(fileNr) + << " from list\n" << endl; + + fitsFileNames.erase(fitsFileNames.begin()+fileNr); + drsFileNames.erase(drsFileNames.begin()+fileNr); + fileNr--; +// continue; + } + } + cout << endl << "size of file vectors: fits: " + << fitsFileNames.size() + << " drs: " + << drsFileNames.size() << endl; + +//---------------------------------------------------------------------------- +// Save-Root-File Settings +//---------------------------------------------------------------------------- + + OutputPath = SetHostsPaths(true, OutputPath ); + OutRootFileName.Prepend(OutputPath); +// drsfilename = SetHostsPaths(false, drsfilename ); +// OutRootFileName = SetHostsPaths(true, OutRootFileName ); + cout << endl; +// cout << "drs filename:\t" << drsfilename << endl; +// cout << "data filename:\t" << datafilename << endl; + cout << "out filename:\t" << OutRootFileName << endl; + + + if (saveResults) + { + CreateRootFile( OutRootFileName, true, verbosityLevel ); + } + + + + + +//---------------------------------------------------------------------------- +// root global Settings +//---------------------------------------------------------------------------- + + gGainMean = GainMean; + gBSLMean = BSLMean; + gPixelOverlayXaxisLeft = OverlayWindowLeft; + gPixelOverlayXaxisRight = OverlayWindowRight; + + gStyle->SetPalette(1,0); + gROOT->SetStyle("Plain"); +// gPad->SetGrid(); + +//---------------------------------------------------------------------------- +// root Canvas Settings +//---------------------------------------------------------------------------- + //Canvas Pad numbering + int pulsesCanvasFrameNrs[4] = { + 1, // Top left + 2, // Top right + 3, // bottom left + 4 // bootom right + }; + + //Canvas Pad numbering + int distributionCanvasFrameNrs[4] = { + 1, // Top left + 2, // Top right + 3, // bottom left + 4 // bootom right + }; + + if (ProduceGraphic) + { + + //Canvases + cgpPixelPulses = new TCanvas*[maxPulseOrder]; + cgpDistributions = new TCanvas*[maxPulseOrder]; + + //TCanvas* gpcDevelopment = NULL; + TString cName = ""; + TString cTitle = ""; + + //naming of pulse canvases + for ( + int pulse_order = maxPulseOrder - 1; + pulse_order >= 0 ; + pulse_order-- + ) + { + cName ="cgpDistributions"; + cName += pulse_order; + cTitle ="Distributions of Pulses with Order of: "; + cTitle += pulse_order; + cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720); + cgpDistributions[pulse_order]->Divide(2, 2); + cName ="cgpPixelPulses"; + cName += pulse_order; + cTitle ="Overlays of Pulses with Order of: "; + cTitle += pulse_order; + cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720); + cgpPixelPulses[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 ... + if (spikeDebug) + { + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1080,420,360,360); + cFiltered->Divide(1, 3); + } + + if (testmode) + { + //additional Test histograms + cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 ); + 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.; + +//----------------------------------------------------------------------------- +// Book the histograms +//----------------------------------------------------------------------------- + + + + if (testmode) + { + BookTestHistos( verbosityLevel ); + } +//---------------------------------------------------------------------------- +// Initialize Pixel +//---------------------------------------------------------------------------- + Pixel** pixel = new Pixel*[NPIX]; + + for (int i = 0 ; i < NPIX; i++) + { + pixel[i] = NULL; + } + + +//----------------------------------------------------------------------------- +// Main Cycle +//----------------------------------------------------------------------------- + + if ( pixelSetSize == -1 ) + { + pixelSetSize = firstpixel +npixel; + } + int lastPixelOfSet = 0; +//------------------------------------- +// Loop over Pixel Sets +//------------------------------------- + for ( int firstPixelOfSet = firstpixel; + firstPixelOfSet < firstpixel + npixel; + firstPixelOfSet += pixelSetSize ) + { + lastPixelOfSet = firstPixelOfSet + pixelSetSize-1; + if (lastPixelOfSet > firstpixel + npixel) + { + lastPixelOfSet = firstpixel + npixel; + } + if (lastPixelOfSet > 1439) + { + lastPixelOfSet = 1439; + } + if (verbosityLevel == 0) + { + cout << "------------------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet << endl; + } + + //------------------------------------- + // Loop over Fits Files + //------------------------------------- + + for( unsigned int fileId = 0; fileId < fitsFileNames.size(); fileId++ ) + { + datafile = NULL; +// int return_val = 0; +// return_val = HandleFitsFile( +// datafile, +// datafiles.at(fileId), +// drsfilename, +// nevents, +// npixel, +// verbosityLevel +// ); +// if (return_val) +// { +// cerr << "ERROR" << endl; +// return return_val; +// } +//================================================== + + // 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. + cout << "\nopening " << fitsFileNames.at(fileId) << endl; + + gNEvents = openDataFits( + fitsFileNames.at(fileId), + &datafile, + gAllPixelDataVector, + gStartCellVector, + gCurrentEventID, + gRegionOfInterest, + gNumberOfPixels, + gPXLxROI, + verbosityLevel + ); + + if (gNEvents == 0) + { + cout << "return code of OpenDataFile:" << fitsFileNames.at(fileId)<< 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: "<< gNumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + gRC = openCalibFits( + drsFileNames.at(fileId), + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gTriggerOffsetROI + ); + + if (gRC == 0) + { + cout << "return code of openCalibFits:" << drsFileNames.at(fileId) << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + //================================================================== + +//-------------------------------------------------------------------- +// Loop over every Event of Pixel Set +//-------------------------------------------------------------------- + for ( int ev = firstevent ; ev < firstevent + nevents; ev++) + { + // Get an Event --> consists of 1440 Pixel ...erm....data + datafile->GetRow( ev - 1 ); + + if (verbosityLevel == 1) + { + cout << "-------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet + << "... Event: " << gCurrentEventID + << "/" << nevents << endl; + } + + cout << "\t\t success!" << endl; + //-------------------------------------------------------------------- + // Loops over every Pixel of a Set of Pixels + //-------------------------------------------------------------------- + for ( int pixelID = firstPixelOfSet; + pixelID < lastPixelOfSet + 1 + && pixelID < firstpixel + npixel; + pixelID++ ) + { + if (verbosityLevel > 1) + { + cout << "-------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet + << "... Event: " << gCurrentEventID + << "/" << nevents << endl + << " Pixel: " << pixelID + << "/" << firstpixel + npixel -1 << endl; + } + + //------------------------------------- + // Create individual Pixel + //------------------------------------- + if (ev == firstevent) + { + pixel[pixelID] = new Pixel( + pixelID, + maxPulseOrder, + verbosityLevel, + stats, + histoOptions, + gPixelOverlayXaxisLeft, + gPixelOverlayXaxisRight, + gBSLMean, + gGainMean + ); + } + + //------------------------------------- + // Apply Calibration + //------------------------------------- + if (verbosityLevel > 2) cout << "...applying DrsCalibration"; + applyDrsCalibration( + gAmeas, + pixelID, + 12, + 12, + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gRegionOfInterest, + gAllPixelDataVector, + gStartCellVector + ); + 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 (gAmeas, gVcorr); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVcorr with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying sliding average filter"; + sliding_avg(gVcorr, gVslide, avg1); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVslide with CFD using FIR filter function + if (verbosityLevel > 2) cout << "...apllying factfir filter"; + factfir(b_cfd , a_cfd, k_cfd, gVslide, gVcfd); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVcfd with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; + sliding_avg(gVcfd, gVcfd2, avg2); + if (verbosityLevel > 2) cout << "...done " << endl; + + //------------------------------------- + // SeagRCh vor Zero crossings + //------------------------------------- + if (verbosityLevel > 2) cout << endl << "...seagRChing zero crossings" ; + // peaks in gAmeas[] are found by seagRChing for zero crossings + // in gVcfd2 + // first Argument 1 means ... *rising* edge + // second Argument 1 means ... seagRCh with stepsize 1 ... 10 is okay as well + vector* pZXings = zerosearch( gVcfd2 , 1 , 8); + // pZXings means "zero cross ings" + EnlargeRegion(*pZXings, 10, 10); + findAbsMaxInRegions(*pZXings, gVslide); + removeMaximaBelow( *pZXings, 3.0); + removeRegionWithMaxOnEdge( *pZXings, 2); + removeRegionOnFallingEdge( *pZXings, 100); + findTimeOfHalfMaxLeft(*pZXings, gVslide, gBSLMean, 5, 10, verbosityLevel ); + removeRegionWithToFlatSlope(*pZXings, 0.5); + if (verbosityLevel > 2) cout << "...done" << endl; + + //----------------------------------------------------------------------------- + // Fill Overlay Histos + //----------------------------------------------------------------------------- + if (verbosityLevel > 2) + { + cout << endl << "Filling Histograms of Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + FillHistograms( + pixel[pixelID], + pZXings, + AmplWindowWidth, + ev, + // histoOptions, + maxPulseOrder, + verbosityLevel + ); + //----------------------------------------------------------------------------- + // Spike Debug + //----------------------------------------------------------------------------- + if ( spikeDebug ) + { + BookDebugHistos(verbosityLevel ); + gBreakout = ProduceDebugHistos( pZXings ); + }// end of if(spikeDebug) + + delete pZXings; + if (gBreakout) break; + //------------------------------------- + // Draw 1. Set of Pixel Histograms + //------------------------------------- + if ((ev % refresh_rate) == 0) + { + if (ProduceGraphic) + { + if (debugPixel) + { + // for ( int order = 0; + // order < maxPulseOrder; + // order++ + // ) + // { + pixel[firstPixelOfSet]->DrawOverlayHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + pixel[firstPixelOfSet]->DrawDistributionHistograms( + cgpDistributions, + distributionCanvasFrameNrs + ); + // } + if (testmode) + { + DrawTestHistograms( + verbosityLevel + ); + } + + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + } + } + } + + if (gBreakout) break; + + if (verbosityLevel > 1) + { + cout << endl << "...End of Pixel" + << endl << "-------------------------------------"<< endl; + } + }//End of Loop over Set + //------------------------------------- + // end of Loops over Set + //------------------------------------- + + }//End of Loop over Events + //------------------------------------- + // end of Loops over Events + //------------------------------------- + + }//End of Loop over Files + //------------------------------------- + // end of Loops over Files + //------------------------------------- + + +//------------------------------------- +// Draw Pixel Histogramms of Overlay Spectra +//------------------------------------- + + for ( int pixelID = firstPixelOfSet; + pixelID < firstPixelOfSet + pixelSetSize + && pixelID < firstpixel + npixel; + pixelID++ ) + { + //here is what happends at the end of each loop over all Events + + if (ProduceGraphic) + { + pixel[pixelID]->DrawOverlayHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + pixel[pixelID]->DrawDistributionHistograms( + cgpDistributions, + distributionCanvasFrameNrs + ); + + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + } + + //Save Histograms of Pixels into Output rootfile + if (verbosityLevel > 2) + { + cout << endl << "Saving Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + pixel[pixelID]->SavePixelHistograms( OutRootFileName, saveResults ); + + + 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; + } + } + + if (verbosityLevel > 2) + { + cout << endl << "Deleting Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + delete pixel[pixelID]; + + if (verbosityLevel > 2) + { + cout << endl << "...End of Pixel-Set" + << endl << "------------------------------------------------" + << endl; + } + } + } + // End of Loop over all Pixels + +//------------------------------------- +// Draw Histograms +//------------------------------------- + TString test = "root"; + SaveHistograms( //save histograms of generell results into output root file + OutRootFileName, + test, + hRootList, + saveResults, + verbosityLevel + ); + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + + //Process gui events asynchronously during input + cout << endl; + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + timer.TurnOn(); + Getline("Press to go on"); + timer.TurnOff(); + + DeletePixelCanvases( + maxPulseOrder, + verbosityLevel + ); + + if (testmode) + delete cgpTestHistos; + + if (spikeDebug) + delete cFiltered; + } + + delete datafile; + delete[] pixel; + + 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; + hRootList = new TList(); + 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]"); +// hRootList->Add( &debugHistos[ type ] ); + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of BookDebugHistos +//---------------------------------------------------------------------------- + + +void +FillHistograms( + Pixel* CurrentPixel, + vector* pZXings, + int AmplWindowWidth, + int eventNumber, + int maxPulseOrder, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; + if (verbosityLevel > 3) cout << endl << "...EventNR " << eventNumber ; + 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; + } + + // 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 (EdgeLeft < 0) EdgeLeft = 0; + if (Right > (int)gAmeas.size() ) Right = gAmeas.size(); + if (EdgeRight > (int)gAmeas.size() ) EdgeRight = gAmeas.size(); + + + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + //determine order of pulse and which histogram shall be filled + for ( int order = 0; order < maxPulseOrder; order++) + { + use_this_peak = UseThisPulse( + reg->maxPos, + order, + AmplWindowWidth, + gGainMean, + maxPulseOrder, + verbosityLevel + ); + if ( use_this_peak) + { + order_of_pulse = order; + break; + } + } + //Fill histograms + if (use_this_peak) + { + //Histograms for Distributions + if (CurrentPixel->mOptions.Contains("S") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling Edge Histogram" ; + CurrentPixel->hSlopeRisingEdge[order_of_pulse]->Fill( reg->slopeOfRisingEdge ) ; + } + + if (CurrentPixel->mOptions.Contains("R") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling RisingEdgeToMax Histogram" ; + CurrentPixel->hRisingEdgeToMax[order_of_pulse]->Fill( reg->distanceEdgeToMax ) ; + } + + if (CurrentPixel->mOptions.Contains("M") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling PosOfMax Histogram" ; + CurrentPixel->hPosOfMax[order_of_pulse]->Fill( reg->maxPos ) ; + } + + //Histograms for Maximum Overlay + for ( int pos = Left; pos < Right; pos++) + { + CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), gAmeas[pos]) ; + CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), gAmeas[pos]) ; + } + //Histograms for Edge Overlay + for ( int pos = EdgeLeft; pos < EdgeRight; pos++) + { + CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gAmeas[pos]) ; + CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gAmeas[pos]) ; + } + + if (verbosityLevel > 2) cout << "...done" << endl; + } + + //gBreakout option + if (gBreakout) + { + break; + } + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of FillHistograms +//---------------------------------------------------------------------------- + +bool UseThisPulse( + int maxPos, + int order, + int AmplWindowWidth, + float gainMean, + int maxPulseOrder, + int verbosityLevel + ) +{ + //determine if pulse is of given order +bool use_this_peak = false; + if (gAmeas[maxPos] >= ((gainMean*(order+1)) - (AmplWindowWidth/2)) + && gAmeas[maxPos] < ((gainMean*(order+1)) + AmplWindowWidth/2)) + { + if (verbosityLevel > 3) cout << "...#" << order ; + use_this_peak = true; + } + else + { + if (verbosityLevel > 3) + { + if (order >= (maxPulseOrder - 1)) + { + cout << "...NONE" << endl ; + } + } + use_this_peak = false; + } + +return use_this_peak; +} + + +void +DrawTestHistograms( + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms TODO CANVAS" ; + +// cgpTestHistos->cd(1); +// hTesthisto->Draw(); +// cgpTestHistos->cd(2); +// hTesthisto2->Draw("BOX"); +} +// end of DrawTestHistograms +//---------------------------------------------------------------------------- + +void +UpdateCanvases( + int verbosityLevel, + int max_pulse_order, + bool testmode + ) +{ + 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(); + cgpDistributions[pulse_order]->Modified(); + cgpDistributions[pulse_order]->Update(); + if ( testmode ) + { + cgpTestHistos->Modified(); + cgpTestHistos->Update(); + } + } +} +// end of UpdateCanvases +//---------------------------------------------------------------------------- + +void +DeletePixelCanvases( + int maxPulseOrder, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ; + for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ ) + { + delete cgpPixelPulses[pulse_order]; + delete cgpDistributions[pulse_order]; + } + delete[] cgpPixelPulses; + delete[] cgpDistributions; +} + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- + +bool +ProduceDebugHistos( + vector* pZXings + ) +{ + // TODO do this correct. The vectors should be the rigt ones... this is just luck + debugHistos[gAmeas_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcorr_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVslide_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcfd_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcfd2_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + + for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) + { + debugHistos[gAmeas_].SetBinContent(sl, gAmeas[sl]); + debugHistos[gVcorr_].SetBinContent(sl, gVcorr[sl]); + debugHistos[gVslide_].SetBinContent( sl, gVslide[sl] ); + debugHistos[gVcfd_].SetBinContent( sl, gVcfd[sl] ); + debugHistos[gVcfd2_].SetBinContent( sl, gVcfd2[sl] ); + } + + bool has_negative_slope = false; + cFiltered->cd(1); + gPad->SetGrid(); + debugHistos[gAmeas_].Draw(); + + TF1 *OneEdge; + TLine *OneLine; + TBox *OneBox; + vector MyEdges; + vector MyBoxes; + vector MyLines; + + int left = debugHistos[gAmeas_].GetXaxis()->GetFirst(); + int right = debugHistos[gAmeas_].GetXaxis()->GetLast(); + vector::iterator region; + for (region = pZXings->begin() ; region < pZXings->end() ; region++) + { + if (region->slopeOfRisingEdge < 0) + { + has_negative_slope = true; + } + } + + if (has_negative_slope) + { + for (unsigned int i=0; isize(); i++){ + Region* reg = &pZXings->at(i); + cout << "Slope is: " << reg->slopeOfRisingEdge << endl; + +// xBegin = (reg->halfRisingEdgePos - reg->distanceEdgeToMax); +// yBegin = reg->integRCeptRisingEdge + reg->slopeOfRisingEdge * (xBegin) + + OneEdge = new TF1( + "OneEdge", + "[0]*x+[1]", + left, + right + ); + + OneEdge->SetParameters( + reg->slopeOfRisingEdge, + reg->interceptRisingEdge + ); + + OneEdge->SetLineColor(kRed); + OneEdge->SetLineWidth(1); + MyEdges.push_back(OneEdge); + + OneEdge->Draw("SAME"); +// delete OneEdge; + + 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(); + + OneLine = new TLine( + pZXings->at(i).halfRisingEdgePos, + gAmeas[pZXings->at(i).halfRisingEdgePos], + pZXings->at(i).halfRisingEdgePos, + pZXings->at(i).halfRisingEdgePos * pZXings->at(i).slopeOfRisingEdge + pZXings->at(i).interceptRisingEdge + ); + OneLine->SetLineColor(kBlue); + OneLine->SetLineWidth(3); + MyLines.push_back(OneBox); + OneLine->Draw(); + } + + cFiltered->cd(2); + gPad->SetGrid(); + debugHistos[gVslide_].Draw(); + debugHistos[gVcorr_].Draw("SAME"); + + for (unsigned int i=0; iDraw("SAME"); + } + for (unsigned int i=0; iDraw(); + } + for (unsigned int i=0; iDraw(); + } + + cFiltered->cd(3); + gPad->SetGrid(); + debugHistos[gVcorr_].SetLineColor(kBlue); + debugHistos[gVcorr_].Draw(); + debugHistos[gAmeas_].Draw("SAME"); + debugHistos[gVslide_].SetLineColor(kRed); + debugHistos[gVslide_].Draw("SAME"); + debugHistos[gVcfd_].SetLineColor(kGreen); + debugHistos[gVcfd_].Draw("SAME"); + debugHistos[gVcfd2_].SetLineColor(kViolet); + debugHistos[gVcfd2_].Draw("SAME"); + for (unsigned int i=0; iDraw("SAME"); + } + for (unsigned int i=0; iDraw(); + } + for (unsigned int i=0; iDraw(); + } + + +// cFiltered->cd(3); +// gPad->SetGrid(); +// debugHistos[gVcfd2_].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") { + return true; + } + + MyLines.clear(); + MyBoxes.clear(); + MyEdges.clear(); + +// for (unsigned int i=0; isize(); i++) +// { +// delete MyEdges.at(i) +// delete MyBoxes.at(i); +// } +// delete OneBox; +// delete OneEdge; +// delete OneLine; + +// MyEdges.erase(); +// MyBoxes.erase(); + } + + //TODO!!!!!!!!! + // do some Garbage collection ... + // all the Objects on the heap should be deleted here. +return false; +}// end of if(spikeDebug) + +int HandleFitsFile( + fits * datafile, + TString datafilename, + TString drsfilename, + int nevents, + int npixel, + int verbosityLevel + ) +{ + // 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. + cout << "\nopening " << datafilename << endl; + + gNEvents = openDataFits( + datafilename, + &datafile, + gAllPixelDataVector, + gStartCellVector, + gCurrentEventID, + gRegionOfInterest, + gNumberOfPixels, + gPXLxROI, + 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: "<< gNumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + gRC = openCalibFits( + drsfilename, + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gTriggerOffsetROI + ); + + if (gRC == 0) + { + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + return 0; +} + +void ReadSequenzFile(TString fileName, vector *fitsFiles, vector *drsFiles) +{ + ifstream seqFile; + seqFile.open(fileName,ios::in); + + TString fitsName, drsName; + + if(!seqFile) + { + cerr << "***** no config-file in the user defined directory.\n" + << "Check the path .\n" << endl; + return; + } + + if(!seqFile.good()) { + cerr << "***** not able to read config-file." << endl; + return; + } + + while(!seqFile.eof()) + { + seqFile >> fitsName >> drsName; + + if (!fitsName.IsWhitespace() && !fitsName.IsNull()) + { + fitsFiles->push_back(fitsName); + + if (!drsName.IsWhitespace() && !drsName.IsNull()) + { + drsFiles->push_back(drsName); + } + else drsFiles->push_back("0"); + } + } + + seqFile.close(); + + return; +} + + +//int main(int argc,char *argv[]) +//{ + +// TString test; +// TString gRCFileName; +// TString processType = "overlay"; +// bool gRCFileNameCmdSet = false; +// int verbLevel = 0; // different verbosity levels can be implemented here +// bool verbLevelCmdSet = false; +// bool save = false; +// bool produceGraphic = false; + +// // decode arguments +// if(argc < 2) +// { +// printf("no arguements given, using standard arguments\n"); +// } + +// // set conditions for functions arguments +// for (int i=1;i +#include +#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 FAD_MAX_SAMPLES 1024 +#define HAVE_ZLIB + +//---------------------------------------------------------------------------- +// rootmacros +//---------------------------------------------------------------------------- + +//#include "../fits.h" + +#include "../openFits.h" +//#include "../openFits.c" + +#include "../discriminator.h" +//#include "../discriminator.C" +#include "../zerosearch.h" +//#include "../zerosearch.h" +# include "../factfir.h" +//#include "../factfir.C" + +//#include "../DrsCalibration.C" +#include "../DrsCalibration.h" + +#include "../SpikeRemoval.h" +//#include "../SpikeRemoval.C" + +#include "rootfilehandler.h" +//#include "rootfilehandler.C" + +#include "pixel.h" +//#include "pixel.C" + +#include "configfile.h" + + +//---------------------------------------------------------------------------- +// Decleration of global variables +//---------------------------------------------------------------------------- + +bool gBreakout=false; + +vector gAllPixelDataVector; +vector gStartCellVector; +unsigned int gCurrentEventID; +size_t gPXLxROI; +UInt_t gRegionOfInterest; +UInt_t gNumberOfPixels; +TString ghistotitle; + +size_t gTriggerOffsetROI; +size_t gRC; +size_t gDrs_n; +vector gDrsBaseMean; +vector gDrsGainMean; +vector gDrsTriggeroffsetMean; + +vector gAmeas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude +vector gVcorr(FAD_MAX_SAMPLES); // corrected Values +vector gVslide(FAD_MAX_SAMPLES); // sliding average result +vector gVcfd(FAD_MAX_SAMPLES); // CDF result +vector gVcfd2(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 + gAmeas_ = 0, + N1mean_ = 1, + gVcorr_ = 2, + Vtest_ = 3, + gVslide_ = 4, + gVcfd_ = 5, + gVcfd2_ = 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[pixelSetSize] = {NULL}; + TList* hRootList = NULL; + + TCanvas** cgpPixelPulses = NULL; + TCanvas** cgpDistributions = NULL; + TCanvas* cFiltered = NULL; + TCanvas* cgpTestHistos = NULL; + +//---------------------------------------------------------------------------- +// Functions +//---------------------------------------------------------------------------- + +void BookDebugHistos(int ); +void BookTestHistos( int ); + +//void DeletePixelHistos(int ); +//void SaveHistograms( const char*, const char*, TObjArray*, int ); + +void FillHistograms(Pixel*, vector*, int, int, int , int ); +void DrawTestHistograms( int); +bool ProduceDebugHistos( vector *pZXings); +bool UseThisPulse( int, int, int, float, int, int); +void UpdateCanvases( int, int, bool); +void DeletePixelCanvases( int, int ); + +int HandleFitsFile( + fits * datafile, + TString datafilename, + TString drsfilename, + int nevents, + int npixel, + int verbosityLevel + ); + +TString gDataDirectory; +TString gHomeDirectory; +TString ghostname; + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +// MAIN - Funtion +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- +int FPulseOverlay_multifile( + TString inputPath = "raw/2012/03/09/", + TString dataFileIdBase = "20120309", + TString dataFileSuffix = ".fits.gz", + int firstRunId = 18, + int nRuns = 50, + TString drsfilename = "20120309_012.drs.fits.gz", + TString OutputPath = "analysis/FPulseTemplate/20120309_018/Overlay/", + TString OutRootFileName = "20120309_018.root", + int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!! + int nevents = 100, + int firstpixel = 0, + int npixel = 2, + int pixelSetSize = 150, + int maxPulseOrder = 1, + int verbosityLevel = 6, // different verbosity levels can be implemented here + TString histoOptions = "SRM", + int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms + float GainMean = 9, + 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, + bool ProduceGraphic = false, + bool spikeDebug = false, + bool debugPixel = false, + bool stats = false, + bool saveResults = false, + bool testmode = false, + int refresh_rate = 500 //refresh rate for canvases + ) +{ +//----------------------------------------------------------------------------- +// prepare datafiles +//----------------------------------------------------------------------------- + + // declare the data file + fits * datafile = NULL; + + // 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. + + //vector of fits files to read + vector datafiles; + + inputPath = SetHostsPaths(false, inputPath ); + drsfilename.Prepend(inputPath); + + //fill vector of fits files to read + cout << "generating file list:" << endl; + for (int RunId = firstRunId; RunId < firstRunId + nRuns; RunId++) + { + TString run_id; + + run_id = ""; + run_id += RunId; + while (run_id.Length() <3) + run_id.Prepend("0"); + + //filename/path construction + TString temp_filename = inputPath; + temp_filename.Append(dataFileIdBase); + temp_filename.Append("_"); + temp_filename.Append(run_id); + temp_filename.Append(dataFileSuffix); + + cout << "\t checking" << temp_filename; + + //check if file is readable + ifstream fits_file(temp_filename); + if (fits_file.good()) + { + //fill vector + cout << "\t...added" << endl; + datafiles.push_back(temp_filename); + } + else { + cout << "\t...failed" << endl; + } + + } + +//---------------------------------------------------------------------------- +// Save-Root-File Settings +//---------------------------------------------------------------------------- + + OutputPath = SetHostsPaths(true, OutputPath ); + OutRootFileName.Prepend(OutputPath); +// drsfilename = SetHostsPaths(false, drsfilename ); +// OutRootFileName = SetHostsPaths(true, OutRootFileName ); + cout << endl; + cout << "drs filename:\t" << drsfilename << endl; +// cout << "data filename:\t" << datafilename << endl; + cout << "out filename:\t" << OutRootFileName << endl; + + + if (saveResults) + { + CreateRootFile( OutRootFileName, true, verbosityLevel ); + } + + + + + +//---------------------------------------------------------------------------- +// root global Settings +//---------------------------------------------------------------------------- + + gGainMean = GainMean; + gBSLMean = BSLMean; + gPixelOverlayXaxisLeft = OverlayWindowLeft; + gPixelOverlayXaxisRight = OverlayWindowRight; + + gStyle->SetPalette(1,0); + gROOT->SetStyle("Plain"); +// gPad->SetGrid(); + +//---------------------------------------------------------------------------- +// root Canvas Settings +//---------------------------------------------------------------------------- + //Canvas Pad numbering + int pulsesCanvasFrameNrs[4] = { + 1, // Top left + 2, // Top right + 3, // bottom left + 4 // bootom right + }; + + //Canvas Pad numbering + int distributionCanvasFrameNrs[4] = { + 1, // Top left + 2, // Top right + 3, // bottom left + 4 // bootom right + }; + + if (ProduceGraphic) + { + + //Canvases + cgpPixelPulses = new TCanvas*[maxPulseOrder]; + cgpDistributions = new TCanvas*[maxPulseOrder]; + + //TCanvas* gpcDevelopment = NULL; + TString cName = ""; + TString cTitle = ""; + + //naming of pulse canvases + for ( + int pulse_order = maxPulseOrder - 1; + pulse_order >= 0 ; + pulse_order-- + ) + { + cName ="cgpDistributions"; + cName += pulse_order; + cTitle ="Distributions of Pulses with Order of: "; + cTitle += pulse_order; + cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720); + cgpDistributions[pulse_order]->Divide(2, 2); + cName ="cgpPixelPulses"; + cName += pulse_order; + cTitle ="Overlays of Pulses with Order of: "; + cTitle += pulse_order; + cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720); + cgpPixelPulses[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 ... + if (spikeDebug) + { + cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1080,420,360,360); + cFiltered->Divide(1, 3); + } + + if (testmode) + { + //additional Test histograms + cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 ); + 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.; + +//----------------------------------------------------------------------------- +// Book the histograms +//----------------------------------------------------------------------------- + + + + if (testmode) + { + BookTestHistos( verbosityLevel ); + } +//---------------------------------------------------------------------------- +// Initialize Pixel +//---------------------------------------------------------------------------- + Pixel** pixel = new Pixel*[NPIX]; + + for (int i = 0 ; i < NPIX; i++) + { + pixel[i] = NULL; + } + + +//----------------------------------------------------------------------------- +// Main Cycle +//----------------------------------------------------------------------------- + + if ( pixelSetSize == -1 ) + { + pixelSetSize = firstpixel +npixel; + } + int lastPixelOfSet = 0; +//------------------------------------- +// Loop over Pixel Sets +//------------------------------------- + for ( int firstPixelOfSet = firstpixel; + firstPixelOfSet < firstpixel + npixel; + firstPixelOfSet += pixelSetSize ) + { + lastPixelOfSet = firstPixelOfSet + pixelSetSize-1; + if (lastPixelOfSet > firstpixel + npixel) + { + lastPixelOfSet = firstpixel + npixel; + } + if (lastPixelOfSet > 1439) + { + lastPixelOfSet = 1439; + } + if (verbosityLevel == 0) + { + cout << "------------------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet << endl; + } + + //------------------------------------- + // Loop over Fits Files + //------------------------------------- + + for( unsigned int fileId = 0; fileId < datafiles.size(); fileId++ ) + { + datafile = NULL; +// int return_val = 0; +// return_val = HandleFitsFile( +// datafile, +// datafiles.at(fileId), +// drsfilename, +// nevents, +// npixel, +// verbosityLevel +// ); +// if (return_val) +// { +// cerr << "ERROR" << endl; +// return return_val; +// } +//================================================== + + // 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. + cout << "\nopening " << datafiles.at(fileId) << endl; + + gNEvents = openDataFits( + datafiles.at(fileId), + &datafile, + gAllPixelDataVector, + gStartCellVector, + gCurrentEventID, + gRegionOfInterest, + gNumberOfPixels, + gPXLxROI, + verbosityLevel + ); + + if (gNEvents == 0) + { + cout << "return code of OpenDataFile:" << datafiles.at(fileId)<< 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: "<< gNumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + gRC = openCalibFits( + drsfilename, + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gTriggerOffsetROI + ); + + if (gRC == 0) + { + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + + //================================================================== + +//-------------------------------------------------------------------- +// Loop over every Event of Pixel Set +//-------------------------------------------------------------------- + for ( int ev = firstevent ; ev < firstevent + nevents; ev++) + { + // Get an Event --> consists of 1440 Pixel ...erm....data + datafile->GetRow( ev - 1 ); + + if (verbosityLevel == 1) + { + cout << "-------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet + << "... Event: " << gCurrentEventID + << "/" << nevents << endl; + } + + cout << "\t\t success!" << endl; + //-------------------------------------------------------------------- + // Loops over every Pixel of a Set of Pixels + //-------------------------------------------------------------------- + for ( int pixelID = firstPixelOfSet; + pixelID < lastPixelOfSet + 1 + && pixelID < firstpixel + npixel; + pixelID++ ) + { + if (verbosityLevel > 1) + { + cout << "-------------------------------------" << endl + << "...processing Set from Pixel " + << firstPixelOfSet + << " to Pixel " + << lastPixelOfSet + << "... Event: " << gCurrentEventID + << "/" << nevents << endl + << " Pixel: " << pixelID + << "/" << firstpixel + npixel -1 << endl; + } + + //------------------------------------- + // Create individual Pixel + //------------------------------------- + if (ev == firstevent) + { + pixel[pixelID] = new Pixel( + pixelID, + maxPulseOrder, + verbosityLevel, + stats, + histoOptions, + gPixelOverlayXaxisLeft, + gPixelOverlayXaxisRight, + gBSLMean, + gGainMean + ); + } + + //------------------------------------- + // Apply Calibration + //------------------------------------- + if (verbosityLevel > 2) cout << "...applying DrsCalibration"; + applyDrsCalibration( + gAmeas, + pixelID, + 12, + 12, + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gRegionOfInterest, + gAllPixelDataVector, + gStartCellVector + ); + 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 (gAmeas, gVcorr); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVcorr with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying sliding average filter"; + sliding_avg(gVcorr, gVslide, avg1); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVslide with CFD using FIR filter function + if (verbosityLevel > 2) cout << "...apllying factfir filter"; + factfir(b_cfd , a_cfd, k_cfd, gVslide, gVcfd); + if (verbosityLevel > 2) cout << "...done " << endl; + + // filter gVcfd with sliding average using FIR filter function + if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; + sliding_avg(gVcfd, gVcfd2, avg2); + if (verbosityLevel > 2) cout << "...done " << endl; + + //------------------------------------- + // SeagRCh vor Zero crossings + //------------------------------------- + if (verbosityLevel > 2) cout << endl << "...seagRChing zero crossings" ; + // peaks in gAmeas[] are found by seagRChing for zero crossings + // in gVcfd2 + // first Argument 1 means ... *rising* edge + // second Argument 1 means ... seagRCh with stepsize 1 ... 10 is okay as well + vector* pZXings = zerosearch( gVcfd2 , 1 , 8); + // pZXings means "zero cross ings" + EnlargeRegion(*pZXings, 10, 10); + findAbsMaxInRegions(*pZXings, gVslide); + removeMaximaBelow( *pZXings, 3.0); + removeRegionWithMaxOnEdge( *pZXings, 2); + removeRegionOnFallingEdge( *pZXings, 100); + findTimeOfHalfMaxLeft(*pZXings, gVslide, gBSLMean, 5, 10, verbosityLevel ); + removeRegionWithToFlatSlope(*pZXings, 0.5); + if (verbosityLevel > 2) cout << "...done" << endl; + + //----------------------------------------------------------------------------- + // Fill Overlay Histos + //----------------------------------------------------------------------------- + if (verbosityLevel > 2) + { + cout << endl << "Filling Histograms of Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + FillHistograms( + pixel[pixelID], + pZXings, + AmplWindowWidth, + ev, + // histoOptions, + maxPulseOrder, + verbosityLevel + ); + //----------------------------------------------------------------------------- + // Spike Debug + //----------------------------------------------------------------------------- + if ( spikeDebug ) + { + BookDebugHistos(verbosityLevel ); + gBreakout = ProduceDebugHistos( pZXings ); + }// end of if(spikeDebug) + + delete pZXings; + if (gBreakout) break; + //------------------------------------- + // Draw 1. Set of Pixel Histograms + //------------------------------------- + if ((ev % refresh_rate) == 0) + { + if (ProduceGraphic) + { + if (debugPixel) + { + // for ( int order = 0; + // order < maxPulseOrder; + // order++ + // ) + // { + pixel[firstPixelOfSet]->DrawOverlayHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + pixel[firstPixelOfSet]->DrawDistributionHistograms( + cgpDistributions, + distributionCanvasFrameNrs + ); + // } + if (testmode) + { + DrawTestHistograms( + verbosityLevel + ); + } + + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + } + } + } + + if (gBreakout) break; + + if (verbosityLevel > 1) + { + cout << endl << "...End of Pixel" + << endl << "-------------------------------------"<< endl; + } + }//End of Loop over Set + //------------------------------------- + // end of Loops over Set + //------------------------------------- + + }//End of Loop over Events + //------------------------------------- + // end of Loops over Events + //------------------------------------- + + }//End of Loop over Files + //------------------------------------- + // end of Loops over Files + //------------------------------------- + + +//------------------------------------- +// Draw Pixel Histogramms of Overlay Spectra +//------------------------------------- + + for ( int pixelID = firstPixelOfSet; + pixelID < firstPixelOfSet + pixelSetSize + && pixelID < firstpixel + npixel; + pixelID++ ) + { + //here is what happends at the end of each loop over all Events + + if (ProduceGraphic) + { + pixel[pixelID]->DrawOverlayHistograms( + cgpPixelPulses, + pulsesCanvasFrameNrs + ); + + pixel[pixelID]->DrawDistributionHistograms( + cgpDistributions, + distributionCanvasFrameNrs + ); + + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + } + + //Save Histograms of Pixels into Output rootfile + if (verbosityLevel > 2) + { + cout << endl << "Saving Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + pixel[pixelID]->SavePixelHistograms( OutRootFileName, saveResults ); + + + 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; + } + } + + if (verbosityLevel > 2) + { + cout << endl << "Deleting Pixel# [Set] " << pixelID + << " Pixel# [Got] " << pixel[pixelID]->mChid; + } + delete pixel[pixelID]; + + if (verbosityLevel > 2) + { + cout << endl << "...End of Pixel-Set" + << endl << "------------------------------------------------" + << endl; + } + } + } + // End of Loop over all Pixels + +//------------------------------------- +// Draw Histograms +//------------------------------------- + TString test = "root"; + SaveHistograms( //save histograms of generell results into output root file + OutRootFileName, + test, + hRootList, + saveResults, + verbosityLevel + ); + + if (ProduceGraphic) + { + UpdateCanvases( + verbosityLevel, + maxPulseOrder, + testmode + ); + + //Process gui events asynchronously during input + cout << endl; + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + timer.TurnOn(); + Getline("Press to go on"); + timer.TurnOff(); + + DeletePixelCanvases( + maxPulseOrder, + verbosityLevel + ); + + if (testmode) + delete cgpTestHistos; + + if (spikeDebug) + delete cFiltered; + } + + delete datafile; + delete[] pixel; + + 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; + hRootList = new TList(); + 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]"); +// hRootList->Add( &debugHistos[ type ] ); + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of BookDebugHistos +//---------------------------------------------------------------------------- + + +void +FillHistograms( + Pixel* CurrentPixel, + vector* pZXings, + int AmplWindowWidth, + int eventNumber, + int maxPulseOrder, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; + if (verbosityLevel > 3) cout << endl << "...EventNR " << eventNumber ; + 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; + } + + // 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 (EdgeLeft < 0) EdgeLeft = 0; + if (Right > (int)gAmeas.size() ) Right = gAmeas.size(); + if (EdgeRight > (int)gAmeas.size() ) EdgeRight = gAmeas.size(); + + + if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; + //determine order of pulse and which histogram shall be filled + for ( int order = 0; order < maxPulseOrder; order++) + { + use_this_peak = UseThisPulse( + reg->maxPos, + order, + AmplWindowWidth, + gGainMean, + maxPulseOrder, + verbosityLevel + ); + if ( use_this_peak) + { + order_of_pulse = order; + break; + } + } + //Fill histograms + if (use_this_peak) + { + //Histograms for Distributions + if (CurrentPixel->mOptions.Contains("S") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling Edge Histogram" ; + CurrentPixel->hSlopeRisingEdge[order_of_pulse]->Fill( reg->slopeOfRisingEdge ) ; + } + + if (CurrentPixel->mOptions.Contains("R") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling RisingEdgeToMax Histogram" ; + CurrentPixel->hRisingEdgeToMax[order_of_pulse]->Fill( reg->distanceEdgeToMax ) ; + } + + if (CurrentPixel->mOptions.Contains("M") ) + { + if (verbosityLevel > 2) + cout << endl << "\t...filling PosOfMax Histogram" ; + CurrentPixel->hPosOfMax[order_of_pulse]->Fill( reg->maxPos ) ; + } + + //Histograms for Maximum Overlay + for ( int pos = Left; pos < Right; pos++) + { + CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), gAmeas[pos]) ; + CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), gAmeas[pos]) ; + } + //Histograms for Edge Overlay + for ( int pos = EdgeLeft; pos < EdgeRight; pos++) + { + CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gAmeas[pos]) ; + CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gAmeas[pos]) ; + } + + if (verbosityLevel > 2) cout << "...done" << endl; + } + + //gBreakout option + if (gBreakout) + { + break; + } + } + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of FillHistograms +//---------------------------------------------------------------------------- + +bool UseThisPulse( + int maxPos, + int order, + int AmplWindowWidth, + float gainMean, + int maxPulseOrder, + int verbosityLevel + ) +{ + //determine if pulse is of given order +bool use_this_peak = false; + if (gAmeas[maxPos] >= ((gainMean*(order+1)) - (AmplWindowWidth/2)) + && gAmeas[maxPos] < ((gainMean*(order+1)) + AmplWindowWidth/2)) + { + if (verbosityLevel > 3) cout << "...#" << order ; + use_this_peak = true; + } + else + { + if (verbosityLevel > 3) + { + if (order >= (maxPulseOrder - 1)) + { + cout << "...NONE" << endl ; + } + } + use_this_peak = false; + } + +return use_this_peak; +} + + +void +DrawTestHistograms( + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms TODO CANVAS" ; + +// cgpTestHistos->cd(1); +// hTesthisto->Draw(); +// cgpTestHistos->cd(2); +// hTesthisto2->Draw("BOX"); +} +// end of DrawTestHistograms +//---------------------------------------------------------------------------- + +void +UpdateCanvases( + int verbosityLevel, + int max_pulse_order, + bool testmode + ) +{ + 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(); + cgpDistributions[pulse_order]->Modified(); + cgpDistributions[pulse_order]->Update(); + if ( testmode ) + { + cgpTestHistos->Modified(); + cgpTestHistos->Update(); + } + } +} +// end of UpdateCanvases +//---------------------------------------------------------------------------- + +void +DeletePixelCanvases( + int maxPulseOrder, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ; + for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ ) + { + delete cgpPixelPulses[pulse_order]; + delete cgpDistributions[pulse_order]; + } + delete[] cgpPixelPulses; + delete[] cgpDistributions; +} + +//---------------------------------------------------------------------------- +//---------------------------------------------------------------------------- + +bool +ProduceDebugHistos( + vector* pZXings + ) +{ + // TODO do this correct. The vectors should be the rigt ones... this is just luck + debugHistos[gAmeas_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcorr_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVslide_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcfd_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + debugHistos[gVcfd2_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5); + + for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) + { + debugHistos[gAmeas_].SetBinContent(sl, gAmeas[sl]); + debugHistos[gVcorr_].SetBinContent(sl, gVcorr[sl]); + debugHistos[gVslide_].SetBinContent( sl, gVslide[sl] ); + debugHistos[gVcfd_].SetBinContent( sl, gVcfd[sl] ); + debugHistos[gVcfd2_].SetBinContent( sl, gVcfd2[sl] ); + } + + bool has_negative_slope = false; + cFiltered->cd(1); + gPad->SetGrid(); + debugHistos[gAmeas_].Draw(); + + TF1 *OneEdge; + TLine *OneLine; + TBox *OneBox; + vector MyEdges; + vector MyBoxes; + vector MyLines; + + int left = debugHistos[gAmeas_].GetXaxis()->GetFirst(); + int right = debugHistos[gAmeas_].GetXaxis()->GetLast(); + vector::iterator region; + for (region = pZXings->begin() ; region < pZXings->end() ; region++) + { + if (region->slopeOfRisingEdge < 0) + { + has_negative_slope = true; + } + } + + if (has_negative_slope) + { + for (unsigned int i=0; isize(); i++){ + Region* reg = &pZXings->at(i); + cout << "Slope is: " << reg->slopeOfRisingEdge << endl; + +// xBegin = (reg->halfRisingEdgePos - reg->distanceEdgeToMax); +// yBegin = reg->integRCeptRisingEdge + reg->slopeOfRisingEdge * (xBegin) + + OneEdge = new TF1( + "OneEdge", + "[0]*x+[1]", + left, + right + ); + + OneEdge->SetParameters( + reg->slopeOfRisingEdge, + reg->interceptRisingEdge + ); + + OneEdge->SetLineColor(kRed); + OneEdge->SetLineWidth(1); + MyEdges.push_back(OneEdge); + + OneEdge->Draw("SAME"); +// delete OneEdge; + + 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(); + + OneLine = new TLine( + pZXings->at(i).halfRisingEdgePos, + gAmeas[pZXings->at(i).halfRisingEdgePos], + pZXings->at(i).halfRisingEdgePos, + pZXings->at(i).halfRisingEdgePos * pZXings->at(i).slopeOfRisingEdge + pZXings->at(i).interceptRisingEdge + ); + OneLine->SetLineColor(kBlue); + OneLine->SetLineWidth(3); + MyLines.push_back(OneBox); + OneLine->Draw(); + } + + cFiltered->cd(2); + gPad->SetGrid(); + debugHistos[gVslide_].Draw(); + debugHistos[gVcorr_].Draw("SAME"); + + for (unsigned int i=0; iDraw("SAME"); + } + for (unsigned int i=0; iDraw(); + } + for (unsigned int i=0; iDraw(); + } + + cFiltered->cd(3); + gPad->SetGrid(); + debugHistos[gVcorr_].SetLineColor(kBlue); + debugHistos[gVcorr_].Draw(); + debugHistos[gAmeas_].Draw("SAME"); + debugHistos[gVslide_].SetLineColor(kRed); + debugHistos[gVslide_].Draw("SAME"); + debugHistos[gVcfd_].SetLineColor(kGreen); + debugHistos[gVcfd_].Draw("SAME"); + debugHistos[gVcfd2_].SetLineColor(kViolet); + debugHistos[gVcfd2_].Draw("SAME"); + for (unsigned int i=0; iDraw("SAME"); + } + for (unsigned int i=0; iDraw(); + } + for (unsigned int i=0; iDraw(); + } + + +// cFiltered->cd(3); +// gPad->SetGrid(); +// debugHistos[gVcfd2_].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") { + return true; + } + + MyLines.clear(); + MyBoxes.clear(); + MyEdges.clear(); + +// for (unsigned int i=0; isize(); i++) +// { +// delete MyEdges.at(i) +// delete MyBoxes.at(i); +// } +// delete OneBox; +// delete OneEdge; +// delete OneLine; + +// MyEdges.erase(); +// MyBoxes.erase(); + } + + //TODO!!!!!!!!! + // do some Garbage collection ... + // all the Objects on the heap should be deleted here. +return false; +}// end of if(spikeDebug) + +int HandleFitsFile( + fits * datafile, + TString datafilename, + TString drsfilename, + int nevents, + int npixel, + int verbosityLevel + ) +{ + // 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. + cout << "\nopening " << datafilename << endl; + + gNEvents = openDataFits( + datafilename, + &datafile, + gAllPixelDataVector, + gStartCellVector, + gCurrentEventID, + gRegionOfInterest, + gNumberOfPixels, + gPXLxROI, + 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: "<< gNumberOfPixels << "\t"; + } + + if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all! + + if (verbosityLevel > 0) + { + cout <<"of, which "<< npixel << " will be processed"<< endl; + } + + //Get the DRS calibration + gRC = openCalibFits( + drsfilename, + gDrsBaseMean, + gDrsGainMean, + gDrsTriggeroffsetMean, + gTriggerOffsetROI + ); + + if (gRC == 0) + { + cout << "return code of openCalibFits:" << drsfilename << endl; + cout << "is zero -> aborting." << endl; + return 1; + } + return 0; +} + +//int main(int argc,char *argv[]) +//{ + +// TString test; +// TString gRCFileName; +// TString processType = "overlay"; +// bool gRCFileNameCmdSet = false; +// int verbLevel = 0; // different verbosity levels can be implemented here +// bool verbLevelCmdSet = false; +// bool save = false; +// bool produceGraphic = false; + +// // decode arguments +// if(argc < 2) +// { +// printf("no arguements given, using standard arguments\n"); +// } + +// // set conditions for functions arguments +// for (int i=1;iHostName(); +if ( hostname.Contains("isdc") ) +{ + gSystem->Load("/usr/lib64/libz.so"); +} +gROOT->ProcessLine(".L ../openFits.c+"); +gROOT->ProcessLine(".L ../discriminator.C+"); +gROOT->ProcessLine(".L ../zerosearch.C+"); +gROOT->ProcessLine(".L ../factfir.C+"); +gROOT->ProcessLine(".L ../DrsCalibration.C+"); +gROOT->ProcessLine(".L ../SpikeRemoval.C+"); +gROOT->ProcessLine(".L rootfilehandler.C+"); +gROOT->ProcessLine(".L pixel.C+"); +gROOT->ProcessLine(".L FPulseOverlay_filelist.C+"); +} Index: fact/tools/rootmacros/PulseTemplates/make_pulseoverlay_multi.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/make_pulseoverlay_multi.C (revision 14782) +++ fact/tools/rootmacros/PulseTemplates/make_pulseoverlay_multi.C (revision 14782) @@ -0,0 +1,16 @@ +{ +TString hostname = gSystem->HostName(); +if ( hostname.Contains("isdc") ) +{ + gSystem->Load("/usr/lib64/libz.so"); +} +gROOT->ProcessLine(".L ../openFits.c+"); +gROOT->ProcessLine(".L ../discriminator.C+"); +gROOT->ProcessLine(".L ../zerosearch.C+"); +gROOT->ProcessLine(".L ../factfir.C+"); +gROOT->ProcessLine(".L ../DrsCalibration.C+"); +gROOT->ProcessLine(".L ../SpikeRemoval.C+"); +gROOT->ProcessLine(".L rootfilehandler.C+"); +gROOT->ProcessLine(".L pixel.C+"); +gROOT->ProcessLine(".L FPulseOverlay_multifile.C+"); +}