DST_MC/old_code/INTT_xyvtx.C

0001 // #include "DAC_Scan_ladder.h"
0002 //#include "InttConversion.h"
0003 #include "/sphenix/user/ChengWei/INTT/INTT_commissioning/INTT_CW/INTT_commissioning/DAC_Scan/InttConversion_new.h"
0004 #include "/sphenix/user/ChengWei/INTT/INTT_commissioning/INTT_CW/INTT_commissioning/DAC_Scan/InttClustering.h"
0005 #include "../sigmaEff.h"
0006 #include "../sPhenixStyle.C"
0007 #include "../INTTDSTchain.C"
0008
0009 // todo : the number of number is given by the adc_setting_run !!!
0010 // todo : also the range of the hist.
0011 // todo : the adc follows the following convention
0012 // todo : 1. the increment has to be 4
0013 // todo : remember to check the "adc_conv"
0014 // vector<vector<int>> adc_setting_run = {
0015 //     // {8  , 12 , 16 , 20 , 24 , 28 , 32 , 36 },
0016 //     // {28 , 32 , 36 , 40 , 44 , 48 , 52 , 56 },
0017 //     {48 , 52 , 56 , 60 , 64 , 68 , 72 , 76 }, // note : 3
0018 //     {68 , 72 , 76 , 80 , 84 , 88 , 92 , 96 }, // note : 4
0019 //     {88 , 92 , 96 , 100, 104, 108, 112, 116}, // note : 5
0020 //     {108, 112, 116, 120, 124, 128, 132, 136}, // note : 6
0021 //     {128, 132, 136, 140, 144, 148, 152, 156}, // note : 7
0022 //     // {148, 152, 156, 160, 164, 168, 172, 176},
0023 //     // {168, 172, 176, 180, 184, 188, 192, 196},
0024 //     // {188, 192, 196, 200, 204, 208, 212, 216}
0025 // };
0026
0027 vector<vector<int>> adc_setting_run = {
0028     {15, 30, 60, 90, 120, 150, 180, 210, 240}
0029     // {15, 30, 50, 70, 90, 110, 130, 150,170}
0030     // {8  , 12 , 16 , 20 , 24 , 28 , 32 , 36 },
0031     // {28 , 32 , 36 , 40 , 44 , 48 , 52 , 56 },
0032     // {48 , 52 , 56 , 60 , 64 , 68 , 72 , 76 }, // note : 3
0033     // {68 , 72 , 76 , 80 , 84 , 88 , 92 , 96 }, // note : 4
0034     // {88 , 92 , 96 , 100, 104, 108, 112, 116}, // note : 5
0035     // {108, 112, 116, 120, 124, 128, 132, 136}, // note : 6
0036     // {128, 132, 136, 140, 144, 148, 152, 156}, // note : 7
0037     // {148, 152, 156, 160, 164, 168, 172, 176},
0038     // {168, 172, 176, 180, 184, 188, 192, 196},
0039     // {188, 192, 196, 200, 204, 208, 212, 216}
0040 };
0041
0042 TString color_code_2[8] = {"#CC768D","#19768D","#DDA573","#009193","#6E9193","#941100","#A08144","#517E66"};
0043
0044 struct full_hit_info {
0045     int FC;
0046     int chip_id;
0047     int chan_id;
0048     int adc;
0049 };
0050
0051
0052 struct ladder_info {
0053     int FC;
0054     TString Port;
0055     int ROC;
0056     int Direction; // note : 0 : south, 1 : north
0057 };
0058
0059 double get_radius(double x, double y)
0060 {
0061     return sqrt(pow(x,2)+pow(y,2));
0062 }
0063
0064 double get_radius_sign(double x, double y)
0065 {
0066     double phi = ((y) < 0) ? atan2((y),(x)) * (180./TMath::Pi()) + 360 : atan2((y),(x)) * (180./TMath::Pi());
0067
0068     return (phi > 180) ? sqrt(pow(x,2)+pow(y,2)) * -1 : sqrt(pow(x,2)+pow(y,2));
0069 }
0070
0071 double sin_func(double *x, double *par)
0072 {
0073     return par[0] * sin(par[1] * x[0] + par[2]) + par[3];
0074 }
0075
0076 double cos_func(double *x, double *par)
0077 {
0078     return -1 * par[0] * cos(par[1] * (x[0] + par[2])) + par[3];
0079 }
0080
0081 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)
0082 {
0083     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
0084     pad -> SetLeftMargin   (left);
0085     pad -> SetRightMargin  (right);
0086     pad -> SetTopMargin    (top);
0087     pad -> SetBottomMargin (bottom);
0088     pad -> SetTicks(1,1);
0089     if (set_logY == true)
0090     {
0091         pad -> SetLogy (1);
0092     }
0093
0094 }
0095
0096 std::vector<double> calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y) {
0097
0098     if (x1 != x2)
0099     {
0100         // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
0101         double a = (y2 - y1) / (x2 - x1);
0102         double b = y1 - a * x1;
0103
0104         // cout<<"slope : y="<<a<<"x+"<<b<<endl;
0105
0106         // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
0107         double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
0108
0109         // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
0110         double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
0111         double Yc = a * Xc + b;
0112
0113         return { closest_distance, Xc, Yc };
0114     }
0115     else
0116     {
0117         double closest_distance = std::abs(x1 - target_x);
0118         double Xc = x1;
0119         double Yc = target_y;
0120
0121         return { closest_distance, Xc, Yc };
0122     }
0123
0124
0125 }
0126
0127 double get_z_vertex(clu_info inner_clu, clu_info outer_clu, double target_x, double target_y)
0128 {
0129     // note : x = z,
0130     // note : y = radius
0131     double inner_clu_r = sqrt(pow(inner_clu.x,2)+ pow(inner_clu.y,2));
0132     double outer_clu_r = sqrt(pow(outer_clu.x,2)+ pow(outer_clu.y,2));
0133     double target_r    = sqrt(pow(target_x,2)   + pow(target_y,2));
0134
0135     // Use the slope equation (y = ax + b) to calculate the x-coordinate for the target y
0136     if ( fabs(outer_clu.z - inner_clu.z) < 0.00001 ){
0137         return outer_clu.z;
0138     }
0139     else {
0140         double slope = (outer_clu_r - inner_clu_r) / (outer_clu.z - inner_clu.z);
0141         double yIntercept = inner_clu_r - slope * inner_clu.z;
0142         double xCoordinate = (target_r - yIntercept) / slope;
0143         return xCoordinate;
0144     }
0145
0146 }
0147
0148 // Function to calculate the angle between two vectors in degrees using the cross product
0149 double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
0150     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
0151     double vector1X = x2 - x1;
0152     double vector1Y = y2 - y1;
0153
0154     double vector2X = targetX - x1;
0155     double vector2Y = targetY - y1;
0156
0157     // Calculate the cross product of vector_1 and vector_2 (z-component)
0158     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
0159
0160     // cout<<" crossProduct : "<<crossProduct<<endl;
0161
0162     // Calculate the magnitudes of vector_1 and vector_2
0163     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
0164     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
0165
0166     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
0167     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
0168
0169     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
0170     // Convert the angle from radians to degrees and return it
0171     double angleInDegrees = angleInRadians * 180.0 / M_PI;
0172
0173     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
0174     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
0175
0176     // cout<<"angle : "<<angleInDegrees_new<<endl;
0177
0178     double DCA_distance = sin(angleInRadians_new) * magnitude2;
0179
0180     return DCA_distance;
0181 }
0182
0183 void temp_bkg(TPad * c1, string conversion_mode, double peek, pair<double,double> beam_origin, InttConversion * ch_pos_DB)
0184 {
0185     c1 -> cd();
0186
0187     int N_ladder[4] = {12, 12, 16, 16};
0188     string ladder_index_string[16] = {"00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15"};
0189
0190     vector<double> x_vec; x_vec.clear();
0191     vector<double> y_vec; y_vec.clear();
0192
0193     vector<double> x_vec_2; x_vec_2.clear();
0194     vector<double> y_vec_2; y_vec_2.clear();
0195
0196     TGraph * bkg = new TGraph();
0197     bkg -> SetTitle("INTT event display X-Y plane");
0198     bkg -> SetMarkerStyle(20);
0199     bkg -> SetMarkerSize(0.1);
0200     bkg -> SetPoint(0,0,0);
0201     bkg -> SetPoint(1,beam_origin.first,beam_origin.second);
0202     bkg -> GetXaxis() -> SetLimits(-150,150);
0203     bkg -> GetYaxis() -> SetRangeUser(-150,150);
0204     bkg -> GetXaxis() -> SetTitle("X [mm]");
0205     bkg -> GetYaxis() -> SetTitle("Y [mm]");
0206
0207     bkg -> Draw("ap");
0208
0209     TLine * ladder_line = new TLine();
0210     ladder_line -> SetLineWidth(1);
0211
0212     for (int server_i = 0; server_i < 4; server_i++)
0213     {
0214         for (int FC_i = 0; FC_i < 14; FC_i++)
0215         {
0216             ladder_line -> DrawLine(
0217                 ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,14,0).x, ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,14,0).y,
0218                 ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,1,0).x, ch_pos_DB -> Get_XY_all(Form("intt%i",server_i),FC_i,1,0).y
0219             );
0220         }
0221     }
0222
0223     ladder_line -> Draw("l same");
0224
0225 }
0226
0227 double grEY_stddev(TGraphErrors * input_grr)
0228 {
0229     vector<double> input_vector; input_vector.clear();
0230     for (int i = 0; i < input_grr -> GetN(); i++)
0231     {
0232         input_vector.push_back( input_grr -> GetPointY(i) );
0233     }
0234
0235     double sum=0;
0236     double average;
0237     double sum_subt = 0;
0238     for (int i=0; i<input_vector.size(); i++)
0239         {
0240             sum+=input_vector[i];
0241
0242         }
0243     average=sum/input_vector.size();
0244     //cout<<"average is : "<<average<<endl;
0245
0246     for (int i=0; i<input_vector.size(); i++)
0247         {
0248             //cout<<input_vector[i]-average<<endl;
0249             sum_subt+=pow((input_vector[i]-average),2);
0250
0251         }
0252     //cout<<"sum_subt : "<<sum_subt<<endl;
0253     double stddev;
0254     stddev=sqrt(sum_subt/(input_vector.size()-1));
0255     return stddev;
0256 }
0257
0258 pair<double, double> mirrorPolynomial(double a, double b) {
0259     // Interchange 'x' and 'y'
0260     double mirroredA = 1.0 / a;
0261     double mirroredB = -b / a;
0262
0263     return {mirroredA, mirroredB};
0264 }
0265
0266 // note : pair<reduced chi2, eta of the track>
0267 // note : vector : {r,z}
0268 // note : p0 : vertex point {r,z,zerror}
0269 // note : p1 : inner layer
0270 // note : p2 : outer layer
0271 pair<double,double> Get_eta(vector<double>p0, vector<double>p1, vector<double>p2)
0272 {
0273     vector<double> r_vec  = {p0[0], p1[0], p2[0]};
0274     vector<double> z_vec  = {p0[1], p1[1], p2[1]};
0275     vector<double> re_vec = {0, 0, 0};
0276     vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10., ( fabs( p2[1] ) < 130 ) ? 8. : 10.};
0277
0278     // note : swap the r and z, to have easier fitting
0279     // note : in principle, Z should be in the x axis, R should be in the Y axis
0280     TGraphErrors * track_gr = new TGraphErrors(3,&r_vec[0],&z_vec[0],&re_vec[0],&ze_vec[0]);
0281
0282     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
0283
0284     if (vertical_line) {return {0,0};}
0285     else
0286     {
0287         TF1 * fit_rz = new TF1("fit_rz","pol1",-500,500);
0288         fit_rz -> SetParameters(0,0);
0289
0290         track_gr -> Fit(fit_rz,"NQ");
0291
0292         pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
0293
0294         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
0295
0296     }
0297
0298 }
0299
0300 // note : vector : {r,z}
0301 // note : p0 : vertex point {r,z,zerror}
0302 // note : p1 : another point from detector
0303 // note : since only two points -> no strip width considered
0304 double Get_eta_2P(vector<double>p0, vector<double>p1)
0305 {
0306     // vector<double> r_vec  = {p0[0], p1[0]};
0307     // vector<double> z_vec  = {p0[1], p1[1]};
0308     // vector<double> re_vec = {0, 0};
0309     // vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10.};
0310
0311     // note : swap the r and z, to have easier fitting
0312     // note : in principle, Z should be in the x axis, R should be in the Y axis
0313     // TGraphErrors * track_gr = new TGraphErrors(3,&r_vec[0],&z_vec[0],&re_vec[0],&ze_vec[0]);
0314
0315     // double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
0316
0317     return  -1 * TMath::Log( fabs( tan( atan2(p1[0] - p0[0], p1[1] - p0[1]) / 2 ) ) );
0318
0319     // if (vertical_line) {return 0;}
0320     // else
0321     // {
0322     //     TF1 * fit_rz = new TF1("fit_rz","pol1",-500,500);
0323     //     fit_rz -> SetParameters(0,0);
0324
0325     //     track_gr -> Fit(fit_rz,"NQ");
0326
0327     //     pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
0328     // }
0329
0330 }
0331
0332 double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
0333 {
0334     if ( fabs(p0x - p1x) < 0.00001 ){ // note : the line is vertical (if z is along the x axis)
0335         return p0x;
0336     }
0337     else {
0338         double slope = (p1y - p0y) / (p1x - p0x);
0339         double yIntercept = p0y - slope * p0x;
0340         double xCoordinate = (given_y - yIntercept) / slope;
0341         return xCoordinate;
0342     }
0343 }
0344
0345 pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) // note : inner p0, outer p1, vector {r,z}, -> {y,x}
0346 {
0347     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}
0348     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}
0349
0350     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
0351     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
0352
0353     double mid_point = (edge_first + edge_second) / 2.;
0354     double possible_width = fabs(edge_first - edge_second) / 2.;
0355
0356     return {mid_point, possible_width}; // note : first : mid point, second : width
0357
0358 }
0359
0360 void TH2F_threshold(TH2F * input_hist, double threshold)
0361 {
0362     for (int xi = 0; xi < input_hist -> GetNbinsX(); xi++){
0363         for(int yi = 0; yi < input_hist -> GetNbinsY(); yi++){
0364             if (input_hist -> GetBinContent(xi + 1, yi +1) < threshold){ input_hist -> SetBinContent(xi + 1, yi +1, 0); }
0365         }
0366     }
0367 }
0368
0369 // note : rotate the ClusX and ClusY
0370 pair<double,double> rotatePoint(double x, double y) {
0371     // Convert the rotation angle from degrees to radians
0372     double rotation = 180.; // todo rotation is here
0373     double angleRad = rotation * M_PI / 180.0;
0374
0375     // Perform the rotation
0376     double xOut = x * cos(angleRad) - y * sin(angleRad);
0377     double yOut = x * sin(angleRad) + y * cos(angleRad);
0378
0379     return {xOut,yOut};
0380 }
0381
0382 // note : use "ls *.root > file_list.txt" to create the list of the file in the folder, full directory in the file_list.txt
0383 // note : set_folder_name = "folder_xxxx"
0384 // note : server_name = "inttx"
0385 void INTT_xyvtx(/*bool full_event, int run_Nevent*/)
0386 {
0387
0388     SetsPhenixStyle();
0389
0390     TCanvas * c4 = new TCanvas("","",1200,800);
0391     c4 -> cd();
0392
0393     TCanvas * c2 = new TCanvas("","",2500,800);
0394     c2 -> cd();
0395     TPad *pad_xy = new TPad(Form("pad_xy"), "", 0.0, 0.0, 0.33, 1.0);
0396     Characterize_Pad(pad_xy, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0397     pad_xy -> Draw();
0398
0399     TPad *pad_rz = new TPad(Form("pad_rz"), "", 0.33, 0.0, 0.66, 1.0);
0400     Characterize_Pad(pad_rz, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0401     pad_rz -> Draw();
0402
0403     TPad *pad_z = new TPad(Form("pad_z"), "", 0.66, 0.0, 1.0, 1.0);
0404     Characterize_Pad(pad_z, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0405     pad_z -> Draw();
0406
0407     TCanvas * c1 = new TCanvas("","",950,800);
0408     c1 -> cd();
0409
0410     //todo : can be more than it, possibly
0411     vector<string> conversion_mode_BD = {"ideal", "survey_1_XYAlpha_Peek_3.32mm", "full_survey_3.32"};
0412     // vector<string> conversion_mode_BD = {"ideal", "survey_1_XYAlpha_Peek", "full_survey_3.32"};
0413
0414     // string mother_folder_directory = "/home/phnxrc/INTT/cwshih/DACscan_data/zero_magnet_Takashi_used";
0415     // string file_name = "beam_inttall-00020869-0000_event_base_ana_cluster_ideal_excludeR1500_100kEvent";
0416     // string file_name = "beam_inttall-00020869-0000_event_base_ana_cluster_survey_1_XYAlpha_Peek_3.32mm_excludeR1500_100kEvent";
0417
0418     // string run_ID = "20869"; string file_event = "150";
0419     // todo : change the file name.
0420     // string mother_folder_directory = "/sphenix/user/ChengWei/INTT/INTT_commissioning/ZeroField/" + run_ID;
0421     // string file_name = "beam_inttall-000" +run_ID+ "-0000_event_base_ana_cluster_survey_1_XYAlpha_Peek_3.32mm_excludeR20000_"+ file_event +"kEvent_3HotCut";
0422     // string file_name = "beam_inttall-000" +run_ID+ "-0000_event_base_ana_cluster_"+ conversion_mode_BD[geo_mode_id] +"_excludeR1500_"+ file_event +"kEvent_3HotCut";
0423
0424     int geo_mode_id = 0;
0425     string mother_folder_directory = "/sphenix/user/ChengWei/sPH_dNdeta/dNdEta_INTT_MC_388_000"; string file_name = "MC_ZF_all";
0426     string MC_list_name = "dst_INTT_dNdEta_388_000.list";
0427
0428     TChain * chain_in = new TChain("EventTree");
0429     INTTDSTchain inttDSTMC(chain_in,mother_folder_directory,MC_list_name);
0430     std::printf("inttDSTMC N event : %lli \n", chain_in->GetEntries());
0431     long long N_event = chain_in->GetEntries();
0432
0433     cout<<"the input file name : "<<file_name<<endl;
0434     sleep(5);
0435
0436     // string mother_folder_directory = "/home/phnxrc/INTT/cwshih/DACscan_data/new_DAC_Scan_0722/AllServer/DAC2";
0437     // string file_name = "beam_inttall-00023058-0000_event_base_ana_cluster_ideal_excludeR2000_100kEvent";
0438
0439
0440
0441     // pair<double,double> beam_origin = {-0.457 + 0.0276, 2.657 - 0.2814}; // note : for run20869
0442     // pair<double,double> beam_origin_U = {0.0004, 0.0035}; // note : for run20869, the uncertainty of the VTXxy, calculated by the fit errors with the error propagation
0443     // pair<double,double> beam_origin = {-0.015, 0.012}; // note : for the MC run /sphenix/user/ChengWei/sPH_dNdeta/dNdEta_INTT_MC
0444
0445     // pair<double,double> beam_origin = {-0.0, -2};
0446     pair<double,double> beam_origin = {-0.4, 2.4};
0447     // pair<double,double> beam_origin = {-0.0, -0.0};
0448
0449     // pair<double,double> beam_origin = {1., -4.};
0450     double temp_Y_align = 0.;
0451     double temp_X_align = -0.;
0452
0453     double zvtx_hist_l = -500;
0454     double zvtx_hist_r = 500;
0455
0456     // todo : MC doesn't keep Nhits information for the moment
0457     // int Nhit_cut = 1500;           // note : if (> Nhit_cut)          -> continue
0458     // int Nhit_cutl = 100;          // note : if (< Nhit_cutl)         -> continue
0459     int clu_size_cut = 4;           // note : if (> clu_size_cut)      -> continue
0460     double clu_sum_adc_cut = -1;    // note : if (< clu_sum_adc_cut)   -> continue
0461     int N_clu_cut = 500;          // note : if (> N_clu_cut)         -> continue  unit number
0462     int N_clu_cutl = 20;           // note : if (< N_clu_cutl)        -> continue  unit number
0463     double phi_diff_cut = 5.72;      // note : if (< phi_diff_cut)      -> pass      unit degree
0464     double DCA_cut = 4;             // note : if (< DCA_cut)           -> pass      unit mm
0465     int zvtx_cal_require = 15;      // note : if (> zvtx_cal_require)  -> pass
0466     // int zvtx_draw_requireL = 15;
0467     // int zvtx_draw_requireR = 40000;   // note : if ( zvtx_draw_requireL < event && event < zvtx_draw_requireR) -> pass
0468     // double Integrate_portion = 0.4; // todo : it was 0.6826, try to require less event, as most of the tracklets are not that useful
0469     // double Integrate_portion_final = 0.68;
0470     bool draw_event_display = false;
0471     int print_rate = 5;
0472
0473     bool rotate_inner_phi = false;
0474
0475     string output_directory = (rotate_inner_phi == true) ? Form("%s/folder_%s_advanced_inner_phi_rotate",mother_folder_directory.c_str(),file_name.c_str()) : Form("%s/folder_%s_advanced",mother_folder_directory.c_str(),file_name.c_str());
0476     system(Form("mkdir %s",output_directory.c_str()));
0477
0478     // todo : change the geo draw mode if needed
0479     int geo_mode_id_draw = 0;
0480     string conversion_mode = (geo_mode_id_draw == 0) ? "ideal" : "survey_1_XYAlpha_Peek";
0481     double peek = 3.32405;
0482
0483     // note : the initiator of the INTT geometry class
0484     InttConversion * ch_pos_DB = new InttConversion(conversion_mode_BD[geo_mode_id], peek);
0485
0486     TLatex *draw_text = new TLatex();
0487     draw_text -> SetNDC();
0488     draw_text -> SetTextSize(0.03);
0489
0490     vector<clu_info> temp_sPH_inner_nocolumn_vec; temp_sPH_inner_nocolumn_vec.clear();
0491     vector<clu_info> temp_sPH_outer_nocolumn_vec; temp_sPH_outer_nocolumn_vec.clear();
0492     vector<vector<double>> temp_sPH_nocolumn_vec(2);
0493     vector<vector<double>> temp_sPH_nocolumn_rz_vec(2);
0494
0495     TH2F * angle_correlation = new TH2F("","angle_correlation",361,0,361,361,0,361);
0496     angle_correlation -> SetStats(0);
0497     angle_correlation -> GetXaxis() -> SetTitle("Inner Phi [degree]");
0498     angle_correlation -> GetYaxis() -> SetTitle("Outer Phi [degree]");
0499
0500     TH2F * angle_diff_DCA_dist = new TH2F("","angle_diff_DCA_dist",100,-1.5,1.5,100,-3.5,3.5);
0501     angle_diff_DCA_dist -> SetStats(0);
0502     angle_diff_DCA_dist -> GetXaxis() -> SetTitle("Inner - Outer [degree]");
0503     angle_diff_DCA_dist -> GetYaxis() -> SetTitle("DCA distance [mm]");
0504
0505     TH1F * angle_diff = new TH1F("angle_diff","angle_diff",100,0,5);
0506     angle_diff -> SetStats(0);
0507     angle_diff -> GetXaxis() -> SetTitle("|Inner - Outer| [degree]");
0508     angle_diff -> GetYaxis() -> SetTitle("Entry");
0509
0510     TH2F * angle_diff_inner_phi = new TH2F("","angle_diff_inner_phi",361,0,361,100,-1.5,1.5);
0511     angle_diff_inner_phi -> SetStats(0);
0512     angle_diff_inner_phi -> GetXaxis() -> SetTitle("Inner phi [degree]");
0513     angle_diff_inner_phi -> GetYaxis() -> SetTitle("Inner - Outer [degree]");
0514
0515     TH2F * angle_diff_outer_phi = new TH2F("","angle_diff_outer_phi",361,0,361,100,-1.5,1.5);
0516     angle_diff_outer_phi -> SetStats(0);
0517     angle_diff_outer_phi -> GetXaxis() -> SetTitle("Outer phi [degree]");
0518     angle_diff_outer_phi -> GetYaxis() -> SetTitle("Inner - Outer [degree]");
0519
0520     TH2F * inner_pos_xy = new TH2F("","inner_pos_xy",360,-100,100,360,-100,100);
0521     inner_pos_xy -> SetStats(0);
0522     inner_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0523     inner_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0524
0525     TH2F * outer_pos_xy = new TH2F("","outer_pos_xy",360,-150,150,360,-150,150);
0526     outer_pos_xy -> SetStats(0);
0527     outer_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0528     outer_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0529
0530     TH2F * inner_outer_pos_xy = new TH2F("","inner_outer_pos_xy",360,-150,150,360,-150,150);
0531     inner_outer_pos_xy -> SetStats(0);
0532     inner_outer_pos_xy -> GetXaxis() -> SetTitle("X axis [mm]");
0533     inner_outer_pos_xy -> GetYaxis() -> SetTitle("Y axis [mm]");
0534
0535     TH1F * z_pos_diff = new TH1F("","z_pos_diff",360,-150,150);
0536     z_pos_diff -> GetXaxis() -> SetTitle("Inner zpos - Outer zpos [mm]");
0537     z_pos_diff -> GetYaxis() -> SetTitle("Eentry");
0538
0539     TH2F * z_pos_diff_angle_diff = new TH2F("","z_pos_diff_angle_diff",100,-25,25,200,-11,11);
0540     z_pos_diff_angle_diff -> SetStats(0);
0541     z_pos_diff_angle_diff -> GetXaxis() -> SetTitle("Inner zpos - Outer zpos [mm]");
0542     z_pos_diff_angle_diff -> GetYaxis() -> SetTitle("Inner phi - Outer phi [degree]");
0543
0544     TH1F * Nhits_good = new TH1F("","Nhits_good",360,0,1000);
0545     Nhits_good -> GetXaxis() -> SetTitle("N hits in one event");
0546     Nhits_good -> GetYaxis() -> SetTitle("Eentry");
0547
0548     TH1F * z_pos_inner = new TH1F("","z_pos_inner",200,-150,150);
0549     z_pos_inner -> GetXaxis() -> SetTitle("Inner zpos [mm]");
0550     z_pos_inner -> GetYaxis() -> SetTitle("Eentry");
0551
0552     TH1F * z_pos_outer = new TH1F("","z_pos_outer",200,-150,150);
0553     z_pos_outer -> GetXaxis() -> SetTitle("Outer zpos [mm]");
0554     z_pos_outer -> GetYaxis() -> SetTitle("Eentry");
0555
0556     TH2F * z_pos_inner_outer = new TH2F("","z_pos_inner_outer",100,-150,150, 100,-150,150);
0557     z_pos_inner_outer -> SetStats(0);
0558     z_pos_inner_outer -> GetXaxis() -> SetTitle("Inner zpos [mm]");
0559     z_pos_inner_outer -> GetYaxis() -> SetTitle("Outer pos [mm]");
0560
0561     TH2F * DCA_point = new TH2F("","DCA_point",100,-10,10,100,-10,10);
0562     DCA_point -> SetStats(0);
0563     DCA_point -> GetXaxis() -> SetTitle("X pos [mm]");
0564     DCA_point -> GetYaxis() -> SetTitle("Y pos [mm]");
0565
0566     TH2F * DCA_distance_inner_phi = new TH2F("","DCA_distance_inner_phi",100,0,360,100,-10,10);
0567     DCA_distance_inner_phi -> SetStats(0);
0568     DCA_distance_inner_phi -> GetXaxis() -> SetTitle("Inner phi [degree]");
0569     DCA_distance_inner_phi -> GetYaxis() -> SetTitle("DCA [mm]");
0570
0571     TH2F * DCA_distance_outer_phi = new TH2F("","DCA_distance_outer_phi",100,0,360,100,-10,10);
0572     DCA_distance_outer_phi -> SetStats(0);
0573     DCA_distance_outer_phi -> GetXaxis() -> SetTitle("Outer phi [degree]");
0574     DCA_distance_outer_phi -> GetYaxis() -> SetTitle("DCA [mm]");
0575
0576     TH1F * N_cluster_outer_pass = new TH1F("","N_cluster_outer_pass",100,0,100);
0577     N_cluster_outer_pass -> GetXaxis() -> SetTitle("N_cluster");
0578     N_cluster_outer_pass -> GetYaxis() -> SetTitle("Eentry");
0579
0580     TH1F * N_cluster_inner_pass = new TH1F("","N_cluster_inner_pass",100,0,100);
0581     N_cluster_inner_pass -> GetXaxis() -> SetTitle("N_cluster");
0582     N_cluster_inner_pass -> GetYaxis() -> SetTitle("Eentry");
0583
0584     TH2F * N_cluster_correlation = new TH2F("","N_cluster_correlation",100,0,500,100,0,500);
0585     N_cluster_correlation -> SetStats(0);
0586     N_cluster_correlation -> GetXaxis() -> SetTitle("inner N_cluster");
0587     N_cluster_correlation -> GetYaxis() -> SetTitle("Outer N_cluster");
0588
0589     TH1F * temp_event_zvtx = new TH1F("","Z vertex dist",125,zvtx_hist_l,zvtx_hist_r);
0590     temp_event_zvtx -> GetXaxis() -> SetTitle("Z vertex position [mm]");
0591     temp_event_zvtx -> GetYaxis() -> SetTitle("Entry");
0592     vector<float> temp_event_zvtx_vec; temp_event_zvtx_vec.clear();
0593     vector<float> temp_event_zvtx_info; temp_event_zvtx_info.clear();
0594     TLine * effi_sig_range_line = new TLine();
0595     effi_sig_range_line -> SetLineWidth(3);
0596     effi_sig_range_line -> SetLineColor(TColor::GetColor("#A08144"));
0597     effi_sig_range_line -> SetLineStyle(2);
0598     TF1 * zvtx_fitting = new TF1("","gaus",-500,500);
0599
0600     double N_good_event = 0;
0601
0602     TH2F * Good_inner_outer_pos_xy_nzB = new TH2F("","Good_inner_outer_pos_xy_nzB",360,-150,150,360,-150,150);
0603     Good_inner_outer_pos_xy_nzB -> SetStats(0);
0604     Good_inner_outer_pos_xy_nzB -> GetXaxis() -> SetTitle("X axis");
0605     Good_inner_outer_pos_xy_nzB -> GetYaxis() -> SetTitle("Y axis");
0606
0607     TH2F * Good_inner_outer_pos_xy_nzA = new TH2F("","Good_inner_outer_pos_xy_nzA",360,-150,150,360,-150,150);
0608     Good_inner_outer_pos_xy_nzA -> SetStats(0);
0609     Good_inner_outer_pos_xy_nzA -> GetXaxis() -> SetTitle("X axis");
0610     Good_inner_outer_pos_xy_nzA -> GetYaxis() -> SetTitle("Y axis");
0611
0612     TH2F * Good_inner_outer_pos_xy_pzA = new TH2F("","Good_inner_outer_pos_xy_pzA",360,-150,150,360,-150,150);
0613     Good_inner_outer_pos_xy_pzA -> SetStats(0);
0614     Good_inner_outer_pos_xy_pzA -> GetXaxis() -> SetTitle("X axis");
0615     Good_inner_outer_pos_xy_pzA -> GetYaxis() -> SetTitle("Y axis");
0616
0617     TH2F * Good_inner_outer_pos_xy_pzB = new TH2F("","Good_inner_outer_pos_xy_pzB",360,-150,150,360,-150,150);
0618     Good_inner_outer_pos_xy_pzB -> SetStats(0);
0619     Good_inner_outer_pos_xy_pzB -> GetXaxis() -> SetTitle("X axis");
0620     Good_inner_outer_pos_xy_pzB -> GetYaxis() -> SetTitle("Y axis");
0621
0622     TH2F * Good_track_space = new TH2F("","Good_track_space",200,-300,300,200,-300,300);
0623     Good_track_space -> SetStats(0);
0624     Good_track_space -> GetXaxis() -> SetTitle("X axis");
0625     Good_track_space -> GetYaxis() -> SetTitle("Y axis");
0626
0627     // dNdeta_hist -> GetXaxis() -> SetLimits(-10,10);
0628     // dNdeta_hist -> GetXaxis() -> SetNdivisions  (505);
0629
0630     TF1 * zvtx_finder = new TF1("zvtx_finder","pol0",-1,3000);
0631
0632
0633     vector<vector<double>> good_track_xy_vec; good_track_xy_vec.clear();
0634     vector<vector<double>> good_track_rz_vec; good_track_rz_vec.clear();
0635     vector<vector<double>> good_comb_rz_vec; good_comb_rz_vec.clear();
0636     vector<vector<double>> good_comb_xy_vec; good_comb_xy_vec.clear();
0637     vector<vector<double>> good_comb_phi_vec; good_comb_phi_vec.clear();
0638     vector<vector<double>> good_tracklet_rz; good_tracklet_rz.clear();
0639     TLine * track_line = new TLine();
0640     track_line -> SetLineWidth(1);
0641     track_line -> SetLineColorAlpha(38,0.5);
0642
0643     TLine * coord_line = new TLine();
0644     coord_line -> SetLineWidth(1);
0645     coord_line -> SetLineColor(16);
0646     coord_line -> SetLineStyle(2);
0647
0648
0649     vector<float> avg_event_zvtx_vec; avg_event_zvtx_vec.clear();
0650     TH1F * avg_event_zvtx = new TH1F("","avg_event_zvtx",125,zvtx_hist_l,zvtx_hist_r);
0651     // avg_event_zvtx -> SetMarkerStyle(20);
0652     // avg_event_zvtx -> SetMarkerSize(0.8);
0653     // avg_event_zvtx -> SetMarkerColor(TColor::GetColor("#1A3947"));
0654     avg_event_zvtx -> SetLineColor(TColor::GetColor("#1A3947"));
0655     avg_event_zvtx -> SetLineWidth(2);
0656     avg_event_zvtx -> GetYaxis() -> SetTitle("Entry");
0657     avg_event_zvtx -> GetXaxis() -> SetTitle("Z vertex position [mm]");
0658     avg_event_zvtx -> GetYaxis() -> SetRangeUser(0,50);
0659     avg_event_zvtx -> SetTitleSize(0.06, "X");
0660     avg_event_zvtx -> SetTitleSize(0.06, "Y");
0661     avg_event_zvtx -> GetXaxis() -> SetTitleOffset (0.82);
0662     avg_event_zvtx -> GetYaxis() -> SetTitleOffset (1.1);
0663     avg_event_zvtx -> GetXaxis() -> CenterTitle(true);
0664     avg_event_zvtx -> GetYaxis() -> CenterTitle(true);
0665     avg_event_zvtx -> GetXaxis() -> SetNdivisions  (505);
0666
0667     TH1F * zvtx_evt_width = new TH1F("","zvtx_evt_width",150,0,800);
0668     zvtx_evt_width -> GetXaxis() -> SetTitle(" mm ");
0669     zvtx_evt_width -> GetYaxis() -> SetTitle("entry");
0670
0671     TH2F * zvtx_evt_fitError_corre = new TH2F("","zvtx_evt_fitError_corre",200,0,10000,200,0,20);
0672     zvtx_evt_fitError_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0673     zvtx_evt_fitError_corre -> GetYaxis() -> SetTitle(" #pm mm ");
0674
0675     TH2F * zvtx_evt_width_corre = new TH2F("","zvtx_evt_width_corre",200,0,10000,200,0,300);
0676     zvtx_evt_width_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0677     zvtx_evt_width_corre -> GetYaxis() -> SetTitle(" mm ");
0678
0679     TH2F * zvtx_evt_nclu_corre = new TH2F("","zvtx_evt_nclu_corre",200,0,10000,200,-500,500);
0680     zvtx_evt_nclu_corre -> GetXaxis() -> SetTitle(" # of clusters ");
0681     zvtx_evt_nclu_corre -> GetYaxis() -> SetTitle(" zvtx [mm] ");
0682     zvtx_evt_nclu_corre -> GetXaxis() -> SetNdivisions  (505);
0683
0684     TH1F * width_density = new TH1F("","width_density",200,0,1); // note : N good hits / width
0685     width_density -> GetXaxis() -> SetTitle(" N good zvtx / width ");
0686     width_density -> GetYaxis() -> SetTitle(" Entry ");
0687
0688     int inner_1_check = 0; int inner_2_check = 0; int inner_3_check = 0; int inner_4_check = 0;
0689     int outer_1_check = 0; int outer_2_check = 0; int outer_3_check = 0; int outer_4_check = 0;
0690     vector<int> used_outer_check(temp_sPH_outer_nocolumn_vec.size(),0);
0691     vector<float> N_comb; vector<float> N_comb_e; vector<float> z_mid; vector<float> z_range;
0692     vector<double> effi_N_comb; vector<double> effi_z_mid; vector<double> effi_N_comb_e; vector<double> effi_z_range;
0693
0694     int N_event_pass_number = 0;
0695     int greatest_N_clu = 0;
0696
0697     vector<double> angle_diff_vec; angle_diff_vec.clear();
0698     TTree * tree_out = new TTree("tree","after ana");
0699     tree_out -> Branch("angle_diff",&angle_diff_vec);
0700
0701     if (draw_event_display) c2 -> Print(Form("%s/temp_event_display.pdf(",output_directory.c_str()));
0702
0703     for (int event_i = 0; event_i < 20000; event_i++)
0704     {
0705         if (event_i % 1000 == 0) cout<<"code running... "<<event_i<<endl;
0706
0707         // cout<<"========================= test : "<<event_i<<" ========================="<<endl;
0708
0709         // tree -> GetEntry(event_i);
0710         inttDSTMC.LoadTree(event_i);
0711         inttDSTMC.GetEntry(event_i);
0712         unsigned int length = inttDSTMC.ClusX->size();
0713
0714         // if (event_i == 13) cout<<"test, eID : "<<event_i<<" Nhits "<<N_hits<<endl;
0715
0716         // todo : MC doesn't keep Nhits information for the moment
0717         // if (N_hits > Nhit_cut) {continue;}
0718         // if (N_hits < Nhit_cutl) {continue;}
0719
0720         N_event_pass_number += 1;
0721
0722         // if (N_cluster_inner == 0 || N_cluster_outer == 0) {continue;}
0723         // if (N_cluster_inner == -1 || N_cluster_outer == -1) {continue;}
0724         if ((length) < zvtx_cal_require) {continue;}
0725         if (inttDSTMC.TruthPV_trig_z * 10. < -100) {continue;}
0726         if (inttDSTMC.NTruthVtx != 1) {continue;}
0727
0728         // cout<<"test point 1"<<endl;
0729         // note : apply some selection to remove the hot channels
0730         // note : and make the inner_clu_vec and outer_clu_vec
0731         for (int clu_i = 0; clu_i < length; clu_i++)
0732         {
0733             if (int(inttDSTMC.ClusPhiSize -> at(clu_i)) > clu_size_cut) continue;
0734             if (int(inttDSTMC.ClusAdc -> at(clu_i)) < clu_sum_adc_cut) continue;
0735
0736             int    clu_layer = (inttDSTMC.ClusLayer -> at(clu_i) == 3 || inttDSTMC.ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0737             double clu_x, clu_y;
0738             if (rotate_inner_phi == true && clu_layer == 0){ // note : inner layer only
0739                 clu_x = rotatePoint(inttDSTMC.ClusX -> at(clu_i) * 10, inttDSTMC.ClusY -> at(clu_i) * 10).first;
0740                 clu_y = rotatePoint(inttDSTMC.ClusX -> at(clu_i) * 10, inttDSTMC.ClusY -> at(clu_i) * 10).second;
0741             }
0742             else {
0743                 clu_x = inttDSTMC.ClusX -> at(clu_i) * 10; // note : change the unit from cm to mm
0744                 clu_y = inttDSTMC.ClusY -> at(clu_i) * 10;
0745             }
0746
0747             double clu_z = inttDSTMC.ClusZ -> at(clu_i) * 10;
0748             // double clu_phi = (clu_y - beam_origin.second < 0) ? atan2(clu_y - beam_origin.second,clu_x + beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(clu_y - beam_origin.second,clu_x + beam_origin.first) * (180./TMath::Pi());
0749             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0750             double clu_radius = get_radius(clu_x, clu_y);
0751
0752             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0753             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0754
0755             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0756             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0757
0758
0759             if (clu_layer == 0) // note : inner
0760                 temp_sPH_inner_nocolumn_vec.push_back({
0761                     -1,
0762                     -1,
0763                     int(inttDSTMC.ClusAdc -> at(clu_i)),
0764                     int(inttDSTMC.ClusAdc -> at(clu_i)),
0765                     int(inttDSTMC.ClusPhiSize -> at(clu_i)),
0766                     clu_x,
0767                     clu_y,
0768                     clu_z,
0769                     clu_layer,
0770                     clu_phi
0771                 });
0772
0773             if (clu_layer == 1) // note : outer
0774                 temp_sPH_outer_nocolumn_vec.push_back({
0775                     -1,
0776                     -1,
0777                     int(inttDSTMC.ClusAdc -> at(clu_i)),
0778                     int(inttDSTMC.ClusAdc -> at(clu_i)),
0779                     int(inttDSTMC.ClusPhiSize -> at(clu_i)),
0780                     clu_x,
0781                     clu_y,
0782                     clu_z,
0783                     clu_layer,
0784                     clu_phi
0785                 });
0786         }
0787         // cout<<"test point 2"<<endl;
0788
0789         inner_1_check = 0;
0790         inner_2_check = 0;
0791         inner_3_check = 0;
0792         inner_4_check = 0;
0793         for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
0794             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 0 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 90)    inner_1_check = 1;
0795             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 90 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 180)  inner_2_check = 1;
0796             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 180 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 270) inner_3_check = 1;
0797             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 270 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 360) inner_4_check = 1;
0798
0799             if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check) == 4 ) break;
0800         }
0801
0802         outer_1_check = 0;
0803         outer_2_check = 0;
0804         outer_3_check = 0;
0805         outer_4_check = 0;
0806         for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
0807             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 0 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 90)    outer_1_check = 1;
0808             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 90 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 180)  outer_2_check = 1;
0809             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 180 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 270) outer_3_check = 1;
0810             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 270 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 360) outer_4_check = 1;
0811
0812             if ( (outer_1_check + outer_2_check + outer_3_check + outer_4_check) == 4 ) break;
0813         }
0814
0815         // cout<<"test point 3"<<endl;
0816
0817         N_cluster_outer_pass -> Fill(temp_sPH_outer_nocolumn_vec.size());
0818         N_cluster_inner_pass -> Fill(temp_sPH_inner_nocolumn_vec.size());
0819         N_cluster_correlation -> Fill( temp_sPH_inner_nocolumn_vec.size(), temp_sPH_outer_nocolumn_vec.size() );
0820
0821         if ( (temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size()) > greatest_N_clu) {greatest_N_clu = (temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size());}
0822
0823         if (temp_sPH_inner_nocolumn_vec.size() < 10 || temp_sPH_outer_nocolumn_vec.size() < 10 || (temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size()) > N_clu_cut || (temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size()) < N_clu_cutl)
0824         {
0825             temp_event_zvtx_info = {-1000,-1000,-1000};
0826             temp_event_zvtx_vec.clear();
0827             temp_event_zvtx -> Reset("ICESM");
0828             good_track_xy_vec.clear();
0829             good_track_rz_vec.clear();
0830             good_comb_rz_vec.clear();
0831             good_comb_xy_vec.clear();
0832             good_comb_phi_vec.clear();
0833             temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
0834             temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
0835             temp_sPH_inner_nocolumn_vec.clear();
0836             temp_sPH_outer_nocolumn_vec.clear();
0837             angle_diff_vec.clear();
0838             continue;
0839         }
0840
0841         if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check + outer_1_check + outer_2_check + outer_3_check + outer_4_check) != 8 )
0842         {
0843             cout<<"some quater of INTT doens't work !! eID : "<<event_i<<endl;
0844             temp_event_zvtx_info = {-1000,-1000,-1000};
0845             temp_event_zvtx_vec.clear();
0846             temp_event_zvtx -> Reset("ICESM");
0847             good_track_xy_vec.clear();
0848             good_track_rz_vec.clear();
0849             good_comb_rz_vec.clear();
0850             good_comb_xy_vec.clear();
0851             good_comb_phi_vec.clear();
0852             temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
0853             temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
0854             temp_sPH_inner_nocolumn_vec.clear();
0855             temp_sPH_outer_nocolumn_vec.clear();
0856             angle_diff_vec.clear();
0857             continue;
0858         }
0859
0860         N_comb.clear();
0861         N_comb_e.clear();
0862         // z_mid.clear();
0863         // z_range.clear();
0864
0865
0866         // note : try to make sure that the clusters not to be used twice or more
0867         // used_outer_check.clear();
0868         // used_outer_check = vector<int>(temp_sPH_outer_nocolumn_vec.size(),0);
0869
0870         // cout<<"test point 4 di loop start"<<endl;
0871
0872         for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
0873         {
0874
0875             for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
0876             {
0877                 // bool DCA_tag = false;
0878                 // if (used_outer_check[outer_i] == 1) continue; // note : this outer cluster was already used, skip the trial of this combination
0879
0880                 vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0881                     temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0882                     temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0883                     beam_origin.first, beam_origin.second
0884                 );
0885
0886                 double DCA_sign = calculateAngleBetweenVectors(
0887                     temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0888                     temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0889                     beam_origin.first, beam_origin.second
0890                 );
0891
0892                 angle_correlation -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi,temp_sPH_outer_nocolumn_vec[outer_i].phi);
0893                 angle_diff -> Fill(fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi));
0894                 angle_diff_vec.push_back(fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi));
0895
0896                 if (DCA_info_vec[0] != fabs(DCA_sign) && fabs( DCA_info_vec[0] - fabs(DCA_sign) ) > 0.1 ){
0897                     cout<<"different DCA : "<<DCA_info_vec[0]<<" "<<DCA_sign<<" diff : "<<DCA_info_vec[0] - fabs(DCA_sign)<<endl;}
0898
0899                 if (fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi) < phi_diff_cut)
0900                 {
0901                     // if (DCA_info_vec[0] < DCA_cut){
0902
0903                     //     used_outer_check[outer_i] = 1; //note : this outer cluster was already used!
0904
0905                     //     pair<double,double> z_range_info = Get_possible_zvtx(
0906                     //         get_radius(beam_origin.first,beam_origin.second),
0907                     //         {get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y), temp_sPH_inner_nocolumn_vec[inner_i].z}, // note : unsign radius
0908                     //         {get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y), temp_sPH_outer_nocolumn_vec[outer_i].z}  // note : unsign radius
0909                     //     );
0910
0911                     //     N_comb.push_back(inner_i);
0912                     //     N_comb_e.push_back(0);
0913                     //     z_mid.push_back(z_range_info.first);
0914                     //     z_range.push_back(z_range_info.second);
0915
0916                     //     // DCA_tag = true;
0917                     // }
0918
0919                     angle_diff_inner_phi -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi);
0920                     angle_diff_outer_phi -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].phi, temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi);
0921
0922                     if (temp_sPH_inner_nocolumn_vec[inner_i].phi > 270)
0923                         angle_diff_DCA_dist -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi, DCA_sign);
0924
0925                     DCA_point -> Fill( DCA_info_vec[1], DCA_info_vec[2] );
0926                     // angle_correlation -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi,temp_sPH_outer_nocolumn_vec[outer_i].phi);
0927                     z_pos_diff -> Fill( temp_sPH_inner_nocolumn_vec[inner_i].z - temp_sPH_outer_nocolumn_vec[outer_i].z );
0928                     z_pos_diff_angle_diff -> Fill( temp_sPH_inner_nocolumn_vec[inner_i].z - temp_sPH_outer_nocolumn_vec[outer_i].z, temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi );
0929                     // Nhits_good -> Fill(N_hits);
0930                     z_pos_inner -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].z);
0931                     z_pos_outer -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].z);
0932                     z_pos_inner_outer -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].z, temp_sPH_outer_nocolumn_vec[outer_i].z);
0933                     // DCA_distance_inner_phi -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, (temp_sPH_inner_nocolumn_vec[inner_i].phi > 90 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 270) ? DCA_sign * -1 : DCA_sign );
0934                     // DCA_distance_outer_phi -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].phi, (temp_sPH_outer_nocolumn_vec[outer_i].phi > 90 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 270) ? DCA_sign * -1 : DCA_sign );
0935                     DCA_distance_inner_phi -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].phi, DCA_sign );
0936                     DCA_distance_outer_phi -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].phi, DCA_sign );
0937
0938                     // if(DCA_tag == true) break; // note : since this combination (one inner cluster, one outer cluster) satisfied the reuqiremet, no reason to ask this inner cluster try with other outer clusters
0939
0940                     // cout<<"good comb : "<<fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi)<<" radius in : "<<get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y)<<" radius out : "<<get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y)<<endl;
0941                 } // note : end of if
0942
0943
0944             } // note : end of outer loop
0945         } // note : end of inner loop
0946
0947         if ( draw_event_display == true && N_comb.size() > zvtx_cal_require)
0948         {
0949             TGraph * temp_event_xy = new TGraph(temp_sPH_nocolumn_vec[0].size(),&temp_sPH_nocolumn_vec[0][0],&temp_sPH_nocolumn_vec[1][0]);
0950             temp_event_xy -> SetTitle("INTT event display X-Y plane");
0951             temp_event_xy -> GetXaxis() -> SetLimits(-150,150);
0952             temp_event_xy -> GetYaxis() -> SetRangeUser(-150,150);
0953             temp_event_xy -> GetXaxis() -> SetTitle("X [mm]");
0954             temp_event_xy -> GetYaxis() -> SetTitle("Y [mm]");
0955             temp_event_xy -> SetMarkerStyle(20);
0956             temp_event_xy -> SetMarkerColor(2);
0957             temp_event_xy -> SetMarkerSize(1);
0958             TGraph * temp_event_rz = new TGraph(temp_sPH_nocolumn_rz_vec[0].size(),&temp_sPH_nocolumn_rz_vec[0][0],&temp_sPH_nocolumn_rz_vec[1][0]);
0959             temp_event_rz -> SetTitle("INTT event display r-Z plane");
0960             temp_event_rz -> GetXaxis() -> SetLimits(-500,500);
0961             temp_event_rz -> GetYaxis() -> SetRangeUser(-150,150);
0962             temp_event_rz -> GetXaxis() -> SetTitle("Z [mm]");
0963             temp_event_rz -> GetYaxis() -> SetTitle("Radius [mm]");
0964             temp_event_rz -> SetMarkerStyle(20);
0965             temp_event_rz -> SetMarkerColor(2);
0966             temp_event_rz -> SetMarkerSize(1);
0967
0968             pad_xy -> cd();
0969             temp_bkg(pad_xy, conversion_mode, peek, beam_origin, ch_pos_DB);
0970             temp_event_xy -> Draw("p same");
0971             for (int track_i = 0; track_i < good_track_xy_vec.size(); track_i++){
0972                 track_line -> DrawLine(good_track_xy_vec[track_i][0],good_track_xy_vec[track_i][1],good_track_xy_vec[track_i][2],good_track_xy_vec[track_i][3]);
0973             }
0974             track_line -> Draw("l same");
0975             draw_text -> DrawLatex(0.2, 0.85, Form("eID : %i, innter Ncluster : %i, outer Ncluster : %i",event_i, temp_sPH_inner_nocolumn_vec.size(), temp_sPH_outer_nocolumn_vec.size()));
0976
0977             pad_rz -> cd();
0978             temp_event_rz -> Draw("ap");
0979             // effi_sig_range_line -> DrawLine(temp_event_zvtx_info[0],-150,temp_event_zvtx_info[0],150);
0980             coord_line -> DrawLine(0,-150,0,150);
0981             coord_line -> DrawLine(-500,0,500,0);
0982             for (int track_i = 0; track_i < good_track_rz_vec.size(); track_i++){
0983                 track_line -> DrawLine(good_track_rz_vec[track_i][0],good_track_rz_vec[track_i][1],good_track_rz_vec[track_i][2],good_track_rz_vec[track_i][3]);
0984             }
0985             draw_text -> DrawLatex(0.2, 0.85, Form("Negative radius : Clu_{outer} > 180^{0}"));
0986             // draw_text -> DrawLatex(0.2, 0.81, Form("EffiSig avg : %.2f mm",temp_event_zvtx_info[0]));
0987
0988             // cout<<"test tag 2-5"<<endl;
0989             pad_z -> cd();
0990
0991             if(draw_event_display && (event_i % print_rate) == 0){c2 -> Print(Form("%s/temp_event_display.pdf",output_directory.c_str()));}
0992             pad_xy -> Clear();
0993             pad_rz -> Clear();
0994             pad_z  -> Clear();
0995
0996             temp_event_xy -> Delete();
0997             temp_event_rz -> Delete();
0998             // z_range_gr_draw -> Delete();
0999
1000         }
1001         // cout<<"test tag 3"<<endl;
1002
1003         for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
1004         {
1005             inner_pos_xy -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y);
1006             inner_outer_pos_xy -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y);
1007         }
1008
1009         for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
1010         {
1011             outer_pos_xy -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y);
1012             inner_outer_pos_xy -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y);
1013         }
1014
1015         // cout<<"test point 7"<<endl;
1016
1017
1018         if (angle_diff_vec.size() > 0) tree_out -> Fill();
1019
1020         angle_diff_vec.clear();
1021         temp_event_zvtx_info = {-1000,-1000,-1000};
1022         temp_event_zvtx_vec.clear();
1023         temp_event_zvtx -> Reset("ICESM");
1024         good_track_xy_vec.clear();
1025         good_track_rz_vec.clear();
1026         good_comb_rz_vec.clear();
1027         good_comb_xy_vec.clear();
1028         good_comb_phi_vec.clear();
1029         temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
1030         temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
1031         temp_sPH_inner_nocolumn_vec.clear();
1032         temp_sPH_outer_nocolumn_vec.clear();
1033         N_comb.clear();
1034         N_comb_e.clear();
1035         z_mid.clear();
1036         z_range.clear();
1037     } // note : end of event
1038
1039     if (draw_event_display) {c2 -> Print(Form("%s/temp_event_display.pdf)",output_directory.c_str()));}
1040     c2 -> Clear();
1041     c1 -> Clear();
1042
1043
1044     c1 -> cd();
1045
1046     TLatex *ltx = new TLatex();
1047     ltx->SetNDC();
1048     ltx->SetTextSize(0.045);
1049     ltx->SetTextAlign(31);
1050
1051     N_cluster_inner_pass -> Draw("hist");
1052     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1053     c1 -> Print(Form("%s/N_cluster_inner_pass.pdf",output_directory.c_str()));
1054     c1 -> Clear();
1055
1056     N_cluster_outer_pass -> Draw("hist");
1057     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1058     c1 -> Print(Form("%s/N_cluster_outer_pass.pdf",output_directory.c_str()));
1059     c1 -> Clear();
1060
1061     N_cluster_correlation -> Draw("colz0");
1062     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1063     c1 -> Print(Form("%s/N_cluster_correlation.pdf",output_directory.c_str()));
1064     c1 -> Clear();
1065
1066     DCA_point -> Draw("colz0");
1067     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1068     c1 -> Print(Form("%s/DCA_point_X%.3fY%.3f_.pdf",output_directory.c_str(),beam_origin.first,beam_origin.second));
1069     c1 -> Clear();
1070
1071     DCA_distance_inner_phi -> Draw("colz0");
1072     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1073     c1 -> Print(Form("%s/DCA_distance_inner_phi_X%.3fY%.3f_.pdf",output_directory.c_str(),beam_origin.first,beam_origin.second));
1074     c1 -> Clear();
1075
1076
1077     TH2F * DCA_distance_inner_phi_peak = (TH2F*)DCA_distance_inner_phi -> Clone();
1078     TF1 * cos_fit = new TF1("cos_fit",cos_func,0,360,4);
1079     cos_fit -> SetParNames("[A]", "[B]", "[C]", "[D]");
1080     cos_fit -> SetParameters(6.5, 0.015,  -185, 0.3);
1081     cos_fit -> SetLineColor(2);
1082     TH2F_threshold(DCA_distance_inner_phi_peak,600); // todo : change the threshold
1083     TProfile * DCA_distance_inner_phi_peak_profile =  DCA_distance_inner_phi_peak->ProfileX();
1084     DCA_distance_inner_phi_peak_profile -> Fit(cos_fit,"N","",50,250);
1085
1086
1087     DCA_distance_inner_phi_peak -> Draw("colz0");
1088     DCA_distance_inner_phi_peak_profile -> Draw("same");
1089     cos_fit -> Draw("l same");
1090     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1091     c1 -> Print(Form("%s/DCA_distance_inner_phi_peak_X%.3fY%.3f_.pdf",output_directory.c_str(),beam_origin.first,beam_origin.second));
1092     c1 -> Clear();
1093
1094
1095     DCA_distance_outer_phi -> Draw("colz0");
1096     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1097     c1 -> Print(Form("%s/DCA_distance_outer_phi_X%.3fY%.3f_.pdf",output_directory.c_str(),beam_origin.first,beam_origin.second));
1098     c1 -> Clear();
1099
1100     z_pos_inner_outer -> Draw("colz0");
1101     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1102     c1 -> Print(Form("%s/z_pos_inner_outer.pdf",output_directory.c_str()));
1103     c1 -> Clear();
1104
1105     z_pos_inner -> Draw("hist");
1106     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1107     c1 -> Print(Form("%s/z_pos_inner.pdf",output_directory.c_str()));
1108     c1 -> Clear();
1109
1110     z_pos_outer -> Draw("hist");
1111     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1112     c1 -> Print(Form("%s/z_pos_outer.pdf",output_directory.c_str()));
1113     c1 -> Clear();
1114
1115     // Nhits_good -> Draw("hist");
1116     // ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1117     // c1 -> Print(Form("%s/Nhits_good.pdf",output_directory.c_str()));
1118     // c1 -> Clear();
1119
1120     angle_correlation -> Draw("colz0");
1121     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1122     c1 -> Print(Form("%s/angle_correlation.pdf",output_directory.c_str()));
1123     c1 -> Clear();
1124
1125     z_pos_diff -> Draw("hist");
1126     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1127     c1 -> Print(Form("%s/z_pos_diff.pdf",output_directory.c_str()));
1128     c1 -> Clear();
1129
1130     z_pos_diff_angle_diff -> Draw("colz0");
1131     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1132     c1 -> Print(Form("%s/z_pos_diff_angle_diff.pdf",output_directory.c_str()));
1133     c1 -> Clear();
1134
1135     angle_diff -> Draw("hist");
1136     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1137     c1 -> Print(Form("%s/angle_diff.pdf",output_directory.c_str()));
1138     c1 -> Clear();
1139
1140
1141     angle_diff_DCA_dist -> Draw("colz0");
1142     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1143     c1 -> Print(Form("%s/angle_diff_DCA_dist.pdf",output_directory.c_str()));
1144     c1 -> Clear();
1145
1146     angle_diff_inner_phi -> Draw("colz0");
1147     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1148     c1 -> Print(Form("%s/angle_diff_inner_phi.pdf",output_directory.c_str()));
1149     c1 -> Clear();
1150
1151     angle_diff_outer_phi -> Draw("colz0");
1152     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1153     c1 -> Print(Form("%s/angle_diff_outer_phi.pdf",output_directory.c_str()));
1154     c1 -> Clear();
1155
1156     inner_pos_xy -> Draw("colz0");
1157     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1158     c1 -> Print(Form("%s/inner_pos_xy.pdf",output_directory.c_str()));
1159     c1 -> Clear();
1160
1161     outer_pos_xy -> Draw("colz0");
1162     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1163     c1 -> Print(Form("%s/outer_pos_xy.pdf",output_directory.c_str()));
1164     c1 -> Clear();
1165
1166     inner_outer_pos_xy -> Draw("colz0");
1167     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Simulation");
1168     c1 -> Print(Form("%s/inner_outer_pos_xy.pdf",output_directory.c_str()));
1169     c1 -> Clear();
1170
1171     TFile * INTT_MC_plot_out = new TFile(Form("%s/MC_angle_diff.root",output_directory.c_str()), "recreate");
1172     angle_diff -> Write();
1173     tree_out -> Write("", TObject::kOverwrite);
1174     INTT_MC_plot_out -> Close();
1175
1176     cout<<"the geatest N clu event under the fNhits cut : "<<greatest_N_clu<<endl;
1177 }