Polarization_RelativeLuminosity/macro/relative_luminosity.cc

0001 #include "polarization.cc"
0002
0003 int relative_luminosity()
0004 {
0005   SetsPhenixStyle();
0006   string golden_run_list = "list/Full_ppGoldenRunList_Version2.list";
0007   vector<int> golden_runs;
0008   GetGoldenRuns(golden_run_list, golden_runs);
0009   cout << "- Golden run list: " << golden_run_list << endl;
0010   cout << "- " << golden_runs.size() << " golden runs taken" << endl;
0011
0012   string output_base = "results/relative_luminosity";
0013   string output_root = output_base + ".root";
0014
0015   TFile* tf = new TFile(output_root.c_str(), "RECREATE");
0016
0017   /////////////////////////////////////////////////////////////////////////////////////
0018   // variables
0019   /////////////////////////////////////////////////////////////////////////////////////
0020   int run = 0;
0021
0022   ////////////////////////
0023   // for Blue ring
0024   double pol_b[kBunch_num_max] = {-9999.9};        // blue beam polarization for each bunch
0025   double pol_error_b[kBunch_num_max] = {-9999.9};  // Error of blue beam polarization for each bunch
0026   double pol_sys_error_b[kBunch_num_max] = {0};
0027   int pattern_b[kBunch_num_max] = {-9999};  // spin pattern of blue beam
0028
0029   ////////////////////////
0030   // for Yellow ring
0031   double pol_y[kBunch_num_max] = {-9999.9};        // yellow beam polarization for each bunch
0032   double pol_error_y[kBunch_num_max] = {-9999.9};  // Error of yellow beam polarization for each bunch
0033   double pol_sys_error_y[kBunch_num_max] = {-9999};
0034   int pattern_y[kBunch_num_max] = {-9999};  // spin pattern of yellow beam
0035
0036   ////////////////////////
0037   // for both beam
0038   Long64_t mbd_scalers_no_cut[kBunch_num_max] = {0};
0039   double mean_mbd_scalers_no_cut = 0, sigma_mbd_scalers_no_cut = 0;
0040
0041   /////////////////////////////////////////////////////////////////////////////////////
0042   // Database setup
0043   /////////////////////////////////////////////////////////////////////////////////////
0044   // Get the spin patterns from the spin DB
0045   //  0xffff is the qa_level from XingShiftCal //////
0046   unsigned int qa_level = 0;  // 0xffff; // <--- old version
0047   // qa_level = 0xffff;
0048
0049   SpinDBContent spin_db_content;
0050   SpinDBOutput spin_db_output = SpinDBOutput("phnxrc");
0051
0052   int run_min = *min_element(golden_runs.begin(), golden_runs.end());
0053   int run_max = *max_element(golden_runs.begin(), golden_runs.end());
0054
0055   cout << "- Retrieving spin infomation from run " << run_min << " to " << run_max << "...... ";
0056   // spin_db_output.StoreDBContent(run, run, qa_level);
0057   spin_db_output.StoreDBContent(run_min, run_max);  // the default QA Lv. is used, <--- use this
0058   cout << "done" << endl;
0059
0060   /////////////////////////////////////////////////////////////////////////////////////
0061   // Histograms setup
0062   /////////////////////////////////////////////////////////////////////////////////////
0063   double min_x = 0.8, max_x = 1.2;
0064   int bin = 80;
0065   TH1D* hist_lumi_b = new TH1D("hist_lumi_b", "Blue beam relative luminosity;Relative luminosity for TSSA;Integrated luminosity", bin, min_x, max_x);
0066   HistSetting(hist_lumi_b, GetSphenixColor());
0067
0068   TH1D* hist_weighted_lumi_b = new TH1D("hist_weighted_lumi_b", "Blue beam relative luminosity with luminosity weight;Relative luminosity = #mathcal{L};Integrated luminosity [a.u.]", bin, min_x, max_x);
0069   HistSetting(hist_weighted_lumi_b, GetSphenixColor());
0070
0071   TH1D* hist_lumi_y = new TH1D("hist_lumi_y", "Yellow beam relative luminosity;Relative luminosity for TSSA;Integrated luminosity", bin, min_x, max_x);
0072   HistSetting(hist_lumi_y, kOrange + 1);
0073
0074   TH1D* hist_weighted_lumi_y = new TH1D("hist_weighted_lumi_y", "Yellow beam relative luminosity with luminosity weight;Relative luminosity for TSSA;Integrated luminosity [a.u.]", bin, min_x, max_x); // #equiv#font[12]{L}^{#uparrow}/#font[12]{L}^{#downarrow}
0075   HistSetting(hist_weighted_lumi_y, kOrange + 1);
0076
0077   Long64_t luminosity_sum = 0;
0078   for (auto& run : golden_runs)
0079   {
0080     auto status = spin_db_output.GetDBContentStore(spin_db_content, run);
0081     if (status == 0)
0082     {
0083       continue;
0084     }
0085
0086     bool is_bad_run = bool(spin_db_content.GetBadRunFlag());  // int 0 -> false, int 1 -> true
0087
0088     if (is_bad_run == true)
0089     {
0090       continue;
0091     }
0092
0093     for (int i = 0; i < kBunch_num_max; i++)
0094     {
0095       spin_db_content.GetPolarizationBlue(i, pol_b[i], pol_error_b[i], pol_sys_error_b[i]);
0096       spin_db_content.GetPolarizationYellow(i, pol_y[i], pol_error_y[i], pol_sys_error_y[i]);
0097       pattern_b[i] = spin_db_content.GetSpinPatternBlue(i);
0098       pattern_y[i] = spin_db_content.GetSpinPatternYellow(i);
0099
0100       long long mbd_scaler = spin_db_content.GetScalerMbdNoCut(i);
0101       mbd_scalers_no_cut[i] = mbd_scaler;
0102        /*
0103       cout << setw(4) << i << " "
0104            << setw(3) << pattern_b[i] << "\t"
0105            << pol_b[i] << "\t"
0106            << setw(3) << pattern_y[i] << "\t"
0107            << pol_y[i] << "\t"
0108            << mbd_scalers_no_cut[i] << endl;
0109            */
0110     }
0111
0112     /////////////////////////////////////////////////////////////////////////////////////
0113     // luminosity calc for blue/yellow depending on polarization direction
0114     // blue beam
0115     Long64_t lumi_pos_b = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_b, 1);
0116     Long64_t lumi_nega_b = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_b, -1);
0117     double lumi_ratio_b = 0;
0118     if (lumi_nega_b != 0)
0119       lumi_ratio_b = 1. * lumi_pos_b / lumi_nega_b;
0120
0121     auto lumi_b = lumi_pos_b + lumi_nega_b;
0122
0123     // yellow beam
0124     Long64_t lumi_pos_y = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_y, 1);
0125     Long64_t lumi_nega_y = GetPolarizedLuminosity(mbd_scalers_no_cut, pattern_y, -1);
0126     double lumi_ratio_y = 0;
0127     if (lumi_nega_y != 0)
0128       lumi_ratio_y = 1. * lumi_pos_y / lumi_nega_y;
0129
0130     auto lumi_y = lumi_pos_y + lumi_nega_y;
0131
0132     luminosity_sum += lumi_b;
0133
0134     hist_lumi_b->Fill( lumi_ratio_b );
0135     hist_weighted_lumi_b->Fill(lumi_ratio_b, lumi_b);
0136     hist_lumi_y->Fill( lumi_ratio_y );
0137     hist_weighted_lumi_y->Fill( lumi_ratio_y , lumi_y);
0138
0139   } // end of loop over runs
0140
0141   hist_weighted_lumi_b->Scale(1.0 / luminosity_sum);
0142   hist_weighted_lumi_y->Scale(1.0 / luminosity_sum);
0143   cout << hist_lumi_b->GetEntries() << " good runs" << endl;
0144
0145   ////////////////////////////////////////////////////////////////////
0146   // Draw!
0147   // string output_pdf = output_base + ".pdf";
0148   TCanvas* c = new TCanvas("canvas", "title", 800, 800);
0149   c->SetFillStyle( 0 );
0150   gPad->SetFrameFillStyle( 0 );
0151
0152   hist_weighted_lumi_y->GetXaxis()->SetNdivisions( 1005 );
0153   hist_weighted_lumi_y->GetXaxis()->SetRangeUser( 0.9, 1.1);
0154   hist_weighted_lumi_y->GetXaxis()->SetLabelOffset( 0.0125 );
0155   hist_weighted_lumi_y->GetXaxis()->SetLabelSize( 0.04 );
0156   hist_weighted_lumi_y->GetYaxis()->SetTitleOffset( 1.5 );
0157   hist_weighted_lumi_y->GetYaxis()->SetLabelSize( 0.04 );
0158   hist_weighted_lumi_y->GetYaxis()->SetRangeUser( 0.0, 0.2 );
0159
0160   //for( int mode=0; mode<4; mode++ ) // loop over modes: preliminary, internal, work in progress
0161   {
0162     int mode = 3;
0163     // polarization with luminosity weight
0164     hist_weighted_lumi_y->Draw("HIST");
0165     hist_weighted_lumi_b->Draw("HIST same");
0166
0167     string output_pdf = output_base;
0168     if( mode == 0 )
0169       output_pdf += "_internal.pdf";
0170     else if( mode == 1 )
0171       output_pdf += "_preliminary.pdf";
0172     else if( mode == 2 )
0173       output_pdf += "_wip.pdf";
0174     else if( mode == 3 )
0175       output_pdf += "_performance.pdf";
0176
0177     int digits = 2;
0178     WriteDate( 0.775, 0.955, 0.04, "9/19/2025" );
0179     WritesPhenix( mode );
0180     double dy = 0.05, y = 0.875 + dy;
0181
0182     y -= dy;
0183     WriteRunCondition( 0.2, y, false, 0.04, true );
0184
0185     y += dy*3/4;
0186     TLegend* leg = new TLegend( 0.675, y - 0.1, 0.875, y );
0187     leg->SetTextSize( 0.04 );
0188     leg->SetBorderSize( 0 );
0189     leg->SetFillStyle( 0 );
0190     leg->AddEntry( hist_weighted_lumi_b, "Blue beam", "l" );
0191     leg->AddEntry( hist_weighted_lumi_y, "Yellow beam", "l" );
0192
0193     leg->Draw();
0194     c->Print( output_pdf.c_str() );
0195   }
0196
0197   tf->WriteTObject(hist_lumi_b, hist_lumi_b->GetName());
0198   tf->WriteTObject(hist_weighted_lumi_b, hist_weighted_lumi_b->GetName());
0199   tf->WriteTObject(hist_lumi_y, hist_lumi_y->GetName());
0200   tf->WriteTObject(hist_weighted_lumi_y, hist_weighted_lumi_y->GetName());
0201   tf->Close();
0202   // tr->Print();
0203   return 0;
0204 }