Newer
Older
std::vector<Double_t> dedx_temp;
fdEdX.push_back(dedx_temp); // Create array of dedx per hit
std::vector<std::vector<GoodFlags> > goodflagstmp1;
fGoodFlags.push_back(goodflagstmp1);
// timeAtFP[itrack] = 0.;
Double_t sumFPTime = 0.; // Line 138
fNScinHit.push_back(0);
//! 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.
//! 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 hasnt defined something reasonable
// Loop over scintillator planes.
// In ENGINE, its loop over good scintillator hits.
fTOFCalc.clear(); // SAW - Can we
fTOFPInfo.clear(); // SAW - combine these two?
Int_t ihhit = 0; // Hit # overall
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
std::vector<GoodFlags> goodflagstmp2;
fGoodFlags[itrack].push_back(goodflagstmp2);
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Double_t zPos = fPlanes[ip]->GetZpos();
Double_t dzPos = fPlanes[ip]->GetDzpos();
// first loop over hits with in a single plane
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
// iphit is hit # within a plane
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
fTOFPInfo.push_back(TOFPInfo());
// Can remove these as we will initialize in the constructor
// fTOFPInfo[ihhit].time_pos = -99.0;
// fTOFPInfo[ihhit].time_neg = -99.0;
// fTOFPInfo[ihhit].keep_pos = kFALSE;
// fTOFPInfo[ihhit].keep_neg = kFALSE;
fTOFPInfo[ihhit].scin_pos_time = 0.0;
fTOFPInfo[ihhit].scin_neg_time = 0.0;
fTOFPInfo[ihhit].hit = hit;
fTOFPInfo[ihhit].planeIndex = ip;
fTOFPInfo[ihhit].hitNumInPlane = iphit;
fTOFPInfo[ihhit].onTrack = kFALSE;
Int_t paddle = hit->GetPaddleNumber()-1;
Double_t zposition = zPos + (paddle%2)*dzPos;
Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
( zposition ); // Line 183
Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
( zposition ); // Line 184
if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 185
scinTrnsCoord = xHitCoord;
scinLongCoord = yHitCoord;
} else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 188
scinTrnsCoord = yHitCoord;
scinLongCoord = xHitCoord;
fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;
Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// Index to access the 2d arrays of paddle/scintillator properties
Int_t fPIndex = GetScinIndex(ip,paddle);
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
fTOFPInfo[ihhit].onTrack = kTRUE;
Double_t zcor = zposition/(29.979*fBetaNominal)*
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());
Double_t tdc_pos = hit->GetPosTDC();
if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) {
Double_t adc_pos = hit->GetPosADC();
Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
fTOFPInfo[ihhit].pathp = pathp;
Double_t timep = tdc_pos*fScinTdcToTime;
if(fTofUsingInvAdc) {
timep -= fHodoPosInvAdcOffset[fPIndex]
+ pathp/fHodoPosInvAdcLinear[fPIndex]
+ fHodoPosInvAdcAdc[fPIndex]
/TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
//Double_t tw_corr_pos = fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
Double_t tw_corr_pos=0.;
pathp=scinLongCoord;
if (adcamp_pos>0) tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
timep += -tw_corr_pos + fHodo_LCoeff[fPIndex]+ pathp/fHodoVelFit[fPIndex];
fTOFPInfo[ihhit].scin_pos_time = timep;
timep -= zcor;
fTOFPInfo[ihhit].time_pos = timep;
Double_t tdc_neg = hit->GetNegTDC();
if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) {
Double_t adc_neg = hit->GetNegADC();
Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
fTOFPInfo[ihhit].pathn = pathn;
Double_t timen = tdc_neg*fScinTdcToTime;
if(fTofUsingInvAdc) {
timen -= fHodoNegInvAdcOffset[fPIndex]
+ pathn/fHodoNegInvAdcLinear[fPIndex]
+ fHodoNegInvAdcAdc[fPIndex]
/TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
pathn=scinLongCoord ;
Double_t tw_corr_neg =0 ;
if (adcamp_neg >0) tw_corr_neg= 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]);
timen += -tw_corr_neg- 2*fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex]- pathn/fHodoVelFit[fPIndex];
fTOFPInfo[ihhit].scin_neg_time = timen;
fTOFPInfo[ihhit].time_neg = timen;
} // condition for cenetr on a paddle
} // First loop over hits in a plane <---------
//-----------------------------------------------------------------------------------------------
//------------- First large loop over scintillator hits ends here --------------------
//-----------------------------------------------------------------------------------------------
Int_t nhits=ihhit;
if(0.5*hTime->GetMaximumBin() > 0) {
Double_t tmin = 0.5*hTime->GetMaximumBin();
for(Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
if ( ( fTOFPInfo[ih].time_pos > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_pos < ( tmin + fTofTolerance ) ) ) {
fTOFPInfo[ih].keep_pos=kTRUE;
if ( ( fTOFPInfo[ih].time_neg > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_neg < ( tmin + fTofTolerance ) ) ){
fTOFPInfo[ih].keep_neg=kTRUE;
//---------------------------------------------------------------------------------------------
// ---------------------- Scond loop over scint. hits in a plane ------------------------------
//---------------------------------------------------------------------------------------------
for(Int_t ih=0; ih < nhits; ih++) {
THcHodoHit *hit = fTOFPInfo[ih].hit;
Int_t iphit = fTOFPInfo[ih].hitNumInPlane;
Int_t ip = fTOFPInfo[ih].planeIndex;
// fDumpOut << " looping over hits = " << ih << " plane = " << ip+1 << endl;
GoodFlags flags;
// Flags are used by THcHodoEff
fGoodFlags[itrack][ip].push_back(flags);
fGoodFlags[itrack][ip][iphit].onTrack = kFALSE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kFALSE;
fGoodFlags[itrack][ip][iphit].goodTdcNeg = kFALSE;
fGoodFlags[itrack][ip][iphit].goodTdcPos = kFALSE;
fTOFCalc.push_back(TOFCalc());
// Do we set back to false for each track, or just once per event?
fTOFCalc[ih].good_scin_time = kFALSE;
// These need a track index too to calculate efficiencies
fTOFCalc[ih].good_tdc_pos = kFALSE;
fTOFCalc[ih].good_tdc_neg = kFALSE;
fTOFCalc[ih].pindex = ip;
Int_t paddle = hit->GetPaddleNumber()-1;
fTOFCalc[ih].hit_paddle = paddle;
fTOFCalc[ih].good_raw_pad = paddle;
// Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// Double_t scinTrnsCoord = fTOFPInfo[ih].scinTrnsCoord;
// Double_t scinLongCoord = fTOFPInfo[ih].scinLongCoord;
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;
if ( fTOFPInfo[ih].keep_neg ) { //
fTOFCalc[ih].good_tdc_neg = kTRUE;
fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
// ** Calculate ave time for scin and error.
if ( fTOFCalc[ih].good_tdc_pos ){
if ( fTOFCalc[ih].good_tdc_neg ){
fTOFCalc[ih].scin_time = ( fTOFPInfo[ih].scin_pos_time +
fTOFCalc[ih].scin_time_fp = ( fTOFPInfo[ih].time_pos +
fTOFPInfo[ih].time_neg ) / 2.;
if (fTofUsingInvAdc){
fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.;
}
else {
fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoSigmaPos[fPIndex] * fHodoSigmaPos[fPIndex] +
fHodoSigmaNeg[fPIndex] * fHodoSigmaNeg[fPIndex] )/2.;
}
fTOFCalc[ih].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
} else{
fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_pos_time;
fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_pos;
if (fTofUsingInvAdc){
fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex];
}
else{
fTOFCalc[ih].scin_sigma = fHodoSigmaPos[fPIndex];
}
fTOFCalc[ih].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
if ( fTOFCalc[ih].good_tdc_neg ){
fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_neg_time;
fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_neg;
if (fTofUsingInvAdc){
fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex];
}
else{
fTOFCalc[ih].scin_sigma = fHodoSigmaNeg[fPIndex];
}
fTOFCalc[ih].good_scin_time = kTRUE;
fGoodFlags[itrack][ip][iphit].goodScinTime = 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 ( fTOFCalc[ih].good_scin_time ){
// scin_time_fp doesn't need to be an array
// Is this any different than the average of time_pos and time_neg?
// Double_t scin_time_fp = ( fTOFPInfo[ih].time_pos +
// fTOFPInfo[ih].time_neg ) / 2.;
Double_t scin_time_fp = fTOFCalc[ih].scin_time_fp;
sumFPTime = sumFPTime + scin_time_fp;
nFPTime ++;
fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
fNPlaneTime[ip] ++;
fNScinHit[itrack] ++;
if ( ( fTOFCalc[ih].good_tdc_pos ) && ( fTOFCalc[ih].good_tdc_neg ) ){
nPmtHit[itrack] = nPmtHit[itrack] + 2;
} else {
nPmtHit[itrack] = nPmtHit[itrack] + 1;
}
fdEdX[itrack].push_back(0.0);
// --------------------------------------------------------------------------------------------
if ( fTOFCalc[ih].good_tdc_pos ){
if ( fTOFCalc[ih].good_tdc_neg ){
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Sqrt( TMath::Max( 0., hit->GetPosADC() * hit->GetNegADC() ) );
} else{
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Max( 0., hit->GetPosADC() );
if ( fTOFCalc[ih].good_tdc_neg ){
fdEdX[itrack][fNScinHit[itrack]-1]=
TMath::Max( 0., hit->GetNegADC() );
} else{
fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
}
// --------------------------------------------------------------------------------------------
} // time at focal plane condition
} // on track condition
// ** See if there are any good time measurements in the plane.
if ( fTOFCalc[ih].good_scin_time ){
fGoodPlaneTime[ip] = kTRUE;
fTOFCalc[ih].dedx = fdEdX[itrack][fNScinHit[itrack]-1];
fTOFCalc[ih].dedx = 0.0;
} // Second loop over hits of a scintillator plane ends here
theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
if (fNumPlanesBetaCalc==4) theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
//
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// * * Fit beta if there are enough time measurements (one upper, one lower)
// From h_tof_fit
if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) &&
( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){
Double_t sumW = 0.;
Double_t sumT = 0.;
Double_t sumZ = 0.;
Double_t sumZZ = 0.;
Double_t sumTZ = 0.;
for(Int_t ih=0; ih < nhits; ih++) {
Int_t ip = fTOFPInfo[ih].planeIndex;
if ( fTOFCalc[ih].good_scin_time ) {
Double_t scinWeight = 1 / ( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );
Double_t zPosition = ( fPlanes[ip]->GetZpos()
+( fTOFCalc[ih].hit_paddle % 2 ) *
fPlanes[ip]->GetDzpos() );
sumW += scinWeight;
sumT += scinWeight * fTOFCalc[ih].scin_time;
sumZ += scinWeight * zPosition;
sumZZ += scinWeight * ( zPosition * zPosition );
sumTZ += scinWeight * zPosition * fTOFCalc[ih].scin_time;
} // condition of good scin time
} // loop over hits
Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
betaChiSq = 0.;
for(Int_t ih=0; ih < nhits; ih++) {
Int_t ip = fTOFPInfo[ih].planeIndex;
if ( fTOFCalc[ih].good_scin_time ){
Double_t zPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ih].hit_paddle % 2 ) *
fPlanes[ip]->GetDzpos() );
Double_t timeDif = ( fTOFCalc[ih].scin_time - t0 );
betaChiSq += ( ( zPosition / beta - timeDif ) *
( zPosition / beta - timeDif ) ) /
( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );
} // condition for good scin time
} // loop over hits
Double_t pathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi() );
// Take angle into account
beta = beta / 29.979; // velocity / c
} // else condition for fTmpDenom
timeAtFP[itrack] = ( sumFPTime / nFPTime );
//
// ---------------------------------------------------------------------------
Double_t FPTimeSum=0.0;
Int_t nFPTimeSum=0;
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
if ( fNPlaneTime[ip] != 0 ){
fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
FPTimeSum += fSumPlaneTime[ip];
nFPTimeSum += fNPlaneTime[ip];
fFPTime[ip] = 1000. * ( ip + 1 );
Double_t fptime=-1000;
if (nFPTimeSum>0) fptime = FPTimeSum/nFPTimeSum;
fFPTimeAll = fptime;
Double_t dedx=0.0;
for(UInt_t ih=0;ih<fTOFCalc.size();ih++) {
if(fTOFCalc[ih].good_scin_time) {
dedx = fTOFCalc[ih].dedx;
break;
}
}
theTrack->SetDedx(dedx);
theTrack->SetFPTime(fptime);
theTrack->SetBeta(beta);
theTrack->SetBetaChi2( betaChiSq );
theTrack->SetNPMT(nPmtHit[itrack]);
theTrack->SetFPTime( timeAtFP[itrack]);
} // Main loop over tracks ends here.
} // If condition for at least one track
return 0;
}
void THcHodoscope::TrackEffTest(void)
{
// assume X planes are 0,2 and Y planes are 1,3
std::array<int,4> PadLow = {fxLoScin[0], fyLoScin[0], fxLoScin[1], fyLoScin[1]};
std::array<int,4> PadHigh = {fxHiScin[0], fyHiScin[0], fxHiScin[1], fyHiScin[1]};
std::array<double,4> PadPosLo;
std::array<double,4> PadPosHi;
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
Double_t lowtemp = fPlanes[ip]->GetPosCenter(PadLow[ip] - 1) + fPlanes[ip]->GetPosOffset();
Double_t hitemp = fPlanes[ip]->GetPosCenter(PadHigh[ip] - 1) + fPlanes[ip]->GetPosOffset();
if (lowtemp < hitemp) {
PadPosLo[ip]=lowtemp;
PadPosHi[ip]=hitemp;
} else {
PadPosLo[ip]=hitemp;
PadPosHi[ip]=lowtemp;
}
}
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
//{
// std::vector<int> test_vector = {1,2,5,6,7, 9,10, 20};
// auto test_res = hcana::find_discontinuity(test_vector);
// std::cout << " find_discontinuity test: " << *test_res << "\n";
// std::cout << " count_discontinuities test: " << hcana::count_discontinuities(test_vector)
// << "\n";
// auto cont_ranges = hcana::get_discontinuities(test_vector);
// int ir = 0;
// std::cout << "ranges: \n";
// for(auto r : cont_ranges) {
// for(auto v = r.first; v< r.second; v++) {
// std::cout << *v << " ";
// }
// std::cout << "\n";
// }
//}
//{
// std::vector<int> test_vector = {0,2,4,8,10,12,13, 16,18,20,23};
// auto test_res = hcana::find_discontinuity(
// test_vector, [](const int& x1, const int& x2) { return x1 + 2 == x2; });
// std::cout << " find_discontinuity test2: " << *test_res << "\n";
// std::cout << " count_discontinuities test2: "
// << hcana::count_discontinuities(
// test_vector, [](const int& x1, const int& x2) { return x1 + 2 == x2; })
// << "\n";
//}
using HitIterator = std::vector<THcHodoHit*>::iterator;
using HitRangeVector = typename std::vector<std::pair<HitIterator, HitIterator>>;
using ClusterPositions = typename std::vector<double>;
std::vector<std::vector<THcHodoHit*>> hit_vecs;
std::vector<HitRangeVector> hit_clusters;
std::vector<ClusterPositions> pos_clusters;
std::vector<std::vector<int>> good_clusters;
std::array<int,4> n_good_clusters;
// Loop over each plane
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
// Get the hit array
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// Create vector for good hits.
hit_vecs.push_back(std::vector<THcHodoHit*>());
for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
// keep only those with "two good times"
hit_vecs[ip].push_back(hit);
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
// Cluster (group adjacent) hits
if( hit_vecs[ip].size() > 0 ) {
auto& hv = hit_vecs[ip];
hcana::sort(hv, [](const THcHodoHit* h1, const THcHodoHit* h2) {
return h1->GetPaddleNumber() < h2->GetPaddleNumber();
});
// get the clusters for this layer
hit_clusters.push_back(
hcana::get_discontinuities(hv, [](const THcHodoHit* h1, const THcHodoHit* h2) {
return h1->GetPaddleNumber() + 1 == h2->GetPaddleNumber();
}));
} else {
hit_clusters.push_back( HitRangeVector{} );
}
// Calculate each cluster's mean position (in one coordinate)
ClusterPositions clust_positions;
std::vector<int> clust_bound;
for(auto& r : hit_clusters.back()) {
double a_pos = 0.0;
int n_hit = 0;
for(auto chit = r.first; chit < r.second; chit++) {
a_pos += fPlanes[ip]->GetPosCenter((*chit)->GetPaddleNumber()-1) + fPlanes[ip]->GetPosOffset();
n_hit++;
}
// take the average position
double avg_pos = a_pos / double(n_hit);
// Calculate if it is the bound by the upper and lower limits where we expect
// full tracks to be reconstructed.
int is_bound = int((avg_pos >= PadPosLo[ip]) && (avg_pos <= PadPosHi[ip]));
clust_bound.push_back(is_bound);
if(is_bound) {
// we are only interested in clusters positioned inthe "sweet spot"
clust_positions.push_back(avg_pos);
n_good_clusters[ip]++;
if (n_hit > 10) {
_logger->warn("cluster in hodoscope track efficiency has {} hits", n_hit);
}
}
pos_clusters.push_back(clust_positions);
good_clusters.push_back(clust_bound);
//_logger->info("{} clusters in Hod plane {}", hit_clusters.back().size(), ip );
bool has_both_x_clusters = (n_good_clusters[0] > 0) && (n_good_clusters[2] > 0);
bool has_both_y_clusters = (n_good_clusters[1] > 0) && (n_good_clusters[3] > 0);
bool at_least_one_good_x_clusters = (n_good_clusters[0] > 0) || (n_good_clusters[2] > 0);
bool at_least_one_good_y_clusters = (n_good_clusters[1] > 0) || (n_good_clusters[3] > 0);
bool good_x_match = false;
bool good_y_match = false;
// Superficial matching. Find out if clusters in front and back could possibly belong to the same
// track. This done by checking to see if the position differences are less than 10 cm;
// ie, |x1-x2|<10 or |y1-y2| <10
// TODO do actual matching with tracks elsewhere
if(has_both_x_clusters) {
for(auto x1_pos : pos_clusters[0] ) {
for(auto x2_pos : pos_clusters[2] ) {
if( std::abs(x1_pos-x2_pos) < 10.0) {
good_x_match = true;
break;
}
}
if(has_both_y_clusters) {
for(auto y1_pos : pos_clusters[1] ) {
for(auto y2_pos : pos_clusters[3] ) {
if( std::abs(y1_pos-y2_pos) < 10.0) {
good_y_match = true;
break;
}
fGoodScinHits = 0;
// (2+2) 4 good hits
if( good_y_match && good_x_match ) {
fGoodScinHits = 1;
}
// (2+1) 3 good hits
if( at_least_one_good_y_clusters && good_x_match ) {
fGoodScinHits = 1;
}
// (1+2) 3 good hits
if( at_least_one_good_x_clusters && good_y_match ) {
fGoodScinHits = 1;
}
}
//
//
//
void THcHodoscope::OriginalTrackEffTest(void)
{
/**
Translation of h_track_tests.f file for tracking efficiency
*/
//************************now look at some hodoscope tests
// *second, we move the scintillators. here we use scintillator cuts to see
// *if a track should have been found.
cout << " enter track eff" << fNumPlanesBetaCalc << endl;
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
cout << " loop over planes " << ip+1 << endl;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
// TClonesArray* scinPosTDC = fPlanes[ip]->GetPosTDC();
// TClonesArray* scinNegTDC = fPlanes[ip]->GetNegTDC();
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
cout << " hits = " << fNScinHits[ip] << endl;
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
fScinHitPaddle[ip][paddle] = 1;
cout << " hit = " << iphit+1 << " " << paddle+1 << endl;
}
}
// *next, look for clusters of hits in a scin plane. a cluster is a group of
// *adjacent scintillator hits separated by a non-firing scintillator.
// *Wwe count the number of three adjacent scintillators too. (A signle track
// *shouldn't fire three adjacent scintillators.
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
// Planes ip = 0 = 1X
// Planes ip = 2 = 2X
fNClust.push_back(0);
fThreeScin.push_back(0);
}
// *look for clusters in x planes... (16 scins) !this assume both x planes have same
// *number of scintillators.
cout << " looking for cluster in x planes" << endl;
Int_t icount;
for (Int_t ip = 0; ip < 3; ip +=2 ) {
icount = 0;
if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
// !look for number of clusters of 1 or more hits
if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
icount ++;
cout << " paddle = " << ipaddle+1 << " " << icount << endl;
cout << "Two cluster in plane = " << ip+1 << " " << icount << endl;
fNClust[ip] = icount;
icount = 0;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
// !look for three or more adjacent hits
if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
icount ++;
} // Second loop over paddles
cout << "Three clusters in plane = " << ip+1 << " " << icount << endl;
if ( icount > 0 )
fThreeScin[ip] = 1;
} // Loop over X plane
// *look for clusters in y planes... (10 scins) !this assume both y planes have same
cout << " looking for cluster in y planes" << endl;
for (Int_t ip = 1; ip < fNumPlanesBetaCalc; ip +=2 ) {
// Planes ip = 1 = 1Y
// Planes ip = 3 = 2Y
icount = 0;
if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
cout << "plane =" << ip << "check if paddle 1 hit " << icount << endl;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
// !look for number of clusters of 1 or more hits
if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
icount ++;
cout << " paddle = " << ipaddle+1 << " " << icount << endl;
cout << "Two cluster in plane = " << ip+1 << " " << icount << endl;
fNClust[ip] = icount;
icount = 0;
for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
// !look for three or more adjacent hits
if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
icount ++;
} // Second loop over Y paddles
cout << "Three clusters in plane = " << ip+1 << " " << icount << endl;
if ( icount > 0 )
fThreeScin[ip] = 1;
}// Loop over Y planes
// *next we mask out the edge scintillators, and look at triggers that happened
// *at the center of the acceptance. To change which scins are in the mask
// *change the values of h*loscin and h*hiscin in htracking.param
// fGoodScinHits = 0;
for (Int_t ifidx = fxLoScin[0]; ifidx < (Int_t) fxHiScin[0]; ifidx ++ ){
fGoodScinHitsX.push_back(0);
}
fHitSweet1X=0;
fHitSweet2X=0;
fHitSweet1Y=0;
fHitSweet2Y=0;
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
// *first x plane. first see if there are hits inside the scin region
for (Int_t ifidx = fxLoScin[0]-1; ifidx < fxHiScin[0]; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){
fHitSweet1X = 1;
fSweet1XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fxLoScin[0]-1; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
}
for (Int_t ifidx = fxHiScin[0]; ifidx < (Int_t) fNPaddle[0]; ifidx ++ ){
if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
}
// *second x plane. first see if there are hits inside the scin region
for (Int_t ifidx = fxLoScin[1]-1; ifidx < fxHiScin[1]; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){
fHitSweet2X = 1;
fSweet2XScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fxLoScin[1]-1; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
}
for (Int_t ifidx = fxHiScin[1]; ifidx < (Int_t) fNPaddle[2]; ifidx ++ ){
if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
}
// *first y plane. first see if there are hits inside the scin region
for (Int_t ifidx = fyLoScin[0]-1; ifidx < fyHiScin[0]; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){
fHitSweet1Y = 1;
fSweet1YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fyLoScin[0]-1; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
}
for (Int_t ifidx = fyHiScin[0]; ifidx < (Int_t) fNPaddle[1]; ifidx ++ ){
if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
}
// *second y plane. first see if there are hits inside the scin region
for (Int_t ifidx = fyLoScin[1]-1; ifidx < fyHiScin[1]; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){
fHitSweet2Y = 1;
fSweet2YScin = ifidx + 1;
}
}
// * next make sure nothing fired outside the good region
for (Int_t ifidx = 0; ifidx < fyLoScin[1]-1; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
}
for (Int_t ifidx = fyHiScin[1]; ifidx < (Int_t) fNPaddle[3]; ifidx ++ ){
if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
}
fTestSum = fHitSweet1X + fHitSweet2X + fHitSweet1Y + fHitSweet2Y;
// * now define a 3/4 or 4/4 trigger of only good scintillators the value
if ( fTestSum >= fTrackEffTestNScinPlanes ){
fGoodScinHits = 1;
for (Int_t ifidx = fxLoScin[0]; ifidx < fxHiScin[0]; ifidx ++ ){
if ( fSweet1XScin == ifidx )
fGoodScinHitsX[ifidx] = 1;
}
}
// * require front/back hodoscopes be close to each other
if ( ( fGoodScinHits == 1 ) && ( fTrackEffTestNScinPlanes == 4 ) ){
if ( TMath::Abs( fSweet1XScin - fSweet2XScin ) > 3 )
fGoodScinHits = 0;
if ( TMath::Abs( fSweet1YScin - fSweet2YScin ) > 2 )
fGoodScinHits = 0;
}
Gabriel Niculescu
committed
//_____________________________________________________________________________
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();
fFPTimeAll=theTrack->GetFPTime();
Double_t shower_track_enorm = theTrack->GetEnergy()/theTrack->GetP();
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; 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 ( fTOFCalc[ih].good_tdc_pos && fTOFCalc[ih].good_tdc_neg ) {
Bool_t sh_pid=(shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi);
Bool_t beta_pid=( fBeta > fTOFCalib_beta_lo && fBeta < fTOFCalib_beta_hi);
// cer_pid is true if there is no Cherenkov detector
Bool_t cer_pid=(fCherenkov?(fCherenkov->GetCerNPE()>fTOFCalib_cer_lo):kTRUE);
if(fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && sh_pid && beta_pid && cer_pid) {
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;
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;
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);
}
if(ih>0&&fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi ) {
fDumpOut << "0 " << endl;
}
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;
}
//_____________________________________________________________________________
Int_t THcHodoscope::End(THaRunBase* run)
{
MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
return 0;
}
ClassImp(THcHodoscope)
////////////////////////////////////////////////////////////////////////////////