diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index 16c30403e4f8a637a2deea47cee76d39fecf6be8..297e51b8b210cdcdb049229e98ce469b2025fe8a 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -259,8 +259,10 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
// Int_t plen=strlen(parname);
cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;
+ fBetaNoTrk = 0.;
+ fBetaNoTrkChiSq = 0.;
+
fNPaddle = new UInt_t [fNPlanes];
-
fFPTime = new Double_t [fNPlanes];
fPlaneCenter = new Double_t[fNPlanes];
fPlaneSpacing = new Double_t[fNPlanes];
@@ -430,6 +432,8 @@ Int_t THcHodoscope::DefineVariables( EMode mode )
RVarDef vars[] = {
// Move these into THcHallCSpectrometer using track fTracks
+ {"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
+ {"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
{"fpHitsTime", "Time at focal plane from all hits", "fFPTime"},
{"starttime", "Hodoscope Start Time", "fStartTime"},
{"goodstarttime", "Hodoscope Good Start Time", "fGoodStartTime"},
@@ -504,12 +508,9 @@ void THcHodoscope::DeleteArrays()
inline
void THcHodoscope::ClearEvent()
{
- // Reset per-event data.
- // for ( Int_t imaxhit = 0; imaxhit < MAXHODHITS; imaxhit++ ){
- // fBeta[imaxhit] = 0.;
- // fBetaChisq[imaxhit] = 0.;
- // }
+ fBetaNoTrk = 0.0;
+ fBetaNoTrkChiSq = 0.0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip]->Clear();
@@ -590,6 +591,7 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
//_____________________________________________________________________________
void THcHodoscope::EstimateFocalPlaneTime( void )
{
+
Int_t timehist[200];
for (Int_t i=0;i<200;i++) {
@@ -629,17 +631,27 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
ihit = 0;
Double_t fpTimeSum = 0.0;
+ Int_t jhit = 0;
fNfptimes=0;
+ fNoTrkPlaneInfo.clear();
+ fNoTrkHitInfo.clear();
for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fNoTrkPlaneInfo.push_back(NoTrkPlaneInfo());
+ fNoTrkPlaneInfo[ip].goodplanetime = kFALSE;
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for(Int_t i=0;i<nphits;i++) {
+ fNoTrkHitInfo.push_back(NoTrkHitInfo());
+ fNoTrkHitInfo[jhit].goodtwotimes = kFALSE;
+ fNoTrkHitInfo[jhit].goodscintime = kFALSE;
Double_t tmin = 0.5*binmax;
Double_t postime=((THcHodoHit*) hodoHits->At(i))->GetPosTOFCorrectedTime();
Double_t negtime=((THcHodoHit*) hodoHits->At(i))->GetNegTOFCorrectedTime();
if ((postime>tmin) && (postime<tmin+fTofTolerance) &&
(negtime>tmin) && (negtime<tmin+fTofTolerance)) {
+ fNoTrkHitInfo[jhit].goodtwotimes = kTRUE;
+ fNoTrkHitInfo[jhit].goodscintime = kTRUE;
// Both tubes fired
Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;
Double_t fptime = ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime()
@@ -649,11 +661,14 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
// Should also fill the all FP times histogram
fpTimeSum += fptime;
fNfptimes++;
+ fNoTrkPlaneInfo[ip].goodplanetime = kTRUE;
}
}
+ jhit++;
}
ihit++;
}
+
if(fNfptimes>0) {
fStartTime = fpTimeSum/fNfptimes;
fGoodStartTime=kTRUE;
@@ -661,6 +676,89 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
fStartTime = fStartTimeCenter;
fGoodStartTime=kFALSE;
}
+
+
+ if ( ( fNoTrkPlaneInfo[0].goodplanetime || fNoTrkPlaneInfo[1].goodplanetime ) &&
+ ( fNoTrkPlaneInfo[2].goodplanetime || fNoTrkPlaneInfo[3].goodplanetime ) ){
+
+ Double_t sumW = 0.;
+ Double_t sumT = 0.;
+ Double_t sumZ = 0.;
+ Double_t sumZZ = 0.;
+ Double_t sumTZ = 0.;
+ Int_t ihhit = 0;
+
+ for(Int_t ip=0;ip<fNPlanes;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 ( fNoTrkHitInfo[ihhit].goodscintime ) {
+
+ Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
+ TMath::Power( fHodoNegSigma[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;
+
+ 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
+
+ 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 ) ) {
+
+ fBetaNoTrk = tmp / tmpDenom;
+ fBetaNoTrkChiSq = 0.;
+ ihhit = 0;
+
+ for (Int_t ip = 0; ip < fNPlanes; 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;
+
+ if ( fNoTrkHitInfo[ihhit].goodscintime ) {
+
+ Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
+ Double_t timeDif = ( ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0 );
+ Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
+ TMath::Power( fHodoNegSigma[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
+
+ Double_t pathNorm = 1.0;
+
+ fBetaNoTrk = fBetaNoTrk * pathNorm;
+ fBetaNoTrk = fBetaNoTrk / 29.979; // velocity / c
+
+ } // condition for fTmpDenom
+ else {
+ fBetaNoTrk = 0.;
+ fBetaNoTrkChiSq = -2.;
+ } // else condition for fTmpDenom
+ }
+
}
//_____________________________________________________________________________
Int_t THcHodoscope::ApplyCorrections( void )
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index efc502bc3f34dd040ce714972d4ebef5f8ddebf9..514d04658fcb22e341ea9e63c1665452e3f3c8c4 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -69,7 +69,7 @@ public:
Double_t GetHodoVelLight(Int_t iii) const {return fHodoVelLight[iii];}
Double_t GetStartTimeCenter() const {return fStartTimeCenter;}
Double_t GetStartTimeSlop() const {return fStartTimeSlop;}
- Double_t GetBetaNotrk() const {return fBetaNotrk;}
+ Double_t GetBetaNotrk() const {return fBetaNoTrk;}
Int_t GetGoodRawPad(Int_t iii){return fTOFCalc[iii].good_raw_pad;}
Int_t GetGoodRawPlane(Int_t iii){return fTOFCalc[iii].pindex;}
@@ -107,7 +107,8 @@ protected:
Double_t fStartTime;
Int_t fNfptimes;
- Double_t fBetaNotrk;
+ Double_t fBetaNoTrk;
+ Double_t fBetaNoTrkChiSq;
// Per-event data
// Potential Hall C parameters. Mostly here for demonstration
@@ -223,6 +224,22 @@ protected:
// Double_t* gain;
// } fDataDest[NDEST]; // Lookup table for decoder
+ // Inforamtion for each plane
+ struct NoTrkPlaneInfo {
+ Bool_t goodplanetime;
+ NoTrkPlaneInfo () : goodplanetime(kFALSE) {}
+ };
+ std::vector<NoTrkPlaneInfo> fNoTrkPlaneInfo;
+
+ // Inforamtion for each plane
+ struct NoTrkHitInfo {
+ Bool_t goodtwotimes;
+ Bool_t goodscintime;
+ NoTrkHitInfo () : goodtwotimes(kFALSE) {}
+ };
+ std::vector<NoTrkHitInfo> fNoTrkHitInfo;
+
+
// Used in TOF calculation (FineProcess) to hold information about hits
// within a given plane
struct TOFPInfo {