diff --git a/.clang-format b/.clang-format
index 0cbccb84d2af16059203b8e1f849b37a0c88e4d0..abc048740674e237d1e012b45148fda7a7592f4d 100644
--- a/.clang-format
+++ b/.clang-format
@@ -91,8 +91,7 @@ PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 60
PointerAlignment: Left
RawStringFormats:
- - Delimiter: pb
- Language: TextProto
+ - Language: TextProto
BasedOnStyle: google
ReflowComments: true
SortIncludes: true
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a0144f0608c7ae8eb579d798e4af69eeeb0db12f..66a089314e07ba10c5daf1d3bdbbf906254b48de 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
cmake_minimum_required(VERSION 3.5)
-project(hcana VERSION 0.90 LANGUAGES CXX)
+project(hcana VERSION 0.95 LANGUAGES CXX)
option(HCANA_BUILTIN_PODD "Use built-in Podd submodule (default: YES)" OFF)
diff --git a/src/THcCoinTime.cxx b/src/THcCoinTime.cxx
index 3e520f3176305643c234d2cb7ab57e2cbbd97e9d..653e6e2d5b69261e866dbe2dfe694488ab5653ca 100644
--- a/src/THcCoinTime.cxx
+++ b/src/THcCoinTime.cxx
@@ -168,18 +168,28 @@ Int_t THcCoinTime::DefineVariables( EMode mode )
const RVarDef vars[] = {
{"epCoinTime_ROC1", "ROC1 Corrected ep Coincidence Time", "fROC1_epCoinTime"},
{"epCoinTime_ROC2", "ROC2 Corrected ep Coincidence Time", "fROC2_epCoinTime"},
+ {"epCoinTime_TRIG1", "TRIG1 Corrected ep Coincidence Time", "fTRIG1_epCoinTime"},
+ {"epCoinTime_TRIG4", "TRIG4 Corrected ep Coincidence Time", "fTRIG4_epCoinTime"},
{"eKCoinTime_ROC1", "ROC1 Corrected eK Coincidence Time", "fROC1_eKCoinTime"},
{"eKCoinTime_ROC2", "ROC2 Corrected eK Coincidence Time", "fROC2_eKCoinTime"},
+ {"eKCoinTime_TRIG1", "TRIG1 Corrected eK Coincidence Time", "fTRIG1_eKCoinTime"},
+ {"eKCoinTime_TRIG4", "TRIG4 Corrected eK Coincidence Time", "fTRIG4_eKCoinTime"},
{"ePiCoinTime_ROC1", "ROC1 Corrected ePi Coincidence Time", "fROC1_ePiCoinTime"},
{"ePiCoinTime_ROC2", "ROC2 Corrected ePi Coincidence Time", "fROC2_ePiCoinTime"},
+ {"ePiCoinTime_TRIG1", "TRIG1 Corrected ePi Coincidence Time", "fTRIG1_ePiCoinTime"},
+ {"ePiCoinTime_TRIG4", "TRIG4 Corrected ePi Coincidence Time", "fTRIG4_ePiCoinTime"},
{"ePositronCoinTime_ROC1", "ROC1 Corrected e-Positorn Coincidence Time", "fROC1_ePosCoinTime"},
{"ePositronCoinTime_ROC2", "ROC2 Corrected e-Positron Coincidence Time", "fROC2_ePosCoinTime"},
+ {"ePositronCoinTime_TRIG1", "TRIG1 Corrected e-Positorn Coincidence Time", "fTRIG1_ePosCoinTime"},
+ {"ePositronCoinTime_TRIG4", "TRIG4 Corrected e-Positron Coincidence Time", "fTRIG4_ePosCoinTime"},
{"CoinTime_RAW_ROC1", "ROC1 RAW Coincidence Time", "fROC1_RAW_CoinTime"},
{"CoinTime_RAW_ROC2", "ROC2 RAW Coincidence Time", "fROC2_RAW_CoinTime"},
+ {"CoinTime_RAW_TRIG1", "TRIG1 RAW Coincidence Time", "fTRIG1_RAW_CoinTime"},
+ {"CoinTime_RAW_TRIG4", "TRIG4 RAW Coincidence Time", "fTRIG4_RAW_CoinTime"},
{"DeltaSHMSPathLength","DeltaSHMSpathLength (cm)","DeltaSHMSpathLength"},
{"DeltaHMSPathLength", "DeltaHMSpathLength (cm)","DeltaHMSpathLength"},
{"had_coinCorr_Positron", "", "had_coinCorr_Positron"},
@@ -294,6 +304,8 @@ Int_t THcCoinTime::Process( const THaEvData& evdata )
//Raw, Uncorrected Coincidence Time
fROC1_RAW_CoinTime = (pTRIG1_TdcTime_ROC1 + SHMS_FPtime) - (pTRIG4_TdcTime_ROC1 + HMS_FPtime);
fROC2_RAW_CoinTime = (pTRIG1_TdcTime_ROC2 + SHMS_FPtime) - (pTRIG4_TdcTime_ROC2 + HMS_FPtime);
+ fTRIG1_RAW_CoinTime = (pTRIG1_TdcTime_ROC1 + SHMS_FPtime) - (pTRIG1_TdcTime_ROC2 + HMS_FPtime);
+ fTRIG4_RAW_CoinTime = (pTRIG4_TdcTime_ROC1 + SHMS_FPtime) - (pTRIG4_TdcTime_ROC2 + HMS_FPtime);
//Corrected Coincidence Time for ROC1/ROC2 (ROC1 Should be identical to ROC2)
@@ -301,19 +313,26 @@ Int_t THcCoinTime::Process( const THaEvData& evdata )
//PROTON
fROC1_epCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset;
fROC2_epCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset;
+ fTRIG1_epCoinTime = fTRIG1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset;
+ fTRIG4_epCoinTime = fTRIG4_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset;
//KAON
fROC1_eKCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
fROC2_eKCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
+ fTRIG1_eKCoinTime = fTRIG1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
+ fTRIG4_eKCoinTime = fTRIG4_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
//PION
fROC1_ePiCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;
fROC2_ePiCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;
+ fTRIG1_ePiCoinTime = fTRIG1_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;
+ fTRIG4_ePiCoinTime = fTRIG4_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;
//POSITRON
fROC1_ePosCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Positron) - eHad_CT_Offset ;
fROC2_ePosCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Positron) - eHad_CT_Offset;
-
+ fTRIG1_ePosCoinTime = fTRIG1_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Positron) - eHad_CT_Offset ;
+ fTRIG4_ePosCoinTime = fTRIG4_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Positron) - eHad_CT_Offset;
diff --git a/src/THcCoinTime.h b/src/THcCoinTime.h
index 8cc33411f8c9db3f133589f60eb884b43cb6612b..ecaf5008c6365d9c317deaa3a28b32ab583826f9 100644
--- a/src/THcCoinTime.h
+++ b/src/THcCoinTime.h
@@ -85,19 +85,29 @@ public:
Double_t fROC1_RAW_CoinTime;
Double_t fROC2_RAW_CoinTime;
+ Double_t fTRIG1_RAW_CoinTime;
+ Double_t fTRIG4_RAW_CoinTime;
Double_t fROC1_epCoinTime;
Double_t fROC2_epCoinTime;
+ Double_t fTRIG1_epCoinTime;
+ Double_t fTRIG4_epCoinTime;
Double_t fROC1_eKCoinTime;
Double_t fROC2_eKCoinTime;
-
+ Double_t fTRIG1_eKCoinTime;
+ Double_t fTRIG4_eKCoinTime;
+
Double_t fROC1_ePiCoinTime;
Double_t fROC2_ePiCoinTime;
+ Double_t fTRIG1_ePiCoinTime;
+ Double_t fTRIG4_ePiCoinTime;
Double_t fROC1_ePosCoinTime; //electron-positron coin time
Double_t fROC2_ePosCoinTime;
+ Double_t fTRIG1_ePosCoinTime; //electron-positron coin time
+ Double_t fTRIG4_ePosCoinTime;
Double_t elec_coinCorr;
Double_t elecArm_BetaCalc;
@@ -128,10 +138,10 @@ public:
Double_t had_FPtime; //hadron focal plane time
// trigger times pTrig1 (SHMS 3/4 trig) and pTrig4 (HMS 3/4 trig)
- Int_t pTRIG1_TdcTime_ROC1;
- Int_t pTRIG4_TdcTime_ROC1;
- Int_t pTRIG1_TdcTime_ROC2;
- Int_t pTRIG4_TdcTime_ROC2;
+ Double_t pTRIG1_TdcTime_ROC1;
+ Double_t pTRIG4_TdcTime_ROC1;
+ Double_t pTRIG1_TdcTime_ROC2;
+ Double_t pTRIG4_TdcTime_ROC2;
//--------------------------------------------------------------------
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 2651bb7f0f37e4c6fc07b1683c5946281d903348..c09ce44dd2ed9d2e6191762c46e9df2d18086f71 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -27,6 +27,7 @@ the number of parameters per plane.
#include "TMath.h"
#include "TVectorD.h"
#include "THaApparatus.h"
+#include "THcHallCSpectrometer.h"
#include <cstring>
#include <cstdio>
@@ -446,6 +447,7 @@ Int_t THcDC::DefineVariables( EMode mode )
{ "chisq", "chisq/dof (golden track) ", "fChisq_best"},
{ "sp1_id", " (golden track) ", "fSp1_ID_best"},
{ "sp2_id", " (golden track) ", "fSp2_ID_best"},
+ { "InsideDipoleExit", " ","fInSideDipoleExit_best"},
{ "gtrack_nsp", " Number of space points in golden track ", "fNsp_best"},
{ "residual", "Residuals", "fResiduals"},
{ "residualExclPlane", "Residuals", "fResidualsExclPlane"},
@@ -540,6 +542,7 @@ void THcDC::ClearEvent()
fYp_fp_best=-10000.;
fChisq_best=kBig;
fNsp_best=0;
+ fInSideDipoleExit_best = kTRUE;
for(UInt_t i=0;i<fNChambers;i++) {
fChambers[i]->Clear();
}
@@ -681,6 +684,8 @@ void THcDC::SetFocalPlaneBestTrack(Int_t golden_track_index)
fY_fp_best=tr1->GetY();
fXp_fp_best=tr1->GetXP();
fYp_fp_best=tr1->GetYP();
+ THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+ fInSideDipoleExit_best = app->InsideDipoleExitWindow(fX_fp_best, fXp_fp_best ,fY_fp_best,fYp_fp_best);
fSp1_ID_best=tr1->GetSp1_ID();
fSp2_ID_best=tr1->GetSp2_ID();
fChisq_best=tr1->GetChisq();
@@ -771,11 +776,11 @@ void THcDC::LinkStubs()
std::vector<THcSpacePoint*> fSp;
fNSp=0;
fSp.clear();
- fSp.reserve(10);
+ fSp.reserve(100);
fNDCTracks=0; // Number of Focal Plane tracks found
fDCTracks->Delete();
// Make a vector of pointers to the SpacePoints
- if (fChambers[0]->GetNSpacePoints()+fChambers[1]->GetNSpacePoints()>10) return;
+ //if (fChambers[0]->GetNSpacePoints()+fChambers[1]->GetNSpacePoints()>10) return;
for(UInt_t ich=0;ich<fNChambers;ich++) {
Int_t nchamber=fChambers[ich]->GetChamberNum();
@@ -785,7 +790,8 @@ void THcDC::LinkStubs()
fSp[fNSp]->fNChamber = nchamber;
fSp[fNSp]->fNChamber_spnum = isp;
fNSp++;
- if (fNSp>10) break;
+ if (ich==0 && fNSp>50) break;
+ if (fNSp>100) break;
}
}
Double_t stubminx = 999999;
@@ -852,7 +858,7 @@ void THcDC::LinkStubs()
fStubTest = 1;
//stubtest=1; Used in h_track_tests.f
// Make a new track if there are not to many
- if(fNDCTracks < MAXTRACKS) {
+ if(fNDCTracks < fNTracksMaxFP) {
sptracks=0; // Number of tracks with this seed
stub_tracks[sptracks++] = fNDCTracks;
THcDCTrack *theDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
diff --git a/src/THcDC.h b/src/THcDC.h
index df5c89f35c733335a3646e82048b5d7d5ce1214a..1e88d8eb09793e23789c05fa68bd1ed6f280bc09 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -209,6 +209,7 @@ protected:
Double_t fChisq_best;
Int_t fSp1_ID_best;
Int_t fSp2_ID_best;
+ Bool_t fInSideDipoleExit_best;
// For accumulating statitics for efficiencies
Int_t fTotEvents;
Int_t* fNChamHits;
@@ -222,7 +223,7 @@ protected:
// double tan_angle, sin_angle, cos_angle;
// Intermediate structure for building
- static const UInt_t MAXTRACKS = 10;
+ static const UInt_t MAXTRACKS = 50;
public:
THcDriftChamberPlane* GetPlane(unsigned int i_plane) {
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 6f2dac9964fdd2e2362f6684e90c53980b6e3c0e..27e1440e410c950edc453021e55cf6cc8a6da0fe 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -1194,7 +1194,7 @@ void THcDriftChamber::LeftRight()
iswhit <<= 1;
}
}
- if (nplaneshit >= fNPlanes-1) {
+ if ( (nplaneshit >= fNPlanes-1) || (nplaneshit >= fNPlanes-2 && !fHMSStyleChambers)) {
Double_t chi2;
chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub);
if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl;
diff --git a/src/THcHallCSpectrometer.cxx b/src/THcHallCSpectrometer.cxx
index a1868918f58246e6fc10bc205499c8fbc3b531e1..2f84f3d435bb92a1bea5b01b4bbe946fe1fd5538 100644
--- a/src/THcHallCSpectrometer.cxx
+++ b/src/THcHallCSpectrometer.cxx
@@ -142,6 +142,7 @@ Int_t THcHallCSpectrometer::DefineVariables( EMode mode )
fIsSetup = ( mode == kDefine );
RVarDef vars[] = {
{ "tr.betachisq", "Chi2 of beta", "fTracks.THaTrack.GetBetaChi2()"},
+ { "tr.PruneSelect", "Prune Select ID", "fPruneSelect"},
{ "present", "Trigger Type includes this spectrometer", "fPresent"},
{ 0 }
};
@@ -277,17 +278,32 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
{"prune_chibeta", &fPruneChiBeta, kDouble, 0, 1},
{"prune_npmt", &fPruneNPMT, kDouble, 0, 1},
{"prune_fptime", &fPruneFpTime, kDouble, 0, 1},
+ {"prune_DipoleExit", &fPruneDipoleExit, kDouble, 0, 1},
{0}
};
// Default values
+ fPruneDipoleExit=0;
fSelUsingScin = 0;
fSelUsingPrune = 0;
+ fPruneXp = .2;
+ fPruneYp = .2;
+ fPruneYtar = 20.;
+ fPruneDelta = 30.;
+ fPruneBeta = 30.;
+ fPruneDf= 1;
+ fPruneChiBeta= 100.;
+ fPruneNPMT= 6;
+ fPruneFpTime= 1000.;
fPhi_lab = 0.;
fSatCorr=0.;
fMispointing_x=999.;
fMispointing_y=999.;
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ fUseHMSDipoleExitWindow=kFALSE;
+ fUseSHMSDipoleExitWindow=kFALSE;
+ if (prefix[0]=='h') fUseHMSDipoleExitWindow=kTRUE;
+ if (prefix[0]=='p') fUseSHMSDipoleExitWindow=kTRUE;
// mispointing in transport system y is horizontal and +x is vertical down
if (fMispointing_y == 999.) {
@@ -344,7 +360,7 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
//
- EnforcePruneLimits();
+ //EnforcePruneLimits();
#ifdef WITH_DEBUG
_logger->debug("Spectrometer {} ", GetName());
@@ -486,7 +502,7 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
TransportToLab(track->GetP(),track->GetTTheta(),track->GetTPhi(),pvect_temp);
track->SetPvect(pvect_temp);
}
-
+ fPruneSelect=-1.;
if (fHodo==0 || (( fSelUsingScin == 0 ) && ( fSelUsingPrune == 0 )) ) {
BestTrackSimple();
} else if (fHodo!=0 && fSelUsingPrune !=0) {
@@ -559,11 +575,15 @@ Int_t THcHallCSpectrometer::TrackCalc()
{
if( fNtracks > 0 ) {
Int_t hit_gold_track=0; // find track with index =0 which is best track
+ Int_t hit_dc_track=1; //
for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
- if (aTrack->GetIndex()==0) hit_gold_track=itrack;
+ if (aTrack->GetIndex()==0) {
+ hit_gold_track=itrack;
+ hit_dc_track = aTrack->GetTrkNum();
+ }
}
- fDC->SetFocalPlaneBestTrack(hit_gold_track);
+ fDC->SetFocalPlaneBestTrack(hit_dc_track-1);
fGoldenTrack = static_cast<THaTrack*>( fTracks->At(hit_gold_track) );
fTrkIfo = *fGoldenTrack;
fTrk = fGoldenTrack;
@@ -780,7 +800,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
testTracks[ptrack] = static_cast<THaTrack*>( fTracks->At(ptrack) );
if (!testTracks[ptrack]) return -1;
}
-
+ fPruneSelect = 0;
+ Double_t PruneSelect=0;
// ! Prune on xptar
nGood = 0;
for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
@@ -796,7 +817,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
-
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on yptar
nGood = 0;
for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
@@ -813,6 +835,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on ytar
nGood = 0;
@@ -829,6 +853,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on delta
nGood = 0;
@@ -845,6 +871,34 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
+
+ // ! Prune on dipole exit
+ nGood = 0;
+ for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
+ Double_t xfp=testTracks[ptrack]->GetX();
+ Double_t yfp=testTracks[ptrack]->GetY();
+ Double_t xpfp=testTracks[ptrack]->GetTheta();
+ Double_t ypfp=testTracks[ptrack]->GetPhi();
+ if ( fPruneDipoleExit==1 && InsideDipoleExitWindow(xfp,xpfp,yfp,ypfp) && ( keep[ptrack] ) ){
+ nGood ++;
+ }
+ }
+ if (nGood > 0 ) {
+ for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
+ Double_t xfp=testTracks[ptrack]->GetX();
+ Double_t yfp=testTracks[ptrack]->GetY();
+ Double_t xpfp=testTracks[ptrack]->GetTheta();
+ Double_t ypfp=testTracks[ptrack]->GetPhi();
+ if (!InsideDipoleExitWindow(xfp,xpfp,yfp,ypfp) ){
+ keep[ptrack] = kFALSE;
+ reject[ptrack] += 30;
+ }
+ }
+ }
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on beta
nGood = 0;
@@ -865,6 +919,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on deg. freedom for track chisq
nGood = 0;
@@ -882,6 +938,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
//! Prune on num pmt hits
nGood = 0;
for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
@@ -897,6 +955,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on beta chisqr
nGood = 0;
@@ -915,6 +975,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on fptime
nGood = 0;
@@ -932,6 +994,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on Y2 being hit
nGood = 0;
@@ -948,6 +1012,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Prune on X2 being hit
nGood = 0;
@@ -964,6 +1030,8 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
}
}
}
+ PruneSelect++;
+ if (nGood==1 && fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
// ! Pick track with best chisq if more than one track passed prune tests
Double_t chi2PerDeg = 0.;
@@ -976,7 +1044,9 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune()
chi2Min = chi2PerDeg;
}
}
- // Set index=0 for fGoodTrack
+ PruneSelect++;
+ if (fPruneSelect ==0 && fNtracks>1) fPruneSelect=PruneSelect;
+ // Set index=0 for fGoodTrack
for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
aTrack->SetIndex(1);
@@ -1042,6 +1112,45 @@ Bool_t THcHallCSpectrometer::IsMyEvent(Int_t evtype) const
return kFALSE;
}
+//
+Bool_t THcHallCSpectrometer::InsideDipoleExitWindow(Double_t x_fp, Double_t xp_fp, Double_t y_fp, Double_t yp_fp) {
+ Bool_t inside=kTRUE;
+ Double_t DipoleExitWindowZpos=0.;
+ if (fUseSHMSDipoleExitWindow) DipoleExitWindowZpos=-307.;
+ if (fUseHMSDipoleExitWindow) DipoleExitWindowZpos=-147.48;
+ Double_t xdip = x_fp + xp_fp*DipoleExitWindowZpos;
+ Double_t ydip = y_fp + yp_fp*DipoleExitWindowZpos;
+ if (fUseSHMSDipoleExitWindow) inside = SHMSDipoleExitWindow(xdip,ydip);
+ if (fUseHMSDipoleExitWindow) inside = HMSDipoleExitWindow(xdip,ydip);
+ return inside;
+}
+//
+Bool_t THcHallCSpectrometer::SHMSDipoleExitWindow(Double_t xdip,Double_t ydip ) {
+ Bool_t insideSHMS=kTRUE;
+ Double_t crad=23.81; // radius of semicircle
+ Double_t voffset= crad-24.035;
+ Double_t hwid=11.549/2.;
+ if ( TMath::Abs(ydip) < hwid) {
+ if (TMath::Abs(xdip) > (crad+voffset)) insideSHMS=kFALSE;
+ } else {
+ if ( ydip >=hwid) {
+ if ( ((xdip-voffset)*(xdip-voffset)+(ydip-hwid)*(ydip-hwid)) > crad*crad) insideSHMS=kFALSE;
+ }
+ if ( ydip <=-hwid) {
+ if ( ((xdip-voffset)*(xdip-voffset)+(ydip+hwid)*(ydip+hwid)) > crad*crad) insideSHMS=kFALSE;
+ }
+ }
+ return insideSHMS;
+ }
+//
+Bool_t THcHallCSpectrometer::HMSDipoleExitWindow(Double_t xdip,Double_t ydip) {
+ Bool_t insideHMS=kTRUE;
+ Double_t xpipe_offset = 2.8;
+ Double_t ypipe_offset = 0.0;
+ Double_t pipe_rad=46.507;
+ if ( ((xdip-xpipe_offset)*(xdip-xpipe_offset)+(ydip-ypipe_offset)*(ydip-ypipe_offset)) > pipe_rad*pipe_rad) insideHMS=kFALSE;
+ return insideHMS ;
+}
//_____________________________________________________________________________
ClassImp(THcHallCSpectrometer)
diff --git a/src/THcHallCSpectrometer.h b/src/THcHallCSpectrometer.h
index bb9a6a1cb4f5ec37947983235e783f0a4bce556b..a0a0061b71f36c1820aeeb319698802908b53b32 100644
--- a/src/THcHallCSpectrometer.h
+++ b/src/THcHallCSpectrometer.h
@@ -75,13 +75,19 @@ public:
virtual Int_t GetNumTypes() { return eventtypes.size(); };
virtual Bool_t IsPresent() {return fPresent;};
+ Bool_t InsideDipoleExitWindow(Double_t x_fp, Double_t xp_fp, Double_t y_fp, Double_t yp_fp);
protected:
void InitializeReconstruction();
- // Bool_t* fKeep;
+ Bool_t SHMSDipoleExitWindow(Double_t x_dip, Double_t y_dip);
+ Bool_t HMSDipoleExitWindow(Double_t x_dip, Double_t y_dip);
+ Bool_t fUseSHMSDipoleExitWindow;
+ Bool_t fUseHMSDipoleExitWindow;
+ // Bool_t* fKeep;
// Int_t* fReject;
+ Double_t fPruneDipoleExit;
Double_t fPartMass;
Double_t fPruneXp;
Double_t fPruneYp;
@@ -93,6 +99,7 @@ protected:
Double_t fPruneFpTime;
Double_t fPruneNPMT;
Double_t fSatCorr;
+ Double_t fPruneSelect;
Int_t fGoodTrack;
Int_t fSelUsingScin;
diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx
index 770148acd575abfe67996ac819bbd3f6e28750db..ae335b1cd918af23fceee1b5283ac78cfd32134d 100644
--- a/src/THcHelicity.cxx
+++ b/src/THcHelicity.cxx
@@ -8,27 +8,25 @@
// +1 = plus, -1 = minus, 0 = unknown
//
// Also supports in-time mode with delay = 0
-//
+//
////////////////////////////////////////////////////////////////////////
#include "THcHelicity.h"
+#include "TH1F.h"
#include "THaApparatus.h"
#include "THaEvData.h"
#include "THcGlobals.h"
#include "THcParmList.h"
-#include "TH1F.h"
#include "TMath.h"
#include <iostream>
using namespace std;
//_____________________________________________________________________________
-THcHelicity::THcHelicity( const char* name, const char* description,
- THaApparatus* app ):
- hcana::ConfigLogging<THaHelicityDet>( name, description, app ),
- fnQrt(-1), fHelDelay(8), fMAXBIT(30)
-{
+THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus* app)
+ : hcana::ConfigLogging<THaHelicityDet>(name, description, app), fnQrt(-1), fHelDelay(8),
+ fMAXBIT(30) {
// for( Int_t i = 0; i < NHIST; ++i )
// fHisto[i] = 0;
// memset(fHbits, 0, sizeof(fHbits));
@@ -36,8 +34,7 @@ THcHelicity::THcHelicity( const char* name, const char* description,
//_____________________________________________________________________________
THcHelicity::THcHelicity()
- : hcana::ConfigLogging<THaHelicityDet>(),fnQrt(-1), fHelDelay(8), fMAXBIT(30)
-{
+ : hcana::ConfigLogging<THaHelicityDet>(), fnQrt(-1), fHelDelay(8), fMAXBIT(30) {
// Default constructor for ROOT I/O
// for( Int_t i = 0; i < NHIST; ++i )
@@ -45,9 +42,8 @@ THcHelicity::THcHelicity()
}
//_____________________________________________________________________________
-THcHelicity::~THcHelicity()
-{
- DefineVariables( kDelete );
+THcHelicity::~THcHelicity() {
+ DefineVariables(kDelete);
// for( Int_t i = 0; i < NHIST; ++i ) {
// delete fHisto[i];
@@ -59,11 +55,11 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
// Call `Setup` before everything else.
Setup(GetName(), GetTitle());
- fFirstEvProcessed = kFALSE;
- fActualHelicity = kUnknown;
+ fFirstEvProcessed = kFALSE;
+ fActualHelicity = kUnknown;
fPredictedHelicity = kUnknown;
- fLastMPSTime = 0;
- fFoundMPS = kFALSE;
+ fLastMPSTime = 0;
+ fFoundMPS = kFALSE;
// Call initializer for base class.
// This also calls `ReadDatabase` and `DefineVariables`.
@@ -73,6 +69,9 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
return fStatus;
}
+ fEvNumCheck = 0;
+ fDisabled = kFALSE;
+
fStatus = kOK;
return fStatus;
}
@@ -86,131 +85,125 @@ void THcHelicity::Setup(const char* name, const char* description) {
}
//_____________________________________________________________________________
-Int_t THcHelicity::ReadDatabase( const TDatime& date )
-{
+Int_t THcHelicity::ReadDatabase(const TDatime& date) {
_logger->info("In THcHelicity::ReadDatabase");
- //cout << "In THcHelicity::ReadDatabase" << endl;
+ // cout << "In THcHelicity::ReadDatabase" << endl;
// Read general HelicityDet database values (e.g. fSign)
// Int_t st = THaHelicityDet::ReadDatabase( date );
// if( st != kOK )
// return st;
// Read readout parameters (ROC addresses etc.)
- Int_t st = THcHelicityReader::ReadDatabase( GetDBFileName(), GetPrefix(),
- date, fQWEAKDebug );
- if( st != kOK )
- return st;
+ Int_t st = THcHelicityReader::ReadDatabase(GetDBFileName(), GetPrefix(), date, fQWEAKDebug);
+ if (st != kOK)
+ return st;
- fSign = 1; // Default helicity sign
+ fSign = 1; // Default helicity sign
fRingSeed_reported_initial = 0; // Initial see that should predict reported
- // helicity of first quartet.
- fFirstCycle = -1; // First Cycle that starts a quad (0 to 3)
- fFreq = 29.5596;
- fHelDelay=8;
-
- DBRequest list[]= {
- // {"_hsign", &fSign, kInt, 0, 1},
- {"helicity_delay", &fHelDelay, kInt, 0, 1},
- {"helicity_freq", &fFreq, kDouble, 0, 1},
- // {"helicity_seed", &fRingSeed_reported_initial, kInt, 0, 1},
- // {"helicity_cycle", &fFirstCycle, kInt, 0, 1},
- {0}
- };
+ // helicity of first quartet.
+ fFirstCycle = -1; // First Cycle that starts a quad (0 to 3)
+ fFreq = 29.5596;
+ fHelDelay = 8;
+
+ DBRequest list[] = {// {"_hsign", &fSign, kInt, 0, 1},
+ {"helicity_delay", &fHelDelay, kInt, 0, 1},
+ {"helicity_freq", &fFreq, kDouble, 0, 1},
+ // {"helicity_seed", &fRingSeed_reported_initial, kInt, 0, 1},
+ // {"helicity_cycle", &fFirstCycle, kInt, 0, 1},
+ {0}};
gHcParms->LoadParmValues(list, "");
- fMAXBIT=30;
+ fMAXBIT = 30;
- fTIPeriod = 250000000.0/fFreq;
+ fTIPeriod = 250000000.0 / fFreq;
// maximum of event in the pattern, for now we are working with quartets
- // Int_t localpattern[4]={1,-1,-1,1};
+ // Int_t localpattern[4]={1,-1,-1,1};
// careful, the first value here should always +1
// for(int i=0;i<fQWEAKNPattern;i++)
// {
// fPatternSequence.push_back(localpattern[i]);
// }
- HWPIN=kTRUE;
+ HWPIN = kTRUE;
- fQuartet[0]=fQuartet[1]=fQuartet[2]=fQuartet[3]=0;
+ fQuartet[0] = fQuartet[1] = fQuartet[2] = fQuartet[3] = 0;
- if (fFirstCycle>=0 && fRingSeed_reported_initial!=0) {
+ if (fFirstCycle >= 0 && fRingSeed_reported_initial != 0) {
// Set the seed for predicted reported and predicted actual
} else {
// Initialize mode to find quartets and then seed
}
+ _logger->info(
+ "Helicity decoder initialized with frequency of {} Hz and reporting delay of {} cycles.",
+ fFreq, fHelDelay);
+ // cout << "Helicity decoder initialized with frequency of " << fFreq
+ // << " Hz and reporting delay of " << fHelDelay << " cycles." << endl;
+
return kOK;
}
//_____________________________________________________________________________
-void THcHelicity::MakePrefix()
-{
- THaDetector::MakePrefix();
-}
+void THcHelicity::MakePrefix() { THaDetector::MakePrefix(); }
//_____________________________________________________________________________
-Int_t THcHelicity::DefineVariables( EMode mode )
-{
+Int_t THcHelicity::DefineVariables(EMode mode) {
// Initialize global variables
_logger->info("Called THcHelicity::DefineVariables with mode == {}", mode);
- //cout << "Called THcHelicity::DefineVariables with mode == "
+ // cout << "Called THcHelicity::DefineVariables with mode == "
// << mode << endl;
- if( mode == kDefine && fIsSetup ) return kOK;
- fIsSetup = ( mode == kDefine );
+ if (mode == kDefine && fIsSetup)
+ return kOK;
+ fIsSetup = (mode == kDefine);
// Define standard variables from base class
- THaHelicityDet::DefineVariables( mode );
-
- const RVarDef var[] = {
- { "nqrt", "position of cycle in quartet", "fnQrt" },
- { "hel", "actual helicity for event", "fActualHelicity" },
- { "helrep", "reported helicity for event", "fReportedHelicity" },
- { "helpred", "predicted reported helicity for event", "fPredictedHelicity" },
- { "mps", "In MPS blanking period", "fMPS"},
- { 0 }
- };
- //cout << "Calling THcHelicity DefineVarsFromList" << endl;
+ THaHelicityDet::DefineVariables(mode);
+
+ const RVarDef var[] = {{"nqrt", "position of cycle in quartet", "fnQrt"},
+ {"hel", "actual helicity for event", "fActualHelicity"},
+ {"helrep", "reported helicity for event", "fReportedHelicity"},
+ {"helpred", "predicted reported helicity for event", "fPredictedHelicity"},
+ {"mps", "In MPS blanking period", "fMPS"},
+ {0}};
+ // cout << "Calling THcHelicity DefineVarsFromList" << endl;
_logger->info("Calling THcHelicity DefineVarsFromList");
- return DefineVarsFromList( var, mode );
+ return DefineVarsFromList(var, mode);
}
//_____________________________________________________________________________
-void THcHelicity::PrintEvent(Int_t evtnum)
-{
+void THcHelicity::PrintEvent(Int_t evtnum) {
- cout<<" ++++++ THcHelicity::Print ++++++\n";
+ cout << " ++++++ THcHelicity::Print ++++++\n";
- cout<<" +++++++++++++++++++++++++++++++++++++\n";
+ cout << " +++++++++++++++++++++++++++++++++++++\n";
return;
}
//_____________________________________________________________________________
-Int_t THcHelicity::Begin( THaRunBase* )
-{
+Int_t THcHelicity::Begin(THaRunBase*) {
THcHelicityReader::Begin();
// fHisto[0] = new TH1F("hel.seed","hel.seed",32,-1.5,30.5);
// fHisto[1] = new TH1F("hel.error.code","hel.error.code",35,-1.5,33.5);
-
+
return 0;
}
//_____________________________________________________________________________
-//void THcHelicity::FillHisto()
+// void THcHelicity::FillHisto()
//{
// fHisto[0]->Fill(fRing_NSeed);
// fHisto[1]->Fill(fErrorCode);
// return;
//}
//_____________________________________________________________________________
-void THcHelicity::SetErrorCode(Int_t error)
-{
- // used as a control for the helciity computation
+void THcHelicity::SetErrorCode(Int_t error) {
+ // used as a control for the helciity computation
// 2^0: if the reported number of events in a pattern is larger than fQWEAKNPattern
// 2^1: if the offset between the ring reported value and TIR value is not fOffsetTIRvsRing
// 2^2: if the reported time in the ring is 0
@@ -218,184 +211,238 @@ void THcHelicity::SetErrorCode(Int_t error)
// 2^4: if the helicity cannot be computed using the SetHelicity routine
// 2^5: if seed is being gathered
- if(fErrorCode==0)
- fErrorCode=(1<<error);
+ if (fErrorCode == 0)
+ fErrorCode = (1 << error);
// only one reported error at the time
return;
}
//_____________________________________________________________________________
-void THcHelicity::Clear( Option_t* opt )
-{
+void THcHelicity::Clear(Option_t* opt) {
// Clear event-by-event data
THaHelicityDet::Clear(opt);
THcHelicityReader::Clear(opt);
fEvtype = 0;
- fQrt=0;
- fErrorCode=0;
+ fQrt = 0;
+ fErrorCode = 0;
return;
}
//_____________________________________________________________________________
-Int_t THcHelicity::Decode( const THaEvData& evdata )
-{
+Int_t THcHelicity::Decode(const THaEvData& evdata) {
// Decode Helicity data.
// Return 1 if helicity was assigned, 0 if not, <0 if error.
- Int_t err = ReadData( evdata ); // from THcHelicityReader class
- if( err ) {
- Error( Here("THcHelicity::Decode"), "Error decoding helicity data." );
+ Int_t err = ReadData(evdata); // from THcHelicityReader class
+ if (err) {
+ Error(Here("THcHelicity::Decode"), "Error decoding helicity data.");
return err;
}
- fReportedHelicity = (fIsHelp?(fIsHelm?kUnknown:kPlus):(fIsHelm?kMinus:kUnknown));
- fMPS = fIsMPS?1:0;
- if(fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
- fActualHelicity = fIsMPS?kUnknown:fReportedHelicity;
+ fReportedHelicity = (fIsHelp ? (fIsHelm ? kUnknown : kPlus) : (fIsHelm ? kMinus : kUnknown));
+ fMPS = fIsMPS ? 1 : 0;
+ if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
+ fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
+ return 0;
+ }
+
+ if (fDisabled) {
+ fActualHelicity = kUnknown;
+ return 0;
+ }
+
+ fEvNumCheck++;
+ Int_t evnum = evdata.GetEvNum();
+ if (fEvNumCheck != evnum) {
+ _logger->info("THcHelicity: Missed {} events at event {}.", evnum - fEvNumCheck, evnum);
+ _logger->info(" Disabling helicity decoding for rest of run.");
+ _logger->info(
+ " Make sure \"RawDecode_master in cuts file accepts all physics events.");
+ fDisabled = kTRUE;
+ fActualHelicity = kUnknown;
return 0;
}
- Long64_t lastlastmpstime = fLastMPSTime;
- if(fFirstEvProcessed) { // Normal processing
+
+ fActualHelicity = -10.0;
+ if (fFirstEvProcessed) { // Normal processing
+ // cout << evnum << " " << fNCycle << " " << fIsMPS << " " << fFoundMPS << " " << fTITime <<
+ // " "
+ // << fLastMPSTime << " " << fNBits << endl;
Int_t missed = 0;
// Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0;
- if(fIsMPS) {
- if(fFoundMPS) {
- missed = TMath::Nint(fTITime/fTIPeriod-fLastMPSTime/fTIPeriod);
- if(missed < 1) { // was <=1
- fLastMPSTime = (fTITime+fLastMPSTime+missed*fTIPeriod)/2;
- fIsNewCycle = kTRUE;
- fActualHelicity = kUnknown;
- fPredictedHelicity = kUnknown;
- } else {
- fLastMPSTime = (fLastMPSTime + fTITime - missed*fTIPeriod)/2;
- }
- // If there is a skip, pass it off to next non MPS event
- // Need to also check here for missed MPS's
- // cout << "Found MPS" << endl;
- // check for Nint((time-last)/period) > 1
+ if (fIsMPS) {
+ fActualHelicity = kUnknown;
+ fPredictedHelicity = kUnknown;
+ if (fFoundMPS) {
+ missed = TMath::Nint(fTITime / fTIPeriod - fLastMPSTime / fTIPeriod);
+ if (missed < 1) { // was <=1
+ fLastMPSTime = (fTITime + fLastMPSTime + missed * fTIPeriod) / 2;
+ fIsNewCycle = kTRUE;
+ fActualHelicity = kUnknown;
+ fPredictedHelicity = kUnknown;
+ } else {
+ fLastMPSTime = (fLastMPSTime + fTITime - missed * fTIPeriod) / 2;
+ }
+ // If there is a skip, pass it off to next non MPS event
+ // Need to also check here for missed MPS's
+ // cout << "Found MPS" << endl;
+ // check for Nint((time-last)/period) > 1
} else {
- fFoundMPS = kTRUE;
- fLastMPSTime = fTITime;
+ fFoundMPS = kTRUE;
+ fLastMPSTime = fTITime;
}
- } else if (fFoundMPS) { //
- if(fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods
- missed = TMath::Nint(floor((fTITime-fLastMPSTime)/fTIPeriod));
- if(missed > 1) {
- // cout << "Missed " << missed << " MPSes" << endl;
- Int_t newNCycle = fNCycle + missed -1; // How many cycles really missed
- Int_t quartets_missed = (newNCycle-fFirstCycle)/4 - (fNCycle-fFirstCycle)/4;
- for(Int_t i=0;i<quartets_missed;i++) { // Advance the seeds.
- fRingSeed_reported = RanBit30(fRingSeed_reported);
- fRingSeed_actual = RanBit30(fRingSeed_actual);
- }
- int quartetphase = (newNCycle-fFirstCycle)%4;
- // cout << " " << fNCycle << " " << newNCycle << " " << fFirstCycle << " " << quartets_missed << " " << quartetphase << endl;
- fQuartetStartHelicity = (fRingSeed_actual&1)?kPlus:kMinus;
- fQuartetStartPredictedHelicity = (fRingSeed_reported&1)?kPlus:kMinus;
- fActualHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartHelicity:-fQuartetStartHelicity;
- fPredictedHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartPredictedHelicity:-fQuartetStartPredictedHelicity;
- // cout << "Cycles " << fNCycle << " " << newNCycle << " " << fFirstCycle
- // << " skipped " << quartets_missed << " quartets" << endl;
- fNCycle = newNCycle;
- // Need to reset fQuartet to reflect where we are based on the current
- // reported helicity. So we don't fail quartet testing.
- // But only do this if we are calibrated.
- if(fNBits >= fMAXBIT) {
- if (((fNCycle - fFirstCycle)%2)==1) {
- fQuartet[0] = fReportedHelicity;
- fQuartet[1] = fQuartet[2] = -fQuartet[0];
- } else {
- fQuartet[0] = fQuartet[1] = -fReportedHelicity;
- fQuartet[2] = -fQuartet[1];
- }
- } else {
- fQuartet[0] = fReportedHelicity;
- fQuartet[1] = 0;
- }
- }
- fLastMPSTime += missed*fTIPeriod;
- fIsNewCycle = kTRUE;
- fLastReportedHelicity = fReportedHelicity;
+ } else if (fFoundMPS) { //
+ if (fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods
+ missed = TMath::Nint(floor((fTITime - fLastMPSTime) / fTIPeriod));
+ if (missed > 1) {
+ // cout << "Missed " << missed << " MPSes" << endl;
+ Int_t newNCycle = fNCycle + missed - 1; // How many cycles really missed
+ Int_t quartets_missed = (newNCycle - fFirstCycle) / 4 - (fNCycle - fFirstCycle) / 4;
+ for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
+ fRingSeed_reported = RanBit30(fRingSeed_reported);
+ fRingSeed_actual = RanBit30(fRingSeed_actual);
+ }
+ int quartetphase = (newNCycle - fFirstCycle) % 4;
+ // cout << " " << fNCycle << " " << newNCycle << " " << fFirstCycle << " " <<
+ // quartets_missed << " " << quartetphase << endl; cout << "Cycles " << fNCycle << "
+ // " << newNCycle << " " << fFirstCycle
+ // << " skipped " << quartets_missed << " quartets" << endl;
+ fNCycle = newNCycle;
+ // Need to reset fQuartet to reflect where we are based on the current
+ // reported helicity. So we don't fail quartet testing.
+ // But only do this if we are calibrated.
+ if (fNBits >= fMAXBIT) {
+ fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
+ fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
+ fActualHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity
+ : -fQuartetStartHelicity;
+ fPredictedHelicity = (quartetphase == 0 || quartetphase == 3)
+ ? fQuartetStartPredictedHelicity
+ : -fQuartetStartPredictedHelicity;
+
+ if (((fNCycle - fFirstCycle) % 2) == 1) {
+ fQuartet[0] = fReportedHelicity;
+ fQuartet[1] = fQuartet[2] = -fQuartet[0];
+ } else {
+ fQuartet[0] = fQuartet[1] = -fReportedHelicity;
+ fQuartet[2] = -fQuartet[1];
+ }
+ } else {
+ fActualHelicity = kUnknown;
+ fQuartet[0] = fReportedHelicity;
+ fQuartet[1] = 0;
+ }
+ }
+ fLastMPSTime += missed * fTIPeriod;
+ fIsNewCycle = kTRUE;
+ fLastReportedHelicity = fReportedHelicity;
+ } else { // No missed periods. Get helicities from rings
+ if (fNBits >= fMAXBIT) {
+ int quartetphase = (fNCycle - fFirstCycle) % 4;
+ fQuartetStartHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
+ fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
+ fActualHelicity = (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity
+ : -fQuartetStartHelicity;
+ fPredictedHelicity = (quartetphase == 0 || quartetphase == 3)
+ ? fQuartetStartPredictedHelicity
+ : -fQuartetStartPredictedHelicity;
+ } else {
+ fActualHelicity = 0;
+ }
}
- if(fIsNewCycle) {
- fQuartet[3]=fQuartet[2]; fQuartet[2]=fQuartet[1]; fQuartet[1]=fQuartet[0];
- fQuartet[0]=fReportedHelicity;
- fNCycle++;
- if((fNCycle-fFirstCycle)%4 == 3) {// Test if last in a quartet
- if((abs(fQuartet[0]+fQuartet[3]-fQuartet[1]-fQuartet[2])==4)) {
- if(!fFoundQuartet) {
- // fFirstCycle = fNCycle - 3;
- _logger->info("Quartet potentially found, starting at cycle {} - event {}", fFirstCycle, evdata.GetEvNum());
- //cout << "Quartet potentially found, starting at cycle " << fFirstCycle << " - event "
+ if (fIsNewCycle) {
+ fQuartet[3] = fQuartet[2];
+ fQuartet[2] = fQuartet[1];
+ fQuartet[1] = fQuartet[0];
+ fQuartet[0] = fReportedHelicity;
+ fNCycle++;
+ if ((fNCycle - fFirstCycle) % 4 == 3) { // Test if last in a quartet
+ if ((abs(fQuartet[0] + fQuartet[3] - fQuartet[1] - fQuartet[2]) == 4)) {
+ if (!fFoundQuartet) {
+ // fFirstCycle = fNCycle - 3;
+ _logger->info("Quartet potentially found, starting at cycle {} - event {}",
+ fFirstCycle, evdata.GetEvNum());
+ // cout << "Quartet potentially found, starting at cycle " << fFirstCycle << " - event
+ // "
// << evdata.GetEvNum() << endl;
fFoundQuartet = kTRUE;
- }
- } else {
- if(fNCycle - fFirstCycle > 4) { // Not at start of run. Reset
- _logger->warn("Lost quartet sync at cycle {} - event {}", fNCycle,evdata.GetEvNum());
- _logger->warn("{} {} {} {}",fQuartet[0],fQuartet[1],fQuartet[2],fQuartet[3]);
- //cout << "Lost quartet sync at cycle " << fNCycle << " - event " << evdata.GetEvNum()
+ }
+ } else {
+ if (fNCycle - fFirstCycle > 4) { // Not at start of run. Reset
+ _logger->warn("Lost quartet sync at cycle {} - event {}", fNCycle, evdata.GetEvNum());
+ _logger->warn("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
+ // cout << "Lost quartet sync at cycle " << fNCycle << " - event " <<
+ // evdata.GetEvNum()
// << endl;
- //cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " << fQuartet[3]
+ // cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " <<
+ // fQuartet[3]
// << endl;
- fFirstCycle += 4*((fNCycle-fFirstCycle)/4); // Update, but don't change phase
- }
- fFoundQuartet = kFALSE;
+ fFirstCycle += 4 * ((fNCycle - fFirstCycle) / 4); // Update, but don't change phase
+ }
+ fFoundQuartet = kFALSE;
fNBits = 0;
- _logger->info("Searching for first of a quartet at cycle {} - event {}", fFirstCycle, evdata.GetEvNum());
- //cout << "Searching for first of a quartet at cycle "
+ _logger->info("Searching for first of a quartet at cycle {} - event {}", fFirstCycle,
+ evdata.GetEvNum());
+ // cout << "Searching for first of a quartet at cycle "
// << " " << fFirstCycle << " - event " << evdata.GetEvNum() << endl;
- //cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " << fQuartet[3]
+ // cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " " << fQuartet[3]
// << endl;
_logger->info("{} {} {} {}", fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
- fFirstCycle++;
- }
- }
- // Load the actual helicity. Calibrate if not calibrated.
- fActualHelicity = kUnknown;
- LoadHelicity(fReportedHelicity, fNCycle, missed);
- fLastReportedHelicity = fReportedHelicity;
- fIsNewCycle = kFALSE;
- // cout << fTITime/250000000.0 << " " << fNCycle << " " << fReportedHelicity << endl;
- // cout << fNCycle << ": " << fReportedHelicity << " "
- // << fPredictedHelicity << " " << fActualHelicity << endl;
+ fFirstCycle++;
+ }
+ }
+ // Load the actual helicity. Calibrate if not calibrated.
+ fActualHelicity = kUnknown;
+ LoadHelicity(fReportedHelicity, fNCycle, missed);
+ fLastReportedHelicity = fReportedHelicity;
+ fIsNewCycle = kFALSE;
+ // cout << fTITime/250000000.0 << " " << fNCycle << " " << fReportedHelicity << endl;
+ // cout << fNCycle << ": " << fReportedHelicity << " "
+ // << fPredictedHelicity << " " << fActualHelicity << endl;
}
// Ignore until a MPS Is found
+ } else { // No MPS found yet
+ fActualHelicity = kUnknown;
}
} else {
- //cout << "Initializing" << endl;
+ // cout << "Initializing" << endl;
_logger->info("Initializing Helicity");
fLastReportedHelicity = fReportedHelicity;
- fActualHelicity = kUnknown;
- fPredictedHelicity = kUnknown;
- fFirstEvTime = fTITime;
- fLastEvTime = fTITime;
- fLastMPSTime = fTITime; // Not necessarily during the MPS
- fNCycle = 0;
- fFirstEvProcessed = kTRUE;
- fFoundMPS = kFALSE;
- fFoundQuartet = kFALSE;
- fIsNewCycle = kFALSE;
- fNBits = 0;
+ fActualHelicity = kUnknown;
+ fPredictedHelicity = kUnknown;
+ fFirstEvTime = fTITime;
+ fLastEvTime = fTITime;
+ fLastMPSTime = fTITime; // Not necessarily during the MPS
+ fNCycle = 0;
+ fFirstEvProcessed = kTRUE;
+ fFoundMPS = kFALSE;
+ fFoundQuartet = kFALSE;
+ fIsNewCycle = kFALSE;
+ fNBits = 0;
}
- // cout << setprecision(9) << "HEL " << fTITime/250000000.0 << " " << fNCycle << "(" << (fNCycle-fFirstCycle)%4 << "): "
- // << fMPS << " " << fReportedHelicity << " "
- // << fPredictedHelicity << " " << fActualHelicity << " " << lastlastmpstime/250000000.0 << endl;
- // bitset<32>(v)
+
+ // Some sanity checks
+ if (fActualHelicity < -5) {
+ _logger->info("Actual Helicity never got defined");
+ }
+ if (fNBits < fMAXBIT) {
+ if (fActualHelicity == -1 || fActualHelicity == 1) {
+ cout << "Helicity of " << fActualHelicity << " reported prematurely at cycle " << fNCycle
+ << endl;
+ }
+ }
+ fLastActualHelicity = fActualHelicity;
return 0;
}
//_____________________________________________________________________________
-Int_t THcHelicity::End( THaRunBase* )
-{
+Int_t THcHelicity::End(THaRunBase*) {
// End of run processing. Write histograms.
THcHelicityReader::End();
@@ -406,158 +453,160 @@ Int_t THcHelicity::End( THaRunBase* )
}
//_____________________________________________________________________________
-void THcHelicity::SetDebug( Int_t level )
-{
- // Set debug level of this detector as well as the THcHelicityReader
+void THcHelicity::SetDebug(Int_t level) {
+ // Set debug level of this detector as well as the THcHelicityReader
// helper class.
- THaHelicityDet::SetDebug( level );
+ THaHelicityDet::SetDebug(level);
fQWEAKDebug = level;
}
//_____________________________________________________________________________
-void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles)
-{
+void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles) {
// static const char* const here = "THcHelicity::LoadHelicity";
- int quartetphase = (cyclecount-fFirstCycle)%4;
- fnQrt = quartetphase;
+ int quartetphase = (cyclecount - fFirstCycle) % 4;
+ fnQrt = quartetphase;
- if(missedcycles > 1) { // If we missed windows
- if(fNBits< fMAXBIT) { // and we haven't gotten the seed, start over
+ if (missedcycles > 1) { // If we missed windows
+ if (fNBits < fMAXBIT) { // and we haven't gotten the seed, start over
fNBits = 0;
return;
}
}
- if(!fFoundQuartet) { // Wait until we have found quad phase before starting
- return; // to calibrate
+ if (!fFoundQuartet) { // Wait until we have found quad phase before starting
+ return; // to calibrate
}
- if(quartetphase == 0) { // Start of a quad
- if(fNBits < fMAXBIT) {
- if(fNBits == 0) {
- _logger->info("Start calibrating at cycle {}" ,cyclecount );
- //cout << "Start calibrating at cycle " << cyclecount << endl;
- fRingSeed_reported = 0;
+ if (quartetphase == 0) { // Start of a quad
+ if (fNBits < fMAXBIT) {
+ if (fNBits == 0) {
+ _logger->info("Start calibrating at cycle {}", cyclecount);
+ // cout << "Start calibrating at cycle " << cyclecount << endl;
+ fRingSeed_reported = 0;
}
- if(fReportedHelicity == kPlus) {
- fRingSeed_reported = ((fRingSeed_reported<<1) | 1) & 0x3FFFFFFF;
+ if (fReportedHelicity == kPlus) {
+ fRingSeed_reported = ((fRingSeed_reported << 1) | 1) & 0x3FFFFFFF;
} else {
- fRingSeed_reported = (fRingSeed_reported<<1) & 0x3FFFFFFF;
+ fRingSeed_reported = (fRingSeed_reported << 1) & 0x3FFFFFFF;
}
fNBits++;
- if(fReportedHelicity == kUnknown) {
- fNBits = 0;
- fRingSeed_reported = 0;
- } else if (fNBits==fMAXBIT) {
-
- _logger->info("Seed Found {} at cycle {} with first cycle {}" , fRingSeed_reported , cyclecount , fFirstCycle );
- //cout << "Seed Found " << hex << fRingSeed_reported << dec << " at cycle " << cyclecount << " with first cycle " << fFirstCycle << endl;
- Int_t backseed = GetSeed30(fRingSeed_reported);
- _logger->info("Seed at cycle {} should be {}", fFirstCycle , backseed);
- //cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec << endl;
+ if (fReportedHelicity == kUnknown) {
+ fNBits = 0;
+ fRingSeed_reported = 0;
+ } else if (fNBits == fMAXBIT) {
+
+ _logger->info("Seed Found {} at cycle {} with first cycle {}", fRingSeed_reported,
+ cyclecount, fFirstCycle);
+ // cout << "Seed Found " << hex << fRingSeed_reported << dec << " at cycle " << cyclecount
+ // << " with first cycle " << fFirstCycle << endl;
+ Int_t backseed = GetSeed30(fRingSeed_reported);
+ _logger->info("Seed at cycle {} should be {}", fFirstCycle, backseed);
+ // cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec <<
+ // endl;
}
fActualHelicity = kUnknown;
} else if (fNBits >= fMAXBIT) {
fRingSeed_reported = RanBit30(fRingSeed_reported);
- if(fNBits==fMAXBIT) {
- fRingSeed_actual = fRingSeed_reported;
- for(Int_t i=0;i<fHelDelay/4; i++) {
- fRingSeed_actual = RanBit30(fRingSeed_actual);
- }
- fNBits++;
+ if (fNBits == fMAXBIT) {
+ fRingSeed_actual = fRingSeed_reported;
+ for (Int_t i = 0; i < fHelDelay / 4; i++) {
+ fRingSeed_actual = RanBit30(fRingSeed_actual);
+ }
+ fNBits++;
} else {
- fRingSeed_actual = RanBit30(fRingSeed_actual);
+ fRingSeed_actual = RanBit30(fRingSeed_actual);
}
- fActualHelicity = (fRingSeed_actual&1)?kPlus:kMinus;
- fPredictedHelicity = (fRingSeed_reported&1)?kPlus:kMinus;
- // if(fTITime/250000000.0 > 380.0) cout << fTITime/250000000.0 << " " << fNCycle << " " << hex <<
- // fRingSeed_reported << " " << fRingSeed_actual << dec << endl;
- if(fReportedHelicity != fPredictedHelicity) {
- _logger->warn("Helicity prediction failed {} {} {}", fReportedHelicity, fPredictedHelicity, fActualHelicity);
- //cout << "Helicity prediction failed " << fReportedHelicity << " "
- // << fPredictedHelicity << " " << fActualHelicity << endl;
- //cout << hex << fRingSeed_reported << " " << fRingSeed_actual << dec << endl;
- fNBits = 0; // Need to reaquire seed
- fActualHelicity = kUnknown;
- fPredictedHelicity = kUnknown;
+ fActualHelicity = (fRingSeed_actual & 1) ? kPlus : kMinus;
+ fPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
+ // if(fTITime/250000000.0 > 380.0) cout << fTITime/250000000.0 << " " << fNCycle << " "
+ // << hex <<
+ // fRingSeed_reported << " " << fRingSeed_actual << dec
+ //<< endl;
+ if (fReportedHelicity != fPredictedHelicity) {
+ _logger->warn("Helicity prediction failed {} {} {}", fReportedHelicity, fPredictedHelicity,
+ fActualHelicity);
+ // cout << "Helicity prediction failed " << fReportedHelicity << " "
+ // << fPredictedHelicity << " " << fActualHelicity << endl;
+ // cout << hex << fRingSeed_reported << " " << fRingSeed_actual << dec << endl;
+ fNBits = 0; // Need to reaquire seed
+ fActualHelicity = kUnknown;
+ fPredictedHelicity = kUnknown;
}
}
- fQuartetStartHelicity = fActualHelicity;
+ fQuartetStartHelicity = fActualHelicity;
fQuartetStartPredictedHelicity = fPredictedHelicity;
- } else { // Not the beginning of a quad
- if(fNBits>=fMAXBIT) {
- fActualHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartHelicity:-fQuartetStartHelicity;
- fPredictedHelicity = (quartetphase==0||quartetphase==3)?
- fQuartetStartPredictedHelicity:-fQuartetStartPredictedHelicity;
+ } else { // Not the beginning of a quad
+ if (fNBits >= fMAXBIT) {
+ fActualHelicity =
+ (quartetphase == 0 || quartetphase == 3) ? fQuartetStartHelicity : -fQuartetStartHelicity;
+ fPredictedHelicity = (quartetphase == 0 || quartetphase == 3)
+ ? fQuartetStartPredictedHelicity
+ : -fQuartetStartPredictedHelicity;
}
}
return;
}
//_____________________________________________________________________________
-Int_t THcHelicity::RanBit30(Int_t ranseed)
-{
-
- UInt_t bit7 = (ranseed & 0x00000040) != 0;
- UInt_t bit28 = (ranseed & 0x08000000) != 0;
- UInt_t bit29 = (ranseed & 0x10000000) != 0;
- UInt_t bit30 = (ranseed & 0x20000000) != 0;
+Int_t THcHelicity::RanBit30(Int_t ranseed) {
+
+ UInt_t bit7 = (ranseed & 0x00000040) != 0;
+ UInt_t bit28 = (ranseed & 0x08000000) != 0;
+ UInt_t bit29 = (ranseed & 0x10000000) != 0;
+ UInt_t bit30 = (ranseed & 0x20000000) != 0;
UInt_t newbit = (bit30 ^ bit29 ^ bit28 ^ bit7) & 0x1;
-
- if(ranseed<=0) {
- if(fQWEAKDebug>1)
- std::cerr<<"ranseed must be greater than zero!"<<"\n";
+ if (ranseed <= 0) {
+ if (fQWEAKDebug > 1)
+ std::cerr << "ranseed must be greater than zero!"
+ << "\n";
newbit = 0;
}
- ranseed = ( (ranseed<<1) | newbit ) & 0x3FFFFFFF;
- //here ranseed is changed
- if(fQWEAKDebug>1)
- {
- cout<< "THcHelicity::RanBit30, newbit="<<newbit<<"\n";
- }
+ ranseed = ((ranseed << 1) | newbit) & 0x3FFFFFFF;
+ // here ranseed is changed
+ if (fQWEAKDebug > 1) {
+ cout << "THcHelicity::RanBit30, newbit=" << newbit << "\n";
+ }
return ranseed;
-
}
//_____________________________________________________________________________
-Int_t THcHelicity::GetSeed30(Int_t currentseed)
+Int_t THcHelicity::GetSeed30(Int_t currentseed)
/* Back track the seed by 30 samples */
{
#if 1
Int_t seed = currentseed;
- for(Int_t i=0;i<30;i++) {
- UInt_t bit1 = (seed & 0x00000001) != 0;
- UInt_t bit8 = (seed & 0x00000080) != 0;
- UInt_t bit29 = (seed & 0x10000000) != 0;
- UInt_t bit30 = (seed & 0x20000000) != 0;
-
+ for (Int_t i = 0; i < 30; i++) {
+ UInt_t bit1 = (seed & 0x00000001) != 0;
+ UInt_t bit8 = (seed & 0x00000080) != 0;
+ UInt_t bit29 = (seed & 0x10000000) != 0;
+ UInt_t bit30 = (seed & 0x20000000) != 0;
+
UInt_t newbit30 = (bit30 ^ bit29 ^ bit8 ^ bit1) & 0x1;
- seed = (seed >> 1) | (newbit30<<29);
+ seed = (seed >> 1) | (newbit30 << 29);
}
#else
Int_t bits = currentseed;
- Int_t seed=0;
- for(Int_t i=0;i<30;i++) {
+ Int_t seed = 0;
+ for (Int_t i = 0; i < 30; i++) {
Int_t val;
// XOR at virtual position 0 and 29
- if(i==0) {
- val = ((bits & (1<<(i)))!=0) ^ ((bits & (1<<(i+29)))!=0);
+ if (i == 0) {
+ val = ((bits & (1 << (i))) != 0) ^ ((bits & (1 << (i + 29))) != 0);
} else {
- val = ((bits & (1<<(i)))!=0) ^ ((seed & (1<<(i-1)))!=0);
+ val = ((bits & (1 << (i))) != 0) ^ ((seed & (1 << (i - 1))) != 0);
}
- if(i<=1) {
- val = ((bits & (1<<(1-i)))!=0) ^ val;
+ if (i <= 1) {
+ val = ((bits & (1 << (1 - i))) != 0) ^ val;
} else {
- val = ((seed & (1<<(i-2)))!=0) ^ val;
+ val = ((seed & (1 << (i - 2))) != 0) ^ val;
}
- if(i<=22) {
- val = ((bits & (1<<(i-22)))!=0) ^ val;
+ if (i <= 22) {
+ val = ((bits & (1 << (i - 22))) != 0) ^ val;
} else {
- val = ((seed & (1<<(i-23)))!=0) ^ val;
+ val = ((seed & (1 << (i - 23))) != 0) ^ val;
}
- seed |= (val<<i);
+ seed |= (val << i);
}
#endif
return seed;
diff --git a/src/THcHelicity.h b/src/THcHelicity.h
index 8eabbffdd52624183e40397db7f04666140c3288..70c651e529ef26e0b829080c8226018f2e47b18d 100644
--- a/src/THcHelicity.h
+++ b/src/THcHelicity.h
@@ -95,6 +95,9 @@ protected:
Double_t fErrorCode;
Int_t fEvtype; // Current CODA event type
+ Int_t fLastActualHelicity;
+ Int_t fEvNumCheck;
+ Bool_t fDisabled;
static const Int_t NHIST = 2;
TH1F* fHisto[NHIST];
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index f1856c64372bc2df1ebf67cac57d991dfd94c2ee..45bc5d4dcb9b78eafbe87aa21963fb979e949694 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -9,29 +9,30 @@ hodoscope array, not just one plane.
*/
-#include "THcSignalHit.h"
-#include "THcShower.h"
+#include "TClass.h"
+#include "THaSubDetector.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcRawShowerHit.h"
-#include "TClass.h"
+#include "THcScintPlaneCluster.h"
+#include "THcShower.h"
+#include "THcSignalHit.h"
#include "math.h"
-#include "THaSubDetector.h"
-#include "THcHodoscope.h"
-#include "THaEvData.h"
+#include "TClonesArray.h"
+#include "THaCutList.h"
#include "THaDetMap.h"
-#include "THcDetectorMap.h"
+#include "THaEvData.h"
#include "THaGlobals.h"
-#include "THaCutList.h"
+#include "THaTrack.h"
+#include "THcDetectorMap.h"
#include "THcGlobals.h"
+#include "THcHodoscope.h"
#include "THcParmList.h"
+#include "TMath.h"
#include "VarDef.h"
#include "VarType.h"
-#include "THaTrack.h"
-#include "TClonesArray.h"
-#include "TMath.h"
#include "THaOutput.h"
#include "TTree.h"
@@ -40,13 +41,13 @@ hodoscope array, not just one plane.
#include <vector>
#include <algorithm>
-#include <cstring>
+#include <array>
#include <cstdio>
#include <cstdlib>
-#include <iostream>
-#include <iomanip>
+#include <cstring>
#include <fstream>
-#include <array>
+#include <iomanip>
+#include <iostream>
#include "hcana/helpers.hxx"
#include <cassert>
@@ -55,98 +56,93 @@ using namespace std;
using std::vector;
//_____________________________________________________________________________
-THcHodoscope::THcHodoscope( const char* name, const char* description,
- THaApparatus* apparatus ) :
- hcana::ConfigLogging<THaNonTrackingDetector>(name,description,apparatus)
-{
+THcHodoscope::THcHodoscope(const char* name, const char* description, THaApparatus* apparatus)
+ : hcana::ConfigLogging<THaNonTrackingDetector>(name, description, apparatus) {
// Constructor
- //fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+ // fTrackProj = new TClonesArray( "THaTrackProj", 5 );
// Construct the planes
- fNPlanes = 0; // No planes until we make them
- fStartTime=-1e5;
- fGoodStartTime=kFALSE;
+ fNPlanes = 0; // No planes until we make them
+ fStartTime = -1e5;
+ fGoodStartTime = kFALSE;
}
//_____________________________________________________________________________
-THcHodoscope::THcHodoscope( ) :
- hcana::ConfigLogging<THaNonTrackingDetector>()
-{
+THcHodoscope::THcHodoscope() : hcana::ConfigLogging<THaNonTrackingDetector>() {
// Constructor
}
//_____________________________________________________________________________
-void THcHodoscope::Setup(const char* name, const char* description)
-{
+void THcHodoscope::Setup(const char* name, const char* description) {
/**
Create the scintillator plane objects for the hodoscope.
-
+
Uses the Xhodo_num_planes and Xhodo_plane_names to get the number of
planes and their names.
Gets a pointer to the Cherenkov named "cer" ("hgcer" in the case of the SHMS.)
- */
- if( IsZombie()) return;
+ */
+ if (IsZombie())
+ return;
// fDebug = 1; // Keep this at one while we're working on the code
char prefix[2];
- prefix[0]=tolower(GetApparatus()->GetName()[0]);
- prefix[1]='\0';
+ prefix[0] = tolower(GetApparatus()->GetName()[0]);
+ prefix[1] = '\0';
TString temp(prefix[0]);
- fSHMS=kFALSE;
- if (temp == "p" ) fSHMS=kTRUE;
- TString histname=temp+"_timehist";
- hTime = new TH1F(histname,"",400,0,200);
+ fSHMS = kFALSE;
+ if (temp == "p")
+ fSHMS = kTRUE;
+ TString histname = temp + "_timehist";
+ hTime = new TH1F(histname, "", 400, 0, 200);
// cout << " fSHMS = " << fSHMS << endl;
- string planenamelist;
- DBRequest listextra[]={
- {"hodo_num_planes", &fNPlanes, kInt},
- {"hodo_plane_names",&planenamelist, kString},
- {"hodo_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
- {"hodo_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
- {0}
- };
- //fNPlanes = 4; // Default if not defined
- fTDC_RefTimeCut = 0; // Minimum allowed reference times
+ string planenamelist;
+ DBRequest listextra[] = {{"hodo_num_planes", &fNPlanes, kInt},
+ {"hodo_plane_names", &planenamelist, kString},
+ {"hodo_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
+ {"hodo_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
+ {0}};
+ // fNPlanes = 4; // Default if not defined
+ fTDC_RefTimeCut = 0; // Minimum allowed reference times
fADC_RefTimeCut = 0;
- gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);
+ gHcParms->LoadParmValues((DBRequest*)&listextra, prefix);
- _det_logger->info("Plane Name List : {}" , planenamelist);
- //cout << "Plane Name List : " << planenamelist << endl;
+ _det_logger->info("Plane Name List : {}", planenamelist);
+ // 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;
+ 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()+1];
+ fPlaneNames = new char*[fNPlanes];
+ for (Int_t i = 0; i < fNPlanes; i++) {
+ fPlaneNames[i] = new char[plane_names[i].length() + 1];
strcpy(fPlaneNames[i], plane_names[i].c_str());
}
// Probably shouldn't assume that description is defined
- char* desc = new char[strlen(description)+100];
- fPlanes = new THcScintillatorPlane* [fNPlanes];
- for(Int_t i=0;i < fNPlanes;i++) {
+ 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, this); // Number planes starting from zero!!
- //cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
+ fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i + 1,
+ this); // Number planes starting from zero!!
+ // cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
}
// Save the nominal particle mass
- THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
- fPartMass = app->GetParticleMass();
- fBetaNominal = app->GetBetaAtPcentral();
-
-
+ THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+ fPartMass = app->GetParticleMass();
+ fBetaNominal = app->GetBetaAtPcentral();
if (fSHMS) {
fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("hgcer"));
@@ -154,21 +150,20 @@ void THcHodoscope::Setup(const char* name, const char* description)
fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("cer"));
}
- delete [] desc;
+ delete[] desc;
}
//_____________________________________________________________________________
-THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
-{
+THaAnalysisObject::EStatus THcHodoscope::Init(const TDatime& date) {
// cout << "In THcHodoscope::Init()" << endl;
Setup(GetName(), GetTitle());
char EngineDID[] = "xSCIN";
- EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
- if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
+ EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
+ if (gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0) {
static const char* const here = "Init()";
- //Error( Here(here), "Error filling detectormap for %s.", EngineDID );
- _det_logger->error("THcHodoscope::Init : Error filling detectormap for {}.",EngineDID);
+ // Error( Here(here), "Error filling detectormap for %s.", EngineDID );
+ _det_logger->error("THcHodoscope::Init : Error filling detectormap for {}.", EngineDID);
return kInitError;
}
@@ -177,27 +172,27 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
// But it needs to happen before the sub detectors are initialized
// so that they can get the pointer to the hitlist.
_det_logger->info("Hodo TDC and ADC ref time cut = {} {}", fTDC_RefTimeCut, fADC_RefTimeCut);
- //cout << " Hodo tdc ref time cut = " << fTDC_RefTimeCut << " " << fADC_RefTimeCut << endl;
+ // cout << " Hodo tdc ref time cut = " << fTDC_RefTimeCut << " " << fADC_RefTimeCut << endl;
- InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan()+1,
- fTDC_RefTimeCut, fADC_RefTimeCut);
+ InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan() + 1, fTDC_RefTimeCut,
+ fADC_RefTimeCut);
EStatus status;
// This triggers call of ReadDatabase and DefineVariables
- if( (status = THaNonTrackingDetector::Init( date )) )
- return fStatus=status;
+ if ((status = THaNonTrackingDetector::Init(date)))
+ return fStatus = status;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- if((status = fPlanes[ip]->Init( date ))) {
- return fStatus=status;
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ if ((status = fPlanes[ip]->Init(date))) {
+ return fStatus = status;
}
}
// Why not std::vector?
- fNScinHits = new Int_t [fNPlanes];
- fGoodPlaneTime = new Bool_t [fNPlanes];
- fNPlaneTime = new Int_t [fNPlanes];
- fSumPlaneTime = new Double_t [fNPlanes];
+ fNScinHits = new Int_t[fNPlanes];
+ fGoodPlaneTime = new Bool_t[fNPlanes];
+ fNPlaneTime = new Int_t[fNPlanes];
+ fSumPlaneTime = new Double_t[fNPlanes];
// Double_t fHitCnt4 = 0., fHitCnt3 = 0.;
@@ -207,21 +202,20 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
// fScinHit[m] = new Double_t[fNPaddle[0]];
// }
- for (int ip=0; ip<fNPlanes; ++ip) {
+ for (int ip = 0; ip < fNPlanes; ++ip) {
fScinHitPaddle.push_back(std::vector<Int_t>(fNPaddle[ip], 0));
}
- fPresentP = 0;
- THaVar* vpresent = gHaVars->Find(Form("%s.present",GetApparatus()->GetName()));
- if(vpresent) {
- fPresentP = (Bool_t *) vpresent->GetValuePointer();
+ fPresentP = 0;
+ THaVar* vpresent = gHaVars->Find(Form("%s.present", GetApparatus()->GetName()));
+ if (vpresent) {
+ fPresentP = (Bool_t*)vpresent->GetValuePointer();
}
return fStatus = kOK;
}
//_____________________________________________________________________________
-Int_t THcHodoscope::ReadDatabase( const TDatime& date )
-{
+Int_t THcHodoscope::ReadDatabase(const TDatime& date) {
/**
Read this detector's parameters from the ThcParmlist.
@@ -235,140 +229,139 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
// 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[0] = tolower(GetApparatus()->GetName()[0]);
//
- prefix[1]='\0';
- strcpy(parname,prefix);
- strcat(parname,"scin_");
+ prefix[1] = '\0';
+ strcpy(parname, prefix);
+ strcat(parname, "scin_");
// Int_t plen=strlen(parname);
// cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;
- CreateMissReportParms(Form("%sscin",prefix));
+ CreateMissReportParms(Form("%sscin", prefix));
- fBetaNoTrk = 0.;
+ fBetaNoTrk = 0.;
fBetaNoTrkChiSq = 0.;
- fNPaddle = new UInt_t [fNPlanes];
- fFPTime = new Double_t [fNPlanes];
- fPlaneCenter = new Double_t[fNPlanes];
+ fNPaddle = new UInt_t[fNPlanes];
+ fFPTime = new Double_t[fNPlanes];
+ fPlaneCenter = new Double_t[fNPlanes];
fPlaneSpacing = new Double_t[fNPlanes];
- prefix[0]=tolower(GetApparatus()->GetName()[0]);
+ prefix[0] = tolower(GetApparatus()->GetName()[0]);
//
- prefix[1]='\0';
+ prefix[1] = '\0';
- for(Int_t i=0;i<fNPlanes;i++) {
+ for (Int_t i = 0; i < fNPlanes; i++) {
- DBRequest list[]={
- {Form("scin_%s_nr",fPlaneNames[i]), &fNPaddle[i], kInt},
- {0}
- };
+ DBRequest list[] = {{Form("scin_%s_nr", fPlaneNames[i]), &fNPaddle[i], kInt}, {0}};
-
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ gHcParms->LoadParmValues((DBRequest*)&list, prefix);
}
// GN added
// reading variables from *hodo.param
- fMaxScinPerPlane=fNPaddle[0];
- for (Int_t i=1;i<fNPlanes;i++) {
- fMaxScinPerPlane=(fMaxScinPerPlane > fNPaddle[i])? fMaxScinPerPlane : fNPaddle[i];
+ fMaxScinPerPlane = fNPaddle[0];
+ for (Int_t i = 1; i < fNPlanes; i++) {
+ fMaxScinPerPlane = (fMaxScinPerPlane > fNPaddle[i]) ? fMaxScinPerPlane : fNPaddle[i];
}
// need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
- fMaxHodoScin=fMaxScinPerPlane*fNPlanes;
- if (fDebug>=1) cout <<"fMaxScinPerPlane = "<<fMaxScinPerPlane<<" fMaxHodoScin = "<<fMaxHodoScin<<endl;
-
- fHodoVelLight=new Double_t [fMaxHodoScin];
- fHodoPosSigma=new Double_t [fMaxHodoScin];
- fHodoNegSigma=new Double_t [fMaxHodoScin];
- fHodoPosMinPh=new Double_t [fMaxHodoScin];
- fHodoNegMinPh=new Double_t [fMaxHodoScin];
- fHodoPosPhcCoeff=new Double_t [fMaxHodoScin];
- fHodoNegPhcCoeff=new Double_t [fMaxHodoScin];
- fHodoPosTimeOffset=new Double_t [fMaxHodoScin];
- fHodoNegTimeOffset=new Double_t [fMaxHodoScin];
- fHodoPosPedLimit=new Int_t [fMaxHodoScin];
- fHodoNegPedLimit=new Int_t [fMaxHodoScin];
- fHodoPosInvAdcOffset=new Double_t [fMaxHodoScin];
- fHodoNegInvAdcOffset=new Double_t [fMaxHodoScin];
- fHodoPosInvAdcLinear=new Double_t [fMaxHodoScin];
- fHodoNegInvAdcLinear=new Double_t [fMaxHodoScin];
- fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin];
- fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin];
-
- //New Time-Walk Calibration Parameters
- fHodoVelFit=new Double_t [fMaxHodoScin];
- fHodoCableFit=new Double_t [fMaxHodoScin];
- fHodo_LCoeff=new Double_t [fMaxHodoScin];
- fHodoPos_c1=new Double_t [fMaxHodoScin];
- fHodoNeg_c1=new Double_t [fMaxHodoScin];
- fHodoPos_c2=new Double_t [fMaxHodoScin];
- fHodoNeg_c2=new Double_t [fMaxHodoScin];
- fHodoSigmaPos=new Double_t [fMaxHodoScin];
- fHodoSigmaNeg=new Double_t [fMaxHodoScin];
-
- fNHodoscopes = 2;
- fxLoScin = new Int_t [fNHodoscopes];
- fxHiScin = new Int_t [fNHodoscopes];
- fyLoScin = new Int_t [fNHodoscopes];
- fyHiScin = new Int_t [fNHodoscopes];
- fHodoSlop = new Double_t [fNPlanes];
- fTdcOffset = new Int_t [fNPlanes];
- fAdcTdcOffset = new Double_t [fNPlanes];
- fHodoPosAdcTimeWindowMin = new Double_t [fMaxHodoScin];
- fHodoPosAdcTimeWindowMax = new Double_t [fMaxHodoScin];
- fHodoNegAdcTimeWindowMin = new Double_t [fMaxHodoScin];
- fHodoNegAdcTimeWindowMax = new Double_t [fMaxHodoScin];
-
-
- for(Int_t ip=0;ip<fNPlanes;ip++) { // Set a large default window
- fTdcOffset[ip] = 0 ;
- fAdcTdcOffset[ip] = 0.0 ;
+ fMaxHodoScin = fMaxScinPerPlane * fNPlanes;
+ if (fDebug >= 1)
+ cout << "fMaxScinPerPlane = " << fMaxScinPerPlane << " fMaxHodoScin = " << fMaxHodoScin << endl;
+
+ fHodoVelLight = new Double_t[fMaxHodoScin];
+ fHodoPosSigma = new Double_t[fMaxHodoScin];
+ fHodoNegSigma = new Double_t[fMaxHodoScin];
+ fHodoPosMinPh = new Double_t[fMaxHodoScin];
+ fHodoNegMinPh = new Double_t[fMaxHodoScin];
+ fHodoPosPhcCoeff = new Double_t[fMaxHodoScin];
+ fHodoNegPhcCoeff = new Double_t[fMaxHodoScin];
+ fHodoPosTimeOffset = new Double_t[fMaxHodoScin];
+ fHodoNegTimeOffset = new Double_t[fMaxHodoScin];
+ fHodoPosPedLimit = new Int_t[fMaxHodoScin];
+ fHodoNegPedLimit = new Int_t[fMaxHodoScin];
+ fHodoPosInvAdcOffset = new Double_t[fMaxHodoScin];
+ fHodoNegInvAdcOffset = new Double_t[fMaxHodoScin];
+ fHodoPosInvAdcLinear = new Double_t[fMaxHodoScin];
+ fHodoNegInvAdcLinear = new Double_t[fMaxHodoScin];
+ fHodoPosInvAdcAdc = new Double_t[fMaxHodoScin];
+ fHodoNegInvAdcAdc = new Double_t[fMaxHodoScin];
+
+ // New Time-Walk Calibration Parameters
+ fHodoVelFit = new Double_t[fMaxHodoScin];
+ fHodoCableFit = new Double_t[fMaxHodoScin];
+ fHodo_LCoeff = new Double_t[fMaxHodoScin];
+ fHodoPos_c1 = new Double_t[fMaxHodoScin];
+ fHodoNeg_c1 = new Double_t[fMaxHodoScin];
+ fHodoPos_c2 = new Double_t[fMaxHodoScin];
+ fHodoNeg_c2 = new Double_t[fMaxHodoScin];
+ fHodoSigmaPos = new Double_t[fMaxHodoScin];
+ fHodoSigmaNeg = new Double_t[fMaxHodoScin];
+
+ fNHodoscopes = 2;
+ fxLoScin = new Int_t[fNHodoscopes];
+ fxHiScin = new Int_t[fNHodoscopes];
+ fyLoScin = new Int_t[fNHodoscopes];
+ fyHiScin = new Int_t[fNHodoscopes];
+ fHodoSlop = new Double_t[fNPlanes];
+ fTdcOffset = new Int_t[fNPlanes];
+ fAdcTdcOffset = new Double_t[fNPlanes];
+ fHodoPosAdcTimeWindowMin = new Double_t[fMaxHodoScin];
+ fHodoPosAdcTimeWindowMax = new Double_t[fMaxHodoScin];
+ fHodoNegAdcTimeWindowMin = new Double_t[fMaxHodoScin];
+ fHodoNegAdcTimeWindowMax = new Double_t[fMaxHodoScin];
+
+ for (Int_t ip = 0; ip < fNPlanes; ip++) { // Set a large default window
+ fTdcOffset[ip] = 0;
+ fAdcTdcOffset[ip] = 0.0;
}
- DBRequest list[]={
- {"cosmicflag", &fCosmicFlag, kInt, 0, 1},
- {"NumPlanesBetaCalc", &fNumPlanesBetaCalc, kInt, 0, 1},
- {"start_time_center", &fStartTimeCenter, kDouble},
- {"start_time_slop", &fStartTimeSlop, kDouble},
- {"scin_tdc_to_time", &fScinTdcToTime, kDouble},
- {"scin_tdc_min", &fScinTdcMin, kDouble},
- {"scin_tdc_max", &fScinTdcMax, kDouble},
- {"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
- {"pathlength_central", &fPathLengthCentral, kDouble},
- {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1},
- {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1},
- {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1},
- {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1},
- {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
- {"xloscin", &fxLoScin[0], kInt, (UInt_t) fNHodoscopes},
- {"xhiscin", &fxHiScin[0], kInt, (UInt_t) fNHodoscopes},
- {"yloscin", &fyLoScin[0], kInt, (UInt_t) fNHodoscopes},
- {"yhiscin", &fyHiScin[0], kInt, (UInt_t) fNHodoscopes},
- {"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
- {"cer_npe", &fNCerNPE, kDouble, 0, 1},
- {"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
- {"hodo_slop", fHodoSlop, kDouble, (UInt_t) fNPlanes},
- {"debugprintscinraw", &fdebugprintscinraw, kInt, 0,1},
- {"hodo_tdc_offset", fTdcOffset, kInt, (UInt_t) fNPlanes, 1},
- {"hodo_adc_tdc_offset", fAdcTdcOffset, kDouble, (UInt_t) fNPlanes, 1},
- {"hodo_PosAdcTimeWindowMin", fHodoPosAdcTimeWindowMin, kDouble, (UInt_t) fMaxHodoScin, 1},
- {"hodo_PosAdcTimeWindowMax", fHodoPosAdcTimeWindowMax, kDouble, (UInt_t) fMaxHodoScin, 1},
- {"hodo_NegAdcTimeWindowMin", fHodoNegAdcTimeWindowMin, kDouble, (UInt_t) fMaxHodoScin, 1},
- {"hodo_NegAdcTimeWindowMax", fHodoNegAdcTimeWindowMax, kDouble, (UInt_t) fMaxHodoScin, 1},
- {"dumptof", &fDumpTOF, kInt, 0, 1},
- {"TOFCalib_shtrk_lo", &fTOFCalib_shtrk_lo, kDouble, 0, 1},
- {"TOFCalib_shtrk_hi", &fTOFCalib_shtrk_hi, kDouble, 0, 1},
- {"TOFCalib_cer_lo", &fTOFCalib_cer_lo, kDouble, 0, 1},
- {"TOFCalib_beta_lo", &fTOFCalib_beta_lo, kDouble, 0, 1},
- {"TOFCalib_beta_hi", &fTOFCalib_beta_hi, kDouble, 0, 1},
- {"dumptof_filename", &fTOFDumpFile, kString, 0, 1},
- {0}
- };
+ DBRequest list[] = {
+ {"cosmicflag", &fCosmicFlag, kInt, 0, 1},
+ {"NumPlanesBetaCalc", &fNumPlanesBetaCalc, kInt, 0, 1},
+ {"start_time_center", &fStartTimeCenter, kDouble},
+ {"start_time_slop", &fStartTimeSlop, kDouble},
+ {"scin_tdc_to_time", &fScinTdcToTime, kDouble},
+ {"scin_tdc_min", &fScinTdcMin, kDouble},
+ {"scin_tdc_max", &fScinTdcMax, kDouble},
+ {"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
+ {"pathlength_central", &fPathLengthCentral, kDouble},
+ {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1},
+ {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1},
+ {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
+ {"xloscin", &fxLoScin[0], kInt, (UInt_t)fNHodoscopes},
+ {"xhiscin", &fxHiScin[0], kInt, (UInt_t)fNHodoscopes},
+ {"yloscin", &fyLoScin[0], kInt, (UInt_t)fNHodoscopes},
+ {"yhiscin", &fyHiScin[0], kInt, (UInt_t)fNHodoscopes},
+ {"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
+ {"trackeff_scint_ydiff_max", &trackeff_scint_ydiff_max, kDouble, 0, 1},
+ {"trackeff_scint_xdiff_max", &trackeff_scint_xdiff_max, kDouble, 0, 1},
+ {"cer_npe", &fNCerNPE, kDouble, 0, 1},
+ {"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
+ {"hodo_slop", fHodoSlop, kDouble, (UInt_t)fNPlanes},
+ {"debugprintscinraw", &fdebugprintscinraw, kInt, 0, 1},
+ {"hodo_tdc_offset", fTdcOffset, kInt, (UInt_t)fNPlanes, 1},
+ {"hodo_adc_tdc_offset", fAdcTdcOffset, kDouble, (UInt_t)fNPlanes, 1},
+ {"hodo_PosAdcTimeWindowMin", fHodoPosAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
+ {"hodo_PosAdcTimeWindowMax", fHodoPosAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
+ {"hodo_NegAdcTimeWindowMin", fHodoNegAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
+ {"hodo_NegAdcTimeWindowMax", fHodoNegAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
+ {"dumptof", &fDumpTOF, kInt, 0, 1},
+ {"TOFCalib_shtrk_lo", &fTOFCalib_shtrk_lo, kDouble, 0, 1},
+ {"TOFCalib_shtrk_hi", &fTOFCalib_shtrk_hi, kDouble, 0, 1},
+ {"TOFCalib_cer_lo", &fTOFCalib_cer_lo, kDouble, 0, 1},
+ {"TOFCalib_beta_lo", &fTOFCalib_beta_lo, kDouble, 0, 1},
+ {"TOFCalib_beta_hi", &fTOFCalib_beta_hi, kDouble, 0, 1},
+ {"dumptof_filename", &fTOFDumpFile, kString, 0, 1},
+ {0}};
// Defaults if not defined in parameter file
- for(UInt_t ip=0;ip<fMaxHodoScin;ip++) {
+ trackeff_scint_ydiff_max = 20.;
+ trackeff_scint_xdiff_max = 20.;
+ for (UInt_t ip = 0; ip < fMaxHodoScin; ip++) {
fHodoPosAdcTimeWindowMin[ip] = -1000.;
fHodoPosAdcTimeWindowMax[ip] = 1000.;
fHodoNegAdcTimeWindowMin[ip] = -1000.;
@@ -376,33 +369,34 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
fHodoPosPedLimit[ip] = 0.0;
fHodoNegPedLimit[ip] = 0.0;
- fHodoPosSigma[ip] = 0.2;
- fHodoNegSigma[ip] = 0.2;
+ fHodoPosSigma[ip] = 0.2;
+ fHodoNegSigma[ip] = 0.2;
}
- fTOFCalib_shtrk_lo=-kBig;
- fTOFCalib_shtrk_hi= kBig;
- fTOFCalib_cer_lo=-kBig;
- fTOFCalib_beta_lo=-kBig;
- fTOFCalib_beta_hi= kBig;
- fdebugprintscinraw=0;
- fDumpTOF = 0;
- fTOFDumpFile="";
- fTofUsingInvAdc = 1;
- fTofTolerance = 3.0;
- fNCerNPE = 2.0;
- fNormETot = 0.7;
- fCosmicFlag=0;
- fNumPlanesBetaCalc=4;
+ fTOFCalib_shtrk_lo = -kBig;
+ fTOFCalib_shtrk_hi = kBig;
+ fTOFCalib_cer_lo = -kBig;
+ fTOFCalib_beta_lo = -kBig;
+ fTOFCalib_beta_hi = kBig;
+ fdebugprintscinraw = 0;
+ fDumpTOF = 0;
+ fTOFDumpFile = "";
+ fTofUsingInvAdc = 1;
+ fTofTolerance = 3.0;
+ fNCerNPE = 2.0;
+ fNormETot = 0.7;
+ fCosmicFlag = 0;
+ fNumPlanesBetaCalc = 4;
// Gets added to each reference time corrected raw TDC value
// to make sure valid range is all positive.
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- if (fCosmicFlag==1) cout << "Setup for cosmics in TOF"<< endl;
+ gHcParms->LoadParmValues((DBRequest*)&list, prefix);
+ if (fCosmicFlag == 1)
+ cout << "Setup for cosmics in TOF" << endl;
// cout << " cosmic flag = " << fCosmicFlag << endl;
- if(fDumpTOF) {
+ if (fDumpTOF) {
fDumpOut.open(fTOFDumpFile.c_str());
- if(fDumpOut.is_open()) {
- //fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
+ if (fDumpOut.is_open()) {
+ // fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
} else {
fDumpTOF = 0;
cout << "WARNING: Unable to open TOF Dump file " << fTOFDumpFile << endl;
@@ -427,24 +421,23 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
// << " number of photo electrons = " << fNCerNPE
// << endl;
- //Set scin Velocity/Cable to default
- for (UInt_t i=0; i<fMaxHodoScin; i++) {
+ // Set scin Velocity/Cable to default
+ for (UInt_t i = 0; i < fMaxHodoScin; i++) {
fHodoVelLight[i] = 15.0;
}
if (fTofUsingInvAdc) {
- DBRequest list2[]={
- {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin, 1},
- {"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin},
- {"hodo_neg_invadc_offset",&fHodoNegInvAdcOffset[0],kDouble,fMaxHodoScin},
- {"hodo_pos_invadc_linear",&fHodoPosInvAdcLinear[0],kDouble,fMaxHodoScin},
- {"hodo_neg_invadc_linear",&fHodoNegInvAdcLinear[0],kDouble,fMaxHodoScin},
- {"hodo_pos_invadc_adc",&fHodoPosInvAdcAdc[0],kDouble,fMaxHodoScin},
- {"hodo_neg_invadc_adc",&fHodoNegInvAdcAdc[0],kDouble,fMaxHodoScin},
- {0}
- };
-
- gHcParms->LoadParmValues((DBRequest*)&list2,prefix);
+ DBRequest list2[] = {
+ {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_pos_invadc_offset", &fHodoPosInvAdcOffset[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_invadc_offset", &fHodoNegInvAdcOffset[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_invadc_linear", &fHodoPosInvAdcLinear[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_invadc_linear", &fHodoNegInvAdcLinear[0], kDouble, fMaxHodoScin},
+ {"hodo_pos_invadc_adc", &fHodoPosInvAdcAdc[0], kDouble, fMaxHodoScin},
+ {"hodo_neg_invadc_adc", &fHodoNegInvAdcAdc[0], kDouble, fMaxHodoScin},
+ {0}};
+
+ gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
};
/* if (!fTofUsingInvAdc) {
DBRequest list3[]={
@@ -455,65 +448,60 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
{"hodo_neg_phc_coeff", &fHodoNegPhcCoeff[0], kDouble, fMaxHodoScin},
{"hodo_pos_time_offset", &fHodoPosTimeOffset[0], kDouble, fMaxHodoScin},
{"hodo_neg_time_offset", &fHodoNegTimeOffset[0], kDouble, fMaxHodoScin},
- {0}
+ {0}
};
-
-
+
+
gHcParms->LoadParmValues((DBRequest*)&list3,prefix);
*/
- DBRequest list4[]={
- {"hodo_velFit", &fHodoVelFit[0], kDouble, fMaxHodoScin, 1},
- {"hodo_cableFit", &fHodoCableFit[0], kDouble, fMaxHodoScin, 1},
- {"hodo_LCoeff", &fHodo_LCoeff[0], kDouble, fMaxHodoScin, 1},
- {"c1_Pos", &fHodoPos_c1[0], kDouble, fMaxHodoScin, 1},
- {"c1_Neg", &fHodoNeg_c1[0], kDouble, fMaxHodoScin, 1},
- {"c2_Pos", &fHodoPos_c2[0], kDouble, fMaxHodoScin, 1},
- {"c2_Neg", &fHodoNeg_c2[0], kDouble, fMaxHodoScin, 1},
- {"TDC_threshold", &fTdc_Thrs, kDouble, 0, 1},
- {"hodo_PosSigma", &fHodoSigmaPos[0], kDouble, fMaxHodoScin, 1},
- {"hodo_NegSigma", &fHodoSigmaNeg[0], kDouble, fMaxHodoScin, 1},
- {0}
- };
-
- fTdc_Thrs = 1.0;
- //Set Default Values if NOT defined in param file
- for (UInt_t i=0; i<fMaxHodoScin; i++)
- {
-
- //Turn OFF Time-Walk Correction if param file NOT found
- fHodoPos_c1[i] = 0.0;
- fHodoPos_c2[i] = 0.0;
- fHodoNeg_c1[i] = 0.0;
- fHodoNeg_c2[i] = 0.0;
- }
- for (UInt_t i=0; i<fMaxHodoScin; i++)
- {
- //Set scin Velocity/Cable to default
- fHodoCableFit[i] = 0.0;
- fHodoVelFit[i] = 15.0;
- //set time coeff between paddles to default
- fHodo_LCoeff[i] = 0.0;
-
- }
-
- gHcParms->LoadParmValues((DBRequest*)&list4,prefix);
-
- if (fDebug >=1) {
- cout <<"******* Testing Hodoscope Parameter Reading ***\n";
- cout<<"StarTimeCenter = "<<fStartTimeCenter<<endl;
- cout<<"StartTimeSlop = "<<fStartTimeSlop<<endl;
- cout <<"ScintTdcToTime = "<<fScinTdcToTime<<endl;
- cout <<"TdcMin = "<<fScinTdcMin<<" TdcMax = "<<fScinTdcMax<<endl;
- cout <<"TofTolerance = "<<fTofTolerance<<endl;
- cout <<"*** VelLight ***\n";
- for (Int_t i1=0;i1<fNPlanes;i1++) {
- cout<<"Plane "<<i1<<endl;
- for (UInt_t i2=0;i2<fMaxScinPerPlane;i2++) {
- cout<<fHodoVelLight[GetScinIndex(i1,i2)]<<" ";
+ DBRequest list4[] = {{"hodo_velFit", &fHodoVelFit[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_cableFit", &fHodoCableFit[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_LCoeff", &fHodo_LCoeff[0], kDouble, fMaxHodoScin, 1},
+ {"c1_Pos", &fHodoPos_c1[0], kDouble, fMaxHodoScin, 1},
+ {"c1_Neg", &fHodoNeg_c1[0], kDouble, fMaxHodoScin, 1},
+ {"c2_Pos", &fHodoPos_c2[0], kDouble, fMaxHodoScin, 1},
+ {"c2_Neg", &fHodoNeg_c2[0], kDouble, fMaxHodoScin, 1},
+ {"TDC_threshold", &fTdc_Thrs, kDouble, 0, 1},
+ {"hodo_PosSigma", &fHodoSigmaPos[0], kDouble, fMaxHodoScin, 1},
+ {"hodo_NegSigma", &fHodoSigmaNeg[0], kDouble, fMaxHodoScin, 1},
+ {0}};
+
+ fTdc_Thrs = 1.0;
+ // Set Default Values if NOT defined in param file
+ for (UInt_t i = 0; i < fMaxHodoScin; i++) {
+
+ // Turn OFF Time-Walk Correction if param file NOT found
+ fHodoPos_c1[i] = 0.0;
+ fHodoPos_c2[i] = 0.0;
+ fHodoNeg_c1[i] = 0.0;
+ fHodoNeg_c2[i] = 0.0;
+ }
+ for (UInt_t i = 0; i < fMaxHodoScin; i++) {
+ // Set scin Velocity/Cable to default
+ fHodoCableFit[i] = 0.0;
+ fHodoVelFit[i] = 15.0;
+ // set time coeff between paddles to default
+ fHodo_LCoeff[i] = 0.0;
+ }
+
+ gHcParms->LoadParmValues((DBRequest*)&list4, prefix);
+
+ if (fDebug >= 1) {
+ cout << "******* Testing Hodoscope Parameter Reading ***\n";
+ cout << "StarTimeCenter = " << fStartTimeCenter << endl;
+ cout << "StartTimeSlop = " << fStartTimeSlop << endl;
+ cout << "ScintTdcToTime = " << fScinTdcToTime << endl;
+ cout << "TdcMin = " << fScinTdcMin << " TdcMax = " << fScinTdcMax << endl;
+ cout << "TofTolerance = " << fTofTolerance << endl;
+ cout << "*** VelLight ***\n";
+ for (Int_t i1 = 0; i1 < fNPlanes; i1++) {
+ cout << "Plane " << i1 << endl;
+ for (UInt_t i2 = 0; i2 < fMaxScinPerPlane; i2++) {
+ cout << fHodoVelLight[GetScinIndex(i1, i2)] << " ";
}
- cout <<endl;
+ cout << endl;
}
- cout <<endl<<endl;
+ cout << endl << endl;
// check fHodoPosPhcCoeff
/*
cout <<"fHodoPosPhcCoeff = ";
@@ -525,87 +513,96 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
}
//
if ((fTofTolerance > 0.5) && (fTofTolerance < 10000.)) {
- //cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
- _det_logger->info("THcHodoscope: Using {} nsec window for fp no_track calculations.",fTofTolerance);
- }
- else {
- fTofTolerance= 3.0;
- //cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
- _det_logger->warn("THcHodoscope: Using default {} nsec window for fp no_track calculations.",fTofTolerance);
+ // cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
+ _det_logger->info("THcHodoscope: Using {} nsec window for fp no_track calculations.",
+ fTofTolerance);
+ } else {
+ fTofTolerance = 3.0;
+ // cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
+ _det_logger->warn("THcHodoscope: Using default {} nsec window for fp no_track calculations.",
+ fTofTolerance);
}
+ fRatio_xpfp_to_xfp = 0.00;
+ TString SHMS = "p";
+ TString HMS = "h";
+ TString test = prefix[0];
+ if (test == SHMS)
+ fRatio_xpfp_to_xfp = 0.0018; // SHMS
+ if (test == HMS)
+ fRatio_xpfp_to_xfp = 0.0011; // HMS
+ cout << " fRatio_xpfp_to_xfp= " << fRatio_xpfp_to_xfp << endl;
+
fIsInit = true;
return kOK;
}
//_____________________________________________________________________________
-Int_t THcHodoscope::DefineVariables( EMode mode )
-{
+Int_t THcHodoscope::DefineVariables(EMode mode) {
/**
Initialize global variables for histograms and Root tree
*/
// cout << "THcHodoscope::DefineVariables called " << GetName() << endl;
- if( mode == kDefine && fIsSetup ) return kOK;
- fIsSetup = ( mode == kDefine );
+ if (mode == kDefine && fIsSetup)
+ return kOK;
+ fIsSetup = (mode == kDefine);
// Register variables in global list
RVarDef vars[] = {
- // Move these into THcHallCSpectrometer using track fTracks
- {"beta" , "Beta including track info" , "fBeta" } ,
- {"betanotrack" , "Beta from scintillator hits" , "fBetaNoTrk" } ,
- {"betachisqnotrack" , "Chi square of beta from scintillator hits" , "fBetaNoTrkChiSq" } ,
- {"fpHitsTime" , "Time at focal plane from all hits" , "fFPTimeAll" } ,
- {"starttime" , "Hodoscope Start Time" , "fStartTime" } ,
- {"goodstarttime" , "Hodoscope Good Start Time (logical flag)" , "fGoodStartTime" } ,
- {"goodscinhit" , "Hit in fid area" , "fGoodScinHits" } ,
- {"TimeHist_Sigma" , "" , "fTimeHist_Sigma" } ,
- {"TimeHist_Peak" , "" , "fTimeHist_Peak" } ,
- {"TimeHist_Hits" , "" , "fTimeHist_Hits" } ,
- { 0 }
- };
-
- return DefineVarsFromList( vars, mode );
+ // Move these into THcHallCSpectrometer using track fTracks
+ {"beta", "Beta including track info", "fBeta"},
+ {"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
+ {"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
+ {"fpHitsTime", "Time at focal plane from all hits", "fFPTimeAll"},
+ {"starttime", "Hodoscope Start Time", "fStartTime"},
+ {"goodstarttime", "Hodoscope Good Start Time (logical flag)", "fGoodStartTime"},
+ {"goodscinhit", "Hit in fid area", "fGoodScinHits"},
+ {"TimeHist_Sigma", "", "fTimeHist_Sigma"},
+ {"TimeHist_Peak", "", "fTimeHist_Peak"},
+ {"TimeHist_Hits", "", "fTimeHist_Hits"},
+ {0}};
+
+ return DefineVarsFromList(vars, mode);
// return kOK;
}
//_____________________________________________________________________________
-Int_t THcHodoscope::ManualInitTree( TTree* t ){
+Int_t THcHodoscope::ManualInitTree(TTree* t) {
// The most direct path to the output tree!!!
std::string app_name = GetApparatus()->GetName();
std::string det_name = GetName();
std::string branch_name = (app_name + "_" + det_name + "_data");
if (t) {
- _det_logger->info("THcHodoscope::ManualInitTree : Adding branch, {}, to output tree", branch_name);
+ _det_logger->info("THcHodoscope::ManualInitTree : Adding branch, {}, to output tree",
+ branch_name);
t->Branch(branch_name.c_str(), &_basic_data, 32000, 99);
}
return 0;
}
//_____________________________________________________________________________
-THcHodoscope::~THcHodoscope()
-{
+THcHodoscope::~THcHodoscope() {
// Destructor. Remove variables from global list.
- delete [] fFPTime;
- delete [] fPlaneCenter;
- delete [] fPlaneSpacing;
+ delete[] fFPTime;
+ delete[] fPlaneCenter;
+ delete[] fPlaneSpacing;
- if( fIsSetup )
+ if (fIsSetup)
RemoveVariables();
- if( fIsInit )
+ if (fIsInit)
DeleteArrays();
- for( int i = 0; i < fNPlanes; ++i ) {
+ for (int i = 0; i < fNPlanes; ++i) {
delete fPlanes[i];
- delete [] fPlaneNames[i];
+ delete[] fPlaneNames[i];
}
- delete [] fPlanes;
- delete [] fPlaneNames;
+ delete[] fPlanes;
+ delete[] fPlaneNames;
}
//_____________________________________________________________________________
-void THcHodoscope::DeleteArrays()
-{
+void THcHodoscope::DeleteArrays() {
// Delete member arrays. Used by destructor.
// Int_t k;
// for( k = 0; k < fNPlanes; k++){
@@ -613,80 +610,121 @@ void THcHodoscope::DeleteArrays()
// }
// delete [] fScinHit;
- delete [] fxLoScin; fxLoScin = NULL;
- delete [] fxHiScin; fxHiScin = NULL;
- delete [] fyLoScin; fyLoScin = NULL;
- delete [] fyHiScin; fyHiScin = NULL;
- delete [] fHodoSlop; fHodoSlop = NULL;
-
- delete [] fNPaddle; fNPaddle = NULL;
- delete [] fHodoVelLight; fHodoVelLight = NULL;
- delete [] fHodoPosSigma; fHodoPosSigma = NULL;
- delete [] fHodoNegSigma; fHodoNegSigma = NULL;
- delete [] fHodoPosMinPh; fHodoPosMinPh = NULL;
- delete [] fHodoNegMinPh; fHodoNegMinPh = NULL;
- delete [] fHodoPosPhcCoeff; fHodoPosPhcCoeff = NULL;
- delete [] fHodoNegPhcCoeff; fHodoNegPhcCoeff = NULL;
- delete [] fHodoPosTimeOffset; fHodoPosTimeOffset = NULL;
- delete [] fHodoNegTimeOffset; fHodoNegTimeOffset = NULL;
- delete [] fHodoPosPedLimit; fHodoPosPedLimit = NULL;
- delete [] fHodoNegPedLimit; fHodoNegPedLimit = NULL;
- delete [] fHodoPosInvAdcOffset; fHodoPosInvAdcOffset = NULL;
- delete [] fHodoNegInvAdcOffset; fHodoNegInvAdcOffset = NULL;
- delete [] fHodoPosInvAdcLinear; fHodoPosInvAdcLinear = NULL;
- delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL;
- delete [] fHodoPosInvAdcAdc; fHodoPosInvAdcAdc = NULL;
- delete [] fHodoNegInvAdcAdc; fHodoNegInvAdcAdc = NULL;
- delete [] fGoodPlaneTime; fGoodPlaneTime = NULL;
- delete [] fNPlaneTime; fNPlaneTime = NULL;
- delete [] fSumPlaneTime; fSumPlaneTime = NULL;
- delete [] fNScinHits; fNScinHits = NULL;
- delete [] fTdcOffset; fTdcOffset = NULL;
- delete [] fAdcTdcOffset; fAdcTdcOffset = NULL;
- delete [] fHodoNegAdcTimeWindowMin; fHodoNegAdcTimeWindowMin = NULL;
- delete [] fHodoNegAdcTimeWindowMax; fHodoNegAdcTimeWindowMax = NULL;
- delete [] fHodoPosAdcTimeWindowMin; fHodoPosAdcTimeWindowMin = NULL;
- delete [] fHodoPosAdcTimeWindowMax; fHodoPosAdcTimeWindowMax = NULL;
-
- delete [] fHodoVelFit; fHodoVelFit = NULL;
- delete [] fHodoCableFit; fHodoCableFit = NULL;
- delete [] fHodo_LCoeff; fHodo_LCoeff = NULL;
- delete [] fHodoPos_c1; fHodoPos_c1 = NULL;
- delete [] fHodoNeg_c1; fHodoNeg_c1 = NULL;
- delete [] fHodoPos_c2; fHodoPos_c2 = NULL;
- delete [] fHodoNeg_c2; fHodoNeg_c2 = NULL;
- delete [] fHodoSigmaPos; fHodoSigmaPos = NULL;
- delete [] fHodoSigmaNeg; fHodoSigmaNeg = NULL;
+ delete[] fxLoScin;
+ fxLoScin = NULL;
+ delete[] fxHiScin;
+ fxHiScin = NULL;
+ delete[] fyLoScin;
+ fyLoScin = NULL;
+ delete[] fyHiScin;
+ fyHiScin = NULL;
+ delete[] fHodoSlop;
+ fHodoSlop = NULL;
+
+ delete[] fNPaddle;
+ fNPaddle = NULL;
+ delete[] fHodoVelLight;
+ fHodoVelLight = NULL;
+ delete[] fHodoPosSigma;
+ fHodoPosSigma = NULL;
+ delete[] fHodoNegSigma;
+ fHodoNegSigma = NULL;
+ delete[] fHodoPosMinPh;
+ fHodoPosMinPh = NULL;
+ delete[] fHodoNegMinPh;
+ fHodoNegMinPh = NULL;
+ delete[] fHodoPosPhcCoeff;
+ fHodoPosPhcCoeff = NULL;
+ delete[] fHodoNegPhcCoeff;
+ fHodoNegPhcCoeff = NULL;
+ delete[] fHodoPosTimeOffset;
+ fHodoPosTimeOffset = NULL;
+ delete[] fHodoNegTimeOffset;
+ fHodoNegTimeOffset = NULL;
+ delete[] fHodoPosPedLimit;
+ fHodoPosPedLimit = NULL;
+ delete[] fHodoNegPedLimit;
+ fHodoNegPedLimit = NULL;
+ delete[] fHodoPosInvAdcOffset;
+ fHodoPosInvAdcOffset = NULL;
+ delete[] fHodoNegInvAdcOffset;
+ fHodoNegInvAdcOffset = NULL;
+ delete[] fHodoPosInvAdcLinear;
+ fHodoPosInvAdcLinear = NULL;
+ delete[] fHodoNegInvAdcLinear;
+ fHodoNegInvAdcLinear = NULL;
+ delete[] fHodoPosInvAdcAdc;
+ fHodoPosInvAdcAdc = NULL;
+ delete[] fHodoNegInvAdcAdc;
+ fHodoNegInvAdcAdc = NULL;
+ delete[] fGoodPlaneTime;
+ fGoodPlaneTime = NULL;
+ delete[] fNPlaneTime;
+ fNPlaneTime = NULL;
+ delete[] fSumPlaneTime;
+ fSumPlaneTime = NULL;
+ delete[] fNScinHits;
+ fNScinHits = NULL;
+ delete[] fTdcOffset;
+ fTdcOffset = NULL;
+ delete[] fAdcTdcOffset;
+ fAdcTdcOffset = NULL;
+ delete[] fHodoNegAdcTimeWindowMin;
+ fHodoNegAdcTimeWindowMin = NULL;
+ delete[] fHodoNegAdcTimeWindowMax;
+ fHodoNegAdcTimeWindowMax = NULL;
+ delete[] fHodoPosAdcTimeWindowMin;
+ fHodoPosAdcTimeWindowMin = NULL;
+ delete[] fHodoPosAdcTimeWindowMax;
+ fHodoPosAdcTimeWindowMax = NULL;
+
+ delete[] fHodoVelFit;
+ fHodoVelFit = NULL;
+ delete[] fHodoCableFit;
+ fHodoCableFit = NULL;
+ delete[] fHodo_LCoeff;
+ fHodo_LCoeff = NULL;
+ delete[] fHodoPos_c1;
+ fHodoPos_c1 = NULL;
+ delete[] fHodoNeg_c1;
+ fHodoNeg_c1 = NULL;
+ delete[] fHodoPos_c2;
+ fHodoPos_c2 = NULL;
+ delete[] fHodoNeg_c2;
+ fHodoNeg_c2 = NULL;
+ delete[] fHodoSigmaPos;
+ fHodoSigmaPos = NULL;
+ delete[] fHodoSigmaNeg;
+ fHodoSigmaNeg = NULL;
}
//_____________________________________________________________________________
-void THcHodoscope::Clear( Option_t* opt )
-{
+void THcHodoscope::Clear(Option_t* opt) {
/*! \brief Clears variables
*
* Called by THcHodoscope::Decode
*
*/
- fTimeHist_Sigma= kBig;
- fTimeHist_Peak= kBig;
- fTimeHist_Hits= kBig;
+ fTimeHist_Sigma = kBig;
+ fTimeHist_Peak = kBig;
+ fTimeHist_Hits = kBig;
- fBeta = 0.0;
- fBetaNoTrk = 0.0;
+ fBeta = 0.0;
+ fBetaNoTrk = 0.0;
fBetaNoTrkChiSq = 0.0;
- fStartTime = -1000.;
- fFPTimeAll= -1000.;
- fGoodStartTime = kFALSE;
- fGoodScinHits = 0;
-
- if( *opt != 'I' ) {
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fStartTime = -1000.;
+ fFPTimeAll = -1000.;
+ fGoodStartTime = kFALSE;
+ fGoodScinHits = 0;
+
+ if (*opt != 'I') {
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
fPlanes[ip]->Clear();
- fFPTime[ip]=0.;
- fPlaneCenter[ip]=0.;
- fPlaneSpacing[ip]=0.;
- for(UInt_t iPaddle=0;iPaddle<fNPaddle[ip]; ++iPaddle) {
- fScinHitPaddle[ip][iPaddle]=0;
+ fFPTime[ip] = 0.;
+ fPlaneCenter[ip] = 0.;
+ fPlaneSpacing[ip] = 0.;
+ for (UInt_t iPaddle = 0; iPaddle < fNPaddle[ip]; ++iPaddle) {
+ fScinHitPaddle[ip][iPaddle] = 0;
}
}
}
@@ -701,23 +739,24 @@ void THcHodoscope::Clear( Option_t* opt )
}
//_____________________________________________________________________________
-Int_t THcHodoscope::Decode( const THaEvData& evdata )
-{
- /*! \brief Decodes raw data and processes raw data into hits for each instance of THcScintillatorPlane
+Int_t THcHodoscope::Decode(const THaEvData& evdata) {
+ /*! \brief Decodes raw data and processes raw data into hits for each instance of
+ * THcScintillatorPlane
*
* - Reads raw data using THcHitList::DecodeToHitList
* - If one wants to subtract pedestals (assumed to be a set of data at beginning of run)
* + Must define "Pedestal_event" cut in the cuts definition file
* + For each "Pedestal_event" calls THcScintillatorPlane::AccumulatePedestals and returns
- * + After First event which is not a "Pedestal_event" calls THcScintillatorPlane::CalculatePedestals
+ * + After First event which is not a "Pedestal_event" calls
+ * THcScintillatorPlane::CalculatePedestals
* - For each scintillator plane THcScintillatorPlane::ProcessHits
* - Calls THcHodoscope::EstimateFocalPlaneTime
*
*
*/
// Get the Hall C style hitlist (fRawHitList) for this event
- Bool_t present = kTRUE; // Suppress reference time warnings
- if(fPresentP) { // if this spectrometer not part of trigger
+ Bool_t present = kTRUE; // Suppress reference time warnings
+ if (fPresentP) { // if this spectrometer not part of trigger
present = *fPresentP;
}
fNHits = DecodeToHitList(evdata, !present);
@@ -728,308 +767,304 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
// cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;
fCheckEvent = evdata.GetEvNum();
- fEventType = evdata.GetEvType();
+ fEventType = evdata.GetEvType();
- if(gHaCuts->Result("Pedestal_event")) {
+ if (gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
}
- fAnalyzePedestals = 1; // Analyze pedestals first normal events
- return(0);
+ fAnalyzePedestals = 1; // Analyze pedestals first normal events
+ return (0);
}
- if(fAnalyzePedestals) {
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if (fAnalyzePedestals) {
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
fPlanes[ip]->CalculatePedestals();
}
- fAnalyzePedestals = 0; // Don't analyze pedestals next event
+ fAnalyzePedestals = 0; // Don't analyze pedestals next event
}
// Let each plane get its hits
Int_t nexthit = 0;
- fNfptimes=0;
+ fNfptimes = 0;
Int_t thits = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
- fPlaneCenter[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
+ fPlaneCenter[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
fPlaneSpacing[ip] = fPlanes[ip]->GetSpacing();
// nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
// GN: select only events that have reasonable TDC values to start with
// as per the Engine h_strip_scin.f
- nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
- thits+=fPlanes[ip]->GetNScinHits();
+ nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
+ thits += fPlanes[ip]->GetNScinHits();
}
- fStartTime=-1000;
- if (thits>0 ) EstimateFocalPlaneTime();
+ fStartTime = -1000;
+ if (thits > 0)
+ EstimateFocalPlaneTime();
if (fdebugprintscinraw == 1) {
- for(UInt_t ihit = 0; ihit < fNRawHits ; ihit++) {
-// THcRawHodoHit* hit = (THcRawHodoHit *) fRawHitList->At(ihit);
-// cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
-// << hit->fADC_pos << " " << hit->fADC_neg << " " << hit->fTDC_pos
-// << " " << hit->fTDC_neg << endl;
+ for (UInt_t ihit = 0; ihit < fNRawHits; ihit++) {
+ // THcRawHodoHit* hit = (THcRawHodoHit *) fRawHitList->At(ihit);
+ // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+ // << hit->fADC_pos << " " << hit->fADC_neg << " " << hit->fTDC_pos
+ // << " " << hit->fTDC_neg << endl;
}
cout << endl;
}
-
return fNHits;
}
//_____________________________________________________________________________
-void THcHodoscope::EstimateFocalPlaneTime()
-{
- /*! \brief Calculates the Drift Chamber start time and fBetaNoTrk (velocity determined without track info)
+void THcHodoscope::EstimateFocalPlaneTime() {
+ /*! \brief Calculates the Drift Chamber start time and fBetaNoTrk (velocity determined without
+ * track info)
*
* - Called by THcHodoscope::Decode
* - selects good scintillator paddle hits
- * + loops through hits in each scintillator plane and fills histogram array, "timehist", with corrected times for positive
- * and negative ends of each paddle
+ * + loops through hits in each scintillator plane and fills histogram array, "timehist", with
+ * corrected times for positive and negative ends of each paddle
* + Determines the peak of "timehist"
*
*/
- Int_t ihit=0;
- Int_t nscinhits=0; // Total # hits with at least one good tdc
+ Int_t ihit = 0;
+ Int_t nscinhits = 0; // Total # hits with at least one good tdc
hTime->Reset();
//
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- Int_t nphits=fPlanes[ip]->GetNScinHits();
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ Int_t nphits = fPlanes[ip]->GetNScinHits();
nscinhits += nphits;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
- for(Int_t i=0;i<nphits;i++) {
- THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
- if(hit->GetHasCorrectedTimes()) {
- Double_t postime=hit->GetPosTOFCorrectedTime();
- Double_t negtime=hit->GetNegTOFCorrectedTime();
- hTime->Fill(postime);
- hTime->Fill(negtime);
+ for (Int_t i = 0; i < nphits; i++) {
+ THcHodoHit* hit = (THcHodoHit*)hodoHits->At(i);
+ if (hit->GetHasCorrectedTimes()) {
+ Double_t postime = hit->GetPosTOFCorrectedTime();
+ Double_t negtime = hit->GetNegTOFCorrectedTime();
+ hTime->Fill(postime);
+ hTime->Fill(negtime);
}
}
}
//
//
- ihit = 0;
- Double_t fpTimeSum = 0.0;
- fNfptimes=0;
- Int_t Ngood_hits_plane=0;
- Double_t Plane_fptime_sum=0.0;
- Bool_t goodplanetime[fNPlanes];
- Bool_t twogoodtimes[nscinhits];
- Double_t tmin = 0.5*hTime->GetMaximumBin();
- fTimeHist_Peak= tmin;
- fTimeHist_Sigma= hTime->GetRMS();
- fTimeHist_Hits= hTime->Integral();
- _basic_data.fTimeHist_Peak = fTimeHist_Peak ;
+ ihit = 0;
+ Double_t fpTimeSum = 0.0;
+ fNfptimes = 0;
+ Int_t Ngood_hits_plane = 0;
+ Double_t Plane_fptime_sum = 0.0;
+ Bool_t goodplanetime[fNPlanes];
+ Bool_t twogoodtimes[nscinhits];
+ Double_t tmin = 0.5 * hTime->GetMaximumBin();
+ fTimeHist_Peak = tmin;
+ fTimeHist_Sigma = hTime->GetRMS();
+ fTimeHist_Hits = hTime->Integral();
+ _basic_data.fTimeHist_Peak = fTimeHist_Peak;
_basic_data.fTimeHist_Sigma = fTimeHist_Sigma;
- _basic_data.fTimeHist_Hits = fTimeHist_Hits ;
- for(Int_t ip=0;ip<fNumPlanesBetaCalc;ip++) {
- goodplanetime[ip] = kFALSE;
- Int_t nphits=fPlanes[ip]->GetNScinHits();
+ _basic_data.fTimeHist_Hits = fTimeHist_Hits;
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ goodplanetime[ip] = kFALSE;
+ Int_t nphits = fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
- Ngood_hits_plane=0;
- Plane_fptime_sum=0.0;
- for(Int_t i=0;i<nphits;i++) {
- THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
+ Ngood_hits_plane = 0;
+ Plane_fptime_sum = 0.0;
+ for (Int_t i = 0; i < nphits; i++) {
+ THcHodoHit* hit = (THcHodoHit*)hodoHits->At(i);
twogoodtimes[ihit] = kFALSE;
- if(hit->GetHasCorrectedTimes()) {
- Double_t postime=hit->GetPosTOFCorrectedTime();
- Double_t negtime=hit->GetNegTOFCorrectedTime();
- if ((postime>(tmin-fTofTolerance)) && (postime<(tmin+fTofTolerance)) &&
- (negtime>(tmin-fTofTolerance)) && (negtime<(tmin+fTofTolerance)) ) {
- hit->SetTwoGoodTimes(kTRUE);
- twogoodtimes[ihit] = kTRUE; // Both tubes fired
- Int_t index=hit->GetPaddleNumber()-1; //
- Double_t fptime;
- if(fCosmicFlag==1) {
- fptime = hit->GetScinCorrectedTime()
- + (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
- / (29.979 * fBetaNominal);
- }else{
- fptime = hit->GetScinCorrectedTime()
- - (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
- / (29.979 * fBetaNominal);
- }
+ if (hit->GetHasCorrectedTimes()) {
+ Double_t postime = hit->GetPosTOFCorrectedTime();
+ Double_t negtime = hit->GetNegTOFCorrectedTime();
+ if ((postime > (tmin - fTofTolerance)) && (postime < (tmin + fTofTolerance)) &&
+ (negtime > (tmin - fTofTolerance)) && (negtime < (tmin + fTofTolerance))) {
+ hit->SetTwoGoodTimes(kTRUE);
+ twogoodtimes[ihit] = kTRUE; // Both tubes fired
+ Int_t index = hit->GetPaddleNumber() - 1; //
+ Double_t fptime;
+ if (fCosmicFlag == 1) {
+ fptime = hit->GetScinCorrectedTime() +
+ (fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
+ (29.979 * fBetaNominal);
+ } else {
+ fptime = hit->GetScinCorrectedTime() -
+ (fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
+ (29.979 * fBetaNominal);
+ }
Ngood_hits_plane++;
- Plane_fptime_sum+=fptime;
- fpTimeSum += fptime;
- fNfptimes++;
- goodplanetime[ip] = kTRUE;
- } else {
- hit->SetTwoGoodTimes(kFALSE);
- }
+ Plane_fptime_sum += fptime;
+ fpTimeSum += fptime;
+ fNfptimes++;
+ goodplanetime[ip] = kTRUE;
+ } else {
+ hit->SetTwoGoodTimes(kFALSE);
+ }
}
ihit++;
}
- if (Ngood_hits_plane) fPlanes[ip]->SetFpTime(Plane_fptime_sum/float(Ngood_hits_plane));
+ if (Ngood_hits_plane)
+ fPlanes[ip]->SetFpTime(Plane_fptime_sum / float(Ngood_hits_plane));
fPlanes[ip]->SetNGoodHits(Ngood_hits_plane);
}
- if(fNfptimes>0) {
- fStartTime = fpTimeSum/fNfptimes;
- fGoodStartTime=kTRUE;
- fFPTimeAll = fStartTime ;
+ if (fNfptimes > 0) {
+ fStartTime = fpTimeSum / fNfptimes;
+ fGoodStartTime = kTRUE;
+ fFPTimeAll = fStartTime;
} else {
- fStartTime = fStartTimeCenter;
- fGoodStartTime=kFALSE;
- fFPTimeAll = fStartTime ;
+ fStartTime = fStartTimeCenter;
+ fGoodStartTime = kFALSE;
+ fFPTimeAll = fStartTime;
}
- //
- //
+ //
+ //
hTime->Reset();
//
- if((goodplanetime[0]||goodplanetime[1]) &&(goodplanetime[2]||goodplanetime[3])) {
+ if ((goodplanetime[0] || goodplanetime[1]) && (goodplanetime[2] || goodplanetime[3])) {
- Double_t sumW = 0.;
- Double_t sumT = 0.;
- Double_t sumZ = 0.;
+ Double_t sumW = 0.;
+ Double_t sumT = 0.;
+ Double_t sumZ = 0.;
Double_t sumZZ = 0.;
Double_t sumTZ = 0.;
- Int_t ihhit = 0;
+ Int_t ihhit = 0;
- for(Int_t ip=0;ip<fNumPlanesBetaCalc;ip++) {
- Int_t nphits=fPlanes[ip]->GetNScinHits();
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ Int_t nphits = fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
-
- for(Int_t i=0;i<nphits;i++) {
- Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
-
- if(twogoodtimes[ihhit]){
-
- Double_t sigma = 0.0;
- if(fTofUsingInvAdc)
- {
- sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
- TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
- }
- else{
- sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) +
- TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) );
- }
-
- Double_t scinWeight = 1 / TMath::Power(sigma,2);
- Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
- // cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " << sigma << endl;
- //cout << "fHodoSigma+ = " << fHodoSigmaPos[GetScinIndex(ip,index)] << endl;
- sumW += scinWeight;
- sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
- sumZ += scinWeight * zPosition;
- sumZZ += scinWeight * ( zPosition * zPosition );
- sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
-
- } // condition of good scin time
- ihhit ++;
+
+ for (Int_t i = 0; i < nphits; i++) {
+ Int_t index = ((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber() - 1;
+
+ if (twogoodtimes[ihhit]) {
+
+ Double_t sigma = 0.0;
+ if (fTofUsingInvAdc) {
+ sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoPosSigma[GetScinIndex(ip, index)], 2) +
+ TMath::Power(fHodoNegSigma[GetScinIndex(ip, index)], 2)));
+ } else {
+ sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoSigmaPos[GetScinIndex(ip, index)], 2) +
+ TMath::Power(fHodoSigmaNeg[GetScinIndex(ip, index)], 2)));
+ }
+
+ Double_t scinWeight = 1 / TMath::Power(sigma, 2);
+ Double_t zPosition = fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos();
+ // cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " <<
+ // sigma << endl; cout << "fHodoSigma+ = " << fHodoSigmaPos[GetScinIndex(ip,index)] <<
+ // endl;
+ sumW += scinWeight;
+ sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
+ sumZ += scinWeight * zPosition;
+ sumZZ += scinWeight * (zPosition * zPosition);
+ sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
+
+ } // condition of good scin time
+ ihhit++;
} // loop over hits of plane
- } // loop over planes
+ } // loop over planes
- Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
- Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
+ Double_t tmp = sumW * sumZZ - sumZ * sumZ;
+ Double_t t0 = (sumT * sumZZ - sumZ * sumTZ) / tmp;
Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
- if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
+ if (TMath::Abs(tmpDenom) > (1 / 10000000000.0)) {
- fBetaNoTrk = tmp / tmpDenom;
+ fBetaNoTrk = tmp / tmpDenom;
fBetaNoTrkChiSq = 0.;
- ihhit = 0;
+ ihhit = 0;
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){ // Loop over planes
- Int_t nphits=fPlanes[ip]->GetNScinHits();
- TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) { // Loop over planes
+ Int_t nphits = fPlanes[ip]->GetNScinHits();
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
- for(Int_t i=0;i<nphits;i++) {
- Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
+ for (Int_t i = 0; i < nphits; i++) {
+ Int_t index = ((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber() - 1;
- if(twogoodtimes[ihhit]) {
+ if (twogoodtimes[ihhit]) {
- Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
- Double_t timeDif = ( ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0 );
-
- Double_t sigma = 0.0;
- if(fTofUsingInvAdc){
- sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
- TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
- }
- else {
- sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) +
- TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) );
- }
+ Double_t zPosition = fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos();
+ Double_t timeDif = (((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0);
- fBetaNoTrkChiSq += ( ( zPosition / fBetaNoTrk - timeDif ) *
- ( zPosition / fBetaNoTrk - timeDif ) ) / ( sigma * sigma );
+ Double_t sigma = 0.0;
+ if (fTofUsingInvAdc) {
+ sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoPosSigma[GetScinIndex(ip, index)], 2) +
+ TMath::Power(fHodoNegSigma[GetScinIndex(ip, index)], 2)));
+ } else {
+ sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoSigmaPos[GetScinIndex(ip, index)], 2) +
+ TMath::Power(fHodoSigmaNeg[GetScinIndex(ip, index)], 2)));
+ }
+ fBetaNoTrkChiSq +=
+ ((zPosition / fBetaNoTrk - timeDif) * (zPosition / fBetaNoTrk - timeDif)) /
+ (sigma * sigma);
- } // condition for good scin time
- ihhit++;
- } // loop over hits of a plane
- } // loop over planes
+ } // condition for good scin time
+ ihhit++;
+ } // loop over hits of a plane
+ } // loop over planes
Double_t pathNorm = 1.0;
fBetaNoTrk = fBetaNoTrk * pathNorm;
- fBetaNoTrk = fBetaNoTrk / 29.979; // velocity / c
+ fBetaNoTrk = fBetaNoTrk / 29.979; // velocity / c
- } // condition for fTmpDenom
+ } // condition for fTmpDenom
else {
- fBetaNoTrk = 0.;
+ fBetaNoTrk = 0.;
fBetaNoTrkChiSq = -2.;
} // else condition for fTmpDenom
//
- fGoodEventTOFCalib=kFALSE;
- if ((fNumPlanesBetaCalc==4)&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&goodplanetime[3]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1&&fPlanes[3]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
- if ((fNumPlanesBetaCalc==3)&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
+ fGoodEventTOFCalib = kFALSE;
+ if ((fNumPlanesBetaCalc == 4) && goodplanetime[0] && goodplanetime[1] && goodplanetime[2] &&
+ goodplanetime[3] && fPlanes[0]->GetNGoodHits() == 1 && fPlanes[1]->GetNGoodHits() == 1 &&
+ fPlanes[2]->GetNGoodHits() == 1 && fPlanes[3]->GetNGoodHits() == 1)
+ fGoodEventTOFCalib = kTRUE;
+ if ((fNumPlanesBetaCalc == 3) && goodplanetime[0] && goodplanetime[1] && goodplanetime[2] &&
+ fPlanes[0]->GetNGoodHits() == 1 && fPlanes[1]->GetNGoodHits() == 1 &&
+ fPlanes[2]->GetNGoodHits() == 1)
+ fGoodEventTOFCalib = kTRUE;
//
//
}
}
//_____________________________________________________________________________
-Int_t THcHodoscope::ApplyCorrections( void )
-{
- return(0);
-}
+Int_t THcHodoscope::ApplyCorrections(void) { return (0); }
//_____________________________________________________________________________
-Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle,
- const ESide side)
-{
- return(0.0);
-}
-
+Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle, const ESide side) { return (0.0); }
//_____________________________________________________________________________
-Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
-{
+Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
-
- Int_t ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+ Int_t ntracks = tracks.GetLast() + 1; // Number of reconstructed tracks
// -------------------------------------------------
// fDumpOut << " ntrack = " << ntracks << endl;
- if (ntracks > 0 ) {
+ if (ntracks > 0) {
// **MAIN LOOP: Loop over all tracks and get corrected time, tof, beta...
vector<Double_t> nPmtHit(ntracks);
vector<Double_t> timeAtFP(ntracks);
fdEdX.reserve(ntracks);
fGoodFlags.reserve(ntracks);
- for ( Int_t itrack = 0; itrack < ntracks; itrack++ ) { // Line 133
- nPmtHit[itrack]=0;
- timeAtFP[itrack]=0;
-
- THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
- if (!theTrack) return -1;
-
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
- fGoodPlaneTime[ip] = kFALSE;
- fNScinHits[ip] = 0;
- fNPlaneTime[ip] = 0;
- fSumPlaneTime[ip] = 0.;
+ for (Int_t itrack = 0; itrack < ntracks; itrack++) { // Line 133
+ nPmtHit[itrack] = 0;
+ timeAtFP[itrack] = 0;
+
+ THaTrack* theTrack = dynamic_cast<THaTrack*>(tracks.At(itrack));
+ if (!theTrack)
+ return -1;
+
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ fGoodPlaneTime[ip] = kFALSE;
+ fNScinHits[ip] = 0;
+ fNPlaneTime[ip] = 0;
+ fSumPlaneTime[ip] = 0.;
}
- std::vector<Double_t> dedx_temp;
- std::vector<std::vector<GoodFlags> > goodflagstmp1;
+ std::vector<Double_t> dedx_temp;
+ std::vector<std::vector<GoodFlags>> goodflagstmp1;
goodflagstmp1.reserve(fNumPlanesBetaCalc);
#if __cplusplus >= 201103L
fdEdX.push_back(std::move(dedx_temp)); // Create array of dedx per hit
@@ -1038,9 +1073,9 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
fdEdX.push_back(dedx_temp); // Create array of dedx per hit
fGoodFlags.push_back(goodflagstmp1);
#endif
- Int_t nFPTime = 0;
+ Int_t nFPTime = 0;
Double_t betaChiSq = -3;
- Double_t beta = 0;
+ Double_t beta = 0;
// timeAtFP[itrack] = 0.;
Double_t sumFPTime = 0.; // Line 138
fNScinHit.push_back(0);
@@ -1057,317 +1092,321 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
//! to accomodate difference in TOF for other particles
//! Default value in case user hasnt defined something reasonable
-
// Loop over scintillator planes.
// In ENGINE, its loop over good scintillator hits.
hTime->Reset();
- fTOFCalc.clear(); // SAW - Can we
- fTOFPInfo.clear(); // SAW - combine these two?
- Int_t ihhit = 0; // Hit # overall
+ fTOFCalc.clear(); // SAW - Can we
+ fTOFPInfo.clear(); // SAW - combine these two?
+ Int_t ihhit = 0; // Hit # overall
- for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
- std::vector<GoodFlags> goodflagstmp2;
- goodflagstmp2.reserve(fNScinHits[ip]);
+ std::vector<GoodFlags> goodflagstmp2;
+ goodflagstmp2.reserve(fNScinHits[ip]);
#if __cplusplus >= 201103L
- fGoodFlags[itrack].push_back(std::move(goodflagstmp2));
+ fGoodFlags[itrack].push_back(std::move(goodflagstmp2));
#else
- fGoodFlags[itrack].push_back(goodflagstmp2);
+ fGoodFlags[itrack].push_back(goodflagstmp2);
#endif
- fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
- TClonesArray* hodoHits = fPlanes[ip]->GetHits();
-
- Double_t zPos = fPlanes[ip]->GetZpos();
- Double_t dzPos = fPlanes[ip]->GetDzpos();
-
- // first loop over hits with in a single plane
- for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
- // iphit is hit # within a plane
- THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
-
- fTOFPInfo.push_back(TOFPInfo());
- // Can remove these as we will initialize in the constructor
- // fTOFPInfo[ihhit].time_pos = -99.0;
- // fTOFPInfo[ihhit].time_neg = -99.0;
- // fTOFPInfo[ihhit].keep_pos = kFALSE;
- // fTOFPInfo[ihhit].keep_neg = kFALSE;
- fTOFPInfo[ihhit].scin_pos_time = 0.0;
- fTOFPInfo[ihhit].scin_neg_time = 0.0;
- fTOFPInfo[ihhit].hit = hit;
- fTOFPInfo[ihhit].planeIndex = ip;
- fTOFPInfo[ihhit].hitNumInPlane = iphit;
- fTOFPInfo[ihhit].onTrack = kFALSE;
-
- Int_t paddle = hit->GetPaddleNumber()-1;
- Double_t zposition = zPos + (paddle%2)*dzPos;
-
- Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
- ( zposition ); // Line 183
-
- Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
- ( zposition ); // Line 184
-
- Double_t scinTrnsCoord, scinLongCoord;
- if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 185
- scinTrnsCoord = xHitCoord;
- scinLongCoord = yHitCoord;
- } else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 188
- scinTrnsCoord = yHitCoord;
- scinLongCoord = xHitCoord;
- } else { return -1; } // Line 195
-
- fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
- fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;
-
- Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
-
- // Index to access the 2d arrays of paddle/scintillator properties
- Int_t fPIndex = GetScinIndex(ip,paddle);
-
- if ( TMath::Abs( scinCenter - scinTrnsCoord ) <
- ( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
-
- fTOFPInfo[ihhit].onTrack = kTRUE;
- Double_t zcor = zposition/(29.979*fBetaNominal)*
- TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
- + theTrack->GetPhi()*theTrack->GetPhi());
- fTOFPInfo[ihhit].zcor = zcor;
- if (fCosmicFlag) {
- Double_t zcor = -zposition/(29.979*1.0)*
- TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
- + theTrack->GetPhi()*theTrack->GetPhi());
- fTOFPInfo[ihhit].zcor = zcor;
- }
- Double_t tdc_pos = hit->GetPosTDC();
- if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) {
- Double_t adc_pos = hit->GetPosADC();
- Double_t adcamp_pos = hit->GetPosADCpeak();
- Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
- fTOFPInfo[ihhit].pathp = pathp;
- Double_t timep = tdc_pos*fScinTdcToTime;
- if(fTofUsingInvAdc) {
- timep -= fHodoPosInvAdcOffset[fPIndex]
- + pathp/fHodoPosInvAdcLinear[fPIndex]
- + fHodoPosInvAdcAdc[fPIndex]
- /TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
- } else {
- //Double_t tw_corr_pos = fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
- Double_t tw_corr_pos=0.;
- pathp=scinLongCoord;
- if (adcamp_pos>0) tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
- timep += -tw_corr_pos + fHodo_LCoeff[fPIndex]+ pathp/fHodoVelFit[fPIndex];
- }
- fTOFPInfo[ihhit].scin_pos_time = timep;
- timep -= zcor;
- fTOFPInfo[ihhit].time_pos = timep;
+ fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+
+ Double_t zPos = fPlanes[ip]->GetZpos();
+ Double_t dzPos = fPlanes[ip]->GetDzpos();
+
+ // first loop over hits with in a single plane
+ for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++) {
+ // iphit is hit # within a plane
+ THcHodoHit* hit = (THcHodoHit*)hodoHits->At(iphit);
+
+ fTOFPInfo.push_back(TOFPInfo());
+ // Can remove these as we will initialize in the constructor
+ // fTOFPInfo[ihhit].time_pos = -99.0;
+ // fTOFPInfo[ihhit].time_neg = -99.0;
+ // fTOFPInfo[ihhit].keep_pos = kFALSE;
+ // fTOFPInfo[ihhit].keep_neg = kFALSE;
+ fTOFPInfo[ihhit].scin_pos_time = 0.0;
+ fTOFPInfo[ihhit].scin_neg_time = 0.0;
+ fTOFPInfo[ihhit].hit = hit;
+ fTOFPInfo[ihhit].planeIndex = ip;
+ fTOFPInfo[ihhit].hitNumInPlane = iphit;
+ fTOFPInfo[ihhit].onTrack = kFALSE;
+
+ Int_t paddle = hit->GetPaddleNumber() - 1;
+ Double_t zposition = zPos + (paddle % 2) * dzPos;
+
+ Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() * (zposition); // Line 183
+
+ Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() * (zposition); // Line 184
+
+ Double_t scinTrnsCoord, scinLongCoord;
+ if ((ip == 0) || (ip == 2)) { // !x plane. Line 185
+ scinTrnsCoord = xHitCoord;
+ scinLongCoord = yHitCoord;
+ } else if ((ip == 1) || (ip == 3)) { // !y plane. Line 188
+ scinTrnsCoord = yHitCoord;
+ scinLongCoord = xHitCoord;
+ } else {
+ return -1;
+ } // Line 195
+
+ fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
+ fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;
+
+ Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
+
+ // Index to access the 2d arrays of paddle/scintillator properties
+ Int_t fPIndex = GetScinIndex(ip, paddle);
+ Double_t betatrack = theTrack->GetP() / TMath::Sqrt(theTrack->GetP() * theTrack->GetP() +
+ fPartMass * fPartMass);
+
+ if (TMath::Abs(scinCenter - scinTrnsCoord) <
+ (fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop())) { // Line 293
+
+ fTOFPInfo[ihhit].onTrack = kTRUE;
+ Double_t zcor = zposition / (29.979 * betatrack) *
+ TMath::Sqrt(1. + theTrack->GetTheta() * theTrack->GetTheta() +
+ theTrack->GetPhi() * theTrack->GetPhi());
+ fTOFPInfo[ihhit].zcor = zcor;
+ if (fCosmicFlag) {
+ Double_t zcor = -zposition / (29.979 * 1.0) *
+ TMath::Sqrt(1. + theTrack->GetTheta() * theTrack->GetTheta() +
+ theTrack->GetPhi() * theTrack->GetPhi());
+ fTOFPInfo[ihhit].zcor = zcor;
+ }
+ Double_t tdc_pos = hit->GetPosTDC();
+ if (tdc_pos >= fScinTdcMin && tdc_pos <= fScinTdcMax) {
+ Double_t adc_pos = hit->GetPosADC();
+ Double_t adcamp_pos = hit->GetPosADCpeak();
+ Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
+ fTOFPInfo[ihhit].pathp = pathp;
+ Double_t timep = tdc_pos * fScinTdcToTime;
+ if (fTofUsingInvAdc) {
+ timep -= fHodoPosInvAdcOffset[fPIndex] + pathp / fHodoPosInvAdcLinear[fPIndex] +
+ fHodoPosInvAdcAdc[fPIndex] / TMath::Sqrt(TMath::Max(20.0 * .020, adc_pos));
+ } else {
+ // Double_t tw_corr_pos =
+ // fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) -
+ // fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
+ Double_t tw_corr_pos = 0.;
+ pathp = scinLongCoord;
+ if (adcamp_pos > 0)
+ tw_corr_pos = 1. / pow(adcamp_pos / fTdc_Thrs, fHodoPos_c2[fPIndex]) -
+ 1. / pow(200. / fTdc_Thrs, fHodoPos_c2[fPIndex]);
+ timep += -tw_corr_pos + fHodo_LCoeff[fPIndex] + pathp / fHodoVelFit[fPIndex];
+ }
+ fTOFPInfo[ihhit].scin_pos_time = timep;
+ timep -= zcor;
+ fTOFPInfo[ihhit].time_pos = timep;
hTime->Fill(timep);
- }
- Double_t tdc_neg = hit->GetNegTDC();
- if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) {
- Double_t adc_neg = hit->GetNegADC();
- Double_t adcamp_neg = hit->GetNegADCpeak();
- Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
- fTOFPInfo[ihhit].pathn = pathn;
- Double_t timen = tdc_neg*fScinTdcToTime;
- if(fTofUsingInvAdc) {
- timen -= fHodoNegInvAdcOffset[fPIndex]
- + pathn/fHodoNegInvAdcLinear[fPIndex]
- + fHodoNegInvAdcAdc[fPIndex]
- /TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
- } else {
- pathn=scinLongCoord ;
- Double_t tw_corr_neg =0 ;
- if (adcamp_neg >0) tw_corr_neg= 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]);
- timen += -tw_corr_neg- 2*fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex]- pathn/fHodoVelFit[fPIndex];
-
- }
- fTOFPInfo[ihhit].scin_neg_time = timen;
- timen -= zcor;
- fTOFPInfo[ihhit].time_neg = timen;
+ }
+ Double_t tdc_neg = hit->GetNegTDC();
+ if (tdc_neg >= fScinTdcMin && tdc_neg <= fScinTdcMax) {
+ Double_t adc_neg = hit->GetNegADC();
+ Double_t adcamp_neg = hit->GetNegADCpeak();
+ Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
+ fTOFPInfo[ihhit].pathn = pathn;
+ Double_t timen = tdc_neg * fScinTdcToTime;
+ if (fTofUsingInvAdc) {
+ timen -= fHodoNegInvAdcOffset[fPIndex] + pathn / fHodoNegInvAdcLinear[fPIndex] +
+ fHodoNegInvAdcAdc[fPIndex] / TMath::Sqrt(TMath::Max(20.0 * .020, adc_neg));
+ } else {
+ pathn = scinLongCoord;
+ Double_t tw_corr_neg = 0;
+ if (adcamp_neg > 0)
+ tw_corr_neg = 1. / pow(adcamp_neg / fTdc_Thrs, fHodoNeg_c2[fPIndex]) -
+ 1. / pow(200. / fTdc_Thrs, fHodoNeg_c2[fPIndex]);
+ timen += -tw_corr_neg - 2 * fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex] -
+ pathn / fHodoVelFit[fPIndex];
+ }
+ fTOFPInfo[ihhit].scin_neg_time = timen;
+ timen -= zcor;
+ fTOFPInfo[ihhit].time_neg = timen;
hTime->Fill(timen);
- }
- } // condition for cenetr on a paddle
- ihhit++;
- } // First loop over hits in a plane <---------
-
- //-----------------------------------------------------------------------------------------------
- //------------- First large loop over scintillator hits ends here --------------------
- //-----------------------------------------------------------------------------------------------
+ }
+ } // condition for cenetr on a paddle
+ ihhit++;
+ } // First loop over hits in a plane <---------
+
+ //-----------------------------------------------------------------------------------------------
+ //------------- First large loop over scintillator hits ends here --------------------
+ //-----------------------------------------------------------------------------------------------
}
- Int_t nhits=ihhit;
-
-
- if(0.5*hTime->GetMaximumBin() > 0) {
- Double_t tmin = 0.5*hTime->GetMaximumBin();
-
- for(Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
- if ( ( fTOFPInfo[ih].time_pos > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_pos < ( tmin + fTofTolerance ) ) ) {
- fTOFPInfo[ih].keep_pos=kTRUE;
- }
- if ( ( fTOFPInfo[ih].time_neg > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_neg < ( tmin + fTofTolerance ) ) ){
- fTOFPInfo[ih].keep_neg=kTRUE;
- }
- }
+ Int_t nhits = ihhit;
+
+ if (0.5 * hTime->GetMaximumBin() > 0) {
+ Double_t tmin = 0.5 * hTime->GetMaximumBin();
+
+ for (Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
+ if ((fTOFPInfo[ih].time_pos > (tmin - fTofTolerance)) &&
+ (fTOFPInfo[ih].time_pos < (tmin + fTofTolerance))) {
+ fTOFPInfo[ih].keep_pos = kTRUE;
+ }
+ if ((fTOFPInfo[ih].time_neg > (tmin - fTofTolerance)) &&
+ (fTOFPInfo[ih].time_neg < (tmin + fTofTolerance))) {
+ fTOFPInfo[ih].keep_neg = kTRUE;
+ }
+ }
}
//---------------------------------------------------------------------------------------------
- // ---------------------- Second loop over scint. hits in a plane -----------------------------
+ // ---------------------- Second loop over scint. hits in a plane
+ // -----------------------------
//---------------------------------------------------------------------------------------------
-
+
fdEdX[itrack].reserve(nhits);
fTOFCalc.reserve(nhits);
- for(Int_t ih=0; ih < nhits; ih++) {
- THcHodoHit *hit = fTOFPInfo[ih].hit;
- Int_t iphit = fTOFPInfo[ih].hitNumInPlane;
- Int_t ip = fTOFPInfo[ih].planeIndex;
- // fDumpOut << " looping over hits = " << ih << " plane = " << ip+1 << endl;
- // Flags are used by THcHodoEff
- fGoodFlags[itrack][ip].reserve(nhits);
- fGoodFlags[itrack][ip].push_back(GoodFlags());
- assert( iphit >= 0 && (size_t)iphit < fGoodFlags[itrack][ip].size() );
- fGoodFlags[itrack][ip][iphit].onTrack = kFALSE;
- fGoodFlags[itrack][ip][iphit].goodScinTime = kFALSE;
- fGoodFlags[itrack][ip][iphit].goodTdcNeg = kFALSE;
- fGoodFlags[itrack][ip][iphit].goodTdcPos = kFALSE;
-
- fTOFCalc.push_back(TOFCalc());
- // Do we set back to false for each track, or just once per event?
- assert( ih >= 0 && (size_t)ih < fTOFCalc.size() );
- fTOFCalc[ih].good_scin_time = kFALSE;
- // These need a track index too to calculate efficiencies
- fTOFCalc[ih].good_tdc_pos = kFALSE;
- fTOFCalc[ih].good_tdc_neg = kFALSE;
- fTOFCalc[ih].pindex = ip;
-
- Int_t paddle = hit->GetPaddleNumber()-1;
- fTOFCalc[ih].hit_paddle = paddle;
- fTOFCalc[ih].good_raw_pad = paddle;
-
- // Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
- // Double_t scinTrnsCoord = fTOFPInfo[ih].scinTrnsCoord;
- // Double_t scinLongCoord = fTOFPInfo[ih].scinLongCoord;
-
- Int_t fPIndex = GetScinIndex(ip,paddle);
-
- if (fTOFPInfo[ih].onTrack) {
- fGoodFlags[itrack][ip][iphit].onTrack = kTRUE;
- if ( fTOFPInfo[ih].keep_pos ) { // 301
- fTOFCalc[ih].good_tdc_pos = kTRUE;
- fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
- }
- if ( fTOFPInfo[ih].keep_neg ) { //
- fTOFCalc[ih].good_tdc_neg = kTRUE;
- fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
- }
- // ** Calculate ave time for scin and error.
- if ( fTOFCalc[ih].good_tdc_pos ){
- if ( fTOFCalc[ih].good_tdc_neg ){
- fTOFCalc[ih].scin_time = ( fTOFPInfo[ih].scin_pos_time +
- fTOFPInfo[ih].scin_neg_time ) / 2.;
- fTOFCalc[ih].scin_time_fp = ( fTOFPInfo[ih].time_pos +
- fTOFPInfo[ih].time_neg ) / 2.;
-
- if (fTofUsingInvAdc){
- fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
- fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.;
- }
- else {
- fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoSigmaPos[fPIndex] * fHodoSigmaPos[fPIndex] +
- fHodoSigmaNeg[fPIndex] * fHodoSigmaNeg[fPIndex] )/2.;
- }
-
- fTOFCalc[ih].good_scin_time = kTRUE;
- fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
- } else{
- fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_pos_time;
- fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_pos;
-
- if (fTofUsingInvAdc){
- fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex];
- }
- else{
- fTOFCalc[ih].scin_sigma = fHodoSigmaPos[fPIndex];
- }
-
- fTOFCalc[ih].good_scin_time = kTRUE;
- fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
- }
- } else {
- if ( fTOFCalc[ih].good_tdc_neg ){
- fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_neg_time;
- fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_neg;
- if (fTofUsingInvAdc){
- fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex];
- }
- else{
- fTOFCalc[ih].scin_sigma = fHodoSigmaNeg[fPIndex];
- }
- fTOFCalc[ih].good_scin_time = kTRUE;
- fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
- }
- } // In h_tof.f this includes the following if condition for time at focal plane
- // // because it is written in FORTRAN code
-
- // c Get time at focal plane
- if ( fTOFCalc[ih].good_scin_time ){
-
- // scin_time_fp doesn't need to be an array
- // Is this any different than the average of time_pos and time_neg?
- // Double_t scin_time_fp = ( fTOFPInfo[ih].time_pos +
- // fTOFPInfo[ih].time_neg ) / 2.;
- Double_t scin_time_fp = fTOFCalc[ih].scin_time_fp;
-
- sumFPTime = sumFPTime + scin_time_fp;
- nFPTime ++;
-
- fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
- fNPlaneTime[ip] ++;
- fNScinHit[itrack] ++;
-
- if ( ( fTOFCalc[ih].good_tdc_pos ) && ( fTOFCalc[ih].good_tdc_neg ) ){
- nPmtHit[itrack] = nPmtHit[itrack] + 2;
- } else {
- nPmtHit[itrack] = nPmtHit[itrack] + 1;
- }
-
- fdEdX[itrack].push_back(0.0);
- assert( fNScinHit[itrack] > 0 && (size_t)fNScinHit[itrack] < fdEdX[itrack].size()+1 );
-
- // --------------------------------------------------------------------------------------------
- if ( fTOFCalc[ih].good_tdc_pos ){
- if ( fTOFCalc[ih].good_tdc_neg ){
- fdEdX[itrack][fNScinHit[itrack]-1]=
- TMath::Sqrt( TMath::Max( 0., hit->GetPosADC() * hit->GetNegADC() ) );
- } else{
- fdEdX[itrack][fNScinHit[itrack]-1]=
- TMath::Max( 0., hit->GetPosADC() );
- }
- } else{
- if ( fTOFCalc[ih].good_tdc_neg ){
- fdEdX[itrack][fNScinHit[itrack]-1]=
- TMath::Max( 0., hit->GetNegADC() );
- } else{
- fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
- }
- }
- // --------------------------------------------------------------------------------------------
-
-
- } // time at focal plane condition
- } // on track condition
-
- // ** See if there are any good time measurements in the plane.
- if ( fTOFCalc[ih].good_scin_time ){
- fGoodPlaneTime[ip] = kTRUE;
- fTOFCalc[ih].dedx = fdEdX[itrack][fNScinHit[itrack]-1];
- } else {
- fTOFCalc[ih].dedx = 0.0;
- }
-
+ for (Int_t ih = 0; ih < nhits; ih++) {
+ THcHodoHit* hit = fTOFPInfo[ih].hit;
+ Int_t iphit = fTOFPInfo[ih].hitNumInPlane;
+ Int_t ip = fTOFPInfo[ih].planeIndex;
+ // fDumpOut << " looping over hits = " << ih << " plane = " << ip+1 << endl;
+ // Flags are used by THcHodoEff
+ fGoodFlags[itrack][ip].reserve(nhits);
+ fGoodFlags[itrack][ip].push_back(GoodFlags());
+ assert(iphit >= 0 && (size_t)iphit < fGoodFlags[itrack][ip].size());
+ fGoodFlags[itrack][ip][iphit].onTrack = kFALSE;
+ fGoodFlags[itrack][ip][iphit].goodScinTime = kFALSE;
+ fGoodFlags[itrack][ip][iphit].goodTdcNeg = kFALSE;
+ fGoodFlags[itrack][ip][iphit].goodTdcPos = kFALSE;
+
+ fTOFCalc.push_back(TOFCalc());
+ // Do we set back to false for each track, or just once per event?
+ assert(ih >= 0 && (size_t)ih < fTOFCalc.size());
+ fTOFCalc[ih].good_scin_time = kFALSE;
+ // These need a track index too to calculate efficiencies
+ fTOFCalc[ih].good_tdc_pos = kFALSE;
+ fTOFCalc[ih].good_tdc_neg = kFALSE;
+ fTOFCalc[ih].pindex = ip;
+
+ Int_t paddle = hit->GetPaddleNumber() - 1;
+ fTOFCalc[ih].hit_paddle = paddle;
+ fTOFCalc[ih].good_raw_pad = paddle;
+
+ // Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) +
+ // fPlanes[ip]->GetPosOffset(); Double_t scinTrnsCoord =
+ // fTOFPInfo[ih].scinTrnsCoord; Double_t scinLongCoord =
+ // fTOFPInfo[ih].scinLongCoord;
+
+ Int_t fPIndex = GetScinIndex(ip, paddle);
+
+ if (fTOFPInfo[ih].onTrack) {
+ fGoodFlags[itrack][ip][iphit].onTrack = kTRUE;
+ if (fTOFPInfo[ih].keep_pos) { // 301
+ fTOFCalc[ih].good_tdc_pos = kTRUE;
+ fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
+ }
+ if (fTOFPInfo[ih].keep_neg) { //
+ fTOFCalc[ih].good_tdc_neg = kTRUE;
+ fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
+ }
+ // ** Calculate ave time for scin and error.
+ if (fTOFCalc[ih].good_tdc_pos) {
+ if (fTOFCalc[ih].good_tdc_neg) {
+ fTOFCalc[ih].scin_time =
+ (fTOFPInfo[ih].scin_pos_time + fTOFPInfo[ih].scin_neg_time) / 2.;
+ fTOFCalc[ih].scin_time_fp = (fTOFPInfo[ih].time_pos + fTOFPInfo[ih].time_neg) / 2.;
+
+ if (fTofUsingInvAdc) {
+ fTOFCalc[ih].scin_sigma =
+ TMath::Sqrt(fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
+ fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex]) /
+ 2.;
+ } else {
+ fTOFCalc[ih].scin_sigma =
+ TMath::Sqrt(fHodoSigmaPos[fPIndex] * fHodoSigmaPos[fPIndex] +
+ fHodoSigmaNeg[fPIndex] * fHodoSigmaNeg[fPIndex]) /
+ 2.;
+ }
+
+ fTOFCalc[ih].good_scin_time = kTRUE;
+ fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
+ } else {
+ fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_pos_time;
+ fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_pos;
+
+ if (fTofUsingInvAdc) {
+ fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex];
+ } else {
+ fTOFCalc[ih].scin_sigma = fHodoSigmaPos[fPIndex];
+ }
+
+ fTOFCalc[ih].good_scin_time = kTRUE;
+ fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
+ }
+ } else {
+ if (fTOFCalc[ih].good_tdc_neg) {
+ fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_neg_time;
+ fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_neg;
+ if (fTofUsingInvAdc) {
+ fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex];
+ } else {
+ fTOFCalc[ih].scin_sigma = fHodoSigmaNeg[fPIndex];
+ }
+ fTOFCalc[ih].good_scin_time = kTRUE;
+ fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
+ }
+ } // In h_tof.f this includes the following if condition for time at focal plane
+ // // because it is written in FORTRAN code
+
+ // c Get time at focal plane
+ if (fTOFCalc[ih].good_scin_time) {
+
+ // scin_time_fp doesn't need to be an array
+ // Is this any different than the average of time_pos and time_neg?
+ // Double_t scin_time_fp = ( fTOFPInfo[ih].time_pos +
+ // fTOFPInfo[ih].time_neg ) / 2.;
+ Double_t scin_time_fp = fTOFCalc[ih].scin_time_fp;
+
+ sumFPTime = sumFPTime + scin_time_fp;
+ nFPTime++;
+
+ fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
+ fNPlaneTime[ip]++;
+ fNScinHit[itrack]++;
+
+ if ((fTOFCalc[ih].good_tdc_pos) && (fTOFCalc[ih].good_tdc_neg)) {
+ nPmtHit[itrack] = nPmtHit[itrack] + 2;
+ } else {
+ nPmtHit[itrack] = nPmtHit[itrack] + 1;
+ }
+
+ fdEdX[itrack].push_back(0.0);
+ assert(fNScinHit[itrack] > 0 && (size_t)fNScinHit[itrack] < fdEdX[itrack].size() + 1);
+
+ // --------------------------------------------------------------------------------------------
+ if (fTOFCalc[ih].good_tdc_pos) {
+ if (fTOFCalc[ih].good_tdc_neg) {
+ fdEdX[itrack][fNScinHit[itrack] - 1] =
+ TMath::Sqrt(TMath::Max(0., hit->GetPosADC() * hit->GetNegADC()));
+ } else {
+ fdEdX[itrack][fNScinHit[itrack] - 1] = TMath::Max(0., hit->GetPosADC());
+ }
+ } else {
+ if (fTOFCalc[ih].good_tdc_neg) {
+ fdEdX[itrack][fNScinHit[itrack] - 1] = TMath::Max(0., hit->GetNegADC());
+ } else {
+ fdEdX[itrack][fNScinHit[itrack] - 1] = 0.0;
+ }
+ }
+ // --------------------------------------------------------------------------------------------
+
+ } // time at focal plane condition
+ } // on track condition
+
+ // ** See if there are any good time measurements in the plane.
+ if (fTOFCalc[ih].good_scin_time) {
+ fGoodPlaneTime[ip] = kTRUE;
+ fTOFCalc[ih].dedx = fdEdX[itrack][fNScinHit[itrack] - 1];
+ } else {
+ fTOFCalc[ih].dedx = 0.0;
+ }
+
} // Second loop over hits of a scintillator plane ends here
- theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
- if (fNumPlanesBetaCalc==4) theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
+ theTrack->SetGoodPlane3(fGoodPlaneTime[2] ? 1 : 0);
+ if (fNumPlanesBetaCalc == 4)
+ theTrack->SetGoodPlane4(fGoodPlaneTime[3] ? 1 : 0);
//
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
@@ -1380,142 +1419,136 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// * * Fit beta if there are enough time measurements (one upper, one lower)
// From h_tof_fit
- if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) &&
- ( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){
+ if (((fGoodPlaneTime[0]) || (fGoodPlaneTime[1])) &&
+ ((fGoodPlaneTime[2]) || (fGoodPlaneTime[3]))) {
- Double_t sumW = 0.;
- Double_t sumT = 0.;
- Double_t sumZ = 0.;
- Double_t sumZZ = 0.;
- Double_t sumTZ = 0.;
+ Double_t sumW = 0.;
+ Double_t sumT = 0.;
+ Double_t sumZ = 0.;
+ Double_t sumZZ = 0.;
+ Double_t sumTZ = 0.;
- for(Int_t ih=0; ih < nhits; ih++) {
- Int_t ip = fTOFPInfo[ih].planeIndex;
+ for (Int_t ih = 0; ih < nhits; ih++) {
+ Int_t ip = fTOFPInfo[ih].planeIndex;
- if ( fTOFCalc[ih].good_scin_time ) {
+ if (fTOFCalc[ih].good_scin_time) {
- Double_t scinWeight = 1 / ( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );
- Double_t zPosition = ( fPlanes[ip]->GetZpos()
- +( fTOFCalc[ih].hit_paddle % 2 ) *
- fPlanes[ip]->GetDzpos() );
+ Double_t scinWeight = 1 / (fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma);
+ Double_t zPosition =
+ (fPlanes[ip]->GetZpos() + (fTOFCalc[ih].hit_paddle % 2) * fPlanes[ip]->GetDzpos());
- sumW += scinWeight;
- sumT += scinWeight * fTOFCalc[ih].scin_time;
- sumZ += scinWeight * zPosition;
- sumZZ += scinWeight * ( zPosition * zPosition );
- sumTZ += scinWeight * zPosition * fTOFCalc[ih].scin_time;
+ sumW += scinWeight;
+ sumT += scinWeight * fTOFCalc[ih].scin_time;
+ sumZ += scinWeight * zPosition;
+ sumZZ += scinWeight * (zPosition * zPosition);
+ sumTZ += scinWeight * zPosition * fTOFCalc[ih].scin_time;
- } // condition of good scin time
- } // loop over hits
+ } // condition of good scin time
+ } // loop over hits
- Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
- Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
- Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
+ Double_t tmp = sumW * sumZZ - sumZ * sumZ;
+ Double_t t0 = (sumT * sumZZ - sumZ * sumTZ) / tmp;
+ Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
- if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
+ if (TMath::Abs(tmpDenom) > (1 / 10000000000.0)) {
- beta = tmp / tmpDenom;
- betaChiSq = 0.;
+ beta = tmp / tmpDenom;
+ betaChiSq = 0.;
- for(Int_t ih=0; ih < nhits; ih++) {
- Int_t ip = fTOFPInfo[ih].planeIndex;
+ for (Int_t ih = 0; ih < nhits; ih++) {
+ Int_t ip = fTOFPInfo[ih].planeIndex;
- if ( fTOFCalc[ih].good_scin_time ){
+ if (fTOFCalc[ih].good_scin_time) {
- Double_t zPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ih].hit_paddle % 2 ) *
- fPlanes[ip]->GetDzpos() );
- Double_t timeDif = ( fTOFCalc[ih].scin_time - t0 );
- betaChiSq += ( ( zPosition / beta - timeDif ) *
- ( zPosition / beta - timeDif ) ) /
- ( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );
+ Double_t zPosition = (fPlanes[ip]->GetZpos() +
+ (fTOFCalc[ih].hit_paddle % 2) * fPlanes[ip]->GetDzpos());
+ Double_t timeDif = (fTOFCalc[ih].scin_time - t0);
+ betaChiSq += ((zPosition / beta - timeDif) * (zPosition / beta - timeDif)) /
+ (fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma);
- } // condition for good scin time
- } // loop over hits
+ } // condition for good scin time
+ } // loop over hits
- Double_t pathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
- theTrack->GetPhi() * theTrack->GetPhi() );
- // Take angle into account
- beta = beta / pathNorm;
- beta = beta / 29.979; // velocity / c
+ Double_t pathNorm = TMath::Sqrt(1. + theTrack->GetTheta() * theTrack->GetTheta() +
+ theTrack->GetPhi() * theTrack->GetPhi());
+ // Take angle into account
+ beta = beta / pathNorm;
+ beta = beta / 29.979; // velocity / c
- } // condition for fTmpDenom
- else {
- beta = 0.;
- betaChiSq = -2.;
- } // else condition for fTmpDenom
+ } // condition for fTmpDenom
+ else {
+ beta = 0.;
+ betaChiSq = -2.;
+ } // else condition for fTmpDenom
} else {
- beta = 0.;
- betaChiSq = -1;
+ beta = 0.;
+ betaChiSq = -1;
}
- if ( nFPTime != 0 ){
- timeAtFP[itrack] = ( sumFPTime / nFPTime );
+ if (nFPTime != 0) {
+ timeAtFP[itrack] = (sumFPTime / nFPTime);
}
//
// ---------------------------------------------------------------------------
- Double_t FPTimeSum=0.0;
- Int_t nFPTimeSum=0;
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
- if ( fNPlaneTime[ip] != 0 ){
- fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
- FPTimeSum += fSumPlaneTime[ip];
- nFPTimeSum += fNPlaneTime[ip];
- } else {
- fFPTime[ip] = 1000. * ( ip + 1 );
- }
+ Double_t FPTimeSum = 0.0;
+ Int_t nFPTimeSum = 0;
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ if (fNPlaneTime[ip] != 0) {
+ fFPTime[ip] = (fSumPlaneTime[ip] / fNPlaneTime[ip]);
+ FPTimeSum += fSumPlaneTime[ip];
+ nFPTimeSum += fNPlaneTime[ip];
+ } else {
+ fFPTime[ip] = 1000. * (ip + 1);
+ }
}
- Double_t fptime=-1000;
- if (nFPTimeSum>0) fptime = FPTimeSum/nFPTimeSum;
- fFPTimeAll = fptime;
- Double_t dedx=0.0;
- for(UInt_t ih=0;ih<fTOFCalc.size();ih++) {
- if(fTOFCalc[ih].good_scin_time) {
- dedx = fTOFCalc[ih].dedx;
- break;
- }
+ Double_t fptime = -1000;
+ if (nFPTimeSum > 0)
+ fptime = FPTimeSum / nFPTimeSum;
+ fFPTimeAll = fptime;
+ Double_t dedx = 0.0;
+ for (UInt_t ih = 0; ih < fTOFCalc.size(); ih++) {
+ if (fTOFCalc[ih].good_scin_time) {
+ dedx = fTOFCalc[ih].dedx;
+ break;
+ }
}
theTrack->SetDedx(dedx);
theTrack->SetFPTime(fptime);
theTrack->SetBeta(beta);
- theTrack->SetBetaChi2( betaChiSq );
+ theTrack->SetBetaChi2(betaChiSq);
theTrack->SetNPMT(nPmtHit[itrack]);
- theTrack->SetFPTime( timeAtFP[itrack]);
-
+ theTrack->SetFPTime(timeAtFP[itrack]);
} // Main loop over tracks ends here.
} // If condition for at least one track
-
- //OriginalTrackEffTest();
+ // OriginalTrackEffTest();
TrackEffTest();
-
return 0;
-
}
-void THcHodoscope::TrackEffTest(void)
-{
+void THcHodoscope::TrackEffTest(void) {
// assume X planes are 0,2 and Y planes are 1,3
- std::array<int,4> PadLow = {fxLoScin[0], fyLoScin[0], fxLoScin[1], fyLoScin[1]};
- std::array<int,4> PadHigh = {fxHiScin[0], fyHiScin[0], fxHiScin[1], fyHiScin[1]};
-
- std::array<double,4> PadPosLo;
- std::array<double,4> PadPosHi;
+ std::array<int, 4> PadLow = {fxLoScin[0], fyLoScin[0], fxLoScin[1], fyLoScin[1]};
+ std::array<int, 4> PadHigh = {fxHiScin[0], fyHiScin[0], fxHiScin[1], fyHiScin[1]};
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
+ std::array<double, 4> PadPosLo;
+ std::array<double, 4> PadPosHi;
+
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
Double_t lowtemp = fPlanes[ip]->GetPosCenter(PadLow[ip] - 1) + fPlanes[ip]->GetPosOffset();
Double_t hitemp = fPlanes[ip]->GetPosCenter(PadHigh[ip] - 1) + fPlanes[ip]->GetPosOffset();
if (lowtemp < hitemp) {
- PadPosLo[ip]=lowtemp;
- PadPosHi[ip]=hitemp;
+ PadPosLo[ip] = lowtemp;
+ PadPosHi[ip] = hitemp;
} else {
- PadPosLo[ip]=hitemp;
- PadPosHi[ip]=lowtemp;
+ PadPosLo[ip] = hitemp;
+ PadPosHi[ip] = lowtemp;
}
- }
+ }
//{
// std::vector<int> test_vector = {1,2,5,6,7, 9,10, 20};
// auto test_res = hcana::find_discontinuity(test_vector);
@@ -1546,32 +1579,32 @@ void THcHodoscope::TrackEffTest(void)
using HitIterator = std::vector<THcHodoHit*>::iterator;
using HitRangeVector = typename std::vector<std::pair<HitIterator, HitIterator>>;
using ClusterPositions = typename std::vector<double>;
- std::map<int,std::vector<THcHodoHit*>> hit_vecs;
- std::vector<HitRangeVector> hit_clusters;
- std::vector<ClusterPositions> pos_clusters;
- std::vector<std::vector<int>> good_clusters;
- std::array<int,4> n_good_clusters;
+ std::map<int, std::vector<THcHodoHit*>> hit_vecs;
+ std::vector<HitRangeVector> hit_clusters;
+ std::vector<ClusterPositions> pos_clusters;
+ std::vector<std::vector<int>> good_clusters;
+ std::array<int, 4> n_good_clusters;
// Loop over each plane
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
// Get the hit array
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// Create vector for good hits.
hit_vecs[ip] = std::vector<THcHodoHit*>();
- for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
- THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
+ for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) {
+ THcHodoHit* hit = (THcHodoHit*)hodoHits->At(iphit);
// keep only those with "two good times"
- if ( hit->GetTwoGoodTimes() ) {
+ if (hit->GetTwoGoodTimes()) {
hit_vecs[ip].push_back(hit);
}
}
- // Cluster (group adjacent) hits
- if( hit_vecs[ip].size() > 0 ) {
+ // Cluster (group adjacent) hits
+ if (hit_vecs[ip].size() > 0) {
auto& hv = hit_vecs[ip];
hcana::sort(hv, [](const THcHodoHit* h1, const THcHodoHit* h2) {
- return h1->GetPaddleNumber() < h2->GetPaddleNumber();
+ return h1->GetPaddleNumber() < h2->GetPaddleNumber();
});
// get the clusters for this layer
@@ -1580,33 +1613,39 @@ void THcHodoscope::TrackEffTest(void)
return h1->GetPaddleNumber() + 1 == h2->GetPaddleNumber();
}));
} else {
- hit_clusters.push_back( HitRangeVector{} );
- }
+ hit_clusters.push_back(HitRangeVector{});
+ }
+ fPlanes[ip]->SetNumberClusters(hit_clusters.back().size());
// Calculate each cluster's mean position (in one coordinate)
ClusterPositions clust_positions;
std::vector<int> clust_bound;
- for(auto& r : hit_clusters.back()) {
- double a_pos = 0.0;
+ int ic = 0;
+ for (auto& r : hit_clusters.back()) {
+ double a_pos = 0.0;
int n_hit = 0;
- for(auto chit = r.first; chit < r.second; chit++) {
- a_pos += fPlanes[ip]->GetPosCenter((*chit)->GetPaddleNumber()-1) + fPlanes[ip]->GetPosOffset();
+ for (auto chit = r.first; chit < r.second; chit++) {
+ a_pos +=
+ fPlanes[ip]->GetPosCenter((*chit)->GetPaddleNumber() - 1) + fPlanes[ip]->GetPosOffset();
n_hit++;
}
// take the average position
double avg_pos = a_pos / double(n_hit);
+ fPlanes[ip]->SetCluster(ic, avg_pos);
+ fPlanes[ip]->SetClusterSize(ic, n_hit);
// Calculate if it is the bound by the upper and lower limits where we expect
// full tracks to be reconstructed.
- int is_bound = int((avg_pos >= PadPosLo[ip]) && (avg_pos <= PadPosHi[ip]));
+ int is_bound = int((avg_pos >= PadPosLo[ip]) && (avg_pos <= PadPosHi[ip]));
clust_bound.push_back(is_bound);
- if(is_bound) {
- // we are only interested in clusters positioned inthe "sweet spot"
+ if (is_bound) {
+ // we are only interested in clusters positioned inthe "sweet spot"
clust_positions.push_back(avg_pos);
n_good_clusters[ip]++;
if (n_hit > 10) {
_det_logger->warn("cluster in hodoscope track efficiency has {} hits", n_hit);
}
}
+ ic += 1;
}
pos_clusters.push_back(clust_positions);
good_clusters.push_back(clust_bound);
@@ -1617,54 +1656,73 @@ void THcHodoscope::TrackEffTest(void)
bool has_both_y_clusters = (n_good_clusters[1] > 0) && (n_good_clusters[3] > 0);
bool at_least_one_good_x_clusters = (n_good_clusters[0] > 0) || (n_good_clusters[2] > 0);
bool at_least_one_good_y_clusters = (n_good_clusters[1] > 0) || (n_good_clusters[3] > 0);
- bool good_x_match = false;
- bool good_y_match = false;
+ bool good_x_match = false;
+ bool good_y_match = false;
// Superficial matching. Find out if clusters in front and back could possibly belong to the same
// track. This done by checking to see if the position differences are less than 10 cm;
// ie, |x1-x2|<10 or |y1-y2| <10
// TODO do actual matching with tracks elsewhere
- if(has_both_x_clusters) {
- for(auto x1_pos : pos_clusters[0] ) {
- for(auto x2_pos : pos_clusters[2] ) {
- if( std::abs(x1_pos-x2_pos) < 10.0) {
+ if (has_both_x_clusters) {
+ int ic0 = 0;
+
+ // factor to project X1 to X2 z position
+ const double x1_projection_factor =
+ (1 + fRatio_xpfp_to_xfp * (fPlanes[2]->GetZpos() - fPlanes[0]->GetZpos()));
+
+ for (auto x1_pos : pos_clusters[0]) {
+ int ic2 = 0;
+ // project X1 to X2 Z position
+ Double_t x1_proj = x1_pos * x1_projection_factor;
+ for (auto x2_pos : pos_clusters[2]) {
+ if (std::abs(x1_proj - x2_pos) < trackeff_scint_xdiff_max) {
good_x_match = true;
+ fPlanes[0]->SetClusterUsedFlag(ic0, 1.);
+ fPlanes[2]->SetClusterUsedFlag(ic2, 1.);
break;
}
+ ic2 += 1;
}
+ ic0 += 1;
}
}
- if(has_both_y_clusters) {
- for(auto y1_pos : pos_clusters[1] ) {
- for(auto y2_pos : pos_clusters[3] ) {
- if( std::abs(y1_pos-y2_pos) < 10.0) {
+ if (has_both_y_clusters) {
+ int ic1 = 0;
+ for (auto y1_pos : pos_clusters[1]) {
+ int ic3 = 0;
+ for (auto y2_pos : pos_clusters[3]) {
+ if (std::abs(y1_pos - y2_pos) < trackeff_scint_ydiff_max) {
good_y_match = true;
+ fPlanes[1]->SetClusterUsedFlag(ic1, 1.);
+ fPlanes[3]->SetClusterUsedFlag(ic3, 1.);
break;
}
+ ic3 += 1;
}
+ ic1 += 1;
}
}
fGoodScinHits = 0;
- // (2+2) 4 good hits
- if( good_y_match && good_x_match ) {
- fGoodScinHits = 1;
- }
- // (2+1) 3 good hits
- if( at_least_one_good_y_clusters && good_x_match ) {
- fGoodScinHits = 1;
- }
- // (1+2) 3 good hits
- if( at_least_one_good_x_clusters && good_y_match ) {
- fGoodScinHits = 1;
+ // (2+2) 4 good hits
+ if (fTrackEffTestNScinPlanes == 4 || fTrackEffTestNScinPlanes == 3) {
+ if (good_y_match && good_x_match) {
+ fGoodScinHits = 1;
+ }
+ // (2+1) 3 good hits
+ if (at_least_one_good_y_clusters && good_x_match) {
+ fGoodScinHits = 1;
+ }
+ // (1+2) 3 good hits
+ if (at_least_one_good_x_clusters && good_y_match) {
+ fGoodScinHits = 1;
+ }
}
-
}
//
//
//
-void THcHodoscope::OriginalTrackEffTest(void)
-{
+void THcHodoscope::OriginalTrackEffTest(void) {
/**
Translation of h_track_tests.f file for tracking efficiency
*/
@@ -1672,21 +1730,21 @@ void THcHodoscope::OriginalTrackEffTest(void)
//************************now look at some hodoscope tests
// *second, we move the scintillators. here we use scintillator cuts to see
// *if a track should have been found.
- cout << " enter track eff" << fNumPlanesBetaCalc << endl;
- for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
- cout << " loop over planes " << ip+1 << endl;
+ cout << " enter track eff" << fNumPlanesBetaCalc << endl;
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ cout << " loop over planes " << ip + 1 << endl;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// TClonesArray* scinPosTDC = fPlanes[ip]->GetPosTDC();
// TClonesArray* scinNegTDC = fPlanes[ip]->GetNegTDC();
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
- cout << " hits = " << fNScinHits[ip] << endl;
- for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
- Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
+ cout << " hits = " << fNScinHits[ip] << endl;
+ for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++) {
+ Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber() - 1;
fScinHitPaddle[ip][paddle] = 1;
- cout << " hit = " << iphit+1 << " " << paddle+1 << endl;
+ cout << " hit = " << iphit + 1 << " " << paddle + 1 << endl;
}
}
@@ -1695,7 +1753,7 @@ void THcHodoscope::OriginalTrackEffTest(void)
// *Wwe count the number of three adjacent scintillators too. (A signle track
// *shouldn't fire three adjacent scintillators.
- for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
// Planes ip = 0 = 1X
// Planes ip = 2 = 2X
fNClust.push_back(0);
@@ -1706,34 +1764,33 @@ void THcHodoscope::OriginalTrackEffTest(void)
// *number of scintillators.
cout << " looking for cluster in x planes" << endl;
Int_t icount;
- for (Int_t ip = 0; ip < 3; ip +=2 ) {
+ for (Int_t ip = 0; ip < 3; ip += 2) {
icount = 0;
- if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
- cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
+ if (fScinHitPaddle[ip][0] == 1)
+ icount++;
+ cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
- for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
+ for (Int_t ipaddle = 0; ipaddle < (Int_t)fNPaddle[ip] - 1; ipaddle++) {
// !look for number of clusters of 1 or more hits
- if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
- ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
- icount ++;
- cout << " paddle = " << ipaddle+1 << " " << icount << endl;
-
+ if ((fScinHitPaddle[ip][ipaddle] == 0) && (fScinHitPaddle[ip][ipaddle + 1] == 1))
+ icount++;
+ cout << " paddle = " << ipaddle + 1 << " " << icount << endl;
+
} // Loop over paddles
- cout << "Two cluster in plane = " << ip+1 << " " << icount << endl;
+ cout << "Two cluster in plane = " << ip + 1 << " " << icount << endl;
fNClust[ip] = icount;
- icount = 0;
+ icount = 0;
- for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
+ for (Int_t ipaddle = 0; ipaddle < (Int_t)fNPaddle[ip] - 2; ipaddle++) {
// !look for three or more adjacent hits
- if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
- ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
- ( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
- icount ++;
+ if ((fScinHitPaddle[ip][ipaddle] == 1) && (fScinHitPaddle[ip][ipaddle + 1] == 1) &&
+ (fScinHitPaddle[ip][ipaddle + 2] == 1))
+ icount++;
} // Second loop over paddles
- cout << "Three clusters in plane = " << ip+1 << " " << icount << endl;
+ cout << "Three clusters in plane = " << ip + 1 << " " << icount << endl;
- if ( icount > 0 )
+ if (icount > 0)
fThreeScin[ip] = 1;
} // Loop over X plane
@@ -1742,224 +1799,248 @@ void THcHodoscope::OriginalTrackEffTest(void)
// *number of scintillators.
cout << " looking for cluster in y planes" << endl;
- for (Int_t ip = 1; ip < fNumPlanesBetaCalc; ip +=2 ) {
+ for (Int_t ip = 1; ip < fNumPlanesBetaCalc; ip += 2) {
// Planes ip = 1 = 1Y
// Planes ip = 3 = 2Y
icount = 0;
- if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
- cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
+ if (fScinHitPaddle[ip][0] == 1)
+ icount++;
+ cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
- for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
+ for (Int_t ipaddle = 0; ipaddle < (Int_t)fNPaddle[ip] - 1; ipaddle++) {
// !look for number of clusters of 1 or more hits
- if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
- ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
- icount ++;
- cout << " paddle = " << ipaddle+1 << " " << icount << endl;
+ if ((fScinHitPaddle[ip][ipaddle] == 0) && (fScinHitPaddle[ip][ipaddle + 1] == 1))
+ icount++;
+ cout << " paddle = " << ipaddle + 1 << " " << icount << endl;
} // Loop over Y paddles
- cout << "Two cluster in plane = " << ip+1 << " " << icount << endl;
+ cout << "Two cluster in plane = " << ip + 1 << " " << icount << endl;
fNClust[ip] = icount;
- icount = 0;
+ icount = 0;
- for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
+ for (Int_t ipaddle = 0; ipaddle < (Int_t)fNPaddle[ip] - 2; ipaddle++) {
// !look for three or more adjacent hits
- if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
- ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
- ( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
- icount ++;
+ if ((fScinHitPaddle[ip][ipaddle] == 1) && (fScinHitPaddle[ip][ipaddle + 1] == 1) &&
+ (fScinHitPaddle[ip][ipaddle + 2] == 1))
+ icount++;
} // Second loop over Y paddles
- cout << "Three clusters in plane = " << ip+1 << " " << icount << endl;
+ cout << "Three clusters in plane = " << ip + 1 << " " << icount << endl;
- if ( icount > 0 )
+ if (icount > 0)
fThreeScin[ip] = 1;
- }// Loop over Y planes
-
+ } // Loop over Y planes
// *next we mask out the edge scintillators, and look at triggers that happened
// *at the center of the acceptance. To change which scins are in the mask
// *change the values of h*loscin and h*hiscin in htracking.param
// fGoodScinHits = 0;
- for (Int_t ifidx = fxLoScin[0]; ifidx < (Int_t) fxHiScin[0]; ifidx ++ ){
+ for (Int_t ifidx = fxLoScin[0]; ifidx < (Int_t)fxHiScin[0]; ifidx++) {
fGoodScinHitsX.push_back(0);
}
- fHitSweet1X=0;
- fHitSweet2X=0;
- fHitSweet1Y=0;
- fHitSweet2Y=0;
+ fHitSweet1X = 0;
+ fHitSweet2X = 0;
+ fHitSweet1Y = 0;
+ fHitSweet2Y = 0;
// *first x plane. first see if there are hits inside the scin region
- for (Int_t ifidx = fxLoScin[0]-1; ifidx < fxHiScin[0]; ifidx ++ ){
- if ( fScinHitPaddle[0][ifidx] == 1 ){
- fHitSweet1X = 1;
+ for (Int_t ifidx = fxLoScin[0] - 1; ifidx < fxHiScin[0]; ifidx++) {
+ if (fScinHitPaddle[0][ifidx] == 1) {
+ fHitSweet1X = 1;
fSweet1XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
- for (Int_t ifidx = 0; ifidx < fxLoScin[0]-1; ifidx ++ ){
- if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
+ for (Int_t ifidx = 0; ifidx < fxLoScin[0] - 1; ifidx++) {
+ if (fScinHitPaddle[0][ifidx] == 1) {
+ fHitSweet1X = -1;
+ }
}
- for (Int_t ifidx = fxHiScin[0]; ifidx < (Int_t) fNPaddle[0]; ifidx ++ ){
- if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
+ for (Int_t ifidx = fxHiScin[0]; ifidx < (Int_t)fNPaddle[0]; ifidx++) {
+ if (fScinHitPaddle[0][ifidx] == 1) {
+ fHitSweet1X = -1;
+ }
}
// *second x plane. first see if there are hits inside the scin region
- for (Int_t ifidx = fxLoScin[1]-1; ifidx < fxHiScin[1]; ifidx ++ ){
- if ( fScinHitPaddle[2][ifidx] == 1 ){
- fHitSweet2X = 1;
+ for (Int_t ifidx = fxLoScin[1] - 1; ifidx < fxHiScin[1]; ifidx++) {
+ if (fScinHitPaddle[2][ifidx] == 1) {
+ fHitSweet2X = 1;
fSweet2XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
- for (Int_t ifidx = 0; ifidx < fxLoScin[1]-1; ifidx ++ ){
- if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
+ for (Int_t ifidx = 0; ifidx < fxLoScin[1] - 1; ifidx++) {
+ if (fScinHitPaddle[2][ifidx] == 1) {
+ fHitSweet2X = -1;
+ }
}
- for (Int_t ifidx = fxHiScin[1]; ifidx < (Int_t) fNPaddle[2]; ifidx ++ ){
- if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
+ for (Int_t ifidx = fxHiScin[1]; ifidx < (Int_t)fNPaddle[2]; ifidx++) {
+ if (fScinHitPaddle[2][ifidx] == 1) {
+ fHitSweet2X = -1;
+ }
}
// *first y plane. first see if there are hits inside the scin region
- for (Int_t ifidx = fyLoScin[0]-1; ifidx < fyHiScin[0]; ifidx ++ ){
- if ( fScinHitPaddle[1][ifidx] == 1 ){
- fHitSweet1Y = 1;
+ for (Int_t ifidx = fyLoScin[0] - 1; ifidx < fyHiScin[0]; ifidx++) {
+ if (fScinHitPaddle[1][ifidx] == 1) {
+ fHitSweet1Y = 1;
fSweet1YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
- for (Int_t ifidx = 0; ifidx < fyLoScin[0]-1; ifidx ++ ){
- if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
+ for (Int_t ifidx = 0; ifidx < fyLoScin[0] - 1; ifidx++) {
+ if (fScinHitPaddle[1][ifidx] == 1) {
+ fHitSweet1Y = -1;
+ }
}
- for (Int_t ifidx = fyHiScin[0]; ifidx < (Int_t) fNPaddle[1]; ifidx ++ ){
- if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
+ for (Int_t ifidx = fyHiScin[0]; ifidx < (Int_t)fNPaddle[1]; ifidx++) {
+ if (fScinHitPaddle[1][ifidx] == 1) {
+ fHitSweet1Y = -1;
+ }
}
// *second y plane. first see if there are hits inside the scin region
- for (Int_t ifidx = fyLoScin[1]-1; ifidx < fyHiScin[1]; ifidx ++ ){
- if ( fScinHitPaddle[3][ifidx] == 1 ){
- fHitSweet2Y = 1;
+ for (Int_t ifidx = fyLoScin[1] - 1; ifidx < fyHiScin[1]; ifidx++) {
+ if (fScinHitPaddle[3][ifidx] == 1) {
+ fHitSweet2Y = 1;
fSweet2YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
- for (Int_t ifidx = 0; ifidx < fyLoScin[1]-1; ifidx ++ ){
- if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
+ for (Int_t ifidx = 0; ifidx < fyLoScin[1] - 1; ifidx++) {
+ if (fScinHitPaddle[3][ifidx] == 1) {
+ fHitSweet2Y = -1;
+ }
}
- for (Int_t ifidx = fyHiScin[1]; ifidx < (Int_t) fNPaddle[3]; ifidx ++ ){
- if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
+ for (Int_t ifidx = fyHiScin[1]; ifidx < (Int_t)fNPaddle[3]; ifidx++) {
+ if (fScinHitPaddle[3][ifidx] == 1) {
+ fHitSweet2Y = -1;
+ }
}
fTestSum = fHitSweet1X + fHitSweet2X + fHitSweet1Y + fHitSweet2Y;
// * now define a 3/4 or 4/4 trigger of only good scintillators the value
// * is specified in htracking
- if ( fTestSum >= fTrackEffTestNScinPlanes ){
+ if (fTestSum >= fTrackEffTestNScinPlanes) {
fGoodScinHits = 1;
- for (Int_t ifidx = fxLoScin[0]; ifidx < fxHiScin[0]; ifidx ++ ){
- if ( fSweet1XScin == ifidx )
- fGoodScinHitsX[ifidx] = 1;
+ for (Int_t ifidx = fxLoScin[0]; ifidx < fxHiScin[0]; ifidx++) {
+ if (fSweet1XScin == ifidx)
+ fGoodScinHitsX[ifidx] = 1;
}
}
// * require front/back hodoscopes be close to each other
- if ( ( fGoodScinHits == 1 ) && ( fTrackEffTestNScinPlanes == 4 ) ){
- if ( TMath::Abs( fSweet1XScin - fSweet2XScin ) > 3 )
+ if ((fGoodScinHits == 1) && (fTrackEffTestNScinPlanes == 4)) {
+ if (TMath::Abs(fSweet1XScin - fSweet2XScin) > 3)
fGoodScinHits = 0;
- if ( TMath::Abs( fSweet1YScin - fSweet2YScin ) > 2 )
+ if (TMath::Abs(fSweet1YScin - fSweet2YScin) > 2)
fGoodScinHits = 0;
}
-//
+ //
}
//_____________________________________________________________________________
-Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
-{
- Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+Int_t THcHodoscope::FineProcess(TClonesArray& tracks) {
+ Int_t Ntracks = tracks.GetLast() + 1; // Number of reconstructed tracks
Double_t hitPos;
Double_t hitDistance;
- Int_t ih=0;
- for (Int_t itrk=0; itrk<Ntracks; itrk++) {
- THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
- if (theTrack->GetIndex()==0) {
- fBeta=theTrack->GetBeta();
- fFPTimeAll=theTrack->GetFPTime();
+ Int_t ih = 0;
+ for (Int_t itrk = 0; itrk < Ntracks; itrk++) {
+ THaTrack* theTrack = static_cast<THaTrack*>(tracks[itrk]);
+ if (theTrack->GetIndex() == 0) {
+ fBeta = theTrack->GetBeta();
+ fFPTimeAll = theTrack->GetFPTime();
_basic_data.fBeta = fBeta;
_basic_data.fFPTimeAll = fFPTimeAll;
- Double_t shower_track_enorm = theTrack->GetEnergy()/theTrack->GetP();
- for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
- Double_t pl_xypos=0;
- Double_t pl_zpos=0;
- Int_t num_good_pad=0;
- TClonesArray* hodoHits = fPlanes[ip]->GetHits();
- for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
- THcHodoHit *hit = fTOFPInfo[ih].hit;
- if ( fTOFCalc[ih].good_tdc_pos && fTOFCalc[ih].good_tdc_neg ) {
- Bool_t sh_pid=(shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi);
- Bool_t beta_pid=( fBeta > fTOFCalib_beta_lo && fBeta < fTOFCalib_beta_hi);
- // cer_pid is true if there is no Cherenkov detector
- Bool_t cer_pid=(fCherenkov?(fCherenkov->GetCerNPE()>fTOFCalib_cer_lo):kTRUE);
- if(fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && sh_pid && beta_pid && cer_pid) {
- fDumpOut << fixed << setprecision(2);
- fDumpOut << showpoint << " 1" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber() << setw(10) << hit->GetPosTDC()*fScinTdcToTime << setw(10) << fTOFPInfo[ih].pathp << setw(10) << fTOFPInfo[ih].zcor << setw(10) << fTOFPInfo[ih].time_pos << setw(10) << hit->GetPosADC() << endl;
- fDumpOut << showpoint << " 2" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber() << setw(10) << hit->GetNegTDC()*fScinTdcToTime << setw(10) << fTOFPInfo[ih].pathn << setw(10) << fTOFPInfo[ih].zcor << setw(10) << fTOFPInfo[ih].time_neg << setw(10) << hit->GetNegADC() << endl;
- }
- Int_t padind = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
- pl_xypos+=fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
- pl_zpos+=fPlanes[ip]->GetZpos()+ (padind%2)*fPlanes[ip]->GetDzpos();
- num_good_pad++;
- }
- ih++;
- // cout << ip << " " << iphit << " " << fGoodFlags[itrk][ip][iphit].goodScinTime << " " << fGoodFlags[itrk][ip][iphit].goodTdcPos << " " << fGoodFlags[itrk][ip][iphit].goodTdcNeg << endl;
- }
- hitDistance=kBig;
- if (num_good_pad !=0 ) {
- pl_xypos=pl_xypos/num_good_pad;
- pl_zpos=pl_zpos/num_good_pad;
- hitPos = theTrack->GetY() + theTrack->GetPhi()*pl_zpos ;
- if (ip%2 == 0) hitPos = theTrack->GetX() + theTrack->GetTheta()*pl_zpos ;
- hitDistance = hitPos - pl_xypos;
- fPlanes[ip]->SetTrackXPosition(theTrack->GetX() + theTrack->GetTheta()*pl_zpos );
- fPlanes[ip]->SetTrackYPosition(theTrack->GetY() + theTrack->GetPhi()*pl_zpos );
- }
- // cout << " ip " << ip << " " << hitPos << " " << pl_xypos << " " << pl_zpos << " " << hitDistance << endl;
- fPlanes[ip]->SetHitDistance(hitDistance);
+ Double_t shower_track_enorm = theTrack->GetEnergy() / theTrack->GetP();
+ for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
+ Double_t pl_xypos = 0;
+ Double_t pl_zpos = 0;
+ Int_t num_good_pad = 0;
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+ for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) {
+ THcHodoHit* hit = fTOFPInfo[ih].hit;
+ if (fTOFCalc[ih].good_tdc_pos && fTOFCalc[ih].good_tdc_neg) {
+ Bool_t sh_pid = (shower_track_enorm > fTOFCalib_shtrk_lo &&
+ shower_track_enorm < fTOFCalib_shtrk_hi);
+ Bool_t beta_pid = (fBeta > fTOFCalib_beta_lo && fBeta < fTOFCalib_beta_hi);
+ // cer_pid is true if there is no Cherenkov detector
+ Bool_t cer_pid = (fCherenkov ? (fCherenkov->GetCerNPE() > fTOFCalib_cer_lo) : kTRUE);
+ if (fDumpTOF && Ntracks == 1 && fGoodEventTOFCalib && sh_pid && beta_pid && cer_pid) {
+ fDumpOut << fixed << setprecision(2);
+ fDumpOut << showpoint << " 1" << setw(3) << ip + 1 << setw(3)
+ << hit->GetPaddleNumber() << setw(10) << hit->GetPosTDC() * fScinTdcToTime
+ << setw(10) << fTOFPInfo[ih].pathp << setw(10) << fTOFPInfo[ih].zcor
+ << setw(10) << fTOFPInfo[ih].time_pos << setw(10) << hit->GetPosADC()
+ << endl;
+ fDumpOut << showpoint << " 2" << setw(3) << ip + 1 << setw(3)
+ << hit->GetPaddleNumber() << setw(10) << hit->GetNegTDC() * fScinTdcToTime
+ << setw(10) << fTOFPInfo[ih].pathn << setw(10) << fTOFPInfo[ih].zcor
+ << setw(10) << fTOFPInfo[ih].time_neg << setw(10) << hit->GetNegADC()
+ << endl;
+ }
+ Int_t padind = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber() - 1;
+ pl_xypos += fPlanes[ip]->GetPosCenter(padind) + fPlanes[ip]->GetPosOffset();
+ pl_zpos += fPlanes[ip]->GetZpos() + (padind % 2) * fPlanes[ip]->GetDzpos();
+ num_good_pad++;
+ }
+ ih++;
+ // cout << ip << " " << iphit << " " << fGoodFlags[itrk][ip][iphit].goodScinTime <<
+ //" " << fGoodFlags[itrk][ip][iphit].goodTdcPos << " " <<
+ // fGoodFlags[itrk][ip][iphit].goodTdcNeg << endl;
+ }
+ hitDistance = kBig;
+ if (num_good_pad != 0) {
+ pl_xypos = pl_xypos / num_good_pad;
+ pl_zpos = pl_zpos / num_good_pad;
+ hitPos = theTrack->GetY() + theTrack->GetPhi() * pl_zpos;
+ if (ip % 2 == 0)
+ hitPos = theTrack->GetX() + theTrack->GetTheta() * pl_zpos;
+ hitDistance = hitPos - pl_xypos;
+ fPlanes[ip]->SetTrackXPosition(theTrack->GetX() + theTrack->GetTheta() * pl_zpos);
+ fPlanes[ip]->SetTrackYPosition(theTrack->GetY() + theTrack->GetPhi() * pl_zpos);
+ }
+ // cout << " ip " << ip << " " << hitPos << " " << pl_xypos << " " << pl_zpos << " "
+ // << hitDistance << endl;
+ fPlanes[ip]->SetHitDistance(hitDistance);
}
- if(ih>0&&fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi ) {
- fDumpOut << "0 " << endl;
+ if (ih > 0 && fDumpTOF && Ntracks == 1 && fGoodEventTOFCalib &&
+ shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi) {
+ fDumpOut << "0 " << endl;
}
}
- } //over tracks
+ } // over tracks
//
return 0;
}
//_____________________________________________________________________________
-Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
+Int_t THcHodoscope::GetScinIndex(Int_t nPlane, Int_t nPaddle) {
// GN: Return the index of a scintillator given the plane # and the paddle #
// This assumes that both planes and
// paddles start counting from 0!
// Result also counts from 0.
- return fNPlanes*nPaddle+nPlane;
+ return fNPlanes * nPaddle + nPlane;
}
//_____________________________________________________________________________
-Int_t THcHodoscope::GetScinIndex( Int_t nSide, Int_t nPlane, Int_t nPaddle ) {
- return nSide*fMaxHodoScin+fNPlanes*nPaddle+nPlane-1;
+Int_t THcHodoscope::GetScinIndex(Int_t nSide, Int_t nPlane, Int_t nPaddle) {
+ return nSide * fMaxHodoScin + fNPlanes * nPaddle + nPlane - 1;
}
//_____________________________________________________________________________
-Double_t THcHodoscope::GetPathLengthCentral() {
- return fPathLengthCentral;
-}
+Double_t THcHodoscope::GetPathLengthCentral() { return fPathLengthCentral; }
//_____________________________________________________________________________
-Int_t THcHodoscope::End(THaRunBase* run)
-{
+Int_t THcHodoscope::End(THaRunBase* run) {
MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
return 0;
}
ClassImp(THcHodoscope)
-////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index 8a0a1c55ef6d021f5a533f04181db84d74c2429c..1d30cb95349609d2a6f45d11770a00e08d9d79a7 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -185,7 +185,9 @@ protected:
TH1F *hTime;
// Calibration
-
+ Double_t fRatio_xpfp_to_xfp;
+ Double_t trackeff_scint_ydiff_max ;
+ Double_t trackeff_scint_xdiff_max ;
// Per-event data
Bool_t fSHMS;
Bool_t fGoodStartTime;
diff --git a/src/THcRawAdcHit.cxx b/src/THcRawAdcHit.cxx
index 6e2bfc7884fdf893e7d205cf54612caf2d042b0d..a5a86f4883514ae0009b3d6d94470b5afe79691e 100644
--- a/src/THcRawAdcHit.cxx
+++ b/src/THcRawAdcHit.cxx
@@ -245,8 +245,8 @@ void THcRawAdcHit::Clear(Option_t* opt) {
void THcRawAdcHit::SetData(Int_t data) {
if (fNPulses >= fMaxNPulses) {
- _hit_logger->error("THcRawAdcHit::SetData: too many pulses! Ignoring pulse {}",fNPulses);
- //throw std::out_of_range("`THcRawAdcHit::SetData`: too many pulses!");
+ _hit_logger->error("THcRawAdcHit::SetData: too many pulses! Ignoring pulse {}", fNPulses);
+ // throw std::out_of_range("`THcRawAdcHit::SetData`: too many pulses!");
return;
}
fPulseInt[fNPulses] = data;
@@ -260,8 +260,8 @@ void THcRawAdcHit::SetRefTime(Int_t refTime) {
void THcRawAdcHit::SetSample(Int_t data) {
if (fNSamples >= fMaxNSamples) {
- //throw std::out_of_range("`THcRawAdcHit::SetSample`: too many samples!");
- _hit_logger->error("THcRawAdcHit::SetSample: too many samples! Ignoring sample {}",fNSamples);
+ // throw std::out_of_range("`THcRawAdcHit::SetSample`: too many samples!");
+ _hit_logger->error("THcRawAdcHit::SetSample: too many samples! Ignoring sample {}", fNSamples);
}
fSample[fNSamples] = data;
++fNSamples;
@@ -269,8 +269,9 @@ void THcRawAdcHit::SetSample(Int_t data) {
void THcRawAdcHit::SetDataTimePedestalPeak(Int_t data, Int_t time, Int_t pedestal, Int_t peak) {
if (fNPulses >= fMaxNPulses) {
- _hit_logger->error("THcRawAdcHit::SetDataTimePedestalPeak: too many pulses! Ignoring pulse {}",fNPulses);
- //throw std::out_of_range("`THcRawAdcHit::SetData`: too many pulses!");
+ _hit_logger->error("THcRawAdcHit::SetDataTimePedestalPeak: too many pulses! Ignoring pulse {}",
+ fNPulses);
+ // throw std::out_of_range("`THcRawAdcHit::SetData`: too many pulses!");
return;
}
fPulseInt[fNPulses] = data;
@@ -283,11 +284,13 @@ void THcRawAdcHit::SetDataTimePedestalPeak(Int_t data, Int_t time, Int_t pedesta
Int_t THcRawAdcHit::GetRawData(UInt_t iPulse) const {
if (iPulse >= fNPulses && iPulse != 0) {
- //TString msg = TString::Format(
+ // TString msg = TString::Format(
// "`THcRawAdcHit::GetRawData`: requested pulse %d where only %d pulses available!", iPulse,
// fNPulses);
- //throw std::out_of_range(msg.Data());
- _hit_logger->error("THcRawAdcHit::GetRawData: requested pulse {} where only {} pulses available!", iPulse, fNPulses);
+ // throw std::out_of_range(msg.Data());
+ _hit_logger->error(
+ "THcRawAdcHit::GetRawData: requested pulse {} where only {} pulses available!", iPulse,
+ fNPulses);
return 0;
} else if (iPulse >= fNPulses && iPulse == 0) {
return 0;
@@ -311,8 +314,8 @@ Double_t THcRawAdcHit::GetAverage(UInt_t iSampleLow, UInt_t iSampleHigh) const {
Int_t THcRawAdcHit::GetIntegral(UInt_t iSampleLow, UInt_t iSampleHigh) const {
if (iSampleHigh >= fNSamples || iSampleLow >= fNSamples) {
- //TString msg = TString::Format("`THcRawAdcHit::GetAverage`: not this many samples available!");
- //throw std::out_of_range(msg.Data());
+ // TString msg = TString::Format("`THcRawAdcHit::GetAverage`: not this many samples
+ // available!"); throw std::out_of_range(msg.Data());
_hit_logger->error("THcRawAdcHit::GetRawData: not this many samples available!");
return 0;
} else {
@@ -343,11 +346,13 @@ Int_t THcRawAdcHit::GetPulseIntRaw(UInt_t iPulse) const {
} else if (iPulse == 0) {
return 0;
} else {
- //TString msg = TString::Format(
+ // TString msg = TString::Format(
// "`THcRawAdcHit::GetPulseIntRaw`: Trying to get pulse %d where only %d pulses available!",
// iPulse, fNPulses);
- //throw std::out_of_range(msg.Data());
- _hit_logger->error("THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!", iPulse, fNPulses);
+ // throw std::out_of_range(msg.Data());
+ _hit_logger->error(
+ "THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!",
+ iPulse, fNPulses);
return 0;
}
}
@@ -358,11 +363,13 @@ Int_t THcRawAdcHit::GetPulseAmpRaw(UInt_t iPulse) const {
} else if (iPulse == 0) {
return 0;
} else {
- //TString msg = TString::Format(
+ // TString msg = TString::Format(
// "`THcRawAdcHit::GetPulseAmpRaw`: Trying to get pulse %d where only %d pulses available!",
// iPulse, fNPulses);
- //throw std::out_of_range(msg.Data());
- _hit_logger->error("THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!", iPulse, fNPulses);
+ // throw std::out_of_range(msg.Data());
+ _hit_logger->error(
+ "THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!",
+ iPulse, fNPulses);
return 0;
}
}
@@ -373,11 +380,13 @@ Int_t THcRawAdcHit::GetPulseTimeRaw(UInt_t iPulse) const {
} else if (iPulse == 0) {
return 0;
} else {
- //TString msg = TString::Format(
+ // TString msg = TString::Format(
// "`THcRawAdcHit::GetPulseTimeRaw`: Trying to get pulse %d where only %d pulses available!",
// iPulse, fNPulses);
- //throw std::out_of_range(msg.Data());
- _hit_logger->error("THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!", iPulse, fNPulses);
+ // throw std::out_of_range(msg.Data());
+ _hit_logger->error(
+ "THcRawAdcHit::GetPulseIntRaw: Trying to get pulse {} where only {} pulses available!",
+ iPulse, fNPulses);
return 0;
}
}
@@ -386,11 +395,12 @@ Int_t THcRawAdcHit::GetSampleRaw(UInt_t iSample) const {
if (iSample < fNSamples) {
return fSample[iSample];
} else {
- //TString msg = TString::Format(
+ // TString msg = TString::Format(
// "`THcRawAdcHit::GetSampleRaw`: Trying to get sample %d where only %d samples available!",
// iSample, fNSamples);
- //throw std::out_of_range(msg.Data());
- _hit_logger->error("THcRawAdcHit::GetSampleRaw: Trying to get sample {} where only {} samples available!",
+ // throw std::out_of_range(msg.Data());
+ _hit_logger->error(
+ "THcRawAdcHit::GetSampleRaw: Trying to get sample {} where only {} samples available!",
iSample, fNSamples);
return 0;
}
diff --git a/src/THcScintPlaneCluster.cxx b/src/THcScintPlaneCluster.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f5d56819800c8915dc341b12dd4456ddebb23e50
--- /dev/null
+++ b/src/THcScintPlaneCluster.cxx
@@ -0,0 +1,8 @@
+/** \class THcScintPlaneCluster
+ \ingroup DetSupport
+
+*/
+
+#include "THcScintPlaneCluster.h"
+
+ClassImp(THcScintPlaneCluster)
diff --git a/src/THcScintPlaneCluster.h b/src/THcScintPlaneCluster.h
new file mode 100644
index 0000000000000000000000000000000000000000..d6c31aabe4075db4fc2c3ba25a5af98c4dcfd651
--- /dev/null
+++ b/src/THcScintPlaneCluster.h
@@ -0,0 +1,43 @@
+#ifndef ROOT_THcScintPlaneCluster
+#define ROOT_THcScintPlaneCluster
+
+/////////////////////////////////////////////////////////////////////////////
+// //
+// THcScintPlaneCluster //
+// //
+/////////////////////////////////////////////////////////////////////////////
+
+#include "TObject.h"
+#include <cstdio>
+
+class THcScintPlaneCluster : public TObject {
+
+ public:
+ THcScintPlaneCluster(Int_t cluster=0, Double_t pos=0.0) :
+ fClusterNumber(cluster), fClusterPosition(pos) {}
+ virtual ~THcScintPlaneCluster() {}
+
+ Int_t GetClusterNumber() {return fClusterNumber;}
+ Double_t GetClusterPosition() {return fClusterPosition;}
+ Double_t GetClusterSize() {return fClusterSize;}
+ Double_t GetClusterFlag() {return fClusterFlag;}
+ Double_t GetClusterUsedFlag() {return fClusterUsedFlag;}
+
+ virtual void Set(Int_t cluster, Double_t pos)
+ { fClusterNumber=cluster; fClusterPosition=pos; fClusterFlag=0.; fClusterUsedFlag=0.;}
+
+ void SetSize(Double_t size) {fClusterSize=size;}
+ void SetFlag(Double_t flag) {fClusterFlag=flag;}
+ void SetUsedFlag(Double_t flag) {fClusterUsedFlag=flag;}
+
+ private:
+ Int_t fClusterNumber;
+ Double_t fClusterPosition;
+ Double_t fClusterSize;
+ Double_t fClusterFlag;
+ Double_t fClusterUsedFlag;
+
+ ClassDef(THcScintPlaneCluster,0); // Single signal value and wire/counter number
+};
+/////////////////////////////////////////////////////////////////
+#endif
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index c9a6b2e65f822258fbae3122bd2006b5f3ef6708..1e4578907372a505e9b7f832e44d1e116353ff17 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -8,6 +8,7 @@ The THcHodoscope class instatiates one of these objects per plane.
*/
#include "TMath.h"
#include "THcScintillatorPlane.h"
+#include "THcScintPlaneCluster.h"
#include "TClonesArray.h"
#include "THcSignalHit.h"
#include "THcHodoHit.h"
@@ -34,7 +35,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
const Int_t planenum,
THaDetectorBase* parent )
: THaSubDetector(name,description,parent),
- fParentHitList(0), frPosAdcErrorFlag(0), frNegAdcErrorFlag(0),
+ fParentHitList(0), fCluster(0),frPosAdcErrorFlag(0), frNegAdcErrorFlag(0),
frPosTDCHits(0), frNegTDCHits(0), frPosADCHits(0), frNegADCHits(0),
frPosADCSums(0), frNegADCSums(0), frPosADCPeds(0), frNegADCPeds(0),
fHodoHits(0), frPosTdcTimeRaw(0), frPosAdcPedRaw(0),
@@ -61,6 +62,8 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
// Normal constructor with name and description
fHodoHits = new TClonesArray("THcHodoHit",16);
+ fCluster = new TClonesArray("THcScintPlaneCluster", 10);
+
frPosAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
frNegAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
@@ -97,6 +100,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
frNegAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
frNegAdcPulseTime = new TClonesArray("THcSignalHit", 16);
+
fPlaneNum = planenum;
fTotPlanes = planenum;
fNScinHits = 0;
@@ -111,6 +115,8 @@ THcScintillatorPlane::~THcScintillatorPlane()
delete frPosAdcErrorFlag; frPosAdcErrorFlag = NULL;
delete frNegAdcErrorFlag; frNegAdcErrorFlag = NULL;
+ delete fCluster; fCluster = NULL;
+
delete fHodoHits;
delete frPosTDCHits;
delete frNegTDCHits;
@@ -384,6 +390,7 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
// Create arrays to hold results here
InitializePedestals();
+
fNumGoodPosAdcHits = vector<Int_t> (fNelem, 0.0);
fNumGoodNegAdcHits = vector<Int_t> (fNelem, 0.0);
fNumGoodPosTdcHits = vector<Int_t> (fNelem, 0.0);
@@ -527,7 +534,12 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode )
{"DiffDisTrackCorr", "TW Corrected Dist Difference between track and scintillator position (cm)", "fGoodDiffDistTrack"},
{"TrackXPos", "Track X position at plane (cm)", "fTrackXPosition"},
{"TrackYPos", "Track Y position at plane (cm)", "fTrackYPosition"},
- //{"ngoodhits", "Number of paddle hits (passed tof tolerance and used to determine the focal plane time )", "GetNGoodHits() "},
+ {"NumClus", "Number of clusters", "fNumberClusters"},
+ {"Clus.Pos", "Position of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterPosition()"},
+ {"Clus.Size", "Size of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterSize()"},
+ {"Clus.Flag", "Flag of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterFlag()"},
+ {"Clus.UsedFlag", "USed Flag of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterUsedFlag()"},
+ //{"ngoodhits", "Number of paddle hits (passed tof tolerance and used to determine the focal plane time )", "GetNGoodHits() "},
{ 0 }
};
@@ -543,6 +555,8 @@ void THcScintillatorPlane::Clear( Option_t* )
*/
//cout << " Calling THcScintillatorPlane::Clear " << GetName() << endl;
// Clears the hit lists
+ fCluster->Clear();
+
frPosAdcErrorFlag->Clear();
frNegAdcErrorFlag->Clear();
@@ -576,6 +590,7 @@ void THcScintillatorPlane::Clear( Option_t* )
frNegAdcPulseAmp->Clear();
frNegAdcPulseTime->Clear();
+
//Clear occupancies
for (UInt_t ielem = 0; ielem < fNumGoodPosAdcHits.size(); ielem++)
fNumGoodPosAdcHits.at(ielem) = 0;
@@ -630,6 +645,7 @@ void THcScintillatorPlane::Clear( Option_t* )
fHitDistance = kBig;
fTrackXPosition = kBig;
fTrackYPosition = kBig;
+ fNumberClusters=0;
}
//_____________________________________________________________________________
@@ -725,7 +741,6 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
frNegAdcPulseInt->Clear();
frNegAdcPulseAmp->Clear();
frNegAdcPulseTime->Clear();
-
//stripped
fNScinHits=0;
@@ -1031,7 +1046,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
// Do corrections if valid TDC on both ends of bar
- if(btdcraw_pos && btdcraw_neg && badcraw_pos && badcraw_neg) {
+ if( (btdcraw_pos && btdcraw_neg) && (badcraw_pos && badcraw_neg) ) {
// Do the pulse height correction to the time. (Position dependent corrections later)
Double_t timec_pos, timec_neg;
if(fTofUsingInvAdc) {
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index 62ef1761337367b7fec9e82446573291fac5fbb1..052ac9512ee1b7706824a5f0eb4562706a57fdbf 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -14,6 +14,7 @@
#include "THaSubDetector.h"
#include "TClonesArray.h"
+#include "THcScintPlaneCluster.h"
using namespace std;
@@ -57,22 +58,31 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t GetPosOffset() {return fPosOffset;};
Double_t GetPosCenter(Int_t PaddleNo) {return fPosCenter[PaddleNo];}; // counting from zero!
Double_t GetFpTime() {return fFptime;};
+ Double_t GetNumberClusters() {return fNumberClusters;}; // number of paddle clusters
void SetFpTime(Double_t f) {fFptime=f;};
void SetNGoodHits(Int_t ng) {fNGoodHits=ng;};
void SetHitDistance(Double_t f) {fHitDistance=f;}; // Distance between track and hit paddle
void SetTrackXPosition(Double_t f) {fTrackXPosition=f;}; // Distance track X position at plane
void SetTrackYPosition(Double_t f) {fTrackYPosition=f;}; // Distance track Y position at plane
+ void SetNumberClusters(Int_t nclus) {fNumberClusters=nclus;}; // number of paddle
+ void SetCluster(Int_t ic,Double_t pos) {((THcScintPlaneCluster*) fCluster->ConstructedAt(ic))->Set(ic,pos);}
+ void SetClusterSize(Int_t ic,Double_t size) {((THcScintPlaneCluster*) fCluster->ConstructedAt(ic))->SetSize(size);}
+ void SetClusterFlag(Int_t ic,Double_t flag) {((THcScintPlaneCluster*) fCluster->ConstructedAt(ic))->SetFlag(flag);}
+ void SetClusterUsedFlag(Int_t ic,Double_t flag) {((THcScintPlaneCluster*) fCluster->ConstructedAt(ic))->SetUsedFlag(flag);}
TClonesArray* fParentHitList;
TClonesArray* GetHits() { return fHodoHits;};
+ TClonesArray* fCluster;
+
protected:
TClonesArray* frPosAdcErrorFlag;
TClonesArray* frNegAdcErrorFlag;
+
TClonesArray* frPosTDCHits;
TClonesArray* frNegTDCHits;
TClonesArray* frPosADCHits;
@@ -113,6 +123,7 @@ class THcScintillatorPlane : public THaSubDetector {
Int_t fTotNumPosAdcHits;
Int_t fTotNumNegAdcHits;
Int_t fTotNumAdcHits;
+ Int_t fNumberClusters;
Int_t fTotNumPosTdcHits;
Int_t fTotNumNegTdcHits;