diff --git a/engine b/engine
new file mode 120000
index 0000000000000000000000000000000000000000..d2882316dab6a6403777d8369a1239d7f2944bdb
--- /dev/null
+++ b/engine
@@ -0,0 +1 @@
+/u/group/shms/zahmed/analysis/hcana_replay
\ No newline at end of file
diff --git a/examples/THcCherenkov.cxx b/examples/THcCherenkov.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..da107c0fc86b2f37cc47840f52aa29fd8db8f1a5
--- /dev/null
+++ b/examples/THcCherenkov.cxx
@@ -0,0 +1,611 @@
+///////////////////////////////////////////////////////////////////////////////////////
+// //
+// THcCherenkov //
+// //
+// Class for an Cherenkov detector consisting of onw pair of PMT's //
+// //
+// Zafar Ahmed. First attempt. September 27 2013. //
+// //
+///////////////////////////////////////////////////////////////////////////////////////
+
+#include "THcCherenkov.h"
+#include "TClonesArray.h"
+#include "THcSignalHit.h"
+#include "THaEvData.h"
+#include "THaDetMap.h"
+#include "THcDetectorMap.h"
+#include "THcGlobals.h"
+#include "THaCutList.h"
+#include "THcParmList.h"
+#include "THcHitList.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;
+
+using std::cout;
+using std::cin;
+using std::endl;
+
+#include <iomanip>
+using std::setw;
+using std::setprecision;
+
+//_____________________________________________________________________________
+THcCherenkov::THcCherenkov( const char* name, const char* description,
+ THaApparatus* apparatus ) :
+ THaNonTrackingDetector(name,description,apparatus)
+{
+ // Normal constructor with name and description
+ fPosTDCHits = new TClonesArray("THcSignalHit",16);
+ // fNegTDCHits = new TClonesArray("THcSignalHit",16); // check
+ fPosADCHits = new TClonesArray("THcSignalHit",16);
+ // fNegADCHits = new TClonesArray("THcSignalHit",16); // check
+
+// fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+}
+
+//_____________________________________________________________________________
+THcCherenkov::THcCherenkov( ) :
+ THaNonTrackingDetector()
+{
+ // Constructor
+}
+
+//_____________________________________________________________________________
+THcCherenkov::~THcCherenkov()
+{
+ // Destructor
+ delete [] fA_Pos;
+ delete [] fA_Pos_1;
+ // delete [] fA_Neg; // check
+ delete [] fA_Pos_p;
+ // delete [] fA_Neg_p; // check
+ delete [] fT_Pos;
+ // delete [] fT_Neg; // check
+ delete [] fPosGain;
+ // delete [] fNegGain; // check
+ delete [] fPosPedLimit;
+ // delete [] fNegPedLimit; // check
+ delete [] fPosPedMean;
+ // delete [] fNegPedMean; // check
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
+{
+ static const char* const here = "Init()";
+
+ cout << "THcCherenkov::Init " << GetName() << endl;
+
+ // Should probably put this in ReadDatabase as we will know the
+ // maximum number of hits after setting up the detector map
+ THcHitList::InitHitList(fDetMap, "THcCherenkovHit", 100); // 100 is max hits
+
+ EStatus status;
+ if( (status = THaNonTrackingDetector::Init( date )) )
+ return fStatus=status;
+
+ // Will need to determine which apparatus it belongs to and use the
+ // appropriate detector ID in the FillMap call
+ if( gHcDetectorMap->FillMap(fDetMap, "HCER") < 0 ) {
+ Error( Here(here), "Error filling detectormap for %s.",
+ "HCER");
+ return kInitError;
+ }
+
+ return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::ReadDatabase( const TDatime& date )
+{
+ // This function is called by THaDetectorBase::Init() once at the beginning
+ // of the analysis.
+
+ cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
+
+ char prefix[2];
+
+ prefix[0]=tolower(GetApparatus()->GetName()[0]);
+ prefix[1]='\0';
+
+ // DBRequest listextra[]={
+ // {"cher_num_pairs", &fNelem, kInt},
+ // {0}
+ // };
+ fNelem = 2; // Ahmed. One pair of PMTs for Cherenkov detctor
+ // gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);
+
+ fA_Pos = new Float_t[fNelem];
+ fA_Pos_1 = new Float_t[fNelem];
+ // fA_Neg = new Float_t[fNelem]; // check
+ fA_Pos_p = new Float_t[fNelem];
+ // fA_Neg_p = new Float_t[fNelem]; // check
+ fT_Pos = new Float_t[fNelem];
+ // fT_Neg = new Float_t[fNelem]; // check
+
+ fPosGain = new Double_t[fNelem];
+ // fNegGain = new Double_t[fNelem]; // check
+ fPosPedLimit = new Int_t[fNelem];
+ // fNegPedLimit = new Int_t[fNelem]; // check
+ fPosPedMean = new Double_t[fNelem];
+ // fNegPedMean = new Double_t[fNelem]; // check
+
+ DBRequest list[]={
+ {"cer_adc_to_npe", fPosGain, kDouble}, // Ahmed
+ {"cer_ped_limit", fPosPedLimit, kInt}, // Ahmed
+ // {"cher_pos_gain", fPosGain, kDouble}, // Ahmed
+ // {"cher_neg_gain", fNegGain, kDouble}, // Ahmed check
+ // {"cher_pos_ped_limit", fPosPedLimit, kInt}, // Ahmed
+ // {"cher_neg_ped_limit", fNegPedLimit, kInt}, // Ahmed check
+ // {"aero_pos_ped_mean", fPosPedMean, kDouble},
+ // {"aero_neg_ped_mean", fNegPedMean, kDouble},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ fIsInit = true;
+
+ // Create arrays to hold pedestal results
+ InitializePedestals();
+
+ return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::DefineVariables( EMode mode )
+{
+ // Initialize global variables for histogramming and tree
+
+ cout << "THcCherenkov::DefineVariables called " << GetName() << endl;
+
+ if( mode == kDefine && fIsSetup ) return kOK;
+ fIsSetup = ( mode == kDefine );
+
+ // Register variables in global list
+
+ // Do we need to put the number of pos/neg TDC/ADC hits into the variables?
+ // No. They show up in tree as Ndata.H.aero.postdchits for example
+
+ RVarDef vars[] = {
+ {"postdchits", "List of Positive TDC hits","fPosTDCHits.THcSignalHit.GetPaddleNumber()"},
+ // {"negtdchits", "List of Negative TDC hits","fNegTDCHits.THcSignalHit.GetPaddleNumber()"}, // check
+ {"posadchits", "List of Positive ADC hits","fPosADCHits.THcSignalHit.GetPaddleNumber()"},
+ // {"negadchits", "List of Negative ADC hits","fNegADCHits.THcSignalHit.GetPaddleNumber()"}, // check
+ {"apos", "Raw Positive ADC Amplitudes", "fA_Pos"},
+ {"adc_1", "Raw First ADC Amplitudes", "fA_Pos_1"},
+ // {"aneg", "Raw Negative ADC Amplitudes", "fA_Neg"}, // check
+ {"apos_p", "Ped-subtracted Positive ADC Amplitudes", "fA_Pos_p"},
+ // {"aneg_p", "Ped-subtracted Negative ADC Amplitudes", "fA_Neg_p"}, // check
+ {"tpos", "Raw Positive TDC", "fT_Pos"},
+ // {"tneg", "Raw Negative TDC", "fT_Neg"}, // check
+ {"pos_npe","PEs Positive Tube","fPosNpe"},
+ // {"neg_npe","PEs Negative Tube","fNegNpe"}, // check
+ {"pos_npe_sum", "Total Positive Tube PEs", "fPosNpeSum"},
+ // {"neg_npe_sum", "Total Negative Tube PEs", "fNegNpeSum"}, // check
+ {"npe_sum", "Total PEs", "fNpeSum"},
+ {"ntdc_pos_hits", "Number of Positive Tube Hits", "fNTDCPosHits"},
+ // {"ntdc_neg_hits", "Number of Negative Tube Hits", "fNTDCNegHits"}, // check
+ {"ngood_hits", "Total number of good hits", "fNGoodHits"},
+ { 0 }
+ };
+
+ return DefineVarsFromList( vars, mode );
+}
+//_____________________________________________________________________________
+inline
+void THcCherenkov::Clear(Option_t* opt)
+{
+ // Clear the hit lists
+ fPosTDCHits->Clear();
+ // fNegTDCHits->Clear(); // check
+ fPosADCHits->Clear();
+ // fNegADCHits->Clear(); // check
+
+ // Clear Cherenkov variables from h_trans_cer.f
+
+ fNhits = 0; // Don't really need to do this. (Be sure this called before Decode)
+
+ fPosNpeSum = 0.0;
+ // fNegNpeSum = 0.0; // check
+ fNpeSum = 0.0;
+
+ fNGoodHits = 0;
+
+ fNADCPosHits = 0;
+ // fNADCNegHits = 0; // check
+ fNTDCPosHits = 0;
+ // fNTDCNegHits = 0; // check
+
+ for(Int_t itube = 0;itube < fNelem;itube++) {
+ fA_Pos[itube] = 0;
+ fA_Pos_1[itube] = 0;
+ // fA_Neg[itube] = 0; // check
+ fA_Pos_p[itube] = 0;
+ // fA_Neg_p[itube] = 0; // check
+ fT_Pos[itube] = 0;
+ // fT_Neg[itube] = 0; // check
+ fPosNpe[itube] = 0.0;
+ // fNegNpe[itube] = 0.0; // check
+ }
+
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::Decode( const THaEvData& evdata )
+{
+ // Get the Hall C style hitlist (fRawHitList) for this event
+ fNhits = THcHitList::DecodeToHitList(evdata);
+
+ if(gHaCuts->Result("Pedestal_event")) {
+
+ 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;
+ Int_t nPosTDCHits=0;
+ // Int_t nNegTDCHits=0; // check
+ Int_t nPosADCHits=0;
+ // Int_t nNegADCHits=0; // check
+ while(ihit < fNhits) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
+
+ // TDC positive hit
+ if(hit->fTDC_pos > 0) {
+ THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
+ sighit->Set(hit->fCounter, hit->fTDC_pos);
+ }
+
+ // TDC negative hit // check
+ // if(hit->fTDC_neg > 0) {
+ // THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
+ // sighit->Set(hit->fCounter, hit->fTDC_neg);
+ // }
+
+ // ADC positive hit
+ if(hit->fADC_pos > 0) {
+ THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ sighit->Set(hit->fCounter, hit->fADC_pos);
+ }
+
+ // ADC negative hit // check
+ // if(hit->fADC_neg > 0) {
+ // THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+ // sighit->Set(hit->fCounter, hit->fADC_neg);
+ // }
+
+ ihit++;
+ }
+ return ihit;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::ApplyCorrections( void )
+{
+ return(0);
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::CoarseProcess( TClonesArray& ) //tracks
+{
+ /*
+ ------------------------------------------------------------------------------------------------------------------
+ Conversion of h_trans_cer.f code:
+
+ hcer_num_hits = 0 <--- clear event
+ do tube=1,hcer_num_mirrors
+ hcer_npe(tube) = 0. <--- clear event
+ hcer_adc(tube) = 0. <--- clear event
+ enddo
+ hcer_npe_sum = 0. <--- clear event
+ do nhit = 1, hcer_tot_hits <--- loop over total hits. Very first line of this method
+ tube = hcer_tube_num(nhit) <--- tube is number of PMT on either side and it is this
+ line: Int_t npmt = hit->fCounter - 1
+ hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube) <--- This is done above:
+ fA_Pos_p[npmt] = hit->fADC_pos - fPosPedMean[npmt];
+ fA_Neg_p[npmt] = hit->fADC_neg - fNegPedMean[npmt];
+ if (hcer_adc(tube) .gt. hcer_width(tube)) then <--- This needs to convert in hcana
+ hcer_num_hits = hcer_num_hits + 1
+ hcer_tube_num(hcer_num_hits) = tube
+ hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
+ hcer_npe_sum = hcer_npe_sum + hcer_npe(tube)
+ endif
+ enddo
+ ------------------------------------------------------------------------------------------------------------------
+ */
+ for(Int_t ihit=0; ihit < fNhits; ihit++) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
+
+ // Pedestal subtraction and gain adjustment
+
+ // An ADC value of less than zero occurs when that particular
+ // channel has been sparsified away and has not been read.
+ // The NPE for that tube will be assigned zero by this code.
+ // An ADC value of greater than 8192 occurs when the ADC overflows on
+ // an input that is too large. Tubes with this characteristic will
+ // be assigned NPE = 100.0.
+
+ Int_t npmt = hit->fCounter - 1; // tube = hcer_tube_num(nhit)
+ // Should probably check that npmt is in range
+
+ cout << "fNhits = " << fNhits << endl;
+ cout << "ihit = " << ihit << endl;
+ cout << "hit->fCounter = " << hit->fCounter << endl;
+ cout << "hit->fADC_pos = " << hit->fADC_pos << endl;
+ cout << "fPosPedMean[npmt] = " << fPosPedMean[0] << endl;
+ if ( ihit == 1 )
+ cout << "--------------------------" << endl;
+
+ fA_Pos[npmt] = hit->fADC_pos;
+ // fA_Neg[npmt] = hit->fADC_neg; // check
+ fA_Pos_p[npmt] = hit->fADC_pos - fPosPedMean[npmt]; // hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube)
+ // fA_Neg_p[npmt] = hit->fADC_neg - fNegPedMean[npmt]; // hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube) // check
+ fT_Pos[npmt] = hit->fTDC_pos;
+ // fT_Neg[npmt] = hit->fTDC_neg; // check
+
+ if( ( fA_Pos_p[npmt] > 50 ) && ( hit->fADC_pos < 8000 ) ) { // if (hcer_adc(tube) .gt. hcer_width(tube)) then
+ fNADCPosHits++; // hcer_num_hits = hcer_num_hits + 1
+ fNGoodHits++;
+ fPosNpe[npmt] = fPosGain[npmt]*fA_Pos_p[npmt]; // hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
+ } else if ( hit->fADC_pos > 8000 ) {
+ fPosNpe[npmt] = 100.0;
+ } else {
+ fPosNpe[npmt] = 0.0;
+ }
+
+ if ( ihit == 0 )
+ {
+
+ fA_Pos_1[npmt] = hit->fADC_pos;
+
+ }
+
+ // check
+ // if( ( fA_Neg_p[npmt] > 50 ) && ( hit->fADC_neg < 8000 ) ) { // if (hcer_adc(tube) .gt. hcer_width(tube)) then
+ // fNADCNegHits++; // hcer_num_hits = hcer_num_hits + 1
+ // fNGoodHits++;
+ // fNegNpe[npmt] = fNegGain[npmt]*fA_Neg_p[npmt]; // hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
+ // } else if ( hit->fADC_neg > 8000 ) {
+ // fNegNpe[npmt] = 100.0;
+ // } else {
+ // fNegNpe[npmt] = 0.0;
+ // }
+
+ fPosNpeSum += fPosNpe[npmt]; // hcer_npe_sum = hcer_npe_sum + hcer_npe(tube)
+ // fNegNpeSum += fNegNpe[npmt]; // hcer_npe_sum = hcer_npe_sum + hcer_npe(tube) // check
+
+ // Sum positive and negative hits to fill tot_good_hits
+ // if(fPosNpe[npmt] > 0.3) {
+ // fNADCPosHits++;
+ // fNGoodHits++;
+ // }
+ // if(fNegNpe[npmt] > 0.3) {
+ // fNADCNegHits++;
+ // fNGoodHits++;
+ // }
+
+ if(hit->fTDC_pos > 0 && hit->fTDC_pos < 8000) {
+ fNTDCPosHits++;
+ }
+ // check
+ // if(hit->fTDC_neg > 0 && hit->fTDC_neg < 8000) {
+ // fNTDCNegHits++;
+ // }
+
+ // Fill raw adc variables with actual tubve values
+ // mainly for diagnostic purposes
+
+ }
+
+ // check
+ // if(fPosNpeSum > 0.5 || fNegNpeSum > 0.5) {
+ // fNpeSum = fPosNpeSum + fNegNpeSum;
+ // } else {
+ // fNpeSum = 0.0;
+ // }
+
+ // Ahmed
+ if(fPosNpeSum > 0.5 ) {
+ fNpeSum = fPosNpeSum;
+ } else {
+ fNpeSum = 0.0;
+ }
+
+ // If total hits are 0, then give a noticable ridiculous NPE
+ if(fNhits < 1) {
+ fNpeSum = 0.0;
+ }
+
+ // The following code is in the fortran. It probably doesn't work
+ // right because the arrays are not cleared first and the aero_ep,
+ // aero_en, ... lines make no sense.
+
+ //* Next, fill the rawadc variables with the actual tube values
+ //* mainly for diagnostic purposes.
+ //
+ // do ihit=1,haero_tot_hits
+ //
+ // npmt=haero_pair_num(ihit)
+ //
+ // haero_rawadc_pos(npmt)=haero_adc_pos(ihit)
+ // aero_ep(npmt)=haero_rawadc_pos(ihit)
+ //
+ // haero_rawadc_neg(npmt)=haero_adc_neg(ihit)
+ // aero_en(npmt)=haero_rawadc_neg(ihit)
+ //
+ // haero_rawtdc_neg(npmt)=haero_tdc_neg(ihit)
+ // aero_tn(npmt)= haero_tdc_neg(ihit)
+ //
+ // haero_rawtdc_pos(npmt)=haero_tdc_pos(ihit)
+ // aero_tp(npmt)= haero_tdc_pos(ihit)
+ //
+ // enddo
+
+ ApplyCorrections();
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::FineProcess( TClonesArray& tracks )
+{
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::InitializePedestals( )
+{
+ fNPedestalEvents = 0;
+ fMinPeds = 500; // In engine, this is set in parameter file
+ fPosPedSum = new Int_t [fNelem];
+ fPosPedSum2 = new Int_t [fNelem];
+ fPosPedLimit = new Int_t [fNelem];
+ fPosPedCount = new Int_t [fNelem];
+ // fNegPedSum = new Int_t [fNelem]; // check
+ // fNegPedSum2 = new Int_t [fNelem]; // check
+ // fNegPedLimit = new Int_t [fNelem]; // check
+ // fNegPedCount = new Int_t [fNelem]; // check
+
+ fPosPed = new Double_t [fNelem];
+ // fNegPed = new Double_t [fNelem]; // check
+ fPosThresh = new Double_t [fNelem];
+ // fNegThresh = new Double_t [fNelem]; // check
+ for(Int_t i=0;i<fNelem;i++) {
+ fPosPedSum[i] = 0;
+ fPosPedSum2[i] = 0;
+ fPosPedLimit[i] = 1000; // In engine, this are set in parameter file
+ fPosPedCount[i] = 0;
+ // fNegPedSum[i] = 0; // check
+ // fNegPedSum2[i] = 0; // check
+ // fNegPedLimit[i] = 1000; // check // In engine, this are set in parameter file
+ // fNegPedCount[i] = 0; // check
+ }
+
+ fPosNpe = new Double_t [fNelem];
+ // fNegNpe = new Double_t [fNelem]; // check
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
+{
+ // Extract data from the hit list, accumulating into arrays for
+ // calculating pedestals
+
+ Int_t nrawhits = rawhits->GetLast()+1;
+
+ Int_t ihit = 0;
+ while(ihit < nrawhits) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) rawhits->At(ihit);
+
+ Int_t element = hit->fCounter - 1;
+ Int_t adcpos = hit->fADC_pos;
+ // Int_t adcneg = hit->fADC_neg; // check
+ if(adcpos <= fPosPedLimit[element]) {
+ fPosPedSum[element] += adcpos;
+ fPosPedSum2[element] += adcpos*adcpos;
+ fPosPedCount[element]++;
+ if(fPosPedCount[element] == fMinPeds/5) {
+ fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
+ }
+ }
+ // check
+ // if(adcneg <= fNegPedLimit[element]) {
+ // fNegPedSum[element] += adcneg;
+ // fNegPedSum2[element] += adcneg*adcneg;
+ // fNegPedCount[element]++;
+ // if(fNegPedCount[element] == fMinPeds/5) {
+ // fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
+ // }
+ // }
+ ihit++;
+ }
+
+ fNPedestalEvents++;
+
+ return;
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::CalculatePedestals( )
+{
+ // Use the accumulated pedestal data to calculate pedestals
+ // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
+ // cout << "Plane: " << fPlaneNum << endl;
+ for(Int_t i=0; i<fNelem;i++) {
+
+ // Positive tubes
+ fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
+ fPosThresh[i] = fPosPed[i] + 15;
+
+ // check // Negative tubes
+ // fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
+ // fNegThresh[i] = fNegPed[i] + 15;
+
+ // cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
+
+ // Just a copy for now, but allow the possibility that fXXXPedMean is set
+ // in a parameter file and only overwritten if there is a sufficient number of
+ // pedestal events. (So that pedestals are sensible even if the pedestal events were
+ // not acquired.)
+ if(fMinPeds > 0) {
+ if(fPosPedCount[i] > fMinPeds) {
+ fPosPedMean[i] = fPosPed[i];
+ }
+ // check
+ // if(fNegPedCount[i] > fMinPeds) {
+ // fNegPedMean[i] = fNegPed[i];
+ // }
+ }
+ }
+ // cout << " " << endl;
+
+}
+void THcCherenkov::Print( const Option_t* opt) const {
+ THaNonTrackingDetector::Print(opt);
+
+ // Print out the pedestals
+
+ cout << endl;
+ cout << "Cherenkov Pedestals" << endl;
+
+ // check
+ // cout << "No. Neg Pos" << endl;
+ // for(Int_t i=0; i<fNelem; i++){
+ // cout << " " << i << " " << fNegPed[i] << " " << fPosPed[i] << endl;
+ // }
+
+ // Ahmed
+ cout << "No. Pos" << endl;
+ for(Int_t i=0; i<fNelem; i++){
+ cout << " " << i << " " << fPosPed[i] << endl;
+ }
+
+ cout << endl;
+}
+
+ClassImp(THcCherenkov)
+////////////////////////////////////////////////////////////////////////////////
+
diff --git a/quick.csh b/quick.csh
new file mode 100755
index 0000000000000000000000000000000000000000..9da9cc21cc474ff9031611cd8a23a1000c4fbfee
--- /dev/null
+++ b/quick.csh
@@ -0,0 +1,11 @@
+make
+
+cd examples
+
+../hcana <<EOF
+.X hodtest.C
+.q
+EOF
+
+cd ..
+
diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c486b2a2108852e7444590d093a0b8a12dd176a3
--- /dev/null
+++ b/src/THcCherenkov.cxx
@@ -0,0 +1,412 @@
+///////////////////////////////////////////////////////////////////////////////////////
+// //
+// THcCherenkov //
+// //
+// Class for an Cherenkov detector consisting of onw pair of PMT's //
+// //
+// Zafar Ahmed. Second attempt. November 14 2013. //
+// Comment:No need to cahnge the Map file but may need to change the parameter file //
+// //
+// This code is written for following variables: //
+// //
+// Variable Name Description //
+// adc_1 Raw value of adc 1 //
+// adc_2 Raw value of adc 2 //
+// adc_p_1 Pedestal substracted value of adc 1 //
+// adc_p_2 Pedestal substracted value of adc 2 //
+// npe_1 Number of photo electrons of adc 1 //
+// npe_2 Number of photo electrons of adc 2 //
+// //
+///////////////////////////////////////////////////////////////////////////////////////
+
+#include "THcCherenkov.h"
+#include "TClonesArray.h"
+#include "THcSignalHit.h"
+#include "THaEvData.h"
+#include "THaDetMap.h"
+#include "THcDetectorMap.h"
+#include "THcGlobals.h"
+#include "THaCutList.h"
+#include "THcParmList.h"
+#include "THcHitList.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;
+
+using std::cout;
+using std::cin;
+using std::endl;
+
+#include <iomanip>
+using std::setw;
+using std::setprecision;
+
+//_____________________________________________________________________________
+THcCherenkov::THcCherenkov( const char* name, const char* description,
+ THaApparatus* apparatus ) :
+ THaNonTrackingDetector(name,description,apparatus)
+{
+ // Normal constructor with name and description
+ fPosADCHits = new TClonesArray("THcSignalHit",16);
+
+// fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+}
+
+//_____________________________________________________________________________
+THcCherenkov::THcCherenkov( ) :
+ THaNonTrackingDetector()
+{
+ // Constructor
+}
+
+//_____________________________________________________________________________
+THcCherenkov::~THcCherenkov()
+{
+ // Destructor
+ delete [] fPosGain;
+ delete [] fPosPedLimit;
+ delete [] fPosPedMean;
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
+{
+ static const char* const here = "Init()";
+
+ cout << "THcCherenkov::Init " << GetName() << endl;
+
+ // Should probably put this in ReadDatabase as we will know the
+ // maximum number of hits after setting up the detector map
+ THcHitList::InitHitList(fDetMap, "THcCherenkovHit", 100); // 100 is max hits
+
+ EStatus status;
+ if( (status = THaNonTrackingDetector::Init( date )) )
+ return fStatus=status;
+
+ // Will need to determine which apparatus it belongs to and use the
+ // appropriate detector ID in the FillMap call
+ if( gHcDetectorMap->FillMap(fDetMap, "HCER") < 0 ) {
+ Error( Here(here), "Error filling detectormap for %s.",
+ "HCER");
+ return kInitError;
+ }
+
+ return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::ReadDatabase( const TDatime& date )
+{
+ // This function is called by THaDetectorBase::Init() once at the beginning
+ // of the analysis.
+
+ cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
+
+ char prefix[2];
+
+ prefix[0]=tolower(GetApparatus()->GetName()[0]);
+ prefix[1]='\0';
+
+ fNelem = 2; // Default if not defined
+
+ fPosGain = new Double_t[fNelem];
+ fPosPedLimit = new Int_t[fNelem];
+ fPosPedMean = new Double_t[fNelem];
+
+ DBRequest list[]={
+ {"cer_adc_to_npe", fPosGain, kDouble, fNelem}, // Ahmed
+ {"cer_ped_limit", fPosPedLimit, kInt, fNelem}, // Ahmed
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ fIsInit = true;
+
+ // Create arrays to hold pedestal results
+ InitializePedestals();
+
+ return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::DefineVariables( EMode mode )
+{
+ // Initialize global variables for histogramming and tree
+
+ cout << "THcCherenkov::DefineVariables called " << GetName() << endl;
+
+ if( mode == kDefine && fIsSetup ) return kOK;
+ fIsSetup = ( mode == kDefine );
+
+ // Register variables in global list
+
+ // Do we need to put the number of pos/neg TDC/ADC hits into the variables?
+ // No. They show up in tree as Ndata.H.aero.postdchits for example
+
+ RVarDef vars[] = {
+ {"adc_1", "Raw First ADC Amplitudes", "fA_1"},
+ {"adc_2", "Raw Second ADC Amplitudes", "fA_2"},
+ {"adc_p_1", "Raw First ADC Amplitudes", "fA_p_1"},
+ {"adc_p_2", "Raw Second ADC Amplitudes", "fA_p_2"},
+ {"npe_1","PEs First Tube", "fNpe_1"},
+ {"npe_2","PEs Second Tube","fNpe_2"},
+ {"posadchits", "List of Positive ADC hits","fPosADCHits.THcSignalHit.GetPaddleNumber()"},
+ { 0 }
+ };
+
+ return DefineVarsFromList( vars, mode );
+}
+//_____________________________________________________________________________
+inline
+void THcCherenkov::Clear(Option_t* opt)
+{
+ // Clear the hit lists
+ fPosADCHits->Clear();
+
+ // Clear Cherenkov variables from h_trans_cer.f
+
+ fNhits = 0; // Don't really need to do this. (Be sure this called before Decode)
+
+ fA_1 = 0;
+ fA_2 = 0;
+ fA_p_1 = 0;
+ fA_p_2 = 0;
+ fNpe_1 = 0;
+ fNpe_2 = 0;
+
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::Decode( const THaEvData& evdata )
+{
+ // Get the Hall C style hitlist (fRawHitList) for this event
+ fNhits = THcHitList::DecodeToHitList(evdata);
+
+ if(gHaCuts->Result("Pedestal_event")) {
+
+ 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;
+ Int_t nPosADCHits=0;
+ while(ihit < fNhits) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
+
+ // ADC positive hit
+ if(hit->fADC_pos > 0) {
+ THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ sighit->Set(hit->fCounter, hit->fADC_pos);
+ }
+
+ ihit++;
+ }
+ return ihit;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::ApplyCorrections( void )
+{
+ return(0);
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::CoarseProcess( TClonesArray& ) //tracks
+{
+ /*
+ ------------------------------------------------------------------------------------------------------------------
+ h_trans_cer.f code:
+
+ hcer_num_hits = 0 <--- clear event
+ do tube=1,hcer_num_mirrors
+ hcer_npe(tube) = 0. <--- clear event
+ hcer_adc(tube) = 0. <--- clear event
+ enddo
+ hcer_npe_sum = 0. <--- clear event
+ do nhit = 1, hcer_tot_hits <--- loop over total hits. Very first line of this method
+ tube = hcer_tube_num(nhit) <--- tube is number of PMT on either side and it is this
+ line: Int_t npmt = hit->fCounter - 1
+ hcer_adc(tube) = hcer_raw_adc(nhit) - hcer_ped(tube) <--- This is done above:
+ fA_Pos_p[npmt] = hit->fADC_pos - fPosPedMean[npmt];
+ fA_Neg_p[npmt] = hit->fADC_neg - fNegPedMean[npmt];
+ if (hcer_adc(tube) .gt. hcer_width(tube)) then <--- This needs to convert in hcana
+ hcer_num_hits = hcer_num_hits + 1
+ hcer_tube_num(hcer_num_hits) = tube
+ hcer_npe(tube) = hcer_adc(tube) * hcer_adc_to_npe(tube)
+ hcer_npe_sum = hcer_npe_sum + hcer_npe(tube)
+ endif
+ enddo
+ ------------------------------------------------------------------------------------------------------------------
+ */
+ for(Int_t ihit=0; ihit < fNhits; ihit++) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
+
+ // Pedestal subtraction and gain adjustment
+
+ // An ADC value of less than zero occurs when that particular
+ // channel has been sparsified away and has not been read.
+ // The NPE for that tube will be assigned zero by this code.
+ // An ADC value of greater than 8192 occurs when the ADC overflows on
+ // an input that is too large. Tubes with this characteristic will
+ // be assigned NPE = 100.0.
+ Int_t npmt = hit->fCounter - 1; // tube = hcer_tube_num(nhit)
+ // Should probably check that npmt is in range
+ if ( ihit != npmt )
+ cout << "ihit != npmt " << endl;
+
+ if ( npmt == 0 ) {
+ fA_1 = hit->fADC_pos;
+ fA_p_1 = hit->fADC_pos - fPosPedMean[npmt];
+ if ( ( fA_p_1 > 50 ) && ( hit->fADC_pos < 8000 ) ) {
+ fNpe_1 = fPosGain[npmt]*fA_p_1;
+ } else if ( hit->fADC_pos > 8000 ) {
+ fNpe_1 = 100.0;
+ } else {
+ fNpe_1 = 0.0;
+ }
+ }
+
+ if ( npmt == 1 ) {
+ fA_2 = hit->fADC_pos;
+ fA_p_2 = hit->fADC_pos - fPosPedMean[npmt];
+ if ( ( fA_p_2 > 50 ) && ( hit->fADC_pos < 8000 ) ) {
+ fNpe_2 = fPosGain[npmt]*fA_p_2;
+ } else if ( hit->fADC_pos > 8000 ) {
+ fNpe_2 = 100.0;
+ } else {
+ fNpe_2 = 0.0;
+ }
+ }
+ }
+
+ ApplyCorrections();
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcCherenkov::FineProcess( TClonesArray& tracks )
+{
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::InitializePedestals( )
+{
+ fNPedestalEvents = 0;
+ fMinPeds = 500; // In engine, this is set in parameter file
+ fPosPedSum = new Int_t [fNelem];
+ fPosPedSum2 = new Int_t [fNelem];
+ fPosPedLimit = new Int_t [fNelem];
+ fPosPedCount = new Int_t [fNelem];
+
+ fPosPed = new Double_t [fNelem];
+ fPosThresh = new Double_t [fNelem];
+ for(Int_t i=0;i<fNelem;i++) {
+ fPosPedSum[i] = 0;
+ fPosPedSum2[i] = 0;
+ fPosPedLimit[i] = 1000; // In engine, this are set in parameter file
+ fPosPedCount[i] = 0;
+ }
+
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
+{
+ // Extract data from the hit list, accumulating into arrays for
+ // calculating pedestals
+
+ Int_t nrawhits = rawhits->GetLast()+1;
+
+ Int_t ihit = 0;
+ while(ihit < nrawhits) {
+ THcCherenkovHit* hit = (THcCherenkovHit *) rawhits->At(ihit);
+
+ Int_t element = hit->fCounter - 1;
+ Int_t adcpos = hit->fADC_pos;
+ if(adcpos <= fPosPedLimit[element]) {
+ fPosPedSum[element] += adcpos;
+ fPosPedSum2[element] += adcpos*adcpos;
+ fPosPedCount[element]++;
+ if(fPosPedCount[element] == fMinPeds/5) {
+ fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
+ }
+ }
+ ihit++;
+ }
+
+ fNPedestalEvents++;
+
+ return;
+}
+
+//_____________________________________________________________________________
+void THcCherenkov::CalculatePedestals( )
+{
+ // Use the accumulated pedestal data to calculate pedestals
+ // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
+ // cout << "Plane: " << fPlaneNum << endl;
+ for(Int_t i=0; i<fNelem;i++) {
+
+ // Positive tubes
+ fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
+ fPosThresh[i] = fPosPed[i] + 15;
+
+ // cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
+
+ // Just a copy for now, but allow the possibility that fXXXPedMean is set
+ // in a parameter file and only overwritten if there is a sufficient number of
+ // pedestal events. (So that pedestals are sensible even if the pedestal events were
+ // not acquired.)
+ if(fMinPeds > 0) {
+ if(fPosPedCount[i] > fMinPeds) {
+ fPosPedMean[i] = fPosPed[i];
+ }
+ }
+ }
+ // cout << " " << endl;
+
+}
+void THcCherenkov::Print( const Option_t* opt) const {
+ THaNonTrackingDetector::Print(opt);
+
+ // Print out the pedestals
+
+ cout << endl;
+ cout << "Cherenkov Pedestals" << endl;
+
+ // Ahmed
+ cout << "No. Pos" << endl;
+ for(Int_t i=0; i<fNelem; i++){
+ cout << " " << i << " " << fPosPed[i] << endl;
+ }
+
+ cout << endl;
+}
+
+ClassImp(THcCherenkov)
+////////////////////////////////////////////////////////////////////////////////
+
diff --git a/src/THcCherenkov.h b/src/THcCherenkov.h
new file mode 100644
index 0000000000000000000000000000000000000000..135a8b6705f26706096b8dc9806f90df460272ed
--- /dev/null
+++ b/src/THcCherenkov.h
@@ -0,0 +1,111 @@
+#ifndef ROOT_THcCherenkov
+#define ROOT_THcCherenkov
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Cherenkov //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "TClonesArray.h"
+#include "THaNonTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcCherenkovHit.h"
+
+class THcCherenkov : public THaNonTrackingDetector, public THcHitList {
+
+ public:
+ THcCherenkov( const char* name, const char* description = "",
+ THaApparatus* a = NULL );
+ virtual ~THcCherenkov();
+
+ virtual void Clear( Option_t* opt="" );
+ virtual Int_t Decode( const THaEvData& );
+ virtual EStatus Init( const TDatime& run_time );
+ virtual Int_t ReadDatabase( const TDatime& date );
+ virtual Int_t DefineVariables( EMode mode = kDefine );
+ virtual Int_t CoarseProcess( TClonesArray& tracks );
+ virtual Int_t FineProcess( TClonesArray& tracks );
+
+ virtual void AccumulatePedestals(TClonesArray* rawhits);
+ virtual void CalculatePedestals();
+
+ virtual Int_t ApplyCorrections( void );
+
+ virtual void Print(const Option_t* opt) const;
+
+ THcCherenkov(); // for ROOT I/O
+ protected:
+ Int_t fAnalyzePedestals;
+
+ // Parameters
+ Double_t* fPosGain;
+ Double_t* fNegGain;
+
+ // Event information
+ Int_t fNhits;
+
+ Float_t* fA_Pos; // [fNelem] Array of ADC amplitudes
+ Float_t* fA_Neg; // [fNelem] Array of ADC amplitudes
+ Float_t* fA_Pos_p; // [fNelem] Array of ped-subtracted ADC amplitudes
+ Float_t* fA_Neg_p; // [fNelem] Array of ped-subtracted ADC amplitudes
+ Float_t* fT_Pos; // [fNelem] Array of TDCs
+ Float_t* fT_Neg; // [fNelem] Array of TDCs
+
+ Float_t fA_1; // Ahmed
+ Float_t fA_2; // Ahmed
+ Float_t fA_p_1; // Ahmed
+ Float_t fA_p_2; // Ahmed
+ Double_t fNpe_1; // Ahmed
+ Double_t fNpe_2; // Ahmed
+ Int_t fNHits_1; // Ahmed
+ Int_t fNHits_2; // Ahmed
+
+ Double_t fPosNpeSum;
+ Double_t fNegNpeSum;
+ Double_t fNpeSum;
+ Int_t fNGoodHits;
+ Int_t fNADCPosHits;
+ Int_t fNADCNegHits;
+ Int_t fNTDCPosHits;
+ Int_t fNTDCNegHits;
+
+ Double_t* fPosNpe; // [fNelem] # Photoelectrons per positive tube
+ Double_t* fNegNpe; // [fNelem] # Photoelectrons per negative tube
+
+
+ // Hits
+ TClonesArray* fPosTDCHits;
+ TClonesArray* fNegTDCHits;
+ TClonesArray* fPosADCHits;
+ TClonesArray* fNegADCHits;
+
+ // Pedestals
+ Int_t fNPedestalEvents;
+ Int_t fMinPeds;
+ Int_t *fPosPedSum; /* Accumulators for pedestals */
+ Int_t *fPosPedSum2;
+ Int_t *fPosPedLimit;
+ Int_t *fPosPedCount;
+ Int_t *fNegPedSum;
+ Int_t *fNegPedSum2;
+ Int_t *fNegPedLimit;
+ Int_t *fNegPedCount;
+
+ Double_t *fPosPed;
+ Double_t *fPosSig;
+ Double_t *fPosThresh;
+ Double_t *fNegPed;
+ Double_t *fNegSig;
+ Double_t *fNegThresh;
+
+ Double_t *fPosPedMean; /* Can be supplied in parameters and then */
+ Double_t *fNegPedMean; /* be overwritten from ped analysis */
+
+ void Setup(const char* name, const char* description);
+ virtual void InitializePedestals( );
+
+ ClassDef(THcCherenkov,0) // Generic cherenkov class
+};
+
+#endif
diff --git a/src/THcCherenkovHit.cxx b/src/THcCherenkovHit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c946328351e36c814fe88dc672b120420ae0e5ec
--- /dev/null
+++ b/src/THcCherenkovHit.cxx
@@ -0,0 +1,19 @@
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcCherenkovHit //
+// //
+// Class representing a single raw hit for a pair of aerogel tubes //
+// //
+// Contains plane, counter and pos/neg adc //
+// //
+// Assumes Aerogel design where Aerogel PMT's are in pairs and only have //
+// ADCs. //
+// //
+// Ahmed. First attempt. September 27 2013. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcCherenkovHit.h"
+
+//////////////////////////////////////////////////////////////////////////
+ClassImp(THcCherenkovHit)
diff --git a/src/THcCherenkovHit.h b/src/THcCherenkovHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..67ef4c13c1e0a344a06803a581d81b869d79e426
--- /dev/null
+++ b/src/THcCherenkovHit.h
@@ -0,0 +1,18 @@
+#ifndef ROOT_THcCherenkovHit
+#define ROOT_THcCherenkovHit
+
+#include "THcHodoscopeHit.h"
+
+class THcCherenkovHit : public THcHodoscopeHit {
+
+ public:
+
+ protected:
+
+ private:
+
+ ClassDef(THcCherenkovHit,0); // Cherenkov hit class
+};
+
+#endif
+