macros/tpc2019/DrawTpcPrototypeGenFitTrkFitter_Summary.C

0001 #include "SaveCanvas.C"
0002 #include "sPhenixStyle.C"
0003
0004 #include <TChain.h>
0005 #include <TCut.h>
0006 #include <TEfficiency.h>
0007 #include <TF1.h>
0008 #include <TGraphAsymmErrors.h>
0009 #include <TGraphErrors.h>
0010 #include <TH2.h>
0011 #include <TH3.h>
0012 #include <TPolyLine.h>
0013 #include <TTree.h>
0014
0015 #include <TFile.h>
0016
0017 #include <TColor.h>
0018 #include <TLatex.h>
0019 #include <TLegend.h>
0020 #include <TLine.h>
0021 #include <TStyle.h>
0022
0023 #include <TDirectory.h>
0024 #include <TMath.h>
0025 #include <TPad.h>
0026 #include <TString.h>
0027 #include <TTree.h>
0028 #include <TVirtualFitter.h>
0029
0030 #include <cmath>
0031 #include <iostream>
0032
0033 using namespace std;
0034
0035 // ROOT6 disabled assert. Well....
0036 #ifdef assert
0037 #undef assert
0038 #endif
0039 #define assert(exp)                                                                             \
0040   {                                                                                             \
0041     if (!(exp))                                                                                 \
0042     {                                                                                           \
0043       cout << "Assert (" << #exp << ") failed at " << __FILE__ << " line " << __LINE__ << endl; \
0044       exit(1);                                                                                  \
0045     }                                                                                           \
0046   }
0047
0048 //! utility function to
0049 void useLogBins(TAxis *axis)
0050 {
0051   assert(axis);
0052   assert(axis->GetXmin() > 0);
0053   assert(axis->GetXmax() > 0);
0054
0055   const int bins = axis->GetNbins();
0056
0057   Axis_t from = log10(axis->GetXmin());
0058   Axis_t to = log10(axis->GetXmax());
0059   Axis_t width = (to - from) / bins;
0060   vector<Axis_t> new_bins(bins + 1);
0061
0062   for (int i = 0; i <= bins; i++)
0063   {
0064     new_bins[i] = TMath::Power(10, from + i * width);
0065   }
0066   axis->Set(bins, new_bins.data());
0067 }
0068
0069 TGraphErrors *
0070 FitProfile(const TH2 *h2)
0071 {
0072   TProfile *p2 = h2->ProfileX();
0073
0074   int n = 0;
0075   double x[1000];
0076   double ex[1000];
0077   double y[1000];
0078   double ey[1000];
0079
0080   for (int i = 1; i <= h2->GetNbinsX(); i++)
0081   {
0082     TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0083
0084     if (h1->GetSum() < 30)
0085     {
0086       cout << "FitProfile - ignore bin " << i << endl;
0087       continue;
0088     }
0089     else
0090     {
0091       //      cout << "FitProfile - fit bin " << i << endl;
0092     }
0093
0094     TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0095               p2->GetBinError(i) * 4);
0096
0097     //    TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0098     //           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0099
0100     fgaus.SetParameter(1, p2->GetBinContent(i));
0101     fgaus.SetParameter(2, p2->GetBinError(i));
0102
0103     h1->Fit(&fgaus, "MQ");
0104
0105     TF1 f2(Form("dgaus"), "gaus  + [3]",
0106            -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
0107
0108     f2.SetParameters(fgaus.GetParameter(0), fgaus.GetParameter(1),
0109                      fgaus.GetParameter(2), 0);
0110
0111     h1->Fit(&f2, "MQR");
0112
0113     //      new TCanvas;
0114     //      h1->Draw();
0115     //      fgaus.Draw("same");
0116     //      break;
0117
0118     x[n] = p2->GetBinCenter(i);
0119     ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0120     y[n] = fgaus.GetParameter(1);
0121     ey[n] = fgaus.GetParError(1);
0122
0123     //      p2->SetBinContent(i, fgaus.GetParameter(1));
0124     //      p2->SetBinError(i, fgaus.GetParameter(2));
0125
0126     n++;
0127     delete h1;
0128   }
0129
0130   return new TGraphErrors(n, x, y, ex, ey);
0131 }
0132
0133 string basePath;
0134 const static double x0 = 6 - (2 + 7. / 16);  // inch
0135
0136 pair<TGraph *, TGraph *> getResolutions(vector<pair<double, string>> datafiles, Color_t color = kRed + 1)
0137 {
0138   vector<double> x_cm;
0139   vector<double> x;
0140   vector<double> res;
0141   vector<double> res_err;
0142   vector<double> res2;
0143   vector<double> res2_err;
0144
0145   for (const auto &data : datafiles)
0146   {
0147     cout << __PRETTY_FUNCTION__ << " "
0148          << " processing x = " << data.first << " from " << data.second << endl;
0149
0150     TFile *f = TFile::Open((TString(basePath + data.second.c_str()) + "_TpcPrototypeGenFitTrkFitter.root_DrawJet_Resolution.root"));
0151     assert(f);
0152     assert(f->IsOpen());
0153     TH1 *hresidual_cor = (TH1 *) f->GetObjectChecked("hresidual_cor", "TH1");
0154     assert(hresidual_cor);
0155     TF1 *hresidual_dgaus = (TF1 *) hresidual_cor->GetListOfFunctions()->At(0);
0156     assert(hresidual_dgaus);
0157
0158     if (hresidual_dgaus->GetParameter(2) < hresidual_dgaus->GetParameter(4))
0159     {
0160       res.push_back(hresidual_dgaus->GetParameter(2) * 1e4);
0161       res_err.push_back(hresidual_dgaus->GetParError(2) * 1e4);
0162     }
0163     else
0164     {
0165       res.push_back(hresidual_dgaus->GetParameter(4) * 1e4);
0166       res_err.push_back(hresidual_dgaus->GetParError(4) * 1e4);
0167     }
0168     x.push_back(data.first);
0169     x_cm.push_back((data.first - x0) * 2.54);
0170     res2.push_back(res[res.size() - 1] * res[res.size() - 1]);
0171     res2_err.push_back(2 * res[res.size() - 1] * res_err[res.size() - 1]);
0172
0173     cout << "x = " << x[res.size() - 1] << " x_cm = " << x_cm[res.size() - 1] << " res = " << res[res.size() - 1] << " res_err = " << res_err[res.size() - 1] << endl;
0174   }
0175
0176   TGraphErrors *g = new TGraphErrors(x.size(), x_cm.data(), res.data(), 0, res_err.data());
0177   g->SetMarkerColor(color);
0178   g->SetLineColor(color);
0179   g->Print();
0180
0181   TGraphErrors *g2 = new TGraphErrors(x.size(), x.data(), res2.data(), 0, res2_err.data());
0182   g2->SetMarkerColor(color);
0183   g2->SetLineColor(color);
0184   g2->Print();
0185
0186   return pair<TGraph *, TGraph *>(g, g2);
0187 }
0188
0189 void DrawTpcPrototypeGenFitTrkFitter_Summary(  //
0190     string basepath = "/phenix/u/jinhuang/links/sPHENIX_work/TPC/fnal_June2019/SimPadPlaneIter5/")
0191 {
0192   //  gSystem->Load("libtpc2019.so");
0193   //
0194   basePath = basepath;
0195   const double max_res = 300;
0196   const double field_off_trans_diffusion = 150;  // um/sqrt(cm), From: Prakhar Garg <prakhar.garg@stonybrook.edu>  Date: Tue, Jun 18, 2019 at 1:04 PM
0197   const double field_on_trans_diffusion = 40;    // um/sqrt(cm), From: Prakhar Garg <prakhar.garg@stonybrook.edu>  Date: Tue, Jun 18, 2019 at 1:04 PM
0198
0199   SetsPhenixStyle();
0200   TVirtualFitter::SetDefaultFitter("Minuit2");
0201   gStyle->SetLegendTextSize(0);
0202
0203   TCanvas *c1 = new TCanvas("DrawTpcPrototypeGenFitTrkFitter_Summary", "DrawTpcPrototypeGenFitTrkFitter_Summary", 1800, 700);
0204   c1->Divide(2, 1);
0205   int idx = 1;
0206   TPad *p = nullptr;
0207
0208   p = (TPad *) c1->cd(idx++);
0209   c1->Update();
0210
0211   TH1 *frame = p->DrawFrame(2, 0, 20, max_res * max_res);
0212   frame->SetTitle(";Horizontal Position [in];Resolution^{2}, Cluster w/ 2+ pad [#mum^{2}]");
0213   frame->GetYaxis()->SetTitleOffset(1.75);
0214   TLine *l = new TLine(x0, 0, x0, max_res * max_res);
0215   l->Draw();
0216
0217   TLegend *leg1 = new TLegend(.35, .7, .9, .9, "2019 TPC Beam Test Preview");
0218
0219   p = (TPad *) c1->cd(idx++);
0220   c1->Update();
0221   p->DrawFrame(0, 50, 125, max_res)->SetTitle(";Drift Length, L [cm];Resolution, Cluster w/ 2+ pad [#mum]");
0222   TLegend *leg2 = new TLegend(.42, .45, .92, .93, "2019 TPC Beam Test Preview");
0223
0224   {
0225     Color_t color = kRed + 2;
0226     const TString name = "Scan2: Gap 300V, GEM 370V";
0227
0228     pair<TGraph *, TGraph *> scan2 = getResolutions(
0229         {{6, "tpc_beam_00000288-0000.evt"},
0230          {6, "tpc_beam_00000292-0000.evt"},
0231          {10, "tpc_beam_00000293-0000.evt"},
0232          {14, "tpc_beam_00000294-0000.evt"},
0233          {18, "tpc_beam_00000295-0000.evt"}},
0234         color);
0235
0236     c1->cd(1);
0237     scan2.second->Draw("p");
0238     TF1 *fline_scan2 = new TF1("fline_scan2", "pol1", x0, 20);
0239     scan2.second->Fit(fline_scan2, "MR0");
0240     fline_scan2->SetLineColor(color);
0241     fline_scan2->Draw("same");
0242
0243     leg1->AddEntry(scan2.second, name, "p");
0244
0245     c1->cd(2);
0246     scan2.first->Draw("p");
0247     TF1 *fdiff_scan2 = new TF1("fline_scan2", "sqrt([0]*[0] + ([1]*[1]/[2])*x)", 0, 40);
0248     fdiff_scan2->SetParameters(70, field_off_trans_diffusion, 20);
0249     fdiff_scan2->FixParameter(1, field_off_trans_diffusion);
0250     scan2.first->Fit(fdiff_scan2, "MR0");
0251     fdiff_scan2->SetLineColor(color);
0252     fdiff_scan2->Draw("same");
0253
0254     TF1 *fdiff_scan2_FieldON = (TF1 *) fdiff_scan2->Clone("fdiff_scan2_FieldON");
0255     fdiff_scan2_FieldON->SetRange(0, 105);
0256     fdiff_scan2_FieldON->SetParameter(1, field_on_trans_diffusion);
0257     fdiff_scan2_FieldON->SetLineStyle(kDashed);
0258     fdiff_scan2_FieldON->Draw("same");
0259     //    fdiff_scan2_FieldON->Print();
0260
0261     leg2->AddEntry(scan2.first, name, "pe");
0262     leg2->AddEntry(fdiff_scan2, Form("#sqrt{(%.1f #mum)^{2} + (%.0f #mum/#sqrt{cm} #times#sqrt{L/%.1f})^{2}}", abs(fdiff_scan2->GetParameter(0)), abs(fdiff_scan2->GetParameter(1)), abs(fdiff_scan2->GetParameter(2))), "l");
0263     leg2->AddEntry(fdiff_scan2_FieldON, Form("Field = 1.4T, #sigma_{T, Diffusion} = %.0f #mum/#sqrt{cm}", field_on_trans_diffusion), "l");
0264
0265     c1->Update();
0266   }
0267
0268   {
0269     Color_t color = kBlue + 2;
0270     const TString name = "Scan3: Gap 150V, GEM 380V";
0271
0272     pair<TGraph *, TGraph *> scan3 = getResolutions(
0273         {{6, "tpc_beam_00000297-0000.evt"},
0274          {6, "tpc_beam_00000298-0000.evt"},
0275          {10, "tpc_beam_00000299-0000.evt"},
0276          {14, "tpc_beam_00000300-0000.evt"},
0277          {18, "tpc_beam_00000301-0000.evt"}},
0278         color);
0279
0280     c1->cd(1);
0281     scan3.second->Draw("p");
0282     TF1 *fline_scan3 = new TF1("fline_scan3", "pol1", x0, 20);
0283     scan3.second->Fit(fline_scan3, "MR0");
0284     fline_scan3->SetLineColor(color);
0285     fline_scan3->Draw("same");
0286     leg1->AddEntry(scan3.second, name, "p");
0287
0288     c1->cd(2);
0289     scan3.first->Draw("p");
0290     TF1 *fdiff_scan3 = new TF1("fline_scan3", "sqrt([0]*[0] + ([1]*[1]/[2])*x)", 0, 40);
0291     fdiff_scan3->SetParameters(70, field_off_trans_diffusion, 20);
0292     fdiff_scan3->FixParameter(1, field_off_trans_diffusion);
0293     scan3.first->Fit(fdiff_scan3, "MR0");
0294     fdiff_scan3->SetLineColor(color);
0295     fdiff_scan3->Draw("same");
0296
0297     TF1 *fdiff_scan3_FieldON = (TF1 *) fdiff_scan3->Clone("fdiff_scan3_FieldON");
0298     fdiff_scan3_FieldON->SetRange(0, 105);
0299     fdiff_scan3_FieldON->SetParameter(1, field_on_trans_diffusion);
0300     fdiff_scan3_FieldON->SetLineStyle(kDashed);
0301     fdiff_scan3_FieldON->Draw("same");
0302     //    fdiff_scan3_FieldON->Print();
0303
0304     leg2->AddEntry(scan3.first, name, "pe");
0305     leg2->AddEntry(fdiff_scan3, Form("#sqrt{(%.1f #mum)^{2} + (%.0f #mum/#sqrt{cm} #times#sqrt{L/%.1f})^{2}}", abs(fdiff_scan3->GetParameter(0)), abs(fdiff_scan3->GetParameter(1)), abs(fdiff_scan3->GetParameter(2))), "l");
0306     leg2->AddEntry(fdiff_scan3_FieldON, Form("Field = 1.4T, #sigma_{T, Diffusion} = %.0f #mum/#sqrt{cm}", field_on_trans_diffusion), "l");
0307
0308     c1->Update();
0309   }
0310
0311   if (0)
0312   {
0313     Color_t color = kGreen + 2;
0314
0315     pair<TGraph *, TGraph *> scan4 = getResolutions(
0316         {{6, "tpc_beam_00000357-0000.evt"},
0317          {6, "tpc_beam_00000358-0000.evt"},
0318          {10, "tpc_beam_00000359-0000.evt"},
0319          {10, "tpc_beam_00000360-0000.evt"},
0320          {13, "tpc_beam_00000356-0000.evt"},
0321          {14, "tpc_beam_00000361-0000.evt"},
0322          {14, "tpc_beam_00000362-0000.evt"},
0323          {18, "tpc_beam_00000363-0000.evt"},
0324          {18, "tpc_beam_00000366-0000.evt"}},
0325         color);
0326
0327     c1->cd(1);
0328     scan4.second->Draw("p");
0329     TF1 *fline_scan4 = new TF1("fline_scan4", "pol1", x0, 16);
0330     scan4.second->Fit(fline_scan4, "MR0");
0331     fline_scan4->SetLineColor(color);
0332     fline_scan4->Draw("same");
0333
0334     c1->cd(2);
0335     scan4.first->Draw("p");
0336     TF1 *fdiff_scan4 = new TF1("fline_scan4", "sqrt([0]*[0] + [1]*[1]*x)", 0, 30);
0337     scan4.first->Fit(fdiff_scan4, "MR0");
0338     fdiff_scan4->SetLineColor(color);
0339     fdiff_scan4->Draw("same");
0340
0341     c1->Update();
0342   }
0343
0344   c1->cd(1);
0345   leg1->Draw();
0346   c1->cd(2);
0347   leg2->Draw();
0348   SaveCanvas(c1,
0349              TString(basePath) + TString(c1->GetName()), true);
0350 }