Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:10 

0001 #include <TFile.h>
0002 #include <TLine.h>
0003 #include <TString.h>
0004 #include <TTree.h>
0005 #include <cassert>
0006 #include <cmath>
0007 #include "SaveCanvas.C"
0008 #include "sPhenixStyle.C"
0009 using namespace std;
0010 
0011 class lin_res
0012 {
0013  public:
0014   TString name;
0015   TGraphErrors *linearity;
0016   TGraphErrors *resolution;
0017   TF1 *f_res;
0018 };
0019 
0020 
0021 
0022 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_nightly/";
0023 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_hcaldigi/";
0024 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_hcaldigi_sampling/";
0025 //TString base_dataset ="/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap";
0026 TString base_dataset = "/sphenix/u/weihuma/analysis/Calorimeter/macros/";
0027 
0028 void DrawHadronShowers(void)
0029 {
0030   SetsPhenixStyle();
0031   TVirtualFitter::SetDefaultFitter("Minuit2");
0032 
0033   ClusterSizeScan("eta0", "0<|#eta|<0.1");
0034   //  ClusterSizeScan("eta0.60","0.6<|#eta|<0.7");
0035  // PIDScan();
0036   //EtaScan();
0037 }
0038 
0039 void EtaScan(void)
0040 {
0041   const TString config = "Single #pi^{-}, Track-based 5x5 clusterizer";
0042 
0043   lin_res *QA_Eta0 =
0044       GetResolution("pi-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0045   QA_Eta0->name = "0<|eta|<0.1";
0046 
0047   lin_res *QA_Eta3 =
0048       GetResolution("pi-", "eta0.30", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0049   QA_Eta3->name = "0.3<|eta|<0.4";
0050 
0051   lin_res *QA_Eta6 =
0052       GetResolution("pi-", "eta0.60", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0053   QA_Eta6->name = "0.6<|eta|<0.7";
0054 
0055   lin_res *QA_Eta9 =
0056       GetResolution("pi-", "eta0.90", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0057   QA_Eta9->name = "0.9<|eta|<1.0";
0058 
0059   // ------------
0060   QA_Eta0->resolution->SetMarkerStyle(kFullCircle);
0061   QA_Eta0->resolution->SetMarkerColor(kRed + 3);
0062   QA_Eta0->resolution->SetLineColor(kRed + 3);
0063 
0064   QA_Eta0->linearity->SetMarkerStyle(kFullCircle);
0065   QA_Eta0->linearity->SetMarkerColor(kRed + 3);
0066   QA_Eta0->linearity->SetLineColor(kRed + 3);
0067 
0068   QA_Eta0->f_res->SetLineColor(kRed + 3);
0069 
0070   // ------------
0071   QA_Eta3->resolution->SetMarkerStyle(kFullSquare);
0072   QA_Eta3->resolution->SetMarkerColor(kBlue + 3);
0073   QA_Eta3->resolution->SetLineColor(kBlue + 3);
0074 
0075   QA_Eta3->linearity->SetMarkerStyle(kFullSquare);
0076   QA_Eta3->linearity->SetMarkerColor(kBlue + 3);
0077   QA_Eta3->linearity->SetLineColor(kBlue + 3);
0078 
0079   QA_Eta3->f_res->SetLineColor(kBlue + 3);
0080 
0081   // ------------
0082   QA_Eta6->resolution->SetMarkerStyle(kFullCross);
0083   QA_Eta6->resolution->SetMarkerColor(kMagenta + 3);
0084   QA_Eta6->resolution->SetLineColor(kMagenta + 3);
0085 
0086   QA_Eta6->linearity->SetMarkerStyle(kFullCross);
0087   QA_Eta6->linearity->SetMarkerColor(kMagenta + 3);
0088   QA_Eta6->linearity->SetLineColor(kMagenta + 3);
0089 
0090   QA_Eta6->f_res->SetLineColor(kMagenta + 3);
0091   QA_Eta6->f_res->SetLineStyle(kDashed);
0092 
0093   // ------------
0094   QA_Eta9->resolution->SetMarkerStyle(kStar);
0095   QA_Eta9->resolution->SetMarkerColor(kCyan + 3);
0096   QA_Eta9->resolution->SetLineColor(kCyan + 3);
0097 
0098   QA_Eta9->linearity->SetMarkerStyle(kStar);
0099   QA_Eta9->linearity->SetMarkerColor(kCyan + 3);
0100   QA_Eta9->linearity->SetLineColor(kCyan + 3);
0101 
0102   QA_Eta9->f_res->SetLineColor(kCyan + 3);
0103   QA_Eta9->f_res->SetLineStyle(kDashed);
0104 
0105   TCanvas *c1 = new TCanvas(Form("EtaScan"),
0106                             Form("EtaScan"), 1300, 600);
0107   c1->Divide(2, 1);
0108   int idx = 1;
0109   TPad *p;
0110 
0111   p = (TPad *) c1->cd(idx++);
0112   c1->Update();
0113   //  p->SetLogy();
0114 
0115   p->DrawFrame(0, 0, 60, 1.2,
0116                Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
0117   TLine *l = new TLine(0, 1, 60, 1);
0118   l->SetLineWidth(2);
0119   l->SetLineColor(kGray);
0120   l->Draw();
0121 
0122   QA_Eta0->linearity->Draw("ep");
0123   QA_Eta3->linearity->Draw("ep");
0124   QA_Eta6->linearity->Draw("ep");
0125   QA_Eta9->linearity->Draw("ep");
0126 
0127   TLegend *leg = new TLegend(.2, .2, .85, .45);
0128   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0129   leg->AddEntry("", "", "");
0130   leg->AddEntry(QA_Eta0->linearity, QA_Eta0->name, "ep");
0131   leg->AddEntry(QA_Eta3->linearity, QA_Eta3->name, "ep");
0132   leg->AddEntry(QA_Eta6->linearity, QA_Eta6->name, "ep");
0133   leg->AddEntry(QA_Eta9->linearity, QA_Eta9->name, "ep");
0134   leg->Draw();
0135 
0136   p = (TPad *) c1->cd(idx++);
0137   c1->Update();
0138 
0139   TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
0140                               32);
0141   f_calo_t1044->SetParameters(13.5, 64.9);
0142   f_calo_t1044->SetLineWidth(3);
0143   f_calo_t1044->SetLineColor(kGreen + 2);
0144 
0145   TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
0146                              Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
0147 
0148   QA_Eta0->f_res->Draw("same");
0149   QA_Eta0->resolution->Draw("ep");
0150 
0151   QA_Eta3->f_res->Draw("same");
0152   QA_Eta3->resolution->Draw("ep");
0153 
0154   QA_Eta6->f_res->Draw("same");
0155   QA_Eta6->resolution->Draw("ep");
0156 
0157   QA_Eta9->f_res->Draw("same");
0158   QA_Eta9->resolution->Draw("ep");
0159 
0160   f_calo_t1044->Draw("same");
0161 
0162   TLegend *leg = new TLegend(.2, .6, .85, .9);
0163   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0164   leg->AddEntry("", "", "");
0165   leg->AddEntry(QA_Eta0->linearity,
0166                 QA_Eta0->name + "      " +
0167                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0168                          QA_Eta0->f_res->GetParameter(0),
0169                          QA_Eta0->f_res->GetParameter(1)),
0170                 "lpe");
0171   leg->AddEntry(f_calo_t1044,
0172                 TString("T-1044 (#eta=0)") + "   " +
0173                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0174                          f_calo_t1044->GetParameter(0),
0175                          f_calo_t1044->GetParameter(1)),
0176                 "l");
0177   leg->AddEntry(QA_Eta3->linearity,
0178                 QA_Eta3->name + "  " +
0179                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0180                          QA_Eta3->f_res->GetParameter(0),
0181                          QA_Eta3->f_res->GetParameter(1)),
0182                 "lpe");
0183   leg->AddEntry(QA_Eta6->linearity,
0184                 QA_Eta6->name + "  " +
0185                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0186                          QA_Eta6->f_res->GetParameter(0),
0187                          QA_Eta6->f_res->GetParameter(1)),
0188                 "lpe");
0189   leg->AddEntry(QA_Eta9->linearity,
0190                 QA_Eta9->name + "  " +
0191                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0192                          QA_Eta9->f_res->GetParameter(0),
0193                          QA_Eta9->f_res->GetParameter(1)),
0194                 "lpe");
0195   leg->Draw();
0196 
0197   SaveCanvas(c1,
0198              base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
0199 }
0200 
0201 void PIDScan(void)
0202 {
0203   const TString config = "Single Particles, 0<|#eta|<0.1, Track-based 5x5 clusterizer";
0204 
0205   lin_res *QA_electron =
0206       GetResolution("e-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0207   QA_electron->name = "Electron";
0208 
0209   lin_res *QA_PiNeg =
0210       GetResolution("pi-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0211   QA_PiNeg->name = "#pi^{-}";
0212 
0213   lin_res *QA_PiPos =
0214       GetResolution("pi\\+", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0215   QA_PiPos->name = "#pi^{+}";
0216 
0217   lin_res *QA_Proton =
0218       GetResolution("proton", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0219   QA_Proton->name = "Proton";
0220 
0221   // ------------
0222   QA_electron->resolution->SetMarkerStyle(kFullCircle);
0223   QA_electron->resolution->SetMarkerColor(kRed + 3);
0224   QA_electron->resolution->SetLineColor(kRed + 3);
0225 
0226   QA_electron->linearity->SetMarkerStyle(kFullCircle);
0227   QA_electron->linearity->SetMarkerColor(kRed + 3);
0228   QA_electron->linearity->SetLineColor(kRed + 3);
0229 
0230   QA_electron->f_res->SetLineColor(kRed + 3);
0231 
0232   // ------------
0233   QA_PiNeg->resolution->SetMarkerStyle(kFullSquare);
0234   QA_PiNeg->resolution->SetMarkerColor(kBlue + 3);
0235   QA_PiNeg->resolution->SetLineColor(kBlue + 3);
0236 
0237   QA_PiNeg->linearity->SetMarkerStyle(kFullSquare);
0238   QA_PiNeg->linearity->SetMarkerColor(kBlue + 3);
0239   QA_PiNeg->linearity->SetLineColor(kBlue + 3);
0240 
0241   QA_PiNeg->f_res->SetLineColor(kBlue + 3);
0242 
0243   // ------------
0244   QA_PiPos->resolution->SetMarkerStyle(kFullCross);
0245   QA_PiPos->resolution->SetMarkerColor(kMagenta + 3);
0246   QA_PiPos->resolution->SetLineColor(kMagenta + 3);
0247 
0248   QA_PiPos->linearity->SetMarkerStyle(kFullCross);
0249   QA_PiPos->linearity->SetMarkerColor(kMagenta + 3);
0250   QA_PiPos->linearity->SetLineColor(kMagenta + 3);
0251 
0252   QA_PiPos->f_res->SetLineColor(kMagenta + 3);
0253   QA_PiPos->f_res->SetLineStyle(kDashed);
0254 
0255   // ------------
0256   QA_Proton->resolution->SetMarkerStyle(kStar);
0257   QA_Proton->resolution->SetMarkerColor(kCyan + 3);
0258   QA_Proton->resolution->SetLineColor(kCyan + 3);
0259 
0260   QA_Proton->linearity->SetMarkerStyle(kStar);
0261   QA_Proton->linearity->SetMarkerColor(kCyan + 3);
0262   QA_Proton->linearity->SetLineColor(kCyan + 3);
0263 
0264   QA_Proton->f_res->SetLineColor(kCyan + 3);
0265   QA_Proton->f_res->SetLineStyle(kDashed);
0266 
0267   TCanvas *c1 = new TCanvas(Form("PIDScan"),
0268                             Form("PIDScan"), 1300, 600);
0269   c1->Divide(2, 1);
0270   int idx = 1;
0271   TPad *p;
0272 
0273   p = (TPad *) c1->cd(idx++);
0274   c1->Update();
0275   //  p->SetLogy();
0276 
0277   p->DrawFrame(0, 0, 60, 1.2,
0278                Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
0279   TLine *l = new TLine(0, 1, 60, 1);
0280   l->SetLineWidth(2);
0281   l->SetLineColor(kGray);
0282 
0283   l->Draw();
0284 
0285   QA_electron->linearity->Draw("ep");
0286   QA_PiNeg->linearity->Draw("ep");
0287   QA_PiPos->linearity->Draw("ep");
0288   QA_Proton->linearity->Draw("ep");
0289 
0290   TLegend *leg = new TLegend(.2, .2, 1, .45);
0291   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0292   leg->AddEntry("", "", "");
0293   leg->AddEntry(QA_electron->linearity, QA_electron->name, "ep");
0294   leg->AddEntry(QA_PiNeg->linearity, QA_PiNeg->name, "ep");
0295   leg->AddEntry(QA_PiPos->linearity, QA_PiPos->name, "ep");
0296   leg->AddEntry(QA_Proton->linearity, QA_Proton->name, "ep");
0297   leg->Draw();
0298 
0299   p = (TPad *) c1->cd(idx++);
0300   c1->Update();
0301 
0302   TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
0303                               32);
0304   f_calo_t1044->SetParameters(13.5, 64.9);
0305   f_calo_t1044->SetLineWidth(3);
0306   f_calo_t1044->SetLineColor(kGreen + 2);
0307 
0308   TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
0309                              Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
0310 
0311   QA_electron->f_res->Draw("same");
0312   QA_electron->resolution->Draw("ep");
0313 
0314   QA_PiNeg->f_res->Draw("same");
0315   QA_PiNeg->resolution->Draw("ep");
0316 
0317   QA_PiPos->f_res->Draw("same");
0318   QA_PiPos->resolution->Draw("ep");
0319 
0320   QA_Proton->f_res->Draw("same");
0321   QA_Proton->resolution->Draw("ep");
0322 
0323   f_calo_t1044->Draw("same");
0324 
0325   TLegend *leg = new TLegend(.2, .6, 1, .9);
0326   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0327   leg->AddEntry("", "", "");
0328   leg->AddEntry(QA_electron->linearity,
0329                 QA_electron->name + "      " +
0330                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0331                          QA_electron->f_res->GetParameter(0),
0332                          QA_electron->f_res->GetParameter(1)),
0333                 "lpe");
0334   leg->AddEntry(f_calo_t1044,
0335                 TString("T-1044 Data, #pi^{-}") + "  " +
0336                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0337                          f_calo_t1044->GetParameter(0),
0338                          f_calo_t1044->GetParameter(1)),
0339                 "l");
0340   leg->AddEntry(QA_PiNeg->linearity,
0341                 QA_PiNeg->name + "      " +
0342                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0343                          QA_PiNeg->f_res->GetParameter(0),
0344                          QA_PiNeg->f_res->GetParameter(1)),
0345                 "lpe");
0346   leg->AddEntry(QA_PiPos->linearity,
0347                 QA_PiPos->name + "  " +
0348                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0349                          QA_PiPos->f_res->GetParameter(0),
0350                          QA_PiPos->f_res->GetParameter(1)),
0351                 "lpe");
0352   leg->AddEntry(QA_Proton->linearity,
0353                 QA_Proton->name + "  " +
0354                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0355                          QA_Proton->f_res->GetParameter(0),
0356                          QA_Proton->f_res->GetParameter(1)),
0357                 "lpe");
0358   leg->Draw();
0359 
0360   SaveCanvas(c1,
0361              base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
0362 }
0363 
0364 void ClusterSizeScan(TString seta = "eta0", TString eta_desc = "0<|#eta|<0.1")
0365 {
0366   const TString config = "Single #pi^{-}, " + eta_desc;
0367 
0368   lin_res *QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP =
0369       GetResolution("pi-", seta, "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP");
0370   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name = "3x3 cluster";
0371 
0372   lin_res *QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP =
0373       GetResolution("pi-", seta, "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
0374   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name = "5x5 cluster";
0375 
0376   lin_res *QAG4SimJet_AntiKt_Tower_r02_Leading_Et =
0377       GetResolution("pi-", seta, "h_QAG4SimJet_AntiKt_Tower_r02_Leading_Et");
0378   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name = "Anti-k_{T} R=0.2";
0379 
0380   lin_res *QAG4SimJet_AntiKt_Tower_r04_Leading_Et =
0381       GetResolution("pi-", seta, "h_QAG4SimJet_AntiKt_Tower_r04_Leading_Et");
0382   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name = "Anti-k_{T} R=0.4";
0383 
0384   // ------------
0385   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetMarkerStyle(kFullCircle);
0386   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetMarkerColor(kRed + 3);
0387   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetLineColor(kRed + 3);
0388 
0389   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetMarkerStyle(kFullCircle);
0390   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetMarkerColor(kRed + 3);
0391   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetLineColor(kRed + 3);
0392 
0393   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->SetLineColor(kRed + 3);
0394 
0395   // ------------
0396   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetMarkerStyle(kFullSquare);
0397   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetMarkerColor(kBlue + 3);
0398   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetLineColor(kBlue + 3);
0399 
0400   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetMarkerStyle(kFullSquare);
0401   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetMarkerColor(kBlue + 3);
0402   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetLineColor(kBlue + 3);
0403 
0404   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->SetLineColor(kBlue + 3);
0405 
0406   // ------------
0407   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetMarkerStyle(kFullCross);
0408   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetMarkerColor(kMagenta + 3);
0409   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetLineColor(kMagenta + 3);
0410 
0411   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetMarkerStyle(kFullCross);
0412   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetMarkerColor(kMagenta + 3);
0413   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetLineColor(kMagenta + 3);
0414 
0415   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->SetLineColor(kMagenta + 3);
0416   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->SetLineStyle(kDashed);
0417 
0418   // ------------
0419   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetMarkerStyle(kStar);
0420   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetMarkerColor(kCyan + 3);
0421   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetLineColor(kCyan + 3);
0422 
0423   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetMarkerStyle(kStar);
0424   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetMarkerColor(kCyan + 3);
0425   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetLineColor(kCyan + 3);
0426 
0427   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->SetLineColor(kCyan + 3);
0428   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->SetLineStyle(kDashed);
0429 
0430   TCanvas *c1 = new TCanvas(Form("ClusterSizeScan_") + seta,
0431                             Form("ClusterSizeScan_") + seta, 1300, 600);
0432   c1->Divide(2, 1);
0433   int idx = 1;
0434   TPad *p;
0435 
0436   p = (TPad *) c1->cd(idx++);
0437   c1->Update();
0438   //  p->SetLogy();
0439 
0440   p->DrawFrame(0, 0, 60, 1.2,
0441                Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
0442   TLine *l = new TLine(0, 1, 60, 1);
0443   l->SetLineWidth(2);
0444   l->SetLineColor(kGray);
0445   l->Draw();
0446 
0447   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->Draw("ep");
0448   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->Draw("ep");
0449   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->Draw("ep");
0450   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->Draw("ep");
0451 
0452   TLegend *leg = new TLegend(.2, .2, .85, .45);
0453   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0454   leg->AddEntry("", "", "");
0455   leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity, QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name, "ep");
0456   leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity, QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name, "ep");
0457   leg->AddEntry(QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity, QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name, "ep");
0458   leg->AddEntry(QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity, QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name, "ep");
0459   leg->Draw();
0460 
0461   p = (TPad *) c1->cd(idx++);
0462   c1->Update();
0463 
0464   TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
0465                               32);
0466   f_calo_t1044->SetParameters(13.5, 64.9);
0467   f_calo_t1044->SetLineWidth(3);
0468   f_calo_t1044->SetLineColor(kGreen + 2);
0469 
0470   TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
0471                              Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
0472 
0473   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->Draw("same");
0474   QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->Draw("ep");
0475 
0476   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->Draw("same");
0477   QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->Draw("ep");
0478 
0479   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->Draw("same");
0480   QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->Draw("ep");
0481 
0482   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->Draw("same");
0483   QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->Draw("ep");
0484 
0485   f_calo_t1044->Draw("same");
0486 
0487   TLegend *leg = new TLegend(.2, .6, .85, .9);
0488   leg->SetHeader("#splitline{#it{#bf{sPHENIX}}  Geant4 Simulation}{" + config + "}");
0489   leg->AddEntry("", "", "");
0490   leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity,
0491                 QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name + "      " +
0492                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0493                          QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->GetParameter(0),
0494                          QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->GetParameter(1)),
0495                 "lpe");
0496   leg->AddEntry(f_calo_t1044,
0497                 TString("T-1044 (4x4, #eta=0)") + " " +
0498                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0499                          f_calo_t1044->GetParameter(0),
0500                          f_calo_t1044->GetParameter(1)),
0501                 "l");
0502   leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity,
0503                 QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name + "      " +
0504                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0505                          QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->GetParameter(0),
0506                          QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->GetParameter(1)),
0507                 "lpe");
0508   leg->AddEntry(QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity,
0509                 QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name + "  " +
0510                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0511                          QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->GetParameter(0),
0512                          QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->GetParameter(1)),
0513                 "lpe");
0514   leg->AddEntry(QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity,
0515                 QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name + "  " +
0516                     Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
0517                          QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->GetParameter(0),
0518                          QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->GetParameter(1)),
0519                 "lpe");
0520   leg->Draw();
0521 
0522   SaveCanvas(c1,
0523              base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
0524 }
0525 
0526 lin_res *
0527 GetResolution(TString PID = "pi-", TString eta = "eta0", TString qa_histo = "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP")
0528 {
0529   TString dataset_name = PID + "_" + eta + "_" + qa_histo;
0530 
0531   const int N = 6;
0532   const double Ps[] = {8, 12, 16, 32, 40, 50, -1, -1};
0533 
0534   vector<double> mean;
0535   vector<double> mean_err;
0536   vector<double> res;
0537   vector<double> res_err;
0538 
0539   TCanvas *c1 = new TCanvas("GetResolution_LineShape_" + dataset_name,
0540                             "GetResolution_LineShape_" + dataset_name, 1800, 900);
0541 
0542   c1->Divide(3, 2);
0543   int idx = 1;
0544   TPad *p;
0545 
0546  // for (int i = 0; i < N; ++i)
0547   for (int i = 2; i < 3; ++i) //for 16GeV
0548   {
0549   // if(i==2){ //for 16GeV
0550 
0551     p = (TPad *) c1->cd(idx++);
0552     c1->Update();
0553 
0554     const double P = Ps[i];
0555 
0556     TString sEnergy = Form("%.0fGeV", P);
0557 
0558     TH1 *h = NULL;
0559 
0560     //if (qa_histo.Length() == 0)
0561       h = GetEvalHisto(PID, eta, sEnergy);
0562   //  else
0563    //   h = GetQAHisto(PID, eta, sEnergy, qa_histo);
0564 
0565     TF1 *fgaus_g = new TF1("fgaus_LG_g_" + dataset_name, "gaus",
0566                            h->GetMean() - 4 * h->GetRMS(), h->GetMean() + 4 * h->GetRMS());
0567     fgaus_g->SetParameters(1, h->GetMean(),
0568                            h->GetRMS());
0569     h->Fit(fgaus_g, "MR0N");
0570     fgaus_g->SetLineColor(kGray);
0571 
0572     TF1 *fgaus = new TF1("fgaus_LG_" + dataset_name, "gaus",
0573                          fgaus_g->GetParameter(1) - 1.5 * fgaus_g->GetParameter(2),
0574                          fgaus_g->GetParameter(1) + 1.5 * fgaus_g->GetParameter(2));
0575     fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0576                          fgaus_g->GetParameter(2));
0577     fgaus->SetLineColor(kBlue);
0578     h->Fit(fgaus, "MR0N");
0579 
0580     h->Draw();
0581     fgaus_g->Draw("same");
0582     fgaus->Draw("same");
0583 
0584     double scale = 1;
0585     if (qa_histo.Contains("Leading_Et")) scale = 1. / P;
0586 
0587     mean.push_back(fgaus->GetParameter(1) * scale);
0588     mean_err.push_back(fgaus->GetParError(1) * scale);
0589     res.push_back(fgaus->GetParameter(2) / fgaus->GetParameter(1));
0590     res_err.push_back(fgaus->GetParError(2) / fgaus->GetParameter(1));
0591 
0592 
0593   //}// for 16GeV
0594   }
0595 
0596   TGraphErrors *ge_linear = new TGraphErrors(N, Ps,
0597                                              &mean[0], 0, &mean_err[0]);
0598 
0599   TGraphErrors *ge_res = new TGraphErrors(N, Ps,
0600                                           &res[0], 0, &res_err[0]);
0601   ge_res->GetHistogram()->SetStats(0);
0602   ge_res->Print();
0603 
0604   lin_res *ret = new lin_res;
0605   ret->name = dataset_name;
0606   ret->linearity = ge_linear;
0607   ret->resolution = ge_res;
0608   ret->f_res = new TF1("f_calo_r_" + dataset_name, "sqrt([0]*[0]+[1]*[1]/x)/100",
0609                        6, 60);
0610   ge_res->Fit(ret->f_res, "RM0N");
0611 
0612   SaveCanvas(c1,
0613              base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
0614 
0615   new TCanvas;
0616 
0617   ret->f_res->Draw();
0618   ge_res->Draw("*e");
0619 
0620   return ret;
0621 }
0622 
0623 TH1 *GetEvalHisto(TString PID = "pi-", TString eta = "eta0", TString energy = "16GeV")
0624 {
0625   TString dataset_name = PID + "_" + eta + "_" + energy;
0626   TH1 *h_E = new TH1F("h_ESumOverp_" + dataset_name, dataset_name + " h_ESumOverp;E/p;Count per bin", 100, 0, 2);
0627   //TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + "_*.root_g4svtx_eval.root";
0628   TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + ".root_g4svtx_eval.root";
0629 
0630   TChain *ntp_track = new TChain("ntp_track");
0631   const int n = ntp_track->Add(infile);
0632 
0633   cout << "GetEvalHisto - loaded " << n << " files from " << infile << endl;
0634 
0635   ntp_track->Draw("(cemce3x3 + hcaloute3x3 + hcaline3x3)/sqrt(gpx**2+gpy**2+gpz**2)>>h_ESumOverp_" + dataset_name,
0636                   "nfromtruth>50 && gtrackID == 1", "goff");
0637 
0638   delete ntp_track;
0639 
0640   return h_E;
0641 }
0642 
0643 TH1 *GetQAHisto(TString PID = "pi-", TString eta = "eta0", TString energy = "16GeV", TString qa_histo = "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP")
0644 {
0645   cout << "Weihu Ma run test " << endl; 
0646  
0647   TString dataset_name = PID + "_" + eta + "_" + energy;
0648   TH1 *h = NULL;
0649   //TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + "_*.root_qa.root";
0650   TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + ".root_qa.root";
0651   cout << "GetQAHisto - searching " << infile << endl;
0652 
0653   TChain *T = new TChain("T");
0654   const int n = T->Add(infile);
0655 
0656   TObjArray *fileElements = T->GetListOfFiles();
0657   TIter next(fileElements);
0658   TChainElement *chEl = 0;
0659   while ((chEl = (TChainElement *) next()))
0660   {
0661     cout << "GetQAHisto - processing " << chEl->GetTitle() << endl;
0662 
0663     TFile f(chEl->GetTitle());
0664 
0665     if (f.IsOpen())
0666     {
0667       TH1 *hqa = (TH1 *) f.GetObjectChecked(qa_histo, "TH1");
0668       cout << "Weihu Ma run test1 " << endl; 
0669       if (hqa)
0670       {
0671         cout << "Weihu Ma run test4 " << endl; 
0672         if (h)
0673         {
0674           cout << "Weihu Ma run test2 " << endl; 
0675           h->Add(hqa);
0676         }
0677         else
0678         {
0679           cout << "Weihu Ma run test3 " << endl; 
0680           h = (TH1 *) (hqa->Clone(TString(hqa->GetName()) + "_" + dataset_name));
0681           h->SetDirectory(NULL);
0682           assert(h);
0683         }
0684       }
0685     }
0686   }
0687   //assert(h);
0688 
0689   return h;
0690 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!