diff --git a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc
index e1f64c7fb546e00ec093bb1ef39327f125af842b..fa341074a75d3399cbf11cf18bccdc356d4487fa 100644
--- a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc
+++ b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.cc
@@ -1170,7 +1170,7 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToPi0PiToKstarKAmp(Psi2STo2K2PiGamData
}
-complex<double> AbsPsi2STo2K2PiGamLh::chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi_2_Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& Pi_2_Dec, double Pi_2_Mass, double Pi_2_Width, double Kstar_Mass, double Kstar_Width){
+complex<double> AbsPsi2STo2K2PiGamLh::chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi_2_Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& Pi_2_Dec, double Pi_2_Mass, double Pi_2_Width, double Kstar_Mass, double Kstar_Width, Spin Kstar_Spin){
Vector4<double> KKPi0(theData->KKPi0_HeliChic0_4V.E(), theData->KKPi0_HeliChic0_4V.Px(),
theData->KKPi0_HeliChic0_4V.Py(), theData->KKPi0_HeliChic0_4V.Pz());
@@ -1194,6 +1194,8 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamDat
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator it;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >::iterator itDec;
+
+ int L_KstarDec=int(Kstar_Spin);
for ( it=ChiToPi_2_Prod.begin(); it!=ChiToPi_2_Prod.end(); ++it){
@@ -1210,18 +1212,22 @@ complex<double> AbsPsi2STo2K2PiGamLh::chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamDat
double thePhiDec=itDec->second.second;
complex<double> expiphiDec(cos(thePhiDec), sin(thePhiDec));
complex<double> tmpDec(0.,0.);
- for (Spin lamKstar=-1; lamKstar<=1; ++lamKstar){
- tmpDec+=sqrt(2.*theDecState->L+1)*Clebsch(theDecState->L, 0., theDecState->S, lamKstar, theDecState->J, lamKstar)*Clebsch(1,lamKstar, 0, 0, theDecState->S, lamKstar)*sqrt(3.)
+
+ for (Spin lamKstar=-Kstar_Spin; lamKstar<=Kstar_Spin; ++lamKstar){
+
+ if( fabs(lamKstar)>theDecState->J || fabs(lamKstar)>theDecState->S) continue;
+
+ tmpDec+=sqrt(2.*theDecState->L+1)*Clebsch(theDecState->L, 0., theDecState->S, lamKstar, theDecState->J, lamKstar)*Clebsch(Kstar_Spin,lamKstar, 0, 0, theDecState->S, lamKstar)*sqrt(2.*L_KstarDec+1.)
*(( BreitWigner(KKPi0, Pi_2_Mass, Pi_2_Width)
- *( BreitWignerBlattW(Kstarp_pi0, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,1)
- *theData->DfPi2ToKstarpK_pi0[theDecState->J][lamPi2][lamKstar]*theData->DfKst1pToKpPi0ViaKKPi0[1][lamKstar][0]
- +BreitWignerBlattW(Kstarm_pi0, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,1)
- *theData->DfPi2ToKstarmK_pi0[theDecState->J][lamPi2][lamKstar]*theData->DfKst1mToKmPi0ViaKKPi0[1][lamKstar][0]))
+ *( BreitWignerBlattW(Kstarp_pi0, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,L_KstarDec)
+ *theData->DfPi2ToKstarpK_pi0[theDecState->J][lamPi2][lamKstar]*theData->DfKst1pToKpPi0ViaKKPi0[Kstar_Spin][lamKstar][0]
+ +BreitWignerBlattW(Kstarm_pi0, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,L_KstarDec)
+ *theData->DfPi2ToKstarmK_pi0[theDecState->J][lamPi2][lamKstar]*theData->DfKst1mToKmPi0ViaKKPi0[Kstar_Spin][lamKstar][0]))
+( BreitWigner(KKPi1, Pi_2_Mass, Pi_2_Width)
- *( BreitWignerBlattW(Kstarp_pi1, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,1)
- *theData->DfPi2ToKstarpK_pi1[theDecState->J][lamPi2][lamKstar]*theData->DfKst1pToKpPi1ViaKKPi1[1][lamKstar][0]
- +BreitWignerBlattW(Kstarm_pi1, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,1)
- *theData->DfPi2ToKstarmK_pi1[theDecState->J][lamPi2][lamKstar]*theData->DfKst1mToKmPi1ViaKKPi1[1][lamKstar][0] ))
+ *( BreitWignerBlattW(Kstarp_pi1, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,L_KstarDec)
+ *theData->DfPi2ToKstarpK_pi1[theDecState->J][lamPi2][lamKstar]*theData->DfKst1pToKpPi1ViaKKPi1[Kstar_Spin][lamKstar][0]
+ +BreitWignerBlattW(Kstarm_pi1, 0.493677, 0.1349766, Kstar_Mass, Kstar_Width,L_KstarDec)
+ *theData->DfPi2ToKstarmK_pi1[theDecState->J][lamPi2][lamKstar]*theData->DfKst1mToKmPi1ViaKKPi1[Kstar_Spin][lamKstar][0] ))
);
}
diff --git a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh
index 0f82074c97922ebdb174b45ebe6380947af05820..9464615cffe2c395351543e7b54c56123f3de4b1 100644
--- a/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh
+++ b/Examples/Psi2STo2K2PiGam/AbsPsi2STo2K2PiGamLh.hh
@@ -118,7 +118,7 @@ protected:
virtual complex<double> chiToPi0PiToKstarKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi_0_Prod, double Pi_0_Mass, double Pi_0_Width, double Kstar_Mass, double Kstar_Width);
- virtual complex<double> chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi_2_Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& Pi_2_Dec, double Pi_2_Mass, double Pi_2_Width, double Kstar_Mass, double Kstar_Width);
+ virtual complex<double> chiToPi2Pi0ToKstarKAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi_2_Prod, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& Pi_2_Dec, double Pi_2_Mass, double Pi_2_Width, double Kstar_Mass, double Kstar_Width, Spin Kstar_Spin=1);
virtual complex<double> chiToPi0Pi0Tof980kkAmp(Psi2STo2K2PiGamData::Psi2STo2K2PiGamEvtData* theData, std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess >& ChiToPi0Pi0Amp, double f980_Mass, double f980_gKK, double f980_gPiPi, double Pi_0_Mass, double Pi_0_Width);
diff --git a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc
index d859ace82991507bb8d1542b8f7f426877dab178..dbe295df688064f41dbb59fa500cf09a9bb810db 100644
--- a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc
+++ b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.cc
@@ -80,6 +80,7 @@ FitParams2K2PiGam::FitParams2K2PiGam()
_jpclsMap[paramEnum2K2PiGam::ChiToPi_2_2285Pi]=theStates.ChiToPi_2PiStates();
_jpclsMap[paramEnum2K2PiGam::Pi_2_2285Tof1700Pi]=theStates.Pi_2Tof0PiStates();
_jpclsMap[paramEnum2K2PiGam::Pi_2_2285ToK892KPi]=theStates.Pi_2ToKst1KStates();
+ _jpclsMap[paramEnum2K2PiGam::Pi_2_2285ToK_2_1430K]=theStates.Pi_2ToKst2KStates();
_jpclsMap[paramEnum2K2PiGam::ChiToK_2_1770K]=theStates.ChiToK2mK0mStates();
_jpclsMap[paramEnum2K2PiGam::K_2_1770ToK_2_1430Pi]=theStates.K2mToK2pPiStates();
_jpclsMap[paramEnum2K2PiGam::ChiToK_0_1430KPi]=theStates.ChiTo2K_0_States();
@@ -173,6 +174,7 @@ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collecti
else if (index==paramEnum2K2PiGam::ChiToPi_2_2285Pi) return params.ChiToPi_2_2285Pi;
else if (index==paramEnum2K2PiGam::Pi_2_2285Tof1700Pi) return params.Pi_2_2285Tof1700Pi;
else if (index==paramEnum2K2PiGam::Pi_2_2285ToK892KPi) return params.Pi_2_2285ToK892KPi;
+ else if (index==paramEnum2K2PiGam::Pi_2_2285ToK_2_1430K) return params.Pi_2_2285ToK_2_1430K;
else if (index==paramEnum2K2PiGam::ChiToK_2_1770K) return params.ChiToK_2_1770K;
else if (index==paramEnum2K2PiGam::K_2_1770ToK_2_1430Pi) return params.K_2_1770ToK_2_1430Pi;
else if (index==paramEnum2K2PiGam::ChiToK_0_1430KPi) return params.ChiToK_0_1430KPi;
diff --git a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh
index 801506bfef1d4093b8258a55644163f42bc473a5..70d6076f7d9bc84865ebc8bae829ba6be5f7c189 100644
--- a/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh
+++ b/Examples/Psi2STo2K2PiGam/FitParams2K2PiGam.hh
@@ -93,6 +93,7 @@ struct param2K2PiGam
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToPi_2_2285Pi;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > Pi_2_2285Tof1700Pi;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > Pi_2_2285ToK892KPi;
+ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > Pi_2_2285ToK_2_1430K;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_2_1770K;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > K_2_1770ToK_2_1430Pi;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToK_0_1430KPi;
@@ -148,7 +149,8 @@ struct paramEnum2K2PiGam{
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,
ChiToPi_2_1670Pi, Pi_2_1670Tof_2_1270Pi, Pi_2_1670ToK892K, ChiToPi1800Pi0Tof980, ChiToPi1800Pi0Tof1370, ChiToPi1800Pi0ToKappa, ChiToPi1800Pi0ToK892K,
- ChiToPi_2_2285Pi,Pi_2_2285Tof1700Pi, Pi_2_2285ToK892KPi, ChiToK_2_1770K, K_2_1770ToK_2_1430Pi, ChiToK_0_1430KPi, nAmps,
+ ChiToPi_2_2285Pi,Pi_2_2285Tof1700Pi, Pi_2_2285ToK892KPi, Pi_2_2285ToK_2_1430K,
+ ChiToK_2_1770K, K_2_1770ToK_2_1430Pi, ChiToK_0_1430KPi, 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,
@@ -173,7 +175,8 @@ struct paramEnum2K2PiGam{
"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",
"ChiToPi_2_1670Pi","Pi_2_1670Tof_2_1270Pi","Pi_2_1670ToK892K","ChiToPi1800Pi0Tof980","ChiToPi1800Pi0Tof1370","ChiToPi1800Pi0ToKappa","ChiToPi1800Pi0ToK892K",
- "ChiToPi_2_2285Pi","Pi_2_2285Tof1700Pi","Pi_2_2285ToK892KPi","ChiToK_2_1770K","K_2_1770ToK_2_1430Pi","ChiToK_0_1430KPi",
+ "ChiToPi_2_2285Pi","Pi_2_2285Tof1700Pi","Pi_2_2285ToK892KPi","Pi_2_2285ToK_2_1430K",
+ "ChiToK_2_1770K","K_2_1770ToK_2_1430Pi","ChiToK_0_1430KPi",
"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",
diff --git a/Examples/Psi2STo2K2PiGam/Hyp9Lh.cc b/Examples/Psi2STo2K2PiGam/Hyp9Lh.cc
index 109e29fcf1c5c2f12ff7957d07de00f1f95cd116..c6493fb9d9ea5537ee399727571c3e4199f69d15 100644
--- a/Examples/Psi2STo2K2PiGam/Hyp9Lh.cc
+++ b/Examples/Psi2STo2K2PiGam/Hyp9Lh.cc
@@ -17,6 +17,7 @@ Hyp9Lh::Hyp9Lh(boost::shared_ptr<const Psi2STo2K2PiGamEvtList> theEvtList, const
,_Pi1800Pi0ToK892KHyp9(true)
,_Pi_2_2285Tof1710PiHyp9(true)
,_Pi_2_2285ToK892KHyp9(true)
+ ,_Pi_2_2285ToK_2_1430KHyp9(true)
,_f980f_2_2300Hyp9(true)
,_f_2_2300sigmaHyp9(true)
,_K_2_1770ToK_2_1430PiHyp9(true)
@@ -36,6 +37,7 @@ Hyp9Lh::Hyp9Lh( boost::shared_ptr<AbsPsi2STo2K2PiGamLh> theLhPtr, const std::map
,_Pi1800Pi0ToK892KHyp9(true)
,_Pi_2_2285Tof1710PiHyp9(true)
,_Pi_2_2285ToK892KHyp9(true)
+ ,_Pi_2_2285ToK_2_1430KHyp9(true)
,_f980f_2_2300Hyp9(true)
,_f_2_2300sigmaHyp9(true)
,_K_2_1770ToK_2_1430PiHyp9(true)
@@ -125,7 +127,7 @@ complex<double> Hyp9Lh::chi0DecAmps(const param2K2PiGam& theParamVal, Psi2STo2K2
}
- if (_Pi_2_2285Tof1710PiHyp9 || _Pi_2_2285ToK892KHyp9){
+ if (_Pi_2_2285Tof1710PiHyp9 || _Pi_2_2285ToK892KHyp9 || _Pi_2_2285ToK_2_1430KHyp9){
double Pi_2_2285Mass=theParamVal.BwPi_2_2285.first;
double Pi_2_2285Width=theParamVal.BwPi_2_2285.second;
std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > ChiToPi_2_2285Pi=theParamVal.ChiToPi_2_2285Pi;
@@ -145,6 +147,15 @@ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collecti
result+=chiToPi2Pi0ToKstarKAmp(theData, ChiToPi_2_2285Pi, Pi_2_2285ToK892KPi, Pi_2_2285Mass, Pi_2_2285Width, K892Mass, K892Width);
}
+
+
+ if(_Pi_2_2285ToK_2_1430KHyp9){
+ std::map< boost::shared_ptr<const JPCLS>, pair<double, double>, pawian::Collection::SharedPtrLess > Pi_2_2285ToK_2_1430K=theParamVal.Pi_2_2285ToK_2_1430K;
+ double K_2_1430Mass=theParamVal.BwK_2_1430.first;
+ double K_2_1430Width=theParamVal.BwK_2_1430.second;
+
+ result+=chiToPi2Pi0ToKstarKAmp(theData, ChiToPi_2_2285Pi, Pi_2_2285ToK_2_1430K, Pi_2_2285Mass, Pi_2_2285Width, K_2_1430Mass, K_2_1430Width,2 );
+ }
}
if (_f980f_2_2300Hyp9 || _f_2_2300sigmaHyp9){
@@ -408,6 +419,19 @@ void Hyp9Lh::setUp(const std::map<const std::string, bool>& hypMap){
exit(0);
}
+ iter= hypMap.find("Pi_2_2285ToK_2_1430KHyp9");
+
+ if (iter !=hypMap.end()){
+ _Pi_2_2285ToK_2_1430KHyp9= iter->second;
+ Info<< "hypothesis " << iter->first << "\t" << _Pi_2_2285ToK_2_1430KHyp9 <<endmsg;
+ _hypMap[iter->first]= iter->second;
+ }
+ else{
+ Alert << "Pi_2_2285ToK_2_1430KHyp9 not set!!!" <<endmsg;
+ exit(0);
+ }
+
+
iter= hypMap.find("f980f_2_2300Hyp9");
@@ -488,7 +512,7 @@ void Hyp9Lh::setUp(const std::map<const std::string, bool>& hypMap){
}
- if (_Pi_2_2285Tof1710PiHyp9 || _Pi_2_2285ToK892KHyp9){
+ if (_Pi_2_2285Tof1710PiHyp9 || _Pi_2_2285ToK892KHyp9 || _Pi_2_2285ToK_2_1430KHyp9){
_ampVec.push_back(paramEnum2K2PiGam::ChiToPi_2_2285Pi);
_massVec.push_back(paramEnum2K2PiGam::Pi_2_2285);
@@ -499,6 +523,11 @@ void Hyp9Lh::setUp(const std::map<const std::string, bool>& hypMap){
if (_Pi_2_2285ToK892KHyp9){
_ampVec.push_back(paramEnum2K2PiGam::Pi_2_2285ToK892KPi);
}
+
+ if(_Pi_2_2285ToK_2_1430KHyp9){
+ _ampVec.push_back(paramEnum2K2PiGam::Pi_2_2285ToK_2_1430K);
+ if(!_K2_1430_K2_1430Hyp && !_K0_1430_K2_1430Hyp) _massVec.push_back(paramEnum2K2PiGam::K_2_1430);
+ }
}
if (_f980f_2_2300Hyp9 || _f_2_2300sigmaHyp9){
diff --git a/Examples/Psi2STo2K2PiGam/Hyp9Lh.hh b/Examples/Psi2STo2K2PiGam/Hyp9Lh.hh
index 94358e5ce9560139943d757636b8febdd62081fc..e79b1cca3f447147c9c401b28fc82200562d0224 100644
--- a/Examples/Psi2STo2K2PiGam/Hyp9Lh.hh
+++ b/Examples/Psi2STo2K2PiGam/Hyp9Lh.hh
@@ -58,6 +58,7 @@ protected:
bool _Pi1800Pi0ToK892KHyp9;
bool _Pi_2_2285Tof1710PiHyp9;
bool _Pi_2_2285ToK892KHyp9;
+ bool _Pi_2_2285ToK_2_1430KHyp9;
bool _f980f_2_2300Hyp9;
bool _f_2_2300sigmaHyp9;
bool _K_2_1770ToK_2_1430PiHyp9;
diff --git a/Examples/Psi2STo2K2PiGam/Jamfile b/Examples/Psi2STo2K2PiGam/Jamfile
index e0e81a396caa66ff5e0d08664e0829b2c0ca3899..f0a330b8a7290fd63dfb68f412e826f18d5e4969 100644
--- a/Examples/Psi2STo2K2PiGam/Jamfile
+++ b/Examples/Psi2STo2K2PiGam/Jamfile
@@ -13,3 +13,5 @@ lib Psi2STo2K2PiGam :
exe ParserTestApp : ParserTestApp.cc Psi2STo2K2PiGam : ;
exe Mpsi2STo2K2PiGamTestApp : Mpsi2STo2K2PiGamTestApp.cc Psi2STo2K2PiGam : ;
+exe qaMixedSampleApp : qaMixedSampleApp.cc Psi2STo2K2PiGam : ;
+exe qaMixedSampleCacheApp : qaMixedSampleCacheApp.cc Psi2STo2K2PiGam : ;
diff --git a/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc b/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc
index 56d321199f72dca5b1626f5d76719f2a3273fd7f..9100b8bacdcd1812d9e4952443959a8cd8fcd3fe 100644
--- a/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc
+++ b/Examples/Psi2STo2K2PiGam/Mpsi2STo2K2PiGamTestApp.cc
@@ -184,6 +184,7 @@ int main(int __argc,char *__argv[]){
hypMap["Pi_2_2285Tof1710PiHyp9"]=true;
hypMap["Pi_2_2285ToK892KHyp9"]=true;
hypMap["Pi1800Pi0ToK892KHyp9"]=true;
+ hypMap["Pi_2_2285ToK_2_1430KHyp9"]=true;
hypMap["f980f_2_2300Hyp9"]=true;
hypMap["f_2_2300sigmaHyp9"]=true;
hypMap["K_2_1770ToK_2_1430PiHyp9"]=true;
diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc
index 36247168e4f3493b17e6e5da23dbccbb00169026..267c13d4f0ee16d52ac0eb0259c9faff48db4481 100644
--- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc
+++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamEvtList.cc
@@ -254,25 +254,23 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int
// calculate and store WignerD functions for KKpi0->KK and KKpi1-> KK
Spin jpi2=2;
- for (Spin M=-2; M<=2; M++){
for (Spin lamf2=-2; lamf2<=2; lamf2++){
- thePsi2STo2K2PiGamEvtData->DfPi2Tof2Pi0[jpi2][M][lamf2]=Wigner_D(KK_HeliKKPi0_4V.Phi(), KK_HeliKKPi0_4V.Theta(),0, jpi2,M,lamf2);
- thePsi2STo2K2PiGamEvtData->DfPi2Tof2Pi1[jpi2][M][lamf2]=Wigner_D(KK_HeliKKPi1_4V.Phi(), KK_HeliKKPi1_4V.Theta(),0, jpi2,M,lamf2);
+ thePsi2STo2K2PiGamEvtData->DfPi2Tof2Pi0[jpi2][0][lamf2]=Wigner_D(KK_HeliKKPi0_4V.Phi(), KK_HeliKKPi0_4V.Theta(),0, jpi2,0,lamf2);
+ thePsi2STo2K2PiGamEvtData->DfPi2Tof2Pi1[jpi2][0][lamf2]=Wigner_D(KK_HeliKKPi1_4V.Phi(), KK_HeliKKPi1_4V.Theta(),0, jpi2,0,lamf2);
}
- }
+
// calculate and store WignerD functions for KKpi0->KK and KKpi1-> KK
jpi2=2;
- for (Spin M=-2; M<=2; M++){
- for (Spin lamKstar=-1; lamKstar<=1; lamKstar++){
- thePsi2STo2K2PiGamEvtData->DfPi2ToKstarpK_pi0[jpi2][M][lamKstar]=Wigner_D(KpPi0_HeliKKPi0_4V.Phi(), KpPi0_HeliKKPi0_4V.Theta(),0, jpi2,M,lamKstar);
- thePsi2STo2K2PiGamEvtData->DfPi2ToKstarpK_pi1[jpi2][M][lamKstar]=Wigner_D(KpPi1_HeliKKPi1_4V.Phi(), KpPi1_HeliKKPi1_4V.Theta(),0, jpi2,M,lamKstar);
- thePsi2STo2K2PiGamEvtData->DfPi2ToKstarmK_pi0[jpi2][M][lamKstar]=Wigner_D(KmPi0_HeliKKPi0_4V.Phi(), KmPi0_HeliKKPi0_4V.Theta(),0, jpi2,M,lamKstar);
- thePsi2STo2K2PiGamEvtData->DfPi2ToKstarmK_pi1[jpi2][M][lamKstar]=Wigner_D(KmPi1_HeliKKPi1_4V.Phi(), KmPi1_HeliKKPi1_4V.Theta(),0, jpi2,M,lamKstar);
+ for (Spin lamKstar=-2; lamKstar<=2; lamKstar++){
+ thePsi2STo2K2PiGamEvtData->DfPi2ToKstarpK_pi0[jpi2][0][lamKstar]=Wigner_D(KpPi0_HeliKKPi0_4V.Phi(), KpPi0_HeliKKPi0_4V.Theta(),0, jpi2,0,lamKstar);
+ thePsi2STo2K2PiGamEvtData->DfPi2ToKstarpK_pi1[jpi2][0][lamKstar]=Wigner_D(KpPi1_HeliKKPi1_4V.Phi(), KpPi1_HeliKKPi1_4V.Theta(),0, jpi2,0,lamKstar);
+ thePsi2STo2K2PiGamEvtData->DfPi2ToKstarmK_pi0[jpi2][0][lamKstar]=Wigner_D(KmPi0_HeliKKPi0_4V.Phi(), KmPi0_HeliKKPi0_4V.Theta(),0, jpi2,0,lamKstar);
+ thePsi2STo2K2PiGamEvtData->DfPi2ToKstarmK_pi1[jpi2][0][lamKstar]=Wigner_D(KmPi1_HeliKKPi1_4V.Phi(), KmPi1_HeliKKPi1_4V.Theta(),0, jpi2,0,lamKstar);
}
- }
+
// 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++){
@@ -285,6 +283,20 @@ Psi2STo2K2PiGamEvtData* Psi2STo2K2PiGamEvtList::fillEvtData(Event* anEvent, int
}
}
+
+// 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);
+
+ }
+ }
+
+
for (Spin lam1=-1; lam1<=1; lam1++){
for (Spin lam2=-1; lam2<=1; lam2++){
diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc
index 9727df675ed283fd29dc5800d5d9f2f4774c1bb6..8ed3a6b8d842cc86fa279706570a5de8e6bc0a5a 100644
--- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc
+++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamHist.cc
@@ -192,9 +192,10 @@ 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:costhetapipi:phipipi:costhetapiViapipi:phipiViapipi:weight");
+ _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");
- _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:costhetapipi:phipipi:costhetapiViapipi:phipiViapipi:weight");
}
@@ -298,6 +299,8 @@ void Psi2STo2K2PiGamHist::fillHistos(const std::string& theKey, const Psi2STo2K2
void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEvtData* theData, double weight)
{
+ float datatype = 1.; // 1 = real data, 2 = pwamc
+
float Phichic0_HeliPsi2S = theData->chic0_HeliPsi2S_4V.Phi();
float CosThetachic0_HeliPsi2S = theData->chic0_HeliPsi2S_4V.CosTheta();
float PhiK1pi1pi2 = theData->KpPiPi_HeliChic0_4V.Phi();
@@ -320,10 +323,10 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv
float mK2pi1pi2 = theData->KmPiPi_HeliChic0_4V.M();
float mK1pi1 = theData->KpPi0_HeliChic0_4V.M();
float mK1pi2 = theData->KpPi1_HeliChic0_4V.M();
- float MK2pi1 = theData->KmPi0_HeliChic0_4V.M();
- float MK2pi2 = theData->KmPi1_HeliChic0_4V.M();
+ float mK2pi1 = theData->KmPi0_HeliChic0_4V.M();
+ float mK2pi2 = theData->KmPi1_HeliChic0_4V.M();
- float Mpipi = theData->PiPi_HeliChic0_4V.M();
+ float mpipi = theData->PiPi_HeliChic0_4V.M();
float CosThetaPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.CosTheta();
float PhiPiPiFromK1PiPi = theData->PiPi_HeliKpPi0Pi0_4V.Phi();
@@ -356,23 +359,23 @@ void Psi2STo2K2PiGamHist::writeNTuple(TNtuple* theTuple, const Psi2STo2K2PiGamEv
fourVectors[19] = mK2pi1pi2;
fourVectors[20] = mK1pi1;
fourVectors[21] = mK1pi2;
- fourVectors[22] = MK2pi1;
- fourVectors[23] = MK2pi2;
+ fourVectors[22] = mK2pi1;
+ fourVectors[23] = mK2pi2;
- fourVectors[24] = Mpipi;
+ fourVectors[24] = mpipi;
fourVectors[25] = CosThetaPiPiFromK1PiPi;
fourVectors[26] = PhiPiPiFromK1PiPi;
fourVectors[27] = CosThetaPiFromPiPi;
fourVectors[28] = PhiPiFromPiPi;
fourVectors[29] = evtweight;
+ fourVectors[30] = datatype;
// cout << evtweight << endl;
theTuple->Fill(fourVectors);
}
-
void Psi2STo2K2PiGamHist::plotCosPsi(const std::string& theKey, const Psi2STo2K2PiGamEvtData* theData, double weight){
std::string cosPsiString="cosPsi"+theKey;
TH1F* theHisto=_hist1DMap[cosPsiString];
diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc
index 07f93f547f4d0f3633e5d45a90fbde095f80975b..6b3117912349908abb4b51f5b23048358d29dc53 100644
--- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc
+++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.cc
@@ -71,6 +71,9 @@ Psi2STo2K2PiGamStates::Psi2STo2K2PiGamStates() :
//LS combinations for the pi_2 decay to K*1 K
fillJPCLS(_pi2JPC, _kJPC, _Kst1JPC,_Pi_2ToKst1KJPCLS);
+ //LS combinations for the pi_2 decay to K*2 K
+ fillJPCLS(_pi2JPC, _kJPC, _Kst2JPC,_Pi_2ToKst2KJPCLS);
+
//LS combinations for the chi_c0 to pi0 pi0
fillJPCLS(_chic0JPC, _pi0JPC, _pi0JPC, _ChiToPi0Pi0JPCLS);
@@ -87,6 +90,7 @@ Psi2STo2K2PiGamStates::Psi2STo2K2PiGamStates() :
fillJPCLS(_K1400JPC, _f0JPC, _kJPC, _K1pTof0KJPCLS);
fillJPCLS(_pi0JPC, _Kst1JPC, _pi0JPC, _Pi0pToKstKJPCLS);
+
}
Psi2STo2K2PiGamStates::~Psi2STo2K2PiGamStates()
@@ -221,6 +225,12 @@ void Psi2STo2K2PiGamStates::print(std::ostream& os) const
os << "\n" << std::endl;
}
+ os << "\n *** pi2: LS combinations for the pi2 decay to K*2 K *** "<< std::endl;
+ for ( itJPCLS=_Pi_2ToKst2KJPCLS.begin(); itJPCLS!=_Pi_2ToKst2KJPCLS.end(); ++itJPCLS){
+ (*itJPCLS)->print(os);
+ os << "\n" << std::endl;
+ }
+
os << "\n *** pi2: LS combinations for the chi_c0 decay to pi0 pi0 *** "<< std::endl;
for ( itJPCLS=_ChiToPi0Pi0JPCLS.begin(); itJPCLS!=_ChiToPi0Pi0JPCLS.end(); ++itJPCLS){
(*itJPCLS)->print(os);
diff --git a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh
index 3c1d6e47008763e1acdcc8190579160adaf36181..ef7cdee7bb048504830a6500d9544114a8c26210 100644
--- a/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh
+++ b/Examples/Psi2STo2K2PiGam/Psi2STo2K2PiGamStates.hh
@@ -39,6 +39,7 @@ public:
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;}
std::vector< boost::shared_ptr<const JPCLS> > Pi_2ToKst1KStates() const {return _Pi_2ToKst1KJPCLS;}
+ std::vector< boost::shared_ptr<const JPCLS> > Pi_2ToKst2KStates() const {return _Pi_2ToKst2KJPCLS;}
std::vector< boost::shared_ptr<const JPCLS> > ChiToPi0Pi0States() const {return _ChiToPi0Pi0JPCLS;}
std::vector< boost::shared_ptr<const JPCLS> > Pi_2Tof0PiStates() const {return _Pi_2Tof0PiJPCLS;}
std::vector< boost::shared_ptr<const JPCLS> > ChiToK2mK0mStates() const {return _ChiToK2mK0mJPCLS;}
@@ -80,6 +81,7 @@ private:
std::vector< boost::shared_ptr<const JPCLS> > _ChiToPi_2PiJPCLS;
std::vector< boost::shared_ptr<const JPCLS> > _Pi_2Tof_2PiJPCLS;
std::vector< boost::shared_ptr<const JPCLS> > _Pi_2ToKst1KJPCLS;
+ std::vector< boost::shared_ptr<const JPCLS> > _Pi_2ToKst2KJPCLS;
std::vector< boost::shared_ptr<const JPCLS> > _ChiToPi0Pi0JPCLS;
std::vector< boost::shared_ptr<const JPCLS> > _Pi_2Tof0PiJPCLS;
std::vector< boost::shared_ptr<const JPCLS> > _ChiToK2mK0mJPCLS;
diff --git a/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc b/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc
new file mode 100644
index 0000000000000000000000000000000000000000..dadecd82984b36e90d387f95718aa2bc48a31a80
--- /dev/null
+++ b/Examples/Psi2STo2K2PiGam/qaMixedSampleApp.cc
@@ -0,0 +1,767 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..b06397e225603b5620f0ebf3c2b0be4a8743ee2c
--- /dev/null
+++ b/Examples/Psi2STo2K2PiGam/qaMixedSampleCacheApp.cc
@@ -0,0 +1,726 @@
+#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<float> vecdatacosThetaKKpipi;
+std::vector<float> vecdataphiKKpipi;
+std::vector<float> vecdatacosThetaK1pi1pi2;
+std::vector<float> vecdataphiK1pi1pi2;
+std::vector<float> vecdatacosThetapi1pi2;
+std::vector<float> vecdataphipi1pi2;
+std::vector<float> vecdatacosThetapi1;
+std::vector<float> vecdataphipi1;
+std::vector<float> vecdatainvmassK1pi1pi2;
+std::vector<float> vecdatainvmasspi1pi2;
+
+std::vector<float> vecpwamccosThetaKKpipi;
+std::vector<float> vecpwamcphiKKpipi;
+std::vector<float> vecpwamccosThetaK1pi1pi2;
+std::vector<float> vecpwamcphiK1pi1pi2;
+std::vector<float> vecpwamccosThetapi1pi2;
+std::vector<float> vecpwamcphipi1pi2;
+std::vector<float> vecpwamccosThetapi1;
+std::vector<float> vecpwamcphipi1;
+std::vector<float> vecpwamcinvmassK1pi1pi2;
+std::vector<float> vecpwamcinvmasspi1pi2;
+
+std::vector<float> vecevtweight;
+
+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);
+
+
+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;
+
+ numOfEntriesData = ntdata->GetEntries();
+ cout << "Events in DATA Set: " << numOfEntriesData << endl;
+ numOfEntriesPwaMc = ntpwamc->GetEntries();
+ cout << "Events in PWAMC Set: " << numOfEntriesPwaMc << 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);
+
+ cout << "Read branches vom ntdata." << endl;
+
+ for(int i = 0; i<numOfEntriesData; i++) {
+ ntdata->GetEntry(i);
+ vecdatacosThetaKKpipi.push_back(datacosThetaKKpipi);
+ vecdataphiKKpipi.push_back(dataphiKKpipi);
+ vecdatacosThetaK1pi1pi2.push_back(datacosThetaK1pi1pi2);
+ vecdataphiK1pi1pi2.push_back(dataphiK1pi1pi2);
+ vecdatacosThetapi1pi2.push_back(datacosThetapi1pi2);
+ vecdataphipi1pi2.push_back(dataphipi1pi2);
+ vecdatacosThetapi1.push_back(datacosThetapi1);
+ vecdataphipi1.push_back(dataphipi1);
+ vecdatainvmassK1pi1pi2.push_back(datainvmassK1pi1pi2);
+ vecdatainvmasspi1pi2.push_back(datainvmasspi1pi2);
+ }
+ cout << "Data Events filled to vectors." << endl;
+
+ 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);
+
+ cout << "Read branches vom ntpwamc." << endl;
+
+ for(int i = 0; i<numOfEntriesPwaMc; i++) {
+ ntpwamc->GetEntry(i);
+ vecpwamccosThetaKKpipi.push_back(pwamccosThetaKKpipi);
+ vecpwamcphiKKpipi.push_back(pwamcphiKKpipi);
+ vecpwamccosThetaK1pi1pi2.push_back(pwamccosThetaK1pi1pi2);
+ vecpwamcphiK1pi1pi2.push_back(pwamcphiK1pi1pi2);
+ vecpwamccosThetapi1pi2.push_back(pwamccosThetapi1pi2);
+ vecpwamcphipi1pi2.push_back(pwamcphipi1pi2);
+ vecpwamccosThetapi1.push_back(pwamccosThetapi1);
+ vecpwamcphipi1.push_back(pwamcphipi1);
+ vecpwamcinvmassK1pi1pi2.push_back(pwamcinvmassK1pi1pi2);
+ vecpwamcinvmasspi1pi2.push_back(pwamcinvmasspi1pi2);
+ }
+ cout << "PwaMC Events filled to vectors." << 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);
+
+ 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 = 3;
+ cout << "Performing " << limit_runs_upper - limit_runs_lower << " runs." << endl;
+ std::vector<float> fpullvector;
+ unsigned int limit_data = 2874; // numOfEntriesData;
+ unsigned int limit_pwamc = 2874; // 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 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) 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);
+ 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;
+
+ float distance = sqrt(
+ pow((vecdatacosThetaKKpipi[dataentry1]-vecdatacosThetaKKpipi[dataentry2])/2.,2) +
+ pow((vecdataphiKKpipi[dataentry1]-vecdataphiKKpipi[dataentry2])/(2*3.1415),2) +
+ pow((vecdatacosThetaK1pi1pi2[dataentry1]-vecdatacosThetaK1pi1pi2[dataentry2])/2. ,2) +
+ pow((vecdataphiK1pi1pi2[dataentry1]-vecdataphiK1pi1pi2[dataentry2])/(2*3.1415) ,2) +
+ pow((vecdatacosThetapi1pi2[dataentry1]-vecdatacosThetapi1pi2[dataentry2])/2. ,2) +
+ pow((vecdataphipi1pi2[dataentry1]-vecdataphipi1pi2[dataentry2])/(2*3.1415) ,2) +
+ pow((vecdatacosThetapi1[dataentry1]-vecdatacosThetapi1[dataentry2]) ,2) +
+ pow((vecdataphipi1[dataentry1]-vecdataphipi1[dataentry2])/(3.1415) ,2) +
+ pow((vecdatainvmassK1pi1pi2[dataentry1]-vecdatainvmassK1pi1pi2[dataentry2] )/2.2 ,2) +
+ pow((vecdatainvmasspi1pi2[dataentry1]-vecdatainvmasspi1pi2[dataentry2])/2.3 ,2)
+ );
+
+ return distance;
+
+}
+
+
+
+float GetDistanceDataToPwaMc(int dataentry, int pwamcentry) {
+ // cout << "Distance from Data Event " << dataentry << " to PwaMc Event " << pwamcentry << endl;
+
+ float distance = sqrt(
+ pow((vecdatacosThetaKKpipi[dataentry]-vecpwamccosThetaKKpipi[pwamcentry])/2.,2) +
+ pow((vecdataphiKKpipi[dataentry]-vecpwamcphiKKpipi[pwamcentry])/(2*3.1415),2) +
+ pow((vecdatacosThetaK1pi1pi2[dataentry]-vecpwamccosThetaK1pi1pi2[pwamcentry])/2. ,2) +
+ pow((vecdataphiK1pi1pi2[dataentry]-vecpwamcphiK1pi1pi2[pwamcentry])/(2*3.1415) ,2) +
+ pow((vecdatacosThetapi1pi2[dataentry]-vecpwamccosThetapi1pi2[pwamcentry])/2. ,2) +
+ pow((vecdataphipi1pi2[dataentry]-vecpwamcphipi1pi2[pwamcentry])/(2*3.1415) ,2) +
+ pow((vecdatacosThetapi1[dataentry]-vecpwamccosThetapi1[pwamcentry]) ,2.) +
+ pow((vecdataphipi1[dataentry]-vecpwamcphipi1[pwamcentry])/(3.1415) ,2) +
+ pow((vecdatainvmassK1pi1pi2[dataentry]-vecpwamcinvmassK1pi1pi2[pwamcentry])/2.2 ,2) +
+ pow((vecdatainvmasspi1pi2[dataentry]-vecpwamcinvmasspi1pi2[pwamcentry])/2.3 ,2)
+ );
+
+ return distance;
+
+}
+
+
+
+
+
+
+
+
+float GetDistancePwaMcToPwaMc(int pwamcentry1, int pwamcentry2) {
+ // cout << "Distance from PwaMC Event " << dataentry1 << " to PwaMC Event " << dataentry2 << endl;
+
+ float distance = sqrt(
+ pow((vecpwamccosThetaKKpipi[pwamcentry1]-vecpwamccosThetaKKpipi[pwamcentry2])/2.,2) +
+ pow((vecpwamcphiKKpipi[pwamcentry1]-vecpwamcphiKKpipi[pwamcentry2])/(2*3.1415),2) +
+ pow((vecpwamccosThetaK1pi1pi2[pwamcentry1]-vecpwamccosThetaK1pi1pi2[pwamcentry2])/2. ,2) +
+ pow((vecpwamcphiK1pi1pi2[pwamcentry1]-vecpwamcphiK1pi1pi2[pwamcentry2])/(2*3.1415) ,2) +
+ pow((vecpwamccosThetapi1pi2[pwamcentry1]-vecpwamccosThetapi1pi2[pwamcentry2])/2. ,2) +
+ pow((vecpwamcphipi1pi2[pwamcentry1]-vecpwamcphipi1pi2[pwamcentry2])/(2*3.1415) ,2) +
+ pow((vecpwamccosThetapi1[pwamcentry1]-vecpwamccosThetapi1[pwamcentry2]) ,2) +
+ pow((vecpwamcphipi1[pwamcentry1]-vecpwamcphipi1[pwamcentry2])/(3.1415) ,2) +
+ pow((vecpwamcinvmassK1pi1pi2[pwamcentry1]-vecpwamcinvmassK1pi1pi2[pwamcentry2])/2.2 ,2) +
+ pow((vecpwamcinvmasspi1pi2[pwamcentry1]-vecpwamcinvmasspi1pi2[pwamcentry2])/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/PspGen/PspGenTestApp.cc b/PspGen/PspGenTestApp.cc
index 297a9bb3f6d42c68cb29c969ae1f3fc137f67214..64358f0cff02b0c8415951f2419db4f585e25ed3 100644
--- a/PspGen/PspGenTestApp.cc
+++ b/PspGen/PspGenTestApp.cc
@@ -19,133 +19,159 @@
#include <iostream>
#include <string>
+#include <sstream>
using namespace std;
int main(int argc, char* argv[])
{
- HepMC::IO_GenEvent hepMC_out("phspEvents_1mio_RndSeed4720.out",std::ios::out);
-
- EvtSimpleRandomEngine myRandom(4720);
- EvtRandom::setRandomEngine(&myRandom);
-
- // variables for the production e+ e- Psi2s
- double mass_e=0.511e-3;
- double mass_psi2s = 3.68609;
- double pe_tot=0.5*sqrt(mass_psi2s*mass_psi2s-2.*mass_e);
-
- // variables for the first decay Psi2s->Chi_c0+gamma
- int firstNdaug = 2;
- double firstMass[8];
- EvtVector4R firstP4[8];
- double firstMp = 3.68609;
-
- firstMass[0] = 3.41475;
- firstMass[10] = 0.;
-
- // variables for the second decay Chi->c0->K+ K- pi0 pi0
- int ndaug = 4;
- double mass[30];
- EvtVector4R p4[30];
-
- mass[0] = .493677;
- mass[1] = .493677;
- mass[2] = .1349766;
- mass[3] = .1349766;
-
- // 5 events for visual control
- for (int count = 0; count < 5; count++) {
- std::cout << "Event #" << count << std::endl;
- EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
- for (int t = 0; t < firstNdaug; t++) {
- std::cout << firstP4[t] << " m = " << sqrt(firstP4[t]*firstP4[t]) << std::endl;
+
+ int set_limit = 200;
+ int eventsperset = 500000;
+ bool silent = true;
+
+ for(int set=100; set < set_limit; set++) {
+
+ std::stringstream out;
+ out << set+1001;
+ string setstring = out.str();
+
+ string filename = "phspEvents_500k_set" + setstring + ".out";
+ cout << "Filename: " << filename << endl;
+ cout << endl;
+
+ HepMC::IO_GenEvent hepMC_out(filename,std::ios::out);
+
+ EvtSimpleRandomEngine myRandom(4701+set);
+ EvtRandom::setRandomEngine(&myRandom);
+
+ // variables for the production e+ e- Psi2s
+ double mass_e=0.511e-3;
+ double mass_psi2s = 3.68609;
+ double pe_tot=0.5*sqrt(mass_psi2s*mass_psi2s-2.*mass_e);
+
+ // variables for the first decay Psi2s->Chi_c0+gamma
+ int firstNdaug = 2;
+ double firstMass[8];
+ EvtVector4R firstP4[8];
+ double firstMp = 3.68609;
+
+ firstMass[0] = 3.41475;
+ firstMass[10] = 0.;
+
+ // variables for the second decay Chi->c0->K+ K- pi0 pi0
+ int ndaug = 4;
+ double mass[30];
+ EvtVector4R p4[30];
+
+ mass[0] = .493677;
+ mass[1] = .493677;
+ mass[2] = .1349766;
+ mass[3] = .1349766;
+
+ // 5 events for visual control
+ if(!silent) {
+ for (int count = 0; count < 5; count++) {
+ std::cout << "Event #" << count << std::endl;
+ EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
+ for (int t = 0; t < firstNdaug; t++) {
+ std::cout << firstP4[t] << " m = " << sqrt(firstP4[t]*firstP4[t]) << std::endl;
+ }
+ // now use the mass of the Chi_c0 from the first decay as input for the second decay
+ EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
+ for (int t = 0; t < ndaug; t++) {
+ std::cout << p4[t] << " m = " << sqrt(p4[t]*p4[t]) << std::endl;
+ p4[t].applyBoostTo(firstP4[0]);
+ std::cout << p4[t] << " boostet: m = " << sqrt(p4[t]*p4[t]) << std::endl;
+ }
+ }
}
- // now use the mass of the Chi_c0 from the first decay as input for the second decay
- EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
- for (int t = 0; t < ndaug; t++) {
- std::cout << p4[t] << " m = " << sqrt(p4[t]*p4[t]) << std::endl;
- p4[t].applyBoostTo(firstP4[0]);
- std::cout << p4[t] << " boostet: m = " << sqrt(p4[t]*p4[t]) << std::endl;
- }
- }
-
- TApplication* rootapp = new TApplication("example",&argc, argv);
- TH2F dalitz("Dalitz plot", "Dalitz Plot", 80,0.,5., 80,0.,6.);
-
- TH1F massKpi("invMassKpi", "invMassKpi", 512, 0., 3.);
- TH1F massKpipi("invMassKpipi", "invMassKpipi", 512, 0., 3.);
- TH1F cosThetapi("cosThetapi", "cosThetapi", 512, -1.2, 1.2);
- TH1F cosThetaKpi("cosThetaKpi", "cosThetaKpi", 512, -1.2, 1.2);
-
- for (int count = 0; count < 1000000; count++) {
- if(count % 100000 == 0) cout << "Event " << count << " generated." << endl;
-
- EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
- EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
- for (int t = 0; t < ndaug; t++) {
- p4[t].applyBoostTo(firstP4[0]);
+ /*
+ TApplication* rootapp = new TApplication("example",&argc, argv);
+ TH2F dalitz("Dalitz plot", "Dalitz Plot", 80,0.,5., 80,0.,6.);
+
+ TH1F massKpi("invMassKpi", "invMassKpi", 512, 0., 3.);
+ TH1F massKpipi("invMassKpipi", "invMassKpipi", 512, 0., 3.);
+ TH1F cosThetapi("cosThetapi", "cosThetapi", 512, -1.2, 1.2);
+ TH1F cosThetaKpi("cosThetaKpi", "cosThetaKpi", 512, -1.2, 1.2);
+ */
+
+ for (int count = 0; count < eventsperset; count++) {
+ if(count % 100000 == 0) cout << "Event " << count << " generated." << endl;
+
+ EvtGenKine::PhaseSpace(firstNdaug, firstMass, firstP4, firstMp);
+ EvtGenKine::PhaseSpace(ndaug, mass, p4, firstP4[0].mass());
+ for (int t = 0; t < ndaug; t++) {
+ p4[t].applyBoostTo(firstP4[0]);
+ }
+
+ // dalitz.Fill((p4[0]+p4[1])*(p4[0]+p4[1]), (p4[1]+p4[2])*(p4[1]+p4[2]));
+
+ HepMC::GenEvent* evt = new HepMC::GenEvent( 20, count );
+ evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
+
+ // create production vertex
+ HepMC::GenVertex* vtx_prod = new HepMC::GenVertex();
+ evt->add_vertex( vtx_prod );
+ HepMC::GenParticle* eplus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), -11, 3 );
+ HepMC::GenParticle* eminus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., -pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), 11, 3 );
+ vtx_prod->add_particle_in( eplus_part );
+ vtx_prod->add_particle_in( eminus_part );
+
+ HepMC::GenParticle* psi2s_part=new HepMC::GenParticle( HepMC::FourVector(0., 0., 0., mass_psi2s), 20443, 999);
+ vtx_prod->add_particle_out(psi2s_part);
+
+ // create Psi(2S) vertex
+ HepMC::GenVertex* vtx_psi2S = new HepMC::GenVertex();
+ evt->add_vertex( vtx_psi2S );
+ vtx_psi2S->add_particle_in( psi2s_part );
+ HepMC::GenParticle* chi_c0_part = new HepMC::GenParticle(HepMC::FourVector(firstP4[0].get(1), firstP4[0].get(2), firstP4[0].get(3), firstP4[0].get(0)), 10441, 999 );
+ vtx_psi2S->add_particle_out(chi_c0_part);
+ vtx_psi2S->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(firstP4[1].get(1), firstP4[1].get(2), firstP4[1].get(3), firstP4[1].get(0)), 22, 1 ) );
+
+
+ // now do the chi_c0 vertex
+ HepMC::GenVertex* vtx_chi_c0 = new HepMC::GenVertex();
+ evt->add_vertex( vtx_chi_c0 );
+ vtx_chi_c0->add_particle_in( chi_c0_part );
+ vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[0].get(1),p4[0].get(2),p4[0].get(3), p4[0].get(0)), 321, 1 ));
+ vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[1].get(1),p4[1].get(2),p4[1].get(3), p4[1].get(0)), -321, 1 ) );
+ vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[2].get(1),p4[2].get(2),p4[2].get(3), p4[2].get(0)), 111, 1 ) );
+ vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[3].get(1),p4[3].get(2),p4[3].get(3), p4[3].get(0)), 111, 1 ) );
+
+ // evt->print(std::cout);
+ hepMC_out << evt;
+ delete evt;
+
+ /*
+ massKpi.Fill( (p4[0]+p4[2]).mass(), 1.);
+ massKpi.Fill( (p4[0]+p4[3]).mass(), 1.);
+ massKpi.Fill( (p4[1]+p4[2]).mass(), 1.);
+ massKpi.Fill( (p4[1]+p4[3]).mass(), 1.);
+
+ massKpipi.Fill( (p4[0]+p4[2]+p4[3]).mass(), 1.);
+ massKpipi.Fill( (p4[1]+p4[2]+p4[3]).mass(), 1.);
+
+ cosThetapi.Fill( p4[2].get(3) /
+ sqrt(p4[2].get(1)*p4[2].get(1) + p4[2].get(2)*p4[2].get(2) + p4[2].get(3)*p4[2].get(3)),
+ 1.);
+
+ cosThetaKpi.Fill( (p4[0]+p4[2]).get(3) /
+ sqrt((p4[0]+p4[2]).get(1)*(p4[0]+p4[2]).get(1) +
+ (p4[0]+p4[2]).get(2)*(p4[0]+p4[2]).get(2) +
+ (p4[0]+p4[2]).get(3)*(p4[0]+p4[2]).get(3)),
+ 1.);
+ */
}
- // dalitz.Fill((p4[0]+p4[1])*(p4[0]+p4[1]), (p4[1]+p4[2])*(p4[1]+p4[2]));
-
- HepMC::GenEvent* evt = new HepMC::GenEvent( 20, count );
- evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
-
- // create production vertex
- HepMC::GenVertex* vtx_prod = new HepMC::GenVertex();
- evt->add_vertex( vtx_prod );
- HepMC::GenParticle* eplus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), -11, 3 );
- HepMC::GenParticle* eminus_part = new HepMC::GenParticle( HepMC::FourVector(0., 0., -pe_tot, sqrt(mass_e*mass_e+pe_tot*pe_tot)), 11, 3 );
- vtx_prod->add_particle_in( eplus_part );
- vtx_prod->add_particle_in( eminus_part );
-
- HepMC::GenParticle* psi2s_part=new HepMC::GenParticle( HepMC::FourVector(0., 0., 0., mass_psi2s), 20443, 999);
- vtx_prod->add_particle_out(psi2s_part);
-
- // create Psi(2S) vertex
- HepMC::GenVertex* vtx_psi2S = new HepMC::GenVertex();
- evt->add_vertex( vtx_psi2S );
- vtx_psi2S->add_particle_in( psi2s_part );
- HepMC::GenParticle* chi_c0_part = new HepMC::GenParticle(HepMC::FourVector(firstP4[0].get(1), firstP4[0].get(2), firstP4[0].get(3), firstP4[0].get(0)), 10441, 999 );
- vtx_psi2S->add_particle_out(chi_c0_part);
- vtx_psi2S->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(firstP4[1].get(1), firstP4[1].get(2), firstP4[1].get(3), firstP4[1].get(0)), 22, 1 ) );
-
-
- // now do the chi_c0 vertex
- HepMC::GenVertex* vtx_chi_c0 = new HepMC::GenVertex();
- evt->add_vertex( vtx_chi_c0 );
- vtx_chi_c0->add_particle_in( chi_c0_part );
- vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[0].get(1),p4[0].get(2),p4[0].get(3), p4[0].get(0)), 321, 1 ));
- vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[1].get(1),p4[1].get(2),p4[1].get(3), p4[1].get(0)), -321, 1 ) );
- vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[2].get(1),p4[2].get(2),p4[2].get(3), p4[2].get(0)), 111, 1 ) );
- vtx_chi_c0->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(p4[3].get(1),p4[3].get(2),p4[3].get(3), p4[3].get(0)), 111, 1 ) );
-
- // evt->print(std::cout);
- hepMC_out << evt;
- delete evt;
-
- massKpi.Fill( (p4[0]+p4[2]).mass(), 1.);
- massKpi.Fill( (p4[0]+p4[3]).mass(), 1.);
- massKpi.Fill( (p4[1]+p4[2]).mass(), 1.);
- massKpi.Fill( (p4[1]+p4[3]).mass(), 1.);
-
- massKpipi.Fill( (p4[0]+p4[2]+p4[3]).mass(), 1.);
- massKpipi.Fill( (p4[1]+p4[2]+p4[3]).mass(), 1.);
-
- cosThetapi.Fill( p4[2].get(3) /
- sqrt(p4[2].get(1)*p4[2].get(1) + p4[2].get(2)*p4[2].get(2) + p4[2].get(3)*p4[2].get(3)),
- 1.);
-
- cosThetaKpi.Fill( (p4[0]+p4[2]).get(3) /
- sqrt((p4[0]+p4[2]).get(1)*(p4[0]+p4[2]).get(1) +
- (p4[0]+p4[2]).get(2)*(p4[0]+p4[2]).get(2) +
- (p4[0]+p4[2]).get(3)*(p4[0]+p4[2]).get(3)),
- 1.);
+ cout << "Generated PHSP set " << set+1 << ", wrote to " << filename << endl;
+
}
-
+
+ cout << endl;
cout << "Finished!" << endl;
-// massKpipi.Draw();
-// rootapp->Run();
-
+ // massKpipi.Draw();
+ // rootapp->Run();
+
return 0;
}