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

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);
}

//***************************************************************************
void EmissionFunctionArray::calculate_Energyflows(int to_order, string flow_differential_filename_in, string flow_integrated_filename_in)
// Calculate flow from order from_order to to_order and store them to files.
{
  Stopwatch sw;
  sw.tic();

  int from_order = 1;

  int number_of_flows = to_order-from_order+1;

  Table vn_diff(3+number_of_flows*3, pT_tab_length); // line format: pT, mT, dN/(pT dpT), flow_1_real, flow_1_imag, flow_1_norm, ...
  Table vn_inte(6, to_order+1); // line format: order# (starting from 0), numerator_real, numerator_imag, flow_real, flow_imag, flow_norm

  double mass = particles[last_particle_idx].mass;

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

  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++)
  {
//.........这里部分代码省略.........

//*********************************************************************************************
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;
//.........这里部分代码省略.........

void EmissionFunctionArray::calculate_dN_ptdptdphidy(int particle_idx)
// Calculate dN_xxx array.
{
  last_particle_idx = particle_idx;
  Stopwatch sw;
  sw.tic();

  double y = particle_y;
  particle_info* particle;
  particle = &particles[particle_idx];

  double mass = particle->mass;
  double sign = particle->sign;
  double degen = particle->gspin;

  double prefactor = 1.0/(8.0*(M_PI*M_PI*M_PI))/hbarC/hbarC/hbarC;

  FO_surf* surf = &FOsurf_ptr[0];

  // for intermediate results
  double dN_ptdptdphidy_tab[pT_tab_length][phi_tab_length];
  double dE_ptdptdphidy_tab[pT_tab_length][phi_tab_length];
  for (int i=0; i<pT_tab_length; i++)
  for (int j=0; j<phi_tab_length; j++)
  {
    dN_ptdptdphidy_tab[i][j] = 0.0;
    dE_ptdptdphidy_tab[i][j] = 0.0;
  }

  // pre-calculated variables //!!!!!!
  double trig_phi_table[phi_tab_length][2]; // 2: 0,1-> cos,sin
  for (int j=0; j<phi_tab_length; j++)
  {
    double phi = phi_tab->get(1,j+1);
    trig_phi_table[j][0] = cos(phi);
    trig_phi_table[j][1] = sin(phi);
  }

  double hypertrig_etas_table[eta_tab_length][2]; // 2: 0,1-> cosh,sinh
  double delta_eta_tab[eta_tab_length]; // cache it
  for (int k=0; k<eta_tab_length; k++)
  {
    double eta_s = eta_tab->get(1,k+1);
    hypertrig_etas_table[k][0] = cosh(y-eta_s);
    hypertrig_etas_table[k][1] = sinh(y-eta_s);
    delta_eta_tab[k] = eta_tab->get(2,k+1);
  }

  //---------------------------
  // THE main summation loop
  //---------------------------
  double progress_total = pT_tab_length*phi_tab_length;
  if (AMOUNT_OF_OUTPUT>0) print_progressbar(-1);

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

      for (int j=0; j<phi_tab_length; j++)
      //for (int j=0; j<1; j++) // for debug
      {
          // double phi = phi_tab->get(1,j+1), phi_weight = phi_tab->get(2,j+1); // unused

          // old way
          //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;
//.........这里部分代码省略.........