Classes for Shower were added, reading the database, readening raw hits, some outputs

Makefile

@@ -13,7 +13,8 @@ SRC  =  src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
 	src/THcDCHit.cxx \
 	src/THcHitList.cxx src/THcDetectorMap.cxx src/THcHodoscope.cxx \
 	src/THcHallCSpectrometer.cxx src/THcDriftChamber.cxx \
-	src/THcScintillatorPlane.cxx src/THcSignalHit.cxx
+	src/THcScintillatorPlane.cxx src/THcSignalHit.cxx \
+	src/THcShower.cxx src/THcShowerHit.cxx 
 
 # Name of your package. 
 # The shared library that will be built will get the name lib$(PACKAGE).so

examples/hodtest.C

@@ -17,6 +17,7 @@
 
   // Add hodoscope
   HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
+  HMS->AddDetector( new THcShower("Cal", "Shower" ));
   //  HMS->AddDetector( new THcDriftChamber("dc", "Drift Chambers" ));
 
   // Set up the analyzer - we use the standard one,

src/HallC_LinkDef.h

@@ -20,5 +20,7 @@
 #pragma link C++ class THcHallCSpectrometer+;
 #pragma link C++ class THcScintillatorPlane+;
 #pragma link C++ class THcSignalHit+;
+#pragma link C++ class THcShowerHit+;
+#pragma link C++ class THcShower+;
 
 #endif

src/THcShower.cxx

+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// THcShower                                                                 //
+//                                                                           //
+// Shower counter class, describing a generic segmented shower detector.     //
+//                                                                           //
+//                                                                           //
+//                                                                           //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcShower.h"
+#include "THaEvData.h"
+#include "THaDetMap.h"
+#include "THcDetectorMap.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "THaTrack.h"
+#include "TClonesArray.h"
+#include "TMath.h"
+
+#include "THaTrackProj.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+using namespace std;
+
+//_____________________________________________________________________________
+THcShower::THcShower( const char* name, const char* description,
+				  THaApparatus* apparatus ) :
+  THaNonTrackingDetector(name,description,apparatus)
+{
+  // Constructor
+
+//  fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+}
+
+//_____________________________________________________________________________
+THcShower::THcShower( ) :
+  THaNonTrackingDetector()
+{
+  // Constructor
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
+{
+  static const char* const here = "Init()";
+
+  if( THaNonTrackingDetector::Init( date ) )
+    return fStatus;
+
+
+  // Should probably put this in ReadDatabase as we will know the
+  // maximum number of hits after setting up the detector map
+
+  THcHitList::InitHitList(fDetMap, "THcShowerHit", 100);
+
+  // Will need to determine which apparatus it belongs to and use the
+  // appropriate detector ID in the FillMap call
+  if( gHcDetectorMap->FillMap(fDetMap, "HCAL") < 0 ) {
+    Error( Here(here), "Error filling detectormap for %s.", 
+	     "HCAL");
+      return kInitError;
+  }
+
+  return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcShower::ReadDatabase( const TDatime& date )
+{
+  // Read this detector's parameters from the database file 'fi'.
+  // This function is called by THaDetectorBase::Init() once at the
+  // beginning of the analysis.
+  // 'date' contains the date/time of the run being analyzed.
+
+  //  static const char* const here = "ReadDatabase()";
+
+  // Read data from database 
+  // Pull values from the THcParmList instead of reading a database
+  // file like Hall A does.
+
+  //  DBRequest list[] = {
+  //    { "TDC_offsetsL", fLOff, kDouble, fNelem },
+  //    { "TDC_offsetsR", fROff, kDouble, fNelem },
+  //    { "ADC_pedsL", fLPed, kDouble, fNelem },
+  //    { "ADC_pedsR", fRPed, kDouble, fNelem },
+  //    { "ADC_coefL", fLGain, kDouble, fNelem },
+  //    { "ADC_coefR", fRGain, kDouble, fNelem },
+  //    { "TDC_res",   &fTdc2T },
+  //    { "TransSpd",  &fCn },
+  //    { "AdcMIP",    &fAdcMIP },
+  //    { "NTWalk",    &fNTWalkPar, kInt },
+  //    { "Timewalk",  fTWalkPar, kDouble, 2*fNelem },
+  //    { "ReTimeOff", fTrigOff, kDouble, fNelem },
+  //    { "AvgRes",    &fResolution },
+  //    { "Atten",     &fAttenuation },
+  //    { 0 }
+  //  };
+
+  // We will probably want to add some kind of method to gHcParms to allow
+  // bulk retrieval of parameters of interest.
+
+  // Will need to determine which spectrometer in order to construct
+  // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)
+
+cout << "THcShower::ReadDatabase called " << GetName() << endl;
+
+  NLayers = 4;			// Hardwire for now
+
+  BlockThick = new Double_t [NLayers];
+
+  BlockThick[0] = *(Double_t *)gHcParms->Find("hcal_1pr_thick")->GetValuePointer();
+  BlockThick[1] = *(Double_t *)gHcParms->Find("hcal_2ta_thick")->GetValuePointer();
+  BlockThick[2] = *(Double_t *)gHcParms->Find("hcal_3ta_thick")->GetValuePointer();
+  BlockThick[3] = *(Double_t *)gHcParms->Find("hcal_4ta_thick")->GetValuePointer();
+
+cout << "Block thickness: " << BlockThick[2]  << endl;
+
+fNBlocks = new Int_t [NLayers];
+
+  fNBlocks[0] = *(Int_t *)gHcParms->Find("hcal_1pr_nr")->GetValuePointer();
+  fNBlocks[1] = *(Int_t *)gHcParms->Find("hcal_2ta_nr")->GetValuePointer();
+  fNBlocks[2] = *(Int_t *)gHcParms->Find("hcal_3ta_nr")->GetValuePointer();
+  fNBlocks[3] = *(Int_t *)gHcParms->Find("hcal_4ta_nr")->GetValuePointer();
+
+  cout << "Number of blocks per layer: " << fNBlocks[2]  << endl;
+
+  fNLayerZPos = new Double_t [NLayers];
+
+  fNLayerZPos[0] = *(Double_t *)gHcParms->Find("hcal_1pr_zpos")->GetValuePointer();
+  fNLayerZPos[1] = *(Double_t *)gHcParms->Find("hcal_2ta_zpos")->GetValuePointer();
+  fNLayerZPos[2] = *(Double_t *)gHcParms->Find("hcal_3ta_zpos")->GetValuePointer();
+  fNLayerZPos[3] = *(Double_t *)gHcParms->Find("hcal_4ta_zpos")->GetValuePointer();
+
+cout << "Z Position: " << fNLayerZPos[2]  << endl;
+
+  XPos = new Double_t [2*NLayers];
+
+  XPos[0] = *(Double_t *)gHcParms->Find("hcal_1pr_left")->GetValuePointer();
+  XPos[1] = *(Double_t *)gHcParms->Find("hcal_1pr_right")->GetValuePointer();
+  XPos[2] = *(Double_t *)gHcParms->Find("hcal_2ta_left")->GetValuePointer();
+  XPos[3] = *(Double_t *)gHcParms->Find("hcal_2ta_right")->GetValuePointer();
+  XPos[4] = *(Double_t *)gHcParms->Find("hcal_3ta_left")->GetValuePointer();
+  XPos[5] = *(Double_t *)gHcParms->Find("hcal_3ta_right")->GetValuePointer();
+  XPos[6] = *(Double_t *)gHcParms->Find("hcal_4ta_left")->GetValuePointer();
+  XPos[7] = *(Double_t *)gHcParms->Find("hcal_4ta_right")->GetValuePointer();
+
+cout << "X Positions: " << XPos[0]  << ", " << XPos[1] << endl;
+
+  YPos = new Double_t* [4];
+  cout << "Y Positions:";
+
+  Double_t* p;
+  Int_t ilayer;
+
+  ilayer = 0; 
+  p = (Double_t *)gHcParms->Find("hcal_1pr_top")->GetValuePointer();
+  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
+  // Print out some parameters just to demonstrate that it works
+  
+  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
+    YPos[ilayer][i] = p[i];
+    cout << " " << YPos[ilayer][i];
+  }
+  cout << endl;
+
+  ilayer = 1; 
+  p = (Double_t *)gHcParms->Find("hcal_2ta_top")->GetValuePointer();
+  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
+  
+   for(Int_t i=0;i<fNBlocks[ilayer];i++) {
+    YPos[ilayer][i] = p[i];
+  }
+ 
+  ilayer = 2; 
+  p = (Double_t *)gHcParms->Find("hcal_3ta_top")->GetValuePointer();
+  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
+
+  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
+    YPos[ilayer][i] = p[i];
+  }
+
+  ilayer = 3; 
+  p = (Double_t *)gHcParms->Find("hcal_4ta_top")->GetValuePointer();
+  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
+
+  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
+    YPos[ilayer][i] = p[i];
+  }
+
+  fIsInit = true;
+
+  return kOK;
+}
+
+
+
+//_____________________________________________________________________________
+Int_t THcShower::DefineVariables( EMode mode )
+{
+  // Initialize global variables and lookup table for decoder
+
+  if( mode == kDefine && fIsSetup ) return kOK;
+  fIsSetup = ( mode == kDefine );
+
+  // Register variables in global list
+
+ //  RVarDef vars[] = {
+ //   { "nhit",   "Number of hits",                     "fNhits" },
+ //   { "a",      "Raw ADC amplitudes",                 "fA" },
+ //   { "a_p",    "Ped-subtracted ADC amplitudes",      "fA_p" },
+ //   { "a_c",    "Calibrated ADC amplitudes",          "fA_c" },
+ //   { "asum_p", "Sum of ped-subtracted ADCs",         "fAsum_p" },
+ //   { "asum_c", "Sum of calibrated ADCs",             "fAsum_c" },
+ //   { "nclust", "Number of clusters",                 "fNclust" },
+ //   { "e",      "Energy (MeV) of largest cluster",    "fE" },
+ //   { "x",      "x-position (cm) of largest cluster", "fX" },
+ //   { "y",      "y-position (cm) of largest cluster", "fY" },
+ //   { "mult",   "Multiplicity of largest cluster",    "fMult" },
+ //   { "nblk",   "Numbers of blocks in main cluster",  "fNblk" },
+ //   { "eblk",   "Energies of blocks in main cluster", "fEblk" },
+ //   { "trx",    "track x-position in det plane",      "fTRX" },
+ //   { "try",    "track y-position in det plane",      "fTRY" },
+ //   { 0 }
+ // };
+ //  return DefineVarsFromList( vars, mode );
+  return kOK;
+}
+
+//_____________________________________________________________________________
+THcShower::~THcShower()
+{
+  // Destructor. Remove variables from global list.
+
+  if( fIsSetup )
+    RemoveVariables();
+  if( fIsInit )
+    DeleteArrays();
+  if (fTrackProj) {
+    fTrackProj->Clear();
+    delete fTrackProj; fTrackProj = 0;
+  }
+}
+
+//_____________________________________________________________________________
+void THcShower::DeleteArrays()
+{
+  // Delete member arrays. Used by destructor.
+
+  delete [] BlockThick;  BlockThick = NULL;
+  delete [] fNBlocks;  fNBlocks = NULL;
+  delete [] fNLayerZPos;  fNLayerZPos = NULL;
+  delete [] XPos;  XPos = NULL;
+  delete [] YPos;  YPos = NULL;
+  //delete [] fSpacing;  fSpacing = NULL;
+  //delete [] fCenter;   fCenter = NULL; // This 2D. What is correct way to delete?
+}
+
+//_____________________________________________________________________________
+inline 
+void THcShower::ClearEvent()
+{
+  // Reset per-event data.
+
+  fTrackProj->Clear();
+}
+
+//_____________________________________________________________________________
+Int_t THcShower::Decode( const THaEvData& evdata )
+{
+
+  // Get the Hall C style hitlist (fRawHitList) for this event
+  Int_t nhits = THcHitList::DecodeToHitList(evdata);
+
+//   fRawHitList is TClones array of THcShowerHit objects
+  for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+    THcShowerHit* hit = (THcShowerHit *) fRawHitList->At(ihit);
+    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+	 << hit->fADC_pos << " " << hit->fADC_neg << " "  << endl;
+  }
+  cout << endl;
+
+  return nhits;
+}
+
+//_____________________________________________________________________________
+Int_t THcShower::ApplyCorrections( void )
+{
+  return(0);
+}
+
+//_____________________________________________________________________________
+Double_t THcShower::TimeWalkCorrection(const Int_t& paddle,
+					     const ESide side)
+{
+  return(0.0);
+}
+
+//_____________________________________________________________________________
+Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
+{
+  // Calculation of coordinates of particle track cross point with shower
+  // plane in the detector coordinate system. For this, parameters of track 
+  // reconstructed in THaVDC::CoarseTrack() are used.
+  //
+  // Apply corrections and reconstruct the complete hits.
+  //
+  //  static const Double_t sqrt2 = TMath::Sqrt(2.);
+  
+  ApplyCorrections();
+
+  return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcShower::FineProcess( TClonesArray& tracks )
+{
+  // Reconstruct coordinates of particle track cross point with shower
+  // plane, and copy the data into the following local data structure:
+  //
+  // Units of measurements are meters.
+
+  // Calculation of coordinates of particle track cross point with shower
+  // plane in the detector coordinate system. For this, parameters of track 
+  // reconstructed in THaVDC::FineTrack() are used.
+
+  return 0;
+}
+
+ClassImp(THcShower)
+////////////////////////////////////////////////////////////////////////////////
+ 

src/THcShower.h

+#ifndef ROOT_THcShower
+#define ROOT_THcShower
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// THcHodoscope                                                              //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "TClonesArray.h"
+#include "THaNonTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcShowerHit.h"
+#include "THcScintillatorPlane.h"
+
+class THaScCalib;
+
+class THcShower : public THaNonTrackingDetector, public THcHitList {
+
+public:
+  THcShower( const char* name, const char* description = "",
+		   THaApparatus* a = NULL );
+  virtual ~THcShower();
+
+  virtual Int_t      Decode( const THaEvData& );
+  virtual EStatus    Init( const TDatime& run_time );
+  virtual Int_t      CoarseProcess( TClonesArray& tracks );
+  virtual Int_t      FineProcess( TClonesArray& tracks );
+  
+  virtual Int_t      ApplyCorrections( void );
+
+  //  Int_t GetNHits() const { return fNhit; }
+  
+  Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
+  const TClonesArray* GetTrackHits() const { return fTrackProj; }
+  
+  friend class THaScCalib;
+
+  THcShower();  // for ROOT I/O
+protected:
+
+  // Calibration
+
+  // Per-event data
+
+
+  // Potential Hall C parameters.  Mostly here for demonstration
+  char** LayerNames;
+  Int_t NLayers;
+  Double_t* fNLayerZPos;		// Z position of front of shower counter layers
+  Double_t* BlockThick;		// Thickness of shower counter blocks, blocks
+  Int_t* fNBlocks;           // Number of shower counter blocks per layer
+  Double_t** YPos;		//X,Y positions of shower counter blocks
+  Double_t* XPos;
+
+  THcScintillatorPlane** fPlane; // List of plane objects
+
+  TClonesArray*  fTrackProj;  // projection of track onto scintillator plane
+                              // and estimated match to TOF paddle
+  // Useful derived quantities
+  // double tan_angle, sin_angle, cos_angle;
+  
+  //  static const char NDEST = 2;
+  //  struct DataDest {
+  //    Int_t*    nthit;
+  //    Int_t*    nahit;
+  //    Double_t*  tdc;
+  //    Double_t*  tdc_c;
+  //    Double_t*  adc;
+  //    Double_t*  adc_p;
+  //    Double_t*  adc_c;
+  //    Double_t*  offset;
+  //    Double_t*  ped;
+  //    Double_t*  gain;
+  //  } fDataDest[NDEST];     // Lookup table for decoder
+
+  void           ClearEvent();
+  void           DeleteArrays();
+  virtual Int_t  ReadDatabase( const TDatime& date );
+  virtual Int_t  DefineVariables( EMode mode = kDefine );
+
+  enum ESide { kLeft = 0, kRight = 1 };
+  
+  virtual  Double_t TimeWalkCorrection(const Int_t& paddle,
+					   const ESide side);
+
+  ClassDef(THcShower,0)   // Generic hodoscope class
+};
+
+////////////////////////////////////////////////////////////////////////////////
+
+#endif
+ 

src/THcShowerHit.cxx

+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// THcShowerHit                                                           //
+//                                                                           //
+// Class representing a single raw hit for a hodoscope paddle                //
+//                                                                           //
+// Contains plane, counter and pos/neg adc and tdc values                    //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcShowerHit.h"
+
+using namespace std;
+
+
+void THcShowerHit::SetData(Int_t signal, Int_t data) {
+  if(signal==0) {
+    fADC_pos = data;
+  } else if (signal==1) {
+    fADC_neg = data;
+  } 
+}
+
+Int_t THcShowerHit::GetData(Int_t signal) {
+  if(signal==0) {
+    return(fADC_pos);
+  } else if (signal==1) {
+    return(fADC_neg);
+  } 
+
+  return(-1); // Actually should throw exception
+}
+
+#if 0
+Int_t THcShowerHit::Compare(const TObject* obj) const
+{
+  // Compare to sort by plane and counter
+
+  const THcShowerHit* hit = dynamic_cast<const THcShowerHit*>(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;
+  else {
+    Int_t c1 = fCounter;
+    Int_t c2 = hit->fCounter;
+    if(c1 < c2) return -1;
+    else if (c1 == c2) return 0;
+    else return 1;
+  }
+}
+#endif
+//_____________________________________________________________________________
+THcShowerHit& THcShowerHit::operator=( const THcShowerHit& rhs )
+{
+  // Assignment operator.
+
+  THcRawHit::operator=(rhs);
+  if ( this != &rhs ) {
+    fPlane = rhs.fPlane;
+    fCounter = rhs.fCounter;
+    fADC_pos = rhs.fADC_pos;
+    fADC_neg = rhs.fADC_neg;
+
+  }
+  return *this;
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+ClassImp(THcShowerHit)
+
+ 

src/THcShowerHit.h

+#ifndef ROOT_THcShowerHit
+#define ROOT_THcShowerHit
+
+#include "THcRawHit.h"
+
+class THcShowerHit : public THcRawHit {
+
+ public:
+
+ THcShowerHit(Int_t plane=0, Int_t counter=0) : THcRawHit(plane, counter), 
+    fADC_pos(-1), fADC_neg(-1){
+  }
+  THcShowerHit& operator=( const THcShowerHit& );
+  virtual ~THcShowerHit() {}
+
+  virtual void Clear( Option_t* opt="" )
+    { fADC_pos = -1; fADC_neg = -1; }
+
+  void SetData(Int_t signal, Int_t data);
+  Int_t GetData(Int_t signal);
+
+  //  virtual Bool_t  IsSortable () const {return kTRUE; }
+  //  virtual Int_t   Compare(const TObject* obj) const;
+
+  Int_t fADC_pos;
+  Int_t fADC_neg;
+
+ protected:
+
+ private:
+
+  ClassDef(THcShowerHit, 0);	// Shower hit class
+};  
+
+#endif
+