From 3ca6f6c09b58e6d5bbc4d88149e2571264daf6ff Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Fri, 14 Jan 2011 23:21:31 +0000
Subject: [PATCH] * Fix bug in TW1Physics.cxx + GetPositionOfInteraction()
method is now implemented
* Add analysis support to W1
+ To be tested
---
NPAnalysis/W1/Makefile | 10 +
NPAnalysis/W1/RunToTreat.txt | 7 +
NPAnalysis/W1/include/ObjectManager.hh | 131 ++++++++++++
.../W1/macros/DisplayInputCrossSection.C | 148 ++++++++++++++
NPAnalysis/W1/src/Analysis.cc | 189 ++++++++++++++++++
NPAnalysis/W1/src/GNUmakefile | 42 ++++
NPLib/W1/TW1Physics.cxx | 11 +
7 files changed, 538 insertions(+)
create mode 100644 NPAnalysis/W1/Makefile
create mode 100644 NPAnalysis/W1/RunToTreat.txt
create mode 100644 NPAnalysis/W1/include/ObjectManager.hh
create mode 100644 NPAnalysis/W1/macros/DisplayInputCrossSection.C
create mode 100644 NPAnalysis/W1/src/Analysis.cc
create mode 100644 NPAnalysis/W1/src/GNUmakefile
diff --git a/NPAnalysis/W1/Makefile b/NPAnalysis/W1/Makefile
new file mode 100644
index 000000000..215794b28
--- /dev/null
+++ b/NPAnalysis/W1/Makefile
@@ -0,0 +1,10 @@
+
+Analyse:
+ make -C ./src
+
+clean:
+ make clean -C ./src
+
+distclean:
+ make clean -C ./src
+ rm Analysis
diff --git a/NPAnalysis/W1/RunToTreat.txt b/NPAnalysis/W1/RunToTreat.txt
new file mode 100644
index 000000000..c4154eda3
--- /dev/null
+++ b/NPAnalysis/W1/RunToTreat.txt
@@ -0,0 +1,7 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ ../../Outputs/Simulation/myResult.root
+% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T0_B0_E0_S2mm.root
+% ../../Outputs/Simulation/132Sndp_3p3_10MeVA_T0_B1_E0_S05mm.root
+% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T1_B1_E0_S05mm.root
diff --git a/NPAnalysis/W1/include/ObjectManager.hh b/NPAnalysis/W1/include/ObjectManager.hh
new file mode 100644
index 000000000..7f91f1ce3
--- /dev/null
+++ b/NPAnalysis/W1/include/ObjectManager.hh
@@ -0,0 +1,131 @@
+// 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"
+#include "NPOptionManager.h"
+#include "TW1Physics.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 <TRandom3.h>
+
+// NPL
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TInteractionCoordinates.h"
+#include "TInitialConditions.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 LightTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 1); // LISE++
+ // For 132Sn(d,p)
+ // CD2
+ EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4
+ EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
+ // solid D2
+// EnergyLoss LightTarget = EnergyLoss("proton_D2_solid.G4table", "G4Table", 1000); // G4
+// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_D2_solid.G4table", "G4Table", 1000); // G4
+ // For 134Sn(p,t)
+ // CH2
+// EnergyLoss LightTarget = EnergyLoss("triton_CH2.G4table", "G4Table", 1000); // G4
+// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_CH2.G4table", "G4Table", 1000); // G4
+ // solid H2
+// EnergyLoss LightTarget = EnergyLoss("triton_H2_solid.G4table", "G4Table", 1000); // G4
+// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_H2_solid.G4table", "G4Table", 1000); // G4
+ // For 132Sn(d,t)
+// EnergyLoss LightTarget = EnergyLoss("triton_CD2.G4table", "G4Table", 1000); // G4
+// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
+}
+
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+
diff --git a/NPAnalysis/W1/macros/DisplayInputCrossSection.C b/NPAnalysis/W1/macros/DisplayInputCrossSection.C
new file mode 100644
index 000000000..8f2725bde
--- /dev/null
+++ b/NPAnalysis/W1/macros/DisplayInputCrossSection.C
@@ -0,0 +1,148 @@
+void DisplayInputCrossSection()
+{
+ // Path to cross-section files
+ TString path = gSystem->Getenv("NPTOOL");
+ path += "/Inputs/CrossSection/";
+
+ // Read cross-sections 132Sn(d,p)
+/* TGraph *gr1 = new TGraph(path + "132Sndp_5A_MeV_3p3_ZR_FRC.lis");
+ TGraph *gr2 = new TGraph(path + "132Sndp_5A_MeV_2f7_ZR_FRC.lis");
+ TGraph *gr3 = new TGraph(path + "132Sndp_10A_MeV_3p3_ZR_FRC.lis");
+ TGraph *gr4 = new TGraph(path + "132Sndp_10A_MeV_2f7_ZR_FRC.lis");*/
+ // Jacques
+ TGraph *gr1 = ReadCrossSection("132Sndp_5A_MeV_3p3_ZR_FRC.lis");
+ TGraph *gr2 = ReadCrossSection("132Sndp_5A_MeV_2f7_ZR_FRC.lis");
+ TGraph *gr3 = ReadCrossSection("132Sndp_10A_MeV_3p3_ZR_FRC.lis");
+ TGraph *gr4 = ReadCrossSection("132Sndp_10A_MeV_2f7_ZR_FRC.lis");
+ // Didier
+ TGraph *gr5 = new TGraph(path + "sn132dp_gs_10AMeV.txt");
+
+ // Read cross-section 134Sn(p,t)
+ // Didier
+ TGraph *grpt1 = ReadCrossSectionPT("CS_Ep10MeV_sn134pt_gs_1h9demi.dat");
+ TGraph *grpt2 = ReadCrossSectionPT("CS_Ep15MeV_sn134pt_gs_1h9demi.dat");
+ TGraph *grpt3 = ReadCrossSectionPT("CS_Ep20MeV_sn134pt_gs_1h9demi.dat");
+
+ // Read cross-section 132Sn(d,d)
+ // Angel
+ TGraph *grdd = new TGraph(path + "132Sndd_10A_MeV_ruth.dat");
+ TGraph *grpp = new TGraph(path + "132Snpp_10A_MeV_ruth.dat");
+
+ // Draw cross-sections
+ TCanvas *can = new TCanvas("can");
+ can->SetLogy();
+ can->Draw();
+// TH2F *hframe = new TH2F("hframe", "^{2}H(^{132}Sn,p)^{133}Sn", 180, 0, 180, 100, 1e-2, 100);
+// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-2, 100);
+// TH2F *hframe = new TH2F("hframe", "^{1}H(^{134}Sn,t)^{132}Sn_{g.s.}", 180, 0, 180, 100, 1e-8, 1e-5);
+// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-8, 1e-5);
+ TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-3, 10);
+ hframe->Draw();
+ hframe->SetXTitle("#Theta_{c.m.} [deg]");
+// hframe->SetYTitle("d#sigma/d#Omega [mb/sr]");
+ hframe->SetYTitle("d#sigma/d#Omega / (d#sigma/d#Omega)_{R}");
+// hframe->SetYTitle("#propto d#sigma/d#Omega");
+/* grpt1->SetLineColor(kRed); grpt1->Draw("l");
+ grpt2->SetLineColor(kMagenta); grpt2->Draw("l");
+ grpt3->SetLineColor(kBlue); grpt3->Draw("l");*/
+/* gr1->SetLineColor(kRed); gr1->SetLineStyle(2); gr1->Draw("l");
+ gr2->SetLineColor(kRed); gr2->Draw("l");
+ gr3->SetLineColor(kBlue); gr3->SetLineStyle(2); gr3->Draw("l");
+ gr4->SetLineColor(kBlue); gr4->Draw("l");*/
+// gr5->Draw("l");
+// grdd->Draw("l");
+ grpp->Draw("l");
+
+ // TLegend
+ TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
+// leg->AddEntry(grdd, "10 MeV/u", "l");
+ leg->AddEntry(grpp, "10 MeV/u", "l");
+/* TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
+ leg->AddEntry(grpt1, "1h9/2 10 MeV/u", "l");
+ leg->AddEntry(grpt2, "1h9/2 15 MeV/u", "l");
+ leg->AddEntry(grpt3, "1h9/2 20 MeV/u", "l");*/
+/* TLegend *leg = new TLegend(0.16, 0.17, 0.48, 0.37);
+ leg->AddEntry(gr1, "3p3/2 5 MeV/u", "l");
+ leg->AddEntry(gr2, "2f7/2 5 MeV/u", "l");
+ leg->AddEntry(gr3, "3p3/2 10 MeV/u", "l");
+ leg->AddEntry(gr4, "2f7/2 10 MeV/u", "l");*/
+ leg->SetBorderSize(1);
+ leg->Draw();
+
+/* TMultiGraph *mgr = new TMultiGraph();
+ mgr->Add(gr1, "lp");
+ mgr->Add(gr2, "lp");
+ mgr->Add(gr3, "lp");
+ mgr->Add(gr4, "lp");
+ mgr->Draw("a*");*/
+// gr1->Draw("alp");
+}
+
+
+
+TGraph* ReadCrossSection(const char* fname)
+{
+ // Path to cross-section files
+ TString path = gSystem->Getenv("NPTOOL");
+ path += "/Inputs/CrossSection/";
+
+ // Open file
+ ifstream fich;
+ fich.open(path + fname);
+ if (!fich) cout << "Probleme d'ouverture dans le fichier " << fname << endl;
+
+ // Read file
+ Double_t angle, sigma;
+ Int_t nlines = 0;
+ TGraph *gr = new TGraph();
+ while (fich >> angle >> sigma) {
+ gr->SetPoint(nlines++, angle, sigma * 15); // 15: fm^2 -> mb + D0^2
+ }
+
+ // Close file
+ fich.close();
+
+ // TGraph name
+ gr->SetTitle(fname);
+
+ // test pour savoir si on a bien rempli le TGraph
+ if (gr->GetN() == 0)
+ cout << "Mauvaise lecture du fichier --> probablement mauvais format" << endl;
+
+ return gr;
+}
+
+
+
+TGraph* ReadCrossSectionPT(const char* fname)
+{
+ // Path to cross-section files
+ TString path = gSystem->Getenv("NPTOOL");
+ path += "/Inputs/CrossSection/";
+
+ // Open file
+ ifstream fich;
+ fich.open(path + fname);
+ if (!fich) cout << "Probleme d'ouverture dans le fichier " << fname << endl;
+
+ // Read file
+ Double_t angle, sigma, dum;
+ Int_t nlines = 0;
+ TGraph *gr = new TGraph();
+// while (fich >> angle >> sigma >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum) {
+ while (fich >> angle >> sigma) {
+ gr->SetPoint(nlines++, angle, sigma);
+ }
+
+ // Close file
+ fich.close();
+
+ // TGraph name
+ gr->SetTitle(fname);
+
+ // test pour savoir si on a bien rempli le TGraph
+ if (gr->GetN() == 0)
+ cout << "Mauvaise lecture du fichier --> probablement mauvais format" << endl;
+
+ return gr;
+}
diff --git a/NPAnalysis/W1/src/Analysis.cc b/NPAnalysis/W1/src/Analysis.cc
new file mode 100644
index 000000000..152c3fc40
--- /dev/null
+++ b/NPAnalysis/W1/src/Analysis.cc
@@ -0,0 +1,189 @@
+#include "ObjectManager.hh"
+
+using namespace std;
+
+
+int main(int argc,char** argv)
+{
+ // command line parsing
+ NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+ string reactionfileName = myOptionManager->GetReactionFilePath();
+ string detectorfileName = myOptionManager->GetDetectorFilePath();
+ string calibrationfileName = myOptionManager->GetCalibrationFilePath();
+ string runToReadfileName = myOptionManager->GetRunToReadFilePath();
+ string OutputfileName = myOptionManager->GetOutputFilePath();
+
+ // Instantiate RootInput and RootOutput singleton classes
+ RootInput:: getInstance(runToReadfileName);
+// RootOutput::getInstance("Analysis/"+OutputfileName, "AnalyzedTree") ;
+ RootOutput::getInstance("Analysis/W1_AnalyzedData", "AnalyzedTree");
+
+ // Initialize the reaction
+ NPL::Reaction* myReaction = new Reaction();
+ myReaction->ReadConfigurationFile(reactionfileName);
+
+ // Initialize the detector
+ NPA::DetectorManager* myDetector = new DetectorManager;
+ myDetector->ReadConfigurationFile(detectorfileName);
+
+ // Calculate beam energy at target middle
+ // Target informations
+ cout << "/////////// Target information ///////////" << endl;
+ cout << "Thickness (um): " << myDetector->GetTargetThickness() << endl;
+ // Get nominal beam energy
+ Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV;
+ cout << "Nominal beam energy (MeV): " << BeamEnergyNominal << endl;
+ // Slow beam at target middle
+ Double_t BeamEnergy = BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
+ cout << "Middle-target beam energy (MeV): " << BeamEnergy << endl << endl;
+ // Set energy beam at target middle
+ myReaction->SetBeamEnergy(BeamEnergy);
+
+ // Attach more branch to the output
+ double Ex = 0 ; double ExNoStrips = 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("ExcitationEnergyNoStrips",&ExNoStrips,"ExNoStrips/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 GaspardTracker pointer
+ TW1Physics* W1 = (TW1Physics*) myDetector->m_Detector["W1"];
+
+ // Get the input TChain and treat it
+ TChain* chain = RootInput:: getInstance() -> GetChain();
+
+ // Connect TInitialConditions branch
+ TInitialConditions *initCond = 0;
+ chain->SetBranchAddress("InitialConditions", &initCond);
+ chain->SetBranchStatus("InitialConditions", 1);
+
+ // Connect TInteractionCoordinates branch
+ TInteractionCoordinates *interCoord = 0;
+ chain->SetBranchAddress("InteractionCoordinates", &interCoord);
+ chain->SetBranchStatus("InteractionCoordinates", 0);
+
+ // Analysis is here!
+ int nentries = chain->GetEntries();
+ cout << "/////////// Loop information ///////////" << endl;
+ cout << "Number of entries to be analysed: " << nentries << endl;
+
+ // Default initialization
+ double XTarget = 0;
+ double YTarget = 0;
+ double BeamTheta = 0;
+ double BeamPhi = 0;
+
+ // random generator
+ TRandom3 *gene = new TRandom3();
+
+ // Loop on all events
+ for (int i = 0; i < nentries; i ++) {
+ if (i%10000 == 0 && i!=0) cout << "\r" << i << " analyzed events" << flush;
+ chain -> GetEntry(i);
+
+ // Treat Gaspard event
+ myDetector->ClearEventPhysics();
+ myDetector->BuildPhysicalEvent();
+
+ // Get total energy
+ double E = W1->GetEnergy(0);
+
+ // if there is a hit in the detector array, treat it.
+ double Theta, ThetaStrip, angle, ThetaCM;
+ double DetecX, DetecY, DetecZ;
+ double r;
+ TVector3 A;
+ if (E > -1000) {
+ // Get c.m. angle
+ ThetaCM = initCond->GetICEmittedAngleThetaCM(0) * deg;
+
+ // Get exact scattering angle from TInteractionCoordinates object
+// interCoord->Dump();
+// cout << i << " mult: " << interCoord->GetDetectedMultiplicity() << endl;
+ DetecX = interCoord->GetDetectedPositionX(0);
+ DetecY = interCoord->GetDetectedPositionY(0);
+ DetecZ = interCoord->GetDetectedPositionZ(0);
+ TVector3 Detec(DetecX, DetecY, DetecZ);
+
+ // Get interaction position in detector
+ // This takes into account the strips
+ A = W1->GetPositionOfInteraction(0);
+
+ // Get beam interaction coordinates on target (from initial condition)
+ XTarget = initCond->GetICPositionX(0);
+ YTarget = initCond->GetICPositionY(0);
+ BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
+ BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
+ TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
+
+ // Hit direction taking into account beam position on target
+ TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
+// cout << "A: " << A.X() << " " << A.Y() << " " << A.Z() << endl;
+// cout << "HitDirection: " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
+
+ // Calculate scattering angle w.r.t. optical beam axis (do not take into account beam position on target)
+ ThetaStrip = ThetaCalculation(A, TVector3(0,0,1));
+ Theta = ThetaCalculation(Detec, TVector3(0, 0, 1));
+ // Calculate scattering angle w.r.t. beam (ideal case)
+// ThetaStrip = ThetaCalculation(HitDirection, BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection);
+ // Calculate scattering angle w.r.t. beam (finite spatial resolution)
+// double resol = 800; // in micrometer
+// angle = gene->Rndm() * 2*3.14;
+// r = fabs(gene->Gaus(0, resol)) * micrometer;
+// ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+// Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
+//
+ // Correct for energy loss in the target
+ E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip);
+
+ // Calculate excitation energy
+// if (Theta/deg > 150 && Theta/deg < 180) {
+// if (Theta/deg < 60 && ThetaCM/deg < 90) {
+// if (Theta/deg > 35 && Theta/deg < 45 && E/MeV < 17) {
+// if (Theta/deg < 45) {
+// if (E/MeV < 38) { // for (p,t) reaction
+ if (Theta/deg > 30) { // for (d,p) reaction
+ ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
+ Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
+ }
+ else {
+ Ex = -200;
+ ExNoStrips = -200;
+ }
+ }
+ else {
+ Ex = -100;
+ ExNoStrips = -100;
+ }
+
+ EE = E ; TT = ThetaStrip/deg;
+ if (E>-1000) {
+ X = A . X();
+ Y = A . Y();
+ }
+ else {
+ X = -1000 ; Y = -1000;
+ }
+
+ // Fill output tree
+ RootOutput::getInstance()->GetTree()->Fill();
+ }
+
+ // delete singleton classes
+ RootOutput::getInstance()->Destroy();
+ RootInput::getInstance()->Destroy();
+ NPOptionManager::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/W1/src/GNUmakefile b/NPAnalysis/W1/src/GNUmakefile
new file mode 100644
index 000000000..f6161b9f8
--- /dev/null
+++ b/NPAnalysis/W1/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: Analysis.o $(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)
diff --git a/NPLib/W1/TW1Physics.cxx b/NPLib/W1/TW1Physics.cxx
index 08df9a528..f53322512 100644
--- a/NPLib/W1/TW1Physics.cxx
+++ b/NPLib/W1/TW1Physics.cxx
@@ -445,6 +445,17 @@ void TW1Physics::AddDetector(double theta, double phi, double distance,
+TVector3 TW1Physics::GetPositionOfInteraction(int i)
+{
+ TVector3 Position = TVector3(GetStripPositionX(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]),
+ GetStripPositionY(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]),
+ GetStripPositionZ(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]));
+
+ return Position;
+}
+
+
+
TVector3 TW1Physics::GetDetectorNormal(int i)
{
TVector3 U = TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, 1),
--
GitLab