diff --git a/macro/mcbm/mcbm_transport.C b/macro/mcbm/mcbm_transport.C index 0f156dfedec9554370e42ab67babfb3189d4d6d3..1ec2d36c1c8a6264eb1f0884dd3fc1a096bac5a4 100644 --- a/macro/mcbm/mcbm_transport.C +++ b/macro/mcbm/mcbm_transport.C @@ -15,9 +15,9 @@ void SetTrack(CbmTransport*, Double_t, Int_t, Double_t, Double_t, Double_t); void mcbm_transport(Int_t nEvents = 10, - const char* setupName = "mcbm_beam_2019_11", + const char* setupName = "mcbm_beam_2020_03", +// const char* setupName = "mcbm_beam_2019_11", // const char* setupName = "mcbm_beam_2019_03", -// const char* setupName = "mcbm_beam_2018_11", // const char* setupName = "sis18_mcbm_25deg_long", const char* output = "test", const char* inputFile = "") diff --git a/macro/psd/build_correlations_tof.C b/macro/psd/build_correlations_tof.C index a6d24d1a893e9442b9a370b7aeed6168e9cd3fd7..81b99c28c654747a20b8fc47165f61fa584694eb 100644 --- a/macro/psd/build_correlations_tof.C +++ b/macro/psd/build_correlations_tof.C @@ -2,12 +2,12 @@ #include <FairRuntimeDb.h> void build_correlations_tof( - const string& parFile = "/home/nikolay/FairRoot/cbmroot_trunk/macro/beamtime/mcbm2019/data/unp_mcbm_params_384.root", - const string& digiFile = "/home/nikolay/FairRoot/cbmroot_trunk/macro/beamtime/mcbm2019/data/unp_mcbm_384.root", - const string& outDir = "/home/nikolay/FairRoot/cbmroot_trunk/macro/beamtime/mcbm2019/data/", - const string& outFile = "reco_mcbm_dec19.root", - const unsigned int uRunId = 384, // used for the output folder - int nEvents = 0 + const string& digiFile = "/home/nikolay/FairRoot/cbmroot_trunk/macro/beamtime/mcbm2020/data/unp_mcbm_582.root", + const string& parFile = "/home/nikolay/FairRoot/cbmroot_trunk/macro/beamtime/mcbm2020/data/unp_mcbm_params_582.root", + const string& outDir = "/home/nikolay/FairRoot/cbmroot_trunk/macro/psd/Results/", + const string& outFile = "reco_mcbm_mar20.root", + const unsigned int uRunId = 582, // used for the output folder + int nEvents = 100 ) { @@ -26,15 +26,16 @@ void build_correlations_tof( TList *parFileList = new TList(); const Double_t eb_fixedTimeWindow {200.}; - const Int_t eb_TriggerMinNumberT0 {1}; + const Int_t eb_TriggerMinNumberT0 {0}; const Int_t eb_TriggerMinNumberSts {0}; const Int_t eb_TriggerMinNumberMuch{0}; - const Int_t eb_TriggerMinNumberTof {10}; + const Int_t eb_TriggerMinNumberTof {0}; const Int_t eb_TriggerMinNumberPsd {1}; //ToF Settings Int_t calMode = 93; Int_t calSel = 1; + Int_t calSelRead = 0; Int_t calSm = 40; Int_t RefSel = 1; TString cFileId = "385.100.-192.0"; @@ -78,9 +79,11 @@ void build_correlations_tof( fRun->SetEventHeaderPersistence(kFALSE); FairFileSource* inputSource = new FairFileSource(digiFile); fRun->SetSource(inputSource); - TString runId = TString::Format("%03u_", uRunId); - TString outFileName = outDir + "/events_" + runId + outFile; + TString runId = TString::Format("%03u", uRunId); + TString QaDir = outDir + "result_run" + runId; + TString outFileName = QaDir + "/events_" + runId + "_" + outFile; //remove(outFileName.c_str()); + gSystem->Exec(Form("mkdir -p %s", QaDir.Data())); FairRootFileSink* outputSink = new FairRootFileSink(outFileName); fRun->SetSink(outputSink); //fRun->SetGenerateRunInfo(kTRUE); @@ -139,9 +142,6 @@ void build_correlations_tof( //tofClust->SetDeadStrips(15,23); // declare dead strip for T0M3,Rpc0,Strip 23 - //Int_t calSelRead = calSel; //TODO CHECK!!!! - Int_t calSelRead = 0; - //if (calSel<0) calSelRead=0; TString cCalibFname=Form("/%s_set%09d_%02d_%01dtofClust.hst.root",cFileId.Data(),iCalSet,calMode,calSelRead); tofClust->SetCalParFileName(TofFileFolder + cCalibFname); TString cOutFname=Form("tofClust_%s_set%09d.hst.root",cFileId.Data(),iCalSet); @@ -346,7 +346,7 @@ switch (calMode) { CbmPsdMCbmQaReal* qaTask = new CbmPsdMCbmQaReal(); - qaTask->SetOutputDir(Form("result_run%d", uRunId)); + qaTask->SetOutputDir(QaDir.Data()); fRun->AddTask(qaTask); diff --git a/macro/psd/standalone_event_build.C b/macro/psd/standalone_event_build.C index 9344603a826ef8df56ed96e1608d2c4018c9ab32..9377b5aa866656494514757b46ebc2ce505d7631 100644 --- a/macro/psd/standalone_event_build.C +++ b/macro/psd/standalone_event_build.C @@ -13,15 +13,13 @@ void standalone_event_build(TString inFileName) TTree* inTree = (TTree*)inFile->Get("cbmsim"); if(inTree == nullptr) { cout<<"tree not found"<<endl; return; } - std::vector<CbmPsdDigi> *fPsdDigis = new std::vector<CbmPsdDigi>(); + std::vector<CbmPsdDigi> *fPsdDigis = new std::vector<CbmPsdDigi>; inTree->SetBranchAddress("PsdDigi", &fPsdDigis); + if(fPsdDigis == nullptr) { cout<<"no psd digis vector in tree"<<endl; return; } - const Float_t event_time_gate = 100.; //in ns - const Int_t total_channels = 16; - const Int_t total_cells = 10; - UInt_t ch_address[total_channels]; - for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) - ch_address[ch_iter] = 8 + ch_iter*1024; + const Float_t event_time_gate = 200.; //in ns + const Int_t total_channels = 9; + const Int_t total_cells = 9; Double_t ch_energy[total_channels]; for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) @@ -31,50 +29,56 @@ void standalone_event_build(TString inFileName) for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) ch_energy_overEv[ch_iter] = 0.; Int_t Ev_counter = 0; - Int_t Fired_channels_in_ev = 0; + Int_t Fired_psd_channels_in_ev = 0; TH1 *th1_hist_ptr = nullptr; TH2 *th2_hist_ptr = nullptr; - th1_hist_ptr = new TH1F("Energy_profile", "PSD energy profile per event; cell # []; Energy [adc counts]", 16, 0, 16); - th1_hist_ptr = new TH1F("Energy_distrib", "PSD energy distribution; Energy over event all channels [adc counts]; Counts []", 1000, 0, 100000); - th1_hist_ptr = new TH1F("Spill_structure", "PSD time spill structure; Time [s]; Counts []", 5000, 0, 5000); - th1_hist_ptr = new TH1F("Event_time_diff", "PSD event time difference; Time [ns]; Counts []", (Int_t)300, 0, 300); - th1_hist_ptr = new TH1F("Fired_chs", "Fired channels in event; fired channels []; Counts []", 16, 0, 16); + th1_hist_ptr = new TH1F("Energy_profile", "PSD energy profile per event; cell # []; Energy [adc counts]", total_channels, 0, total_channels); + th1_hist_ptr = new TH1F("Energy_distrib", "PSD energy distribution; Energy over event all channels [adc counts]; Counts []", 1000, 0, 10000); + th1_hist_ptr = new TH1F("Spill_structure", "PSD time spill structure; Time [s]; Counts []", 600, 20, 40); + + th1_hist_ptr = new TH1F("Digi_time_diff", "PSD digi time difference; Time [ns]; Counts []", (Int_t)500, 0, 500); + th1_hist_ptr = new TH1F("Fired_psd_chs", "Fired PSD channels in event; fired channels []; Counts []", total_cells, 0, total_cells); + th1_hist_ptr = new TH1F("Cosmic_check", "Half a sum from two hodo channels; Energy []; Counts []", 100, 0, 100); + + th2_hist_ptr = new TH2F("PsdVsHodo", "Energy in Psd vs energy in Hodo; Psd energy [adc counts]; Hodo energy [adc counts]", 1000, 0, 2000, 140, 0, 70); for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) - th1_hist_ptr = new TH1F( Form("Energy_distrib_ch%i", ch_iter), Form("Energy_distrib_ch%i; Energy [adc counts]; Counts []", ch_iter), 1000, 0, 10000); + th1_hist_ptr = new TH1F( Form("Energy_distrib_ch%i", ch_iter), Form("Energy_distrib_ch%i; Energy [adc counts]; Counts []", ch_iter), 100, 0, 2000); + for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) + th1_hist_ptr = new TH1F( Form("Energy_distrib_good_ev_ch%i", ch_iter), Form("Energy_distrib_good_ev_ch%i; Energy [adc counts]; Counts []", ch_iter), 100, 0, 2000); + + Double_t event_time = -1.; + Double_t prev_event_time = -1.; + Bool_t IsFirstEvent = true; for(UInt_t ev = 0; ev < inTree->GetEntries(); ev++) - //for(UInt_t ev = 0; ev < 901; ev++) { inTree->GetEntry(ev); Int_t nPsdDigis = fPsdDigis->size(); - if (ev % 1000 == 0 ) std::cout << "Entry " << ev << " / " << inTree->GetEntries() << " Psd Digis: " << nPsdDigis << std::endl; - - //if(nPsdDigis > 100) inTree->Show(); + if (ev % 10 == 0 ) std::cout << "Entry " << ev << " / " << inTree->GetEntries() << " Psd Digis: " << nPsdDigis << std::endl; + if(nPsdDigis == 0) continue; CbmPsdDigi digi; - Double_t event_time = 0.; - Double_t prev_event_time = 0.; Bool_t IsNewEvent = false; for (Int_t iDigi = 0; iDigi < nPsdDigis; iDigi++) { digi = (*fPsdDigis)[iDigi]; - - //if(! (digi.GetSectionID() == 0 && digi.GetTime()>4e12) ) continue; - if(! ( digi.GetTime()>4e12) ) continue; + //if(! (digi.GetSectionID() == 0 ) ) continue; + //if(digi.GetSectionID() > 8) continue; if( digi.GetTime()-event_time > event_time_gate ) { - prev_event_time = event_time; + prev_event_time = event_time; event_time = digi.GetTime(); - if( prev_event_time > 1. ) IsNewEvent = true; + IsNewEvent = true; } else IsNewEvent = false; if(IsNewEvent) { + if(IsFirstEvent) { IsFirstEvent = false; continue; } //skip info before first event Bool_t GoodEvent = true; if(GoodEvent) @@ -82,28 +86,45 @@ void standalone_event_build(TString inFileName) th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Spill_structure" ))); th1_hist_ptr->Fill(event_time/1e9); - th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Event_time_diff" ))); + + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Digi_time_diff" ))); th1_hist_ptr->Fill(event_time - prev_event_time); - Double_t total_energy = 0.; - Fired_channels_in_ev = 0; + Double_t psd_energy = 0.; + Fired_psd_channels_in_ev = 0; - for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) + for(UInt_t ch_iter = 0; ch_iter < total_cells; ch_iter++) { ch_energy_overEv[ch_iter] += ch_energy[ch_iter]; - total_energy += ch_energy[ch_iter]; + psd_energy += ch_energy[ch_iter]; + + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( Form("Energy_distrib_good_ev_ch%i", ch_iter) ))); + if(ch_energy[ch_iter] > 1.) th1_hist_ptr->Fill(ch_energy[ch_iter]); - if(ch_energy[ch_iter] > 1.) Fired_channels_in_ev++; + if(ch_energy[ch_iter] > 1.) Fired_psd_channels_in_ev++; ch_energy[ch_iter] = 0.; } - th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Energy_distrib" ))); - th1_hist_ptr->Fill(total_energy); - th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Fired_chs" ))); - th1_hist_ptr->Fill(Fired_channels_in_ev); - - //cout<<worked_cells_counter<<endl; + //if(Fired_psd_channels_in_ev < 2) continue; + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Energy_distrib" ))); + th1_hist_ptr->Fill(psd_energy); + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Fired_psd_chs" ))); + th1_hist_ptr->Fill(Fired_psd_channels_in_ev); + +/* + Double_t hodo_energy = 0.5*(ch_energy[17] + ch_energy[19]); + ch_energy_overEv[17] += ch_energy[17]; + ch_energy_overEv[19] += ch_energy[19]; + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Cosmic_check" ))); + if(ch_energy[17] > 1. && ch_energy[19] > 1.) th1_hist_ptr->Fill(0.5*(ch_energy[17] + ch_energy[19])); + + ch_energy[17] = 0.; + ch_energy[19] = 0.; + + th2_hist_ptr = ((TH2*)(gDirectory->FindObjectAny( "PsdVsHodo" ))); + th2_hist_ptr->Fill(psd_energy, hodo_energy); +*/ Ev_counter++; }//is good event @@ -118,12 +139,12 @@ void standalone_event_build(TString inFileName) Double_t Edep = digi.GetEdep(); UInt_t SectionID = digi.GetSectionID(); + //if(ch_energy[SectionID] > 1.) cout<<"PileUp candidate in channel "<<SectionID<<endl; ch_energy[SectionID] = Edep; th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( Form("Energy_distrib_ch%i", SectionID) ))); th1_hist_ptr->Fill(Edep); - }//digi loop }//event loop @@ -137,7 +158,7 @@ void standalone_event_build(TString inFileName) th1_hist_ptr->Draw(); TCanvas *canvas_time_diff_ev = new TCanvas("c_time_diff_ev", "c_time_diff_ev"); - th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Event_time_diff" ))); + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Digi_time_diff" ))); th1_hist_ptr->Draw(); TCanvas *canvas_energy_pr = new TCanvas("c_energy_pr", "c_energy_pr"); @@ -149,7 +170,14 @@ void standalone_event_build(TString inFileName) th1_hist_ptr->Draw(); TCanvas *canvas_fired_chs = new TCanvas("c_fired_ch", "c_fired_ch"); - th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Fired_chs" ))); + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Fired_psd_chs" ))); + th1_hist_ptr->Draw(); + + TCanvas *canvas_psd_hodo = new TCanvas("c_psd_hodo", "c_psd_hodo"); + th2_hist_ptr = ((TH2*)(gDirectory->FindObjectAny( "PsdVsHodo" ))); + th2_hist_ptr->Draw("colz"); + + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Cosmic_check" ))); th1_hist_ptr->Draw(); TCanvas *canvas_energy_dstr_ch = new TCanvas("c_energy_dstr_ch", "c_energy_dstr_ch"); @@ -161,4 +189,13 @@ void standalone_event_build(TString inFileName) th1_hist_ptr->Draw(); } + TCanvas *canvas_energy_dstr_gEv_ch = new TCanvas("c_energy_dstr_gEv_ch", "c_energy_dstr_gEv_ch"); + canvas_energy_dstr_gEv_ch->DivideSquare(total_channels); + for(UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) + { + canvas_energy_dstr_gEv_ch->cd(ch_iter+1); + th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( Form("Energy_distrib_good_ev_ch%i", ch_iter) ))); + th1_hist_ptr->Draw(); + } + }