updating calib. macros in Vendeta project

Projects/Vendeta/macro/FillTOFHisto.C

+// 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();
+
+}

Projects/Vendeta/macro/FitTofGammaPeak.C

+// 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;
+  }
+}