Back to home page

sPhenix code displayed by LXR

 
 

    


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

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