PrepareAllDist/RestDist/RestDist.h

0001 #ifndef RESTDIST_H
0002 #define RESTDIST_H
0003
0004 #include <iostream>
0005 #include <vector>
0006 #include <string>
0007 #include <numeric>
0008
0009 #include <TFile.h>
0010 #include <TTree.h>
0011 #include <TH2D.h>
0012 #include <TH1D.h>
0013 #include <TMath.h>
0014 #include <TF1.h>
0015 #include <TProfile.h>
0016 #include <TCanvas.h>
0017 #include <TPad.h>
0018 #include <TGraphErrors.h>
0019 #include <TLatex.h>
0020
0021 #include <TCanvas.h> // note : for the combined case
0022 #include <TGraph.h>  // note : for the combined case
0023
0024 #include <TRandom.h> // note : for the offset
0025 #include <TRandom3.h> // note : for the offset
0026
0027 #include <TColor.h>
0028
0029 #include <TObjArray.h>
0030
0031 #include "../../Constants.h"
0032 #include "../../EvtVtxZTracklet/structure.h"
0033
0034 class RestDist{
0035     public:
0036         RestDist(
0037             int process_id_in,
0038             int runnumber_in,
0039             int run_nEvents_in,
0040             std::string input_directory_in,
0041             std::string input_file_name_in,
0042             std::string output_directory_in,
0043
0044             std::string output_file_name_suffix_in = "",
0045             std::pair<double, double> vertexXYIncm_in = {std::nan(""),std::nan("")},
0046
0047             bool Apply_cut_in = false,
0048             bool ApplyVtxZReWeighting_in = false,
0049             std::pair<bool, int> ApplyEvtBcoFullDiffCut_in = {true, 61},
0050             std::pair<bool, std::pair<double,double>> RequireVtxZRange_in = {true, {-10, 10}},
0051             std::pair<bool, std::pair<double,double>> isClusQA_in = {false, {0, 0}},
0052             bool isRotated_in = false
0053         );
0054
0055         std::string GetOutputFileName() {return output_filename;}
0056         void SetINTTvtxZReweighting(TH1D * h1D_INTT_vtxZ_reweighting_in) {h1D_INTT_vtxZ_reweighting = (TH1D*)h1D_INTT_vtxZ_reweighting_in->Clone("INTT_vtxZ_reweighting");}
0057         void PrepareEvent();
0058         void SetInnerBarrelRotation(double angle_in) {rotate_phi_angle = angle_in;}
0059         void EndRun();
0060
0061     protected:
0062         // division : ---For constructor---------------------------------------------------------------------------------------------------------------------------------------------
0063         int process_id;
0064         int runnumber;
0065         int run_nEvents;
0066         std::string input_directory;
0067         std::string input_file_name;
0068         std::string output_directory;
0069         std::string output_file_name_suffix;
0070         std::pair<double, double> vertexXYIncm;
0071         bool Apply_cut;
0072         bool ApplyVtxZReWeighting;
0073         std::pair<bool, int> ApplyEvtBcoFullDiffCut;
0074         std::pair<bool, std::pair<double,double>> RequireVtxZRange;
0075         std::pair<bool, std::pair<double,double>> isClusQA;
0076         bool isRotated;
0077
0078         // division : ---For input file---------------------------------------------------------------------------------------------------------------------------------------------
0079         void PrepareInputFile();
0080         std::map<std::string, int> GetInputTreeBranches(TTree * m_tree_in);
0081         TFile * file_in;
0082         TTree * tree_in;
0083         std::string tree_name = "EventTree";
0084
0085         bool is_min_bias;
0086         float MBD_centrality;
0087         float MBD_z_vtx;
0088         float MBD_charge_sum;
0089         float MBD_south_charge_sum;
0090         float MBD_north_charge_sum;
0091         float MBD_charge_asymm;
0092
0093         double INTTvtxZ;
0094         double INTTvtxZError;
0095         double NgroupTrapezoidal;
0096         double NgroupCoarse;
0097         double TrapezoidalFitWidth;
0098         double TrapezoidalFWHM;
0099
0100         int NClus;
0101         int NClus_Layer1;
0102
0103         std::vector<float> * ClusX;
0104         std::vector<float> * ClusY;
0105         std::vector<float> * ClusZ;
0106
0107         std::vector<unsigned int> * ClusAdc;
0108         std::vector<float> * ClusPhiSize;
0109         std::vector<float> * ClusZSize;
0110
0111         std::vector<int> * ClusLayer;
0112         std::vector<unsigned char> * ClusLadderZId;
0113         std::vector<unsigned char> * ClusLadderPhiId;
0114
0115         std::vector<double> * ClusEta_INTTz;
0116         std::vector<double> * ClusEta_MBDz;
0117         std::vector<double> * ClusPhi_AvgPV;
0118
0119         // note : for MC
0120         std::vector<double> * ClusEta_TrueXYZ;
0121         std::vector<double> * ClusPhi_TrueXY;
0122
0123
0124         // note : for data?
0125         int MBDNSg2;
0126         int MBDNSg2_vtxZ10cm;
0127         int MBDNSg2_vtxZ30cm;
0128
0129         int InttBcoFullDiff_next; // note : for data
0130
0131         // note : for MC
0132         int NTruthVtx;
0133         float TruthPV_trig_x;
0134         float TruthPV_trig_y;
0135         float TruthPV_trig_z;
0136
0137
0138         // note : the flag
0139         int m_withTrig = false;
0140
0141         // division : ---For Analysis---------------------------------------------------------------------------------------------------------------------------------------------
0142         struct clu_info{
0143             double x;
0144             double y;
0145             double z;
0146
0147             int index;
0148
0149             int adc;
0150             int phi_size;
0151             int sensorZID;
0152             int ladderPhiID;
0153             int layerID;
0154
0155             double eta_INTTz;
0156             double eta_MBDz;
0157             double phi_avgXY;
0158
0159             double eta_TrueXYZ;
0160             double phi_TrueXY;
0161         };
0162
0163         int NClus_inner_typeA;
0164         int NClus_inner_typeB;
0165         int NClus_outer_typeA;
0166         int NClus_outer_typeB;
0167
0168         void PrepareClusterVec();
0169         void EvtCleanUp();
0170         std::vector<clu_info> evt_sPH_inner_nocolumn_vec;
0171         std::vector<clu_info> evt_sPH_outer_nocolumn_vec;
0172         std::vector<std::vector<std::pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0173         std::vector<std::vector<std::pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0174
0175         double shortestDistanceLineToPoint(const std::vector<double> point1, const std::vector<double> point2, const std::vector<double> target);
0176         void GetRotatedClusterVec(std::vector<RestDist::clu_info> &input_cluster_vec);
0177         std::pair<double,double> rotatePoint(double x, double y);
0178         double get_clu_eta(std::vector<double> vertex, std::vector<double> clu_pos);
0179         double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0180
0181         // note : for the best pair quick check
0182         std::vector<pair_str> temp_all_pair_vec;
0183         std::vector<double> temp_all_pair_deltaR_vec;
0184         std::map<int, TH2D*> BestPair_h2D_map;
0185         std::map<int, TGraph*> BestPair_gr_map;
0186
0187         // division : ---For output file---------------------------------------------------------------------------------------------------------------------------------------------
0188         void PrepareOutPutFileName();
0189         void PrepareOutputFile();
0190         std::string output_filename;
0191
0192         TFile * file_out;
0193
0194         // division : ---For Prepare Hist---------------------------------------------------------------------------------------------------------------------------------------------
0195         void PrepareHist();
0196         std::map<std::string, TH1D*> h1D_map;
0197         std::map<std::string, TH2D*> h2D_map;
0198         TH1D * h1D_INTT_vtxZ_reweighting = nullptr; // note : should be long range
0199         TH1D * h1D_centrality_bin;
0200
0201         // division : ---For constants---------------------------------------------------------------------------------------------------------------------------------------------
0202         const std::vector<int> typeA_sensorZID = {0,2}; // note : sensor Z ID for type A // note -> 1, 0, 2, 3
0203
0204         int HighNClus = Constants::HighNClus;
0205
0206         double rotate_phi_angle = 180.;
0207
0208         std::vector<double> centrality_edges = Constants::centrality_edges;
0209         int nCentrality_bin;
0210
0211         // int nVtxZ_bin = 60;
0212         // std::pair<double, double> vtxZ_range = {-60, 60};
0213
0214         // int nVtxZ_bin_narrow = 20;
0215         // std::pair<double, double> vtxZ_range_narrow = {-20, 20};
0216
0217         // note : for deta_phi
0218         double DeltaPhiEdge_min = -0.07; // note : rad ~ -4 degree
0219         double DeltaPhiEdge_max = 0.07;  // note : rad ~ 4 degree
0220         int    nDeltaPhiBin = 140;
0221
0222         // note : for deta_eta
0223         double DeltaEtaEdge_min = -1.; // note : rad ~ -4 degree
0224         double DeltaEtaEdge_max = 1.;  // note : rad ~ 4 degree
0225         int    nDeltaEtaBin = 100;
0226
0227         double EtaEdge_min = Constants::EtaEdge_min;
0228         double EtaEdge_max = Constants::EtaEdge_max;
0229         int nEtaBin = Constants::nEtaBin;
0230
0231
0232         std::pair<double, double> cut_vtxZDiff = Constants::cut_vtxZDiff;
0233         std::pair<double, double> cut_TrapezoidalFitWidth = Constants::cut_TrapezoidalFitWidth;
0234         std::pair<double, double> cut_TrapezoidalFWHM = Constants::cut_TrapezoidalFWHM;
0235         std::pair<double, double> cut_INTTvtxZError = Constants::cut_INTTvtxZError;
0236         std::pair<double, double> cut_GlobalMBDvtxZ = Constants::cut_GlobalMBDvtxZ;
0237         int cut_InttBcoFullDIff_next = Constants::cut_InttBcoFullDIff_next;
0238         int Semi_inclusive_bin = Constants::Semi_inclusive_bin;
0239
0240         std::vector<double> data_edep_bins_vec = {
0241             0, 35, 45, 60, 90, 120, 150, 180, 210, 300
0242         };
0243
0244         std::vector<double> MC_edep_bins_vec = {
0245             0, 15, 30, 60, 90, 120, 150, 180, 210, 300
0246         };
0247
0248 };
0249
0250 #endif