Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:12:35 

0001 #include "EvtVtxZProtoTracklet.h"
0002 
0003 EvtVtxZProtoTracklet::EvtVtxZProtoTracklet(
0004     int process_id_in,
0005     int runnumber_in,
0006     int run_nEvents_in,
0007     std::string input_directory_in,
0008     std::string input_file_name_in,
0009     std::string output_directory_in,
0010     std::string output_file_name_suffix_in,
0011 
0012     std::pair<double, double> vertexXYIncm_in, // note : in cm
0013     bool IsFieldOn_in,
0014     bool IsDCACutApplied_in,
0015     std::pair<std::pair<double,double>,std::pair<double,double>> DeltaPhiCutInDegree_in, // note : in degree
0016     std::pair<std::pair<double,double>,std::pair<double,double>> DCAcutIncm_in, // note : in cm
0017     int ClusAdcCut_in,
0018     int ClusPhiSizeCut_in,
0019     
0020     bool PrintRecoDetails_in,
0021     bool DrawEvtVtxZ_in,
0022 
0023     bool RunInttBcoFullDiff_in,
0024     bool RunVtxZReco_in,
0025     bool RunTrackletPair_in,
0026     bool RunTrackletPairRotate_in,
0027     bool HaveGeoOffsetTag_in
0028 ):
0029     process_id(process_id_in),
0030     runnumber(runnumber_in),
0031     run_nEvents(run_nEvents_in),
0032     input_directory(input_directory_in),
0033     input_file_name(input_file_name_in),
0034     output_directory(output_directory_in),
0035     output_file_name_suffix(output_file_name_suffix_in),
0036 
0037     vertexXYIncm(vertexXYIncm_in),
0038     IsFieldOn(IsFieldOn_in),
0039     IsDCACutApplied(IsDCACutApplied_in),
0040     DeltaPhiCutInDegree(DeltaPhiCutInDegree_in),
0041     DCAcutIncm(DCAcutIncm_in),
0042     ClusAdcCut(ClusAdcCut_in),
0043     ClusPhiSizeCut(ClusPhiSizeCut_in),
0044 
0045     PrintRecoDetails(PrintRecoDetails_in),
0046     DrawEvtVtxZ(DrawEvtVtxZ_in),
0047     
0048     RunInttBcoFullDiff(RunInttBcoFullDiff_in),
0049     RunVtxZReco(RunVtxZReco_in),
0050     RunTrackletPair(RunTrackletPair_in),
0051     RunTrackletPairRotate(RunTrackletPairRotate_in),
0052     HaveGeoOffsetTag(HaveGeoOffsetTag_in)
0053 {
0054     
0055     PrepareOutPutFileName();
0056 
0057     for (int layer_i = B0L0_index; layer_i <= B1L1_index; layer_i++){
0058         double N_layer_ladder = (layer_i == B0L0_index || layer_i == B0L1_index) ? nLadder_inner : nLadder_outer;
0059         for (int phi_i = 0; phi_i < N_layer_ladder; phi_i++){
0060             geo_offset_map[Form("%i_%i", layer_i, phi_i)] = {0, 0, 0};
0061         }
0062     }
0063 
0064     // system(Form("if [ -f %s/completed/%s ]; then rm %s/completed/%s; fi;", output_directory.c_str(), output_filename.c_str(), output_directory.c_str(), output_filename.c_str()));
0065     // system(Form("cp %s/%s %s/%s", input_directory.c_str(), input_file_name.c_str(), output_directory.c_str(), output_filename.c_str()));
0066 
0067     PrepareRootFile();
0068     PrepareINTTvtxZ();
0069 
0070     track_gr = new TGraphErrors();
0071     fit_rz = new TF1("fit_rz","pol1",-1000,1000);
0072 
0073     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0074     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0075 
0076     out_eID_count = 0;
0077 
0078     if (DrawEvtVtxZ){
0079         PrepareEvtCanvas();
0080     }
0081 
0082 }
0083 
0084 void EvtVtxZProtoTracklet::PrepareOutPutFileName()
0085 {
0086     std::string job_index = std::to_string( process_id );
0087     int job_index_len = 5;
0088     job_index.insert(0, job_index_len - job_index.size(), '0');
0089 
0090     std::string runnumber_str = std::to_string( runnumber );
0091     if (runnumber != -1){
0092         int runnumber_str_len = 8;
0093         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0094     }
0095 
0096     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0097         output_file_name_suffix = "_" + output_file_name_suffix;
0098     }
0099 
0100     output_filename = "EvtVtxZProtoTracklet";
0101     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0102     output_filename += (IsFieldOn) ? "_FieldOn" : "_FieldOff";
0103     output_filename += (RunInttBcoFullDiff && runnumber != -1) ? "_BcoFullDiff" : "";
0104     output_filename += (RunVtxZReco) ? "_VtxZReco" : ""; 
0105     output_filename += (RunTrackletPair) ? "_Tracklet" : ""; 
0106     output_filename += (RunTrackletPairRotate) ? "_TrackletRotate" : ""; 
0107     output_filename += (HaveGeoOffsetTag) ? "_GeoOffset" : "";
0108     output_filename += output_file_name_suffix;
0109     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0110 }
0111 
0112 void EvtVtxZProtoTracklet::SetGeoOffset(std::map<std::string, std::vector<double>> input_geo_offset_map)
0113 {
0114     std::cout<<"In EvtVtxZProtoTracklet, reading geo offset map:"<<std::endl;
0115 
0116     geo_offset_map.clear();
0117     geo_offset_map = input_geo_offset_map;
0118 
0119     for (auto &pair : geo_offset_map)
0120     {
0121         std::cout<<pair.first<<" : {"<<pair.second[0]<<", "<<pair.second[1]<<", "<<pair.second[2]<<"}"<<std::endl;
0122     }
0123 }
0124 
0125 double EvtVtxZProtoTracklet::CheckGeoOffsetMap()
0126 {
0127     double sum = 0;
0128     for (auto &pair : geo_offset_map)
0129     {
0130         sum += fabs(pair.second[0]) + fabs(pair.second[1]) + fabs(pair.second[2]);
0131     }
0132     return sum;
0133 }
0134 
0135 std::map<std::string, int> EvtVtxZProtoTracklet::GetInputTreeBranches(TTree * m_tree_in)
0136 {
0137     std::map<std::string, int> branch_map;
0138     TObjArray * branch_list = m_tree_in -> GetListOfBranches();
0139     for (int i = 0; i < branch_list -> GetEntries(); i++)
0140     {
0141         TBranch * branch = dynamic_cast<TBranch*>(branch_list->At(i));
0142         branch_map[branch -> GetName()] = 1;
0143     }
0144     return branch_map;
0145 }
0146 
0147 void EvtVtxZProtoTracklet::PrepareRootFile()
0148 {
0149     file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0150     if (!file_in || file_in -> IsZombie() || file_in == nullptr) {
0151         std::cout << "Error: cannot open file: " << input_file_name << std::endl;
0152         exit(1);
0153     }
0154 
0155     tree_in = (TTree*)file_in -> Get("EventTree");
0156     // tree_in -> SetBranchStatus("*", 0);
0157     std::map<std::string, int> branch_map = GetInputTreeBranches(tree_in);
0158 
0159     // note : -------------------------------------------------------
0160     if(branch_map.find("TrkrHitRow") != branch_map.end())          {tree_in -> SetBranchStatus("TrkrHitRow", 0);}
0161     if(branch_map.find("TrkrHitColumn") != branch_map.end())       {tree_in -> SetBranchStatus("TrkrHitColumn", 0);}
0162     if(branch_map.find("TrkrHitLadderZId") != branch_map.end())    {tree_in -> SetBranchStatus("TrkrHitLadderZId", 0);}
0163     if(branch_map.find("TrkrHitLadderPhiId") != branch_map.end())  {tree_in -> SetBranchStatus("TrkrHitLadderPhiId", 0);}
0164     if(branch_map.find("TrkrHitLayer") != branch_map.end())        {tree_in -> SetBranchStatus("TrkrHitLayer", 0);}
0165     if(branch_map.find("TrkrHitADC") != branch_map.end())          {tree_in -> SetBranchStatus("TrkrHitADC", 0);}
0166     if(branch_map.find("TrkrHitX") != branch_map.end())            {tree_in -> SetBranchStatus("TrkrHitX", 0);}
0167     if(branch_map.find("TrkrHitY") != branch_map.end())            {tree_in -> SetBranchStatus("TrkrHitY", 0);}
0168     if(branch_map.find("TrkrHitZ") != branch_map.end())            {tree_in -> SetBranchStatus("TrkrHitZ", 0);}
0169     if(branch_map.find("TrkrHit_truthHit_x0") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_x0", 0);}
0170     if(branch_map.find("TrkrHit_truthHit_y0") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_y0", 0);}
0171     if(branch_map.find("TrkrHit_truthHit_z0") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_z0", 0);}
0172     if(branch_map.find("TrkrHit_truthHit_x1") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_x1", 0);}
0173     if(branch_map.find("TrkrHit_truthHit_y1") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_y1", 0);}
0174     if(branch_map.find("TrkrHit_truthHit_z1") != branch_map.end()) {tree_in -> SetBranchStatus("TrkrHit_truthHit_z1", 0);}
0175 
0176     if(branch_map.find("MBD_PMT_charge") != branch_map.end())      {tree_in -> SetBranchStatus("MBD_PMT_charge", 0);}
0177 
0178     if(branch_map.find("triggervec") != branch_map.end())                {tree_in -> SetBranchStatus("triggervec", 0);}
0179     // if(branch_map.find("firedTriggers") != branch_map.end())             {tree_in -> SetBranchStatus("firedTriggers", 0);}
0180     if(branch_map.find("firedTriggers_name") != branch_map.end())        {tree_in -> SetBranchStatus("firedTriggers_name", 0);}
0181     if(branch_map.find("firedTriggers_checkraw") != branch_map.end())    {tree_in -> SetBranchStatus("firedTriggers_checkraw", 0);}
0182     if(branch_map.find("firedTriggers_prescale") != branch_map.end())    {tree_in -> SetBranchStatus("firedTriggers_prescale", 0);}
0183     if(branch_map.find("firedTriggers_scalers") != branch_map.end())     {tree_in -> SetBranchStatus("firedTriggers_scalers", 0);}
0184     if(branch_map.find("firedTriggers_livescalers") != branch_map.end()) {tree_in -> SetBranchStatus("firedTriggers_livescalers", 0);}
0185     if(branch_map.find("firedTriggers_rawscalers") != branch_map.end())  {tree_in -> SetBranchStatus("firedTriggers_rawscalers", 0);}
0186 
0187     if (runnumber == -1){
0188         if(branch_map.find("PHG4Hit_x0") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_x0", 0);}
0189         if(branch_map.find("PHG4Hit_y0") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_y0", 0);}
0190         if(branch_map.find("PHG4Hit_z0") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_z0", 0);}
0191         if(branch_map.find("PHG4Hit_x1") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_x1", 0);}
0192         if(branch_map.find("PHG4Hit_y1") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_y1", 0);}
0193         if(branch_map.find("PHG4Hit_z1") != branch_map.end())     {tree_in -> SetBranchStatus("PHG4Hit_z1", 0);}
0194         if(branch_map.find("PHG4Hit_edep") != branch_map.end())   {tree_in -> SetBranchStatus("PHG4Hit_edep", 0);}
0195 
0196         if(branch_map.find("HepMCFSPrtl_Pt") != branch_map.end())    {tree_in -> SetBranchStatus("HepMCFSPrtl_Pt", 0);}
0197         if(branch_map.find("HepMCFSPrtl_Eta") != branch_map.end())   {tree_in -> SetBranchStatus("HepMCFSPrtl_Eta", 0);}
0198         if(branch_map.find("HepMCFSPrtl_Phi") != branch_map.end())   {tree_in -> SetBranchStatus("HepMCFSPrtl_Phi", 0);}
0199         if(branch_map.find("HepMCFSPrtl_E") != branch_map.end())     {tree_in -> SetBranchStatus("HepMCFSPrtl_E", 0);}
0200         if(branch_map.find("HepMCFSPrtl_PID") != branch_map.end())   {tree_in -> SetBranchStatus("HepMCFSPrtl_PID", 0);}
0201         if(branch_map.find("HepMCFSPrtl_prodx") != branch_map.end()) {tree_in -> SetBranchStatus("HepMCFSPrtl_prodx", 0);}
0202         if(branch_map.find("HepMCFSPrtl_prody") != branch_map.end()) {tree_in -> SetBranchStatus("HepMCFSPrtl_prody", 0);}
0203         if(branch_map.find("HepMCFSPrtl_prodz") != branch_map.end()) {tree_in -> SetBranchStatus("HepMCFSPrtl_prodz", 0);}
0204     }
0205     // note : -------------------------------------------------------
0206     
0207     
0208     // tree_in -> SetBranchStatus("event", 1);
0209     // tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0210     // tree_in -> SetBranchStatus("INTT_BCO", 1);
0211     // tree_in -> SetBranchStatus("ClusX", 1);
0212     // tree_in -> SetBranchStatus("ClusY", 1);
0213     // tree_in -> SetBranchStatus("ClusZ", 1);
0214     // tree_in -> SetBranchStatus("ClusLayer", 1);
0215     // tree_in -> SetBranchStatus("ClusLadderZId", 1);
0216     // tree_in -> SetBranchStatus("ClusLadderPhiId", 1);
0217     // tree_in -> SetBranchStatus("ClusAdc", 1);
0218     // tree_in -> SetBranchStatus("ClusPhiSize", 1);
0219     // if (!RunVtxZReco){
0220     //     tree_in -> SetBranchStatus("INTTvtxZ",1);
0221     //     tree_in -> SetBranchStatus("INTTvtxZError",1);
0222     // }
0223 
0224     // if (runnumber == -1) {
0225     //     tree_in -> SetBranchStatus("TruthPV_trig_z",1);
0226     // }
0227     
0228     ClusX = 0;
0229     ClusY = 0;
0230     ClusZ = 0;
0231     ClusLayer = 0;
0232     ClusLadderZId = 0;
0233     ClusLadderPhiId = 0;
0234     ClusAdc = 0;
0235     ClusPhiSize = 0;
0236     firedTriggers = 0;
0237 
0238     // tree_in -> SetBranchAddress("event", &event);
0239     if(branch_map.find("firedTriggers") != branch_map.end()) {
0240         tree_in -> SetBranchAddress("firedTriggers", &firedTriggers);
0241         m_withTrig = true;
0242     }
0243 
0244     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0245     tree_in -> SetBranchAddress("INTT_BCO", &INTT_BCO);
0246     tree_in -> SetBranchAddress("ClusX", &ClusX);
0247     tree_in -> SetBranchAddress("ClusY", &ClusY);
0248     tree_in -> SetBranchAddress("ClusZ", &ClusZ);
0249     tree_in -> SetBranchAddress("ClusLayer", &ClusLayer);
0250     tree_in -> SetBranchAddress("ClusLadderZId", &ClusLadderZId);
0251     tree_in -> SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
0252     tree_in -> SetBranchAddress("ClusAdc", &ClusAdc);
0253     tree_in -> SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0254     if (!RunVtxZReco){
0255         tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);
0256         tree_in -> SetBranchAddress("INTTvtxZError",&INTTvtxZError);
0257     }
0258     
0259     if (runnumber == -1) {
0260         tree_in -> SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
0261         tree_in -> SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
0262         tree_in -> SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0263         tree_in -> SetBranchAddress("NTruthVtx", &NTruthVtx);
0264     }
0265 
0266     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0267     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0268 
0269 
0270     // note : for the output file
0271     file_out = new TFile(Form("%s/%s",output_directory.c_str(),output_filename.c_str()),"recreate");
0272     tree_out = tree_in->CloneTree(0);
0273 
0274     // note : add new branches
0275     if (RunInttBcoFullDiff && runnumber != -1) {b_InttBcoFullDiff_next = tree_out -> Branch("InttBcoFullDiff_next", &out_InttBcoFullDiff_next);}
0276 
0277     if (RunVtxZReco){
0278         b_INTTvtxZ = tree_out -> Branch("INTTvtxZ", &out_INTTvtxZ);
0279         b_INTTvtxZError = tree_out -> Branch("INTTvtxZError", &out_INTTvtxZError);
0280         b_NgroupTrapezoidal = tree_out -> Branch("NgroupTrapezoidal", &out_NgroupTrapezoidal);
0281         b_NgroupCoarse = tree_out -> Branch("NgroupCoarse", &out_NgroupCoarse);
0282         b_TrapezoidalFitWidth = tree_out -> Branch("TrapezoidalFitWidth", &out_TrapezoidalFitWidth);
0283         b_TrapezoidalFWHM = tree_out -> Branch("TrapezoidalFWHM", &out_TrapezoidalFWHM);
0284 
0285         b_ClusEta_INTTz = tree_out -> Branch("ClusEta_INTTz", &out_ClusEta_INTTz);
0286         b_ClusEta_MBDz = tree_out -> Branch("ClusEta_MBDz", &out_ClusEta_MBDz);
0287         b_ClusPhi_AvgPV = tree_out -> Branch("ClusPhi_AvgPV", &out_ClusPhi_AvgPV);
0288         
0289         if (runnumber == -1){
0290             b_ClusEta_TrueXYZ = tree_out -> Branch("ClusEta_TrueXYZ", &out_ClusEta_TrueXYZ);
0291             b_ClusPhi_TrueXY = tree_out -> Branch("ClusPhi_TrueXY", &out_ClusPhi_TrueXY);
0292         }
0293 
0294         if (m_withTrig){
0295             b_MBDNSg2 = tree_out -> Branch("MBDNSg2", &out_MBDNSg2);
0296             b_MBDNSg2_vtxZ10cm = tree_out -> Branch("MBDNSg2_vtxZ10cm", &out_MBDNSg2_vtxZ10cm);
0297             b_MBDNSg2_vtxZ30cm = tree_out -> Branch("MBDNSg2_vtxZ30cm", &out_MBDNSg2_vtxZ30cm);
0298             b_MBDNSg2_vtxZ60cm = tree_out -> Branch("MBDNSg2_vtxZ60cm", &out_MBDNSg2_vtxZ60cm);
0299         }
0300 
0301         b_eID_count = tree_out -> Branch("eID_count", &out_eID_count);
0302     }
0303      
0304     if (RunTrackletPair) {b_evt_TrackletPair_vec = tree_out -> Branch("TrackletPair", &out_evt_TrackletPair_vec);}
0305     if (RunTrackletPairRotate) {b_evt_TrackletPairRotate_vec = tree_out -> Branch("TrackletPairRotate", &out_evt_TrackletPairRotate_vec);}
0306 
0307 }
0308 
0309 void EvtVtxZProtoTracklet::GetTriggerInfo()
0310 {
0311     std::vector<int> firedTriggers_vec =  *(firedTriggers);
0312     std::map<int, int> firedTriggers_map; firedTriggers_map.clear();
0313 
0314     for (int trig : firedTriggers_vec){
0315         firedTriggers_map[trig] = 1;
0316     }
0317 
0318     out_MBDNSg2 = (firedTriggers_map.find(index_MBDNSg2) != firedTriggers_map.end()) ? 1 : 0;
0319     out_MBDNSg2_vtxZ10cm = (firedTriggers_map.find(index_MBDNSg2_vtxZ10cm) != firedTriggers_map.end()) ? 1 : 0;
0320     out_MBDNSg2_vtxZ30cm = (firedTriggers_map.find(index_MBDNSg2_vtxZ30cm) != firedTriggers_map.end()) ? 1 : 0;
0321     out_MBDNSg2_vtxZ60cm = (firedTriggers_map.find(index_MBDNSg2_vtxZ60cm) != firedTriggers_map.end()) ? 1 : 0;
0322 }
0323 
0324 void EvtVtxZProtoTracklet::PrepareINTTvtxZ()
0325 {
0326     evt_possible_z = new TH1D("","evt_possible_z",50, evt_possible_z_range.first, evt_possible_z_range.second);
0327     evt_possible_z -> SetLineWidth(1);
0328     evt_possible_z -> GetXaxis() -> SetTitle("Z [cm]");
0329     evt_possible_z -> GetYaxis() -> SetTitle("Entries");
0330 
0331 
0332     line_breakdown_hist = new TH1D("line_breakdown_hist", "line_breakdown_hist", 2*zvtx_hist_Nbins+1, -1*(fine_bin_width*zvtx_hist_Nbins + fine_bin_width/2.), fine_bin_width*zvtx_hist_Nbins + fine_bin_width/2.);
0333     line_breakdown_hist -> SetLineWidth(1);
0334     line_breakdown_hist -> GetXaxis() -> SetTitle("Z [cm]");
0335     line_breakdown_hist -> GetYaxis() -> SetTitle("Entries [A.U.]");
0336     line_breakdown_hist -> GetXaxis() -> SetNdivisions(705);
0337 
0338     combination_zvtx_range_shape = new TH1D(
0339         "",
0340         "",
0341         line_breakdown_hist -> GetNbinsX(),
0342         line_breakdown_hist -> GetXaxis() -> GetXmin(),
0343         line_breakdown_hist -> GetXaxis() -> GetXmax()
0344     );
0345 
0346     gaus_fit = new TF1("gaus_fit",gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4);
0347     gaus_fit -> SetLineColor(2);
0348     gaus_fit -> SetLineWidth(1);
0349     gaus_fit -> SetNpx(1000);
0350 
0351     gaus_fit_vec.clear();
0352     for (int i = 0; i < number_of_gaus; i++) 
0353     {
0354         gaus_fit_vec.push_back(new TF1(Form("gaus_fit_vec_%i",i),gaus_func,evt_possible_z_range.first,evt_possible_z_range.second,4));
0355         // gaus_fit_vec[i] -> SetLineColor(TColor::GetColor(color_code[i].c_str()));
0356         gaus_fit_vec[i] -> SetLineColor(ROOT_color_code[i]);
0357         gaus_fit_vec[i] -> SetLineWidth(1);
0358         gaus_fit_vec[i] -> SetNpx(1000);
0359     }
0360 
0361     reco_INTTvtxZ = new TH1D("reco_INTTvtxZ","reco_INTTvtxZ;INTT vtxZ [cm];Entries",nVtxZ_bin,vtxZ_range.first,vtxZ_range.second);
0362 }
0363 
0364 void EvtVtxZProtoTracklet::PrepareClusterVec()
0365 {
0366     for (int clu_i = 0; clu_i < ClusX -> size(); clu_i++)
0367     {
0368         EvtVtxZProtoTracklet::clu_info this_clu;
0369 
0370         this_clu.adc = ClusAdc -> at(clu_i);
0371         this_clu.phi_size = ClusPhiSize -> at(clu_i);
0372         this_clu.sensorZID = ClusLadderZId -> at(clu_i);
0373         this_clu.ladderPhiID = ClusLadderPhiId -> at(clu_i);
0374         this_clu.layerID = ClusLayer -> at(clu_i);
0375 
0376         this_clu.index = clu_i;
0377 
0378         this_clu.x = ClusX -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][0];
0379         this_clu.y = ClusY -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][1];
0380         this_clu.z = ClusZ -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][2];
0381 
0382 
0383         std::vector<EvtVtxZProtoTracklet::clu_info>* p_evt_sPH_nocolumn_vec =
0384         (this_clu.layerID == 3 || this_clu.layerID == 4) ? (&evt_sPH_inner_nocolumn_vec) : (&evt_sPH_outer_nocolumn_vec);
0385         p_evt_sPH_nocolumn_vec -> push_back(this_clu);
0386 
0387         if (this_clu.adc <= ClusAdcCut || this_clu.phi_size > ClusPhiSizeCut) {continue;}
0388 
0389         std::vector<EvtVtxZProtoTracklet::clu_info>* p_evt_sPH_nocolumn_vec_PostCut =
0390         (this_clu.layerID == 3 || this_clu.layerID == 4) ? (&evt_sPH_inner_nocolumn_vec_PostCut) : (&evt_sPH_outer_nocolumn_vec_PostCut);
0391         p_evt_sPH_nocolumn_vec_PostCut -> push_back(this_clu);
0392     }
0393 }
0394 
0395 void EvtVtxZProtoTracklet::GetINTTvtxZ()
0396 {
0397 
0398     inner_clu_phi_map_PostCut.clear();
0399     outer_clu_phi_map_PostCut.clear();
0400     inner_clu_phi_map_PostCut = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
0401     outer_clu_phi_map_PostCut = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
0402 
0403     for (int inner_i = 0; inner_i < int(evt_sPH_inner_nocolumn_vec_PostCut.size()); inner_i++) {
0404         double Clus_InnerPhi_Offset = (evt_sPH_inner_nocolumn_vec_PostCut[inner_i].y - vertexXYIncm.second < 0) ? atan2(evt_sPH_inner_nocolumn_vec_PostCut[inner_i].y - vertexXYIncm.second, evt_sPH_inner_nocolumn_vec_PostCut[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(evt_sPH_inner_nocolumn_vec_PostCut[inner_i].y - vertexXYIncm.second, evt_sPH_inner_nocolumn_vec_PostCut[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0405         inner_clu_phi_map_PostCut[ int(Clus_InnerPhi_Offset) ].push_back({false,evt_sPH_inner_nocolumn_vec_PostCut[inner_i]});
0406     }
0407     for (int outer_i = 0; outer_i < int(evt_sPH_outer_nocolumn_vec_PostCut.size()); outer_i++) {
0408         double Clus_OuterPhi_Offset = (evt_sPH_outer_nocolumn_vec_PostCut[outer_i].y - vertexXYIncm.second < 0) ? atan2(evt_sPH_outer_nocolumn_vec_PostCut[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec_PostCut[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(evt_sPH_outer_nocolumn_vec_PostCut[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec_PostCut[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0409         outer_clu_phi_map_PostCut[ int(Clus_OuterPhi_Offset) ].push_back({false,evt_sPH_outer_nocolumn_vec_PostCut[outer_i]});
0410     }
0411 
0412   // note : prepare good tracklet
0413   for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0414   {
0415       // note : N cluster in this phi cell
0416       for (int inner_phi_clu_i = 0; inner_phi_clu_i < int(inner_clu_phi_map_PostCut[inner_phi_i].size()); inner_phi_clu_i++)
0417       {
0418           if (inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
0419 
0420           double Clus_InnerPhi_Offset = (inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second < 0) ? atan2(inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second, inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second, inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first) * (180./TMath::Pi());
0421 
0422           // todo: change the outer phi scan range
0423           // note : the outer phi index, -1, 0, 1
0424           // note : the outer phi index, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 for the scan test
0425           for (int scan_i = -5; scan_i < 6; scan_i++)
0426           {
0427               int true_scan_i = ((inner_phi_i + scan_i) < 0) ? 360 + (inner_phi_i + scan_i) : ((inner_phi_i + scan_i) > 359) ? (inner_phi_i + scan_i)-360 : inner_phi_i + scan_i;
0428 
0429               // note : N clusters in that outer phi cell
0430               for (int outer_phi_clu_i = 0; outer_phi_clu_i < int(outer_clu_phi_map_PostCut[true_scan_i].size()); outer_phi_clu_i++)
0431               {
0432                   if (outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].first == true) {continue;}
0433 
0434                   double Clus_OuterPhi_Offset = (outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second < 0) ? atan2(outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second, outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second, outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first) * (180./TMath::Pi());
0435                   double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0436                   
0437                   double delta_phi_correct = delta_phi; // note : this is for the further correction, potentially 
0438 
0439                   // note : ----------------------------------------------------------------------------------------------------------------------------------
0440                   // note : delta phi cut
0441                   if (!IsFieldOn){ // note : field off
0442                     if (delta_phi_correct <= DeltaPhiCutInDegree.first.first || delta_phi_correct >= DeltaPhiCutInDegree.first.second) {continue;}
0443                   }
0444                   // note : field ON
0445                   else{
0446                     if (
0447                           (delta_phi_correct > DeltaPhiCutInDegree.first.first && delta_phi_correct < DeltaPhiCutInDegree.first.second) == false && 
0448                           (delta_phi_correct > DeltaPhiCutInDegree.second.first && delta_phi_correct < DeltaPhiCutInDegree.second.second) == false
0449                        ) 
0450                        {continue;}
0451                   }
0452 
0453                   double DCA_sign = calculateAngleBetweenVectors(
0454                       outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.x, outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y,
0455                       inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.x, inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y,
0456                       vertexXYIncm.first, vertexXYIncm.second
0457                   );   
0458 
0459                   // note : ----------------------------------------------------------------------------------------------------------------------------------
0460                   // note : DCA cut
0461                   if (IsDCACutApplied)
0462                   {
0463                     if (!IsFieldOn){ // note : field off
0464                       if (DCA_sign <= DCAcutIncm.first.first || DCA_sign >= DCAcutIncm.first.second) {continue;}
0465                     }
0466                     else{
0467                       if (
0468                           (DCA_sign > DCAcutIncm.first.first && DCA_sign < DCAcutIncm.first.second) == false && 
0469                           (DCA_sign > DCAcutIncm.second.first && DCA_sign < DCAcutIncm.second.second) == false
0470                        ) 
0471                        {continue;}
0472                     }
0473                   }
0474 
0475                   // note : ----------------------------------------------------------------------------------------------------------------------------------
0476                   // note : good proto-tracklet
0477                   std::pair<double,double> z_range_info = Get_possible_zvtx( 
0478                       0., // get_radius(vertexXYIncm.first,vertexXYIncm.second), 
0479                       {
0480                         get_radius(inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first, inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second), 
0481                         inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.z,
0482                         double(inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.sensorZID)
0483                       }, // note : unsign radius
0484                       
0485                       {
0486                         get_radius(outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first, outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second), 
0487                         outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.z,
0488                         double(outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.sensorZID)
0489                       }  // note : unsign radius
0490                   );
0491 
0492                   // note : ----------------------------------------------------------------------------------------------------------------------------------
0493                   // note : check the coverage
0494                   if (evt_possible_z_range.first < z_range_info.first && z_range_info.first < evt_possible_z_range.second) {
0495                       evt_possible_z -> Fill(z_range_info.first);
0496 
0497                       // note : fill the line_breakdown histogram as well as a vector for the width determination
0498                       trapezoidal_line_breakdown( 
0499                           line_breakdown_hist, 
0500                           
0501                           // note : inner_r and inner_z
0502                           get_radius(inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first, inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second), 
0503                           inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.z,
0504                           inner_clu_phi_map_PostCut[inner_phi_i][inner_phi_clu_i].second.sensorZID,
0505 
0506                           // note : outer_r and outer_z
0507                           get_radius(outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first, outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second), 
0508                           outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.z,
0509                           outer_clu_phi_map_PostCut[true_scan_i][outer_phi_clu_i].second.sensorZID
0510                       );
0511                   }
0512                   
0513               }
0514           } // note : end of outer clu loop
0515       } 
0516 
0517   } // note : end of inner clu loop
0518 
0519   // todo : some of the events have the special behavior which is having the high entry in the both edges
0520   // note : prepare the vertex Z by INTT
0521   for (int bin_i = 0; bin_i < line_breakdown_hist->GetNbinsX(); bin_i++){
0522       if (line_breakdown_hist -> GetBinCenter(bin_i+1) < edge_rejection.first || line_breakdown_hist -> GetBinCenter(bin_i+1) > edge_rejection.second){
0523           line_breakdown_hist -> SetBinContent(bin_i+1,0);
0524       }
0525   }
0526 
0527   // note : no good proto-tracklet
0528   if (line_breakdown_hist->Integral() == 0) {
0529     return;
0530   }
0531 
0532   std::vector<double> N_group_info = find_Ngroup(evt_possible_z);
0533 
0534 
0535   std::vector<double> N_group_info_detail = find_Ngroup(line_breakdown_hist);
0536   double Half_N_group_half_width = fabs(N_group_info_detail[3] - N_group_info_detail[2]) / 2.;
0537   double line_breakdown_hist_max_content = line_breakdown_hist -> GetBinContent( line_breakdown_hist -> GetMaximumBin() );
0538   double line_breakdown_hist_max_center = line_breakdown_hist -> GetBinCenter( line_breakdown_hist -> GetMaximumBin() );
0539   
0540   // note : first fit is for the width, so apply the constraints on the Gaussian offset
0541   gaus_fit -> SetParameters(line_breakdown_hist_max_content, line_breakdown_hist_max_center, Half_N_group_half_width, 0);
0542   gaus_fit -> SetParLimits(0,0,100000);  // note : size 
0543   gaus_fit -> SetParLimits(2,0.5,1000);   // note : Width
0544   gaus_fit -> SetParLimits(3,-200,10000);   // note : offset
0545   // todo : the width fit range is 60% of the peak width, // todo : try to use single gaus to fit the distribution
0546   double width_fit_range_l = line_breakdown_hist_max_center - Half_N_group_half_width * 0.65;
0547   double width_fit_range_r = line_breakdown_hist_max_center + Half_N_group_half_width * 0.65;
0548   line_breakdown_hist -> Fit(gaus_fit, "NQ", "", width_fit_range_l, width_fit_range_r); 
0549   
0550   // std::cout<<"evnet_i: "<<", width fitting, SetParameters: "<<line_breakdown_hist_max_content<<" "<<line_breakdown_hist_max_center<<" "<<Half_N_group_half_width<<std::endl;
0551   // std::cout<<"evnet_i: "<<", width fitting, FitRange: "<<width_fit_range_l<<" "<<width_fit_range_r<<std::endl;
0552   // std::cout<<"evnet_i: "<<", width fitting, result: "<<gaus_fit -> GetParameter(0)<<" "<<gaus_fit -> GetParameter(1)<<" "<<gaus_fit -> GetParameter(2)<<" "<<gaus_fit -> GetParameter(3)<<std::endl;
0553   
0554 
0555   for (int fit_i = 0; fit_i < int(gaus_fit_vec.size()); fit_i++) // note : special_tag, uses a vector of gaus fit to determine the vertex Z
0556   {
0557       gaus_fit_vec[fit_i] -> SetParameters(line_breakdown_hist_max_content, line_breakdown_hist_max_center, Half_N_group_half_width, 0);
0558       gaus_fit_vec[fit_i] -> SetParLimits(0,0,100000);  // note : size 
0559       gaus_fit_vec[fit_i] -> SetParLimits(2,0.5,1000);   // note : Width
0560       gaus_fit_vec[fit_i] -> SetParLimits(3,-200,10000);   // note : offset
0561       double N_width_ratio = 0.2 + 0.15 * fit_i;
0562       double mean_fit_range_l = line_breakdown_hist_max_center - Half_N_group_half_width * N_width_ratio;
0563       double mean_fit_range_r = line_breakdown_hist_max_center + Half_N_group_half_width * N_width_ratio;
0564 
0565       line_breakdown_hist -> Fit(gaus_fit_vec[fit_i], "NQ", "", mean_fit_range_l, mean_fit_range_r);
0566       gaus_fit_vec[fit_i] -> SetRange(mean_fit_range_l, mean_fit_range_r);
0567 
0568       fit_mean_mean_vec.push_back(gaus_fit_vec[fit_i] -> GetParameter(1));
0569       fit_mean_reducedChi2_vec.push_back(gaus_fit_vec[fit_i] -> GetChisquare() / double(gaus_fit_vec[fit_i] -> GetNDF()));
0570       fit_mean_width_vec.push_back(fabs(gaus_fit_vec[fit_i] -> GetParameter(2)));
0571   }
0572 
0573   out_INTTvtxZ = vector_average(fit_mean_mean_vec);
0574   out_INTTvtxZError = vector_stddev(fit_mean_mean_vec);
0575 
0576   out_NgroupTrapezoidal = N_group_info_detail[0];
0577   out_NgroupCoarse = N_group_info[0];
0578   out_TrapezoidalFitWidth = gaus_fit -> GetParameter(2);
0579   out_TrapezoidalFWHM = (fabs(N_group_info_detail[3] - N_group_info_detail[2]) / 2.);
0580 
0581   // for (int hist_i = 0; hist_i < line_breakdown_hist->GetNbinsX(); hist_i++){
0582   //     // std::cout<<"("<<hist_i+1<<", "<<line_breakdown_hist->GetBinContent(hist_i+1)<<"), ";
0583   //     std::cout<<"{"<<hist_i+1<<", "<<line_breakdown_hist->GetBinContent(hist_i+1)<<"}, "<<std::endl;
0584   // }
0585   // std::cout<<std::endl;
0586 
0587   if (
0588         DrawEvtVtxZ && 
0589         (
0590             out_eID_count % 50 == 0 || 
0591             (out_INTTvtxZ - MBD_z_vtx) < -4.5 || 
0592             (out_INTTvtxZ - MBD_z_vtx) > 3.0 || 
0593             (out_TrapezoidalFitWidth < 0.9 && out_TrapezoidalFWHM < 6)
0594         )
0595     )
0596   {
0597     line_breakdown_hist_zoomin = (TH1D*)line_breakdown_hist -> Clone("line_breakdown_hist_zoomin");
0598 
0599     line_breakdown_hist_zoomin -> GetXaxis() -> SetRangeUser(
0600         gaus_fit_vec.front() -> GetParameter(1) - Half_N_group_half_width,
0601         gaus_fit_vec.front() -> GetParameter(1) + Half_N_group_half_width
0602     );
0603     line_breakdown_hist_zoomin -> SetMinimum( gaus_fit_vec.front() -> GetParameter(0) * 0.5);
0604     line_breakdown_hist_zoomin -> SetMaximum( gaus_fit_vec.front() -> GetParameter(0) * 1.5);
0605 
0606     c1 -> cd();
0607     pad_EvtZDist -> Draw();
0608     pad_EvtZDist -> cd();
0609 
0610     line_breakdown_hist -> SetMinimum(0);
0611     line_breakdown_hist -> SetMaximum(line_breakdown_hist -> GetBinContent( line_breakdown_hist->GetMaximumBin() ) * 1.65);
0612     line_breakdown_hist -> Draw("hist");
0613     for (int fit_i = 0; fit_i < gaus_fit_vec.size(); fit_i++){
0614         gaus_fit_vec[gaus_fit_vec.size() - fit_i - 1] -> Draw("l same");
0615     }
0616 
0617     draw_text -> DrawLatex(0.2, 0.9, Form("Event ID: %i", out_eID_count));
0618     draw_text -> DrawLatex(0.2, 0.86, Form("NClusGood: %i", evt_sPH_inner_nocolumn_vec_PostCut.size() + evt_sPH_outer_nocolumn_vec_PostCut.size()));
0619     draw_text -> DrawLatex(0.2, 0.82, Form("Reco. vtx Z: %.3f cm, StdDev: %.3f cm", out_INTTvtxZ, out_INTTvtxZError));
0620     draw_text -> DrawLatex(0.2, 0.78, Form("MBD_z_vtx: %.3f cm", MBD_z_vtx));
0621     if (runnumber == -1) {
0622         draw_text -> DrawLatex(0.2, 0.74, Form("True vtx Z: %.3f cm",TruthPV_trig_z));
0623         draw_text -> DrawLatex(0.2, 0.70, Form("Diff w.r.t TrueZ: %.3f cm", out_INTTvtxZ - TruthPV_trig_z));
0624     }
0625 
0626     INTTvtxZ_EvtDisplay_file_out -> cd();
0627     c1 -> Write(Form("eID_%d_EvtZDist", out_eID_count));
0628 
0629     c1 -> cd();
0630     pad_ZoomIn_EvtZDist -> Draw();
0631     pad_ZoomIn_EvtZDist -> cd();
0632     line_breakdown_hist_zoomin -> Draw("hist");
0633     for (int fit_i = 0; fit_i < gaus_fit_vec.size(); fit_i++){
0634         gaus_fit_vec[gaus_fit_vec.size() - fit_i - 1] -> Draw("l same");
0635     }
0636 
0637     INTTvtxZ_EvtDisplay_file_out -> cd();
0638     c1 -> Write(Form("eID_%d_EvtZDist_ZoomIn", out_eID_count));
0639 
0640     std::cout<<"Printed: eID: "<<out_eID_count<<" NClusGood: "<<evt_sPH_inner_nocolumn_vec_PostCut.size() + evt_sPH_outer_nocolumn_vec_PostCut.size()<<", NClusAll: "<<evt_sPH_inner_nocolumn_vec.size() + evt_sPH_outer_nocolumn_vec.size()<<", INTTvtxZ : "<<out_INTTvtxZ<<" INTTvtxZError : "<<out_INTTvtxZError<<" NgroupTrapezoidal : "<<out_NgroupTrapezoidal<<" NgroupCoarse : "<<out_NgroupCoarse<<" TrapezoidalFitWidth : "<<out_TrapezoidalFitWidth<<" TrapezoidalFWHM : "<<out_TrapezoidalFWHM<<std::endl;
0641 
0642     pad_ZoomIn_EvtZDist -> Clear();
0643     pad_EvtZDist -> Clear();
0644   }
0645 
0646   // note : print everything 
0647   if (PrintRecoDetails && out_eID_count % 200 == 0){
0648     std::cout<<"eID: "<<out_eID_count<<" NClusGood: "<<evt_sPH_inner_nocolumn_vec_PostCut.size() + evt_sPH_outer_nocolumn_vec_PostCut.size()<<", NClusAll: "<<evt_sPH_inner_nocolumn_vec.size() + evt_sPH_outer_nocolumn_vec.size()<<", INTTvtxZ : "<<out_INTTvtxZ<<" INTTvtxZError : "<<out_INTTvtxZError<<" NgroupTrapezoidal : "<<out_NgroupTrapezoidal<<" NgroupCoarse : "<<out_NgroupCoarse<<" TrapezoidalFitWidth : "<<out_TrapezoidalFitWidth<<" TrapezoidalFWHM : "<<out_TrapezoidalFWHM<<std::endl;
0649   }
0650 
0651   return;
0652 }
0653 
0654 std::pair<double,double> EvtVtxZProtoTracklet::rotatePoint(double x, double y)
0655 {
0656     // Convert the rotation angle from degrees to radians
0657     double rotation = rotate_phi_angle; // todo rotation is here
0658     double angleRad = rotation * M_PI / 180.0;
0659 
0660     // Perform the rotation
0661     double xOut = x * cos(angleRad) - y * sin(angleRad);
0662     double yOut = x * sin(angleRad) + y * cos(angleRad);
0663 
0664     xOut = (fabs(xOut) < 0.0000001) ? 0 : xOut;
0665     yOut = (fabs(yOut) < 0.0000001) ? 0 : yOut;
0666 
0667     // cout<<"Post rotation: "<<xOut<<" "<<yOut<<endl;
0668 
0669     return {xOut,yOut};
0670 }
0671 
0672 std::vector<EvtVtxZProtoTracklet::clu_info> EvtVtxZProtoTracklet::GetRoatedClusterVec(std::vector<EvtVtxZProtoTracklet::clu_info> input_cluster_vec)
0673 {
0674     std::vector<EvtVtxZProtoTracklet::clu_info> output_cluster_vec; output_cluster_vec.clear();
0675 
0676     for (EvtVtxZProtoTracklet::clu_info this_clu : input_cluster_vec)
0677     {
0678         std::pair<double,double> rotated_xy = rotatePoint(this_clu.x, this_clu.y);
0679      
0680         EvtVtxZProtoTracklet::clu_info rotated_clu = this_clu;
0681         rotated_clu.x = rotated_xy.first;
0682         rotated_clu.y = rotated_xy.second;
0683         output_cluster_vec.push_back(rotated_clu);
0684     }
0685 
0686     return output_cluster_vec;
0687 }
0688 
0689 void EvtVtxZProtoTracklet::GetClusUpdated()
0690 {
0691 
0692     for (int clu_i = 0; clu_i < ClusX -> size(); clu_i++)
0693     {
0694         if (MBD_z_vtx == MBD_z_vtx)
0695         {
0696             out_ClusEta_MBDz.push_back(
0697                 get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, MBD_z_vtx}, {ClusX -> at(clu_i), ClusY -> at(clu_i), ClusZ -> at(clu_i)})
0698             );
0699         }
0700 
0701         if (temp_INTTvtxZ == temp_INTTvtxZ && temp_INTTvtxZError == temp_INTTvtxZError)
0702         {
0703             out_ClusEta_INTTz.push_back(
0704                 get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, temp_INTTvtxZ}, {ClusX -> at(clu_i), ClusY -> at(clu_i), ClusZ -> at(clu_i)})
0705             );
0706         }
0707 
0708         if (runnumber == -1)
0709         {
0710             out_ClusEta_TrueXYZ.push_back(
0711                 get_clu_eta({TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z}, {ClusX -> at(clu_i), ClusY -> at(clu_i), ClusZ -> at(clu_i)})
0712             );
0713 
0714             out_ClusPhi_TrueXY.push_back(
0715                 atan2(ClusY -> at(clu_i) - TruthPV_trig_y, ClusX -> at(clu_i) - TruthPV_trig_x)
0716             );
0717         }
0718 
0719         out_ClusPhi_AvgPV.push_back(
0720             atan2(ClusY -> at(clu_i) - vertexXYIncm.second, ClusX -> at(clu_i) - vertexXYIncm.first)
0721         );
0722     }
0723     
0724 }
0725 
0726 void EvtVtxZProtoTracklet::GetTrackletPair(std::vector<pair_str> &input_TrackletPair_vec, bool isRotated)
0727 {
0728 
0729     inner_clu_phi_map.clear();
0730     outer_clu_phi_map.clear();
0731     inner_clu_phi_map = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
0732     outer_clu_phi_map = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
0733 
0734     if (temp_INTTvtxZ != temp_INTTvtxZ || temp_INTTvtxZError != temp_INTTvtxZError) {return;}
0735 
0736     std::vector<EvtVtxZProtoTracklet::clu_info> temp_evt_sPH_inner_nocolumn_vec = (isRotated) ? GetRoatedClusterVec(evt_sPH_inner_nocolumn_vec) : evt_sPH_inner_nocolumn_vec;
0737 
0738     for (int inner_i = 0; inner_i < int(temp_evt_sPH_inner_nocolumn_vec.size()); inner_i++) {
0739       double Clus_InnerPhi_Offset = (temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second < 0) ? atan2(temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second, temp_evt_sPH_inner_nocolumn_vec[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second, temp_evt_sPH_inner_nocolumn_vec[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0740       inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_evt_sPH_inner_nocolumn_vec[inner_i]});
0741     }
0742     for (int outer_i = 0; outer_i < int(evt_sPH_outer_nocolumn_vec.size()); outer_i++) {
0743         double Clus_OuterPhi_Offset = (evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second < 0) ? atan2(evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0744         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,evt_sPH_outer_nocolumn_vec[outer_i]});
0745     }
0746 
0747 
0748     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0749     {
0750         // note : N cluster in this phi cell
0751         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
0752         {
0753             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
0754 
0755             double Clus_InnerPhi_Offset_radian = atan2(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first);
0756             double Clus_InnerPhi_Offset = (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second < 0) ? Clus_InnerPhi_Offset_radian * (180./TMath::Pi()) + 360 : Clus_InnerPhi_Offset_radian * (180./TMath::Pi());
0757 
0758             // todo: change the outer phi scan range
0759             // note : the outer phi index, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5
0760             for (int scan_i = -5; scan_i < 6; scan_i++)
0761             {
0762                 int true_scan_i = ((inner_phi_i + scan_i) < 0) ? 360 + (inner_phi_i + scan_i) : ((inner_phi_i + scan_i) > 359) ? (inner_phi_i + scan_i)-360 : inner_phi_i + scan_i;
0763 
0764                 // note : N clusters in that outer phi cell
0765                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
0766                 {
0767                     if (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].first == true) {continue;}
0768 
0769                     double Clus_OuterPhi_Offset_radian = atan2(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first);
0770                     double Clus_OuterPhi_Offset = (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second < 0) ? Clus_OuterPhi_Offset_radian * (180./TMath::Pi()) + 360 : Clus_OuterPhi_Offset_radian * (180./TMath::Pi());
0771                     double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0772                     // note : degree to radian 
0773                     double delta_phi_radian = delta_phi * (TMath::Pi() / 180.);
0774                     
0775                     // if (fabs(delta_phi) > 5.72) {continue;}
0776                     // if (fabs(delta_phi) > 3.5) {continue;}
0777 
0778                     double inner_clu_eta = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, temp_INTTvtxZ},{inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z});
0779                     double outer_clu_eta = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, temp_INTTvtxZ},{outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z});
0780                     double delta_eta = inner_clu_eta - outer_clu_eta;
0781 
0782                     std::pair<double,double> z_range_info = Get_possible_zvtx( 
0783                         0., // get_radius(vertexXYIncm.first,vertexXYIncm.second), 
0784                         {
0785                             get_radius(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second), 
0786                             inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z,
0787                             double(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.sensorZID)
0788                         }, // note : unsign radius
0789                         
0790                         {
0791                             get_radius(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second), 
0792                             outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z,
0793                             double(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.sensorZID)
0794                         }  // note : unsign radius
0795                     );
0796 
0797                     // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0798                     if (z_range_info.first - z_range_info.second > temp_INTTvtxZ + temp_INTTvtxZError || z_range_info.first + z_range_info.second < temp_INTTvtxZ - temp_INTTvtxZError) {continue;}
0799 
0800                     // note : do the fill here (find the best match outer cluster with the inner cluster )
0801                     std::pair<double,double> Get_eta_pair = Get_eta(
0802                         {0., temp_INTTvtxZ,temp_INTTvtxZError},
0803                         {
0804                             get_radius(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x - vertexXYIncm.first, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y - vertexXYIncm.second), 
0805                             inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z,
0806                             double(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.sensorZID)
0807                         },
0808                         {
0809                             get_radius(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x - vertexXYIncm.first, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y - vertexXYIncm.second), 
0810                             outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z,
0811                             double(outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.sensorZID)
0812                         }
0813                     );  
0814 
0815                     pair_str temp_pair_str;
0816                     temp_pair_str.delta_phi = delta_phi_radian;
0817                     temp_pair_str.delta_eta = delta_eta;
0818                     temp_pair_str.pair_eta_num = (inner_clu_eta + outer_clu_eta) / 2.;
0819                     temp_pair_str.pair_eta_fit = Get_eta_pair.second;
0820 
0821                     // temp_pair_str.inner_phi_id = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.ladderPhiID;
0822                     // temp_pair_str.inner_layer_id = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.layerID;
0823                     temp_pair_str.inner_zid = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.sensorZID;
0824                     temp_pair_str.inner_phi_size = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.phi_size;
0825                     temp_pair_str.inner_adc = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.adc;
0826                     temp_pair_str.inner_x = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.x;
0827                     temp_pair_str.inner_y = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.y;
0828                     // temp_pair_str.inner_z = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z;
0829                     // temp_pair_str.inner_phi = Clus_InnerPhi_Offset_radian;
0830                     // temp_pair_str.inner_eta = inner_clu_eta;
0831                     temp_pair_str.inner_index = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.index;
0832                     
0833                     // temp_pair_str.outer_phi_id = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.ladderPhiID;
0834                     // temp_pair_str.outer_layer_id = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.layerID;
0835                     temp_pair_str.outer_zid = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.sensorZID;
0836                     temp_pair_str.outer_phi_size = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.phi_size;
0837                     temp_pair_str.outer_adc = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.adc;
0838                     temp_pair_str.outer_x = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.x;
0839                     temp_pair_str.outer_y = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.y;
0840                     // temp_pair_str.outer_z = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z;
0841                     // temp_pair_str.outer_phi = Clus_OuterPhi_Offset_radian;
0842                     // temp_pair_str.outer_eta = outer_clu_eta;
0843                     temp_pair_str.outer_index = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.index;
0844 
0845                     input_TrackletPair_vec.push_back(temp_pair_str);
0846 
0847                 }
0848             }
0849         }
0850     }
0851 
0852 
0853 }
0854 
0855 void EvtVtxZProtoTracklet::FillRecoINTTVtxZH1D(int event_index)
0856 {
0857 
0858     if (event_index % 20 == 0)
0859     {
0860         std::cout<<"In Filling RecoZ: "<<event_index<<", "<<MBD_z_vtx<<", "<<temp_INTTvtxZ<<", diff: "<<MBD_z_vtx - temp_INTTvtxZ<<std::endl;
0861     }
0862 
0863     // =======================================================================================================================================================
0864     // note : hard cut
0865 
0866     // note : for data
0867     if (runnumber != -1 && out_InttBcoFullDiff_next <= cut_InttBcoFullDIff_next) {return;}
0868     if (runnumber != -1 && out_MBDNSg2 != 1) {return;} // todo: assume MC no trigger
0869 
0870     // note : for MC
0871     // if (runnumber == -1 && NTruthVtx != 1) {return;}
0872 
0873     // note : both data and MC
0874     if (is_min_bias != 1) {return;}
0875     if (MBD_z_vtx != MBD_z_vtx) {return;}
0876     if (MBD_centrality != MBD_centrality) {return;}
0877     if (MBD_centrality < 0 || MBD_centrality > 100) {return;}
0878     if (temp_INTTvtxZ != temp_INTTvtxZ) {return;}
0879     if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {return;} // todo: the hard cut 60 cm 
0880 
0881     if (MBD_centrality > cut_MBD_centrality) {return;} // note : 1 - 70, for example 
0882 
0883     // =======================================================================================================================================================
0884     // note : optional cut
0885     if ((MBD_z_vtx - temp_INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - temp_INTTvtxZ > cut_vtxZDiff.second)) {return;}
0886     if ((out_TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || out_TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){return;}
0887     if ((out_TrapezoidalFWHM < cut_TrapezoidalFWHM.first || out_TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){return;}
0888 
0889     // =======================================================================================================================================================
0890 
0891     reco_INTTvtxZ -> Fill(temp_INTTvtxZ);
0892 }
0893 
0894 void EvtVtxZProtoTracklet::MainProcess()
0895 {
0896     if (HaveGeoOffsetTag == true && CheckGeoOffsetMap() <= 0) {
0897         std::cout<<"The HaveGeoOffsetTag is true, but the GeoOffsetMap is empty, please check the GeoOffsetMap"<<std::endl;
0898         exit(1);
0899     }
0900 
0901     if (HaveGeoOffsetTag == false && CheckGeoOffsetMap() > 0.0001) {
0902         std::cout<<"The HaveGeoOffsetTag is false, but the GeoOffsetMap has some non-zero numbers, please check the GeoOffsetMap"<<std::endl;
0903         exit(1);
0904     }
0905 
0906     for (int i = 0; i < run_nEvents; i++){
0907         tree_in -> GetEntry(i);
0908 
0909         if (i % 20 == 0) {
0910             std::cout<<"In EvtVtxZProtoTracklet, Processing event: "<<i<<std::endl;
0911         }
0912 
0913         EvtCleanUp();
0914 
0915         PrepareClusterVec();
0916 
0917         if (RunVtxZReco){
0918             GetINTTvtxZ();
0919             temp_INTTvtxZ = out_INTTvtxZ;
0920             temp_INTTvtxZError = out_INTTvtxZError;
0921 
0922             GetClusUpdated();
0923 
0924             if (m_withTrig) {GetTriggerInfo();}
0925         }
0926         else {
0927             temp_INTTvtxZ = INTTvtxZ;
0928             temp_INTTvtxZError = INTTvtxZError;
0929         }
0930 
0931         if (RunTrackletPair){
0932             GetTrackletPair(out_evt_TrackletPair_vec, false);
0933         }
0934 
0935         if (RunTrackletPairRotate){
0936             GetTrackletPair(out_evt_TrackletPairRotate_vec, true);
0937         }
0938 
0939         if (RunInttBcoFullDiff && runnumber != -1){
0940             if (i != run_nEvents -1){
0941                 ULong_t this_InttBcoFull = INTT_BCO;
0942 
0943                 tree_in -> GetEntry(i+1);
0944                 ULong_t next_InttBcoFull = INTT_BCO;
0945                 
0946                 out_InttBcoFullDiff_next = next_InttBcoFull - this_InttBcoFull;
0947 
0948                 tree_in -> GetEntry(i);
0949             }
0950             else {
0951                 out_InttBcoFullDiff_next = -1;
0952             }
0953         }
0954         
0955         FillRecoINTTVtxZH1D(i);
0956 
0957         tree_out -> Fill();
0958         // if (RunInttBcoFullDiff) {b_InttBcoFullDiff_next -> Fill();}
0959         // if (RunVtxZReco){
0960         //     b_INTTvtxZ -> Fill();
0961         //     b_INTTvtxZError -> Fill();
0962         //     b_NgroupTrapezoidal -> Fill();
0963         //     b_NgroupCoarse -> Fill();
0964         //     b_TrapezoidalFitWidth -> Fill();
0965         //     b_TrapezoidalFWHM -> Fill();
0966         // }
0967         // if (RunTrackletPair) {b_evt_TrackletPair_vec -> Fill();}
0968         // if (RunTrackletPairRotate) {b_evt_TrackletPairRotate_vec -> Fill();}
0969 
0970         out_eID_count++;
0971     }
0972 
0973     // note : if you use the update function, then you need to fill the rest of the events
0974     // // note : if the number of events is less than the run_nEvents, then fill the rest of the events
0975     // EvtCleanUp();
0976     // for (int i = run_nEvents; i < tree_in -> GetEntries(); i++){
0977     //     tree_in -> GetEntry(i);
0978     //     if (RunInttBcoFullDiff) {b_InttBcoFullDiff_next -> Fill();}
0979     //     if (RunVtxZReco){
0980     //         b_INTTvtxZ -> Fill();
0981     //         b_INTTvtxZError -> Fill();
0982     //         b_NgroupTrapezoidal -> Fill();
0983     //         b_NgroupCoarse -> Fill();
0984     //         b_TrapezoidalFitWidth -> Fill();
0985     //         b_TrapezoidalFWHM -> Fill();
0986     //     }
0987     //     if (RunTrackletPair) {b_evt_TrackletPair_vec -> Fill();}
0988     //     if (RunTrackletPairRotate) {b_evt_TrackletPairRotate_vec -> Fill();}
0989 
0990     //     // out_eID_count++;
0991     // }
0992 }
0993 
0994 void EvtVtxZProtoTracklet::EvtCleanUp()
0995 {
0996     evt_sPH_inner_nocolumn_vec_PostCut.clear();
0997     evt_sPH_outer_nocolumn_vec_PostCut.clear();
0998     inner_clu_phi_map_PostCut.clear();
0999     outer_clu_phi_map_PostCut.clear();
1000     inner_clu_phi_map_PostCut = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
1001     outer_clu_phi_map_PostCut = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
1002     
1003     evt_possible_z -> Reset("ICESM");
1004     line_breakdown_hist -> Reset("ICESM");
1005 
1006     out_INTTvtxZ = std::nan("");
1007     out_INTTvtxZError = std::nan("");
1008     out_NgroupTrapezoidal = std::nan("");
1009     out_NgroupCoarse = std::nan("");
1010     out_TrapezoidalFitWidth = std::nan("");
1011     out_TrapezoidalFWHM = std::nan("");
1012 
1013     temp_INTTvtxZ = std::nan("");
1014     temp_INTTvtxZError = std::nan("");
1015 
1016     fit_mean_mean_vec.clear();
1017     fit_mean_reducedChi2_vec.clear();
1018     fit_mean_width_vec.clear();
1019 
1020     out_ClusEta_INTTz.clear();
1021     out_ClusEta_MBDz.clear();
1022     out_ClusEta_TrueXYZ.clear();
1023     out_ClusPhi_AvgPV.clear();
1024     out_ClusPhi_TrueXY.clear();
1025 
1026 
1027     evt_sPH_inner_nocolumn_vec.clear();
1028     evt_sPH_outer_nocolumn_vec.clear();
1029     inner_clu_phi_map.clear();
1030     outer_clu_phi_map.clear();
1031     inner_clu_phi_map = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
1032     outer_clu_phi_map = std::vector<std::vector<std::pair<bool,EvtVtxZProtoTracklet::clu_info>>>(360);
1033 
1034     out_evt_TrackletPair_vec.clear();
1035     out_evt_TrackletPairRotate_vec.clear();
1036 
1037     out_InttBcoFullDiff_next = -1;
1038 
1039 }
1040 
1041 void EvtVtxZProtoTracklet::EndRun(int close_file_in)
1042 {   
1043     file_out -> cd();
1044     tree_out -> Write();
1045     reco_INTTvtxZ -> Write();
1046     file_out -> Close();
1047     
1048     if (DrawEvtVtxZ && INTTvtxZ_EvtDisplay_file_out != nullptr)
1049     {
1050         INTTvtxZ_EvtDisplay_file_out -> cd();
1051         INTTvtxZ_EvtDisplay_file_out -> Close();
1052     }
1053 
1054 
1055     if (close_file_in == 1){
1056         file_in -> cd();
1057         file_in -> Close();
1058     }
1059 }
1060 
1061 double EvtVtxZProtoTracklet::get_radius(double x, double y)
1062 {
1063     return sqrt(pow(x,2)+pow(y,2));
1064 }
1065 
1066 // note : angle_1 = inner clu phi
1067 // note: angle_2 = outer clu phi
1068 double EvtVtxZProtoTracklet::get_delta_phi(double angle_1, double angle_2)
1069 {
1070     std::vector<double> vec_abs = {std::fabs(angle_1 - angle_2), std::fabs(angle_1 - angle_2 + 360), std::fabs(angle_1 - angle_2 - 360)};
1071     std::vector<double> vec = {(angle_1 - angle_2), (angle_1 - angle_2 + 360), (angle_1 - angle_2 - 360)};
1072     return vec[std::distance(vec_abs.begin(), std::min_element(vec_abs.begin(),vec_abs.end()))];
1073 }
1074 
1075 double EvtVtxZProtoTracklet::calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
1076     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
1077     double vector1X = x2 - x1;
1078     double vector1Y = y2 - y1;
1079 
1080     double vector2X = targetX - x1;
1081     double vector2Y = targetY - y1;
1082 
1083     // Calculate the cross product of vector_1 and vector_2 (z-component)
1084     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
1085     
1086     // cout<<" crossProduct : "<<crossProduct<<endl;
1087 
1088     // Calculate the magnitudes of vector_1 and vector_2
1089     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
1090     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
1091 
1092     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
1093     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
1094 
1095     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
1096     // Convert the angle from radians to degrees and return it
1097     double angleInDegrees __attribute__((unused)) = angleInRadians * 180.0 / M_PI;
1098     
1099     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
1100     double angleInDegrees_new __attribute__((unused)) = angleInRadians_new * 180.0 / M_PI;
1101     
1102     // cout<<"angle : "<<angleInDegrees_new<<endl;
1103 
1104     double DCA_distance = sin(angleInRadians_new) * magnitude2;
1105 
1106     return DCA_distance;
1107 }
1108 
1109 double EvtVtxZProtoTracklet::Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
1110 {
1111     if ( fabs(p0x - p1x) < 0.000001 ){ // note : the line is vertical (if z is along the x axis)
1112         return p0x;
1113     }
1114     else {
1115         double slope = (p1y - p0y) / (p1x - p0x);
1116         double yIntercept = p0y - slope * p0x;
1117         double xCoordinate = (given_y - yIntercept) / slope;
1118         return xCoordinate;
1119     }
1120 }
1121 
1122 std::pair<double,double> EvtVtxZProtoTracklet::Get_possible_zvtx(double rvtx, std::vector<double> p0, std::vector<double> p1) // note : inner p0, outer p1, vector {r,z, zid}, -> {y,x}
1123 {
1124     std::vector<double> p0_z_edge = { ( p0[2] == typeA_sensorZID[0] || p0[2] == typeA_sensorZID[1] ) ? p0[1] - typeA_sensor_half_length_incm : p0[1] - typeB_sensor_half_length_incm, ( p0[2] == typeA_sensorZID[0] || p0[2] == typeA_sensorZID[1] ) ? p0[1] + typeA_sensor_half_length_incm : p0[1] + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1125     std::vector<double> p1_z_edge = { ( p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1] ) ? p1[1] - typeA_sensor_half_length_incm : p1[1] - typeB_sensor_half_length_incm, ( p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1] ) ? p1[1] + typeA_sensor_half_length_incm : p1[1] + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1126 
1127     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
1128     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
1129 
1130     double mid_point = (edge_first + edge_second) / 2.;
1131     double possible_width = fabs(edge_first - edge_second) / 2.;
1132 
1133     return {mid_point, possible_width}; // note : first : mid point, second : width
1134 
1135 }
1136 
1137 void EvtVtxZProtoTracklet::line_breakdown(TH1D* hist_in, std::pair<double,double> line_range)
1138 { 
1139     int first_bin = int((line_range.first  - hist_in->GetXaxis()->GetXmin()) /  hist_in->GetBinWidth(1)) + 1;
1140     int last_bin  = int((line_range.second - hist_in->GetXaxis()->GetXmin()) /  hist_in->GetBinWidth(1)) + 1;
1141     
1142     if (first_bin < 1) {first_bin = 0;}
1143     else if (first_bin > hist_in -> GetNbinsX()) {first_bin = hist_in -> GetNbinsX() + 1;}
1144 
1145     if (last_bin < 1) {last_bin = 0;}
1146     else if (last_bin > hist_in -> GetNbinsX()) {last_bin = hist_in -> GetNbinsX() + 1;}
1147 
1148     // double fill_weight = 1./fabs(line_range.second - line_range.first); // note : the entry is weitghted by the width of the line, by Akiba san's suggestion // special_tag cancel
1149     double fill_weight = 1.; // note : the weight is 1, it's for testing the trapezoidal method // special_tag
1150 
1151     // cout<<"Digitize the bin : "<<first_bin<<" "<<last_bin<<endl;
1152 
1153     // note : if first:last = (0:0) or (N+1:N+1) -> the subtraction of them euqals to zero.
1154     for (int i = 0; i < (last_bin - first_bin) + 1; i++){
1155         hist_in -> SetBinContent(first_bin + i,  hist_in -> GetBinContent(first_bin + i) + fill_weight ); // note : special_tag
1156         // line_breakdown_vec.push_back(hist_in -> GetBinCenter(first_bin + i));
1157     }
1158 
1159 }
1160 
1161 // note : for each combination
1162 // note : corrected inner_r and outer_r
1163 void EvtVtxZProtoTracklet::trapezoidal_line_breakdown(TH1D * hist_in, double inner_r, double inner_z, int inner_zid, double outer_r, double outer_z, int outer_zid)
1164 {
1165     combination_zvtx_range_shape -> Reset("ICESM");
1166 
1167     std::vector<double> inner_edge = { ( inner_zid == typeA_sensorZID[0] || inner_zid == typeA_sensorZID[1] ) ? inner_z - typeA_sensor_half_length_incm : inner_z - typeB_sensor_half_length_incm, ( inner_zid == typeA_sensorZID[0] || inner_zid == typeA_sensorZID[1] ) ? inner_z + typeA_sensor_half_length_incm : inner_z + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1168     std::vector<double> outer_edge = { ( outer_zid == typeA_sensorZID[0] || outer_zid == typeA_sensorZID[1] ) ? outer_z - typeA_sensor_half_length_incm : outer_z - typeB_sensor_half_length_incm, ( outer_zid == typeA_sensorZID[0] || outer_zid == typeA_sensorZID[1] ) ? outer_z + typeA_sensor_half_length_incm : outer_z + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1169 
1170     for (int possible_i = 0; possible_i < 2001; possible_i++) // todo : the segmentation of the strip
1171     {
1172         // double random_inner_z = rand_gen -> Uniform(inner_edge[0], inner_edge[1]);
1173         double random_inner_z = inner_edge[0] + ((inner_edge[1] - inner_edge[0]) / 2000.) * possible_i; 
1174         double edge_first  = Get_extrapolation(0, random_inner_z, inner_r, outer_edge[1], outer_r);
1175         double edge_second = Get_extrapolation(0, random_inner_z, inner_r, outer_edge[0], outer_r);
1176 
1177         double mid_point = (edge_first + edge_second) / 2.;
1178         double possible_width = fabs(edge_first - edge_second) / 2.;
1179 
1180         line_breakdown(combination_zvtx_range_shape, {mid_point - possible_width, mid_point + possible_width});
1181     }
1182 
1183     // note : this integral should take the overflow bins into account, since we normalize each of the trapezoidal shape,
1184     // note : even if some part of the shape is outside the range, it should be taken into account
1185     combination_zvtx_range_shape -> Scale(1./ combination_zvtx_range_shape -> Integral(-1,-1));
1186 
1187     for (int bin_i = 1; bin_i <= combination_zvtx_range_shape -> GetNbinsX(); bin_i++)
1188     {
1189         hist_in -> SetBinContent(bin_i, hist_in -> GetBinContent(bin_i) + combination_zvtx_range_shape -> GetBinContent(bin_i));
1190     }
1191 }
1192 
1193 std::vector<double> EvtVtxZProtoTracklet::find_Ngroup(TH1D * hist_in)
1194 {
1195     double Highest_bin_Content  = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
1196     double Highest_bin_Center   = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
1197 
1198     int group_Nbin = 0;
1199     int peak_group_ID = -9999;
1200     double group_entry = 0;
1201     double peak_group_ratio;
1202     std::vector<int> group_Nbin_vec; group_Nbin_vec.clear();
1203     std::vector<double> group_entry_vec; group_entry_vec.clear();
1204     std::vector<double> group_widthL_vec; group_widthL_vec.clear();
1205     std::vector<double> group_widthR_vec; group_widthR_vec.clear();
1206 
1207     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
1208         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
1209         double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
1210 
1211         if (bin_content != 0){
1212             
1213             if (group_Nbin == 0) {
1214                 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1215             }
1216 
1217             group_Nbin += 1; 
1218             group_entry += bin_content;
1219         }
1220         else if (bin_content == 0 && group_Nbin != 0){
1221             group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
1222             group_Nbin_vec.push_back(group_Nbin);
1223             group_entry_vec.push_back(group_entry);
1224             group_Nbin = 0;
1225             group_entry = 0;
1226         }
1227     }
1228     if (group_Nbin != 0) {
1229         group_Nbin_vec.push_back(group_Nbin);
1230         group_entry_vec.push_back(group_entry);
1231         group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
1232     } // note : the last group at the edge
1233 
1234     // note : find the peak group
1235     for (int i = 0; i < int(group_Nbin_vec.size()); i++){
1236         if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
1237             peak_group_ID = i;
1238             break;
1239         }
1240     }
1241     
1242     // note : On Nov 6, 2024, we no longer need to calculate the ratio of the peak group
1243     // double denominator = (accumulate( group_entry_vec.begin(), group_entry_vec.end(), 0.0 ));
1244     // denominator = (denominator <= 0) ? 1. : denominator;
1245     // peak_group_ratio = group_entry_vec[peak_group_ID] / denominator;
1246     peak_group_ratio = -999;
1247 
1248     double peak_group_left  = (double(group_Nbin_vec.size()) == 0) ? -999 : group_widthL_vec[peak_group_ID];
1249     double peak_group_right = (double(group_Nbin_vec.size()) == 0) ? 999  : group_widthR_vec[peak_group_ID];
1250 
1251     // for (int i = 0; i < group_Nbin_vec.size(); i++)
1252     // {
1253     //     cout<<" "<<endl;
1254     //     cout<<"group size : "<<group_Nbin_vec[i]<<endl;
1255     //     cout<<"group entry : "<<group_entry_vec[i]<<endl;
1256     //     cout<<group_widthL_vec[i]<<" "<<group_widthR_vec[i]<<endl;
1257     // }
1258 
1259     // cout<<" "<<endl;
1260     // cout<<"N group : "<<group_Nbin_vec.size()<<endl;
1261     // cout<<"Peak group ID : "<<peak_group_ID<<endl;
1262     // cout<<"peak group width : "<<group_widthL_vec[peak_group_ID]<<" "<<group_widthR_vec[peak_group_ID]<<endl;
1263     // cout<<"ratio : "<<peak_group_ratio<<endl;
1264     
1265     // note : {N_group, ratio (if two), peak widthL, peak widthR}
1266     return {double(group_Nbin_vec.size()), peak_group_ratio, peak_group_left, peak_group_right};
1267 }
1268 
1269 double  EvtVtxZProtoTracklet::vector_average (std::vector <double> input_vector) {
1270     return accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
1271 }
1272 
1273 double EvtVtxZProtoTracklet::vector_stddev (std::vector <double> input_vector){
1274     
1275     double sum_subt = 0;
1276     double average  = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
1277     
1278     // cout<<"average is : "<<average<<endl;
1279 
1280     for (int i=0; i<int(input_vector.size()); i++){ sum_subt += pow((input_vector[i] - average),2); }
1281 
1282     //cout<<"sum_subt : "<<sum_subt<<endl;
1283     // cout<<"print from the function, average : "<<average<<" std : "<<stddev<<endl;
1284 
1285     return sqrt( sum_subt / double(input_vector.size()-1) );
1286 }   
1287 
1288 
1289 double EvtVtxZProtoTracklet::get_clu_eta(std::vector<double> vertex, std::vector<double> clu_pos)
1290 {
1291     double correct_x = clu_pos[0] - vertex[0];
1292     double correct_y = clu_pos[1] - vertex[1];
1293     double correct_z = clu_pos[2] - vertex[2];
1294     double clu_r = sqrt(pow(correct_x,2) + pow(correct_y,2));
1295     // cout<<"correct info : "<<correct_x<<" "<<correct_y<<" "<<correct_z<<" "<<clu_r<<endl;
1296 
1297     return -0.5 * TMath::Log((sqrt(pow(correct_z,2)+pow(clu_r,2))-(correct_z)) / (sqrt(pow(correct_z,2)+pow(clu_r,2))+(correct_z)));
1298 }
1299 
1300 // note : pair<reduced chi2, eta of the track>
1301 // note : vector : {r,z}
1302 // note : p0 : vertex point {r,z,zerror}
1303 // note : p1 : inner layer
1304 // note : p2 : outer layer
1305 std::pair<double,double> EvtVtxZProtoTracklet::Get_eta(std::vector<double>p0, std::vector<double>p1, std::vector<double>p2)
1306 {
1307     // if (track_gr -> GetN() != 0) {cout<<"In EvtVtxZProtoTracklet class, track_gr is not empty, track_gr size : "<<track_gr -> GetN()<<endl; exit(0);}
1308     
1309     if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
1310     
1311     std::vector<double> r_vec  = {p0[0], p1[0], p2[0]}; 
1312     std::vector<double> z_vec  = {p0[1], p1[1], p2[1]}; 
1313     std::vector<double> re_vec = {0, 0, 0}; 
1314     std::vector<double> ze_vec = {
1315         p0[2],
1316         (p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1]) ? typeA_sensor_half_length_incm : typeB_sensor_half_length_incm,
1317         (p2[2] == typeA_sensorZID[0] || p2[2] == typeA_sensorZID[1]) ? typeA_sensor_half_length_incm : typeB_sensor_half_length_incm
1318     };
1319 
1320     // std::vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < sensor_type_segment ) ? sensor_length_typeA : sensor_length_typeB, ( fabs( p2[1] ) < sensor_type_segment ) ? sensor_length_typeA : sensor_length_typeB}; 
1321 
1322     // note : swap the r and z, to have easier fitting 
1323     // note : in principle, Z should be in the x axis, R should be in the Y axis
1324     for (int i = 0; i < 3; i++)
1325     {
1326         track_gr -> SetPoint(i,r_vec[i],z_vec[i]);
1327         track_gr -> SetPointError(i,re_vec[i],ze_vec[i]);
1328 
1329         // cout<<"In EvtVtxZProtoTracklet class, point : "<<i<<" r : "<<r_vec[i]<<" z : "<<z_vec[i]<<" re : "<<re_vec[i]<<" ze : "<<ze_vec[i]<<endl;
1330     }    
1331     
1332     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
1333     
1334     if (vertical_line) {return {0,0};}
1335     else 
1336     {
1337         fit_rz -> SetParameters(0,0);
1338 
1339         track_gr -> Fit(fit_rz,"NQ");
1340 
1341         std::pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
1342 
1343         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
1344     }
1345 
1346 }
1347 
1348 double EvtVtxZProtoTracklet::grEY_stddev(TGraphErrors * input_grr)
1349 {
1350     std::vector<double> input_vector; input_vector.clear();
1351     for (int i = 0; i < input_grr -> GetN(); i++) {input_vector.push_back( input_grr -> GetPointY(i) );}
1352 
1353     double average = std::accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
1354 
1355     double sum_subt = 0;
1356     for (int ele : input_vector) {sum_subt+=pow((ele-average),2);}
1357     
1358     return sqrt(sum_subt/(input_vector.size()-1));
1359 }   
1360 
1361 std::pair<double, double> EvtVtxZProtoTracklet::mirrorPolynomial(double a, double b) {
1362     // Interchange 'x' and 'y'
1363     double mirroredA = 1.0 / a;
1364     double mirroredB = -b / a;
1365 
1366     return {mirroredA, mirroredB};
1367 }
1368 
1369 void EvtVtxZProtoTracklet::Characterize_Pad (TPad *pad, float left, float right, float top, float bottom, bool set_logY, int setgrid_bool)
1370 {
1371     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
1372     pad -> SetLeftMargin   (left);
1373     pad -> SetRightMargin  (right);
1374     pad -> SetTopMargin    (top);
1375     pad -> SetBottomMargin (bottom);
1376     pad -> SetTicks(1,1);
1377     if (set_logY == true)
1378     {
1379         pad -> SetLogy (1);
1380     }
1381     
1382 }
1383 
1384 void EvtVtxZProtoTracklet::PrepareEvtCanvas()
1385 {
1386     std::string PlotOut_filename = "Plot_" + output_filename;
1387     INTTvtxZ_EvtDisplay_file_out = new TFile(Form("%s/%s",output_directory.c_str(),PlotOut_filename.c_str()),"RECREATE");
1388 
1389     c1 = new TCanvas("","",950,800);
1390     c1 -> cd();
1391 
1392     pad_EvtZDist = new TPad("pad_EvtZDist", "pad_EvtZDist", 0.0, 0.0, 1.0, 1.0);
1393     // Characterize_Pad(pad_EvtZDist, 0.15, 0.1, 0.1, 0.2, 0, 0);
1394     // pad_EvtZDist -> Draw();
1395 
1396     pad_ZoomIn_EvtZDist = new TPad("pad_ZoomIn_EvtZDist", "pad_ZoomIn_EvtZDist", 0.52, 0.15+0.1, 0.82, 0.5+0.1);
1397     Characterize_Pad(pad_ZoomIn_EvtZDist, 0.15, 0.1, 0.1, 0.2, 0, 0);
1398     pad_ZoomIn_EvtZDist -> SetFillColor(0); // Set background color of the pad to white
1399     // pad_ZoomIn_EvtZDist -> Draw();
1400 
1401     draw_text = new TLatex();
1402     draw_text -> SetName("draw_text");
1403     draw_text -> SetNDC();
1404     draw_text -> SetTextSize(0.03);
1405 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!