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 estimateLuminosityByRun(const std::string& runListName = "runs.txt", const std::string& outputfilename = "luminosity.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_run[nBeams][nScalers];
0034   TGraph *graph_relative_luminosity_run[nBeams][nScalers];
0035   TGraph *graph_luminosity_up_run[nBeams][nScalers];
0036   TGraph *graph_luminosity_down_run[nBeams][nScalers];
0037   TGraph *graph_crossangle_run[nBeams][nScalers];
0038   int runnumber;
0039   double Rmbd;
0040   double Rzdc;
0041   TTree *tree_relative_luminosity[nBeams];
0042   for (int iB = 0; iB < nBeams; iB++)
0043   {
0044     std::stringstream treename; treename << "tree_relative_luminosity_" << beams[iB]; 
0045     tree_relative_luminosity[iB] = new TTree(treename.str().c_str(), "relative luminosity tree");
0046     tree_relative_luminosity[iB]->Branch("runnumber", &runnumber, "runnumber/I"); 
0047     tree_relative_luminosity[iB]->Branch("mbd", &Rmbd, "mbd/D");
0048     tree_relative_luminosity[iB]->Branch("zdc", &Rzdc, "zdc/D");
0049     for (int iScal = 0; iScal < nScalers; iScal++) 
0050     {
0051       std::stringstream hname; hname << "h_polarization_" << beams[iB] << "_" << scalers[iScal];
0052       std::stringstream htitle; htitle << ";" << beams[iB] << " (" << scalers[iScal] << ") " << " beam polarization [%]; luminosity weighted distribution";
0053       h_polarization[iB][iScal] = new TH1F(hname.str().c_str(), htitle.str().c_str(), 50, 30, 65);
0054       hname.str(""); hname << "h_relative_luminosity_" << beams[iB] << "_" << scalers[iScal];
0055       htitle.str(""); htitle << ";" << beams[iB] << " (" << scalers[iScal] << ") " << " beam relative luminosity; luminosity weighted distribution";
0056       h_relative_luminosity[iB][iScal] = new TH1F(hname.str().c_str(), htitle.str().c_str(), 50, 0.9, 1.1);
0057       hname.str(""); hname << "graph_polarization_run_" << beams[iB] << "_" << scalers[iScal];
0058       htitle.str(""); htitle << ";run number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam polarization [%]";
0059       graph_polarization_run[iB][iScal] = new TGraph();
0060       graph_polarization_run[iB][iScal]->SetName(hname.str().c_str());
0061       graph_polarization_run[iB][iScal]->SetTitle(htitle.str().c_str());
0062       hname.str(""); hname << "graph_relative_luminosity_run_" << beams[iB] << "_" << scalers[iScal];
0063       htitle.str(""); htitle << ";run number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam relative luminosity";
0064       graph_relative_luminosity_run[iB][iScal] = new TGraph();
0065       graph_relative_luminosity_run[iB][iScal]->SetName(hname.str().c_str());
0066       graph_relative_luminosity_run[iB][iScal]->SetTitle(htitle.str().c_str());
0067       hname.str(""); hname << "graph_luminosity_up_run_" << beams[iB] << "_" << scalers[iScal];
0068       htitle.str(""); htitle << ";run number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity up";
0069       graph_luminosity_up_run[iB][iScal] = new TGraph();
0070       graph_luminosity_up_run[iB][iScal]->SetName(hname.str().c_str());
0071       graph_luminosity_up_run[iB][iScal]->SetTitle(htitle.str().c_str());
0072       hname.str(""); hname << "graph_luminosity_down_run_" << beams[iB] << "_" << scalers[iScal];
0073       htitle.str(""); htitle << ";run number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity down";
0074       graph_luminosity_down_run[iB][iScal] = new TGraph();
0075       graph_luminosity_down_run[iB][iScal]->SetName(hname.str().c_str());
0076       graph_luminosity_down_run[iB][iScal]->SetTitle(htitle.str().c_str());
0077 
0078       hname.str(""); hname << "graph_crossangle_run_" << beams[iB] << "_" << scalers[iScal];
0079       htitle.str(""); htitle << ";run number;" << beams[iB] << " (" << scalers[iScal] << ") " << " beam luminosity down";
0080       graph_crossangle_run[iB][iScal] = new TGraph();
0081       graph_crossangle_run[iB][iScal]->SetName(hname.str().c_str());
0082       graph_crossangle_run[iB][iScal]->SetTitle(htitle.str().c_str());
0083     }
0084   }
0085   
0086   unsigned long normalization[nBeams][nScalers] = {0};
0087   long double average_polarization[nBeams][nScalers] = {0};
0088   long double average_relative_luminosity[nBeams][nScalers] = {0};
0089   
0090   //std::cout << "runnumber beam scaler LuminosityUp LuminosityDown RelativeLuminosity" << std::endl;
0091   int count = -1;
0092   while (std::getline(fin, line))
0093   {
0094     runnumber = std::stoi(line);
0095     std::cout << "run = " << runnumber << std::endl;
0096     count++;
0097     for (int iB = 0; iB < nBeams; iB++)
0098     {
0099       for (int iScal = 0; iScal < nScalers; iScal++)
0100       {
0101         unsigned long LuminosityUp = 0;
0102         unsigned long LuminosityDown = 0;
0103         double pol = 0;
0104         double polerr = 0;
0105         double crossangle;
0106         double crossangle_std;
0107         std::stringstream GL1Pscalers;
0108         unsigned long sumScalers = 0;
0109         estimateLuminosityRun(iB, iScal, GL1Pscalers, sumScalers, runnumber, LuminosityUp, LuminosityDown, pol, polerr, crossangle, crossangle_std);
0110 
0111         unsigned long Luminosity = LuminosityUp + LuminosityDown;
0112         /*double blueRelativeLuminosity = (bluepol * (double) blueUpLuminosity + (1 - bluepol) * 0.5 * ((double) blueUpLuminosity + (double) blueDownLuminosity)) /
0113           (bluepol * (double) blueDownLuminosity + (1 - bluepol) * 0.5 * ((double) blueUpLuminosity + (double) blueDownLuminosity));
0114         */
0115         long double RelativeLuminosity = (long double) LuminosityUp / (long double) LuminosityDown;
0116         if (iScal == 0) Rmbd = RelativeLuminosity;
0117         else if (iScal == 1) Rzdc = RelativeLuminosity;
0118         
0119         if (false) {
0120           average_polarization[iB][iScal] += Luminosity * pol;
0121           average_relative_luminosity[iB][iScal] += (long double) Luminosity * (long double) LuminosityDown / (long double) LuminosityUp;
0122         } else {
0123           average_polarization[iB][iScal] += Luminosity * pol;
0124           average_relative_luminosity[iB][iScal] += (long double) Luminosity * (long double) LuminosityUp / (long double) LuminosityDown;
0125         }
0126 
0127         //average_relative_luminosity_blue += (float) blueLuminosity * ((float) blueUpLuminosity / (float) blueDownLuminosity);
0128         normalization[iB][iScal] += Luminosity;
0129 
0130         //std::cout << runnumber << " " << beams[iB] << " " << scalers[iScal] << " " << Luminosity << std::endl;
0131         //std::cout << runnumber << " " << beams[iB] << " " << scalers[iScal] << " " << LuminosityUp << " " << LuminosityDown << " " << RelativeLuminosity << std::endl;
0132         
0133         // fill histograms
0134         h_polarization[iB][iScal]->Fill(pol * 100, Luminosity);
0135         graph_polarization_run[iB][iScal]->SetPoint(count, runnumber, pol);
0136         h_relative_luminosity[iB][iScal]->Fill(RelativeLuminosity, Luminosity);
0137         graph_relative_luminosity_run[iB][iScal]->SetPoint(count, runnumber, RelativeLuminosity);
0138         graph_luminosity_up_run[iB][iScal]->SetPoint(count, runnumber, LuminosityUp);
0139         graph_luminosity_down_run[iB][iScal]->SetPoint(count, runnumber, LuminosityDown);
0140         graph_crossangle_run[iB][iScal]->SetPoint(count, runnumber, crossangle);
0141         //graph_crossangle_run[iB][iScal]->SetPointError(count, 0, crossangle_std);
0142       }
0143       //std::cout << runnumber << " " << beams
0144       tree_relative_luminosity[iB]->Fill();
0145     }
0146   }
0147   for (int iB = 0; iB < nBeams; iB++)
0148   {
0149     for (int iScal = 0; iScal < nScalers; iScal++)
0150     {
0151       average_polarization[iB][iScal] /= (long double) normalization[iB][iScal];
0152       average_relative_luminosity[iB][iScal] /= (long double) normalization[iB][iScal];
0153       
0154       std::cout << std::fixed << std::setprecision(10)
0155                 << "Average values " << beams[iB] << " " << scalers[iScal] << ":\n"
0156                 << "Pol = " << average_polarization[iB][iScal] << "  "
0157                 << "RL = " << average_relative_luminosity[iB][iScal] << "\n";
0158       
0159       graph_polarization_run[iB][iScal]->Write();
0160       graph_relative_luminosity_run[iB][iScal]->Write();
0161       graph_luminosity_up_run[iB][iScal]->Write();
0162       graph_luminosity_down_run[iB][iScal]->Write();
0163       graph_crossangle_run[iB][iScal]->Write();
0164       fin.close();
0165     }
0166   }
0167   
0168   outfile->Write();
0169   outfile->Close();
0170   delete outfile;
0171 
0172   gSystem->Exit(0);
0173 
0174 }
0175 
0176 void estimateLuminosityRun(int iB, int iScal, std::stringstream& GL1Pscalers, unsigned long& sumScalers, unsigned long runnumber, 
0177                            unsigned long& LuminosityUp, unsigned long& LuminosityDown,
0178                            double& pol, double& polerr, double& crossangle, double& crossangle_std)
0179 {
0180   SpinDBOutput spin_out("phnxrc");
0181   SpinDBContent *spin_cont = new SpinDBContentv1();
0182   spin_out.StoreDBContent(runnumber, runnumber);
0183   spin_out.GetDBContentStore(spin_cont, runnumber);
0184 
0185   // General information
0186   //std::cout << "Run number: " << runnumber << std::endl;
0187 
0188   // Spin patterns
0189   int spinpat[120] = {0};
0190   if (iB == 0)
0191   {
0192     for (int i = 0; i < 120; ++i)
0193     {
0194       spinpat[i] = -1 * spin_cont->GetSpinPatternBlue(i);
0195     }
0196   }
0197   else if (iB == 1)
0198   {
0199     for (int i = 0; i < 120; ++i)
0200     {
0201       spinpat[i] = -1 * spin_cont->GetSpinPatternYellow(i);
0202     }
0203   }
0204 
0205   crossangle = spin_cont->GetCrossAngle();
0206   crossangle_std = spin_cont->GetCrossAngleStd();
0207 
0208   // Beam polarization
0209   //double bluepol, yellpol, bluepolerr, yellpolerr;
0210   if (iB == 0) spin_cont->GetPolarizationBlue(0, pol, polerr);
0211   else if (iB == 1) spin_cont->GetPolarizationYellow(0, pol, polerr);
0212   pol /= 100; polerr /= 100;
0213 
0214   // GL1p Scalers
0215   // I don't know which one is the most relevant
0216   unsigned long scaler[120] = {0};
0217   //GL1Pscalers << "\"";
0218   sumScalers = 0;
0219   if (iScal == 0)
0220   {
0221     for (int i = 0; i < 111; ++i)
0222     {
0223       scaler[i] = spin_cont->GetScalerMbdNoCut(i); sumScalers += scaler[i];
0224     }
0225   } else if (iScal == 1)
0226   {
0227     for (int i = 0; i < 111; ++i)
0228     {
0229       scaler[i] = spin_cont->GetScalerZdcNoCut(i); sumScalers += scaler[i];
0230     }
0231   } else if (iScal == 2)
0232   {
0233     std::cout << "This should never happen!\n";
0234     for (int i = 0; i < 111; ++i)
0235     {
0236       scaler[i] = spin_cont->GetScalerMbdVertexCut(i); sumScalers += scaler[i];
0237     }
0238   }
0239 
0240   // Estimate luminosity:
0241   LuminosityUp = 0;
0242   LuminosityDown = 0;
0243   for (int i = 0; i < 120; i++)
0244   {
0245     if (spinpat[i] == 1)
0246       LuminosityUp += scaler[i];
0247     else if (spinpat[i] == -1)
0248       LuminosityDown += scaler[i];
0249   }
0250 }
0251
[ 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!