title: Multithreaded code uses only one CPU thread while the number of threads should be more
body: <p>I am trying to implement multithreading to a CPU intensive simulation that normally takes around 10hrs to run. I wrote a dummy code which works fine for multithreading and uses the specified number of CPU threads. But when I implement the same thing in my main program, it uses only one thread. What could be the reason for this?</p>
  <p>My code:
  <strong>photongen()</strong></p>
  <pre><code>void photongen() {
  
      // Iterating over Muon list
      for (auto m = muonList.begin(); m != muonList.end(); ++m) {
          
          nP = 0; // Initialise the number of photons to 0
  
          vector&lt;Photon&gt; tempP;   // Temporary vector to store generated photons
          
          // Muon propagation
          bool stopLoop = false;
          while (!stopLoop) {
              
              // Photon generation
              if (nPhotons()) {
                  for (int i = 0; i &lt; nPhotons(); i++){
  
                      // Create photon instance
                      Photon *P = new Photon(m-&gt;x, m-&gt;y, m-&gt;z, m-&gt;run, m-&gt;event, i);
                      tempP.push_back(*P);
  
                      delete P;   // Free memory by photon object
                  }
              }
  
              // Check if muon has left the detector
              if (m-&gt;x &gt; S_LENGTH/2. - PATH_STEP*sin(m-&gt;theta)*cos(m-&gt;phi) || m-&gt;x &lt; -S_LENGTH/2. - PATH_STEP*sin(m-&gt;theta)*cos(m-&gt;phi)) {
                  stopLoop = true;
              } else if (m-&gt;y &gt; S_WIDTH/2. - PATH_STEP*sin(m-&gt;theta)*sin(m-&gt;phi) || m-&gt;y &lt; -S_WIDTH/2. - PATH_STEP*sin(m-&gt;theta)*sin(m-&gt;phi)) {
                  stopLoop = true;
              } else if (m-&gt;z &gt; S_HEIGHT - PATH_STEP*cos(m-&gt;theta)) {
                  stopLoop = true;
              } else {
                  // Update Muon path
                  m-&gt;x += PATH_STEP*sin(m-&gt;theta)*cos(m-&gt;phi);
                  m-&gt;y += PATH_STEP*sin(m-&gt;theta)*sin(m-&gt;phi);
                  m-&gt;z += PATH_STEP*cos(m-&gt;theta);
              }
          }
  
          multiprocess(8, tempP, &amp;propagation);
          tempP.clear();
      }
  }
  </code></pre>
  <p><strong>propagation()</strong></p>
  <pre><code>void propagation(Photon *P) {
      if (photonProp(P)) {
          pLock.lock();
          photonList.push_back(*P);
          nP++;
          pLock.unlock();
      }
  }
  </code></pre>
  <p><strong>photonProp()</strong></p>
  <pre><code>bool photonProp(Photon*&amp; P) {
  
      float dist[3];
      for (int i = 0; i &lt; MAX_REFL; i++) {
          
          // Precalculation for optimisation
          float theta = P-&gt;theta;
          float phi = P-&gt;phi;
  
          float sintheta = sin(theta);
          float costheta = cos(theta);
          float sinphi = sin(phi);
          float cosphi = cos(phi);
  
          float P_x = P-&gt;x;
          float P_y = P-&gt;y;
          float P_z = P-&gt;z;
  
          // Find plane of intersection of particle
          dist[0] = (sintheta*cosphi &gt; 0) ? (S_LENGTH/2. - P_x)/(sintheta*cosphi) : (-S_LENGTH/2. - P_x)/(sintheta*cosphi);
          dist[1] = (sintheta*sinphi &gt; 0) ? (S_WIDTH/2. - P_y)/(sintheta*sinphi) : (-S_WIDTH/2. - P_y)/(sintheta*sinphi);
          dist[2] = (costheta &gt; 0) ? (S_HEIGHT - P_z)/costheta : -P_z/costheta;
  
          // Calculate the plane of intersection
          int min_index = distance(dist,min_element(dist, dist + 3));
          float d = dist[min_index]; // distance from original point to plane
  
          // Check for attenuation
          float lambda_sc = 100.;
          float x_surv = -lambda_sc*log(r-&gt;Uniform());
          if (d &gt; x_surv) return false;
          
          // Update particle coords
          P-&gt;x = P_x + d*sintheta*cosphi;
          P-&gt;y = P_y + d*sintheta*sinphi;
          P-&gt;z = P_z + d*costheta;
  
          // check if particle hit photodetector
          if (min_index  == 0 &amp;&amp; sintheta*cosphi &gt; 0 &amp;&amp; P-&gt;y &lt;= -S_WIDTH/2. + 1) {
              return true;
          }
  
          // Check for reflective loss
          if (r-&gt;Uniform() &lt; 0.05) return false;
  
          // reflect particle
          if (min_index  == 0) {
              P-&gt;phi = M_PI - phi;
          } else if (min_index == 1) {
              P-&gt;phi = -phi;
          } else if (min_index == 2) {
              P-&gt;theta = -theta;
          }
      }
      
      return false;
  }
  </code></pre>
  <p><strong>multiprocess()</strong></p>
  <pre><code>template &lt;typename T&gt;
  void multiprocess(int nThreads, vector&lt;T&gt; arg, void (*func)(T*)) {
  
      size_t i;
  
      for (i=0; i &lt; arg.size() - nThreads; i += nThreads) {
          
          vector&lt;thread&gt; t(nThreads);
  
          for (int j = 0; j &lt; nThreads; j++){
              t[j] = thread(func, &amp;arg[i+j]);
          }
  
          for (int j = 0; j &lt; nThreads; j++){
              t[j].join();
          }
      }
  
      vector&lt;thread&gt; t(arg.size()%nThreads);
  
      for (std::size_t j = 0; j &lt; arg.size()%nThreads; j++){ 
          t[j] = thread(func, &amp;arg[i+j]);        
      }
  
      for (std::size_t j = 0; j &lt; arg.size()%nThreads; j++){
          t[j].join();
      } 
  }
  </code></pre>