Line data Source code
1 : #include "Tauola.h"
2 : #include "TauolaParticlePair.h"
3 : #include "TauolaLog.h"
4 : #include <stdlib.h>
5 : #include <math.h>
6 : #include <iostream>
7 :
8 : namespace Tauolapp
9 : {
10 :
11 : /** constructor. Get the mothers, grandmothers and siblings of the tau */
12 0 : TauolaParticlePair::TauolaParticlePair(std::vector<TauolaParticle*> &particle_list){
13 0 : TauolaParticle *particle=particle_list.back();
14 0 : particle_list.pop_back();
15 0 : setBornKinematics(0,0,-1.,0.); //set the default born variables
16 :
17 : // In case of decayOne() we need only the tau information
18 0 : if(Tauola::isUsingDecayOne())
19 : {
20 0 : m_mother = 0;
21 0 : m_mother_exists = false;
22 0 : m_production_particles.push_back(particle);
23 0 : m_final_particles.push_back(particle);
24 0 : initializeDensityMatrix();
25 0 : Log::AddDecay(0);
26 0 : return;
27 : }
28 :
29 : //See if there are any mothers
30 0 : std::vector<TauolaParticle *> temp_mothers = particle->findProductionMothers();
31 :
32 : // Case that there are no mothers or grandmothers
33 0 : if(temp_mothers.size()==0){
34 : // NOTE: Not executed by release examples
35 : // However, such cases were present if tests used older Pythia8.1 or so.
36 0 : Log::Warning()<< "WARNING: Could not find taus mother or grandmothers. "
37 0 : << "Ignoring spin effects" << std::endl;
38 0 : m_mother = 0;
39 0 : m_mother_exists = false;
40 0 : m_production_particles.push_back(particle);
41 0 : m_final_particles.push_back(particle->findLastSelf());
42 0 : initializeDensityMatrix();
43 0 : Log::AddDecay(1);
44 0 : return;
45 : }
46 :
47 : //fill the sibling pointers
48 0 : std::vector<TauolaParticle*> temp_daughters;
49 0 : temp_daughters = temp_mothers.at(0)->getDaughters();
50 0 : if(temp_daughters.size()==0)
51 0 : Log::Fatal("WARNING: Something wrong with event structure or there is a bug in the TAUOLA interface.",6);
52 :
53 0 : m_production_particles=temp_daughters;
54 0 : m_final_particles.push_back(particle->findLastSelf());
55 :
56 : //Second tau-like particle selection should match properties of the hard (exotic) process. At present, solution
57 : //match one of the possible diagrams contributing to SM process. Rare case like tau+ tau+ nutau nutau will not be treted
58 : //accordingly to any diagram of SM
59 0 : for(signed int i=0; i < (int) m_production_particles.size(); i++)
60 : {
61 : //check if it has opposite PDGID sign
62 0 : if(m_final_particles.at(0)->getPdgID()*m_production_particles.at(i)->getPdgID()>=0) continue;
63 :
64 0 : if(m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_PLUS||
65 0 : m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_MINUS)
66 : {
67 : //tau+ or tau- - check if it's on the particle_list
68 : int j=-1;
69 0 : for(j=0;j<(int)particle_list.size();j++)
70 0 : if(m_production_particles.at(i)->getBarcode()==particle_list.at(j)->getBarcode()) break;
71 0 : if(j>=(int)particle_list.size()) continue;
72 :
73 : //exists on the list - add to m_final_particles and delete from particle_list
74 0 : m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
75 0 : particle_list.erase(particle_list.begin()+j);
76 0 : }
77 0 : else if(m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_NEUTRINO||
78 0 : m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_ANTINEUTRINO)
79 : {
80 : //neutrino - for now - just add to m_final_particles
81 0 : m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
82 : }
83 : }
84 :
85 :
86 : //fill the mother and grandmother pointers
87 0 : if(temp_mothers.size()==1){ //one mother
88 0 : m_mother_exists = true;
89 0 : m_mother = temp_mothers.at(0);
90 0 : m_grandmothers = m_mother->findProductionMothers();
91 0 : Log::AddDecay(3);
92 : }
93 : else{ //no mother, but grandparents exist
94 0 : m_mother_exists = false;
95 0 : m_grandmothers = temp_mothers;
96 0 : m_mother = makeTemporaryMother(m_production_particles);
97 0 : Log::AddDecay(2);
98 : }
99 :
100 0 : initializeDensityMatrix();
101 :
102 : return;
103 0 : }
104 :
105 :
106 : /** The axis is defined by the boosting routine but our standard convention
107 : is:
108 : - Axis 3 is along the direction of the +ve tau,
109 : - Axis 1 is perpendicular to the reaction plane (sign=??)
110 : - Axis 2 is defined through the vector product so the system
111 : is right handed (?? check). Axis 1,2 and 3 are parrellel for the 1st
112 : and second tau. */
113 : void TauolaParticlePair::initializeDensityMatrix(){
114 0 : int incoming_pdg_id=0;
115 0 : int outgoing_pdg_id=0;
116 0 : double invariant_mass_squared=-5.0;
117 0 : double cosTheta=3.0;
118 :
119 : //initialize all elements of the density matrix to zero
120 0 : for(int x = 0; x < 4; x ++) {
121 0 : for(int y = 0; y < 4; y ++)
122 0 : m_R[x][y] = 0;
123 : }
124 :
125 0 : m_R[0][0]=1;
126 :
127 0 : if(Tauola::isUsingDecayOne())
128 : {
129 0 : const double *pol = Tauola::getDecayOnePolarization();
130 :
131 0 : m_R[0][1]=pol[0];
132 0 : m_R[0][2]=pol[1];
133 0 : m_R[0][3]=pol[2];
134 :
135 0 : m_R[1][0]=pol[0];
136 0 : m_R[2][0]=pol[1];
137 0 : m_R[3][0]=pol[2];
138 0 : }
139 :
140 0 : if(!m_mother) return;
141 : // fill the matrix depending on mother
142 :
143 :
144 : // do scalar-pseudoscalar mixed higgs case separately since
145 : // switch doesn't allow non-constants in a case statement.
146 : // very annoying!
147 :
148 0 : if(m_mother->getPdgID()==Tauola::getHiggsScalarPseudoscalarPDG()){
149 0 : if(!Tauola::spin_correlation.HIGGS_H) return;
150 :
151 0 : double phi = Tauola::getHiggsScalarPseudoscalarMixingAngle();
152 0 : double mass_ratio = Tauola::getTauMass()/m_mother->getMass();
153 :
154 0 : double beta = sqrt(1-4*mass_ratio*mass_ratio);
155 0 : double denominator = pow(cos(phi)*beta,2)+pow(sin(phi),2);
156 :
157 0 : m_R[0][0]= 1;
158 0 : m_R[1][1]= (pow(cos(phi)*beta,2)-pow(sin(phi),2))/denominator;
159 0 : m_R[1][2]= 2*cos(phi)*sin(phi)*beta/denominator;
160 0 : m_R[2][1]= -m_R[1][2];
161 0 : m_R[2][2]= m_R[1][1];
162 0 : m_R[3][3]= -1;
163 0 : }
164 : else {
165 :
166 : double pz = 0.0;
167 :
168 0 : switch(m_mother->getPdgID()){
169 :
170 : case TauolaParticle::Z0:
171 0 : if(!Tauola::spin_correlation.Z0) break;
172 : // Here we calculate SVAR and COSTHE as well as IDE and IDF
173 : // ANGULU(&IDE,&IDF,&SVAR,&COSTHE);
174 : // this is ++ configuration of Fig 5 from HEPPH0101311:
175 : //pol_tau_minus[2]=-1; pol_tau_plus[2]=1;
176 :
177 0 : pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
178 0 : m_R[0][0]=1;
179 0 : m_R[0][3]=2*pz-1;
180 0 : m_R[3][0]=2*pz-1;
181 0 : m_R[3][3]=1;
182 : // alternatively this may be overwritten if better solution exist
183 0 : recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
184 0 : break;
185 :
186 : case TauolaParticle::GAMMA:
187 0 : if(!Tauola::spin_correlation.GAMMA) break;
188 : // Here we calculate SVAR and COSTHE as well as IDE and IDF
189 : // ANGULU(&IDE,&IDF,&SVAR,&COSTHE);
190 : // this is ++ configuration of Fig 5 from HEPPH0101311:
191 : //pol_tau_minus[2]=-1; pol_tau_plus[2]=1;
192 :
193 0 : pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
194 0 : m_R[0][0]=1;
195 0 : m_R[0][3]=2*pz-1;
196 0 : m_R[3][0]=2*pz-1;
197 0 : m_R[3][3]=1;
198 : // alternatively this may be overwritten if better solution exist
199 0 : recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
200 0 : break;
201 :
202 : //scalar higgs
203 : case TauolaParticle::HIGGS:
204 0 : if(!Tauola::spin_correlation.HIGGS) break;
205 0 : m_R[0][0]=1;
206 0 : m_R[1][1]=1;
207 0 : m_R[2][2]=1;
208 0 : m_R[3][3]=-1;
209 0 : break;
210 :
211 : //pseudoscalar higgs case
212 : case TauolaParticle::HIGGS_A:
213 0 : if(!Tauola::spin_correlation.HIGGS_A) break;
214 0 : m_R[0][0]=1;
215 0 : m_R[1][1]=-1;
216 0 : m_R[2][2]=-1;
217 0 : m_R[3][3]=-1;
218 0 : break;
219 :
220 : case TauolaParticle::HIGGS_PLUS:
221 0 : if(!Tauola::spin_correlation.HIGGS_PLUS) break;
222 0 : m_R[0][0]=1;
223 0 : m_R[0][3]=-1;
224 0 : m_R[3][0]=-1;
225 0 : break;
226 :
227 : case TauolaParticle::HIGGS_MINUS:
228 0 : if(!Tauola::spin_correlation.HIGGS_MINUS) break;
229 0 : m_R[0][0]=1;
230 0 : m_R[0][3]=-1;
231 0 : m_R[3][0]=-1;
232 0 : break;
233 :
234 : case TauolaParticle::W_PLUS:
235 0 : if(!Tauola::spin_correlation.W_PLUS) break;
236 0 : m_R[0][0]=1;
237 0 : m_R[0][3]=1; //tau minus (tau minus is on the -ve Z axis)
238 0 : m_R[3][0]=1; //tau plus
239 0 : break;
240 :
241 : case TauolaParticle::W_MINUS:
242 0 : if(!Tauola::spin_correlation.W_MINUS) break;
243 0 : m_R[0][0]=1;
244 0 : m_R[0][3]=1; //tau minus (tau minus is on the -ve Z axis)
245 0 : m_R[3][0]=1; //tau plus
246 0 : break;
247 :
248 : //ignore spin effects when mother is unknown
249 : default:
250 0 : m_R[0][0]=1;
251 0 : break;
252 : }
253 : }
254 :
255 0 : }
256 :
257 : /**************************************************************/
258 : void TauolaParticlePair::setBornKinematics(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared,double cosTheta){
259 0 : Tauola::buf_incoming_pdg_id=incoming_pdg_id;
260 0 : Tauola::buf_outgoing_pdg_id=outgoing_pdg_id;
261 0 : Tauola::buf_invariant_mass_squared=invariant_mass_squared;
262 0 : Tauola::buf_cosTheta=cosTheta;
263 : //cout<<"(TauolaParticlePair::Just to be sure:) "<<buf_incoming_pdg_id<<" "<<buf_outgoing_pdg_id<<" "<<buf_invariant_mass_squared<<" "<<buf_cosTheta<<endl;
264 : // m_R[0][0] to be added in next step;
265 0 : }
266 :
267 : /**************************************************************/
268 : double TauolaParticlePair::getZPolarization(int *incoming_pdg_id, int *outgoing_pdg_id, double *invariant_mass_squared,double *cosTheta){
269 :
270 : //defaults
271 0 : *incoming_pdg_id = TauolaParticle::ELECTRON;
272 0 : *cosTheta = 0 ;
273 0 : *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
274 0 : *invariant_mass_squared = pow(m_mother->getMass(),2);
275 0 : setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
276 :
277 : //TRIVIAL CASE:
278 : //if we don't know the incoming beams then
279 : //return the average z polarisation
280 0 : if(m_grandmothers.size()<2){
281 0 : Log::Warning()<<"Not enough mothers of Z to "
282 0 : <<"calculate cos(theta) between "
283 0 : <<"incoming about outgoing beam"
284 0 : << endl;
285 0 : return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
286 : invariant_mass_squared, cosTheta);
287 : }
288 :
289 0 : if(!getTauPlus(m_production_particles)||!getTauMinus(m_production_particles)){
290 0 : Log::Error()<<"tau+ or tau- not found in Z decay"<< endl;
291 0 : return 0;
292 : }
293 :
294 : //NOW CHECK FOR THE DIFFICULT EVENTS:
295 : //case f1 + f2 + f3 -> Z -> tau tau
296 : //case f1 + f2 -> Z + f3, Z-> tau tau or f1 + f2 -> tau tau f3
297 : //case f1 + f2 -> Z -> tau tau gamma
298 0 : if(m_grandmothers.size()>2 ||
299 0 : (m_grandmothers.at(0)->getDaughters().size()>1 && m_mother_exists==true) ||
300 0 : m_production_particles.size() > 2){
301 :
302 : //make a vector of the extra grandmother particles
303 0 : vector<TauolaParticle*> extra_grandmothers;
304 0 : for(int i=0; i<(int) m_grandmothers.size(); i++){
305 : // temp_grandmothers.push_back(m_grandmothers.at(i));
306 0 : if(m_grandmothers.at(i)!=getGrandmotherPlus(m_grandmothers)&&
307 0 : m_grandmothers.at(i)!=getGrandmotherMinus(m_grandmothers))
308 0 : extra_grandmothers.push_back(m_grandmothers.at(i));
309 : }
310 :
311 : //make a vector of the tau siblings
312 : //and copy all the production particle vector.
313 0 : vector<TauolaParticle*> extra_tau_siblings;
314 0 : vector<TauolaParticle*> temp_production_particles;
315 0 : for(int i=0; i<(int) m_production_particles.size(); i++){
316 0 : if(m_production_particles.at(i)!=getTauPlus(m_production_particles)&&
317 0 : m_production_particles.at(i)!=getTauMinus(m_production_particles))
318 0 : extra_tau_siblings.push_back(m_production_particles.at(i));
319 : }
320 :
321 : //make a vector of the Z's sibling
322 0 : vector<TauolaParticle*> extra_Z_siblings;
323 0 : for(int i=0; m_mother_exists && i<(int) m_grandmothers.at(0)->getDaughters().size(); i++){
324 0 : if(m_grandmothers.at(0)->getDaughters().at(i)->getPdgID()!=TauolaParticle::Z0)
325 0 : extra_Z_siblings.push_back(m_grandmothers.at(0)->getDaughters().at(i));
326 : }
327 :
328 : //make temporary particles for the effect beams
329 : //and copy into a vector
330 0 : std::vector<TauolaParticle*> effective_taus;
331 0 : effective_taus.push_back(getTauPlus(m_production_particles)->clone());
332 0 : effective_taus.push_back(getTauMinus(m_production_particles)->clone());
333 :
334 : //copy grandmothers into the m_grandmothers vector since we want this to be perminante.
335 0 : TauolaParticle * g1 = getGrandmotherPlus(m_grandmothers)->clone();
336 0 : TauolaParticle * g2 = getGrandmotherMinus(m_grandmothers)->clone();
337 0 : m_grandmothers.clear();
338 0 : m_grandmothers.push_back(g1);
339 0 : m_grandmothers.push_back(g2);
340 :
341 : //loop over extra grandmothers
342 0 : for(int i=0; i<(int) extra_grandmothers.size(); i++)
343 0 : addToBeam(extra_grandmothers.at(i),&m_grandmothers,0);
344 :
345 : //loop over siblings to the Z
346 0 : for(int i=0; i<(int) extra_Z_siblings.size(); i++)
347 0 : addToBeam(extra_Z_siblings.at(i),0,&m_grandmothers);
348 :
349 : //loop over siblings of the taus
350 0 : for(int i=0; i<(int) extra_tau_siblings.size() ; i++){
351 0 : if(m_mother_exists)
352 0 : addToBeam(extra_tau_siblings.at(i),&effective_taus,0);
353 : else
354 0 : addToBeam(extra_tau_siblings.at(i),&effective_taus,&m_grandmothers);
355 : }
356 : //And we are finish making the effective income and outgoing beams
357 :
358 0 : TauolaParticle * temp_mother = makeTemporaryMother(effective_taus);
359 0 : *invariant_mass_squared = pow(temp_mother->getMass(),2);
360 :
361 : //now we are ready to do the boosting, calculate the angle
362 : //between incoming and outgoing, and get the polarisation of the Z
363 :
364 0 : double angle1,angle2,angle3;
365 0 : boostFromLabToTauPairFrame(&angle1, &angle2, &angle3, temp_mother,
366 0 : m_grandmothers, effective_taus);
367 :
368 : /*double theta1 = -getGrandmotherPlus(m_grandmothers)->getRotationAngle(TauolaParticle::Y_AXIS);
369 : double theta2 = -(M_PI+getGrandmotherMinus(m_grandmothers)->getRotationAngle(TauolaParticle::Y_AXIS));
370 : *cosTheta = cos((theta1+theta2)/2.0); //just average the angles for now.*/
371 :
372 0 : TauolaParticle *tM=getTauPlus(effective_taus);
373 0 : TauolaParticle *gM=getGrandmotherPlus(m_grandmothers);
374 :
375 0 : double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
376 0 : sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
377 0 : sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
378 0 : tM=getTauMinus(effective_taus);
379 0 : gM=getGrandmotherMinus(m_grandmothers);
380 :
381 0 : double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
382 0 : sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
383 0 : sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
384 0 : double sintheta1 = sqrt(1-costheta1*costheta1);
385 0 : double sintheta2 = sqrt(1-costheta2*costheta2);
386 0 : double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
387 :
388 0 : *cosTheta = avg;
389 :
390 0 : *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
391 :
392 0 : if(abs(*incoming_pdg_id)==21 || abs(*incoming_pdg_id)==22)
393 : {
394 0 : *incoming_pdg_id = -getGrandmotherMinus(m_grandmothers)->getPdgID();
395 : //cout<<"INFO:\tgluon pdg id changed!"<<endl;
396 0 : }
397 :
398 0 : if( abs(*incoming_pdg_id)> 8 &&
399 0 : abs(*incoming_pdg_id)!=11 &&
400 0 : abs(*incoming_pdg_id)!=13 &&
401 0 : abs(*incoming_pdg_id)!=15 )
402 : {
403 0 : Log::Error()<<"Second class disaster: incoming_pdg_id = "<<*incoming_pdg_id<<endl;
404 :
405 : // Return avarage Z polarization
406 0 : *incoming_pdg_id = TauolaParticle::ELECTRON;
407 0 : *cosTheta = 0 ;
408 0 : *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
409 :
410 0 : return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
411 : invariant_mass_squared, cosTheta);
412 : }
413 :
414 0 : boostFromTauPairToLabFrame(angle1, angle2, angle3, temp_mother,
415 0 : m_grandmothers, effective_taus);
416 0 : setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
417 0 : return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
418 :
419 0 : }
420 : //REGULAR CASE:
421 : //f1 + f2 -> Z -> tau+ tau - or f1 f2 -> tau+ tau-
422 : //This includes Z -> tau tau, tau -> gamma tau?
423 :
424 0 : double angle1,angle2,angle3;
425 0 : boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
426 0 : m_mother,m_grandmothers,m_production_particles);
427 :
428 0 : TauolaParticle *tM=getTauPlus(m_production_particles);
429 0 : TauolaParticle *gM=getGrandmotherPlus(m_grandmothers);
430 0 : double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
431 0 : sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
432 0 : sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
433 :
434 0 : tM=getTauMinus(m_production_particles);
435 0 : gM=getGrandmotherMinus(m_grandmothers);
436 :
437 0 : double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
438 0 : sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
439 0 : sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
440 :
441 0 : double sintheta1 = sqrt(1-costheta1*costheta1);
442 0 : double sintheta2 = sqrt(1-costheta2*costheta2);
443 :
444 0 : double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
445 :
446 0 : *cosTheta = avg;
447 :
448 0 : *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
449 :
450 0 : boostFromTauPairToLabFrame(angle1, angle2, angle3,
451 0 : m_mother,m_grandmothers,m_production_particles);
452 :
453 0 : setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
454 0 : return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
455 : // return 0.5 - (-0.12 + 0.12*cosTheta)/2;
456 0 : }
457 :
458 : /** WHERE WE CALCULATE THE EFFECTIVE BEAMS **/
459 : /** This is where we decide which particle should be added into which beam,
460 : add it and change the flavour if necessary. candidates_same
461 : are on the same side of the vertex as the particle. This is needed
462 : for negative the particle 4-momentum. **/
463 : void TauolaParticlePair::addToBeam(TauolaParticle * pcle,
464 : std::vector<TauolaParticle*> *candidates_same,
465 : std::vector<TauolaParticle*> *candidates_opp){
466 :
467 :
468 : double s=0.,o=-10.;
469 : TauolaParticle * s_beam = NULL;
470 : TauolaParticle * o_beam = NULL;
471 :
472 0 : if(candidates_same){
473 0 : double s0 =1.0/getVirtuality(pcle,candidates_same->at(0),false);
474 0 : double s1 = 1.0/getVirtuality(pcle,candidates_same->at(1),false);
475 0 : if(s0>s1){
476 : s=s0;
477 0 : s_beam = candidates_same->at(0);
478 0 : }
479 : else{
480 : s=s1;
481 0 : s_beam = candidates_same->at(1);
482 : }
483 0 : }
484 0 : if(candidates_opp){
485 :
486 0 : double o0 =1.0/getVirtuality(pcle,candidates_opp->at(0),true);
487 0 : double o1 = 1.0/getVirtuality(pcle,candidates_opp->at(1),true);
488 0 : if(o0>o1){
489 : o=o0;
490 0 : o_beam = candidates_opp->at(0);
491 0 : }
492 : else{
493 : o=o1;
494 0 : o_beam = candidates_opp->at(1);
495 : }
496 0 : }
497 :
498 0 : if(s>o)
499 : {
500 0 : s_beam->add(pcle);
501 0 : int pdg1 = pcle->getPdgID();
502 0 : int pdg2 = s_beam->getPdgID();
503 0 : if(abs(pdg2)==15) return;
504 0 : if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) s_beam->setPdgID( pdg1);
505 0 : }
506 : else
507 : {
508 0 : o_beam->subtract(pcle);
509 0 : int pdg1 = pcle->getPdgID();
510 0 : int pdg2 = o_beam->getPdgID();
511 0 : if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) o_beam->setPdgID(-pdg1);
512 : }
513 :
514 : //should we also do something with the flavours??
515 :
516 0 : }
517 :
518 : double TauolaParticlePair::getVirtuality(TauolaParticle * p1, TauolaParticle*p2, bool flip){
519 :
520 : //if one particle in a gluon and the other a lepton
521 0 : if((p1->getPdgID()==TauolaParticle::GLUON&&abs(p2->getPdgID())>10&&abs(p2->getPdgID())<20) ||
522 0 : (p2->getPdgID()==TauolaParticle::GLUON&&abs(p1->getPdgID())>10&&abs(p1->getPdgID())<20))
523 0 : return -1;
524 :
525 : //add some others...
526 :
527 : //otherwise we calculate the virtuality:
528 0 : double kp = p1->getE()*p2->getE() - p1->getPx()*p2->getPx()
529 0 : - p1->getPy()*p2->getPy() - p1->getPz()*p2->getPz();
530 0 : if(flip) // Note that we have 1/(p1+p2)^2 or 1/(p1-p2)^2 and p2^2=0
531 0 : kp = kp - p1->getMass()*p1->getMass();
532 :
533 : double q = 1;
534 0 : if(p1->getPdgID()==TauolaParticle::GAMMA){
535 0 : if(abs(p2->getPdgID())==TauolaParticle::UP ||
536 0 : abs(p2->getPdgID())==TauolaParticle::CHARM ||
537 0 : abs(p2->getPdgID())==TauolaParticle::TOP)
538 0 : q=2.0/3.0;
539 0 : if(abs(p2->getPdgID())==TauolaParticle::DOWN ||
540 0 : abs(p2->getPdgID())==TauolaParticle::STRANGE ||
541 0 : abs(p2->getPdgID())==TauolaParticle::BOTTOM)
542 0 : q=1.0/3.0;
543 : }
544 0 : if(p2->getPdgID()==TauolaParticle::GAMMA){
545 0 : if(abs(p1->getPdgID())==TauolaParticle::UP ||
546 0 : abs(p1->getPdgID())==TauolaParticle::CHARM ||
547 0 : abs(p1->getPdgID())==TauolaParticle::TOP)
548 0 : q=2.0/3.0;
549 0 : if(abs(p1->getPdgID())==TauolaParticle::DOWN ||
550 0 : abs(p1->getPdgID())==TauolaParticle::STRANGE ||
551 0 : abs(p1->getPdgID())==TauolaParticle::BOTTOM)
552 0 : q=1.0/3.0;
553 : }
554 :
555 0 : return fabs(2*kp)/q;
556 0 : }
557 :
558 : /***********************************************************
559 : ** TauolaParticlePair::decayTauPair().
560 : ** Generate tau decay
561 : ************************************************************/
562 : void TauolaParticlePair::decayTauPair(){
563 :
564 : //initalize h vectors in case the partner is not a tau
565 0 : double h_tau_minus[4]={2,0,0,0}; //2 used to compensate for maximum weight
566 0 : double h_tau_plus[4]={2,0,0,0}; //in the case when there is only 1 tau
567 :
568 0 : TauolaParticle * tau_minus = getTauMinus(m_final_particles);
569 0 : TauolaParticle * tau_plus = getTauPlus(m_final_particles);
570 :
571 : //now calculate the spin weight
572 0 : for(double weight=0; weight < Tauola::randomDouble();){
573 : //call tauola decay and get polarimetric vectors
574 0 : if(tau_minus){
575 0 : Tauola::redefineTauMinusProperties(tau_minus);
576 0 : tau_minus->decay();
577 0 : h_tau_minus[0]=1;
578 0 : h_tau_minus[1]=tau_minus->getPolarimetricX();
579 0 : h_tau_minus[2]=tau_minus->getPolarimetricY();
580 0 : h_tau_minus[3]=tau_minus->getPolarimetricZ();
581 0 : }
582 0 : if(tau_plus){
583 0 : Tauola::redefineTauPlusProperties(tau_plus);
584 0 : tau_plus->decay();
585 0 : h_tau_plus[0]=1;
586 0 : h_tau_plus[1]=tau_plus->getPolarimetricX();
587 0 : h_tau_plus[2]=tau_plus->getPolarimetricY();
588 0 : h_tau_plus[3]=tau_plus->getPolarimetricZ();
589 0 : }
590 : // cout<<" tau? tau? "<< tau_plus->getPdgID()<<" "<< tau_minus->getPdgID()<<endl;
591 : // cout<<" przy uzyciu rrrRRRrrrRRR "<< m_R[0][0]<<" " <<m_R[3][3]<<" " << m_R[0][3]<<" " <<m_R[3][0] <<endl;
592 : //calculate the weight
593 : weight=0;
594 0 : for(int i =0; i<4; i++){
595 0 : for(int j=0; j<4; j++)
596 0 : weight+=m_R[i][j]*h_tau_minus[i]*h_tau_plus[j];
597 : }
598 0 : weight = weight/4.0;
599 : }
600 : double wthel[4]={0};
601 0 : wthel[0]=(h_tau_minus[0]+h_tau_minus[3])*(h_tau_plus[0]+h_tau_plus[3])*(m_R[0][0]+m_R[0][3]+m_R[3][0]+m_R[3][3]);
602 0 : wthel[1]=(h_tau_minus[0]+h_tau_minus[3])*(h_tau_plus[0]-h_tau_plus[3])*(m_R[0][0]-m_R[0][3]+m_R[3][0]-m_R[3][3]);
603 0 : wthel[2]=(h_tau_minus[0]-h_tau_minus[3])*(h_tau_plus[0]+h_tau_plus[3])*(m_R[0][0]+m_R[0][3]-m_R[3][0]-m_R[3][3]);
604 0 : wthel[3]=(h_tau_minus[0]-h_tau_minus[3])*(h_tau_plus[0]-h_tau_plus[3])*(m_R[0][0]-m_R[0][3]-m_R[3][0]+m_R[3][3]);
605 :
606 0 : double rn=Tauola::randomDouble();
607 0 : if (rn>(wthel[0]+wthel[1]+wthel[2])/(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,-1);
608 0 : else if (rn>(wthel[0]+wthel[1]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,+1);
609 0 : else if (rn>(wthel[0]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities( 1,-1);
610 0 : else Tauola::setHelicities( 1, 1);
611 :
612 :
613 :
614 :
615 : //boost into the frame used to define the density matrices
616 : //and add the tau's new daughters.
617 :
618 0 : double angle1,angle2,angle3;
619 0 : TauolaParticle * mum = makeTemporaryMother(m_final_particles);
620 :
621 0 : if(!Tauola::isUsingDecayOneBoost())
622 0 : boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
623 0 : mum,m_grandmothers,m_final_particles);
624 :
625 : //add the accepted decay into the event record
626 0 : if(tau_plus)
627 0 : tau_plus->addDecayToEventRecord();
628 0 : if(tau_minus)
629 0 : tau_minus->addDecayToEventRecord();
630 :
631 0 : if(!Tauola::isUsingDecayOneBoost())
632 0 : boostFromTauPairToLabFrame(angle1,angle2,angle3,
633 0 : mum,m_grandmothers,m_final_particles);
634 :
635 : // Apply final decay modification, once taus are in lab frame.
636 : // At present 28.02.11, vertex position shift (in HepMC) is implemented.
637 : // Thanks Sho Iwamoto for feedback
638 0 : if(tau_plus)
639 0 : tau_plus->decayEndgame();
640 0 : if(tau_minus)
641 0 : tau_minus->decayEndgame();
642 :
643 0 : }
644 :
645 : /***********************************************************
646 : ** Below are the methods used for boosting and rotating
647 : ** into another reference frame.
648 : ************************************************************/
649 :
650 : /** Step 1. (Transformation A). Any modification to this method also requires a modification
651 : to the inverse method boostFromTauPairFrameToLab (transformation A^-1).*/
652 : void TauolaParticlePair::boostFromLabToTauPairFrame(double * rotation_angle1,
653 : double * rotation_angle2,
654 : double * rotation_angle3,
655 : TauolaParticle * mother,
656 : vector<TauolaParticle *> grandmothers,
657 : vector<TauolaParticle *> taus
658 : ){ //in positive z axis (+1) //or negative z axis (-1)
659 :
660 : /** boost all gradmothers and daughters (taus, neutrinos, etc,)
661 : to the mothers rest frame */
662 :
663 0 : for(int i=0; i< (int) grandmothers.size(); i++)
664 0 : grandmothers.at(i)->boostToRestFrame(mother);
665 : //If taus.size()==1 then taus[1] is the same as mum. Disaster to be avoided.
666 0 : if(taus.size()!=1)
667 0 : for(int i=0; i< (int) taus.size(); i++){
668 0 : taus.at(i)->boostToRestFrame(mother);
669 : // NOTE: Following line has no effect in release examples
670 : // because taus don't have daughters at this step
671 0 : taus.at(i)->boostDaughtersToRestFrame(mother);
672 0 : }
673 :
674 : /** rotate all particles so taus are on the z axis */
675 0 : if(getTauPlus(taus)){ //if there's a tau+ use this to get the rotation angle
676 0 : *rotation_angle1 = getTauPlus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
677 0 : rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
678 0 : *rotation_angle2 = getTauPlus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
679 0 : rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
680 0 : }
681 : else{ //otherwise use the tau-
682 0 : *rotation_angle1 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
683 0 : rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
684 0 : *rotation_angle2 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
685 0 : rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
686 : }
687 :
688 : //now rotate incoming beams so they are aligned in the y-z plane
689 0 : if(grandmothers.size()<1){ //if there are no grandmothers
690 0 : *rotation_angle3 = 0;
691 0 : return;
692 : }
693 :
694 : //the first grandmother will have a positive y component
695 0 : if(getGrandmotherPlus(grandmothers))
696 0 : *rotation_angle3 = getGrandmotherPlus(grandmothers)->getRotationAngle(TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
697 : else
698 0 : *rotation_angle3 = M_PI+getGrandmotherMinus(grandmothers)->getRotationAngle(TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
699 :
700 0 : rotateSystem(grandmothers,taus,*rotation_angle3,TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
701 :
702 0 : }
703 :
704 : /** Reverses boostFromLabtoMotherFrame. The three rotation angle must be provided.
705 : Any modification to this would require a modification to boostFromLabToTauPairFrame
706 : since this is the inverse transformation (Step 2: A^-1). */
707 : void TauolaParticlePair::boostFromTauPairToLabFrame(double rotation_angle1,
708 : double rotation_angle2,
709 : double rotation_angle3,
710 : TauolaParticle * mother,
711 : vector<TauolaParticle *> grandmothers,
712 : vector<TauolaParticle *> taus){
713 :
714 :
715 0 : if(mother==0) //if there are no mothers
716 : return;
717 :
718 :
719 : //rotate grand mothers back away from the y-z plane
720 0 : rotateSystem(grandmothers,taus,-rotation_angle3, TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
721 :
722 : //rotate all so that taus are no longer on the z axis
723 0 : rotateSystem(grandmothers,taus,-rotation_angle2,TauolaParticle::X_AXIS);
724 0 : rotateSystem(grandmothers,taus,-rotation_angle1,TauolaParticle::Y_AXIS);
725 :
726 :
727 : /*** boost grandmothers and daughters (taus) back into the lab frame */
728 0 : for(int i=0; i< (int) grandmothers.size(); i++)
729 0 : grandmothers.at(i)->boostFromRestFrame(mother);
730 : //If taus.size()==1 then taus[1] is the same as mum. Disaster to be avoided.
731 0 : if(taus.size()!=1)
732 0 : for(int i=0; i< (int) taus.size(); i++){
733 0 : taus.at(i)->boostFromRestFrame(mother);
734 0 : taus.at(i)->boostDaughtersFromRestFrame(mother);
735 0 : }
736 0 : }
737 :
738 : void TauolaParticlePair::rotateSystem(vector<TauolaParticle *> grandmothers,
739 : vector<TauolaParticle *> taus,
740 : double theta, int axis, int axis2){
741 0 : for(int i=0; i< (int) grandmothers.size(); i++)
742 0 : grandmothers.at(i)->rotate(axis,theta,axis2);
743 0 : for(int i=0; i< (int) taus.size(); i++){
744 0 : taus.at(i)->rotate(axis,theta,axis2);
745 0 : taus.at(i)->rotateDaughters(axis,theta,axis2);
746 : }
747 0 : }
748 :
749 :
750 :
751 : /***********************************************************
752 : Some extra useful methods
753 : ************************************************************/
754 : TauolaParticle * TauolaParticlePair::makeTemporaryMother(vector<TauolaParticle *> taus){
755 :
756 : //calculate mothers 4-momentum based on the daughters.
757 : double e=0;
758 : double px=0;
759 : double py=0;
760 : double pz=0;
761 :
762 : bool tau_plus = false;
763 : bool tau_minus = false;
764 : bool tau_neutrino = false;
765 : bool tau_antineutrino = false;
766 :
767 0 : for(int d=0; d < (int) taus.size(); d++){
768 0 : TauolaParticle * daughter = taus.at(d);
769 0 : e+=daughter->getE();
770 0 : px+=daughter->getPx();
771 0 : py+=daughter->getPy();
772 0 : pz+=daughter->getPz();
773 0 : switch(daughter->getPdgID()){
774 : case TauolaParticle::TAU_PLUS:
775 : tau_plus=true;
776 0 : break;
777 : case TauolaParticle::TAU_MINUS:
778 : tau_minus=true;
779 0 : break;
780 : case TauolaParticle::TAU_NEUTRINO:
781 : tau_neutrino=true;
782 0 : break;
783 : case TauolaParticle::TAU_ANTINEUTRINO:
784 : tau_antineutrino=true;
785 0 : break;
786 : }
787 : }
788 0 : double m = e*e-px*px-py*py-pz*pz;
789 0 : if(m<0)
790 0 : m= -sqrt(-m);
791 : else
792 0 : m=sqrt(m);
793 :
794 : //look for mothers type:
795 : int pdg=0;
796 :
797 : //Assume Z
798 0 : if(tau_plus&&tau_minus)
799 0 : pdg=TauolaParticle::Z0 ;
800 :
801 : //Assume W+
802 0 : if(tau_plus&&tau_neutrino)
803 0 : pdg=TauolaParticle::W_PLUS;
804 :
805 : //Assume W-
806 0 : if(tau_minus&&tau_antineutrino)
807 0 : pdg=TauolaParticle::W_MINUS;
808 :
809 : // now make the mother
810 0 : return taus.at(0)->createNewParticle(pdg,2,m,px,py,pz,e);
811 : }
812 :
813 : // see if "particle" is one of the final taus in this tau pair
814 : // (based on particle barcode)
815 : // NOTE: Not executed by release examples
816 : bool TauolaParticlePair::contains(TauolaParticle * particle){
817 :
818 0 : for(int i=0; i< (int) m_final_particles.size(); i++){
819 0 : if(m_final_particles.at(i)->getBarcode()==particle->getBarcode())
820 0 : return true;
821 : }
822 0 : return false;
823 0 : }
824 :
825 : TauolaParticle * TauolaParticlePair::getTauMinus(vector<TauolaParticle*> taus){
826 0 : for(int i=0; i< (int) taus.size(); i++){
827 0 : if(taus.at(i)->getPdgID()==TauolaParticle::TAU_MINUS)
828 0 : return taus.at(i);
829 : }
830 0 : return 0;
831 0 : }
832 :
833 : TauolaParticle * TauolaParticlePair::getTauPlus(vector<TauolaParticle*> taus){
834 0 : for(int i=0; i< (int) taus.size(); i++){
835 0 : if(taus.at(i)->getPdgID()==TauolaParticle::TAU_PLUS)
836 0 : return taus.at(i);
837 : }
838 0 : return 0;
839 0 : }
840 :
841 : TauolaParticle * TauolaParticlePair::getGrandmotherPlus(vector<TauolaParticle*> grandmothers){
842 : //choose based on energy if there are more than one??
843 : double e = -1e30;
844 : int index = -1;
845 0 : for(int i=0; i< (int) grandmothers.size(); i++){
846 0 : if( (grandmothers.at(i)->getPdgID()<0 && grandmothers.at(i)->getPdgID()>-9) ||
847 0 : grandmothers.at(i)->getPdgID()== TauolaParticle::POSITRON||
848 0 : grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_PLUS){
849 0 : if(e<grandmothers.at(i)->getE()){
850 0 : e=grandmothers.at(i)->getE();
851 : index=i;
852 0 : }
853 : }
854 : }
855 0 : if(index==-1)
856 : {
857 0 : for(int i=0; i< (int) grandmothers.size(); i++)
858 : {
859 0 : if(grandmothers.at(i)->getPdgID()==21 || grandmothers.at(i)->getPdgID()==22)
860 : {
861 : index=i;
862 0 : e=grandmothers.at(i)->getE();
863 0 : break;
864 : }
865 : }
866 0 : }
867 0 : if(index==-1) index=0;
868 0 : if(e==0)
869 : {
870 0 : grandmothers.at(index)->print();
871 0 : return 0;
872 : }
873 : else
874 0 : return grandmothers.at(index);
875 0 : }
876 :
877 : TauolaParticle * TauolaParticlePair::getGrandmotherMinus(vector<TauolaParticle*> grandmothers){
878 : //choose based on energy if there are more than one??
879 : double e = -1e30;
880 : int index = -1;
881 0 : for(int i=0; i< (int) grandmothers.size(); i++){
882 0 : if( (grandmothers.at(i)->getPdgID()>0 && grandmothers.at(i)->getPdgID()<9) ||
883 0 : grandmothers.at(i)->getPdgID()== TauolaParticle::ELECTRON||
884 0 : grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_MINUS){
885 0 : if(e<grandmothers.at(i)->getE()){
886 0 : e=grandmothers.at(i)->getE();
887 : index=i;
888 0 : }
889 : }
890 : }
891 0 : if(index==-1)
892 : {
893 0 : for(int i=(int) grandmothers.size()-1; i>=0 ; i--)
894 : {
895 0 : if(grandmothers.at(i)->getPdgID()==21||grandmothers.at(i)->getPdgID()==22)
896 : {
897 : index=i;
898 0 : e=grandmothers.at(i)->getE();
899 0 : break;
900 : }
901 : }
902 0 : }
903 0 : if(index==-1) index=0;
904 0 : if(e==0)
905 0 : return 0;
906 : else
907 0 : return grandmothers.at(index);
908 0 : }
909 :
910 :
911 :
912 : void TauolaParticlePair::print(){
913 :
914 0 : Log::RedirectOutput(Log::Info());
915 0 : std::cout << "Daughters final:" << std::endl;
916 0 : for(int i=0; i< (int) m_final_particles.size(); i++)
917 0 : m_final_particles.at(i)->print();
918 :
919 :
920 0 : std::cout << "Daughters at production:" << std::endl;
921 0 : for(int i=0; i< (int) m_production_particles.size(); i++)
922 0 : m_production_particles.at(i)->print();
923 :
924 :
925 0 : std::cout << "Mother particle: " << std::endl;
926 0 : if(m_mother)
927 0 : m_mother->print();
928 :
929 0 : std::cout << "Grandmother particles: " << std::endl;
930 0 : for(int i=0; i< (int) m_grandmothers.size(); i++)
931 0 : m_grandmothers.at(i)->print();
932 :
933 0 : Log::RevertOutput();
934 0 : }
935 :
936 :
937 : void TauolaParticlePair::checkMomentumConservation(){
938 :
939 0 : for(int i=0; i< (int) m_final_particles.size(); i++)
940 0 : m_final_particles.at(i)->checkMomentumConservation();
941 :
942 0 : for(int i=0; i< (int) m_production_particles.size(); i++)
943 0 : m_production_particles.at(i)->checkMomentumConservation();
944 :
945 0 : if(m_mother)
946 0 : m_mother->checkMomentumConservation();
947 :
948 0 : for(int i=0; i< (int) m_grandmothers.size(); i++)
949 0 : m_grandmothers.at(i)->checkMomentumConservation();
950 :
951 0 : }
952 :
953 : void TauolaParticlePair::recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta){
954 :
955 0 : if (abs(outgoing_pdg_id)!=15)
956 : {
957 0 : Log::Warning()<<"interface was not used for taus pdg id="<<outgoing_pdg_id<<endl;
958 0 : return;
959 : }
960 :
961 : double (*T) [Tauola::NCOS][4][4] = NULL;
962 : double (*Tw) [Tauola::NCOS] = NULL;
963 : double (*Tw0)[Tauola::NCOS] = NULL;
964 : double smin = 0.0; // Tauola::sminA;
965 : double smax = 0.0; // Tauola::smaxA;
966 : double step = 0.0; // (smaxl-sminl)/(Tauola::NS1-1);
967 :
968 : // Select table for appropriate incoming particles
969 0 : switch(abs(incoming_pdg_id)){
970 : case 11:
971 0 : if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
972 : {
973 : T = Tauola::table11B;
974 : Tw = Tauola::wtable11B;
975 : Tw0 = Tauola::w0table11B;
976 : smin = Tauola::sminB;
977 : smax = Tauola::smaxB;
978 0 : step = (smax-smin)/(Tauola::NS2-1);
979 0 : }
980 0 : else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
981 : {
982 : T = Tauola::table11C;
983 : Tw = Tauola::wtable11C;
984 : Tw0 = Tauola::w0table11C;
985 : smin = Tauola::sminC;
986 : smax = Tauola::smaxC;
987 0 : step = (smax-smin)/(Tauola::NS3-1);
988 0 : }
989 : else
990 : {
991 : T = Tauola::table11A;
992 : Tw = Tauola::wtable11A;
993 : Tw0 = Tauola::w0table11A;
994 0 : smin = Tauola::sminA;
995 0 : smax = Tauola::smaxA;
996 0 : step = (smax-smin)/(Tauola::NS1-1);
997 : }
998 : break;
999 :
1000 : // NOTE: Use the same tables for 1, 3 and 5
1001 : case 1:
1002 : case 3:
1003 : case 5:
1004 0 : if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
1005 : {
1006 : T = Tauola::table1B;
1007 : Tw = Tauola::wtable1B;
1008 : Tw0 = Tauola::w0table1B;
1009 : smin = Tauola::sminB;
1010 : smax = Tauola::smaxB;
1011 0 : step = (smax-smin)/(Tauola::NS2-1);
1012 0 : }
1013 0 : else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
1014 : {
1015 : T = Tauola::table1C;
1016 : Tw = Tauola::wtable1C;
1017 : Tw0 = Tauola::w0table1C;
1018 : smin = Tauola::sminC;
1019 : smax = Tauola::smaxC;
1020 0 : step = (smax-smin)/(Tauola::NS3-1);
1021 0 : }
1022 : else
1023 : {
1024 : T = Tauola::table1A;
1025 : Tw = Tauola::wtable1A;
1026 : Tw0 = Tauola::w0table1A;
1027 0 : smin = Tauola::sminA;
1028 0 : smax = Tauola::smaxA;
1029 0 : step = (smax-smin)/(Tauola::NS1-1);
1030 : }
1031 : break;
1032 :
1033 : // NOTE: Use the same tables for 2 and 4
1034 : case 2:
1035 : case 4:
1036 0 : if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
1037 : {
1038 : T = Tauola::table2B;
1039 : Tw = Tauola::wtable2B;
1040 : Tw0 = Tauola::w0table2B;
1041 : smin = Tauola::sminB;
1042 : smax = Tauola::smaxB;
1043 0 : step = (smax-smin)/(Tauola::NS2-1);
1044 0 : }
1045 0 : else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
1046 : {
1047 : T = Tauola::table2C;
1048 : Tw = Tauola::wtable2C;
1049 : Tw0 = Tauola::w0table2C;
1050 : smin = Tauola::sminC;
1051 : smax = Tauola::smaxC;
1052 0 : step = (smax-smin)/(Tauola::NS3-1);
1053 0 : }
1054 : else
1055 : {
1056 : T = Tauola::table2A;
1057 : Tw = Tauola::wtable2A;
1058 : Tw0 = Tauola::w0table2A;
1059 0 : smin = Tauola::sminA;
1060 0 : smax = Tauola::smaxA;
1061 0 : step = (smax-smin)/(Tauola::NS1-1);
1062 : }
1063 : break;
1064 :
1065 : default:
1066 0 : Log::Warning()<<"interface was not used for proper beams pdg id="<<incoming_pdg_id<<endl;
1067 0 : return;
1068 : }
1069 :
1070 : // Check if we have found a table for matrix recalculation
1071 0 : if (T[0][0][0][0]<0.5) return;
1072 :
1073 : // If mass is too small
1074 0 : if (invariant_mass_squared <= exp(Tauola::sminA)){
1075 : double mta = 1.777;
1076 : double cosf = cosTheta;
1077 : double s = invariant_mass_squared;
1078 0 : double sinf = sqrt(1-cosf*cosf);
1079 0 : double MM = sqrt(4e0*mta*mta/s);
1080 0 : double xnorm = 1+cosf*cosf + MM*MM*sinf*sinf;
1081 :
1082 : // Zero matrix
1083 0 : for (int i=0;i<4;i++)
1084 0 : for (int j=0;j<4;j++)
1085 0 : m_R[i][j]=0.0;
1086 :
1087 0 : m_R[0][0] = (1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
1088 0 : m_R[1][1] = (-(1- MM*MM)*sinf*sinf)/xnorm;
1089 0 : m_R[2][2] = ( (1+ MM*MM)*sinf*sinf)/xnorm;
1090 0 : m_R[3][3] = (1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
1091 0 : m_R[2][3] = (2*MM*sinf*cosf)/xnorm;
1092 0 : m_R[3][2] = (2*MM*sinf*cosf)/xnorm;
1093 :
1094 : // Get weights
1095 : double w = 1.;
1096 : double w0 = 1.;
1097 :
1098 0 : if (Tauola::wtable2A[0][0]>0 ) w=Tauola::wtable2A[0][0];
1099 0 : else if(Tauola::wtable1A[0][0]>0 ) w=Tauola::wtable1A[0][0];
1100 0 : else if(Tauola::wtable11A[0][0]>0) w=Tauola::wtable11A[0][0];
1101 :
1102 0 : if (Tauola::wtable2A[0][0]>0 ) w0=Tauola::w0table2A[0][0];
1103 0 : else if(Tauola::wtable1A[0][0]>0 ) w0=Tauola::w0table1A[0][0];
1104 0 : else if(Tauola::wtable11A[0][0]>0) w0=Tauola::w0table11A[0][0];
1105 :
1106 : // Tauola::setEWwt(w/w0);
1107 0 : Tauola::setEWwt(w,w0);
1108 : return;
1109 : } // if(mass too small)
1110 :
1111 : int x = 0;
1112 : double buf = smin;
1113 : double remnant = 0.0;
1114 :
1115 : // Interpolate s
1116 0 : if(T==Tauola::table11A || T==Tauola::table1A || T== Tauola::table2A)
1117 : {
1118 0 : while(log(invariant_mass_squared) > buf){
1119 0 : x++;
1120 0 : buf += (step);
1121 : }
1122 0 : remnant = (log(invariant_mass_squared)-(buf-step))/step;
1123 0 : }
1124 : else
1125 : {
1126 0 : while(invariant_mass_squared > buf){
1127 0 : x++;
1128 0 : buf += step;
1129 : }
1130 0 : remnant = (invariant_mass_squared-(buf-step))/step;
1131 : }
1132 :
1133 : double cmin = -1.;
1134 : int y = 0;
1135 : double remnantc = 0.0;
1136 :
1137 : // Interpolate cosTheta
1138 : buf = cmin+1./Tauola::NCOS;
1139 0 : while(cosTheta > buf){
1140 0 : y++;
1141 0 : buf+=2./Tauola::NCOS;
1142 : }
1143 :
1144 0 : remnantc = (cosTheta-(buf-2./Tauola::NCOS))/(2./Tauola::NCOS);
1145 :
1146 : // Special cases at edges - extrapolation
1147 : bool isUsingExtrapolation = false;
1148 0 : if(y >= Tauola::NCOS) { isUsingExtrapolation = true; y = Tauola::NCOS-1; }
1149 0 : if(y == 0) { isUsingExtrapolation = true; y = 1; }
1150 :
1151 : // Bilinear interpolation
1152 : double b11,b21,b12,b22;
1153 :
1154 0 : for (int i=0; i<4; i++){
1155 0 : for (int j=0; j<4; j++){
1156 :
1157 0 : if(isUsingExtrapolation){
1158 0 : if(y == 1){
1159 0 : b11 = 2*T[x-1][0][i][j] - T[x-1][1][i][j];
1160 0 : b21 = 2*T[x][0][i][j] - T[x][1][i][j];
1161 : b12 = T[x-1][0][i][j];
1162 : b22 = T[x][0][i][j];
1163 0 : }
1164 : else{
1165 0 : b11 = T[x-1][y][i][j];
1166 0 : b21 = T[x][y][i][j];
1167 0 : b12 = 2*T[x-1][y][i][j] - T[x-1][y-1][i][j];
1168 0 : b22 = 2*T[x][y][i][j] - T[x][y-1][i][j];
1169 : }
1170 : } // if(isUsingExtrapolation)
1171 : else{
1172 0 : b11 = T[x-1][y-1][i][j];
1173 0 : b21 = T[x][y-1][i][j];
1174 0 : b12 = T[x-1][y][i][j];
1175 0 : b22 = T[x][y][i][j];
1176 : }
1177 0 : m_R[i][j] = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1178 0 : + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1179 : } // for(j)
1180 : } // for(i)
1181 :
1182 : // Calculate electroweak weights
1183 : double w,w0;
1184 :
1185 0 : if(isUsingExtrapolation){
1186 0 : if(y == 1){
1187 0 : b11 = 2*Tw[x-1][0] - Tw[x-1][1];
1188 0 : b21 = 2*Tw[x][0] - Tw[x][1];
1189 : b12 = Tw[x-1][0];
1190 : b22 = Tw[x][0];
1191 0 : }
1192 : else
1193 : {
1194 0 : b11 = Tw[x-1][y];
1195 0 : b21 = Tw[x][y];
1196 0 : b12 = 2*Tw[x-1][y] - Tw[x-1][y-1];
1197 0 : b22 = 2*Tw[x][y] - Tw[x][y-1];
1198 : }
1199 : } // if(isUsingExtrapolation)
1200 : else
1201 : {
1202 0 : b11 = Tw[x-1][y-1];
1203 0 : b21 = Tw[x][y-1];
1204 0 : b12 = Tw[x-1][y];
1205 0 : b22 = Tw[x][y];
1206 : }
1207 :
1208 0 : w = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1209 0 : + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1210 :
1211 0 : if(isUsingExtrapolation){
1212 0 : if(y == 1){
1213 0 : b11 = 2*Tw0[x-1][0] - Tw0[x-1][1];
1214 0 : b21 = 2*Tw0[x][0] - Tw0[x][1];
1215 : b12 = Tw0[x-1][0];
1216 : b22 = Tw0[x][0];
1217 0 : }
1218 : else{
1219 0 : b11 = Tw0[x-1][y];
1220 0 : b21 = Tw0[x][y];
1221 0 : b12 = 2*Tw0[x-1][y] - Tw0[x-1][y-1];
1222 0 : b22 = 2*Tw0[x][y] - Tw0[x][y-1];
1223 : }
1224 : } // if (isUsingExtrapolation)
1225 : else{
1226 0 : b11 = Tw0[x-1][y-1];
1227 0 : b21 = Tw0[x][y-1];
1228 0 : b12 = Tw0[x-1][y];
1229 0 : b22 = Tw0[x][y];
1230 : }
1231 :
1232 0 : w0 = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1233 0 : + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1234 :
1235 0 : Tauola::setEWwt(w,w0);
1236 0 : }
1237 :
1238 : } // namespace Tauolapp
|
