DST_MC/ReadMCselfgen/calculate_MBin.C

0001 #include "MCSelfGenReader.C"
0002 #include "../ana_map_folder/ana_map_v1.h"
0003
0004 // note : this code is mainly for the selfgen MC, to calculate the MBin cut on NClus
0005
0006 void calculate_MBin(string input_directory, string input_file_name, string output_directory)
0007 {
0008     TFile * file_in =  TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0009     TTree * tree = (TTree*)file_in->Get("tree");
0010     long long N_event = tree->GetEntries();
0011
0012     std::cout << "N_event = " << N_event << std::endl;
0013
0014     MCSelfGenReader * inttMCselfgen = new MCSelfGenReader(tree);
0015
0016     vector<int> NClus_event; NClus_event.clear();
0017     vector<int> NClus_bin_cut; NClus_bin_cut.clear();
0018
0019     TH1F * NClus_hist = new TH1F("",";N clusters;Entry",100,0,4500);
0020
0021     TLine * coord_line = new TLine();
0022     coord_line -> SetLineWidth(1);
0023     coord_line -> SetLineColor(2);
0024     // coord_line -> SetLineStyle(2);
0025
0026     vector<int> MBin_NClus_cut = {
0027         5,
0028         15,
0029         25,
0030         35,
0031         45,
0032         55,
0033         65,
0034         75,
0035         85,
0036         95,
0037     };
0038
0039     vector<string> centrality_region = ANA_MAP_V1::centrality_region;
0040
0041
0042     for (long long event_i = 0; event_i < N_event; event_i++)
0043     {
0044         inttMCselfgen->LoadTree(event_i);
0045         inttMCselfgen->GetEntry(event_i);
0046
0047         if (event_i % 3000 == 0) {cout<<"eID : "<<event_i<<" NClus : "<<inttMCselfgen->NClus<<endl;}
0048
0049         NClus_hist->Fill(inttMCselfgen->NClus);
0050         NClus_event.push_back(inttMCselfgen->NClus);
0051     }
0052
0053     int total_entry = NClus_event.size();
0054
0055     int sort_clu_index[NClus_event.size()];
0056     cout<<"NClus_event: "<<NClus_event.size()<<endl;
0057     TMath::Sort(int(NClus_event.size()), &NClus_event[0], sort_clu_index);
0058
0059
0060
0061     TFile * file_out = new TFile(Form("%s/%s", output_directory.c_str(), "MBin_NClus_cut.root"), "RECREATE");
0062     file_out -> cd();
0063     TTree * tree_out = new TTree("tree", "MBin cut on NClus");
0064     tree_out -> Branch("MBinCut", &NClus_bin_cut);
0065
0066     TCanvas * c1 = new TCanvas("","",900,800);
0067     NClus_hist -> Draw("hist");
0068     c1 -> SetLogy(1);
0069
0070     for (int i = 0; i < MBin_NClus_cut.size(); i++)
0071     {
0072         cout<<"centrality upper cut of portion : "<<centrality_region[i]<<", the N event ID : "<< int(total_entry * (double(MBin_NClus_cut[i]) /100.))<<" entry_line : "<<NClus_event[ sort_clu_index[ int(total_entry * (double(MBin_NClus_cut[i]) /100.)) ] ]<<endl;
0073
0074         NClus_bin_cut.push_back(NClus_event[ sort_clu_index[ int(total_entry * (double(MBin_NClus_cut[i]) /100.)) ] ]);
0075
0076         coord_line -> DrawLine(NClus_event[ sort_clu_index[ int(total_entry * (double(MBin_NClus_cut[i]) /100.)) ] ], 0, NClus_event[ sort_clu_index[ int(total_entry * (double(MBin_NClus_cut[i]) /100.)) ] ], NClus_hist -> GetBinContent( int( NClus_event[ sort_clu_index[ int(total_entry * (double(MBin_NClus_cut[i]) /100.)) ] ] / double(NClus_hist -> GetBinWidth(1)) ) + 1 ));
0077     }
0078
0079     tree_out -> Fill();
0080
0081     c1 -> cd();
0082     c1 -> Print(Form("%s/%s", output_directory.c_str(), "NClus_hist.pdf"));
0083
0084     file_out -> cd();
0085     tree_out -> SetDirectory(file_out);
0086     tree_out -> Write("", TObject::kOverwrite);
0087
0088     file_out -> Close();
0089
0090 }