diff --git a/Projects/Vendeta/macro/FillTOFHisto.C b/Projects/Vendeta/macro/FillTOFHisto.C new file mode 100644 index 0000000000000000000000000000000000000000..70f18372e2370a19f5241a5090f025a5aa563f6d --- /dev/null +++ b/Projects/Vendeta/macro/FillTOFHisto.C @@ -0,0 +1,194 @@ +// c++ includes +#include <vector> +#include <iostream> +#include <iomanip> +#include <stdlib.h> + +// ROOT includes +#include "TChain.h" +#include "TFile.h" +#include "TObject.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TGraph.h" + +using namespace std; + +int NumberOfDetectors= 72; +int NumberOfAnodes= 11; +int nentries=1e6; +<<<<<<< HEAD + +///////////////////////////////////// +int main(int argc, char** argv){ + int run = atof(argv[1]); + + auto chain = new TChain("PhysicsTree"); + chain->Add(Form("/home/cyril/Workflow/nptool/Outputs/Analysis/run%i/run%i.root",run,run)); + + nentries = chain->GetEntries(); + /* nentries = 50000000; */ + cout << "Number of entries: " << nentries << endl; + + TFile* ofile = new TFile(Form("histos_ToF/histo_tof_file_run%i.root",run),"recreate"); + TH2F* hLG_LSci[11]; + TH2F* hHG_LSci[11]; + TH1F* hLG[791]; + TH1F* hHG[791]; + + + TGraph *gFlyPath = new TGraph(); + gFlyPath->SetTitle(" ; Vendeta-Anode Index ; fly length (mm);"); + +======= +int run=86; +///////////////////////////////////// +void LoadRootFile(){ + chain = new TChain("PhysicsTree"); + chain->Add(Form("/home/faster/nptool/Outputs/Analysis/run%i.root",run)); + //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_1.root"); + //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_2.root"); + //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_3.root"); + //chain->Add("/home/faster/nptool/Outputs/Analysis/test_sampler_qdc_cf_4.root"); +} + +///////////////////////////////////// +void FillTOFHisto(){ + + LoadRootFile(); + nentries = chain->GetEntries(); + cout << "Number of entries: " << nentries << endl; + + TFile* ofile = new TFile(Form("histo_tof_file_run%i.root",run),"recreate"); + TH1F* hLG[791]; + TH1F* hHG[791]; + +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + vector<double>* FC_Q1 = new vector<double>(); + vector<int>* FC_Anode_ID = new vector<int>(); + vector<bool>* FC_FakeFission = new vector<bool>(); + + vector<double>* LG_Tof = new vector<double>(); + vector<int>* LG_ID = new vector<int>(); + vector<double>* LG_Q1 = new vector<double>(); + vector<double>* LG_Q2 = new vector<double>(); +<<<<<<< HEAD + vector<double>* LG_FlyPath = new vector<double>(); +======= +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + + vector<double>* HG_Tof = new vector<double>(); + vector<int>* HG_ID = new vector<int>(); + vector<double>* HG_Q1 = new vector<double>(); + vector<double>* HG_Q2 = new vector<double>(); +<<<<<<< HEAD + vector<double>* HG_FlyPath = new vector<double>(); +======= + + TFissionChamberPhysics* FC = new TFissionChamberPhysics(); + chain->SetBranchAddress("FissionChamber",&FC); +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + + chain->SetBranchAddress("FC_Q1",&FC_Q1); + chain->SetBranchAddress("FC_Anode_ID",&FC_Anode_ID); + chain->SetBranchAddress("FC_FakeFission",&FC_FakeFission); + chain->SetBranchAddress("LG_Tof",&LG_Tof); + chain->SetBranchAddress("LG_ID",&LG_ID); + chain->SetBranchAddress("LG_Q1",&LG_Q1); + chain->SetBranchAddress("LG_Q2",&LG_Q2); +<<<<<<< HEAD + chain->SetBranchAddress("LG_FlyPath",&LG_FlyPath); +======= +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + chain->SetBranchAddress("HG_Tof",&HG_Tof); + chain->SetBranchAddress("HG_ID",&HG_ID); + chain->SetBranchAddress("HG_Q1",&HG_Q1); + chain->SetBranchAddress("HG_Q2",&HG_Q2); +<<<<<<< HEAD + chain->SetBranchAddress("HG_FlyPath",&HG_FlyPath); + + for(int j=0; j<NumberOfAnodes; j++){ + + TString histo_name = Form("hLG_LSci_Anode%i",j+1); + hLG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300); + + histo_name = Form("hHG_LSci_Anode%i",j+1); + hHG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300); + + for(int i=0; i<NumberOfDetectors; i++){ + int index = j + i*NumberOfAnodes; + histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1); + hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300); + + histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1); + hHG[index] = new TH1F(histo_name,histo_name,2000,-100,300); +======= + + for(int i=0; i<NumberOfDetectors; i++){ + for(int j=0; j<NumberOfAnodes; j++){ + int index = j + i*NumberOfAnodes; + TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1); + hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300); + + histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1); + hHG[index] = new TH1F(histo_name,histo_name,2500,0,500); +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + } + } + for(int i=0; i<nentries; i++){ + chain->GetEntry(i); + + if(i%100000==0){ +<<<<<<< HEAD + cout << setprecision(2) << "\033[34m\r Processing run "<< run <<" ..." << (double)i/nentries*100 << "\% done" << flush; +======= + cout << setprecision(2) << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush; +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + } + + if(FC_Anode_ID->size()>0){ + bool Fake = FC_FakeFission->at(0); + int anode = FC_Anode_ID->at(0); + + + int mysize = LG_Tof->size(); + for(int j=0; j<mysize; j++){ + // LG // + int index_LG = (anode-1) + (LG_ID->at(j)-1)*NumberOfAnodes; + double PSD = LG_Q2->at(j)/LG_Q1->at(j); + if(LG_ID->at(j)>0 && anode>0 && Fake==0 && LG_Q1->at(j)>2000){ +<<<<<<< HEAD + hLG[index_LG]->Fill(LG_Tof->at(j)); + hLG_LSci[anode-1]->Fill(LG_ID->at(j),LG_Tof->at(j)); + gFlyPath->SetPoint(index_LG,index_LG,LG_FlyPath->at(j)); +======= + hLG[index_LG]->Fill(LG_Tof->at(j)); +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + } + } + + mysize = HG_Tof->size(); + for(int j=0; j<mysize; j++){ + // HG // + int index_HG = (anode-1) + (HG_ID->at(j)-1)*NumberOfAnodes; + double PSD = HG_Q2->at(j)/HG_Q1->at(j); + if(HG_ID->at(j)>0 && anode>0 && Fake==0 && HG_ID->size()==1){ +<<<<<<< HEAD + hHG[index_HG]->Fill(HG_Tof->at(j)); + hHG_LSci[anode-1]->Fill(HG_ID->at(j),HG_Tof->at(j)); +======= + hHG[index_HG]->Fill(HG_Tof->at(j)); +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + } + } + } + } + +<<<<<<< HEAD + ofile->WriteObject(gFlyPath,"gFlyPath"); +======= +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + ofile->Write(); + ofile->Close(); + +} diff --git a/Projects/Vendeta/macro/FitTofGammaPeak.C b/Projects/Vendeta/macro/FitTofGammaPeak.C new file mode 100644 index 0000000000000000000000000000000000000000..e4074e33e4d5b349b7f76c39e1e3fcf63c7101d4 --- /dev/null +++ b/Projects/Vendeta/macro/FitTofGammaPeak.C @@ -0,0 +1,205 @@ +// c++ includes +#include <vector> +#include <iostream> +#include <iomanip> +#include <fstream> +#include <stdlib.h> + +// ROOT includes +#include "TChain.h" +#include "TCanvas.h" +#include "TFile.h" +#include "TObject.h" +#include "TSpectrum.h" +#include "TGraph.h" +#include "TF1.h" +#include "TH1F.h" +#include "TH2F.h" + +using namespace std; + +int NumberOfDetectors=72; +int NumberOfAnodes=11; +int m_NumberOfGammaPeak=1; +double PosGammaPeak = 3.2; +double c = 299.792458; // mm/ns + +double sigma_anode[11], sigma_anode_HG[11]; + +bool Finder(TH1F* h, Double_t *mean, Double_t *sigma); + +///////////////////////////////////////////////////// +int main(int argc, char** argv){ + + int run = atof(argv[1]); + TString Append = argv[2]; + + auto ifile = new TFile(Form("histos_ToF/histo_tof_file_run%d"+Append+".root",run)); + auto gFlyPath = (TGraph*)ifile->Get("gFlyPath"); + + TCanvas *cLG[11]; + TCanvas *cHG[11]; + + ofstream ofile; +<<<<<<< HEAD + ofile.open(Form("ToF_calibs/Vendeta_Time_run%d.cal",run)); + TFile* orootfile = new TFile(Form("histos_ToF/histo_tof_fitted_run%d"+Append+".root",run),"recreate"); + + for(int i=0; i<NumberOfAnodes; i++){ + sigma_anode[i]=0; + sigma_anode_HG[i]=0; +======= + ofile.open(Form("Vendeta_Time_run%i.cal",run)); + + TFile* orootfile = new TFile(Form("histo_tof_fitted_run%i.root",run),"recreate"); + + for(int i=0; i<NumberOfAnodes; i++){ + sigma_anode[i]=0; +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + } + + Double_t* mean = new Double_t[1]; + Double_t* sigma = new Double_t[1]; + TGraph* gSigma_LG = new TGraph(); + TGraph* gSigma_HG = new TGraph(); + + int countLG=0; + int countHG=0; + for(int i=0; i<NumberOfDetectors; i++){ + for(int j=0; j<NumberOfAnodes; j++){ + // LG // + TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1); + TH1F* h = (TH1F*) ifile->FindObjectAny(histo_name); + + //cLG[j]->cd(i+1); + h->Draw(); + + int index = j + i*11; + double FlyPath = gFlyPath->GetPointY(index); + PosGammaPeak = FlyPath/c ; + + mean[0] = 0; + sigma[0] = 0; + TString LG_token = Form("Vendeta_DET%i_LG_ANODE%i_TIMEOFFSET",i+1,j+1); + if(Finder(h, mean, sigma)){ + ofile << LG_token << " " << -mean[0]+PosGammaPeak << endl; + gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]); + + sigma_anode[j] += sigma[0]; + + countLG++; + h->Write(); + } + else{ + ofile << LG_token << " 0" << endl; + } + + // HG // + histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1); + h = (TH1F*) ifile->FindObjectAny(histo_name); + + //cHG[j]->cd(i+1); + h->Draw(); + + mean[0] = 0; + sigma[0] = 0; + TString HG_token = Form("Vendeta_DET%i_HG_ANODE%i_TIMEOFFSET",i+1,j+1); + if(Finder(h, mean, sigma)){ + ofile << HG_token << " " << -mean[0]+PosGammaPeak << endl; + gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]); +<<<<<<< HEAD + sigma_anode_HG[j] += sigma[0]; +======= +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + countHG++; + h->Write(); + } + else{ + ofile << HG_token << " 0" << endl; + } + } + } + + TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600); + c1->Divide(2,1); + + gSigma_LG->SetMarkerStyle(8); + gSigma_HG->SetMarkerStyle(8); + + c1->cd(1); + gSigma_LG->Draw(); + c1->cd(2); + gSigma_HG->Draw(); + + gSigma_LG->SetName("sigma_LG"); + gSigma_HG->SetName("sigma_HG"); + + gSigma_LG->Write(); + gSigma_HG->Write(); + + orootfile->Write(); + orootfile->Close(); +<<<<<<< HEAD +======= + + for(int i=0; i<NumberOfAnodes; i++){ + cout << "Anode= " << i+1 << " / " << sigma_anode[i]/NumberOfDetectors << endl; + } +>>>>>>> 1f5a4966c35bb110ec3e906e38ea506f92dad29d + + double totLG = 0, totHG=0; + for(int i=0; i<NumberOfAnodes; i++){ + cout << "Anode= " << i+1 << " | LG: " << sigma_anode[i]/NumberOfDetectors << " HG: " << sigma_anode_HG[i]/NumberOfDetectors << endl; + totLG += sigma_anode[i]; + totHG += sigma_anode_HG[i]; + } + totLG = totLG/NumberOfDetectors/11.; + totHG = totHG/NumberOfDetectors/11.; + cout << "<sigma> LG: "<< totLG << " HG: "<<totHG << endl; +} + + +///////////////////////////////////////////////////// +bool Finder(TH1F* h, Double_t *mean, Double_t *sigma){ + Double_t resolsig = 2; + Double_t resolsigTSpec = 1; + Double_t seuil = 0.3; + + TSpectrum* s = new TSpectrum(m_NumberOfGammaPeak,resolsigTSpec); + + Int_t nfound = s->Search(h,resolsig,"new",seuil); + Double_t *xpeaks = s->GetPositionX(); + + Double_t linf=0; + Double_t lsup=0; + if(nfound == m_NumberOfGammaPeak){ + cout << "Gamma Peak Found" << endl; + for(int i=0; i<nfound; i++){ + linf = xpeaks[i]-1.1; + lsup = xpeaks[i]+0.8; + /* if(anode==9){ */ + /* lsup =xpeaks[i]+0.8 */ + /* } */ + + TF1* gauss = new TF1("gaus","gaus",linf,lsup); + h->Fit(gauss,"RQ"); + mean[i] = gauss->GetParameter(1); + sigma[i] = gauss->GetParameter(2); + } + for(int i=0; i<nfound; i++){ + linf = xpeaks[i]-1.3*sigma[i]; + lsup = xpeaks[i]+1*sigma[i]; + TF1* gauss = new TF1("gaus","gaus",linf,lsup); + h->Fit(gauss,"RQ"); + mean[i] = gauss->GetParameter(1); + sigma[i] = gauss->GetParameter(2); + } + + return true; + } + + else if(nfound != m_NumberOfGammaPeak){ + cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl; + return false; + } +}