[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtVSSBMixCPT.cxx

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 2002      INFN-Pisa
//
// Module: EvtVSSBMixCPT.cc
//
// Description:
//    Routine to decay vector-> scalar scalar with coherent BB-like mixing
//                              including CPT effects
//    Based on VSSBMIX
//
// Modification history:
//
//    F. Sandrelli, Fernando M-V March 03, 2002 
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenModels/EvtVSSBMixCPT.hh"
#include "EvtGenBase/EvtId.hh"
#include <string>
#include "EvtGenBase/EvtRandom.hh"
using std::endl;

EvtVSSBMixCPT::~EvtVSSBMixCPT() {}

std::string EvtVSSBMixCPT::getName(){
  return "VSS_BMIX";
}


EvtDecayBase* EvtVSSBMixCPT::clone(){
  return new EvtVSSBMixCPT;
}

void EvtVSSBMixCPT::init(){

  if ( getNArg()>4) checkNArg(14,12,8);

  if (getNArg()<1) {
    report(ERROR,"EvtGen") << "EvtVSSBMix generator expected "
                           << " at least 1 argument (deltam) but found:"<<getNArg()<<endl;
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  }
  // check that we are asked to produced exactly 2 daughters
  //4 are allowed if they are aliased..
  checkNDaug(2,4);

  if ( getNDaug()==4) {
    if ( getDaug(0)!=getDaug(2)||getDaug(1)!=getDaug(3)){
      report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator allows "
			     << " 4 daughters only if 1=3 and 2=4"
			     << " (but 3 and 4 are aliased "<<endl;
      report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
      ::abort();
    }
  }
  // check that we are asked to decay a vector particle into a pair
  // of scalar particles

  checkSpinParent(EvtSpinType::VECTOR);

  checkSpinDaughter(0,EvtSpinType::SCALAR);
  checkSpinDaughter(1,EvtSpinType::SCALAR);

  // check that our daughter particles are charge conjugates of each other
  if(!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
    report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
			   << "to be charge conjugate." << endl
			   << "  Found " << EvtPDL::name(getDaug(0)).c_str() << " and "
			   << EvtPDL::name(getDaug(1)).c_str() << endl;
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  }
  // check that both daughter particles have the same lifetime
  if(EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
    report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
			   << "to have the same lifetime." << endl
			   << "  Found ctau = "
			   << EvtPDL::getctau(getDaug(0)) << " mm and "
			   << EvtPDL::getctau(getDaug(1)) << " mm" << endl;
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  }
  // precompute quantities that will be used to generate events
  // and print out a summary of parameters for this decay

  // mixing frequency in hbar/mm
  _freq= getArg(0)/EvtConst::c; 

  // deltaG 
  double gamma= 1/EvtPDL::getctau(getDaug(0));  // gamma/c (1/mm)
  _dGamma=0.0;
  double dgog=0.0;
  if ( getNArg() > 1 ) {
    dgog=getArg(1);
    _dGamma=dgog*gamma;
  }
  // q/p
  _qoverp = EvtComplex(1.0,0.0);
  if ( getNArg() > 2){
    _qoverp = EvtComplex(getArg(2),0.0); 
  } 
  if ( getNArg() > 3) {
    _qoverp = getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
  }
  _poverq=1.0/_qoverp;

  // decay amplitudes
  _A_f=EvtComplex(1.0,0.0);
  _Abar_f=EvtComplex(0.0,0.0);
  _A_fbar=_Abar_f;  // CPT conservation
  _Abar_fbar=_A_f;  // CPT conservation
  if ( getNArg() > 4){
    _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));    // this allows for DCSD
    _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7))); // this allows for DCSD 
    if ( getNArg() > 8 ){
      // CPT violation in decay 
      _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
      _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
    } else {
      // CPT conservation in decay
      _A_fbar=_Abar_f;
      _Abar_fbar=_A_f;
    }
  }

  // CPT violation in mixing
  _z = EvtComplex(0.0,0.0);
  if ( getNArg() > 12 ){
    _z = EvtComplex(getArg(12),getArg(13)); 
  }


  // some printout
  double tau= 1e12*EvtPDL::getctau(getDaug(0))/EvtConst::c; // in ps
  double dm= 1e-12*getArg(0); // B0/anti-B0 mass difference in hbar/ps
  double x= dm*tau; 
  double y= dgog*0.5; //y=dgamma/(2*gamma) 
  double qop2 = abs(_qoverp*_qoverp);
  _chib0_b0bar=qop2*(x*x+y*y)/(qop2*(x*x+y*y)+2+x*x-y*y);  // does not include CPT in mixing
  _chib0bar_b0=(1/qop2)*(x*x+y*y)/((1/qop2)*(x*x+y*y)+2+x*x-y*y); // does not include CPT in mixing 

  if ( verbose() ) {
    report(INFO,"EvtGen") << "VSS_BMIXCPT will generate mixing and CPT/CP effects in mixing:"
			  << endl << endl
			  << "    " << EvtPDL::name(getParentId()).c_str() << " --> "
			  << EvtPDL::name(getDaug(0)).c_str() << " + "
			  << EvtPDL::name(getDaug(1)).c_str() << endl << endl
			  << "using parameters:" << endl << endl
			  << "  delta(m)  = " << dm << " hbar/ps" << endl
			  << "  _freq     = " << _freq << " hbar/mm" << endl
			  << "  dgog      = "  << dgog <<endl
			  << "  dGamma    = "  << _dGamma <<" hbar/mm" <<endl
			  << "  q/p       = " << _qoverp << endl  
			  << "  z         = " << _z << endl  
			  << "  tau       = " << tau << " ps" << endl
			  << "  x         = " << x << endl
			  << " chi(B0->B0bar) = " << _chib0_b0bar << endl
			  << " chi(B0bar->B0) = " << _chib0bar_b0 << endl 
			  << " Af         = " << _A_f << endl
			  << " Abarf      = " << _Abar_f << endl
			  << " Afbar      = " << _A_fbar << endl
			  << " Abarfbar   = " << _Abar_fbar << endl
			  << endl;
  }
}

void EvtVSSBMixCPT::initProbMax(){
  // this value is ok for reasonable values of all the parameters
  setProbMax(4.0);
}

void EvtVSSBMixCPT::decay( EvtParticle *p ){

  static EvtId B0=EvtPDL::getId("B0");
  static EvtId B0B=EvtPDL::getId("anti-B0");

  // generate a final state according to phase space

  double rndm= EvtRandom::random();

  if ( getNDaug()==4) {
    EvtId tempDaug[2];

    if ( rndm < 0.5 ) { tempDaug[0]=getDaug(0);  tempDaug[1]=getDaug(3); }
    else{ tempDaug[0]=getDaug(2);  tempDaug[1]=getDaug(1); }

    p->initializePhaseSpace(2,tempDaug);
  }
  else{ //nominal case.
    p->initializePhaseSpace(2,getDaugs());
  }

  EvtParticle *s1,*s2;

  s1 = p->getDaug(0);
  s2 = p->getDaug(1);
  //delete any daughters - if there are daughters, they
  //are from the initialization and will be redone later
  if ( s1->getNDaug() > 0 ) { s1->deleteDaughters();}
  if ( s2->getNDaug() > 0 ) { s2->deleteDaughters();}
  
  EvtVector4R p1= s1->getP4();
  EvtVector4R p2= s2->getP4();

  // throw a random number to decide if this final state should be mixed
  rndm= EvtRandom::random();
  int mixed= (rndm < 0.5) ? 1 : 0;

  // if this decay is mixed, choose one of the 2 possible final states
  // with equal probability (re-using the same random number)
  if(mixed==1) {
    EvtId mixedId= (rndm < 0.25) ? getDaug(0) : getDaug(1);
    EvtId mixedId2= mixedId;
    if (getNDaug()==4&&rndm<0.25)  mixedId2=getDaug(2);
    if (getNDaug()==4&&rndm>0.25)  mixedId2=getDaug(3);
    s1->init(mixedId, p1);
    s2->init(mixedId2, p2);
  }


  // if this decay is unmixed, choose one of the 2 possible final states
  // with equal probability (re-using the same random number)
  if(mixed==0) {
    EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
    EvtId unmixedId2= (rndm < 0.75) ? getDaug(1) : getDaug(0);
    if (getNDaug()==4&&rndm<0.75)  unmixedId2=getDaug(3);
    if (getNDaug()==4&&rndm>0.75)  unmixedId2=getDaug(2);
    s1->init(unmixedId, p1);
    s2->init(unmixedId2, p2);
  }

  // choose a decay time for each final state particle using the
  // lifetime (which must be the same for both particles) in pdt.table
  // and calculate the lifetime difference for this event
  s1->setLifetime();
  s2->setLifetime();
  double dct= s1->getLifetime() - s2->getLifetime(); // in mm

  // Convention: _dGamma=GammaLight-GammaHeavy
  EvtComplex exp1(-0.25*_dGamma*dct,0.5*_freq*dct);

  /*
  //Find the flavor of the B that decayed first.
  EvtId firstDec = (dct > 0 ) ? s2->getId() : s1->getId();
 
  //This tags the flavor of the other particle at that time. 
  EvtId stateAtDeltaTeq0 = ( firstDec==B0 ) ? B0B : B0; 
  */
  EvtId stateAtDeltaTeq0 = (s2->getId()==B0) ? B0B : B0;

  // calculate the oscillation amplitude, based on wether this event is mixed or not
  EvtComplex osc_amp;

  //define some useful functions: (see BAD #188 eq. 39 for ref.) 
  EvtComplex gp=0.5*(exp(-1.0*exp1)+exp(exp1)); 
  EvtComplex gm=0.5*(exp(-1.0*exp1)-exp(exp1));
  EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
  
  EvtComplex BB=gp+_z*gm;                // <B0|B0(t)> 
  EvtComplex barBB=-sqz*_qoverp*gm;      // <B0bar|B0(t)>
  EvtComplex BbarB=-sqz*_poverq*gm;      // <B0|B0bar(t)>
  EvtComplex barBbarB=gp-_z*gm;          // <B0bar|B0bar(t)>

  //
  if ( !mixed&&stateAtDeltaTeq0==B0 ) {
    osc_amp= BB*_A_f+barBB*_Abar_f;
  }
  if ( !mixed&&stateAtDeltaTeq0==B0B ) {
    osc_amp= barBbarB*_Abar_fbar+BbarB*_A_fbar;
  }

  if ( mixed&&stateAtDeltaTeq0==B0 ) {
    osc_amp=barBB*_Abar_fbar+BB*_A_fbar;  
  }
  if ( mixed&&stateAtDeltaTeq0==B0B ) {
    osc_amp=BbarB*_A_f+barBbarB*_Abar_f;
  }

  // store the amplitudes for each parent spin basis state
  double norm=1.0/p1.d3mag();
  vertex(0,norm*osc_amp*p1*(p->eps(0)));
  vertex(1,norm*osc_amp*p1*(p->eps(1)));
  vertex(2,norm*osc_amp*p1*(p->eps(2)));

  return ;
}
Commit	Line	Data
da0e9ce3	1	//--------------------------------------------------------------------------
	2	//
	3	// Environment:
	4	// This software is part of the EvtGen package developed jointly
	5	// for the BaBar and CLEO collaborations. If you use all or part
	6	// of it, please give an appropriate acknowledgement.
	7	//
	8	// Copyright Information: See EvtGen/COPYRIGHT
	9	// Copyright (C) 2002 INFN-Pisa
	10	//
	11	// Module: EvtVSSBMixCPT.cc
	12	//
	13	// Description:
	14	// Routine to decay vector-> scalar scalar with coherent BB-like mixing
	15	// including CPT effects
	16	// Based on VSSBMIX
	17	//
	18	// Modification history:
	19	//
	20	// F. Sandrelli, Fernando M-V March 03, 2002
	21	//
	22	//------------------------------------------------------------------------
	23	//
	24	#include "EvtGenBase/EvtPatches.hh"
	25	#include <stdlib.h>
	26	#include "EvtGenBase/EvtConst.hh"
	27	#include "EvtGenBase/EvtParticle.hh"
	28	#include "EvtGenBase/EvtGenKine.hh"
	29	#include "EvtGenBase/EvtPDL.hh"
	30	#include "EvtGenBase/EvtReport.hh"
	31	#include "EvtGenBase/EvtVector4C.hh"
	32	#include "EvtGenModels/EvtVSSBMixCPT.hh"
	33	#include "EvtGenBase/EvtId.hh"
	34	#include <string>
	35	#include "EvtGenBase/EvtRandom.hh"
	36	using std::endl;
	37
	38	EvtVSSBMixCPT::~EvtVSSBMixCPT() {}
	39
	40	std::string EvtVSSBMixCPT::getName(){
	41	return "VSS_BMIX";
	42	}
	43
	44
	45	EvtDecayBase* EvtVSSBMixCPT::clone(){
	46	return new EvtVSSBMixCPT;
	47	}
	48
	49	void EvtVSSBMixCPT::init(){
	50
	51	if ( getNArg()>4) checkNArg(14,12,8);
	52
	53	if (getNArg()<1) {
	54	report(ERROR,"EvtGen") << "EvtVSSBMix generator expected "
	55	<< " at least 1 argument (deltam) but found:"<<getNArg()<<endl;
	56	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
	57	::abort();
	58	}
	59	// check that we are asked to produced exactly 2 daughters
	60	//4 are allowed if they are aliased..
	61	checkNDaug(2,4);
	62
	63	if ( getNDaug()==4) {
	64	if ( getDaug(0)!=getDaug(2)\|\|getDaug(1)!=getDaug(3)){
65	report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator allows "
66	<< " 4 daughters only if 1=3 and 2=4"
67	<< " (but 3 and 4 are aliased "<<endl;
68	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
69	::abort();
70	}
71	}
72	// check that we are asked to decay a vector particle into a pair
73	// of scalar particles
74
75	checkSpinParent(EvtSpinType::VECTOR);
76
77	checkSpinDaughter(0,EvtSpinType::SCALAR);
78	checkSpinDaughter(1,EvtSpinType::SCALAR);
79
80	// check that our daughter particles are charge conjugates of each other
81	if(!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
82	report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
83	<< "to be charge conjugate." << endl
84	<< " Found " << EvtPDL::name(getDaug(0)).c_str() << " and "
85	<< EvtPDL::name(getDaug(1)).c_str() << endl;
86	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
87	::abort();
88	}
89	// check that both daughter particles have the same lifetime
90	if(EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
91	report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
92	<< "to have the same lifetime." << endl
93	<< " Found ctau = "
94	<< EvtPDL::getctau(getDaug(0)) << " mm and "
95	<< EvtPDL::getctau(getDaug(1)) << " mm" << endl;
96	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
97	::abort();
98	}
99	// precompute quantities that will be used to generate events
100	// and print out a summary of parameters for this decay
101
102	// mixing frequency in hbar/mm
103	_freq= getArg(0)/EvtConst::c;
104
105	// deltaG
106	double gamma= 1/EvtPDL::getctau(getDaug(0)); // gamma/c (1/mm)
107	_dGamma=0.0;
108	double dgog=0.0;
109	if ( getNArg() > 1 ) {
110	dgog=getArg(1);
111	_dGamma=dgog*gamma;
112	}
113	// q/p
114	_qoverp = EvtComplex(1.0,0.0);
115	if ( getNArg() > 2){
116	_qoverp = EvtComplex(getArg(2),0.0);
117	}
118	if ( getNArg() > 3) {
119	_qoverp = getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
120	}
121	_poverq=1.0/_qoverp;
122
123	// decay amplitudes
124	_A_f=EvtComplex(1.0,0.0);
125	_Abar_f=EvtComplex(0.0,0.0);
126	_A_fbar=_Abar_f; // CPT conservation
127	_Abar_fbar=_A_f; // CPT conservation
128	if ( getNArg() > 4){
129	_A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5))); // this allows for DCSD
130	_Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7))); // this allows for DCSD
131	if ( getNArg() > 8 ){
132	// CPT violation in decay
133	_A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
134	_Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
135	} else {
136	// CPT conservation in decay
137	_A_fbar=_Abar_f;
138	_Abar_fbar=_A_f;
139	}
140	}
141
142	// CPT violation in mixing
143	_z = EvtComplex(0.0,0.0);
144	if ( getNArg() > 12 ){
145	_z = EvtComplex(getArg(12),getArg(13));
146	}
147
148
149	// some printout
150	double tau= 1e12*EvtPDL::getctau(getDaug(0))/EvtConst::c; // in ps
151	double dm= 1e-12*getArg(0); // B0/anti-B0 mass difference in hbar/ps
152	double x= dm*tau;
153	double y= dgog0.5; //y=dgamma/(2gamma)
154	double qop2 = abs(_qoverp*_qoverp);
155	_chib0_b0bar=qop2(xx+yy)/(qop2(xx+yy)+2+xx-yy); // does not include CPT in mixing
156	_chib0bar_b0=(1/qop2)(xx+yy)/((1/qop2)(xx+yy)+2+xx-yy); // does not include CPT in mixing
157
158	if ( verbose() ) {
159	report(INFO,"EvtGen") << "VSS_BMIXCPT will generate mixing and CPT/CP effects in mixing:"
160	<< endl << endl
161	<< " " << EvtPDL::name(getParentId()).c_str() << " --> "
162	<< EvtPDL::name(getDaug(0)).c_str() << " + "
163	<< EvtPDL::name(getDaug(1)).c_str() << endl << endl
164	<< "using parameters:" << endl << endl
165	<< " delta(m) = " << dm << " hbar/ps" << endl
166	<< " _freq = " << _freq << " hbar/mm" << endl
167	<< " dgog = " << dgog <<endl
168	<< " dGamma = " << _dGamma <<" hbar/mm" <<endl
169	<< " q/p = " << _qoverp << endl
170	<< " z = " << _z << endl
171	<< " tau = " << tau << " ps" << endl
172	<< " x = " << x << endl
173	<< " chi(B0->B0bar) = " << _chib0_b0bar << endl
174	<< " chi(B0bar->B0) = " << _chib0bar_b0 << endl
175	<< " Af = " << _A_f << endl
176	<< " Abarf = " << _Abar_f << endl
177	<< " Afbar = " << _A_fbar << endl
178	<< " Abarfbar = " << _Abar_fbar << endl
179	<< endl;
180	}
181	}
182
183	void EvtVSSBMixCPT::initProbMax(){
184	// this value is ok for reasonable values of all the parameters
185	setProbMax(4.0);
186	}
187
188	void EvtVSSBMixCPT::decay( EvtParticle *p ){
189
190	static EvtId B0=EvtPDL::getId("B0");
191	static EvtId B0B=EvtPDL::getId("anti-B0");
192
193	// generate a final state according to phase space
194
195	double rndm= EvtRandom::random();
196
197	if ( getNDaug()==4) {
198	EvtId tempDaug[2];
199
200	if ( rndm < 0.5 ) { tempDaug[0]=getDaug(0); tempDaug[1]=getDaug(3); }
201	else{ tempDaug[0]=getDaug(2); tempDaug[1]=getDaug(1); }
202
203	p->initializePhaseSpace(2,tempDaug);
204	}
205	else{ //nominal case.
206	p->initializePhaseSpace(2,getDaugs());
207	}
208
209	EvtParticle s1,s2;
210
211	s1 = p->getDaug(0);
212	s2 = p->getDaug(1);
213	//delete any daughters - if there are daughters, they
214	//are from the initialization and will be redone later
215	if ( s1->getNDaug() > 0 ) { s1->deleteDaughters();}
216	if ( s2->getNDaug() > 0 ) { s2->deleteDaughters();}
217
218	EvtVector4R p1= s1->getP4();
219	EvtVector4R p2= s2->getP4();
220
221	// throw a random number to decide if this final state should be mixed
222	rndm= EvtRandom::random();
223	int mixed= (rndm < 0.5) ? 1 : 0;
224
225	// if this decay is mixed, choose one of the 2 possible final states
226	// with equal probability (re-using the same random number)
227	if(mixed==1) {
228	EvtId mixedId= (rndm < 0.25) ? getDaug(0) : getDaug(1);
229	EvtId mixedId2= mixedId;
230	if (getNDaug()==4&&rndm<0.25) mixedId2=getDaug(2);
231	if (getNDaug()==4&&rndm>0.25) mixedId2=getDaug(3);
232	s1->init(mixedId, p1);
233	s2->init(mixedId2, p2);
234	}
235
236
237	// if this decay is unmixed, choose one of the 2 possible final states
238	// with equal probability (re-using the same random number)
239	if(mixed==0) {
240	EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
241	EvtId unmixedId2= (rndm < 0.75) ? getDaug(1) : getDaug(0);
242	if (getNDaug()==4&&rndm<0.75) unmixedId2=getDaug(3);
243	if (getNDaug()==4&&rndm>0.75) unmixedId2=getDaug(2);
244	s1->init(unmixedId, p1);
245	s2->init(unmixedId2, p2);
246	}
247
248	// choose a decay time for each final state particle using the
249	// lifetime (which must be the same for both particles) in pdt.table
250	// and calculate the lifetime difference for this event
251	s1->setLifetime();
252	s2->setLifetime();
253	double dct= s1->getLifetime() - s2->getLifetime(); // in mm
254
255	// Convention: _dGamma=GammaLight-GammaHeavy
256	EvtComplex exp1(-0.25_dGammadct,0.5_freqdct);
257
258	/*
259	//Find the flavor of the B that decayed first.
260	EvtId firstDec = (dct > 0 ) ? s2->getId() : s1->getId();
261
262	//This tags the flavor of the other particle at that time.
263	EvtId stateAtDeltaTeq0 = ( firstDec==B0 ) ? B0B : B0;
264	*/
265	EvtId stateAtDeltaTeq0 = (s2->getId()==B0) ? B0B : B0;
266
267	// calculate the oscillation amplitude, based on wether this event is mixed or not
268	EvtComplex osc_amp;
269
270	//define some useful functions: (see BAD #188 eq. 39 for ref.)
271	EvtComplex gp=0.5(exp(-1.0exp1)+exp(exp1));
272	EvtComplex gm=0.5(exp(-1.0exp1)-exp(exp1));
273	EvtComplex sqz=sqrt(abs(1-_z_z))exp(EvtComplex(0,arg(1-_z*_z)/2));
274
275	EvtComplex BB=gp+_z*gm; // <B0\|B0(t)>
276	EvtComplex barBB=-sqz_qoverpgm; // <B0bar\|B0(t)>
277	EvtComplex BbarB=-sqz_poverqgm; // <B0\|B0bar(t)>
278	EvtComplex barBbarB=gp-_z*gm; // <B0bar\|B0bar(t)>
279
280	//
281	if ( !mixed&&stateAtDeltaTeq0==B0 ) {
282	osc_amp= BB_A_f+barBB_Abar_f;
283	}
284	if ( !mixed&&stateAtDeltaTeq0==B0B ) {
285	osc_amp= barBbarB_Abar_fbar+BbarB_A_fbar;
286	}
287
288	if ( mixed&&stateAtDeltaTeq0==B0 ) {
289	osc_amp=barBB_Abar_fbar+BB_A_fbar;
290	}
291	if ( mixed&&stateAtDeltaTeq0==B0B ) {
292	osc_amp=BbarB_A_f+barBbarB_Abar_f;
293	}
294
295	// store the amplitudes for each parent spin basis state
296	double norm=1.0/p1.d3mag();
297	vertex(0,normosc_ampp1*(p->eps(0)));
298	vertex(1,normosc_ampp1*(p->eps(1)));
299	vertex(2,normosc_ampp1*(p->eps(2)));
300
301	return ;
302	}
303
304
305
306
307
308
309