diff --git a/NPAnalysis/e530/Makefile b/NPAnalysis/e530/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..f5bfeb69ffd25b28c9ccb8a81eebb1b66664dfa2 --- /dev/null +++ b/NPAnalysis/e530/Makefile @@ -0,0 +1,6 @@ + +Analyse: + make -C ./src + +clean: + make clean -C ./src diff --git a/NPAnalysis/e530/RunToTreat.txt b/NPAnalysis/e530/RunToTreat.txt new file mode 100644 index 0000000000000000000000000000000000000000..1a9cc7e2870b0895210c65798057d40d0dcfe585 --- /dev/null +++ b/NPAnalysis/e530/RunToTreat.txt @@ -0,0 +1,5 @@ +TTreeName + AutoTree +RootFileName + ../../../../root/run_0003.root + diff --git a/NPAnalysis/e530/include/ObjectManager.hh b/NPAnalysis/e530/include/ObjectManager.hh new file mode 100644 index 0000000000000000000000000000000000000000..036e1a34f65df03c02e04424554c6185972b913c --- /dev/null +++ b/NPAnalysis/e530/include/ObjectManager.hh @@ -0,0 +1,194 @@ +// You can use this file to declare your spectra, file, energy loss , ... and whatever you want. +// This way you can remove all unnecessary declaration in the main programm. +// In order to help debugging and organizing we use Name Space. + +///////////////////////////////////////////////////////////////////////////////////////////////// +// -------------------------------------- VARIOUS INCLUDE --------------------------------------- + +// NPA +#include "DetectorManager.h" + +// STL C++ +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +#include <cmath> +#include <cstdlib> +#include <time.h> + +// ROOT +#include <TROOT.h> +#include <TChain.h> +#include <TFile.h> +#include <TLeaf.h> +#include <TVector3.h> +#include <TRandom.h> + +// NPL +#include "TPlasticData.h" +#include "NPReaction.h" +#include "RootInput.h" +#include "RootOutput.h" +#include "TInitialConditions.h" +#include "TMust2Physics.h" +#include "TSSSDPhysics.h" +#include "TPlasticPhysics.h" + +// Use CLHEP System of unit and Physical Constant +#include "CLHEP/Units/GlobalSystemOfUnits.h" +#include "CLHEP/Units/PhysicalConstants.h" + + +// ---------------------------------------------------------------------------------------------- +double ThetaCalculation (TVector3 A , TVector3 B) ; +///////////////////////////////////////////////////////////////////////////////////////////////// +// ----------------------------------- DOUBLE, INT, BOOL AND MORE ------------------------------- +namespace VARIABLE + { + // Declare your Variable here: + + double X1,Y1,Z1 ; + int N1,N2 = 0 ; + bool check= false ; + + // A Usefull Simple Random Generator + TRandom Rand; + } + +using namespace VARIABLE ; +// ---------------------------------------------------------------------------------------------- + + + +///////////////////////////////////////////////////////////////////////////////////////////////// +// -----------------------------------GRAPH------------------------------------------------------ +#include <TObject.h> +#include <TH1.h> +#include <TH1F.h> +#include <TH2.h> +#include <TH2F.h> +#include <TGraph2D.h> + +namespace GRAPH + { + // Declare your Spectra here: + + TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ; + + TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ; + + } + +using namespace GRAPH ; +// -------------------------------------------------------------------------------------------- + + + +/////////////////////////////////////////////////////////////////////////////////////////////// +// -----------------------------------CUT------------------------------------------------------ +#include <TCutG.h> +namespace CUT + { + // Declare your Cut here: + + } + +using namespace CUT ; +// -------------------------------------------------------------------------------------------- + + + +//////////////////////////////////////////////////////////////////////////////////////////////// +// -----------------------------------ENERGY LOSS---------------------------------------------- +#include "NPEnergyLoss.h" +using namespace NPL ; +namespace ENERGYLOSS + { + + + // 3He Energy Loss + /* + EnergyLoss He3TargetWind = EnergyLoss ("He3_Mylar.G4table", + "G4Table" , + 10000 ); + + EnergyLoss He3TargetGaz = EnergyLoss ("He3_D2.G4table" , + "G4Table" , + 10000 ); + + EnergyLoss He3StripAl = EnergyLoss ("3He_Al.txt" , + "LISE" , + 10000 , + 1 , + 3 ); + + EnergyLoss He3StripSi = EnergyLoss ("3He_Si.txt" , + "LISE" , + 10000 , + 1 , + 3 ); + + */ + + // // 3He Energy Loss + // EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" , + // 10000 , + // 1 , + // 3 ); + // + // EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" , + // 10000 , + // 1 , + // 3 ); + // + // EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" , + // 10000 , + // 1 , + // 3 ); + // + // EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" , + // 10000 , + // 1 , + // 3 ); + + + // proton Energy Loss + EnergyLoss protonTargetWind = EnergyLoss ("proton_Mylar.txt", + "LISE" , + 1000 , + 1 , + 1 ); + + EnergyLoss protonTargetGaz = EnergyLoss ("proton_D2gaz_1b_26K.txt" , + "LISE" , + 1000 , + 1 , + 1 ); + + EnergyLoss protonStripAl = EnergyLoss ("proton_Al.txt" , + "LISE" , + 100 , + 1 , + 1 ); + + + } + +using namespace ENERGYLOSS ; +// ---------------------------------------------------------------------------------------------- + +///////////////////////////////////////////////////////////////////////////////////////////////// +// -----------------------------------Random Engine---------------------------------------------- +#include "TRandom3.h" +namespace RANDOMENGINE + { + + TRandom3 RandomEngine = TRandom3(); + + } + +using namespace RANDOMENGINE ; +// ---------------------------------------------------------------------------------------------- +///////////////////////////////////////////////////////////////////////////////////////////////// + diff --git a/NPAnalysis/e530/macro/CrossSection.c b/NPAnalysis/e530/macro/CrossSection.c new file mode 100644 index 0000000000000000000000000000000000000000..5640efb446f5464a98882b5a40ea015ef2745171 --- /dev/null +++ b/NPAnalysis/e530/macro/CrossSection.c @@ -0,0 +1,121 @@ +{ + gROOT->Reset(); + gStyle->SetOptStat(1); + gStyle->SetPalette(51,0); + /////////////////////// + /////////////////////// + + double Pi = 3.141592654 ; + +///Usefull myAngleInRad * RadToDeg = myAngle In deg :p + double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180 + double RadToDeg = 180./Pi ; // 360/2Pi = 180/Pi + +TFile *file0 = TFile::Open("./Result/myResult.root"); + + cEA = new TCanvas("cEA","Kinematic Line" ,100,100,900,900); + hEA->Draw("COLZ"); + cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600); + hEx->Draw(); + + cEHexa = new TCanvas("cEHexa","Hexaneutron bound Energy" ,100,100,600,600); + hEHexa->Draw(); + + cCM = new TCanvas("cCm" , "Cross Section (CM)" , 100 , 100 , 900, 900) ; + hThetaCM->Draw(); + + vector<double> relative_error; + double NumberOfHitX=0; + ifstream efficFile; + efficFile.open("./17cm.efficiency"); + string EffString, ThetaString; + double xxx = hThetaLabCM->GetNbinsY(); + + for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++) + { + // Calculate Relative error + NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx ); + if(NumberOfHitX!=0)relative_error.push_back( 1. / ( sqrt(NumberOfHitX) ) ); + else {relative_error.push_back(0) ;} + } + + for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++) + { + //Efficacity correction + efficFile >> ThetaString >> EffString; + double Eff = atoi(EffString.c_str()); + for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++) + { + if (Eff!=0) + { + double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff) ; + hThetaLabCM->SetCellContent( i, j, pipo ); + } + else hThetaLabCM->SetCellContent( i, j, 0 ); + } + } + + TH1D* hCrossSection = new TH1D(); + hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ; + + for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++) + { + // Calculate Error + hCrossSection->SetBinContent(i, + NumberOfHitX=hCrossSection->GetBinContent(i)/(sin(hCrossSection->GetBinCenter(i)*DegToRad) )); + double error=NumberOfHitX*relative_error[i]; + hCrossSection->SetBinError(i,error); + } + + hCrossSection->Draw(); + + string Path="../Data/CrossSection/11Li(d,3He)10He.txt"; + ifstream CrossSectionFile; + CrossSectionFile.open(Path.c_str()); + if( CrossSectionFile.is_open() ) + cout << " file " << Path << " loading " << endl; + else{ + cout << " Error, no file " << Path << " found" << endl; return;} + + + //Normalisation: + //Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ; + Int_t Maximum_Bin = 3 ; + Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ; + Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ; + Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ; + + Double_t Normalisation_Factor = 0 ; + Double_t ThetaVal = 0 ; + Double_t CSVal = 0 ; + TMarker marker ; + string theta,CS ; + + while( !CrossSectionFile.eof() && ThetaVal<60) + { + CrossSectionFile >> theta >> CS ; + + ThetaVal=atof(theta.c_str()) ; + + if(ThetaVal>(Maximum_theta-Bin_Width/1000) && ThetaVal<(Maximum_theta+Bin_Width/1000) ) + { + CSVal=atof(CS.c_str()) ; + Normalisation_Factor = Maximum/CSVal; + } + } + CrossSectionFile.close(); + CrossSectionFile.open(Path.c_str()); + + //Reading Cross Section + + ThetaVal=0; + while( !CrossSectionFile.eof() && ThetaVal<60) + { + + CrossSectionFile >> theta >> CS ; + CSVal=atof(CS.c_str()) ; + ThetaVal=atof(theta.c_str()) ; + marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor ) ; + } + +} diff --git a/NPAnalysis/e530/macro/TimeOfFlight.c b/NPAnalysis/e530/macro/TimeOfFlight.c new file mode 100644 index 0000000000000000000000000000000000000000..916932fc04c540832f14c0c4cf14b4889d8972be --- /dev/null +++ b/NPAnalysis/e530/macro/TimeOfFlight.c @@ -0,0 +1,14 @@ +{ + + gROOT->Reset(); + gStyle->SetOptStat(1); + gStyle->SetPalette(51,0); + //gStyle->SetPalette(1); + /////////////////////// + /////////////////////// +TFile *file0 = TFile::Open("./Result/myResult.root"); + + TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800); + hTOF1234->Draw("COLZ"); + +} diff --git a/NPAnalysis/e530/macro/affich.c b/NPAnalysis/e530/macro/affich.c new file mode 100644 index 0000000000000000000000000000000000000000..2c8537ea62be3b91ec797d490b630dd8f4ecc6eb --- /dev/null +++ b/NPAnalysis/e530/macro/affich.c @@ -0,0 +1,291 @@ +{ + + gROOT->Reset(); + gStyle->SetOptStat(1); + //gStyle->SetPalette(51,0); + gStyle->SetPalette(1); + /////////////////////// + /////////////////////// +TFile *file0 = TFile::Open("./Result/myResult.root"); + +/* TCanvas* cTOF = new TCanvas("cTOF","Time of Flight" ,10,10,800,800); + cTOF->Divide(3,2); + + cTOF->cd(1); + hTOF1->Draw("COLZ"); + + cTOF->cd(2); + hTOF2->Draw("COLZ"); + + cTOF->cd(3); + hTOF3->Draw("COLZ"); + + cTOF->cd(4); + hTOF4->Draw("COLZ"); + + cTOF->cd(5); + hTOF5->Draw("COLZ"); */ + + double Pi = 3.141592654 ; + double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180 + +/* vector<double> relative_error; + double NumberOfHitX=0; + cTheta2D = new TCanvas("cTheta2d","Cross Section" ,100,100,900,900); + hThetaLabCM->Draw("COLZ");*/ + cTheta = new TCanvas("cTheta","Cross Section" ,100,100,900,900); + hTheta->Draw(); +/* ifstream efficFile; + + efficFile.open("./efficiency/20cm2.efficiency"); + string EffString, ThetaString; + double xxx = hThetaLabCM->GetNbinsY(); + + for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++) + { + // Calculate Relative error + NumberOfHitX=hThetaLabCM->Integral(i, i, 1, xxx ); + if(NumberOfHitX!=0)relative_error.push_back( ( sqrt(NumberOfHitX) ) / (NumberOfHitX) ); + else {relative_error.push_back(0) ;} + } + + for(Int_t j = 1 ; j < hThetaLabCM->GetNbinsY() ; j++) + { + //Efficacity correction + efficFile >> ThetaString >> EffString; + double Eff = atoi(EffString.c_str()); + for(Int_t i = 1 ; i < hThetaLabCM->GetNbinsX() ; i++) + { + if (Eff!=0) + { + double pipo = (hThetaLabCM->GetCellContent(i,j))/(Eff/100) ; + hThetaLabCM->SetCellContent( i, j, pipo ); + } + else hThetaLabCM->SetCellContent( i, j, 0 ); + } + } + + + TH1D* hCrossSection = new TH1D(); + hCrossSection = hThetaLabCM->ProjectionX("Cross Section", 0, -1, "") ; + + + for(Int_t i = 1 ; i < hCrossSection->GetNbinsX() ; i++) + { + // Calculate Error + NumberOfHitX=hCrossSection->GetBinContent(i); + double error=NumberOfHitX*relative_error[i]; + hCrossSection->SetBinError(i,error); + } + + hCrossSection->Draw(); + +/* ofstream efficFile; + efficFile.open("20cm.efficiency"); + + for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++) + { + + Double_t Flux = 100000*sin(hTheta->GetBinCenter(i)*DegToRad)*2*DegToRad/2; + + if(Flux!=0) + { + Double_t Efficiency = ( ( Flux - hTheta->GetBinContent(i) ) / Flux ) * 100 ; + + hTheta->SetBinContent(i, 100-Efficiency ) ; + + efficFile << hTheta->GetBinCenter(i) << " " << 100-Efficiency << endl; + } + + else hTheta->SetBinContent(i, 0 ); + } + + + //hTheta->Draw();*/ + + //efficiency correction and error bar calculation +/* ifstream efficFile; + efficFile.open("./efficiency/20cmCM.efficiency"); + string EffString, ThetaString; + double Eff, relative_error, error, toto ; + + for(Int_t i = 1 ; i < hThetaCM->GetNbinsX() ; i++) + { + + efficFile >> ThetaString >> EffString; + Eff = atoi(EffString.c_str()); + if(hThetaCM->GetBinContent(i)!=0)relative_error = sqrt(hThetaCM->GetBinContent(i))/(hThetaCM->GetBinContent(i)); + + toto=hThetaCM->GetBinContent(i)*(100-Eff) ; + hThetaCM->SetBinContent(i, toto) ; + + error=hThetaCM->GetBinContent(i)*relative_error; + hThetaCM->SetBinError(i,error); + }*/ + + + + + //hThetaCM->Draw(); + +/* string Path="/home/Adrien/Desktop/geant/Simulation/Data/CrossSection/11Li(d,3He)10He.txt"; + ifstream CrossSectionFile; + CrossSectionFile.open(Path.c_str()); + if( CrossSectionFile.is_open() ) + cout << " file " << Path << " loading " << endl; + else{ + cout << " Error, no file " << Path << " found" << endl; return;} + + + //Normalisation: + Int_t Maximum_Bin = hCrossSection->GetMaximumBin() ; + Double_t Maximum_theta = hCrossSection->GetBinCenter(Maximum_Bin) ; + Double_t Bin_Width = hCrossSection->GetBinWidth(Maximum_Bin) ; + Double_t Maximum = hCrossSection->GetBinContent(Maximum_Bin) ; + + Double_t Normalisation_Factor=0; + + + + Double_t ThetaVal=0 ; + Double_t CSVal=0 ; + TMarker marker ; + string theta,CS ; + + while( !CrossSectionFile.eof() && ThetaVal<60) + { + CrossSectionFile >> theta >> CS ; + + ThetaVal=atof(theta.c_str()) ; + + if(ThetaVal>(Maximum_theta-Bin_Width/4) && ThetaVal<(Maximum_theta+Bin_Width/4) ) + { + CSVal=atof(CS.c_str()) ; + Normalisation_Factor = Maximum/CSVal; + } + } + CrossSectionFile.close(); + CrossSectionFile.open(Path.c_str()); + + //Reading Cross Section + + ThetaVal=0; + while( !CrossSectionFile.eof() && ThetaVal<60) + { + + CrossSectionFile >> theta >> CS ; + CSVal=atof(CS.c_str()) ; + ThetaVal=atof(theta.c_str()) ; + marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ; + } + + + +/* + //Normalisation: + Double_t Maximum_Bin = hTheta->GetMaximumBin() ; + Double_t Maximum = hTheta->GetBinContent(Maximum_Bin) ; + + Double_t Normalisation_Factor = Maximum/3.; + + + + //Reading Cross Section + Double_t ThetaVal=0 ; + Double_t CSVal=0 ; + TMarker marker ; + string theta,CS ; + + while( !CrossSectionFile.eof() && ThetaVal<60) + { + + CrossSectionFile >> theta >> CS ; + CSVal=atof(CS.c_str()) ; + ThetaVal=atof(theta.c_str()) ; + + ThetaVal=0.000479597*pow(ThetaVal,3)-0.0670771*pow(ThetaVal,2)+3.08149*ThetaVal-0.339958 ; + + marker->DrawMarker(ThetaVal,CSVal*Normalisation_Factor) ; + } +// cTheta->SetLogy(1); + + + + + /* + cKine = new TCanvas("cKine","Kinematics lines" ,10,10,800,800); + hKine->Draw("COLZ");*/ + +/* cKineInit = new TCanvas("cKineInit","Initial Kinematics lines" ,100,100,600,600); + hKineInit->Draw("COLZ"); + + cEDE = new TCanvas("EDE","EDE add Strip" ,100,100,600,600); + hEDEAddStrip->Draw("COLZ"); +/* cG = new TCanvas("cG","Strip position",500,100,800,600); + cG->Divide(3,2); + + cG->cd(1); + Agraph2D->Draw("P0") ; + cG->cd(2); + Bgraph2D->Draw("P0") ; + cG->cd(3); + Cgraph2D->Draw("P0") ; + cG->cd(4); + Dgraph2D->Draw("P0") ; + cG->cd(5); + Egraph2D->Draw("P0") ; + + cG2 = new TCanvas("cG2","all Strip position",500,100,800,600); + TOTgraph2D->Draw("P0") ; + */ + cH = new TCanvas("cH","Hit density",500,100,1000,800); + cH->Divide(3,2); + + cH->cd(1); + hHIT4->Draw("COLZ"); + cH->cd(2); + hHIT2->Draw("COLZ"); + cH->cd(3); + hHIT5->Draw("COLZ"); + cH->cd(4); + hHIT1->Draw("COLZ"); + cH->cd(5); + hHIT3->Draw("COLZ"); + cH->cd(6); + hThetaHeavy->Draw(); + +/* cH->cd(5); + hXZ->Draw("COLZ"); + cH->cd(6); + hXY->Draw("COLZ"); + + cEDE = new TCanvas("cEDE","EDE indentification Spectra",500,100,1000,800); + cEDE->Divide(3,2); + + cEDE->cd(1); + hEDE1->Draw("COLZ"); + cEDE->cd(2); + hEDE2->Draw("COLZ"); + cEDE->cd(3); + hEDE3->Draw("COLZ"); + cEDE->cd(4); + hEDE4->Draw("COLZ"); + cEDE->cd(5); + hEDE5->Draw("COLZ"); + cEDE->cd(6); + TH2F* hEDET= new TH2F("hEDET","",4000,-1,600, 400, -1, 25) ; + hEDET->Add(hEDE1); + hEDET->Add(hEDE2); + hEDET->Add(hEDE3); + hEDET->Add(hEDE4); + hEDET->Draw("CLOZ");*/ + + cEx = new TCanvas("cEx","Excitation Energy" ,100,100,600,600); + hEx->Draw(); + /*cE = new TCanvas("cE","Light Energy" ,500,100,1000,800); + cE->Divide(2); + cE->cd(1); + hE1234->Draw(); + cE->cd(2); + hE5->Draw();*/ +} diff --git a/NPAnalysis/e530/macro/efficiency.c b/NPAnalysis/e530/macro/efficiency.c new file mode 100644 index 0000000000000000000000000000000000000000..9a04af185f75853d0bce5bae2009cd8cd01cae1a --- /dev/null +++ b/NPAnalysis/e530/macro/efficiency.c @@ -0,0 +1,58 @@ +{ + gROOT->Reset(); + gStyle->SetOptStat(1); + gStyle->SetPalette(51,0); + /////////////////////// + /////////////////////// +TFile *file0 = TFile::Open("./Result/myResult.root"); + + double Pi = 3.141592654 ; + double DegToRad = Pi/180. ; // 2Pi/360 = Pi/180 + + + ofstream efficFile; + efficFile.open("12cm.efficiency"); + + for(Int_t i = 1 ; i < hTheta->GetNbinsX() ; i++) + { + Double_t Flux = 1000000./2. * sin(hTheta->GetBinCenter(i)*DegToRad) * (2*DegToRad) ; + + if(Flux!=0) + { + Double_t Efficiency = 100*( (hTheta->GetBinContent(i)) / Flux ); + hTheta->SetBinContent(i, Efficiency ) ; + } + + else hTheta->SetBinContent(i, 0 ); + + efficFile << hTheta->GetBinCenter(i) << " " << Efficiency << endl; + } + + cEff = new TCanvas("cEff","Efficiency" ,100,100,600,600); + hTheta->Draw(); + + cHit = new TCanvas("cHit","Hit" ,100,100,600,600); + hXY->Draw("COLZ"); + + cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600); + cHit2->Divide(3,2); + + cHit2->cd(1); + hXY1->Draw("COLZ"); + + cHit2->cd(2); + hXY2->Draw("COLZ"); + + cHit2->cd(3); + hXY3->Draw("COLZ"); + + cHit2->cd(4); + hXY4->Draw("COLZ"); + + cHit2->cd(5); + hXY5->Draw("COLZ"); + + cHit2->cd(6); + hXY6->Draw("COLZ"); + +} diff --git a/NPAnalysis/e530/macro/hit.c b/NPAnalysis/e530/macro/hit.c new file mode 100644 index 0000000000000000000000000000000000000000..50db19da5e7c54b1cba977ffe34232fb328dcdd0 --- /dev/null +++ b/NPAnalysis/e530/macro/hit.c @@ -0,0 +1,44 @@ +{ + gROOT->Reset(); + gStyle->SetOptStat(1); + gStyle->SetPalette(51,0); + /////////////////////// + /////////////////////// +TFile *file0 = TFile::Open("./Result/myResult.root"); + + cEff = new TCanvas("cEff","Theta Distribution" ,100,100,600,600); + hTheta->Draw(); + + + cCM = new TCanvas("cCm" , "ThetaCM" , 100 , 100 , 600, 600) ; + hThetaCM->Draw(); + + cHit = new TCanvas("cHit","Hit" ,100,100,600,600); + hXY->Draw("COLZ"); + + + cHit2 = new TCanvas("cHit2","Hit" ,100,100,600,600); + cHit2->Divide(3,2); + + cHit2->cd(1); + hXY1->Draw("COLZ"); + + cHit2->cd(2); + hXY2->Draw("COLZ"); + + cHit2->cd(3); + hXY3->Draw("COLZ"); + + cHit2->cd(4); + hXY4->Draw("COLZ"); + + cHit2->cd(5); + hXY5->Draw("COLZ"); + + cHit2->cd(6); + hXY6->Draw("COLZ"); + + cEx = new TCanvas("cEx","Excitation Energy" ,100,100,300,300); + hEx->Draw(); + +} diff --git a/NPAnalysis/e530/pipo.txt b/NPAnalysis/e530/pipo.txt new file mode 100644 index 0000000000000000000000000000000000000000..7b4859f4567953d6a47cd82929be6fa61bda9f41 --- /dev/null +++ b/NPAnalysis/e530/pipo.txt @@ -0,0 +1,2 @@ +CalibrationFilePath + dummy.txt diff --git a/NPAnalysis/e530/src/Analysis.cc b/NPAnalysis/e530/src/Analysis.cc new file mode 100644 index 0000000000000000000000000000000000000000..555c8163e55ee1d72133708a212f64e76133592b --- /dev/null +++ b/NPAnalysis/e530/src/Analysis.cc @@ -0,0 +1,275 @@ +#include "ObjectManager.hh" + +using namespace std; + +int main(int argc,char** argv) +{ + + if(argc!=5) + { + cout << + "you need to specify both a Reaction file and a Detector file such as : "<< endl; + cout << + "Analysis myReaction.reaction myDetector.detector runToRead.run" << endl; + return 0; + } + + string reactionfileName = argv[1] ; + string detectorfileName = argv[2] ; + string calibrationfileName = argv[3] ; + string runToReadfileName = argv[4] ; + + // First of All instantiate RootInput and Output + // Detector will be attached later + RootInput:: getInstance(runToReadfileName) ; + //RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree") ; + RootOutput::getInstance("Analysis/Fe60", "Analysed_Data") ; + + // Instantiate some Reaction + + NPL::Reaction* Fe60Reaction = new Reaction; + Fe60Reaction -> ReadConfigurationFile(reactionfileName); + + // Instantiate the detector using a file + NPA::DetectorManager* myDetector = new DetectorManager ; + myDetector -> ReadConfigurationFile(detectorfileName) ; + + // Instantiate the Calibration Manger using a file + // CalibrationManager* myCalibration = new CalibrationManager(calibrationfileName) ; + CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName); + + // Attach more branch to the output + + double ELab[2], ExcitationEnergy[2] ; + double ThetaLab[2] , ThetaCM[2] ; + double X[2] , Y[2] ; + + // Exitation Energy + RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",ExcitationEnergy,"ExcitationEnergy[2]/D") ; + + //E lab et Theta lab + RootOutput::getInstance()->GetTree()->Branch("ThetaCM",ThetaCM,"ThetaCM[2]/D") ; + RootOutput::getInstance()->GetTree()->Branch("ELab",ELab,"ELab[2]/D") ; + RootOutput::getInstance()->GetTree()->Branch("ThetaLab",ThetaLab,"ThetaLab[2]/D") ; + RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") ; + RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]") ; + + + // Get the formed Chained Tree and Treat it + TChain* Chain = RootInput:: getInstance() -> GetChain() ; + + // Open the ThinSi Branch + // Chain->SetBranchStatus("InitialConditions",true) ; + // Chain->SetBranchStatus("fIC_*",true) ; + + // TInitialConditions* Init = new TInitialConditions(); + // Chain->SetBranchAddress("InitialConditions" ,&Init ); + + double XTarget=0 ; double YTarget=0; double BeamTheta = 0 ; double BeamPhi = 0 ; double E=-1000; + + // Get Detector Pointer: + TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ; + // TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ; + // TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ; + cout << " ///////// Starting Analysis ///////// "<< endl << endl ; + + int i ,N=Chain -> GetEntries(); + + N = 1000; + cout << " Number of Event to be treated : " << N << endl ; + + clock_t begin=clock(); + clock_t end=begin; + for ( i = 0 ; i < N ; i ++ ) + { + // Clear local branch + for(int hh = 0 ; hh <2 ; hh++) + { + ELab[hh] = -1 ; ExcitationEnergy[hh] = -1 ; ThetaLab[hh] = -1 ; + X[hh] = -1000 ; Y[hh] = -1000 ; ThetaCM[hh] = -1 ; + } + + // Minimum code + if( i%10000 == 0 && i!=0) { + cout.precision(5); + end=clock(); + double TimeElapsed = (end-begin)/CLOCKS_PER_SEC; + double percent = (double)i/N ; + double TimeToWait = (TimeElapsed/percent) - TimeElapsed ; + cout << "\r Progression:" << percent*100 + << " % \t | \t Remaining time : ~" + << TimeToWait <<"s"<< flush; + } + + else if (i==N-1) cout << "\r Progression:" + << " 100% " <<endl; + + Chain -> GetEntry(i); + // Clear Previous Event + myDetector -> ClearEventPhysics() ; + // Build the new one + myDetector -> BuildPhysicalEvent() ; + //// + + + // Target (from initial condition) + //XTarget = Init->GetICPositionX(0); + // YTarget = Init->GetICPositionY(0); + // XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1); + // YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1); + // BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; + // BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ; + TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ; + //// + + // Must 2 And ThinSi // + //for(int hit = 0; hit < M2 -> GetEventMultiplicity() ; hit ++) + for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++) + { + ELab[hit] = -1 ; ThetaLab[hit] = -1; + // Get Hit Direction + TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0); + + // Angle between beam and particle + ThetaLab[hit] = ThetaCalculation ( HitDirection , BeamDirection ) ; + // Angle between particule and z axis (target Normal) + double ThetaN = ThetaCalculation ( HitDirection , TVector3(0,0,1) ) ; + // ANgle between particule and Must2 Si surface + double ThetaMM2Surface = ThetaCalculation ( HitDirection , M2 -> GetTelescopeNormal(hit) ); + + if(M2 -> GetPositionOfInteraction(hit).Z() > 0) + { + if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0) + { + ELab[hit] = M2 -> Si_E[hit] ; + + /* + ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 2*0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + */ + // ELab[hit]= He3StripSi.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + // 20*micrometer , // Target Thickness at 0 degree + // ThetaMM2Surface ); + + //if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit]; + + /* + ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 15*micrometer , // Target Thickness at 0 degree + ThetaN ); + + ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 1.5*mm , // Target Thickness at 0 degree + ThetaN ); + */ + + ThetaCM[hit] = Fe60Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ; + ExcitationEnergy[hit] = Fe60Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ; + X[hit] = HitDirection . X(); + Y[hit] = HitDirection . Y(); + ThetaLab[hit] = ThetaLab[hit] / deg ; + } + + /* + else if(M2 ->CsI_E[hit] > 0 && M2 -> Si_E[hit] > 0) + { + + ELab[hit]= M2 ->CsI_E[hit] ; + + ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 3*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + ELab[hit]+= M2 ->Si_E[hit]; + + ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + //if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit]; + + + ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 15*micrometer , // Target Thickness at 0 degree + ThetaN ); + + ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 1.5*mm , // Target Thickness at 0 degree + ThetaN ); + + ThetaCM[hit]= Fe60Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ; + ExcitationEnergy[hit] = Fe60Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ; + X[hit] = HitDirection . X(); + Y[hit] = HitDirection . Y(); + ThetaLab[hit] = ThetaLab[hit] / deg ; + } + */ + else {ExcitationEnergy[hit]=-100 ; X[hit] = -100 ; Y[hit]= -100 ; ThetaLab[hit]=-100;ThetaCM[hit]=-100;} + + } + + /*else if(M2 -> GetPositionOfInteraction(hit).Z()<0) + { + + if(ELab[hit]>-1000 ) + { + if(ELab[hit]>18) + { + ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + } + + ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 0.4*micrometer , // Target Thickness at 0 degree + ThetaMM2Surface ); + + ELab[hit]= protonTargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 15*micrometer , // Target Thickness at 0 degree + ThetaN ); + + ELab[hit]= protonTargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle + 1.5*mm , // Target Thickness at 0 degree + ThetaN ); + ThetaCM[hit] = myReaction -> EnergyLabToThetaCM( ELab[hit] , 1 ) /deg ; + ExcitationEnergy[hit] = myReaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ; + X[hit] = HitDirection . X(); + Y[hit] = HitDirection . Y(); + ThetaLab[hit] = ThetaLab[hit] / deg ; + } + + else {ExcitationEnergy[hit]=-100 ; X[hit] = -100 ; Y[hit] = -100 ;ThetaLab[hit]=-100; ThetaCM[hit]=-100 ;} + + } */ + + + } + + RootOutput::getInstance()->GetTree()->Fill() ; + } + + cout << " A total of " << i << " event has been annalysed " << endl ; + cout << endl << " ///////////////////////////////////// "<< endl<< endl ; + RootOutput::getInstance()->Destroy(); + return 0 ; +} + +double ThetaCalculation (TVector3 A , TVector3 B) +{ + double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ; + return Theta*rad ; +} + diff --git a/NPAnalysis/e530/src/GNUmakefile b/NPAnalysis/e530/src/GNUmakefile new file mode 100644 index 0000000000000000000000000000000000000000..28c404622cec4d40b2e93aa0009fc4fc4fa97d26 --- /dev/null +++ b/NPAnalysis/e530/src/GNUmakefile @@ -0,0 +1,42 @@ +###Make file made by Adrien MATTA/ Institut de Physique Nucleaire d'Orsay IPNO### +# Made to compile the ROOT Analyser for MUST2 experiment + +CPP=g++ +EXEC=Analysis + +# local includes +NPAINCLUDES = ../include + +# ROOT includes +CXXFLAGS += `root-config --cflags` + +# CLHEP includes +CXXFLAGS += -I$(CLHEP_INCLUDE_DIR) +CXXFLAGS += -I$(NPAINCLUDES) +CXXFLAGS += -I$(NPLIB)/include + +LDFLAGS = `root-config --libs` -lMathMore +LDFLAGS+= `$(NPLIB)/liblist` +LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB) + +SRC= $(wildcard *.cc) +INC= $(wildcard $(NPAINCLUDES)/*.hh) +OBJ=$(SRC:.cc=.o) + +#all:$(EXEC) +# @$(CPP) -o $@ -c $< $(CXXFLAGS) + +Analysis:$(OBJ) $(INC) + @$(CPP) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) + mv Analysis ../Analysis + +%.o: %.cc + @$(CPP) $(CXXFLAGS) -o $@ -c $< + +.PHONY: clean mrproper + +clean: + rm -rf *.o + +mrproper: clean + rm -rf $(EXEC)