LuminosityCounterGoodRuns/run/get_rate_map.C

0001 #include "/sphenix/user/jocl/projects/chi2checker/src/dlUtility.h"
0002
0003 unsigned long long int factorial(int n)
0004 {
0005   if(n < 2) return 1;
0006   return n*factorial(n-1);
0007 }
0008
0009 void increment_in_base(vector<unsigned int> &counter, int base)
0010 {
0011
0012   for(int i = 0; i<counter.size(); ++i)
0013     {
0014       ++counter[i];
0015       if(counter[i] > base)
0016     {
0017       cout << "SOMETHING WRONG - ENTRY > BASE" << endl;
0018       exit(1);
0019     }
0020       if(counter[i] == base)
0021     {
0022       counter[i] = 0;
0023     }
0024       else
0025     {
0026       return;
0027     }
0028     }
0029 }
0030
0031 float get_prob_given_ncol(int n, float mbdprob)
0032 {
0033   float onehitprob = 0.175; //old: sqrt(mbdprob) - mbdprob; //probability of firing only one side
0034   float nohitsprob = 1 - mbdprob - 2*onehitprob; //
0035   vector<unsigned int> counter = {};
0036   float prob = 0;
0037   for(int i=0; i<n; ++i)
0038     {
0039       counter.push_back(0);
0040     }
0041   for(int i=0; i<pow(4,n); ++i)
0042     {
0043       int nnone = std::count(counter.begin(), counter.end(), 0);
0044       int nnorth = std::count(counter.begin(), counter.end(), 1);
0045       int nsouth = std::count(counter.begin(), counter.end(), 2);
0046       int nboth = std::count(counter.begin(), counter.end(), 3);
0047       if((nnorth && nsouth) || nboth)
0048     {
0049       prob += pow(nohitsprob,nnone)*pow(onehitprob,nnorth)*pow(onehitprob,nsouth)*pow(mbdprob,nboth);
0050     }
0051       increment_in_base(counter, 4);
0052     }
0053   return prob;
0054 }
0055
0056 float p_from_rate(float colrate, float beamrate, int n)
0057 {
0058   float num = (1-colrate/beamrate)*pow(-log(1-colrate/beamrate),n);
0059   float den = factorial(n);
0060   return num/den;
0061 }
0062
0063 float p_from_truerate(float truerate, float beamrate, int n)
0064 {
0065   float ratio = truerate/beamrate;
0066   return exp(-ratio)*pow(ratio,n)/factorial(n);
0067 }
0068
0069 float get_true_rate(float colrate, float beamrate)
0070 {
0071   float sum = 0;
0072   for(int i=1; i<10; ++i)
0073     {
0074       sum += pow(-log(1-colrate/beamrate),i)/factorial(i-1);
0075     }
0076   sum *= (beamrate-colrate);
0077   return sum;
0078 }
0079 int get_rate_map()
0080 {
0081   gStyle->SetOptStat(0);
0082   gStyle->SetOptTitle(0);
0083   gStyle->SetPadTickX(1);
0084   gStyle->SetPadTickY(1);
0085   // for a single collision...
0086   float mbdprob = 0.562; //from PYTHIA or use 25.2/42; //MBD xsec / pp inel xsec
0087   const int ncount = 5;
0088   float beamrate = 111*78.2e3;
0089   float p_mbd_given_ncol[ncount];
0090   float p_of_n[ncount];
0091   const int nrate = 400000;
0092   float mbd_rate[nrate];
0093   float mbd_r2[nrate];
0094   float col_rate[nrate];
0095   float true_rate[nrate];
0096   float true_rate_sparse[nrate/100];
0097   float mbd_rate_sparse[nrate/100];
0098   for(int i=0; i<ncount; ++i)
0099     {
0100       p_mbd_given_ncol[i] = get_prob_given_ncol(i+1, mbdprob);
0101       cout << p_mbd_given_ncol[i] << endl;
0102     }
0103   for(int r=0; r<nrate; ++r)
0104     {
0105       col_rate[r] = 10*r;
0106       float total = 0;
0107       float total2 = 0;
0108       true_rate[r] = get_true_rate(col_rate[r],beamrate);
0109       for(int i=0; i<ncount; ++i)
0110     {
0111       total += p_mbd_given_ncol[i]*p_from_rate(col_rate[r], beamrate, i+1);
0112       total2 += p_mbd_given_ncol[i]*p_from_truerate(true_rate[r], beamrate, i+1);
0113     }
0114       if(total > 0.5) cout << total << endl;
0115       mbd_rate[r] = total*beamrate;
0116       mbd_r2[r] = total*beamrate;
0117       if(r%100==0)
0118     {
0119       mbd_rate_sparse[r/100]=mbd_rate[r];
0120       true_rate_sparse[r/100]=true_rate[r];
0121       cout << "total p: " << total - total2 << endl;
0122       cout << "rate diff [Hz]: "  << mbd_rate[r] - mbd_r2[r] << endl;
0123     }
0124     }
0125   TGraph* g = new TGraph(nrate, mbd_rate, col_rate);
0126   g->SetName("mbdtogr1");
0127   g->GetHistogram()->GetXaxis()->SetTitle("R_{MBD} [Hz]");
0128   g->GetHistogram()->GetYaxis()->SetTitle("R_{#geq1} [Hz]");
0129   TGraph* gr = new TGraph(nrate, mbd_rate, true_rate);
0130   gr->SetName("mbdtotrue");
0131   gr->GetHistogram()->GetXaxis()->SetTitle("R_{MBD} [Hz]");
0132   gr->GetHistogram()->GetYaxis()->SetTitle("R_{true} [Hz]");
0133   TGraph* sparse = new TGraph(nrate/100,mbd_rate_sparse,true_rate_sparse);
0134   TH1D* temphist = new TH1D("temphist",";R_{MBD} [Hz];R_{true} [Hz]",1,0,3000e3);
0135   temphist->GetYaxis()->SetRangeUser(0,6000e3);
0136   temphist->GetXaxis()->SetRangeUser(0,3000e3);
0137
0138   sparse->SetMarkerStyle(20);
0139   sparse->SetMarkerSize(1);
0140   sparse->SetMarkerColor(kRed+2);
0141   gr->SetMarkerStyle(20);
0142   gr->SetMarkerSize(1);
0143   gr->SetMarkerColor(kBlue+2);
0144   TCanvas* c = new TCanvas("","",1000,1000);
0145   c->SetLeftMargin(0.15);
0146   c->SetBottomMargin(0.15);
0147   temphist->Draw();
0148   sparse->Draw("SAME P");
0149   //gr->Draw("SAME P");
0150   sphenixtext();
0151   TLegend* leg = new TLegend(0.5,0.2,0.8,0.3);
0152   leg->SetFillStyle(0);
0153   leg->SetFillColor(0);
0154   leg->SetBorderSize(0);
0155   leg->AddEntry(g,"R_{#geq1}","p");
0156   leg->AddEntry(gr,"R_{true}","p");
0157   //leg->Draw();
0158   c->SaveAs("ratevsrate.pdf");
0159   TFile* outf = new TFile("mbd_to_col_map.root","RECREATE");
0160   outf->cd();
0161   g->Write();
0162   gr->Write();
0163   outf->Close();
0164   return 0;
0165 }