analysis/INTTClusZSize_Tomoya/chfit_read_tapple.cc

0001 #include <iostream>
0002 #include <random>
0003
0004 TH1D* nCp(const std::vector<Double_t>& pro_list,int Event);
0005 int chip_reaction(Double_t probability);
0006
0007 int chfit_read_tapple()
0008 {
0009
0010 // TFile *file = TFile::Open("/sphenix/u/tomoya/work/24.05/own_dst_ana/macro/data/all/test602_v14_41620_all.root");//test602_v6_43409.root");
0011 TFile *file = TFile::Open("/sphenix/u/tomoya/work/24.05/own_dst_ana/macro/data/draw_macor/result/out20869_all_v2.root");
0012
0013 // TFile *file = TFile::Open("/sphenix/u/tomoya/work/24.05/own_dst_ana/macro/data/test602_v4.root");
0014 TH1D* hist = (TH1D*)file->Get("hist_z_layyer_0_ladder_0");
0015
0016  TH1D* hist_2 = (TH1D*)file->Get("hist_Zsize_Scale_layyer_0_ladder_0");
0017
0018   if (!hist)
0019   {
0020     std::cerr << "Unable to get ntuple" << std::endl;
0021     file->Close();
0022      return 1;
0023   }
0024
0025 int num_chip =0;
0026 std::vector<Double_t>pro_list;
0027  for(int i=0 ;i <401;i++)
0028 {
0029
0030   Double_t entry=hist->GetBinContent(i);
0031   if(entry>300000)
0032   {
0033
0034   cout<<"entry="<<entry<<endl;
0035   Double_t  pro = entry/550108/128;
0036 pro_list.push_back(pro);
0037
0038    cout<<"pro="<<pro<<endl;
0039   }
0040 }
0041
0042
0043  TCanvas *c1 = new TCanvas("c1", "Canvas Example", 800, 600);
0044 //  c1->Divide(3,2);
0045 //  c1->cd(1);
0046 // hist_2->Scale(1/5);
0047   hist_2->SetLineColor(kRed);
0048   hist_2->Draw("HIST");
0049
0050  TFile *td_out = new TFile("G609Ghypo_run20869_100000.root","RECREATE") ;
0051 //  cout<<"write"<<endl;
0052 //  histZ_all_layer->Write();
0053 //  outfile->Close();
0054
0055 std::vector<Double_t>pro_list_1(26,0.008);
0056
0057 TH1D* hist_m = nCp(pro_list_1,100000);
0058 hist_m->SetLineColor(kGreen);
0059 td_out->WriteTObject( hist_m, "hist_m" );
0060 hist_m->Draw("SAME");
0061
0062 //td_out->WriteTObject( hist_m, "hist_m" );
0063 hist_m->Scale(1.0/100000);
0064
0065 td_out->WriteTObject( hist_m, "hist_m_scale" );
0066 td_out->WriteTObject( hist_2, "hist_scale_Zsize_ladder0_layyer0" );
0067
0068
0069 return 0;
0070 }
0071
0072
0073 int chip_reaction(double probability)
0074 {
0075
0076         std::random_device seed_gen;
0077         std::default_random_engine engine(seed_gen());
0078
0079         std::binomial_distribution<> dist(1, probability);
0080
0081         int result = dist(engine);
0082
0083 return result;
0084 }
0085
0086
0087
0088 TH1D* nCp(const std::vector<Double_t>& pro_list, int Event)
0089 {
0090     std::vector<int> chip_list;
0091     std::vector<int> name_fire_chip;
0092     TH1D* hist_hypo = new TH1D("hist_hypo", "hist_hypo", 25, 1, 26);
0093
0094     for (int i = 0; i < Event; i++)
0095     {
0096 //        cout<<"i="<<i<<endl;
0097         int b_first = 0;
0098         int b_count = 0;
0099         std::vector<std::vector<float>> name_fire_channel(28);
0100
0101         chip_list.clear();
0102         name_fire_chip.clear();
0103
0104         for (int chip = 0; chip < 26; chip++)
0105         {
0106             for (int chp = 0; chp < 128; chp++)
0107             {
0108                 int b = chip_reaction(pro_list[chip]);
0109 //      cout<<"chip"<<chip<<"chp"<<chp<<"b"<<b<<endl;
0110                 chip_list.push_back(b);
0111                 if (b == 1)
0112                 {
0113                   //  cout<<"chip"<<chip<<"chp"<<chp<<"b"<<b<<endl;
0114                     b_first = b;
0115                     b_count++;
0116                     name_fire_channel[chip].push_back(chp);
0117                 }
0118             }
0119
0120             if (b_first == 1)
0121             {
0122                 name_fire_chip.push_back(chip);
0123                 b_first = 0;
0124             }
0125         }
0126
0127 int consecutiveCount = 1;
0128         int time_fill = 0;
0129
0130         for (int t = 0; t < name_fire_chip.size() ; t++)
0131         {
0132             bool shouldBreak = false;
0133
0134             if (name_fire_chip[t] == name_fire_chip[t + 1] - 1)
0135             {
0136                 for (size_t l = 0; l < name_fire_channel[name_fire_chip[t]].size(); l++)
0137                 {
0138                     for (size_t p = 0; p < name_fire_channel[name_fire_chip[t + 1]].size(); p++)
0139                     {
0140                         if (name_fire_channel[name_fire_chip[t]][l] == name_fire_channel[name_fire_chip[t + 1]][p] ||
0141                             name_fire_channel[name_fire_chip[t]][l] - 1 == name_fire_channel[name_fire_chip[t + 1]][p] ||
0142                             name_fire_channel[name_fire_chip[t]][l] + 1 == name_fire_channel[name_fire_chip[t + 1]][p])
0143                         {
0144                             consecutiveCount++;
0145                             shouldBreak = true;
0146                             break;
0147                         }
0148                         else if (name_fire_chip[t] == name_fire_chip[t - 1] + 1 && consecutiveCount > 1)
0149                         {
0150                             hist_hypo->Fill(consecutiveCount);
0151                     //        cout<<"consecutiveCount"<<consecutiveCount<<endl;
0152                             time_fill += consecutiveCount;
0153                             b_count -= consecutiveCount;
0154                             consecutiveCount = 1;
0155                             shouldBreak = true;
0156                             break;
0157                         }
0158                     }
0159                     if (shouldBreak) break;
0160                 }
0161             }
0162             else if (t==name_fire_chip.size()-1)
0163             {
0164                hist_hypo->Fill(consecutiveCount);
0165 //     cout<<"consecutiveCount_2"<<consecutiveCount<<endl;
0166
0167                time_fill= time_fill+consecutiveCount;
0168                consecutiveCount =1;
0169             }
0170             else
0171             {
0172                        hist_hypo->Fill(1, name_fire_channel[name_fire_chip[t]].size());
0173   //                cout<<"consecutiveCount_3,1"<<endl;
0174
0175                 time_fill += name_fire_channel[name_fire_chip[t]].size();
0176             }
0177         }
0178     }
0179     return hist_hypo;
0180 }
0181