Skip to content
Snippets Groups Projects
Commit 968cd667 authored by Cyril Lenain's avatar Cyril Lenain :surfer_tone3:
Browse files

Updating calib. macros in Vendeta Project

parent d44869d9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Projects/Vendeta/macro/FillTOFHisto.C
+ 122
− 96
View file @ 968cd667
TChain* chain; // 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 NumberOfDetectors= 72;
int NumberOfAnodes= 11; int NumberOfAnodes= 11;
int nentries=1e6; int nentries=1e6;
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(){ int main(int argc, char** argv){
int run = atof(argv[1]);
LoadRootFile();
nentries = chain->GetEntries(); auto chain = new TChain("PhysicsTree");
cout << "Number of entries: " << nentries << endl; chain->Add(Form("/home/cyril/Workflow/nptool/Outputs/Analysis/run%i/run%i.root",run,run));
TFile* ofile = new TFile(Form("histo_tof_file_run%i.root",run),"recreate"); nentries = chain->GetEntries();
TH1F* hLG[791]; /* nentries = 50000000; */
TH1F* hHG[791]; cout << "Number of entries: " << nentries << endl;
vector<double>* FC_Q1 = new vector<double>(); TFile* ofile = new TFile(Form("histos_ToF/histo_tof_file_run%i.root",run),"recreate");
vector<int>* FC_Anode_ID = new vector<int>(); TH2F* hLG_LSci[11];
vector<bool>* FC_FakeFission = new vector<bool>(); TH2F* hHG_LSci[11];
TH1F* hLG[791];
vector<double>* LG_Tof = new vector<double>(); TH1F* hHG[791];
vector<int>* LG_ID = new vector<int>();
vector<double>* LG_Q1 = new vector<double>();
vector<double>* LG_Q2 = new vector<double>(); TGraph *gFlyPath = new TGraph();
gFlyPath->SetTitle(" ; Vendeta-Anode Index ; fly length (mm);");
vector<double>* HG_Tof = new vector<double>();
vector<int>* HG_ID = new vector<int>(); vector<double>* FC_Q1 = new vector<double>();
vector<double>* HG_Q1 = new vector<double>(); vector<int>* FC_Anode_ID = new vector<int>();
vector<double>* HG_Q2 = new vector<double>(); vector<bool>* FC_FakeFission = new vector<bool>();
TFissionChamberPhysics* FC = new TFissionChamberPhysics(); vector<double>* LG_Tof = new vector<double>();
chain->SetBranchAddress("FissionChamber",&FC); vector<int>* LG_ID = new vector<int>();
vector<double>* LG_Q1 = new vector<double>();
chain->SetBranchAddress("FC_Q1",&FC_Q1); vector<double>* LG_Q2 = new vector<double>();
chain->SetBranchAddress("FC_Anode_ID",&FC_Anode_ID); vector<double>* LG_FlyPath = new vector<double>();
chain->SetBranchAddress("FC_FakeFission",&FC_FakeFission);
chain->SetBranchAddress("LG_Tof",&LG_Tof); vector<double>* HG_Tof = new vector<double>();
chain->SetBranchAddress("LG_ID",&LG_ID); vector<int>* HG_ID = new vector<int>();
chain->SetBranchAddress("LG_Q1",&LG_Q1); vector<double>* HG_Q1 = new vector<double>();
chain->SetBranchAddress("LG_Q2",&LG_Q2); vector<double>* HG_Q2 = new vector<double>();
chain->SetBranchAddress("HG_Tof",&HG_Tof); vector<double>* HG_FlyPath = new vector<double>();
chain->SetBranchAddress("HG_ID",&HG_ID);
chain->SetBranchAddress("HG_Q1",&HG_Q1); chain->SetBranchAddress("FC_Q1",&FC_Q1);
chain->SetBranchAddress("HG_Q2",&HG_Q2); chain->SetBranchAddress("FC_Anode_ID",&FC_Anode_ID);
chain->SetBranchAddress("FC_FakeFission",&FC_FakeFission);
for(int i=0; i<NumberOfDetectors; i++){ chain->SetBranchAddress("LG_Tof",&LG_Tof);
for(int j=0; j<NumberOfAnodes; j++){ chain->SetBranchAddress("LG_ID",&LG_ID);
int index = j + i*NumberOfAnodes; chain->SetBranchAddress("LG_Q1",&LG_Q1);
TString histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1); chain->SetBranchAddress("LG_Q2",&LG_Q2);
hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300); chain->SetBranchAddress("LG_FlyPath",&LG_FlyPath);
chain->SetBranchAddress("HG_Tof",&HG_Tof);
histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1); chain->SetBranchAddress("HG_ID",&HG_ID);
hHG[index] = new TH1F(histo_name,histo_name,2500,0,500); chain->SetBranchAddress("HG_Q1",&HG_Q1);
} chain->SetBranchAddress("HG_Q2",&HG_Q2);
} chain->SetBranchAddress("HG_FlyPath",&HG_FlyPath);
for(int i=0; i<nentries; i++){
chain->GetEntry(i); for(int j=0; j<NumberOfAnodes; j++){
if(i%100000==0){ TString histo_name = Form("hLG_LSci_Anode%i",j+1);
cout << setprecision(2) << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush; hLG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300);
}
histo_name = Form("hHG_LSci_Anode%i",j+1);
if(FC_Anode_ID->size()>0){ hHG_LSci[j] = new TH2F(histo_name,histo_name,72,1,73,2000,-100,300);
bool Fake = FC_FakeFission->at(0);
int anode = FC_Anode_ID->at(0); for(int i=0; i<NumberOfDetectors; i++){
int index = j + i*NumberOfAnodes;
histo_name = Form("hLG_Det%i_Anode%i",i+1,j+1);
int mysize = LG_Tof->size(); hLG[index] = new TH1F(histo_name,histo_name,2000,-100,300);
for(int j=0; j<mysize; j++){
// LG // histo_name = Form("hHG_Det%i_Anode%i",i+1,j+1);
int index_LG = (anode-1) + (LG_ID->at(j)-1)*NumberOfAnodes; hHG[index] = new TH1F(histo_name,histo_name,2000,-100,300);
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){ }
hLG[index_LG]->Fill(LG_Tof->at(j)); for(int i=0; i<nentries; i++){
} chain->GetEntry(i);
}
if(i%100000==0){
mysize = HG_Tof->size(); cout << setprecision(2) << "\033[34m\r Processing run "<< run <<" ..." << (double)i/nentries*100 << "\% done" << flush;
for(int j=0; j<mysize; j++){ }
// HG //
int index_HG = (anode-1) + (HG_ID->at(j)-1)*NumberOfAnodes; if(FC_Anode_ID->size()>0){
double PSD = HG_Q2->at(j)/HG_Q1->at(j); bool Fake = FC_FakeFission->at(0);
if(HG_ID->at(j)>0 && anode>0 && Fake==0 && HG_ID->size()==1){ int anode = FC_Anode_ID->at(0);
hHG[index_HG]->Fill(HG_Tof->at(j));
}
} int mysize = LG_Tof->size();
} for(int j=0; j<mysize; j++){
} // LG //
int index_LG = (anode-1) + (LG_ID->at(j)-1)*NumberOfAnodes;
ofile->Write(); double PSD = LG_Q2->at(j)/LG_Q1->at(j);
ofile->Close(); if(LG_ID->at(j)>0 && anode>0 && Fake==0 && LG_Q1->at(j)>2000){
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));
}
}
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){
hHG[index_HG]->Fill(HG_Tof->at(j));
hHG_LSci[anode-1]->Fill(HG_ID->at(j),HG_Tof->at(j));
}
}
}
}
ofile->WriteObject(gFlyPath,"gFlyPath");
ofile->Write();
ofile->Close();
} }
Projects/Vendeta/macro/FitTofGammaPeak.C
+ 87
− 47
View file @ 968cd667
TFile* ifile; // 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 NumberOfDetectors=72;
int NumberOfAnodes=11; int NumberOfAnodes=11;
int m_NumberOfGammaPeak=1; int m_NumberOfGammaPeak=1;
int run=86;
double PosGammaPeak = 3.2; double PosGammaPeak = 3.2;
double c = 299.792458; // mm/ns
double sigma_anode[11]; double sigma_anode[11], sigma_anode_HG[11];
bool Finder(TH1F* h, Double_t *mean, Double_t *sigma); bool Finder(TH1F* h, Double_t *mean, Double_t *sigma);
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
void OpenRootFile(){ int main(int argc, char** argv){
ifile = new TFile(Form("histo_tof_file_run%i.root",run));
} int run = atof(argv[1]);
TString Append = argv[2];
/////////////////////////////////////////////////////
void FitTofGammaPeak(){ auto ifile = new TFile(Form("histos_ToF/histo_tof_file_run%d"+Append+".root",run));
OpenRootFile(); auto gFlyPath = (TGraph*)ifile->Get("gFlyPath");
TCanvas *cLG[11]; TCanvas *cLG[11];
TCanvas *cHG[11]; TCanvas *cHG[11];
ofstream ofile; ofstream ofile;
ofile.open(Form("Vendeta_Time_run%i.cal",run)); 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");
TFile* orootfile = new TFile(Form("histo_tof_fitted_run%i.root",run),"recreate");
for(int i=0; i<NumberOfAnodes; i++){ for(int i=0; i<NumberOfAnodes; i++){
sigma_anode[i]=0; sigma_anode[i]=0;
sigma_anode_HG[i]=0;
} }
Double_t* mean = new Double_t[1]; Double_t* mean = new Double_t[1];
Double_t* sigma = new Double_t[1]; Double_t* sigma = new Double_t[1];
TGraph* gSigma_LG = new TGraph(); TGraph* gSigma_LG = new TGraph();
TGraph* gSigma_HG = new TGraph(); TGraph* gSigma_HG = new TGraph();
int countLG=0; int countLG=0;
int countHG=0; int countHG=0;
for(int i=0; i<NumberOfDetectors; i++){ for(int i=0; i<NumberOfDetectors; i++){
for(int j=0; j<NumberOfAnodes; j++){ for(int j=0; j<NumberOfAnodes; j++){
// LG // // LG //
...@@ -46,18 +64,22 @@ void FitTofGammaPeak(){ ...@@ -46,18 +64,22 @@ void FitTofGammaPeak(){
//cLG[j]->cd(i+1); //cLG[j]->cd(i+1);
h->Draw(); h->Draw();
int index = j + i*11;
double FlyPath = gFlyPath->GetPointY(index);
PosGammaPeak = FlyPath/c ;
mean[0] = 0; mean[0] = 0;
sigma[0] = 0; sigma[0] = 0;
TString LG_token = Form("Vendeta_DET%i_LG_ANODE%i_TIMEOFFSET",i+1,j+1); TString LG_token = Form("Vendeta_DET%i_LG_ANODE%i_TIMEOFFSET",i+1,j+1);
if(Finder(h, mean, sigma)){ if(Finder(h, mean, sigma)){
ofile << LG_token << " " << -mean[0]+PosGammaPeak << endl; ofile << LG_token << " " << -mean[0]+PosGammaPeak << endl;
gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]); gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]);
sigma_anode[j] += sigma[0]; sigma_anode[j] += sigma[0];
countLG++; countLG++;
h->Write(); h->Write();
} }
else{ else{
ofile << LG_token << " 0" << endl; ofile << LG_token << " 0" << endl;
...@@ -75,9 +97,10 @@ void FitTofGammaPeak(){ ...@@ -75,9 +97,10 @@ void FitTofGammaPeak(){
TString HG_token = Form("Vendeta_DET%i_HG_ANODE%i_TIMEOFFSET",i+1,j+1); TString HG_token = Form("Vendeta_DET%i_HG_ANODE%i_TIMEOFFSET",i+1,j+1);
if(Finder(h, mean, sigma)){ if(Finder(h, mean, sigma)){
ofile << HG_token << " " << -mean[0]+PosGammaPeak << endl; ofile << HG_token << " " << -mean[0]+PosGammaPeak << endl;
gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]); gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]);
countHG++; sigma_anode_HG[j] += sigma[0];
h->Write(); countHG++;
h->Write();
} }
else{ else{
ofile << HG_token << " 0" << endl; ofile << HG_token << " 0" << endl;
...@@ -85,30 +108,35 @@ void FitTofGammaPeak(){ ...@@ -85,30 +108,35 @@ void FitTofGammaPeak(){
} }
} }
TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600); TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
c1->Divide(2,1); c1->Divide(2,1);
gSigma_LG->SetMarkerStyle(8); gSigma_LG->SetMarkerStyle(8);
gSigma_HG->SetMarkerStyle(8); gSigma_HG->SetMarkerStyle(8);
c1->cd(1); c1->cd(1);
gSigma_LG->Draw(); gSigma_LG->Draw();
c1->cd(2); c1->cd(2);
gSigma_HG->Draw(); gSigma_HG->Draw();
gSigma_LG->SetName("sigma_LG"); gSigma_LG->SetName("sigma_LG");
gSigma_HG->SetName("sigma_HG"); gSigma_HG->SetName("sigma_HG");
gSigma_LG->Write();
gSigma_HG->Write();
orootfile->Write(); gSigma_LG->Write();
orootfile->Close(); gSigma_HG->Write();
for(int i=0; i<NumberOfAnodes; i++){ orootfile->Write();
cout << "Anode= " << i+1 << " / " << sigma_anode[i]/NumberOfDetectors << endl; orootfile->Close();
}
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;
} }
...@@ -128,18 +156,30 @@ bool Finder(TH1F* h, Double_t *mean, Double_t *sigma){ ...@@ -128,18 +156,30 @@ bool Finder(TH1F* h, Double_t *mean, Double_t *sigma){
if(nfound == m_NumberOfGammaPeak){ if(nfound == m_NumberOfGammaPeak){
cout << "Gamma Peak Found" << endl; cout << "Gamma Peak Found" << endl;
for(int i=0; i<nfound; i++){ for(int i=0; i<nfound; i++){
linf = xpeaks[i]-2; linf = xpeaks[i]-1.1;
lsup = xpeaks[i]+1; lsup = xpeaks[i]+0.8;
/* if(anode==9){ */
/* lsup =xpeaks[i]+0.8 */
/* } */
TF1* gauss = new TF1("gaus","gaus",linf,lsup); TF1* gauss = new TF1("gaus","gaus",linf,lsup);
h->Fit(gauss,"RQ"); h->Fit(gauss,"RQ");
mean[i] = gauss->GetParameter(1); mean[i] = gauss->GetParameter(1);
sigma[i] = gauss->GetParameter(2); 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; return true;
} }
else if(nfound != m_NumberOfGammaPeak){ else if(nfound != m_NumberOfGammaPeak){
cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl; cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl;
return false; return false;
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment