diff --git a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc index a4b0cb9a75fb574ec281637d2b4a6a295ae900e9..d1b1ed576360272480578f7b598d868962c8b237 100644 --- a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc +++ b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc @@ -333,9 +333,6 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToK0K2Amp(Psi2STo2K2PiGamData::Psi2STo2 - // result+=theMag*expiphi* - // (ampK0ToKpPi0*ampK2ToKmPi1+ampK0ToKpPi1*ampK2ToKmPi0+ampK0ToKmPi0+ampK2ToKpPi1+ampK0ToKmPi1*ampK2ToKpPi0); - result+=theMag*expiphi* (ampK0ToKpPi0*ampK2ToKmPi1+ampK0ToKpPi1*ampK2ToKmPi0+ampK0ToKmPi0*ampK2ToKpPi1+ampK0ToKmPi1*ampK2ToKpPi0); @@ -345,6 +342,106 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToK0K2Amp(Psi2STo2K2PiGamData::Psi2STo2 } +complex<double> AbsPsi2STo2K2PiGamLh::chiToK0K1Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToK0K1, double K0Mass, double K0Width, double K1Mass, double K1Width){ + complex<double> result(0.,0.); + + Vector4<double> KpPi0(theData->KpPi0_HeliChic0_4V.E(), theData->KpPi0_HeliChic0_4V.Px(), + theData->KpPi0_HeliChic0_4V.Py(), theData->KpPi0_HeliChic0_4V.Pz()); + + Vector4<double> KpPi1(theData->KpPi1_HeliChic0_4V.E(), theData->KpPi1_HeliChic0_4V.Px(), + theData->KpPi1_HeliChic0_4V.Py(), theData->KpPi1_HeliChic0_4V.Pz()); + + Vector4<double> KmPi0(theData->KmPi0_HeliChic0_4V.E(), theData->KmPi0_HeliChic0_4V.Px(), + theData->KmPi0_HeliChic0_4V.Py(), theData->KmPi0_HeliChic0_4V.Pz()); + + Vector4<double> KmPi1(theData->KmPi1_HeliChic0_4V.E(), theData->KmPi1_HeliChic0_4V.Px(), + theData->KmPi1_HeliChic0_4V.Py(), theData->KmPi1_HeliChic0_4V.Pz()); + + complex<double> ampK0ToKpPi0= BreitWigner(KpPi0, K0Mass, K0Width); + complex<double> ampK0ToKpPi1= BreitWigner(KpPi1, K0Mass, K0Width); + complex<double> ampK0ToKmPi0= BreitWigner(KmPi0, K0Mass, K0Width); + complex<double> ampK0ToKmPi1= BreitWigner(KmPi1, K0Mass, K0Width); + complex<double> ampK1ToKpPi0= sqrt(3.)*(conj(theData->DfKst1pToKpPi0[1][0][0])*BreitWignerBlattW(KpPi0, 0.493677, 0.1349766, K1Mass, K1Width, 1)); + complex<double> ampK1ToKpPi1= sqrt(3.)*(conj(theData->DfKst1pToKpPi1[1][0][0])*BreitWignerBlattW(KpPi1, 0.493677, 0.1349766, K1Mass, K1Width, 1)); + complex<double> ampK1ToKmPi0= sqrt(3.)*(conj(theData->DfKst1mToKmPi0[1][0][0])*BreitWignerBlattW(KmPi0, 0.493677, 0.1349766, K1Mass, K1Width, 1)); + complex<double> ampK1ToKmPi1= sqrt(3.)*(conj(theData->DfKst1mToKmPi1[1][0][0])*BreitWignerBlattW(KmPi1, 0.493677, 0.1349766, K1Mass, K1Width, 1)); + + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator it; + for ( it=ChiToK0K1.begin(); it!=ChiToK0K1.end(); ++it){ + + boost::shared_ptr<const JPCLS> ChiToK0K1_State=it->first; + double theMag=it->second.first; + double thePhi=it->second.second; + complex<double> expiphi(cos(thePhi), sin(thePhi)); + + + result+=theMag*expiphi*sqrt(2.*ChiToK0K1_State->L+1.)* + (ampK0ToKpPi0*ampK1ToKmPi1+ampK0ToKpPi1*ampK1ToKmPi0+ampK0ToKmPi0*ampK1ToKpPi1+ampK0ToKmPi1*ampK1ToKpPi0); + + } + + return result; + +} + + + +complex<double> AbsPsi2STo2K2PiGamLh::chiToK1K2Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToK1K2, double K1Mass, double K1Width, double K2Mass, double K2Width){ + complex<double> result(0.,0.); + + Vector4<double> KpPi0(theData->KpPi0_HeliChic0_4V.E(), theData->KpPi0_HeliChic0_4V.Px(), + theData->KpPi0_HeliChic0_4V.Py(), theData->KpPi0_HeliChic0_4V.Pz()); + + Vector4<double> KpPi1(theData->KpPi1_HeliChic0_4V.E(), theData->KpPi1_HeliChic0_4V.Px(), + theData->KpPi1_HeliChic0_4V.Py(), theData->KpPi1_HeliChic0_4V.Pz()); + + Vector4<double> KmPi0(theData->KmPi0_HeliChic0_4V.E(), theData->KmPi0_HeliChic0_4V.Px(), + theData->KmPi0_HeliChic0_4V.Py(), theData->KmPi0_HeliChic0_4V.Pz()); + + Vector4<double> KmPi1(theData->KmPi1_HeliChic0_4V.E(), theData->KmPi1_HeliChic0_4V.Px(), + theData->KmPi1_HeliChic0_4V.Py(), theData->KmPi1_HeliChic0_4V.Pz()); + + + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator it; + for ( it=ChiToK1K2.begin(); it!=ChiToK1K2.end(); ++it){ + + boost::shared_ptr<const JPCLS> ChiToK1K2_State=it->first; + double theMag=it->second.first; + double thePhi=it->second.second; + complex<double> expiphi(cos(thePhi), sin(thePhi)); + + complex<double> tmpResult(0.,0.); + for (Spin lamK1=-1; lamK1<=1; ++lamK1){ + Spin lamK2=lamK1; + Spin lamKK=lamK2-lamK1; + + complex<double> ampK1ToKpPi0= conj(theData->DfKst1pToKpPi0[1][lamK1][0])*BreitWignerBlattW(KpPi0, 0.493677, 0.1349766, K1Mass, K1Width, 1); + complex<double> ampK1ToKpPi1= conj(theData->DfKst1pToKpPi1[1][lamK1][0])*BreitWignerBlattW(KpPi1, 0.493677, 0.1349766, K1Mass, K1Width, 1); + complex<double> ampK1ToKmPi0= conj(theData->DfKst1mToKmPi0[1][lamK1][0])*BreitWignerBlattW(KmPi0, 0.493677, 0.1349766, K1Mass, K1Width, 1); + complex<double> ampK1ToKmPi1= conj(theData->DfKst1mToKmPi1[1][lamK1][0])*BreitWignerBlattW(KmPi1, 0.493677, 0.1349766, K1Mass, K1Width, 1); + + complex<double> ampK2ToKpPi0= conj(theData->DfKst2pToKpPi0[2][lamK2][0])*BreitWignerBlattW(KpPi0, 0.493677, 0.1349766, K2Mass, K2Width, 2); + complex<double> ampK2ToKpPi1= conj(theData->DfKst2pToKpPi1[2][lamK2][0])*BreitWignerBlattW(KpPi1, 0.493677, 0.1349766, K2Mass, K2Width, 2); + complex<double> ampK2ToKmPi0= conj(theData->DfKst2mToKmPi0[2][lamK2][0])*BreitWignerBlattW(KmPi0, 0.493677, 0.1349766, K2Mass, K2Width, 2); + complex<double> ampK2ToKmPi1= conj(theData->DfKst2mToKmPi1[2][lamK2][0])*BreitWignerBlattW(KmPi1, 0.493677, 0.1349766, K2Mass, K2Width, 2); + + tmpResult+=Clebsch(ChiToK1K2_State->L, 0, ChiToK1K2_State->S, lamKK, ChiToK1K2_State->J, lamKK)*Clebsch(2,lamK2, 1, -lamK1, ChiToK1K2_State->S, lamKK)* + ( + ampK1ToKpPi0*ampK2ToKmPi1 + +ampK1ToKpPi1*ampK2ToKmPi0 + +ampK1ToKmPi0*ampK2ToKpPi1 + +ampK1ToKmPi1*ampK2ToKpPi0 ); + + } + result+=theMag*expiphi*sqrt(2.*ChiToK1K2_State->L+1.)*sqrt(3.)*sqrt(5.)*tmpResult; + } + + return result; + +} + + + complex<double> AbsPsi2STo2K2PiGamLh::chiTo2K_0_Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiTo2K_0, double K_0_Mass0, double K_0_Width0, double K_0_Mass1, double K_0_Width1){ complex<double> result(0.,0.); @@ -1526,7 +1623,7 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToK2KToK2PiAmp(Psi2STo2K2PiGamData::Psi return result; } -complex<double> AbsPsi2STo2K2PiGamLh::chiToK1Tof980_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Dec, double K1Mass, double K1Width, double f980_Mass, double f980_gKK, double f980_gPiPi){ +complex<double> AbsPsi2STo2K2PiGamLh::chiToKjTof980_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjProd, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjDec, Spin Kj, double KjMass, double KjWidth, double f980_Mass, double f980_gKK, double f980_gPiPi){ const double massPi0 = 0.1349766; std::pair <const double, const double> decPair1=make_pair(massPi0, massPi0); @@ -1547,13 +1644,13 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToK1Tof980_piKAmp(Psi2STo2K2PiGamData:: Vector4<double> PiPi(theData->PiPi_HeliChic0_4V.E(), theData->PiPi_HeliChic0_4V.Px(), theData->PiPi_HeliChic0_4V.Py(), theData->PiPi_HeliChic0_4V.Pz()); - Spin lamK1=0; + Spin lamKj=0; complex<double> result(0.,0.); std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator itProd; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator itDec; - for ( itProd=K1Prod.begin(); itProd!=K1Prod.end(); ++itProd){ + for ( itProd=KjProd.begin(); itProd!=KjProd.end(); ++itProd){ boost::shared_ptr<const JPCLS> theState=itProd->first; double theMagProd=itProd->second.first; @@ -1562,26 +1659,71 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToK1Tof980_piKAmp(Psi2STo2K2PiGamData:: complex<double> currentResultDec(0.,0.); - for ( itDec=K1Dec.begin(); itDec!=K1Dec.end(); ++itDec){ + for ( itDec=KjDec.begin(); itDec!=KjDec.end(); ++itDec){ boost::shared_ptr<const JPCLS> theDecState=itDec->first; double theMagDec=itDec->second.first; double thePhiDec=itDec->second.second; complex<double> expiphiDec(cos(thePhiDec), sin(thePhiDec)); - complex<double> ampKpPiPi=conj(theData->DfK1pTof0Kp[1][0][0])*BreitWigner(KpPiPi, K1Mass, K1Width); - complex<double> ampKmPiPi=conj(theData->DfK1mTof0Km[1][0][0])*BreitWigner(KmPiPi, K1Mass, K1Width); + complex<double> ampKpPiPi=conj(theData->DfKjpTof0Kp[Kj][0][0])*BreitWigner(KpPiPi, KjMass, KjWidth); + complex<double> ampKmPiPi=conj(theData->DfKjmTof0Km[Kj][0][0])*BreitWigner(KmPiPi, KjMass, KjWidth); - currentResultDec+=theMagDec*expiphiDec*sqrt(2.*theState->L+1.)*Clebsch(theState->L, 0., theState->S, 0, theState->J, 0)*Clebsch(1,0, 0, 0, theState->S, 0)*Flatte(PiPi, decPair1, decPair2, f980_Mass, f980_gPiPi, f980_gKK)*(ampKpPiPi+ampKmPiPi); + currentResultDec+=theMagDec*expiphiDec*sqrt(2.*theDecState->L+1.) + *Flatte(PiPi, decPair1, decPair2, f980_Mass, f980_gPiPi, f980_gKK)*(ampKpPiPi+ampKmPiPi); } - result+=theMagProd*expiphiProd*currentResultDec; + result+=theMagProd*expiphiProd*sqrt(2.*theState->L+1.)*Clebsch(theState->L, 0., theState->S, 0, theState->J, 0)*Clebsch(Kj,0, 0, 0, theState->S, 0)*currentResultDec; } return result; } +complex<double> AbsPsi2STo2K2PiGamLh::chiToKjTof0_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjProd, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjDec, Spin Kj, double KjMass, double KjWidth, double f0Mass, double f0Width){ + Vector4<double> KpPiPi(theData->KpPiPi_HeliChic0_4V.E(), theData->KpPiPi_HeliChic0_4V.Px(), + theData->KpPiPi_HeliChic0_4V.Py(), theData->KpPiPi_HeliChic0_4V.Pz()); + Vector4<double> KmPiPi(theData->KmPiPi_HeliChic0_4V.E(), theData->KmPiPi_HeliChic0_4V.Px(), + theData->KmPiPi_HeliChic0_4V.Py(), theData->KmPiPi_HeliChic0_4V.Pz()); + + Vector4<double> PiPi(theData->PiPi_HeliChic0_4V.E(), theData->PiPi_HeliChic0_4V.Px(), + theData->PiPi_HeliChic0_4V.Py(), theData->PiPi_HeliChic0_4V.Pz()); + + + Spin lamKj=0; + complex<double> result(0.,0.); + + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator itProd; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator itDec; + + for ( itProd=KjProd.begin(); itProd!=KjProd.end(); ++itProd){ + + boost::shared_ptr<const JPCLS> theState=itProd->first; + double theMagProd=itProd->second.first; + double thePhiProd=itProd->second.second; + complex<double> expiphiProd(cos(thePhiProd), sin(thePhiProd)); + + complex<double> currentResultDec(0.,0.); + + for ( itDec=KjDec.begin(); itDec!=KjDec.end(); ++itDec){ + + boost::shared_ptr<const JPCLS> theDecState=itDec->first; + double theMagDec=itDec->second.first; + double thePhiDec=itDec->second.second; + complex<double> expiphiDec(cos(thePhiDec), sin(thePhiDec)); + + complex<double> ampKpPiPi=conj(theData->DfKjpTof0Kp[Kj][0][0])*BreitWigner(KpPiPi, KjMass, KjWidth); + complex<double> ampKmPiPi=conj(theData->DfKjmTof0Km[Kj][0][0])*BreitWigner(KmPiPi, KjMass, KjWidth); + + + currentResultDec+=theMagDec*expiphiDec*sqrt(2.*theDecState->L+1.)*BreitWigner(PiPi, f0Mass, f0Width)*(ampKpPiPi+ampKmPiPi); + } + + result+=theMagProd*expiphiProd*sqrt(2.*theState->L+1.)*Clebsch(theState->L, 0., theState->S, 0, theState->J, 0)*Clebsch(Kj,0, 0, 0, theState->S, 0)*currentResultDec; + } + return result; + +} void AbsPsi2STo2K2PiGamLh::print(std::ostream& os) const{ os << "AbsPsi2STo2K2PiGamLh::print\n"; diff --git a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh index 7f8335790e29a49d04363be566d6ac9f2131b341..d88521630a1549c2fc68e018ddc086d7ab78c4c0 100644 --- a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh +++ b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh @@ -77,7 +77,9 @@ protected: virtual complex<double> chiToK0K2Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToK0K2, double K0Mass, double K0Width, double K2Mass, double K2Width); -// virtual complex<double> chiToK1ToK1piAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1ToK1dPi, double K1Mass, double K1Width, double K1dMass, double K1dWidth); + virtual complex<double> chiToK0K1Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToK0K1, double K0Mass, double K0Width, double K1Mass, double K1Width); + + virtual complex<double> chiToK1K2Amp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToK1K2, double K1Mass, double K1Width, double K2Mass, double K2Width); virtual complex<double> chiToK1ToK1piAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Dec, double K1Mass, double K1Width, double K1dMass, double K1dWidth); @@ -133,7 +135,10 @@ protected: virtual complex<double> chiToPi0Pi0ToK0KAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi0Pi0Amp, double Pi_0_Mass, double Pi_0_Width, double K0Mass, double K0Width); - virtual complex<double> chiToK1Tof980_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& K1Dec, double K1Mass, double K1Width, double f980_Mass, double f980_gKK, double f980_gPiPi); + virtual complex<double> chiToKjTof980_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjProd, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjDec, Spin Kj, double KjMass, double KjWidth, double f980_Mass, double f980_gKK, double f980_gPiPi); + + virtual complex<double> chiToKjTof0_piKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjProd, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& KjDec, Spin Kj, double KjMass, double KjWidth, double f0Mass, double f0Width); + private: diff --git a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc index 40b55a97a63754d1aaf30df91212a88d3dbefb3b..6d32421f93f2f3009a980d6c419d441b0990ff37 100644 --- a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc +++ b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc @@ -18,7 +18,10 @@ FitParams2K2PiGam::FitParams2K2PiGam() _jpclsMap[paramEnum2K2PiGam::K_0_1430K_0_1430]=theStates.ChiTo2K_0_States(); _jpclsMap[paramEnum2K2PiGam::K_2_1430K_2_1430]=theStates.ChiTo2K_2_1430States(); _jpclsMap[paramEnum2K2PiGam::K_0_1430K_2_1430]=theStates.ChiToK0K2_States(); + _jpclsMap[paramEnum2K2PiGam::ChiToK_0_1430_K892]=theStates.ChiToKst0Kst1States(); _jpclsMap[paramEnum2K2PiGam::K_1_1410K_1_1410]=theStates.ChiTo2K892States(); + _jpclsMap[paramEnum2K2PiGam::ChiToK_1_1410_K892]=theStates.ChiTo2K892States(); + _jpclsMap[paramEnum2K2PiGam::ChiToK_2_1430_K892]=theStates.ChiToKst1Kst2States(); _jpclsMap[paramEnum2K2PiGam::ChiToK_1_1400K]=theStates.ChiToK1400ToK892piStates(); _jpclsMap[paramEnum2K2PiGam::K_1_1400ToK892Pi]=theStates.K1400ToKst1PiStates(); _jpclsMap[paramEnum2K2PiGam::ChiToK_1_1270_K]=theStates.ChiToK1400ToK892piStates(); @@ -32,6 +35,7 @@ FitParams2K2PiGam::FitParams2K2PiGam() _jpclsMap[paramEnum2K2PiGam::KappaK_0_1950]=theStates.ChiTo2K_0_States(); _jpclsMap[paramEnum2K2PiGam::f980_pif1710_k]=theStates.ChiTof0f0States(); _jpclsMap[paramEnum2K2PiGam::f980_kf1710_pi]=theStates.ChiTof0f0States(); + _jpclsMap[paramEnum2K2PiGam::ChiTof1710f1710]=theStates.ChiTof0f0States(); _jpclsMap[paramEnum2K2PiGam::f980f980]=theStates.ChiTof0f0States(); _jpclsMap[paramEnum2K2PiGam::f980f2200]=theStates.ChiTof0f0States(); _jpclsMap[paramEnum2K2PiGam::ChiTof980f_2_2300]=theStates.ChiTof0f2States(); @@ -60,7 +64,14 @@ FitParams2K2PiGam::FitParams2K2PiGam() _jpclsMap[paramEnum2K2PiGam::ChiToK_0_2400ToKf_0_1710]=theStates.ChiToK0K0States(); _jpclsMap[paramEnum2K2PiGam::ChiToK_1_2400K]=theStates.ChiToK1400ToK892piStates(); _jpclsMap[paramEnum2K2PiGam::K_1_2400Tof980K]=theStates.K1pTof0KStates(); + _jpclsMap[paramEnum2K2PiGam::K_1_2400Tof1710K]=theStates.K1pTof0KStates(); + _jpclsMap[paramEnum2K2PiGam::K_1_2400ToK_0_1430Pi]=theStates.K1ToK0PiStates(); + _jpclsMap[paramEnum2K2PiGam::ChiToK_2_2400K]=theStates.ChiToK0K2_States(); + _jpclsMap[paramEnum2K2PiGam::K_2_2400Tof980K]=theStates.K2Tof0KStates(); + _jpclsMap[paramEnum2K2PiGam::K_2_2400Tof1710K]=theStates.K2Tof0KStates(); _jpclsMap[paramEnum2K2PiGam::K892K_1_1680]=theStates.ChiTo2K892States(); + _jpclsMap[paramEnum2K2PiGam::ChiToK1680K1680]=theStates.ChiTo2K892States(); + _jpclsMap[paramEnum2K2PiGam::ChiToK1680K_0_1430]=theStates.ChiToKst0Kst1States(); _jpclsMap[paramEnum2K2PiGam::K892K_1_2300]=theStates.ChiTo2K892States(); _jpclsMap[paramEnum2K2PiGam::sigmaf980]=theStates.ChiTof0f0States(); _jpclsMap[paramEnum2K2PiGam::ChiToSigmaf1370]=theStates.ChiTof0f0States(); @@ -116,7 +127,10 @@ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collecti else if (index==paramEnum2K2PiGam::K_0_1430K_0_1430) return params.ChiTo2K_0_1430; else if (index==paramEnum2K2PiGam::K_2_1430K_2_1430) return params.ChiTo2K_2_1430; else if (index==paramEnum2K2PiGam::K_0_1430K_2_1430) return params.ChiToK_0_1430_K_2_1430; + else if (index==paramEnum2K2PiGam::ChiToK_0_1430_K892) return params.ChiToK_0_1430_K892; else if (index==paramEnum2K2PiGam::K_1_1410K_1_1410) return params.ChiToK_1_1410_K_1_1410; + else if (index==paramEnum2K2PiGam::ChiToK_1_1410_K892) return params.ChiToK_1_1410_K892; + else if (index==paramEnum2K2PiGam::ChiToK_2_1430_K892) return params.ChiToK_2_1430_K892; else if (index==paramEnum2K2PiGam::ChiToK_1_1400K) return params.ChiToK_1_1400K; else if (index==paramEnum2K2PiGam::K_1_1400ToK892Pi) return params.K1400ToK892Pi; else if (index==paramEnum2K2PiGam::ChiToK_1_1270_K) return params.ChiToK_1_1270_K; @@ -130,6 +144,7 @@ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collecti else if (index==paramEnum2K2PiGam::KappaK_0_1950) return params.ChiToKappaK_0_1950; else if (index==paramEnum2K2PiGam::f980_pif1710_k) return params.f980_pif1710_k; else if (index==paramEnum2K2PiGam::f980_kf1710_pi) return params.f980_kf1710_pi; + else if (index==paramEnum2K2PiGam::ChiTof1710f1710) return params.ChiTof1710f1710; else if (index==paramEnum2K2PiGam::f980f980) return params.ChiTof980f980; else if (index==paramEnum2K2PiGam::f980f2200) return params.ChiTof980f2200; else if (index==paramEnum2K2PiGam::ChiTof980f_2_2300) return params.ChiTof980f_2_2300; @@ -158,7 +173,14 @@ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collecti else if (index==paramEnum2K2PiGam::ChiToK_0_2400ToKf_0_1710) return params.ChiToK_0_2400ToKf_0_1710; else if (index==paramEnum2K2PiGam::ChiToK_1_2400K) return params.ChiToK_1_2400K; else if (index==paramEnum2K2PiGam::K_1_2400Tof980K) return params.K_1_2400Tof980K; + else if (index==paramEnum2K2PiGam::K_1_2400Tof1710K) return params.K_1_2400Tof1710K; + else if (index==paramEnum2K2PiGam::K_1_2400ToK_0_1430Pi) return params.K_1_2400ToK_0_1430Pi; + else if (index==paramEnum2K2PiGam::ChiToK_2_2400K) return params.ChiToK_2_2400K; + else if (index==paramEnum2K2PiGam::K_2_2400Tof980K) return params.K_2_2400Tof980K; + else if (index==paramEnum2K2PiGam::K_2_2400Tof1710K) return params.K_2_2400Tof1710K; else if (index==paramEnum2K2PiGam::K892K_1_1680) return params.ChiToK892K1680; + else if (index==paramEnum2K2PiGam::ChiToK1680K1680) return params.ChiToK1680K1680; + else if (index==paramEnum2K2PiGam::ChiToK1680K_0_1430) return params.ChiToK1680K_0_1430; else if (index==paramEnum2K2PiGam::K892K_1_2300) return params.ChiToK892K2300; else if (index==paramEnum2K2PiGam::sigmaf980) return params.ChiToSigmaf980; else if (index==paramEnum2K2PiGam::ChiToSigmaf1370) return params.ChiToSigmaf1370; @@ -221,6 +243,7 @@ pair<double, double>& FitParams2K2PiGam::massPair(param2K2PiGam& params, unsigne else if (index==paramEnum2K2PiGam::f_2_2300) return params.Bwf_2_2300; else if (index==paramEnum2K2PiGam::K_0_2400) return params.BwK_0_2400; else if (index==paramEnum2K2PiGam::K_1_2400) return params.BwK_1_2400; + else if (index==paramEnum2K2PiGam::K_2_2400) return params.BwK_2_2400; else if (index==paramEnum2K2PiGam::K_0_1950) return params.BwK_0_1950; else if (index==paramEnum2K2PiGam::K_1_1680) return params.BwK_1_1680; else if (index==paramEnum2K2PiGam::K_1_2300) return params.BwK_1_2300; @@ -346,11 +369,11 @@ void FitParams2K2PiGam::setMnUsrParamsDec(MnUserParameters& upar, param2K2PiGam& double magErr=errPair.first; double phiErr=errPair.second; - double magMin=magVal-magErr; + double magMin=magVal-2.*magErr; if (magMin<0.) magMin=0.; - upar.Add(magStr, magVal, magErr, magMin, magVal+magErr); - upar.Add(phiStr, phiVal, phiErr, -3*M_PI, 3*M_PI); + upar.Add(magStr, magVal, magErr, magMin, magVal+2.*magErr); + upar.Add(phiStr, phiVal, phiErr, -3.*M_PI, 3.*M_PI); counter++; } diff --git a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh index b4f0566b45a79031d3c73715d0524609147db5a2..0ebc98d28fd35567c8c907cd7cd406f75611fc52 100644 --- a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh +++ b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh @@ -25,56 +25,67 @@ using namespace ROOT::Minuit2; struct param2K2PiGam { - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > PsiToChiGam; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K892; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2Kappa; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K_2_1430; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K_0_1430; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430_K_2_1430; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K_1_1410; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1400K; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K1400ToK892Pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1270_K; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_1270ToK892Pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_1270ToK_0_1430Pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_2400K; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof980K; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1710_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1710_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f980; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f2200; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f_2_2300; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > PsiToChiGam; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K892; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2Kappa; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K_2_1430; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTo2K_0_1430; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430_K_2_1430; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430_K892; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K_1_1410; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K892; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_2_1430_K892; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1400K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K1400ToK892Pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1270_K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_1270ToK892Pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_1270ToK_0_1430Pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_2400K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof980K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof1710K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400ToK_0_1430Pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_2_2400K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_2_2400Tof980K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_2_2400Tof1710K; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1710_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1710_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof1710f1710; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f980; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f2200; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof980f_2_2300; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof_2_2300sigma; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf980; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf1370; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf1710; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf2200; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1370_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1370_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1500_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1500_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif1370_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf1370_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1270_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1270_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1430_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1430_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1525_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1525_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1950_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1950_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1500_pif_2_1525_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1500_kf_2_1525_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif_2_1430_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf_2_1430_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif_2_1950_k; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf_2_1950_pi; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_2400ToKf980; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_2400ToKf_0_1710; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430K_0_1950; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf980; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf1370; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf1710; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToSigmaf2200; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1370_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1370_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif1500_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf1500_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif1370_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf1370_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1270_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1270_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1430_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1430_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1525_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1525_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_pif_2_1950_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f980_kf_2_1950_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1500_pif_2_1525_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1500_kf_2_1525_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif_2_1430_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf_2_1430_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_pif_2_1950_k; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > f1710_kf_2_1950_pi; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_2400ToKf980; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_2400ToKf_0_1710; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430K_0_1950; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToKappaK_0_1430; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToKappaK_0_1950; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK892K1680; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK1680K1680; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK1680K_0_1430; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK892K2300; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_0_1460ToK892Pi; std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_0_1460ToK_0_1430Pi; @@ -122,6 +133,7 @@ struct param2K2PiGam pair<double, double> Bwf_2_2300; pair<double, double> BwK_0_2400; pair<double, double> BwK_1_2400; + pair<double, double> BwK_2_2400; pair<double, double> BwK_2_1950; pair<double, double> BwK_0_1950; pair<double, double> BwK_1_1680; @@ -143,13 +155,15 @@ struct param2K2PiGam }; struct paramEnum2K2PiGam{ - enum {ChiGam=0, K892K892, KappaKappa, K_0_1430K_0_1430, K_2_1430K_2_1430, K_0_1430K_2_1430, K_1_1410K_1_1410, + enum {ChiGam=0, K892K892, KappaKappa, K_0_1430K_0_1430, ChiToK_0_1430_K892, + K_2_1430K_2_1430, K_0_1430K_2_1430, K_1_1410K_1_1410, ChiToK_1_1410_K892, ChiToK_2_1430_K892, ChiToK_1_1400K, K_1_1400ToK892Pi, ChiToK_1_1270_K, K_1_1270ToK892Pi, K_1_1270ToK_0_1430Pi, ChiToK_1_1650K, K_1_1650ToK892Pi, K_1_1650ToK_0_1430Pi, K_0_1430K_0_1950, KappaK_0_1950, ChiToKappaK_0_1430, - ChiToK_1_2400K, K_1_2400Tof980K, - f980_pif1710_k, f980_kf1710_pi, f980f980, f980f2200, ChiTof980f_2_2300, ChiTof_2_2300sigma, + ChiToK_1_2400K, K_1_2400Tof980K, K_1_2400Tof1710K, K_1_2400ToK_0_1430Pi, + ChiToK_2_2400K, K_2_2400Tof980K, K_2_2400Tof1710K, + f980_pif1710_k, f980_kf1710_pi, ChiTof1710f1710, f980f980, f980f2200, ChiTof980f_2_2300, ChiTof_2_2300sigma, f980_pif1370_k, f980_kf1370_pi, f980_pif1500_k, f980_kf1500_pi, f1710_pif1370_k, f1710_kf1370_pi, - K_0_2400KToKf980, ChiToK_0_2400ToKf_0_1710, K892K_1_1680, K892K_1_2300, + K_0_2400KToKf980, ChiToK_0_2400ToKf_0_1710, K892K_1_1680, ChiToK1680K1680, ChiToK1680K_0_1430, K892K_1_2300, sigmaf980, ChiToSigmaf1370, sigmaf1710, sigmaf2200, K_0_1460ToK892Pi, K_0_1460ToK_0_1430Pi, K_0_1830ToK892Pi, K_0_1830ToK_0_1430Pi, f980_pif_2_1270_k, f980_kf_2_1270_pi, f980_pif_2_1430_k, f980_kf_2_1430_pi, f980_pif_2_1525_k, f980_kf_2_1525_pi, f980_pif_2_1950_k, f980_kf_2_1950_pi, f1500_pif_2_1525_k, f1500_kf_2_1525_pi, f1710_pif_2_1430_k, f1710_kf_2_1430_pi, f1710_pif_2_1950_k, f1710_kf_2_1950_pi, @@ -159,7 +173,7 @@ struct paramEnum2K2PiGam{ ChiToK_2_1770K, K_2_1770ToK_2_1430Pi, ChiToK_0_1430KPi, ChiToK892KPi, nAmps, K892=nAmps,Kappa, K_0_1430, K_1_1400, K_1_1410, K_2_1430, K_1_1270,K_1_1650, - f1500, f1710, f2200, sigma, f1360, f1370, K_0_2400, K_1_2400, K_0_1950, K_1_1680, K_1_2300, K_0_1460, K_0_1830, + f1500, f1710, f2200, sigma, f1360, f1370, K_0_2400, K_1_2400, K_2_2400, K_0_1950, K_1_1680, K_1_2300, K_0_1460, K_0_1830, f_2_1270, f_2_1430, f_2_1525, f_2_1950, f_2_2300, Pi_2_1670, Pi1800, Pi3000, Pi_2_2285, K_2_1770, nMasses, @@ -170,13 +184,15 @@ struct paramEnum2K2PiGam{ static const std::string& name(unsigned int t) { static std::string fitName[paramEnum2K2PiGam::nPhaseSpace] - ={"ChiGam", "K892K892", "KappaKappa", "K_0_1430K_0_1430", "K_2_1430K_2_1430", "K_0_1430K_2_1430", "K_1_1410K_1_1410", + ={"ChiGam", "K892K892", "KappaKappa", "K_0_1430K_0_1430", "ChiToK_0_1430_K892", + "K_2_1430K_2_1430", "K_0_1430K_2_1430", "K_1_1410K_1_1410","ChiToK_1_1410_K892","ChiToK_2_1430_K892", "ChiToK_1_1400K", "K_1_1400ToK892Pi", "ChiToK_1_1270_K", "K_1_1270ToK892Pi", "K_1_1270ToK_0_1430Pi", "ChiToK_1_1650K", "K_1_1650ToK892Pi", "K_1_1650ToK_0_1430Pi","K_0_1430K_0_1950", "KappaK_0_1950","ChiToKappaK_0_1430", - "ChiToK_1_2400K", "K_1_2400Tof980K", - "f980_pif1710_k", "f980_kf1710_pi", "f980f980", "f980f2200", "ChiTof980f_2_2300","ChiTof_2_2300sigma", + "ChiToK_1_2400K", "K_1_2400Tof980K","K_1_2400Tof1710K", "K_1_2400ToK_0_1430Pi", + "ChiToK_2_2400K", "K_2_2400Tof980K","K_2_2400Tof1710K", + "f980_pif1710_k", "f980_kf1710_pi", "ChiTof1710f1710", "f980f980", "f980f2200", "ChiTof980f_2_2300","ChiTof_2_2300sigma", "f980_pif1370_k", "f980_kf1370_pi", "f980_pif1500_k", "f980_kf1500_pi", "f1710_pif1370_k", "f1710_kf1370_pi", - "K_0_2400KToKf980", "ChiToK_0_2400ToKf_0_1710", "K892K_1_1680", "K892K_1_2300", + "K_0_2400KToKf980", "ChiToK_0_2400ToKf_0_1710", "K892K_1_1680", "ChiToK1680K1680", "ChiToK1680K_0_1430", "K892K_1_2300", "sigmaf980", "ChiToSigmaf1370", "sigmaf1710", "sigmaf2200", "K_0_1460ToK892Pi", "K_0_1460ToK_0_1430Pi","K_0_1830ToK892Pi", "K_0_1830ToK_0_1430Pi", "f980_pif_2_1270_k", "f980_kf_2_1270_pi","f980_pif_2_1430_k", "f980_kf_2_1430_pi","f980_pif_2_1525_k", "f980_kf_2_1525_pi", "f980_pif_2_1950_k", "f980_kf_2_1950_pi", "f1500_pif_2_1525_k", "f1500_kf_2_1525_pi", "f1710_pif_2_1430_k", "f1710_kf_2_1430_pi", "f1710_pif_2_1950_k", "f1710_kf_2_1950_pi", @@ -184,9 +200,8 @@ struct paramEnum2K2PiGam{ "ChiToPi1800Pi0Tof980","ChiToPi1800Pi0Tof1370","ChiToPi1800Pi0ToKappa","ChiToPi1800Pi0ToK892K","ChiToPi3000Pi0ToK892K","ChiToPi3000Pi0ToK_0_1950K", "ChiToPi_2_2285Pi","Pi_2_2285Tof1700Pi","Pi_2_2285ToK892KPi","Pi_2_2285ToK_0_1430K","Pi_2_2285ToK_2_1430K", "ChiToK_2_1770K","K_2_1770ToK_2_1430Pi","ChiToK_0_1430KPi","ChiToK892KPi", - "K892", "Kappa", "K_0_1430", "K_1_1400", "K_1_1410", "K_2_1430", "K_1_1270", "K_1_1650", - "f1500", "f1710", "f2200", "sigma", "f1360", "f1370", "K_0_2400", "K_1_2400", "K_0_1950", "K_1_1680", "K_1_2300", "K_0_1460", "K_0_1830", + "f1500", "f1710", "f2200", "sigma", "f1360", "f1370", "K_0_2400", "K_1_2400", "K_2_2400", "K_0_1950", "K_1_1680", "K_1_2300", "K_0_1460", "K_0_1830", "f_2_1270", "f_2_1430", "f_2_1525","f_2_1950","f_2_2300", "Pi_2_1670","Pi1800","Pi3000","Pi_2_2285","K_2_1770", diff --git a/Examples/Psi2STo2K2PiGam/Hyp1Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp1Lh.cc index 4b0205e697d4040b53ce371259ed111ab7ae422a..33f2e8fd9ab5f28a7e346b1ff66037b908efcf2c 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp1Lh.cc +++ b/Examples/Psi2STo2K2PiGam/Hyp1Lh.cc @@ -12,7 +12,11 @@ Hyp1Lh::Hyp1Lh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList> theEvtList, const ,_K0_1430_K0_1430Hyp(true) ,_K2_1430_K2_1430Hyp(false) ,_K0_1430_K2_1430Hyp(false) + ,_K0_1430_K892Hyp1(true) + ,_K2_1430_K892Hyp1(true) ,_K1_1410_K1_1410Hyp(false) + ,_K1_1410_K892Hyp1(true) + ,_f1710_f1710Hyp1(true) ,_nFitParams(0) { setUp(hypMap); @@ -25,7 +29,11 @@ Hyp1Lh::Hyp1Lh( boost::shared_ptr<AbsPsi2STo2K2PiGamLh> theLhPtr, const std::map ,_K0_1430_K0_1430Hyp(true) ,_K2_1430_K2_1430Hyp(false) ,_K0_1430_K2_1430Hyp(false) + ,_K0_1430_K892Hyp1(true) + ,_K2_1430_K892Hyp1(true) ,_K1_1410_K1_1410Hyp(false) + ,_K1_1410_K892Hyp1(true) + ,_f1710_f1710Hyp1(true) ,_nFitParams(0) { setUp(hypMap); @@ -79,11 +87,29 @@ complex<double> Hyp1Lh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2 if (_K0_1430_K0_1430Hyp) result+=chiTo2K_0_Amp(theData, ChiTo2K_0_1430, K_0_1430Mass, K_0_1430Width, K_0_1430Mass, K_0_1430Width); if (_K2_1430_K2_1430Hyp) result+=chiTo2K_2_Amp(theData, ChiTo2K_2_1430, K_2_1430Mass, K_2_1430Width); if (_K0_1430_K2_1430Hyp) result+=chiToK0K2Amp(theData, ChiToK_0_1430_K_2_1430, K_0_1430Mass, K_0_1430Width, K_2_1430Mass, K_2_1430Width); - if (_K1_1410_K1_1410Hyp){ + if(_K0_1430_K892Hyp1){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430_K892=theParamVal.ChiToK_0_1430_K892; + result+=chiToK0K1Amp(theData, ChiToK_0_1430_K892, K_0_1430Mass, K_0_1430Width, K892Mass, K892Width); + } + + if(_K2_1430_K892Hyp1){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_2_1430_K892=theParamVal.ChiToK_2_1430_K892; + result+=chiToK1K2Amp(theData, ChiToK_2_1430_K892, K892Mass, K892Width, K_2_1430Mass, K_2_1430Width); + } + + if (_K1_1410_K1_1410Hyp || _K1_1410_K892Hyp1){ double K_1_1410Mass=theParamVal.BwK_1_1410.first; double K_1_1410Width=theParamVal.BwK_1_1410.second; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K_1_1410=theParamVal.ChiToK_1_1410_K_1_1410; - result+=chiToK1K1Amp(theData, ChiToK_1_1410_K_1_1410, K_1_1410Mass, K_1_1410Width, K_1_1410Mass, K_1_1410Width); + + if(_K1_1410_K1_1410Hyp){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K_1_1410=theParamVal.ChiToK_1_1410_K_1_1410; + result+=chiToK1K1Amp(theData, ChiToK_1_1410_K_1_1410, K_1_1410Mass, K_1_1410Width, K_1_1410Mass, K_1_1410Width); + } + + if(_K1_1410_K892Hyp1){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_1410_K892=theParamVal.ChiToK_1_1410_K892; + result+=chiToK1K1Amp(theData, ChiToK_1_1410_K892, K892Mass, K892Width, K_1_1410Mass, K_1_1410Width); + } } //Chi_c0 decay to K1*(1400) -> K1*(892) pi0 -> (K pi0) pi0 @@ -113,6 +139,11 @@ complex<double> Hyp1Lh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2 //Chi_c0 decay to f0(980) f0(980) -> (pi0 pi0) (K K) result+=chiTof980f980Amp(theData, ChiTof980f980, f980_Mass, f980_gPiPi, f980_gKK); + if(_f1710_f1710Hyp1){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiTof1710f1710=theParamVal.ChiTof1710f1710; + result+=chiTof0_pif0_kAmp(theData, ChiTof1710f1710, f1710Mass, f1710Width, f1710Mass, f1710Width); + } + return result; } @@ -304,6 +335,24 @@ void Hyp1Lh::setUp(const std::map<const std::string, bool>& hypMap){ _hypMap[iter->first]= iter->second; Info<< "hypothesis " << iter->first << "\t" << _K0_1430_K2_1430Hyp <<endmsg; } + else Alert << "hypothesis K0_1430_K2_1430Hyp not set!!!" <<endmsg; + + + iter= hypMap.find("K0_1430_K892Hyp1"); + if (iter !=hypMap.end()){ + _K0_1430_K892Hyp1= iter->second; + _hypMap[iter->first]= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K0_1430_K892Hyp1 <<endmsg; + } + else Alert << "hypothesis K0_1430_K892Hyp1 not set!!!" <<endmsg; + + iter= hypMap.find("K2_1430_K892Hyp1"); + if (iter !=hypMap.end()){ + _K2_1430_K892Hyp1= iter->second; + _hypMap[iter->first]= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K2_1430_K892Hyp1 <<endmsg; + } + else Alert << "hypothesis K2_1430_K892Hyp1 not set!!!" <<endmsg; iter= hypMap.find("K1_1410_K1_1410Hyp"); if (iter !=hypMap.end()){ @@ -311,9 +360,25 @@ void Hyp1Lh::setUp(const std::map<const std::string, bool>& hypMap){ _hypMap[iter->first]= iter->second; Info<< "hypothesis " << iter->first << "\t" << _K1_1410_K1_1410Hyp <<endmsg; } - else Alert << "hypothesis K1_1410_K1_1410Hyp not set!!!" <<endmsg; + iter= hypMap.find("K1_1410_K892Hyp1"); + if (iter !=hypMap.end()){ + _K1_1410_K892Hyp1= iter->second; + _hypMap[iter->first]= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K1_1410_K892Hyp1 <<endmsg; + } + else Alert << "hypothesis K1_1410_K892Hyp1 not set!!!" <<endmsg; + + iter= hypMap.find("f1710_f1710Hyp1"); + if (iter !=hypMap.end()){ + _f1710_f1710Hyp1= iter->second; + _hypMap[iter->first]= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _f1710_f1710Hyp1 <<endmsg; + } + else Alert << "hypothesis f1710_f1710Hyp1 not set!!!" <<endmsg; + + _ampVec.push_back(paramEnum2K2PiGam::ChiGam); _ampVec.push_back(paramEnum2K2PiGam::K892K892); _ampVec.push_back(paramEnum2K2PiGam::f980_pif1710_k); @@ -332,17 +397,31 @@ void Hyp1Lh::setUp(const std::map<const std::string, bool>& hypMap){ if(_K2_1430_K2_1430Hyp) _ampVec.push_back(paramEnum2K2PiGam::K_2_1430K_2_1430); if(_K0_1430_K0_1430Hyp) _ampVec.push_back(paramEnum2K2PiGam::K_0_1430K_0_1430); if(_K0_1430_K2_1430Hyp) _ampVec.push_back(paramEnum2K2PiGam::K_0_1430K_2_1430); - if(_K1_1410_K1_1410Hyp){ - _ampVec.push_back(paramEnum2K2PiGam::K_1_1410K_1_1410); + if(_K0_1430_K892Hyp1) _ampVec.push_back(paramEnum2K2PiGam::ChiToK_0_1430_K892); + if(_K2_1430_K892Hyp1) _ampVec.push_back(paramEnum2K2PiGam::ChiToK_2_1430_K892); + + if(_K1_1410_K1_1410Hyp || _K1_1410_K892Hyp1){ _massVec.push_back(paramEnum2K2PiGam::K_1_1410); + + if(_K1_1410_K1_1410Hyp){ + _ampVec.push_back(paramEnum2K2PiGam::K_1_1410K_1_1410); + } + + if(_K1_1410_K892Hyp1){ + _ampVec.push_back(paramEnum2K2PiGam::ChiToK_1_1410_K892); + } } + if(_f1710_f1710Hyp1) _ampVec.push_back(paramEnum2K2PiGam::ChiTof1710f1710); + _massVec.push_back(paramEnum2K2PiGam::K892); _massVec.push_back(paramEnum2K2PiGam::f1710); if(_K1_1270Hyp) _massVec.push_back(paramEnum2K2PiGam::K_1_1270); if(_K1_1400Hyp) _massVec.push_back(paramEnum2K2PiGam::K_1_1400); - if(_K2_1430_K2_1430Hyp || _K0_1430_K2_1430Hyp) _massVec.push_back(paramEnum2K2PiGam::K_2_1430); - if(_K0_1430_K0_1430Hyp || _K0_1430_K2_1430Hyp || _K1_1270Hyp) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); + if(_K2_1430_K2_1430Hyp || _K0_1430_K2_1430Hyp || _K2_1430_K892Hyp1) _massVec.push_back(paramEnum2K2PiGam::K_2_1430); + if(_K0_1430_K0_1430Hyp || _K0_1430_K2_1430Hyp || _K1_1270Hyp || _K0_1430_K892Hyp1) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); + + std::vector<unsigned int>::iterator ampIt; for (ampIt=_ampVec.begin(); ampIt!=_ampVec.end(); ++ampIt){ diff --git a/Examples/Psi2STo2K2PiGam/Hyp1Lh.hh b/Examples/Psi2STo2K2PiGam/Hyp1Lh.hh index 792ce88a81641b6331155e6a4663de4c698ac851..616c2a4f8b8e44a00d76c7fc3095b11254387621 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp1Lh.hh +++ b/Examples/Psi2STo2K2PiGam/Hyp1Lh.hh @@ -58,7 +58,11 @@ protected: bool _K0_1430_K0_1430Hyp; bool _K2_1430_K2_1430Hyp; bool _K0_1430_K2_1430Hyp; + bool _K0_1430_K892Hyp1; + bool _K2_1430_K892Hyp1; bool _K1_1410_K1_1410Hyp; + bool _K1_1410_K892Hyp1; + bool _f1710_f1710Hyp1; std::map<const std::string, bool> _hypMap; virtual complex<double> chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData); diff --git a/Examples/Psi2STo2K2PiGam/Hyp5Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp5Lh.cc index 9187426673f3c2f6649806ad7b48d43418992016..c23ca579134b6fab7f2d5e5ae2a0dc9b98e4d66f 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp5Lh.cc +++ b/Examples/Psi2STo2K2PiGam/Hyp5Lh.cc @@ -13,6 +13,10 @@ Hyp5Lh::Hyp5Lh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList> theEvtList, const , _K_0_2400KHyp5(true) ,_K_0_2400KTof_0_1710Hyp5(true) ,_K_1_2400KHyp5(true) + ,_K_1_2400KTof_0_1710Hyp5(true) + ,_K_1_2400KToK_0_1430Hyp5(true) + ,_K_2_2400KTof980Hyp5(true) + ,_K_2_2400KTof_0_1710Hyp5(true) ,_ChiToK_0_1430KPiHyp5(true) ,_ChiToK892KPiHyp5(true) , _nFitParams(0) @@ -26,6 +30,10 @@ Hyp5Lh::Hyp5Lh( boost::shared_ptr<AbsPsi2STo2K2PiGamLh> theLhPtr, const std::map , _K_0_2400KHyp5(true) ,_K_0_2400KTof_0_1710Hyp5(true) ,_K_1_2400KHyp5(true) + ,_K_1_2400KTof_0_1710Hyp5(true) + ,_K_1_2400KToK_0_1430Hyp5(true) + ,_K_2_2400KTof980Hyp5(true) + ,_K_2_2400KTof_0_1710Hyp5(true) ,_ChiToK_0_1430KPiHyp5(true) ,_ChiToK892KPiHyp5(true) , _nFitParams(0) @@ -68,12 +76,48 @@ complex<double> Hyp5Lh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2 } - if (_K_1_2400KHyp5){ + if (_K_1_2400KHyp5 || _K_1_2400KTof_0_1710Hyp5 || _K_1_2400KToK_0_1430Hyp5){ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_1_2400K=theParamVal.ChiToK_1_2400K; - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof980K=theParamVal.K_1_2400Tof980K; double K_1_2400Mass=theParamVal.BwK_1_2400.first; double K_1_2400Width=theParamVal.BwK_1_2400.second; - result+=chiToK1Tof980_piKAmp(theData, ChiToK_1_2400K, K_1_2400Tof980K, K_1_2400Mass, K_1_2400Width, f980_Mass, f980_gKK, f980_gPiPi); + + if(_K_1_2400KHyp5){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof980K=theParamVal.K_1_2400Tof980K; + result+=chiToKjTof980_piKAmp(theData, ChiToK_1_2400K, K_1_2400Tof980K, 1, K_1_2400Mass, K_1_2400Width, f980_Mass, f980_gKK, f980_gPiPi); + } + + if(_K_1_2400KTof_0_1710Hyp5){ + double f1710Mass=theParamVal.Bwf1710.first; + double f1710Width=theParamVal.Bwf1710.second; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400Tof1710K=theParamVal.K_1_2400Tof1710K; + result+=chiToKjTof0_piKAmp(theData, ChiToK_1_2400K, K_1_2400Tof1710K, 1, K_1_2400Mass, K_1_2400Width, f1710Mass, f1710Width); + } + + if(_K_1_2400KToK_0_1430Hyp5){ + double K_0_1430Mass=theParamVal.BwK_0_1430.first; + double K_0_1430Width=theParamVal.BwK_0_1430.second; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_1_2400ToK_0_1430Pi=theParamVal.K_1_2400ToK_0_1430Pi; + result+=chiToK1ToK0piAmp(theData, ChiToK_1_2400K, K_1_2400ToK_0_1430Pi, K_1_2400Mass, K_1_2400Width, K_0_1430Mass, K_0_1430Width); + } + } + + + if(_K_2_2400KTof980Hyp5 || _K_2_2400KTof_0_1710Hyp5){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_2_2400K=theParamVal.ChiToK_2_2400K; + double K_2_2400Mass=theParamVal.BwK_2_2400.first; + double K_2_2400Width=theParamVal.BwK_2_2400.second; + + if(_K_2_2400KTof980Hyp5){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_2_2400Tof980K=theParamVal.K_2_2400Tof980K; + result+=chiToKjTof980_piKAmp(theData, ChiToK_2_2400K, K_2_2400Tof980K, 2, K_2_2400Mass, K_2_2400Width, f980_Mass, f980_gKK, f980_gPiPi); + } + + if (_K_2_2400KTof_0_1710Hyp5){ + double f1710Mass=theParamVal.Bwf1710.first; + double f1710Width=theParamVal.Bwf1710.second; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_2_2400Tof1710K=theParamVal.K_2_2400Tof1710K; + result+=chiToKjTof0_piKAmp(theData, ChiToK_2_2400K, K_2_2400Tof1710K, 2, K_2_2400Mass, K_2_2400Width, f1710Mass, f1710Width); + } } if(_ChiToK_0_1430KPiHyp5){ @@ -251,6 +295,43 @@ void Hyp5Lh::setUp(const std::map<const std::string, bool>& hypMap){ } else Alert << "hypothesis K_1_2400KHyp5 not set!!!" <<endmsg; + + iter= hypMap.find("K_1_2400KTof_0_1710Hyp5"); + if (iter !=hypMap.end()){ + _K_1_2400KTof_0_1710Hyp5= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K_1_2400KTof_0_1710Hyp5 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K_1_2400KTof_0_1710Hyp5 not set!!!" <<endmsg; + + iter= hypMap.find("K_1_2400KToK_0_1430Hyp5"); + if (iter !=hypMap.end()){ + _K_1_2400KToK_0_1430Hyp5= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K_1_2400KToK_0_1430Hyp5 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K_1_2400KToK_0_1430Hyp5 not set!!!" <<endmsg; + + + + iter= hypMap.find("K_2_2400KTof980Hyp5"); + if (iter !=hypMap.end()){ + _K_2_2400KTof980Hyp5= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K_2_2400KTof980Hyp5 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K_2_2400KTof980Hyp5 not set!!!" <<endmsg; + + + iter= hypMap.find("K_2_2400KTof_0_1710Hyp5"); + if (iter !=hypMap.end()){ + _K_2_2400KTof_0_1710Hyp5= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K_2_2400KTof_0_1710Hyp5 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K_2_2400KTof_0_1710Hyp5 not set!!!" <<endmsg; + + iter= hypMap.find("K_0_2400KTof_0_1710Hyp5"); if (iter !=hypMap.end()){ @@ -291,12 +372,34 @@ void Hyp5Lh::setUp(const std::map<const std::string, bool>& hypMap){ } - if(_K_1_2400KHyp5){ + if(_K_1_2400KHyp5 || _K_1_2400KTof_0_1710Hyp5 || _K_1_2400KToK_0_1430Hyp5){ _ampVec.push_back(paramEnum2K2PiGam::ChiToK_1_2400K); - _ampVec.push_back(paramEnum2K2PiGam::K_1_2400Tof980K); _massVec.push_back(paramEnum2K2PiGam::K_1_2400); + if(_K_1_2400KHyp5){ + _ampVec.push_back(paramEnum2K2PiGam::K_1_2400Tof980K); + } + if(_K_1_2400KTof_0_1710Hyp5){ + _ampVec.push_back(paramEnum2K2PiGam::K_1_2400Tof1710K); + } + if(_K_1_2400KToK_0_1430Hyp5){ + _ampVec.push_back(paramEnum2K2PiGam::K_1_2400ToK_0_1430Pi); + } + } + + if(_K_2_2400KTof980Hyp5 || _K_2_2400KTof_0_1710Hyp5){ + _ampVec.push_back(paramEnum2K2PiGam::ChiToK_2_2400K); + _massVec.push_back(paramEnum2K2PiGam::K_2_2400); + + if(_K_2_2400KTof980Hyp5){ + _ampVec.push_back(paramEnum2K2PiGam::K_2_2400Tof980K); + } + + if(_K_2_2400KTof_0_1710Hyp5){ + _ampVec.push_back(paramEnum2K2PiGam::K_2_2400Tof1710K); + } } + if(_ChiToK_0_1430KPiHyp5){ _ampVec.push_back(paramEnum2K2PiGam::ChiToK_0_1430KPi); } diff --git a/Examples/Psi2STo2K2PiGam/Hyp5Lh.hh b/Examples/Psi2STo2K2PiGam/Hyp5Lh.hh index 8229526b0096f944533e2b55591704db0ca265cc..5507ea0d0a0afb3656df1b2bbd0c1a1a2d4de1f1 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp5Lh.hh +++ b/Examples/Psi2STo2K2PiGam/Hyp5Lh.hh @@ -55,6 +55,10 @@ protected: bool _K_0_2400KHyp5; bool _K_0_2400KTof_0_1710Hyp5; bool _K_1_2400KHyp5; + bool _K_1_2400KTof_0_1710Hyp5; + bool _K_1_2400KToK_0_1430Hyp5; + bool _K_2_2400KTof980Hyp5; + bool _K_2_2400KTof_0_1710Hyp5; bool _ChiToK_0_1430KPiHyp5; bool _ChiToK892KPiHyp5; diff --git a/Examples/Psi2STo2K2PiGam/Hyp6Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp6Lh.cc index ddcb9a627a31c14eb44624da145097fda0491239..2c78cb8a4a7b42fe4b786580ba37d9f564d17ad2 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp6Lh.cc +++ b/Examples/Psi2STo2K2PiGam/Hyp6Lh.cc @@ -243,7 +243,7 @@ void Hyp6Lh::setUp(const std::map<const std::string, bool>& hypMap){ } if (_K_0_1430K_0_1950Hyp6 || _KappaK_0_1430Hyp6){ - if(!_K0_1430_K0_1430Hyp && !_K0_1430_K2_1430Hyp && !_K1_1270Hyp) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); + if(!_K0_1430_K0_1430Hyp && !_K0_1430_K2_1430Hyp && !_K1_1270Hyp && !_K0_1430_K892Hyp1) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); } diff --git a/Examples/Psi2STo2K2PiGam/Hyp7Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp7Lh.cc index ee8686ddbe94a2c0ab95145ea0b060a4ea80a2da..b8f993dfc4b986a645071598c5926745108c7881 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp7Lh.cc +++ b/Examples/Psi2STo2K2PiGam/Hyp7Lh.cc @@ -11,6 +11,8 @@ Hyp7Lh::Hyp7Lh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList> theEvtList, const Hyp6Lh(theEvtList, hypMap ) ,_KappaHyp(true) ,_K1_1680Hyp(true) + ,_K1_1680K1_1680Hyp7(true) + ,_K1_1680K0_1430Hyp7(true) ,_K1_2300Hyp(true) ,_nFitParams(0) { @@ -21,6 +23,8 @@ Hyp7Lh::Hyp7Lh( boost::shared_ptr<AbsPsi2STo2K2PiGamLh> theLhPtr, const std::map Hyp6Lh(theLhPtr->getEventList(), hypMap) ,_KappaHyp(true) ,_K1_1680Hyp(true) + ,_K1_1680K1_1680Hyp7(true) + ,_K1_1680K0_1430Hyp7(true) ,_K1_2300Hyp(true) ,_nFitParams(0) { @@ -47,11 +51,29 @@ complex<double> Hyp7Lh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2 result+=chiTo2K_0_Amp(theData, ChiTo2Kappa, KappaMass, KappaWidth, KappaMass, KappaWidth); } - if(_K1_1680Hyp){ - std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK892K1680=theParamVal.ChiToK892K1680; + if(_K1_1680Hyp || _K1_1680K1_1680Hyp7 || _K1_1680K0_1430Hyp7){ + double K_1_1680Mass=theParamVal.BwK_1_1680.first; - double K_1_1680Width=theParamVal.BwK_1_1680.second; - result+=chiToK1K1Amp(theData, ChiToK892K1680, K892Mass, K892Width, K_1_1680Mass, K_1_1680Width); + double K_1_1680Width=theParamVal.BwK_1_1680.second; + + if (_K1_1680Hyp){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK892K1680=theParamVal.ChiToK892K1680; + result+=chiToK1K1Amp(theData, ChiToK892K1680, K892Mass, K892Width, K_1_1680Mass, K_1_1680Width); + } + + if (_K1_1680K1_1680Hyp7){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK1680K1680=theParamVal.ChiToK1680K1680; + result+=chiToK1K1Amp(theData, ChiToK1680K1680, K_1_1680Mass, K_1_1680Width, K_1_1680Mass, K_1_1680Width); + } + + if(_K1_1680K0_1430Hyp7){ + + double K_0_1430Mass=theParamVal.BwK_0_1430.first; + double K_0_1430Width=theParamVal.BwK_0_1430.second; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK1680K_0_1430=theParamVal.ChiToK1680K_0_1430; + result+=chiToK0K1Amp(theData, ChiToK1680K_0_1430, K_0_1430Mass, K_0_1430Width, K_1_1680Mass, K_1_1680Width); + } + } if(_K1_2300Hyp){ @@ -71,19 +93,15 @@ void Hyp7Lh::setMnUsrParams(MnUserParameters& upar, param2K2PiGam& startVal, pa Hyp6Lh::setMnUsrParams(upar, startVal, errVal); - if(_KappaHyp){ - _fitParams2K2PiGam.setMnUsrParamsDec(upar, startVal, errVal, paramEnum2K2PiGam::KappaKappa); - if (!_KappaK_0_1950Hyp6) _fitParams2K2PiGam.setMnUsrParamsMass(upar, startVal, errVal, paramEnum2K2PiGam::Kappa); - } + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ - if(_K1_1680Hyp){ - _fitParams2K2PiGam.setMnUsrParamsDec(upar, startVal, errVal, paramEnum2K2PiGam::K892K_1_1680); - _fitParams2K2PiGam.setMnUsrParamsMass(upar, startVal, errVal, paramEnum2K2PiGam::K_1_1680); + _fitParams2K2PiGam.setMnUsrParamsDec(upar, startVal, errVal, (*itAmps)); } - if(_K1_2300Hyp){ - _fitParams2K2PiGam.setMnUsrParamsDec(upar, startVal, errVal, paramEnum2K2PiGam::K892K_1_2300); - _fitParams2K2PiGam.setMnUsrParamsMass(upar, startVal, errVal, paramEnum2K2PiGam::K_1_2300); + std::vector<unsigned int>::const_iterator itMasses; + for ( itMasses=_massVec.begin(); itMasses!=_massVec.end(); ++itMasses){ + _fitParams2K2PiGam.setMnUsrParamsMass(upar, startVal, errVal, (*itMasses) ); } } @@ -99,21 +117,19 @@ int Hyp7Lh::setFitParamVal(param2K2PiGam& theParamVal, const std::vector<double> } int counter=Hyp6Lh::setFitParamVal(theParamVal, par); - if(_KappaHyp){ - counter=_fitParams2K2PiGam.setFitParamValDec(theParamVal, par, counter, paramEnum2K2PiGam::KappaKappa); - if (!_KappaK_0_1950Hyp6) counter=_fitParams2K2PiGam.setFitParamValMass(theParamVal, par, counter, paramEnum2K2PiGam::Kappa); - } - if(_K1_1680Hyp){ - counter=_fitParams2K2PiGam.setFitParamValDec(theParamVal, par, counter, paramEnum2K2PiGam::K892K_1_1680); - counter=_fitParams2K2PiGam.setFitParamValMass(theParamVal, par, counter, paramEnum2K2PiGam::K_1_1680); - } - if(_K1_2300Hyp){ - counter=_fitParams2K2PiGam.setFitParamValDec(theParamVal, par, counter, paramEnum2K2PiGam::K892K_1_2300); - counter=_fitParams2K2PiGam.setFitParamValMass(theParamVal, par, counter, paramEnum2K2PiGam::K_1_2300); - } + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ + counter=_fitParams2K2PiGam.setFitParamValDec(theParamVal, par, counter, (*itAmps)); + } + + std::vector<unsigned int>::const_iterator itMasses; + for ( itMasses=_massVec.begin(); itMasses!=_massVec.end(); ++itMasses){ + counter=_fitParams2K2PiGam.setFitParamValMass(theParamVal, par, counter, (*itMasses) ); + } return counter; + } unsigned int Hyp7Lh::nFitParams(){ @@ -204,7 +220,26 @@ void Hyp7Lh::setUp(const std::map<const std::string, bool>& hypMap){ _hypMap[iter->first]= iter->second; } else Alert << "hypothesis K1_1680Hyp7 not set!!!" <<endmsg; - + + + iter= hypMap.find("K1_1680K1_1680Hyp7"); + if (iter !=hypMap.end()){ + _K1_1680K1_1680Hyp7= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K1_1680K1_1680Hyp7 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K1_1680K1_1680Hyp7 not set!!!" <<endmsg; + + + iter= hypMap.find("K1_1680K0_1430Hyp7"); + if (iter !=hypMap.end()){ + _K1_1680K0_1430Hyp7= iter->second; + Info<< "hypothesis " << iter->first << "\t" << _K1_1680K0_1430Hyp7 <<endmsg; + _hypMap[iter->first]= iter->second; + } + else Alert << "hypothesis K1_1680K0_1430Hyp7 not set!!!" <<endmsg; + + iter= hypMap.find("K1_2300Hyp7"); if (iter !=hypMap.end()){ @@ -221,9 +256,20 @@ void Hyp7Lh::setUp(const std::map<const std::string, bool>& hypMap){ if(!_KappaK_0_1950Hyp6) _massVec.push_back(paramEnum2K2PiGam::Kappa); } - if(_K1_1680Hyp){ - _ampVec.push_back(paramEnum2K2PiGam::K892K_1_1680); + if(_K1_1680Hyp || _K1_1680K1_1680Hyp7 || _K1_1680K0_1430Hyp7){ + _massVec.push_back(paramEnum2K2PiGam::K_1_1680); + if(_K1_1680Hyp){ + _ampVec.push_back(paramEnum2K2PiGam::K892K_1_1680); + } + if(_K1_1680K1_1680Hyp7){ + _ampVec.push_back(paramEnum2K2PiGam::ChiToK1680K1680); + } + + if(_K1_1680K0_1430Hyp7){ + _ampVec.push_back(paramEnum2K2PiGam::ChiToK1680K_0_1430); + } + } if(_K1_2300Hyp){ diff --git a/Examples/Psi2STo2K2PiGam/Hyp7Lh.hh b/Examples/Psi2STo2K2PiGam/Hyp7Lh.hh index 015234c01246d0ed915820a18efe0801b2f62f7c..27ea4ebc2b26836c671f5f87d263a2b1ded11f69 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp7Lh.hh +++ b/Examples/Psi2STo2K2PiGam/Hyp7Lh.hh @@ -52,6 +52,8 @@ public: protected: bool _KappaHyp; bool _K1_1680Hyp; + bool _K1_1680K1_1680Hyp7; + bool _K1_1680K0_1430Hyp7; bool _K1_2300Hyp; virtual complex<double> chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData); diff --git a/Examples/Psi2STo2K2PiGam/Hyp8Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp8Lh.cc index 728923e854021cb3a96bc9fc5d7dc75ba08cd1e6..84b99ec73314425998bfd2f21f27f5ccd1c89ff1 100644 --- a/Examples/Psi2STo2K2PiGam/Hyp8Lh.cc +++ b/Examples/Psi2STo2K2PiGam/Hyp8Lh.cc @@ -267,7 +267,7 @@ void Hyp8Lh::setUp(const std::map<const std::string, bool>& hypMap){ _massVec.push_back(paramEnum2K2PiGam::K_1_1650); } - if(!_K0_1430_K0_1430Hyp && !_K0_1430_K0_1430Hyp && !_K1_1270Hyp) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); + if(!_K0_1430_K0_1430Hyp && !_K0_1430_K0_1430Hyp && !_K1_1270Hyp && !_K0_1430_K892Hyp1) _massVec.push_back(paramEnum2K2PiGam::K_0_1430); std::vector<unsigned int>::iterator ampIt; diff --git a/Examples/Psi2STo2K2PiGam/HypProdLh.cc b/Examples/Psi2STo2K2PiGam/HypProdLh.cc new file mode 100644 index 0000000000000000000000000000000000000000..c6e70112b7feec737ee6590ea698f5d7ea964276 --- /dev/null +++ b/Examples/Psi2STo2K2PiGam/HypProdLh.cc @@ -0,0 +1,145 @@ +#include <getopt.h> +#include <fstream> + +#include "Examples/Psi2STo2K2PiGam/HypProdLh.hh" +#include "Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.hh" +#include "ErrLogger/ErrLogger.hh" + +HypProdLh::HypProdLh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList> theEvtList) : + AbsPsi2STo2K2PiGamLh(theEvtList) + ,_nFitParams(0) +{ + setUp(); +} + +HypProdLh::HypProdLh( boost::shared_ptr<AbsPsi2STo2K2PiGamLh> theLhPtr) : + AbsPsi2STo2K2PiGamLh(theLhPtr->getEventList()) + ,_nFitParams(0) +{ + setUp(); +} + +HypProdLh::~HypProdLh() +{; +} + + +complex<double> HypProdLh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData){ + + complex<double> result(1.,0.); + + return result; + +} + + + +void HypProdLh::setMnUsrParams(MnUserParameters& upar, param2K2PiGam& startVal, param2K2PiGam& errVal){ + + _fitParams2K2PiGam.setMnUsrParamsFlatteMass(upar, startVal, errVal, paramEnum2K2PiGam::name(paramEnum2K2PiGam::phaseSpace)); + + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ + + _fitParams2K2PiGam.setMnUsrParamsDec(upar, startVal, errVal, (*itAmps)); + } + + + +} + + + +int HypProdLh::setFitParamVal(param2K2PiGam& theParamVal, const std::vector<double>& par){ + + if (par.size() != nFitParams() ) { + Alert << "size of parameters wrong!!! par.size()=" << par.size() << + "\t it should be" << nFitParams() << endmsg; + exit(1); + } + + std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; + int counter=0; + + counter=_fitParams2K2PiGam.setFitParamFlatteMass(theParamVal, par, counter, paramEnum2K2PiGam::name(paramEnum2K2PiGam::phaseSpace)); + + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ + counter=_fitParams2K2PiGam.setFitParamValDec(theParamVal, par, counter, (*itAmps)); + } + + return counter; +} + +unsigned int HypProdLh::nFitParams(){ + + return _nFitParams; +} + + +void HypProdLh::print(std::ostream& os) const{ + os << "HypProdLh::print\n"; +} + +void HypProdLh::printCurrentFitResult(param2K2PiGam& theParamVal){ + + DebugMsg<< "phaseSpace: " << theParamVal.phaseSpace << endmsg; + + + std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; + + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ + std::vector< boost::shared_ptr<const JPCLS> > JPCLSs=_fitParams2K2PiGam.jpclsVec(*itAmps); + + for ( itJPCLS=JPCLSs.begin(); itJPCLS!=JPCLSs.end(); ++itJPCLS){ + DebugMsg<< (*itJPCLS)->name()<< paramEnum2K2PiGam::name(*itAmps) << endmsg; + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > currentMap=_fitParams2K2PiGam.ampMap(theParamVal, *itAmps); + std::pair<double, double> tmpParam=currentMap[(*itJPCLS)]; + DebugMsg <<"\t mag:" << tmpParam.first <<"\t phi:" << tmpParam.second << endmsg; + } + } + +} + +void HypProdLh::dumpCurrentResult(std::ostream& os, param2K2PiGam& theParamVal, std::string& suffix){ + + if ( suffix.compare("Val") != 0 && suffix.compare("Err") !=0 ){ + Warning << "suffix " << suffix << " not supported!!! Use Val or Err" << endmsg; + return; + } + + std::string tmpStringPs="phaseSpace"+suffix; + os << tmpStringPs << "\t" << theParamVal.phaseSpace << "\t" << 0. << std::endl; + + std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; + + std::vector<unsigned int>::const_iterator itAmps; + for ( itAmps=_ampVec.begin(); itAmps!=_ampVec.end(); ++itAmps){ + std::vector< boost::shared_ptr<const JPCLS> > JPCLSs=_fitParams2K2PiGam.jpclsVec(*itAmps); + + for ( itJPCLS=JPCLSs.begin(); itJPCLS!=JPCLSs.end(); ++itJPCLS){ + std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > currentMap=_fitParams2K2PiGam.ampMap(theParamVal, *itAmps); + std::pair<double, double> tmpParam=currentMap[(*itJPCLS)]; + + std::string tmpStringDec=(*itJPCLS)->name()+paramEnum2K2PiGam::name(*itAmps)+suffix; + os << tmpStringDec << "\t" << tmpParam.first << "\t" << tmpParam.second << std::endl; + } + } + +} + + +void HypProdLh::setUp(){ + + _ampVec.push_back(paramEnum2K2PiGam::ChiGam); + + std::vector<unsigned int>::iterator ampIt; + for (ampIt=_ampVec.begin(); ampIt!=_ampVec.end(); ++ampIt){ + std::vector< boost::shared_ptr<const JPCLS> > JPCLSs=_fitParams2K2PiGam.jpclsVec(*ampIt); + _nFitParams+=2*JPCLSs.size(); + } + + _nFitParams+=1; //phase space + +} diff --git a/Examples/Psi2STo2K2PiGam/HypProdLh.hh b/Examples/Psi2STo2K2PiGam/HypProdLh.hh new file mode 100644 index 0000000000000000000000000000000000000000..ffa9c4d5fb82ee3548acb2076e7c82e7abcdf71c --- /dev/null +++ b/Examples/Psi2STo2K2PiGam/HypProdLh.hh @@ -0,0 +1,64 @@ +#ifndef _HypProdLh_H +#define _HypProdLh_H + +#include <iostream> +#include <fstream> +#include <string> +#include <vector> +#include <complex> +//#include <map> + +#include <cassert> +#include <boost/shared_ptr.hpp> + +#include "TROOT.h" +// #include <TSystem.h> +#include "qft++/topincludes/relativistic-quantum-mechanics.hh" + +#include "Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh" +#include "Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamData.hh" +#include "PwaUtils/DataUtils.hh" + +#include "Minuit2/MnUserParameters.h" + + + +class HypProdLh : public AbsPsi2STo2K2PiGamLh{ + +public: + + // create/copy/destroy: + + ///Constructor + HypProdLh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList>); + HypProdLh(boost::shared_ptr<AbsPsi2STo2K2PiGamLh>); + + /** Destructor */ + virtual ~HypProdLh(); + + virtual AbsPsi2STo2K2PiGamLh* clone_() const{ + return new HypProdLh(_Psi2STo2K2PiGamEvtListPtr); + } + + + // Getters: + virtual void setMnUsrParams(MnUserParameters& upar, param2K2PiGam& startVal, param2K2PiGam& errVal); + virtual int setFitParamVal(param2K2PiGam& theParamVal, const std::vector<double>& par); + virtual unsigned int nFitParams(); + + virtual void print(std::ostream& os) const; + virtual void printCurrentFitResult(param2K2PiGam& theParamVal); + virtual void dumpCurrentResult(std::ostream& os, param2K2PiGam& theParamVal, std::string& suffix); + + +protected: + + virtual complex<double> chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData); + +private: + unsigned int _nFitParams; + std::vector<unsigned int> _ampVec; + void setUp(); +}; + +#endif diff --git a/Examples/Psi2STo2K2PiGam/Jamfile b/Examples/Psi2STo2K2PiGam/Jamfile index f0a330b8a7290fd63dfb68f412e826f18d5e4969..35c51f5fb87c22268f0779b84d1746d0ac9f576d 100644 --- a/Examples/Psi2STo2K2PiGam/Jamfile +++ b/Examples/Psi2STo2K2PiGam/Jamfile @@ -13,5 +13,5 @@ lib Psi2STo2K2PiGam : exe ParserTestApp : ParserTestApp.cc Psi2STo2K2PiGam : ; exe Mpsi2STo2K2PiGamTestApp : Mpsi2STo2K2PiGamTestApp.cc Psi2STo2K2PiGam : ; -exe qaMixedSampleApp : qaMixedSampleApp.cc Psi2STo2K2PiGam : ; exe qaMixedSampleCacheApp : qaMixedSampleCacheApp.cc Psi2STo2K2PiGam : ; +exe qaMixedSampleFastApp : qaMixedSampleFastApp.cc Psi2STo2K2PiGam : ; diff --git a/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc b/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc index c0c82926e942a947fd6cfdfb4ec379cbe7c26941..15d1c3e464a479245ccf95f59f04fd5c01e9b5e5 100644 --- a/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc +++ b/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc @@ -13,6 +13,7 @@ #include "Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh" #include "Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.hh" #include "Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh" +#include "Examples/Psi2STo2K2PiGam/HypProdLh.hh" #include "Examples/Psi2STo2K2PiGam/Hyp1Lh.hh" #include "Examples/Psi2STo2K2PiGam/Hyp2Lh.hh" #include "Examples/Psi2STo2K2PiGam/Hyp3Lh.hh" @@ -147,7 +148,11 @@ int main(int __argc,char *__argv[]){ hypMap["K0_1430_K0_1430Hyp"]=true; hypMap["K2_1430_K2_1430Hyp"]=true; hypMap["K0_1430_K2_1430Hyp"]=true; + hypMap["K0_1430_K892Hyp1"]=true; + hypMap["K2_1430_K892Hyp1"]=true; hypMap["K1_1410_K1_1410Hyp"]=true; + hypMap["K1_1410_K892Hyp1"]=true; + hypMap["f1710_f1710Hyp1"]=true; hypMap["doHyp2"]=true; hypMap["sigmaf980Hyp3"]=true; hypMap["sigmaf1710Hyp3"]=true; @@ -167,12 +172,18 @@ int main(int __argc,char *__argv[]){ hypMap["K_0_2400KHyp5"]=true; hypMap["K_0_2400KTof_0_1710Hyp5"]=true; hypMap["K_1_2400KHyp5"]=true; + hypMap["K_1_2400KTof_0_1710Hyp5"]=true; + hypMap["K_1_2400KToK_0_1430Hyp5"]=true; + hypMap["K_2_2400KTof980Hyp5"]=true; + hypMap["K_2_2400KTof_0_1710Hyp5"]=true; hypMap["ChiToK_0_1430KPiHyp5"]=true; hypMap["ChiToK892KPiHyp5"]=true; hypMap["K_0_1430K_0_1950Hyp6"]=true; hypMap["KappaK_0_1430Hyp6"]=true; hypMap["KappaK_0_1950Hyp6"]=true; hypMap["K1_1680Hyp7"]=true; + hypMap["K1_1680K1_1680Hyp7"]=true; + hypMap["K1_1680K0_1430Hyp7"]=true; hypMap["K1_2300Hyp7"]=true; hypMap["KappaHyp7"]=true; hypMap["K_0_1460ToKstPiHyp8"]=true; @@ -201,7 +212,7 @@ int main(int __argc,char *__argv[]){ std::map<const std::string, bool>::const_iterator iter= hypMap.find( (*itStr) ); if (iter !=hypMap.end()){ hypMap[iter->first]= false; - Info<< "hypothesis " << iter->first << " disabed" <<endmsg; + Info<< "hypothesis " << iter->first << " disabled" <<endmsg; } else { Alert << "hypothesis " << (*itStr) << " can not be disabled"<<endmsg; exit(0); @@ -211,8 +222,8 @@ int main(int __argc,char *__argv[]){ boost::shared_ptr<AbsPsi2STo2K2PiGamLh> thePsi2STo2K2PiGamLhPtr; std::string startWithHyp=theAppParams.startHypo(); - - if (startWithHyp=="hyp1") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp1Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); + if (startWithHyp=="prod") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new HypProdLh(thePsi2STo2K2PiGamEvtListPtr)); + else if (startWithHyp=="hyp1") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp1Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); else if (startWithHyp=="hyp2") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp2Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); else if (startWithHyp=="hyp3") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp3Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); else if (startWithHyp=="hyp4") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp4Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); @@ -221,7 +232,7 @@ int main(int __argc,char *__argv[]){ else if (startWithHyp=="hyp7") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp7Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); else if (startWithHyp=="hyp8") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp8Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); else if (startWithHyp=="hyp9") thePsi2STo2K2PiGamLhPtr= boost::shared_ptr<AbsPsi2STo2K2PiGamLh>(new Hyp9Lh(thePsi2STo2K2PiGamEvtListPtr, hypMap)); - else { Alert << "start with hypthesis " << startWithHyp << " not supported!!!!" << endmsg; + else { Alert << "start with hypothesis " << startWithHyp << " not supported!!!!" << endmsg; exit(1); } @@ -249,7 +260,7 @@ int main(int __argc,char *__argv[]){ if (qaMode){ thePsi2STo2K2PiGamLhPtr->printCurrentFitResult(theStartparams); double theLh=thePsi2STo2K2PiGamLhPtr->calcLogLh(theStartparams); - Info <<"theLh = "<< theLh << endmsg; + Info << "theLh = " << theLh << endmsg; Psi2STo2K2PiGamHist Psi2STo2K2PiGamHist(thePsi2STo2K2PiGamLhPtr, theStartparams); return 0; diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamData.hh b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamData.hh index 5c4c3fecf937994512e01be641d865e0387632a5..0d07e6be009af6d7c58bd09ae207659a7b3dc9a1 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamData.hh +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamData.hh @@ -112,9 +112,8 @@ namespace Psi2STo2K2PiGamData { map<Spin,map<Spin,map<Spin,complex<double> > > > Dff2ToKKviaKKpi0; map<Spin,map<Spin,map<Spin,complex<double> > > > Dff2ToKKviaKKpi1; - map<Spin,map<Spin,map<Spin,complex<double> > > > DfK1pTof0Kp; - map<Spin,map<Spin,map<Spin,complex<double> > > > DfK1mTof0Km; - + map<Spin,map<Spin,map<Spin,complex<double> > > > DfKjpTof0Kp; + map<Spin,map<Spin,map<Spin,complex<double> > > > DfKjmTof0Km; bool operator<(const Psi2STo2K2PiGamEvtData& compare) const{ bool result=false; diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc index 621734fa315d117d057d3b578feb21cd995f41ec..71cca09ba4fb161dc6b6943748dde7c8cfe3660e 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc @@ -170,15 +170,14 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int // calculate and store WignerD functions for K*+(1400)->K*+(892) pi1; K*+(892)->K+ pi0 + cc Spin jKst1400=1; - for (Spin M=-1; M<=1; M++){ - for (Spin lam=-1; lam<=1; lam++){ - thePsi2STo2K2PiGamEvtData->DfK1400pToKstpPi1[jKst1400][M][lam]=Wigner_D(KpPi0_HeliKpPi0Pi0_4V.Phi(),KpPi0_HeliKpPi0Pi0_4V.Theta(),0,jKst1400,M,lam); - thePsi2STo2K2PiGamEvtData->DfK1400pToKstpPi0[jKst1400][M][lam]=Wigner_D(KpPi1_HeliKpPi0Pi0_4V.Phi(),KpPi1_HeliKpPi0Pi0_4V.Theta(),0,jKst1400,M,lam); - - thePsi2STo2K2PiGamEvtData->DfK1400mToKstmPi1[jKst1400][M][lam]=Wigner_D(KmPi0_HeliKmPi0Pi0_4V.Phi(),KmPi0_HeliKmPi0Pi0_4V.Theta(),0,jKst1400,M,lam); - thePsi2STo2K2PiGamEvtData->DfK1400mToKstmPi0[jKst1400][M][lam]=Wigner_D(KmPi1_HeliKmPi0Pi0_4V.Phi(),KmPi1_HeliKmPi0Pi0_4V.Theta(),0,jKst1400,M,lam); - } + for (Spin lam=-1; lam<=1; lam++){ + thePsi2STo2K2PiGamEvtData->DfK1400pToKstpPi1[jKst1400][0][lam]=Wigner_D(KpPi0_HeliKpPi0Pi0_4V.Phi(),KpPi0_HeliKpPi0Pi0_4V.Theta(),0,jKst1400,0,lam); + thePsi2STo2K2PiGamEvtData->DfK1400pToKstpPi0[jKst1400][0][lam]=Wigner_D(KpPi1_HeliKpPi0Pi0_4V.Phi(),KpPi1_HeliKpPi0Pi0_4V.Theta(),0,jKst1400,0,lam); + + thePsi2STo2K2PiGamEvtData->DfK1400mToKstmPi1[jKst1400][0][lam]=Wigner_D(KmPi0_HeliKmPi0Pi0_4V.Phi(),KmPi0_HeliKmPi0Pi0_4V.Theta(),0,jKst1400,0,lam); + thePsi2STo2K2PiGamEvtData->DfK1400mToKstmPi0[jKst1400][0][lam]=Wigner_D(KmPi1_HeliKmPi0Pi0_4V.Phi(),KmPi1_HeliKmPi0Pi0_4V.Theta(),0,jKst1400,0,lam); } + // calculate and store WignerD functions for K2->K*2 pi; Spin jK2=2; @@ -192,12 +191,10 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int } // calculate and store WignerD functions for K*2->K pi0 via K pi0 pi0 for (Spin M=-2; M<=2; M++){ - for (Spin lam2=-2; lam2<=2; lam2++){ - thePsi2STo2K2PiGamEvtData->DfK2pToKpPi0ViaKpPiPi[jK2][M][lam2]=Wigner_D(Kp_HeliKpPi0_ViaKpPiPi_4V.Phi(), Kp_HeliKpPi0_ViaKpPiPi_4V.Theta(),0,jK2,M,lam2); - thePsi2STo2K2PiGamEvtData->DfK2pToKpPi1ViaKpPiPi[jK2][M][lam2]=Wigner_D(Kp_HeliKpPi1_ViaKpPiPi_4V.Phi(),Kp_HeliKpPi1_ViaKpPiPi_4V.Theta(),0,jK2,M,lam2); - thePsi2STo2K2PiGamEvtData->DfK2mToKmPi0ViaKmPiPi[jK2][M][lam2]=Wigner_D(Km_HeliKmPi0_ViaKmPiPi_4V.Phi(),Km_HeliKmPi0_ViaKmPiPi_4V.Theta(),0,jK2,M,lam2); - thePsi2STo2K2PiGamEvtData->DfK2mToKmPi1ViaKmPiPi[jK2][M][lam2]=Wigner_D(Km_HeliKmPi1_ViaKmPiPi_4V.Phi(),Km_HeliKmPi1_ViaKmPiPi_4V.Theta(),0,jK2,M,lam2); - } + thePsi2STo2K2PiGamEvtData->DfK2pToKpPi0ViaKpPiPi[jK2][M][0]=Wigner_D(Kp_HeliKpPi0_ViaKpPiPi_4V.Phi(), Kp_HeliKpPi0_ViaKpPiPi_4V.Theta(),0,jK2,M,0); + thePsi2STo2K2PiGamEvtData->DfK2pToKpPi1ViaKpPiPi[jK2][M][0]=Wigner_D(Kp_HeliKpPi1_ViaKpPiPi_4V.Phi(),Kp_HeliKpPi1_ViaKpPiPi_4V.Theta(),0,jK2,M,0); + thePsi2STo2K2PiGamEvtData->DfK2mToKmPi0ViaKmPiPi[jK2][M][0]=Wigner_D(Km_HeliKmPi0_ViaKmPiPi_4V.Phi(),Km_HeliKmPi0_ViaKmPiPi_4V.Theta(),0,jK2,M,0); + thePsi2STo2K2PiGamEvtData->DfK2mToKmPi1ViaKmPiPi[jK2][M][0]=Wigner_D(Km_HeliKmPi1_ViaKmPiPi_4V.Phi(),Km_HeliKmPi1_ViaKmPiPi_4V.Theta(),0,jK2,M,0); } // calculate and store WignerD functions for K*1+->K+ pi0, K*1+->K+ pi1+ cc Spin jKst892=1; @@ -211,19 +208,11 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int // calculate and store WignerD functions for K*1+->K+ pi0, K*1+->K+ pi1+ cc Via K+ pi0 pi0 for (Spin M=-1; M<=1; M++){ - for (Spin lam2=-1; lam2<=1; lam2++){ - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKpPiPi[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi0_ViaKpPiPi_4V.Phi(), Kp_HeliKpPi0_ViaKpPiPi_4V.Theta(),0,jKst892,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKpPiPi[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi1_ViaKpPiPi_4V.Phi(),Kp_HeliKpPi1_ViaKpPiPi_4V.Theta(),0,jKst892,M,lam2); - - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKmPiPi[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi0_ViaKmPiPi_4V.Phi(),Km_HeliKmPi0_ViaKmPiPi_4V.Theta(),0,jKst892,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKmPiPi[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi1_ViaKmPiPi_4V.Phi(),Km_HeliKmPi1_ViaKmPiPi_4V.Theta(),0,jKst892,M,lam2); + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKpPiPi[jKst892][M][0]=Wigner_D(Kp_HeliKpPi0_ViaKpPiPi_4V.Phi(), Kp_HeliKpPi0_ViaKpPiPi_4V.Theta(),0,jKst892,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKpPiPi[jKst892][M][0]=Wigner_D(Kp_HeliKpPi1_ViaKpPiPi_4V.Phi(),Kp_HeliKpPi1_ViaKpPiPi_4V.Theta(),0,jKst892,M,0); -// thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKpPiPi[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi0_4V.Phi(), Kp_HeliKpPi0_4V.Theta(),0,jKst892,M,lam2); -// thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKpPiPi[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi1_4V.Phi(), Kp_HeliKpPi1_4V.Theta(),0,jKst892,M,lam2); - -// thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKmPiPi[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi0_4V.Phi(), Km_HeliKmPi0_4V.Theta(),0,jKst892,M,lam2); -// thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKmPiPi[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi1_4V.Phi(), Km_HeliKmPi1_4V.Theta(),0,jKst892,M,lam2); - } + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKmPiPi[jKst892][M][0]=Wigner_D(Km_HeliKmPi0_ViaKmPiPi_4V.Phi(),Km_HeliKmPi0_ViaKmPiPi_4V.Theta(),0,jKst892,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKmPiPi[jKst892][M][0]=Wigner_D(Km_HeliKmPi1_ViaKmPiPi_4V.Phi(),Km_HeliKmPi1_ViaKmPiPi_4V.Theta(),0,jKst892,M,0); } @@ -273,28 +262,24 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int // calculate and store WignerD functions for K*1+->K+ pi0, K*1+->K+ pi1+ cc Via K+ K- pi0 - for (Spin M=-1; M<=1; M++){ - for (Spin lam2=-1; lam2<=1; lam2++){ - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKKPi0[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi0_ViaKKPi0_4V.Phi(), Kp_HeliKpPi0_ViaKKPi0_4V.Theta(),0,jKst892,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKKPi1[jKst892][M][lam2]=Wigner_D(Kp_HeliKpPi1_ViaKKPi1_4V.Phi(),Kp_HeliKpPi1_ViaKKPi1_4V.Theta(),0,jKst892,M,lam2); - - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKKPi0[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi0_ViaKKPi0_4V.Phi(),Km_HeliKmPi0_ViaKKPi0_4V.Theta(),0,jKst892,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKKPi1[jKst892][M][lam2]=Wigner_D(Km_HeliKmPi1_ViaKKPi1_4V.Phi(),Km_HeliKmPi1_ViaKKPi1_4V.Theta(),0,jKst892,M,lam2); - + for (Spin M=-1; M<=1; M++){ + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKKPi0[jKst892][M][0]=Wigner_D(Kp_HeliKpPi0_ViaKKPi0_4V.Phi(), Kp_HeliKpPi0_ViaKKPi0_4V.Theta(),0,jKst892,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKKPi1[jKst892][M][0]=Wigner_D(Kp_HeliKpPi1_ViaKKPi1_4V.Phi(),Kp_HeliKpPi1_ViaKKPi1_4V.Theta(),0,jKst892,M,0); + + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKKPi0[jKst892][M][0]=Wigner_D(Km_HeliKmPi0_ViaKKPi0_4V.Phi(),Km_HeliKmPi0_ViaKKPi0_4V.Theta(),0,jKst892,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKKPi1[jKst892][M][0]=Wigner_D(Km_HeliKmPi1_ViaKKPi1_4V.Phi(),Km_HeliKmPi1_ViaKKPi1_4V.Theta(),0,jKst892,M,0); + } - } - -// calculate and store WignerD functions for K*2+->K+ pi0, K*2+->K+ pi1+ cc Via K+ K- pi0 - for (Spin M=-2; M<=2; M++){ - for (Spin lam2=-2; lam2<=2; lam2++){ - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKKPi0[2][M][lam2]=Wigner_D(Kp_HeliKpPi0_ViaKKPi0_4V.Phi(), Kp_HeliKpPi0_ViaKKPi0_4V.Theta(),0,2,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKKPi1[2][M][lam2]=Wigner_D(Kp_HeliKpPi1_ViaKKPi1_4V.Phi(),Kp_HeliKpPi1_ViaKKPi1_4V.Theta(),0,2,M,lam2); - - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKKPi0[2][M][lam2]=Wigner_D(Km_HeliKmPi0_ViaKKPi0_4V.Phi(),Km_HeliKmPi0_ViaKKPi0_4V.Theta(),0,2,M,lam2); - thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKKPi1[2][M][lam2]=Wigner_D(Km_HeliKmPi1_ViaKKPi1_4V.Phi(),Km_HeliKmPi1_ViaKKPi1_4V.Theta(),0,2,M,lam2); - + + // calculate and store WignerD functions for K*2+->K+ pi0, K*2+->K+ pi1+ cc Via K+ K- pi0 + for (Spin M=-2; M<=2; M++){ + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKKPi0[2][M][0]=Wigner_D(Kp_HeliKpPi0_ViaKKPi0_4V.Phi(), Kp_HeliKpPi0_ViaKKPi0_4V.Theta(),0,2,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi1ViaKKPi1[2][M][0]=Wigner_D(Kp_HeliKpPi1_ViaKKPi1_4V.Phi(),Kp_HeliKpPi1_ViaKKPi1_4V.Theta(),0,2,M,0); + + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKKPi0[2][M][0]=Wigner_D(Km_HeliKmPi0_ViaKKPi0_4V.Phi(),Km_HeliKmPi0_ViaKKPi0_4V.Theta(),0,2,M,0); + thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKKPi1[2][M][0]=Wigner_D(Km_HeliKmPi1_ViaKKPi1_4V.Phi(),Km_HeliKmPi1_ViaKKPi1_4V.Theta(),0,2,M,0); + } - } // calculate and store WignerD functions for K*0+->K+ pi0, K*0+->K+ pi1+ cc Via K+ K- pi0 thePsi2STo2K2PiGamEvtData->DfKst1pToKpPi0ViaKKPi0[0][0][0]=Wigner_D(Kp_HeliKpPi0_ViaKKPi0_4V.Phi(), Kp_HeliKpPi0_ViaKKPi0_4V.Theta(),0,0,0,0); @@ -303,12 +288,12 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi0ViaKKPi0[0][0][0]=Wigner_D(Km_HeliKmPi0_ViaKKPi0_4V.Phi(),Km_HeliKmPi0_ViaKKPi0_4V.Theta(),0,0,0,0); thePsi2STo2K2PiGamEvtData->DfKst1mToKmPi1ViaKKPi1[0][0][0]=Wigner_D(Km_HeliKmPi1_ViaKKPi1_4V.Phi(),Km_HeliKmPi1_ViaKKPi1_4V.Theta(),0,0,0,0); - for (Spin lam1=-1; lam1<=1; lam1++){ - for (Spin lam2=-1; lam2<=1; lam2++){ - thePsi2STo2K2PiGamEvtData->DfK1pTof0Kp[1][lam1][lam2]=Wigner_D(PiPi_HeliKpPi0Pi0_4V.Phi(), PiPi_HeliKpPi0Pi0_4V.Theta(),0, 1, lam1,lam2); - thePsi2STo2K2PiGamEvtData->DfK1mTof0Km[1][lam1][lam2]=Wigner_D(PiPi_HeliKmPi0Pi0_4V.Phi(), PiPi_HeliKmPi0Pi0_4V.Theta(),0, 1, lam1,lam2); - } + + for (Spin Kj=0; Kj<=3; Kj++){ + thePsi2STo2K2PiGamEvtData->DfKjpTof0Kp[Kj][0][0]=Wigner_D(PiPi_HeliKpPi0Pi0_4V.Phi(), PiPi_HeliKpPi0Pi0_4V.Theta(),0, Kj, 0,0); + thePsi2STo2K2PiGamEvtData->DfKjmTof0Km[Kj][0][0]=Wigner_D(PiPi_HeliKmPi0Pi0_4V.Phi(), PiPi_HeliKmPi0Pi0_4V.Theta(),0, Kj, 0,0); } + return thePsi2STo2K2PiGamEvtData; } diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc index 8ed3a6b8d842cc86fa279706570a5de8e6bc0a5a..9bb2bf0c85370eab67f514b861058e39a1313695 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc @@ -192,11 +192,11 @@ void Psi2STo2K2PiGamHist::initRootStuff(const std::string& fileName) initHistMap(); - _dataTuple = new TNtuple("dataTuple","dataTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipiViaK1pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype"); - - _mcTuple = new TNtuple("mcTuple","mcTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipiViaK1pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype"); + _dataTuple = new TNtuple("dataTuple","dataTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1k2pi1:mk1k2pi2:mk1k2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetaK1ViaK1K2:costhetaKKViaK1K2pi1:costhetaKKViaK1K2pi2:costhetapipiViaK1pipi:costhetapipiViaK2pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype"); + // _mcTuple = new TNtuple("mcTuple","mcTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mK1K2pi1:mK1K2pi2:mK1K2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetapipiViaK1pipi:costhetapipiViaK2pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype"); + _mcTuple = new TNtuple("mcTuple","mcTuple","phiKKpipi:costhetaKKpipi:phiK1pi1pi2:costhetaK1pi1pi2:phiK2pi1pi2:costhetaK2pi1pi2:phiK1pi1:costhetaK1pi1:phiK1pi2:costhetaK1pi2:phiK2pi1:costhetaK2pi1:phiK2pi2:costhetaK2pi2:phipi1:costhetapi1:phipi2:costhetapi2:mk1k2pi1:mk1k2pi2:mk1k2:mk1pi1pi2:mk2pi1pi2:mk1pi1:mk1pi2:mk2pi1:mk2pi2:mpipi:costhetaK1ViaK1K2:costhetaKKFromK1K2pi1:costhetaKKFromK1K2pi2:costhetapipiViaK1pipi:costhetapipiViaK2pipi:phipipiViaK1pipi:costhetapiViapipi:phipiViapipi:weight:datatype"); } void Psi2STo2K2PiGamHist::initHistMap(){ @@ -319,6 +319,9 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv float CosThetapi1 = -theData->Kp_HeliKpPi0_4V.CosTheta(); float Phipi2 = -theData->Kp_HeliKpPi1_4V.Phi(); float CosThetapi2 = -theData->Kp_HeliKpPi1_4V.CosTheta(); + float mK1K2pi1 = theData->KKPi0_HeliChic0_4V.M(); + float mK1K2pi2 = theData->KKPi1_HeliChic0_4V.M(); + float mK1K2 = theData->KpKm_HeliChic0_4V.M(); float mK1pi1pi2 = theData->KpPiPi_HeliChic0_4V.M(); float mK2pi1pi2 = theData->KmPiPi_HeliChic0_4V.M(); float mK1pi1 = theData->KpPi0_HeliChic0_4V.M(); @@ -327,7 +330,11 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv float mK2pi2 = theData->KmPi1_HeliChic0_4V.M(); float mpipi = theData->PiPi_HeliChic0_4V.M(); + float CosThetaK1FromK1K2 = theData->Km_HeliKmKp_4V.CosTheta(); + float CosThetaKKFromK1K2Pi1 = theData->KK_HeliKKPi0_4V.CosTheta(); + float CosThetaKKFromK1K2Pi2 = theData->KK_HeliKKPi1_4V.CosTheta(); float CosThetaPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.CosTheta(); + float CosThetaPiPiFromK2PiPi = theData->PiPi_HeliKmPi0Pi0_4V.CosTheta(); float PhiPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.Phi(); float CosThetaPiFromPiPi = theData->Pi0_HeliPi0Pi0_4V.CosTheta(); @@ -355,21 +362,28 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv fourVectors[15] = CosThetapi1; fourVectors[16] = Phipi2; fourVectors[17] = CosThetapi2; - fourVectors[18] = mK1pi1pi2; - fourVectors[19] = mK2pi1pi2; - fourVectors[20] = mK1pi1; - fourVectors[21] = mK1pi2; - fourVectors[22] = mK2pi1; - fourVectors[23] = mK2pi2; - - fourVectors[24] = mpipi; - fourVectors[25] = CosThetaPiPiFromK1PiPi; - fourVectors[26] = PhiPiPiFromK1PiPi; - fourVectors[27] = CosThetaPiFromPiPi; - fourVectors[28] = PhiPiFromPiPi; - - fourVectors[29] = evtweight; - fourVectors[30] = datatype; + fourVectors[18] = mK1K2pi1; + fourVectors[19] = mK1K2pi2; + fourVectors[20] = mK1K2; + fourVectors[21] = mK1pi1pi2; + fourVectors[22] = mK2pi1pi2; + fourVectors[23] = mK1pi1; + fourVectors[24] = mK1pi2; + fourVectors[25] = mK2pi1; + fourVectors[26] = mK2pi2; + fourVectors[27] = mpipi; + fourVectors[28] = CosThetaK1FromK1K2; + fourVectors[29] = CosThetaKKFromK1K2Pi1; + fourVectors[30] = CosThetaKKFromK1K2Pi2; + fourVectors[31] = CosThetaPiPiFromK1PiPi; + fourVectors[32] = CosThetaPiPiFromK2PiPi; + fourVectors[33] = PhiPiPiFromK1PiPi; + fourVectors[34] = CosThetaPiFromPiPi; + fourVectors[35] = PhiPiFromPiPi; + + fourVectors[36] = evtweight; + fourVectors[37] = datatype; + // cout << evtweight << endl; theTuple->Fill(fourVectors); diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamParser.hh b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamParser.hh index 426d405d4d6f2916647259853b6ca367aa3113e5..1e54ead466be8795231308b1c1407d8cb4c858a7 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamParser.hh +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamParser.hh @@ -79,7 +79,7 @@ class Psi2STo2K2PiGamParser , _qaMode(false) , _genMode(false) { - _disabledHyps.push_back("blainit"); +// _disabledHyps.push_back("blainit"); if (!parseCommandLine(argc, argv)) throw false; } diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc index 6f90d3383d7d8129d2178ec19b6168368a5cb331..82c2f07113c8e026db7640ac3ebc1a2a6dc62022 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc @@ -38,6 +38,12 @@ Psi2STo2K2PiGamStates::Psi2STo2K2PiGamStates() : //LS combinations for the Chi_c0 decay to K*0(1430) K*2(1430) fillJPCLS(_chic0JPC, _Kst0JPC, _Kst2JPC, _ChiToK0K2_JPCLS); + //LS combinations for the Chi_c0 decay to K*0 K*1 + fillJPCLS(_chic0JPC, _Kst0JPC, _Kst1JPC, _ChiToKst0Kst1JPCLS); + + //LS combinations for the Chi_c0 decay to K*1 K*2 + fillJPCLS(_chic0JPC, _Kst1JPC, _Kst2JPC, _ChiToKst1Kst2JPCLS); + //LS combinations for the Chi_c0 decay to f0 f0 fillJPCLS(_chic0JPC, _f0JPC, _f0JPC,_ChiTof0f0JPCLS); @@ -94,6 +100,9 @@ Psi2STo2K2PiGamStates::Psi2STo2K2PiGamStates() : fillJPCLS(_pi0JPC, _Kst1JPC, _pi0JPC, _Pi0pToKstKJPCLS); + + fillJPCLS(_K2mJPC, _f0JPC, _kJPC, _K2Tof0KJPCLS); + } Psi2STo2K2PiGamStates::~Psi2STo2K2PiGamStates() @@ -168,6 +177,18 @@ void Psi2STo2K2PiGamStates::print(std::ostream& os) const os << "\n" << std::endl; } + os << "*** Chi_c0: LS combinations for the decay to K*0 K*1 *** "<< std::endl; + for ( itJPCLS=_ChiToKst0Kst1JPCLS.begin(); itJPCLS!=_ChiToKst0Kst1JPCLS.end(); ++itJPCLS){ + (*itJPCLS)->print(os); + os << "\n" << std::endl; + } + + os << "*** Chi_c0: LS combinations for the decay to K*1 K*2 *** "<< std::endl; + for ( itJPCLS=_ChiToKst1Kst2JPCLS.begin(); itJPCLS!=_ChiToKst1Kst2JPCLS.end(); ++itJPCLS){ + (*itJPCLS)->print(os); + os << "\n" << std::endl; + } + os << "*** Chi_c0: LS combinations for the decay to f0 f0 *** "<< std::endl; for ( itJPCLS=_ChiTof0f0JPCLS.begin(); itJPCLS!=_ChiTof0f0JPCLS.end(); ++itJPCLS){ (*itJPCLS)->print(os); @@ -277,6 +298,13 @@ void Psi2STo2K2PiGamStates::print(std::ostream& os) const os << "\n" << std::endl; } + os << "\n *** LS combinations for the decay K2 to f0 K *** "<< std::endl; + for ( itJPCLS=_K2Tof0KJPCLS.begin(); itJPCLS!=_K2Tof0KJPCLS.end(); ++itJPCLS){ + (*itJPCLS)->print(os); + os << "\n" << std::endl; + } + + os << "\n *** K+- *** "<< std::endl; _kJPC->print(os); diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh index b85172d2ae17b33f4d83c0012b66c871ae0a309c..e488f9615585d0ab3d64a4d5a7af5d16c9c6cb62 100644 --- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh +++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh @@ -35,6 +35,8 @@ public: std::vector< boost::shared_ptr<const JPCLS> > ChiTof0f2States() const {return _ChiTof0f2JPCLS;} std::vector< boost::shared_ptr<const JPCLS> > K1400ToKst1PiStates() const {return _K1400ToKst1PiJPCLS;} std::vector< boost::shared_ptr<const JPCLS> > ChiToK0K0States() const {return _ChiToK0K0JPCLS;} + std::vector< boost::shared_ptr<const JPCLS> > ChiToKst0Kst1States() const {return _ChiToKst0Kst1JPCLS;} + std::vector< boost::shared_ptr<const JPCLS> > ChiToKst1Kst2States() const {return _ChiToKst1Kst2JPCLS;} std::vector< boost::shared_ptr<const JPCLS> > K1ToK0PiStates() const {return _K1ToK0PiJPCLS;} std::vector< boost::shared_ptr<const JPCLS> > ChiToPi_2PiStates() const {return _ChiToPi_2PiJPCLS;} std::vector< boost::shared_ptr<const JPCLS> > Pi_2Tof_2PiStates() const {return _Pi_2Tof_2PiJPCLS;} @@ -47,6 +49,7 @@ public: std::vector< boost::shared_ptr<const JPCLS> > K2mToK2pPiStates() const {return _K2mToK2pPiJPCLS;} std::vector< boost::shared_ptr<const JPCLS> > K1pTof0KStates() const {return _K1pTof0KJPCLS;} std::vector< boost::shared_ptr<const JPCLS> > Pi0pToKstKStates() const {return _Pi0pToKstKJPCLS;} + std::vector< boost::shared_ptr<const JPCLS> > K2Tof0KStates() const {return _K2Tof0KJPCLS;} void print(std::ostream& os) const; protected: @@ -77,6 +80,8 @@ private: std::vector< boost::shared_ptr<const JPCLS> > _ChiTof0f0JPCLS; std::vector< boost::shared_ptr<const JPCLS> > _ChiTof0f2JPCLS; std::vector< boost::shared_ptr<const JPCLS> > _ChiToK0K0JPCLS; + std::vector< boost::shared_ptr<const JPCLS> > _ChiToKst0Kst1JPCLS; + std::vector< boost::shared_ptr<const JPCLS> > _ChiToKst1Kst2JPCLS; std::vector< boost::shared_ptr<const JPCLS> > _K1ToK0PiJPCLS; std::vector< boost::shared_ptr<const JPCLS> > _K1400ToKst1PiJPCLS; std::vector< boost::shared_ptr<const JPCLS> > _ChiToPi_2PiJPCLS; @@ -93,6 +98,7 @@ private: std::vector< boost::shared_ptr<const JPCLS> > _Kst2JPCLS; std::vector< boost::shared_ptr<const JPCLS> > _f2JPCLS; std::vector< boost::shared_ptr<const JPCLS> > _Pi0pToKstKJPCLS; + std::vector< boost::shared_ptr<const JPCLS> > _K2Tof0KJPCLS; void fillJPCLS(boost::shared_ptr<jpcRes>, boost::shared_ptr<jpcRes>, boost::shared_ptr<jpcRes>, std::vector< boost::shared_ptr<const JPCLS> >& ); }; diff --git a/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc b/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc deleted file mode 100644 index dadecd82984b36e90d387f95718aa2bc48a31a80..0000000000000000000000000000000000000000 --- a/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc +++ /dev/null @@ -1,767 +0,0 @@ -#include <vector> -#include <map> -#include <string.h> -#include <math.h> -#include <stdlib.h> -#include "TFile.h" -#include "TH1F.h" -#include "TNtuple.h" -#include "TMath.h" -#include <iostream> -#include "TCanvas.h" -#include "TF1.h" - -TFile* fdata; -TNtuple* ntdata; -TFile* fpwamc; -TNtuple* ntpwamc; - -TH1F* dataCosThetaKKpipi; -TH1F* pwamcCosThetaKKpipi; -TH1F* dataPhiKKpipi; -TH1F* pwamcPhiKKpipi; -TH1F* dataCosThetaKpipi; -TH1F* pwamcCosThetaKpipi; -TH1F* dataPhiKpipi; -TH1F* pwamcPhiKpipi; -TH1F* dataCosThetapipi; -TH1F* pwamcCosThetapipi; -TH1F* dataPhipipi; -TH1F* pwamcPhipipi; -TH1F* dataCosThetapi; -TH1F* pwamcCosThetapi; -TH1F* dataPhipi; -TH1F* pwamcPhipi; -TH1F* datainvMassKpipi; -TH1F* pwamcinvMassKpipi; -TH1F* datainvMasspipi; -TH1F* pwamcinvMasspipi; - -TH1F* histofpull; - -Int_t evtNo = 0; -float datacosThetaKKpipi; -float dataphiKKpipi; -float datacosThetaK1pi1pi2; -float dataphiK1pi1pi2; -float datacosThetapi1pi2; -float dataphipi1pi2; -float datacosThetapi1; -float dataphipi1; -float datainvmassK1pi1pi2; -float datainvmasspi1pi2; - -float pwamccosThetaKKpipi; -float pwamcphiKKpipi; -float pwamccosThetaK1pi1pi2; -float pwamcphiK1pi1pi2; -float pwamccosThetapi1pi2; -float pwamcphipi1pi2; -float pwamccosThetapi1; -float pwamcphipi1; -float pwamcinvmassK1pi1pi2; -float pwamcinvmasspi1pi2; - -float evtweight; - -std::vector<TH1F*> histVectData; -std::vector<TH1F*> histVectPwaMc; - -Int_t numOfEntriesData; -Int_t numOfEntriesPwaMc; - -TString type = "mc"; -bool drawFlag = true; -bool output = false; -bool longoutput = false; -bool silent = false; -TString method = "both"; -using namespace std; - -struct nextneighbours -{ - nextneighbours( unsigned int theevtnumber1, unsigned int theevtnumber2, float thedistance, bool thedatatype){ - evtnumber1 = theevtnumber1; - evtnumber2 = theevtnumber2; - distance = thedistance; - datatype = thedatatype; - } - - virtual ~nextneighbours(){}; - - unsigned int evtnumber1; - unsigned int evtnumber2; - float distance; - bool datatype; -}; - -void fillHistograms(); -void drawHistograms(); -void chisq_method(); -void mixed_sample(); -float GetDistanceDataToData(int dataentry1, int dataentry2); -float GetDistanceDataToPwaMc(int dataentry, int pwamcentry); -float GetDistancePwaMcToPwaMc(int pwamcentry1, int pwamcentry2); -bool Compare(const nextneighbours a, const nextneighbours b); - -/* -void initNtData(){ - - if (0 != fdata){ - delete fdata; - fdata=0; - } - - TString fileData="./Psi2STo2K2PiGamPWA_data.root"; - fdata = new TFile(fileData,"READ"); - ntdata = (TNtuple*)fdata->Get("dataTuple"); - ntdata->SetBranchAddress("costhetaKKpipi", &datacosThetaKKpipi); - ntdata->SetBranchAddress("phiKKpipi", &dataphiKKpipi); - ntdata->SetBranchAddress("costhetaK1pi1pi2", &datacosThetaK1pi1pi2); - ntdata->SetBranchAddress("phiK1pi1pi2", &dataphiK1pi1pi2); - ntdata->SetBranchAddress("costhetapipiViaK1pipi", &datacosThetapi1pi2); - ntdata->SetBranchAddress("phipipiViaK1pipi", &dataphipi1pi2); - ntdata->SetBranchAddress("costhetapiViapipi", &datacosThetapi1); - ntdata->SetBranchAddress("phipiViapipi", &dataphipi1); - ntdata->SetBranchAddress("mk1pi1pi2", &datainvmassK1pi1pi2); - ntdata->SetBranchAddress("mpipi", &datainvmasspi1pi2); -} - - -void initNtPwamc(){ - - if (0 != fpwamc){ - delete fpwamc; - fpwamc=0; - } - - TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc__FitParamHyp8_base_K0K2_newflatte__10000events.root"; - fpwamc = new TFile(filePwaMc,"READ"); - ntpwamc = (TNtuple*)fpwamc->Get("dataTuple"); - ntpwamc->SetBranchAddress("costhetaKKpipi", &pwamccosThetaKKpipi); - ntpwamc->SetBranchAddress("phiKKpipi", &pwamcphiKKpipi); - ntpwamc->SetBranchAddress("costhetaK1pi1pi2", &pwamccosThetaK1pi1pi2); - ntpwamc->SetBranchAddress("phiK1pi1pi2", &pwamcphiK1pi1pi2); - ntpwamc->SetBranchAddress("costhetapipiViaK1pipi", &pwamccosThetapi1pi2); - ntpwamc->SetBranchAddress("phipipiViaK1pipi", &pwamcphipi1pi2); - ntpwamc->SetBranchAddress("costhetapiViapipi", &pwamccosThetapi1); - ntpwamc->SetBranchAddress("phipiViapipi", &pwamcphipi1); - ntpwamc->SetBranchAddress("mk1pi1pi2", &pwamcinvmassK1pi1pi2); - ntpwamc->SetBranchAddress("mpipi", &pwamcinvmasspi1pi2); -} -*/ - -int main(){ - - TString fileData="./Psi2STo2K2PiGamPWA_data.root"; - //TString fileData="./Psi2STo2K2PiGamPWA_pwamc__FitParamHyp8_base_K0K2_newflatte__10000events.root"; - TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc__FitParamHyp8_base_K0K2_newflatte__10000events.root"; - TString option="mixedsample,pwamc,nodraw,silent"; - - cout << "--- ENTERING QAMIXEDSAMPLE" << endl; - // gROOT->SetStyle("Plain"); - // gStyle->SetCanvasColor(0); - // gStyle->SetStatBorderSize(0); - // gStyle->SetPalette(1); - // gStyle->SetOptStat(1111); - - - if (option.Contains("pwamc")) {type = "pwamc";} - if (option.Contains("nodraw")) {drawFlag = false;} - if (option.Contains("chisq")) {method = "chisq";} - if (option.Contains("mixedsample")) {method = "mixedsample";} - if (option.Contains("output")) {output = true;} - if (option.Contains("longoutput")) {longoutput = true;} - if (option.Contains("silent")) {silent = true;} - - fdata = new TFile(fileData,"READ"); - ntdata = (TNtuple*)fdata->Get("dataTuple"); - - fpwamc = new TFile(filePwaMc,"READ"); - if(type=="pwamc") { - cout << "Vergleich mit PWA-MC" << endl; - ntpwamc = (TNtuple*)fpwamc->Get("dataTuple"); - } - else if(type=="mc") { - cout << "Vergleich mit MC" << endl; - ntpwamc = (TNtuple*)fdata->Get("mcTuple"); - } - cout << endl; - - ntdata->SetBranchAddress("costhetaKKpipi", &datacosThetaKKpipi); - ntdata->SetBranchAddress("phiKKpipi", &dataphiKKpipi); - ntdata->SetBranchAddress("costhetaK1pi1pi2", &datacosThetaK1pi1pi2); - ntdata->SetBranchAddress("phiK1pi1pi2", &dataphiK1pi1pi2); - ntdata->SetBranchAddress("costhetapipiViaK1pipi", &datacosThetapi1pi2); - ntdata->SetBranchAddress("phipipiViaK1pipi", &dataphipi1pi2); - ntdata->SetBranchAddress("costhetapiViapipi", &datacosThetapi1); - ntdata->SetBranchAddress("phipiViapipi", &dataphipi1); - ntdata->SetBranchAddress("mk1pi1pi2", &datainvmassK1pi1pi2); - ntdata->SetBranchAddress("mpipi", &datainvmasspi1pi2); - - - ntpwamc->SetBranchAddress("costhetaKKpipi", &pwamccosThetaKKpipi); - ntpwamc->SetBranchAddress("phiKKpipi", &pwamcphiKKpipi); - ntpwamc->SetBranchAddress("costhetaK1pi1pi2", &pwamccosThetaK1pi1pi2); - ntpwamc->SetBranchAddress("phiK1pi1pi2", &pwamcphiK1pi1pi2); - ntpwamc->SetBranchAddress("costhetapipiViaK1pipi", &pwamccosThetapi1pi2); - ntpwamc->SetBranchAddress("phipipiViaK1pipi", &pwamcphipi1pi2); - ntpwamc->SetBranchAddress("costhetapiViapipi", &pwamccosThetapi1); - ntpwamc->SetBranchAddress("phipiViapipi", &pwamcphipi1); - ntpwamc->SetBranchAddress("mk1pi1pi2", &pwamcinvmassK1pi1pi2); - ntpwamc->SetBranchAddress("mpipi", &pwamcinvmasspi1pi2); - ntpwamc->SetBranchAddress("weight",&evtweight); - - numOfEntriesData = ntdata->GetEntries(); - numOfEntriesPwaMc = ntpwamc->GetEntries(); - - - dataCosThetaKKpipi = new TH1F("dataCosThetaKKpipi","dataCosThetaKKpipi",60,-1.,1.); - dataPhiKKpipi = new TH1F("dataPhiKKpipi","dataPhiKKpipi",60,-TMath::Pi(),TMath::Pi()); - dataCosThetaKpipi = new TH1F("dataCosThetaKpipi","dataCosThetaKpipi",60,-1.,1.); - dataPhiKpipi = new TH1F("dataPhiKpipi","dataPhiKpipi",60,-TMath::Pi(),TMath::Pi()); - dataCosThetapipi = new TH1F("dataCosThetapipi","dataCosThetapipi",60,-1.,1.); - dataPhipipi = new TH1F("dataPhipipi","dataPhipipi",60,-TMath::Pi(),TMath::Pi()); - dataCosThetapi = new TH1F("dataCosThetapi","dataCosThetapi",60,0.,1.); - dataPhipi = new TH1F("dataPhipi","dataPhipi",60,0.,TMath::Pi()); - datainvMassKpipi = new TH1F("datainvMassKpipi", "datainvMassKpipi", 60, 0.7, 3.0); - datainvMasspipi = new TH1F("datainvMasspipi", "datainvMasspipi", 60, 0.2, 2.8); - - pwamcCosThetaKKpipi = new TH1F("pwamcCosThetaKKpipi","pwamcCosThetaKKpipi",60,-1.,1.); - pwamcPhiKKpipi = new TH1F("pwamcPhiKKpipi","pwamcPhiKKpipi",60,-TMath::Pi(),TMath::Pi()); - pwamcCosThetaKpipi = new TH1F("pwamcCosThetaKpipi","pwamcCosThetaKpipi",60,-1.,1.); - pwamcPhiKpipi = new TH1F("pwamcPhiKpipi","pwamcPhiKpipi",60,-TMath::Pi(),TMath::Pi()); - pwamcCosThetapipi = new TH1F("pwamcCosThetapipi","pwamcCosThetapipi",60,-1.,1.); - pwamcPhipipi = new TH1F("pwamcPhipipi","pwamcPhipipi",60,-TMath::Pi(),TMath::Pi()); - pwamcCosThetapi = new TH1F("pwamcCosThetapi","pwamcCosThetapi",60,0.,1.); - pwamcPhipi = new TH1F("pwamcPhipi","pwamcPhipi",60,0.,TMath::Pi()); - pwamcinvMassKpipi = new TH1F("pwamcinvMassKpipi", "pwamcinvMassKpipi", 60, 0.7, 3.0); - pwamcinvMasspipi = new TH1F("pwamcinvMasspipi", "pwamcinvMasspipi", 60, 0.2, 2.8); - - histofpull = new TH1F("histofpull", "histofpull", 100, -10, 10); - - histVectData.push_back(dataCosThetaKKpipi); - histVectData.push_back(dataPhiKKpipi); - histVectData.push_back(dataCosThetaKpipi); - histVectData.push_back(dataPhiKpipi); - histVectData.push_back(dataCosThetapipi); - histVectData.push_back(dataPhipipi); - histVectData.push_back(dataCosThetapi); - histVectData.push_back(dataPhipi); - histVectData.push_back(datainvMassKpipi); - histVectData.push_back(datainvMasspipi); - - histVectPwaMc.push_back(pwamcCosThetaKKpipi); - histVectPwaMc.push_back(pwamcPhiKKpipi); - histVectPwaMc.push_back(pwamcCosThetaKpipi); - histVectPwaMc.push_back(pwamcPhiKpipi); - histVectPwaMc.push_back(pwamcCosThetapipi); - histVectPwaMc.push_back(pwamcPhipipi); - histVectPwaMc.push_back(pwamcCosThetapi); - histVectPwaMc.push_back(pwamcPhipi); - histVectPwaMc.push_back(pwamcinvMassKpipi); - histVectPwaMc.push_back(pwamcinvMasspipi); - - fillHistograms(); - - if(drawFlag) { - cout << "--- Histogramme gezeichnet." << endl; - drawHistograms(); - cout << endl; - cout << endl; - cout << "-----------------------------------------" << endl; - } - if(method=="chisq" || method=="both") { - cout << "--- CHISQ-METHOD-ERGEBNISSE:" << endl; - cout << endl; - chisq_method(); - cout << "-----------------------------------------" << endl; - cout << endl; - } - if(method=="mixedsample" || method=="both") { - cout << "--- MIXED-SAMPLE-ERGEBNISSE:" << endl; - cout << endl; - mixed_sample(); - cout << "-----------------------------------------" << endl; - cout << endl; - } - cout << "--- EXITING QAMIXEDSAMPLE" << endl; - - return 1; - -} - - -bool Compare(const nextneighbours a, const nextneighbours b) { - return a.distance < b.distance; -} - - -void mixed_sample() { - - if(!silent) cout << "Entering mixed-sample method" << endl; - if(!silent) cout << endl; - - unsigned int limit_runs_lower = 0; - unsigned int limit_runs_upper = 10; - cout << "Performing " << limit_runs_upper - limit_runs_lower << " runs." << endl; - std::vector<float> fpullvector; - unsigned int limit_data = 1000; // numOfEntriesData; - unsigned int limit_pwamc = 1000; // numOfEntriesPwaMc; - unsigned int numberofneighbours = 10; - - - for(unsigned int run = limit_runs_lower; run < limit_runs_upper; ++run) { - - std::vector<struct nextneighbours> neighbourList; - cout << "Starting run " << run+1 << endl; - - if(!silent) cout << endl; - - float distance = 1E6; - float T=0; - // unsigned int limit_data = 2500; // numOfEntriesData; - // unsigned int limit_pwamc = 5000; // numOfEntriesPwaMc; - // unsigned int numberofneighbours = 10; - unsigned int pwaoffset = 100000; - - if(!silent) cout << "--- Calculating distances DATA <-> DATA" << endl; - for(unsigned int i = 0; i<limit_data; i++) { - for(unsigned int j = i; j<limit_data; j++) { - //cout << run << ", loop 1:\t" << i << "\t" << j << endl; - if(i!=j) { - distance = GetDistanceDataToData(i,j); - nextneighbours tmpNeighbours(i,j,distance,true); - neighbourList.push_back(tmpNeighbours); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; - } - } - } - - if(!silent) cout << "--- Calculating distances DATA <-> PWAMC" << endl; - for(unsigned int i = 0; i<limit_data; i++) { - for(unsigned int j = 0; j<limit_pwamc; j++) { - //cout << run << ", loop 2:\t" << i << "\t" << j << endl; - distance = GetDistanceDataToPwaMc(i,j+run*limit_pwamc); - nextneighbours tmpNeighbours(i,j+pwaoffset,distance,false); - neighbourList.push_back(tmpNeighbours); - // neighbourList.push_back((struct nextneighbours) {i,j+pwaoffset,distance,false}); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; - } - } - - if(!silent) cout << "--- Calculating distances PWAMC <-> PWAMC" << endl; - for(unsigned int i = 0; i<limit_pwamc; i++) { - for(unsigned int j = i; j<limit_pwamc; j++) { - //cout << run << ", loop 3:\t" << i << "\t" << j << endl; - if(i!=j) { - distance = GetDistancePwaMcToPwaMc(i+run*limit_pwamc,j+run*limit_pwamc); - nextneighbours tmpNeighbours(i+pwaoffset,j+pwaoffset,distance,true); - neighbourList.push_back(tmpNeighbours); - // neighbourList.push_back((struct nextneighbours) {i+pwaoffset,j+pwaoffset,distance,true}); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; - } - } - } - - if(!silent) { - cout << endl; - cout << "How many distances have been calculated?\t" << neighbourList.size() << endl; - cout << endl; - cout << endl; - } - - if(longoutput) { - cout << "--- Distances:" << endl; - for(unsigned int i = 0; i < neighbourList.size(); i++) { - cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; - } - } - if(!silent) cout << endl; - - std::sort(neighbourList.begin(),neighbourList.end(), Compare); - if(!silent) cout << "neighbourList has been sorted!" << endl; - if(!silent) cout << endl; - - if(longoutput) { - cout << "--- Distances after sorting:" << endl; - for(unsigned int i = 0; i < neighbourList.size(); i++) { - cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; - } - } - if(!silent) cout << endl; - - float numOfEntriesTotal = float(limit_data + limit_pwamc); - - unsigned int m; - - for(unsigned int i = 0; i < limit_data; i++) { - if(output) cout << "--- Checking neighbours of event " << i << endl; - m = 0; - for(unsigned int j = 0; j < neighbourList.size(); j++) - { - if (m >= numberofneighbours) continue; - if(longoutput) cout << "Checking Event in neighbourList for data: " << j << endl; - if(neighbourList[j].evtnumber1 == i || neighbourList[j].evtnumber2 == i) - { - m++; - if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i; - if(neighbourList[j].datatype==true) - { - if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; - T=T+1.; - } - else - { - if(longoutput) cout << ", but doesn't belong to the same sample." << endl; - } - } - else - { - if(longoutput) cout << "Not a neighbour of event " << i << endl; - } - } - if(output) cout << "Event " << i << " has been analysed with " << m << " neighbours." << endl; - if(output) cout << "T: " << T << endl; - if(output) cout << endl; - } - - for(unsigned int i = 0; i < limit_pwamc; i++) { - if(output) cout << "--- Checking neighbours of event " << i+pwaoffset << endl; - m = 0; - - for(unsigned int j = 0; j < neighbourList.size(); j++) - { - if (m >= numberofneighbours) continue; - if(longoutput) cout << "Checking Event in neighbourList for pwamc: " << j << endl; - if(neighbourList[j].evtnumber1 == i+pwaoffset || neighbourList[j].evtnumber2 == i+pwaoffset) - { - m++; - if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i+pwaoffset; - if(neighbourList[j].datatype==true) - { - if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; - T=T+1; - } - else - { - if(longoutput) cout << ", but doesn't belong to the same sample." << endl; - } - } - else - { - if(longoutput) cout << "Not a neighbour of event " << i+pwaoffset << endl; - } - - } - - if(output) cout << "Event " << i+pwaoffset << " has been analysed." << endl; - if(output) cout << "T: " << T << endl; - if(output) cout << endl; - } - - - if(!silent) { - cout << endl; - cout << "Number of Entries in DATA: " << numOfEntriesData << " events available (used " << limit_data << ")"<< endl; - cout << "Number of Entries in PWAMC: " << numOfEntriesPwaMc << " events available (used " << limit_pwamc << ")"<< endl; - cout << "Total Number of Entries: " << numOfEntriesTotal << endl; - cout << endl; - } - - numOfEntriesData = limit_data; - numOfEntriesPwaMc = limit_pwamc; - - float Tnorm = 0; - Tnorm = T/(float(numberofneighbours)*float(numOfEntriesTotal)); - if(!silent) cout << "Tnorm: " << Tnorm << endl; - - - float mu = ((float(numOfEntriesData)*(float(numOfEntriesData)-1)) + (float(numOfEntriesPwaMc)*(float(numOfEntriesPwaMc)-1))) / (float(numOfEntriesTotal)*(float(numOfEntriesTotal)-1)); - if(!silent) cout << "Erwartungswert mu: " << mu << endl; - - - float variancesquared = - 1 / (float(numOfEntriesTotal)*float(numberofneighbours)) - * - ( - (float(numOfEntriesData)*float(numOfEntriesPwaMc))/(pow(float(numOfEntriesTotal), 2)) - + - 4*(pow(float(numOfEntriesData), 2)*pow(float(numOfEntriesPwaMc), 2))/(pow(float(numOfEntriesTotal),4)) - ); - - float fpull = (Tnorm-mu)/sqrt(variancesquared); - - if(!silent) { - cout << "Varianz: " << variancesquared << endl; - cout << endl; - cout << "Differenz zwischen Resultat und Erwartungswert: " << Tnorm-mu << endl; - } - cout << "f_pull in run " << run+1 << ": " << fpull << endl; - cout << endl; - fpullvector.push_back(fpull); - } - - cout << "--------------------------" << endl; - cout << fpullvector.size() << endl; - for(unsigned int i = 0; i < fpullvector.size(); i++) { - cout << "Calculated fpull:\t" << fpullvector[i] << endl; - histofpull->Fill( fpullvector[i] ); - } - - - TCanvas* chist1 = new TCanvas("chist1","chist1",1200,800); - histofpull->SetLineColor(2); - histofpull->Draw(); - - TF1* fit1 = new TF1("fit1","gaus(0)",-10.,10.); - fit1->SetParameters(5,histofpull->GetMean(),histofpull->GetRMS()); - fit1->SetLineColor(4); - fit1->SetLineWidth(3); - histofpull->Fit("fit1","LMRS","",-10.,10.); - chist1->Print("histofpull.pdf"); - -} - - - - -float GetDistanceDataToData(int dataentry1, int dataentry2) { - // cout << "Distance from Data Event " << dataentry1 << " to Data Event " << dataentry2 << endl; - - ntdata->GetEntry(dataentry1); - float datacosThetaKKpipi_1 = datacosThetaKKpipi; - float dataphiKKpipi_1 = dataphiKKpipi; - float datacosThetaK1pi1pi2_1 = datacosThetaK1pi1pi2; - float dataphiK1pi1pi2_1 = dataphiK1pi1pi2; - float datacosThetapi1pi2_1 = datacosThetapi1pi2; - float dataphipi1pi2_1 = dataphipi1pi2; - float datacosThetapi1_1 = datacosThetapi1; - float dataphipi1_1 = dataphipi1; - float datainvmassK1pi1pi2_1 = datainvmassK1pi1pi2; - float datainvmasspi1pi2_1 = datainvmasspi1pi2; - ntdata->GetEntry(dataentry2); - float datacosThetaKKpipi_2 = datacosThetaKKpipi; - float dataphiKKpipi_2 = dataphiKKpipi; - float datacosThetaK1pi1pi2_2 = datacosThetaK1pi1pi2; - float dataphiK1pi1pi2_2 = dataphiK1pi1pi2; - float datacosThetapi1pi2_2 = datacosThetapi1pi2; - float dataphipi1pi2_2 = dataphipi1pi2; - float datacosThetapi1_2 = datacosThetapi1; - float dataphipi1_2 = dataphipi1; - float datainvmassK1pi1pi2_2 = datainvmassK1pi1pi2; - float datainvmasspi1pi2_2 = datainvmasspi1pi2; - - float distance = sqrt( - pow((datacosThetaKKpipi_1-datacosThetaKKpipi_2)/2.,2) + - pow((dataphiKKpipi_1-dataphiKKpipi_2)/(2*3.1415),2) + - pow((datacosThetaK1pi1pi2_1-datacosThetaK1pi1pi2_2)/2. ,2) + - pow((dataphiK1pi1pi2_1-dataphiK1pi1pi2_2)/(2*3.1415) ,2) + - pow((datacosThetapi1pi2_1-datacosThetapi1pi2_2)/2. ,2) + - pow((dataphipi1pi2_1-dataphipi1pi2_2)/(2*3.1415) ,2) + - pow((datacosThetapi1_1-datacosThetapi1_2) ,2) + - pow((dataphipi1_1-dataphipi1_2)/(3.1415) ,2) + - pow((datainvmassK1pi1pi2_1-datainvmassK1pi1pi2_2 )/2.2 ,2) + - pow((datainvmasspi1pi2_1-datainvmasspi1pi2_2)/2.3 ,2) - ); - - return distance; - -} - - - -float GetDistanceDataToPwaMc(int dataentry, int pwamcentry) { - // cout << "Distance from Data Event " << dataentry << " to PwaMc Event " << pwamcentry << endl; - - ntdata->GetEntry(dataentry); - ntpwamc->GetEntry(pwamcentry); - - float distance = sqrt( - pow((datacosThetaKKpipi-pwamccosThetaKKpipi)/2.,2) + - pow((dataphiKKpipi-pwamcphiKKpipi)/(2*3.1415),2) + - pow((datacosThetaK1pi1pi2-pwamccosThetaK1pi1pi2)/2. ,2) + - pow((dataphiK1pi1pi2-pwamcphiK1pi1pi2)/(2*3.1415) ,2) + - pow((datacosThetapi1pi2-pwamccosThetapi1pi2)/2. ,2) + - pow((dataphipi1pi2-pwamcphipi1pi2)/(2*3.1415) ,2) + - pow((datacosThetapi1-pwamccosThetapi1) ,2.) + - pow((dataphipi1-pwamcphipi1)/(3.1415) ,2) + - pow((datainvmassK1pi1pi2-pwamcinvmassK1pi1pi2)/2.2 ,2) + - pow((datainvmasspi1pi2-pwamcinvmasspi1pi2)/2.3 ,2) - ); - - return distance; - -} - - - - - - - - -float GetDistancePwaMcToPwaMc(int pwamcentry1, int pwamcentry2) { - // cout << "Distance from PwaMC Event " << dataentry1 << " to PwaMC Event " << dataentry2 << endl; - - ntpwamc->GetEntry(pwamcentry1); - float pwamccosThetaKKpipi_1 = pwamccosThetaKKpipi; - float pwamcphiKKpipi_1 = pwamcphiKKpipi; - float pwamccosThetaK1pi1pi2_1 = pwamccosThetaK1pi1pi2; - float pwamcphiK1pi1pi2_1 = pwamcphiK1pi1pi2; - float pwamccosThetapi1pi2_1 = pwamccosThetapi1pi2; - float pwamcphipi1pi2_1 = pwamcphipi1pi2; - float pwamccosThetapi1_1 = pwamccosThetapi1; - float pwamcphipi1_1 = pwamcphipi1; - float pwamcinvmassK1pi1pi2_1 = pwamcinvmassK1pi1pi2; - float pwamcinvmasspi1pi2_1 = pwamcinvmasspi1pi2; - ntpwamc->GetEntry(pwamcentry2); - float pwamccosThetaKKpipi_2 = pwamccosThetaKKpipi; - float pwamcphiKKpipi_2 = pwamcphiKKpipi; - float pwamccosThetaK1pi1pi2_2 = pwamccosThetaK1pi1pi2; - float pwamcphiK1pi1pi2_2 = pwamcphiK1pi1pi2; - float pwamccosThetapi1pi2_2 = pwamccosThetapi1pi2; - float pwamcphipi1pi2_2 = pwamcphipi1pi2; - float pwamccosThetapi1_2 = pwamccosThetapi1; - float pwamcphipi1_2 = pwamcphipi1; - float pwamcinvmassK1pi1pi2_2 = pwamcinvmassK1pi1pi2; - float pwamcinvmasspi1pi2_2 = pwamcinvmasspi1pi2; - - float distance = sqrt( - pow((pwamccosThetaKKpipi_1-pwamccosThetaKKpipi_2)/2.,2) + - pow((pwamcphiKKpipi_1-pwamcphiKKpipi_2)/(2*3.1415),2) + - pow((pwamccosThetaK1pi1pi2_1-pwamccosThetaK1pi1pi2_2)/2. ,2) + - pow((pwamcphiK1pi1pi2_1-pwamcphiK1pi1pi2_2)/(2*3.1415) ,2) + - pow((pwamccosThetapi1pi2_1-pwamccosThetapi1pi2_2)/2. ,2) + - pow((pwamcphipi1pi2_1-pwamcphipi1pi2_2)/(2*3.1415) ,2) + - pow((pwamccosThetapi1_1-pwamccosThetapi1_2) ,2) + - pow((pwamcphipi1_1-pwamcphipi1_2)/(3.1415) ,2) + - pow((pwamcinvmassK1pi1pi2_1-pwamcinvmassK1pi1pi2_2)/2.2 ,2) + - pow((pwamcinvmasspi1pi2_1-pwamcinvmasspi1pi2_2)/2.3 ,2) - ); - - return distance; - -} - - - - - - - -void chisq_method() { - - cout << "Entering chisq-method" << endl; - - Double_t chisq = 0; - - Int_t constraints = 50; // 12-01-06_FitParamHyp8_base_K0K2_Hyp5/rerun_correctedampK0K2/Psi2STo2K2PiGamPWA_noK1_1680Hyp7.root - // float constraints = 94; // 12-01-06_FitParamHyp8_base_K0K2_Hyp5/rerun_correctedampK0K2/Psi2STo2K2PiGamPWA_noK1_1680Hyp7_noK1_2300Hyp7.root - Int_t ndf; - int numberofbins = 0; - - for(unsigned int i=0; i<histVectData.size(); i++) { - int currentNoOfBins=0; - - for (int j=0; j<(histVectData[i]->GetSize()); j++) { - float binContentData=float(histVectData[i]->GetBinContent(j)); - if (binContentData>0.){ - chisq += pow( binContentData - (histVectPwaMc[i]->GetBinContent(j)),2)/binContentData; - currentNoOfBins++; - } - // cout << "chi^2 after bin " << j << " of histogram " << histVectData[i]->GetName() << ": " << chisq << endl; - } - cout << "chi^2 after " << histVectData[i]->GetName() << ": " << chisq << endl; - - numberofbins += currentNoOfBins; - } - - ndf = numberofbins - constraints; - - cout << "Total chi^2:\t" << chisq << endl; - cout << "Total number of bins:\t" << numberofbins << endl; - cout << "Number of constraints:\t" << constraints << endl; - cout << "Degrees of freedom:\t" << ndf << endl; - cout << "chi^2/ndf:\t" << chisq/double(ndf) << endl; - cout << "CL:\t" << TMath::Prob(chisq, ndf) << endl; -} - - - - - - -void fillHistograms() { - - Int_t maxEvtNo = 0; - for(Int_t i=0; i<max(numOfEntriesPwaMc,numOfEntriesData); i++) - { - if(i<numOfEntriesData){ - ntdata->GetEntry(i); - - dataCosThetaKKpipi->Fill(datacosThetaKKpipi); - dataPhiKKpipi->Fill(dataphiKKpipi); - - dataCosThetaKpipi->Fill(datacosThetaK1pi1pi2); - dataPhiKpipi->Fill(dataphiK1pi1pi2); - - dataCosThetapipi->Fill(datacosThetapi1pi2); - dataPhipipi->Fill(dataphipi1pi2); - - dataCosThetapi->Fill(fabs(datacosThetapi1)); - dataPhipi->Fill(fabs(dataphipi1)); - - datainvMassKpipi->Fill(datainvmassK1pi1pi2); - datainvMasspipi->Fill(datainvmasspi1pi2); - } - - ntpwamc->GetEntry(i); - - pwamcCosThetaKKpipi->Fill(pwamccosThetaKKpipi, evtweight); - pwamcPhiKKpipi->Fill(pwamcphiKKpipi, evtweight); - - pwamcCosThetaKpipi->Fill(pwamccosThetaK1pi1pi2, evtweight); - pwamcPhiKpipi->Fill(pwamcphiK1pi1pi2, evtweight); - - pwamcCosThetapipi->Fill(pwamccosThetapi1pi2, evtweight); - pwamcPhipipi->Fill(pwamcphipi1pi2, evtweight); - - pwamcCosThetapi->Fill(fabs(pwamccosThetapi1), evtweight); - pwamcPhipi->Fill(fabs(pwamcphipi1), evtweight); - - pwamcinvMassKpipi->Fill(pwamcinvmassK1pi1pi2, evtweight); - pwamcinvMasspipi->Fill(pwamcinvmasspi1pi2, evtweight); - - if(maxEvtNo<evtNo) maxEvtNo=evtNo; - } - - -} - - - - -void drawHistograms() { - - double scalefactor = double(ntdata->GetEntries())/double(ntpwamc->GetEntries()); - - Int_t rebin = 1; - - TCanvas* cmain = new TCanvas("cmain","cmain",1400,600); - cmain->Divide(5,2); - for(unsigned int i=0; i<histVectData.size(); i++) { - cmain->cd(i+1); - histVectData[i]->Rebin(rebin); - histVectData[i]->SetLineColor(2); - histVectData[i]->Draw("E"); - histVectData[i]->SetMinimum(0); - histVectPwaMc[i]->Rebin(rebin); - histVectPwaMc[i]->Scale(scalefactor); - histVectPwaMc[i]->Draw("same"); - } - -} - - - diff --git a/Examples/Psi2STo2K2PiGam/qaMixedSampleCacheApp.cc b/Examples/Psi2STo2K2PiGam/qaMixedSampleCacheApp.cc index 23a983ee2a0cd2b5386205b2b3924c6fc909d9de..9bddcaaced46dbc5d8a1a84c7d45cb65f59c628c 100644 --- a/Examples/Psi2STo2K2PiGam/qaMixedSampleCacheApp.cc +++ b/Examples/Psi2STo2K2PiGam/qaMixedSampleCacheApp.cc @@ -10,6 +10,8 @@ #include <iostream> #include "TCanvas.h" #include "TF1.h" +#include <sstream> + TFile* fdata; TNtuple* ntdata; @@ -134,10 +136,10 @@ bool Compare(const nextneighbours a, const nextneighbours b); int main(){ //TString fileData="./Psi2STo2K2PiGamPWA_data.root"; - //TString fileData="./Psi2STo2K2PiGamPWA_pwamc__FitParamHyp8_base_K0K2_newflatte__10000events.root"; - //TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc__FitParamHyp8_base_K0K2_newflatte__10000events.root"; + //TString fileData="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp1_76695events.root"; TString fileData="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp9_43791events.root"; - TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp9_43791events.root"; + TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp1_76695events.root"; + //TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp9_43791events.root"; TString option="mixedsample,pwamc,nodraw,silent"; cout << "--- ENTERING QAMIXEDSAMPLE" << endl; @@ -254,8 +256,6 @@ int main(){ pwamcinvMassKpipi = new TH1F("pwamcinvMassKpipi", "pwamcinvMassKpipi", 60, 0.7, 3.0); pwamcinvMasspipi = new TH1F("pwamcinvMasspipi", "pwamcinvMasspipi", 60, 0.2, 2.8); - histofpull = new TH1F("histofpull", "histofpull", 100, -10, 10); - histVectData.push_back(dataCosThetaKKpipi); histVectData.push_back(dataPhiKKpipi); histVectData.push_back(dataCosThetaKpipi); @@ -318,218 +318,265 @@ void mixed_sample() { if(!silent) cout << "Entering mixed-sample method" << endl; if(!silent) cout << endl; - unsigned int limit_runs_lower = 1; - unsigned int limit_runs_upper = 8; - cout << "Performing " << limit_runs_upper - limit_runs_lower << " runs." << endl; - std::vector<float> fpullvector; - unsigned int limit_data = 1000; // numOfEntriesData; - unsigned int limit_pwamc = 5000; // numOfEntriesPwaMc; - unsigned int numberofneighbours = 10; - - - for(unsigned int run = limit_runs_lower; run < limit_runs_upper; ++run) { - - std::vector<struct nextneighbours> neighbourList; - cout << "Starting run " << run+1 << endl; + std::vector<int> vec_limit_data; + vec_limit_data.push_back(100); + // vec_limit_data.push_back(200); + // vec_limit_data.push_back(500); + // vec_limit_data.push_back(1000); + // vec_limit_data.push_back(2000); + // vec_limit_data.push_back(3000); + // vec_limit_data.push_back(5000); + // vec_limit_data.push_back(7500); + // vec_limit_data.push_back(10000); + - if(!silent) cout << endl; + for(unsigned int i = 0; i<vec_limit_data.size(); i++) { + cout << "Run " << i << " will have " << vec_limit_data[i] << " data events." << endl; + } - float distance = 1E6; - float T=0; - unsigned int pwaoffset = 100000; + for(unsigned int k = 0; k<vec_limit_data.size(); k++) { + cout << "Cycle: " << k << endl; + //unsigned int limit_data = 2874; // numOfEntriesData (usually 2874); + unsigned int limit_data = vec_limit_data[k]; + unsigned int limit_pwamc = 10*limit_data; // numOfEntriesPwaMc; + unsigned int total_events_in_pwamc = 76695; + string sourcename1 = "Hyp9pwamc"; + string sourcename2 = "Hyp1pwamc"; + cout << "Using " << limit_data << " DATA events." << endl; + cout << "Using " << limit_pwamc << " PWAMC events." << endl; + + unsigned int limit_runs_lower = 0; + unsigned int limit_runs_upper = total_events_in_pwamc/limit_pwamc; + cout << "Performing " << limit_runs_upper - limit_runs_lower << " runs." << endl; + unsigned int numberofneighbours = 10; + + std::vector<float> fpullvector; + + for(unsigned int run = limit_runs_lower; run < limit_runs_upper; ++run) { - if(!silent) cout << "--- Calculating distances DATA <-> DATA" << endl; - for(unsigned int i = 0; i<limit_data; i++) { - for(unsigned int j = i; j<limit_data; j++) { - if(i!=j) { - distance = GetDistanceDataToData(i,j); - nextneighbours tmpNeighbours(i,j,distance,true); - neighbourList.push_back(tmpNeighbours); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; + std::vector<struct nextneighbours> neighbourList; + cout << "Starting run " << run+1 << endl; + + if(!silent) cout << endl; + + float distance = 1E6; + float T=0; + unsigned int pwaoffset = 100000; + + if(!silent) cout << "--- Calculating distances DATA <-> DATA" << endl; + for(unsigned int i = 0; i<limit_data; i++) { + for(unsigned int j = i; j<limit_data; j++) { + if(i!=j) { + distance = GetDistanceDataToData(i,j); + nextneighbours tmpNeighbours(i,j,distance,true); + neighbourList.push_back(tmpNeighbours); + if(neighbourList.size() % 500000 == 0 && !silent) cout << neighbourList.size() << "th distance has been calculated." << endl; + } } } - } - - if(!silent) cout << "--- Calculating distances DATA <-> PWAMC" << endl; - for(unsigned int i = 0; i<limit_data; i++) { - for(unsigned int j = 0; j<limit_pwamc; j++) { - distance = GetDistanceDataToPwaMc(i,j+run*limit_pwamc); - nextneighbours tmpNeighbours(i,j+pwaoffset,distance,false); - neighbourList.push_back(tmpNeighbours); - // neighbourList.push_back((struct nextneighbours) {i,j+pwaoffset,distance,false}); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; - } - } - - if(!silent) cout << "--- Calculating distances PWAMC <-> PWAMC" << endl; - for(unsigned int i = 0; i<limit_pwamc; i++) { - for(unsigned int j = i; j<limit_pwamc; j++) { - if(i!=j) { - distance = GetDistancePwaMcToPwaMc(i+run*limit_pwamc,j+run*limit_pwamc); - nextneighbours tmpNeighbours(i+pwaoffset,j+pwaoffset,distance,true); + + if(!silent) cout << "--- Calculating distances DATA <-> PWAMC" << endl; + for(unsigned int i = 0; i<limit_data; i++) { + for(unsigned int j = 0; j<limit_pwamc; j++) { + distance = GetDistanceDataToPwaMc(i,j+run*limit_pwamc); + nextneighbours tmpNeighbours(i,j+pwaoffset,distance,false); neighbourList.push_back(tmpNeighbours); - // neighbourList.push_back((struct nextneighbours) {i+pwaoffset,j+pwaoffset,distance,true}); - if(neighbourList.size() % 500000 == 0) cout << neighbourList.size() << "th distance has been calculated." << endl; + // neighbourList.push_back((struct nextneighbours) {i,j+pwaoffset,distance,false}); + if(neighbourList.size() % 500000 == 0 && !silent) cout << neighbourList.size() << "th distance has been calculated." << endl; } } - } - - if(!silent) { - cout << endl; - cout << "How many distances have been calculated?\t" << neighbourList.size() << endl; - cout << endl; - cout << endl; - } - - if(longoutput) { - cout << "--- Distances:" << endl; - for(unsigned int i = 0; i < neighbourList.size(); i++) { - cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; + + if(!silent) cout << "--- Calculating distances PWAMC <-> PWAMC" << endl; + for(unsigned int i = 0; i<limit_pwamc; i++) { + for(unsigned int j = i; j<limit_pwamc; j++) { + if(i!=j) { + distance = GetDistancePwaMcToPwaMc(i+run*limit_pwamc,j+run*limit_pwamc); + nextneighbours tmpNeighbours(i+pwaoffset,j+pwaoffset,distance,true); + neighbourList.push_back(tmpNeighbours); + // neighbourList.push_back((struct nextneighbours) {i+pwaoffset,j+pwaoffset,distance,true}); + if(neighbourList.size() % 500000 == 0 && !silent) cout << neighbourList.size() << "th distance has been calculated." << endl; + } + } } - } - if(!silent) cout << endl; - - std::sort(neighbourList.begin(),neighbourList.end(), Compare); - if(!silent) cout << "neighbourList has been sorted!" << endl; - if(!silent) cout << endl; - - if(longoutput) { - cout << "--- Distances after sorting:" << endl; - for(unsigned int i = 0; i < neighbourList.size(); i++) { - cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; + + if(!silent) { + cout << endl; + cout << "How many distances have been calculated?\t" << neighbourList.size() << endl; + cout << endl; + cout << endl; } - } - if(!silent) cout << endl; - - float numOfEntriesTotal = float(limit_data + limit_pwamc); - - unsigned int m; - - for(unsigned int i = 0; i < limit_data; i++) { - if(output) cout << "--- Checking neighbours of event " << i << endl; - m = 0; - for(unsigned int j = 0; j < neighbourList.size(); j++) - { - if (m >= numberofneighbours) continue; - if(longoutput) cout << "Checking Event in neighbourList for data: " << j << endl; - if(neighbourList[j].evtnumber1 == i || neighbourList[j].evtnumber2 == i) - { - m++; - if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i; - if(neighbourList[j].datatype==true) - { - if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; - T=T+1.; - } - else - { - if(longoutput) cout << ", but doesn't belong to the same sample." << endl; - } - } - else - { - if(longoutput) cout << "Not a neighbour of event " << i << endl; - } + + if(longoutput) { + cout << "--- Distances:" << endl; + for(unsigned int i = 0; i < neighbourList.size(); i++) { + cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; } - if(output) cout << "Event " << i << " has been analysed with " << m << " neighbours." << endl; - if(output) cout << "T: " << T << endl; - if(output) cout << endl; - } - - for(unsigned int i = 0; i < limit_pwamc; i++) { - if(output) cout << "--- Checking neighbours of event " << i+pwaoffset << endl; - m = 0; + } + if(!silent) cout << endl; - for(unsigned int j = 0; j < neighbourList.size(); j++) - { - if (m >= numberofneighbours) continue; - if(longoutput) cout << "Checking Event in neighbourList for pwamc: " << j << endl; - if(neighbourList[j].evtnumber1 == i+pwaoffset || neighbourList[j].evtnumber2 == i+pwaoffset) - { - m++; - if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i+pwaoffset; - if(neighbourList[j].datatype==true) - { - if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; - T=T+1; - } - else - { - if(longoutput) cout << ", but doesn't belong to the same sample." << endl; - } - } - else - { - if(longoutput) cout << "Not a neighbour of event " << i+pwaoffset << endl; - } - + std::sort(neighbourList.begin(),neighbourList.end(), Compare); + if(!silent) cout << "neighbourList has been sorted!" << endl; + if(!silent) cout << endl; + + if(longoutput) { + cout << "--- Distances after sorting:" << endl; + for(unsigned int i = 0; i < neighbourList.size(); i++) { + cout << neighbourList[i].evtnumber1 << "\t" << neighbourList[i].evtnumber2 << "\t" << neighbourList[i].distance << " \t\t" << neighbourList[i].datatype << endl; } - - if(output) cout << "Event " << i+pwaoffset << " has been analysed." << endl; - if(output) cout << "T: " << T << endl; - if(output) cout << endl; - } - - - if(!silent) { - cout << endl; - cout << "Number of Entries in DATA: " << numOfEntriesData << " events available (used " << limit_data << ")"<< endl; - cout << "Number of Entries in PWAMC: " << numOfEntriesPwaMc << " events available (used " << limit_pwamc << ")"<< endl; - cout << "Total Number of Entries: " << numOfEntriesTotal << endl; + } + if(!silent) cout << endl; + + float numOfEntriesTotal = float(limit_data + limit_pwamc); + + unsigned int m; + + for(unsigned int i = 0; i < limit_data; i++) { + if(output) cout << "--- Checking neighbours of event " << i << endl; + m = 0; + for(unsigned int j = 0; j < neighbourList.size(); j++) + { + if (m >= numberofneighbours) continue; + if(longoutput) cout << "Checking Event in neighbourList for data: " << j << endl; + if(neighbourList[j].evtnumber1 == i || neighbourList[j].evtnumber2 == i) + { + m++; + if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i; + if(neighbourList[j].datatype==true) + { + if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; + T=T+1.; + } + else + { + if(longoutput) cout << ", but doesn't belong to the same sample." << endl; + } + } + else + { + if(longoutput) cout << "Not a neighbour of event " << i << endl; + } + } + if(output) cout << "Event " << i << " has been analysed with " << m << " neighbours." << endl; + if(output) cout << "T: " << T << endl; + if(output) cout << endl; + } + + for(unsigned int i = 0; i < limit_pwamc; i++) { + if(output) cout << "--- Checking neighbours of event " << i+pwaoffset << endl; + m = 0; + + for(unsigned int j = 0; j < neighbourList.size(); j++) + { + if (m >= numberofneighbours) continue; + if(longoutput) cout << "Checking Event in neighbourList for pwamc: " << j << endl; + if(neighbourList[j].evtnumber1 == i+pwaoffset || neighbourList[j].evtnumber2 == i+pwaoffset) + { + m++; + if(longoutput) cout << j << "th event is the " << m << "th neighbour of event " << i+pwaoffset; + if(neighbourList[j].datatype==true) + { + if(longoutput) cout << " and belongs to the same sample! Adding +1 to T" << endl; + T=T+1; + } + else + { + if(longoutput) cout << ", but doesn't belong to the same sample." << endl; + } + } + else + { + if(longoutput) cout << "Not a neighbour of event " << i+pwaoffset << endl; + } + + } + + if(output) cout << "Event " << i+pwaoffset << " has been analysed." << endl; + if(output) cout << "T: " << T << endl; + if(output) cout << endl; + } + + + if(!silent) { + cout << endl; + cout << "Number of Entries in DATA: " << numOfEntriesData << " events available (used " << limit_data << ")"<< endl; + cout << "Number of Entries in PWAMC: " << numOfEntriesPwaMc << " events available (used " << limit_pwamc << ")"<< endl; + cout << "Total Number of Entries: " << numOfEntriesTotal << endl; + cout << endl; + } + + numOfEntriesData = limit_data; + numOfEntriesPwaMc = limit_pwamc; + + float Tnorm = 0; + Tnorm = T/(float(numberofneighbours)*float(numOfEntriesTotal)); + if(!silent) cout << "Tnorm: " << Tnorm << endl; + + + float mu = ((float(numOfEntriesData)*(float(numOfEntriesData)-1)) + (float(numOfEntriesPwaMc)*(float(numOfEntriesPwaMc)-1))) / (float(numOfEntriesTotal)*(float(numOfEntriesTotal)-1)); + if(!silent) cout << "Erwartungswert mu: " << mu << endl; + + + float variancesquared = + 1 / (float(numOfEntriesTotal)*float(numberofneighbours)) + * + ( + (float(numOfEntriesData)*float(numOfEntriesPwaMc))/(pow(float(numOfEntriesTotal), 2)) + + + 4*(pow(float(numOfEntriesData), 2)*pow(float(numOfEntriesPwaMc), 2))/(pow(float(numOfEntriesTotal),4)) + ); + + float fpull = (Tnorm-mu)/sqrt(variancesquared); + + if(!silent) { + cout << "Varianz: " << variancesquared << endl; + cout << endl; + cout << "Differenz zwischen Resultat und Erwartungswert: " << Tnorm-mu << endl; + } + cout << "f_pull in run " << run+1 << ": " << fpull << endl; cout << endl; + fpullvector.push_back(fpull); } + + std::stringstream out; + out << limit_data; + string limit_data_string = out.str(); + out.str(""); + out.clear(); + out << limit_pwamc; + string limit_pwamc_string = out.str(); + out.str(""); + out.clear(); + out << limit_runs_upper - limit_runs_lower; + string runs_string = out.str(); - numOfEntriesData = limit_data; - numOfEntriesPwaMc = limit_pwamc; - - float Tnorm = 0; - Tnorm = T/(float(numberofneighbours)*float(numOfEntriesTotal)); - if(!silent) cout << "Tnorm: " << Tnorm << endl; - - - float mu = ((float(numOfEntriesData)*(float(numOfEntriesData)-1)) + (float(numOfEntriesPwaMc)*(float(numOfEntriesPwaMc)-1))) / (float(numOfEntriesTotal)*(float(numOfEntriesTotal)-1)); - if(!silent) cout << "Erwartungswert mu: " << mu << endl; - + TString histname = "histofpull_"+limit_data_string+sourcename1+"_"+limit_pwamc_string+sourcename2+"_"+runs_string+"runs.pdf"; + cout << histname << endl; + histofpull = new TH1F("histofpull", histname, 100, -10, 10); + + cout << "--------------------------" << endl; + cout << fpullvector.size() << " fpulls have been calculated." << endl; + for(unsigned int i = 0; i < fpullvector.size(); i++) { + cout << "Calculated fpull:\t" << fpullvector[i] << endl; + histofpull->Fill( fpullvector[i] ); + } + + TCanvas* chist1 = new TCanvas("chist1","chist1",1200,800); + histofpull->SetLineColor(2); + histofpull->Draw(); - float variancesquared = - 1 / (float(numOfEntriesTotal)*float(numberofneighbours)) - * - ( - (float(numOfEntriesData)*float(numOfEntriesPwaMc))/(pow(float(numOfEntriesTotal), 2)) - + - 4*(pow(float(numOfEntriesData), 2)*pow(float(numOfEntriesPwaMc), 2))/(pow(float(numOfEntriesTotal),4)) - ); - - float fpull = (Tnorm-mu)/sqrt(variancesquared); + TF1* fit1 = new TF1("fit1","gaus(0)",-10.,10.); + fit1->SetParameters(5,histofpull->GetMean(),histofpull->GetRMS()); + fit1->SetLineColor(4); + fit1->SetLineWidth(3); + histofpull->Fit("fit1","LMRS","",-10.,10.); + + chist1->Print(histname); + histofpull->Reset(); + fpullvector.clear(); - if(!silent) { - cout << "Varianz: " << variancesquared << endl; - cout << endl; - cout << "Differenz zwischen Resultat und Erwartungswert: " << Tnorm-mu << endl; - } - cout << "f_pull in run " << run+1 << ": " << fpull << endl; - cout << endl; - fpullvector.push_back(fpull); - } - - cout << "--------------------------" << endl; - cout << fpullvector.size() << endl; - for(unsigned int i = 0; i < fpullvector.size(); i++) { - cout << "Calculated fpull:\t" << fpullvector[i] << endl; - histofpull->Fill( fpullvector[i] ); } - TCanvas* chist1 = new TCanvas("chist1","chist1",1200,800); - histofpull->SetLineColor(2); - histofpull->Draw(); - - TF1* fit1 = new TF1("fit1","gaus(0)",-10.,10.); - fit1->SetParameters(5,histofpull->GetMean(),histofpull->GetRMS()); - fit1->SetLineColor(4); - fit1->SetLineWidth(3); - histofpull->Fit("fit1","LMRS","",-10.,10.); - chist1->Print("histofpull.pdf"); } diff --git a/Examples/Psi2STo2K2PiGam/qaMixedSampleFastApp.cc b/Examples/Psi2STo2K2PiGam/qaMixedSampleFastApp.cc new file mode 100644 index 0000000000000000000000000000000000000000..1b3fb986a151e5c75a15802ce88319714a7876e5 --- /dev/null +++ b/Examples/Psi2STo2K2PiGam/qaMixedSampleFastApp.cc @@ -0,0 +1,662 @@ +#include <vector> +#include <map> +#include <string.h> +#include <math.h> +#include <stdlib.h> +#include "TFile.h" +#include "TH1F.h" +#include "TNtuple.h" +#include "TMath.h" +#include <iostream> +#include "TCanvas.h" +#include "TF1.h" +#include <sstream> + + +TFile* fdata; +TNtuple* ntdata; +TFile* fpwamc; +TNtuple* ntpwamc; + +TFile* fhistos; + +TH1F* dataCosThetaKKpipi; +TH1F* pwamcCosThetaKKpipi; +TH1F* dataPhiKKpipi; +TH1F* pwamcPhiKKpipi; +TH1F* dataCosThetaKpipi; +TH1F* pwamcCosThetaKpipi; +TH1F* dataPhiKpipi; +TH1F* pwamcPhiKpipi; +TH1F* dataCosThetapipi; +TH1F* pwamcCosThetapipi; +TH1F* dataPhipipi; +TH1F* pwamcPhipipi; +TH1F* dataCosThetapi; +TH1F* pwamcCosThetapi; +TH1F* dataPhipi; +TH1F* pwamcPhipi; +TH1F* datainvMassKpipi; +TH1F* pwamcinvMassKpipi; +TH1F* datainvMasspipi; +TH1F* pwamcinvMasspipi; + +TH1F* histofpull; + +Int_t evtNo = 0; + +float datacosThetaKKpipi; +float dataphiKKpipi; +float datacosThetaK1pi1pi2; +float dataphiK1pi1pi2; +float datacosThetapi1pi2; +float dataphipi1pi2; +float datacosThetapi1; +float dataphipi1; +float datainvmassK1pi1pi2; +float datainvmasspi1pi2; + +float pwamccosThetaKKpipi; +float pwamcphiKKpipi; +float pwamccosThetaK1pi1pi2; +float pwamcphiK1pi1pi2; +float pwamccosThetapi1pi2; +float pwamcphipi1pi2; +float pwamccosThetapi1; +float pwamcphipi1; +float pwamcinvmassK1pi1pi2; +float pwamcinvmasspi1pi2; + +float evtweight; + + +std::vector<float> veccosThetaKKpipi; +std::vector<float> vecphiKKpipi; +std::vector<float> veccosThetaK1pi1pi2; +std::vector<float> vecphiK1pi1pi2; +std::vector<float> veccosThetapi1pi2; +std::vector<float> vecphipi1pi2; +std::vector<float> veccosThetapi1; +std::vector<float> vecphipi1; +std::vector<float> vecinvmassK1pi1pi2; +std::vector<float> vecinvmasspi1pi2; +std::vector<float> vecevtweight; + +std::vector<TH1F*> histVectData; +std::vector<TH1F*> histVectPwaMc; +std::vector<TH1F*> histVecthistofpull; + +Int_t numOfEntriesData; +Int_t numOfEntriesPwaMc; + +TString type = "mc"; +bool drawFlag = true; +bool output = false; +bool longoutput = false; +bool silent = false; +TString method = "both"; +using namespace std; + +struct nextneighbours +{ + nextneighbours( unsigned int theevtnumber1, unsigned int theevtnumber2, float thedistance, int thedatatype){ + evtnumber1 = theevtnumber1; + evtnumber2 = theevtnumber2; + distance = thedistance; + datatype = thedatatype; + } + + virtual ~nextneighbours(){}; + + unsigned int evtnumber1; + unsigned int evtnumber2; + float distance; + int datatype; +}; + +void fillHistograms(); +void fillVectors(int start1, int start2, int end1, int end2); +void drawHistograms(); +void chisq_method(); +void mixed_sample(); +float GetDistance(int entry1, int entry2); +bool Compare(const nextneighbours a, const nextneighbours b); + + +int main(){ + + TString option="mixedsample,pwamc,nodraw,output"; + + //TString fileData="./Psi2STo2K2PiGamPWA_data.root"; + //TString fileData="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp1_76695events.root"; + TString fileData="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp9_43791events.root"; + TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp1_76695events.root"; + //TString filePwaMc="./Psi2STo2K2PiGamPWA_pwamc_FitParamHyp9_43791events.root"; + + cout << "--- ENTERING QAMIXEDSAMPLE" << endl; + // gROOT->SetStyle("Plain"); + // gStyle->SetCanvasColor(0); + // gStyle->SetStatBorderSize(0); + // gStyle->SetPalette(1); + // gStyle->SetOptStat(1111); + + + if (option.Contains("pwamc")) {type = "pwamc";} + if (option.Contains("nodraw")) {drawFlag = false;} + if (option.Contains("chisq")) {method = "chisq";} + if (option.Contains("mixedsample")) {method = "mixedsample";} + if (option.Contains("output")) {output = true;} + if (option.Contains("longoutput")) {longoutput = true;} + if (option.Contains("silent")) {silent = true;} + + fdata = new TFile(fileData,"READ"); + ntdata = (TNtuple*)fdata->Get("dataTuple"); + + fpwamc = new TFile(filePwaMc,"READ"); + if(type=="pwamc") { + cout << "Vergleich mit PWA-MC" << endl; + ntpwamc = (TNtuple*)fpwamc->Get("dataTuple"); + } + else if(type=="mc") { + cout << "Vergleich mit MC" << endl; + ntpwamc = (TNtuple*)fdata->Get("mcTuple"); + } + cout << endl; + + numOfEntriesData = ntdata->GetEntries(); + cout << "Events in DATA Set: " << numOfEntriesData << endl; + numOfEntriesPwaMc = ntpwamc->GetEntries(); + cout << "Events in PWAMC Set: " << numOfEntriesPwaMc << endl; + + dataCosThetaKKpipi = new TH1F("dataCosThetaKKpipi","dataCosThetaKKpipi",60,-1.,1.); + dataPhiKKpipi = new TH1F("dataPhiKKpipi","dataPhiKKpipi",60,-TMath::Pi(),TMath::Pi()); + dataCosThetaKpipi = new TH1F("dataCosThetaKpipi","dataCosThetaKpipi",60,-1.,1.); + dataPhiKpipi = new TH1F("dataPhiKpipi","dataPhiKpipi",60,-TMath::Pi(),TMath::Pi()); + dataCosThetapipi = new TH1F("dataCosThetapipi","dataCosThetapipi",60,-1.,1.); + dataPhipipi = new TH1F("dataPhipipi","dataPhipipi",60,-TMath::Pi(),TMath::Pi()); + dataCosThetapi = new TH1F("dataCosThetapi","dataCosThetapi",60,0.,1.); + dataPhipi = new TH1F("dataPhipi","dataPhipi",60,0.,TMath::Pi()); + datainvMassKpipi = new TH1F("datainvMassKpipi", "datainvMassKpipi", 60, 0.7, 3.0); + datainvMasspipi = new TH1F("datainvMasspipi", "datainvMasspipi", 60, 0.2, 2.8); + + pwamcCosThetaKKpipi = new TH1F("pwamcCosThetaKKpipi","pwamcCosThetaKKpipi",60,-1.,1.); + pwamcPhiKKpipi = new TH1F("pwamcPhiKKpipi","pwamcPhiKKpipi",60,-TMath::Pi(),TMath::Pi()); + pwamcCosThetaKpipi = new TH1F("pwamcCosThetaKpipi","pwamcCosThetaKpipi",60,-1.,1.); + pwamcPhiKpipi = new TH1F("pwamcPhiKpipi","pwamcPhiKpipi",60,-TMath::Pi(),TMath::Pi()); + pwamcCosThetapipi = new TH1F("pwamcCosThetapipi","pwamcCosThetapipi",60,-1.,1.); + pwamcPhipipi = new TH1F("pwamcPhipipi","pwamcPhipipi",60,-TMath::Pi(),TMath::Pi()); + pwamcCosThetapi = new TH1F("pwamcCosThetapi","pwamcCosThetapi",60,0.,1.); + pwamcPhipi = new TH1F("pwamcPhipi","pwamcPhipi",60,0.,TMath::Pi()); + pwamcinvMassKpipi = new TH1F("pwamcinvMassKpipi", "pwamcinvMassKpipi", 60, 0.7, 3.0); + pwamcinvMasspipi = new TH1F("pwamcinvMasspipi", "pwamcinvMasspipi", 60, 0.2, 2.8); + + histVectData.push_back(dataCosThetaKKpipi); + histVectData.push_back(dataPhiKKpipi); + histVectData.push_back(dataCosThetaKpipi); + histVectData.push_back(dataPhiKpipi); + histVectData.push_back(dataCosThetapipi); + histVectData.push_back(dataPhipipi); + histVectData.push_back(dataCosThetapi); + histVectData.push_back(dataPhipi); + histVectData.push_back(datainvMassKpipi); + histVectData.push_back(datainvMasspipi); + + histVectPwaMc.push_back(pwamcCosThetaKKpipi); + histVectPwaMc.push_back(pwamcPhiKKpipi); + histVectPwaMc.push_back(pwamcCosThetaKpipi); + histVectPwaMc.push_back(pwamcPhiKpipi); + histVectPwaMc.push_back(pwamcCosThetapipi); + histVectPwaMc.push_back(pwamcPhipipi); + histVectPwaMc.push_back(pwamcCosThetapi); + histVectPwaMc.push_back(pwamcPhipi); + histVectPwaMc.push_back(pwamcinvMassKpipi); + histVectPwaMc.push_back(pwamcinvMasspipi); + + fillHistograms(); + + if(drawFlag) { + cout << "--- Histogramme gezeichnet." << endl; + drawHistograms(); + cout << endl; + cout << endl; + cout << "-----------------------------------------" << endl; + } + if(method=="chisq" || method=="both") { + cout << "--- CHISQ-METHOD-ERGEBNISSE:" << endl; + cout << endl; + chisq_method(); + cout << "-----------------------------------------" << endl; + cout << endl; + } + if(method=="mixedsample" || method=="both") { + cout << "--- MIXED-SAMPLE-ERGEBNISSE:" << endl; + cout << endl; + mixed_sample(); + cout << "-----------------------------------------" << endl; + cout << endl; + } + cout << "--- EXITING QAMIXEDSAMPLE" << endl; + + return 1; + +} + + +bool Compare(const nextneighbours a, const nextneighbours b) { + return a.distance < b.distance; +} + + +void fillVectors(int start1, int end1, int start2, int end2) { + + veccosThetaKKpipi.clear(); + vecphiKKpipi.clear(); + veccosThetaK1pi1pi2.clear(); + vecphiK1pi1pi2.clear(); + veccosThetapi1pi2.clear(); + vecphipi1pi2.clear(); + veccosThetapi1.clear(); + vecphipi1.clear(); + vecinvmassK1pi1pi2.clear(); + vecinvmasspi1pi2.clear(); + cout << "Vectors cleared." << endl; + cout << endl; + cout << "Read branches vom ntdata." << endl; + + ntdata->SetBranchAddress("costhetaKKpipi", &datacosThetaKKpipi); + ntdata->SetBranchAddress("phiKKpipi", &dataphiKKpipi); + ntdata->SetBranchAddress("costhetaK1pi1pi2", &datacosThetaK1pi1pi2); + ntdata->SetBranchAddress("phiK1pi1pi2", &dataphiK1pi1pi2); + ntdata->SetBranchAddress("costhetapipiViaK1pipi", &datacosThetapi1pi2); + ntdata->SetBranchAddress("phipipiViaK1pipi", &dataphipi1pi2); + ntdata->SetBranchAddress("costhetapiViapipi", &datacosThetapi1); + ntdata->SetBranchAddress("phipiViapipi", &dataphipi1); + ntdata->SetBranchAddress("mk1pi1pi2", &datainvmassK1pi1pi2); + ntdata->SetBranchAddress("mpipi", &datainvmasspi1pi2); + ntpwamc->SetBranchAddress("costhetaKKpipi", &pwamccosThetaKKpipi); + ntpwamc->SetBranchAddress("phiKKpipi", &pwamcphiKKpipi); + ntpwamc->SetBranchAddress("costhetaK1pi1pi2", &pwamccosThetaK1pi1pi2); + ntpwamc->SetBranchAddress("phiK1pi1pi2", &pwamcphiK1pi1pi2); + ntpwamc->SetBranchAddress("costhetapipiViaK1pipi", &pwamccosThetapi1pi2); + ntpwamc->SetBranchAddress("phipipiViaK1pipi", &pwamcphipi1pi2); + ntpwamc->SetBranchAddress("costhetapiViapipi", &pwamccosThetapi1); + ntpwamc->SetBranchAddress("phipiViapipi", &pwamcphipi1); + ntpwamc->SetBranchAddress("mk1pi1pi2", &pwamcinvmassK1pi1pi2); + ntpwamc->SetBranchAddress("mpipi", &pwamcinvmasspi1pi2); + ntpwamc->SetBranchAddress("weight",&evtweight); + + for(int i = start1; i<end1; i++) { + ntdata->GetEntry(i); + veccosThetaKKpipi.push_back(datacosThetaKKpipi); + vecphiKKpipi.push_back(dataphiKKpipi); + veccosThetaK1pi1pi2.push_back(datacosThetaK1pi1pi2); + vecphiK1pi1pi2.push_back(dataphiK1pi1pi2); + veccosThetapi1pi2.push_back(datacosThetapi1pi2); + vecphipi1pi2.push_back(dataphipi1pi2); + veccosThetapi1.push_back(datacosThetapi1); + vecphipi1.push_back(dataphipi1); + vecinvmassK1pi1pi2.push_back(datainvmassK1pi1pi2); + vecinvmasspi1pi2.push_back(datainvmasspi1pi2); + } + cout << "Data Events filled to vectors." << endl; + + cout << "Read branches vom ntpwamc." << endl; + + for(int i = start2; i<end2; i++) { + ntpwamc->GetEntry(i); + veccosThetaKKpipi.push_back(pwamccosThetaKKpipi); + vecphiKKpipi.push_back(pwamcphiKKpipi); + veccosThetaK1pi1pi2.push_back(pwamccosThetaK1pi1pi2); + vecphiK1pi1pi2.push_back(pwamcphiK1pi1pi2); + veccosThetapi1pi2.push_back(pwamccosThetapi1pi2); + vecphipi1pi2.push_back(pwamcphipi1pi2); + veccosThetapi1.push_back(pwamccosThetapi1); + vecphipi1.push_back(pwamcphipi1); + vecinvmassK1pi1pi2.push_back(pwamcinvmassK1pi1pi2); + vecinvmasspi1pi2.push_back(pwamcinvmasspi1pi2); + } + cout << "PwaMC Events filled to vectors." << endl; + cout << "Filled " << veccosThetaKKpipi.size() << " events in total." << endl; + +} + + + +void mixed_sample() { + + if(!silent) cout << "Entering mixed-sample method" << endl; + if(!silent) cout << endl; + + std::vector<int> vec_limit_data; + //vec_limit_data.push_back(100); + // vec_limit_data.push_back(200); + // vec_limit_data.push_back(500); + vec_limit_data.push_back(1000); + // vec_limit_data.push_back(2000); + // vec_limit_data.push_back(3000); + // vec_limit_data.push_back(5000); + // vec_limit_data.push_back(7500); + // vec_limit_data.push_back(10000); + + + for(unsigned int i = 0; i<vec_limit_data.size(); i++) { + cout << "Cycle " << i << " will have " << vec_limit_data[i] << " data events." << endl; + } + + for(unsigned int k = 0; k<vec_limit_data.size(); k++) { + + //unsigned int limit_data = 2874; // numOfEntriesData (usually 2874); + unsigned int limit_data = vec_limit_data[k]; + unsigned int limit_pwamc = 1000; //10*limit_data; // numOfEntriesPwaMc; + string sourcename1 = "Hyp9pwamc"; + string sourcename2 = "Hyp1pwamc"; + unsigned int limit_runs_lower = 1; + unsigned int limit_runs_upper = 4; //numOfEntriesPwaMc/limit_pwamc; + unsigned int numberofneighbours = 10; + std::vector<float> fpullvector; + + cout << "Cycle: " << k << endl; + cout << "Using " << limit_data << " DATA events per run." << endl; + cout << "Using " << limit_pwamc << " PWAMC events per run." << endl; + cout << "Performing " << limit_runs_upper - limit_runs_lower << " runs." << endl; + + + for(unsigned int run = limit_runs_lower; run < limit_runs_upper; run++) { + + cout << "Starting run " << run << endl; + if(!silent) cout << endl; + + float distance = 1E6; + int type = 0; + float T=0; + nextneighbours maxNeighbours(0,0,1000000.,0); + nextneighbours tmpNeighbours(0,0,100000,99); + std::vector<struct nextneighbours> neighbourList; + unsigned int numOfEntriesTotalUsed = limit_data + limit_pwamc; + + fillVectors(0, limit_data, run*limit_pwamc, (run+1)*limit_pwamc); + cout << "Vectors have been filled." << endl; + + if(!silent) { + cout << endl; + cout << "Number of Entries in DATA: " << numOfEntriesData << " events available (used " << limit_data << ")"<< endl; + cout << "Number of Entries in PWAMC: " << numOfEntriesPwaMc << " events available (used " << limit_pwamc << " starting at " << run*limit_pwamc << " to " << ((run+1)*limit_pwamc)-1 << ")" << endl; + cout << "Total Number of used Entries: " << numOfEntriesTotalUsed << endl; + cout << endl; + } + + for(unsigned int i = 0; i < numOfEntriesTotalUsed; i++) { + neighbourList.clear(); + for(unsigned int nn = 0; nn < numberofneighbours; nn++) { + neighbourList.push_back(maxNeighbours); + } + for(unsigned int j = 0; j < numOfEntriesTotalUsed; j++) { + if(i!=j) { + std::sort(neighbourList.begin(),neighbourList.end(), Compare); + if(i<limit_data && j<limit_data) + { + distance = GetDistance(i,j); + type = 1; + tmpNeighbours.evtnumber1 = i; + tmpNeighbours.evtnumber2 = j; + tmpNeighbours.distance = distance; + tmpNeighbours.datatype = type; + } + else if (i<limit_data && j>=limit_data) + { + distance = GetDistance(i,j); + type = 0; + tmpNeighbours.evtnumber1 = i; + tmpNeighbours.evtnumber2 = j; + tmpNeighbours.distance = distance; + tmpNeighbours.datatype = type; + } + else if (i>=limit_data && j<limit_data) + { + distance = GetDistance(i,j); + type = 0; + tmpNeighbours.evtnumber1 = i; + tmpNeighbours.evtnumber2 = j; + tmpNeighbours.distance = distance; + tmpNeighbours.datatype = type; + } + else if (i>=limit_data && j>=limit_data) + { + distance = GetDistance(i,j); + type = 1; + tmpNeighbours.evtnumber1 = i; + tmpNeighbours.evtnumber2 = j; + tmpNeighbours.distance = distance; + tmpNeighbours.datatype = type; + } + if(distance < neighbourList[numberofneighbours-1].distance) { + neighbourList[numberofneighbours-1] = tmpNeighbours; + } + } + } + std::sort(neighbourList.begin(),neighbourList.end(), Compare); + if(longoutput) cout << neighbourList.size() << " next neighbours of event " << i << " are" << endl; + for(unsigned int k = 0; k < neighbourList.size(); k++) { + if(longoutput) cout << neighbourList[k].evtnumber1 << "\t" << neighbourList[k].evtnumber2 << "\t" << neighbourList[k].distance << " \t\t" << neighbourList[k].datatype << endl; + T=T+neighbourList[k].datatype; + } + neighbourList.clear(); + if(((i+1) % 1000) == 0 && !silent) cout << "Finished " << i+1 << "th event." << endl; + } + + float Tnorm = T/(float(numberofneighbours)*float(numOfEntriesTotalUsed)); + float mu = ((float(limit_data)*(float(limit_data)-1)) + (float(limit_pwamc)*(float(limit_pwamc)-1))) / (float(numOfEntriesTotalUsed)*(float(numOfEntriesTotalUsed)-1)); + float variancesquared = + 1 / (float(numOfEntriesTotalUsed)*float(numberofneighbours)) + * + ( + (float(limit_data)*float(limit_pwamc))/(pow(float(numOfEntriesTotalUsed), 2)) + + + 4*(pow(float(limit_data), 2)*pow(float(limit_pwamc), 2))/(pow(float(numOfEntriesTotalUsed),4)) + ); + float fpull = (Tnorm-mu)/sqrt(variancesquared); + + if(!silent) cout << "T: " << T << endl; + if(!silent) cout << "Tnorm: " << Tnorm << endl; + if(!silent) cout << "Erwartungswert mu: " << mu << endl; + if(!silent) { + cout << "Varianz: " << variancesquared << endl; + cout << endl; + cout << "Differenz zwischen Resultat und Erwartungswert: " << Tnorm-mu << endl; + } + cout << "f_pull in run " << run << ": " << fpull << endl; + cout << endl; + fpullvector.push_back(fpull); + } + + std::stringstream out; + out << limit_data; + string limit_data_string = out.str(); + out.str(""); + out.clear(); + out << limit_pwamc; + string limit_pwamc_string = out.str(); + out.str(""); + out.clear(); + out << limit_runs_upper - limit_runs_lower; + string runs_string = out.str(); + + TString histname = "histofpull_"+limit_data_string+sourcename1+"_"+limit_pwamc_string+sourcename2+"_"+runs_string+"runs.pdf"; + //cout << histname << endl; + histofpull = new TH1F("histofpull", histname, 100, -10, 10); + + cout << "--------------------------" << endl; + cout << fpullvector.size() << " fpulls have been calculated." << endl; + for(unsigned int i = 0; i < fpullvector.size(); i++) { + cout << "Calculated fpull:\t" << fpullvector[i] << endl; + histofpull->Fill( fpullvector[i] ); + } + + TCanvas* chist1 = new TCanvas("chist1","chist1",1200,800); + histofpull->SetLineColor(2); + histofpull->Draw(); + + TF1* fit1 = new TF1("fit1","gaus(0)",-10.,10.); + fit1->SetParameters(5,histofpull->GetMean(),histofpull->GetRMS()); + fit1->SetLineColor(4); + fit1->SetLineWidth(3); + histofpull->Fit("fit1","LMRS","",-10.,10.); + + chist1->Print(histname); + + histVecthistofpull.push_back(histofpull); + + TString rootfile = "histofpull_"+limit_data_string+sourcename1+"_"+limit_pwamc_string+sourcename2+"_"+runs_string+"runs.root"; + fhistos = new TFile(rootfile,"RECREATE"); + for(unsigned int i = 0; i < histVecthistofpull.size(); i++) + { + histVecthistofpull[i]->Write(); + } + + histofpull->Reset(); + fpullvector.clear(); + + + } + +} + + + + +float GetDistance(int entry1, int entry2) { + + float distance = sqrt( + pow((veccosThetaKKpipi[entry1]-veccosThetaKKpipi[entry2])/2.,2) + + pow((vecphiKKpipi[entry1]-vecphiKKpipi[entry2])/(2*3.1415),2) + + pow((veccosThetaK1pi1pi2[entry1]-veccosThetaK1pi1pi2[entry2])/2. ,2) + + pow((vecphiK1pi1pi2[entry1]-vecphiK1pi1pi2[entry2])/(2*3.1415) ,2) + + pow((veccosThetapi1pi2[entry1]-veccosThetapi1pi2[entry2])/2. ,2) + + pow((vecphipi1pi2[entry1]-vecphipi1pi2[entry2])/(2*3.1415) ,2) + + pow((veccosThetapi1[entry1]-veccosThetapi1[entry2]) ,2) + + pow((vecphipi1[entry1]-vecphipi1[entry2])/(3.1415) ,2) + + pow((vecinvmassK1pi1pi2[entry1]-vecinvmassK1pi1pi2[entry2] )/2.2 ,2) + + pow((vecinvmasspi1pi2[entry1]-vecinvmasspi1pi2[entry2])/2.3 ,2) + ); + + return distance; + +} + + + +void chisq_method() { + + cout << "Entering chisq-method" << endl; + + Double_t chisq = 0; + + Int_t constraints = 50; // 12-01-06_FitParamHyp8_base_K0K2_Hyp5/rerun_correctedampK0K2/Psi2STo2K2PiGamPWA_noK1_1680Hyp7.root + // float constraints = 94; // 12-01-06_FitParamHyp8_base_K0K2_Hyp5/rerun_correctedampK0K2/Psi2STo2K2PiGamPWA_noK1_1680Hyp7_noK1_2300Hyp7.root + Int_t ndf; + int numberofbins = 0; + + for(unsigned int i=0; i<histVectData.size(); i++) { + int currentNoOfBins=0; + + for (int j=0; j<(histVectData[i]->GetSize()); j++) { + float binContentData=float(histVectData[i]->GetBinContent(j)); + if (binContentData>0.){ + chisq += pow( binContentData - (histVectPwaMc[i]->GetBinContent(j)),2)/binContentData; + currentNoOfBins++; + } + // cout << "chi^2 after bin " << j << " of histogram " << histVectData[i]->GetName() << ": " << chisq << endl; + } + cout << "chi^2 after " << histVectData[i]->GetName() << ": " << chisq << endl; + + numberofbins += currentNoOfBins; + } + + ndf = numberofbins - constraints; + + cout << "Total chi^2:\t" << chisq << endl; + cout << "Total number of bins:\t" << numberofbins << endl; + cout << "Number of constraints:\t" << constraints << endl; + cout << "Degrees of freedom:\t" << ndf << endl; + cout << "chi^2/ndf:\t" << chisq/double(ndf) << endl; + cout << "CL:\t" << TMath::Prob(chisq, ndf) << endl; +} + + + + + + +void fillHistograms() { + + Int_t maxEvtNo = 0; + for(Int_t i=0; i<max(numOfEntriesPwaMc,numOfEntriesData); i++) + { + if(i<numOfEntriesData){ + ntdata->GetEntry(i); + + dataCosThetaKKpipi->Fill(datacosThetaKKpipi); + dataPhiKKpipi->Fill(dataphiKKpipi); + + dataCosThetaKpipi->Fill(datacosThetaK1pi1pi2); + dataPhiKpipi->Fill(dataphiK1pi1pi2); + + dataCosThetapipi->Fill(datacosThetapi1pi2); + dataPhipipi->Fill(dataphipi1pi2); + + dataCosThetapi->Fill(fabs(datacosThetapi1)); + dataPhipi->Fill(fabs(dataphipi1)); + + datainvMassKpipi->Fill(datainvmassK1pi1pi2); + datainvMasspipi->Fill(datainvmasspi1pi2); + } + + ntpwamc->GetEntry(i); + + pwamcCosThetaKKpipi->Fill(pwamccosThetaKKpipi, evtweight); + pwamcPhiKKpipi->Fill(pwamcphiKKpipi, evtweight); + + pwamcCosThetaKpipi->Fill(pwamccosThetaK1pi1pi2, evtweight); + pwamcPhiKpipi->Fill(pwamcphiK1pi1pi2, evtweight); + + pwamcCosThetapipi->Fill(pwamccosThetapi1pi2, evtweight); + pwamcPhipipi->Fill(pwamcphipi1pi2, evtweight); + + pwamcCosThetapi->Fill(fabs(pwamccosThetapi1), evtweight); + pwamcPhipi->Fill(fabs(pwamcphipi1), evtweight); + + pwamcinvMassKpipi->Fill(pwamcinvmassK1pi1pi2, evtweight); + pwamcinvMasspipi->Fill(pwamcinvmasspi1pi2, evtweight); + + if(maxEvtNo<evtNo) maxEvtNo=evtNo; + } + + +} + + + + +void drawHistograms() { + + double scalefactor = double(ntdata->GetEntries())/double(ntpwamc->GetEntries()); + + Int_t rebin = 1; + + TCanvas* cmain = new TCanvas("cmain","cmain",1400,600); + cmain->Divide(5,2); + for(unsigned int i=0; i<histVectData.size(); i++) { + cmain->cd(i+1); + histVectData[i]->Rebin(rebin); + histVectData[i]->SetLineColor(2); + histVectData[i]->Draw("E"); + histVectData[i]->SetMinimum(0); + histVectPwaMc[i]->Rebin(rebin); + histVectPwaMc[i]->Scale(scalefactor); + histVectPwaMc[i]->Draw("same"); + } + +} + + + diff --git a/Examples/Psi2STo2K2PiGam/viewHistograms.C b/Examples/Psi2STo2K2PiGam/viewHistograms.C index 2dc1b25ebdcbf2b6885a278cbfbd5311749c5bb2..a9ea41cca697bc58a4f7a41e0a039ab5b30eab93 100644 --- a/Examples/Psi2STo2K2PiGam/viewHistograms.C +++ b/Examples/Psi2STo2K2PiGam/viewHistograms.C @@ -1,12 +1,86 @@ -#include <vector> -#include <map> +#include <vector.h> +#include <map.h> #include <string.h> #include <math.h> #include <stdlib.h> -#include <TFile> -#include <TH1F> -#include <TNtuple> -#include <TCanvas> +#include <TFile.h> +#include <TH1F.h> +#include <TH2F.h> +#include <TNtuple.h> +#include <TCanvas.h> +#include <TF1.h> + +TNtuple* ntdata; + +float datainvmassK1pi1; +float datainvmassK1pi2; +float datainvmassK2pi1; +float datainvmassK2pi2; +float datainvmassK1pi1pi2; +float datainvmassK2pi1pi2; +float datainvmasspi1pi2; +float datainvmassK1K2pi1; +float datainvmassK1K2pi2; +float datainvmassK1K2; + +float datacosthetapi1; +float datacosthetapi2; +float datacosthetapipiViaK1pi1pi2; +float datacosthetapipiViaK2pi1pi2; +float datacosthetaK1K2ViaK1K2pi1; +float datacosthetaK1K2ViaK1K2pi2; + +float datacosthetaK1pi1ViaK1pi1pi2; +float datacosthetaK1pi2ViaK1pi1pi2; +float datacosthetaK2pi1ViaK2pi1pi2; +float datacosthetaK2pi2ViaK2pi1pi2; + +float datacosthetaK1ViaK1K2; +float datacosthetaK2ViaK1K2; + +std::vector<float> vecdatainvmassK1pi1; +std::vector<float> vecdatainvmassK1pi2; +std::vector<float> vecdatainvmassK2pi1; +std::vector<float> vecdatainvmassK2pi2; +std::vector<float> vecdatacosthetapi1; +std::vector<float> vecdatacosthetapi2; + +std::vector<float> vecdatainvmassK1pi1pi2; +std::vector<float> vecdatainvmassK2pi1pi2; +std::vector<float> vecdatainvmasspi1pi2; +std::vector<float> vecdatacosthetapipiViaK1pi1pi2; +std::vector<float> vecdatacosthetapipiViaK2pi1pi2; + +std::vector<float> vecdatainvmassK1K2pi1; +std::vector<float> vecdatainvmassK1K2pi2; +std::vector<float> vecdatainvmassK1K2; +std::vector<float> vecdatacosthetaK1K2ViaK1K2pi1; +std::vector<float> vecdatacosthetaK1K2ViaK1K2pi2; + +std::vector<float> vecdatacosthetaK1pi1ViaK1pi1pi2; +std::vector<float> vecdatacosthetaK1pi2ViaK1pi1pi2; +std::vector<float> vecdatacosthetaK2pi1ViaK2pi1pi2; +std::vector<float> vecdatacosthetaK2pi2ViaK2pi1pi2; + +std::vector<float> vecdatacosthetaK1ViaK1K2; +std::vector<float> vecdatacosthetaK2ViaK1K2; + +float datacosThetapi1; +float datacosThetapi2; + +Int_t numOfEntriesData; + +TH2F* datainvMassKpivsKpi; +TH1F* datainvMassKpiAngular; + +TH2F* datainvMassKpipivspipi; +TH1F* datainvMassKpipiAngular; + +TH2F* datainvMassKKpivsKK; +TH1F* datainvMassKKpiAngular; + +TH2F* datainvMassKpivsKpipi; +TH1F* datainvMassKpivsKpipiAngular; bool printToPDF = false; @@ -14,23 +88,80 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG using namespace std; gROOT->SetStyle("Plain"); - + gStyle->SetCanvasColor(0); gStyle->SetStatBorderSize(0); gStyle->SetPalette(1); gStyle->SetOptStat(1111); - + TFile* f1 = new TFile(fname,"READ"); + ntdata = (TNtuple*)f1->Get("dataTuple"); + numOfEntriesData = ntdata->GetEntries(); + + ntdata->SetBranchAddress("mk1pi1", &datainvmassK1pi1); + ntdata->SetBranchAddress("mk1pi2", &datainvmassK1pi2); + ntdata->SetBranchAddress("mk2pi1", &datainvmassK2pi1); + ntdata->SetBranchAddress("mk2pi2", &datainvmassK2pi2); + ntdata->SetBranchAddress("costhetapi1", &datacosthetapi1); + ntdata->SetBranchAddress("costhetapi2", &datacosthetapi2); + + ntdata->SetBranchAddress("mk1pi1pi2", &datainvmassK1pi1pi2); + ntdata->SetBranchAddress("mk2pi1pi2", &datainvmassK2pi1pi2); + ntdata->SetBranchAddress("mpipi", &datainvmasspi1pi2); + ntdata->SetBranchAddress("costhetapipiViaK1pipi", &datacosthetapipiViaK1pi1pi2); + ntdata->SetBranchAddress("costhetapipiViaK2pipi", &datacosthetapipiViaK2pi1pi2); + + ntdata->SetBranchAddress("mk1k2pi1", &datainvmassK1K2pi1); + ntdata->SetBranchAddress("mk1k2pi2", &datainvmassK1K2pi2); + ntdata->SetBranchAddress("mk1k2", &datainvmassK1K2); + ntdata->SetBranchAddress("costhetaKKViaK1K2pi1", &datacosthetaK1K2ViaK1K2pi1); + ntdata->SetBranchAddress("costhetaKKViaK1K2pi2", &datacosthetaK1K2ViaK1K2pi2); + + ntdata->SetBranchAddress("costhetaK1pi1", &datacosthetaK1pi1ViaK1pi1pi2); + ntdata->SetBranchAddress("costhetaK1pi2", &datacosthetaK1pi2ViaK1pi1pi2); + ntdata->SetBranchAddress("costhetaK2pi1", &datacosthetaK2pi1ViaK2pi1pi2); + ntdata->SetBranchAddress("costhetaK2pi2", &datacosthetaK2pi2ViaK2pi1pi2); + + ntdata->SetBranchAddress("costhetaK1ViaK1K2", &datacosthetaK1ViaK1K2); + + for(int i = 0; i<numOfEntriesData; i++) { + ntdata->GetEntry(i); + vecdatainvmassK1pi1.push_back(datainvmassK1pi1); + vecdatainvmassK1pi2.push_back(datainvmassK1pi2); + vecdatainvmassK2pi1.push_back(datainvmassK2pi1); + vecdatainvmassK2pi2.push_back(datainvmassK2pi2); + vecdatacosthetapi1.push_back(datacosthetapi1); + vecdatacosthetapi2.push_back(datacosthetapi2); + + vecdatainvmassK1pi1pi2.push_back(datainvmassK1pi1pi2); + vecdatainvmassK2pi1pi2.push_back(datainvmassK2pi1pi2); + vecdatainvmasspi1pi2.push_back(datainvmasspi1pi2); + vecdatacosthetapipiViaK1pi1pi2.push_back(datacosthetapipiViaK1pi1pi2); + vecdatacosthetapipiViaK2pi1pi2.push_back(datacosthetapipiViaK2pi1pi2); + + vecdatainvmassK1K2pi1.push_back(datainvmassK1K2pi1); + vecdatainvmassK1K2pi2.push_back(datainvmassK1K2pi2); + vecdatainvmassK1K2.push_back(datainvmassK1K2); + vecdatacosthetaK1K2ViaK1K2pi1.push_back(datacosthetaK1K2ViaK1K2pi1); + vecdatacosthetaK1K2ViaK1K2pi2.push_back(datacosthetaK1K2ViaK1K2pi2); + + vecdatacosthetaK1pi1ViaK1pi1pi2.push_back(datacosthetaK1pi1ViaK1pi1pi2); + vecdatacosthetaK1pi2ViaK1pi1pi2.push_back(datacosthetaK1pi2ViaK1pi1pi2); + vecdatacosthetaK2pi1ViaK2pi1pi2.push_back(datacosthetaK2pi1ViaK2pi1pi2); + vecdatacosthetaK2pi2ViaK2pi1pi2.push_back(datacosthetaK2pi2ViaK2pi1pi2); + + vecdatacosthetaK1ViaK1K2.push_back(datacosthetaK1ViaK1K2); + } std::vector<TH1F*> histVectData; std::vector<TH1F*> histVectMc; std::vector<TH2F*> histVectData2d; std::vector<TH2F*> histVectMc2d; - + bool drawswitch = false; - + if (option.Contains("makepdf")) {printToPDF = true;} - + histVectData.push_back(invKKDataHist); histVectData.push_back(invKPiDataHist); histVectData.push_back(invPiPiDataHist); @@ -39,7 +170,7 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG histVectData.push_back(cosK892DataHist); histVectData.push_back(cosK1430DataHist); histVectData.push_back(cosK1430ViaK892DataHist); - + histVectMc.push_back(invKKFittedHist); histVectMc.push_back(invKPiFittedHist); histVectMc.push_back(invPiPiFittedHist); @@ -49,6 +180,7 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG histVectMc.push_back(cosK1430FittedHist); histVectMc.push_back(cosK1430ViaK892FittedHist); + /* TCanvas* cmain = new TCanvas("cmain","cmain",1400,600); cmain->Divide(4,2); for(int i=0; i<histVectData.size(); i++) { @@ -58,20 +190,20 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG histVectData[i]->Draw("E"); histVectMc[i]->SetLineWidth(3); histVectMc[i]->Draw("same"); - } + } + + cout << endl; - cout << endl; - histVectData2d.push_back(KPivsKPiDataHist); histVectData2d.push_back(KKvsPiPiDataHist); histVectData2d.push_back(KPiPivsPiPiDataHist); histVectData2d.push_back(KKPivsKPiDataHist); - + histVectMc2d.push_back(KPivsKPiFittedHist); histVectMc2d.push_back(KKvsPiPiFittedHist); histVectMc2d.push_back(KPiPivsPiPiFittedHist); histVectMc2d.push_back(KKPivsKPiFittedHist); - + TCanvas* cmain2d = new TCanvas("cmain2d","cmain2d",1400,600); cmain2d->Divide(4,2); for(int i=0; i<histVectData2d.size(); i++) { @@ -80,16 +212,16 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG cmain2d->cd(i+1+histVectData2d.size()); histVectMc2d[i]->SetMaximum(histVectData2d[i]->GetMaximum()); histVectMc2d[i]->Draw("colz"); - } + } + + cout << endl; + - cout << endl; - - TCanvas* c1 = new TCanvas("c1","c1",1200,800); TCanvas* c2 = new TCanvas("c2","c2",1200,800); TCanvas* c3 = new TCanvas("c3","c3",1200,800); TCanvas* c4 = new TCanvas("c4","c4",1200,800); - + for(int i=0; i<4; i++) { if(i==0) c1->cd(); if(i==1) c2->cd(); @@ -113,8 +245,214 @@ void viewHistograms(TString fname="bin/gcc-4.1.2/debug/link-static/Psi2STo2K2PiG if(i==2) c3->Print("pdfplots/pwa_"+hypname+"_"+histname+".pdf"); if(i==3) c4->Print("pdfplots/pwa_"+hypname+"_"+histname+".pdf"); } + } + + */ + + datainvMassKpivsKpi = new TH2F("datainvMassKpivsKpi", "datainvMassKpivsKpi", 48, 0.5, 2.9, 48, 0.5, 2.9); + datainvMassKpiAngular = new TH1F("datainvMassKpiAngular", "datainvMassKpiAngular", 20, -1., 1.); + + datainvMassKpivsKpi_2 = new TH2F("datainvMassKpivsKpi_2", "datainvMassKpivsKpi_2", 48, 0.5, 2.9, 48, 0.5, 2.9); + datainvMassKpiAngular_2 = new TH1F("datainvMassKpiAngular_2", "datainvMassKpiAngular_2", 20, -1., 1.); + + datainvMassKpipivspipi = new TH2F("datainvMassKpipivspipi", "datainvMassKpipivspipi", 46, 0.2, 2.5, 48, 0.8, 3.2); + datainvMassKpipiAngular = new TH1F("datainvMassKpipiAngular", "datainvMassKpipiAngular", 20, -1., 1.); + + datainvMassKKpivsKK = new TH2F("datainvMassKKpivsKK", "datainvMassKKpivsKK", 48, 0.9, 3.3, 48, 1.1, 3.5); + datainvMassKKpiAngular = new TH1F("datainvMassKKpivsKKAngular", "datainvMassKKpivsKKAngular", 20, -1., 1.); + + datainvMassKpivsKpipi = new TH2F("datainvMassKpivsKpipi", "datainvMassKpivsKpipi", 48, 0.5, 2.9, 48, 0.5, 2.9); + datainvMassKpivsKpipiAngular = new TH1F("datainvMassKpivsKpipiAngular", "datainvMassKpivsKpipiAngular", 20, -1., 1.); + + datainvMasspipivsKK = new TH2F("datainvMasspipivsKK", "datainvMasspipivsKK", 48, 0.9, 3.3, 46, 0.2, 2.5); + datainvMasspipivsKKAngular = new TH1F("datainvMasspipivsKKAngular", "datainvMasspipivsKKAngular", 20, -1., 1.); + + + for(i=0; i<numOfEntriesData; i++) { + if( (vecdatainvmassK1pi1[i]>0.870 && vecdatainvmassK1pi1[i]<0.930) && (vecdatainvmassK2pi2[i]>1.8 && vecdatainvmassK2pi2[i]<2.3) ) + { + datainvMassKpivsKpi->Fill(vecdatainvmassK1pi1[i],vecdatainvmassK2pi2[i]); + datainvMassKpivsKpi->Fill(vecdatainvmassK2pi2[i],vecdatainvmassK1pi1[i]); + datainvMassKpiAngular->Fill(vecdatacosthetapi2[i]); + } + if( (vecdatainvmassK1pi2[i]>0.870 && vecdatainvmassK1pi2[i]<0.930) && (vecdatainvmassK2pi1[i]>1.8 && vecdatainvmassK2pi1[i]<2.3) ) + { + datainvMassKpivsKpi->Fill(vecdatainvmassK1pi2[i],vecdatainvmassK2pi1[i]); + datainvMassKpivsKpi->Fill(vecdatainvmassK2pi1[i],vecdatainvmassK1pi2[i]); + datainvMassKpiAngular->Fill(vecdatacosthetapi1[i]); + } + if( (vecdatainvmassK2pi1[i]>0.870 && vecdatainvmassK2pi1[i]<0.930) && (vecdatainvmassK1pi2[i]>1.8 && vecdatainvmassK1pi2[i]<2.3) ) + { + datainvMassKpivsKpi->Fill(vecdatainvmassK2pi1[i],vecdatainvmassK1pi2[i]); + datainvMassKpivsKpi->Fill(vecdatainvmassK1pi2[i],vecdatainvmassK2pi1[i]); + datainvMassKpiAngular->Fill(vecdatacosthetapi2[i]); + } + if( (vecdatainvmassK2pi2[i]>0.870 && vecdatainvmassK2pi2[i]<0.930) && (vecdatainvmassK1pi1[i]>1.8 && vecdatainvmassK1pi1[i]<2.3) ) + { + datainvMassKpivsKpi->Fill(vecdatainvmassK2pi2[i],vecdatainvmassK1pi1[i]); + datainvMassKpivsKpi->Fill(vecdatainvmassK1pi1[i],vecdatainvmassK2pi2[i]); + datainvMassKpiAngular->Fill(vecdatacosthetapi1[i]); + } + } + + TCanvas* cangular1 = new TCanvas("cangular1","cangular1",1200,600); + cangular1->Divide(2,1); + cangular1->cd(1); + datainvMassKpivsKpi->Draw("colz"); + cangular1->cd(2); + datainvMassKpiAngular->Draw(); + + // TF1* fit1 = new TF1("fit1","gaus(0)",-10.,10.); +// fit1->SetParameters(5,histofpull->GetMean(),histofpull->GetRMS()); +// fit1->SetLineColor(4); +// fit1->SetLineWidth(3); +// datainvMassKpiAngular->Fit("fit1","LMRS","",-10.,10.); + + float lower_limit = 2.50; // 2.3; + float upper_limit = 2.60; // 2.6; + + for(i=0; i<numOfEntriesData; i++) { + if( (vecdatainvmasspi1pi2[i] > 0.94 && vecdatainvmasspi1pi2[i] < 1.01 ) && (vecdatainvmassK1pi1pi2[i] > lower_limit && vecdatainvmassK1pi1pi2[i] < upper_limit) ) + { + datainvMassKpipivspipi->Fill(vecdatainvmasspi1pi2[i], vecdatainvmassK1pi1pi2[i]); + datainvMassKpipiAngular->Fill(vecdatacosthetapipiViaK1pi1pi2[i]); + } + if( (vecdatainvmasspi1pi2[i] > 0.94 && vecdatainvmasspi1pi2[i] < 1.01 ) && (vecdatainvmassK2pi1pi2[i] > lower_limit && vecdatainvmassK2pi1pi2[i] < upper_limit) ) + { + datainvMassKpipivspipi->Fill(vecdatainvmasspi1pi2[i], vecdatainvmassK2pi1pi2[i]); + datainvMassKpipiAngular->Fill(vecdatacosthetapipiViaK2pi1pi2[i]); + } + } + TCanvas* cangular2 = new TCanvas("cangular2","cangular2",1200,600); + cangular2->Divide(2,1); + cangular2->cd(1); + datainvMassKpipivspipi->Draw("colz"); + cangular2->cd(2); + datainvMassKpipiAngular->Draw(); + + TF1* fit2 = new TF1("fit2","[0]+pow([1]*(3/2*x*x-1/2)+[2],2)",-1.,1.); + fit2->SetParameters(3,5,0.3); + fit2->SetLineColor(4); + fit2->SetLineWidth(3); + datainvMassKpipiAngular->Fit("fit2","LMRS","",-1.,1.); + cangular2->Print("KpipivsKpi_2500-2600.pdf"); + + + + for(i=0; i<numOfEntriesData; i++) { + if( (vecdatainvmassK1K2[i] > 1.7 && vecdatainvmassK1K2[i] < 1.77 ) && (vecdatainvmassK1K2pi1[i] > 2.25 && vecdatainvmassK1K2pi1[i] < 2.42) ) + { + datainvMassKKpivsKK->Fill(vecdatainvmassK1K2[i], vecdatainvmassK1K2pi1[i]); + datainvMassKKpiAngular->Fill(vecdatacosthetaK1K2ViaK1K2pi1[i]); + } + if( (vecdatainvmassK1K2[i] > 1.7 && vecdatainvmassK1K2[i] < 1.77 ) && (vecdatainvmassK1K2pi2[i] > 2.25 && vecdatainvmassK1K2pi2[i] < 2.42) ) + { + datainvMassKKpivsKK->Fill(vecdatainvmassK1K2[i], vecdatainvmassK1K2pi2[i]); + datainvMassKKpiAngular->Fill(vecdatacosthetaK1K2ViaK1K2pi2[i]); + } + } + + TCanvas* cangular3 = new TCanvas("cangular3","cangular3",1200,600); + cangular3->Divide(2,1); + cangular3->cd(1); + datainvMassKKpivsKK->Draw("colz"); + cangular3->cd(2); + datainvMassKKpiAngular->Draw(); + + + + for(int i=0; i<numOfEntriesData; i++) { + if( (vecdatainvmassK1pi1[i] > 0.80 && vecdatainvmassK1pi1[i] < 0.95 ) && (vecdatainvmassK1pi1pi2[i] > 1.35 && vecdatainvmassK1pi1pi2[i] < 1.5) ) + { + datainvMassKpivsKpipi->Fill(vecdatainvmassK1pi1[i], vecdatainvmassK1pi1pi2[i]); + datainvMassKpivsKpipiAngular->Fill(vecdatacosthetaK1pi1ViaK1pi1pi2[i]); + } + if( (vecdatainvmassK1pi2[i] > 0.80 && vecdatainvmassK1pi2[i] < 0.95 ) && (vecdatainvmassK1pi1pi2[i] > 1.35 && vecdatainvmassK1pi1pi2[i] < 1.5) ) + { + datainvMassKpivsKpipi->Fill(vecdatainvmassK1pi2[i], vecdatainvmassK1pi1pi2[i]); + datainvMassKpivsKpipiAngular->Fill(vecdatacosthetaK1pi2ViaK1pi1pi2[i]); + } + if( (vecdatainvmassK2pi1[i] > 0.80 && vecdatainvmassK2pi1[i] < 0.95 ) && (vecdatainvmassK2pi1pi2[i] > 1.35 && vecdatainvmassK2pi1pi2[i] < 1.5) ) + { + datainvMassKpivsKpipi->Fill(vecdatainvmassK2pi1[i], vecdatainvmassK2pi1pi2[i]); + datainvMassKpivsKpipiAngular->Fill(vecdatacosthetaK2pi1ViaK2pi1pi2[i]); + } + if( (vecdatainvmassK2pi2[i] > 0.80 && vecdatainvmassK2pi2[i] < 0.95 ) && (vecdatainvmassK2pi1pi2[i] > 1.35 && vecdatainvmassK2pi1pi2[i] < 1.5) ) + { + datainvMassKpivsKpipi->Fill(vecdatainvmassK2pi2[i], vecdatainvmassK2pi1pi2[i]); + datainvMassKpivsKpipiAngular->Fill(vecdatacosthetaK2pi2ViaK2pi1pi2[i]); + } + + } + + + TCanvas* cangular4 = new TCanvas("cangular4","cangular4",1200,600); + cangular4->Divide(2,1); + cangular4->cd(1); + datainvMassKpivsKpipi->Draw("colz"); + cangular4->cd(2); + datainvMassKpivsKpipiAngular->Draw(); + + + + + for(i=0; i<numOfEntriesData; i++) { + if( !(vecdatainvmasspi1pi2[i]>0.9 && vecdatainvmasspi1pi2[i]<1.01)) { + if( (vecdatainvmassK1pi1[i]>1.35 && vecdatainvmassK1pi1[i]<1.5) && (vecdatainvmassK2pi2[i]>1.7 && vecdatainvmassK2pi2[i]<1.9) ) + { + datainvMassKpivsKpi_2->Fill(vecdatainvmassK1pi1[i],vecdatainvmassK2pi2[i]); + datainvMassKpivsKpi_2->Fill(vecdatainvmassK2pi2[i],vecdatainvmassK1pi1[i]); + datainvMassKpiAngular_2->Fill(vecdatacosthetapi2[i]); + } + if( (vecdatainvmassK1pi2[i]>1.35 && vecdatainvmassK1pi2[i]<1.5) && (vecdatainvmassK2pi1[i]>1.7 && vecdatainvmassK2pi1[i]<1.9) ) + { + datainvMassKpivsKpi_2->Fill(vecdatainvmassK1pi2[i],vecdatainvmassK2pi1[i]); + datainvMassKpivsKpi_2->Fill(vecdatainvmassK2pi1[i],vecdatainvmassK1pi2[i]); + datainvMassKpiAngular_2->Fill(vecdatacosthetapi1[i]); + } + if( (vecdatainvmassK2pi1[i]>1.35 && vecdatainvmassK2pi1[i]<1.5) && (vecdatainvmassK1pi2[i]>1.7 && vecdatainvmassK1pi2[i]<1.9) ) + { + datainvMassKpivsKpi_2->Fill(vecdatainvmassK2pi1[i],vecdatainvmassK1pi2[i]); + datainvMassKpivsKpi_2->Fill(vecdatainvmassK1pi2[i],vecdatainvmassK2pi1[i]); + datainvMassKpiAngular_2->Fill(vecdatacosthetapi2[i]); + } + if( (vecdatainvmassK2pi2[i]>1.35 && vecdatainvmassK2pi2[i]<1.5) && (vecdatainvmassK1pi1[i]>1.7 && vecdatainvmassK1pi1[i]<1.9) ) + { + datainvMassKpivsKpi_2->Fill(vecdatainvmassK2pi2[i],vecdatainvmassK1pi1[i]); + datainvMassKpivsKpi_2->Fill(vecdatainvmassK1pi1[i],vecdatainvmassK2pi2[i]); + datainvMassKpiAngular_2->Fill(vecdatacosthetapi1[i]); + } + } + } + + TCanvas* cangular5 = new TCanvas("cangular5","cangular5",1200,600); + cangular5->Divide(2,1); + cangular5->cd(1); + datainvMassKpivsKpi_2->Draw("colz"); + cangular5->cd(2); + datainvMassKpiAngular_2->Draw(); + + + + + for(i=0; i<numOfEntriesData; i++) { + if( (vecdatainvmasspi1pi2[i]>0.94 && vecdatainvmasspi1pi2[i]<1.01) && (vecdatainvmassK1K2[i]>2.33 && vecdatainvmassK1K2[i]<2.42) ) + { + datainvMasspipivsKK->Fill(vecdatainvmassK1K2[i],vecdatainvmasspi1pi2[i]); + datainvMasspipivsKKAngular->Fill(vecdatacosthetaK1ViaK1K2[i]); + } + } + + TCanvas* cangular6 = new TCanvas("cangular6","cangular6",1200,600); + cangular6->Divide(2,1); + cangular6->cd(1); + datainvMasspipivsKK->Draw("colz"); + cangular6->cd(2); + datainvMasspipivsKKAngular->Draw(); + + } diff --git a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc index cf105ee62a0b5e8d0ef313007548bea525b02aca..5c6f96be7db0de751ce24bd889aaa29143d82d2c 100644 --- a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc +++ b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc @@ -79,6 +79,8 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) po::options_description common("Common Options"); common.add_options() + ("datFile",po::value<string>(&_dataFile), "full path of data file") + ("mcFile",po::value<string>(&_mcFile), "full path of Monte Carlo file") ("lmax", po::value<unsigned>(&lMax)->default_value(lMax),"choose lmax.") ("pbarmom", po::value<unsigned>(&pbarMom)->default_value(pbarMom),"choose pbar momentum.") ("errLogMode,e", po::value<string>(&strErrLogMode)->default_value(strErrLogMode), @@ -88,7 +90,7 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) ("name,n", po::value<string>(&strName)->default_value("myApp"), "Name that is attached to all otuput file names to be able to run multiple fits in parallel.") ("LhMode", po::value<std::string>(&theLhMode)->default_value(theLhMode), - "Specifies the likelihood mode: OmegaPiLhGamma or OmegaTo3PiLhGamma") + "Specifies the likelihood mode: OmegaPiLhGamma, OmegaTo3PiLhGamma or OmegaTo3PiLhProd") ; po::options_description config("Configuration file options"); @@ -121,8 +123,6 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) "Specifies whether results should be returned even if they are not better than before") ("waitFactor", po::value<boost::uint32_t>(&waitFactor)->default_value(waitFactor), "Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individuum") - ("SourcePath", po::value<std::string>(&theSourcePath)->default_value(theSourcePath), - "Specifies the path to root directory of the source") ; po::options_description cmdline_options; @@ -334,23 +334,24 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) if(verbose){ - std::cout << "\nRunning with the following options using " << configFile << ":\n" - << "nProducerThreads = " << (boost::uint16_t)nProducerThreads << "\n" // boost::uint8_t not printable on gcc ??? - << "populationSize = " << populationSize << "\n" - << "nParents = " << nParents << "\n" - << "maxIterations = " << maxIterations << "\n" - << "maxMinutes = " << maxMinutes << "\n" - << "reportIteration = " << reportIteration << "\n" - << "rScheme = " << (boost::uint16_t)rScheme << "\n" - << "sortingScheme = " << smode << "\n" - << "arraySize = " << arraySize << "\n" - << "processingCycles = " << processingCycles << "\n" - << "returnRegardless = " << (returnRegardless?"true":"false") << "\n" - << "lMax = " << lMax << "\n" - << "pbarMom = " << pbarMom << "\n" - << "errLogMode = " << strErrLogMode << "\n" - << "SourcePath = " << theSourcePath << "\n" - << "LhMode = " << theLhMode << "\n" + std::cout << "\nRunning with the following options using " << configFile << ":\n" + << "data file: " << _dataFile <<"\n" + << "mc file: " << _mcFile <<"\n" + << "nProducerThreads = " << (boost::uint16_t)nProducerThreads << "\n" // boost::uint8_t not printable on gcc ??? + << "populationSize = " << populationSize << "\n" + << "nParents = " << nParents << "\n" + << "maxIterations = " << maxIterations << "\n" + << "maxMinutes = " << maxMinutes << "\n" + << "reportIteration = " << reportIteration << "\n" + << "rScheme = " << (boost::uint16_t)rScheme << "\n" + << "sortingScheme = " << smode << "\n" + << "arraySize = " << arraySize << "\n" + << "processingCycles = " << processingCycles << "\n" + << "returnRegardless = " << (returnRegardless?"true":"false") << "\n" + << "lMax = " << lMax << "\n" + << "pbarMom = " << pbarMom << "\n" + << "errLogMode = " << strErrLogMode << "\n" + << "LhMode = " << theLhMode << "\n" << endl; } diff --git a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh index 26d642b4ecbc2e2e16faafbe6b338d92f57c341d..71659d23d6f713770256a54654ce3fe843a28a93 100644 --- a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh +++ b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh @@ -71,6 +71,8 @@ class ApplicationParameterLS public: ApplicationParameterLS(int argc,char **argv) : configFile("./GOmegaPiProject.cfg") + , _dataFile("/data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/PWA_0900.dat") + , _mcFile("/data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/mcPWA_0900.dat") , parallelizationMode(1) , ip("localhost") , port(10000) @@ -97,13 +99,14 @@ class ApplicationParameterLS , lMax(3) , pbarMom(600) , errLogMode(debug) - , theSourcePath("./") - , theLhMode("OmegaPiLhGamma") + , theLhMode("OmegaTo3PiLhGamma") { if (!parseCommandLine(argc, argv)) throw false; } const std::string& getConfigFile() const { return configFile;} + const std::string dataFile() const {return _dataFile;} + const std::string mcFile() const {return _mcFile;} const boost::uint16_t& getParallelizationMode() const {return parallelizationMode;} const bool& getServerMode() const { return serverMode; } const std::string& getIp() const { return ip; } @@ -126,7 +129,6 @@ class ApplicationParameterLS const enErrLogMode& getErrLogMode() const { return errLogMode; } const Gem::Common::serializationMode& getSerMode() const { return serMode; } const enExecMode& getAppExecMode() const { return enAppExecMode; } - const std::string& getSourcePath() const { return theSourcePath; } const std::string& getLhMode() const { return theLhMode; } const std::string& getPathStartParamFile() const { return strPathStartParamFile; } const std::string& getMinuitFixParamFile() const { return strMinuitFixParamFile; } @@ -146,7 +148,9 @@ protected: std::string strName; std::string strPathStartParamFile; std::string strMinuitFixParamFile; - std::string configFile; + std::string configFile; + std::string _dataFile; + std::string _mcFile; boost::uint16_t parallelizationMode; bool serverMode; std::string ip; @@ -176,7 +180,6 @@ protected: unsigned lMax; unsigned pbarMom; enErrLogMode errLogMode; - std::string theSourcePath; std::string theLhMode; }; diff --git a/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg b/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg index bc0f456aa6442c74c22a64c3d088c0017004e47e..c9b7e8f406a840de2b518b15c39914e5eb2380e6 100644 --- a/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg +++ b/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg @@ -3,7 +3,8 @@ # # # Authors: Ruediger Berlich # ################################################################################# - +datFile = /data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/PWA_0900.dat +mcFile = /data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/mcPWA_0900.dat # The amount of random number producer threads nProducerThreads = 16 #nProducerThreads = 4 @@ -38,7 +39,8 @@ returnRegardless = 1 # Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individual waitFactor = 2 # Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individual -lmax = 2 -pbarmom = 600 +lmax = 4 +pbarmom = 900 errLogMode = debug -SourcePath = /data/sleipnir1/bertram/PawianGit1107/Pawian/ +appMode = 1 +LhMode = OmegaTo3PiLhProd \ No newline at end of file diff --git a/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc b/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc index 733849ae035b95ac51d19687587f94c81add8596..509bddaa577bcb09f387023f2ed021000852ea17 100644 --- a/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc +++ b/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc @@ -18,6 +18,7 @@ #include "Examples/pbarpToOmegaPiLS/AbsOmegaPiLhLS.hh" #include "Examples/pbarpToOmegaPiLS/OmegaPiLhPi0GammaLS.hh" #include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhPi0GammaLS.hh" +#include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh" #include "Examples/pbarpToOmegaPiLS/MOmegaPiFcnLS.hh" #include "Examples/pbarpToOmegaPiLS/SpinDensityHistLS.hh" @@ -614,21 +615,20 @@ int main(int __argc,char *__argv[]){ exit(1); } - std::string piomegaDatFile; - std::string piomegaMcFile; + std::string piomegaDatFile=theAppParams.dataFile(); + std::string piomegaMcFile=theAppParams.mcFile(); + Info << "data file: " << piomegaDatFile << endmsg; + Info << "mc file: " << piomegaMcFile << endmsg; + int nParticlesPerEvt=0; bool readWeightData=true; if (theAppParams.getLhMode()=="OmegaPiLhGamma" || (theAppParams.getLhMode()=="OmegaPiLhGammaBw") ){ - nParticlesPerEvt=3; - constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile); - constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile); - } - else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ - nParticlesPerEvt=4; - constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/om3pi_",theAppParams.getPbarMom(),piomegaDatFile); - constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mcom3pi_",theAppParams.getPbarMom(),piomegaMcFile); - } + nParticlesPerEvt=3; + } + else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhProd") || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ + nParticlesPerEvt=4; + } else { Alert << "LhMode: " << theAppParams.getLhMode() << " does not exist!!!" << endmsg; @@ -636,7 +636,6 @@ int main(int __argc,char *__argv[]){ } - // constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile); if (checkFileExist(piomegaDatFile)){ Info << "Using Data file " << piomegaDatFile << endmsg; } @@ -646,7 +645,6 @@ int main(int __argc,char *__argv[]){ exit(1); } - // constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile); if (checkFileExist(piomegaMcFile)){ Info << "Using Monte Carlo file " << piomegaMcFile << endmsg; } @@ -661,7 +659,9 @@ int main(int __argc,char *__argv[]){ ParticleTable pTable; PdtParser parser; - std::string pdtFile(theAppParams.getSourcePath()+"/Particle/pdt.table"); + + std::string theSourcePath=getenv("CMAKE_SOURCE_DIR"); + std::string pdtFile(theSourcePath+"/Particle/pdt.table"); if (!parser.parse(pdtFile, pTable)) { Alert << "Error: could not parse " << pdtFile << endmsg; exit(1); @@ -709,11 +709,10 @@ int main(int __argc,char *__argv[]){ theOmegaPiEventPtr = boost::shared_ptr<const AbsOmegaPiEventListLS>(new OmegaPiEventListLS (theDataEventList, theMcEventList, theAppParams.getLMax()+1, theAppParams.getPbarMom() ) ); theRootFilePath << "./" << theAppParams.getName() << "OmegaPi0Fit_Lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root"; } - else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ + else if ( (theAppParams.getLhMode()=="OmegaTo3PiLhGamma") || (theAppParams.getLhMode()=="OmegaTo3PiLhProd") || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw") ){ theOmegaPiEventPtr = boost::shared_ptr<const AbsOmegaPiEventListLS>(new OmegaTo3PiEventListLS (theDataEventList, theMcEventList, theAppParams.getLMax()+1, theAppParams.getPbarMom() ) ); theRootFilePath << "./" << theAppParams.getName() << "OmegaTo3PiFit_Lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root"; } - if (theAppParams.getAppExecMode()==ApplicationParameterLS::HistTest){ histTest(theOmegaPiEventPtr, theRootFilePath.str()); exit(1); @@ -737,6 +736,9 @@ int main(int __argc,char *__argv[]){ else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma"){ theOmegaPiLh = boost::shared_ptr<AbsOmegaPiLhLS>(new OmegaTo3PiLhPi0GammaLS(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr)); } + else if (theAppParams.getLhMode()=="OmegaTo3PiLhProd"){ + theOmegaPiLh = boost::shared_ptr<AbsOmegaPiLhLS>(new OmegaTo3PiLhProdLS(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr)); + } else{ Alert <<"LhMode " << theAppParams.getLhMode() << " doesn't exist !!! " << endmsg; exit(1); diff --git a/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc new file mode 100644 index 0000000000000000000000000000000000000000..0a425aba9e092b8527bd9ecafe89a5741a67d870 --- /dev/null +++ b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc @@ -0,0 +1,84 @@ +#include <getopt.h> + +#include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh" +#include "Examples/pbarpToOmegaPiLS/AbsOmegaPiEventListLS.hh" +#include "PwaUtils/pbarpStatesLS.hh" +#include "Examples/pbarpToOmegaPiLS/pbarpToOmegaPi0StatesLS.hh" +#include "ErrLogger/ErrLogger.hh" + +// #include <geneva/GConstrainedDoubleObject.hpp> + + +OmegaTo3PiLhProdLS::OmegaTo3PiLhProdLS(boost::shared_ptr<const AbsOmegaPiEventListLS> theEvtList, boost::shared_ptr<const pbarpToOmegaPi0StatesLS> theStates) : + AbsOmegaPiLhLS(theEvtList, theStates, "OmegaTo3PiLhProd") +{ +} + +OmegaTo3PiLhProdLS::OmegaTo3PiLhProdLS(boost::shared_ptr<OmegaTo3PiLhProdLS> theOmegaTo3PiLhProdLSPtr): + AbsOmegaPiLhLS(theOmegaTo3PiLhProdLSPtr, "OmegaTo3PiLhProd") +{ +} + +OmegaTo3PiLhProdLS::~OmegaTo3PiLhProdLS() +{ +} + +double OmegaTo3PiLhProdLS::calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal){ + + double result=AbsOmegaPiLhLS::calcLogLh(theParamVal); + if (_globalItCounter%10 == 0) printFitParams(std::cout, theParamVal); + +// if (_globalItCounter%1000 == 0){ +// std::ofstream theStream ("currentResult.dat"); +// dumpCurrentResult(theStream, theParamVal); +// theStream.close(); +// } + + return result; +} + + +double OmegaTo3PiLhProdLS::calcEvtIntensity(OmegaPiDataLS::OmPiEvtDataLS* theData, const OmegaPiDataLS::fitParamVal& theParamVal){ + + double result=0.; + + for (Spin lamOmega=-1; lamOmega<=1; ++lamOmega){ + result+= norm(calcCoherentAmp(lamOmega,0, theParamVal, _singlet_JPCLS_States, theData)); //singlet amp + result+= norm(calcCoherentAmp(lamOmega,0, theParamVal, _triplet0_JPCLS_States, theData)); //triplet0 + result+= norm(calcCoherentAmp(lamOmega,1, theParamVal, _tripletp1_JPCLS_States, theData)); //triplet+1 + result+= norm(calcCoherentAmp(lamOmega,-1, theParamVal, _tripletm1_JPCLS_States, theData)); //triplet-1 + } + + return result; +} + + +complex<double> OmegaTo3PiLhProdLS::calcCoherentAmp(Spin lamOm, Spin Minit, const OmegaPiDataLS::fitParamVal& theParamVal, std::vector< boost::shared_ptr<const JPCLSls> >& theJPCLSlsStates, OmegaPiDataLS::OmPiEvtDataLS* theData){ + + complex<double> result(0.,0.); + + std::vector< boost::shared_ptr<const JPCLSls> >::const_iterator it; + for ( it=theJPCLSlsStates.begin(); it!=theJPCLSlsStates.end(); ++it){ + + + std::map< boost::shared_ptr<const JPCLSls>, pair<double, double>, pawian::Collection::SharedPtrLess >::const_iterator itmap; + itmap = theParamVal.lsParam.find((*it)); + pair<double, double> fitPair= (*itmap).second; + + if (fabs(lamOm)> (*it)->J ) continue; + + complex<double> omegaProdAmp=calcOmegaProdPartAmp(Minit, lamOm, (*it), fitPair, theData); + + result += omegaProdAmp; + + } + + + + return result; + } + + +void OmegaTo3PiLhProdLS::print(std::ostream& os) const{ + os << "OmegaTo3PiLhProdLS\n"; +} diff --git a/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh new file mode 100644 index 0000000000000000000000000000000000000000..f0f0aad8e81d01f5de67aaec056713969ca63df3 --- /dev/null +++ b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh @@ -0,0 +1,63 @@ +#ifndef _OmegaTo3PiLhProdLS_H +#define _OmegaTo3PiLhProdLS_H + +#include <iostream> +#include <fstream> +#include <string> +#include <vector> +#include <complex> + +#include <cassert> +#include <boost/shared_ptr.hpp> + +#include "TROOT.h" +// #include <TSystem.h> +#include "qft++/topincludes/relativistic-quantum-mechanics.hh" +#include "Examples/pbarpToOmegaPiLS/OmegaPiDataLS.hh" +#include "Examples/pbarpToOmegaPiLS/AbsOmegaPiLhLS.hh" + + +#include "PwaUtils/DataUtils.hh" + +// using namespace std; + + +class AbsOmegaPiEventListLS; +class pbarpToOmegaPi0StatesLS; + +class OmegaTo3PiLhProdLS : public AbsOmegaPiLhLS{ + +public: + + // create/copy/destroy: + + ///Constructor + OmegaTo3PiLhProdLS(boost::shared_ptr<const AbsOmegaPiEventListLS>, boost::shared_ptr<const pbarpToOmegaPi0StatesLS>); + OmegaTo3PiLhProdLS(boost::shared_ptr<OmegaTo3PiLhProdLS>); + + /** Destructor */ + virtual ~OmegaTo3PiLhProdLS(); + + virtual AbsOmegaPiLhLS* clone_() const{ + return new OmegaTo3PiLhProdLS(_omegaPiEventListPtr, _omegaPi0StatesPtr); + } + + + // Getters: + virtual double calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal); + virtual double calcEvtIntensity(OmegaPiDataLS::OmPiEvtDataLS* theData, const OmegaPiDataLS::fitParamVal& theParamVal); + + + virtual void print(std::ostream& os) const; + +protected: + virtual complex<double> calcCoherentAmp(Spin lamOm, Spin Minit, const OmegaPiDataLS::fitParamVal& theParamVal, std::vector< boost::shared_ptr<const JPCLSls> >& theJPCLSlsStates, OmegaPiDataLS::OmPiEvtDataLS* theData); + + + +private: + + +}; + +#endif diff --git a/PspGen/PspGenTestApp.cc b/PspGen/PspGenTestApp.cc index 64358f0cff02b0c8415951f2419db4f585e25ed3..04c9c9f8528f86d34a814e1bd09d8e5bbdbbf6d8 100644 --- a/PspGen/PspGenTestApp.cc +++ b/PspGen/PspGenTestApp.cc @@ -30,7 +30,7 @@ int main(int argc, char* argv[]) int eventsperset = 500000; bool silent = true; - for(int set=100; set < set_limit; set++) { + for(int set=0; set < set_limit; set++) { std::stringstream out; out << set+1001; diff --git a/qft++/relativistic-quantum-mechanics/BlattWeisskopf.cc b/qft++/relativistic-quantum-mechanics/BlattWeisskopf.cc index 0d9082eae76307fc3ad23ec1002b6e22b58d4ecb..7d24affc175ea87a17fa00f6f527e7bd55268349 100644 --- a/qft++/relativistic-quantum-mechanics/BlattWeisskopf.cc +++ b/qft++/relativistic-quantum-mechanics/BlattWeisskopf.cc @@ -78,7 +78,7 @@ double BlattWeisskopf::compute(double p) const if(0 == _LL) result=1.; else - if(1 == _LL) result=sqrt(x/(1.0+x)); + if(1 == _LL) result=sqrt(2.*x/(1.0+x)); else if(2 == _LL) result=sqrt((13.*x*x)/((x-3.)*(x-3.)+9.*x)); else