calo/calo_syst/genCaloSyst.C

0001 #include "cdbHistConv.C"
0002
0003 #include <TCanvas.h>
0004 #include <TRandom3.h>
0005
0006 #include <format>
0007
0008 void genCaloSyst()
0009 {
0010   // std::string calo = "CEMC";
0011   // std::string calo = "HCALOUT";
0012   // std::string calo = "HCALIN";
0013   std::vector<std::string> calos = {"CEMC", "HCALOUT", "HCALIN"};
0014
0015   for (auto calo : calos)
0016   {
0017     // int begin = 54908;
0018     // int end   = 54921;
0019     int begin = 46000;
0020     int end   = 53862;
0021
0022     float tbyt_stat;
0023     float calib_unc;
0024     float sys_v1;
0025 //    float sys_v2;
0026     float had_frac;
0027     int eta;
0028     int phi;
0029     bool embool; // embool = 0 for EMCal shape, embool = 1 for HCal shape
0030     float had_resp_unc = 0.117;  // value from test beam
0031
0032     if (calo == "CEMC")
0033     {
0034       tbyt_stat = 0.03;  // tower by tower stat uncertainty
0035       calib_unc = 0.026;     // global scale shift
0036       sys_v1 = 0.02;
0037       //      sys_v2 = 0.02;
0038       had_frac = 0.61;
0039       embool = false;
0040       eta = 96;
0041       phi = 256;
0042     }
0043     if (calo == "HCALOUT")
0044     {
0045       tbyt_stat = 0.02;  // tower by tower stat uncertainty
0046       calib_unc = 0.02;      // global scale shift
0047       sys_v1 = 0.02;
0048       //      sys_v2 = 0.02;
0049       had_frac = 0.972;
0050       embool = true;
0051       eta = 24;
0052       phi = 64;
0053     }
0054     if (calo == "HCALIN")
0055     {
0056       tbyt_stat = 0.05;  // tower by tower stat uncertainty
0057       calib_unc = 0.05;       // global scale shift
0058       sys_v1 = 0.02;
0059       //      sys_v2 = 0.05;
0060       had_frac = 0.961;
0061       embool = true;
0062       eta = 24;
0063       phi = 64;
0064     }
0065
0066     TH2F* h_sys = new TH2F("h_sys","",eta,0,eta,phi,0,phi);
0067     TH2F* h_calib_stat      = (TH2F*) h_sys->Clone("h_calib_stat");
0068     TH2F* h_calib_unc_up     = (TH2F*) h_sys->Clone("h_calib_unc_up");
0069     TH2F* h_v1        = (TH2F*) h_sys->Clone("h_v1");
0070     TH2F* h_v2        = (TH2F*) h_sys->Clone("h_v2");
0071     TH2F* h_had_resp_up  = (TH2F*) h_sys->Clone("h_had_resp_up");
0072     TH2F* h_no_calib  = (TH2F*) h_sys->Clone("h_no_calib");
0073     TH2F* h_calib_unc_down = (TH2F*) h_sys->Clone("h_calib_unc_down");
0074     TH2F* h_had_resp_down = (TH2F*) h_sys->Clone("h_had_resp_down");
0075
0076     TRandom3* rnd = new TRandom3();
0077
0078     for (int ie = 1; ie < h_sys->GetXaxis()->GetNbins() + 1; ie++)
0079     {
0080       for (int ip = 1; ip < h_sys->GetYaxis()->GetNbins() + 1; ip++)
0081       {
0082         float val = rnd->Gaus(1.0, tbyt_stat);
0083         h_calib_stat->SetBinContent(ie, ip, val);
0084         h_calib_unc_up->SetBinContent(ie, ip, 1.0 + calib_unc);
0085         float v1 = 1.0 + sys_v1 * std::cos((float) ip / phi * 2. * M_PI);
0086         h_v1->SetBinContent(ie, ip, v1);
0087         float v2 = 1.0 + sys_v1 * std::cos(2.0 * (float) ip / phi * 2. * M_PI);
0088         h_v2->SetBinContent(ie, ip, v2);
0089         float had_resp_up = 1.0 + 2.0 * had_resp_unc * had_frac / (sqrt(12));
0090         h_had_resp_up->SetBinContent(ie,ip,had_resp_up);
0091         h_no_calib->SetBinContent(ie,ip,1.0);
0092         h_calib_unc_down->SetBinContent(ie, ip, 1.0 - calib_unc);
0093         float had_resp_down = 1.0 - 2.0 * had_resp_unc * had_frac / (sqrt(12));
0094         h_had_resp_down->SetBinContent(ie, ip, had_resp_down);
0095       }
0096     }
0097
0098     histToCaloCDBTree(std::format("cdbFiles/{}_calib_stat_syst_{}_{}.root"        ,calo,begin,end),"calo_sys", embool, h_calib_stat);
0099     histToCaloCDBTree(std::format("cdbFiles/{}_calib_unc_up_syst_{}_{}.root"       ,calo,begin,end),"calo_sys", embool, h_calib_unc_up);
0100     histToCaloCDBTree(std::format("cdbFiles/{}_v1Modulation_syst_{}_{}.root",calo,begin,end),"calo_sys", embool, h_v1);
0101   //  histToCaloCDBTree(std::format("cdbFilesAuAu/{}_v2Modulation_syst_{}_{}.root",calo.c_str(),begin,end),"calo_sys", embool, h_v2);
0102     histToCaloCDBTree(std::format("cdbFiles/{}_had_resp_up_syst_{}_{}.root"    ,calo,begin,end),"calo_sys", embool, h_had_resp_up);
0103     histToCaloCDBTree(std::format("cdbFiles/{}_no_calib_syst_{}_{}.root"    ,calo,begin,end),"calo_sys", embool, h_no_calib);
0104     histToCaloCDBTree(std::format("cdbFiles/{}_calib_unc_down_syst_{}_{}.root"    ,calo,begin,end),"calo_sys", embool, h_calib_unc_down);
0105     histToCaloCDBTree(std::format("cdbFiles/{}_had_resp_down_syst_{}_{}.root" ,calo,begin,end),"calo_sys", embool, h_had_resp_down);
0106
0107     TCanvas* c1 = new TCanvas("c1","c1",600,600);
0108     h_had_resp_up->Draw("COLZ");
0109     c1->Draw();
0110     c1->SaveAs(std::format("h_had_resp_up_{}.png", calo).c_str());
0111
0112     TCanvas* c2 = new TCanvas("c2", "c2", 600, 600);
0113     h_calib_unc_up->Draw("COLZ");
0114     c2->Draw();
0115     c2->SaveAs(std::format("h_calib_unc_up_{}.png", calo).c_str());
0116
0117     TCanvas* c3 = new TCanvas("c3", "c3", 600, 600);
0118     h_calib_stat->Draw("COLZ");
0119     c3->Draw();
0120     c3->SaveAs(std::format("h_calib_stat_{}.png", calo).c_str());
0121
0122     TCanvas* c4 = new TCanvas("c4", "c4", 600, 600);
0123     h_v1->Draw("COLZ");
0124     c4->Draw();
0125     c4->SaveAs(std::format("h_v1_{}.png",calo).c_str());
0126
0127     TCanvas* c5 = new TCanvas("c5","c5",600,600);
0128     h_no_calib->Draw("COLZ");
0129     c5->Draw();
0130     c5->SaveAs(std::format("h_no_calib_{}.png",calo).c_str());
0131
0132     TCanvas* c6 = new TCanvas("c6","c6",600,600);
0133     h_calib_unc_down->Draw("COLZ");
0134     c6->Draw();
0135     c6->SaveAs(std::format("h_calib_unc_down_{}.png",calo).c_str());
0136
0137     TCanvas* c7 = new TCanvas("c7","c7",600,600);
0138     h_had_resp_down->Draw("COLZ");
0139     c7->Draw();
0140     c7->SaveAs(std::format("h_had_resp_down_{}.png",calo).c_str());
0141
0142   };
0143 }