[u/mrichter/AliRoot.git] / TEvtGen / Tauola / TauolaParticle.cxx

#include "TauolaParticle.h"
#include "Log.h"
#include <stdexcept> 

namespace Tauolapp
{

double TauolaParticle::getPolarimetricX(){
  return m_pol_x;
}

double TauolaParticle::getPolarimetricY(){
  return m_pol_y;
}

double TauolaParticle::getPolarimetricZ(){
  return m_pol_z;
}


TauolaParticle * TauolaParticle::clone(){
  
  return createNewParticle(getPdgID(),getStatus(),getMass(),
                           getPx(),getPy(),getPz(),getE());

}

// NOTE: Not executed by release examples
double TauolaParticle::getAngle(TauolaParticle * other_particle){

  //use the dot product
  double x1 = getPx();
  double y1 = getPy();
  double z1 = getPz();
  double x2 = other_particle->getPx();
  double y2 = other_particle->getPy();
  double z2 = other_particle->getPz();

  return acos( (x1*x2+y1*y2+z1*z2) / sqrt((x1*x1+y1*y1+z1*z1)*(x2*x2+y2*y2+z2*z2)) );

}

// NOTE: Not executed by release examples
void TauolaParticle::add(TauolaParticle * other_particle){

  setPx(getPx() + other_particle->getPx());
  setPy(getPy() + other_particle->getPy());
  setPz(getPz() + other_particle->getPz());
  setE(getE() + other_particle->getE());
  setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() ));
}

void TauolaParticle::subtract(TauolaParticle * other_particle){

  setPx(getPx() - other_particle->getPx());
  setPy(getPy() - other_particle->getPy());
  setPz(getPz() - other_particle->getPz());
  setE(getE() - other_particle->getE());
  setMass( sqrt( getE()*getE()-getPx()*getPx()-getPy()*getPy()-getPz()*getPz() ));
}

int TauolaParticle::getSign(){
  if(getPdgID()== Tauola::getDecayingParticle())
    return SAME_SIGN;
  else if(getPdgID()== -1 * Tauola::getDecayingParticle())
    return OPPOSITE_SIGN;
  else
    return NA_SIGN;
}

bool TauolaParticle::hasDaughters(){
  if(getDaughters().size()==0)
    return 0;
  else
    return 1;
}

TauolaParticle * TauolaParticle::findLastSelf(){
  vector<TauolaParticle*> daughters = getDaughters();
  vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
  
  //get all daughters and look for stable with same pgd id
  for(;pcl_itr != daughters.end();pcl_itr++){
    if((*pcl_itr)->getPdgID()==this->getPdgID())
      return (*pcl_itr)->findLastSelf();
  }

  return this;
}

std::vector<TauolaParticle*> TauolaParticle::findProductionMothers(){
  vector<TauolaParticle*> mothers = getMothers();
  vector<TauolaParticle*>::iterator pcl_itr = mothers.begin();
  
    //get all mothers and check none have pdg id of this one
    for(;pcl_itr != mothers.end();pcl_itr++){
      if((*pcl_itr)->getPdgID()==this->getPdgID())
        return (*pcl_itr)->findProductionMothers();
    }
    return mothers;
}

void TauolaParticle::decay(){

  //Do the decay and set the polarimetric vectors
  TauolaDecay(getSign(),&m_pol_x, &m_pol_y, &m_pol_z, &m_pol_n);
}

void TauolaParticle::addDecayToEventRecord(){

  //Add to decay list used by f_filhep.c
  DecayList::addToEnd(this);
  TauolaWriteDecayToEventRecord(getSign());

  double xmom[4]={0};
  double *pp=xmom;

  if (Tauola::ion[2])
    for(int i=1;1;i++)
    {
      TauolaParticle *x;
      try{ x=DecayList::getParticle(i); }
      catch(std::out_of_range d) {break;}
      if(x->getPdgID()==221){ 

        // Fix 28.04.2011 The pp vector must have boost for eta undone
        TauolaParticle *x_copy = x->clone();
        if(getP(TauolaParticle::Z_AXIS)>0)
          x_copy->boostAlongZ(-getP(),getE());
        else
          x_copy->boostAlongZ(getP(),getE());

        pp[3]=x_copy->getE();
        pp[0]=x_copy->getPx();
        pp[1]=x_copy->getPy();
        pp[2]=x_copy->getPz();
        taueta_(pp,&i);
      } 
    }

  if (Tauola::ion[1])
    for(int i=1;1;i++)
    {
      TauolaParticle *x;
      try{ x=DecayList::getParticle(i); }
      catch(std::out_of_range d) {break;}
      if(x->getPdgID()==310){

        // Fix 28.04.2011 The pp vector must have boost for k0 undone
        TauolaParticle *x_copy = x->clone();
        if(getP(TauolaParticle::Z_AXIS)>0)
          x_copy->boostAlongZ(-getP(),getE());
        else
          x_copy->boostAlongZ(getP(),getE());

        pp[3]=x_copy->getE();
        pp[0]=x_copy->getPx();
        pp[1]=x_copy->getPy();
        pp[2]=x_copy->getPz();
        tauk0s_(pp,&i);
      } 
    }

  if (Tauola::ion[0])
    for(int i=1;1;i++)
    {
      TauolaParticle *x;
      try{ x=DecayList::getParticle(i); }
      catch(std::out_of_range d) {break;}
      if(x->getPdgID()==111){ 

        // Fix 28.04.2011 The pp vector must have boost for pi0 undone
        TauolaParticle *x_copy = x->clone();
        if(getP(TauolaParticle::Z_AXIS)>0)
          x_copy->boostAlongZ(-getP(),getE());
        else
          x_copy->boostAlongZ(getP(),getE());

        pp[3]=x_copy->getE();
        pp[0]=x_copy->getPx();
        pp[1]=x_copy->getPy();
        pp[2]=x_copy->getPz();
        taupi0_(pp,&i);
      } 
    }
  DecayList::clear();

  if(!hasDaughters())
    Log::Fatal("TAUOLA failed. No decay was created",5);
  //  checkMomentumConservation();
  //  decayEndgame(); // vertex shift was wrongly calculated, 
  //                     used 4-momenta should be in lab frame,
  //                     thanks to Sho Iwamoto for debug info.

}


void TauolaParticle::boostDaughtersFromRestFrame(TauolaParticle * tau_momentum){

  if(!hasDaughters()) //if there are no daughters
    return;

  vector<TauolaParticle*> daughters = getDaughters();
  vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
  
  //get all daughters then rotate and boost them.
  for(;pcl_itr != daughters.end();pcl_itr++){
 
    (*pcl_itr)->boostFromRestFrame(tau_momentum);
    (*pcl_itr)->boostDaughtersFromRestFrame(tau_momentum);
  }
  //checkMomentumConservation();
}

void TauolaParticle::boostDaughtersToRestFrame(TauolaParticle * tau_momentum){

  if(!hasDaughters()) //if there are no daughters
    return;
  // NOTE: Not executed by release examples
  //       because !hasDaughters() is always true
  vector<TauolaParticle*> daughters = getDaughters();
  vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
  
  //get all daughters then rotate and boost them.
  for(;pcl_itr != daughters.end();pcl_itr++){
 
    (*pcl_itr)->boostToRestFrame(tau_momentum);
    (*pcl_itr)->boostDaughtersToRestFrame(tau_momentum);
  }
  //checkMomentumConservation();
}


void TauolaParticle::boostToRestFrame(TauolaParticle * tau_momentum){

  double theta = tau_momentum->getRotationAngle(Y_AXIS);
  tau_momentum->rotate(Y_AXIS,theta);
  double phi = tau_momentum->getRotationAngle(X_AXIS);
  tau_momentum->rotate(Y_AXIS,-theta);

  //Now rotate coordinates to get boost in Z direction.
  rotate(Y_AXIS,theta);
  rotate(X_AXIS,phi);
  boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE());
  rotate(X_AXIS,-phi);
  rotate(Y_AXIS,-theta);

}

void TauolaParticle::boostFromRestFrame(TauolaParticle * tau_momentum){
  //get the rotation angles
  //and boost z

  double theta = tau_momentum->getRotationAngle(Y_AXIS);
  tau_momentum->rotate(Y_AXIS,theta);
  double phi = tau_momentum->getRotationAngle(X_AXIS);
  tau_momentum->rotate(Y_AXIS,-theta);

  //Now rotate coordinates to get boost in Z direction.
  rotate(Y_AXIS,theta);
  rotate(X_AXIS,phi);
  boostAlongZ(tau_momentum->getP(),tau_momentum->getE());
  rotate(X_AXIS,-phi);
  rotate(Y_AXIS,-theta);
}

/** Get the angle needed to rotate the 4 momentum vector so that
    the x (y) component disapears. (and the Z component is > 0) */
double TauolaParticle::getRotationAngle(int axis, int second_axis){

  /**if(getP(axis)==0){
    if(getPz()>0)
      return 0; //no rotaion required
    else
      return M_PI;
      }**/
  if(getP(second_axis)==0){
    if(getP(axis)>0)
      return -M_PI/2.0;
    else
      return M_PI/2.0;
  }
  if(getP(second_axis)>0)
    return -atan(getP(axis)/getP(second_axis));
  else
    return M_PI-atan(getP(axis)/getP(second_axis));

}

/** Boost this vector along the Z direction.
    Assume no momentum components in the X or Y directions. */
void TauolaParticle::boostAlongZ(double boost_pz, double boost_e){

  // Boost along the Z axis
  double m=sqrt(boost_e*boost_e-boost_pz*boost_pz);

  double p=getPz();
  double e=getE();

  setPz((boost_e*p + boost_pz*e)/m);
  setE((boost_pz*p + boost_e*e )/m);
}

/** Rotation around an axis X or Y */
void TauolaParticle::rotate(int axis,double theta, int second_axis){
  
  double temp_px=getP(axis);
  double temp_pz=getP(second_axis);
  setP(axis,cos(theta)*temp_px + sin(theta)*temp_pz);
  setP(second_axis,-sin(theta)*temp_px + cos(theta)*temp_pz);
}

void TauolaParticle::rotateDaughters(int axis,double theta, int second_axis){
  if(!hasDaughters()) //if there are no daughters
    return;

  vector<TauolaParticle*> daughters = getDaughters();
  vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
  
  //get all daughters then rotate and boost them.
  for(;pcl_itr != daughters.end();pcl_itr++){
 
    (*pcl_itr)->rotate(axis,theta,second_axis);
    (*pcl_itr)->rotateDaughters(axis,theta,second_axis);
  }
  //checkMomentumConservation();
}

double TauolaParticle::getMass(){
  double e_sq=getE()*getE();
  double p_sq=getP()*getP();

  if(e_sq>p_sq)
    return sqrt(e_sq-p_sq);
  else
    return -1*sqrt(p_sq-e_sq); //if it's negative
}

double TauolaParticle::getP(){
  return sqrt(getPx()*getPx()+getPy()*getPy()+getPz()*getPz());
}

double TauolaParticle::getP(int axis){
  if(axis==X_AXIS)
    return getPx();

  if(axis==Y_AXIS)
    return getPy();

  if(axis==Z_AXIS)
    return getPz();

  return 0;
}

void TauolaParticle::setP(int axis, double p_component){
  if(axis==X_AXIS)
    setPx(p_component);
  if(axis==Y_AXIS)
    setPy(p_component);
  if(axis==Z_AXIS)
    setPz(p_component);
}

} // namespace Tauolapp
Commit	Line	Data
0ca57c2f	1	#include "TauolaParticle.h"
	2	#include "Log.h"
	3	#include <stdexcept>
	4
	5	namespace Tauolapp
	6	{
	7
	8	double TauolaParticle::getPolarimetricX(){
	9	return m_pol_x;
	10	}
	11
	12	double TauolaParticle::getPolarimetricY(){
	13	return m_pol_y;
	14	}
	15
	16	double TauolaParticle::getPolarimetricZ(){
	17	return m_pol_z;
	18	}
	19
	20
	21	TauolaParticle * TauolaParticle::clone(){
	22
	23	return createNewParticle(getPdgID(),getStatus(),getMass(),
	24	getPx(),getPy(),getPz(),getE());
	25
	26	}
	27
	28	// NOTE: Not executed by release examples
	29	double TauolaParticle::getAngle(TauolaParticle * other_particle){
	30
	31	//use the dot product
	32	double x1 = getPx();
	33	double y1 = getPy();
	34	double z1 = getPz();
	35	double x2 = other_particle->getPx();
	36	double y2 = other_particle->getPy();
	37	double z2 = other_particle->getPz();
	38
	39	return acos( (x1x2+y1y2+z1z2) / sqrt((x1x1+y1y1+z1z1)(x2x2+y2y2+z2z2)) );
	40
	41	}
	42
	43	// NOTE: Not executed by release examples
	44	void TauolaParticle::add(TauolaParticle * other_particle){
	45
	46	setPx(getPx() + other_particle->getPx());
	47	setPy(getPy() + other_particle->getPy());
	48	setPz(getPz() + other_particle->getPz());
	49	setE(getE() + other_particle->getE());
	50	setMass( sqrt( getE()getE()-getPx()getPx()-getPy()getPy()-getPz()getPz() ));
	51	}
	52
	53	void TauolaParticle::subtract(TauolaParticle * other_particle){
	54
	55	setPx(getPx() - other_particle->getPx());
	56	setPy(getPy() - other_particle->getPy());
	57	setPz(getPz() - other_particle->getPz());
	58	setE(getE() - other_particle->getE());
	59	setMass( sqrt( getE()getE()-getPx()getPx()-getPy()getPy()-getPz()getPz() ));
	60	}
	61
	62	int TauolaParticle::getSign(){
	63	if(getPdgID()== Tauola::getDecayingParticle())
	64	return SAME_SIGN;
65	else if(getPdgID()== -1 * Tauola::getDecayingParticle())
66	return OPPOSITE_SIGN;
67	else
68	return NA_SIGN;
69	}
70
71	bool TauolaParticle::hasDaughters(){
72	if(getDaughters().size()==0)
73	return 0;
74	else
75	return 1;
76	}
77
78	TauolaParticle * TauolaParticle::findLastSelf(){
79	vector<TauolaParticle*> daughters = getDaughters();
80	vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
81
82	//get all daughters and look for stable with same pgd id
83	for(;pcl_itr != daughters.end();pcl_itr++){
84	if((*pcl_itr)->getPdgID()==this->getPdgID())
85	return (*pcl_itr)->findLastSelf();
86	}
87
88	return this;
89	}
90
91	std::vector<TauolaParticle*> TauolaParticle::findProductionMothers(){
92	vector<TauolaParticle*> mothers = getMothers();
93	vector<TauolaParticle*>::iterator pcl_itr = mothers.begin();
94
95	//get all mothers and check none have pdg id of this one
96	for(;pcl_itr != mothers.end();pcl_itr++){
97	if((*pcl_itr)->getPdgID()==this->getPdgID())
98	return (*pcl_itr)->findProductionMothers();
99	}
100	return mothers;
101	}
102
103	void TauolaParticle::decay(){
104
105	//Do the decay and set the polarimetric vectors
106	TauolaDecay(getSign(),&m_pol_x, &m_pol_y, &m_pol_z, &m_pol_n);
107	}
108
109	void TauolaParticle::addDecayToEventRecord(){
110
111	//Add to decay list used by f_filhep.c
112	DecayList::addToEnd(this);
113	TauolaWriteDecayToEventRecord(getSign());
114
115	double xmom[4]={0};
116	double *pp=xmom;
117
118	if (Tauola::ion[2])
119	for(int i=1;1;i++)
120	{
121	TauolaParticle *x;
122	try{ x=DecayList::getParticle(i); }
123	catch(std::out_of_range d) {break;}
124	if(x->getPdgID()==221){
125
126	// Fix 28.04.2011 The pp vector must have boost for eta undone
127	TauolaParticle *x_copy = x->clone();
128	if(getP(TauolaParticle::Z_AXIS)>0)
129	x_copy->boostAlongZ(-getP(),getE());
130	else
131	x_copy->boostAlongZ(getP(),getE());
132
133	pp[3]=x_copy->getE();
134	pp[0]=x_copy->getPx();
135	pp[1]=x_copy->getPy();
136	pp[2]=x_copy->getPz();
137	taueta_(pp,&i);
138	}
139	}
140
141	if (Tauola::ion[1])
142	for(int i=1;1;i++)
143	{
144	TauolaParticle *x;
145	try{ x=DecayList::getParticle(i); }
146	catch(std::out_of_range d) {break;}
147	if(x->getPdgID()==310){
148
149	// Fix 28.04.2011 The pp vector must have boost for k0 undone
150	TauolaParticle *x_copy = x->clone();
151	if(getP(TauolaParticle::Z_AXIS)>0)
152	x_copy->boostAlongZ(-getP(),getE());
153	else
154	x_copy->boostAlongZ(getP(),getE());
155
156	pp[3]=x_copy->getE();
157	pp[0]=x_copy->getPx();
158	pp[1]=x_copy->getPy();
159	pp[2]=x_copy->getPz();
160	tauk0s_(pp,&i);
161	}
162	}
163
164	if (Tauola::ion[0])
165	for(int i=1;1;i++)
166	{
167	TauolaParticle *x;
168	try{ x=DecayList::getParticle(i); }
169	catch(std::out_of_range d) {break;}
170	if(x->getPdgID()==111){
171
172	// Fix 28.04.2011 The pp vector must have boost for pi0 undone
173	TauolaParticle *x_copy = x->clone();
174	if(getP(TauolaParticle::Z_AXIS)>0)
175	x_copy->boostAlongZ(-getP(),getE());
176	else
177	x_copy->boostAlongZ(getP(),getE());
178
179	pp[3]=x_copy->getE();
180	pp[0]=x_copy->getPx();
181	pp[1]=x_copy->getPy();
182	pp[2]=x_copy->getPz();
183	taupi0_(pp,&i);
184	}
185	}
186	DecayList::clear();
187
188	if(!hasDaughters())
189	Log::Fatal("TAUOLA failed. No decay was created",5);
190	// checkMomentumConservation();
191	// decayEndgame(); // vertex shift was wrongly calculated,
192	// used 4-momenta should be in lab frame,
193	// thanks to Sho Iwamoto for debug info.
194
195	}
196
197
198	void TauolaParticle::boostDaughtersFromRestFrame(TauolaParticle * tau_momentum){
199
200	if(!hasDaughters()) //if there are no daughters
201	return;
202
203	vector<TauolaParticle*> daughters = getDaughters();
204	vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
205
206	//get all daughters then rotate and boost them.
207	for(;pcl_itr != daughters.end();pcl_itr++){
208
209	(*pcl_itr)->boostFromRestFrame(tau_momentum);
210	(*pcl_itr)->boostDaughtersFromRestFrame(tau_momentum);
211	}
212	//checkMomentumConservation();
213	}
214
215	void TauolaParticle::boostDaughtersToRestFrame(TauolaParticle * tau_momentum){
216
217	if(!hasDaughters()) //if there are no daughters
218	return;
219	// NOTE: Not executed by release examples
220	// because !hasDaughters() is always true
221	vector<TauolaParticle*> daughters = getDaughters();
222	vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
223
224	//get all daughters then rotate and boost them.
225	for(;pcl_itr != daughters.end();pcl_itr++){
226
227	(*pcl_itr)->boostToRestFrame(tau_momentum);
228	(*pcl_itr)->boostDaughtersToRestFrame(tau_momentum);
229	}
230	//checkMomentumConservation();
231	}
232
233
234	void TauolaParticle::boostToRestFrame(TauolaParticle * tau_momentum){
235
236	double theta = tau_momentum->getRotationAngle(Y_AXIS);
237	tau_momentum->rotate(Y_AXIS,theta);
238	double phi = tau_momentum->getRotationAngle(X_AXIS);
239	tau_momentum->rotate(Y_AXIS,-theta);
240
241	//Now rotate coordinates to get boost in Z direction.
242	rotate(Y_AXIS,theta);
243	rotate(X_AXIS,phi);
244	boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE());
245	rotate(X_AXIS,-phi);
246	rotate(Y_AXIS,-theta);
247
248	}
249
250	void TauolaParticle::boostFromRestFrame(TauolaParticle * tau_momentum){
251	//get the rotation angles
252	//and boost z
253
254	double theta = tau_momentum->getRotationAngle(Y_AXIS);
255	tau_momentum->rotate(Y_AXIS,theta);
256	double phi = tau_momentum->getRotationAngle(X_AXIS);
257	tau_momentum->rotate(Y_AXIS,-theta);
258
259	//Now rotate coordinates to get boost in Z direction.
260	rotate(Y_AXIS,theta);
261	rotate(X_AXIS,phi);
262	boostAlongZ(tau_momentum->getP(),tau_momentum->getE());
263	rotate(X_AXIS,-phi);
264	rotate(Y_AXIS,-theta);
265	}
266
267	/** Get the angle needed to rotate the 4 momentum vector so that
268	the x (y) component disapears. (and the Z component is > 0) */
269	double TauolaParticle::getRotationAngle(int axis, int second_axis){
270
271	/**if(getP(axis)==0){
272	if(getPz()>0)
273	return 0; //no rotaion required
274	else
275	return M_PI;
276	}**/
277	if(getP(second_axis)==0){
278	if(getP(axis)>0)
279	return -M_PI/2.0;
280	else
281	return M_PI/2.0;
282	}
283	if(getP(second_axis)>0)
284	return -atan(getP(axis)/getP(second_axis));
285	else
286	return M_PI-atan(getP(axis)/getP(second_axis));
287
288	}
289
290	/** Boost this vector along the Z direction.
291	Assume no momentum components in the X or Y directions. */
292	void TauolaParticle::boostAlongZ(double boost_pz, double boost_e){
293
294	// Boost along the Z axis
295	double m=sqrt(boost_eboost_e-boost_pzboost_pz);
296
297	double p=getPz();
298	double e=getE();
299
300	setPz((boost_ep + boost_pze)/m);
301	setE((boost_pzp + boost_ee )/m);
302	}
303
304	/** Rotation around an axis X or Y */
305	void TauolaParticle::rotate(int axis,double theta, int second_axis){
306
307	double temp_px=getP(axis);
308	double temp_pz=getP(second_axis);
309	setP(axis,cos(theta)temp_px + sin(theta)temp_pz);
310	setP(second_axis,-sin(theta)temp_px + cos(theta)temp_pz);
311	}
312
313	void TauolaParticle::rotateDaughters(int axis,double theta, int second_axis){
314	if(!hasDaughters()) //if there are no daughters
315	return;
316
317	vector<TauolaParticle*> daughters = getDaughters();
318	vector<TauolaParticle*>::iterator pcl_itr = daughters.begin();
319
320	//get all daughters then rotate and boost them.
321	for(;pcl_itr != daughters.end();pcl_itr++){
322
323	(*pcl_itr)->rotate(axis,theta,second_axis);
324	(*pcl_itr)->rotateDaughters(axis,theta,second_axis);
325	}
326	//checkMomentumConservation();
327	}
328
329	double TauolaParticle::getMass(){
330	double e_sq=getE()*getE();
331	double p_sq=getP()*getP();
332
333	if(e_sq>p_sq)
334	return sqrt(e_sq-p_sq);
335	else
336	return -1*sqrt(p_sq-e_sq); //if it's negative
337	}
338
339	double TauolaParticle::getP(){
340	return sqrt(getPx()getPx()+getPy()getPy()+getPz()*getPz());
341	}
342
343	double TauolaParticle::getP(int axis){
344	if(axis==X_AXIS)
345	return getPx();
346
347	if(axis==Y_AXIS)
348	return getPy();
349
350	if(axis==Z_AXIS)
351	return getPz();
352
353	return 0;
354	}
355
356	void TauolaParticle::setP(int axis, double p_component){
357	if(axis==X_AXIS)
358	setPx(p_component);
359	if(axis==Y_AXIS)
360	setPy(p_component);
361	if(axis==Z_AXIS)
362	setPz(p_component);
363	}
364
365	} // namespace Tauolapp