Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-09 08:12:17

0001 #include "PreparedNdEta.h"
0002 
0003 PreparedNdEta::PreparedNdEta(
0004     int process_id_in,
0005     int runnumber_in,
0006     std::string input_directory_in,
0007     std::string input_file_name_in,
0008     std::string output_directory_in,
0009     std::string output_file_name_suffix_in,
0010     bool ApplyAlphaCorr_in
0011 ) : process_id(process_id_in),
0012     runnumber(runnumber_in),
0013     input_directory(input_directory_in),
0014     input_file_name(input_file_name_in),
0015     output_directory(output_directory_in),
0016     output_file_name_suffix(output_file_name_suffix_in),
0017     ApplyAlphaCorr(ApplyAlphaCorr_in)
0018 {
0019 
0020     PrepareInputRootFie();
0021 
0022     PrepareOutPutFileName();
0023     PrepareOutPutRootFile();
0024     PrepareHistFits();
0025 
0026     vtxZ_index_map = GetVtxZIndexMap();
0027 
0028     c1 = new TCanvas("c1", "c1", 950, 800);
0029 
0030     h2D_alpha_correction_map_in.clear();
0031     h2D_alpha_correction_map_out.clear();
0032 }
0033 
0034 void PreparedNdEta::PrepareInputRootFie()
0035 {
0036     file_in = TFile::Open(Form("%s/%s",input_directory.c_str(),input_file_name.c_str()));
0037     if (file_in == nullptr)
0038     {
0039         std::cout << "Error: file_in can not be opened" << std::endl;
0040         exit(1);
0041     }
0042 
0043     for (TObject* keyAsObj : *file_in->GetListOfKeys())
0044     {
0045         auto key = dynamic_cast<TKey*>(keyAsObj);
0046         std::string hist_name  = key->GetName();
0047         std::string class_name = key->GetClassName();
0048 
0049         if (class_name == "TH2D")
0050         {
0051             h2D_input_map[hist_name.c_str()] = (TH2D*) file_in -> Get( hist_name.c_str() );
0052         }
0053         else if (class_name == "TH1D")
0054         {
0055             h1D_input_map[hist_name.c_str()] = (TH1D*) file_in -> Get( hist_name.c_str() );
0056             h1D_input_map[hist_name.c_str()] -> SetLineColor(kBlack);
0057             h1D_input_map[hist_name.c_str()] -> SetLineWidth(2);
0058         }
0059     }
0060     std::cout<<"In PrepareInputRootFie(), number of input TH1D: "<< h1D_input_map.size() <<std::endl;
0061     std::cout<<"In PrepareInputRootFie(), number of input TH2D: "<< h2D_input_map.size() <<std::endl;
0062 
0063     tree_in = (TTree*) file_in -> Get("tree_par");
0064     
0065     centrality_edges = 0;
0066     cut_GoodProtoTracklet_DeltaPhi = 0;
0067 
0068     tree_in -> SetBranchAddress("centrality_edges", &centrality_edges);
0069     tree_in -> SetBranchAddress("nCentrality_bin", &nCentrality_bin);
0070 
0071     tree_in -> SetBranchAddress("CentralityFineEdge_min", &CentralityFineEdge_min);
0072     tree_in -> SetBranchAddress("CentralityFineEdge_max", &CentralityFineEdge_max);
0073     tree_in -> SetBranchAddress("nCentralityFineBin", &nCentralityFineBin);
0074 
0075     tree_in -> SetBranchAddress("EtaEdge_min", &EtaEdge_min);
0076     tree_in -> SetBranchAddress("EtaEdge_max", &EtaEdge_max);
0077     tree_in -> SetBranchAddress("nEtaBin", &nEtaBin);
0078 
0079     tree_in -> SetBranchAddress("VtxZEdge_min", &VtxZEdge_min);
0080     tree_in -> SetBranchAddress("VtxZEdge_max", &VtxZEdge_max);
0081     tree_in -> SetBranchAddress("nVtxZBin", &nVtxZBin);
0082 
0083     tree_in -> SetBranchAddress("DeltaPhiEdge_min", &DeltaPhiEdge_min);
0084     tree_in -> SetBranchAddress("DeltaPhiEdge_max", &DeltaPhiEdge_max);
0085     tree_in -> SetBranchAddress("nDeltaPhiBin", &nDeltaPhiBin);    
0086 
0087     tree_in -> SetBranchAddress("cut_GoodProtoTracklet_DeltaPhi", &cut_GoodProtoTracklet_DeltaPhi);
0088 
0089     tree_in -> GetEntry(0);
0090 
0091     std::cout<<"------------------------------------------------------------"<<std::endl;
0092     std::cout<<"In PrepareInputRootFile(), nEntries : "<<tree_in->GetEntries()<<std::endl;
0093     std::cout<<"nCentrality_bin : "<<nCentrality_bin<<std::endl;
0094     for(int i = 0; i < nCentrality_bin; i++)
0095     {
0096         std::cout<<"centrality_edges["<<i<<"] : "<<centrality_edges->at(i)<<std::endl;
0097     }
0098     std::cout<<"CentralityFineBin : "<<nCentralityFineBin<<", "<<CentralityFineEdge_min<<", "<<CentralityFineEdge_max<<std::endl;
0099     std::cout<<"EtaBin : "<<nEtaBin<<", "<<EtaEdge_min<<", "<<EtaEdge_max<<std::endl;
0100     std::cout<<"VtxZBin : "<<nVtxZBin<<", "<<VtxZEdge_min<<", "<<VtxZEdge_max<<std::endl;
0101     std::cout<<"DeltaPhiBin : "<<nDeltaPhiBin<<", "<<DeltaPhiEdge_min<<", "<<DeltaPhiEdge_max<<std::endl;
0102     std::cout<<"DeltaPhi signal region: "<<cut_GoodProtoTracklet_DeltaPhi->first<<" ~ "<<cut_GoodProtoTracklet_DeltaPhi->second<<std::endl;
0103     std::cout<<"------------------------------------------------------------"<<std::endl;
0104 
0105 }
0106 
0107 void PreparedNdEta::PrepareOutPutFileName()
0108 {
0109     std::string job_index = std::to_string( process_id );
0110     int job_index_len = 5;
0111     job_index.insert(0, job_index_len - job_index.size(), '0');
0112 
0113     std::string runnumber_str = std::to_string( runnumber );
0114     if (runnumber != -1){
0115         int runnumber_str_len = 8;
0116         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0117     }
0118 
0119     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0120         output_file_name_suffix = "_" + output_file_name_suffix;
0121     }
0122 
0123     output_filename = "PreparedNdEta";
0124     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0125     output_filename += (ApplyAlphaCorr) ? "_ApplyAlphaCorr" : "";
0126     output_filename += output_file_name_suffix;
0127     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0128     // output_filename += (runnumber != -1) ? Form("_%s.root",runnumber_str.c_str()) : Form(".root");
0129 }
0130 
0131 void PreparedNdEta::PrepareOutPutRootFile()
0132 {
0133     file_out = new TFile(Form("%s/%s",output_directory.c_str(),output_filename.c_str()), "RECREATE");
0134     file_out_dNdEta = new TFile(Form("%s/%s_dNdEta.root",output_directory.c_str(),output_filename.c_str()), "RECREATE");
0135 }
0136 
0137 void PreparedNdEta::PrepareTrackletEtaHist(std::map<std::string, TH1D*> &map_in, std::string name_in)
0138 {
0139     map_in.clear();
0140 
0141     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0142         map_in.insert( std::make_pair(
0143                 Form("h1D_Mbin%d", Mbin),
0144                 new TH1D(Form("h1D_%s_Mbin%d", name_in.c_str(), Mbin), Form("h1D_%s_Mbin%d;Tracklet #eta;Entries", name_in.c_str(), Mbin), nEtaBin, EtaEdge_min, EtaEdge_max)
0145             )
0146         ).second;
0147 
0148 
0149         map_in.insert( std::make_pair(
0150                 Form("h1D_typeA_Mbin%d", Mbin),
0151                 new TH1D(Form("h1D_typeA_%s_Mbin%d", name_in.c_str(), Mbin), Form("h1D_typeA_%s_Mbin%d;Tracklet #eta;Entries", name_in.c_str(), Mbin), nEtaBin, EtaEdge_min, EtaEdge_max)
0152             )
0153         ).second;
0154     }
0155 
0156     map_in.insert( std::make_pair(
0157             Form("h1D_Inclusive100"),
0158             new TH1D(Form("h1D_%s_Inclusive100", name_in.c_str()), Form("h1D_%s_Inclusive100;Tracklet #eta;Entries", name_in.c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0159         )
0160     ).second;
0161     
0162     map_in.insert( std::make_pair(
0163             Form("h1D_Inclusive70"),
0164             new TH1D(Form("h1D_%s_Inclusive70", name_in.c_str()), Form("h1D_%s_Inclusive70;Tracklet #eta;Entries", name_in.c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0165         )
0166     ).second;
0167 
0168     map_in.insert( std::make_pair(
0169             Form("h1D_typeA_Inclusive100"),
0170             new TH1D(Form("h1D_typeA_%s_Inclusive100", name_in.c_str()), Form("h1D_typeA_%s_Inclusive100;Tracklet #eta;Entries", name_in.c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0171         )
0172     ).second;
0173     
0174     map_in.insert( std::make_pair(
0175             Form("h1D_typeA_Inclusive70"),
0176             new TH1D(Form("h1D_typeA_%s_Inclusive70", name_in.c_str()), Form("h1D_typeA_%s_Inclusive70;Tracklet #eta;Entries", name_in.c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0177         )
0178     ).second;
0179 
0180 }
0181 
0182 void PreparedNdEta::PrepareHistFits()
0183 {
0184     std::vector<int> insert_check; insert_check.clear();
0185     bool isInserted = false;
0186 
0187     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0188     // note : multiplicity, dNdeta, type A, inclusive[centrality segemntation, 70%, 100%]
0189     PrepareTrackletEtaHist(h1D_FitBkg_RecoTrackletEta_map, "FitBkg_RecoTrackletEta");
0190     // PrepareTrackletEtaHist(h1D_FitBkg_RecoTrackletEtaPerEvt_map);
0191     // PrepareTrackletEtaHist(h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map);
0192     PrepareTrackletEtaHist(h1D_RotatedBkg_RecoTrackletEta_map, "RotatedBkg_RecoTrackletEta");
0193     // PrepareTrackletEtaHist(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map);
0194     // PrepareTrackletEtaHist(h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map);
0195     
0196 
0197 
0198     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0199     hstack1D_map.clear();
0200     if (runnumber == -1){
0201         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0202             isInserted = hstack1D_map.insert( std::make_pair(
0203                 Form("hstack1D_TrueEta_Mbin%d", Mbin),
0204                 new THStack(Form("hstack1D_TrueEta_Mbin%d", Mbin), Form("hstack1D_TrueEta_Mbin%d;PHG4Particle #eta;Entries", Mbin))
0205             )
0206             ).second; insert_check.push_back(isInserted);
0207         }
0208 
0209         isInserted = hstack1D_map.insert( std::make_pair(
0210                 Form("hstack1D_TrueEta_Inclusive100"),
0211                 new THStack(Form("hstack1D_TrueEta_Inclusive100"), Form("hstack1D_TrueEta_Inclusive100;PHG4Particle #eta;Entries"))
0212         )
0213         ).second; insert_check.push_back(isInserted);
0214 
0215         isInserted = hstack1D_map.insert( std::make_pair(
0216                 Form("hstack1D_TrueEta_Inclusive70"),
0217                 new THStack(Form("hstack1D_TrueEta_Inclusive70"), Form("hstack1D_TrueEta_Inclusive70;PHG4Particle #eta;Entries"))
0218         )
0219         ).second; insert_check.push_back(isInserted);
0220     }
0221 
0222     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0223     // note : for the delta phi, {typeA, rotated} x {eta} x {Mbin [0-100%], Inclusive70, Inclusive100}
0224     for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++){
0225         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0226             isInserted = hstack1D_map.insert( std::make_pair(
0227                     Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", Mbin, eta_bin),
0228                     new THStack(Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", Mbin, eta_bin), Form("hstack1D_DeltaPhi_Mbin%d_Eta%d;Pair #Delta#phi [radian];Entries", Mbin, eta_bin))
0229                 )
0230             ).second; insert_check.push_back(isInserted);
0231 
0232             
0233             isInserted = hstack1D_map.insert( std::make_pair(
0234                     Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", Mbin, eta_bin),
0235                     new THStack(Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", Mbin, eta_bin), Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated;Pair #Delta#phi [radian];Entries", Mbin, eta_bin))
0236                 )
0237             ).second; insert_check.push_back(isInserted);
0238             
0239 
0240 
0241             isInserted = hstack1D_map.insert( std::make_pair(
0242                     Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", Mbin, eta_bin),
0243                     new THStack(Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", Mbin, eta_bin), Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d;Pair #Delta#phi [radian] (type A);Entries", Mbin, eta_bin))
0244                 )
0245             ).second; insert_check.push_back(isInserted);
0246 
0247             isInserted = hstack1D_map.insert( std::make_pair(
0248                     Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", Mbin, eta_bin),
0249                     new THStack(Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", Mbin, eta_bin), Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated;Pair #Delta#phi [radian] (type A);Entries", Mbin, eta_bin))
0250                 )
0251             ).second; insert_check.push_back(isInserted);
0252         }
0253 
0254         // note : inclusive 100
0255         isInserted = hstack1D_map.insert( std::make_pair(
0256                 Form("hstack1D_DeltaPhi_Eta%d_Inclusive100", eta_bin),
0257                 new THStack(Form("hstack1D_DeltaPhi_Eta%d_Inclusive100", eta_bin), Form("hstack1D_DeltaPhi_Eta%d_Inclusive100;Pair #Delta#phi [radian];Entries", eta_bin))
0258             )
0259         ).second; insert_check.push_back(isInserted);
0260 
0261         
0262         isInserted = hstack1D_map.insert( std::make_pair(
0263                 Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated", eta_bin),
0264                 new THStack(Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated", eta_bin), Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated;Pair #Delta#phi [radian];Entries", eta_bin))
0265             )
0266         ).second; insert_check.push_back(isInserted);
0267         
0268 
0269         isInserted = hstack1D_map.insert( std::make_pair(
0270                 Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100", eta_bin),
0271                 new THStack(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100", eta_bin), Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100;Pair #Delta#phi [radian] (type A);Entries", eta_bin))
0272             )
0273         ).second; insert_check.push_back(isInserted);
0274 
0275         isInserted = hstack1D_map.insert( std::make_pair(
0276                 Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated", eta_bin),
0277                 new THStack(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated", eta_bin), Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated;Pair #Delta#phi [radian] (type A);Entries", eta_bin))
0278             )
0279         ).second; insert_check.push_back(isInserted);
0280 
0281         // note : inclusive 70
0282         isInserted = hstack1D_map.insert( std::make_pair(
0283                 Form("hstack1D_DeltaPhi_Eta%d_Inclusive70", eta_bin),
0284                 new THStack(Form("hstack1D_DeltaPhi_Eta%d_Inclusive70", eta_bin), Form("hstack1D_DeltaPhi_Eta%d_Inclusive70;Pair #Delta#phi [radian];Entries", eta_bin))
0285             )
0286         ).second; insert_check.push_back(isInserted);
0287 
0288         
0289         isInserted = hstack1D_map.insert( std::make_pair(
0290                 Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated", eta_bin),
0291                 new THStack(Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated", eta_bin), Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated;Pair #Delta#phi [radian];Entries", eta_bin))
0292             )
0293         ).second; insert_check.push_back(isInserted);
0294         
0295 
0296         isInserted = hstack1D_map.insert( std::make_pair(
0297                 Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70", eta_bin),
0298                 new THStack(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70", eta_bin), Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70;Pair #Delta#phi [radian] (type A);Entries", eta_bin))
0299             )
0300         ).second; insert_check.push_back(isInserted);
0301 
0302         isInserted = hstack1D_map.insert( std::make_pair(
0303                 Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated", eta_bin),
0304                 new THStack(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated", eta_bin), Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated;Pair #Delta#phi [radian] (type A);Entries", eta_bin))
0305             )
0306         ).second; insert_check.push_back(isInserted);
0307     } // note : end of eta_bin loop
0308 
0309     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0310     // note : for the fittings, {typeA, rotated} x {eta} x {Mbin [0-100%], Inclusive70, Inclusive100}
0311     f1_BkgPol2_Fit_map.clear();
0312     f1_BkgPol2_Draw_map.clear();
0313 
0314     f1_SigBkgPol2_Fit_map.clear();
0315     f1_SigBkgPol2_DrawSig_map.clear();
0316     f1_SigBkgPol2_DrawBkgPol2_map.clear();
0317 
0318     for (auto &pair : hstack1D_map){
0319         if (pair.first.find("hstack1D_DeltaPhi") != std::string::npos || pair.first.find("hstack1D_typeA_DeltaPhi") != std::string::npos){
0320 
0321             std::string f1_name = pair.first + "_BkgPol2_Fit";
0322             isInserted = f1_BkgPol2_Fit_map.insert(
0323                 std::make_pair(
0324                     f1_name,
0325                     new TF1(f1_name.c_str(), bkg_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 6)
0326                 )
0327             ).second; insert_check.push_back(isInserted);
0328             f1_BkgPol2_Fit_map[f1_name] -> SetNpx(1000);
0329 
0330 
0331             f1_name = pair.first + "_BkgPol2_Draw";
0332             isInserted = f1_BkgPol2_Draw_map.insert(
0333                 std::make_pair(
0334                     f1_name,
0335                     new TF1(f1_name.c_str(), full_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 4)
0336                 )
0337             ).second; insert_check.push_back(isInserted);
0338             f1_BkgPol2_Draw_map[f1_name] -> SetLineColor(6);
0339             f1_BkgPol2_Draw_map[f1_name] -> SetLineStyle(9);
0340             f1_BkgPol2_Draw_map[f1_name] -> SetLineWidth(3);
0341             f1_BkgPol2_Draw_map[f1_name] -> SetNpx(1000);
0342 
0343             
0344             f1_name = pair.first + "_SigBkgPol2_Fit";
0345             isInserted = f1_SigBkgPol2_Fit_map.insert(
0346                 std::make_pair(
0347                     f1_name,
0348                     new TF1(f1_name.c_str(), gaus_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 7)
0349                 )
0350             ).second; insert_check.push_back(isInserted);
0351             f1_SigBkgPol2_Fit_map[f1_name] -> SetLineColor(46);
0352             f1_SigBkgPol2_Fit_map[f1_name] -> SetLineStyle(1);
0353             f1_SigBkgPol2_Fit_map[f1_name] -> SetNpx(1000);
0354 
0355 
0356             f1_name = pair.first + "_SigBkgPol2_DrawSig";
0357             isInserted = f1_SigBkgPol2_DrawSig_map.insert(
0358                 std::make_pair(
0359                     f1_name,
0360                     new TF1(f1_name.c_str(), gaus_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 4)
0361                 )
0362             ).second; insert_check.push_back(isInserted);
0363             f1_SigBkgPol2_DrawSig_map[f1_name] -> SetLineColor(46);
0364             f1_SigBkgPol2_DrawSig_map[f1_name] -> SetLineStyle(2);
0365             f1_SigBkgPol2_DrawSig_map[f1_name] -> SetNpx(1000);
0366 
0367 
0368             f1_name = pair.first + "_SigBkgPol2_DrawBkgPol2";
0369             isInserted = f1_SigBkgPol2_DrawBkgPol2_map.insert(
0370                 std::make_pair(
0371                     f1_name,
0372                     new TF1(f1_name.c_str(), full_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 4)
0373                 )
0374             ).second; insert_check.push_back(isInserted);
0375             f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetLineColor(8);
0376             f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetLineStyle(2);
0377             f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetNpx(1000);
0378 
0379         }
0380     }
0381 
0382     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0383     // note : signal (eta-vtxZ) [centrality segemntation, 70% and 100%]
0384     h2D_map.clear();
0385 
0386 
0387     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0388     // note : truth tracks (eta-vtxZ) [70% and 100%], {coarse bin, fine bin}
0389 
0390     hstack2D_map.clear();
0391     if (runnumber == -1){
0392 
0393         // note : centrality inclusive 
0394         isInserted = hstack2D_map.insert( std::make_pair(
0395                 Form("hstack2D_TrueEtaVtxZ_Inclusive100"),
0396                 new THStack(Form("hstack2D_TrueEtaVtxZ_Inclusive100"), Form("hstack2D_TrueEtaVtxZ_Inclusive100;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0397             )
0398         ).second; insert_check.push_back(isInserted);
0399 
0400         isInserted = hstack2D_map.insert( std::make_pair(
0401                 Form("hstack2D_TrueEtaVtxZ_Inclusive100_FineBin"),
0402                 new THStack(Form("hstack2D_TrueEtaVtxZ_Inclusive100_FineBin"), Form("hstack2D_TrueEtaVtxZ_Inclusive100_FineBin;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0403             )
0404         ).second; insert_check.push_back(isInserted);
0405 
0406         // note : centrality inclusive 
0407         isInserted = hstack2D_map.insert( std::make_pair(
0408                 Form("hstack2D_TrueEtaVtxZ_Inclusive70"),
0409                 new THStack(Form("hstack2D_TrueEtaVtxZ_Inclusive70"), Form("hstack2D_TrueEtaVtxZ_Inclusive70;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0410             )
0411         ).second; insert_check.push_back(isInserted);
0412 
0413         isInserted = hstack2D_map.insert( std::make_pair(
0414                 Form("hstack2D_TrueEtaVtxZ_Inclusive70_FineBin"),
0415                 new THStack(Form("hstack2D_TrueEtaVtxZ_Inclusive70_FineBin"), Form("hstack2D_TrueEtaVtxZ_Inclusive70_FineBin;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0416             )
0417         ).second; insert_check.push_back(isInserted);   
0418     }
0419     
0420     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0421     // note : good pair (eta-vtxZ) [70% and 100%], {typeA, rotated, normal} {coarse bin, fine bin}
0422 
0423     // note : inclusive 100
0424     isInserted = hstack2D_map.insert( std::make_pair(
0425             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100"),
0426             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100;Pair #eta;INTT vtxZ [cm]")) 
0427         )
0428     ).second; insert_check.push_back(isInserted);
0429 
0430     isInserted = hstack2D_map.insert( std::make_pair(
0431             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100"),
0432             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100;Pair #eta;INTT vtxZ [cm]")) 
0433         )
0434     ).second; insert_check.push_back(isInserted);
0435 
0436     // note : normal fine bin
0437     isInserted = hstack2D_map.insert( std::make_pair(
0438             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"),
0439             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin;Pair #eta;INTT vtxZ [cm]")) 
0440         )
0441     ).second; insert_check.push_back(isInserted);
0442 
0443     isInserted = hstack2D_map.insert( std::make_pair(
0444             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"),
0445             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin;Pair #eta;INTT vtxZ [cm]")) 
0446         )
0447     ).second; insert_check.push_back(isInserted);
0448     
0449 
0450     // note : rotated
0451     isInserted = hstack2D_map.insert( std::make_pair(
0452             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"),
0453             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated;Pair #eta;INTT vtxZ [cm]")) 
0454         )
0455     ).second; insert_check.push_back(isInserted);
0456 
0457     isInserted = hstack2D_map.insert( std::make_pair(
0458             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"),
0459             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated;Pair #eta;INTT vtxZ [cm]")) 
0460         )
0461     ).second; insert_check.push_back(isInserted);
0462 
0463     // note : inclusive 70
0464     isInserted = hstack2D_map.insert( std::make_pair(
0465             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70"),
0466             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70;Pair #eta;INTT vtxZ [cm]")) 
0467         )
0468     ).second; insert_check.push_back(isInserted);
0469 
0470     isInserted = hstack2D_map.insert( std::make_pair(
0471             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70"),
0472             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70;Pair #eta;INTT vtxZ [cm]")) 
0473         )
0474     ).second; insert_check.push_back(isInserted);
0475 
0476     // note : normal fine bin
0477     isInserted = hstack2D_map.insert( std::make_pair(
0478             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"),
0479             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin;Pair #eta;INTT vtxZ [cm]")) 
0480         )
0481     ).second; insert_check.push_back(isInserted);
0482 
0483     isInserted = hstack2D_map.insert( std::make_pair(
0484             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"),
0485             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin;Pair #eta;INTT vtxZ [cm]")) 
0486         )
0487     ).second; insert_check.push_back(isInserted);
0488     
0489     // note : rotated
0490     isInserted = hstack2D_map.insert( std::make_pair(
0491             Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"),
0492             new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated;Pair #eta;INTT vtxZ [cm]")) 
0493         )
0494     ).second; insert_check.push_back(isInserted);
0495 
0496     isInserted = hstack2D_map.insert( std::make_pair(
0497             Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"),
0498             new THStack(Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"), Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated;Pair #eta;INTT vtxZ [cm]")) 
0499         )
0500     ).second; insert_check.push_back(isInserted);
0501 
0502 
0503     for (int isInserted_tag : insert_check){
0504         if (isInserted_tag == 0){
0505             std::cout<<"Histogram insertion failed"<<std::endl;
0506             exit(1);
0507         }
0508     }
0509 
0510 }
0511 
0512 std::map<int, int> PreparedNdEta::GetVtxZIndexMap()
0513 {
0514     if(h2D_input_map.find("h2D_RecoEvtCount_vtxZCentrality") == h2D_input_map.end()){
0515         std::cout<<"Error: h2D_RecoEvtCount_vtxZCentrality not found"<<std::endl;
0516         exit(1);
0517     }
0518 
0519     std::map<int, int> vtxZIndexMap; vtxZIndexMap.clear();
0520 
0521     for (int y_i = 1; y_i <= h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetNbinsY(); y_i++){
0522 
0523         if (h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinLowEdge(y_i) < cut_INTTvtxZ.first || h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinLowEdge(y_i) > cut_INTTvtxZ.second){
0524             continue;
0525         }
0526         if (h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinUpEdge(y_i) < cut_INTTvtxZ.first || h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinUpEdge(y_i) > cut_INTTvtxZ.second){
0527             continue;
0528         }
0529 
0530         vtxZIndexMap.insert(
0531             std::make_pair(
0532                 y_i - 1,
0533                 y_i - 1
0534             )
0535         );
0536     }
0537 
0538     std::cout<<"The selected INTT vtxZ bin : [";
0539     for (auto &pair : vtxZIndexMap){
0540         std::cout<<pair.first<<", ";
0541     }
0542     std::cout<<"]"<<std::endl;
0543 
0544     for (auto &pair : vtxZIndexMap){
0545         std::cout<<"vtxZ index : "<<pair.first<<", range : {"<<h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinLowEdge(pair.first + 1)<<", "<<h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetXaxis()->GetBinUpEdge(pair.first + 1)<<"}"<<std::endl;
0546     }
0547 
0548     return vtxZIndexMap;
0549 }
0550 
0551 // note : return the bin index of eta, vtxZ, Mbin, typeA, rotated, finebin
0552 std::tuple<int, int, int, int, int, int> PreparedNdEta::GetHistStringInfo(std::string hist_name)
0553 {
0554 
0555     std::cout<<"In GetHistStringInfo(), input name: "<<hist_name<<std::endl;
0556     
0557     int eta_bin;
0558     int vtxz_bin;
0559     int Mbin;
0560 
0561     std::string eta_bin_str = "";
0562     std::string vtxz_bin_str = "";
0563     std::string Mbin_str = "";
0564 
0565     // note : eta bin
0566     if (hist_name.find("_Eta") == std::string::npos) {eta_bin = -1;}
0567     else {
0568         for (int i = 0; i < hist_name.size(); i++){
0569             int start_index = hist_name.find("_Eta") + 4 + i;
0570             if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0571                 eta_bin_str += hist_name[start_index];
0572             }
0573             else {
0574                 break;
0575             }
0576 
0577             if (start_index == hist_name.size() - 1){
0578                 break;
0579             }
0580         }
0581 
0582         eta_bin = (eta_bin_str.length() == 0) ? -1 : std::stoi(eta_bin_str);
0583     }
0584 
0585     
0586 
0587     // note : vtxZ bin
0588     if (hist_name.find("_VtxZ") == std::string::npos) {vtxz_bin = -1;}
0589     else {
0590         for (int i = 0; i < hist_name.size(); i++){
0591             int start_index = hist_name.find("_VtxZ") + 5 + i;
0592             if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0593                 vtxz_bin_str += hist_name[start_index];
0594             }
0595             else {
0596                 break;
0597             }
0598 
0599             if (start_index == hist_name.size() - 1){
0600                 break;
0601             }
0602         }
0603 
0604         vtxz_bin = (vtxz_bin_str.length() == 0) ? -1 : std::stoi(vtxz_bin_str);
0605     }
0606 
0607     
0608 
0609     // note : Mbin
0610     if (hist_name.find("_Mbin") == std::string::npos) {Mbin = -1;}
0611     else {
0612         for (int i = 0; i < hist_name.size(); i++){
0613             int start_index = hist_name.find("_Mbin") + 5 + i;
0614             if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0615                 Mbin_str += hist_name[start_index];
0616             }
0617             else {
0618                 break;
0619             }
0620 
0621             if (start_index == hist_name.size() - 1){
0622                 break;
0623             }
0624         }
0625 
0626         Mbin = (Mbin_str.length() == 0) ? -1 : std::stoi(Mbin_str);
0627     }
0628 
0629     // note : typeA, rotated, finebin
0630     int typeA = (hist_name.find("_typeA") != std::string::npos) ? 1 : 0;
0631     int rotated = (hist_name.find("_rotated") != std::string::npos) ? 1 : 0;
0632     int finebin = (hist_name.find("_FineBin") != std::string::npos) ? 1 : 0;
0633 
0634 
0635     return std::make_tuple(eta_bin, vtxz_bin, Mbin, typeA, rotated, finebin);
0636 }
0637 
0638 void PreparedNdEta::PrepareStacks()
0639 {   
0640 
0641     std::cout<<"In PrepareStacks()"<<std::endl;
0642 
0643     int h1D_eta_bin, h1D_vtxz_bin, h1D_Mbin;
0644     int h2D_eta_bin, h2D_vtxz_bin, h2D_Mbin;
0645     int typeA, rotated, finebin;
0646 
0647     for (auto &pair : h1D_input_map)
0648     {
0649         std::tie(h1D_eta_bin, h1D_vtxz_bin, h1D_Mbin, typeA, rotated, finebin) = GetHistStringInfo(pair.first);
0650 
0651         // note : the vtxz_bin of the hist is not in the interest
0652         if (vtxZ_index_map.find(h1D_vtxz_bin) == vtxZ_index_map.end()){
0653             continue;
0654         }
0655 
0656         // note : reco, pair DeltaPhi for each Mbin and eta_bin, stacked up by vertex z bin
0657         // note : with typeA, rotate, inclusive100, inclusive70
0658         if (pair.first.find("h1D") != std::string::npos && pair.first.find("_DeltaPhi") != std::string::npos && h1D_Mbin != -1 && h1D_eta_bin != -1){
0659             std::cout<<"----- "<<pair.first<<", eta_bin: "<<h1D_eta_bin<<", vtxz_bin: "<<h1D_vtxz_bin<<", Mbin: "<<h1D_Mbin<<", typeA: "<<typeA<<", rotated: "<<rotated<<", finebin: "<<finebin<<std::endl;
0660 
0661 
0662             if (typeA == 0 && rotated == 0){
0663                 pair.second->SetFillColor(ROOT_color_code[hstack1D_map[Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0664                 if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)) == hstack1D_map.end() ) {
0665                     std::cout<<Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0666                     exit(1);
0667                 } 
0668                 hstack1D_map[Form("hstack1D_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)]->Add(pair.second);
0669 
0670                 if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)) == hstack1D_map.end() ) {
0671                     std::cout<<Form("hstack1D_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0672                     exit(1);
0673                 } 
0674                 hstack1D_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)]->Add(pair.second);
0675                 
0676                 if (h1D_Mbin <= Semi_inclusive_Mbin){
0677                     if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)) == hstack1D_map.end() ) {
0678                         std::cout<<Form("hstack1D_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0679                         exit(1);
0680                     } 
0681                     hstack1D_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)]->Add(pair.second);
0682                 }
0683             }
0684 
0685             if (typeA == 0 && rotated == 1){
0686                 pair.second->SetFillColor(ROOT_color_code[hstack1D_map[Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0687                 if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)) == hstack1D_map.end() ) {
0688                     std::cout<<Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0689                     exit(1);
0690                 } 
0691                 hstack1D_map[Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)]->Add(pair.second);
0692 
0693                 if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)) == hstack1D_map.end() ) {
0694                     std::cout<<Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0695                     exit(1);
0696                 } 
0697                 hstack1D_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)]->Add(pair.second);
0698 
0699                 if (h1D_Mbin <= Semi_inclusive_Mbin){
0700                     if ( hstack1D_map.find(Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)) == hstack1D_map.end() ) {
0701                         std::cout<<Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0702                         exit(1);
0703                     } 
0704                     hstack1D_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)]->Add(pair.second);
0705                 }
0706             }
0707 
0708             // note : type A
0709             if (typeA == 1 && rotated == 0){
0710                 pair.second->SetFillColor(ROOT_color_code[hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0711                 if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)) == hstack1D_map.end() ) {
0712                     std::cout<<Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0713                     exit(1);
0714                 } 
0715                 hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d", h1D_Mbin, h1D_eta_bin)]->Add(pair.second);
0716 
0717                 if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)) == hstack1D_map.end() ) {
0718                     std::cout<<Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0719                     exit(1);
0720                 } 
0721                 hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100", h1D_eta_bin)]->Add(pair.second);
0722 
0723                 if (h1D_Mbin <= Semi_inclusive_Mbin){
0724                     if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)) == hstack1D_map.end() ) {
0725                         std::cout<<Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0726                         exit(1);
0727                     } 
0728                     hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70", h1D_eta_bin)]->Add(pair.second);
0729                 }
0730             }
0731 
0732             if (typeA == 1 && rotated == 1){
0733                 pair.second->SetFillColor(ROOT_color_code[hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0734                 if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)) == hstack1D_map.end() ) {
0735                     std::cout<<Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0736                     exit(1);
0737                 } 
0738                 hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_rotated", h1D_Mbin, h1D_eta_bin)]->Add(pair.second);
0739 
0740                 if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)) == hstack1D_map.end() ) {
0741                     std::cout<<Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0742                     exit(1);
0743                 } 
0744                 hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_rotated", h1D_eta_bin)]->Add(pair.second);
0745 
0746                 if (h1D_Mbin <= Semi_inclusive_Mbin){
0747                     if ( hstack1D_map.find(Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)) == hstack1D_map.end() ) {
0748                         std::cout<<Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)<<" not found in hstack1D_map !!"<<std::endl;
0749                         exit(1);
0750                     } 
0751                     hstack1D_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_rotated", h1D_eta_bin)]->Add(pair.second);
0752                 }
0753             }
0754 
0755 
0756 
0757         }
0758 
0759 
0760         // note : truth eta distribution, stacked up by vertex Z bin
0761         // note : different Mbin, inclusive100, inclusive70
0762         if (runnumber == -1 && pair.first.find("h1D") != std::string::npos && pair.first.find("_TrueEta") != std::string::npos && h1D_Mbin != -1){
0763             
0764             pair.second->SetFillColor(ROOT_color_code[hstack1D_map[Form("hstack1D_TrueEta_Mbin%d", h1D_Mbin)]->GetNhists() % ROOT_color_code.size()]);
0765 
0766             if (hstack1D_map.find(Form("hstack1D_TrueEta_Mbin%d", h1D_Mbin)) == hstack1D_map.end()){
0767                 std::cout<<Form("hstack1D_TrueEta_Mbin%d", h1D_Mbin)<<" not found in hstack1D_map !!"<<std::endl;
0768                 exit(1);
0769             }
0770             hstack1D_map[Form("hstack1D_TrueEta_Mbin%d", h1D_Mbin)]->Add(pair.second);
0771             
0772             if (hstack1D_map.find(Form("hstack1D_TrueEta_Inclusive100")) == hstack1D_map.end()){
0773                 std::cout<<Form("hstack1D_TrueEta_Inclusive100")<<" not found in hstack1D_map !!"<<std::endl;
0774                 exit(1);
0775             }
0776             hstack1D_map["hstack1D_TrueEta_Inclusive100"] -> Add(pair.second);
0777 
0778             if (h1D_Mbin <= Semi_inclusive_Mbin){
0779                 if (hstack1D_map.find(Form("hstack1D_TrueEta_Inclusive70")) == hstack1D_map.end()){
0780                     std::cout<<Form("hstack1D_TrueEta_Inclusive70")<<" not found in hstack1D_map !!"<<std::endl;
0781                     exit(1);
0782                 }
0783                 hstack1D_map["hstack1D_TrueEta_Inclusive70"] -> Add(pair.second);
0784             }
0785 
0786         }
0787     }
0788 
0789     // for (auto &pair : h2D_input_map){
0790     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100             V
0791     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated     V
0792     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70              V
0793     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated      V
0794     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100         V
0795     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated V
0796     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70          V
0797     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated  V
0798 
0799     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin        V
0800     //     // hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin         V
0801     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin  V
0802     //     // hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin   V
0803 
0804     //     std::tie(h2D_eta_bin, h2D_vtxz_bin, h2D_Mbin, typeA, rotated, finebin) = GetHistStringInfo(pair.first);
0805     //     std::string HS2D = "hstack2D";
0806     //     std::string GPTEVz = "GoodProtoTracklet_EtaVtxZ";
0807     //     if (pair.first.find("h2D_") != std::string::npos && pair.first.find(GPTEVz) != std::string::npos && h2D_Mbin != -1){
0808     //         if (typeA == 0 && rotated == 0 && finebin == 0){
0809     //             hstack2D_map[Form("%s_%s_Inclusive100", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0810 
0811     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0812     //                 hstack2D_map[Form("%s_%s_Inclusive70", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0813     //             }
0814     //         } 
0815 
0816     //         if (typeA == 0 && rotated == 1 && finebin == 0){
0817     //             hstack2D_map[Form("%s_%s_Inclusive100_rotated", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0818 
0819     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0820     //                 hstack2D_map[Form("%s_%s_Inclusive70_rotated", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0821     //             }
0822     //         }
0823 
0824     //         if (typeA == 1 && rotated == 0 && finebin == 0){
0825     //             hstack2D_map[Form("%s_typeA_%s_Inclusive100", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0826 
0827     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0828     //                 hstack2D_map[Form("%s_typeA_%s_Inclusive70", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0829     //             }
0830     //         }
0831 
0832     //         if (typeA == 1 && rotated == 1 && finebin == 0){
0833     //             hstack2D_map[Form("%s_typeA_%s_Inclusive100_rotated", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0834 
0835     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0836     //                 hstack2D_map[Form("%s_typeA_%s_Inclusive70_rotated", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0837     //             }
0838     //         }
0839 
0840     //         if (typeA == 0 && rotated == 0 && finebin == 1){
0841     //             hstack2D_map[Form("%s_%s_Inclusive100_FineBin", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0842 
0843     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0844     //                 hstack2D_map[Form("%s_%s_Inclusive70_FineBin", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0845     //             }
0846     //         }
0847 
0848     //         if (typeA == 1 && rotated == 0 && finebin == 1){
0849     //             hstack2D_map[Form("%s_typeA_%s_Inclusive100_FineBin", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0850 
0851     //             if (h2D_Mbin <= Semi_inclusive_Mbin){
0852     //                 hstack2D_map[Form("%s_typeA_%s_Inclusive70_FineBin", HS2D.c_str(), GPTEVz.c_str())]->Add(pair.second);
0853     //             }
0854     //         }
0855     //     }
0856     // }
0857 
0858     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0859         hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0860         hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0861         hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0862 
0863         hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0864         hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0865         hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0866 
0867         if (Mbin <= Semi_inclusive_Mbin){
0868             hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0869             hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0870             hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0871 
0872             hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0873             hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0874             hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0875         }
0876     }
0877 
0878     // note : for MC
0879     // hstack2D_TrueEtaVtxZ_Inclusive100         V
0880     // hstack2D_TrueEtaVtxZ_Inclusive100_FineBin V
0881     // hstack2D_TrueEtaVtxZ_Inclusive70          V
0882     // hstack2D_TrueEtaVtxZ_Inclusive70_FineBin  V
0883     // hstack1D_TrueEta_Mbin%d", Mbin),          V
0884     // hstack1D_TrueEta_Inclusive100"),          V
0885     // hstack1D_TrueEta_Inclusive70"),           V
0886     if (runnumber == -1){
0887         for (auto &pair : h2D_input_map)
0888         {
0889             std::tie(h2D_eta_bin, h2D_vtxz_bin, h2D_Mbin, typeA, rotated, finebin) = GetHistStringInfo(pair.first);
0890 
0891             if (pair.first.find("h2D_TrueEtaVtxZ") != std::string::npos && h2D_Mbin != -1){
0892                 
0893                 std::cout<<" MCMCMC ----- "<<pair.first<<", eta_bin: "<<h2D_eta_bin<<", vtxz_bin: "<<h2D_vtxz_bin<<", Mbin: "<<h2D_Mbin<<", typeA: "<<typeA<<", rotated: "<<rotated<<", finebin: "<<finebin<<std::endl;
0894 
0895                 if(finebin == 1){
0896                     hstack2D_map[Form("hstack2D_TrueEtaVtxZ_Inclusive100_FineBin")] -> Add(pair.second);
0897 
0898                     if (h2D_Mbin <= Semi_inclusive_Mbin){
0899                         hstack2D_map[Form("hstack2D_TrueEtaVtxZ_Inclusive70_FineBin")] -> Add(pair.second);
0900                     }
0901                 }
0902                 else {
0903                     hstack2D_map[Form("hstack2D_TrueEtaVtxZ_Inclusive100")] -> Add(pair.second);
0904 
0905                     if (h2D_Mbin <= Semi_inclusive_Mbin){
0906                         hstack2D_map[Form("hstack2D_TrueEtaVtxZ_Inclusive70")] -> Add(pair.second);
0907                     }
0908                 }
0909             }
0910             
0911             
0912         }
0913     }
0914         
0915 
0916 
0917 
0918     std::cout<<"End the PrepareStacks()"<<std::endl;
0919 }
0920 
0921 void PreparedNdEta::DoFittings()
0922 {
0923     std::cout<<"In DoFittings()"<<std::endl;
0924 
0925     int stack_count = 0;
0926 
0927     for (auto &pair : hstack1D_map){
0928         if (pair.first.find("hstack1D") != std::string::npos && pair.first.find("_DeltaPhi") != std::string::npos){
0929             
0930             if (stack_count % 20 == 0){
0931                 std::cout<<"Fitting stack : "<<stack_count<<", "<<pair.first<<std::endl;
0932             }
0933 
0934             auto temp_hist = (TH1D*) pair.second -> GetStack() -> Last();
0935             std::vector<double> N_group_info = find_Ngroup(temp_hist);
0936 
0937             std::string f1_BkgPol2_Fit_map_key  = pair.first + "_BkgPol2_Fit";
0938             std::string f1_BkgPol2_Draw_map_key = pair.first + "_BkgPol2_Draw";
0939 
0940             if (f1_BkgPol2_Fit_map.find(f1_BkgPol2_Fit_map_key) == f1_BkgPol2_Fit_map.end()){
0941                 std::cout<<f1_BkgPol2_Fit_map_key<<" not found in f1_BkgPol2_Fit_map !!"<<std::endl;
0942                 exit(1);
0943             }
0944 
0945             std::string f1_SigBkgPol2_Fit_map_key         = pair.first + "_SigBkgPol2_Fit";
0946             std::string f1_SigBkgPol2_DrawSig_map_key     = pair.first + "_SigBkgPol2_DrawSig";
0947             std::string f1_SigBkgPol2_DrawBkgPol2_map_key = pair.first + "_SigBkgPol2_DrawBkgPol2";
0948 
0949             if (f1_SigBkgPol2_Fit_map.find(f1_SigBkgPol2_Fit_map_key) == f1_SigBkgPol2_Fit_map.end()){
0950                 std::cout<<f1_SigBkgPol2_Fit_map_key<<" not found in f1_SigBkgPol2_Fit_map !!"<<std::endl;
0951                 exit(1);
0952             }
0953 
0954             double hist_offset = get_dist_offset(temp_hist, 15); // todo:
0955             double signal_region_l = cut_GoodProtoTracklet_DeltaPhi->first;
0956             double signal_region_r = cut_GoodProtoTracklet_DeltaPhi->second; 
0957 
0958             // todo:
0959             // note : par[0] + par[1]* (x[0]-par[3]) + par[2] * pow((x[0]-par[3]),2);
0960             // std::cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<std::endl;
0961             f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> SetParameters(hist_offset, 0, -0.2, 0, signal_region_l, signal_region_r);
0962             f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> FixParameter(4, signal_region_l);
0963             f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> FixParameter(5, signal_region_r);
0964             f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> SetParLimits(2, -100, 0);
0965 
0966             // std::cout<<"for "<<f1_BkgPol2_Fit_map_key<<"fit parameters: "<<hist_offset<<", 0, -0.2, 0, "<<signal_region_l<<", "<<signal_region_r<<std::endl;
0967 
0968             temp_hist -> Fit(f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key], "NQ");
0969 
0970             // std::cout<<"for "<<f1_BkgPol2_Fit_map_key<<"fit parameters: "<<f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(0)<<", "<<f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(1)<<", "<<f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(2)<<", "<<f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(3)<<std::endl;
0971 
0972             // std::cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<std::endl;
0973 
0974             f1_BkgPol2_Draw_map[f1_BkgPol2_Draw_map_key] -> SetParameters(
0975                 f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(0),
0976                 f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(1),
0977                 f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(2),
0978                 f1_BkgPol2_Fit_map[f1_BkgPol2_Fit_map_key] -> GetParameter(3)
0979             );
0980 
0981             // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0982 
0983 
0984             // note : the normal cases, with the peak
0985             if (pair.first.find("_rotated") == std::string::npos)
0986             {
0987                 // note : for the signal+bkg fit
0988                 double gaus_height = temp_hist->GetBinContent(temp_hist->GetMaximumBin()) - hist_offset;
0989                 double gaus_width = fabs(N_group_info[3]-N_group_info[2]) / 2.;
0990 
0991                 // note : 
0992                 // note : par[0] : size
0993                 // note : par[1] : mean
0994                 // note : par[2] : width
0995                 // note : double gaus_func = par[0] * TMath::Gaus(x[0],par[1],par[2]);
0996                 // note : double pol2_func = par[3] + par[4]* (x[0]-par[6]) + par[5] * pow((x[0]-par[6]),2);
0997 
0998                 // std::cout<<"=================================================================================================================================="<<std::endl;
0999 
1000                 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> SetParameters(
1001                     gaus_height, 0, gaus_width, 
1002                     hist_offset, 0, -0.2, 0
1003                 );
1004                 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> SetParLimits(5, -100, 0);
1005 
1006                 // std::cout<<"for "<<f1_SigBkgPol2_Fit_map_key<<", all pars: "<<gaus_height<<", 0, "<<gaus_width<<", "<<hist_offset<<", 0, -0.2, 0"<<std::endl;
1007 
1008                 temp_hist -> Fit(f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key], "NQ");
1009 
1010                 // std::cout<<"for "<<f1_SigBkgPol2_Fit_map_key<<", all pars: "<<
1011                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(0)<<", "<<
1012                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(1)<<", "<<
1013                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(2)<<", "<<
1014                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(3)<<", "<<
1015                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(4)<<", "<<
1016                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(5)<<", "<<
1017                 // f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(6)<<std::endl;
1018                 
1019                 // std::cout<<"=================================================================================================================================="<<std::endl;
1020 
1021                 f1_SigBkgPol2_DrawSig_map[f1_SigBkgPol2_DrawSig_map_key] -> SetParameters(
1022                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(0),
1023                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(1),
1024                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(2)
1025                 );
1026 
1027                 f1_SigBkgPol2_DrawBkgPol2_map[f1_SigBkgPol2_DrawBkgPol2_map_key] -> SetParameters(
1028                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(3),
1029                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(4),
1030                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(5),
1031                     f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(6)
1032                 );
1033             } // note : end of the normal cases, with the peak
1034             
1035             // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1036 
1037             // note : save the plot into root files
1038             c1 -> cd();
1039             hstack1D_map[pair.first] -> Draw(""); // note : inclusive
1040 
1041             if (pair.first.find("_rotated") == std::string::npos) // note : normal
1042             {
1043                 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> Draw("l same");
1044                 f1_SigBkgPol2_DrawSig_map[f1_SigBkgPol2_DrawSig_map_key] -> Draw("l same");
1045                 f1_SigBkgPol2_DrawBkgPol2_map[f1_SigBkgPol2_DrawBkgPol2_map_key] -> Draw("l same");
1046             }
1047 
1048             f1_BkgPol2_Draw_map[f1_BkgPol2_Draw_map_key] -> Draw("l same");
1049 
1050             file_out -> cd();
1051             c1 -> Write(Form("c1_%s", pair.first.c_str()));
1052             c1 -> Clear();
1053 
1054             if (pair.first.find("_rotated") == std::string::npos){ // note : normal
1055                 temp_hist -> SetFillColor(0);
1056                 temp_hist -> SetLineColor(9);
1057                 temp_hist -> SetMaximum(temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.6);
1058                 
1059                 auto temp_hist_rotate = (TH1D*) hstack1D_map[pair.first + "_rotated"] -> GetStack() -> Last();
1060                 temp_hist_rotate -> SetFillColor(0);
1061                 temp_hist_rotate -> SetLineColor(46);
1062 
1063                 c1 -> cd();
1064                 temp_hist -> Draw();
1065                 temp_hist_rotate -> Draw("same");
1066 
1067                 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> Draw("l same");
1068                 f1_SigBkgPol2_DrawSig_map[f1_SigBkgPol2_DrawSig_map_key] -> Draw("l same");
1069                 f1_SigBkgPol2_DrawBkgPol2_map[f1_SigBkgPol2_DrawBkgPol2_map_key] -> Draw("l same");
1070 
1071                 f1_BkgPol2_Draw_map[f1_BkgPol2_Draw_map_key] -> Draw("l same");
1072 
1073                 file_out -> cd();
1074                 c1 -> Write(Form("c1_%s_h1DFit", pair.first.c_str()));
1075                 c1 -> Clear();
1076             }
1077 
1078         }
1079 
1080         stack_count++;
1081     }
1082 
1083     std::cout<<"End the DoFittings()"<<std::endl;
1084 }
1085 
1086 double PreparedNdEta::get_h2D_GoodProtoTracklet_count(TH2D * hist_in, int eta_bin_in)
1087 {
1088     double count = 0;
1089     
1090     // note : hist
1091     // note : X: eta
1092     // note : Y: VtxZ
1093 
1094     for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){ // note : vtxZ bin
1095             
1096         if (vtxZ_index_map.find(yi - 1) == vtxZ_index_map.end()){
1097             continue;
1098         }
1099 
1100         count += hist_in -> GetBinContent(eta_bin_in + 1, yi);
1101     }
1102 
1103     return count;
1104 
1105 }
1106 
1107 double PreparedNdEta::get_EvtCount(TH2D * hist_in, int centrality_bin_in)
1108 {
1109     double count = 0;
1110 
1111     // note : X : vtxZ
1112     // note : Y : centrality
1113 
1114     if (centrality_bin_in == 100){
1115         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
1116 
1117             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1118                 continue;
1119             }
1120 
1121             for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){ // note : centrality bin
1122                 count += hist_in -> GetBinContent(xi, yi);
1123             }
1124         }
1125     }
1126     else if (centrality_bin_in == 70){
1127         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
1128 
1129             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1130                 continue;
1131             }
1132 
1133             for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){ // note : centrality bin
1134                 if (yi-1 <= Semi_inclusive_Mbin){
1135                     count += hist_in -> GetBinContent(xi, yi);
1136                 }
1137             }
1138         }
1139     }
1140     else {
1141         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
1142 
1143             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1144                 continue;
1145             }
1146 
1147             count += hist_in -> GetBinContent(xi, centrality_bin_in + 1);
1148         }
1149     }
1150 
1151     return count;
1152 
1153 }
1154 
1155 void PreparedNdEta::Convert_to_PerEvt(TH1D * hist_in, double Nevent)
1156 {
1157     for (int i = 1; i <= hist_in->GetNbinsX(); i++){
1158         hist_in -> SetBinContent(i, hist_in -> GetBinContent(i) / Nevent);
1159         hist_in -> SetBinError(i, hist_in -> GetBinError(i) / Nevent);
1160     }
1161 
1162 }
1163 
1164 void PreparedNdEta::PrepareMultiplicity()
1165 {
1166     // note : data {typeA, inclusive} {fit bkg, rotated bkg} {centrality bin, inclusive 100, inclusive 70}
1167 
1168     // note : "h2D_RecoEvtCount_vtxZCentrality", the event counting 
1169     
1170     // note : h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d
1171     // note : h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated
1172     // note :   V hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100
1173     // note :   V hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated
1174     // note :   V hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70
1175     // note :   V hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated
1176 
1177     // note : h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d
1178     // note : h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated
1179     // note :   V hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100
1180     // note :   V hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated
1181     // note :   V hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70
1182     // note :   V hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated
1183 
1184     // note : V h1D_FitBkg_RecoTrackletEta_map
1185     // note : h1D_FitBkg_RecoTrackletEtaPerEvt_map
1186     // note : h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map
1187     // note : V h1D_RotatedBkg_RecoTrackletEta_map
1188     // note : h1D_RotatedBkg_RecoTrackletEtaPerEvt_map
1189     // note : h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map
1190     // note :   V h1D_Mbin%d
1191     // note :   V h1D_typeA_Mbin%d
1192     // note :   V h1D_Inclusive100")
1193     // note :   V h1D_Inclusive70")
1194     // note :   V h1D_typeA_Inclusive100")
1195     // note :   V h1D_typeA_Inclusive70")
1196 
1197     // note : f1_BkgPol2_Draw_map
1198 
1199     
1200 
1201     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
1202         for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++){
1203             
1204             // note : normal
1205             double pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_DeltaPhi_Mbin%d_Eta%d_BkgPol2_Draw", Mbin, eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1206             h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Mbin%d",Mbin)] -> SetBinContent(
1207                 eta_bin + 1,
1208                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)], eta_bin) - pol2_bkg_integral
1209             );
1210             h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Mbin%d",Mbin)] -> SetBinContent(
1211                 eta_bin + 1,
1212                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)], eta_bin) - 
1213                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)], eta_bin)
1214             );
1215 
1216             // note : type A
1217             pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_typeA_DeltaPhi_Mbin%d_Eta%d_BkgPol2_Draw", Mbin, eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1218             h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinContent(
1219                 eta_bin + 1,
1220                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)], eta_bin) - pol2_bkg_integral
1221             );
1222             h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinContent(
1223                 eta_bin + 1,
1224                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)], eta_bin) - 
1225                 get_h2D_GoodProtoTracklet_count(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)], eta_bin)
1226             );
1227         }   
1228     }
1229 
1230     // note : inclusive 100, inclusive 70
1231     for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++){
1232         
1233         // note : normal
1234         double pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive100_BkgPol2_Draw", eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1235         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Inclusive100")] -> SetBinContent(
1236             eta_bin + 1,
1237             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100")]->GetStack()->Last()), eta_bin) - pol2_bkg_integral
1238         );
1239         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Inclusive100")] -> SetBinContent(
1240             eta_bin + 1,
1241             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100")]->GetStack()->Last()), eta_bin) - 
1242             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated")]->GetStack()->Last()), eta_bin)
1243         );
1244 
1245         pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_DeltaPhi_Eta%d_Inclusive70_BkgPol2_Draw", eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1246         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Inclusive70")] -> SetBinContent(
1247             eta_bin + 1,
1248             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70")]->GetStack()->Last()), eta_bin) - pol2_bkg_integral
1249         );
1250         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Inclusive70")] -> SetBinContent(
1251             eta_bin + 1,
1252             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70")]->GetStack()->Last()), eta_bin) - 
1253             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated")]->GetStack()->Last()), eta_bin)
1254         );
1255 
1256 
1257 
1258         // note : type A
1259         pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive100_BkgPol2_Draw", eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1260         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive100")] -> SetBinContent(
1261             eta_bin + 1,
1262             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100")]->GetStack()->Last()), eta_bin) - pol2_bkg_integral
1263         );
1264         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive100")] -> SetBinContent(
1265             eta_bin + 1,
1266             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100")]->GetStack()->Last()), eta_bin) - 
1267             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated")]->GetStack()->Last()), eta_bin)
1268         );
1269 
1270         pol2_bkg_integral = fabs(f1_BkgPol2_Draw_map[Form("hstack1D_typeA_DeltaPhi_Eta%d_Inclusive70_BkgPol2_Draw", eta_bin)] -> Integral( cut_GoodProtoTracklet_DeltaPhi->first, cut_GoodProtoTracklet_DeltaPhi->second )) / ((EtaEdge_max - EtaEdge_min)/double(nEtaBin));
1271         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive70")] -> SetBinContent(
1272             eta_bin + 1,
1273             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70")]->GetStack()->Last()), eta_bin) - pol2_bkg_integral
1274         );
1275         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive70")] -> SetBinContent(
1276             eta_bin + 1,
1277             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70")]->GetStack()->Last()), eta_bin) - 
1278             get_h2D_GoodProtoTracklet_count(((TH2D*)hstack2D_map[Form("hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated")]->GetStack()->Last()), eta_bin)
1279         );
1280 
1281     }
1282 
1283     
1284 
1285     // note : MC {centrality bin, inclusive 100, inclusive 70}
1286     // note : "h2D_TrueEvtCount_vtxZCentrality", the event counting
1287 
1288     // note : hstack1D_TrueEta_Mbin%d", Mbin)
1289     // note : V hstack1D_TrueEta_Inclusive100")
1290     // note : V hstack1D_TrueEta_Inclusive70")
1291     // note : V hstack2D_TrueEtaVtxZ_Inclusive100
1292     // note : V hstack2D_TrueEtaVtxZ_Inclusive100_FineBin
1293     // note : V hstack2D_TrueEtaVtxZ_Inclusive70
1294     // note : V hstack2D_TrueEtaVtxZ_Inclusive70_FineBin
1295 
1296     // note : h1D_TruedNdEta_map
1297 
1298 }
1299 
1300 void PreparedNdEta::PreparedNdEtaHist()
1301 {
1302     std::cout<<1111111<<std::endl;
1303     if (runnumber == -1){
1304         h1D_TruedNdEta_map.clear();
1305 
1306         h1D_TruedNdEta_map.insert( 
1307             std::make_pair(
1308                 "h1D_TruedNdEta_Inclusive100",
1309                 (TH1D*)hstack1D_map["hstack1D_TrueEta_Inclusive100"] -> GetStack() -> Last()
1310             )
1311         );
1312         
1313         std::cout<<22222222<<std::endl;
1314 
1315         Convert_to_PerEvt(h1D_TruedNdEta_map["h1D_TruedNdEta_Inclusive100"], get_EvtCount(h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"], 100));
1316         
1317         std::cout<<33333333<<std::endl;
1318 
1319         h1D_TruedNdEta_map.insert( 
1320             std::make_pair(
1321                 "h1D_TruedNdEta_Inclusive70",
1322                 (TH1D*)hstack1D_map["hstack1D_TrueEta_Inclusive70"] -> GetStack() -> Last()
1323             )
1324         );
1325         Convert_to_PerEvt(h1D_TruedNdEta_map["h1D_TruedNdEta_Inclusive70"], get_EvtCount(h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"], 70));
1326 
1327         std::cout<<44444<<std::endl;
1328 
1329         for (int i = 0; i < nCentrality_bin; i++){
1330             h1D_TruedNdEta_map.insert( 
1331                 std::make_pair(
1332                     Form("h1D_TruedNdEta_Mbin%d", i),
1333                     (TH1D*)hstack1D_map[Form("hstack1D_TrueEta_Mbin%d", i)] -> GetStack() -> Last()
1334                 )
1335             );
1336             Convert_to_PerEvt(h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", i)], get_EvtCount(h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"], i));
1337         }
1338     } // note : end of truth
1339 
1340     std::cout<<5555<<std::endl;
1341 
1342     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1343     // note : FitBkg, {TypeA, normal} x {Mbin, Inclusive100, Inclusive70}
1344     h1D_FitBkg_RecoTrackletEtaPerEvt_map.clear();
1345     h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1346         std::make_pair(            
1347             "h1D_Inclusive100",
1348             (TH1D*)h1D_FitBkg_RecoTrackletEta_map["h1D_Inclusive100"] -> Clone("h1D_FitBkg_RecoTrackletEtaPerEvt_Inclusive100")
1349         )
1350     );
1351 
1352     std::cout<<"h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1353 
1354     Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 100));
1355     std::cout<<"h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1356 
1357     std::cout<<6666<<std::endl;
1358 
1359     h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1360         std::make_pair(            
1361             "h1D_Inclusive70",
1362             (TH1D*)h1D_FitBkg_RecoTrackletEta_map["h1D_Inclusive70"] -> Clone("h1D_FitBkg_RecoTrackletEtaPerEvt_Inclusive70")
1363         )
1364     );
1365     Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 70));
1366 
1367     h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1368         std::make_pair(            
1369             "h1D_typeA_Inclusive100",
1370             (TH1D*)h1D_FitBkg_RecoTrackletEta_map["h1D_typeA_Inclusive100"] -> Clone("h1D_typeA_FitBkg_RecoTrackletEtaPerEvt_Inclusive100")
1371         )
1372     );
1373     Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 100));
1374 
1375     h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1376         std::make_pair(            
1377             "h1D_typeA_Inclusive70",
1378             (TH1D*)h1D_FitBkg_RecoTrackletEta_map["h1D_typeA_Inclusive70"] -> Clone("h1D_typeA_FitBkg_RecoTrackletEtaPerEvt_Inclusive70")
1379         )
1380     );
1381     Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 70));
1382 
1383     for (int i = 0; i < nCentrality_bin; i++){
1384         h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1385             std::make_pair(                
1386                 Form("h1D_Mbin%d", i),
1387                 (TH1D*)h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Mbin%d", i)] -> Clone(Form("h1D_FitBkg_RecoTrackletEtaPerEvt_Mbin%d", i))
1388             )
1389         );
1390         Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], i));
1391 
1392         h1D_FitBkg_RecoTrackletEtaPerEvt_map.insert( 
1393             std::make_pair(                
1394                 Form("h1D_typeA_Mbin%d", i),
1395                 (TH1D*)h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d", i)] -> Clone(Form("h1D_typeA_FitBkg_RecoTrackletEtaPerEvt_Mbin%d", i))
1396             )
1397         );
1398         Convert_to_PerEvt(h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], i));
1399     }
1400 
1401     std::cout<<7777<<std::endl;
1402 
1403     // note : copy to the next map
1404     // note : FitBkg, {TypeA, normal} x {Mbin, Inclusive100, Inclusive70}
1405     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1406 
1407     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1408         std::make_pair(            
1409             "h1D_Inclusive100",
1410             (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"] -> Clone("h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_Inclusive100")
1411         )
1412     );
1413 
1414     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1415         std::make_pair(            
1416             "h1D_Inclusive70",
1417             (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"] -> Clone("h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_Inclusive70")
1418         )
1419     );
1420 
1421     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1422         std::make_pair(            
1423             "h1D_typeA_Inclusive100",
1424             (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"] -> Clone("h1D_typeA_FitBkg_RecoTrackletEtaPerEvtPostAC_Inclusive100")
1425         )
1426     );
1427 
1428     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1429         std::make_pair(            
1430             "h1D_typeA_Inclusive70",
1431             (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"] -> Clone("h1D_typeA_FitBkg_RecoTrackletEtaPerEvtPostAC_Inclusive70")
1432         )
1433     );
1434 
1435     std::cout<<88888<<std::endl;
1436 
1437     for (int i = 0; i < nCentrality_bin; i++){
1438         h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1439             std::make_pair(                
1440                 Form("h1D_Mbin%d", i),
1441                 (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)] -> Clone(Form("h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_Mbin%d", i))
1442             )
1443         );
1444 
1445         h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1446             std::make_pair(                
1447                 Form("h1D_typeA_Mbin%d", i),
1448                 (TH1D*)h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)] -> Clone(Form("h1D_typeA_FitBkg_RecoTrackletEtaPerEvtPostAC_Mbin%d", i))
1449             )
1450         );
1451     }
1452 
1453     std::cout<<9999<<std::endl;
1454 
1455     // note : RotatedBkg, {TypeA, normal} x {Mbin, Inclusive100, Inclusive70}
1456     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1457     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.clear();
1458     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1459         std::make_pair(            
1460             "h1D_Inclusive100",
1461             (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map["h1D_Inclusive100"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_Inclusive100")
1462         )
1463     );
1464     Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 100));
1465 
1466     std::cout<<10101010<<std::endl;
1467 
1468     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1469         std::make_pair(            
1470             "h1D_Inclusive70",
1471             (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map["h1D_Inclusive70"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt_Inclusive70")
1472         )
1473     );
1474     Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 70));
1475 
1476     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1477         std::make_pair(            
1478             "h1D_typeA_Inclusive100",
1479             (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map["h1D_typeA_Inclusive100"] -> Clone("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvt_Inclusive100")
1480         )
1481     );
1482     Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 100));
1483 
1484     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1485         std::make_pair(            
1486             "h1D_typeA_Inclusive70",
1487             (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map["h1D_typeA_Inclusive70"] -> Clone("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvt_Inclusive70")
1488         )
1489     );
1490     Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], 70));
1491 
1492     for (int i = 0; i < nCentrality_bin; i++){
1493         h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1494             std::make_pair(                
1495                 Form("h1D_Mbin%d", i),
1496                 (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Mbin%d", i)] -> Clone(Form("h1D_RotatedBkg_RecoTrackletEtaPerEvt_Mbin%d", i))
1497             )
1498         );
1499         Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], i));
1500 
1501         h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert( 
1502             std::make_pair(                
1503                 Form("h1D_typeA_Mbin%d", i),
1504                 (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d", i)] -> Clone(Form("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvt_Mbin%d", i))
1505             )
1506         );
1507         Convert_to_PerEvt(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)], get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], i));
1508     }
1509 
1510     std::cout<<10101111<<std::endl;
1511 
1512     // note : copy to the next map
1513     // note : RotatedBkg, {TypeA, normal} x {Mbin, Inclusive100, Inclusive70}
1514     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1515 
1516     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1517         std::make_pair(            
1518             "h1D_Inclusive100",
1519             (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Inclusive100")
1520         )
1521     );
1522 
1523 
1524     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1525         std::make_pair(            
1526             "h1D_Inclusive70",
1527             (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Inclusive70")
1528         )
1529     );
1530 
1531     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1532         std::make_pair(            
1533             "h1D_typeA_Inclusive100",
1534             (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"] -> Clone("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Inclusive100")
1535         )
1536     );
1537 
1538     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1539         std::make_pair(            
1540             "h1D_typeA_Inclusive70",
1541             (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"] -> Clone("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Inclusive70")
1542         )
1543     );
1544 
1545     for (int i = 0; i < nCentrality_bin; i++){
1546         h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1547             std::make_pair(                
1548                 Form("h1D_Mbin%d", i),
1549                 (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)] -> Clone(Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Mbin%d", i))
1550             )
1551         );
1552 
1553         h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert( 
1554             std::make_pair(                
1555                 Form("h1D_typeA_Mbin%d", i),
1556                 (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)] -> Clone(Form("h1D_typeA_RotatedBkg_RecoTrackletEtaPerEvtPostAC_Mbin%d", i))
1557             )
1558         );
1559     }
1560 
1561     std::cout<<101022222<<std::endl;
1562 
1563     std::cout<<"a h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1564     // note : after the copy, make sure that the title is also changed
1565     for (auto pair : h1D_FitBkg_RecoTrackletEta_map){
1566         h1D_FitBkg_RecoTrackletEtaPerEvt_map[pair.first] -> SetTitle(h1D_FitBkg_RecoTrackletEtaPerEvt_map[pair.first]->GetName());
1567         h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[pair.first] -> SetTitle(h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[pair.first]->GetName());
1568 
1569         h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[pair.first] -> SetTitle(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[pair.first]->GetName());
1570         h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[pair.first] -> SetTitle(h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[pair.first]->GetName());
1571     }
1572     std::cout<<"b h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1573     
1574 
1575     // note : apply the alpha correction from the "h2D_alpha_correction_map_in"
1576     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1577 
1578     if (h2D_alpha_correction_map_in.size() != 0 && ApplyAlphaCorr == true){
1579         // h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map
1580         // h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map
1581 
1582         // h2D_FitBkg_alpha_correction
1583         // h2D_RotatedBkg_alpha_correction
1584         // h2D_typeA_FitBkg_alpha_correction
1585         // h2D_typeA_RotatedBkg_alpha_correction
1586 
1587 
1588         // note : for [0-100%]
1589         // note : X: eta, Y: Mbin
1590         // todo: need to fix the zero division issue
1591         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
1592             for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1593             {   
1594                 // note : FitBkg
1595                 h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1596                 h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1597                 h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1598                 h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1599 
1600                 // note : RotatedBkg
1601                 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1602                 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1603 
1604                 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1605                 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, Mbin + 1));
1606             }
1607         }
1608 
1609         // note : for Inclusive 70
1610         // todo : the index in the Y axis (for Mbin) is given by (nCentrality_bin + 1)
1611         for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1612         {   
1613             // note : FitBkg
1614             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1615             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1616             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1617             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1618 
1619             // note : RotatedBkg
1620             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1621             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1622 
1623             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1624             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 1)));
1625         }
1626 
1627         // note : for Inclusive 100
1628         // todo : the index in the Y axis (for Mbin) is given by (nCentrality_bin + 2)
1629         for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1630         {   
1631             // note : FitBkg
1632             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1633             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1634             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> SetBinContent(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1635             h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> SetBinError(eta_bin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_FitBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1636 
1637             // note : RotatedBkg
1638             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1639             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1640 
1641             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> SetBinContent(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinContent(eta_bin + 1) / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1642             h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> SetBinError(eta_bin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinError(eta_bin + 1)     / h2D_alpha_correction_map_in["h2D_typeA_RotatedBkg_alpha_correction"] -> GetBinContent(eta_bin + 1, (nCentrality_bin + 2)));
1643         }
1644 
1645     }
1646 
1647     std::cout<<10103333<<std::endl;
1648 
1649 }
1650 
1651 void PreparedNdEta::DeriveAlphaCorrection()
1652 {
1653     h2D_alpha_correction_map_out.insert(
1654         std::make_pair(
1655             "h2D_FitBkg_alpha_correction",
1656             new TH2D("h2D_FitBkg_alpha_correction", "h2D_FitBkg_alpha_correction", nEtaBin, 0, nEtaBin, nCentrality_bin+2, 0, nCentrality_bin+2) // todo: additional 2 bins for inclusive 70, 100
1657         )
1658     );   
1659 
1660     h2D_alpha_correction_map_out.insert(
1661         std::make_pair(
1662             "h2D_RotatedBkg_alpha_correction",
1663             new TH2D("h2D_RotatedBkg_alpha_correction", "h2D_RotatedBkg_alpha_correction", nEtaBin, 0, nEtaBin, nCentrality_bin+2, 0, nCentrality_bin+2) // todo: additional 2 bins for inclusive 70, 100
1664         )
1665     );
1666 
1667     h2D_alpha_correction_map_out.insert(
1668         std::make_pair(
1669             "h2D_typeA_FitBkg_alpha_correction",
1670             new TH2D("h2D_typeA_FitBkg_alpha_correction", "h2D_typeA_FitBkg_alpha_correction", nEtaBin, 0, nEtaBin, nCentrality_bin+2, 0, nCentrality_bin+2) // todo: additional 2 bins for inclusive 70, 100
1671         )
1672     );   
1673 
1674     h2D_alpha_correction_map_out.insert(
1675         std::make_pair(
1676             "h2D_typeA_RotatedBkg_alpha_correction",
1677             new TH2D("h2D_typeA_RotatedBkg_alpha_correction", "h2D_typeA_RotatedBkg_alpha_correction", nEtaBin, 0, nEtaBin, nCentrality_bin+2, 0, nCentrality_bin+2) // todo: additional 2 bins for inclusive 70, 100
1678         )
1679     );   
1680 
1681 
1682     // note : derive the alpha correction regardless whether you have the input alpha correction map or not
1683     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1684     // h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map
1685     // h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map
1686 
1687     // h2D_FitBkg_alpha_correction
1688     // h2D_RotatedBkg_alpha_correction
1689     // h2D_typeA_FitBkg_alpha_correction
1690     // h2D_typeA_RotatedBkg_alpha_correction
1691 
1692     // note : for [0-100%]
1693     // note : X: eta, Y: Mbin
1694     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
1695         for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1696         {   
1697             // note : FitBkg
1698             // todo : how to properly handle the errors
1699             h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, Mbin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1700             // h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, Mbin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1701             
1702             h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, Mbin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1703             // h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, Mbin + 1, h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1704 
1705             // note : RotatedBkg
1706             h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, Mbin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1707             // h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, Mbin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1708 
1709             h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, Mbin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1710             // h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, Mbin + 1, h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d",Mbin)] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", Mbin)] -> GetBinContent(eta_bin + 1));
1711         }
1712     }
1713 
1714     // note : for Inclusive 70
1715     // todo : the index in the Y axis (for Mbin) is given by (nCentrality_bin + 1)
1716     for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1717     {   
1718         // note : FitBkg
1719         // todo : how to properly handle the errors
1720         h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 1), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1721         // h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 1), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1722         
1723         h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 1), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1724         // h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 1), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1725 
1726         // note : RotatedBkg
1727         h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 1), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1728         // h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 1), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive70")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1729 
1730         h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 1), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1731         // h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 1), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive70")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive70" )] -> GetBinContent(eta_bin + 1));
1732     }
1733 
1734     // note : for Inclusive 100
1735     // todo : the index in the Y axis (for Mbin) is given by (nCentrality_bin + 2)
1736     for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1737     {   
1738         // note : FitBkg
1739         // todo : how to properly handle the errors
1740         h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 2), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1741         // h2D_alpha_correction_map_out["h2D_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 2), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1742         
1743         h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 2), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1744         // h2D_alpha_correction_map_out["h2D_typeA_FitBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 2), h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1745 
1746         // note : RotatedBkg
1747         h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 2), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1748         // h2D_alpha_correction_map_out["h2D_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 2), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Inclusive100")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1749 
1750         h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinContent(eta_bin + 1, (nCentrality_bin + 2), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinContent(eta_bin + 1) / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1751         // h2D_alpha_correction_map_out["h2D_typeA_RotatedBkg_alpha_correction"] -> SetBinError(eta_bin + 1, (nCentrality_bin + 2), h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Inclusive100")] -> GetBinError(eta_bin + 1)     / h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Inclusive100" )] -> GetBinContent(eta_bin + 1));
1752     }
1753 
1754     std::cout<<"cccc h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1755 
1756 }
1757 
1758 void PreparedNdEta::EndRun()
1759 {
1760     
1761 
1762     std::cout<<"aa"<<111<<std::endl;
1763 
1764     std::cout<<"c h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1765 
1766     file_out_dNdEta -> cd();
1767     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
1768         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin)] -> Write();
1769         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin)] -> Write();
1770         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin)] -> Write();
1771         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin)] -> Write();
1772 
1773         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin)] -> Write();
1774         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin)] -> Write();
1775         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin)] -> Write();
1776         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin)] -> Write();
1777 
1778         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin)] -> Write();
1779         h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin)] -> Write();
1780         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin)] -> Write();
1781         h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin)] -> Write();
1782     }
1783 
1784     std::cout<<"d h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1785 
1786     std::cout<<"aa"<<222<<std::endl;
1787 
1788     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100"] -> Write();
1789     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"] -> Write();
1790     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70"] -> Write();
1791     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"] -> Write();
1792 
1793     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100"] -> Write();
1794     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_rotated"] -> Write();
1795     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70"] -> Write();
1796     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_rotated"] -> Write();
1797 
1798     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"] -> Write();
1799     hstack2D_map["hstack2D_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"] -> Write();
1800     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive100_FineBin"] -> Write();
1801     hstack2D_map["hstack2D_typeA_GoodProtoTracklet_EtaVtxZ_Inclusive70_FineBin"] -> Write();
1802 
1803     std::cout<<"aa"<<333<<std::endl;
1804 
1805     std::cout<<"e h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1806 
1807     // note : reco, multiplicity
1808     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
1809         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Mbin%d",Mbin)] -> Write();
1810         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Mbin%d",Mbin)] -> Write();
1811 
1812         h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d",Mbin)] -> Write();
1813         h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Mbin%d",Mbin)] -> Write();
1814     }
1815 
1816     h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Inclusive100")] -> Write();
1817     h1D_FitBkg_RecoTrackletEta_map[Form("h1D_Inclusive70")]  -> Write();
1818     h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive100")] -> Write();
1819     h1D_FitBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive70")]  -> Write();
1820 
1821     h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Inclusive100")] -> Write();
1822     h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_Inclusive70")]  -> Write();
1823     h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive100")] -> Write();
1824     h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_typeA_Inclusive70")]  -> Write();
1825     
1826     std::cout<<"aa"<<444<<std::endl;
1827 
1828     std::cout<<"f h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1829 
1830     // note : multiplicity
1831     if (runnumber == -1){
1832 
1833         hstack2D_map["hstack2D_TrueEtaVtxZ_Inclusive100"] -> Write();
1834         hstack2D_map["hstack2D_TrueEtaVtxZ_Inclusive100_FineBin"] -> Write();
1835         hstack2D_map["hstack2D_TrueEtaVtxZ_Inclusive70"] -> Write();
1836         hstack2D_map["hstack2D_TrueEtaVtxZ_Inclusive70_FineBin"] -> Write();
1837 
1838         hstack1D_map["hstack1D_TrueEta_Inclusive100"] -> Write();
1839         hstack1D_map["hstack1D_TrueEta_Inclusive70"] -> Write();
1840 
1841         for (int i = 0; i < nCentrality_bin; i++){
1842             hstack1D_map[Form("hstack1D_TrueEta_Mbin%d", i)] -> Write();
1843         }
1844 
1845         std::cout<<"aa"<<555<<std::endl;
1846     }
1847 
1848     // note : N things for per event 
1849     if (runnumber == -1){
1850         h1D_TruedNdEta_map["h1D_TruedNdEta_Inclusive100"] -> Write("h1D_TruedNdEta_Inclusive100");
1851         h1D_TruedNdEta_map["h1D_TruedNdEta_Inclusive70"] -> Write("h1D_TruedNdEta_Inclusive70");
1852         for (int i = 0; i < nCentrality_bin; i++){
1853             h1D_TruedNdEta_map[Form("h1D_TruedNdEta_Mbin%d", i)] -> Write(Form("h1D_TruedNdEta_Mbin%d", i));
1854         }
1855     }
1856 
1857     std::cout<<"g h1D_FitBkg_RecoTrackletEtaPerEvt_map[\"h1D_Inclusive100\"]->GetEntries(): "<<h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"]->GetEntries()<<std::endl;
1858 
1859     std::cout<<"aa"<<666<<std::endl;
1860 
1861     for (auto pair : h1D_FitBkg_RecoTrackletEtaPerEvt_map){
1862         std::cout<<"maps in h1D_FitBkg_RecoTrackletEtaPerEvt_map: "<<pair.first<<", hist entries: "<<pair.second->GetEntries()<<std::endl;
1863     }
1864 
1865 
1866     h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"] -> Write("a");
1867     h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"] -> Write("b");
1868     h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"] -> Write("c");
1869     h1D_FitBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"] -> Write("d");
1870 
1871     std::cout<<"aa"<<66666<<std::endl;
1872 
1873     for (int i = 0; i < nCentrality_bin; i++){
1874         h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)] -> Write();
1875         h1D_FitBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)] -> Write();
1876     }
1877 
1878     std::cout<<"aa"<<777<<std::endl;
1879 
1880     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive100"] -> Write();
1881     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_Inclusive70"] -> Write();
1882     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive100"] -> Write();
1883     h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_typeA_Inclusive70"] -> Write();
1884 
1885     for (int i = 0; i < nCentrality_bin; i++){
1886         h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_Mbin%d", i)] -> Write();
1887         h1D_RotatedBkg_RecoTrackletEtaPerEvt_map[Form("h1D_typeA_Mbin%d", i)] -> Write();
1888     }
1889 
1890     std::cout<<"aa"<<888<<std::endl;
1891 
1892     // note : post the alpha correction {if it's applied}
1893     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_Inclusive100"] -> Write();
1894     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_Inclusive70"] -> Write();
1895     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_typeA_Inclusive100"] -> Write();
1896     h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_typeA_Inclusive70"] -> Write();
1897 
1898     for (int i = 0; i < nCentrality_bin; i++){
1899         h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d", i)] -> Write();
1900         h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d", i)] -> Write();
1901     }
1902 
1903     std::cout<<"aa"<<999<<std::endl;
1904 
1905     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_Inclusive100"] -> Write();
1906     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_Inclusive70"] -> Write();
1907     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_typeA_Inclusive100"] -> Write();
1908     h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_typeA_Inclusive70"] -> Write();
1909 
1910     for (int i = 0; i < nCentrality_bin; i++){
1911         h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_Mbin%d", i)] -> Write();
1912         h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_typeA_Mbin%d", i)] -> Write();
1913     }
1914 
1915     std::cout<<"aa"<<111101010<<std::endl;
1916 
1917     // note : the alpha correction map
1918     for (auto &pair : h2D_alpha_correction_map_out){
1919         pair.second -> Write();
1920     }
1921 
1922     std::cout<<"aa"<<111122222<<std::endl;
1923 
1924     file_out_dNdEta -> Close();
1925 
1926     std::cout<<"aa"<<11113333<<std::endl;
1927 
1928     file_out -> cd();
1929     file_out -> Close();
1930 }
1931 
1932 std::vector<double> PreparedNdEta::find_Ngroup(TH1D * hist_in)
1933 {
1934     double Highest_bin_Content  = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
1935     double Highest_bin_Center   = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
1936 
1937     int group_Nbin = 0;
1938     int peak_group_ID = -9999;
1939     double group_entry = 0;
1940     double peak_group_ratio;
1941     std::vector<int> group_Nbin_vec; group_Nbin_vec.clear();
1942     std::vector<double> group_entry_vec; group_entry_vec.clear();
1943     std::vector<double> group_widthL_vec; group_widthL_vec.clear();
1944     std::vector<double> group_widthR_vec; group_widthR_vec.clear();
1945 
1946     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
1947         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
1948         double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
1949 
1950         if (bin_content != 0){
1951             
1952             if (group_Nbin == 0) {
1953                 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1954             }
1955 
1956             group_Nbin += 1; 
1957             group_entry += bin_content;
1958         }
1959         else if (bin_content == 0 && group_Nbin != 0){
1960             group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1961             group_Nbin_vec.push_back(group_Nbin);
1962             group_entry_vec.push_back(group_entry);
1963             group_Nbin = 0;
1964             group_entry = 0;
1965         }
1966     }
1967     if (group_Nbin != 0) {
1968         group_Nbin_vec.push_back(group_Nbin);
1969         group_entry_vec.push_back(group_entry);
1970         group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
1971     } // note : the last group at the edge
1972 
1973     // note : find the peak group
1974     for (int i = 0; i < int(group_Nbin_vec.size()); i++){
1975         if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
1976             peak_group_ID = i;
1977             break;
1978         }
1979     }
1980     
1981     // note : On Nov 6, 2024, we no longer need to calculate the ratio of the peak group
1982     // double denominator = (accumulate( group_entry_vec.begin(), group_entry_vec.end(), 0.0 ));
1983     // denominator = (denominator <= 0) ? 1. : denominator;
1984     // peak_group_ratio = group_entry_vec[peak_group_ID] / denominator;
1985     peak_group_ratio = -999;
1986 
1987     double peak_group_left  = (double(group_Nbin_vec.size()) == 0) ? -999 : group_widthL_vec[peak_group_ID];
1988     double peak_group_right = (double(group_Nbin_vec.size()) == 0) ? 999  : group_widthR_vec[peak_group_ID];
1989 
1990     // for (int i = 0; i < group_Nbin_vec.size(); i++)
1991     // {
1992     //     cout<<" "<<endl;
1993     //     cout<<"group size : "<<group_Nbin_vec[i]<<endl;
1994     //     cout<<"group entry : "<<group_entry_vec[i]<<endl;
1995     //     cout<<group_widthL_vec[i]<<" "<<group_widthR_vec[i]<<endl;
1996     // }
1997 
1998     // cout<<" "<<endl;
1999     // cout<<"N group : "<<group_Nbin_vec.size()<<endl;
2000     // cout<<"Peak group ID : "<<peak_group_ID<<endl;
2001     // cout<<"peak group width : "<<group_widthL_vec[peak_group_ID]<<" "<<group_widthR_vec[peak_group_ID]<<endl;
2002     // cout<<"ratio : "<<peak_group_ratio<<endl;
2003     
2004     // note : {N_group, ratio (if two), peak widthL, peak widthR}
2005     return {double(group_Nbin_vec.size()), peak_group_ratio, peak_group_left, peak_group_right};
2006 }
2007 
2008 double PreparedNdEta::get_dist_offset(TH1D * hist_in, int check_N_bin) // note : check_N_bin 1 to N bins of hist
2009 {
2010     if (check_N_bin < 0 || check_N_bin > hist_in -> GetNbinsX()) {std::cout<<" wrong check_N_bin "<<std::endl; exit(1);}
2011     double total_entry = 0;
2012     for (int i = 0; i < check_N_bin; i++)
2013     {
2014         total_entry += hist_in -> GetBinContent(i+1); // note : 1, 2, 3.....
2015         total_entry += hist_in -> GetBinContent(hist_in -> GetNbinsX() - i);
2016     }
2017 
2018     return total_entry/double(2. * check_N_bin);
2019 }