Newer
Older
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;
}
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
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==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 (fhdebugflagpr) cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
} else if (nhits == 2) {
if (fhdebugflagpr) 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(chi2 < minchi2) {
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
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)
////////////////////////////////////////////////////////////////////////////////