Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-07 08:08:32 

0001 #include "AnNeutralMeson_nano.h"
0002 #include <fun4all/Fun4AllReturnCodes.h>
0003 
0004 // Spin DB
0005 #include <uspin/SpinDBContent.h>
0006 #include <uspin/SpinDBContentv1.h>
0007 #include <uspin/SpinDBOutput.h>
0008 
0009 #include <algorithm>
0010 #include <fstream>
0011 #include <iostream>
0012 #include <iomanip> // for setprecision
0013 #include <random>
0014 #include <sstream>
0015 
0016 AnNeutralMeson_nano::AnNeutralMeson_nano(const std::string &name,
0017                                          const std::string &inputlistname,
0018                                          const std::string &inputfiletemplate,
0019                                          const std::string &outputfiletemplate)
0020   : SubsysReco(name)
0021   , inlistname(inputlistname)
0022   , infiletemplate(inputfiletemplate)
0023   , treename("tree_diphoton_compact")
0024   , outfiletemplate(outputfiletemplate)
0025 {
0026 }
0027 
0028 AnNeutralMeson_nano::~AnNeutralMeson_nano()
0029 {
0030 }
0031 
0032 int AnNeutralMeson_nano::Init(PHCompositeNode *)
0033 {
0034   std::ifstream inlistfile;
0035   inlistfile.open(inlistname.c_str());
0036 
0037   std::string runLine;
0038   while (std::getline(inlistfile, runLine))
0039   {
0040     int runnumber = std::stoi(runLine);
0041     runList.push_back(runnumber);
0042   }
0043   nRuns = runList.size();
0044   
0045   return Fun4AllReturnCodes::EVENT_OK;
0046 }
0047 
0048 int AnNeutralMeson_nano::process_event(PHCompositeNode *)
0049 {
0050   std::stringstream infilename;
0051   std::stringstream outfilename;
0052   std::cout << "nRuns = " << nRuns << std::endl;
0053   for (int irun = 0; irun < nRuns; irun++)
0054   {
0055     std::cout << "run " << runList[irun] << std::endl;
0056     infilename.str("");
0057     infilename << infiletemplate << runList[irun] << ".root";
0058     outfilename.str("");
0059     outfilename << outfiletemplate << runList[irun] << ".root";
0060 
0061     // Book histograms
0062     BookHistos(outfilename.str());
0063 
0064     // Open TTree
0065     TFile *infile = TFile::Open(infilename.str().c_str(), "READ");
0066 
0067     if (infile == nullptr)
0068     {
0069       std::cerr << "Error. Could not open file " << infilename.str()
0070                 << std::endl;
0071       return Fun4AllReturnCodes::ABORTEVENT;
0072     }
0073 
0074     TTree *nanoDST = nullptr;
0075     infile->GetObject(treename.c_str(), nanoDST);
0076     if (nanoDST == nullptr)
0077     {
0078       std::cerr << "Error: tree " << treename << " was not found in the file "
0079                 << infilename.str() << std::endl;
0080       return Fun4AllReturnCodes::ABORTEVENT;
0081     }
0082 
0083     // Optimization for reading
0084     nanoDST->SetBranchStatus("*", 0);
0085     nanoDST->SetBranchStatus("diphoton_vertex_z", 1);
0086     nanoDST->SetBranchStatus("diphoton_bunchnumber", 1);
0087     nanoDST->SetBranchStatus("diphoton_mass", 1);
0088     nanoDST->SetBranchStatus("diphoton_eta", 1);
0089     nanoDST->SetBranchStatus("diphoton_phi", 1);
0090     nanoDST->SetBranchStatus("diphoton_pt", 1);
0091     nanoDST->SetBranchStatus("diphoton_xf", 1);
0092 
0093     nanoDST->SetBranchAddress(
0094         "diphoton_vertex_z", &diphoton_vertex_z);
0095     nanoDST->SetBranchAddress(
0096         "diphoton_bunchnumber", &diphoton_bunchnumber);
0097     nanoDST->SetBranchAddress(
0098         "diphoton_mass", &diphoton_mass);
0099     nanoDST->SetBranchAddress(
0100         "diphoton_eta", &diphoton_eta);
0101     nanoDST->SetBranchAddress(
0102         "diphoton_phi", &diphoton_phi);
0103     nanoDST->SetBranchAddress(
0104         "diphoton_pt", &diphoton_pt);
0105     nanoDST->SetBranchAddress(
0106         "diphoton_xf", &diphoton_xf);
0107 
0108     // Spin pattern from SpinDB
0109     // Get spin pattern from SpinDB
0110 
0111     // Use the default QA level
0112     SpinDBOutput spin_out(
0113         "phnxrc");
0114     SpinDBContent *spin_cont = new SpinDBContentv1();
0115     spin_out.StoreDBContent(runList[irun], runList[irun]);
0116     spin_out.GetDBContentStore(spin_cont, runList[irun]);
0117 
0118     crossingshift = spin_cont->GetCrossingShift();
0119 
0120     for (int i = 0; i < nBunches; i++)
0121     {
0122       beamspinpat[0][i] = -1 * spin_cont->GetSpinPatternYellow(i);
0123       beamspinpat[1][i] = -1 * spin_cont->GetSpinPatternBlue(i);
0124       // spin pattern at sPHENIX is -1*(CDEV pattern)
0125       // The spin pattern corresponds to the expected spin at IP12, before the
0126       // siberian snake
0127     }
0128 
0129     if (seednb != 0) shuffle_spin_pattern(irun);
0130 
0131     // Loop over all recorded events;
0132     Long64_t nentries = nanoDST->GetEntries();
0133     std::cout << "nentries = " << nentries << std::endl;
0134     for (Long64_t ientry = 0; ientry < nentries; ientry++)
0135     {
0136       nanoDST->GetEntry(ientry);
0137 
0138       if (diphoton_pt < pTCutMin || diphoton_pt > pTCutMax) continue;
0139 
0140       if (require_low_vtx_cut) {
0141         if (std::abs(diphoton_vertex_z) > 30) continue;
0142       }
0143       else if (require_high_vtx_cut) {
0144         if (std::abs(diphoton_vertex_z) <= 30) continue;
0145       }
0146 
0147       if (require_phenix_cut) {
0148         if (std::abs(diphoton_eta) > 0.35) continue;
0149       }
0150       
0151       // Select the right pt bin:
0152       int iPt = FindBinBinary(diphoton_pt, pTBins, nPtBins);
0153       if (!((iPt < nPtBins) && (iPt >= 0)))
0154         continue;
0155 
0156       // Select the right eta bin:
0157       int ietaBin = FindBinBinary(diphoton_eta, etaBins, nEtaBins);
0158       if (!((ietaBin < nEtaBins) && (ietaBin >= 0)))
0159         continue;
0160 
0161       // Select the right xf bin:
0162       int ixfBin = FindBinBinary(diphoton_xf, xfBins, nXfBins);
0163       if (!((ixfBin < nXfBins) && (ixfBin >= 0)))
0164         continue;
0165 
0166       // Store invariant mass distributions
0167       h_pair_mass->Fill(diphoton_mass);
0168       if (iPt >= 0 && iPt < nPtBins)
0169         h_pair_mass_pt[iPt]->Fill(diphoton_mass);
0170       if (ietaBin >= 0 && ietaBin < nEtaBins)
0171         h_pair_mass_eta[ietaBin]->Fill(diphoton_mass);
0172       if (ixfBin >= 0 && ixfBin < nXfBins)
0173         h_pair_mass_xf[ixfBin]->Fill(diphoton_mass);
0174 
0175       // Select particle and region index;
0176       int ival = FindBinBinary(diphoton_mass, band_limits,
0177                                nParticles * (nRegions + 1) * 2);
0178       
0179       // Ignore any diphoton which is not within the side band or the peak band
0180       if (ival < 0 || ival % 2 == 1)
0181       {
0182         continue;
0183       }
0184       int iP = ival / 2 / (nRegions + 1);            // particle index
0185       int iR = (ival / 2 % (nRegions + 1) + 1) % 2;  // region index
0186 
0187       // Store kinematic correlations
0188       if (iR == 0) {
0189         h_pair_meson_zvtx[iP]->Fill(diphoton_vertex_z);
0190         h_pair_meson_pt_eta[iP][iPt]->Fill(diphoton_eta);
0191         h_pair_meson_pt_xf[iP][iPt]->Fill(diphoton_xf);
0192         h_pair_meson_eta_pt[iP][ietaBin]->Fill(diphoton_pt);
0193         h_pair_meson_eta_xf[iP][ietaBin]->Fill(diphoton_xf);
0194         h_pair_meson_xf_pt[iP][ixfBin]->Fill(diphoton_pt);
0195         h_pair_meson_xf_eta[iP][ixfBin]->Fill(diphoton_eta);
0196       }
0197 
0198       // Store pT, eta and xF values in order to compute the average value of
0199       // each bin
0200       if (iP != -1) 
0201       {
0202         if (iPt >= 0 && iPt < nPtBins) 
0203         {
0204           h_average_pt[iP]->Fill(iPt, diphoton_pt);
0205           h_norm_pt[iP]->Fill(iPt, 1);
0206         }
0207         if (ietaBin >= 0 && ietaBin < nEtaBins)
0208         {
0209           h_average_eta[iP]->Fill(ietaBin, diphoton_eta);
0210           h_norm_eta[iP]->Fill(ietaBin, 1);
0211         }
0212         if (ixfBin >= 0 && ixfBin < nXfBins)
0213         {
0214           h_average_xf[iP]->Fill(ixfBin, diphoton_xf);
0215           h_norm_xf[iP]->Fill(ixfBin, 1);
0216         }
0217       }
0218       
0219       // Compute yield
0220       float phi_beam = 0;
0221       for (int iB = 0; iB < nBeams; iB++)
0222       {
0223         // Forward vs backward asymmetry
0224         int iDir = (diphoton_eta * beamDirection[iB] > 0 ? 0 : 1);
0225         int ietaBinRelative =
0226             (beamDirection[iB] > 0 ? ietaBin : nEtaBins - 1 - ietaBin);
0227         int ixfBinRelative =
0228             (beamDirection[iB] > 0 ? ixfBin : nXfBins - 1 - ixfBin);
0229 
0230         if (require_high_xf_cut) {
0231           if (ixfBinRelative < 6) continue; // i.e. last most forward pT bins 6 and 7 -> xF > 0.035
0232         }
0233 
0234         if (require_low_xf_cut) {
0235           if (ixfBinRelative >= 6) continue; // i.e. last most forward pT bins 6 and 7 -> xF > 0.035
0236         }
0237 
0238         // phi global to phi yellow/blue
0239         phi_beam = WrapAngle(diphoton_phi + phi_shift[iB]); // i.e. phi_beam between - M_PI and +M_PI
0240 
0241         // select spin value;
0242         int iBunch = (crossingshift + diphoton_bunchnumber) % nBunches;
0243         int iS = -1;
0244         if (beamspinpat[iB][iBunch] == 1)
0245           iS = 0;
0246         else if (beamspinpat[iB][iBunch] == -1)
0247           iS = 1;
0248         if (iS == -1)
0249           continue;
0250 
0251         // Determine phi index
0252         if ((phi_beam < phiMin) || (phi_beam > phiMax))
0253           continue;
0254 
0255         // Increment kinematic-dependent yields
0256         h_yield_pt[iB][iP][iR][iPt][iDir][iS]->Fill(phi_beam);
0257         h_yield_eta[iB][iP][iR][ietaBinRelative][iS]->Fill(phi_beam);
0258         h_yield_xf[iB][iP][iR][ixfBinRelative][iS]->Fill(phi_beam);
0259 
0260         if (store_bunch_yields) {
0261           // Different pT thresholds for the pi0 and eta trigger
0262           if ((trigger_mbd && iP == 0 && iPt <= 1) ||
0263               (trigger_mbd && iP == 1 && iPt <= 2) ||
0264               (trigger_photon && iP == 0 && iPt >= 2) ||
0265               (trigger_photon && iP == 1 && iPt >= 3))
0266           {
0267             int iPhi = FindBinDirect(phi_beam, -M_PI, M_PI, nPhiBins);
0268             array_yield_pt[iB][iP][iR][iPt][iDir][iS][iPhi][iBunch]++;
0269             array_yield_eta[iB][iP][iR][ietaBinRelative][iS][iPhi][iBunch]++;
0270             array_yield_xf[iB][iP][iR][ixfBinRelative][iS][iPhi][iBunch]++;
0271           }
0272         }
0273       }
0274     }
0275     infile->Close();
0276 
0277     // Store histogram yields
0278     outfile->cd();
0279     outfile->Write();
0280     outfile->Close();
0281     delete outfile;
0282     outfile = nullptr; // To prevent double deletion with the destructor
0283 
0284     if (store_bunch_yields)
0285     {
0286       outfilename.str("");
0287       outfilename << outbunchtemplate << runList[irun] << ".bin";
0288       SaveBunchYields(outfilename.str());
0289     }
0290   }
0291   return Fun4AllReturnCodes::EVENT_OK;
0292 }
0293 
0294 int AnNeutralMeson_nano::End(PHCompositeNode *)
0295 {
0296   return 0;
0297 }
0298 
0299 void AnNeutralMeson_nano::BookHistos(const std::string &outputfilename)
0300 {
0301   // Book histograms
0302   outfile = new TFile(outputfilename.c_str(), "RECREATE");
0303   outfile->cd();
0304 
0305   // diphoton invariant mass
0306   h_pair_mass = new TH1F(
0307       "h_pair_mass",
0308       ";M_{#gamma} [GeV];counts", 500, 0, 1);
0309 
0310   for (int iPt = 0; iPt < nPtBins; iPt++)
0311   {
0312     std::stringstream h_pair_mass_name;
0313     h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_pt_"
0314                      << iPt;
0315     std::stringstream h_pair_mass_title;
0316     h_pair_mass_title << std::fixed << std::setprecision(1) << "diphoton mass ["
0317                       << pTBins[iPt] << " < p_{T} [GeV/c] < "
0318                       << pTBins[iPt + 1]
0319                       << "];M_{#gamma#gamma} [GeV/c^{2}];counts";
0320     h_pair_mass_pt[iPt] = new TH1F(h_pair_mass_name.str().c_str(),
0321                                    h_pair_mass_title.str().c_str(), 500, 0, 1);
0322   }
0323 
0324   
0325 
0326   for (int ietaBin = 0; ietaBin < nEtaBins; ietaBin++)
0327   {
0328     std::stringstream h_pair_mass_name;
0329     h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_eta_"
0330                      << ietaBin;
0331     std::stringstream h_pair_mass_title;
0332     h_pair_mass_title << std::fixed << std::setprecision(2) << "diphoton mass ["
0333                       << etaBins[ietaBin] << " < #eta < "
0334                       << etaBins[ietaBin + 1] << "];M_{#gamma#gamma};counts";
0335     h_pair_mass_eta[ietaBin] =
0336         new TH1F(h_pair_mass_name.str().c_str(),
0337                  h_pair_mass_title.str().c_str(), 500, 0, 1);
0338   }
0339 
0340   for (int ixfBin = 0; ixfBin < nXfBins; ixfBin++)
0341   {
0342     std::stringstream h_pair_mass_name;
0343     h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_xf_"
0344                      << ixfBin;
0345     std::stringstream h_pair_mass_title;
0346     h_pair_mass_title << std::fixed << std::setprecision(2) << "diphoton mass ["
0347                       << xfBins[ixfBin] << " < x_{F} < " << xfBins[ixfBin + 1]
0348                       << "];M_{#gamma#gamma};counts";
0349     h_pair_mass_xf[ixfBin] =
0350         new TH1F(h_pair_mass_name.str().c_str(),
0351                  h_pair_mass_title.str().c_str(), 500, 0, 1);
0352   }
0353 
0354   // Histogram for the average bin value
0355   h_average_pt[0] = new TH1F(
0356       "h_average_pt_pi0",
0357       ";iPt; pT average [GeV/c]", nPtBins,
0358       0, nPtBins);
0359   h_average_eta[0] = new TH1F(
0360       "h_average_eta_pi0",
0361       ";ieta; #eta average", nEtaBins, 0,
0362       nEtaBins);
0363   h_average_xf[0] = new TH1F(
0364       "h_average_xf_pi0",
0365       ";ixf; pT average", nXfBins, 0, nXfBins);
0366   h_norm_pt[0] = new TH1F(
0367       "h_norm_pt_pi0",
0368       ";iPt; pT norm [GeV/c]", nPtBins,
0369       0, nPtBins);
0370   h_norm_eta[0] = new TH1F(
0371       "h_norm_eta_pi0",
0372       ";ieta; #eta norm", nEtaBins, 0,
0373       nEtaBins);
0374   h_norm_xf[0] = new TH1F(
0375       "h_norm_xf_pi0",
0376       ";ixf; pT norm", nXfBins, 0, nXfBins);
0377   h_average_pt[1] = new TH1F(
0378       "h_average_pt_eta",
0379       ";iPt; pT average [GeV/c]", nPtBins,
0380       0, nPtBins);
0381   h_average_eta[1] = new TH1F(
0382       "h_average_eta_eta",
0383       ";ieta; #eta average", nEtaBins, 0,
0384       nEtaBins);
0385   h_average_xf[1] = new TH1F(
0386       "h_average_xf_eta",
0387       ";ixf; pT average", nXfBins, 0, nXfBins);
0388   h_norm_pt[1] = new TH1F(
0389       "h_norm_pt_eta",
0390       ";iPt; pT norm [GeV/c]", nPtBins,
0391       0, nPtBins);
0392   h_norm_eta[1] = new TH1F(
0393       "h_norm_eta_eta",
0394       ";ieta; #eta norm", nEtaBins, 0,
0395       nEtaBins);
0396   h_norm_xf[1] = new TH1F(
0397       "h_norm_xf_eta",
0398       ";ixf; pT norm", nXfBins, 0, nXfBins);
0399 
0400   for (int iP = 0; iP < 2; iP++) {
0401     std::stringstream h_name;
0402     h_name << "h_pair_" << particle[iP] << "_zvtx";
0403     h_pair_meson_zvtx[iP] = new TH1F(
0404       h_name.str().c_str(),
0405       ";z_{vtx} [cm]; Counts / [2 cm]",
0406       200, -200, 200);
0407   }
0408   
0409   // Kinematic relation pT vs eta vs xF
0410   // diphoton eta pT-bin by pT-bin
0411   for (int iP = 0; iP < 2; iP++) {
0412     for (int iPt = 1; iPt <= 9; iPt++) {
0413       std::stringstream h_name;
0414       h_name << "h_pair_" << particle[iP] << "_eta_pt_" << iPt;
0415       h_pair_meson_pt_eta[iP][iPt-1] = new TH1F(
0416         h_name.str().c_str(),
0417         ";#eta; Counts / [20 mrad]",
0418         200, -2.0, 2.0);
0419     }
0420   }
0421 
0422   // Kinematic relation pT vs eta vs xF
0423   // diphoton xF pT-bin by pT-bin
0424   for (int iP = 0; iP < 2; iP++) {
0425     for (int iPt = 1; iPt <= 9; iPt++) {
0426       std::stringstream h_name;
0427       h_name << "h_pair_" << particle[iP] << "_xf_pt_" << iPt;
0428       h_pair_meson_pt_xf[iP][iPt-1] = new TH1F(
0429         h_name.str().c_str(),
0430         ";x_{F}; Counts / [2 mrad]",
0431         200, -0.2, 0.2);
0432     }
0433   }
0434 
0435   // Kinematic relation pT vs eta vs xF
0436   // diphoton pT eta-bin by eta-bin
0437   for (int iP = 0; iP < 2; iP++) {
0438     for (int iEta = 1; iEta <= 8; iEta++) {
0439       std::stringstream h_name;
0440       h_name << "h_pair_" << particle[iP] << "_pt_eta_" << iEta;
0441       h_pair_meson_eta_pt[iP][iEta-1] = new TH1F(
0442         h_name.str().c_str(),
0443         ";p_{T} [GeV]; Counts / [100 MeV]",
0444         200, 0.0, 20);
0445     }
0446   }
0447 
0448   // Kinematic relation pT vs eta vs xF
0449   // diphoton xF eta-bin by eta-bin
0450   for (int iP = 0; iP < 2; iP++) {
0451     for (int iEta = 1; iEta <= 8; iEta++) {
0452       std::stringstream h_name;
0453       h_name << "h_pair_" << particle[iP] << "_xf_eta_" << iEta;
0454       h_pair_meson_eta_xf[iP][iEta-1] = new TH1F(
0455         h_name.str().c_str(),
0456         ";x_{F}; Counts / [2 mrad]",
0457         200, -0.2, 0.2);
0458     }
0459   }
0460 
0461   // Kinematic relation pT vs eta vs xF
0462   // diphoton pT xf-bin by xf-bin
0463   for (int iP = 0; iP < 2; iP++) {
0464     for (int iXf = 1; iXf <= 8; iXf++) {
0465       std::stringstream h_name;
0466       h_name << "h_pair_" << particle[iP] << "_pt_xf_" << iXf;
0467       h_pair_meson_xf_pt[iP][iXf-1] = new TH1F(
0468         h_name.str().c_str(),
0469         ";p_{T} [GeV]; Counts / [100 MeV]",
0470         200, 0.0, 20);
0471     }
0472   }
0473 
0474   // Kinematic relation pT vs eta vs xF
0475   // diphoton eta xf-bin by xf-bin
0476   for (int iP = 0; iP < 2; iP++) {
0477     for (int iXf = 1; iXf <= 8; iXf++) {
0478       std::stringstream h_name;
0479       h_name << "h_pair_" << particle[iP] << "_eta_xf_" << iXf;
0480       h_pair_meson_xf_eta[iP][iXf-1] = new TH1F(
0481         h_name.str().c_str(),
0482         ";#eta; Counts / [20 mrad]",
0483         200, -2.0, 2.0);
0484     }
0485   }
0486   
0487   // Beam- spin- and kinematic-dependent yields -> pT dependent
0488   for (int iB = 0; iB < nBeams; iB++)
0489   {
0490     for (int iP = 0; iP < nParticles; iP++)
0491     {
0492       for (int iR = 0; iR < nRegions; iR++)
0493       {
0494         for (int iS = 0; iS < nSpins; iS++)
0495         {
0496           for (int iPt = 0; iPt < nPtBins; iPt++)
0497           {
0498             for (int iDir = 0; iDir < nDirections; iDir++)
0499             {
0500               std::stringstream h_yield_name;
0501               h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0502                            << "_" << regions[iR] << "_"
0503                            << "pT_" << iPt << "_"
0504                            << directions[iDir]
0505                            << "_" << spins[iS];
0506               std::stringstream h_yield_title;
0507               h_yield_title << std::fixed << std::setprecision(1)
0508                             << "P_{T} #in [" << pTBins[iPt] << ", "
0509                             << pTBins[iPt + 1] << "] " << directions[iDir]
0510                             << " " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0511                             << regions[iR] << " range);"
0512                             << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0513                             << "};counts";
0514               h_yield_pt[iB][iP][iR][iPt][iDir][iS] =
0515                 new TH1F(h_yield_name.str().c_str(),
0516                          h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0517             }
0518           }
0519           for (int ietaBin = 0; ietaBin < nEtaBins; ietaBin++)
0520           {
0521             std::stringstream h_yield_name;
0522             h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0523                          << "_" << regions[iR] << "_"
0524                          << "eta_" << (int) ietaBin << "_" << spins[iS];
0525             std::stringstream h_yield_title;
0526             h_yield_title << std::fixed << std::setprecision(1) << "#eta #in ["
0527                           << etaBins[ietaBin] << ", " << etaBins[ietaBin]
0528                           << "] " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0529                           << regions[iR] << " range);"
0530                           << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0531                           << "};counts";
0532             h_yield_eta[iB][iP][iR][ietaBin][iS] =
0533                 new TH1F(h_yield_name.str().c_str(),
0534                          h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0535           }
0536           for (int ixfBin = 0; ixfBin < nXfBins; ixfBin++)
0537           {
0538             std::stringstream h_yield_name;
0539             h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0540                          << "_" << regions[iR] << "_"
0541                          << "xf_" << (int) ixfBin << "_" << spins[iS];
0542             std::stringstream h_yield_title;
0543             h_yield_title << std::fixed << std::setprecision(1) << "x_{F} #in ["
0544                           << xfBins[ixfBin] << ", " << xfBins[ixfBin + 1]
0545                           << "] " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0546                           << regions[iR] << " range);"
0547                           << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0548                           << "};counts";
0549             h_yield_xf[iB][iP][iR][ixfBin][iS] =
0550                 new TH1F(h_yield_name.str().c_str(),
0551                          h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0552           }
0553         }
0554       }
0555     }
0556   }
0557 }
0558 
0559 void AnNeutralMeson_nano::shuffle_spin_pattern(const int irun)
0560 {
0561   // Shuffle the spin pattern if the seed number is non zero
0562   // Only consider the first 111/(nBunches = 120) bunches, the last 9 are
0563   // slots are empty.
0564   const int nFilledBunches = 111;
0565   int shuffled_indices[nBunches] = {0};
0566   std::iota(shuffled_indices, shuffled_indices + nFilledBunches, 0);
0567 
0568   std::mt19937 rng(seednb + 100000 * irun);
0569   // Pseudo-random generator: Mersenne twister with a period of 2^19937 - 1
0570   std::shuffle(shuffled_indices, shuffled_indices + nFilledBunches, rng);
0571 
0572   // With a probability 1/2, flip all the bunch spins
0573   std::uniform_int_distribution<std::mt19937::result_type> dist2(0, 1);
0574   int factor = 2 * (int) dist2(rng) - 1;
0575   int beamspinpat_tmp[2][nBunches] = {0};
0576   for (int ibunch = 0; ibunch < nBunches; ibunch++)
0577   {
0578     int ishuffle = shuffled_indices[ibunch];
0579     beamspinpat_tmp[0][ibunch] = factor * beamspinpat[0][ishuffle];
0580     beamspinpat_tmp[1][ibunch] = factor * beamspinpat[1][ishuffle];
0581   }
0582   for (int ibunch = 0; ibunch < nBunches; ibunch++)
0583   {
0584     beamspinpat[0][ibunch] = beamspinpat_tmp[0][ibunch];
0585     beamspinpat[1][ibunch] = beamspinpat_tmp[1][ibunch];
0586   }
0587 }
0588 
0589 float AnNeutralMeson_nano::WrapAngle(const float phi)
0590 {
0591   float phi_out = phi;
0592   while (phi_out < -M_PI)
0593   {
0594     phi_out += 2 * M_PI;
0595   }
0596   while (phi_out > M_PI)
0597   {
0598     phi_out -= 2 * M_PI;
0599   }
0600   return phi_out;
0601 }
0602 
0603 void AnNeutralMeson_nano::SaveBunchYields(const std::string &outputfilename)
0604 {
0605   // Binary output
0606   std::ofstream outbinary(outputfilename, std::ios::binary);
0607   outbinary.write(reinterpret_cast<const char *>(array_yield_pt),
0608                   sizeof(array_yield_pt));
0609   outbinary.write(reinterpret_cast<const char *>(array_yield_eta),
0610                   sizeof(array_yield_eta));
0611   outbinary.write(reinterpret_cast<const char *>(array_yield_xf),
0612                   sizeof(array_yield_xf));
0613   outbinary.close();
0614 
0615   // Reset arrays to zero
0616   std::memset(array_yield_pt, 0, sizeof(array_yield_pt));
0617   std::memset(array_yield_eta, 0, sizeof(array_yield_eta));
0618   std::memset(array_yield_xf, 0, sizeof(array_yield_xf));
0619 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!