Loading...
Searching...
No Matches
QCMCFilter Class Reference
#include <QCMCFilter.h>
Inheritance diagram for QCMCFilter:
Public Member Functions | |
| QCMCFilter (const std::string &name, ISvcLocator *pSvcLocator) | |
| std::vector< double > | findD0Decay (int charm) |
| StatusCode | initialize () |
| StatusCode | execute () |
| StatusCode | finalize () |
Detailed Description
Definition at line 17 of file QCMCFilter.h.
Constructor & Destructor Documentation
◆ QCMCFilter()
| QCMCFilter::QCMCFilter | ( | const std::string & | name, |
| ISvcLocator * | pSvcLocator | ||
| ) |
Definition at line 235 of file QCMCFilter.cxx.
235 :
236 Algorithm(name, pSvcLocator){
237 //Declare the properties
238 declareProperty("Debug", m_debug= true );
239 declareProperty("x", m_x= 0.00407);
240 declareProperty("y", m_y= 0.00647); //For values outside of (-.11 , .06) measured y value may vary from input y
241
242 declareProperty("rForCFModes", m_rCF= 0.0621); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
243 declareProperty("zForCFModes", m_zCF= 2.0);
244 declareProperty("wSignForCFModes", m_wCFSign= true);
245 declareProperty("rForCFMode_0", m_rCF_0= 0.0587); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
246 declareProperty("RForCFMode_0", m_RCF_0= 1.); //should be kept near this.
247 declareProperty("zForCFMode_0", m_zCF_0= 1.9585);
248 declareProperty("dForCFMode_0", m_dCF_0= 191.7*3.14159/180.);
249 declareProperty("wSignForCFMode_0", m_wCFSign_0 = true);
250 declareProperty("rForCFModes_1", m_rCF_1= 0.044); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
251 declareProperty("RForCFMode_1", m_RCF_1= 0.79); //should be kept near this.
252 declareProperty("zForCFModes_1", m_zCF_1= 1.9225);
253 declareProperty("dForCFMode_1", m_dCF_1= 196.*3.14159/180.);
254 declareProperty("wSignForCFModes_1", m_wCFSign_1= true);
255 declareProperty("rForCFModes_3", m_rCF_3= 0.0549); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
256 declareProperty("RForCFMode_3", m_RCF_3= 0.43); //should be kept near this.
257 declareProperty("zForCFModes_3", m_zCF_3= 1.2313);
258 declareProperty("dForCFMode_3", m_dCF_3= 161.*3.14159/180.);
259 declareProperty("wSignForCFModes_3", m_wCFSign_3= true);
260
261 declareProperty("rForCSModes", m_rCS= 1.0);
262 declareProperty("zForCSModes", m_zCS= -2.0);
263 declareProperty("wSignForCSModes", m_wCSSign= true);
264 declareProperty("DplusDminusWeight", m_dplusDminusWeight= 0.5);
265 declareProperty("dalitzModel", m_dalitzModel= 0); // 0 CLEO-c, 1 Babar, 2 Belle
266 declareProperty("UseNewWeights", m_useNewWeights= true);
267 declareProperty("InvertSelection", m_invertSelection= false);
268 }
Member Function Documentation
◆ execute()
| StatusCode QCMCFilter::execute | ( | ) |
Definition at line 826 of file QCMCFilter.cxx.
826 {
827
828 //interface to event data service
829 ISvcLocator* svcLocator = Gaudi::svcLocator();
830 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
831 if (sc.isFailure())
832 std::cout<<"Could not access EventDataSvc!"<<std::endl;
833
834 //setFilterPassed(false);
835
836 MsgStream log(msgSvc(), name());
837 log<< MSG::INFO << "in execute()" << endmsg;
838
839 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
840
841 if(!eventHeader)
842 {
843 cout<<" eventHeader "<<endl;
844 return StatusCode::FAILURE;
845 }
847 long event = eventHeader->eventNumber();
848 log << MSG::DEBUG << "run, evtnum = "
850 << event << endmsg;
851
852 unsigned int QC_status = 0;
853 //eventHeader->setEventTag(0);
854
855 // Get McParticle collection
856 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
857
858 //Check if event is Psi(3770)
859 int psipp_daughter = 0;
860 int d0_count = 0;
861 int chd_count = 0;
862 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
863 {
864 int pdg_code = (*it)->particleProperty();
865 if ((*it)->primaryParticle()) continue;
867 int mother_pdg = 0;
868 //finds if the particle is daughter of Psi(3770)
869 mother_pdg = it2.particleProperty();
870 if (mother_pdg == kPsi3770ID)
871 {
872 psipp_daughter++; //Found daughter
876 }
877 }
878
879 if( psipp_daughter == 0 ) //didn't find Psi(3770)
880 {
881 ++m_nUnknownEvents ;
882
883 if( m_invertSelection )
884 {
885 //eventHeader->setEventTag(0);
886 QC_status = 0;
887 eventHeader->setEventTag(QC_status);
888
889 log << MSG::INFO << "QC_status " << QC_status << endmsg;
890 log<< MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg ;
891 return StatusCode::SUCCESS; //discard event
892 }
893 else
894 {
895 // -------- Write to root file
896 //setFilterPassed(true);
897 //eventHeader->setEventTag(1);
898 QC_status = 1;
899 eventHeader->setEventTag(QC_status);
900
901 log << MSG::INFO << "QC_status " << QC_status << endmsg;
902
903 return StatusCode::SUCCESS; //kept event
904 }
905 }
906
907 log<< MSG::DEBUG << "Found psi(3770) -->" << endmsg ;
908
909 // Check that top particle has only two daughters + possible FSR photon
910 if ( psipp_daughter > 3 )
911 {
912 ++m_nUnknownEvents ;
913
914 if( m_invertSelection )
915 {
916 //eventHeader->setEventTag(0);
917 QC_status = 0;
918 eventHeader->setEventTag(QC_status);
919
920 log << MSG::INFO << "QC_status " << QC_status << endmsg;
921 log<< MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg ;
922 return StatusCode::SUCCESS; //discard event
923 }
924 else
925 {
926 // -------- Write to root file
927 //eventHeader->setEventTag(1);
928 QC_status = 1;
929 eventHeader->setEventTag(QC_status);
930
931 log << MSG::INFO << "QC_status " << QC_status << endmsg;
932 //setFilterPassed(true);
933
934
935 return StatusCode::SUCCESS; //kept event
936 }
937 }
938
939 //Check if D+D- state
940 if (chd_count == 2) { // D+D- event
941 ++m_nDpDmEvents ;
942
943 // D+D- must be rewieghted otherwise too many D+D- events relative to D0D0bar.
944 double random = RandFlat::shoot() ;
945
946 if( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
947 ( random > m_dplusDminusWeight && m_invertSelection ) )
948 {
949 // -------- Write to root file
950 //setFilterPassed(true);
951 //eventHeader->setEventTag(1);
952 QC_status = 1;
953 eventHeader->setEventTag(QC_status);
954
955 log << MSG::INFO << "QC_status " << QC_status << endmsg;
956
957
958 return StatusCode::SUCCESS; //event kept
959 }
960 else
961 {
962 //eventHeader->setEventTag(0);
963 QC_status = 0;
964 eventHeader->setEventTag(QC_status);
965
966 log << MSG::INFO << "QC_status " << QC_status << endmsg;
967 log<< MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg ;
968 ++m_nDpDmDiscarded ;
969 return StatusCode::SUCCESS; //event discarded
970 }
971 }
972
973 //Check if D0D0Bar
974 if (d0_count != 2) {
975 if( !m_invertSelection )
976 {
977 //eventHeader->setEventTag(0);
978 QC_status = 0;
979 eventHeader->setEventTag(QC_status);
980
981 log << MSG::INFO << "QC_status " << QC_status << endmsg;
982 log<< MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg ;
983 return StatusCode::SUCCESS; //discard event
984 }
985 else
986 {
987 // -------- Write to root file
988 //eventHeader->setEventTag(1);
989 QC_status = 1;
990 eventHeader->setEventTag(QC_status);
991
992 log << MSG::INFO << "QC_status " << QC_status << endmsg;
993 //setFilterPassed(true);
994
995
996 return StatusCode::SUCCESS; //kept event
997 }
998 }
999
1000 log<< MSG::INFO << "D0-D0bar event." << endmsg ;
1001 m_nD0D0barEvents++ ;
1002
1003 m_rho_flag = 0;
1004
1005 //Get D0-D0bar modes and info
1006 vector<double> d0_decay = findD0Decay(1);
1007 vector<double> d0bar_decay = findD0Decay(-1);
1008
1009 int d0mode = d0_decay[0];
1010 int d0barmode = d0bar_decay[0];
1011 if(m_debug)cout<<"D0 "<<d0mode<<endl;
1012 if(m_debug)cout<<"D0bar "<<d0barmode<<endl;
1013 m_modeCounter[d0mode][d0barmode]++;
1014
1015
1016
1017 //if any event unidentified remove
1018 if (d0_decay[0] == 18 || d0bar_decay[0] == 18 ) {
1019 if( !m_invertSelection )
1020 {
1021 //eventHeader->setEventTag(0);
1022 QC_status = 0;
1023 eventHeader->setEventTag(QC_status);
1024
1025 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1026 log<< MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg ;
1027 return StatusCode::SUCCESS; //discard event
1028 }
1029 else
1030 {
1031 // -------- Write to root file
1032 //eventHeader->setEventTag(1);
1033 QC_status = 1;
1034 eventHeader->setEventTag(QC_status);
1035
1036 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1037 //setFilterPassed(true);
1038
1039
1040
1041 return StatusCode::SUCCESS; //kept event
1042 }
1043 }
1044
1045 // Loop over D's
1046 double r2D0 = d0_decay[1] ;
1047 double deltaD0 = d0_decay[2] ;
1048 double RD0 = d0_decay[3] ;
1049 double FCP = d0_decay[4] ;
1050 double r2D0bar = d0bar_decay[1] ;
1051 double deltaD0bar = d0bar_decay[2] ;
1052 double RD0bar = d0bar_decay[3] ;
1053 double FCPbar = d0bar_decay[4] ;
1054
1055
1056 //Weight based on DDbar combination
1058 if( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz || d0_decay[0] == k2Pip2Pim || d0bar_decay[0] == k2Pip2Pim || d0_decay[0] == kPipPim2Pi0 || d0bar_decay[0] == kPipPim2Pi0 || d0_decay[0] == kK0PiPiPi0 || d0bar_decay[0] == kK0PiPiPi0
1059 || d0_decay[0] == kKKPiPi || d0bar_decay[0] == kKKPiPi ||d0bar_decay[0] == kFlavoredCF_3 || d0bar_decay[0] == kFlavoredCFBar_3 )
1060 {
1061 double r2prod = r2D0 * r2D0bar ;
1064
1065 //double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
1066 //double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
1067 //double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
1068 //double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
1069
1070 //weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
1071 //weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
1072 double bD0 = 1. - (2*FCP - 1.) * RD0 * sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
1073 double rwsD0 = ( r2D0 - (2*FCP - 1.) * RD0 * sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
1074 double bD0bar = 1. - (2*FCPbar - 1.) * RD0bar * sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
1075 double rwsD0bar = ( r2D0bar - (2*FCPbar - 1.) * RD0bar * sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
1076
1077 log << MSG::INFO << "bD0 == " << bD0<< endmsg;
1078 log << MSG::INFO << "bD0bar == " << bD0bar<< endmsg;
1079 log << MSG::INFO << "rwsD0 == " << rwsD0<< endmsg;
1080 log << MSG::INFO << "rwsD0bar == " << rwsD0bar<< endmsg;
1081 weight = 1. + r2prod - 2. * RD0 * (2*FCP - 1.) * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
1083 weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
1085 weight /= m_largestWeight ;
1087
1088 //get dalitz graph information before cut
1089 int k = -10;
1090 double m2i= 0;
1091 double m2j= 0;
1092 double m2k= 0;
1093 if( d0_decay[0] == kDalitz) {
1094 k = d0bar_decay[0];
1095 m2i= d0_decay[4];
1096 m2j= d0_decay[5];
1097 m2k= d0_decay[7];}
1098 if( d0bar_decay[0] == kDalitz) {
1099 k = d0_decay[0];
1100 m2i= d0bar_decay[4];
1101 m2j= d0bar_decay[5];
1102 m2k= d0bar_decay[7];}
1103
1104 }
1105 else
1106 {
1108 }
1109
1110 double random = RandFlat::shoot() ;
1111 log<< MSG::DEBUG << "type: " << d0_decay[0] << " " << d0bar_decay[0] << ", weight " <<
1113
1115 {
1116 // -------- Write to root file
1117 QC_status = 1;
1118 eventHeader->setEventTag(QC_status);
1119
1120 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1121 //eventHeader->setEventTag(1);
1122 //setFilterPassed(true);
1123
1124 //get dalitz graph information after filter
1125 /*int k = -10;
1126 double m2i= 0;
1127 double m2j= 0;
1128 double m2k= 0;
1129 double cosd = cos( deltaD0 ) ;
1130 double sind = sin( deltaD0 ) ;
1131 if( d0_decay[0] == kDalitz) {
1132 k = d0bar_decay[0];
1133 m2i= d0_decay[4];
1134 m2j= d0_decay[5];
1135 m2k= d0_decay[7];
1136 if ( m2j - m2i > 0. ) { //avoids doublecounting the branching fraction
1137 dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
1138 dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
1139 dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
1140 }
1141 if( d0bar_decay[0] == kDalitz) {
1142 k = d0_decay[0];
1143 m2i= d0bar_decay[4];
1144 m2j= d0bar_decay[5];
1145 m2k= d0bar_decay[7];
1146 cosd = cos( deltaD0bar ) ;
1147 sind = sin( deltaD0bar ) ;
1148 if( m2j - m2i < 0. ) { //avoids doublecounting the branching fraction
1149 dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
1150 dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
1151 dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
1152 }
1153
1154 m_keptModeCounter[d0mode][d0barmode]++;*/
1155 return StatusCode::SUCCESS; //event kept
1156 }
1157 else
1158 {
1159 QC_status = 0;
1160 eventHeader->setEventTag(QC_status);
1161
1162 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1163 //eventHeader->setEventTag(0);
1164 log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg ;
1165
1166 ++m_nD0D0barDiscarded ;
1167 return StatusCode::SUCCESS; //event discarded
1168 }
1169
1170}
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
HepMatrix m_modeCounter(19, 19, 0)
HepSymMatrix m_weights(18, 0)
Definition: McParticle.h:28
◆ finalize()
| StatusCode QCMCFilter::finalize | ( | ) |
Definition at line 1174 of file QCMCFilter.cxx.
1174 {
1175 MsgStream log(msgSvc(), name());
1176 log << MSG::INFO << "in finalize()" << endmsg;
1181 {
1182 log << MSG::INFO <<"Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents ) << endmsg ;
1183 }
1184
1185 log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar << endmsg ;
1187 /*if( m_nD0D0barEvents > 0 && m_nCPPlus > 0 && m_nCPMinus > 0)
1188 {
1189 log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double( m_nD0D0barEvents ) << endl <<
1190 "Fraction unknown decays: " << 0.5 * double( m_nUnknownDecays ) / double( m_nD0D0barEvents ) <<endmsg ;
1191
1192 double nD0 = 2. * m_nD0D0barEvents ;
1193 double brSL = m_nSL / nD0 ;
1194 double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
1195 double brCF_0 = m_nFlavoredCFD0_0 / nD0 ;
1196 double dBrCF_0 = sqrt( brCF_0 * ( 1. - brCF_0 ) / nD0 ) ;
1197 double brCF_1 = m_nFlavoredCFD0_1 / nD0 ;
1198 double dBrCF_1 = sqrt( brCF_1 * ( 1. - brCF_1 ) / nD0 ) ;
1199 double brCF_3 = m_nFlavoredCFD0_3 / nD0 ;
1200 double dBrCF_3 = sqrt( brCF_3 * ( 1. - brCF_3 ) / nD0 ) ;
1201 double brCS = m_nFlavoredCSD0 / nD0 ;
1202 double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
1203 double brDCS_0 = m_nFlavoredDCSD0_0 / nD0 ;
1204 double dBrDCS_0 = sqrt( brDCS_0 * ( 1. - brDCS_0 ) / nD0 ) ;
1205 double brDCS_1 = m_nFlavoredDCSD0_1 / nD0 ;
1206 double dBrDCS_1 = sqrt( brDCS_1 * ( 1. - brDCS_1 ) / nD0 ) ;
1207 double brDCS_3 = m_nFlavoredDCSD0_3 / nD0 ;
1208 double dBrDCS_3 = sqrt( brDCS_3 * ( 1. - brDCS_3 ) / nD0 ) ;
1209 double brCPPlus = m_nCPPlus / nD0 ;
1210 double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
1211 double brCPMinus = m_nCPMinus / nD0 ;
1212 double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
1213 double brDalitz = m_nDalitz / nD0 ;
1214 double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
1215 double fdBrDalitz = dBrDalitz / brDalitz ;
1216 double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
1217 double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
1218 double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
1219 double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
1220 double dalitzDenomNorm = m_dalitzDenom / nD0 ;
1221 double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
1222
1223 double br2Pip2Pim = m_n2Pip2Pim / nD0 ;
1224 double dBr2Pip2Pim = sqrt( br2Pip2Pim * ( 1. - br2Pip2Pim ) / nD0 ) ;
1225 double fdBr2Pip2Pim= dBr2Pip2Pim / br2Pip2Pim ;
1226 double Pip2Pim2Numer1Norm = m_2Pip2PimNumer1 / nD0 ;
1227 double d2Pip2PimNumer1Norm = Pip2Pim2Numer1Norm * fdBr2Pip2Pim ;
1228 double Pip2Pim2Numer2Norm = m_2Pip2PimNumer2 / nD0 ;
1229 double d2Pip2PimNumer2Norm = Pip2Pim2Numer2Norm * fdBr2Pip2Pim ;
1230 double Pip2Pim2DenomNorm = m_2Pip2PimDenom / nD0 ;
1231 double d2Pip2PimDenomNorm = Pip2Pim2DenomNorm * fdBr2Pip2Pim ;
1232
1233 double brPipPim2Pi0 = m_nPipPim2Pi0 / nD0 ;
1234 double dBrPipPim2Pi0 = sqrt( brPipPim2Pi0 * ( 1. - brPipPim2Pi0 ) / nD0 ) ;
1235 double fdBrPipPim2Pi0= dBrPipPim2Pi0 / brPipPim2Pi0 ;
1236 double Pip2Pim2Numer1Norm = m_PipPim2Pi0Numer1 / nD0 ;
1237 double dPipPim2Pi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrPipPim2Pi0 ;
1238 double Pip2Pim2Numer2Norm = m_PipPim2Pi0Numer2 / nD0 ;
1239 double dPipPim2Pi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrPipPim2Pi0 ;
1240 double Pip2Pim2DenomNorm = m_PipPim2Pi0Denom / nD0 ;
1241 double dPipPim2Pi0DenomNorm = Pip2Pim2DenomNorm * fdBrPipPim2Pi0 ;
1242
1243 double brKSPiPiPi0 = m_nKSPiPiPi0 / nD0 ;
1244 double dBrKSPiPiPi0 = sqrt( brKSPiPiPi0 * ( 1. - brKSPiPiPi0 ) / nD0 ) ;
1245 double fdBrKSPiPiPi0= dBrKSPiPiPi0 / brKSPiPiPi0 ;
1246 double Pip2Pim2Numer1Norm = m_KSPiPiPi0Numer1 / nD0 ;
1247 double dKSPiPiPi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrKSPiPiPi0 ;
1248 double Pip2Pim2Numer2Norm = m_KSPiPiPi0Numer2 / nD0 ;
1249 double dKSPiPiPi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrKSPiPiPi0 ;
1250 double Pip2Pim2DenomNorm = m_KSPiPiPi0Denom / nD0 ;
1251 double dKSPiPiPi0DenomNorm = Pip2Pim2DenomNorm * fdBrKSPiPiPi0 ;
1252
1253 //after the filter was applied
1254 double fil_SL = 0, fil_FlavoredCFD0_0 = 0, fil_FlavoredCFD0_1 = 0, fil_FlavoredCFD0_3 = 0, fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0_0 =0, fil_FlavoredDCSD0_1 =0, fil_FlavoredDCSD0_3 =0, fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0, fil_2Pip2Pim = 0;
1255 for( int i = 0 ; i < 14 ; ++i )
1256 {
1257 fil_SL += m_keptModeCounter[i][9] + m_keptModeCounter[9][i] + m_keptModeCounter[i][8] + m_keptModeCounter[8][i] ;
1258 fil_FlavoredCFD0_0 += m_keptModeCounter[i][1] + m_keptModeCounter[0][i];
1259 fil_FlavoredCFD0_1 += m_keptModeCounter[i][3] + m_keptModeCounter[2][i];
1260 fil_FlavoredCFD0_3 += m_keptModeCounter[i][5] + m_keptModeCounter[4][i];
1261 fil_FlavoredCSD0 += m_keptModeCounter[i][7] + m_keptModeCounter[7][i] + m_keptModeCounter[i][6] + m_keptModeCounter[6][i];
1262 fil_FlavoredDCSD0_0 += m_keptModeCounter[i][0] + m_keptModeCounter[1][i];
1263 fil_FlavoredDCSD0_1 += m_keptModeCounter[i][2] + m_keptModeCounter[3][i];
1264 fil_FlavoredDCSD0_3 += m_keptModeCounter[i][4] + m_keptModeCounter[5][i];
1265 fil_CPPlus += m_keptModeCounter[i][10] + m_keptModeCounter[10][i];
1266 fil_CPMinus += m_keptModeCounter[i][11] + m_keptModeCounter[11][i];
1267 fil_Dalitz += m_keptModeCounter[i][12] + m_keptModeCounter[12][i];
1268 fil_2Pip2Pim += m_keptModeCounter[i][13] + m_keptModeCounter[13][i];
1269 fil_PipPim2Pi0 += m_keptModeCounter[i][14] + m_keptModeCounter[14][i];
1270 fil_KSPiPiPi0 += m_keptModeCounter[i][15] + m_keptModeCounter[15][i];
1271}
1272double nD0_fil = 2. * ( m_nD0D0barEvents - m_nD0D0barDiscarded ) ;
1273double brSL_fil = fil_SL / nD0_fil ;
1274double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
1275double brCF_0_fil = fil_FlavoredCFD0_0 / nD0_fil ;
1276double dBrCF_0_fil = sqrt( brCF_0_fil * ( 1. - brCF_0_fil ) / nD0_fil ) ;
1277double brCF_1_fil = fil_FlavoredCFD0_1 / nD0_fil ;
1278double dBrCF_1_fil = sqrt( brCF_1_fil * ( 1. - brCF_1_fil ) / nD0_fil ) ;
1279double brCF_3_fil = fil_FlavoredCFD0_3 / nD0_fil ;
1280double dBrCF_3_fil = sqrt( brCF_3_fil * ( 1. - brCF_3_fil ) / nD0_fil ) ;
1281double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
1282double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
1283double brDCS_0_fil = fil_FlavoredDCSD0_0 / nD0_fil ;
1284double brDCS_1_fil = fil_FlavoredDCSD0_1 / nD0_fil ;
1285double brDCS_3_fil = fil_FlavoredDCSD0_3 / nD0_fil ;
1286double dBrDCS_0_fil = sqrt( brDCS_0_fil * ( 1. - brDCS_0_fil ) / nD0_fil ) ;
1287double dBrDCS_1_fil = sqrt( brDCS_1_fil * ( 1. - brDCS_1_fil ) / nD0_fil ) ;
1288double dBrDCS_3_fil = sqrt( brDCS_3_fil * ( 1. - brDCS_3_fil ) / nD0_fil ) ;
1289double brCPPlus_fil = fil_CPPlus / nD0_fil ;
1290double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
1291double brCPMinus_fil = fil_CPMinus / nD0_fil ;
1292double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
1293double brDalitz_fil = fil_Dalitz / nD0_fil ;
1294double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
1295double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
1296double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
1297double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
1298double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
1299double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
1300double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
1301double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
1302
1303double br2Pip2Pim_fil = fil_2Pip2Pim / nD0_fil ;
1304double dBr2Pip2Pim_fil = sqrt( br2Pip2Pim_fil * ( 1. - br2Pip2Pim_fil ) / nD0_fil ) ;
1305double fdBr2Pip2Pim_fil = dBr2Pip2Pim_fil / br2Pip2Pim_fil ;
1306double Pip2Pim2Numer1Norm_fil = Pip2Pim2Numer1_fil / nD0_fil ;
1307double d2Pip2PimNumer1Norm_fil = Pip2Pim2Numer1Norm_fil * fdBr2Pip2Pim_fil ;
1308double Pip2Pim2Numer2Norm_fil = Pip2Pim2Numer2_fil / nD0_fil ;
1309double d2Pip2PimNumer2Norm_fil = Pip2Pim2Numer2Norm_fil * fdBr2Pip2Pim_fil ;
1310double Pip2Pim2DenomNorm_fil = Pip2Pim2Denom_fil / nD0_fil ;
1311double d2Pip2PimDenomNorm_fil = Pip2Pim2DenomNorm_fil * fdBr2Pip2Pim_fil ;
1312
1313log << MSG::INFO <<
1314"Original number of SL decays: " << m_nSL <<" ==> Br(SL) = " << brSL << " +- " << dBrSL << endl <<
1315"Filtered number of SL decays: " << fil_SL <<" ==> Br(SL) = " << brSL_fil << " +- " << dBrSL_fil << endl <<
1316"Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<m_nFlavoredCFD0_0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0_0 << endl <<
1317" ==> Br(CF) = " << brCF_0 << " +- " << dBrCF_0 << endl <<
1318"Original number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_1 << "/" << m_nFlavoredDCSD0_1 << endl <<
1319" ==> Br(CF) = " << brCF_1 << " +- " << dBrCF_1 << endl <<
1320"Original number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_3 << "/" << m_nFlavoredDCSD0_3 << endl <<
1321" ==> Br(CF) = " << brCF_3 << " +- " << dBrCF_3 << endl <<
1322" ==> Br(CS) = " << brCS << " +- " << dBrCS << endl <<
1323" ==> Br(DCS) = " << brDCS_0 << " +- " << dBrDCS_0 << endl <<
1324"Filtered number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0_0 << "/" << fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0_0 << endl <<
1325" ==> Br(DCS) = " << brDCS_1 << " +- " << dBrDCS_1 << endl <<
1326"Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<fil_FlavoredCFD0_1 << "/" << fil_FlavoredDCSD0_1 << endl <<
1327" ==> Br(DCS) = " << brDCS_3 << " +- " << dBrDCS_3 << endl <<
1328"Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<fil_FlavoredCFD0_3 << "/" << fil_FlavoredDCSD0_3 << endl <<
1329" ==> Br(CF) = " << brCF_0_fil << " +- " << dBrCF_0_fil << endl <<
1330" ==> Br(CF) = " << brCF_1_fil << " +- " << dBrCF_1_fil << endl <<
1331" ==> Br(CF) = " << brCF_3_fil << " +- " << dBrCF_3_fil << endl <<
1332" ==> Br(CS) = " << brCS_fil << " +- " << dBrCS_fil << endl <<
1333" ==> Br(DCS) = " << brDCS_0_fil << " +- " << dBrDCS_0_fil << endl <<
1334" ==> Br(DCS) = " << brDCS_1_fil << " +- " << dBrDCS_1_fil << endl <<
1335" ==> Br(DCS) = " << brDCS_3_fil << " +- " << dBrDCS_3_fil << endl <<
1336
1337"Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl <<
1338" ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl <<
1339" ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl <<
1340"Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl <<
1341" ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl <<
1342" ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl <<
1343
1344"Original number of Dalitz: " << m_nDalitz << endl <<
1345" ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
1346" y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
1347" numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
1348" denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
1349"Filtered number of Dalitz: " << fil_Dalitz << endl <<
1350" ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
1351" y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
1352" numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
1353" denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg <<
1354
1355"Original number of 2Pip2Pim: " << m_n2Pip2Pim << endl <<
1356" ==> Br(2Pip2Pim) = " << br2Pip2Pim << " +- " << dBr2Pip2Pim <<endl <<
1357" y contrib. numer1 " << Pip2Pim2Numer1Norm << " +- " << d2Pip2PimNumer1Norm << endl <<
1358" numer2 " << Pip2Pim2Numer2Norm << " +- " << d2Pip2PimNumer2Norm << endl <<
1359" denom " << Pip2Pim2DenomNorm << " +- " << d2Pip2PimDenomNorm << endmsg <<
1360"Filtered number of 2Pip2Pim: " << fil_2Pip2Pim << endl <<
1361" ==> Br(2Pip2Pim) = " << br2Pip2Pim_fil << " +- " << dBr2Pip2Pim_fil <<endl <<
1362" y contrib. numer1 " << Pip2Pim2Numer1Norm_fil << " +- " << d2Pip2PimNumer1Norm_fil << endl <<
1363" numer2 " << Pip2Pim2Numer2Norm_fil << " +- " << d2Pip2PimNumer2Norm_fil << endl <<
1364" denom " << Pip2Pim2DenomNorm_fil << " +- " << d2Pip2PimDenomNorm_fil << endmsg ;
1365
1366
1367HepMatrix differencetoWeight(17,17,0);
1368for( int i = 0 ; i < 16 ; ++i ) {
1369 for( int j = 0 ; j < 16 ; ++j ) {
1370 if (m_modeCounter[i][j] != 0) {
1371 differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] - m_weights[i][j]; } } }
1372
1373 log << MSG::INFO << "Weight matrix" << m_weights << endmsg ;
1374 log << MSG::INFO << "D0 Modes before filter" << m_modeCounter << endmsg ;
1375 log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter << endmsg ;
1376 log << MSG::INFO << "Compare percentage difference to weights " << differencetoWeight << endmsg ;
1377
1378 //Calculating and comparing Y before/after the filter
1379 // To calculate y, we want half the CS BF
1380 brCS /= 2. ;
1381 dBrCS /= 2. ;
1382 brCS_fil /= 2 ;
1383 dBrCS_fil /= 2 ;
1384
1385 double y, dy, y_fil, dy_fil ;
1386 double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
1387
1388 //Original MC
1389 double rCF_0 = m_rCF_0 ;
1390 double zCF_0 = m_zCF_0 ;
1391 double rCF2_0 = rCF_0 * rCF_0 ;
1392 double zCF2_0 = zCF_0 * zCF_0 ;
1393 double rCF_1 = m_rCF_1 ;
1394 double zCF_1 = m_zCF_1 ;
1395 double rCF2_1 = rCF_1 * rCF_1 ;
1396 double zCF2_1 = zCF_1 * zCF_1 ;
1397 double rCF_3 = m_rCF_3 ;
1398 double zCF_3 = m_zCF_3 ;
1399 double rCF2_3 = rCF_3 * rCF_3 ;
1400 double zCF2_3 = zCF_3 * zCF_3 ;
1401 double rCS = m_rCS ;
1402 double zCS = m_zCS ;
1403 double rCS2 = rCS * rCS ;
1404 double zCS2 = zCS * zCS ;
1405
1406 double wCF_0 = ( m_wCFSign_0 ? 1. : -1. ) * sqrt( 4. - zCF2_0 ) ;
1407 double wCF_1 = ( m_wCFSign_1 ? 1. : -1. ) * sqrt( 4. - zCF2_1 ) ;
1408 double wCF_3 = ( m_wCFSign_3 ? 1. : -1. ) * sqrt( 4. - zCF2_3 ) ;
1409 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
1410
1411 double numFactCF_0 =
1412 -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x ) ;
1413 double numFactCF_1 =
1414 -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x ) ;
1415 double numFactCF_3 =
1416 -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x ) ;
1417 double numFactCS =
1418 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
1419 double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0 ;
1420 double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1 ;
1421 double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3 ;
1422 double denFactCS = 0.5 * rCS2 * zCS2 ;
1423 double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ;
1424 double Pip2Pim2NumerNorm = Pip2Pim2Numer1Norm + Pip2Pim2Numer2Norm ;
1425
1426 double num =
1427 brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 + brCF_3 * numFactCF_3 + brCS * numFactCS +
1428 dalitzNumerNorm + Pip2Pim2NumerNorm;
1429 double den =
1430 1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 - brCF_3 * denFactCF_3 - brCS * denFactCS +
1431 dalitzDenomNorm + Pip2Pim2DenomNorm;
1432 y = num / den ;
1433 dy = sqrt(
1434 ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
1435 ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
1436 ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 *
1437 ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 +
1438 ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 *
1439 ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 +
1440 ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 *
1441 ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 +
1442 ( numFactCS * den + denFactCS * num ) * dBrCS *
1443 ( numFactCS * den + denFactCS * num ) * dBrCS +
1444 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz *
1445 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz +
1446 ( Pip2Pim2NumerNorm * den + Pip2Pim2DenomNorm * num ) * fdBrDalitz *
1447 ( Pip2Pim2NumerNorm * den + Pip2Pim2DenomNorm * num ) * fdBrDalitz
1448 ) / ( den * den ) ;
1449
1450 double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] + m_modeCounter[4][6] + m_modeCounter[5][6] ;
1451 double n_cpminussl = m_modeCounter[7][4] + m_modeCounter[7][5] + m_modeCounter[4][7] + m_modeCounter[5][7] ;
1452 double a_cpsl = ( n_cpplussl * m_nCPMinus ) / ( n_cpminussl * m_nCPPlus ) ;
1453 double b_cpsl = ( n_cpminussl * m_nCPPlus ) / ( n_cpplussl * m_nCPMinus ) ;
1454 y_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
1455 dy_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) ) ;
1456
1457
1458 //Filtered MC
1459 double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
1460 double Pip2Pim2NumerNorm_fil = Pip2Pim2Numer1Norm_fil + Pip2Pim2Numer2Norm_fil ;
1461 double num_fil =
1462 brCPPlus_fil - brCPMinus_fil + brCF_0_fil * numFactCF_0 + brCF_1_fil * numFactCF_1 + brCF_3_fil * numFactCF_3 + brCS_fil * numFactCS +
1463 dalitzNumerNorm_fil + Pip2Pim2NumerNorm_fil;
1464 double den_fil =
1465 1. - brCPPlus_fil - brCPMinus_fil - brCF_0_fil * denFactCF_0 - brCF_1_fil * denFactCF_1 - brCF_3_fil * denFactCF_3 - brCS_fil * denFactCS +
1466 dalitzDenomNorm_fil + Pip2Pim2DenomNorm_fil;
1467 y_fil = num_fil / den_fil ;
1468 dy_fil = sqrt(
1469 ( num_fil + den_fil ) * ( num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil +
1470 ( num_fil - den_fil ) * ( num_fil - den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil +
1471 ( numFactCF_0 * den_fil + denFactCF_0 * num_fil ) * dBrCF_0_fil *
1472 ( numFactCF_0 * den_fil + denFactCF_0 * num_fil ) * dBrCF_0_fil +
1473 ( numFactCF_1 * den_fil + denFactCF_1 * num_fil ) * dBrCF_1_fil *
1474 ( numFactCF_1 * den_fil + denFactCF_1 * num_fil ) * dBrCF_1_fil +
1475 ( numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil *
1476 ( numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil +
1477 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil *
1478 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
1479 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil *
1480 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil +
1481 ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) * fdBr2Pip2Pim_fil *
1482 ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) * fdBr2Pip2Pim_fil
1483 ) / ( den_fil * den_fil ) ;
1484
1485 double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] + m_keptModeCounter[4][6] + m_keptModeCounter[5][6] ;
1486 double fil_cpminussl = m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] + m_keptModeCounter[5][7] ;
1487 a_cpsl = ( fil_cpplussl * fil_CPMinus ) / ( fil_cpminussl * fil_CPPlus ) ;
1488 b_cpsl = ( fil_cpminussl * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus ) ;
1489 y_fil_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
1490 dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ;
1491
1492
1493
1494 log << MSG::INFO <<"BR Method : Parent MC ==> y = " << y << " +- " << dy << endmsg ;
1495 log << MSG::INFO <<"BR Method : Filtered MC ==> y = " << y_fil << " +- " << dy_fil << endmsg ;
1496 log << MSG::INFO <<"CP-SL doubletag : Parent MC ==> y = " << y_cpsl << " +- " << dy_cpsl <<endl<<"Previous line should be equivalent with 0 as parent MC not correalated"<< endmsg ;
1497 log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " << dy_fil_cpsl << endmsg ;
1498
1499 }*/
1500
1501return StatusCode::SUCCESS;
1502}
◆ findD0Decay()
| vector< double > QCMCFilter::findD0Decay | ( | int | charm | ) |
Definition at line 1507 of file QCMCFilter.cxx.
1508{
1509 MsgStream log(msgSvc(), name());
1510 log << MSG::INFO <<" In findD0Decay " << endmsg ;
1511 vector<double> d0_info(9);
1512 for(int i=0; i< 8; i++) d0_info[i]=0;
1513 double decay_mode = 18; //d0_info[0] = decay mode;
1514 double r2D0 = -1; //d0_info[1] = r2D0;
1515 double deltaD0 = -1; //d0_info[2] = deltaD0;
1516 double RD0 = 1; //d0_info[3] = RD0;
1517 double FCP = 1; //d0_info[4] = FCP;
1518
1521 // Number of each partile type
1522 int kPiPlus = 0;// num[0] = Pi+, a0+
1523 int kPiMinus = 1;// num[1] = Pi-, a0-
1524 int kPi0 = 2;// num[2] = Pi0
1525 int kEta = 3;// num[3] = Eta
1526 int kEtaPrime = 4;// num[4] = Eta Prime
1527 int kNeutralScalar = 5;// num[5] = Neutral Scalar : f0, f'0, f0(1500), a0, (f_2 can be included here because it is only used in f_2 pi0 decay)
1528 int kUFVPlus = 6;// num[6] = unflavored positive (axial) vectors: rho+, a1+, rho(2S)+
1529 int kUFVMinus = 7;// num[7] = unflavored negative (axial) vectors: rho-, a1-, rho(2S)-
1530 int kRho0 = 8;// num[8] = Rho0, Rho(2S)0
1531 int kOmega = 9;// num[9] = Omega
1532 int kPhi = 10;// num[10] = Phi
1533 int kKPlus = 11;// num[11] = K+
1534 int kKMinus = 12;// num[12] = K-
1535 int kK0S = 13;// num[13] = K short
1536 int kK0L = 14;// num[14] = K long
1537 int kFVPlus = 15;// num[15] = flavored positive (axial) vectors: K*+, K_1+, K_2*+
1538 int kFVMinus = 16;// num[16] = flavored negative (axial) vectors: K*-, K_1-, K_2*-
1539 int kKStar0 = 17;// num[17] = K*0
1540 int kKStar0Bar = 18;// num[18] = anti-K*0
1541 int kK10 = 19;// num[19] = K1_0
1542 int kK10Bar = 20;// num[20] = anti-K1_0
1543 int kLeptonPlus = 21;// num[21] = LeptonPlus: e+, mu+
1544 int kLeptonMinus = 22;// num[22] = LeptonMinus: e-, mu-
1545 int kNeutrino = 23;// num[23] = Neutrino: nu, nubar
1546 int kGamma = 24;// num[24] = gamma
1547 int kUnknown = 25;// num[25] = Not identified particle
1548 //int kGammaFSR GammaFSR is ignored as we are interested in the PreFSR quantum #
1549
1550 // Get McParticle collection
1551 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
1552 int d0_pdg = charm == 1 ? kD0ID : kD0BarID;
1553
1554 HepLorentzVector piPlus, piMinus, k0;
1555 vector<HepLorentzVector> piPlus_4pi, piMinus_4pi, pi0_p4;
1556 vector<HepLorentzVector> kPlus_kkpipi, kMinus_kkpipi;
1557 piPlus_4pi.clear();piMinus_4pi.clear(), pi0_p4.clear();
1558 kPlus_kkpipi.clear();kMinus_kkpipi.clear();
1559
1560 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
1561 {
1562 int pdg_code = (*it)->particleProperty();
1563 if ((*it)->primaryParticle()) continue;
1565 int mother_pdg = 0;
1566
1567 //finds if the particle is related to the d0 we need
1568 mother_pdg = it2.particleProperty();
1569
1570 //special cases
1571 if (pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID || pdg_code == kK0BarID || pdg_code == kKDPStar0ID || pdg_code == kKDPStar0BarID
1572 || pdg_code == kK10ID || pdg_code == kK10BarID
1573 || pdg_code == kK0Star0ID || pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID ) continue; //don't count special intermediate states
1574
1575 if (mother_pdg == kK0ID || mother_pdg == kK0BarID) {//Looks another step up if mother was k0 or k0bar
1576 it2 = it2.mother();
1577 mother_pdg = it2.particleProperty(); }
1578
1579 if (mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID) {//Looks another step up if mother was k*0 or k*0bar
1580 //if code is found to be K*0 -> K+Pi- we save it as a KStar0 and K*0Bar -> K-Pi+ as KStar0Bar
1581 //if is it a neutral decay K*0(Bar)->K0(bar) Pi0 or Gamma. We save them as K0(bar) and Pi0 or Gamma.
1582 if (pdg_code == kPiPlusID || pdg_code == kPiMinusID) continue;
1583 if (mother_pdg == kKStar0ID && pdg_code == kKPlusID) pdg_code = kKStar0ID;
1584 if (mother_pdg == kKStar0BarID && pdg_code == kKMinusID) pdg_code = kKStar0BarID;
1585 it2 = it2.mother();
1586 mother_pdg = it2.particleProperty();
1587 }
1588
1589 if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID ) { // Looks another step up if mother was K_0
1590 it2 = it2.mother();
1591 mother_pdg = it2.particleProperty();
1592 }
1593
1594 if (mother_pdg == kK10ID || mother_pdg == kK10BarID) {//Looks another step up if mother was k_10 or k_10bar
1595 it2 = it2.mother();
1596 mother_pdg = it2.particleProperty();
1597 }
1598
1599
1600
1601 //Identifies what particles are the daughters
1602 if (mother_pdg == d0_pdg) //Found daughter
1603 {
1609 else if (pdg_code == kF0ID || pdg_code == kA00ID
1610 || pdg_code == kFPrime0ID|| pdg_code == kF0m1500ID
1611 || pdg_code == kF2ID) num[5]++;
1612 else if (pdg_code == kRhoPlusID || pdg_code == kRho2SPlusID
1613 || pdg_code == kA1PlusID ) num[6]++;
1614 else if (pdg_code == kRhoMinusID || pdg_code == kRho2SMinusID
1615 || pdg_code == kA1MinusID) num[7]++;
1623 else if (pdg_code == kKStarPlusID || pdg_code == kK1PlusID
1624 || pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID
1625 || pdg_code == kK0StarPlusID) num[15]++;
1626 else if (pdg_code == kKStarMinusID || pdg_code == kK1MinusID
1627 || pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID
1628 || pdg_code == kK0StarMinusID ) num[16]++;
1635 else if (pdg_code == kNuEID || pdg_code == kNuEBarID
1636 || pdg_code == kNuMuID || pdg_code == kNuMuBarID ) num[23]++;
1638 else if (pdg_code == kGammaFSRID ) continue;
1639 else {
1640 num[25]++;
1641 cout <<"Unknown particle: "<< pdg_code << endl;
1642 }
1643
1644 //if (pdg_code == kA10ID ) num[0]++; // not present
1645
1646 //Enter Momentums for Dalitz canidates
1647 //It would be prefered if we could get PreFSR Momentum
1648 //m_D0daughter_P4.push_back( (*it)->initialFourMomentum().vect() );
1649 //m_D0daughter_PID.push_back( pdg_code );
1650
1655 HepLorentzVector daughterP4;
1656 if(piPlus_4pi.size()>2||piMinus_4pi.size()>2)continue;
1657 if(kPlus_kkpipi.size()>1||kMinus_kkpipi.size()>1)continue;
1658 //log << MSG::INFO << "in D0->pipipipi;" << endmsg;
1659 if (pdg_code == kPiPlusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpion);piPlus_4pi.push_back(daughterP4);}
1660 if (pdg_code == kPiMinusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpion);piMinus_4pi.push_back(daughterP4);}
1661 if (pdg_code == kKPlusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmkaon);kPlus_kkpipi.push_back(daughterP4);}
1662 if (pdg_code == kKMinusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmkaon);kMinus_kkpipi.push_back(daughterP4);}
1663 if (pdg_code == kPi0ID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpi0);pi0_p4.push_back(daughterP4);}
1664 }
1665 }
1666
1667 // D decay daughters have been tabulated. Now, classify decay.
1668 int nDaughters = 0 ;
1669 for( int i = 0 ; i < 26 ; i++ )
1670 {
1671 nDaughters += num[ i ] ;
1672 }
1673
1680 int nCPPlusEig = num[ kNeutralScalar ] + num[ kRho0 ] + num[ kOmega ] + num[ kPhi ] + num[ kK0S ] + num[ kGamma ] ;
1682 int nCPEig = nCPPlusEig + nCPMinusEig ;
1685
1686
1687 //check Dalitz modes: KsPi+Pi- or KlPi+Pi-
1688 double mnm_gen = 0. ;
1689 double mpn_gen = 0. ;
1690 double mpm_gen = 0. ;
1693
1694 //check K3Pi modes:
1695
1696
1698 {
1699 decay_mode = kDalitz ;
1700 //k0.boost(-0.011, 0, 0);
1701 //piMinus.boost(-0.011, 0, 0);
1702 //piPlus.boost(-0.011, 0, 0);
1703 //mnm_gen = (k0 + piMinus).m2();
1704 //mpn_gen = (k0 + piPlus).m2();
1705 //mpm_gen = (piPlus + piMinus).m2();
1707 }
1709 {
1710 decay_mode = k2Pip2Pim;
1711 }
1713 {
1714 decay_mode = kPipPim2Pi0;
1715 }
1717 {
1718 decay_mode = kK0PiPiPi0;
1720 }
1721 if( num[ kPiPlus ] == 1&&num[ kPiMinus ]==1 &&num[ kKPlus ] == 1&&num[ kKMinus ]==1 && nDaughters == 4 )
1722 {
1723 decay_mode = kKKPiPi;
1724 }
1725
1726 // The order of if-statements below matters.
1727 if( decay_mode == kDalitz )
1728 {
1729 ++m_nDalitz ;
1730
1731 // Dalitz t(8) ;
1732 // TComplex D0, D0bar;
1733 // if (m_dalitzModel == 0) { //Cleo model
1734 // D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1735 // D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1736 // if (m_dalitzModel == 1) { //Babar model
1737 // D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1738 // D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1739 // if (m_dalitzModel == 2) { //Belle model
1740 // D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1741 // D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1742
1743 complex<double> D0, D0bar;
1744 vector<double> k0l;vector<double> pip;vector<double> pim;
1745 k0l.clear();pip.clear();pim.clear();
1746 k0l.push_back(k0[0]);k0l.push_back(k0[1]);k0l.push_back(k0[2]);k0l.push_back(k0[3]);
1747 pip.push_back(piPlus[0]);pip.push_back(piPlus[1]);pip.push_back(piPlus[2]);pip.push_back(piPlus[3]);
1748 pim.push_back(piMinus[0]);pim.push_back(piMinus[1]);pim.push_back(piMinus[2]);pim.push_back(piMinus[3]);
1749
1750 if(isKsPiPi){
1752 D0 = tKSpipi.Amp_PFT(k0l,pip,pim);
1753
1754 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
1755 cpk0l.clear();pip.clear();pim.clear();
1756 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
1757 cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
1758 cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
1759 D0bar = tKSpipi.Amp_PFT(cpk0l, cppip, cppim);
1760 }else{
1761 D0ToKLpipi2018 tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx is not the same as used in other analyses
1762 D0 = tKLpipi.Amp_PFT(k0l,pip,pim);
1763 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
1764 cpk0l.clear();pip.clear();pim.clear();
1765 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
1766 cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
1767 cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
1768 D0bar = tKLpipi.Amp_PFT(cpk0l, cppip, cppim);
1769 }
1770
1771 if( charm == 1){
1772
1773 r2D0 = std::norm(D0bar)/std::norm(D0);
1774 double temp = std::arg(D0)-std::arg(D0bar);
1775 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
1776 while (temp < -TMath::Pi())
1777 {
1778 temp += 2.0 * TMath::Pi();
1779 }
1780 while (temp > TMath::Pi())
1781 {
1782 temp -= 2.0 * TMath::Pi();
1783 }
1784 deltaD0 = temp;
1785 deltaD0 = isKsPiPi ? deltaD0 : (deltaD0+PI);
1786
1789 if( mpn_gen - mnm_gen > 0. )
1790 {
1791 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1792 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1793 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1794 }
1795 }
1796 else {
1797 r2D0 = std::norm(D0)/std::norm(D0bar);
1798 double temp = std::arg(D0bar)-std::arg(D0);
1799 while (temp < -TMath::Pi())
1800 {
1801 temp += 2.0 * TMath::Pi();
1802 }
1803 while (temp > TMath::Pi())
1804 {
1805 temp -= 2.0 * TMath::Pi();
1806 }
1807 deltaD0 = temp;
1808 deltaD0 = isKsPiPi ? deltaD0 : (deltaD0+PI);
1809
1812 if( mpn_gen - mnm_gen < 0. )
1813 {
1814 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1815 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1816 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1817 }
1818 }
1819 }else if( decay_mode == k2Pip2Pim ){
1820
1821 ++m_n2Pip2Pim ;
1822
1823 HepLorentzVector pip1_4pi = (piPlus_4pi.at(0));
1824 HepLorentzVector pim1_4pi = (piMinus_4pi.at(0));
1825 HepLorentzVector pip2_4pi = (piPlus_4pi.at(1));
1826 HepLorentzVector pim2_4pi = (piMinus_4pi.at(1));
1827 pip1_4pi.boost(-0.011, 0, 0);
1828 pip2_4pi.boost(-0.011, 0, 0);
1829 pim1_4pi.boost(-0.011, 0, 0);
1830 pim2_4pi.boost(-0.011, 0, 0);
1831 std::complex<double> D0, D0bar;
1832 //EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
1833 //EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
1834 //EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
1835 //EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
1836 std::vector<double> pip1; pip1.clear(); pip1.push_back(pip1_4pi[0]);pip1.push_back(pip1_4pi[1]);pip1.push_back(pip1_4pi[2]);pip1.push_back(pip1_4pi[3]);
1837 std::vector<double> pim1; pim1.clear(); pim1.push_back(pim1_4pi[0]);pim1.push_back(pim1_4pi[1]);pim1.push_back(pim1_4pi[2]);pim1.push_back(pim1_4pi[3]);
1838 std::vector<double> pip2; pip2.clear(); pip2.push_back(pip2_4pi[0]);pip2.push_back(pip2_4pi[1]);pip2.push_back(pip2_4pi[2]);pip2.push_back(pip2_4pi[3]);
1839 std::vector<double> pim2; pim2.clear(); pim2.push_back(pim2_4pi[0]);pim2.push_back(pim2_4pi[1]);pim2.push_back(pim2_4pi[2]);pim2.push_back(pim2_4pi[3]);
1840
1842 D0 = t2pip2pim.Amp(pip1,pim1,pip2,pim2);
1843 //EvtVector4R cppim1(pip1_4pi[3], -1*pip1_4pi[0], -1*pip1_4pi[1], -1*pip1_4pi[2]);
1844 //EvtVector4R cppip1(pim1_4pi[3], -1*pim1_4pi[0], -1*pim1_4pi[1], -1*pim1_4pi[2]);
1845 //EvtVector4R cppim2(pip2_4pi[3], -1*pip2_4pi[0], -1*pip2_4pi[1], -1*pip2_4pi[2]);
1846 //EvtVector4R cppip2(pim2_4pi[3], -1*pim2_4pi[0], -1*pim2_4pi[1], -1*pim2_4pi[2]);
1847 std::vector<double> cppip1; cppip1.clear(); cppip1.push_back(-pim1_4pi[0]);cppip1.push_back(-pim1_4pi[1]);cppip1.push_back(-pim1_4pi[2]);cppip1.push_back(pim1_4pi[3]);
1848 std::vector<double> cppim1; cppim1.clear(); cppim1.push_back(-pip1_4pi[0]);cppim1.push_back(-pip1_4pi[1]);cppim1.push_back(-pip1_4pi[2]);cppim1.push_back(pip1_4pi[3]);
1849 std::vector<double> cppip2; cppip2.clear(); cppip2.push_back(-pim2_4pi[0]);cppip2.push_back(-pim2_4pi[1]);cppip2.push_back(-pim2_4pi[2]);cppip2.push_back(pim2_4pi[3]);
1850 std::vector<double> cppim2; cppim2.clear(); cppim2.push_back(-pip2_4pi[0]);cppim2.push_back(-pip2_4pi[1]);cppim2.push_back(-pip2_4pi[2]);cppim2.push_back(pip2_4pi[3]);
1851
1852 D0bar = t2pip2pim.Amp(cppip1, cppim1, cppip2, cppim2);
1853
1854 if( charm == 1){
1855
1856 r2D0 = std::norm(D0bar)/std::norm(D0);
1857 double temp = std::arg(D0)-std::arg(D0bar);
1858 while (temp < -TMath::Pi())
1859 {
1860 temp += 2.0 * TMath::Pi();
1861 }
1862 while (temp > TMath::Pi())
1863 {
1864 temp -= 2.0 * TMath::Pi();
1865 }
1866 deltaD0 = temp;
1867
1870 //if( mpn_gen - mnm_gen > 0. )
1871 //{
1872 m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1873 m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1874 m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd ;
1875 //}
1876 }
1877 else {
1878 r2D0 = std::norm(D0)/std::norm(D0bar);
1879 double temp = std::arg(D0bar)-std::arg(D0);
1880 while (temp < -TMath::Pi())
1881 {
1882 temp += 2.0 * TMath::Pi();
1883 }
1884 while (temp > TMath::Pi())
1885 {
1886 temp -= 2.0 * TMath::Pi();
1887 }
1888 deltaD0 = temp;
1889
1892 //if( mpn_gen - mnm_gen < 0. )
1893 //{
1894 m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1895 m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1896 m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd ;
1897 //}
1898 }
1899 }else if( decay_mode == kPipPim2Pi0 ){
1900
1901 ++m_nPipPim2Pi0 ;
1902
1903
1904 HepLorentzVector pip1_4pi = (piPlus_4pi.at(0));
1905 HepLorentzVector pim1_4pi = (piMinus_4pi.at(0));
1906 HepLorentzVector pi01_p4 = (pi0_p4.at(0));
1907 HepLorentzVector pi02_p4 = (pi0_p4.at(1));
1908 pip1_4pi.boost(-0.011, 0, 0);
1909 pim1_4pi.boost(-0.011, 0, 0);
1910 pi01_p4.boost(-0.011, 0, 0);
1911 pi02_p4.boost(-0.011, 0, 0);
1912 std::complex<double> D0, D0bar;
1913 //EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
1914 //EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
1915 //EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
1916 //EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
1917 std::vector<double> pip1; pip1.clear(); pip1.push_back(pip1_4pi[0]);pip1.push_back(pip1_4pi[1]);pip1.push_back(pip1_4pi[2]);pip1.push_back(pip1_4pi[3]);
1918 std::vector<double> pim1; pim1.clear(); pim1.push_back(pim1_4pi[0]);pim1.push_back(pim1_4pi[1]);pim1.push_back(pim1_4pi[2]);pim1.push_back(pim1_4pi[3]);
1919 std::vector<double> pi01; pi01.clear(); pi01.push_back(pi01_p4[0]);pi01.push_back(pi01_p4[1]);pi01.push_back(pi01_p4[2]);pi01.push_back(pi01_p4[3]);
1920 std::vector<double> pi02; pi02.clear(); pi02.push_back(pi02_p4[0]);pi02.push_back(pi02_p4[1]);pi02.push_back(pi02_p4[2]);pi02.push_back(pi02_p4[3]);
1921
1923 D0 = tpippim2pi0.Amp(pip1,pim1,pi01,pi02);
1924 //EvtVector4R cppim1(pip1_4pi[3], -1*pip1_4pi[0], -1*pip1_4pi[1], -1*pip1_4pi[2]);
1925 //EvtVector4R cppip1(pim1_4pi[3], -1*pim1_4pi[0], -1*pim1_4pi[1], -1*pim1_4pi[2]);
1926 //EvtVector4R cppim2(pip2_4pi[3], -1*pip2_4pi[0], -1*pip2_4pi[1], -1*pip2_4pi[2]);
1927 //EvtVector4R cppip2(pim2_4pi[3], -1*pim2_4pi[0], -1*pim2_4pi[1], -1*pim2_4pi[2]);
1928 std::vector<double> cppip1; cppip1.clear(); cppip1.push_back(-pim1_4pi[0]);cppip1.push_back(-pim1_4pi[1]);cppip1.push_back(-pim1_4pi[2]);cppip1.push_back(pim1_4pi[3]);
1929 std::vector<double> cppim1; cppim1.clear(); cppim1.push_back(-pip1_4pi[0]);cppim1.push_back(-pip1_4pi[1]);cppim1.push_back(-pip1_4pi[2]);cppim1.push_back(pip1_4pi[3]);
1930 std::vector<double> cppi01; cppi01.clear(); cppi01.push_back(-pi01_p4[0]);cppi01.push_back(-pi01_p4[1]);cppi01.push_back(-pi01_p4[2]);cppi01.push_back(pi01_p4[3]);
1931 std::vector<double> cppi02; cppi02.clear(); cppi02.push_back(-pi02_p4[0]);cppi02.push_back(-pi02_p4[1]);cppi02.push_back(-pi02_p4[2]);cppi02.push_back(pi02_p4[3]);
1932
1933 D0bar = tpippim2pi0.Amp(cppip1, cppim1, cppi01, cppi02);
1934
1935 if( charm == 1){
1936
1937 r2D0 = std::norm(D0bar)/std::norm(D0);
1938 double temp = std::arg(D0)-std::arg(D0bar);
1939 while (temp < -TMath::Pi())
1940 {
1941 temp += 2.0 * TMath::Pi();
1942 }
1943 while (temp > TMath::Pi())
1944 {
1945 temp -= 2.0 * TMath::Pi();
1946 }
1947 deltaD0 = temp;
1948
1951 //if( mpn_gen - mnm_gen > 0. )
1952 //{
1953 m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
1954 m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
1955 m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd ;
1956 //}
1957 }
1958 else {
1959 r2D0 = std::norm(D0)/std::norm(D0bar);
1960 double temp = std::arg(D0bar)-std::arg(D0);
1961 while (temp < -TMath::Pi())
1962 {
1963 temp += 2.0 * TMath::Pi();
1964 }
1965 while (temp > TMath::Pi())
1966 {
1967 temp -= 2.0 * TMath::Pi();
1968 }
1969 deltaD0 = temp;
1970
1973
1974 m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
1975 m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
1976 m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd ;
1977 }
1978 }else if( decay_mode == kK0PiPiPi0 )
1979 {
1980 ++m_nK0PiPiPi0 ;
1981
1982 // Dalitz t(8) ;
1983 // TComplex D0, D0bar;
1984 // if (m_dalitzModel == 0) { //Cleo model
1985 // D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1986 // D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1987 // if (m_dalitzModel == 1) { //Babar model
1988 // D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1989 // D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1990 // if (m_dalitzModel == 2) { //Belle model
1991 // D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1992 // D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1993 HepLorentzVector pip_kspipipi0 = (piPlus_4pi.at(0));
1994 HepLorentzVector pim_kspipipi0 = (piMinus_4pi.at(0));
1995 HepLorentzVector k0_kspipipi0 = k0;
1996 HepLorentzVector pi0_kspipipi0 = (pi0_p4.at(0));
1997 pip_kspipipi0.boost(-0.011, 0, 0);
1998 pim_kspipipi0.boost(-0.011, 0, 0);
1999 k0_kspipipi0.boost(-0.011, 0, 0);
2000 pi0_kspipipi0.boost(-0.011, 0, 0);
2001
2002 complex<double> D0, D0bar;
2003 vector<double> k0l;vector<double> pip;vector<double> pim; vector<double> pi0;
2004 k0l.clear();pip.clear();pim.clear();pi0.clear();
2005 k0l.push_back(k0_kspipipi0[0]);k0l.push_back(k0_kspipipi0[1]);k0l.push_back(k0_kspipipi0[2]);k0l.push_back(k0_kspipipi0[3]);
2006 pip.push_back(pip_kspipipi0[0]);pip.push_back(pip_kspipipi0[1]);pip.push_back(pip_kspipipi0[2]);pip.push_back(pip_kspipipi0[3]);
2007 pim.push_back(pim_kspipipi0[0]);pim.push_back(pim_kspipipi0[1]);pim.push_back(pim_kspipipi0[2]);pim.push_back(pim_kspipipi0[3]);
2008 pi0.push_back(pi0_kspipipi0[0]);pi0.push_back(pi0_kspipipi0[1]);pi0.push_back(pi0_kspipipi0[2]);pi0.push_back(pi0_kspipipi0[3]);
2009
2010 if(isKsPiPiPi0){
2012 D0 = tKSpipipi0.Amp(k0l,pip,pim,pi0);
2013
2014 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;vector<double> cppi0;
2015 cpk0l.clear();cppip.clear();cppim.clear();cppi0.clear();
2016 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2017 cppim.push_back(-pip[0]);cppim.push_back(-pip[1]);cppim.push_back(-pip[2]);cppim.push_back(pip[3]);
2018 cppip.push_back(-pim[0]);cppip.push_back(-pim[1]);cppip.push_back(-pim[2]);cppip.push_back(pim[3]);
2019 cppi0.push_back(-pi0[0]);cppi0.push_back(-pi0[1]);cppi0.push_back(-pi0[2]);cppi0.push_back(pi0[3]);
2020 D0bar = tKSpipipi0.Amp(cpk0l, cppip, cppim, cppi0);
2021 }else{
2023 D0 = tKSpipipi0.Amp(k0l,pip,pim,pi0);
2024
2025 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;vector<double> cppi0;
2026 cpk0l.clear();cppip.clear();cppim.clear();cppi0.clear();
2027 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2028 cppim.push_back(-pip[0]);cppim.push_back(-pip[1]);cppim.push_back(-pip[2]);cppim.push_back(pip[3]);
2029 cppip.push_back(-pim[0]);cppip.push_back(-pim[1]);cppip.push_back(-pim[2]);cppip.push_back(pim[3]);
2030 cppi0.push_back(-pi0[0]);cppi0.push_back(-pi0[1]);cppi0.push_back(-pi0[2]);cppi0.push_back(pi0[3]);
2031 D0bar = tKSpipipi0.Amp(cpk0l, cppip, cppim, cppi0);
2032 //D0ToKLpipi tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx is not the same as used in other analyses
2033 //D0 = tKLpipi.Amp_PFT(k0l,pip,pim);
2034 //vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
2035 //cpk0l.clear();pip.clear();pim.clear();
2036 //cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2037 //cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
2038 //cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
2039 //D0bar = tKLpipi.Amp_PFT(cpk0l, cppip, cppim);
2040 }
2041
2042 if( charm == 1){
2043
2044 r2D0 = std::norm(D0bar)/std::norm(D0);
2045 double temp = std::arg(D0)-std::arg(D0bar);
2046 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
2047 while (temp < -TMath::Pi())
2048 {
2049 temp += 2.0 * TMath::Pi();
2050 }
2051 while (temp > TMath::Pi())
2052 {
2053 temp -= 2.0 * TMath::Pi();
2054 }
2055 deltaD0 = temp;
2056 deltaD0 = isKsPiPiPi0 ? (deltaD0+PI) : (deltaD0);
2057
2060 if( mpn_gen - mnm_gen > 0. )
2061 {
2062 m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
2063 m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
2064 m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd ;
2065 }
2066 }
2067 else {
2068 r2D0 = std::norm(D0)/std::norm(D0bar);
2069 double temp = std::arg(D0bar)-std::arg(D0);
2070 while (temp < -TMath::Pi())
2071 {
2072 temp += 2.0 * TMath::Pi();
2073 }
2074 while (temp > TMath::Pi())
2075 {
2076 temp -= 2.0 * TMath::Pi();
2077 }
2078 deltaD0 = temp;
2079 deltaD0 = isKsPiPiPi0 ? (deltaD0+PI) : (deltaD0);
2080
2083 if( mpn_gen - mnm_gen < 0. )
2084 {
2085 m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
2086 m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
2087 m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd ;
2088 }
2089 }
2090 }else if( decay_mode == kKKPiPi ){
2091
2092 ++m_nKKPiPi ;
2093
2094 complex<double> D0, D0bar;
2095 vector<double> kp; vector<double> km; vector<double> pip;vector<double> pim;
2096 double pdau[16];
2097 pip.clear();pim.clear();
2098 kp.clear();km.clear();
2099
2100 HepLorentzVector pip_kkpipi = (piPlus_4pi.at(0));
2101 HepLorentzVector pim_kkpipi = (piMinus_4pi.at(0));
2102 HepLorentzVector kp_kkpipi = (kPlus_kkpipi.at(0));
2103 HepLorentzVector km_kkpipi = (kMinus_kkpipi.at(0));
2104
2105 pip.push_back(pip_kkpipi[0]);pip.push_back(pip_kkpipi[1]);pip.push_back(pip_kkpipi[2]);pip.push_back(pip_kkpipi[3]);
2106 pim.push_back(pim_kkpipi[0]);pim.push_back(pim_kkpipi[1]);pim.push_back(pim_kkpipi[2]);pim.push_back(pim_kkpipi[3]);
2107 kp.push_back(kp_kkpipi[0]);kp.push_back(kp_kkpipi[1]);kp.push_back(kp_kkpipi[2]);kp.push_back(kp_kkpipi[3]);
2108 km.push_back(km_kkpipi[0]);km.push_back(km_kkpipi[1]);km.push_back(km_kkpipi[2]);km.push_back(km_kkpipi[3]);
2109
2110 pdau[0]=kp_kkpipi[0]; pdau[1]=kp_kkpipi[1]; pdau[2]=kp_kkpipi[2]; pdau[3]=kp_kkpipi[3];
2111 pdau[4]=km_kkpipi[0]; pdau[5]=km_kkpipi[1]; pdau[6]=km_kkpipi[2]; pdau[7]=km_kkpipi[3];
2112 pdau[8]=pip_kkpipi[0];pdau[9]=pip_kkpipi[1];pdau[10]=pip_kkpipi[2];pdau[11]=pip_kkpipi[3];
2113 pdau[12]=pim_kkpipi[0];pdau[13]=pim_kkpipi[1];pdau[14]=pim_kkpipi[2];pdau[15]=pim_kkpipi[3];
2114
2116 D0 = tKKpipi.AMP(pdau,1);
2117
2118 pdau[0]=-km_kkpipi[0]; pdau[1]=-km_kkpipi[1]; pdau[2]=-km_kkpipi[2]; pdau[3]=km_kkpipi[3];
2119 pdau[4]=-kp_kkpipi[0]; pdau[5]=-kp_kkpipi[1]; pdau[6]=-kp_kkpipi[2]; pdau[7]=kp_kkpipi[3];
2120 pdau[8]=-pim_kkpipi[0]; pdau[9]=-pim_kkpipi[1]; pdau[10]=-pim_kkpipi[2];pdau[11]=pim_kkpipi[3];
2121 pdau[12]=-pip_kkpipi[0];pdau[13]=-pip_kkpipi[1];pdau[14]=-pip_kkpipi[2];pdau[15]=pip_kkpipi[3];
2122 D0bar = tKKpipi.AMP(pdau,1);
2123
2124 if( charm == 1){
2125
2126 r2D0 = std::norm(D0bar)/std::norm(D0);
2127 double temp = std::arg(D0)-std::arg(D0bar);
2128 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
2129 while (temp < -TMath::Pi())
2130 {
2131 temp += 2.0 * TMath::Pi();
2132 }
2133 while (temp > TMath::Pi())
2134 {
2135 temp -= 2.0 * TMath::Pi();
2136 }
2137 deltaD0 = temp;
2138 }
2139 else {
2140 r2D0 = std::norm(D0)/std::norm(D0bar);
2141 double temp = std::arg(D0bar)-std::arg(D0);
2142 while (temp < -TMath::Pi())
2143 {
2144 temp += 2.0 * TMath::Pi();
2145 }
2146 while (temp > TMath::Pi())
2147 {
2148 temp -= 2.0 * TMath::Pi();
2149 }
2150 deltaD0 = temp;
2151 }
2152
2154 {
2155 decay_mode = kSLPlus ;
2156 ++m_nSL ;
2157
2158 r2D0 = 0. ;
2159 deltaD0 = 0. ;
2160 }
2162 {
2163 decay_mode = kSLMinus ;
2164 ++m_nSL ;
2165
2166 r2D0 = 0. ;
2167 deltaD0 = 0. ;
2168 }
2169 //Check CP even modes
2170 else if(
2171 ( nDaughters == 2 &&
2173 nUFP0 == 2 ||
2174 num[ kNeutralScalar ] == 2 ||
2175 ( nUFP0 == 1 && nUFV0 == 1 ) ||
2177 num[ kK0S ] == 2 ||
2178 num[ kK0L ] == 2 ||
2179 ( num[ kK0L ] == 1 && nUFP0 == 1 ) ||
2181 ( num[ kK0L ] == 1 && nUFV0 == 1 ) ||
2184 ( nDaughters == 3 &&
2186 ( nCPPlusEig == 1 && num[ kPi0 ] == 2 ) ||
2189 ( num[ kK0S ] == 1 && nUFP0 == 2 )
2190 ) ) ||
2191 ( nDaughters == 4 &&
2192 ( num[ kK0L ] == 1 && nUFP0 == 3 )||
2194 )
2195 {
2196 decay_mode = kCPPlus ;
2197 ++m_nCPPlus ;
2198
2199 r2D0 = 1. ;
2200 deltaD0 = 0.;
2201 }
2203
2204 decay_mode = kCPPlus ;
2205 ++m_nCPPlus ;
2206
2207 FCP = 0.973;
2208 r2D0 = 1. ;
2209 deltaD0 = 0.;
2210 }
2211 //Check CP odd modes
2212 else if(
2213 ( nDaughters == 2 &&
2214 ( ( num[ kK0S ] == 1 && nUFP0 == 1 ) ||
2216 ( num[ kK0S ] == 1 && nUFV0 == 1 ) ||
2217 ( num[ kNeutralScalar ] == 1 && nUFV0 == 1 ) ||
2218 ( nUFP0 == 1 && num[ kNeutralScalar ] == 1 ) ) ) ||
2219 ( nDaughters == 3 &&
2221 ( num[ kPi0 ] == 3 ) ||
2224 ( num[ kK0L ] == 1 && nUFP0 == 2 ) ) ) ||
2225 ( nDaughters == 4 &&
2226 ( num[ kK0S ] == 1 && num[ kPi0 ] == 3 )|| ( ( num[ kK0S ] == 1 && num[ kK0L ] == 1 ) && nUFP0 == 2 ) )
2227 )
2228 {
2229 decay_mode = kCPMinus ;
2230 ++m_nCPMinus ;
2231
2232 r2D0 = 1. ;
2233 deltaD0 = PI ;
2234 }
2235 else if( nStrange == nAntiStrange + 1 )
2236 {
2237 if(m_debug)cout<< nDaughters <<endl;
2241 if( charm == 1 )
2242 {
2243 if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 2 && num[ kPiMinus ] == 1 ) && nDaughters == 4 ){
2244 if(m_debug)cout<< "True K-3Pi "<<endl;
2245 decay_mode = kFlavoredCF_3 ;
2246 ++m_nFlavoredCFD0_3 ;
2247
2248 complex<double> D0, D0bar;
2249 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2250 double pdau[16];
2251 pi1.clear();pi2.clear();pi3.clear();
2252 kp.clear();
2253
2254 HepLorentzVector kp_kpipipi = (kMinus_kkpipi.at(0));
2255 HepLorentzVector pi1_kpipipi = (piPlus_4pi.at(0));
2256 HepLorentzVector pi2_kpipipi = (piPlus_4pi.at(1));
2257 HepLorentzVector pi3_kpipipi = (piMinus_4pi.at(0));
2258
2259 kp.push_back( kp_kpipipi[0]); kp.push_back(kp_kpipipi[1]); kp.push_back(kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2260 pi1.push_back(pi1_kpipipi[0]);pi1.push_back(pi1_kpipipi[1]);pi1.push_back(pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2261 pi2.push_back(pi2_kpipipi[0]);pi2.push_back(pi2_kpipipi[1]);pi2.push_back(pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2262 pi3.push_back(pi3_kpipipi[0]);pi3.push_back(pi3_kpipipi[1]);pi3.push_back(pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2263
2264 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2265 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2266 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2267 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2268
2270 D0 = tKpipipi.AMP(pdau,1);
2271
2273 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2274 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2275 deltaD0 = ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2276 r2D0 = std::norm(D0bar)/std::norm(D0);
2277
2278 }else if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2279 decay_mode = kFlavoredCF_1 ;
2280 ++m_nFlavoredCFD0_1 ;
2281 r2D0 = pow(m_rCF_1,2) ;
2282 deltaD0 = m_deltaCF_1 ;
2283 RD0 = m_RCF_1;
2284 }else{
2285 //if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2286 if(m_debug)cout<< "True K-Pi "<<endl;
2287 decay_mode = kFlavoredCF_0 ;
2288 ++m_nFlavoredCFD0_0 ;
2289 r2D0 = pow(m_rCF_0,2) ;
2290 deltaD0 = m_deltaCF_0 ;
2291 RD0 = m_RCF_0;
2292 }
2293 }
2294 else
2295 {
2296 if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 2 && num[ kPiMinus ] == 1 ) && nDaughters == 4 ){
2297 decay_mode = kFlavoredCF_3 ;
2298 ++m_nFlavoredDCSD0_3 ;
2299
2300 complex<double> D0, D0bar;
2301 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2302 double pdau[16];
2303 pi1.clear();pi2.clear();pi3.clear();
2304 kp.clear();
2305
2306 HepLorentzVector kp_kpipipi = (kMinus_kkpipi.at(0));
2307 HepLorentzVector pi1_kpipipi = (piPlus_4pi.at(0));
2308 HepLorentzVector pi2_kpipipi = (piPlus_4pi.at(1));
2309 HepLorentzVector pi3_kpipipi = (piMinus_4pi.at(0));
2310
2311 kp.push_back( kp_kpipipi[0]); kp.push_back(kp_kpipipi[1]); kp.push_back(kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2312 pi1.push_back(pi1_kpipipi[0]);pi1.push_back(pi1_kpipipi[1]);pi1.push_back(pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2313 pi2.push_back(pi2_kpipipi[0]);pi2.push_back(pi2_kpipipi[1]);pi2.push_back(pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2314 pi3.push_back(pi3_kpipipi[0]);pi3.push_back(pi3_kpipipi[1]);pi3.push_back(pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2315
2316 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2317 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2318 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2319 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2320
2322 D0 = tKpipipi.AMP(pdau,1);
2323
2325 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2326 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2327 deltaD0 = - ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2328 r2D0 = std::norm(D0)/std::norm(D0bar);
2329
2330 }else if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2331 decay_mode = kFlavoredCF_1 ;
2332 ++m_nFlavoredDCSD0_1 ;
2333 r2D0 = 1. / pow( m_rCF_1,2 ) ;
2334 deltaD0 = -m_deltaCF_1 ;
2335 RD0 = m_RCF_1;
2336 }else{
2337 //if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 ) && nDaughters == 2 )
2338 decay_mode = kFlavoredCF_0 ;
2339 ++m_nFlavoredDCSD0_0 ;
2340 r2D0 = 1. / pow( m_rCF_0,2 ) ;
2341 deltaD0 = -m_deltaCF_0 ;
2342 RD0 = m_RCF_0;
2343 }
2344 }
2345 }
2346 else if( nAntiStrange == nStrange + 1 )
2347 {
2348 if( charm == -1 )
2349 {
2350 if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 2 && num[ kPiPlus ] == 1 ) && nDaughters == 4 ){
2351 decay_mode = kFlavoredCFBar_3 ;
2352 ++m_nFlavoredCFD0_3 ;
2353
2354 complex<double> D0, D0bar;
2355 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2356 double pdau[16];
2357 pi1.clear();pi2.clear();pi3.clear();
2358 kp.clear();
2359
2360 HepLorentzVector kp_kpipipi = (kPlus_kkpipi.at(0));
2361 HepLorentzVector pi1_kpipipi = (piMinus_4pi.at(0));
2362 HepLorentzVector pi2_kpipipi = (piMinus_4pi.at(1));
2363 HepLorentzVector pi3_kpipipi = (piPlus_4pi.at(0));
2364
2365 kp.push_back( -kp_kpipipi[0]); kp.push_back(-kp_kpipipi[1]); kp.push_back(-kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2366 pi1.push_back(-pi1_kpipipi[0]);pi1.push_back(-pi1_kpipipi[1]);pi1.push_back(-pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2367 pi2.push_back(-pi2_kpipipi[0]);pi2.push_back(-pi2_kpipipi[1]);pi2.push_back(-pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2368 pi3.push_back(-pi3_kpipipi[0]);pi3.push_back(-pi3_kpipipi[1]);pi3.push_back(-pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2369
2370 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2371 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2372 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2373 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2374
2376 D0 = tKpipipi.AMP(pdau,1);
2377
2379 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2380 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2381 deltaD0 = ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2382 r2D0 = std::norm(D0bar)/std::norm(D0);
2383
2384
2385 }else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2386 decay_mode = kFlavoredCFBar_1 ;
2387 ++m_nFlavoredCFD0_1 ;
2388 r2D0 = pow(m_rCF_1,2) ;
2389 deltaD0 = m_deltaCF_1 ;
2390 RD0 = m_RCF_1;
2391 }else{
2392 //else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2393 decay_mode = kFlavoredCFBar_0 ;
2394 ++m_nFlavoredCFD0_0 ;
2395 r2D0 = pow(m_rCF_0,2) ;
2396 deltaD0 = m_deltaCF_0 ;
2397 RD0 = m_RCF_0;
2398 }
2399 }
2400 else
2401 {
2402 if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 2 && num[ kPiPlus ] == 1 ) && nDaughters == 4 ){
2403 decay_mode = kFlavoredCFBar_3 ;
2404 ++m_nFlavoredDCSD0_3 ;
2405
2406 complex<double> D0, D0bar;
2407 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2408 double pdau[16];
2409 pi1.clear();pi2.clear();pi3.clear();
2410 kp.clear();
2411
2412 HepLorentzVector kp_kpipipi = (kPlus_kkpipi.at(0));
2413 HepLorentzVector pi1_kpipipi = (piMinus_4pi.at(0));
2414 HepLorentzVector pi2_kpipipi = (piMinus_4pi.at(1));
2415 HepLorentzVector pi3_kpipipi = (piPlus_4pi.at(0));
2416
2417 kp.push_back( -kp_kpipipi[0]); kp.push_back(-kp_kpipipi[1]); kp.push_back(-kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]);//k-
2418 pi1.push_back(-pi1_kpipipi[0]);pi1.push_back(-pi1_kpipipi[1]);pi1.push_back(-pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2419 pi2.push_back(-pi2_kpipipi[0]);pi2.push_back(-pi2_kpipipi[1]);pi2.push_back(-pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2420 pi3.push_back(-pi3_kpipipi[0]);pi3.push_back(-pi3_kpipipi[1]);pi3.push_back(-pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2421
2422 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2423 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2424 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2425 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2426
2428 D0 = tKpipipi.AMP(pdau,1);
2429
2431 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2432 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2433 deltaD0 = - ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2434 r2D0 = std::norm(D0)/std::norm(D0bar);
2435
2436 }else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2437 decay_mode = kFlavoredCFBar_1 ;
2438 ++m_nFlavoredDCSD0_1 ;
2439 r2D0 = 1. / pow( m_rCF_1,2 ) ;
2440 deltaD0 = -m_deltaCF_1 ;
2441 RD0 = m_RCF_1;
2442 }else{
2443 //if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2444 decay_mode = kFlavoredCFBar_0 ;
2445 ++m_nFlavoredDCSD0_0 ;
2446 r2D0 = 1. / pow( m_rCF_0,2 ) ;
2447 deltaD0 = -m_deltaCF_0 ;
2448 RD0 = m_RCF_0;
2449 }
2450 }
2451 }
2452 else if( nStrange == nAntiStrange ) // Unflavored or Cabibbo-suppressed
2453 {
2457 ( nK0 > 0 &&
2458 nFV0 != nFV0Bar &&
2459 nK0 == nFV0Bar ) ||
2460 ( num[ kPiPlus ] > 0 &&
2462 {
2463 decay_mode = kFlavoredCS ;
2464 ++m_nFlavoredCSD0 ;
2465
2466 if( charm == 1 )
2467 {
2468 r2D0 = pow( m_rCS, 2 ) ;
2469 deltaD0 = m_deltaCS ;
2470 }
2471 else
2472 {
2473 r2D0 = 1. / pow( m_rCS, 2 ) ;
2474 deltaD0 = -m_deltaCS ;
2475 }
2476 }
2477 else if( ( num[ kKMinus ] > 0 && ( num[ kKMinus ] == num[ kFVPlus ] || num[ kKMinus ] == nFV0 ) ) || // for K*0 K- pi+
2478 ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0 ) ||
2480 {
2481 decay_mode = kFlavoredCSBar ;
2482 ++m_nFlavoredCSD0 ;
2483
2484 if( charm == -1 )
2485 {
2486 r2D0 = pow( m_rCS, 2 ) ;
2487 deltaD0 = m_deltaCS ;
2488 }
2489 else
2490 {
2491 r2D0 = 1. / pow( m_rCS, 2 ) ;
2492 deltaD0 = -m_deltaCS ;
2493 }
2494 }
2495
2496 // Self-conjugate mixed-CP final states -- pick CP or non-CP at
2497 // random. If CP, pick + or - at random. If non-CP, pick
2498 // flavored or flavoredBar according to the appropriate value of r.
2501 nUFV0 == nDaughters ||
2505 // for rho+rho-; no a1+a1- and rhoa1 final states in DECAY.DEC
2506 )
2507 {
2508 log << MSG::DEBUG <<" [ Self-conjugate mixed-CP ]" << endmsg ;
2509
2510 if( RandFlat::shoot() > 0.5 )
2511 {
2512 if( RandFlat::shoot() > 0.5 )
2513 {
2514 decay_mode = kCPPlus ;
2515 ++m_nCPPlus ;
2516
2517 r2D0 = 1. ;
2518 deltaD0 = PI;
2519 }
2520 else
2521 {
2522 decay_mode = kCPMinus ;
2523 ++m_nCPMinus ;
2524
2525 r2D0 = 1. ;
2526 deltaD0 = 0. ;
2527 }
2528 }
2529 else
2530 {
2531 // Odd # of K0S or K0L = CF; even # of K0S or K0L = SC.
2532 if( nK0 % 2 )
2533 {
2534 // Nov 2007: correct for r2 -> RWS and BR != A^2
2536
2537 if( charm == -1 )
2538 {
2539 cutoff = 1. - cutoff ;
2540 }
2541
2542 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
2543
2544 decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF_0 : kFlavoredCFBar_0 ;
2545
2546 if( ( decay_mode == kFlavoredCF_0 && charm == 1 ) ||
2547 ( decay_mode == kFlavoredCFBar_0 && charm == -1 ) )
2548 {
2549 ++m_nFlavoredCFD0_0 ;
2550
2551 r2D0 = sqrt( m_rCF_0 ) ;
2552 deltaD0 = m_deltaCF_0 ;
2553 }
2554 else
2555 {
2556 ++m_nFlavoredDCSD0_0 ;
2557
2558 r2D0 = 1. / sqrt( m_rCF_0 ) ;
2559 deltaD0 = -m_deltaCF_0 ;
2560 }
2561 }
2562 else
2563 {
2564 // Nov 2007: correct for r2 -> RWS and BR != A^2
2566
2567 if( charm == -1 )
2568 {
2569 cutoff = 1. - cutoff ;
2570 }
2571
2572 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
2573
2574 decay_mode = ( RandFlat::shoot() > cutoff ) ?
2575 kFlavoredCS : kFlavoredCSBar ;
2576 ++m_nFlavoredCSD0 ;
2577
2578 if( ( decay_mode == kFlavoredCS && charm == 1 ) ||
2579 ( decay_mode == kFlavoredCSBar && charm == -1 ) )
2580 {
2581 r2D0 = sqrt( m_rCS ) ;
2582 deltaD0 = m_deltaCS ;
2583 }
2584 else
2585 {
2586 r2D0 = 1. / sqrt( m_rCS ) ;
2587 deltaD0 = -m_deltaCS ;
2588 }
2589 }
2590 }
2591 }
2592 }
2593
2595 {
2596 cout << "decay mode " << decay_mode << endl ;
2597 cout << "D #" << charm << endl ;
2624 }
2625
2626
2627 d0_info[0]=decay_mode;
2628 d0_info[1]=r2D0;
2629 d0_info[2]=deltaD0;
2630 d0_info[3]=RD0;
2631 d0_info[4]=FCP;
2632 d0_info[5]=mnm_gen;
2633 d0_info[6]=mpn_gen;
2634 d0_info[7]= double (isKsPiPi);
2635 d0_info[8]=mpm_gen;
2636 return d0_info;
2637}
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmpi0
Definition: RRes.h:32
Definition: D0To2pip2pim.h:9
complex< double > Amp(vector< double > Pip1, vector< double > Pim1, vector< double > Pip2, vector< double > Pim2)
Definition: D0To2pip2pim.cxx:1426
Definition: D0ToKKpipi.h:9
complex< double > AMP(double const *x0, const int &x1)
Definition: D0ToKKpipi.cxx:2656
Definition: D0ToKLpipi2018.h:9
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition: D0ToKLpipi2018.cxx:92
Definition: D0ToKSpipi2018.h:9
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition: D0ToKSpipi2018.cxx:94
Definition: D0ToKSpipipi0.h:10
complex< double > Amp(vector< double > ks0, vector< double > pip, vector< double > pim, vector< double > pi0)
Definition: D0ToKSpipipi0.cxx:1615
Definition: D0ToKpipipiLHCb.h:8
std::complex< double > AMP(const double *x0, const int &x1)
Definition: D0ToKpipipiLHCb.cxx:2976
Definition: D0TopiKpipi.h:9
complex< double > AMP(double const *x0, const int &x1)
Definition: D0TopiKpipi.cxx:1086
Definition: D0Topippim2pi0.h:9
complex< double > Amp(vector< double > Pip, vector< double > Pim, vector< double > Pi01, vector< double > Pi02)
Definition: D0Topippim2pi0.cxx:1384
Referenced by execute().
◆ initialize()
| StatusCode QCMCFilter::initialize | ( | ) |
Definition at line 272 of file QCMCFilter.cxx.
272 {
273 MsgStream log(msgSvc(), name());
274
275 log << MSG::INFO << "in initialize()" << endmsg;
276
277 // Get the Particle Properties Service
278 IPartPropSvc* p_PartPropSvc;
279 static const bool CREATEIFNOTTHERE(true);
280 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
281 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc)
282 {
283 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
284 return PartPropStatus;
285 }
286 m_particleTable = p_PartPropSvc->PDT();
287
288
289 //Calculates parameters here based on Parameter Input by User
290 // (e.g. run expensive algorithms that are based on parameters
291 // specified by user at run-time)
292
293 //m_y = (m_y - 0.008)/0.292;
294 //corrects the y input so that we get the desired measurable y from CP:Semi-leptonic doubletag method
295 // This equation is a result of solving equations for y dependence.
296 //If parent MC Branching ratio's this equation will need resolved. January 2013 -DA
299 double rCF_0 = m_rCF_0 ;
300 double dCF_0 = m_dCF_0 ;
301 double RCF_0 = m_RCF_0 ;
302 //double zCF_0 = m_zCF_0 ;
304 double rCF2_0 = rCF_0 * rCF_0 ;
305 double zCF2_0 = zCF_0 * zCF_0 ;
306 double rCF_1 = m_rCF_1 ;
307 double dCF_1 = m_dCF_1 ;
308 double RCF_1 = m_RCF_1 ;
309 //double zCF_1 = m_zCF_1 ;
311 double rCF2_1 = rCF_1 * rCF_1 ;
312 double zCF2_1 = zCF_1 * zCF_1 ;
313 double rCF_3 = m_rCF_3 ;
314 double dCF_3 = m_dCF_3 ;
315 double RCF_3 = m_RCF_3 ;
316 //double zCF_3 = m_zCF_3 ;
318 double rCF2_3 = rCF_3 * rCF_3 ;
319 double zCF2_3 = zCF_3 * zCF_3 ;
320 double rCS = m_rCS ;
321 double zCS = m_zCS ;
322 double rCS2 = rCS * rCS ;
323 double zCS2 = zCS * zCS ;
324
329 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
330 double vCF_0CSPlus = ( zCF_0 * zCS + wCF_0 * wCS ) / 2. ;
331 double vCF_1CSPlus = ( zCF_1 * zCS + wCF_1 * wCS ) / 2. ;
332 double vCF_3CSPlus = ( zCF_3 * zCS + wCF_3 * wCS ) / 2. ;
333 double vCF_0CSMinus = ( zCF_0 * zCS - wCF_0 * wCS ) / 2. ;
334 double vCF_1CSMinus = ( zCF_1 * zCS - wCF_1 * wCS ) / 2. ;
335 double vCF_3CSMinus = ( zCF_3 * zCS - wCF_3 * wCS ) / 2. ;
336
337 m_deltaCF_0 = acos( zCF_0 / 2. ) * ( m_wCFSign_0 ? 1. : -1. ) ;
338 m_deltaCF_1 = acos( zCF_1 / 2. ) * ( m_wCFSign_1 ? 1. : -1. ) ;
339 m_deltaCF_3 = acos( zCF_3 / 2. ) * ( m_wCFSign_3 ? 1. : -1. ) ;
340 m_deltaCS = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. ) ;
341
342 double rCF = m_rCF ;
343 double zCF = m_zCF ;
344 double rCF2 = rCF * rCF ;
345 double zCF2 = zCF * zCF ;
346 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
347 double vCFCSPlus = ( zCF * zCS + wCF * wCS ) / 2. ;
348 double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2. ;
349 m_deltaCF = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. ) ;
350
351 double bCF_0 = 1. + 0.5 * rCF_0 * ( zCF_0 * y + wCF_0 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
352 double bCF_1 = 1. + 0.5 * rCF_1 * ( zCF_1 * y + wCF_1 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
353 double bCF_3 = 1. + 0.5 * rCF_3 * ( zCF_3 * y + wCF_3 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
357
358 if( !m_useNewWeights )
359 {
360 bCF_0 = 1. ;
361 bCF_1 = 1. ;
362 bCF_3 = 1. ;
363 m_rwsCF_0 = rCF2_0 ;
364 m_rwsCF_1 = rCF2_1 ;
365 m_rwsCF_3 = rCF2_3 ;
366 bCS = 1. ;
367 m_rwsCS = rCS2 ;
368 bCPPlus = 1. ;
369 bCPMinus = 1. ;
370 }
371
372 m_rwsCF_0 = ( rCF2_0 + 0.5 * rCF_0 * ( zCF_0 * y - wCF_0 * x ) + rmix ) / bCF_0;//though division by b not mentioned in memo it is correct
373 m_rwsCF_1 = ( rCF2_1 + 0.5 * rCF_1 * ( zCF_1 * y - wCF_1 * x ) + rmix ) / bCF_1;//though division by b not mentioned in memo it is correct
374 m_rwsCF_3 = ( rCF2_3 + 0.5 * rCF_3 * ( zCF_3 * y - wCF_3 * x ) + rmix ) / bCF_3;//though division by b not mentioned in memo it is correct
376
377 double bfCF_0 = ( 1. - RCF_0 * rCF_0 * ( 0.5*zCF_0*y + 0.5*wCF_0*x ) + 0.5 * ( -1*x*x + y*y) );//B(D0->F)=AF*AF*brCF_0
380 double bfCFbar_0 = ( rCF2_0 - RCF_0 * rCF_0 * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + 0.5 * ( x*x + y*y) );//B(D0->Fbar)=AF*AF*brCFbar_0
381 double bfCFbar_1 = ( rCF2_1 - RCF_1 * rCF_1 * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) + 0.5 * ( x*x + y*y) );
382 double bfCFbar_3 = ( rCF2_3 - RCF_3 * rCF_3 * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) + 0.5 * ( x*x + y*y) );
383
384
385 // Calculate reweighting factors, normalized to the largest factor.
386 // Note:A fake Dalitz weight is added in to just to make initial predictions
387 // These Dalitz weights are not used in actual code
388
389 // FlavoredCF
390
391 // First calculate rate factors.
392 // m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] = ( ( 1 - RCF_0 * RCF_0 ) / ( 1 - RCF_0 * (1./rCF_0) * (0.5*zCF_0*y - 0.5*wCF_0*x) + rmix/(2*rCF2_0) ) ) ;
393 // rmix * m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;
394 // m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_0 ] = ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( brCF_0 * brCF_0 );
395 // 1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
396
397 // Nov 2007: correct for r2 -> RWS and BR != A^2
398 // m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
399 // m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
400 // ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
401 // m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
402 // ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
403 // m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
404 // ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
405 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] = ( ( 1 - RCF_0 * RCF_0 ) / ( 1 - RCF_0 * (1./rCF_0) * (0.5*zCF_0*y - 0.5*wCF_0*x) + rmix/(2*rCF2_0) ) ) ;
406 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_0 ] = ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( bfCF_0 * bfCF_0 );
407 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_1 ] = ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - 2. * ( rCF2_0 / rCF2_1 ) * RCF_0 * RCF_1 * cos( dCF_0 - dCF_1 ) ) / ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - RCF_1 * (1./rCF_1) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) - (RCF_0*rCF_1/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) );
408 cout<<( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - RCF_1 * (1./rCF_1) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) - (RCF_0*rCF_1/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) )<<endl;
409 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_1 ] = ( 1. + rCF2_0 * rCF2_1 - 2. * RCF_0 * RCF_1 * rCF_0 * rCF_1 ) / ( bfCF_0 * bfCF_1 );
410 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_3 ] = ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) - 2. * ( rCF2_0 / rCF2_3 ) * RCF_0 * RCF_3 * cos( dCF_0 - dCF_3 ) ) / ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) - RCF_3 * (1./rCF_3) * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) - (RCF_0/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) );
411 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_0 * rCF2_3 - 2. * RCF_0 * RCF_3 * rCF_0 * rCF_3 ) / ( bfCF_0 * bfCF_3 );
412 m_weights[ kFlavoredCF_0 ][ kFlavoredCSBar ] =
413 1. + rCF2_0 * rCS2 - rCF_0 * rCS * vCFCSMinus ;
414 m_weights[ kFlavoredCF_0 ][ kFlavoredCS ] =
415 rCF2_0 + rCS2 - rCF_0 * rCS * vCF_0CSPlus +
416 rmix * m_weights[ kFlavoredCF_0 ][ kFlavoredCSBar ] ;
417
420
421 //m_weights[ kFlavoredCF_0 ][ kCPPlus ] =
422 // ( 1. + rCF2_0 + rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPPlus ;
423 //m_weights[ kFlavoredCF_0 ][ kCPMinus ] =
424 // ( 1. + rCF2_0 - rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPMinus ;
425 m_weights[ kFlavoredCF_0 ][ kCPPlus ] =
426 ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) ) / ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) * y + y * y ) * bCPPlus );
427 m_weights[ kFlavoredCF_0 ][ kCPMinus ] =
428 ( 1. + rCF2_0 + 2. * rCF_0 * RCF_0 * cos(dCF_0) ) / ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) * y + y * y ) * bCPMinus );
429 m_weights[ kFlavoredCF_0 ][ kDalitz ] = 1. / bfCF_0;
430 m_weights[ kFlavoredCF_0 ][ k2Pip2Pim ] = 1. / bfCF_0;
431 m_weights[ kFlavoredCF_0 ][ kPipPim2Pi0 ] = 1. / bfCF_0;
432 m_weights[ kFlavoredCF_0 ][ kK0PiPiPi0 ] = 1. / bfCF_0;
433 m_weights[ kFlavoredCF_0 ][ kKKPiPi ] = 1. / bfCF_0;
434
435 // FlavoredCFBar
436 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_0 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] ;
437 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCF_1 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_1 ] ;
438 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_1 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_1 ] ;
439 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCF_3 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_3 ] ;
440 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_3 ] ;
443 m_weights[ kFlavoredCFBar_0 ][ kSLPlus ] = 1. / bCF_0 ;
444 m_weights[ kFlavoredCFBar_0 ][ kSLMinus ] = 1. / bCF_0 ;
447 m_weights[ kFlavoredCFBar_0 ][ kDalitz ] = 1. / bCF_0;
448 m_weights[ kFlavoredCFBar_0 ][ k2Pip2Pim ] = 1. / bCF_0;
449 m_weights[ kFlavoredCFBar_0 ][ kPipPim2Pi0 ] = 1. / bCF_0;
450 m_weights[ kFlavoredCFBar_0 ][ kK0PiPiPi0 ] = 1. / bCF_0;
451 m_weights[ kFlavoredCFBar_0 ][ kKKPiPi ] = 1. / bCF_0;
452
453 m_weights[ kFlavoredCF_1 ][ kFlavoredCF_1 ] = ( ( 1 - RCF_1 * RCF_1 ) / ( 1 - RCF_1 * (1./rCF_1) * (0.5*zCF_1*y - 0.5*wCF_1*x) + rmix/(2*rCF2_1) ) ) ;
454 m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_1 ] = ( 1. + rCF2_1 * rCF2_1 - 2. * RCF_1 * RCF_1 * rCF2_1 ) / ( bfCF_1 * bfCF_1 );
455 m_weights[ kFlavoredCF_1 ][ kFlavoredCF_3 ] = ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) - 2. * ( rCF2_1 / rCF2_3 ) * RCF_1 * RCF_3 * cos( dCF_1 - dCF_3 ) ) / ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) - RCF_3 * (1./rCF_3) * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) - (RCF_1/(rCF_1*rCF_1)) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) + rmix/(rCF_1*rCF_1) );;
456 m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_1 * rCF2_3 - 2. * RCF_1 * RCF_3 * rCF_1 * rCF_3 ) / ( bfCF_1 * bfCF_3 );;
457 m_weights[ kFlavoredCF_1 ][ kFlavoredCSBar ] =
458 1. + rCF2_1 * rCS2 - rCF_1 * rCS * vCFCSMinus ;
459 m_weights[ kFlavoredCF_1 ][ kFlavoredCS ] =
460 rCF2_1 + rCS2 - rCF_1 * rCS * vCF_1CSPlus +
461 rmix * m_weights[ kFlavoredCF_1 ][ kFlavoredCSBar ] ;
462
463 m_weights[ kFlavoredCF_1 ][ kSLPlus ] = 1. / bfCF_1 ;
464 m_weights[ kFlavoredCF_1 ][ kSLMinus ] = 1. / bfCF_1 ;
465
466 // m_weights[ kFlavoredCF_1 ][ kCPPlus ] =
467 // ( 1. + rCF2_1 + rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPPlus ;
468 // m_weights[ kFlavoredCF_1 ][ kCPMinus ] =
469 // ( 1. + rCF2_1 - rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPMinus ;
470 m_weights[ kFlavoredCF_1 ][ kCPPlus ] =
471 ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) ) / ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) *y + y * y ) * bCPPlus );
472 m_weights[ kFlavoredCF_1 ][ kCPMinus ] =
473 ( 1. + rCF2_1 + 2. * rCF_1 * RCF_1 * cos(dCF_1) ) / ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) *y + y * y ) * bCPMinus );
474 m_weights[ kFlavoredCF_1 ][ kDalitz ] = 1. / bCF_1;
475 m_weights[ kFlavoredCF_1 ][ k2Pip2Pim ] = 1. / bCF_1;
476 m_weights[ kFlavoredCF_1 ][ kPipPim2Pi0 ] = 1. / bCF_1;
477 m_weights[ kFlavoredCF_1 ][ kK0PiPiPi0 ] = 1. / bCF_1;
478 m_weights[ kFlavoredCF_1 ][ kKKPiPi ] = 1. / bCF_1;
479
480 // FlavoredCFBar
481 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCFBar_1 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCF_1 ] ;
482 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCF_3 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_3 ] ;
483 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCF_3 ] ;
486 m_weights[ kFlavoredCFBar_1 ][ kSLPlus ] = 1. / bCF_1 ;
487 m_weights[ kFlavoredCFBar_1 ][ kSLMinus ] = 1. / bCF_1 ;
490 m_weights[ kFlavoredCFBar_1 ][ kDalitz ] = 1. / bCF_1;
491 m_weights[ kFlavoredCFBar_1 ][ k2Pip2Pim ] = 1. / bCF_1;
492 m_weights[ kFlavoredCFBar_1 ][ kPipPim2Pi0 ] = 1. / bCF_1;
493 m_weights[ kFlavoredCFBar_1 ][ kK0PiPiPi0 ] = 1. / bCF_1;
494 m_weights[ kFlavoredCFBar_1 ][ kKKPiPi ] = 1. / bCF_1;
495
496 m_weights[ kFlavoredCF_3 ][ kFlavoredCF_3 ] = ( ( 1 - RCF_3 * RCF_3 ) / ( 1 - RCF_3 * (1./rCF_3) * (0.5*zCF_3*y - 0.5*wCF_3*x) + rmix/(2*rCF2_3) ) ) ;
497 m_weights[ kFlavoredCF_3 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_3 * rCF2_3 - 2. * RCF_3 * RCF_3 * rCF2_3 ) / ( bfCF_3 * bfCF_3 );
498 m_weights[ kFlavoredCF_3 ][ kFlavoredCSBar ] =
499 1. + rCF2_3 * rCS2 - rCF_3 * rCS * vCFCSMinus ;
500 m_weights[ kFlavoredCF_3 ][ kFlavoredCS ] =
501 rCF2_3 + rCS2 - rCF_3 * rCS * vCF_3CSPlus +
502 rmix * m_weights[ kFlavoredCF_3 ][ kFlavoredCSBar ] ;
503
504 m_weights[ kFlavoredCF_3 ][ kSLPlus ] = 1. / bfCF_3 ;
505 m_weights[ kFlavoredCF_3 ][ kSLMinus ] = 1. / bfCF_3 ;
506
507 //m_weights[ kFlavoredCF_3 ][ kCPPlus ] =
508 // ( 1. + rCF2_3 + rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPPlus ;
509 //m_weights[ kFlavoredCF_3 ][ kCPMinus ] =
510 // ( 1. + rCF2_3 - rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPMinus ;
511 m_weights[ kFlavoredCF_3 ][ kCPPlus ] =
512 ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) ) / ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) * y + y * y ) * bCPPlus );
513 m_weights[ kFlavoredCF_3 ][ kCPMinus ] =
514 ( 1. + rCF2_3 + 2. * rCF_3 * RCF_3 * cos(dCF_3) ) / ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) * y + y * y ) * bCPMinus );
515 m_weights[ kFlavoredCF_3 ][ kDalitz ] = 1. / bCF_3;
516 m_weights[ kFlavoredCF_3 ][ k2Pip2Pim ] = 1. / bCF_3;
517 m_weights[ kFlavoredCF_3 ][ kPipPim2Pi0 ] = 1. / bCF_3;
518 m_weights[ kFlavoredCF_3 ][ kK0PiPiPi0 ] = 1. / bCF_3;
519 m_weights[ kFlavoredCF_3 ][ kKKPiPi ] = 1. / bCF_3;
520
521 // FlavoredCFBar
522 m_weights[ kFlavoredCFBar_3 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_3 ][ kFlavoredCF_3 ] ;
525 m_weights[ kFlavoredCFBar_3 ][ kSLPlus ] = 1. / bCF_3 ;
526 m_weights[ kFlavoredCFBar_3 ][ kSLMinus ] = 1. / bCF_3 ;
529 m_weights[ kFlavoredCFBar_3 ][ kDalitz ] = 1. / bCF_3;
530 m_weights[ kFlavoredCFBar_3 ][ k2Pip2Pim ] = 1. / bCF_3;
531 m_weights[ kFlavoredCFBar_3 ][ kPipPim2Pi0 ] = 1. / bCF_3;
532 m_weights[ kFlavoredCFBar_3 ][ kK0PiPiPi0 ] = 1. / bCF_3;
533 m_weights[ kFlavoredCFBar_3 ][ kKKPiPi ] = 1. / bCF_3;
534
535 // FlavoredCS
536
537 // First calculate rate factors.
538 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2 ;
540
541 // Nov 2007: correct for r2 -> RWS and BR != A^2
544
545 m_weights[ kFlavoredCS ][ kSLPlus ] = 1. / bCS ;
546 m_weights[ kFlavoredCS ][ kSLMinus ] = 1. / bCS ;
547
548 m_weights[ kFlavoredCS ][ kCPPlus ] = ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus ;
549 m_weights[ kFlavoredCS ][ kCPMinus ] = ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus ;
550
551 m_weights[ kFlavoredCS ][ kDalitz ] = 1. / bCS;
552 m_weights[ kFlavoredCS ][ k2Pip2Pim ] = 1. / bCS;
553 m_weights[ kFlavoredCS ][ kPipPim2Pi0 ] = 1. / bCS;
554 m_weights[ kFlavoredCS ][ kK0PiPiPi0 ] = 1. / bCS;
555 m_weights[ kFlavoredCS ][ kKKPiPi ] = 1. / bCS;
556
557 // FlavoredCSBar
558
560 m_weights[ kFlavoredCSBar ][ kSLPlus ] = 1. / bCS ;
561 m_weights[ kFlavoredCSBar ][ kSLMinus ] = 1. / bCS ;
564 m_weights[ kFlavoredCSBar ][ kDalitz ] = 1. / bCS;
565 m_weights[ kFlavoredCSBar ][ k2Pip2Pim ] = 1. / bCS;
566 m_weights[ kFlavoredCSBar ][ kPipPim2Pi0 ] = 1. / bCS;
567 m_weights[ kFlavoredCSBar ][ kK0PiPiPi0 ] = 1. / bCS;
568 m_weights[ kFlavoredCSBar ][ kKKPiPi ] = 1. / bCS;
569
570 // SLPlus
571
572 // SLPlus/SLPlus should have rate factor = Rmix, but none of these are
573 // generated in uncorrelated MC.
574 m_weights[ kSLPlus ][ kSLPlus ] = 0. ;
575 m_weights[ kSLPlus ][ kSLMinus ] = 1. ;
576 m_weights[ kSLPlus ][ kCPPlus ] = 1. / bCPPlus ;
577 m_weights[ kSLPlus ][ kCPMinus ] = 1. / bCPMinus ;
578 m_weights[ kSLPlus ][ kDalitz ] = 1. ;
579 m_weights[ kSLPlus ][ k2Pip2Pim ] = 1. ;
580 m_weights[ kSLPlus ][ kPipPim2Pi0 ] = 1. ;
581 m_weights[ kSLPlus ][ kK0PiPiPi0 ] = 1. ;
582 m_weights[ kSLPlus ][ kKKPiPi ] = 1. ;
583
584 // SLMinus
585
586 m_weights[ kSLMinus ][ kSLMinus ] = 0. ;
587 m_weights[ kSLMinus ][ kCPPlus ] = 1. / bCPPlus ;
588 m_weights[ kSLMinus ][ kCPMinus ] = 1. / bCPMinus ;
589 m_weights[ kSLMinus ][ kDalitz ] = 1. ;
590 m_weights[ kSLMinus ][ k2Pip2Pim ] = 1. ;
591 m_weights[ kSLMinus ][ kPipPim2Pi0 ] = 1. ;
592 m_weights[ kSLMinus ][ kK0PiPiPi0 ] = 1. ;
593 m_weights[ kSLMinus ][ kKKPiPi ] = 1. ;
594
595 // CPPlus
596
597 m_weights[ kCPPlus ][ kCPPlus ] = 0. ;
598 m_weights[ kCPPlus ][ kCPMinus ] = 2. / bCPPlus / bCPMinus ;
599 m_weights[ kCPPlus ][ kDalitz ] = 1. / bCPPlus;
600 m_weights[ kCPPlus ][ k2Pip2Pim ] = 1. / bCPPlus;
601 m_weights[ kCPPlus ][ kPipPim2Pi0 ] = 1. / bCPPlus;
602 m_weights[ kCPPlus ][ kK0PiPiPi0 ] = 1. / bCPPlus;
603 m_weights[ kCPPlus ][ kKKPiPi ] = 1. / bCPPlus;
604
605 // CPMinus
606
607 m_weights[ kCPMinus ][ kCPMinus ] = 0. ;
608 m_weights[ kCPMinus ][ kDalitz ] = 1. / bCPMinus;
609 m_weights[ kCPMinus ][ k2Pip2Pim ] = 1. / bCPMinus;
610 m_weights[ kCPMinus ][ kPipPim2Pi0 ] = 1. / bCPMinus;
611 m_weights[ kCPMinus ][ kK0PiPiPi0 ] = 1. / bCPMinus;
612 m_weights[ kCPMinus ][ kKKPiPi ] = 1. / bCPMinus;
613
614 //Dalitz estimate
615 m_weights[ kDalitz ][ kDalitz ] = 1;
616 m_weights[ kDalitz ][ k2Pip2Pim ] = 1;
617 m_weights[ kDalitz ][ kPipPim2Pi0 ] = 1;
618 m_weights[ kDalitz ][ kK0PiPiPi0 ] = 1;
619 m_weights[ kDalitz ][ kKKPiPi ] = 1;
620
621 m_weights[ k2Pip2Pim ][ k2Pip2Pim ] = 1;
622 m_weights[ k2Pip2Pim ][ kPipPim2Pi0 ] = 1;
623 m_weights[ k2Pip2Pim ][ kK0PiPiPi0 ] = 1;
624 m_weights[ k2Pip2Pim ][ kKKPiPi ] = 1;
625
626 m_weights[ kPipPim2Pi0 ][ kPipPim2Pi0 ] = 1;
627 m_weights[ kPipPim2Pi0 ][ kK0PiPiPi0 ] = 1;
628 m_weights[ kPipPim2Pi0 ][ kKKPiPi ] = 1;
629
630 m_weights[ kK0PiPiPi0 ][ kK0PiPiPi0 ] = 1;
631 m_weights[ kK0PiPiPi0 ][ kKKPiPi ] = 1;
632
633 m_weights[ kKKPiPi ][ kKKPiPi ] = 1;
635
636 //Boss 6.6.2 MC Branching fractions
637 HepVector fractionsD0( kNDecayTypes+1, 0 ) ;
638 fractionsD0[ kFlavoredCF_0 ] = 0.305 ;
639 fractionsD0[ kFlavoredCFBar_0 ] = 0.0076 ;
640 fractionsD0[ kFlavoredCF_1 ] = 0.144 ;
641 fractionsD0[ kFlavoredCFBar_1 ] = 0.00031 ;
642 fractionsD0[ kFlavoredCF_3 ] = 0.0822 ;
643 fractionsD0[ kFlavoredCFBar_3 ] = 0.00027 ;
644 fractionsD0[ kFlavoredCS ] = 0.031 ;
645 fractionsD0[ kFlavoredCSBar ] = 0.031 ;
646 fractionsD0[ kSLPlus ] = 0.125 ;
647 fractionsD0[ kSLMinus ] = 0. ;
648 fractionsD0[ kCPPlus ] = 0.113 ;
649 fractionsD0[ kCPMinus ] = 0.102 ;
650 fractionsD0[ kDalitz ] = 0.056 ;
651 fractionsD0[ k2Pip2Pim ] = 0.00720 ;
652 fractionsD0[ kPipPim2Pi0 ] = 0.0102 ;
653 fractionsD0[ kK0PiPiPi0 ] = 0.104 ;
654 fractionsD0[ kKKPiPi ] = 0.00273 ;
655 HepVector scales = m_weights * fractionsD0 ;
656 log << MSG::INFO << "Integrated scale factors " <<scales << endmsg ;
657
658
659 HepVector fractionsD0bar( kNDecayTypes+1, 0 ) ;
660 fractionsD0bar[ kFlavoredCF_0 ] = 0.0076 ;
661 fractionsD0bar[ kFlavoredCFBar_0 ] = 0.305 ;
662 fractionsD0bar[ kFlavoredCF_1 ] = 0.00031 ;
663 fractionsD0bar[ kFlavoredCFBar_1 ] = 0.144 ;
664 fractionsD0bar[ kFlavoredCF_3 ] = 0.00027 ;
665 fractionsD0bar[ kFlavoredCFBar_3 ] = 0.0822 ;
666 fractionsD0bar[ kFlavoredCS ] = 0.031 ;
667 fractionsD0bar[ kFlavoredCSBar ] = 0.031 ;
668 fractionsD0bar[ kSLPlus ] = 0. ;
669 fractionsD0bar[ kSLMinus ] = 0.125 ;
670 fractionsD0bar[ kCPPlus ] = 0.113 ;
671 fractionsD0bar[ kCPMinus ] = 0.102 ;
672 fractionsD0bar[ kDalitz ] = 0.056 ;
673 fractionsD0bar[ k2Pip2Pim ] = 0.00720 ;
674 fractionsD0bar[ kPipPim2Pi0 ] = 0.0102 ;
675 fractionsD0bar[ kK0PiPiPi0 ] = 0.104 ;
676 fractionsD0bar[ kKKPiPi ] = 0.00273 ;
677 HepVector weight_norm = fractionsD0bar.T() * scales;
678 log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg ;
679
680
681
682 //Original MC
683 //Branching fraction predictions after weight
684 double brCF_0 = ( fractionsD0bar[ kFlavoredCFBar_0 ] * scales[ kFlavoredCFBar_0 ] ) / weight_norm[0] ;
685 double brCF_1 = ( fractionsD0bar[ kFlavoredCFBar_1 ] * scales[ kFlavoredCFBar_1 ] ) / weight_norm[0] ;
686 double brCF_3 = ( fractionsD0bar[ kFlavoredCFBar_3 ] * scales[ kFlavoredCFBar_3 ] ) / weight_norm[0] ;
687 double brCS = ( fractionsD0bar[ kFlavoredCS ] * scales[ kFlavoredCS ] + fractionsD0bar[ kFlavoredCSBar ] * scales[ kFlavoredCSBar ] ) / weight_norm[0] ;
688 double brDCS_0 = ( fractionsD0bar[ kFlavoredCF_0 ] * scales[ kFlavoredCF_0 ] ) / weight_norm[0] ;
689 double brDCS_1 = ( fractionsD0bar[ kFlavoredCF_1 ] * scales[ kFlavoredCF_1 ] ) / weight_norm[0] ;
690 double brDCS_3 = ( fractionsD0bar[ kFlavoredCF_3 ] * scales[ kFlavoredCF_3 ] ) / weight_norm[0] ;
691 double brCPPlus = ( fractionsD0bar[ kCPPlus ] * scales[ kCPPlus ] ) / weight_norm[0] ;
692 double brCPMinus = ( fractionsD0bar[ kCPMinus ] * scales[ kCPMinus ] ) / weight_norm[0] ;
693
694 //y=0.0 dalitz example used for prediction
695 double dalitz_y_num = -0.000177536;
696 double dalitz_y_den = -0.0150846;
697
698 double numFactCF_0 =
699 -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x ) ;
700 double numFactCF_1 =
701 -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x ) ;
702 double numFactCF_3 =
703 -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x ) ;
704 double numFactCS =
705 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
706 double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0 ;
707 double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1 ;
708 double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3 ;
709 double denFactCS = 0.5 * rCS2 * zCS2 ;
710
711 double num_pre =
712 brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 + brCF_3 * numFactCF_3 + 0.5 * brCS * numFactCS + dalitz_y_num;
713 double den_pre =
714 1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 - brCF_3 * denFactCF_3 - 0.5 * brCS * denFactCS + dalitz_y_den ;
715 double y_pre = num_pre / den_pre ;
717 fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] + fractionsD0bar[ kFlavoredCFBar_0 ] * numFactCF_0 + fractionsD0bar[ kFlavoredCFBar_1 ] * numFactCF_1 + fractionsD0bar[ kFlavoredCFBar_3 ] * numFactCF_3 + fractionsD0bar[ kFlavoredCS ] * numFactCS + dalitz_y_num;
718 double den =
719 1. - fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] - fractionsD0bar[ kFlavoredCFBar_0 ] * denFactCF_0 - fractionsD0bar[ kFlavoredCFBar_1 ] * denFactCF_1 - fractionsD0bar[ kFlavoredCFBar_3 ] * denFactCF_3 - fractionsD0bar[ kFlavoredCS ] * denFactCS + dalitz_y_den;
721 double y_test1 = 0.25 * ( ( ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) / ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) ) -
722 ( ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) / ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) ) ) ;
723
724 log << MSG::INFO << "An Input Y of " << m_y << endmsg ;
725 log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg ;
726 log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg ;
727 log << MSG::INFO << "The CP-SL double tag method predicts y of " <<y_test1 << endmsg ;
728
729 // Renormalize
730 m_largestWeight = 0. ;
731 for( int i = 0 ; i < kNDecayTypes ; ++i )
732 {
733 for( int j = 0 ; j < kNDecayTypes ; ++j )
734 {
735 if( m_weights[ i ][ j ] < 0 ) m_weights[ i ][ j ] = 0; //gets rid of negative weights in early calculations
737 {
738 m_largestWeight = m_weights[ i ][ j ] ;
739 }
740 }
741 }
742 m_weights /= m_largestWeight ;
743
745
746
747 //Set up weight parameters
749 double D0Sum[ 10 ], D0barSum[ 10 ] ; // last index is for sum over bins
750 TComplex D0D0barSum[ 10 ] ;
751 int points[ 10 ] ;
752
753 for( int i = 0 ; i < 10 ; ++i )
754 {
755 D0Sum[ i ] = 0. ;
756 D0barSum[ i ] = 0. ;
757 D0D0barSum[ i ] = TComplex( 0., 0. ) ;
758 points[ i ] = 0 ;
759 }
760
761 double m2min = 0. ;
762 double m2max = 3. ;
763 int nsteps = 1000 ;
764 double stepsize = ( m2max - m2min ) / ( double ) nsteps ;
765
766 for( int i = 0 ; i < nsteps ; ++i )
767 {
768 double m2i = m2min + ( ( double ) i + 0.5 ) * stepsize ;
769
770 for( int j = 0 ; j < nsteps ; ++j )
771 {
772 double m2j = m2min + ( ( double ) j + 0.5 ) * stepsize ;
773
775 {
776 double m2k = 1.86450*1.86450 + 0.497648*0.497648 +
777 0.139570*0.139570 + 0.139570*0.139570 - m2i - m2j ;
778
779 TComplex D0, D0bar;
780 if (m_dalitzModel == 0) { //Cleo model
781 D0 = t.CLEO_Amplitude( m2i, m2j, m2k ) ;
782 D0bar = t.CLEO_Amplitude( m2j, m2i, m2k ) ;}
783 if (m_dalitzModel == 1) { //Babar model
784 D0 = t.Babar_Amplitude( m2i, m2j, m2k ) ;
785 D0bar = t.Babar_Amplitude( m2j, m2i, m2k ) ;}
786 if (m_dalitzModel == 2) { //Belle model
787 D0 = t.Amplitude( m2i, m2j, m2k ) ;
788 D0bar = t.Amplitude( m2j, m2i, m2k ) ;}
789
790 if ( j <= i ) {// lower half only
793
794 ++points[ bin ] ;
795 D0Sum[ bin ] += norm( D0 ) ;
796 D0barSum[ bin ] += norm( D0bar ) ;
798
799 ++points[ 9 ] ;
800 D0Sum[ 9 ] += norm( D0 ) ;
801 D0barSum[ 9 ] += norm( D0bar ) ;
802 D0D0barSum[ 9 ] += D0 * conj( D0bar ) ;
803 }
804
805 }
806 }
807 }
808
809 for( int i = 0 ; i < 10 ; ++i )
810 {
811 double r2 = D0barSum[ i ] / D0Sum[ i ] ;
814
815 cout << "BIN " << i << ": "
816 << points[ i ] << " pts "
818 << endmsg ;
819 }
820
821 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
822 return StatusCode::SUCCESS;
823}
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
Definition: Dalitz.h:30
The documentation for this class was generated from the following files:
