Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:12:26 

0001 #include "../INTTDSTchain.C"
0002 vector<string> read_list (string MC_directory, string MC_list_name)
0003 {
0004     string list_buffer;
0005     ifstream data_list;
0006     data_list.open((MC_directory + "/" + MC_list_name).c_str());
0007 
0008     vector<string> output_vec; output_vec.clear();
0009 
0010     while (1)
0011     {
0012         data_list >> list_buffer;
0013         if (!data_list.good())
0014         {
0015             break;
0016         }
0017         output_vec.push_back(list_buffer);
0018     }
0019     cout<<"size : "<<output_vec.size()<<endl;
0020 
0021     return output_vec;
0022 }
0023 
0024 void sanity_check()
0025 {
0026     string folder_direction = "/sphenix/user/ChengWei/sPH_dNdeta/HIJING_ana398_xvtx-0p04cm_yvtx0p24cm_zvtx-20cm_dummyAlignParams";
0027     string MC_list_name = "file_list.txt";
0028     vector <string> file_list = read_list(folder_direction,MC_list_name);
0029 
0030     TChain * chain_in = new TChain("EventTree");
0031     INTTDSTchain inttDSTMC(chain_in,folder_direction,file_list);
0032     std::printf("inttDSTMC N event : %lli \n", chain_in->GetEntries());
0033 
0034     TH2F * MC_vtx_xy = new TH2F("MC_vtx_xy","",100,-0.2,0.2,100,-0.,0.4);
0035     MC_vtx_xy -> GetXaxis() -> SetTitle("X axis [cm]");
0036     MC_vtx_xy -> GetYaxis() -> SetTitle("Y axis [cm]");
0037 
0038     TH1F * MC_vtx_z = new TH1F("MC_vtx_z","",100,-35,35);
0039     MC_vtx_z -> GetXaxis() -> SetTitle("Z axis [cm]");
0040     MC_vtx_z -> GetYaxis() -> SetTitle("Entry");
0041 
0042     TH1F * MC_vtx_Npart = new TH1F("Npart","",100,0,10000);
0043     MC_vtx_Npart -> GetXaxis() -> SetTitle("Npart");
0044     MC_vtx_Npart -> GetYaxis() -> SetTitle("Entry");
0045 
0046     TH2F * MC_Npart_Ntrack = new TH2F("MC_Npart_Ntrack","",200,0,10000, 200, 0, 4000);
0047     MC_Npart_Ntrack -> GetXaxis() -> SetTitle("Npart");
0048     MC_Npart_Ntrack -> GetYaxis() -> SetTitle("Ntrack");
0049 
0050     TH1F * MC_INTT_multiplicity = new TH1F("MC_INTT_multiplicity","",100,0,8000);
0051     MC_INTT_multiplicity -> GetXaxis() -> SetTitle("INTT N cluster");
0052     MC_INTT_multiplicity -> GetYaxis() -> SetTitle("Entry");
0053 
0054     TH1F * MC_INTT_CluPhiSize = new TH1F("MC_INTT_CluPhiSize","",25,0,25);
0055     MC_INTT_CluPhiSize -> GetXaxis() -> SetTitle("Size");
0056     MC_INTT_CluPhiSize -> GetYaxis() -> SetTitle("Entry");
0057 
0058     TH1F * MC_INTT_CluZSize = new TH1F("MC_INTT_CluZSize","",25,0,25);
0059     MC_INTT_CluZSize -> GetXaxis() -> SetTitle("Size");
0060     MC_INTT_CluZSize -> GetYaxis() -> SetTitle("Entry");
0061 
0062     TH1F * MC_INTT_CluADC = new TH1F("MC_INTT_CluADC","",150,0,3000);
0063     MC_INTT_CluADC -> GetXaxis() -> SetTitle("ADC");
0064     MC_INTT_CluADC -> GetYaxis() -> SetTitle("Entry");
0065 
0066     TH2F * MC_INTT_clu_nhit_corr = new TH2F("MC_INTT_clu_nhit_corr","",200,0,45000,200,0,9000);
0067     MC_INTT_clu_nhit_corr -> GetXaxis() -> SetTitle("NTrkrhits");
0068     MC_INTT_clu_nhit_corr -> GetYaxis() -> SetTitle("NClu");
0069 
0070     TH2F * MC_Ntrack_NClu = new TH2F("MC_Ntrack_NClu","",200,0,9000, 200, 0, 4000);
0071     MC_Ntrack_NClu -> GetXaxis() -> SetTitle("NClus");
0072     MC_Ntrack_NClu -> GetYaxis() -> SetTitle("Ntrack");
0073 
0074     TH1F * MC_Ntrack_1D = new TH1F("MC_Ntrack_1D","",200,0,4000);
0075     MC_Ntrack_1D -> GetXaxis() -> SetTitle("Ntrack");
0076     MC_Ntrack_1D -> GetYaxis() -> SetTitle("Entry");
0077     
0078     TH1F * MC_true_track_eta = new TH1F("MC_true_track_eta","",29,-2.9,2.9);
0079     MC_true_track_eta -> GetXaxis() -> SetTitle("#eta");
0080     MC_true_track_eta -> GetYaxis() -> SetTitle("dN/d#eta");
0081 
0082     TH2F * MC_NClu_Zvtx = new TH2F("MC_NClu_Zvtx","",200,0,9000, 200, -500, 500);
0083     MC_NClu_Zvtx -> GetXaxis() -> SetTitle("NClus");
0084     MC_NClu_Zvtx -> GetYaxis() -> SetTitle("Z vertex [mm]");
0085 
0086     TH1F * MC_inner_phi_3 = new TH1F("MC_inner_phi_3","",100,0,360);
0087     MC_inner_phi_3 -> GetXaxis() -> SetTitle("Inner phi");
0088     MC_inner_phi_3 -> GetYaxis() -> SetTitle("Entry");
0089 
0090     TH1F * MC_inner_phi_4 = new TH1F("MC_inner_phi_4","",100,0,360);
0091     MC_inner_phi_4 -> GetXaxis() -> SetTitle("Inner phi");
0092     MC_inner_phi_4 -> GetYaxis() -> SetTitle("Entry");
0093 
0094     TH1F * MC_outer_phi_5 = new TH1F("MC_outer_phi_5","",100,0,360);
0095     MC_outer_phi_5 -> GetXaxis() -> SetTitle("Outer phi");
0096     MC_outer_phi_5 -> GetYaxis() -> SetTitle("Entry");
0097 
0098     TH1F * MC_outer_phi_6 = new TH1F("MC_outer_phi_6","",100,0,360);
0099     MC_outer_phi_6 -> GetXaxis() -> SetTitle("Outer phi");
0100     MC_outer_phi_6 -> GetYaxis() -> SetTitle("Entry");
0101 
0102     int good_evt_counter = 0;
0103     int eta_one_track_counter = 0;
0104     int eta_one_track_counter_evt = 0;
0105     int eta_two_track_counter_evt = 0;
0106 
0107     for (int i = 0; i < chain_in->GetEntries(); i++) {
0108         inttDSTMC.LoadTree(i);
0109         inttDSTMC.GetEntry(i);
0110 
0111         if (i % 1000 == 0){
0112             cout<<"evnet running : "<<i<<endl;
0113             cout<<"----------------------------------------------- -----------------------------------------------"<<endl;
0114             printf(" event : %i, NTruthVtx : %i, NClus : %i, NTrkrhits : %i \n", inttDSTMC.event, inttDSTMC.NTruthVtx, inttDSTMC.NClus, inttDSTMC.NTrkrhits);
0115             printf(" size of ClusX : %zu \n", inttDSTMC.ClusX -> size());
0116             printf("triggered VTX ? : %.2f, %.2f, %.2f \n", inttDSTMC.TruthPV_trig_x, inttDSTMC.TruthPV_trig_y, inttDSTMC.TruthPV_trig_z);
0117         }
0118 
0119         // MC_Ntrack_1D   -> Fill(inttDSTMC.UniqueAncG4P_TrackID -> size());
0120 
0121         if (inttDSTMC.NTruthVtx == 1){
0122 
0123             // for (int track_i = 0; track_i < inttDSTMC.UniqueAncG4P_TrackID -> size(); track_i++){
0124             //     if (fabs(inttDSTMC.UniqueAncG4P_Eta -> at(track_i)) <= 1.){
0125             //         eta_one_track_counter_evt += 1;
0126             //     }
0127 
0128             //     if (fabs(inttDSTMC.UniqueAncG4P_Eta -> at(track_i)) <= 2.){
0129             //         eta_two_track_counter_evt += 1;
0130             //     }
0131 
0132             // }
0133             MC_NClu_Zvtx -> Fill(inttDSTMC.NClus, (inttDSTMC.TruthPV_trig_z) * 10.);
0134             // MC_Ntrack_1D   -> Fill(inttDSTMC.UniqueAncG4P_TrackID -> size());
0135             MC_Ntrack_NClu -> Fill(inttDSTMC.NClus, eta_two_track_counter_evt);
0136             eta_one_track_counter_evt = 0;
0137             eta_two_track_counter_evt = 0;
0138             
0139             // if (0/*good_evt_counter*/ == 0) {
0140             //     if (inttDSTMC.NClus > 4813){
0141             //         good_evt_counter += 1;
0142             //         // cout<<"test : "<<inttDSTMC.UniqueAncG4P_TrackID -> size()<<" "<<inttDSTMC.NClus<<endl;
0143             //         for (int track_i = 0; track_i < inttDSTMC.UniqueAncG4P_TrackID -> size(); track_i++){
0144             //             MC_true_track_eta -> Fill(inttDSTMC.UniqueAncG4P_Eta -> at(track_i));
0145             //             if (fabs(inttDSTMC.UniqueAncG4P_Eta -> at(track_i)) <= 1.){
0146             //                 eta_one_track_counter += 1;
0147             //             }
0148             //         }
0149             //     }
0150             // }
0151             
0152             
0153             // cout<<i<<" z vertex : "<<inttDSTMC.TruthPV_trig_z<<endl;
0154             MC_vtx_xy -> Fill(inttDSTMC.TruthPV_trig_x, inttDSTMC.TruthPV_trig_y);
0155             MC_vtx_z  -> Fill( inttDSTMC.TruthPV_trig_z );
0156             MC_vtx_Npart -> Fill( inttDSTMC.TruthPV_trig_Npart );
0157             // MC_Npart_Ntrack -> Fill(inttDSTMC.TruthPV_trig_Npart, inttDSTMC.UniqueAncG4P_TrackID -> size());
0158             
0159             
0160             MC_INTT_multiplicity -> Fill( inttDSTMC.ClusX -> size() );
0161             MC_INTT_clu_nhit_corr -> Fill(inttDSTMC.NTrkrhits, inttDSTMC.ClusX -> size());
0162             
0163 
0164             for (int clu_i = 0; clu_i < inttDSTMC.ClusX -> size(); clu_i++){
0165                 MC_INTT_CluPhiSize -> Fill( inttDSTMC.ClusPhiSize -> at(clu_i) );
0166                 MC_INTT_CluZSize -> Fill( inttDSTMC.ClusZSize -> at(clu_i) );
0167                 MC_INTT_CluADC -> Fill( inttDSTMC.ClusAdc -> at(clu_i) );
0168                 
0169                 if (inttDSTMC.ClusLayer -> at(clu_i) == 3){
0170                     MC_inner_phi_3 -> Fill( (inttDSTMC.ClusY -> at(clu_i) < 0) ? atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) + 360 : atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) );
0171                 }
0172                 else if (inttDSTMC.ClusLayer -> at(clu_i) == 4){
0173                     MC_inner_phi_4 -> Fill( (inttDSTMC.ClusY -> at(clu_i) < 0) ? atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) + 360 : atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) );
0174                 }
0175                 else if (inttDSTMC.ClusLayer -> at(clu_i) == 5){
0176                     MC_outer_phi_5 -> Fill( (inttDSTMC.ClusY -> at(clu_i) < 0) ? atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) + 360 : atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) );
0177                 }
0178                 else if (inttDSTMC.ClusLayer -> at(clu_i) == 6){
0179                     MC_outer_phi_6 -> Fill( (inttDSTMC.ClusY -> at(clu_i) < 0) ? atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) + 360 : atan2(inttDSTMC.ClusY -> at(clu_i), inttDSTMC.ClusX -> at(clu_i) ) * (180./TMath::Pi()) );
0180                 }
0181                 
0182                
0183             }
0184 
0185 
0186         
0187         }
0188 
0189         // for (int vtx_i = 0; vtx_i < inttDSTMC.TruthPV_x -> size(); vtx_i++){
0190         //     printf("VTX vec ID-%i : %.2f, %.2f, %.2f \n", vtx_i, inttDSTMC.TruthPV_x -> at(vtx_i), inttDSTMC.TruthPV_y -> at(vtx_i), inttDSTMC.TruthPV_z -> at(vtx_i));
0191         // }
0192         // cout<<"~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ "<<endl;
0193         // for (int clu_i = 0; clu_i < ClusX -> size(); clu_i++){
0194         //     printf("cluster %zu : X %.2f, Y %.2f, Z %.2f, sizePhi : %f, sizeZ : %f,  \n", clu_i, ClusX -> at(clu_i), ClusY -> at(clu_i), ClusZ -> at(clu_i), ClusPhiSize -> at(clu_i), ClusZSize -> at(clu_i));
0195         // }
0196     }
0197 
0198     MC_true_track_eta -> Scale(1./double(good_evt_counter));
0199     cout<<"N track in eta +- 1 : "<<eta_one_track_counter<<" good_evt_counter "<<good_evt_counter<<endl;
0200     // MC_true_track_eta -> Scale(1/0.2);
0201 
0202     TFile * INTT_MC_sanity_out = new TFile(Form("%s/MC_sanity_intt.root",folder_direction.c_str()), "recreate");
0203     MC_vtx_xy -> Write();
0204     MC_vtx_z -> Write();
0205     MC_NClu_Zvtx -> Write();
0206     MC_vtx_Npart -> Write();
0207     MC_Npart_Ntrack -> Write();
0208     MC_INTT_multiplicity -> Write();
0209     MC_INTT_CluPhiSize -> Write();
0210     MC_INTT_CluZSize -> Write();
0211     MC_INTT_CluADC -> Write();
0212     MC_INTT_clu_nhit_corr -> Write();
0213     MC_Ntrack_NClu -> Write();
0214     MC_Ntrack_1D -> Write();
0215     MC_true_track_eta -> Write();
0216     MC_inner_phi_3 -> Write();
0217     MC_inner_phi_4 -> Write();
0218     MC_outer_phi_5 -> Write();
0219     MC_outer_phi_6 -> Write();
0220     INTT_MC_sanity_out -> Close();
0221 
0222 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!