bunchAym/macro/DrawBunchAysm_UD.C

0001 #include <TFile.h>
0002 #include <TH1I.h>
0003 #include <TGraphErrors.h>
0004 #include <TCanvas.h>
0005 #include <TMath.h>
0006 #include <iostream>
0007 #include <string>
0008
0009 #include <uspin/SpinDBContent.h>
0010 #include <uspin/SpinDBOutput.h>
0011
0012 R__LOAD_LIBRARY(libuspin.so)
0013 static const int NBUNCHES = 120;
0014 float zdc_adc_cut_[2][3] = {0.}; // ADC cut for each ZDC module ::array[side][module]
0015 float smd_adc_cut_[2][2] = {0.}; // ADC cut for each smd  ::array[side][direction]
0016 int spinPatternBlue[NBUNCHES] = {0};
0017 int spinPatternYellow[NBUNCHES] = {0};
0018 int crossingShift = 1;
0019
0020 void GetSpinPatterns(int runNum)
0021 {
0022   // Get the spin patterns from the spin DB
0023
0024   //  0xffff is the qa_level from XingShiftCal //////
0025   unsigned int qa_level = 0xffff;
0026   SpinDBContent spin_cont;
0027   SpinDBOutput spin_out("phnxrc");
0028
0029   spin_out.StoreDBContent(runNum, runNum, qa_level);
0030   spin_out.GetDBContentStore(spin_cont, runNum);
0031
0032   // Get crossing shift
0033   crossingShift = spin_cont.GetCrossingShift();
0034   std::cout << "Crossing shift: " << crossingShift << std::endl;
0035
0036   std::cout << "Blue spin pattern: [";
0037   for (int i = 0; i < 120; i++)
0038   {
0039     spinPatternBlue[i] = spin_cont.GetSpinPatternBlue(i);
0040     std::cout << spinPatternBlue[i];
0041     if (i < 119)
0042       std::cout << ", ";
0043   }
0044   std::cout << "]" << std::endl;
0045
0046   std::cout << "Yellow spin pattern: [";
0047   for (int i = 0; i < 120; i++)
0048   {
0049     spinPatternYellow[i] = spin_cont.GetSpinPatternYellow(i);
0050     std::cout << spinPatternYellow[i];
0051     if (i < 119)
0052       std::cout << ", ";
0053   }
0054   std::cout << "]" << std::endl;
0055
0056   if (false)
0057   {
0058     // for testing -- if we couldn't find the spin pattern, fill it with a dummy pattern
0059     // +-+-+-+- ...
0060     std::cout << "Could not find spin pattern packet for blue beam! Using dummy pattern" << std::endl;
0061     for (int i = 0; i < NBUNCHES; i++)
0062     {
0063       int mod = i % 2;
0064       if (mod == 0)
0065         spinPatternBlue[i] = 1;
0066       else
0067         spinPatternBlue[i] = -1;
0068     }
0069     // for testing -- if we couldn't find the spin pattern, fill it with a dummy pattern
0070     // ++--++-- ,,,
0071     std::cout << "Could not find spin pattern packet for yellow beam! Using dummy pattern" << std::endl;
0072     for (int i = 0; i < NBUNCHES; i++)
0073     {
0074       int mod = i % 4;
0075       if (mod == 0 || mod == 1)
0076         spinPatternYellow[i] = 1;
0077       else
0078         spinPatternYellow[i] = -1;
0079     }
0080   }
0081 }
0082 void DrawBunchAysm_UD(int runnumber = 42797)
0083 {
0084 //    std::string filename = "result_count_000" + std::to_string(runnumber) + "-"+std::to_string(cutvalue)+".root";
0085     std::string filename = "result_count_00042797-5-30M-merged.root";
0086     // Open the ROOT file
0087     TFile* file = TFile::Open(filename.c_str());
0088     if (!file || file->IsZombie()) {
0089         std::cerr << "Error opening file: " << filename << std::endl;
0090         return;
0091     }
0092     GetSpinPatterns(runnumber);
0093     crossingShift = 0;
0094     // Get the histograms
0095     TH1I *n_smd_v12 = dynamic_cast<TH1I *>(file->Get("h1_NSMD_h_12"));
0096     TH1I *n_smd_v78 = dynamic_cast<TH1I *>(file->Get("h1_NSMD_h_78"));
0097     TH1I *s_smd_v12 = dynamic_cast<TH1I *>(file->Get("h1_SSMD_h_12"));
0098     TH1I *s_smd_v78 = dynamic_cast<TH1I *>(file->Get("h1_SSMD_h_78"));
0099     if (!n_smd_v12 || !n_smd_v78 || !s_smd_v12 || !s_smd_v78)
0100     {
0101         std::cerr << "Error: could not retrieve histograms" << std::endl;
0102         file->Close();
0103         return;
0104     }
0105
0106     // Number of bins
0107     int nBins = n_smd_v12->GetNbinsX();
0108     if (nBins != n_smd_v78->GetNbinsX())
0109     {
0110         std::cerr << "Error: histograms have different number of bins" << std::endl;
0111         file->Close();
0112         return;
0113     }
0114
0115     // Prepare arrays for TGraphErrors
0116     std::vector<double> x_spinup_blue, y_spinup_blue, ex_spinup_blue, ey_spinup_blue;
0117     std::vector<double> x_spindown_blue, y_spindown_blue, ex_spindown_blue, ey_spindown_blue;
0118     std::vector<double> x_spinup_yellow, y_spinup_yellow, ex_spinup_yellow, ey_spinup_yellow;
0119     std::vector<double> x_spindown_yellow, y_spindown_yellow, ex_spindown_yellow, ey_spindown_yellow;
0120
0121     // Calculate the ratio and errors
0122     for (int i = 1; i <= nBins; ++i)
0123     {
0124       double A = n_smd_v12->GetBinContent(i);
0125       double B = n_smd_v78->GetBinContent(i);
0126       double A_y = s_smd_v12->GetBinContent(i);
0127       double B_y = s_smd_v78->GetBinContent(i);
0128       std::cout << "A : " << A << " //// B : " << B << std::endl;
0129       int nbunch = i - 1;
0130       nbunch = (nbunch + crossingShift) % 120;
0131       int b_spin = spinPatternBlue[nbunch];
0132       int y_spin = spinPatternYellow[nbunch];
0133       if (b_spin!=1 && b_spin!=-1)
0134         continue;
0135       if (A != 0 || B != 0)
0136       {
0137         double aym = (B - A) / (A + B);
0138         double error = TMath::Sqrt(4 * (A * B) / ((A + B) * (A + B) * (A + B)));
0139         if (b_spin == 1)
0140         { // spinup_blue bin
0141           x_spinup_blue.push_back(nbunch);
0142           y_spinup_blue.push_back(aym);
0143           ex_spinup_blue.push_back(0); // No error on x
0144           ey_spinup_blue.push_back(error);
0145         }
0146         else if (b_spin == -1)
0147         { // spindown_blue bin
0148           x_spindown_blue.push_back(nbunch);
0149           y_spindown_blue.push_back(aym);
0150           ex_spindown_blue.push_back(0); // No error on x
0151           ey_spindown_blue.push_back(error);
0152         }
0153       }
0154       if (A_y != 0 || B_y != 0)
0155       {
0156         double aym = (B_y - A_y) / (A_y + B_y);
0157         double error = TMath::Sqrt(4 * (A_y * B_y) / ((A_y + B_y) * (A_y + B_y) * (A_y + B_y)));
0158         if (y_spin == 1)
0159         { // spinup_blue bin
0160           x_spinup_yellow.push_back(nbunch);
0161           y_spinup_yellow.push_back(aym);
0162           ex_spinup_yellow.push_back(0); // No error on x
0163           ey_spinup_yellow.push_back(error);
0164         }
0165         else if (y_spin == -1)
0166         { // spindown_yellow bin
0167           x_spindown_yellow.push_back(nbunch);
0168           y_spindown_yellow.push_back(aym);
0169           ex_spindown_yellow.push_back(0); // No error on x
0170           ey_spindown_yellow.push_back(error);
0171         }
0172       }
0173     }
0174
0175     // Create the TGraphErrors for spinup_blue bins
0176     TGraphErrors *graph_spinup_blue = new TGraphErrors(x_spinup_blue.size(), x_spinup_blue.data(), y_spinup_blue.data(), ex_spinup_blue.data(), ey_spinup_blue.data());
0177     graph_spinup_blue->SetMarkerStyle(21);
0178     graph_spinup_blue->SetMarkerColor(kRed);
0179
0180     // Create the TGraphErrors for spindown_blue bins
0181     TGraphErrors *graph_spindown_blue = new TGraphErrors(x_spindown_blue.size(), x_spindown_blue.data(), y_spindown_blue.data(), ex_spindown_blue.data(), ey_spindown_blue.data());
0182     graph_spindown_blue->SetMarkerStyle(21);
0183     graph_spindown_blue->SetMarkerColor(kBlue);
0184
0185     // Create the TGraphErrors for spinup_yellow bins
0186     TGraphErrors *graph_spinup_yellow = new TGraphErrors(x_spinup_yellow.size(), x_spinup_yellow.data(), y_spinup_yellow.data(), ex_spinup_yellow.data(), ey_spinup_yellow.data());
0187     graph_spinup_yellow->SetMarkerStyle(21);
0188     graph_spinup_yellow->SetMarkerColor(kRed);
0189
0190     // Create the TGraphErrors for spindown_yellow bins
0191     TGraphErrors *graph_spindown_yellow = new TGraphErrors(x_spindown_yellow.size(), x_spindown_yellow.data(), y_spindown_yellow.data(), ex_spindown_yellow.data(), ey_spindown_yellow.data());
0192     graph_spindown_yellow->SetMarkerStyle(21);
0193     graph_spindown_yellow->SetMarkerColor(kBlue);
0194
0195
0196
0197
0198     // Draw the graphs
0199     TCanvas *canvas = new TCanvas("canvas_blue", "North SMD", 800, 600);
0200     graph_spinup_blue->SetTitle(Form("NSMD North Crossing Shift : %d;Bunch number;(U-D)/(U+D)",crossingShift));
0201     double x_min = 0;
0202     double x_max = 128;
0203     double y_min = -0.05;
0204     double y_max = 0.05; // Adjust as necessary
0205
0206     graph_spinup_blue->GetXaxis()->SetLimits(x_min, 110);
0207     graph_spinup_blue->SetMinimum(y_min);
0208     graph_spinup_blue->SetMaximum(y_max);
0209     graph_spinup_blue->Draw("AP");
0210     graph_spindown_blue->Draw("P SAME");
0211
0212     // Add legend
0213     TLegend *legend = new TLegend(0.6, 0.8, 0.9, 0.9);
0214     legend->AddEntry(graph_spinup_blue, "blue spin up Bunches", "p");
0215     legend->AddEntry(graph_spindown_blue, "blue spin down Bunches", "p");
0216     legend->Draw();
0217
0218 TCanvas *canvas2 = new TCanvas("canvas_yellow", "South SMD", 800, 600);
0219     graph_spinup_yellow->SetTitle(Form("SSMD South Crossing Shift : %d;Bunch number;(U-D)/(U+D)",crossingShift));
0220     double x_min2 = 0;
0221     double x_max2 = 128;
0222     double y_min2 = -0.05;
0223     double y_max2 = 0.05; // Adjust as necessary
0224
0225     graph_spinup_yellow->GetXaxis()->SetLimits(x_min2, 110);
0226     graph_spinup_yellow->SetMinimum(y_min2);
0227     graph_spinup_yellow->SetMaximum(y_max2);
0228     graph_spinup_yellow->Draw("AP");
0229     graph_spindown_yellow->Draw("P SAME");
0230
0231     // Add legend
0232     TLegend *legend2 = new TLegend(0.6, 0.8, 0.9, 0.9);
0233     legend2->AddEntry(graph_spinup_yellow, "yellow spin up Bunches", "p");
0234     legend2->AddEntry(graph_spindown_yellow, "yellow spin down Bunches", "p");
0235     legend2->Draw();
0236
0237
0238     std::size_t pos = filename.find(".root");
0239
0240     if (pos == std::string::npos) {
0241         return filename;
0242     }
0243     std::string result = filename;
0244     result.insert(pos, "-plot-UD-crossing"+std::to_string(crossingShift));
0245     canvas->SaveAs(Form("run_%d_crossing_%d_blue.png",runnumber,crossingShift));
0246     canvas2->SaveAs(Form("run_%d_crossing_%d_yellow.png",runnumber,crossingShift));
0247     TFile* sf = new TFile(result.c_str(),"RECREATE");
0248     sf->cd();
0249     canvas->Write();
0250     canvas2->Write();
0251     sf->Close();
0252
0253
0254     // Cleanup
0255     file->Close();
0256     delete file;
0257     delete canvas;
0258     delete graph_spinup_yellow;
0259     delete graph_spindown_yellow;
0260 }