diff --git a/src/libraries/AMPTOOLS_AMPS/DeltaAngles.cc b/src/libraries/AMPTOOLS_AMPS/DeltaAngles.cc new file mode 100644 index 000000000..1b1a67463 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/DeltaAngles.cc @@ -0,0 +1,180 @@ + +#include +#include +#include +#include +#include + +#include "TLorentzVector.h" +#include "TLorentzRotation.h" + +#include "IUAmpTools/Kinematics.h" +#include "AMPTOOLS_AMPS/DeltaAngles.h" +#include "AMPTOOLS_AMPS/clebschGordan.h" +#include "AMPTOOLS_AMPS/wignerD.h" +#include "AMPTOOLS_AMPS/decayAngles.h" + +DeltaAngles::DeltaAngles( const vector< string >& args ) : +UserAmplitude< DeltaAngles >( args ) +{ + assert( args.size() == 13 || args.size() == 11 ); + + rho011 = AmpParameter( args[0] ); + rho031 = AmpParameter( args[1] ); + rho03m1 = AmpParameter( args[2] ); + + rho111 = AmpParameter( args[3] ); + rho133 = AmpParameter( args[4] ); + rho131 = AmpParameter( args[5] ); + rho13m1 = AmpParameter( args[6] ); + + rho231 = AmpParameter( args[7] ); + rho23m1 = AmpParameter( args[8] ); + + lowerVertex = args[9].c_str(); + upperVertex = args[10].c_str(); + + // need to register any free parameters so the framework knows about them + registerParameter( rho011 ); + registerParameter( rho031 ); + registerParameter( rho03m1 ); + + registerParameter( rho111 ); + registerParameter( rho133 ); + registerParameter( rho131 ); + registerParameter( rho13m1 ); + + registerParameter( rho231 ); + registerParameter( rho23m1 ); + + if(args.size() == 13){ + polAngle = atof(args[11].c_str() ); // azimuthal angle of the photon polarization vector in the lab. + polFraction = AmpParameter( args[12] ); // fraction of polarization (0-1) +// std::cout << "Fixed polarisation of " << polFraction << " and angle of " << polAngle << " degrees." << std::endl; + } + else + assert(0); +} + + +complex< GDouble > +DeltaAngles::calcAmplitude( GDouble** pKin, GDouble* userVars ) const { + + GDouble cosTheta = userVars[kCosTheta]; + GDouble sinSqTheta = userVars[kSinSqTheta]; + GDouble sin2Theta = userVars[kSin2Theta]; + + GDouble cosPhi = userVars[kCosPhi]; + GDouble cos2Phi = userVars[kCos2Phi]; + GDouble sinPhi = userVars[kSinPhi]; + GDouble sin2Phi = userVars[kSin2Phi]; + + GDouble sin2BigPhi = userVars[kSin2BigPhi]; + GDouble cos2BigPhi = userVars[kCos2BigPhi]; + GDouble Pgamma = userVars[kPgamma]; + GDouble sqrt3 = TMath::Sqrt(3); + + // SDMEs for 3/2- -> 1/2+ + 0- (doi.org/10.1103/PhysRevC.96.025208) + GDouble W = 3.*(0.5 - rho011)*sinSqTheta + rho011*(1.+3.*cosTheta*cosTheta) - 2.*sqrt3*rho031*cosPhi*sin2Theta - 2.*sqrt3*rho03m1*cos2Phi*sinSqTheta; + + W -= Pgamma*cos2BigPhi * (3.*rho133*sinSqTheta + rho111*(1.+3.*cosTheta*cosTheta) - 2.*sqrt3*rho131*cosPhi*sin2Theta - 2.*sqrt3*rho13m1*cos2Phi*sinSqTheta); + + W -= Pgamma*sin2BigPhi * (2.*sqrt3*rho231*sinPhi*sin2Theta + 2.*sqrt3*rho23m1*sin2Phi*sinSqTheta); + + W *= 1./(4.*PI); + +// return W; + return complex< GDouble > ( sqrt(fabs(W)) ); +} + +void +DeltaAngles::calcUserVars( GDouble** pKin, GDouble* userVars ) const { + + TLorentzVector target ( 0, 0, 0, 0.9382720813 ); + TLorentzVector beam ( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); +// TLorentzVector p1, p2, p3, ptot, ptemp; //p1 and p2 from decaying lower vertex, p2 used to calculate angles for SDME calculation, p3 = upper vertex resonance (b1 in this case) + TLorentzVector p1, p2, pDelta; //p1 and p2 from decaying lower vertex, p2 used to calculate angles for SDME calculation, p3 = upper vertex resonance (b1 in this case) + + string lv1; lv1 += lowerVertex[0]; + string lv2; lv2 += lowerVertex[1]; + + int index1 = atoi( lv1.c_str() ); + int index2 = atoi( lv2.c_str() ); + + p1.SetPxPyPzE( pKin[index1][1], pKin[index1][2], pKin[index1][3], pKin[index1][0] ); + p2.SetPxPyPzE( pKin[index2][1], pKin[index2][2], pKin[index2][3], pKin[index2][0] ); + + pDelta = p1 + p2; + + vector< double > thetaPhi = getOneStepAngles( pDelta, p1, beam, target, 2, false ); + + userVars[kCosTheta] = TMath::Cos( thetaPhi[0] ); + userVars[kSinSqTheta] = TMath::Sin( thetaPhi[0] ) * TMath::Sin( thetaPhi[0] ); + userVars[kSin2Theta] = TMath::Sin( 2.*thetaPhi[0] ); + userVars[kCosPhi] = cos(thetaPhi[1]); + userVars[kCos2Phi] = cos(2*thetaPhi[1]); + userVars[kSinPhi] = sin(thetaPhi[1]); + userVars[kSin2Phi] = sin(2*thetaPhi[1]); + + double phiProd = getPhiProd( polAngle, pDelta, beam, target, 2, false ); + userVars[kCos2BigPhi] = cos(2*phiProd); + userVars[kSin2BigPhi] = sin(2*phiProd); + + +/* + //p3 is sum of all particles in upper vertex + for( unsigned int i = 0; i < upperVertex.size(); ++i ){ + string num; num += upperVertex[i]; + int index = atoi(num.c_str()); + ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], pKin[index][3], pKin[index][0] ); + p3 += ptemp; + ptot += ptemp; + } + + TLorentzVector lowerVertexResonance = p1 + p2; + TLorentzRotation lowerVertexBoost( -lowerVertexResonance.BoostVector() ); + + TLorentzVector target_lowerVertexRF = lowerVertexBoost * target; + TLorentzVector beam_lowerVertexRF = lowerVertexBoost * beam; + TLorentzVector upperVertexResonance_lowerVertexRF = lowerVertexBoost * p3; + TLorentzVector p2_lowerVertexRF = lowerVertexBoost * p2; + + // normal to the production plane + TVector3 y = (target_lowerVertexRF.Vect().Unit().Cross(upperVertexResonance_lowerVertexRF.Vect().Unit())).Unit(); + // choose Gottfried-Jackson frame: z-axis along -target direction in baryon rest frame + TVector3 z = target_lowerVertexRF.Vect().Unit(); + TVector3 x = y.Cross(z).Unit(); + + TVector3 angles( (p2_lowerVertexRF.Vect()).Dot(x), + (p2_lowerVertexRF.Vect()).Dot(y), + (p2_lowerVertexRF.Vect()).Dot(z) ); + + userVars[kCosTheta] = angles.CosTheta(); + userVars[kSinSqTheta] = sin(angles.Theta())*sin(angles.Theta()); +// userVars[kCosSqTheta] = cos(angles.Theta())*cos(angles.Theta()); + userVars[kSin2Theta] = sin(2.*angles.Theta()); + userVars[kPhi] = angles.Phi(); + + TVector3 eps(cos(polAngle*TMath::DegToRad()), sin(polAngle*TMath::DegToRad()), 0.0); // beam polarization vector in lab + userVars[kBigPhi] = atan2(y.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(y))); +// Phi = Phi > 0? Phi : Phi + 3.14159; +*/ + // polarization BeamProperties + GDouble Pgamma = polFraction; + + if(polAngle == -1) + Pgamma = 0.; + + userVars[kPgamma] = Pgamma; +} + + +#ifdef GPU_ACCELERATION +void +DeltaAngles::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const { + + GPUDeltaAngles_exec( dimGrid, dimBlock, GPU_AMP_ARGS, + rho011, rho031, rho03m1, rho111, rho133, rho131, rho13m1, rho231, rho23m1, polAngle ); +} + +#endif // GPU_ACCELERATION diff --git a/src/libraries/AMPTOOLS_AMPS/DeltaAngles.h b/src/libraries/AMPTOOLS_AMPS/DeltaAngles.h new file mode 100644 index 000000000..56b0e48fc --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/DeltaAngles.h @@ -0,0 +1,79 @@ +#if !defined(DELTAANGLES) +#define DELTAANGLES + +#include "IUAmpTools/Amplitude.h" +#include "IUAmpTools/UserAmplitude.h" +#include "IUAmpTools/AmpParameter.h" +#include "GPUManager/GPUCustomTypes.h" + +#include "TH1D.h" + +#include +#include +#include + + + +using std::complex; +using namespace std; + +#ifdef GPU_ACCELERATION +void GPUDeltaAngles_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble rho011, GDouble rho031, GDouble rho03m1, + GDouble rho111, GDouble rho133, GDouble rho131, GDouble rho13m1, + GDouble rho231, GDouble rho23m1, GDouble polAngle ); +#endif // GPU_ACCELERATION + + +class Kinematics; + +class DeltaAngles : public UserAmplitude< DeltaAngles > +{ + +public: + + DeltaAngles() : UserAmplitude< DeltaAngles >() { }; + DeltaAngles( const vector< string >& args ); + + enum UserVars { kPgamma = 0, kCosTheta, kSinSqTheta, kSin2Theta, + kCosPhi, kCos2Phi, kSinPhi, kSin2Phi, kCos2BigPhi, kSin2BigPhi, kNumUserVars }; + unsigned int numUserVars() const { return kNumUserVars; } + + string name() const { return "DeltaAngles"; } + + complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const; + void calcUserVars( GDouble** pKin, GDouble* userVars ) const; + + bool needsUserVarsOnly() const { return true; } + +#ifdef GPU_ACCELERATION + + void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const; + bool isGPUEnabled() const { return true; } + +#endif // GPU_ACCELERATION + +private: + + AmpParameter rho011; + AmpParameter rho031; + AmpParameter rho03m1; + + AmpParameter rho111; + AmpParameter rho133; + AmpParameter rho131; + AmpParameter rho13m1; + + AmpParameter rho231; + AmpParameter rho23m1; + + GDouble polFraction=0.; + GDouble polAngle=-1; + TH1D *polFrac_vs_E; + + string lowerVertex; + string upperVertex; + +}; + +#endif diff --git a/src/libraries/AMPTOOLS_AMPS/GPUDeltaAngles_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPUDeltaAngles_kernel.cu new file mode 100644 index 000000000..fcdc18c8f --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/GPUDeltaAngles_kernel.cu @@ -0,0 +1,58 @@ +#include + +#include "GPUManager/GPUCustomTypes.h" +#include "GPUManager/CUDA-Complex.cuh" + +#include "AMPTOOLS_AMPS/breakupMomentum.cuh" +#include "AMPTOOLS_AMPS/barrierFactor.cuh" + +#include "GPUUtils/wignerD.cuh" +#include "GPUUtils/clebsch.cuh" + +/////////////////////////////////////////////////////////////////////////////// +__global__ void +GPUDeltaAngles_kernel( GPU_AMP_PROTO, GDouble rho011, GDouble rho031, GDouble rho03m1, + GDouble rho111, GDouble rho133, GDouble rho131, GDouble rho13m1, + GDouble rho231, GDouble rho23m1, GDouble polAngle ) +{ + int iEvent = GPU_THIS_EVENT; + + GDouble Pgamma = GPU_UVARS(0); + GDouble cosTheta = GPU_UVARS(1); + GDouble sinSqTheta = GPU_UVARS(2); + GDouble sin2Theta = GPU_UVARS(3); + GDouble cosPhi = GPU_UVARS(4); + GDouble cos2Phi = GPU_UVARS(5); + GDouble sinPhi = GPU_UVARS(6); + GDouble sin2Phi = GPU_UVARS(7); + GDouble cos2BigPhi = GPU_UVARS(8); + GDouble sin2BigPhi = GPU_UVARS(9); + GDouble sqrt3 = sqrt(3.0); + + GDouble W = 3.*(0.5 - rho011)*sinSqTheta + rho011*(1.+3.*cosTheta*cosTheta) - 2.*sqrt3*rho031*cosPhi*sin2Theta - 2.*sqrt3*rho03m1*cos2Phi*sinSqTheta; + + W -= Pgamma*cos2BigPhi * (3.*rho133*sinSqTheta + rho111*(1.+3.*cosTheta*cosTheta) - 2.*sqrt3*rho131*cosPhi*sin2Theta - 2.*sqrt3*rho13m1*cos2Phi*sinSqTheta); + + W -= Pgamma*sin2BigPhi * (2.*sqrt3*rho231*sinPhi*sin2Theta + 2.*sqrt3*rho23m1*sin2Phi*sinSqTheta); + + W *= 1.0/(4.0*PI); + + WCUComplex amp = { sqrt( fabs( W ) ), 0 }; + + pcDevAmp[iEvent] = amp; + +} + +void +GPUDeltaAngles_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble rho011, GDouble rho031, GDouble rho03m1, + GDouble rho111, GDouble rho133, GDouble rho131, GDouble rho13m1, + GDouble rho231, GDouble rho23m1, GDouble polAngle ) + +{ + + GPUDeltaAngles_kernel<<< dimGrid, dimBlock >>> + ( GPU_AMP_ARGS, rho011, rho031, rho03m1, rho111, rho133, rho131, rho13m1, rho231, rho23m1, polAngle ); + +} + diff --git a/src/libraries/AMPTOOLS_AMPS/GPULinear_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPULinear_kernel.cu new file mode 100644 index 000000000..1c37c6d33 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/GPULinear_kernel.cu @@ -0,0 +1,30 @@ +#include + +#include "GPUManager/GPUCustomTypes.h" +#include "GPUManager/CUDA-Complex.cuh" + + +__global__ void +GPULinear_kernel( GPU_AMP_PROTO, GDouble real_p0, GDouble real_p1, + GDouble imag_p0, GDouble imag_p1 ){ + + int iEvent = GPU_THIS_EVENT; + + GDouble mass = GPU_UVARS(0); + + GDouble real_tot = real_p0 + real_p1 * mass; + GDouble imag_tot = imag_p0 + imag_p1 * mass; + + WCUComplex amp = { real_tot, imag_tot }; + pcDevAmp[iEvent] = amp; +} + +void +GPULinear_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble real_p0, GDouble real_p1, + GDouble imag_p0, GDouble imag_p1 ){ + + GPULinear_kernel<<< dimGrid, dimBlock >>> + ( GPU_AMP_ARGS, real_p0, real_p1, imag_p0, imag_p1 ); +} + diff --git a/src/libraries/AMPTOOLS_AMPS/Linear.cc b/src/libraries/AMPTOOLS_AMPS/Linear.cc new file mode 100644 index 000000000..badc1c9c4 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/Linear.cc @@ -0,0 +1,87 @@ +#include +#include +#include +#include +#include + +#include "IUAmpTools/Kinematics.h" +#include "AMPTOOLS_AMPS/Linear.h" + +// this amplitude returns a function that is linear in +// invariant mass and has a constant phase. It takes +// the form of e^(i*phi) ( a + b*M ), where M is an +// invaraiant mass and phi, a, and b, are all real +// numbers -- the computation below is equivalent but +// done in terms of real and imaginary parts + +Linear::Linear( const vector< string >& args ) : +UserAmplitude< Linear >( args ) +{ + assert( args.size() == 5 ); + + m_daughters = pair< string, string >( args[0], args[1] ); + + m_real_p0 = AmpParameter( args[2] ); + m_real_p1 = AmpParameter( args[3] ); + m_imag_p0 = AmpParameter( args[4] ); + + registerParameter( m_real_p0 ); + registerParameter( m_real_p1 ); + registerParameter( m_imag_p0 ); +} + +complex< GDouble > +Linear::calcAmplitude( GDouble** pKin, GDouble* userVars ) const +{ + GDouble mass = userVars[kMass]; + + double real_tot = m_real_p0 + m_real_p1 * mass; + double imag_tot = m_imag_p0 + m_imag_p1 * mass; + + complex< GDouble > ans( real_tot, imag_tot ); + return ans; +} + +void Linear::calcUserVars( GDouble** pKin, GDouble* userVars ) const +{ + TLorentzVector P1, P2, Ptot, Ptemp; + + for( unsigned int i = 0; i < m_daughters.first.size(); ++i ){ + + string num; num += m_daughters.first[i]; + int index = atoi(num.c_str()); + Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], + pKin[index][3], pKin[index][0] ); + P1 += Ptemp; + Ptot += Ptemp; + } + + for( unsigned int i = 0; i < m_daughters.second.size(); ++i ){ + + string num; num += m_daughters.second[i]; + int index = atoi(num.c_str()); + Ptemp.SetPxPyPzE( pKin[index][1], pKin[index][2], + pKin[index][3], pKin[index][0] ); + P2 += Ptemp; + Ptot += Ptemp; + } + userVars[kMass] = Ptot.M(); +} + +void +Linear::updatePar( const AmpParameter& par ) +{ + m_imag_p1 = m_real_p1 * m_imag_p0 / m_real_p0; +} + + +#ifdef GPU_ACCELERATION +void +Linear::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const { + + GPULinear_exec( dimGrid, dimBlock, GPU_AMP_ARGS, + m_real_p0, m_real_p1, m_imag_p0, m_imag_p1 ); + +} +#endif //GPU_ACCELERATION + diff --git a/src/libraries/AMPTOOLS_AMPS/Linear.h b/src/libraries/AMPTOOLS_AMPS/Linear.h new file mode 100644 index 000000000..c7b4b6b16 --- /dev/null +++ b/src/libraries/AMPTOOLS_AMPS/Linear.h @@ -0,0 +1,67 @@ +#if !defined(LINEAR) +#define LINEAR + +#include "IUAmpTools/UserAmplitude.h" +#include "IUAmpTools/AmpParameter.h" +#include "GPUManager/GPUCustomTypes.h" + +#include +#include +#include +#include + +#ifdef GPU_ACCELERATION +void GPULinear_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO, + GDouble real_p0, GDouble real_p1, + GDouble imag_p0, GDouble imag_p1 ); +#endif // GPU_ACCELERATION + +using std::complex; +using namespace std; + +class Kinematics; + +class Linear : public UserAmplitude< Linear > +{ +public: + + Linear() : UserAmplitude< Linear >() { } + + Linear( const vector< string >& args ); + + ~Linear(){} + + string name() const { return "Linear"; } + + complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const; + + enum UserVars { kMass = 0, kNumUserVars }; + unsigned int numUserVars() const { return kNumUserVars; } + + void calcUserVars( GDouble** pKin, GDouble* userVars ) const; + + bool needsUserVarsOnly() const { return true; } + bool areUserVarsStatic() const { return false; } + + void updatePar( const AmpParameter& par ); + +#ifdef GPU_ACCELERATION + + void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const; + + bool isGPUEnabled() const { return true; } + +#endif // GPU_ACCELERATION + +private: + + pair< string, string > m_daughters; + AmpParameter m_real_p0; + AmpParameter m_real_p1; + AmpParameter m_imag_p0; + + GDouble m_imag_p1; + +}; + +#endif diff --git a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.cc b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.cc index 7e944b6bf..ecd8fd04b 100644 --- a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.cc +++ b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.cc @@ -11,6 +11,11 @@ #include "AMPTOOLS_DATAIO/FSRootDataReader.h" #include "TSystem.h" +#include "IUAmpTools/report.h" + +const char* FSRootDataReader::kModule = "FSRootDataReader"; + + using namespace std; // Constructor expects one of the following argument patterns: @@ -77,7 +82,7 @@ FSRootDataReader::FSRootDataReader( const vector< string >& args ) : fileexists.close(); m_inFile = new TFile( inFileName.c_str() ); if (!m_inFile || m_inFile->IsZombie()) { - cout << "FSRootDataReader WARNING: Cannot open file... " << inFileName << endl; + report( ERROR, kModule ) << "FSRootDataReader ERROR: Cannot open file... " << inFileName << endl; m_inFile = NULL; m_inTree = NULL; return; @@ -85,7 +90,7 @@ FSRootDataReader::FSRootDataReader( const vector< string >& args ) : m_inTree = static_cast( m_inFile->Get( inTreeName.c_str() ) ); if (!m_inTree) { - cout << "FSRootDataReader WARNING: Cannot open tree... " << inTreeName << endl; + report( ERROR, kModule ) << "FSRootDataReader ERROR: Cannot open tree... " << inTreeName << endl; m_inFile->Close(); delete m_inFile; m_inFile = NULL; @@ -97,17 +102,18 @@ FSRootDataReader::FSRootDataReader( const vector< string >& args ) : m_inTree->AddFriend(friendTreeName, friendFileName); } else{ - cout << "FSRootDataReader WARNING: Cannot find file... " << inFileName << endl; + report( ERROR, kModule ) << "FSRootDataReader ERROR: Cannot find file... " << inFileName << endl; m_inFile = NULL; m_inTree = NULL; return; } - cout << "Opening Tree: " << inFileName << " " << inTreeName << " (numParticles=" << m_numParticles << ")"; - if (fourMomentumPrefix != "") cout << " fourMomentumPrefix=" << fourMomentumPrefix; - if (friendFileName != "") cout << " friendFile=" << friendFileName << " friendTree=" << friendTreeName; - if (weightBranchName != "weight") cout << " weightBranch=" << weightBranchName; - cout << endl; + report( DEBUG, kModule ) << "Opening Tree: " << inFileName << " " << inTreeName << " (numParticles=" << m_numParticles << ")"; + if (fourMomentumPrefix != "") report( DEBUG, kModule ) << " fourMomentumPrefix=" << fourMomentumPrefix; + if (friendFileName != "") report( DEBUG, kModule ) << " friendFile=" << friendFileName << " friendTree=" << friendTreeName; + if (weightBranchName != "weight") report( DEBUG, kModule ) << " weightBranch=" << weightBranchName; + report( DEBUG, kModule ) << endl; + if (m_inTree){ TString sEnPB = fourMomentumPrefix+"EnPB"; TString sPxPB = fourMomentumPrefix+"PxPB"; diff --git a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.h b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.h index 4738f661e..14e69cb51 100644 --- a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.h +++ b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReader.h @@ -48,6 +48,7 @@ class FSRootDataReader : public UserDataReader< FSRootDataReader >{ double m_weight; + static const char* kModule; }; #endif diff --git a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.cc b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.cc index 745062862..180a6e9b6 100644 --- a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.cc +++ b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include // Added for multiset #include "TH1.h" #include "TFile.h" @@ -12,9 +13,9 @@ #include "AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.h" #include "TSystem.h" -using namespace std; +#include "IUAmpTools/report.h" -#include // Include for printing +const char* FSRootDataReaderBootstrap::kModule = "FSRootDataReaderBootstrap"; // Constructor expects one of the following argument patterns: // @@ -73,37 +74,47 @@ FSRootDataReaderBootstrap::FSRootDataReaderBootstrap( const vector< string >& ar fourMomentumPrefix = args[3]; randSeed = atoi(args[4].c_str()); } - else if (args.size() == 6) { - if (isInteger(args[4])) { - fourMomentumPrefix = args[3]; - randSeed = atoi(args[4].c_str()); - weightBranchName = args[5]; - } else { - friendFileName = args[3]; - friendTreeName = args[4]; - weightBranchName = args[5]; - } - } - else if (args.size() >= 7) { + else if (args.size() >= 6) { friendFileName = args[3]; friendTreeName = args[4]; weightBranchName = args[5]; + + if (args.size() == 7) { + if (isInteger(args[6])) { + randSeed = atoi(args[6].c_str()); + } else { + fourMomentumPrefix = args[6]; + } + } else if (args.size() == 8) { fourMomentumPrefix = args[6]; - } - else if (args.size() == 8) { randSeed = atoi(args[7].c_str()); + } } // Initialize random number generator with provided seed m_randGenerator = new TRandom3(randSeed); - cout << "******************** WARNING ***********************" << endl; - cout << "* You are using the boostrap data reader, which *" << endl; - cout << "* should only be used for evaluating errors. *" << endl; - cout << "* The results with different seeds will be random *" << endl; - cout << "* due to random oversampling of the input file. *" << endl; - cout << "****************************************************" << endl; - cout << endl; + // Compact printing of parsed arguments + report( DEBUG, kModule ) << "FSRootDataReaderBootstrap initialized with:\n" + << " inFileName: " << inFileName << "\n" + << " inTreeName: " << inTreeName << "\n" + << " numParticles: " << m_numParticles << "\n" + << " friendFileName: " << (friendFileName.Length() == 0 ? "N/A" : friendFileName.Data()) << "\n" + << " friendTreeName: " << (friendTreeName.Length() == 0 ? "N/A" : friendTreeName.Data()) << "\n" + << " weightBranchName: " << (weightBranchName.Length() == 0 ? "N/A" : weightBranchName.Data()) << "\n" + << " fourMomentumPrefix:" << (fourMomentumPrefix.Length() == 0 ? "N/A" : fourMomentumPrefix.Data()) << "\n" + << " randSeed: " << randSeed << "\n" + << std::endl; + + + report( NOTICE, kModule ) << "******************** NOTICE ************************" << endl; + report( NOTICE, kModule ) << "* You are using the boostrap data reader, which *" << endl; + report( NOTICE, kModule ) << "* should only be used for evaluating errors. *" << endl; + report( NOTICE, kModule ) << "* The results with different seeds will be random *" << endl; + report( NOTICE, kModule ) << "* due to random oversampling of the input file. *" << endl; + report( NOTICE, kModule ) << "* Random Seed: " << std::setw(7) << randSeed << " *" << endl; + report( NOTICE, kModule ) << "****************************************************" << endl; + report( NOTICE, kModule ) << endl; TH1::AddDirectory( kFALSE ); gSystem->Load( "libTree" ); @@ -114,7 +125,7 @@ FSRootDataReaderBootstrap::FSRootDataReaderBootstrap( const vector< string >& ar fileexists.close(); m_inFile = new TFile( inFileName.c_str() ); if (!m_inFile || m_inFile->IsZombie()) { - cout << "FSRootDataReaderBootstrap WARNING: Cannot open file... " << inFileName << endl; + report( ERROR, kModule ) << "FSRootDataReaderBootstrap WARNING: Cannot open file... " << inFileName << endl; m_inFile = NULL; m_inTree = NULL; return; @@ -122,7 +133,7 @@ FSRootDataReaderBootstrap::FSRootDataReaderBootstrap( const vector< string >& ar m_inTree = static_cast( m_inFile->Get( inTreeName.c_str() ) ); if (!m_inTree) { - cout << "FSRootDataReaderBootstrap WARNING: Cannot open tree... " << inTreeName << endl; + report( ERROR, kModule ) << "FSRootDataReaderBootstrap WARNING: Cannot open tree... " << inTreeName << endl; m_inFile->Close(); delete m_inFile; m_inFile = NULL; @@ -134,19 +145,13 @@ FSRootDataReaderBootstrap::FSRootDataReaderBootstrap( const vector< string >& ar m_inTree->AddFriend(friendTreeName, friendFileName); } else{ - cout << "FSRootDataReaderBootstrap WARNING: Cannot find file... " << inFileName << endl; + report( ERROR, kModule ) << "FSRootDataReaderBootstrap WARNING: Cannot find file... " << inFileName << endl; m_inFile = NULL; m_inTree = NULL; return; } - - cout << "Opening Tree: " << inFileName << " " << inTreeName << " (numParticles=" << m_numParticles << ")"; - if (fourMomentumPrefix != "") cout << " fourMomentumPrefix=" << fourMomentumPrefix; - if (friendFileName != "") cout << " friendFile=" << friendFileName << " friendTree=" << friendTreeName; - if (weightBranchName != "weight") cout << " weightBranch=" << weightBranchName; - cout << " randSeed=" << randSeed << endl << endl; - if (m_inTree){ + if (m_inTree){ TString sEnPB = fourMomentumPrefix+"EnPB"; TString sPxPB = fourMomentumPrefix+"PxPB"; TString sPyPB = fourMomentumPrefix+"PyPB"; diff --git a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.h b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.h index edfe5dbee..3d02fbc7e 100644 --- a/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.h +++ b/src/libraries/AMPTOOLS_DATAIO/FSRootDataReaderBootstrap.h @@ -55,6 +55,7 @@ class FSRootDataReaderBootstrap : public UserDataReader< FSRootDataReaderBootstr std::multiset m_entryOrder; // Stores bootstrap sampled event indices mutable std::multiset::const_iterator m_nextEntry; // Iterator for sampling + static const char* kModule; }; #endif