Newer
Older
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) { // Adjacent plane
if(hit2->GetWireNum() <= hit1->GetWireNum()) {
plusminusknown[ihit1] = -1;
plusminusknown[ihit2] = 1;
} else {
plusminusknown[ihit1] = 1;
plusminusknown[ihit2] = -1;
}
npaired+=2;
}
}
}
}
nplusminus = 1 << (nhits-npaired);
}
}
if(nhits < 2) {
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) {
Double_t chi2 = FindStub(nhits, sp,
plane_list, bitpat, plusminus, stub);
if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl;
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
if(fHMSStyleChambers) { // Perhaps a different flag here
// 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]]);
// Tune depdendent. Definitely HMS specific
Double_t xp_expect = sp->GetX()/875.0;
if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit];
}
Double_t *spstub = sp->GetStubP();
for(Int_t i=0;i<4;i++) {
spstub[i] = stub[i];
}
} 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];
}
}
} else { // Not HMS specific
minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit];
}
Double_t *spstub = sp->GetStubP();
for(Int_t i=0;i<4;i++) {
spstub[i] = stub[i];
} else if (nplaneshit >= fNPlanes-2 && !fHMSStyleChambers) { // Two planes missing
Double_t chi2 = FindStub(nhits, sp,
plane_list, bitpat, plusminus, stub);
//if(debugging)
//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];
}
Double_t *spstub = sp->GetStubP();
for(Int_t i=0;i<4;i++) {
spstub[i] = stub[i];
}
}
} else {
if (fhdebugflagpr) cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
}
} // End loop of pm combinations
Double_t *spstub = sp->GetStubP();
if(minchi2 > 9.9e9) { // No track passed angle cut
minchi2 = tmp_minchi2;
for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = tmp_plusminus[ihit];
}
for(Int_t i=0;i<4;i++) {
spstub[i] = tmp_stub[i];
}
}
// 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 space 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];
for(Int_t i=0;i<4;i++) {
spstub[i] = stub[i];
//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->GetHit(ihit)->GetDist()
- 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.
//_____________________________________________________________________________
inline
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)
////////////////////////////////////////////////////////////////////////////////