Newer
Older
for( ihit = 0; ihit < fNScinHits[ip]; ihit++) { // Loop over sinc. hits. in plane
if ( ( fTimePos[ihit] > fTMin ) && ( fTimePos[ihit] < ( fTMin + fTofTolerance ) ) ) {
fKeepPos[ihit]=kTRUE;
if ( ( fTimeNeg[ihit] > fTMin ) && ( fTimeNeg[ihit] < ( fTMin + fTofTolerance ) ) ){
fKeepNeg[ihit] = kTRUE;
} // fJMax > 0 condition
// fScinPosTime = new Double_t [MAXHODHITS];
// fScinNegTime = new Double_t [MAXHODHITS];
for ( k = 0; k < MAXHODHITS; k++ ){
fScinPosTime[k] = 0;
fScinNegTime[k] = 0;
//---------------------------------------------------------------------------------------------
// ---------------------- Scond loop over scint. hits in a plane ------------------------------
//---------------------------------------------------------------------------------------------
fPaddle = ((THcSignalHit*)scinPosTDC->At(ihit))->GetPaddleNumber()-1;
fXcoord = theTrack->GetX() + theTrack->GetTheta() *
( fPlanes[ip]->GetZpos() + ( fPaddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 277
fYcoord = theTrack->GetY() + theTrack->GetPhi() *
( fPlanes[ip]->GetZpos() + ( fPaddle % 2 ) * fPlanes[ip]->GetDzpos() ); // Line 278
if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 278
fScinTrnsCoord = fXcoord;
fScinLongCoord = fYcoord;
}
else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 281
fScinTrnsCoord = fYcoord;
fScinLongCoord = fXcoord;
}
else { return -1; } // Line 288
fScinCenter = fPlanes[ip]->GetPosCenter(fPaddle) + fPlanes[ip]->GetPosOffset();
fIndex = fNPlanes * fPaddle + ip;
// ** Check if scin is on track
if ( TMath::Abs( fScinCenter - fScinTrnsCoord ) >
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
// scinOnTrack[itrack][ihit] = kFALSE;
// scinOnTrack[itrack][ihit] = kTRUE;
// * * Check for good TDC
if ( ( ((THcSignalHit*)scinPosTDC->At(ihit))->GetData() > fScinTdcMin ) &&
( ((THcSignalHit*)scinPosTDC->At(ihit))->GetData() < fScinTdcMax ) &&
( fKeepPos[ihit] ) ) { // 301
// ** Calculate time for each tube with a good tdc. 'pos' side first.
fGoodTDCPos[fGoodTimeIndex] = kTRUE;
// fNtof ++;
adcPh[ihit] = ((THcSignalHit*)scinPosADC->At(ihit))->GetData();
fPath[ihit] = fPlanes[ip]->GetPosLeft() - fScinLongCoord;
// * Convert TDC value to time, do pulse height correction, correction for
// * propogation of light thru scintillator, and offset.
fTime[ihit] = ((THcSignalHit*)scinPosTDC->At(ihit))->GetData() * fScinTdcToTime;
fTime[ihit] = fTime[ihit] - ( fHodoPosPhcCoeff[fIndex] * TMath::Sqrt( TMath::Max( 0. ,
( ( adcPh[ihit] / fHodoPosMinPh[fIndex] ) - 1 ) ) ) );
fTime[ihit] = fTime[ihit] - ( fPath[ihit] / fHodoVelLight[fIndex] );
fScinPosTime[ihit] = fTime[ihit] - fHodoPosTimeOffset[fIndex];
} // check for good pos TDC condition
// ** Repeat for pmts on 'negative' side
if ( ( ((THcSignalHit*)scinNegTDC->At(ihit))->GetData() > fScinTdcMin ) &&
( ((THcSignalHit*)scinNegTDC->At(ihit))->GetData() < fScinTdcMax ) &&
( fKeepNeg[ihit] ) ) { //
// ** Calculate time for each tube with a good tdc. 'pos' side first.
fGoodTDCNeg[fGoodTimeIndex] = kTRUE;
// fNtof ++;
adcPh[ihit] = ((THcSignalHit*)scinNegADC->At(ihit))->GetData();
fPath[ihit] = fScinLongCoord - fPlanes[ip]->GetPosRight(); // pos = left, neg = right
// * Convert TDC value to time, do pulse height correction, correction for
// * propogation of light thru scintillator, and offset.
fTime[ihit] = ((THcSignalHit*)scinNegTDC->At(ihit))->GetData() * fScinTdcToTime;
fTime[ihit] = fTime[ihit] - ( fHodoNegPhcCoeff[fIndex] *
TMath::Sqrt( TMath::Max( 0. , ( ( adcPh[ihit] / fHodoNegMinPh[fIndex] ) - 1 ) ) ) );
fTime[ihit] = fTime[ihit] - ( fPath[ihit] / fHodoVelLight[fIndex] );
fScinNegTime[ihit] = fTime[ihit] - fHodoNegTimeOffset[fIndex];
} // check for good neg TDC condition
// ** Calculate ave time for scin and error.
if ( fGoodTDCPos[fGoodTimeIndex] ){
if ( fGoodTDCNeg[fGoodTimeIndex] ){
fScinTime[fGoodTimeIndex] = ( fScinPosTime[ihit] + fScinNegTime[ihit] ) / 2.;
fScinSigma[fGoodTimeIndex] = TMath::Sqrt( fHodoPosSigma[fIndex] * fHodoPosSigma[fIndex] +
fHodoNegSigma[fIndex] * fHodoNegSigma[fIndex] )/2.;
fGoodScinTime[fGoodTimeIndex] = kTRUE;
// fScinGoodTime[fGoodTimeIndex] = kTRUE;
// fNtofPairs ++;
fScinTime[fGoodTimeIndex] = fScinPosTime[ihit];
fScinSigma[fGoodTimeIndex] = fHodoPosSigma[fIndex];
fGoodScinTime[fGoodTimeIndex] = kTRUE;
// fScinGoodTime[fGoodTimeIndex] = kTRUE;
if ( fGoodTDCNeg[fGoodTimeIndex] ){
fScinTime[fGoodTimeIndex] = fScinNegTime[ihit];
fScinSigma[fGoodTimeIndex] = fHodoNegSigma[fIndex];
fGoodScinTime[fGoodTimeIndex] = kTRUE;
// fScinGoodTime[fGoodTimeIndex] = kTRUE;
} // In h_tof.f this includes the following if condition for time at focal plane
// // because it is written in FORTRAN code
// c Get time at focal plane
if ( fGoodScinTime[fGoodTimeIndex] ){
fScinTimefp[ihit] = fScinTime[fGoodTimeIndex] -
( fPlanes[ip]->GetZpos() + ( fPaddle % 2 ) * fPlanes[ip]->GetDzpos() ) /
( 29.979 * fBetaP ) *
TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
// ---------------------------------------------------------------------------
// Date: July 8 2014
//
// Right now we do not need this code for beta and chisquare
//
//
// fSumfpTime = fSumfpTime + fScinTimefp[ihit];
// fNfpTime ++;
//
// ---------------------------------------------------------------------------
fSumPlaneTime[ip] = fSumPlaneTime[ip] + fScinTimefp[ihit];
fNPlaneTime[ip] ++;
fNScinHit[itrack] ++;
// scinHit[itrack][fNScinHit[itrack]] = ihit;
// scinfFPTime[itrack][fNScinHit[itrack]] = fScinTimefp[ihit];
// ---------------------------------------------------------------------------
// Date: July 8 2014
// This counts the pmt hits. Right now we don't need it so it is commentd off
//
// if ( ( fGoodTDCPos[fGoodTimeIndex] ) && ( fGoodTDCNeg[fGoodTimeIndex] ) ){
// fNPmtHit[itrack] = fNPmtHit[itrack] + 2;
// }
// else {
// fNPmtHit[itrack] = fNPmtHit[itrack] + 1;
// }
// ---------------------------------------------------------------------------
// fdEdX[itrack][ihit] = 5.0;
// --------------------------------------------------------------------------------------------
// Date: July 8 201 May be we need this, not sure.
//
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
if ( fGoodTDCPos[fGoodTimeIndex] ){
if ( fGoodTDCNeg[fGoodTimeIndex] ){
fdEdX[itrack][fNScinHit[itrack]-1] =
TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() *
((THcSignalHit*)scinNegADC->At(ihit))->GetData() ) );
}
else{
fdEdX[itrack][fNScinHit[itrack]-1] =
TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(ihit))->GetData() );
}
}
else{
if ( fGoodTDCNeg[fGoodTimeIndex] ){
fdEdX[itrack][fNScinHit[itrack]-1] =
TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(ihit))->GetData() );
}
else{
fdEdX[itrack][fNScinHit[itrack]-1] = 0.;
}
}
// --------------------------------------------------------------------------------------------
} // time at focal plane condition
} // on track else condition
// ** See if there are any good time measurements in the plane.
if ( fGoodScinTime[fGoodTimeIndex] ){
fGoodPlaneTime[ip] = kTRUE;
} // Second loop over hits of a scintillator plane ends here
} // Loop over scintillator planes ends here
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// * * Fit beta if there are enough time measurements (one upper, one lower)
if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) &&
( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){
Double_t sumw, sumt, sumz, sumzz, sumtz;
Double_t scinWeight, tmp, t0, tmpDenom, pathNorm, zPosition, timeDif;
sumw = 0.; sumt = 0.; sumz = 0.; sumzz = 0.; sumtz = 0.;
fGoodTimeIndex = -1;
for ( ip = 0; ip < fNPlanes; ip++ ){
if (!fPlanes[ip])
return -1;
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
fGoodTimeIndex ++;
if ( fGoodScinTime[fGoodTimeIndex] ) {
scinWeight = 1 / ( fScinSigma[fGoodTimeIndex] * fScinSigma[fGoodTimeIndex] );
zPosition = ( fPlanes[ip]->GetZpos() + ( fHitPaddle[fGoodTimeIndex] % 2 ) * fPlanes[ip]->GetDzpos() );
sumw += scinWeight;
sumt += scinWeight * fScinTime[fGoodTimeIndex];
sumz += scinWeight * zPosition;
sumzz += scinWeight * ( zPosition * zPosition );
sumtz += scinWeight * zPosition * fScinTime[fGoodTimeIndex];
} // condition of good scin time
} // loop over hits of plane
} // loop over planes
tmp = sumw * sumzz - sumz * sumz ;
t0 = ( sumt * sumzz - sumz * sumtz ) / tmp ;
tmpDenom = sumw * sumtz - sumz * sumt;
if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
fBeta[itrack] = tmp / tmpDenom;
fBetaChisq[itrack] = 0.;
fGoodTimeIndex = -1;
for ( ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes
if (!fPlanes[ip])
return -1;
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
for (ihit = 0; ihit < fNScinHits[ip]; ihit++ ){ // Loop over hits of a plane
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
fGoodTimeIndex ++;
if ( fGoodScinTime[fGoodTimeIndex] ){
zPosition = ( fPlanes[ip]->GetZpos() + ( fHitPaddle[fGoodTimeIndex] % 2 ) * fPlanes[ip]->GetDzpos() );
timeDif = ( fScinTime[fGoodTimeIndex] - t0 );
fBetaChisq[itrack] += ( ( zPosition / fBeta[itrack] - timeDif ) * ( zPosition / fBeta[itrack] - timeDif ) ) /
( fScinSigma[fGoodTimeIndex] * fScinSigma[fGoodTimeIndex] );
} // condition for good scin time
} // loop over hits of a plane
} // loop over planes
pathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() + theTrack->GetPhi() * theTrack->GetPhi() );
fBeta[itrack] = fBeta[itrack] / pathNorm;
fBeta[itrack] = fBeta[itrack] / 29.979; // velocity / c
// cout << "track = " << itrack + 1
// << " beta = " << fBeta[itrack] << endl;
} // condition for tmpDenom
else {
fBeta[itrack] = 0.;
fBetaChisq[itrack] = -2.;
} // else condition for tmpDenom
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
// --------------------------------------------------------------------
fBeta[itrack] = 0.;
fBetaChisq[itrack] = -1;
// ---------------------------------------------------------------------------
// Date: July 8 2014
//
// Right now we do not need this code for beta and chisquare
//
// if ( fNfpTime != 0 ){
// // fTimeAtFP[itrack] = ( fSumfpTime / fNfpTime );
// }
//
// ---------------------------------------------------------------------------
for ( Int_t ind = 0; ind < fNPlanes; ind++ ){
if ( fNPlaneTime[ind] != 0 ){
fFPTime[ind] = ( fSumPlaneTime[ind] / fNPlaneTime[ind] );
fFPTime[ind] = 1000. * ( ind + 1 );
theTrack->SetDedx(fdEdX[itrack][0]);
theTrack->SetBeta(fBeta[itrack]);
theTrack->SetBetaChi2( fBetaChisq[itrack] );
} // Main loop over tracks ends here.
} // If condition for at least one track
return 0;
Gabriel Niculescu
committed
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
// 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;
Gabriel Niculescu
committed
}
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nSide, Int_t nPlane, Int_t nPaddle ) {
return nSide*fMaxHodoScin+fNPlanes*nPaddle+nPlane-1;
}
//_____________________________________________________________________________
Double_t THcHodoscope::GetPathLengthCentral() {
return fPathLengthCentral;
}
ClassImp(THcHodoscope)
////////////////////////////////////////////////////////////////////////////////