ColumnCheck/GetMultiplicityMap/GetMultiplicityMap.cpp

0001 #include "GetMultiplicityMap.h"
0002
0003 GetMultiplicityMap::GetMultiplicityMap(
0004     int runnumber_in,
0005     std::string data_directory_in,
0006     std::string data_file_name_in,
0007     std::string MC_directory_in,
0008     std::string MC_file_name_in,
0009     std::string output_directory_in,
0010
0011     std::string output_file_name_suffix_in,
0012
0013     double SetMbinFloat_in,
0014     std::pair<double, double> VtxZRange_in,
0015     bool IsZClustering_in,
0016     bool BcoFullDiffCut_in,
0017     std::pair<bool, std::pair<double, double>> isClusQA_in
0018 ):
0019     runnumber(runnumber_in),
0020     data_directory(data_directory_in),
0021     data_file_name(data_file_name_in),
0022     MC_directory(MC_directory_in),
0023     MC_file_name(MC_file_name_in),
0024     output_directory(output_directory_in),
0025
0026     output_file_name_suffix(output_file_name_suffix_in),
0027
0028     SetMbinFloat(SetMbinFloat_in),
0029     VtxZRange(VtxZRange_in),
0030     IsZClustering(IsZClustering_in),
0031     BcoFullDiffCut(BcoFullDiffCut_in),
0032     isClusQA(isClusQA_in)
0033
0034 {
0035     std::cout<<1111<<std::endl;
0036     h1D_target_map.clear();
0037     h2D_target_map.clear();
0038     std::cout<<222<<std::endl;
0039     h1D_target_map = {
0040         {"h1D_ClusZ", {nullptr,nullptr}}
0041     };
0042
0043     for (int layer_i = 3; layer_i < 7; layer_i++)
0044     {
0045         int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0046
0047         for (int phi_i = 0; phi_i < phi_max; phi_i++)
0048         {
0049             h2D_target_map[Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",layer_i,phi_i)] = {nullptr,nullptr}; // note : X: Clus radius, Y: ClusZID
0050         }
0051     }
0052
0053     for (int i = 0; i < nZbin; i++)
0054     {
0055         h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)] = {nullptr,nullptr}; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0056     }
0057
0058     std::cout<<555<<std::endl;
0059
0060     PrepareInputRootFile();
0061
0062     std::cout<<666<<std::endl;
0063
0064     PrepareOutPutFileName();
0065
0066     std::cout<<777<<std::endl;
0067     PrepareOutPutRootFile();
0068
0069     std::cout<<888<<std::endl;
0070 }
0071
0072 void GetMultiplicityMap::PrepareInputRootFile()
0073 {
0074     file_in_data = new TFile(Form("%s/%s", data_directory.c_str(), data_file_name.c_str()), "READ");
0075     file_in_MC = new TFile(Form("%s/%s", MC_directory.c_str(), MC_file_name.c_str()), "READ");
0076
0077     if (file_in_data -> IsZombie() || file_in_MC -> IsZombie() || !file_in_data || !file_in_MC)
0078     {
0079         std::cerr<<"Error : Cannot open the input file"<<std::endl;
0080         exit(1);
0081     }
0082
0083     TFile * temp;
0084
0085     for (int i = 0; i < 2; i++){ // note : data -> MC
0086         temp = (i == 0) ? file_in_data : file_in_MC;
0087
0088         for (TObject* keyAsObj : *temp->GetListOfKeys())
0089         {
0090             auto key = dynamic_cast<TKey*>(keyAsObj);
0091             std::string hist_name  = key->GetName();
0092             std::string class_name = key->GetClassName();
0093
0094             if (class_name == "TH1D" && h1D_target_map.find(hist_name) != h1D_target_map.end())
0095             {
0096                 if (i == 0) {h1D_target_map[hist_name].first = (TH1D*) temp -> Get( hist_name.c_str() );}
0097                 else {h1D_target_map[hist_name].second = (TH1D*) temp -> Get( hist_name.c_str() );}
0098             }
0099
0100             if (class_name == "TH2D" && h2D_target_map.find(hist_name) != h2D_target_map.end())
0101             {
0102                 if (i == 0) {h2D_target_map[hist_name].first = (TH2D*) temp -> Get( hist_name.c_str() );}
0103                 else {h2D_target_map[hist_name].second = (TH2D*) temp -> Get( hist_name.c_str() );}
0104             }
0105         }
0106
0107     }
0108 }
0109
0110 void GetMultiplicityMap::PrepareOutPutFileName()
0111 {
0112     std::string runnumber_str = std::to_string( runnumber );
0113     if (runnumber != -1){
0114         int runnumber_str_len = 8;
0115         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0116     }
0117
0118     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0119         output_file_name_suffix = "_" + output_file_name_suffix;
0120     }
0121
0122     output_filename = "MulMap";
0123
0124     output_filename += (BcoFullDiffCut) ? "_BcoFullDiffCut" : "";
0125     output_filename += (IsZClustering) ? "_ZClustering" : "";
0126     output_filename += Form("_Mbin%.0f",SetMbinFloat);
0127     output_filename += Form("_VtxZ%.0fto%.0fcm",VtxZRange.first,VtxZRange.second);
0128     output_filename += (isClusQA.first) ? Form("_ClusQAAdc%.0fPhiSize%.0f",isClusQA.second.first,isClusQA.second.second) : "";
0129
0130     output_filename += output_file_name_suffix;
0131     output_filename += Form("_%s.root",runnumber_str.c_str());
0132 }
0133
0134 void GetMultiplicityMap::PrepareOutPutRootFile()
0135 {
0136     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0137 }
0138
0139 void GetMultiplicityMap::GetUsedZIDVec()
0140 {
0141     used_zid_data_vec.clear();
0142     used_zid_MC_vec.clear();
0143
0144     // note : data
0145     double HalfMaxEntries_data = h1D_target_map["h1D_ClusZ"].first -> GetBinContent( h1D_target_map["h1D_ClusZ"].first -> GetMaximumBin() ) / 2;
0146     for (int i = 1; i <= h1D_target_map["h1D_ClusZ"].first -> GetNbinsX(); i++)
0147     {
0148         if (h1D_target_map["h1D_ClusZ"].first -> GetBinContent(i) > HalfMaxEntries_data)
0149         {
0150             used_zid_data_vec.push_back(i-1); // note : start from zero
0151         }
0152     }
0153
0154     // note : MC
0155     double HalfMaxEntries_MC = h1D_target_map["h1D_ClusZ"].second -> GetBinContent( h1D_target_map["h1D_ClusZ"].second -> GetMaximumBin() ) / 2;
0156     for (int i = 1; i <= h1D_target_map["h1D_ClusZ"].second -> GetNbinsX(); i++)
0157     {
0158         if (h1D_target_map["h1D_ClusZ"].second -> GetBinContent(i) > HalfMaxEntries_MC)
0159         {
0160             used_zid_MC_vec.push_back(i-1); // note : start from zero
0161         }
0162     }
0163
0164     if (used_zid_data_vec.size() != used_zid_MC_vec.size())
0165     {
0166         std::cerr<<"Error : The number of ZID is different between data and MC"<<std::endl;
0167         exit(1);
0168     }
0169
0170     std::cout<<"used_zid_data_vec.size(): "<<used_zid_data_vec.size()<<std::endl;
0171
0172     for (int i = 0; i < used_zid_data_vec.size(); i++)
0173     {
0174         std::cout<<"Data: ZID : "<<used_zid_data_vec[i]<<". MC: "<<used_zid_MC_vec[i]<<std::endl;
0175
0176         if (used_zid_data_vec[i] != used_zid_MC_vec[i])
0177         {
0178             std::cerr<<"Error : The ZID is different between data and MC"<<std::endl;
0179             exit(1);
0180         }
0181     }
0182
0183     for (int i = 0; i < nZbin; i++)
0184     {
0185         if (std::find(used_zid_data_vec.begin(), used_zid_data_vec.end(), i) == used_zid_data_vec.end())
0186         {
0187             std::cout<<"ZID: "<<i<<" is dropped. Data Entries: "
0188             <<h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].first -> GetEntries()<<" MC Entries: "
0189             <<h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].second -> GetEntries()<<std::endl;
0190
0191             h2D_target_map.erase(Form("h2D_ClusCountLayerPhiId_ZId%d",i));
0192         }
0193     }
0194
0195 }
0196
0197 void GetMultiplicityMap::h2DNormalizedByRadius()
0198 {
0199     for (int i : used_zid_data_vec)
0200     {
0201         TH2D * temp_h2D_data = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].first; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0202
0203         for (int layer_i = 3; layer_i < 7; layer_i++)
0204         {
0205             int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0206
0207             for (int phi_i = 0; phi_i < phi_max; phi_i++)
0208             {
0209                 double radius = h2D_target_map[Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",layer_i,phi_i)].first -> GetMean(1); // note : X axis
0210                 double cover_angle = 2 * atan2((sensor_width/2), radius) * (180./M_PI); // note : in radian
0211
0212                 std::cout<<"Data: Layer: "<<layer_i<<", PhiId: "<<phi_i<<", Radius: "<<radius<<", CoverAngle: "<<cover_angle<<std::endl;
0213
0214                 int correspond_bin = temp_h2D_data -> FindBin(layer_i, phi_i);
0215                 double bin_content = temp_h2D_data -> GetBinContent(correspond_bin);
0216
0217                 temp_h2D_data -> SetBinContent(correspond_bin, bin_content / cover_angle);
0218                 temp_h2D_data -> SetBinError(correspond_bin, temp_h2D_data -> GetBinError(correspond_bin) / cover_angle);
0219             }
0220         }
0221
0222         // Division: ------------------------------------------------------------------------------------------------------------------------------------------------
0223         std::cout<<std::endl;
0224
0225         TH2D * temp_h2D_MC   = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].second; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0226         for (int layer_i = 3; layer_i < 7; layer_i++)
0227         {
0228             int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0229
0230             for (int phi_i = 0; phi_i < phi_max; phi_i++)
0231             {
0232                 double radius = h2D_target_map[Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",layer_i,phi_i)].second -> GetMean(1); // note : X axis
0233                 double cover_angle = 2 * atan2((sensor_width/2), radius) * (180./M_PI); // note : in radian
0234
0235                 std::cout<<"MC: Layer: "<<layer_i<<", PhiId: "<<phi_i<<", Radius: "<<radius<<", CoverAngle: "<<cover_angle<<std::endl;
0236
0237                 int correspond_bin = temp_h2D_MC -> FindBin(layer_i, phi_i);
0238                 double bin_content = temp_h2D_MC -> GetBinContent(correspond_bin);
0239
0240                 temp_h2D_MC -> SetBinContent(correspond_bin, bin_content / cover_angle);
0241                 temp_h2D_MC -> SetBinError(correspond_bin, temp_h2D_MC -> GetBinError(correspond_bin) / cover_angle);
0242             }
0243         }
0244     }
0245 }
0246
0247 void GetMultiplicityMap::h2DNormalized()
0248 {
0249     for (int i : used_zid_data_vec)
0250     {
0251         TH2D * temp_h2D_data = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].first; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0252         TH2D * temp_h2D_MC = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].second; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0253
0254         temp_h2D_data -> Sumw2(true);
0255         temp_h2D_MC -> Sumw2(true);
0256
0257         temp_h2D_data -> Scale(1. / temp_h2D_data -> GetBinContent(temp_h2D_data -> GetMaximumBin()));
0258         temp_h2D_MC -> Scale(1. / temp_h2D_MC -> GetBinContent(temp_h2D_MC -> GetMaximumBin()));
0259     }
0260 }
0261
0262 void GetMultiplicityMap::DataMCDivision()
0263 {
0264     h2D_map.clear();
0265     h1D_map.clear();
0266
0267     h1D_map["h1D_Ratio_All"] = new TH1D("h1D_Ratio_All","h1D_Ratio_All;Multiplicity Ratio (data/MC);Entries", 100, 0, 2);
0268
0269     for (int i : used_zid_data_vec)
0270     {
0271         TH2D * temp_h2D_data = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].first; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0272         TH2D * temp_h2D_MC = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].second; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0273
0274         h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] = (TH2D*) temp_h2D_data -> Clone(Form("h2D_RatioLayerPhiId_ZId%d",i));
0275         h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> SetTitle(Form("h2D_RatioLayerPhiId_ZId%d",i));
0276         h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> Divide(temp_h2D_data, temp_h2D_MC, 1, 1);
0277
0278         h1D_map[Form("h1D_Ratio_ZId%d",i)] = new TH1D(Form("h1D_Ratio_ZId%d",i),Form("h1D_Ratio_ZId%d;Multiplicity Ratio (data/MC);Entries",i), 100, 0, 2);
0279
0280         for (int layer_i = 3; layer_i < 7; layer_i++)
0281         {
0282             int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0283
0284             for (int phi_i = 0; phi_i < phi_max; phi_i++)
0285             {
0286                 int correspond_bin = h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> FindBin(layer_i, phi_i);
0287                 double bin_content = h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> GetBinContent(correspond_bin);
0288
0289                 if (bin_content == 0)
0290                 {
0291                     std::cout<<"Warning: bin_content == 0, ZID: "<<i<<", Layer: "<<layer_i<<", PhiId: "<<phi_i<<", data multiplicity: "<<temp_h2D_data->GetBinContent(correspond_bin)<<", MC multiplicity: "<<temp_h2D_MC->GetBinContent(correspond_bin)
0292                     <<std::endl;
0293                 }
0294
0295                 h1D_map[Form("h1D_Ratio_ZId%d",i)] -> Fill(bin_content);
0296                 h1D_map["h1D_Ratio_All"] -> Fill(bin_content);
0297             }
0298         }
0299     }
0300 }
0301
0302 void GetMultiplicityMap::PrepareMulMap()
0303 {
0304     h2D_MulMap = new TH2D("h2D_MulMap","h2D_MulMap;(LayerID-3) * 20 + PhiId;ZID", 80, 0, 80, nZbin, Zmin, Zmax);
0305     h2D_MaskedMap = new TH2D("h2D_MaskedMap","h2D_MaskedMap;(LayerID-3) * 20 + PhiId;ZID", 80, 0, 80, nZbin, Zmin, Zmax);
0306     h2D_RatioMap = new TH2D("h2D_RatioMap","h2D_RatioMap;(LayerID-3) * 20 + PhiId;ZID", 80, 0, 80, nZbin, Zmin, Zmax);
0307
0308
0309     for (int i : used_zid_data_vec) // note : only good ZID are here
0310     {
0311
0312         for (int layer_i = 3; layer_i < 7; layer_i++)
0313         {
0314             int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0315
0316             for (int phi_i = 0; phi_i < phi_max; phi_i++)
0317             {
0318                 int correspond_bin = h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> FindBin(layer_i, phi_i);
0319                 double column_ratio = h2D_map[Form("h2D_RatioLayerPhiId_ZId%d",i)] -> GetBinContent(correspond_bin);
0320                 column_ratio = (column_ratio != column_ratio) ? 0 : column_ratio;
0321
0322                 TH2D * temp_h2D_data = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].first; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0323                 TH2D * temp_h2D_MC = h2D_target_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)].second; // note : X: Layer {3 - 7}, Y: PhiId {0 -16}
0324
0325                 int index_x = (layer_i - 3) * 20 + phi_i + 1; // note : bin index, start at 1
0326                 int index_y = i + 1; // note : bin index, start at 1
0327
0328                 // std::cout<<"ZID: "<<i<<", Layer: "<<layer_i<<", PhiId: "<<phi_i<<", index_x: "<<index_x<<", index_y: "<<index_y<<", Ratio: "<<column_ratio<<std::endl;
0329
0330                 h2D_RatioMap -> SetBinContent(index_x, index_y, column_ratio);
0331
0332                 if (column_ratio >= Ratio_cut_pair.first && column_ratio <= Ratio_cut_pair.second)
0333                 {
0334                     h2D_MulMap -> SetBinContent(index_x, index_y, 1);
0335                 }
0336                 else
0337                 {
0338                     if (temp_h2D_data->GetBinContent(correspond_bin) > 0.000001 || temp_h2D_MC->GetBinContent(correspond_bin) > 0.000001)
0339                     {
0340                         std::cout<<"Masked: ZID: "<<i<<", Layer: "<<layer_i<<", PhiId: "<<phi_i<<", index_x: "<<index_x<<", index_y: "<<index_y<<", Ratio: "<<column_ratio
0341                         <<", Data: "<<temp_h2D_data->GetBinContent(correspond_bin)<<", MC: "<<temp_h2D_MC->GetBinContent(correspond_bin)
0342                         <<std::endl;
0343                         h2D_MaskedMap -> SetBinContent(index_x, index_y, 1);
0344                     }
0345
0346                 }
0347             }
0348         }
0349     }
0350 }
0351
0352 void GetMultiplicityMap::EndRun()
0353 {
0354     file_out -> cd();
0355
0356     h2D_RatioMap -> Write();
0357     h2D_MulMap -> Write();
0358     h2D_MaskedMap -> Write();
0359
0360
0361     for (auto pair : h2D_target_map){
0362         pair.second.first -> Write(Form("%s_data",pair.first.c_str()));
0363         pair.second.second -> Write(Form("%s_MC",pair.first.c_str()));
0364     }
0365     for (auto pair : h2D_map){
0366         pair.second -> Write();
0367     }
0368
0369     for (auto pair : h1D_map){
0370         pair.second -> Write();
0371     }
0372
0373     for (auto pair : h1D_target_map){
0374         pair.second.first -> Write(Form("%s_data",pair.first.c_str()));
0375         pair.second.second -> Write(Form("%s_MC",pair.first.c_str()));
0376     }
0377
0378
0379     file_out -> Close();
0380 }