/********************** 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 , bool , int verbosityLevel ); void DrawTestHistograms( int); bool ProduceDebugHistos( vector *pZXings); bool UseThisPulse( int, int, int, float, float bslMean, int maxPulseOrder, int verbosityLevel); 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( //---------------{Process Files}--------------- TString inputPath = "raw/2012/08/02/", TString sequenzFileName = "/home_nfs/isdc/jbbuss/analysis/pulseShape/lists/20120802templ_vs_100mV.csv", TString OutputPath = "analysis/PulsTemplateAnalysis/20120802/", TString OutRootFileName = "20120802_100mV_Stacking_G12B-2AWW8.root", //---------------{Process parameters}--------------- int nRuns = -1, int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!! int nevents = -1, int firstpixel = 0, int npixel = -1, int pixelSetSize = 150, int maxPulseOrder = 6, //---------------{stacking parameters}--------------- TString histoOptions = "SRMTILDAP", int AmplWindowWidth = 8, //Width of Window for selection of pulses to histograms float GainMean = 12, float BSLMean = -2, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother int avg1 = 8, int avg2 = 8, int OverlayWindowLeft = 100, int OverlayWindowRight = 1024, // pulses within this distance behind pulse max will be discarted int beginRisingEdge = 5, int endRisingEdge = 10, // pulses within this distance behind pulse max will be discarted float risingEdgeMinSlope = 0.5, //---------------{graphical options}--------------- bool select_pulses = false, bool ProduceGraphic = false, bool spikeDebug = false, bool debugPixel = false, bool stats = true, bool saveResults = true, bool testmode = false, //---------------{verbosity}--------------- int refresh_rate = 500, //refresh rate for canvases int verbosityLevel = 0 // different verbosity levels can be implemented here ) { //----------------------------------------------------------------------------- // configuration //----------------------------------------------------------------------------- // Amplitude below which pulses will be discarted float minPulseAmplitude = 0.3*(GainMean-BSLMean); // Amplitude above which pulses will be discarted float maxPulseAmplitude = 2*(GainMean-BSLMean)*maxPulseOrder; // pulses within this distance to timeline edge be discarted int timeLineEdgeWidth = 2; // pulses within this distance behind pulse max will be discarted int FallingEdgeWidth = 200; //----------------------------------------------------------------------------- // prepare datafiles //----------------------------------------------------------------------------- inputPath = SetHostsPaths(false, inputPath ); TString fitsFileSuffix = ".fits.gz"; TString drsFileSuffix = ".drs.fits.gz"; // 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; cout << "...reading sequnezfile" << endl; ReadSequenzFile( sequenzFileName, &fitsFileNames, &drsFileNames ); cout << "filelist\n" ; //fill vector of fits files to read cout << "generating file list:" << endl; if(nRuns < fitsFileNames.size() && nRuns != -1) { fitsFileNames.resize( nRuns); drsFileNames.resize( nRuns); } 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)); fits_filename.Append(fitsFileSuffix); TString drs_filename = inputPath; drs_filename.Append(drsFileNames.at(fileNr)); drs_filename.Append(drsFileSuffix); 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; if (fitsFileNames.size() < 1) return 1; //---------------------------------------------------------------------------- // 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 //----------------------------------------------------------------------------- cout << "------------------Main cylcle------------------------------" << endl; if(npixel == -1) npixel = 1439; if ( pixelSetSize == -1 ) { pixelSetSize = firstpixel +npixel; } int lastPixelOfSet = 0; cout << "check " << npixel << endl; //------------------------------------- // 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 //------------------------------------- cout << "loop over files" << endl; 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; } //-------------------------------------------------------------------- // 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 && fileId == 0) { 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; //------------------------------------- // Search 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" int init_reg_size = pZXings->size(); int curr_reg_size = initial_region_size; int step = 0; pixel[pixelID]->hDiscartedPulses->Fill(step, init_reg_size - curr_reg_size); //enlare the search window curr_reg_size = EnlargeRegion(*pZXings, 10, 10); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); // find the the max position of a peak curr_reg_size = findAbsMaxInRegions(*pZXings, gVslide); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); //do a small after pulse analysis, fill in delayed pulses // +++++++++++++++++++++++++++++++++++++++++++++++++++++++ vector::iterator reg_it = pZXings->begin(); int first_pos = pZXings->begin()->maxPos; while( reg_it != pZXings->end() ) { float minimum = 0.0; for (int sl = reg_it->maxPos; sl < 1; i--) { if (gVcfd2[sl] > 0 && gVcfd2[sl-2] <= 0 && gVcfd2[sl-4] < 0) { minimum = gVcfd2[sl-2]; break; } } pixel[pixelID]->hAfterPulses->Fill( reg_it->maxPos - first_pos, reg_it->maxVal - minimum ); } // +++++++++++++++++++++++++++++++++++++++++++++++++++++++ // remove peaks below 1 p.e. curr_reg_size = removeMaximaBelow( *pZXings, minPulseAmplitude); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); // remove peaks above that are larger tha the largest amplitude to be studied curr_reg_size = removeMaximaAbove( *pZXings, maxPulseAmplitude); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); // remove pulses at the pipeline fringes curr_reg_size = removeRegionWithMaxOnEdge( *pZXings, timeLineEdgeWidth); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); // remove "Afterpulses" from pulse list curr_reg_size = removeRegionOnFallingEdge( *pZXings, FallingEdgeWidth); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); // Calculate the position (sample) of a pulses half maximum curr_reg_size = findTimeOfHalfMaxLeft(*pZXings, gVslide, gBSLMean, beginRisingEdge, endRisingEdge, verbosityLevel ); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); //sometimes pulses where found with a very flat slope, DISCARD them curr_reg_size = removeRegionWithToFlatSlope(*pZXings, risingEdgeMinSlope); pixel[pixelID]->hDiscartedPulses->Fill(step++, init_reg_size - curr_reg_size); if (verbosityLevel > 2) cout << "...done" << endl; // Fill arrival times into histograms vector::iterator it = pZXings->begin(); while( it != pZXings->end() ) { pixel[pixelID]->hMaxPos->Fill( it->maxPos); pixel[pixelID]->hEdgePos->Fill( it->halfRisingEdgePos); pixel[pixelID]->hEdgeSlope->Fill( it->slopeOfRisingEdge); pixel[pixelID]->hIntercept->Fill( it->interceptRisingEdge); pixel[pixelID]->hEdgeLength->Fill( it->lengthOfRisingEdge); pixel[pixelID]->hMaxAmpl->Fill( it->maxVal); ++it; } //----------------------------------------------------------------------------- // 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, select_pulses, 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"; SaveList( //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, bool select_pulses, int verbosityLevel ) { //entry counters int maxOverlayEntries = 0; int edgeOverlayEntries = 0; int maxProfileEntries = 0; int edgeProfileEntries = 0; 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)gVcorr.size() ) Right = gVcorr.size(); if (EdgeRight > (int)gVcorr.size() ) EdgeRight = gVcorr.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, gBSLMean, maxPulseOrder, verbosityLevel ); if ( use_this_peak) { order_of_pulse = order; break; } } if (select_pulses){ order_of_pulse = 0; use_this_peak = true; } //Fill histograms if (use_this_peak) { //get number of entries of used histograms maxOverlayEntries = CurrentPixel->hMaxOverlay[order_of_pulse]->GetEntries(); edgeOverlayEntries = CurrentPixel->hEdgeOverlay[order_of_pulse]->GetEntries(); maxProfileEntries = CurrentPixel->hMaxProfile[order_of_pulse]->GetEntries(); edgeProfileEntries = CurrentPixel->hEdgeProfile[order_of_pulse]->GetEntries(); //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 ) ; } float average_of_max = gVcorr[reg->maxPos]; average_of_max += gVcorr[reg->maxPos + 1]; average_of_max += gVcorr[reg->maxPos - 1]; average_of_max /= 3; float slope_of_follow = 0.0; //Histograms for Maximum Overlay for ( int pos = Left; pos < Right; pos++) { if (pos > reg->maxPos + 1) { if ( pos+5 >= Right ) { break; } if (gVcorr[pos+1] > average_of_max) { break; } slope_of_follow = (gVcorr[pos+5] - gVcorr[pos+1])/5; if (slope_of_follow > 1) { break; } } CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), gVcorr[pos]) ; CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), gVcorr[pos]) ; } //Histograms for Edge Overlay int last_pos = 0; bool ap_cutoff = false; slope_of_follow = 0.0; for ( int pos = EdgeLeft; pos < EdgeRight; pos++) { if (pos > reg->maxPos + 1) { // check if the current slice is close to the end of the overlay window if ( pos+5 >= Right ) { ap_cutoff = false; break; } // check if the next slice is larger than the pulse maximum // i guess it makes more sense to compare the average amplitude of the slices before //to the following pulse if (gVcorr[pos+1] > average_of_max) { ap_cutoff = true; break; } // check if the next 5 slices have a steep slope slope_of_follow = (gVcorr[pos+5] - gVcorr[pos+1])/5; if (slope_of_follow > 1) { ap_cutoff = true; break; } } // if (pos > reg->maxPos // && gVcorr[pos] < average_of_max // && gVcorr[pos+3] > average_of_max // && gVcorr[pos+10] > average_of_max // ) // { // break; // } CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gVcorr[pos]) ; CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gVcorr[pos]) ; } CurrentPixel->hPulseLength->Fill(last_pos - (reg->halfRisingEdgePos)); if (ap_cutoff) { CurrentPixel->hPulseLengthAPcutoff->Fill(last_pos - (reg->halfRisingEdgePos)); } else { CurrentPixel->hPulseLengthTLcutoff->Fill(last_pos - (reg->halfRisingEdgePos)); } //set histogram entries CurrentPixel->hMaxOverlay[order_of_pulse]->SetEntries( maxOverlayEntries +1 ); CurrentPixel->hEdgeOverlay[order_of_pulse]->SetEntries( edgeOverlayEntries +1 ); CurrentPixel->hMaxProfile[order_of_pulse]->SetEntries( maxProfileEntries +1 ); CurrentPixel->hEdgeProfile[order_of_pulse]->SetEntries( edgeProfileEntries +1 ); 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, float bslMean, int maxPulseOrder, int verbosityLevel ) { //determine if pulse is of given order bool use_this_peak = false; if (gVcorr[maxPos] >= ((gainMean*(order+1)) - bslMean - (AmplWindowWidth/2)) && gVcorr[maxPos] < ((gainMean*(order+1)) -bslMean + 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