NewCode2024/GeoOffset/Run_GeoOffset.C

0001 #include "../AvgVtxXY/AvgVtxXY.h"
0002 #include "../EvtVtxZTracklet/EvtVtxZProtoTracklet.h"
0003 #include "../TrackletHistogram/TrackletHistogramNew.h"
0004 #include "../PreparedNdEta/PreparedNdEtaEach.h"
0005
0006 R__LOAD_LIBRARY(../AvgVtxXY/libavgvtxxy.so)
0007 R__LOAD_LIBRARY(../EvtVtxZTracklet/libEvtVtxZProtoTracklet.so)
0008 R__LOAD_LIBRARY(../TrackletHistogram/libTrackletHistogramNew.so)
0009 R__LOAD_LIBRARY(../PreparedNdEta/libPreparedNdEtaEach.so)
0010
0011 TH2D * GetGoodColMap (std::string ColMulMask_map_dir_in, std::string ColMulMask_map_file_in, std::string map_name_in)
0012 {
0013   TFile * f = TFile::Open(Form("%s/%s", ColMulMask_map_dir_in.c_str(), ColMulMask_map_file_in.c_str()));
0014   TH2D * h = (TH2D*)f->Get(map_name_in.c_str());
0015   return h;
0016 }
0017
0018 TH1D * GetDataVtxZH1D (std::string data_vtxZ_directory_file_in, std::string data_vtxZ_hist_name_in)
0019 {
0020   TFile * f2 = TFile::Open(Form("%s", data_vtxZ_directory_file_in.c_str()));
0021   TH1D * h2 = (TH1D*)f2->Get(data_vtxZ_hist_name_in.c_str());
0022   return h2;
0023 }
0024
0025 TH1D * GetVtxZReweightH1D (TH1D * data_hist, TH1D * mc_hist)
0026 {
0027     if (data_hist->GetNbinsX() != mc_hist->GetNbinsX()){
0028         cout<<"GetVtxZReweightH1D: data_hist->GetNbinsX() != mc_hist->GetNbinsX()"<<endl;
0029         exit(1);
0030         return nullptr;
0031     }
0032     if (data_hist->GetXaxis()->GetXmin() != mc_hist->GetXaxis()->GetXmin()){
0033         cout<<"GetVtxZReweightH1D: data_hist->GetXaxis()->GetXmin() != mc_hist->GetXaxis()->GetXmin()"<<endl;
0034         exit(1);
0035         return nullptr;
0036     }
0037     if (data_hist->GetXaxis()->GetXmax() != mc_hist->GetXaxis()->GetXmax()){
0038         cout<<"GetVtxZReweightH1D: data_hist->GetXaxis()->GetXmax() != mc_hist->GetXaxis()->GetXmax()"<<endl;
0039         exit(1);
0040         return nullptr;
0041     }
0042
0043     TH1D * h3 = (TH1D*)data_hist->Clone("reweight_hist");
0044     h3 -> Scale(1.0 / h3->Integral());
0045     mc_hist -> Scale(1.0 / mc_hist->Integral());
0046     h3->Divide(mc_hist); // note : data / mc
0047     return h3;
0048 }
0049
0050 int Run_GeoOffset(
0051     int process_id = 0,
0052     int run_num = -1,
0053     int nevents = -1,
0054     string input_directory = "/sphenix/user/ChengWei/INTT/INTT/general_codes/CWShih/INTTBcoResolution/macro", // note :
0055     string input_filename = "file_list_54280_intt.txt", // note : the ntuple
0056     string output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280/completed/BCO_check",
0057
0058     // todo : modify here
0059     std::string output_file_name_suffix = "_test1",
0060
0061     // todo : the data vertex Z distribution should be updated
0062     // std::string data_vtxZ_directory_file = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Jan172025/Run4/EvtVtxZ/completed/VtxZDist/completed/Data_vtxZDist_VtxZQA_EvtBcoFullDiffCut61_00054280_merged.root", // note : full directory
0063     // std::string data_vtxZ_hist_name = "h1D_INTTz_Inclusive70",
0064
0065     // todo: the map file should be updated
0066     std::string ColMulMask_map_dir = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/Run7/EvtVtxZ/ColumnCheck/baseline/completed/MulMap/completed",
0067     std::string ColMulMask_map_file = "MulMap_BcoFullDiffCut_Mbin70_VtxZ-30to30cm_ClusQAAdc35PhiSize40_00054280.root"
0068 )
0069 {
0070
0071     bool Global_GeoOffsetTag = (process_id == 0) ? false : true;
0072     int Global_ClusAdcCut = 35;
0073     int Global_ClusPhiSizeCut = 40;
0074
0075     bool is_tracklet_DeltaPhiCut_changed = true;
0076     std::pair<double,double> tracklet_DeltaPhiCut = {-0.04, 0.04};
0077
0078     std::string GeoOffset_string = (Global_GeoOffsetTag) ? "" : "_NoGeoOffset";
0079
0080     std::string job_index = std::to_string( process_id );
0081     int job_index_len = 5;
0082     job_index.insert(0, job_index_len - job_index.size(), '0');
0083
0084     std::string final_output_directory = output_directory + Form("/Run%s_%s", GeoOffset_string.c_str(), job_index.c_str());
0085     system(Form("if [ ! -d %s ]; then mkdir -p %s; fi", final_output_directory.c_str(), final_output_directory.c_str()));
0086
0087     std::pair<double,double> AvgVtxXY_MBD_vtxZ_cut = {-20, 20};
0088     std::pair<int,int> AvgVtxXY_INTTNClus_cut = {20, 350};
0089     int AvgVtxXY_ClusAdc_cut = Global_ClusAdcCut;
0090     int AvgVtxXY_ClusPhiSize_cut = Global_ClusPhiSizeCut;
0091
0092     bool AvgVtxXY_HaveGeoOffsetTag = Global_GeoOffsetTag;
0093     double AvgVtxXY_random_range_XYZ = 0.025; // note : unit : cm [250 um]
0094
0095     AvgVtxXY * avgXY = new AvgVtxXY(
0096         process_id,
0097         run_num,
0098         nevents,
0099         input_directory,
0100         input_filename,
0101         final_output_directory,
0102         output_file_name_suffix,
0103
0104         AvgVtxXY_MBD_vtxZ_cut,
0105         AvgVtxXY_INTTNClus_cut,
0106         AvgVtxXY_ClusAdc_cut,
0107         AvgVtxXY_ClusPhiSize_cut,
0108
0109         AvgVtxXY_HaveGeoOffsetTag,
0110         AvgVtxXY_random_range_XYZ
0111     );
0112
0113     std::string AvgXY_output_filename = avgXY->GetOutputFileName();
0114     cout<<"AvgXY: final_output_file_name: "<<AvgXY_output_filename<<endl;
0115     std::map<std::string, std::vector<double>> GeoOffset_map = avgXY->GetGeoOffsetMap();
0116
0117     avgXY -> PreparePairs();
0118
0119     // avgXY -> FindVertexQuadrant( 1, 0.4, {0, 0} ); // note : unit : cm
0120     // avgXY -> FindVertexLineFill({-0.022,0.2229}, 100, 0.25, 0.0001); // note : unit : cm
0121     // avgXY -> EndRun();
0122
0123     avgXY -> FindVertexQuadrant( 9, 0.4, {0, 0} ); // note : unit : cm
0124     avgXY -> FindVertexLineFill(avgXY -> GetVertexQuadrant(), 100, 0.25, 0.0001); // note : unit : cm
0125
0126     std::pair<double, double> vertexXYIncm = avgXY -> GetFinalAvgVtxXY();
0127     cout<<"vertexXYIncm: "<<vertexXYIncm.first<<" "<<vertexXYIncm.second<<endl;
0128
0129     avgXY -> EndRun();
0130
0131     std::cout<<"test test: done with AvgVtxXY"<<std::endl;
0132
0133     // Division : --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0134
0135     bool evtZ_IsFieldOn = false;
0136     bool evtZ_IsDCACutApplied = true;
0137     std::pair<std::pair<double,double>,std::pair<double,double>> evtZ_DeltaPhiCutInDegree = {{-0.6,0.6},{-1000.,1000.}}; // note : in degree
0138     std::pair<std::pair<double,double>,std::pair<double,double>> evtZ_DCAcutIncm = {{-0.1,0.1},{-1000.,1000.}}; // note : in cm
0139     int evtZ_ClusAdcCut = Global_ClusAdcCut;
0140     int evtZ_ClusPhiSizeCut = Global_ClusPhiSizeCut;
0141
0142     bool evtZ_PrintRecoDetails = false;
0143     bool evtZ_DrawEvtVtxZ = (process_id < 10) ? true : false;
0144
0145     bool evtZ_RunInttBcoFullDiff = false;
0146     bool evtZ_RunVtxZReco = true;
0147     bool evtZ_RunTrackletPair = false;
0148     bool evtZ_RunTrackletPairRotate = false;
0149
0150     bool evtZ_HaveGeoOffsetTag = Global_GeoOffsetTag;
0151
0152     EvtVtxZProtoTracklet * evtVtxZ = new EvtVtxZProtoTracklet(
0153         process_id,
0154         run_num,
0155         nevents,
0156         input_directory,
0157         input_filename,
0158         final_output_directory,
0159         output_file_name_suffix,
0160
0161         vertexXYIncm,
0162         evtZ_IsFieldOn,
0163         evtZ_IsDCACutApplied,
0164         evtZ_DeltaPhiCutInDegree,
0165         evtZ_DCAcutIncm,
0166         evtZ_ClusAdcCut,
0167         evtZ_ClusPhiSizeCut,
0168
0169         evtZ_PrintRecoDetails,
0170         evtZ_DrawEvtVtxZ,
0171
0172         evtZ_RunInttBcoFullDiff,
0173         evtZ_RunVtxZReco,
0174         evtZ_RunTrackletPair,
0175         evtZ_RunTrackletPairRotate,
0176
0177         evtZ_HaveGeoOffsetTag
0178     );
0179
0180     std::string evtVtxZ_output_filename = evtVtxZ->GetOutputFileName();
0181     cout<<"evtVtxZ: final_output_file_name: "<<evtVtxZ_output_filename<<endl;
0182
0183     evtVtxZ -> SetGeoOffset(GeoOffset_map);
0184     evtVtxZ -> MainProcess();
0185
0186     // TH1D * h1D_MC_INTTRecoVtxZ = (TH1D*) (evtVtxZ -> GetINTTRecovtxZH1D())->Clone("h1D_MC_INTTRecoVtxZ");
0187     // TH1D * h1D_data_INTTRecoVtxZ = GetDataVtxZH1D(data_vtxZ_directory_file, data_vtxZ_hist_name);
0188     // TH1D * h1D_vtxZReweight = GetVtxZReweightH1D(h1D_data_INTTRecoVtxZ, h1D_MC_INTTRecoVtxZ);
0189
0190     // TFile * file_out = new TFile(Form("%s/%s", final_output_directory.c_str(), "VtxZ_hist.root"), "recreate");
0191     // std::cout<<"test1: h1D_MC_INTTRecoVtxZ->GetName(): "<<h1D_MC_INTTRecoVtxZ->GetName()<<std::endl;
0192     // std::cout<<"test1: h1D_MC_INTTRecoVtxZ->GetNbinsX(): "<<h1D_MC_INTTRecoVtxZ->GetNbinsX()<<std::endl;
0193     // h1D_MC_INTTRecoVtxZ -> Write("h1D_MC_INTTRecoVtxZ");
0194     // h1D_data_INTTRecoVtxZ -> Write("h1D_data_INTTRecoVtxZ");
0195     // h1D_vtxZReweight -> Write("h1D_vtxZReweight");
0196     // file_out -> Close();
0197
0198     // for (int i = 0; i < h1D_vtxZReweight->GetNbinsX(); i++){
0199     //     cout<<"h1D_vtxZReweight, index: "<<i+1
0200     //     <<", center: "<<h1D_vtxZReweight -> GetBinCenter(i+1)
0201     //     <<", data: "  <<Form("%.3f",h1D_data_INTTRecoVtxZ->GetBinContent(i+1))
0202     //     <<", MC: "    <<Form("%.3f",h1D_MC_INTTRecoVtxZ->GetBinContent(i+1))
0203     //     <<", ratio: " <<Form("%.3f",h1D_vtxZReweight->GetBinContent(i+1))<<endl;
0204     // }
0205
0206     std::cout<<"test test: closing in EvtVtxZProtoTracklet "<<std::endl;
0207     evtVtxZ -> EndRun(0); // note : do not close the file_in
0208
0209     std::cout<<"test test: done with EvtVtxZProtoTracklet"<<std::endl;
0210     sleep(30);
0211
0212     // Division : --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0213
0214     // std::pair<bool, TH1D*> Track_vtxZReweight = {true, h1D_vtxZReweight};
0215     std::pair<bool, TH1D*> Track_vtxZReweight = {false, nullptr};
0216     bool Track_BcoFullDiffCut = false;
0217     bool Track_INTT_vtxZ_QA = true;
0218     std::pair<bool, std::pair<double, double>> Track_isClusQA = {true, {Global_ClusAdcCut, Global_ClusPhiSizeCut}}; // note : {adc, phi size}
0219     bool Track_HaveGeoOffsetTag = Global_GeoOffsetTag;
0220     std::pair<bool, int> Track_SetRandomHits = {false, 0};
0221     bool Track_RandInttZ = false;
0222
0223     bool Track_ColMulMask = true;
0224
0225
0226     TrackletHistogramNew * THN = new TrackletHistogramNew(
0227         process_id,
0228         run_num,
0229         nevents,
0230         final_output_directory,
0231         evtVtxZ_output_filename,
0232         final_output_directory,
0233
0234         output_file_name_suffix,
0235         vertexXYIncm,
0236
0237         Track_vtxZReweight,
0238         Track_BcoFullDiffCut,
0239         Track_INTT_vtxZ_QA,
0240         Track_isClusQA,
0241         Track_HaveGeoOffsetTag,
0242         Track_SetRandomHits,
0243         Track_RandInttZ,
0244         Track_ColMulMask
0245     );
0246
0247     std::string Track_output_filename = THN->GetOutputFileName();
0248     cout<<"Track: final_output_file_name: "<<Track_output_filename<<endl;
0249
0250     THN -> SetGeoOffset(GeoOffset_map);
0251
0252     if (Track_ColMulMask){
0253         THN -> SetGoodColMap(
0254         GetGoodColMap(ColMulMask_map_dir, ColMulMask_map_file, THN->GetGoodColMapName())
0255         );
0256     }
0257
0258     THN -> MainProcess();
0259
0260     std::cout<<"test test: closing in TrackletHistogramNew "<<std::endl;
0261     THN -> EndRun();
0262
0263     std::cout<<"test test: done with TrackletHistogramNew"<<std::endl;
0264     sleep(30);
0265
0266     // Division : --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0267
0268     bool dNdeta_ApplyAlphaCorr = false;
0269     bool dNdeta_isTypeA = false;
0270     std::pair<double,double> dNdeta_cut_INTTvtxZ = {-10, 10};
0271     int dNdeta_SelectedMbin = 70; // note : 0, 1, ---- 10, 70, 100
0272
0273     bool dNdeta_ApplyGeoAccCorr = false;
0274
0275     PreparedNdEtaEach * pdndeta = new PreparedNdEtaEach(
0276         process_id,
0277         run_num,
0278         final_output_directory,
0279         Track_output_filename,
0280         final_output_directory,
0281         output_file_name_suffix,
0282
0283         dNdeta_ApplyAlphaCorr,
0284         dNdeta_ApplyGeoAccCorr,
0285
0286         dNdeta_isTypeA,
0287         dNdeta_cut_INTTvtxZ,
0288         dNdeta_SelectedMbin
0289     );
0290
0291     // todo: if change cut
0292     if (is_tracklet_DeltaPhiCut_changed) {pdndeta->SetBkgRotated_DeltaPhi_Signal_range(tracklet_DeltaPhiCut);}
0293
0294     std::vector<std::string> dNdeta_output_filename = pdndeta->GetOutputFileName();
0295     cout<<"dNdeta: final_output_file_name: "<<dNdeta_output_filename[0]<<endl;
0296     cout<<"dNdeta: final_output_file_name: "<<dNdeta_output_filename[1]<<endl;
0297
0298     pdndeta -> PrepareStacks();
0299     pdndeta -> DoFittings();
0300     pdndeta -> PrepareMultiplicity();
0301     pdndeta -> PreparedNdEtaHist();
0302     pdndeta -> DeriveAlphaCorrection();
0303
0304     std::cout<<"test test: closing in PreparedNdEtaEach "<<std::endl;
0305     pdndeta -> EndRun();
0306
0307     std::cout<<"test test: done with PreparedNdEtaEach"<<std::endl;
0308
0309     // Division : --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0310
0311     std::cout<<"test test: closing in AvgVtxXY "<<std::endl;
0312
0313     std::cout<<"test test: done with in AvgVtxXY "<<std::endl;
0314
0315     sleep(30);
0316
0317     system(Form("mv %s %s/completed", final_output_directory.c_str(), output_directory.c_str()));
0318
0319     return 0;
0320 }