diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index a2c401064566f04d9b2771af18d799d2f856e1c0..7fa88fa8149d920c6ee0c0fd1c59284adb56c31e 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -597,6 +597,9 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
// Get the Hall C style hitlist (fRawHitList) for this event
Int_t nhits = THcHitList::DecodeToHitList(evdata);
+ //
+ // GN: print event number so we can cross-check with engine
+ // if (evdata.GetEvNum()>1000) cout <<"hcana event no = "<<evdata.GetEvNum()<<endl;
if(gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
@@ -615,9 +618,9 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
// Let each plane get its hits
Int_t nexthit = 0;
- Int_t nfptimes=0;
fStartTime=0;
+ fNfptimes=0;
for(Int_t ip=0;ip<fNPlanes;ip++) {
// nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
// GN: select only events that have reasonable TDC values to start with
@@ -625,21 +628,27 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
if (fPlanes[ip]->GetNScinHits()>0) {
fPlanes[ip]->PulseHeightCorrection();
- if (TMath::Abs(fPlanes[ip]->GetFpTime()-fStartTimeCenter)<=fStartTimeSlop) {
- fStartTime=fStartTime+fPlanes[ip]->GetFpTime();
- nfptimes++;
+ // GN: allow for more than one fptime per plane!!
+ for (Int_t i=0;i<fPlanes[ip]->GetNScinGoodHits();i++) {
+ if (TMath::Abs(fPlanes[ip]->GetFpTime(i)-fStartTimeCenter)<=fStartTimeSlop) {
+ fStartTime=fStartTime+fPlanes[ip]->GetFpTime(i);
+ // GN write stuff out so I can compare with engine
+ /// cout<<"hcana event= "<<evdata.GetEvNum()<<" fNfptimes= "<<fNfptimes<<" fptime= "<<fPlanes[ip]->GetFpTime(i)<<endl;
+ fNfptimes++;
+ }
}
}
}
- if (nfptimes>0) {
- fStartTime=fStartTime/nfptimes;
+ if (fNfptimes>0) {
+ fStartTime=fStartTime/fNfptimes;
fGoodStartTime=kTRUE;
} else {
fGoodStartTime=kFALSE;
fStartTime=fStartTimeCenter;
}
- fStartTime=32.; // mkj force to constant
- // cout <<" stats = "<<fGoodStartTime<<" "<<nfptimes<<" fStartTime = "<<fStartTime<<endl;
+
+/// fStartTime=32.; // mkj force to constant
+/// if (fGoodStartTime) cout <<"hcana event= "<<evdata.GetEvNum()<<" fNfptimes= "<<fNfptimes<<" fStartTime= "<<fStartTime<<endl<<endl;
// fRawHitList is TClones array of THcHodoscopeHit objects
#if 0
for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
@@ -650,7 +659,6 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
}
cout << endl;
#endif
-
/// fStartTime = 500; // Drift Chamber will need this
return nhits;
@@ -720,4 +728,66 @@ Double_t THcHodoscope::GetPathLengthCentral() {
return fPathLengthCentral;
}
ClassImp(THcHodoscope)
+//_____________________________________________________________________________
+Double_t THcHodoscope::CalcBeta(Int_t trackno) {
+ // clone of Engine's h_tof_fit. Original comments follow:
+ /*
+*-------------------------------------------------------------------
+* author: John Arrington
+* created: 2/22/94
+*
+* h_tof_fit fits the velocity of the paritcle from the corrected
+* times generated by h_tof.
+*
+* modifications:
+* $Log: h_tof_fit.f,v $
+* Revision 1.10 1996/09/04 13:36:24 saw
+* (JRA) Include actual beta in calculation of focal plane time.
+*
+* Revision 1.9 1995/05/22 19:39:29 cdaq
+* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
+*
+* Revision 1.8 1995/02/23 13:29:15 cdaq
+* (SAW) Cosmetic Changes
+*
+* Revision 1.7 1995/02/10 18:49:57 cdaq
+* (JRA) Add track index to hgood_scin_time
+*
+* Revision 1.6 1994/09/13 21:26:53 cdaq
+* (JRA) fix chisq calculation
+*
+* Revision 1.5 1994/07/13 15:05:08 cdaq
+* (SAW) Add abs around tmpdenom that I left out last update
+*
+* Revision 1.4 1994/07/11 18:34:35 cdaq
+* (JRA) Increase comparison of tmpdenom from 1e-15 to 1e-10
+*
+* Revision 1.3 1994/07/08 19:42:31 cdaq
+* (JRA) Change fit from velocity to beta. Bad fits give beta=0
+*
+* Revision 1.2 1994/06/14 04:53:41 cdaq
+* (DFG) Protect against divide by 0 in beta calc
+*
+* Revision 1.1 1994/04/13 16:29:15 cdaq
+* Initial revision
+*
+*-------------------------------------------------------------------
+ */
+
+ Double_t SumW, SumT, SumZ, SumZZ, SumTZ;
+ Double_t ScinWeight;
+ Double_t tmp, t0 ,tmpdenom;
+ Double_t PathNorm;
+ Double_t tmpbeta;
+ Int_t hit;
+
+ SumW = 0.;
+ SumT = 0.;
+ SumZ = 0.;
+ SumZZ = 0.;
+ SumTZ = 0.;
+ // for (int i=0;i<
+
+ return tmpbeta;
+}
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index 7d18bc8d808b8ea6cead145c8c615e79544b0504..11fe399e9d286bebfb58bf547a0fe640759133c3 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -32,6 +32,7 @@ public:
virtual Int_t ApplyCorrections( void );
Double_t GetStartTime() const { return fStartTime; }
Bool_t IsStartTimeGood() const {return fGoodStartTime;};
+ Int_t GetNfptimes() const {return fNfptimes;};
Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle);
Int_t GetScinIndex(Int_t nSide, Int_t nPlane, Int_t nPaddle);
Double_t GetPathLengthCentral();
@@ -47,6 +48,13 @@ public:
Double_t GetHodoPosTimeOffset(Int_t iii) const {return fHodoPosTimeOffset[iii];}
Double_t GetHodoNegTimeOffset(Int_t iii) const {return fHodoNegTimeOffset[iii];}
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 GetBeta() const {return fBeta;}
+ Double_t GetHodoPosSigma(Int_t iii) const {return fHodoPosSigma[iii];}
+ Double_t GetHodoNegSigma(Int_t iii) const {return fHodoNegSigma[iii];}
+ Double_t CalcBeta(Int_t trackno);
const TClonesArray* GetTrackHits() const { return fTrackProj; }
@@ -61,8 +69,11 @@ protected:
// Per-event data
Bool_t fGoodStartTime;
- Double_t fStartTime;
-
+ Double_t fStartTime;
+ Int_t fNfptimes;
+
+ Double_t fBetaNotrk;
+ Double_t fBeta;
// Per-event data
// Potential Hall C parameters. Mostly here for demonstration
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index c326fe0a179917e068f8cb951033a06e8d383211..b39ca7adb3c69750b323e86978db9fa2585ecddd 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -46,6 +46,10 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
//
fMaxHits=53;
fpTime = -1.e5;
+ fpTimes = new Double_t [fMaxHits];
+ fScinTime = new Double_t [fMaxHits];
+ fScinSigma = new Double_t [fMaxHits];
+ fScinZpos = new Double_t [fMaxHits];
}
//______________________________________________________________________________
THcScintillatorPlane::THcScintillatorPlane( const char* name,
@@ -70,7 +74,10 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
//
fMaxHits=53;
fpTime = -1.e5;
-
+ fpTimes = new Double_t [fMaxHits];
+ fScinTime = new Double_t [fMaxHits];
+ fScinSigma = new Double_t [fMaxHits];
+ fScinZpos = new Double_t [fMaxHits];
}
//______________________________________________________________________________
@@ -85,6 +92,10 @@ THcScintillatorPlane::~THcScintillatorPlane()
delete frNegTDCHits;
delete frPosADCHits;
delete frNegADCHits;
+ delete fpTimes;
+ delete fScinTime;
+ delete fScinSigma;
+ delete fScinZpos;
}
//______________________________________________________________________________
@@ -508,14 +519,29 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
scin_corrected_time[i]=0.5*(postime[i]+negtime[i]);
- fpTime=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
}
else { // only one tube fired
scin_corrected_time[i]=0.0; // not a very good "flag" but there is the logical two_good_hits...
+ // no fpTimes for U!
+ }
+ }
+ //start time calculation. assume xp=yp=0 radians. project all
+ //time values to focal plane. use average for start time.
+
+ fNScinGoodHits=0;
+ for (i=0;i<fNScinHits;i++) {
+ j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber()-1;
+ index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
+ if (two_good_times[i]) { // both tubes fired
+ fpTimes[fNScinGoodHits]=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
+ fScinTime[fNScinGoodHits]=scin_corrected_time[i];
+ fScinSigma[fNScinGoodHits]=TMath::Sqrt(((THcHodoscope *)GetParent())->GetHodoPosSigma(index)*((THcHodoscope *)GetParent())->GetHodoPosSigma(index)+((THcHodoscope *)GetParent())->GetHodoNegSigma(index)*((THcHodoscope *)GetParent())->GetHodoNegSigma(index)); // not ideal by any stretch!!!
+ fScinZpos[fNScinGoodHits]=fZpos+(j%2)*fDzpos; // see comment above
+ // h_rfptime(hscin_plane_num(ihit))=fptime
+ fNScinGoodHits++; // increment the number of good hits
}
}
- // Start time calculation, assume xp=yp=0 radians. project all
- // time values to focal plane. use average for start time.
+ CalcFpTime();
return 0;
}
//_____________________________________________________________________________
@@ -616,7 +642,8 @@ void THcScintillatorPlane::InitializePedestals( )
Int_t THcScintillatorPlane::GetNelem()
{
return fNelem;
-}//____________________________________________________________________________
+}
+//____________________________________________________________________________
Int_t THcScintillatorPlane::GetNScinHits()
{
return fNScinHits;
@@ -663,6 +690,24 @@ Double_t THcScintillatorPlane::GetPosCenter(Int_t PaddleNo) {
return fPosCenter[PaddleNo];
}
//____________________________________________________________________________
+Double_t THcScintillatorPlane::CalcFpTime()
+{
+ Double_t tmp=0;
+ Int_t i,counter=0;
+ for (i=0;i<fNScinHits;i++) {
+ if (TMath::Abs(fpTimes[i]-((THcHodoscope *)GetParent())->GetStartTimeCenter())<=((THcHodoscope *)GetParent())->GetStartTimeSlop()) {
+ tmp+=fpTimes[i];
+ counter++;
+ }
+ }
+ if (counter>0) {
+ fpTime=tmp/counter;
+ } else {
+ fpTime=((THcHodoscope *)GetParent())->GetStartTimeCenter();
+ }
+ return fpTime;
+}
+//____________________________________________________________________________
ClassImp(THcScintillatorPlane)
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index cb832fa7dd0975bf1eb0d064e6e6cb91a2e60cfb..2a5cec516dfe899045b00539f5787fb571d12c8b 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -53,7 +53,13 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t GetPosRight();
Double_t GetPosOffset();
Double_t GetPosCenter(Int_t PaddleNo); // here we're counting from zero!
- Double_t GetFpTime() { return fpTime;};
+ Double_t CalcFpTime();
+ Double_t GetFpTime() {return fpTime;};
+ Double_t GetFpTime(Int_t index) { return fpTimes[index];};
+ Double_t GetScinTime(Int_t index) { return fScinTime[index];};
+ Double_t GetScinSigma(Int_t index) { return fScinSigma[index];};
+ Double_t GetScinZpos(Int_t index) { return fScinZpos[index];};
+ Int_t GetNScinGoodHits() const {return fNScinGoodHits;};
TClonesArray* fParentHitList;
@@ -106,8 +112,12 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t *fNegThresh;
//
- Double_t fpTime; // the original code only has one fpTime per plane!
-
+ Int_t fNScinGoodHits; // number of hits for which both ends of the paddle fired in time!
+ Double_t fpTime; // the original code only has one fpTime per plane!
+ Double_t *fpTimes; // ... but also allows for more than one hit per plane
+ Double_t *fScinTime; // array of scintillator times (only filled for goodhits)
+ Double_t *fScinSigma; // errors for the above
+ Double_t *fScinZpos; // zpositions for the above
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );