From ddd3cf9ad0cadd07b4970aa70bdbcb0e67f2e8b7 Mon Sep 17 00:00:00 2001
From: Gabriel Niculescu <gabriel@jalb.org>
Date: Tue, 5 Mar 2013 11:59:41 -0500
Subject: [PATCH] G.N. 2013 Calculating pulse height correction for the
hodoscope. Calculating focal plane times for each scintillator plane - (in
THcScintillatorPlane) Calculating hodoscope start time (average of fp times)
- (in THcHodoscope)
---
src/THcHodoscope.cxx | 24 +++++--
src/THcHodoscope.h | 4 +-
src/THcScintillatorPlane.cxx | 132 +++++++++++++++++++++++++++++++----
src/THcScintillatorPlane.h | 9 ++-
4 files changed, 147 insertions(+), 22 deletions(-)
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index c5c1fd6..de6278f 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -43,7 +43,8 @@ THcHodoscope::THcHodoscope( const char* name, const char* description,
//fTrackProj = new TClonesArray( "THaTrackProj", 5 );
// Construct the planes
fNPlanes = 0; // No planes until we make them
-
+ fStartTime=-1e5;
+ fGoodStartTime=kFALSE;
}
//_____________________________________________________________________________
@@ -416,7 +417,6 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
}
cout <<endl;
}
-
cout <<endl<<endl;
// check fHodoPosPhcCoeff
/*
@@ -592,16 +592,30 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
// Let each plane get its hits
Int_t nexthit = 0;
+ Int_t nfptimes=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
// as per the Engine h_strip_scin.f
nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
- if (fPlanes[ip]->GetNScinHits()>0) {
+ if (fPlanes[ip]->GetNScinHits()>0) {
fPlanes[ip]->PulseHeightCorrection();
+ ///cout <<"Plane "<<ip<<" number of fpTimes = "<<fPlanes[ip]->GetFpTimeHits()<<endl;
+ for (Int_t ifptimes=0;ifptimes<fPlanes[ip]->GetFpTimeHits();ifptimes++) {
+ fStartTime=fStartTime+fPlanes[ip]->GetFpTime(ifptimes);
+ nfptimes++;
}
+ }
}
-
+ if (nfptimes>0) {
+ fStartTime=fStartTime/nfptimes;
+ fGoodStartTime=kTRUE;
+ } else {
+ fGoodStartTime=kFALSE;
+ fStartTime=fStartTimeCenter;
+ }
+ cout <<" stats = "<<fGoodStartTime<<" "<<nfptimes<<" fStartTime = "<<fStartTime<<endl;
// fRawHitList is TClones array of THcHodoscopeHit objects
#if 0
for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
@@ -613,7 +627,7 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
cout << endl;
#endif
- fStartTime = 500; // Drift Chamber will need this
+ /// fStartTime = 500; // Drift Chamber will need this
return nhits;
}
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index ab24f0e..7d18bc8 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -30,8 +30,8 @@ public:
virtual Int_t FineProcess( TClonesArray& tracks );
virtual Int_t ApplyCorrections( void );
- // Int_t GetNHits() const { return fNhit; }
Double_t GetStartTime() const { return fStartTime; }
+ Bool_t IsStartTimeGood() const {return fGoodStartTime;};
Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle);
Int_t GetScinIndex(Int_t nSide, Int_t nPlane, Int_t nPaddle);
Double_t GetPathLengthCentral();
@@ -60,7 +60,7 @@ protected:
// Calibration
// Per-event data
-
+ Bool_t fGoodStartTime;
Double_t fStartTime;
// Per-event data
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index b5fe6e3..f1ccda7 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -38,7 +38,11 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fNegADCHits = new TClonesArray("THcSignalHit",16);
fPlaneNum = planenum;
fTotPlanes = planenum;
- fNScinHits = 0;
+ fNScinHits = 0;
+ //
+ fMaxHits=53;
+ fpTimeHits=0;
+ fpTime = new Double_t [fMaxHits];
}
//______________________________________________________________________________
THcScintillatorPlane::THcScintillatorPlane( const char* name,
@@ -56,6 +60,11 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
fPlaneNum = planenum;
fTotPlanes = totplanes;
fNScinHits = 0;
+ //
+ fMaxHits=53;
+ fpTimeHits=0;
+ fpTime = new Double_t [fMaxHits];
+
}
//______________________________________________________________________________
@@ -66,8 +75,10 @@ THcScintillatorPlane::~THcScintillatorPlane()
delete fNegTDCHits;
delete fPosADCHits;
delete fNegADCHits;
-
+ delete fpTime;
}
+
+//______________________________________________________________________________
THaAnalysisObject::EStatus THcScintillatorPlane::Init( const TDatime& date )
{
// Extra initialization for scintillator plane: set up DataDest map
@@ -336,18 +347,37 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
! reference particle, need to make sure this is big enough
! to accomodate difference in TOF for other particles
! Default value in case user hasn't defined something reasonable */
- Int_t i,j,index,nfound=0;
- Double_t mintdc, maxtdc,tdctotime;
- Double_t pos_ph[53],neg_ph[53],postime[53],negtime[53]; // the 53 should go in a param file (was hmax_scin_hits originally)
- Bool_t keep_pos[53],keep_neg[53]; // are these all really needed?
+ Int_t i,j,index;
+ Double_t mintdc, maxtdc,tdctotime,toftolerance,tmin;
+ Double_t pos_ph[53],neg_ph[53],postime[53],negtime[53],scin_corrected_time[53]; // the 53 should go in a param file (was hmax_scin_hits originally)
+ Bool_t keep_pos[53],keep_neg[53],two_good_times[53]; // are these all really needed?
Double_t dist_from_center,scint_center,hit_position,time_pos[100],time_neg[100],hbeta_pcent;
+ Int_t timehist[200],jmax,maxhit,nfound=0; // This seems as a pretty old-fashioned way of doing things. Is there a better way?
+
+
+ // protect against spam events
+ if (fNScinHits>1000) {
+ cout <<"Too many hits "<<fNScinHits<<" in this event! Skipping!\n";
+ return -1;
+ }
+ // zero out histogram
+ for (i=0;i<200;i++) {
+ timehist[i]=0;
+ }
+ for (i=0;i<53;i++) {
+ keep_pos[i]=kFALSE;
+ keep_neg[i]=kFALSE;
+ two_good_times[i]=kFALSE;
+ }
mintdc=((THcHodoscope *)GetParent())->GetTdcMin();
maxtdc=((THcHodoscope *)GetParent())->GetTdcMax();
+ tdctotime=((THcHodoscope *)GetParent())->GetTdcToTime();
+ toftolerance=((THcHodoscope *)GetParent())->GetTofTolerance();
// hbeta_pcent=(TH((THcHodoscope *)GetParent())->GetParent()
// Horrible hack until I find out where to get the central beta from momentum!! GN
hbeta_pcent=0.99776;
- tdctotime=((THcHodoscope *)GetParent())->GetTdcToTime();
+
for (i=0;i<fNScinHits;i++) {
if ((((THcSignalHit*) fPosTDCHits->At(i))->GetData()>=mintdc) &&
(((THcSignalHit*) fPosTDCHits->At(i))->GetData()<=maxtdc) &&
@@ -360,7 +390,6 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1)));
postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index);
- // cout <<postime[i]<<endl;
neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData();
negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime;
j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber();
@@ -368,25 +397,100 @@ Int_t THcScintillatorPlane::PulseHeightCorrection()
negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)*
TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1)));
negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegTimeOffset(index);
- // cout <<"postime = "<<postime[i]<<" negtime = "<<negtime[i]<<endl;
// Find hit position. If postime larger, then hit was nearer negative side.
dist_from_center=0.5*(negtime[i]-postime[i])*((THcHodoscope *)GetParent())->GetHodoVelLight(index);
scint_center=0.5*(fPosLeft+fPosRight);
hit_position=scint_center+dist_from_center;
- // cout <<fPosLeft<<" "<<fPosRight<<" hit position = "<<hit_position<<" ";
hit_position=TMath::Min(hit_position,fPosLeft);
- // cout<<hit_position<<" ";
hit_position=TMath::Max(hit_position,fPosRight);
- //cout<<hit_position<<endl;
postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
time_pos[i]=postime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
time_neg[i]=negtime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
nfound++;
- cout <<"time pos/neg = "<<time_pos[i]<<" "<<time_neg[i]<<endl;
- // for (int k=0;k<200;k++) {
+ for (int k=0;k<200;k++) {
+ tmin=0.5*k;
+ if ((time_pos[i]> tmin) && (time_pos[i] < tmin+toftolerance)) {
+ timehist[k]++;
+ }
+ if ((time_neg[i]> tmin) && (time_neg[i] < tmin+toftolerance)) {
+ timehist[k]++;
+ }
+ }
+ nfound++;
+ }
+ }
+ // Find the bin with most hits
+ jmax=-1;
+ maxhit=0;
+ for (i=0;i<200;i++) {
+ if (timehist[i]>maxhit) {
+ jmax=i;
+ maxhit=timehist[i];
+ }
+ }
+ if (jmax>=0) {
+ tmin=0.5*jmax;
+ for (i=0;i<fNScinHits;i++) {
+ if ((time_pos[i]>tmin) && (time_pos[i]<tmin+toftolerance)) {
+ keep_pos[i]=kTRUE;
+ }
+ if ((time_neg[i]>tmin) && (time_neg[i]<tmin+toftolerance)) {
+ keep_neg[i]=kTRUE;
+ }
+ }
+ }
+ // Resume regular tof code, now using time filer(?) from above
+ // Check for TWO good TDC hits
+ for (i=0;i<fNScinHits;i++) {
+ if ((((THcSignalHit*) fPosTDCHits->At(i))->GetData()>=mintdc) &&
+ (((THcSignalHit*) fPosTDCHits->At(i))->GetData()<=maxtdc) &&
+ (((THcSignalHit*) fNegTDCHits->At(i))->GetData()>=mintdc) &&
+ (((THcSignalHit*) fNegTDCHits->At(i))->GetData()<=maxtdc) &&
+ keep_pos[i] && keep_neg[i]) {
+ two_good_times[i]=kTRUE;
+ }
+ } // end of loop that finds tube setting time
+ fpTimeHits=0;
+ for (i=0;i<fNScinHits;i++) {
+ if (two_good_times[i]) { // both tubes fired
+ // correct time for everything except veloc. correction in order
+ // to find hit location from difference in TDC.
+ pos_ph[i]=((THcSignalHit*) fPosADCHits->At(i))->GetData();
+ postime[i]=((THcSignalHit*) fPosTDCHits->At(i))->GetData()*tdctotime;
+ j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber();
+ index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
+ postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)*
+ TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1)));
+ postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index);
+ //
+ neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData();
+ negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime;
+ j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber();
+ index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j);
+ negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)*
+ TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1)));
+ negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegTimeOffset(index);
+ // find hit position. If postime larger, then hit was nearer negative side
+ dist_from_center=0.5*(negtime[i]-postime[i])*((THcHodoscope *)GetParent())->GetHodoVelLight(index);
+ scint_center=0.5*(fPosLeft+fPosRight);
+ hit_position=scint_center+dist_from_center;
+ hit_position=TMath::Min(hit_position,fPosLeft);
+ hit_position=TMath::Max(hit_position,fPosRight);
+ postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
+ negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index);
+ scin_corrected_time[i]=0.5*(postime[i]+negtime[i]);
+ fpTime[i]=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent);
+ /// cout <<"fptime ["<<i<<"]= "<<fpTime[i]<<endl;
+ fpTimeHits++;
+ }
+ else { // only one tube fired
+ scin_corrected_time[i]=0.0; // not a very good "flag" but there is the logical two_good_hits...
+ fpTime[i]=0.0; // NOTE: we're not incrementing fpTimeHits. In principle these two lines might be deleted
}
}
+ // Start time calculation, assume xp=yp=0 radians. project all
+ // time values to focal plane. use average for start time.
return 0;
}
//_____________________________________________________________________________
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index e3c5f4d..f47b1eb 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -53,6 +53,8 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t GetPosRight();
Double_t GetPosOffset();
Double_t GetPosCenter(Int_t PaddleNo); // here we're counting from zero!
+ Double_t GetFpTime(Int_t index) { return fpTime[index];};
+ Int_t GetFpTimeHits() { return fpTimeHits;};
TClonesArray* fParentHitList;
@@ -68,6 +70,7 @@ class THcScintillatorPlane : public THaSubDetector {
Int_t fNelem; /* Need since we don't inherit from
detector base class */
Int_t fNScinHits; /* Number of hits in this plane */
+ Int_t fMaxHits; /* maximum number of hits to be considered - useful for dimensioning arrays */
Double_t fSpacing; /* paddle spacing */
Double_t fSize; /* paddle size */
Double_t fZpos; /* z position */
@@ -98,7 +101,11 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t *fNegPed;
Double_t *fNegSig;
Double_t *fNegThresh;
-
+
+ //
+ Int_t fpTimeHits; // number of entries in *fpTime
+ Double_t *fpTime; // array with fpTimes for this scintillator plane
+
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
virtual void InitializePedestals( );
--
GitLab