diff --git a/Projects/Vendeta/macro/FillTOFHisto.C b/Projects/Vendeta/macro/FillTOFHisto.C
index 558383cc70048659cedcb3ecf36b9062859766d7..144a783605a0c9a1c6c70cd6e00a39d061717654 100644
--- a/Projects/Vendeta/macro/FillTOFHisto.C
+++ b/Projects/Vendeta/macro/FillTOFHisto.C
@@ -1,104 +1,130 @@
-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 NumberOfAnodes= 11;
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(){
-
- 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];
-
- 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>();
-
- 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>();
-
- TFissionChamberPhysics* FC = new TFissionChamberPhysics();
- chain->SetBranchAddress("FissionChamber",&FC);
-
- 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);
- chain->SetBranchAddress("HG_Tof",&HG_Tof);
- chain->SetBranchAddress("HG_ID",&HG_ID);
- chain->SetBranchAddress("HG_Q1",&HG_Q1);
- chain->SetBranchAddress("HG_Q2",&HG_Q2);
-
- 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);
- }
- }
- for(int i=0; i<nentries; i++){
- chain->GetEntry(i);
-
- if(i%100000==0){
- cout << setprecision(2) << "\033[34m\r Processing tree..." << (double)i/nentries*100 << "\% done" << flush;
- }
-
- 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){
- hLG[index_LG]->Fill(LG_Tof->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));
- }
- }
- }
- }
-
- ofile->Write();
- ofile->Close();
+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);");
+
+ 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>();
+ vector<double>* LG_FlyPath = new vector<double>();
+
+ 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>();
+ vector<double>* HG_FlyPath = new vector<double>();
+
+ 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);
+ chain->SetBranchAddress("LG_FlyPath",&LG_FlyPath);
+ chain->SetBranchAddress("HG_Tof",&HG_Tof);
+ chain->SetBranchAddress("HG_ID",&HG_ID);
+ chain->SetBranchAddress("HG_Q1",&HG_Q1);
+ chain->SetBranchAddress("HG_Q2",&HG_Q2);
+ 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<nentries; i++){
+ chain->GetEntry(i);
+
+ if(i%100000==0){
+ cout << setprecision(2) << "\033[34m\r Processing run "<< run <<" ..." << (double)i/nentries*100 << "\% done" << flush;
+ }
+
+ 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){
+ 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();
}
diff --git a/Projects/Vendeta/macro/FitTofGammaPeak.C b/Projects/Vendeta/macro/FitTofGammaPeak.C
index 62d4967acbb704150e0ec3efccbc886ea0751966..d008591bcfdfa60d4e72323eac059efe1abc1e9c 100644
--- a/Projects/Vendeta/macro/FitTofGammaPeak.C
+++ b/Projects/Vendeta/macro/FitTofGammaPeak.C
@@ -1,43 +1,61 @@
-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 NumberOfAnodes=11;
int m_NumberOfGammaPeak=1;
-int run=86;
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);
/////////////////////////////////////////////////////
-void OpenRootFile(){
- ifile = new TFile(Form("histo_tof_file_run%i.root",run));
-}
-
-/////////////////////////////////////////////////////
-void FitTofGammaPeak(){
- OpenRootFile();
+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;
- ofile.open(Form("Vendeta_Time_run%i.cal",run));
-
- TFile* orootfile = new TFile(Form("histo_tof_fitted_run%i.root",run),"recreate");
+ 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[i]=0;
+ sigma_anode_HG[i]=0;
}
Double_t* mean = new Double_t[1];
Double_t* sigma = new Double_t[1];
- TGraph* gSigma_LG = new TGraph();
- TGraph* gSigma_HG = new TGraph();
+ TGraph* gSigma_LG = new TGraph();
+ TGraph* gSigma_HG = new TGraph();
- int countLG=0;
- int countHG=0;
+ int countLG=0;
+ int countHG=0;
for(int i=0; i<NumberOfDetectors; i++){
for(int j=0; j<NumberOfAnodes; j++){
// LG //
@@ -46,18 +64,22 @@ void FitTofGammaPeak(){
//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]);
+ gSigma_LG->SetPoint(countLG,countLG+1,sigma[0]);
- sigma_anode[j] += sigma[0];
+ sigma_anode[j] += sigma[0];
- countLG++;
- h->Write();
+ countLG++;
+ h->Write();
}
else{
ofile << LG_token << " 0" << endl;
@@ -75,9 +97,10 @@ void FitTofGammaPeak(){
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]);
- countHG++;
- h->Write();
+ gSigma_HG->SetPoint(countHG,countHG+1,sigma[0]);
+ sigma_anode_HG[j] += sigma[0];
+ countHG++;
+ h->Write();
}
else{
ofile << HG_token << " 0" << endl;
@@ -85,30 +108,35 @@ void FitTofGammaPeak(){
}
}
- TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
- c1->Divide(2,1);
+ TCanvas* c1 = new TCanvas("Sigma","Sigma",1200,600);
+ c1->Divide(2,1);
- gSigma_LG->SetMarkerStyle(8);
- gSigma_HG->SetMarkerStyle(8);
+ gSigma_LG->SetMarkerStyle(8);
+ gSigma_HG->SetMarkerStyle(8);
- c1->cd(1);
- gSigma_LG->Draw();
- c1->cd(2);
- gSigma_HG->Draw();
+ 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();
+ gSigma_LG->SetName("sigma_LG");
+ gSigma_HG->SetName("sigma_HG");
- orootfile->Write();
- orootfile->Close();
+ gSigma_LG->Write();
+ gSigma_HG->Write();
- for(int i=0; i<NumberOfAnodes; i++){
- cout << "Anode= " << i+1 << " / " << sigma_anode[i]/NumberOfDetectors << endl;
- }
+ orootfile->Write();
+ 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){
if(nfound == m_NumberOfGammaPeak){
cout << "Gamma Peak Found" << endl;
for(int i=0; i<nfound; i++){
- linf = xpeaks[i]-2;
- lsup = xpeaks[i]+1;
+ 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){
+ else if(nfound != m_NumberOfGammaPeak){
cout << "Warning. Number of peak different of " << m_NumberOfGammaPeak << " !! / nfound = " << nfound << endl;
return false;
}