diff --git a/Examples/JpsiGamKsKlKK/Jamfile b/Examples/JpsiGamKsKlKK/Jamfile index 5aaad9d466b6e9886685ba67ac674ef7fb49a674..768a4245d34f4cfffa8e650e3dd66d0fbdc7df26 100644 --- a/Examples/JpsiGamKsKlKK/Jamfile +++ b/Examples/JpsiGamKsKlKK/Jamfile @@ -12,5 +12,6 @@ lib JpsiGamKlKsKK : exe MJpsiGamKlKsKKApp : MJpsiGamKlKsKKApp.cc JpsiGamKlKsKK : ; exe MJpsiGamKsKlKKParserTestApp : JpsiGamKsKlKKParserTestApp.cc JpsiGamKlKsKK : ; - +exe MMassIndepResultApp : MMassIndepResultApp.cc JpsiGamKlKsKK : ; +exe MTestFcnApp : MTestFcnApp.cc JpsiGamKlKsKK : ; diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.cc b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.cc index 2bacc00bd7b1e930c6812b40a364bb8240d2d0fa..2e60161cfa3c55337b528087419219119fd2e20e 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.cc +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.cc @@ -141,6 +141,16 @@ void JpsiGamKsKlKKEventList::read4Vecs(EventList& evtList, std::vector<EvtData*> } } + //WignerD function for 1+ -> phi phi + Spin jAxialVec =1; + for(Spin M=-jAxialVec; M<=jAxialVec; M++){ + for (Spin lam=-1; lam<=1; lam++){ + evtData->WignerDs[enumJpsiGamKsKlKKData::Df_Spin2][jAxialVec][M][lam] = Wigner_D(V4_KsKl_HeliKsKlKpKm.Phi(),V4_KsKl_HeliKsKlKpKm.Theta(), 0,jAxialVec,M,lam); //use Df_Spin2 (!) to be able to use same function for fJ(J=1,2) decay amplitudes + } + } + + + //WignerD function for 0-/0+ -> phi phi Spin jScalar =0; Spin M =0; diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh index 735258c1acb00af6f92ba38f3d31cad9fcc560bd..67852b5b7cf5dd4a3ad9ae690afbecf8c53d6b22 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh @@ -22,6 +22,9 @@ #include "Minuit2/MnUserParameters.h" + + + // using namespace std; using namespace ROOT::Minuit2; diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.cc b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.cc index 0d66542e46b57c5cff1abce12c5bbb90ef51736e..52d2e8d8dbd635beb3b882a2c537b8f1324964ae 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.cc +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.cc @@ -2,10 +2,13 @@ #include <fstream> #include <sstream> #include <string> + #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh" #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.hh" #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKData.hh" +#include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh" #include "PwaUtils/KinUtils.hh" +#include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh" #include "TFile.h" #include "TH1F.h" @@ -15,74 +18,79 @@ #include "TLorentzVector.h" #include "ErrLogger/ErrLogger.hh" -JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<const EvtDataBaseList> theEvtList) : - _theTFile(0), - _PhiPhiMassDataHist(0), - _PhiPhiMassMcHist(0), - _PhiPhiMassFittedHist(0), - _KpKmMassDataHist(0), - _KpKmMassMcHist(0), - _KpKmMassFittedHist(0), - _KsKlMassDataHist(0), - _KsKlMassMcHist(0), - _KsKlMassFittedHist(0), - _costPhi_KpKmDataHist(0), - _costPhi_KpKmMcHist(0), - _costPhi_KpKmFittedHist(0), - _phiPhi_KpKmDataHist(0), - _phiPhi_KpKmMcHist(0), - _phiPhi_KpKmFittedHist(0), - _chiDataHist(0), - _chiMcHist(0), - _chiFittedHist(0), - _dataTuple(0), - _mcTuple(0) -{ - if(0==theEvtList){ - Alert <<"JpsiGamKsKlKKEventList* theEvtList is a 0 pointer !!!!" ; // << endmsg; - exit(1); - } - - initRootStuff(); - - const std::vector<EvtData*> dataList=theEvtList->getDataVecs(); - std::vector<EvtData*>::const_iterator it=dataList.begin(); - while(it!=dataList.end()) - { - plotDalitz(_dalitzDataHist, (*it), 1.); - plotPhiPhi(_PhiPhiMassDataHist, (*it), 1. ); - plotKsKl( _KsKlMassDataHist, (*it), 1. ); - plotKpKm( _KpKmMassDataHist, (*it), 1. ); - plotCostPhiKs( _costKs_KsKlHeliDataHist, _phiKs_KsKlHeliDataHist,(*it), 1. ); - plotCostPhiKp( _costKp_KpKmHeliDataHist, _phiKp_KpKmHeliDataHist,(*it), 1. ); - plotCostGam( _costGamCmDataHist,(*it), 1. ); - plotCostPhi_PhiPhiHeli(_costPhi_KpKmDataHist, _phiPhi_KpKmDataHist, (*it)->FourVecs[enumJpsiGamKsKlKKData::V4_KpKm_HeliKsKlKpKm], 1.); - fillTuple(_dataTuple, (*it), 1.); +// JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<const EvtDataBaseList> theEvtList) : +// _theTFile(0), +// _PhiPhiMassDataHist(0), +// _PhiPhiMassMcHist(0), +// _PhiPhiMassFittedHist(0), +// _KpKmMassDataHist(0), +// _KpKmMassMcHist(0), +// _KpKmMassFittedHist(0), +// _KsKlMassDataHist(0), +// _KsKlMassMcHist(0), +// _KsKlMassFittedHist(0), +// _costPhi_KpKmDataHist(0), +// _costPhi_KpKmMcHist(0), +// _costPhi_KpKmFittedHist(0), +// _phiPhi_KpKmDataHist(0), +// _phiPhi_KpKmMcHist(0), +// _phiPhi_KpKmFittedHist(0), +// _chiDataHist(0), +// _chiMcHist(0), +// _chiFittedHist(0), +// _dataTuple(0), +// _mcTuple(0), +// _massIndepTuple(0), +// _massRange(make_pair(0,100)) +// { +// if(0==theEvtList){ +// Alert <<"JpsiGamKsKlKKEventList* theEvtList is a 0 pointer !!!!" ; // << endmsg; +// exit(1); +// } + +// initRootStuff(); + +// const std::vector<EvtData*> dataList=theEvtList->getDataVecs(); + +// std::vector<EvtData*>::const_iterator it=dataList.begin(); +// while(it!=dataList.end()) +// { +// plotDalitz(_dalitzDataHist, (*it), 1.); +// plotPhiPhi(_PhiPhiMassDataHist, (*it), 1. ); +// plotKsKl( _KsKlMassDataHist, (*it), 1. ); +// plotKpKm( _KpKmMassDataHist, (*it), 1. ); +// plotCostPhiKs( _costKs_KsKlHeliDataHist, _phiKs_KsKlHeliDataHist,(*it), 1. ); +// plotCostPhiKp( _costKp_KpKmHeliDataHist, _phiKp_KpKmHeliDataHist,(*it), 1. ); +// plotCostGam( _costGamCmDataHist,(*it), 1. ); +// plotCostPhi_PhiPhiHeli(_costPhi_KpKmDataHist, _phiPhi_KpKmDataHist, (*it)->FourVecs[enumJpsiGamKsKlKKData::V4_KpKm_HeliKsKlKpKm], 1.); +// fillTuple(_dataTuple, (*it), 1.); - ++it; - } - - const std::vector<EvtData*> mcList=theEvtList->getMcVecs(); - it=mcList.begin(); - while(it!=mcList.end()) - { - plotDalitz(_dalitzMcHist, (*it), 1.); - plotDalitz(_dalitzMcHist, (*it), 1.); - plotPhiPhi(_PhiPhiMassMcHist, (*it), 1. ); - plotKsKl( _KsKlMassMcHist, (*it), 1. ); - plotKpKm( _KpKmMassMcHist, (*it), 1. ); - plotCostPhiKs( _costKs_KsKlHeliMcHist, _phiKs_KsKlHeliMcHist,(*it), 1. ); - plotCostPhiKp( _costKp_KpKmHeliMcHist, _phiKp_KpKmHeliMcHist, (*it), 1. ); - plotCostGam( _costGamCmMcHist, (*it), 1. ); - plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecs[enumJpsiGamKsKlKKData::V4_KpKm_HeliKsKlKpKm], 1.); - fillTuple(_mcTuple, (*it), 1.); - - ++it; - } -} - -JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh> theJpsiGamKsKlKKLh, fitParams& fitParam) : +// ++it; +// } + +// const std::vector<EvtData*> mcList=theEvtList->getMcVecs(); +// it=mcList.begin(); +// while(it!=mcList.end()) +// { +// plotDalitz(_dalitzMcHist, (*it), 1.); +// plotDalitz(_dalitzMcHist, (*it), 1.); +// plotPhiPhi(_PhiPhiMassMcHist, (*it), 1. ); +// plotKsKl( _KsKlMassMcHist, (*it), 1. ); +// plotKpKm( _KpKmMassMcHist, (*it), 1. ); +// plotCostPhiKs( _costKs_KsKlHeliMcHist, _phiKs_KsKlHeliMcHist,(*it), 1. ); +// plotCostPhiKp( _costKp_KpKmHeliMcHist, _phiKp_KpKmHeliMcHist, (*it), 1. ); +// plotCostGam( _costGamCmMcHist, (*it), 1. ); +// plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecs[enumJpsiGamKsKlKKData::V4_KpKm_HeliKsKlKpKm], 1.); +// fillTuple(_mcTuple, (*it), 1.); + +// ++it; +// } + +// } + + +JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(JpsiGamKsKlKKProdLh* theJpsiGamKsKlKKLh, fitParams& fitParam, FitParamErrorMatrix* theErrorMatrix) : _theTFile(0), _dalitzDataHist(0), _dalitzMcHist(0), @@ -106,16 +114,27 @@ JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh> theJpsiGamKsKlKKLh _chiMcHist(0), _chiFittedHist(0), _dataTuple(0), - _mcTuple(0) + _mcTuple(0), + _massIndepTuple(0), + _massRange(make_pair(0.,0.) ) { if(0==theJpsiGamKsKlKKLh){ Alert <<"JpsiGamKsKlKKLh* theJpsiGamKsKlKKLh is a 0 pointer !!!!" ; // << endmsg; exit(1); } - + initRootStuff(); + _theJpsiGamKsKlKKLh = theJpsiGamKsKlKKLh; + _fitParam = fitParam; + std::vector<double> data; + _errMatrix = theErrorMatrix; + +} + + +void JpsiGamKsKlKKHist::fill(){ - boost::shared_ptr<const EvtDataBaseList> theEvtList=theJpsiGamKsKlKKLh->getEventList(); + boost::shared_ptr<const EvtDataBaseList> theEvtList=_theJpsiGamKsKlKKLh->getEventList(); const std::vector<EvtData*> dataList=theEvtList->getDataVecs(); std::vector<EvtData*>::const_iterator it=dataList.begin(); @@ -149,7 +168,7 @@ JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh> theJpsiGamKsKlKKLh plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecs[enumJpsiGamKsKlKKData::V4_KpKm_HeliKsKlKpKm], 1.); plotChi(_chiMcHist, (*it), 1. ); - double evtWeight= theJpsiGamKsKlKKLh->calcEvtIntensity((*it), fitParam); + double evtWeight= _theJpsiGamKsKlKKLh->calcEvtIntensity((*it), _fitParam); plotDalitz(_dalitzFittedHist, (*it),evtWeight ); plotPhiPhi(_PhiPhiMassFittedHist, (*it), evtWeight ); plotKsKl( _KsKlMassFittedHist, (*it), evtWeight ); @@ -187,6 +206,50 @@ JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh> theJpsiGamKsKlKKLh _costPhi_KpKmFittedHist->Scale(scaleFactor); _phiPhi_KpKmFittedHist->Scale(scaleFactor); _chiFittedHist->Scale(scaleFactor); + + + double iEta(0.), iEtaErr(0.), iF0(0.), iF0Err(0.), iEta2(0.), iEta2Err(0.), iF1(0.), iF1Err(0.), iF2(0.), iF2Err(0.); + + + it=mcList.begin(); + while(it!=mcList.end()){ + + std::pair<double, double> intensityEvent = make_pair(0.,0.); + _theJpsiGamKsKlKKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "etacHyp", intensityEvent); + iEta+= intensityEvent.first*scaleFactor; + iEtaErr+= intensityEvent.second*scaleFactor; + + _theJpsiGamKsKlKKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f02020Hyp",intensityEvent); + iF0 += intensityEvent.first*scaleFactor; + iF0Err += intensityEvent.second*scaleFactor; + + _theJpsiGamKsKlKKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f1Hyp", intensityEvent); + iF1 += intensityEvent.first*scaleFactor; + iF1Err += intensityEvent.second*scaleFactor; + + _theJpsiGamKsKlKKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f22300Hyp",intensityEvent); + iF2 += intensityEvent.first*scaleFactor; + iF2Err += intensityEvent.second*scaleFactor; + + _theJpsiGamKsKlKKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "eta21870Hyp",intensityEvent); + iEta2+= intensityEvent.first*scaleFactor; + iEta2Err += intensityEvent.second*scaleFactor; + + it++; + } + double meanMassRange = _massRange.first + 0.5*(_massRange.second-_massRange.first); + + + Info << "Events for component eta : " << iEta << " +/- " << iEtaErr ; + Info << "Events for component f0: " << iF0 << " +/- " << iF0Err ; + Info << "Events for component f1: " << iF1 << " +/- " << iF1Err ; + Info << "Events for component f2: " << iF2 << " +/- " << iF2Err ; + Info << "Events for component eta2: " << iEta2 << " +/- " << iEta2Err ; + + + + _massIndepTuple->Fill(meanMassRange, iEta, iEtaErr, iF0, iF0Err, iF1, iF1Err, iF2, iF2Err, iEta2, iEta2Err ); + } JpsiGamKsKlKKHist::~JpsiGamKsKlKKHist() @@ -265,6 +328,8 @@ void JpsiGamKsKlKKHist::initRootStuff() _dataTuple=new TNtuple("_dataTuple", "data ntuple", tupleVariables.data()); _mcTuple=new TNtuple("_mcTuple", "mc ntuple", tupleVariables.data()); + + _massIndepTuple = new TNtuple("_massIndepTuple","results from mass indep. fits", ("mass:eta:etaErr:f0:f0Err:f1:f1Err:f2:f2Err:eta2:eta2Err") ); } void JpsiGamKsKlKKHist::plotDalitz(TH2F* theHisto, EvtData* theData, double weight) @@ -430,3 +495,5 @@ double JpsiGamKsKlKKHist::decayAngleChi(const Vector4<double>& v4_p,const Vector } + + diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh index 08b06625477a268829d5e2a8f6bafff64a44f0cf..ec94a19c1a13a5065a2cf6d09a32704494ba4d25 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh @@ -5,6 +5,8 @@ #include <fstream> #include <string> #include <vector> +#include <utility> + #include <cassert> @@ -17,13 +19,15 @@ #include "PwaUtils/EvtDataBaseList.hh" #include "PwaUtils/FitParamsBase.hh" #include "PwaUtils/AbsLh.hh" - +#include "Examples/JpsiGamKsKlKK/FitParamErrorMatrix.hh" class TFile; class TH2F; class TH1F; class TNtuple; - +class JpsiGamKsKlKKProdLh; +class EvtDataBaseList; +class FitParamErrorMatrix; class JpsiGamKsKlKKHist { @@ -33,7 +37,11 @@ public: ///Constructor JpsiGamKsKlKKHist(boost::shared_ptr<const EvtDataBaseList>); - JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh>, fitParams&); + //JpsiGamKsKlKKHist(boost::shared_ptr<AbsLh>, fitParams&); + JpsiGamKsKlKKHist(JpsiGamKsKlKKProdLh* theJpsiGamKsKlKKLh, fitParams& fitParam, FitParamErrorMatrix* theErrorMatrix ); + void fill(); + void setMassRange(std::pair<double, double> theMassRange){ _massRange = theMassRange; } + /** Destructor */ virtual ~JpsiGamKsKlKKHist(); @@ -45,7 +53,7 @@ protected: private: - + TFile* _theTFile; @@ -97,7 +105,8 @@ private: TNtuple* _dataTuple; TNtuple* _mcTuple; - + TNtuple* _massIndepTuple; + std::pair<double, double> _massRange; void initRootStuff(); void plotDalitz(TH2F* theHisto, EvtData* theData, double weight); @@ -118,6 +127,10 @@ private: const Vector4<double>& v4_h2 ) ; + JpsiGamKsKlKKProdLh* _theJpsiGamKsKlKKLh; + fitParams _fitParam; + FitParamErrorMatrix* _errMatrix; + }; #endif diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.cc b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.cc index 549d25ea113c77d364bfb109d190e80dd8540229..f9a7f2a95b177e15fdb3d2f2b951677691541540 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.cc +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.cc @@ -5,8 +5,16 @@ #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh" #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.hh" #include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKFitParams.hh" +#include "Examples/JpsiGamKsKlKK/FitParamErrorMatrix.hh" +#include "Examples/JpsiGamKsKlKK/FitParamIndex.hh" + #include "ErrLogger/ErrLogger.hh" +#include <boost/bind.hpp> +#include <boost/numeric/ublas/matrix.hpp> +#include <boost/numeric/ublas/io.hpp> + + JpsiGamKsKlKKProdLh::JpsiGamKsKlKKProdLh(boost::shared_ptr<const EvtDataBaseList> theEvtList, const std::map<const std::string, bool>& hypMap) : AbsLh(theEvtList) ,_etacHyp(false) @@ -22,13 +30,8 @@ JpsiGamKsKlKKProdLh::JpsiGamKsKlKKProdLh(boost::shared_ptr<const EvtDataBaseList ,_phiMass( 1.019455) ,_kaonMass(0.493677) { - - initializeHypothesisMap( hypMap); - - - - + } JpsiGamKsKlKKProdLh::JpsiGamKsKlKKProdLh( boost::shared_ptr<AbsLh> theLhPtr, const std::map<const std::string, bool>& hypMap ) : @@ -118,7 +121,6 @@ complex<double> JpsiGamKsKlKKProdLh::etaGammaAmp(Spin Minit, Spin Metac, Spin Mg boost::shared_ptr<const JPCLS> PsiState=itPsi->first; double thePsiMag=itPsi->second; double thePsiPhi=PsiToEtaPhi[PsiState]; - complex<double> expiphiPsi(cos(thePsiPhi), sin(thePsiPhi)); Spin lambda = Metac-Mgamma; if( fabs(lambda)>PsiState->J || fabs(lambda)>PsiState->S) continue; @@ -602,24 +604,29 @@ JpsiGamKsKlKKProdLh::initializeHypothesisMap( const std::map<const std::string, } - -double -JpsiGamKsKlKKProdLh::calcComponentIntensity( EvtData* theData, fitParams& theParamVal, std::string component ){ +bool +JpsiGamKsKlKKProdLh::calcComponentIntensity( EvtData* theData, fitParams& theParamVal, FitParamErrorMatrix& theErrMatrix, std::string component, std::pair<double, double> &intensity ){ - double result=0.0; complex<double> JmpGmp(0.0,0.0); complex<double> JmpGmm(0.0,0.0); complex<double> JmmGmp(0.0,0.0); complex<double> JmmGmm(0.0,0.0); + intensity = make_pair(0.,0.); + + std::map<const std::string, bool>::const_iterator iter= _hypMap.find(component); if (iter ==_hypMap.end()){ - Alert << "Component " << component << " was not included in fit hypothesis" <<endmsg ; - return 0; + static int alerts=0; + if(alerts<5) Alert << "Component " << component << " was not included in fit hypothesis"; + alerts++; + return true; } if (iter !=_hypMap.end() && !iter->second ){ - Alert << "Component " << component << " was disabled in fit hypothesis" <<endmsg ; - return 0; + static int alerts=0; + if(alerts<5) Alert << "Component " << component << " was disabled in fit hypothesis" ; + alerts++; + return true; } @@ -627,22 +634,69 @@ JpsiGamKsKlKKProdLh::calcComponentIntensity( EvtData* theData, fitParams& thePa if(component == "etacHyp" ){ calcEtacGammaAmp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm,dynamicModelParams::MassIndependent ); + intensity.first =norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + + boost::function<void(EvtData* , fitParams&, complex<double>&, complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel) > ampFunction = boost::bind(&JpsiGamKsKlKKProdLh::calcEtacGammaAmp, this, _1, _2, _3, _4, _5, _6,_7 ); + double intensityError=0.0; + std::vector< int > theAmps; + theAmps.push_back( paramEnumJpsiGamKsKlKK::PsiToEtacGamma ); + calcAmpError( theData, theParamVal, theErrMatrix, dynamicModelParams::MassIndependent, ampFunction, theAmps, intensityError ); + intensity.second=intensityError; } + + else if( component == "f02020Hyp" ){ calcF02020Amp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::MassIndependent ); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + boost::function<void(EvtData* , fitParams&, complex<double>&, complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel) > ampFunction = boost::bind(&JpsiGamKsKlKKProdLh::calcF02020Amp, this, _1, _2, _3, _4, _5, _6,_7 ); + double intensityError=0.0; + std::vector< int > theAmps; + theAmps.push_back( paramEnumJpsiGamKsKlKK::PsiToF02020Gamma ); + theAmps.push_back( paramEnumJpsiGamKsKlKK::F02020ToPhiPhi ); + calcAmpError( theData, theParamVal, theErrMatrix, dynamicModelParams::MassIndependent, ampFunction, theAmps, intensityError ); + intensity.second=intensityError; } + else if(component == "f22300Hyp"){ calcF22300Amp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::MassIndependent); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + boost::function<void(EvtData* , fitParams&, complex<double>&, complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel) > ampFunction = boost::bind(&JpsiGamKsKlKKProdLh::calcF22300Amp, this, _1, _2, _3, _4, _5, _6,_7 ); + double intensityError=0.0; + std::vector< int > theAmps; + theAmps.push_back( paramEnumJpsiGamKsKlKK::PsiToF22300Gamma ); + theAmps.push_back( paramEnumJpsiGamKsKlKK::F22300ToPhiPhi ); + calcAmpError( theData, theParamVal, theErrMatrix, dynamicModelParams::MassIndependent, ampFunction, theAmps, intensityError ); + intensity.second=intensityError; } + else if(component == "eta21870Hyp"){ calcE21870Amp(theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::MassIndependent); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + boost::function<void(EvtData* , fitParams&, complex<double>&, complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel) > ampFunction = boost::bind(&JpsiGamKsKlKKProdLh::calcE21870Amp, this, _1, _2, _3, _4, _5, _6,_7 ); + double intensityError=0.0; + std::vector< int > theAmps; + theAmps.push_back( paramEnumJpsiGamKsKlKK::PsiToEta21870Gamma ); + theAmps.push_back( paramEnumJpsiGamKsKlKK::Eta21870ToPhiPhi ); + calcAmpError( theData, theParamVal, theErrMatrix, dynamicModelParams::MassIndependent, ampFunction, theAmps, intensityError ); + intensity.second=intensityError; + } + else if(component == "f1Hyp"){ calcF1Amp(theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::MassIndependent); - } + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); - result=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); - return result; + boost::function<void(EvtData* , fitParams&, complex<double>&, complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel) > ampFunction = boost::bind(&JpsiGamKsKlKKProdLh::calcF1Amp, this, _1, _2, _3, _4, _5, _6,_7 ); + double intensityError=0.0; + std::vector< int > theAmps; + theAmps.push_back( paramEnumJpsiGamKsKlKK::PsiToF1Gamma ); + theAmps.push_back( paramEnumJpsiGamKsKlKK::F1ToPhiPhi ); + calcAmpError( theData, theParamVal, theErrMatrix, dynamicModelParams::MassIndependent, ampFunction, theAmps, intensityError ); + intensity.second=intensityError; + + + } + return true; } @@ -653,39 +707,41 @@ JpsiGamKsKlKKProdLh::calcComponentIntensity( EvtData* theData, fitParams& thePa if(component == "etacHyp" ){ calcEtacGammaAmp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::BreitWigner ); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } else if(component == "eta2225Hyp"){ calcEta2225Amp(theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm,dynamicModelParams::BreitWigner ); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } else if( component == "f02020Hyp" ){ calcF02020Amp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::BreitWigner ); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } else if(component == "f02020FlatteHyp"){ calcF02020Amp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::Flatte ); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } else if(component == "f22300Hyp"){ calcF22300Amp( theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::BreitWigner); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } else if(component == "eta21870Hyp"){ calcE21870Amp(theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::BreitWigner); - } + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); +} else if(component == "f1Hyp"){ calcF1Amp(theData, theParamVal, JmpGmp,JmpGmm, JmmGmp,JmmGmm, dynamicModelParams::BreitWigner); + intensity.first=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); } - return norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); - - - - - return 0.0; + return true; } @@ -870,3 +926,111 @@ JpsiGamKsKlKKProdLh::calcE21870Amp(EvtData* theData, fitParams& theParamVal, JmmGmp+=f2GammaAmp(-1, 1, f1Spin, theData, PsiTof1GamMag, PsiTof1GamPhi,F1ToPhiPhiMag,F1ToPhiPhiPhi,mass,width, _useCommonProductionPhase, theDynModel ); JmmGmm+=f2GammaAmp(-1, -1, f1Spin, theData, PsiTof1GamMag, PsiTof1GamPhi,F1ToPhiPhiMag,F1ToPhiPhiPhi,mass,width, _useCommonProductionPhase, theDynModel ); } + + + + void + JpsiGamKsKlKKProdLh::dumpComponentIntensity( std::ostream &os, fitParams& theParams, FitParamErrorMatrix &theErrMatrix ){ + + const std::vector<EvtData*> mcList=getEventList()->getMcVecs(); + const std::vector<EvtData*> dataList=getEventList()->getDataVecs(); + + std::map<const std::string, bool>::const_iterator hypo= _hypMap.begin(); + while(hypo !=_hypMap.end()){ + if( hypo->second ){ + + std::vector<EvtData*>::const_iterator it=mcList.begin(); + double integral=0.0; + double integralError=0.0; + while(it!=mcList.end()){ + + std::pair<double, double> intensityEvent=make_pair(0.,0.); + calcComponentIntensity( *it, theParams, theErrMatrix, hypo->first, intensityEvent ); + integral+=intensityEvent.first; + integralError+=intensityEvent.second; + it++; + } + Info << "Unscaled intensity for component " << hypo->first << ": " << integral << " +/- " << integralError << endmsg; + double scale = (1.*dataList.size()) / (1.*mcList.size()) ; + Info << "Scale factor: " << scale << endmsg; + integral*=scale; + integralError*=scale; + Info << "Events for component " << hypo->first << ": " << integral << " +/- " << integralError << endmsg; + + os << hypo->first << "\t" << integral << endl; + } + hypo++; + } + } + + + + + +void JpsiGamKsKlKKProdLh::calcAmpError( EvtData* theData, fitParams& theParamVal, FitParamErrorMatrix& theErrMatrix, + dynamicModelParams::enumDynamicModel theModel, + boost::function<void(EvtData* , fitParams&, complex<double>&, + complex<double>&, complex<double>&, complex<double>&, + dynamicModelParams::enumDynamicModel)> calcAmp, + std::vector< int > theAmpsEnum, + double& theAmpError) +{ + + double epsilon=0.5e-5; + double error(0.); + std::vector<double> derivatives; + std::vector<int> parameterIndices; + + std::vector<int>::iterator ampIter; + for( ampIter=theAmpsEnum.begin(); ampIter!=theAmpsEnum.end(); ampIter++){ + + + std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > theMag = theParamVal.Mags[*ampIter]; + + std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::iterator it; + for ( it=theMag.begin(); it!=theMag.end(); ++it){ + + fitParams newParamsPlus = theParamVal; + fitParams newParamsMinus = theParamVal; + boost::shared_ptr<const JPCLS> theState=it->first; + newParamsPlus.Mags[*ampIter][theState]=theParamVal.Mags[*ampIter][theState]+epsilon; + newParamsMinus.Mags[*ampIter][theState]=theParamVal.Mags[*ampIter][theState]-epsilon; + + complex<double> JmpGmp(0.0,0.0); + complex<double> JmpGmm(0.0,0.0); + complex<double> JmmGmp(0.0,0.0); + complex<double> JmmGmm(0.0,0.0); + + calcAmp( theData, newParamsPlus, JmpGmp,JmpGmm, JmmGmp,JmmGmm,dynamicModelParams::MassIndependent ); + double resultPlus=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + + JmpGmp = complex<double>(0.0,0.0); + JmpGmm = complex<double>(0.0,0.0); + JmmGmp = complex<double>(0.0,0.0); + JmmGmm = complex<double>(0.0,0.0); + + calcAmp( theData, newParamsMinus, JmpGmp,JmpGmm, JmmGmp,JmmGmm,dynamicModelParams::MassIndependent ); + double resultMinus=norm(JmpGmp)+norm(JmpGmm)+norm(JmmGmp)+norm(JmmGmm); + + derivatives.push_back((resultPlus-resultMinus)/(2*epsilon)); + static FitParamIndex theParamIndex(theParamVal); + parameterIndices.push_back( theParamIndex.Mag(*ampIter, theState ) ); + } + + int elements = derivatives.size(); + for(int i=0; i<elements; i++){ + for(int j=i;j<elements;j++){ + if(i==j){ + error+=derivatives[i]*derivatives[j]*theErrMatrix(parameterIndices[i],parameterIndices[j] ); + }else{ + error+=derivatives[i]*derivatives[j]*theErrMatrix(parameterIndices[i],parameterIndices[j] ); + } + } + } + } + theAmpError=sqrt(error); + +} + + + diff --git a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh index e40f3ed5b89aaa162d04a61eec3c5c22bc8ae9b5..ee92d49ce2fb5e5e9030626436438b1ad05ac7e1 100644 --- a/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh +++ b/Examples/JpsiGamKsKlKK/JpsiGamKsKlKKProdLh.hh @@ -9,6 +9,7 @@ #include <cassert> #include <boost/shared_ptr.hpp> +#include <boost/function.hpp> #include "TROOT.h" // #include <TSystem.h> @@ -22,7 +23,7 @@ #include "Minuit2/MnUserParameters.h" - +class FitParamErrorMatrix; class JpsiGamKsKlKKProdLh : public AbsLh{ @@ -41,8 +42,6 @@ public: return new JpsiGamKsKlKKProdLh(_evtListPtr, _hypMap); } - - virtual double calcEvtIntensity( EvtData* theData, fitParams& theParamVal); //Getters: @@ -52,8 +51,9 @@ public: void useCommonProductionPhase( bool commonPhase ){_useCommonProductionPhase=commonPhase;} void massIndependentFit( bool massIndep){ _massIndependentFit=massIndep; } - double calcComponentIntensity( EvtData* theData, fitParams& theParamVal, std::string component ); - + bool calcComponentIntensity( EvtData* theData, fitParams& theParamVal, FitParamErrorMatrix& theErrMatrix, std::string component, std::pair<double, double> &intensity ); + void dumpComponentIntensity( std::ostream &os, fitParams& theParams, FitParamErrorMatrix& theErrMatrix ); + protected: @@ -113,6 +113,8 @@ protected: void calcF1Amp(EvtData* theData, fitParams& theParamVal, complex<double> &JmpGmp, complex<double> &JmpGmm, complex<double> &JmmGmp, complex<double> &JmmGmm, dynamicModelParams::enumDynamicModel theModel); + + bool _etacHyp; bool _eta2225Hyp; @@ -134,6 +136,13 @@ protected: private: + void calcAmpError( EvtData* theData, fitParams& theParamVal, FitParamErrorMatrix& theErrMatrix, + dynamicModelParams::enumDynamicModel theModel, + boost::function<void(EvtData* , fitParams&, complex<double>&, + complex<double>&, complex<double>&, complex<double>&, dynamicModelParams::enumDynamicModel)>, + std::vector< int >, double& ); + + bool initializeHypothesisMap( const std::map<const std::string, bool>& hypMap ); const double _phiMass;// = 1.019455; @@ -143,3 +152,5 @@ private: }; #endif + + diff --git a/Examples/JpsiGamKsKlKK/MJpsiGamKlKsKKApp.cc b/Examples/JpsiGamKsKlKK/MJpsiGamKlKsKKApp.cc index 2669fd8ff821b09ca916e57ef40c14f392d95057..a057df93f44dc3092fb6edb9434ca55942e8414c 100644 --- a/Examples/JpsiGamKsKlKK/MJpsiGamKlKsKKApp.cc +++ b/Examples/JpsiGamKsKlKK/MJpsiGamKlKsKKApp.cc @@ -37,7 +37,12 @@ #include "Minuit2/FunctionMinimum.h" #include "Minuit2/MnMinos.h" #include "Minuit2/MnStrategy.h" +#include "Minuit2/MnPrint.h" +#include "Minuit2/MnScan.h" +#include "Examples/JpsiGamKsKlKK/FitParamErrorMatrix.hh" +#include "Examples/JpsiGamKsKlKK/FitParamIndex.hh" +#include "Examples/JpsiGamKsKlKK/FitParamErrorMatrixStreamer.hh" using namespace ROOT::Minuit2; @@ -171,23 +176,21 @@ int main(int __argc,char *__argv[]){ } boost::shared_ptr<AbsLh> theLhPtr; + JpsiGamKsKlKKProdLh* theProdLh = new JpsiGamKsKlKKProdLh(theJpsiGamKsKlKKEventListPtr, hypMap); + theProdLh->massIndependentFit( theAppParams.massIndependentFit() ); + theProdLh->useCommonProductionPhase( theAppParams.useCommonProductionPhases() ); + std::string startWithHyp=theAppParams.startHypo(); if (startWithHyp=="production"){ - JpsiGamKsKlKKProdLh* tmplh = new JpsiGamKsKlKKProdLh(theJpsiGamKsKlKKEventListPtr, hypMap); - tmplh->massIndependentFit( theAppParams.massIndependentFit() ); - tmplh->useCommonProductionPhase( theAppParams.useCommonProductionPhases() ); - theLhPtr = boost::shared_ptr<AbsLh> (tmplh); - //theLhPtr = boost::shared_ptr<AbsLh> (new JpsiGamKsKlKKProdLh(theJpsiGamKsKlKKEventListPtr, hypMap )); + theLhPtr = boost::shared_ptr<AbsLh> (theProdLh); } else { Alert << "start with hypothesis " << startWithHyp << " not supported!!!!" ; // << endmsg; exit(1); } - - - + std::string mode=theAppParams.mode(); if (mode=="dumpDefaultParams"){ fitParams defaultVal; @@ -220,37 +223,25 @@ int main(int __argc,char *__argv[]){ double theLh=theLhPtr->calcLogLh(theStartparams); Info <<"theLh = "<< theLh << endmsg; - JpsiGamKsKlKKHist theHist(theLhPtr, theStartparams); + + std::string errFile = "finalErrorMatrix.dat"; + FitParamErrorMatrixStreamer theErrStreamer( errFile ); + std::vector<double> theErrData; + int ncols(0); + theErrStreamer.matrixData( theErrData, ncols ); + FitParamErrorMatrix theErrorMatrix(theErrData, ncols ); + + JpsiGamKsKlKKHist theHist(theProdLh, theStartparams, &theErrorMatrix); + theHist.setMassRange(theAppParams.massRange() ); + theHist.fill(); if(theAppParams.massIndependentFit()){ //calculate intensity contributions - JpsiGamKsKlKKProdLh* contrLh = new JpsiGamKsKlKKProdLh(theJpsiGamKsKlKKEventListPtr, hypMap); - contrLh->massIndependentFit( theAppParams.massIndependentFit() ); - contrLh->useCommonProductionPhase( theAppParams.useCommonProductionPhases() ); - - boost::shared_ptr<const EvtDataBaseList> theEvtList=contrLh->getEventList(); - const std::vector<EvtData*> mcList=theEvtList->getMcVecs(); - const std::vector<EvtData*> dataList=theEvtList->getDataVecs(); - - std::map<const std::string, bool>::const_iterator hypo= hypMap.begin(); - while(hypo !=hypMap.end()){ - if( hypo->second ){ - - std::vector<EvtData*>::const_iterator it=mcList.begin(); - double integral=0.0; - while(it!=mcList.end()){ - integral+=contrLh->calcComponentIntensity( *it, theStartparams, hypo->first ); - it++; - } - double scale = (1.*dataList.size()) / (1.*mcList.size()) ; - Info << "Scale factor: " << scale << endmsg; - integral*=scale; - Info << "Events for component " << hypo->first << ": " << integral << endmsg; - } - hypo++; - } + //std::ofstream theStream ( "componentIntensity.dat"); + //theProdLh->dumpComponentIntensity( theStream, theStartparams, theErrorMatrix ); } - + + end= clock(); double cpuTime= (end-start)/ (CLOCKS_PER_SEC); Info << "cpuTime:\t" << cpuTime << "\tsec" << endmsg; @@ -274,12 +265,15 @@ int main(int __argc,char *__argv[]){ for (itFix=fixedParams.begin(); itFix!=fixedParams.end(); ++itFix){ upar.Fix( (*itFix) ); } - + + //MnScan theScan(theFcn, upar); + //theScan(); + //cout << upar << endl; MnMigrad migrad(theFcn, upar); Info <<"start migrad "<< endmsg; FunctionMinimum min = migrad(); - + if(!min.IsValid()) { //try with higher strategy Info <<"FM is invalid, try with strategy = 2."<< endmsg; @@ -289,11 +283,20 @@ int main(int __argc,char *__argv[]){ MnUserParameters finalUsrParameters=min.UserParameters(); const std::vector<double> finalParamVec=finalUsrParameters.Params(); - + fitParams finalFitParams=theFitParamBase->getFitParamVal(finalParamVec); - JpsiGamKsKlKKHist theHist(theLhPtr, finalFitParams); + //MnUserCovariance theCov = min.UserCovariance() ; + //cout << "User vov : "<< endl; + //cout << theCov << endl; + + MnUserParameterState theState = min.UserState(); + cout << "User state " << theState << endl; + + + + theFitParamBase->printParams(finalFitParams); double theLh=theLhPtr->calcLogLh(finalFitParams); Info <<"theLh = "<< theLh << endmsg; @@ -307,9 +310,23 @@ int main(int __argc,char *__argv[]){ fitParams finalFitErrs=theFitParamBase->getFitParamVal(finalParamErrorVec); std::ofstream theStream ( "finalResult.dat"); theFitParamBase->dumpParams(theStream, finalFitParams, finalFitErrs); + + + + MnUserCovariance theCovMatrix = min.UserCovariance(); + std::cout << min << std::endl; + std::ofstream theErrMatStream ( "finalErrorMatrix.dat"); + FitParamErrorMatrix theErrMatrix(theCovMatrix, finalUsrParameters ); + theErrMatrix.Write(theErrMatStream); + + //std::ofstream theCompStream ( "componentIntensity.dat"); + //theProdLh->dumpComponentIntensity( theCompStream, finalFitParams, theErrMatrix ); + JpsiGamKsKlKKHist theHist(theProdLh, theStartparams, &theErrMatrix); + theHist.fill(); + return 0; } - return 0; + return 0; } diff --git a/Examples/JpsiGamKsKlKK/MassIndep.cfg b/Examples/JpsiGamKsKlKK/MassIndep.cfg index c7870a82f435844e4fa7878eadcbe0cee4bcb122..77ce777afcf01ace81cceda2146eb80cb3faa5fc 100644 --- a/Examples/JpsiGamKsKlKK/MassIndep.cfg +++ b/Examples/JpsiGamKsKlKK/MassIndep.cfg @@ -7,23 +7,31 @@ errLogMode = debug datFile = /data/sleipnir1/marc/PwaJpsiGamKsKlKK/FS4Vectors-Data-JpsiGamKsKlKpKm.dat mcFile = /data/sleipnir1/marc/PwaJpsiGamKsKlKK/FS4Vectors-PHSPMC-JpsiGamKsKlKpKm.dat -paramFile = ./startParams_MassIndependent.dat +paramFile = ./startParams.dat startHypo = production massIndependentFit = 1 -massRangeMin = 2.9 -massRangeMax = 3.1 +massRangeMin = 2.95 +massRangeMax = 3.05 enableHyp = etacHyp enableHyp = f02020Hyp -enableHyp = f22300Hyp + mnParFix = etacMass mnParFix = etacWidth mnParFix = f02020Mass mnParFix = f02020Width -mnParFix = f22300Mass -mnParFix = f22300Width + + +mnParFix = J1P-1C-1L2S1PsiToF02020GammaMag + +mnParFix = J1P-1C-1L1S1PsiToEtacGamPhi +mnParFix = J1P-1C-1L0S1PsiToF02020GammaPhi +mnParFix = J1P-1C-1L2S1PsiToF02020GammaPhi +mnParFix = J0P1C1L0S0F02020ToPhiPhiPhi +mnParFix = J0P1C1L2S2F02020ToPhiPhiPhi + mode = dumpDefaultParams