analysis_INTT_CW/old_code/truth_eta_dist.C

0001 // #include "../DST_MC/INTTReadTree.h"
0002 #include "../../INTTDSTchain.C"
0003
0004 void Characterize_Pad (TPad *pad, float left = 0.15, float right = 0.1, float top = 0.1, float bottom = 0.12, bool set_logY = false, int setgrid_bool = 0)
0005 {
0006     if (setgrid_bool == true) {pad -> SetGrid (1, 1);}
0007     pad -> SetLeftMargin   (left);
0008     pad -> SetRightMargin  (right);
0009     pad -> SetTopMargin    (top);
0010     pad -> SetBottomMargin (bottom);
0011     pad -> SetTicks(1,1);
0012     if (set_logY == true)
0013     {
0014         pad -> SetLogy (1);
0015     }
0016
0017 }
0018
0019 void truth_eta_dist()
0020 {
0021     string input_directory = "/sphenix/user/ChengWei/sPH_dNdeta/HIJING_ana398_xvtx-0p04cm_yvtx0p24cm_zvtx-20cm_dummyAlignParams";
0022     string file_name = "MC_ZF_xyvtx";
0023     string out_folder_directory = input_directory + "/Eta_test_" + file_name;
0024     string MC_list_name = "file_list.txt";
0025     string tree_name = "EventTree";
0026
0027     double INTT_layer_R = 7.188; // note: the innermost, B0L0
0028     double runEvent = 80000;
0029
0030     // INTTReadTree * INTTClu = new INTTReadTree(data_type, input_directory, MC_list_name, tree_name, clu_size_cut, clu_sum_adc_cut);
0031
0032     TChain * chain_in = new TChain(tree_name.c_str());
0033     INTTDSTchain * inttDSTMC = new INTTDSTchain(chain_in, input_directory, MC_list_name);
0034     long long N_event = chain_in->GetEntries();
0035     std::printf("inttDSTMC N event : %lli \n", N_event);
0036
0037     // cout<<"Total event : "<<INTTClu -> GetNEvt()<<endl;
0038
0039     TCanvas * c3 = new TCanvas("c3","c3",900,800); c3 -> cd();
0040     // TPad *pad_obj = new TPad(Form("pad_obj"), "", 0.0, 0.0, 1.0, 1.0);
0041     // Characterize_Pad(pad_obj, 0.18,  0.1,  0.1,  0.12, 0, 0);
0042     // pad_obj -> SetTicks(1,1);
0043     // pad_obj -> Draw();
0044     // pad_obj -> cd();
0045
0046     TH1F * true_track_Eta_1D = new TH1F("","",145,-2.9,2.9);
0047     true_track_Eta_1D -> SetMarkerStyle(20);
0048     true_track_Eta_1D -> SetMarkerSize(0.8);
0049     true_track_Eta_1D -> SetMarkerColor(TColor::GetColor("#1A3947"));
0050     true_track_Eta_1D -> SetLineColor(TColor::GetColor("#1A3947"));
0051     true_track_Eta_1D -> SetLineWidth(2);
0052     true_track_Eta_1D -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0053     true_track_Eta_1D -> GetXaxis() -> SetTitle("#eta");
0054     true_track_Eta_1D -> GetYaxis() -> SetRangeUser(0,50);
0055     true_track_Eta_1D -> SetTitleSize(0.06, "X");
0056     true_track_Eta_1D -> SetTitleSize(0.06, "Y");
0057     true_track_Eta_1D -> GetXaxis() -> SetTitleOffset (0.71);
0058     true_track_Eta_1D -> GetYaxis() -> SetTitleOffset (1.1);
0059     true_track_Eta_1D -> GetXaxis() -> CenterTitle(true);
0060     true_track_Eta_1D -> GetYaxis() -> CenterTitle(true);
0061
0062     TH2F * true_track_EtaPhi_2D = new TH2F("","",1450,-2.9,2.9, 300, -4, 4);
0063     true_track_EtaPhi_2D -> GetYaxis() -> SetTitle("#phi");
0064     true_track_EtaPhi_2D -> GetXaxis() -> SetTitle("#eta");
0065     // true_track_EtaPhi_2D -> GetYaxis() -> SetRangeUser(0,50);
0066     // true_track_EtaPhi_2D -> SetTitleSize(0.06, "X");
0067     // true_track_EtaPhi_2D -> SetTitleSize(0.06, "Y");
0068     // true_track_EtaPhi_2D -> GetXaxis() -> SetTitleOffset (0.71);
0069     // true_track_EtaPhi_2D -> GetYaxis() -> SetTitleOffset (1.1);
0070     // true_track_EtaPhi_2D -> GetXaxis() -> CenterTitle(true);
0071     // true_track_EtaPhi_2D -> GetYaxis() -> CenterTitle(true);
0072
0073     TH2F * true_track_EtaZ_2D = new TH2F("","",1450,-2.9,2.9, 300, -23, -17);
0074     true_track_EtaZ_2D -> GetYaxis() -> SetTitle("Z");
0075     true_track_EtaZ_2D -> GetXaxis() -> SetTitle("#eta");
0076
0077     TH1F * true_track_PID = new TH1F("","",20000,-10000,10000);
0078     true_track_PID -> GetYaxis() -> SetTitle("PID");
0079     true_track_PID -> GetXaxis() -> SetTitle("Entry");
0080
0081
0082     cout<<"true_track_Eta_1D, binwidth : "<<true_track_Eta_1D->GetBinWidth(1)<<endl;
0083     double good_event_counting = 0;
0084
0085     double z_range_l = -23;
0086     double z_range_r = -17;
0087
0088     std::unordered_map<int, int> PID_countMap;
0089
0090     for (int event_i = 0; event_i < 80000; event_i ++)
0091     {
0092         inttDSTMC->LoadTree(event_i);
0093         inttDSTMC->GetEntry(event_i);
0094
0095         double z = inttDSTMC -> TruthPV_trig_z;
0096
0097         if (inttDSTMC -> NTruthVtx != 1) {continue;}
0098         // if (event_i % 100 == 0){cout<<inttDSTMC -> centrality_bimp<<endl;}
0099         // if (inttDSTMC -> centrality_bimp != 5) {continue;}
0100         // if (z < z_range_l || z > z_range_r) {continue;}
0101
0102         double INTT_eta_acceptance_l = -0.5 * TMath::Log((sqrt(pow(-23.-z_range_l,2)+pow(INTT_layer_R,2))-(-23.-z_range_l)) / (sqrt(pow(-23.-z_range_l,2)+pow(INTT_layer_R,2))+(-23.-z_range_l))); // note : left
0103         double INTT_eta_acceptance_r =  -0.5 * TMath::Log((sqrt(pow(23.-z_range_r,2)+pow(INTT_layer_R,2))-(23.-z_range_r)) / (sqrt(pow(23.-z_range_r,2)+pow(INTT_layer_R,2))+(23.-z_range_r))); // note : right
0104
0105
0106         good_event_counting += 1;
0107
0108         // if (event_i % 100 == 0){cout<<"z : "<<z<<" eta : "<<INTT_eta_acceptance_l<<" "<<INTT_eta_acceptance_r<<" N tracks : "<<inttDSTMC->UniqueAncG4P_Eta->size()<<" NClus : "<<inttDSTMC->NClus<<endl;}
0109
0110         for (int track_i = 0; track_i < inttDSTMC->UniqueAncG4P_Eta->size(); track_i++)
0111         {
0112             if (inttDSTMC->UniqueAncG4P_Eta->at(track_i) > INTT_eta_acceptance_l && inttDSTMC->UniqueAncG4P_Eta->at(track_i) < INTT_eta_acceptance_r)
0113             {
0114                 true_track_Eta_1D -> Fill(inttDSTMC->UniqueAncG4P_Eta->at(track_i));
0115                 // double fillID = true_track_PID -> Fill(inttDSTMC->UniqueAncG4P_PID->at(track_i));
0116                 // if (fillID == -1) {cout<<"PID out of the set range : "<<inttDSTMC->UniqueAncG4P_PID->at(track_i)<<endl;}
0117                 int value = inttDSTMC->UniqueAncG4P_PID->at(track_i);
0118                 PID_countMap[value]++;
0119
0120                 // true_track_EtaZ_2D -> Fill(inttDSTMC->UniqueAncG4P_Eta->at(track_i),z);
0121             }
0122
0123             // true_track_Eta_1D -> Fill(inttDSTMC->UniqueAncG4P_Eta->at(track_i));
0124             // true_track_EtaPhi_2D -> Fill(inttDSTMC->UniqueAncG4P_Eta->at(track_i),inttDSTMC->UniqueAncG4P_Phi->at(track_i));
0125
0126             true_track_EtaZ_2D -> Fill(inttDSTMC->UniqueAncG4P_Eta->at(track_i),z);
0127         }
0128
0129         cout<<"z : "<<z<<" eta : "<<INTT_eta_acceptance_l<<" "<<INTT_eta_acceptance_r<<" N tracks : "<<inttDSTMC->UniqueAncG4P_Eta->size()<<" NClus : "<<inttDSTMC->NClus<<endl;
0130         // break;
0131
0132     }
0133
0134     // true_track_Eta_1D -> Scale(1./double(runEvent * true_track_Eta_1D -> GetBinWidth(1) ));
0135     // true_track_Eta_1D -> Scale(1./double(runEvent));
0136     // true_track_Eta_1D -> GetYaxis() -> SetRangeUser(0,800);
0137
0138     cout<<" good_event_counting : "<<good_event_counting<<endl;
0139
0140     true_track_Eta_1D -> Scale(1./double(true_track_Eta_1D -> GetBinWidth(1) ));
0141     true_track_Eta_1D -> Scale(1./double(good_event_counting));
0142     true_track_Eta_1D -> GetYaxis() -> SetRangeUser(0,1000);
0143
0144     true_track_Eta_1D -> Draw("hist");
0145     c3->Print(Form("%s/true_track_Eta_1D.pdf",out_folder_directory.c_str()));
0146
0147     // pad_obj -> Clear();
0148
0149     true_track_EtaPhi_2D -> Draw("colz0");
0150     c3->Print(Form("%s/true_track_EtaPhi_2D.pdf",out_folder_directory.c_str()));
0151
0152     true_track_EtaZ_2D -> Draw("colz0");
0153     c3->Print(Form("%s/true_track_EtaZ_2D.pdf",out_folder_directory.c_str()));
0154     c3 -> Clear();
0155
0156     std::cout << "PID   Count" << std::endl;
0157     for (const auto& pair : PID_countMap) {
0158         std::cout << pair.first << "       " << pair.second << std::endl;
0159     }
0160
0161     // c3 -> SetLogy(1);
0162     // true_track_PID -> Draw("hist");
0163     // c3->Print(Form("%s/true_track_PID.pdf",out_folder_directory.c_str()));
0164     // c3 -> Clear();
0165
0166 }