tpc/dEdx/test_sample_size.C

0001 #include "GlobaldEdxFitter.h"
0002
0003 #include <TCanvas.h>
0004 #include <TH1F.h>
0005 #include <TH2F.h>
0006 #include <TGraphErrors.h>
0007 #include <TFile.h>
0008
0009 #include <algorithm>
0010 #include <numeric>
0011
0012 void test_sample_size(const std::string& infile="/sphenix/tg/tg01/hf/mjpeters/run53877_tracks/track_output_53877_*.root")
0013 {
0014   std::vector<float> samplesizes = {1000,2000,5000,10000,20000};//,50000,100000,200000,500000,1000000};//,2000000,5000000};
0015
0016   const int n_samples = 20;
0017   const float fluctuation_ymin = 5.;
0018   const float fluctuation_ymax = 26.;
0019   const int distribution_nbins = 30;
0020   const float distribution_xmin = 5.;
0021   const float distribution_xmax = 26.;
0022
0023   std::vector<std::vector<float>> fitvalues_all;
0024   std::vector<float> fitvalues_avg;
0025   std::vector<float> fitvalues_stdev;
0026
0027   std::vector<TCanvas*> fluctuations;
0028   std::vector<TCanvas*> distributions;
0029
0030   std::vector<TH1F*> dist_h;
0031
0032   std::vector<std::unique_ptr<GlobaldEdxFitter>> gfs;
0033
0034   for(int i=0; i<samplesizes.size(); i++)
0035   {
0036     fitvalues_all.emplace_back();
0037     gfs.push_back(std::make_unique<GlobaldEdxFitter>());
0038
0039     const std::string& fluctuation_canvasname = "fluctuations_"+std::to_string((int)floor(samplesizes[i]));
0040     const std::string& distribution_canvasname = "distributions_"+std::to_string((int)floor(samplesizes[i]));
0041     fluctuations.push_back(new TCanvas(fluctuation_canvasname.c_str(),fluctuation_canvasname.c_str(),600,600));
0042     distributions.push_back(new TCanvas(distribution_canvasname.c_str(),distribution_canvasname.c_str(),600,600));
0043
0044     for(int j=0;j<n_samples;j++)
0045     {
0046       gfs[i]->processResidualData(infile,floor(samplesizes[i]),j*samplesizes[i]);
0047       double min = gfs[i]->get_minimum();
0048       std::cout << "minimum: " << min << std::endl;
0049       fitvalues_all[i].push_back(min);
0050       if(j<n_samples-1)
0051       {
0052         gfs[i]->reset();
0053       }
0054 /*
0055       tf1s[i]->cd();
0056       TF1* tf1copy = gfs[i]->create_TF1(("ntrk_"+std::to_string(samplesizes[i])).c_str());
0057       tf1copy->SetLineColor(base_color);
0058       tf1copy->GetYaxis()->SetRangeUser(1.,tf1copy->GetMaximum());
0059       if(i==0) tf1copy->Draw();
0060       else tf1copy->Draw("SAME");
0061 */
0062     }
0063   }
0064
0065   std::vector<float> sample_index(n_samples);
0066   std::iota(sample_index.begin(),sample_index.end(),0.);
0067
0068   for(int i=0; i<samplesizes.size(); i++)
0069   {
0070     fluctuations[i]->cd();
0071     TGraph* g = new TGraph(n_samples,sample_index.data(),fitvalues_all[i].data());
0072     g->GetYaxis()->SetRangeUser(fluctuation_ymin,fluctuation_ymax);
0073     g->SetMarkerStyle(kFullCircle);
0074     g->SetMarkerSize(1.);
0075     g->Draw("APL");
0076
0077     distributions[i]->cd();
0078     std::string hname = "h_"+std::to_string(floor(samplesizes[i]));
0079     std::string htitle = "Distribution of fit results for sample size "+std::to_string(floor(samplesizes[i]));
0080 /*
0081     auto bounds = std::minmax_element(fitvalues_all[i].begin(),fitvalues_all[i].end());
0082     float lowerbound = floor(*bounds.first);
0083     float upperbound = ceil(*bounds.second);
0084 */
0085     TH1F* h = new TH1F(hname.c_str(),htitle.c_str(),distribution_nbins,distribution_xmin,distribution_xmax);
0086     for(int j=0; j<n_samples; j++)
0087     {
0088       h->Fill(fitvalues_all[i][j]);
0089     }
0090     h->Draw();
0091   }
0092
0093   for(int i=0; i<samplesizes.size(); i++)
0094   {
0095     float avg = 0.;
0096     float stdev = 0.;
0097     for(int j=0; j<n_samples; j++)
0098     {
0099       avg += fitvalues_all[i][j];
0100     }
0101     avg /= (float)n_samples;
0102
0103     for(int j=0; j<n_samples; j++)
0104     {
0105       stdev += pow(fitvalues_all[i][j]-avg,2.);
0106     }
0107     stdev /= (float)n_samples;
0108     stdev = std::sqrt(stdev);
0109
0110     fitvalues_avg.push_back(avg);
0111     fitvalues_stdev.push_back(stdev);
0112   }
0113
0114   std::vector<float> errx(n_samples,0.);
0115
0116   TCanvas* cg = new TCanvas("cg","sizes",600,600);
0117   TGraph* g = new TGraphErrors(samplesizes.size(),samplesizes.data(),fitvalues_avg.data(),errx.data(),fitvalues_stdev.data());
0118   g->SetMarkerStyle(kFullCircle);
0119   g->SetMarkerSize(1);
0120   g->Draw("APL");
0121   cg->SetLogx();
0122
0123   TCanvas* cbg = new TCanvas("vsbetagamma","vsbetagamma",600,600);
0124   TGraph* gbg = gfs.back()->graph_vsbetagamma(fitvalues_avg.back());
0125   gbg->SetMarkerStyle(kFullCircle);
0126   gbg->SetMarkerSize(0.2);
0127   gbg->Draw("AP");
0128   cbg->SetLogx();
0129
0130   double best_A = fitvalues_avg.back();
0131
0132   TF1* bethe = new TF1("bethebloch_vslnbg",bethe_bloch_new_1D_wrapper,0.,100.,2,1);
0133   bethe->SetParameter(0,best_A);
0134   bethe->SetNpx(1000);
0135   bethe->Draw("SAME");
0136
0137   TF1* bethe_directfit = new TF1("bethebloch_directfit",bethe_bloch_new_1D_wrapper,0.,10.,1,1);
0138   bethe_directfit->SetParameter(0,best_A);
0139   bethe_directfit->SetLineColor(kBlue);
0140   gbg->Fit(bethe_directfit);
0141   double newbest_A = bethe_directfit->GetParameter(0);
0142   std::cout << "new best: " << newbest_A << std::endl;
0143
0144   TCanvas* cbands = new TCanvas("bands","bands",600,600);
0145   TGraph* gp = gfs.back()->graph_vsp();
0146   gp->SetMarkerStyle(kFullCircle);
0147   gp->SetMarkerSize(0.1);
0148   gp->Draw("AP");
0149   cbands->SetLogx();
0150
0151   for(double mass : {dedx_constants::m_pi, dedx_constants::m_K, dedx_constants::m_p, dedx_constants::m_d})
0152   {
0153     TF1* band = new TF1(("band_"+std::to_string(mass)).c_str(),bethe_bloch_vs_p_wrapper_new_1D,0.,10.,2,1);
0154     band->SetParameter(0,best_A);
0155     band->SetParameter(1,mass);
0156     band->SetNpx(1000);
0157     band->Draw("SAME");
0158
0159     TF1* directband = new TF1(("directband_"+std::to_string(mass)).c_str(),bethe_bloch_vs_p_wrapper_new_1D,0.,10.,2,1);
0160     directband->SetLineColor(kBlue);
0161     directband->SetParameters(best_A,mass);
0162     directband->SetNpx(1000);
0163     directband->Draw("SAME");
0164   }
0165
0166   TCanvas* cb = new TCanvas("fullbands","fullbands",600,600);
0167   TFile* f_h = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/run53877_tracks/dedx/merged_dedx.root");
0168   TH2F* dedx_h = (TH2F*)f_h->Get("dedx_log_30");
0169   dedx_h->Draw("COLZ");
0170   cb->SetLogz();
0171
0172   for(double mass : {dedx_constants::m_pi, dedx_constants::m_K, dedx_constants::m_p, dedx_constants::m_d})
0173   {
0174     TF1* band = new TF1(("band_"+std::to_string(mass)).c_str(),bethe_bloch_vs_logp_wrapper_new_1D,-1.,5.,2,1);
0175     band->SetParameter(0,best_A);
0176     band->SetParameter(1,mass);
0177     band->Draw("SAME");
0178
0179     TF1* directband = new TF1(("directband_"+std::to_string(mass)).c_str(),bethe_bloch_vs_logp_wrapper_new_1D,-1.,5.,2,1);
0180     directband->SetLineColor(kBlue);
0181     directband->SetParameters(newbest_A,mass);
0182     directband->SetNpx(1000);
0183     directband->Draw("SAME");
0184   }
0185
0186   TFile* fout = new TFile("dedxfitvals.root","RECREATE");
0187   for(auto& c : fluctuations)
0188   {
0189     c->Write();
0190   }
0191   for(auto& c : distributions)
0192   {
0193     c->Write();
0194   }
0195   cg->Write();
0196   cbg->Write();
0197   cbands->Write();
0198   cb->Write();
0199   fout->Close();
0200 }