Skip to content
Snippets Groups Projects
Commit b7ce84cf authored by Gabriel Niculescu's avatar Gabriel Niculescu
Browse files

hodoscope start time calculation. verified against the engine output.

parent 4300132a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/THcHodoscope.cxx
+ 79
− 9
View file @ b7ce84cf
......@@ -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;
}
////////////////////////////////////////////////////////////////////////////////
src/THcHodoscope.h
+ 13
− 2
View file @ b7ce84cf
......@@ -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
......
src/THcScintillatorPlane.cxx
+ 50
− 5
View file @ b7ce84cf
......@@ -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)
////////////////////////////////////////////////////////////////////////////////
src/THcScintillatorPlane.h
+ 13
− 3
View file @ b7ce84cf
......@@ -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 );
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment