macros/scattered-electron_analysis/plot_true_false_positive.C

0001 #include "TGraph.h"
0002 vector< TString > v_cuts_CEMC;
0003 vector< TString > v_cuts_CEMC_short;
0004 vector< TString > v_cuts_EEMC;
0005 vector< TString > v_cuts_EEMC_short;
0006 vector< TString > v_cuts_FEMC;
0007 vector< TString > v_cuts_FEMC_short;
0008
0009 void add_CEMC_Cut(double, double);
0010
0011 void add_EEMC_Cut(double);
0012
0013 void add_FEMC_Cut(double);
0014
0015 int plot_true_false_positive()
0016 {
0017
0018   /*----------------------------------------- */
0019   /* Add Cuts Here                            */
0020   /*                                          */
0021   /* add_CEMC_cut(float_1, float_2)           */
0022   /*   Cut 1 : E/p  > float_1                 */
0023   /*   Cut 2 : prob > float_2                 */
0024   /*   (if float_#=0, Cut # is not applied    */
0025   /* add_EEMC_cut(float_1)                    */
0026   /*   Cut 1 : E/True E > float_1             */
0027   /* add_FEMC_cut(float_1)                    */
0028   /*   Cut 1 : E/True E > float_1             */
0029   /* ---------------------------------------- */
0030
0031   add_CEMC_Cut(0.807,0.02);
0032   add_CEMC_Cut(0.703,0.02);
0033
0034   add_EEMC_Cut(0.703);
0035   add_EEMC_Cut(0.807);
0036   add_EEMC_Cut(0.95);
0037
0038   add_FEMC_Cut(0.703);
0039   add_FEMC_Cut(0.807);
0040   add_FEMC_Cut(0.95);
0041
0042   /* ---------------------------------------- */
0043
0044   // Momentums //
0045   vector< float > v_momenta;
0046   v_momenta.push_back(1.0);
0047   v_momenta.push_back(2.0);
0048   v_momenta.push_back(5.0);
0049   v_momenta.push_back(10.0);
0050   v_momenta.push_back(20.0);
0051
0052   // Marker Styles //
0053   vector< int > v_styles;
0054   v_styles.push_back(20);
0055   v_styles.push_back(24);
0056   v_styles.push_back(21);
0057   v_styles.push_back(25);
0058   v_styles.push_back(22);
0059   v_styles.push_back(26);
0060   v_styles.push_back(23);
0061   v_styles.push_back(27);
0062   v_styles.push_back(34);
0063   v_styles.push_back(28);
0064   v_styles.push_back(29);
0065   v_styles.push_back(30);
0066
0067   /* Loading Macros */
0068
0069   gROOT->LoadMacro("efficiency_purity_analysisMacro.C");
0070   gROOT->LoadMacro("/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C");
0071   SetsPhenixStyle();
0072
0073   /* Useful local variables */
0074   float v_rate_tp; //true positive rate value holder
0075   float v_rate_fp; //false positive rate value holder
0076
0077   //------------------------------------------------------------//
0078   //--------                CEMC Plots                  --------//
0079   //------------------------------------------------------------//
0080
0081   //CEMC Graph Initialization//
0082   const int size_C=2*v_cuts_CEMC.size(); //T.P & F.P graph per cut
0083   TCanvas *c_CEMC = new TCanvas("cc","cc",1000,600);
0084   TGraph *graphs_CEMC[size_C]; //Array of TGraph pointers
0085   auto legend_CEMC_t = new TLegend(0.51,0.2,1.0,0.9); //Legend
0086   auto legend_CEMC_f = new TLegend(0.51,0.2,1.0,0.9); //Legend
0087
0088   legend_CEMC_t->SetTextSize(0.06);
0089   legend_CEMC_f->SetTextSize(0.06);
0090   // Create the CEMC TGraphs ( 2 per cut )//
0091   for(unsigned idx_cut=0;idx_cut < v_cuts_CEMC.size(); idx_cut++)
0092     /* loop over cuts */
0093     {
0094       // Index for each graph in TGraph vector //
0095       int idx_tp=idx_cut*2;
0096       int idx_fp=idx_cut*2+1;
0097
0098       // Create new TGraphs_CEMC //
0099       graphs_CEMC[idx_tp]=new TGraph(v_momenta.size());
0100       graphs_CEMC[idx_fp]=new TGraph(v_momenta.size());
0101
0102       // Clear true positive and false positive //
0103       v_rate_tp = v_rate_fp = 0.0;
0104
0105       // Set Points in TGraph for True Positive and False Positive per momenta//
0106       for ( unsigned idx_p = 0; idx_p < v_momenta.size(); idx_p++ )
0107     {
0108       v_rate_tp=get_true_positive("C",v_cuts_CEMC.at(idx_cut),v_momenta.at(idx_p));
0109       v_rate_fp=get_false_positive("C",v_cuts_CEMC.at(idx_cut),v_momenta.at(idx_p));
0110       graphs_CEMC[idx_tp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_tp);
0111       graphs_CEMC[idx_fp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_fp);
0112     }
0113
0114       // Tweak Graph Style //
0115       graphs_CEMC[idx_tp]->SetMarkerColor(kGreen+2);
0116       graphs_CEMC[idx_fp]->SetMarkerColor(kRed+2);
0117
0118       graphs_CEMC[idx_tp]->SetLineColor(graphs_CEMC[idx_tp]->GetMarkerColor());
0119       graphs_CEMC[idx_fp]->SetLineColor(graphs_CEMC[idx_fp]->GetMarkerColor());
0120
0121       graphs_CEMC[idx_tp]->SetMarkerStyle(v_styles.at(idx_cut));
0122       graphs_CEMC[idx_fp]->SetMarkerStyle(v_styles.at(idx_cut));
0123
0124       graphs_CEMC[idx_tp]->SetMarkerSize(1); // Possibly too large? //
0125       graphs_CEMC[idx_fp]->SetMarkerSize(1);
0126
0127       graphs_CEMC[idx_tp]->SetTitle("");
0128       graphs_CEMC[idx_fp]->SetTitle("");
0129
0130
0131
0132       legend_CEMC_t->AddEntry(graphs_CEMC[idx_tp],v_cuts_CEMC_short.at(idx_cut)+" True Positive","P");
0133       legend_CEMC_f->AddEntry(graphs_CEMC[idx_fp],v_cuts_CEMC_short.at(idx_cut)+ " False Positive","P");
0134     }
0135
0136   // Draw all graphs on same canvas (CEMC) //
0137   c_CEMC->Divide(1,2);
0138   for(int i = 0; i<size_C; i++)
0139     {
0140       if(i==0)
0141     {
0142       c_CEMC->cd(1);
0143       //gPad->SetBottomMargin(0.2);
0144       gPad->SetRightMargin(0.5);
0145       graphs_CEMC[0]->Draw("PA");
0146       graphs_CEMC[0]->GetXaxis()->SetTitle("True Momentum (GeV)");
0147       graphs_CEMC[0]->GetYaxis()->SetRangeUser(0.7,1.0);
0148       graphs_CEMC[0]->GetYaxis()->SetTitle("Rate");
0149     }
0150       else if(i==1)
0151     {
0152       c_CEMC->cd(2);
0153       //gPad->SetTopMargin(0.2);
0154       gPad->SetRightMargin(0.5);
0155       graphs_CEMC[1]->GetXaxis()->SetTitle("True Momentum (GeV)");
0156       graphs_CEMC[1]->GetYaxis()->SetRangeUser(-0.01,0.11);
0157       graphs_CEMC[1]->GetYaxis()->SetTitle("Rate");
0158       graphs_CEMC[1]->Draw("PA");
0159     }
0160       else if (i%2==0)
0161     {
0162       c_CEMC->cd(1);
0163       graphs_CEMC[i]->Draw("PSame");
0164     }
0165       else
0166     {
0167       c_CEMC->cd(2);
0168       graphs_CEMC[i]->Draw("P");
0169     }
0170
0171     }
0172   c_CEMC->cd(1);
0173   legend_CEMC_t->Draw();
0174   c_CEMC->cd(2);
0175   legend_CEMC_f->Draw();
0176
0177   //------------------------------------------------------------//
0178   //--------                EEMC Plots                  --------//
0179   //------------------------------------------------------------//
0180
0181
0182   // Create EEMC Graphs (2 per cut) //
0183    const int size_E=v_cuts_EEMC.size(); //T.P & F.P graph per cut
0184    TCanvas *c_EEMC = new TCanvas("ce","ce",1000,600);
0185    TGraph *graphs_EEMC[size_E]; //Array of TGraph pointers
0186    auto legend_EEMC_t = new TLegend(0.51,0.2,1,0.9); //Legend
0187    auto legend_EEMC_f = new TLegend(0.51,0.2,1,0.9); //Legend
0188
0189    legend_EEMC_t->SetTextSize(0.06);
0190    legend_EEMC_f->SetTextSize(0.06);
0191    for(unsigned idx_cut=0;idx_cut < v_cuts_EEMC.size()/2; idx_cut++)
0192     // loop over cuts //
0193     {
0194       // Index for each graph in TGraph vector //
0195       int idx_tp=idx_cut*2;
0196       int idx_fp=idx_cut*2+1;
0197
0198       // Create new TGraphs //
0199       graphs_EEMC[idx_tp]=new TGraph(v_momenta.size());
0200       graphs_EEMC[idx_fp]=new TGraph(v_momenta.size());
0201
0202       // Clear true positive and false positive //
0203       v_rate_tp = v_rate_fp = 0.0;
0204
0205       // Fill data points for True Positve and False Positive Plot //
0206       for ( unsigned idx_p = 0; idx_p < v_momenta.size(); idx_p++ )
0207     {
0208       TString cut = v_cuts_EEMC.at(idx_cut*2)+ getCut(v_momenta.at(idx_p)) + v_cuts_EEMC.at(idx_cut*2+1);
0209       v_rate_tp=get_true_positive("E",cut,v_momenta.at(idx_p));
0210       v_rate_fp=get_false_positive("E",cut,v_momenta.at(idx_p));
0211       graphs_EEMC[idx_tp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_tp);
0212       graphs_EEMC[idx_fp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_fp);
0213     }
0214
0215       graphs_EEMC[idx_tp]->SetMarkerColor(kGreen+2);
0216       graphs_EEMC[idx_fp]->SetMarkerColor(kRed+2);
0217
0218       graphs_EEMC[idx_tp]->SetLineColor(graphs_EEMC[idx_tp]->GetMarkerColor());
0219       graphs_EEMC[idx_fp]->SetLineColor(graphs_EEMC[idx_fp]->GetMarkerColor());
0220
0221       graphs_EEMC[idx_tp]->SetMarkerStyle(v_styles.at(idx_cut));
0222       graphs_EEMC[idx_fp]->SetMarkerStyle(v_styles.at(idx_cut));
0223
0224       graphs_EEMC[idx_tp]->SetMarkerSize(1); // Possibly too large? //
0225       graphs_EEMC[idx_fp]->SetMarkerSize(1);
0226
0227       graphs_EEMC[idx_tp]->SetTitle("");
0228       graphs_EEMC[idx_fp]->SetTitle("");
0229
0230       legend_EEMC_t->AddEntry(graphs_EEMC[idx_tp],v_cuts_EEMC_short.at(idx_cut)+" True Positive","P");
0231       legend_EEMC_f->AddEntry(graphs_EEMC[idx_tp],v_cuts_EEMC_short.at(idx_cut)+" False Positive","P");
0232     }
0233     // Draw all graphs on same canvas (EEMC) //
0234    c_EEMC->Divide(1,2);
0235    for(int i = 0; i<size_E; i++)
0236      {
0237        if(i==0)
0238      {
0239        c_EEMC->cd(1);
0240        gPad->SetRightMargin(0.5);
0241        graphs_EEMC[0]->Draw("PA");
0242        graphs_EEMC[0]->GetXaxis()->SetTitle("True Momentum (GeV)");
0243        graphs_EEMC[0]->GetYaxis()->SetRangeUser(0.7,1.0);
0244        graphs_EEMC[0]->GetYaxis()->SetTitle("Rate");
0245      }
0246        else if (i==1)
0247      {
0248        c_EEMC->cd(2);
0249        gPad->SetRightMargin(0.5);
0250        graphs_EEMC[1]->Draw("PA");
0251        graphs_EEMC[1]->GetXaxis()->SetTitle("True Momentum (GeV)");
0252        graphs_EEMC[1]->GetYaxis()->SetRangeUser(-0.01,0.11);
0253        graphs_EEMC[1]->GetYaxis()->SetTitle("Rate");
0254      }
0255        else if (i%2==0)
0256      {
0257        c_EEMC->cd(1);
0258        graphs_EEMC[i]->Draw("PSame");
0259      }
0260        else
0261      {
0262        c_EEMC->cd(2);
0263        graphs_EEMC[i]->Draw("P");
0264      }
0265      }
0266    c_EEMC->cd(1);
0267    legend_EEMC_t->Draw();
0268    c_EEMC->cd(2);
0269    legend_EEMC_f->Draw();
0270
0271    //------------------------------------------------------------//
0272    //--------                FEMC Plots                  --------//
0273    //------------------------------------------------------------//
0274
0275
0276    // Create FEMC Graphs (2 per cut) //
0277    const int size_F=v_cuts_FEMC.size(); //T.P & F.P graph per cut
0278    TCanvas *c_FEMC = new TCanvas("cf","cf",1000,600);
0279    TGraph *graphs_FEMC[size_F]; //Array of TGraph pointers
0280    auto legend_FEMC_t = new TLegend(0.51,0.2,1.0,0.9); //Legend
0281    auto legend_FEMC_f = new TLegend(0.51,0.2,1.0,0.9); //Legend
0282
0283    legend_FEMC_t->SetTextSize(0.06);
0284    legend_FEMC_f->SetTextSize(0.06);
0285
0286    for(unsigned idx_cut=0;idx_cut < v_cuts_FEMC.size()/2; idx_cut++)
0287     // loop over cuts //
0288      {
0289        // Index for each graph in TGraph vector //
0290       int idx_tp=idx_cut*2;
0291       int idx_fp=idx_cut*2+1;
0292
0293       // Create new TGraphs //
0294       graphs_FEMC[idx_tp]=new TGraph(v_momenta.size());
0295       graphs_FEMC[idx_fp]=new TGraph(v_momenta.size());
0296
0297       // Clear true positive and false positive //
0298       v_rate_tp = v_rate_fp = 0.0;
0299
0300       // Fill data points for True Positve and False Positive Plot //
0301       for ( unsigned idx_p = 0; idx_p < v_momenta.size(); idx_p++ )
0302     {
0303       TString cut = v_cuts_FEMC.at(idx_cut*2)+ getCut(v_momenta.at(idx_p)) + v_cuts_FEMC.at(idx_cut*2+1);
0304       v_rate_tp=get_true_positive("F",cut,v_momenta.at(idx_p));
0305       v_rate_fp=get_false_positive("F",cut,v_momenta.at(idx_p));
0306       graphs_FEMC[idx_tp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_tp);
0307       graphs_FEMC[idx_fp]->SetPoint(idx_p,v_momenta[idx_p],v_rate_fp);
0308     }
0309
0310       graphs_FEMC[idx_tp]->SetMarkerColor(kGreen+2);
0311       graphs_FEMC[idx_fp]->SetMarkerColor(kRed+2);
0312
0313       graphs_FEMC[idx_tp]->SetLineColor(graphs_FEMC[idx_tp]->GetMarkerColor());
0314       graphs_FEMC[idx_fp]->SetLineColor(graphs_FEMC[idx_fp]->GetMarkerColor());
0315
0316       graphs_FEMC[idx_tp]->SetMarkerStyle(v_styles.at(idx_cut));
0317       graphs_FEMC[idx_fp]->SetMarkerStyle(v_styles.at(idx_cut));
0318
0319       graphs_FEMC[idx_tp]->SetMarkerSize(1); // Possibly too large? //
0320       graphs_FEMC[idx_fp]->SetMarkerSize(1);
0321
0322       graphs_FEMC[idx_tp]->SetTitle("");
0323       graphs_FEMC[idx_fp]->SetTitle("");
0324
0325       legend_FEMC_t->AddEntry(graphs_FEMC[idx_tp],v_cuts_FEMC_short.at(idx_cut)+" True Positive","P");
0326       legend_FEMC_f->AddEntry(graphs_FEMC[idx_fp],v_cuts_FEMC_short.at(idx_cut)+ " False Positive","P");
0327     }
0328
0329    // Draw all graphs on same canvas (FEMC) //
0330    c_FEMC->Divide(1,2);
0331    for(int i = 0; i<size_F; i++)
0332      {
0333         if(i==0)
0334     {
0335       c_FEMC->cd(1);
0336       //gPad->SetBottomMargin(0.2);
0337       gPad->SetRightMargin(0.5);
0338       graphs_FEMC[0]->Draw("PA");
0339       graphs_FEMC[0]->GetXaxis()->SetTitle("True Momentum (GeV)");
0340       graphs_FEMC[0]->GetYaxis()->SetRangeUser(0.7,1.0);
0341       graphs_FEMC[0]->GetYaxis()->SetTitle("Rate");
0342     }
0343       else if(i==1)
0344     {
0345       c_FEMC->cd(2);
0346       //gPad->SetTopMargin(0.2);
0347       gPad->SetRightMargin(0.5);
0348       graphs_FEMC[1]->GetXaxis()->SetTitle("True Momentum (GeV)");
0349       graphs_FEMC[1]->GetYaxis()->SetRangeUser(-0.01,0.11);
0350       graphs_FEMC[1]->GetYaxis()->SetTitle("Rate");
0351       graphs_FEMC[1]->Draw("PA");
0352     }
0353       else if (i%2==0)
0354     {
0355       c_FEMC->cd(1);
0356       graphs_FEMC[i]->Draw("PSame");
0357     }
0358       else
0359     {
0360       c_FEMC->cd(2);
0361       graphs_FEMC[i]->Draw("PSame");
0362     }
0363      }
0364    c_FEMC->cd(1);
0365    legend_FEMC_t->Draw();
0366    c_FEMC->cd(2);
0367    legend_FEMC_f->Draw();
0368
0369   return 0;
0370 }
0371
0372 void add_CEMC_Cut(double ep, double prob)
0373 {
0374   //TStrings necessary to represent cuts//
0375
0376   TString base_ep = "em_cluster_e/em_track_ptotal > ";
0377   TString base_prob = "em_cluster_prob >= ";
0378
0379   //TString necessary to print legible legends//
0380
0381   TString base_ep_short = "E/P > ";
0382   TString base_prob_short = "Prob >= ";
0383   char str_1[40];
0384   char str_2[40];
0385   sprintf(str_1,"%.2f", ep);
0386   sprintf(str_2,"%.3f", prob);
0387
0388   //Test which cuts to use//
0389
0390   bool do_ep=true; bool do_prob=true;
0391   if(ep==0)
0392     do_ep=false;
0393   if(prob==0)
0394     do_prob=false;
0395
0396   //Add cuts to vector//
0397
0398   if(do_prob&&!do_ep)
0399     {
0400       v_cuts_CEMC.push_back(TString(base_prob) += prob);
0401       v_cuts_CEMC_short.push_back((TString(base_prob_short)+=str_2)+" True Positive");
0402     }
0403   else if(do_ep&&!do_prob)
0404     {
0405       v_cuts_CEMC.push_back(TString(base_ep) += ep);
0406       v_cuts_CEMC_short.push_back((TString(base_ep_short)+=str_1)+" True Positive");
0407     }
0408   else
0409     {
0410       v_cuts_CEMC.push_back((TString(base_ep) += ep) + "&&" + (TString(base_prob) += prob));
0411       v_cuts_CEMC_short.push_back((TString(base_ep_short) += str_1) + "&&" +(TString(base_prob_short) += str_2));
0412     }
0413 }
0414
0415 void add_EEMC_Cut(double ep)
0416 {
0417   TString base_ep="em_cluster_e/";
0418   TString base_ep_short = "E/P > ";
0419   char str_1[40];
0420   sprintf(str_1,"%.2f", ep);
0421
0422   v_cuts_EEMC.push_back(base_ep); // 1st element : 'em_cluster_e/'
0423   v_cuts_EEMC.push_back(TString(">")+= str_1); //: 'em_cluster_e/*GeV*>cut'
0424
0425   v_cuts_EEMC_short.push_back(((TString(base_ep_short)) += str_1));
0426 }
0427
0428 void add_FEMC_Cut(double ep)
0429 {
0430   TString base_ep="em_cluster_e/";
0431   TString base_ep_short = "E/P > ";
0432   char str_1[40];
0433   sprintf(str_1,"%.2f", ep);
0434
0435   v_cuts_FEMC.push_back(base_ep); // 1st element : 'em_cluster_e/'
0436   v_cuts_FEMC.push_back(TString(">")+= str_1); //: 'em_cluster_e/*GeV*>cut'
0437
0438   v_cuts_FEMC_short.push_back(((TString(base_ep_short)) += str_1));
0439 }