Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:12:27

0001 #ifndef EtaDistReader_h
0002 #define EtaDistReader_h
0003 
0004 #include <filesystem>
0005 #include <cstdlib>
0006 
0007 #include "INTTEta.h"
0008 
0009 
0010 // note : the code works with the histograms, so the centrality bins are already defined in the histograms, ranging from 0 to something big number
0011 // todo : the signal_region and centrality_region is hard_written in the INTTEta.h, which is written in the "ana_map_v1.h", so modify there
0012 
0013 
0014 class EtaDistReader : public INTTEta
0015 {
0016     public:
0017         EtaDistReader(string run_type, string out_folder_directory, vector<int> fulleta_MC_included_z_bin, vector<pair<int,vector<int>>> included_eta_z_map, string input_file_directory, bool centrality_Z_map_bool = 0) : 
0018         INTTEta(run_type, out_folder_directory), fulleta_MC_included_z_bin(fulleta_MC_included_z_bin), included_eta_z_map(included_eta_z_map), input_file_directory(input_file_directory), centrality_Z_map_bool(centrality_Z_map_bool)
0019         {
0020             gErrorIgnoreLevel = kError;
0021 
0022             cout<<"In EtaDistReader, centrality region size : "<<centrality_region.size()<<endl;    
0023             cout<<"In EtaDistReader, eta region size : "<<eta_region.size() - 1<<endl;
0024             cout<<"In EtaDistReader, plot_text: "<<plot_text<<endl;
0025             cout<<"In EtaDistReader, signal_region: "<<signal_region<<endl;
0026 
0027             if(std::filesystem::exists(Form("%s",out_folder_directory.c_str())) == false){
0028                 system(Form("mkdir %s", out_folder_directory.c_str()));
0029                 cout<<"----------- check folder exists -----------"<<endl;
0030                 system(Form("ls %s", out_folder_directory.c_str()));
0031             }
0032             else 
0033             {
0034                 cout<<"----------- folder exists already -----------"<<endl;
0035                 system(Form("ls %s", out_folder_directory.c_str()));
0036             }
0037 
0038 
0039             file_in = TFile::Open(input_file_directory.c_str());
0040 
0041             N_centrality_bin = centrality_region.size();
0042             eta_correction = included_eta_z_map[0].first - 1; // todo : the eta bin starts from 15, so the correction is 14
0043             
0044             SignalNTrack_eta_z_Multi_2D.clear();
0045             SignalNTrack_eta_z_Single_2D.clear();
0046             TrueNtrack_eta_z_MC_2D.clear();
0047             TrueNtrack_eta_z_fulleta_MC_2D.clear();
0048             dNdeta_1D_MC.clear();
0049             dNdeta_1D_fulleta_MC.clear();
0050             dNdeta_1D_reco_single_original.clear();
0051             dNdeta_1D_reco_multi_original.clear();
0052             dNdeta_1D_reco_single_postalpha = vector<TH1F*>(N_centrality_bin, nullptr);
0053             dNdeta_1D_reco_multi_postalpha = vector<TH1F*>(N_centrality_bin, nullptr);
0054             DeltaPhi_Multi_1D.clear();
0055             DeltaPhi_Multi_Stack_hist_out.clear();
0056             DeltaPhi_Multi_Stack.clear();
0057             SignalNTrack_Single.clear(); 
0058             SignalNTrack_Multi.clear();
0059             N_event_counting            = vector<vector<double>>(N_centrality_bin, vector<double>(included_eta_z_map.size(),0));
0060             N_event_counting_MC         = vector<vector<int>>(N_centrality_bin, vector<int>(included_eta_z_map.size(),0));
0061             N_event_counting_fulleta_MC = vector<vector<int>>(N_centrality_bin, vector<int>(eta_region.size() - 1, 0)); // note : N centrality_bin x N eta_bin
0062 
0063             DeltaPhi_Multi_Stack_hist_out_radian.clear();
0064             DeltaPhi_Multi_Stack_radian.clear();
0065             DeltaPhi_Multi_1D_radian.clear();
0066             
0067             
0068             solid_line = new TLine();
0069             solid_line -> SetLineWidth(1);
0070             solid_line -> SetLineColor(2);
0071             solid_line -> SetLineStyle(1);
0072 
0073             // file_out = new TFile(Form("%s/alpha_correction_map.root",out_folder_directory.c_str()),"RECREATE");
0074 
0075             legend = new TLegend(0.45,0.8,0.70,0.9);
0076             // legend -> SetMargin(0);
0077             legend->SetTextSize(0.03);
0078 
0079             post_alpha_tag = 0;
0080 
0081             ReadFileHist();
0082             InitHist();
0083             MainPreparation();
0084             HistDivision();
0085             PrintPlots_exclusive();
0086             // FinaldNdEta();
0087 
0088             return;
0089         };
0090 
0091         vector<TH1F *> GetDeltaPhi_Multi_stack_1D();
0092         vector<TH1F *> GetDeltaPhi_Multi_stack_1D_radian();
0093         vector<TH1F *> GetdNdeta_1D_MC();
0094         vector<TH1F *> GetdNdeta_1D_fulleta_MC();
0095         vector<TH1F *> GetdNdeta_1D_reco_single_original();
0096         vector<TH1F *> GetdNdeta_1D_reco_multi_original();
0097         vector<TH1F *> GetdNdeta_1D_reco_single_postalpha();
0098         vector<TH1F *> GetdNdeta_1D_reco_multi_postalpha();
0099         vector<string> Get_centrality_region() {return centrality_region;};
0100         string         Get_eta_region_txt(int bin_x_id);
0101         vector<TH2F *> Get_alpha_corr_map();
0102         void FindCoveredRegion();
0103         void DeriveAlphaCorrection(int correction_step);
0104         void Set_alpha_corr_map(TH2F * hist_in_loose, TH2F * hist_in_tight);
0105         void ApplyAlphaCorrection();
0106         void PrintPlots();
0107         void PrintPlots_exclusive();
0108         // void FinaldNdEta();
0109         void EndRun();
0110 
0111         ~EtaDistReader()
0112         {
0113             file_in -> Close();
0114             cout<<"EtaDistReader done, goodbye"<<endl;
0115         };
0116 
0117     protected:
0118 
0119         bool centrality_Z_map_bool;
0120         vector<pair<int,vector<int>>> included_eta_z_map;
0121         vector<int> fulleta_MC_included_z_bin;
0122         string input_file_directory;
0123         int eta_correction; // todo : the eta bin starts from 15, so the correction is 14
0124         int N_centrality_bin;
0125 
0126         TFile * file_in;
0127         string color_code[5] = {"#fedc97", "#b5b682", "#7c9885", "#28666e", "#033f63"};
0128         // vector<string> centrality_region;
0129 
0130         vector<TH2F *> SignalNTrack_eta_z_Multi_2D;
0131         vector<TH2F *> SignalNTrack_eta_z_Single_2D;
0132         map<string,TH1F *> DeltaPhi_Multi_1D;
0133         vector<TH2F *> TrueNtrack_eta_z_MC_2D;
0134         vector<TH2F *> TrueNtrack_eta_z_fulleta_MC_2D;
0135         TH2F * centrality_Z_map;
0136         TH2F * centrality_Z_map_MC;
0137         TH2F * eta_z_ref;
0138         TH2F * eta_z_ref_used_check;
0139         TH2F * eta_z_ref_used_check_binID;
0140 
0141         TH2F * eta_z_ref_cm;
0142         TH2F * eta_z_ref_used_check_cm;
0143         TH2F * eta_z_ref_used_check_binID_cm;
0144         TH2F * exclusive_cluster_inner_eta_phi_2D;
0145         TH2F * exclusive_cluster_outer_eta_phi_2D;
0146 
0147         TH2F * exclusive_cluster_inner_eta_adc_2D;
0148         TH2F * exclusive_cluster_outer_eta_adc_2D;
0149         TH2F * exclusive_cluster_inner_eta_adc_2D_bkgrm;
0150         TH2F * exclusive_cluster_outer_eta_adc_2D_bkgrm;
0151         TProfile * exclusive_cluster_inner_eta_adc_2D_bkgrm_profile;
0152         TProfile * exclusive_cluster_outer_eta_adc_2D_bkgrm_profile;
0153         TGraphErrors * exclusive_cluster_inner_eta_adc_2D_bkgrm_profile_graph;
0154         TGraphErrors * exclusive_cluster_outer_eta_adc_2D_bkgrm_profile_graph;
0155 
0156 
0157         TH1F * exclusive_cluster_inner_adc;
0158         TH1F * exclusive_cluster_outer_adc;
0159 
0160         TH2F * eta_Mbin_correction_tight;
0161         TH2F * eta_Mbin_correction_loose;
0162         TH2F * eta_Mbin_correction_loose_noUPC; // note : no ultra peripheral collision bin
0163         TH2F * input_eta_Mbin_correction_tight;
0164         TH2F * input_eta_Mbin_correction_loose;
0165 
0166 
0167         vector<TH1F *> dNdeta_1D_MC;
0168         vector<TH1F *> dNdeta_1D_fulleta_MC;
0169         vector<TH1F *> dNdeta_1D_reco_single_original;
0170         vector<TH1F *> dNdeta_1D_reco_multi_original;
0171         vector<TH1F *> dNdeta_1D_reco_single_postalpha;
0172         vector<TH1F *> dNdeta_1D_reco_multi_postalpha;
0173         
0174 
0175         map<string, THStack *> DeltaPhi_Multi_Stack;
0176         vector<TH1F *> DeltaPhi_Multi_Stack_hist_out;
0177         map<string, THStack *> DeltaPhi_Multi_Stack_radian;
0178         vector<TH1F *> DeltaPhi_Multi_Stack_hist_out_radian;
0179         map<string,TH1F *> DeltaPhi_Multi_1D_radian;
0180 
0181         map<string, double> SignalNTrack_Single;
0182         map<string, double> SignalNTrack_Multi;
0183         vector<vector<double>> N_event_counting;
0184         vector<vector<int>> N_event_counting_MC;
0185         vector<vector<int>> N_event_counting_fulleta_MC;
0186 
0187         TLegend * legend;
0188         TLine * solid_line;
0189         TH1F * temp_hist;
0190 
0191         int post_alpha_tag;
0192 
0193         double radian_conversion = 180./M_PI;
0194 
0195         // TFile * file_out; 
0196 
0197         void ReadFileHist();
0198         void InitHist();
0199         void MainPreparation();
0200         void DrawCoordLine(TH2F * hist_in);
0201         void DrawEtaCoverBox(TH2F * hist_in, string unit_string);
0202         void HistDivision();
0203         void Hist1DSetting(TH1F * hist_in, string color_code);
0204         
0205 
0206 };
0207 
0208 string EtaDistReader::Get_eta_region_txt(int bin_x_id)
0209 {
0210     return Form("%.2f ~ %.2f",eta_z_ref->GetXaxis()->GetBinLowEdge(bin_x_id), eta_z_ref->GetXaxis()->GetBinLowEdge(bin_x_id) + eta_z_ref->GetXaxis()->GetBinWidth(bin_x_id));
0211 }
0212 
0213 void EtaDistReader::DrawCoordLine(TH2F * hist_in)
0214 {
0215     // note : draw the vertical line, to segment the eta region
0216     for (int i = 0; i < hist_in -> GetNbinsX(); i++) { 
0217         coord_line -> DrawLine(
0218             hist_in->GetXaxis()->GetBinLowEdge(i+1), 
0219             hist_in->GetYaxis()->GetBinLowEdge(1), 
0220             hist_in->GetXaxis()->GetBinLowEdge(i+1), 
0221             hist_in->GetYaxis()->GetBinLowEdge(hist_in -> GetNbinsY()) + hist_in->GetYaxis()->GetBinWidth(hist_in -> GetNbinsY())
0222         ); 
0223     }
0224     coord_line -> DrawLine(
0225         hist_in->GetXaxis()->GetBinLowEdge(hist_in->GetNbinsX()) + hist_in->GetXaxis()->GetBinWidth(hist_in->GetNbinsX()), 
0226         hist_in->GetYaxis()->GetBinLowEdge(1), 
0227         hist_in->GetXaxis()->GetBinLowEdge(hist_in->GetNbinsX()) + hist_in->GetXaxis()->GetBinWidth(hist_in->GetNbinsX()), 
0228         hist_in->GetYaxis()->GetBinLowEdge(hist_in -> GetNbinsY()) + hist_in->GetYaxis()->GetBinWidth(hist_in -> GetNbinsY())
0229     );
0230 
0231     // note : draw the horizontal line, to segment the z region
0232     for (int i = 0; i < hist_in -> GetNbinsY(); i++) { 
0233         coord_line -> DrawLine(
0234             hist_in->GetXaxis()->GetBinLowEdge(1), 
0235             hist_in->GetYaxis()->GetBinLowEdge(i+1), 
0236             hist_in->GetXaxis()->GetBinLowEdge(hist_in->GetNbinsX()) + hist_in->GetXaxis()->GetBinWidth(hist_in->GetNbinsX()), 
0237             hist_in->GetYaxis()->GetBinLowEdge(i+1)
0238         ); 
0239     }
0240     coord_line -> DrawLine(
0241         hist_in->GetXaxis()->GetBinLowEdge(1), 
0242         hist_in->GetYaxis()->GetBinLowEdge(hist_in -> GetNbinsY()) + hist_in->GetYaxis()->GetBinWidth(hist_in -> GetNbinsY()), 
0243         hist_in->GetXaxis()->GetBinLowEdge(hist_in->GetNbinsX()) + hist_in->GetXaxis()->GetBinWidth(hist_in->GetNbinsX()), 
0244         hist_in->GetYaxis()->GetBinLowEdge(hist_in -> GetNbinsY()) + hist_in->GetYaxis()->GetBinWidth(hist_in -> GetNbinsY())
0245     ); 
0246 }
0247 
0248 void EtaDistReader::DrawEtaCoverBox(TH2F * hist_in, string unit_string)
0249 {
0250     double unit_correction;
0251     if (unit_string == "mm") { unit_correction = 1.; }
0252     else if (unit_string == "cm") { unit_correction = 0.1; }
0253     else {cout<<"In EtaDistReader, wrong unit_string input"<<endl;}
0254 
0255     double INTT_layer_R[4] = {
0256         71.88 * unit_correction, 
0257         77.32 * unit_correction, 
0258         96.8 * unit_correction, 
0259         102.62 * unit_correction
0260     }; // note : the radii of the INTT layers
0261 
0262     double INTT_left_edge = -230.; INTT_left_edge = INTT_left_edge * unit_correction;
0263     double INTT_right_edge = 230.; INTT_right_edge = INTT_right_edge * unit_correction;
0264 
0265     // note : low  Z edge for left  eta
0266     // note : high Z edge for right eta
0267     for (int i = 0; i < hist_in -> GetNbinsY(); i++)
0268     {
0269         double bin_low_zvtx  = hist_in->GetYaxis()->GetBinLowEdge(i+1);
0270         double bin_high_zvtx = hist_in->GetYaxis()->GetBinLowEdge(i+1) + hist_in->GetYaxis()->GetBinWidth(i+1);
0271 
0272         double INTT_eta_acceptance_l = -0.5 * TMath::Log((sqrt(pow(INTT_left_edge - bin_low_zvtx,2)+pow(INTT_layer_R[0],2))-(INTT_left_edge - bin_low_zvtx)) / (sqrt(pow(INTT_left_edge - bin_low_zvtx,2)+pow(INTT_layer_R[0],2))+(INTT_left_edge - bin_low_zvtx))); // note : left
0273         double INTT_eta_acceptance_r = -0.5 * TMath::Log((sqrt(pow(INTT_right_edge - bin_high_zvtx,2)+pow(INTT_layer_R[3],2))-(INTT_right_edge - bin_high_zvtx)) / (sqrt(pow(INTT_right_edge - bin_high_zvtx,2)+pow(INTT_layer_R[3],2))+(INTT_right_edge - bin_high_zvtx))); // note : right
0274         
0275         solid_line -> DrawLine(INTT_eta_acceptance_l, bin_low_zvtx, INTT_eta_acceptance_l, bin_high_zvtx); // note : vertical
0276         solid_line -> DrawLine(INTT_eta_acceptance_r, bin_low_zvtx, INTT_eta_acceptance_r, bin_high_zvtx); // note : vertical
0277         solid_line -> DrawLine(INTT_eta_acceptance_l, bin_low_zvtx, INTT_eta_acceptance_r, bin_low_zvtx); // note : horizontal
0278         solid_line -> DrawLine(INTT_eta_acceptance_l, bin_high_zvtx, INTT_eta_acceptance_r, bin_high_zvtx); // note : horizontal
0279         
0280         cout<<"bin Yaxis ID: "<< i+1 <<" low_edge "<<bin_low_zvtx<<" high edge: "<<bin_high_zvtx<<" eta region: "<<INTT_eta_acceptance_l<<"~"<<INTT_eta_acceptance_r<<endl;   
0281     }
0282 }
0283 
0284 void EtaDistReader::ReadFileHist()
0285 {
0286 
0287     // note : the histograms of the number of tracklets in the signal region for both methods
0288     // note : the 1D vector is for different centrality bins
0289     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0290     {
0291         // note : the entries with the deltaphi within 1 degree, which means the background entries are also included. 
0292         // note : the loose method
0293         SignalNTrack_eta_z_Multi_2D.push_back((TH2F *) file_in -> Get(Form("Reco_SignalNTracklet_Multi_MBin%i",Mbin)));
0294         // note : the tight method
0295         SignalNTrack_eta_z_Single_2D.push_back((TH2F *) file_in -> Get(Form("Reco_SignalNTracklet_Single_MBin%i",Mbin)));
0296     }
0297 
0298     // note : the 1D histograms of the DeltaPhi distributions of each centrality-eta-z bin (the multi method, for the background subtraction)
0299     // note : the format of the key is "MBin_EtaBin_ZBin"
0300     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++) // note : centrality bin
0301     {
0302         for (auto eta_z : included_eta_z_map) // note : eta bin
0303         {
0304             for(int zbin = 0; zbin < eta_z.second.size(); zbin++) // note :  Z bin
0305             {
0306                 DeltaPhi_Multi_1D[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] = (TH1F *) file_in -> Get(Form("Reco_DeltaPhi1D_Multi_MBin%i_Eta%i_Z%i",Mbin,eta_z.first,eta_z.second[zbin]));
0307                 DeltaPhi_Multi_1D[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] -> SetFillColor(TColor::GetColor(color_code[zbin%5].c_str()));
0308 
0309                 // note : change the unit, for the presentation
0310                 DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] = (TH1F *) file_in -> Get(Form("Reco_DeltaPhi1D_Multi_MBin%i_Eta%i_Z%i",Mbin,eta_z.first,eta_z.second[zbin]));
0311                 DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] -> SetFillColor(TColor::GetColor(color_code[zbin%5].c_str()));
0312                 DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] -> GetXaxis() -> SetLimits(
0313                     DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] -> GetXaxis() -> GetXmin() / radian_conversion,
0314                     DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,eta_z.second[zbin])] -> GetXaxis() -> GetXmax() / radian_conversion
0315                 );
0316             }
0317         }
0318     }
0319 
0320     // note : the histograms of the number of True track from the MC
0321     // note : the 1D vector is for different centrality bins
0322     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0323     {
0324         TrueNtrack_eta_z_MC_2D.push_back((TH2F *) file_in -> Get(Form("NTrueTrack_MBin%i",Mbin)));
0325     }
0326 
0327     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0328     {
0329         TrueNtrack_eta_z_fulleta_MC_2D.push_back((TH2F *) file_in -> Get(Form("NTrueTrack_FullEta_MBin%i",Mbin)));
0330     }
0331 
0332     // note : the centrality Z map, to keep the N event counting
0333     string centrality_Z_map_name = (centrality_Z_map_bool == 0) ? "MBin_Z_evt_map" : "Z_MBin_evt_map"; // todo : the correct name should be "MBin_Z_evt_map", correct next time
0334     centrality_Z_map    = (TH2F *) file_in -> Get(centrality_Z_map_name.c_str());
0335     centrality_Z_map_MC = (TH2F *) file_in -> Get((centrality_Z_map_name+"_MC").c_str()); 
0336 
0337     // note : the reference map of Eta-Z 2D histogram
0338     eta_z_ref = (TH2F *) file_in -> Get("Eta_Z_reference");
0339 
0340     exclusive_cluster_inner_eta_phi_2D = (TH2F *) file_in -> Get("exclusive_cluster_inner_eta_phi_2D");
0341     exclusive_cluster_outer_eta_phi_2D = (TH2F *) file_in -> Get("exclusive_cluster_outer_eta_phi_2D");
0342     exclusive_cluster_inner_eta_adc_2D = (TH2F *) file_in -> Get("exclusive_cluster_inner_eta_adc_2D");
0343     exclusive_cluster_outer_eta_adc_2D = (TH2F *) file_in -> Get("exclusive_cluster_outer_eta_adc_2D");
0344     
0345     exclusive_cluster_inner_eta_adc_2D_bkgrm = (TH2F *) file_in -> Get("exclusive_cluster_inner_eta_adc_2D");
0346     // TH2F_threshold_advanced_2(exclusive_cluster_inner_eta_adc_2D_bkgrm, 0.5);
0347     vector<double> hist_row_content_inner = SumTH2FColumnContent_row(exclusive_cluster_inner_eta_adc_2D_bkgrm);
0348     
0349     exclusive_cluster_inner_eta_adc_2D_bkgrm_profile = (TProfile *) exclusive_cluster_inner_eta_adc_2D_bkgrm -> ProfileY("exclusive_cluster_inner_eta_adc_2D_bkgrm_profile");
0350     exclusive_cluster_inner_eta_adc_2D_bkgrm_profile_graph = new TGraphErrors(); 
0351     int point_index_inner = 0;
0352     for (int i = 0; i < hist_row_content_inner.size(); i++){
0353         if (hist_row_content_inner[i] < 5){continue;} // note : in order to remove some remaining background
0354         exclusive_cluster_inner_eta_adc_2D_bkgrm_profile_graph -> SetPoint(
0355             point_index_inner, 
0356             exclusive_cluster_inner_eta_adc_2D_bkgrm_profile->GetBinContent(i+1), 
0357             exclusive_cluster_inner_eta_adc_2D_bkgrm_profile->GetBinCenter(i+1)
0358         );
0359         
0360         exclusive_cluster_inner_eta_adc_2D_bkgrm_profile_graph -> SetPointError(
0361             point_index_inner, 
0362             exclusive_cluster_inner_eta_adc_2D_bkgrm_profile->GetBinError(i+1), 
0363             exclusive_cluster_inner_eta_adc_2D_bkgrm_profile->GetBinWidth(i+1)/2.
0364         );
0365         
0366         point_index_inner += 1;   
0367     }
0368 
0369     exclusive_cluster_outer_eta_adc_2D_bkgrm = (TH2F *) file_in -> Get("exclusive_cluster_outer_eta_adc_2D");
0370     // TH2F_threshold_advanced_2(exclusive_cluster_outer_eta_adc_2D_bkgrm, 0.5);
0371     vector<double> hist_row_content_outer = SumTH2FColumnContent_row(exclusive_cluster_outer_eta_adc_2D_bkgrm);
0372     
0373     exclusive_cluster_outer_eta_adc_2D_bkgrm_profile = (TProfile *) exclusive_cluster_outer_eta_adc_2D_bkgrm -> ProfileY("exclusive_cluster_outer_eta_adc_2D_bkgrm_profile");
0374     exclusive_cluster_outer_eta_adc_2D_bkgrm_profile_graph = new TGraphErrors(); 
0375     int point_index_outer = 0;
0376     for (int i = 0; i < hist_row_content_outer.size(); i++){
0377         if (hist_row_content_outer[i] < 5){continue;} // note : in order to remove some remaining background
0378         exclusive_cluster_outer_eta_adc_2D_bkgrm_profile_graph -> SetPoint(
0379             point_index_outer, 
0380             exclusive_cluster_outer_eta_adc_2D_bkgrm_profile->GetBinContent(i+1), 
0381             exclusive_cluster_outer_eta_adc_2D_bkgrm_profile->GetBinCenter(i+1)
0382         );
0383         exclusive_cluster_outer_eta_adc_2D_bkgrm_profile_graph -> SetPointError(
0384             point_index_outer, 
0385             exclusive_cluster_outer_eta_adc_2D_bkgrm_profile->GetBinError(i+1), 
0386             exclusive_cluster_outer_eta_adc_2D_bkgrm_profile->GetBinWidth(i+1)/2.
0387         );
0388         
0389         point_index_outer += 1;   
0390     }
0391 
0392     
0393 
0394     exclusive_cluster_inner_adc = (TH1F *) file_in -> Get("exclusive_cluster_inner_adc");
0395     exclusive_cluster_outer_adc = (TH1F *) file_in -> Get("exclusive_cluster_outer_adc");
0396 
0397     eta_z_ref_used_check = (TH2F*) eta_z_ref -> Clone("eta_z_ref_used_check");
0398     eta_z_ref_used_check_binID = (TH2F*) eta_z_ref -> Clone("eta_z_ref_used_check");
0399     
0400     // todo, the unit_correction is hard written here, the 0.1, I meant. 
0401     eta_z_ref_cm = new TH2F (
0402         "",
0403         "eta_z_ref_cm;reco #eta;reco zvtx [cm]",
0404         eta_z_ref -> GetNbinsX(),
0405         eta_z_ref -> GetXaxis() -> GetXmin(),
0406         eta_z_ref -> GetXaxis() -> GetXmax(),
0407         eta_z_ref -> GetNbinsY(),
0408         eta_z_ref -> GetYaxis() -> GetXmin() * 0.1,
0409         eta_z_ref -> GetYaxis() -> GetXmax() * 0.1
0410     );
0411     eta_z_ref_used_check_cm = (TH2F *) eta_z_ref_cm -> Clone("eta_z_ref_used_check_cm");
0412     eta_z_ref_used_check_binID_cm = (TH2F *) eta_z_ref_cm -> Clone("eta_z_ref_used_check_binID_cm");
0413 
0414     cout<<"ReadFileHist done"<<endl;
0415 
0416     return;
0417 
0418 }
0419 
0420 void EtaDistReader::FindCoveredRegion()
0421 {
0422     for (int eta_bin = 0; eta_bin < eta_z_ref_used_check_binID->GetNbinsX(); eta_bin++)
0423     {
0424         for (int z_bin = 0; z_bin < eta_z_ref_used_check_binID->GetNbinsY(); z_bin++)
0425         {
0426             eta_z_ref_used_check_binID    -> SetBinContent(eta_bin+1,z_bin+1, 100 * (eta_bin+1) + (z_bin+1));
0427             eta_z_ref_used_check_binID_cm -> SetBinContent(eta_bin+1,z_bin+1, 100 * (eta_bin+1) + (z_bin+1));
0428         }
0429     }
0430 
0431     c1 -> cd();
0432     eta_z_ref -> GetXaxis() -> SetTitle("reco #eta");
0433     eta_z_ref -> GetYaxis() -> SetTitle("reco zvtx [mm]");
0434     eta_z_ref -> Draw("colz0");
0435     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0436     DrawCoordLine(eta_z_ref);    
0437     DrawEtaCoverBox(eta_z_ref, "mm");
0438     c1 -> Print(Form("%s/eta_z_ref.pdf",out_folder_directory.c_str()));
0439     c1 -> Clear();
0440 
0441     // note : for the used bins check
0442     c1 -> cd();
0443     double eta_z_ref_used_check_max_entry = eta_z_ref_used_check->GetMaximum();
0444     for (auto eta_z : included_eta_z_map){
0445         for (auto zbin : eta_z.second){
0446             eta_z_ref_used_check -> SetBinContent(eta_z.first, zbin, eta_z_ref_used_check_max_entry * 1000);        
0447         }
0448     }
0449     eta_z_ref_used_check -> GetXaxis() -> SetTitle("reco #eta");
0450     eta_z_ref_used_check -> GetYaxis() -> SetTitle("reco zvtx [mm]");
0451     eta_z_ref_used_check -> Draw("colz0");
0452     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0453     DrawCoordLine(eta_z_ref_used_check);    
0454     DrawEtaCoverBox(eta_z_ref_used_check, "mm");
0455 
0456     gStyle->SetPaintTextFormat("1.0f");
0457     eta_z_ref_used_check_binID -> SetMarkerSize(0.7);
0458     eta_z_ref_used_check_binID -> Draw("HIST TEXT45 SAME");
0459 
0460     c1 -> Print(Form("%s/eta_z_ref_used_check.pdf",out_folder_directory.c_str()));
0461     c1 -> Clear();
0462 
0463     // note : unit in cm, for the presentation
0464     TH2F_FakeClone(eta_z_ref, eta_z_ref_cm);
0465     TH2F_FakeClone(eta_z_ref_used_check, eta_z_ref_used_check_cm);
0466 
0467     c1 -> cd();
0468     eta_z_ref_cm -> GetXaxis() -> SetTitle("reco #eta");
0469     eta_z_ref_cm -> GetYaxis() -> SetTitle("reco zvtx [cm]");
0470     eta_z_ref_cm -> Draw("colz0");
0471     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0472     DrawCoordLine(eta_z_ref_cm);    
0473     DrawEtaCoverBox(eta_z_ref_cm, "cm");
0474     c1 -> Print(Form("%s/eta_z_ref_cm.pdf",out_folder_directory.c_str()));
0475     c1 -> Clear();
0476 
0477     // note : for the used bins check
0478     c1 -> cd();
0479     eta_z_ref_used_check_cm -> GetXaxis() -> SetTitle("reco #eta");
0480     eta_z_ref_used_check_cm -> GetYaxis() -> SetTitle("reco zvtx [cm]");
0481     eta_z_ref_used_check_cm -> Draw("colz0");
0482     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0483     DrawCoordLine(eta_z_ref_used_check_cm);    
0484     DrawEtaCoverBox(eta_z_ref_used_check_cm, "cm");
0485 
0486     gStyle->SetPaintTextFormat("1.0f");
0487     eta_z_ref_used_check_binID_cm -> SetMarkerSize(0.7);
0488     eta_z_ref_used_check_binID_cm -> Draw("HIST TEXT45 SAME");
0489 
0490     c1 -> Print(Form("%s/eta_z_ref_used_check_cm.pdf",out_folder_directory.c_str()));
0491     c1 -> Clear();
0492 
0493 }
0494 
0495 void EtaDistReader::Hist1DSetting(TH1F * hist_in, string color_code)
0496 {
0497     hist_in -> SetMarkerStyle(20);
0498     hist_in -> SetMarkerSize(0.8);
0499     hist_in -> SetMarkerColor(TColor::GetColor(color_code.c_str()));
0500     hist_in -> SetLineColor(TColor::GetColor(color_code.c_str()));
0501     hist_in -> SetLineWidth(2);
0502     hist_in -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0503     hist_in -> GetXaxis() -> SetTitle("#eta");
0504     // hist_in -> GetYaxis() -> SetRangeUser(0,50);
0505     hist_in -> SetTitleSize(0.06, "X");
0506     hist_in -> SetTitleSize(0.06, "Y");
0507     hist_in -> GetXaxis() -> SetTitleOffset (0.71);
0508     hist_in -> GetYaxis() -> SetTitleOffset (1.1);
0509     hist_in -> GetXaxis() -> CenterTitle(true);
0510     hist_in -> GetYaxis() -> CenterTitle(true);
0511 }
0512 
0513 void EtaDistReader::InitHist()
0514 {
0515     for (int i = 0; i < N_centrality_bin; i++)
0516     {
0517         dNdeta_1D_MC.push_back(new TH1F("","",included_eta_z_map.size(), eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.));
0518         Hist1DSetting(dNdeta_1D_MC[i], "#1A3947");
0519 
0520         dNdeta_1D_fulleta_MC.push_back(
0521             new TH1F(
0522                 "",
0523                 "",
0524                 eta_region.size() - 1,
0525                 &eta_region[0]
0526             )
0527         );
0528         Hist1DSetting(dNdeta_1D_fulleta_MC[i], "#D8364D");
0529 
0530         dNdeta_1D_reco_single_original.push_back(new TH1F("","",included_eta_z_map.size(), eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.));
0531         Hist1DSetting(dNdeta_1D_reco_single_original[i], "#1A3947");
0532 
0533         dNdeta_1D_reco_multi_original.push_back(new TH1F("","",included_eta_z_map.size(), eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.));
0534         Hist1DSetting(dNdeta_1D_reco_multi_original[i], "#c48045");
0535     }
0536 
0537     eta_Mbin_correction_tight = new TH2F(
0538         "",
0539         "eta_Mbin_correction_tight;#eta;Mbin;#alpha factor",
0540         included_eta_z_map.size(), 
0541         eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), 
0542         eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.,
0543         N_centrality_bin,
0544         0, 
0545         N_centrality_bin
0546     );
0547 
0548     eta_Mbin_correction_loose = new TH2F(
0549         "",
0550         "eta_Mbin_correction_loose;#eta;Mbin;#alpha factor",
0551         included_eta_z_map.size(), 
0552         eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), 
0553         eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.,
0554         N_centrality_bin,
0555         0, 
0556         N_centrality_bin
0557     );
0558 
0559     eta_Mbin_correction_loose_noUPC = new TH2F(
0560         "",
0561         "eta_Mbin_correction_loose_noUPC;#eta;Mbin;#alpha factor",
0562         included_eta_z_map.size(), 
0563         eta_z_ref -> GetXaxis() ->GetBinLowEdge(included_eta_z_map[0].first), 
0564         eta_z_ref -> GetXaxis() ->GetBinCenter(included_eta_z_map[included_eta_z_map.size() - 1].first) + eta_z_ref -> GetXaxis() ->GetBinWidth(included_eta_z_map[included_eta_z_map.size() - 1].first) / 2.,
0565         N_centrality_bin,
0566         0, 
0567         N_centrality_bin
0568     );
0569 
0570     input_eta_Mbin_correction_tight = nullptr;
0571     input_eta_Mbin_correction_loose = nullptr;
0572 
0573     cout<<"InitHist done"<<endl;
0574 }
0575 
0576 // note : to count the number of reco tracklets in the signal region with the selected eta bin and z bin, for each centrality bin
0577 // note : for both methods
0578 // note : to increase the statistic, try to stack the DeltaPhi distributions with same eta bin but different z bin
0579 // note : the key is "MBin_EtaBin"
0580 // note : for each eta bin and each centrality bin, count the number of event. (Each eta bin has different selected z region, so the number of event is different for each eta bin)
0581 void EtaDistReader::MainPreparation()
0582 {
0583     c1 -> cd();
0584     c1 -> Print(Form("%s/Reco_DeltaPhi_Multi_Stack.pdf(", out_folder_directory.c_str())); 
0585     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++) // note : centrality bin
0586     {
0587         // note : for the MC
0588         // note: this is fro the full dN/deta distribution in the MC, the full coverage
0589         for (int eta_i = 0; eta_i < TrueNtrack_eta_z_fulleta_MC_2D[Mbin]->GetNbinsX(); eta_i++)
0590         {
0591             for (int z_bin : fulleta_MC_included_z_bin) // note : the bin ID is used.
0592             {
0593                 dNdeta_1D_fulleta_MC[Mbin]->SetBinContent(eta_i+1, dNdeta_1D_fulleta_MC[Mbin] -> GetBinContent(eta_i+1) + TrueNtrack_eta_z_fulleta_MC_2D[Mbin]->GetBinContent(eta_i+1, z_bin));   
0594                 N_event_counting_fulleta_MC[Mbin][eta_i] += centrality_Z_map_MC -> GetBinContent(Mbin + 1, z_bin);
0595             }
0596         }
0597 
0598         // note : for the MC
0599         //todo: test, which means that, for the true dN/deta distribution, we don't use the selected region only. Instead, we use the full range (kinda directly from the physics)
0600         //todo: in principle, the two methods should work, but turned out that at the edge bins, the entries are a bit fewer than the expectation. 
0601         for (int eta_i = 0; eta_i < dNdeta_1D_MC[Mbin]->GetNbinsX(); eta_i++) {
0602             dNdeta_1D_MC[Mbin] -> SetBinContent(eta_i+1, dNdeta_1D_fulleta_MC[Mbin]->GetBinContent( dNdeta_1D_fulleta_MC[Mbin] -> FindBin(dNdeta_1D_MC[Mbin] -> GetBinCenter(eta_i+1)) ));   
0603             for (int z_bin : fulleta_MC_included_z_bin){
0604                 //todo: test
0605                 N_event_counting_MC[Mbin][eta_i] += centrality_Z_map_MC -> GetBinContent(Mbin + 1, z_bin);
0606             }
0607         }
0608 
0609         // note : only for the data, now 
0610         for (auto eta_z : included_eta_z_map) // note : eta bin
0611         {
0612             DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] = new THStack(Form("DeltaPhi_Multi_Stack_MBin%i_Eta%i",Mbin,eta_z.first),Form("DeltaPhi_Multi_Stack_MBin%i_Eta%i;#Delta#phi (Inner - Outer) [degree];Entry",Mbin,eta_z.first));
0613             for(int zbin : eta_z.second) // note :  Z bin
0614             {
0615                 // note : the key to the map: "MBin_EtaBin"
0616                 DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] -> Add(DeltaPhi_Multi_1D[Form("%i_%i_%i",Mbin,eta_z.first,zbin)]);
0617                 SignalNTrack_Single[Form("%i_%i",Mbin,eta_z.first)] += SignalNTrack_eta_z_Single_2D[Mbin] -> GetBinContent(eta_z.first, zbin); // note : the signal counting for the single method
0618                 SignalNTrack_Multi[Form("%i_%i",Mbin,eta_z.first)]  +=  SignalNTrack_eta_z_Multi_2D[Mbin] -> GetBinContent(eta_z.first, zbin); // note : the signal counting for the multi method
0619 
0620                 // note : to count the TrueNTtrack in the true level information
0621                 //todo: test
0622                 // dNdeta_1D_MC[Mbin]->SetBinContent(eta_z.first - eta_correction, TrueNtrack_eta_z_MC_2D[Mbin] -> GetBinContent(eta_z.first, zbin) + dNdeta_1D_MC[Mbin]->GetBinContent(eta_z.first - eta_correction));
0623 
0624                 N_event_counting[Mbin][eta_z.first - eta_correction - 1]    += centrality_Z_map -> GetBinContent(Mbin + 1, zbin);
0625                 //todo: test
0626                 // N_event_counting_MC[Mbin][eta_z.first - eta_correction - 1] += centrality_Z_map_MC -> GetBinContent(Mbin + 1, zbin);
0627             }
0628 
0629             // note : Background fit on the DeltaPhi distributions, to know the background level
0630             // note : print the DeltaPhi distributions for each centrality bin adn each eta bin, after stacking
0631             temp_hist = (TH1F *) DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] -> GetStack() -> Last();
0632             
0633             DeltaPhi_Multi_Stack_hist_out.push_back(temp_hist);
0634             DeltaPhi_Multi_Stack_hist_out[DeltaPhi_Multi_Stack_hist_out.size() - 1] -> SetTitle(Form("MBin%i_EtaBin%i",Mbin,eta_z.first));
0635 
0636             double hist_offset = get_dist_offset(temp_hist, 15);
0637 
0638             bkg_fit_pol2 -> SetParameters(hist_offset, 0, -0.2, 0, signal_region);
0639             bkg_fit_pol2 -> FixParameter(4, signal_region);
0640             bkg_fit_pol2 -> SetParLimits(2, -100, 0);
0641             temp_hist -> Fit(bkg_fit_pol2,"NQ");
0642             // note : extract the background region (which includes the signal region also)
0643             draw_pol2_line -> SetParameters(bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3));
0644 
0645             DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] -> SetMinimum(0);
0646             DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] -> SetMaximum( temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0647             DeltaPhi_Multi_Stack[Form("%i_%i",Mbin,eta_z.first)] -> Draw(); 
0648 
0649             draw_pol2_line -> Draw("lsame");
0650             
0651             double pol2_bkg_integral = fabs(draw_pol2_line -> Integral( -1. * signal_region, signal_region )) / temp_hist -> GetBinWidth(1);
0652             // cout<<i<<" "<<i1<<" pol2_bkg integral: "<<pol2_bkg_integral<<endl;
0653 
0654             // note : for the multi method, count the number of tracklets after the back subtraction
0655             dNdeta_1D_reco_multi_original[Mbin]  -> SetBinContent(eta_z.first - eta_correction, SignalNTrack_Multi[Form("%i_%i",Mbin,eta_z.first)] - pol2_bkg_integral);
0656             dNdeta_1D_reco_single_original[Mbin] -> SetBinContent(eta_z.first - eta_correction, SignalNTrack_Single[Form("%i_%i",Mbin,eta_z.first)]);
0657 
0658             coord_line -> DrawLine(-1*signal_region, 0, -1 * signal_region, temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0659             coord_line -> DrawLine(signal_region, 0, signal_region, temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0660             
0661             ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0662             
0663             draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, #eta: %.2f ~ %.2f",centrality_region[Mbin].c_str(), 
0664                 eta_z_ref -> GetXaxis() -> GetBinCenter(eta_z.first) - eta_z_ref -> GetXaxis() -> GetBinWidth(eta_z.first)/2., 
0665                 eta_z_ref -> GetXaxis() -> GetBinCenter(eta_z.first) + eta_z_ref -> GetXaxis() -> GetBinWidth(eta_z.first)/2.)
0666             );
0667             
0668             draw_text -> DrawLatex(0.21, 0.86, Form("MBin: %i, #eta_bin: %i, #Delta#Phi_bin width: %.2f", Mbin, eta_z.first, temp_hist -> GetBinWidth(1)));
0669             draw_text -> DrawLatex(0.21, 0.82, Form("pol2: %.2f + %.2f(x-%.2f) + %.2f(x-%.2f)^{2}", bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(3), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3)));
0670             draw_text -> DrawLatex(0.21, 0.78, Form("Signal size: %.3f, pol2 bkg size: %.2f", SignalNTrack_Multi[Form("%i_%i",Mbin,eta_z.first)], pol2_bkg_integral));
0671 
0672             c1 -> Print(Form("%s/Reco_DeltaPhi_Multi_Stack.pdf", out_folder_directory.c_str()));
0673             c1 -> Clear();
0674 
0675 
0676             // note : ============================================================================================================================================================================================================
0677             // note : ============================================================================================================================================================================================================
0678             // note : ============================================================================================================================================================================================================
0679             // note : ============================================================================================================================================================================================================
0680             // note : change the unit for the presentation
0681             DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] = new THStack(Form("DeltaPhi_Multi_Stack_MBin%i_Eta%i",Mbin,eta_z.first),Form("DeltaPhi_Multi_Stack_MBin%i_Eta%i;#Delta#phi (Inner - Outer) [radian];Entry",Mbin,eta_z.first));
0682             for(int zbin : eta_z.second) { // note :  Z bin
0683                 // note : the key to the map: "MBin_EtaBin"
0684                 DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] -> Add(DeltaPhi_Multi_1D_radian[Form("%i_%i_%i",Mbin,eta_z.first,zbin)]);
0685             }
0686 
0687             // note : Background fit on the DeltaPhi distributions, to know the background level
0688             // note : print the DeltaPhi distributions for each centrality bin adn each eta bin, after stacking
0689             temp_hist = (TH1F *) DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] -> GetStack() -> Last();
0690             
0691             DeltaPhi_Multi_Stack_hist_out_radian.push_back(temp_hist);
0692             DeltaPhi_Multi_Stack_hist_out_radian[DeltaPhi_Multi_Stack_hist_out_radian.size() - 1] -> SetTitle(Form("MBin%i_EtaBin%i",Mbin,eta_z.first));
0693 
0694             // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2
0695             draw_pol2_line -> SetParameters(bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1) * radian_conversion, bkg_fit_pol2 -> GetParameter(2) * pow(radian_conversion,2), bkg_fit_pol2 -> GetParameter(3)/radian_conversion);
0696 
0697             DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] -> SetMinimum(0);
0698             DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] -> SetMaximum( temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0699             DeltaPhi_Multi_Stack_radian[Form("%i_%i",Mbin,eta_z.first)] -> Draw(); 
0700 
0701             draw_pol2_line -> Draw("lsame");
0702 
0703             coord_line -> DrawLine(-1 * signal_region / radian_conversion, 0, -1 * signal_region / radian_conversion, temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0704             coord_line -> DrawLine(     signal_region / radian_conversion, 0,      signal_region / radian_conversion, temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
0705             
0706             ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0707             
0708             draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, #eta: %.2f ~ %.2f",centrality_region[Mbin].c_str(), 
0709                 eta_z_ref -> GetXaxis() -> GetBinCenter(eta_z.first) - eta_z_ref -> GetXaxis() -> GetBinWidth(eta_z.first)/2., 
0710                 eta_z_ref -> GetXaxis() -> GetBinCenter(eta_z.first) + eta_z_ref -> GetXaxis() -> GetBinWidth(eta_z.first)/2.)
0711             );
0712             
0713             draw_text -> DrawLatex(0.21, 0.86, Form("#Delta#Phi_bin width: %.5f", temp_hist -> GetBinWidth(1)));
0714             draw_text -> DrawLatex(0.21, 0.82, Form("pol2: %.2f + %.2f(x-%.2f) + %.2f(x-%.2f)^{2}", draw_pol2_line -> GetParameter(0), draw_pol2_line -> GetParameter(1), draw_pol2_line -> GetParameter(3), draw_pol2_line -> GetParameter(2), draw_pol2_line -> GetParameter(3)));
0715 
0716             c1 -> Print(Form("%s/Reco_DeltaPhi_Multi_Stack.pdf", out_folder_directory.c_str()));
0717             c1 -> Clear();
0718         }
0719     }
0720     c1 -> Print(Form("%s/Reco_DeltaPhi_Multi_Stack.pdf)", out_folder_directory.c_str()));
0721     c1 -> Clear();
0722 
0723     cout<<"MainPreparation done"<<endl;
0724     return;
0725 }
0726 
0727 void EtaDistReader::HistDivision()
0728 {
0729     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0730     {   
0731         dNdeta_1D_MC[Mbin]          -> Scale(1./double(dNdeta_1D_MC[Mbin]          -> GetBinWidth(1)));
0732         dNdeta_1D_fulleta_MC[Mbin]  -> Scale(1./double(dNdeta_1D_fulleta_MC[Mbin]  -> GetBinWidth(1)));
0733         dNdeta_1D_reco_single_original[Mbin] -> Scale(1./double(dNdeta_1D_reco_single_original[Mbin] -> GetBinWidth(1)));
0734         dNdeta_1D_reco_multi_original[Mbin]  -> Scale(1./double(dNdeta_1D_reco_multi_original[Mbin]  -> GetBinWidth(1)));
0735 
0736         // note : only for the fulleta region
0737         for (int i = 0; i < dNdeta_1D_fulleta_MC[Mbin] -> GetNbinsX(); i++)
0738         {
0739             dNdeta_1D_fulleta_MC[Mbin] -> SetBinContent(i+1, dNdeta_1D_fulleta_MC[Mbin] -> GetBinContent(i+1) / double(N_event_counting_fulleta_MC[Mbin][i]));
0740             dNdeta_1D_fulleta_MC[Mbin] -> SetBinError(i+1, dNdeta_1D_fulleta_MC[Mbin] -> GetBinError(i+1) / double(N_event_counting_fulleta_MC[Mbin][i]));
0741         }
0742 
0743         for (int i = 0; i < included_eta_z_map.size(); i++)
0744         {
0745             dNdeta_1D_MC[Mbin] -> SetBinContent(i+1, dNdeta_1D_MC[Mbin] -> GetBinContent(i+1) / double(N_event_counting_MC[Mbin][i]));
0746             dNdeta_1D_MC[Mbin] -> SetBinError(i+1, dNdeta_1D_MC[Mbin] -> GetBinError(i+1) / double(N_event_counting_MC[Mbin][i]));
0747 
0748             dNdeta_1D_reco_single_original[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1) / double(N_event_counting[Mbin][i]));
0749             dNdeta_1D_reco_single_original[Mbin] -> SetBinError(i+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1) / double(N_event_counting[Mbin][i]));
0750 
0751             dNdeta_1D_reco_multi_original[Mbin]  -> SetBinContent(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1) / double(N_event_counting[Mbin][i]));
0752             dNdeta_1D_reco_multi_original[Mbin]  -> SetBinError(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) / double(N_event_counting[Mbin][i]));
0753 
0754             cout<<"check denominators, Mbin : "<<Mbin<<" eta bin : "<<i+1<<" truth_count: "<<N_event_counting_MC[Mbin][i]<<" reco_count: "<<N_event_counting[Mbin][i]<<endl;
0755 
0756         }
0757 
0758         //note : clone the histogram for the alpha correction
0759         dNdeta_1D_reco_multi_postalpha[Mbin] = (TH1F *) dNdeta_1D_reco_multi_original[Mbin] -> Clone(Form("dNdeta_1D_reco_multi_postalpha_Mbin%i",Mbin));
0760         dNdeta_1D_reco_single_postalpha[Mbin] = (TH1F *) dNdeta_1D_reco_single_original[Mbin] -> Clone(Form("dNdeta_1D_reco_single_postalpha_Mbin%i",Mbin));
0761     }
0762 
0763 
0764 }
0765 
0766 // note : correction_step 0: correction b/w the original distribution and true dist. ---> this is for deriving 
0767 // note : correction_step 1: b/w post_alpha dist. and true dist. ---> this is for the test
0768 void EtaDistReader::DeriveAlphaCorrection(int correction_step)
0769 {
0770     double MC_hist_counting = 0;
0771     double data_tight_counting = 0;
0772     double data_loose_counting = 0;
0773 
0774     vector<TH1F *> temp_hist_single_vec;
0775     vector<TH1F *> temp_hist_multi_vec;
0776 
0777     if (correction_step == 0) {
0778         temp_hist_single_vec = dNdeta_1D_reco_single_original;
0779         temp_hist_multi_vec = dNdeta_1D_reco_multi_original;
0780     }
0781     else if (correction_step == 1){
0782         temp_hist_single_vec = dNdeta_1D_reco_single_postalpha;
0783         temp_hist_multi_vec = dNdeta_1D_reco_multi_postalpha;
0784     }
0785     else {
0786         std::cout<<"wrong correction_step input !!!"<<std::endl;
0787         return;
0788     }
0789 
0790     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0791     {   
0792         std::cout<<" "<<std::endl;
0793         
0794         cout<<"----------- for the case of method tight ----------------" <<endl;
0795         // note : to print the ratio between reco track and MC track
0796         std::cout << "centrality bin : "<<Mbin<<", ";
0797         // note : the eta bin
0798         for (int bin_i = 0; bin_i < dNdeta_1D_MC[Mbin] -> GetNbinsX(); bin_i++) {
0799             MC_hist_counting += dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1);
0800             data_tight_counting += temp_hist_single_vec[Mbin] -> GetBinContent(bin_i+1);
0801             std::cout <<"~~"<<temp_hist_single_vec[Mbin] -> GetBinContent(bin_i+1) <<", "<< dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)<< ", " << Form("%.3f",temp_hist_single_vec[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)) <<"~~, ";
0802 
0803             eta_Mbin_correction_tight -> SetBinContent(bin_i + 1, Mbin+1, temp_hist_single_vec[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
0804         }
0805         std::cout << std::endl;
0806         cout<<"Mbin : "<<Mbin<<" MC_hist_counting: "<<MC_hist_counting<<" data_tight_counting: "<<data_tight_counting<<endl;
0807 
0808         MC_hist_counting = 0;
0809         data_tight_counting = 0;
0810     }
0811 
0812     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0813     {
0814         cout<<"----------- for the case of method inclusive ----------------" <<endl;
0815         // note : to print the ratio between reco track and MC track
0816         std::cout << "centrality bin : "<<Mbin<<", ";
0817         for (int bin_i = 0; bin_i < dNdeta_1D_MC[Mbin] -> GetNbinsX(); bin_i++) {
0818             MC_hist_counting += dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1);
0819             data_loose_counting += temp_hist_multi_vec[Mbin] -> GetBinContent(bin_i+1);
0820             std::cout <<"~~"<<temp_hist_multi_vec[Mbin] -> GetBinContent(bin_i+1) <<", "<< dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)<< ", " << Form( "%.3f", temp_hist_multi_vec[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)) <<"~~, ";
0821 
0822             eta_Mbin_correction_loose       -> SetBinContent(bin_i + 1, Mbin + 1, temp_hist_multi_vec[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
0823             eta_Mbin_correction_loose_noUPC -> SetBinContent(bin_i + 1, Mbin + 1, temp_hist_multi_vec[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
0824         }
0825         std::cout << std::endl;
0826         cout<<"Mbin : "<<Mbin<<" MC_hist_counting: "<<MC_hist_counting<<" data_loose_counting: "<<data_loose_counting<<endl;
0827 
0828         MC_hist_counting = 0;
0829         data_loose_counting = 0;
0830     }
0831 }
0832 
0833 vector<TH2F *> EtaDistReader::Get_alpha_corr_map() {return {eta_Mbin_correction_loose, eta_Mbin_correction_tight};}
0834 
0835 void EtaDistReader::Set_alpha_corr_map(TH2F * hist_in_loose, TH2F * hist_in_tight)
0836 {
0837     input_eta_Mbin_correction_loose = hist_in_loose;
0838     input_eta_Mbin_correction_tight = hist_in_tight;
0839 
0840     std::cout<<"in EtaDistReader, input map loose check bin (3,3): "<<input_eta_Mbin_correction_loose -> GetBinContent(3,3)<<endl;
0841     std::cout<<"in EtaDistReader, input map tight check bin (3,3): "<<input_eta_Mbin_correction_tight -> GetBinContent(3,3)<<endl;
0842 
0843     return;
0844 }
0845 
0846 void EtaDistReader::ApplyAlphaCorrection()
0847 {
0848     if (input_eta_Mbin_correction_loose == nullptr || input_eta_Mbin_correction_tight == nullptr)
0849     {
0850         std::cout<<" no alpha correction map included !!! "<<endl;
0851         return;
0852     }
0853 
0854     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0855     {   
0856         for (int i = 0; i < included_eta_z_map.size(); i++)
0857         {        
0858             // note : for the loose method 
0859             dNdeta_1D_reco_multi_postalpha[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_multi_postalpha[Mbin] -> GetBinContent(i+1) / input_eta_Mbin_correction_loose -> GetBinContent(i+1, Mbin+1));
0860             dNdeta_1D_reco_multi_postalpha[Mbin] -> SetBinError(i+1,   dNdeta_1D_reco_multi_postalpha[Mbin] -> GetBinError(i+1)   /  input_eta_Mbin_correction_loose -> GetBinContent(i+1, Mbin+1));
0861             
0862             // note : for the tight method 
0863             dNdeta_1D_reco_single_postalpha[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_single_postalpha[Mbin] -> GetBinContent(i+1) / input_eta_Mbin_correction_tight -> GetBinContent(i+1, Mbin+1));
0864             dNdeta_1D_reco_single_postalpha[Mbin] -> SetBinError(i+1,   dNdeta_1D_reco_single_postalpha[Mbin] -> GetBinError(i+1)   / input_eta_Mbin_correction_tight -> GetBinContent(i+1, Mbin+1));
0865         }
0866     }
0867 
0868     post_alpha_tag = 1;
0869 
0870     std::cout<<"finish the alpha correction application"<<std::endl;
0871 
0872     return;
0873 
0874 }
0875 
0876 void EtaDistReader::PrintPlots_exclusive()
0877 {
0878     c1 ->cd();
0879     exclusive_cluster_inner_adc -> SetMinimum(0);
0880     exclusive_cluster_inner_adc -> Draw("hist");
0881     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0882     c1 -> Print(Form("%s/exclusive_cluster_inner_adc.pdf", out_folder_directory.c_str()));
0883     c1 -> Clear();
0884 
0885     c1 ->cd();
0886     exclusive_cluster_outer_adc -> SetMinimum(0);
0887     exclusive_cluster_outer_adc -> Draw("hist");
0888     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0889     c1 -> Print(Form("%s/exclusive_cluster_outer_adc.pdf", out_folder_directory.c_str()));
0890     c1 -> Clear();
0891 
0892     c1 ->cd();
0893     exclusive_cluster_inner_eta_phi_2D -> Draw("colz0");
0894     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0895     c1 -> Print(Form("%s/exclusive_cluster_inner_eta_phi_2D.pdf", out_folder_directory.c_str()));
0896     c1 -> Clear();
0897 
0898     c1 ->cd();
0899     exclusive_cluster_outer_eta_phi_2D -> Draw("colz0");
0900     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0901     c1 -> Print(Form("%s/exclusive_cluster_outer_eta_phi_2D.pdf", out_folder_directory.c_str()));
0902     c1 -> Clear();
0903 
0904     c1 ->cd();
0905     exclusive_cluster_inner_eta_adc_2D -> Draw("colz0");
0906     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0907     c1 -> Print(Form("%s/exclusive_cluster_inner_eta_adc_2D.pdf", out_folder_directory.c_str()));
0908     c1 -> Clear();
0909 
0910     c1 ->cd();
0911     exclusive_cluster_outer_eta_adc_2D -> Draw("colz0");
0912     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0913     c1 -> Print(Form("%s/exclusive_cluster_outer_eta_adc_2D.pdf", out_folder_directory.c_str()));
0914     c1 -> Clear();
0915 
0916     c1 ->cd();
0917     exclusive_cluster_inner_eta_adc_2D_bkgrm -> Draw("colz0");
0918     exclusive_cluster_inner_eta_adc_2D_bkgrm_profile_graph -> Draw("p same");
0919     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0920     c1 -> Print(Form("%s/exclusive_cluster_inner_eta_adc_2D_bkgrm.pdf", out_folder_directory.c_str()));
0921     c1 -> Clear();
0922 
0923     c1 ->cd();
0924     exclusive_cluster_outer_eta_adc_2D_bkgrm -> Draw("colz0");
0925     exclusive_cluster_outer_eta_adc_2D_bkgrm_profile_graph -> Draw("p same");
0926     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0927     c1 -> Print(Form("%s/exclusive_cluster_outer_eta_adc_2D_bkgrm.pdf", out_folder_directory.c_str()));
0928     c1 -> Clear();
0929 
0930 }
0931 
0932 void EtaDistReader::PrintPlots()
0933 {
0934     gStyle->SetPaintTextFormat("1.3f");
0935 
0936     c1 -> cd();
0937     eta_Mbin_correction_loose -> Draw("colz0");
0938     eta_Mbin_correction_loose -> SetMarkerSize(0.7);
0939     eta_Mbin_correction_loose -> Draw("HIST TEXT45 SAME");
0940     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0941     c1 -> Print(Form("%s/eta_Mbin_correction_loose.pdf", out_folder_directory.c_str()));
0942     c1 -> Clear();
0943 
0944     c1 -> cd();
0945     for (int i = 0; i < eta_Mbin_correction_loose_noUPC -> GetNbinsX(); i++) {eta_Mbin_correction_loose_noUPC -> SetBinContent(i+1, eta_Mbin_correction_loose_noUPC -> GetNbinsY()-1, 0);}
0946     eta_Mbin_correction_loose_noUPC -> Draw("colz0");
0947     eta_Mbin_correction_loose_noUPC -> SetMarkerSize(0.7);
0948     eta_Mbin_correction_loose_noUPC -> Draw("HIST TEXT45 SAME");
0949     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0950     c1 -> Print(Form("%s/eta_Mbin_correction_loose_noUPC.pdf", out_folder_directory.c_str()));
0951     c1 -> Clear();
0952 
0953     c1 -> cd();
0954     eta_Mbin_correction_tight -> Draw("colz0");
0955     eta_Mbin_correction_tight -> SetMarkerSize(0.7);
0956     eta_Mbin_correction_tight -> Draw("HIST TEXT45 SAME");
0957     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0958     c1 -> Print(Form("%s/eta_Mbin_correction_tight.pdf", out_folder_directory.c_str()));
0959     c1 -> Clear();
0960 
0961 
0962     c1 -> cd();
0963     c1 -> Print(Form("%s/dNdeta_combine_final_original.pdf(", out_folder_directory.c_str()));
0964     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
0965     {   
0966         if (Mbin == 0)
0967         {
0968             legend -> AddEntry(dNdeta_1D_fulleta_MC[Mbin], "Truth dist. (full)","f");
0969             legend -> AddEntry(dNdeta_1D_MC[Mbin], "Truth dist. (selected)","f");
0970             legend -> AddEntry(dNdeta_1D_reco_single_original[Mbin], "Reco. Method tight","lep");
0971             legend -> AddEntry(dNdeta_1D_reco_multi_original[Mbin], "Reco. Method loose","lep");
0972         }
0973 
0974         // todo: this is for the test, to full the fulleta dNdeta together with others
0975         dNdeta_1D_fulleta_MC[Mbin] -> SetMinimum(0);
0976         dNdeta_1D_fulleta_MC[Mbin] -> SetMaximum( dNdeta_1D_fulleta_MC[Mbin] -> GetBinContent( dNdeta_1D_fulleta_MC[Mbin] -> GetMaximumBin() ) * 1.5);
0977         dNdeta_1D_fulleta_MC[Mbin] -> Draw("hist");
0978 
0979         // dNdeta_1D_MC[Mbin] -> GetYaxis() -> SetRangeUser(0, dNdeta_1D_MC[Mbin] -> GetMaximum() * 1.5);
0980 
0981         dNdeta_1D_MC[Mbin] -> Draw("hist same");
0982         dNdeta_1D_reco_single_original[Mbin] -> Draw("p same");
0983         dNdeta_1D_reco_multi_original[Mbin]  -> Draw("p same");
0984 
0985         legend -> Draw("same");
0986 
0987         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0988         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[Mbin].c_str()));
0989        
0990         c1 -> Print(Form("%s/dNdeta_combine_final_original.pdf", out_folder_directory.c_str()));
0991         c1 -> Clear();
0992     }
0993     c1 -> Print(Form("%s/dNdeta_combine_final_original.pdf)", out_folder_directory.c_str()));
0994     c1 -> Clear();
0995     legend -> Clear();
0996 
0997     
0998 
0999     c1 -> cd();
1000     c1 -> Print(Form("%s/dNdeta_combine_final_postalpha.pdf(", out_folder_directory.c_str()));
1001     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1002     {   
1003         if (Mbin == 0)
1004         {
1005             legend -> AddEntry(dNdeta_1D_fulleta_MC[Mbin], "Truth dist. (full)","f");
1006             legend -> AddEntry(dNdeta_1D_MC[Mbin], "Truth dist. (selected)","f");
1007             legend -> AddEntry(dNdeta_1D_reco_single_postalpha[Mbin], "Reco. Method tight","lep");
1008             legend -> AddEntry(dNdeta_1D_reco_multi_postalpha[Mbin], "Reco. Method loose","lep");
1009         }
1010 
1011         // todo: this is for the test, to full the fulleta dNdeta together with others
1012         dNdeta_1D_fulleta_MC[Mbin] -> SetMinimum(0);
1013         dNdeta_1D_fulleta_MC[Mbin] -> SetMaximum( dNdeta_1D_fulleta_MC[Mbin] -> GetBinContent( dNdeta_1D_fulleta_MC[Mbin] -> GetMaximumBin() ) * 1.5);
1014         dNdeta_1D_fulleta_MC[Mbin] -> Draw("hist");
1015 
1016         // dNdeta_1D_MC[Mbin] -> GetYaxis() -> SetRangeUser(0, dNdeta_1D_MC[Mbin] -> GetMaximum() * 1.5);
1017 
1018         dNdeta_1D_MC[Mbin] -> Draw("hist same");
1019         Hist1DSetting(dNdeta_1D_reco_single_postalpha[Mbin], "#1A3947");
1020         dNdeta_1D_reco_single_postalpha[Mbin] -> Draw("p same");
1021         Hist1DSetting(dNdeta_1D_reco_multi_postalpha[Mbin], "#c48045");
1022         dNdeta_1D_reco_multi_postalpha[Mbin]  -> Draw("p same");
1023 
1024         legend -> Draw("same");
1025 
1026         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1027         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[Mbin].c_str()));
1028        
1029         c1 -> Print(Form("%s/dNdeta_combine_final_postalpha.pdf", out_folder_directory.c_str()));
1030         c1 -> Clear();
1031     }
1032     c1 -> Print(Form("%s/dNdeta_combine_final_postalpha.pdf)", out_folder_directory.c_str()));
1033     c1 -> Clear();
1034     legend -> Clear();
1035 
1036 
1037 }
1038 
1039 void EtaDistReader::EndRun()
1040 {
1041     TFile * file_out = new TFile(Form("%s/alpha_correction_map.root",out_folder_directory.c_str()),"RECREATE");
1042     // file_out -> cd();
1043     eta_Mbin_correction_loose -> Write("eta_Mbin_correction_loose");
1044     eta_Mbin_correction_loose_noUPC -> Write("eta_Mbin_correction_loose_noUPC");
1045     eta_Mbin_correction_tight -> Write("eta_Mbin_correction_tight");
1046     
1047     // note : the dNdeta 1D for fulleta region, MC -> this is truth 
1048     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1049     {
1050         dNdeta_1D_fulleta_MC[Mbin] -> Write(Form("dNdeta_1D_fulleta_MC_Mbin%i",Mbin));
1051     }
1052     
1053     // note : the dNdeta 1D for selected eta region, MC
1054     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1055     {
1056         dNdeta_1D_MC[Mbin] -> Write(Form("dNdeta_1D_MC_Mbin%i",Mbin));
1057     }
1058 
1059     // note : the reconstructed dNdeta
1060     // note : loose method
1061     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1062     {
1063         dNdeta_1D_reco_multi_original[Mbin] -> Write(Form("dNdeta_1D_reco_loose_original_Mbin%i", Mbin));
1064     }
1065 
1066     // note : the reconstructed dNdeta 
1067     // note : tight method 
1068     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1069     {
1070         dNdeta_1D_reco_single_original[Mbin] -> Write(Form("dNdeta_1D_reco_tight_original_Mbin%i", Mbin));
1071     }
1072 
1073     // note : the reconstructed dNdeta
1074     // note : loose method
1075     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1076     {
1077         dNdeta_1D_reco_multi_postalpha[Mbin] -> Write(Form("dNdeta_1D_reco_loose_postalpha_Mbin%i", Mbin));
1078     }
1079 
1080     // note : the reconstructed dNdeta 
1081     // note : tight method 
1082     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1083     {
1084         dNdeta_1D_reco_single_postalpha[Mbin] -> Write(Form("dNdeta_1D_reco_tight_postalpha_Mbin%i", Mbin));
1085     }
1086 
1087     file_out -> Close();
1088 }
1089 
1090 vector<TH1F *> EtaDistReader::GetDeltaPhi_Multi_stack_1D()
1091 {
1092     return DeltaPhi_Multi_Stack_hist_out;
1093 }
1094 vector<TH1F *> EtaDistReader::GetDeltaPhi_Multi_stack_1D_radian()
1095 {
1096     return DeltaPhi_Multi_Stack_hist_out_radian;
1097 }
1098 vector<TH1F *> EtaDistReader::GetdNdeta_1D_MC()
1099 {
1100     return dNdeta_1D_MC;    
1101 }
1102 
1103 vector<TH1F *> EtaDistReader::GetdNdeta_1D_fulleta_MC()
1104 {
1105     return dNdeta_1D_fulleta_MC;
1106 }
1107 
1108 vector<TH1F *> EtaDistReader::GetdNdeta_1D_reco_single_original()
1109 {
1110     return dNdeta_1D_reco_single_original;    
1111 }
1112 vector<TH1F *> EtaDistReader::GetdNdeta_1D_reco_multi_original()
1113 {
1114     return dNdeta_1D_reco_multi_original;    
1115 }
1116 
1117 vector<TH1F *> EtaDistReader::GetdNdeta_1D_reco_single_postalpha()
1118 {
1119     return dNdeta_1D_reco_single_postalpha;    
1120 }
1121 vector<TH1F *> EtaDistReader::GetdNdeta_1D_reco_multi_postalpha()
1122 {
1123     return dNdeta_1D_reco_multi_postalpha;    
1124 }
1125 
1126 #endif
1127 
1128 
1129 // void EtaDistReader::FinaldNdEta()
1130 // {
1131 //     double MC_hist_counting = 0;
1132 //     double data_tight_counting = 0;
1133 //     double data_loose_counting = 0;
1134 
1135 //     c1 -> cd();
1136 //     c1 -> Print(Form("%s/dNdeta_combine_final_no_correction.pdf(", out_folder_directory.c_str()));
1137 //     for (int Mbin = 0; Mbin < N_centrality_bin; Mbin++)
1138 //     {   
1139 //         // // note : check the ratio of the bin error to the bin content of the three histograms
1140 //         // for (int i = 0; i < included_eta_z_map.size(); i++)
1141 //         // {
1142 //         //     cout<<" "<<endl;
1143 //         //     cout<<"1MBin: "<<Mbin<<" eta bin: "<<i<<" MC: "<<dNdeta_1D_MC[Mbin] -> GetBinError(i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(i+1)<<" reco_single: "<<dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1)<<" reco_multi: "<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1)<<endl;
1144 //         // }
1145 
1146 //         std::cout<<" "<<std::endl;
1147 
1148 //         dNdeta_1D_MC[Mbin]          -> Scale(1./double(dNdeta_1D_MC[Mbin]          -> GetBinWidth(1) ));
1149 //         dNdeta_1D_fulleta_MC[Mbin]  -> Scale(1./double(dNdeta_1D_fulleta_MC[Mbin]  -> GetBinWidth(1) ));
1150 //         dNdeta_1D_reco_single_original[Mbin] -> Scale(1./double(dNdeta_1D_reco_single_original[Mbin] -> GetBinWidth(1) ));
1151 //         dNdeta_1D_reco_multi_original[Mbin]  -> Scale(1./double(dNdeta_1D_reco_multi_original[Mbin]  -> GetBinWidth(1) ));
1152 
1153 //         // // note : check the ratio of the bin error to the bin content of the three histograms
1154 //         // for (int i = 0; i < included_eta_z_map.size(); i++)
1155 //         // {
1156 //         //     cout<<" "<<endl;
1157 //         //     cout<<"2MBin: "<<Mbin<<" eta bin: "<<i<<" MC: "<<dNdeta_1D_MC[Mbin] -> GetBinError(i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(i+1)<<" reco_single: "<<dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1)<<" reco_multi: "<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1)<<endl;
1158 //         // }
1159 
1160 //         // note : only for the fulleta region
1161 //         for (int i = 0; i < dNdeta_1D_fulleta_MC[Mbin] -> GetNbinsX(); i++)
1162 //         {
1163 //             dNdeta_1D_fulleta_MC[Mbin] -> SetBinContent(i+1, dNdeta_1D_fulleta_MC[Mbin] -> GetBinContent(i+1) / double(N_event_counting_fulleta_MC[Mbin][i]));
1164 //             dNdeta_1D_fulleta_MC[Mbin] -> SetBinError(i+1, dNdeta_1D_fulleta_MC[Mbin] -> GetBinError(i+1) / double(N_event_counting_fulleta_MC[Mbin][i]));
1165 //         }
1166 
1167 //         for (int i = 0; i < included_eta_z_map.size(); i++)
1168 //         {
1169 //             dNdeta_1D_MC[Mbin] -> SetBinContent(i+1, dNdeta_1D_MC[Mbin] -> GetBinContent(i+1) / double(N_event_counting_MC[Mbin][i]));
1170 //             dNdeta_1D_MC[Mbin] -> SetBinError(i+1, dNdeta_1D_MC[Mbin] -> GetBinError(i+1) / double(N_event_counting_MC[Mbin][i]));
1171 
1172 //             dNdeta_1D_reco_single_original[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1) / double(N_event_counting[Mbin][i]));
1173 //             dNdeta_1D_reco_single_original[Mbin] -> SetBinError(i+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1) / double(N_event_counting[Mbin][i]));
1174 
1175 //             dNdeta_1D_reco_multi_original[Mbin]  -> SetBinContent(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1) / double(N_event_counting[Mbin][i]));
1176 //             dNdeta_1D_reco_multi_original[Mbin]  -> SetBinError(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) / double(N_event_counting[Mbin][i]));
1177 
1178 //             cout<<"check denominators, Mbin : "<<Mbin<<" eta bin : "<<i+1<<" MC_count: "<<N_event_counting_MC[Mbin][i]<<" data_count: "<<N_event_counting[Mbin][i]<<endl;
1179         
1180 //             if (input_eta_Mbin_correction_loose != nullptr)
1181 //             {
1182 //                 dNdeta_1D_reco_multi_original[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1) / input_eta_Mbin_correction_loose -> GetBinContent(i+1, Mbin+1));
1183 //                 dNdeta_1D_reco_multi_original[Mbin] -> SetBinError(i+1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) /     input_eta_Mbin_correction_loose -> GetBinContent(i+1, Mbin+1));
1184 //             }
1185 
1186 //             if (input_eta_Mbin_correction_tight != nullptr)
1187 //             {
1188 //                 dNdeta_1D_reco_single_original[Mbin] -> SetBinContent(i+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1) / input_eta_Mbin_correction_tight -> GetBinContent(i+1, Mbin+1));
1189 //                 dNdeta_1D_reco_single_original[Mbin] -> SetBinError(i+1,   dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1)   / input_eta_Mbin_correction_tight -> GetBinContent(i+1, Mbin+1));
1190 //             }
1191 
1192 //         }
1193         
1194 //         cout<<"----------- for the case of method tight ----------------" <<endl;
1195 //         // note : to print the ratio between reco track and MC track
1196 //         std::cout << "centrality bin : "<<Mbin<<", ";
1197 //         for (int bin_i = 0; bin_i < dNdeta_1D_MC[Mbin] -> GetNbinsX(); bin_i++) {
1198 //             MC_hist_counting += dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1);
1199 //             data_tight_counting += dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(bin_i+1);
1200 //             std::cout <<"~~"<<dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(bin_i+1) <<", "<< dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)<< ", " << Form("%.3f",dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)) <<"~~, ";
1201 
1202 //             eta_Mbin_correction_tight -> SetBinContent(bin_i + 1, Mbin+1, dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
1203 //         }
1204 //         std::cout << std::endl;
1205 
1206 //         cout<<"----------- for the case of method inclusive ----------------" <<endl;
1207 //         // note : to print the ratio between reco track and MC track
1208 //         std::cout << "centrality bin : "<<Mbin<<", ";
1209 //         for (int bin_i = 0; bin_i < dNdeta_1D_MC[Mbin] -> GetNbinsX(); bin_i++) {
1210 //             data_loose_counting += dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(bin_i+1);
1211 //             std::cout <<"~~"<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(bin_i+1) <<", "<< dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)<< ", " << Form( "%.3f", dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1)) <<"~~, ";
1212 
1213 //             eta_Mbin_correction_loose       -> SetBinContent(bin_i + 1, Mbin + 1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
1214 //             eta_Mbin_correction_loose_noUPC -> SetBinContent(bin_i + 1, Mbin + 1, dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(bin_i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(bin_i+1));
1215 //         }
1216 //         std::cout << std::endl;
1217 
1218 //         // note : check the bin content of the three histograms
1219 //         // note : and check the bin error of the three histograms
1220 //         // for (int i = 0; i < included_eta_z_map.size(); i++)
1221 //         // {
1222 //         //     cout<<" "<<endl;
1223 //         //     cout<<"3MBin: "<<Mbin<<" eta bin: "<<i<<" MC: "<<dNdeta_1D_MC[Mbin] -> GetBinContent(i+1)<<" reco_single: "<<dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1)<<" reco_multi: "<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1)<<endl;
1224 //         //     cout<<"3MBin: "<<Mbin<<" eta bin: "<<i<<" MC: "<<dNdeta_1D_MC[Mbin] -> GetBinError(i+1)<<" reco_single: "<<dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1)<<" reco_multi: "<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1)<<endl;
1225 //         // }
1226 
1227 //         // // note : check the ratio of the bin error to the bin content of the three histograms
1228 //         // for (int i = 0; i < included_eta_z_map.size(); i++)
1229 //         // {
1230 //         //     cout<<" "<<endl;
1231 //         //     cout<<"3MBin: "<<Mbin<<" eta bin: "<<i<<" MC: "<<dNdeta_1D_MC[Mbin] -> GetBinError(i+1) / dNdeta_1D_MC[Mbin] -> GetBinContent(i+1)<<" reco_single: "<<dNdeta_1D_reco_single_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_single_original[Mbin] -> GetBinContent(i+1)<<" reco_multi: "<<dNdeta_1D_reco_multi_original[Mbin] -> GetBinError(i+1) / dNdeta_1D_reco_multi_original[Mbin] -> GetBinContent(i+1)<<endl;
1232 //         // }
1233 
1234 //         if (Mbin == 0)
1235 //         {
1236 //             legend -> AddEntry(dNdeta_1D_fulleta_MC[Mbin], "Truth dist. (full)","f");
1237 //             legend -> AddEntry(dNdeta_1D_MC[Mbin], "Truth dist. (selected)","f");
1238 //             legend -> AddEntry(dNdeta_1D_reco_single_original[Mbin], "Reco. Method tight","lep");
1239 //             legend -> AddEntry(dNdeta_1D_reco_multi_original[Mbin], "Reco. Method loose","lep");
1240 //         }
1241 
1242 //         // todo: this is for the test, to full the fulleta dNdeta together with others
1243 //         dNdeta_1D_fulleta_MC[Mbin] -> GetYaxis() -> SetRangeUser(0, dNdeta_1D_MC[Mbin] -> GetMaximum() * 1.5);
1244 //         dNdeta_1D_fulleta_MC[Mbin] -> Draw("hist");
1245 
1246 //         // dNdeta_1D_MC[Mbin] -> GetYaxis() -> SetRangeUser(0, dNdeta_1D_MC[Mbin] -> GetMaximum() * 1.5);
1247 
1248 //         dNdeta_1D_MC[Mbin] -> Draw("hist same");
1249 //         dNdeta_1D_reco_single_original[Mbin] -> Draw("p same");
1250 //         dNdeta_1D_reco_multi_original[Mbin]  -> Draw("p same");
1251 
1252 //         legend -> Draw("same");
1253 
1254 //         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1255 //         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[Mbin].c_str()));
1256        
1257 //         c1 -> Print(Form("%s/dNdeta_combine_final_no_correction.pdf", out_folder_directory.c_str()));
1258 //         c1 -> Clear();
1259 
1260 //         cout<<"Mbin : "<<Mbin<<" MC_hist_counting: "<<MC_hist_counting<<" data_tight_counting: "<<data_tight_counting<<" data_loose_counting: "<<data_loose_counting<<endl;
1261 
1262 //         MC_hist_counting = 0;
1263 //         data_tight_counting = 0;
1264 //         data_loose_counting = 0;
1265 //     }
1266 //     c1 -> Print(Form("%s/dNdeta_combine_final_no_correction.pdf)", out_folder_directory.c_str()));
1267 //     c1 -> Clear();
1268 
1269 //     gStyle->SetPaintTextFormat("1.3f");
1270 
1271 //     c1 -> cd();
1272 //     eta_Mbin_correction_loose -> Draw("colz0");
1273 //     eta_Mbin_correction_loose -> SetMarkerSize(0.7);
1274 //     eta_Mbin_correction_loose -> Draw("HIST TEXT45 SAME");
1275 //     c1 -> Print(Form("%s/eta_Mbin_correction_loose.pdf", out_folder_directory.c_str()));
1276 //     c1 -> Clear();
1277 
1278 //     c1 -> cd();
1279 //     for (int i = 0; i < eta_Mbin_correction_loose_noUPC -> GetNbinsX(); i++) {eta_Mbin_correction_loose_noUPC -> SetBinContent(i+1, eta_Mbin_correction_loose_noUPC -> GetNbinsY()-1, 0);}
1280 //     eta_Mbin_correction_loose_noUPC -> Draw("colz0");
1281 //     eta_Mbin_correction_loose_noUPC -> SetMarkerSize(0.7);
1282 //     eta_Mbin_correction_loose_noUPC -> Draw("HIST TEXT45 SAME");
1283 //     c1 -> Print(Form("%s/eta_Mbin_correction_loose_noUPC.pdf", out_folder_directory.c_str()));
1284 //     c1 -> Clear();
1285 
1286 //     c1 -> cd();
1287 //     eta_Mbin_correction_tight -> Draw("colz0");
1288 //     eta_Mbin_correction_tight -> SetMarkerSize(0.7);
1289 //     eta_Mbin_correction_tight -> Draw("HIST TEXT45 SAME");
1290 //     c1 -> Print(Form("%s/eta_Mbin_correction_tight.pdf", out_folder_directory.c_str()));
1291 //     c1 -> Clear();
1292 
1293 //     cout<<"FinaldNdEta done"<<endl;
1294 //     return;
1295 // }