From 95331d70db55d3ef533856b1ba82fcd7e576c559 Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Tue, 10 Sep 2019 17:12:57 -0500
Subject: [PATCH] fix for THcDC
---
src/THcDC.cxx | 1789 +++++++++++++++++++++++++------------------------
1 file changed, 912 insertions(+), 877 deletions(-)
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index c09ce44..cad4ce9 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -13,132 +13,125 @@ the number of parameters per plane.
*/
#include "THcDC.h"
-#include "THaEvData.h"
+#include "TClonesArray.h"
+#include "THaApparatus.h"
+#include "THaCutList.h"
#include "THaDetMap.h"
+#include "THaEvData.h"
+#include "THaTrack.h"
+#include "THcDCTrack.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
-#include "THaCutList.h"
+#include "THcHallCSpectrometer.h"
#include "THcParmList.h"
-#include "THcDCTrack.h"
-#include "VarDef.h"
-#include "VarType.h"
-#include "THaTrack.h"
-#include "TClonesArray.h"
#include "TMath.h"
#include "TVectorD.h"
-#include "THaApparatus.h"
-#include "THcHallCSpectrometer.h"
+#include "VarDef.h"
+#include "VarType.h"
-#include <cstring>
#include <cstdio>
#include <cstdlib>
+#include <cstring>
#include <iostream>
-#include "spdlog/spdlog.h"
+#include "spdlog/sinks/basic_file_sink.h" // support for basic file logging
#include "spdlog/sinks/stdout_color_sinks.h" //support for stdout logging
-#include "spdlog/sinks/basic_file_sink.h" // support for basic file logging
-
+#include "spdlog/spdlog.h"
using namespace std;
//_____________________________________________________________________________
-THcDC::THcDC(
- const char* name, const char* description,
- THaApparatus* apparatus ) :
- THaTrackingDetector(name,description,apparatus)
-{
+THcDC::THcDC(const char* name, const char* description, THaApparatus* apparatus)
+ : THaTrackingDetector(name, description, apparatus) {
// Constructor
- fNPlanes = 0; // No planes until we make them
- _sub_logger = _det_logger->clone("THcDC");
- fXCenter = NULL;
- fYCenter = NULL;
- fMinHits = NULL;
- fMaxHits = NULL;
- fMinCombos = NULL;
+ fNPlanes = 0; // No planes until we make them
+ _sub_logger = _det_logger->clone("THcDC");
+ fXCenter = NULL;
+ fYCenter = NULL;
+ fMinHits = NULL;
+ fMaxHits = NULL;
+ fMinCombos = NULL;
fSpace_Point_Criterion = NULL;
- fTdcWinMin = NULL;
- fTdcWinMax = NULL;
- fCentralTime = NULL;
- fNWires = NULL;
- fNChamber = NULL;
- fWireOrder = NULL;
+ fTdcWinMin = NULL;
+ fTdcWinMax = NULL;
+ fCentralTime = NULL;
+ fNWires = NULL;
+ fNChamber = NULL;
+ fWireOrder = NULL;
fDriftTimeSign = NULL;
- fReadoutTB = NULL;
- fReadoutLR = NULL;
-
- fXPos = NULL;
- fYPos = NULL;
- fZPos = NULL;
- fAlphaAngle = NULL;
- fBetaAngle = NULL;
- fGammaAngle = NULL;
- fPitch = NULL;
- fCentralWire = NULL;
+ fReadoutTB = NULL;
+ fReadoutLR = NULL;
+
+ fXPos = NULL;
+ fYPos = NULL;
+ fZPos = NULL;
+ fAlphaAngle = NULL;
+ fBetaAngle = NULL;
+ fGammaAngle = NULL;
+ fPitch = NULL;
+ fCentralWire = NULL;
fPlaneTimeZero = NULL;
- fSigma = NULL;
+ fSigma = NULL;
// These should be set to zero (in a parameter file) in order to
// replicate historical ENGINE behavior
- fFixLR = 1;
+ fFixLR = 1;
fFixPropagationCorrection = 1;
- fProjectToChamber = 0; // Use 1 for SOS chambers
+ fProjectToChamber = 0; // Use 1 for SOS chambers
- fDCTracks = new TClonesArray( "THcDCTrack", 20 );
+ fDCTracks = new TClonesArray("THcDCTrack", 20);
- fNChamHits = 0;
+ fNChamHits = 0;
fPlaneEvents = 0;
- //The version defaults to 0 (old HMS style). 1 is new HMS style and 2 is SHMS style.
+ // The version defaults to 0 (old HMS style). 1 is new HMS style and 2 is SHMS style.
fVersion = 0;
- //Create and return a shared_ptr to a multithreaded console logger.
+ // Create and return a shared_ptr to a multithreaded console logger.
//_logger = spdlog::get("config");
- //if(!_logger) {
+ // if(!_logger) {
// _logger = spdlog::stdout_color_mt("config");
//}
_sub_logger->debug("THcDC constructor");
}
//_____________________________________________________________________________
-void THcDC::Setup(const char* name, const char* description)
-{
+void THcDC::Setup(const char* name, const char* description) {
- Bool_t optional = true;
+ Bool_t optional = true;
// Create the chamber and plane objects using parameters.
static const char* const here = "Setup";
- THaApparatus *app = GetApparatus();
- if(app) {
+ THaApparatus* app = GetApparatus();
+ if (app) {
// cout << app->GetName() << endl;
- fPrefix[0]=tolower(app->GetName()[0]);
- fPrefix[1]='\0';
+ fPrefix[0] = tolower(app->GetName()[0]);
+ fPrefix[1] = '\0';
} else {
cout << "No apparatus found" << endl;
- fPrefix[0]='\0';
+ fPrefix[0] = '\0';
}
// For now, decide chamber style from the spectrometer name.
// Should override with a paramter
- //cout<<"HMS Style??\t"<<fHMSStyleChambers<<endl;
- string planenamelist;
- DBRequest list[]={
- {"dc_num_planes",&fNPlanes, kInt},
- {"dc_num_chambers",&fNChambers, kInt},
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
- {"dc_plane_names",&planenamelist, kString},
- {"dc_version", &fVersion, kInt, 0, optional},
- {"dc_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
- {0}
- };
-
- fTDC_RefTimeCut = 0; // Minimum allowed reference times
- gHcParms->LoadParmValues((DBRequest*)&list,fPrefix);
-
- if(fVersion==0) {
+ // cout<<"HMS Style??\t"<<fHMSStyleChambers<<endl;
+ string planenamelist;
+ DBRequest list[] = {{"dc_num_planes", &fNPlanes, kInt},
+ {"dc_num_chambers", &fNChambers, kInt},
+ {"dc_tdc_time_per_channel", &fNSperChan, kDouble},
+ {"dc_wire_velocity", &fWireVelocity, kDouble},
+ {"dc_plane_names", &planenamelist, kString},
+ {"dc_version", &fVersion, kInt, 0, optional},
+ {"dc_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
+ {0}};
+
+ fTDC_RefTimeCut = 0; // Minimum allowed reference times
+ gHcParms->LoadParmValues((DBRequest*)&list, fPrefix);
+
+ if (fVersion == 0) {
fHMSStyleChambers = 1;
} else {
fHMSStyleChambers = 0;
@@ -146,123 +139,121 @@ void THcDC::Setup(const char* name, const char* description)
_det_logger->info("Plane Name List: {}", planenamelist);
_det_logger->info("Drift Chambers: {} planes in {} chambers", fNPlanes, fNChambers);
- //cout << "Plane Name List: " << planenamelist << endl;
- //cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
+ // cout << "Plane Name List: " << planenamelist << endl;
+ // cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
vector<string> plane_names = vsplit(planenamelist);
- if(plane_names.size() != (UInt_t) fNPlanes) {
- //cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
+ if (plane_names.size() != (UInt_t)fNPlanes) {
+ // cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names
+ // " << plane_names.size() << endl;
_det_logger->error("ERROR: Number of planes {} doesn't agree with number of plane names {}",
- fNPlanes, plane_names.size());
+ fNPlanes, plane_names.size());
// Should quit. Is there an official way to quit?
}
- fPlaneNames = new char* [fNPlanes];
- for(Int_t i=0;i<fNPlanes;i++) {
- fPlaneNames[i] = new char[plane_names[i].length()+1];
+ fPlaneNames = new char*[fNPlanes];
+ for (Int_t i = 0; i < fNPlanes; i++) {
+ fPlaneNames[i] = new char[plane_names[i].length() + 1];
strcpy(fPlaneNames[i], plane_names[i].c_str());
}
- char *desc = new char[strlen(description)+100];
- char *desc1= new char[strlen(description)+100];
+ char* desc = new char[strlen(description) + 100];
+ char* desc1 = new char[strlen(description) + 100];
fPlanes.clear();
- for(Int_t i=0;i<fNPlanes;i++) {
+ for (Int_t i = 0; i < fNPlanes; i++) {
strcpy(desc, description);
strcat(desc, " Plane ");
strcat(desc, fPlaneNames[i]);
- THcDriftChamberPlane* newplane = new THcDriftChamberPlane(fPlaneNames[i], desc, i+1, this);
- if( !newplane or newplane->IsZombie() ) {
- _det_logger->error( "{} Error creating Drift Chamber plane {}. Call expert.", Here(here), name);
+ THcDriftChamberPlane* newplane = new THcDriftChamberPlane(fPlaneNames[i], desc, i + 1, this);
+ if (!newplane or newplane->IsZombie()) {
+ _det_logger->error("{} Error creating Drift Chamber plane {}. Call expert.", Here(here),
+ name);
MakeZombie();
- delete [] desc;
- delete [] desc1;
+ delete[] desc;
+ delete[] desc1;
return;
}
fPlanes.push_back(newplane);
newplane->SetDebug(fDebug);
// cout << "Created Drift Chamber Plane " << fPlaneNames[i] << ", " << desc << endl;
-
}
fChambers.clear();
- for(UInt_t i=0;i<fNChambers;i++) {
- sprintf(desc1,"Ch%d",i+1);
+ for (UInt_t i = 0; i < fNChambers; i++) {
+ sprintf(desc1, "Ch%d", i + 1);
// Should construct a better chamber name
- THcDriftChamber* newchamber = new THcDriftChamber(desc1, desc, i+1, this);
+ THcDriftChamber* newchamber = new THcDriftChamber(desc1, desc, i + 1, this);
fChambers.push_back(newchamber);
- //cout << "Created Drift Chamber " << i+1 << ", " << desc1 << endl;
- _det_logger->info("Created Drift Chamber {}, {}" , i+1 , desc1);
+ // cout << "Created Drift Chamber " << i+1 << ", " << desc1 << endl;
+ _det_logger->info("Created Drift Chamber {}, {}", i + 1, desc1);
newchamber->SetHMSStyleFlag(fHMSStyleChambers); // Tell the chamber its style
}
- delete [] desc;
- delete [] desc1;
+ delete[] desc;
+ delete[] desc1;
}
//_____________________________________________________________________________
-THcDC::THcDC( )
-{
+THcDC::THcDC() {
// Constructor
}
//_____________________________________________________________________________
-THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
-{
+THaAnalysisObject::EStatus THcDC::Init(const TDatime& date) {
// Register the plane objects with the appropriate chambers.
// Trigger ReadDatabase to load the remaining parameters
- Setup(GetName(), GetTitle()); // Create the subdetectors here
+ Setup(GetName(), GetTitle()); // Create the subdetectors here
EffInit();
char EngineDID[] = "xDC";
- EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
- if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
+ EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
+ if (gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0) {
static const char* const here = "Init()";
- _det_logger->error("{} Error filling detectormap for {}.", Here(here), EngineDID );
+ _det_logger->error("{} Error filling detectormap for {}.", Here(here), EngineDID);
return kInitError;
}
// Should probably put this in ReadDatabase as we will know the
// maximum number of hits after setting up the detector map
_det_logger->info("DC tdc ref time cut = {} ", fTDC_RefTimeCut);
- //cout << " DC tdc ref time cut = " << fTDC_RefTimeCut << endl;
- InitHitList(fDetMap, "THcRawDCHit", fDetMap->GetTotNumChan()+1, fTDC_RefTimeCut, 0);
+ // cout << " DC tdc ref time cut = " << fTDC_RefTimeCut << endl;
+ InitHitList(fDetMap, "THcRawDCHit", fDetMap->GetTotNumChan() + 1, fTDC_RefTimeCut, 0);
- CreateMissReportParms(Form("%sdc",fPrefix));
+ CreateMissReportParms(Form("%sdc", fPrefix));
EStatus status;
// This triggers call of ReadDatabase and DefineVariables
- if( (status = THaTrackingDetector::Init( date )) )
- return fStatus=status;
+ if ((status = THaTrackingDetector::Init(date)))
+ return fStatus = status;
// Initialize planes and add them to chambers
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- if((status = fPlanes[ip]->Init( date ))) {
- return fStatus=status;
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ if ((status = fPlanes[ip]->Init(date))) {
+ return fStatus = status;
} else {
- Int_t chamber=fNChamber[ip];
- fChambers[chamber-1]->AddPlane(fPlanes[ip]);
+ Int_t chamber = fNChamber[ip];
+ fChambers[chamber - 1]->AddPlane(fPlanes[ip]);
}
}
// Initialize chambers
- for(UInt_t ic=0;ic<fNChambers;ic++) {
- if((status = fChambers[ic]->Init ( date ))) {
- return fStatus=status;
+ for (UInt_t ic = 0; ic < fNChambers; ic++) {
+ if ((status = fChambers[ic]->Init(date))) {
+ return fStatus = status;
}
}
// Retrieve the fiting coefficients
- fPlaneCoeffs = new Double_t* [fNPlanes];
- for(Int_t ip=0; ip<fNPlanes;ip++) {
+ fPlaneCoeffs = new Double_t*[fNPlanes];
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
fPlaneCoeffs[ip] = fPlanes[ip]->GetPlaneCoef();
}
- fResiduals = new Double_t [fNPlanes];
- fResidualsExclPlane = new Double_t [fNPlanes];
- fWire_hit_did = new Double_t [fNPlanes];
- fWire_hit_should = new Double_t [fNPlanes];
-
+ fResiduals = new Double_t[fNPlanes];
+ fResidualsExclPlane = new Double_t[fNPlanes];
+ fWire_hit_did = new Double_t[fNPlanes];
+ fWire_hit_should = new Double_t[fNPlanes];
// Replace with what we need for Hall C
// const DataDest tmp[NDEST] = {
@@ -270,17 +261,17 @@ THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
// { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
// };
// memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
-
- fPresentP = 0;
- THaVar* vpresent = gHaVars->Find(Form("%s.present",GetApparatus()->GetName()));
- if(vpresent) {
- fPresentP = (Bool_t *) vpresent->GetValuePointer();
+
+ fPresentP = 0;
+ THaVar* vpresent = gHaVars->Find(Form("%s.present", GetApparatus()->GetName()));
+ if (vpresent) {
+ fPresentP = (Bool_t*)vpresent->GetValuePointer();
}
return fStatus = kOK;
}
//_____________________________________________________________________________
-Int_t THcDC::ManualInitTree( TTree* t ){
+Int_t THcDC::ManualInitTree(TTree* t) {
std::string app_name = GetApparatus()->GetName();
std::string det_name = GetName();
std::string branch_name = (app_name + "_" + det_name + "_data");
@@ -292,8 +283,7 @@ Int_t THcDC::ManualInitTree( TTree* t ){
}
//_____________________________________________________________________________
-Int_t THcDC::ReadDatabase( const TDatime& date )
-{
+Int_t THcDC::ReadDatabase(const TDatime& date) {
/**
Read this detector's parameters from the ThcParmList
This function is called by THaDetectorBase::Init() once at the
@@ -302,114 +292,135 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
*/
// static const char* const here = "ReadDatabase()";
- delete [] fXCenter; fXCenter = new Double_t [fNChambers];
- delete [] fYCenter; fYCenter = new Double_t [fNChambers];
- delete [] fMinHits; fMinHits = new Int_t [fNChambers];
- delete [] fMaxHits; fMaxHits = new Int_t [fNChambers];
- delete [] fMinCombos; fMinCombos = new Int_t [fNChambers];
- delete [] fSpace_Point_Criterion; fSpace_Point_Criterion = new Double_t [fNChambers];
-
- delete [] fTdcWinMin; fTdcWinMin = new Int_t [fNPlanes];
- delete [] fTdcWinMax; fTdcWinMax = new Int_t [fNPlanes];
- delete [] fCentralTime; fCentralTime = new Double_t [fNPlanes];
- delete [] fNWires; fNWires = new Int_t [fNPlanes];
- delete [] fNChamber; fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
- delete [] fWireOrder; fWireOrder = new Int_t [fNPlanes]; // Wire readout order
- delete [] fDriftTimeSign; fDriftTimeSign = new Int_t [fNPlanes];
- delete [] fReadoutLR; fReadoutLR = new Int_t [fNPlanes];
- delete [] fReadoutTB; fReadoutTB = new Int_t [fNPlanes];
-
- delete [] fXPos; fXPos = new Double_t [fNPlanes];
- delete [] fYPos; fYPos = new Double_t [fNPlanes];
- delete [] fZPos; fZPos = new Double_t [fNPlanes];
- delete [] fAlphaAngle; fAlphaAngle = new Double_t [fNPlanes];
- delete [] fBetaAngle; fBetaAngle = new Double_t [fNPlanes];
- delete [] fGammaAngle; fGammaAngle = new Double_t [fNPlanes];
- delete [] fPitch; fPitch = new Double_t [fNPlanes];
- delete [] fCentralWire; fCentralWire = new Double_t [fNPlanes];
- delete [] fPlaneTimeZero; fPlaneTimeZero = new Double_t [fNPlanes];
- delete [] fSigma; fSigma = new Double_t [fNPlanes];
+ delete[] fXCenter;
+ fXCenter = new Double_t[fNChambers];
+ delete[] fYCenter;
+ fYCenter = new Double_t[fNChambers];
+ delete[] fMinHits;
+ fMinHits = new Int_t[fNChambers];
+ delete[] fMaxHits;
+ fMaxHits = new Int_t[fNChambers];
+ delete[] fMinCombos;
+ fMinCombos = new Int_t[fNChambers];
+ delete[] fSpace_Point_Criterion;
+ fSpace_Point_Criterion = new Double_t[fNChambers];
+
+ delete[] fTdcWinMin;
+ fTdcWinMin = new Int_t[fNPlanes];
+ delete[] fTdcWinMax;
+ fTdcWinMax = new Int_t[fNPlanes];
+ delete[] fCentralTime;
+ fCentralTime = new Double_t[fNPlanes];
+ delete[] fNWires;
+ fNWires = new Int_t[fNPlanes];
+ delete[] fNChamber;
+ fNChamber = new Int_t[fNPlanes]; // Which chamber is this plane
+ delete[] fWireOrder;
+ fWireOrder = new Int_t[fNPlanes]; // Wire readout order
+ delete[] fDriftTimeSign;
+ fDriftTimeSign = new Int_t[fNPlanes];
+ delete[] fReadoutLR;
+ fReadoutLR = new Int_t[fNPlanes];
+ delete[] fReadoutTB;
+ fReadoutTB = new Int_t[fNPlanes];
+
+ delete[] fXPos;
+ fXPos = new Double_t[fNPlanes];
+ delete[] fYPos;
+ fYPos = new Double_t[fNPlanes];
+ delete[] fZPos;
+ fZPos = new Double_t[fNPlanes];
+ delete[] fAlphaAngle;
+ fAlphaAngle = new Double_t[fNPlanes];
+ delete[] fBetaAngle;
+ fBetaAngle = new Double_t[fNPlanes];
+ delete[] fGammaAngle;
+ fGammaAngle = new Double_t[fNPlanes];
+ delete[] fPitch;
+ fPitch = new Double_t[fNPlanes];
+ delete[] fCentralWire;
+ fCentralWire = new Double_t[fNPlanes];
+ delete[] fPlaneTimeZero;
+ fPlaneTimeZero = new Double_t[fNPlanes];
+ delete[] fSigma;
+ fSigma = new Double_t[fNPlanes];
Bool_t optional = true;
- DBRequest list[]={
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
-
- {"dc_xcenter", fXCenter, kDouble, fNChambers},
- {"dc_ycenter", fYCenter, kDouble, fNChambers},
- {"min_hit", fMinHits, kInt, fNChambers},
- {"max_pr_hits", fMaxHits, kInt, fNChambers},
- {"min_combos", fMinCombos, kInt, fNChambers},
- {"space_point_criterion", fSpace_Point_Criterion, kDouble, fNChambers},
-
- {"dc_tdc_min_win", fTdcWinMin, kInt, (UInt_t)fNPlanes},
- {"dc_tdc_max_win", fTdcWinMax, kInt, (UInt_t)fNPlanes},
- {"dc_central_time", fCentralTime, kDouble, (UInt_t)fNPlanes},
- {"dc_nrwire", fNWires, kInt, (UInt_t)fNPlanes},
- {"dc_chamber_planes", fNChamber, kInt, (UInt_t)fNPlanes},
- {"dc_wire_counting", fWireOrder, kInt, (UInt_t)fNPlanes},
- {"dc_drifttime_sign", fDriftTimeSign, kInt, (UInt_t)fNPlanes},
- {"dc_readoutLR", fReadoutLR, kInt, (UInt_t)fNPlanes, optional},
- {"dc_readoutTB", fReadoutTB, kInt, (UInt_t)fNPlanes, optional},
-
- {"dc_zpos", fZPos, kDouble, (UInt_t)fNPlanes},
- {"dc_alpha_angle", fAlphaAngle, kDouble, (UInt_t)fNPlanes},
- {"dc_beta_angle", fBetaAngle, kDouble, (UInt_t)fNPlanes},
- {"dc_gamma_angle", fGammaAngle, kDouble, (UInt_t)fNPlanes},
- {"dc_pitch", fPitch, kDouble, (UInt_t)fNPlanes},
- {"dc_central_wire", fCentralWire, kDouble, (UInt_t)fNPlanes},
- {"dc_plane_time_zero", fPlaneTimeZero, kDouble, (UInt_t)fNPlanes},
- {"dc_sigma", fSigma, kDouble, (UInt_t)fNPlanes},
- {"single_stub",&fSingleStub, kInt,0,1},
- {"ntracks_max_fp", &fNTracksMaxFP, kInt},
- {"xt_track_criterion", &fXtTrCriterion, kDouble},
- {"yt_track_criterion", &fYtTrCriterion, kDouble},
- {"xpt_track_criterion", &fXptTrCriterion, kDouble},
- {"ypt_track_criterion", &fYptTrCriterion, kDouble},
- {"dc_fix_lr", &fFixLR, kInt},
- {"dc_fix_propcorr", &fFixPropagationCorrection, kInt},
- {"debuglinkstubs", &fdebuglinkstubs, kInt},
- {"debugprintrawdc", &fdebugprintrawdc, kInt},
- {"debugprintdecodeddc", &fdebugprintdecodeddc, kInt},
- {"debugflagpr", &fdebugflagpr, kInt},
- {"debugflagstubs", &fdebugflagstubs, kInt},
- {"debugtrackprint", &fdebugtrackprint , kInt},
- {0}
- };
- fSingleStub=0;
- for(Int_t ip=0; ip<fNPlanes;ip++) {
+ DBRequest list[] = {{"dc_tdc_time_per_channel", &fNSperChan, kDouble},
+ {"dc_wire_velocity", &fWireVelocity, kDouble},
+
+ {"dc_xcenter", fXCenter, kDouble, fNChambers},
+ {"dc_ycenter", fYCenter, kDouble, fNChambers},
+ {"min_hit", fMinHits, kInt, fNChambers},
+ {"max_pr_hits", fMaxHits, kInt, fNChambers},
+ {"min_combos", fMinCombos, kInt, fNChambers},
+ {"space_point_criterion", fSpace_Point_Criterion, kDouble, fNChambers},
+
+ {"dc_tdc_min_win", fTdcWinMin, kInt, (UInt_t)fNPlanes},
+ {"dc_tdc_max_win", fTdcWinMax, kInt, (UInt_t)fNPlanes},
+ {"dc_central_time", fCentralTime, kDouble, (UInt_t)fNPlanes},
+ {"dc_nrwire", fNWires, kInt, (UInt_t)fNPlanes},
+ {"dc_chamber_planes", fNChamber, kInt, (UInt_t)fNPlanes},
+ {"dc_wire_counting", fWireOrder, kInt, (UInt_t)fNPlanes},
+ {"dc_drifttime_sign", fDriftTimeSign, kInt, (UInt_t)fNPlanes},
+ {"dc_readoutLR", fReadoutLR, kInt, (UInt_t)fNPlanes, optional},
+ {"dc_readoutTB", fReadoutTB, kInt, (UInt_t)fNPlanes, optional},
+
+ {"dc_zpos", fZPos, kDouble, (UInt_t)fNPlanes},
+ {"dc_alpha_angle", fAlphaAngle, kDouble, (UInt_t)fNPlanes},
+ {"dc_beta_angle", fBetaAngle, kDouble, (UInt_t)fNPlanes},
+ {"dc_gamma_angle", fGammaAngle, kDouble, (UInt_t)fNPlanes},
+ {"dc_pitch", fPitch, kDouble, (UInt_t)fNPlanes},
+ {"dc_central_wire", fCentralWire, kDouble, (UInt_t)fNPlanes},
+ {"dc_plane_time_zero", fPlaneTimeZero, kDouble, (UInt_t)fNPlanes},
+ {"dc_sigma", fSigma, kDouble, (UInt_t)fNPlanes},
+ {"single_stub", &fSingleStub, kInt, 0, 1},
+ {"ntracks_max_fp", &fNTracksMaxFP, kInt},
+ {"xt_track_criterion", &fXtTrCriterion, kDouble},
+ {"yt_track_criterion", &fYtTrCriterion, kDouble},
+ {"xpt_track_criterion", &fXptTrCriterion, kDouble},
+ {"ypt_track_criterion", &fYptTrCriterion, kDouble},
+ {"dc_fix_lr", &fFixLR, kInt},
+ {"dc_fix_propcorr", &fFixPropagationCorrection, kInt},
+ {"debuglinkstubs", &fdebuglinkstubs, kInt},
+ {"debugprintrawdc", &fdebugprintrawdc, kInt},
+ {"debugprintdecodeddc", &fdebugprintdecodeddc, kInt},
+ {"debugflagpr", &fdebugflagpr, kInt},
+ {"debugflagstubs", &fdebugflagstubs, kInt},
+ {"debugtrackprint", &fdebugtrackprint, kInt},
+ {0}};
+ fSingleStub = 0;
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
fReadoutLR[ip] = 0.0;
fReadoutTB[ip] = 0.0;
- }
-
+ }
- gHcParms->LoadParmValues((DBRequest*)&list,fPrefix);
+ gHcParms->LoadParmValues((DBRequest*)&list, fPrefix);
- //Set the default plane x,y positions to those of the chamber
- for(Int_t ip=0; ip<fNPlanes;ip++) {
- fXPos[ip] = fXCenter[GetNChamber(ip+1)-1];
- fYPos[ip] = fYCenter[GetNChamber(ip+1)-1];
- }
+ // Set the default plane x,y positions to those of the chamber
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ fXPos[ip] = fXCenter[GetNChamber(ip + 1) - 1];
+ fYPos[ip] = fYCenter[GetNChamber(ip + 1) - 1];
+ }
- //Load the x,y positions of the planes if they exist (overwrites defaults)
- DBRequest listOpt[]={
- {"dc_xpos", fXPos, kDouble, (UInt_t)fNPlanes, optional},
- {"dc_ypos", fYPos, kDouble, (UInt_t)fNPlanes, optional},
- {0}
- };
- gHcParms->LoadParmValues((DBRequest*)&listOpt,fPrefix);
- if(fNTracksMaxFP <= 0) fNTracksMaxFP = 10;
+ // Load the x,y positions of the planes if they exist (overwrites defaults)
+ DBRequest listOpt[] = {{"dc_xpos", fXPos, kDouble, (UInt_t)fNPlanes, optional},
+ {"dc_ypos", fYPos, kDouble, (UInt_t)fNPlanes, optional},
+ {0}};
+ gHcParms->LoadParmValues((DBRequest*)&listOpt, fPrefix);
+ if (fNTracksMaxFP <= 0)
+ fNTracksMaxFP = 10;
// if(fNTracksMaxFP > HNRACKS_MAX) fNTracksMaxFP = NHTRACKS_MAX;
-
+
std::string plane_counts_string;
- //cout << "Plane counts:";
- for(Int_t i=0;i<fNPlanes;i++) {
- //cout << " " << fNWires[i];
+ // cout << "Plane counts:";
+ for (Int_t i = 0; i < fNPlanes; i++) {
+ // cout << " " << fNWires[i];
plane_counts_string += std::string(" ");
plane_counts_string += std::to_string(fNWires[i]);
}
- //cout << endl;
+ // cout << endl;
_det_logger->info("Plane counts: {}", plane_counts_string);
fIsInit = true;
@@ -418,148 +429,173 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
}
//_____________________________________________________________________________
-Int_t THcDC::DefineVariables( EMode mode )
-{
+Int_t THcDC::DefineVariables(EMode mode) {
/**
Initialize global variables for histograms and Root tree
*/
- if( mode == kDefine && fIsSetup ) return kOK;
- fIsSetup = ( mode == kDefine );
+ if (mode == kDefine && fIsSetup)
+ return kOK;
+ fIsSetup = (mode == kDefine);
// Register variables in global list
RVarDef vars[] = {
- { "stubtest", "stub test", "fStubTest" },
- { "nhit", "Number of DC hits", "fNhits" },
- { "tnhit", "Number of good DC hits", "fNthits" },
- { "trawhit", "Number of true raw DC hits", "fN_True_RawHits" },
- { "ntrack", "Number of Tracks", "fNDCTracks" },
- { "nsp", "Number of Space Points", "fNSp" },
- { "track_nsp", "Number of spacepoints in track", "fDCTracks.THcDCTrack.GetNSpacePoints()"},
- { "x", "X at focal plane", "fDCTracks.THcDCTrack.GetX()"},
- { "y", "Y at focal plane", "fDCTracks.THcDCTrack.GetY()"},
- { "xp", "XP at focal plane", "fDCTracks.THcDCTrack.GetXP()"},
- { "yp", "YP at focal plane", "fDCTracks.THcDCTrack.GetYP()"},
- { "x_fp", "X at focal plane (golden track)", "fX_fp_best"},
- { "y_fp", "Y at focal plane( golden track)", "fY_fp_best"},
- { "xp_fp", "XP at focal plane (golden track)", "fXp_fp_best"},
- { "yp_fp", "YP at focal plane(golden track) ", "fYp_fp_best"},
- { "chisq", "chisq/dof (golden track) ", "fChisq_best"},
- { "sp1_id", " (golden track) ", "fSp1_ID_best"},
- { "sp2_id", " (golden track) ", "fSp2_ID_best"},
- { "InsideDipoleExit", " ","fInSideDipoleExit_best"},
- { "gtrack_nsp", " Number of space points in golden track ", "fNsp_best"},
- { "residual", "Residuals", "fResiduals"},
- { "residualExclPlane", "Residuals", "fResidualsExclPlane"},
- { "wireHitDid","Wire did have matched track hit", "fWire_hit_did"},
- { "wireHitShould", "Wire should have matched track hit", "fWire_hit_should"},
- { 0 }
- };
- return DefineVarsFromList( vars, mode );
-
+ {"stubtest", "stub test", "fStubTest"},
+ {"nhit", "Number of DC hits", "fNhits"},
+ {"tnhit", "Number of good DC hits", "fNthits"},
+ {"trawhit", "Number of true raw DC hits", "fN_True_RawHits"},
+ {"ntrack", "Number of Tracks", "fNDCTracks"},
+ {"nsp", "Number of Space Points", "fNSp"},
+ {"track_nsp", "Number of spacepoints in track", "fDCTracks.THcDCTrack.GetNSpacePoints()"},
+ {"x", "X at focal plane", "fDCTracks.THcDCTrack.GetX()"},
+ {"y", "Y at focal plane", "fDCTracks.THcDCTrack.GetY()"},
+ {"xp", "XP at focal plane", "fDCTracks.THcDCTrack.GetXP()"},
+ {"yp", "YP at focal plane", "fDCTracks.THcDCTrack.GetYP()"},
+ {"x_fp", "X at focal plane (golden track)", "fX_fp_best"},
+ {"y_fp", "Y at focal plane( golden track)", "fY_fp_best"},
+ {"xp_fp", "XP at focal plane (golden track)", "fXp_fp_best"},
+ {"yp_fp", "YP at focal plane(golden track) ", "fYp_fp_best"},
+ {"chisq", "chisq/dof (golden track) ", "fChisq_best"},
+ {"sp1_id", " (golden track) ", "fSp1_ID_best"},
+ {"sp2_id", " (golden track) ", "fSp2_ID_best"},
+ {"InsideDipoleExit", " ", "fInSideDipoleExit_best"},
+ {"gtrack_nsp", " Number of space points in golden track ", "fNsp_best"},
+ {"residual", "Residuals", "fResiduals"},
+ {"residualExclPlane", "Residuals", "fResidualsExclPlane"},
+ {"wireHitDid", "Wire did have matched track hit", "fWire_hit_did"},
+ {"wireHitShould", "Wire should have matched track hit", "fWire_hit_should"},
+ {0}};
+ return DefineVarsFromList(vars, mode);
}
//_____________________________________________________________________________
-THcDC::~THcDC()
-{
+THcDC::~THcDC() {
// Destructor. Remove variables from global list.
- if( fIsSetup )
+ if (fIsSetup)
RemoveVariables();
- if( fIsInit )
+ if (fIsInit)
DeleteArrays();
// Delete the plane objects
- for (vector<THcDriftChamberPlane*>::iterator ip = fPlanes.begin();
- ip != fPlanes.end(); ++ip) delete *ip;
+ for (vector<THcDriftChamberPlane*>::iterator ip = fPlanes.begin(); ip != fPlanes.end(); ++ip)
+ delete *ip;
// Delete the chamber objects
- for (vector<THcDriftChamber*>::iterator ip = fChambers.begin();
- ip != fChambers.end(); ++ip) delete *ip;
+ for (vector<THcDriftChamber*>::iterator ip = fChambers.begin(); ip != fChambers.end(); ++ip)
+ delete *ip;
delete fDCTracks;
}
//_____________________________________________________________________________
-void THcDC::DeleteArrays()
-{
+void THcDC::DeleteArrays() {
// Delete member arrays. Used by destructor.
- delete [] fXCenter; fXCenter = NULL;
- delete [] fYCenter; fYCenter = NULL;
- delete [] fMinHits; fMinHits = NULL;
- delete [] fMaxHits; fMaxHits = NULL;
- delete [] fMinCombos; fMinCombos = NULL;
- delete [] fSpace_Point_Criterion; fSpace_Point_Criterion = NULL;
-
- delete [] fTdcWinMin; fTdcWinMin = NULL;
- delete [] fTdcWinMax; fTdcWinMax = NULL;
- delete [] fCentralTime; fCentralTime = NULL;
- delete [] fNWires; fNWires = NULL;
- delete [] fNChamber; fNChamber = NULL;
- delete [] fWireOrder; fWireOrder = NULL;
- delete [] fDriftTimeSign; fDriftTimeSign = NULL;
- delete [] fReadoutLR; fReadoutLR = NULL;
- delete [] fReadoutTB; fReadoutTB = NULL;
-
- delete [] fXPos; fXPos = NULL;
- delete [] fYPos; fYPos = NULL;
- delete [] fZPos; fZPos = NULL;
- delete [] fAlphaAngle; fAlphaAngle = NULL;
- delete [] fBetaAngle; fBetaAngle = NULL;
- delete [] fGammaAngle; fGammaAngle = NULL;
- delete [] fPitch; fPitch = NULL;
- delete [] fCentralWire; fCentralWire = NULL;
- delete [] fPlaneTimeZero; fPlaneTimeZero = NULL;
- delete [] fSigma; fSigma = NULL;
+ delete[] fXCenter;
+ fXCenter = NULL;
+ delete[] fYCenter;
+ fYCenter = NULL;
+ delete[] fMinHits;
+ fMinHits = NULL;
+ delete[] fMaxHits;
+ fMaxHits = NULL;
+ delete[] fMinCombos;
+ fMinCombos = NULL;
+ delete[] fSpace_Point_Criterion;
+ fSpace_Point_Criterion = NULL;
+
+ delete[] fTdcWinMin;
+ fTdcWinMin = NULL;
+ delete[] fTdcWinMax;
+ fTdcWinMax = NULL;
+ delete[] fCentralTime;
+ fCentralTime = NULL;
+ delete[] fNWires;
+ fNWires = NULL;
+ delete[] fNChamber;
+ fNChamber = NULL;
+ delete[] fWireOrder;
+ fWireOrder = NULL;
+ delete[] fDriftTimeSign;
+ fDriftTimeSign = NULL;
+ delete[] fReadoutLR;
+ fReadoutLR = NULL;
+ delete[] fReadoutTB;
+ fReadoutTB = NULL;
+
+ delete[] fXPos;
+ fXPos = NULL;
+ delete[] fYPos;
+ fYPos = NULL;
+ delete[] fZPos;
+ fZPos = NULL;
+ delete[] fAlphaAngle;
+ fAlphaAngle = NULL;
+ delete[] fBetaAngle;
+ fBetaAngle = NULL;
+ delete[] fGammaAngle;
+ fGammaAngle = NULL;
+ delete[] fPitch;
+ fPitch = NULL;
+ delete[] fCentralWire;
+ fCentralWire = NULL;
+ delete[] fPlaneTimeZero;
+ fPlaneTimeZero = NULL;
+ delete[] fSigma;
+ fSigma = NULL;
// Efficiency arrays
- delete [] fNChamHits; fNChamHits = NULL;
- delete [] fPlaneEvents; fPlaneEvents = NULL;
-
- for( Int_t i = 0; i<fNPlanes; ++i )
- delete [] fPlaneNames[i];
- delete [] fPlaneNames;
-
- delete [] fPlaneCoeffs; fPlaneCoeffs = 0;
- delete [] fResiduals; fResiduals = 0;
- delete [] fResidualsExclPlane; fResidualsExclPlane = 0;
- delete [] fWire_hit_did; fWire_hit_did = 0;
- delete [] fWire_hit_should; fWire_hit_should = 0;
+ delete[] fNChamHits;
+ fNChamHits = NULL;
+ delete[] fPlaneEvents;
+ fPlaneEvents = NULL;
+
+ for (Int_t i = 0; i < fNPlanes; ++i)
+ delete[] fPlaneNames[i];
+ delete[] fPlaneNames;
+
+ delete[] fPlaneCoeffs;
+ fPlaneCoeffs = 0;
+ delete[] fResiduals;
+ fResiduals = 0;
+ delete[] fResidualsExclPlane;
+ fResidualsExclPlane = 0;
+ delete[] fWire_hit_did;
+ fWire_hit_did = 0;
+ delete[] fWire_hit_should;
+ fWire_hit_should = 0;
}
//_____________________________________________________________________________
-inline
-void THcDC::ClearEvent()
-{
+inline void THcDC::ClearEvent() {
// Reset per-event data.
- fStubTest = 0;
- fNhits = 0;
- fNthits = 0;
- fN_True_RawHits=0;
- fX_fp_best=-10000.;
- fY_fp_best=-10000.;
- fXp_fp_best=-10000.;
- fYp_fp_best=-10000.;
- fChisq_best=kBig;
- fNsp_best=0;
+ fStubTest = 0;
+ fNhits = 0;
+ fNthits = 0;
+ fN_True_RawHits = 0;
+ fX_fp_best = -10000.;
+ fY_fp_best = -10000.;
+ fXp_fp_best = -10000.;
+ fYp_fp_best = -10000.;
+ fChisq_best = kBig;
+ fNsp_best = 0;
fInSideDipoleExit_best = kTRUE;
- for(UInt_t i=0;i<fNChambers;i++) {
+ for (UInt_t i = 0; i < fNChambers; i++) {
fChambers[i]->Clear();
}
- for(Int_t i=0;i<fNPlanes;i++) {
- fResiduals[i] = 1000.0;
+ for (Int_t i = 0; i < fNPlanes; i++) {
+ fResiduals[i] = 1000.0;
fResidualsExclPlane[i] = 1000.0;
- fWire_hit_did[i] = 1000.0;
- fWire_hit_should[i] = 1000.0;
+ fWire_hit_did[i] = 1000.0;
+ fWire_hit_should[i] = 1000.0;
}
// fTrackProj->Clear();
}
//_____________________________________________________________________________
-Int_t THcDC::Decode( const THaEvData& evdata )
-{
+Int_t THcDC::Decode(const THaEvData& evdata) {
/**
Decode event into hit list.
Pass hit list to the planes.
@@ -567,54 +603,52 @@ Int_t THcDC::Decode( const THaEvData& evdata )
*/
ClearEvent();
Int_t num_event = evdata.GetEvNum();
- if (fdebugprintrawdc ||fdebugprintdecodeddc || fdebuglinkstubs || fdebugtrackprint) cout << " event num = " << num_event << endl;
+ if (fdebugprintrawdc || fdebugprintdecodeddc || fdebuglinkstubs || fdebugtrackprint)
+ cout << " event num = " << num_event << endl;
// Get the Hall C style hitlist (fRawHitList) for this event
- Bool_t present = kTRUE; // Suppress reference time warnings
- if(fPresentP) { // if this spectrometer not part of trigger
+ Bool_t present = kTRUE; // Suppress reference time warnings
+ if (fPresentP) { // if this spectrometer not part of trigger
present = *fPresentP;
}
fNhits = DecodeToHitList(evdata, !present);
-
-
- if(!gHaCuts->Result("Pedestal_event")) {
+ if (!gHaCuts->Result("Pedestal_event")) {
// Let each plane get its hits
Int_t nexthit = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
fN_True_RawHits += fPlanes[ip]->GetNRawhits();
-
}
// fRawHitList is TClones array of THcRawDCHit objects
- Int_t counter=0;
+ Int_t counter = 0;
if (fdebugprintrawdc) {
- cout << " RAW_TOT_HITS = " << fNRawHits << endl;
- cout << " Hit # " << "Plane " << " Wire " << " Raw TDC " << endl;
- for(UInt_t ihit = 0; ihit < fNRawHits ; ihit++) {
- THcRawDCHit* hit = (THcRawDCHit *) fRawHitList->At(ihit);
- for(UInt_t imhit = 0; imhit < hit->GetRawTdcHit().GetNHits(); imhit++) {
- counter++;
- cout << counter << " " << hit->fPlane << " " << hit->fCounter << " " << hit->GetRawTdcHit().GetTimeRaw(imhit) << endl;
- }
+ cout << " RAW_TOT_HITS = " << fNRawHits << endl;
+ cout << " Hit # "
+ << "Plane "
+ << " Wire "
+ << " Raw TDC " << endl;
+ for (UInt_t ihit = 0; ihit < fNRawHits; ihit++) {
+ THcRawDCHit* hit = (THcRawDCHit*)fRawHitList->At(ihit);
+ for (UInt_t imhit = 0; imhit < hit->GetRawTdcHit().GetNHits(); imhit++) {
+ counter++;
+ cout << counter << " " << hit->fPlane << " " << hit->fCounter << " "
+ << hit->GetRawTdcHit().GetTimeRaw(imhit) << endl;
+ }
}
cout << endl;
}
- Eff(); // Accumlate statistics
+ Eff(); // Accumlate statistics
}
return fNhits;
}
//_____________________________________________________________________________
-Int_t THcDC::ApplyCorrections( void )
-{
- return(0);
-}
+Int_t THcDC::ApplyCorrections(void) { return (0); }
//_____________________________________________________________________________
-Int_t THcDC::CoarseTrack( TClonesArray& tracks )
-{
+Int_t THcDC::CoarseTrack(TClonesArray& tracks) {
/**
Find a set of tracks through the drift chambers and put them
into the tracks TClonesArray.
@@ -622,48 +656,50 @@ Int_t THcDC::CoarseTrack( TClonesArray& tracks )
*/
// Subtract starttimes from each plane hit
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- fPlanes[ip]->SubtractStartTime();
- }
+ for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ fPlanes[ip]->SubtractStartTime();
+ }
//
- // Let each chamber get its hits
- for(UInt_t ic=0;ic<fNChambers;ic++) {
- fChambers[ic]->ProcessHits();
- fNthits += fChambers[ic]->GetNHits();
- if (fdebugprintdecodeddc)fChambers[ic]->PrintDecode();
- }
- //
- for(UInt_t i=0;i<fNChambers;i++) {
+ // Let each chamber get its hits
+ for (UInt_t ic = 0; ic < fNChambers; ic++) {
+ fChambers[ic]->ProcessHits();
+ fNthits += fChambers[ic]->GetNHits();
+ if (fdebugprintdecodeddc)
+ fChambers[ic]->PrintDecode();
+ }
+ //
+ for (UInt_t i = 0; i < fNChambers; i++) {
fChambers[i]->FindSpacePoints();
fChambers[i]->CorrectHitTimes();
fChambers[i]->LeftRight();
}
- if (fdebugflagstubs) PrintSpacePoints();
- if (fdebugflagstubs) PrintStubs();
+ if (fdebugflagstubs)
+ PrintSpacePoints();
+ if (fdebugflagstubs)
+ PrintStubs();
// Now link the stubs between chambers
LinkStubs();
- if(fNDCTracks > 0) {
- TrackFit();
+ if (fNDCTracks > 0) {
+ TrackFit();
// Copy tracks into podd tracks list
- for(UInt_t itrack=0;itrack<fNDCTracks;itrack++) {
+ for (UInt_t itrack = 0; itrack < fNDCTracks; itrack++) {
THaTrack* theTrack = NULL;
- theTrack = AddTrack(tracks, 0.0, 0.0, 0.0, 0.0); // Leaving off trackID
+ theTrack = AddTrack(tracks, 0.0, 0.0, 0.0, 0.0); // Leaving off trackID
// Should we add stubs with AddCluster? Could we do this
// by having stubs inherit from cluster
- THcDCTrack *tr = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
+ THcDCTrack* tr = static_cast<THcDCTrack*>(fDCTracks->At(itrack));
theTrack->Set(tr->GetX(), tr->GetY(), tr->GetXP(), tr->GetYP());
- theTrack->SetFlag((UInt_t) 0);
+ theTrack->SetFlag((UInt_t)0);
// Need to look at how engine does chi2 and track selection. Reduced?
- theTrack->SetChi2(tr->GetChisq(),tr->GetNFree());
+ theTrack->SetChi2(tr->GetChisq(), tr->GetNFree());
// CalcFocalPlaneCoords. Aren't our tracks already in focal plane coords
// We should have some kind of track ID so that the THaTrack can be
// associate back with the DC track
// Assign the track number
- theTrack->SetTrkNum(itrack+1);
+ theTrack->SetTrkNum(itrack + 1);
}
- }
-
+ }
ApplyCorrections();
@@ -671,68 +707,66 @@ Int_t THcDC::CoarseTrack( TClonesArray& tracks )
}
//_____________________________________________________________________________
-Int_t THcDC::FineTrack( TClonesArray& tracks )
-{
-
- return 0;
-}
+Int_t THcDC::FineTrack(TClonesArray& tracks) { return 0; }
//
-void THcDC::SetFocalPlaneBestTrack(Int_t golden_track_index)
-{
- THcDCTrack *tr1 = static_cast<THcDCTrack*>( fDCTracks->At(golden_track_index));
- fX_fp_best=tr1->GetX();
- fY_fp_best=tr1->GetY();
- fXp_fp_best=tr1->GetXP();
- fYp_fp_best=tr1->GetYP();
- THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
- fInSideDipoleExit_best = app->InsideDipoleExitWindow(fX_fp_best, fXp_fp_best ,fY_fp_best,fYp_fp_best);
- fSp1_ID_best=tr1->GetSp1_ID();
- fSp2_ID_best=tr1->GetSp2_ID();
- fChisq_best=tr1->GetChisq();
- fNsp_best=tr1->GetNSpacePoints();
- for (UInt_t ihit = 0; ihit < UInt_t (tr1->GetNHits()); ihit++) {
- THcDCHit *hit = tr1->GetHit(ihit);
- Int_t plane = hit->GetPlaneNum() - 1;
- _basic_data._Residuals[plane] = tr1->GetResidual(plane);
- fResiduals[plane] = tr1->GetResidual(plane);
- _basic_data._ResidualsExclPlane[plane] = tr1->GetResidualExclPlane(plane);
- fResidualsExclPlane[plane] = tr1->GetResidualExclPlane(plane);
- }
- EfficiencyPerWire(golden_track_index);
+void THcDC::SetFocalPlaneBestTrack(Int_t golden_track_index) {
+ THcDCTrack* tr1 = static_cast<THcDCTrack*>(fDCTracks->At(golden_track_index));
+ fX_fp_best = tr1->GetX();
+ fY_fp_best = tr1->GetY();
+ fXp_fp_best = tr1->GetXP();
+ fYp_fp_best = tr1->GetYP();
+ THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+ fInSideDipoleExit_best =
+ app->InsideDipoleExitWindow(fX_fp_best, fXp_fp_best, fY_fp_best, fYp_fp_best);
+ fSp1_ID_best = tr1->GetSp1_ID();
+ fSp2_ID_best = tr1->GetSp2_ID();
+ fChisq_best = tr1->GetChisq();
+ fNsp_best = tr1->GetNSpacePoints();
+ for (UInt_t ihit = 0; ihit < UInt_t(tr1->GetNHits()); ihit++) {
+ THcDCHit* hit = tr1->GetHit(ihit);
+ Int_t plane = hit->GetPlaneNum() - 1;
+ _basic_data._Residuals[plane] = tr1->GetResidual(plane);
+ fResiduals[plane] = tr1->GetResidual(plane);
+ _basic_data._ResidualsExclPlane[plane] = tr1->GetResidualExclPlane(plane);
+ fResidualsExclPlane[plane] = tr1->GetResidualExclPlane(plane);
+ }
+ EfficiencyPerWire(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;
+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;
}
}
//
-void THcDC::PrintSpacePoints()
-{
- for(UInt_t ich=0;ich<fNChambers;ich++) {
- printf("%s %2d %s %3d %s %3d \n"," chamber = ",fChambers[ich]->GetChamberNum()," number of hits = ",fChambers[ich]->GetNHits()," number of spacepoints = ",fChambers[ich]->GetNSpacePoints());
- printf("%6s %-8s %-8s %6s %6s %10s \n"," "," "," ","Number","Number","Plane Wire");
- printf("%6s %-8s %-8s %6s %6s %10s \n","Point","x","y"," hits ","combos"," for each hit");
+void THcDC::PrintSpacePoints() {
+ for (UInt_t ich = 0; ich < fNChambers; ich++) {
+ printf("%s %2d %s %3d %s %3d \n", " chamber = ", fChambers[ich]->GetChamberNum(),
+ " number of hits = ", fChambers[ich]->GetNHits(),
+ " number of spacepoints = ", fChambers[ich]->GetNSpacePoints());
+ printf("%6s %-8s %-8s %6s %6s %10s \n", " ", " ", " ", "Number", "Number", "Plane Wire");
+ 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++) {
+ 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("%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");
}
@@ -740,24 +774,23 @@ void THcDC::PrintSpacePoints()
}
//
//
-void THcDC::PrintStubs()
-{
- for(UInt_t ich=0;ich<fNChambers;ich++) {
- printf("%s %3d \n"," Stub fit results Chamber = ",ich+1);
- printf("%-5s %-18s %-18s %-18s %-18s\n","point","x_t","y_t","xp_t","yp_t");
- printf("%-5s %-18s %-18s %-18s %-18s\n"," ","[cm]","[cm]","[cm]","[cm]");
+void THcDC::PrintStubs() {
+ for (UInt_t ich = 0; ich < fNChambers; ich++) {
+ printf("%s %3d \n", " Stub fit results Chamber = ", ich + 1);
+ printf("%-5s %-18s %-18s %-18s %-18s\n", "point", "x_t", "y_t", "xp_t", "yp_t");
+ 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]);
+ 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]);
}
}
}
//
//_____________________________________________________________________________
-void THcDC::LinkStubs()
-{
+void THcDC::LinkStubs() {
/**
The logic is
0) Put all space points in a single list
@@ -774,194 +807,215 @@ void THcDC::LinkStubs()
*/
std::vector<THcSpacePoint*> fSp;
- fNSp=0;
+ fNSp = 0;
fSp.clear();
fSp.reserve(100);
- fNDCTracks=0; // Number of Focal Plane tracks found
+ fNDCTracks = 0; // Number of Focal Plane tracks found
fDCTracks->Delete();
// Make a vector of pointers to the SpacePoints
- //if (fChambers[0]->GetNSpacePoints()+fChambers[1]->GetNSpacePoints()>10) return;
+ // if (fChambers[0]->GetNSpacePoints()+fChambers[1]->GetNSpacePoints()>10) return;
- for(UInt_t ich=0;ich<fNChambers;ich++) {
- Int_t nchamber=fChambers[ich]->GetChamberNum();
+ for (UInt_t ich = 0; ich < fNChambers; ich++) {
+ Int_t nchamber = fChambers[ich]->GetChamberNum();
TClonesArray* spacepointarray = fChambers[ich]->GetSpacePointsP();
- for(Int_t isp=0;isp<fChambers[ich]->GetNSpacePoints();isp++) {
+ for (Int_t isp = 0; isp < fChambers[ich]->GetNSpacePoints(); isp++) {
fSp.push_back(static_cast<THcSpacePoint*>(spacepointarray->At(isp)));
- fSp[fNSp]->fNChamber = nchamber;
+ fSp[fNSp]->fNChamber = nchamber;
fSp[fNSp]->fNChamber_spnum = isp;
fNSp++;
- if (ich==0 && fNSp>50) break;
- if (fNSp>100) break;
+ if (ich == 0 && fNSp > 50)
+ break;
+ if (fNSp > 100)
+ break;
}
}
- Double_t stubminx = 999999;
- Double_t stubminy = 999999;
+ Double_t stubminx = 999999;
+ Double_t stubminy = 999999;
Double_t stubminxp = 999999;
Double_t stubminyp = 999999;
- Int_t stub_tracks[MAXTRACKS];
- if(fSingleStub==0) {
- for(Int_t isp1=0;isp1<fNSp-1;isp1++) { // isp1 is index/id in total list of space points
- THcSpacePoint* sp1 = fSp[isp1];
- Int_t sptracks=0;
+ Int_t stub_tracks[MAXTRACKS];
+ if (fSingleStub == 0) {
+ for (Int_t isp1 = 0; isp1 < fNSp - 1;
+ isp1++) { // isp1 is index/id in total list of space points
+ THcSpacePoint* sp1 = fSp[isp1];
+ Int_t sptracks = 0;
// Now make sure this sp is not already used in a sp.
// Could this be done by having a sp point to the track it is in?
- Int_t tryflag=1;
- for(UInt_t itrack=0;itrack<fNDCTracks;itrack++) {
- THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
- for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
- // isp is index into list of space points attached to theDCTrack
- if(theDCTrack->GetSpacePoint(isp) == sp1) {
- tryflag=0;
- }
- }
+ Int_t tryflag = 1;
+ for (UInt_t itrack = 0; itrack < fNDCTracks; itrack++) {
+ THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(itrack));
+ for (Int_t isp = 0; isp < theDCTrack->GetNSpacePoints(); isp++) {
+ // isp is index into list of space points attached to theDCTrack
+ if (theDCTrack->GetSpacePoint(isp) == sp1) {
+ tryflag = 0;
+ }
+ }
}
- if(tryflag) { // SP not already part of a track
- Int_t newtrack=1;
- for(Int_t isp2=isp1+1;isp2<fNSp;isp2++) {
- THcSpacePoint* sp2=fSp[isp2];
- if(sp1->fNChamber!=sp2->fNChamber&&sp1->GetSetStubFlag()&&sp2->GetSetStubFlag()) {
- Double_t *spstub1=sp1->GetStubP();
- Double_t *spstub2=sp2->GetStubP();
- Double_t dposx = spstub1[0] - spstub2[0];
- Double_t dposy;
- if(fProjectToChamber) { // From SOS s_link_stubs
- // Since single chamber resolution is ~50mr, and the maximum y`
- // angle is about 30mr, use differenece between y AT CHAMBERS, rather
- // than at focal plane. (Project back to chamber, to take out y' uncertainty)
- // (Should this be done for SHMS and HMS too?)
- Double_t y1=spstub1[1]+fChambers[sp1->fNChamber]->GetZPos()*spstub1[3];
- Double_t y2=spstub2[1]+fChambers[sp2->fNChamber]->GetZPos()*spstub2[3];
- dposy = y1-y2;
- } else {
- dposy = spstub1[1] - spstub2[1];
- }
- Double_t dposxp = spstub1[2] - spstub2[2];
- Double_t dposyp = spstub1[3] - spstub2[3];
-
- // What is the point of saving these stubmin values. They
- // Don't seem to be used anywhere except that they can be
- // printed out if hbypass_track_eff_files is zero.
- if(TMath::Abs(dposx)<TMath::Abs(stubminx)) stubminx = dposx;
- if(TMath::Abs(dposy)<TMath::Abs(stubminy)) stubminy = dposy;
- if(TMath::Abs(dposxp)<TMath::Abs(stubminxp)) stubminxp = dposxp;
- if(TMath::Abs(dposyp)<TMath::Abs(stubminyp)) stubminyp = dposyp;
-
- // if hbypass_track_eff_files == 0 then
- // Print out each stubminX that is less that its criterion
-
- if((TMath::Abs(dposx) < fXtTrCriterion)
- && (TMath::Abs(dposy) < fYtTrCriterion)
- && (TMath::Abs(dposxp) < fXptTrCriterion)
- && (TMath::Abs(dposyp) < fYptTrCriterion)) {
- if(newtrack) {
- assert(sptracks==0);
- fStubTest = 1;
- //stubtest=1; Used in h_track_tests.f
- // Make a new track if there are not to many
- if(fNDCTracks < fNTracksMaxFP) {
- sptracks=0; // Number of tracks with this seed
- stub_tracks[sptracks++] = fNDCTracks;
- THcDCTrack *theDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
- theDCTrack->AddSpacePoint(sp1);
- theDCTrack->AddSpacePoint(sp2);
- if (sp1->fNChamber==1) theDCTrack->SetSp1_ID(sp1->fNChamber_spnum);
- if (sp1->fNChamber==2) theDCTrack->SetSp2_ID(sp1->fNChamber_spnum);
- if (sp2->fNChamber==1) theDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
- if (sp2->fNChamber==2) theDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
- newtrack = 0; // Make no more tracks in this loop
- // (But could replace a SP?)
- } else {
+ if (tryflag) { // SP not already part of a track
+ Int_t newtrack = 1;
+ for (Int_t isp2 = isp1 + 1; isp2 < fNSp; isp2++) {
+ THcSpacePoint* sp2 = fSp[isp2];
+ if (sp1->fNChamber != sp2->fNChamber && sp1->GetSetStubFlag() && sp2->GetSetStubFlag()) {
+ Double_t* spstub1 = sp1->GetStubP();
+ Double_t* spstub2 = sp2->GetStubP();
+ Double_t dposx = spstub1[0] - spstub2[0];
+ Double_t dposy;
+ if (fProjectToChamber) { // From SOS s_link_stubs
+ // Since single chamber resolution is ~50mr, and the maximum y`
+ // angle is about 30mr, use differenece between y AT CHAMBERS, rather
+ // than at focal plane. (Project back to chamber, to take out y' uncertainty)
+ // (Should this be done for SHMS and HMS too?)
+ Double_t y1 = spstub1[1] + fChambers[sp1->fNChamber]->GetZPos() * spstub1[3];
+ Double_t y2 = spstub2[1] + fChambers[sp2->fNChamber]->GetZPos() * spstub2[3];
+ dposy = y1 - y2;
+ } else {
+ dposy = spstub1[1] - spstub2[1];
+ }
+ Double_t dposxp = spstub1[2] - spstub2[2];
+ Double_t dposyp = spstub1[3] - spstub2[3];
+
+ // What is the point of saving these stubmin values. They
+ // Don't seem to be used anywhere except that they can be
+ // printed out if hbypass_track_eff_files is zero.
+ if (TMath::Abs(dposx) < TMath::Abs(stubminx))
+ stubminx = dposx;
+ if (TMath::Abs(dposy) < TMath::Abs(stubminy))
+ stubminy = dposy;
+ if (TMath::Abs(dposxp) < TMath::Abs(stubminxp))
+ stubminxp = dposxp;
+ if (TMath::Abs(dposyp) < TMath::Abs(stubminyp))
+ stubminyp = dposyp;
+
+ // if hbypass_track_eff_files == 0 then
+ // Print out each stubminX that is less that its criterion
+
+ if ((TMath::Abs(dposx) < fXtTrCriterion) && (TMath::Abs(dposy) < fYtTrCriterion) &&
+ (TMath::Abs(dposxp) < fXptTrCriterion) && (TMath::Abs(dposyp) < fYptTrCriterion)) {
+ if (newtrack) {
+ assert(sptracks == 0);
+ fStubTest = 1;
+ // stubtest=1; Used in h_track_tests.f
+ // Make a new track if there are not to many
+ if (fNDCTracks < fNTracksMaxFP) {
+ sptracks = 0; // Number of tracks with this seed
+ stub_tracks[sptracks++] = fNDCTracks;
+ THcDCTrack* theDCTrack = new ((*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
+ theDCTrack->AddSpacePoint(sp1);
+ theDCTrack->AddSpacePoint(sp2);
+ if (sp1->fNChamber == 1)
+ theDCTrack->SetSp1_ID(sp1->fNChamber_spnum);
+ if (sp1->fNChamber == 2)
+ theDCTrack->SetSp2_ID(sp1->fNChamber_spnum);
+ if (sp2->fNChamber == 1)
+ theDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
+ if (sp2->fNChamber == 2)
+ theDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
+ newtrack = 0; // Make no more tracks in this loop
+ // (But could replace a SP?)
+ } else {
if (fHMSStyleChambers) {
- fNDCTracks=0;
- return;
+ fNDCTracks = 0;
+ return;
}
- }
- } else {
- // Check if there is another space point in the same chamber
- for(Int_t itrack=0;itrack<sptracks;itrack++) {
- Int_t track=stub_tracks[itrack];
- THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(track));
- Int_t spoint=-1;
- Int_t duppoint=0;
- for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
- // isp is index of space points in theDCTrack
- if(sp2->fNChamber ==
- theDCTrack->GetSpacePoint(isp)->fNChamber) {
- spoint=isp;
- }
- if(sp2==theDCTrack->GetSpacePoint(isp)) {
- duppoint=1;
- }
- } // End loop over sp in tracks with isp1
- // If there is no other space point in this chamber
- // add this space point to current track(2)
- if(!duppoint) {
- if(spoint<0) {
- theDCTrack->AddSpacePoint(sp2);
- if (sp2->fNChamber==1) theDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
- if (sp2->fNChamber==2) theDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
- } else {
- // If there is another point in the same chamber
- // in this track create a new track with all the
- // same space points except spoint
- if(fNDCTracks < MAXTRACKS) {
- stub_tracks[sptracks++] = fNDCTracks;
- THcDCTrack *newDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
- for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
- if(isp!=spoint) {
- newDCTrack->AddSpacePoint(theDCTrack->GetSpacePoint(isp));
- if (theDCTrack->GetSpacePoint(isp)->fNChamber==1) newDCTrack->SetSp1_ID(theDCTrack->GetSpacePoint(isp)->fNChamber_spnum);
- if (theDCTrack->GetSpacePoint(isp)->fNChamber==2) newDCTrack->SetSp2_ID(theDCTrack->GetSpacePoint(isp)->fNChamber_spnum);
- } else {
- newDCTrack->AddSpacePoint(sp2);
- if (sp2->fNChamber==1) newDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
- if (sp2->fNChamber==2) newDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
- } // End check for dup on copy
- } // End copy of track
- } else {
- if (fHMSStyleChambers) {
- if (fdebuglinkstubs) cout << "EPIC FAIL 2: Too many tracks found in THcDC::LinkStubs maxtracks = " << MAXTRACKS << endl;
- fNDCTracks=0;
- return; // Max # of allowed tracks
+ }
+ } else {
+ // Check if there is another space point in the same chamber
+ for (Int_t itrack = 0; itrack < sptracks; itrack++) {
+ Int_t track = stub_tracks[itrack];
+ THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(track));
+ Int_t spoint = -1;
+ Int_t duppoint = 0;
+ for (Int_t isp = 0; isp < theDCTrack->GetNSpacePoints(); isp++) {
+ // isp is index of space points in theDCTrack
+ if (sp2->fNChamber == theDCTrack->GetSpacePoint(isp)->fNChamber) {
+ spoint = isp;
+ }
+ if (sp2 == theDCTrack->GetSpacePoint(isp)) {
+ duppoint = 1;
+ }
+ } // End loop over sp in tracks with isp1
+ // If there is no other space point in this chamber
+ // add this space point to current track(2)
+ if (!duppoint) {
+ if (spoint < 0) {
+ theDCTrack->AddSpacePoint(sp2);
+ if (sp2->fNChamber == 1)
+ theDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
+ if (sp2->fNChamber == 2)
+ theDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
+ } else {
+ // If there is another point in the same chamber
+ // in this track create a new track with all the
+ // same space points except spoint
+ if (fNDCTracks < MAXTRACKS) {
+ stub_tracks[sptracks++] = fNDCTracks;
+ THcDCTrack* newDCTrack =
+ new ((*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
+ for (Int_t isp = 0; isp < theDCTrack->GetNSpacePoints(); isp++) {
+ if (isp != spoint) {
+ newDCTrack->AddSpacePoint(theDCTrack->GetSpacePoint(isp));
+ if (theDCTrack->GetSpacePoint(isp)->fNChamber == 1)
+ newDCTrack->SetSp1_ID(
+ theDCTrack->GetSpacePoint(isp)->fNChamber_spnum);
+ if (theDCTrack->GetSpacePoint(isp)->fNChamber == 2)
+ newDCTrack->SetSp2_ID(
+ theDCTrack->GetSpacePoint(isp)->fNChamber_spnum);
+ } else {
+ newDCTrack->AddSpacePoint(sp2);
+ if (sp2->fNChamber == 1)
+ newDCTrack->SetSp1_ID(sp2->fNChamber_spnum);
+ if (sp2->fNChamber == 2)
+ newDCTrack->SetSp2_ID(sp2->fNChamber_spnum);
+ } // End check for dup on copy
+ } // End copy of track
+ } else {
+ if (fHMSStyleChambers) {
+ if (fdebuglinkstubs)
+ cout << "EPIC FAIL 2: Too many tracks found in THcDC::LinkStubs "
+ "maxtracks = "
+ << MAXTRACKS << endl;
+ fNDCTracks = 0;
+ return; // Max # of allowed tracks
}
- }
- } // end if on same chamber
- } // end if on duplicate point
- } // end for over tracks with isp1
- } // else newtrack
- } // criterion
- } // end test on same chamber
- } // end isp2 loop over new space points
- } // end test on tryflag
- } // end isp1 outer loop over space points
+ }
+ } // end if on same chamber
+ } // end if on duplicate point
+ } // end for over tracks with isp1
+ } // else newtrack
+ } // criterion
+ } // end test on same chamber
+ } // end isp2 loop over new space points
+ } // end test on tryflag
+ } // end isp1 outer loop over space points
+ //
//
- //
} else { // Make track out of each single space point
- for(Int_t isp=0;isp<fNSp;isp++) {
- if(fNDCTracks<MAXTRACKS) {
- // Need some constructed t thingy
+ for (Int_t isp = 0; isp < fNSp; isp++) {
+ 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;
- fNDCTracks=0;
- // Do something here to fail this event
- return; // Max # of allowed tracks
+ if (fdebuglinkstubs)
+ cout << "EPIC FAIL 3: Too many tracks found in THcDC::LinkStubs" << endl;
+ fNDCTracks = 0;
+ // Do something here to fail this event
+ return; // Max # of allowed tracks
}
}
}
///
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("%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");
}
@@ -969,123 +1023,116 @@ void THcDC::LinkStubs()
}
//_____________________________________________________________________________
-void THcDC::TrackFit()
-{
+void THcDC::TrackFit() {
/**
Primary track fitting routine
*/
// Number of ray parameters in focal plane.
- const Int_t raycoeffmap[]={4,5,2,3};
+ const Int_t raycoeffmap[] = {4, 5, 2, 3};
Double_t dummychi2 = 1.0E4;
- for(UInt_t itrack=0;itrack<fNDCTracks;itrack++) {
+ for (UInt_t itrack = 0; itrack < fNDCTracks; itrack++) {
// Double_t chi2 = dummychi2;
// Int_t htrack_fit_num = itrack;
- THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
+ THcDCTrack* theDCTrack = static_cast<THcDCTrack*>(fDCTracks->At(itrack));
Double_t coords[theDCTrack->GetNHits()];
- Int_t planes[theDCTrack->GetNHits()];
- for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
- 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();
- }
+ Int_t planes[theDCTrack->GetNHits()];
+ for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {
+ 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();
+ }
} 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();
+ }
}
-
- } //end loop over hits
+ } // end loop over hits
theDCTrack->SetNFree(theDCTrack->GetNHits() - NUM_FPRAY);
Double_t chi2 = dummychi2;
- if(theDCTrack->GetNFree() > 0) {
+ 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);
-
- 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
+ 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;
+ // }
+
+ } // 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);
-
-
- AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
- fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
- pow(hit->GetWireSigma(),2);
-
- } //end ihit 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);
+
+ AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]] *
+ fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]] /
+ pow(hit->GetWireSigma(), 2);
+
+ } // end ihit loop
+ }
+ }
}
// Solve 4x4 equations
TVectorD dray(NUM_FPRAY);
// Should check that it is invertable
AA.Invert();
- dray = AA*TT;
- // cout << "DRAY: " << dray[0] << " "<< dray[1] << " "<< dray[2] << " "<< dray[3] << " " << endl;
- // if(bad_determinant) {
+ dray = AA * TT;
+ // cout << "DRAY: " << dray[0] << " "<< dray[1] << " "<< dray[2] << " "<< dray[3] << " "
+ // << endl; if(bad_determinant) {
// dray[0] = dray[1] = 10000.; dray[2] = dray[3] = 2.0;
// }
// Calculate hit coordinate for each plane for chi2 and efficiency
// calculations
// 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);
+ 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);
}
// Compute Chi2 and residuals
chi2 = 0.0;
- for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
-
- THcDCHit* hit=theDCTrack->GetHit(ihit);
-
+ for (Int_t ihit = 0; ihit < theDCTrack->GetNHits(); ihit++) {
- Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
- theDCTrack->SetResidual(planes[ihit], residual);
+ THcDCHit* hit = theDCTrack->GetHit(ihit);
- // double track_coord = theDCTrack->GetCoord(planes[ihit]);
-//cout<<planes[ihit]<<"\t"<<track_coord<<"\t"<<coords[ihit]<<"\t"<<residual<<endl;
- chi2 += pow(residual/hit->GetWireSigma(),2);
+ Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
+ theDCTrack->SetResidual(planes[ihit], residual);
+ // double track_coord = theDCTrack->GetCoord(planes[ihit]);
+ // cout<<planes[ihit]<<"\t"<<track_coord<<"\t"<<coords[ihit]<<"\t"<<residual<<endl;
+ chi2 += pow(residual / hit->GetWireSigma(), 2);
}
theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
@@ -1093,181 +1140,171 @@ void THcDC::TrackFit()
theDCTrack->SetChisq(chi2);
// 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++) {
-
-
- THcDCHit* hit=theDCTrack->GetHit(ihit);
-
- 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 ihit=0;ihit < theDCTrack->GetNHits();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);
-
- }
- }
- }
- }
- }
- //
- // 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];
- // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
- }
- Double_t residual = coords[ipl_hit] - coord;
- theDCTrack->SetResidualExclPlane(planes[ipl_hit], residual);
+ 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++) {
+
+ THcDCHit* hit = theDCTrack->GetHit(ihit);
+
+ 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 ihit = 0; ihit < theDCTrack->GetNHits(); 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);
+ }
+ }
+ }
+ }
+ }
+ //
+ // 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];
+ // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
+ }
+ 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
// Specific for two chambers. Can/should it be generalized?
- if(fSingleStub != 0) {
- if(fNDCTracks == 2) {
- THcDCTrack *theDCTrack1 = static_cast<THcDCTrack*>( fDCTracks->At(0));
- THcDCTrack *theDCTrack2 = static_cast<THcDCTrack*>( fDCTracks->At(1));
+ if (fSingleStub != 0) {
+ if (fNDCTracks == 2) {
+ THcDCTrack* theDCTrack1 = static_cast<THcDCTrack*>(fDCTracks->At(0));
+ THcDCTrack* theDCTrack2 = static_cast<THcDCTrack*>(fDCTracks->At(1));
// Int_t itrack=0;
- Int_t ihit=0;
- THcDCHit* hit=theDCTrack1->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- Int_t chamber=fNChamber[plane];
- if(chamber==1) {
- // itrack=1;
- hit=theDCTrack2->GetHit(ihit);
- plane=hit->GetPlaneNum()-1;
- chamber=fNChamber[plane];
- if(chamber==2) {
- Double_t ray1[4];
- Double_t ray2[4];
- theDCTrack1->GetRay(ray1);
- theDCTrack2->GetRay(ray2);
- // itrack = 1;
- // Loop over hits in second chamber
- for(Int_t ihit=0;ihit < theDCTrack2->GetNHits();ihit++) {
- // Calculate residual in second chamber from first chamber track
- THcDCHit* hit=theDCTrack2->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- Double_t pos = DpsiFun(ray1,plane);
- Double_t coord;
- if(fFixLR) {
- if(fFixPropagationCorrection) {
- coord = hit->GetPos()
- + theDCTrack2->GetHitLR(ihit)*theDCTrack2->GetHitDist(ihit);
- } else {
- coord = hit->GetPos()
- + theDCTrack2->GetHitLR(ihit)*hit->GetDist();
- }
- } else {
- if(fFixPropagationCorrection) {
- coord = hit->GetPos()
- + hit->GetLR()*theDCTrack2->GetHitDist(ihit);
- } else {
- coord = hit->GetCoord();
- }
- }
- theDCTrack1->SetDoubleResidual(plane,coord - pos);
- // hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
- }
- // itrack=0;
- // Loop over hits in first chamber
- for(Int_t ihit=0;ihit < theDCTrack1->GetNHits();ihit++) {
- // Calculate residual in first chamber from second chamber track
- THcDCHit* hit=theDCTrack1->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- Double_t pos = DpsiFun(ray1,plane);
- Double_t coord;
- if(fFixLR) {
- if(fFixPropagationCorrection) {
- coord = hit->GetPos()
- + theDCTrack1->GetHitLR(ihit)*theDCTrack1->GetHitDist(ihit);
- } else {
- coord = hit->GetPos()
- + theDCTrack1->GetHitLR(ihit)*hit->GetDist();
- }
- } else {
- if(fFixPropagationCorrection) {
- coord = hit->GetPos()
- + hit->GetLR()*theDCTrack1->GetHitDist(ihit);
- } else {
- coord = hit->GetCoord();
- }
- }
- theDCTrack2->SetDoubleResidual(plane,coord - pos);
- // hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
- }
- }
+ Int_t ihit = 0;
+ THcDCHit* hit = theDCTrack1->GetHit(ihit);
+ Int_t plane = hit->GetPlaneNum() - 1;
+ Int_t chamber = fNChamber[plane];
+ if (chamber == 1) {
+ // itrack=1;
+ hit = theDCTrack2->GetHit(ihit);
+ plane = hit->GetPlaneNum() - 1;
+ chamber = fNChamber[plane];
+ if (chamber == 2) {
+ Double_t ray1[4];
+ Double_t ray2[4];
+ theDCTrack1->GetRay(ray1);
+ theDCTrack2->GetRay(ray2);
+ // itrack = 1;
+ // Loop over hits in second chamber
+ for (Int_t ihit = 0; ihit < theDCTrack2->GetNHits(); ihit++) {
+ // Calculate residual in second chamber from first chamber track
+ THcDCHit* hit = theDCTrack2->GetHit(ihit);
+ Int_t plane = hit->GetPlaneNum() - 1;
+ Double_t pos = DpsiFun(ray1, plane);
+ Double_t coord;
+ if (fFixLR) {
+ if (fFixPropagationCorrection) {
+ coord = hit->GetPos() + theDCTrack2->GetHitLR(ihit) * theDCTrack2->GetHitDist(ihit);
+ } else {
+ coord = hit->GetPos() + theDCTrack2->GetHitLR(ihit) * hit->GetDist();
+ }
+ } else {
+ if (fFixPropagationCorrection) {
+ coord = hit->GetPos() + hit->GetLR() * theDCTrack2->GetHitDist(ihit);
+ } else {
+ coord = hit->GetCoord();
+ }
+ }
+ theDCTrack1->SetDoubleResidual(plane, coord - pos);
+ // hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
+ }
+ // itrack=0;
+ // Loop over hits in first chamber
+ for (Int_t ihit = 0; ihit < theDCTrack1->GetNHits(); ihit++) {
+ // Calculate residual in first chamber from second chamber track
+ THcDCHit* hit = theDCTrack1->GetHit(ihit);
+ Int_t plane = hit->GetPlaneNum() - 1;
+ Double_t pos = DpsiFun(ray1, plane);
+ Double_t coord;
+ if (fFixLR) {
+ if (fFixPropagationCorrection) {
+ coord = hit->GetPos() + theDCTrack1->GetHitLR(ihit) * theDCTrack1->GetHitDist(ihit);
+ } else {
+ coord = hit->GetPos() + theDCTrack1->GetHitLR(ihit) * hit->GetDist();
+ }
+ } else {
+ if (fFixPropagationCorrection) {
+ coord = hit->GetPos() + hit->GetLR() * theDCTrack1->GetHitDist(ihit);
+ } else {
+ coord = hit->GetCoord();
+ }
+ }
+ theDCTrack2->SetDoubleResidual(plane, coord - pos);
+ // hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
+ }
+ }
}
}
}
//
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());
+ 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);
- } else {
- coords_temp = hit->GetCoord();
- }
- }
- printf("%-5d %15.7e %15.7e %15.7e \n",plane+1,coords_temp,theDCTrack->GetCoord(plane),theDCTrack->GetResidual(plane));
+ 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);
+ } else {
+ coords_temp = hit->GetCoord();
+ }
+ }
+ printf("%-5d %15.7e %15.7e %15.7e \n", plane + 1, coords_temp, theDCTrack->GetCoord(plane),
+ theDCTrack->GetResidual(plane));
}
}
}
@@ -1276,8 +1313,7 @@ void THcDC::TrackFit()
}
//
//
-Double_t THcDC::DpsiFun(Double_t ray[4], Int_t plane)
-{
+Double_t THcDC::DpsiFun(Double_t ray[4], Int_t plane) {
/**
this function calculates the psi coordinate of the intersection
of a ray (defined by ray) with a hms wire chamber plane. the geometry
@@ -1296,74 +1332,73 @@ Double_t THcDC::DpsiFun(Double_t ray[4], Int_t plane)
ray(4) = tan(yp)
*/
- Double_t infinity = 1.0E+20;
- Double_t cinfinity = 1/infinity;
- Double_t DpsiFun =
- ray[2]*ray[1]*fPlaneCoeffs[plane][0] +
- ray[3]*ray[0]*fPlaneCoeffs[plane][1] +
- ray[2]*fPlaneCoeffs[plane][2] +
- ray[3]*fPlaneCoeffs[plane][3] +
- ray[0]*fPlaneCoeffs[plane][4] +
- ray[1]*fPlaneCoeffs[plane][5];
- Double_t denom = ray[2]*fPlaneCoeffs[plane][6]
- + ray[3]*fPlaneCoeffs[plane][7]
- + fPlaneCoeffs[plane][8];
- if(TMath::Abs(denom) < cinfinity) {
+ Double_t infinity = 1.0E+20;
+ Double_t cinfinity = 1 / infinity;
+ Double_t DpsiFun = ray[2] * ray[1] * fPlaneCoeffs[plane][0] +
+ ray[3] * ray[0] * fPlaneCoeffs[plane][1] + ray[2] * fPlaneCoeffs[plane][2] +
+ ray[3] * fPlaneCoeffs[plane][3] + ray[0] * fPlaneCoeffs[plane][4] +
+ ray[1] * fPlaneCoeffs[plane][5];
+ Double_t denom =
+ ray[2] * fPlaneCoeffs[plane][6] + ray[3] * fPlaneCoeffs[plane][7] + fPlaneCoeffs[plane][8];
+ if (TMath::Abs(denom) < cinfinity) {
DpsiFun = infinity;
} else {
- DpsiFun = DpsiFun/denom;
+ DpsiFun = DpsiFun / denom;
}
- return(DpsiFun);
+ return (DpsiFun);
}
//_____________________________________________________________________________
-Int_t THcDC::End(THaRunBase* run)
-{
+Int_t THcDC::End(THaRunBase* run) {
// EffCalc();
MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
return 0;
}
//_____________________________________________________________________________
-void THcDC::EffInit()
-{
+void THcDC::EffInit() {
/**
Create, and initialize counters used to calculate
efficiencies. Register the counters in gHcParms so that the
variables can be used in end of run reports.
*/
- delete [] fNChamHits; fNChamHits = new Int_t [fNChambers];
- delete [] fPlaneEvents; fPlaneEvents = new Int_t [fNPlanes];
+ delete[] fNChamHits;
+ fNChamHits = new Int_t[fNChambers];
+ delete[] fPlaneEvents;
+ fPlaneEvents = new Int_t[fNPlanes];
fTotEvents = 0;
- for(UInt_t i=0;i<fNChambers;i++) {
+ for (UInt_t i = 0; i < fNChambers; i++) {
fNChamHits[i] = 0;
}
- for(Int_t i=0;i<fNPlanes;i++) {
+ for (Int_t i = 0; i < fNPlanes; i++) {
fPlaneEvents[i] = 0;
}
- gHcParms->Define(Form("%sdc_tot_events",fPrefix),"Total DC Events",fTotEvents);
- gHcParms->Define(Form("%sdc_cham_hits[%d]",fPrefix,fNChambers),"N events with hits per chamber",*fNChamHits);
- gHcParms->Define(Form("%sdc_events[%d]",fPrefix,fNPlanes),"N events with hits per plane",*fPlaneEvents);
+ gHcParms->Define(Form("%sdc_tot_events", fPrefix), "Total DC Events", fTotEvents);
+ gHcParms->Define(Form("%sdc_cham_hits[%d]", fPrefix, fNChambers),
+ "N events with hits per chamber", *fNChamHits);
+ gHcParms->Define(Form("%sdc_events[%d]", fPrefix, fNPlanes), "N events with hits per plane",
+ *fPlaneEvents);
}
//_____________________________________________________________________________
-void THcDC::Eff()
-{
+void THcDC::Eff() {
/**
Accumulate statistics for efficiency calculations
*/
fTotEvents++;
- for(UInt_t i=0;i<fNChambers;i++) {
- if(fChambers[i]->GetNHits()>0) fNChamHits[i]++;
+ for (UInt_t i = 0; i < fNChambers; i++) {
+ if (fChambers[i]->GetNHits() > 0)
+ fNChamHits[i]++;
}
- for(Int_t i=0;i<fNPlanes;i++) {
- if(fPlanes[i]->GetNHits() > 0) fPlaneEvents[i]++;
+ for (Int_t i = 0; i < fNPlanes; i++) {
+ if (fPlanes[i]->GetNHits() > 0)
+ fPlaneEvents[i]++;
}
return;
}
ClassImp(THcDC)
-////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
--
GitLab