QA/Intt/InttMixupQA.cc

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