NewCode2024/AvgVtxXY/AvgVtxXY.h

0001 #ifndef AVGVTXXY_H
0002 #define AVGVTXXY_H
0003
0004 #include <iostream>
0005 #include <vector>
0006 #include <string>
0007
0008 #include <TFile.h>
0009 #include <TTree.h>
0010 #include <TH2D.h>
0011 #include <TH1D.h>
0012 #include <TMath.h>
0013 #include <TF1.h>
0014 #include <TProfile.h>
0015
0016 #include <TCanvas.h> // note : for the combined case
0017 #include <TGraph.h>  // note : for the combined case
0018
0019 #include <TRandom.h> // note : for the offset
0020 #include <TRandom3.h> // note : for the offset
0021
0022
0023
0024 class AvgVtxXY {
0025     public:
0026         AvgVtxXY(
0027             int process_id_in,
0028             int runnumber_in,
0029             int run_nEvents_in,
0030             std::string input_directory_in,
0031             std::string input_file_name_in,
0032             std::string output_directory_in,
0033             std::string output_file_name_suffix_in,
0034
0035             // note : for the event selection
0036             std::pair<double,double> MBD_vtxZ_cut_in,
0037             std::pair<int,int> INTTNClus_cut_in,
0038             int ClusAdc_cut_in,
0039             int ClusPhiSize_cut_in,
0040
0041             bool HaveGeoOffsetTag_in = false,
0042             double random_range_XYZ_in = 0.02, // note : unit : cm // note : for the generation of the offset
0043             int random_seed_in = -999,  // note : for the generation of the offset
0044             std::string input_offset_map_in = "no_map" // note : full map
0045         );
0046
0047         std::map<std::string, std::vector<double>> GetGeoOffsetMap(){
0048             return geo_offset_map;
0049         };
0050         void PreparePairs();
0051         void FindVertexQuadrant(int nIteration, double search_Window_half, std::pair<double,double> set_center);
0052         std::pair<double,double> GetVertexQuadrant(){
0053             return std::make_pair(out_vtxX_quadrant, out_vtxY_quadrant);
0054         };
0055         void FindVertexLineFill(std::pair<double,double> set_center, int Nbins = 100, double search_Window_half = 0.25 /* unit : cm*/, double segmentation = 0.0001 /*unit: cm*/);
0056         void EndRun();
0057
0058         std::pair<double,double> GetFinalAvgVtxXY(){
0059             return std::make_pair((out_vtxX_quadrant + out_vtxX_LineFill) / 2., (out_vtxY_quadrant + out_vtxY_LineFill) / 2.);
0060         };
0061
0062         std::string GetOutputFileName() {return output_filename;}
0063
0064     protected:
0065         // note : --------------------- the structure ---------------------
0066         struct pos_xy{
0067             double x;
0068             double y;
0069         };
0070
0071         void ReadRootFile();
0072         void InitOutRootFile();
0073
0074         // note : --------------------- for the constructor ---------------------
0075         int process_id;
0076         int runnumber;
0077         int run_nEvents;
0078         std::string input_directory;
0079         std::string input_file_name;
0080         std::string output_directory;
0081         std::string output_file_name_suffix;
0082         std::pair<double,double> MBD_vtxZ_cut;
0083         std::pair<int,int> INTTNClus_cut;
0084         int ClusAdc_cut;
0085         int ClusPhiSize_cut;
0086         bool HaveGeoOffsetTag;
0087         double random_range_XYZ;
0088         int random_seed;
0089         std::string input_offset_map;
0090
0091         // note : --------------------- for reading the file ---------------------
0092         TFile * file_in;
0093         TTree * tree_in;
0094         std::string input_tree_name = "EventTree";
0095         float MBD_z_vtx;
0096         bool is_min_bias;
0097         float MBD_centrality;
0098         bool InttBco_IsToBeRemoved;
0099
0100         int NClus;
0101         std::vector<float> *ClusX;
0102         std::vector<float> *ClusY;
0103         std::vector<int> *ClusLayer;
0104         std::vector<unsigned char> *ClusLadderPhiId;
0105         std::vector<int> *ClusAdc;
0106         std::vector<float> *ClusPhiSize;
0107
0108         // note : --------------------- for output file ---------------------
0109         TFile * file_out;
0110         std::string output_filename;
0111
0112         TTree * tree_vtxXY;
0113         double out_vtxX_quadrant = std::nan("");
0114         double out_vtxXE_quadrant = std::nan("");
0115         double out_vtxY_quadrant = std::nan("");
0116         double out_vtxYE_quadrant = std::nan("");
0117         double out_vtxX_LineFill = std::nan("");
0118         double out_vtxXE_LineFill = std::nan("");
0119         double out_vtxXStdDev_LineFill = std::nan("");
0120         double out_vtxY_LineFill = std::nan("");
0121         double out_vtxYE_LineFill = std::nan("");
0122         double out_vtxYStdDev_LineFill = std::nan("");
0123         // int out_run_nEvents;
0124         int out_job_index;
0125         int out_file_total_event;
0126
0127         TTree * tree_geooffset;
0128         int out_LayerID;
0129         int out_LadderPhiID;
0130         double out_offset_x;
0131         double out_offset_y;
0132         double out_offset_z;
0133
0134         TTree * tree_quadrant_detail;
0135         int out_iteration;
0136         int out_quadrant;
0137         double out_assume_center_x;
0138         double out_assume_center_y;
0139         double out_Fit_InnerPhi_DeltaPhi_RChi2;
0140         double out_Fit_InnerPhi_DeltaPhi_Err0;
0141         double out_Fit_InnerPhi_DeltaPhi_LineY;
0142         double out_Fit_InnerPhi_DCA_RChi2;
0143         double out_Fit_InnerPhi_DCA_Err0;
0144         double out_Fit_InnerPhi_DCA_LineY;
0145
0146         TTree * tree_parameters;
0147         int out_run_nEvents;
0148         double out_MBD_vtxZ_cut_l;
0149         double out_MBD_vtxZ_cut_r;
0150         int out_INTTNClus_cut_l;
0151         int out_INTTNClus_cut_r;
0152         int out_ClusAdc_cut;
0153         int out_ClusPhiSize_cut;
0154         int out_HaveGeoOffsetTag;
0155         double out_random_range_XYZ;
0156         int out_random_seed;
0157         int out_input_offset_map;
0158         double out_hopeless_pair_angle_cut;
0159         double out_rough_DeltaPhi_cut_degree;
0160         double out_InnerOuterPhi_DeltaPhi_DCA_Threshold;
0161         double out_LineFill_Hist_Threshold;
0162         double out_InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold;
0163         int out_TH2D_threshold_advanced_2_CheckNBins;
0164
0165         // note : --------------------- for the geo offset ---------------------
0166         TRandom3 * geo_rand;
0167         std::map<std::string, std::vector<double>> geo_offset_map;
0168         std::vector<double> offset_correction(std::map<std::string,std::vector<double>> input_map);
0169         void GenGeoOffset(); // note : gen, save them into a map, and also in the tree
0170         void SetGeoOffset(); // note : read the map
0171         void geo_offset_correction_fill();
0172
0173         // note : --------------------- for preapre the pairs ---------------------
0174         std::vector<pos_xy> inner_cluster_vec;
0175         std::vector<pos_xy> outer_cluster_vec;
0176         std::vector<std::pair<pos_xy,pos_xy>> cluster_pair_vec;
0177         double hopeless_pair_angle_cut;
0178
0179         // note : --------------------- for the common ---------------------
0180         void TH2D_threshold_advanced_2(TH2D * hist, double threshold);
0181         int TH2D_threshold_advanced_2_CheckNBins;
0182
0183         // note : --------------------- run the iteration ---------------------
0184         double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0185         void SetUpHistograms(int iteration, int quadrant);
0186         std::pair<double,double> RunSingleCandidate(int iteration, int quadrant, std::pair<double,double> assume_center);
0187         std::vector<double> SumTH2FColumnContent(TH2D * hist_in);
0188         std::vector<std::pair<double,double>> Get4vtx(std::pair<double,double> origin, double length);
0189         std::map<std::string, TGraph*> gr_map; // note : for each iteration
0190         std::map<std::string, TProfile*> profile_map; // note : for each iteration
0191         TF1 * fit_innerphi_DeltaPhi;
0192         TF1 * fit_innerphi_DCA;
0193         double InnerOuterPhi_DeltaPhi_DCA_Threshold;
0194         double InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold;
0195
0196         // note : --------------------- for the histograms ---------------------
0197         std::map<std::string, TH2D*> h2D_map; // note : for each iteration
0198         std::map<std::string, TH1D*> h1D_map; // note : for each iteration
0199
0200         TH2D * final_InnerPhi_DeltaPhi;
0201         TH2D * final_InnerPhi_DCA;
0202         TH2D * final_InnerPhi_DeltaPhi_bkgrm;
0203         TH2D * final_InnerPhi_DCA_bkgrm;
0204
0205         TH2D * final_OuterPhi_DeltaPhi;
0206         TH2D * final_OuterPhi_DCA;
0207         TH2D * final_OuterPhi_DeltaPhi_bkgrm;
0208         TH2D * final_OuterPhi_DCA_bkgrm;
0209
0210         TH2D * final_DeltaPhi_DCA;
0211         TH1D * final_DeltaPhi_1D;
0212         TH1D * final_DCA_1D;
0213
0214         TH2D * final_LineFill_BeamSpot;
0215         TH2D * final_LineFill_BeamSpot_bkgrm;
0216
0217         // note : --------------------- for the line filled method ---------------------
0218         double rough_DeltaPhi_cut_degree;
0219         void TH2DSampleLineFill(TH2D * hist_in, double segmentation, std::pair<double,double> inner_clu, std::pair<double,double> outer_clu);
0220         std::vector<double> calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y);
0221         double LineFill_Hist_Threshold;
0222
0223         // note : --------------------- constants ---------------------
0224         const int h2_nbins = 100;
0225
0226         const double Hist2D_angle_xmin = -M_PI;
0227         const double Hist2D_angle_xmax =  M_PI;
0228
0229         const double Hist2D_DeltaPhi_min = -0.03;
0230         const double Hist2D_DeltaPhi_max =  0.03;
0231
0232         const double Hist2D_DCA_min = -1.0;
0233         const double Hist2D_DCA_max =  1.0;
0234
0235         const int B0L0_index = 3;
0236         const int B0L1_index = 4;
0237         const int B1L0_index = 5;
0238         const int B1L1_index = 6;
0239         const int nLadder_inner = 12;
0240         const int nLadder_outer = 16;
0241
0242
0243 };
0244
0245
0246
0247 #endif