diff --git a/NPLib/Detectors/GeTAMU/CMakeLists.txt b/NPLib/Detectors/GeTAMU/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f903e7c54184130b4847ddf163ef18529e89dd84
--- /dev/null
+++ b/NPLib/Detectors/GeTAMU/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_custom_command(OUTPUT TGeTAMUPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TGeTAMUPhysics.h TGeTAMUPhysicsDict.cxx TGeTAMUPhysics.rootmap libNPGeTAMU.dylib DEPENDS TGeTAMUPhysics.h)
+add_custom_command(OUTPUT TGeTAMUDataDict.cxx COMMAND ../../scripts/build_dict.sh TGeTAMUData.h TGeTAMUDataDict.cxx TGeTAMUData.rootmap libNPGeTAMU.dylib DEPENDS TGeTAMUData.h)
+add_library(NPGeTAMU SHARED TGeTAMUData.cxx TGeTAMUPhysics.cxx TGeTAMUDataDict.cxx TGeTAMUPhysicsDict.cxx )
+target_link_libraries(NPGeTAMU ${ROOT_LIBRARIES} NPCore)
+install(FILES TGeTAMUData.h TGeTAMUPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..16915a316343261f6358a30f448b8a75e07d37e0
--- /dev/null
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUData.cxx
@@ -0,0 +1,68 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold the GeTAMU raw data (Made for TIG10 card) *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include <iostream>
+#include <fstream>
+#include <sstream>
+using namespace std;
+
+#include "TGeTAMUData.h"
+
+ClassImp(TGeTAMUData)
+
+/////////////////////////
+TGeTAMUData::TGeTAMUData(){
+}
+
+/////////////////////////
+TGeTAMUData::~TGeTAMUData(){
+}
+
+/////////////////////////
+void TGeTAMUData::Clear(){
+ fTIG_Ge_CloverNbr.clear();
+ fTIG_Ge_CrystalNbr.clear();
+ fTIG_Ge_SegmentNbr.clear();
+ fTIG_Ge_Energy.clear();
+ fTIG_Ge_TimeCFD.clear();
+ fTIG_Ge_TimeLED.clear();
+
+ fTIG_BGO_CloverNbr.clear();
+ fTIG_BGO_CrystalNbr.clear();
+ fTIG_BGO_PmNbr.clear();
+ fTIG_BGO_Energy.clear();
+ fTIG_BGO_TimeCFD.clear();
+ fTIG_BGO_TimeLED.clear();
+}
+
+/////////////////////////
+void TGeTAMUData::Dump() const{
+ // Energy
+ // cout << "GeTAMU_Mult = " << fTIG_CloverNbr.size() << endl;
+
+ // Front
+ // for (UShort_t i = 0; i < fTIG_CloverNbr.size(); i++){
+ // cout << "Clover: " << fTIG_CloverNbr[i]
+ // << " Crystal: " << fTIG_CrystalNbr[i]
+ // << " Energy: " << fTIG_Energy[i]
+ // << " Time: " << fTIG_Time[i] << endl;
+ // }
+}
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUData.h b/NPLib/Detectors/GeTAMU/TGeTAMUData.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3bd30ad9baac1d6038efb6084cebdcf65ccde01
--- /dev/null
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUData.h
@@ -0,0 +1,114 @@
+#ifndef __GeTAMUDATA__
+#define __GeTAMUDATA__
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold the GeTAMU raw data (Made for TIG10 card) *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// STL
+#include<stdlib.h>
+#include <vector>
+#include <map>
+using namespace std ;
+
+// ROOT
+#include "TObject.h"
+
+class TGeTAMUData : public TObject {
+private:
+ // GeTAMU
+ // Energy
+ vector<UShort_t> fTIG_Ge_CloverNbr;
+ vector<UShort_t> fTIG_Ge_CrystalNbr;
+ vector<UShort_t> fTIG_Ge_SegmentNbr;
+ vector<Double_t> fTIG_Ge_Energy;
+ vector<Double_t> fTIG_Ge_TimeCFD;
+ vector<Double_t> fTIG_Ge_TimeLED;
+
+ vector<UShort_t> fTIG_BGO_CloverNbr;
+ vector<UShort_t> fTIG_BGO_CrystalNbr;
+ vector<UShort_t> fTIG_BGO_PmNbr;
+ vector<Double_t> fTIG_BGO_Energy;
+ vector<Double_t> fTIG_BGO_TimeCFD;
+ vector<Double_t> fTIG_BGO_TimeLED;
+
+
+public:
+ TGeTAMUData();
+ virtual ~TGeTAMUData();
+
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+ ///////////////////// SETTERS ////////////////////////
+ inline void SetGeCloverNbr(const UShort_t &GeCloverNbr){fTIG_Ge_CloverNbr.push_back(GeCloverNbr); }
+ inline void SetGeCrystalNbr(const UShort_t &GeCrystalNbr){fTIG_Ge_CrystalNbr.push_back(GeCrystalNbr);}
+ inline void SetGeSegmentNbr(const UShort_t &GeSegmentNbr){fTIG_Ge_SegmentNbr.push_back(GeSegmentNbr);}
+ inline void SetGeEnergy(const Double_t &GeEnergy){fTIG_Ge_Energy.push_back(GeEnergy);}
+ inline void SetGeTimeCFD(const Double_t &GeTimeCFD){fTIG_Ge_TimeCFD.push_back(GeTimeCFD);}
+ inline void SetGeTimeLED(const Double_t &GeTimeLED){fTIG_Ge_TimeLED.push_back(GeTimeLED);}
+
+ inline void SetBGOCloverNbr(const UShort_t &BGOCloverNbr){fTIG_BGO_CloverNbr.push_back(BGOCloverNbr); }
+ inline void SetBGOCrystalNbr(const UShort_t &BGOCrystalNbr){fTIG_BGO_CrystalNbr.push_back(BGOCrystalNbr);}
+ inline void SetBGOPmNbr(const UShort_t &BGOPmNbr){fTIG_BGO_PmNbr.push_back(BGOPmNbr);}
+ inline void SetBGOEnergy(const Double_t &BGOEnergy){fTIG_BGO_Energy.push_back(BGOEnergy);}
+ inline void SetBGOTimeCFD(const Double_t &BGOTimeCFD){fTIG_BGO_TimeCFD.push_back(BGOTimeCFD);}
+ inline void SetBGOTimeLED(const Double_t &BGOTimeLED){fTIG_BGO_TimeLED.push_back(BGOTimeLED);}
+
+ ///////////////////// GETTERS ////////////////////////
+ inline UShort_t GetGeCloverNbr(const unsigned int &i) {return fTIG_Ge_CloverNbr[i]; }
+ inline UShort_t GetGeCrystalNbr(const unsigned int &i) {return fTIG_Ge_CrystalNbr[i]; }
+ inline UShort_t GetGeSegmentNbr(const unsigned int &i) {return fTIG_Ge_SegmentNbr[i]; }
+
+ inline Double_t GetGeEnergy(const unsigned int &i) {return fTIG_Ge_Energy[i];}
+ inline Double_t GetGeTimeCFD(const unsigned int &i) {return fTIG_Ge_TimeCFD[i];}
+ inline Double_t GetGeTimeLED(const unsigned int &i) {return fTIG_Ge_TimeLED[i];}
+
+ inline UShort_t GetBGOCloverNbr(const unsigned int &i) {return fTIG_BGO_CloverNbr[i]; }
+ inline UShort_t GetBGOCrystalNbr(const unsigned int &i) {return fTIG_BGO_CrystalNbr[i]; }
+ inline UShort_t GetBGOPmNbr(const unsigned int &i) {return fTIG_BGO_PmNbr[i]; }
+ inline Double_t GetBGOEnergy(const unsigned int &i) {return fTIG_BGO_Energy[i];}
+ inline Double_t GetBGOTimeCFD(const unsigned int &i) {return fTIG_BGO_TimeCFD[i];}
+ inline Double_t GetBGOTimeLED(const unsigned int &i) {return fTIG_BGO_TimeLED[i];}
+
+ inline unsigned int GetMultiplicityGe() {return fTIG_Ge_CloverNbr.size();}
+ inline unsigned int GetMultiplicityBGO() {return fTIG_BGO_CloverNbr.size();}
+
+ ClassDef(TGeTAMUData,1) // GeTAMUData structure
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e6829fbc5399c9dc2bff4dbddb9af9938850b808
--- /dev/null
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.cxx
@@ -0,0 +1,416 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold GeTAMU treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+using namespace std;
+
+#include "TGeTAMUPhysics.h"
+
+// NPL
+#include "RootInput.h"
+#include "NPDetectorFactory.h"
+#include "RootOutput.h"
+#include "TAsciiFile.h"
+#include "NPSystemOfUnits.h"
+using namespace NPUNITS;
+
+// ROOT
+#include "TChain.h"
+///////////////////////////////////////////////////////////////////////////
+
+ClassImp(TGeTAMUPhysics)
+ ///////////////////////////////////////////////////////////////////////////
+ TGeTAMUPhysics::TGeTAMUPhysics() {
+ m_EventData = new TGeTAMUData ;
+ m_PreTreatedData = new TGeTAMUData ;
+ m_EventPhysics = this ;
+ };
+
+
+/////////////////////////////////////////////////
+void TGeTAMUPhysics::BuildPhysicalEvent(){
+ PreTreat();
+ // Addback Map
+ unsigned int mysize = Gamma_Energy.size();
+ for(unsigned int i = 0 ; i < 16 ; i ++) {
+ for(unsigned int g = 0 ; g < mysize ; g++){
+ if(Clover_Number[g] == i+1 && Segment_Number[g]==0){
+ m_map_E[i] += Gamma_Energy[g];
+ if( Gamma_Energy[g]> m_map_Core_MaxE[i] ){
+ m_map_Core_MaxE[i] = Gamma_Energy[g];
+ m_map_Core_Crystal[i] = Crystal_Number[g];
+ }
+ }
+ if(Clover_Number[g] == i+1 && Segment_Number[g]>0 && Segment_Number[g]<9){
+ if( Gamma_Energy[g]>m_map_Segment_MaxE[i]){
+ m_map_Segment_MaxE[i] = Gamma_Energy[g];
+ m_map_Segment_Crystal[i] = Crystal_Number[g];
+ m_map_Segment[i] = Segment_Number[g];
+ }
+ }
+ }
+ }
+
+ // Final Addback and Doppler Correction
+ int zero = 0;
+ for(int i = 0 ; i < 16 ; i++) {
+ if(m_map_E.find(i)!=m_map_E.end()){
+ int clover = i+1;
+ TVector3 Pos;
+ if(m_map_Segment_MaxE[i]>0)
+ Pos = GetSegmentPosition(clover,m_map_Segment_Crystal[i],m_map_Segment[i]);
+ else if(m_map_Core_MaxE[i]>0)
+ Pos = GetSegmentPosition(clover,m_map_Core_Crystal[i],zero);
+
+ if(Pos.Mag()!=0){
+ static TVector3 Beta = TVector3(0,0,0.10);
+ double E = GetDopplerCorrectedEnergy(m_map_E[i],Pos,Beta);
+ AddBack_DC.push_back(E);
+ AddBack_E.push_back(m_map_E[i]);
+ AddBack_Theta.push_back(Pos.Angle(Beta)*180./3.141592653589793);
+ AddBack_Clover.push_back(clover);
+ if(m_map_Segment_MaxE[i]>0){
+ AddBack_Crystal.push_back(m_map_Segment_Crystal[i]);
+ AddBack_Segment.push_back(m_map_Segment[i]);
+ }
+ else{
+ AddBack_Crystal.push_back(m_map_Core_Crystal[i]);
+ AddBack_Segment.push_back(0);
+ }
+ AddBack_X.push_back(Pos.X());
+ AddBack_Y.push_back(Pos.Y());
+ AddBack_Z.push_back(Pos.Z());
+ }
+ }
+ }
+}
+
+/////////////////////////////////////////////////
+void TGeTAMUPhysics::PreTreat(){
+ static CalibrationManager* cal = CalibrationManager::getInstance();
+ static string name;
+ unsigned int mysize = m_EventData->GetMultiplicityGe();
+ double Eraw,Energy;
+ int clover, crystal, segment;
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ Eraw = m_EventData->GetGeEnergy(i);
+ if( Eraw>20000){
+ clover = m_EventData->GetGeCloverNbr(i);
+ crystal = m_EventData->GetGeCrystalNbr(i);
+ segment = m_EventData->GetGeSegmentNbr(i);
+ name = "GETAMU/D"+ NPL::itoa(clover)+"_CRY"+ NPL::itoa(crystal)+"_SEG"+ NPL::itoa(segment)+"_E";
+ Energy = cal->ApplyCalibration(name, Eraw);
+ Gamma_Energy.push_back(Energy);
+ Clover_Number.push_back(clover);
+ Crystal_Number.push_back(crystal);
+ Segment_Number.push_back(segment);
+ Gamma_Time.push_back(m_EventData->GetGeTimeCFD(i));
+
+ // Look for Associate BGO
+ bool BGOcheck = false ;
+ for(unsigned j = 0 ; j < m_EventData->GetMultiplicityBGO() ; j++){
+
+ if( m_EventData->GetBGOCloverNbr(j)== m_EventData->GetGeCloverNbr(i) && m_EventData->GetBGOEnergy(j)>20 )
+ BGOcheck = true ;
+ }
+ BGO.push_back(BGOcheck);
+ }
+ }
+}
+
+/////////////////////////////////////////////////
+TVector3 TGeTAMUPhysics::GetPositionOfInteraction(unsigned int& i){
+ return GetSegmentPosition(Clover_Number[i],Crystal_Number[i],Segment_Number[i]);
+}
+/////////////////////////////////////////////////
+// original energy, position, beta
+double TGeTAMUPhysics::GetDopplerCorrectedEnergy(double& energy , TVector3 position, TVector3& beta){
+ // renorm pos vector
+ position.SetMag(1);
+ m_GammaLV.SetPx(energy*position.X());
+ m_GammaLV.SetPy(energy*position.Y());
+ m_GammaLV.SetPz(energy*position.Z());
+ m_GammaLV.SetE(energy);
+ m_GammaLV.Boost(-beta);
+ return m_GammaLV.Energy();
+}
+
+/////////////////////////////////////////////////
+void TGeTAMUPhysics::AddClover(unsigned int ID ,double R, double Theta, double Phi){
+ TVector3 Pos(0,0,1);
+ Pos.SetTheta(Theta);
+ Pos.SetPhi(Phi);
+ Pos.SetMag(R);
+ m_CloverPosition[ID] = Pos;
+ return;
+}
+/////////////////////////////////////////////////
+TVector3 TGeTAMUPhysics::GetCloverPosition(int& CloverNbr){
+ return m_CloverPosition[CloverNbr];
+}
+/////////////////////////////////////////////////
+TVector3 TGeTAMUPhysics::GetCorePosition(int& CloverNbr,int& CoreNbr){
+ static double offset = 33.4; // mm
+ static double depth = 45;
+ static TVector3 Pos;
+ TVector3 CloverPos = m_CloverPosition[CloverNbr];
+
+ if(CoreNbr==1)
+ Pos.SetXYZ(-offset,offset,depth);
+ else if(CoreNbr==2)
+ Pos.SetXYZ(offset,offset,depth);
+ else if(CoreNbr==3)
+ Pos.SetXYZ(offset,-offset,depth);
+ else if(CoreNbr==4)
+ Pos.SetXYZ(-offset,-offset,depth);
+ else
+ cout << "Warning: GeTAMU crystal number " << CoreNbr << " is out of range (1 to 4)" << endl;
+
+ // Define reference axis as the clover direction
+ Pos.RotateUz(CloverPos.Unit());
+ Pos+=CloverPos;
+ return (Pos);
+}
+/////////////////////////////////////////////////
+TVector3 TGeTAMUPhysics::GetSegmentPosition(int& CloverNbr,int& CoreNbr, int& SegmentNbr){
+ static double offsetXY1 = 10.4; // mm
+ static double offsetXY2 = 16.7; // mm
+ static double offsetZ1 = 15.5; // mm
+ static double offsetZ2 = 60.5; // mm
+ TVector3 CorePos = GetCorePosition(CloverNbr,CoreNbr);
+ TVector3 CloverPos = GetCloverPosition(CloverNbr);
+ static TVector3 Pos;
+
+ if(SegmentNbr == 0 || SegmentNbr == 9)
+ return GetCorePosition(CloverNbr,CoreNbr);
+ else if(SegmentNbr==1)
+ Pos.SetXYZ(-offsetXY1,offsetXY1,offsetZ1);
+ else if(SegmentNbr==2)
+ Pos.SetXYZ(offsetXY1,offsetXY1,offsetZ1);
+ else if(SegmentNbr==3)
+ Pos.SetXYZ(offsetXY1,-offsetXY1,offsetZ1);
+ else if(SegmentNbr==4)
+ Pos.SetXYZ(-offsetXY1,-offsetXY1,offsetZ1);
+ else if(SegmentNbr==5)
+ Pos.SetXYZ(-offsetXY2,offsetXY2,offsetZ2);
+ else if(SegmentNbr==6)
+ Pos.SetXYZ(offsetXY2,offsetXY2,offsetZ2);
+ else if(SegmentNbr==7)
+ Pos.SetXYZ(offsetXY2,-offsetXY2,offsetZ2);
+ else if(SegmentNbr==8)
+ Pos.SetXYZ(-offsetXY2,-offsetXY2,offsetZ2);
+ else
+ cout << "Warning: GeTAMU segment number " << SegmentNbr << " is out of range (0 to 9)" << endl;
+
+
+ // Each crystal is a rotation of the previous one
+ if (CoreNbr == 2 )
+ Pos.RotateZ(90*deg);
+ else if (CoreNbr == 3 )
+ Pos.RotateZ(180*deg);
+ else if (CoreNbr == 4)
+ Pos.RotateZ(270*deg);
+
+ // Define reference axis as the core position
+ Pos.RotateUz(CorePos.Unit());
+ Pos+=CorePos;
+ return (Pos);
+
+}
+
+
+/////////////////////////////////////////////////
+void TGeTAMUPhysics::ReadConfiguration(string Path) {
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+
+ if(!ConfigFile.is_open()) cout << "Config File not Found" << endl ;
+
+ string LineBuffer;
+ string DataBuffer;
+
+ bool check_CloverId= false;
+ bool check_R= false;
+ bool check_Theta= false;
+ bool check_Phi= false;
+ bool ReadingStatus = true;
+
+ int CloverId=0;
+ double R=0;
+ double Theta=0;
+ double Phi=0;
+
+ while (!ConfigFile.eof()) {
+
+ getline(ConfigFile, LineBuffer);
+ // If line is a Start Up GeTAMU bloc, Reading toggle to true
+ if (LineBuffer.compare(0, 13, "GeTAMUClover") == 0) {
+ cout << "///" << endl ;
+ cout << "GeTAMU Clover found: " << endl ;
+ ReadingStatus = true ;
+ }
+ // Else don't toggle to Reading Block Status
+ else ReadingStatus = false ;
+
+ // Reading Block
+ while(ReadingStatus) {
+ // Pickup Next Word
+
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+
+ if (DataBuffer.compare(0, 1, "%") == 0) {ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ // Finding another Clover toggle out (safety)
+ else if (DataBuffer.compare(0, 13, "GeTAMUClover") == 0) {
+ cout << "WARNING: Another Detector is find before standard sequence of Token, Error may occured in Clover definition" << endl ;
+ ReadingStatus = false ;
+ }
+
+ else if (DataBuffer=="CloverId=") {
+ check_CloverId = true;
+ ConfigFile >> DataBuffer ;
+ CloverId=atoi(DataBuffer.c_str());
+ cout << "CloverId: " << CloverId << endl;
+ }
+
+ else if (DataBuffer=="R=") {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str());
+ cout << "R: " << R << "mm" << endl;
+ }
+
+ else if (DataBuffer=="Theta=") {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str());
+ cout << "Theta: " << Theta << "deg" << endl;
+ }
+
+ else if (DataBuffer=="Phi=") {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str());
+ cout << "Phi: " << Phi << "deg" << endl;
+ }
+
+ ///////////////////////////////////////////////////
+ // If no Detector Token and no comment, toggle out
+ else {
+ ReadingStatus = false; cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl ;}
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ( check_Theta && check_Phi && check_R && check_CloverId) {
+ ReadingStatus = false;
+ check_CloverId= false;
+ check_R= false;
+ check_Theta= false;
+ check_Phi= false;
+ AddClover(CloverId,R,Theta*3.141592653589793/180.,Phi*3.141592653589793/180.);
+ }
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TGeTAMUPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus( "GeTAMU" , true );
+ if(inputChain->FindBranch( "fTIG_*" ))
+ inputChain->SetBranchStatus( "fTIG_*" , true );
+ inputChain->SetBranchAddress( "GeTAMU" , &m_EventData );
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TGeTAMUPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch( "GeTAMU" , "TGeTAMUPhysics" , &m_EventPhysics );
+}
+///////////////////////////////////////////////////////////////////////////
+void TGeTAMUPhysics::Clear() {
+ Gamma_Energy.clear();
+ Gamma_Time.clear();
+ Crystal_Number.clear();
+ Clover_Number.clear();
+ Segment_Number.clear();
+ BGO.clear();
+ AddBack_E.clear();
+ AddBack_DC.clear();
+ AddBack_Theta.clear();
+ AddBack_X.clear();
+ AddBack_Y.clear();
+ AddBack_Z.clear();
+ m_map_E.clear();
+ m_map_Core_Crystal.clear();
+ m_map_Core_MaxE.clear();
+ m_map_Segment_Crystal.clear();
+ m_map_Segment.clear();
+ m_map_Segment_MaxE.clear();
+ AddBack_Clover.clear();
+ AddBack_Crystal.clear();
+ AddBack_Segment.clear();
+}
+///////////////////////////////////////////////////////////////////////////
+void TGeTAMUPhysics::ClearEventData() {
+ m_EventData->Clear();
+ m_PreTreatedData->Clear();
+}
+///////////////////////////////////////////////////////////////////////////
+void TGeTAMUPhysics::AddParameterToCalibrationManager(){
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ for(int i = 0 ; i < 16; ++i){
+ for(int cry = 0 ; cry < 4 ; cry++){
+ // core are 0 and 9 , segment 1 to 8
+ for( int j = 0 ; j < 10 ; ++j){
+ Cal->AddParameter("GETAMU", "D"+ NPL::itoa(i+1)+"_CRY"+NPL::itoa(cry+1)+"_SEG"+ NPL::itoa(j)+"_E","GETAMU_D"+ NPL::itoa(i+1)+"_CRY"+NPL::itoa(cry+1)+"_SEG"+NPL::itoa(j)+"_E");
+ }
+ }
+ }
+ return;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TGeTAMUPhysics::Construct(){
+ return (NPL::VDetector*) new TGeTAMUPhysics();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+class proxy_getamu{
+ public:
+ proxy_getamu(){
+ NPL::DetectorFactory::getInstance()->AddToken("GeTAMU","GeTAMU");
+ NPL::DetectorFactory::getInstance()->AddDetector("GeTAMU",TGeTAMUPhysics::Construct);
+ }
+};
+
+proxy_getamu p;
+}
+
diff --git a/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..84c5d82e2f478d776befbffe3fa1ef9b3cbb174d
--- /dev/null
+++ b/NPLib/Detectors/GeTAMU/TGeTAMUPhysics.h
@@ -0,0 +1,138 @@
+#ifndef TGETAMUPHYSICS_H
+#define TGETAMUPHYSICS_H
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * Peter Bender contact address: bender@triumf.ca *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold GeTAMU treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+#include <vector>
+#include <map>
+#include <string>
+using namespace std;
+
+// NPL
+#include "TGeTAMUData.h"
+#include "NPCalibrationManager.h"
+#include "NPVDetector.h"
+
+// ROOT
+#include "TObject.h"
+#include "TVector3.h"
+#include "TLorentzVector.h"
+
+class TGeTAMUPhysics : public TObject, public NPL::VDetector{
+
+ public:
+ TGeTAMUPhysics();
+ ~TGeTAMUPhysics() { };
+
+ void Clear();
+ void Clear(const Option_t*){Clear();}
+
+ public: // inherited from VDetector
+ // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+ void ReadConfiguration(string);
+
+ // Add Parameter to the CalibrationManger
+ void AddParameterToCalibrationManager();
+
+ // Activated associated Branches and link it to the private member DetectorData address
+ // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
+ void InitializeRootInputRaw() ;
+
+ // Activated associated Branches and link it to the private member DetectorPhysics address
+ // In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
+ void InitializeRootInputPhysics() {};
+
+ // Create associated branches and associated private member DetectorPhysics address
+ void InitializeRootOutput() ;
+
+ // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
+ void BuildPhysicalEvent() ;
+
+ // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
+ // This method aimed to be used for analysis performed during experiment, when speed is requiered.
+ // NB: This method can eventually be the same as BuildPhysicalEvent.
+ void BuildSimplePhysicalEvent(){BuildPhysicalEvent();} ;
+
+ // Clear the Event Physics
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() ;
+
+ public:
+ void PreTreat();
+
+ private: // Root Input and Output tree classes
+
+ TGeTAMUData* m_EventData;//!
+ TGeTAMUData* m_PreTreatedData;//!
+ TGeTAMUPhysics* m_EventPhysics;//!
+
+ public: // Data Member
+ vector<double> Gamma_Energy;
+ vector<int> Crystal_Number;
+ vector<int> Clover_Number;
+ vector<int> Segment_Number;
+ vector<bool> BGO;
+ vector<double> Gamma_Time;
+
+ // add back by clover
+ vector<double> AddBack_E;
+ vector<double> AddBack_DC;
+ vector<double> AddBack_Theta;
+ vector<double> AddBack_X;
+ vector<double> AddBack_Y;
+ vector<double> AddBack_Z;
+ vector<int> AddBack_Clover;
+ vector<int> AddBack_Crystal;
+ vector<int> AddBack_Segment;
+
+ private: // use for anlysis
+ // Keep track of the core
+ map<int,double> m_map_E; //!
+ map<int,int> m_map_Core_Crystal; //!
+ map<int,double> m_map_Core_MaxE; //!
+
+ // Keep track of the segment
+ map<int,int> m_map_Segment_Crystal; //!
+ map<int,int> m_map_Segment; //!
+ map<int,double> m_map_Segment_MaxE; //!
+
+
+ TLorentzVector m_GammaLV; //!
+ public:
+ TVector3 GetPositionOfInteraction(unsigned int& i);
+ double GetDopplerCorrectedEnergy(double& energy , TVector3 position, TVector3& beta);
+ // Add a detector and computes its coordinate
+ void AddClover(unsigned int ID, double R, double Theta, double Phi);
+ TVector3 GetCloverPosition(int& CloverNbr);
+ TVector3 GetCorePosition(int& CloverNbr, int& CoreNbr);
+ TVector3 GetSegmentPosition(int& CloverNbr, int& CoreNbr, int& SegmentNbr);
+ inline TVector3 GetCrystalPosition(int& CloverNbr, int& CoreNbr){return GetCorePosition(CloverNbr,CoreNbr);};
+
+ private:
+ map<unsigned int,TVector3> m_CloverPosition;//!
+
+ public: // Static constructor to be passed to the Detector Factory
+ static NPL::VDetector* Construct();
+ ClassDef(TGeTAMUPhysics,1) // GeTAMUPhysics structure
+};
+
+
+
+#endif
diff --git a/NPSimulation/Detectors/GeTAMU/CMakeLists.txt b/NPSimulation/Detectors/GeTAMU/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8dca6410672fbf7db864b5a1c91a9fbfbbae5d97
--- /dev/null
+++ b/NPSimulation/Detectors/GeTAMU/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_library(NPSGeTAMU SHARED GeTAMU.cc )
+target_link_libraries(NPSGeTAMU NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPGeTAMU)
diff --git a/NPSimulation/Detectors/GeTAMU/GeTAMU.cc b/NPSimulation/Detectors/GeTAMU/GeTAMU.cc
new file mode 100644
index 0000000000000000000000000000000000000000..739885d582d6f6ae3cb772f5b8d8bb3dd637ea4e
--- /dev/null
+++ b/NPSimulation/Detectors/GeTAMU/GeTAMU.cc
@@ -0,0 +1,992 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe the GeTAMU Germanium array *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <sstream>
+#include <cmath>
+#include <limits>
+//G4 Geometry object
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4Trd.hh"
+#include "G4Trap.hh"
+#include "G4Cons.hh"
+
+//G4 sensitive
+#include "G4SDManager.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+//G4 various object
+#include "G4Material.hh"
+#include "G4Polycone.hh"
+#include "G4Polyhedra.hh"
+#include "G4LogicalVolume.hh"
+#include "G4ThreeVector.hh"
+#include "G4Transform3D.hh"
+#include "G4RotationMatrix.hh"
+#include "G4PVPlacement.hh"
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+#include "G4RunManager.hh"
+#include "G4ios.hh"
+#include "G4SubtractionSolid.hh"
+#include "G4IntersectionSolid.hh"
+#include "G4UnionSolid.hh"
+#include "G4ThreeVector.hh"
+
+// NPS
+#include "GeTAMU.hh"
+//#include "GeTAMUScorers.hh"
+#include "MaterialManager.hh"
+#include "NPSDetectorFactory.hh"
+
+// NPL
+#include "NPOptionManager.h"
+#include "RootOutput.h"
+using namespace GETAMU;
+
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+
+using namespace std;
+using namespace CLHEP;
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+namespace {
+
+ // Ge crystal
+ // Cylindrical part
+ const G4double CrystalOuterRadius = 30.0*mm; // outer radius for crystal
+ const G4double CrystalInnerRadius = 5.0*mm; // inner radius for hole in crystal
+ const G4double CrystalLength = 90.0*mm; // crystal length
+ const G4double CrystalHoleDepth = 15.0*mm; // depth at which starts the hole
+ //const G4double CrystaHoleRadius = 0*cm;
+ //const G4double CrystalInterDistance = 0.6*mm; // Distance between two crystal
+
+ // Squared part
+ const G4double CrystalWidth = 56.5*mm; // Width of one crystal
+
+ // Exogam Stuff
+ const G4double CrystalEdgeOffset1 = 26.0*mm; // distance of the edge from the center of the crystal
+ const G4double CrystalEdgeOffset2 = 28.5*mm; // distance of the edge from the center of the crystal
+
+ const G4double CapsuleWidth = 1.5*mm; // capsule width
+ const G4double CapsuleLength = 110.*mm; // capsule length
+ const G4double CapsuleEdgeDepth = 3.3*cm; // same as crystal !
+ const G4double CrystalToCapsule = 7*mm; // to be adjusted ..
+
+ //const G4double BGOLength = 120.0*mm;
+ //const G4double BGOWidth = 25.0*mm;
+
+ //const G4double CsILength = 20.0*mm;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// GeTAMU Specific Method
+GeTAMU::GeTAMU(){
+ InitializeMaterial();
+ m_Event = new TGeTAMUData();
+
+ BlueVisAtt = new G4VisAttributes(G4Colour(0, 0, 1)) ;
+ GreenVisAtt = new G4VisAttributes(G4Colour(0, 1, 0)) ;
+ RedVisAtt = new G4VisAttributes(G4Colour(1, 0, 0)) ;
+ WhiteVisAtt = new G4VisAttributes(G4Colour(1, 1, 1)) ;
+ TrGreyVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5)) ;
+
+ m_LogicClover = 0;
+
+}
+
+GeTAMU::~GeTAMU(){
+ delete m_MaterialVacuum;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of NPS::VDetector class
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// Called in DetecorConstruction::ReadDetextorConfiguration Method
+void GeTAMU::ReadConfiguration(string Path){
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+ string LineBuffer ;
+ string DataBuffer ;
+ istringstream LineStream ;
+ // Standard Case:
+ bool check_CloverId = false;
+
+ vector<int> CloverId;
+ int CloverId_Free = 0;
+ double R = 0;
+ double Theta = 0;
+ double Phi = 0;
+ double BetaX=0;
+ double BetaY=0;
+ double BetaZ=0;
+
+ // Free postion case:
+ bool check_R = false ;
+ bool check_Theta = false ;
+ bool check_Phi = false ;
+ bool check_Beta = false ;
+
+ // Frame Case
+ bool check_RightFrame = false ;
+ bool check_LeftFrame = false ;
+
+ bool ReadingStatusStandard = false ;
+ bool ReadingStatusFree = false ;
+ bool ReadingStatusFrame = false ;
+ bool ReadingStatus = false ;
+
+ while (!ConfigFile.eof()){
+ int VerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+
+ getline(ConfigFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 7, "GeTAMU") == 0)
+ ReadingStatus = true;
+
+ while (ReadingStatus && !ConfigFile.eof()) {
+ getline(ConfigFile, LineBuffer);
+
+ // Comment Line
+ while (LineBuffer.compare(0, 1, "%") == 0) {
+ // Take the next line
+ getline(ConfigFile, LineBuffer);
+ }
+
+ // Standard case
+ if (LineBuffer.compare(0, 15, "GeTAMUStandard") == 0){
+ if(VerboseLevel==1)
+ G4cout << "/// Clovers in Standard Configuration : ///" << G4endl ;
+ ReadingStatusStandard = true ;
+ }
+
+ // Free placing case
+ else if (LineBuffer.compare(0, 13, "GeTAMUClover") == 0){
+ if(VerboseLevel==1)
+ G4cout << "/// Free placed clover : ///" << G4endl ;
+ ReadingStatusFree = true ;
+ }
+
+ // Frame case
+ else if (LineBuffer.compare(0, 12, "GeTAMUFrame") == 0){
+ if(VerboseLevel==1)
+ G4cout << "/// Support Frame : ///" << G4endl ;
+ ReadingStatusFrame = true ;
+ }
+
+ // Reading Block
+ while(ReadingStatusStandard){
+ // Pickup Next Line
+ getline(ConfigFile, LineBuffer);
+ // Comment Line
+ while (LineBuffer.compare(0, 1, "%") == 0) {
+ // Take the next line
+ getline(ConfigFile, LineBuffer);
+ }
+
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+
+ if ( DataBuffer == "CloverId=" ) {
+ check_CloverId = true;
+
+ if(VerboseLevel==1) G4cout << "CloverId: " ;
+ while(LineStream >> DataBuffer){
+ CloverId.push_back(atoi(DataBuffer.c_str()));
+ if(VerboseLevel==1) G4cout << atoi(DataBuffer.c_str()) << " ";
+ }
+ if(VerboseLevel==1) G4cout << G4endl << G4endl;
+ }
+
+ ///////////////////////////////////////////////////
+ // If no Detector Token and no comment, toggle out
+ else{
+ ReadingStatusStandard = false;
+ G4cout << "Error: Wrong Token Sequence: Getting out " << DataBuffer << G4endl ;
+ exit(1);
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+
+ if (check_CloverId){
+ ReadingStatusStandard = false;
+ AddCloverStandard(CloverId);
+ CloverId.clear();
+ check_CloverId = false ;
+ }
+ }
+
+ // Reading Block
+ while(ReadingStatusFree){
+ // Pickup Next Line
+ getline(ConfigFile, LineBuffer);
+ // Comment Line
+ while (LineBuffer.compare(0, 1, "%") == 0) {
+ // Take the next line
+ getline(ConfigFile, LineBuffer);
+ }
+
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+
+ if ( DataBuffer == "CloverId=" ) {
+ check_CloverId = true;
+ LineStream >> DataBuffer;
+ CloverId_Free = atoi(DataBuffer.c_str());
+ if(VerboseLevel==1)
+ G4cout << "CloverId: " << atoi(DataBuffer.c_str()) << " " << G4endl ;
+ }
+
+ else if ( DataBuffer == "R=" ) {
+ check_R = true;
+ LineStream >> DataBuffer;
+ R = atof(DataBuffer.c_str())*mm;
+ if(VerboseLevel==1)
+ G4cout << "R: " << R/mm << " " << G4endl ;
+ }
+
+ else if ( DataBuffer == "Theta=" ) {
+ check_Theta = true;
+ LineStream >> DataBuffer;
+ Theta = atof(DataBuffer.c_str())*deg;
+ if(VerboseLevel==1)
+ G4cout << "Theta: " << Theta/deg << " " << G4endl ;
+ }
+
+ else if ( DataBuffer == "Phi=" ) {
+ check_Phi = true;
+ LineStream >> DataBuffer;
+ Phi = atof(DataBuffer.c_str())*deg;
+ if(VerboseLevel==1)
+ G4cout << "Phi: " << Phi/deg << " " << G4endl ;
+ }
+
+ else if ( DataBuffer == "Beta=" ) {
+ check_Beta = true;
+ LineStream >> DataBuffer;
+ BetaX = atof(DataBuffer.c_str())*deg;
+ if(VerboseLevel==1)
+ G4cout << "BetaX: " << BetaX/deg << " " << G4endl ;
+ LineStream >> DataBuffer;
+ BetaY = atof(DataBuffer.c_str())*deg;
+ if(VerboseLevel==1)
+ G4cout << "BetaY: " << BetaY/deg << " " << G4endl ;
+ LineStream >> DataBuffer;
+ BetaZ = atof(DataBuffer.c_str())*deg;
+ if(VerboseLevel==1)
+ G4cout << "BetaZ: " << BetaZ/deg << " " << G4endl ;
+ }
+
+
+ ///////////////////////////////////////////////////
+ // If no Detector Token and no comment, toggle out
+ else{
+ ReadingStatusStandard = false;
+ G4cout << "Error: Wrong Token Sequence: Getting out " << DataBuffer << G4endl ;
+ exit(1);
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+
+ if (check_CloverId && check_R && check_Theta && check_Phi && check_Beta){
+ ReadingStatusFree = false;
+ AddCloverFreePosition(CloverId_Free,R,Theta,Phi,BetaX,BetaY,BetaZ);
+ check_CloverId = false ;
+ check_R = false ;
+ check_Theta = false ;
+ check_Phi = false ;
+ check_Beta = false ;
+ }
+ }
+
+ // Reading Block
+ while(ReadingStatusFrame){
+ // Pickup Next Line
+ getline(ConfigFile, LineBuffer);
+ // Comment Line
+ while (LineBuffer.compare(0, 1, "%") == 0) {
+ // Take the next line
+ getline(ConfigFile, LineBuffer);
+ }
+
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+
+ if ( DataBuffer == "RightFrame=" ) {
+ check_RightFrame = true;
+
+ LineStream >> DataBuffer;
+ m_RightFrame=atoi(DataBuffer.c_str());
+
+ if (VerboseLevel==1) {
+ if(m_RightFrame)
+ G4cout << "Right frame: yes" << G4endl;
+ else
+ G4cout << "Right frame: no" << G4endl;
+ }
+ }
+
+ else if ( DataBuffer == "LeftFrame=" ) {
+ check_LeftFrame = true;
+
+ LineStream >> DataBuffer;
+ m_LeftFrame=atoi(DataBuffer.c_str());
+
+ if (VerboseLevel==1) {
+ if(m_LeftFrame)
+ G4cout << "Left frame: yes" << G4endl;
+ else
+ G4cout << "Left frame: no" << G4endl;
+ }
+ }
+
+ ///////////////////////////////////////////////////
+ // If no Detector Token and no comment, toggle out
+ else{
+ ReadingStatusStandard = false;
+ G4cout << "Error: Wrong Token Sequence: Getting out " << DataBuffer << G4endl ;
+ exit(1);
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+
+ if (check_RightFrame && check_LeftFrame){
+ ReadingStatusFrame = false;
+ AddCloverStandard(CloverId);
+ CloverId.clear();
+ check_RightFrame = false;
+ check_LeftFrame = false;
+ }
+ }
+
+ }
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return a G4VSolid modeling the Crystal
+G4LogicalVolume* GeTAMU::ConstructCrystal(){
+ G4Tubs* Crystal_Cylinder = new G4Tubs("Crystal_Cylinder", 0, CrystalOuterRadius, CrystalLength*0.5, 0, 2*M_PI);
+
+ // Central Hole for cold finger
+ G4RotationMatrix* BoxRotation = new G4RotationMatrix(0,0,0);
+ G4Tubs* Crystal_Hole = new G4Tubs("Crystal_Hole", 0, CrystalInnerRadius, (CrystalLength-CrystalHoleDepth)*0.5, 0, 2*M_PI);
+ G4SubtractionSolid* Crystal_Stage1 = new G4SubtractionSolid("Crystal_Stage1",Crystal_Cylinder,Crystal_Hole,BoxRotation,G4ThreeVector(0,0,CrystalHoleDepth));
+
+ // Flat surface on the side
+ G4Box* Crystal_Box1 = new G4Box("Crystal_Box1", CrystalWidth*0.6, CrystalWidth*0.6,CrystalLength*0.6);
+ G4SubtractionSolid* Crystal_Stage2 = new G4SubtractionSolid("Crystal_Stage2",Crystal_Stage1,Crystal_Box1,BoxRotation,G4ThreeVector(24.5+CrystalWidth*0.6,0,0));
+ G4SubtractionSolid* Crystal_Stage3 = new G4SubtractionSolid("Crystal_Stage3",Crystal_Stage2,Crystal_Box1,BoxRotation,G4ThreeVector(-29-CrystalWidth*0.6,0,0));
+ G4SubtractionSolid* Crystal_Stage4 = new G4SubtractionSolid("Crystal_Stage4",Crystal_Stage3,Crystal_Box1,BoxRotation,G4ThreeVector(0,29+CrystalWidth*0.6,0));
+ G4SubtractionSolid* Crystal_Stage5 = new G4SubtractionSolid("Crystal_Stage5",Crystal_Stage4,Crystal_Box1,BoxRotation,G4ThreeVector(0,-24.5-CrystalWidth*0.6,0));
+
+ // Bezel
+ G4RotationMatrix* BoxRotation1 = new G4RotationMatrix(0,0,0);
+ BoxRotation1->rotate(22.5*deg,G4ThreeVector(1,0,0));
+ G4SubtractionSolid* Crystal_Stage6= new G4SubtractionSolid("Crystal_Stage6",Crystal_Stage5,Crystal_Box1,BoxRotation1,G4ThreeVector(0,20.54*mm+CrystalWidth*0.6,-45*mm));
+
+ G4RotationMatrix* BoxRotation2 = new G4RotationMatrix(0,0,0);
+ BoxRotation2->rotate(22.5*deg,G4ThreeVector(0,1,0));
+ G4SubtractionSolid* Crystal_Stage7= new G4SubtractionSolid("Crystal_Stage7",Crystal_Stage6,Crystal_Box1,BoxRotation2,G4ThreeVector(-20.54*mm-CrystalWidth*0.6,0,-45*mm));
+
+ G4LogicalVolume* logicCrystal =
+ new G4LogicalVolume(Crystal_Stage7,m_MaterialGe,"LogicCrystal", 0, 0, 0);
+
+ return logicCrystal;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return a G4VSolid modeling the Capsule
+G4LogicalVolume* GeTAMU::ConstructCapsule(){
+
+ G4int nbslice = 7;
+ const G4double widthface = 45.5*mm;
+ G4double zSlice[7] = { 0.0*mm,
+ CapsuleWidth-0.1*mm,
+ CapsuleWidth,
+ CapsuleEdgeDepth,
+ CapsuleLength-CapsuleWidth,
+ CapsuleLength-CapsuleWidth-0.1*mm,
+ CapsuleLength };
+
+ G4double InnNullRad[7] = {0,0,0,0,0,0,0};
+
+ G4double InnRad[7] = { 0.0*mm,
+ 0.0*mm,
+ widthface-CapsuleWidth,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule - CapsuleWidth,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule - CapsuleWidth,
+ 0.0*mm,
+ 0.0*mm};
+
+ G4double OutRad[7] = { widthface-1.5*mm,
+ widthface,
+ widthface,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule,
+ CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule};
+
+ // The whole volume of the Capsule, made of N2
+ G4Polyhedra* caps = new G4Polyhedra(G4String("Capsule"), 45.*deg, 360.*deg, 4, nbslice, zSlice, InnNullRad, OutRad);
+ G4LogicalVolume* LogicCapsule=
+ new G4LogicalVolume(caps,m_MaterialN2,"LogicCapsule", 0, 0, 0);
+ LogicCapsule->SetVisAttributes(G4VisAttributes::Invisible);
+
+ // The wall of the Capsule made of Al
+ G4Polyhedra* capsWall = new G4Polyhedra(G4String("CapsuleWall"), 45.*deg, 360.*deg, 4, nbslice, zSlice, InnRad, OutRad);
+ G4LogicalVolume* logicCapsuleWall =
+ new G4LogicalVolume(capsWall,m_MaterialAl,"LogicCapsuleWall", 0, 0, 0);
+
+ new G4PVPlacement(G4Transform3D(*(new G4RotationMatrix()), G4ThreeVector(0,0,0)),
+ logicCapsuleWall,"CapsuleWall",LogicCapsule,false,1);
+ logicCapsuleWall->SetVisAttributes(TrGreyVisAtt);
+
+ return LogicCapsule;
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* GeTAMU::ConstructDewar(){
+ G4Tubs* DewarSolid = new G4Tubs("DewarSolid",0,90*mm*0.5,90*mm*0.5,0,2*M_PI);
+ G4Tubs* DewarCFSolid = new G4Tubs("DewarCFSolid",0,45*mm*0.5,145*mm*0.5,0,2*M_PI);
+
+ G4UnionSolid* DewarFull =
+ new G4UnionSolid("Dewarfull", DewarSolid, DewarCFSolid, new G4RotationMatrix(),G4ThreeVector(0,0,-90*mm-(145-90)*0.5*mm));
+
+ G4LogicalVolume* LogicDewar = new G4LogicalVolume(DewarFull,m_MaterialAl,"LogicDewar",0,0,0);
+
+ G4Tubs* N2Solid = new G4Tubs("N2Solid",0,90*mm*0.5-1*mm,90*mm*0.5-1*mm,0,2*M_PI);
+ G4Tubs* N2CFSolid = new G4Tubs("N2CFSolid",0,45*mm*0.5-1*mm,145*mm*0.5-1*mm,0,2*M_PI);
+
+ G4LogicalVolume* LogicN2 = new G4LogicalVolume(N2Solid,m_MaterialN2,"LogicN2",0,0,0);
+ G4LogicalVolume* LogicN2CF = new G4LogicalVolume(N2CFSolid,m_MaterialN2,"LogicN2CF",0,0,0);
+
+ LogicN2->SetVisAttributes(GreenVisAtt);
+ LogicN2CF->SetVisAttributes(GreenVisAtt);
+ new G4PVPlacement(G4Transform3D(*(new G4RotationMatrix()), G4ThreeVector(0,0,0)),
+ LogicN2,"N2 Deware",LogicDewar,false,1);
+
+ //new G4PVPlacement(G4Transform3D(*(new G4RotationMatrix()), G4ThreeVector(0,0,-90*mm)),
+ // LogicN2CF,"N2 Deware",LogicDewar,false,1);
+
+ LogicDewar->SetVisAttributes(TrGreyVisAtt);
+
+ return LogicDewar;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return a G4VSolid modeling the BGO
+G4LogicalVolume* GeTAMU::ConstructBGO(){
+
+ return NULL;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return a clover in the configuration given by option (not use a the moment)
+void GeTAMU::ConstructClover(){
+
+ if(m_LogicClover==0){
+ // Construct the clover itself
+ m_LogicClover = ConstructCapsule();
+
+ // Place the cristal in the clover
+ double CrystalOffset = (24.5*mm+0.5*mm);
+
+ G4LogicalVolume* logicCrystal = ConstructCrystal();
+
+ G4RotationMatrix* CrystalRotation = new G4RotationMatrix(0,0,0);
+ G4ThreeVector CrystalPositionB = G4ThreeVector(-CrystalOffset,+CrystalOffset,0.5*CrystalLength+7*mm);
+ new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPositionB),
+ logicCrystal,"LogicCrystalB",m_LogicClover,false,1);
+ logicCrystal->SetVisAttributes(BlueVisAtt);
+
+ CrystalRotation->rotate(-90*deg, G4ThreeVector(0,0,1));
+ G4ThreeVector CrystalPositionG = G4ThreeVector(+CrystalOffset,+CrystalOffset,0.5*CrystalLength+7*mm);
+ new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPositionG),
+ logicCrystal,"LogicCrystalG",m_LogicClover,false,2);
+
+ CrystalRotation->rotate(-90*deg, G4ThreeVector(0,0,1));
+ G4ThreeVector CrystalPositionR = G4ThreeVector(+CrystalOffset,-CrystalOffset,0.5*CrystalLength+7*mm);
+ new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPositionR),
+ logicCrystal,"LogicCrystalR",m_LogicClover,false,3);
+
+ CrystalRotation->rotate(-90*deg, G4ThreeVector(0,0,1));
+ G4ThreeVector CrystalPositionW = G4ThreeVector(-CrystalOffset,-CrystalOffset,0.5*CrystalLength+7*mm);
+ new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPositionW),
+ logicCrystal,"LogicCrystalW",m_LogicClover,false,4);
+
+ }
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Construct detector and inialise sensitive part.
+// Called After DetecorConstruction::AddDetector Method
+void GeTAMU::ConstructDetector(G4LogicalVolume* world){
+ ConstructClover();
+
+ G4RotationMatrix* DetectorRotation = new G4RotationMatrix(0,0,0);
+ for (unsigned int i = 0 ; i < m_CloverId.size(); i++) {
+
+ // Constructing the Detector referential and the transition matrix
+ G4ThreeVector U,V,W;
+ G4double wX = sin(m_Theta[i]) * cos(m_Phi[i]) ;
+ G4double wY = sin(m_Theta[i]) * sin(m_Phi[i]) ;
+ G4double wZ = cos(m_Theta[i]);
+ W = G4ThreeVector(wX, wY, wZ) ;
+
+ // vector parallel to one axis of the entrance plane
+ G4double vX = cos(m_Theta[i]) * cos(m_Phi[i]);
+ G4double vY = cos(m_Theta[i]) * sin(m_Phi[i]);
+ G4double vZ = -sin(m_Theta[i]);
+ V = G4ThreeVector(vX, vY, vZ);
+
+ W = W.unit();
+ U = V.cross(W);
+ U = U.unit();
+ V = W.cross(U);
+ V = V.unit();
+ // Passage Matrix from Lab Referential to Clover Referential
+ delete DetectorRotation;
+ DetectorRotation = new G4RotationMatrix(U, V, W);
+
+ DetectorRotation->rotate(m_BetaX[i], U);
+ DetectorRotation->rotate(m_BetaY[i], V);
+ DetectorRotation->rotate(m_BetaZ[i], W);
+ G4ThreeVector DetectorPosition = m_R[i]*W;
+
+ new G4PVPlacement(G4Transform3D(*DetectorRotation, DetectorPosition),
+ m_LogicClover,"Clover",world,false,m_CloverId[i]);
+
+ G4LogicalVolume* LogicDewar = ConstructDewar();
+
+ new G4PVPlacement(G4Transform3D(*DetectorRotation, DetectorPosition+W*((90*mm+(145)*mm)+CapsuleLength*0.5+90*0.5*mm)),
+ LogicDewar,"Dewar",world,false,m_CloverId[i]);
+
+
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add clover at the standard position of the array
+// Take as argument the standard clover Id.
+void GeTAMU::AddCloverStandard(vector<int> CloverId){
+
+ for (unsigned int i = 0 ; i < CloverId.size(); i++) {
+ if(CloverId[i] == 1 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(45*deg);
+ m_Phi.push_back(22.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(0);
+ }
+
+ else if(CloverId[i] == 2 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(45*deg);
+ m_Phi.push_back(112.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(0);
+ }
+
+ else if(CloverId[i] == 3 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(45*deg);
+ m_Phi.push_back(202.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(0);
+ }
+
+ else if(CloverId[i] == 4 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(45*deg);
+ m_Phi.push_back(292.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(0);
+ }
+
+ else if(CloverId[i] == 5 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(22.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 6 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(67.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 7 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(112.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 8 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(157.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 9 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(202.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 10 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(247.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 11 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(292.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 12 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(90*deg);
+ m_Phi.push_back(337.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 13 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(135*deg);
+ m_Phi.push_back(22.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 14 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(135*deg);
+ m_Phi.push_back(112.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 15 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(135*deg);
+ m_Phi.push_back(202.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+
+ else if(CloverId[i] == 16 ){
+ m_CloverId.push_back(CloverId[i]);
+ m_R.push_back(145*mm);
+ m_Theta.push_back(135*deg);
+ m_Phi.push_back(292.5*deg);
+ m_BetaX.push_back(0);
+ m_BetaY.push_back(0);
+ m_BetaZ.push_back(180*deg);
+ }
+ }
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add clover at a free position in space with coordinate
+// in spherical coordinate
+// Beta are the three angles of rotation in the Clover frame
+void GeTAMU::AddCloverFreePosition(int CloverId,double R,double Theta,double Phi,double BetaX,double BetaY,double BetaZ){
+
+ m_CloverId.push_back(CloverId);
+ m_R.push_back(R);
+ m_Theta.push_back(Theta);
+ m_Phi.push_back(Phi);
+ m_BetaX.push_back(BetaX);
+ m_BetaY.push_back(BetaY);
+ m_BetaZ.push_back(BetaZ);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void GeTAMU::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ pTree->Branch("GeTAMU", "TGeTAMUData", &m_Event) ;
+ pTree->SetBranchAddress("GeTAMU", &m_Event) ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void GeTAMU::ReadSensitive(const G4Event* event){
+ event->GetHCofThisEvent(); // event should be used to remove compilation warning
+ /*m_Event->Clear();
+
+ ///////////
+ // BOX
+ G4THitsMap<G4double*>* BOXHitMap;
+ std::map<G4int, G4double**>::iterator BOX_itr;
+
+ G4int BOXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("GeTAMU_BOXScorer/GeTAMUBOX");
+ BOXHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(BOXCollectionID));
+
+ // Loop on the BOX map
+ for (BOX_itr = BOXHitMap->GetMap()->begin() ; BOX_itr != BOXHitMap->GetMap()->end() ; BOX_itr++){
+
+ G4double* Info = *(BOX_itr->second);
+
+ double Energy = Info[0];
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+ int StripFront = (int) Info[3];
+ int StripBack = (int) Info[4];
+
+ m_Event->SetFront_DetectorNbr(DetNbr);
+ m_Event->SetFront_StripNbr(StripFront);
+ m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ m_Event->SetBack_DetectorNbr(DetNbr);
+ m_Event->SetBack_StripNbr(StripBack);
+ m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+
+ // Interraction Coordinates
+ ms_InterCoord->SetDetectedPositionX(Info[5]) ;
+ ms_InterCoord->SetDetectedPositionY(Info[6]) ;
+ ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
+ ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
+ ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
+
+ }
+
+ // clear map for next event
+ BOXHitMap->clear();
+
+ ///////////
+ // PAD
+ G4THitsMap<G4double*>* PADHitMap;
+ std::map<G4int, G4double**>::iterator PAD_itr;
+
+ G4int PADCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("GeTAMU_PADScorer/GeTAMUPAD");
+ PADHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(PADCollectionID));
+
+ // Loop on the BOX map
+ for (PAD_itr = PADHitMap->GetMap()->begin() ; PAD_itr != PADHitMap->GetMap()->end() ; PAD_itr++){
+
+ G4double* Info = *(PAD_itr->second);
+
+ double Energy = Info[0];
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+
+ m_Event->SetPAD_DetectorNbr(DetNbr);
+ m_Event->SetPAD_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetPAD_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetPAD_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ }
+
+ // clear map for next event
+ PADHitMap->clear();
+
+ ///////////
+ // QQQ
+ G4THitsMap<G4double*>* QQQHitMap;
+ std::map<G4int, G4double**>::iterator QQQ_itr;
+
+ G4int QQQCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("GeTAMU_QQQScorer/GeTAMUQQQ");
+ QQQHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(QQQCollectionID));
+
+ // Loop on the BOX map
+ for (QQQ_itr = QQQHitMap->GetMap()->begin() ; QQQ_itr != QQQHitMap->GetMap()->end() ; QQQ_itr++){
+
+ G4double* Info = *(QQQ_itr->second);
+
+ double Energy = Info[0];
+ double Time = Info[1];
+ int DetNbr = (int) Info[2];
+ int StripFront = (int) Info[3];
+ int StripBack = (int) Info[4];
+
+ m_Event->SetFront_DetectorNbr(DetNbr);
+ m_Event->SetFront_StripNbr(StripFront);
+ m_Event->SetFront_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetFront_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetFront_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ m_Event->SetBack_DetectorNbr(DetNbr);
+ m_Event->SetBack_StripNbr(StripBack);
+ m_Event->SetBack_Energy(RandGauss::shoot(Energy, ResoEnergy));
+ m_Event->SetBack_TimeCFD(RandGauss::shoot(Time, ResoTime));
+ m_Event->SetBack_TimeLED(RandGauss::shoot(Time, ResoTime));
+
+ // Interraction Coordinates
+ ms_InterCoord->SetDetectedPositionX(Info[5]) ;
+ ms_InterCoord->SetDetectedPositionY(Info[6]) ;
+ ms_InterCoord->SetDetectedPositionZ(Info[7]) ;
+ ms_InterCoord->SetDetectedAngleTheta(Info[8]/deg) ;
+ ms_InterCoord->SetDetectedAnglePhi(Info[9]/deg) ;
+
+ }
+
+ // clear map for next event
+ QQQHitMap->clear();
+ */
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void GeTAMU::InitializeScorers(){
+ /*
+ // Silicon Associate Scorer
+ m_BOXScorer = new G4MultiFunctionalDetector("GeTAMU_BOXScorer");
+ m_PADScorer = new G4MultiFunctionalDetector("GeTAMU_PADScorer");
+ m_QQQScorer = new G4MultiFunctionalDetector("GeTAMU_QQQScorer");
+
+ G4VPrimitiveScorer* BOXScorer =
+ new GeTAMU::PS_Silicon_Rectangle("GeTAMUBOX",
+ BOX_Wafer_Length,
+ BOX_Wafer_Width,
+ BOX_Wafer_Back_NumberOfStrip ,
+ BOX_Wafer_Front_NumberOfStrip,
+ EnergyThreshold);
+
+ G4VPrimitiveScorer* PADScorer =
+ new GeTAMU::PS_Silicon_Rectangle("GeTAMUPAD",
+ PAD_Wafer_Length,
+ PAD_Wafer_Width,
+ 1 ,
+ 1,
+ EnergyThreshold);
+
+ G4VPrimitiveScorer* QQQScorer =
+ new GeTAMU::PS_Silicon_Annular("GeTAMUQQQ",
+ QQQ_Wafer_Inner_Radius,
+ QQQ_Wafer_Outer_Radius,
+ QQQ_Wafer_Stopping_Phi-QQQ_Wafer_Starting_Phi,
+ QQQ_Wafer_NumberOf_RadialStrip,
+ QQQ_Wafer_NumberOf_AnnularStrip,
+ EnergyThreshold);
+
+
+
+ //and register it to the multifunctionnal detector
+ m_BOXScorer->RegisterPrimitive(BOXScorer);
+ m_PADScorer->RegisterPrimitive(PADScorer);
+ m_QQQScorer->RegisterPrimitive(QQQScorer);
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_BOXScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_PADScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_QQQScorer) ;*/
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+/////////////////Material Definition ///////////////////////////
+////////////////////////////////////////////////////////////////
+void GeTAMU::InitializeMaterial(){
+ m_MaterialVacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
+ m_MaterialGe= MaterialManager::getInstance()->GetMaterialFromLibrary("Ge");
+ m_MaterialAl= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ m_MaterialN2= MaterialManager::getInstance()->GetMaterialFromLibrary("N2_liquid");
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPS::VDetector* GeTAMU::Construct(){
+ return (NPS::VDetector*) new GeTAMU();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern"C" {
+ class proxy_nps_getamu{
+ public:
+ proxy_nps_getamu(){
+ NPS::DetectorFactory::getInstance()->AddToken("GeTAMU","GeTAMU");
+ NPS::DetectorFactory::getInstance()->AddDetector("GeTAMU",GeTAMU::Construct);
+ }
+ };
+
+ proxy_nps_getamu p_nps_getamu;
+}
diff --git a/NPSimulation/Detectors/GeTAMU/GeTAMU.hh b/NPSimulation/Detectors/GeTAMU/GeTAMU.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c184f50f2839de5ee14674e3b5314d7ddc540e7
--- /dev/null
+++ b/NPSimulation/Detectors/GeTAMU/GeTAMU.hh
@@ -0,0 +1,174 @@
+#ifndef GeTAMU_h
+#define GeTAMU_h 1
+/*****************************************************************************
+ * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : November 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe the GeTAMU Silicon detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+#include <vector>
+
+// G4 header defining G4 types
+#include "globals.hh"
+
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4VisAttributes.hh"
+#include "G4VSolid.hh"
+
+// NPSimulation header
+#include "NPSVDetector.hh"
+
+// NPLib
+#include "TGeTAMUData.h"
+using namespace std;
+using namespace CLHEP;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+namespace GETAMU{
+ // Energy and time Resolution
+ const G4double ResoTime = 0 ;
+ const G4double ResoEnergy = 0.035*MeV ;// = zzkeV of Resolution // Unit is MeV/2.35
+ const G4double EnergyThreshold = 0.4*MeV;
+}
+
+using namespace GETAMU ;
+class GeTAMU : public NPS::VDetector
+{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+public:
+ GeTAMU() ;
+ ~GeTAMU() ;
+
+ ////////////////////////////////////////////////////
+ //////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // Add clover at the standard position of the array
+ // Take as argument the standard clover Id.
+ void AddCloverStandard(vector<int> CloverId);
+
+ // Add clover at a free position in space with coordinate
+ // in spherical coordinate
+ // Beta are the three angles of rotation in the Clover frame
+ void AddCloverFreePosition(int CloverId,
+ double R,
+ double Theta,
+ double Phi,
+ double BetaX,
+ double BetaY,
+ double BetaZ);
+
+ // Return a clover in the configuration given by option (not use a the moment)
+ void ConstructClover();
+
+ // Return a modeling of the Crystal
+ G4LogicalVolume* ConstructCrystal();
+
+ // Return a modeling of the Capsule
+ G4LogicalVolume* ConstructCapsule();
+
+ // Return a modeling of the Dewar
+ G4LogicalVolume* ConstructDewar();
+
+ // Return a G4VSolid modeling the BGO
+ G4LogicalVolume* ConstructBGO();
+
+ ////////////////////////////////////////////////////
+ ///////// Inherite from NPS::VDetector class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path) ;
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world) ;
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event) ;
+
+
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+private:
+ TGeTAMUData* m_Event ;
+
+ ////////////////////////////////////////////////////
+ ///////////////// Scorer Related ///////////////////
+ ////////////////////////////////////////////////////
+
+private:
+ // Geometry related
+ G4LogicalVolume* m_LogicClover;
+
+
+ // Initialize all Scorer
+ void InitializeScorers() ;
+
+ // Scorer Associate to the Silicon
+ //G4MultiFunctionalDetector* m_GermaniumScorer ;
+ //G4MultiFunctionalDetector* m_BGOScorer ;
+private:
+ // Initialize material used in detector definition
+ void InitializeMaterial();
+
+ // List of material
+ G4Material* m_MaterialVacuum ;
+ G4Material* m_MaterialGe;
+ G4Material* m_MaterialAl;
+ G4Material* m_MaterialN2;
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private:
+ // Clover Position
+ vector<int> m_CloverId;
+ vector<double> m_R;
+ vector<double> m_Theta;
+ vector<double> m_Phi;
+ vector<double> m_BetaX;
+ vector<double> m_BetaY;
+ vector<double> m_BetaZ;
+
+ // Frame: true if the frame is to be done
+ bool m_LeftFrame;
+ bool m_RightFrame;
+
+private:/// Visualisation Attribute:
+ G4VisAttributes* BlueVisAtt;
+ G4VisAttributes* GreenVisAtt;
+ G4VisAttributes* RedVisAtt;
+ G4VisAttributes* WhiteVisAtt;
+ G4VisAttributes* TrGreyVisAtt;
+public:
+ static NPS::VDetector* Construct();
+};
+#endif
diff --git a/Projects/S1554/macro/Mg.cxx b/Projects/S1554/macro/Mg.cxx
deleted file mode 100644
index e6a8a87770a649cb01d65f06286a038a06a2d085..0000000000000000000000000000000000000000
--- a/Projects/S1554/macro/Mg.cxx
+++ /dev/null
@@ -1,121 +0,0 @@
-void Mg(){
-
- TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mgdp_Phy.root");
- TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
-
- tree->SetAlias("GoodTrifoil","(Trifoil.Time>150 && Trifoil.Time<160 && RunNumber>35110) && ELab > 0 ") ;
- tree->SetAlias("GoodProton","ThetaLab>90 && Sharc.Strip_E>0.8");
- TCanvas* c = new TCanvas("Mg","Mg");
- c->Divide(3,3);
-
- // Setting Trifoil entry list for fast drawing
- TFile* fileEL = new TFile("macro/GoodTrifoil.root");
- TEntryList* EL = (TEntryList*) fileEL->FindObjectAny("TrifoilEL");
- if(!EL){
- cout << "Trifoil Entry List does not exist, generating : ";
- tree->Draw(">>TrifoilEL","GoodTrifoil","entrylist");
- EL = (TEntryList*) gDirectory->FindObjectAny("TrifoilEL");
- cout << " done" << endl;
- EL->SaveAs("macro/GoodTrifoil.root");
- }
-
- cout <<"Setting Trifoil entry list :" ;
- tree->SetEntryList(EL);
- cout << " Done" << endl;
-
-
- c->cd(1);
- tree->Draw("Ex:AddBack_DC/1000.>>h2(2000,0,8,200,-1,8)","ThetaLab>90 && GoodProton","colz");
- TH1* h2 = (TH1*) gDirectory->FindObjectAny("h2");
- h2->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
- h2->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
- gPad->SetLogz();
- TLine* line = new TLine(0,0,8,8);
- line -> Draw("c");
-
- c->cd(2);
-
- new TCanvas();
- tree->Draw("Ex>>h(200,-8,8)","GoodTrifoil && GoodProton ","colz");
- TH1* h = (TH1*) gDirectory->FindObjectAny("h");
- h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
- h->GetYaxis()->SetTitle("counts/40 keV");
- return ;
- int maximum = h->GetMaximum();
- TLine* state1= new TLine(0,0,0,maximum);
- state1->Draw("l");
-
- TLine* state2= new TLine(0.054,0,0.054,maximum);
- state2->Draw("l");
-
- TLine* state3= new TLine(1.095,0,1.095,maximum);
- state3->Draw("l");
-
- TLine* state4= new TLine(1.431,0,1.431,maximum);
- state4->Draw("l");
-
- TLine* state5= new TLine(2.266,0,2.266,maximum);
- state5->Draw("l");
-
- TLine* state6= new TLine(2.5,0,2.5,maximum);
- state6->Draw("l");
-
- TLine* state7= new TLine(3.2,0,3.2,maximum);
- state7->Draw("l");
-
- TLine* state8= new TLine(4.43,0,4.43,maximum);
- state8->Draw("l");
-
- c->cd(3);
- tree->Draw("ELab:ThetaLab>>hK(180,90,180,900,0,14)"," GoodProton","colz");
- NPL::Reaction r("28Mg(d,p)29Mg@222");
- r.GetKinematicLine3()->Draw("c");
-
- NPL::Reaction r2("28Mg(d,d)28Mg@222");
- r2.GetKinematicLine3()->Draw("c");
-
- NPL::Reaction r3("28Mg(p,p)28Mg@222");
- r3.GetKinematicLine3()->Draw("c");
-
-
- TH2* hk = (TH2*) gDirectory->FindObjectAny("hK");
- hk->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
- hk->GetYaxis()->SetTitle("E_{Lab} (MeV)");
-
- c->cd(4);
- tree->Draw("AddBack_DC/1000>>hDC(300,0,3)","ELab>0&& ThetaLab>90 && GoodProton && Ex > 0.5 && Ex<2","");
- gPad->SetLogz();
- h = (TH1*) gDirectory->FindObjectAny("hDC");
- h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
- h->GetYaxis()->SetTitle("counts / 10 keV");
-
- c->cd(5);
- tree->Draw("Ex>>h11(50,-1,8)","GoodProton && AddBack_DC/1000. > 1.0 && AddBack_DC/1000. < 1.06 ","");
- h = (TH1*) gDirectory->FindObjectAny("h11");
- h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
- tree->Draw("Ex>>h12(50,-1,8)","GoodProton && AddBack_DC/1000. > 0.3 && AddBack_DC/1000. < 0.36 ","same");
- h = (TH1*) gDirectory->FindObjectAny("h12");
- h->SetFillColor(kBlack);
- h->SetLineColor(kBlack);
- h->SetFillStyle(3004);
-
- c->cd(6);
- tree->Draw("AddBack_DC/1000>>hDC2(300,0,3)","ELab>0&& ThetaLab>90 && GoodProton && Ex > 2 && Ex<3","");
- gPad->SetLogz();
- h = (TH1*) gDirectory->FindObjectAny("hDC2");
- h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
- h->GetYaxis()->SetTitle("counts / 10 keV");
-
-
- c->cd(7);
- tree->Draw("Ex>>h13(50,-1+0.1,8+0.1)","GoodProton && AddBack_DC/1000. > 2.2 && AddBack_DC/1000. < 2.28 ","");
- h = (TH1*) gDirectory->FindObjectAny("h11");
- h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
-
- tree->Draw("Ex>>h14(50,-1+0.1,8+0.1)","GoodProton && AddBack_DC/1000. > 2.35 && AddBack_DC/1000. < 2.65 ","same");
- h = (TH1*) gDirectory->FindObjectAny("h14");
- h->SetFillColor(kBlack);
- h->SetLineColor(kBlack);
- h->SetFillStyle(3004);
-}
-
diff --git a/Projects/S1554/macro/MgCS.cxx b/Projects/S1554/macro/MgCS.cxx
deleted file mode 100644
index adf76990bbe0205f3afb42316a2c22ed92741f0f..0000000000000000000000000000000000000000
--- a/Projects/S1554/macro/MgCS.cxx
+++ /dev/null
@@ -1,344 +0,0 @@
-void Scale(TGraph* g , TGraphErrors* ex);
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
-
-////////////////////////////////////////////////////////////////////////////////
-void MgCS(){
-
-
- file = new TFile("SharcEfficiency_29Mg.root");
- TH2* Eff = (TH2*) file->FindObjectAny("SolidAngleCM_2D");
-
- //TCanvas* c0 = new TCanvas("Mg0","Mg0");
-
- Goodfile = new TFile("GoodExThetaCM.root");
- TH2D* goodEx = (TH2D*) Goodfile->FindObjectAny("hexcmG");
-
- Badfile = new TFile("BadExThetaCM.root");
- TH2D* badEx = (TH2D*) Badfile->FindObjectAny("hexcmB");
-
- TH2D* h = goodEx;
- badEx->Scale(7./(145-1+300-160));
- // h->Sumw2();
- h->SetMarkerSize(0.5);
- h->SetMarkerColor(kAzure+7);
- // badEx->Sumw2();
-// h->Draw("colz");
- h->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
- h->GetXaxis()->SetTitle("#theta_{CM} (deg)");
- h->GetXaxis()->SetRangeUser(-1,50);
-
- vector<double> E,W;
- E.push_back(0.0);
- E.push_back(0.054);
-// E.push_back(0.590);
-
- E.push_back(1.094);
- E.push_back(1.431);
-// E.push_back(1.638);
-
- E.push_back(2.266);
-// E.push_back(2.500);
- E.push_back(2.8);
- E.push_back(3.223);
- E.push_back(3.674);
- E.push_back(3.985);
- E.push_back(4.280);
- E.push_back(4.43);
-
-
- string formula="[0]+[1]*x+";
-
- for (unsigned int i = 0 ; i < E.size() ; i++){
- //formula+=Form("gaus(%i)",i*3);
- formula+=Form("[%i]*exp( -(x-[5]-[%i])*(x-[5]-[%i])/(([2]+[3]*x+[4]*x*x)*([2]+[3]*x+[4]*x*x)))",i*2+6,i*2+6+1,i*2+6+1);
- if(i!=E.size()-1)
- formula+="+";
- W.push_back(1.1*(0.144206+0.00644644*E[i]+0.0006305*E[i]*E[i]));
- }
-
- TGraph* g_shift = new TGraph();
- vector<TGraph*> g_width;
- for (unsigned int i = 0 ; i < E.size() ; i++){
- g_width.push_back(new TGraph);
- }
- TF1* bck = new TF1("bck","[0]+[1]*x",-6,-1);
- bck->SetLineColor(kBlack);
- TF1* func = new TF1("func",formula.c_str(),-2,5);
- func->FixParameter(0,0);
- func->FixParameter(1,0);
-
- func->FixParameter(2,0.226217);
- func->FixParameter(3,0.00294652);
- func->FixParameter(4,0);
- func->FixParameter(5,-0.055);
-
- // func->SetParLimits(2,0.1,0.3);
- // func->SetParLimits(3,0,1);
- // func->SetParLimits(4,0,1e-31);
- // func->SetParLimits(5,-0.06,0.06);
-
- func->SetNpx(500);
-
- double Energy;
- for (unsigned int i = 0 ; i < E.size() ; i++){
- func->SetParameter(i*2+6,10);
- func->SetParLimits(i*2+6,0,100);
- func->FixParameter(i*2+6+1,E[i]);
- }
-
- //////////
- vector<TGraphErrors*> CS;
- CS.resize(E.size(),NULL);
-
- for(unsigned int i = 0 ; i < CS.size() ; i++){
- CS[i]=new TGraphErrors();
- }
-
- const double* par;
- unsigned int nbin = 10;
- double minCM = 2;
- double maxCM = 42;
- double step = (maxCM-minCM)/nbin;
- Eff->RebinX(step);
-
-
- TH1D* p; TH1D* bbb;
- double ThetaCM;
- int pad = 1;
-
- TCanvas* c4 = new TCanvas("Fit","Fit",1800,1000);
- c4->Divide(4,2,1e-4,1e-4);
- c4->cd(1);
- for(unsigned int i = 0 ; i < nbin ; i++){
- ThetaCM = minCM+i*step;
- p = h->ProjectionY(Form("p%i",i),h->GetXaxis()->FindBin(ThetaCM),h->GetXaxis()->FindBin(ThetaCM+step));
- p->Rebin(2);
- // Background
- bbb = badEx->ProjectionY(Form("b%i",i),badEx->GetXaxis()->FindBin(ThetaCM),badEx->GetXaxis()->FindBin(ThetaCM+step));
- bbb->Rebin(2);
- // Substract the Background
- p->Add(bbb,-1);
- // Presentation
- p->SetMarkerSize(0.5);
- p->SetMarkerColor(kAzure+7);
- bbb->SetMarkerSize(0.1);
- bbb->SetLineColor(kOrange+7);
- bbb->SetFillColor(kOrange+7);
- bbb->SetFillStyle(3004);
-
- p->GetYaxis()->SetRangeUser(0,30);
- p->GetXaxis()->SetRangeUser(-1,6);
-
- // Limit the heigh of the gaussian to the bin at the state position
- for(unsigned int n = 0 ; n < CS.size() ; n++){
- func->SetParameter(n*2+6,p->GetBinContent(p->FindBin(func->GetParameter(n*2+6+1))));
- double integral = p->Integral(p->FindBin(E[n]-2*W[n]),p->FindBin(E[n]+2*W[n]));
- if(integral >0)
- func->SetParLimits(n*2+6,0, integral+1e-20);
- else
- func->SetParLimits(n*2+6,0, 10);
- }
- c4->cd(pad++);
- p->Draw();
- // p->Fit(bck,"RBFQ");
- bbb->Draw("same");
-
- // func->FixParameter(0,bck->GetParameter(0));
- // func->FixParameter(1,bck->GetParameter(1));
- p->Fit(func,"RBFQ+");
- par = func->GetParameters();
- TLatex* tex = new TLatex(2,25,Form("%.1f deg",ThetaCM+0.5*step));
- tex->Draw();
- // bck->Draw("same");
-
- g_shift->SetPoint(g_shift->GetN(),ThetaCM,func->GetParameter(5));
- if(p->Integral(p->FindBin(-1),p->FindBin(6))>-1000 ){
- for(unsigned int n = 0 ; n < CS.size() ; n++){
- double SolidAngle = Eff->ProjectionX("px",Eff->GetYaxis()->FindBin(E[n]-W[n]),Eff->GetYaxis()->FindBin(E[n]+W[n]))->Interpolate(ThetaCM);
- double Width = func->GetParameter(2)+func->GetParameter(3)*E[n]+func->GetParameter(4)*E[n]*E[n];
- g_width[n]->SetPoint(g_width[n]->GetN(),ThetaCM,Width);
- CS[n]->SetPoint(CS[n]->GetN(),ThetaCM+0.5*step,par[n*2+6]/SolidAngle);
- CS[n]->SetPointError(CS[n]->GetN()-1,0,func->GetParError(n*2+6)/SolidAngle);
- }
- }
- else{
- pad= pad-1;
- }
- }
-
- // Draw Width of each state vs ThetaCM
- c4->cd(pad++);
- const UInt_t Number = 2;
- Double_t Red[Number] = { 0, 0 };
- Double_t Green[Number] = { 0, 0.8 };
- Double_t Blue[Number] = { 0, 1.00 };
-
- Double_t Length[Number] = { 0,1.00 };
- Int_t nb=g_width.size() ;
- Int_t color = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
-
- for(unsigned int n = 0 ; n < g_width.size() ; n++){
- g_width[n]->SetMarkerColor(color++) ;
- g_width[n]->SetMarkerSize(0.4);
- if(n==0){
- g_width[n]->Draw("ap");
- g_width[n]->GetYaxis()->SetRangeUser(0.3,1);
- }
- else
- g_width[n]->Draw("p");
- }
-
- // Draw the shift of Ex spectrum for each ThetaCM
- c4->cd(pad++);
- g_shift->SetLineColor(kAzure+7);
- g_shift->Draw("acp");
-
-
-
- TCanvas* c = new TCanvas("Mg","Mg",2000,2000*2./5.);
- c->Divide(5,2,0.00001,0.00001);
-
-
- for(unsigned int i = 0 ; i < CS.size() ; i++){
- c->cd(i+1);
-
- CS[i]->Draw("ap");
- gPad->SetLogy();
-
- TGraph* gs = TWOFNR(E[i], 0,0.5,1,0, 0.5);
- TGraph* gp = TWOFNR(E[i], 0,1.5,1,1, 1.5);
- TGraph* gd = TWOFNR(E[i], 0,1.5,0,2, 1.5);
- TGraph* gf = TWOFNR(E[i], 0,3.5,0,3, 3.5);
- TGraph* gg = TWOFNR(E[i], 0,4.5,0,4, 4.5);
-
- Scale(gs,CS[i]);
- Scale(gp,CS[i]);
- Scale(gd,CS[i]);
- Scale(gf,CS[i]);
- Scale(gg,CS[i]);
-
- gs->SetLineColor(kGreen);
- gp->SetLineColor(kBlue);
- gd->SetLineColor(kRed);
- gf->SetLineColor(kBlack);
- gg->SetLineColor(kAzure+7);
-
- if(i!=0){
- gp->Draw("c");
- gd->Draw("c");
- gf->Draw("c");
- gs->Draw("c");
- gg->Draw("c");
- TLegend* leg = new TLegend(0.7,0.6,1,0.9);
- leg->SetBorderSize(0);
- leg->SetFillStyle(0);
- leg->AddEntry(gs,"L=0","l");
- leg->AddEntry(gp,"L=1","l");
- leg->AddEntry(gd,"L=2","l");
- leg->AddEntry(gf,"L=3","l");
- leg->AddEntry(gg,"L=4","l");
- leg->Draw();
- }
-
- else{
- gs->Draw("c");
- gd->Draw("c");
- TGraph* gm = new TGraph(*gs);
- double* XX = gm->GetX();
- double* YY = gm->GetY();
- double* ZZ = gd->GetY();
- for(int n = 0 ; n < gs->GetN() ; n++)
- gm->SetPoint(n,XX[n],YY[n]+ZZ[n]);
- Scale(gm,CS[i]);
- gm->SetLineStyle(2);
- gm->SetLineColor(kGreen-3);
- gm->Draw("c");
- TLegend* leg = new TLegend(0.7,0.6,1,0.9);
- leg->SetBorderSize(0);
- leg->SetFillStyle(0);
- leg->AddEntry(gs,"L=0","l");
- leg->AddEntry(gd,"L=2","l");
- leg->AddEntry(gm,"L=2+0","l");
- leg->Draw();
- }
-
- TLatex* tex = new TLatex(10,1.7e3,Form("%.2f MeV",E[i]));
- tex->Draw();
- CS[i]->GetYaxis()->SetRangeUser(1,1e4);
- CS[i]->GetXaxis()->SetRangeUser(0,40);
- }
-}
-////////////////////////////////////////////////////////////////////////////////
-void Scale(TGraph* g , TGraphErrors* ex){
- double scale;
- double mean = 0 ;
- double* eX = ex->GetX();
- double* eY = ex->GetY();
- double totalW = 0;
- double W = 0;
- for(int i = 0 ; i < ex->GetN() ; i++){
- if(eY[i]>1 && eY[i] <1e4){
- // Incremental Error weighted average
- W = 1./ex->GetErrorY(i);
- scale = eY[i]/g->Eval(eX[i]);
- totalW +=W;
- mean = mean + W*(scale - mean)/(totalW);
- }
- }
-
- double* x = g->GetX();
- double* y = g->GetY();
- for(unsigned int i = 0 ; i < g->GetN() ; i++)
- g->SetPoint(i,x[i],y[i]*mean);
-}
-////////////////////////////////////////////////////////////////////////////////
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
- double BeamEnergy = 8;
-
- NPL::Reaction r("28Mg(d,p)29Mg@224");
- r.SetExcitationHeavy(E);
- double QValue = r.GetQValue();
-
- std::ofstream Front_Input("in.front");
- Front_Input << "jjj" << std::endl;
- Front_Input << "pipo" << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << BeamEnergy << std::endl;
- Front_Input << 28 << " " << 12 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << "0 0 0" << std::endl;
- Front_Input << l << " " << j << std::endl;
- Front_Input << n << std::endl;
- Front_Input << 2 << std::endl;
-
- // unbound case:
- if( QValue < 0 )
- Front_Input << 0.01 << std::endl;
- else
- Front_Input << QValue << std::endl;
-
- Front_Input << 1 << std::endl;
- Front_Input << J0 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 5 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << J << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1.25 << " " << 0.65 << std::endl;
- Front_Input << 6.0 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << 0 << std::endl;
-
- Front_Input.close() ;
-
- system("(exec FRONT < in.front &> /dev/null)");
- system("(exec echo tran.jjj | TWOFNR &> /dev/null)");
-
- TGraph* CS = new TGraph("22.jjj");
- return CS;
-}
diff --git a/Projects/S1554/macro/MgCS2.cxx b/Projects/S1554/macro/MgCS2.cxx
deleted file mode 100644
index c19bfecd911bc03bfd3fb88a2357ef92067ef87a..0000000000000000000000000000000000000000
--- a/Projects/S1554/macro/MgCS2.cxx
+++ /dev/null
@@ -1,273 +0,0 @@
-void Scale(TGraph* g , TGraphErrors* ex);
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
-TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad, TH2* Eff);
-double ToMininize(const double* parameter);
-TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp, double& s1, double& s2);
-
-TGraph* g1;
-TGraph* g2;
-TH1* current;
-
-////////////////////////////////////////////////////////////////////////////////
-void MgCS2(){
- file = new TFile("SharcEfficiency_29Mg.root");
- TH2* Eff = (TH2*) file->FindObjectAny("SolidAngleCM_2D");
-
- Goodfile = new TFile("GoodExThetaCM.root");
- TH2D* goodEx = (TH2D*) Goodfile->FindObjectAny("hexcmG");
-
- Badfile = new TFile("BadExThetaCM.root");
- TH2D* badEx = (TH2D*) Badfile->FindObjectAny("hexcmB");
-
- TH2D* h = goodEx;
- badEx->Scale(7./(145-1+300-160));
- h->Sumw2();
- h->SetMarkerSize(0.5);
- h->SetMarkerColor(kAzure+7);
- badEx->Sumw2();
- h->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
- h->GetXaxis()->SetTitle("#theta_{CM} (deg)");
- h->Add(badEx,-1);
- Eff->Sumw2();
-// //h->Scale(0.010048*10/7.37);
- h->Scale(0.015791*10/7.37);
-
- h->Divide(Eff);
-
- TH1D* p;
- int nn=0;
- // Full Ex
- p = h->ProjectionY(Form("p%i",nn++),h->GetXaxis()->FindBin(0.),h->GetXaxis()->FindBin(50.));
- p->Draw();
-
- new TCanvas();
- h->Draw("colz");
- // CS
- TCanvas* c = new TCanvas("CS","CS",800,800);
- c->Divide(2,2);
- TGraph* CSth;
- TLegend* leg;
- double s1,s2;
-
- // Region 1
- current = Region(-0.1,0.4,h,c,nn++,Eff);
- g1 = TWOFNR(0.000, 0, 0.5, 1, 0, 0.5);
- g2 = TWOFNR(0.054, 0, 1.5, 0, 2, 1.5);
- g1->SetLineStyle(2);g1->Draw("c");
- g2->SetLineStyle(3);g2->Draw("c");
- FindNormalisation(g1,g2,current,s1,s2)->Draw("c");
-
-
- CSth = FindNormalisation(g1,g2,current,s1,s2);
- CSth->Draw("c");
- leg = new TLegend(0.6,0.5,0.9,0.9);
- leg->SetBorderSize(0);
- leg->SetFillStyle(0);
- leg->AddEntry(g1,Form("L=0 S=%.2f",s1),"l");
- leg->AddEntry(g2,Form("L=2 S=%.2f",s2),"l");
- leg->AddEntry(CSth,"L=0 + L=2 ","l");
- leg->Draw();
-
- // Region 2
- current = Region(0.6,1.8,h,c,nn++,Eff);
- g1 = TWOFNR(1.09, 0, 1.5, 1, 1, 1.5);
- g2 = TWOFNR(1.43, 0, 3.5, 0, 3, 3.5);
- g1->SetLineStyle(2);g1->Draw("c");
- g2->SetLineStyle(3);g2->Draw("c");
- FindNormalisation(g1,g2,current,s1,s2)->Draw("c");
-
- CSth = FindNormalisation(g1,g2,current,s1,s2);
- CSth->Draw("c");
- leg = new TLegend(0.6,0.5,0.9,0.9);
- leg->SetBorderSize(0);
- leg->SetFillStyle(0);
- leg->AddEntry(g1,Form("L=1 S=%.2f",s1),"l");
- leg->AddEntry(g2,Form("L=3 S=%.2f",s2),"l");
- leg->AddEntry(CSth,"L=1 + L=3 ","l");
- leg->Draw();
-
- // Region 3
- current = Region(1.8,3.4,h,c,nn++,Eff);
- g1 = TWOFNR(2.27, 0, 1.5, 0, 2, 1.5);
- g2 = TWOFNR(2.50, 0, 0.5, 1, 1, 0.5);
- g1->SetLineStyle(2);g1->Draw("c");
- g2->SetLineStyle(3);g2->Draw("c");
- FindNormalisation(g1,g2,current,s1,s2)->Draw("c");
-
- CSth = FindNormalisation(g1,g2,current,s1,s2);
- CSth->Draw("c");
- leg = new TLegend(0.6,0.5,0.9,0.9);
- leg->SetBorderSize(0);
- leg->SetFillStyle(0);
- leg->AddEntry(g1,Form("L=2 S=%.2f",s1),"l");
- leg->AddEntry(g2,Form("L=1 S=%.2f",s2),"l");
- leg->AddEntry(CSth,"L=2 + L=1 ","l");
- leg->Draw();
-
-
-}
-////////////////////////////////////////////////////////////////////////////////
-TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad , TH2* Eff){
- c->cd(pad);
- TH1D* p = h->ProjectionX(Form("p%i",pad),h->GetYaxis()->FindBin(Emin),h->GetYaxis()->FindBin(Emax));
- p->Draw();
- p->GetYaxis()->SetTitle("d#sigma/d#Omega (mb/sr)");
- gPad->SetLogy();
- p->GetYaxis()->SetRangeUser(1e-1,1e3);
- return p;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double Chi2(TGraph* g , TH1* h){
- double Chi2 = 0 ;
- double chi;
-
- for(int i = 1 ; i < h->GetNbinsX() ; i++){
- if(h->GetBinContent(i)>0) {
- chi = (h->GetBinContent(i) - g->Eval(h->GetBinCenter(i)) ) / (h->GetBinError(i));
- Chi2 +=chi*chi;
- }
- }
-
- return Chi2;
-}
-////////////////////////////////////////////////////////////////////////////////
-double ToMininize(const double* parameter){
-static int f = 0 ;
- TGraph* g = new TGraph();
- double* X = g1->GetX();
- double* Y = g1->GetY();
- for(int i = 0 ; i < g1->GetN() ; i++){
- g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i] + parameter[1]*g2->Eval(X[i]) );
- }
- double chi2 = Chi2(g,current);
- return chi2;
-}
-////////////////////////////////////////////////////////////////////////////////
-void Scale(TGraph* g , TGraphErrors* ex){
- double scale;
- double mean = 0 ;
- double* eX = ex->GetX();
- double* eY = ex->GetY();
- double totalW = 0;
- double W = 0;
- for(int i = 0 ; i < ex->GetN() ; i++){
- if(eY[i]>1 && eY[i] <1e4){
- // Incremental Error weighted average
- W = 1./ex->GetErrorY(i);
- scale = eY[i]/g->Eval(eX[i]);
- totalW +=W;
- mean = mean + W*(scale - mean)/(totalW);
- }
- }
-
- double* x = g->GetX();
- double* y = g->GetY();
- for(unsigned int i = 0 ; i < g->GetN() ; i++)
- g->SetPoint(i,x[i],y[i]*mean);
-}
-////////////////////////////////////////////////////////////////////////////////
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
- double BeamEnergy = 8;
-
- NPL::Reaction r("28Mg(d,p)29Mg@224");
- r.SetExcitationHeavy(E);
- double QValue = r.GetQValue();
-
- std::ofstream Front_Input("in.front");
- Front_Input << "jjj" << std::endl;
- Front_Input << "pipo" << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << BeamEnergy << std::endl;
- Front_Input << 28 << " " << 12 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << "0 0 0" << std::endl;
- Front_Input << l << " " << j << std::endl;
- Front_Input << n << std::endl;
- Front_Input << 2 << std::endl;
-
- // unbound case:
-// if( QValue < 0 )
-// Front_Input << 0.01 << std::endl;
-// else
- Front_Input << QValue << std::endl;
-
- Front_Input << 1 << std::endl;
- Front_Input << J0 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 5 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << J << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1.25 << " " << 0.65 << std::endl;
- Front_Input << 6.0 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << 0 << std::endl;
-
- Front_Input.close() ;
-
- system("(exec FRONT < in.front &> /dev/null)");
- system("(exec echo tran.jjj | TWOFNR &> /dev/null)");
- // system("exec FRONT < in.front");
- // system("(exec echo tran.jjj | TWOFNR)");
-
-
- TGraph* CS = new TGraph("22.jjj");
- return CS;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp, double& s1, double& s2){
- // Set global variable
- g1 = cs1;
- g2 = cs2;
- current = hexp;
-
- const char* minName ="Minuit";const char* algoName="Fumili2";
- ROOT::Math::Minimizer* min =
- ROOT::Math::Factory::CreateMinimizer(minName, algoName);
- min->SetValidError(true);
-
- // Number of parameter
- mysize = 2;
-
- // create funciton wrapper for minimizer
- // a IMultiGenFunction type
- ROOT::Math::Functor f(&ToMininize,mysize);
- min->SetFunction(f);
-
- double* parameter = new double[mysize];
- for(unsigned int i = 0 ; i < mysize ; i++){
- parameter[i] = 0.5;
- char name[4];
- sprintf(name,"S%d",i+1);
- min->SetLimitedVariable(i,name,parameter[i],0.1,0,1000);
- }
-
- // do the minimization
- min->Minimize();
- const double *xs = min->X();
- const double *err = min->Errors();
-
- for(int i = 0 ; i < mysize ; i++){
- cout << Form("S%d=",i+1) << xs[i] <<"("<<err[i] << ") ";
- }
- cout << endl;
- // Return the Fitted CS
- TGraph* g = new TGraph();
- double* X = cs1->GetX();
- double* Y = cs1->GetY();
- for(int i = 0 ; i < cs1->GetN() ; i++){
- g->SetPoint(g->GetN(),X[i],xs[0]*Y[i] + xs[1]*cs2->Eval(X[i]) );
- }
- s1 = xs[0];
- s2 = xs[1];
-
- return g;
-}
-
diff --git a/Projects/S1554/macro/SM.cxx b/Projects/S1554/macro/SM.cxx
index 752a1580da02e7935da2bc51c519df9333edbcba..63908d94d7d87b6d1a0a664815b2b2f0ee8407e9 100644
--- a/Projects/S1554/macro/SM.cxx
+++ b/Projects/S1554/macro/SM.cxx
@@ -4,7 +4,7 @@
void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collection, vector<int> matched);
vector<TPad*> thePad;
double Emin = -0.2;
-double Emax = 10;
+double Emax = 3;
double ScaleMin = 0;
double ScaleMax = 7;
double LabelSize = 0.04;
@@ -75,15 +75,15 @@ void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collecti
string final_label;
if(state.GetOrder()!=0){
if(state.GetParity()>0)
- final_label = Form("%d^{+}_{%d} %s %.2f" ,(int)state.GetJ() ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ final_label = Form("%d/2^{+}_{%d} %s %.2f" ,(int)(2*state.GetJ()) ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
else
- final_label = Form("%d^{-}_{%d} %s %.2f" ,(int)state.GetJ() ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ final_label = Form("%d/2^{-}_{%d} %s %.2f" ,(int)(2*state.GetJ()) ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
}
else{
if(state.GetParity()>0)
- final_label = Form("%d^{+} %s %.2f" ,(int)state.GetJ() ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ final_label = Form("%d/2^{+} %s %.2f" ,(int)(2*state.GetJ()) ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
else
- final_label = Form("%d^{-} %s %.2f" ,(int)state.GetJ() ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ final_label = Form("%d/2^{-} %s %.2f" ,(int)(2*state.GetJ()) ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
}
@@ -99,7 +99,7 @@ void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collecti
unsigned int orbsize = state.GetNumberOfOrbital();
// Where the last one stopped
- int thickness = 4;
+ int thickness = 5;
for(unsigned int orb = 0 ; orb < orbsize ; orb++){
TLine* Level;
double SF = state.GetOrbitalS(orb);
@@ -111,11 +111,14 @@ void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collecti
}
else if(orbital.find("p")!=string::npos){
color = kBlue;
- }if(orbital.find("d")!=string::npos){
+ }
+ else if(orbital.find("d")!=string::npos){
color = kGreen+2;
- }if(orbital.find("f")!=string::npos){
- color = kOrange+7;
- }if(orbital.find("g")!=string::npos){
+ }
+ else if(orbital.find("f")!=string::npos){
+ color = kBlack;
+ }
+ else if(orbital.find("g")!=string::npos){
color = kMagenta+2;
}
diff --git a/Projects/SharcEfficiency/Analysis.cxx b/Projects/SharcEfficiency/Analysis.cxx
index 18ef52b81e0570790cec7bd57012e4d50772fe0d..eb6c07d23572ae7db4ad91469fd70390d5175f23 100644
--- a/Projects/SharcEfficiency/Analysis.cxx
+++ b/Projects/SharcEfficiency/Analysis.cxx
@@ -73,15 +73,14 @@ void Analysis::Init(){
X_Trifoil = 0;
Y_Trifoil = 0 ;
-
Si_E_Sharc = 0 ;
E_Sharc = 0;
ThetaDetector = 0 ;
BeamDirection = TVector3(0,0,1);
// S1554
- // TargetPosition = TVector3(0.1635909,0.910980,m_DetectorManager->GetTargetZ() );
+ TargetPosition = TVector3(0.1635909,0.910980,m_DetectorManager->GetTargetZ() );
// S1107
- TargetPosition = TVector3(0.0808323,0.177073,m_DetectorManager->GetTargetZ() );
+ //TargetPosition = TVector3(0.0808323,0.177073,m_DetectorManager->GetTargetZ() );
double finalEnergy = BeamCD2.Slow(myReaction->GetBeamEnergy(),TargetThickness*0.5,0);
myReaction->SetBeamEnergy(finalEnergy);
cout << "Set Beam energy to: " << finalEnergy << " MeV" << endl;
@@ -100,10 +99,10 @@ void Analysis::Init(){
ThetaLab_emmitted_2D = new TH2F("ThetaLab_emmitted_2D","ThetaLab_emmitted_2D",72,0,180,400,-8,8);
ThetaLab_detected_2D = new TH2F("ThetaLab_detected_2D","ThetaLab_detected_2D",72,0,180,400,-8,8);
*/
- ThetaCM_emmitted_2D = new TH2F("ThetaCM_emmitted_2D","ThetaCM_emmitted_2D",180,0,180,250,-1,10);
- ThetaCM_detected_2D = new TH2F("ThetaCM_detected_2D","ThetaCM_detected_2D",180,0,180,250,-1,10);
- ThetaLab_emmitted_2D = new TH2F("ThetaLab_emmitted_2D","ThetaLab_emmitted_2D",180,0,180,250,-1,10);
- ThetaLab_detected_2D = new TH2F("ThetaLab_detected_2D","ThetaLab_detected_2D",180,0,180,250,-1,10);
+ ThetaCM_emmitted_2D = new TH2F("ThetaCM_emmitted_2D","ThetaCM_emmitted_2D",180,0,180,1100,-1,10);
+ ThetaCM_detected_2D = new TH2F("ThetaCM_detected_2D","ThetaCM_detected_2D",180,0,180,1100,-1,10);
+ ThetaLab_emmitted_2D = new TH2F("ThetaLab_emmitted_2D","ThetaLab_emmitted_2D",180,0,180,1100,-1,10);
+ ThetaLab_detected_2D = new TH2F("ThetaLab_detected_2D","ThetaLab_detected_2D",180,0,180,1100,-1,10);
Kine_2D = new TH2F("Kine_2D","Kine_2D",1800,0,180,600,0,60);
@@ -120,7 +119,7 @@ void Analysis::TreatEvent(){
ThetaCM_emmitted->Fill(myInit->GetThetaCM(0));
ThetaLab_emmitted->Fill(myInit->GetThetaLab_WorldFrame(0));
- myReaction->SetBeamEnergy(myInit->GetIncidentFinalKineticEnergy());
+ myReaction->SetBeamEnergy(myInit->GetIncidentInitialKineticEnergy());
double EXD = myReaction->ReconstructRelativistic(myInit->GetKineticEnergy(0),myInit->GetThetaLab_IncidentFrame(0)*deg);
ThetaCM_emmitted_2D->Fill(myInit->GetThetaCM(0),EXD);
ThetaLab_emmitted_2D->Fill(myInit->GetThetaLab_WorldFrame(0),EXD);
@@ -232,7 +231,7 @@ void Analysis::TreatEvent(){
*/
check =true;
if(Trifoil->Energy.size()>0){
- if( check && abs(Ex-EXD)<0.5 && Trifoil->Energy.size()>0 && Trifoil->Energy[0]>0) { // for S1107
+ if( abs(Ex-EXD)<0.5 && Trifoil->Energy.size()>0 && Trifoil->Energy[0]>0) {
ThetaCM_detected->Fill(myInit->GetThetaCM(0));
ThetaLab_detected->Fill(myInit->GetThetaLab_WorldFrame(0));
ThetaCM_detected_2D->Fill(myInit->GetThetaCM(0),EXD);