analysis_INTT_CW/old_code/dNdeta_Z.C

0001 // #include "DAC_Scan_ladder.h"
0002 //#include "InttConversion.h"
0003 #include "/sphenix/user/ChengWei/INTT/INTT_commissioning/INTT_CW/INTT_commissioning/DAC_Scan/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 struct z_str {
0058     int nclu_inner;
0059     int nclu_outer;
0060     double zvtx;
0061     double zvtxE;
0062     double rangeL;
0063     double rangeR;
0064     int N_good;
0065     double width_density;
0066 };
0067
0068 double get_radius(double x, double y)
0069 {
0070     return sqrt(pow(x,2)+pow(y,2));
0071 }
0072
0073 double get_radius_sign(double x, double y)
0074 {
0075     double phi = ((y) < 0) ? atan2((y),(x)) * (180./TMath::Pi()) + 360 : atan2((y),(x)) * (180./TMath::Pi());
0076
0077     return (phi > 180) ? sqrt(pow(x,2)+pow(y,2)) * -1 : sqrt(pow(x,2)+pow(y,2));
0078 }
0079
0080 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)
0081 {
0082     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
0083     pad -> SetLeftMargin   (left);
0084     pad -> SetRightMargin  (right);
0085     pad -> SetTopMargin    (top);
0086     pad -> SetBottomMargin (bottom);
0087     pad -> SetTicks(1,1);
0088     if (set_logY == true)
0089     {
0090         pad -> SetLogy (1);
0091     }
0092
0093 }
0094
0095 std::vector<double> calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y) {
0096
0097     if (x1 != x2)
0098     {
0099         // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
0100         double a = (y2 - y1) / (x2 - x1);
0101         double b = y1 - a * x1;
0102
0103         // cout<<"slope : y="<<a<<"x+"<<b<<endl;
0104
0105         // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
0106         double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
0107
0108         // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
0109         double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
0110         double Yc = a * Xc + b;
0111
0112         return { closest_distance, Xc, Yc };
0113     }
0114     else
0115     {
0116         double closest_distance = std::abs(x1 - target_x);
0117         double Xc = x1;
0118         double Yc = target_y;
0119
0120         return { closest_distance, Xc, Yc };
0121     }
0122
0123
0124 }
0125
0126 double get_z_vertex(clu_info inner_clu, clu_info outer_clu, double target_x, double target_y)
0127 {
0128     // note : x = z,
0129     // note : y = radius
0130     double inner_clu_r = sqrt(pow(inner_clu.x,2)+ pow(inner_clu.y,2));
0131     double outer_clu_r = sqrt(pow(outer_clu.x,2)+ pow(outer_clu.y,2));
0132     double target_r    = sqrt(pow(target_x,2)   + pow(target_y,2));
0133
0134     // Use the slope equation (y = ax + b) to calculate the x-coordinate for the target y
0135     if ( fabs(outer_clu.z - inner_clu.z) < 0.00001 ){
0136         return outer_clu.z;
0137     }
0138     else {
0139         double slope = (outer_clu_r - inner_clu_r) / (outer_clu.z - inner_clu.z);
0140         double yIntercept = inner_clu_r - slope * inner_clu.z;
0141         double xCoordinate = (target_r - yIntercept) / slope;
0142         return xCoordinate;
0143     }
0144
0145 }
0146
0147 // Function to calculate the angle between two vectors in degrees using the cross product
0148 double calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
0149     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
0150     double vector1X = x2 - x1;
0151     double vector1Y = y2 - y1;
0152
0153     double vector2X = targetX - x1;
0154     double vector2Y = targetY - y1;
0155
0156     // Calculate the cross product of vector_1 and vector_2 (z-component)
0157     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
0158
0159     // cout<<" crossProduct : "<<crossProduct<<endl;
0160
0161     // Calculate the magnitudes of vector_1 and vector_2
0162     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
0163     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
0164
0165     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
0166     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
0167
0168     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
0169     // Convert the angle from radians to degrees and return it
0170     double angleInDegrees = angleInRadians * 180.0 / M_PI;
0171
0172     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
0173     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
0174
0175     // cout<<"angle : "<<angleInDegrees_new<<endl;
0176
0177     double DCA_distance = sin(angleInRadians_new) * magnitude2;
0178
0179     return DCA_distance;
0180 }
0181
0182 double grEY_stddev (TGraphErrors * input_grr)
0183 {
0184     vector<double> input_vector; input_vector.clear();
0185     for (int i = 0; i < input_grr -> GetN(); i++)
0186     {
0187         input_vector.push_back( input_grr -> GetPointY(i) );
0188     }
0189
0190     double sum=0;
0191     double average;
0192     double sum_subt = 0;
0193     for (int i=0; i<input_vector.size(); i++)
0194         {
0195             sum+=input_vector[i];
0196
0197         }
0198     average=sum/input_vector.size();
0199     //cout<<"average is : "<<average<<endl;
0200
0201     for (int i=0; i<input_vector.size(); i++)
0202         {
0203             //cout<<input_vector[i]-average<<endl;
0204             sum_subt+=pow((input_vector[i]-average),2);
0205
0206         }
0207     //cout<<"sum_subt : "<<sum_subt<<endl;
0208     double stddev;
0209     stddev=sqrt(sum_subt/(input_vector.size()-1));
0210     return stddev;
0211 }
0212
0213 pair<double, double> mirrorPolynomial(double a, double b) {
0214     // Interchange 'x' and 'y'
0215     double mirroredA = 1.0 / a;
0216     double mirroredB = -b / a;
0217
0218     return {mirroredA, mirroredB};
0219 }
0220
0221 // note : pair<reduced chi2, eta of the track>
0222 // note : vector : {r,z}
0223 // note : p0 : vertex point {r,z,zerror}
0224 // note : p1 : inner layer
0225 // note : p2 : outer layer
0226 pair<double,double> Get_eta(vector<double>p0, vector<double>p1, vector<double>p2)
0227 {
0228     vector<double> r_vec  = {p0[0], p1[0], p2[0]};
0229     vector<double> z_vec  = {p0[1], p1[1], p2[1]};
0230     vector<double> re_vec = {0, 0, 0};
0231     vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10., ( fabs( p2[1] ) < 130 ) ? 8. : 10.};
0232
0233     // note : swap the r and z, to have easier fitting
0234     // note : in principle, Z should be in the x axis, R should be in the Y axis
0235     TGraphErrors * track_gr = new TGraphErrors(3,&r_vec[0],&z_vec[0],&re_vec[0],&ze_vec[0]);
0236
0237     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
0238
0239     if (vertical_line) {return {0,0};}
0240     else
0241     {
0242         TF1 * fit_rz = new TF1("fit_rz","pol1",-500,500);
0243         fit_rz -> SetParameters(0,0);
0244
0245         track_gr -> Fit(fit_rz,"NQ");
0246
0247         pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
0248
0249         track_gr -> Delete();
0250         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
0251
0252     }
0253
0254 }
0255
0256 // note : vector : {r,z}
0257 // note : p0 : vertex point {r,z,zerror}
0258 // note : p1 : another point from detector
0259 // note : since only two points -> no strip width considered
0260 double Get_eta_2P(vector<double>p0, vector<double>p1)
0261 {
0262     return  -1 * TMath::Log( fabs( tan( atan2(p1[0] - p0[0], p1[1] - p0[1]) / 2 ) ) );
0263 }
0264
0265 double Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
0266 {
0267     if ( fabs(p0x - p1x) < 0.00001 ){ // note : the line is vertical (if z is along the x axis)
0268         return p0x;
0269     }
0270     else {
0271         double slope = (p1y - p0y) / (p1x - p0x);
0272         double yIntercept = p0y - slope * p0x;
0273         double xCoordinate = (given_y - yIntercept) / slope;
0274         return xCoordinate;
0275     }
0276 }
0277
0278 pair<double,double> Get_possible_zvtx(double rvtx, vector<double> p0, vector<double> p1) // note : inner p0, outer p1, vector {r,z}, -> {y,x}
0279 {
0280     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}
0281     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}
0282
0283     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
0284     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
0285
0286     double mid_point = (edge_first + edge_second) / 2.;
0287     double possible_width = fabs(edge_first - edge_second) / 2.;
0288
0289     return {mid_point, possible_width}; // note : first : mid point, second : width
0290
0291 }
0292
0293 vector<z_str> Get_INTT_Z (string full_directory, int z_capacity)
0294 {
0295     cout<<"----------------- loading the Z information ----------------- "<<endl;
0296     int z_eID, z_nclu_inner, z_nclu_outer, z_N_good;
0297     double z_zvtx, z_zvtxE, z_rangeL, z_rangeR, z_width_density;
0298
0299     TFile * z_file_in = new TFile( Form("%s",full_directory.c_str()) ,"read");
0300     TTree * z_tree =  (TTree *)z_file_in->Get("tree_Z");
0301     z_tree -> SetBranchAddress("eID",&z_eID);
0302     z_tree -> SetBranchAddress("nclu_inner",&z_nclu_inner);
0303     z_tree -> SetBranchAddress("nclu_outer",&z_nclu_outer);
0304     z_tree -> SetBranchAddress("zvtx",&z_zvtx);
0305     z_tree -> SetBranchAddress("zvtxE",&z_zvtxE);
0306     z_tree -> SetBranchAddress("rangeL",&z_rangeL);
0307     z_tree -> SetBranchAddress("rangeR",&z_rangeR);
0308     z_tree -> SetBranchAddress("N_good",&z_N_good);
0309     z_tree -> SetBranchAddress("Width_density",&z_width_density);
0310
0311     vector<z_str> out_vec(z_capacity,{-1,-1,-1,-1,-1,-1,-1,-1});
0312     // cout<<"test1"<<endl;
0313
0314     for (int i = 0; i < z_tree -> GetEntries(); i++){
0315         z_tree -> GetEntry(i);
0316         // cout<<"test 0 : "<<i<<"  "<<z_eID<<endl;
0317         out_vec[z_eID] = {z_nclu_inner,z_nclu_outer,z_zvtx,z_zvtxE,z_rangeL,z_rangeR,z_N_good,z_width_density};
0318     }
0319
0320
0321     cout<<"-----------------  INTT Z info. loader done -----------------"<<endl;
0322
0323     return out_vec;
0324 }
0325
0326 // note : use "ls *.root > file_list.txt" to create the list of the file in the folder, full directory in the file_list.txt
0327 // note : set_folder_name = "folder_xxxx"
0328 // note : server_name = "inttx"
0329 void dNdeta_Z(/*pair<double,double>beam_origin*/)
0330 {
0331
0332     TCanvas * c4 = new TCanvas("","",850,800);
0333     c4 -> cd();
0334
0335     TCanvas * c2 = new TCanvas("","",2500,800);
0336     c2 -> cd();
0337     TPad *pad_xy = new TPad(Form("pad_xy"), "", 0.0, 0.0, 0.33, 1.0);
0338     Characterize_Pad(pad_xy, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0339     pad_xy -> Draw();
0340
0341     TPad *pad_rz = new TPad(Form("pad_rz"), "", 0.33, 0.0, 0.66, 1.0);
0342     Characterize_Pad(pad_rz, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0343     pad_rz -> Draw();
0344
0345     TPad *pad_z = new TPad(Form("pad_z"), "", 0.66, 0.0, 1.0, 1.0);
0346     Characterize_Pad(pad_z, 0.15, 0.1, 0.1, 0.1 , 0, 0);
0347     pad_z -> Draw();
0348
0349     TCanvas * c1 = new TCanvas("","",1000,800);
0350     c1 -> cd();
0351
0352
0353
0354     // string mother_folder_directory = "/home/phnxrc/INTT/cwshih/DACscan_data/zero_magnet_Takashi_used";
0355     // string file_name = "beam_inttall-00020869-0000_event_base_ana_cluster_ideal_excludeR1500_100kEvent";
0356     // string file_name = "beam_inttall-00020869-0000_event_base_ana_cluster_survey_1_XYAlpha_Peek_3.32mm_excludeR1500_100kEvent";
0357
0358     string mother_folder_directory = "/sphenix/user/ChengWei/INTT/INTT_commissioning/ZeroField/20869";
0359     string file_name = "beam_inttall-00020869-0000_event_base_ana_cluster_survey_1_XYAlpha_Peek_3.32mm_excludeR20000_150kEvent_3HotCut";
0360
0361     // string mother_folder_directory = "/home/phnxrc/INTT/cwshih/DACscan_data/new_DAC_Scan_0722/AllServer/DAC2";
0362     // string file_name = "beam_inttall-00023058-0000_event_base_ana_cluster_ideal_excludeR2000_100kEvent";
0363
0364     system(Form("mkdir %s/folder_%s_advanced",mother_folder_directory.c_str(),file_name.c_str()));
0365     system(Form("mkdir %s/folder_%s_advanced/good_track",mother_folder_directory.c_str(),file_name.c_str()));
0366
0367     // note : new updated centrality with the fNhits cut 40k.
0368     // note :                   0      5   10     15    20    25    30    35    40   45   50
0369     int centrality_line[20] = {8000, 4822, 3897, 3117, 2476, 1964, 1564, 1242, 985, 772, 598, 456, 347, 263, 195, 139, 91, 47, 13, 3};
0370     pair<int,int>centrality_interval = {10,15};
0371     pair<double,double> beam_origin = {-0,2};
0372     double temp_Y_align = 0.;
0373     double temp_X_align = -0.;
0374
0375     double zvtx_hist_l = -500;
0376     double zvtx_hist_r = 500;
0377
0378     int Nhit_cut = 15000;           // note : if (> Nhit_cut)          -> continue
0379     int Nhit_cutl = 500;          // note : if (< Nhit_cutl)         -> continue
0380     int clu_size_cut = 4;           // note : if (> clu_size_cut)      -> continue
0381     double clu_sum_adc_cut = 31;    // note : if (< clu_sum_adc_cut)   -> continue
0382     int N_clu_cut = centrality_line[ centrality_interval.first/5 ];          // note : if (> N_clu_cut)         -> continue  unit number
0383     int N_clu_cutl = centrality_line[ centrality_interval.second/5 ];;           // note : if (< N_clu_cutl)        -> continue  unit number
0384     double phi_diff_cut = 3.5;      // note : if (< phi_diff_cut)      -> pass      unit degree
0385     double DCA_cut = 4;             // note : if (< DCA_cut)           -> pass      unit mm
0386     int zvtx_cal_require = 15;      // note : if (> zvtx_cal_require)  -> pass
0387     int zvtx_draw_requireL = 15;
0388     int zvtx_draw_requireR = 40000;   // note : if ( zvtx_draw_requireL < event && event < zvtx_draw_requireR) -> pass
0389     double Integrate_portion = 0.4; // todo : it was 0.6826, try to require less event, as most of the tracklets are not that useful
0390     double Integrate_portion_final = 0.68;
0391     bool draw_event_display = false;
0392     int print_rate = 5;
0393     int z_capacity = 20000; // todo : filled the Z capacity
0394
0395     double mean_zvtx = -198.38; // note : unit : mm
0396
0397     cout<<" Centality line : "<<N_clu_cutl<<" "<<N_clu_cut<<endl;
0398
0399     bool full_event = false;
0400     long long used_event = 10000;
0401
0402     double dNdeta_upper_range = 650;
0403
0404     int geo_mode_id = 1;
0405     string conversion_mode = (geo_mode_id == 0) ? "ideal" : "survey_1_XYAlpha_Peek";
0406     double peek = 3.32405;
0407
0408
0409     vector<z_str> INTT_Z = Get_INTT_Z (Form("%s/folder_%s_advanced/INTT_zvtx.root",mother_folder_directory.c_str(),file_name.c_str()), z_capacity);
0410
0411     TFile * file_in = new TFile(Form("%s/%s.root",mother_folder_directory.c_str(),file_name.c_str()),"read");
0412     TTree * tree = (TTree *)file_in->Get("tree_clu");
0413
0414     long long N_event = (full_event == true) ? tree -> GetEntries() : used_event;
0415     cout<<Form("N_event in file %s : %lli",file_name.c_str(), N_event)<<endl;
0416
0417     int N_hits;
0418     int N_cluster_inner;
0419     int N_cluster_outer;
0420     vector<int>* column_vec = new vector<int>();
0421     vector<double>* avg_chan_vec = new vector<double>();
0422     vector<int>* sum_adc_vec = new vector<int>();
0423     vector<int>* sum_adc_conv_vec = new vector<int>();
0424     vector<int>* size_vec = new vector<int>();
0425     vector<double>* x_vec = new vector<double>();
0426     vector<double>* y_vec = new vector<double>();
0427     vector<double>* z_vec = new vector<double>();
0428     vector<int>* layer_vec = new vector<int>();
0429     vector<double>* phi_vec = new vector<double>();
0430
0431     tree -> SetBranchStatus("*",0);
0432     tree -> SetBranchStatus("nhits",1);
0433     tree -> SetBranchStatus("nclu_inner",1);
0434     tree -> SetBranchStatus("nclu_outer",1);
0435     tree -> SetBranchStatus("column",1);
0436     tree -> SetBranchStatus("avg_chan",1);
0437     tree -> SetBranchStatus("sum_adc",1);
0438     tree -> SetBranchStatus("sum_adc_conv",1);
0439     tree -> SetBranchStatus("size",1);
0440     tree -> SetBranchStatus("x",1);
0441     tree -> SetBranchStatus("y",1);
0442     tree -> SetBranchStatus("z",1);
0443     tree -> SetBranchStatus("layer",1);
0444     tree -> SetBranchStatus("phi",1);
0445
0446     tree -> SetBranchAddress("nhits",&N_hits);
0447     tree -> SetBranchAddress("nclu_inner",&N_cluster_inner);
0448     tree -> SetBranchAddress("nclu_outer",&N_cluster_outer);
0449
0450     tree -> SetBranchAddress("column", &column_vec);
0451     tree -> SetBranchAddress("avg_chan", &avg_chan_vec);
0452     tree -> SetBranchAddress("sum_adc", &sum_adc_vec);
0453     tree -> SetBranchAddress("sum_adc_conv", &sum_adc_conv_vec);
0454     tree -> SetBranchAddress("size", &size_vec);
0455     tree -> SetBranchAddress("x", &x_vec);
0456     tree -> SetBranchAddress("y", &y_vec);
0457     tree -> SetBranchAddress("z", &z_vec);
0458     tree -> SetBranchAddress("layer", &layer_vec);
0459     tree -> SetBranchAddress("phi", &phi_vec);
0460
0461     TFile * out_file = new TFile(Form("%s/folder_%s_advanced/INTT_dNdeta.root",mother_folder_directory.c_str(),file_name.c_str()),"RECREATE");
0462
0463     int N_hits_out, N_cluster_inner_out, N_cluster_outer_out, ntrack_out;
0464     int eID_out;
0465     double out_xvtx, out_yvtx, out_zvtx, out_zvtxE;
0466     vector<double> out_inner_x;
0467     vector<double> out_inner_y;
0468     vector<double> out_inner_z;
0469     vector<double> out_inner_r;
0470     vector<double> out_inner_phi;
0471     vector<double> out_inner_zE;
0472     vector<double> out_outer_x;
0473     vector<double> out_outer_y;
0474     vector<double> out_outer_z;
0475     vector<double> out_outer_r;
0476     vector<double> out_outer_phi;
0477     vector<double> out_outer_zE;
0478     vector<double> out_eta;
0479     vector<double> out_eta_rchi2;
0480
0481     TTree * tree_out =  new TTree ("tree_eta", "eta info.");
0482     tree_out -> Branch("eID",&eID_out);
0483     tree_out -> Branch("nhits",&N_hits_out);
0484     tree_out -> Branch("nclu_inner",&N_cluster_inner_out);
0485     tree_out -> Branch("nclu_outer",&N_cluster_outer_out);
0486     tree_out -> Branch("ntrack",&ntrack_out);
0487     tree_out -> Branch("xvtx",&out_xvtx);
0488     tree_out -> Branch("yvtx",&out_yvtx);
0489     tree_out -> Branch("zvtx",&out_zvtx);
0490     tree_out -> Branch("zvtxE",&out_zvtxE);
0491
0492     tree_out -> Branch("inner_x",&out_inner_x); // note : inner
0493     tree_out -> Branch("inner_y",&out_inner_y);
0494     tree_out -> Branch("inner_z",&out_inner_z);
0495     tree_out -> Branch("inner_r",&out_inner_r);
0496     tree_out -> Branch("inner_phi",&out_inner_phi);
0497     tree_out -> Branch("inner_zE",&out_inner_zE);
0498     tree_out -> Branch("outer_x",&out_outer_x); // note : outer
0499     tree_out -> Branch("outer_y",&out_outer_y);
0500     tree_out -> Branch("outer_z",&out_outer_z);
0501     tree_out -> Branch("outer_r",&out_outer_r);
0502     tree_out -> Branch("outer_phi",&out_outer_phi);
0503     tree_out -> Branch("outer_zE",&out_outer_zE);
0504     tree_out -> Branch("eta",&out_eta); // note : track
0505     tree_out -> Branch("eta_rchi",&out_eta_rchi2);
0506
0507     TLatex *draw_text = new TLatex();
0508     draw_text -> SetNDC();
0509     draw_text -> SetTextSize(0.02);
0510
0511     vector<clu_info> temp_sPH_inner_nocolumn_vec; temp_sPH_inner_nocolumn_vec.clear();
0512     vector<clu_info> temp_sPH_outer_nocolumn_vec; temp_sPH_outer_nocolumn_vec.clear();
0513     vector<vector<double>> temp_sPH_nocolumn_vec(2);
0514     vector<vector<double>> temp_sPH_nocolumn_rz_vec(2);
0515
0516     TH1F * dNdeta_hist = new TH1F("","",29,-2.9,2.9);
0517     // TH1F * dNdeta_hist = new TH1F("","",145,-2.9,2.9);
0518     dNdeta_hist -> SetMarkerStyle(20);
0519     dNdeta_hist -> SetMarkerSize(0.8);
0520     dNdeta_hist -> SetMarkerColor(TColor::GetColor("#1A3947"));
0521     dNdeta_hist -> SetLineColor(TColor::GetColor("#1A3947"));
0522     dNdeta_hist -> SetLineWidth(2);
0523     dNdeta_hist -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0524     dNdeta_hist -> GetXaxis() -> SetTitle("#eta");
0525     dNdeta_hist -> GetYaxis() -> SetRangeUser(0,50);
0526     dNdeta_hist -> SetTitleSize(0.06, "X");
0527     dNdeta_hist -> SetTitleSize(0.06, "Y");
0528     dNdeta_hist -> GetXaxis() -> SetTitleOffset (0.71);
0529     dNdeta_hist -> GetYaxis() -> SetTitleOffset (1.1);
0530     dNdeta_hist -> GetXaxis() -> CenterTitle(true);
0531     dNdeta_hist -> GetYaxis() -> CenterTitle(true);
0532
0533     // TH1F * dNdeta_2P_inner_hist = new TH1F("","",29,-2.9,2.9);
0534     TH1F * dNdeta_2P_inner_hist = new TH1F("","",145,-2.9,2.9);
0535     dNdeta_2P_inner_hist -> SetMarkerStyle(20);
0536     dNdeta_2P_inner_hist -> SetMarkerSize(0.8);
0537     dNdeta_2P_inner_hist -> SetMarkerColor(TColor::GetColor("#1A3947"));
0538     dNdeta_2P_inner_hist -> SetLineColor(TColor::GetColor("#1A3947"));
0539     dNdeta_2P_inner_hist -> SetLineWidth(2);
0540     dNdeta_2P_inner_hist -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0541     dNdeta_2P_inner_hist -> GetXaxis() -> SetTitle("#eta");
0542     dNdeta_2P_inner_hist -> GetYaxis() -> SetRangeUser(0,50);
0543     dNdeta_2P_inner_hist -> SetTitleSize(0.06, "X");
0544     dNdeta_2P_inner_hist -> SetTitleSize(0.06, "Y");
0545     dNdeta_2P_inner_hist -> GetXaxis() -> SetTitleOffset (0.71);
0546     dNdeta_2P_inner_hist -> GetYaxis() -> SetTitleOffset (1.1);
0547     dNdeta_2P_inner_hist -> GetXaxis() -> CenterTitle(true);
0548     dNdeta_2P_inner_hist -> GetYaxis() -> CenterTitle(true);
0549
0550     TH1F * dNdeta_2P_outer_hist = new TH1F("","",29,-2.9,2.9);
0551     dNdeta_2P_outer_hist -> SetMarkerStyle(20);
0552     dNdeta_2P_outer_hist -> SetMarkerSize(0.8);
0553     dNdeta_2P_outer_hist -> SetMarkerColor(TColor::GetColor("#1A3947"));
0554     dNdeta_2P_outer_hist -> SetLineColor(TColor::GetColor("#1A3947"));
0555     dNdeta_2P_outer_hist -> SetLineWidth(2);
0556     dNdeta_2P_outer_hist -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0557     dNdeta_2P_outer_hist -> GetXaxis() -> SetTitle("#eta");
0558     dNdeta_2P_outer_hist -> GetYaxis() -> SetRangeUser(0,50);
0559     dNdeta_2P_outer_hist -> SetTitleSize(0.06, "X");
0560     dNdeta_2P_outer_hist -> SetTitleSize(0.06, "Y");
0561     dNdeta_2P_outer_hist -> GetXaxis() -> SetTitleOffset (0.71);
0562     dNdeta_2P_outer_hist -> GetYaxis() -> SetTitleOffset (1.1);
0563     dNdeta_2P_outer_hist -> GetXaxis() -> CenterTitle(true);
0564     dNdeta_2P_outer_hist -> GetYaxis() -> CenterTitle(true);
0565
0566     double N_good_event = 0;
0567
0568     TH2F * Good_inner_outer_pos_xy_nzB = new TH2F("","Good_inner_outer_pos_xy_nzB",360,-150,150,360,-150,150);
0569     Good_inner_outer_pos_xy_nzB -> SetStats(0);
0570     Good_inner_outer_pos_xy_nzB -> GetXaxis() -> SetTitle("X axis");
0571     Good_inner_outer_pos_xy_nzB -> GetYaxis() -> SetTitle("Y axis");
0572
0573     TH2F * Good_inner_outer_pos_xy_nzA = new TH2F("","Good_inner_outer_pos_xy_nzA",360,-150,150,360,-150,150);
0574     Good_inner_outer_pos_xy_nzA -> SetStats(0);
0575     Good_inner_outer_pos_xy_nzA -> GetXaxis() -> SetTitle("X axis");
0576     Good_inner_outer_pos_xy_nzA -> GetYaxis() -> SetTitle("Y axis");
0577
0578     TH2F * Good_inner_outer_pos_xy_pzA = new TH2F("","Good_inner_outer_pos_xy_pzA",360,-150,150,360,-150,150);
0579     Good_inner_outer_pos_xy_pzA -> SetStats(0);
0580     Good_inner_outer_pos_xy_pzA -> GetXaxis() -> SetTitle("X axis");
0581     Good_inner_outer_pos_xy_pzA -> GetYaxis() -> SetTitle("Y axis");
0582
0583     TH2F * Good_inner_outer_pos_xy_pzB = new TH2F("","Good_inner_outer_pos_xy_pzB",360,-150,150,360,-150,150);
0584     Good_inner_outer_pos_xy_pzB -> SetStats(0);
0585     Good_inner_outer_pos_xy_pzB -> GetXaxis() -> SetTitle("X axis");
0586     Good_inner_outer_pos_xy_pzB -> GetYaxis() -> SetTitle("Y axis");
0587
0588     TH2F * Good_track_space = new TH2F("","Good_track_space",200,-300,300,200,-300,300);
0589     Good_track_space -> SetStats(0);
0590     Good_track_space -> GetXaxis() -> SetTitle("X axis");
0591     Good_track_space -> GetYaxis() -> SetTitle("Y axis");
0592
0593     // dNdeta_hist -> GetXaxis() -> SetLimits(-10,10);
0594     // dNdeta_hist -> GetXaxis() -> SetNdivisions  (505);
0595
0596     TF1 * zvtx_finder = new TF1("zvtx_finder","pol0",-1,3000);
0597
0598
0599     vector<vector<double>> good_track_xy_vec; good_track_xy_vec.clear();
0600     vector<vector<double>> good_track_rz_vec; good_track_rz_vec.clear();
0601     vector<vector<double>> good_comb_rz_vec; good_comb_rz_vec.clear();
0602     vector<vector<double>> good_comb_xy_vec; good_comb_xy_vec.clear();
0603     vector<vector<double>> good_comb_phi_vec; good_comb_phi_vec.clear();
0604     vector<vector<double>> good_tracklet_rz; good_tracklet_rz.clear();
0605     TLine * track_line = new TLine();
0606     track_line -> SetLineWidth(1);
0607     track_line -> SetLineColorAlpha(38,0.5);
0608
0609     TLine * coord_line = new TLine();
0610     coord_line -> SetLineWidth(1);
0611     coord_line -> SetLineColor(16);
0612     coord_line -> SetLineStyle(2);
0613
0614
0615     int N_event_pass_number = 0;
0616
0617     vector<double> DCA_info_vec; DCA_info_vec.clear();
0618     pair<double,double> Get_eta_pair;
0619
0620     int inner_1_check = 0;  int inner_2_check = 0;  int inner_3_check = 0;  int inner_4_check = 0;
0621     int outer_1_check = 0;  int outer_2_check = 0;  int outer_3_check = 0;  int outer_4_check = 0;
0622     unsigned int length;
0623     vector<int> used_outer_check_eta; used_outer_check_eta.clear();
0624
0625     for (int event_i = 0; event_i < N_event; event_i++)
0626     {
0627         if (event_i % 100 == 0) cout<<"code running... "<<event_i<<endl;
0628
0629         tree -> GetEntry(event_i);
0630         length = column_vec -> size();
0631
0632         if (event_i == 13) cout<<"test, eID : "<<event_i<<" Nhits "<<N_hits<<endl;
0633
0634         if (N_hits > Nhit_cut) continue;
0635         if (N_hits < Nhit_cutl) continue;
0636
0637         N_event_pass_number += 1;
0638
0639         if (N_cluster_inner == 0 || N_cluster_outer == 0) continue;
0640         if (N_cluster_inner == -1 || N_cluster_outer == -1) continue;
0641         if ((N_cluster_inner + N_cluster_outer) < zvtx_cal_require) continue;
0642         if (N_cluster_inner < 30) continue;
0643         if (N_cluster_outer < 30) continue;
0644
0645         // cout<<"pass : "<<INTT_Z[event_i].zvtx<<" "<<INTT_Z[event_i].width_density<<endl;
0646
0647         if ( INTT_Z[event_i].N_good == -1) continue;
0648         if ( fabs(INTT_Z[event_i].zvtx - mean_zvtx) > 20 ) continue;
0649         if ( INTT_Z[event_i].width_density <= 0.3 ) continue;
0650
0651
0652
0653
0654         // note : apply some selection to remove the hot channels
0655         // note : and make the inner_clu_vec and outer_clu_vec
0656         for (int clu_i = 0; clu_i < length; clu_i++)
0657         {
0658             if (size_vec -> at(clu_i) > clu_size_cut) continue;
0659             // if (size_vec -> at(clu_i) < 2) continue;
0660             if (sum_adc_conv_vec -> at(clu_i) < clu_sum_adc_cut) continue;
0661             // if (z_vec -> at(clu_i) < 0) continue;
0662
0663             // note : inner
0664             // 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;
0665             // // 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;
0666             // if (layer_vec -> at(clu_i) == 0 && phi_vec -> at(clu_i) > 295 && phi_vec -> at(clu_i) < 302) continue;
0667             // if (layer_vec -> at(clu_i) == 0 && phi_vec -> at(clu_i) > 210 && phi_vec -> at(clu_i) < 213) continue;
0668             // if (layer_vec -> at(clu_i) == 0 && phi_vec -> at(clu_i) > 55 && phi_vec -> at(clu_i) < 65) continue;
0669             // if (layer_vec -> at(clu_i) == 0 && phi_vec -> at(clu_i) > 348 && phi_vec -> at(clu_i) < 353) continue;
0670             // if (layer_vec -> at(clu_i) == 0 && phi_vec -> at(clu_i) > 265 && phi_vec -> at(clu_i) < 270) continue;
0671
0672             // note : outer
0673             // 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;
0674             // // 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;
0675             // // 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;
0676             // 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;
0677             // if (layer_vec -> at(clu_i) == 1 && phi_vec -> at(clu_i) > 335 && phi_vec -> at(clu_i) < 340) continue;
0678             // if (layer_vec -> at(clu_i) == 1 && phi_vec -> at(clu_i) > 105 && phi_vec -> at(clu_i) < 115) continue;
0679             // 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"
0680
0681
0682             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) );
0683             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) );
0684
0685             double clu_radius = get_radius(
0686                 (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),
0687                 (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)
0688             );
0689             temp_sPH_nocolumn_rz_vec[0].push_back(z_vec -> at(clu_i));
0690             temp_sPH_nocolumn_rz_vec[1].push_back( ( phi_vec -> at(clu_i) > 180 ) ? clu_radius * -1 : clu_radius );
0691
0692
0693             if (layer_vec -> at(clu_i) == 0) // note : inner
0694                 temp_sPH_inner_nocolumn_vec.push_back({
0695                     column_vec -> at(clu_i),
0696                     avg_chan_vec -> at(clu_i),
0697                     sum_adc_vec -> at(clu_i),
0698                     sum_adc_conv_vec -> at(clu_i),
0699                     size_vec -> at(clu_i),
0700                     (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),
0701                     (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),
0702                     z_vec -> at(clu_i),
0703                     layer_vec -> at(clu_i),
0704                     phi_vec -> at(clu_i)
0705                 });
0706
0707             if (layer_vec -> at(clu_i) == 1) // note : outer
0708                 temp_sPH_outer_nocolumn_vec.push_back({
0709                     column_vec -> at(clu_i),
0710                     avg_chan_vec -> at(clu_i),
0711                     sum_adc_vec -> at(clu_i),
0712                     sum_adc_conv_vec -> at(clu_i),
0713                     size_vec -> at(clu_i),
0714                     (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),
0715                     (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),
0716                     z_vec -> at(clu_i),
0717                     layer_vec -> at(clu_i),
0718                     phi_vec -> at(clu_i)
0719                 });
0720         }
0721
0722         inner_1_check = 0;
0723         inner_2_check = 0;
0724         inner_3_check = 0;
0725         inner_4_check = 0;
0726         for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
0727             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 0 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 90)    inner_1_check = 1;
0728             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 90 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 180)  inner_2_check = 1;
0729             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 180 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 270) inner_3_check = 1;
0730             if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 270 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 360) inner_4_check = 1;
0731
0732             if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check) == 4 ) break;
0733         }
0734
0735         outer_1_check = 0;
0736         outer_2_check = 0;
0737         outer_3_check = 0;
0738         outer_4_check = 0;
0739         for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
0740             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 0 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 90)    outer_1_check = 1;
0741             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 90 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 180)  outer_2_check = 1;
0742             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 180 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 270) outer_3_check = 1;
0743             if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 270 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 360) outer_4_check = 1;
0744
0745             if ( (outer_1_check + outer_2_check + outer_3_check + outer_4_check) == 4 ) break;
0746         }
0747
0748         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)
0749         {
0750             good_track_xy_vec.clear();
0751             good_track_rz_vec.clear();
0752             good_comb_rz_vec.clear();
0753             good_comb_xy_vec.clear();
0754             good_comb_phi_vec.clear();
0755             temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
0756             temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
0757             temp_sPH_inner_nocolumn_vec.clear();
0758             temp_sPH_outer_nocolumn_vec.clear();
0759             continue;
0760         }
0761
0762         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 )
0763         {
0764             cout<<"some quater of INTT doens't work !! eID : "<<event_i<<endl;
0765
0766             good_track_xy_vec.clear();
0767             good_track_rz_vec.clear();
0768             good_comb_rz_vec.clear();
0769             good_comb_xy_vec.clear();
0770             good_comb_phi_vec.clear();
0771             temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
0772             temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
0773             temp_sPH_inner_nocolumn_vec.clear();
0774             temp_sPH_outer_nocolumn_vec.clear();
0775             continue;
0776         }
0777
0778         used_outer_check_eta.clear();
0779         used_outer_check_eta = vector<int>(temp_sPH_outer_nocolumn_vec.size(),0);
0780
0781         for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ ) {
0782             for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ ) {
0783
0784                 if (used_outer_check_eta[outer_i] == 1) continue;
0785                 if (fabs(temp_sPH_inner_nocolumn_vec[inner_i].phi - temp_sPH_outer_nocolumn_vec[outer_i].phi) > phi_diff_cut) continue;
0786
0787                 DCA_info_vec = calculateDistanceAndClosestPoint(
0788                     temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0789                     temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0790                     beam_origin.first, beam_origin.second
0791                 );
0792
0793                 if (DCA_info_vec[0] > DCA_cut) continue;
0794
0795                 Get_eta_pair = Get_eta(
0796                     {get_radius(beam_origin.first,beam_origin.second), INTT_Z[event_i].zvtx, INTT_Z[event_i].zvtxE},
0797                     {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},
0798                     {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}
0799                 );
0800
0801
0802                 // note : following cut means that the three points have good connection
0803                 if (fabs(Get_eta_pair.first) > 7) continue; // todo the chi2/NDF cut is here
0804
0805                 eID_out = event_i;
0806                 N_hits_out = N_hits;
0807                 N_cluster_inner_out = temp_sPH_inner_nocolumn_vec.size();
0808                 N_cluster_outer_out = temp_sPH_outer_nocolumn_vec.size();
0809                 out_xvtx = beam_origin.first;
0810                 out_yvtx = beam_origin.second;
0811                 out_zvtx = INTT_Z[event_i].zvtx;
0812                 out_zvtxE = INTT_Z[event_i].zvtxE;
0813
0814                 out_inner_x.push_back( temp_sPH_inner_nocolumn_vec[inner_i].x ); // note : inner
0815                 out_inner_y.push_back( temp_sPH_inner_nocolumn_vec[inner_i].y );
0816                 out_inner_z.push_back( temp_sPH_inner_nocolumn_vec[inner_i].z );
0817                 out_inner_r.push_back( get_radius(temp_sPH_inner_nocolumn_vec[inner_i].x,temp_sPH_inner_nocolumn_vec[inner_i].y) );
0818                 out_inner_phi.push_back( temp_sPH_inner_nocolumn_vec[inner_i].phi );
0819                 out_inner_zE.push_back( ( fabs( temp_sPH_inner_nocolumn_vec[inner_i].z ) < 130 ) ? 8. : 10. );
0820                 out_outer_x.push_back( temp_sPH_outer_nocolumn_vec[outer_i].x ); // note : outer
0821                 out_outer_y.push_back( temp_sPH_outer_nocolumn_vec[outer_i].y );
0822                 out_outer_z.push_back( temp_sPH_outer_nocolumn_vec[outer_i].z );
0823                 out_outer_r.push_back( get_radius(temp_sPH_outer_nocolumn_vec[outer_i].x,temp_sPH_outer_nocolumn_vec[outer_i].y) );
0824                 out_outer_phi.push_back( temp_sPH_outer_nocolumn_vec[outer_i].phi );
0825                 out_outer_zE.push_back( ( fabs( temp_sPH_outer_nocolumn_vec[outer_i].z ) < 130 ) ? 8. : 10. );
0826                 out_eta.push_back( Get_eta_pair.second );
0827                 out_eta_rchi2.push_back( Get_eta_pair.first );
0828
0829                 dNdeta_hist -> Fill(Get_eta_pair.second);
0830                 used_outer_check_eta[outer_i] = 1;
0831
0832                 DCA_info_vec.clear();
0833
0834                 break;
0835             } // note : for loop outer
0836         } // note : for loop inner
0837
0838         ntrack_out = out_eta.size();
0839         if (out_eta.size() > 0) {
0840             tree_out -> Fill();
0841             N_good_event += 1;
0842         }
0843
0844         out_inner_x.clear();
0845         out_inner_y.clear();
0846         out_inner_z.clear();
0847         out_inner_r.clear();
0848         out_inner_phi.clear();
0849         out_inner_zE.clear();
0850         out_outer_x.clear();
0851         out_outer_y.clear();
0852         out_outer_z.clear();
0853         out_outer_r.clear();
0854         out_outer_phi.clear();
0855         out_outer_zE.clear();
0856         out_eta.clear();
0857         out_eta_rchi2.clear();
0858
0859         good_track_xy_vec.clear();
0860         good_track_rz_vec.clear();
0861         good_comb_rz_vec.clear();
0862         good_comb_xy_vec.clear();
0863         good_comb_phi_vec.clear();
0864         temp_sPH_nocolumn_rz_vec.clear(); temp_sPH_nocolumn_rz_vec = vector<vector<double>>(2);
0865         temp_sPH_nocolumn_vec.clear(); temp_sPH_nocolumn_vec = vector<vector<double>>(2);
0866         temp_sPH_inner_nocolumn_vec.clear();
0867         temp_sPH_outer_nocolumn_vec.clear();
0868     } // note : end of event
0869
0870
0871     c1 -> Clear();
0872
0873     tree_out->SetDirectory(out_file);
0874     tree_out -> Write("", TObject::kOverwrite);
0875
0876
0877     c1 -> cd();
0878     // Good_inner_outer_pos_xy_nzB-> Draw("colz0");
0879     // c1 -> Print(Form("%s/folder_%s_advanced/Good_inner_outer_pos_xy_nzB.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0880     // c1 -> Clear();
0881     // Good_inner_outer_pos_xy_nzA-> Draw("colz0");
0882     // c1 -> Print(Form("%s/folder_%s_advanced/Good_inner_outer_pos_xy_nzA.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0883     // c1 -> Clear();
0884     // Good_inner_outer_pos_xy_pzA-> Draw("colz0");
0885     // c1 -> Print(Form("%s/folder_%s_advanced/Good_inner_outer_pos_xy_pzA.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0886     // c1 -> Clear();
0887     // Good_inner_outer_pos_xy_pzB-> Draw("colz0");
0888     // c1 -> Print(Form("%s/folder_%s_advanced/Good_inner_outer_pos_xy_pzB.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0889     // c1 -> Clear();
0890
0891     // Good_track_space -> Draw("");
0892     // for (int i = 0; i < good_tracklet_rz.size(); i++){
0893     //     track_line -> DrawLine(good_tracklet_rz[i][1],good_tracklet_rz[i][0],good_tracklet_rz[i][3],good_tracklet_rz[i][2]);
0894     // }
0895     // c1 -> Print(Form("%s/folder_%s_advanced/good_tracklet_rz.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0896     // c1 -> Clear();
0897
0898     //   note : normalized
0899     dNdeta_hist -> Scale(1./double(N_good_event * dNdeta_hist -> GetBinWidth(1) ));
0900     dNdeta_hist -> GetYaxis() -> SetRangeUser(0,dNdeta_upper_range);
0901     dNdeta_2P_inner_hist -> Scale(1./double(N_good_event));
0902     dNdeta_2P_inner_hist -> GetYaxis() -> SetRangeUser(0,dNdeta_upper_range);
0903     dNdeta_2P_outer_hist -> Scale(1./double(N_good_event));
0904     dNdeta_2P_outer_hist -> GetYaxis() -> SetRangeUser(0,dNdeta_upper_range);
0905
0906     cout<<" final N good event : "<<N_good_event<<" N event with correct hit : "<<N_event_pass_number<<endl;
0907
0908
0909     TCanvas * c3 = new TCanvas("c3","c3",900,800); c3 -> cd();
0910     TPad *pad_obj = new TPad(Form("pad_obj"), "", 0.0, 0.0, 1.0, 1.0);
0911     Characterize_Pad(pad_obj, 0.18,  0.1,  0.1,  0.12, 0, 0);
0912     pad_obj -> SetTicks(1,1);
0913     pad_obj -> Draw();
0914     pad_obj -> cd();
0915
0916     SetsPhenixStyle();
0917     dNdeta_hist -> Draw("ep");
0918
0919
0920     TLatex *ltx = new TLatex();
0921     ltx->SetNDC();
0922     ltx->SetTextSize(0.045);
0923     ltx->DrawLatex(gPad->GetLeftMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Work-in-progress");
0924     ltx->DrawLatex(0.48, 0.835, "Run 20869");
0925     ltx->DrawLatex(0.48, 0.785, "Au+Au #sqrt{s_{NN}} = 200 GeV");
0926     ltx->DrawLatex(0.48, 0.735, Form("%i%% - %i%%",centrality_interval.first,centrality_interval.second));
0927     c3 -> Print(Form("%s/folder_%s_advanced/dNdeta_%i_%i.pdf",mother_folder_directory.c_str(),file_name.c_str(),centrality_interval.first,centrality_interval.second));
0928     pad_obj -> Clear();
0929
0930     // pad_obj -> cd();
0931     // dNdeta_2P_inner_hist -> Draw("ep");
0932     // ltx->DrawLatex(gPad->GetLeftMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Work-in-progress");
0933     // ltx->DrawLatex(0.38, 0.835, "Run 20869");
0934     // ltx->DrawLatex(0.38, 0.785, "Au+Au #sqrt{s_{NN}} = 200 GeV");
0935     // ltx->DrawLatex(0.38, 0.735, "Tracklet : vertex + inner cluster");
0936     // c3 -> Print(Form("%s/folder_%s_advanced/dNdeta_2P_inner.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0937     // pad_obj -> Clear();
0938
0939     // pad_obj -> cd();
0940     // dNdeta_2P_outer_hist -> Draw("ep");
0941     // ltx->DrawLatex(gPad->GetLeftMargin(), 1 - gPad->GetTopMargin() + 0.01, "#it{#bf{sPHENIX INTT}} Work-in-progress");
0942     // ltx->DrawLatex(0.38, 0.835, "Run 20869");
0943     // ltx->DrawLatex(0.38, 0.785, "Au+Au #sqrt{s_{NN}} = 200 GeV");
0944     // ltx->DrawLatex(0.38, 0.735, "Tracklet : vertex + outer cluster");
0945     // c3 -> Print(Form("%s/folder_%s_advanced/dNdeta_2P_outer.pdf",mother_folder_directory.c_str(),file_name.c_str()));
0946     // pad_obj -> Clear();
0947
0948
0949     MemInfo_t test_aaa;
0950     gSystem->GetMemInfo(&test_aaa);
0951
0952     cout<<"Memory usage ? "<<test_aaa.fMemUsed<<endl;
0953
0954 }