NewCode2024/EvtVtxZTracklet/EvtVtxZProtoTracklet.h

0001 #ifndef EVTVTXZPROTOTRACKLET_H
0002 #define EVTVTXZPROTOTRACKLET_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 "structure.h"
0032 #include "../Constants.h"
0033
0034 class EvtVtxZProtoTracklet{
0035     public:
0036         EvtVtxZProtoTracklet(
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             std::string output_file_name_suffix_in,
0044
0045             std::pair<double, double> vertexXYIncm_in = {0, 0}, // note : in cm
0046             bool IsFieldOn_in = false,
0047             bool IsDCACutApplied_in = true,
0048             std::pair<std::pair<double,double>,std::pair<double,double>> DeltaPhiCutInDegree_in = {{-1,1},{-1000.,1000.}}, // note : in degree
0049             std::pair<std::pair<double,double>,std::pair<double,double>> DCAcutIncm_in = {{-1,1},{-1000.,1000.}}, // note : in cm
0050             int ClusAdcCut_in = 35,
0051             int ClusPhiSizeCut_in = 8,
0052
0053             bool PrintRecoDetails_in = false,
0054             bool DrawEvtVtxZ_in = true,
0055
0056             bool RunInttBcoFullDiff_in = true,
0057             bool RunVtxZReco_in = true,
0058             bool RunTrackletPair_in = true,
0059             bool RunTrackletPairRotate_in = true,
0060
0061             bool HaveGeoOffsetTag_in = false
0062         );
0063
0064         std::string GetOutputFileName() {return output_filename;}
0065         void SetGeoOffset(std::map<std::string, std::vector<double>> input_geo_offset_map);
0066         void MainProcess();
0067         void EndRun(int close_file_in = 1);
0068         TH1D * GetINTTRecovtxZH1D(){
0069             return reco_INTTvtxZ;
0070         }
0071
0072     protected:
0073         struct clu_info{
0074             double x;
0075             double y;
0076             double z;
0077
0078             int index;
0079
0080             int adc;
0081             int phi_size;
0082             int sensorZID;
0083             int ladderPhiID;
0084             int layerID;
0085         };
0086
0087     // note : ---------------------- For the constructor ----------------------
0088     int process_id;
0089     int runnumber;
0090     int run_nEvents;
0091     std::string input_directory;
0092     std::string input_file_name;
0093     std::string output_directory;
0094     std::string output_file_name_suffix;
0095
0096     std::pair<double, double> vertexXYIncm;
0097     bool IsFieldOn;
0098     bool IsDCACutApplied;
0099     std::pair<std::pair<double,double>,std::pair<double,double>> DeltaPhiCutInDegree;
0100     std::pair<std::pair<double,double>,std::pair<double,double>> DCAcutIncm;
0101     int ClusAdcCut;
0102     int ClusPhiSizeCut;
0103     bool PrintRecoDetails;
0104     bool DrawEvtVtxZ;
0105     bool RunInttBcoFullDiff;
0106     bool RunVtxZReco;
0107     bool RunTrackletPair;
0108     bool RunTrackletPairRotate;
0109     bool HaveGeoOffsetTag;
0110
0111
0112     // note : ---------------------- For read & outputTFile ----------------------
0113     TFile * file_in;
0114     TTree * tree_in;
0115     TFile * file_out;
0116     TTree * tree_out;
0117
0118     void PrepareOutPutFileName();
0119     std::string output_filename;
0120     void PrepareRootFile();
0121     std::map<std::string, int> GetInputTreeBranches(TTree * m_tree_in);
0122
0123     // int event;
0124     float MBD_z_vtx;
0125     ULong_t INTT_BCO;
0126     std::vector<float> *ClusX;
0127     std::vector<float> *ClusY;
0128     std::vector<float> *ClusZ;
0129     std::vector<int> *ClusLayer;
0130     std::vector<unsigned char> *ClusLadderZId;
0131     std::vector<unsigned char> *ClusLadderPhiId;
0132     std::vector<int> *ClusAdc;
0133     std::vector<float> *ClusPhiSize;
0134     std::vector<int> *firedTriggers;
0135
0136     float TruthPV_trig_x;
0137     float TruthPV_trig_y;
0138     float TruthPV_trig_z;
0139
0140     double INTTvtxZ;
0141     double INTTvtxZError;
0142
0143     // note : additional branches used for making the vertex Z distribution
0144     int NTruthVtx;
0145     bool is_min_bias;
0146     float MBD_centrality;
0147
0148     // note : the flag
0149     int m_withTrig = false;
0150
0151     // note : ---------------------- For update tree vertex Z ----------------------
0152     double out_INTTvtxZ;
0153     double out_INTTvtxZError;
0154     double out_NgroupTrapezoidal;
0155     double out_NgroupCoarse;
0156     double out_TrapezoidalFitWidth;
0157     double out_TrapezoidalFWHM;
0158
0159     TBranch * b_INTTvtxZ;
0160     TBranch * b_INTTvtxZError;
0161     TBranch * b_NgroupTrapezoidal;
0162     TBranch * b_NgroupCoarse;
0163     TBranch * b_TrapezoidalFitWidth;
0164     TBranch * b_TrapezoidalFWHM;
0165
0166     // note : ---------------------- For the cluster eta and phi ----------------------
0167     std::vector<double> out_ClusEta_INTTz;
0168     std::vector<double> out_ClusEta_MBDz;
0169     std::vector<double> out_ClusEta_TrueXYZ;
0170     std::vector<double> out_ClusPhi_AvgPV;
0171     std::vector<double> out_ClusPhi_TrueXY;
0172
0173     TBranch * b_ClusEta_INTTz;
0174     TBranch * b_ClusEta_MBDz;
0175     TBranch * b_ClusEta_TrueXYZ;
0176
0177     TBranch * b_ClusPhi_AvgPV;
0178     TBranch * b_ClusPhi_TrueXY;
0179
0180     void GetClusUpdated();
0181
0182     // note : ---------------------- For InttBcoFullDiff w.r.t next ----------------------
0183     int out_InttBcoFullDiff_next;
0184     TBranch * b_InttBcoFullDiff_next;
0185
0186     // note : ---------------------- For common ----------------------
0187     double CheckGeoOffsetMap();
0188     void PrepareClusterVec();
0189     void EvtCleanUp();
0190     std::vector<clu_info> evt_sPH_inner_nocolumn_vec;
0191     std::vector<clu_info> evt_sPH_outer_nocolumn_vec;
0192     std::vector<std::vector<std::pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0193     std::vector<std::vector<std::pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0194     std::map<std::string, std::vector<double>> geo_offset_map;
0195     double temp_INTTvtxZ;
0196     double temp_INTTvtxZError;
0197
0198     int out_MBDNSg2;
0199     int out_MBDNSg2_vtxZ10cm;
0200     int out_MBDNSg2_vtxZ30cm;
0201     int out_MBDNSg2_vtxZ60cm;
0202     int out_eID_count;
0203
0204     TBranch * b_MBDNSg2;
0205     TBranch * b_MBDNSg2_vtxZ10cm;
0206     TBranch * b_MBDNSg2_vtxZ30cm;
0207     TBranch * b_MBDNSg2_vtxZ60cm;
0208     TBranch * b_eID_count;
0209
0210     int index_MBDNSg2 = 10;
0211     int index_MBDNSg2_vtxZ10cm = 12;
0212     int index_MBDNSg2_vtxZ30cm = 13;
0213     int index_MBDNSg2_vtxZ60cm = 14;
0214
0215     void GetTriggerInfo();
0216
0217     // note : ---------------------- For vertex Z calculation ----------------------
0218     void PrepareINTTvtxZ();
0219     void GetINTTvtxZ();
0220
0221     double get_radius(double x, double y);
0222     double get_delta_phi(double angle_1, double angle_2);
0223     double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0224     double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y);
0225     std::pair<double,double> Get_possible_zvtx(double rvtx, std::vector<double> p0, std::vector<double> p1);
0226     void line_breakdown(TH1D* hist_in, std::pair<double,double> line_range);
0227     void trapezoidal_line_breakdown(TH1D * hist_in, double inner_r, double inner_z, int inner_zid, double outer_r, double outer_z, int outer_zid);
0228     std::vector<double> find_Ngroup(TH1D * hist_in);
0229     double vector_average (std::vector <double> input_vector);
0230     double vector_stddev (std::vector <double> input_vector);
0231
0232     static double gaus_func(double *x, double *par)
0233     {
0234         // note : par[0] : size
0235         // note : par[1] : mean
0236         // note : par[2] : width
0237         // note : par[3] : offset
0238         return par[0] * TMath::Gaus(x[0],par[1],par[2]) + par[3];
0239     };
0240
0241     // note : objects for the vertex Z reco
0242     TH1D * evt_possible_z;
0243     TH1D * line_breakdown_hist;
0244     TH1D * combination_zvtx_range_shape;
0245     std::vector<clu_info> evt_sPH_inner_nocolumn_vec_PostCut;
0246     std::vector<clu_info> evt_sPH_outer_nocolumn_vec_PostCut;
0247     std::vector<std::vector<std::pair<bool,clu_info>>> inner_clu_phi_map_PostCut; // note: phi
0248     std::vector<std::vector<std::pair<bool,clu_info>>> outer_clu_phi_map_PostCut; // note: phi
0249     // std::pair<double, double> m_vertexXYInmm;
0250     // std::pair<std::pair<double,double>,std::pair<double,double>> m_DCAcutInmm;
0251
0252     TF1 * gaus_fit; // note : cm
0253     std::vector<TF1 *> gaus_fit_vec; // note : cm
0254     std::vector<double> fit_mean_mean_vec;
0255     std::vector<double> fit_mean_reducedChi2_vec;
0256     std::vector<double> fit_mean_width_vec;
0257
0258     // note : ---------------------- For DrawEvtVtxZ ----------------------
0259     TFile * INTTvtxZ_EvtDisplay_file_out = nullptr;
0260     TH1D * line_breakdown_hist_zoomin;
0261     TCanvas * c1;
0262     TPad * pad_EvtZDist;
0263     TPad * pad_ZoomIn_EvtZDist;
0264     TLatex * draw_text;
0265     void PrepareEvtCanvas();
0266     void Characterize_Pad(TPad *pad, float left = 0.15, float right = 0.1, float top = 0.1, float bottom = 0.12, bool set_logY = false, int setgrid_bool = 0);
0267
0268     // note : ---------------------- For update tree tracklet pair ----------------------
0269     std::vector<pair_str> out_evt_TrackletPair_vec;
0270     std::vector<pair_str> out_evt_TrackletPairRotate_vec;
0271
0272     TBranch * b_evt_TrackletPair_vec;
0273     TBranch * b_evt_TrackletPairRotate_vec;
0274
0275     double rotate_phi_angle = 180.; // note : unit : degree
0276
0277     TGraphErrors * track_gr;
0278     TF1 * fit_rz;
0279
0280     void GetTrackletPair(std::vector<pair_str> &input_TrackletPair_vec, bool isRotated);
0281     std::vector<clu_info> GetRoatedClusterVec(std::vector<clu_info> input_cluster_vec);
0282     std::pair<double,double> rotatePoint(double x, double y);
0283     double get_clu_eta(std::vector<double> vertex, std::vector<double> clu_pos);
0284     std::pair<double,double> Get_eta(std::vector<double>p0, std::vector<double>p1, std::vector<double>p2);
0285     double grEY_stddev(TGraphErrors * input_grr);
0286     std::pair<double, double> mirrorPolynomial(double a, double b);
0287
0288     // note : ---------------------- It's to have a vertex Z distribution ----------------------
0289     TH1D * reco_INTTvtxZ;
0290     void FillRecoINTTVtxZH1D(int event_index);
0291     int nVtxZ_bin = 60;
0292     std::pair<double, double> vtxZ_range = {-60, 60};
0293
0294     std::pair<double, double> cut_vtxZDiff = Constants::cut_vtxZDiff;
0295     std::pair<double, double> cut_TrapezoidalFitWidth = Constants::cut_TrapezoidalFitWidth;
0296     std::pair<double, double> cut_TrapezoidalFWHM = Constants::cut_TrapezoidalFWHM;
0297     std::pair<double, double> cut_GlobalMBDvtxZ = Constants::cut_GlobalMBDvtxZ;
0298
0299     int cut_InttBcoFullDIff_next = Constants::cut_InttBcoFullDIff_next;
0300     int cut_MBD_centrality = 70;
0301
0302     // note : ---------------------- For constants ----------------------
0303     const int B0L0_index = 3;
0304     const int B0L1_index = 4;
0305     const int B1L0_index = 5;
0306     const int B1L1_index = 6;
0307     const int nLadder_inner = 12;
0308     const int nLadder_outer = 16;
0309
0310     const int zvtx_hist_Nbins = 1200;   // note : N bins for each side, regardless the bin at zero
0311     const double fine_bin_width = 0.05;  // note : bin width with the unit [cm]
0312     const std::pair<double,double> evt_possible_z_range = {-70, 70}; // note : [cm]
0313     const std::pair<double,double> edge_rejection = {-50, 50}; // note : [cm]
0314     const double typeA_sensor_half_length_incm = 0.8; // note : [cm]
0315     const double typeB_sensor_half_length_incm = 1.0; // note : [cm]
0316     const std::vector<int> typeA_sensorZID = {0,2}; // note : sensor Z ID for type A // note -> 1, 0, 2, 3
0317     const std::vector<std::string> color_code = {
0318         "#9e0142",
0319         "#d53e4f",
0320         "#f46d43",
0321         "#fdae61",
0322         "#fee08b",
0323         "#e6f598",
0324         "#abdda4",
0325         "#66c2a5",
0326         "#3288bd",
0327         "#5e4fa2"
0328     };
0329
0330     const std::vector<int> ROOT_color_code = {
0331         51, 61, 70, 79, 88, 98, 100, 113, 120, 129
0332     };
0333
0334     // note : constants for the INTT vtxZ
0335     const int number_of_gaus = 7;
0336
0337 };
0338
0339 #endif