6464#include < thrust/extrema.h>
6565#include < thrust/count.h>
6666#include < thrust/fill.h>
67+ #include < thrust/sort.h>
68+ #include < thrust/iterator/counting_iterator.h>
6769
6870#if !(THRUST_DEVICE_SYSTEM==THRUST_DEVICE_BACKEND_OMP || THRUST_DEVICE_SYSTEM==THRUST_DEVICE_BACKEND_TBB)
6971#include < thrust/system/cuda/execution_policy.h>
@@ -171,7 +173,8 @@ class PhaseSpace {
171173 */
172174 PhaseSpace (GReal_t _MotherMass, vector<GReal_t> _Masses, GLong_t _NEvents) :
173175 fNDaughters (_Masses.size()), fNEvents (_NEvents),
174- fRandNumbers ((3 * _Masses.size() - 2 ) * _NEvents, 0.0 ), fMaxWeight (0.0 ),
176+ fSeed (1 ),
177+ fMaxWeight (0.0 ),
175178 RND_Time (0.0 ), EVT_Time(0.0 ), EXP_Time(0.0 ), fNAccepted (0 )
176179 {
177180
@@ -224,8 +227,8 @@ class PhaseSpace {
224227
225228 fWeights .clear ();
226229 fWeights .shrink_to_fit ();
227- fRandNumbers .clear ();
228- fRandNumbers .shrink_to_fit ();
230+ // fRandNumbers.clear();
231+ // fRandNumbers.shrink_to_fit();
229232 fAccRejFlags .clear ();
230233 fAccRejFlags .shrink_to_fit ();
231234 }
@@ -294,6 +297,14 @@ class PhaseSpace {
294297 return fMaxWeight ;
295298 }
296299
300+ inline GInt_t GetSeed () const {
301+ return fSeed ;
302+ }
303+
304+ inline void SetSeed (GInt_t _seed) {
305+ fSeed =_seed;
306+ }
307+
297308 /* *
298309 * Export the events and all related information to host.
299310 */
@@ -334,6 +345,7 @@ class PhaseSpace {
334345
335346 GLong_t fNEvents ; // /< Number of events.
336347 GInt_t fNDaughters ;// /< Number of daughters.
348+ GInt_t fSeed ;// /< seed.
337349 GLong_t fNAccepted ;
338350 GReal_t RND_Time;// /< Random number generation time interval seconds.
339351 GReal_t EVT_Time;// /< Event generation time interval in seconds.
@@ -344,7 +356,6 @@ class PhaseSpace {
344356 RealVector_d fWeights ;// /< Device vector of weights.
345357 BoolVector_d fAccRejFlags ;// /< Device vector of Accept/reject flags
346358 Particles_d fDaughters [kMAXP ];// /< Array of device vectors with the daughter four-vectors
347- RealVector_d fRandNumbers ;// /<
348359
349360
350361};
@@ -363,10 +374,13 @@ GULong_t PhaseSpace::Unweight()
363374 }
364375 else
365376 {
366- mc_device_vector<GLong_t> Evt (fNEvents );
367- thrust::sequence (Evt.begin (), Evt.end ());
368377
369- thrust::transform (Evt.begin (), Evt.end (), fWeights .begin (),
378+ // create iterators
379+ thrust::counting_iterator<GLong_t> first (0 );
380+ thrust::counting_iterator<GLong_t> last = first + fNEvents ;
381+
382+
383+ thrust::transform (first, last, fWeights .begin (),
370384 fAccRejFlags .begin (), FlagAcceptReject (fMaxWeight ));
371385
372386 count = thrust::count (fAccRejFlags .begin (), fAccRejFlags .end (),
@@ -534,20 +548,11 @@ void PhaseSpace::Generate(const Vector4R fMother) {
534548 cudaDeviceSetCacheConfig (cudaFuncCachePreferL1);
535549#endif
536550 /* random number generation */
537- timespec time_rnd_start, time_rnd_end;
538- clock_gettime (TIMER, &time_rnd_start);
539-
540- mc_device_vector<GLong_t> Evt (fNEvents );
541- thrust::sequence (Evt.begin (), Evt.end ());
542551
543- thrust::for_each (Evt.begin (), Evt.end (),
544- RandGen (fNDaughters ,
545- (thrust::raw_pointer_cast (fRandNumbers .data ()))));
546-
547- clock_gettime (TIMER, &time_rnd_end);
548-
549- RND_Time = ((GReal_t) (time_diff (time_rnd_start, time_rnd_end).tv_sec
550- + time_diff (time_rnd_start, time_rnd_end).tv_nsec * 1.0e-9 ));
552+ RND_Time = 0 ;
553+ // create iterators
554+ thrust::counting_iterator<GLong_t> first (0 );
555+ thrust::counting_iterator<GLong_t> last = first + fNEvents ;
551556
552557 // Vai!!!
553558
@@ -559,92 +564,92 @@ void PhaseSpace::Generate(const Vector4R fMother) {
559564
560565 case 2 :
561566
562- thrust::transform (Evt. begin (), Evt. end () ,
567+ thrust::transform (first, last ,
563568 thrust::make_zip_iterator (
564569 thrust::make_tuple (fDaughters [0 ].begin (),
565570 fDaughters [1 ].begin ())), fWeights .begin (),
566- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
571+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
567572
568573 break ;
569574
570575 case 3 :
571576
572- thrust::transform (Evt. begin (), Evt. end () ,
577+ thrust::transform (first, last ,
573578 thrust::make_zip_iterator (
574579 thrust::make_tuple (fDaughters [0 ].begin (),
575580 fDaughters [1 ].begin (), fDaughters [2 ].begin ())),
576581 fWeights .begin (),
577- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
582+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
578583
579584 break ;
580585 case 4 :
581586
582- thrust::transform (Evt. begin (), Evt. end () ,
587+ thrust::transform (first, last ,
583588 thrust::make_zip_iterator (
584589 thrust::make_tuple (fDaughters [0 ].begin (),
585590 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
586591 fDaughters [3 ].begin ())), fWeights .begin (),
587- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
592+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
588593
589594 break ;
590595 case 5 :
591596
592- thrust::transform (Evt. begin (), Evt. end () ,
597+ thrust::transform (first, last ,
593598 thrust::make_zip_iterator (
594599 thrust::make_tuple (fDaughters [0 ].begin (),
595600 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
596601 fDaughters [3 ].begin (), fDaughters [4 ].begin ())),
597602 fWeights .begin (),
598- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
603+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
599604
600605 // }
601606 break ;
602607 case 6 :
603608
604- thrust::transform (Evt. begin (), Evt. end () ,
609+ thrust::transform (first, last ,
605610 thrust::make_zip_iterator (
606611 thrust::make_tuple (fDaughters [0 ].begin (),
607612 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
608613 fDaughters [3 ].begin (), fDaughters [4 ].begin (),
609614 fDaughters [5 ].begin ())), fWeights .begin (),
610- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
615+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
611616
612617 break ;
613618 case 7 :
614619
615- thrust::transform (Evt. begin (), Evt. end () ,
620+ thrust::transform (first, last ,
616621 thrust::make_zip_iterator (
617622 thrust::make_tuple (fDaughters [0 ].begin (),
618623 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
619624 fDaughters [3 ].begin (), fDaughters [4 ].begin (),
620625 fDaughters [5 ].begin (), fDaughters [6 ].begin ())),
621626 fWeights .begin (),
622- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
627+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
623628
624629 break ;
625630 case 8 :
626631
627- thrust::transform (Evt. begin (), Evt. end () ,
632+ thrust::transform (first, last ,
628633 thrust::make_zip_iterator (
629634 thrust::make_tuple (fDaughters [0 ].begin (),
630635 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
631636 fDaughters [3 ].begin (), fDaughters [4 ].begin (),
632637 fDaughters [5 ].begin (), fDaughters [6 ].begin (),
633638 fDaughters [7 ].begin ())), fWeights .begin (),
634- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
639+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
635640
636641 break ;
637642 case 9 :
638643
639- thrust::transform (Evt. begin (), Evt. end () ,
644+ thrust::transform (first, last ,
640645 thrust::make_zip_iterator (
641646 thrust::make_tuple (fDaughters [0 ].begin (),
642647 fDaughters [1 ].begin (), fDaughters [2 ].begin (),
643648 fDaughters [3 ].begin (), fDaughters [4 ].begin (),
644649 fDaughters [5 ].begin (), fDaughters [6 ].begin (),
645650 fDaughters [7 ].begin (), fDaughters [8 ].begin ())),
646651 fWeights .begin (),
647- DecayMother (fMother , fMasses , fNDaughters , fRandNumbers ));
652+ DecayMother (fMother , fMasses , fNDaughters , fSeed ));
648653
649654 break ;
650655
@@ -674,118 +679,120 @@ void PhaseSpace::Generate(Particles_d fMothers) {
674679 cout << " fNEvents != fMothers.size()" << endl;
675680
676681 /* random number generation */
682+ /*
677683 timespec time_rnd_start, time_rnd_end;
678684 clock_gettime(TIMER, &time_rnd_start);
679685
680- mc_device_vector<GLong_t> Evt (fNEvents );
681- thrust::sequence (Evt.begin (), Evt.end ());
682-
683- thrust::for_each (Evt.begin (), Evt.end (),
684- RandGen (fNDaughters ,
685- (thrust::raw_pointer_cast (fRandNumbers .data ()))));
686686
687687 clock_gettime(TIMER, &time_rnd_end);
688688
689689 RND_Time = ((GReal_t) (time_diff(time_rnd_start, time_rnd_end).tv_sec
690690 + time_diff(time_rnd_start, time_rnd_end).tv_nsec * 1.0e-9));
691+ */
691692
692693 /* event generation */
693694 timespec time_event_start, time_event_end;
694695 clock_gettime (TIMER, &time_event_start);
695696
697+ RND_Time = 0.0 ;
698+ // create iterators
699+ thrust::counting_iterator<GLong_t> first (0 );
700+ thrust::counting_iterator<GLong_t> last = first + fNEvents ;
701+
702+
696703 switch (fNDaughters ) {
697704
698705 case 2 :
699706
700- thrust::transform (Evt. begin (), Evt. end () ,
707+ thrust::transform (first, last ,
701708 thrust::make_zip_iterator (
702709 thrust::make_tuple (fMothers .begin (),
703710 fDaughters [0 ].begin (), fDaughters [1 ].begin ())),
704711 fWeights .begin (),
705- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
712+ DecayMothers (fMasses , fNDaughters , fSeed ));
706713
707714 break ;
708715
709716 case 3 :
710717
711- thrust::transform (Evt. begin (), Evt. end () ,
718+ thrust::transform (first, last ,
712719 thrust::make_zip_iterator (
713720 thrust::make_tuple (fMothers .begin (),
714721 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
715722 fDaughters [2 ].begin ())), fWeights .begin (),
716- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
723+ DecayMothers (fMasses , fNDaughters , fSeed ));
717724
718725 break ;
719726 case 4 :
720727
721- thrust::transform (Evt. begin (), Evt. end () ,
728+ thrust::transform (first, last ,
722729 thrust::make_zip_iterator (
723730 thrust::make_tuple (fMothers .begin (),
724731 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
725732 fDaughters [2 ].begin (), fDaughters [3 ].begin ())),
726733 fWeights .begin (),
727- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
734+ DecayMothers (fMasses , fNDaughters , fSeed ));
728735
729736 break ;
730737 case 5 :
731738
732- thrust::transform (Evt. begin (), Evt. end () ,
739+ thrust::transform (first, last ,
733740 thrust::make_zip_iterator (
734741 thrust::make_tuple (fMothers .begin (),
735742 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
736743 fDaughters [2 ].begin (), fDaughters [3 ].begin (),
737744 fDaughters [4 ].begin ())), fWeights .begin (),
738- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
745+ DecayMothers (fMasses , fNDaughters , fSeed ));
739746
740747 break ;
741748 case 6 :
742749
743- thrust::transform (Evt. begin (), Evt. end () ,
750+ thrust::transform (first, last ,
744751 thrust::make_zip_iterator (
745752 thrust::make_tuple (fMothers .begin (),
746753 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
747754 fDaughters [2 ].begin (), fDaughters [3 ].begin (),
748755 fDaughters [4 ].begin (), fDaughters [5 ].begin ())),
749756 fWeights .begin (),
750- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
757+ DecayMothers (fMasses , fNDaughters , fSeed ));
751758
752759 break ;
753760 case 7 :
754761
755- thrust::transform (Evt. begin (), Evt. end () ,
762+ thrust::transform (first, last ,
756763 thrust::make_zip_iterator (
757764 thrust::make_tuple (fMothers .begin (),
758765 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
759766 fDaughters [2 ].begin (), fDaughters [3 ].begin (),
760767 fDaughters [4 ].begin (), fDaughters [5 ].begin (),
761768 fDaughters [6 ].begin ())), fWeights .begin (),
762- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
769+ DecayMothers (fMasses , fNDaughters , fSeed ));
763770
764771 break ;
765772 case 8 :
766773
767- thrust::transform (Evt. begin (), Evt. end () ,
774+ thrust::transform (first, last ,
768775 thrust::make_zip_iterator (
769776 thrust::make_tuple (fMothers .begin (),
770777 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
771778 fDaughters [2 ].begin (), fDaughters [3 ].begin (),
772779 fDaughters [4 ].begin (), fDaughters [5 ].begin (),
773780 fDaughters [6 ].begin (), fDaughters [7 ].begin ())),
774781 fWeights .begin (),
775- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
782+ DecayMothers (fMasses , fNDaughters , fSeed ));
776783
777784 break ;
778785 case 9 :
779786
780- thrust::transform (Evt. begin (), Evt. end () ,
787+ thrust::transform (first, last ,
781788 thrust::make_zip_iterator (
782789 thrust::make_tuple (fMothers .begin (),
783790 fDaughters [0 ].begin (), fDaughters [1 ].begin (),
784791 fDaughters [2 ].begin (), fDaughters [3 ].begin (),
785792 fDaughters [4 ].begin (), fDaughters [5 ].begin (),
786793 fDaughters [6 ].begin (), fDaughters [7 ].begin (),
787794 fDaughters [8 ].begin ())), fWeights .begin (),
788- DecayMothers (fMasses , fNDaughters , fRandNumbers ));
795+ DecayMothers (fMasses , fNDaughters , fSeed ));
789796
790797 break ;
791798
0 commit comments