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;
0086 int is_min_bias_wozdc;
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;
0103 InttConversion * inttConv;
0104 MCSelfGenReader * inttMCselfgen;
0105 TFile * file_in = nullptr;
0106 TTree * tree;
0107
0108
0109
0110 ReadF4ANtupleMC * inttDSTMC;
0111 ReadF4ANtupleData * inttDSTData;
0112
0113 vector<int> MBin_NClus_cut_vec;
0114
0115 string run_type_out;
0116
0117
0118
0119
0120
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();
0134 int get_centrality_bin(int NClus);
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
0183 TChainInit_MC_geo_test();
0184
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
0248
0249
0250
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
0257
0258
0259
0260 }
0261
0262 void INTTReadTree::TChainInit_MC_geo_test()
0263 {
0264
0265
0266
0267
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
0274
0275
0276
0277
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
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
0294 tree = (TTree *)file_in->Get(tree_name.c_str());
0295 inttCluData = new PrivateCluReader(tree);
0296 N_event = tree -> GetEntries();
0297
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.;
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;
0323 Centrality_mbd = inttDSTMC->MBD_centrality;
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
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.;
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
0441
0442
0443 int clu_layer = (inttDSTMC -> ClusLayer -> at(clu_i) == 3 || inttDSTMC -> ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0444
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;
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) {
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),
0469 clu_phi,
0470 int(inttDSTMC -> ClusLayer -> at(clu_i)),
0471 int(inttDSTMC -> ClusLadderPhiId -> at(clu_i)),
0472 int(inttDSTMC -> ClusLadderZId -> at(clu_i)),
0473 int(inttDSTMC -> ClusPhiSize -> at(clu_i)),
0474 int(inttDSTMC -> ClusZSize -> at(clu_i))
0475 });
0476 }
0477
0478 if (clu_layer == 1) {
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),
0489 clu_phi,
0490 int(inttDSTMC -> ClusLayer -> at(clu_i)),
0491 int(inttDSTMC -> ClusLadderPhiId -> at(clu_i)),
0492 int(inttDSTMC -> ClusLadderZId -> at(clu_i)),
0493 int(inttDSTMC -> ClusPhiSize -> at(clu_i)),
0494 int(inttDSTMC -> ClusZSize -> at(clu_i))
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);
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) {
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) {
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
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
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
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;
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) {
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) {
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
0709
0710
0711 int clu_layer = (inttMCselfgen -> layer -> at(clu_i) == 3 || inttMCselfgen -> layer -> at(clu_i) == 4) ? 0 : 1;
0712
0713 double clu_x = inttMCselfgen -> X -> at(clu_i) * 10;
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) {
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),
0737 clu_phi
0738 });
0739 }
0740
0741 if (clu_layer == 1) {
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),
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
0765
0766
0767 int clu_layer = (inttDSTData -> ClusLayer -> at(clu_i) == 3 || inttDSTData -> ClusLayer -> at(clu_i) == 4) ? 0 : 1;
0768
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;
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) {
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),
0793 clu_phi,
0794 int(inttDSTData -> ClusLayer -> at(clu_i)),
0795 int(inttDSTData -> ClusLadderPhiId -> at(clu_i)),
0796 int(inttDSTData -> ClusLadderZId -> at(clu_i)),
0797 int(inttDSTData -> ClusPhiSize -> at(clu_i)),
0798 int(inttDSTData -> ClusZSize -> at(clu_i))
0799 });
0800 }
0801
0802 if (clu_layer == 1) {
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),
0813 clu_phi,
0814 int(inttDSTData -> ClusLayer -> at(clu_i)),
0815 int(inttDSTData -> ClusLadderPhiId -> at(clu_i)),
0816 int(inttDSTData -> ClusLadderZId -> at(clu_i)),
0817 int(inttDSTData -> ClusPhiSize -> at(clu_i)),
0818 int(inttDSTData -> ClusZSize -> at(clu_i))
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)
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)
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)
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
0862 }
0863 }
0864
0865
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
0946 pair<double,double> INTTReadTree::rotatePoint(double x, double y)
0947 {
0948
0949 double rotation = 180.;
0950 double angleRad = rotation * M_PI / 180.0;
0951
0952
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
0960
0961 return {xOut,yOut};
0962 }
0963
0964 #endif