analysis_INTT_CW/DST_MC/MegaTrackFinder.h

0001 #ifndef MegaTrackFinder_h
0002 #define MegaTrackFinder_h
0003
0004 #include "INTTXYvtxEvt.h"
0005
0006 class MegaTrackFinder : public INTTXYvtxEvt
0007 {
0008     public:
0009
0010         // note : the N_centrality_bin here means the number without the "inclusive"
0011         MegaTrackFinder(string run_type, string out_folder_directory, int N_centrality_bin, pair<double,double> beam_origin_in, bool mark_used = true) :
0012         INTTXYvtxEvt(run_type, out_folder_directory), N_centrality_bin(N_centrality_bin), mark_used(mark_used)
0013         {
0014             inner_used_mega_clu.clear();
0015             outer_used_mega_clu.clear();
0016
0017             clu4_track_performance.clear();
0018             clu4_track_index.clear();
0019             clu4_track_clu_phi.clear();
0020
0021             clu3_track_performance.clear();
0022             clu3_track_index.clear();
0023             clu3_track_clu_phi.clear();
0024             clu3_track_type.clear();
0025
0026             mega_track_eta.clear();
0027
0028             INTTXYvtxEvt::beam_origin = beam_origin_in;
0029
0030             NClu3_track_count = 0;
0031             NClu4_track_count = 0;
0032
0033             Init();
0034         };
0035
0036         void FindMegaTracks(vector<vector<pair<bool,clu_info>>>& inner_clu_phi_map, vector<vector<pair<bool,clu_info>>>& outer_clu_phi_map, pair<double,double> evt_z, int final_centrality_bin);
0037         void ClearEvt();
0038         int Get_NClu3_track_count() {return NClu3_track_count;};
0039         int Get_NClu4_track_count() {return NClu4_track_count;};
0040         double Get_performance_cut() {return performance_cut;};
0041         vector<double> Get_mega_track_eta() {return mega_track_eta;};
0042
0043         TH2F * Get_NClu3_track_centrality_2D() {return NClu3_track_centrality_2D;};
0044         TH2F * Get_NClu4_track_centrality_2D() {return NClu4_track_centrality_2D;};
0045         TH1F * Get_cluster4_track_phi_1D() {return cluster4_track_phi_1D;};
0046         TH1F * Get_cluster3_track_phi_1D() {return cluster3_track_phi_1D;};
0047         TH1F * Get_cluster3_inner_track_phi_1D() {return cluster3_inner_track_phi_1D;};
0048         TH1F * Get_cluster3_outer_track_phi_1D() {return cluster3_outer_track_phi_1D;};
0049         TH1F * Get_clu4_track_ReducedChi2_1D() {return clu4_track_ReducedChi2_1D;};
0050         TH1F * Get_clu3_track_ReducedChi2_1D() {return clu3_track_ReducedChi2_1D;};
0051         TH1F * Get_mega_track_eta_1D() {return mega_track_eta_1D;};
0052         TH1F * Get_cluster3_inner_track_eta_1D() {return cluster3_inner_track_eta_1D;};
0053         TH1F * Get_cluster3_outer_track_eta_1D() {return cluster3_outer_track_eta_1D;};
0054
0055
0056     protected:
0057
0058         int N_centrality_bin;
0059         bool mark_used;
0060
0061         TH2F * NClu3_track_centrality_2D;
0062         TH2F * NClu4_track_centrality_2D;
0063         TH1F * cluster4_track_phi_1D;
0064         TH1F * cluster3_track_phi_1D;
0065         TH1F * cluster3_inner_track_eta_1D;
0066         TH1F * cluster3_inner_track_phi_1D;
0067         TH1F * cluster3_outer_track_eta_1D;
0068         TH1F * cluster3_outer_track_phi_1D;
0069         TH1F * clu4_track_ReducedChi2_1D;
0070         TH1F * clu3_track_ReducedChi2_1D;
0071         TH1F * mega_track_eta_1D;
0072
0073         int NClu3_track_count;
0074         int NClu4_track_count;
0075
0076         double mega_cluster_deltaphi_cut = 0.1;
0077         double performance_cut = 0.0005;
0078
0079
0080         TF1 * pol0_fit;
0081         TGraph * r_phi_gr;
0082         TGraphErrors * track_gr;
0083
0084         map<string,int> inner_used_mega_clu;
0085         map<string,int> outer_used_mega_clu;
0086
0087         // note : to keep the cluster-4 tracklet fit result and the final 1D array
0088         vector<double>              clu4_track_performance;
0089         vector<vector<int>>         clu4_track_index;
0090         vector<pair<double,double>> clu4_track_clu_phi;    // note : 8 index, 4 inner, 4 outer
0091
0092         vector<double>              clu3_track_performance;
0093         vector<vector<int>>         clu3_track_index;
0094         vector<pair<double,double>> clu3_track_clu_phi;    // note : 6 index
0095         vector<int>                 clu3_track_type;
0096
0097         pair<double,double> Get_eta_pair;
0098         vector<double> mega_track_eta; // note : keep the eta of the mega tracklet (mega tracklet means: 3/4-cluster tracks)
0099
0100         void Init();
0101         double get_delta_phi(double angle_1, double angle_2);
0102         double get_track_phi(double inner_clu_phi_in, double delta_phi_in);
0103         double get_offset_clu_phi(double x_in, double y_in);
0104         double grEY_stddev(TGraphErrors * input_grr);
0105         pair<double, double> mirrorPolynomial(double a, double b);
0106         pair<double,double> Get_eta3(vector<double>p0, vector<double>p1, vector<double>p2, vector<double>p3);
0107         pair<double,double> Get_eta4(vector<double>p0, vector<double>p1, vector<double>p2, vector<double>p3, vector<double>p4);
0108         void self_check(vector<vector<pair<bool,clu_info>>> inner_clu_phi_map, vector<vector<pair<bool,clu_info>>> outer_clu_phi_map);
0109         // double get_radius(double x, double y);
0110         // pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) ;// note : inner p0, outer p1, vector {r,z}, -> {y,x}
0111
0112 };
0113
0114
0115
0116
0117 void MegaTrackFinder::Init()
0118 {
0119     cluster4_track_phi_1D = new TH1F("","cluster4_track_phi_1D", 200, -40, 400);
0120     cluster4_track_phi_1D -> GetXaxis() -> SetTitle("Tracklet #phi");
0121     cluster4_track_phi_1D -> GetYaxis() -> SetTitle("Entry");
0122     cluster4_track_phi_1D -> GetXaxis() -> SetNdivisions(505);
0123
0124     cluster3_track_phi_1D = new TH1F("","cluster3_track_phi_1D", 200, -40, 400);
0125     cluster3_track_phi_1D -> GetXaxis() -> SetTitle("Tracklet #phi");
0126     cluster3_track_phi_1D -> GetYaxis() -> SetTitle("Entry");
0127     cluster3_track_phi_1D -> GetXaxis() -> SetNdivisions(505);
0128
0129     cluster3_inner_track_phi_1D = new TH1F("","cluster3_inner_track_phi_1D", 200, -40, 400);
0130     cluster3_inner_track_phi_1D -> GetXaxis() -> SetTitle("Tracklet #phi");
0131     cluster3_inner_track_phi_1D -> GetYaxis() -> SetTitle("Entry");
0132     cluster3_inner_track_phi_1D -> GetXaxis() -> SetNdivisions(505);
0133
0134     cluster3_outer_track_phi_1D = new TH1F("","cluster3_outer_track_phi_1D", 200, -40, 400);
0135     cluster3_outer_track_phi_1D -> GetXaxis() -> SetTitle("Tracklet #phi");
0136     cluster3_outer_track_phi_1D -> GetYaxis() -> SetTitle("Entry");
0137     cluster3_outer_track_phi_1D -> GetXaxis() -> SetNdivisions(505);
0138
0139     NClu3_track_centrality_2D = new TH2F("","NClu3_track_centrality_2D", N_centrality_bin, 0, N_centrality_bin, 50, 0, 50);
0140     NClu3_track_centrality_2D -> GetXaxis() -> SetTitle("Centrality bin");
0141     NClu3_track_centrality_2D -> GetYaxis() -> SetTitle("N Clu3 tracks");
0142     NClu3_track_centrality_2D -> GetXaxis() -> SetNdivisions(505);
0143
0144     NClu4_track_centrality_2D = new TH2F("","NClu4_track_centrality_2D", N_centrality_bin, 0, N_centrality_bin, 20, 0, 20);
0145     NClu4_track_centrality_2D -> GetXaxis() -> SetTitle("Centrality bin");
0146     NClu4_track_centrality_2D -> GetYaxis() -> SetTitle("N Clu4 tracks");
0147     NClu4_track_centrality_2D -> GetXaxis() -> SetNdivisions(505);
0148
0149     clu4_track_ReducedChi2_1D = new TH1F("","clu4_track_ReducedChi2_1D", 100, 0, 0.01);
0150     clu4_track_ReducedChi2_1D -> GetXaxis() -> SetTitle("Reduced #chi^{2}");
0151     clu4_track_ReducedChi2_1D -> GetYaxis() -> SetTitle("Entry");
0152     clu4_track_ReducedChi2_1D -> GetXaxis() -> SetNdivisions(505);
0153
0154     clu3_track_ReducedChi2_1D = new TH1F("","clu3_track_ReducedChi2_1D", 100, 0, 0.01);
0155     clu3_track_ReducedChi2_1D -> GetXaxis() -> SetTitle("Reduced #chi^{2}");
0156     clu3_track_ReducedChi2_1D -> GetYaxis() -> SetTitle("Entry");
0157     clu3_track_ReducedChi2_1D -> GetXaxis() -> SetNdivisions(505);
0158
0159     double Eta_NBin = 61;
0160     double Eta_Min = -3.05;
0161     double Eta_Max = 3.05;
0162
0163     mega_track_eta_1D = new TH1F("","mega_track_eta_1D", Eta_NBin, Eta_Min, Eta_Max);
0164     mega_track_eta_1D -> GetXaxis() -> SetTitle("Mega track #eta");
0165     mega_track_eta_1D -> GetYaxis() -> SetTitle("Entry");
0166     mega_track_eta_1D -> GetXaxis() -> SetNdivisions(505);
0167
0168     cluster3_inner_track_eta_1D = new TH1F("","cluster3_inner_track_eta_1D", Eta_NBin, Eta_Min, Eta_Max);
0169     cluster3_inner_track_eta_1D -> GetXaxis() -> SetTitle("clu3 inner track #eta");
0170     cluster3_inner_track_eta_1D -> GetYaxis() -> SetTitle("Entry");
0171     cluster3_inner_track_eta_1D -> GetXaxis() -> SetNdivisions(505);
0172
0173     cluster3_outer_track_eta_1D = new TH1F("","cluster3_outer_track_eta_1D", Eta_NBin, Eta_Min, Eta_Max);
0174     cluster3_outer_track_eta_1D -> GetXaxis() -> SetTitle("clu3 outer track #eta");
0175     cluster3_outer_track_eta_1D -> GetYaxis() -> SetTitle("Entry");
0176     cluster3_outer_track_eta_1D -> GetXaxis() -> SetNdivisions(505);
0177
0178     pol0_fit = new TF1("pol0_fit","pol0",0,360);
0179     r_phi_gr = new TGraph();
0180     r_phi_gr -> Set(0);
0181
0182     track_gr = new TGraphErrors();
0183     track_gr -> Set(0);
0184 }
0185
0186 void MegaTrackFinder::FindMegaTracks(vector<vector<pair<bool,clu_info>>>& inner_clu_phi_map, vector<vector<pair<bool,clu_info>>>& outer_clu_phi_map, pair<double,double> evt_z, int final_centrality_bin)
0187 {
0188     // cout<<" "<<endl;
0189
0190     // note : find the proto "mega-inner + single outer" 3-cluster track
0191     // note : if such a 3-cluster track is found, then keep this one, and try to seek for the 4-cluster track as it should fit to tighter selection
0192     for (int inner_clu1_first = 0; inner_clu1_first < 360; inner_clu1_first++)
0193     {
0194         for (int inner_clu1_second = 0; inner_clu1_second < inner_clu_phi_map[inner_clu1_first].size(); inner_clu1_second++) // note : the first inner clu
0195         {
0196             if (inner_clu_phi_map[inner_clu1_first][inner_clu1_second].first == true) {continue;} // note : if the cluster is used, skip it
0197
0198             double inner_clu1_radius = get_radius(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.x - beam_origin.first, inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.y - beam_origin.second);
0199             double inner_clu1_phi    = get_offset_clu_phi(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.x, inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.y);
0200
0201             for (int scan_inner_i2 = -1; scan_inner_i2 < 2; scan_inner_i2++) // note : -1, 0, 1
0202             {
0203                 int inner_clu2_first = ((inner_clu1_first + scan_inner_i2) < 0) ? 360 + (inner_clu1_first + scan_inner_i2) : ((inner_clu1_first + scan_inner_i2) > 359) ? (inner_clu1_first + scan_inner_i2)-360 : inner_clu1_first + scan_inner_i2;
0204
0205                 for (int inner_clu2_second = 0; inner_clu2_second < inner_clu_phi_map[inner_clu2_first].size(); inner_clu2_second++) // note : the second inner clu
0206                 {
0207                     if (inner_clu_phi_map[inner_clu2_first][inner_clu2_second].first == true) {continue;} // note : if the cluster is used, skip it
0208
0209                     // note : the cluster itself
0210                     if (inner_clu1_first == inner_clu2_first && inner_clu1_second == inner_clu2_second) {continue;}
0211                     // note : in the same sub-layer, skip
0212                     // todo : todo : it may not work for the data
0213                     if (inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.layer == inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.layer) {continue;}
0214                     // note : I expect the two cluster have to have the same z position
0215                     // note : but in case of data, the z position of the same strip may be fluctuated a little bit
0216                     // todo : currently, set the Z position flutuation to be 4 mm
0217                     if (fabs(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.z - inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.z) > 4 ) {continue;}
0218
0219                     double inner_clu2_radius = get_radius(inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.x - beam_origin.first, inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.y - beam_origin.second);
0220                     double inner_clu2_phi    = get_offset_clu_phi(inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.x, inner_clu_phi_map[inner_clu2_first][inner_clu2_second].second.y);
0221                     // double delta_phi_in1_in2 = get_delta_phi(inner_clu1_phi, inner_clu2_phi);
0222
0223                     int inner_comb_index[4] = {inner_clu1_first, inner_clu1_second, inner_clu2_first, inner_clu2_second};
0224                     int id_inner_small_r = ( inner_clu1_radius < inner_clu2_radius ) ? 0 : 2;
0225
0226                     for (int scan_outer_i1 = -1; scan_outer_i1 < 2; scan_outer_i1++) // note : -1, 0, 1
0227                      {
0228                         int outer_clu1_first = ((inner_clu2_first + scan_outer_i1) < 0) ? 360 + (inner_clu2_first + scan_outer_i1) : ((inner_clu2_first + scan_outer_i1) > 359) ? (inner_clu2_first + scan_outer_i1)-360 : inner_clu2_first + scan_outer_i1;
0229
0230                         for (int outer_clu1_second = 0; outer_clu1_second < outer_clu_phi_map[outer_clu1_first].size(); outer_clu1_second++)
0231                         {
0232                             if (outer_clu_phi_map[outer_clu1_first][outer_clu1_second].first == true) {continue;} // note : if the cluster is used, skip it
0233
0234                             double outer_clu1_radius = get_radius(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.x - beam_origin.first, outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.y - beam_origin.second);
0235                             double outer_clu1_phi    = get_offset_clu_phi(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.x, outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.y);
0236
0237                             // double delta_phi_in1_outout= in_delta_phi(inner_clu1_phi, outer_clu1_phi);
0238                             // double delta_phi_in2_outout= in_delta_phi(inner_clu2_phi, outer_clu1_phi);
0239
0240                             pair<double,double> z_range_info = Get_possible_zvtx(
0241                                 0., // get_radius(beam_origin.first,beam_origin.second),
0242                                 {get_radius(inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.x - beam_origin.first, inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.y - beam_origin.second), inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.z}, // note : unsign radius
0243                                 {get_radius(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.x - beam_origin.first, outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.y - beam_origin.second), outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.z}  // note : unsign radius
0244                             );
0245
0246                             // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0247                             if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0248
0249                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), inner_clu1_radius, inner_clu1_phi);
0250                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), inner_clu2_radius, inner_clu2_phi);
0251                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), outer_clu1_radius, outer_clu1_phi);
0252                             r_phi_gr -> Fit(pol0_fit, "NQ");
0253
0254                             clu3_track_performance.push_back(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF())); // note: the reduced chi2 is used as the performance of the tracklet
0255                             clu3_track_index.push_back({inner_clu1_first, inner_clu1_second, inner_clu2_first, inner_clu2_second, outer_clu1_first, outer_clu1_second});
0256                             clu3_track_clu_phi.push_back({pol0_fit->GetParameter(0), (inner_clu1_phi + inner_clu2_phi + outer_clu1_phi)/3. });
0257                             clu3_track_type.push_back(0); // note : inner-mega and outer single
0258                             clu3_track_ReducedChi2_1D -> Fill(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF()));
0259
0260
0261                             r_phi_gr -> Set(0);
0262
0263                             // note : here we start to work on the 4-cluster tracklet instead of spin out the code, in order to make it more faster
0264                             // note : the reason is that, the 4-cluster tracklet should have more tight selections, so once we found a 3-cluster track, then we try to find the 4-cluster track
0265                             for (int scan_outer_i2 = -1; scan_outer_i2 < 2; scan_outer_i2++) // note : -1, 0, 1
0266                             {
0267                                 int outer_clu2_first = ((outer_clu1_first + scan_outer_i2) < 0) ? 360 + (outer_clu1_first + scan_outer_i2) : ((outer_clu1_first + scan_outer_i2) > 359) ? (outer_clu1_first + scan_outer_i2)-360 : outer_clu1_first + scan_outer_i2;
0268
0269                                 for (int outer_clu2_second = 0; outer_clu2_second < outer_clu_phi_map[outer_clu2_first].size(); outer_clu2_second++)
0270                                 {
0271                                     if (outer_clu_phi_map[outer_clu2_first][outer_clu2_second].first == true) {continue;}
0272
0273                                     // note : the cluster itself
0274                                     if (outer_clu1_first == outer_clu2_first && outer_clu1_second == outer_clu2_second) {continue;}
0275                                     // note : in the same sub-layer, skip
0276                                     // todo : todo : it may not work for the data
0277                                     if (outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.layer == outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.layer) {continue;}
0278                                     // note : I expect the two cluster have to have the same z position
0279                                     // note : but in case of data, the z position of the same strip may be fluctuated a little bit
0280                                     // todo : currently, set the Z position flutuation to be 4 mm
0281                                     if (fabs(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.z - outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.z) > 4 ) {continue;}
0282
0283                                     double outer_clu2_radius = get_radius(outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.x - beam_origin.first, outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.y - beam_origin.second);
0284                                     double outer_clu2_phi    = get_offset_clu_phi(outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.x, outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.y);
0285                                     // double delta_phi_out2_out2 = get_delta_phi(outer_clu1_phi, outer_clu2_phi);
0286
0287                                     int outer_comb_index[4] = {outer_clu1_first, outer_clu1_second, outer_clu2_first, outer_clu2_second};
0288                                     int id_outer_large_r = ( outer_clu1_radius > outer_clu2_radius ) ? 0 : 2;
0289
0290                                     pair<double,double> z_range_info = Get_possible_zvtx(
0291                                         0., // get_radius(beam_origin.first,beam_origin.second),
0292                                         {get_radius(inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.x - beam_origin.first, inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.y - beam_origin.second), inner_clu_phi_map[inner_comb_index[id_inner_small_r]][inner_comb_index[id_inner_small_r+1]].second.z}, // note : unsign radius
0293                                         {get_radius(outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.x - beam_origin.first, outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.y - beam_origin.second), outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.z}  // note : unsign radius
0294                                     );
0295
0296                                     // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0297                                     if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0298
0299                                     r_phi_gr -> SetPoint(r_phi_gr -> GetN(), inner_clu1_radius, inner_clu1_phi);
0300                                     r_phi_gr -> SetPoint(r_phi_gr -> GetN(), inner_clu2_radius, inner_clu2_phi);
0301                                     r_phi_gr -> SetPoint(r_phi_gr -> GetN(), outer_clu1_radius, outer_clu1_phi);
0302                                     r_phi_gr -> SetPoint(r_phi_gr -> GetN(), outer_clu2_radius, outer_clu2_phi);
0303                                     r_phi_gr -> Fit(pol0_fit, "NQ");
0304
0305                                     clu4_track_performance.push_back(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF())); // note: the reduced chi2 is used as the performance of the tracklet
0306                                     clu4_track_index.push_back({inner_clu1_first, inner_clu1_second, inner_clu2_first, inner_clu2_second, outer_clu1_first, outer_clu1_second, outer_clu2_first, outer_clu2_second});
0307                                     clu4_track_clu_phi.push_back({pol0_fit->GetParameter(0), (inner_clu1_phi + inner_clu2_phi + outer_clu1_phi + outer_clu2_phi)/4. });
0308                                     clu4_track_ReducedChi2_1D -> Fill(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF()));
0309
0310                                     r_phi_gr -> Set(0);
0311                                 }
0312
0313                             } // note : outer_clu2 loop
0314                         }
0315                     } // note : outer_clu1_loop
0316
0317
0318
0319
0320                 }
0321             }// note : inner clu2 loop
0322
0323         }
0324     } // note : inner clu1 loop
0325
0326
0327     // note : try to check another case of the 3-cluster tracklet, single inner and mega outer
0328     for (int outer_clu1_first = 0; outer_clu1_first < 360; outer_clu1_first++)
0329     {
0330         for (int outer_clu1_second = 0; outer_clu1_second < outer_clu_phi_map[outer_clu1_first].size(); outer_clu1_second++) // note : the first outer clu
0331         {
0332             if (outer_clu_phi_map[outer_clu1_first][outer_clu1_second].first == true) {continue;} // note : if the cluster is used, skip it
0333
0334             double outer_clu1_radius = get_radius(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.x - beam_origin.first, outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.y - beam_origin.second);
0335             double outer_clu1_phi    = get_offset_clu_phi(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.x, outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.y);
0336
0337             for (int scan_outer_i2 = -1; scan_outer_i2 < 2; scan_outer_i2++) // note : -1, 0, 1
0338             {
0339                 int outer_clu2_first = ((outer_clu1_first + scan_outer_i2) < 0) ? 360 + (outer_clu1_first + scan_outer_i2) : ((outer_clu1_first + scan_outer_i2) > 359) ? (outer_clu1_first + scan_outer_i2)-360 : outer_clu1_first + scan_outer_i2;
0340
0341                 for (int outer_clu2_second = 0; outer_clu2_second < outer_clu_phi_map[outer_clu2_first].size(); outer_clu2_second++) // note : the second outer clu
0342                 {
0343                     if (outer_clu_phi_map[outer_clu2_first][outer_clu2_second].first == true) {continue;} // note : if the cluster is used, skip it
0344
0345                     // note : the cluster itself
0346                     if (outer_clu1_first == outer_clu2_first && outer_clu1_second == outer_clu2_second) {continue;}
0347                     // note : in the same sub-layer, skip
0348                     // todo : todo : it may not work for the data
0349                     if (outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.layer == outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.layer) {continue;}
0350                     // note : I expect the two cluster have to have the same z position
0351                     // note : but in case of data, the z position of the same strip may be fluctuated a little bit
0352                     // todo : currently, set the Z position flutuation to be 4 mm
0353                     if (fabs(outer_clu_phi_map[outer_clu1_first][outer_clu1_second].second.z - outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.z) > 4 ) {continue;}
0354
0355                     double outer_clu2_radius = get_radius(outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.x - beam_origin.first, outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.y - beam_origin.second);
0356                     double outer_clu2_phi    = get_offset_clu_phi(outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.x, outer_clu_phi_map[outer_clu2_first][outer_clu2_second].second.y);
0357                     // double delta_phi_out1_out2 = get_delta_phi(outer_clu1_phi, outer_clu2_phi);
0358
0359                     int outer_comb_index[4] = {outer_clu1_first, outer_clu1_second, outer_clu2_first, outer_clu2_second};
0360                     int id_outer_large_r = ( outer_clu1_radius > outer_clu2_radius ) ? 0 : 2;
0361
0362                     for (int scan_inner_i1 = -1; scan_inner_i1 < 2; scan_inner_i1++) // note : -1, 0, 1
0363                      {
0364                         int inner_clu1_first = ((outer_clu2_first + scan_inner_i1) < 0) ? 360 + (outer_clu2_first + scan_inner_i1) : ((outer_clu2_first + scan_inner_i1) > 359) ? (outer_clu2_first + scan_inner_i1)-360 : outer_clu2_first + scan_inner_i1;
0365
0366                         for (int inner_clu1_second = 0; inner_clu1_second < inner_clu_phi_map[inner_clu1_first].size(); inner_clu1_second++)
0367                         {
0368                             if (inner_clu_phi_map[inner_clu1_first][inner_clu1_second].first == true) {continue;} // note : if the cluster is used, skip it
0369
0370                             double inner_clu1_radius = get_radius(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.x - beam_origin.first, inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.y - beam_origin.second);
0371                             double inner_clu1_phi    = get_offset_clu_phi(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.x, inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.y);
0372
0373                             // double delta_phi_out1_in1 = get_delta_phi(outer_clu1_phi, inner_clu1_phi);
0374                             // double delta_phi_out2_in1 = get_delta_phi(outer_clu2_phi, inner_clu1_phi);
0375
0376                             pair<double,double> z_range_info = Get_possible_zvtx(
0377                                 0., // get_radius(beam_origin.first,beam_origin.second),
0378                                 {get_radius(inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.x - beam_origin.first, inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.y - beam_origin.second), inner_clu_phi_map[inner_clu1_first][inner_clu1_second].second.z},  // note : unsign radius
0379                                 {get_radius(outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.x - beam_origin.first, outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.y - beam_origin.second), outer_clu_phi_map[outer_comb_index[id_outer_large_r]][outer_comb_index[id_outer_large_r+1]].second.z} // note : unsign radius
0380                             );
0381
0382                             // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0383                             if (z_range_info.first - z_range_info.second > evt_z.first + evt_z.second || z_range_info.first + z_range_info.second < evt_z.first - evt_z.second) {continue;}
0384
0385                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), outer_clu1_radius, outer_clu1_phi);
0386                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), outer_clu2_radius, outer_clu2_phi);
0387                             r_phi_gr -> SetPoint(r_phi_gr -> GetN(), inner_clu1_radius, inner_clu1_phi);
0388                             r_phi_gr -> Fit(pol0_fit, "NQ");
0389
0390                             clu3_track_performance.push_back(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF())); // note: the reduced chi2 is used as the performance of the tracklet
0391                             clu3_track_index.push_back({outer_clu1_first, outer_clu1_second, outer_clu2_first, outer_clu2_second, inner_clu1_first, inner_clu1_second});
0392                             clu3_track_clu_phi.push_back({pol0_fit->GetParameter(0), (outer_clu1_phi + outer_clu2_phi + inner_clu1_phi)/3. });
0393                             clu3_track_type.push_back(1); // note : outer-mega and inner single
0394                             clu3_track_ReducedChi2_1D -> Fill(pol0_fit->GetChisquare()/double(pol0_fit->GetNDF()));
0395
0396
0397                             r_phi_gr -> Set(0);
0398                         }
0399                     } // note : inner_clu1_loop
0400
0401
0402
0403
0404                 }
0405             }// note : outer clu2 loop
0406
0407         }
0408     } // note : outer clu1 loop
0409
0410
0411
0412     // note : try to sort all the proto 4-cluster tracklets
0413     long long N_proto_clu4_track = clu4_track_performance.size();
0414     long long ind_clu4_track[clu4_track_performance.size()];
0415     TMath::Sort(N_proto_clu4_track, &clu4_track_performance[0], ind_clu4_track, false);
0416
0417     for (int i = 0; i < clu4_track_index.size(); i++)
0418     {
0419         int better_track_index = ind_clu4_track[i];
0420
0421         if (clu4_track_performance[better_track_index] > performance_cut) {continue;}
0422         if (inner_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][0], clu4_track_index[better_track_index][1])] != 0) {continue;}
0423         if (inner_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][2], clu4_track_index[better_track_index][3])] != 0) {continue;}
0424         if (outer_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][4], clu4_track_index[better_track_index][5])] != 0) {continue;}
0425         if (outer_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][6], clu4_track_index[better_track_index][7])] != 0) {continue;}
0426
0427
0428         // note : the new trackleet is found!
0429         double inner_clu1_radius = get_radius(inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.x - beam_origin.first, inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y - beam_origin.second);
0430         double inner_clu2_radius = get_radius(inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.x - beam_origin.first, inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y - beam_origin.second);
0431         double outer_clu1_radius = get_radius(outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.x - beam_origin.first, outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y - beam_origin.second);
0432         double outer_clu2_radius = get_radius(outer_clu_phi_map[clu3_track_index[better_track_index][6]][clu3_track_index[better_track_index][7]].second.x - beam_origin.first, outer_clu_phi_map[clu3_track_index[better_track_index][6]][clu3_track_index[better_track_index][7]].second.y - beam_origin.second);
0433
0434         Get_eta_pair = Get_eta4(
0435             {0., evt_z.first, evt_z.second},
0436             {inner_clu1_radius, inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.z},
0437             {inner_clu2_radius, inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.z},
0438             {outer_clu1_radius, outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.z},
0439             {outer_clu2_radius, outer_clu_phi_map[clu3_track_index[better_track_index][6]][clu3_track_index[better_track_index][7]].second.z}
0440         );
0441
0442         mega_track_eta_1D -> Fill(Get_eta_pair.second);
0443         mega_track_eta.push_back(Get_eta_pair.second);
0444         cluster4_track_phi_1D -> Fill(clu4_track_clu_phi[better_track_index].first);
0445         NClu4_track_count += 1;
0446
0447
0448         // note : mark the clusters as used
0449         inner_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][0], clu4_track_index[better_track_index][1])] += 1;
0450         inner_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][2], clu4_track_index[better_track_index][3])] += 1;
0451         outer_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][4], clu4_track_index[better_track_index][5])] += 1;
0452         outer_used_mega_clu[Form("%i_%i", clu4_track_index[better_track_index][6], clu4_track_index[better_track_index][7])] += 1;
0453
0454         if (mark_used == false) {continue;}
0455         inner_clu_phi_map[clu4_track_index[better_track_index][0]][clu4_track_index[better_track_index][1]].first = true;
0456         inner_clu_phi_map[clu4_track_index[better_track_index][2]][clu4_track_index[better_track_index][3]].first = true;
0457         outer_clu_phi_map[clu4_track_index[better_track_index][4]][clu4_track_index[better_track_index][5]].first = true;
0458         outer_clu_phi_map[clu4_track_index[better_track_index][6]][clu4_track_index[better_track_index][7]].first = true;
0459     }
0460
0461
0462     // note : try to sort all the proto 3-cluster tracklets with mega inner cluster by the fit performance,
0463     long long N_proto_clu3_track = clu3_track_performance.size();
0464     long long ind_clu3_track[clu3_track_performance.size()];
0465     TMath::Sort(N_proto_clu3_track, &clu3_track_performance[0], ind_clu3_track, false);
0466
0467     for (int i = 0; i < clu3_track_index.size(); i++)
0468     {
0469         int better_track_index = ind_clu3_track[i];
0470         if (clu3_track_performance[better_track_index] > performance_cut) {continue;}
0471
0472         // note : inner mega and outer single
0473         if (clu3_track_type[better_track_index] == 0){
0474
0475             if (inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][0], clu3_track_index[better_track_index][1])] != 0) {continue;}
0476             if (inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][2], clu3_track_index[better_track_index][3])] != 0) {continue;}
0477             if (outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][4], clu3_track_index[better_track_index][5])] != 0) {continue;}
0478
0479             // note : the new trackleet is found!
0480             double inner_clu1_radius = get_radius(inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.x - beam_origin.first, inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y - beam_origin.second);
0481             double inner_clu2_radius = get_radius(inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.x - beam_origin.first, inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y - beam_origin.second);
0482             double outer_clu1_radius = get_radius(outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.x - beam_origin.first, outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y - beam_origin.second);
0483
0484             Get_eta_pair = Get_eta3(
0485                 {0., evt_z.first, evt_z.second},
0486                 {inner_clu1_radius, inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.z},
0487                 {inner_clu2_radius, inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.z},
0488                 {outer_clu1_radius, outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.z}
0489             );
0490
0491             mega_track_eta_1D -> Fill(Get_eta_pair.second);
0492             mega_track_eta.push_back(Get_eta_pair.second);
0493             cluster3_inner_track_eta_1D -> Fill(Get_eta_pair.second);
0494             cluster3_inner_track_phi_1D -> Fill(clu3_track_clu_phi[better_track_index].first);
0495             cluster3_track_phi_1D       -> Fill(clu3_track_clu_phi[better_track_index].first);
0496             NClu3_track_count += 1;
0497
0498             // cout<<"the clu indics : "<<clu3_track_index[better_track_index][0]<<", "<<clu3_track_index[better_track_index][1]<<", "<<clu3_track_index[better_track_index][2]<<", "<<clu3_track_index[better_track_index][3]<<", "<<clu3_track_index[better_track_index][4]<<", "<<clu3_track_index[better_track_index][5]
0499             // <<"\t trackPhi: "<<clu3_track_clu_phi[better_track_index].first<<" radius: "<< inner_clu1_radius<<", "<<inner_clu2_radius<<", "<<outer_clu1_radius<<"\t ---- clu3 proto track performance : "<<clu3_track_performance[better_track_index]<<endl;
0500
0501             // cout<<"inner posXY ("
0502             // <<inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.x<<", "<<inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y<<"), ("
0503             // <<inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.x<<", "<<inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y<<"), ("
0504             // <<outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.x<<", "<<outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y<<")"<<endl;
0505
0506             // cout<<"inner PosRZ ("<<evt_z.first<<", 0),("
0507             // <<inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.z<<", "<<inner_clu1_radius<<"), ("
0508             // <<inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.z<<", "<<inner_clu2_radius<<"), ("
0509             // <<outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.z<<", "<<outer_clu1_radius<<")"<<endl;
0510
0511             // ", "<<<inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y
0512             // ", "<<<inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y
0513             // ", "<<<outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y
0514             cout<<"("<<clu3_track_clu_phi[better_track_index].first<<", "<<Get_eta_pair.second<<"), ";
0515
0516             // note : mark the clusters as used
0517             inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][0], clu3_track_index[better_track_index][1])] += 1;
0518             inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][2], clu3_track_index[better_track_index][3])] += 1;
0519             outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][4], clu3_track_index[better_track_index][5])] += 1;
0520
0521             if (mark_used == false) {continue;}
0522             inner_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].first = true;
0523             inner_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].first = true;
0524             outer_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].first = true;
0525         }
0526         else if (clu3_track_type[better_track_index] == 1)
0527         {
0528             if (outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][0], clu3_track_index[better_track_index][1])] != 0) {continue;}
0529             if (outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][2], clu3_track_index[better_track_index][3])] != 0) {continue;}
0530             if (inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][4], clu3_track_index[better_track_index][5])] != 0) {continue;}
0531
0532             // note : the new trackleet is found!
0533             double inner_clu1_radius = get_radius(inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.x - beam_origin.first, inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y - beam_origin.second);
0534             double outer_clu1_radius = get_radius(outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.x - beam_origin.first, outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y - beam_origin.second);
0535             double outer_clu2_radius = get_radius(outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.x - beam_origin.first, outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y - beam_origin.second);
0536
0537             Get_eta_pair = Get_eta3(
0538                 {0., evt_z.first, evt_z.second},
0539                 {inner_clu1_radius, inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.z},
0540                 {outer_clu1_radius, outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.z},
0541                 {outer_clu2_radius, outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.z}
0542             );
0543
0544             mega_track_eta_1D -> Fill(Get_eta_pair.second);
0545             mega_track_eta.push_back(Get_eta_pair.second);
0546             cluster3_outer_track_eta_1D -> Fill(Get_eta_pair.second);
0547             cluster3_outer_track_phi_1D -> Fill(clu3_track_clu_phi[better_track_index].first);
0548             cluster3_track_phi_1D       -> Fill(clu3_track_clu_phi[better_track_index].first);
0549             NClu3_track_count += 1;
0550
0551             // cout<<"the clu indics : "<<clu3_track_index[better_track_index][0]<<", "<<clu3_track_index[better_track_index][1]<<", "<<clu3_track_index[better_track_index][2]<<", "<<clu3_track_index[better_track_index][3]<<", "<<clu3_track_index[better_track_index][4]<<", "<<clu3_track_index[better_track_index][5]
0552             // <<"\t trackPhi: "<<clu3_track_clu_phi[better_track_index].first<<" radius: "<< outer_clu1_radius<<", "<<outer_clu2_radius<<", "<<inner_clu1_radius<<"\t ---- clu3 proto track performance : "<<clu3_track_performance[better_track_index]<<endl;
0553
0554             // cout<<"outer posXY ("
0555             // <<outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.x<<", "<<outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y<<"), ("
0556             // <<outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.x<<", "<<outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y<<"), ("
0557             // <<inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.x<<", "<<inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y<<")"<<endl;
0558
0559             // cout<<"outer PosRZ ("<<evt_z.first<<", 0),("
0560             // <<outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.z<<", "<<-1 * outer_clu1_radius<<"), ("
0561             // <<outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.z<<", "<<-1 * outer_clu2_radius<<"), ("
0562             // <<inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.z<<", "<<-1 * inner_clu1_radius<<")"<<endl;
0563
0564             // ", "<<<outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].second.y
0565             // ", "<<<outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].second.y
0566             // ", "<<<inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].second.y
0567
0568             cout<<"("<<clu3_track_clu_phi[better_track_index].first<<", "<<Get_eta_pair.second<<"), ";
0569
0570             // note : mark the clusters as used
0571             outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][0], clu3_track_index[better_track_index][1])] += 1;
0572             outer_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][2], clu3_track_index[better_track_index][3])] += 1;
0573             inner_used_mega_clu[Form("%i_%i", clu3_track_index[better_track_index][4], clu3_track_index[better_track_index][5])] += 1;
0574
0575             if (mark_used == false) {continue;}
0576             outer_clu_phi_map[clu3_track_index[better_track_index][0]][clu3_track_index[better_track_index][1]].first = true;
0577             outer_clu_phi_map[clu3_track_index[better_track_index][2]][clu3_track_index[better_track_index][3]].first = true;
0578             inner_clu_phi_map[clu3_track_index[better_track_index][4]][clu3_track_index[better_track_index][5]].first = true;
0579         }
0580
0581     }
0582
0583
0584     // note : to check discrepancy between the number of mega tracks found and number of clusters marked as used
0585     // self_check(inner_clu_phi_map, outer_clu_phi_map);
0586
0587     NClu4_track_centrality_2D -> Fill(final_centrality_bin, NClu4_track_count);
0588     NClu3_track_centrality_2D -> Fill(final_centrality_bin, NClu3_track_count);
0589 }
0590
0591 void MegaTrackFinder::ClearEvt()
0592 {
0593     inner_used_mega_clu.clear();
0594     outer_used_mega_clu.clear();
0595
0596     clu4_track_performance.clear();
0597     clu4_track_index.clear();
0598     clu4_track_clu_phi.clear();
0599
0600     clu3_track_performance.clear();
0601     clu3_track_index.clear();
0602     clu3_track_clu_phi.clear();
0603     clu3_track_type.clear();
0604
0605     mega_track_eta.clear();
0606
0607     NClu3_track_count = 0;
0608     NClu4_track_count = 0;
0609
0610 }
0611
0612 // note : angle_1 = inner clu phi
0613 // note: angle_2 = outer clu phi
0614 double MegaTrackFinder::get_delta_phi(double angle_1, double angle_2)
0615 {
0616     vector<double> vec_abs = {fabs(angle_1 - angle_2), fabs(angle_1 - angle_2 + 360), fabs(angle_1 - angle_2 - 360)};
0617     vector<double> vec = {(angle_1 - angle_2), (angle_1 - angle_2 + 360), (angle_1 - angle_2 - 360)};
0618     return vec[std::distance(vec_abs.begin(), std::min_element(vec_abs.begin(),vec_abs.end()))];
0619 }
0620
0621 double MegaTrackFinder::get_track_phi(double inner_clu_phi_in, double delta_phi_in)
0622 {
0623     double track_phi = inner_clu_phi_in - (delta_phi_in/2.);
0624     if (track_phi < 0) {track_phi += 360;}
0625     else if (track_phi > 360) {track_phi -= 360;}
0626     else if (track_phi == 360) {track_phi = 0;}
0627     else {track_phi = track_phi;}
0628     return track_phi;
0629 }
0630
0631 double MegaTrackFinder::get_offset_clu_phi(double x_in, double y_in)
0632 {
0633     // if (beam_origin.second == 0) {cout<<"wrong beam position"<<endl; exit(0);}
0634     return (y_in - beam_origin.second < 0) ? atan2(y_in - beam_origin.second, x_in - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(y_in - beam_origin.second, x_in - beam_origin.first) * (180./TMath::Pi());
0635 }
0636
0637 pair<double, double> MegaTrackFinder::mirrorPolynomial(double a, double b) {
0638     // Interchange 'x' and 'y'
0639     double mirroredA = 1.0 / a;
0640     double mirroredB = -b / a;
0641
0642     return {mirroredA, mirroredB};
0643 }
0644
0645 double MegaTrackFinder::grEY_stddev(TGraphErrors * input_grr)
0646 {
0647     vector<double> input_vector; input_vector.clear();
0648     for (int i = 0; i < input_grr -> GetN(); i++) {input_vector.push_back( input_grr -> GetPointY(i) );}
0649
0650     double average = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0651
0652     double sum_subt = 0;
0653     for (int ele : input_vector) {sum_subt+=pow((ele-average),2);}
0654
0655     return sqrt(sum_subt/(input_vector.size()-1));
0656 }
0657
0658 pair<double,double> MegaTrackFinder::Get_eta3(vector<double>p0, vector<double>p1, vector<double>p2, vector<double>p3)
0659 {
0660     // if (track_gr -> GetN() != 0) {cout<<"In INTTEta class, track_gr is not empty, track_gr size : "<<track_gr -> GetN()<<endl; exit(0);}
0661
0662     // if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
0663
0664     track_gr -> Set(0);
0665
0666     vector<double> r_vec  = {p0[0], p1[0], p2[0], p3[0]};
0667     vector<double> z_vec  = {p0[1], p1[1], p2[1], p3[1]};
0668     vector<double> re_vec = {0, 0, 0, 0};
0669     vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10., ( fabs( p2[1] ) < 130 ) ? 8. : 10., ( fabs( p3[1] ) < 130 ) ? 8. : 10.};
0670
0671     // note : swap the r and z, to have easier fitting
0672     // note : in principle, Z should be in the x axis, R should be in the Y axis
0673     for (int i = 0; i < r_vec.size(); i++)
0674     {
0675         track_gr -> SetPoint(i,r_vec[i],z_vec[i]);
0676         track_gr -> SetPointError(i,re_vec[i],ze_vec[i]);
0677
0678         // cout<<"In INTTEta class, point : "<<i<<" r : "<<r_vec[i]<<" z : "<<z_vec[i]<<" re : "<<re_vec[i]<<" ze : "<<ze_vec[i]<<endl;
0679     }
0680
0681     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
0682
0683     if (vertical_line) {return {0,0};}
0684     else
0685     {
0686         fit_rz -> SetParameters(0,0);
0687
0688         track_gr -> Fit(fit_rz,"NQ");
0689
0690         pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
0691
0692         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
0693     }
0694
0695 }
0696
0697
0698 pair<double,double> MegaTrackFinder::Get_eta4(vector<double>p0, vector<double>p1, vector<double>p2, vector<double>p3, vector<double>p4)
0699 {
0700     // if (track_gr -> GetN() != 0) {cout<<"In INTTEta class, track_gr is not empty, track_gr size : "<<track_gr -> GetN()<<endl; exit(0);}
0701
0702     // if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
0703     track_gr -> Set(0);
0704
0705     vector<double> r_vec  = {p0[0], p1[0], p2[0], p3[0], p4[0]};
0706     vector<double> z_vec  = {p0[1], p1[1], p2[1], p3[1], p4[1]};
0707     vector<double> re_vec = {0, 0, 0, 0, 0};
0708     vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10., ( fabs( p2[1] ) < 130 ) ? 8. : 10., ( fabs( p3[1] ) < 130 ) ? 8. : 10., ( fabs( p4[1] ) < 130 ) ? 8. : 10.};
0709
0710     // note : swap the r and z, to have easier fitting
0711     // note : in principle, Z should be in the x axis, R should be in the Y axis
0712     for (int i = 0; i < r_vec.size(); i++)
0713     {
0714         track_gr -> SetPoint(i,r_vec[i],z_vec[i]);
0715         track_gr -> SetPointError(i,re_vec[i],ze_vec[i]);
0716
0717         // cout<<"In INTTEta class, point : "<<i<<" r : "<<r_vec[i]<<" z : "<<z_vec[i]<<" re : "<<re_vec[i]<<" ze : "<<ze_vec[i]<<endl;
0718     }
0719
0720     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
0721
0722     if (vertical_line) {return {0,0};}
0723     else
0724     {
0725         fit_rz -> SetParameters(0,0);
0726
0727         track_gr -> Fit(fit_rz,"NQ");
0728
0729         pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
0730
0731         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
0732     }
0733
0734 }
0735
0736 void MegaTrackFinder::self_check(vector<vector<pair<bool,clu_info>>> inner_clu_phi_map, vector<vector<pair<bool,clu_info>>> outer_clu_phi_map)
0737 {
0738     int N_used_clu_inner_check = 0;
0739     int N_used_clu_outer_check = 0;
0740
0741     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0742     {
0743         // note : N cluster in this phi cell
0744         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
0745         {
0746             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {N_used_clu_inner_check += 1;}
0747         }
0748     }
0749
0750     for (int outer_phi_i = 0; outer_phi_i < 360; outer_phi_i++) // note : each phi cell (1 degree)
0751     {
0752         // note : N cluster in this phi cell
0753         for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[outer_phi_i].size(); outer_phi_clu_i++)
0754         {
0755             if (outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].first == true) {N_used_clu_outer_check += 1;}
0756         }
0757     }
0758     if (NClu4_track_count * 4 + NClu3_track_count * 3 != (N_used_clu_inner_check + N_used_clu_outer_check))
0759     {
0760         cout<<"---- N 4-cluster tracklet: "<<NClu4_track_count<<" N 3-cluster tracklet: "<<NClu3_track_count<<endl;
0761         cout<<"---- N used inner clu : "<<N_used_clu_inner_check<<" N used outer clu : "<<N_used_clu_outer_check<<endl;
0762         cout<<"---- diff : "<<(N_used_clu_inner_check + N_used_clu_outer_check) - (NClu4_track_count * 4 + NClu3_track_count * 3)<<endl;
0763     }
0764 }
0765
0766 // double MegaTrackFinder::get_radius(double x, double y)
0767 // {
0768 //     return sqrt(pow(x,2)+pow(y,2));
0769 // }
0770
0771 // pair<double,double> MegaTrackFinder::Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) // note : inner p0, outer p1, vector {r,z}, -> {y,x}
0772 // {
0773 //     vector<double> p0_z_edge = { ( fabs( p0[1] ) < 130 ) ? p0[1] - 8. : p0[1] - 10., ( fabs( p0[1] ) < 130 ) ? p0[1] + 8. : p0[1] + 10.}; // note : vector {left edge, right edge}
0774 //     vector<double> p1_z_edge = { ( fabs( p1[1] ) < 130 ) ? p1[1] - 8. : p1[1] - 10., ( fabs( p1[1] ) < 130 ) ? p1[1] + 8. : p1[1] + 10.}; // note : vector {left edge, right edge}
0775
0776 //     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
0777 //     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
0778
0779 //     double mid_point = (edge_first + edge_second) / 2.;
0780 //     double possible_width = fabs(edge_first - edge_second) / 2.;
0781
0782 //     return {mid_point, possible_width}; // note : first : mid point, second : width
0783
0784 // }
0785
0786 #endif