Skip to content
Snippets Groups Projects
Commit d9aabc09 authored by Bertram Kopf's avatar Bertram Kopf
Browse files

test application added: PWA with covariant tensors and Minuit2

parent 4e3cf138
No related branches found
No related tags found
No related merge requests found
...@@ -63,8 +63,10 @@ INSTALL ( TARGETS ${EXECUTABLENAME1} DESTINATION bin/${EXECUTABLENAME}) ...@@ -63,8 +63,10 @@ INSTALL ( TARGETS ${EXECUTABLENAME1} DESTINATION bin/${EXECUTABLENAME})
SET (EXECUTABLENAME EtacToa1320pi0fitApp) SET (EXECUTABLENAME EtacToa1320pi0fitApp)
SET (${EXECUTABLENAME}_SRCS SET (${EXECUTABLENAME}_SRCS
EtacToa1320pi0Fcn.cc
EtacToa1320pi0fit.cc EtacToa1320pi0fit.cc
EtacToa1320pi0fitMain.cc
) )
ADD_EXECUTABLE(${EXECUTABLENAME} ADD_EXECUTABLE(${EXECUTABLENAME}
......
#ifndef _etacToa1320pi0Data_H
#define _etacToa1320pi0Data_H
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
struct fitParamVal
{
double a1320Mass;
double a1320Width;
double cont0spin;
double cont1spin;
double cont2spin;
};
struct evt4Vec
{
Vector4<double> pi0_0_4Vec;
Vector4<double> pi0_1_4Vec;
Vector4<double> eta_4Vec;
Vector4<double> a2_0_4Vec;
Vector4<double> a2_1_4Vec;
Vector4<double> cm_4Vec;
Tensor<complex<double> > spin2_0_amp;
Tensor<complex<double> > spin2_1_amp;
Tensor<complex<double> > spin1_0_amp;
Tensor<complex<double> > spin1_1_amp;
};
#endif /* _etacToa1320pi0Data_H */
//#include <getopt.h>
//#include <fstream>
//#include <string>
#include "Examples/qft++/EtacToa1320pi0Fcn.hh"
#include "Examples/qft++/EtacToa1320pi0Data.hh"
#include "Examples/qft++/EtacToa1320pi0fit.hh"
#include <cassert>
using namespace ROOT::Minuit2;
EtacToa1320pi0Fcn::EtacToa1320pi0Fcn(EtacToa1320pi0fit* etacToa1320pi0fit) :
_etacToa1320pi0fit(etacToa1320pi0fit)
{
if (0==_etacToa1320pi0fit){
std::cout << "EtacToa1320pi0fit pointer is 0 !!!!";
assert(0);
}
}
EtacToa1320pi0Fcn::~EtacToa1320pi0Fcn()
{
}
double EtacToa1320pi0Fcn::operator()(const std::vector<double>& par) const
{
fitParamVal theFitParmValTmp;
assert(_etacToa1320pi0fit->setFitParamVal(theFitParmValTmp, par));
double result=_etacToa1320pi0fit->calcLogLh(theFitParmValTmp);
std::cout << "InterMassFit= " << theFitParmValTmp.a1320Mass
<< " InterWidthFit= " << theFitParmValTmp.a1320Width
<< " spin0= " << theFitParmValTmp.cont0spin
<< " spin1= " << theFitParmValTmp.cont1spin
<< " spin2= " << theFitParmValTmp.cont2spin
<< " logLH= " << result << std::endl;
return result;
}
double EtacToa1320pi0Fcn::Up() const
{
return .5;
}
#ifndef _etacToa1320pi0fcn_H
#define _etacToa1320pi0fcn_H
#include <iostream>
#include <fstream>
//#include <string>
#include <vector>
//#include <cassert>
#include "Minuit2/FCNBase.h"
class EtacToa1320pi0fit;
namespace ROOT {
namespace Minuit2 {
class EtacToa1320pi0Fcn : public FCNBase {
public:
EtacToa1320pi0Fcn(EtacToa1320pi0fit* etacToa1320pi0fit);
virtual ~EtacToa1320pi0Fcn();
double operator()(const std::vector<double>& par) const;
double Up() const;
private:
EtacToa1320pi0fit* _etacToa1320pi0fit;
};
} // namespace Minuit2
} // namespace ROOT
#endif /* _etacToa1320pi0fcn_H */
#include <getopt.h>
#include <fstream>
#include <string>
#include "Examples/qft++/EtacToa1320pi0fit.hh"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
EtacToa1320pi0fit::EtacToa1320pi0fit(int kindOfData) :
_kindOfData(kindOfData),
_nOfData(10000),
_nOfMc(30000),
_fitHistsFilled(false),
_theTFile(0),
_dalitzDataHist(0),
_dalitzMcHist(0),
_dalitzFittedHist(0),
_invpietaDataHist(0),
_invpietaMcHist(0),
_invpietaFittedHist(0),
_invpipiDataHist(0),
_invpipiMcHist(0),
_invpipiFittedHist(0)
{
std::string theSourcePath=getenv("CMAKE_SOURCE_DIR");
_mcPath=theSourcePath+std::string("/Examples/qft++/data/mc100Mgev.dat");
if (kindOfData==0) _dataPath=theSourcePath+std::string("/Examples/qft++/data/dataEtacToA0Pi_100Mgev.dat");
else if (kindOfData==2) _dataPath=theSourcePath+std::string("/Examples/qft++/data/dataSpin2100MgevNew.dat");
else{
std::cout <<"this kind of data: " << kindOfData << " is not supported!!!!\n"
<<"initialze EtacToa1320pi0fit either with 0 or with 2 !!!!"
<< std::endl;
assert(0);
}
initRootStuff();
read4Vecs(_dataPath, _nOfData, _data4Vecs);
read4Vecs(_mcPath,_nOfMc, _mc4Vecs);
for (size_t it=0; it<_data4Vecs.size(); it++)
{
plotDalitz(_dalitzDataHist, _data4Vecs[it], 1);
plotInvPiEta(_invpietaDataHist, _data4Vecs[it], 1);
plotInvPiPi(_invpipiDataHist, _data4Vecs[it], 1);
}
for (size_t it=0; it<_mc4Vecs.size(); it++)
{
plotDalitz(_dalitzMcHist, _mc4Vecs[it], 1);
plotInvPiEta(_invpietaMcHist, _mc4Vecs[it], 1);
plotInvPiPi(_invpipiMcHist, _mc4Vecs[it], 1);
}
}
EtacToa1320pi0fit::~EtacToa1320pi0fit()
{
_theTFile->Write();
_theTFile->Close();
// if (0!=_theTFile) delete _theTFile;
// if (0!=_dalitzDataHist) delete _dalitzDataHist;
// if (0!=_dalitzMcHist) delete _dalitzMcHist;
// if (0!=_dalitzFittedHist) delete _dalitzFittedHist;
// if (0!=_invpietaDataHist) delete _invpietaDataHist;
// if (0!=_invpietaMcHist) delete _invpietaMcHist;
// if (0!=_invpietaFittedHist) delete _invpietaFittedHist;
// if (0!=_invpipiDataHist) delete _invpipiDataHist;
// if (0!=_invpipiMcHist) delete _invpipiMcHist;
// if (0!=_invpipiFittedHist) delete _invpipiFittedHist;
}
void EtacToa1320pi0fit::initRootStuff()
{
_theTFile=new TFile("./a1320Fit.root","recreate");
_dalitzDataHist= new TH2F("_dalitzDataHist","dalitzplot data",50, 0., 8.,50, 0., 8.);
_dalitzMcHist= new TH2F("_dalitzMcHist","dalitzplot mc",50, 0., 8.,50, 0., 8.);
_dalitzFittedHist= new TH2F("_dalitzFittedHist","dalitzplot fit",50, 0., 8.,50, 0., 8.);
_invpietaDataHist= new TH1F("_invpietaDataHist","M(#pi #eta) data",120, 0.6, 3.0);
_invpietaMcHist= new TH1F("_invpietaMcHist","M(#pi #eta) mc",120, 0.6, 3.0);
_invpietaFittedHist= new TH1F("_invpietaFittedHist","M(#pi #eta) fit",120, 0.6, 3.0);
_invpipiDataHist= new TH1F("_invpipiDataHist","M(#pi #pi) data",120, 0.2, 2.6);
_invpipiMcHist= new TH1F("_invpipiMcHist","M(#pi #pi) mc",120, 0.2, 2.6);
_invpipiFittedHist= new TH1F("_invpipiFittedHist","M(#pi #pi) fit",120, 0.2, 2.6);
}
void EtacToa1320pi0fit::read4Vecs(std::string& path, int nEvts, std::vector<evt4Vec>& the4Vecs)
{
std::cout << "calculate 4Vecs and amplitudes for " << path << std::endl;
std::ifstream inputStream(path.c_str(), std::ios::in);
if (!inputStream.good())
{
std::cout << "Input file doesn't exsits !!! ";
assert(0);
}
int counter=0;
while (!inputStream.eof() && counter<nEvts )
{
if ( counter%1000 == 0 ) std::cout << "event " << counter << std::endl;
Vector4<double> pi14V,pi24V,eta4V; // 4-momenta
get4Vecs(inputStream, pi14V);
get4Vecs(inputStream, pi24V);
get4Vecs(inputStream, eta4V);
Vector4<double> cm_4V(pi14V+pi24V+eta4V);
Vector4<double> pi1_cm_4V(pi14V);
pi1_cm_4V.Boost(cm_4V);
Vector4<double> pi2_cm_4V(pi24V);
pi2_cm_4V.Boost(cm_4V);
Vector4<double> eta_cm_4V(eta4V);
eta_cm_4V.Boost(cm_4V);
Vector4<double> a2_0_cm_4Vec(pi1_cm_4V+eta_cm_4V);
Vector4<double> a2_1_cm_4Vec(pi2_cm_4V+eta_cm_4V);
Vector4<double> cm_cm_4V(cm_4V);
cm_cm_4V.Boost(cm_4V);
evt4Vec theEvt4Vec;
theEvt4Vec.pi0_0_4Vec=pi1_cm_4V;
theEvt4Vec.pi0_1_4Vec=pi2_cm_4V;
theEvt4Vec.eta_4Vec=eta_cm_4V;
theEvt4Vec.a2_0_4Vec=a2_0_cm_4Vec;
theEvt4Vec.a2_1_4Vec=a2_1_cm_4Vec;
theEvt4Vec.cm_4Vec=cm_cm_4V;
theEvt4Vec.spin2_0_amp = calcSpin2Amp(theEvt4Vec.pi0_0_4Vec, theEvt4Vec.eta_4Vec, theEvt4Vec.pi0_1_4Vec);
theEvt4Vec.spin2_1_amp = calcSpin2Amp(theEvt4Vec.pi0_1_4Vec, theEvt4Vec.eta_4Vec, theEvt4Vec.pi0_0_4Vec);
theEvt4Vec.spin1_0_amp = calcSpin1Amp(theEvt4Vec.pi0_0_4Vec, theEvt4Vec.eta_4Vec, theEvt4Vec.pi0_1_4Vec);
theEvt4Vec.spin1_1_amp = calcSpin1Amp(theEvt4Vec.pi0_1_4Vec, theEvt4Vec.eta_4Vec, theEvt4Vec.pi0_0_4Vec);
the4Vecs.push_back(theEvt4Vec);
counter++;
}
}
void EtacToa1320pi0fit::get4Vecs (std::ifstream& inStream, Vector4<double>& the4Vec)
{
double tmpPx,tmpPy, tmpPz, tmpE;
inStream >> tmpPx;
inStream >> tmpPy;
inStream >> tmpPz;
inStream >> tmpE;
the4Vec.SetP4(tmpE, tmpPx, tmpPy, tmpPz);
}
void EtacToa1320pi0fit::plotDalitz(TH2F* theHisto, evt4Vec& theEvt4Vecs, double weight)
{
Vector4<double> pi0eta4V=theEvt4Vecs.pi0_0_4Vec + theEvt4Vecs.eta_4Vec;
Vector4<double> pi1eta4V=theEvt4Vecs.pi0_1_4Vec + theEvt4Vecs.eta_4Vec;
theHisto->Fill(pi0eta4V.M()*pi0eta4V.M(),pi1eta4V.M()*pi1eta4V.M(), weight);
theHisto->Fill(pi1eta4V.M()*pi1eta4V.M(),pi0eta4V.M()*pi0eta4V.M(), weight);
}
void EtacToa1320pi0fit::plotInvPiEta (TH1F* theHisto, evt4Vec& theEvt4Vecs, double weight)
{
Vector4<double> pi0eta4V=theEvt4Vecs.pi0_0_4Vec + theEvt4Vecs.eta_4Vec;
Vector4<double> pi1eta4V=theEvt4Vecs.pi0_1_4Vec + theEvt4Vecs.eta_4Vec;
theHisto->Fill(pi0eta4V.M(), weight);
theHisto->Fill(pi1eta4V.M(), weight);
}
void EtacToa1320pi0fit::plotInvPiPi (TH1F* theHisto, evt4Vec& theEvt4Vecs, double weight)
{
Vector4<double> pi0pi0=theEvt4Vecs.pi0_0_4Vec + theEvt4Vecs.pi0_1_4Vec;
theHisto->Fill(pi0pi0.M(), weight);
}
Tensor<complex<double> > EtacToa1320pi0fit::calcSpin2Amp(Vector4<double>& pi_cm, Vector4<double>& eta_cm, Vector4<double>& pi_recoil_cm)
{
Tensor<complex<double> > result(0); //rank-0 Tensor for the amplitude
Vector4<double> a1320_cm=eta_cm + pi_cm;
PolVector a1320_Pol_cm(2);
a1320_Pol_cm.SetP4(a1320_cm,a1320_cm.M());
OrbitalTensor orb_eta_To_a2_piRecoil(2);
orb_eta_To_a2_piRecoil.SetP4(a1320_cm, pi_recoil_cm);
OrbitalTensor orb_a2_To_pieta(2);
orb_a2_To_pieta.SetP4(pi_cm, eta_cm);
result = (orb_a2_To_pieta) | ( a1320_Pol_cm.Projector() | orb_eta_To_a2_piRecoil );
return result;
}
Tensor<complex<double> > EtacToa1320pi0fit::calcSpin1Amp(Vector4<double>& pi_cm, Vector4<double>& eta_cm, Vector4<double>& pi_recoil_cm)
{
Tensor<complex<double> > result(0); // rank-0 Tensor for the amplitude
Vector4<double> a1320_cm=eta_cm + pi_cm;
PolVector a1320_Pol_cm(1);
a1320_Pol_cm.SetP4(a1320_cm,a1320_cm.M());
OrbitalTensor orb_eta_To_a2_piRecoil(1);
orb_eta_To_a2_piRecoil.SetP4(a1320_cm, pi_recoil_cm);
OrbitalTensor orb_a2_To_pieta(1);
orb_a2_To_pieta.SetP4(pi_cm, eta_cm);
result = (orb_a2_To_pieta) | ( a1320_Pol_cm.Projector() | orb_eta_To_a2_piRecoil );
return result;
}
double EtacToa1320pi0fit::calcIntensityCache(evt4Vec& theEvtVec, const fitParamVal& theParamVal)
{
double spin2cont=theParamVal.cont2spin;
double spin1cont=theParamVal.cont1spin;
double spin0cont=theParamVal.cont0spin;
double intensity=0;
Tensor<complex<double> > amp(0); // rank-0 Tensor for the amplitude
amp = theEvtVec.spin2_0_amp*BreitWigner(theEvtVec.a2_0_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin2cont*norm(amp(0,0));
amp = theEvtVec.spin2_1_amp*BreitWigner(theEvtVec.a2_1_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin2cont*norm(amp(0,0));
amp = theEvtVec.spin1_0_amp*BreitWigner(theEvtVec.a2_0_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin1cont*norm(amp(0,0));
amp = theEvtVec.spin1_1_amp*BreitWigner(theEvtVec.a2_1_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin1cont*norm(amp(0,0));
//dont know how to normalize
amp = 0.1*BreitWigner(theEvtVec.a2_0_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin0cont*norm(amp(0,0));
amp = 0.01*BreitWigner(theEvtVec.a2_1_4Vec, theParamVal.a1320Mass, theParamVal.a1320Width);
intensity += spin0cont*norm(amp(0,0));
return intensity;
}
double EtacToa1320pi0fit::calcLogLh(const fitParamVal& theParamVal)
{
double logLH=0.;
double logLH_data=0.;
for (size_t it=0; it<_data4Vecs.size(); it++)
{
double intensity=calcIntensityCache(_data4Vecs[it], theParamVal);
if (intensity>0.) logLH_data+=log10(intensity);
}
std::cout << "logLH_data= " << logLH_data << std::endl;
double LH_mc=0.;
for (size_t it=0; it<_mc4Vecs.size(); it++)
{
double intensity=calcIntensityCache(_mc4Vecs[it], theParamVal);
LH_mc+=intensity;
}
std::cout << "LH_mc= " << LH_mc << std::endl;
double logLH_mc_Norm=0.;
if (LH_mc>0.) logLH_mc_Norm=log10(LH_mc/_mc4Vecs.size());
logLH=_data4Vecs.size()/2.*(LH_mc/_mc4Vecs.size()-1)*(LH_mc/_mc4Vecs.size()-1)
-logLH_data
+_data4Vecs.size()*logLH_mc_Norm;
return logLH;
}
bool EtacToa1320pi0fit::initFitParameters(MnUserParameters& minuitParams)
{
bool result=true;
if (_kindOfData==2)
{
minuitParams.Add("InterMass", 1.6, .1, 2.3, 0.6);
minuitParams.Add("InterWidth", 0.04, .01, 0.8, 0.01);
minuitParams.Add("spin0", 0.3, .1, 1., 0.);
minuitParams.Add("spin1", 0.3, .1, 1., 0.);
minuitParams.Add("spin2", 0.3, .1, 1., 0.);
}
else //_kindOfData==0
{
minuitParams.Add("InterMass", 1.1, .1, 2.3, 0.6);
minuitParams.Add("InterWidth", 0.04, .01, 0.8, 0.01);
minuitParams.Add("spin0", 0.3, .1, 1., 0.);
minuitParams.Add("spin1", 0.3, .1, 1., 0.);
minuitParams.Add("spin2", 0.3, .1, 1., 0.);
}
return result;
}
bool EtacToa1320pi0fit::setFitParamVal(fitParamVal& fitParamVal, const std::vector<double>& par)
{
bool result=true;
fitParamVal.a1320Mass=par[0];
fitParamVal.a1320Width=par[1];
fitParamVal.cont0spin=par[2];
fitParamVal.cont1spin=par[3];
fitParamVal.cont2spin=par[4];
return result;
}
bool EtacToa1320pi0fit::fillFitHists(const fitParamVal& fitParamVal)
{
if(_fitHistsFilled)
{
std::cout <<"Not possible to fill the histos with fit results!!!!\n"
<<"They are already filled!!!!"
<<std::endl;
return false;
}
for (size_t it=0; it<_mc4Vecs.size(); it++)
{
double evtWeight=calcIntensityCache(_mc4Vecs[it], fitParamVal);
plotDalitz(_dalitzFittedHist, _mc4Vecs[it], evtWeight);
plotInvPiEta(_invpietaFittedHist, _mc4Vecs[it], evtWeight);
plotInvPiPi(_invpipiFittedHist, _mc4Vecs[it], evtWeight);
}
//normalize fitted histos
double integralData=_invpipiDataHist->Integral();
std::cout << "integralData= " << integralData << std::endl;
double integralFittedMc=_invpipiFittedHist->Integral();
std::cout << "integralFittedMc= " << integralFittedMc << std::endl;
_invpipiFittedHist->Scale(integralData/integralFittedMc);
_invpietaFittedHist->Scale(integralData/integralFittedMc);
_dalitzFittedHist->Scale(integralData/integralFittedMc);
return true;
}
#ifndef _etacToa1320pi0fit_H
#define _etacToa1320pi0fit_H
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cassert>
#include "TROOT.h"
// #include <TSystem.h>
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "Examples/qft++/EtacToa1320pi0Data.hh"
#include "Minuit2/MnUserParameters.h"
//#include "Minuit2/FCNBase.h"
// using namespace std;
using namespace ROOT::Minuit2;
class TFile;
class TH2F;
class TH1F;
class EtacToa1320pi0fit {
public:
// create/copy/destroy:
///Constructor
EtacToa1320pi0fit(int);
/** Destructor */
virtual ~EtacToa1320pi0fit();
// Getters:
double calcLogLh(const fitParamVal& theParamVal);
bool initFitParameters(MnUserParameters& minuitParams);
bool setFitParamVal(fitParamVal& fitParamVal, const std::vector<double>& par);
bool fillFitHists(const fitParamVal& fitParamVal);
protected:
private:
std::string _dataPath;
std::string _mcPath;
int _kindOfData;
int _nOfData;
int _nOfMc;
bool _fitHistsFilled;
TFile* _theTFile;
TH2F* _dalitzDataHist;
TH2F* _dalitzMcHist;
TH2F* _dalitzFittedHist;
TH1F* _invpietaDataHist;
TH1F* _invpietaMcHist;
TH1F* _invpietaFittedHist;
TH1F* _invpipiDataHist;
TH1F* _invpipiMcHist;
TH1F* _invpipiFittedHist;
std::vector<evt4Vec> _data4Vecs;
std::vector<evt4Vec> _mc4Vecs;
void initRootStuff();
void read4Vecs(std::string& path, int nEvts, std::vector<evt4Vec>& the4Vecs);
void get4Vecs (std::ifstream& inStream, Vector4<double>& the4Vec);
void plotDalitz (TH2F* theHisto, evt4Vec& theEvt4Vecs, double weight);
void plotInvPiEta (TH1F* theHisto, evt4Vec& theEvt4Vecs, double weight);
void plotInvPiPi (TH1F* theHisto, evt4Vec& theEvt4Vecs, double weight);
double calcIntensityCache(evt4Vec& theEvtVec, const fitParamVal& theParamVal);
Tensor<complex<double> > calcSpin2Amp(Vector4<double>& pi_cm, Vector4<double>& eta_cm, Vector4<double>& pi_recoil_cm);
Tensor<complex<double> > calcSpin1Amp(Vector4<double>& pi_cm, Vector4<double>& eta_cm, Vector4<double>& pi_recoil_cm);
};
#endif /* _EtacToa1320pi0fit_H */
#include <iostream>
#include <cstring>
#include <string>
#include <sstream>
#include <vector>
#include "Examples/qft++/EtacToa1320pi0Data.hh"
#include "Examples/qft++/EtacToa1320pi0fit.hh"
#include "Examples/qft++/EtacToa1320pi0Fcn.hh"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnStrategy.h"
//#include "Minuit2/MnUserTransformation.h"
using namespace ROOT::Minuit2;
int main(int __argc,char *__argv[]){
if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
|| strcmp( __argv[1], "--help" ) == 0 ) ){
std::cout << "This application is a simple PWA fit for the decay chain\n"
<< "eta_c -> intermediate + pi0 ; intermediate -> pi0 + eta\n"
<< "It makes use of the Covariant Tensor Formalism\n"
<< "The fit determines whether the intermediate resonance is a spin 0,1 or 2 particle\n"
<< "In addition the mass and width of the intermediate resonance will be fitted with a simple Breit-Wigner\n"
<< "To start the application with data containing an intermediate resonance with Spin=0 and mass 0.98 GeV, type: ./EtacToa1320pi0fitApp -d 0\n"
<< "To start the application with data containing an intermediate resonance with Spin=2 and mass 1.32 GeV, type: ./EtacToa1320pi0fitApp -d 2\n"
<< std::endl;
return 0;
}
int optind=1;
std::string dataSpinStr="";
// decode arguments
while ((optind < __argc) && (__argv[optind][0]=='-')) {
std::string sw = __argv[optind];
if (sw=="-d") {
optind++;
dataSpinStr = __argv[optind];
}
else{
cout << "Unknown switch: "
<< __argv[optind] << endl;
optind++;
}
}
std::stringstream dataSpinStrStr(dataSpinStr);
int dataSpin=2;
dataSpinStrStr >> dataSpin ;
std::cout << "dataSpin: " << dataSpin << std::endl;
EtacToa1320pi0fit* etacToa1320pi0fit=new EtacToa1320pi0fit(dataSpin);
EtacToa1320pi0Fcn fcn(etacToa1320pi0fit);
MnUserParameters upar;
if (! etacToa1320pi0fit->initFitParameters(upar)) assert(0);
MnMigrad migrad(fcn, upar);
std::cout<<"start migrad "<<std::endl;
FunctionMinimum min = migrad();
if(!min.IsValid()) {
//try with higher strategy
std::cout<<"FM is invalid, try with strategy = 2."<<std::endl;
MnMigrad migrad2(fcn, min.UserState(), MnStrategy(2));
min = migrad2();
}
// std::cout << "minimum: " << min << std::endl;
std::cout << "migrad.Fval(): " << min.Fval() << std::endl;
// std::cout<<"start Minos"<<std::endl;
// MnMinos Minos(fcn, min);
// std::pair<double,double> e0 = Minos(0);
// std::pair<double,double> e1 = Minos(1);
// std::pair<double,double> e2 = Minos(2);
// std::pair<double,double> e3 = Minos(3);
// std::pair<double,double> e4 = Minos(4);
// std::cout<<"a1320mass: "<<min.UserState().Value("a1320mass")<<" "<<e0.first<<" "<<e0.second<<std::endl;
// std::cout<<"a1320width: "<<min.UserState().Value("a1320width")<<" "<<e1.first<<" "<<e1.second<<std::endl;
// std::cout<<"spin0 content: "<<min.UserState().Value("spin0")<<" "<<e2.first<<" "<<e2.second<<std::endl;
// std::cout<<"spin1 content: "<<min.UserState().Value("spin1")<<" "<<e3.first<<" "<<e3.second<<std::endl;
// std::cout<<"spin2 content: "<<min.UserState().Value("spin2")<<" "<<e4.first<<" "<<e4.second<<std::endl;
fitParamVal theFitResult;
theFitResult.a1320Mass=min.UserState().Value("InterMass");
theFitResult.a1320Width=min.UserState().Value("InterWidth");
theFitResult.cont0spin=min.UserState().Value("spin0");
theFitResult.cont1spin=min.UserState().Value("spin1");
theFitResult.cont2spin=min.UserState().Value("spin2");
etacToa1320pi0fit->fillFitHists(theFitResult);
delete etacToa1320pi0fit;
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment