From b625a5a7e977a4760bf3b728e5471aa965bba786 Mon Sep 17 00:00:00 2001
From: Mark Jones <jones@jlab.org>
Date: Thu, 2 Nov 2017 16:22:50 -0400
Subject: [PATCH] Update THcHodoscope, THcScintillatorPlane
THcHodoscope.cxx
1. Remove from the use of "temp_planes" for SHMS which
used only 3 planes for beta and determination of start time.
2. With change of Integral ADC to pC, change TMath::Sqrt(TMath::Max(20.0,adc_pos))
to Math::Sqrt(TMath::Max(20.0*.020,adc_pos)) to convert 20 channels to pC
3. Set fGoodFlags[itrack][ip][iphit].onTrack = kTRUE; when fTOFPInfo[ih].onTrack
Was always FALSE before.
4. Modify FineProcess
a. Add calculation of the track X and Y track position at each plane.
b. Calculate the difference between the track position and position
measured by hodoscope center. If multiple paddles hit then take
average of the paddles.
THcScintillatorPlane.h and cxx
1. add variables fHitDistance,fTrackXPosition and fTrackYPosition
2. add methods GetHitDistance(),GetTrackXPosition(),GetTrackYPosition()
3. add methods SetHitDistance,SetTrackXPosition,SetTrackYPosition
4. add variables to tree DiffDisTrack,TrackXPos and TrackYPos
5. set GoodPosAdcPulseTime to frPosAdcPulseTime
6. set fGoodNegAdcPulseTime to frNegAdcPulseTime
7. in HodoHit set PosADCtime and NegADCtime to PulseTime
---
src/THcHodoscope.cxx | 66 ++++++++++++++++++++++++------------
src/THcScintillatorPlane.cxx | 22 +++++++-----
src/THcScintillatorPlane.h | 10 +++++-
3 files changed, 68 insertions(+), 30 deletions(-)
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index 8d0fc57..164e95d 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -173,7 +173,7 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
// memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
fNScinHits = new Int_t [fNPlanes];
- fGoodPlaneTime = new Bool_t [fNPlanes];
+ fGoodPlaneTime = new Bool_t [fNPlanes];
fNPlaneTime = new Int_t [fNPlanes];
fSumPlaneTime = new Double_t [fNPlanes];
@@ -681,9 +681,7 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
Bool_t goodplanetime[fNPlanes];
Bool_t twogoodtimes[nscinhits];
- Int_t temp_planes=fNPlanes;
- if (fSHMS) temp_planes=3;
- for(Int_t ip=0;ip<temp_planes;ip++) {
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
goodplanetime[ip] = kFALSE;
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -747,9 +745,7 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
Double_t sumTZ = 0.;
Int_t ihhit = 0;
- temp_planes=fNPlanes;
- if (fSHMS) temp_planes=3;
- for(Int_t ip=0;ip<temp_planes;ip++) {
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -785,7 +781,7 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
fBetaNoTrkChiSq = 0.;
ihhit = 0;
- for (Int_t ip = 0; ip < temp_planes; ip++ ){ // Loop over planes
+ for (Int_t ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -845,8 +841,6 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// -------------------------------------------------
// fDumpOut << " ntrack = " << ntracks << endl;
- Int_t temp_planes=fNPlanes;
- if (fSHMS) temp_planes=3;
if (tracks.GetLast()+1 > 0 ) {
@@ -860,7 +854,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
if (!theTrack) return -1;
- for (Int_t ip = 0; ip < temp_planes; ip++ ){
+ for (Int_t ip = 0; ip < fNPlanes; ip++ ){
fGoodPlaneTime[ip] = kFALSE;
fNScinHits[ip] = 0;
fNPlaneTime[ip] = 0;
@@ -899,7 +893,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
fTOFPInfo.clear(); // SAW - combine these two?
Int_t ihhit = 0; // Hit # overall
- for(Int_t ip = 0; ip < temp_planes; ip++ ) {
+ for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
std::vector<GoodFlags> goodflagstmp2;
fGoodFlags[itrack].push_back(goodflagstmp2);
@@ -980,7 +974,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
timep -= fHodoPosInvAdcOffset[fPIndex]
+ pathp/fHodoPosInvAdcLinear[fPIndex]
+ fHodoPosInvAdcAdc[fPIndex]
- /TMath::Sqrt(TMath::Max(20.0,adc_pos));
+ /TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
} else {
timep -= fHodoPosPhcCoeff[fPIndex]*
TMath::Sqrt(TMath::Max(0.0,adc_pos/fHodoPosMinPh[fPIndex]-1.0))
@@ -1008,7 +1002,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
timen -= fHodoNegInvAdcOffset[fPIndex]
+ pathn/fHodoNegInvAdcLinear[fPIndex]
+ fHodoNegInvAdcAdc[fPIndex]
- /TMath::Sqrt(TMath::Max(20.0,adc_neg));
+ /TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
} else {
timen -= fHodoNegPhcCoeff[fPIndex]*
TMath::Sqrt(TMath::Max(0.0,adc_neg/fHodoNegMinPh[fPIndex]-1.0))
@@ -1094,6 +1088,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
Int_t fPIndex = GetScinIndex(ip,paddle);
if (fTOFPInfo[ih].onTrack) {
+ fGoodFlags[itrack][ip][iphit].onTrack = kTRUE;
if ( fTOFPInfo[ih].keep_pos ) { // 301
fTOFCalc[ih].good_tdc_pos = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
@@ -1197,7 +1192,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
} // Second loop over hits of a scintillator plane ends here
theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
- theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
+ if (fNPlanes==4) theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
//
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
@@ -1287,7 +1282,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
Double_t FPTimeSum=0.0;
Int_t nFPTimeSum=0;
- for (Int_t ip = 0; ip < temp_planes; ip++ ){
+ for (Int_t ip = 0; ip < fNPlanes; ip++ ){
if ( fNPlaneTime[ip] != 0 ){
fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
FPTimeSum += fSumPlaneTime[ip];
@@ -1332,7 +1327,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// *second, we move the scintillators. here we use scintillator cuts to see
// *if a track should have been found.
- for(Int_t ip = 0; ip < temp_planes; ip++ ) {
+ for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
if (!fPlanes[ip])
return -1;
@@ -1353,7 +1348,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// *Wwe count the number of three adjacent scintillators too. (A signle track
// *shouldn't fire three adjacent scintillators.
- for(Int_t ip = 0; ip < temp_planes; ip++ ) {
+ for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
// Planes ip = 0 = 1X
// Planes ip = 2 = 2X
if (!fPlanes[ip]) return -1;
@@ -1398,7 +1393,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
// *look for clusters in y planes... (10 scins) !this assume both y planes have same
// *number of scintillators.
- for (Int_t ip = 1; ip < temp_planes; ip +=2 ) {
+ for (Int_t ip = 1; ip < fNPlanes; ip +=2 ) {
// Planes ip = 1 = 1Y
// Planes ip = 3 = 2Y
if (!fPlanes[ip]) return -1;
@@ -1539,13 +1534,42 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
{
Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+ Double_t hitPos;
+ Double_t hitDistance;
for (Int_t itrk=0; itrk<Ntracks; itrk++) {
THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
if (theTrack->GetIndex()==0) {
fBeta=theTrack->GetBeta();
+ for (Int_t ip = 0; ip < fNPlanes; ip++ ){
+ Double_t pl_xypos=0;
+ Double_t pl_zpos=0;
+ Int_t num_good_pad=0;
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+ for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
+ if (fGoodFlags[itrk][ip][iphit].goodScinTime) {
+ Int_t padind = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
+ pl_xypos+=fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
+ pl_zpos+=fPlanes[ip]->GetZpos()+ (padind%2)*fPlanes[ip]->GetDzpos();
+ num_good_pad++;
+ }
+ // cout << ip << " " << iphit << " " << fGoodFlags[itrk][ip][iphit].goodScinTime << " " << fGoodFlags[itrk][ip][iphit].goodTdcPos << " " << fGoodFlags[itrk][ip][iphit].goodTdcNeg << endl;
+ }
+ hitDistance=kBig;
+ if (num_good_pad !=0 ) {
+ pl_xypos=pl_xypos/num_good_pad;
+ pl_zpos=pl_zpos/num_good_pad;
+ hitPos = theTrack->GetY() + theTrack->GetPhi()*pl_zpos ;
+ if (ip%2 == 0) hitPos = theTrack->GetX() + theTrack->GetTheta()*pl_zpos ;
+ hitDistance = hitPos - pl_xypos;
+ fPlanes[ip]->SetTrackXPosition(theTrack->GetX() + theTrack->GetTheta()*pl_zpos );
+ fPlanes[ip]->SetTrackYPosition(theTrack->GetY() + theTrack->GetPhi()*pl_zpos );
+ }
+ // cout << " ip " << ip << " " << hitPos << " " << pl_xypos << " " << pl_zpos << " " << hitDistance << endl;
+ fPlanes[ip]->SetHitDistance(hitDistance);
}
- } //over tracks
-
+ }
+ } //over tracks
+ //
return 0;
}
//_____________________________________________________________________________
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index 7f621d2..61f3dea 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -425,6 +425,9 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode )
{"GoodNegAdcPulseAmp", "List of Negative ADC peak amp (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegAdcPulseAmp"},
{"GoodPosAdcPulseTime", "List of positive ADC time (passed TDC && ADC Min and Max cuts for either end)", "fGoodPosAdcPulseTime"},
{"GoodNegAdcPulseTime", "List of Negative ADC time (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegAdcPulseTime"},
+ {"DiffDisTrack", "Difference between track and scintillator position (cm)", "fHitDistance"},
+ {"TrackXPos", "Track X position at plane (cm)", "fTrackXPosition"},
+ {"TrackYPos", "Track Y position at plane (cm)", "fTrackYPosition"},
//{"ngoodhits", "Number of paddle hits (passed tof tolerance and used to determine the focal plane time )", "GetNGoodHits() "},
{ 0 }
};
@@ -514,6 +517,9 @@ void THcScintillatorPlane::Clear( Option_t* )
fpTime = -1.e4;
+ fHitDistance = kBig;
+ fTrackXPosition = kBig;
+ fTrackYPosition = kBig;
}
//_____________________________________________________________________________
@@ -861,7 +867,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fGoodPosAdcPed.at(padnum-1) = ((THcSignalHit*) frPosAdcPed->ConstructedAt(good_ielem_posadc))->GetData();
fGoodPosAdcPulseInt.at(padnum-1) = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(good_ielem_posadc))->GetData();
fGoodPosAdcPulseAmp.at(padnum-1) = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(good_ielem_posadc))->GetData();
- fGoodPosAdcPulseTime.at(padnum-1) = ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(good_ielem_posadc))->GetData();
+ fGoodPosAdcPulseTime.at(padnum-1) = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(good_ielem_posadc))->GetData();
}
@@ -878,7 +884,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fGoodNegAdcPed.at(padnum-1) = ((THcSignalHit*) frNegAdcPed->ConstructedAt(good_ielem_negadc))->GetData();
fGoodNegAdcPulseInt.at(padnum-1) = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(good_ielem_negadc))->GetData();
fGoodNegAdcPulseAmp.at(padnum-1) = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(good_ielem_negadc))->GetData();
- fGoodNegAdcPulseTime.at(padnum-1) = ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(good_ielem_negadc))->GetData();
+ fGoodNegAdcPulseTime.at(padnum-1) = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(good_ielem_negadc))->GetData();
}
//DEFINE THE "GOOD +TDC Multiplicities and Occupancies"
@@ -912,9 +918,9 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adc_peak);
adc_peak=hit->GetRawAdcHitNeg().GetPulseAmp();
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adc_peak);
- Double_t time_peak=hit->GetRawAdcHitPos().GetPulseTimeRaw();
+ Double_t time_peak=hit->GetRawAdcHitPos().GetPulseTime();
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCtime(time_peak);
- time_peak=hit->GetRawAdcHitNeg().GetPulseTimeRaw();
+ time_peak=hit->GetRawAdcHitNeg().GetPulseTime();
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCtime(time_peak);
//Define GoodTdcUnCorrTime
@@ -935,10 +941,10 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(fTofUsingInvAdc) {
timec_pos = tdc_pos*fScinTdcToTime
- fHodoPosInvAdcOffset[index]
- - fHodoPosInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0,adc_pos));
+ - fHodoPosInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
timec_neg = tdc_neg*fScinTdcToTime
- fHodoNegInvAdcOffset[index]
- - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0,adc_neg));
+ - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
} else { // Old style
timec_pos = tdc_pos*fScinTdcToTime - fHodoPosPhcCoeff[index]*
TMath::Sqrt(TMath::Max(0.0,adc_pos/fHodoPosMinPh[index]-1.0))
@@ -1005,7 +1011,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(fTofUsingInvAdc) {
timec_pos = tdc_pos*fScinTdcToTime
- fHodoPosInvAdcOffset[index]
- - fHodoPosInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0,adc_pos));
+ - fHodoPosInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
} else { // Old style
timec_pos = tdc_pos*fScinTdcToTime - fHodoPosPhcCoeff[index]*
TMath::Sqrt(TMath::Max(0.0,adc_pos/fHodoPosMinPh[index]-1.0))
@@ -1016,7 +1022,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(fTofUsingInvAdc) {
timec_neg = tdc_neg*fScinTdcToTime
- fHodoNegInvAdcOffset[index]
- - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0,adc_neg));
+ - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
} else { // Old style
timec_neg = tdc_neg*fScinTdcToTime - fHodoNegPhcCoeff[index]*
TMath::Sqrt(TMath::Max(0.0,adc_neg/fHodoNegMinPh[index]-1.0))
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index bc900ed..866c9f0 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -44,6 +44,9 @@ class THcScintillatorPlane : public THaSubDetector {
Int_t GetNelem() {return fNelem;}; // return number of paddles in this plane
Int_t GetNScinHits() {return fNScinHits;}; // Get # hits in plane (that pass min/max TDC cuts)
Int_t GetNGoodHits() {return fNGoodHits;}; // Get # hits in plane (used in determining focal plane time)
+ Double_t GetHitDistance() {return fHitDistance;}; // Distance between track and hit paddle
+ Double_t GetTrackXPosition() {return fTrackXPosition;}; // Distance track X position at plane
+ Double_t GetTrackYPosition() {return fTrackYPosition;}; // Distance track Y position at plane
Double_t GetSpacing() {return fSpacing;}; // spacing of paddles
Double_t GetSize() {return fSize;}; // paddle size
Double_t GetHodoSlop() {return fHodoSlop;}; // hodo slop
@@ -59,6 +62,9 @@ class THcScintillatorPlane : public THaSubDetector {
void SetFpTime(Double_t f) {fFptime=f;};
void SetNGoodHits(Int_t ng) {fNGoodHits=ng;};
+ void SetHitDistance(Double_t f) {fHitDistance=f;}; // Distance between track and hit paddle
+ void SetTrackXPosition(Double_t f) {fTrackXPosition=f;}; // Distance track X position at plane
+ void SetTrackYPosition(Double_t f) {fTrackYPosition=f;}; // Distance track Y position at plane
TClonesArray* fParentHitList;
@@ -152,6 +158,9 @@ class THcScintillatorPlane : public THaSubDetector {
Int_t fDebugAdc;
+ Double_t fHitDistance;
+ Double_t fTrackXPosition;
+ Double_t fTrackYPosition;
Int_t fCosmicFlag; //
Int_t fPlaneNum; /* Which plane am I 1-4 */
UInt_t fTotPlanes; /* so we can read variables that are not indexed by plane id */
@@ -204,7 +213,6 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t *fHodoPosInvAdcAdc;
Double_t *fHodoNegInvAdcAdc;
Double_t *fHodoSigma;
-
Double_t fTolerance; /* need this for Focal Plane Time estimation */
Double_t fFptime;
/* Pedestal Quantities */
--
GitLab