Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-17 09:21:17

0001 #include "InttMixupQA.h"
0002 
0003 #include <ffarawobjects/InttRawHitContainer.h>
0004 
0005 #include <qautils/QAHistManagerDef.h>
0006 
0007 #include <ffarawobjects/InttRawHit.h>
0008 
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>
0012 
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>  // PHWHERE
0015 
0016 #include <TCanvas.h>
0017 #include <TFile.h>
0018 #include <TGraph.h>
0019 #include <TH1.h>
0020 #include <TH2.h>
0021 #include <TLegend.h>
0022 #include <TLine.h>
0023 #include <TROOT.h>
0024 #include <TStyle.h>
0025 #include <TSystem.h>
0026 
0027 #include <algorithm>
0028 #include <cstdint>
0029 #include <cstdlib>
0030 #include <format>
0031 #include <fstream>
0032 #include <iostream>
0033 #include <limits>
0034 #include <map>
0035 #include <sstream>
0036 #include <string>
0037 #include <utility>
0038 
0039 InttMixupQA::InttMixupQA(const std::string &name, const int run_num, const int felix_num)
0040   : SubsysReco(name)
0041   , run_num_(run_num)
0042   , felix_num_(felix_num)
0043 {
0044   std::cout << "felix_num_=" << felix_num_ << "felix_num=" << felix_num << std::endl;
0045 }
0046 
0047 int InttMixupQA::Init(PHCompositeNode * /*topNode*/)
0048 {
0049   if (Verbosity() > 5)
0050   {
0051     std::cout << "Beginning Init in InttMixupQA" << std::endl;
0052   }
0053 
0054   return 0;
0055 }
0056 
0057 int InttMixupQA::InitRun(PHCompositeNode *topNode)
0058 {
0059   if (!topNode)
0060   {
0061     std::cout << "InttMixupQA::InitRun(PHCompositeNode* topNode)" << std::endl;
0062     std::cout << "\tCould not retrieve topNode; doing nothing" << std::endl;
0063 
0064     return 1;
0065   }
0066 
0067   std::cout << "felix_num_=" << felix_num_ << std::endl;
0068 
0069   // Initialize histograms
0070   for (int felix = 0; felix < kFelix_num_; felix++)
0071   {
0072     std::string name = std::format("{}_allmulti_intt{}", getHistoPrefix(), felix);
0073     std::string title = name + std::format("_Run{}", run_num_) + ": with clone cut";
0074     h_allmulti_[felix] = new TH1F(name.c_str(), title.c_str(), 200, 0, bin);
0075     h_allmulti_[felix]->SetXTitle("Multiplicity");
0076     // h_allmulti_[felix]->SetLineColor(felix+1);
0077 
0078     name = std::format("{}_allclone_intt{}", getHistoPrefix(), felix);
0079     title = name + std::format("_Run{}", run_num_);
0080     h_allclone_[felix] = new TH1F(name.c_str(), title.c_str(), 200, 0, 200);
0081     h_allclone_[felix]->SetXTitle("Clone multiplicity");
0082     h_allclone_[felix]->SetLineColor(felix + 1);
0083 
0084     // These next 4 will be plotted
0085     auto *hm = QAHistManagerDef::getHistoManager();
0086     if (!hm)
0087     {
0088       std::cerr << PHWHERE << "\n"
0089                 << "\tQAHistManagerDef::getHistoManager() returned null\n"
0090                 << "\texit(1)" << std::endl;
0091       gSystem->Exit(1);
0092       exit(1);
0093     }
0094 
0095     name = std::format("{}_bco_full&0x7F_prev_vs_bco_intt{}", getHistoPrefix(), felix);
0096     title = name + std::format("_Run{}", run_num_);
0097     h_vsprefull_bco[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0098     h_vsprefull_bco[felix]->SetXTitle("BCO");
0099     h_vsprefull_bco[felix]->SetYTitle("BCO_FULL previous event &0x7F");
0100     hm->registerHisto(h_vsprefull_bco[felix]);
0101 
0102     name = std::format("{}_bco_full&0x7F_vs_bco_intt{}", getHistoPrefix(), felix);
0103     title = name + std::format("_Run{}", run_num_);
0104     h_vsfull_bco[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0105     h_vsfull_bco[felix]->SetXTitle("BCO");
0106     h_vsfull_bco[felix]->SetYTitle("BCO_FULL &0x7F");
0107     hm->registerHisto(h_vsfull_bco[felix]);
0108 
0109     name = std::format("{}_bco_full&0x7F_vs_bco_intt{}", getHistoPrefix(), felix);
0110     title = name + std::format("_Run{}", run_num_);
0111     h_prefull_bco[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0112     h_prefull_bco[felix]->SetXTitle("BCO_FULL previous event - BCO");
0113     h_prefull_bco[felix]->SetMinimum(0);
0114     hm->registerHisto(h_prefull_bco[felix]);
0115 
0116     name = std::format("{}_bco_full&0x7F_bco_intt{}", getHistoPrefix(), felix);
0117     title = name + std::format("_Run{}", run_num_);
0118     h_full_bco[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0119     h_full_bco[felix]->SetXTitle("BCO_FULL - BCO");
0120     h_full_bco[felix]->SetMinimum(0);
0121     hm->registerHisto(h_full_bco[felix]);
0122 
0123     name = std::format("{}_bco_full&0x7F_prev_bco_intt_all{}", getHistoPrefix(), felix);
0124     title = name + std::format("_Run{}", run_num_);
0125     h_prefull_bco_all[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0126     h_prefull_bco_all[felix]->SetXTitle("BCO_FULL previous event - BCO");
0127     h_prefull_bco_all[felix]->SetMinimum(0);
0128 
0129     name = std::format("{}_bco_full&0x7F_prev_vs_bco_intt_all{}", getHistoPrefix(), felix);
0130     title = name + std::format("_Run{}", run_num_);
0131     h_vsprefull_bco_all[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0132     h_vsprefull_bco_all[felix]->SetXTitle("BCO");
0133     h_vsprefull_bco_all[felix]->SetYTitle("BCO_FULL previous event &0x7F");
0134 
0135     std::string const name1 = "Mixup Multiplicity" + std::to_string(felix);
0136     std::string const name2 = "Mixup/all Multiplcity" + std::to_string(felix);
0137     std::string const title1 = name1 + std::format("_Run{}", run_num_);
0138     std::string const title2 = name2 + std::format("_Run{}", run_num_);
0139     for (int p = 0; p < divimul; p++)
0140     {
0141       h_mixupmulti[felix][p] = new TH1F((name1 + std::format("_{}", p)).c_str(), (title1 + std::format("_{}", p)).c_str(), 200, 0, bin);
0142       h_divmul[felix][p] = new TH1F((title1 + std::format("_{}", p)).c_str(), (title2 + std::format("_{}", p)).c_str(), 200, 0, bin);
0143       ;
0144     }
0145 
0146     name = std::format("{}_Number of mixup hit{}", getHistoPrefix(), felix);
0147     title = name + std::format("_Run{}", run_num_);
0148     h_mixup[felix] = new TH1F(name.c_str(), title.c_str(), 400, 0, 400);
0149 
0150     name = std::format("{}_prev_allhit vs Nmixup{}", getHistoPrefix(), felix);
0151     title = name + std::format("_Run{}", run_num_);
0152     h_prevsNmix[felix] = new TH2F(name.c_str(), title.c_str(), divimul + 10, 0, divimul + 10, bin, 0, bin);
0153     h_prevsNmix[felix]->GetXaxis()->SetNdivisions(405);
0154     h_prevsNmix[felix]->GetYaxis()->SetNdivisions(405);
0155     h_prevsNmix[felix]->SetXTitle("Number of Mixup hits");
0156     h_prevsNmix[felix]->SetYTitle("Number of previous event hits");
0157 
0158     name = std::format("{}_bco_full_prev-bco No copyhit{}", getHistoPrefix(), felix);
0159     title = name + std::format("_Run{}", run_num_);
0160     h_nocopyhit[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0161     h_nocopyhit[felix]->SetXTitle("bco_full_prev-bco");
0162 
0163     name = std::format("{}_bco_full_prev-bco copyhit{}", getHistoPrefix(), felix);
0164     title = name + std::format("_Run{}", run_num_);
0165     h_copyhit[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0166     h_copyhit[felix]->SetXTitle("bco_full_prev-bco");
0167 
0168     name = std::format("{}_Mixup & copy hit {}", getHistoPrefix(), felix);
0169     title = name + std::format("_Run{}", run_num_);
0170     h_mixcopy[felix] = new TH1F(name.c_str(), title.c_str(), 50, 0, 50);
0171     h_mixcopy[felix]->SetXTitle("Number of Mixup copy hit");
0172 
0173     name = std::format("{}_Nmixup vs Nclone{}", getHistoPrefix(), felix);
0174     title = name + std::format("_Run{}", run_num_);
0175     h_mixvscopy[felix] = new TH2F(name.c_str(), title.c_str(), 100, 0, 100, 50, 0, 50);
0176     h_mixvscopy[felix]->SetXTitle("Nmixup");
0177     h_mixvscopy[felix]->SetYTitle("Ncopy");
0178 
0179     name = std::format("{}_Nmixup vs pre_Nhits{}", getHistoPrefix(), felix);
0180     title = name + std::format("_Run{}", run_num_);
0181     h_hitfra[felix] = new TH2F(name.c_str(), title.c_str(), bin, 0, bin, bin, 0, bin);
0182 
0183     name = std::format("{}_Nmixup vs pre_Nhits others 4bin{}", getHistoPrefix(), felix);
0184     title = name + std::format("_Run{}", run_num_);
0185     h_bghit[felix] = new TH2F(name.c_str(), title.c_str(), bin, 0, bin, bin, 0, bin);
0186 
0187     fNhit[felix].open(std::format("./txtfile/NHit_{}_intt{}.txt", run_num_, felix));
0188 
0189     bcopeak_file[felix] = bcopeak_dir_ + std::format("bco_000{}_intt{}.root", run_num_, felix);
0190   }
0191 
0192   std::string name = "bco_full - prev_bco_full";
0193   std::string title = name + std::format("_Run{}", run_num_);
0194   h_interval = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0195   h_interval->SetXTitle("bco_full - prev_bco_full");
0196 
0197   name = "bco_full - prev_bco_full Mixup";
0198   title = name + std::format("_Run{}", run_num_);
0199   h_mixinterval = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0200 
0201   name = "Interval normarize";
0202   title = name + std::format("_Run{}", run_num_);
0203   h_divinter = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0204   h_divinter->SetXTitle("bco_full - prev_bco_full");
0205 
0206   name = "bco_full";
0207   title = name + std::format("_Run{}", run_num_);
0208   h_bcofull_7 = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0209   h_bcofull_7->SetXTitle("bco_full &0x7F");
0210 
0211   name = "bco";
0212   title = name + std::format("_Run{}", run_num_);
0213   h_bco = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0214   h_bco->SetXTitle("bco");
0215 
0216   name = "felix vs NmixupEv";
0217   title = name + std::format("_Run{}", run_num_);
0218   // h_NmixEv=new TH2F(name.c_str(), title.c_str(),8,0,8,100000,0,100000);
0219   h_NmixEv = new TH1F(name.c_str(), title.c_str(), 8, 0, 8);
0220 
0221   name = "All Event";
0222   title = name + std::format("_Run{}", run_num_);
0223   h_AllEv = new TH1F(name.c_str(), title.c_str(), 8, 0, 8);
0224 
0225   g_evfraction = new TGraph(n);
0226 
0227   g_cloevfraction = new TGraph(n);
0228 
0229   g_hitfraction = new TGraph(n);
0230 
0231   g_copyfraction = new TGraph(n);
0232 
0233   // bcopeak file set
0234   // bcopeak_file=bcopeak_dir_+std::format("bco_000%d",run_num_)+this->GetFileSuffix()+".root";
0235 
0236   // hotchan file set
0237   // hotchan_file=hotchan_dir_+std::format("InttHotDeadMap_000%d_50",run_num_)+this->GetFileSuffix()+".txt";
0238 
0239   tf_output_ = new TFile(output_root_.c_str(), "RECREATE");
0240 
0241   GetBcopeak();
0242 
0243   Readpeak();
0244 
0245   // Hotchancut();
0246 
0247   return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249 
0250 int InttMixupQA::process_event(PHCompositeNode *topNode)
0251 {
0252   if (Verbosity() > 5)
0253   {
0254     std::cout << "prev_bcofull=" << prev_bcofull << " "
0255               << "pre_allhit=" << pre_allhit << std::endl;
0256   }
0257 
0258   // get event bco_full
0259   /* std::string node_name_intteventheader = "INTTEVENTHEADER";
0260    InttEventInfo *node_intteventheader_map_ = findNode::getClass<InttEventInfo>(topNode, node_name_intteventheader);
0261 
0262    if (!node_intteventheader_map_){
0263      cerr << PHWHERE << node_name_intteventheader << " node is missing." << std::endl;
0264      return Fun4AllReturnCodes::ABORTEVENT;
0265    }
0266    uint64_t longbco_full = node_intteventheader_map_->get_bco_full();
0267    uint64_t bco_full =longbco_full &0x7FU;
0268    uint64_t difevent_bcofull = (longbco_full &bit )-(long_prev_bcofull &bit);
0269    h_interval->Fill(difevent_bcofull);
0270    h_bcofull_7->Fill(bco_full);
0271 
0272    if(Verbosity()>5){
0273    std::cout<<"bco_full="<<bco_full<<std::endl;
0274    }*/
0275 
0276   // get raw hit
0277   std::string const m_InttRawNodeName = "INTTRAWHIT";
0278   InttRawHitContainer *inttcont = findNode::getClass<InttRawHitContainer>(topNode, m_InttRawNodeName);
0279   if (!inttcont)
0280   {
0281     std::cout << PHWHERE << std::endl;
0282     std::cout << "InttMixupQA::process_event(PHCompositeNode* topNode)" << std::endl;
0283     std::cout << "Could not get \"" << m_InttRawNodeName << "\" from Node Tree" << std::endl;
0284     std::cout << "Exiting" << std::endl;
0285     gSystem->Exit(1);
0286     exit(1);
0287   }
0288 
0289   uint64_t const longbco_full = (0 < inttcont->get_nhits()) ? inttcont->get_hit(0)->get_bco() : std::numeric_limits<uint64_t>::max();
0290   uint64_t const difevent_bcofull = (longbco_full & bit) - (long_prev_bcofull & bit);
0291   h_interval->Fill(difevent_bcofull);
0292   uint64_t const bco_full = longbco_full & 0x7FU;
0293   // std::cout<<"longbco_full="<<longbco_full<<std::endl;
0294 
0295   ievent_++;
0296 
0297   if (ievent_ == 1)
0298   {
0299     first_bcofull = longbco_full;
0300   }
0301 
0302   /*if(ievent_<500000){
0303     return Fun4AllReturnCodes::EVENT_OK;
0304   }*/
0305 
0306   if ((ievent_ % 100 == 0 && ievent_ < 1000) || ievent_ % 1000 == 0)
0307   {
0308     std::cout << "Process event #" << ievent_ << std::endl;
0309   }
0310 
0311   int const nhits = inttcont->get_nhits();
0312   if (Verbosity() > 5)
0313   {
0314     std::cout << "Nhits = " << nhits << std::endl;
0315   }
0316 
0317   std::map<int, int> map_hit;
0318   int nhit_fx[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0319   int ncln_fx[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0320 
0321   int Nmixup[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0322   int Nother[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0323   // int Nclone[kFelix_num_] ={0,0,0,0,0,0,0,0};
0324   int Ncopy[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0325   // Loop over all INTTRAWHIT
0326   for (int i = 0; i < nhits; i++)
0327   {
0328     InttRawHit *intthit = inttcont->get_hit(i);
0329 
0330     int const fnum = intthit->get_packetid() - 3001;  // packet id
0331     int const fchn = intthit->get_fee();              // module
0332     int const adc = intthit->get_adc();               // adc
0333     int const chip = intthit->get_chip_id();          // chip
0334     int const chan = intthit->get_channel_id();       // channel
0335     int const bco = intthit->get_FPHX_BCO();          // FPHX bco
0336     int const hitID = 100000000 * (fnum + 1) + 1000000 * fchn + 10000 * chip + 10 * chan + adc;
0337     // int hitID_h = 10000000 * (fnum+1) + 100000 * fchn + 1000 * chip + chan;
0338 
0339     // clone hit cut
0340     auto itrHit = map_hit.find(hitID);
0341     if (itrHit == map_hit.end())
0342     {
0343       map_hit.insert(std::make_pair(hitID, 0));
0344       if (Verbosity() > 5)
0345       {
0346         std::cout << hitID << " " << fnum + 1 << " " << fchn << " " << chip << " " << chan << " " << adc << " " << bco << std::endl;
0347       }
0348 
0349       // hot channel cut
0350       // auto itrHit_h=hotmap.find(hitID_h);
0351       // if (itrHit_h==hotmap.end())
0352       {
0353         nhit_fx[fnum]++;
0354 
0355         h_bco->Fill(bco);
0356 
0357         int bco_diff = bco_full - bco;
0358         // int bco_diff = bco - bco_full;
0359         if (bco_diff < 0)
0360         {
0361           bco_diff = bco_diff + 128;
0362         }
0363 
0364         int pre_bco_diff = prev_bcofull - bco;
0365         // int pre_bco_diff =bco - prev_bcofull ;
0366         if (pre_bco_diff < 0)
0367         {
0368           pre_bco_diff = pre_bco_diff + 128;
0369         }
0370 
0371         /*if(bco_diff>21 && bco_diff<80)//for Run46681
0372         {
0373           std::cout<<"bco_diff="<<bco_diff<<std::endl;
0374         }
0375         else*/
0376         {
0377           if (!bcopar_[fnum].empty() && !bcopar_[fnum].contains(bco_diff))
0378           {
0379             if (bcopar_[fnum].contains(pre_bco_diff))
0380             {
0381               Nmixup[fnum]++;
0382 
0383               auto pre_itrHit = premap_hit.find(hitID);
0384               if (pre_itrHit != premap_hit.end())
0385               {
0386                 h_copyhit[fnum]->Fill(pre_bco_diff);
0387                 Ncopy[fnum]++;  // mixup & copy hit
0388               }
0389               else
0390               {
0391                 // mixup & not copy hit
0392                 h_nocopyhit[fnum]->Fill(pre_bco_diff);
0393               }
0394             }
0395             else if (otbcopar_[fnum].contains(pre_bco_diff))
0396             {
0397               Nother[fnum]++;
0398             }
0399 
0400             h_prefull_bco[fnum]->Fill(pre_bco_diff);  // without this event collision hit
0401             h_vsprefull_bco[fnum]->Fill(bco, prev_bcofull);
0402           }
0403           h_prefull_bco_all[fnum]->Fill(pre_bco_diff);  // with this event collision hit
0404           h_vsprefull_bco_all[fnum]->Fill(bco, prev_bcofull);
0405 
0406           h_full_bco[fnum]->Fill(bco_diff);
0407           h_vsfull_bco[fnum]->Fill(bco, bco_full);
0408         }
0409       }
0410       /*else
0411       {
0412         std::cout<<"Hot channel"<<std::endl;
0413       }*/
0414     }
0415     else
0416     {
0417       ncln_fx[fnum]++;  // Number of clone hit
0418     }
0419 
0420     if (Verbosity() > 5)
0421     {
0422       std::cout << "bco_full=" << bco_full << " "
0423                 << "bco=" << bco << " " << std::endl;
0424     }
0425   }
0426 
0427   if (Verbosity() > 5)
0428   {
0429     for (int id = 0; id < 8; id++)
0430     {
0431       std::cout << "Felix#" << id << " "
0432                 << "Nmixup=" << Nmixup[id] << " "
0433                 << "Ncopy=" << Ncopy[id] << " "
0434                 << "Nclone=" << ncln_fx[id] << std::endl;
0435     }
0436   }
0437 
0438   for (int felix = 0; felix < kFelix_num_; felix++)
0439   {
0440     // h_allmulti_[felix]->Fill(nhit_fx[felix]);
0441     h_allmulti_[felix]->Fill(pre_allhit[felix]);
0442     h_allclone_[felix]->Fill(ncln_fx[felix]);
0443     h_mixup[felix]->Fill(Nmixup[felix]);
0444     h_mixcopy[felix]->Fill(Ncopy[felix]);
0445     h_hitfra[felix]->Fill(Nmixup[felix], pre_allhit[felix]);
0446     h_bghit[felix]->Fill(Nother[felix], pre_allhit[felix]);
0447     for (int p = 0; p < divimul; p++)
0448     {
0449       if (Nmixup[felix] == p + 1)
0450       // if (Nmixup[felix] > p * 4 && Nmixup[felix] <= (p + 1) * 4)
0451       {
0452         h_mixupmulti[felix][p]->Fill(pre_allhit[felix]);
0453 
0454         /*if (1500 < pre_allhit && pre_allhit < 2500)
0455         {
0456           // h_fullmixup[fnumber]->Fill(difevent_bcofull);
0457         }*/
0458       }
0459       else
0460       {
0461         continue;
0462       }
0463 
0464       h_divmul[felix][p]->Divide(h_mixupmulti[felix][p], h_allmulti_[felix]);
0465     }
0466 
0467     mixupfraction[felix] = (double(Nmixup[felix]) / double(pre_allhit[felix] + Nmixup[felix]));
0468     // copyfraction[felix] = (double(Ncopy[felix]) / double(Nmixup[felix])) * 100;
0469     // thisclonefraction[felix]=(double(Nthisclone)/double(map_hit.size()))*100;
0470 
0471     h_AllEv->Fill(felix);
0472     if (Nmixup[felix] > 0 && pre_allhit[felix] > 0)
0473     // if (Nmixup[felix] > 0 && pre_allhit[felix]>0 )
0474     {
0475       // std::cout<<"Nmixup="<<Nmixup[felix]<<" "<<"mixupfraction="<<mixupfraction[felix]<<std::endl;
0476       h_mixinterval->Fill(difevent_bcofull);
0477       h_divinter->Divide(h_mixinterval, h_interval);
0478 
0479       h_prevsNmix[felix]->Fill(Nmixup[felix], pre_allhit[felix]);
0480       NmixupEv[felix]++;
0481       h_NmixEv->Fill(felix);
0482 
0483       // mixupfraction[felix]= (double(Nmixup[felix]) / double(pre_allhit[felix] + Nmixup[felix])) * 100;
0484       // copyfraction[felix] = (double(Ncopy[felix]) / double(Nmixup[felix])) * 100;
0485       mixupfraction_sum[felix] = mixupfraction_sum[felix] + mixupfraction[felix];
0486       Nmixup_sum[felix] = Nmixup_sum[felix] + Nmixup[felix];
0487       pre_allhit_sum[felix] = pre_allhit_sum[felix] + pre_allhit[felix];
0488       if (Ncopy[felix] > 0)
0489       {
0490         NmixcopyEv[felix]++;
0491         copyfraction_sum[felix] = copyfraction_sum[felix] + copyfraction[felix];
0492       }
0493       if (Verbosity() > 5)
0494       {
0495         std::cout << "Felix=" << felix << " "
0496                   << "mixupfraction=" << mixupfraction[felix] << "Nmixup=" << Nmixup[felix] << "pre_allhit=" << pre_allhit[felix] << "Ncopy=" << Ncopy[felix] << "copyfraction=" << copyfraction[felix] << std::endl;
0497       }
0498     }
0499     // h_divinter->Divide(h_mixinterval,h_interval);
0500     if (felix == felix_num_)
0501     {
0502       // std::cout<<nhit_fx[felix]<<std::endl;
0503 
0504       fNhit[felix] << nhit_fx[felix] << std::endl;
0505     }
0506     pre_allhit[felix] = nhit_fx[felix];
0507     // Initialize pre_allhit array
0508 
0509     /* pre_allhit[felix] = 0;
0510 
0511      pre_allhit[felix] = std::count_if(map_hit.begin(), map_hit.end(),
0512      [felix](const std::pair<int, int>& hit) {
0513      int felix_id = hit.first / 100000000 - 1;
0514      return felix_id == felix;
0515      });*/
0516   }
0517 
0518   long_prev_bcofull = longbco_full;
0519   prev_bcofull = bco_full;
0520   premap_hit = map_hit;
0521   //  h_allmulti_[fnum]->Fill(nhits);
0522   return Fun4AllReturnCodes::EVENT_OK;
0523 }
0524 
0525 int InttMixupQA::End(PHCompositeNode * /*topNode*/)
0526 {
0527   // Mixupfraction();
0528 
0529   if (Verbosity() > 1)
0530   {
0531     std::cout << "Processing InttMixupQA done" << std::endl;
0532   }
0533 
0534   /*for (int felix=0;felix<8;felix++){
0535   std::cout<<"felix="<< NmixupEv[felix]<<std::endl;
0536   }*/
0537   this->DrawHists();
0538 
0539   /*for( auto& hist : h_allmulti_ ) {
0540     tf_output_->WriteTObject( hist, hist->GetName() );
0541   }
0542 
0543   tf_output_->Close();*/
0544 
0545   return Fun4AllReturnCodes::EVENT_OK;
0546 }
0547 
0548 int InttMixupQA::SetOutputDir(std::string const &dir)
0549 {
0550   output_dir_ = dir;
0551 
0552   std::string const run_num_str = std::string(8 - std::to_string(run_num_).size(), '0') + std::to_string(run_num_);
0553   std::string const fnum_str = std::to_string(felix_num_);
0554 
0555   output_root_ = output_dir_ + "root/" + output_basename_ + run_num_str + "intt" + fnum_str + ".root";
0556   output_pdf_ = output_dir_ + "plots/" + output_basename_ + run_num_str + "intt" + fnum_str + ".pdf";
0557 
0558   // output_root_ = output_dir_ + "root/" + output_basename_ + run_num_str +".root";
0559   // output_pdf_  = output_dir_ + "plots/" + output_basename_ + run_num_str +".pdf";
0560 
0561   output_txt_ = output_dir_ + "txtfile/";
0562 
0563   return Fun4AllReturnCodes::EVENT_OK;
0564 }
0565 
0566 /*int InttMixupQA::SetHistBin(std::string type)
0567 {
0568   if (type==p)
0569   {
0570     bin=400;
0571     bin2=200000;
0572     bit=0x1FFF;//13bit
0573   }
0574   else if(type==Au)
0575   {
0576     bin=6000;
0577     bin2=200000;
0578     bit=0x7FFFFF;//21bit
0579   }
0580   else
0581   {
0582       std::cout<<"No set bin type"<<std::endl;
0583 
0584       return 1;
0585   }
0586   return Fun4AllReturnCodes::EVENT_OK;
0587 }*/
0588 
0589 // Private functions
0590 void InttMixupQA::GetBcopeak()
0591 {
0592   // get bco finder file
0593   /*tf_bcopeak_=new TFile(bcopeak_file.c_str(),"READ");
0594   std::cout<<tf_bcopeak_<<std::endl;*/
0595 
0596   // get hist from INTT QA
0597   /*for(int i=0; i<8; i++){
0598    h2_bco_felix[i]=(TH2D*)gROOT->FindObject(std::format("h2_bco_felix_%d",i));
0599    hbco[i]=h2_bco_felix[i]->ProjectionY(std::format("hbco_%d",i));
0600 
0601    //default peakbin for can't found peak bin
0602    h2_bco_felix_sub[i]=(TH2D*)gROOT->FindObject(std::format("h2_bco_felix_cut%d",i));
0603    hbco_sub[i]=h2_bco_felix_sub[i]->ProjectionY(std::format("hbco_sub%d",i));
0604    int maxbin_sub=hbco_sub[i]->GetMaximumBin();
0605    DEFAULT_BCO_VALUE[i]=hbco_sub[i]->GetBinLowEdge(maxbin_sub);
0606    }*/
0607 
0608   // get hist from made by myself
0609   for (int i = 0; i < 8; i++)
0610   {
0611     tf_bcopeak_[i] = new TFile(bcopeak_file[i].c_str(), "READ");
0612     // std::cout<<tf_bcopeak_[i]<<std::endl;
0613     hbco[i] = dynamic_cast<TH1D *>(gROOT->FindObject(std::format("h2_bco_felix_{}", i).c_str()));
0614     if (!hbco[i])
0615     {
0616       if (Verbosity())
0617       {
0618         std::cerr << PHWHERE << std::endl;
0619       }
0620       continue;
0621     }
0622 
0623     // default peakbin for can't found peak bin
0624     int const maxbin_sub = hbco[i]->GetMaximumBin();
0625     DEFAULT_BCO_VALUE[i] = hbco[i]->GetBinLowEdge(maxbin_sub);
0626   }
0627 
0628   // check_bcohist
0629   int peakbco[8];
0630   double bg_mean[8];
0631   double bg_rms[8];
0632   /////////////////////
0633   // find peak and BG area
0634   for (int id = 0; id < 8; id++)
0635   {
0636     double const max = hbco[id]->GetMaximum();
0637     int const maxbin = hbco[id]->GetMaximumBin();
0638     peakbco[id] = hbco[id]->GetBinLowEdge(maxbin);
0639     std::cout << id << " " << max << " " << maxbin << " " << peakbco[id] << std::endl;
0640 
0641     hbcohist[id] = new TH1F(std::format("hbcohist_{}", id).c_str(), std::format("hbcohist_{}", id).c_str(), 200, 0, max * 1.1);
0642     int const N = hbco[id]->GetNbinsX();
0643     for (int ibin = 0; ibin < N; ibin++)
0644     {
0645       double const cont = hbco[id]->GetBinContent(ibin + 1);
0646       hbcohist[id]->Fill(cont);
0647     }
0648   }
0649 
0650   for (int id = 0; id < 8; id++)
0651   {
0652     double const max = hbco[id]->GetMaximum();
0653     hbcohist2[id] = new TH1F(std::format("hbcohist2_{}", id).c_str(), std::format("hbcohist2_{}", id).c_str(), 200, 0, max * 1.1);
0654 
0655     int const N = hbco[id]->GetNbinsX();
0656     double const mean = hbcohist[id]->GetMean();
0657     double const rms = hbcohist[id]->GetRMS();
0658 
0659     for (int ibin = 0; ibin < N; ibin++)
0660     {
0661       double const cont = hbco[id]->GetBinContent(ibin + 1);
0662       if (mean + 3 * rms > cont)
0663       {
0664         hbcohist2[id]->Fill(cont);
0665       }
0666     }
0667     bg_mean[id] = hbcohist2[id]->GetMean();
0668     bg_rms[id] = hbcohist2[id]->GetRMS();
0669   }
0670 
0671   ///////
0672   // check peak area w/ peak +-10bin
0673   {
0674     f_felixpeak.open(output_txt_ + std::format("bco_{}.txt", run_num_));
0675 
0676     for (int id = 0; id < 8; id++)
0677     {
0678       f_felixpeak << id << " : ";
0679       // double thre   = bg_mean[id] + 6*bg_rms[id];
0680       double const thre = 5 * bg_mean[id];
0681       ;
0682       int const maxbin = hbco[id]->GetMaximumBin();
0683 
0684       int const N = hbco[id]->GetNbinsX();
0685       bool peak_found = false;
0686 
0687       for (int ibin = 0; ibin < 10; ibin++)
0688       {
0689         int binid = ibin + maxbin - 10;
0690         if (binid < 1)
0691         {
0692           binid += N;  // binid+=128;
0693         }
0694         double const cont = hbco[id]->GetBinContent(binid + 1);
0695         std::cout << id << " :  " << ibin << " " << binid << " " << maxbin;
0696         if (cont > thre)
0697         {
0698           int const peakbin = hbco[id]->GetBinLowEdge(binid + 1);
0699           std::cout << " exceed : " << peakbin;
0700           f_felixpeak << peakbin << " ";
0701           peak_found = true;
0702         }
0703         std::cout << std::endl;
0704       }
0705       for (int ibin = 0; ibin < 10; ibin++)
0706       {
0707         int binid = ibin + maxbin;
0708         if (binid >= N)
0709         {
0710           binid -= N;  // binid+=128;
0711         }
0712         double const cont = hbco[id]->GetBinContent(binid + 1);
0713         std::cout << id << " :  " << ibin << " " << binid << " " << maxbin;
0714         if (cont > thre)
0715         {
0716           int const peakbin = hbco[id]->GetBinLowEdge(binid + 1);
0717           std::cout << " exceed : " << peakbin;
0718           f_felixpeak << peakbin << " ";
0719           peak_found = true;
0720         }
0721         std::cout << std::endl;
0722       }
0723       if (!peak_found)
0724       {
0725         f_felixpeak << DEFAULT_BCO_VALUE[id] << " ";
0726       }
0727 
0728       f_felixpeak << std::endl;
0729     }
0730     f_felixpeak.close();
0731   }
0732 
0733   // get min
0734   float minimum = 1000000;
0735   for (auto &id : hbco)
0736   {
0737     std::cout << id->GetMinimumBin() << std::endl;
0738     std::cout << id->GetBinContent(id->GetMinimumBin()) << std::endl;
0739     double const min = id->GetBinContent(id->GetMinimumBin());
0740     std::cout << min << std::endl;
0741     minimum = std::min<double>(min, minimum);
0742   }
0743 
0744   // Draw hist
0745 
0746   TCanvas *c1 = new TCanvas("c1", "c1", 1200, 700);
0747   c1->Divide(4, 4);
0748   for (int id = 0; id < 8; id++)
0749   {
0750     c1->cd(id + 1);
0751     gPad->SetLogy();
0752     hbco[id]->SetMinimum(minimum * 0.2);
0753     hbco[id]->Draw();
0754 
0755     TLine *la = new TLine(0, 5 * bg_mean[id], 128, 5 * bg_mean[id]);
0756     la->Draw("same");
0757     TLine *lb = new TLine(0, bg_mean[id] + 6 * bg_rms[id], 128, bg_mean[id] + 6 * bg_rms[id]);
0758     lb->SetLineColor(2);
0759     lb->Draw("same");
0760   }
0761 
0762   for (int id = 0; id < 8; id++)
0763   {
0764     c1->cd(id + 9);
0765     gPad->SetLogy();
0766     double const max = hbcohist[id]->GetMaximum();
0767 
0768     hbcohist[id]->Draw();
0769     hbcohist2[id]->SetLineColor(2);
0770     hbcohist2[id]->Draw("same");
0771 
0772     TLine *lm = new TLine(bg_mean[id], 0.001, bg_mean[id], max * 1.1);
0773     lm->Draw();
0774     TLine *l5m = new TLine(5 * bg_mean[id], 0.001, 5 * bg_mean[id], max * 1.1);
0775     l5m->Draw();
0776     TLine *l = new TLine(bg_mean[id] + 6 * bg_rms[id], 0.001, bg_mean[id] + 6 * bg_rms[id], max * 1.1);
0777     l->SetLineColor(2);
0778     l->Draw();
0779   }
0780 
0781   c1->Print((output_dir_ + "plots/" + std::format("check_bco_file_{}.pdf", run_num_)).c_str());
0782 }
0783 
0784 void InttMixupQA::Readpeak()
0785 {
0786   std::string const bcofile_ = (output_txt_ + std::format("bco_{}.txt", run_num_));
0787 
0788   if (!bcofile_.empty())
0789   {
0790     std::cout << "BCO file : " << bcofile_ << std::endl;
0791     std::ifstream file(bcofile_);
0792     if (!file.is_open())
0793     {
0794       std::cerr << "Failed to open file: " << bcofile_ << std::endl;
0795       return;
0796     }
0797 
0798     std::string sline;
0799 
0800     while (getline(file, sline))
0801     {
0802       // std::cout<<sline<<std::endl;
0803       std::istringstream ss(sline);
0804       std::string sbuf;
0805 
0806       int felix = -1;
0807       int idx = 0;
0808       while (getline(ss, sbuf, ':'))
0809       {
0810         // 1st for felix id
0811         if (idx == 0)
0812         {
0813           // std::cout<<idx<<", id "<<sbuf<<std::endl;
0814           felix = std::stoi(sbuf);
0815           if (felix < 0 || felix > 7)
0816           {
0817             std::cout << "felixid out of range. id=" << felix << std::endl;
0818             break;
0819           }
0820         }
0821         else
0822         {
0823           // std::cout<<idx<<" "<<felix<<", bco "<<sbuf<<std::endl;
0824           std::istringstream ss2(sbuf);
0825           std::string sbuf2;
0826           std::string sbuf3;
0827           std::string sbuf4;
0828           std::string sbuf5;
0829           std::string sbuf6;
0830           while (ss2 >> sbuf2)
0831           {
0832             // std::cout<<"   "<<sbuf2<<std::endl;
0833             bcopar_[felix].insert(std::stoi(sbuf2));
0834             int const peak = std::stoi(sbuf2);
0835             int const peak1 = peak + 1;
0836             int const peak2 = peak + 2;
0837             int const peak_1 = peak - 1;
0838             int const peak_2 = peak - 2;
0839             sbuf3 = std::to_string(peak1);
0840             sbuf4 = std::to_string(peak2);
0841             sbuf5 = std::to_string(peak_1);
0842             sbuf6 = std::to_string(peak_2);
0843             otbcopar_[felix].insert(std::stoi(sbuf3));
0844             otbcopar_[felix].insert(std::stoi(sbuf4));
0845             otbcopar_[felix].insert(std::stoi(sbuf5));
0846             otbcopar_[felix].insert(std::stoi(sbuf6));
0847           }
0848         }
0849         idx++;
0850       }
0851     }
0852     file.close();
0853   }
0854 
0855   std::cout << "BCO peaks " << std::endl;
0856   for (int i = 0; i < 8; i++)
0857   {
0858     std::cout << "    felix " << i << " : ";
0859     for (int const itr : bcopar_[i])
0860     {
0861       std::cout << itr << " ";
0862     }
0863     std::cout << std::endl;
0864   }
0865 
0866   std::cout << "Around BCO peaks " << std::endl;
0867   for (int i = 0; i < 8; i++)
0868   {
0869     std::cout << "    felix " << i << " : ";
0870     for (int const itr : otbcopar_[i])
0871     {
0872       std::cout << itr << " ";
0873     }
0874     std::cout << std::endl;
0875   }
0876 }
0877 
0878 /*void InttMixupQA::Hotchancut()
0879 {
0880   //tf_hotchan_=new TFile(hotchan_file.c_str(),"READ");
0881   std::ifstream tf_hotchan_(hotchan_file.c_str());
0882  // std::cout <<tf_hotchan_<< std::endl;
0883 
0884   //std::map<int,int> hotmap;
0885   std::string sline;
0886   int lineNumber=0;
0887 
0888   getline(tf_hotchan_, sline);
0889 
0890   while (getline(tf_hotchan_,sline))
0891   {
0892     std::istringstream ss(sline);
0893     int felix,felix_ch,chip,channel;
0894 
0895     ss>>felix>>felix_ch>>chip>>channel;
0896 
0897     int hotchanID=10000000 * (felix+1) + 100000 * felix_ch + 1000 * chip + channel;
0898 
0899     hotmap[lineNumber++]=hotchanID;
0900 
0901   }
0902 
0903   tf_hotchan_.close();
0904 
0905   for(const auto&[key,value]:hotmap){
0906     std::cout<<"Line"<<key+1<<":"<<value<<std::endl;
0907   }
0908 }*/
0909 
0910 void InttMixupQA::DrawHists()
0911 {
0912   std::cout << "output pdf: " << output_pdf_ << std::endl;
0913 
0914   /////////////////////////Mixup plot Draw//////////////////////////////
0915   TCanvas *c0 = new TCanvas("c0", "c0");
0916   c0->Print((output_pdf_ + "[").c_str());
0917   TCanvas *c1[8];
0918   TCanvas *c2[8];
0919   TCanvas *c3[8];
0920   // TCanvas *c4[8];
0921   TCanvas *c5[8];
0922   TCanvas *c6[8];
0923   TCanvas *c10[8];
0924   TCanvas *c7 = new TCanvas("canvas7", "interval", 1600, 1200);
0925   TCanvas *c8 = new TCanvas("canvas8", "fraction", 1600, 1200);
0926   TCanvas *c9 = new TCanvas("canvas9", "bco bco_full", 1600, 1200);
0927 
0928   gStyle->SetOptStat(0);
0929 
0930   for (int felix = 0; felix < kFelix_num_; felix++)
0931   {
0932     if (felix != felix_num_)
0933     {
0934       continue;
0935     }
0936     c1[felix] = new TCanvas(std::format("canvas1_{}", felix).c_str(), "bcofull&bco", 1600, 1200);
0937     c1[felix]->cd();
0938     gStyle->SetOptStat(0);
0939     c1[felix]->Divide(2, 2);
0940 
0941     c1[felix]->cd(1);
0942     h_vsfull_bco[felix]->Draw("colz");
0943     c1[felix]->cd(2);
0944     h_vsprefull_bco[felix]->Draw("colz");
0945     c1[felix]->cd(3);
0946     h_full_bco[felix]->SetMinimum(0);
0947     h_full_bco[felix]->Draw();
0948     c1[felix]->cd(4);
0949     h_prefull_bco[felix]->SetMinimum(0);
0950     h_prefull_bco[felix]->Draw();
0951 
0952     c2[felix] = new TCanvas(std::format("canvas2_{}", felix).c_str(), "Multiplicity", 1600, 1200);
0953     c2[felix]->cd();
0954     gStyle->SetOptStat(0);
0955     gPad->SetRightMargin(0.15);
0956     TLegend *legend = new TLegend(0.7, 0.45, 0.9, 0.9);
0957     gPad->SetLogy();
0958     h_allmulti_[felix]->GetXaxis()->SetNdivisions(405);
0959     h_allmulti_[felix]->GetYaxis()->SetNdivisions(405);
0960     // h_allmulti[i]->SetXTitle("Hit Multiplcity/event");
0961     int const a = h_allmulti_[felix]->GetMaximum();
0962     h_allmulti_[felix]->GetYaxis()->SetRangeUser(1, a * 1.6);
0963     h_allmulti_[felix]->SetXTitle("Number of hits of previous event");
0964     h_allmulti_[felix]->SetYTitle("Entry");
0965     h_allmulti_[felix]->Draw();
0966     h_allmulti_[felix]->SetLineColor(1);
0967     legend->SetTextSize(0.04);
0968     legend->AddEntry(h_allmulti_[felix], "Number of", "");
0969     legend->AddEntry(h_allmulti_[felix], "Mixup hits", "");
0970     legend->AddEntry(h_allmulti_[felix], "All event", "l");
0971 
0972     for (int p = 0; p < 10; p++)
0973     {
0974       h_mixupmulti[felix][p]->Draw("same");
0975       int const color = p + 2;
0976       if (color == 10)
0977       {
0978         h_mixupmulti[felix][p]->SetLineColor(46);
0979       }
0980       else
0981       {
0982         h_mixupmulti[felix][p]->SetLineColor(color);
0983       }
0984       legend->AddEntry(h_mixupmulti[felix][p], std::format("{}", p + 1).c_str(), "l");
0985     }
0986     legend->Draw();
0987 
0988     c3[felix] = new TCanvas(std::format("canvas3_{}", felix).c_str(), "divide multiplicity", 1600, 1200);
0989     c3[felix]->cd();
0990     gStyle->SetOptStat(0);
0991     gPad->SetRightMargin(0.15);
0992     TLegend *leg = new TLegend(0.75, 0.5, 0.9, 0.9);
0993     gPad->SetLogy();
0994     for (int p = 0; p < 10; p++)
0995     {
0996       h_divmul[felix][p]->SetXTitle("Mixup/all hit Multiplicity");
0997       if (p == 0)
0998       {
0999         h_divmul[felix][p]->Draw();
1000       }
1001       else
1002       {
1003         h_divmul[felix][p]->Draw("same");
1004       }
1005       int const color2 = p + 2;
1006       if (color2 == 10)
1007       {
1008         h_divmul[felix][p]->SetLineColor(46);
1009       }
1010       else
1011       {
1012         h_divmul[felix][p]->SetLineColor(color2);
1013       }
1014       leg->AddEntry(h_divmul[felix][p], std::format("{}", p + 1).c_str(), "l");
1015     }
1016 
1017     /*c4[felix]= new TCanvas( std::format("canvas4_{}",felix), "copy hit",1600, 1200 );
1018     c4[felix]->cd();
1019     gStyle->SetOptStat( 0 );
1020     c4[felix]->Divide(2,2);
1021     c4[felix]->cd(1);
1022     h_nocopyhit[felix]->Draw();
1023     c4[felix]->cd(2);
1024     h_copyhit[felix]->Draw();
1025     c4[felix]->cd(3);
1026     h_mixcopy[felix]->Draw();
1027     c4[felix]->cd(4);
1028     h_mixvscopy[felix]->Draw();*/
1029 
1030     c5[felix] = new TCanvas(std::format("canvas5_{}", felix).c_str(), "previous event hit vs Mixup hit", 1600, 1200);
1031     c5[felix]->cd();
1032     gStyle->SetOptStat(0);
1033     gPad->SetRightMargin(0.15);
1034     h_prevsNmix[felix]->Draw("colz");
1035 
1036     c6[felix] = new TCanvas(std::format("canvas6_{}", felix).c_str(), "Mixup", 1600, 1200);
1037     c6[felix]->cd();
1038     c6[felix]->Divide(2, 2);
1039     c6[felix]->cd(1);
1040     gStyle->SetOptStat(0);
1041     gPad->SetLogy();
1042     h_mixup[felix]->Draw();
1043 
1044     c6[felix]->cd(4);
1045     h_prefull_bco_all[felix]->Draw();
1046     c6[felix]->cd(2);
1047     h_vsprefull_bco_all[felix]->Draw();
1048 
1049     c10[felix] = new TCanvas(std::format("canvas10_{}", felix).c_str(), "for fraction");
1050     c10[felix]->cd();
1051     c10[felix]->Divide(2, 1);
1052     c10[felix]->cd(1);
1053     h_hitfra[felix]->Draw("colz");
1054     c10[felix]->cd(2);
1055     h_bghit[felix]->Draw("colz");
1056 
1057     c1[felix]->Print(output_pdf_.c_str());
1058     c2[felix]->Print(output_pdf_.c_str());
1059     c3[felix]->Print(output_pdf_.c_str());
1060     // c4[felix]->Print(output_pdf_.c_str());
1061     c5[felix]->Print(output_pdf_.c_str());
1062     c6[felix]->Print(output_pdf_.c_str());
1063     c10[felix]->Print(output_pdf_.c_str());
1064     tf_output_->cd();
1065     /*c1[felix]->Write();
1066     c2[felix]->Write();
1067     c3[felix]->Write();
1068     c4[felix]->Write();
1069     c5[felix]->Write();
1070     c6[felix]->Write();*/
1071 
1072     h_vsprefull_bco[felix]->Write();
1073     h_vsfull_bco[felix]->Write();
1074     h_prefull_bco[felix]->Write();
1075     h_full_bco[felix]->Write();
1076     /*for (int p = 0; p < 10; p++)
1077     {
1078       h_allmulti_[felix]->Add(h_mixupmulti[felix][p]);
1079     }
1080     h_allmulti_[felix]->Write();*/
1081     c2[felix]->Write();
1082     h_mixup[felix]->Write();
1083     h_prevsNmix[felix]->Write();
1084     h_hitfra[felix]->Write();
1085     h_bghit[felix]->Write();
1086     h_NmixEv->Write();
1087     h_AllEv->Write();
1088     h_divinter->Write();
1089     // tf_output_->Close();
1090   }
1091   // h_interval->Write();
1092   // h_mixinterval->Write();
1093 
1094   c7->cd();
1095   gStyle->SetOptStat(0);
1096   TLegend *legsph3 = new TLegend(.55, .4, .65, .6);
1097   gPad->SetRightMargin(0.15);
1098   h_interval->Draw();
1099   h_mixinterval->Draw("same");
1100   h_mixinterval->SetLineColor(2);
1101   legsph3->AddEntry(h_interval, "All event", "l");
1102   legsph3->AddEntry(h_mixinterval, "Mixup event", "l");
1103 
1104   /*c8->cd();
1105   gStyle->SetOptStat( 0 );
1106   c8->Divide(2, 2);
1107   c8->cd(1);
1108   g_evfraction->SetMarkerStyle(20);
1109   g_evfraction->SetMarkerSize(1.1);
1110   g_evfraction->SetMinimum(0);
1111   g_evfraction->SetTitle(std::format("Mixup event fraction Run{};Felix;Mixup Event fraction ",run_num_).c_str());
1112   g_evfraction->Draw("AP");
1113 
1114   c8->cd(2);
1115   g_hitfraction->SetMarkerStyle(20);
1116   g_hitfraction->SetMarkerSize(1.1);
1117   g_hitfraction->SetMinimum(0);
1118   g_hitfraction->SetTitle(std::format("Mixup Hit fraction Run{};Felix;Mixup Hit fraction ",run_num_).c_str());
1119   g_hitfraction->Draw("AP");
1120 
1121   c8->cd(3);
1122   g_cloevfraction->SetMarkerStyle(20);
1123   g_cloevfraction->SetMarkerSize(1.1);
1124   g_cloevfraction->SetMinimum(0);
1125   g_cloevfraction->SetTitle(std::format("Mixclone event fraction Run{}/Mixup;Felix;Mixclone event fraction ",run_num_).c_str());
1126   g_cloevfraction->Draw("AP");
1127 
1128   c8->cd(4);
1129   g_copyfraction->SetMarkerStyle(20);
1130   g_copyfraction->SetMarkerSize(1.1);
1131   g_copyfraction->SetMinimum(0);
1132   g_copyfraction->SetTitle(std::format("Mixclone hit fraction Run{}/Mixup;Felix;clone hit from prev/Mixup ",run_num_).c_str());
1133   g_copyfraction->Draw("AP");
1134 
1135   fgraph->cd();
1136   g_hitfraction->Write("g_hitfraction");
1137   g_evfraction->Write("g_evfraction");
1138   g_copyfraction->Write("g_copyfraction");
1139   g_cloevfraction->Write("g_cloevfraction");
1140   fgraph->Close();*/
1141 
1142   c9->cd();
1143   gStyle->SetOptStat(0);
1144   c9->Divide(2, 1);
1145   c9->cd(1);
1146   h_bcofull_7->Draw();
1147   c9->cd(2);
1148   h_bco->Draw();
1149 
1150   tf_output_->cd();
1151   c7->Write();
1152   /*c8->Write();
1153   c9->Write();*/
1154 
1155   /*g_hitfraction->Write("g_hitfraction");
1156   g_evfraction->Write("g_evfraction");
1157   g_copyfraction->Write("g_copyfraction");
1158   g_cloevfraction->Write("g_cloevfraction");*/
1159   tf_output_->Close();
1160 
1161   c7->Print((output_pdf_).c_str());
1162   c8->Print((output_pdf_).c_str());
1163   c9->Print((output_pdf_).c_str());
1164 
1165   c9->Print((output_pdf_ + "]").c_str());
1166 }