diff --git a/hc_hodo_calib/tofcal.f b/hc_hodo_calib/tofcal.f
index 39c7917caefecc245e8c2d3ba9c82ae896d69416..5fe197d5b13077a475d28e153079eb8bab20ed41 100644
--- a/hc_hodo_calib/tofcal.f
+++ b/hc_hodo_calib/tofcal.f
@@ -50,7 +50,7 @@
open(unit=9,file=filename,err=998)
if(filename(1:1).ne.'h'.and.filename(1:1).ne.'s') then
write(6,'(1x,''error, input file name '',
- > '' should start with h or s'')')
+ > '' sould start with h or s'')')
goto 999
endif
write(*,*) ' filling adc histograms'
@@ -102,7 +102,7 @@ c write(13,'(1x,''adcmin='',f6.1)') adcmin(i)
! Do everything twice, 2nd time with tighter tolerance (ttol) on
! time differences, based on first pass
- ttol=5.0
+ ttol=20.0
do iloop=1,2
write(*,*) ' loop = ',iloop
! Initialize the fitting arrays
@@ -128,7 +128,6 @@ c write(13,'(1x,''adcmin='',f6.1)') adcmin(i)
n=n+1
read(string,'(1x,i1,2i3,5f10.3)') ipn,ipl,idt,
> tr(n),p(n),zc(n),tc1(n),adc(n)
-
! linearize the detector numbers
idet(n) = 100*(ipn-1) + 20*(ipl-1) + idt
if(idet(n).lt.1.or.idet(n).gt.200) write(6,'(1x,
@@ -149,7 +148,6 @@ c adc(n) = min(500, max(0., adc(n)))
! Loop over all pairs, if at least 6
11 if(n.ge.6) then
-
! see if this is first time a detector is used
do j=1,n
nhit(idet(j))=nhit(idet(j))+1
@@ -158,7 +156,7 @@ c adc(n) = min(500, max(0., adc(n)))
! Note that detector had has a fixed time offset (ip1) of zero
! since all times are relative)
if(nhit(idet(j)).eq.1) then
- if(idet(j).eq.4) then
+ if(idet(j).eq.10) then
ip1(idet(j))=0
else
@@ -334,17 +332,17 @@ c adc(n) = min(500, max(0., adc(n)))
enddo
write(10+iloop,'(/a,''hodo_pos_invadc_linear ='',3(f8.2,'',''),
- > f8.2)')filename(1:1),( -1./min(-0.02,vel(i)),i=1,80,20)
+ > f8.2)')filename(1:1),( -1./min(-1./15.,vel(i)),i=1,80,20)
do j=2,16
write(10+iloop,'(1x,'' '',3(f8.2,'',''),
- > f8.2)')(-1./min(-0.02,vel(i)),i=j,79+j,20)
+ > f8.2)')(-1./min(-1./15.,vel(i)),i=j,79+j,20)
enddo
write(10+iloop,'(/a,''hodo_neg_invadc_linear ='',3(f8.2,'',''),
- > f8.2)')filename(1:1),( -1./min(-0.02,vel(i)),i=101,180,20)
+ > f8.2)')filename(1:1),( -1./min(-1./15.,vel(i)),i=101,180,20)
do j=2,16
write(10+iloop,'(1x,'' '',3(f8.2,'',''),
- > f8.2)')(-1./min(-0.02,vel(i)),i=100+j,179+j,20)
+ > f8.2)')(-1./min(-1./15.,vel(i)),i=100+j,179+j,20)
enddo
write(10+iloop,'(/a,''hodo_pos_invadc_adc='',3(f9.2,'',''),
diff --git a/src/THcHitList.cxx b/src/THcHitList.cxx
index eafad160db1562747325c21badb75d12f4f91083..41f7815ad4b5136f63e89e28e6032afe5ae3f7ae 100644
--- a/src/THcHitList.cxx
+++ b/src/THcHitList.cxx
@@ -203,7 +203,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata ) {
Int_t reftime = evdata.GetData(d->crate, d->slot, d->refchan, 0);
rawhit->SetReference(signal, reftime);
} else {
- cout << "HitList: refchan " << d->refchan <<
+ cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
" missing for (" << d->crate << ", " << d->slot <<
", " << chan << ")" << endl;
}
@@ -212,7 +212,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata ) {
if(fRefIndexMaps[d->refindex].hashit) {
rawhit->SetReference(signal, fRefIndexMaps[d->refindex].reftime);
} else {
- cout << "HitList: refindex " << d->refindex <<
+ cout << "HitList(event=" << evdata.GetEvNum() << "): refindex " << d->refindex <<
" (" << fRefIndexMaps[d->refindex].crate <<
", " << fRefIndexMaps[d->refindex].slot <<
", " << fRefIndexMaps[d->refindex].channel << ")" <<
@@ -223,7 +223,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata ) {
}
} else { // This is a Flash ADC
- if (fPSE125 && fPSE125->IsPresent(d->crate)) { // Set F250 parameters.
+ if (fPSE125) { // Set F250 parameters.
rawhit->SetF250Params(
fPSE125->GetNSA(d->crate),
fPSE125->GetNSB(d->crate),
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index cb4f531acad08a5cd5dd24deea7f680a785c85fd..d82373e1734c389ac961dbdfe1e738e3026068b5 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -40,6 +40,7 @@ hodoscope array, not just one plane.
#include <cstdio>
#include <cstdlib>
#include <iostream>
+#include <iomanip>
#include <fstream>
using namespace std;
@@ -85,6 +86,10 @@ void THcHodoscope::Setup(const char* name, const char* description)
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
+ TString temp(prefix[0]);
+ fSHMS=kFALSE;
+ if (temp == "p" ) fSHMS=kTRUE;
+ cout << " fSHMS = " << fSHMS << endl;
string planenamelist;
DBRequest listextra[]={
{"hodo_num_planes", &fNPlanes, kInt},
@@ -360,9 +365,9 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date )
if (fCosmicFlag==1) cout << " setup for cosmics in TOF"<< endl;
cout << " cosmic flag = " << fCosmicFlag << endl;
if(fDumpTOF) {
- fDumpOut.open(fTOFDumpFile.c_str());
- if(fDumpOut.is_open()) {
- fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
+ fDumpOut.open(fTOFDumpFile.c_str());
+ if(fDumpOut.is_open()) {
+ //fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
} else {
fDumpTOF = 0;
cout << "WARNING: Unable to open TOF Dump file " << fTOFDumpFile << endl;
@@ -561,7 +566,7 @@ void THcHodoscope::ClearEvent()
fFPTime[ip]=0.;
fPlaneCenter[ip]=0.;
fPlaneSpacing[ip]=0.;
- for(Int_t iPaddle=0;iPaddle<fNPaddle[ip]; ++iPaddle) {
+ for(UInt_t iPaddle=0;iPaddle<fNPaddle[ip]; ++iPaddle) {
fScinHitPaddle[ip][iPaddle]=0;
}
}
@@ -706,7 +711,9 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
Bool_t goodplanetime[fNPlanes];
Bool_t twogoodtimes[nscinhits];
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t temp_planes=fNPlanes;
+ if (fSHMS) temp_planes=3;
+ for(Int_t ip=0;ip<temp_planes;ip++) {
goodplanetime[ip] = kFALSE;
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -722,10 +729,8 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
if ((postime>tmin) && (postime<tmin+fTofTolerance) &&
(negtime>tmin) && (negtime<tmin+fTofTolerance)) {
hit->SetTwoGoodTimes(kTRUE);
- twogoodtimes[ihit] = kTRUE;
- // Both tubes fired
- Int_t index=hit->GetPaddleNumber()-1;
- // Need to put this in a multihit histo
+ twogoodtimes[ihit] = kTRUE; // Both tubes fired
+ Int_t index=hit->GetPaddleNumber()-1; //
Double_t fptime;
if(fCosmicFlag==1) {
fptime = hit->GetScinCorrectedTime()
@@ -772,7 +777,9 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
Double_t sumTZ = 0.;
Int_t ihhit = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ temp_planes=fNPlanes;
+ if (fSHMS) temp_planes=3;
+ for(Int_t ip=0;ip<temp_planes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -785,8 +792,7 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
Double_t scinWeight = 1 / TMath::Power(sigma,2);
Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
-
- // cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " << sigma << endl;
+ // cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " << sigma << endl;
sumW += scinWeight;
sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
@@ -809,7 +815,7 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
fBetaNoTrkChiSq = 0.;
ihhit = 0;
- for (Int_t ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes
+ for (Int_t ip = 0; ip < temp_planes; ip++ ){ // Loop over planes
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
@@ -842,6 +848,9 @@ void THcHodoscope::EstimateFocalPlaneTime( void )
fBetaNoTrkChiSq = -2.;
} // else condition for fTmpDenom
}
+ fGoodEventTOFCalib=kFALSE;
+ if (!fSHMS&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&goodplanetime[3]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1&&fPlanes[3]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
+ if (fSHMS&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
}
//_____________________________________________________________________________
@@ -874,6 +883,8 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
// -------------------------------------------------
// fDumpOut << " ntrack = " << ntracks << endl;
+ Int_t temp_planes=fNPlanes;
+ if (fSHMS) temp_planes=3;
if (tracks.GetLast()+1 > 0 ) {
@@ -887,7 +898,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
if (!theTrack) return -1;
- for (Int_t ip = 0; ip < fNPlanes; ip++ ){
+ for (Int_t ip = 0; ip < temp_planes; ip++ ){
fGoodPlaneTime[ip] = kFALSE;
fNScinHits[ip] = 0;
fNPlaneTime[ip] = 0;
@@ -928,7 +939,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
fTOFPInfo.clear(); // SAW - combine these two?
Int_t ihhit = 0; // Hit # overall
- for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
+ for(Int_t ip = 0; ip < temp_planes; ip++ ) {
std::vector<GoodFlags> goodflagstmp2;
fGoodFlags[itrack].push_back(goodflagstmp2);
@@ -989,12 +1000,16 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
fTOFPInfo[ihhit].onTrack = kTRUE;
-
Double_t zcor = zposition/(29.979*fBetaP)*
TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
+ theTrack->GetPhi()*theTrack->GetPhi());
fTOFPInfo[ihhit].zcor = zcor;
-
+ if (fCosmicFlag) {
+ Double_t zcor = -zposition/(29.979*1.0)*
+ TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
+ + theTrack->GetPhi()*theTrack->GetPhi());
+ fTOFPInfo[ihhit].zcor = zcor;
+ }
Double_t tdc_pos = hit->GetPosTDC();
if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) {
Double_t adc_pos = hit->GetPosADC();
@@ -1013,7 +1028,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
+ fHodoPosTimeOffset[fPIndex];
}
fTOFPInfo[ihhit].scin_pos_time = timep;
- timep -= zcor;
+ timep -= zcor;
fTOFPInfo[ihhit].time_pos = timep;
for ( Int_t k = 0; k < 200; k++ ){ // Line 211
@@ -1122,21 +1137,17 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
if ( fTOFPInfo[ih].keep_pos ) { // 301
fTOFCalc[ih].good_tdc_pos = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
- if(fDumpTOF && ntracks==1) {
- fDumpOut << "1 " << ip+1 << " " << paddle+1 << " " <<
- hit->GetPosTDC()*fScinTdcToTime << " " << fTOFPInfo[ih].pathp <<
- " " << fTOFPInfo[ih].zcor << " " << fTOFPInfo[ih].time_pos <<
- " " << hit->GetPosADC() << endl;
+ if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
+ fDumpOut << fixed << setprecision(2);
+ fDumpOut << showpoint << " 1" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber() << setw(10) << hit->GetPosTDC()*fScinTdcToTime << setw(10) << fTOFPInfo[ih].pathp << setw(10) << fTOFPInfo[ih].zcor << setw(10) << fTOFPInfo[ih].time_pos << setw(10) << hit->GetPosADC() << endl;
}
}
if ( fTOFPInfo[ih].keep_neg ) { //
fTOFCalc[ih].good_tdc_neg = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
- if(fDumpTOF && ntracks==1) {
- fDumpOut << "2 " << ip+1 << " " << paddle+1 << " " <<
- hit->GetNegTDC()*fScinTdcToTime << " " << fTOFPInfo[ih].pathn <<
- " " << fTOFPInfo[ih].zcor << " " << fTOFPInfo[ih].time_neg <<
- " " << hit->GetNegADC() << endl;
+ if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
+ fDumpOut << fixed << setprecision(2);
+ fDumpOut << showpoint << " 2" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber() << setw(10) << hit->GetNegTDC()*fScinTdcToTime << setw(10) << fTOFPInfo[ih].pathn << setw(10) << fTOFPInfo[ih].zcor << setw(10) << fTOFPInfo[ih].time_neg << setw(10) << hit->GetNegADC() << endl;
}
}
@@ -1225,11 +1236,10 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
}
} // Second loop over hits of a scintillator plane ends here
-
- theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
+ theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
-
- //------------------------------------------------------------------------------
+ //
+ //------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
@@ -1317,7 +1327,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
Double_t FPTimeSum=0.0;
Int_t nFPTimeSum=0;
- for (Int_t ip = 0; ip < fNPlanes; ip++ ){
+ for (Int_t ip = 0; ip < temp_planes; ip++ ){
if ( fNPlaneTime[ip] != 0 ){
fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
FPTimeSum += fSumPlaneTime[ip];
@@ -1348,8 +1358,8 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks )
} // If condition for at least one track
- if(fDumpTOF && ntracks==1 ) {
- fDumpOut << "0 " << endl;
+ if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
+ fDumpOut << "0 " << endl;
}
//-----------------------------------------------------------------------
@@ -1362,7 +1372,7 @@ Int_t THcHodoscope::FineProcess( 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 < fNPlanes; ip++ ) {
+ for(Int_t ip = 0; ip < temp_planes; ip++ ) {
if (!fPlanes[ip])
return -1;
@@ -1383,7 +1393,7 @@ Int_t THcHodoscope::FineProcess( 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 < fNPlanes; ip++ ) {
+ for(Int_t ip = 0; ip < temp_planes; ip++ ) {
// Planes ip = 0 = 1X
// Planes ip = 2 = 2X
if (!fPlanes[ip]) return -1;
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index 27a9d6b3111c1c2cd3d3b80276fdc700ae029616..4e3cd4e0d5367c01591c5ab143768b2d1bd65b1b 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -126,6 +126,7 @@ protected:
// Calibration
// Per-event data
+ Bool_t fSHMS;
Bool_t fGoodStartTime;
Double_t fStartTime;
Double_t fFPTimeAll;
@@ -212,8 +213,9 @@ protected:
Int_t fNHodoscopes;
Int_t fDumpTOF;
- ofstream fDumpOut;
+ ofstream fDumpOut;
string fTOFDumpFile;
+ Bool_t fGoodEventTOFCalib;
Int_t fHitSweet1X;
diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx
index ffa690bc7c1cca2a98f8d657476528a5a55ee1fc..29fec53b5f5f992917bc29ca61860a3f4ec506f5 100644
--- a/src/THcScintillatorPlane.cxx
+++ b/src/THcScintillatorPlane.cxx
@@ -214,15 +214,19 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
{"hodo_adc_mode", &fADCMode, kInt, 0, 1},
{"hodo_pedestal_scale", &fADCPedScaleFactor, kDouble, 0, 1},
{"hodo_adc_diag_cut", &fADCDiagCut, kInt, 0, 1},
- {0}
+ {"cosmicflag", &fCosmicFlag, kInt, 0, 1},
+ {0}
};
fTofUsingInvAdc = 1;
fADCMode = kADCStandard;
fADCPedScaleFactor = 1.0;
fADCDiagCut = 50.0;
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- // fetch the parameter from the temporary list
+ fCosmicFlag=0;
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ if (fCosmicFlag==1) cout << " setup for cosmics in scint plane"<< endl;
+ cout << " cosmic flag = " << fCosmicFlag << endl;
+ // fetch the parameter from the temporary list
// Retrieve parameters we need from parent class
// Common for all planes
@@ -675,14 +679,24 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
timec_neg -= (hit_position-fPosRight)/
fHodoNegInvAdcLinear[index];
scin_corrected_time = 0.5*(timec_pos+timec_neg);
+ if (fCosmicFlag) {
+ postime = timec_pos + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ negtime = timec_neg + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ } else {
postime = timec_pos - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
negtime = timec_neg - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ }
} else {
postime=timec_pos-(fPosLeft-hit_position)/fHodoVelLight[index];
negtime=timec_neg-(hit_position-fPosRight)/fHodoVelLight[index];
scin_corrected_time = 0.5*(postime+negtime);
- postime = postime-(fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
- negtime = negtime-(fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ if (fCosmicFlag) {
+ postime = timec_pos + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ negtime = timec_neg + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ } else {
+ postime = timec_pos - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ negtime = timec_neg - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
+ }
}
// cout << fNScinHits<< " " << timec_pos << " " << timec_neg << endl;
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPaddleCenter(fPosCenter[index]);
diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h
index 123eaf5348ea3bf2451cc2efdd5765c17af89002..e14fa397c9da0ed1486483bbb9f8541c3738fc84 100644
--- a/src/THcScintillatorPlane.h
+++ b/src/THcScintillatorPlane.h
@@ -94,6 +94,7 @@ class THcScintillatorPlane : public THaSubDetector {
TClonesArray* frNegAdcPulseInt;
TClonesArray* frNegAdcPulseAmp;
+ 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 */
UInt_t fNelem; /* Need since we don't inherit from