LuminosityCounterGoodRuns/run/get_luminosity_182630.C

0001 /*
0002 int get_scaler(int runnumber, int bit)
0003 {
0004   TSQLServer *db = TSQLServer::Connect("pgsql://sphnxdaqdbreplica:5432/daq","","");
0005
0006   if (db)
0007     {
0008       printf("Server info: %s\n", db->ServerInfo());
0009     }
0010   else
0011     {
0012       printf("bad\n");
0013     }
0014
0015   TSQLRow *row;
0016   TSQLResult *res;
0017   TString cmd = "";
0018   char sql[1000];
0019
0020   sprintf(sql, "select scaled from gl1_scalers where runnumber = %d and index = %d;", runnumber, bit);
0021   res = db->Query(sql);
0022   row = res->Next();
0023   try
0024     {
0025       int retval = stoi(row->GetField(0));
0026       delete row;
0027       delete res;
0028       delete db;
0029       return retval;
0030     }
0031   catch(...)
0032     {
0033       delete row;
0034       delete res;
0035       delete db;
0036       return 0;
0037     }
0038 }
0039 */
0040 int get_luminosity_182630(string rnlist, int zsel)
0041 {
0042   gStyle->SetOptStat(0);
0043
0044   if(zsel > 5 || zsel < 0)
0045     {
0046       cout << "z selection must be between 0 and 5 (inclusive) for |zvtx| < 30/60/200/10/1000/none, respectively! Exiting." << endl;
0047       return 2;
0048     }
0049   const int ntrig = 3;
0050   float mbsig = 25.2;
0051   int upperrun = 53000;
0052   int lowerrun = 47200;
0053   ifstream is(rnlist);
0054   string rnstr;
0055   TH1D* mbdrate = new TH1D("mbdrate","",upperrun-lowerrun,lowerrun-.5,upperrun-.5);
0056   TH1D* nSB10 = new TH1D("nSB10","",upperrun-lowerrun,lowerrun-0.5,upperrun-0.5);
0057   TH1D* nSB18 = new TH1D("nSB18","",upperrun-lowerrun,lowerrun-0.5,upperrun-0.5);
0058   TH1D* fsb10l30 = new TH1D("fsb10l30","",upperrun-lowerrun,lowerrun-0.5,upperrun-0.5);
0059   TH1D* fsb18l30 = new TH1D("fsb18l30","",upperrun-lowerrun,lowerrun-0.5,upperrun-0.5);
0060   //int it = 0;
0061   float lumi[ntrig] = {0};
0062   float uclumi[ntrig] = {0};
0063   int trigs[ntrig] = {18,22,30};
0064   cout << "RN   Bit"<<to_string(trigs[0])<<"Corr  Bit"<<to_string(trigs[0])<<"UC   Bit"<<to_string(trigs[1])<<"Corr  Bit"<<to_string(trigs[1])<<"UC   Bit"<<to_string(trigs[2])<<"Corr  Bit"<<to_string(trigs[2])<<"UC (units pb^-1)" << endl << endl;
0065   while(getline(is, rnstr))
0066     {
0067       //cerr << "RN: " << rnstr << endl;
0068       /*
0069       try
0070     {
0071       if(stoi(rnstr) > 48862 && stoi(rnstr) < 48868) continue;
0072     }
0073       catch(...)
0074     {
0075       continue;
0076     }
0077       */
0078       TFile* file;
0079       file = TFile::Open(("output/added_nblair/triggeroutput_nblair_"+rnstr+".root").c_str(),"READ");
0080       if(!file)
0081     {
0082       cerr << "RN " << rnstr << " not found! Continuing." << endl;
0083       cout << rnstr << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << endl;
0084       continue;
0085     }
0086       long long unsigned int tottrigcounts[6][64];
0087       double avgPS[64];
0088       long long unsigned int sumgoodscaled[64];
0089       long long unsigned int sumgoodraw[64];
0090       long long unsigned int sumgoodlive[64];
0091       float nmbdc;
0092       float ttseg;
0093       float nblair;
0094       TTree* tree = (TTree*)file->Get("outt");
0095
0096       if(!tree)
0097     {
0098       cerr << "RN " << rnstr << " does not contain tree outt! Continuing." << endl;
0099       cout << rnstr << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << endl;
0100       continue;
0101     }
0102       tree->SetBranchAddress("avgPS",avgPS);
0103       tree->SetBranchAddress("sumgoodscaled",sumgoodscaled);
0104       tree->SetBranchAddress("sumgoodlive",sumgoodlive);
0105       tree->SetBranchAddress("sumgoodraw",sumgoodraw);
0106       tree->SetBranchAddress("tottrigcounts",tottrigcounts);
0107       tree->SetBranchAddress("nmbdc",&nmbdc);
0108       tree->SetBranchAddress("ttseg",&ttseg);
0109       tree->SetBranchAddress("nblair",&nblair);
0110       tree->GetEntry(0);
0111
0112
0113       /*
0114       try
0115     {
0116       int scaleddb18 = get_scaler(stoi(rnstr),18);
0117       /*
0118       if(sumgoodscaled[18] < 0.99*scaleddb18)
0119         {
0120           cerr << sumgoodscaled[18] << " " << scaleddb18 << endl;
0121           cerr << "< 99% of db" << endl;
0122           cout << "Run " << rnstr << " " << 0 << " " <<0 << endl;
0123           continue;
0124         }
0125     }
0126
0127       catch(...)
0128     {
0129       cout << "Run " << rnstr << " " << 0 << " " << 0 << endl;
0130       cerr << "error with getting db counts" << endl;
0131       continue;
0132     }
0133       */
0134       //tree->SetBranchAddress("runLumi",runLumi);
0135       /*
0136     if(sumgoodlive[10]/sumgoodscaled[10] > avgPS[10]+0.0000000001)
0137     {
0138       cout << rnstr << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " " << 0 << endl;
0139       continue;
0140     }
0141       */
0142       //cout << rnstr << " " << ((double)tottrigcounts[zsel][10])*avgPS[10]/avgPS[trigger] << endl;
0143       if(nmbdc/((double)sumgoodraw[10]) < 1)
0144     {
0145       cerr << rnstr << " " << nmbdc << " " << sumgoodraw[10] << " " << nmbdc / (double)sumgoodraw[10] << endl;
0146       continue;
0147     }
0148       cout << rnstr;
0149
0150
0151
0152       for(int i=0; i<3; ++i)
0153     {
0154       int trigger = trigs[i];
0155       if(avgPS[trigger] > 0 && !isnan(avgPS[trigger]) && !isnan(avgPS[10]) && !isnan((tottrigcounts[zsel][10])) && !isinf(avgPS[10]) && !isinf(avgPS[trigger]) && !isinf((tottrigcounts[zsel][10])) && avgPS[10] > 0)
0156         {
0157
0158           float clumiseg = (tottrigcounts[zsel][10])*avgPS[10]*((((nmbdc/((double)sumgoodraw[10])))))/avgPS[trigger]/(mbsig*1e9);
0159           float uclumiseg = (tottrigcounts[zsel][10])*((double)avgPS[10])/avgPS[trigger]/(mbsig*1e9);//((double)(trigger==22?sumgoodscaled[10]:tottrigcounts[zsel][10]))*avgPS[10]/avgPS[trigger]/(mbdsig*1e9;
0160           float blumiseg = (tottrigcounts[zsel][10])*avgPS[10]*((((nblair/((double)sumgoodraw[10]))-1)/2)+1)/avgPS[trigger]/(mbsig*1e9);
0161           uclumi[i] += uclumiseg;
0162           lumi[i] += clumiseg;
0163           cout << " " << fixed << setprecision(7) << clumiseg << " " << uclumiseg;// << " " << blumiseg ;
0164         }
0165       else
0166         {
0167           cout << " " << 0 << " " << 0;// << " " << 0 ;
0168         }
0169       //if(runLumi > 0) lumi += runLumi;
0170       //cout << lumi << " " << (trigger==22?sumgoodscaled[10]:tottrigcounts[zsel][10]) << " " << avgPS[10] << " " << avgPS[trigger] << endl;
0171     }
0172       cout << endl;
0173     }
0174   /*
0175   lumi[0] /= ((mbdsig*1e9);
0176   lumi[1] /= (mbdsig*1e9;
0177   lumi[2] /= (mbdsig*1e9;
0178   uclumi[0] /= ((mbdsig*1e9);
0179   uclumi[1] /= (mbdsig*1e9;
0180   uclumi[2] /= (mbdsig*1e9;
0181   */
0182   /*
0183   formatHist(nSB10,0,nSB10->GetMaximum()*1.1,"Run Number","N_{event} Firing MBDNS>=1 Trigger",kBlack,1.5);
0184   formatHist(nSB18,0,nSB18->GetMaximum()*1.1,"Run Number","N_{event} Firing Jet10+MBDNS>=1 Trigger",kBlack,1.5);
0185   formatHist(fsb10l30,0,1,"Run Number","N_{ev}MBDNS>=1 w/ |z_{vtx}|<30cm / N_{ev}MBDNS>=1",kBlack,1.5);
0186   formatHist(fsb18l30,0,1,"Run Number","N_{ev}Jet10+MBDNS>=1 w/ |z_{vtx}|<30cm / N_{ev}Jet10+MBDNS>=1",kBlack,1.5);
0187   */
0188   cout << "RN   Bit"<<to_string(trigs[0])<<"Corr  Bit"<<to_string(trigs[0])<<"UC   Bit"<<to_string(trigs[1])<<"Corr  Bit"<<to_string(trigs[1])<<"UC   Bit"<<to_string(trigs[2])<<"Corr  Bit"<<to_string(trigs[2])<<"UC (units pb^-1)" << endl << endl;
0189   cout << "Totals: " << fixed << setprecision(5)<<  lumi[0] << " " << uclumi[0] << " " << lumi[1] << " " << uclumi[1] << " " << lumi[2] << " " << uclumi[2] << " " << " pb^-1" << endl;
0190   /*
0191   mbdrate->SetMarkerStyle(20);
0192   mbdrate->SetMarkerSize(2);
0193   TCanvas* c1 = new TCanvas("","",2000,1000);
0194   gPad->SetLeftMargin(0.15);
0195   gPad->SetRightMargin(0.01);
0196   gPad->SetTopMargin(0.05);
0197   gPad->SetBottomMargin(0.15);
0198   mbdrate->GetYaxis()->SetTitle("MBD Live Rate [Hz]");
0199   mbdrate->GetXaxis()->SetTitle("Run Number");
0200   mbdrate->Draw("P HIST");
0201   c1->SaveAs("mbdrate_full.png");
0202   nSB10->Draw("P HIST");
0203   c1->SaveAs("nsb10.png");
0204   nSB18->Draw("P HIST");
0205   c1->SaveAs("nsb18.png");
0206   fsb10l30->Draw("P HIST");
0207   c1->SaveAs("fsb10l30.png");
0208   fsb18l30->Draw("P HIST");
0209   c1->SaveAs("fsb18l30.png");
0210   */
0211   return 0;
0212 }