macros/MvtxProduction/Fun4All_Mvtx_hot_pixel_masking.C

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