Skip to content
Snippets Groups Projects
Commit bd295b0b authored by Buddhini Waidyawansa's avatar Buddhini Waidyawansa Committed by Stephen A. Wood
Browse files

Add THcRasterRawHit, THcRasterBeam classes. Add raster_test.C script.

  THcRasterRawHit class to helps with the decoding of the raster
     signals via THcHitList class.
  THcRasteredBeam class derives from THaApparatus to handle rastered beam.
  Modify Makefile and HallC_LinkDef.h to compile new classes
parent 642bccd8
No related branches found
No related tags found
No related merge requests found
......@@ -24,7 +24,10 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
src/THcRawShowerHit.cxx \
src/THcAerogel.cxx src/THcAerogelHit.cxx \
src/THcCherenkov.cxx src/THcCherenkovHit.cxx \
src/THcFormula.cxx
src/THcFormula.cxx\
src/THcRaster.cxx\
src/THcRasteredBeam.cxx\
src/THcRasterRawHit.cxx
# Name of your package.
# The shared library that will be built will get the name lib$(PACKAGE).so
......
......@@ -8,6 +8,8 @@ hhodo_plane_names = "1x 1y 2x 2y"
hcal_num_layers = 4
rraster_num_signals = 4
# Exclusion band width for the calorimeter's fiducial volume.
hcal_fv_delta = 5.
......@@ -83,4 +85,5 @@ hcer_tot_pmts = 2
# parameter list for the raster
# Maybe we can add this to geabm.param later ?
graster_channels = "xsync xsig ysync ysig"
rraster_signal_names = "xsync xsig ysync ysig"
......@@ -2,82 +2,103 @@
{
//
// Steering script to test hodoscope decoding
// Steering script to test raster signal decoding
//
Int_t RunNumber=50017;
char* RunFileNamePattern="daq04_%d.log.0";
// Open the database
//
gHcParms->Define("gen_run_number", "Run Number", RunNumber);
gHcParms->AddString("g_ctp_database_filename", "DBASE/test.database");
gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
// g_ctp_parm_filename and g_decode_map_filename should now be defined
// Open and load parameter files
//
gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
// Constants not in ENGINE PARAM files that we want to be
// configurable
// parameters not found in usual engine parameter files
gHcParms->Load("PARAM/hcana.param");
// Generate db_cratemap to correspond to map file contents
// Generate db_cratemap to correspond to map file contents
// make_cratemap.pl scripts reads a Hall C style MAP file and output a
// Hall A style crate map DB file
//
char command[100];
sprintf(command,"./make_cratemap.pl < %s > db_cratemap.dat",gHcParms->GetString("g_decode_map_filename"));
system(command);
// Load the Hall C style detector map
//
gHcDetectorMap=new THcDetectorMap();
gHcDetectorMap->Load(gHcParms->GetString("g_decode_map_filename"));
// Set up the equipment to be analyzed.
// Set up the equipment to be analyzed.
//
// HMS and its detectors
THaApparatus* HMS = new THcHallCSpectrometer("H","HMS");
gHaApps->Add( HMS );
// Add hodoscope
HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
HMS->AddDetector( new THcShower("cal", "Shower" ));
HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
THcAerogel* aerogel = new THcAerogel("aero", "Aerogel Cerenkov" );
HMS->AddDetector( aerogel );
THcCherenkov* cherenkov = new THcCherenkov("cher", "Gas Cerenkov" );
HMS->AddDetector( cherenkov );
HMS->AddDetector( new THcAerogel("aero", "Aerogel Cerenkov" ));
HMS->AddDetector( new THcCherenkov("cher", "Gas Cerenkov" ));
// Beamline and its detectors
THaApparatus * BEAM = new THcRasteredBeam("rb","Beamline");
gHaApps->Add( BEAM );
// Set up the analyzer - we use the standard one,
// but this could be an experiment-specific one as well.
// The Analyzer controls the reading of the data, executes
// tests/cuts, loops over Acpparatus's and PhysicsModules,
// and executes the output routines.
//
THcAnalyzer* analyzer = new THcAnalyzer;
// A simple event class to be output to the resulting tree.
// Creating your own descendant of THaEvent is one way of
// defining and controlling the output.
//
THaEvent* event = new THaEvent;
// Define the run(s) that we want to analyze.
// We just set up one, but this could be many.
//
char RunFileName[100];
sprintf(RunFileName,RunFileNamePattern,RunNumber);
THaRun* run = new THaRun(RunFileName);
// Eventually need to learn to skip over, or properly analyze
// the pedestal events
run->SetEventRange(1,5000);// Physics Event number, does not
// include scaler or control events
//
run->SetEventRange(1,2);// Physics Event number, does not
// include scaler or control events
// Define the analysis parameters
analyzer->SetEvent( event );
analyzer->SetOutFile( "hodtest.root" );
//
analyzer->SetEvent(event);
analyzer->SetOutFile("raster_test.root");
analyzer->SetOdefFile("output.def");
analyzer->SetCutFile("hodtest_cuts.def"); // optional
// analyzer->SetCutFile("rastertest_cuts.def"); // optional
analyzer->SetCountMode(2);// Counter event number same as gen_event_ID_number
// File to record cuts accounting information
// analyzer->SetSummaryFile("summary_example.log"); // optional
analyzer->Process(run); // start the actual analysis
// start the actual analysis
//
analyzer->Process(run);
analyzer->PrintReport("report.template","report.out");
}
......@@ -48,5 +48,8 @@
#pragma link C++ class THcCherenkov+;
#pragma link C++ class THcCherenkovHit+;
#pragma link C++ class THcFormula+;
#pragma link C++ class THcRaster+;
#pragma link C++ class THcRasteredBeam+;
#pragma link C++ class THcRasterRawHit+;
#endif
......@@ -12,12 +12,22 @@
///////////////////////////////////////////////////////////////////////////////
#include "THcRaster.h"
#include "THcHitList.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcParmList.h"
#include "THcGlobals.h"
#include "THaGlobals.h"
//#include "THcHitList.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
using namespace std;
......@@ -27,83 +37,239 @@ THcRaster::THcRaster( const char* name, const char* description,
THaApparatus* apparatus ) :
THaBeamDet(name,description,apparatus)
{
// // Initializing channels
fAnalyzePedestals = 0;
fRawPos[0]=0;
fRawPos[1]=0;
// fRasterXADC = new TClonesArray("THcSignalHit",16);
// fRasterYADC = new TClonesArray("THcSignalHit",16);
}
//_____________________________________________________________________________
THcRaster::THcRaster( ) :
THaBeamDet()
{
// Default constructor
}
//_____________________________________________________________________________
THcRaster::~THRaster()
THcRaster::~THcRaster()
{
// Distructor
}
//_____________________________________________________________________________
void THcRaster::ClearEvent(){
// nothing for now
}
//_____________________________________________________________________________
void THcRaster::InitializeReconstruction()
{
// do nothing for now
}
//_____________________________________________________________________________
// ReadDatabase: if detectors cant be added to detmap
Int_t THcRaster::ReadDatabase( const TDatime& date )
{
static const char* const here = "THcRaster::ReadDatabase";
InitializeReconstruction();
// Read parameters such as calibration factor, of this detector from the database.
// static const char* const here = "THcRaster::ReadDatabase";
// InitializeReconstruction();
// char prefix[2];
// cout << " THcRaster::ReadDatabase GetName() called " << GetName() << endl;
// prefix[0]=tolower(GetName()[0]);
// prefix[1]='\0';
char prefix[2];
// string names;
// DBRequest list[]={
// {"raster_channel_name",&names,kString},
// {"raster_channel_number",&names,kString},
// {0}
// };
cout << " GetName() " << GetName() << endl;
// // get the channel assignments from the parameter file
// gHcParms->LoadParmValues((DBRequest*)&list,prefix);
prefix[0]=tolower(GetName()[0]);
prefix[1]='\0';
// std::cout <<"Raster channel no : " << std::endl;
// std::cout <<" X signal channel = "<<fChNO[0]<<std::endl;
// std::cout <<" Y signal channel = "<<fChNO[1]<<std::endl;
string channelNames;
DBRequest channellist[]={
{"graster_channels",&channelNames,kString},
{0}
};
// get the channel list from the parameter file
gHcParms->LoadParmValues((DBRequest*)&channellist,prefix);
std::cout << "Raster channel list : " << channellist << std::endl;
// std::vector<string> channel_names = vsplit(names);
// Int_t NChannels = channel_names.size();
// for(Int_t i=0;i<NChannels;i++) {
// // fPlaneNames[i] = new char[plane_names[i].length()];
// // strcpy(fPlaneNames[i], plane_names[i].c_str());
// //std::cout<<" channel = "<<channel_names.at(i)<<std::endl;
// std::cout<<" Channel numbers = "<<
// }
vector<string> channel_names = vsplit(channellist);
Int_t NChannels = channel_names.size_of();
for(Int_t i=0;i<NChannels;i++) {
// fPlaneNames[i] = new char[plane_names[i].length()];
// strcpy(fPlaneNames[i], plane_names[i].c_str());
std::cout<<" channel = "<<channel_names.at(i)<<std::endl;
// Read the detector map file
// Based on the ungly method used in THaRaster class. I have no idea
// how to do this more elegantly.
// Buddhini
// const int LEN=100;
// char buf[LEN];
// FILE* fi = OpenFile( date );
// if( !fi ) return kFileError;
// fDetMap->Clear();
// int first_chan, crate, dummy, slot, first, last, modulid;
// std::cout<<"$$$$$$$$$$$$ THcRaster detmap"<<std::endl;
// do {
// fgets( buf, LEN, fi);
// sscanf(buf,"%d %d %d %d %d %d %d",&first_chan, &crate, &dummy, &slot, &first, &last, &modulid);
// std::cout<<"$$$$$$$$$$$$ "<<crate<<std::endl;
// // if (first_chan>=0) {
// // if ( fDetMap->AddModule (crate, slot, first, last, first_chan )<0) {
// // Error( Here(here), "Couldnt add Raster to DetMap!");
// // fclose(fi);
// // return kInitError;
// // }
// // }
// } while (first_chan>=0);
return kOK;
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcRaster::Init( const TDatime& date )
{
cout << "THcRaster::Init()" << endl;
THcHitList::InitHitList(fDetMap,"THcRasterRawHit",4);
EStatus status;
if( (status = THaBeamDet::Init( date )) )
return fStatus=status;
// Fill detector map with RASTER type channels
if( gHcDetectorMap->FillMap(fDetMap, "RASTER") < 0 ) {
static const char* const here = "Init()";
Error( Here(here), "Error filling detectormap for %s.",
"RASTER");
return kInitError;
}
return kOK
return fStatus = kOK;
}
//_____________________________________________________________________________
void THcRaster::AccumulatePedestals(TClonesArray* rawhits)
{
// Extract data from the hit list, accumulating into arrays for
// calculating pedestals
}
//_____________________________________________________________________________
void THcRaster::CalculatePedestals( )
{
// Use the accumulated pedestal data to calculate pedestals
}
//_____________________________________________________________________________
Int_t THcRaster::Decode( const THaEvData& evdata )
{
cout << "THcRaster::Decode()" << endl;
// clears the event structure
// loops over all channels defined in the detector map
// copies raw data into local variables
// performs pedestal subtraction
ClearEvent();
// Get the Hall C style hitlist (fRawHitList) for this event
// empty for now while I work on reading the channel map.
Int_t fNhits = THcHitList::DecodeToHitList(evdata);
std::cout<<"$$$$ Hits.."<<fNhits<<std::endl;
// // Get the pedestals
// if(gHaCuts->Result("Pedestal_event")) {
// Int_t nexthit = 0;
// AccumulatePedestals(fRawHitList);
// fAnalyzePedestals = 1; // Analyze pedestals first normal events
// return(0);
// }
// if(fAnalyzePedestals) {
// CalculatePedestals();
// fAnalyzePedestals = 0; // Don't analyze pedestals next event
// }
Int_t ihit = 0;
std::cout<<"$$$$ Hits.."<<fNhits<<std::endl;
while(ihit < fNhits) {
//std::cout<<"Get raw hit.."<<std::endl;
THcRasterRawHit* hit = (THcRasterRawHit *) fRawHitList->At(ihit);
fRawADC[0] = hit->fADC_xsig;
fRawADC[1] = hit->fADC_ysig;
for(Int_t i=0;i<2;i++){
fRawPos[i] = fRawADC[i]-fPedADC[i];
// std::cout<<" Raw ADC = "<<fRawADC[i]
// <<" ped = "<<fPedADC[i]
// <<" corrected ="<<fRawPos[i]<<std::endl;
}
}
for(Int_t ihit = 0; ihit < fNhits ; ihit++) {
THcRasterRawHit* hit = (THcRasterRawHit *) fRawHitList->At(ihit);
cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
<< hit->fADC_xsig << " " << hit->fADC_ysig << " " << hit->fADC_xsync
<< " " << hit->fADC_ysync << endl;
}
cout << endl;
//Get X and Y ADC values
return 0;
}
//_____________________________________________________________________________
Int_t THaRaster::Process( ){
Int_t THcRaster::Process( ){
/* From ENGINE/g_analyze_misc.f
Fast Raster Signals:
===================
gfrx_raw_adc = gmisc_dec_data(14,2)
gfry_raw_adc = gmisc_dec_data(16,2)
calculate raster position from ADC value.
gfrx_adc = gfrx_raw_adc - gfrx_adc_ped
gfry_adc = gfry_raw_adc - gfry_adc_ped
*/
// empty for now while I work on reading the channel map.
return o;
return 0;
}
......
#ifndef ROOT_THcRaster
#define ROOT_THcRaster
///////////////////////////////////////////////////////////////////////////////
// //
// THcRaster //
......@@ -9,36 +8,53 @@
///////////////////////////////////////////////////////////////////////////////
#include "THaBeamDet.h"
#include "TVector.h"
#include "TClonesArray.h"
#include "THcHitList.h"
#include "THcDetectorMap.h"
#include "THcRasterRawHit.h"
#include "THaCutList.h"
class THcRaster : public THaBeamDet {
class THcRaster : public THaBeamDet, public THcHitList {
public:
THcRaster( const char* name, const char* description ="" ,
THaApparatus* a = NULL );
THcRaster(const char* name, const char* description ="",THaApparatus* a = NULL );
virtual ~THaRaster();
~THcRaster();
EStatus Init( const TDatime& run_time );
virtual Int_t Decode( const THaEvData& );
virtual Int_t Process();
Int_t Decode( const THaEvData& );
Int_t Process();
virtual TVector3 GetPosition() const { return fPosition[2]; }
virtual TVector3 GetDirection() const { return NULL; } // Hall C we don't use raster direction yet.
TVector3 GetPosition() const { return fPosition[2]; }
TVector3 GetDirection() const { return fDirection; } // Hall C we don't use raster direction yet.
Double_t GetCurrentX() { return fRawPos(0); }
Double_t GetCurrentY() { return fRawPos(1); }
Double_t GetCurrentX() { return fRawPos[0]; }
Double_t GetCurrentY() { return fRawPos[1]; }
protected:
void InitializeReconstruction();
void ClearEvent();
virtual Int_t ReadDatabase( const TDatime& date );
Int_t ReadDatabase( const TDatime& date );
Double_t fRawADC[2]; // raw ADC values
Double_t fPedADC[2]; // ADC poedestals
TVector fRawPos[2]; // current in Raster ADCs for position
Double_t fRawPos[2]; // current in Raster ADCs for position
TVector3 fPosition[2]; // Beam position at 1st, 2nd BPM or at the target (meters)
TVector3 fPosOff[2]; // pedestals
TVector3 fPosPed[2]; // pedestals
TVector3 fDirection; // Beam angle at the target (meters)
// Array to store channel numbers corresponding to X and Ysignals. 0 - X, 1 - Y.
Double_t fChN0[2];
private:
Bool_t fAnalyzePedestals;
void CalculatePedestals();
void AccumulatePedestals(TClonesArray* rawhits);
ClassDef(THcRaster, 0); // add THcRaster to ROOT library
};
......
// Author : Buddhini Waidyawansa
// Date : 23-01-2014
///////////////////////////////////////////////////////////////////////////////
// //
// THcRasterRawHit //
// //
// Class representing a single raw hit for the raster //
// //
// Contains the X, Y raster voltage signals and sync signals //
// //
///////////////////////////////////////////////////////////////////////////////
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include "THcRasterRawHit.h"
using namespace std;
void THcRasterRawHit::SetData(Int_t signal, Int_t data) {
if(signal==0) {
fADC_xsync = data;
} else if (signal==1) {
fADC_xsig = data;
} else if(signal==2) {
fADC_ysync = data;
} else if (signal==3) {
fADC_ysig = data;
}
}
Int_t THcRasterRawHit::GetData(Int_t signal) {
if(signal==1) {
return(fADC_xsync);
} else if (signal==2) {
return(fADC_xsig);
} else if(signal==3) {
return(fADC_ysync);
} else if (signal==4) {
return(fADC_ysig);
}
return(-1); // Actually should throw exception
}
// Int_t THcRasterRawHit::Compare(const TObject* obj) const
// {
// // Compare to sort by the plane
// // There is only one raster so no need for an additional check on the counter
// const THcRasterRawHit* hit = dynamic_cast<const THcRasterRawHit*>(obj);
// if(!hit) return -1;
// Int_t p1 = fPlane;
// Int_t p2 = hit->fPlane;
// if(p1 < p2) return -1;
// else if(p1 > p2) return 1;
// }
//_____________________________________________________________________________
THcRasterRawHit& THcRasterRawHit::operator=( const THcRasterRawHit& rhs )
{
// Assignment operator.
THcRawHit::operator=(rhs);
if ( this != &rhs ) {
fPlane = rhs.fPlane;
fCounter = rhs.fCounter;
fADC_xsync = rhs.fADC_xsync;
fADC_xsig = rhs.fADC_xsig;
fADC_ysync = rhs.fADC_ysync;
fADC_ysig = rhs.fADC_ysig;
}
return *this;
}
//////////////////////////////////////////////////////////////////////////
ClassImp(THcRasterRawHit)
#ifndef ROOT_THcRasterRawHit
#define ROOT_THcRasterRawHit
///////////////////////////////////////////////////////////////////////////////
// //
// THcRasterRawHit //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THcRawHit.h"
class THcRasterRawHit : public THcRawHit {
public:
THcRasterRawHit(Int_t plane=0, Int_t counter=0) : THcRawHit(plane, counter),
fADC_xsig(-1), fADC_ysig(-1),
fADC_xsync(-1), fADC_ysync(-1) {
}
THcRasterRawHit& operator=( const THcRasterRawHit& );
~THcRasterRawHit() {}
void Clear( Option_t* opt="" )
{ fADC_xsig = -1; fADC_ysig = -1; fADC_xsync = -1; fADC_ysync = -1; }
void SetData(Int_t signal, Int_t data);
Int_t GetData(Int_t signal);
//Int_t Compare(TObject* obj);
// signals
Int_t fADC_xsig;
Int_t fADC_ysig;
Int_t fADC_xsync;
Int_t fADC_ysync;
protected:
private:
ClassDef(THcRasterRawHit, 0); // Raw hit class for raster data
};
#endif
// Author : Buddhini Waidyawansa
// Date : 01-08-2014
///////////////////////////////////////////////////////////////////////////////
// //
// THcRasteredBeam //
// //
// A class to handle the raster processing tasks //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THcRasteredBeam.h"
#include "THcRaster.h"
#include "TMath.h"
#include "TDatime.h"
#include "TList.h"
#include "VarDef.h"
//_____________________________________________________________________________
THcRasteredBeam::THcRasteredBeam( const char* name, const char* description ) :
THaBeam( name, description )
{
AddDetector( new THcRaster("Raster","raster",this) );
}
//_____________________________________________________________________________
Int_t THcRasteredBeam::Reconstruct()
{
TIter nextDet( fDetectors );
nextDet.Reset();
// This apparatus assumes that there is only one detector
// in the list. If someone adds detectors by hand, the first
// detector in the list will be used to get the beam position
// the others will be processed
// -- Iam not sure why the code is written like this. But for now, I am
// going to work with this code as it is since all I need for is to
// decode the raster - Buddhini
if (THaBeamDet* theBeamDet=
static_cast<THaBeamDet*>( nextDet() )) {
theBeamDet->Process();
fPosition = theBeamDet->GetPosition();
fDirection = theBeamDet->GetDirection();
}
else {
Error( Here("THcRasteredBeam::Reconstruct"),
"Beamline Detectors Missing in Detector List" );
}
// Process any other detectors that may have been added (by default none)
while (THaBeamDet * theBeamDet=
static_cast<THaBeamDet*>( nextDet() )) {
theBeamDet->Process();
}
Update();
return 0;
}
//_____________________________________________________________________________
ClassImp(THcRasteredBeam)
////////////////////////////////////////////////////////////////////////////////
#ifndef ROOT_THcRasteredBeam
#define ROOT_THcRasteredBeam
///////////////////////////////////////////////////////////////////////////////
// //
// THcRasteredBeam //
// A beam with rastered beam, analyzed event by event using raster currents //
// This is identical to THaRasteredBeam except that we need to use THcRaster //
// For raster signal processing. //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THaBeam.h"
class THcRasteredBeam : public THaBeam {
public:
THcRasteredBeam( const char* name, const char* description) ;
virtual ~THcRasteredBeam(){};
virtual Int_t Reconstruct() ;
ClassDef(THcRasteredBeam, 0); // add THcRasteredBeam to ROOT library
};
////////////////////////////////////////////////////////////////////////////////
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment