hcal_calib_year2/fileQA/plot_hits.C

0001 #include "BlairUtils.C"
0002 #include "sPhenixStyle.C"
0003
0004 std::vector<int> runNumbers;    // Stores the run numbers
0005 std::vector<double> timestamps;  // Stores the timestamps
0006 std::vector<double> values;     // Stores the corresponding values
0007 std::vector<int> hours;     // Stores the corresponding values
0008
0009 std::vector<int> runNumbers_cosmic;    // Stores the run numbers
0010 std::vector<double> timestamps_cosmic;  // Stores the timestamps
0011
0012 void getInfo();
0013 void getCosmicInfo();
0014 float corrDouble(float hits_per_event,double rate);
0015
0016 void plot_hits() {
0017
0018
0019    SetsPhenixStyle();
0020
0021   getInfo();
0022   getCosmicInfo();
0023
0024
0025   TFile* fin = new TFile("output/anaOut.root");
0026   TH1F* h_events_run = (TH1F*) fin->Get("h_events_run");
0027   TH1F* h_hits= (TH1F*) fin->Get("h_hits");
0028   TH1F* h_runnumber = (TH1F*) fin->Get("runNumberHist");
0029   TH1F* h_temp_run = (TH1F*) fin->Get("h_temp_run");
0030
0031   TH1F * h_hits_per_event = (TH1F*) h_hits->Clone("h_hits_per_event");
0032   h_hits_per_event->Sumw2();
0033   h_events_run->Sumw2();
0034   h_hits_per_event->Divide(h_events_run);
0035
0036    int n_runs = h_runnumber->GetNbinsX();
0037
0038   //TH1F* h_hits_time = new TH1F("h_hits_time","",n_runs
0039   TGraphErrors *gr_hits_time = new TGraphErrors();
0040   TGraphErrors *gr_hits_run  = new TGraphErrors();
0041   TGraphErrors *gr_rate_time = new TGraphErrors();
0042   TGraphErrors *gr_hits_time_corr = new TGraphErrors();
0043   TProfile* pr_temp_hour = new TProfile("pr_temp_hour","",24,-0.5,23.5);
0044   TGraphErrors* gr_temp_time = new TGraphErrors();
0045   TH1F* h_hits_corr_dis = new TH1F("h_hits_corr_dis","",100,0.002,0.003);
0046
0047    float mean_corr_hits = 0;
0048    int   mean_corr_hitsC = 0;
0049
0050   for(int it=0; it<n_runs; it++){
0051     int runnumber = h_runnumber->GetBinContent(1+it);
0052     for (int ir=0; ir<runNumbers.size(); ir++){
0053       if(runnumber == runNumbers[ir]){
0054         float hits = h_hits_per_event->GetBinContent(it+1);
0055         float hits_err = h_hits_per_event->GetBinError(it+1);
0056         gr_hits_time->SetPoint(it,timestamps[ir],hits);
0057         gr_hits_time->SetPointError(it,0,hits_err);
0058
0059         gr_hits_run->SetPoint(it,runnumber,hits);
0060         gr_hits_run->SetPointError(it,0,hits_err);
0061
0062         gr_rate_time->SetPoint(it,timestamps[ir],values[ir]);
0063         cout << "timestamps[ir] " << timestamps[ir] << " values[ir]] " << values[ir] << endl;
0064
0065         pr_temp_hour->Fill(hours[ir],h_temp_run->GetBinContent(1+it));
0066         gr_temp_time->SetPoint(it,timestamps[ir],h_temp_run->GetBinContent(1+it));
0067
0068         if (values[ir] != 0){
0069            float corr_hits = corrDouble(hits,values[ir]);
0070            gr_hits_time_corr->SetPoint(it,timestamps[ir],corr_hits);
0071            gr_hits_time_corr->SetPointError(it,0,hits_err);
0072            mean_corr_hits += corr_hits; mean_corr_hitsC++;
0073            h_hits_corr_dis->Fill(corr_hits);
0074         }
0075       }
0076
0077     }
0078   }
0079
0080
0081    mean_corr_hits /= mean_corr_hitsC;
0082    TGraph* gr_mean_corr_hits = new TGraph();
0083    gr_mean_corr_hits->SetPoint(0,0,mean_corr_hits);
0084    gr_mean_corr_hits->SetPoint(1,1e10,mean_corr_hits);
0085    gr_mean_corr_hits->SetLineColor(kBlue);
0086
0087
0088    TGraph* gr_cosmic_times = new TGraph();
0089
0090
0091   for(int ir=0; ir<runNumbers_cosmic.size(); ir++){
0092     cout << timestamps_cosmic[ir] << endl;
0093     gr_cosmic_times->SetPoint(ir*2,timestamps_cosmic[ir],10000*(ir%2));
0094     gr_cosmic_times->SetPoint(ir*2+1,timestamps_cosmic[ir],10000*((ir+1)%2));
0095   }
0096
0097    TCanvas *c1 = new TCanvas("c1","c1",600,600);
0098    //h_hits_per_event->Draw("ex0");
0099    gr_hits_run->Draw("ap");
0100
0101     TCanvas* c2 = new TCanvas("c2","c2",2000,600);
0102     TH1F* h_frame = new TH1F("h_frame","",100,1.7183e+09,1.729e+09);
0103     h_frame->GetXaxis()->SetTimeFormat("%m-%d %H:%M");
0104     h_frame->GetXaxis()->SetTimeDisplay(1);
0105     h_frame->GetXaxis()->SetTitle("run time");
0106     h_frame->GetYaxis()->SetTitle("mean hits / event / run");
0107     h_frame->GetXaxis()->SetRangeUser(1.7183e+09,1.729e+09);
0108     //h_frame->GetXaxis()->SetRangeUser(1.72e+09,1.721e+09);
0109     h_frame->GetYaxis()->SetRangeUser(0,0.005);
0110     h_frame->GetYaxis()->SetTitleOffset(0.7);
0111     h_frame->Draw("axis");
0112     gPad->SetLeftMargin(0.08);
0113
0114     gr_hits_time->Draw("p");
0115     gr_hits_time->SetMarkerStyle(20);
0116     gr_hits_time->SetMarkerSize(1.5);
0117     gr_hits_time->GetXaxis()->SetTimeDisplay(1);
0118     gr_hits_time->GetYaxis()->SetRangeUser(0,0.005);
0119     gr_hits_time->GetXaxis()->SetNdivisions(-503);
0120     gr_hits_time->GetXaxis()->SetTimeFormat("%m-%d %H:%M");
0121     gr_hits_time->GetXaxis()->SetTitle("run time");
0122     gr_hits_time->GetYaxis()->SetTitle("mean hits / event / run");
0123     //gr_hits_time->GetXaxis()->SetRangeUser(1.7185e+09,1.729e+09);
0124     gr_hits_time->GetXaxis()->SetRangeUser(1.72e+09,1.721e+09);
0125
0126     gr_hits_time_corr->Draw("p");
0127     gr_hits_time_corr->SetMarkerStyle(33);;
0128     gr_hits_time_corr->SetMarkerColor(kBlue);;
0129     gr_hits_time_corr->SetLineColor(kBlue);;
0130 /*
0131     gr_cosmic_times->Draw("l");
0132     gr_cosmic_times->SetLineStyle(7);
0133     gr_cosmic_times->SetLineWidth(2);
0134     gr_cosmic_times->SetLineColor(kRed);
0135     gr_cosmic_times->GetXaxis()->SetTimeDisplay(1);
0136     gr_cosmic_times->GetXaxis()->SetNdivisions(-503);
0137     gr_cosmic_times->GetXaxis()->SetTimeFormat("%m-%d %H:%M");
0138
0139    gr_mean_corr_hits->Draw("l");
0140 */
0141    myText         (0.11,0.90-0.5,1,"#bf{#it{sPHENIX}} Internal",0.05);
0142    myMarkerLineText(0.14,0.85-0.5,1.5,kBlack ,20,kBlack,1,"data", 0.05, true);
0143    myMarkerLineText(0.14,0.80-0.5,1.5,kBlue ,33,kBlue,1,"pileup corrected", 0.05, true);
0144    //myMarkerLineText(0.14,0.75-0.5,0  ,kRed,1,kRed,7,"cosmic calib runs", 0.05, true);
0145      myText         (0.5,0.90-0.5,1,"Not final calib - no temp correction or cosmics",0.05);
0146
0147    c2->SaveAs("figures/hit_rate_vs_time.pdf");
0148
0149
0150    TCanvas* c11 = new TCanvas("c1","c1",600,600);
0151    h_hits_corr_dis->Draw();
0152    h_hits_corr_dis->Fit("gaus");
0153    h_hits_corr_dis->SetXTitle("Hits/run");
0154    h_hits_corr_dis->SetYTitle("runs");
0155    h_hits_corr_dis->GetXaxis()->SetNdivisions(505);
0156
0157    TF1* fitFunc = h_hits_corr_dis->GetFunction("gaus");
0158    double fmean = fitFunc->GetParameter(1);      // Mean (mu)
0159    double fsigma = fitFunc->GetParameter(2);     // Standard deviation (sigma)
0160
0161    myText(0.20, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.05);
0162    myText(0.20, 0.85, 1, "OHCal", 0.05);
0163    myText(0.20, 0.8, 1, Form("Fit #mu/#sigma = %.3f", fsigma/fmean), 0.05);
0164    myText(0.20, 0.75, 1, Form("RMS/Mean = %.3f", h_hits_corr_dis->GetRMS()/h_hits_corr_dis->GetMean()), 0.05);
0165
0166    c11->SaveAs("figures/hit_rate_allRuns.pdf");
0167
0168   ///////////////////////
0169   // temp plots
0170   /////////////////////////
0171
0172   TProfile2D* pr2d_temp     = (TProfile2D*) fin->Get("pr2d_temp");
0173   //TH2F* pr2d_temp_var = (TH2F*) pr2d_temp->Clone("pr2d_temp_var");
0174   TH2F* pr2d_temp_var = new TH2F("pr2d_temp_var","",24,0,24,64,0,64);
0175   pr2d_temp_var->Reset();
0176   for(int ix=0; ix<24; ix++){
0177     for(int iy=0; iy<64; iy++){
0178       float var = pr2d_temp->GetBinError(ix+1,iy+1);
0179       cout << var << endl;
0180       pr2d_temp_var->SetBinContent(ix+1,iy+1,var);
0181       pr2d_temp_var->SetBinError(ix+1,iy+1,0);
0182     }
0183   }
0184
0185    TCanvas *c3 = new TCanvas("c3","c3",600,600);
0186    //pr_temp_hour->Draw();
0187    pr2d_temp->Draw("COLZ");
0188    pr2d_temp->SetXTitle("OHCal #it{#eta}^{i}");
0189    pr2d_temp->SetYTitle("OHCal #it{#phi}^{i}");
0190    gPad->SetRightMargin(0.15);
0191
0192   TCanvas *c4 = new TCanvas("c4","c4",2000,600);
0193   TH1F* h_frame2 = new TH1F("h_frame2","",100,1.7183e+09,1.729e+09);
0194   h_frame2->GetXaxis()->SetTimeFormat("%m-%d %H:%M");
0195   h_frame2->GetXaxis()->SetTimeDisplay(1);
0196   h_frame2->GetXaxis()->SetTitle("run time");
0197   h_frame2->GetYaxis()->SetTitle("mean OHCal temperature / run");
0198   h_frame2->GetXaxis()->SetRangeUser(1.7183e+09,1.729e+09);
0199   h_frame2->GetYaxis()->SetRangeUser(0,25);
0200   h_frame2->GetYaxis()->SetTitleOffset(0.7);
0201   h_frame2->Draw("axis");
0202   gPad->SetLeftMargin(0.08);
0203
0204   gr_temp_time->Draw("p");
0205   gr_temp_time->SetMarkerStyle(20);
0206   gr_temp_time->SetMarkerSize(1.0);
0207   gr_temp_time->GetXaxis()->SetTimeDisplay(1);
0208   gr_temp_time->GetYaxis()->SetRangeUser(0,25);
0209   gr_temp_time->GetXaxis()->SetNdivisions(-503);
0210   gr_temp_time->GetXaxis()->SetTimeFormat("%m-%d %H:%M");
0211   gr_temp_time->GetXaxis()->SetTitle("run time");
0212   gr_temp_time->GetYaxis()->SetTitle("ohcal tower avg temp [C]");
0213   gr_temp_time->GetXaxis()->SetRangeUser(1.7185e+09,1.729e+09);
0214
0215   myText         (0.11,0.90-0.5,1,"#bf{#it{sPHENIX}} Internal",0.05);
0216
0217   c4->SaveAs("figures/OHCal_mean_temp.pdf");
0218
0219 }
0220
0221
0222 float corrDouble(float hits_per_event,double rate){
0223   //double newnumber = TMath::Abs(lambda * (double)(endtime-starttime)*(111*78e3));
0224   //double prob = newnumber / 78e3/111;
0225   double prob1orMore = rate / 78.e3/111.;
0226   double lambda = -1.0*TMath::Log(1-prob1orMore);
0227   // labda = prob of a single collision
0228   // p_double = prob of an additional collsion which = lambda
0229   // x = true probability of a hit per a single mb event
0230   // obs_x is the observed probability (prob per bunch crossing)
0231   // (x*events +  x*event*p_double ) / events = ob_x
0232   // x + x*p_double = ob_x
0233   // x = ob_x / (1+p_double)
0234   return hits_per_event / (1 + lambda);
0235 }
0236
0237
0238
0239
0240 void getInfo(){
0241     // Define vectors to store the data
0242
0243     // Open the CSV file
0244     std::ifstream file("run_timestamps_with_rate.csv");  // Replace with your actual CSV file path
0245     if (!file.is_open()) {
0246         std::cerr << "Error opening file!" << std::endl;
0247         return;
0248     }
0249
0250     std::string line;
0251     while (std::getline(file, line)) {
0252         std::istringstream ss(line);
0253         std::string runNumberStr, timestampStr, valueStr;
0254
0255         // Read the run number, timestamp, and value from each line
0256         std::getline(ss, runNumberStr, ',');
0257         std::getline(ss, timestampStr, ',');
0258         std::getline(ss, valueStr, ',');
0259
0260         // Convert the run number to integer
0261         int runNumber = std::stoi(runNumberStr);
0262
0263
0264
0265         int year, month, day, hour, minute, second;
0266         sscanf(timestampStr.c_str(), "%d-%d-%d %d:%d:%d", &year, &month, &day, &hour, &minute, &second);
0267         TDatime t(year, month, day, hour, minute, second);
0268         double timeStamp = t.Convert();
0269
0270         // Convert the value to double
0271         double value = valueStr.empty() ? 0.0 : std::stod(valueStr);
0272
0273        // cout << runNumber << "  "   << timeStamp << "  "  << value << endl;
0274         // Store the values in the vectors
0275         runNumbers.push_back(runNumber);
0276         timestamps.push_back(timeStamp);
0277         values.push_back(value);
0278         hours.push_back(t.GetHour());
0279     }
0280
0281     file.close();
0282
0283
0284 }
0285
0286
0287
0288
0289 void getCosmicInfo(){
0290     std::ifstream file("cosmic_times.csv");  // Replace with your actual CSV file path
0291     if (!file.is_open()) {
0292         std::cerr << "Error opening file!" << std::endl;
0293         return;
0294     }
0295
0296     std::string line;
0297     while (std::getline(file, line)) {
0298         std::istringstream ss(line);
0299         std::string runNumberStr, timestampStr, valueStr;
0300
0301         std::getline(ss, runNumberStr, ',');
0302         std::getline(ss, timestampStr, ',');
0303         std::getline(ss, valueStr, ',');
0304
0305         int runNumber = std::stoi(runNumberStr);
0306
0307         int year, month, day, hour, minute, second;
0308         sscanf(timestampStr.c_str(), "%d-%d-%d %d:%d:%d", &year, &month, &day, &hour, &minute, &second);
0309         TDatime t(year, month, day, hour, minute, second);
0310         double timeStamp = t.Convert();
0311
0312          cout << runNumber << "  "   << timeStamp << endl;
0313
0314         runNumbers_cosmic.push_back(runNumber);
0315         timestamps_cosmic.push_back(timeStamp);
0316     }
0317
0318     file.close();
0319
0320
0321 }
0322
0323
0324
0325
0326
0327
0328
0329
0330