MdEdxPlot.cpp File Reference

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <map>
#include <utility>
#include <sstream>
#include <cmath>
#include "midas.h"
#include <TH1.h>
#include <TH2.h>
#include "TGlobalData.h"
#include "TSetupData.h"
#include "TPulseIsland.h"

Go to the source code of this file.

Typedefs
typedef std::vector < std::vector< TPulseIsland * > ::iterator >	vecOfIterToPulses

Functions
INT	MdEdxPlot_init (void)
INT	MdEdxPlot (EVENT_HEADER *, void *)
std::map< double, double >	MakePulseTimeHeightMap (std::vector< TPulseIsland * > thePulses)
std::map< double, double >	MakePulseHeightSums (std::vector< std::map< double, double > > timeHeightMapsVector, double time_difference)

Variables
HNDLE	hDB
TGlobalData *	gData
TSetupData *	gSetup
static TH2D *	hdEdx_left
static TH2D *	hdEdx_right
std::vector< std::string >	left_thin_banknames
std::vector< std::string >	left_thick_banknames
std::vector< std::string >	right_thin_banknames
std::vector< std::string >	right_thick_banknames
ANA_MODULE	MdEdxPlot_module

Typedef Documentation

typedef std::vector<std::vector<TPulseIsland*>::iterator> vecOfIterToPulses

Definition at line 36 of file MdEdxPlot.cpp.

Function Documentation

std::map< double, double > MakePulseHeightSums	(	std::vector< std::map< double, double > >	timeHeightMapsVector,
		double	time_difference
	)

Definition at line 235 of file MdEdxPlot.cpp.

References time.

Referenced by MdEdxPlot().

                                                                                                                            {

  // Create the iterators
  std::vector<std::map<double, double>::iterator > pulseIters;
  std::vector<std::map<double, double>::iterator > backwardIters;

  for (int b = 0; b < timeHeightMapsVector.size(); b++) {
    pulseIters.push_back((timeHeightMapsVector.at(b)).begin()); 
    backwardIters.push_back((timeHeightMapsVector.at(b)).end());
  }

  // Prepare the output map:
  std::map<double, double> output;

  while ( pulseIters.size() > 0 ) {

    //Iterators/pointers into pulseIters and backwardIters must not be held past
    //one loop cycle as they may become invalidated.

    //Get the time of the next (unused) TPulseIslands (one from each bank)
    //work out which is the earliest
    double min_time = 99999999999.;  //something huge.
    for (int b = 0; b < pulseIters.size(); ++b){
      double time = (pulseIters.at(b))->first;
      min_time = std::min(min_time , time);
      //      printf("Bank #%d: time = %f ms\n", b, time * 1e-6);
      //Or use GetTimeStamp()?
    }
    //    printf("Min time = %f ms\n", min_time*1e-6);

    //Sum all the pulseHeights at a given time 
    double pulseHeight = 0;
    double pulseTime = 0;
    int nSummed = 0;
    for (int b = 0; b < pulseIters.size(); ++b){
      std::pair<double, double> timeHeightPair = *(pulseIters.at(b));
      //Need to define a comparison for CONSISTENT_TIME
      if ( std::abs(timeHeightPair.first - min_time) < time_difference ) {
        pulseHeight += timeHeightPair.second;
        pulseTime += timeHeightPair.first;
        ++nSummed;
        //and increment the iterator because we used the island
        ++(pulseIters.at(b));
        //      printf("Bank #%d: height = %f, time = %f ms\n", b, timeHeightPair.second, timeHeightPair.first * 1e-6);
      }
    } // for (int b)

    //append to output
    //    printf("Summed time = %f ms, Summed height = %f\n", (pulseTime/nSummed)*1e-6, pulseHeight);
    output[pulseTime/nSummed] = pulseHeight;


    //Delete the iterators to finished banks. Go through in reverse to
    //avoid invalidation problems
    for (int b = pulseIters.size()-1; b >= 0 ; --b){
      if (pulseIters.at(b) == backwardIters.at(b)){
        pulseIters.erase(pulseIters.begin() + b);
        backwardIters.erase(backwardIters.begin() + b);
      }  
    } // for (int b -reversed)
  } // while (size > 0)

  return output;
}

std::map< double, double > MakePulseTimeHeightMap

(

std::vector< TPulseIsland * >

thePulses)

Definition at line 219 of file MdEdxPlot.cpp.

Referenced by MdEdxPlot().

                                                                                  {

  std::map<double, double> output;

  for (std::vector<TPulseIsland*>::iterator pulseIter = thePulses.begin(); 
       pulseIter != thePulses.end(); pulseIter++) {

    output[(*pulseIter)->GetPulseTime()] = (*pulseIter)->GetPulseHeight();
  }

  return output;  
}

INT MdEdxPlot	(	EVENT_HEADER *	pheader,
		void *	pevent
	)

This method processes one MIDAS block, producing a vector of TOctalFADCIsland objects from the raw Octal FADC data.

Definition at line 115 of file MdEdxPlot.cpp.

References TGlobalData::fPulseIslandToChannelMap, hdEdx_left, hdEdx_right, MakePulseHeightSums(), MakePulseTimeHeightMap(), and SUCCESS.

{
        // Get the event number
        int midas_event_number = pheader->serial_number;

        // Some typedefs
        typedef map<string, vector<TPulseIsland*> > TStringPulseIslandMap;
        typedef pair<string, vector<TPulseIsland*> > TStringPulseIslandPair;
        typedef map<string, vector<TPulseIsland*> >::iterator map_iterator;

        // Fetch a reference to the gData structure that stores a map
        // of (bank_name, vector<TPulseIsland*>) pairs
        TStringPulseIslandMap& pulse_islands_map =
                gData->fPulseIslandToChannelMap;

        // SiL
        // Just loop over the banks we know are connected to these detectors
        std::vector<std::map<double, double> > thinSiLeft_TimeHeightMaps;
        for (std::vector<std::string>::iterator bankNameIter = left_thin_banknames.begin();
             bankNameIter != left_thin_banknames.end(); bankNameIter++) {

          std::vector<TPulseIsland*> temp_vector = gData->fPulseIslandToChannelMap[*bankNameIter];
          //skip bank if it doesen't exist
          if (temp_vector.size() == 0) continue;
          
          thinSiLeft_TimeHeightMaps.push_back(MakePulseTimeHeightMap(temp_vector));
        }
        std::map<double, double> thinLeftSum;
        if (thinSiLeft_TimeHeightMaps.size() != 0)
          thinLeftSum = MakePulseHeightSums(thinSiLeft_TimeHeightMaps, 1000);

        std::vector<std::map<double, double> > totalLeft_TimeHeightMaps;
        for (std::vector<std::string>::iterator leftThickIter = left_thick_banknames.begin();
             leftThickIter != left_thick_banknames.end(); leftThickIter++) {

          std::vector<TPulseIsland*> temp_vector = gData->fPulseIslandToChannelMap[*leftThickIter];
          //skip bank if it doesen't exist
          if (temp_vector.size() == 0) continue;
          
          totalLeft_TimeHeightMaps.push_back(MakePulseTimeHeightMap(temp_vector));
        }
        if (thinLeftSum.size() != 0)
          totalLeft_TimeHeightMaps.push_back(thinLeftSum);

        std::map<double, double> totalLeftSum;
        if (totalLeft_TimeHeightMaps.size() != 0)
          totalLeftSum = MakePulseHeightSums(totalLeft_TimeHeightMaps, 1000);

        for (std::map<double, double>::iterator thinHit = thinLeftSum.begin(), totalHit = totalLeftSum.begin(); 
             thinHit != thinLeftSum.end() || totalHit != totalLeftSum.end(); thinHit++, totalHit++) {

          hdEdx_left->Fill((*thinHit).second, (*totalHit).second);
        }       
        //


        // SiR
        // Just loop over the banks we know are connected to these detectors
        std::vector<std::map<double, double> > thinSiRight_TimeHeightMaps;
        for (std::vector<std::string>::iterator bankNameIter = right_thin_banknames.begin();
             bankNameIter != right_thin_banknames.end(); bankNameIter++) {

          std::vector<TPulseIsland*> temp_vector = gData->fPulseIslandToChannelMap[*bankNameIter];
          //skip bank if it doesen't exist
          if (temp_vector.size() == 0) continue;
          
          thinSiRight_TimeHeightMaps.push_back(MakePulseTimeHeightMap(temp_vector));
        }
        std::map<double, double> thinRightSum;
        if (thinSiRight_TimeHeightMaps.size() != 0)
          thinRightSum = MakePulseHeightSums(thinSiRight_TimeHeightMaps, 1000);

        std::vector<std::map<double, double> > totalRight_TimeHeightMaps;
        for (std::vector<std::string>::iterator rightThickIter = right_thick_banknames.begin();
             rightThickIter != right_thick_banknames.end(); rightThickIter++) {

          std::vector<TPulseIsland*> temp_vector = gData->fPulseIslandToChannelMap[*rightThickIter];
          //skip bank if it doesen't exist
          if (temp_vector.size() == 0) continue;
          
          totalRight_TimeHeightMaps.push_back(MakePulseTimeHeightMap(temp_vector));
        }
        if (thinRightSum.size() != 0)
          totalRight_TimeHeightMaps.push_back(thinRightSum);

        std::map<double, double> totalRightSum;
        if (totalRight_TimeHeightMaps.size() != 0)
          totalRightSum = MakePulseHeightSums(totalRight_TimeHeightMaps, 1000);

        for (std::map<double, double>::iterator thinHit = thinRightSum.begin(), totalHit = totalRightSum.begin(); 
             thinHit != thinRightSum.end() || totalHit != totalRightSum.end(); thinHit++, totalHit++) {

          hdEdx_right->Fill((*thinHit).second, (*totalHit).second);
        }       

        return SUCCESS;
}

INT MdEdxPlot_init

(

)

This method initializes histograms.

Definition at line 72 of file MdEdxPlot.cpp.

References TSetupData::fBankToDetectorMap, TSetupData::GetDetectorName(), hdEdx_left, hdEdx_right, and SUCCESS.

{
  // The dE/dx histogram is created for the left and right arms of the detector:
  // energy in Si1 (x-axis) vs total energy in Si1 + Si2 (y-axis)

  std::map<std::string, std::string> bank_to_detector_map = gSetup->fBankToDetectorMap;
  for(std::map<std::string, std::string>::iterator mapIter = bank_to_detector_map.begin(); 
      mapIter != bank_to_detector_map.end(); mapIter++) { 

    std::string bankname = mapIter->first;
    std::string detname = gSetup->GetDetectorName(bankname);

    if (detname.substr(0,4) == "SiL1" && detname.substr(7,4) == "slow")
      left_thin_banknames.push_back(bankname);
    else if (detname.substr(0,4) == "SiL2" && detname.substr(5,4) == "slow")
      left_thick_banknames.push_back(bankname);
    else if (detname.substr(0,4) == "SiR1" && detname.substr(7,4) == "slow")
      right_thin_banknames.push_back(bankname);
    else if (detname.substr(0,4) == "SiR2" && detname.substr(5,4) == "slow")
      right_thick_banknames.push_back(bankname);
  }

  int max_adc_value = 4096; // probably want to check which channel we're in before setting this
  max_adc_value = 7500;
  // hdEdx
  std::string histname = "hdEdx_left";
  std::string histtitle = "dE/dx plot for the left arm";
  hdEdx_left = new TH2D(histname.c_str(), histtitle.c_str(), max_adc_value,0,max_adc_value, max_adc_value,0,max_adc_value);
  hdEdx_left->GetXaxis()->SetTitle("Pulse Height in Si1");
  hdEdx_left->GetYaxis()->SetTitle("Total Pulse Height in Si1 and Si2");

  histname = "hdEdx_right";
  histtitle = "dE/dx plot for the right arm";
  hdEdx_right = new TH2D(histname.c_str(), histtitle.c_str(), max_adc_value,0,max_adc_value, max_adc_value,0,max_adc_value);
  hdEdx_right->GetXaxis()->SetTitle("Pulse Height in Si1");
  hdEdx_right->GetYaxis()->SetTitle("Total Pulse Height in Si1 and Si2");

  return SUCCESS;
}

Variable Documentation

TGlobalData* gData

Definition at line 54 of file analyzer.cpp.

TSetupData* gSetup

Definition at line 55 of file analyzer.cpp.

HNDLE hDB

Definition at line 11 of file mucap_compress.cpp.

TH2D* hdEdx_left

static

Definition at line 48 of file MdEdxPlot.cpp.

Referenced by MdEdxPlot(), and MdEdxPlot_init().

TH2D* hdEdx_right

static

Definition at line 49 of file MdEdxPlot.cpp.

Referenced by MdEdxPlot(), and MdEdxPlot_init().

std::vector<std::string> left_thick_banknames

Definition at line 52 of file MdEdxPlot.cpp.

std::vector<std::string> left_thin_banknames

Definition at line 51 of file MdEdxPlot.cpp.

ANA_MODULE MdEdxPlot_module

Initial value:

=
{
        "MdEdxPlot",                    
        "Andrew Edmonds",              
        MdEdxPlot,                      
        NULL,                          
        NULL,                          
        MdEdxPlot_init,                 
        NULL,                          
        NULL,                          
        0,                             
        NULL,                          
}

Definition at line 56 of file MdEdxPlot.cpp.

std::vector<std::string> right_thick_banknames

Definition at line 54 of file MdEdxPlot.cpp.

std::vector<std::string> right_thin_banknames

Definition at line 53 of file MdEdxPlot.cpp.