PreparedNdEta/old_codes/PreparedNdEta.h

0001 #ifndef PREPAREDNDETA_H
0002 #define PREPAREDNDETA_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
0034 // #include "sPhenixStyle.h"
0035
0036 class PreparedNdEta{
0037     public:
0038         PreparedNdEta(
0039             int process_id_in,
0040             int runnumber_in,
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         );
0048
0049         std::map<std::string, int> GetAlphaCorrectionNameMap() {return alpha_correction_name_map;}
0050
0051         void SetAlphaCorrectionH2DMap(std::map<std::string, TH2D*> map_in) {
0052             h2D_alpha_correction_map_in = map_in;
0053
0054             for (auto &pair : h2D_alpha_correction_map_in){
0055                 if (alpha_correction_name_map.find(pair.first) == alpha_correction_name_map.end()){
0056                     std::cout << "Error : the alpha correction name is not found in the map" << std::endl;
0057                     exit(1);
0058                 }
0059             }
0060         }
0061         std::string GetOutputFileName() {return output_filename;}
0062
0063         void PrepareStacks();
0064         void DoFittings();
0065         void PrepareMultiplicity();
0066         void PreparedNdEtaHist();
0067         void DeriveAlphaCorrection();
0068         void EndRun();
0069
0070     protected:
0071
0072         // note : ----------------- for the constructor -----------------
0073         int process_id;
0074         int runnumber;
0075         std::string input_directory;
0076         std::string input_file_name;
0077         std::string output_directory;
0078         std::string output_file_name_suffix;
0079         bool ApplyAlphaCorr;
0080
0081         // note : ----------------- for the input file -----------------
0082         TFile * file_in;
0083         TTree * tree_in;
0084         void PrepareInputRootFie();
0085         std::map<std::string, TH1D*> h1D_input_map;
0086         std::map<std::string, TH2D*> h2D_input_map;
0087
0088         std::vector<double> *centrality_edges;
0089         int nCentrality_bin;
0090
0091         double CentralityFineEdge_min;
0092         double CentralityFineEdge_max;
0093         int nCentralityFineBin;
0094
0095         double EtaEdge_min;
0096         double EtaEdge_max;
0097         int nEtaBin;
0098
0099         double VtxZEdge_min;
0100         double VtxZEdge_max;
0101         int nVtxZBin;
0102
0103         double DeltaPhiEdge_min;
0104         double DeltaPhiEdge_max;
0105         int nDeltaPhiBin;
0106
0107         std::pair<double, double> *cut_GoodProtoTracklet_DeltaPhi;
0108
0109         // note : ----------------- for the histogram -----------------
0110         void PrepareHistFits();
0111         void PrepareTrackletEtaHist(std::map<std::string, TH1D*> &map_in, std::string name_in);
0112
0113         // note : multiplicity, dNdeta
0114         std::map<std::string, TH1D*> h1D_FitBkg_RecoTrackletEta_map;
0115         std::map<std::string, TH1D*> h1D_FitBkg_RecoTrackletEtaPerEvt_map;
0116         std::map<std::string, TH1D*> h1D_FitBkg_RecoTrackletEtaPerEvtPostAC_map;
0117         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEta_map;
0118         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEtaPerEvt_map;
0119         std::map<std::string, TH1D*> h1D_RotatedBkg_RecoTrackletEtaPerEvtPostAC_map;
0120
0121         std::map<std::string, THStack*> hstack1D_map;
0122
0123         // note : signal (eta-vtxZ) centrality segemntation
0124         std::map<std::string, TH2D*> h2D_map;
0125         std::map<std::string, THStack*> hstack2D_map;
0126
0127         // note : for the final Truth
0128         std::map<std::string, TH1D*> h1D_TruedNdEta_map;
0129
0130         // note : for the alpha correction
0131         // note : X: eta, Y: Mbin
0132         std::map<std::string, TH2D*> h2D_alpha_correction_map_in;
0133         std::map<std::string, TH2D*> h2D_alpha_correction_map_out;
0134         std::map<std::string, int> alpha_correction_name_map = {
0135             {"h2D_FitBkg_alpha_correction", 1},
0136             {"h2D_RotatedBkg_alpha_correction", 1},
0137             {"h2D_typeA_FitBkg_alpha_correction", 1},
0138             {"h2D_typeA_RotatedBkg_alpha_correction", 1}
0139         };
0140
0141         // note : ----------------- for the fitting -----------------
0142         std::map<std::string, TF1*> f1_BkgPol2_Fit_map;
0143         std::map<std::string, TF1*> f1_BkgPol2_Draw_map;
0144         std::map<std::string, TF1*> f1_SigBkgPol2_Fit_map;
0145         std::map<std::string, TF1*> f1_SigBkgPol2_DrawSig_map;
0146         std::map<std::string, TF1*> f1_SigBkgPol2_DrawBkgPol2_map;
0147
0148         // note : ----------------- for the output file -----------------
0149         TFile * file_out;
0150         TFile * file_out_dNdEta;
0151         std::string output_filename;
0152         void PrepareOutPutFileName();
0153         void PrepareOutPutRootFile();
0154
0155         TCanvas * c1;
0156
0157         // note : ----------------- for the analysis -----------------
0158         std::map<int, int> GetVtxZIndexMap();
0159         std::tuple<int, int, int, int, int, int> GetHistStringInfo(std::string hist_name);
0160         std::vector<double> find_Ngroup(TH1D * hist_in);
0161         double get_dist_offset(TH1D * hist_in, int check_N_bin);
0162         double get_h2D_GoodProtoTracklet_count(TH2D * hist_in, int eta_bin_in);
0163         double get_EvtCount(TH2D * hist_in, int centrality_bin_in);
0164         void Convert_to_PerEvt(TH1D * hist_in, double Nevent);
0165
0166         std::map<int, int> vtxZ_index_map;
0167
0168         // note : ----------------- for some constant -----------------
0169         std::pair<double, double> cut_INTTvtxZ = {-10, 10};
0170         const std::vector<int> ROOT_color_code = {
0171             51, 61, 70, 79, 88, 98, 100, 113, 120, 129
0172         };
0173
0174         int Semi_inclusive_Mbin = 7;
0175
0176         // note : ----------------- for the functions -----------------
0177         static double gaus_func(double *x, double *par)
0178         {
0179             // note : par[0] : size
0180             // note : par[1] : mean
0181             // note : par[2] : width
0182             // note : par[3] : offset
0183             return par[0] * TMath::Gaus(x[0],par[1],par[2]); //+ par[3];
0184         }
0185
0186         static  double gaus_pol2_func(double *x, double *par)
0187         {
0188             // note : par[0] : size
0189             // note : par[1] : mean
0190             // note : par[2] : width
0191
0192             double gaus_func = par[0] * TMath::Gaus(x[0],par[1],par[2]);
0193             double pol2_func = par[3] + par[4]* (x[0]-par[6]) + par[5] * pow((x[0]-par[6]),2);
0194
0195             return gaus_func + pol2_func;
0196         }
0197
0198         static  double bkg_pol2_func(double *x, double *par)
0199         {
0200             if (x[0] > (-1 * par[4]) && x[0] < par[5]) {
0201                 TF1::RejectPoint();
0202                 return 0;
0203             }
0204             return par[0] + par[1]* (x[0]-par[3]) + par[2] * pow((x[0]-par[3]),2);
0205
0206             // 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
0207         }
0208
0209         static  double full_pol2_func(double *x, double *par)
0210         {
0211             return par[0] + par[1]* (x[0]-par[3]) + par[2] * pow((x[0]-par[3]),2);
0212
0213             // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2
0214         }
0215
0216
0217
0218 };
0219
0220 #endif // PREPAREDNDETA_H