diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis new file mode 100755 index 0000000000000000000000000000000000000000..e4e550e76ef53467904d0205641e15791e9579a0 Binary files /dev/null and b/NPAnalysis/10He_Riken/Analysis differ diff --git a/NPAnalysis/10He_Riken/Makefile b/NPAnalysis/10He_Riken/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..f5bfeb69ffd25b28c9ccb8a81eebb1b66664dfa2 --- /dev/null +++ b/NPAnalysis/10He_Riken/Makefile @@ -0,0 +1,6 @@ + +Analyse: + make -C ./src + +clean: + make clean -C ./src diff --git a/NPAnalysis/10He_Riken/RunToTreat.txt b/NPAnalysis/10He_Riken/RunToTreat.txt new file mode 100644 index 0000000000000000000000000000000000000000..008354fd1b9fa2d1c15477a463be1e2bf76a6b7e --- /dev/null +++ b/NPAnalysis/10He_Riken/RunToTreat.txt @@ -0,0 +1,4 @@ +TTreeName + SimulatedTree +RootFileName + ../../Outputs/Simulation/mySimul.root diff --git a/NPAnalysis/10He_Riken/include/DetectorManager.hh b/NPAnalysis/10He_Riken/include/DetectorManager.hh new file mode 100644 index 0000000000000000000000000000000000000000..eb49f01ef92509cd0531f4f12c7eace833047e61 --- /dev/null +++ b/NPAnalysis/10He_Riken/include/DetectorManager.hh @@ -0,0 +1,39 @@ +#ifndef DetectorManager_h +#define DetectorManager_h + +// NPL +#include "VDetector.h" + +// STL +#include <string> +#include <map> + +using namespace std ; +using namespace NPA ; + +// This class manage a map of virtual detector + +class DetectorManager + { + public: + DetectorManager() ; + ~DetectorManager() ; + + public: + // Read stream at Path and pick-up Token declaration of Detector + void ReadConfigurationFile(string Path) ; + void BuildPhysicalEvent() ; + void BuildSimplePhysicalEvent() ; + void InitializeRootInput() ; + void InitializeRootOutput() ; + void AddDetector(string,VDetector*) ; + void ClearEventPhysics() ; + void ClearEventData() ; + + public: // The map containning all detectors + // Using a Map one can access to any detector using its name + map<string,VDetector*> m_Detector ; + + }; + +#endif diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh new file mode 100644 index 0000000000000000000000000000000000000000..e0a66004be5b8c354e8e3f7175613fbd2c625d7a --- /dev/null +++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh @@ -0,0 +1,117 @@ +// 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.hh" +#include "Must2Array.h" + +// STL C++ +#include <iostream> +#include <fstream> +#include <sstream> +#include <string> +#include <cmath> +#include <cstdlib> + +// ROOT +#include <TROOT.h> +#include <TChain.h> +#include <TFile.h> +#include <TLeaf.h> +#include <TVector3.h> +#include <TRandom.h> + +// NPL +#include "TMust2Data.h" +#include "TMust2Physics.h" +#include "NPReaction.h" +#include "RootInput.h" +#include "RootOutput.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 + { + + // Declare your Energy loss here : + /* EnergyLoss ProtonTarget = EnergyLoss ( "CD2.txt" , + 100 , + 1 ); + */ + } + +using namespace ENERGYLOSS ; +// ---------------------------------------------------------------------------------------------- +///////////////////////////////////////////////////////////////////////////////////////////////// + + diff --git a/NPAnalysis/10He_Riken/macro/CrossSection.c b/NPAnalysis/10He_Riken/macro/CrossSection.c new file mode 100644 index 0000000000000000000000000000000000000000..5640efb446f5464a98882b5a40ea015ef2745171 --- /dev/null +++ b/NPAnalysis/10He_Riken/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/10He_Riken/macro/TimeOfFlight.c b/NPAnalysis/10He_Riken/macro/TimeOfFlight.c new file mode 100644 index 0000000000000000000000000000000000000000..916932fc04c540832f14c0c4cf14b4889d8972be --- /dev/null +++ b/NPAnalysis/10He_Riken/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/10He_Riken/macro/affich.c b/NPAnalysis/10He_Riken/macro/affich.c new file mode 100644 index 0000000000000000000000000000000000000000..2c8537ea62be3b91ec797d490b630dd8f4ecc6eb --- /dev/null +++ b/NPAnalysis/10He_Riken/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/10He_Riken/macro/efficiency.c b/NPAnalysis/10He_Riken/macro/efficiency.c new file mode 100644 index 0000000000000000000000000000000000000000..9a04af185f75853d0bce5bae2009cd8cd01cae1a --- /dev/null +++ b/NPAnalysis/10He_Riken/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/10He_Riken/macro/hit.c b/NPAnalysis/10He_Riken/macro/hit.c new file mode 100644 index 0000000000000000000000000000000000000000..50db19da5e7c54b1cba977ffe34232fb328dcdd0 --- /dev/null +++ b/NPAnalysis/10He_Riken/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/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc new file mode 100644 index 0000000000000000000000000000000000000000..c06a0e94a6eb74807be7e8d399d5dc8fe30e9150 --- /dev/null +++ b/NPAnalysis/10He_Riken/src/Analysis.cc @@ -0,0 +1,80 @@ +#include "ObjectManager.hh" + +using namespace std; + +int main(int argc,char** argv) +{ + + if(argc!=4) + { + cout << + "you need to specify both a Reaction file and a Detector file such as : Analysis myReaction.reaction myDetector.detector runToRead.run" + << endl; + return 0; + } + + string reactionfileName = argv[1] ; + string detectorfileName = argv[2] ; + string runToReadfileName = argv[3] ; + + // First of All instantiate RootInput and Output + // Detector will be attached later + RootInput:: getInstance(runToReadfileName) ; + RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree") ; + + // Instantiate a Reaction + NPL::Reaction* myReaction = new Reaction ; + myReaction -> ReadConfigurationFile(reactionfileName) ; + + // Instantiate the detector using a file + DetectorManager* myDetector = new DetectorManager ; + myDetector -> ReadConfigurationFile(detectorfileName) ; + + // Attach more branch to the output + double Ex = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ; + RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ; + RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ; + RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ; + RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ; + RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ; + + // Get Must2 Pointer: + MUST2Array* M2 = (MUST2Array*) myDetector -> m_Detector["MUST2"] ; + // Get the formed Chained Tree and Treat it + TChain* Chain = RootInput:: getInstance() -> GetChain() ; + int i; + for ( i = 0 ; i < Chain -> GetEntries() ; i ++ ) + { + if( i%10000 == 0 && i!=0) cout << i << " Event annalysed " << endl ; + Chain -> GetEntry(i); + + myDetector -> ClearEventPhysics() ; + myDetector -> BuildPhysicalEvent() ; + + + double E = M2 -> GetEnergyDeposit(); + TVector3 A = M2 -> GetPositionOfInteraction(); + double Theta = ThetaCalculation ( A , TVector3(0,0,1) ) ; + if(E>-1000) Ex = myReaction -> ReconstructRelativistic( E , Theta ) ; + else Ex = -100 ; + EE = E ; TT = Theta/deg ; + if(E>-1000){ + X = A . X(); + Y = A . Y();} + else{X = -1000 ; Y = -1000;} + + RootOutput::getInstance()->GetTree()->Fill() ; + } + cout << "A total of " << i << " event has been annalysed " << endl ; + + RootOutput::getInstance()->Destroy(); + return 0 ; +} + + +double ThetaCalculation (TVector3 A , TVector3 B) + { + double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ; + return Theta ; + } + diff --git a/NPAnalysis/10He_Riken/src/DetectorManager.cc b/NPAnalysis/10He_Riken/src/DetectorManager.cc new file mode 100644 index 0000000000000000000000000000000000000000..000d61fe3827077359d77aa4c4f2af460364fbfd --- /dev/null +++ b/NPAnalysis/10He_Riken/src/DetectorManager.cc @@ -0,0 +1,228 @@ +#include "DetectorManager.hh" + +// STL +#include <iostream> +#include <fstream> +#include <cstdlib> + +// Detector +#include "Must2Array.h" + +///////////////////////////////////////////////////////////////////////////////////////////////// +// Default Constructor and Destructor +DetectorManager::DetectorManager() + {} + +///////////////////////////////////////////////////////////////////////////////////////////////// +DetectorManager::~DetectorManager() + {} + + +///////////////////////////////////////////////////////////////////////////////////////////////// +// Read stream at ConfigFile and pick-up Token declaration of Detector +void DetectorManager::ReadConfigurationFile(string Path) + { + ////////General Reading needs//////// + string LineBuffer; + string DataBuffer; + + /////////Boolean//////////////////// + bool MUST2 = false ; + bool AddThinSi = false ; + bool GeneralTarget = false ; + bool GPDTracker = false ; // Gaspard Tracker + ////////////////////////////////////////////////////////////////////////////////////////// + // added by Nicolas [07/05/09] + string GlobalPath = getenv("NPTOOL"); + Path = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; + ifstream ConfigFile; + ConfigFile.open(Path.c_str()); + + if (ConfigFile.is_open()) + { + cout << "/////////////////////////////" << endl; + cout << " Configuration file " << Path << " loading " << endl; + } + + else + { + cout << " Error, no configuration file" << Path << " found" << endl; + return; + } + + + while (!ConfigFile.eof()) { + //Pick-up next line + getline(ConfigFile, LineBuffer); + //Search for comment Symbol: % + if (LineBuffer.compare(0, 1, "%") == 0) { /*Do Nothing*/ ;} + + /* //////////////////////////////////////////// + //////////// Search for Gaspard //////////// + //////////////////////////////////////////// + else if (LineBuffer.compare(0, 14, "GaspardTracker") == 0 && GPDTracker == false) { + GPDTracker = true ; + cout << "//////// Gaspard Tracker ////////" << endl ; + + // Instantiate the new array as a VDetector Object + VDetector* myDetector = new GaspardTracker() ; + + // Read Position of Telescope + ConfigFile.close() ; + myDetector->ReadConfiguration(Path) ; + ConfigFile.open(Path.c_str()) ; + + // Add array to the VDetector Vector + AddDetector(myDetector) ; + }*/ + + //////////////////////////////////////////// + //////// Search for MUST2 Array //////// + //////////////////////////////////////////// + else if (LineBuffer.compare(0, 10, "MUST2Array") == 0 && MUST2 == false) { + MUST2 = true ; + cout << "//////// MUST2 Array ////////" << endl << endl ; + + // Instantiate the new array as a VDetector Object + VDetector* myDetector = new MUST2Array() ; + + // Read Position of Telescope + ConfigFile.close() ; + myDetector->ReadConfiguration(Path) ; + ConfigFile.open(Path.c_str()) ; + + // Add array to the VDetector Vector + AddDetector("MUST2" , myDetector) ; + } + + /* //////////////////////////////////////////// + ////////// Search for Add.ThinSi /////////// + //////////////////////////////////////////// + else if (LineBuffer.compare(0, 9, "AddThinSi") == 0 && AddThinSi == false) { + AddThinSi = true ; + cout << "//////// Thin Si ////////" << endl << endl ; + + // Instantiate the new array as a VDetector Object + VDetector* myDetector = new ThinSi() ; + + // Read Position of Telescope + ConfigFile.close() ; + myDetector->ReadConfiguration(Path) ; + ConfigFile.open(Path.c_str()) ; + + // Add array to the VDetector Vector + AddDetector(myDetector) ; + } + + //////////////////////////////////////////// + //////////// Search for Target ///////////// + //////////////////////////////////////////// + + else if (LineBuffer.compare(0, 13, "GeneralTarget") == 0 && GeneralTarget == false) { + GeneralTarget = true ; + cout << "////////// Target ///////////" << endl << endl ; + + // Instantiate the new array as a VDetector Objects + VDetector* myDetector = new Target() ; + + // Read Position and target specification + ConfigFile.close() ; + myDetector->ReadConfiguration(Path) ; + ConfigFile.open(Path.c_str()) ; + + m_TargetThickness = ((Target*)myDetector)->GetTargetThickness() ; + m_TargetRadius = ((Target*)myDetector)->GetTargetRadius() ; + m_TargetX = ((Target*)myDetector)->GetTargetX() ; + m_TargetY = ((Target*)myDetector)->GetTargetY() ; + m_TargetZ = ((Target*)myDetector)->GetTargetZ() ; + + // Add target to the VDetector Vector + AddDetector(myDetector) ; + }*/ + + //Nothing understandable + //else ; + } + + ConfigFile.close(); + + InitializeRootInput(); + InitializeRootOutput(); + + return ; + } + +///////////////////////////////////////////////////////////////////////////////////////////////// + +void DetectorManager::BuildPhysicalEvent() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->BuildPhysicalEvent() ; + } + } +///////////////////////////////////////////////////////////////////////////////////////////////// + +void DetectorManager::BuildSimplePhysicalEvent() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->BuildSimplePhysicalEvent() ; + } + } +///////////////////////////////////////////////////////////////////////////////////////////////// + +void DetectorManager::InitializeRootInput() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->InitializeRootInput() ; + } + } +///////////////////////////////////////////////////////////////////////////////////////////////// + +void DetectorManager::InitializeRootOutput() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->InitializeRootOutput() ; + } + + } +///////////////////////////////////////////////////////////////////////////////////////////////// + +void DetectorManager::AddDetector(string DetectorName , VDetector* newDetector) + { + m_Detector["MUST2"] = newDetector ; + } +///////////////////////////////////////////////////////////////////////////////////////////////// +void DetectorManager::ClearEventPhysics() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->ClearEventPhysics() ; + } + + } +///////////////////////////////////////////////////////////////////////////////////////////////// +void DetectorManager::ClearEventData() + { + map<string,VDetector*>::iterator it ; + + for( it = m_Detector.begin() ; it != m_Detector.end() ; ++it) + { + it->second->ClearEventData() ; + } + + } + diff --git a/NPAnalysis/10He_Riken/src/GNUmakefile b/NPAnalysis/10He_Riken/src/GNUmakefile new file mode 100644 index 0000000000000000000000000000000000000000..742a90a491df079d882318b4948c91f27adfe056 --- /dev/null +++ b/NPAnalysis/10He_Riken/src/GNUmakefile @@ -0,0 +1,44 @@ +###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+= -L$(NPLIB)/lib -lVDetector -lIORoot -lReaction -lEnergyLoss \ + -lMust2Data -lMust2Physics \ + -lAnnularS1Data -lGaspardData \ + -lInitialConditions -lInteractionCoordinates +LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB) + +SRC= $(wildcard *.cc) +OBJ=$(SRC:.cc=.o) + +#all:$(EXEC) +# @$(CPP) -o $@ -c $< $(CXXFLAGS) + +Analysis:$(OBJ) + @$(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) diff --git a/NPSimulation/src/ThinSi.cc b/NPSimulation/src/ThinSi.cc index a6680545a43a6acea4a4a5b6091e31d55ec80670..ff0c030aa3c7a38cb0a50375527044c4f3a688ee 100644 --- a/NPSimulation/src/ThinSi.cc +++ b/NPSimulation/src/ThinSi.cc @@ -44,15 +44,15 @@ namespace THINSI const G4double ResoEnergy = 0.042 ;// = 100keV of Resolution // Unit is MeV/2.35 // Geometry - const G4double DetectorSize = 70*mm ; + const G4double DetectorSize = 70*mm ; const G4double SiliconThickness = 20*micrometer ; - const G4double SiliconSize = 50*mm ; + const G4double SiliconSize = 50*mm ; const G4double AluThickness = 0.4*micrometer ; - const G4int NumberOfStrip = 32 ; + const G4int NumberOfStrip = 32 ; - const G4double AluStripFront_PosZ = -0.5*SiliconThickness - 0.5*AluThickness ; - const G4double Si_PosZ = 0 ; - const G4double AluStripBack_PosZ = 0.5*SiliconThickness + 0.5*AluThickness ; + const G4double AluStripFront_PosZ = -0.5*SiliconThickness - 0.5*AluThickness ; + const G4double Si_PosZ = 0 ; + const G4double AluStripBack_PosZ = 0.5*SiliconThickness + 0.5*AluThickness ; const G4double Si_PosX_Shift = 4*mm ; const G4double Si_PosY_Shift = 2*mm ;