Newer
Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
case 2://readout from right of chamber
if (y>yc){wireDistance = -readhoriz*wireDistance;}
else{wireDistance = readhoriz*wireDistance;}
break;
case 3: //readout from bottom of chamber
if (xc>x){wireDistance= -readvert*wireDistance;}
else{wireDistance = readvert*wireDistance;}
break;
case 4: //readout from left of chamber
if(yc>y){wireDistance = -readhoriz*wireDistance;}
else{wireDistance = readhoriz*wireDistance;}
break;
default:
wireDistance = 0.0;
}
Double_t timeCorrection = wireDistance/fWireVelocity;
if(fFixPropagationCorrection==0) { // ENGINE behavior
Double_t time=hit->GetTime();
hit->SetTime(time - timeCorrection);
hit->ConvertTimeToDist();
} else {
Double_t time=hit->GetTime();
Double_t dist=hit->GetDist();
hit->SetTime(time - timeCorrection);
//double usingOldTime = time-time_corr;
//hit->SetTime(time- time_corr);
hit->ConvertTimeToDist();
sp->SetHitDist(ihit, hit->GetDist());
hit->SetTime(time); // Restore time
hit->SetDist(dist); // Restore distance
}
} else {
/////////////////////////////////////////////////////////////
// if (fhdebugflagpr) cout << "Correcting hit " << hit << " " << plane->GetPlaneNum() << " " << isp << "/" << ihit << " " << x << "," << y << endl;
// Fortran ENGINE does not do this check, so hits can get "corrected"
// multiple times if they belong to multiple space points.
// To reproduce the precise ENGINE behavior, remove this if condition.
if(fFixPropagationCorrection==0) { // ENGINE behavior
hit->SetTime(hit->GetTime() - plane->GetCentralTime()
+ plane->GetDriftTimeSign()*time_corr);
hit->ConvertTimeToDist();
// hit->SetCorrectedStatus(1);
} else {
// New behavior: Save corrected distance with the hit in the space point
// so that the same hit can have a different correction depending on
// which space point it is in.
//
// This is a hack now because the converttimetodist method is connected to the hit
// so I compute the corrected time and distance, and then restore the original
// time and distance. Can probably add a method to hit that does a conversion on a time
// but does not modify the hit data.
Double_t time=hit->GetTime();
Double_t dist=hit->GetDist();
hit->SetTime(time - plane->GetCentralTime()
+ plane->GetDriftTimeSign()*time_corr);
hit->ConvertTimeToDist();
sp->SetHitDist(ihit, hit->GetDist());
hit->SetTime(time); // Restore time
hit->SetDist(dist); // Restore distance
}
}
}
Stephen A. Wood
committed
}
UInt_t THcDriftChamber::Count1Bits(UInt_t x)
// From http://graphics.stanford.edu/~seander/bithacks.html
x = x - ((x >> 1) & 0x55555555);
x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
return (((x + (x >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
}
//_____________________________________________________________________________
void THcDriftChamber::LeftRight()
{
/**
For each space point,
Fit stubs to all possible left-right combinations of drift distances
and choose the set with the minimum chi**2.
*/
for(Int_t isp=0; isp<fNSpacePoints; isp++) {
// Build a bit pattern of which planes are hit
THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
Int_t nhits = sp->GetNHits();
UInt_t bitpat = 0; // Bit pattern of which planes are hit
Double_t maxchi2= 1.0e10;
Double_t minchi2 = maxchi2;
Double_t tmp_minchi2=maxchi2;
Double_t minxp = 0.25;
Int_t hasy1 = -1;
Int_t hasy2 = -1;
Int_t plusminusknown[nhits];
Int_t plusminusbest[nhits];
Int_t plusminus[nhits]; // ENGINE makes this array float. Why?
Int_t tmp_plusminus[nhits];
Int_t plane_list[nhits];
Double_t stub[4];
Double_t tmp_stub[4];
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() nhits < 0" << endl;
} else if (nhits==0) {
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() nhits = 0" << endl;
}
for(Int_t ihit=0;ihit < nhits;ihit++) {
THcDCHit* hit = sp->GetHit(ihit);
Int_t pindex = hit->GetPlaneIndex();
plane_list[ihit] = pindex;
bitpat |= 1<<pindex;
plusminusknown[ihit] = 0;
if(pindex == YPlaneInd) hasy1 = ihit;
if(pindex == YPlanePInd) hasy2 = ihit;
}
if(fHMSStyleChambers) {
Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes
if(sp->GetHit(hasy2)->GetPos() <=
sp->GetHit(hasy1)->GetPos()) {
plusminusknown[hasy1] = -1;
plusminusknown[hasy2] = 1;
} else {
plusminusknown[hasy1] = 1;
plusminusknown[hasy2] = -1;
}
nplusminus = 1<<(nhits-2);
// if (fhdebugflagpr) cout << " Small angle approx = " << smallAngOK << " " << plusminusknown[hasy1] << endl;
//if (fhdebugflagpr) cout << "pm = " << plusminusknown[hasy2] << " " << hasy1 << " " << hasy2 << endl;
//if (fhdebugflagpr) cout << " Plane index " << YPlaneInd << " " << YPlanePInd << endl;
}
} else { // SOS Style
if(fSmallAngleApprox !=0) {
// Brookhaven style chamber L/R code
Int_t npaired=0;
for(Int_t ihit1=0;ihit1 < nhits;ihit1++) {
THcDCHit* hit1 = sp->GetHit(ihit1);
Int_t pindex1=hit1->GetPlaneIndex();
if((pindex1%2)==0) { // Odd plane (or even index)
for(Int_t ihit2=0;ihit2<nhits;ihit2++) {
THcDCHit* hit2 = sp->GetHit(ihit2);
if(hit2->GetPlaneIndex()-pindex1 == 1 && TMath::Abs(hit2->GetPos()-hit1->GetPos())<0.51) { // Adjacent plane
if(hit2->GetPos() <= hit1->GetPos() ) {
plusminusknown[ihit1] = -1;
plusminusknown[ihit2] = 1;
} else {
plusminusknown[ihit1] = 1;
plusminusknown[ihit2] = -1;
}
npaired+=2;
}
}
}
}
nplusminus = 1 << (nhits-npaired);
}//end if fSmallAngleApprox!=0
}//end else SOS style
if (fdebugstubchisq) cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
} else if (nhits == 2) {
if (fdebugstubchisq) cout << "THcDriftChamber::LeftRight: numhits-2 = 0" << endl;
}
Int_t nplaneshit = Count1Bits(bitpat);
//if (fhdebugflagpr) cout << " num of pm = " << nplusminus << " num of hits =" << nhits << endl;
// Use bit value of integer word to set + or -
// Loop over all combinations of left right.
for(Int_t pmloop=0;pmloop<nplusminus;pmloop++) {
Int_t iswhit = 1;
for(Int_t ihit=0;ihit<nhits;ihit++) {
if(plusminusknown[ihit]!=0) {
plusminus[ihit] = plusminusknown[ihit];
} else {
// Max hits per point has to be less than 32.
if(pmloop & iswhit) {
plusminus[ihit] = 1;
} else {
plusminus[ihit] = -1;
}
iswhit <<= 1;
}
}
if ( (nplaneshit >= fNPlanes-1) || (nplaneshit >= fNPlanes-2 && !fHMSStyleChambers)) {
chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub);
if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl;
// Take best chi2 IF x' of the stub agrees with x' as expected from x.
// Sometimes an incorrect x' gives a good chi2 for the stub, even though it is
// not the correct left/right combination for the real track.
// Rotate x'(=stub(3)) to hut coordinates and compare to x' expected from x.
// THIS ASSUMES STANDARD HMS TUNE!!!!, for which x' is approx. x/875.
if(stub[2]*fTanBeta[plane_list[0]]==-1.0) {
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() Error 3" << endl;
}
Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
/(1+stub[2]*fTanBeta[plane_list[0]]);
Double_t xp_expect = sp->GetX()*fRatio_xpfp_to_xfp;
if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit];
}
sp->SetStub(stub);
} else { // Record best stub failing angle cut
tmp_minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
tmp_plusminus[ihit] = plusminus[ihit];
}
for(Int_t i=0;i<4;i++) {
tmp_stub[i] = stub[i];
}
minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit];
}
sp->SetStub(stub);
///////////////
} else if (nplaneshit >= fNPlanes-2 && fHMSStyleChambers) { // Two planes missing
Double_t chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub);
//if (fhdebugflagpr) cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
// Isn't this a bad idea, doing == with reals
if(stub[2]*fTanBeta[plane_list[0]] == -1.0) {
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight() Error 3" << endl;
}
Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
/(1+stub[2]*fTanBeta[plane_list[0]]);
if(TMath::Abs(xp_fit) <= minxp) {
minxp = TMath::Abs(xp_fit);
minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit];
}
sp->SetStub(stub);
if (fhdebugflagpr) cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
}
} // End loop of pm combinations
if (minchi2 == maxchi2 && tmp_minchi2 == maxchi2) {
} else {
if(minchi2 == maxchi2 ) { // No track passed angle cut
minchi2 = tmp_minchi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = tmp_plusminus[ihit];
}
sp->SetStub(tmp_stub);
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
1302
Double_t *spstub = sp->GetStubP();
// Calculate final coordinate based on plusminusbest
// Update the hit positions in the space points
for(Int_t ihit=0; ihit<nhits; ihit++) {
// Save left/right status with the hit and in the spaleftce point
// In THcDC will decide which to used based on fix_lr flag
sp->GetHit(ihit)->SetLeftRight(plusminusbest[ihit]);
sp->SetHitLR(ihit, plusminusbest[ihit]);
}
// Stubs are calculated in rotated coordinate system
// (I think this rotates in case chambers not perpendicular to central ray)
Int_t pindex=plane_list[0];
if(spstub[2] - fTanBeta[pindex] == -1.0) {
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight(): stub3 error" << endl;
}
stub[2] = (spstub[2] - fTanBeta[pindex])
/ (1.0 + spstub[2]*fTanBeta[pindex]);
if(spstub[2]*fSinBeta[pindex] == -fCosBeta[pindex]) {
if (fhdebugflagpr) cout << "THcDriftChamber::LeftRight(): stub4 error" << endl;
}
stub[3] = spstub[3]
/ (spstub[2]*fSinBeta[pindex]+fCosBeta[pindex]);
stub[0] = spstub[0]*fCosBeta[pindex]
- spstub[0]*stub[2]*fSinBeta[pindex];
stub[1] = spstub[1]
- spstub[1]*stub[3]*fSinBeta[pindex];
sp->SetStub(stub);
//if (fhdebugflagpr) cout << " Left/Right space pt " << isp+1 << " " << stub[0]<< " " << stub[1] << " " << stub[2]<< " " << stub[3] << endl;
// Option to print stubs
}
//_____________________________________________________________________________
Double_t THcDriftChamber::FindStub(Int_t nhits, THcSpacePoint *sp,
Int_t* plane_list, UInt_t bitpat,
Int_t* plusminus, Double_t* stub)
{
// For a given combination of L/R, fit a stub to the space point
// This method does a linear least squares fit of a line to the
// hits in an individual chamber. It assumes that the y slope is 0
// The wire coordinate is calculated by
// wire center + plusminus*(drift distance).
// Method is called in a loop over all combinations of plusminus
Double_t zeros[] = {0.0,0.0,0.0};
TVectorD TT; TT.Use(3, zeros); // X, X', Y
Double_t dpos[nhits];
for(Int_t ihit=0;ihit<nhits; ihit++) {
dpos[ihit] = sp->GetHit(ihit)->GetPos()
+ plusminus[ihit]*sp->GetHitDist(ihit)
- fPsi0[plane_list[ihit]];
for(Int_t index=0;index<3;index++) {
TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
/fSigma[plane_list[ihit]];
}
}
// fAA3Inv[bitpat].Print();
// if (fhdebugflagpr) cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
// TT.Print();
TT *= fAA3Inv[bitpat];
// if (fhdebugflagpr) cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
// Calculate Chi2. Remember one power of sigma is in fStubCoefs
stub[0] = TT[0];
stub[1] = TT[1];
stub[2] = TT[2];
stub[3] = 0.0;
Double_t chi2=0.0;
for(Int_t ihit=0;ihit<nhits; ihit++) {
chi2 += pow( dpos[ihit]/fSigma[plane_list[ihit]]
- fStubCoefs[plane_list[ihit]][0]*stub[0]
- fStubCoefs[plane_list[ihit]][1]*stub[1]
- fStubCoefs[plane_list[ihit]][2]*stub[2]
, 2);
}
return(chi2);
}
//_____________________________________________________________________________
{
// Destructor. Remove variables from global list.
if (fhdebugflagpr) cout << fChamberNum << ": delete fSpacePoints: " << hex << fSpacePoints << endl;
delete fSpacePoints;
// Should delete the matrices
if( fIsSetup )
RemoveVariables();
if( fIsInit )
DeleteArrays();
if (fTrackProj) {
fTrackProj->Clear();
delete fTrackProj; fTrackProj = 0;
}
}
//_____________________________________________________________________________
void THcDriftChamber::DeleteArrays()
{
// Delete member arrays. Used by destructor.
delete [] fCosBeta; fCosBeta = NULL;
delete [] fSinBeta; fSinBeta = NULL;
delete [] fTanBeta; fTanBeta = NULL;
delete [] fSigma; fSigma = NULL;
delete [] fPsi0; fPsi0 = NULL;
delete [] fStubCoefs; fStubCoefs = NULL;
//_____________________________________________________________________________
void THcDriftChamber::Clear( const Option_t* )
{
// Reset per-event data.
// fNhits = 0;
// fTrackProj->Clear();
fNhits = 0;
}
//_____________________________________________________________________________
Int_t THcDriftChamber::ApplyCorrections( void )
{
return(0);
}
ClassImp(THcDriftChamber)
////////////////////////////////////////////////////////////////////////////////