Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtSSDCP.cpp
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) 2001      Caltech
10 //
11 // Module: EvtSSDCP.cc
12 //
13 // Description: See EvtSSDCP.hh
14 //
15 // Modification history:
16 //
17 //    RYD       August 12, 2001        Module created
18 //    F. Sandrelli, Fernando M-V  March 1, 2002     Debugged and added z parameter (CPT violation)
19 //------------------------------------------------------------------------
20 // 
21 #include "EvtGenBase/EvtPatches.hh"
22 #include <stdlib.h>
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtRandom.hh"
25 #include "EvtGenBase/EvtGenKine.hh"
26 #include "EvtGenBase/EvtCPUtil.hh"
27 #include "EvtGenBase/EvtPDL.hh"
28 #include "EvtGenBase/EvtReport.hh"
29 #include "EvtGenBase/EvtVector4C.hh"
30 #include "EvtGenBase/EvtTensor4C.hh"
31 #include "EvtGenModels/EvtSSDCP.hh"
32 #include "EvtGenBase/EvtIncoherentMixing.hh"
33 #include <string>
34 #include "EvtGenBase/EvtConst.hh"
35 using std::endl;
36
37 EvtSSDCP::~EvtSSDCP() {}
38
39 std::string EvtSSDCP::getName(){
40
41   return "SSD_CP";     
42
43 }
44
45
46 EvtDecayBase* EvtSSDCP::clone(){
47
48   return new EvtSSDCP;
49
50 }
51
52 void EvtSSDCP::init(){
53
54   // check that there are 8 or 12 or 14 arguments
55
56   checkNArg(14,12,8);
57   checkNDaug(2);
58
59   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
60   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
61
62   // Check it is a B0 or B0s
63   if ( ( getParentId() != EvtPDL::getId( "B0" ) )
64        && ( getParentId() != EvtPDL::getId( "anti-B0" ) ) 
65        && ( getParentId() != EvtPDL::getId( "B_s0" ) )
66        && ( getParentId() != EvtPDL::getId( "anti-B_s0" ) ) ) {
67     report( ERROR , "EvtGen" ) << "EvtSSDCP only decays B0 and B0s" 
68                                << std::endl ;
69     ::abort() ;
70   }
71
72   
73   if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
74        (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||(d2type==EvtSpinType::TENSOR)))||
75        (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||(d1type==EvtSpinType::TENSOR)))
76        ) {
77     report(ERROR,"EvtGen") << "EvtSSDCP generator expected "
78                            << "one of the daugters to be a scalar, the other either scalar, vector, or tensor, found:"
79                            << EvtPDL::name(getDaug(0)).c_str()<<" and "<<EvtPDL::name(getDaug(1)).c_str()<<endl;
80     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
81     ::abort();
82   }
83
84   _dm=getArg(0)/EvtConst::c; //units of 1/mm
85
86   _dgog=getArg(1);
87
88   
89   _qoverp=getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
90   _poverq=1.0/_qoverp;
91
92   _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));
93
94   _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7)));
95   
96   if (getNArg()>=12){
97     _eigenstate=false;
98     _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
99     _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
100   }
101   else{
102     //I'm somewhat confused about this. For a CP eigenstate set the
103     //amplitudes to the same. For a non CP eigenstate CPT invariance
104     //is enforced. (ryd)
105     if (
106         (getDaug(0)==EvtPDL::chargeConj(getDaug(0))&&
107          getDaug(1)==EvtPDL::chargeConj(getDaug(1)))||
108         (getDaug(0)==EvtPDL::chargeConj(getDaug(1))&&
109          getDaug(1)==EvtPDL::chargeConj(getDaug(0)))){
110       _eigenstate=true;
111     }else{
112       _eigenstate=false;
113       _A_fbar=conj(_Abar_f);
114       _Abar_fbar=conj(_A_f);
115     }
116   }
117
118   //FS: new check for z 
119   if (getNArg()==14){ //FS Set _z parameter if provided else set it 0
120     _z=EvtComplex(getArg(12),getArg(13));
121   }
122   else{
123     _z=EvtComplex(0.0,0.0);
124   }
125
126   // FS substituted next 2 lines...
127
128
129   //
130   //  _gamma=EvtPDL::getctau(EvtPDL::getId("B0"));  //units of 1/mm
131   //_dgamma=_gamma*0.5*_dgog;
132   //
133   // ...with:
134
135   if ( ( getParentId() == EvtPDL::getId("B0") ) || 
136        ( getParentId() == EvtPDL::getId("anti-B0") ) ) {
137     _gamma=1./EvtPDL::getctau(EvtPDL::getId("B0")); //gamma/c (1/mm)
138   }
139   else{
140     _gamma=1./EvtPDL::getctau(EvtPDL::getId("B_s0")) ;
141   }
142
143   _dgamma=_gamma*_dgog;  //dgamma/c (1/mm) 
144
145   if (verbose()){
146     report(INFO,"EvtGen") << "SSD_CP will generate CP/CPT violation:"
147                           << endl << endl
148                           << "    " << EvtPDL::name(getParentId()).c_str() << " --> "
149                           << EvtPDL::name(getDaug(0)).c_str() << " + "
150                           << EvtPDL::name(getDaug(1)).c_str() << endl << endl
151                           << "using parameters:" << endl << endl
152                           << "  delta(m)  = " << _dm << " hbar/ps" << endl
153                           << "dGamma      = "  << _dgamma <<" ps-1" <<endl
154                           << "       q/p  = " << _qoverp << endl  
155                           << "        z  = " << _z << endl  
156                           << "       tau  = " << 1./_gamma << " ps" << endl;
157   }
158
159 }
160
161 void EvtSSDCP::initProbMax() {
162   double theProbMax = 
163     abs(_A_f) * abs(_A_f) +
164     abs(_Abar_f) * abs(_Abar_f) +
165     abs(_A_fbar) * abs(_A_fbar) +
166     abs(_Abar_fbar) * abs(_Abar_fbar);
167
168   if (_eigenstate) theProbMax*=2;
169
170   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
171   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
172   if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
173
174   setProbMax(theProbMax);
175 }
176
177 void EvtSSDCP::decay( EvtParticle *p){
178
179   static EvtId B0=EvtPDL::getId("B0");
180   static EvtId B0B=EvtPDL::getId("anti-B0");
181   
182   static EvtId B0s = EvtPDL::getId("B_s0");
183   static EvtId B0Bs = EvtPDL::getId("anti-B_s0");
184
185   double t;
186   EvtId other_b;
187   EvtId daugs[2];
188
189   int flip=0;
190   if (!_eigenstate){
191     if (EvtRandom::Flat(0.0,1.0)<0.5) flip=1;
192   }
193
194   if (!flip) {
195     daugs[0]=getDaug(0);
196     daugs[1]=getDaug(1);
197   }
198   else{
199     daugs[0]=EvtPDL::chargeConj(getDaug(0));
200     daugs[1]=EvtPDL::chargeConj(getDaug(1));
201   }
202
203   EvtParticle *d;
204   p->initializePhaseSpace(2, daugs);
205
206   EvtComplex amp;
207
208
209   EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5); // t is c*Dt (mm)
210 //  EvtIncoherentMixing::OtherB( p , t , other_b , 0.5 ) ;
211   
212   //if (flip) t=-t;
213
214   //FS We assume DGamma=GammaLow-GammaHeavy and Dm=mHeavy-mLow
215   EvtComplex expH=exp(-EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
216   EvtComplex expL=exp(EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
217   //FS Definition of gp and gm
218   EvtComplex gp=0.5*(expL+expH);
219   EvtComplex gm=0.5*(expL-expH);
220   //FS Calculation os sqrt(1-z^2) 
221   EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
222   
223   //EvtComplex BB=0.5*(expL+expH);                  // <B0|B0(t)>
224   //EvtComplex barBB=_qoverp*0.5*(expL-expH);       // <B0bar|B0(t)>
225   //EvtComplex BbarB=_poverq*0.5*(expL-expH);       // <B0|B0bar(t)>
226   //EvtComplex barBbarB=BB;                         // <B0bar|B0bar(t)>
227   //  FS redefinition of these guys... (See BAD #188 eq.35 for ref.)
228   //  q/p is taken as in the BaBar Phys. Book (opposite sign wrt ref.)
229   EvtComplex BB=gp+_z*gm;                 // <B0|B0(t)>
230   EvtComplex barBB=sqz*_qoverp*gm;       // <B0bar|B0(t)>
231   EvtComplex BbarB=sqz*_poverq*gm;       // <B0|B0bar(t)>
232   EvtComplex barBbarB=gp-_z*gm;           // <B0bar|B0bar(t)>
233
234   if (!flip){
235     if (other_b==B0B||other_b==B0Bs){
236       //at t=0 we have a B0
237       //report(INFO,"EvtGen") << "B0B"<<endl;
238       amp=BB*_A_f+barBB*_Abar_f;
239       //std::cout << "noflip B0B tag:"<<amp<<std::endl;
240       //amp=0.0;
241     }
242     if (other_b==B0||other_b==B0s){
243       //report(INFO,"EvtGen") << "B0"<<endl;
244       amp=BbarB*_A_f+barBbarB*_Abar_f;
245     }
246   }else{
247     if (other_b==B0||other_b==B0s){
248       amp=BbarB*_A_fbar+barBbarB*_Abar_fbar;
249       //std::cout << "flip B0 tag:"<<amp<<std::endl;
250       //amp=0.0;
251     }
252     if (other_b==B0B||other_b==B0Bs){
253       amp=BB*_A_fbar+barBB*_Abar_fbar;
254     }
255   }
256
257
258   EvtVector4R p4_parent=p->getP4Restframe();
259   double m_parent=p4_parent.mass();
260
261   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
262
263   EvtVector4R momv;
264   EvtVector4R moms;
265
266   if (d2type==EvtSpinType::SCALAR){
267     d2type=EvtPDL::getSpinType(getDaug(0));
268     d= p->getDaug(0);
269     momv = d->getP4();
270     moms = p->getDaug(1)->getP4();
271   }
272   else{
273     d= p->getDaug(1);
274     momv = d->getP4();
275     moms = p->getDaug(0)->getP4();
276   }
277
278
279
280   if (d2type==EvtSpinType::SCALAR) {
281     vertex(amp);
282   }
283   
284   if (d2type==EvtSpinType::VECTOR) {
285     
286     double norm=momv.mass()/(momv.d3mag()*p->mass());
287
288     //std::cout << amp << " " << norm << " " << p4_parent << d->getP4()<< std::endl;
289     //    std::cout << EvtPDL::name(d->getId()) << " " << EvtPDL::name(p->getDaug(0)->getId()) << 
290     //  " 1and2 " << EvtPDL::name(p->getDaug(1)->getId()) << std::endl;
291     //std::cout << d->eps(0) << std::endl;
292     //std::cout << d->epsParent(0) << std::endl;
293     vertex(0,amp*norm*p4_parent*(d->epsParent(0)));
294     vertex(1,amp*norm*p4_parent*(d->epsParent(1)));
295     vertex(2,amp*norm*p4_parent*(d->epsParent(2)));
296   
297   }
298
299   if (d2type==EvtSpinType::TENSOR) {
300
301     double norm=d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
302  
303    
304    vertex(0,amp*norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
305    vertex(1,amp*norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
306    vertex(2,amp*norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
307    vertex(3,amp*norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
308    vertex(4,amp*norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);
309   
310   }
311
312
313   return ;
314 }
315
316 std::string EvtSSDCP::getParamName(int i) {
317   switch(i) {
318   case 0:
319     return "deltaM";
320   case 1:
321     return "deltaGammaOverGamma";
322   case 2:
323     return "qOverP";
324   case 3:
325     return "qOverPPhase";
326   case 4:
327     return "Af";
328   case 5:
329     return "AfPhase";
330   case 6:
331     return "Abarf";
332   case 7:
333     return "AbarfPhase";
334   case 8:
335     return "Afbar";
336   case 9:
337     return "AfbarPhase";
338   case 10:
339     return "Abarfbar";
340   case 11:
341     return "AbarfbarPhase";
342   case 12:
343     return "Z";
344   case 13:
345     return "ZPhase";
346   default:
347     return "";
348   }
349 }
350
351 std::string EvtSSDCP::getParamDefault(int i) {
352   switch(i) {
353   case 3:
354     return "0.0";
355   case 4:
356     return "1.0";
357   case 5:
358     return "0.0";
359   case 6:
360     return "1.0";
361   case 7:
362     return "0.0";
363   default:
364     return "";
365   }
366 }