File indexing completed on 2025-08-06 08:12:42
0001 #include "PreparedNdEtaEach.h"
0002
0003 PreparedNdEtaEach::PreparedNdEtaEach(
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 ApplyGeoAccCorr_in,
0013
0014 bool isTypeA_in,
0015 std::pair<double,double> cut_INTTvtxZ_in,
0016 int SelectedMbin_in
0017 ) : process_id(process_id_in),
0018 runnumber(runnumber_in),
0019 input_directory(input_directory_in),
0020 input_file_name(input_file_name_in),
0021 output_directory(output_directory_in),
0022 output_file_name_suffix(output_file_name_suffix_in),
0023
0024 ApplyAlphaCorr(ApplyAlphaCorr_in),
0025 ApplyGeoAccCorr(ApplyGeoAccCorr_in),
0026
0027 isTypeA(isTypeA_in),
0028 cut_INTTvtxZ(cut_INTTvtxZ_in),
0029 SelectedMbin(SelectedMbin_in)
0030 {
0031 PrepareInputRootFie();
0032
0033 if (SelectedMbin >= nCentrality_bin && SelectedMbin != Semi_inclusive_interval && SelectedMbin != 100){
0034 std::cout << "Error : SelectedMbin is out of range" << std::endl;
0035 std::cout << "SelectedMbin : " << SelectedMbin <<", but : nCentrality_bin: "<< nCentrality_bin << std::endl;
0036 exit(1);
0037 }
0038
0039 vtxZ_index_map = GetVtxZIndexMap();
0040
0041 PrepareOutPutFileName();
0042
0043 PrepareOutPutRootFile();
0044
0045 PrepareHistFits();
0046
0047 c1 = new TCanvas("c1", "c1", 950, 800);
0048 c1 -> Print(Form("%s/%s(", output_directory.c_str(), output_filename_pdf.c_str()));
0049
0050 h1D_alpha_correction_map_in.clear();
0051 h1D_alpha_correction_map_out.clear();
0052
0053 h2D_GeoAccCorr_map.clear();
0054
0055 std::string isTypeA_str = (isTypeA) ? "_typeA" : "";
0056 GeoAccCorr_name_map = {
0057 Form("h2D%s_BestPair_EtaVtxZ_Rebin_Selected",isTypeA_str.c_str()),
0058 Form("h2D%s_BestPair_EtaVtxZ_Rebin_Selected",isTypeA_str.c_str())
0059 };
0060 }
0061
0062 void PreparedNdEtaEach::PrepareInputRootFie()
0063 {
0064 file_in = TFile::Open(Form("%s/%s",input_directory.c_str(),input_file_name.c_str()));
0065 if (file_in == nullptr)
0066 {
0067 std::cout << "Error: file_in can not be opened" << std::endl;
0068 exit(1);
0069 }
0070
0071 for (TObject* keyAsObj : *file_in->GetListOfKeys())
0072 {
0073 auto key = dynamic_cast<TKey*>(keyAsObj);
0074 std::string hist_name = key->GetName();
0075 std::string class_name = key->GetClassName();
0076
0077 if (class_name == "TH2D")
0078 {
0079 h2D_input_map[hist_name.c_str()] = (TH2D*) file_in -> Get( hist_name.c_str() );
0080 }
0081 else if (class_name == "TH1D")
0082 {
0083 h1D_input_map[hist_name.c_str()] = (TH1D*) file_in -> Get( hist_name.c_str() );
0084 h1D_input_map[hist_name.c_str()] -> SetLineColor(kBlack);
0085 h1D_input_map[hist_name.c_str()] -> SetLineWidth(2);
0086 }
0087 }
0088 std::cout<<"In PrepareInputRootFie(), number of input TH1D: "<< h1D_input_map.size() <<std::endl;
0089 std::cout<<"In PrepareInputRootFie(), number of input TH2D: "<< h2D_input_map.size() <<std::endl;
0090
0091 tree_in = (TTree*) file_in -> Get("tree_par");
0092
0093 centrality_edges = 0;
0094 cut_GoodProtoTracklet_DeltaPhi = 0;
0095
0096 tree_in -> SetBranchAddress("centrality_edges", ¢rality_edges);
0097 tree_in -> SetBranchAddress("nCentrality_bin", &nCentrality_bin);
0098
0099 tree_in -> SetBranchAddress("CentralityFineEdge_min", &CentralityFineEdge_min);
0100 tree_in -> SetBranchAddress("CentralityFineEdge_max", &CentralityFineEdge_max);
0101 tree_in -> SetBranchAddress("nCentralityFineBin", &nCentralityFineBin);
0102
0103 tree_in -> SetBranchAddress("EtaEdge_min", &EtaEdge_min);
0104 tree_in -> SetBranchAddress("EtaEdge_max", &EtaEdge_max);
0105 tree_in -> SetBranchAddress("nEtaBin", &nEtaBin);
0106
0107 tree_in -> SetBranchAddress("VtxZEdge_min", &VtxZEdge_min);
0108 tree_in -> SetBranchAddress("VtxZEdge_max", &VtxZEdge_max);
0109 tree_in -> SetBranchAddress("nVtxZBin", &nVtxZBin);
0110
0111 tree_in -> SetBranchAddress("DeltaPhiEdge_min", &DeltaPhiEdge_min);
0112 tree_in -> SetBranchAddress("DeltaPhiEdge_max", &DeltaPhiEdge_max);
0113 tree_in -> SetBranchAddress("nDeltaPhiBin", &nDeltaPhiBin);
0114
0115 tree_in -> SetBranchAddress("cut_GoodProtoTracklet_DeltaPhi", &cut_GoodProtoTracklet_DeltaPhi);
0116
0117 tree_in -> GetEntry(0);
0118
0119 std::cout<<"------------------------------------------------------------"<<std::endl;
0120 std::cout<<"In PrepareInputRootFile(), nEntries : "<<tree_in->GetEntries()<<std::endl;
0121 std::cout<<"nCentrality_bin : "<<nCentrality_bin<<std::endl;
0122 for(int i = 0; i < nCentrality_bin; i++)
0123 {
0124 std::cout<<"centrality_edges["<<i<<"] : "<<centrality_edges->at(i)<<std::endl;
0125 }
0126 std::cout<<"CentralityFineBin : "<<nCentralityFineBin<<", "<<CentralityFineEdge_min<<", "<<CentralityFineEdge_max<<std::endl;
0127 std::cout<<"EtaBin : "<<nEtaBin<<", "<<EtaEdge_min<<", "<<EtaEdge_max<<std::endl;
0128 std::cout<<"VtxZBin : "<<nVtxZBin<<", "<<VtxZEdge_min<<", "<<VtxZEdge_max<<std::endl;
0129 std::cout<<"DeltaPhiBin : "<<nDeltaPhiBin<<", "<<DeltaPhiEdge_min<<", "<<DeltaPhiEdge_max<<std::endl;
0130 std::cout<<"DeltaPhi signal region: "<<cut_GoodProtoTracklet_DeltaPhi->first<<" ~ "<<cut_GoodProtoTracklet_DeltaPhi->second<<std::endl;
0131 std::cout<<"------------------------------------------------------------"<<std::endl;
0132
0133 }
0134
0135 std::map<int, int> PreparedNdEtaEach::GetVtxZIndexMap()
0136 {
0137 if(h2D_input_map.find("h2D_RecoEvtCount_vtxZCentrality") == h2D_input_map.end()){
0138 std::cout<<"Error: h2D_RecoEvtCount_vtxZCentrality not found"<<std::endl;
0139 exit(1);
0140 }
0141
0142 std::map<int, int> vtxZIndexMap; vtxZIndexMap.clear();
0143
0144 for (int y_i = 1; y_i <= h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"]->GetNbinsY(); y_i++){
0145
0146 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){
0147 continue;
0148 }
0149 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){
0150 continue;
0151 }
0152
0153 vtxZIndexMap.insert(
0154 std::make_pair(
0155 y_i - 1,
0156 y_i - 1
0157 )
0158 );
0159 }
0160
0161 std::cout<<"The selected INTT vtxZ bin : [";
0162 for (auto &pair : vtxZIndexMap){
0163 std::cout<<pair.first<<", ";
0164 }
0165 std::cout<<"]"<<std::endl;
0166
0167 for (auto &pair : vtxZIndexMap){
0168 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;
0169 }
0170
0171 return vtxZIndexMap;
0172 }
0173
0174 void PreparedNdEtaEach::PrepareOutPutFileName()
0175 {
0176 std::string job_index = std::to_string( process_id );
0177 int job_index_len = 5;
0178 job_index.insert(0, job_index_len - job_index.size(), '0');
0179
0180 std::string runnumber_str = std::to_string( runnumber );
0181 if (runnumber != -1){
0182 int runnumber_str_len = 8;
0183 runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0184 }
0185
0186 if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0187 output_file_name_suffix = "_" + output_file_name_suffix;
0188 }
0189
0190 output_filename = "PreparedNdEtaEach";
0191 output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0192
0193 output_filename += (ApplyAlphaCorr) ? "_AlphaCorr" : "";
0194 output_filename += (ApplyGeoAccCorr) ? "_GeoAccCorr" : "";
0195
0196 output_filename += (isTypeA) ? "_TypeA" : "_AllSensor";
0197 output_filename += Form("_VtxZ%.0f", fabs(cut_INTTvtxZ.first - cut_INTTvtxZ.second)/2.);
0198 output_filename += Form("_Mbin%d",SelectedMbin);
0199
0200 output_filename += output_file_name_suffix;
0201
0202 output_filename_DeltaPhi = output_filename;
0203 output_filename_dNdEta = output_filename;
0204 output_filename_pdf = output_filename;
0205
0206
0207 output_filename_DeltaPhi += (runnumber != -1) ? Form("_%s_%s_DeltaPhi.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s_DeltaPhi.root",job_index.c_str());
0208 output_filename_dNdEta += (runnumber != -1) ? Form("_%s_%s_dNdEta.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s_dNdEta.root",job_index.c_str());
0209 output_filename_pdf += (runnumber != -1) ? Form("_%s_%s.pdf",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.pdf",job_index.c_str());
0210 }
0211
0212 void PreparedNdEtaEach::PrepareOutPutRootFile()
0213 {
0214 file_out_DeltaPhi = new TFile(Form("%s/%s",output_directory.c_str(),output_filename_DeltaPhi.c_str()), "RECREATE");
0215 file_out_dNdEta = new TFile(Form("%s/%s",output_directory.c_str(),output_filename_dNdEta.c_str()), "RECREATE");
0216 }
0217
0218 void PreparedNdEtaEach::PrepareHistFits()
0219 {
0220 hstack1D_BestPair_ClusPhiSize = new THStack("hstack1D_BestPair_ClusPhiSize","hstack1D_BestPair_ClusPhiSize");
0221 hstack1D_BestPair_ClusAdc = new THStack("hstack1D_BestPair_ClusAdc","hstack1D_BestPair_ClusAdc");
0222 hstack1D_BestPair_DeltaPhi = new THStack("hstack1D_BestPair_DeltaPhi","hstack1D_BestPair_DeltaPhi");
0223 hstack1D_BestPair_DeltaEta = new THStack("hstack1D_BestPair_DeltaEta","hstack1D_BestPair_DeltaEta");
0224 hstack2D_BestPairEtaVtxZ = new THStack("hstack2D_BestPairEtaVtxZ","hstack2D_BestPairEtaVtxZ");
0225 hstack2D_BestPairEtaVtxZ_FineBin = new THStack("hstack2D_BestPairEtaVtxZ_FineBin","hstack2D_BestPairEtaVtxZ_FineBin");
0226
0227
0228 hstack2D_GoodProtoTracklet_map.clear();
0229
0230 hstack2D_GoodProtoTracklet_map.insert( std::make_pair(
0231 Form("hstack2D_GoodProtoTracklet_EtaVtxZ"),
0232 new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ;Pair #eta;INTT vtxZ [cm]"))
0233 )
0234 ).second;
0235
0236
0237 hstack2D_GoodProtoTracklet_map.insert( std::make_pair(
0238 Form("hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"),
0239 new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin;Pair #eta;INTT vtxZ [cm]"))
0240 )
0241 ).second;
0242
0243
0244 hstack2D_GoodProtoTracklet_map.insert( std::make_pair(
0245 Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"),
0246 new THStack(Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"), Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated;Pair #eta;INTT vtxZ [cm]"))
0247 )
0248 ).second;
0249
0250
0251
0252 h1D_BestPair_RecoTrackletEta_map.clear();
0253 h1D_BestPair_RecoTrackletEta_map.insert( std::make_pair(
0254 Form("h1D_BestPair_RecoTrackletEta"),
0255 new TH1D(Form("h1D_BestPair_RecoTrackletEta"), Form("h1D_BestPair_RecoTrackletEta;Tracklet #eta;Entries"), nEtaBin, EtaEdge_min, EtaEdge_max)
0256 )
0257 ).second;
0258
0259 h1D_RotatedBkg_RecoTrackletEta_map.clear();
0260 h1D_RotatedBkg_RecoTrackletEta_map.insert( std::make_pair(
0261 Form("h1D_RotatedBkg_RecoTrackletEta"),
0262 new TH1D(Form("h1D_RotatedBkg_RecoTrackletEta"), Form("h1D_RotatedBkg_RecoTrackletEta;Tracklet #eta;Entries"), nEtaBin, EtaEdge_min, EtaEdge_max)
0263 )
0264 ).second;
0265
0266 h1D_RotatedBkg_DeltaPhi_Signal_map.clear();
0267 for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++){
0268 h1D_RotatedBkg_DeltaPhi_Signal_map.insert(
0269 std::make_pair(
0270 Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin),
0271 new TH1D(Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin), Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d;Pair #Delta#phi [radian];Entries", eta_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0272 )
0273 );
0274
0275 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin)] -> Sumw2(true);
0276
0277 if (isTypeA) {
0278 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin)] -> GetXaxis() -> SetTitle("Pair #Delta#phi (type A) [radian]");
0279 }
0280
0281 }
0282
0283
0284
0285 hstack1D_DeltaPhi_map.clear();
0286 for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
0287 {
0288 hstack1D_DeltaPhi_map.insert( std::make_pair(
0289 Form("hstack1D_DeltaPhi_Eta%d", eta_bin),
0290 new THStack(Form("hstack1D_DeltaPhi_Eta%d", eta_bin), Form("hstack1D_DeltaPhi_Eta%d;Pair #Delta#phi [radian];Entries", eta_bin))
0291 )
0292 ).second;
0293
0294
0295 hstack1D_DeltaPhi_map.insert( std::make_pair(
0296 Form("hstack1D_DeltaPhi_Eta%d_rotated", eta_bin),
0297 new THStack(Form("hstack1D_DeltaPhi_Eta%d_rotated", eta_bin), Form("hstack1D_DeltaPhi_Eta%d_rotated;Pair #Delta#phi [radian];Entries", eta_bin))
0298 )
0299 ).second;
0300 }
0301
0302
0303
0304
0305
0306 f1_SigBkgPol2_Fit_map.clear();
0307 f1_SigBkgPol2_DrawSig_map.clear();
0308 f1_SigBkgPol2_DrawBkgPol2_map.clear();
0309
0310 f1_BkgPolTwo_Fit_map.clear();
0311
0312 for (auto &pair : hstack1D_DeltaPhi_map){
0313 if (pair.first.find("hstack1D_DeltaPhi") != std::string::npos){
0314
0315 std::string f1_name = pair.first + "_BkgPolTwo_Fit";
0316 f1_BkgPolTwo_Fit_map.insert(
0317 std::make_pair(
0318 f1_name,
0319 new TF1(f1_name.c_str(), "pol2", DeltaPhiEdge_min, DeltaPhiEdge_max)
0320 )
0321 ).second;
0322 f1_BkgPolTwo_Fit_map[f1_name] -> SetLineColor(6);
0323 f1_BkgPolTwo_Fit_map[f1_name] -> SetLineStyle(9);
0324 f1_BkgPolTwo_Fit_map[f1_name] -> SetLineWidth(3);
0325 f1_BkgPolTwo_Fit_map[f1_name] -> SetNpx(1000);
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352 f1_name = pair.first + "_SigBkgPol2_Fit";
0353 f1_SigBkgPol2_Fit_map.insert(
0354 std::make_pair(
0355 f1_name,
0356 new TF1(f1_name.c_str(), gaus_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 7)
0357 )
0358 ).second;
0359 f1_SigBkgPol2_Fit_map[f1_name] -> SetLineColor(46);
0360 f1_SigBkgPol2_Fit_map[f1_name] -> SetLineStyle(1);
0361 f1_SigBkgPol2_Fit_map[f1_name] -> SetNpx(1000);
0362
0363
0364 f1_name = pair.first + "_SigBkgPol2_DrawSig";
0365 f1_SigBkgPol2_DrawSig_map.insert(
0366 std::make_pair(
0367 f1_name,
0368 new TF1(f1_name.c_str(), gaus_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 4)
0369 )
0370 ).second;
0371 f1_SigBkgPol2_DrawSig_map[f1_name] -> SetLineColor(46);
0372 f1_SigBkgPol2_DrawSig_map[f1_name] -> SetLineStyle(2);
0373 f1_SigBkgPol2_DrawSig_map[f1_name] -> SetNpx(1000);
0374
0375
0376 f1_name = pair.first + "_SigBkgPol2_DrawBkgPol2";
0377 f1_SigBkgPol2_DrawBkgPol2_map.insert(
0378 std::make_pair(
0379 f1_name,
0380 new TF1(f1_name.c_str(), full_pol2_func, DeltaPhiEdge_min, DeltaPhiEdge_max, 4)
0381 )
0382 ).second;
0383 f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetLineColor(8);
0384 f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetLineStyle(2);
0385 f1_SigBkgPol2_DrawBkgPol2_map[f1_name] -> SetNpx(1000);
0386
0387 }
0388 }
0389
0390
0391
0392 if (runnumber == -1){
0393 hstack2D_TrueEtaVtxZ_map.clear();
0394
0395
0396 hstack2D_TrueEtaVtxZ_map.insert( std::make_pair(
0397 Form("hstack2D_TrueEtaVtxZ"),
0398 new THStack(Form("hstack2D_TrueEtaVtxZ"), Form("hstack2D_TrueEtaVtxZ;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0399 )
0400 ).second;
0401
0402 hstack2D_TrueEtaVtxZ_map.insert( std::make_pair(
0403 Form("hstack2D_TrueEtaVtxZ_FineBin"),
0404 new THStack(Form("hstack2D_TrueEtaVtxZ_FineBin"), Form("hstack2D_TrueEtaVtxZ_FineBin;PHG4Particle #eta;TruthPV_trig_z [cm]"))
0405 )
0406 ).second;
0407 }
0408
0409
0410
0411 if (runnumber == -1){
0412 hstack1D_TrueEta_map.clear();
0413
0414 hstack1D_TrueEta_map.insert( std::make_pair(
0415 Form("hstack1D_TrueEta"),
0416 new THStack(Form("hstack1D_TrueEta"), Form("hstack1D_TrueEta;PHG4Particle #eta;Entries"))
0417 )
0418 ).second;
0419
0420 }
0421
0422 }
0423
0424
0425 std::tuple<int, int, int, int, int, int> PreparedNdEtaEach::GetHistStringInfo(std::string hist_name)
0426 {
0427
0428
0429
0430 int eta_bin;
0431 int vtxz_bin;
0432 int Mbin;
0433
0434 std::string eta_bin_str = "";
0435 std::string vtxz_bin_str = "";
0436 std::string Mbin_str = "";
0437
0438
0439 if (hist_name.find("_Eta") == std::string::npos) {eta_bin = -1;}
0440 else {
0441 for (int i = 0; i < hist_name.size(); i++){
0442 int start_index = hist_name.find("_Eta") + 4 + i;
0443 if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0444 eta_bin_str += hist_name[start_index];
0445 }
0446 else {
0447 break;
0448 }
0449
0450 if (start_index == hist_name.size() - 1){
0451 break;
0452 }
0453 }
0454
0455 eta_bin = (eta_bin_str.length() == 0) ? -1 : std::stoi(eta_bin_str);
0456 }
0457
0458
0459
0460
0461 if (hist_name.find("_VtxZ") == std::string::npos) {vtxz_bin = -1;}
0462 else {
0463 for (int i = 0; i < hist_name.size(); i++){
0464 int start_index = hist_name.find("_VtxZ") + 5 + i;
0465 if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0466 vtxz_bin_str += hist_name[start_index];
0467 }
0468 else {
0469 break;
0470 }
0471
0472 if (start_index == hist_name.size() - 1){
0473 break;
0474 }
0475 }
0476
0477 vtxz_bin = (vtxz_bin_str.length() == 0) ? -1 : std::stoi(vtxz_bin_str);
0478 }
0479
0480
0481
0482
0483 if (hist_name.find("_Mbin") == std::string::npos) {Mbin = -1;}
0484 else {
0485 for (int i = 0; i < hist_name.size(); i++){
0486 int start_index = hist_name.find("_Mbin") + 5 + i;
0487 if (hist_name[start_index] != '_' && isdigit(hist_name[start_index])){
0488 Mbin_str += hist_name[start_index];
0489 }
0490 else {
0491 break;
0492 }
0493
0494 if (start_index == hist_name.size() - 1){
0495 break;
0496 }
0497 }
0498
0499 Mbin = (Mbin_str.length() == 0) ? -1 : std::stoi(Mbin_str);
0500 }
0501
0502
0503 int typeA = (hist_name.find("_typeA") != std::string::npos) ? 1 : 0;
0504 int rotated = (hist_name.find("_rotated") != std::string::npos) ? 1 : 0;
0505 int finebin = (hist_name.find("_FineBin") != std::string::npos) ? 1 : 0;
0506
0507
0508 return std::make_tuple(eta_bin, vtxz_bin, Mbin, typeA, rotated, finebin);
0509 }
0510
0511 void PreparedNdEtaEach::PrepareStacks()
0512 {
0513 std::cout<<"In PrepareStacks()"<<std::endl;
0514
0515 int l_eta_bin, l_vtxz_bin, l_Mbin;
0516 int l_typeA, l_rotated, l_finebin;
0517
0518 if (ApplyGeoAccCorr && h2D_GeoAccCorr_map.size() == 1){
0519 if (
0520 cut_INTTvtxZ.first < h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmin() ||
0521 cut_INTTvtxZ.second > h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmax()
0522 )
0523 {
0524 std::cout<<"Error: cut_INTTvtxZ is outside the h2D_GeoAccr map"<<std::endl;
0525 exit(1);
0526 }
0527 }
0528 else if (ApplyGeoAccCorr && h2D_GeoAccCorr_map.size() != 1){
0529 std::cout<<"Error: h2D_GeoAccCorr_map is empty"<<std::endl;
0530 exit(1);
0531 }
0532
0533 for (auto &pair : h1D_input_map)
0534 {
0535 std::tie(l_eta_bin, l_vtxz_bin, l_Mbin, l_typeA, l_rotated, l_finebin) = GetHistStringInfo(pair.first);
0536
0537
0538 if (vtxZ_index_map.find(l_vtxz_bin) == vtxZ_index_map.end()){
0539 continue;
0540 }
0541
0542 if (pair.first.find("h1D") != std::string::npos && pair.first.find("_DeltaPhi") != std::string::npos && l_Mbin != -1 && l_eta_bin != -1){
0543
0544 if (!ApplyGeoAccCorr) {std::cout<<"----- "<<pair.first<<", eta_bin: "<<l_eta_bin<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", rotated: "<<l_rotated<<", finebin: "<<l_finebin<<std::endl;}
0545
0546 double GeoAcc_correction;
0547
0548 if (ApplyGeoAccCorr && h2D_GeoAccCorr_map.size() != 0)
0549 {
0550
0551 double h1D_eta_center = h2D_input_map["h2D_GoodProtoTracklet_EtaVtxZ_Mbin0"] -> GetXaxis() -> GetBinCenter(l_eta_bin + 1);
0552 double h1D_Z_center = h2D_input_map["h2D_GoodProtoTracklet_EtaVtxZ_Mbin0"] -> GetYaxis() -> GetBinCenter(l_vtxz_bin + 1);
0553 double h1D_Z_BinWidth = h2D_input_map["h2D_GoodProtoTracklet_EtaVtxZ_Mbin0"] -> GetYaxis() -> GetBinWidth(l_vtxz_bin + 1);
0554
0555 auto temp_RotatedBkg_GeoAccCorr_h2D = h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]];
0556
0557 if (
0558 h1D_Z_center < temp_RotatedBkg_GeoAccCorr_h2D -> GetYaxis() -> GetXmin() ||
0559 h1D_Z_center > temp_RotatedBkg_GeoAccCorr_h2D -> GetYaxis() -> GetXmax() ||
0560 h1D_Z_BinWidth != temp_RotatedBkg_GeoAccCorr_h2D -> GetYaxis() -> GetBinWidth(1)
0561 )
0562 {
0563 std::cout<<"Error: h1D_Z_center is outside the h2D_GeoAccCorr_map or the binwidths are not matched"<<std::endl;
0564 exit(1);
0565 }
0566
0567 GeoAcc_correction = temp_RotatedBkg_GeoAccCorr_h2D -> GetBinContent( temp_RotatedBkg_GeoAccCorr_h2D -> FindBin(h1D_eta_center, h1D_Z_center) );
0568 GeoAcc_correction = (GeoAcc_correction <= 0) ? 0. : 1.;
0569
0570 std::cout<<std::endl;
0571 std::cout<<"----- "<<pair.first<<", eta_bin: "<<l_eta_bin<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", rotated: "<<l_rotated<<", finebin: "<<l_finebin<<std::endl;
0572 std::cout<<"----- eta_center: "<<h1D_eta_center<<", Z_center: "<<h1D_Z_center<<", GeoAcc_correction: "<<GeoAcc_correction<<std::endl;
0573 std::cout<<"----- h1D_DeltaPhi Integral: "<<pair.second -> Integral()<<std::endl;
0574 pair.second -> Scale(GeoAcc_correction);
0575 std::cout<<"----- h1D_DeltaPhi Integral post : "<<pair.second -> Integral()<<std::endl;
0576
0577 }
0578 else if (ApplyGeoAccCorr && h2D_GeoAccCorr_map.size() == 0)
0579 {
0580 std::cout<<"Error: h2D_GeoAccCorr_map is empty"<<std::endl;
0581 exit(1);
0582 }
0583 else
0584 {
0585 GeoAcc_correction = 1;
0586 }
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599 if (isTypeA == l_typeA){
0600 if (
0601 SelectedMbin == 100 ||
0602 (SelectedMbin == Semi_inclusive_interval && l_Mbin <= Semi_inclusive_Mbin) ||
0603 SelectedMbin == l_Mbin
0604 )
0605 {
0606 if (l_rotated == 0){
0607 pair.second->SetFillColor(ROOT_color_code[hstack1D_DeltaPhi_map[Form("hstack1D_DeltaPhi_Eta%d", l_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0608
0609 hstack1D_DeltaPhi_map[Form("hstack1D_DeltaPhi_Eta%d", l_eta_bin)] -> Add(pair.second);
0610 }
0611 else if (l_rotated == 1){
0612 pair.second->SetFillColor(ROOT_color_code[hstack1D_DeltaPhi_map[Form("hstack1D_DeltaPhi_Eta%d_rotated", l_eta_bin)]->GetNhists() % ROOT_color_code.size()]);
0613
0614 hstack1D_DeltaPhi_map[Form("hstack1D_DeltaPhi_Eta%d_rotated", l_eta_bin)] -> Add(pair.second);
0615 }
0616 else {
0617 std::cout<<"wtf_Inclusive100: "<<pair.first<<std::endl;
0618 }
0619 }
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652
0653
0654 else
0655 {
0656 std::cout<<"wtf_inclusive_to_Mbin: "<<pair.first<<std::endl;
0657 }
0658 }
0659 }
0660
0661
0662 if (runnumber == -1 && pair.first.find("h1D") != std::string::npos && pair.first.find("_TrueEta") != std::string::npos && l_Mbin != -1){
0663 if (SelectedMbin == 100){
0664 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0665 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0666 }
0667
0668 else if (SelectedMbin == Semi_inclusive_interval){
0669 if (l_Mbin > Semi_inclusive_Mbin) {continue;}
0670
0671 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0672 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0673 }
0674
0675 else if (SelectedMbin == l_Mbin){
0676
0677 pair.second->SetFillColor(ROOT_color_code[hstack1D_TrueEta_map[Form("hstack1D_TrueEta")]->GetNhists() % ROOT_color_code.size()]);
0678 hstack1D_TrueEta_map[Form("hstack1D_TrueEta")] -> Add(pair.second);
0679 }
0680 else {
0681 std::cout<<"wtf_TruthEta_inclusive_to_Mbin: "<<pair.first<<std::endl;
0682 }
0683 }
0684
0685
0686 if (pair.first.find("h1D") != std::string::npos && pair.first.find("BestPair") != std::string::npos && isTypeA == l_typeA){
0687 if (
0688 SelectedMbin == 100 ||
0689 (SelectedMbin == Semi_inclusive_interval && l_Mbin <= Semi_inclusive_Mbin) ||
0690 SelectedMbin == l_Mbin
0691 ){
0692 if (pair.first.find("ClusPhiSize") != std::string::npos){
0693 pair.second->SetFillColor(ROOT_color_code[hstack1D_BestPair_ClusPhiSize->GetNhists() % ROOT_color_code.size()]);
0694 hstack1D_BestPair_ClusPhiSize -> Add(pair.second);
0695 }
0696 else if (pair.first.find("ClusAdc") != std::string::npos){
0697 pair.second->SetFillColor(ROOT_color_code[hstack1D_BestPair_ClusAdc->GetNhists() % ROOT_color_code.size()]);
0698 hstack1D_BestPair_ClusAdc -> Add(pair.second);
0699 }
0700 else if (pair.first.find("DeltaPhi") != std::string::npos){
0701 pair.second->SetFillColor(ROOT_color_code[hstack1D_BestPair_DeltaPhi->GetNhists() % ROOT_color_code.size()]);
0702 hstack1D_BestPair_DeltaPhi -> Add(pair.second);
0703 }
0704 else if (pair.first.find("DeltaEta") != std::string::npos){
0705 pair.second->SetFillColor(ROOT_color_code[hstack1D_BestPair_DeltaEta->GetNhists() % ROOT_color_code.size()]);
0706 hstack1D_BestPair_DeltaEta -> Add(pair.second);
0707 }
0708 }
0709 }
0710 }
0711
0712
0713
0714 for (auto &pair : hstack1D_DeltaPhi_map)
0715 {
0716 std::tie(l_eta_bin, l_vtxz_bin, l_Mbin, l_typeA, l_rotated, l_finebin) = GetHistStringInfo(pair.first);
0717
0718 if (pair.first.find("hstack1D") != std::string::npos && pair.first.find("_DeltaPhi") != std::string::npos && pair.first.find("_rotated") == std::string::npos && l_eta_bin != -1){
0719
0720 std::cout<<"----- "<<pair.first<<", eta_bin: "<<l_eta_bin<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", rotated: "<<l_rotated<<", finebin: "<<l_finebin<<std::endl;
0721
0722
0723 auto temp_hist = (TH1D*) pair.second -> GetStack() -> Last();
0724 auto temp_hist_rotate = (TH1D*) hstack1D_DeltaPhi_map[pair.first + "_rotated"] -> GetStack() -> Last();
0725
0726 double hist_offset_rotate = get_dist_offset(temp_hist_rotate, 15);
0727
0728 std::string f1_BkgPolTwo_Fit_map_key = pair.first + "_rotated" + "_BkgPolTwo_Fit";
0729 f1_BkgPolTwo_Fit_map[f1_BkgPolTwo_Fit_map_key] -> SetParameters(hist_offset_rotate, 0, 0);
0730 temp_hist_rotate -> Fit(f1_BkgPolTwo_Fit_map[f1_BkgPolTwo_Fit_map_key], "N");
0731
0732
0733
0734
0735
0736 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)] = (TH1D*) temp_hist -> Clone(Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin));
0737 TH1D * temp_signal_hist = h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)];
0738 for (int hist_i = 1; hist_i <= temp_signal_hist -> GetNbinsX(); hist_i++){
0739
0740 temp_signal_hist -> SetBinContent(hist_i, temp_hist -> GetBinContent(hist_i) - f1_BkgPolTwo_Fit_map[f1_BkgPolTwo_Fit_map_key] -> Eval(temp_hist -> GetBinCenter(hist_i)));
0741 }
0742
0743
0744
0745 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)] -> SetFillColorAlpha(1,0);
0746 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)] -> SetLineColor(8);
0747
0748 }
0749 }
0750
0751
0752
0753 for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0754
0755 if (isTypeA == 0 && SelectedMbin == 100){
0756 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0757 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0758 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0759 }
0760 else if (isTypeA == 0 && SelectedMbin == Semi_inclusive_interval){
0761 if (Mbin > Semi_inclusive_Mbin) {continue;}
0762
0763 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0764 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0765 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0766 }
0767 else if (isTypeA == 0 && SelectedMbin == Mbin){
0768 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0769 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0770 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0771 }
0772
0773
0774 else if (isTypeA == 1 && SelectedMbin == 100){
0775 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0776 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0777 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0778 }
0779 else if (isTypeA == 1 && SelectedMbin == Semi_inclusive_interval){
0780 if (Mbin > Semi_inclusive_Mbin) {continue;}
0781
0782 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0783 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0784 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0785 }
0786 else if (isTypeA == 1 && SelectedMbin == Mbin){
0787 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d",Mbin)]);
0788 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin",Mbin)]);
0789 hstack2D_GoodProtoTracklet_map["hstack2D_GoodProtoTracklet_EtaVtxZ_rotated"] -> Add(h2D_input_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated",Mbin)]);
0790 }
0791
0792 }
0793
0794
0795
0796
0797
0798
0799 if (runnumber == -1)
0800 {
0801 for (auto &pair : h2D_input_map)
0802 {
0803 std::tie(l_eta_bin, l_vtxz_bin, l_Mbin, l_typeA, l_rotated, l_finebin) = GetHistStringInfo(pair.first);
0804
0805 if (pair.first.find("h2D_TrueEtaVtxZ") != std::string::npos && l_Mbin != -1){
0806
0807 std::cout<<" Truth for TrueEtaVtxZ ----- "<<pair.first<<", eta_bin: "<<l_eta_bin<<", vtxz_bin: "<<l_vtxz_bin<<", Mbin: "<<l_Mbin<<", typeA: "<<l_typeA<<", rotated: "<<l_rotated<<", finebin: "<<l_finebin<<std::endl;
0808
0809 if (SelectedMbin == 100){
0810 if (l_finebin == 0){
0811 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"] -> Add(pair.second);
0812 }
0813 else if (l_finebin == 1){
0814 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"] -> Add(pair.second);
0815 }
0816 else {
0817 std::cout<<"wtf_TrueEtaVtxZ_inclusive100: "<<pair.first<<std::endl;
0818 }
0819 }
0820 else if (SelectedMbin == Semi_inclusive_interval){
0821 if (l_Mbin > Semi_inclusive_Mbin) {continue;}
0822
0823 if (l_finebin == 0){
0824 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"] -> Add(pair.second);
0825 }
0826 else if (l_finebin == 1){
0827 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"] -> Add(pair.second);
0828 }
0829 else {
0830 std::cout<<"wtf_TrueEtaVtxZ_inclusive"<<Semi_inclusive_interval<<" : "<<pair.first<<std::endl;
0831 }
0832 }
0833 else if (SelectedMbin == l_Mbin){
0834 if (l_finebin == 0){
0835 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"] -> Add(pair.second);
0836 }
0837 else if (l_finebin == 1){
0838 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"] -> Add(pair.second);
0839 }
0840 else {
0841 std::cout<<"wtf_TrueEtaVtxZ_singlebin: "<<pair.first<<std::endl;
0842 }
0843 }
0844 else {
0845 std::cout<<"wtf_TrueEtaVtxZ_inclusive_to_Mbin: "<<pair.first<<std::endl;
0846 }
0847
0848 }
0849
0850
0851 }
0852 }
0853
0854
0855
0856 for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0857 {
0858 std::string isTypeA_str = (isTypeA == 0) ? "" : "_typeA";
0859
0860 if (
0861 SelectedMbin == 100 ||
0862 (SelectedMbin == Semi_inclusive_interval && Mbin <= Semi_inclusive_Mbin) ||
0863 SelectedMbin == Mbin
0864 )
0865 {
0866 hstack2D_BestPairEtaVtxZ -> Add(h2D_input_map[Form("h2D%s_BestPairEtaVtxZ_Mbin%d",isTypeA_str.c_str(), Mbin)]);
0867 hstack2D_BestPairEtaVtxZ_FineBin -> Add(h2D_input_map[Form("h2D%s_BestPairEtaVtxZ_Mbin%d_FineBin",isTypeA_str.c_str(), Mbin)]);
0868 }
0869 }
0870
0871 }
0872
0873 void PreparedNdEtaEach::DoFittings()
0874 {
0875 std::cout<<"In DoFittings()"<<std::endl;
0876
0877 int stack_count = 0;
0878
0879 int l_eta_bin, l_vtxz_bin, l_Mbin;
0880 int l_typeA, l_rotated, l_finebin;
0881
0882 for (auto &pair : hstack1D_DeltaPhi_map)
0883 {
0884 std::tie(l_eta_bin, l_vtxz_bin, l_Mbin, l_typeA, l_rotated, l_finebin) = GetHistStringInfo(pair.first);
0885
0886 if (pair.first.find("hstack1D") != std::string::npos && pair.first.find("_DeltaPhi") != std::string::npos){
0887 if (stack_count % 20 == 0){
0888 std::cout<<"Fitting stack : "<<stack_count<<", "<<pair.first<<std::endl;
0889 }
0890
0891 auto temp_hist = (TH1D*) pair.second -> GetStack() -> Last();
0892
0893 std::vector<double> N_group_info = find_Ngroup(temp_hist);
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903 std::string f1_SigBkgPol2_Fit_map_key = pair.first + "_SigBkgPol2_Fit";
0904 std::string f1_SigBkgPol2_DrawSig_map_key = pair.first + "_SigBkgPol2_DrawSig";
0905 std::string f1_SigBkgPol2_DrawBkgPol2_map_key = pair.first + "_SigBkgPol2_DrawBkgPol2";
0906
0907 if (f1_SigBkgPol2_Fit_map.find(f1_SigBkgPol2_Fit_map_key) == f1_SigBkgPol2_Fit_map.end()){
0908 std::cout<<f1_SigBkgPol2_Fit_map_key<<" not found in f1_SigBkgPol2_Fit_map !!"<<std::endl;
0909 exit(1);
0910 }
0911
0912 double hist_offset = get_dist_offset(temp_hist, 15);
0913 double signal_region_l = cut_GoodProtoTracklet_DeltaPhi->first;
0914 double signal_region_r = cut_GoodProtoTracklet_DeltaPhi->second;
0915
0916
0917
0918
0919
0920
0921
0922
0923
0924
0925
0926
0927
0928
0929
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939
0940
0941
0942 if (pair.first.find("_rotated") == std::string::npos)
0943 {
0944
0945 double gaus_height = temp_hist->GetBinContent(temp_hist->GetMaximumBin()) - hist_offset;
0946 double gaus_width = fabs(N_group_info[3]-N_group_info[2]) / 2.;
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> SetParameters(
0958 gaus_height, 0, gaus_width,
0959 hist_offset, 0, 500000, 0
0960 );
0961
0962
0963
0964
0965 temp_hist -> Fit(f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key], "NQ");
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978 f1_SigBkgPol2_DrawSig_map[f1_SigBkgPol2_DrawSig_map_key] -> SetParameters(
0979 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(0),
0980 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(1),
0981 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(2)
0982 );
0983
0984 f1_SigBkgPol2_DrawBkgPol2_map[f1_SigBkgPol2_DrawBkgPol2_map_key] -> SetParameters(
0985 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(3),
0986 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(4),
0987 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(5),
0988 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> GetParameter(6)
0989 );
0990 }
0991
0992
0993
0994 file_out_DeltaPhi -> cd();
0995 c1 -> cd();
0996 if (pair.first.find("_rotated") == std::string::npos){
0997 auto temp_hist_rotate = (TH1D*) hstack1D_DeltaPhi_map[pair.first + "_rotated"] -> GetStack() -> Last();
0998
0999 std::string f1_BkgPolTwo_Fit_map_key = pair.first + "_rotated" + "_BkgPolTwo_Fit";
1000
1001 temp_hist -> SetFillColor(0);
1002 temp_hist -> SetLineColor(9);
1003 temp_hist -> SetMinimum(0);
1004 temp_hist -> SetMaximum(temp_hist -> GetBinContent(temp_hist -> GetMaximumBin()) * 1.5);
1005 temp_hist_rotate -> SetFillColor(0);
1006 temp_hist_rotate -> SetLineColor(46);
1007
1008 temp_hist -> Draw("hist");
1009 temp_hist_rotate -> Draw("ep same");
1010 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)] -> Draw("hist same");
1011 f1_BkgPolTwo_Fit_map[f1_BkgPolTwo_Fit_map_key] -> Draw("l same");
1012 c1 -> Write(Form("c1_%s", pair.first.c_str()));
1013
1014
1015
1016
1017
1018 c1 -> Clear();
1019
1020 temp_hist -> Draw("hist");
1021 f1_SigBkgPol2_Fit_map[f1_SigBkgPol2_Fit_map_key] -> Draw("l same");
1022 f1_SigBkgPol2_DrawSig_map[f1_SigBkgPol2_DrawSig_map_key] -> Draw("l same");
1023 f1_SigBkgPol2_DrawBkgPol2_map[f1_SigBkgPol2_DrawBkgPol2_map_key] -> Draw("l same");
1024 c1 -> Write(Form("c1_%s", f1_SigBkgPol2_Fit_map_key.c_str()));
1025 c1 -> Clear();
1026
1027
1028 temp_hist->SetMinimum(hist_offset * 0.9);
1029 temp_hist->SetMaximum(hist_offset * 1.1);
1030 temp_hist -> Draw("hist");
1031 temp_hist_rotate -> Draw("ep same");
1032 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", l_eta_bin)] -> Draw("hist same");
1033 f1_BkgPolTwo_Fit_map[f1_BkgPolTwo_Fit_map_key] -> Draw("l same");
1034 c1 -> Print(Form("%s/%s", output_directory.c_str(), output_filename_pdf.c_str()));
1035 c1 -> Clear();
1036 }
1037
1038
1039 }
1040
1041 else {
1042 std::cout<<"wtf_DoFittings: "<<pair.first<<std::endl;
1043 }
1044
1045 stack_count++;
1046 }
1047 }
1048
1049 void PreparedNdEtaEach::PrepareMultiplicity()
1050 {
1051 std::cout<<"In PrepareMultiplicity()"<<std::endl;
1052
1053 h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")] -> Sumw2(false);
1054
1055 for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++){
1056
1057
1058 double rotated_bkg_count = get_hstack2D_GoodProtoTracklet_count(hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated")], eta_bin);
1059
1060
1061
1062 auto temp_h2D = (TH2D*) hstack2D_BestPairEtaVtxZ -> GetStack() -> Last();
1063
1064 if (false )
1065 {
1066 SetH2DNonZeroBins(h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]], 1);
1067 }
1068
1069 if (false )
1070 {
1071 double GelAcc_BestPair_BinWidthY = h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetBinWidth(1);
1072 double temp_h2D_BinWidthY = temp_h2D -> GetYaxis() -> GetBinWidth(1);
1073 int Ngroup = temp_h2D_BinWidthY / GelAcc_BestPair_BinWidthY;
1074
1075 std::cout<<"GelAcc_BestPair_BinWidthY: "<<GelAcc_BestPair_BinWidthY<<", temp_h2D_BinWidthY: "<<temp_h2D_BinWidthY<<", Ngroup: "<<Ngroup<<std::endl;
1076
1077 h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> Rebin2D(1, Ngroup);
1078 }
1079
1080 if (false ){
1081
1082 if (
1083 temp_h2D -> GetNbinsX() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetNbinsX() ||
1084 temp_h2D -> GetXaxis() -> GetXmin() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetXaxis() -> GetXmin() ||
1085 temp_h2D -> GetXaxis() -> GetXmax() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetXaxis() -> GetXmax() ||
1086
1087 temp_h2D -> GetNbinsY() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetNbinsY() ||
1088 temp_h2D -> GetYaxis() -> GetXmin() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmin() ||
1089 temp_h2D -> GetYaxis() -> GetXmax() != h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmax()
1090 )
1091 {
1092 std::cout<<"h2D_BestPairEtaVtxZ and h2D_GeoAccCorr_map[For_BestPair] have different binning!!"<<std::endl;
1093 std::cout<<"temp_h2D X: "<<temp_h2D -> GetNbinsX()<<", "<<temp_h2D -> GetXaxis() -> GetXmin()<<", "<<temp_h2D -> GetXaxis() -> GetXmax()<<std::endl;
1094 std::cout<<"temp_h2D Y: "<<temp_h2D -> GetNbinsY()<<", "<<temp_h2D -> GetYaxis() -> GetXmin()<<", "<<temp_h2D -> GetYaxis() -> GetXmax()<<std::endl;
1095
1096 std::cout<<"h2D_GeoAccCorr_map[For_BestPair] X: "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetNbinsX()<<", "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetXaxis() -> GetXmin()<<", "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetXaxis() -> GetXmax()<<std::endl;
1097 std::cout<<"h2D_GeoAccCorr_map[For_BestPair] Y: "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetNbinsY()<<", "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmin()<<", "<<h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> GetYaxis() -> GetXmax()<<std::endl;
1098 exit(1);
1099 }
1100
1101 temp_h2D -> Multiply(h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]]);
1102 }
1103
1104 h1D_BestPair_RecoTrackletEta_map[Form("h1D_BestPair_RecoTrackletEta")] -> SetBinContent(
1105 eta_bin + 1,
1106
1107 get_h2D_GoodProtoTracklet_count(temp_h2D, eta_bin)
1108 );
1109
1110 std::pair<int, int> DeltaPhi_signal_bin = get_DeltaPhi_SingleBin(h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin)], BkgRotated_DeltaPhi_Signal_range);
1111
1112
1113 h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")] -> SetBinContent(
1114 eta_bin + 1,
1115
1116 h1D_RotatedBkg_DeltaPhi_Signal_map[Form("h1D_RotatedBkg_DeltaPhi_Signal_Eta%d", eta_bin)] -> Integral(DeltaPhi_signal_bin.first, DeltaPhi_signal_bin.second)
1117
1118
1119
1120 );
1121
1122 std::cout<<"h1D_RotatedBkg_RecoTrackletEta, bin: "<<eta_bin + 1
1123 <<", Eta: ["<<h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")]->GetXaxis()->GetBinLowEdge(eta_bin + 1)<<", "<<h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")]->GetXaxis()->GetBinUpEdge(eta_bin + 1)<<"], "
1124 <<"Content&Errors: {"<<h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")]->GetBinContent(eta_bin + 1)<<", "<<h1D_RotatedBkg_RecoTrackletEta_map[Form("h1D_RotatedBkg_RecoTrackletEta")]->GetBinError(eta_bin + 1)<<"}"<<std::endl;
1125 }
1126 }
1127
1128 void PreparedNdEtaEach::PreparedNdEtaHist()
1129 {
1130 std::cout<<"In PreparedNdEtaHist()"<<std::endl;
1131
1132 if (runnumber == -1){
1133 h1D_TruedNdEta_map.clear();
1134
1135 h1D_TruedNdEta_map.insert(
1136 std::make_pair(
1137 "h1D_TruedNdEta",
1138 (TH1D*)((TH1D*)hstack1D_TrueEta_map["hstack1D_TrueEta"] -> GetStack() -> Last()) -> Clone("h1D_TruedNdEta")
1139 )
1140 );
1141 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetName("h1D_TruedNdEta");
1142 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetTitle("h1D_TruedNdEta");
1143 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Scale(1. / get_EvtCount(h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"], SelectedMbin));
1144 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Scale(1. / double(h1D_TruedNdEta_map["h1D_TruedNdEta"]->GetBinWidth(1)));
1145
1146 }
1147
1148
1149
1150 h1D_BestPair_RecoTrackletEtaPerEvt_map.clear();
1151 h1D_BestPair_RecoTrackletEtaPerEvt_map.insert(
1152 std::make_pair(
1153 "h1D_BestPair_RecoTrackletEtaPerEvt",
1154 (TH1D*)h1D_BestPair_RecoTrackletEta_map["h1D_BestPair_RecoTrackletEta"] -> Clone("h1D_BestPair_RecoTrackletEtaPerEvt")
1155 )
1156 );
1157 h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> SetName("h1D_BestPair_RecoTrackletEtaPerEvt");
1158 h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> SetTitle("h1D_BestPair_RecoTrackletEtaPerEvt");
1159 h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> Scale(1. / get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], SelectedMbin));
1160 h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> Scale(1. / double(h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"]->GetBinWidth(1)));
1161
1162
1163
1164 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.clear();
1165 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map.insert(
1166 std::make_pair(
1167 "h1D_RotatedBkg_RecoTrackletEtaPerEvt",
1168 (TH1D*)h1D_RotatedBkg_RecoTrackletEta_map["h1D_RotatedBkg_RecoTrackletEta"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvt")
1169 )
1170 );
1171 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> SetName("h1D_RotatedBkg_RecoTrackletEtaPerEvt");
1172 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> SetTitle("h1D_RotatedBkg_RecoTrackletEtaPerEvt");
1173 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> Scale(1. / get_EvtCount(h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"], SelectedMbin));
1174 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> Scale(1. / double(h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"]->GetBinWidth(1)));
1175
1176
1177
1178
1179 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map.clear();
1180 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map.insert(
1181 std::make_pair(
1182 "h1D_BestPair_RecoTrackletEtaPerEvtPostAC",
1183 (TH1D*)h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> Clone("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")
1184 )
1185 );
1186 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map["h1D_BestPair_RecoTrackletEtaPerEvtPostAC"] -> SetName("h1D_BestPair_RecoTrackletEtaPerEvtPostAC");
1187 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map["h1D_BestPair_RecoTrackletEtaPerEvtPostAC"] -> SetTitle("h1D_BestPair_RecoTrackletEtaPerEvtPostAC");
1188
1189 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.clear();
1190 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map.insert(
1191 std::make_pair(
1192 "h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC",
1193 (TH1D*)h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> Clone("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")
1194 )
1195 );
1196 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC"] -> SetName("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC");
1197 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC"] -> SetTitle("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC");
1198
1199
1200
1201
1202
1203 if (h1D_alpha_correction_map_in.size() != 0 && ApplyAlphaCorr == true)
1204 {
1205
1206 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> Sumw2(true);
1207 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> Divide(h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")], h1D_alpha_correction_map_in[alpha_correction_name_map[0]], 1, 1);
1208
1209
1210 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")] -> Sumw2(true);
1211 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")] -> Divide(h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")], h1D_alpha_correction_map_in[alpha_correction_name_map[1]], 1, 1);
1212
1213
1214 for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
1215 {
1216 if (
1217 cut_EtaRange_pair.first == cut_EtaRange_pair.first &&
1218 (
1219 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> GetXaxis() -> GetBinUpEdge(eta_bin + 1) < cut_EtaRange_pair.first ||
1220 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> GetXaxis() -> GetBinLowEdge(eta_bin + 1) > cut_EtaRange_pair.second
1221 )
1222 ){
1223 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> SetBinContent(eta_bin + 1, 0);
1224 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")] -> SetBinError(eta_bin + 1, 0);
1225
1226 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")] -> SetBinContent(eta_bin + 1, 0);
1227 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")] -> SetBinError(eta_bin + 1, 0);
1228
1229 continue;
1230 }
1231 }
1232
1233 }
1234 else if (h1D_alpha_correction_map_in.size() == 0 && ApplyAlphaCorr == true){
1235 std::cout<<"wtf_PreparedNdEtaHist: h1D_alpha_correction_map_in.size() == 0"<<std::endl;
1236 exit(1);
1237 }
1238
1239 }
1240
1241 void PreparedNdEtaEach::DeriveAlphaCorrection()
1242 {
1243 if (runnumber != -1) {
1244 std::cout<<"wtf_DeriveAlphaCorrection: runnumber != -1"<<std::endl;
1245 return;
1246 }
1247
1248 std::cout<<"In DeriveAlphaCorrection()"<<std::endl;
1249
1250 h1D_alpha_correction_map_out.insert(
1251 std::make_pair(
1252 alpha_correction_name_map[0],
1253 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)
1254 )
1255 );
1256
1257 h1D_alpha_correction_map_out.insert(
1258 std::make_pair(
1259 alpha_correction_name_map[1],
1260 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)
1261 )
1262 );
1263
1264
1265 h1D_alpha_correction_map_out[alpha_correction_name_map[0]] -> Sumw2(true);
1266 h1D_alpha_correction_map_out[alpha_correction_name_map[1]] -> Sumw2(true);
1267
1268 h1D_alpha_correction_map_out[alpha_correction_name_map[0]] -> Divide(h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_BestPair_RecoTrackletEtaPerEvtPostAC")], h1D_TruedNdEta_map[Form("h1D_TruedNdEta")]);
1269 h1D_alpha_correction_map_out[alpha_correction_name_map[1]] -> Divide(h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map[Form("h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC")], h1D_TruedNdEta_map[Form("h1D_TruedNdEta")]);
1270 }
1271
1272 void PreparedNdEtaEach::EndRun()
1273 {
1274 c1 -> Print(Form("%s/%s)", output_directory.c_str(), output_filename_pdf.c_str()));
1275
1276 std::cout<<111<<std::endl;
1277
1278 file_out_DeltaPhi -> cd();
1279
1280 for (auto &pair : hstack1D_DeltaPhi_map)
1281 {
1282 pair.second -> Write();
1283 }
1284
1285 std::cout<<2222<<std::endl;
1286
1287 for (auto &pair : h1D_RotatedBkg_DeltaPhi_Signal_map){
1288 pair.second -> Write();
1289 }
1290
1291 std::cout<<3333<<std::endl;
1292
1293 file_out_dNdEta -> cd();
1294 ((TH2D*)hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ")]->GetStack()->Last()) -> Write("h2D_GoodProtoTracklet_EtaVtxZ");
1295 ((TH2D*)hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin")]->GetStack()->Last()) -> Write("h2D_GoodProtoTracklet_EtaVtxZ_FineBin");
1296 ((TH2D*)hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated")]->GetStack()->Last()) -> Write("h2D_GoodProtoTracklet_EtaVtxZ_rotated");
1297 h2D_input_map["h2D_RecoEvtCount_vtxZCentrality"] -> Write();
1298
1299 std::cout<<4444<<std::endl;
1300
1301 if (runnumber == -1){
1302 ((TH2D*)hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"]->GetStack()->Last()) -> Write("h2D_TrueEtaVtxZ");
1303 ((TH2D*)hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"]->GetStack()->Last()) -> Write("h2D_TrueEtaVtxZ_FineBin");
1304
1305 h2D_input_map["h2D_TrueEvtCount_vtxZCentrality"] -> Write();
1306 }
1307
1308 std::cout<<5555<<std::endl;
1309 std::cout<<"test: hstack1D_BestPair_ClusPhiSize size(): "<<hstack1D_BestPair_ClusPhiSize->GetNhists()<<std::endl;
1310 ((TH1D*) hstack1D_BestPair_ClusPhiSize -> GetStack() -> Last()) -> Write("h1D_BestPair_ClusPhiSize");
1311 std::cout<<55551<<std::endl;
1312 ((TH1D*) hstack1D_BestPair_ClusAdc -> GetStack() -> Last()) -> Write("h1D_BestPair_ClusAdc");
1313 std::cout<<55552<<std::endl;
1314 ((TH1D*) hstack1D_BestPair_DeltaPhi -> GetStack() -> Last()) -> Write("h1D_BestPair_DeltaPhi");
1315 std::cout<<55553<<std::endl;
1316 ((TH1D*) hstack1D_BestPair_DeltaEta -> GetStack() -> Last()) -> Write("h1D_BestPair_DeltaEta");
1317 std::cout<<55554<<std::endl;
1318 ((TH2D*) hstack2D_BestPairEtaVtxZ -> GetStack() -> Last()) -> Write("h2D_BestPairEtaVtxZ");
1319 std::cout<<55555<<std::endl;
1320 ((TH2D*) hstack2D_BestPairEtaVtxZ_FineBin -> GetStack() -> Last()) -> Write("h2D_BestPairEtaVtxZ_FineBin");
1321 std::cout<<55556<<std::endl;
1322
1323 std::cout<<6666<<std::endl;
1324
1325 h1D_BestPair_RecoTrackletEta_map["h1D_BestPair_RecoTrackletEta"] -> Write();
1326 h1D_BestPair_RecoTrackletEtaPerEvt_map["h1D_BestPair_RecoTrackletEtaPerEvt"] -> Write();
1327 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map["h1D_BestPair_RecoTrackletEtaPerEvtPostAC"] -> Write();
1328 h1D_RotatedBkg_RecoTrackletEta_map["h1D_RotatedBkg_RecoTrackletEta"] -> Write();
1329 h1D_RotatedBkg_RecoTrackletEtaPerEvt_map["h1D_RotatedBkg_RecoTrackletEtaPerEvt"] -> Write();
1330 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC"] -> Write();
1331
1332 std::cout<<7777<<std::endl;
1333
1334 if (ApplyGeoAccCorr && h2D_GeoAccCorr_map.size() != 0){
1335 h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> Write("GeoAccMap_For_BestPair");
1336 h2D_GeoAccCorr_map[GeoAccCorr_name_map[0]] -> Write("GeoAccMap_For_RotatedBkg");
1337 }
1338
1339
1340 if (runnumber == -1){
1341 ((TH1D*)hstack1D_TrueEta_map["hstack1D_TrueEta"]->GetStack()->Last()) -> Write("h1D_TrueEta");
1342 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Write();
1343 }
1344
1345 std::cout<<8888<<std::endl;
1346
1347 if (h1D_alpha_correction_map_out.size() != 0){
1348 h1D_alpha_correction_map_out[alpha_correction_name_map[0]] -> Write();
1349 h1D_alpha_correction_map_out[alpha_correction_name_map[1]] -> Write();
1350 }
1351
1352 std::cout<<9999<<std::endl;
1353
1354 hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ")] -> Write();
1355 hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_FineBin")] -> Write();
1356 hstack2D_GoodProtoTracklet_map[Form("hstack2D_GoodProtoTracklet_EtaVtxZ_rotated")] -> Write();
1357 if (runnumber == -1){
1358 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ"] -> Write();
1359 hstack2D_TrueEtaVtxZ_map["hstack2D_TrueEtaVtxZ_FineBin"] -> Write();
1360 hstack1D_TrueEta_map["hstack1D_TrueEta"] -> Write();
1361 }
1362
1363 std::cout<<1010101<<std::endl;
1364
1365 hstack1D_BestPair_ClusPhiSize -> Write();
1366 hstack1D_BestPair_ClusAdc -> Write();
1367 hstack1D_BestPair_DeltaPhi -> Write();
1368 hstack1D_BestPair_DeltaEta -> Write();
1369 hstack2D_BestPairEtaVtxZ -> Write();
1370 hstack2D_BestPairEtaVtxZ_FineBin -> Write();
1371 std::cout<<101011111<<std::endl;
1372
1373 if (runnumber == -1){
1374 c1 -> Clear();
1375 c1 -> cd();
1376
1377 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetFillColorAlpha(0,0);
1378 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> SetMinimum(0.);
1379 h1D_TruedNdEta_map["h1D_TruedNdEta"] -> Draw("hist");
1380 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map["h1D_BestPair_RecoTrackletEtaPerEvtPostAC"]->SetLineColor(2);
1381 h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map["h1D_BestPair_RecoTrackletEtaPerEvtPostAC"]->Draw("ep same");
1382 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC"]->SetLineColor(8);
1383 h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map["h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC"]->Draw("ep same");
1384
1385 c1 -> Write("c1_comp_h1D_dNdEta");
1386 c1 -> Clear();
1387 }
1388
1389 std::cout<<10102222<<std::endl;
1390
1391 file_out_dNdEta -> Close();
1392
1393 file_out_DeltaPhi -> Close();
1394 }
1395
1396 std::vector<double> PreparedNdEtaEach::find_Ngroup(TH1D * hist_in)
1397 {
1398 double Highest_bin_Content = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
1399 double Highest_bin_Center = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
1400
1401 int group_Nbin = 0;
1402 int peak_group_ID = -9999;
1403 double group_entry = 0;
1404 double peak_group_ratio;
1405 std::vector<int> group_Nbin_vec; group_Nbin_vec.clear();
1406 std::vector<double> group_entry_vec; group_entry_vec.clear();
1407 std::vector<double> group_widthL_vec; group_widthL_vec.clear();
1408 std::vector<double> group_widthR_vec; group_widthR_vec.clear();
1409
1410 for (int i = 0; i < hist_in -> GetNbinsX(); i++){
1411
1412 double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
1413
1414 if (bin_content != 0){
1415
1416 if (group_Nbin == 0) {
1417 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1418 }
1419
1420 group_Nbin += 1;
1421 group_entry += bin_content;
1422 }
1423 else if (bin_content == 0 && group_Nbin != 0){
1424 group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1425 group_Nbin_vec.push_back(group_Nbin);
1426 group_entry_vec.push_back(group_entry);
1427 group_Nbin = 0;
1428 group_entry = 0;
1429 }
1430 }
1431 if (group_Nbin != 0) {
1432 group_Nbin_vec.push_back(group_Nbin);
1433 group_entry_vec.push_back(group_entry);
1434 group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
1435 }
1436
1437
1438 for (int i = 0; i < int(group_Nbin_vec.size()); i++){
1439 if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
1440 peak_group_ID = i;
1441 break;
1442 }
1443 }
1444
1445
1446
1447
1448
1449 peak_group_ratio = -999;
1450
1451 double peak_group_left = (double(group_Nbin_vec.size()) == 0) ? -999 : group_widthL_vec[peak_group_ID];
1452 double peak_group_right = (double(group_Nbin_vec.size()) == 0) ? 999 : group_widthR_vec[peak_group_ID];
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469 return {double(group_Nbin_vec.size()), peak_group_ratio, peak_group_left, peak_group_right};
1470 }
1471
1472 double PreparedNdEtaEach::get_dist_offset(TH1D * hist_in, int check_N_bin)
1473 {
1474 if (check_N_bin < 0 || check_N_bin > hist_in -> GetNbinsX()) {std::cout<<" wrong check_N_bin "<<std::endl; exit(1);}
1475 double total_entry = 0;
1476 for (int i = 0; i < check_N_bin; i++)
1477 {
1478 total_entry += hist_in -> GetBinContent(i+1);
1479 total_entry += hist_in -> GetBinContent(hist_in -> GetNbinsX() - i);
1480 }
1481
1482 return total_entry/double(2. * check_N_bin);
1483 }
1484
1485 double PreparedNdEtaEach::get_hstack2D_GoodProtoTracklet_count(THStack * stack_in, int eta_bin_in)
1486 {
1487 auto temp_hist = (TH2D*) stack_in -> GetStack() -> Last();
1488
1489 double count = 0;
1490
1491
1492
1493
1494
1495 for (int yi = 1; yi <= temp_hist->GetNbinsY(); yi++){
1496
1497 if (vtxZ_index_map.find(yi - 1) == vtxZ_index_map.end()){
1498 continue;
1499 }
1500
1501 count += temp_hist -> GetBinContent(eta_bin_in + 1, yi);
1502 }
1503
1504 return count;
1505
1506 }
1507
1508 double PreparedNdEtaEach::get_h2D_GoodProtoTracklet_count(TH2D * h2D_in, int eta_bin_in)
1509 {
1510 double count = 0;
1511
1512
1513
1514
1515
1516 for (int yi = 1; yi <= h2D_in->GetNbinsY(); yi++){
1517
1518 if (vtxZ_index_map.find(yi - 1) == vtxZ_index_map.end()){
1519 continue;
1520 }
1521
1522 count += h2D_in -> GetBinContent(eta_bin_in + 1, yi);
1523 }
1524
1525 return count;
1526 }
1527
1528 double PreparedNdEtaEach::get_EvtCount(TH2D * hist_in, int centrality_bin_in)
1529 {
1530 double count = 0;
1531
1532
1533
1534
1535 if (centrality_bin_in == 100){
1536 for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){
1537
1538 if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1539 continue;
1540 }
1541
1542 for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){
1543 count += hist_in -> GetBinContent(xi, yi);
1544 }
1545 }
1546 }
1547
1548 else if (centrality_bin_in == Semi_inclusive_interval){
1549 for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){
1550
1551 if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1552 continue;
1553 }
1554
1555 for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){
1556 if (yi-1 <= Semi_inclusive_Mbin){
1557 count += hist_in -> GetBinContent(xi, yi);
1558 }
1559 }
1560 }
1561 }
1562 else {
1563 for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){
1564
1565 if (vtxZ_index_map.find(xi - 1) == vtxZ_index_map.end()){
1566 continue;
1567 }
1568
1569 count += hist_in -> GetBinContent(xi, centrality_bin_in + 1);
1570 }
1571 }
1572
1573 return count;
1574
1575 }
1576
1577 void PreparedNdEtaEach::Convert_to_PerEvt(TH1D * hist_in, double Nevent)
1578 {
1579 for (int i = 1; i <= hist_in->GetNbinsX(); i++){
1580 hist_in -> SetBinContent(i, hist_in -> GetBinContent(i) / Nevent);
1581 hist_in -> SetBinError(i, hist_in -> GetBinError(i) / Nevent);
1582 }
1583
1584 }
1585
1586 std::pair<int,int> PreparedNdEtaEach::get_DeltaPhi_SingleBin(TH1D * hist_in, std::pair<double, double> range_in)
1587 {
1588 int bin_min = hist_in -> FindBin(range_in.first);
1589
1590 int bin_max = hist_in -> FindBin(range_in.second);
1591 bin_max = (fabs(hist_in -> GetXaxis() -> GetBinLowEdge(bin_max) - range_in.second) < 1e-9) ? bin_max - 1 : bin_max;
1592
1593 std::cout<<std::endl;
1594 std::cout<<"Input range : "<<range_in.first<<", "<<range_in.second<<std::endl;
1595 std::cout<<"Final selected bins: "<<bin_min<<", "<<bin_max<<std::endl;
1596 std::cout<<"Output range : "<<hist_in -> GetXaxis() -> GetBinLowEdge(bin_min)<<", "<<hist_in -> GetXaxis() -> GetBinUpEdge(bin_max)<<std::endl;
1597
1598 return std::make_pair(bin_min, bin_max);
1599 }
1600
1601 void PreparedNdEtaEach::SetH2DNonZeroBins(TH2D* hist_in, double value_in)
1602 {
1603 for (int xi = 1; xi <= hist_in->GetNbinsX(); xi++){
1604 for (int yi = 1; yi <= hist_in->GetNbinsY(); yi++){
1605 if (hist_in -> GetBinContent(xi, yi) != 0){
1606 hist_in -> SetBinContent(xi, yi, value_in);
1607 }
1608 }
1609 }
1610 }
