Ignore:
Timestamp:
05/10/12 16:44:09 (13 years ago)
Author:
Jens Buss
Message:
All funktions implemented to compute template
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C

    r13465 r13641  
    11/********************** FPulseTemplate ***********************
    2 * read FACT raw data
    3 * remove spikes
    4 * calculate baseline
    5 * find peaks (CFD and normal discriminator)
     2* open a root-file with pulse overlay histograms
     3* calculate most propable Pulse-Shape
    64* compute pulse height and pulse integral spektrum of the peaks
    7 + search for Peaks in data
    8 + put peaks into different histos depending und Amplitude
    95+ draw histos for single, double, tripple, ....above 100mV Photonpulses
    106+ draw Parameterdevelopment depending on photon quantity
    … …  
    1915#include <TROOT.h>
    2016#include <TCanvas.h>
    21 #include <TProfile.h>
    2217#include <TTimer.h>
    23 #include <TH1F.h>
    24 #include <TH2F.h>
    2518#include <Getline.h>
    26 #include <TLine.h>
    27 #include <TBox.h>
    2819#include <TMath.h>
    2920#include <TFile.h>
    3021#include <TStyle.h>
    3122#include <TString.h>
    32 #include <TProfile.h>
    3323
    3424#include <stdio.h>
    … …  
    4838// rootmacros
    4939//----------------------------------------------------------------------------
    50 
    51 #include "../fits.h"
    52 
    53 #include "../openFits.h"
    54 #include "../openFits.c"
    55 
    56 #include "../discriminator.h"
    57 #include "../discriminator.C"
    58 #include "../zerosearch.h"
    59 #include "../zerosearch.C"
    60 #include "../factfir.C"
    61 
    62 #include "../DrsCalibration.C"
    63 #include "../DrsCalibration.h"
    64 
    65 #include "../SpikeRemoval.h"
    66 #include "../SpikeRemoval.C"
     40#include "rootfilehandler.h"
     41#include "pixel.h"
     42#include "templateextractors.h"
    6743
    6844//----------------------------------------------------------------------------
    6945// Decleration of global variables
    7046//----------------------------------------------------------------------------
    71 
    7247bool breakout=false;
    7348
    74 int gNEvents;
    75 vector<int16_t> AllPixelDataVector;
    76 vector<int16_t> StartCellVector;
    77 unsigned int CurrentEventID;
    78 size_t PXLxROI;
    79 UInt_t gRegionOfInterest;
    80 UInt_t NumberOfPixels;
    81 TString histotitle;
    82 
    83 size_t TriggerOffsetROI, RC;
    84 size_t drs_n;
    85 vector<float> drs_basemean;
    86 vector<float> drs_gainmean;
    87 vector<float> drs_triggeroffsetmean;
    88 
    89 vector<float> Ameas(FAD_MAX_SAMPLES);   // copy of the data (measured amplitude
    90 vector<float> Vcorr(FAD_MAX_SAMPLES);   // corrected Values
    91 vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
    92 vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
    93 vector<float> Vcfd2(FAD_MAX_SAMPLES);   // CDF result + 2nd sliding average
    94 
    95 float gGainMean = 9;
    96 float gBSLMean = 0;
    97 
    98 typedef struct{
    99   double maxAmpl;
    100   int countOfMax;
    101   } OverlayMaximum;
    102 
    103 //typedef struct{
    104 //  TString name;
    105 //  TString title;
    106 //  TString xTitle;
    107 //  TString yTitle;
    108 //  float yMax;
    109 //  float yMin;
    110 //} histAttrib_t;
    111 
    112 //histAttrib_t* gHistoAttrib[MAX_PULS_ORDER];
    113 //histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER];
    114 //histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER];
    115 //histAttrib_t* gProfileAttrib[MAX_PULS_ORDER];
    116 
    117 
    118 // histograms
    119 const int Number_Of_Debug_Histo_Types = 7;
    120 
    121 const unsigned int
    122         Ameas_  = 0,
    123         N1mean_ = 1,
    124         Vcorr_  = 2,
    125         Vtest_  = 3,
    126         Vslide_ = 4,
    127         Vcfd_   = 5,
    128         Vcfd2_  = 6;
    129 
    130 //const char* gHistoTypes[8] = {
    131 //    "hPixelOverlay",
    132 //    "hPixelMax",
    133 //    "hPixelEdgeOverlay",
    134 //    "hPixelProfile",
    135 //    "hAllPixelOverlay",
    136 //    "hAllPixelMax",
    137 //    "hAllPixelMaxGaus",
    138 //    "hAllPixelProfile"
    139 //};
    140 
    141 //----------------------------------------------------------------------------
    142 // Initialisation of histograms
     49//----------------------------------------------------------------------------
     50// Initialisation of Root Objects
    14351//----------------------------------------------------------------------------
    14452
    14553// Temporary Objects
    146 TH1F*       debugHistos = NULL;
    147 TH2F*       hOverlayTemp = NULL;
    148 TH1F*       hMaximumTemp = NULL;
    149 TH1D*       hProjPeak = NULL;
    150 TH1F*       hTesthisto = NULL;
    151 TH2F*       hTesthisto2 = NULL;
    15254
    15355//Pixelwise Histograms
    154 TH2F*       hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks
    155 TH2F*       hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
    156 TProfile*   hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
    157 TProfile*   hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
    158 TH1F*       hPixelMax[MAX_PULS_ORDER] = {NULL};
    159 TH1F*       hPixelMedian[MAX_PULS_ORDER] = {NULL};
    160 TH1F*       hPixelMean[MAX_PULS_ORDER] = {NULL};
    16156
    16257//All Pixel Histograms
    163 TH2F*       hAllPixelOverlay[MAX_PULS_ORDER] = {NULL};
    164 TProfile*   hAllPixelProfile[MAX_PULS_ORDER] = {NULL};
    165 TProfile*   hAllPixelProfile2[MAX_PULS_ORDER] = {NULL};
    166 TH1F*       hAllPixelMax[MAX_PULS_ORDER]     = {NULL};
    167 TH1F*       hAllPixelMedian[MAX_PULS_ORDER] = {NULL};
    168 TH1F*       hAllPixelMean[MAX_PULS_ORDER] = {NULL};
    169 TH2F*       hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
    170 
    171 //Histogram Parameters
    172 Int_t       gPixelOverlayXaxisLeft = 0;
    173 Int_t       gPixelOverlayXaxisRight = 0;
    17458
    17559//Root-File Objects
    176 TObjArray*  hList = NULL;
    177 TObjArray*  hAllPixelList = NULL;
    178 TObjArray*  hRootList = NULL;
     60TObjArray*  hAllPixelList       = NULL;
     61TObjArray*  hRootList           = NULL;
    17962
    18063//Canvases
    181 TCanvas*    cgpPixelPulses[MAX_PULS_ORDER] = {NULL};
    182 TCanvas*    cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL};
    183 TCanvas*    cgpTestHistos = NULL;
    184 //TCanvas*    gpcDevelopment = NULL;
     64TCanvas**   cgpPixelPulses      = NULL;
     65TCanvas**   cgpAllPixelPulses   = NULL;
     66TCanvas**   cgpDistributions    = NULL;
     67//TCanvas*    cgpTestHistos      = NULL;
    18568
    18669//----------------------------------------------------------------------------
    18770// Functions
    18871//----------------------------------------------------------------------------
    189 
    190 void BookHistos(int );
    191 void BookPixelHistos(int );
    192 void BookTestHistos( int );
    193 
    194 void DeletePixelHistos(int );
    195 void SaveHistograms( const char*, const char*, TObjArray*, int );
    196 
    197 void FillHistograms(vector<Region>*, int, int, int);
    198 void DrawPulseHistograms(int, int);
    199 void DrawMaximumHistograms( int, int );
    200 void DrawTestHistograms( int);
    201 
    202 void PlotPulseShapeOfMaxPropability(unsigned int, int, int);
    203 void FitMaxPropabilityPuls( TH1F* , int );
    204 
    205 Double_t MedianOfH1( TH1 *h1 );
    206 void PlotMedianEachSliceOfPulse(unsigned int, int, int);
    207 
    208 void UpdateCanvases( int, int);
    209 void DeletePixelCanvases( int );
    210 
    211 void CreateRootFile(const char*, int );
    212 TFile *OpenRootFile(const char*, int );
    213 void CloseRootFile(TFile *tf);
    214 
    215 TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file
    216 TString CreateSubDirName(const char* );
    217 
    218 void WritePixelTemplateToCsv( TString, const char*, const char*, int, int);
    219 void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int);
    220 //void SetHistogrammSettings( const char*, int, int , int);
    221 //void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int);
    222 //void DeletePulseCanvas( TCanvas* , unsigned int);
    223 
    224 
     72void DeletePixelCanvases( int, int );
     73void
     74UpdateCanvases(
     75        int verbosityLevel,
     76        int max_pulse_order,
     77        bool testmode
     78        );
    22579//----------------------------------------------------------------------------
    22680//----------------------------------------------------------------------------
    … …  
    22983//----------------------------------------------------------------------------
    23084int FCalcPulseTemplate(
    231   char*         datafilename        = "../data/2011/11/09/20111109_006.fits.gz",
    232   const char*   drsfilename         = "../data/2011/11/09/20111109_003.drs.fits.gz",
    233   const char*   OutRootFileName     = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root",
    234   const char*   OutPutPath          = "../analysis/FPulseTemplate/20111109_006/Templates/",
    235   bool          ProduceGraphic      = false,
    236   int           refresh_rate        = 500,      //refresh rate for canvases
    237   bool          spikeDebug          = false,
    238   bool          debugPixel          = false,
    239   bool          fitdata             = false,
    240   int           verbosityLevel      = 0,        // different verbosity levels can be implemented here
    241   int           firstevent          = 0,
    242   int           nevents             = -1,
    243   int           firstpixel          = 0,
    244   int           npixel              = -1,
    245   int           sampleSetSize       = 36,
    246   int           AmplWindowWidth     = 14,       //Width of Window for selection of pulses to histograms
    247   float         GainMean           = 8,
    248   float         BSLMean            = -1,        //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
    249   int           avg1                = 8,
    250   int           avg2                = 8,
    251   int           OverlayWindowLeft   = 70,
    252   int           OverlayWindowRight  = 230
     85    TString     InRootFileName        = "20111109_006.pulses.root",
     86    TString     InputPath           = "../analysis/FPulseTemplate/20111109_006/",
     87    TString     OutputRootFileName     = "20111109_006.pulses.root",
     88    TString     OutPutPath          = "../analysis/FPulseTemplate/20111109_006/Templates/",
     89    bool        ProduceGraphic      = false,
     90    bool        fitdata             = false,
     91    bool        stats               = false,
     92    bool        debugPixel          = false,
     93    int         pixelSetSize        = 40,
     94    int         maxPulseOrder       = 3,
     95    int         refresh_rate        = 500,      //refresh rate for canvases
     96    int         verbosityLevel      = 0,        // different verbosity levels can be implemented here
     97    int         firstpixel          = 0,
     98    int         npixel              = -1
    25399        )
    254100{
    255     gGainMean = GainMean;
    256     gBSLMean = BSLMean;
    257 //----------------------------------------------------------------------------
    258 //      Save-Root-File Settings
    259 //----------------------------------------------------------------------------
    260     CreateRootFile( OutRootFileName, verbosityLevel );
    261 
     101
     102//----------------------------------------------------------------------------
     103//      Open-Root-File Settings
     104//----------------------------------------------------------------------------
     105    TFile * inputRootFile = OpenRootFile( InputPath, InRootFileName, verbosityLevel );
     106
     107//----------------------------------------------------------------------------
     108//      global variable Settings
     109//----------------------------------------------------------------------------
     110    if ( npixel == -1 )
     111    {
     112        npixel = NPIX - firstpixel;
     113    }
     114
     115
     116    if ( pixelSetSize == -1 )
     117    {
     118        pixelSetSize = firstpixel +npixel;
     119    }
     120//    float GainMean  = GainMean;  // this has to be extracted from root files
     121//    float BSLMean   = BSLMean;    // this has to be extracted from root files
    262122
    263123//----------------------------------------------------------------------------
    264124//      root global Settings
    265125//----------------------------------------------------------------------------
    266 
    267     gPixelOverlayXaxisLeft = OverlayWindowLeft;
    268     gPixelOverlayXaxisRight = OverlayWindowRight;
    269 
    270126    gStyle->SetPalette(1,0);
    271127    gROOT->SetStyle("Plain");
    272 //    gPad->SetGrid();
    273 
     128    gPad->SetGrid();
     129
     130
     131
     132//----------------------------------------------------------------------------
     133//      root Canvas Settings
     134//----------------------------------------------------------------------------
     135
     136    //Canvas Pad numbering
     137    int pulsesCanvasFrameNrs[4] = {
     138        1,  // Top left
     139        2,  // Top right
     140        3,  // bottom left
     141        4   // bootom right
     142    };
     143
     144    //Canvas Pad numbering
     145    int distributionCanvasFrameNrs[4] = {
     146        1,  // Top left
     147        2,  // Top right
     148        3,  // bottom left
     149        4   // bootom right
     150    };
     151
     152    if (ProduceGraphic)
     153    {
     154
     155        //Canvases
     156        cgpPixelPulses      = new TCanvas*[maxPulseOrder];
     157        cgpDistributions    = new TCanvas*[maxPulseOrder];
     158
     159        //TCanvas*    gpcDevelopment = NULL;
     160        TString cName   = "";
     161        TString cTitle  = "";
     162
     163        //naming of pulse canvases
    274164        for (
    275              int pulse_order = MAX_PULS_ORDER - 1;
     165             int pulse_order = maxPulseOrder - 1;
    276166             pulse_order >= 0 ;
    277167             pulse_order--
    278168              )
    279169        {
    280             TString cName   ="cgpPixelPulses";
    281                     cName   += pulse_order;
    282             TString cTitle   ="Pulses of Order: ";
    283                     cTitle   += pulse_order;
    284             cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800);
     170            cName   ="cgpDistributions";
     171            cName   += pulse_order;
     172            cTitle   ="Distributions of Pulses with Order of: ";
     173            cTitle   += pulse_order;
     174            cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720);
     175            cgpDistributions[pulse_order]->Divide(2, 2);
     176            cName   ="cgpPixelPulses";
     177            cName   += pulse_order;
     178            cTitle   ="Overlays of Pulses with Order of: ";
     179            cTitle   += pulse_order;
     180            cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720);
    285181            cgpPixelPulses[pulse_order]->Divide(2, 2);
    286                     cName   = "cgpAllPixelPulses";
    287                     cName   += pulse_order;
    288                     cTitle   ="AllPixel average of puls shapes of Order: ";
    289                     cTitle   += pulse_order;
    290             cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 );
    291             cgpAllPixelPulses[pulse_order]->Divide(2, 2);
    292182        }
     183
    293184
    294185        // Create (pointer to) Canvases, which are used in every run,
    … …  
    297188                // the pointers anyway ...
    298189
    299         TCanvas *cFiltered = NULL;
    300         if (spikeDebug)
    301         {
    302           cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300);
    303           cFiltered->Divide(1, 3);
    304         }
    305 
    306         cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 );
    307         cgpTestHistos->Divide(2,0);
     190//        if (testmode)
     191//        {
     192//        //additional Test histograms
     193//        cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 );
     194//        cgpTestHistos->Divide(2,0);
     195//        }
     196    }
     197
    308198//-----------------------------------------------------------------------------
    309199// Filter-Settings
    310200//-----------------------------------------------------------------------------
    311 // CFD filter settings
    312     int k_cfd = 10;
    313     vector<double> a_cfd(k_cfd, 0);
    314     double b_cfd = 1.;
    315     a_cfd[0]=-0.75;
    316     a_cfd[k_cfd-1]=1.;
     201
    317202
    318203//-----------------------------------------------------------------------------
    … …  
    320205//-----------------------------------------------------------------------------
    321206
    322 // Open the data file
    323 
    324     fits * datafile;
    325     // Opens the raw data file and 'binds' the variables given as
    326     // Parameters to the data file. So they are filled with
    327     // raw data as soon as datafile->GetRow(int) is called.
    328     gNEvents = openDataFits(
    329                     datafilename,
    330                     &datafile,
    331                     AllPixelDataVector,
    332                     StartCellVector,
    333                     CurrentEventID,
    334                     gRegionOfInterest,
    335                     NumberOfPixels,
    336                     PXLxROI,
    337                     verbosityLevel
    338                 );
    339 
    340     if (gNEvents == 0)
    341     {
    342         cout << "return code of OpenDataFile:" << datafilename<< endl;
    343         cout << "is zero -> aborting." << endl;
    344         return 1;
    345     }
    346 
    347     if (verbosityLevel > 0)
    348     {
    349         cout << endl <<"number of events in file: "<< gNEvents << "\t";
    350     }
    351 
    352     if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!
    353 
    354     if (verbosityLevel > 0)
    355     {
    356         cout <<"of, which "<< nevents << " will be processed"<< endl;
    357         cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
    358     }
    359 
    360     if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
    361 
    362     if (verbosityLevel > 0)
    363     {
    364         cout <<"of, which "<< npixel << " will be processed"<< endl;
    365     }
    366 
    367     //Get the DRS calibration
    368     RC = openCalibFits(
    369                 drsfilename,
    370                 drs_basemean,
    371                 drs_gainmean,
    372                 drs_triggeroffsetmean,
    373                 TriggerOffsetROI
    374                 );
    375 
    376     if (RC == 0)
    377     {
    378      cout << "return code of openCalibFits:" << drsfilename << endl;
    379      cout << "is zero -> aborting." << endl;
    380      return 1;
    381     }
    382 
    383     // Book the histograms
    384     BookHistos(verbosityLevel );
    385     BookTestHistos( verbosityLevel );
     207//----------------------------------------------------------------------------
     208// Initialize Pixel
     209//----------------------------------------------------------------------------
     210    Pixel** pixel = new Pixel*[NPIX];
     211
     212    for (int i = 0 ; i < NPIX; i++)
     213    {
     214        pixel[i] = NULL;
     215    }
     216
    386217//-------------------------------------
    387218// Loop over Pixel Sets
    388219//-------------------------------------
    389     for ( int firstPixelOfSample = firstpixel;
    390           firstPixelOfSample < firstpixel + npixel;
    391           firstPixelOfSample += sampleSetSize )
    392     {
    393 
     220    for ( int firstPixelOfSet = firstpixel;
     221          firstPixelOfSet < firstpixel + npixel;
     222          firstPixelOfSet += pixelSetSize )
     223    {
    394224        if (verbosityLevel >= 0)
    395225        {
    396226            cout << "------------------------------------------------" << endl
    397                  << "...processing Sample from Pixel: "
    398                  << firstPixelOfSample;
     227                 << "...processing Pixel: "
     228                 << firstPixelOfSet
    399229                 << "to Pixel: "
    400                  << firstPixelOfSample+sampleSetSize-1 << endl;
     230                 << firstPixelOfSet+pixelSetSize-1 << endl;
    401231        }
    402         BookPixelHistos(verbosityLevel );
    403 
    404 //--------------------------------------------------------------------
    405 // Loops over Every Event of Pixel
    406 //--------------------------------------------------------------------
    407         for ( int ev = firstevent; ev < firstevent + nevents; ev++)
     232
     233        //--------------------------------------------------------------------
     234        // Loops over every Pixel of a Set of Pixels
     235        //--------------------------------------------------------------------
     236        for ( int pixelID = firstPixelOfSet;
     237                pixelID < firstPixelOfSet + pixelSetSize
     238                && pixelID < firstpixel + npixel;
     239                pixelID++ )
    408240        {
    409             // Get an Event --> consists of 1440 Pixel ...erm....data
    410             datafile->GetRow( ev );
    411 
    412             if (verbosityLevel > 0)
     241
     242            if (verbosityLevel > 1)
    413243            {
    414244             cout << "-------------------------------------" << endl
    415                   << "...processing Event: " << CurrentEventID
    416                   << "/" << nevents << " of Pixel: " << pixel << endl;
     245                  << "...processing Set from Pixel "
     246                  << firstPixelOfSet
     247                  << " to Pixel "
     248                  << firstPixelOfSet+pixelSetSize-1
     249                  << " Pixel: " << pixelID
     250                  << "/" << firstpixel + npixel -1 << endl;
    417251            }
    418252
    419 //--------------------------------------------------------------------
    420 // Loops over every Pixel of a Set of Pixels
    421 //--------------------------------------------------------------------
    422         for ( int pixel = firstPixelOfSample;
    423             pixel < sampleSetSize || pixel < firstpixel + npixel;
    424             pixel++ )
    425         {
    426             if (verbosityLevel > 0)
    427            {
    428              cout << "-------------------------------------" << endl
    429                   << "...processing Pixel: " << pixel
    430                   << "/" << nPixelSample << endl;
    431            }
    432 
    433 //-------------------------------------
    434 // Apply Calibration
    435 //-------------------------------------
    436             if (verbosityLevel > 2) cout << "...applying DrsCalibration";
    437             applyDrsCalibration(
    438                         Ameas,
    439                         pixel,
    440                         12,
    441                         12,
    442                         drs_basemean,
    443                         drs_gainmean,
    444                         drs_triggeroffsetmean,
    445                         gRegionOfInterest,
    446                         AllPixelDataVector,
    447                         StartCellVector
     253            //-------------------------------------
     254            // Create individual Pixel
     255            //-------------------------------------
     256           pixel[pixelID] = new Pixel(
     257                       pixelID,
     258                       maxPulseOrder,  ///One should get this from the file
     259                       verbosityLevel,
     260                       stats,
     261                       inputRootFile
     262                    );
     263
     264            if (breakout)   break;
     265
     266            //-------------------------------------
     267            // Histogramms of Maximas in Overlay Spectra
     268            //-------------------------------------
     269
     270            for ( int pulse_order = 0 ;
     271                  pulse_order < pixel[pixelID]->mMaxPulseOrder ;
     272                  pulse_order ++)
     273            {
     274                if (verbosityLevel > 2)
     275                {
     276                    cout << "-------------------------------------" << endl
     277                         << "...processing Set from Pixel "
     278                         << firstPixelOfSet
     279                         << " to Pixel "
     280                         << firstPixelOfSet+pixelSetSize-1
     281                         << " Pixel: " << pixelID
     282                         << "/" << firstpixel + npixel -1 << endl
     283                         << " Pulse-Order: " << pulse_order;
     284                }
     285
     286                // Calculate Max Prop. Value of each slice
     287                //-------------------------------------
     288
     289                //from Maximum Overlay
     290                ExtractPulseTemplate(
     291                            pixel[pixelID],
     292                            "Maximum",
     293                            pulse_order,
     294                            verbosityLevel
     295                            );
     296
     297                //from Edge Overlay
     298                ExtractPulseTemplate(
     299                            pixel[pixelID],
     300                            "Edge",
     301                            pulse_order,
     302                            verbosityLevel
     303                            );
     304
     305                WritePixelTemplateToCsv(
     306                            pixel[pixelID],
     307                            OutPutPath,
     308                            "Maximum",
     309                            pulse_order,
     310                            verbosityLevel
     311                            );
     312
     313                WritePixelTemplateToCsv(
     314                            pixel[pixelID],
     315                            OutPutPath,
     316                            "Edge",
     317                            pulse_order,
     318                            verbosityLevel
     319                            );
     320
     321                pixel[pixelID]->DrawTemplateHistograms(
     322                            cgpPixelPulses,
     323                            pulsesCanvasFrameNrs
     324                            );
     325
     326                pixel[pixelID]->DrawEdgeTemplateHistograms(
     327                            cgpPixelPulses,
     328                            pulsesCanvasFrameNrs
     329                            );
     330
     331
     332
     333            }
     334            // End of Loop over pulsorder of current Pixel
     335
     336            if (ProduceGraphic)
     337            {
     338            UpdateCanvases(
     339                        verbosityLevel,
     340                        MAX_PULS_ORDER,
     341                        false
    448342                        );
    449             if (verbosityLevel > 2) cout << "...done " << endl;
    450 
    451 //-------------------------------------
    452 // Apply Filters
    453 //-------------------------------------
    454             // finds spikes in the raw data, and interpolates the value
    455             // spikes are: 1 or 2 slice wide, positive non physical artifacts
    456             if (verbosityLevel > 2) cout << "...removeing Spikes";
    457             removeSpikes (Ameas, Vcorr);
    458             if (verbosityLevel > 2) cout << "...done " << endl;
    459 
    460             // filter Vcorr with sliding average using FIR filter function
    461             if (verbosityLevel > 2) cout << "...applying sliding average filter";
    462             sliding_avg(Vcorr, Vslide, avg1);
    463             if (verbosityLevel > 2) cout << "...done " << endl;
    464 
    465             // filter Vslide with CFD using FIR filter function
    466             if (verbosityLevel > 2) cout << "...apllying factfir filter";
    467             factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
    468             if (verbosityLevel > 2) cout << "...done " << endl;
    469 
    470             // filter Vcfd with sliding average using FIR filter function
    471             if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
    472             sliding_avg(Vcfd, Vcfd2, avg2);
    473             if (verbosityLevel > 2) cout << "...done " << endl;
    474 
    475 //-------------------------------------
    476 // Search vor Zero crossings
    477 //-------------------------------------
    478             if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ;
    479             // peaks in Ameas[] are found by searching for zero crossings
    480             // in Vcfd2
    481             // first Argument 1 means ... *rising* edge
    482             // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
    483             vector<Region>* pZXings = zerosearch( Vcfd2 , 1 , 8);
    484             // pZXings means "zero cross ings"
    485             EnlargeRegion(*pZXings, 10, 10);
    486             findAbsMaxInRegions(*pZXings, Vslide);
    487             removeMaximaBelow( *pZXings, 3.0);
    488             removeRegionWithMaxOnEdge( *pZXings, 2);
    489             removeRegionOnFallingEdge( *pZXings, 100);
    490             findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel );
    491             if (verbosityLevel > 2) cout << "...done" << endl;
    492 
    493 //-----------------------------------------------------------------------------
    494 // Fill Overlay Histos
    495 //-----------------------------------------------------------------------------
    496         FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel );
    497 
    498 //-----------------------------------------------------------------------------
    499 // Spike Debug
    500 //-----------------------------------------------------------------------------
    501         if ( spikeDebug )
    502         {
    503             // TODO do this correct. The vectors should be the rigt ones... this is just luck
    504             debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
    505             debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
    506             debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
    507             debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
    508             debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
    509 
    510             for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++)
     343            }
     344
     345            if (debugPixel)
    511346            {
    512                debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
    513                debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
    514                debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
    515                debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
    516                debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
     347                //Process gui events asynchronously during input
     348                cout << endl;
     349                TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
     350                timer.TurnOn();
     351                TString input = Getline("Type 'q' to exit, <return> to go on: ");
     352                timer.TurnOff();
     353                if (input=="q\n")
     354                {
     355                       break;
     356                }
    517357            }
    518358
    519 
    520             cFiltered->cd(1);
    521             gPad->SetGrid();
    522             debugHistos[Ameas_].Draw();
    523 
    524             cFiltered->cd(2);
    525             gPad->SetGrid();
    526             debugHistos[Vcorr_].Draw();
    527 
    528             cFiltered->cd(3);
    529             gPad->SetGrid();
    530             debugHistos[Vslide_].Draw();
    531 
    532             TBox *OneBox;
    533             vector<TBox*> MyBoxes;
    534             for (unsigned int i=0; i<pZXings->size(); i++){
    535                     OneBox = new TBox(
    536                             pZXings->at(i).maxPos -10 ,
    537                             pZXings->at(i).maxVal -0.5,
    538                             pZXings->at(i).maxPos +10 ,
    539                             pZXings->at(i).maxVal +0.5);
    540                     OneBox->SetLineColor(kBlue);
    541                     OneBox->SetLineWidth(1);
    542                     OneBox->SetFillStyle(0);
    543                     OneBox->SetFillColor(kRed);
    544                     MyBoxes.push_back(OneBox);
    545                     OneBox->Draw();
    546             }
    547 
    548 //            cFiltered->cd(3);
    549 //            gPad->SetGrid();
    550 //            debugHistos[Vcfd2_].Draw();
    551 //            TLine *zeroline = new TLine(0, 0, 1024, 0);
    552 //            zeroline->SetLineColor(kBlue);
    553 //            zeroline->Draw();
    554 
    555             cFiltered->Update();
    556 
    557 
    558             //Process gui events asynchronously during input
    559             TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
    560             timer.TurnOn();
    561             TString input = Getline("Type 'q' to exit, <return> to go on: ");
    562             timer.TurnOff();
    563             if (input=="q\n") {
    564                     breakout=true;
    565             }
    566 
    567             //TODO!!!!!!!!!
    568             // do some Garbage collection ...
    569             // all the Objects on the heap should be deleted here.
    570 
    571         }// end of if(spikeDebug)
    572 
    573         delete pZXings;
    574         if (breakout)   break;
    575 //-------------------------------------
    576 // Draw 1. Set of Pixel Histograms
    577 //-------------------------------------
    578             if ((ev % refresh_rate) == 0)
     359            if (verbosityLevel > 2)
    579360            {
    580                 DrawPulseHistograms(
    581                         verbosityLevel,
    582                         MAX_PULS_ORDER
    583                         );
    584                 DrawTestHistograms(
    585                         verbosityLevel
    586                         );
    587               if (ProduceGraphic)
    588               {
    589                 UpdateCanvases(
    590                                 verbosityLevel,
    591                                 MAX_PULS_ORDER
    592                                 );
    593               }
    594             }
    595 
    596         if (breakout) break;
    597 
    598         if (verbosityLevel > 2)
    599         {
    600             cout << endl << "-------------------------------------"<< endl;
    601         }
    602 
    603 
    604 
    605 
    606         }//End of Loop over Sample
    607 //-------------------------------------
    608 // end of Loops over Sample
    609 //-------------------------------------
    610 
    611         }//End of Loop over Events
    612 //-------------------------------------
    613 // end of Loops over Events
    614 //-------------------------------------
    615 
    616 
    617 
    618 
    619 //-------------------------------------
    620 // Histogramms of Maximas in Overlay Spectra
    621 //-------------------------------------
    622 
    623     PlotPulseShapeOfMaxPropability(
    624             MAX_PULS_ORDER,
    625             fitdata,
    626             verbosityLevel
    627          );
    628 
    629     for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    630     {
    631     if (verbosityLevel > 2) cout << "\t...plotting most probable pulse of"
    632                                  << "pulse order " << pulse_order;
    633     PlotMedianEachSliceOfPulse(
    634          hPixelOverlay[pulse_order],
    635          hPixelMedian[pulse_order],
    636          fitdata,
    637          verbosityLevel
    638          );
    639 
    640     PlotMeanEachSliceOfPulse(
    641          hPixelOverlay[pulse_order],
    642          hPixelMean[pulse_order],
    643          fitdata,
    644          verbosityLevel
    645          );
    646     }
    647          WritePixelTemplateToCsv(
    648                      OutPutPath,
    649                      "PulseTemplate_PointSet",
    650                      "Maximum",
    651                      pixel,
    652                      verbosityLevel
    653                      );
    654 
    655        if (verbosityLevel > 2) cout << "...done" << endl;
    656        //here is what happends at the end of each loop over all pixels
    657        //Save Histograms of Pixel into Output rootfile
    658        DrawPulseHistograms(
    659                verbosityLevel,
    660                MAX_PULS_ORDER
    661                );
    662 
    663        DrawMaximumHistograms(
    664                verbosityLevel,
    665                MAX_PULS_ORDER
    666                );
    667 
    668         if (ProduceGraphic)
    669         {
    670         UpdateCanvases(
    671                        verbosityLevel,
    672                        MAX_PULS_ORDER
    673                        );
    674         }
    675 
    676 //       SaveHistograms(
    677 //                   OutRootFileName,
    678 //                   CreateSubDirName(pixel),
    679 //                   hList,
    680 //                   verbosityLevel
    681 //                   );
    682 
    683        if (debugPixel)
    684        {
    685            //Process gui events asynchronously during input
    686            cout << endl;
    687            TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
    688            timer.TurnOn();
    689            TString input = Getline("Type 'q' to exit, <return> to go on: ");
    690            timer.TurnOff();
    691            if (input=="q\n") {
    692                    break;
    693            }
    694        }
    695 
    696        DeletePixelHistos(verbosityLevel);
    697 
    698        if (verbosityLevel > 2)
    699        {
    700            cout << endl << "...End of Pixel"
     361            cout << endl << "...End of Pixel"
    701362                << endl << "------------------------------------------------"
    702363                << endl;
    703        }
    704     }
    705     // End of Loop over all Pixels
     364            }
     365        }
     366        // End of Loop over all Pixels of set
     367
     368        if (verbosityLevel > 1)
     369        {
     370        cout << endl << "...End of Loop over all Pixels of set"
     371            << endl << "------------------------------------------------"
     372            << endl;
     373        }
     374    }
     375    // End of Loop over all Pixelsets
     376
     377    if (verbosityLevel > 0)
     378    {
     379    cout << endl << "...End of Loop over all Pixelsets"
     380        << endl << "------------------------------------------------"
     381        << endl;
     382    }
    706383
    707384//-------------------------------------
    708 // Draw Histograms
     385// Draw All Pixel Histograms
    709386//-------------------------------------
    710387
    711     SaveHistograms(     //save histograms of all pixel into output root file
    712                 OutRootFileName,
    713                 CreateSubDirName("All"),
    714                 hAllPixelList,
    715                 verbosityLevel
    716                 );
     388//    SaveHistograms(     //save histograms of all pixel into output root file
     389//                OutInRootFileName,
     390//                CreateSubDirName("All"),
     391//                hAllPixelList,
     392//                verbosityLevel
     393//                );
    717394
    718395//    SaveHistograms(     //save histograms of generell results into output root file
    719 //                OutRootFileName,
     396//                OutInRootFileName,
    720397//                "root",
    721398//                hRootList,
    … …  
    725402
    726403
    727     if (ProduceGraphic)
    728     {
    729         UpdateCanvases(
    730                     verbosityLevel,
    731                     MAX_PULS_ORDER
    732                     );
    733     }
    734 
    735     WriteAllPixelTemplateToCsv(
    736                 OutPutPath,
    737                 "PulseTemplate_PointSet",
    738                 "Maximum",
    739                 verbosityLevel
    740                 );
    741 
    742     DeletePixelCanvases( verbosityLevel );
     404//    if (ProduceGraphic)
     405//    {
     406//        UpdateCanvases(
     407//                    verbosityLevel,
     408//                    MAX_PULS_ORDER
     409//                    );
     410//    }
     411
     412//    WriteAllPixelTemplateToCsv(
     413//                OutPutPath,
     414//                "PulseTemplate_PointSet",
     415//                "Maximum",
     416//                verbosityLevel
     417//                );
     418
     419    DeletePixelCanvases( maxPulseOrder ,verbosityLevel );
    743420    return( 0 );
    744421}
    … …  
    753430// Funktions
    754431//-----------------------------------------------------------------------------
    755 
    756 /*
    757 ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel
    758 wie viel speicher braucht das?
    759 Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren
    760 am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute.
    761 ne klasse wäre super
    762 
    763 */
    764432void
    765 BookPixelHistos(
    766         int     sampleSize,
    767 
    768         int     verbosityLevel
    769         )
    770 {
    771     if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
    772     hList = new TObjArray;
    773 
    774     TString histoname;
    775     for (int order = 0; order < MAX_PULS_ORDER; order ++)
    776     {
    777 //        SetHistogramAttributes(
    778 //                    "PeakOverlay",
    779 //                    order,
    780 //                    gOverlayAttrib,
    781 //                    MaxAmplOfFirstPulse
    782 //                    );
    783         histoname = "hPixelOverlay";
    784         histoname += order;
    785         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    786         hPixelOverlay[order] = new TH2F(
    787                     histoname,
    788                     "Overlay of detected pulses of one pulse order for one Pixel",
    789                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    790                     (-1*gPixelOverlayXaxisLeft)-0.5,
    791                     gPixelOverlayXaxisRight-0.5 ,
    792                     512,
    793                     -55.5,
    794                     200.5
    795                     );
    796         hPixelOverlay[order]->SetAxisRange(
    797                     gBSLMean - 5,
    798                     (gGainMean*(order+1)) + 10,
    799                     "Y");
    800         hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    801         hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    802         //hPixelProfile->SetBit(TH2F::kCanRebin);
    803         hList->Add( hPixelOverlay[order] );
    804 
    805     //------------------------------------------------------------------------
    806 //        SetHistogramAttributes(
    807 //                    "PeakMaximum",
    808 //                    order,
    809 //                    gMaximumAttrib,
    810 //                    MaxAmplOfFirstPulse
    811 //                    );
    812         histoname = "hPixelMax";
    813         histoname += order;
    814         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    815         hPixelMax[order] = new TH1F(
    816                     histoname,
    817                     "Maximum value of each slice in overlay plot",
    818                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    819                     (-1*gPixelOverlayXaxisLeft)-0.5,
    820                     gPixelOverlayXaxisRight-0.5
    821                     );
    822         hPixelMax[order]->SetAxisRange(
    823                     gBSLMean - 5,
    824                     (gGainMean*(order+1)) + 10,
    825                     "Y");
    826         hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    827         hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    828         hList->Add( hPixelMax[order] );
    829 
    830 //------------------------------------------------------------------------
    831     //        SetHistogramAttributes(
    832     //                    "PeakMaximum",
    833     //                    order,
    834     //                    gMaximumAttrib,
    835     //                    MaxAmplOfFirstPulse
    836     //                    );
    837             histoname = "hPixelMedian";
    838             histoname += order;
    839             if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    840             hPixelMedian[order] = new TH1F(
    841                         histoname,
    842                         "Median of each slice in overlay plot",
    843                         gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    844                         (-1*gPixelOverlayXaxisLeft)-0.5,
    845                         gPixelOverlayXaxisRight-0.5
    846                         );
    847             hPixelMedian[order]->SetAxisRange(
    848                         gBSLMean - 5,
    849                         (gGainMean*(order+1)) + 10,
    850                         "Y");
    851             hPixelMedian[order]->SetLineColor(kRed);
    852             hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    853             hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    854             hList->Add( hPixelMedian[order] );
    855 
    856 //------------------------------------------------------------------------
    857     //        SetHistogramAttributes(
    858     //                    "PeakMaximum",
    859     //                    order,
    860     //                    gMaximumAttrib,
    861     //                    MaxAmplOfFirstPulse
    862     //                    );
    863             histoname = "hPixelMean";
    864             histoname += order;
    865             if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    866             hPixelMean[order] = new TH1F(
    867                         histoname,
    868                         "Mean of each slice in overlay plot",
    869                         gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    870                         (-1*gPixelOverlayXaxisLeft)-0.5,
    871                         gPixelOverlayXaxisRight-0.5
    872                         );
    873             hPixelMean[order]->SetAxisRange(
    874                         gBSLMean - 5,
    875                         (gGainMean*(order+1)) + 10,
    876                         "Y");
    877             hPixelMean[order]->SetLineColor(kBlue);
    878             hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    879             hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    880             hList->Add( hPixelMean[order] );
    881 
    882 //------------------------------------------------------------------------
    883 //        SetHistogramAttributes(
    884 //                    "PeakMaxGaus",
    885 //                    order,
    886 //                    gMaxGausAttrib,
    887 //                    MaxAmplOfFirstPulse
    888 //                    );
    889         histoname = "hPixelEdgeOverlay";
    890         histoname += order;
    891         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    892         hPixelEdgeOverlay[order] = new TH2F(
    893                     histoname,
    894                     "Overlay at rising edge of detected pulses of one pulse order for one Pixel ",
    895                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    896                     (-1*gPixelOverlayXaxisLeft)-0.5,
    897                     gPixelOverlayXaxisRight-0.5 ,
    898                     512,
    899                     -55.5,
    900                     200.5
    901                     );
    902 
    903         hPixelEdgeOverlay[order]->SetAxisRange(
    904                     gBSLMean - 5,
    905                     (gGainMean*(order+1)) + 10,
    906                     "Y");
    907         hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    908         hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    909         hList->Add( hPixelEdgeOverlay[order] );
    910 
    911     //------------------------------------------------------------------------
    912 //        SetHistogramAttributes(
    913 //                    "PeakProfile",
    914 //                    order,
    915 //                    gProfileAttrib,
    916 //                    MaxAmplOfFirstPulse
    917 //                    );
    918         histoname = "hPixelProfile";
    919         histoname += order;
    920         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    921         hPixelProfile[order] = new TProfile(
    922                     histoname,
    923                     "Mean value of each slice in overlay plot (Tprofile)",
    924                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
    925                     (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
    926                     gPixelOverlayXaxisRight-0.5 ,                    //xup
    927     //                512,                                            //nbinsy
    928                     -55.5,                                          //ylow
    929                     300.5,                                          //yup
    930                     "s");                                           //option
    931         hPixelProfile[order]->SetAxisRange(
    932                     gBSLMean - 5,
    933                     (gGainMean*(order+1)) + 10,
    934                     "Y");
    935         hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    936         hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    937         //hPixelProfile->SetBit(TH2F::kCanRebin);
    938         hList->Add( hPixelProfile[order] );
    939 
    940 
    941     //------------------------------------------------------------------------
    942     //        SetHistogramAttributes(
    943     //                    "PeakProfile",
    944     //                    order,
    945     //                    gProfileAttrib,
    946     //                    MaxAmplOfFirstPulse
    947     //                    );
    948             histoname = "hPixelProfile2";
    949             histoname += order;
    950             if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    951             hPixelProfile2[order] = new TProfile(
    952                         histoname,
    953                         "Mean value of each slice in overlay plot (Tprofile)",
    954                         gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
    955                         (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
    956                         gPixelOverlayXaxisRight-0.5 ,                    //xup
    957         //                512,                                            //nbinsy
    958                         -55.5,                                          //ylow
    959                         300.5,                                          //yup
    960                         "s");                                           //option
    961             hPixelProfile2[order]->SetLineColor(kRed);
    962             hPixelProfile2[order]->SetAxisRange(
    963                         gBSLMean - 5,
    964                         (gGainMean*(order+1)) + 10,
    965                         "Y");
    966             hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    967             hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    968             //hPixelProfile->SetBit(TH2F::kCanRebin);
    969 
    970 
    971 
    972 
    973     }
    974     if (verbosityLevel > 2) cout << "...done" << endl;
    975 }
    976 //end of BookPixelHistos
    977 //----------------------------------------------------------------------------
    978 
    979 
    980 void DeletePixelHistos( int verbosityLevel )
    981 {
    982     if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ;
    983     for (int order = 0; order < MAX_PULS_ORDER; order ++)
    984     {
    985         if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order;
    986         delete hPixelOverlay[order];
    987         if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order;
    988         delete hPixelMax[order];
    989         if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order;
    990         delete hPixelMedian[order];
    991         if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order;
    992         delete hPixelMean[order];
    993         if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order;
    994         delete hPixelEdgeOverlay[order];
    995         if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order;
    996         delete hPixelProfile[order];
    997         if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order;
    998         delete hPixelProfile2[order];
    999     }
    1000     if (verbosityLevel > 3) cout << endl << "...deleting hList";
    1001     delete hList;
    1002     if (verbosityLevel > 2) cout << endl << "...done" << endl;
    1003 }
    1004 //end of DeletePixelHistos
    1005 //----------------------------------------------------------------------------
    1006 
    1007 void BookTestHistos( int verbosityLevel )
    1008 {
    1009     if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
    1010 
    1011     hTesthisto = new TH1F (
    1012                 "hTesthisto",
    1013                 "Deviation of rising edge and maximum",
    1014                 600,
    1015                 -10.1,
    1016                 10.1
    1017                 );
    1018     hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1019     hTesthisto->GetYaxis()->SetTitle( "counts" );
    1020 
    1021     hTesthisto2 = new TH2F (
    1022                 "hTesthisto2",
    1023                 "Deviation of rising edge and maximum by event #",
    1024                 gNEvents,
    1025                 250,
    1026                 800,
    1027                 600,
    1028                 -10.1,
    1029                 10.1
    1030                 );
    1031 //    hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1032 //    hTesthisto2->GetYaxis()->SetTitle( "counts" );
    1033 //    hTesthisto2->SetMarkerStyle( 2 );
    1034 
    1035     if (verbosityLevel > 2) cout << "...done" << endl;
    1036 }
    1037 //end of BookTestHistos
    1038 //----------------------------------------------------------------------------
    1039 
    1040 void BookHistos( int verbosityLevel )
    1041 {
    1042         if (verbosityLevel > 2) cout << endl << "...book histograms" << endl;
    1043     hAllPixelList = new TObjArray;
    1044     hRootList = new TObjArray;
    1045     debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ];
    1046     for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){
    1047         debugHistos[ type ].SetBins(1024, 0, 1024);
    1048         debugHistos[ type ].SetLineColor(1);
    1049         debugHistos[ type ].SetLineWidth(2);
    1050         debugHistos[ type ].SetStats(false);
    1051 
    1052         // set X axis paras
    1053         debugHistos[ type ].GetXaxis()->SetLabelSize(0.08);
    1054         debugHistos[ type ].GetXaxis()->SetTitleSize(0.08);
    1055         debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0);
    1056         debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.));
    1057 
    1058         // set Y axis paras
    1059         debugHistos[ type ].GetYaxis()->SetLabelSize(0.08);
    1060         debugHistos[ type ].GetYaxis()->SetTitleSize(0.08);
    1061         debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
    1062         debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]");
    1063     }
    1064 
    1065     // All Pixel Histograms
    1066     //------------------------------------------------------------------------
    1067     TString histoname;
    1068     for (int order = 0; order < MAX_PULS_ORDER; order ++)
    1069     {
    1070 //        SetHistogramAttributes(
    1071 //                    "AllPixelPeakOverlay",
    1072 //                    order,
    1073 //                    gOverlayAttrib,
    1074 //                    MaxAmplOfFirstPulse
    1075 //                    );
    1076         histoname = "hAllPixelOverlay";
    1077         histoname += order;
    1078         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1079         hAllPixelOverlay[order] = new TH2F(
    1080                     histoname,
    1081                     "Overlay of detected Pulses of one Order for All Pixels",
    1082                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    1083                     (-1*gPixelOverlayXaxisLeft)-0.5,
    1084                     gPixelOverlayXaxisRight-0.5 ,
    1085                     512,
    1086                     -55.5,
    1087                     300.5
    1088                     );
    1089         hAllPixelOverlay[order]->SetAxisRange(
    1090                     gBSLMean - 5,
    1091                     (gGainMean*(order+1)) + 10,
    1092                     "Y");
    1093         hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1094         hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1095         //hAllPixelOverlay->SetBit(TH2F::kCanRebin);
    1096         hAllPixelList->Add( hAllPixelOverlay[order] );
    1097 
    1098     //------------------------------------------------------------------------
    1099 
    1100 //        SetHistogramAttributes(
    1101 //                    "AllPixelPeakMax",
    1102 //                    order,
    1103 //                    gMaximumAttrib,
    1104 //                    MaxAmplOfFirstPulse
    1105 //                    );
    1106         histoname = "hAllPixelMax";
    1107         histoname += order;
    1108         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1109         hAllPixelMax[order] = new TH1F(
    1110                     histoname,
    1111                     "Maximum value of each slice in overlay plot for All Pixels",
    1112                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    1113                     (-1*gPixelOverlayXaxisLeft)-0.5,
    1114                     gPixelOverlayXaxisRight-0.5
    1115                     );
    1116         hAllPixelMax[order]->SetAxisRange(
    1117                     gBSLMean - 5,
    1118                     (gGainMean*(order+1)) + 10,
    1119                     "Y");
    1120         hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1121         hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1122         hAllPixelList->Add( hAllPixelMax[order] );
    1123 
    1124 //------------------------------------------------------------------------
    1125 //        SetHistogramAttributes(
    1126 //                    "PeakMaximum",
    1127 //                    order,
    1128 //                    gMaximumAttrib,
    1129 //                    MaxAmplOfFirstPulse
    1130 //                    );
    1131         histoname = "hAllPixelMedian";
    1132         histoname += order;
    1133         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1134         hAllPixelMedian[order] = new TH1F(
    1135                     histoname,
    1136                     "Median of each slice in overlay plot for All Pixels",
    1137                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    1138                     (-1*gPixelOverlayXaxisLeft)-0.5,
    1139                     gPixelOverlayXaxisRight-0.5
    1140                     );
    1141         hAllPixelMedian[order]->SetAxisRange(
    1142                     gBSLMean - 5,
    1143                     (gGainMean*(order+1)) + 10,
    1144                     "Y");
    1145         hAllPixelMedian[order]->SetLineColor(kRed);
    1146         hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1147         hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1148         hAllPixelList->Add( hAllPixelMedian[order] );
    1149         if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl;
    1150 //------------------------------------------------------------------------
    1151 //        SetHistogramAttributes(
    1152 //                    "PeakMaximum",
    1153 //                    order,
    1154 //                    gMaximumAttrib,
    1155 //                    MaxAmplOfFirstPulse
    1156 //                    );
    1157         histoname = "hAllPixelMean";
    1158         histoname += order;
    1159         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1160         hAllPixelMean[order] = new TH1F(
    1161                     histoname,
    1162                     "Mean of each slice in overlay plot for All Pixels",
    1163                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    1164                     (-1*gPixelOverlayXaxisLeft)-0.5,
    1165                     gPixelOverlayXaxisRight-0.5
    1166                     );
    1167         hAllPixelMean[order]->SetAxisRange(
    1168                     gBSLMean - 5,
    1169                     (gGainMean*(order+1)) + 10,
    1170                     "Y");
    1171         hAllPixelMean[order]->SetLineColor(kBlue);
    1172         hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1173         hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1174         hAllPixelList->Add( hAllPixelMean[order] );
    1175 
    1176 //------------------------------------------------------------------------
    1177 //        SetHistogramAttributes(
    1178 //                    "PeakMaxGaus",
    1179 //                    order,
    1180 //                    gMaxGausAttrib,
    1181 //                    MaxAmplOfFirstPulse
    1182 //                    );
    1183         histoname = "hAllPixelEdgeOverlay";
    1184         histoname += order;
    1185         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1186         hAllPixelEdgeOverlay[order] = new TH2F(
    1187                     histoname,
    1188                     "Overlay at rising edge of detected pulses of one pulse order for All Pixels ",
    1189                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    1190                     (-1*gPixelOverlayXaxisLeft)-0.5,
    1191                     gPixelOverlayXaxisRight-0.5 ,
    1192                     512,
    1193                     -55.5,
    1194                     200.5
    1195                     );
    1196 
    1197         hAllPixelEdgeOverlay[order]->SetAxisRange(
    1198                     gBSLMean - 5,
    1199                     (gGainMean*(order+1)) + 10,
    1200                     "Y");
    1201         hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1202         hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1203         hAllPixelList->Add( hAllPixelEdgeOverlay[order] );
    1204 
    1205 //------------------------------------------------------------------------
    1206 //        SetHistogramAttributes(
    1207 //                    "AllPixelPeakProfile",
    1208 //                    order,
    1209 //                    gProfileAttrib,
    1210 //                    MaxAmplOfFirstPulse
    1211 //                    );
    1212         histoname = "hAllPixelProfile";
    1213         histoname += order;
    1214         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1215         hAllPixelProfile[order] = new TProfile(
    1216                     histoname,
    1217                     "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
    1218                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
    1219                     (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
    1220                     gPixelOverlayXaxisRight-0.5 ,                    //xup
    1221     //                512,                                            //nbinsy
    1222                     -55.5,                                          //ylow
    1223                     300.5,                                          //yup
    1224                     "s");                                           //option
    1225         hAllPixelProfile[order]->SetAxisRange(
    1226                     gBSLMean - 5,
    1227                     (gGainMean*(order+1)) + 10,
    1228                     "Y");
    1229         hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1230         hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1231         //hPixelProfile->SetBit(TH2F::kCanRebin);
    1232         hAllPixelList->Add( hAllPixelProfile[order] );
    1233 
    1234 //------------------------------------------------------------------------
    1235 //        SetHistogramAttributes(
    1236 //                    "AllPixelPeakProfile",
    1237 //                    order,
    1238 //                    gProfileAttrib,
    1239 //                    MaxAmplOfFirstPulse
    1240 //                    );
    1241         histoname = "hAllPixelProfile2";
    1242         histoname += order;
    1243         if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    1244         hAllPixelProfile2[order] = new TProfile(
    1245                     histoname,
    1246                     "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
    1247                     gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
    1248                     (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
    1249                     gPixelOverlayXaxisRight-0.5 ,                    //xup
    1250     //                512,                                            //nbinsy
    1251                     -55.5,                                          //ylow
    1252                     300.5,                                          //yup
    1253                     "s");                                           //option
    1254         hAllPixelProfile2[order]->SetAxisRange(
    1255                     gBSLMean - 5,
    1256                     (gGainMean*(order+1)) + 10,
    1257                     "Y");
    1258         hAllPixelProfile2[order]->SetLineColor(kRed);
    1259         hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    1260         hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    1261         //hPixelProfile->SetBit(TH2F::kCanRebin);
    1262         hAllPixelList->Add( hAllPixelProfile2[order] );
    1263 
    1264 
    1265 
    1266     }//end of for over order
    1267     if (verbosityLevel > 2) cout << "...done" << endl;
    1268 }
    1269 // end of BookHistos
    1270 //----------------------------------------------------------------------------
    1271 
    1272 
    1273 void CreateRootFile(const char * loc_fname, int verbosityLevel)
    1274 {
    1275     TFile* tf = new TFile(loc_fname,"UPDATE");
    1276     if (tf->IsZombie())
    1277     {
    1278         cout << "Error opening file" << endl;
    1279         exit(-1);
    1280     }
    1281     else {
    1282         if (verbosityLevel > 1) cout << "creating root-file successfull" << endl;
    1283         tf->Close();
    1284     }
    1285 }
    1286 //end of CreateRootFile
    1287 //----------------------------------------------------------------------------
    1288 
    1289 
    1290 TFile* OpenRootFile(const char * loc_fname, int verbosityLevel)
    1291 {
    1292     TFile* tf = new TFile(loc_fname,"UPDATE");
    1293     if (tf->IsZombie())
    1294     {
    1295         cout << "Error opening file" << endl;
    1296         exit(-1);
    1297     }
    1298     else {
    1299         if (verbosityLevel > 1) cout << "...opening root-file:"
    1300                                      << loc_fname
    1301                                      << " successfull" << endl;
    1302     }
    1303     return tf;
    1304 }
    1305 //end of OpenRootFile
    1306 //----------------------------------------------------------------------------
    1307 
    1308 
    1309 void CloseRootFile(TFile* tf)
    1310 {
    1311     if (tf->IsZombie())
    1312     {
    1313     cout << "Error closing file" << endl;
    1314     exit(-1);
    1315     } else {
    1316     tf->Close();
    1317     }
    1318 }
    1319 //end of CloseRootFile
    1320 //----------------------------------------------------------------------------
    1321 
    1322 
    1323 void SaveHistograms(
    1324         const char* loc_fname,
    1325         const char* subdirectory,
    1326         TObjArray* histList,
     433DeletePixelCanvases(
     434        int maxPulseOrder,
    1327435        int verbosityLevel
    1328436        )
    1329437{
    1330     if (verbosityLevel > 2) cout << endl
    1331                                  << "...saving pixel histograms to subdirectory: "
    1332                                  << subdirectory << endl ;
    1333 
    1334     TFile* tf = OpenRootFile(loc_fname, verbosityLevel);
    1335     if ( subdirectory != "root")
    1336     {
    1337         tf->cd();
    1338         tf->mkdir(subdirectory,subdirectory);
    1339         tf->cd(subdirectory);
    1340     }
    1341     else
    1342     {
    1343         tf->cd();
    1344     }
    1345     histList->Write(); // write the major histograms into the top level directory
    1346     if (verbosityLevel > 3) tf->ls();
    1347     CloseRootFile( tf ); // close the file
    1348     if (verbosityLevel > 2) cout << "...done" << endl;
     438    if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ;
     439    for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ )
     440    {
     441        delete cgpPixelPulses[pulse_order];
     442        delete cgpDistributions[pulse_order];
     443    }
     444    delete[] cgpPixelPulses;
     445    delete[] cgpDistributions;
    1349446}
    1350 // end of SaveHistograms
    1351 //----------------------------------------------------------------------------
    1352 
    1353 
    1354 TString CreateSubDirName(int pixel)
    1355 {
    1356     TString path = "Pixel_";
    1357     path += pixel;
    1358 //    cout << "path " << path << endl ;
    1359     return path;
    1360 }
    1361 // end of CreateSubDirName
    1362 //----------------------------------------------------------------------------
    1363 
    1364 TString CreateSubDirName(const char* title)
    1365 {
    1366     TString path = title;
    1367     path += "_Pixel";
    1368 //    cout << "path " << path << endl ;
    1369     return path;
    1370 }
    1371 // end of CreateSubDirName
    1372 //----------------------------------------------------------------------------
    1373 
    1374 void FillHistograms(
    1375         vector<Region>* pZXings,
    1376         int AmplWindowWidth,
    1377         int eventNumber,
    1378         int verbosityLevel
    1379         )
    1380 {
    1381     if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ;
    1382     bool use_this_peak=false;
    1383     int order_of_pulse=0;
    1384     vector<Region>::iterator reg;
    1385 
    1386     //Loop over all found zerocrossings in Region
    1387     for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
    1388     {
    1389         //skip those who are at the Rim of the ROI
    1390         if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12)
    1391         {
    1392             if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl;
    1393             continue;
    1394         }
    1395         //domstest->Fill(reg->maxPos);
    1396 
    1397         // Define axis range of Histogram
    1398         int Left = reg->maxPos - gPixelOverlayXaxisLeft;
    1399         int Right = reg->maxPos + gPixelOverlayXaxisRight;
    1400         int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft;
    1401         int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight;
    1402 
    1403         if (Left < 0) Left =0;
    1404         if (Right > (int)Ameas.size() ) Right=Ameas.size() ;
    1405         if (EdgeLeft < 0) EdgeLeft =0;
    1406         if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ;
    1407 
    1408 
    1409         hTesthisto->Fill( reg->slopeOfRisingEdge ) ;
    1410         hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ;
    1411 
    1412         if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
    1413 
    1414         //determine order of pulse and which histogram shall be filled
    1415         if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
    1416         for(int order = 0; order < MAX_PULS_ORDER; order++)
    1417         {
    1418             if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2))
    1419                 && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2))
    1420             {
    1421                 //hOverlayTemp = hPixelOverlay[order];
    1422                 if (verbosityLevel > 2) cout << "...#" << order ;
    1423                 order_of_pulse = order;
    1424                 use_this_peak = true;
    1425                 break;
    1426             }
    1427          else if (order >= (MAX_PULS_ORDER - 1)){
    1428                 use_this_peak = false;
    1429                 if (verbosityLevel > 2) cout << "...NONE" << endl ;
    1430          }
    1431         }
    1432 
    1433         //Fill overlay und profile histograms
    1434         if (use_this_peak){
    1435             if (verbosityLevel > 2) cout << "...filling" ;
    1436             for ( int pos = Left; pos < Right; pos++){
    1437 //                if ();
    1438                 hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    1439                 hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    1440                 hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    1441                 hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    1442             }
    1443             for ( int pos = EdgeLeft; pos < EdgeRight; pos++){
    1444 //                cout << endl << "###filling edge histo###" << endl;
    1445                 hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
    1446                 hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
    1447                 hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
    1448                 hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
    1449 //                cout << endl << "######" << endl;
    1450             }
    1451             if (verbosityLevel > 2) cout << "...done" << endl;
    1452         }
    1453         //Breakout option
    1454         if (breakout)
    1455         break;
    1456    }
    1457     if (verbosityLevel > 2) cout << "...done" << endl;
    1458 }
    1459 // end of FillHistograms
    1460 //----------------------------------------------------------------------------
    1461 
    1462 void DrawTestHistograms(
    1463         int verbosityLevel
    1464         )
    1465 {
    1466     if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;
    1467 
    1468     cgpTestHistos->cd(1);
    1469     hTesthisto->Draw();
    1470     cgpTestHistos->cd(2);
    1471     hTesthisto2->Draw("BOX");
    1472 }
    1473 // end of DrawTestHistograms
    1474 //----------------------------------------------------------------------------
    1475 
    1476 void DrawPulseHistograms(
     447
     448void
     449UpdateCanvases(
    1477450        int verbosityLevel,
    1478         int max_pulse_order
    1479         )
    1480 {
    1481     if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;
    1482     for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
    1483     {
    1484         cgpPixelPulses[pulse_order]->cd(1);
    1485         hPixelOverlay[pulse_order]->Draw("COLZ");
    1486         cgpPixelPulses[pulse_order]->cd(3);
    1487         hPixelProfile[pulse_order]->Draw("COLZ");
    1488         hPixelProfile2[pulse_order]->Draw("SAME");
    1489 
    1490         cgpPixelPulses[pulse_order]->cd(2);
    1491         hPixelEdgeOverlay[pulse_order]->Draw("COLZ");
    1492 //        cgpPixelPulses[pulse_order]->cd(4);
    1493 //        hTesthisto->Draw();
    1494 
    1495         cgpAllPixelPulses[pulse_order]->cd(1);
    1496         hAllPixelOverlay[pulse_order]->Draw("COLZ");
    1497         cgpAllPixelPulses[pulse_order]->cd(2);
    1498         hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ");
    1499 
    1500         cgpAllPixelPulses[pulse_order]->cd(3);
    1501         hAllPixelProfile[pulse_order]->Draw("COLZ");
    1502         hAllPixelProfile2[pulse_order]->Draw("SAME");
    1503 
    1504     }
    1505 }
    1506 // end of DrawPulseHistograms
    1507 //----------------------------------------------------------------------------
    1508 
    1509 
    1510 void DrawMaximumHistograms(
    1511         int verbosityLevel,
    1512         int max_pulse_order
    1513         )
    1514 {
    1515     if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ;
    1516     for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
    1517     {
    1518         cgpPixelPulses[pulse_order]->cd(4);
    1519         hPixelMax[pulse_order]->Draw();
    1520         hPixelMedian[pulse_order]->Draw("SAME");
    1521         hPixelMean[pulse_order]->Draw("SAME");
    1522 
    1523         cgpAllPixelPulses[pulse_order]->cd(4);
    1524         hAllPixelMax[pulse_order]->Draw();
    1525         hAllPixelMedian[pulse_order]->Draw("SAME");
    1526         hAllPixelMean[pulse_order]->Draw("SAME");
    1527 //        cgpAllPixelPulses[pulse_order]->cd(3);
    1528      //   hAllPixelMaxGaus[pulse_order]->Draw("COLZ");
    1529     }
    1530 }
    1531 // end of DrawMaximumHistograms
    1532 //----------------------------------------------------------------------------
    1533 
    1534 
    1535 void UpdateCanvases(
    1536         int verbosityLevel,
    1537         int max_pulse_order
     451        int max_pulse_order,
     452        bool testmode
    1538453        )
    1539454{
    … …  
    1543458        cgpPixelPulses[pulse_order]->Modified();
    1544459        cgpPixelPulses[pulse_order]->Update();
    1545         cgpAllPixelPulses[pulse_order]->Modified();
    1546         cgpAllPixelPulses[pulse_order]->Update();
    1547         cgpTestHistos->Modified();
    1548         cgpTestHistos->Update();
     460        cgpDistributions[pulse_order]->Modified();
     461        cgpDistributions[pulse_order]->Update();
     462//        if ( testmode )
     463//        {
     464//            cgpTestHistos->Modified();
     465//            cgpTestHistos->Update();
     466//        }
    1549467    }
    1550468}
    1551 // end of UpdateCanvases
    1552 //----------------------------------------------------------------------------
    1553 
    1554 
    1555 
    1556 void
    1557 PlotMedianEachSliceOfPulse(
    1558         unsigned int    max_pulse_order,
    1559         int             fitdata,
    1560         int             verbosityLevel
    1561         )
     469//-----------------------------------------------------------------------------
     470// MAIN Funktion for C compilation
     471//-----------------------------------------------------------------------------
     472
     473int main()
    1562474{
    1563     if (verbosityLevel > 2) cout << endl
    1564                                  << "...calculating pulse shape of slice's Median"
    1565                                  << endl;
    1566     for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    1567     {
    1568         if (verbosityLevel > 2) cout << "\t...calculation of "
    1569                                      << "pulse order " << pulse_order;
    1570 
    1571         Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
    1572 
    1573         for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
    1574            TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
    1575            float median = MedianOfH1(h1);
    1576            float mean = h1->GetMean();
    1577            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
    1578            delete h1;
    1579            hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
    1580            hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
    1581            hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
    1582            hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
    1583         }
    1584 
    1585         if (verbosityLevel > 2) cout << "\t...done" << endl;
    1586 
    1587         if (fitdata)
    1588         {
    1589            FitMaxPropabilityPuls(
    1590                     hPixelMedian[pulse_order],
    1591                     verbosityLevel
    1592                     );
    1593         }
    1594     }
     475
     476FCalcPulseTemplate();
     477return 0;
     478
    1595479}
    1596 // end of PlotMedianEachSliceOfPulse
    1597 //----------------------------------------------------------------------------
    1598 
    1599 Double_t MedianOfH1 (
    1600         TH1*            h1
    1601         )
    1602 {
    1603    //compute the median for 1-d histogram h1
    1604    Int_t nbins = h1->GetXaxis()->GetNbins();
    1605    Double_t *x = new Double_t[nbins];
    1606    Double_t *y = new Double_t[nbins];
    1607    for (Int_t i=0;i<nbins;i++) {
    1608       x[i] = h1->GetXaxis()->GetBinCenter(i+1);
    1609       y[i] = h1->GetBinContent(i+1);
    1610    }
    1611    Double_t median = TMath::Median(nbins,x,y);
    1612    delete [] x;
    1613    delete [] y;
    1614    return median;
    1615 }
    1616 // end of PlotMedianEachSliceOfPulse
    1617 //----------------------------------------------------------------------------
    1618 
    1619 void
    1620 PlotPulseShapeOfMaxPropability(
    1621         unsigned int    max_pulse_order,
    1622         int             fitdata,
    1623         int             verbosityLevel
    1624         )
    1625 {
    1626     if (verbosityLevel > 2) cout << endl
    1627                                  << "...calculating pulse shape of max. propability"
    1628                                  << endl;
    1629     for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    1630     {
    1631         if (verbosityLevel > 2) cout << "\t...calculation of "
    1632                                      << "pulse order " << pulse_order;
    1633         //  vector max_value_of to number of timeslices in Overlay Spectra
    1634         float max_value_of_slice;
    1635 
    1636         Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
    1637         if (verbosityLevel > 2) cout << "...generating projections" << endl;
    1638         for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    1639         {
    1640             if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
    1641             //put maximumvalue of every Y-projection of every timeslice into vector
    1642             TH1D *h1 =  hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice);
    1643 
    1644             max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() );
    1645 
    1646             if (verbosityLevel > 2) cout << " with max value "
    1647                                          << max_value_of_slice << endl;
    1648             hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
    1649             hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
    1650             delete h1;
    1651         }
    1652 
    1653         if (verbosityLevel > 2) cout << "\t...done" << endl;
    1654 
    1655         if (fitdata)
    1656         {
    1657            FitMaxPropabilityPuls(
    1658                     hPixelMax[pulse_order],
    1659                     verbosityLevel
    1660                     );
    1661         }
    1662     }
    1663 }
    1664 
    1665 
    1666 void
    1667 FitMaxPropabilityPuls(
    1668         TH1F*       hMaximumTemp,
    1669         int         verbosityLevel
    1670         )
    1671     {
    1672       if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
    1673       if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
    1674       hMaximumTemp->Fit("landau", "", "", -50, 250);
    1675       if (verbosityLevel > 2) cout << "...done" << endl;
    1676     }
    1677 
    1678 
    1679 //void SetHistogrammSettings(
    1680 //        const char*     histoType,
    1681 //        int             max_puls_order,
    1682 //        int             MaxAmplOfFirstPulse,
    1683 //        int             verbosityLevel
    1684 //        )
    1685 //{
    1686 //    if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl;
    1687 //    for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++)
    1688 //    {
    1689 //        qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order);
    1690 //        cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl;
    1691 //        gHistSetting[pulse_order].name = "h";
    1692 //        gHistSetting[pulse_order].name += histoType;
    1693 //        gHistSetting[pulse_order].name += "PeakOverlay";
    1694 //        gHistSetting[pulse_order].name += pulse_order;
    1695 //        gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around ";
    1696 //        gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax;
    1697 //        gHistSetting[pulse_order].title += " mV";
    1698 //        gHistSetting[pulse_order].Mname = "h";
    1699 //        gHistSetting[pulse_order].Mname += histoType;
    1700 //        gHistSetting[pulse_order].Mname += "PeakMaximum";
    1701 //        gHistSetting[pulse_order].Mname += pulse_order;
    1702 //        gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";
    1703 //        gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax;
    1704 //        gHistSetting[pulse_order].Mtitle += " mV";
    1705 //    }
    1706 //    if (verbosityLevel > 2) cout << "...done" << endl;
    1707 //}
    1708 
    1709 
    1710 //void CreatePulseCanvas(
    1711 //        TCanvas*    array_of_canvases,
    1712 //        unsigned int        max_pulse_order,
    1713 //        const char*      canvas_name,
    1714 //        const char*      canvas_title,
    1715 //        const char*      pixel,
    1716 //        int verbosityLevel
    1717 //        )
    1718 //{
    1719 //    if (verbosityLevel > 2) cout << "...building Canvases" ;
    1720 
    1721 //    for (
    1722 //         unsigned int pulse_order = 1;
    1723 //         pulse_order <= max_pulse_order;
    1724 //         pulse_order++
    1725 //          )
    1726 //    {
    1727 //        TString name = canvas_name;
    1728 //                name += " ";
    1729 //                name += pulse_order;
    1730 //                name += " ";
    1731 
    1732 //                TString title = "Order ";
    1733 //                title += pulse_order;
    1734 //                title += " ";
    1735 //                title += canvas_title;
    1736 //                title += " of ";
    1737 //                title += pixel;
    1738 //                title += ". Pixel";
    1739 
    1740 //        array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700);
    1741 //        array_of_canvases[pulse_order].Divide(2,2);
    1742 //    }
    1743 //    if (verbosityLevel > 2) cout << "...done" << endl;
    1744 //}
    1745 
    1746 
    1747 //void
    1748 //DeletePulseCanvas(
    1749 //        TCanvas* array_of_canvases,
    1750 //        unsigned int     max_pulse_order
    1751 //        )
    1752 //{
    1753 //    for (
    1754 //         unsigned int pulse_order = 1;
    1755 //         pulse_order <= max_pulse_order;
    1756 //         pulse_order++
    1757 //          )
    1758 //    {
    1759 //        delete array_of_canvases[pulse_order];
    1760 //    }
    1761 //}
    1762 
    1763 
    1764 
    1765 
    1766 void
    1767 DeletePixelCanvases( int verbosityLevel )
    1768 {
    1769     if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ;
    1770     for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++)
    1771     {
    1772         delete cgpPixelPulses[pulse_order];
    1773     }
    1774 }
    1775 
    1776 //----------------------------------------------------------------------------
    1777 //----------------------------------------------------------------------------
    1778 
    1779 //void SetHistogramAttributes(
    1780 //        const char*     histoType,
    1781 //        int             pulse_order,
    1782 //        histAttrib_t*   histoAttrib,
    1783 //        int             max_ampl_of_first_pulse
    1784 //        )
    1785 //{
    1786 //    histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order);
    1787 //    histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2);
    1788 //    histoAttrib[pulse_order].name = "h";
    1789 //    histoAttrib[pulse_order].name += histoType;
    1790 //    histoAttrib[pulse_order].name = "_";
    1791 //    histoAttrib[pulse_order].name += pulse_order;
    1792 
    1793 //    histoAttrib[pulse_order].title += histoType;
    1794 //    histoAttrib[pulse_order].title += " with Amplitude around ";
    1795 //    histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax;
    1796 //    histoAttrib[pulse_order].title += " mV";
    1797 
    1798 //    histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]";
    1799 //    histoAttrib[pulse_order].yTitle += "Amplitude [mV]";
    1800 //}
    1801 
    1802 void WritePixelTemplateToCsv(
    1803         TString         path,
    1804         const char*     csv_file_name,
    1805         const char*     overlay_method,
    1806         int             pixel,
    1807         int             verbosityLevel
    1808         )
    1809 {
    1810     path += csv_file_name;
    1811     path += "_";
    1812     path += pixel;
    1813     path += ".csv";
    1814 
    1815     Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins();
    1816     if (verbosityLevel > 0)
    1817     {
    1818         cout << "writing point-set to csv file: " ;
    1819         cout << path << endl;
    1820         cout << "...opening file" << endl;
    1821     }
    1822     if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
    1823     ofstream out;
    1824     out.open( path );
    1825     out << "### point-set of a single photon pulse template" << endl;
    1826     out << "### template determined with pulse overlay at: "
    1827         << overlay_method << endl;
    1828     out << "### Slice's Amplitude determined by calculating the " << endl
    1829         << "### value of maximum propability of slice -> Amplitude1 " << endl
    1830         << "### mean of slice -> Amplitude2 " << endl
    1831         << "### median of slice -> Amplitude3 " << endl
    1832         << "### for each slice" << endl;
    1833     out << "### Pixel number (CHid): " << pixel << endl
    1834         << endl;
    1835 
    1836     out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
    1837 
    1838     for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    1839     {
    1840         out << TimeSlice << "," ;
    1841         out << hPixelMax[0]->GetBinContent(TimeSlice) << ",";
    1842         out << hPixelMean[0]->GetBinContent(TimeSlice) << ",";
    1843         out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl;
    1844     }
    1845     out.close();
    1846     if (verbosityLevel > 0) cout << "...file closed" << endl;
    1847 }
    1848 
    1849 void WriteAllPixelTemplateToCsv(
    1850         TString         path,
    1851         const char*     csv_file_name,
    1852         const char*     overlay_method,
    1853         int             verbosityLevel
    1854         )
    1855 {
    1856     path += csv_file_name;
    1857     path += "_AllPixel.csv";
    1858 
    1859     Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins();
    1860     if (verbosityLevel > 0)
    1861     {
    1862         cout << "writing point-set to csv file: " ;
    1863         cout << path << endl;
    1864         cout << "...opening file" << endl;
    1865     }
    1866     if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
    1867     ofstream out;
    1868     out.open( path );
    1869     out << "### point-set of a single photon pulse template" << endl;
    1870     out << "### template determined with pulse overlay at: "
    1871         << overlay_method << "of all Pixels";
    1872     out << "### Slice's Amplitude determined by calculating the " << endl
    1873         << "### value of maximum propability of slice -> Amplitude1 " << endl
    1874         << "### mean of slice -> Amplitude2 " << endl
    1875         << "### median of slice -> Amplitude3 " << endl
    1876         << "### for each slice" << endl
    1877         << endl;
    1878 
    1879     out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
    1880 
    1881     for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    1882     {
    1883         out << TimeSlice << "," ;
    1884         out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ",";
    1885         out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ",";
    1886         out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl;
    1887     }
    1888     out.close();
    1889     if (verbosityLevel > 0) cout << "...file closed" << endl;
    1890 }
Note: See TracChangeset for help on using the changeset viewer.