macros/tpc2019/DrawTpcPrototypeGenFitTrkFitter.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 TFile *_file0 = NULL;
0134 TString description;
0135 TTree *T(nullptr);
0136
0137 void Resolution(const TCut &cut = "Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=12 && TPCTrack.clusterSizePhi > 3.5",
0138                 const double phi_start = -2.905, const double phi_end = -2.885
0139
0140 )
0141 {
0142   TH2 *hresidual_phi = new TH2F("hresidual_phi", "hresidual_phi", 60, phi_start, phi_end, 60, -.2, .2);
0143   T->Draw("TPCTrack.clusterResidualPhi:TPCTrack.clusterProjectionPhi>>hresidual_phi",
0144           cut, "goff");
0145   hresidual_phi->SetTitle(";Global Phi [rad];Phi Residual [cm]");
0146
0147   TH2 *hresidual_phi_highres = new TH2F("hresidual_phi_highres", "hresidual_phi", 500, phi_start, phi_end, 1000, -.2, .2);
0148   T->Draw("TPCTrack.clusterResidualPhi:TPCTrack.clusterProjectionPhi>>hresidual_phi_highres",
0149           cut, "goff");
0150   hresidual_phi_highres->SetTitle(";Global Phi [rad];Phi Residual [cm]");
0151
0152   TCanvas *c1 = new TCanvas("Resolution", "Resolution", 1800, 860);
0153   c1->Divide(3, 1);
0154   int idx = 1;
0155   TPad *p = nullptr;
0156
0157   p = (TPad *) c1->cd(idx++);
0158   c1->Update();
0159
0160   TGraphErrors *gcenter = FitProfile(hresidual_phi);
0161
0162   hresidual_phi->Draw("colz");
0163   gcenter->Draw("p");
0164   //  hresidual_phi
0165
0166   //  TF1 *fPeriod = new TF1("fPeriod",
0167   //                         "[0] * sin(2*pi*x/(2*pi/12/8/16)) + [1] * sin(2*pi*x/(2*pi/12/8/16)) + [2] + [3]*x + [4]*x*x + [5] * sin(4*pi*x/(2*pi/12/8/16)) + [6] * sin(4*pi*x/(2*pi/12/8/16))+ [7] * sin(3*pi*x/(2*pi/12/8/16)) + [8] * sin(3*pi*x/(2*pi/12/8/16))+ [9] * sin(pi*x/(5*pi/12/8/16)) + [10] * sin(pi*x/(5*pi/12/8/16)) + [11] * x*x*x", -3, 0);
0168   TF1 *fPeriod = new TF1("fPeriod",
0169                          "([0]+[1]*x) * sin(2*pi*x/(2*pi/12/8/16)) + ([1]+[2]*x) * cos(2*pi*x/(2*pi/12/8/16)) + [3] + [4]*x + [5]*x*x + ([6]+[7]*x) * sin(4*pi*x/(2*pi/12/8/16)) + ([8]+[9]*x) * cos(4*pi*x/(2*pi/12/8/16))+ ([10] + [11]*x) * sin(3*pi*x/(2*pi/12/8/16)) + ([12]+[13]*x) * cos(3*pi*x/(2*pi/12/8/16))+ ([14]+[15]*x) * sin(pi*x/(5*pi/12/8/16)) + ([16]+[17]*x) * cos(pi*x/(5*pi/12/8/16)) + ([18]+[19]*x) * x*x*x", -3, 0);
0170   gcenter->Fit(fPeriod, "M");
0171
0172   p = (TPad *) c1->cd(idx++);
0173   c1->Update();
0174
0175   TH2 *hresidual_phi_cor = (TH2 *) hresidual_phi->Clone("hresidual_phi_cor");
0176   hresidual_phi_cor->Reset();
0177   hresidual_phi_cor->SetTitle(";Global Phi [rad];Modulation-Corrected Phi Residual [cm]");
0178
0179   double sum = 0;
0180   for (int binx = 1; binx <= hresidual_phi_highres->GetNbinsX(); ++binx)
0181     for (int biny = 1; biny <= hresidual_phi_highres->GetNbinsY(); ++biny)
0182     {
0183       const double x = hresidual_phi_highres->GetXaxis()->GetBinCenter(binx);
0184       const double y = hresidual_phi_highres->GetYaxis()->GetBinCenter(biny);
0185       const double value = hresidual_phi_highres->GetBinContent(binx, biny);
0186
0187       sum += value;
0188       //      const double y_cor = y;
0189       const double y_cor = y - fPeriod->Eval(x);
0190
0191       hresidual_phi_cor->Fill(x, y_cor, value);
0192     }
0193   cout << __PRETTY_FUNCTION__ << " sum = " << sum << endl;
0194
0195   hresidual_phi_cor->Draw("colz");
0196
0197   p = (TPad *) c1->cd(idx++);
0198   c1->Update();
0199   TH1 *hresidual_cor = hresidual_phi_cor->ProjectionY("hresidual_cor", 1, hresidual_phi_cor->GetNbinsX());
0200
0201   hresidual_cor->Draw();
0202
0203   TH1 *hresidual = hresidual_phi->ProjectionY("hresidual", 1, hresidual_phi->GetNbinsX());
0204   //  hresidual->SetMarkerColor(kRed - 3);
0205   hresidual->SetLineColor(kRed - 3);
0206   hresidual->Draw("same");
0207
0208   TF1 *hresidual_fgaus = new TF1("hresidual_fgaus", "gaus", -.2, .2);
0209   //
0210   TF1 *hresidual_dgaus = new TF1(Form("hresidual_dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0211                                  -.2, .2);
0212   //
0213   hresidual_fgaus->SetParameter(1, hresidual->GetMean());
0214   hresidual_fgaus->SetParameter(2, hresidual->GetRMS());
0215   //
0216   hresidual_cor->Fit(hresidual_fgaus, "MQ");
0217   //
0218   //  TF1 f2(Form("dgaus"), "gaus  + [3]",
0219   //         -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
0220   //
0221   hresidual_dgaus->SetParameters(hresidual_fgaus->GetParameter(0),
0222                                  hresidual_fgaus->GetParameter(1),
0223                                  hresidual_fgaus->GetParameter(2) / 2,
0224                                  hresidual_fgaus->GetParameter(0) / 20,
0225                                  hresidual_fgaus->GetParameter(2) * 10, 0);
0226   //
0227   hresidual_cor->Fit(hresidual_dgaus, "M");
0228   hresidual_dgaus->SetLineColor(kBlue + 2);
0229   const double resolution = TMath::Min(hresidual_dgaus->GetParameter(2), hresidual_dgaus->GetParameter(4));
0230
0231   cout << __PRETTY_FUNCTION__ << " hresidual_fgaus->GetParameter(2) = " << hresidual_dgaus->GetParameter(2)
0232        << "  hresidual_fgaus->GetParameter(4) = " << hresidual_dgaus->GetParameter(4) << " resolution = " << resolution << endl;
0233   //
0234   //  //      new TCanvas;
0235   //  //      h1->Draw();
0236   //  //      fgaus.Draw("same");
0237   //  //      break;
0238   //
0239   //  x[n] = p2->GetBinCenter(i);
0240   //  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0241   //  y[n] = fgaus.GetParameter(1);
0242   //  ey[n] = fgaus.GetParError(1);
0243
0244   gPad->SetTopMargin(.3);
0245   TLegend *leg = new TLegend(.03, .7, .98, .99, description + ": 1-removed residual");
0246   leg->AddEntry(hresidual, "Raw residual", "l");
0247   leg->AddEntry(hresidual_cor, "Position dependence corected residual", "p");
0248   leg->AddEntry(hresidual_dgaus, Form("resolution = %.0f #mum", resolution * 1e4), "l");
0249   leg->Draw();
0250
0251   SaveCanvas(c1,
0252              TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
0253 }
0254
0255 void TrackQA()
0256 {
0257   TCanvas *c1 = new TCanvas("TrackQA", "TrackQA", 1800, 860);
0258   c1->Divide(3, 2);
0259   int idx = 1;
0260   TPad *p = nullptr;
0261
0262   p = (TPad *) c1->cd(idx++);
0263   c1->Update();
0264
0265   TH1 *hnTrack = new TH1F("hnTrack", ";# track per event", 17, -.5, 16.5);
0266   T->Draw("nTrack>>hnTrack");
0267
0268   p = (TPad *) c1->cd(idx++);
0269   c1->Update();
0270
0271   TH1 *hCluster = new TH1F("hCluster", ";# cluster per track", 16, .5, 16.5);
0272   T->Draw("TPCTrack.nCluster>>hCluster");
0273
0274   p = (TPad *) c1->cd(idx++);
0275   c1->Update();
0276
0277   TH1 *hclusterSizePhi = new TH1F("hclusterSizePhi", ";Cluster phi size", 16, .5, 16.5);
0278   T->Draw("TPCTrack.clusterSizePhi>>hclusterSizePhi");
0279
0280   p = (TPad *) c1->cd(idx++);
0281   c1->Update();
0282
0283   TH1 *hclusterProjectionPhi = new TH1F("hclusterProjectionPhi", ";Cluster phi [rad]", 100, -3.5, -2.5);
0284   T->Draw("TPCTrack.clusterProjectionPhi>>hclusterProjectionPhi");
0285
0286   p = (TPad *) c1->cd(idx++);
0287   c1->Update();
0288
0289   TH1 *hAngle = new TH1F("hAngle", ";Horizontal angle [degree]", 100, -30, 30);
0290   T->Draw("atan(TPCTrack.pz/ TPCTrack.px)/pi*180>>hAngle");
0291
0292   p = (TPad *) c1->cd(idx++);
0293   c1->Update();
0294
0295   TH1 *hAngleV = new TH1F("hAngleV", ";Vertical angle [degree]", 100, -30, 30);
0296   T->Draw("(atan(TPCTrack.py/ TPCTrack.px) - 2*pi/12/2 )/pi*180>>hAngleV");
0297   //  TH1 *hresidualRough = new TH1F("hresidualRough", ";Rought phi residual [cm]", 1000, -1, 1);
0298   //  T->Draw("TPCTrack.clusterResidualPhi>>hresidualRough", "Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=10");
0299
0300   SaveCanvas(c1,
0301              TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
0302 }
0303
0304 void Track3D()
0305 {
0306   TCanvas *c1 = new TCanvas("Track3D", "Track3D", 1800, 900);
0307   c1->Divide(2, 1);
0308   int idx = 1;
0309   TPad *p = nullptr;
0310
0311   p = (TPad *) c1->cd(idx++);
0312   c1->Update();
0313   //  p->SetLogx();
0314   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0315   //
0316   TH3F *h3ClusterOverlay = new TH3F("h3ClusterOverlay", "h3ClusterOverlay", 128, -40, 10, 64, 35, 65, 128, -15, 15);
0317   //
0318   T->SetAlias("PhiCenter", "pi/12 + pi");  // center line azimuthal angle for TPC sector 0
0319   T->Draw("TPCTrack.clusterX*cos(PhiCenter + pi/2) + TPCTrack.clusterY*sin(PhiCenter + pi/2):TPCTrack.clusterX*cos(PhiCenter) + TPCTrack.clusterY*sin(PhiCenter):TPCTrack.clusterZ>>h3ClusterOverlay", "", "BOX2");
0320   //               "TPCTrack.clusterE", "BOX2");
0321   h3ClusterOverlay->SetTitle(";Drift Direction [cm];Pad Row Direction [cm];Azimuth Direction [cm]");
0322   h3ClusterOverlay->SetLineWidth(0);
0323   //
0324   p = (TPad *) c1->cd(idx++);
0325   c1->Update();
0326   //  p->SetLogx();
0327   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0328
0329   TH1 *h2ClusterOverlay = h3ClusterOverlay->Project3D("zx");
0330   h2ClusterOverlay->Draw("colz");
0331   //    eventT->Draw("Clusters.avg_pady:Clusters.min_sample+Clusters.peak_sample>>h2ClusterOverlay",
0332   //                 "Clusters.peak", "colz");
0333   //    h2ClusterOverlay->SetTitle(";Time [0-127*50ns];Pads [0-127]");
0334   //  h2ClusterOverlay->SetLineWidth(0);
0335
0336   p->SetTopMargin(.9);
0337   TLegend *leg = new TLegend(.15, .9, .95, .99, description + ": accumulated clusters on tracks");
0338   leg->Draw();
0339
0340   SaveCanvas(c1,
0341              TString(_file0->GetName()) + TString(c1->GetName()), false);
0342 }
0343
0344 void TrackDistortion(const TCut &cut = "TPCTrack.nCluster>=10")
0345 {
0346   TCanvas *c1 = new TCanvas("TrackDistortion", "TrackDistortion", 1800, 860);
0347   c1->Divide(2, 1);
0348   int idx = 1;
0349   TPad *p = nullptr;
0350
0351   p = (TPad *) c1->cd(idx++);
0352   c1->Update();
0353
0354   TH2 *hPhiDistortion = new TH2F("hPhiDistortion", ";Pad Layers;Phi Residual [cm]", 16, -.5, 15.5, 200, -3, 3);
0355   T->Draw("TPCTrack.clusterResidualPhi:Iteration$>>hPhiDistortion", cut, "goff");
0356   TGraph *gPhiDistortion = FitProfile(hPhiDistortion);
0357   hPhiDistortion->Draw("colz");
0358   gPhiDistortion->Draw("p");
0359
0360   p = (TPad *) c1->cd(idx++);
0361   c1->Update();
0362
0363   TH2 *hZDistortion = new TH2F("hZDistortion", ";Pad Layers;Z Residual [cm]", 16, -.5, 15.5, 200, -3, 3);
0364   T->Draw("TPCTrack.clusterResidualZ:Iteration$>>hZDistortion", cut, "goff");
0365   TGraph *gZDistortion = FitProfile(hZDistortion);
0366   hZDistortion->Draw("colz");
0367   gZDistortion->Draw("p");
0368
0369   SaveCanvas(c1,
0370              TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
0371 }
0372
0373 Double_t langaufun(Double_t *x, Double_t *par)
0374 {
0375   //Fit parameters:
0376   //par[0]=Width (scale) parameter of Landau density
0377   //par[1]=Most Probable (MP, location) parameter of Landau density
0378   //par[2]=Total area (integral -inf to inf, normalization constant)
0379   //par[3]=Width (sigma) of convoluted Gaussian function
0380   //
0381   //In the Landau distribution (represented by the CERNLIB approximation),
0382   //the maximum is located at x=-0.22278298 with the location parameter=0.
0383   //This shift is corrected within this function, so that the actual
0384   //maximum is identical to the MP parameter.
0385
0386   // Numeric constants
0387   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0388   Double_t mpshift = -0.22278298;       // Landau maximum location
0389
0390   // Control constants
0391   Double_t np = 100.0;  // number of convolution steps
0392   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0393
0394   // Variables
0395   Double_t xx;
0396   Double_t mpc;
0397   Double_t fland;
0398   Double_t sum = 0.0;
0399   Double_t xlow, xupp;
0400   Double_t step;
0401   Double_t i;
0402
0403   // MP shift correction
0404   mpc = par[1] - mpshift * par[0];
0405
0406   // Range of convolution integral
0407   xlow = x[0] - sc * par[3];
0408   xupp = x[0] + sc * par[3];
0409
0410   step = (xupp - xlow) / np;
0411
0412   // Convolution integral of Landau and Gaussian by sum
0413   for (i = 1.0; i <= np / 2; i++)
0414   {
0415     xx = xlow + (i - .5) * step;
0416     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0417     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0418
0419     xx = xupp - (i - .5) * step;
0420     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0421     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0422   }
0423
0424   return (par[2] * step * sum * invsq2pi / par[3]);
0425 }
0426
0427 TF1 *langaufit(TH1 *his,
0428                const vector<double> &fitrange,
0429                const vector<double> &startvalues,
0430                const vector<double> &parlimitslo, const vector<double> &parlimitshi
0431                //    ,
0432                //    const vector<double> &fitparams, const vector<double> &fiterrors,
0433                //    const vector<double> &ChiSqr, Int_t *NDF
0434 )
0435 {
0436   // Once again, here are the Landau * Gaussian parameters:
0437   //   par[0]=Width (scale) parameter of Landau density
0438   //   par[1]=Most Probable (MP, location) parameter of Landau density
0439   //   par[2]=Total area (integral -inf to inf, normalization constant)
0440   //   par[3]=Width (sigma) of convoluted Gaussian function
0441   //
0442   // Variables for langaufit call:
0443   //   his             histogram to fit
0444   //   fitrange[2]     lo and hi boundaries of fit range
0445   //   startvalues[4]  reasonable start values for the fit
0446   //   parlimitslo[4]  lower parameter limits
0447   //   parlimitshi[4]  upper parameter limits
0448   //   fitparams[4]    returns the final fit parameters
0449   //   fiterrors[4]    returns the final fit errors
0450   //   ChiSqr          returns the chi square
0451   //   NDF             returns ndf
0452
0453   Int_t i;
0454   Char_t FunName[100];
0455
0456   sprintf(FunName, "Fitfcn_%s", his->GetName());
0457
0458   TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(FunName);
0459   if (ffitold) delete ffitold;
0460
0461   TF1 *ffit = new TF1(FunName, langaufun, fitrange[0], fitrange[1], 4);
0462   ffit->SetParameters(startvalues.data());
0463   ffit->SetParNames("Width", "MP", "Area", "GSigma");
0464
0465   for (i = 0; i < 4; i++)
0466   {
0467     ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0468   }
0469
0470   his->Fit(FunName, "RB0");  // fit within specified range, use ParLimits, do not plot
0471
0472   //   ffit->GetParameters(fitparams);    // obtain fit parameters
0473   //   for (i=0; i<4; i++) {
0474   //      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
0475   //   }
0476   //   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
0477   //   NDF[0] = ffit->GetNDF();           // obtain ndf
0478
0479   return (ffit);  // return fit function
0480 }
0481
0482 void TrackClusterEnergy(const TCut &cut = "TPCTrack.nCluster>=10")
0483 {
0484   TCanvas *c1 = new TCanvas("TrackClusterEnergy", "TrackClusterEnergy", 1200, 860);
0485   //  c1->Divide(2, 1);
0486   int idx = 1;
0487   TPad *p = nullptr;
0488
0489   p = (TPad *) c1->cd(idx++);
0490   c1->Update();
0491   p->SetLogy();
0492
0493   TH1 *hClusterEnergy = new TH1F("hClusterEnergy", ";Cluster Energy on good track [ADU];Count / bin", 2500, 0, 5000);
0494   T->Draw("TPCTrack.clusterE>>hClusterEnergy", cut, "goff");
0495   //  TGraph *gPhiDistortion = FitProfile(hPhiDistortion);
0496
0497   TF1 *fit = langaufit(hClusterEnergy,
0498                        {0, 5000},
0499                        {200, 500, (double) T->GetEntries(cut) * 10, 200},
0500                        {0, 0, 0, 0},
0501                        {5000, 5000, 1e10, 5000});
0502   fit->SetLineColor(kBlue + 2);
0503
0504   hClusterEnergy->Draw();
0505   fit->Draw("same");
0506   //    gPhiDistortion->Draw("p");
0507
0508   TLegend *leg = new TLegend(.3, .7, .95, .9, description + ": Cluster Energy on good track");
0509   leg->AddEntry(hClusterEnergy, TString("Data: ") + cut.GetTitle(), "l");
0510   leg->AddEntry(fit,
0511                 Form("Langdau * Gauss Fit, #mu= %.0f ADU", fit->GetParameter(1)), "l");
0512   leg->Draw();
0513
0514   SaveCanvas(c1,
0515              TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
0516 }
0517
0518 void DrawTpcPrototypeGenFitTrkFitter(
0519     const TString infile = "data/tpc_beam/tpc_beam_00000292-0000.evt_TpcPrototypeGenFitTrkFitter.root", const TString desc = "Run 292"  //
0520 )
0521 {
0522   gSystem->Load("libtpc2019.so");
0523   //
0524   SetsPhenixStyle();
0525   TVirtualFitter::SetDefaultFitter("Minuit2");
0526   gStyle->SetLegendTextSize(0);
0527   //
0528   if (!_file0)
0529   {
0530     TString chian_str = infile;
0531     chian_str.ReplaceAll("ALL", "*");
0532
0533     TChain *t = new TChain("T");
0534     const int n = t->Add(chian_str);
0535
0536     cout << "Loaded " << n << " root files with eventT in " << chian_str << endl;
0537     assert(n > 0);
0538
0539     T = t;
0540
0541     _file0 = new TFile;
0542     _file0->SetName(infile);
0543   }
0544   //
0545   description = desc;
0546
0547   TrackQA();
0548   TrackDistortion();
0549   Track3D();
0550   TrackClusterEnergy();
0551   //  Resolution();
0552   Resolution("Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=12 && TPCTrack.clusterSizePhi > 7");
0553 }