[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtDecayAmp.cpp

//--------------------------------------------------------------------------
//
// 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) 1998      Caltech, UCSB
//
// Module: EvtGen/EvtDecayAmp.cc
//
// Description: Baseclass for models that calculates amplitudes
//
// Modification history:
//
//    DJL/RYD     August 11, 1998         Module created
//
//------------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"


#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtDecayAmp.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtRadCorr.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtReport.hh"
using std::endl;


void EvtDecayAmp::makeDecay(EvtParticle* p, bool recursive){

  //original default value
  int ntimes=10000;
  
  int more;

  EvtSpinDensity rho;
  double prob,prob_max;

  _amp2.init(p->getId(),getNDaug(),getDaugs());

  do{

    _daugsDecayedByParentModel=false;
    _weight = 1.0;
    decay(p);

    rho=_amp2.getSpinDensity();

    prob=p->getSpinDensityForward().normalizedProb(rho);

    if (prob<0.0) {
      report(ERROR,"EvtGen")<<"Negative prob:"<<p->getId().getId()
			    <<" "<<p->getChannel()<<endl;

      report(ERROR,"EvtGen") << "rho_forward:"<<endl;
      report(ERROR,"EvtGen") << p->getSpinDensityForward();
      report(ERROR,"EvtGen") << "rho decay:"<<endl;
      report(ERROR,"EvtGen") << rho <<endl;
    }

    if (prob!=prob) {

      report(DEBUG,"EvtGen") << "Forward density matrix:"<<endl;
      report(DEBUG,"EvtGen") << p->getSpinDensityForward();

      report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl;
      report(DEBUG,"EvtGen") << rho;

      report(DEBUG,"EvtGen") << "prob:"<<prob<<endl;
      
      report(DEBUG,"EvtGen") << "Particle:"
			     <<EvtPDL::name(p->getId()).c_str()<<endl;
      report(DEBUG,"EvtGen") << "channel        :"<<p->getChannel()<<endl;
      report(DEBUG,"EvtGen") << "Momentum:" << p->getP4() << " " << p->mass() << endl;
      if( p->getParent()!=0){
	report(DEBUG,"EvtGen") << "parent:"
			       <<EvtPDL::name(
				p->getParent()->getId()).c_str()<<endl;
	report(DEBUG,"EvtGen") << "parent channel        :"
			       <<p->getParent()->getChannel()<<endl;

        size_t i;
	report(DEBUG,"EvtGen") << "parent daughters  :";
        for (i=0;i<p->getParent()->getNDaug();i++){
	  report(DEBUG,"") << EvtPDL::name(
			    p->getParent()->getDaug(i)->getId()).c_str()
				 << " ";
        }
	report(DEBUG,"") << endl;

	report(DEBUG,"EvtGen") << "daughters  :";
        for (size_t i=0;i<p->getNDaug();i++){
	  report(DEBUG,"") << EvtPDL::name(
			    p->getDaug(i)->getId()).c_str()
				 << " ";
        }
	report(DEBUG,"") << endl;

	report(DEBUG,"EvtGen") << "daughter momenta  :" << endl;;
        for (size_t i=0;i<p->getNDaug();i++){
	  report(DEBUG,"") << p->getDaug(i)->getP4() << " " << p->getDaug(i)->mass();
	  report(DEBUG,"") << endl;
        }

      }
    }


    prob/=_weight;

    prob_max = getProbMax(prob);
    p->setDecayProb(prob/prob_max);

    more=prob<EvtRandom::Flat(prob_max);
    
    ntimes--;

  }while(ntimes&&more);

  if (ntimes==0){
    report(DEBUG,"EvtGen") << "Tried accept/reject: 10000" 
			   <<" times, and rejected all the times!"<<endl;
   
    report(DEBUG,"EvtGen")<<p->getSpinDensityForward()<<endl;
    report(DEBUG,"EvtGen") << "Is therefore accepting the last event!"<<endl;
    report(DEBUG,"EvtGen") << "Decay of particle:"<<
      EvtPDL::name(p->getId()).c_str()<<"(channel:"<<
      p->getChannel()<<") with mass "<<p->mass()<<endl;
    
    for(size_t ii=0;ii<p->getNDaug();ii++){
      report(DEBUG,"EvtGen") <<"Daughter "<<ii<<":"<<
	EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<<
	p->getDaug(ii)->mass()<<endl;
    }				   
  }

  EvtSpinDensity rho_list[10];

  rho_list[0]=p->getSpinDensityForward();

  EvtAmp ampcont;

  if (_amp2._pstates!=1){
    ampcont=_amp2.contract(0,p->getSpinDensityForward());
  }
  else{
    ampcont=_amp2;
  }


  // it may be that the parent decay model has already
  // done the decay - this should be rare and the
  // model better know what it is doing..
  
  if ( !daugsDecayedByParentModel() ){

    if(recursive) {   
    
      for(size_t i=0;i<p->getNDaug();i++){

	rho.setDim(_amp2.dstates[i]);
      
	if (_amp2.dstates[i]==1) {
	  rho.set(0,0,EvtComplex(1.0,0.0));
	}
	else{
	  rho=ampcont.contract(_amp2._dnontrivial[i],_amp2);
	}
	
	if (!rho.check()) {
	  
	  report(ERROR,"EvtGen") << "-------start error-------"<<endl;
	  report(ERROR,"EvtGen")<<"forward rho failed Check:"<<
	    EvtPDL::name(p->getId()).c_str()<<" "<<p->getChannel()<<" "<<i<<endl;
	  
	  p->printTree();

	  for (size_t idaug = 0; idaug < p->getNDaug(); idaug++) {
	    EvtParticle* daughter = p->getDaug(idaug);
	    if (daughter != 0) {daughter->printTree();}
	  }

	  EvtParticle* pParent = p->getParent();
	  if (pParent != 0) {
	    report(ERROR,"EvtGen")<<"Parent:"<<EvtPDL::name(pParent->getId()).c_str()<<endl;

	    EvtParticle* grandParent = pParent->getParent();

	    if (grandParent != 0) {
	      report(ERROR,"EvtGen")<<"GrandParent:"<<EvtPDL::name(grandParent->getId()).c_str()<<endl;
	    }
	  }

	  report(ERROR,"EvtGen") << " EvtSpinDensity rho: " << rho;
	  
	  _amp2.dump();

	  for(size_t ii=0;ii<i+1;ii++){
	    report(ERROR,"EvtGen") << "rho_list[" << ii << "] = " << rho_list[ii];
	  }

	  report(ERROR,"EvtGen") << "-------Done with error-------"<<endl;  

	}
      
	p->getDaug(i)->setSpinDensityForward(rho);
	p->getDaug(i)->decay();
	
	rho_list[i+1]=p->getDaug(i)->getSpinDensityBackward();
	
	if (_amp2.dstates[i]!=1){
	ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]);
	}
      
	
      }
      
      p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list));
      
      
      if (!p->getSpinDensityBackward().check()) {
	
	report(ERROR,"EvtGen")<<"rho_backward failed Check"<<
	  p->getId().getId()<<" "<<p->getChannel()<<endl;
      
	report(ERROR,"EvtGen") << p->getSpinDensityBackward();
      
      }
    }    
  }


  if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) {
    int n_daug_orig=p->getNDaug();
    EvtRadCorr::doRadCorr(p);
    int n_daug_new=p->getNDaug();
    for (int i=n_daug_orig;i<n_daug_new;i++){
      p->getDaug(i)->decay();
    }
  }

}
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) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtGen/EvtDecayAmp.cc
	12	//
	13	// Description: Baseclass for models that calculates amplitudes
	14	//
	15	// Modification history:
	16	//
	17	// DJL/RYD August 11, 1998 Module created
	18	//
	19	//------------------------------------------------------------------------
	20	#include "EvtGenBase/EvtPatches.hh"
	21
	22
	23
	24	#include "EvtGenBase/EvtDecayBase.hh"
	25	#include "EvtGenBase/EvtDecayAmp.hh"
	26	#include "EvtGenBase/EvtParticle.hh"
	27	#include "EvtGenBase/EvtPDL.hh"
	28	#include "EvtGenBase/EvtRandom.hh"
	29	#include "EvtGenBase/EvtRadCorr.hh"
	30	#include "EvtGenBase/EvtAmp.hh"
	31	#include "EvtGenBase/EvtReport.hh"
	32	using std::endl;
	33
	34
	35	void EvtDecayAmp::makeDecay(EvtParticle* p, bool recursive){
	36
	37	//original default value
	38	int ntimes=10000;
	39
	40	int more;
	41
	42	EvtSpinDensity rho;
	43	double prob,prob_max;
	44
	45	_amp2.init(p->getId(),getNDaug(),getDaugs());
	46
	47	do{
	48
	49	_daugsDecayedByParentModel=false;
	50	_weight = 1.0;
	51	decay(p);
	52
	53	rho=_amp2.getSpinDensity();
	54
	55	prob=p->getSpinDensityForward().normalizedProb(rho);
	56
	57	if (prob<0.0) {
	58	report(ERROR,"EvtGen")<<"Negative prob:"<<p->getId().getId()
	59	<<" "<<p->getChannel()<<endl;
	60
	61	report(ERROR,"EvtGen") << "rho_forward:"<<endl;
	62	report(ERROR,"EvtGen") << p->getSpinDensityForward();
	63	report(ERROR,"EvtGen") << "rho decay:"<<endl;
	64	report(ERROR,"EvtGen") << rho <<endl;
65	}
66
67	if (prob!=prob) {
68
69	report(DEBUG,"EvtGen") << "Forward density matrix:"<<endl;
0ca57c2f	70	report(DEBUG,"EvtGen") << p->getSpinDensityForward();
da0e9ce3	71
da0e9ce3	72	report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl;
0ca57c2f	73	report(DEBUG,"EvtGen") << rho;
da0e9ce3	74
	75	report(DEBUG,"EvtGen") << "prob:"<<prob<<endl;
	76
	77	report(DEBUG,"EvtGen") << "Particle:"
	78	<<EvtPDL::name(p->getId()).c_str()<<endl;
	79	report(DEBUG,"EvtGen") << "channel :"<<p->getChannel()<<endl;
	80	report(DEBUG,"EvtGen") << "Momentum:" << p->getP4() << " " << p->mass() << endl;
	81	if( p->getParent()!=0){
	82	report(DEBUG,"EvtGen") << "parent:"
	83	<<EvtPDL::name(
	84	p->getParent()->getId()).c_str()<<endl;
	85	report(DEBUG,"EvtGen") << "parent channel :"
	86	<<p->getParent()->getChannel()<<endl;
	87
	88	size_t i;
	89	report(DEBUG,"EvtGen") << "parent daughters :";
	90	for (i=0;i<p->getParent()->getNDaug();i++){
	91	report(DEBUG,"") << EvtPDL::name(
	92	p->getParent()->getDaug(i)->getId()).c_str()
	93	<< " ";
	94	}
	95	report(DEBUG,"") << endl;
	96
	97	report(DEBUG,"EvtGen") << "daughters :";
	98	for (size_t i=0;i<p->getNDaug();i++){
	99	report(DEBUG,"") << EvtPDL::name(
	100	p->getDaug(i)->getId()).c_str()
	101	<< " ";
	102	}
	103	report(DEBUG,"") << endl;
	104
	105	report(DEBUG,"EvtGen") << "daughter momenta :" << endl;;
	106	for (size_t i=0;i<p->getNDaug();i++){
	107	report(DEBUG,"") << p->getDaug(i)->getP4() << " " << p->getDaug(i)->mass();
	108	report(DEBUG,"") << endl;
	109	}
	110
	111	}
	112	}
	113
	114
	115	prob/=_weight;
	116
	117	prob_max = getProbMax(prob);
	118	p->setDecayProb(prob/prob_max);
	119
	120	more=prob<EvtRandom::Flat(prob_max);
	121
	122	ntimes--;
	123
	124	}while(ntimes&&more);
	125
	126	if (ntimes==0){
	127	report(DEBUG,"EvtGen") << "Tried accept/reject: 10000"
	128	<<" times, and rejected all the times!"<<endl;
	129
	130	report(DEBUG,"EvtGen")<<p->getSpinDensityForward()<<endl;
	131	report(DEBUG,"EvtGen") << "Is therefore accepting the last event!"<<endl;
	132	report(DEBUG,"EvtGen") << "Decay of particle:"<<
	133	EvtPDL::name(p->getId()).c_str()<<"(channel:"<<
	134	p->getChannel()<<") with mass "<<p->mass()<<endl;
	135
	136	for(size_t ii=0;ii<p->getNDaug();ii++){
	137	report(DEBUG,"EvtGen") <<"Daughter "<<ii<<":"<<
138	EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<<
139	p->getDaug(ii)->mass()<<endl;
140	}
141	}
142
143	EvtSpinDensity rho_list[10];
144
145	rho_list[0]=p->getSpinDensityForward();
146
147	EvtAmp ampcont;
148
149	if (_amp2._pstates!=1){
150	ampcont=_amp2.contract(0,p->getSpinDensityForward());
151	}
152	else{
153	ampcont=_amp2;
154	}
155
156
157	// it may be that the parent decay model has already
158	// done the decay - this should be rare and the
159	// model better know what it is doing..
160
161	if ( !daugsDecayedByParentModel() ){
162
163	if(recursive) {
164
165	for(size_t i=0;i<p->getNDaug();i++){
166
167	rho.setDim(_amp2.dstates[i]);
168
169	if (_amp2.dstates[i]==1) {
170	rho.set(0,0,EvtComplex(1.0,0.0));
171	}
172	else{
173	rho=ampcont.contract(_amp2._dnontrivial[i],_amp2);
174	}
175
176	if (!rho.check()) {
177
178	report(ERROR,"EvtGen") << "-------start error-------"<<endl;
179	report(ERROR,"EvtGen")<<"forward rho failed Check:"<<
180	EvtPDL::name(p->getId()).c_str()<<" "<<p->getChannel()<<" "<<i<<endl;
181
0ca57c2f	182	p->printTree();
	183
	184	for (size_t idaug = 0; idaug < p->getNDaug(); idaug++) {
	185	EvtParticle* daughter = p->getDaug(idaug);
	186	if (daughter != 0) {daughter->printTree();}
	187	}
	188
	189	EvtParticle* pParent = p->getParent();
	190	if (pParent != 0) {
	191	report(ERROR,"EvtGen")<<"Parent:"<<EvtPDL::name(pParent->getId()).c_str()<<endl;
	192
	193	EvtParticle* grandParent = pParent->getParent();
	194
	195	if (grandParent != 0) {
	196	report(ERROR,"EvtGen")<<"GrandParent:"<<EvtPDL::name(grandParent->getId()).c_str()<<endl;
	197	}
	198	}
	199
	200	report(ERROR,"EvtGen") << " EvtSpinDensity rho: " << rho;
da0e9ce3	201
da0e9ce3	202	_amp2.dump();
0ca57c2f	203
da0e9ce3	204	for(size_t ii=0;ii<i+1;ii++){
0ca57c2f	205	report(ERROR,"EvtGen") << "rho_list[" << ii << "] = " << rho_list[ii];
da0e9ce3	206	}
0ca57c2f	207
da0e9ce3	208	report(ERROR,"EvtGen") << "-------Done with error-------"<<endl;
0ca57c2f	209
da0e9ce3	210	}
	211
	212	p->getDaug(i)->setSpinDensityForward(rho);
	213	p->getDaug(i)->decay();
	214
	215	rho_list[i+1]=p->getDaug(i)->getSpinDensityBackward();
	216
	217	if (_amp2.dstates[i]!=1){
	218	ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]);
	219	}
	220
	221
	222	}
	223
	224	p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list));
	225
	226
	227	if (!p->getSpinDensityBackward().check()) {
	228
	229	report(ERROR,"EvtGen")<<"rho_backward failed Check"<<
	230	p->getId().getId()<<" "<<p->getChannel()<<endl;
	231
	232	report(ERROR,"EvtGen") << p->getSpinDensityBackward();
	233
	234	}
	235	}
	236	}
	237
	238
	239	if (getPHOTOS() \|\| EvtRadCorr::alwaysRadCorr()) {
	240	int n_daug_orig=p->getNDaug();
	241	EvtRadCorr::doRadCorr(p);
	242	int n_daug_new=p->getNDaug();
	243	for (int i=n_daug_orig;i<n_daug_new;i++){
	244	p->getDaug(i)->decay();
	245	}
	246	}
	247
	248	}
	249
	250
	251
	252
	253
	254
	255
	256
	257
	258