Line data Source code
1 : #include "EvtGenBase/EvtPatches.hh"
2 : #include "EvtGenBase/EvtReport.hh"
3 : #include "EvtGenBase/EvtSpinAmp.hh"
4 : #include <stdlib.h>
5 :
6 : using std::endl;
7 :
8 : std::ostream&
9 : operator<<( std::ostream& os, const EvtSpinAmp& amp )
10 : {
11 0 : vector<int> index = amp.iterinit();
12 :
13 0 : os << ":";
14 : do {
15 0 : os<<"<";
16 0 : for(size_t i=0; i<index.size()-1; ++i) {
17 0 : os<<index[i];
18 : }
19 0 : os<<index[index.size()-1]<<">"<<amp(index)<<":";
20 0 : } while( amp.iterate( index ) );
21 :
22 : return os;
23 0 : }
24 :
25 : EvtSpinAmp operator*( const EvtComplex& real, const EvtSpinAmp& cont )
26 : {
27 0 : EvtSpinAmp ret( cont );
28 :
29 0 : for( size_t i=0; i<ret._elem.size(); ++i ) {
30 0 : ret._elem[i] *= real;
31 : }
32 :
33 : return ret;
34 0 : }
35 :
36 : EvtSpinAmp operator*( const EvtSpinAmp& cont, const EvtComplex& real )
37 : {
38 0 : return real*cont;
39 : }
40 :
41 : EvtSpinAmp operator/( const EvtSpinAmp& cont, const EvtComplex& real )
42 : {
43 0 : EvtSpinAmp ret( cont );
44 :
45 0 : for( size_t i=0; i<ret._elem.size(); ++i ) {
46 0 : ret._elem[i] /= real;
47 : }
48 :
49 : return ret;
50 0 : }
51 :
52 : vector<unsigned int> EvtSpinAmp::calctwospin( const vector<EvtSpinType::spintype>& type ) const
53 : {
54 0 : vector<unsigned int> twospin;
55 :
56 0 : for( size_t i=0; i<type.size(); ++i ) {
57 0 : twospin.push_back( EvtSpinType::getSpin2( type[i] ) );
58 : }
59 :
60 : return twospin;
61 0 : }
62 :
63 0 : EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type )
64 0 : {
65 : int num = 1;
66 0 : _type = type;
67 0 : _twospin=calctwospin( type );
68 :
69 0 : for( size_t i=0; i<_twospin.size(); ++i )
70 0 : num*=_twospin[i]+1;
71 :
72 0 : _elem=vector<EvtComplex>( num );
73 0 : }
74 :
75 0 : EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type, const EvtComplex & val )
76 0 : {
77 : int num = 1;
78 0 : _type = type;
79 0 : _twospin=calctwospin( type );
80 :
81 0 : for( size_t i=0; i<_twospin.size(); ++i )
82 0 : num*=_twospin[i]+1;
83 :
84 0 : _elem=vector<EvtComplex>( num, val );
85 0 : }
86 :
87 0 : EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type,
88 : const vector<EvtComplex>& elem )
89 0 : {
90 : size_t num = 1;
91 :
92 0 : _type = type;
93 0 : _twospin=calctwospin( type );
94 0 : _elem=elem;
95 :
96 0 : for( size_t i=0; i<_twospin.size(); ++i ){
97 0 : num*=(_twospin[i]+1);
98 : }
99 :
100 0 : if(_elem.size() != num ) {
101 0 : report(Severity::Error,"EvtGen")<<"Wrong number of elements input:"
102 0 : <<_elem.size()<<" vs. "<<num<<endl;
103 0 : ::abort();
104 : }
105 :
106 0 : }
107 :
108 0 : EvtSpinAmp::EvtSpinAmp( const EvtSpinAmp& copy )
109 0 : {
110 0 : _twospin = copy._twospin;
111 0 : _elem = copy._elem;
112 0 : _type = copy._type;
113 0 : }
114 :
115 : void EvtSpinAmp::checktwospin( const vector<unsigned int>& twospin ) const
116 : {
117 0 : if( _twospin == twospin )
118 0 : return;
119 :
120 0 : report( Severity::Error, "EvtGen" )
121 0 : <<"Dimension or order of tensors being operated on does not match"
122 0 : <<endl;
123 0 : ::abort();
124 : }
125 :
126 : void EvtSpinAmp::checkindexargs( const vector<int>& index ) const
127 : {
128 0 : if( index.size()==0 ) {
129 0 : report(Severity::Error,"EvtGen") << "EvtSpinAmp can't handle no indices" << endl;
130 0 : ::abort();
131 : }
132 :
133 0 : if( index.size() != _twospin.size() ) {
134 0 : report( Severity::Error, "EvtGen" ) << "Rank of EvtSpinAmp index does not match: "
135 0 : <<_twospin.size()<<" expected "<<index.size()<<" input."<<endl;
136 0 : ::abort();
137 : }
138 :
139 0 : for( size_t i=0; i<_twospin.size(); ++i ) {
140 0 : if( static_cast<int>(_twospin[i])>=abs(index[i]) && static_cast<int>(_twospin[i])%2==index[i]%2 )
141 : continue;
142 0 : report(Severity::Error,"EvtGen")<<"EvtSpinAmp index out of range" << endl;
143 0 : report(Severity::Error,"EvtGen")<<" Index: ";
144 0 : for(size_t j=0; j<_twospin.size(); ++j )
145 0 : report(Severity::Error," ")<<_twospin[j];
146 :
147 0 : report(Severity::Error, " ")<<endl<<" Index: ";
148 0 : for(size_t j=0; j<index.size(); ++j )
149 0 : report(Severity::Error," ")<<index[j];
150 0 : report(Severity::Error, " ")<<endl;
151 0 : ::abort();
152 : }
153 0 : }
154 :
155 : int EvtSpinAmp::findtrueindex( const vector<int>& index ) const
156 : {
157 : int trueindex = 0;
158 :
159 0 : for( size_t i = index.size()-1; i>0; --i ) {
160 0 : trueindex += (index[i]+_twospin[i])/2;
161 0 : trueindex *= _twospin[i-1]+1;
162 : }
163 :
164 0 : trueindex += (index[0]+_twospin[0])/2;
165 :
166 0 : return trueindex;
167 : }
168 :
169 : EvtComplex & EvtSpinAmp::operator()( const vector<int>& index )
170 : {
171 0 : checkindexargs( index );
172 :
173 0 : size_t trueindex = findtrueindex(index);
174 0 : if(trueindex >= _elem.size()) {
175 0 : report(Severity::Error,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
176 0 : for(size_t i=0; i<_twospin.size(); ++i) {
177 0 : report(Severity::Error,"")<<_twospin[i]<<" ";
178 : }
179 0 : report(Severity::Error,"")<<endl;
180 :
181 0 : for(size_t i=0; i<index.size(); ++i) {
182 0 : report(Severity::Error,"")<<index[i]<<" ";
183 : }
184 0 : report(Severity::Error,"")<<endl;
185 :
186 0 : ::abort();
187 : }
188 :
189 0 : return _elem[trueindex];
190 : }
191 :
192 : const EvtComplex & EvtSpinAmp::operator()( const vector<int>& index ) const
193 : {
194 0 : checkindexargs( index );
195 :
196 0 : size_t trueindex = findtrueindex(index);
197 0 : if(trueindex >= _elem.size()) {
198 0 : report(Severity::Error,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
199 0 : for(size_t i=0; i<_twospin.size(); ++i) {
200 0 : report(Severity::Error,"")<<_twospin[i]<<" ";
201 : }
202 0 : report(Severity::Error,"")<<endl;
203 :
204 0 : for(size_t i=0; i<index.size(); ++i) {
205 0 : report(Severity::Error,"")<<index[i]<<" ";
206 : }
207 0 : report(Severity::Error,"")<<endl;
208 :
209 0 : ::abort();
210 : }
211 :
212 0 : return _elem[trueindex];
213 : }
214 :
215 : EvtComplex & EvtSpinAmp::operator()( int i, ... )
216 : {
217 0 : va_list ap;
218 0 : vector<int> index( _twospin.size() );
219 :
220 0 : va_start(ap, i);
221 :
222 0 : index[0]=i;
223 0 : for(size_t n=1; n<_twospin.size(); ++n )
224 0 : index[n]=va_arg( ap, int );
225 :
226 0 : va_end(ap);
227 :
228 0 : return (*this)( index );
229 0 : }
230 :
231 : const EvtComplex & EvtSpinAmp::operator()( int i, ... ) const
232 : {
233 0 : vector<int> index( _twospin.size() );
234 0 : va_list ap;
235 :
236 0 : va_start(ap, i);
237 :
238 0 : index[0]=i;
239 0 : for(size_t n=1; n<_twospin.size(); ++n )
240 0 : index[n]=va_arg( ap, int );
241 :
242 0 : va_end(ap);
243 :
244 0 : return (*this)( index );
245 0 : }
246 :
247 : EvtSpinAmp& EvtSpinAmp::operator=( const EvtSpinAmp& cont )
248 : {
249 0 : _twospin=cont._twospin;
250 0 : _elem=cont._elem;
251 0 : _type=cont._type;
252 :
253 0 : return *this;
254 : }
255 :
256 : EvtSpinAmp EvtSpinAmp::operator+( const EvtSpinAmp & cont ) const
257 : {
258 0 : checktwospin( cont._twospin );
259 :
260 0 : EvtSpinAmp ret( cont );
261 0 : for( size_t i=0; i<ret._elem.size(); ++i ) {
262 0 : ret._elem[i]+=_elem[i];
263 : }
264 :
265 : return ret;
266 0 : }
267 :
268 : EvtSpinAmp& EvtSpinAmp::operator+=( const EvtSpinAmp&
269 : cont )
270 : {
271 0 : checktwospin( cont._twospin );
272 :
273 0 : for( size_t i=0; i<_elem.size(); ++i )
274 0 : _elem[i]+=cont._elem[i];
275 :
276 0 : return *this;
277 : }
278 :
279 : EvtSpinAmp EvtSpinAmp::operator-( const EvtSpinAmp & cont ) const
280 : {
281 0 : checktwospin( cont._twospin );
282 :
283 0 : EvtSpinAmp ret( *this );
284 0 : for( size_t i=0; i<ret._elem.size(); ++i )
285 0 : ret._elem[i]-=cont._elem[i];
286 :
287 : return ret;
288 0 : }
289 :
290 : EvtSpinAmp& EvtSpinAmp::operator-=( const EvtSpinAmp& cont )
291 : {
292 0 : checktwospin( cont._twospin );
293 :
294 0 : for( size_t i=0; i<_elem.size(); ++i )
295 0 : _elem[i]-=cont._elem[i];
296 :
297 0 : return *this;
298 : }
299 :
300 : // amp = amp1 * amp2
301 : EvtSpinAmp EvtSpinAmp::operator*( const EvtSpinAmp & amp2 ) const
302 : {
303 0 : vector<int> index(rank()+amp2.rank());
304 0 : vector<int> index1(rank()), index2(amp2.rank());
305 0 : EvtSpinAmp amp;
306 :
307 0 : amp._twospin=_twospin;
308 0 : amp._type=_type;
309 :
310 0 : for( size_t i=0; i<amp2._twospin.size(); ++i ) {
311 0 : amp._twospin.push_back( amp2._twospin[i] );
312 0 : amp._type.push_back( amp2._type[i] );
313 : }
314 :
315 0 : amp._elem = vector<EvtComplex>( _elem.size() * amp2._elem.size() );
316 :
317 0 : for( size_t i=0; i<index1.size(); ++i )
318 0 : index[i]=index1[i]=-_twospin[i];
319 :
320 0 : for( size_t i=0; i<index2.size(); ++i )
321 0 : index[i+rank()]=index2[i]=-amp2._twospin[i];
322 :
323 0 : while(true) {
324 0 : amp( index ) = (*this)( index1 )*amp2( index2 );
325 0 : if(!amp.iterate( index )) break;
326 :
327 0 : for( size_t i=0; i<index1.size(); ++i )
328 0 : index1[i]=index[i];
329 :
330 0 : for( size_t i=0; i<index2.size(); ++i )
331 0 : index2[i]=index[i+rank()];
332 : }
333 :
334 : return amp;
335 0 : }
336 :
337 : EvtSpinAmp& EvtSpinAmp::operator*=( const EvtSpinAmp& cont )
338 : {
339 0 : EvtSpinAmp ret = (*this)*cont;
340 0 : *this=ret;
341 : return *this;
342 0 : }
343 :
344 : EvtSpinAmp& EvtSpinAmp::operator*=( const EvtComplex& real )
345 : {
346 0 : for( size_t i=0; i<_elem.size(); ++i )
347 0 : _elem[i] *= real;
348 :
349 0 : return *this;
350 : }
351 :
352 : EvtSpinAmp& EvtSpinAmp::operator/=( const EvtComplex& real )
353 : {
354 0 : for( size_t i=0; i<_elem.size(); ++i )
355 0 : _elem[i] /= real;
356 :
357 0 : return *this;
358 : }
359 :
360 : vector<int> EvtSpinAmp::iterinit() const
361 : {
362 0 : vector<int> init( _twospin.size() );
363 :
364 0 : for( size_t i=0; i<_twospin.size(); ++i )
365 0 : init[i]=-_twospin[i];
366 :
367 : return init;
368 0 : }
369 :
370 : bool EvtSpinAmp::iterate( vector<int>& index ) const
371 : {
372 0 : int last = _twospin.size() - 1;
373 :
374 0 : index[0]+=2;
375 0 : for( size_t j=0; static_cast<int>(j)<last; ++j ) {
376 0 : if( index[j] > static_cast<int>(_twospin[j]) ) {
377 0 : index[j] = -_twospin[j];
378 0 : index[j+1]+=2;
379 0 : }
380 : }
381 :
382 0 : return (abs(index[last]))<=((int)_twospin[last]);
383 : }
384 :
385 : // Test whether a particular index is an allowed one (specifically to deal with
386 : // photons and possibly neutrinos)
387 : bool EvtSpinAmp::allowed( const vector<int>& index ) const
388 : {
389 0 : if( index.size() != _type.size() ) {
390 0 : report(Severity::Error,"EvtGen")
391 0 : <<"Wrong dimensino index input to allowed."<<endl;
392 0 : ::abort();
393 : }
394 :
395 0 : for(size_t i=0; i<index.size(); ++i) {
396 0 : switch(_type[i]) {
397 : case EvtSpinType::PHOTON:
398 0 : if(abs(index[i])!=2) return false;
399 : break;
400 : case EvtSpinType::NEUTRINO:
401 0 : report(Severity::Error,"EvtGen")
402 0 : <<"EvMultibody currently cannot handle neutrinos."<<endl;
403 0 : ::abort();
404 : default:
405 : break;
406 : }
407 : }
408 :
409 0 : return true;
410 0 : }
411 :
412 : bool EvtSpinAmp::iterateallowed( vector<int>& index ) const
413 : {
414 0 : while(true) {
415 0 : if(!iterate( index ))
416 0 : return false;
417 0 : if(allowed( index ))
418 0 : return true;
419 : }
420 0 : }
421 :
422 : vector<int> EvtSpinAmp::iterallowedinit() const
423 : {
424 0 : vector<int> init = iterinit();
425 0 : while(!allowed(init)) {
426 0 : iterate(init);
427 : }
428 :
429 : return init;
430 0 : }
431 :
432 : void EvtSpinAmp::intcont( size_t a, size_t b )
433 : {
434 :
435 0 : if(rank()<=2) {
436 0 : report(Severity::Error,"EvtGen")<<"EvtSpinAmp can't handle no indices" << endl;
437 0 : ::abort();
438 : }
439 :
440 0 : size_t newrank=rank()-2;
441 :
442 0 : if(_twospin[a]!=_twospin[b]) {
443 0 : report(Severity::Error,"EvtGen")
444 0 : <<"Contaction called on indices of different dimension"
445 0 : <<endl;
446 0 : report(Severity::Error,"EvtGen")
447 0 : <<"Called on "<<_twospin[a]<<" and "<<_twospin[b]
448 0 : <<endl;
449 0 : ::abort();
450 : }
451 :
452 0 : vector<int> newtwospin( newrank );
453 0 : vector<EvtSpinType::spintype> newtype( newrank );
454 :
455 0 : for( size_t i=0, j=0; i<_twospin.size(); ++i ){
456 0 : if(i==a || i==b) continue;
457 :
458 0 : newtwospin[j] = _twospin[i];
459 0 : newtype[j] = _type[i];
460 0 : ++j;
461 0 : }
462 :
463 0 : EvtSpinAmp newamp( newtype );
464 0 : vector<int> index( rank() ), newindex = newamp.iterinit();
465 :
466 0 : for( size_t i=0; i<newrank; ++i )
467 0 : newindex[i] = -newtwospin[i];
468 :
469 0 : while(true) {
470 :
471 0 : for( size_t i=0, j=0; i<rank(); ++i ) {
472 0 : if(i==a || i==b) continue;
473 0 : index[i]=newindex[j];
474 0 : ++j;
475 0 : }
476 :
477 0 : index[b]=index[a]=-_twospin[a];
478 0 : newamp(newindex) = (*this)(index);
479 0 : for( size_t i=-_twospin[a]+2; i<=_twospin[a]; i+=2 ) {
480 0 : index[b]=index[a]=i;
481 0 : newamp(newindex) += (*this)(index);
482 : }
483 :
484 0 : if(!newamp.iterate(newindex)) break;
485 : }
486 :
487 0 : *this=newamp;
488 0 : }
489 :
490 : // In A.extcont(B), a is the index in A and b is the index in B - note that this
491 : // routine can be extremely improved!
492 : void EvtSpinAmp::extcont( const EvtSpinAmp& cont, int a, int b )
493 : {
494 0 : EvtSpinAmp ret = (*this)*cont;
495 0 : ret.intcont( a, rank()+b );
496 :
497 0 : *this=ret;
498 0 : }
|
