Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:11:56

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