analysis/INTTvtxZTrapezoidal/InttZVertexFinderTrapezoidal.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef INTTZVERTEXFINDERTRAPEZOIDAL_H
0004 #define INTTZVERTEXFINDERTRAPEZOIDAL_H
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <trackbase/ActsGeometry.h>
0010 #include <trackbase/InttDefs.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrHit.h>
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainer.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021
0022 #include <string>
0023 #include <cmath>
0024 #include <cstdlib>
0025 #include <filesystem>
0026 #include <iterator>
0027 #include <algorithm>
0028
0029
0030 #include <TH1.h>
0031 #include <TF1.h>
0032 #include <TMath.h>
0033 #include <TColor.h>
0034
0035 #include "INTTvtxZF4AObj.h"
0036
0037 class PHCompositeNode;
0038 class INTTvtxZF4AObj;
0039
0040
0041 class InttZVertexFinderTrapezoidal : public SubsysReco
0042 {
0043  public:
0044
0045   InttZVertexFinderTrapezoidal(
0046     const std::string &name = "InttZVertexFinderTrapezoidal",
0047     std::pair<double, double> vertexXYIncm_in = {0, 0}, // note : in cm
0048     bool IsFieldOn_in = false,
0049     bool IsDCACutApplied_in = true,
0050     std::pair<std::pair<double,double>,std::pair<double,double>> DeltaPhiCutInDegree_in = {{-1,1},{-1000.,1000.}}, // note : in degree
0051     std::pair<std::pair<double,double>,std::pair<double,double>> DCAcutIncm_in = {{-1,1},{-1000.,1000.}}, // note : in cm
0052     int ClusAdcCut_in = 35,
0053     int ClusPhiSizeCut_in = 8,
0054     bool PrintRecoDetails_in = false
0055   );
0056
0057   ~InttZVertexFinderTrapezoidal() override;
0058
0059   /** Called during initialization.
0060       Typically this is where you can book histograms, and e.g.
0061       register them to Fun4AllServer (so they can be output to file
0062       using Fun4AllServer::dumpHistos() method).
0063    */
0064   int Init(PHCompositeNode *topNode) override;
0065
0066   /** Called for first event when run number is known.
0067       Typically this is where you may want to fetch data from
0068       database, because you know the run number. A place
0069       to book histograms which have to know the run number.
0070    */
0071   int InitRun(PHCompositeNode *topNode) override;
0072
0073   /** Called for each event.
0074       This is where you do the real work.
0075    */
0076   int process_event(PHCompositeNode *topNode) override;
0077
0078   /// Clean up internals after each event.
0079   int ResetEvent(PHCompositeNode *topNode) override;
0080
0081   /// Called at the end of each run.
0082   int EndRun(const int runnumber) override;
0083
0084   /// Called at the end of all processing.
0085   int End(PHCompositeNode *topNode) override;
0086
0087   /// Reset
0088   int Reset(PHCompositeNode * /*topNode*/) override;
0089
0090   void Print(const std::string &what = "ALL") const override;
0091
0092  private:
0093     std::pair<double, double> vertexXYIncm;
0094     bool IsFieldOn;
0095     bool IsDCACutApplied;
0096     std::pair<std::pair<double,double>,std::pair<double,double>> DeltaPhiCutInDegree; // note : this is prepared for the field on also
0097     std::pair<std::pair<double,double>,std::pair<double,double>> DCAcutIncm; // note : this is prepared for the field on also
0098     int ClusAdcCut;
0099     int ClusPhiSizeCut;
0100     bool PrintRecoDetails;
0101
0102     // note : prepare clsuter vec
0103     int PrepareClusterVec(PHCompositeNode *topNode);
0104     int PrepareINTTvtxZ();
0105
0106     // note : the private member functions
0107     double get_radius(double x, double y);
0108     double get_delta_phi(double angle_1, double angle_2);
0109     double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY);
0110     double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y);
0111     std::pair<double,double> Get_possible_zvtx(double rvtx, std::vector<double> p0, std::vector<double> p1);
0112     void line_breakdown(TH1D* hist_in, std::pair<double,double> line_range);
0113     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);
0114     std::vector<double> find_Ngroup(TH1D * hist_in);
0115     double vector_average (std::vector <double> input_vector);
0116     double vector_stddev (std::vector <double> input_vector);
0117
0118     // note : create node tree
0119     int createNodes(PHCompositeNode *topNode);
0120     INTTvtxZF4AObj * inttzobj;
0121
0122     // note : the cluster structure
0123     struct clu_info{
0124         double x;
0125         double y;
0126         double z;
0127
0128         int adc;
0129         int phi_size;
0130         int sensorZID;
0131     };
0132
0133     static double gaus_func(double *x, double *par)
0134     {
0135         // note : par[0] : size
0136         // note : par[1] : mean
0137         // note : par[2] : width
0138         // note : par[3] : offset
0139         return par[0] * TMath::Gaus(x[0],par[1],par[2]) + par[3];
0140     };
0141
0142
0143     // note : objects for the vertex Z reco
0144     TH1D * evt_possible_z;
0145     TH1D * line_breakdown_hist;
0146     TH1D * combination_zvtx_range_shape;
0147     std::vector<clu_info> temp_sPH_inner_nocolumn_vec;
0148     std::vector<clu_info> temp_sPH_outer_nocolumn_vec;
0149     std::vector<std::vector<std::pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0150     std::vector<std::vector<std::pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0151     // std::pair<double, double> m_vertexXYInmm;
0152     // std::pair<std::pair<double,double>,std::pair<double,double>> m_DCAcutInmm;
0153
0154     TF1 * gaus_fit; // note : cm
0155     std::vector<TF1 *> gaus_fit_vec; // note : cm
0156     std::vector<double> fit_mean_mean_vec;
0157     std::vector<double> fit_mean_reducedChi2_vec;
0158     std::vector<double> fit_mean_width_vec;
0159
0160
0161     // note : for the INTT vertex Z out
0162     // note : unit [cm]
0163     double INTTvtxZ;
0164     double INTTvtxZError;
0165     double NgroupTrapezoidal;
0166     double NgroupCoarse;
0167     double TrapezoidalFitWidth;
0168     double TrapezoidalFWHM;
0169     long long NClusAll;
0170     long long NClusAllInner;
0171
0172
0173     // note: some constants
0174     const int zvtx_hist_Nbins = 1200;   // note : N bins for each side, regardless the bin at zero
0175     const double fine_bin_width = 0.05;  // note : bin width with the unit [cm]
0176     const std::pair<double,double> evt_possible_z_range = {-70, 70}; // note : [cm]
0177     const std::pair<double,double> edge_rejection = {-50, 50}; // note : [cm]
0178     const double typeA_sensor_half_length_incm = 0.8; // note : [cm]
0179     const double typeB_sensor_half_length_incm = 1.0; // note : [cm]
0180     const std::vector<int> typeA_sensorZID = {0,2}; // note : sensor Z ID for type A // note -> 1, 0, 2, 3
0181     const std::vector<std::string> color_code = {
0182         "#9e0142",
0183         "#d53e4f",
0184         "#f46d43",
0185         "#fdae61",
0186         "#fee08b",
0187         "#e6f598",
0188         "#abdda4",
0189         "#66c2a5",
0190         "#3288bd",
0191         "#5e4fa2"
0192     };
0193
0194     // note : constants for the INTT vtxZ
0195     const int number_of_gaus = 7;
0196
0197 };
0198
0199 #endif // INTTZVERTEXFINDERTRAPEZOIDAL_H