From c1cbe3d55308cefece12cdbd418833ba159e78ab Mon Sep 17 00:00:00 2001
From: "theodore.efremov" <theodore.efremov@cea.fr>
Date: Thu, 31 Oct 2024 14:10:02 +0100
Subject: [PATCH] Added two macro, one to make cut on the DE_E and the other to
apply a decreasing expo fit on each segment of the IC
---
.../2024/macro/chio/CalibrationChio/CutZ.cxx | 128 +++++++++++++
.../2024/macro/chio/CalibrationChio/ICCorr.C | 13 +-
.../IC_Dampen_Recombination.cxx | 168 ++++++++++++++++++
3 files changed, 306 insertions(+), 3 deletions(-)
create mode 100644 Projects/AlPhaPha/2024/macro/chio/CalibrationChio/CutZ.cxx
create mode 100644 Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/CutZ.cxx b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/CutZ.cxx
new file mode 100644
index 000000000..103aa7acb
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/CutZ.cxx
@@ -0,0 +1,128 @@
+#include "IC_Dampen_Recombination.cxx"
+
+
+void HistoFillerDE_E();
+/**
+ * * Draw the DE_E plot and save the histo into a file
+ * */
+void CutZ(bool GenDE=true){
+
+
+ if (GenDE==true) HistoFillerDE_E() ;
+
+
+ TFile* infile = TFile::Open("Output/DE_E.root");
+
+
+ //===========================================================================================================
+ // Load Histo
+ //===========================================================================================================
+ // Get number of TH2F
+ int HistoCount = 0 ;
+ TIter next(infile->GetListOfKeys());
+ TKey* key;
+
+ while ((key=(TKey*)next())){
+ if (std::string(key->GetClassName()) == "TH2F"){
+ HistoCount ++;
+ }
+ }
+
+ vector<TH2F*> hDE_E;
+ //for each histo in the file
+ for (int i=0 ; i< HistoCount ; i++){
+ if (i==0) hDE_E.push_back((TH2F*)infile->Get("DE_E_online"));
+ if (i==1) hDE_E.push_back((TH2F*)infile->Get("DE_Etot_online"));
+ }
+
+ //===========================================================================================================
+ // Draw histo
+ //===========================================================================================================
+
+ TCanvas *cDE = new TCanvas("cDE","cDE",1000,500);
+ cDE->Divide(2);
+ cDE->cd(1);
+ hDE_E[0]->Draw();
+ cDE->cd(2);
+ hDE_E[1]->Draw();
+}
+
+void HistoFillerDE_E(){
+
+
+ //===========================================================================================================
+ // Loading var
+ //===========================================================================================================
+ TChain* chain = new TChain("PhysicsTree");
+ chain->Add("../../../root/analysis/VamosCalib241.root");
+
+
+ TICPhysics* IC = new TICPhysics() ;
+ double FF_IC_X, FF_IC_Y;
+
+ chain->SetBranchStatus("FF_IC_X", true);
+ chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
+ chain->SetBranchStatus("FF_IC_Y", true);
+ chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+ chain->SetBranchStatus("IC", true);
+ chain->SetBranchAddress("IC", &IC);
+
+ //===========================================================================================================
+ // Defining histo
+ //===========================================================================================================
+
+ TH2F* hDE_E_nocorr = new
+ TH2F("DE_E_online","DE_E_online",1000,0,20000,1000,0,25000);
+ TH2F* hDE_Etot = new
+ TH2F("DE_Etot_online","DE_Etot_online",1000,0,50000,1000,0,25000);
+ //===========================================================================================================
+ // FILL HISTO
+ //===========================================================================================================
+
+ int Nentries = chain->GetEntries();
+
+ for (int e = 0; e < Nentries; e++) {
+ // Printing status
+ if (e % 100000 == 0){
+ cout << "Please wait, we are " << e * 100. / Nentries << "% done ..\r" << flush;
+ }
+ chain->GetEntry(e);
+
+ double DE = 0 ,E = 0, Etot = 0 ;
+
+ //===========================================================================================================
+ // Fill De from seg 1 to 3
+ //===========================================================================================================
+ for (int seg = 0 ; seg < sizeof(IC->fIC_raw)/sizeof(IC->fIC_raw[0]) ; seg++ ){
+ if(seg > 1 && seg < 4) DE += IC->fIC_raw[seg] ;
+ else if (seg >4) E+= IC->fIC_raw[seg] ;
+ }
+
+
+
+ //===========================================================================================================
+ // Filling DE with ICO and 4
+ //===========================================================================================================
+
+ DE += IC->fIC_raw[0] ;
+
+ DE = DE *0.5 + IC->fIC_raw[4];
+
+ double IC1corr = (IC->fIC_raw[1]*(1-0.000686068*FF_IC_Y))*(1-4.88238e-05*FF_IC_Y+7.40395e-06*FF_IC_Y*FF_IC_Y);
+ double DE_Bis = 0.5*(IC1corr+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
+
+ Etot = E+DE_Bis;
+ //===========================================================================================================
+ // Fill histo
+ //===========================================================================================================
+
+ hDE_E_nocorr->Fill(E,DE_Bis);
+ hDE_Etot->Fill(Etot,DE_Bis);
+ }
+
+ TFile* outfile = new TFile("Output/DE_E.root","recreate");
+ hDE_E_nocorr->Write();
+ hDE_Etot->Write();
+ outfile->Close();
+
+}
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
index 5521ceb21..80f5b4ab3 100644
--- a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/ICCorr.C
@@ -51,7 +51,9 @@ void ICCorr() {
TH2F* hDE_E_IC0_Corr = new
TH2F("hDE_E_IC0_Corr","hDE_E_IC0_Corr",1000,0,20000,1000,0,25000);
-
+ TH2F* hDEcut=
+ new TH2F("DEcut","DEcut",1000,0,20000,1000,0,25000);
+
hICorr.push_back(hICorr_X);
hICorr.push_back(hICorr_Y);
@@ -169,6 +171,7 @@ void ICCorr() {
// Chio online
//===========================================================================================================
+
double IC1corr = (IC->fIC_raw[1]*(1-0.000686068*FF_IC_Y))*(1-4.88238e-05*FF_IC_Y+7.40395e-06*FF_IC_Y*FF_IC_Y);
double DE_Bis = 0.5*(IC1corr+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
//double DE_Bis = 0.5*(IC->fIC_raw[1]+IC->fIC_raw[2]+IC->fIC_raw[3])+IC->fIC_raw[4];
@@ -178,7 +181,11 @@ void ICCorr() {
hDE_E->Fill(E,DE);
hDE_Ecorr01->Fill(E,DEcorr);
hDE_E_IC1_Corr->Fill(E,DEIC1corr);
- hDE_E_IC0_Corr->Fill(E,DE_IC0_Corr_Y);
+
+ double DEcut;
+ if ((FF_IC_Y>-60 && FF_IC_Y<-55)|| (FF_IC_Y>55 && FF_IC_Y<60)){ DEcut=DE_Bis;
+ hDEcut->Fill(E,DEcut);
+ }
} //End loop on event
//===========================================================================================================
@@ -205,6 +212,6 @@ void ICCorr() {
c2->cd(3);
hDEBis_Ecorr01->Draw();
c2->cd(4);
- hDE_E_IC0_Corr->Draw();
+ hDEcut->Draw();
} // End spline chio XY
diff --git a/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
new file mode 100644
index 000000000..05ce6484b
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/chio/CalibrationChio/IC_Dampen_Recombination.cxx
@@ -0,0 +1,168 @@
+#include <TFile.h>
+#include <TICPhysics.h>
+#include <TChain.h>
+#include <TF1.h>
+#include <TH2.h>
+#include <Math/Functor.h>
+#include <TTreeReaderValue.h>
+#include <vector>
+using namespace std;
+
+
+void HistoFillerIC();
+
+void IC_Dampen_Recombination(bool Fill = false){
+
+
+ if (Fill==true) HistoFillerIC();
+
+
+ TFile *infile = TFile::Open("Output/HistoDampen.root");
+ vector<TH2F*> hICX, hICY ;
+ vector<TH1F*> hICX_Profile, hICY_Profile ;
+
+ for (int seg = 0 ; seg < 11; seg++ ){
+
+ hICX.push_back((TH2F*)infile->Get(Form("hChio_X_%d",seg)));
+ hICY.push_back((TH2F*)infile->Get(Form("hChio_Y_%d",seg)));
+ hICX_Profile.push_back((TH1F*)infile->Get(Form("hChio_X_%d_pfx",seg)));
+ hICY_Profile.push_back((TH1F*)infile->Get(Form("hChio_Y_%d_pfx",seg)));
+
+ }
+
+ //===========================================================================================================
+ // Fit Histo
+ //===========================================================================================================
+
+
+
+ vector<TF1*> Fit_Y(11);
+ for (int seg=0; seg<11 ; seg++){
+
+ Fit_Y[seg] = new
+ TF1(Form("Fit_Y_Expo_IC%d",seg),[](double *x, double*p) {
+
+ if (x[0]>-50 && x[0]<50) {
+ TF1::RejectPoint();
+ }
+ //Define here the function
+ return p[0] * exp(p[1]* (x[0]-p[3])) + p[2] ;
+ },-60,63,4);
+
+
+ Fit_Y[seg]->SetParameters(1,-10,4800,-60);
+
+
+ cout << "****************" << endl << "IC"<<seg << endl;
+ hICY_Profile[seg]->Fit(Fit_Y[seg] ,"Rnone"); vector<TH2F*> hICX, hICY ;
+ vector<TH1F*> hICX_Profile, hICY_Profile ;
+
+
+
+ }//end loop on histo
+
+
+
+ //===========================================================================================================
+ // Draw histo
+ //===========================================================================================================
+
+ vector<TCanvas*> cIC(11);
+ for (int seg=0; seg<1; seg++){
+ cIC[seg] = new TCanvas
+ (Form("cIC%d",seg),Form("cIC%d",seg),100,100,1000,600);
+ cIC[seg]->Divide(2);
+
+ cIC[seg]->cd(1);
+ hICY[seg]->Draw();
+
+ cIC[seg]->cd(2);
+ hICY_Profile[seg]->Draw();
+ Fit_Y[seg]->Draw("same");
+
+
+ }
+}
+
+
+void HistoFillerIC(){
+ //===========================================================================================================
+ // INPUT
+ //===========================================================================================================
+
+ TChain* chain = new TChain("PhysicsTree");
+ chain->Add("../../../root/analysis/VamosCalib241.root");
+
+ double FF_IC_X, FF_IC_Y;
+ TICPhysics* IC = new TICPhysics();
+
+ chain->SetBranchStatus("FF_IC_X", true);
+ chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
+ chain->SetBranchStatus("FF_IC_Y", true);
+ chain->SetBranchAddress("FF_IC_Y", &FF_IC_Y);
+ chain->SetBranchStatus("IC", true);
+ chain->SetBranchAddress("IC", &IC);
+
+
+ //===========================================================================================================
+ // Declare and fill histo
+ //===========================================================================================================
+
+ vector<TH2F*> hICX, hICY ;
+ vector<TH1F*> hICX_Profile, hICY_Profile ;
+
+ int NSegmentChio = 11 ;
+
+ for (int seg = 0 ; seg < 11 ; seg++){
+
+ TH2F* hChio_X = new TH2F(Form("hChio_X_%d",seg) ,Form("hChio_X_%d",seg) ,1000 ,-800 ,800 ,500 ,2000,9500 );
+ TH2F* hChio_Y = new TH2F(Form("hChio_Y_%d",seg) ,Form("hChio_Y_%d",seg) ,1000 ,-160 ,160 ,500 ,2000,9500 );
+
+ TH1F* hProfile_X = (TH1F*)hChio_X->ProfileX();
+ TH1F* hProfile_Y = (TH1F*)hChio_Y->ProfileX();
+
+ hICX.push_back(hChio_X);
+ hICY.push_back(hChio_Y);
+ hICX_Profile.push_back(hProfile_X);
+ hICY_Profile.push_back(hProfile_Y);
+
+ } //end loop on segment
+
+
+ // Loop on event
+ int Nentries = chain->GetEntries();
+
+ for (int e = 0 ; e<Nentries ; e++){ // to be parallelised
+
+ if ((e/Nentries)*100 %100 ==0){
+ cout << "Please wait, we are " << e * 100. / Nentries << "% done ..\r" << flush;
+ }
+
+ chain->GetEntry(e);
+
+ for (int seg=0 ; seg<11 ; seg++ ){
+
+ hICX[seg]->Fill(FF_IC_X,IC->fIC_raw[seg]);
+ hICY[seg]->Fill(FF_IC_Y,IC->fIC_raw[seg]);
+
+ }//end loop seg
+
+ } //end loop entries
+
+
+ //ReProfile
+
+ for (int seg=0; seg<11 ;seg++){
+ hICX_Profile[seg] = (TH1F*)hICX[seg]->ProfileX();
+ hICY_Profile[seg] = (TH1F*)hICY[seg]->ProfileX();
+ }
+
+ TFile* outfile = new TFile("Output/HistoDampen.root","recreate");
+ for (int seg=0 ; seg<11 ;seg++){
+ hICX[seg]->Write();
+ hICY[seg]->Write();
+ hICX_Profile[seg]->Write();
+ hICY_Profile[seg]->Write();
+ }
+ outfile->Close();
+}
--
GitLab