[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtPHOTOS.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) 1998      Caltech, UCSB
//
// Module: EvtPHOTOS.cc
//
// Description: This routine takes the particle *p and applies
//              the PHOTOS package to generate final state radiation
//              on the produced mesons.
//
// Modification history:
//
//    RYD     October 1, 1997        Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPhotonParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenModels/EvtPHOTOS.hh"
#include "EvtGenBase/EvtReport.hh"
#include "TRandom.h"
#include <stdlib.h>

extern "C" void begevtgenstorex_(int *,int *,int *,int *,
                                int *,int *,int *,int *,
                                double *,double *,double *, 
                                double *,double *,double *, 
                                double *,double *,double *);

extern "C" void begevtgengetx_(int *,int *,int *,int *,
			      int *,int *,int *,int *,
			      double *,double *,double *, 
			      double *,double *,double *, 
			      double *,double *,double *);

extern "C" void heplst_(int *);

extern "C" void photos_(int *);

extern "C" void phoini_(int *, int *);


EvtPHOTOS::EvtPHOTOS(std::string photontype){

    _photontype=photontype;
  
}

void EvtPHOTOS::doRadCorr( EvtParticle *p){

  static int first=1;

  //added by Lange Jan4,2000
  //allow to set photon tupe
  static EvtId GAMM=EvtPDL::getId(_photontype);

  if (GAMM==EvtId(-1,-1)) {
    report(ERROR,"EvtGen") << "In EvtPHOTOS::doRadCorr():Particle:"<<
      _photontype<<" is not in EvtPDL"<<std::endl;
     ::abort();
  }

  if (first) {
    first=0;
    int iseed1=(int)gRandom->Integer(31327);
    int iseed2=(int)gRandom->Integer(30080);
    phoini_(&iseed1,&iseed2);
  }

  double mpho=EvtPDL::getMeanMass(GAMM);

  int entry,eventnum,numparticle,istat,partnum,mother;
  int daugfirst,dauglast;

  int numparticlephotos;

  double px,py,pz,e,m,x,y,z,t;

  static EvtId dq=EvtPDL::getId("d");
  static EvtId adq=EvtPDL::getId("anti-d");
  static EvtId uq=EvtPDL::getId("u");
  static EvtId auq=EvtPDL::getId("anti-u");
  static EvtId sq=EvtPDL::getId("s");
  static EvtId asq=EvtPDL::getId("anti-s");
  static EvtId cq=EvtPDL::getId("c");
  static EvtId acq=EvtPDL::getId("anti-c");
  static EvtId bq=EvtPDL::getId("b");
  static EvtId abq=EvtPDL::getId("anti-b");
  static EvtId tq=EvtPDL::getId("t");
  static EvtId atq=EvtPDL::getId("anti-t");
  static EvtId vpho=EvtPDL::getId("vpho");

  static EvtIdSet quarks(dq,adq,uq,auq,sq,asq,cq,acq,bq,abq,tq,atq);

  if ( p->getId() == vpho ) return;
  if ( p->getNDaug() > 10 ) return;

  px=0.0;
  py=0.0;
  pz=0.0;
  e=p->mass();
  m=p->mass();
  x=0.0;
  y=0.0;
  z=0.0;
  t=0.0;
  
  entry=1;
  eventnum=1;
  numparticle=1;
  istat=2;
  partnum=EvtPDL::getStdHep(p->getId());
  mother=0;
  daugfirst=2;
  dauglast=1+p->getNDaug();

  begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
                  &mother,&daugfirst,&dauglast,
		  &px,&py,&pz,&e,&m,&x,&y,&z,&t);

  //  std::cout << EvtPDL::name(p->getId()) <<" " ;
  for(size_t i=0;i<p->getNDaug();i++){

    //No quarks to photos
    if (quarks.contains(p->getDaug(i)->getId())==1) continue;

    px=p->getDaug(i)->getP4().get(1);
    py=p->getDaug(i)->getP4().get(2);
    pz=p->getDaug(i)->getP4().get(3);
    e=p->getDaug(i)->getP4().get(0);
    m=p->getDaug(i)->mass();
    x=0.0;
    y=0.0;
    z=0.0;
    t=0.0;

    //    std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " ;
    entry+=1;
    eventnum=1;
    numparticle+=1;
    istat=1;
    partnum=EvtPDL::getStdHep(p->getDaug(i)->getId());
    mother=1;
    daugfirst=0;
    dauglast=0;    

    begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
		    &mother,&daugfirst,&dauglast,
		    &px,&py,&pz,&e,&m,&x,&y,&z,&t);
    

  }
  //  std::cout << std::endl;
  //can't use heplst since the common block used by the BaBar
  //implementation of PHOTOS  is renamed due to real*4 vs real*8
  //problems.

  //int mlst=1;

  //heplst_(&mlst);

  entry=1;

  //  report(INFO,"EvtGen") << "Doing photos " << EvtPDL::name(p->getId()) << endl;
  photos_(&entry);
  //  report(INFO,"EvtGen") << "done\n";
  begevtgengetx_(&entry,&eventnum,&numparticlephotos,&istat,&partnum,
		    &mother,&daugfirst,&dauglast,
		    &px,&py,&pz,&e,&m,&x,&y,&z,&t);
    

  //report(INFO,"EvtGen") << "numparticlephotos:"<<numparticlephotos<<endl;
  
  if (numparticle==numparticlephotos) return;

  EvtVector4R new4mom;

  int np;

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

    entry=i+2;

    begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
		    &mother,&daugfirst,&dauglast,
		    &px,&py,&pz,&e,&m,&x,&y,&z,&t);

    //this is needed to ensure that photos does not
    //change the masses. But it will violate energy conservation!
    double mp=p->getDaug(i)->mass();
    e=sqrt(mp*mp+px*px+py*py+pz*pz);
        
    new4mom.set(e,px,py,pz);

    p->getDaug(i)->setP4WithFSR(new4mom);

    
  }

  for(entry=numparticle+1;entry<=numparticlephotos;entry++){

    begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
		    &mother,&daugfirst,&dauglast,
		    &px,&py,&pz,&e,&m,&x,&y,&z,&t);
        
    //Hack here to give the photon the mass of the generated particle
    e=sqrt(mpho*mpho+px*px+py*py+pz*pz);
    new4mom.set(e,px,py,pz);

    //new4mom.dump();

    EvtPhotonParticle* gamma;
    gamma=new EvtPhotonParticle;
    gamma->init(GAMM,new4mom);
    gamma->setFSRP4toZero();
    //    report(INFO,"EvtGen") << gamma << " " << p << " "<< px << " " << py << " " << pz << " " << p->getNDaug() << " " << EvtPDL::name(p->getId())<<" " << entry << " " <<numparticlephotos<<endl;
    gamma->addDaug(p);

//    p->getDaug(i)->set_type(EvtSpinType::PHOTON);

  }
  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) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtPHOTOS.cc
	12	//
	13	// Description: This routine takes the particle *p and applies
	14	// the PHOTOS package to generate final state radiation
	15	// on the produced mesons.
	16	//
	17	// Modification history:
	18	//
	19	// RYD October 1, 1997 Module created
	20	//
	21	//------------------------------------------------------------------------
	22	//
	23	#include "EvtGenBase/EvtPatches.hh"
	24	#include "EvtGenBase/EvtPatches.hh"
	25	#include "EvtGenBase/EvtIdSet.hh"
	26	#include "EvtGenBase/EvtParticle.hh"
	27	#include "EvtGenBase/EvtPhotonParticle.hh"
	28	#include "EvtGenBase/EvtPDL.hh"
	29	#include "EvtGenModels/EvtPHOTOS.hh"
	30	#include "EvtGenBase/EvtReport.hh"
362af318	31	#include "TRandom.h"
da0e9ce3	32	#include <stdlib.h>
	33
	34	extern "C" void begevtgenstorex_(int ,int ,int ,int ,
	35	int ,int ,int ,int ,
	36	double ,double ,double *,
	37	double ,double ,double *,
	38	double ,double ,double *);
	39
	40	extern "C" void begevtgengetx_(int ,int ,int ,int ,
	41	int ,int ,int ,int ,
	42	double ,double ,double *,
	43	double ,double ,double *,
	44	double ,double ,double *);
	45
	46	extern "C" void heplst_(int *);
	47
	48	extern "C" void photos_(int *);
	49
362af318	50	extern "C" void phoini_(int , int );
da0e9ce3	51
	52
	53	EvtPHOTOS::EvtPHOTOS(std::string photontype){
	54
	55	_photontype=photontype;
	56
	57	}
	58
	59	void EvtPHOTOS::doRadCorr( EvtParticle *p){
	60
	61	static int first=1;
	62
	63	//added by Lange Jan4,2000
	64	//allow to set photon tupe
	65	static EvtId GAMM=EvtPDL::getId(_photontype);
	66
	67	if (GAMM==EvtId(-1,-1)) {
	68	report(ERROR,"EvtGen") << "In EvtPHOTOS::doRadCorr():Particle:"<<
	69	_photontype<<" is not in EvtPDL"<<std::endl;
	70	::abort();
	71	}
	72
	73	if (first) {
	74	first=0;
362af318	75	int iseed1=(int)gRandom->Integer(31327);
	76	int iseed2=(int)gRandom->Integer(30080);
	77	phoini_(&iseed1,&iseed2);
da0e9ce3	78	}
	79
	80	double mpho=EvtPDL::getMeanMass(GAMM);
	81
	82	int entry,eventnum,numparticle,istat,partnum,mother;
	83	int daugfirst,dauglast;
	84
	85	int numparticlephotos;
	86
	87	double px,py,pz,e,m,x,y,z,t;
	88
	89	static EvtId dq=EvtPDL::getId("d");
	90	static EvtId adq=EvtPDL::getId("anti-d");
	91	static EvtId uq=EvtPDL::getId("u");
	92	static EvtId auq=EvtPDL::getId("anti-u");
	93	static EvtId sq=EvtPDL::getId("s");
	94	static EvtId asq=EvtPDL::getId("anti-s");
	95	static EvtId cq=EvtPDL::getId("c");
	96	static EvtId acq=EvtPDL::getId("anti-c");
	97	static EvtId bq=EvtPDL::getId("b");
	98	static EvtId abq=EvtPDL::getId("anti-b");
	99	static EvtId tq=EvtPDL::getId("t");
	100	static EvtId atq=EvtPDL::getId("anti-t");
	101	static EvtId vpho=EvtPDL::getId("vpho");
	102
	103	static EvtIdSet quarks(dq,adq,uq,auq,sq,asq,cq,acq,bq,abq,tq,atq);
	104
	105	if ( p->getId() == vpho ) return;
	106	if ( p->getNDaug() > 10 ) return;
	107
	108	px=0.0;
	109	py=0.0;
	110	pz=0.0;
	111	e=p->mass();
	112	m=p->mass();
	113	x=0.0;
	114	y=0.0;
	115	z=0.0;
	116	t=0.0;
	117
	118	entry=1;
	119	eventnum=1;
	120	numparticle=1;
	121	istat=2;
	122	partnum=EvtPDL::getStdHep(p->getId());
	123	mother=0;
	124	daugfirst=2;
	125	dauglast=1+p->getNDaug();
	126
	127	begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
	128	&mother,&daugfirst,&dauglast,
	129	&px,&py,&pz,&e,&m,&x,&y,&z,&t);
	130
	131	// std::cout << EvtPDL::name(p->getId()) <<" " ;
	132	for(size_t i=0;i<p->getNDaug();i++){
	133
	134	//No quarks to photos
	135	if (quarks.contains(p->getDaug(i)->getId())==1) continue;
	136
	137	px=p->getDaug(i)->getP4().get(1);
	138	py=p->getDaug(i)->getP4().get(2);
	139	pz=p->getDaug(i)->getP4().get(3);
	140	e=p->getDaug(i)->getP4().get(0);
	141	m=p->getDaug(i)->mass();
142	x=0.0;
143	y=0.0;
144	z=0.0;
145	t=0.0;
146
147	// std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " ;
148	entry+=1;
149	eventnum=1;
150	numparticle+=1;
151	istat=1;
152	partnum=EvtPDL::getStdHep(p->getDaug(i)->getId());
153	mother=1;
154	daugfirst=0;
155	dauglast=0;
156
157	begevtgenstorex_(&entry,&eventnum,&numparticle,&istat,&partnum,
158	&mother,&daugfirst,&dauglast,
159	&px,&py,&pz,&e,&m,&x,&y,&z,&t);
160
161
162	}
163	// std::cout << std::endl;
164	//can't use heplst since the common block used by the BaBar
165	//implementation of PHOTOS is renamed due to real4 vs real8
166	//problems.
167
168	//int mlst=1;
169
170	//heplst_(&mlst);
171
172	entry=1;
173
174	// report(INFO,"EvtGen") << "Doing photos " << EvtPDL::name(p->getId()) << endl;
175	photos_(&entry);
176	// report(INFO,"EvtGen") << "done\n";
177	begevtgengetx_(&entry,&eventnum,&numparticlephotos,&istat,&partnum,
178	&mother,&daugfirst,&dauglast,
179	&px,&py,&pz,&e,&m,&x,&y,&z,&t);
180
181
182	//report(INFO,"EvtGen") << "numparticlephotos:"<<numparticlephotos<<endl;
183
184	if (numparticle==numparticlephotos) return;
185
186	EvtVector4R new4mom;
187
188	int np;
189
190	for(size_t i=0;i<p->getNDaug();i++){
191
192	entry=i+2;
193
194	begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
195	&mother,&daugfirst,&dauglast,
196	&px,&py,&pz,&e,&m,&x,&y,&z,&t);
197
198	//this is needed to ensure that photos does not
199	//change the masses. But it will violate energy conservation!
200	double mp=p->getDaug(i)->mass();
201	e=sqrt(mpmp+pxpx+pypy+pzpz);
202
203	new4mom.set(e,px,py,pz);
204
205	p->getDaug(i)->setP4WithFSR(new4mom);
206
207
208
209	}
210
211	for(entry=numparticle+1;entry<=numparticlephotos;entry++){
212
213	begevtgengetx_(&entry,&eventnum,&np,&istat,&partnum,
214	&mother,&daugfirst,&dauglast,
215	&px,&py,&pz,&e,&m,&x,&y,&z,&t);
216
217	//Hack here to give the photon the mass of the generated particle
218	e=sqrt(mphompho+pxpx+pypy+pzpz);
219	new4mom.set(e,px,py,pz);
220
221	//new4mom.dump();
222
223	EvtPhotonParticle* gamma;
224	gamma=new EvtPhotonParticle;
225	gamma->init(GAMM,new4mom);
226	gamma->setFSRP4toZero();
227	// report(INFO,"EvtGen") << gamma << " " << p << " "<< px << " " << py << " " << pz << " " << p->getNDaug() << " " << EvtPDL::name(p->getId())<<" " << entry << " " <<numparticlephotos<<endl;
228	gamma->addDaug(p);
229
230	// p->getDaug(i)->set_type(EvtSpinType::PHOTON);
231
232	}
233	return ;
234	}
235