#include <stdio.h>
#include <stdlib.h>
#include "midas.h"
#include "../../compress/mucap_compress.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"

Functions
void	process_file (char *input_filename, int addr, int fadc)

int	process_event (void *pevent, int addr, int fadc)

bool	get_input_event (FILE fp, char pevent, int event_number)

double	findPulse (int nsamples, int *samples)

int	processIsland (int t0, int nsamples, int *samples, int channel)

Variables
TCanvas *	c = 0

TH2 *	hSamples

TH1 *	hChannelNSample

TH1 *	hChannelNIsland

TH2 *	hMinSample

TH2 *	hMaxSample

TH2 *	hPulseEnergy

TH2 *	hPed0

TH1 *	hPulseIF

TH2 *	hPulseIvsEnergy

TH2 *	hPulseFvsEnergy

TH2 *	hPulseAvsEnergy

Function Documentation

double findPulse	(	int	nsamples,
		int *	samples
	)

Definition at line 122 of file res_analysis_new.cxx.

References hPulseAvsEnergy, hPulseFvsEnergy, hPulseIF, hPulseIvsEnergy, i, nsamples, and pedestalsigma.

                                             {
   int sig_length_counter = 0;
   double sig_max_height = 0.;
   int sig_first_index = 0;
   double sig_integral = 0.;
   int sig_over_thr = 0;
   
   bool hasOvershoot = false;
   bool hasUndershoot = false;
   const int kLengthThreshold = 4;
 
   // (Serdar) i increased nsamples to 40 (was 16)
 
   if(nsamples < 40 ) return 0.;
 
   // Let's first find the pedestal 
 
   // (Serdar) here i check whether selecting the first 5 sample or 10 sample do change anything
   // on the pedestal determination, see histograms histPeds0first05bins or 10bins for more info
   // and i find both can be used, so i select first 10 bins to get more info.
 
   const int pedSamples = 10;
   double pedestalOffset = 0;
   double sumpedestal2 = 0;
   double pedestalsigma = 0 ;
 
   if(pedSamples > nsamples) return 0.;
 
   for(int i = 0; i < pedSamples; i++) {
     pedestalOffset += samples[i];
   }
   pedestalOffset /= pedSamples;
 
   for(int i = 0; i < pedSamples; i++) {
     sumpedestal2 += (pedestalOffset-samples[i])*(pedestalOffset-samples[i]);
   }
   pedestalsigma= sqrt(sumpedestal2/pedSamples);
 
   // printf("Pedestal = %f\n", pedestalOffset);
   // Pedestal found, let's look for a pulse now
 
   // (Serdar) here i calculated the sigma for pedestal calculation, then i will start
   // searching the pulse if sample is bigger the 4 pedestal sigma,
 
   const int kThreshold = 4*pedestalsigma;
 
   for(int i = 0; i < nsamples; i++) {
     double bin_content = -1.*((double)samples[i] - pedestalOffset);
     
     double next_bin_content = 0.;
     if(i<nsamples-1) next_bin_content = -1.*((double)samples[i+1] - pedestalOffset);
 
     if(bin_content >= kThreshold || next_bin_content >= kThreshold){
       sig_length_counter++;
       sig_over_thr++;
       
       if (sig_length_counter == 1){
         sig_first_index = i;
         // inspect left shoulder of this pulse
         int back;
         for(back=i-1; back>=0; back--){
           double cont = -1.*((double)samples[back] - pedestalOffset);
 
           // (Serdar) here i made 2 changes, first, instead of 1/4 of the kThreshold
           // i have 3/4, which equals 3 sigma at the left shoulder
 
           if(cont > 3*kThreshold/4 && (i-back)<4){ 
             // accept it still as part of the left shoulder
             sig_length_counter++;
             sig_first_index = back;             
           }
           else break;
         }
       }
     }
     else{
       if (((bin_content < kThreshold) || (i == nsamples - 1)) && 
           (sig_over_thr >= kLengthThreshold)){
         // It's the end of the pulse, so let's include part of the right 
         // shoulder
         int forw;  // index for samples right of momentary position
         for(forw=i+1;forw<nsamples; forw++){
           double cont = -1.*((double)samples[forw] - pedestalOffset);
           if(cont > 3*kThreshold/4 && (forw-i)<4){
             // accept it still as part of the left shoulder
             sig_length_counter++;
           }
           else break;  // probably pedestal so let's stop
           i = forw;
         }
         
         sig_integral = 0.;
         sig_max_height = 0.;
         
         for(int i = sig_first_index; i < sig_first_index+sig_length_counter; i++){
           double bin_content = -1.*((double)samples[i] - pedestalOffset);
           sig_integral += bin_content;
           if (bin_content > sig_max_height){
             sig_max_height = bin_content;
           }
         }
         
         // Let's only deal with one pulse at this time!!!
         break;
       }
     }
   }
 
   // printf("Pulse accepted from sample %d to %d with energy %f\n", sig_first_index,
   //     sig_first_index+sig_length_counter, sig_integral);
   hPulseIF->Fill(sig_first_index);
   hPulseIF->Fill(sig_first_index+sig_length_counter);
   hPulseAvsEnergy->Fill(sig_integral, sig_max_height);
   hPulseIvsEnergy->Fill(sig_integral, sig_first_index);
   hPulseFvsEnergy->Fill(sig_integral, sig_first_index+sig_length_counter);
   return sig_integral;
 }

bool get_input_event	(	FILE *	fp,
		char *	pevent,
		int	event_number
	)

Definition at line 27 of file res_analysis_new.cxx.

References SUCCESS.

 {
     int event_size, size_read;
     EVENT_HEADER *header;
 
     size_read = fread(pevent, 1, sizeof(EVENT_HEADER), fp);
     if (size_read != sizeof(EVENT_HEADER))
         return !SUCCESS;
 
     header = (EVENT_HEADER *) pevent;
 
     if (event_number == 0 && header->trigger_mask != MIDAS_MAGIC) {
         fprintf(stderr, "Input is not a valid MIDAS file\n");
         exit(1);
     }
 
     event_size = header->data_size;
 
     if ((unsigned int) event_size >= MAX_EVENT_SIZE - sizeof(EVENT_HEADER)) {
         fprintf(stderr, "Invalid event size %d\n", event_size);
         exit(1);
     }
 
     size_read = fread(pevent + sizeof(EVENT_HEADER), 1, event_size, fp);
     if (size_read != event_size)
         return !SUCCESS;
 
     return SUCCESS;
 }

int process_event	(	void *	pevent,
		int	addr,
		int	fadc
	)

void process_file	(	char *	input_filename,
		int	addr,
		int	fadc
	)

int processIsland	(	int	t0,
		int	nsamples,
		int *	samples,
		int	channel
	)

Definition at line 240 of file res_analysis_new.cxx.

References findPulse(), hChannelNIsland, hChannelNSample, hMaxSample, hMinSample, hPed0, hPulseEnergy, hSamples, i, and nsamples.

 {
   hChannelNIsland->Fill(channel);
 
   if(nsamples < 40) return 0;
 
   int minI = -1;
   int maxI = -1;
   int minSample = 4097;;
   int maxSample = -1;
 
   double energy = findPulse(nsamples, samples);
 
   for(int i = 0; i < nsamples; i++) {
     hSamples->Fill(channel, samples[i]);
     hChannelNSample->Fill(channel);
     if(samples[i] > maxSample) {
       maxI = i;
       maxSample = samples[i];
     }
     if(samples[i] < minSample) {
       minI = i;
       minSample = samples[i];
     }
     if(i<10){
       hPed0->Fill(channel, samples[i]);
     }
   }
 
   hMinSample->Fill(channel, minSample);
   hMaxSample->Fill(channel, maxSample);
   hPulseEnergy->Fill(channel, energy);
 
   return 0;
 }

Variable Documentation

TCanvas* c = 0

Definition at line 22 of file res_analysis_new.cxx.

TH1* hChannelNIsland

Definition at line 59 of file res_analysis_new.cxx.

TH1* hChannelNSample

Definition at line 58 of file res_analysis_new.cxx.

TH2* hMaxSample

Definition at line 61 of file res_analysis_new.cxx.

TH2* hMinSample

Definition at line 60 of file res_analysis_new.cxx.

TH2* hPed0

Definition at line 63 of file res_analysis_new.cxx.

TH2* hPulseAvsEnergy

Definition at line 67 of file res_analysis_new.cxx.

Referenced by findPulse().

TH2* hPulseEnergy

Definition at line 62 of file res_analysis_new.cxx.

TH2* hPulseFvsEnergy

Definition at line 66 of file res_analysis_new.cxx.

Referenced by findPulse().

TH1* hPulseIF

Definition at line 64 of file res_analysis_new.cxx.

Referenced by findPulse().

TH2* hPulseIvsEnergy

Definition at line 65 of file res_analysis_new.cxx.

Referenced by findPulse().

TH2* hSamples

Definition at line 57 of file res_analysis_new.cxx.

Functions

Variables

Function Documentation

Variable Documentation