Skip to content
Snippets Groups Projects
Commit a48ac69c authored by Stephen A. Wood's avatar Stephen A. Wood Committed by Stephen A. Wood
Browse files

Change space point from a structure to a full class 

  THcDriftChamber::SpacePoint -> THcSpacePoint
  That way space points can be put in a TClonesArray and be made known
  to THcDC when stub information is added.
parent 8ecd5887
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Makefile
+ 1
− 0
View file @ a48ac69c
...@@ -20,6 +20,7 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \ ...@@ -20,6 +20,7 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
src/THcRawDCHit.cxx src/THcDCHit.cxx \ src/THcRawDCHit.cxx src/THcDCHit.cxx \
src/THcDCWire.cxx \ src/THcDCWire.cxx \
src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \ src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \
src/THcSpacePoint.cxx \
src/THcShower.cxx src/THcShowerPlane.cxx \ src/THcShower.cxx src/THcShowerPlane.cxx \
src/THcShowerHit.cxx \ src/THcShowerHit.cxx \
src/THcAerogel.cxx src/THcAerogelHit.cxx src/THcAerogel.cxx src/THcAerogelHit.cxx
......
src/HallC_LinkDef.h
+ 1
− 0
View file @ a48ac69c
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#pragma link C++ class THcDCWire+; #pragma link C++ class THcDCWire+;
#pragma link C++ class THcDCLookupTTDConv+; #pragma link C++ class THcDCLookupTTDConv+;
#pragma link C++ class THcDCTimeToDistConv+; #pragma link C++ class THcDCTimeToDistConv+;
#pragma link C++ class THcSpacePoint+;
#pragma link C++ class THcShower+; #pragma link C++ class THcShower+;
#pragma link C++ class THcShowerPlane+; #pragma link C++ class THcShowerPlane+;
#pragma link C++ class THcShowerHit+; #pragma link C++ class THcShowerHit+;
......
src/THcDriftChamber.cxx
+ 109
− 97
View file @ a48ac69c
...@@ -20,6 +20,7 @@ ...@@ -20,6 +20,7 @@
#include "TClonesArray.h" #include "TClonesArray.h"
#include "TMath.h" #include "TMath.h"
#include "TVectorD.h" #include "TVectorD.h"
#include "THcSpacePoint.h"
#include "THaTrackProj.h" #include "THaTrackProj.h"
...@@ -44,6 +45,8 @@ THcDriftChamber::THcDriftChamber( ...@@ -44,6 +45,8 @@ THcDriftChamber::THcDriftChamber(
fChamberNum = chambernum; fChamberNum = chambernum;
fSpacePoints = new TClonesArray("THcSpacePoint",10);
} }
//_____________________________________________________________________________ //_____________________________________________________________________________
...@@ -253,6 +256,8 @@ Int_t THcDriftChamber::FindSpacePoints( void ) ...@@ -253,6 +256,8 @@ Int_t THcDriftChamber::FindSpacePoints( void )
// Code below is specifically for HMS chambers with Y and Y' planes // Code below is specifically for HMS chambers with Y and Y' planes
fSpacePoints->Clear();
Int_t yplane_hitind=0; Int_t yplane_hitind=0;
Int_t yplanep_hitind=0; Int_t yplanep_hitind=0;
...@@ -372,14 +377,11 @@ Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hiti ...@@ -372,14 +377,11 @@ Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hiti
} }
if(easy_space_point) { // Register the space point if(easy_space_point) { // Register the space point
cout << "Easy Space Point " << xt << " " << yt << endl; cout << "Easy Space Point " << xt << " " << yt << endl;
fNSpacePoints = 1; THcSpacePoint* sp = (THcSpacePoint*)fSpacePoints->ConstructedAt(fNSpacePoints++);
fSpacePoints[0].x = xt; sp->Clear();
fSpacePoints[0].y = yt; sp->SetXY(xt, yt);
fSpacePoints[0].nhits = fNhits;
fSpacePoints[0].ncombos = 0; // No combos
fSpacePoints[0].hits.resize(fNhits);
for(Int_t ihit=0;ihit<fNhits;ihit++) { for(Int_t ihit=0;ihit<fNhits;ihit++) {
fSpacePoints[0].hits[ihit] = fHits[ihit]; sp->AddHit(fHits[ihit]);
} }
} }
return(easy_space_point); return(easy_space_point);
...@@ -454,9 +456,9 @@ Int_t THcDriftChamber::FindHardSpacePoints() ...@@ -454,9 +456,9 @@ Int_t THcDriftChamber::FindHardSpacePoints()
if(fNSpacePoints > 0) { if(fNSpacePoints > 0) {
Int_t add_flag=1; Int_t add_flag=1;
for(Int_t ispace=0;ispace<fNSpacePoints;ispace++) { for(Int_t ispace=0;ispace<fNSpacePoints;ispace++) {
if(fSpacePoints[ispace].nhits > 0) { THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[ispace];
Double_t sqdist_test = pow(xt - fSpacePoints[ispace].x,2) + if(sp->GetNHits() > 0) {
pow(yt - fSpacePoints[ispace].y,2); Double_t sqdist_test = pow(xt - sp->GetX(),2) + pow(yt - sp->GetY(),2);
// I (who is I) want to be careful if sqdist_test is bvetween 1 and // I (who is I) want to be careful if sqdist_test is bvetween 1 and
// 3 fSpacePointCriterion2. Let me ignore not add a new point the // 3 fSpacePointCriterion2. Let me ignore not add a new point the
if(sqdist_test < 3*fSpacePointCriterion2) { if(sqdist_test < 3*fSpacePointCriterion2) {
...@@ -470,9 +472,9 @@ Int_t THcDriftChamber::FindHardSpacePoints() ...@@ -470,9 +472,9 @@ Int_t THcDriftChamber::FindHardSpacePoints()
// Find out which of the four hits in the combo are already // Find out which of the four hits in the combo are already
// in the space point under consideration so that we don't // in the space point under consideration so that we don't
// add duplicate hits to the space point // add duplicate hits to the space point
for(Int_t isp_hit=0;isp_hit<fSpacePoints[ispace].nhits;isp_hit++) { for(Int_t isp_hit=0;isp_hit<sp->GetNHits();isp_hit++) {
for(Int_t icm_hit=0;icm_hit<4;icm_hit++) { // Loop over combo hits for(Int_t icm_hit=0;icm_hit<4;icm_hit++) { // Loop over combo hits
if(fSpacePoints[ispace].hits[isp_hit]==hits[icm_hit]) { if(sp->GetHit(isp_hit)==hits[icm_hit]) {
iflag[icm_hit] = 1; iflag[icm_hit] = 1;
} }
} }
...@@ -489,11 +491,10 @@ Int_t THcDriftChamber::FindHardSpacePoints() ...@@ -489,11 +491,10 @@ Int_t THcDriftChamber::FindHardSpacePoints()
// Add the unique combo hits to the space point // Add the unique combo hits to the space point
for(Int_t icm=0;icm<4;icm++) { for(Int_t icm=0;icm<4;icm++) {
if(iflag[icm]==0) { if(iflag[icm]==0) {
fSpacePoints[ispace].hits.push_back(hits[icm]); sp->AddHit(hits[icm]);
fSpacePoints[ispace].nhits++;
} }
} }
fSpacePoints[ispace].ncombos++; sp->IncCombos();
// Terminate loop since this combo can only belong to one space point // Terminate loop since this combo can only belong to one space point
break; break;
} }
...@@ -502,43 +503,35 @@ Int_t THcDriftChamber::FindHardSpacePoints() ...@@ -502,43 +503,35 @@ Int_t THcDriftChamber::FindHardSpacePoints()
// Create a new space point if more than 2*space_point_criteria // Create a new space point if more than 2*space_point_criteria
if(fNSpacePoints < MAX_SPACE_POINTS) { if(fNSpacePoints < MAX_SPACE_POINTS) {
if(add_flag) { if(add_flag) {
fSpacePoints[fNSpacePoints].nhits=2; THcSpacePoint* sp = (THcSpacePoint*)fSpacePoints->ConstructedAt(fNSpacePoints++);
fSpacePoints[fNSpacePoints].ncombos=1; sp->Clear();
fSpacePoints[fNSpacePoints].hits.resize(fSpacePoints[fNSpacePoints].nhits); sp->SetXY(xt, yt);
fSpacePoints[fNSpacePoints].hits[0]=hits[0]; sp->SetCombos(1);
fSpacePoints[fNSpacePoints].hits[1]=hits[1]; sp->AddHit(hits[0]);
fSpacePoints[fNSpacePoints].x = xt; sp->AddHit(hits[1]);
fSpacePoints[fNSpacePoints].y = yt;
if(hits[0] != hits[2] && hits[1] != hits[2]) { if(hits[0] != hits[2] && hits[1] != hits[2]) {
fSpacePoints[fNSpacePoints].hits.push_back(hits[2]); sp->AddHit(hits[2]);
fSpacePoints[fNSpacePoints].nhits++;
} }
if(hits[0] != hits[3] && hits[1] != hits[3]) { if(hits[0] != hits[3] && hits[1] != hits[3]) {
fSpacePoints[fNSpacePoints].hits.push_back(hits[3]); sp->AddHit(hits[3]);
fSpacePoints[fNSpacePoints].nhits++;
} }
fNSpacePoints++;
} }
} }
} else {// Create first space point } else {// Create first space point
// This duplicates code above. Need to see if we can restructure // This duplicates code above. Need to see if we can restructure
// to avoid // to avoid
fSpacePoints[fNSpacePoints].nhits=2; THcSpacePoint* sp = (THcSpacePoint*)fSpacePoints->ConstructedAt(fNSpacePoints++);
fSpacePoints[fNSpacePoints].ncombos=1; sp->Clear();
fSpacePoints[fNSpacePoints].hits.resize(fSpacePoints[fNSpacePoints].nhits); sp->SetXY(xt, yt);
fSpacePoints[fNSpacePoints].hits[0] = hits[0]; sp->SetCombos(1);
fSpacePoints[fNSpacePoints].hits[1] = hits[1]; sp->AddHit(hits[0]);
fSpacePoints[fNSpacePoints].x = xt; sp->AddHit(hits[1]);
fSpacePoints[fNSpacePoints].y = yt;
if(hits[0] != hits[2] && hits[1] != hits[2]) { if(hits[0] != hits[2] && hits[1] != hits[2]) {
fSpacePoints[fNSpacePoints].hits.push_back(hits[2]); sp->AddHit(hits[2]);
fSpacePoints[fNSpacePoints].nhits++;
} }
if(hits[0] != hits[3] && hits[1] != hits[3]) { if(hits[0] != hits[3] && hits[1] != hits[3]) {
fSpacePoints[fNSpacePoints].hits.push_back(hits[3]); sp->AddHit(hits[3]);
fSpacePoints[fNSpacePoints].nhits++;
} }
fNSpacePoints++;
}//End check on 0 space points }//End check on 0 space points
}//End loop over combos }//End loop over combos
return(fNSpacePoints); return(fNSpacePoints);
...@@ -562,8 +555,9 @@ Int_t THcDriftChamber::DestroyPoorSpacePoints() ...@@ -562,8 +555,9 @@ Int_t THcDriftChamber::DestroyPoorSpacePoints()
nhitsperplane[ip] = 0; nhitsperplane[ip] = 0;
} }
// Count # hits in each plane for this space point // Count # hits in each plane for this space point
for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) { THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
THcDCHit* hit=fSpacePoints[isp].hits[ihit]; for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) {
THcDCHit* hit=sp->GetHit(ihit);
// hit_order(hit) = ihit; // hit_order(hit) = ihit;
Int_t ip = hit->GetPlaneIndex(); Int_t ip = hit->GetPlaneIndex();
nhitsperplane[ip]++; nhitsperplane[ip]++;
...@@ -590,7 +584,7 @@ Int_t THcDriftChamber::DestroyPoorSpacePoints() ...@@ -590,7 +584,7 @@ Int_t THcDriftChamber::DestroyPoorSpacePoints()
if(osp > isp) { if(osp > isp) {
// Does this work, or do we have to copy each member? // Does this work, or do we have to copy each member?
// If it doesn't we should overload the = operator // If it doesn't we should overload the = operator
fSpacePoints[isp] = fSpacePoints[osp]; (*fSpacePoints)[isp] = (*fSpacePoints)[osp];
} }
} }
return nremoved; return nremoved;
...@@ -628,8 +622,9 @@ Int_t THcDriftChamber::SpacePointMultiWire() ...@@ -628,8 +622,9 @@ Int_t THcDriftChamber::SpacePointMultiWire()
nhitsperplane[ip] = 0; nhitsperplane[ip] = 0;
} }
// Sort Space Points hits by plane // Sort Space Points hits by plane
for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) { // All hits in SP THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
THcDCHit* hit=fSpacePoints[isp].hits[ihit]; for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) { // All hits in SP
THcDCHit* hit=sp->GetHit(ihit);
// hit_order Make a hash // hit_order Make a hash
// hash(hit) = ihit; // hash(hit) = ihit;
Int_t ip = hit->GetPlaneIndex(); Int_t ip = hit->GetPlaneIndex();
...@@ -674,20 +669,25 @@ Int_t THcDriftChamber::SpacePointMultiWire() ...@@ -674,20 +669,25 @@ Int_t THcDriftChamber::SpacePointMultiWire()
for(Int_t n3=0;n3<nhitsperplane[maxplane[2]];n3++) { for(Int_t n3=0;n3<nhitsperplane[maxplane[2]];n3++) {
ntot++; ntot++;
newsp_num = nsp_tot + ntot - 2; // ntot will be 2 for first new newsp_num = nsp_tot + ntot - 2; // ntot will be 2 for first new
if(n1==0 && n2==0 && n3==0) newsp_num = isp; // Copy over original SP THcSpacePoint* newsp;
fSpacePoints[newsp_num].x = fSpacePoints[isp].x; if(n1==0 && n2==0 && n3==0) {
fSpacePoints[newsp_num].y = fSpacePoints[isp].y; newsp_num = isp; // Copy over original SP
fSpacePoints[newsp_num].nhits = nplanes_hit; newsp = sp;
fSpacePoints[newsp_num].ncombos = fSpacePoints[isp].ncombos; newsp->Clear(); // Clear doesn't clear X, Y
fSpacePoints[newsp_num].hits.resize(0); } else {
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[0]][n1]); THcSpacePoint* newsp = (THcSpacePoint*)fSpacePoints->ConstructedAt(fNSpacePoints++);
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[1]][n2]); newsp->Clear();
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[2]][n3]); newsp->SetXY(sp->GetX(), sp->GetY());
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[3]][0]); }
newsp->SetCombos(sp->GetCombos());
newsp->AddHit(hits_plane[maxplane[0]][n1]);
newsp->AddHit(hits_plane[maxplane[1]][n2]);
newsp->AddHit(hits_plane[maxplane[2]][n3]);
newsp->AddHit(hits_plane[maxplane[3]][0]);
if(nhitsperplane[maxplane[4]] == 1) { if(nhitsperplane[maxplane[4]] == 1) {
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[4]][0]); newsp->AddHit(hits_plane[maxplane[4]][0]);
if(nhitsperplane[maxplane[5]] == 1) if(nhitsperplane[maxplane[5]] == 1)
fSpacePoints[newsp_num].hits.push_back(hits_plane[maxplane[5]][0]); newsp->AddHit(hits_plane[maxplane[5]][0]);
} }
} }
} }
...@@ -733,21 +733,22 @@ void THcDriftChamber::ChooseSingleHit() ...@@ -733,21 +733,22 @@ void THcDriftChamber::ChooseSingleHit()
{ {
// Look at all hits in a space point. If two hits are in the same plane, // Look at all hits in a space point. If two hits are in the same plane,
// reject the one with the longer drift time. // reject the one with the longer drift time.
Int_t goodhit[MAX_HITS_PER_POINT];
for(Int_t isp=0;isp<fNSpacePoints;isp++) { for(Int_t isp=0;isp<fNSpacePoints;isp++) {
Int_t startnum = fSpacePoints[isp].nhits; THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
Int_t startnum = sp->GetNHits();
Int_t goodhit[startnum];
for(Int_t ihit=0;ihit<startnum;ihit++) { for(Int_t ihit=0;ihit<startnum;ihit++) {
goodhit[ihit] = 1; goodhit[ihit] = 1;
} }
// For each plane, mark all hits longer than the shortest drift time // For each plane, mark all hits longer than the shortest drift time
for(Int_t ihit1=0;ihit1<startnum-1;ihit1++) { for(Int_t ihit1=0;ihit1<startnum-1;ihit1++) {
THcDCHit* hit1 = fSpacePoints[isp].hits[ihit1]; THcDCHit* hit1 = sp->GetHit(ihit1);
Int_t plane1=hit1->GetPlaneIndex(); Int_t plane1=hit1->GetPlaneIndex();
Double_t tdrift1 = hit1->GetTime(); Double_t tdrift1 = hit1->GetTime();
for(Int_t ihit2=ihit1+1;ihit2<startnum;ihit2++) { for(Int_t ihit2=ihit1+1;ihit2<startnum;ihit2++) {
THcDCHit* hit2 = fSpacePoints[isp].hits[ihit2]; THcDCHit* hit2 = sp->GetHit(ihit2);
Int_t plane2=hit2->GetPlaneIndex(); Int_t plane2=hit2->GetPlaneIndex();
Double_t tdrift2 = hit2->GetTime(); Double_t tdrift2 = hit2->GetTime();
if(plane1 == plane2) { if(plane1 == plane2) {
...@@ -764,13 +765,11 @@ void THcDriftChamber::ChooseSingleHit() ...@@ -764,13 +765,11 @@ void THcDriftChamber::ChooseSingleHit()
for(Int_t ihit=0;ihit<startnum;ihit++) { for(Int_t ihit=0;ihit<startnum;ihit++) {
if(goodhit[ihit] > 0) { // Keep this hit if(goodhit[ihit] > 0) { // Keep this hit
if (ihit > finalnum) { // Move hit if (ihit > finalnum) { // Move hit
fSpacePoints[isp].hits[finalnum++] = fSpacePoints[isp].hits[ihit]; sp->ReplaceHit(finalnum++, sp->GetHit(ihit));
} else {
finalnum++;
} }
} }
} }
fSpacePoints[isp].nhits = finalnum; sp->SetNHits(finalnum);
} }
} }
//_____________________________________________________________________________ //_____________________________________________________________________________
...@@ -784,12 +783,17 @@ void THcDriftChamber::SelectSpacePoints() ...@@ -784,12 +783,17 @@ void THcDriftChamber::SelectSpacePoints()
Int_t sp_count=0; Int_t sp_count=0;
for(Int_t isp=0;isp<fNSpacePoints;isp++) { for(Int_t isp=0;isp<fNSpacePoints;isp++) {
// Include fEasySpacePoint because ncombos not filled in // Include fEasySpacePoint because ncombos not filled in
if(fSpacePoints[isp].ncombos >= fMinCombos || fEasySpacePoint) { THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
if(fSpacePoints[isp].nhits >= fMinHits) { if(sp->GetCombos() >= fMinCombos || fEasySpacePoint) {
fSpacePoints[sp_count++] = fSpacePoints[isp]; if(sp->GetNHits() >= fMinHits) {
if(isp > sp_count)
(*fSpacePoints)[sp_count] = (*fSpacePoints)[isp];
sp_count++;
} }
} }
} }
if(sp_count < fNSpacePoints)
cout << "Reduced from " << fNSpacePoints << " to " << sp_count << " space points" << endl;
fNSpacePoints = sp_count; fNSpacePoints = sp_count;
} }
...@@ -805,10 +809,11 @@ void THcDriftChamber::CorrectHitTimes() ...@@ -805,10 +809,11 @@ void THcDriftChamber::CorrectHitTimes()
// time offset per card will cancel much of the error caused by this. The // time offset per card will cancel much of the error caused by this. The
// alternative is to check by card, rather than by plane and this is harder. // alternative is to check by card, rather than by plane and this is harder.
for(Int_t isp=0;isp<fNSpacePoints;isp++) { for(Int_t isp=0;isp<fNSpacePoints;isp++) {
Double_t x = fSpacePoints[isp].x; THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
Double_t y = fSpacePoints[isp].y; Double_t x = sp->GetX();
for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) { Double_t y = sp->GetY();
THcDCHit* hit = fSpacePoints[isp].hits[ihit]; for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) {
THcDCHit* hit = sp->GetHit(ihit);
THcDriftChamberPlane* plane=hit->GetWirePlane(); THcDriftChamberPlane* plane=hit->GetWirePlane();
// How do we know this correction only gets applied once? Is // How do we know this correction only gets applied once? Is
...@@ -850,7 +855,8 @@ void THcDriftChamber::LeftRight() ...@@ -850,7 +855,8 @@ void THcDriftChamber::LeftRight()
for(Int_t isp=0; isp<fNSpacePoints; isp++) { for(Int_t isp=0; isp<fNSpacePoints; isp++) {
// Build a bit pattern of which planes are hit // Build a bit pattern of which planes are hit
Int_t nhits = fSpacePoints[isp].nhits; THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
Int_t nhits = sp->GetNHits();
UInt_t bitpat = 0; // Bit pattern of which planes are hit UInt_t bitpat = 0; // Bit pattern of which planes are hit
Double_t minchi2 = 1.0e10; Double_t minchi2 = 1.0e10;
Double_t tmp_minchi2; Double_t tmp_minchi2;
...@@ -871,7 +877,7 @@ void THcDriftChamber::LeftRight() ...@@ -871,7 +877,7 @@ void THcDriftChamber::LeftRight()
cout << "THcDriftChamber::LeftRight() nhits = 0" << endl; cout << "THcDriftChamber::LeftRight() nhits = 0" << endl;
} }
for(Int_t ihit=0;ihit < nhits;ihit++) { for(Int_t ihit=0;ihit < nhits;ihit++) {
THcDCHit* hit = fSpacePoints[isp].hits[ihit]; THcDCHit* hit = sp->GetHit(ihit);
Int_t pindex = hit->GetPlaneIndex(); Int_t pindex = hit->GetPlaneIndex();
plane_list[ihit] = pindex; plane_list[ihit] = pindex;
...@@ -884,8 +890,8 @@ void THcDriftChamber::LeftRight() ...@@ -884,8 +890,8 @@ void THcDriftChamber::LeftRight()
} }
Int_t smallAngOK = (hasy1>=0) && (hasy2>=0); Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes
if(fSpacePoints[isp].hits[hasy1]->GetPos() <= if(sp->GetHit(hasy1)->GetPos() <=
fSpacePoints[isp].hits[hasy2]->GetPos()) { sp->GetHit(hasy2)->GetPos()) {
plusminusknown[hasy1] = -1; plusminusknown[hasy1] = -1;
plusminusknown[hasy2] = 1; plusminusknown[hasy2] = 1;
} else { } else {
...@@ -919,7 +925,7 @@ void THcDriftChamber::LeftRight() ...@@ -919,7 +925,7 @@ void THcDriftChamber::LeftRight()
} }
} }
if (nplaneshit >= fNPlanes-1) { if (nplaneshit >= fNPlanes-1) {
Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits, Double_t chi2 = FindStub(nhits, sp->GetHitVectorP(),
plane_list, bitpat, plusminus, stub); plane_list, bitpat, plusminus, stub);
//if(debugging) //if(debugging)
...@@ -936,14 +942,15 @@ void THcDriftChamber::LeftRight() ...@@ -936,14 +942,15 @@ void THcDriftChamber::LeftRight()
Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]] Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
/(1+stub[2]*fTanBeta[plane_list[0]]); /(1+stub[2]*fTanBeta[plane_list[0]]);
// Tune depdendent. Definitely HMS specific // Tune depdendent. Definitely HMS specific
Double_t xp_expect = fSpacePoints[isp].x/875.0; Double_t xp_expect = sp->GetX()/875.0;
if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) { if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
minchi2 = chi2; minchi2 = chi2;
for(Int_t ihit=0;ihit<nhits;ihit++) { for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit]; plusminusbest[ihit] = plusminus[ihit];
} }
Double_t *spstub = sp->GetStubP();
for(Int_t i=0;i<4;i++) { for(Int_t i=0;i<4;i++) {
fSpacePoints[isp].stub[i] = stub[i]; spstub[i] = stub[i];
} }
} else { // Record best stub failing angle cut } else { // Record best stub failing angle cut
tmp_minchi2 = chi2; tmp_minchi2 = chi2;
...@@ -956,7 +963,7 @@ void THcDriftChamber::LeftRight() ...@@ -956,7 +963,7 @@ void THcDriftChamber::LeftRight()
} }
} }
} else if (nplaneshit >= fNPlanes-2) { // Two planes missing } else if (nplaneshit >= fNPlanes-2) { // Two planes missing
Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits, Double_t chi2 = FindStub(nhits, sp->GetHitVectorP(),
plane_list, bitpat, plusminus, stub); plane_list, bitpat, plusminus, stub);
//if(debugging) //if(debugging)
//cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl; //cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
...@@ -972,8 +979,9 @@ void THcDriftChamber::LeftRight() ...@@ -972,8 +979,9 @@ void THcDriftChamber::LeftRight()
for(Int_t ihit=0;ihit<nhits;ihit++) { for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = plusminus[ihit]; plusminusbest[ihit] = plusminus[ihit];
} }
Double_t *spstub = sp->GetStubP();
for(Int_t i=0;i<4;i++) { for(Int_t i=0;i<4;i++) {
fSpacePoints[isp].stub[i] = stub[i]; spstub[i] = stub[i];
} }
} }
} else { } else {
...@@ -981,13 +989,14 @@ void THcDriftChamber::LeftRight() ...@@ -981,13 +989,14 @@ void THcDriftChamber::LeftRight()
} }
} // End loop of pm combinations } // End loop of pm combinations
Double_t *spstub = sp->GetStubP();
if(minchi2 > 9.9e9) { // No track passed angle cut if(minchi2 > 9.9e9) { // No track passed angle cut
minchi2 = tmp_minchi2; minchi2 = tmp_minchi2;
for(Int_t ihit=0;ihit<nhits;ihit++) { for(Int_t ihit=0;ihit<nhits;ihit++) {
plusminusbest[ihit] = tmp_plusminus[ihit]; plusminusbest[ihit] = tmp_plusminus[ihit];
} }
for(Int_t i=0;i<4;i++) { for(Int_t i=0;i<4;i++) {
fSpacePoints[isp].stub[i] = tmp_stub[i]; spstub[i] = tmp_stub[i];
} }
} }
...@@ -995,35 +1004,35 @@ void THcDriftChamber::LeftRight() ...@@ -995,35 +1004,35 @@ void THcDriftChamber::LeftRight()
// Calculate final coordinate based on plusminusbest // Calculate final coordinate based on plusminusbest
// Update the hit positions in the space points // Update the hit positions in the space points
for(Int_t ihit; ihit<nhits; ihit++) { for(Int_t ihit; ihit<nhits; ihit++) {
fSpacePoints[isp].hits[ihit]->SetLeftRight(plusminusbest[ihit]); sp->GetHit(ihit)->SetLeftRight(plusminusbest[ihit]);
} }
// Stubs are calculated in rotated coordinate system // Stubs are calculated in rotated coordinate system
// (I think this rotates in case chambers not perpendicular to central ray) // (I think this rotates in case chambers not perpendicular to central ray)
Int_t pindex=plane_list[0]; Int_t pindex=plane_list[0];
if(fSpacePoints[isp].stub[2] - fTanBeta[pindex] == -1.0) { if(spstub[2] - fTanBeta[pindex] == -1.0) {
cout << "THcDriftChamber::LeftRight(): stub3 error" << endl; cout << "THcDriftChamber::LeftRight(): stub3 error" << endl;
} }
stub[2] = (fSpacePoints[isp].stub[2] - fTanBeta[pindex]) stub[2] = (spstub[2] - fTanBeta[pindex])
/ (1.0 + fSpacePoints[isp].stub[2]*fTanBeta[pindex]); / (1.0 + spstub[2]*fTanBeta[pindex]);
if(fSpacePoints[isp].stub[2]*fSinBeta[pindex] == -fCosBeta[pindex]) { if(spstub[2]*fSinBeta[pindex] == -fCosBeta[pindex]) {
cout << "THcDriftChamber::LeftRight(): stub4 error" << endl; cout << "THcDriftChamber::LeftRight(): stub4 error" << endl;
} }
stub[3] = fSpacePoints[isp].stub[3] stub[3] = spstub[3]
/ (fSpacePoints[isp].stub[2]*fSinBeta[pindex]+fCosBeta[pindex]); / (spstub[2]*fSinBeta[pindex]+fCosBeta[pindex]);
stub[0] = fSpacePoints[isp].stub[0]*fCosBeta[pindex] stub[0] = spstub[0]*fCosBeta[pindex]
- fSpacePoints[isp].stub[0]*stub[2]*fSinBeta[pindex]; - spstub[0]*stub[2]*fSinBeta[pindex];
stub[1] = fSpacePoints[isp].stub[1] stub[1] = spstub[1]
- fSpacePoints[isp].stub[1]*stub[3]*fSinBeta[pindex]; - spstub[1]*stub[3]*fSinBeta[pindex];
for(Int_t i=0;i<4;i++) { for(Int_t i=0;i<4;i++) {
fSpacePoints[isp].stub[i] = stub[i]; spstub[i] = stub[i];
} }
} }
// Option to print stubs // Option to print stubs
} }
// if(fAA3Inv.find(bitpat) != fAAInv.end()) { // Valid hit combination // if(fAA3Inv.find(bitpat) != fAAInv.end()) { // Valid hit combination
//_____________________________________________________________________________ //_____________________________________________________________________________
Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>& hits, Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>* hits,
Int_t* plane_list, UInt_t bitpat, Int_t* plane_list, UInt_t bitpat,
Int_t* plusminus, Double_t* stub) Int_t* plusminus, Double_t* stub)
{ {
...@@ -1031,10 +1040,11 @@ Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>& hi ...@@ -1031,10 +1040,11 @@ Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>& hi
Double_t zeros[] = {0.0,0.0,0.0}; Double_t zeros[] = {0.0,0.0,0.0};
TVectorD TT; TT.Use(3, zeros); TVectorD TT; TT.Use(3, zeros);
TVectorD dstub; TVectorD dstub;
Double_t dpos[3]; Double_t dpos[nhits];
// This isn't right. dpos only goes up to 2.
for(Int_t ihit=0;ihit<nhits; ihit++) { for(Int_t ihit=0;ihit<nhits; ihit++) {
dpos[ihit] = hits[ihit]->GetPos() + plusminus[ihit]*hits[ihit]->GetdDist() dpos[ihit] = (*hits)[ihit]->GetPos() + plusminus[ihit]*(*hits)[ihit]->GetdDist()
- fPsi0[plane_list[ihit]]; - fPsi0[plane_list[ihit]];
for(Int_t index=0;index<3;index++) { for(Int_t index=0;index<3;index++) {
TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index] TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
...@@ -1064,6 +1074,8 @@ THcDriftChamber::~THcDriftChamber() ...@@ -1064,6 +1074,8 @@ THcDriftChamber::~THcDriftChamber()
{ {
// Destructor. Remove variables from global list. // Destructor. Remove variables from global list.
delete fSpacePoints;
if( fIsSetup ) if( fIsSetup )
RemoveVariables(); RemoveVariables();
if( fIsInit ) if( fIsInit )
......
src/THcDriftChamber.h
+ 7
− 2
View file @ a48ac69c
...@@ -105,12 +105,13 @@ protected: ...@@ -105,12 +105,13 @@ protected:
void SelectSpacePoints(void); void SelectSpacePoints(void);
UInt_t Count1Bits(UInt_t x); UInt_t Count1Bits(UInt_t x);
void LeftRight(void); void LeftRight(void);
Double_t FindStub(Int_t nhits, const std::vector<THcDCHit*>& hits, Double_t FindStub(Int_t nhits, const std::vector<THcDCHit*>* hits,
Int_t* plane_list, UInt_t bitpat, Int_t* plane_list, UInt_t bitpat,
Int_t* plusminus, Double_t* stub); Int_t* plusminus, Double_t* stub);
std::vector<THcDCHit*> fHits; /* All hits for this chamber */ std::vector<THcDCHit*> fHits; /* All hits for this chamber */
// A simple structure until we figure out what we are doing. // A simple structure until we figure out what we are doing.
#if 0
struct SpacePoint { struct SpacePoint {
Double_t x; Double_t x;
Double_t y; Double_t y;
...@@ -119,8 +120,12 @@ protected: ...@@ -119,8 +120,12 @@ protected:
// THcDCHit* hits[MAX_HITS_PER_POINT]; // THcDCHit* hits[MAX_HITS_PER_POINT];
std::vector<THcDCHit*> hits; std::vector<THcDCHit*> hits;
Double_t stub[4]; Double_t stub[4];
void Clear(Option_t* opt="") {nhits=0;ncombos=0;hits.clear();};
SpacePoint() {nhits=0;ncombos=0;hits.clear();};
void SetXY(Double_t xa, Double_t ya) {x = xa; y = ya;};
}; };
SpacePoint fSpacePoints[MAX_SPACE_POINTS]; #endif
TClonesArray *fSpacePoints;
Int_t fNSpacePoints; Int_t fNSpacePoints;
Int_t fEasySpacePoint; /* This event is an easy space point */ Int_t fEasySpacePoint; /* This event is an easy space point */
......
src/THcDriftChamberPlane.cxx
+ 6
− 0
View file @ a48ac69c
...@@ -41,6 +41,12 @@ THcDriftChamberPlane::THcDriftChamberPlane( const char* name, ...@@ -41,6 +41,12 @@ THcDriftChamberPlane::THcDriftChamberPlane( const char* name,
fPlaneNum = planenum; fPlaneNum = planenum;
} }
//_____________________________________________________________________________
THcDriftChamberPlane::THcDriftChamberPlane() :
THaSubDetector()
{
// Constructor
}
//______________________________________________________________________________ //______________________________________________________________________________
THcDriftChamberPlane::~THcDriftChamberPlane() THcDriftChamberPlane::~THcDriftChamberPlane()
{ {
......
src/THcDriftChamberPlane.h
+ 1
− 0
View file @ a48ac69c
...@@ -66,6 +66,7 @@ class THcDriftChamberPlane : public THaSubDetector { ...@@ -66,6 +66,7 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t GetPsi0() { return fPsi0; } Double_t GetPsi0() { return fPsi0; }
Double_t* GetStubCoef() { return fStubCoef; } Double_t* GetStubCoef() { return fStubCoef; }
THcDriftChamberPlane(); // for ROOT I/O
protected: protected:
TClonesArray* fParentHitList; TClonesArray* fParentHitList;
......
src/THcSpacePoint.cxx 0 → 100644
+ 13
− 0
View file @ a48ac69c
///////////////////////////////////////////////////////////////////////////////
// //
// THcSpacePoint //
// //
// Class representing a single hit for the VDC //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THcSpacePoint.h"
ClassImp(THcSpacePoint)
///////////////////////////////////////////////////////////////////////////////
src/THcSpacePoint.h 0 → 100644
+ 60
− 0
View file @ a48ac69c
#ifndef ROOT_THcSpacePoint
#define ROOT_THcSpacePoint
///////////////////////////////////////////////////////////////////////////////
// //
// THcSpacePoint //
// //
///////////////////////////////////////////////////////////////////////////////
#include "TObject.h"
#include "THcDCHit.h"
class THcSpacePoint : public TObject {
public:
THcSpacePoint(Int_t nhits=0, Int_t ncombos=0) :
fNHits(nhits), fNCombos(ncombos) {
fHits.clear();
}
virtual ~THcSpacePoint() {}
void SetXY(Double_t x, Double_t y) {fX = x; fY = y;};
void Clear(Option_t* opt="") {fNHits=0; fNCombos=0; fHits.clear();};
void AddHit(THcDCHit* hit) {
fHits.push_back(hit);
fNHits++;
}
Int_t GetNHits() {return fNHits;};
void SetNHits(Int_t nhits) {fNHits = nhits;};
Double_t GetX() {return fX;};
Double_t GetY() {return fY;};
THcDCHit* GetHit(Int_t ihit) {return fHits[ihit];};
std::vector<THcDCHit*>* GetHitVectorP() {return &fHits;};
void ReplaceHit(Int_t ihit, THcDCHit *hit) {fHits[ihit] = hit;};
void SetStub(Double_t stub[4]) {
for(Int_t i=0;i<4;i++) {
fStub[i] = stub[i];
}
};
Double_t* GetStubP() { return fStub; };
void IncCombos() { fNCombos++; };
void SetCombos(Int_t ncombos) { fNCombos=ncombos; };
Int_t GetCombos() { return fNCombos; };
protected:
Double_t fX;
Double_t fY;
Int_t fNHits;
Int_t fNCombos;
std::vector<THcDCHit*> fHits;
Double_t fStub[4];
ClassDef(THcSpacePoint,0) // Drift Chamber class
};
////////////////////////////////////////////////////////////////////////////////
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment