macros/tpc2019/DrawTpcPrototypeUnpacker.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 TFile *_file0 = NULL;
0049 TString description;
0050 TTree *eventT(nullptr);
0051 TTree *chanT(nullptr);
0052
0053 //! utility function to
0054 void useLogBins(TAxis *axis)
0055 {
0056   assert(axis);
0057   assert(axis->GetXmin() > 0);
0058   assert(axis->GetXmax() > 0);
0059
0060   const int bins = axis->GetNbins();
0061
0062   Axis_t from = log10(axis->GetXmin());
0063   Axis_t to = log10(axis->GetXmax());
0064   Axis_t width = (to - from) / bins;
0065   vector<Axis_t> new_bins(bins + 1);
0066
0067   for (int i = 0; i <= bins; i++)
0068   {
0069     new_bins[i] = TMath::Power(10, from + i * width);
0070   }
0071   axis->Set(bins, new_bins.data());
0072 }
0073
0074 TGraphErrors *
0075 FitProfile(const TH2 *h2)
0076 {
0077   TProfile *p2 = h2->ProfileX();
0078
0079   int n = 0;
0080   double x[1000];
0081   double ex[1000];
0082   double y[1000];
0083   double ey[1000];
0084
0085   for (int i = 1; i <= h2->GetNbinsX(); i++)
0086   {
0087     TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0088
0089     if (h1->GetSum() < 30)
0090     {
0091       cout << "FitProfile - ignore bin " << i << endl;
0092       continue;
0093     }
0094     else
0095     {
0096       //      cout << "FitProfile - fit bin " << i << endl;
0097     }
0098
0099     TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0100               p2->GetBinError(i) * 4);
0101
0102     //    TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0103     //           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0104
0105     fgaus.SetParameter(1, p2->GetBinContent(i));
0106     fgaus.SetParameter(2, p2->GetBinError(i));
0107
0108     h1->Fit(&fgaus, "MQ");
0109
0110     TF1 f2(Form("dgaus"), "gaus  + [3]",
0111            -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
0112
0113     f2.SetParameters(fgaus.GetParameter(0), fgaus.GetParameter(1),
0114                      fgaus.GetParameter(2), 0);
0115
0116     h1->Fit(&f2, "MQR");
0117
0118     //      new TCanvas;
0119     //      h1->Draw();
0120     //      fgaus.Draw("same");
0121     //      break;
0122
0123     x[n] = p2->GetBinCenter(i);
0124     ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0125     y[n] = fgaus.GetParameter(2);
0126     ey[n] = fgaus.GetParError(2);
0127
0128     //      p2->SetBinContent(i, fgaus.GetParameter(1));
0129     //      p2->SetBinError(i, fgaus.GetParameter(2));
0130
0131     n++;
0132     delete h1;
0133   }
0134
0135   return new TGraphErrors(n, x, y, ex, ey);
0136 }
0137
0138 void Clusters3D()
0139 {
0140   TCanvas *c1 = new TCanvas("Clusters3D", "Clusters3D", 1800, 900);
0141   c1->Divide(2, 1);
0142   int idx = 1;
0143   TPad *p = nullptr;
0144
0145   p = (TPad *) c1->cd(idx++);
0146   c1->Update();
0147   //  p->SetLogx();
0148   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0149
0150   TH3F *h3ClusterOverlay = new TH3F("h3ClusterOverlay", "h3ClusterOverlay", 128, -.5, 127.5, 16, -.5, 15.5, 128, -.5, 127.5);
0151
0152   eventT->Draw("Clusters.avg_pady:Clusters.avg_padx:Clusters.min_sample+Clusters.peak_sample>>h3ClusterOverlay",
0153                "Clusters.peak", "BOX2");
0154   h3ClusterOverlay->SetTitle(";Time [0-127*50ns];Rows [0-15];Pads [0-127]");
0155   h3ClusterOverlay->SetLineWidth(0);
0156
0157   p = (TPad *) c1->cd(idx++);
0158   c1->Update();
0159   //  p->SetLogx();
0160   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0161
0162   TH2 *h2ClusterOverlay = new TH2F("h2ClusterOverlay", "h2ClusterOverlay", 128, -.5, 127.5, 128, -.5, 127.5);
0163
0164   eventT->Draw("Clusters.avg_pady:Clusters.min_sample+Clusters.peak_sample>>h2ClusterOverlay",
0165                "Clusters.peak", "colz");
0166   h2ClusterOverlay->SetTitle(";Time [0-127*50ns];Pads [0-127]");
0167   //  h2ClusterOverlay->SetLineWidth(0);
0168
0169   p->SetTopMargin(.9);
0170   TLegend *leg = new TLegend(.05, .9, .95, .99, description + ": accumulated clusters");
0171   leg->Draw();
0172
0173   SaveCanvas(c1,
0174              TString(_file0->GetName()) + TString("_") + TString("Clusters3D") + TString("_") + TString(c1->GetName()), false);
0175 }
0176
0177 void ClusterQA()
0178 {
0179   TCanvas *c1 = new TCanvas("ClusterQA", "ClusterQA", 1800, 860);
0180   c1->Divide(3, 1);
0181   int idx = 1;
0182   TPad *p = nullptr;
0183
0184   p = (TPad *) c1->cd(idx++);
0185   c1->Update();
0186   //  p->SetLogx();
0187   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0188
0189   TH1 *CLusterEnergyAvg = new TH1F("CLusterEnergyAvg", ";<cluster energy> [ADU];Count", 100, -0, 2000);
0190   eventT->Draw("Sum$(Clusters[].peak)/nClusters>>CLusterEnergyAvg");
0191   CLusterEnergyAvg->Fit("gaus");
0192   TF1 *fgauss = (TF1 *) CLusterEnergyAvg->GetListOfFunctions()->At(0);
0193   assert(fgauss);
0194   fgauss->SetLineColor(kBlue + 1);
0195
0196   TLegend *leg = new TLegend(.2, .8, .99, .99, description);
0197   //    leg->AddEntry("", "Au+Au #sqrt{s_{NN}}=200GeV, 24B 0-10%C", "");
0198   leg->AddEntry(fgauss, Form("Energy = %.1f +/- %.1f", fgauss->GetParameter(1), fgauss->GetParError(1)), "l");
0199   leg->Draw();
0200
0201   p = (TPad *) c1->cd(idx++);
0202   c1->Update();
0203   //  p->SetLogx();
0204   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0205
0206   TH1 *CLusterEnergyAll = new TH1F("CLusterEnergyAll", ";cluster energy [ADU];Count", 100, -0, 2000);
0207   eventT->Draw("(Clusters[].peak)>>CLusterEnergyAll");
0208   //  CLusterEnergyAvg->Fit("gaus");
0209
0210   p = (TPad *) c1->cd(idx++);
0211   c1->Update();
0212   //  p->SetLogx();
0213   //  p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
0214
0215   TH1 *nClusters = new TH1F("nClusters", ";Cluster count;Count", 151, -.5, 150.5);
0216   eventT->Draw("nClusters>>nClusters");
0217   //  CLusterEnergyAvg->Fit("gaus");
0218
0219   SaveCanvas(c1,
0220              TString(_file0->GetName()) + TString("_") + TString("Clusters3D") + TString("_") + TString(c1->GetName()), false);
0221 }
0222
0223 void DrawTpcPrototypeUnpacker(
0224     const TString infile = "data/tpc_beam/tpc_beam_00000292-0000.evt_TpcPrototypeUnpacker.root", const TString desc = "Run 292"  //
0225 )
0226 {
0227   gSystem->Load("libtpc2019.so");
0228   //
0229   SetsPhenixStyle();
0230   TVirtualFitter::SetDefaultFitter("Minuit2");
0231   gStyle->SetLegendTextSize(0);
0232   //
0233   if (!_file0)
0234   {
0235     TString chian_str = infile;
0236     chian_str.ReplaceAll("ALL", "*");
0237
0238     TChain *t = new TChain("eventT");
0239     const int n = t->Add(chian_str);
0240
0241     cout << "Loaded " << n << " root files with eventT in " << chian_str << endl;
0242     assert(n > 0);
0243
0244     eventT = t;
0245
0246     t = new TChain("chanT");
0247     const int n1 = t->Add(chian_str);
0248
0249     cout << "Loaded " << n1 << " root files with chanT in " << chian_str << endl;
0250     assert(n1 > 0);
0251
0252     chanT = t;
0253
0254     _file0 = new TFile;
0255     _file0->SetName(infile);
0256   }
0257   //
0258   description = desc;
0259
0260   Clusters3D();
0261   ClusterQA();
0262 }