Back to home page

sPhenix code displayed by LXR

 
 

    


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

0001 #ifndef INTTReadTree_h
0002 #define INTTReadTree_h
0003 
0004 #include "../INTTDSTchain.C"
0005 #include "PrivateCluReader.C"
0006 #include "../private_cluster_gen/InttConversion_new.h"
0007 #include "../private_cluster_gen/InttClustering.h"
0008 #include "ReadMCselfgen/MCSelfGenReader.C"
0009 #include "ReadF4ANtupleMC/ReadF4ANtupleMC.C"
0010 #include "ReadF4ANtupleData/ReadF4ANtupleData.C"
0011 
0012 class INTTReadTree
0013 {
0014     public : 
0015         vector<clu_info> temp_sPH_inner_nocolumn_vec; 
0016         vector<clu_info> temp_sPH_outer_nocolumn_vec; 
0017         vector<vector<double>> temp_sPH_nocolumn_vec;
0018         vector<vector<double>> temp_sPH_nocolumn_rz_vec;
0019         
0020         INTTReadTree(
0021             int data_type, 
0022             string input_directory, 
0023             string MC_list_name, 
0024             string tree_name, 
0025             int clu_size_cut, 
0026             int clu_sum_adc_cut, 
0027             int random_seed = 65539, 
0028             int N_ladder = 14, 
0029             double offset_range = 0.2, 
0030             vector<string> included_ladder_vec = {}, 
0031             int apply_geo_offset = 1, 
0032             vector<vector<double>> ladder_offset_vec = {}
0033         );
0034         void EvtInit(long long event_i);
0035         void EvtSetCluGroup();
0036         long long GetNEvt();
0037         unsigned long GetEvtNClus();
0038         unsigned long GetEvtNClusPost();
0039         double GetTrigZvtxMC();
0040         vector<double> GetTrigvtxMC();
0041         vector<vector<float>> GetTrueTrackInfo();
0042         bool GetPhiCheckTag();
0043         float GetCentralityBin();
0044         int GetNvtxMC();
0045         string GetRunType();
0046         Long64_t GetBCOFull();
0047         double GetMBDRecoZ();
0048         int GetIsMinBias();
0049         int GetIsMinBiasWoZDC();
0050         double GetMBDNorthChargeSum();
0051         double GetMBDSouthChargeSum();
0052         void EvtClear();
0053         map<string, vector<double>> GetLadderOffsetMap();
0054         void EndRun();
0055         void gen_ladder_offset();
0056         void clear_ladder_offset_map() {ladder_offset_map.clear();}
0057         void set_random_seed(int random_seed_in) {random_seed = random_seed_in;}
0058         void set_ladder_offset(map<string, vector<double>> input_map) {ladder_offset_map = input_map;}
0059 
0060     private : 
0061         string data_type_list[9] = {
0062             "MC",
0063             "data_DST",
0064             "data_private", 
0065             "data_private_geo1", 
0066             "MC_geo_test", 
0067             "MC_inner_phi_rotation", 
0068             "MC_selfgen", 
0069             "data_F4A", 
0070             "data_F4A_inner_phi_rotation"
0071         };
0072         long long N_event;
0073 
0074         string input_directory; 
0075         string MC_list_name;
0076         string tree_name;
0077         double clu_sum_adc_cut;
0078         double offset_range;
0079         int apply_geo_offset;
0080         int random_seed;
0081         int clu_size_cut;
0082         int data_type;   
0083         int N_ladder;
0084         int NvtxMC; 
0085         int is_min_bias;  // note : for data and MC
0086         int is_min_bias_wozdc; // note : only for data
0087         double MBD_south_charge_sum, MBD_north_charge_sum;
0088         double MBD_reco_z;
0089         double TrigXvtxMC;
0090         double TrigYvtxMC;
0091         double TrigZvtxMC;
0092         float  Centrality_bimp;
0093         float  Centrality_mbd;
0094         Long64_t bco_full;
0095 
0096         vector<vector<float>>  true_track_info;
0097         vector<string> included_ladder_vec;
0098         vector<vector<double>> ladder_offset_vec;
0099         map<string, vector<double>> ladder_offset_map;
0100 
0101         TChain * chain_in;
0102         PrivateCluReader * inttCluData; // note : the class to read the private gen cluster file
0103         InttConversion * inttConv; // note : the class to conver the server_module into ladder name
0104         MCSelfGenReader * inttMCselfgen; // note : the class to read the self generated MC sample
0105         TFile * file_in = nullptr;
0106         TTree * tree;
0107 
0108         // note : the class to read the MC sample 
0109         // INTTDSTchain * inttDSTMC;    // note : read with TChain
0110         ReadF4ANtupleMC * inttDSTMC; // note : read with TTree, this is the latest
0111         ReadF4ANtupleData * inttDSTData; // note : read with TTree, this is the latest, for the data produced by F4A
0112 
0113         vector<int> MBin_NClus_cut_vec;
0114 
0115         string run_type_out;
0116 
0117         // vector<string> included_ladder_vec = {
0118         //     "B0L000S","B0L002S","B0L003S", "B0L005S", "B0L006S", "B0L008S", "B0L009S", 
0119         //     "B1L000S","B1L003S","B1L004S", "B1L007S", "B1L008S", "B1L011S", "B1L012S"
0120         // };// todo : should be updated in the future when the progress is made and we moved to the next geoemtry correction step 
0121         
0122         unsigned long evt_length;
0123 
0124         void TChainInit_MC();
0125         void TTreeInit_private();
0126         void TTreeInit_private_geo1();
0127         void TChainInit_MC_geo_test();
0128         void TTreeInit_MC_selfgen();
0129         void TTreeInit_data_F4A();
0130         double get_radius(double x, double y);
0131         pair<double,double> rotatePoint(double x, double y);
0132         vector<double> offset_correction(map<string,vector<double>> input_map);
0133         void read_centrality_cut();        // note : for the selfgen MC
0134         int get_centrality_bin(int NClus); // note : for the selfgen MC
0135 
0136 };
0137 
0138 INTTReadTree::INTTReadTree(int data_type, string input_directory, string MC_list_name, string tree_name, int clu_size_cut, int clu_sum_adc_cut, int random_seed, int N_ladder, double offset_range, vector<string> included_ladder_vec, int apply_geo_offset, vector<vector<double>> ladder_offset_vec)
0139 :data_type(data_type), input_directory(input_directory), MC_list_name(MC_list_name), tree_name(tree_name), clu_size_cut(clu_size_cut), clu_sum_adc_cut(clu_sum_adc_cut), random_seed(random_seed), N_ladder(N_ladder), offset_range(offset_range), included_ladder_vec(included_ladder_vec), apply_geo_offset(apply_geo_offset), ladder_offset_vec(ladder_offset_vec)
0140 {
0141     temp_sPH_inner_nocolumn_vec.clear();
0142     temp_sPH_outer_nocolumn_vec.clear();
0143     
0144     temp_sPH_nocolumn_vec.clear();
0145     temp_sPH_nocolumn_vec = vector<vector<double>> (2);
0146 
0147     temp_sPH_nocolumn_rz_vec.clear();
0148     temp_sPH_nocolumn_rz_vec = vector<vector<double>> (2);
0149 
0150     ladder_offset_map.clear();
0151 
0152     true_track_info.clear();
0153 
0154     MBin_NClus_cut_vec.clear();
0155 
0156     if (data_type_list[data_type] == "MC" || data_type_list[data_type] == "MC_inner_phi_rotation") 
0157     {
0158         run_type_out = "MC";
0159         std::cout<<"--- INTTReadTree -> input data type: MC ---"<<std::endl;
0160         TChainInit_MC();
0161         std::cout<<"--- INTTReadTree -> Initialization done ---"<<std::endl;
0162 
0163         if (data_type_list[data_type] == "MC_inner_phi_rotation") {std::cout<<"--- INTTReadTree -> note, the clusters in the inner barrel will be rotated ---"<<std::endl;}
0164     }
0165     else if (data_type_list[data_type] == "data_private")
0166     {
0167         run_type_out = "data_private";
0168         std::cout<<"--- INTTReadTree -> input data type: private gen cluster ---"<<std::endl;
0169         TTreeInit_private();
0170         std::cout<<"--- INTTReadTree -> Initialization done ---"<<std::endl;
0171     }
0172     else if (data_type_list[data_type] == "data_private_geo1")
0173     {
0174         run_type_out = "data_private_geo1";
0175         std::cout<<"--- INTTReadTree -> input data type: private gen cluster with geo correction 1 ---"<<std::endl;
0176         TTreeInit_private_geo1();
0177         std::cout<<"--- INTTReadTree -> Initialization done ---"<<std::endl;
0178     }
0179     else if (data_type_list[data_type] == "MC_geo_test")
0180     {
0181         run_type_out = "MC";
0182         // std::cout<<"--- INTTReadTree -> input data type: MC geometry test ---"<<std::endl; // note : was used
0183         TChainInit_MC_geo_test();
0184         // std::cout<<"--- INTTReadTree -> Initialization done ---"<<std::endl; // note : was used
0185     }
0186     else if (data_type_list[data_type] == "MC_selfgen")
0187     {
0188         run_type_out = "MC";
0189         read_centrality_cut();
0190         TTreeInit_MC_selfgen();
0191     }
0192     else if (data_type_list[data_type] == "data_F4A" || data_type_list[data_type] == "data_F4A_inner_phi_rotation")
0193     {
0194         run_type_out = "data";
0195         std::cout<<"--- INTTReadTree -> input data type: data F4A ---"<<std::endl;
0196         TTreeInit_data_F4A();
0197         std::cout<<"--- INTTReadTree -> Initialization done ---"<<std::endl;
0198     }
0199 
0200 }
0201 
0202 void INTTReadTree::read_centrality_cut()
0203 {
0204     TFile * file_MBin_cut = TFile::Open(Form("%s/MBin_NClus_cut.root", input_directory.c_str()));
0205     TTree * tree_MBin_cut = (TTree*)file_MBin_cut->Get("tree");
0206     vector<int> *MBinCut_vec; MBinCut_vec = 0;
0207     tree_MBin_cut -> SetBranchAddress("MBinCut", &MBinCut_vec);
0208     tree_MBin_cut -> GetEntry(0);
0209     cout<<"check MBinCut_vec size : "<<MBinCut_vec->size()<<endl;
0210     for (int i = 0; i < MBinCut_vec->size(); i++) {MBin_NClus_cut_vec.push_back(MBinCut_vec->at(i));}
0211 
0212     file_MBin_cut -> Close();
0213     return;
0214 }
0215 
0216 int INTTReadTree::get_centrality_bin(int NClus)
0217 {
0218     int MBin_id = -1;
0219     for (int i = 0; i < MBin_NClus_cut_vec.size(); i++)
0220     {
0221         if (NClus > MBin_NClus_cut_vec[i]) {MBin_id = i; break;}
0222     }
0223 
0224     if (MBin_id == -1) { MBin_id = MBin_NClus_cut_vec.size(); }
0225 
0226     return MBin_id;
0227 }
0228 
0229 void INTTReadTree::TTreeInit_MC_selfgen()
0230 {
0231     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0232     tree = (TTree *)file_in->Get(tree_name.c_str());   
0233     N_event = tree -> GetEntries();
0234     inttMCselfgen = new MCSelfGenReader(tree);
0235 }
0236 
0237 void INTTReadTree::TTreeInit_data_F4A()
0238 {
0239     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0240     tree = (TTree *)file_in->Get(tree_name.c_str());   
0241     N_event = tree -> GetEntries();
0242     inttDSTData = new ReadF4ANtupleData(tree);
0243 }
0244 
0245 void INTTReadTree::TChainInit_MC()
0246 {
0247     // note : this is the latest format, probably fixed
0248     // note : the idea is that, for the avgVTXxy, we are going to read the merged file
0249     // note : for the evtVTXz, and evtTracklet, we are going to read the separated files for each condor job
0250     // note : which means, we might not necessarily need to use TChain
0251     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0252     tree = (TTree *)file_in->Get(tree_name.c_str());   
0253     N_event = tree -> GetEntries();
0254     inttDSTMC = new ReadF4ANtupleMC(tree);
0255 
0256     // chain_in = new TChain(tree_name.c_str());
0257     // inttDSTMC = new INTTDSTchain(chain_in, input_directory, MC_list_name);
0258     // N_event = chain_in->GetEntries();
0259     // std::printf("inttDSTMC N event : %lli \n", N_event); // note : was used
0260 }
0261 
0262 void INTTReadTree::TChainInit_MC_geo_test()
0263 {
0264     // note : this is the latest format, probably fixed
0265     // note : the idea is that, for the avgVTXxy, we are going to read the merged file
0266     // note : for the evtVTXz, and evtTracklet, we are going to read the separated files for each condor job
0267     // note : which means, we might not necessarily need to use TChain
0268     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0269     tree = (TTree *)file_in->Get(tree_name.c_str());   
0270     N_event = tree -> GetEntries();
0271     inttDSTMC = new ReadF4ANtupleMC(tree);
0272 
0273     // chain_in = new TChain(tree_name.c_str());
0274     // inttDSTMC = new INTTDSTchain(chain_in, input_directory, MC_list_name);
0275     // N_event = chain_in->GetEntries();
0276     // std::printf("inttDSTMC N event : %lli \n", N_event); // note : was used
0277     // gen_ladder_offset();
0278 }
0279 
0280 void INTTReadTree::TTreeInit_private()
0281 {
0282     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0283     tree = (TTree *)file_in->Get(tree_name.c_str());
0284     inttCluData = new PrivateCluReader(tree);
0285     N_event = tree -> GetEntries();
0286     // std::printf("private gen tree, N event : %lli \n", N_event); // note : was used
0287 
0288 }
0289 
0290 void INTTReadTree::TTreeInit_private_geo1()
0291 {
0292     file_in = TFile::Open(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()));
0293     // file_in = new TFile(Form("%s/%s.root", input_directory.c_str(), MC_list_name.c_str()),"read");
0294     tree = (TTree *)file_in->Get(tree_name.c_str());
0295     inttCluData = new PrivateCluReader(tree);
0296     N_event = tree -> GetEntries();
0297     // std::printf("private gen tree, N event : %lli \n", N_event); // note : was used
0298     inttConv = new InttConversion();
0299 
0300     if (N_ladder != included_ladder_vec.size()) {cout<<"In INTTReadTree, there is an error in ladder offset generation"<<endl; exit(1);}    
0301     gen_ladder_offset();
0302 }
0303 
0304 void INTTReadTree::EvtInit(long long event_i)
0305 {
0306     if (data_type_list[data_type] == "MC" || data_type_list[data_type] == "MC_geo_test" || data_type_list[data_type] == "MC_inner_phi_rotation")
0307     {
0308         inttDSTMC->LoadTree(event_i);
0309         inttDSTMC->GetEntry(event_i);
0310 
0311         evt_length           = inttDSTMC->NClus;
0312         bco_full             = -1;
0313         is_min_bias          = inttDSTMC->is_min_bias;
0314         is_min_bias_wozdc    = inttDSTMC->is_min_bias_wozdc; 
0315         MBD_south_charge_sum = inttDSTMC->MBD_south_charge_sum;
0316         MBD_north_charge_sum = inttDSTMC->MBD_north_charge_sum;
0317         MBD_reco_z           = inttDSTMC->MBD_z_vtx * 10.; // note : from cm to mm
0318         NvtxMC               = inttDSTMC->NTruthVtx;
0319         TrigXvtxMC           = inttDSTMC->TruthPV_trig_x;
0320         TrigYvtxMC           = inttDSTMC->TruthPV_trig_y;
0321         TrigZvtxMC           = inttDSTMC->TruthPV_trig_z;
0322         Centrality_bimp      = inttDSTMC->MBD_centrality; // note : same thing for both
0323         Centrality_mbd       = inttDSTMC->MBD_centrality; // note : same thing for both
0324         for (int track_i = 0; track_i < inttDSTMC->PrimaryG4P_Eta->size(); track_i++) {true_track_info.push_back({inttDSTMC->PrimaryG4P_Eta->at(track_i),inttDSTMC->PrimaryG4P_Phi->at(track_i),float(inttDSTMC->PrimaryG4P_PID->at(track_i))});}
0325     }
0326     else if (data_type_list[data_type] == "data_private" || data_type_list[data_type] == "data_private_geo1")
0327     {
0328         inttCluData->LoadTree(event_i);
0329         inttCluData->GetEntry(event_i);
0330 
0331         evt_length           = inttCluData->x->size();
0332         bco_full             = inttCluData->bco_full;
0333         is_min_bias          = 1;
0334         is_min_bias_wozdc    = 1;
0335         MBD_south_charge_sum = -1000;
0336         MBD_north_charge_sum = -1000;
0337         MBD_reco_z           = -10000;
0338         NvtxMC               = -1;
0339         TrigXvtxMC           = -1000.;
0340         TrigYvtxMC           = -1000.;
0341         TrigZvtxMC           = -1000.;
0342     }
0343     else if (data_type_list[data_type] == "MC_selfgen")
0344     {
0345         inttMCselfgen->LoadTree(event_i);
0346         inttMCselfgen->GetEntry(event_i);
0347 
0348         evt_length           = inttMCselfgen->NClus;
0349         bco_full             = -1;
0350         NvtxMC               = 1;
0351         is_min_bias          = 1;
0352         is_min_bias_wozdc    = 1;
0353         MBD_south_charge_sum = -1000;
0354         MBD_north_charge_sum = -1000;
0355         MBD_reco_z           = -10000;
0356         TrigXvtxMC           = inttMCselfgen->TruthPV_x;
0357         TrigYvtxMC           = inttMCselfgen->TruthPV_y;
0358         TrigZvtxMC           = inttMCselfgen->TruthPV_z;
0359         Centrality_bimp = get_centrality_bin(inttMCselfgen->NClus);
0360         Centrality_mbd = get_centrality_bin(inttMCselfgen->NClus);
0361 
0362         // cout<<"inttMCselfgen->NClus : "<<inttMCselfgen->NClus<<" Centrality_bimp : "<<Centrality_bimp<<endl;
0363         
0364         for (int track_i = 0; track_i < inttMCselfgen->PrimaryG4P_Eta->size(); track_i++) {
0365             true_track_info.push_back({
0366                 inttMCselfgen->PrimaryG4P_Eta->at(track_i),
0367                 inttMCselfgen->PrimaryG4P_Phi->at(track_i),
0368                 float(inttMCselfgen->PrimaryG4P_PID->at(track_i))
0369             });
0370         }
0371     }
0372     else if (data_type_list[data_type] == "data_F4A" || data_type_list[data_type] == "data_F4A_inner_phi_rotation") 
0373     {
0374         inttDSTData->LoadTree(event_i);
0375         inttDSTData->GetEntry(event_i);
0376 
0377         evt_length           = inttDSTData->NClus;
0378         bco_full             = inttDSTData->INTT_BCO;
0379         NvtxMC               = 1;
0380         is_min_bias          = inttDSTData->is_min_bias;
0381         is_min_bias_wozdc    = inttDSTData->is_min_bias_wozdc;
0382         MBD_south_charge_sum = inttDSTData->MBD_south_charge_sum;
0383         MBD_north_charge_sum = inttDSTData->MBD_north_charge_sum;
0384         MBD_reco_z           = inttDSTData->MBD_z_vtx * 10.; // note : cm to mm
0385         TrigXvtxMC           = -1000;
0386         TrigYvtxMC           = -1000;
0387         TrigZvtxMC           = -1000;
0388         Centrality_bimp      = inttDSTData->MBD_centrality;
0389         Centrality_mbd       = inttDSTData->MBD_centrality;
0390     }
0391 }
0392 
0393 long long INTTReadTree::GetNEvt() { return N_event; }
0394 unsigned long INTTReadTree::GetEvtNClus() { return evt_length; }
0395 double INTTReadTree::GetTrigZvtxMC() {return TrigZvtxMC;}
0396 vector<double> INTTReadTree::GetTrigvtxMC() {return {TrigXvtxMC,TrigYvtxMC,TrigZvtxMC};}
0397 int INTTReadTree::GetNvtxMC() {return NvtxMC;}
0398 string INTTReadTree::GetRunType() { return run_type_out; }
0399 Long64_t INTTReadTree::GetBCOFull() {return bco_full;}
0400 map<string, vector<double>> INTTReadTree::GetLadderOffsetMap() {return ladder_offset_map;}
0401 float INTTReadTree::GetCentralityBin() {return Centrality_bimp;}
0402 vector<vector<float>> INTTReadTree::GetTrueTrackInfo() {return true_track_info;}
0403 double INTTReadTree::GetMBDRecoZ() {return MBD_reco_z;}
0404 int INTTReadTree::GetIsMinBias() {return is_min_bias;}
0405 int INTTReadTree::GetIsMinBiasWoZDC() {return is_min_bias_wozdc;}
0406 double INTTReadTree::GetMBDNorthChargeSum() {return MBD_north_charge_sum;}
0407 double INTTReadTree::GetMBDSouthChargeSum() {return MBD_south_charge_sum;}
0408 
0409 unsigned long INTTReadTree::GetEvtNClusPost() 
0410 { 
0411     return temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size(); 
0412 }
0413 
0414 vector<double> INTTReadTree::offset_correction(map<string,vector<double>> input_map)
0415 {
0416     double N_pair = 0;
0417     double sum_x = 0;
0418     double sum_y = 0;
0419     double sum_z = 0;
0420 
0421     for (const auto& pair : input_map)
0422     {
0423         N_pair += 1;
0424         sum_x += pair.second[0];
0425         sum_y += pair.second[1];
0426         sum_z += pair.second[2];
0427     }
0428 
0429     return {sum_x/N_pair, sum_y/N_pair, sum_z/N_pair};
0430 }
0431 
0432 void INTTReadTree::EvtSetCluGroup()
0433 {
0434     if (data_type_list[data_type] == "MC" || data_type_list[data_type] == "MC_inner_phi_rotation"){
0435         for (int clu_i = 0; clu_i < evt_length; clu_i++)
0436         {
0437             if (int(inttDSTMC -> ClusPhiSize -> at(clu_i)) > clu_size_cut) continue; 
0438             if (int(inttDSTMC -> ClusAdc -> at(clu_i)) <= clu_sum_adc_cut) continue;
0439 
0440             // double clu_x = inttDSTMC -> ClusX -> at(clu_i) * 10; // note : change the unit from cm to mm
0441             // double clu_y = inttDSTMC -> ClusY -> at(clu_i) * 10;
0442 
0443             int    clu_layer = (inttDSTMC -> ClusLayer -> at(clu_i) == 3 || inttDSTMC -> ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0444             // note : this is for rotating the clusters in the inner barrel by 180 degrees.
0445             double clu_x = (data_type_list[data_type] == "MC_inner_phi_rotation" && clu_layer == 0)? rotatePoint(inttDSTMC -> ClusX -> at(clu_i) * 10, inttDSTMC -> ClusY -> at(clu_i) * 10).first  : inttDSTMC -> ClusX -> at(clu_i) * 10; // note : change the unit from cm to mm
0446             double clu_y = (data_type_list[data_type] == "MC_inner_phi_rotation" && clu_layer == 0)? rotatePoint(inttDSTMC -> ClusX -> at(clu_i) * 10, inttDSTMC -> ClusY -> at(clu_i) * 10).second : inttDSTMC -> ClusY -> at(clu_i) * 10;
0447             double clu_z = inttDSTMC -> ClusZ -> at(clu_i) * 10;
0448             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0449             double clu_radius = get_radius(clu_x, clu_y);
0450 
0451             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0452             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0453             
0454             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0455             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0456             
0457 
0458             if (clu_layer == 0) {// note : inner
0459                 temp_sPH_inner_nocolumn_vec.push_back({
0460                     -1, 
0461                     -1, 
0462                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0463                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0464                     int(inttDSTMC -> ClusPhiSize -> at(clu_i)), 
0465                     clu_x, 
0466                     clu_y, 
0467                     clu_z, 
0468                     inttDSTMC -> ClusLayer -> at(clu_i), // note : should be 3 or 4, for the inner layer, this is for the mega cluster search
0469                     clu_phi,
0470                     int(inttDSTMC -> ClusLayer       -> at(clu_i)), // note : raw_layer_id
0471                     int(inttDSTMC -> ClusLadderPhiId -> at(clu_i)), // note : raw_ladder_id
0472                     int(inttDSTMC -> ClusLadderZId   -> at(clu_i)), // note : raw_Z_id
0473                     int(inttDSTMC -> ClusPhiSize     -> at(clu_i)), // note : raw_phi_size
0474                     int(inttDSTMC -> ClusZSize       -> at(clu_i))  // note : raw_z_size
0475                 });
0476             }
0477             
0478             if (clu_layer == 1) {// note : outer
0479                 temp_sPH_outer_nocolumn_vec.push_back({
0480                     -1, 
0481                     -1, 
0482                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0483                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0484                     int(inttDSTMC -> ClusPhiSize -> at(clu_i)), 
0485                     clu_x, 
0486                     clu_y, 
0487                     clu_z, 
0488                     inttDSTMC -> ClusLayer -> at(clu_i), // note : should be 5 or 6, for the outer layer, this is for the mega cluster search
0489                     clu_phi,
0490                     int(inttDSTMC -> ClusLayer       -> at(clu_i)), // note : raw_layer_id
0491                     int(inttDSTMC -> ClusLadderPhiId -> at(clu_i)), // note : raw_ladder_id
0492                     int(inttDSTMC -> ClusLadderZId   -> at(clu_i)), // note : raw_Z_id
0493                     int(inttDSTMC -> ClusPhiSize     -> at(clu_i)), // note : raw_phi_size
0494                     int(inttDSTMC -> ClusZSize       -> at(clu_i))  // note : raw_z_size
0495                 });            
0496             }        
0497         }
0498     }
0499 
0500     else if (data_type_list[data_type] == "data_private"){
0501         for (int clu_i = 0; clu_i < evt_length; clu_i++)
0502         {
0503             if (int(inttCluData -> size -> at(clu_i)) > clu_size_cut) continue; 
0504             if (int(inttCluData -> sum_adc_conv -> at(clu_i)) <= clu_sum_adc_cut) continue;
0505 
0506             double clu_x = inttCluData -> x -> at(clu_i);
0507             double clu_y = inttCluData -> y -> at(clu_i);
0508             double clu_z = inttCluData -> z -> at(clu_i);
0509             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0510             int    clu_layer = inttCluData -> layer -> at(clu_i); // note : should be 0 or 1
0511             double clu_radius = get_radius(clu_x, clu_y);
0512 
0513             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0514             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0515             
0516             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0517             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0518             
0519 
0520             if (clu_layer == 0) {// note : inner
0521                 temp_sPH_inner_nocolumn_vec.push_back({
0522                     -1, 
0523                     -1, 
0524                     int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0525                     int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0526                     int(inttCluData -> size -> at(clu_i)), 
0527                     clu_x, 
0528                     clu_y, 
0529                     clu_z, 
0530                     clu_layer, 
0531                     clu_phi
0532                 });
0533             }
0534             
0535             if (clu_layer == 1) {// note : outer
0536                 temp_sPH_outer_nocolumn_vec.push_back({
0537                     -1, 
0538                     -1, 
0539                     int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0540                     int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0541                     int(inttCluData -> size -> at(clu_i)), 
0542                     clu_x, 
0543                     clu_y, 
0544                     clu_z, 
0545                     clu_layer, 
0546                     clu_phi
0547                 });            
0548             }        
0549         }
0550     }
0551 
0552     // else if (data_type_list[data_type] == "data_private_geo1"){
0553     //     for (int clu_i = 0; clu_i < evt_length; clu_i++)
0554     //     {
0555     //         if (int(inttCluData -> size -> at(clu_i)) > clu_size_cut) {continue;} 
0556     //         if (int(inttCluData -> sum_adc_conv -> at(clu_i)) <= clu_sum_adc_cut) {continue;}
0557 
0558     //         string ladder_name = inttConv -> GetLadderName(Form("intt%i_%i", inttCluData -> server -> at(clu_i) - 3001, inttCluData -> module -> at(clu_i)));
0559 
0560     //         // note : the following selections only keeps a small portion of INTT clusters
0561     //         // if (inttCluData -> server -> at(clu_i) - 3001 > 3) {continue;} // note : only the 3001, 3002, 3003, 3004 are included in the study (south)
0562     //         if ((ladder_name).substr(0,4) == "B0L1") {continue;}
0563     //         if ((ladder_name).substr(0,4) == "B1L1") {continue;}
0564     //         // if (inttCluData -> column -> at(clu_i) > 5) {continue;} // note : only the column 1, 2, 3, 4, 5 are included in the study (south type B sensor only)
0565 
0566     //         if ( (ladder_name).substr(0,6) == "B0L001") {continue;}
0567     //         if ( (ladder_name).substr(0,6) == "B0L004") {continue;}
0568     //         if ( (ladder_name).substr(0,6) == "B0L007") {continue;}
0569     //         if ( (ladder_name).substr(0,6) == "B0L010") {continue;}
0570 
0571     //         if ( (ladder_name).substr(0,6) == "B1L001") {continue;}
0572     //         if ( (ladder_name).substr(0,6) == "B1L002") {continue;}
0573     //         if ( (ladder_name).substr(0,6) == "B1L005") {continue;}
0574     //         if ( (ladder_name).substr(0,6) == "B1L006") {continue;}
0575     //         if ( (ladder_name).substr(0,6) == "B1L009") {continue;}
0576     //         if ( (ladder_name).substr(0,6) == "B1L010") {continue;}
0577     //         if ( (ladder_name).substr(0,6) == "B1L013") {continue;}
0578     //         if ( (ladder_name).substr(0,6) == "B1L014") {continue;}
0579 
0580     //         // note : because the corresponding hald-ladder of this one, the B1L015S, has no data
0581     //         if ( (ladder_name).substr(0,6) == "B0L011") {continue;}
0582     //         if ( (ladder_name).substr(0,6) == "B1L015") {continue;}
0583 
0584     //         // note : for the hypothesis test
0585     //         // if ( (ladder_name).substr(0,6) == "B1L004") {continue;}
0586     //         // if ( (ladder_name).substr(0,6) == "B0L003") {continue;} // note : pair
0587     //         // if ( (ladder_name).substr(0,6) == "B1L007") {continue;}
0588     //         // if ( (ladder_name).substr(0,6) == "B0L005") {continue;} // note : pair
0589     //         // if ( (ladder_name).substr(0,6) == "B1L011") {continue;}
0590     //         // if ( (ladder_name).substr(0,6) == "B0L008") {continue;} // note : pair
0591     //         // if ( (ladder_name).substr(0,6) == "B1L012") {continue;}
0592     //         // if ( (ladder_name).substr(0,6) == "B0L009") {continue;} // note : pair
0593     //         // if ( (ladder_name).substr(0,6) == "B0L000") {continue;}
0594     //         // if ( (ladder_name).substr(0,6) == "B1L000") {continue;} // note : pair
0595 
0596     //         double clu_x_offset = ladder_offset_map[(ladder_name).substr(0,6)].first;
0597     //         double clu_y_offset = ladder_offset_map[(ladder_name).substr(0,6)].second;
0598 
0599     //         double clu_x = inttCluData -> x -> at(clu_i) + clu_x_offset;
0600     //         double clu_y = inttCluData -> y -> at(clu_i) + clu_y_offset;
0601     //         double clu_z = inttCluData -> z -> at(clu_i);
0602     //         double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0603     //         int    clu_layer = inttCluData -> layer -> at(clu_i); // note : should be 0 or 1
0604     //         double clu_radius = get_radius(clu_x, clu_y);
0605 
0606     //         temp_sPH_nocolumn_vec[0].push_back( clu_x );
0607     //         temp_sPH_nocolumn_vec[1].push_back( clu_y );
0608             
0609     //         temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0610     //         temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0611             
0612 
0613     //         if (clu_layer == 0) {// note : inner
0614     //             temp_sPH_inner_nocolumn_vec.push_back({
0615     //                 -1, 
0616     //                 -1, 
0617     //                 int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0618     //                 int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0619     //                 int(inttCluData -> size -> at(clu_i)), 
0620     //                 clu_x, 
0621     //                 clu_y, 
0622     //                 clu_z, 
0623     //                 clu_layer, 
0624     //                 clu_phi
0625     //             });
0626     //         }
0627             
0628     //         if (clu_layer == 1) {// note : outer
0629     //             temp_sPH_outer_nocolumn_vec.push_back({
0630     //                 -1, 
0631     //                 -1, 
0632     //                 int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0633     //                 int(inttCluData -> sum_adc_conv -> at(clu_i)), 
0634     //                 int(inttCluData -> size -> at(clu_i)), 
0635     //                 clu_x, 
0636     //                 clu_y, 
0637     //                 clu_z, 
0638     //                 clu_layer, 
0639     //                 clu_phi
0640     //             });            
0641     //         }        
0642     //     }
0643     // }
0644 
0645     else if (data_type_list[data_type] == "MC_geo_test"){
0646         for (int clu_i = 0; clu_i < evt_length; clu_i++)
0647         {
0648             // if (inttDSTMC -> ClusLadderZId -> at(clu_i) != 0) continue; // note : only the south-side sensor B is included in the study
0649             if (int(inttDSTMC -> ClusPhiSize -> at(clu_i)) > clu_size_cut) continue; 
0650             if (int(inttDSTMC -> ClusAdc -> at(clu_i)) <= clu_sum_adc_cut) continue;
0651 
0652             double clu_x_offset = ladder_offset_map[Form("%i_%i",inttDSTMC -> ClusLayer -> at(clu_i), inttDSTMC -> ClusLadderPhiId -> at(clu_i))][0];
0653             double clu_y_offset = ladder_offset_map[Form("%i_%i",inttDSTMC -> ClusLayer -> at(clu_i), inttDSTMC -> ClusLadderPhiId -> at(clu_i))][1];
0654             double clu_z_offset = ladder_offset_map[Form("%i_%i",inttDSTMC -> ClusLayer -> at(clu_i), inttDSTMC -> ClusLadderPhiId -> at(clu_i))][2];
0655 
0656             double clu_x = inttDSTMC -> ClusX -> at(clu_i) * 10 + clu_x_offset; // note : change the unit from cm to mm
0657             double clu_y = inttDSTMC -> ClusY -> at(clu_i) * 10 + clu_y_offset;
0658             double clu_z = inttDSTMC -> ClusZ -> at(clu_i) * 10 + clu_z_offset;
0659             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0660             int    clu_layer = (inttDSTMC -> ClusLayer -> at(clu_i) == 3 || inttDSTMC -> ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0661             double clu_radius = get_radius(clu_x, clu_y);
0662 
0663             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0664             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0665             
0666             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0667             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0668             
0669 
0670             if (clu_layer == 0) {// note : inner
0671                 temp_sPH_inner_nocolumn_vec.push_back({
0672                     -1, 
0673                     -1, 
0674                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0675                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0676                     int(inttDSTMC -> ClusPhiSize -> at(clu_i)), 
0677                     clu_x, 
0678                     clu_y, 
0679                     clu_z, 
0680                     clu_layer, 
0681                     clu_phi
0682                 });
0683             }
0684             
0685             if (clu_layer == 1) {// note : outer
0686                 temp_sPH_outer_nocolumn_vec.push_back({
0687                     -1, 
0688                     -1, 
0689                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0690                     int(inttDSTMC -> ClusAdc -> at(clu_i)), 
0691                     int(inttDSTMC -> ClusPhiSize -> at(clu_i)), 
0692                     clu_x, 
0693                     clu_y, 
0694                     clu_z, 
0695                     clu_layer, 
0696                     clu_phi
0697                 });            
0698             }        
0699         }
0700     }
0701 
0702     else if (data_type_list[data_type] == "MC_selfgen"){
0703         for (int clu_i = 0; clu_i < evt_length; clu_i++)
0704         {
0705             if (int(inttMCselfgen -> phi_size -> at(clu_i)) > clu_size_cut) continue; 
0706             if (int(inttMCselfgen -> adc -> at(clu_i)) <= clu_sum_adc_cut) continue;
0707 
0708             // double clu_x = inttMCselfgen -> X -> at(clu_i) * 10; // note : change the unit from cm to mm
0709             // double clu_y = inttMCselfgen -> Y -> at(clu_i) * 10;
0710 
0711             int    clu_layer = (inttMCselfgen -> layer -> at(clu_i) == 3 || inttMCselfgen -> layer -> at(clu_i) == 4) ? 0 : 1;
0712             // note : this is for rotating the clusters in the inner barrel by 180 degrees.
0713             double clu_x = inttMCselfgen -> X -> at(clu_i) * 10; // note : change the unit from cm to mm
0714             double clu_y = inttMCselfgen -> Y -> at(clu_i) * 10;
0715             double clu_z = inttMCselfgen -> Z -> at(clu_i) * 10;
0716             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0717             double clu_radius = get_radius(clu_x, clu_y);
0718 
0719             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0720             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0721             
0722             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0723             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0724             
0725 
0726             if (clu_layer == 0) {// note : inner
0727                 temp_sPH_inner_nocolumn_vec.push_back({
0728                     -1, 
0729                     -1, 
0730                     int(inttMCselfgen -> adc -> at(clu_i)), 
0731                     int(inttMCselfgen -> adc -> at(clu_i)), 
0732                     int(inttMCselfgen -> phi_size -> at(clu_i)), 
0733                     clu_x, 
0734                     clu_y, 
0735                     clu_z, 
0736                     inttMCselfgen -> layer -> at(clu_i), // note : should be 3 or 4, for the inner layer, this is for the mega cluster search
0737                     clu_phi
0738                 });
0739             }
0740             
0741             if (clu_layer == 1) {// note : outer
0742                 temp_sPH_outer_nocolumn_vec.push_back({
0743                     -1, 
0744                     -1, 
0745                     int(inttMCselfgen -> adc -> at(clu_i)), 
0746                     int(inttMCselfgen -> adc -> at(clu_i)), 
0747                     int(inttMCselfgen -> phi_size -> at(clu_i)), 
0748                     clu_x, 
0749                     clu_y, 
0750                     clu_z, 
0751                     inttMCselfgen -> layer -> at(clu_i), // note : should be 5 or 6, for the outer layer, this is for the mega cluster search
0752                     clu_phi
0753                 });            
0754             }        
0755         }
0756     }
0757 
0758     else if (data_type_list[data_type] == "data_F4A" || data_type_list[data_type] == "data_F4A_inner_phi_rotation"){
0759         for (int clu_i = 0; clu_i < evt_length; clu_i++)
0760         {
0761             if (int(inttDSTData -> ClusPhiSize -> at(clu_i)) > clu_size_cut) continue; 
0762             if (int(inttDSTData -> ClusAdc -> at(clu_i)) <= clu_sum_adc_cut) continue;
0763 
0764             // double clu_x = inttDSTData -> ClusX -> at(clu_i) * 10; // note : change the unit from cm to mm
0765             // double clu_y = inttDSTData -> ClusY -> at(clu_i) * 10;
0766 
0767             int    clu_layer = (inttDSTData -> ClusLayer -> at(clu_i) == 3 || inttDSTData -> ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0768             // note : this is for rotating the clusters in the inner barrel by 180 degrees.
0769             double clu_x = (data_type_list[data_type] == "data_F4A_inner_phi_rotation" && clu_layer == 0)? rotatePoint(inttDSTData -> ClusX -> at(clu_i) * 10, inttDSTData -> ClusY -> at(clu_i) * 10).first  : inttDSTData -> ClusX -> at(clu_i) * 10; // note : change the unit from cm to mm
0770             double clu_y = (data_type_list[data_type] == "data_F4A_inner_phi_rotation" && clu_layer == 0)? rotatePoint(inttDSTData -> ClusX -> at(clu_i) * 10, inttDSTData -> ClusY -> at(clu_i) * 10).second : inttDSTData -> ClusY -> at(clu_i) * 10;
0771             double clu_z = inttDSTData -> ClusZ -> at(clu_i) * 10;
0772             double clu_phi = (clu_y < 0) ? atan2(clu_y,clu_x) * (180./TMath::Pi()) + 360 : atan2(clu_y,clu_x) * (180./TMath::Pi());
0773             double clu_radius = get_radius(clu_x, clu_y);
0774 
0775             temp_sPH_nocolumn_vec[0].push_back( clu_x );
0776             temp_sPH_nocolumn_vec[1].push_back( clu_y );
0777             
0778             temp_sPH_nocolumn_rz_vec[0].push_back( clu_z );
0779             temp_sPH_nocolumn_rz_vec[1].push_back( ( clu_phi > 180 ) ? clu_radius * -1 : clu_radius );
0780             
0781 
0782             if (clu_layer == 0) {// note : inner
0783                 temp_sPH_inner_nocolumn_vec.push_back({
0784                     -1, 
0785                     -1, 
0786                     int(inttDSTData -> ClusAdc -> at(clu_i)), 
0787                     int(inttDSTData -> ClusAdc -> at(clu_i)), 
0788                     int(inttDSTData -> ClusPhiSize -> at(clu_i)), 
0789                     clu_x, 
0790                     clu_y, 
0791                     clu_z, 
0792                     inttDSTData -> ClusLayer -> at(clu_i), // note : should be 3 or 4, for the inner layer, this is for the mega cluster search
0793                     clu_phi,
0794                     int(inttDSTData -> ClusLayer       -> at(clu_i)), // note : raw_layer_id
0795                     int(inttDSTData -> ClusLadderPhiId -> at(clu_i)), // note : raw_ladder_id
0796                     int(inttDSTData -> ClusLadderZId   -> at(clu_i)), // note : raw_Z_id
0797                     int(inttDSTData -> ClusPhiSize     -> at(clu_i)), // note : raw_phi_size
0798                     int(inttDSTData -> ClusZSize       -> at(clu_i))  // note : raw_z_size
0799                 });
0800             }
0801             
0802             if (clu_layer == 1) {// note : outer
0803                 temp_sPH_outer_nocolumn_vec.push_back({
0804                     -1, 
0805                     -1, 
0806                     int(inttDSTData -> ClusAdc -> at(clu_i)), 
0807                     int(inttDSTData -> ClusAdc -> at(clu_i)), 
0808                     int(inttDSTData -> ClusPhiSize -> at(clu_i)), 
0809                     clu_x, 
0810                     clu_y, 
0811                     clu_z, 
0812                     inttDSTData -> ClusLayer -> at(clu_i), // note : should be 5 or 6, for the outer layer, this is for the mega cluster search
0813                     clu_phi,
0814                     int(inttDSTData -> ClusLayer       -> at(clu_i)), // note : raw_layer_id
0815                     int(inttDSTData -> ClusLadderPhiId -> at(clu_i)), // note : raw_ladder_id
0816                     int(inttDSTData -> ClusLadderZId   -> at(clu_i)), // note : raw_Z_id
0817                     int(inttDSTData -> ClusPhiSize     -> at(clu_i)), // note : raw_phi_size
0818                     int(inttDSTData -> ClusZSize       -> at(clu_i))  // note : raw_z_size
0819                 });            
0820             }        
0821         }
0822     }
0823     
0824 }
0825 
0826 void INTTReadTree::gen_ladder_offset()
0827 {
0828     TRandom * rand = new TRandom3();
0829     rand -> SetSeed(random_seed);
0830 
0831     if (apply_geo_offset == 1)
0832     {
0833         for (int ladder_i = 0; ladder_i < N_ladder; ladder_i++)
0834         {
0835             ladder_offset_map[included_ladder_vec[ladder_i].c_str()] = {rand -> Uniform(-offset_range, offset_range), rand -> Uniform(-offset_range, offset_range)};    
0836         }
0837     }
0838     else if (apply_geo_offset == 0) // note : no geometry offset applied
0839     {
0840         for (int ladder_i = 0; ladder_i < N_ladder; ladder_i++)
0841         {
0842             ladder_offset_map[included_ladder_vec[ladder_i].c_str()] = {0, 0};    
0843         }
0844     }
0845     else if (apply_geo_offset == 2) // note : the geometry offset is given from outside
0846     {
0847         for (int ladder_i = 0; ladder_i < N_ladder; ladder_i++)
0848         {
0849             ladder_offset_map[included_ladder_vec[ladder_i].c_str()] = ladder_offset_vec[ladder_i];    
0850         }
0851     }
0852     else if (apply_geo_offset == 3) // note : the geometry offset for the MC
0853     {
0854         for (int layer_i = 3; layer_i < 7; layer_i++)
0855         {
0856             double N_layer_ladder = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0857 
0858             for (int phi_i = 0; phi_i < N_layer_ladder; phi_i++)
0859             {
0860                 ladder_offset_map[Form("%i_%i", layer_i, phi_i)] = {rand -> Uniform(-offset_range, offset_range), rand -> Uniform(-offset_range, offset_range), rand -> Uniform(-offset_range, offset_range)};
0861                 // ladder_offset_map[Form("%i_%i", layer_i, phi_i)] = {100., 100.};
0862             }
0863         }
0864 
0865         // note : it's possible that the whole system has a systematic offset, so we need to correct it
0866         vector<double> XYZ_correction = offset_correction(ladder_offset_map);
0867         for (const auto& pair : ladder_offset_map)
0868         {
0869             ladder_offset_map[pair.first] = {pair.second[0] - XYZ_correction[0], pair.second[1] - XYZ_correction[1], pair.second[2] - XYZ_correction[2]};
0870         }
0871         
0872     }
0873     
0874 
0875 
0876     rand -> Delete();
0877 }
0878 
0879 double INTTReadTree::get_radius(double x, double y)
0880 {
0881     return sqrt(pow(x,2)+pow(y,2));
0882 }
0883 
0884 bool INTTReadTree::GetPhiCheckTag()
0885 {
0886     int inner_1_check = 0;
0887     int inner_2_check = 0;
0888     int inner_3_check = 0;
0889     int inner_4_check = 0;
0890     for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
0891         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 0 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 90)    inner_1_check = 1;
0892         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 90 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 180)  inner_2_check = 1;
0893         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 180 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 270) inner_3_check = 1;
0894         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 270 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 360) inner_4_check = 1;
0895 
0896         if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check) == 4 ) break;
0897     }
0898 
0899     int outer_1_check = 0;
0900     int outer_2_check = 0;
0901     int outer_3_check = 0;
0902     int outer_4_check = 0;
0903     for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
0904         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 0 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 90)    outer_1_check = 1;
0905         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 90 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 180)  outer_2_check = 1;
0906         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 180 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 270) outer_3_check = 1;
0907         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 270 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 360) outer_4_check = 1;
0908 
0909         if ( (outer_1_check + outer_2_check + outer_3_check + outer_4_check) == 4 ) break;
0910     }
0911 
0912     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 ) {return false;}
0913     else { return true; }
0914 }
0915 
0916 void INTTReadTree::EvtClear()
0917 {
0918     temp_sPH_inner_nocolumn_vec.clear();
0919     temp_sPH_outer_nocolumn_vec.clear();
0920 
0921     temp_sPH_nocolumn_vec.clear();
0922     temp_sPH_nocolumn_vec = vector<vector<double>> (2);
0923     
0924     temp_sPH_nocolumn_rz_vec.clear();
0925     temp_sPH_nocolumn_rz_vec = vector<vector<double>> (2);   
0926 
0927     true_track_info.clear();
0928 
0929     evt_length = 0;
0930     NvtxMC     = 0;
0931 }
0932 
0933 void INTTReadTree::EndRun()
0934 {
0935     ladder_offset_map.clear();
0936     if (file_in != nullptr){
0937         cout<<"closing the file_in"<<endl;
0938         file_in -> Close();
0939         delete file_in;
0940     }
0941 
0942 
0943 }
0944 
0945 // note : rotate the ClusX and ClusY
0946 pair<double,double> INTTReadTree::rotatePoint(double x, double y)
0947 {
0948     // Convert the rotation angle from degrees to radians
0949     double rotation = 180.; // todo rotation is here
0950     double angleRad = rotation * M_PI / 180.0;
0951 
0952     // Perform the rotation
0953     double xOut = x * cos(angleRad) - y * sin(angleRad);
0954     double yOut = x * sin(angleRad) + y * cos(angleRad);
0955 
0956     xOut = (fabs(xOut) < 0.00000001) ? 0 : xOut;
0957     yOut = (fabs(yOut) < 0.00000001) ? 0 : yOut;
0958 
0959     // cout<<"Post rotation: "<<xOut<<" "<<yOut<<endl;
0960 
0961     return {xOut,yOut};
0962 }
0963 
0964 #endif