NewCode2024/PreparedNdEtaClus/PreparedNdEtaClusEach.cpp

0001 #include "PreparedNdEtaClusEach.h"
0002
0003 PreparedNdEtaClusEach::PreparedNdEtaClusEach(
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
0011     bool ApplyAlphaCorr_in,
0012     bool isTypeA_in,
0013     std::pair<double,double> cut_INTTvtxZ_in,
0014     int SelectedMbin_in
0015 ) : process_id(process_id_in),
0016     runnumber(runnumber_in),
0017     input_directory(input_directory_in),
0018     input_file_name(input_file_name_in),
0019     output_directory(output_directory_in),
0020     output_file_name_suffix(output_file_name_suffix_in),
0021
0022     ApplyAlphaCorr(ApplyAlphaCorr_in),
0023     isTypeA(isTypeA_in),
0024     cut_INTTvtxZ(cut_INTTvtxZ_in),
0025     SelectedMbin(SelectedMbin_in)
0026 {
0027     PrepareInputRootFie();
0028
0029     if (SelectedMbin >= nCentrality_bin && SelectedMbin != Semi_inclusive_Mbin*10 && SelectedMbin != 100){
0030         std::cout << "Error : SelectedMbin is out of range" << std::endl;
0031         std::cout << "SelectedMbin : " << SelectedMbin <<", but : nCentrality_bin: "<< nCentrality_bin << std::endl;
0032         exit(1);
0033     }
0034
0035     vtxZ_index_map = GetVtxZIndexMap();
0036
0037     PrepareOutPutFileName();
0038
0039     PrepareOutPutRootFile();
0040
0041     PrepareHist();
0042
0043     c1 = new TCanvas("c1", "c1", 950, 800);
0044
0045     h1D_alpha_correction_map_in.clear();
0046     h1D_alpha_correction_map_out.clear();
0047
0048 }
0049
0050 void PreparedNdEtaClusEach::PrepareInputRootFie()
0051 {
0052     file_in = TFile::Open(Form("%s/%s",input_directory.c_str(),input_file_name.c_str()));
0053     if (file_in == nullptr)
0054     {
0055         std::cout << "Error: file_in can not be opened" << std::endl;
0056         exit(1);
0057     }
0058
0059     for (TObject* keyAsObj : *file_in->GetListOfKeys())
0060     {
0061         auto key = dynamic_cast<TKey*>(keyAsObj);
0062         std::string hist_name  = key->GetName();
0063         std::string class_name = key->GetClassName();
0064
0065         if (class_name == "TH2D")
0066         {
0067             h2D_input_map[hist_name.c_str()] = (TH2D*) file_in -> Get( hist_name.c_str() );
0068         }
0069         else if (class_name == "TH1D")
0070         {
0071             h1D_input_map[hist_name.c_str()] = (TH1D*) file_in -> Get( hist_name.c_str() );
0072             h1D_input_map[hist_name.c_str()] -> SetLineColor(kBlack);
0073             h1D_input_map[hist_name.c_str()] -> SetLineWidth(2);
0074         }
0075     }
0076     std::cout<<"In PrepareInputRootFie(), number of input TH1D: "<< h1D_input_map.size() <<std::endl;
0077     std::cout<<"In PrepareInputRootFie(), number of input TH2D: "<< h2D_input_map.size() <<std::endl;
0078
0079     tree_in = (TTree*) file_in -> Get("tree_par");
0080
0081     centrality_edges = 0;
0082
0083     tree_in -> SetBranchAddress("centrality_edges", &centrality_edges);
0084     tree_in -> SetBranchAddress("nCentrality_bin", &nCentrality_bin);
0085
0086     tree_in -> SetBranchAddress("CentralityFineEdge_min", &CentralityFineEdge_min);
0087     tree_in -> SetBranchAddress("CentralityFineEdge_max", &CentralityFineEdge_max);
0088     tree_in -> SetBranchAddress("nCentralityFineBin", &nCentralityFineBin);
0089
0090     tree_in -> SetBranchAddress("EtaEdge_min", &EtaEdge_min);
0091     tree_in -> SetBranchAddress("EtaEdge_max", &EtaEdge_max);
0092     tree_in -> SetBranchAddress("nEtaBin", &nEtaBin);
0093
0094     tree_in -> SetBranchAddress("VtxZEdge_min", &VtxZEdge_min);
0095     tree_in -> SetBranchAddress("VtxZEdge_max", &VtxZEdge_max);
0096     tree_in -> SetBranchAddress("nVtxZBin", &nVtxZBin);
0097
0098     tree_in -> GetEntry(0);
0099
0100     std::cout<<"------------------------------------------------------------"<<std::endl;
0101     std::cout<<"In PrepareInputRootFile(), nEntries : "<<tree_in->GetEntries()<<std::endl;
0102     std::cout<<"nCentrality_bin : "<<nCentrality_bin<<std::endl;
0103     for(int i = 0; i < nCentrality_bin; i++)
0104     {
0105         std::cout<<"centrality_edges["<<i<<"] : "<<centrality_edges->at(i)<<std::endl;
0106     }
0107     std::cout<<"CentralityFineBin : "<<nCentralityFineBin<<", "<<CentralityFineEdge_min<<", "<<CentralityFineEdge_max<<std::endl;
0108     std::cout<<"EtaBin : "<<nEtaBin<<", "<<EtaEdge_min<<", "<<EtaEdge_max<<std::endl;
0109     std::cout<<"VtxZBin : "<<nVtxZBin<<", "<<VtxZEdge_min<<", "<<VtxZEdge_max<<std::endl;
0110     std::cout<<"------------------------------------------------------------"<<std::endl;
0111 }
0112
0113 std::map<int, int> PreparedNdEtaClusEach::GetVtxZIndexMap()
0114 {
0115     if(h2D_input_map.find("h2D_RecoEvtCount_vtxZCentrality") == h2D_input_map.end()){
0116         std::cout<<"Error: h2D_RecoEvtCount_vtxZCentrality not found"<<std::endl;
0117         exit(1);
0118     }
0119
0120     std::map<int, int> vtxZIndexMap; vtxZIndexMap.clear();
0121
0122     for (int y_i = 1; y_i <= h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetNbinsY(); y_i++){
0123
0124         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){
0125             continue;
0126         }
0127         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){
0128             continue;
0129         }
0130
0131         vtxZIndexMap.insert(
0132             std::make_pair(
0133                 y_i - 1,
0134                 y_i - 1
0135             )
0136         );
0137     }
0138
0139     std::cout<<"The selected INTT vtxZ bin : [";
0140     for (auto &pair : vtxZIndexMap){
0141         std::cout<<pair.first<<", ";
0142     }
0143     std::cout<<"]"<<std::endl;
0144
0145     for (auto &pair : vtxZIndexMap){
0146         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;
0147     }
0148
0149     return vtxZIndexMap;
0150 }
0151
0152 void PreparedNdEtaClusEach::PrepareOutPutFileName()
0153 {
0154     std::string job_index = std::to_string( process_id );
0155     int job_index_len = 5;
0156     job_index.insert(0, job_index_len - job_index.size(), '0');
0157
0158     std::string runnumber_str = std::to_string( runnumber );
0159     if (runnumber != -1){
0160         int runnumber_str_len = 8;
0161         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0162     }
0163
0164     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0165         output_file_name_suffix = "_" + output_file_name_suffix;
0166     }
0167
0168     output_filename = "PreparedNdEtaClus";
0169     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0170
0171     output_filename += (ApplyAlphaCorr) ? "_ApplyAlphaCorr" : "";
0172     output_filename += (isTypeA) ? "_TypeA" : "_AllSensor";
0173     output_filename += Form("_VtxZ%.0f", fabs(cut_INTTvtxZ.first - cut_INTTvtxZ.second)/2.); // todo : check this
0174     output_filename += Form("_Mbin%d",SelectedMbin);
0175
0176     output_filename += output_file_name_suffix;
0177
0178     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0179 }
0180
0181 void PreparedNdEtaClusEach::PrepareOutPutRootFile()
0182 {
0183     file_out = new TFile(Form("%s/%s",output_directory.c_str(),output_filename.c_str()), "RECREATE");
0184 }
0185
0186 void PreparedNdEtaClusEach::PrepareHist()
0187 {
0188     hstack2D_NClusEtaVtxZ_map.clear();
0189
0190     hstack2D_NClusEtaVtxZ_map.insert( std::make_pair(
0191             Form("hstack2D_inner_NClusEtaVtxZ"),
0192             new THStack(Form("hstack2D_inner_NClusEtaVtxZ"), Form("hstack2D_inner_NClusEtaVtxZ;Cluster #eta (inner);INTT vtxZ [cm]"))
0193         )
0194     );
0195
0196     // note : normal fine bin
0197     hstack2D_NClusEtaVtxZ_map.insert( std::make_pair(
0198             Form("hstack2D_inner_NClusEtaVtxZ_FineBin"),
0199             new THStack(Form("hstack2D_inner_NClusEtaVtxZ_FineBin"), Form("hstack2D_inner_NClusEtaVtxZ_FineBin;Cluster #eta (inner);INTT vtxZ [cm]"))
0200         )
0201     );
0202
0203     hstack2D_NClusEtaVtxZ_map.insert( std::make_pair(
0204             Form("hstack2D_outer_NClusEtaVtxZ"),
0205             new THStack(Form("hstack2D_outer_NClusEtaVtxZ"), Form("hstack2D_outer_NClusEtaVtxZ;Cluster #eta (outer);INTT vtxZ [cm]"))
0206         )
0207     );
0208
0209     // note : normal fine bin
0210     hstack2D_NClusEtaVtxZ_map.insert( std::make_pair(
0211             Form("hstack2D_outer_NClusEtaVtxZ_FineBin"),
0212             new THStack(Form("hstack2D_outer_NClusEtaVtxZ_FineBin"), Form("hstack2D_outer_NClusEtaVtxZ_FineBin;Cluster #eta (outer);INTT vtxZ [cm]"))
0213         )
0214     );
0215
0216     // if (isTypeA){
0217     //     std::cout<<111<<std::endl;
0218     //     for (auto &pair : hstack2D_NClusEtaVtxZ_map){
0219     //         std::string X_title = pair.second->GetXaxis()->GetTitle();
0220     //         X_title = X_title.substr(0, X_title.find_last_of(")"));
0221     //         X_title += ", TypeA)";
0222
0223     //         pair.second->GetXaxis()->SetTitle(X_title.c_str());
0224     //     }
0225     //     std::cout<<222<<std::endl;
0226     // }
0227
0228     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0229
0230     hstack1D_NClusEta_map.clear();
0231
0232     hstack1D_NClusEta_map.insert(
0233         std::make_pair(
0234             "hstack1D_inner_NClusEta",
0235             new THStack("hstack1D_inner_NClusEta", "hstack1D_inner_NClusEta;Cluster #eta (inner);Entries")
0236         )
0237     );
0238
0239     hstack1D_NClusEta_map.insert(
0240         std::make_pair(
0241             "hstack1D_outer_NClusEta",
0242             new THStack("hstack1D_outer_NClusEta", "hstack1D_outer_NClusEta;Cluster #eta (outer);Entries")
0243         )
0244     );
0245
0246     // if (isTypeA){
0247     //     for (auto &pair : hstack1D_NClusEta_map){
0248     //         std::string X_title = pair.second->GetXaxis()->GetTitle();
0249     //         X_title = X_title.substr(0, X_title.find_last_of(")"));
0250     //         X_title += ", TypeA)";
0251
0252     //         pair.second->GetXaxis()->SetTitle(X_title.c_str());
0253     //     }
0254     // }
0255
0256     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0257
0258     if (runnumber == -1){
0259         hstack2D_TrueEtaVtxZ_map.clear();
0260
0261         // note : centrality inclusive
0262         hstack2D_TrueEtaVtxZ_map.insert( std::make_pair(
0263                 Form("hstack2D_TrueEtaVtxZ"),
0264                 new THStack(Form("hstack2D_TrueEtaVtxZ"), Form("hstack2D_TrueEtaVtxZ;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0265             )
0266         );
0267
0268         hstack2D_TrueEtaVtxZ_map.insert( std::make_pair(
0269                 Form("hstack2D_TrueEtaVtxZ_FineBin"),
0270                 new THStack(Form("hstack2D_TrueEtaVtxZ_FineBin"), Form("hstack2D_TrueEtaVtxZ_FineBin;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0271             )
0272         );
0273     }
0274
0275     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0276
0277     if (runnumber == -1){
0278         hstack1D_TrueEta_map.clear();
0279
0280         hstack1D_TrueEta_map.insert( std::make_pair(
0281             Form("hstack1D_TrueEta"),
0282             new THStack(Form("hstack1D_TrueEta"), Form("hstack1D_TrueEta;PHG4Particle #eta;Entries"))
0283         )
0284         );
0285
0286     }
0287
0288 }
0289
0290 // note : return the bin index of vtxZ, Mbin, typeA, inner, finebin
0291 std::tuple<int, int, int, int, int> PreparedNdEtaClusEach::GetHistStringInfo(std::string hist_name)
0292 {
0293
0294     // std::cout<<"In GetHistStringInfo(), input name: "<<hist_name<<std::endl;
0295
0296     int vtxz_bin;
0297     int Mbin;
0298
0299     std::string vtxz_bin_str = "";
0300     std::string Mbin_str = "";
0301
0302     // note : vtxZ bin
0303     if (hist_name.find("_VtxZ") == std::string::npos) {vtxz_bin = -1;}
0304     else {
0305         for (int i = 0; i < hist_name.size(); i++){
0306             int start_index = hist_name.find("_VtxZ") + 5 + i;
0307             if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0308                 vtxz_bin_str += hist_name[start_index];
0309             }
0310             else {
0311                 break;
0312             }
0313
0314             if (start_index == hist_name.size() - 1){
0315                 break;
0316             }
0317         }
0318
0319         vtxz_bin = (vtxz_bin_str.length() == 0) ? -1 : std::stoi(vtxz_bin_str);
0320     }
0321
0322
0323
0324     // note : Mbin
0325     if (hist_name.find("_Mbin") == std::string::npos) {Mbin = -1;}
0326     else {
0327         for (int i = 0; i < hist_name.size(); i++){
0328             int start_index = hist_name.find("_Mbin") + 5 + i;
0329             if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0330                 Mbin_str += hist_name[start_index];
0331             }
0332             else {
0333                 break;
0334             }
0335
0336             if (start_index == hist_name.size() - 1){
0337                 break;
0338             }
0339         }
0340
0341         Mbin = (Mbin_str.length() == 0) ? -1 : std::stoi(Mbin_str);
0342     }
0343
0344     // note : typeA, inner, finebin
0345     int typeA = (hist_name.find("_typeA") != std::string::npos) ? 1 : 0;
0346     int isInner = (hist_name.find("_inner") != std::string::npos) ? 1 : 0;
0347     int finebin = (hist_name.find("_FineBin") != std::string::npos) ? 1 : 0;
0348
0349
0350     return std::make_tuple(vtxz_bin, Mbin, typeA, isInner, finebin);
0351 }
0352
0353 void PreparedNdEtaClusEach::PrepareStacks()
0354 {
0355     std::cout<<"In PrepareStacks()"<<std::endl;
0356
0357     int l_vtxz_bin, l_Mbin;
0358     int l_typeA, l_inner, l_finebin;
0359
0360     for (auto &pair : h1D_input_map)
0361     {
0362         std::tie(l_vtxz_bin, l_Mbin, l_typeA, l_inner, l_finebin) = GetHistStringInfo(pair.first);
0363
0364         // note : the vtxz_bin of the hist is not in the interest
0365         if (vtxZ_index_map.find(l_vtxz_bin) == vtxZ_index_map.end()){
0366             continue;
0367         }
0368
0369         if (pair.first.find("h1D") != std::string::npos && pair.first.find("_NClusEta") != std::string::npos && l_Mbin != -1){
0370             std::cout<<"---in hstack1D--- "<<pair.first<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", inner: "<<l_inner<<", finebin: "<<l_finebin<<std::endl;
0371             std::string inner_outer_str = (l_inner) ? "inner" : "outer";
0372
0373
0374             if (isTypeA == l_typeA){ // note : {isTypeA == 0 -> inclusive, l_typeA should be zero}, {isTypeA == 1 -> typeA, l_typeA should be 1}
0375                 if (SelectedMbin == 100) // note : inclusive centrality
0376                 {
0377                     pair.second->SetFillColor(ROOT_color_code[hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->GetNhists() % ROOT_color_code.size()]);
0378                     hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->Add(pair.second);
0379                 }
0380                 // todo: need to check if the centrality bin is changed
0381                 else if (SelectedMbin == Semi_inclusive_Mbin * 10){
0382                     if (l_Mbin > Semi_inclusive_Mbin) {continue;}
0383                     pair.second->SetFillColor(ROOT_color_code[hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->GetNhists() % ROOT_color_code.size()]);
0384                     hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->Add(pair.second);
0385                 }
0386                 else if (SelectedMbin == l_Mbin)
0387                 {
0388                     pair.second->SetFillColor(ROOT_color_code[hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->GetNhists() % ROOT_color_code.size()]);
0389                     hstack1D_NClusEta_map[Form("hstack1D_%s_NClusEta",inner_outer_str.c_str())]->Add(pair.second);
0390                 }
0391                 else
0392                 {
0393                     std::cout<<"wtf the Mbin: "<<pair.first<<std::endl;
0394                 }
0395             }
0396         }
0397
0398         if (pair.first.find("h1D") != std::string::npos && pair.first.find("_TrueEta") != std::string::npos && l_Mbin != -1){
0399
0400             if (SelectedMbin == 100){
0401                 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0402                 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0403             }
0404
0405             else if (SelectedMbin == Semi_inclusive_Mbin * 10){
0406                 if (l_Mbin > Semi_inclusive_Mbin) {continue;}
0407
0408                 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0409                 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0410             }
0411
0412             else if (SelectedMbin == l_Mbin){
0413
0414                 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0415                 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0416             }
0417             else {
0418                 std::cout<<"wtf_TruthEta_inclusive_to_Mbin: "<<pair.first<<std::endl;
0419             }
0420         }
0421     }
0422
0423     // Division:--------------------------------------------------------------------------------------------------------------------------------------------------------------------
0424
0425     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0426     {
0427         std::string isTypeA_str = (isTypeA) ? "_typeA" : "";
0428
0429         if (SelectedMbin == 100)
0430         {
0431             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0432             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0433             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0434             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0435         }
0436         else if (SelectedMbin == Semi_inclusive_Mbin * 10){
0437             if (Mbin > Semi_inclusive_Mbin) {continue;}
0438
0439             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0440             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0441             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0442             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0443         }
0444         else if (SelectedMbin == Mbin)
0445         {
0446             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0447             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d", isTypeA_str.c_str(), Mbin)]);
0448             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_inner_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_inner_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0449             hstack2D_NClusEtaVtxZ_map[Form("hstack2D_outer_NClusEtaVtxZ_FineBin")] -> Add(h2D_input_map[Form("h2D%s_outer_NClusEta_Mbin%d_FineBin", isTypeA_str.c_str(), Mbin)]);
0450         }
0451     }
0452
0453     // Division:--------------------------------------------------------------------------------------------------------------------------------------------------------------------
0454     // note : for hstack2D_TrueEtaVtxZ_map
0455     // Form("hstack2D_TrueEtaVtxZ"),
0456     // Form("hstack2D_TrueEtaVtxZ_FineBin"),
0457     if (runnumber == -1)
0458     {
0459         for (auto &pair : h2D_input_map)
0460         {
0461             std::tie(l_vtxz_bin, l_Mbin, l_typeA, l_inner, l_finebin) = GetHistStringInfo(pair.first);
0462
0463             if (pair.first.find("h2D_TrueEtaVtxZ") != std::string::npos && l_Mbin != -1){
0464
0465                 std::cout<<" Truth for TrueEtaVtxZ ----- "<<pair.first<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", inner: "<<l_inner<<", finebin: "<<l_finebin<<std::endl;
0466
0467                 std::string finebin_str = (l_finebin) ? "_FineBin" : "";
0468
0469                 if (SelectedMbin == 100){
0470                     hstack2D_TrueEtaVtxZ_map[Form("hstack2D_TrueEtaVtxZ%s",finebin_str.c_str())] -> Add(pair.second);
0471                 }
0472                 else if (SelectedMbin == Semi_inclusive_Mbin * 10){
0473                     if (l_Mbin > Semi_inclusive_Mbin) {continue;}
0474
0475                     hstack2D_TrueEtaVtxZ_map[Form("hstack2D_TrueEtaVtxZ%s",finebin_str.c_str())] -> Add(pair.second);
0476                 }
0477                 else if (SelectedMbin == l_Mbin){
0478                     hstack2D_TrueEtaVtxZ_map[Form("hstack2D_TrueEtaVtxZ%s",finebin_str.c_str())] -> Add(pair.second);
0479                 }
0480                 else {
0481                     std::cout<<"wtf_TrueEtaVtxZ_inclusive_to_Mbin: "<<pair.first<<std::endl;
0482                 }
0483
0484             }
0485
0486
0487         }
0488     }
0489
0490     h1D_inner_NClusEta = (TH1D*)((TH1D*) hstack1D_NClusEta_map["hstack1D_inner_NClusEta"] -> GetStack() -> Last())->Clone("h1D_inner_NClusEta");
0491     h1D_outer_NClusEta = (TH1D*)((TH1D*) hstack1D_NClusEta_map["hstack1D_outer_NClusEta"] -> GetStack() -> Last())->Clone("h1D_outer_NClusEta");
0492
0493     h1D_inner_NClusEta -> SetFillColorAlpha(2,0);
0494     h1D_outer_NClusEta -> SetFillColorAlpha(2,0);
0495 }
0496
0497 void PreparedNdEtaClusEach::PreparedNdEtaHist()
0498 {
0499     if (runnumber == -1){
0500         h1D_TruedNdEta_map.clear();
0501
0502         h1D_TruedNdEta_map.insert(
0503             std::make_pair(
0504                 "h1D_TruedNdEta",
0505                 (TH1D*)((TH1D*)hstack1D_TrueEta_map["hstack1D_TrueEta"] -> GetStack() -> Last()) -> Clone("h1D_TruedNdEta")
0506             )
0507         );
0508         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetFillColorAlpha(3,0);
0509         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetName("h1D_TruedNdEta");
0510         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetTitle("h1D_TruedNdEta");
0511         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Scale(1. / get_EvtCount(h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"], SelectedMbin));
0512         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Scale(1. / double(h1D_TruedNdEta_map["h1D_TruedNdEta"]->GetBinWidth(1)));
0513
0514     } // note : end of truth
0515
0516     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0517
0518
0519     h1D_inner_NClusEtaPerEvt = (TH1D*) h1D_inner_NClusEta -> Clone("h1D_inner_NClusEtaPerEvt");
0520     h1D_inner_NClusEtaPerEvt -> SetName("h1D_inner_NClusEtaPerEvt");
0521     h1D_inner_NClusEtaPerEvt -> SetTitle("h1D_inner_NClusEtaPerEvt");
0522     h1D_inner_NClusEtaPerEvt -> Scale(1. / get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], SelectedMbin));
0523     h1D_inner_NClusEtaPerEvt -> Scale(1. / double(h1D_inner_NClusEtaPerEvt->GetBinWidth(1)));
0524
0525
0526     h1D_outer_NClusEtaPerEvt = (TH1D*) h1D_outer_NClusEta -> Clone("h1D_outer_NClusEtaPerEvt");
0527     h1D_outer_NClusEtaPerEvt -> SetName("h1D_outer_NClusEtaPerEvt");
0528     h1D_outer_NClusEtaPerEvt -> SetTitle("h1D_outer_NClusEtaPerEvt");
0529     h1D_outer_NClusEtaPerEvt -> Scale(1. / get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], SelectedMbin));
0530     h1D_outer_NClusEtaPerEvt -> Scale(1. / double(h1D_outer_NClusEtaPerEvt->GetBinWidth(1)));
0531
0532     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0533
0534     h1D_inner_NClusEtaPerEvtPostAC = (TH1D*) h1D_inner_NClusEtaPerEvt -> Clone("h1D_inner_NClusEtaPerEvtPostAC");
0535     h1D_inner_NClusEtaPerEvtPostAC -> SetName("h1D_inner_NClusEtaPerEvtPostAC");
0536     h1D_inner_NClusEtaPerEvtPostAC -> SetTitle("h1D_inner_NClusEtaPerEvtPostAC");
0537
0538
0539     h1D_outer_NClusEtaPerEvtPostAC = (TH1D*) h1D_outer_NClusEtaPerEvt -> Clone("h1D_outer_NClusEtaPerEvtPostAC");
0540     h1D_outer_NClusEtaPerEvtPostAC -> SetName("h1D_outer_NClusEtaPerEvtPostAC");
0541     h1D_outer_NClusEtaPerEvtPostAC -> SetTitle("h1D_outer_NClusEtaPerEvtPostAC");
0542
0543     // Division:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0544
0545     if (h1D_alpha_correction_map_in.size() != 0 && ApplyAlphaCorr == true)
0546     {
0547         h1D_inner_NClusEtaPerEvtPostAC -> Sumw2(true);
0548         h1D_inner_NClusEtaPerEvtPostAC -> Divide(h1D_alpha_correction_map_in[alpha_correction_name_map[0]]);
0549
0550         h1D_outer_NClusEtaPerEvtPostAC -> Sumw2(true);
0551         h1D_outer_NClusEtaPerEvtPostAC -> Divide(h1D_alpha_correction_map_in[alpha_correction_name_map[1]]);
0552
0553         if (cut_EtaRange_pair.first == cut_EtaRange_pair.first){
0554             for (int eta_bin = 0; eta_bin < h1D_inner_NClusEtaPerEvtPostAC -> GetNbinsX(); eta_bin++)
0555             {
0556                 if (
0557                     h1D_inner_NClusEtaPerEvtPostAC -> GetXaxis() -> GetBinUpEdge(eta_bin + 1) < cut_EtaRange_pair.first ||
0558                     h1D_inner_NClusEtaPerEvtPostAC -> GetXaxis() -> GetBinLowEdge(eta_bin + 1) > cut_EtaRange_pair.second
0559                 ){
0560                     h1D_inner_NClusEtaPerEvtPostAC -> SetBinContent(eta_bin + 1, 0);
0561                     h1D_inner_NClusEtaPerEvtPostAC -> SetBinError(eta_bin + 1, 0);
0562
0563                     h1D_outer_NClusEtaPerEvtPostAC -> SetBinContent(eta_bin + 1, 0);
0564                     h1D_outer_NClusEtaPerEvtPostAC -> SetBinError(eta_bin + 1, 0);
0565                 }
0566             }
0567         }
0568
0569
0570     }
0571     else if (h1D_alpha_correction_map_in.size() == 0 && ApplyAlphaCorr == true){
0572         std::cout<<"wtf_PreparedNdEtaHist: h1D_alpha_correction_map_in.size() == 0"<<std::endl;
0573         exit(1);
0574     }
0575
0576 }
0577
0578 void PreparedNdEtaClusEach::DeriveAlphaCorrection()
0579 {
0580     if (runnumber != -1) {
0581         std::cout<<"wtf_DeriveAlphaCorrection: runnumber != -1"<<std::endl;
0582         return;
0583     }
0584
0585     std::cout<<"In DeriveAlphaCorrection()"<<std::endl;
0586
0587     h1D_alpha_correction_map_out.insert(
0588         std::make_pair(
0589             alpha_correction_name_map[0],
0590             new TH1D(alpha_correction_name_map[0].c_str(), Form("%s;#eta;#alpha corrections", alpha_correction_name_map[0].c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0591         )
0592     );
0593
0594     h1D_alpha_correction_map_out.insert(
0595         std::make_pair(
0596             alpha_correction_name_map[1],
0597             new TH1D(alpha_correction_name_map[1].c_str(), Form("%s;#eta;#alpha corrections",alpha_correction_name_map[1].c_str()), nEtaBin, EtaEdge_min, EtaEdge_max)
0598         )
0599     );
0600
0601
0602     h1D_alpha_correction_map_out[alpha_correction_name_map[0]] -> Sumw2(true);
0603     h1D_alpha_correction_map_out[alpha_correction_name_map[1]] -> Sumw2(true);
0604
0605     h1D_alpha_correction_map_out[alpha_correction_name_map[0]] -> Divide(h1D_inner_NClusEtaPerEvtPostAC, h1D_TruedNdEta_map[Form("h1D_TruedNdEta")]);
0606     h1D_alpha_correction_map_out[alpha_correction_name_map[1]] -> Divide(h1D_outer_NClusEtaPerEvtPostAC, h1D_TruedNdEta_map[Form("h1D_TruedNdEta")]);
0607 }
0608
0609 void PreparedNdEtaClusEach::EndRun()
0610 {
0611     file_out -> cd();
0612
0613     ((TH2D*) hstack2D_NClusEtaVtxZ_map["hstack2D_inner_NClusEtaVtxZ"]->GetStack()->Last()) -> Write("h2D_inner_NClusEtaVtxZ");
0614     ((TH2D*) hstack2D_NClusEtaVtxZ_map["hstack2D_inner_NClusEtaVtxZ_FineBin"]->GetStack()->Last()) -> Write("h2D_inner_NClusEtaVtxZ_FineBin");
0615     ((TH2D*) hstack2D_NClusEtaVtxZ_map["hstack2D_outer_NClusEtaVtxZ"]->GetStack()->Last()) -> Write("h2D_outer_NClusEtaVtxZ");
0616     ((TH2D*) hstack2D_NClusEtaVtxZ_map["hstack2D_outer_NClusEtaVtxZ_FineBin"]->GetStack()->Last()) -> Write("h2D_outer_NClusEtaVtxZ_FineBin");
0617     h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"] -> Write();
0618
0619     if (runnumber == -1){
0620         h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"] -> Write();
0621
0622         ((TH2D*)hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"]->GetStack()->Last()) -> Write("h2D_TrueEtaVtxZ");
0623         ((TH2D*)hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"]->GetStack()->Last()) -> Write("h2D_TrueEtaVtxZ_FineBin");
0624
0625
0626         ((TH1D*) hstack1D_TrueEta_map["hstack1D_TrueEta"]->GetStack()->Last()) -> Write("h1D_TrueEta");
0627         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Write();
0628     }
0629
0630     h1D_inner_NClusEta -> Write();
0631     h1D_inner_NClusEtaPerEvt -> Write();
0632     h1D_inner_NClusEtaPerEvtPostAC -> Write();
0633
0634     h1D_outer_NClusEta -> Write();
0635     h1D_outer_NClusEtaPerEvt -> Write();
0636     h1D_outer_NClusEtaPerEvtPostAC -> Write();
0637
0638
0639     for (auto &pair : h1D_alpha_correction_map_out){
0640         pair.second -> Write();
0641     }
0642
0643     if (h1D_alpha_correction_map_in.size() != 0){
0644         for (auto &pair : h1D_alpha_correction_map_in){
0645             pair.second -> Write(("Used_" + pair.first).c_str());
0646         }
0647     }
0648
0649
0650
0651     for (auto &pair : hstack2D_NClusEtaVtxZ_map){
0652         pair.second -> Write();
0653     }
0654     for (auto &pair : hstack1D_NClusEta_map){
0655         pair.second -> Write();
0656     }
0657
0658     if (runnumber == -1){
0659         for (auto &pair : hstack2D_TrueEtaVtxZ_map){
0660             pair.second -> Write();
0661         }
0662         for (auto &pair : hstack1D_TrueEta_map){
0663             pair.second -> Write();
0664         }
0665     }
0666
0667     if (runnumber == -1){
0668         c1 -> Clear();
0669         c1 -> cd();
0670
0671         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetFillColorAlpha(0,0);
0672         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetMinimum(0.);
0673         h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Draw("hist");
0674         h1D_inner_NClusEtaPerEvtPostAC->SetLineColor(2);
0675         h1D_inner_NClusEtaPerEvtPostAC->Draw("ep same");
0676         h1D_outer_NClusEtaPerEvtPostAC->SetLineColor(8);
0677         h1D_outer_NClusEtaPerEvtPostAC->Draw("ep same");
0678
0679         c1 -> Write("c1_comp_h1D_dNdEta");
0680         c1 -> Clear();
0681     }
0682
0683 }
0684
0685
0686 double PreparedNdEtaClusEach::get_EvtCount(TH2D * hist_in, int centrality_bin_in)
0687 {
0688     double count = 0;
0689
0690     // note : X : vtxZ
0691     // note : Y : centrality
0692
0693     if (centrality_bin_in == 100){
0694         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
0695
0696             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
0697                 continue;
0698             }
0699
0700             for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){ // note : centrality bin
0701                 count += hist_in -> GetBinContent(xi, yi);
0702             }
0703         }
0704     }
0705     // todo:
0706     else if (centrality_bin_in == Semi_inclusive_Mbin * 10){
0707         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
0708
0709             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
0710                 continue;
0711             }
0712
0713             for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){ // note : centrality bin
0714                 if (yi-1 <= Semi_inclusive_Mbin){
0715                     count += hist_in -> GetBinContent(xi, yi);
0716                 }
0717             }
0718         }
0719     }
0720     else {
0721         for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){ // note : vtxZ bin
0722
0723             if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
0724                 continue;
0725             }
0726
0727             count += hist_in -> GetBinContent(xi, centrality_bin_in + 1);
0728         }
0729     }
0730
0731     return count;
0732
0733 }