Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,14 @@ ROOT_GENERATE_DICTIONARY(G__Pluto ${INCL_FILES} LINKDEF PlutoLinkDef.h)
#---Create a shared library with geneated dictionary
add_library(Pluto SHARED ${SOURCES} G__Pluto.cxx ${CMAKE_CURRENT_BINARY_DIR}/Compiled.h)
target_link_libraries(Pluto ${ROOT_LIBRARIES})
target_include_directories(Pluto PRIVATE
inc
)

SET(ROOT_VERSION_LEGACY 6.22.00)
if (${ROOT_VERSION} VERSION_LESS ${ROOT_VERSION_LEGACY})
target_compile_definitions(LEGACY_ROOT)
endif()

#---Create a main program using the library
#add_executable(Main MainEvent.cxx)
Expand Down
17 changes: 17 additions & 0 deletions inc/pluto_private.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
#ifndef PLUTO_PRIVATE_H
#define PLUTO_PRIVATE_H

#include <memory>

#if __cplusplus < 201402L
namespace std {
template <typename T, typename... Args> std::unique_ptr<T> make_unique(Args&&... args)
{
return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
}
}
#else
// using std::make_unique;
#endif

#endif /* PLUTO_PRIVATE_H */
98 changes: 52 additions & 46 deletions src/PAngularDistribution.cc
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
//
//
// If "align" is set, the reference will be the composite particle
// of the primary + align
// of the primary + align
// This can be used for the alignment of two particles
//
// The 2nd parameter for TF2 is the mass of the reference
Expand All @@ -30,15 +30,18 @@
// Author: I. Froehlich
// Written: 01.10.2006
// Revised: 24.10.2010
//
//
//
//
/////////////////////////////////////////////////////////////////////

#include "pluto_private.h"

using namespace std;
#include <sstream>
#include <iostream>
#include <iomanip>
#include <memory>

using namespace std;

#include "PAngularDistribution.h"

Expand Down Expand Up @@ -77,12 +80,15 @@ PAngularDistribution::PAngularDistribution(const Char_t *id,const Char_t *de) :
fParErrors.resize(fNpar);
fParMin.resize(fNpar);
fParMax.resize(fNpar);
#ifdef LEGACY_ROOT
fParams = new TF1Parameters(fNpar);

#else
fParams = std::move(std::make_unique<TF1Parameters>(fNpar));
#endif
fParErrors[0] = 0;
fParMin[0] = 0;
fParMax[0] = 0;

fParams->SetParName(0, "Mass reference for 2-dim angular distributions");
fParams->SetParameter(0, 0);

Expand All @@ -105,7 +111,7 @@ PAngularDistribution::PAngularDistribution(const Char_t *id,const Char_t *de) :
fParams[0]=0;
#endif

fNpx = 1000;
fNpx = 1000;
fXmin = -1;
fXmax = 1;

Expand Down Expand Up @@ -141,7 +147,7 @@ Bool_t PAngularDistribution::Init(void) {
mass_reference = GetParticle("mass_reference");
ang_reference = GetParticle("ang_reference");
align = GetParticle("align");

if (align) {//check if align is a daughter or parent
if ((current_flag == PARTICLE_LIST_PARENT) ||
(current_flag == PARTICLE_LIST_DAUGHTER))
Expand All @@ -166,11 +172,11 @@ Bool_t PAngularDistribution::Init(void) {
}

if (!reference) reference = parent;

for (int i=0; i<position; i++) {
if ((particle_flag[i] & PARTICLE_LIST_GRANDPARENT)
| (particle_flag[i] & PARTICLE_LIST_GRANDGRANDPARENT)
| (particle_flag[i] & PARTICLE_LIST_SIBLING))
if ((particle_flag[i] & PARTICLE_LIST_GRANDPARENT)
| (particle_flag[i] & PARTICLE_LIST_GRANDGRANDPARENT)
| (particle_flag[i] & PARTICLE_LIST_SIBLING))
check_abort = kTRUE;
}

Expand All @@ -180,11 +186,11 @@ Bool_t PAngularDistribution::Init(void) {
daughter[n_daughters] = GetParticle("daughter");
if (!daughter[n_daughters]) {
myloop = 0;
}
}
}
n_daughters--;

direct_sampling_possible = kFALSE;
direct_sampling_possible = kFALSE;

if (!align && // no composite
!mass_reference && // too complicated
Expand All @@ -195,7 +201,7 @@ Bool_t PAngularDistribution::Init(void) {
direct_sampling_possible=kTRUE;
//N.B. this is overwritten by inherited models (e.g.PDeltaAngularDistribution)
//cout << direct_sampling_possible << endl;
return kTRUE;
return kTRUE;
};

Bool_t PAngularDistribution::Prepare(void) {
Expand All @@ -220,25 +226,25 @@ Bool_t PAngularDistribution::SampleAngle(void) {
//the rejection method. This works
//in the simple cases, and only if there
//are no correlated angular distributions

if (!direct_sampling_possible) return kFALSE;

// primary->Print();
// for (int i=0;i<n_daughters;i++) daughter[i]->Print();
if (!Rotate(1)) return kFALSE;

if (!Rotate(1)) return kFALSE;

// cout << "---" << endl;
// primary_tmp.Print();
// for (int i=0;i<n_daughters;i++) daughter[i]->Print();
// cout << "*********************" << endl;

//now we rotate all daughters such that primary_tmp is pointing to the z-axis
double prim_phi = primary_tmp.Phi();
double prim_theta = primary_tmp.Theta();

#if 1

//primary->Print();
// primary_tmp.Print();
// for (int i=0;i<n_daughters;i++) daughter[i]->Print();
Expand All @@ -251,7 +257,7 @@ Bool_t PAngularDistribution::SampleAngle(void) {
}

//sample the polar angle
double polar_angle = acos(SamplePolarAngle(0.5*(cos(prim_theta)+1)));
double polar_angle = acos(SamplePolarAngle(0.5*(cos(prim_theta)+1)));
//2r=cos_theta-1 was the definition in the baryon_cos alg of the original PChannel
//polar_angle=prim_theta;
//cout << polar_angle << endl;
Expand All @@ -272,13 +278,13 @@ Bool_t PAngularDistribution::SampleAngle(void) {
daughter[i]->RotateZ(prim_phi);
#endif
//rotate back to the old frame
if (!RotateBack(1)) return kFALSE;
if (!RotateBack(1)) return kFALSE;

//if all worked out we can copy the primary :
*primary = primary_tmp;

direct_sampling_done = kTRUE;

return kTRUE;
}

Expand All @@ -289,10 +295,10 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
//are rotated as well. This is needed for direct samplings

Int_t loc_n_daughters = n_daughters * rotate_daughters;

primary_tmp = primary; //particle under investigation. Make better a copy
primary_tmp.Boost(parent->BoostVector()); // go back to lab frame
for (int i=0; i<loc_n_daughters; i++)
for (int i=0; i<loc_n_daughters; i++)
daughter[i]->Boost(parent->BoostVector());
PParticle compound(primary_tmp);
PParticle atmp;
Expand Down Expand Up @@ -323,15 +329,15 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
} else {
primary_tmp.Boost(-reference->BoostVector());
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->Boost(-reference->BoostVector());
daughter[i]->Boost(-reference->BoostVector());
}
}
} else { // first go to 2nd reference (e.g. base_reference=eta), then
// rotate and boost to 1st reference (e.g. reference=dilepton)
reference_tmp = reference;

primary_tmp.Boost(-base_reference->BoostVector()); // daughter in base_reference
ang_tmp.Boost(-base_reference->BoostVector());
ang_tmp.Boost(-base_reference->BoostVector());
reference_tmp.Boost(-base_reference->BoostVector()); // reference in base_reference
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->Boost(-base_reference->BoostVector());
Expand All @@ -341,7 +347,7 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
primary_tmp.RotateY(-reference_tmp.Theta());
primary_tmp.Boost(0, 0, -reference_tmp.Beta());
ang_tmp.RotateZ(-reference_tmp.Phi());
ang_tmp.RotateY(-reference_tmp.Theta());
ang_tmp.RotateY(-reference_tmp.Theta());
ang_tmp.Boost(0, 0, -reference_tmp.Beta());
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->RotateZ(-reference_tmp.Phi());
Expand All @@ -350,19 +356,19 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
}
} else {
primary_tmp.Boost(-reference_tmp.BoostVector()); // go to reference
for (int i=0; i<loc_n_daughters; i++)
for (int i=0; i<loc_n_daughters; i++)
daughter[i]->Boost(-reference_tmp.BoostVector());
}
}

//finally take ang_reference into account
if (ang_reference && rotate) {
primary_tmp.RotateZ(-ang_tmp.Phi()); // rotate daughter reference onto z-axis
primary_tmp.RotateY(-ang_tmp.Theta());
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->RotateZ(-ang_tmp.Phi());
daughter[i]->RotateY(-ang_tmp.Theta());
}
}
}

return kTRUE;
Expand All @@ -379,12 +385,12 @@ Bool_t PAngularDistribution::RotateBack(Int_t rotate_daughters) {
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->RotateY(ang_tmp.Theta());
daughter[i]->RotateZ(ang_tmp.Phi());
}
}
}

if (base_reference==NULL) { // no 2nd reference frame
if (rotate) {
primary_tmp.Boost(0, 0, reference->Beta()); // go to 1st reference frame
primary_tmp.Boost(0, 0, reference->Beta()); // go to 1st reference frame
primary_tmp.RotateY(reference->Theta());
primary_tmp.RotateZ(reference->Phi());
for (int i=0; i<loc_n_daughters; i++) {
Expand All @@ -408,16 +414,16 @@ Bool_t PAngularDistribution::RotateBack(Int_t rotate_daughters) {
daughter[i]->Boost(0, 0, reference_tmp.Beta());
daughter[i]->RotateY(reference_tmp.Theta());
daughter[i]->RotateZ(reference_tmp.Phi());
}
}

primary_tmp.Boost(base_reference->BoostVector()); // daughter in base_reference
reference_tmp.Boost(base_reference->BoostVector()); // reference in base_reference
for (int i=0; i<loc_n_daughters; i++) {
daughter[i]->Boost(base_reference->BoostVector());
daughter[i]->Boost(base_reference->BoostVector());
}
}
primary_tmp.Boost(-parent->BoostVector()); // go back to parent frame
for (int i=0; i<loc_n_daughters; i++)
for (int i=0; i<loc_n_daughters; i++)
daughter[i]->Boost(-parent->BoostVector());
return kTRUE;
}
Expand All @@ -427,11 +433,11 @@ Bool_t PAngularDistribution::IsNotRejected(void) {
if (direct_sampling_done) return kTRUE;

Double_t tmp_c0, f;

if (!Rotate(0)) return kFALSE;

tmp_c0 = cos(primary_tmp.Theta());

if ((reflection_symmetry) && (tmp_c0<0.)) {
tmp_c0 = -tmp_c0;
}
Expand All @@ -453,9 +459,9 @@ Bool_t PAngularDistribution::IsNotRejected(void) {
Warning("IsNotRejected", "[%s] Weight > max, new max is %lf", GetName(), weight_max);
}

if ((f/weight_max) > PUtils::sampleFlat())
if ((f/weight_max) > PUtils::sampleFlat())
return kTRUE; // sample now distribution

return kFALSE;
};

Expand Down Expand Up @@ -486,7 +492,7 @@ Double_t PAngularDistribution::EvalPar(const Double_t *x, const Double_t *params
}
return Eval(x[0]);
}

Double_t PAngularDistribution::Eval(Double_t x, Double_t, Double_t, Double_t) const
{
if ((reflection_symmetry) && (x<0.)) x = -x;
Expand All @@ -501,7 +507,7 @@ void PAngularDistribution::Print(const Option_t *) const {
BasePrint();

cout << " Formula used: ";
if (angles1) {
if (angles1) {
if (angles1->GetExpFormula() != TString("")) cout << angles1->GetExpFormula();
else cout << "<compiled>";
} else if (angles2) {
Expand All @@ -510,7 +516,7 @@ void PAngularDistribution::Print(const Option_t *) const {
} else if (anglesg) {
cout << "<TGraph>";
} else cout << "NONE";

cout << endl;
}

Expand Down
Loading