GN: Halfway through the implementation of h_trans_scin and getting a

hodoscope start time.

src/THcHodoscope.cxx

@@ -349,26 +349,25 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
 // need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
   fMaxHodoScin=fMaxScinPerPlane*fNPlanes; 
   if (fDebug>=1)  cout <<"fMaxScinPerPlane = "<<fMaxScinPerPlane<<" fMaxHodoScin = "<<fMaxHodoScin<<endl;
-
-  Double_t fHodoVelLight[fMaxHodoScin];
-  Double_t fHodoPosSigma[fMaxHodoScin];
-  Double_t fHodoNegSigma[fMaxHodoScin];
-  Double_t fHodoPosMinPh[fMaxHodoScin];
-  Double_t fHodoNegMinPh[fMaxHodoScin];
-  Double_t fHodoPosPhcCoeff[fMaxHodoScin];
-  Double_t fHodoNegPhcCoeff[fMaxHodoScin];
-  Double_t fHodoPosTimeOffset[fMaxHodoScin];
-  Double_t fHodoNegTimeOffset[fMaxHodoScin];
-  Int_t fHodoPosPedLimit[fMaxHodoScin];
-  Int_t fHodoNegPedLimit[fMaxHodoScin];
-  Int_t fTofUsingInvAdc;
-  Double_t fHodoPosInvAdcOffset[fMaxHodoScin];
-  Double_t fHodoNegInvAdcOffset[fMaxHodoScin];
-  Double_t fHodoPosInvAdcLinear[fMaxHodoScin];
-  Double_t fHodoNegInvAdcLinear[fMaxHodoScin];
-  Double_t fHodoPosInvAdcAdc[fMaxHodoScin];
-  Double_t fHodoNegInvAdcAdc[fMaxHodoScin];
-
+  
+  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];
+  
   prefix[1]='\0';
   DBRequest list[]={
     {"start_time_center", &fStartTimeCenter, kDouble},
@@ -407,7 +406,7 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
     cout <<"TdcMin = "<<fScinTdcMin<<" TdcMax = "<<fScinTdcMax<<endl;
     cout <<"TofTolerance = "<<fTofTolerance<<endl;
     cout <<"*** VelLight ***\n";
-    for (int i1=1;i1<=fNPlanes;i1++) {
+    for (int i1=0;i1<fNPlanes;i1++) {
       cout<<"Plane "<<i1<<endl;
       for (int i2=0;i2<fMaxScinPerPlane;i2++) {
 	cout<<fHodoVelLight[GetScinIndex(i1,i2)]<<" ";
@@ -416,6 +415,22 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
     }
  
     cout <<endl<<endl;
+    // check fHodoPosPhcCoeff
+    /*
+    cout <<"fHodoPosPhcCoeff = ";
+    for (int i1=0;i1<fMaxHodoScin;i1++) {
+      cout<<this->GetHodoPosPhcCoeff(i1)<<" ";
+    }
+    cout<<endl;
+    */
+  }
+  //
+  if ((fTofTolerance > 0.5) && (fTofTolerance < 10000.)) {
+    cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
+  }
+  else {
+    fTofTolerance= 3.0;
+    cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
   }
   fIsInit = true;
   return kOK;
@@ -492,7 +507,24 @@ void THcHodoscope::DeleteArrays()
 {
   // Delete member arrays. Used by destructor.
 
-  delete [] fNPaddle; fNPaddle = 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 [] fSpacing; fSpacing = NULL;
   //delete [] fCenter;  fCenter = NULL; // This 2D. What is correct way to delete?
 
@@ -561,7 +593,10 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
     //    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,fScinTdcMin,fScinTdcMax,nexthit);
+    nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
+        if (fPlanes[ip]->GetNScinHits()>0) {
+      fPlanes[ip]->PulseHeightCorrection();
+      }
   }
 
   // fRawHitList is TClones array of THcHodoscopeHit objects
@@ -585,7 +620,6 @@ Int_t THcHodoscope::ApplyCorrections( void )
 {
   return(0);
 }
-
 //_____________________________________________________________________________
 Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle,
 					     const ESide side)
@@ -603,7 +637,8 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray&  tracks  )
   // Apply corrections and reconstruct the complete hits.
   //
   //  static const Double_t sqrt2 = TMath::Sqrt(2.);
-  /*  cout <<"**** in THcHodoscope CoarseProcess ********\n";  
+  //  cout <<"**** in THcHodoscope CoarseProcess ********\n"; 
+  /*  
   for(Int_t i=0;i<fNPlanes;i++) {
     cout<<i<<" ";
     fPlanes[i]->CoarseProcess(tracks);
@@ -629,10 +664,11 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
 }
 //_____________________________________________________________________________
 Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
-  // Return the index of a scintillator given the plane # and the paddle #
-  // This assumes that planes start counting from 1
-  // and paddles start counting from 0! (not ideal but that's what I have for now)
-  return fNPlanes*nPaddle+nPlane-1;
+  // 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;
 }
 //_____________________________________________________________________________
 Int_t THcHodoscope::GetScinIndex( Int_t nSide, Int_t nPlane, Int_t nPaddle ) {

src/THcHodoscope.h

@@ -30,14 +30,24 @@ public:
   virtual Int_t      FineProcess( TClonesArray& tracks );
   
   virtual Int_t      ApplyCorrections( void );
-
   //  Int_t GetNHits() const { return fNhit; }
   Double_t GetStartTime() const { return fStartTime; }
   Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle);
   Int_t GetScinIndex(Int_t nSide, Int_t nPlane, Int_t nPaddle);
   Double_t GetPathLengthCentral();
-
   Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
+  Double_t GetTdcMin() const {return fScinTdcMin;}
+  Double_t GetTdcMax() const {return fScinTdcMax;}
+  Double_t GetTofTolerance() const {return fTofTolerance;}
+  Double_t GetTdcToTime() const {return fScinTdcToTime;}
+  Double_t GetHodoPosPhcCoeff(Int_t iii) const {return fHodoPosPhcCoeff[iii];}
+  Double_t GetHodoNegPhcCoeff(Int_t iii) const {return fHodoNegPhcCoeff[iii];}
+  Double_t GetHodoPosMinPh(Int_t iii) const {return fHodoPosMinPh[iii];}
+  Double_t GetHodoNegMinPh(Int_t iii) const {return fHodoNegMinPh[iii];}
+  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];}
+
   const TClonesArray* GetTrackHits() const { return fTrackProj; }
   
   friend class THaScCalib;
@@ -53,6 +63,8 @@ protected:
 
   Double_t fStartTime;
   
+  // Per-event data
+
   // Potential Hall C parameters.  Mostly here for demonstration
   Int_t fNPlanes,fMaxScinPerPlane,fMaxHodoScin; // number of planes; max number of scin/plane; product of the first two 
   Double_t fStartTimeCenter, fStartTimeSlop, fScinTdcToTime;

src/THcScintillatorPlane.cxx

@@ -5,7 +5,7 @@
 // THcScintillatorPlane
 //
 //////////////////////////////////////////////////////////////////////////
-
+#include "TMath.h"
 #include "THcScintillatorPlane.h"
 #include "TClonesArray.h"
 #include "THcSignalHit.h"
@@ -38,6 +38,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
   fNegADCHits = new TClonesArray("THcSignalHit",16);
   fPlaneNum = planenum;
   fTotPlanes = planenum;
+  fNScinHits = 0;
 }
 //______________________________________________________________________________
 THcScintillatorPlane::THcScintillatorPlane( const char* name, 
@@ -54,6 +55,7 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name,
   fNegADCHits = new TClonesArray("THcSignalHit",16);
   fPlaneNum = planenum;
   fTotPlanes = totplanes;
+  fNScinHits = 0;
 }
 
 //______________________________________________________________________________
@@ -129,14 +131,16 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
   strcat(parname,"_nr");
   fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer();
   //
+  // Based on the signs of these quantities in the .pos file the correspondence 
+  // should be bot=>left  and top=>right when comparing x and y-type scintillators
   char *tmpleft, *tmpright;
   if (fPlaneNum==1 || fPlaneNum==3) {
     tmpleft="left";
     tmpright="right";
   } 
   else {
-    tmpleft="top";
-    tmpright="bot";
+    tmpleft="bot";
+    tmpright="top";
   }
 
   Double_t tmpdouble[fTotPlanes];
@@ -168,24 +172,6 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
    cout <<"PosOffset = "<<fPosOffset<<endl;
    cout <<"PosCenter[0] = "<<fPosCenter[0]<<endl;
 
-  /* for(Int_t i=0;i<fNPlanes;i++) {
-    fCenter[i] = new Double_t[fNPaddle[i]];
-    DBRequest list[]={
-      {Form("scin_%s_spacing",fPlaneNames[i]), &fSpacing[i], kDouble},
-      {Form("scin_%s_center",fPlaneNames[i]), fCenter[i], kDouble, fNPaddle[i]},
-      {0}
-    };
-    gHcParms->LoadParmValues((DBRequest*)&list,prefix);
-  }
-  */
-  ///  Int_t tmpint = new Int_t [GetParent()->GetNPlanes()]; 
-  /// cout<<"Nplanes in scintplane!! = "<<fTotPlanes <<endl;
-
-  // First letter of GetParent()->GetPrefix() tells us what prefix to
-  // use on parameter names.  
-
-  //  Find the number of elements
-  
   // Think we will make special methods to pass most
   // How generic do we want to make this class?  
   // The way we get parameter data is going to be pretty specific to
@@ -262,78 +248,23 @@ Int_t THcScintillatorPlane::FineProcess( TClonesArray& tracks )
 {
   return 0;
 }
-
 //_____________________________________________________________________________
 Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
-{
-  // Extract the data for this plane from hit list
-  // Assumes that the hit list is sorted by plane, so we stop when the
-  // plane doesn't agree and return the index for the next hit.
-
-  Int_t nPosTDCHits=0;
-  Int_t nNegTDCHits=0;
-  Int_t nPosADCHits=0;
-  Int_t nNegADCHits=0;
-  fPosTDCHits->Clear();
-  fNegTDCHits->Clear();
-  fPosADCHits->Clear();
-  fNegADCHits->Clear();
-
-  Int_t nrawhits = rawhits->GetLast()+1;
-  // cout << "THcScintillatorPlane::ProcessHits " << fPlaneNum << " " << nexthit << "/" << nrawhits << endl;
-
-  Int_t ihit = nexthit;
-  while(ihit < nrawhits) {
-    THcHodoscopeHit* hit = (THcHodoscopeHit *) rawhits->At(ihit);
-    if(hit->fPlane > fPlaneNum) {
-      break;
-    }
-
-    // TDC positive hit
-    if(hit->fTDC_pos >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
-      sighit->Set(hit->fCounter, hit->fTDC_pos);
-    }
-
-    // TDC negative hit
-    if(hit->fTDC_neg >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
-      sighit->Set(hit->fCounter, hit->fTDC_neg);
-    }
-
-    // ADC positive hit
-    if(hit->fADC_pos >  0) {
-      //     cout <<"adc pos hit!!\n";
-      THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
-      sighit->Set(hit->fCounter, hit->fADC_pos);
-    }
-
-    // ADC negative hit
-    if(hit->fADC_neg >  0) {   
-      // cout <<"adc neg hit!!\n";
-      THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
-      sighit->Set(hit->fCounter, hit->fADC_neg);
-    }
-
-    ihit++;
-  }
-  return(ihit);
-}
-
-//_____________________________________________________________________________
-Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Double_t mintdc, Double_t maxtdc,Int_t nexthit)
 {
   // Extract the data for this plane from hit list
   // Assumes that the hit list is sorted by plane, so we stop when the
   // plane doesn't agree and return the index for the next hit.
   // GN: Select only events that have at least one of their TDC signals in the
   // right range. 
-  // NOTE: 01/24/2013: Not quite subtracting pedestals for now! Also subtract the pedestals from the adc signal
+  // Also subtract the pedestals from the adc signal (as per the ENGINE) we're not checking here
+  // if the actual ADC is larger than the pedestal value we subtract!!
 
+  Double_t mintdc, maxtdc;
   Int_t nPosTDCHits=0;
   Int_t nNegTDCHits=0;
   Int_t nPosADCHits=0;
   Int_t nNegADCHits=0;
+  fNScinHits=0;
   fPosTDCHits->Clear();
   fNegTDCHits->Clear();
   fPosADCHits->Clear();
@@ -341,7 +272,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Double_t mintdc,
 
   Int_t nrawhits = rawhits->GetLast()+1;
   // cout << "THcScintillatorPlane::ProcessHits " << fPlaneNum << " " << nexthit << "/" << nrawhits << endl;
-
+  mintdc=((THcHodoscope *)GetParent())->GetTdcMin();
+  maxtdc=((THcHodoscope *)GetParent())->GetTdcMax();
   Int_t ihit = nexthit;
   while(ihit < nrawhits) {
     THcHodoscopeHit* hit = (THcHodoscopeHit *) rawhits->At(ihit);
@@ -352,6 +284,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Double_t mintdc,
     // check TDC values
     if (((hit->fTDC_pos >= mintdc) && (hit->fTDC_pos <= maxtdc)) ||
 	((hit->fTDC_neg >= mintdc) && (hit->fTDC_neg <= maxtdc))) {
+
       //TDC positive hit
       THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
       sighit->Set(hit->fCounter, hit->fTDC_pos);
@@ -359,25 +292,103 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Double_t mintdc,
       THcSignalHit *sighit2 = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
       sighit2->Set(hit->fCounter, hit->fTDC_neg);
       // ADC positive hit
-      if(hit->fADC_pos >  0) {
-	THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
-	sighit->Set(hit->fCounter, hit->fADC_pos);
-      }
+      /// if(hit->fADC_pos >  0) {
+	THcSignalHit *sighit3 = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+	sighit3->Set(hit->fCounter, hit->fADC_pos-fPosPed[ihit]);
+	///} else {
+	///  cout<<"Skipping ADC_pos "<<hit->fADC_pos<<endl;
+	/// }
       // ADC negative hit
-      if(hit->fADC_neg >  0) {   
+	///      if(hit->fADC_neg >  0) {   
 	// cout <<"adc neg hit!!\n";
-	THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
-	sighit->Set(hit->fCounter, hit->fADC_neg);
-      }
-      else {
-	cout <<"skipping BAD tdc event\n";
-      }
+	THcSignalHit *sighit4 = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+	sighit4->Set(hit->fCounter, hit->fADC_neg-fNegPed[ihit]);
+	///} else {
+	///cout<<"Skipping ADC_neg "<<hit->fADC_neg<<endl;
+	///}
+      //      cout <<"test "<<fNHits<<endl;
+      fNScinHits=fNScinHits++;
+    }
+    else {
+      //cout <<"pos TDC "<<hit->fTDC_pos<<" "<<mintdc<<" "<<maxtdc<<endl;
+      //cout <<"neg TDC "<<hit->fTDC_neg<<" "<<mintdc<<" "<<maxtdc<<endl;
+      //cout <<"skipping BAD tdc event\n";
     }
     ihit++;
   }
+
   return(ihit);
 }
+//________________________________________________________________________________
 
+Int_t THcScintillatorPlane::PulseHeightCorrection()
+{
+  // Perform pulse height correction of the TDC values as in the original h_trans_scin
+  // see original comments below
+  /*! Calculate all corrected hit times and histogram
+    ! This uses a copy of code below. Results are save in time_pos,neg
+    ! including the z-pos. correction assuming nominal value of betap
+    ! Code is currently hard-wired to look for a peak in the
+    ! range of 0 to 100 nsec, with a group of times that all
+    ! agree withing a time_tolerance of time_tolerance nsec. The normal
+    ! peak position appears to be around 35 nsec (SOS0 or 31 nsec (HMS)
+    ! NOTE: if want to find particles with beta different than
+    !       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?
+  Double_t dist_from_center,scint_center,hit_position,time_pos[100],time_neg[100],hbeta_pcent;
+
+  mintdc=((THcHodoscope *)GetParent())->GetTdcMin();
+  maxtdc=((THcHodoscope *)GetParent())->GetTdcMax();
+  //  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) &&
+	(((THcSignalHit*) fNegTDCHits->At(i))->GetData()>=mintdc) &&
+	(((THcSignalHit*) fNegTDCHits->At(i))->GetData()<=maxtdc)) {
+	  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);
+	  //	  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();
+	  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);
+	  //	  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++) {
+    }
+  }
+  return 0;
+}
 //_____________________________________________________________________________
 Int_t THcScintillatorPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
 {
@@ -476,6 +487,10 @@ void THcScintillatorPlane::InitializePedestals( )
 Int_t THcScintillatorPlane::GetNelem() 
 {
   return fNelem;
+}//____________________________________________________________________________
+Int_t THcScintillatorPlane::GetNScinHits() 
+{
+  return fNScinHits;
 }
 //____________________________________________________________________________
 Double_t THcScintillatorPlane::GetSpacing()

src/THcScintillatorPlane.h

@@ -37,11 +37,13 @@ class THcScintillatorPlane : public THaSubDetector {
   virtual Bool_t   IsPid()      { return kFALSE; }
 
   virtual Int_t ProcessHits(TClonesArray* rawhits, Int_t nexthit);
-  virtual Int_t ProcessHits(TClonesArray* rawhits, Double_t mintdc, Double_t maxtdc, Int_t nexthit);
+  virtual Int_t PulseHeightCorrection();
+
   virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit);
   virtual void  CalculatePedestals( );
 
   Int_t GetNelem(); // return number of paddles in this plane
+  Int_t GetNScinHits(); // return the number of hits in this plane (that pass min/max TDC cuts)
   Double_t GetSpacing(); // return the spacing of paddles 
   Double_t GetSize();    // return paddle size
   Double_t GetHodoSlop(); // return the hodo slop
@@ -62,12 +64,13 @@ class THcScintillatorPlane : public THaSubDetector {
   TClonesArray* fNegADCHits;
 
   Int_t fPlaneNum;		/* Which plane am I 1-4 */
-  Int_t fTotPlanes;              /* so we can read variables that are not indexed by plane id */
+  Int_t fTotPlanes;             /* so we can read variables that are not indexed by plane id */
   Int_t fNelem;			/* Need since we don't inherit from 
 				 detector base class */
+  Int_t fNScinHits;                 /* Number of hits in this plane */
   Double_t fSpacing;            /* paddle spacing */
   Double_t fSize;               /* paddle size */
-  Double_t fZpos;                /* z position */
+  Double_t fZpos;               /* z position */
   Double_t fDzpos;
   Double_t fHodoSlop;           /* max allowed slop for this plane */
   Double_t fPosLeft;            /* NOTE: "left" = "top" for a Y scintillator */
@@ -75,6 +78,8 @@ class THcScintillatorPlane : public THaSubDetector {
   Double_t fPosOffset;
   Double_t fPosCenter[16];         /* array with centers for all scintillators in the plane */
 
+
+  Double_t fTolerance; /* need this for PulseHeightCorrection */
   /* Pedestal Quantities */
   Int_t fNPedestalEvents;	/* Number of pedestal events */
   Int_t fMinPeds;		/* Only analyze/update if num events > */