diff --git a/CMakeLists.txt b/CMakeLists.txt
index fc8c37f..c17f110 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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)
diff --git a/inc/pluto_private.h b/inc/pluto_private.h
new file mode 100644
index 0000000..772bd86
--- /dev/null
+++ b/inc/pluto_private.h
@@ -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 */
diff --git a/src/PAngularDistribution.cc b/src/PAngularDistribution.cc
index 822a78c..51c953e 100644
--- a/src/PAngularDistribution.cc
+++ b/src/PAngularDistribution.cc
@@ -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
@@ -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"
 
@@ -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);
 
@@ -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;
 
@@ -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))
@@ -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;
     }
 
@@ -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
@@ -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) {
@@ -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();
@@ -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;
@@ -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;
 }
 
@@ -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;
@@ -323,7 +329,7 @@ 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
@@ -331,7 +337,7 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
 	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());
@@ -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());
@@ -350,11 +356,11 @@ 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
@@ -362,7 +368,7 @@ Bool_t PAngularDistribution::Rotate(Int_t rotate_daughters) {
 	for (int i=0; i<loc_n_daughters; i++) {
 	    daughter[i]->RotateZ(-ang_tmp.Phi());
 	    daughter[i]->RotateY(-ang_tmp.Theta());
-	}	    
+	}
     }
 
     return kTRUE;
@@ -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++) {
@@ -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;
 }
@@ -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;
     }
@@ -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;
 };
 
@@ -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;
@@ -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) {
@@ -510,7 +516,7 @@ void PAngularDistribution::Print(const Option_t *) const {
     } else if (anglesg) {
         cout << "<TGraph>";
     } else cout << "NONE";
-    
+
     cout << endl;
 }
 
diff --git a/src/PChannelModel.cc b/src/PChannelModel.cc
index 4acb04f..e4ec5b9 100644
--- a/src/PChannelModel.cc
+++ b/src/PChannelModel.cc
@@ -5,14 +5,15 @@
 //                                  Author:  I. Froehlich
 /////////////////////////////////////////////////////////////////////
 
+#include "pluto_private.h"
 
-using namespace std;
 #include <sstream>
 #include <iostream>
 #include <iomanip>
 
 #include "PChannelModel.h"
 
+using namespace std;
 
 ClassImp(PChannelModel)
 
@@ -46,11 +47,11 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 	Int_t alt_position = 0,
 	    path_position  = 0;
 	for (Int_t i=strlen(id); i>=0; i--) {
-	    if (id[i] == '@') 
+	    if (id[i] == '@')
 		alt_position = i;
 	}
 
-	if (alt_position) { 
+	if (alt_position) {
 	    alt_position++; //skip "@"
 	}
 
@@ -58,11 +59,11 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 	strcpy(db_id, &(id[alt_position]));
 
 	for (Int_t i=strlen(id); i>=alt_position; i--) {
-	    if (id[i] == path_delim) 
+	    if (id[i] == path_delim)
 		path_position = i;
 	}
 
-	if (path_position) { 
+	if (path_position) {
 	    path_position++; //skip "/"
 	    model_def_key = makeStaticData()->
 		MakeDirectoryEntry("modeldef", NMODEL_NAME, LMODEL_NAME, &(id[path_position]));
@@ -74,7 +75,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 
 	PUtils::Tokenize(&(id[alt_position]), parser_delim, arr1, &arr1_s);
 
-	if (path_position) { 
+	if (path_position) {
 	    char *nid = new char[strlen(id)+1];
 	    strcpy(nid,id);
 	    nid[path_position-1] = path_delim;
@@ -101,8 +102,8 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 
 	    Int_t arr1_id[10+2];
 	    arr1_id[0] = makeStaticData()->GetParticleID(arr1[0]);
-	    	   
-	    for (int i=1; i<=(arr1_s-2); i++) 
+
+	    for (int i=1; i<=(arr1_s-2); i++)
 		arr1_id[i] = makeStaticData()->GetParticleID(arr1[i+1]);
 
 	    key = makeStaticData()->GetDecayKey(arr1_id, arr1_s-2);
@@ -117,7 +118,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 		sec_key = makeStaticData()->GetSecondaryKey(key, model_def_key);
 		if (sec_key > 0) {
 		    key = sec_key;
-		} else {		    
+		} else {
 		    key = makeStaticData()->AddAlias(makeStaticData()->GetDecayNameByKey(key), db_id);
 		    //defkey has to be stored for future use!
 		    makeDataBase()->SetParamInt (key, "defkey", &model_def_key);
@@ -128,12 +129,12 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 
 	    //To make life more easy, add "default" track template
  	    Add(arr1[0], "parent");
- 	    for (int i=1; i<=(arr1_s-2); i++) 
+ 	    for (int i=1; i<=(arr1_s-2); i++)
  		Add(arr1[i+1], "daughter");
  	    preliminary_particles = 1;
 	} // end "decay"
 
-	if (alt_position) { 
+	if (alt_position) {
 	    //reset ID:
 	    char *newid = new char[alt_position];
 	    strncpy(newid, id, alt_position);
@@ -154,7 +155,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
     if (sec_key < 0) { //No sec. model
 	is_channel = makeStaticData()->GetDecayIdxByKey(key);
 	is_pid     = makeStaticData()->GetParticleIDByKey(key);
-    } 
+    }
 
     if (is_pid<0 && !is_channel) {
 	Warning("PChannelModel", "Undetermined particle for key %i", key);
@@ -170,7 +171,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 
     draw_option = 0;
     didx_option = -1;
-    
+
     //some TF1 defaults...
 
     //ROOT6
@@ -180,7 +181,11 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
     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
     for (int i = 0; i < fNpar; i++) {
 	fParErrors[i]  = 0;
 	fParMin[i]     = 0;
@@ -225,7 +230,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
 
     mc_max = 1000; //for the integration methods in Width-mesh calculation
 
-    if (is_pid >= 0) 
+    if (is_pid >= 0)
 	SetDynamicRange(PData::LMass(is_pid), PData::UMass(is_pid));
 
     if ((is_channel>=0) && (sec_key<0)) {
@@ -237,7 +242,7 @@ PChannelModel::PChannelModel(const Char_t *id, const Char_t *de, Int_t key) :
     didx_param     = makeDataBase()->GetParamInt("didx");
     scfactor_param = makeDataBase()->GetParamDouble("scfactor");
     unstable_width = makeStaticData()->GetBatchValue("_system_unstable_width");
-    
+
 };
 
 PDistribution *PChannelModel::Clone(const char *) const {
@@ -267,7 +272,7 @@ PChannelModel *PChannelModel::GetSecondaryModel(const char *name) {
 Bool_t PChannelModel::SampleMass(Double_t*,  Int_t*) {
     //Samples the masses of the known decay products (or single particle)
     //Order of particles is the same as in data base
-    
+
     return kFALSE;
 }
 
@@ -279,7 +284,7 @@ Double_t PChannelModel::EvalPar(const Double_t *x, const Double_t *params) {
 
     return Eval(x[0]);
 }
- 
+
 Double_t PChannelModel::Eval(Double_t x, Double_t, Double_t, Double_t) const {
     Double_t res;
     Double_t my_x[11];
@@ -309,7 +314,7 @@ Double_t PChannelModel::Eval(Double_t x, Double_t, Double_t, Double_t) const {
 }
 
 Double_t PChannelModel::GetWeight(Double_t *mass, Int_t *didx) {
-    //Get the weight for the masses of the known decay products 
+    //Get the weight for the masses of the known decay products
     //(or single particle)
     //Decay: Order of particles is the same as in data base,
     //mass[0] is the parent mass
@@ -327,7 +332,7 @@ Double_t PChannelModel::GetWeight(Double_t *mass, Int_t *didx) {
 }
 
 TComplex PChannelModel::GetAmplitude(Double_t *mass, Int_t *didx) {
-    //Get the amplitude for the masses of the known decay products 
+    //Get the amplitude for the masses of the known decay products
     //(or single particle)
     //Decay: Order of particles is the same as in data base,
     //mass[0] is the parent mass
@@ -374,7 +379,7 @@ Bool_t PChannelModel::GetBR(Double_t mass, Double_t *br, Double_t totalwidth) {
 
     if (is_pid >= 0)
 	return kFALSE;
-    
+
     Double_t lwidth;
 
     //test if the GetWidth was overloaded
@@ -383,7 +388,7 @@ Bool_t PChannelModel::GetBR(Double_t mass, Double_t *br, Double_t totalwidth) {
  	*br = makeStaticData()->GetDecayBR(is_channel);
  	return kTRUE;
     }
-    
+
     //Make everything coherent to Dyn.Data
     Double_t sc = 1.;
     Double_t *scfactor = &sc;
@@ -395,13 +400,13 @@ Bool_t PChannelModel::GetBR(Double_t mass, Double_t *br, Double_t totalwidth) {
     if (totalwidth < 0.)
 	totalwidth = makeStaticData()->GetParticleTotalWidth(makeStaticData()->GetDecayParent(is_channel));
 
-    if (totalwidth <= (*unstable_width)) 
+    if (totalwidth <= (*unstable_width))
 	*br = makeStaticData()->GetDecayBR(is_channel);
     else  {
 	//if we have a local width, the BR is the partial width divided by total
 	*br = lwidth/totalwidth;
     }
-    return kTRUE; 
+    return kTRUE;
 }
 
 int PChannelModel::GetDepth(int) {