NewCode2024/PreparedNdEta/PreparedNdEtaEach.h

0001 #ifndef PREPAREDNDETAEACH_H
0002 #define PREPAREDNDETAEACH_H
0003
0004 #include <iostream>
0005 #include <vector>
0006 #include <string>
0007 #include <numeric>
0008 #include <cctype> // For isdigit
0009
0010 #include <TFile.h>
0011 #include <TTree.h>
0012 #include <TH2D.h>
0013 #include <TH1D.h>
0014 #include <TMath.h>
0015 #include <TF1.h>
0016 #include <TProfile.h>
0017 #include <TCanvas.h>
0018 #include <TPad.h>
0019 #include <TGraphErrors.h>
0020 #include <TLatex.h>
0021 #include <THStack.h>
0022 #include <TCanvas.h> // note : for the combined case
0023 #include <TGraph.h>  // note : for the combined case
0024
0025 #include <TKey.h>
0026 #include <TRandom.h> // note : for the offset
0027 #include <TRandom3.h> // note : for the offset
0028
0029 #include <TColor.h>
0030
0031 #include <TObjArray.h>
0032
0033 #include "../Constants.h"
0034
0035
0036 class PreparedNdEtaEach{
0037     public:
0038         PreparedNdEtaEach(
0039             int process_id_in, // note : 1 or 2
0040             int runnumber_in, // note : still, (54280 or -1)
0041             std::string input_directory_in,
0042             std::string input_file_name_in,
0043             std::string output_directory_in,
0044             std::string output_file_name_suffix_in,
0045
0046             bool ApplyAlphaCorr_in,
0047             bool ApplyGeoAccCorr_in,
0048
0049             bool isTypeA_in,
0050             std::pair<double,double> cut_INTTvtxZ_in,
0051             int SelectedMbin_in // note : 0, 1, ---- 10, 70, 100
0052         );
0053
0054         std::vector<std::string> GetOutputFileName() {
0055             return {output_filename_DeltaPhi, output_filename_dNdEta, output_filename_pdf};
0056         }
0057
0058         std::vector<std::string> GetAlphaCorrectionNameMap() {return alpha_correction_name_map;}
0059
0060         void SetAlphaCorrectionH1DMap(std::map<std::string, TH1D*> map_in) {
0061             h1D_alpha_correction_map_in = map_in;
0062
0063             for (auto &pair : h1D_alpha_correction_map_in){
0064                 if (std::find(alpha_correction_name_map.begin(), alpha_correction_name_map.end(), pair.first) == alpha_correction_name_map.end()){
0065                     std::cout << "Error : the alpha correction name is not found in the map" << std::endl;
0066                     exit(1);
0067                 }
0068
0069                 // if (alpha_correction_name_map.find(pair.first) == alpha_correction_name_map.end()){
0070                 //     std::cout << "Error : the alpha correction name is not found in the map" << std::endl;
0071                 //     exit(1);
0072                 // }
0073             }
0074         }
0075
0076         std::vector<std::string> GetGeoAccCorrNameMap() {return GeoAccCorr_name_map;}
0077
0078         void SetGeoAccCorrH2DMap(std::map<std::string, TH2D*> map_in) {
0079             h2D_GeoAccCorr_map = map_in;
0080
0081             for (auto &pair : h2D_GeoAccCorr_map){
0082                 if (std::find(GeoAccCorr_name_map.begin(), GeoAccCorr_name_map.end(), pair.first) == GeoAccCorr_name_map.end()){
0083                     std::cout << "Error : the GeoAccCorr name is not found in the map" << std::endl;
0084                     exit(1);
0085                 }
0086             }
0087
0088         }
0089
0090
0091         void SetSelectedEtaRange(std::pair<double, double> cut_EtaRange_in) {
0092             cut_EtaRange_pair = cut_EtaRange_in;
0093         }
0094
0095         void SetBkgRotated_DeltaPhi_Signal_range(std::pair<double, double> range_in) {
0096             BkgRotated_DeltaPhi_Signal_range = range_in;
0097         }
0098
0099         void PrepareStacks();
0100         void DoFittings();
0101         void PrepareMultiplicity();
0102         void PreparedNdEtaHist();
0103         void DeriveAlphaCorrection();
0104         void EndRun();
0105
0106     protected:
0107
0108         // Division:---for the constructor------------------------------------------------------------------------------------------------
0109         int process_id;
0110         int runnumber;
0111         std::string input_directory;
0112         std::string input_file_name;
0113         std::string output_directory;
0114         std::string output_file_name_suffix;
0115         bool ApplyAlphaCorr;
0116         bool ApplyGeoAccCorr;
0117
0118         bool isTypeA;
0119         std::pair<double,double> cut_INTTvtxZ;
0120         int SelectedMbin;
0121
0122         // Division:---input root file------------------------------------------------------------------------------------------------
0123         TFile * file_in;
0124         TTree * tree_in;
0125         void PrepareInputRootFie();
0126         std::map<std::string, TH1D*> h1D_input_map;
0127         std::map<std::string, TH2D*> h2D_input_map;
0128
0129         std::vector<double> *centrality_edges;
0130         int nCentrality_bin;
0131
0132         double CentralityFineEdge_min;
0133         double CentralityFineEdge_max;
0134         int nCentralityFineBin;
0135
0136         double EtaEdge_min;
0137         double EtaEdge_max;
0138         int nEtaBin;
0139
0140         double VtxZEdge_min;
0141         double VtxZEdge_max;
0142         int nVtxZBin;
0143
0144         double DeltaPhiEdge_min;
0145         double DeltaPhiEdge_max;
0146         int nDeltaPhiBin;
0147
0148         std::pair<double, double> *cut_GoodProtoTracklet_DeltaPhi;
0149
0150
0151         // Division:---output root file------------------------------------------------------------------------------------------------
0152         TFile * file_out_DeltaPhi;
0153         TFile * file_out_dNdEta;
0154         std::string output_filename;
0155         std::string output_filename_DeltaPhi;
0156         std::string output_filename_dNdEta;
0157         std::string output_filename_pdf;
0158         void PrepareOutPutFileName();
0159         void PrepareOutPutRootFile();
0160
0161         TCanvas * c1;
0162
0163         // Division:---analysis------------------------------------------------------------------------------------------------
0164         std::map<int, int> GetVtxZIndexMap();
0165         void PrepareHistFits();
0166         std::string GetObjectName();
0167         std::tuple<int, int, int, int, int, int> GetHistStringInfo(std::string hist_name);
0168         std::vector<double> find_Ngroup(TH1D * hist_in);
0169         double get_dist_offset(TH1D * hist_in, int check_N_bin);
0170         double get_hstack2D_GoodProtoTracklet_count(THStack * stack_in, int eta_bin_in);
0171         double get_h2D_GoodProtoTracklet_count(TH2D * h2D_in, int eta_bin_in);
0172         double get_EvtCount(TH2D * hist_in, int centrality_bin_in);
0173         void Convert_to_PerEvt(TH1D * hist_in, double Nevent);
0174         void SetH2DNonZeroBins(TH2D* hist_in, double value_in);
0175         std::pair<int,int> get_DeltaPhi_SingleBin(TH1D * hist_in, std::pair<double, double> range_in);
0176
0177         std::map<int, int> vtxZ_index_map;
0178         std::pair<double, double> cut_EtaRange_pair = {std::nan(""), std::nan("")};
0179
0180         // Division:---histogram------------------------------------------------------------------------------------------------
0181         // note : signal (eta-vtxZ) centrality segemntation
0182         std::map<std::string, TH1D*> h1D_BestPair_RecoTrackletEta_map;    // note : for the total count
0183         std::map<std::string, TH1D*> h1D_BestPair_RecoTrackletEtaPerEvt_map;
0184         std::map<std::string, TH1D*> h1D_BestPair_RecoTrackletEtaPerEvtPostAC_map;
0185
0186         std::map<std::string, THStack*> hstack1D_DeltaPhi_map;
0187         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEta_map; // note : for the total count
0188         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEtaPerEvt_map;
0189         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map;
0190
0191         std::map<std::string, TH1D*> h1D_RotatedBkg_DeltaPhi_Signal_map;
0192
0193         std::map<std::string, THStack*> hstack2D_GoodProtoTracklet_map; // note : for the good pair
0194
0195         // note : for the best pair
0196         THStack * hstack1D_BestPair_ClusPhiSize;
0197         THStack * hstack1D_BestPair_ClusAdc;
0198         THStack * hstack1D_BestPair_DeltaPhi;
0199         THStack * hstack1D_BestPair_DeltaEta;
0200         THStack * hstack2D_BestPairEtaVtxZ;
0201         THStack * hstack2D_BestPairEtaVtxZ_FineBin;
0202
0203         // note : for the final Truth
0204         std::map<std::string, THStack*> hstack2D_TrueEtaVtxZ_map;
0205         std::map<std::string, THStack*> hstack1D_TrueEta_map; // note : for the total count
0206         std::map<std::string, TH1D*> h1D_TruedNdEta_map;
0207
0208         // note : for the alpha correction
0209         std::map<std::string, TH1D*> h1D_alpha_correction_map_in;
0210         std::map<std::string, TH1D*> h1D_alpha_correction_map_out;
0211         std::vector<std::string> alpha_correction_name_map = {
0212             "h1D_BestPair_alpha_correction",
0213             "h1D_RotatedBkg_alpha_correction"
0214         };
0215
0216         // note : for Geo Acceprance correction
0217         std::map<std::string, TH2D*> h2D_GeoAccCorr_map;
0218         std::vector<std::string> GeoAccCorr_name_map;
0219
0220         // Division:---fitting------------------------------------------------------------------------------------------------
0221         // std::map<std::string, TF1*> f1_BkgPol2_Fit_map;
0222         // std::map<std::string, TF1*> f1_BkgPol2_Draw_map;
0223         std::map<std::string, TF1*> f1_SigBkgPol2_Fit_map;
0224         std::map<std::string, TF1*> f1_SigBkgPol2_DrawSig_map;
0225         std::map<std::string, TF1*> f1_SigBkgPol2_DrawBkgPol2_map;
0226
0227         std::map<std::string, TF1*> f1_BkgPolTwo_Fit_map;
0228
0229         // Division:---Constants------------------------------------------------------------------------------------------------
0230         int Semi_inclusive_Mbin = Constants::Semi_inclusive_bin;
0231         int Semi_inclusive_interval = Constants::Semi_inclusive_interval;
0232
0233         const std::vector<int> ROOT_color_code = {
0234             51, 61, 70, 79, 88, 98
0235         };
0236
0237         std::pair<double,double> BkgRotated_DeltaPhi_Signal_range = {-0.026, 0.026};
0238
0239
0240         // Division:---the functions------------------------------------------------------------------------------------------------
0241         static double gaus_func(double *x, double *par)
0242         {
0243             // note : par[0] : size
0244             // note : par[1] : mean
0245             // note : par[2] : width
0246             // note : par[3] : offset
0247             return par[0] * TMath::Gaus(x[0],par[1],par[2]); //+ par[3];
0248         }
0249
0250         static  double gaus_pol2_func(double *x, double *par)
0251         {
0252             // note : par[0] : size
0253             // note : par[1] : mean
0254             // note : par[2] : width
0255
0256             double gaus_func = par[0] * TMath::Gaus(x[0],par[1],par[2]);
0257             double pol2_func = par[3] + par[4]* (x[0]-par[6]) - fabs(par[5]) * pow((x[0]-par[6]),2);
0258
0259             return gaus_func + pol2_func;
0260         }
0261
0262         static  double bkg_pol2_func(double *x, double *par)
0263         {
0264             if (x[0] > par[4] && x[0] < par[5]) {
0265                 TF1::RejectPoint();
0266                 return 0;
0267             }
0268             return par[0] + par[1]* (x[0]-par[3]) - fabs(par[2]) * pow((x[0]-par[3]),2);
0269
0270             // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2, p[4] sets the signal range that should be excluded in the fit
0271         }
0272
0273         static  double full_pol2_func(double *x, double *par)
0274         {
0275             return par[0] + par[1]* (x[0]-par[3]) - fabs(par[2]) * pow((x[0]-par[3]),2);
0276
0277             // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2
0278         }
0279 };
0280
0281 #endif