Problem 1101

Summary:	G4PrimaryTransformer::GimmePrimaries acceleration
Product:	Geant4	Reporter:	Pierre <pgorel>
Component:	event	Assignee:	asai
Status:	RESOLVED FIXED
Severity:	enhancement
Priority:	P5
Version:	9.3
Hardware:	All
OS:	All

Description Pierre 2010-01-27 01:25:08 CET

Hello,
in G4.9.3 (and before), the function G4PrimaryTransformer::GimmePrimaries use the following loop:

  for( G4int i=0; i<n_vertex; i++ )
  { GenerateTracks( anEvent->GetPrimaryVertex(i) ); }

with 
      inline G4PrimaryVertex* GetPrimaryVertex(G4int i=0)  const
      {
        if( i == 0 )
        { return thePrimaryVertex; }
        else if( i > 0 && i < numberOfPrimaryVertex )
        {
          G4PrimaryVertex* primaryVertex = thePrimaryVertex;
          for( G4int j=0; j<i; j++ )
          {
            if( primaryVertex == 0 ) return 0;
            primaryVertex = primaryVertex->GetNext();
          }
          return primaryVertex;
        }
        else
        { return 0; }
      }

So the loop inside GetPrimaryVertex goes through every vertex each time it is called, up to the requested one. It is not that bad with a low mulitplicity, but when it gets high, the time grows too fast.

I propose the following modification, in order to not have to use two loops :


G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent,G4int trackIDCounter)
{
  trackID = trackIDCounter;

  //TV.clearAndDestroy();
  for( size_t ii=0; ii<TV.size();ii++)
  { delete TV[ii]; }
  TV.clear();
  G4int n_vertex = anEvent->GetNumberOfPrimaryVertex();
  if(n_vertex==0) return 0;
  G4PrimaryVertex *previousVertex=0, *newVertex=0;
  for( G4int i=0; i<n_vertex; i++ )
  {
    newVertex= anEvent->GetPrimaryVertex(i, previousVertex);
    GenerateTracks( newVertex );
    previousVertex = newVertex;
  }
  return &TV;
}

and 

inline G4PrimaryVertex* GetPrimaryVertex(G4int i=0, G4PrimaryVertex* lastVertex=0)  const
{
 if( i == 0)
 { return thePrimaryVertex; }
 else if( i > 0 && i < numberOfPrimaryVertex && lastVertex != 0)
 {
   return lastVertex->GetNext();
 }
 else
 { return 0; }
}

I checked that it gives the same result (hopefully, I did not oversee anything). For a multiplicity of 100000 photons, it increased the speed of computation by a factor 15 !

Regards,
Pierre

Comment 1 Pierre 2010-01-30 00:43:03 CET

I noticed now that it may be useful to keep the possibility to have no "lastVertex" as an argument (at least for backward compatibility). So I would propose the following :

inline G4PrimaryVertex* GetPrimaryVertex(G4int i=0, G4PrimaryVertex* lastVertex=0)  const
      {
        if( i == 0)
        { return thePrimaryVertex; }
        else
          if( i > 0 && i < numberOfPrimaryVertex) {
            if (lastVertex != 0){
              return lastVertex->GetNext();
            }
            else {
              G4PrimaryVertex* primaryVertex = thePrimaryVertex;
              for( G4int j=0; j<i; j++ )
              {
                if( primaryVertex == 0 ) return 0;
                primaryVertex = primaryVertex->GetNext();
              }
              return primaryVertex;
            }
        }
        else
        { return 0; }
      }

Comment 2 asai 2010-04-19 21:46:19 CEST

Dear Pierre,

Thanks for pointing this out.
The implementation of G4PrimaryTransporfer will be changed as following,
which I believe has the same effect as you proposed, without changes in
other class.

G4TrackVector* G4PrimaryTransformer::GimmePrimaries
                          (G4Event* anEvent,G4int trackIDCounter)
{
  trackID = trackIDCounter;

  //TV.clearAndDestroy();
  for( size_t ii=0; ii<TV.size();ii++)
  { delete TV[ii]; }
  TV.clear();
  G4PrimaryVertex* nextVertex = anEvent->GetPrimaryVertex();
  while(nextVertex)
  { 
    GenerateTracks(nextVertex);
    nextVertex = nextVertex->GetNext();
  }
  return &TV;
}

This fix will show up in the next release.

Kind regards,
Makoto