diff --git a/Makefile b/Makefile
index fadfb2944800f96524642667e412b4977efa1e24..422e44b8456bc4d1847309026d5fb6da4dcf2d4b 100644
--- a/Makefile
+++ b/Makefile
@@ -15,7 +15,8 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
src/THcSignalHit.cxx \
src/THcHodoscope.cxx src/THcScintillatorPlane.cxx \
src/THcHodoscopeHit.cxx \
- src/THcDriftChamber.cxx src/THcDriftChamberPlane.cxx \
+ src/THcDC.cxx src/THcDriftChamberPlane.cxx \
+ src/THcDriftChamber.cxx \
src/THcRawDCHit.cxx src/THcDCHit.cxx \
src/THcDCWire.cxx \
src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \
diff --git a/examples/PARAM/hcana.param b/examples/PARAM/hcana.param
index d7f6f97b403c0daeace562384cbacb2285482767..61d60c67416cd2b92a6e69416e20b1281f30689b 100644
--- a/examples/PARAM/hcana.param
+++ b/examples/PARAM/hcana.param
@@ -14,3 +14,4 @@ hdc_plane_names = "1x1 1y1 1u1 1v1 1y2 1x2 2x1 2y1 2u1 2v1 2y2 2x2"
hcal_layer_names = "1pr 2ta 3ta 4ta"
+hhodo_plane_names = "1x 1y 2x 2y"
diff --git a/examples/PARAM/hhodo.param b/examples/PARAM/hhodo.param
index 8a4ba606916588424553151b44d4ccb199bea081..4b7fcf89a862e88c022473c458f697227bee964c 100644
--- a/examples/PARAM/hhodo.param
+++ b/examples/PARAM/hhodo.param
@@ -11,7 +11,7 @@
; new variable for picking good hits for tof fitting
; this should not be set tight until you are ready to fit
; tof and you figured out good values
-htof_tolerance = 50.0
+ htof_tolerance = 3.0
;
; hms_tof_params
; hnum_scin_counters, hhodo_zpos, hhodo_center_coord, hhodo_width
@@ -127,38 +127,38 @@ htof_tolerance = 50.0
; csa 9/8/98 -- I had to hand-twaddle a few of the values
; based on (relative) offsets of older hhodo.param
-hhodo_pos_time_offset = 5.5391, 0.0000, -6.8586, 5.1830
- -1.2985, -0.2551, 0.3953, 8.9369
- 5.0531, 1.4333, -6.5520, 8.0628
- 5.5481, -0.3965, -3.6470, 12.8933
- 1.2416, -7.0707, 9.6246, 10.0208
- 1.9282, -0.9275, -4.1448, 8.2112
- 6.4910, -1.5780, 1.8980, 9.6928
- 4.4770, 1.7009, -3.8385, 4.7545
- 2.7533, -2.7182, -5.8864, 6.3882
- 4.3398, 6.1158, -0.3572, -0.4308
- 6.0782, 0.0000, -7.5343, 0.0000
- 5.4665, 0.0000, 0.2169, 0.0000
- 4.1334, 0.0000, 1.3767, 0.0000
- 1.6088, 0.0000, 2.5930, 0.0000
- 3.9776, 0.0000, -5.0340, 0.0000
- 1.6534, 0.0000, 1.5043, 0.0000
-hhodo_neg_time_offset = -2.5728, 0.0000, 2.8982, 10.8670
- -1.9187, -2.6479, 10.6272, 13.8790
- -4.1126, -4.5084, 3.9705, 15.5799
- 0.7699, -2.3908, 11.7183, 15.1612
- -1.2568, -5.0343, 2.9473, 10.6625
- -2.8197, 0.7670, 10.3919, 7.8739
- 1.2798, -3.9185, 4.3248, 11.3533
- -4.8009, -0.2453, 9.2837, 11.6355
- -6.3004, -5.7362, 0.8352, 14.9451
- 1.8476, 6.1209, 11.9751, 15.7375
- 0.3913, 0.0000, 8.9105, 0.0000
- -1.0702, 0.0000, 9.8926, 0.0000
- -2.3617, 0.0000, 5.7061, 0.0000
- -5.2931, 0.0000, 10.7318, 0.0000
- 0.0632, 0.0000, 6.6962, 0.0000
- 6.7208, 0.0000, 13.4108, 0.0000
+hhodo_pos_time_offset = 6.5821, 0.0000, -6.9727, 4.4189
+ -1.1011, 0.1270, 0.1930, 8.2907
+ 5.5928, 1.2312, -7.0087, 8.3756
+ 5.3494, -0.2264, -3.9023, 13.4700
+ 1.6938, -7.1708, 9.2976, 9.9767
+ 1.7531, -0.9044, -4.6540, 7.8608
+ 6.9700, -1.5083, 1.7866, 8.6117
+ 3.3819, 1.8281, -4.2414, 4.3756
+ 3.4984, -2.5534, -6.0882, 5.4043
+ 4.5793, 6.1158, -0.5954, -1.3631
+ 6.8272, 0.0000, -7.6493, 32.1266
+ 5.4571, 2.0212, 0.2085, 0.0000
+ 5.0147, 2.0219, 1.1072, 32.0001
+ 1.4959, 0.0000, 2.4523, 0.0000
+ 6.2602, 0.0000, -5.1102, -1.9990
+ 1.6534, 0.0000, 1.2513, 0.0000
+hhodo_neg_time_offset = -2.4664, 0.0000, 4.1095, 11.9016
+ -2.0589, -2.5527, 10.4536, 13.4114
+ -4.2564, -3.8944, 3.8712, 15.4475
+ 0.3935, -2.2300, 11.5590, 14.8406
+ -1.7111, -4.8155, 2.8606, 10.2379
+ -3.0792, 0.6972, 10.1887, 7.8083
+ 1.0046, -3.6839, 4.2680, 10.9969
+ -5.2801, 0.9309, 9.1779, 11.4677
+ -6.4931, -5.5225, 0.7246, 14.8212
+ 1.6218, 6.1209, 11.8919, 15.5699
+ 0.0703, -1.9990, 8.8625, 2.0216
+ -1.3146, 0.0000, 9.8650, 0.0000
+ -2.9002, 0.0000, 5.6735, 0.0000
+ -5.6444, 0.0000, 10.7529, 0.0000
+ 0.5662, 0.0000, 6.6559, 2.0115
+ 6.7208, 0.0000, 13.4179, 36.7325
;
hhodo_pos_ped_limit = 1000,1000,1000,1000,1000,1000,1000,1000
1000,1000,1000,1000,1000,1000,1000,1000
@@ -178,11 +178,9 @@ hhodo_neg_ped_limit = 1000,1000,1000,1000,1000,1000,1000,1000
1000,1000,1000,1000,1000,1000,1000,1000
1000,1000,1000,1000,1000,1000,1000,1000
-
; use htofusinginvadc=1 if want invadc_offset
; invadc_linear, and invadc_adc to be used
htofusinginvadc=1
-
hhodo_pos_invadc_offset = 0.00, 0.00, -7.49, 0.00
0.00, -3.93, -1.33, 6.35
-0.47, -1.79, -7.92, 6.66
@@ -251,6 +249,7 @@ hhodo_neg_invadc_linear = 50.00, 50.00, 11.14, 50.00
26.89, 50.00, 11.54, 50.00
50.00, 50.00, 10.51, 50.00
+
hhodo_pos_invadc_adc= 0.00, 0.00, 24.53, 0.00
0.00, 85.08, 40.16, 46.66
104.72, 90.39, 48.39, 40.45
@@ -284,48 +283,3 @@ hhodo_neg_invadc_adc= 0.00, 0.00, 33.27, 0.00
78.19, 0.00, 51.72, 0.00
75.03, 0.00, 49.34, 0.00
0.00, 0.00, 53.40, 0.00
-
-hhodo_pos_sigma = 100.00, 100.00, 0.92, 100.00
- 100.00, 0.55, 0.85, 0.80
- 0.79, 0.58, 0.80, 0.75
- 0.68, 0.58, 0.65, 0.79
- 0.73, 0.58, 0.66, 0.79
- 0.59, 0.58, 0.65, 0.74
- 0.59, 0.58, 0.53, 0.76
- 0.61, 0.57, 0.67, 0.80
- 0.57, 0.59, 0.64, 0.81
- 0.61, 100.00, 0.64, 0.76
- 0.67, 100.00, 0.67, 100.00
- 0.52, 100.00, 0.61, 100.00
- 0.56, 100.00, 0.62, 100.00
- 0.61, 100.00, 0.71, 100.00
- 0.70, 100.00, 0.74, 100.00
- 100.00, 100.00, 0.66, 100.00
-
-hhodo_neg_sigma = 100.00, 100.00, 0.92, 100.00
- 100.00, 0.58, 0.82, 0.80
- 0.74, 0.57, 0.81, 0.81
- 0.68, 0.58, 0.72, 0.81
- 100.00, 0.56, 0.68, 0.88
- 0.61, 0.54, 0.61, 0.76
- 0.53, 0.55, 0.53, 0.73
- 0.58, 0.55, 0.65, 0.80
- 0.61, 0.56, 0.71, 0.87
- 0.59, 100.00, 0.69, 0.88
- 100.00, 100.00, 0.65, 100.00
- 0.57, 100.00, 0.62, 100.00
- 0.58, 100.00, 0.57, 100.00
- 0.59, 100.00, 0.64, 100.00
- 0.57, 100.00, 0.84, 100.00
- 100.00, 100.00, 0.73, 100.00
-
-
-
-
-
-
-
-
-
-
-
diff --git a/examples/hodtest.C b/examples/hodtest.C
index 0a7370289f9f293c9faa1c76b3d9a7bb7ac8a568..982e284adffdc261d98dec97cbf39ec5f105b017 100644
--- a/examples/hodtest.C
+++ b/examples/hodtest.C
@@ -38,7 +38,7 @@
// Add hodoscope
HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
HMS->AddDetector( new THcShower("cal", "Shower" ));
- HMS->AddDetector( new THcDriftChamber("dc", "Drift Chambers" ));
+ HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
THcAerogel* aerogel = new THcAerogel("aero", "Aerogel Cerenkov" );
HMS->AddDetector( aerogel );
@@ -63,7 +63,7 @@
// Eventually need to learn to skip over, or properly analyze
// the pedestal events
- run->SetEventRange(1,1000000);// Physics Event number, does not
+ run->SetEventRange(1,2000);// Physics Event number, does not
// include scaler or control events
// Define the analysis parameters
diff --git a/examples/hodtest_cuts.def b/examples/hodtest_cuts.def
index be6d4297d79edd4e786eed2d543a2f3c064a2342..f6003c2200a18b2f4d7156e51c0250c8b5ce2b06 100644
--- a/examples/hodtest_cuts.def
+++ b/examples/hodtest_cuts.def
@@ -9,3 +9,9 @@ RawDecode_master 1
Block: Decode
Decode_master !Pedestal_event
+Block: CoarseTracking
+CoarseTracking_master !Pedestal_event
+
+Block: CoarseReconstruct
+RawCoarseReconstruct !Pedestal_event
+
diff --git a/examples/output.def b/examples/output.def
index 0e282d441984a9b4f648e231783988f941afd4cc..a8fb8536136ab159c7df20d0efdde3af357b6e76 100644
--- a/examples/output.def
+++ b/examples/output.def
@@ -164,3 +164,9 @@ TH1F hdc2v1_dd 'HDC 2V1 Drift Distance' H.dc.2v1.dist 300 -0.1 0.6
TH1F hdc2y2_dd 'HDC 2Y2 Drift Distance' H.dc.2y2.dist 300 -0.1 0.6
TH1F hdc2x2_dd 'HDC 2X2 Drift Distance' H.dc.2x2.dist 300 -0.1 0.6
+# Focal Plane times
+TH1F hs1xfptime 'HODO s1x fptime' H.hod.1x.fptime 80 0 80 H.hod.hgoodstarttime
+TH1F hs1yfptime 'HODO s1y fptime' H.hod.1y.fptime 80 0 80 H.hod.hgoodstarttime
+TH1F hs2xfptime 'HODO s2x fptime' H.hod.2x.fptime 80 0 80 H.hod.hgoodstarttime
+TH1F hs2yfptime 'HODO s2y fptime' H.hod.2y.fptime 80 0 80 H.hod.hgoodstarttime
+TH1F starttime 'HODO start time' H.hod.starttime 80 0 80 H.hod.hgoodstarttime
diff --git a/src/HallC_LinkDef.h b/src/HallC_LinkDef.h
index 608c0041beef3679cbdce6befbbf801ca61abd6d..ce520fa067dda33b10aebdad9bebadeb2336f719 100644
--- a/src/HallC_LinkDef.h
+++ b/src/HallC_LinkDef.h
@@ -18,6 +18,7 @@
#pragma link C++ class THcHodoscope+;
#pragma link C++ class THcScintillatorPlane+;
#pragma link C++ class THcHodoscopeHit+;
+#pragma link C++ class THcDC+;
#pragma link C++ class THcDriftChamber+;
#pragma link C++ class THcDriftChamberPlane+;
#pragma link C++ class THcRawDCHit+;
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e79868cd26090cc462a36fdb37b40d3942064365
--- /dev/null
+++ b/src/THcDC.cxx
@@ -0,0 +1,432 @@
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcDC //
+// //
+// Class for a generic hodoscope consisting of multiple //
+// planes with multiple paddles with phototubes on both ends. //
+// This differs from Hall A scintillator class in that it is the whole //
+// hodoscope array, not just one plane. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcDC.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;
+
+//_____________________________________________________________________________
+THcDC::THcDC(
+ const char* name, const char* description,
+ THaApparatus* apparatus ) :
+ THaTrackingDetector(name,description,apparatus)
+{
+ // Constructor
+
+ // fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+ fNPlanes = 0; // No planes until we make them
+
+}
+
+//_____________________________________________________________________________
+void THcDC::Setup(const char* name, const char* description)
+{
+
+ char prefix[2];
+ char parname[100];
+
+ THaApparatus *app = GetApparatus();
+ if(app) {
+ cout << app->GetName() << endl;
+ } else {
+ cout << "No apparatus found" << endl;
+ }
+
+ prefix[0]=tolower(app->GetName()[0]);
+ prefix[1]='\0';
+
+ string planenamelist;
+ DBRequest list[]={
+ {"dc_num_planes",&fNPlanes, kInt},
+ {"dc_num_chambers",&fNChambers, kInt},
+ {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
+ {"dc_wire_velocity",&fWireVelocity,kDouble},
+ {"dc_plane_names",&planenamelist, kString},
+ {0}
+ };
+
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ cout << planenamelist << endl;
+ cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
+
+ vector<string> plane_names = vsplit(planenamelist);
+
+ if(plane_names.size() != (UInt_t) fNPlanes) {
+ cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
+ // Should quit. Is there an official way to quit?
+ }
+ fPlaneNames = new char* [fNPlanes];
+ for(Int_t i=0;i<fNPlanes;i++) {
+ fPlaneNames[i] = new char[plane_names[i].length()];
+ strcpy(fPlaneNames[i], plane_names[i].c_str());
+ }
+
+ char *desc = new char[strlen(description)+100];
+ fPlanes = new THcDriftChamberPlane* [fNPlanes];
+
+ for(Int_t i=0;i<fNPlanes;i++) {
+ strcpy(desc, description);
+ strcat(desc, " Plane ");
+ strcat(desc, fPlaneNames[i]);
+
+ fPlanes[i] = new THcDriftChamberPlane(fPlaneNames[i], desc, i+1, this);
+ cout << "Created Drift Chamber Plane " << fPlaneNames[i] << ", " << desc << endl;
+
+ }
+ fChambers = new THcDriftChamber* [fNChambers];
+ for(Int_t i=0;i<fNChambers;i++) {
+ sprintf(desc,"%s Chamber %d",description, i+1);
+
+ // Should construct a better chamber name
+ fChambers[i] = new THcDriftChamber(desc, desc, i+1, this);
+ cout << "Created Drift Chamber " << i+1 << ", " << desc << endl;
+
+
+ }
+}
+
+//_____________________________________________________________________________
+THcDC::THcDC( ) :
+ THaTrackingDetector()
+{
+ // Constructor
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
+{
+ static const char* const here = "Init()";
+
+ Setup(GetName(), GetTitle()); // Create the subdetectors here
+
+ // Should probably put this in ReadDatabase as we will know the
+ // maximum number of hits after setting up the detector map
+ THcHitList::InitHitList(fDetMap, "THcRawDCHit", 1000);
+
+ EStatus status;
+ // This triggers call of ReadDatabase and DefineVariables
+ if( (status = THaTrackingDetector::Init( date )) )
+ return fStatus=status;
+
+ // Initialize planes and add them to chambers
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if((status = fPlanes[ip]->Init( date ))) {
+ return fStatus=status;
+ } else {
+ Int_t chamber=fNChamber[ip];
+ fChambers[chamber-1]->AddPlane(fPlanes[ip]);
+ }
+ }
+ // Initialize chambers
+ for(Int_t ic=0;ic<fNChambers;ic++) {
+ if((status = fChambers[ic]->Init ( date ))) {
+ return fStatus=status;
+ }
+ }
+
+ // Replace with what we need for Hall C
+ // const DataDest tmp[NDEST] = {
+ // { &fRTNhit, &fRANhit, fRT, fRT_c, fRA, fRA_p, fRA_c, fROff, fRPed, fRGain },
+ // { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
+ // };
+ // memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
+
+ // Will need to determine which apparatus it belongs to and use the
+ // appropriate detector ID in the FillMap call
+ char EngineDID[4];
+
+ EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
+ EngineDID[1] = 'D';
+ EngineDID[2] = 'C';
+ EngineDID[3] = '\0';
+
+ if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
+ Error( Here(here), "Error filling detectormap for %s.",
+ EngineDID);
+ return kInitError;
+ }
+
+ return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::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()";
+ char prefix[2];
+ char parname[100];
+
+ // Read data from database
+ // Pull values from the THcParmList instead of reading a database
+ // file like Hall A does.
+
+ // 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)
+
+ prefix[0]=tolower(GetApparatus()->GetName()[0]);
+
+ prefix[1]='\0';
+
+ fXCenter = new Double_t [fNChambers];
+ fYCenter = new Double_t [fNChambers];
+ fMinHits = new Int_t [fNChambers];
+ fMaxHits = new Int_t [fNChambers];
+ fMinCombos = new Int_t [fNChambers];
+ fSpace_Point_Criterion2 = new Double_t [fNChambers];
+
+ fTdcWinMin = new Int_t [fNPlanes];
+ fTdcWinMax = new Int_t [fNPlanes];
+ fCentralTime = new Int_t [fNPlanes];
+ fNWires = new Int_t [fNPlanes];
+ fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
+ fWireOrder = new Int_t [fNPlanes]; // Wire readout order
+ fDriftTimeSign = new Int_t [fNPlanes];
+
+ fZPos = new Double_t [fNPlanes];
+ fAlphaAngle = new Double_t [fNPlanes];
+ fBetaAngle = new Double_t [fNPlanes];
+ fGammaAngle = new Double_t [fNPlanes];
+ fPitch = new Double_t [fNPlanes];
+ fCentralWire = new Double_t [fNPlanes];
+ fPlaneTimeZero = new Double_t [fNPlanes];
+ fSigma = new Double_t [fNPlanes];
+
+ DBRequest list[]={
+ {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
+ {"dc_wire_velocity",&fWireVelocity,kDouble},
+
+ {"dc_xcenter", fXCenter, kDouble, fNChambers},
+ {"dc_ycenter", fYCenter, kDouble, fNChambers},
+ {"min_hit", fMinHits, kInt, fNChambers},
+ {"max_pr_hits", fMaxHits, kInt, fNChambers},
+ {"min_combos", fMinCombos, kInt, fNChambers},
+ {"space_point_criterion", fSpace_Point_Criterion2, kDouble, fNChambers},
+
+ {"dc_tdc_min_win", fTdcWinMin, kInt, fNPlanes},
+ {"dc_tdc_max_win", fTdcWinMax, kInt, fNPlanes},
+ {"dc_central_time", fCentralTime, kInt, fNPlanes},
+ {"dc_nrwire", fNWires, kInt, fNPlanes},
+ {"dc_chamber_planes", fNChamber, kInt, fNPlanes},
+ {"dc_wire_counting", fWireOrder, kInt, fNPlanes},
+ {"dc_drifttime_sign", fDriftTimeSign, kInt, fNPlanes},
+
+ {"dc_zpos", fZPos, kDouble, fNPlanes},
+ {"dc_alpha_angle", fAlphaAngle, kDouble, fNPlanes},
+ {"dc_beta_angle", fBetaAngle, kDouble, fNPlanes},
+ {"dc_gamma_angle", fGammaAngle, kDouble, fNPlanes},
+ {"dc_pitch", fPitch, kDouble, fNPlanes},
+ {"dc_central_wire", fCentralWire, kDouble, fNPlanes},
+ {"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
+ {"dc_sigma", fSigma, kDouble, fNPlanes},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ cout << "Plane counts:";
+ for(Int_t i=0;i<fNPlanes;i++) {
+ cout << " " << fNWires[i];
+ }
+ cout << endl;
+
+ fIsInit = true;
+
+ return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::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 DC hits", "fNhits" },
+ { 0 }
+ };
+ return DefineVarsFromList( vars, mode );
+
+}
+
+//_____________________________________________________________________________
+THcDC::~THcDC()
+{
+ // Destructor. Remove variables from global list.
+
+ if( fIsSetup )
+ RemoveVariables();
+ if( fIsInit )
+ DeleteArrays();
+ if (fTrackProj) {
+ fTrackProj->Clear();
+ delete fTrackProj; fTrackProj = 0;
+ }
+}
+
+//_____________________________________________________________________________
+void THcDC::DeleteArrays()
+{
+ // Delete member arrays. Used by destructor.
+
+ delete [] fNWires; fNWires = NULL;
+ // delete [] fSpacing; fSpacing = NULL;
+ // delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete?
+
+ // delete [] fRA_c; fRA_c = NULL;
+ // delete [] fRA_p; fRA_p = NULL;
+ // delete [] fRA; fRA = NULL;
+ // delete [] fLA_c; fLA_c = NULL;
+ // delete [] fLA_p; fLA_p = NULL;
+ // delete [] fLA; fLA = NULL;
+ // delete [] fRT_c; fRT_c = NULL;
+ // delete [] fRT; fRT = NULL;
+ // delete [] fLT_c; fLT_c = NULL;
+ // delete [] fLT; fLT = NULL;
+
+ // delete [] fRGain; fRGain = NULL;
+ // delete [] fLGain; fLGain = NULL;
+ // delete [] fRPed; fRPed = NULL;
+ // delete [] fLPed; fLPed = NULL;
+ // delete [] fROff; fROff = NULL;
+ // delete [] fLOff; fLOff = NULL;
+ // delete [] fTWalkPar; fTWalkPar = NULL;
+ // delete [] fTrigOff; fTrigOff = NULL;
+
+ // delete [] fHitPad; fHitPad = NULL;
+ // delete [] fTime; fTime = NULL;
+ // delete [] fdTime; fdTime = NULL;
+ // delete [] fYt; fYt = NULL;
+ // delete [] fYa; fYa = NULL;
+}
+
+//_____________________________________________________________________________
+inline
+void THcDC::ClearEvent()
+{
+ // Reset per-event data.
+ fNhits = 0;
+
+ for(Int_t i=0;i<fNChambers;i++) {
+ fChambers[i]->Clear();
+ }
+
+
+ // fTrackProj->Clear();
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::Decode( const THaEvData& evdata )
+{
+
+ ClearEvent();
+
+ // Get the Hall C style hitlist (fRawHitList) for this event
+ fNhits = THcHitList::DecodeToHitList(evdata);
+
+ // Let each plane get its hits
+ Int_t nexthit = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
+ }
+
+ // Let each chamber get its hits
+ for(Int_t ic=0;ic<fNChambers;ic++) {
+ fChambers[ic]->ProcessHits();
+ }
+#if 0
+ // fRawHitList is TClones array of THcRawDCHit objects
+ for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+ THcRawDCHit* hit = (THcRawDCHit *) fRawHitList->At(ihit);
+ // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+ // << endl;
+ for(Int_t imhit = 0; imhit < hit->fNHits; imhit++) {
+ // cout << " " << imhit << " " << hit->fTDC[imhit]
+ // << endl;
+ }
+ }
+ // cout << endl;
+#endif
+
+ return fNhits;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::ApplyCorrections( void )
+{
+ return(0);
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ )
+{
+ // Calculation of coordinates of particle track cross point with scint
+ // 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.);
+ for(Int_t i=0;i<fNChambers;i++) {
+
+ fChambers[i]->FindSpacePoints();
+ fChambers[i]->CorrectHitTimes();
+ }
+
+ ApplyCorrections();
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::FineTrack( TClonesArray& tracks )
+{
+ // Reconstruct coordinates of particle track cross point with scintillator
+ // 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 scint
+ // plane in the detector coordinate system. For this, parameters of track
+ // reconstructed in THaVDC::FineTrack() are used.
+
+ return 0;
+}
+
+ClassImp(THcDC)
+////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcDC.h b/src/THcDC.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a69711388ef1c5685a7dd47174e4060e39868ce
--- /dev/null
+++ b/src/THcDC.h
@@ -0,0 +1,150 @@
+#ifndef ROOT_THcDC
+#define ROOT_THcDC
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcDC //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THaTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcRawDCHit.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
+#include "TMath.h"
+
+//class THaScCalib;
+class TClonesArray;
+
+class THcDC : public THaTrackingDetector, public THcHitList {
+
+public:
+ THcDC( const char* name, const char* description = "",
+ THaApparatus* a = NULL );
+ virtual ~THcDC();
+
+ virtual Int_t Decode( const THaEvData& );
+ virtual EStatus Init( const TDatime& run_time );
+ virtual Int_t CoarseTrack( TClonesArray& tracks );
+ virtual Int_t FineTrack( 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; }
+
+ Int_t GetNWires(Int_t plane) const { return fNWires[plane-1];}
+ Int_t GetNChamber(Int_t plane) const { return fNChamber[plane-1];}
+ Int_t GetWireOrder(Int_t plane) const { return fWireOrder[plane-1];}
+ Int_t GetPitch(Int_t plane) const { return fPitch[plane-1];}
+ Int_t GetCentralWire(Int_t plane) const { return fCentralWire[plane-1];}
+ Int_t GetTdcWinMin(Int_t plane) const { return fTdcWinMin[plane-1];}
+ Int_t GetTdcWinMax(Int_t plane) const { return fTdcWinMax[plane-1];}
+
+ Double_t GetZPos(Int_t plane) const { return fZPos[plane-1];}
+ Double_t GetAlphaAngle(Int_t plane) const { return fAlphaAngle[plane-1];}
+ Double_t GetBetaAngle(Int_t plane) const { return fBetaAngle[plane-1];}
+ Double_t GetGammaAngle(Int_t plane) const { return fGammaAngle[plane-1];}
+
+ Int_t GetMinHits(Int_t chamber) const { return fMinHits[chamber-1];}
+ Int_t GetMaxHits(Int_t chamber) const { return fMaxHits[chamber-1];}
+ Int_t GetMinCombos(Int_t chamber) const { return fMinCombos[chamber-1];}
+ Double_t GetSpacePointCriterion(Int_t chamber) const { return TMath::Sqrt(fSpace_Point_Criterion2[chamber-1]);}
+ Double_t GetCentralTime(Int_t plane) const { return fCentralTime[plane-1];}
+ Int_t GetDriftTimeSign(Int_t plane) const { return fDriftTimeSign[plane-1];}
+
+ Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
+ Double_t GetSigma(Int_t plane) const { return fSigma[plane-1];}
+
+ Double_t GetNSperChan() const { return fNSperChan;}
+
+ Double_t GetCenter(Int_t plane) const {
+ Int_t chamber = GetNChamber(plane)-1;
+ return
+ fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
+ fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
+ }
+ // friend class THaScCalib;
+
+ THcDC(); // for ROOT I/O
+protected:
+
+ // Calibration
+
+ // Per-event data
+ Int_t fNhits;
+
+ // Potential Hall C parameters. Mostly here for demonstration
+ Int_t fNPlanes;
+ char** fPlaneNames;
+ Int_t fNChambers;
+
+ Double_t fNSperChan; /* TDC bin size */
+ Double_t fWireVelocity;
+
+ // Each of these will be dimensioned with the number of chambers
+ Double_t* fXCenter;
+ Double_t* fYCenter;
+ Int_t* fMinHits;
+ Int_t* fMaxHits;
+ Int_t* fMinCombos;
+ Double_t* fSpace_Point_Criterion2;
+
+ // Each of these will be dimensioned with the number of planes
+ // A THcDCPlane class object will need to access the value for
+ // its plane number. Should we have a Get method for each or
+ Int_t* fTdcWinMin;
+ Int_t* fTdcWinMax;
+ Int_t* fCentralTime;
+ Int_t* fNWires; // Number of wires per plane
+ Int_t* fNChamber;
+ Int_t* fWireOrder;
+ Int_t* fDriftTimeSign;
+
+ Double_t* fZPos;
+ Double_t* fAlphaAngle;
+ Double_t* fBetaAngle;
+ Double_t* fGammaAngle;
+ Double_t* fPitch;
+ Double_t* fCentralWire;
+ Double_t* fPlaneTimeZero;
+ Double_t* fSigma;
+
+ THcDriftChamberPlane** fPlanes; // List of plane objects
+ THcDriftChamber** fChambers; // List of chamber 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 );
+
+ void Setup(const char* name, const char* description);
+
+ ClassDef(THcDC,0) // Drift Chamber class
+};
+
+////////////////////////////////////////////////////////////////////////////////
+
+#endif
diff --git a/src/THcDCHit.cxx b/src/THcDCHit.cxx
index 0484b1ad39d1f2d1ee0171dfe0d2f5194f6a0ce0..3b1523c444ce0d4088ab5b422359d7e3d1cfd9ac 100644
--- a/src/THcDCHit.cxx
+++ b/src/THcDCHit.cxx
@@ -9,10 +9,31 @@
#include "THcDCHit.h"
#include "THcDCTimeToDistConv.h"
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
ClassImp(THcDCHit)
const Double_t THcDCHit::kBig = 1.e38; // Arbitrary large value
+//_____________________________________________________________________________
+void THcDCHit::Print( Option_t* opt ) const
+{
+ // Print hit info
+
+ cout << "Hit: wire=" << GetWireNum()
+ << "/" << (fWirePlane ? fWirePlane->GetName() : "??")
+ << " wpos=" << GetPos()
+ << " time=" << GetTime()
+ << " drift=" << GetDist();
+ // << " res=" << GetResolution()
+ // << " z=" << GetZ()
+ if( *opt != 'C' )
+ cout << endl;
+}
+
//_____________________________________________________________________________
Double_t THcDCHit::ConvertTimeToDist()
{
@@ -51,15 +72,17 @@ Int_t THcDCHit::Compare( const TObject* obj ) const
Int_t myWireNum = fWire->GetNum();
Int_t hitWireNum = hit->GetWire()->GetNum();
- // Compare wire numbers
+ Int_t myPlaneNum = GetPlaneNum();
+ Int_t hitPlaneNum = hit->GetPlaneNum();
+ if (myPlaneNum < hitPlaneNum) return -1;
+ if (myPlaneNum > hitPlaneNum) return 1;
+ // If planes are the same, compare wire numbers
if (myWireNum < hitWireNum) return -1;
if (myWireNum > hitWireNum) return 1;
- if (myWireNum == hitWireNum) {
- // If wire numbers are the same, compare times
- Double_t hitTime = hit->GetTime();
- if (fTime < hitTime) return -1;
- if (fTime > hitTime) return 1;
- }
+ // If wire numbers are the same, compare times
+ Double_t hitTime = hit->GetTime();
+ if (fTime < hitTime) return -1;
+ if (fTime > hitTime) return 1;
return 0;
}
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index 79f88e7123f7a8d39d805719e7ce2b74ec831316..e32cc6159bd9590170df39d4ebb15a175aebcec3 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -9,20 +9,26 @@
#include "TObject.h"
#include "THcDCWire.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
#include <cstdio>
class THcDCHit : public TObject {
public:
- THcDCHit( THcDCWire* wire=NULL, Int_t rawtime=0, Double_t time=0.0 ) :
- fWire(wire), fRawTime(rawtime), fTime(time), fDist(0.0), ftrDist(kBig) {
+ THcDCHit( THcDCWire* wire=NULL, Int_t rawtime=0, Double_t time=0.0,
+ THcDriftChamberPlane* wp=0) :
+ fWire(wire), fRawTime(rawtime), fTime(time), fWirePlane(wp),
+ fDist(0.0), ftrDist(kBig) {
ConvertTimeToDist();
+ fCorrected = 0;
}
virtual ~THcDCHit() {}
virtual Double_t ConvertTimeToDist();
Int_t Compare ( const TObject* obj ) const;
Bool_t IsSortable () const { return kTRUE; }
+ virtual void Print( Option_t* opt="" ) const;
// Get and Set Functions
THcDCWire* GetWire() const { return fWire; }
@@ -31,24 +37,40 @@ public:
Double_t GetTime() const { return fTime; }
Double_t GetDist() const { return fDist; }
Double_t GetPos() const { return fWire->GetPos(); } //Position of hit wire
+ Double_t GetCoord() const { return fCoord; }
Double_t GetdDist() const { return fdDist; }
+ Int_t GetCorrectedStatus() const { return fCorrected; }
+
+ THcDriftChamberPlane* GetWirePlane() const { return fWirePlane; }
+
void SetWire(THcDCWire * wire) { fWire = wire; }
void SetRawTime(Int_t time) { fRawTime = time; }
void SetTime(Double_t time) { fTime = time; }
void SetDist(Double_t dist) { fDist = dist; }
+ void SetLeftRight(Int_t lr) { fCoord = GetPos() + lr*fDist;}
void SetdDist(Double_t ddist) { fdDist = ddist; }
void SetFitDist(Double_t dist) { ftrDist = dist; }
+ Int_t GetPlaneNum() const { return fWirePlane->GetPlaneNum(); }
+ Int_t GetPlaneIndex() const { return fWirePlane->GetPlaneIndex(); }
+ Int_t GetChamberNum() const { return fWirePlane->GetChamberNum(); }
+ void SetCorrectedStatus(Int_t c) { fCorrected = c; }
protected:
static const Double_t kBig; //!
- THcDCWire* fWire; // Wire on which the hit occurred
+ THcDCWire* fWire; // Wire on which the hit occurred
Int_t fRawTime; // TDC value (channels)
Double_t fTime; // Time corrected for time offset of wire (s)
+ THcDriftChamberPlane* fWirePlane; //! Pointer to parent wire plane
Double_t fDist; // (Perpendicular) Drift Distance
+ Double_t fCoord; // Actual coordinate of hit
Double_t fdDist; // uncertainty in fDist (for chi2 calc)
Double_t ftrDist; // (Perpendicular) distance from the track
+ Int_t fCorrected; // Has this hit been corrected?
+
+ THcDriftChamber* fChamber; //! Pointer to parent wire plane
+
private:
THcDCHit( const THcDCHit& );
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 1d53b20c58f52e7465596ea9240a58bd1755eb6b..59312b75dd42c826c949ed2d2169f33224831529 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -1,18 +1,17 @@
///////////////////////////////////////////////////////////////////////////////
// //
-// THcDriftChamber //
+// THcDriftChamber //
// //
-// Class for a generic hodoscope consisting of multiple //
-// planes with multiple paddles with phototubes on both ends. //
-// This differs from Hall A scintillator class in that it is the whole //
-// hodoscope array, not just one plane. //
+// Subdetector class to hold a bunch of planes constituting a chamber //
+// This class will be created by the THcDC class which will also create //
+// the plane objects. //
+// The THcDC class will then pass this class a list of the planes. //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THcDriftChamber.h"
-#include "THaEvData.h"
-#include "THaDetMap.h"
-#include "THcDetectorMap.h"
+#include "THcDC.h"
+#include "THcDCHit.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "VarDef.h"
@@ -20,6 +19,7 @@
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"
+#include "TVectorD.h"
#include "THaTrackProj.h"
@@ -27,20 +27,23 @@
#include <cstdio>
#include <cstdlib>
#include <iostream>
+#include <iomanip>
using namespace std;
//_____________________________________________________________________________
THcDriftChamber::THcDriftChamber(
const char* name, const char* description,
- THaApparatus* apparatus ) :
- THaTrackingDetector(name,description,apparatus)
+ const Int_t chambernum, THaDetectorBase* parent ) :
+ THaSubDetector(name,description,parent)
{
// Constructor
// fTrackProj = new TClonesArray( "THaTrackProj", 5 );
fNPlanes = 0; // No planes until we make them
+ fChamberNum = chambernum;
+
}
//_____________________________________________________________________________
@@ -48,7 +51,6 @@ void THcDriftChamber::Setup(const char* name, const char* description)
{
char prefix[2];
- char parname[100];
THaApparatus *app = GetApparatus();
if(app) {
@@ -60,181 +62,153 @@ void THcDriftChamber::Setup(const char* name, const char* description)
prefix[0]=tolower(app->GetName()[0]);
prefix[1]='\0';
- string planenamelist;
- DBRequest list[]={
- {"dc_num_planes",&fNPlanes, kInt},
- {"dc_num_chambers",&fNChambers, kInt},
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
- {"dc_plane_names",&planenamelist, kString},
- {0}
- };
-
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- cout << planenamelist << endl;
- cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
-
- vector<string> plane_names = vsplit(planenamelist);
-
- if(plane_names.size() != (UInt_t) fNPlanes) {
- cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
- // Should quit. Is there an official way to quit?
- }
- fPlaneNames = new char* [fNPlanes];
- for(Int_t i=0;i<fNPlanes;i++) {
- fPlaneNames[i] = new char[plane_names[i].length()];
- strcpy(fPlaneNames[i], plane_names[i].c_str());
- }
-
- char *desc = new char[strlen(description)+100];
- fPlanes = new THcDriftChamberPlane* [fNPlanes];
-
- for(Int_t i=0;i<fNPlanes;i++) {
- strcpy(desc, description);
- strcat(desc, " Plane ");
- strcat(desc, fPlaneNames[i]);
-
- fPlanes[i] = new THcDriftChamberPlane(fPlaneNames[i], desc, i+1, this);
- cout << "Created Drift Chamber Plane " << fPlaneNames[i] << ", " << desc << endl;
-
- }
-
}
//_____________________________________________________________________________
THcDriftChamber::THcDriftChamber( ) :
- THaTrackingDetector()
+ THaSubDetector()
{
// Constructor
}
+//_____________________________________________________________________________
+Int_t THcDriftChamber::Decode( const THaEvData& evdata )
+{
+ cout << "THcDriftChamber::Decode called" << endl;
+ return 0;
+}
+
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcDriftChamber::Init( const TDatime& date )
{
- static const char* const here = "Init()";
+ // static const char* const here = "Init()";
- Setup(GetName(), GetTitle()); // Create the subdetectors here
+ Setup(GetName(), GetTitle());
- // Should probably put this in ReadDatabase as we will know the
- // maximum number of hits after setting up the detector map
- THcHitList::InitHitList(fDetMap, "THcRawDCHit", 1000);
-
EStatus status;
// This triggers call of ReadDatabase and DefineVariables
- if( (status = THaTrackingDetector::Init( date )) )
+ if( (status = THaSubDetector::Init( date )) )
return fStatus=status;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- if((status = fPlanes[ip]->Init( date ))) {
- return fStatus=status;
- }
- }
-
- // Replace with what we need for Hall C
- // const DataDest tmp[NDEST] = {
- // { &fRTNhit, &fRANhit, fRT, fRT_c, fRA, fRA_p, fRA_c, fROff, fRPed, fRGain },
- // { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
- // };
- // memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
-
- // Will need to determine which apparatus it belongs to and use the
- // appropriate detector ID in the FillMap call
- char EngineDID[4];
+ return fStatus = kOK;
+}
- EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
- EngineDID[1] = 'D';
- EngineDID[2] = 'C';
- EngineDID[3] = '\0';
-
- if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
- Error( Here(here), "Error filling detectormap for %s.",
- EngineDID);
- return kInitError;
+void THcDriftChamber::AddPlane(THcDriftChamberPlane *plane)
+{
+ cout << "Added plane " << plane->GetPlaneNum() << " to chamber " << fChamberNum << endl;
+ plane->SetPlaneIndex(fNPlanes);
+ fPlanes[fNPlanes] = plane;
+
+ // HMS Specific
+ // Check if this is a Y Plane
+ if(plane->GetPlaneNum() == YPlaneNum) {
+ YPlaneInd = fNPlanes;
+ } else if (plane->GetPlaneNum() == YPlanePNum) {
+ YPlanePInd = fNPlanes;
}
+ fNPlanes++;
- return fStatus = kOK;
+ return;
}
//_____________________________________________________________________________
Int_t THcDriftChamber::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()";
+ cout << "THcDriftChamber::ReadDatabase()" << endl;
char prefix[2];
- char parname[100];
-
- // Read data from database
- // Pull values from the THcParmList instead of reading a database
- // file like Hall A does.
-
- // 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)
-
prefix[0]=tolower(GetApparatus()->GetName()[0]);
-
prefix[1]='\0';
-
- fXCenter = new Double_t [fNChambers];
- fYCenter = new Double_t [fNChambers];
-
- fTdcWinMin = new Int_t [fNPlanes];
- fTdcWinMax = new Int_t [fNPlanes];
- fCentralTime = new Int_t [fNPlanes];
- fNWires = new Int_t [fNPlanes];
- fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
- fWireOrder = new Int_t [fNPlanes]; // Wire readout order
- fDriftTimeSign = new Int_t [fNPlanes];
-
- fZPos = new Double_t [fNPlanes];
- fAlphaAngle = new Double_t [fNPlanes];
- fBetaAngle = new Double_t [fNPlanes];
- fGammaAngle = new Double_t [fNPlanes];
- fPitch = new Double_t [fNPlanes];
- fCentralWire = new Double_t [fNPlanes];
- fPlaneTimeZero = new Double_t [fNPlanes];
-
-
DBRequest list[]={
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
-
- {"dc_xcenter", fXCenter, kDouble, fNChambers},
- {"dc_ycenter", fYCenter, kDouble, fNChambers},
-
- {"dc_tdc_min_win", fTdcWinMin, kInt, fNPlanes},
- {"dc_tdc_max_win", fTdcWinMax, kInt, fNPlanes},
- {"dc_central_time", fCentralTime, kInt, fNPlanes},
- {"dc_nrwire", fNWires, kInt, fNPlanes},
- {"dc_chamber_planes", fNChamber, kInt, fNPlanes},
- {"dc_wire_counting", fWireOrder, kInt, fNPlanes},
- {"dc_drifttime_sign", fDriftTimeSign, kInt, fNPlanes},
-
- {"dc_zpos", fZPos, kDouble, fNPlanes},
- {"dc_alpha_angle", fAlphaAngle, kDouble, fNPlanes},
- {"dc_beta_angle", fBetaAngle, kDouble, fNPlanes},
- {"dc_gamma_angle", fGammaAngle, kDouble, fNPlanes},
- {"dc_pitch", fPitch, kDouble, fNPlanes},
- {"dc_central_wire", fCentralWire, kDouble, fNPlanes},
- {"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
+ {"_remove_sppt_if_one_y_plane",&fRemove_Sppt_If_One_YPlane, kInt},
+ {"dc_wire_velocity", &fWireVelocity, kDouble},
+ {"SmallAngleApprox", &fSmallAngleApprox, kInt},
+ {"stub_max_xpdiff", &fStubMaxXPDiff, kDouble},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- cout << "Plane counts:";
- for(Int_t i=0;i<fNPlanes;i++) {
- cout << " " << fNWires[i];
+ // Get parameters parent knows about
+ THcDC* fParent;
+ fParent = (THcDC*) GetParent();
+ fMinHits = fParent->GetMinHits(fChamberNum);
+ fMaxHits = fParent->GetMaxHits(fChamberNum);
+ fMinCombos = fParent->GetMinCombos(fChamberNum);
+
+ cout << "Chamber " << fChamberNum << " Min/Max: " << fMinHits << "/" << fMaxHits << endl;
+
+ fSpacePointCriterion = fParent->GetSpacePointCriterion(fChamberNum);
+ fSpacePointCriterion2 = fSpacePointCriterion*fSpacePointCriterion;
+
+ // HMS Specific
+ // Hard code Y plane numbers. Should be able to get from wire angle
+ if(fChamberNum == 1) {
+ YPlaneNum = 2;
+ YPlanePNum = 5;
+ } else {
+ YPlaneNum = 8;
+ YPlanePNum = 11;
+ }
+
+ // Generate the HAA3INV matrix for all the acceptable combinations
+ // of hit planes. Try to make it as generic as possible
+ // pindex=0 -> Plane 1 missing, pindex5 -> plane 6 missing. Won't
+ // replicate the exact values used in the ENGINE, because the engine
+ // had one big list of matrices for both chambers, while here we will
+ // have a list just for one chamber. Also, call pindex, pmindex as
+ // we tend to use pindex as a plane index.
+ fCosBeta = new Double_t [fNPlanes];
+ fSinBeta = new Double_t [fNPlanes];
+ fTanBeta = new Double_t [fNPlanes];
+ fSigma = new Double_t [fNPlanes];
+ fPsi0 = new Double_t [fNPlanes];
+ fStubCoefs = new Double_t* [fNPlanes];
+ Int_t allplanes=0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fCosBeta[ip] = TMath::Cos(fPlanes[ip]->GetBeta());
+ fSinBeta[ip] = TMath::Sin(fPlanes[ip]->GetBeta());
+ fTanBeta[ip] = fSinBeta[ip]/fCosBeta[ip];
+ fSigma[ip] = fPlanes[ip]->GetSigma();
+ fPsi0[ip] = fPlanes[ip]->GetPsi0();
+ fStubCoefs[ip] = fPlanes[ip]->GetStubCoef();
+ allplanes |= 1<<ip;
+ }
+ // Unordered map introduced in C++-11
+ // Can use unordered_map if using C++-11
+ // May not want to use map a all for performance, but using it now
+ // for code clarity
+ for(Int_t ipm1=0;ipm1<fNPlanes+1;ipm1++) { // Loop over missing plane1
+ for(Int_t ipm2=ipm1;ipm2<fNPlanes+1;ipm2++) {
+ if(ipm1==ipm2 && ipm1<fNPlanes) continue;
+ TMatrixDSym* AA3 = new TMatrixDSym(3);
+ for(Int_t i=0;i<3;i++) {
+ for(Int_t j=i;j<3;j++) {
+ (*AA3)[j][i] = 0.0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(ipm1 != ip && ipm2 != ip) {
+ (*AA3)[j][i] += fStubCoefs[ip][i]*fStubCoefs[ip][j];
+ }
+ }
+ }
+ }
+ // Use the bit pattern of missing planes as an hash index.
+ // Perhaps it is more efficient to just use a regular array instead of map
+ Int_t bitpat = allplanes & ~(1<<ipm1) & ~(1<<ipm2);
+ // Should check that it is invertable
+ AA3->Invert();
+ fAA3Inv[bitpat] = AA3;
+ }
+ }
+
+#if 0
+ for(map<int,TMatrixDSym*>::iterator pm=fHaa3map.begin();
+ pm != fHaa3map.end(); pm++) {
+ cout << setbase(8) << pm->first << endl;
+ fAA3Inv[pm->first]->Print();
}
- cout << endl;
+#endif
fIsInit = true;
-
return kOK;
}
@@ -249,165 +223,863 @@ Int_t THcDriftChamber::DefineVariables( EMode mode )
// Register variables in global list
// RVarDef vars[] = {
- // { "nlthit", "Number of Left paddles TDC times", "fLTNhit" },
- // { "nrthit", "Number of Right paddles TDC times", "fRTNhit" },
- // { "nlahit", "Number of Left paddles ADCs amps", "fLANhit" },
- // { "nrahit", "Number of Right paddles ADCs amps", "fRANhit" },
- // { "lt", "TDC values left side", "fLT" },
- // { "lt_c", "Corrected times left side", "fLT_c" },
- // { "rt", "TDC values right side", "fRT" },
- // { "rt_c", "Corrected times right side", "fRT_c" },
- // { "la", "ADC values left side", "fLA" },
- // { "la_p", "Corrected ADC values left side", "fLA_p" },
- // { "la_c", "Corrected ADC values left side", "fLA_c" },
- // { "ra", "ADC values right side", "fRA" },
- // { "ra_p", "Corrected ADC values right side", "fRA_p" },
- // { "ra_c", "Corrected ADC values right side", "fRA_c" },
- // { "nthit", "Number of paddles with l&r TDCs", "fNhit" },
- // { "t_pads", "Paddles with l&r coincidence TDCs", "fHitPad" },
- // { "y_t", "y-position from timing (m)", "fYt" },
- // { "y_adc", "y-position from amplitudes (m)", "fYa" },
- // { "time", "Time of hit at plane (s)", "fTime" },
- // { "dtime", "Est. uncertainty of time (s)", "fdTime" },
- // { "dedx", "dEdX-like deposited in paddle", "fAmpl" },
- // { "troff", "Trigger offset for paddles", "fTrigOff"},
- // { "trn", "Number of tracks for hits", "GetNTracks()" },
- // { "trx", "x-position of track in det plane", "fTrackProj.THaTrackProj.fX" },
- // { "try", "y-position of track in det plane", "fTrackProj.THaTrackProj.fY" },
- // { "trpath", "TRCS pathlen of track to det plane","fTrackProj.THaTrackProj.fPathl" },
- // { "trdx", "track deviation in x-position (m)", "fTrackProj.THaTrackProj.fdX" },
- // { "trpad", "paddle-hit associated with track", "fTrackProj.THaTrackProj.fChannel" },
+ // { "nhit", "Number of DC hits", "fNhits" },
// { 0 }
// };
// return DefineVarsFromList( vars, mode );
return kOK;
-}
-//_____________________________________________________________________________
-THcDriftChamber::~THcDriftChamber()
+}
+void THcDriftChamber::ProcessHits( void)
{
- // Destructor. Remove variables from global list.
+ // Make a list of hits for whole chamber
+ fNhits = 0;
- if( fIsSetup )
- RemoveVariables();
- if( fIsInit )
- DeleteArrays();
- if (fTrackProj) {
- fTrackProj->Clear();
- delete fTrackProj; fTrackProj = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ TClonesArray* hitsarray = fPlanes[ip]->GetHits();
+ for(Int_t ihit=0;ihit<fPlanes[ip]->GetNHits();ihit++) {
+ fHits[fNhits++] = static_cast<THcDCHit*>(hitsarray->At(ihit));
+ }
}
+ // cout << "ThcDriftChamber::ProcessHits() " << fNhits << " hits" << endl;
}
-//_____________________________________________________________________________
-void THcDriftChamber::DeleteArrays()
+
+Int_t THcDriftChamber::FindSpacePoints( void )
{
- // Delete member arrays. Used by destructor.
+ // Following h_pattern_recognition.f, but just for one chamber
- delete [] fNWires; fNWires = NULL;
- // delete [] fSpacing; fSpacing = NULL;
- // delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete?
-
- // delete [] fRA_c; fRA_c = NULL;
- // delete [] fRA_p; fRA_p = NULL;
- // delete [] fRA; fRA = NULL;
- // delete [] fLA_c; fLA_c = NULL;
- // delete [] fLA_p; fLA_p = NULL;
- // delete [] fLA; fLA = NULL;
- // delete [] fRT_c; fRT_c = NULL;
- // delete [] fRT; fRT = NULL;
- // delete [] fLT_c; fLT_c = NULL;
- // delete [] fLT; fLT = NULL;
-
- // delete [] fRGain; fRGain = NULL;
- // delete [] fLGain; fLGain = NULL;
- // delete [] fRPed; fRPed = NULL;
- // delete [] fLPed; fLPed = NULL;
- // delete [] fROff; fROff = NULL;
- // delete [] fLOff; fLOff = NULL;
- // delete [] fTWalkPar; fTWalkPar = NULL;
- // delete [] fTrigOff; fTrigOff = NULL;
-
- // delete [] fHitPad; fHitPad = NULL;
- // delete [] fTime; fTime = NULL;
- // delete [] fdTime; fdTime = NULL;
- // delete [] fYt; fYt = NULL;
- // delete [] fYa; fYa = NULL;
+ // Code below is specifically for HMS chambers with Y and Y' planes
+
+ Int_t yplane_hitind=0;
+ Int_t yplanep_hitind=0;
+
+ fNSpacePoints=0;
+ fEasySpacePoint = 0;
+ // Should really build up array of hits
+ if(fNhits >= fMinHits && fNhits < fMaxHits) {
+ /* Has this list of hits already been cut the way it should have at this point? */
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ Int_t ip=thishit->GetPlaneNum(); // This is the absolute plane mumber
+ if(ip==YPlaneNum) yplane_hitind = ihit;
+ if(ip==YPlanePNum) yplanep_hitind = ihit;
+ }
+ // fPlanes is the array of planes for this chamber.
+ // cout << fPlanes[YPlaneInd]->GetNHits() << " "
+ // << fPlanes[YPlanePInd]->GetNHits() << " " << endl;
+ if (fPlanes[YPlaneInd]->GetNHits() == 1 && fPlanes[YPlanePInd]->GetNHits() == 1) {
+ cout << fHits[yplane_hitind]->GetWireNum() << " "
+ << fHits[yplane_hitind]->GetPos() << " "
+ << fHits[yplanep_hitind]->GetWireNum() << " "
+ << fHits[yplanep_hitind]->GetPos() << " "
+ << fSpacePointCriterion << endl;
+ }
+ if(fPlanes[YPlaneInd]->GetNHits() == 1 && fPlanes[YPlanePInd]->GetNHits() == 1
+ && TMath::Abs(fHits[yplane_hitind]->GetPos() - fHits[yplanep_hitind]->GetPos())
+ < fSpacePointCriterion
+ && fNhits <= 6) { // An easy case, probably one hit per plane
+ fEasySpacePoint = FindEasySpacePoint(yplane_hitind, yplanep_hitind);
+ }
+ if(!fEasySpacePoint) { // It's not easy
+ FindHardSpacePoints();
+ }
+
+ // We have our space points for this chamber
+ cout << fNSpacePoints << " Space Points found" << endl;
+ if(fNSpacePoints > 0) {
+ if(fRemove_Sppt_If_One_YPlane == 1) {
+ // The routine is specific to HMS
+ Int_t ndest=DestroyPoorSpacePoints();
+ cout << ndest << " Poor space points destroyed" << endl;
+ // Loop over space points and remove those with less than 4 planes
+ // hit and missing hits in Y,Y' planes
+ }
+ if(fNSpacePoints == 0) cout << "DestroyPoorSpacePoints() killed SP" << endl;
+ Int_t nadded=SpacePointMultiWire();
+ if (nadded) cout << nadded << " Space Points added with SpacePointMultiWire()" << endl;
+ ChooseSingleHit();
+ SelectSpacePoints();
+ if(fNSpacePoints == 0) cout << "SelectSpacePoints() killed SP" << endl;
+ }
+ cout << fNSpacePoints << " Space Points remain" << endl;
+ // Add these space points to the total list of space points for the
+ // the DC package. Do this in THcDC.cxx.
+#if 0
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t np = fPlanes[ip]->GetPlaneNum(); // Actuall plane number of this plane
+ // (0-11) or (1-12)?
+ TClonesArray* fHits = fPlanes[ip]->GetHits();
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) { // Looping over all hits in all planes of the chamber
+ THcDCHit* hit = static_cast<THcDCHit*>(fHits->At(ihit));
+ //
+
+ }
+ }
+#endif
+ }
+ return(fNSpacePoints);
}
//_____________________________________________________________________________
-inline
-void THcDriftChamber::ClearEvent()
+// HMS Specific
+Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hitind)
{
- // Reset per-event data.
+ // Simplified HMS find_space_point routing. It is given all y hits and
+ // checks to see if all x-like hits are close enough together to make
+ // a space point.
+
+ Int_t easy_space_point=0;
+ cout << "Doing Find Easy Space Point" << endl;
+ Double_t yt = (fHits[yplane_hitind]->GetPos() + fHits[yplanep_hitind]->GetPos())/2.0;
+ Double_t xt = 0.0;
+ Int_t num_xhits = 0;
+ Double_t x_pos[MAX_HITS_PER_POINT];
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
+ // ysp and xsp are from h_generate_geometry
+ x_pos[ihit] = (thishit->GetPos()
+ -yt*thishit->GetWirePlane()->GetYsp())
+ /thishit->GetWirePlane()->GetXsp();
+ xt += x_pos[ihit];
+ num_xhits++;
+ } else {
+ x_pos[ihit] = 0.0;
+ }
+ }
+ xt = (num_xhits>0?xt/num_xhits:0.0);
+ cout << "xt=" << xt << endl;
+ easy_space_point = 1; // Assume we have an easy space point
+ // Rule it out if x points don't cluster well enough
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ cout << "Hit " << ihit << " " << x_pos[ihit] << " " << xt << endl;
+ if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
+ if(TMath::Abs(xt-x_pos[ihit]) >= fSpacePointCriterion)
+ { easy_space_point=0; break;}
+ }
+ }
+ if(easy_space_point) { // Register the space point
+ cout << "Easy Space Point " << xt << " " << yt << endl;
+ fNSpacePoints = 1;
+ fSpacePoints[0].x = xt;
+ fSpacePoints[0].y = yt;
+ fSpacePoints[0].nhits = fNhits;
+ fSpacePoints[0].ncombos = 0; // No combos
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ fSpacePoints[0].hits[ihit] = thishit;
+ }
+ }
+ return(easy_space_point);
+}
- fTrackProj->Clear();
+//_____________________________________________________________________________
+// Generic
+Int_t THcDriftChamber::FindHardSpacePoints()
+{
+ Int_t MAX_NUMBER_PAIRS=1000; // Where does this get set?
+ struct Pair {
+ THcDCHit* hit1;
+ THcDCHit* hit2;
+ Double_t x, y;
+ };
+ Pair pairs[MAX_NUMBER_PAIRS];
+ //
+ Int_t ntest_points=0;
+ for(Int_t ihit1=0;ihit1<fNhits-1;ihit1++) {
+ THcDCHit* hit1=fHits[ihit1];
+ THcDriftChamberPlane* plane1 = hit1->GetWirePlane();
+ for(Int_t ihit2=ihit1+1;ihit2<fNhits;ihit2++) {
+ if(ntest_points < MAX_NUMBER_PAIRS) {
+ THcDCHit* hit2=fHits[ihit2];
+ THcDriftChamberPlane* plane2 = hit2->GetWirePlane();
+ Double_t determinate = plane1->GetXsp()*plane2->GetYsp()
+ -plane1->GetYsp()*plane2->GetXsp();
+ if(TMath::Abs(determinate) > 0.3) { // 0.3 is sin(alpha1-alpha2)=sin(17.5)
+ pairs[ntest_points].hit1 = hit1;
+ pairs[ntest_points].hit2 = hit2;
+ pairs[ntest_points].x = (hit1->GetPos()*plane2->GetYsp()
+ - hit2->GetPos()*plane1->GetYsp())
+ /determinate;
+ pairs[ntest_points].y = (hit2->GetPos()*plane1->GetXsp()
+ - hit1->GetPos()*plane2->GetXsp())
+ /determinate;
+ ntest_points++;
+ }
+ }
+ }
+ }
+ Int_t ncombos=0;
+ struct Combo {
+ Pair* pair1;
+ Pair* pair2;
+ };
+ Combo combos[10*MAX_NUMBER_PAIRS];
+ for(Int_t ipair1=0;ipair1<ntest_points-1;ipair1++) {
+ for(Int_t ipair2=ipair1+1;ipair2<ntest_points;ipair2++) {
+ if(ncombos < 10*MAX_NUMBER_PAIRS) {
+ Double_t dist2 = pow(pairs[ipair1].x - pairs[ipair2].x,2)
+ + pow(pairs[ipair1].y - pairs[ipair2].y,2);
+ if(dist2 <= fSpacePointCriterion2) {
+ combos[ncombos].pair1 = &pairs[ipair1];
+ combos[ncombos].pair2 = &pairs[ipair2];
+ ncombos++;
+ }
+ }
+ }
+ }
+ // Loop over all valid combinations and build space points
+ for(Int_t icombo=0;icombo<ncombos;icombo++) {
+ THcDCHit* hits[4];
+ hits[0]=combos[icombo].pair1->hit1;
+ hits[1]=combos[icombo].pair1->hit2;
+ hits[2]=combos[icombo].pair2->hit1;
+ hits[3]=combos[icombo].pair2->hit2;
+ // Get Average Space point xt, yt
+ Double_t xt = (combos[icombo].pair1->x + combos[icombo].pair2->x)/2.0;
+ Double_t yt = (combos[icombo].pair1->y + combos[icombo].pair2->y)/2.0;
+ // Loop over space points
+ if(fNSpacePoints > 0) {
+ Int_t add_flag=1;
+ for(Int_t ispace=0;ispace<fNSpacePoints;ispace++) {
+ if(fSpacePoints[ispace].nhits > 0) {
+ Double_t sqdist_test = pow(xt - fSpacePoints[ispace].x,2) +
+ pow(yt - fSpacePoints[ispace].y,2);
+ // I (who is I) want to be careful if sqdist_test is bvetween 1 and
+ // 3 fSpacePointCriterion2. Let me ignore not add a new point the
+ if(sqdist_test < 3*fSpacePointCriterion2) {
+ add_flag = 0; // do not add a new space point
+ }
+ if(sqdist_test < fSpacePointCriterion2) {
+ // This is a real match
+ // Add the new hits to the existing space point
+ Int_t iflag[4];
+ iflag[0]=0;iflag[1]=0;iflag[2]=0;iflag[3]=0;
+ // Find out which of the four hits in the combo are already
+ // in the space point under consideration so that we don't
+ // add duplicate hits to the space point
+ for(Int_t isp_hit=0;isp_hit<fSpacePoints[ispace].nhits;isp_hit++) {
+ for(Int_t icm_hit=0;icm_hit<4;icm_hit++) { // Loop over combo hits
+ if(fSpacePoints[ispace].hits[isp_hit]==hits[icm_hit]) {
+ iflag[icm_hit] = 1;
+ }
+ }
+ }
+ // Remove duplicated pionts in the combo so we don't add
+ // duplicate hits to the space point
+ for(Int_t icm1=0;icm1<3;icm1++) {
+ for(Int_t icm2=icm1+1;icm2<4;icm2++) {
+ if(hits[icm1]==hits[icm2]) {
+ iflag[icm2] = 1;
+ }
+ }
+ }
+ // Add the unique combo hits to the space point
+ for(Int_t icm=0;icm<4;icm++) {
+ if(iflag[icm]==0) {
+ fSpacePoints[ispace].hits[fSpacePoints[ispace].nhits++] = hits[icm];
+ }
+ }
+ fSpacePoints[ispace].ncombos++;
+ // Terminate loop since this combo can only belong to one space point
+ break;
+ }
+ }
+ }// End of loop over existing space points
+ // Create a new space point if more than 2*space_point_criteria
+ if(fNSpacePoints < MAX_SPACE_POINTS) {
+ if(add_flag) {
+ fSpacePoints[fNSpacePoints].nhits=2;
+ fSpacePoints[fNSpacePoints].ncombos=1;
+ fSpacePoints[fNSpacePoints].hits[0]=hits[0];
+ fSpacePoints[fNSpacePoints].hits[1]=hits[1];
+ fSpacePoints[fNSpacePoints].x = xt;
+ fSpacePoints[fNSpacePoints].y = yt;
+ if(hits[0] != hits[2] && hits[1] != hits[2]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[2];
+ }
+ if(hits[0] != hits[3] && hits[1] != hits[3]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[3];
+ }
+ fNSpacePoints++;
+ }
+ }
+ } else {// Create first space point
+ // This duplicates code above. Need to see if we can restructure
+ // to avoid
+ fSpacePoints[fNSpacePoints].nhits=2;
+ fSpacePoints[fNSpacePoints].ncombos=1;
+ fSpacePoints[fNSpacePoints].hits[0]=hits[0];
+ fSpacePoints[fNSpacePoints].hits[1]=hits[1];
+ fSpacePoints[fNSpacePoints].x = xt;
+ fSpacePoints[fNSpacePoints].y = yt;
+ if(hits[0] != hits[2] && hits[1] != hits[2]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[2];
+ }
+ if(hits[0] != hits[3] && hits[1] != hits[3]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[3];
+ }
+ fNSpacePoints++;
+ }//End check on 0 space points
+ }//End loop over combos
+ return(fNSpacePoints);
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::Decode( const THaEvData& evdata )
+// HMS Specific?
+Int_t THcDriftChamber::DestroyPoorSpacePoints()
{
+ Int_t nhitsperplane[fNPlanes];
- // Get the Hall C style hitlist (fRawHitList) for this event
- Int_t nhits = THcHitList::DecodeToHitList(evdata);
+ Int_t spacepointsgood[fNSpacePoints];
+ Int_t ngood=0;
- // Let each plane get its hits
- Int_t nexthit = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
+ for(Int_t i=0;i<fNSpacePoints;i++) {
+ spacepointsgood[i] = 0;
+ }
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t nplanes_hit = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nhitsperplane[ip] = 0;
+ }
+ // Count # hits in each plane for this space point
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) {
+ THcDCHit* hit=fSpacePoints[isp].hits[ihit];
+ // hit_order(hit) = ihit;
+ Int_t ip = hit->GetPlaneIndex();
+ nhitsperplane[ip]++;
+ }
+ // Count # planes that have hits
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(nhitsperplane[ip] > 0) {
+ nplanes_hit++;
+ }
+ }
+ if(nplanes_hit >= fMinHits && nhitsperplane[YPlaneInd]>0
+ && nhitsperplane[YPlanePInd] > 0) {
+ spacepointsgood[ngood++] = isp; // Build list of good points
+ } else {
+ // cout << "Missing Y-hit!!";
+ }
}
-#if 0
- // fRawHitList is TClones array of THcRawDCHit objects
- for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
- THcRawDCHit* hit = (THcRawDCHit *) fRawHitList->At(ihit);
- // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
- // << endl;
- for(Int_t imhit = 0; imhit < hit->fNHits; imhit++) {
- // cout << " " << imhit << " " << hit->fTDC[imhit]
- // << endl;
+ // Remove the bad space points
+ Int_t nremoved=fNSpacePoints-ngood;
+ fNSpacePoints = ngood;
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) { // New index num ber
+ Int_t osp=spacepointsgood[isp]; // Original index number
+ if(osp > isp) {
+ // Does this work, or do we have to copy each member?
+ // If it doesn't we should overload the = operator
+ fSpacePoints[isp] = fSpacePoints[osp];
}
}
- // cout << endl;
+ return nremoved;
+}
+
+//_____________________________________________________________________________
+// HMS Specific?
+ /*
+ Purpose and Methods : This routine loops over space points and
+ looks at all hits in the space
+ point. If more than 1 hit is in the same
+ plane then the space point is cloned with
+ all combinations of 1 wire per plane. The
+ requirements for cloning are: 1) at least
+ 4 planes fire, and 2) no more than 6 planes
+ have multiple hits.
+ */
+Int_t THcDriftChamber::SpacePointMultiWire()
+{
+ Int_t nhitsperplane[fNPlanes];
+ THcDCHit* hits_plane[fNPlanes][MAX_HITS_PER_POINT];
+
+ Int_t nsp_check;
+ Int_t nplanes_single;
+
+ Int_t nsp_tot=fNSpacePoints;
+
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t nplanes_hit = 0; // Number of planes with hits
+ Int_t nplanes_mult = 0; // Number of planes with multiple hits
+ Int_t nsp_new = 1;
+ Int_t newsp_num=0;
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nhitsperplane[ip] = 0;
+ }
+ // Sort Space Points hits by plane
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) { // All hits in SP
+ THcDCHit* hit=fSpacePoints[isp].hits[ihit];
+ // hit_order Make a hash
+ // hash(hit) = ihit;
+ Int_t ip = hit->GetPlaneIndex();
+ hits_plane[ip][nhitsperplane[ip]++] = hit;
+ }
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(nhitsperplane[ip] > 0) {
+ nplanes_hit++;
+ nsp_new *= nhitsperplane[ip];
+ if(nhitsperplane[ip] > 1) nplanes_mult++;
+ }
+ }
+ --nsp_new;
+ nsp_check=nsp_tot + nsp_new;
+ nplanes_single = nplanes_hit - nplanes_mult;
+
+ // Check if cloning conditions are met
+ Int_t ntot = 0;
+ if(nplanes_hit >= 4 && nplanes_mult < 4 && nplanes_mult >0
+ && nsp_check < 20) {
+
+ // Order planes by decreasing # of hits
+
+ Int_t maxplane[fNPlanes];
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ maxplane[ip] = ip;
+ }
+ // Sort by decreasing # of hits
+ for(Int_t ip1=0;ip1<fNPlanes-1;ip1++) {
+ for(Int_t ip2=ip1+1;ip2<fNPlanes;ip2++) {
+ if(nhitsperplane[maxplane[ip2]] > nhitsperplane[maxplane[ip1]]) {
+ Int_t temp = maxplane[ip1];
+ maxplane[ip1] = maxplane[ip2];
+ maxplane[ip2] = temp;
+ }
+ }
+ }
+
+ // First fill clones with 1 hit each from the 3 planes with the most hits
+ for(Int_t n1=0;n1<nhitsperplane[maxplane[0]];n1++) {
+ for(Int_t n2=0;n2<nhitsperplane[maxplane[1]];n2++) {
+ for(Int_t n3=0;n3<nhitsperplane[maxplane[2]];n3++) {
+ ntot++;
+ newsp_num = nsp_tot + ntot - 2; // ntot will be 2 for first new
+ if(n1==0 && n2==0 && n3==0) newsp_num = isp; // Copy over original SP
+ fSpacePoints[newsp_num].x = fSpacePoints[isp].x;
+ fSpacePoints[newsp_num].y = fSpacePoints[isp].y;
+ fSpacePoints[newsp_num].nhits = nplanes_hit;
+ fSpacePoints[newsp_num].ncombos = fSpacePoints[isp].ncombos;
+ fSpacePoints[newsp_num].hits[0] = hits_plane[maxplane[0]][n1];
+ fSpacePoints[newsp_num].hits[1] = hits_plane[maxplane[1]][n2];
+ fSpacePoints[newsp_num].hits[2] = hits_plane[maxplane[2]][n3];
+ fSpacePoints[newsp_num].hits[3] = hits_plane[maxplane[3]][0];
+ if(nhitsperplane[maxplane[4]] == 1)
+ fSpacePoints[newsp_num].hits[4] = hits_plane[maxplane[4]][0];
+ if(nhitsperplane[maxplane[5]] == 1)
+ fSpacePoints[newsp_num].hits[5] = hits_plane[maxplane[5]][0];
+ }
+ }
+ }
+#if 0
+ // Loop over clones and order hits in the same way as parent SP
+ // Why do we have to order the hits.
+ for(Int_t i=0;i<ntot;i++) {
+ Int_t newsp_num= nsp_tot + i;
+ if(i == 1) newsp_num = isp;
+ for(Int_t j=0;j<nplanes_hit;j++) {
+ for(Int_t k=0;k<nplanes_hit;k++) {
+ THcDCHit* hit1 = fSpacePointHits[newsp_num][j];
+ THcDCHit* hit2 = fSpacePointHits[newsp_num][k];
+ if(hit_order(hit1) > hit_order(hit2)) {
+ THcDCHit* temp = fSpacePoints[newsp_num].hits[k];
+ fSpacePoints[newsp_num].hits[k] = fSpacePoints[newsp_num].hits[j];
+ fSpacePoints[newsp_num].hits[j] = temp;
+ }
+ }
+ }
+ }
#endif
+ nsp_tot += (ntot-1);
+ } else {
+ ntot=1;
+ }
+ }
+ assert (nsp_tot > 0);
+ Int_t nadded=0;
+ if(nsp_tot <= 20) {
+ nadded = nsp_tot - fNSpacePoints;
+ fNSpacePoints = nsp_tot;
+ }
- return nhits;
+ // In Fortran, fill in zeros.
+ return(nadded);
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::ApplyCorrections( void )
+// HMS Specific?
+void THcDriftChamber::ChooseSingleHit()
{
- return(0);
+ // Look at all hits in a space point. If two hits are in the same plane,
+ // reject the one with the longer drift time.
+ Int_t goodhit[MAX_HITS_PER_POINT];
+
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t startnum = fSpacePoints[isp].nhits;
+
+ for(Int_t ihit=0;ihit<startnum;ihit++) {
+ goodhit[ihit] = 1;
+ }
+ // For each plane, mark all hits longer than the shortest drift time
+ for(Int_t ihit1=0;ihit1<startnum-1;ihit1++) {
+ THcDCHit* hit1 = fSpacePoints[isp].hits[ihit1];
+ Int_t plane1=hit1->GetPlaneIndex();
+ Double_t tdrift1 = hit1->GetTime();
+ for(Int_t ihit2=ihit1+1;ihit2<startnum;ihit2++) {
+ THcDCHit* hit2 = fSpacePoints[isp].hits[ihit2];
+ Int_t plane2=hit2->GetPlaneIndex();
+ Double_t tdrift2 = hit2->GetTime();
+ if(plane1 == plane2) {
+ if(tdrift1 > tdrift2) {
+ goodhit[ihit1] = 0;
+ } else {
+ goodhit[ihit2] = 0;
+ }
+ }
+ }
+ }
+ // Gather the remaining hits
+ Int_t finalnum = 0;
+ for(Int_t ihit=0;ihit<startnum;ihit++) {
+ if(goodhit[ihit] > 0) { // Keep this hit
+ if (ihit > finalnum) { // Move hit
+ fSpacePoints[isp].hits[finalnum++] = fSpacePoints[isp].hits[ihit];
+ } else {
+ finalnum++;
+ }
+ }
+ }
+ fSpacePoints[isp].nhits = finalnum;
+ }
}
+//_____________________________________________________________________________
+// Generic
+void THcDriftChamber::SelectSpacePoints()
+// This routine goes through the list of space_points and space_point_hits
+// found by find_space_points and only accepts those with
+// number of hits > min_hits
+// number of combinations > min_combos
+{
+ Int_t sp_count=0;
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ // Include fEasySpacePoint because ncombos not filled in
+ if(fSpacePoints[isp].ncombos >= fMinCombos || fEasySpacePoint) {
+ if(fSpacePoints[isp].nhits >= fMinHits) {
+ fSpacePoints[sp_count++] = fSpacePoints[isp];
+ }
+ }
+ }
+ fNSpacePoints = sp_count;
+}
+
+void THcDriftChamber::CorrectHitTimes()
+{
+ // Use the rough hit positions in the chambers to correct the drift time
+ // for hits in the space points.
+
+ // Assume all wires for a plane are read out on the same side (l/r or t/b).
+ // If the wire is closer to horizontal, read out left/right. If nearer
+ // vertical, assume top/bottom. (Note, this is not always true for the
+ // SOS u and v planes. They have 1 card each on the side, but the overall
+ // time offset per card will cancel much of the error caused by this. The
+ // alternative is to check by card, rather than by plane and this is harder.
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Double_t x = fSpacePoints[isp].x;
+ Double_t y = fSpacePoints[isp].y;
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) {
+ THcDCHit* hit = fSpacePoints[isp].hits[ihit];
+ THcDriftChamberPlane* plane=hit->GetWirePlane();
+
+ // How do we know this correction only gets applied once? Is
+ // it determined that a given hit can only belong to one space point?
+ Double_t time_corr = plane->GetReadoutX() ?
+ y*plane->GetReadoutCorr()/fWireVelocity :
+ x*plane->GetReadoutCorr()/fWireVelocity;
+
+ // cout << "Correcting hit " << hit << " " << plane->GetPlaneNum() << " " << isp << "/" << ihit << " " << x << "," << y << endl;
+ // Fortran ENGINE does not do this check, so hits can get "corrected"
+ // multiple times if they belong to multiple space points.
+ // To reproduce the precise ENGINE behavior, remove this if condition.
+ if(! hit->GetCorrectedStatus()) {
+ hit->SetTime(hit->GetTime() - plane->GetCentralTime()
+ + plane->GetDriftTimeSign()*time_corr);
+ hit->ConvertTimeToDist();
+ hit->SetCorrectedStatus(1);
+ }
+ }
+ }
+}
+UInt_t THcDriftChamber::Count1Bits(UInt_t x)
+// From "Hacker's Delight"
+{
+ x = x - ((x >> 1) & 0x55555555);
+ x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
+ x = x + (x >> 8);
+ x = x + (x >> 16);
+ return x & 0x0000003F;
+}
+
+//_____________________________________________________________________________
+// HMS Specific
+void THcDriftChamber::LeftRight()
+{
+ // For each space point,
+ // Fit stubs to all possible left-right combinations of drift distances
+ // and choose the set with the minimum chi**2.
+
+ for(Int_t isp=0; isp<fNSpacePoints; isp++) {
+ // Build a bit pattern of which planes are hit
+ Int_t nhits = fSpacePoints[isp].nhits;
+ UInt_t bitpat = 0; // Bit pattern of which planes are hit
+ Double_t minchi2 = 1.0e10;
+ Double_t tmp_minchi2;
+ Double_t minxp = 0.25;
+ Int_t hasy1 = -1;
+ Int_t hasy2 = -1;
+ Int_t plusminusknown[nhits];
+ Int_t plusminusbest[nhits];
+ Int_t plusminus[nhits]; // ENGINE makes this array float. Why?
+ Int_t tmp_plusminus[nhits];
+ Int_t plane_list[nhits];
+ Double_t stub[4];
+ Double_t tmp_stub[4];
+
+ if(nhits < 0) {
+ cout << "THcDriftChamber::LeftRight() nhits < 0" << endl;
+ } else if (nhits==0) {
+ cout << "THcDriftChamber::LeftRight() nhits = 0" << endl;
+ }
+ for(Int_t ihit=0;ihit < nhits;ihit++) {
+ THcDCHit* hit = fSpacePoints[isp].hits[ihit];
+ Int_t pindex = hit->GetPlaneIndex();
+ plane_list[ihit] = pindex;
+
+ bitpat |= 1<<pindex;
+
+ plusminusknown[ihit] = 0;
+ if(pindex == YPlaneInd) hasy1 = ihit;
+ if(pindex == YPlanePInd) hasy2 = ihit;
+ }
+ Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
+ if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes
+ if(fSpacePoints[isp].hits[hasy1]->GetPos() <=
+ fSpacePoints[isp].hits[hasy2]->GetPos()) {
+ plusminusknown[hasy1] = -1;
+ plusminusknown[hasy2] = 1;
+ } else {
+ plusminusknown[hasy1] = 1;
+ plusminusknown[hasy2] = -1;
+ }
+ }
+ if(nhits < 2) {
+ cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
+ } else if (nhits == 2) {
+ cout << "THcDriftChamber::LeftRight: numhits-2 = 0" << endl;
+ }
+ Int_t nplaneshit = Count1Bits(bitpat);
+ Int_t nplusminus = pow(2,nhits-2);
+
+ // Use bit value of integer word to set + or -
+ // Loop over all combinations of left right.
+ for(Int_t pmloop=0;pmloop<nplusminus;pmloop++) {
+ Int_t iswhit = 1;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ if(plusminusknown[ihit]!=0) {
+ plusminus[ihit] = plusminusknown[ihit];
+ } else {
+ // Max hits per point has to be less than 32.
+ if(pmloop & iswhit) {
+ plusminus[ihit] = 1;
+ } else {
+ plusminus[ihit] = -1;
+ }
+ iswhit <<= 1;
+ }
+ }
+ if (nplaneshit >= fNPlanes-1) {
+ Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits,
+ plane_list, bitpat, plusminus, stub);
+
+ //if(debugging)
+ //cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
+ // Take best chi2 IF x' of the stub agrees with x' as expected from x.
+ // Sometimes an incorrect x' gives a good chi2 for the stub, even though it is
+ // not the correct left/right combination for the real track.
+ // Rotate x'(=stub(3)) to hut coordinates and compare to x' expected from x.
+ // THIS ASSUMES STANDARD HMS TUNE!!!!, for which x' is approx. x/875.
+ if(chi2 < minchi2) {
+ if(stub[2]*fTanBeta[plane_list[0]]==-1.0) {
+ cout << "THcDriftChamber::LeftRight() Error 3" << endl;
+ }
+ Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
+ /(1+stub[2]*fTanBeta[plane_list[0]]);
+ // Tune depdendent. Definitely HMS specific
+ Double_t xp_expect = fSpacePoints[isp].x/875.0;
+ if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
+ minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ } else { // Record best stub failing angle cut
+ tmp_minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ tmp_plusminus[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ tmp_stub[i] = stub[i];
+ }
+ }
+ }
+ } else if (nplaneshit >= fNPlanes-2) { // Two planes missing
+ Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits,
+ plane_list, bitpat, plusminus, stub);
+ //if(debugging)
+ //cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
+ // Isn't this a bad idea, doing == with reals
+ if(stub[2]*fTanBeta[plane_list[0]] == -1.0) {
+ cout << "THcDriftChamber::LeftRight() Error 3" << endl;
+ }
+ Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
+ /(1+stub[2]*fTanBeta[plane_list[0]]);
+ if(TMath::Abs(xp_fit) <= minxp) {
+ minxp = TMath::Abs(xp_fit);
+ minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ }
+ } else {
+ cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
+ }
+ } // End loop of pm combinations
+
+ if(minchi2 > 9.9e9) { // No track passed angle cut
+ minchi2 = tmp_minchi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = tmp_plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = tmp_stub[i];
+ }
+
+ }
+
+ // Calculate final coordinate based on plusminusbest
+ // Update the hit positions in the space points
+ for(Int_t ihit; ihit<nhits; ihit++) {
+ fSpacePoints[isp].hits[ihit]->SetLeftRight(plusminusbest[ihit]);
+ }
+
+ // Stubs are calculated in rotated coordinate system
+ // (I think this rotates in case chambers not perpendicular to central ray)
+ Int_t pindex=plane_list[0];
+ if(fSpacePoints[isp].stub[2] - fTanBeta[pindex] == -1.0) {
+ cout << "THcDriftChamber::LeftRight(): stub3 error" << endl;
+ }
+ stub[2] = (fSpacePoints[isp].stub[2] - fTanBeta[pindex])
+ / (1.0 + fSpacePoints[isp].stub[2]*fTanBeta[pindex]);
+ if(fSpacePoints[isp].stub[2]*fSinBeta[pindex] == -fCosBeta[pindex]) {
+ cout << "THcDriftChamber::LeftRight(): stub4 error" << endl;
+ }
+ stub[3] = fSpacePoints[isp].stub[3]
+ / (fSpacePoints[isp].stub[2]*fSinBeta[pindex]+fCosBeta[pindex]);
+ stub[0] = fSpacePoints[isp].stub[0]*fCosBeta[pindex]
+ - fSpacePoints[isp].stub[0]*stub[2]*fSinBeta[pindex];
+ stub[1] = fSpacePoints[isp].stub[1]
+ - fSpacePoints[isp].stub[1]*stub[3]*fSinBeta[pindex];
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ }
+ // Option to print stubs
+}
+ // if(fAA3Inv.find(bitpat) != fAAInv.end()) { // Valid hit combination
//_____________________________________________________________________________
-Int_t THcDriftChamber::CoarseTrack( TClonesArray& /* tracks */ )
+Double_t THcDriftChamber::FindStub(Int_t nhits, THcDCHit** hits,
+ Int_t* plane_list, UInt_t bitpat,
+ Int_t* plusminus, Double_t* stub)
{
- // Calculation of coordinates of particle track cross point with scint
- // 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.);
+ // For a given combination of L/R, fit a stub to the space point
+ Double_t zeros[] = {0.0,0.0,0.0};
+ TVectorD TT; TT.Use(3, zeros);
+ TVectorD dstub;
+ Double_t dpos[3];
- ApplyCorrections();
+ for(Int_t ihit=0;ihit<nhits; ihit++) {
+ dpos[ihit] = hits[ihit]->GetPos() + plusminus[ihit]*hits[ihit]->GetdDist()
+ - fPsi0[plane_list[ihit]];
+ for(Int_t index=0;index<3;index++) {
+ TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
+ /fSigma[plane_list[ihit]];
+ }
+ }
+ dstub = (*fAA3Inv[bitpat]) * TT;
+
+ // Calculate Chi2. Remember one power of sigma is in fStubCoefs
+ stub[0] = dstub[0];
+ stub[1] = dstub[1];
+ stub[2] = dstub[2];
+ stub[3] = 0.0;
+ Double_t chi2=0.0;
+ for(Int_t ihit=0;ihit<nhits; ihit++) {
+ chi2 += pow( dpos[ihit]/fSigma[plane_list[ihit]]
+ - fStubCoefs[plane_list[ihit]][0]*stub[0]
+ - fStubCoefs[plane_list[ihit]][1]*stub[1]
+ - fStubCoefs[plane_list[ihit]][2]*stub[2]
+ , 2);
+ }
+ return(chi2);
+}
- return 0;
+//_____________________________________________________________________________
+THcDriftChamber::~THcDriftChamber()
+{
+ // Destructor. Remove variables from global list.
+
+ if( fIsSetup )
+ RemoveVariables();
+ if( fIsInit )
+ DeleteArrays();
+ if (fTrackProj) {
+ fTrackProj->Clear();
+ delete fTrackProj; fTrackProj = 0;
+ }
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::FineTrack( TClonesArray& tracks )
+void THcDriftChamber::DeleteArrays()
{
- // Reconstruct coordinates of particle track cross point with scintillator
- // plane, and copy the data into the following local data structure:
- //
- // Units of measurements are meters.
+ // Delete member arrays. Used by destructor.
- // Calculation of coordinates of particle track cross point with scint
- // plane in the detector coordinate system. For this, parameters of track
- // reconstructed in THaVDC::FineTrack() are used.
+}
- return 0;
+//_____________________________________________________________________________
+inline
+void THcDriftChamber::Clear( const Option_t* )
+{
+ // Reset per-event data.
+ // fNhits = 0;
+
+ // fTrackProj->Clear();
+ fNhits = 0;
+
+}
+
+//_____________________________________________________________________________
+Int_t THcDriftChamber::ApplyCorrections( void )
+{
+ return(0);
}
ClassImp(THcDriftChamber)
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index a973c466a25804993d2e8c3576d05cefd9c1eb56..db307dd06bfab3b204bf482e8c3820351d03da0b 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -7,52 +7,44 @@
// //
///////////////////////////////////////////////////////////////////////////////
-#include "THaTrackingDetector.h"
-#include "THcHitList.h"
-#include "THcRawDCHit.h"
+#include "THaSubDetector.h"
#include "THcDriftChamberPlane.h"
-#include "TMath.h"
+#include "TClonesArray.h"
+#include "TMatrixDSym.h"
+
+#include <map>
+
+#define MAX_SPACE_POINTS 50
+#define MAX_HITS_PER_POINT 20
+
+//#include "TMath.h"
//class THaScCalib;
class TClonesArray;
-class THcDriftChamber : public THaTrackingDetector, public THcHitList {
+class THcDriftChamber : public THaSubDetector {
public:
- THcDriftChamber( const char* name, const char* description = "",
- THaApparatus* a = NULL );
+ THcDriftChamber( const char* name, const char* description, Int_t chambernum,
+ THaDetectorBase* parent = NULL );
virtual ~THcDriftChamber();
- virtual Int_t Decode( const THaEvData& );
+ virtual Int_t Decode( const THaEvData& );
virtual EStatus Init( const TDatime& run_time );
- virtual Int_t CoarseTrack( TClonesArray& tracks );
- virtual Int_t FineTrack( TClonesArray& tracks );
+ virtual void AddPlane(THcDriftChamberPlane *plane);
virtual Int_t ApplyCorrections( void );
+ virtual void ProcessHits( void );
+ virtual Int_t FindSpacePoints( void ) ;
+ virtual void CorrectHitTimes( void ) ;
+
+ virtual void Clear( Option_t* opt="" );
// Int_t GetNHits() const { return fNhit; }
Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
const TClonesArray* GetTrackHits() const { return fTrackProj; }
- Int_t GetNWires(Int_t plane) const { return fNWires[plane-1];}
- Int_t GetNChamber(Int_t plane) const { return fNChamber[plane-1];}
- Int_t GetWireOrder(Int_t plane) const { return fWireOrder[plane-1];}
- Int_t GetPitch(Int_t plane) const { return fPitch[plane-1];}
- Int_t GetCentralWire(Int_t plane) const { return fCentralWire[plane-1];}
- Int_t GetTdcWinMin(Int_t plane) const { return fTdcWinMin[plane-1];}
- Int_t GetTdcWinMax(Int_t plane) const { return fTdcWinMax[plane-1];}
-
- Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
-
- Double_t GetNSperChan() const { return fNSperChan;}
-
- Double_t GetCenter(Int_t plane) const {
- Int_t chamber = GetNChamber(plane)-1;
- return
- fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
- fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
- }
// friend class THaScCalib;
THcDriftChamber(); // for ROOT I/O
@@ -61,65 +53,75 @@ protected:
// Calibration
// Per-event data
+ Int_t fNhits;
- // Potential Hall C parameters. Mostly here for demonstration
- Int_t fNPlanes;
- char** fPlaneNames;
- Int_t fNChambers;
+ Int_t fNPlanes; // Number of planes in the chamber
- Double_t fNSperChan; /* TDC bin size */
- Double_t fWireVelocity;
+ Int_t fChamberNum;
- // Each of these will be dimensioned with the number of chambers
- Double_t* fXCenter;
- Double_t* fYCenter;
-
- // Each of these will be dimensioned with the number of planes
- // A THcDriftChamberPlane class object will need to access the value for
- // its plane number. Should we have a Get method for each or
- Int_t* fTdcWinMin;
- Int_t* fTdcWinMax;
- Int_t* fCentralTime;
- Int_t* fNWires; // Number of wires per plane
- Int_t* fNChamber;
- Int_t* fWireOrder;
- Int_t* fDriftTimeSign;
-
- Double_t* fZPos;
- Double_t* fAlphaAngle;
- Double_t* fBetaAngle;
- Double_t* fGammaAngle;
- Double_t* fPitch;
- Double_t* fCentralWire;
- Double_t* fPlaneTimeZero;
-
- THcDriftChamberPlane** fPlanes; // List of plane objects
+ // HMS Specific
+ Int_t YPlaneInd; // Index of Yplane for this chamber
+ Int_t YPlanePInd; // Index of Yplanep for this chamber
+ Int_t YPlaneNum; // Absolute plane number of Yplane
+ Int_t YPlanePNum; // Absolute plane number of Yplanep
+
+ // Parameters
+ Int_t fMinHits; // Minimum hits required to do something
+ Int_t fMaxHits; // Maximum required to do something
+ Int_t fMinCombos; // Minimum # pairs in a space point
+ Int_t fRemove_Sppt_If_One_YPlane;
+ Double_t fWireVelocity;
+ Int_t fSmallAngleApprox;
+ Double_t fStubMaxXPDiff;
+
+ Double_t fXCenter;
+ Double_t fYCenter;
+ Double_t fSpacePointCriterion;
+ Double_t fSpacePointCriterion2;
+ Double_t* fSinBeta;
+ Double_t* fCosBeta;
+ Double_t* fTanBeta;
+ Double_t* fSigma;
+ Double_t* fPsi0;
+ Double_t** fStubCoefs;
+
+ THcDriftChamberPlane* fPlanes[20]; // 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 ClearEvent();
void DeleteArrays();
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
void Setup(const char* name, const char* description);
+ Int_t FindEasySpacePoint(Int_t yplane_hitind, Int_t yplanep_hitind);
+ Int_t FindHardSpacePoints(void);
+ Int_t DestroyPoorSpacePoints(void);
+ Int_t SpacePointMultiWire(void);
+ void ChooseSingleHit(void);
+ void SelectSpacePoints(void);
+ UInt_t Count1Bits(UInt_t x);
+ void LeftRight(void);
+ Double_t FindStub(Int_t nhits, THcDCHit** hits, Int_t* plane_list,
+ UInt_t bitpat, Int_t* plusminus, Double_t* stub);
+
+ THcDCHit* fHits[10000]; /* All hits for this chamber */
+ // A simple structure until we figure out what we are doing.
+ struct SpacePoint {
+ Double_t x;
+ Double_t y;
+ Int_t nhits;
+ Int_t ncombos;
+ THcDCHit* hits[MAX_HITS_PER_POINT];
+ Double_t stub[4];
+ };
+ SpacePoint fSpacePoints[MAX_SPACE_POINTS];
+ Int_t fNSpacePoints;
+ Int_t fEasySpacePoint; /* This event is an easy space point */
+
+ Double_t* stubcoef[4];
+ std::map<int,TMatrixDSym*> fAA3Inv;
ClassDef(THcDriftChamber,0) // Drift Chamber class
};
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 816d27a650ea46865c1bd42f0655f5751de1ae1d..bd708bf6c51a7557ef3d09880390f5ab2c637586 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -14,7 +14,7 @@
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THcHitList.h"
-#include "THcDriftChamber.h"
+#include "THcDC.h"
#include "THcHodoscope.h"
#include "TClass.h"
@@ -97,11 +97,11 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
// Retrieve parameters we need from parent class
- THcDriftChamber* fParent;
+ THcDC* fParent;
- fParent = (THcDriftChamber*) GetParent();
+ fParent = (THcDC*) GetParent();
// These are single variables here, but arrays in THcDriftChamber.
- fNChamber = fParent->GetNChamber(fPlaneNum);
+ fChamberNum = fParent->GetNChamber(fPlaneNum);
fNWires = fParent->GetNWires(fPlaneNum);
fWireOrder = fParent->GetWireOrder(fPlaneNum);
fPitch = fParent->GetPitch(fPlaneNum);
@@ -110,10 +110,72 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fTdcWinMax = fParent->GetTdcWinMax(fPlaneNum);
fPlaneTimeZero = fParent->GetPlaneTimeZero(fPlaneNum);
fCenter = fParent->GetCenter(fPlaneNum);
+ fCentralTime = fParent->GetCentralTime(fPlaneNum);
+ fDriftTimeSign = fParent->GetDriftTimeSign(fPlaneNum);
+ fSigma = fParent->GetSigma(fPlaneNum);
fNSperChan = fParent->GetNSperChan();
- cout << fPlaneNum << " " << fNWires << endl;
+ // Calculate Geometry Constants
+ // Do we want to move all this to the Chamber of DC Package leve
+ // as that is where these things will be needed?
+ Double_t z0 = fParent->GetZPos(fPlaneNum);
+ Double_t alpha = fParent->GetAlphaAngle(fPlaneNum);
+ Double_t beta = fParent->GetBetaAngle(fPlaneNum);
+ fBeta = beta;
+ Double_t gamma = fParent->GetGammaAngle(fPlaneNum);
+ Double_t cosalpha = TMath::Cos(alpha);
+ Double_t sinalpha = TMath::Sin(alpha);
+ Double_t cosbeta = TMath::Cos(beta);
+ Double_t sinbeta = TMath::Sin(beta);
+ Double_t cosgamma = TMath::Cos(gamma);
+ Double_t singamma = TMath::Sin(gamma);
+
+ Double_t hzchi = -cosalpha*sinbeta + sinalpha*cosbeta*singamma;
+ Double_t hzpsi = sinalpha*sinbeta + cosalpha*cosbeta*singamma;
+ Double_t hxchi = -cosalpha*cosbeta - sinalpha*sinbeta*singamma;
+ Double_t hxpsi = sinalpha*cosbeta - cosalpha*sinbeta*singamma;
+ Double_t hychi = sinalpha*cosgamma;
+ Double_t hypsi = cosalpha*cosgamma;
+ Double_t stubxchi = -cosalpha;
+ Double_t stubxpsi = sinalpha;
+ Double_t stubychi = sinalpha;
+ Double_t stubypsi = cosalpha;
+
+ if(cosalpha <= 0.707) { // x-like wire, need dist from x=0 line
+ fReadoutX = 1;
+ fReadoutCorr = 1/sinalpha;
+ } else {
+ fReadoutX = 0;
+ fReadoutCorr = 1/cosalpha;
+ }
+
+ Double_t sumsqupsi = hzpsi*hzpsi+hxpsi*hxpsi+hypsi*hypsi;
+ Double_t sumsquchi = hzchi*hzchi+hxchi*hxchi+hychi*hychi;
+ Double_t sumcross = hzpsi*hzchi + hxpsi*hxchi + hypsi*hychi;
+ Double_t denom1 = sumsqupsi*sumsquchi-sumcross*sumcross;
+ fPsi0 = (-z0*hzpsi*sumsquchi
+ +z0*hzchi*sumcross) / denom1;
+ Double_t hchi0 = (-z0*hzchi*sumsqupsi
+ +z0*hzpsi*sumcross) / denom1;
+ Double_t hphi0 = TMath::Sqrt(pow(z0+hzpsi*fPsi0+hzchi*hchi0,2)
+ + pow(hxpsi*fPsi0+hxchi*hchi0,2)
+ + pow(hypsi*fPsi0+hychi*hchi0,2) );
+ if(z0 < 0.0) hphi0 = -hphi0;
+
+ Double_t denom2 = stubxpsi*stubychi - stubxchi*stubypsi;
+
+ // Why are there 4, but only 3 used?
+ fStubCoef[0] = stubychi/(fSigma*denom2); // sin(a)/sigma
+ fStubCoef[1] = -stubxchi/(fSigma*denom2); // cos(a)/sigma
+ fStubCoef[2] = hphi0*fStubCoef[0]; // z0*sin(a)/sig
+ fStubCoef[3] = hphi0*fStubCoef[1]; // z0*cos(a)/sig
+
+ fXsp = hychi/denom2; // sin(a)
+ fYsp = -hxchi/denom2; // cos(a)
+
+
+ cout << fPlaneNum << " " << fNWires << " " << fWireOrder << endl;
fTTDConv = new THcDCLookupTTDConv(DriftMapFirstBin,fPitch/2,DriftMapBinSize,
NumDriftMapBins,DriftMap);
@@ -242,7 +304,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
- rawtdc*fNSperChan + fPlaneTimeZero;
// How do we get this start time from the hodoscope to here
// (or at least have it ready by coarse process)
- new( (*fHits)[nextHit++] ) THcDCHit(wire, rawtdc, time);
+ new( (*fHits)[nextHit++] ) THcDCHit(wire, rawtdc, time, this);
}
wire_last = wireNum;
}
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index 7307a4acf649f75679223878d27d2607747c3798..475e9bc9d1d4038a4f1f8bede46e535f41d940b0 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -48,6 +48,24 @@ class THcDriftChamberPlane : public THaSubDetector {
{ assert( i>=1 && i<=GetNWires() );
return (THcDCWire*)fWires->UncheckedAt(i-1); }
+ Int_t GetNHits() const { return fHits->GetLast()+1; }
+ TClonesArray* GetHits() const { return fHits; }
+
+ Int_t GetPlaneNum() const { return fPlaneNum; }
+ Int_t GetChamberNum() const { return fChamberNum; }
+ void SetPlaneIndex(Int_t index) { fPlaneIndex = index; }
+ Int_t GetPlaneIndex() { return fPlaneIndex; }
+ Double_t GetXsp() const { return fXsp; }
+ Double_t GetYsp() const { return fYsp; }
+ Int_t GetReadoutX() { return fReadoutX; }
+ Double_t GetReadoutCorr() { return fReadoutCorr; }
+ Double_t GetCentralTime() { return fCentralTime; }
+ Int_t GetDriftTimeSign() { return fDriftTimeSign; }
+ Double_t GetBeta() { return fBeta; }
+ Double_t GetSigma() { return fSigma; }
+ Double_t GetPsi0() { return fPsi0; }
+ Double_t* GetStubCoef() { return fStubCoef; }
+
protected:
TClonesArray* fParentHitList;
@@ -56,7 +74,8 @@ class THcDriftChamberPlane : public THaSubDetector {
TClonesArray* fWires;
Int_t fPlaneNum;
- Int_t fNChamber;
+ Int_t fPlaneIndex; /* Index of this plane within it's chamber */
+ Int_t fChamberNum;
Int_t fNWires;
Int_t fWireOrder;
Int_t fTdcWinMin;
@@ -64,6 +83,17 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t fPitch;
Double_t fCentralWire;
Double_t fPlaneTimeZero;
+ Double_t fXsp;
+ Double_t fYsp;
+ Double_t fSigma;
+ Double_t fPsi0;
+ Double_t fStubCoef[4];
+ Double_t fBeta;
+
+ Int_t fReadoutX;
+ Double_t fReadoutCorr;
+ Double_t fCentralTime;
+ Int_t fDriftTimeSign;
Double_t fCenter;
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index 4d09ccf2560968b566b79fcccd8ae438bc95ae92..4a7f365f3cb78d3258b051d56060fe427bcde181 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -73,15 +73,29 @@ void THcHodoscope::Setup(const char* name, const char* description)
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
+ string planenamelist;
DBRequest listextra[]={
{"hodo_num_planes", &fNPlanes, kInt},
+ {"hodo_plane_names",&planenamelist, kString},
{0}
};
- fNPlanes = 4; // Default if not defined
+ //fNPlanes = 4; // Default if not defined
gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);
- // Plane names
+ cout << "Plane Name List : " << planenamelist << endl;
+
+ vector<string> plane_names = vsplit(planenamelist);
+ // Plane names
+ if(plane_names.size() != (UInt_t) fNPlanes) {
+ cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
+ // Should quit. Is there an official way to quit?
+ }
fPlaneNames = new char* [fNPlanes];
+ for(Int_t i=0;i<fNPlanes;i++) {
+ fPlaneNames[i] = new char[plane_names[i].length()];
+ strcpy(fPlaneNames[i], plane_names[i].c_str());
+ }
+ /* fPlaneNames = new char* [fNPlanes];
for(Int_t i=0;i<fNPlanes;i++) {fPlaneNames[i] = new char[3];}
// Should get the plane names from parameters.
// could try this: grep _zpos PARAM/hhodo.pos | sed 's/\_/\ /g' | awk '{print $2}'
@@ -89,16 +103,15 @@ void THcHodoscope::Setup(const char* name, const char* description)
strcpy(fPlaneNames[1],"1y");
strcpy(fPlaneNames[2],"2x");
strcpy(fPlaneNames[3],"2y");
-
+ */
// Probably shouldn't assume that description is defined
- char* desc = new char[strlen(description)+50];
+ char* desc = new char[strlen(description)+100];
fPlanes = new THcScintillatorPlane* [fNPlanes];
for(Int_t i=0;i < fNPlanes;i++) {
strcpy(desc, description);
strcat(desc, " Plane ");
strcat(desc, fPlaneNames[i]);
- fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i+1,fNPlanes, this);
- //fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i+1, this);
+ fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i+1,fNPlanes,this); // Number planes starting from zero!!
cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
}
}
@@ -304,7 +317,7 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
strcpy(parname,prefix);
strcat(parname,"scin_");
Int_t plen=strlen(parname);
-
+ cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;
fNPaddle = new Int_t [fNPlanes];
// fSpacing = new Double_t [fNPlanes];
//fCenter = new Double_t* [fNPlanes];
@@ -393,6 +406,11 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
{"hodo_pos_ped_limit",&fHodoPosPedLimit[0],kInt,fMaxHodoScin},
{"hodo_neg_ped_limit",&fHodoNegPedLimit[0],kInt,fMaxHodoScin},
{"tofusinginvadc",&fTofUsingInvAdc,kInt},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ if (fTofUsingInvAdc) {
+ DBRequest list2[]={
{"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin},
{"hodo_neg_invadc_offset",&fHodoNegInvAdcOffset[0],kDouble,fMaxHodoScin},
{"hodo_pos_invadc_linear",&fHodoPosInvAdcLinear[0],kDouble,fMaxHodoScin},
@@ -401,7 +419,8 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
{"hodo_neg_invadc_adc",&fHodoNegInvAdcAdc[0],kDouble,fMaxHodoScin},
{0}
};
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ gHcParms->LoadParmValues((DBRequest*)&list2,prefix);
+ };
if (fDebug >=1) {
cout <<"******* Testing Hodoscope Parameter Reading ***\n";
cout<<"StarTimeCenter = "<<fStartTimeCenter<<endl;
@@ -455,7 +474,9 @@ Int_t THcHodoscope::DefineVariables( EMode mode )
//...
// hnegtdc4 HMS s2y- TDC hits
- // RVarDef vars[] = {
+ RVarDef vars[] = {
+ {"starttime","Hodoscope Start Time","fStartTime"},
+ {"hgoodstarttime","Hodoscope Good Start Time","fGoodStartTime"},
// { "nlthit", "Number of Left paddles TDC times", "fLTNhit" },
// { "nrthit", "Number of Right paddles TDC times", "fRTNhit" },
// { "nlahit", "Number of Left paddles ADCs amps", "fLANhit" },
@@ -484,10 +505,10 @@ Int_t THcHodoscope::DefineVariables( EMode mode )
// { "trpath", "TRCS pathlen of track to det plane","fTrackProj.THaTrackProj.fPathl" },
// { "trdx", "track deviation in x-position (m)", "fTrackProj.THaTrackProj.fdX" },
// { "trpad", "paddle-hit associated with track", "fTrackProj.THaTrackProj.fChannel" },
- // { 0 }
- // };
- // return DefineVarsFromList( vars, mode );
- return kOK;
+ { 0 }
+ };
+ return DefineVarsFromList( vars, mode );
+ // return kOK;
}
//_____________________________________________________________________________
@@ -594,6 +615,7 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
Int_t nexthit = 0;
Int_t nfptimes=0;
+ fStartTime=0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
// nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
// GN: select only events that have reasonable TDC values to start with
@@ -601,9 +623,8 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
if (fPlanes[ip]->GetNScinHits()>0) {
fPlanes[ip]->PulseHeightCorrection();
- ///cout <<"Plane "<<ip<<" number of fpTimes = "<<fPlanes[ip]->GetFpTimeHits()<<endl;
- for (Int_t ifptimes=0;ifptimes<fPlanes[ip]->GetFpTimeHits();ifptimes++) {
- fStartTime=fStartTime+fPlanes[ip]->GetFpTime(ifptimes);
+ if (TMath::Abs(fPlanes[ip]->GetFpTime()-fStartTimeCenter)<=fStartTimeSlop) {
+ fStartTime=fStartTime+fPlanes[ip]->GetFpTime();
nfptimes++;
}
}
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index f1ccda77f1890d99e1453e9b6928a2275943c675..6595540b5371e11414e1b3a713e0194211f108c6 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -36,13 +36,16 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fNegTDCHits = new TClonesArray("THcSignalHit",16);
fPosADCHits = new TClonesArray("THcSignalHit",16);
fNegADCHits = new TClonesArray("THcSignalHit",16);
+ frPosTDCHits = new TClonesArray("THcSignalHit",16);
+ frNegTDCHits = new TClonesArray("THcSignalHit",16);
+ frPosADCHits = new TClonesArray("THcSignalHit",16);
+ frNegADCHits = new TClonesArray("THcSignalHit",16);
fPlaneNum = planenum;
fTotPlanes = planenum;
fNScinHits = 0;
//
fMaxHits=53;
- fpTimeHits=0;
- fpTime = new Double_t [fMaxHits];
+ fpTime = -1.e5;
}
//______________________________________________________________________________
THcScintillatorPlane::THcScintillatorPlane( const char* name,
@@ -57,13 +60,16 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fNegTDCHits = new TClonesArray("THcSignalHit",16);
fPosADCHits = new TClonesArray("THcSignalHit",16);
fNegADCHits = new TClonesArray("THcSignalHit",16);
+ frPosTDCHits = new TClonesArray("THcSignalHit",16);
+ frNegTDCHits = new TClonesArray("THcSignalHit",16);
+ frPosADCHits = new TClonesArray("THcSignalHit",16);
+ frNegADCHits = new TClonesArray("THcSignalHit",16);
fPlaneNum = planenum;
fTotPlanes = totplanes;
fNScinHits = 0;
//
fMaxHits=53;
- fpTimeHits=0;
- fpTime = new Double_t [fMaxHits];
+ fpTime = -1.e5;
}
@@ -75,7 +81,10 @@ THcScintillatorPlane::~THcScintillatorPlane()
delete fNegTDCHits;
delete fPosADCHits;
delete fNegADCHits;
- delete fpTime;
+ delete frPosTDCHits;
+ delete frNegTDCHits;
+ delete frPosADCHits;
+ delete frNegADCHits;
}
//______________________________________________________________________________
@@ -144,14 +153,14 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
//
// Based on the signs of these quantities in the .pos file the correspondence
// should be bot=>left and top=>right when comparing x and y-type scintillators
- char *tmpleft, *tmpright;
+ char tmpleft[6], tmpright[6];
if (fPlaneNum==1 || fPlaneNum==3) {
- tmpleft="left";
- tmpright="right";
+ strcpy(tmpleft,"left");
+ strcpy(tmpright,"right");
}
else {
- tmpleft="bot";
- tmpright="top";
+ strcpy(tmpleft,"bot");
+ strcpy(tmpright,"top");
}
Double_t tmpdouble[fTotPlanes];
@@ -210,13 +219,15 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode )
// Register variables in global list
RVarDef vars[] = {
{"postdchits", "List of Positive TDC hits",
- "fPosTDCHits.THcSignalHit.GetPaddleNumber()"},
+ "frPosTDCHits.THcSignalHit.GetPaddleNumber()"},
{"negtdchits", "List of Negative TDC hits",
- "fNegTDCHits.THcSignalHit.GetPaddleNumber()"},
+ "frNegTDCHits.THcSignalHit.GetPaddleNumber()"},
{"posadchits", "List of Positive ADC hits",
- "fPosADCHits.THcSignalHit.GetPaddleNumber()"},
+ "frPosADCHits.THcSignalHit.GetPaddleNumber()"},
{"negadchits", "List of Negative ADC hits",
- "fNegADCHits.THcSignalHit.GetPaddleNumber()"},
+ "frNegADCHits.THcSignalHit.GetPaddleNumber()"},
+ {"fptime", "Time at focal plane",
+ "GetFpTime()"},
{ 0 }
};
@@ -232,6 +243,10 @@ void THcScintillatorPlane::Clear( Option_t* )
fNegTDCHits->Clear();
fPosADCHits->Clear();
fNegADCHits->Clear();
+ frPosTDCHits->Clear();
+ frNegTDCHits->Clear();
+ frPosADCHits->Clear();
+ frNegADCHits->Clear();
}
//_____________________________________________________________________________
@@ -271,6 +286,16 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
// if the actual ADC is larger than the pedestal value we subtract!!
Double_t mintdc, maxtdc;
+ //raw
+ Int_t nrPosTDCHits=0;
+ Int_t nrNegTDCHits=0;
+ Int_t nrPosADCHits=0;
+ Int_t nrNegADCHits=0;
+ frPosTDCHits->Clear();
+ frNegTDCHits->Clear();
+ frPosADCHits->Clear();
+ frNegADCHits->Clear();
+ //stripped
Int_t nPosTDCHits=0;
Int_t nNegTDCHits=0;
Int_t nPosADCHits=0;
@@ -280,7 +305,6 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fNegTDCHits->Clear();
fPosADCHits->Clear();
fNegADCHits->Clear();
-
Int_t nrawhits = rawhits->GetLast()+1;
// cout << "THcScintillatorPlane::ProcessHits " << fPlaneNum << " " << nexthit << "/" << nrawhits << endl;
mintdc=((THcHodoscope *)GetParent())->GetTdcMin();
@@ -291,34 +315,28 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(hit->fPlane > fPlaneNum) {
break;
}
-
+ Int_t padnum=hit->fCounter;
+ Int_t index=padnum-1;
+ if (hit->fTDC_pos > 0) ((THcSignalHit*) frPosTDCHits->ConstructedAt(nrPosTDCHits++))->Set(padnum, hit->fTDC_pos);
+ if (hit->fTDC_neg > 0) ((THcSignalHit*) frNegTDCHits->ConstructedAt(nrNegTDCHits++))->Set(padnum, hit->fTDC_neg);
+ if ((hit->fADC_pos-fPosPed[index]) >= 50) ((THcSignalHit*) frPosADCHits->ConstructedAt(nrPosADCHits++))->Set(padnum, hit->fADC_pos-fPosPed[index]);
+ if ((hit->fADC_neg-fNegPed[index]) >= 50) ((THcSignalHit*) frNegADCHits->ConstructedAt(nrNegADCHits++))->Set(padnum, hit->fADC_neg-fNegPed[index]);
// check TDC values
if (((hit->fTDC_pos >= mintdc) && (hit->fTDC_pos <= maxtdc)) ||
((hit->fTDC_neg >= mintdc) && (hit->fTDC_neg <= maxtdc))) {
-
//TDC positive hit
THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
- sighit->Set(hit->fCounter, hit->fTDC_pos);
+ sighit->Set(padnum, hit->fTDC_pos);
// TDC negative hit
THcSignalHit *sighit2 = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
- sighit2->Set(hit->fCounter, hit->fTDC_neg);
+ sighit2->Set(padnum, hit->fTDC_neg);
// ADC positive hit
- /// if(hit->fADC_pos > 0) {
- THcSignalHit *sighit3 = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
- sighit3->Set(hit->fCounter, hit->fADC_pos-fPosPed[ihit]);
- ///} else {
- /// cout<<"Skipping ADC_pos "<<hit->fADC_pos<<endl;
- /// }
+ THcSignalHit *sighit3 = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ sighit3->Set(padnum, hit->fADC_pos-fPosPed[index]);
// ADC negative hit
- /// if(hit->fADC_neg > 0) {
- // cout <<"adc neg hit!!\n";
- THcSignalHit *sighit4 = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
- sighit4->Set(hit->fCounter, hit->fADC_neg-fNegPed[ihit]);
- ///} else {
- ///cout<<"Skipping ADC_neg "<<hit->fADC_neg<<endl;
- ///}
- // cout <<"test "<<fNHits<<endl;
- fNScinHits=fNScinHits++;
+ THcSignalHit *sighit4 = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+ sighit4->Set(padnum, hit->fADC_neg-fNegPed[index]);
+ fNScinHits++;
}
else {
//cout <<"pos TDC "<<hit->fTDC_pos<<" "<<mintdc<<" "<<maxtdc<<endl;
@@ -364,7 +382,7 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
for (i=0;i<200;i++) {
timehist[i]=0;
}
- for (i=0;i<53;i++) {
+ for (i=0;i<fMaxHits;i++) {
keep_pos[i]=kFALSE;
keep_neg[i]=kFALSE;
two_good_times[i]=kFALSE;
@@ -376,8 +394,8 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
toftolerance=((THcHodoscope *)GetParent())->GetTofTolerance();
// hbeta_pcent=(TH((THcHodoscope *)GetParent())->GetParent()
// Horrible hack until I find out where to get the central beta from momentum!! GN
- hbeta_pcent=0.99776;
-
+ hbeta_pcent=1.0;
+ fpTime=-1e5;
for (i=0;i<fNScinHits;i++) {
if ((((THcSignalHit*) fPosTDCHits->At(i))->GetData()>=mintdc) &&
(((THcSignalHit*) fPosTDCHits->At(i))->GetData()<=maxtdc) &&
@@ -385,18 +403,29 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
(((THcSignalHit*) fNegTDCHits->At(i))->GetData()<=maxtdc)) {
pos_ph[i]=((THcSignalHit*) fPosADCHits->At(i))->GetData();
postime[i]=((THcSignalHit*) fPosTDCHits->At(i))->GetData()*tdctotime;
- j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber();
+ j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber()-1;
index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1)));
postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index);
+
neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData();
negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime;
- j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber();
+ j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber()-1;
index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1)));
negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegTimeOffset(index);
+
+ // ***************
+ /// cout <<"hcana i = "<<i<<endl;
+ /// cout <<"hcana tdc_pos = "<<((THcSignalHit*) fPosTDCHits->At(i))->GetData()<<endl;
+ ///cout <<"hcana pos_ph = "<<pos_ph[i]<<endl;
+ // cout <<"hcana pos_phc_coeff = "<<((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)<<endl;
+ // cout <<"hcana postime = "<<postime[i]<<endl;
+ //cout <<"hcana negtime = "<<negtime[i]<<endl;
+ //*************
+
// Find hit position. If postime larger, then hit was nearer negative side.
dist_from_center=0.5*(negtime[i]-postime[i])*((THcHodoscope *)GetParent())->GetHodoVelLight(index);
scint_center=0.5*(fPosLeft+fPosRight);
@@ -405,6 +434,7 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
hit_position=TMath::Max(hit_position,fPosRight);
postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
+
time_pos[i]=postime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
time_neg[i]=negtime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
nfound++;
@@ -421,7 +451,7 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
}
}
// Find the bin with most hits
- jmax=-1;
+ jmax=0;
maxhit=0;
for (i=0;i<200;i++) {
if (timehist[i]>maxhit) {
@@ -451,14 +481,13 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
two_good_times[i]=kTRUE;
}
} // end of loop that finds tube setting time
- fpTimeHits=0;
for (i=0;i<fNScinHits;i++) {
if (two_good_times[i]) { // both tubes fired
// correct time for everything except veloc. correction in order
// to find hit location from difference in TDC.
pos_ph[i]=((THcSignalHit*) fPosADCHits->At(i))->GetData();
postime[i]=((THcSignalHit*) fPosTDCHits->At(i))->GetData()*tdctotime;
- j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber();
+ j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber()-1;
index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1)));
@@ -466,7 +495,7 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
//
neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData();
negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime;
- j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber();
+ j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber()-1;
index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1)));
@@ -480,13 +509,10 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
scin_corrected_time[i]=0.5*(postime[i]+negtime[i]);
- fpTime[i]=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
- /// cout <<"fptime ["<<i<<"]= "<<fpTime[i]<<endl;
- fpTimeHits++;
+ fpTime=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
}
else { // only one tube fired
scin_corrected_time[i]=0.0; // not a very good "flag" but there is the logical two_good_hits...
- fpTime[i]=0.0; // NOTE: we're not incrementing fpTimeHits. In principle these two lines might be deleted
}
}
// Start time calculation, assume xp=yp=0 radians. project all
@@ -552,7 +578,7 @@ void THcScintillatorPlane::CalculatePedestals( )
fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
fNegThresh[i] = fNegPed[i] + 15;
- // cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
+ // cout <<"Pedestals "<< i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
}
// cout << " " << endl;
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index f47b1ebdd721d7c954a84663b42abb492e349230..cb832fa7dd0975bf1eb0d064e6e6cb91a2e60cfb 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -53,13 +53,16 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t GetPosRight();
Double_t GetPosOffset();
Double_t GetPosCenter(Int_t PaddleNo); // here we're counting from zero!
- Double_t GetFpTime(Int_t index) { return fpTime[index];};
- Int_t GetFpTimeHits() { return fpTimeHits;};
+ Double_t GetFpTime() { return fpTime;};
TClonesArray* fParentHitList;
protected:
+ TClonesArray* frPosTDCHits;
+ TClonesArray* frNegTDCHits;
+ TClonesArray* frPosADCHits;
+ TClonesArray* frNegADCHits;
TClonesArray* fPosTDCHits;
TClonesArray* fNegTDCHits;
TClonesArray* fPosADCHits;
@@ -103,8 +106,8 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t *fNegThresh;
//
- Int_t fpTimeHits; // number of entries in *fpTime
- Double_t *fpTime; // array with fpTimes for this scintillator plane
+ Double_t fpTime; // the original code only has one fpTime per plane!
+
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
diff --git a/src/THcSignalHit.h b/src/THcSignalHit.h
index da9d27bbacafefd96e8c1d401a2fd860adcfcb2e..6af80a6e68a203261d306336022d015a237a71cf 100644
--- a/src/THcSignalHit.h
+++ b/src/THcSignalHit.h
@@ -22,6 +22,8 @@ class THcSignalHit : public TObject {
virtual void Set(Int_t paddle, Int_t data)
{ fPaddleNumber=paddle; fData=data; }
+ virtual void Set(Int_t paddle, Double_t data)
+ { fPaddleNumber=paddle; fData=data; }
private:
Int_t fPaddleNumber;