macros/MvtxProduction/Fun4All_Mvtx_hot_pixel_masking.C

0001 // leave the GlobalVariables.C at the beginning, an empty line afterwards
0002 // protects its position against reshuffling by clang-format
0003 #include <GlobalVariables.C>
0004
0005 #include <mvtxcalib/MvtxFakeHitRate.h>
0006
0007 #include <mvtx/MvtxPixelDefs.h>
0008
0009 #include <trackbase/MvtxDefs.h>
0010 #include <trackbase/TrkrDefs.h>
0011
0012 #include <cdbobjects/CDBTTree.h>
0013
0014 #include <fun4all/Fun4AllInputManager.h>
0015 #include <fun4all/Fun4AllServer.h>
0016
0017 #include <fun4allraw/Fun4AllStreamingInputManager.h>
0018 #include <fun4allraw/InputManagerType.h>
0019 #include <fun4allraw/SingleMvtxPoolInput.h>
0020
0021 #include <phool/recoConsts.h>
0022
0023 #include <TFile.h>
0024 #include <TSystem.h>
0025 #include <TTree.h>
0026
0027 #include <nlohmann/json.hpp>
0028
0029 #include <iomanip>
0030 #include <fstream>
0031 #include <string>
0032
0033 R__LOAD_LIBRARY(libfun4all.so)
0034 R__LOAD_LIBRARY(libfun4allraw.so)
0035 R__LOAD_LIBRARY(libffamodules.so)
0036 R__LOAD_LIBRARY(libffarawmodules.so)
0037 R__LOAD_LIBRARY(libmvtx.so)
0038 R__LOAD_LIBRARY(libmvtxcalib.so)
0039 R__LOAD_LIBRARY(libcdbobjects.so)
0040
0041 using json = nlohmann::json;
0042
0043 std::map<std::string, int> felix_map{
0044 {"L0_03", 0},
0045 {"L0_04", 0},
0046 {"L2_01", 0},
0047 {"L2_02", 0},
0048 {"L0_05", 0},
0049 {"L2_03", 0},
0050 {"L2_04", 0},
0051 {"L2_15", 0},
0052 {"L1_00", 1},
0053 {"L1_01", 1},
0054 {"L2_06", 1},
0055 {"L2_07", 1},
0056 {"L1_02", 1},
0057 {"L1_03", 1},
0058 {"L2_08", 1},
0059 {"L2_09", 1},
0060 {"L0_00", 2},
0061 {"L0_01", 2},
0062 {"L1_04", 2},
0063 {"L2_10", 2},
0064 {"L0_02", 2},
0065 {"L1_05", 2},
0066 {"L1_06", 2},
0067 {"L1_07", 2},
0068 {"L0_09", 3},
0069 {"L0_10", 3},
0070 {"L2_05", 3},
0071 {"L2_11", 3},
0072 {"L0_11", 3},
0073 {"L2_12", 3},
0074 {"L2_13", 3},
0075 {"L2_14", 3},
0076 {"L0_06", 4},
0077 {"L0_07", 4},
0078 {"L1_12", 4},
0079 {"L2_00", 4},
0080 {"L0_08", 4},
0081 {"L1_13", 4},
0082 {"L1_14", 4},
0083 {"L1_15", 4},
0084 {"L1_08", 5},
0085 {"L1_09", 5},
0086 {"L2_16", 5},
0087 {"L2_17", 5},
0088 {"L1_10", 5},
0089 {"L1_11", 5},
0090 {"L2_18", 5},
0091 {"L2_19", 5}
0092 };
0093
0094 //https://stackoverflow.com/questions/997946/how-to-get-current-time-and-date-in-c
0095 std::string getDate()
0096 {
0097     std::time_t t = std::time(nullptr);   // get time now
0098     std::tm* now = std::localtime(&t);
0099
0100     std::stringstream date;
0101     date << (now->tm_year + 1900) << (now->tm_mon + 1) << now->tm_mday;
0102     return date.str();
0103 }
0104
0105 // https://stackoverflow.com/questions/478898/how-do-i-execute-a-command-and-get-the-output-of-the-command-within-c-using-po
0106 std::string exec(const char *cmd)
0107 {
0108   std::array<char, 128> buffer{};
0109   std::string result;
0110   std::unique_ptr<FILE, decltype(&pclose)> pipe(popen(cmd, "r"), pclose);
0111   if (!pipe)
0112   {
0113     throw std::runtime_error("popen() failed!");
0114   }
0115   while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr)
0116   {
0117     result += buffer.data();
0118   }
0119   return result;
0120 }
0121
0122 void calculate_FHR(const int run_number = 64639, const std::string& output_name = "output.root")
0123 {
0124   std::string calib_files = "/sphenix/lustre01/sphnxpro/physics/MVTX/calib/calib_mvtx";
0125   std::stringstream nice_run_number;
0126   nice_run_number << std::setw(8) << std::setfill('0') << std::to_string(run_number);
0127
0128   //Lets see how many FELIX servers took data for this run
0129   std::string felix_count_executable_command = "ls -1 ";
0130   felix_count_executable_command += calib_files;
0131   felix_count_executable_command += "*";
0132   felix_count_executable_command += nice_run_number.str();
0133   felix_count_executable_command += "-0000.evt | wc -l";
0134   int felix_count = std::stoi(exec(felix_count_executable_command.c_str()));
0135
0136 /*
0137   std::string strobe_length_executable_command = "psql -h sphnxdaqdbreplica daq --csv -c \"SELECT strobe FROM mvtx_strobe WHERE hostname = 'mvtx0' AND runnumber = ";
0138   strobe_length_executable_command += std::to_string(run_number);
0139   strobe_length_executable_command += ";\" | tail -n 1";
0140   float strobe_length = std::stof(exec(strobe_length_executable_command.c_str()));
0141 */
0142
0143   Fun4AllServer *se = Fun4AllServer::instance();
0144
0145   Fun4AllStreamingInputManager *in = new Fun4AllStreamingInputManager("Comb");
0146   in->Verbosity(1);
0147
0148   int counter = 0;
0149   for(int iflx = 0; iflx < felix_count; iflx++)
0150   {
0151     SingleMvtxPoolInput* mvtx_sngl = new SingleMvtxPoolInput("MVTX_FLX" + std::to_string(iflx));
0152     mvtx_sngl->Verbosity(1);
0153     mvtx_sngl->SetBcoRange(10);
0154     mvtx_sngl->SetNegativeBco(10);
0155     mvtx_sngl->runMVTXstandalone();
0156
0157     bool foundGoodFile = false;
0158     while (!foundGoodFile)
0159     {
0160       if (counter == 6) break;
0161
0162       std::string infile = calib_files + std::to_string(counter) + "-" +  nice_run_number.str() + "-0000.evt";
0163       std::ifstream file(infile.c_str());
0164       if (file.good())
0165       {
0166         mvtx_sngl->AddFile(infile);
0167         foundGoodFile = true;
0168       }
0169       else
0170       {
0171         ++counter;
0172       }
0173     }
0174
0175     in->registerStreamingInput(mvtx_sngl, InputManagerType::MVTX);
0176     ++counter;
0177   }
0178   se->registerInputManager(in);
0179
0180   MvtxFakeHitRate *mfhr = new MvtxFakeHitRate();
0181   mfhr->SetOutputfile(output_name);
0182   mfhr->SetMaxMaskedPixels(100000);
0183   se->registerSubsystem(mfhr);
0184 }
0185
0186 void generate_JSON_mask(const std::string &calibration_file = "output.root",
0187                         const std::string& output_file = "masked_pixels",
0188                         const double target_threshold = 3.2e-9)
0189 {
0190   TFile* f = new TFile(calibration_file.c_str(), "READ");
0191   TTree* t = dynamic_cast<TTree*>(f->Get("masked_pixels"));
0192   int num_strobes = 0;
0193   int num_masked_pixels = 0;
0194   double noise_threshold = 0.0;
0195   std::vector<uint64_t> * masked_pixels = nullptr;
0196   t->SetBranchAddress("num_strobes", &num_strobes);
0197   t->SetBranchAddress("num_masked_pixels", &num_masked_pixels);
0198   t->SetBranchAddress("noise_threshold", &noise_threshold);
0199   t->SetBranchAddress("masked_pixels", &masked_pixels);
0200
0201   int nentries = t->GetEntries();
0202   int selected_entry = -1;
0203   for (int i = 0; i < nentries; i++)
0204   {
0205     t->GetEntry(i);
0206     if(noise_threshold < target_threshold)
0207     {
0208       selected_entry = i;
0209       break;
0210     }
0211   }
0212
0213   t->GetEntry(selected_entry);
0214
0215   json masked_pixels_file;
0216   json masked_pixels_file_by_felix[6];
0217
0218   for (int i = 0; i < num_masked_pixels; i++)
0219   {
0220     uint64_t pixel = masked_pixels->at(i);
0221     uint32_t hitsetkey = MvtxPixelDefs::get_hitsetkey(pixel);
0222     uint32_t hitkey = MvtxPixelDefs::get_hitkey(pixel);
0223
0224     uint8_t layer = TrkrDefs::getLayer(hitsetkey);
0225     uint8_t stave = MvtxDefs::getStaveId(hitsetkey);
0226     uint8_t chip = MvtxDefs::getChipId(hitsetkey);
0227     uint16_t column = MvtxDefs::getCol(hitkey);
0228     uint16_t row = MvtxDefs::getRow(hitkey);
0229
0230     std::stringstream stave_sw_name;
0231     stave_sw_name << "L" << std::to_string(layer) << "_" << std::setw(2) << std::setfill('0') << std::to_string(stave);
0232     json pixel_info = json::array({column, row});
0233     masked_pixels_file[stave_sw_name.str()][std::to_string(chip)] += pixel_info;
0234     masked_pixels_file_by_felix[felix_map[stave_sw_name.str()]][stave_sw_name.str()][std::to_string(chip)] += pixel_info;
0235   }
0236
0237   f->Close();
0238   delete f;
0239
0240   std::cout << "Number of masked pixels: " << num_masked_pixels << std::endl;
0241   std::cout << "Noise threshold: " << noise_threshold << std::endl;
0242
0243   std::cout << "Writing masked pixels to json file" << std::endl;
0244
0245   std::ofstream json_output(output_file);
0246   json_output << masked_pixels_file.dump(2);
0247   json_output.close();
0248
0249   std::string output_file_no_trailer = output_file;
0250   std::string trailer = ".json";
0251   size_t pos = output_file_no_trailer.find(trailer);
0252   if (pos != std::string::npos)
0253   {
0254     output_file_no_trailer.erase(pos, trailer.length());
0255   }
0256
0257   for (int i = 0; i < 6; ++i)
0258   {
0259     std::string individual_mask_name = output_file_no_trailer + "_mvtx" + std::to_string(i) + ".json";
0260     std::ofstream individual_json_output(individual_mask_name);
0261     individual_json_output << masked_pixels_file_by_felix[i].dump(2);
0262     individual_json_output.close();
0263   }
0264 }
0265
0266 void generate_CDB_mask(const std::string& input_file = "masked_pixels.json")
0267 {
0268   std::string output_file = input_file;
0269   std::string trailer = ".json";
0270   size_t pos = output_file.find(trailer);
0271   if (pos != std::string::npos)
0272   {
0273     output_file.erase(pos, trailer.length());
0274   }
0275
0276   output_file += ".root";
0277
0278   CDBTTree *cdbttree = new CDBTTree(output_file);
0279
0280   unsigned int nStave[3] = {12, 16, 20};
0281
0282 //  char *calibrationsroot = getenv("CALIBRATIONROOT");
0283
0284   std::ifstream f(input_file);
0285   json data = json::parse(f);
0286   int PixelIndex = 0;
0287
0288   for (unsigned int layer = 0; layer < 3; ++layer)
0289   {
0290     for (unsigned int stave = 0; stave < nStave[layer]; ++stave)
0291     {
0292       std::stringstream staveStream;
0293       staveStream << "L" << layer << "_" << std::setw(2) << std::setfill('0') << stave;
0294       std::string staveName = staveStream.str();
0295
0296       bool doesStaveExist = data.contains(staveName);
0297       if (!doesStaveExist) continue;
0298
0299       for (unsigned int chip = 0; chip < 9; ++chip)
0300       {
0301         bool doesChipExist = data.at(staveName).contains(std::to_string(chip));
0302         if (!doesChipExist) continue;
0303
0304         std::vector<std::pair<int, int>> dead_pixels = data.at(staveName).at(std::to_string(chip));
0305
0306         for (auto &pixel : dead_pixels)
0307         {
0308           cdbttree->SetIntValue(PixelIndex,"layer",layer);
0309           cdbttree->SetIntValue(PixelIndex,"stave",stave);
0310           cdbttree->SetIntValue(PixelIndex,"chip",chip);
0311           cdbttree->SetIntValue(PixelIndex,"col",pixel.first);
0312           cdbttree->SetIntValue(PixelIndex,"row",pixel.second);
0313           PixelIndex = PixelIndex + 1;
0314         }
0315       }
0316     }
0317   }
0318
0319   cdbttree->Commit();
0320   cdbttree->SetSingleIntValue("TotalHotPixels",PixelIndex);
0321   cdbttree->CommitSingle();
0322   cdbttree->WriteCDBTTree();
0323   delete cdbttree;
0324 }
0325
0326 void Fun4All_Mvtx_hot_pixel_masking(const int nEvents = 5e5, const int run_number = 64639, const float threshold = 3.2e-9, bool doFHR = true, bool doJSONmask = true, bool doCDBmask = true)
0327 {
0328     std::string date = getDate();
0329     std::string run_string = std::to_string(run_number);
0330     std::string outputDir = "FHR_ana_run_" + run_string + "_" + date;
0331     std::string makeDirectory = "mkdir -p " + outputDir;
0332     system(makeDirectory.c_str());
0333
0334     Fun4AllServer *se = Fun4AllServer::instance();
0335     se->Verbosity(1);
0336
0337     recoConsts *rc = recoConsts::instance();
0338     Enable::CDB = true;
0339     rc->set_StringFlag("CDB_GLOBALTAG", CDB::global_tag);
0340     rc->set_uint64Flag("TIMESTAMP", run_number);
0341     rc->set_IntFlag("RUNNUMBER", run_number);
0342
0343     std::string run_info = "nEvents_" + std::to_string(nEvents) + "_run_" + run_string + "_date_" + date;
0344     std::string FHR_output_name = outputDir + "/FHR_output_" + run_info + ".root";
0345     std::string JSON_output_name = outputDir + "/masked_pixels_" + run_info + ".json";
0346
0347     if (doFHR)
0348     {
0349       calculate_FHR(run_number, FHR_output_name);
0350       se->run(nEvents);
0351       se->End();
0352       delete se;
0353     }
0354
0355     if (doJSONmask) generate_JSON_mask(FHR_output_name, JSON_output_name, threshold);
0356     if (doCDBmask) generate_CDB_mask(JSON_output_name);
0357     gSystem->Exit(0);
0358 }