bunchAym/macro/DrawBunchAysm.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(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     crossingShift = 0;
0095     // Get the histograms
0096     TH1I *n_smd_v12 = dynamic_cast<TH1I *>(file->Get("h1_NSMD_v_12"));
0097     TH1I *n_smd_v67 = dynamic_cast<TH1I *>(file->Get("h1_NSMD_v_67"));
0098     TH1I *s_smd_v12 = dynamic_cast<TH1I *>(file->Get("h1_SSMD_v_12"));
0099     TH1I *s_smd_v67 = dynamic_cast<TH1I *>(file->Get("h1_SSMD_v_67"));
0100     if (!n_smd_v12 || !n_smd_v67 || !s_smd_v12 || !s_smd_v67)
0101     {
0102         std::cerr << "Error: could not retrieve histograms" << std::endl;
0103         file->Close();
0104         return;
0105     }
0106
0107     // Number of bins
0108     int nBins = n_smd_v12->GetNbinsX();
0109     if (nBins != n_smd_v67->GetNbinsX())
0110     {
0111         std::cerr << "Error: histograms have different number of bins" << std::endl;
0112         file->Close();
0113         return;
0114     }
0115
0116     // Prepare arrays for TGraphErrors
0117     std::vector<double> x_spinup_blue, y_spinup_blue, ex_spinup_blue, ey_spinup_blue;
0118     std::vector<double> x_spindown_blue, y_spindown_blue, ex_spindown_blue, ey_spindown_blue;
0119     std::vector<double> x_spinup_yellow, y_spinup_yellow, ex_spinup_yellow, ey_spinup_yellow;
0120     std::vector<double> x_spindown_yellow, y_spindown_yellow, ex_spindown_yellow, ey_spindown_yellow;
0121
0122     // Calculate the ratio and errors
0123     for (int i = 1; i <= nBins; ++i)
0124     {
0125       double A = n_smd_v12->GetBinContent(i);
0126       double B = n_smd_v67->GetBinContent(i);
0127       double A_y = s_smd_v12->GetBinContent(i);
0128       double B_y = s_smd_v67->GetBinContent(i);
0129       std::cout << "A : " << A << " //// B : " << B << std::endl;
0130       int nbunch = i - 1;
0131       nbunch = (nbunch + crossingShift) % 120;
0132       int b_spin = spinPatternBlue[nbunch];
0133       int y_spin = spinPatternYellow[nbunch];
0134       if (b_spin!=1 && b_spin!=-1)
0135         continue;
0136       if (A != 0 || B != 0)
0137       {
0138         double aym = (A - B) / (A + B);
0139         double error = TMath::Sqrt(4 * (A * B) / ((A + B) * (A + B) * (A + B)));
0140         if (b_spin == 1)
0141         { // spinup_blue bin
0142           x_spinup_blue.push_back(nbunch);
0143           y_spinup_blue.push_back(aym);
0144           ex_spinup_blue.push_back(0); // No error on x
0145           ey_spinup_blue.push_back(error);
0146         }
0147         else if (b_spin == -1)
0148         { // spindown_blue bin
0149           x_spindown_blue.push_back(nbunch);
0150           y_spindown_blue.push_back(aym);
0151           ex_spindown_blue.push_back(0); // No error on x
0152           ey_spindown_blue.push_back(error);
0153         }
0154       }
0155       if (A_y != 0 || B_y != 0)
0156       {
0157         double aym = (A_y - B_y) / (A_y + B_y);
0158         double error = TMath::Sqrt(4 * (A_y * B_y) / ((A_y + B_y) * (A_y + B_y) * (A_y + B_y)));
0159         if (y_spin == 1)
0160         { // spinup_blue bin
0161           x_spinup_yellow.push_back(nbunch);
0162           y_spinup_yellow.push_back(aym);
0163           ex_spinup_yellow.push_back(0); // No error on x
0164           ey_spinup_yellow.push_back(error);
0165         }
0166         else if (y_spin == -1)
0167         { // spindown_yellow bin
0168           x_spindown_yellow.push_back(nbunch);
0169           y_spindown_yellow.push_back(aym);
0170           ex_spindown_yellow.push_back(0); // No error on x
0171           ey_spindown_yellow.push_back(error);
0172         }
0173       }
0174     }
0175
0176     // Create the TGraphErrors for spinup_blue bins
0177     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());
0178     graph_spinup_blue->SetMarkerStyle(21);
0179     graph_spinup_blue->SetMarkerColor(kRed);
0180
0181     // Create the TGraphErrors for spindown_blue bins
0182     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());
0183     graph_spindown_blue->SetMarkerStyle(21);
0184     graph_spindown_blue->SetMarkerColor(kBlue);
0185
0186     // Create the TGraphErrors for spinup_yellow bins
0187     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());
0188     graph_spinup_yellow->SetMarkerStyle(21);
0189     graph_spinup_yellow->SetMarkerColor(kRed);
0190
0191     // Create the TGraphErrors for spindown_yellow bins
0192     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());
0193     graph_spindown_yellow->SetMarkerStyle(21);
0194     graph_spindown_yellow->SetMarkerColor(kBlue);
0195
0196
0197
0198
0199     // Draw the graphs
0200     TCanvas *canvas = new TCanvas("canvas_blue", "North SMD", 800, 600);
0201     graph_spinup_blue->SetTitle(Form("NSMD North Crossing Shift : %d;Bunch number;(L-R)/(L+R)",crossingShift));
0202     double x_min = 0;
0203     double x_max = 128;
0204     double y_min = -0.05;
0205     double y_max = 0.05; // Adjust as necessary
0206
0207     graph_spinup_blue->GetXaxis()->SetLimits(x_min, 110);
0208     graph_spinup_blue->SetMinimum(y_min);
0209     graph_spinup_blue->SetMaximum(y_max);
0210     graph_spinup_blue->Draw("AP");
0211     graph_spindown_blue->Draw("P SAME");
0212
0213   TF1 *fit1 = new TF1("fit1", "[0]", 0, 120);
0214     graph_spinup_blue->Fit(fit1, "R");
0215
0216     TF1 *fit2 = new TF1("fit2", "[0]", 0, 120);
0217     graph_spindown_blue->Fit(fit2, "R");
0218
0219     // Draw the fits
0220     fit1->SetLineColor(kRed);
0221     fit1->Draw("SAME");
0222
0223     fit2->SetLineColor(kBlue);
0224     fit2->Draw("SAME");
0225
0226     // Add legend
0227     TLegend *legend = new TLegend(0.6, 0.8, 0.9, 0.9);
0228     legend->AddEntry(graph_spinup_blue, "blue spin up Bunches", "p");
0229     legend->AddEntry(graph_spindown_blue, "blue spin down Bunches", "p");
0230     legend->Draw();
0231
0232 TCanvas *canvas2 = new TCanvas("canvas_yellow", "South SMD", 800, 600);
0233     graph_spinup_yellow->SetTitle(Form("SSMD South Crossing Shift : %d;Bunch number;(L-R)/(L+R)",crossingShift));
0234     double x_min2 = 0;
0235     double x_max2 = 128;
0236     double y_min2 = -0.05;
0237     double y_max2 = 0.05; // Adjust as necessary
0238
0239     graph_spinup_yellow->GetXaxis()->SetLimits(x_min2, 110);
0240     graph_spinup_yellow->SetMinimum(y_min2);
0241     graph_spinup_yellow->SetMaximum(y_max2);
0242     graph_spinup_yellow->Draw("AP");
0243     graph_spindown_yellow->Draw("P SAME");
0244
0245     // Add legend
0246     TLegend *legend2 = new TLegend(0.6, 0.8, 0.9, 0.9);
0247     legend2->AddEntry(graph_spinup_yellow, "yellow spin up Bunches", "p");
0248     legend2->AddEntry(graph_spindown_yellow, "yellow spin down Bunches", "p");
0249     legend2->Draw();
0250
0251
0252     std::size_t pos = filename.find(".root");
0253
0254     if (pos == std::string::npos) {
0255         return filename;
0256     }
0257     std::string result = filename;
0258     result.insert(pos, "-plot-crossing"+std::to_string(crossingShift));
0259     canvas->SaveAs(Form("run_%d_crossing_%d_blue.png",runnumber,crossingShift));
0260     canvas2->SaveAs(Form("run_%d_crossing_%d_yellow.png",runnumber,crossingShift));
0261     TFile* sf = new TFile(result.c_str(),"RECREATE");
0262     sf->cd();
0263     canvas->Write();
0264     canvas2->Write();
0265     sf->Close();
0266
0267
0268     // Cleanup
0269     file->Close();
0270     delete file;
0271     delete canvas;
0272     delete graph_spinup_yellow;
0273     delete graph_spindown_yellow;
0274 }