Tracking/dEdx/test_sample_size.C

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