Polarization_RelativeLuminosity/macro/polarization.cc

0001 //#include <uspin/SpinDBContentv1.h> // for ana.505 or later
0002 #include <uspin/SpinDBContent.h>  // for ana.504 or earlier
0003 #include <uspin/SpinDBOutput.h>
0004
0005 #include "functions.hh"
0006 #include "sPhenixStyle.C"
0007
0008 const int kBunch_num_max = 120;
0009
0010 // It writes, for example,
0011 //     <Pb> = 50%
0012 // on the canvas
0013 void WriteAveragePolarization( double value, string beam, double pos_x, double pos_y, int digits=2 )
0014 {
0015   TLatex* tex = new TLatex();
0016   tex->SetTextSize( 0.04 );
0017   stringstream ss;
0018   ss << "#LTP_{" << beam << "}#GT = " << setprecision( digits ) << value << "%";
0019
0020   tex->DrawLatexNDC( pos_x, pos_y, ss.str().c_str() );
0021   delete tex;
0022 }
0023
0024 void WriteAdditionalMessage( double pos_x, double pos_y )
0025 {
0026   TLatex* tex = new TLatex();
0027   tex->SetTextSize( 0.04 );
0028   stringstream ss;
0029   ss << "RHIC Polarimetry Group Preliminary Results";
0030   //ss << "Preliminary CNI Polarimetry Group Results";
0031
0032   tex->DrawLatexNDC( pos_x, pos_y, ss.str().c_str() );
0033   delete tex;
0034 }
0035
0036 long long GetPolarizedLuminosity(const long long scalers[kBunch_num_max], const int pattern[kBunch_num_max], const int polarity)
0037 {
0038   long long sum = 0;
0039   for (int i = 0; i < kBunch_num_max; i++)
0040   {
0041     if (pattern[i] == polarity)
0042     {
0043       sum += scalers[i];
0044     }
0045   }
0046
0047   return sum;
0048 }
0049
0050 bool GetGoldenRuns(string list, vector<int>& runs)
0051 {
0052   ifstream ifs(list.c_str());
0053   if (ifs.fail())
0054   {
0055     cerr << "Fail to open " << list << endl;
0056     return false;
0057   }
0058
0059   int temp;
0060   while (ifs >> temp)
0061     runs.push_back(temp);
0062
0063   return true;
0064 }
0065
0066 int polarization()
0067 {
0068   SetsPhenixStyle();
0069   string golden_run_list = "list/Full_ppGoldenRunList_Version2.list";
0070   vector<int> golden_runs;
0071   GetGoldenRuns(golden_run_list, golden_runs);
0072   cout << "- Golden run list: " << golden_run_list << endl;
0073   cout << "- " << golden_runs.size() << " golden runs taken" << endl;
0074
0075   string output_base = "results/polarization";
0076   string output_root = output_base + ".root";
0077
0078   TFile* tf = new TFile(output_root.c_str(), "RECREATE");
0079
0080   /////////////////////////////////////////////////////////////////////////////////////
0081   // variables
0082   /////////////////////////////////////////////////////////////////////////////////////
0083   int run = 150045;  // 0;
0084   bool is_bad_run = true;
0085
0086   ////////////////////////
0087   // for Blue ring
0088   double pol_b[kBunch_num_max] = {-9999.9};        // blue beam polarization for each bunch
0089   double pol_error_b[kBunch_num_max] = {-9999.9};  // Error of blue beam polarization for each bunch
0090   double pol_sys_error_b[kBunch_num_max] = {0};
0091   int pattern_b[kBunch_num_max] = {-9999};  // spin pattern of blue beam
0092   Long64_t lumi_b = 0, lumi_pos_b = 0, lumi_nega_b = 0;
0093   double lumi_ratio_b = 0;
0094
0095   ////////////////////////
0096   // for Yellow ring
0097   double pol_y[kBunch_num_max] = {-9999.9};        // yellow beam polarization for each bunch
0098   double pol_error_y[kBunch_num_max] = {-9999.9};  // Error of yellow beam polarization for each bunch
0099   double pol_sys_error_y[kBunch_num_max] = {-9999};
0100   int pattern_y[kBunch_num_max] = {-9999};  // spin pattern of yellow beam
0101   Long64_t lumi_y = 0, lumi_pos_y = 0, lumi_nega_y = 0;
0102   double lumi_ratio_y = 0;
0103
0104   ////////////////////////
0105   // for both beam
0106   Long64_t mbd_scalers_no_cut[kBunch_num_max] = {0};
0107   double mean_mbd_scalers_no_cut = 0, sigma_mbd_scalers_no_cut = 0;
0108
0109   /////////////////////////////////////////////////////////////////////////////////////
0110   // Database setup
0111   /////////////////////////////////////////////////////////////////////////////////////
0112   int spin_db_status = 0;
0113
0114   // Get the spin patterns from the spin DB
0115   //  0xffff is the qa_level from XingShiftCal //////
0116   unsigned int qa_level = 0;  // 0xffff; // <--- old version
0117   // qa_level = 0xffff;
0118
0119   SpinDBContent spin_db_content;
0120   SpinDBOutput spin_db_output = SpinDBOutput("phnxrc");
0121
0122   int run_min = *min_element(golden_runs.begin(), golden_runs.end());
0123   int run_max = *max_element(golden_runs.begin(), golden_runs.end());
0124
0125   cout << "- Retrieving spin infomation from run " << run_min << " to " << run_max << "...... ";
0126   // spin_db_output.StoreDBContent(run, run, qa_level);
0127   spin_db_output.StoreDBContent(run_min, run_max);  // the default QA Lv. is used, <--- use this
0128   cout << "done" << endl;
0129
0130   /////////////////////////////////////////////////////////////////////////////////////
0131   // Histograms setup
0132   /////////////////////////////////////////////////////////////////////////////////////
0133   TH1D* hist_pol_b = new TH1D("hist_pol_b", "Blue beam polarization;Polarization [%];Integrated luminosity", 100, 0, 100);
0134   HistSetting(hist_pol_b, GetSphenixColor());
0135
0136   TH1D* hist_weighted_pol_b = new TH1D("hist_weighted_pol_b", "Blue beam polarization with luminosity weight;Polarization [%];Integrated luminosity [a.u.]", 100, 0, 100);
0137   HistSetting(hist_weighted_pol_b, GetSphenixColor());
0138
0139   TH1D* hist_pol_y = new TH1D("hist_pol_y", "Yellow beam polarization;Polarization [%];Integrated luminosity", 100, 0, 100);
0140   HistSetting(hist_pol_y, kOrange + 1);
0141
0142   TH1D* hist_weighted_pol_y = new TH1D("hist_weighted_pol_y", "Yellow beam polarization with luminosity weight;Polarization [%];Integrated luminosity [a.u.]", 100, 0, 100);
0143   HistSetting(hist_weighted_pol_y, kOrange + 1);
0144
0145   Long64_t luminosity_sum = 0;
0146   for (auto& run : golden_runs)
0147   {
0148     auto status = spin_db_output.GetDBContentStore(spin_db_content, run);
0149     if (status == 0)
0150     {
0151       continue;
0152     }
0153
0154     is_bad_run = bool(spin_db_content.GetBadRunFlag());  // int 0 -> false, int 1 -> true
0155
0156     if (is_bad_run == true)
0157       continue;
0158
0159     for (int i = 0; i < kBunch_num_max; i++)
0160     {
0161       spin_db_content.GetPolarizationBlue(i, pol_b[i], pol_error_b[i], pol_sys_error_b[i]);
0162       spin_db_content.GetPolarizationYellow(i, pol_y[i], pol_error_y[i], pol_sys_error_y[i]);
0163       pattern_b[i] = spin_db_content.GetSpinPatternBlue(i);
0164       pattern_y[i] = spin_db_content.GetSpinPatternYellow(i);
0165
0166       long long mbd_scaler = spin_db_content.GetScalerMbdNoCut(i);
0167       mbd_scalers_no_cut[i] = mbd_scaler;
0168
0169       /*
0170       cout << setw(4) << i << " "
0171            << setw(3) << pattern_b[i] << "\t"
0172            << pol_b[i] << "\t"
0173            << setw(3) << pattern_y[i] << "\t"
0174            << pol_y[i] << "\t"
0175            << mbd_scalers_no_cut[i] << endl;
0176            */
0177     }
0178
0179     /////////////////////////////////////////////////////////////////////////////////////
0180     // luminosity calc for blue/yellow depending on polarization direction
0181     // blue beam
0182     lumi_pos_b = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_b, 1);
0183     lumi_nega_b = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_b, -1);
0184     if (lumi_nega_b != 0)
0185       lumi_ratio_b = 1. * lumi_pos_b / lumi_nega_b;
0186     else
0187       lumi_ratio_b = 0;
0188
0189     lumi_b = lumi_pos_b + lumi_nega_b;
0190
0191     // yellow beam
0192     lumi_pos_y = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_y, 1);
0193     lumi_nega_y = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_y, -1);
0194     if (lumi_nega_y != 0)
0195       lumi_ratio_y = 1. * lumi_pos_y / lumi_nega_y;
0196     else
0197       lumi_ratio_y = 0;
0198
0199     lumi_y = lumi_pos_y + lumi_nega_y;
0200
0201     luminosity_sum += lumi_b;
0202     hist_pol_b->Fill(pol_b[0]);
0203     hist_weighted_pol_b->Fill(pol_b[0], lumi_b);
0204     hist_pol_y->Fill(pol_y[0]);
0205     hist_weighted_pol_y->Fill(pol_y[0], lumi_y);
0206
0207   }
0208
0209   hist_weighted_pol_b->Scale(1.0 / luminosity_sum);
0210   hist_weighted_pol_y->Scale(1.0 / luminosity_sum);
0211   cout << hist_pol_b->GetEntries() << " good runs" << endl;
0212
0213   ////////////////////////////////////////////////////////////////////
0214   // Draw!
0215   // string output_pdf = output_base + ".pdf";
0216   TCanvas* c = new TCanvas("canvas", "title", 800, 800);
0217   c->SetFillStyle( 0 );
0218   gPad->SetFrameFillStyle( 0 );
0219
0220   hist_weighted_pol_b->GetXaxis()->SetRangeUser(35, 65);
0221   hist_weighted_pol_b->GetXaxis()->SetLabelOffset( 0.0125 );
0222   //hist_weighted_pol_b->GetXaxis()->SetLabelSize( 0.03 );
0223
0224   hist_weighted_pol_b->GetYaxis()->SetTitleOffset( 1.5 );
0225   hist_weighted_pol_b->GetYaxis()->SetRangeUser( 0, 0.3 );
0226
0227   //for( int mode=0; mode<4; mode++ ) // loop over modes: preliminary, internal, work in progress
0228   {
0229     int mode = 3;
0230     // polarization with luminosity weight
0231     hist_weighted_pol_b->Draw("HIST");
0232     hist_weighted_pol_y->Draw("HIST same");
0233
0234     string output_pdf = output_base;
0235     if( mode == 0 )
0236       output_pdf += "_internal.pdf";
0237     else if( mode == 1 )
0238       output_pdf += "_preliminary.pdf";
0239     else if( mode == 2 )
0240       output_pdf += "_wip.pdf";
0241     else if( mode == 3 )
0242       output_pdf += "_performance.pdf";
0243
0244     int digits = 2;
0245     WriteDate( 0.775, 0.955, 0.04, "9/19/2025" );
0246     WritesPhenix( mode );
0247     double dy = 0.05, y = 0.89 + dy;
0248
0249     y -= dy;
0250     WriteRunCondition( 0.2, y); // p+p 2024 sqrt{s} = 200 GeV
0251
0252     y -= dy;//*1.25;
0253     WriteAdditionalMessage( 0.2, y); // Preliminary RHIC Polarimetry Group Results
0254
0255     y -= dy/2/2;
0256     TLegend* leg = new TLegend( 0.2, y-0.1, 0.4, y );
0257     leg->SetTextSize( 0.04 );
0258     leg->SetBorderSize( 0 );
0259     leg->SetFillStyle( 0 );
0260
0261     stringstream ss_legend_b;
0262     ss_legend_b << "Blue beam,  "
0263         << "#LTPol.#GT = " << setprecision( digits ) << hist_weighted_pol_b->GetMean() << "%";
0264     leg->AddEntry( hist_weighted_pol_b, ss_legend_b.str().c_str(), "l" );
0265
0266     stringstream ss_legend_y;
0267     ss_legend_y << "Yellow beam, "
0268         << "#LTPol.#GT = " << setprecision( digits ) << hist_weighted_pol_y->GetMean() << "%";
0269     leg->AddEntry( hist_weighted_pol_y, ss_legend_y.str().c_str(), "l" );
0270
0271     leg->Draw();
0272     c->Print( output_pdf.c_str() );
0273   }
0274
0275   tf->WriteTObject(hist_pol_b, hist_pol_b->GetName());
0276   tf->WriteTObject(hist_weighted_pol_b, hist_weighted_pol_b->GetName());
0277   tf->WriteTObject(hist_pol_y, hist_pol_y->GetName());
0278   tf->WriteTObject(hist_weighted_pol_y, hist_weighted_pol_y->GetName());
0279   tf->Close();
0280   // tr->Print();
0281   return 0;
0282 }