Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-04 08:08:21 

0001 #include <uspin/SpinDBContent.h>
0002 #include <uspin/SpinDBOutput.h>
0003 
0004 R__LOAD_LIBRARY(libuspin.so)
0005 
0006 void estimateLuminosityRun(int iB, int iScal, std::stringstream& GL1Pscalers, unsigned long& sumScalers, unsigned long runnumber, 
0007                             unsigned long& LuminosityUp, unsigned long& LuminosityDown,
0008                             double& pol, double& polerr, double& crossangle, double& crossangle_std);
0009 
0010 void estimateLuminosityByFill(const std::string& runListName = "run_to_fill.txt", const std::string& outputfilename = "luminosity_fill.root")
0011 {
0012   std::ifstream fin;
0013   fin.open(runListName);
0014   if (!(fin.is_open()))
0015   {
0016     std::cerr << "Error: File " << runListName << " could not be opened.\n";
0017   }
0018   std::string line;
0019 
0020   // Polarizations and luminosities are defined differently for each beam and scaler
0021   const int nBeams = 2;
0022   const std::string beams[2] = {"blue", "yellow"};
0023   const int nScalers = 2;
0024   const std::string scalers[2] = {"mbdns", "zdcns"};
0025 
0026   // Book histograms
0027   TFile *outfile = new TFile(outputfilename.c_str(), "RECREATE");
0028   outfile->cd();
0029 
0030   // Polarization graphs
0031   TH1F *h_polarization[nBeams][nScalers];
0032   TH1F *h_relative_luminosity[nBeams][nScalers];
0033   TGraph *graph_polarization_fill[nBeams][nScalers];
0034   TGraph *graph_relative_luminosity_fill[nBeams][nScalers];
0035   TGraph *graph_luminosity_up_fill[nBeams][nScalers];
0036   TGraph *graph_luminosity_down_fill[nBeams][nScalers];
0037   TGraph *graph_crossangle_fill[nBeams][nScalers];
0038   int runnumber;
0039   int fillnumber;
0040   double Rmbd;
0041   double Rzdc;
0042   TTree *tree_relative_luminosity[nBeams];
0043   for (int iB = 0; iB < nBeams; iB++)
0044   {
0045     std::stringstream treename; treename << "tree_relative_luminosity_" << beams[iB]; 
0046     tree_relative_luminosity[iB] = new TTree(treename.str().c_str(), "relative luminosity tree");
0047     tree_relative_luminosity[iB]->Branch("fillnumber", &fillnumber, "fillnumber/I"); 
0048     tree_relative_luminosity[iB]->Branch("mbd", &Rmbd, "mbd/D");
0049     tree_relative_luminosity[iB]->Branch("zdc", &Rzdc, "zdc/D");
0050     for (int iScal = 0; iScal < nScalers; iScal++) 
0051     {
0052       std::stringstream hname; hname << "h_polarization_" << beams[iB] << "_" << scalers[iScal];
0053       std::stringstream htitle; htitle << ";" << beams[iB] << " (" << scalers[iScal] << ") " << " beam polarization [%]; luminosity weighted distribution";
0054       h_polarization[iB][iScal] = new TH1F(hname.str().c_str(), htitle.str().c_str(), 50, 30, 65);
0055       hname.str(""); hname << "h_relative_luminosity_" << beams[iB] << "_" << scalers[iScal];
0056       htitle.str(""); htitle << ";" << beams[iB] << " (" << scalers[iScal] << ") " << " beam relative luminosity; luminosity weighted distribution";
0057       h_relative_luminosity[iB][iScal] = new TH1F(hname.str().c_str(), htitle.str().c_str(), 50, 0.9, 1.1);
0058       hname.str(""); hname << "graph_polarization_fill_" << beams[iB] << "_" << scalers[iScal];
0059       htitle.str(""); htitle << ";fill number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam polarization [%]";
0060       graph_polarization_fill[iB][iScal] = new TGraph();
0061       graph_polarization_fill[iB][iScal]->SetName(hname.str().c_str());
0062       graph_polarization_fill[iB][iScal]->SetTitle(htitle.str().c_str());
0063       hname.str(""); hname << "graph_relative_luminosity_fill_" << beams[iB] << "_" << scalers[iScal];
0064       htitle.str(""); htitle << ";fill number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam relative luminosity";
0065       graph_relative_luminosity_fill[iB][iScal] = new TGraph();
0066       graph_relative_luminosity_fill[iB][iScal]->SetName(hname.str().c_str());
0067       graph_relative_luminosity_fill[iB][iScal]->SetTitle(htitle.str().c_str());
0068       hname.str(""); hname << "graph_luminosity_up_fill_" << beams[iB] << "_" << scalers[iScal];
0069       htitle.str(""); htitle << ";fill number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity up";
0070       graph_luminosity_up_fill[iB][iScal] = new TGraph();
0071       graph_luminosity_up_fill[iB][iScal]->SetName(hname.str().c_str());
0072       graph_luminosity_up_fill[iB][iScal]->SetTitle(htitle.str().c_str());
0073       hname.str(""); hname << "graph_luminosity_down_fill_" << beams[iB] << "_" << scalers[iScal];
0074       htitle.str(""); htitle << ";fill number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity down";
0075       graph_luminosity_down_fill[iB][iScal] = new TGraph();
0076       graph_luminosity_down_fill[iB][iScal]->SetName(hname.str().c_str());
0077       graph_luminosity_down_fill[iB][iScal]->SetTitle(htitle.str().c_str());
0078 
0079       hname.str(""); hname << "graph_crossangle_fill_" << beams[iB] << "_" << scalers[iScal];
0080       htitle.str(""); htitle << ";fill number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity down";
0081       graph_crossangle_fill[iB][iScal] = new TGraph();
0082       graph_crossangle_fill[iB][iScal]->SetName(hname.str().c_str());
0083       graph_crossangle_fill[iB][iScal]->SetTitle(htitle.str().c_str());
0084     }
0085   }
0086   
0087   unsigned long normalization[nBeams][nScalers] = {0};
0088   long double average_polarization[nBeams][nScalers] = {0};
0089   long double total_luminosity_up[nBeams][nScalers] = {0};
0090   long double total_luminosity_down[nBeams][nScalers] = {0};
0091   long double average_relative_luminosity[nBeams][nScalers] = {0};
0092   
0093   //std::cout << "runnumber beam scaler LuminosityUp LuminosityDown RelativeLuminosity" << std::endl;
0094   std::string entry_content;
0095   // Read all fill and run numbers
0096   std::map<int, std::vector<int>> fill2runs;
0097   while (std::getline(fin, line)) // Run number loop
0098   {
0099     std::istringstream entry(line);
0100     std::getline(entry, entry_content, ' ');
0101     fillnumber = std::stoi(entry_content);
0102     std::getline(entry, entry_content, ' ');
0103     runnumber = std::stoi(entry_content);
0104     fill2runs[fillnumber].push_back(runnumber);
0105   }
0106   
0107   for (int iB = 0; iB < nBeams; iB++)
0108   {
0109     for (int iScal = 0; iScal < nScalers; iScal++)
0110     {
0111       int count = -1;
0112       for (const auto& [fill, runs] : fill2runs)
0113       {
0114         count++;
0115         fillnumber = fill;
0116         std::cout << "fill number = " << fillnumber << std::endl;
0117         unsigned long LuminosityUp = 0;
0118         unsigned long LuminosityDown = 0;
0119         double pol = 0;
0120         double polerr = 0;
0121         double crossangle;
0122         double crossangle_std;
0123         std::stringstream GL1Pscalers;
0124         unsigned long sumScalers = 0;
0125         for (int irun = 0; irun < (int) runs.size(); irun++)
0126         {
0127           runnumber = runs[irun];
0128           estimateLuminosityRun(iB, iScal, GL1Pscalers, sumScalers,
0129                                 runnumber, LuminosityUp, LuminosityDown,
0130                                 pol, polerr, crossangle, crossangle_std);
0131         }
0132         
0133         unsigned long Luminosity = LuminosityUp + LuminosityDown;
0134         long double RelativeLuminosity = (long double) LuminosityUp / (long double) LuminosityDown;
0135         if (iScal == 0) Rmbd = RelativeLuminosity;
0136         else if (iScal == 1) Rzdc = RelativeLuminosity;
0137         
0138 
0139         average_polarization[iB][iScal] += Luminosity * pol;
0140         total_luminosity_up[iB][iScal] += (long double) LuminosityUp;
0141         total_luminosity_down[iB][iScal] += (long double) LuminosityDown;
0142 
0143         //average_relative_luminosity_blue += (float) blueLuminosity * ((float) blueUpLuminosity / (float) blueDownLuminosity);
0144         normalization[iB][iScal] += Luminosity;
0145 
0146         //std::cout << runnumber << " " << beams[iB] << " " << scalers[iScal] << " " << Luminosity << std::endl;
0147         //std::cout << runnumber << " " << beams[iB] << " " << scalers[iScal] << " " << LuminosityUp << " " << LuminosityDown << " " << RelativeLuminosity << std::endl;
0148         
0149         // fill histograms
0150         h_polarization[iB][iScal]->Fill(pol * 100, Luminosity);
0151         graph_polarization_fill[iB][iScal]->SetPoint(count, fillnumber, pol);
0152         h_relative_luminosity[iB][iScal]->Fill(RelativeLuminosity, Luminosity);
0153         graph_relative_luminosity_fill[iB][iScal]->SetPoint(count, fillnumber, RelativeLuminosity);
0154         graph_luminosity_up_fill[iB][iScal]->SetPoint(count, fillnumber, LuminosityUp);
0155         graph_luminosity_down_fill[iB][iScal]->SetPoint(count, fillnumber, LuminosityDown);
0156         graph_crossangle_fill[iB][iScal]->SetPoint(count, fillnumber, crossangle);
0157         //graph_crossangle_run[iB][iScal]->SetPointError(count, 0, crossangle_std);
0158       }
0159       //std::cout << runnumber << " " << beams
0160       tree_relative_luminosity[iB]->Fill();
0161     }
0162   }
0163   for (int iB = 0; iB < nBeams; iB++)
0164   {
0165     for (int iScal = 0; iScal < nScalers; iScal++)
0166     {
0167       average_polarization[iB][iScal] /= (long double) normalization[iB][iScal];
0168       average_relative_luminosity[iB][iScal] = total_luminosity_up[iB][iScal] / total_luminosity_down[iB][iScal];
0169       
0170       std::cout << std::fixed << std::setprecision(10)
0171                 << "Average values " << beams[iB] << " " << scalers[iScal] << ":\n"
0172                 << "Pol = " << average_polarization[iB][iScal] << "  "
0173                 << "RL = " << average_relative_luminosity[iB][iScal] << "\n";
0174       
0175       graph_polarization_fill[iB][iScal]->Write();
0176       graph_relative_luminosity_fill[iB][iScal]->Write();
0177       graph_luminosity_up_fill[iB][iScal]->Write();
0178       graph_luminosity_down_fill[iB][iScal]->Write();
0179       graph_crossangle_fill[iB][iScal]->Write();
0180       fin.close();
0181     }
0182   }
0183   
0184   outfile->Write();
0185   outfile->Close();
0186   delete outfile;
0187 
0188   gSystem->Exit(0);
0189 
0190 }
0191 
0192 void estimateLuminosityRun(int iB, int iScal, std::stringstream& GL1Pscalers, unsigned long& sumScalers, unsigned long runnumber, 
0193                               unsigned long& LuminosityUp, unsigned long& LuminosityDown,
0194                               double& pol, double& polerr, double& crossangle, double& crossangle_std)
0195 {
0196   SpinDBOutput spin_out("phnxrc");
0197   SpinDBContent *spin_cont = new SpinDBContentv1();
0198 
0199   spin_out.StoreDBContent(runnumber, runnumber);
0200   spin_out.GetDBContentStore(spin_cont, runnumber);
0201 
0202   // General information
0203   //std::cout << "Run number: " << runnumber << std::endl;
0204 
0205   // Spin patterns
0206   int spinpat[120] = {0};
0207   if (iB == 0)
0208   {
0209     for (int i = 0; i < 120; ++i)
0210     {
0211       spinpat[i] = -1 * spin_cont->GetSpinPatternBlue(i);
0212     }
0213   }
0214   else if (iB == 1)
0215   {
0216     for (int i = 0; i < 120; ++i)
0217     {
0218       spinpat[i] = -1 * spin_cont->GetSpinPatternYellow(i);
0219     }
0220   }
0221 
0222   crossangle = spin_cont->GetCrossAngle();
0223   crossangle_std = spin_cont->GetCrossAngleStd();
0224 
0225   // Beam polarization
0226   //double bluepol, yellpol, bluepolerr, yellpolerr;
0227   if (iB == 0) spin_cont->GetPolarizationBlue(0, pol, polerr);
0228   else if (iB == 1) spin_cont->GetPolarizationYellow(0, pol, polerr);
0229   pol /= 100; polerr /= 100;
0230 
0231   // GL1p Scalers
0232   // I don't know which one is the most relevant
0233   unsigned long scaler[120] = {0};
0234   //GL1Pscalers << "\"";
0235   sumScalers = 0;
0236   if (iScal == 0)
0237   {
0238     for (int i = 0; i < 120; ++i)
0239     {
0240       scaler[i] = spin_cont->GetScalerMbdNoCut(i); sumScalers += scaler[i];
0241     }
0242   } else if (iScal == 1)
0243   {
0244     for (int i = 0; i < 120; ++i)
0245     {
0246       scaler[i] = spin_cont->GetScalerZdcNoCut(i); sumScalers += scaler[i];
0247     }
0248   } else if (iScal == 2)
0249   {
0250     std::cout << "This should never happen!\n";
0251     for (int i = 0; i < 120; ++i)
0252     {
0253       scaler[i] = spin_cont->GetScalerMbdVertexCut(i); sumScalers += scaler[i];
0254     }
0255   }
0256 
0257   // Estimate luminosity:
0258   // LuminosityUp = 0;
0259   // LuminosityDown = 0;
0260   for (int i = 0; i < 120; i++)
0261   {
0262     if (spinpat[i] == 1)
0263       LuminosityUp += scaler[i];
0264     else if (spinpat[i] == -1)
0265       LuminosityDown += scaler[i];
0266   }
0267 }
0268
[ 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!