Skip to content
Snippets Groups Projects
Commit 6b03462f authored by matta's avatar matta
Browse files

* Modified NPToolLogon to avoid error at root startup

parent 3d071d6d
No related branches found
No related tags found
No related merge requests found
......@@ -19,8 +19,8 @@
GeneralTarget
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
Target
THICKNESS= 18.87
RADIUS= 1.25
THICKNESS= 0
RADIUS= 0.1
MATERIAL= CD2
ANGLE= 0
X= 0
......
......@@ -20,14 +20,14 @@ Transfert
ShootHeavy= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ParticleDecay 10He
ParticleDecay 10He
Daughter= 9He n
ExcitationEnergy= 0.5 0
DifferentialCrossSection= flat.txt
shoot= 1 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GammaDecay 9He
GammaDecay 9He
Cascade
BranchingRatio= 30
Energies= 0.1
......@@ -37,7 +37,7 @@ Transfert
Energies= 0.2 0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ParticleDecay 9He
ParticleDecay 9He
Daughter= 8He n
DifferentialCrossSection= flat.txt
shoot= 1 1
#include "Analysis.h"
using namespace std;
#define MM_DEADLAYER 0.4 // Must2 aluminium dead layer in um
double ThetaCalculation (TVector3 A , TVector3 B);
double PhiCalculation (TVector3 A , TVector3 B);
int main(int argc,char** argv)
{
// command line parsing
NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
// Instantiate RootInput
string runToReadfileName = myOptionManager->GetRunToReadFile();
RootInput:: getInstance(runToReadfileName);
// if input files are not given, use those from TAsciiFile
if (myOptionManager->IsDefault("EventGenerator")) {
string name = RootInput::getInstance()->DumpAsciiFile("EventGenerator");
myOptionManager->SetReactionFile(name);
}
if (myOptionManager->IsDefault("DetectorConfiguration")) {
string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
myOptionManager->SetDetectorFile(name);
}
// Instantiate RootOutput
RootOutput::getInstance("Analysis/Must2_AnalysedData", "AnalysedTree");
// get input files from NPOptionManager
......@@ -57,38 +34,13 @@ int main(int argc,char** argv)
myReaction->SetBeamEnergy(BeamEnergy);
myReaction->Print();
// Attach more branch to the output
const int mult = 2;
double ELab[mult], ThetaLab[mult], ThetaCM[mult], PhiLab[mult], ExcitationEnergy[mult];
// some initializations
for (int i = 0; i < mult; i++) {
ELab[i] = -10000; ThetaLab[i] = -10000; ThetaCM[i] = -10000; ExcitationEnergy[i] = -10000; PhiLab[i] = -10000;
}
RootOutput::getInstance()->GetTree()->Branch("ELab", &ELab, "ELab[2]/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab", &ThetaLab, "ThetaLab[2]/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM[2]/D");
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy, "ExcitationEnergy[2]/D");
// Get the formed Chained Tree and Treat it
TChain* Chain = RootInput::getInstance()->GetChain();
// Connect TInitialConditions branch
TInitialConditions *initCond = 0;
Chain->SetBranchAddress("InitialConditions", &initCond);
Chain->SetBranchStatus("InitialConditions", 1);
// Get TMust2Physics pointer
TMust2Physics *M2 = (TMust2Physics*) myDetector -> GetDetector("MUST2") ;
// define user spectra
TH2F* DE_E_protons = new TH2F("DE_E_protons", "DE_E et cut protons", 1000, 0, 25000, 1000, 0, 25000);
TH2F* DE_E = new TH2F("DE_E", "DE_E", 1000, 0, 25000, 1000, 0, 25000);
TH2F* E_Theta = new TH2F("E_Theta", "E_Theta", 60, 100, 180, 500, 0, 20);
TH2F* Position_M2 = new TH2F("Must2Positions", "Must2Positions", 2000, -200, 200, 2000, -200, 200);
TH2F* Position_M2_tot = new TH2F("Must2Positions_tot", "Must2Positions_tot", 130, -1, 129, 130, -1, 130);
cout << endl << "///////// Starting Analysis ///////// "<< endl;
int nentries = Chain->GetEntries();
cout << " Number of Event to be treated : " << nentries << endl;
......@@ -108,72 +60,13 @@ int main(int argc,char** argv)
}
else if (i==nentries-1) cout << "\r Progression:" << " 100% " <<endl;
// Get raw data
// Clear previous event
// Get new raw data
Chain->GetEntry(i);
// Clear Previous Event
myDetector->ClearEventPhysics();
// Build the new event
myDetector->BuildPhysicalEvent();
// Get target interaction position from initial conditions
double XTarget = initCond->GetICPositionX(0);
double YTarget = initCond->GetICPositionY(0);
double ZTarget = initCond->GetICPositionZ(0);
// Calculate beam direction
double BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg; double BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta));
// Loop over hit
for (int hit = 0; hit < M2->Si_E.size(); hit++) {
ExcitationEnergy[hit] = -10000;
ELab[hit] = -10000 ; ThetaLab[hit] = -10000;
// Get hit direction
TVector3 HitDirection = M2->GetPositionOfInteraction(hit) - TVector3(XTarget, YTarget, myDetector->GetTargetZ());
// Angle between beam and particle
ThetaLab[hit] = ThetaCalculation(HitDirection , BeamDirection) * rad;
PhiLab[hit] = (M2->GetPositionOfInteraction(hit)).Phi();
// 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->Si_E[hit] > 0)) {
Position_M2->Fill((M2->GetPositionOfInteraction(hit)).X(), (M2->GetPositionOfInteraction(hit)).Y());
}
if ((M2->Si_E[hit] > 0) && (M2->SiLi_E[hit]) < 0) {
ELab[hit] = M2->Si_E[hit] * MeV;
ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit], MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
ELab[hit] = protonTargetCD2.EvaluateInitialEnergy(ELab[hit], myDetector -> GetTargetThickness() /2 *micrometer , ThetaLab[hit]) /MeV;
ExcitationEnergy[hit] = myReaction->ReconstructRelativistic(ELab[hit]/MeV, ThetaLab[hit]) /MeV;
ThetaLab[hit] = ThetaLab[hit] /deg;
PhiLab[hit] = PhiLab[hit] /deg;
}
else if ((M2->Si_E[hit] > 0) && (M2->SiLi_E[hit]) > 0) {
Position_M2_tot->Fill(M2->Si_X[hit], M2->Si_Y[hit]);
DE_E ->Fill(M2->SiLi_E[hit], M2->Si_E[hit]);
ELab[hit] = M2->SiLi_E[hit]*MeV;
ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit]/MeV, MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
ELab[hit] += M2->Si_E[hit]*MeV;
ELab[hit] = protonStripAl.EvaluateInitialEnergy(ELab[hit]/MeV, MM_DEADLAYER*micrometer, ThetaMM2Surface) /MeV;
ELab[hit] = protonTargetCD2.EvaluateInitialEnergy(ELab[hit]/MeV,myDetector -> GetTargetThickness() /2 *micrometer,ThetaLab[hit]) /MeV;
ExcitationEnergy[hit] = myReaction -> ReconstructRelativistic( ELab[hit]/MeV, ThetaLab[hit]) /MeV ;
ThetaLab[hit] = ThetaLab[hit] /deg;
PhiLab[hit] = PhiLab[hit] /deg;
}
E_Theta->Fill(ThetaLab[hit],ELab[hit]/MeV);
} // end of for loop over hit
RootOutput::getInstance()->GetTree()->Fill();
......@@ -189,26 +82,3 @@ int main(int argc,char** argv)
return 0;
}
double ThetaCalculation (TVector3 A , TVector3 B)
{
double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
return Theta*rad;
}
double PhiCalculation (TVector3 A, TVector3 B)
{
double Theta = acos((A.Dot(B)) / (A.Mag()*B.Mag()));
double R = A.Mag();
double Phi = 0;
double test = A.X() / R * sin(Theta);
Phi = acos(test);
return Phi*rad;
}
......@@ -5,11 +5,6 @@
/////////////////////////////////////////////////////////////////////////////////////////////////
// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
// NPA
#include "DetectorManager.h"
#include "NPOptionManager.h"
#include "TMust2Physics.h"
// STL C++
#include <iostream>
#include <fstream>
......@@ -18,200 +13,37 @@
#include <cmath>
#include <cstdlib>
#include <time.h>
// ROOT
#include <TROOT.h>
#include <TCutG.h>
#include <TChain.h>
#include <TFile.h>
#include <TLeaf.h>
#include <TVector3.h>
#include <TRandom.h>
using namespace std;
// NPL
#include "TPlasticData.h"
#include "NPReaction.h"
#include "RootInput.h"
#include "RootOutput.h"
#include "TInitialConditions.h"
#include "TSharcData.h"
#include "TTigressData.h"
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
// Global variable of the input tree
// ----------------------------------------------------------------------------------------------
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------------------------------------------------------
// Attached what is needed in the input tree
//int tig_event_number;
int tig_num_chan;
int tig_event_id;
#include <TCutG.h>
namespace CUT
{
//open the ROOT file for cuts
// TFile *FCuts = new TFile("cut_protons.root");
// TCutG *cut_protons = (TCutG*) FCuts->Get("protons");
}
int* tig_midas_id;
int* tig_type; //0 for tig10 and 1 for tig 64, sound ok? good.
using namespace CUT ;
// --------------------------------------------------------------------------------------------
int* channel_number;
int* channel_raw;
int* cfd_value;
int* led_value;
int* charge_raw;
float* charge_cal;
////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------ENERGY LOSS----------------------------------------------
#include "NPEnergyLoss.h"
using namespace NPL ;
namespace ENERGYLOSS
{
// Beam Energy Loss
EnergyLoss Fe60TargetCD2 = EnergyLoss ("Fe60[0.0]_CD2.G4table",
"G4Table" ,
1000 );
/*
EnergyLoss Fe60TargetCD2 = EnergyLoss ("60Fe_CD2_new.txt",
"LISE" ,
1000 ,
3 ,
59.934 );
EnergyLoss Fe60TargetCarbone= EnergyLoss ("60Fe_carbone.txt",
"LISE" ,
1000 ,
3 ,
59.934 );
*/
/*
EnergyLoss S34TargetCD2 = EnergyLoss ("S34[0.0]_CD2.G4table",
"G4Table" ,
1000 );
*/
/*
EnergyLoss Si34TargetCD2 = EnergyLoss ("34Si_CD2.txt",
"LISE" ,
1000 ,
1 ,
34 );
EnergyLoss Si34TargetCarbone= EnergyLoss ("34Si_carbone.txt",
"LISE" ,
1000 ,
1 ,
34 );
EnergyLoss S36TargetCD2 = EnergyLoss ("36S_CD2.txt",
"LISE" ,
1000 ,
1 ,
36 );
EnergyLoss S36TargetCarbone= EnergyLoss ("36S_carbone.txt",
"LISE" ,
1000 ,
1 ,
36 );
*/
// proton Energy Loss
EnergyLoss protonTargetCD2 = EnergyLoss ("proton_CD2.G4table",
"G4Table" ,
1000 );
EnergyLoss protonTargetWind = EnergyLoss ("proton_Mylar.txt",
"LISE" ,
1000 ,
3,
1.008 );
EnergyLoss protonTargetC = EnergyLoss ("proton_carbone.txt" ,
"LISE" ,
1000 ,
3 ,
1.008 );
/*
EnergyLoss protonTargetCD2 = EnergyLoss ("proton_cd2.txt" ,
"LISE" ,
1000 ,
3 ,
1.008 );
*/
EnergyLoss protonStripAl = EnergyLoss ("proton_Al.txt" ,
"LISE" ,
1000 ,
3 ,
1.008 );
EnergyLoss protonPadSi = EnergyLoss ("proton_Si.txt",
"LISE" ,
1000 ,
3 ,
1.008 );
}
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------Random Engine----------------------------------------------
#include "TRandom3.h"
namespace RANDOMENGINE
{
TRandom3 RandomEngine = TRandom3();
}
using namespace RANDOMENGINE ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
int* timestamp_low;
int* timestamp_high;
int* timestamp_live;
int* timestamp_tr; // triggers requested
in*t timestamp_ta; // triggers accepted
//
// File generated by rootcint at Sat Nov 17 14:34:33 2012
// File generated by rootcint at Sat Nov 17 17:59:19 2012
// Do NOT change. Changes will be lost next time file is generated
//
......
......@@ -55,20 +55,40 @@ void NPToolLogon(bool verbosemode = false)
// Test if the lib directory is empty or not
if (listfile->GetEntries() > 2) gSystem->Load(libpath+"/libVDetector.so");
// Since the libMust2Physics.so library uses TVector2
// objects, the libPhysics.so ROOT library is loaded.
gSystem->Load("libPhysics.so");
// Loop on all libraries
gSystem->Load("libPhysics.so"); // Needed by Must2 and Sharc
gSystem->Load("libHist.so"); // Needed by TSpectra Class
gSystem->Load("libCore.so"); // Need by Maya
// Loop on Data libraries
Int_t i = 0;
while (listfile->At(i)) {
TString libname = listfile->At(i++)->GetName();
if (libname.Contains("so") && !libname.Contains("libVDetector.so")) {
if (libname.Contains("so") && libname.Contains("Data") && !libname.Contains("libVDetector.so")) {
TString lib = libpath + "/" + libname;
gSystem->Load(lib);
}
}
// Loop on Physics Library
i = 0;
while (listfile->At(i)) {
TString libname = listfile->At(i++)->GetName();
if (libname.Contains("so") && libname.Contains("Physics") &&!libname.Contains("libVDetector.so")) {
TString lib = libpath + "/" + libname;
gSystem->Load(lib);
}
}
// Loop on the Reset of the Library
i = 0;
while (listfile->At(i)) {
TString libname = listfile->At(i++)->GetName();
if (libname.Contains("so") && !libname.Contains("Physics") && !libname.Contains("Data") &&!libname.Contains("libVDetector.so")) {
TString lib = libpath + "/" + libname;
gSystem->Load(lib);
}
}
gROOT->ProcessLine(".L $NPTOOL/NPLib/include/RootInput.h+");
// Since the libdir.GetListOfFiles() commands cds to the
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment