diff --git a/src/THcDC.cxx b/src/THcDC.cxx index 9ec08cf5317a1bffe394515bea224a8bed0f36e0..67a9039831a11de00efb5d40ac25fcd218186bf9 100644 --- a/src/THcDC.cxx +++ b/src/THcDC.cxx @@ -657,21 +657,21 @@ void THcDC::SetFocalPlaneBestTrack(Int_t golden_track_index) // void THcDC::EfficiencyPerWire(Int_t golden_track_index) { - THcDCTrack *tr1 = static_cast<THcDCTrack*>( fDCTracks->At(golden_track_index)); - Double_t track_pos; - for (UInt_t ihit = 0; ihit < UInt_t (tr1->GetNHits()); ihit++) { - THcDCHit *hit = tr1->GetHit(ihit); - Int_t plane = hit->GetPlaneNum() - 1; - track_pos=tr1->GetCoord(plane); - Int_t wire_num = hit->GetWireNum(); - Int_t wire_track_num=round(fPlanes[plane]->CalcWireFromPos(track_pos)); - if ( (wire_num-wire_track_num) ==0) fWire_hit_did[plane]=wire_num; - } - for(Int_t ip=0; ip<fNPlanes;ip++) { - track_pos=tr1->GetCoord(ip); - Int_t wire_should = round(fPlanes[ip]->CalcWireFromPos(track_pos)); - fWire_hit_should[ip]=wire_should; - } + THcDCTrack *tr1 = static_cast<THcDCTrack*>( fDCTracks->At(golden_track_index)); + Double_t track_pos; + for (UInt_t ihit = 0; ihit < UInt_t (tr1->GetNHits()); ihit++) { + THcDCHit *hit = tr1->GetHit(ihit); + Int_t plane = hit->GetPlaneNum() - 1; + track_pos=tr1->GetCoord(plane); + Int_t wire_num = hit->GetWireNum(); + Int_t wire_track_num=round(fPlanes[plane]->CalcWireFromPos(track_pos)); + if ( (wire_num-wire_track_num) ==0) fWire_hit_did[plane]=wire_num; + } + for(Int_t ip=0; ip<fNPlanes;ip++) { + track_pos=tr1->GetCoord(ip); + Int_t wire_should = round(fPlanes[ip]->CalcWireFromPos(track_pos)); + fWire_hit_should[ip]=wire_should; + } } // void THcDC::PrintSpacePoints() @@ -682,13 +682,13 @@ void THcDC::PrintSpacePoints() printf("%6s %-8s %-8s %6s %6s %10s \n","Point","x","y"," hits ","combos"," for each hit"); TClonesArray* spacepointarray = fChambers[ich]->GetSpacePointsP(); for(Int_t isp=0;isp<fChambers[ich]->GetNSpacePoints();isp++) { - THcSpacePoint* sp = (THcSpacePoint*)(spacepointarray->At(isp)); - printf("%5d %8.5f %8.5f %5d %5d ",isp+1,sp->GetX(),sp->GetY(),sp->GetNHits(),sp->GetCombos()) ; - for (Int_t ii=0;ii<sp->GetNHits();ii++) { - THcDCHit* hittemp = (THcDCHit*)(sp->GetHit(ii)); - printf("%3d %3d %3d",hittemp->GetPlaneNum(),hittemp->GetWireNum(),hittemp->GetLR()); - } - printf("\n"); + THcSpacePoint* sp = (THcSpacePoint*)(spacepointarray->At(isp)); + printf("%5d %8.5f %8.5f %5d %5d ",isp+1,sp->GetX(),sp->GetY(),sp->GetNHits(),sp->GetCombos()) ; + for (Int_t ii=0;ii<sp->GetNHits();ii++) { + THcDCHit* hittemp = (THcDCHit*)(sp->GetHit(ii)); + printf("%3d %3d %3d",hittemp->GetPlaneNum(),hittemp->GetWireNum(),hittemp->GetLR()); + } + printf("\n"); } } } @@ -702,9 +702,9 @@ void THcDC::PrintStubs() printf("%-5s %-18s %-18s %-18s %-18s\n"," ","[cm]","[cm]","[cm]","[cm]"); TClonesArray* spacepointarray = fChambers[ich]->GetSpacePointsP(); for(Int_t isp=0;isp<fChambers[ich]->GetNSpacePoints();isp++) { - THcSpacePoint* sp = (THcSpacePoint*)(spacepointarray->At(isp)); - Double_t *spstubt=sp->GetStubP(); - printf("%-5d % 15.10e % 15.10e % 15.10e % 15.10e \n",isp+1,spstubt[0],spstubt[1],spstubt[2],spstubt[3]); + THcSpacePoint* sp = (THcSpacePoint*)(spacepointarray->At(isp)); + Double_t *spstubt=sp->GetStubP(); + printf("%-5d % 15.10e % 15.10e % 15.10e % 15.10e \n",isp+1,spstubt[0],spstubt[1],spstubt[2],spstubt[3]); } } } @@ -894,8 +894,8 @@ void THcDC::LinkStubs() if(fNDCTracks<MAXTRACKS) { // Need some constructed t thingy if (fSp[isp]->GetSetStubFlag()) { - THcDCTrack *newDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes); - newDCTrack->AddSpacePoint(fSp[isp]); + THcDCTrack *newDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes); + newDCTrack->AddSpacePoint(fSp[isp]); } } else { if (fdebuglinkstubs) cout << "EPIC FAIL 3: Too many tracks found in THcDC::LinkStubs" << endl; @@ -906,17 +906,17 @@ void THcDC::LinkStubs() } } /// -if (fdebuglinkstubs) { - cout << " Number of tracks from link stubs = " << fNDCTracks << endl; - printf("%s %s \n","Track","Plane Wire "); - for (UInt_t itrack=0;itrack<fNDCTracks;itrack++) { - THcDCTrack *tempTrack = (THcDCTrack*)( fDCTracks->At(itrack)); - printf("%-5d ",itrack+1); - for (Int_t ihit=0;ihit<tempTrack->GetNHits();ihit++) { + if (fdebuglinkstubs) { + cout << " Number of tracks from link stubs = " << fNDCTracks << endl; + printf("%s %s \n","Track","Plane Wire "); + for (UInt_t itrack=0;itrack<fNDCTracks;itrack++) { + THcDCTrack *tempTrack = (THcDCTrack*)( fDCTracks->At(itrack)); + printf("%-5d ",itrack+1); + for (Int_t ihit=0;ihit<tempTrack->GetNHits();ihit++) { THcDCHit* hit=(THcDCHit*)(tempTrack->GetHit(ihit)); printf("%3d %3d",hit->GetPlaneNum(),hit->GetWireNum()); - } - printf("\n"); + } + printf("\n"); } } } @@ -943,21 +943,21 @@ void THcDC::TrackFit() THcDCHit* hit=theDCTrack->GetHit(ihit); planes[ihit]=hit->GetPlaneNum()-1; if(fFixLR) { - if(fFixPropagationCorrection) { - coords[ihit] = hit->GetPos() - + theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit); - } else { - coords[ihit] = hit->GetPos() - + theDCTrack->GetHitLR(ihit)*hit->GetDist(); - } + if(fFixPropagationCorrection) { + coords[ihit] = hit->GetPos() + + theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit); + } else { + coords[ihit] = hit->GetPos() + + theDCTrack->GetHitLR(ihit)*hit->GetDist(); + } } else { - if(fFixPropagationCorrection) { - coords[ihit] = hit->GetPos() - + hit->GetLR()*theDCTrack->GetHitDist(ihit); - } else { - coords[ihit] = hit->GetCoord(); - } + if(fFixPropagationCorrection) { + coords[ihit] = hit->GetPos() + + hit->GetLR()*theDCTrack->GetHitDist(ihit); + } else { + coords[ihit] = hit->GetCoord(); + } } @@ -969,39 +969,39 @@ void THcDC::TrackFit() TVectorD TT(NUM_FPRAY); TMatrixD AA(NUM_FPRAY,NUM_FPRAY); for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { - TT[irayp] = 0.0; - for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - - THcDCHit* hit=theDCTrack->GetHit(ihit); - - TT[irayp] += (coords[ihit]*fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])/pow(hit->GetWireSigma(),2); - if (hit->GetPlaneNum()==5) - { - // cout << "Plane: " << hit->GetPlaneNum() << endl; - //cout << "Wire: " <<hit->GetWireNum() << endl; - //cout << "Sigma: " << hit->GetWireSigma() << endl; - } + TT[irayp] = 0.0; + for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - } //end hit loop + THcDCHit* hit=theDCTrack->GetHit(ihit); + + TT[irayp] += (coords[ihit]*fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])/pow(hit->GetWireSigma(),2); + // if (hit->GetPlaneNum()==5) + // { + // // cout << "Plane: " << hit->GetPlaneNum() << endl; + // //cout << "Wire: " <<hit->GetWireNum() << endl; + // //cout << "Sigma: " << hit->GetWireSigma() << endl; + // } + + } //end hit loop } for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { - for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) { - AA[irayp][jrayp] = 0.0; - if(jrayp<irayp) { // Symmetric - AA[irayp][jrayp] = AA[jrayp][irayp]; - } else { - for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - - THcDCHit* hit=theDCTrack->GetHit(ihit); + for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) { + AA[irayp][jrayp] = 0.0; + if(jrayp<irayp) { // Symmetric + AA[irayp][jrayp] = AA[jrayp][irayp]; + } else { + for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { + + THcDCHit* hit=theDCTrack->GetHit(ihit); - AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]* - fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/ - pow(hit->GetWireSigma(),2); + AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]* + fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/ + pow(hit->GetWireSigma(),2); - } //end ihit loop - } - } + } //end ihit loop + } + } } // Solve 4x4 equations @@ -1018,12 +1018,12 @@ void THcDC::TrackFit() // Make sure fCoords, fResiduals, and fDoubleResiduals are clear for(Int_t iplane=0;iplane < fNPlanes; iplane++) { - Double_t coord=0.0; - for(Int_t ir=0;ir<NUM_FPRAY;ir++) { - coord += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir]; - // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl; - } - theDCTrack->SetCoord(iplane,coord); + Double_t coord=0.0; + for(Int_t ir=0;ir<NUM_FPRAY;ir++) { + coord += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir]; + // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl; + } + theDCTrack->SetCoord(iplane,coord); } // Compute Chi2 and residuals chi2 = 0.0; @@ -1045,56 +1045,56 @@ void THcDC::TrackFit() } theDCTrack->SetChisq(chi2); - // calculate ray without a plane in track + // calculate ray without a plane in track for(Int_t ipl_hit=0;ipl_hit < theDCTrack->GetNHits();ipl_hit++) { - if(theDCTrack->GetNFree() > 0) { - TVectorD TT(NUM_FPRAY); - TMatrixD AA(NUM_FPRAY,NUM_FPRAY); - for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { - TT[irayp] = 0.0; - for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { + if(theDCTrack->GetNFree() > 0) { + TVectorD TT(NUM_FPRAY); + TMatrixD AA(NUM_FPRAY,NUM_FPRAY); + for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { + TT[irayp] = 0.0; + for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - THcDCHit* hit=theDCTrack->GetHit(ihit); + THcDCHit* hit=theDCTrack->GetHit(ihit); - if (ihit != ipl_hit) { - TT[irayp] += (coords[ihit]* - fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]) - /pow(hit->GetWireSigma(),2); + if (ihit != ipl_hit) { + TT[irayp] += (coords[ihit]* + fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]) + /pow(hit->GetWireSigma(),2); - } + } + } } - } - for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { - for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) { - AA[irayp][jrayp] = 0.0; - if(jrayp<irayp) { // Symmetric - AA[irayp][jrayp] = AA[jrayp][irayp]; - } else { + for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) { + for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) { + AA[irayp][jrayp] = 0.0; + if(jrayp<irayp) { // Symmetric + AA[irayp][jrayp] = AA[jrayp][irayp]; + } else { - for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { + for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - THcDCHit* hit=theDCTrack->GetHit(ihit); + THcDCHit* hit=theDCTrack->GetHit(ihit); - if (ihit != ipl_hit) { - AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]* - fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/ - pow(hit->GetWireSigma(),2); + if (ihit != ipl_hit) { + AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]* + fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/ + pow(hit->GetWireSigma(),2); + } } } } } - } - // - // Solve 4x4 equations - // Should check that it is invertable - TVectorD dray(NUM_FPRAY); - AA.Invert(); - dray = AA*TT; + // + // Solve 4x4 equations + // Should check that it is invertable + TVectorD dray(NUM_FPRAY); + AA.Invert(); + dray = AA*TT; Double_t coord=0.0; for(Int_t ir=0;ir<NUM_FPRAY;ir++) { coord += fPlaneCoeffs[planes[ipl_hit]][raycoeffmap[ir]]*dray[ir]; @@ -1102,10 +1102,10 @@ void THcDC::TrackFit() } Double_t residual = coords[ipl_hit] - coord; theDCTrack->SetResidualExclPlane(planes[ipl_hit], residual); - } + } } } - //Calculate residual without plane + //Calculate residual without plane // Calculate residuals for each chamber if in single stub mode // and there was a track found in each chamber @@ -1189,41 +1189,41 @@ void THcDC::TrackFit() } } // - if (fdebugtrackprint) { - printf("%5s %-14s %-14s %-14s %-14s %-10s %-10s \n","Track","x_t","y_t","xp_t","yp_t","chi2","DOF"); - printf("%5s %-14s %-14s %-14s %-14s %-10s %-10s \n"," ","[cm]","[cm]","[rad]","[rad]"," "," "); - for(UInt_t itr=0;itr < fNDCTracks;itr++) { - THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itr)); - printf("%-5d %14.6e %14.6e %14.6e %14.6e %10.3e %3d \n", itr+1,theDCTrack->GetX(),theDCTrack->GetY(),theDCTrack->GetXP(),theDCTrack->GetYP(),theDCTrack->GetChisq(),theDCTrack->GetNFree()); - } - for(UInt_t itr=0;itr < fNDCTracks;itr++) { - printf("%s %5d \n","Hit info for track number = ",itr+1); - printf("%5s %-15s %-15s %-15s \n","Plane","WIRE_COORD","Fit postiion","Residual"); - THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itr)); - for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { - THcDCHit* hit=theDCTrack->GetHit(ihit); - Int_t plane=hit->GetPlaneNum()-1; - Double_t coords_temp; - if(fFixLR) { - if(fFixPropagationCorrection) { - coords_temp = hit->GetPos() - + theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit); - } else { - coords_temp = hit->GetPos() - + theDCTrack->GetHitLR(ihit)*hit->GetDist(); - } - } else { - if(fFixPropagationCorrection) { - coords_temp = hit->GetPos() - + hit->GetLR()*theDCTrack->GetHitDist(ihit); + if (fdebugtrackprint) { + printf("%5s %-14s %-14s %-14s %-14s %-10s %-10s \n","Track","x_t","y_t","xp_t","yp_t","chi2","DOF"); + printf("%5s %-14s %-14s %-14s %-14s %-10s %-10s \n"," ","[cm]","[cm]","[rad]","[rad]"," "," "); + for(UInt_t itr=0;itr < fNDCTracks;itr++) { + THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itr)); + printf("%-5d %14.6e %14.6e %14.6e %14.6e %10.3e %3d \n", itr+1,theDCTrack->GetX(),theDCTrack->GetY(),theDCTrack->GetXP(),theDCTrack->GetYP(),theDCTrack->GetChisq(),theDCTrack->GetNFree()); + } + for(UInt_t itr=0;itr < fNDCTracks;itr++) { + printf("%s %5d \n","Hit info for track number = ",itr+1); + printf("%5s %-15s %-15s %-15s \n","Plane","WIRE_COORD","Fit postiion","Residual"); + THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itr)); + for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) { + THcDCHit* hit=theDCTrack->GetHit(ihit); + Int_t plane=hit->GetPlaneNum()-1; + Double_t coords_temp; + if(fFixLR) { + if(fFixPropagationCorrection) { + coords_temp = hit->GetPos() + + theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit); + } else { + coords_temp = hit->GetPos() + + theDCTrack->GetHitLR(ihit)*hit->GetDist(); + } } else { - coords_temp = hit->GetCoord(); + if(fFixPropagationCorrection) { + coords_temp = hit->GetPos() + + hit->GetLR()*theDCTrack->GetHitDist(ihit); + } else { + coords_temp = hit->GetCoord(); + } } + printf("%-5d %15.7e %15.7e %15.7e \n",plane+1,coords_temp,theDCTrack->GetCoord(plane),theDCTrack->GetResidual(plane)); } - printf("%-5d %15.7e %15.7e %15.7e \n",plane+1,coords_temp,theDCTrack->GetCoord(plane),theDCTrack->GetResidual(plane)); - } - } - } + } + } // } diff --git a/src/THcDCHit.cxx b/src/THcDCHit.cxx index c28c146a5ef0094d5cd2e9af6e13799e34018df6..b9e0bad67b6121cba2c4cf975e30be2f491c9221 100644 --- a/src/THcDCHit.cxx +++ b/src/THcDCHit.cxx @@ -27,7 +27,7 @@ void THcDCHit::Print( Option_t* opt ) const << " time=" << GetTime() << " drift=" << GetDist(); // << " res=" << GetResolution() - // << " z=" << GetZ() + // << " z=" << GetZ() if( *opt != 'C' ) cout << endl; } diff --git a/src/THcDCHit.h b/src/THcDCHit.h index 25536be1676b6c53c5bde0608ab9e848aff458c8..6636e9e4f1a465acb780b4249c627e1760644930 100644 --- a/src/THcDCHit.h +++ b/src/THcDCHit.h @@ -59,110 +59,87 @@ public: void SetCorrectedStatus(Int_t c) { fCorrected = c; } Int_t GetReadoutSide() { - Int_t pln = fWirePlane->GetPlaneNum(); - Int_t tb = fWirePlane->GetReadoutTB(); - Int_t wn = fWire->GetNum(); - Int_t version = fWirePlane->GetVersion(); - - //if new HMS - if (version == 1){ - if ((pln>=3 && pln<=4) || (pln>=9 && pln<=10)){ - if (tb>0){ - if (wn < 60){ - fReadoutSide = 2; - } - else { - fReadoutSide = 4; - } - } - else { - if (wn < 44){ - fReadoutSide = 4; - } - else { - fReadoutSide = 2; - } - } - } - else{ - if (tb>0){ - if (wn < 51){ - fReadoutSide = 2; - } - else if (wn >= 51 && wn <= 64){ - fReadoutSide = 1; - } - else { - fReadoutSide =4; - } - } - else { - if (wn < 33){ - fReadoutSide = 4; - } - else if (wn >=33 && wn<=46){ - fReadoutSide = 1; - } - else { - fReadoutSide = 2; - } - } - } + Int_t pln = fWirePlane->GetPlaneNum(); + Int_t tb = fWirePlane->GetReadoutTB(); + Int_t wn = fWire->GetNum(); + Int_t version = fWirePlane->GetVersion(); + + //if new HMS + if (version == 1) { + if ((pln>=3 && pln<=4) || (pln>=9 && pln<=10)) { + if (tb>0) { + if (wn < 60) { + fReadoutSide = 2; + } else { + fReadoutSide = 4; } - - else{//appplies SHMS DC configuration - //check if x board - if ((pln>=3 && pln<=4) || (pln>=9 && pln<=10)){ - if (tb>0){ - if (wn < 49){ - fReadoutSide = 4; - } - else { - fReadoutSide = 2; - } - } - else { - if (wn < 33){ - fReadoutSide = 2; - } - else { - fReadoutSide = 4; - } - } - } - //else is u board - else{ - if (tb>0){ - if (wn < 41){ - fReadoutSide = 4; - } - else if (wn >= 41 && wn <= 63){ - fReadoutSide = 3; - } - else if (wn >=64 && wn <=69){ - fReadoutSide = 1; - } - else { - fReadoutSide = 2; - } - } - else { - if (wn < 39){ - fReadoutSide = 2; - } - else if (wn >=39 && wn<=44){ - fReadoutSide = 1; - } - else if (wn>=45 && wn<=67){ - fReadoutSide = 3; - } - else { - fReadoutSide = 4; - } - } - } + } else { + if (wn < 44) { + fReadoutSide = 4; + } else { + fReadoutSide = 2; + } + } + } else { + if (tb>0) { + if (wn < 51) { + fReadoutSide = 2; + } else if (wn >= 51 && wn <= 64) { + fReadoutSide = 1; + } else { + fReadoutSide =4; + } + } else { + if (wn < 33) { + fReadoutSide = 4; + } else if (wn >=33 && wn<=46) { + fReadoutSide = 1; + } else { + fReadoutSide = 2; + } + } + } + } else{//appplies SHMS DC configuration + //check if x board + if ((pln>=3 && pln<=4) || (pln>=9 && pln<=10)) { + if (tb>0) { + if (wn < 49) { + fReadoutSide = 4; + } else { + fReadoutSide = 2; } - return fReadoutSide; + } else { + if (wn < 33) { + fReadoutSide = 2; + } else { + fReadoutSide = 4; + } + } + } else { //else is u board + if (tb>0) { + if (wn < 41) { + fReadoutSide = 4; + } else if (wn >= 41 && wn <= 63) { + fReadoutSide = 3; + } else if (wn >=64 && wn <=69) { + fReadoutSide = 1; + } else { + fReadoutSide = 2; + } + } else { + if (wn < 39) { + fReadoutSide = 2; + } else if (wn >=39 && wn<=44) { + fReadoutSide = 1; + } else if (wn>=45 && wn<=67) { + fReadoutSide = 3; + } else { + fReadoutSide = 4; + } + } + } + } + return fReadoutSide; } diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx index 18a880b3b7c4f48e750a2b297f3a2ed767c13cca..a47dea84f6143fb494246dc86190543c1d4167cc 100644 --- a/src/THcDriftChamber.cxx +++ b/src/THcDriftChamber.cxx @@ -94,34 +94,34 @@ void THcDriftChamber::AddPlane(THcDriftChamberPlane *plane) { plane->SetPlaneIndex(fNPlanes); fPlanes.push_back(plane); - // HMS Specific + // HMS Specific // Hard code Y plane numbers. Eventually need to get from database if (fHMSStyleChambers) { - if(fChamberNum == 1) { - YPlaneNum=2; - YPlanePNum=5; - if(plane->GetPlaneNum() == 2) YPlaneInd = fNPlanes; - if(plane->GetPlaneNum() == 5) YPlanePInd = fNPlanes; - } else { - YPlaneNum=8; - YPlanePNum=11; - if(plane->GetPlaneNum() == 8) YPlaneInd = fNPlanes; - if(plane->GetPlaneNum() == 11) YPlanePInd = fNPlanes; - } + if(fChamberNum == 1) { + YPlaneNum=2; + YPlanePNum=5; + if(plane->GetPlaneNum() == 2) YPlaneInd = fNPlanes; + if(plane->GetPlaneNum() == 5) YPlanePInd = fNPlanes; + } else { + YPlaneNum=8; + YPlanePNum=11; + if(plane->GetPlaneNum() == 8) YPlaneInd = fNPlanes; + if(plane->GetPlaneNum() == 11) YPlanePInd = fNPlanes; + } } else { - // SOS Specific - // Hard code X plane numbers. Eventually need to get from database - if(fChamberNum == 1) { - XPlaneNum=3; - XPlanePNum=4; - if(plane->GetPlaneNum() == 3) XPlaneInd = fNPlanes; - if(plane->GetPlaneNum() == 4) XPlanePInd = fNPlanes; - } else { - XPlaneNum=9; - XPlanePNum=10; - if(plane->GetPlaneNum() == 9) XPlaneInd = fNPlanes; - if(plane->GetPlaneNum() == 10) XPlanePInd = fNPlanes; - } + // SOS Specific + // Hard code X plane numbers. Eventually need to get from database + if(fChamberNum == 1) { + XPlaneNum=3; + XPlanePNum=4; + if(plane->GetPlaneNum() == 3) XPlaneInd = fNPlanes; + if(plane->GetPlaneNum() == 4) XPlanePInd = fNPlanes; + } else { + XPlaneNum=9; + XPlanePNum=10; + if(plane->GetPlaneNum() == 9) XPlaneInd = fNPlanes; + if(plane->GetPlaneNum() == 10) XPlanePInd = fNPlanes; + } } fNPlanes++; return; @@ -165,7 +165,7 @@ Int_t THcDriftChamber::ReadDatabase( const TDatime& date ) fSpacePointCriterion = static_cast<THcDC*>(fParent)->GetSpacePointCriterion(fChamberNum); fMaxDist = TMath::Sqrt(fSpacePointCriterion/2.0); // For easy space points - if (fhdebugflagpr) cout << " cham = " << fChamberNum << " Set yplane num " << YPlaneNum << " "<< YPlanePNum << endl; + if (fhdebugflagpr) cout << " cham = " << fChamberNum << " Set yplane num " << YPlaneNum << " "<< YPlanePNum << endl; // Generate the HAA3INV matrix for all the acceptable combinations // of hit planes. Try to make it as generic as possible // pindex=0 -> Plane 1 missing, pindex5 -> plane 6 missing. Won't @@ -229,19 +229,19 @@ Int_t THcDriftChamber::DefineVariables( EMode mode ) fIsSetup = ( mode == kDefine ); // Register variables in global list - RVarDef vars[] = { - { "maxhits", "Maximum hits allowed", "fMaxHits" }, - { "spacepoints", "Space points of DC", "fNSpacePoints" }, - { "nhit", "Number of DC hits", "fNhits" }, - { "trawhit", "Number of True Raw hits", "fN_True_RawHits" }, - { "stub_x", "", "fSpacePoints.THcSpacePoint.GetStubX()" }, - { "stub_xp", "", "fSpacePoints.THcSpacePoint.GetStubXP()" }, - { "stub_y", "", "fSpacePoints.THcSpacePoint.GetStubY()" }, - { "stub_yp", "", "fSpacePoints.THcSpacePoint.GetStubYP()" }, - { "ncombos", "", "fSpacePoints.THcSpacePoint.GetCombos()" }, - { 0 } - }; - return DefineVarsFromList( vars, mode ); + RVarDef vars[] = { + { "maxhits", "Maximum hits allowed", "fMaxHits" }, + { "spacepoints", "Space points of DC", "fNSpacePoints" }, + { "nhit", "Number of DC hits", "fNhits" }, + { "trawhit", "Number of True Raw hits", "fN_True_RawHits" }, + { "stub_x", "", "fSpacePoints.THcSpacePoint.GetStubX()" }, + { "stub_xp", "", "fSpacePoints.THcSpacePoint.GetStubXP()" }, + { "stub_y", "", "fSpacePoints.THcSpacePoint.GetStubY()" }, + { "stub_yp", "", "fSpacePoints.THcSpacePoint.GetStubYP()" }, + { "ncombos", "", "fSpacePoints.THcSpacePoint.GetCombos()" }, + { 0 } + }; + return DefineVarsFromList( vars, mode ); //return kOK; } @@ -267,10 +267,10 @@ void THcDriftChamber::PrintDecode( void ) { cout << " Num of nits = " << fNhits << endl; cout << " Num " << " Plane " << " Wire " << " Wire-Center " << " RAW TDC " << " Drift time" << endl; - for(Int_t ihit=0;ihit<fNhits;ihit++) { + for(Int_t ihit=0;ihit<fNhits;ihit++) { THcDCHit* thishit = fHits[ihit]; cout << ihit << " " <<thishit->GetPlaneNum() << " " << thishit->GetWireNum() << " " << thishit->GetPos() << " " << thishit->GetRawTime() << " " << thishit->GetTime() << endl; - } + } } @@ -357,7 +357,7 @@ Int_t THcDriftChamber::FindSpacePoints( void ) if(fRemove_Sppt_If_One_YPlane == 1) { // The routine is specific to HMS //Int_t ndest= - DestroyPoorSpacePoints(); // Only for HMS? + DestroyPoorSpacePoints(); // Only for HMS? // Loop over space points and remove those with less than 4 planes // hit and missing hits in Y,Y' planes } @@ -681,17 +681,17 @@ Int_t THcDriftChamber::DestroyPoorSpacePoints() // Does this work, or do we have to copy each member? // If it doesn't we should overload the = operator //(*fSpacePoints)[isp] = (*fSpacePoints)[osp]; - THcSpacePoint* spi = (THcSpacePoint*)(*fSpacePoints)[isp]; - THcSpacePoint* spo = (THcSpacePoint*)(*fSpacePoints)[osp]; - spi->Clear(); - Double_t xt,yt; - xt=spo->GetX(); - yt=spo->GetY(); - spi->SetXY(xt, yt); - for(Int_t ihit=0;ihit<spo->GetNHits();ihit++) { - THcDCHit* hit = spo->GetHit(ihit); - spi->AddHit(hit); - } + THcSpacePoint* spi = (THcSpacePoint*)(*fSpacePoints)[isp]; + THcSpacePoint* spo = (THcSpacePoint*)(*fSpacePoints)[osp]; + spi->Clear(); + Double_t xt,yt; + xt=spo->GetX(); + yt=spo->GetY(); + spi->SetXY(xt, yt); + for(Int_t ihit=0;ihit<spo->GetNHits();ihit++) { + THcDCHit* hit = spo->GetHit(ihit); + spi->AddHit(hit); + } } } return nremoved; @@ -732,9 +732,9 @@ Int_t THcDriftChamber::SpacePointMultiWire() for(Int_t ip=0;ip<fNPlanes;ip++) { nhitsperplane[ip] = 0; - for(Int_t ih=0;ih<MAX_HITS_PER_POINT;ih++) { - hits_plane[ip][ih] = 0; - } + for(Int_t ih=0;ih<MAX_HITS_PER_POINT;ih++) { + hits_plane[ip][ih] = 0; + } } // Sort Space Points hits by plane THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp]; @@ -913,17 +913,17 @@ void THcDriftChamber::SelectSpacePoints() if(sp->GetNHits() >= fMinHits) { if(isp > sp_count) { // (*fSpacePoints)[sp_count] = (*fSpacePoints)[isp]; - THcSpacePoint* sp1 = (THcSpacePoint*)(*fSpacePoints)[sp_count]; - sp1->Clear(); - Double_t xt,yt; - xt=sp->GetX(); - yt=sp->GetY(); - sp1->SetXY(xt, yt); - sp1->SetCombos(sp->GetCombos()); - for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) { + THcSpacePoint* sp1 = (THcSpacePoint*)(*fSpacePoints)[sp_count]; + sp1->Clear(); + Double_t xt,yt; + xt=sp->GetX(); + yt=sp->GetY(); + sp1->SetXY(xt, yt); + sp1->SetCombos(sp->GetCombos()); + for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) { THcDCHit* hit = sp->GetHit(ihit); - sp1->AddHit(hit); - } + sp1->AddHit(hit); + } } sp_count++; } @@ -958,119 +958,117 @@ void THcDriftChamber::CorrectHitTimes() //if (fhdebugflagpr) cout << "In correcthittimes fNSpacePoints = " << fNSpacePoints << endl; for(Int_t isp=0;isp<fNSpacePoints;isp++) { - THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp]; - Double_t x = sp->GetX(); - Double_t y = sp->GetY(); - for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) { - THcDCHit* hit = sp->GetHit(ihit); - THcDriftChamberPlane* plane=hit->GetWirePlane(); - - // How do we know this correction only gets applied once? Is - // it determined that a given hit can only belong to one space point? - Double_t time_corr = plane->GetReadoutX() ? - y*plane->GetReadoutCorr()/fWireVelocity : - x*plane->GetReadoutCorr()/fWireVelocity; - - // This applies the wire velocity correction for new SHMS chambers --hszumila, SEP17 - if (!fHMSStyleChambers){ - Int_t pln = hit->GetPlaneNum(); - Int_t readoutSide = hit->GetReadoutSide(); - - Double_t posn = hit->GetPos(); - //The following values are determined from param file as permutations on planes 5 and 10 - Int_t readhoriz = plane->GetReadoutLR(); - Int_t readvert = plane->GetReadoutTB(); - - //+x is up and +y is beam right! - double alpha = static_cast<THcDC*>(fParent)->GetAlphaAngle(pln); - double xc = posn*TMath::Sin(alpha); - double yc = posn*TMath::Cos(alpha); - - Double_t wireDistance = plane->GetReadoutX() ? - (abs(y-yc))*abs(plane->GetReadoutCorr()) : - (abs(x-xc))*abs(plane->GetReadoutCorr()); - - //Readout side is based off wiring diagrams - switch (readoutSide){ - case 1: //readout from top of chamber - if (x>xc){wireDistance = -readvert*wireDistance;} - else{wireDistance = readvert*wireDistance;} - - break; - 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 - } - - } - ///////////////////////////////////////////////////////////// - - // 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. - else{ - 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 - } - } - } + THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp]; + Double_t x = sp->GetX(); + Double_t y = sp->GetY(); + for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) { + THcDCHit* hit = sp->GetHit(ihit); + THcDriftChamberPlane* plane=hit->GetWirePlane(); + + // How do we know this correction only gets applied once? Is + // it determined that a given hit can only belong to one space point? + Double_t time_corr = plane->GetReadoutX() ? + y*plane->GetReadoutCorr()/fWireVelocity : + x*plane->GetReadoutCorr()/fWireVelocity; + + // This applies the wire velocity correction for new SHMS chambers --hszumila, SEP17 + if (!fHMSStyleChambers){ + Int_t pln = hit->GetPlaneNum(); + Int_t readoutSide = hit->GetReadoutSide(); + + Double_t posn = hit->GetPos(); + //The following values are determined from param file as permutations on planes 5 and 10 + Int_t readhoriz = plane->GetReadoutLR(); + Int_t readvert = plane->GetReadoutTB(); + + //+x is up and +y is beam right! + double alpha = static_cast<THcDC*>(fParent)->GetAlphaAngle(pln); + double xc = posn*TMath::Sin(alpha); + double yc = posn*TMath::Cos(alpha); + + Double_t wireDistance = plane->GetReadoutX() ? + (abs(y-yc))*abs(plane->GetReadoutCorr()) : + (abs(x-xc))*abs(plane->GetReadoutCorr()); + + //Readout side is based off wiring diagrams + switch (readoutSide){ + case 1: //readout from top of chamber + if (x>xc){wireDistance = -readvert*wireDistance;} + else{wireDistance = readvert*wireDistance;} + + break; + 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 + } + } + } } } UInt_t THcDriftChamber::Count1Bits(UInt_t x) @@ -1197,7 +1195,7 @@ void THcDriftChamber::LeftRight() } if (nplaneshit >= fNPlanes-1) { Double_t chi2; - chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub); + chi2 = FindStub(nhits, sp,plane_list, bitpat, plusminus, stub); if (fdebugstubchisq) cout << " pmloop = " << pmloop << " chi2 = " << chi2 << endl; if(chi2 < minchi2) { if (fStubMaxXPDiff<100. ) { @@ -1220,13 +1218,13 @@ void THcDriftChamber::LeftRight() sp->SetStub(stub); } else { // Record best stub failing angle cut if (chi2 < tmp_minchi2) { - 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]; - } + 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 @@ -1256,52 +1254,51 @@ void THcDriftChamber::LeftRight() } sp->SetStub(stub); } - }////////// - else { + } else { 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]; + } 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); } - sp->SetStub(tmp_stub); - } - 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; + 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 diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx index bc09d4fb699644834f4798f60093fecd5a743fd0..4acfbacbdbb888bdbc6805ff4cc7373d3da86c24 100644 --- a/src/THcDriftChamberPlane.cxx +++ b/src/THcDriftChamberPlane.cxx @@ -262,7 +262,7 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date ) THaApparatus* app = GetApparatus(); const char* nm = "hod"; if( !app || - !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector(nm))) ) { + !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector(nm))) ) { static const char* const here = "ReadDatabase()"; Warning(Here(here),"Hodoscope \"%s\" not found. " "Event-by-event time offsets will NOT be used!!",nm); @@ -360,7 +360,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) fNRawhits++; /* Sort into early, late and ontime */ Int_t rawnorefcorrtdc = hit->GetRawTdcHit().GetTimeRaw(mhit); // Get the ref time subtracted time - Int_t rawtdc = hit->GetRawTdcHit().GetTime(mhit); // Get the ref time subtracted time + Int_t rawtdc = hit->GetRawTdcHit().GetTime(mhit); // Get the ref time subtracted time if(rawtdc < fTdcWinMin) { // Increment early counter (Actually late because TDC is backward) } else if (rawtdc > fTdcWinMax) { @@ -380,11 +380,11 @@ Int_t THcDriftChamberPlane::SubtractStartTime() { Double_t StartTime = 0.0; if( fglHod ) StartTime = fglHod->GetStartTime(); - for(Int_t ihit=0;ihit<GetNHits();ihit++) { - THcDCHit *thishit = (THcDCHit*) fHits->At(ihit); - Double_t temptime= thishit->GetTime()-StartTime; - thishit->SetTime(temptime); - thishit->ConvertTimeToDist(); - } + for(Int_t ihit=0;ihit<GetNHits();ihit++) { + THcDCHit *thishit = (THcDCHit*) fHits->At(ihit); + Double_t temptime= thishit->GetTime()-StartTime; + thishit->SetTime(temptime); + thishit->ConvertTimeToDist(); + } return 0; } diff --git a/src/THcHallCSpectrometer.cxx b/src/THcHallCSpectrometer.cxx index b499688fb1d9f7f162f3190bc5b60e529eeb7d3e..e1a83478bedca920f42a135ed5948e568cb2c9bb 100644 --- a/src/THcHallCSpectrometer.cxx +++ b/src/THcHallCSpectrometer.cxx @@ -288,23 +288,23 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date ) // mispointing in transport system y is horizontal and +x is vertical down if (fMispointing_y == 999.) { if (prefix[0]=='h') { - fMispointing_y = 0.1*(0.52-0.012*40.+0.002*40.*40.); - if (fTheta_lab < 40) fMispointing_y = 0.1*(0.52-0.012*TMath::Abs(fTheta_lab)+0.002*fTheta_lab*fTheta_lab); + fMispointing_y = 0.1*(0.52-0.012*40.+0.002*40.*40.); + if (fTheta_lab < 40) fMispointing_y = 0.1*(0.52-0.012*TMath::Abs(fTheta_lab)+0.002*fTheta_lab*fTheta_lab); } if (prefix[0]=='p') fMispointing_y = 0.1*(-0.6); - cout << prefix[0] << " From Formula Mispointing_y = " << fMispointing_y << endl; + cout << prefix[0] << " From Formula Mispointing_y = " << fMispointing_y << endl; } else { - cout << prefix[0] << " From Parameter Set Mispointing_y = " << fMispointing_y << endl; + cout << prefix[0] << " From Parameter Set Mispointing_y = " << fMispointing_y << endl; } if (fMispointing_x == 999.) { if (prefix[0]=='h') { - fMispointing_x = 0.1*(2.37-0.086*50+0.0012*50.*50.); - if (fTheta_lab < 50)fMispointing_x = 0.1*(2.37-0.086*TMath::Abs(fTheta_lab)+0.0012*fTheta_lab*fTheta_lab); + fMispointing_x = 0.1*(2.37-0.086*50+0.0012*50.*50.); + if (fTheta_lab < 50)fMispointing_x = 0.1*(2.37-0.086*TMath::Abs(fTheta_lab)+0.0012*fTheta_lab*fTheta_lab); } if (prefix[0]=='p') fMispointing_x = 0.1*(-1.26); - cout << prefix[0] << " From Formula Mispointing_x = " << fMispointing_x << endl; + cout << prefix[0] << " From Formula Mispointing_x = " << fMispointing_x << endl; } else { - cout << prefix[0] << " From Parameter Set Mispointing_x = " << fMispointing_x << endl; + cout << prefix[0] << " From Parameter Set Mispointing_x = " << fMispointing_x << endl; } // @@ -385,7 +385,7 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date ) ,&fReconTerms[fNReconTerms].Exp[2] ,&fReconTerms[fNReconTerms].Exp[3] ,&fReconTerms[fNReconTerms].Exp[4]); - fNReconTerms++; + fNReconTerms++; good = getline(ifile,line).good(); } cout << "Read " << fNReconTerms << " matrix element terms" << endl; @@ -433,7 +433,7 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks ) Double_t xptar=kBig,yptar=kBig,ytar=kBig,delta=kBig; Double_t xtar=0; CalculateTargetQuantities(track,xtar,xptar,ytar,yptar,delta); - // Transfer results to track + // Transfer results to track // No beam raster yet //; In transport coordinates phi = hyptar = dy/dz and theta = hxptar = dx/dz //; but for unknown reasons the yp offset is named htheta_offset @@ -442,7 +442,7 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks ) track->SetTarget(0.0, ytar*100.0, xptar, yptar); track->SetDp(delta*100.0); // Percent. Double_t ptemp = fPcentral*(1+track->GetDp()/100.0); - track->SetMomentum(ptemp); + track->SetMomentum(ptemp); TVector3 pvect_temp; TransportToLab(track->GetP(),track->GetTTheta(),track->GetTPhi(),pvect_temp); track->SetPvect(pvect_temp); @@ -528,9 +528,9 @@ Int_t THcHallCSpectrometer::TrackCalc() fGoldenTrack = static_cast<THaTrack*>( fTracks->At(hit_gold_track) ); fTrkIfo = *fGoldenTrack; fTrk = fGoldenTrack; - } else + } else { fGoldenTrack = NULL; - + } return TrackTimes( fTracks ); } @@ -546,11 +546,11 @@ Int_t THcHallCSpectrometer::BestTrackSimple() if( GetTrSorting() ) fTracks->Sort(); // Assign index=0 to the best track, - for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){ - THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) ); - aTrack->SetIndex(1); - if (itrack==0) aTrack->SetIndex(0); - } + for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){ + THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) ); + aTrack->SetIndex(1); + if (itrack==0) aTrack->SetIndex(0); + } return(0); } @@ -607,75 +607,74 @@ Int_t THcHallCSpectrometer::BestTrackUsingScin() ( aTrack->GetBeta() > fSelBetaMin ) && ( aTrack->GetBeta() < fSelBetaMax ) && ( aTrack->GetEnergy() > fSelEtMin ) && - ( aTrack->GetEnergy() < fSelEtMax ) ) - { - Int_t x2D, y2D; - - if ( fNtracks > 1 ){ - Double_t hitpos3 = aTrack->GetX() + aTrack->GetTheta() * ( fScin2XZpos + 0.5 * fScin2XdZpos ); - Int_t icounter3 = TMath::Nint( ( hitpos3 - fHodo->GetPlaneCenter(2) ) / fHodo->GetPlaneSpacing(2) ) + 1; - Int_t hitCnt3 = TMath::Max( TMath::Min(icounter3, (Int_t) fHodo->GetNPaddles(2) ) , 1); // scin_2x_nr = 16 - // fHitDist3 = fHitPos3 - ( fHodo->GetPlaneSpacing(2) * ( hitCnt3 - 1 ) + fHodo->GetPlaneCenter(2) ); - Double_t hitpos4 = aTrack->GetY() + aTrack->GetPhi() * ( fScin2YZpos + 0.5 * fScin2YdZpos ); - Int_t icounter4 = TMath::Nint( ( fHodo->GetPlaneCenter(3) - hitpos4 ) / fHodo->GetPlaneSpacing(3) ) + 1; - Int_t hitCnt4 = TMath::Max( TMath::Min(icounter4, (Int_t) fHodo->GetNPaddles(3) ) , 1); // scin_2y_nr = 10 - // fHitDist4 = fHitPos4 - ( fHodo->GetPlaneCenter(3) - fHodo->GetPlaneSpacing(3) * ( hitCnt4 - 1 ) ); - // Plane 3 - Int_t mindiff=1000; - for (UInt_t i = 0; i < fHodo->GetNPaddles(2); i++ ){ - if ( x2Hits[i] ) { - Int_t diff = TMath::Abs((Int_t)hitCnt3-(Int_t)i-1); - if (diff < mindiff) mindiff = diff; - } - } - if(mindiff < 1000) { - x2D = mindiff; - } else { - x2D = 0; // Is this what we really want if there were no hits on this plane? + ( aTrack->GetEnergy() < fSelEtMax ) ) { + Int_t x2D, y2D; + + if ( fNtracks > 1 ){ + Double_t hitpos3 = aTrack->GetX() + aTrack->GetTheta() * ( fScin2XZpos + 0.5 * fScin2XdZpos ); + Int_t icounter3 = TMath::Nint( ( hitpos3 - fHodo->GetPlaneCenter(2) ) / fHodo->GetPlaneSpacing(2) ) + 1; + Int_t hitCnt3 = TMath::Max( TMath::Min(icounter3, (Int_t) fHodo->GetNPaddles(2) ) , 1); // scin_2x_nr = 16 + // fHitDist3 = fHitPos3 - ( fHodo->GetPlaneSpacing(2) * ( hitCnt3 - 1 ) + fHodo->GetPlaneCenter(2) ); + Double_t hitpos4 = aTrack->GetY() + aTrack->GetPhi() * ( fScin2YZpos + 0.5 * fScin2YdZpos ); + Int_t icounter4 = TMath::Nint( ( fHodo->GetPlaneCenter(3) - hitpos4 ) / fHodo->GetPlaneSpacing(3) ) + 1; + Int_t hitCnt4 = TMath::Max( TMath::Min(icounter4, (Int_t) fHodo->GetNPaddles(3) ) , 1); // scin_2y_nr = 10 + // fHitDist4 = fHitPos4 - ( fHodo->GetPlaneCenter(3) - fHodo->GetPlaneSpacing(3) * ( hitCnt4 - 1 ) ); + // Plane 3 + Int_t mindiff=1000; + for (UInt_t i = 0; i < fHodo->GetNPaddles(2); i++ ){ + if ( x2Hits[i] ) { + Int_t diff = TMath::Abs((Int_t)hitCnt3-(Int_t)i-1); + if (diff < mindiff) mindiff = diff; } + } + if(mindiff < 1000) { + x2D = mindiff; + } else { + x2D = 0; // Is this what we really want if there were no hits on this plane? + } - // Plane 4 - mindiff = 1000; - for (UInt_t i = 0; i < fHodo->GetNPaddles(3); i++ ){ - if ( y2Hits[i] ) { - Int_t diff = TMath::Abs((Int_t)hitCnt4-(Int_t)i-1); - if (diff < mindiff) mindiff = diff; - } - } - if(mindiff < 1000) { - y2D = mindiff; - } else { - y2D = 0; // Is this what we really want if there were no hits on this plane? + // Plane 4 + mindiff = 1000; + for (UInt_t i = 0; i < fHodo->GetNPaddles(3); i++ ){ + if ( y2Hits[i] ) { + Int_t diff = TMath::Abs((Int_t)hitCnt4-(Int_t)i-1); + if (diff < mindiff) mindiff = diff; } - } else { // Only a single track - x2D = 0.; - y2D = 0.; } + if(mindiff < 1000) { + y2D = mindiff; + } else { + y2D = 0; // Is this what we really want if there were no hits on this plane? + } + } else { // Only a single track + x2D = 0.; + y2D = 0.; + } + + if ( y2D <= y2Dmin ) { + if ( y2D < y2Dmin ) { + x2Dmin = 100; + chi2Min = 10000000000.; + } // y2D min - if ( y2D <= y2Dmin ) { - if ( y2D < y2Dmin ) { - x2Dmin = 100; - chi2Min = 10000000000.; - } // y2D min - - if ( x2D <= x2Dmin ){ - if ( x2D < x2Dmin ){ - chi2Min = 10000000000.0; - } // condition x2D - if ( chi2PerDeg < chi2Min ){ - - fGoodTrack = itrack; // fGoodTrack = itrack - y2Dmin = y2D; - x2Dmin = x2D; - chi2Min = chi2PerDeg; - - fGoldenTrack = static_cast<THaTrack*>( fTracks->At( fGoodTrack ) ); - fTrkIfo = *fGoldenTrack; - fTrk = fGoldenTrack; - } + if ( x2D <= x2Dmin ){ + if ( x2D < x2Dmin ){ + chi2Min = 10000000000.0; } // condition x2D - } // condition y2D - } // conditions for dedx, beta and trac energy + if ( chi2PerDeg < chi2Min ){ + + fGoodTrack = itrack; // fGoodTrack = itrack + y2Dmin = y2D; + x2Dmin = x2D; + chi2Min = chi2PerDeg; + + fGoldenTrack = static_cast<THaTrack*>( fTracks->At( fGoodTrack ) ); + fTrkIfo = *fGoldenTrack; + fTrk = fGoldenTrack; + } + } // condition x2D + } // condition y2D + } // conditions for dedx, beta and trac energy } // confition for fNFreeFP greater than fSelNDegreesMin } // loop over tracks @@ -705,7 +704,7 @@ Int_t THcHallCSpectrometer::BestTrackUsingScin() aTrack->SetIndex(1); if (iitrack==fGoodTrack) aTrack->SetIndex(0); } - // + // } } @@ -938,16 +937,16 @@ Int_t THcHallCSpectrometer::BestTrackUsingPrune() chi2Min = chi2PerDeg; } } - // Set index=0 for fGoodTrack - for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){ - THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) ); - aTrack->SetIndex(1); - if (iitrack==fGoodTrack) aTrack->SetIndex(0); - } - // + // Set index=0 for fGoodTrack + for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){ + THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) ); + aTrack->SetIndex(1); + if (iitrack==fGoodTrack) aTrack->SetIndex(0); + } + // - } + } return(0); } diff --git a/src/THcHodoEff.cxx b/src/THcHodoEff.cxx index c2617624fc42d6fd7f0a18a32ec981da2fd64e9a..9bb9d5fe40f8e41d6ac5d5e4de8ee2bb7e6fe756 100644 --- a/src/THcHodoEff.cxx +++ b/src/THcHodoEff.cxx @@ -82,34 +82,30 @@ Int_t THcHodoEff::Begin( THaRunBase* ) Int_t THcHodoEff::End( THaRunBase* ) { // End of analysis - for(Int_t ip=0;ip<fNPlanes;ip++) { - fStatAndEff[ip]=0; - for(Int_t ic=0;ic<fNCounters[ip];ic++) { - fStatTrkSum[ip]+=fStatTrk[fHod->GetScinIndex(ip,ic)]; - fStatAndSum[ip]+=fHodoAndEffi[fHod->GetScinIndex(ip,ic)]; - } - if (fStatTrkSum[ip] !=0) fStatAndEff[ip]=float(fStatAndSum[ip])/float(fStatTrkSum[ip]); - } + for(Int_t ip=0;ip<fNPlanes;ip++) { + fStatAndEff[ip]=0; + for(Int_t ic=0;ic<fNCounters[ip];ic++) { + fStatTrkSum[ip]+=fStatTrk[fHod->GetScinIndex(ip,ic)]; + fStatAndSum[ip]+=fHodoAndEffi[fHod->GetScinIndex(ip,ic)]; + } + if (fStatTrkSum[ip] !=0) fStatAndEff[ip]=float(fStatAndSum[ip])/float(fStatTrkSum[ip]); + } // - Double_t p1; - Double_t p2; - Double_t p3; - Double_t p4; - p1=fStatAndEff[0]; - p2=fStatAndEff[1]; - p3=fStatAndEff[2]; - p4=fStatAndEff[3]; -// probability that ONLY the listed planes had triggers - Double_t p1234= p1*p2*p3*p4; - Double_t p123 = p1*p2*p3*(1.-p4); - Double_t p124 = p1*p2*(1.-p3)*p4; - Double_t p134 = p1*(1.-p2)*p3*p4; - Double_t p234 = (1.-p1)*p2*p3*p4; - fHodoEff_s1 = 1. - ((1.-p1)*(1.-p2)); - fHodoEff_s2 = 1. - ((1.-p3)*(1.-p4)); - fHodoEff_tof=fHodoEff_s1 * fHodoEff_s2; - fHodoEff_3_of_4=p1234+p123+p124+p134+p234; - fHodoEff_4_of_4=p1234; + Double_t p1=fStatAndEff[0]; + Double_t p2=fStatAndEff[1]; + Double_t p3=fStatAndEff[2]; + Double_t p4=fStatAndEff[3]; + // probability that ONLY the listed planes had triggers + Double_t p1234= p1*p2*p3*p4; + Double_t p123 = p1*p2*p3*(1.-p4); + Double_t p124 = p1*p2*(1.-p3)*p4; + Double_t p134 = p1*(1.-p2)*p3*p4; + Double_t p234 = (1.-p1)*p2*p3*p4; + fHodoEff_s1 = 1. - ((1.-p1)*(1.-p2)); + fHodoEff_s2 = 1. - ((1.-p3)*(1.-p4)); + fHodoEff_tof=fHodoEff_s1 * fHodoEff_s2; + fHodoEff_3_of_4=p1234+p123+p124+p134+p234; + fHodoEff_4_of_4=p1234; return 0; } @@ -135,7 +131,7 @@ THaAnalysisObject::EStatus THcHodoEff::Init( const TDatime& run_time ) cout << "THcHodoEff::Init nplanes=" << fHod->GetNPlanes() << endl; cout << "THcHodoEff::Init Apparatus = " << fHod->GetName() << " " << - (fHod->GetApparatus())->GetName() << endl; + (fHod->GetApparatus())->GetName() << endl; return fStatus = kOK; } @@ -220,13 +216,13 @@ Int_t THcHodoEff::ReadDatabase( const TDatime& date ) for(Int_t ic=0;ic<fNCounters[ip];ic++) { fStatTrkDel[ip][ic].resize(20); // Max this settable fStatAndHitDel[ip][ic].resize(20); // Max this settable - + fHodoPosEffi[fHod->GetScinIndex(ip,ic)] = 0; fHodoNegEffi[fHod->GetScinIndex(ip,ic)] = 0; fHodoOrEffi[fHod->GetScinIndex(ip,ic)] = 0; fHodoAndEffi[fHod->GetScinIndex(ip,ic)] = 0; fStatTrk[fHod->GetScinIndex(ip,ic)] = 0; - } + } } // Int_t fHodPaddles = fNCounters[0]; @@ -303,17 +299,17 @@ Int_t THcHodoEff::Process( const THaEvData& evdata ) if(ip%2 == 0) { // X Plane hitPos[ip] = theTrack->GetX() + theTrack->GetTheta()*fPosZ[ip]; hitCounter[ip] = TMath::Max( - TMath::Min( - TMath::Nint((hitPos[ip]-fCenterFirst[ip])/ - fSpacing[ip]+1),fNCounters[ip] ),1); + TMath::Min( + TMath::Nint((hitPos[ip]-fCenterFirst[ip])/ + fSpacing[ip]+1),fNCounters[ip] ),1); hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) + fCenterFirst[ip]); } else { // Y Plane hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip]; hitCounter[ip] = TMath::Max( - TMath::Min( - TMath::Nint((fCenterFirst[ip]-hitPos[ip])/ - fSpacing[ip]+1), fNCounters[ip] ),1); + TMath::Min( + TMath::Nint((fCenterFirst[ip]-hitPos[ip])/ + fSpacing[ip]+1), fNCounters[ip] ),1); hitDistance[ip] = hitPos[ip] -(fCenterFirst[ip] - fSpacing[ip]*(hitCounter[ip]-1)); } @@ -365,7 +361,7 @@ Int_t THcHodoEff::Process( const THaEvData& evdata ) // lookat[ip] = TRUE; } fHitPlane[ip] = 0; - } + } // Is there a hit on or adjacent to paddle that track // passes through? diff --git a/src/THcReactionPoint.cxx b/src/THcReactionPoint.cxx index 65032aa86dca7a38ad301e4fbcdbc6a1a4ecd9f5..4c1aa19d091e2d0f89e64a2a9383faee4b05d782 100644 --- a/src/THcReactionPoint.cxx +++ b/src/THcReactionPoint.cxx @@ -75,9 +75,9 @@ Int_t THcReactionPoint::Process( const THaEvData& ) theTrack->SetVertex(v); if( theTrack == fSpectro->GetGoldenTrack() ) { fVertex = theTrack->GetVertex(); - fVertexOK = kTRUE; + fVertexOK = kTRUE; } - } + } return 0; }