C++ Stopwatch::toc方法代码示例

int main(int argc, char** argv) {
    Stopwatch sw;

    sw.tic();
    ParameterReader* paraRdr = new ParameterReader();
    paraRdr->readFromFile("parameters.dat");
    paraRdr->readFromArguments(argc, argv);

    // create integration grid along eta direction for boost-invariant medium
    int neta = paraRdr->getVal("neta");
    double eta_i = paraRdr->getVal("eta_i");
    double eta_f = paraRdr->getVal("eta_f");
    double* eta_ptr = new double[neta];
    double* etaweight_ptr = new double[neta];
    gauss_quadrature(neta, 1, 0.0, 0.0, eta_i, eta_f, eta_ptr, etaweight_ptr);

    PhotonEmission thermalPhotons(paraRdr);

    // initialize hydro medium
    int hydro_flag = paraRdr->getVal("hydro_flag");
    if (hydro_flag == 0) {
        int bufferSize = paraRdr->getVal("HydroinfoBuffersize");
        int hydroInfoVisflag = paraRdr->getVal("HydroinfoVisflag");
        // hydro data file pointer
        HydroinfoH5* hydroinfo_ptr = new HydroinfoH5(
                        "results/JetData.h5", bufferSize, hydroInfoVisflag);
        // calculate thermal photons from the hydro medium
        thermalPhotons.calPhotonemission(hydroinfo_ptr, eta_ptr, 
                                         etaweight_ptr);
        delete hydroinfo_ptr;
    } else if (hydro_flag == 1) {
        Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
        int hydro_mode = 8;
        int nskip_tau = paraRdr->getVal("hydro_nskip_tau");
        hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau);
        // calculate thermal photons from the hydro medium
        thermalPhotons.calPhotonemission(hydroinfo_ptr, eta_ptr, 
                                         etaweight_ptr);
        delete hydroinfo_ptr;
    } else if (hydro_flag == 3) {
        Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
        int hydro_mode = 9;
        int nskip_tau = paraRdr->getVal("hydro_nskip_tau");
        hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau);
        // calculate thermal photons from the hydro medium
        thermalPhotons.calPhotonemission(hydroinfo_ptr, eta_ptr, 
                                         etaweight_ptr);
        delete hydroinfo_ptr;
    } else if (hydro_flag == 2) {
        Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
        int hydro_mode = 10;
        int nskip_tau = 1;
        hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau);
        // calculate thermal photons from the hydro medium
        thermalPhotons.calPhotonemission_3d(hydroinfo_ptr);
        delete hydroinfo_ptr;
    } else {
        cout << "main: unrecognized hydro_flag = " << hydro_flag << endl;
        exit(1);
    }

    // sum up all channels and compute thermal photon spectra and vn
    thermalPhotons.calPhoton_SpvnpT_individualchannel();
    thermalPhotons.calPhoton_total_SpMatrix();
    thermalPhotons.calPhoton_total_Spvn();

    // output results
    thermalPhotons.outputPhotonSpvn();

    sw.toc();
    cout << "totally takes : " << sw.takeTime() << " seconds." << endl;

    // clean up
    delete [] eta_ptr;
    delete [] etaweight_ptr;

    return(0);
}

//.........这里部分代码省略.........
  double normalization[pT_tab_length]; // normalization factor
  for (int i=0; i<pT_tab_length; i++) normalization[i] = 0.0;

  double vn[pT_tab_length][number_of_flows][2]; // diff_flow numerators; 2: 0,1->real,imag
  for (int i=0; i<pT_tab_length; i++)
  for (int t=0; t<number_of_flows; t++)
    {vn[i][t][0]=0; vn[i][t][1]=0;}

  for (int i=0; i<pT_tab_length; i++)
  //for (int i=0; i<1; i++) // for debugging
  {
    double pT = pT_tab->get(1,i+1); // pT_weight = pT_tab->get(2,i+1);
    double mT = sqrt(mass*mass + pT*pT);

    // phi integration
    for(int j=0; j<phi_tab_length; j++)
    //for(int j=0; j<1; j++) // for debugging
    {
      double phi = phi_tab->get(1,j+1), phi_weight = phi_tab->get(2,j+1);
      double dE = dE_ptdptdphidy->get(i+1,j+1);

      normalization[i] += dE*phi_weight;
      for (int order=from_order; order<=to_order; order++)
      {
        vn[i][order-from_order][0] += dE*phi_weight*cos(order*phi);
        vn[i][order-from_order][1] += dE*phi_weight*sin(order*phi);
      }
    }

    normalization[i] = normalization[i] + 1e-30;
    // store values
    vn_diff.set(1, i+1, pT);
    vn_diff.set(2, i+1, mT-mass);
    vn_diff.set(3, i+1, normalization[i]/(2.0*M_PI)); // 2*pi: azimuthal angle averaged
    for (int t=0; t<number_of_flows; t++)
    {
      vn_diff.set(4+t*3, i+1, vn[i][t][0]/normalization[i]);
      vn_diff.set(5+t*3, i+1, vn[i][t][1]/normalization[i]);
      vn_diff.set(6+t*3, i+1, sqrt(vn[i][t][0]*vn[i][t][0]+vn[i][t][1]*vn[i][t][1])/normalization[i]);
    }

  }
  //cout << "done." << endl;


  //---------------------
  // integrated flow
  //---------------------
  //cout << "Calculating integrated flows... ";

  double normalizationi = 0;

  double vni[number_of_flows][2]; // integrated_flow numerators; 2: 0,1->real,imag
  for (int t=0; t<number_of_flows; t++) {vni[t][0]=0; vni[t][1]=0;}

  for (int i=0; i<pT_tab_length; i++)
  //for (int i=0; i<1; i++) // for debugging
  {
    double pT = pT_tab->get(1,i+1), pT_weight = pT_tab->get(2,i+1);

    normalizationi += normalization[i]*pT*pT_weight;

    for (int order=from_order; order<=to_order; order++)
    {
      vni[order-from_order][0] += vn[i][order-from_order][0]*pT*pT_weight;
      vni[order-from_order][1] += vn[i][order-from_order][1]*pT*pT_weight;
    }

  }

  vn_inte.set(1, 1, 0);
  vn_inte.set(2, 1, normalizationi);
  vn_inte.set(3, 1, 0);
  vn_inte.set(4, 1, 1);
  vn_inte.set(5, 1, 0);
  vn_inte.set(6, 1, 1);

  for (int t=0; t<number_of_flows; t++)
  {
    vn_inte.set(1, t+2, from_order+t);
    vn_inte.set(2, t+2, vni[from_order+t-1][0]);
    vn_inte.set(3, t+2, vni[from_order+t-1][1]);
    vn_inte.set(4, t+2, vni[from_order+t-1][0]/normalizationi);
    vn_inte.set(5, t+2, vni[from_order+t-1][1]/normalizationi);
    vn_inte.set(6, t+2, sqrt(vni[from_order+t-1][0]*vni[from_order+t-1][0]+vni[from_order+t-1][1]*vni[from_order+t-1][1])/normalizationi);
  }

  // save to files
  ofstream of1(flow_differential_filename_in.c_str(), ios_base::app);
  vn_diff.printTable(of1);
  of1.close();

  ofstream of2(flow_integrated_filename_in.c_str(), ios_base::app);
  vn_inte.printTable(of2);
  of2.close();

  sw.toc();
  //cout << "calculate_flows finishes " << sw.takeTime() << " seconds." << endl;

}

//*********************************************************************************************
void EmissionFunctionArray::calculate_dN_ptdptdphidy_and_flows_4all(int to_order)
// Calculate dNArrays and flows for all particles given in chosen_particle file.
{
    if (USE_HISTORIC_FORMAT)
    {
          cout << endl
               <<"**********************************************************"
               << endl
               << "Function calculate_dN_ptdptdphidy_and_flows_4all(old) started... " << endl;
          Stopwatch sw;
          sw.tic();

          remove(dN_ptdptdphidy_filename.c_str());
          remove(flow_differential_filename_old.c_str());
          remove(flow_integrated_filename_old.c_str());
          if(CALCULATEDED3P)
          {
             remove(dE_ptdptdphidy_filename.c_str());
             remove(energyflow_differential_filename_old.c_str());
             remove(energyflow_integrated_filename_old.c_str());
          }

          particle_info* particle = NULL;

          for (int n=0; n<Nparticles; n++)
          {
            particle = &particles[n];
            cout << "Index: " << n << ", Name: " << particle->name << ", Monte-carlo index: " << particle->monval;

            // first, dN_xxx arrays:
            if (chosen_particles_01_table[n]==0)
            {
              cout << " -- Skipped." << endl;
              dN_ptdptdphidy->setAll(0.0);
              if(CALCULATEDED3P) dE_ptdptdphidy->setAll(0.0);
              last_particle_idx = n; // fake a "calculation"
            }
            else
            {
              cout << " -- Processing... " << endl;
              calculate_dN_ptdptdphidy(n);
            }
            write_dN_ptdptdphidy_toFile();

            // next flows:
            ofstream of1(flow_differential_filename_old.c_str(), ios_base::app);
            of1 << "# Output for particle: " << particle->name << endl;
            of1 << "#                 " << particle->monval << endl;
            of1.close();

            ofstream of2(flow_integrated_filename_old.c_str(), ios_base::app);
            of2 << "# For: " << particle->name << endl;
            of2.close();
            calculate_flows(to_order, flow_differential_filename_old, flow_integrated_filename_old);
            if(CALCULATEDED3P)
            {
               ofstream of1(energyflow_differential_filename_old.c_str(), ios_base::app);
               of1 << "# Output for particle: " << particle->name << endl;
               of1 << "#                 " << particle->monval << endl;
               of1.close();
               ofstream of2(energyflow_integrated_filename_old.c_str(), ios_base::app);
               of2 << "# For: " << particle->name << endl;
               of2.close();
               calculate_Energyflows(to_order, energyflow_differential_filename_old, energyflow_integrated_filename_old);
            }
          }


          sw.toc();
          cout << "calculate_dN_ptdptdphidy_and_flows_4all finishes " << sw.takeTime() << " seconds." << endl;
    }
    else
    {
            cout << endl
            << "****************************************************************"
            << endl
            << "Function calculate_dN_ptdptdphidy_and_flows_4all(new) started... " << endl;
            Stopwatch sw;
            sw.tic();

            // prepare a huge array to store calculated dN_ptdptdphidy
            Table* dNs[Nparticles];
            for (int n=0; n<Nparticles; n++) dNs[n]=NULL;
            
            Table* dEs[Nparticles];
            for (int n=0; n<Nparticles; n++) dEs[n]=NULL;

            // loop over chosen particles
            particle_info* particle = NULL;
            for (int m=0; m<number_of_chosen_particles; m++)
            {
                int particle_idx = chosen_particles_sampling_table[m];
                particle = &particles[particle_idx];
                int monval = particle->monval;
                cout << "Index: " << m << ", Name: " << particle->name << ", Monte-carlo index: " << monval << endl;
                // Calculate dN / (ptdpt dphi dy)
                if (m>0 && particles_are_the_same(particle_idx, chosen_particles_sampling_table[m-1]))
                {
                   cout << " -- Using dN_ptdptdphidy from previous calculation... " << endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          //double cos_phi = cos(phi);
          //double sin_phi = sin(phi);
          //double px = pT*cos_phi;
          //double py = pT*sin_phi;
          // new way
          double px = pT*trig_phi_table[j][0];
          double py = pT*trig_phi_table[j][1];

          double dN_ptdptdphidy_tmp = 0.0;
          double dE_ptdptdphidy_tmp = 0.0;

          for (long l=0; l<FO_length; l++)
          {
              surf = &FOsurf_ptr[l];

              double Tdec = surf->Tdec;
              double Pdec = surf->Pdec;
              double Edec = surf->Edec;
              double deltaf_prefactor = 1.0/(2.0*Tdec*Tdec*(Edec+Pdec))*INCLUDE_DELTAF;
              double mu = surf->particle_mu[last_particle_idx];

              double tau = surf->tau;
              double vx = surf->vx;
              double vy = surf->vy;

              double da0 = surf->da0;
              double da1 = surf->da1;
              double da2 = surf->da2;
              double pi00 = surf->pi00;
              double pi01 = surf->pi01;
              double pi02 = surf->pi02;
              double pi11 = surf->pi11;
              double pi12 = surf->pi12;
              double pi22 = surf->pi22;
              double pi33 = surf->pi33;

              double v2 = vx*vx + vy*vy;
              double gammaT = 1.0/sqrt(1.0 - v2);


              for (int k=0; k<eta_tab_length; k++)
              {
                  //double eta_s = eta_tab->get(1,k+1); // unused

                  // old way
                  //double delta_eta = eta_tab->get(2,k+1);
                  // new way
                  double delta_eta = delta_eta_tab[k];

                  // old way
                  //double cosh_y_eta = cosh(y-eta_s);
                  //double sinh_y_eta = sinh(y-eta_s);

                  // old way
                  //double pt = mT*cosh_y_eta;
                  //double pz = mT*sinh_y_eta;

                  // new way
                  double pt = mT*hypertrig_etas_table[k][0];
                  double pz = mT*hypertrig_etas_table[k][1];

                  double expon = (gammaT*(pt*1 - px*vx - py*vy) - mu) / Tdec;
                  double f0 = 1./(exp(expon)+sign);       //thermal equilibrium distributions

                  // Must adjust this to be correct for the p*del \tau term.  The plus sign is
                  // due to the fact that the DA# variables are for the covariant surface integration
                  double pdsigma = pt*da0 + px*da1 + py*da2;

                  //viscous corrections
                  double Wfactor = pt*pt*pi00 - 2.0*pt*px*pi01 - 2.0*pt*py*pi02 + px*px*pi11 + 2.0*px*py*pi12 + py*py*pi22 + pz*pz*pi33;
                  double deltaf = (1 - F0_IS_NOT_SMALL*sign*f0)*Wfactor*deltaf_prefactor;
                  double result;
                  if(1+deltaf < 0.0) //set results to zero when delta f turns whole expression to negative
                     result = 0.0;
                  else
                     result = prefactor*degen*f0*(1+deltaf)*pdsigma*tau;

                  dN_ptdptdphidy_tmp += result*delta_eta;
                  if(CALCULATEDED3P) dE_ptdptdphidy_tmp += result*delta_eta*pt;
              } // k
          } // l

          dN_ptdptdphidy_tab[i][j] = dN_ptdptdphidy_tmp;
          if (CALCULATEDED3P) dE_ptdptdphidy_tab[i][j] = dE_ptdptdphidy_tmp;
          if (AMOUNT_OF_OUTPUT>0) print_progressbar((i*phi_tab_length+j)/progress_total);
      }
  //cout << int(100.0*(i+1)/pT_tab_length) << "% completed" << endl;
  }
  if (AMOUNT_OF_OUTPUT>0) print_progressbar(1);

  for (int i=0; i<pT_tab_length; i++)
  for (int j=0; j<phi_tab_length; j++)
  {
    dN_ptdptdphidy->set(i+1,j+1,dN_ptdptdphidy_tab[i][j]);
    if(CALCULATEDED3P) dE_ptdptdphidy->set(i+1,j+1,dE_ptdptdphidy_tab[i][j]);
  }

  sw.toc();
  cout << endl << "Finished " << sw.takeTime() << " seconds." << endl;
}