Newer
Older
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));
} else {
timep -= fHodoPosPhcCoeff[fPIndex]*
TMath::Sqrt(TMath::Max(0.0,adc_pos/fHodoPosMinPh[fPIndex]-1.0))
+ pathp/fHodoVelLight[fPIndex]
+ fHodoPosTimeOffset[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));
} else {
timen -= fHodoNegPhcCoeff[fPIndex]*
TMath::Sqrt(TMath::Max(0.0,adc_neg/fHodoNegMinPh[fPIndex]-1.0))
+ pathn/fHodoVelLight[fPIndex]
+ fHodoNegTimeOffset[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 +
fTOFPInfo[ih].scin_neg_time ) / 2.;
fTOFCalc[ih].scin_time_fp = ( fTOFPInfo[ih].time_pos +
fTOFPInfo[ih].time_neg ) / 2.;
fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
fHodoNegSigma[fPIndex] * fHodoNegSigma[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;
fTOFCalc[ih].scin_sigma = fHodoPosSigma[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;
fTOFCalc[ih].scin_sigma = fHodoNegSigma[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
} // condition for fTmpDenom
} // 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 );
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
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
//OriginalTrackEffTest();
TrackEffTest();
return 0;
}
//
void THcHodoscope::TrackEffTest(void)
{
Double_t PadLow[4];
Double_t PadHigh[4];
// assume X planes are 0,2 and Y planes are 1,3
PadLow[0]=fxLoScin[0];
PadLow[2]=fxLoScin[1];
PadLow[1]=fyLoScin[0];
PadLow[3]=fyLoScin[1];
PadHigh[0]=fxHiScin[0];
PadHigh[2]=fxHiScin[1];
PadHigh[1]=fyHiScin[0];
PadHigh[3]=fyHiScin[1];
//
Double_t PadPosLo[4];
Double_t PadPosHi[4];
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;
}
}
//
const Int_t MaxNClus=5;
std::vector<Int_t > iw(MaxNClus,0);
std::vector<Double_t > dw(MaxNClus,0);
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
fNClust.push_back(0);
fClustSize.push_back(iw);
fClustPos.push_back(dw);
}
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Int_t prev_padnum=-100;
for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
Int_t padnum = hit->GetPaddleNumber();
if ( hit->GetTwoGoodTimes() ) {
if ( padnum==prev_padnum+1 ) {
fClustSize[ip][fNClust[ip]-1]=fClustSize[ip][fNClust[ip]-1]+1;
fClustPos[ip][fNClust[ip]-1]=fClustPos[ip][fNClust[ip]-1]+fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset();
// cout << "Add to cluster pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus = " << fNClust[ip] << " cl size = " << fClustSize[ip][fNClust[ip]-1] << " pos " << fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset() << endl;
} else {
if (fNClust[ip]<MaxNClus) fNClust[ip]++;
fClustSize[ip][fNClust[ip]-1]=1;
fClustPos[ip][fNClust[ip]-1]=fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset();
// cout << " New clus pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus = " << fNClust[ip] << " cl size = " << fClustSize[ip][fNClust[ip]] << " pos " << fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset() << endl;
}
prev_padnum=padnum;
}
}
}
//
Bool_t inside_bound[4]={kFALSE,kFALSE,kFALSE,kFALSE};
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
for(Int_t ic = 0; ic <fNClust[ip] ; ic++ ) {
fClustPos[ip][ic]=fClustPos[ip][ic]/fClustSize[ip][ic];
inside_bound[ip] = fClustPos[ip][ic]>=PadPosLo[ip] && fClustPos[ip][ic]<=PadPosHi[ip];
//cout << "plane = " << ip+1 << " Cluster = " << ic+1 << " size = " << fClustSize[ip][ic]<< " pos = " << fClustPos[ip][ic] << " inside = " << inside_bound[ip] << " lo = " << PadPosLo[ip]<< " hi = " << PadPosHi[ip]<< endl;
}
}
//
Int_t good_for_track_test[4]={0,0,0,0};
Int_t sum_good_track_test=0;
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
if (fNClust[ip]==1 && inside_bound[ip] && fClustSize[ip][0]<=2) good_for_track_test[ip]=1;
//cout << " good for track = " << good_for_track_test[ip] << endl;
sum_good_track_test+=good_for_track_test[ip];
}
//
Double_t trackeff_scint_ydiff_max= 10. ;
Double_t trackeff_scint_xdiff_max= 10. ;
Bool_t xdiffTest=kFALSE;
Bool_t ydiffTest=kFALSE;
fGoodScinHits = 0;
if (fTrackEffTestNScinPlanes == 4) {
if (fTrackEffTestNScinPlanes==sum_good_track_test) {
xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
if (xdiffTest && ydiffTest) fGoodScinHits = 1;
}
}
//
if (fTrackEffTestNScinPlanes == 3) {
if (fTrackEffTestNScinPlanes==sum_good_track_test) {
if(good_for_track_test[0]==1&&good_for_track_test[2]==1) {
xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
ydiffTest=kTRUE;
}
if (good_for_track_test[1]==1&&good_for_track_test[3]==1) {
xdiffTest=kTRUE;
ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
}
if (xdiffTest && ydiffTest) fGoodScinHits = 1;
}
if (sum_good_track_test==4) {
xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
if (xdiffTest && ydiffTest) fGoodScinHits = 1;
}
}
//
// cout << " good scin = " << fGoodScinHits << " " << sum_good_track_test << " " << xdiffTest << " " << ydiffTest<< endl;
//cout << " ************" << endl;
//
}
//
//
//
void THcHodoscope::OriginalTrackEffTest(void)
{
//-----------------------------------------------------------------------
//
// Trnslation 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;
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
// *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();
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);
Bool_t cer_pid=( fCherenkov->GetCerNPE() > fTOFCalib_cer_lo);
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)
////////////////////////////////////////////////////////////////////////////////