diff --git a/NPLib/Charissa/Makefile b/NPLib/Charissa/Makefile
index b0dd77cd7d857d4b2abdd7a5458bbb85c0256c72..b251010c9da427044978c1d296f6039ff7469aa8 100644
--- a/NPLib/Charissa/Makefile
+++ b/NPLib/Charissa/Makefile
@@ -1,16 +1,14 @@
include ../Makefile.arch
#------------------------------------------------------------------------------
-#SHARELIB = libCharissaData.so libCharissaPhysics.so libCharissaSpectra.so
-SHARELIB = libCharissa.so
+SHARELIB = libCharissa.so
all: $(SHARELIB)
#------------------------------------------------------------------------------
############### Detector ##############
## MUST2 ##
-#libCharissa.so: TCharissaData.o TCharissaDataDict.o TCharissaPhysics.o TCharissaPhysicsDict.o TCharissaSpectra.o
-libCharissa.so: TCharissaData.o TCharissaDataDict.o TCharissaSpectra.o
+libCharissa.so: TCharissaData.o TCharissaDataDict.o TCharissaPhysics.o TCharissaPhysicsDict.o TCharissaSpectra.o
$(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
TCharissaDataDict.cxx: TCharissaData.h
diff --git a/NPLib/Charissa/TCharissaPhysics.cxx b/NPLib/Charissa/TCharissaPhysics.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..64a922ba1e02f9a80783db4b4fe0c5ba2f8e4e90
--- /dev/null
+++ b/NPLib/Charissa/TCharissaPhysics.cxx
@@ -0,0 +1,1409 @@
+/*****************************************************************************
+ * 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 : December 2013 *
+ * *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold Charissa treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+#include "TCharissaPhysics.h"
+using namespace CHARISSA_LOCAL;
+
+// STL
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+
+// NPL
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TAsciiFile.h"
+#include "NPOptionManager.h"
+
+// ROOT
+#include "TChain.h"
+///////////////////////////////////////////////////////////////////////////
+
+ClassImp(TCharissaPhysics)
+///////////////////////////////////////////////////////////////////////////
+TCharissaPhysics::TCharissaPhysics(){
+ EventMultiplicity = 0 ;
+ m_EventData = new TCharissaData ;
+ m_PreTreatedData = new TCharissaData ;
+ m_EventPhysics = this ;
+ m_Spectra = NULL;
+ m_NumberOfTelescope = 0 ;
+ m_MaximumStripMultiplicityAllowed = 10;
+ m_StripEnergyMatchingSigma = 0.020 ;
+ m_StripEnergyMatchingNumberOfSigma = 3;
+ // Raw Threshold
+ m_Si_X_E_RAW_Threshold = 8200 ;
+ m_Si_Y_E_RAW_Threshold = 8200 ;
+ m_CsI_E_RAW_Threshold = 8200 ;
+ // Calibrated Threshold
+ m_Si_X_E_Threshold = 0 ;
+ m_Si_Y_E_Threshold = 0 ;
+ m_CsI_E_Threshold = 0;
+ m_NumberOfStrip = 16;
+
+ m_Take_E_Y=false;
+ m_Take_T_Y=true;
+
+ }
+
+
+
+///////////////////////////////////////////////////////////////////////////
+TCharissaPhysics::~TCharissaPhysics(){
+
+}
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::BuildSimplePhysicalEvent(){
+ BuildPhysicalEvent();
+}
+
+///////////////////////////////////////////////////////////////////////////
+
+void TCharissaPhysics::BuildPhysicalEvent(){
+ PreTreat();
+
+ bool check_CSI = false ;
+
+ m_Layer1_StripXEMult = m_PreTreatedData->GetCharissaLayer1StripXEMult();
+ m_Layer1_StripYEMult = m_PreTreatedData->GetCharissaLayer1StripYEMult();
+ m_Layer1_StripXTMult = m_PreTreatedData->GetCharissaLayer1StripXTMult();
+ m_Layer1_StripYTMult = m_PreTreatedData->GetCharissaLayer1StripYTMult();
+ m_Layer2_StripXEMult = m_PreTreatedData->GetCharissaLayer2StripXEMult();
+ m_Layer2_StripYEMult = m_PreTreatedData->GetCharissaLayer2StripYEMult();
+ m_Layer2_StripXTMult = m_PreTreatedData->GetCharissaLayer2StripXTMult();
+ m_Layer2_StripYTMult = m_PreTreatedData->GetCharissaLayer2StripYTMult();
+ m_CsIEMult = m_PreTreatedData->GetCharissaCsIEMult();
+ m_CsITMult = m_PreTreatedData->GetCharissaCsITMult();
+
+ // Layer 1
+ if( Layer1_CheckEvent() == 1 ){
+ vector< TVector2 > Layer1_couple = Layer1_Match_X_Y() ;
+ EventMultiplicity = Layer1_couple.size();
+
+ for(unsigned int i = 0 ; i < Layer1_couple.size() ; ++i){
+ check_CSI = false ;
+
+ int Layer1_N = m_PreTreatedData->GetCharissaLayer1StripXEDetectorNbr(Layer1_couple[i].X()) ;
+
+ int Layer1_X = m_PreTreatedData->GetCharissaLayer1StripXEStripNbr(Layer1_couple[i].X()) ;
+ int Layer1_Y = m_PreTreatedData->GetCharissaLayer1StripYEStripNbr(Layer1_couple[i].Y()) ;
+
+ double Layer1_Si_X_E = m_PreTreatedData->GetCharissaLayer1StripXEEnergy( Layer1_couple[i].X() ) ;
+ double Layer1_Si_Y_E = m_PreTreatedData->GetCharissaLayer1StripYEEnergy( Layer1_couple[i].Y() ) ;
+
+ // Search for associate Time
+ double Layer1_Si_X_T = -1000 ;
+
+ for(unsigned int t = 0 ; t < m_Layer1_StripXTMult ; ++t )
+ {
+ if( m_PreTreatedData->GetCharissaLayer1StripXTStripNbr( Layer1_couple[i].X() ) == m_PreTreatedData->GetCharissaLayer1StripXTStripNbr(t)
+ ||m_PreTreatedData->GetCharissaLayer1StripXTDetectorNbr( Layer1_couple[i].X() ) == m_PreTreatedData->GetCharissaLayer1StripXTDetectorNbr(t))
+ Layer1_Si_X_T = m_PreTreatedData->GetCharissaLayer1StripXTTime(t);
+ }
+
+ double Layer1_Si_Y_T = -1000 ;
+
+ for(unsigned int t = 0 ; t < m_Layer1_StripYTMult ; ++t )
+ {
+ if( m_PreTreatedData->GetCharissaLayer1StripYTStripNbr( Layer1_couple[i].Y() ) == m_PreTreatedData->GetCharissaLayer1StripYTStripNbr(t)
+ ||m_PreTreatedData->GetCharissaLayer1StripYTDetectorNbr( Layer1_couple[i].Y() ) == m_PreTreatedData->GetCharissaLayer1StripYTDetectorNbr(t))
+ Layer1_Si_Y_T = m_PreTreatedData->GetCharissaLayer1StripYTTime(t);
+ }
+
+ Layer1_Si_X.push_back(Layer1_X) ; Layer1_Si_Y.push_back(Layer1_Y) ; Layer1_TelescopeNumber.push_back(Layer1_N) ;
+
+ if(m_Take_E_Y) Layer1_Si_E.push_back(Layer1_Si_Y_E);
+ else Layer1_Si_E.push_back(Layer1_Si_X_E);
+
+ if(m_Take_T_Y) Layer1_Si_T.push_back(Layer1_Si_Y_T);
+ else Layer1_Si_T.push_back(Layer1_Si_X_T);
+
+ // Store the other value for checking purpose
+ Layer1_Si_EX.push_back(Layer1_Si_X_E);
+ Layer1_Si_TX.push_back(Layer1_Si_X_T);
+ Layer1_Si_EY.push_back(Layer1_Si_Y_E);
+ Layer1_Si_TY.push_back(Layer1_Si_Y_T);
+ }
+ }
+
+
+ // Layer 2
+ if( Layer2_CheckEvent() == 1 ){
+ vector< TVector2 > Layer2_couple = Layer2_Match_X_Y() ;
+ EventMultiplicity = Layer2_couple.size();
+
+ for(unsigned int i = 0 ; i < Layer2_couple.size() ; ++i){
+
+ check_CSI = false ;
+
+ int Layer2_N = m_PreTreatedData->GetCharissaLayer2StripXEDetectorNbr(Layer2_couple[i].X()) ;
+
+ int Layer2_X = m_PreTreatedData->GetCharissaLayer2StripXEStripNbr(Layer2_couple[i].X()) ;
+ int Layer2_Y = m_PreTreatedData->GetCharissaLayer2StripYEStripNbr(Layer2_couple[i].Y()) ;
+
+ double Layer2_Si_X_E = m_PreTreatedData->GetCharissaLayer1StripXEEnergy( Layer2_couple[i].X() ) ;
+ double Layer2_Si_Y_E = m_PreTreatedData->GetCharissaLayer1StripYEEnergy( Layer2_couple[i].Y() ) ;
+
+ // Search for associate Time
+ double Layer2_Si_X_T = -1000 ;
+
+ for(unsigned int t = 0 ; t < m_Layer2_StripXTMult ; ++t )
+ {
+ if( m_PreTreatedData->GetCharissaLayer2StripXTStripNbr( Layer2_couple[i].X() ) == m_PreTreatedData->GetCharissaLayer2StripXTStripNbr(t)
+ ||m_PreTreatedData->GetCharissaLayer2StripXTDetectorNbr( Layer2_couple[i].X() ) == m_PreTreatedData->GetCharissaLayer2StripXTDetectorNbr(t))
+ Layer2_Si_X_T = m_PreTreatedData->GetCharissaLayer2StripXTTime(t);
+ }
+
+ double Layer2_Si_Y_T = -1000 ;
+
+ for(unsigned int t = 0 ; t < m_Layer2_StripYTMult ; ++t )
+ {
+ if( m_PreTreatedData->GetCharissaLayer2StripYTStripNbr( Layer2_couple[i].Y() ) == m_PreTreatedData->GetCharissaLayer2StripYTStripNbr(t)
+ ||m_PreTreatedData->GetCharissaLayer2StripYTDetectorNbr( Layer2_couple[i].Y() ) == m_PreTreatedData->GetCharissaLayer2StripYTDetectorNbr(t))
+ Layer2_Si_Y_T = m_PreTreatedData->GetCharissaLayer2StripYTTime(t);
+ }
+
+ Layer2_Si_X.push_back(Layer2_X) ; Layer2_Si_Y.push_back(Layer2_Y) ; Layer2_TelescopeNumber.push_back(Layer2_N) ;
+
+ if(m_Take_E_Y) Layer2_Si_E.push_back(Layer2_Si_Y_E);
+ else Layer2_Si_E.push_back(Layer2_Si_X_E);
+
+ if(m_Take_T_Y) Layer2_Si_T.push_back(Layer2_Si_Y_T);
+ else Layer2_Si_T.push_back(Layer2_Si_X_T);
+
+ // Store the other value for checking purpose
+ Layer2_Si_EX.push_back(Layer2_Si_X_E);
+ Layer2_Si_TX.push_back(Layer2_Si_X_T);
+ Layer2_Si_EY.push_back(Layer2_Si_Y_E);
+ Layer2_Si_TY.push_back(Layer2_Si_Y_T);
+
+ for(unsigned int j = 0 ; j < m_CsIEMult ; ++j){
+ if(m_PreTreatedData->GetCharissaCsIEDetectorNbr(j)==Layer2_N){
+ if(Match_Si_CsI( Layer2_X, Layer2_Y , m_PreTreatedData->GetCharissaCsIECristalNbr(j) ) ){
+ CsI_N.push_back( m_PreTreatedData->GetCharissaCsIECristalNbr(j) ) ;
+ CsI_E.push_back( m_PreTreatedData->GetCharissaCsIEEnergy(j) ) ;
+
+ // Look for associate Time
+ for(unsigned int k =0 ; k < m_CsITMult ; ++k){
+ // Same Cristal; Same Detector
+ if( m_PreTreatedData->GetCharissaCsIECristalNbr(j)==m_PreTreatedData->GetCharissaCsITCristalNbr(k)
+ && m_PreTreatedData->GetCharissaCsIEDetectorNbr(j)==m_PreTreatedData->GetCharissaCsITDetectorNbr(k) )
+ CsI_T.push_back( m_PreTreatedData->GetCharissaCsITTime(j) ) ; break ;
+ }
+ check_CSI = true ;
+ }
+ }
+ }
+
+ if(!check_CSI){
+ CsI_N.push_back(0) ;
+ CsI_E.push_back(-1000) ;
+ CsI_T.push_back(-1000) ;
+ }
+ }
+ }
+
+
+
+
+
+ return;
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::PreTreat(){
+
+ ClearPreTreatedData();
+ m_Layer1_StripXEMult = m_EventData->GetCharissaLayer1StripXEMult();
+ m_Layer1_StripYEMult = m_EventData->GetCharissaLayer1StripYEMult();
+ m_Layer1_StripXTMult = m_EventData->GetCharissaLayer1StripXTMult();
+ m_Layer1_StripYTMult = m_EventData->GetCharissaLayer1StripYTMult();
+ m_Layer2_StripXEMult = m_PreTreatedData->GetCharissaLayer2StripXEMult();
+ m_Layer2_StripYEMult = m_PreTreatedData->GetCharissaLayer2StripYEMult();
+ m_Layer2_StripXTMult = m_PreTreatedData->GetCharissaLayer2StripXTMult();
+ m_Layer2_StripYTMult = m_PreTreatedData->GetCharissaLayer2StripYTMult();
+
+ m_CsIEMult = m_EventData->GetCharissaCsIEMult();
+ m_CsITMult = m_EventData->GetCharissaCsITMult();
+
+ // X Layer1
+ // E
+ for(unsigned int i = 0 ; i < m_Layer1_StripXEMult ; ++i){
+ if( m_EventData->GetCharissaLayer1StripXEEnergy(i)>m_Si_X_E_RAW_Threshold && IsValidChannel("X", m_EventData->GetCharissaLayer1StripXEDetectorNbr(i), m_EventData->GetCharissaLayer1StripXEStripNbr(i)) ){
+ double Layer1_EX = fSi_X_E(m_EventData , i);
+ if( Layer1_EX > m_Si_X_E_Threshold ){
+ m_PreTreatedData->SetCharissaLayer1StripXEDetectorNbr( m_EventData->GetCharissaLayer1StripXEDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripXEStripNbr( m_EventData->GetCharissaLayer1StripXEStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripXEEnergy( Layer1_EX );
+ }
+
+ }
+ }
+
+ // T
+ for(unsigned int i = 0 ; i < m_Layer1_StripXTMult ; ++i){
+ if(IsValidChannel("X", m_EventData->GetCharissaLayer1StripXTDetectorNbr(i), m_EventData->GetCharissaLayer1StripXTStripNbr(i))){
+ m_PreTreatedData->SetCharissaLayer1StripXTDetectorNbr( m_EventData->GetCharissaLayer1StripXTDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripXTStripNbr( m_EventData->GetCharissaLayer1StripXTStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripXTTime( fSi_X_T(m_EventData , i) );
+ }
+ }
+
+
+ // Y Layer1
+ // E
+ for(unsigned int i = 0 ; i < m_Layer1_StripYEMult ; ++i){
+ if( m_EventData->GetCharissaLayer1StripYEEnergy(i)<m_Si_Y_E_RAW_Threshold && IsValidChannel("Y", m_EventData->GetCharissaLayer1StripYEDetectorNbr(i), m_EventData->GetCharissaLayer1StripYEStripNbr(i))){
+ double Layer1_EY = fSi_Y_E(m_EventData , i);
+ if( Layer1_EY >m_Si_Y_E_Threshold ){
+ m_PreTreatedData->SetCharissaLayer1StripYEDetectorNbr( m_EventData->GetCharissaLayer1StripYEDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripYEStripNbr( m_EventData->GetCharissaLayer1StripYEStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripYEEnergy( Layer1_EY );
+ }
+ }
+ }
+ // T
+ for(unsigned int i = 0 ; i < m_Layer1_StripYTMult ; ++i){
+ if(IsValidChannel("Y", m_EventData->GetCharissaLayer1StripYTDetectorNbr(i), m_EventData->GetCharissaLayer1StripYTStripNbr(i))){
+ m_PreTreatedData->SetCharissaLayer1StripYTDetectorNbr( m_EventData->GetCharissaLayer1StripYTDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripYTStripNbr( m_EventData->GetCharissaLayer1StripYTStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer1StripYTTime( fSi_Y_T(m_EventData , i) );
+ }
+ }
+
+ // X Layer2
+ // E
+ for(unsigned int i = 0 ; i < m_Layer2_StripXEMult ; ++i){
+ if( m_EventData->GetCharissaLayer2StripXEEnergy(i)>m_Si_X_E_RAW_Threshold && IsValidChannel("X", m_EventData->GetCharissaLayer2StripXEDetectorNbr(i), m_EventData->GetCharissaLayer1StripXEStripNbr(i)) ){
+ double Layer2_EX = fSi_X_E(m_EventData , i);
+ if( Layer2_EX > m_Si_X_E_Threshold ){
+ m_PreTreatedData->SetCharissaLayer2StripXEDetectorNbr( m_EventData->GetCharissaLayer2StripXEDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripXEStripNbr( m_EventData->GetCharissaLayer2StripXEStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripXEEnergy( Layer2_EX );
+ }
+
+ }
+ }
+
+ // T
+ for(unsigned int i = 0 ; i < m_Layer2_StripXTMult ; ++i){
+ if(IsValidChannel("X", m_EventData->GetCharissaLayer2StripXTDetectorNbr(i), m_EventData->GetCharissaLayer2StripXTStripNbr(i))){
+ m_PreTreatedData->SetCharissaLayer2StripXTDetectorNbr( m_EventData->GetCharissaLayer2StripXTDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripXTStripNbr( m_EventData->GetCharissaLayer2StripXTStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripXTTime( fSi_X_T(m_EventData , i) );
+ }
+ }
+
+
+ // Y Layer2
+ // E
+ for(unsigned int i = 0 ; i < m_Layer2_StripYEMult ; ++i){
+ if( m_EventData->GetCharissaLayer2StripYEEnergy(i)<m_Si_Y_E_RAW_Threshold && IsValidChannel("Y", m_EventData->GetCharissaLayer2StripYEDetectorNbr(i), m_EventData->GetCharissaLayer1StripYEStripNbr(i))){
+ double Layer2_EY = fSi_Y_E(m_EventData , i);
+ if( Layer2_EY >m_Si_Y_E_Threshold ){
+ m_PreTreatedData->SetCharissaLayer2StripYEDetectorNbr( m_EventData->GetCharissaLayer2StripYEDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripYEStripNbr( m_EventData->GetCharissaLayer2StripYEStripNbr(i) );
+ m_PreTreatedData->SetCharissaLayer2StripYEEnergy( Layer2_EY );
+ }
+ }
+ }
+
+
+
+
+ // CsI
+ // E
+ for(unsigned int i = 0 ; i < m_CsIEMult ; ++i){
+ if( m_EventData->GetCharissaCsIEEnergy(i)>m_CsI_E_RAW_Threshold && IsValidChannel("CsI", m_EventData->GetCharissaCsIEDetectorNbr(i), m_EventData->GetCharissaCsIECristalNbr(i))){
+ double ECsI = fCsI_E(m_EventData , i);
+ if( ECsI > m_CsI_E_Threshold ){
+ m_PreTreatedData->SetCharissaCsIEDetectorNbr( m_EventData->GetCharissaCsIEDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaCsIECristalNbr( m_EventData->GetCharissaCsIECristalNbr(i) );
+ m_PreTreatedData->SetCharissaCsIEEnergy( ECsI );
+ }
+ }
+ }
+
+ // T
+ for(unsigned int i = 0 ; i < m_CsITMult ; ++i){
+ if(IsValidChannel("CsI", m_EventData->GetCharissaCsITDetectorNbr(i), m_EventData->GetCharissaCsITCristalNbr(i))){
+ m_PreTreatedData->SetCharissaCsITDetectorNbr( m_EventData->GetCharissaCsITDetectorNbr(i) );
+ m_PreTreatedData->SetCharissaCsITCristalNbr( m_EventData->GetCharissaCsITCristalNbr(i) );
+ m_PreTreatedData->SetCharissaCsITTime( fCsI_T(m_EventData , i) );
+ }
+ }
+
+
+
+ return;
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+int TCharissaPhysics :: Layer1_CheckEvent(){
+ // Check the size of the different elements
+ if( m_PreTreatedData->GetCharissaLayer1StripXEMult() == m_PreTreatedData->GetCharissaLayer1StripYEMult() /*&& m_PreTreatedData->GetCharissaLayer1StripYEMult() == m_PreTreatedData->GetCharissaLayer1StripXTMult() && m_PreTreatedData->GetCharissaLayer1StripXTMult() == m_PreTreatedData->GetCharissaLayer1StripYTMult()*/ )
+ return 1 ; // Regular Event
+
+ else if( m_PreTreatedData->GetCharissaLayer1StripXEMult() == m_PreTreatedData->GetCharissaLayer1StripYEMult()+1 || m_PreTreatedData->GetCharissaLayer1StripXEMult() == m_PreTreatedData->GetCharissaLayer1StripYEMult()-1 )
+ return 2 ; // Pseudo Event, potentially interstrip
+
+ else
+ return -1 ; // Rejected Event
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+int TCharissaPhysics :: Layer2_CheckEvent(){
+ // Check the size of the different elements
+ if( m_PreTreatedData->GetCharissaLayer2StripXEMult() == m_PreTreatedData->GetCharissaLayer2StripYEMult() /*&& m_PreTreatedData->GetCharissaLayer1StripYEMult() == m_PreTreatedData->GetCharissaLayer1StripXTMult() && m_PreTreatedData->GetCharissaLayer1StripXTMult() == m_PreTreatedData->GetCharissaLayer1StripYTMult()*/ )
+ return 1 ; // Regular Event
+
+ else if( m_PreTreatedData->GetCharissaLayer2StripXEMult() == m_PreTreatedData->GetCharissaLayer2StripYEMult()+1 || m_PreTreatedData->GetCharissaLayer2StripXEMult() == m_PreTreatedData->GetCharissaLayer2StripYEMult()-1 )
+ return 2 ; // Pseudo Event, potentially interstrip
+
+ else
+ return -1 ; // Rejected Event
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+bool TCharissaPhysics :: ResolvePseudoEvent(){
+ return false;
+}
+
+///////////////////////////////////////////////////////////////////////////
+vector < TVector2 > TCharissaPhysics :: Layer1_Match_X_Y(){
+ vector < TVector2 > ArrayOfGoodCouple ;
+ m_Layer1_StripXEMult = m_PreTreatedData->GetCharissaLayer1StripXEMult();
+ m_Layer1_StripYEMult = m_PreTreatedData->GetCharissaLayer1StripXEMult();
+
+ // Prevent code from treating very high multiplicity Event
+ // Those event are not physical anyway and that improve speed.
+ if( m_Layer1_StripXEMult > m_MaximumStripMultiplicityAllowed || m_Layer1_StripYEMult > m_MaximumStripMultiplicityAllowed )
+ return ArrayOfGoodCouple;
+
+ for(unsigned int i = 0 ; i < m_Layer1_StripXEMult ; ++i){
+ for(unsigned int j = 0 ; j < m_Layer1_StripYEMult ; j++){
+ // if same detector check energy
+ if ( m_PreTreatedData->GetCharissaLayer1StripXEDetectorNbr(i) == m_PreTreatedData->GetCharissaLayer1StripYEDetectorNbr(j) ){
+ // Look if energy match
+ if( abs( (m_PreTreatedData->GetCharissaLayer1StripXEEnergy(i)-m_PreTreatedData->GetCharissaLayer1StripYEEnergy(j))/2. ) < m_StripEnergyMatchingNumberOfSigma*m_StripEnergyMatchingSigma ){
+ // Special Option, if the event is between two CsI cristal, it is rejected.
+ if(m_Ignore_not_matching_CsI){
+ bool check_validity=false;
+ for (unsigned int hh = 0 ; hh<m_NumberOfStrip ; ++hh ){
+ if( Match_Si_CsI(m_PreTreatedData->GetCharissaLayer1StripXEStripNbr(i), m_PreTreatedData->GetCharissaLayer1StripYEStripNbr(j) , hh+1) )
+ check_validity=true;
+ }
+
+ if(check_validity)
+ ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ }
+
+
+
+ // Regular case, keep the event
+ else ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ }
+ }
+ }
+ }
+
+ // Prevent to treat event with ambiguous matchin beetween X and Y
+ if( ArrayOfGoodCouple.size() > m_Layer1_StripXEMult ) ArrayOfGoodCouple.clear() ;
+
+ return ArrayOfGoodCouple;
+}
+
+///////////////////////////////////////////////////////////////////////////
+vector < TVector2 > TCharissaPhysics :: Layer2_Match_X_Y(){
+ vector < TVector2 > ArrayOfGoodCouple ;
+ m_Layer2_StripXEMult = m_PreTreatedData->GetCharissaLayer2StripXEMult();
+ m_Layer2_StripYEMult = m_PreTreatedData->GetCharissaLayer2StripXEMult();
+
+ // Prevent code from treating very high multiplicity Event
+ // Those event are not physical anyway and that improve speed.
+ if( m_Layer2_StripXEMult > m_MaximumStripMultiplicityAllowed || m_Layer2_StripYEMult > m_MaximumStripMultiplicityAllowed )
+ return ArrayOfGoodCouple;
+
+ for(unsigned int i = 0 ; i < m_Layer2_StripXEMult ; ++i){
+ for(unsigned int j = 0 ; j < m_Layer2_StripYEMult ; j++){
+ // if same detector check energy
+ if ( m_PreTreatedData->GetCharissaLayer2StripXEDetectorNbr(i) == m_PreTreatedData->GetCharissaLayer2StripYEDetectorNbr(j) ){
+ // Look if energy match
+ if( abs( (m_PreTreatedData->GetCharissaLayer2StripXEEnergy(i)-m_PreTreatedData->GetCharissaLayer2StripYEEnergy(j))/2. ) < m_StripEnergyMatchingNumberOfSigma*m_StripEnergyMatchingSigma ){
+ // Special Option, if the event is between two CsI cristal, it is rejected.
+ if(m_Ignore_not_matching_CsI){
+ bool check_validity=false;
+ for (unsigned int hh = 0 ; hh<m_NumberOfStrip ; ++hh ){
+ if( Match_Si_CsI(m_PreTreatedData->GetCharissaLayer2StripXEStripNbr(i), m_PreTreatedData->GetCharissaLayer2StripYEStripNbr(j) , hh+1) )
+ check_validity=true;
+ }
+
+ if(check_validity)
+ ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ }
+
+
+
+ // Regular case, keep the event
+ else ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ }
+ }
+ }
+ }
+
+ // Prevent to treat event with ambiguous matchin beetween X and Y
+ if( ArrayOfGoodCouple.size() > m_Layer2_StripXEMult ) ArrayOfGoodCouple.clear() ;
+
+ return ArrayOfGoodCouple;
+}
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////
+bool TCharissaPhysics :: IsValidChannel(const string DetectorType, const int telescope , const int channel){
+ if(DetectorType == "CsI")
+ return *(m_CsIChannelStatus[telescope-1].begin()+channel-1);
+
+ else if(DetectorType == "X")
+ return *(m_XChannelStatus[telescope-1].begin()+channel-1);
+
+ else if(DetectorType == "Y")
+ return *(m_YChannelStatus[telescope-1].begin()+channel-1);
+
+ else return false;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::ReadAnalysisConfig(){
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigCharissa.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigCharissa.dat found: Default parameters loaded for Analysis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameters for Analysis from ConfigCharissa.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigCharissaCHARISSA.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer,DataBuffer,whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ if (LineBuffer.compare(0, 11, "ConfigCharissa") == 0) ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus ) {
+
+ whatToDo="";
+ AnalysisConfigFile >> whatToDo;
+
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+
+ else if (whatToDo=="MAX_STRIP_MULTIPLICITY") {
+ AnalysisConfigFile >> DataBuffer;
+ m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() );
+ cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
+ }
+
+ else if (whatToDo=="STRIP_ENERGY_MATCHING_SIGMA") {
+ AnalysisConfigFile >> DataBuffer;
+ m_StripEnergyMatchingSigma = atof(DataBuffer.c_str() );
+ cout << "STRIP ENERGY MATCHING SIGMA " << m_StripEnergyMatchingSigma <<endl;
+ }
+
+ else if (whatToDo=="STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA") {
+ AnalysisConfigFile >> DataBuffer;
+ m_StripEnergyMatchingNumberOfSigma = atoi(DataBuffer.c_str() );
+ cout << "STRIP ENERGY MATCHING NUMBER OF SIGMA " << m_StripEnergyMatchingNumberOfSigma << endl;
+ }
+
+ else if (whatToDo== "DISABLE_ALL") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ int telescope = atoi(DataBuffer.substr(2,1).c_str());
+ vector< bool > ChannelStatus;
+ ChannelStatus.resize( m_NumberOfStrip,false);
+ m_XChannelStatus[telescope-1] = ChannelStatus;
+ m_YChannelStatus[telescope-1] = ChannelStatus;
+ ChannelStatus.resize(m_NumberOfStrip,false);
+ m_CsIChannelStatus[telescope-1] = ChannelStatus;
+ }
+
+ else if (whatToDo == "DISABLE_CHANNEL") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ int telescope = atoi(DataBuffer.substr(2,1).c_str());
+ int channel = -1;
+ if (DataBuffer.compare(3,4,"STRX") == 0) {
+ channel = atoi(DataBuffer.substr(7).c_str());
+ *(m_XChannelStatus[telescope-1].begin()+channel-1) = false;
+ }
+
+ else if (DataBuffer.compare(3,4,"STRY") == 0) {
+ channel = atoi(DataBuffer.substr(7).c_str());
+ *(m_YChannelStatus[telescope-1].begin()+channel-1) = false;
+ }
+
+
+ else if (DataBuffer.compare(3,3,"CSI") == 0) {
+ channel = atoi(DataBuffer.substr(6).c_str());
+ *(m_CsIChannelStatus[telescope-1].begin()+channel-1) = false;
+ }
+
+ else cout << "Warning: detector type for CharissaCHARISSA unknown!" << endl;
+
+ }
+
+ else if (whatToDo=="TAKE_E_Y") {
+ m_Take_E_Y = true;
+ cout << whatToDo << endl;
+ }
+
+ else if (whatToDo=="TAKE_T_Y") {
+ m_Take_T_Y = true;
+ cout << whatToDo << endl;
+ }
+
+ else if (whatToDo=="TAKE_E_X") {
+ m_Take_E_Y = false;
+ cout << whatToDo << endl;
+ }
+
+ else if (whatToDo=="TAKE_T_X") {
+ m_Take_T_Y = false;
+ cout << whatToDo << endl;
+ }
+
+ else if (whatToDo== "IGNORE_NOT_MATCHING_CSI") {
+ m_Ignore_not_matching_CsI = true;
+ cout << whatToDo << endl;
+ }
+
+ else if (whatToDo=="CSI_SIZE") {
+ AnalysisConfigFile >> DataBuffer;
+ m_CsI_Size = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_CsI_Size << endl;
+ }
+
+
+ else if (whatToDo=="SI_X_E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_Si_X_E_RAW_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_Si_X_E_RAW_Threshold << endl;
+ }
+
+ else if (whatToDo=="SI_Y_E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_Si_Y_E_RAW_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_Si_Y_E_RAW_Threshold << endl;
+ }
+
+
+ else if (whatToDo== "CSI_E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_CsI_E_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_CsI_E_Threshold << endl;
+ }
+
+ else if (whatToDo=="SI_X_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_Si_X_E_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_Si_X_E_Threshold << endl;
+ }
+
+ else if (whatToDo== "SI_Y_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_Si_Y_E_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_Si_Y_E_Threshold << endl;
+ }
+
+
+ else if (whatToDo=="CSI_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_CsI_E_Threshold = atoi(DataBuffer.c_str());
+ cout << whatToDo << " " << m_CsI_E_Threshold << endl;
+ }
+
+ else {
+ ReadingStatus = false;
+ }
+
+ }
+ }
+}
+
+
+
+
+///////////////////////////////////////////////////////////////////////////
+bool TCharissaPhysics :: Match_Si_CsI(int X, int Y , int CristalNbr){
+
+ if( abs(m_CsI_MatchingX[CristalNbr-1] - X) < m_CsI_Size/2.&&
+ abs(m_CsI_MatchingY[CristalNbr-1] - Y) < m_CsI_Size/2.)
+
+ return true ;
+
+ else return false;
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::Clear(){
+ EventMultiplicity= 0 ;
+
+ Layer1_TelescopeNumber.clear();
+ EventType.clear();
+ TotalEnergy.clear();
+
+ // Si X
+ Layer1_Si_E.clear();
+ Layer1_Si_T.clear();
+ Layer1_Si_X.clear();
+ Layer1_Si_Y.clear();
+ Layer2_Si_E.clear();
+ Layer2_Si_T.clear();
+ Layer2_Si_X.clear();
+ Layer2_Si_Y.clear();
+
+
+ // CsI
+ CsI_E.clear();
+ CsI_T.clear();
+ CsI_N.clear();
+
+ Layer1_Si_EX.clear();
+ Layer1_Si_TX.clear();
+ Layer1_Si_EY.clear();
+ Layer1_Si_TY.clear();
+ Layer1_TelescopeNumber_X.clear();
+ Layer1_TelescopeNumber_Y.clear();
+ Layer2_Si_EX.clear();
+ Layer2_Si_TX.clear();
+ Layer2_Si_EY.clear();
+ Layer2_Si_TY.clear();
+ Layer2_TelescopeNumber_X.clear();
+ Layer2_TelescopeNumber_Y.clear();
+}
+///////////////////////////////////////////////////////////////////////////
+
+void TCharissaPhysics::ReadCalibrationRun(){
+ m_Layer1_StripXEMult = m_EventData->GetCharissaLayer1StripXEMult();
+ m_Layer1_StripYEMult = m_EventData->GetCharissaLayer1StripYEMult();
+ m_Layer1_StripXTMult = m_EventData->GetCharissaLayer1StripXTMult();
+ m_Layer1_StripYTMult = m_EventData->GetCharissaLayer1StripYTMult();
+ m_Layer2_StripXEMult = m_EventData->GetCharissaLayer2StripXEMult();
+ m_Layer2_StripYEMult = m_EventData->GetCharissaLayer2StripYEMult();
+ m_Layer2_StripXTMult = m_EventData->GetCharissaLayer2StripXTMult();
+ m_Layer2_StripYTMult = m_EventData->GetCharissaLayer2StripYTMult();
+ m_CsIEMult = m_EventData->GetCharissaCsIEMult();
+ m_CsITMult = m_EventData->GetCharissaCsITMult();
+
+ // X Layer 1
+ // E
+ for(unsigned int i = 0 ; i < m_Layer1_StripXEMult;++i){
+ Layer1_TelescopeNumber_X.push_back(m_EventData->GetCharissaLayer1StripXEDetectorNbr(i));
+ Layer1_Si_EX.push_back( fSi_X_E( m_EventData , i) );
+ Layer1_Si_X.push_back(m_EventData->GetCharissaLayer1StripXEStripNbr(i));
+ }
+ // T
+ for(unsigned int i = 0 ; i < m_Layer1_StripXTMult;++i){
+ Layer1_TelescopeNumber_X.push_back(m_EventData->GetCharissaLayer1StripXTDetectorNbr(i));
+ Layer1_Si_TX.push_back( fSi_X_T( m_EventData , i) );
+ Layer1_Si_X.push_back(m_EventData->GetCharissaLayer1StripXTStripNbr(i));
+ }
+
+ // Y Layer 1
+ // E
+ for(unsigned int i = 0 ; i < m_Layer1_StripYEMult;++i){
+ Layer1_TelescopeNumber_Y.push_back(m_EventData->GetCharissaLayer1StripYEDetectorNbr(i));
+ Layer1_Si_EY.push_back( fSi_Y_E( m_EventData , i) );
+ Layer1_Si_Y.push_back(m_EventData->GetCharissaLayer1StripYEStripNbr(i));
+ }
+
+ // T
+ for(unsigned int i = 0 ; i < m_Layer1_StripYTMult;++i){
+ Layer1_TelescopeNumber_Y.push_back(m_EventData->GetCharissaLayer1StripYTDetectorNbr(i));
+ Layer1_Si_TY.push_back( fSi_Y_T( m_EventData , i) );
+ Layer1_Si_Y.push_back(m_EventData->GetCharissaLayer1StripYTStripNbr(i));
+ }
+
+
+ // X Layer 2
+ // E
+ for(unsigned int i = 0 ; i < m_Layer2_StripXEMult;++i){
+ Layer2_TelescopeNumber_X.push_back(m_EventData->GetCharissaLayer2StripXEDetectorNbr(i));
+ Layer2_Si_EX.push_back( fSi_X_E( m_EventData , i) );
+ Layer2_Si_X.push_back(m_EventData->GetCharissaLayer2StripXEStripNbr(i));
+ }
+ // T
+ for(unsigned int i = 0 ; i < m_Layer2_StripXTMult;++i){
+ Layer2_TelescopeNumber_X.push_back(m_EventData->GetCharissaLayer2StripXTDetectorNbr(i));
+ Layer2_Si_TX.push_back( fSi_X_T( m_EventData , i) );
+ Layer2_Si_X.push_back(m_EventData->GetCharissaLayer2StripXTStripNbr(i));
+ }
+
+ // Y Layer 2
+ // E
+ for(unsigned int i = 0 ; i < m_Layer2_StripYEMult;++i){
+ Layer2_TelescopeNumber_Y.push_back(m_EventData->GetCharissaLayer2StripYEDetectorNbr(i));
+ Layer2_Si_EY.push_back( fSi_Y_E( m_EventData , i) );
+ Layer2_Si_Y.push_back(m_EventData->GetCharissaLayer2StripYEStripNbr(i));
+ }
+
+ // T
+ for(unsigned int i = 0 ; i < m_Layer2_StripYTMult;++i){
+ Layer2_TelescopeNumber_Y.push_back(m_EventData->GetCharissaLayer2StripYTDetectorNbr(i));
+ Layer2_Si_TY.push_back( fSi_Y_T( m_EventData , i) );
+ Layer2_Si_Y.push_back(m_EventData->GetCharissaLayer2StripYTStripNbr(i));
+ }
+
+
+ //CsI Energy
+ for(unsigned int j = 0 ; j < m_CsIEMult ; ++j){
+ CsI_N.push_back(m_EventData->GetCharissaCsIECristalNbr(j)) ;
+ CsI_E.push_back(fCsI_E(m_EventData , j)) ;
+ }
+
+ //CsI Time
+ for(unsigned int j = 0 ; j < m_CsITMult ; ++j){
+ CsI_T.push_back(fCsI_T(m_EventData , j)) ;
+ }
+
+}
+
+//// Innherited from VDetector Class ////
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::ReadConfiguration(string Path){
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+ string LineBuffer ;
+ string DataBuffer ;
+
+ // A:X1_Y1 --> X:1 Y:1
+ // B:X m_NumberOfStrip_Y1 --> X: m_NumberOfStrip Y:1
+ // C:X1_Y m_NumberOfStrip --> X:1 Y: m_NumberOfStrip
+ // D:X m_NumberOfStrip_Y m_NumberOfStrip --> X: m_NumberOfStrip Y: m_NumberOfStrip
+
+ double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz;
+ TVector3 A , B , C , D;
+ double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0;
+
+ bool check_A = false ;
+ bool check_C = false ;
+ bool check_B = false ;
+ bool check_D = false ;
+
+ bool check_Theta = false ;
+ bool check_Phi = false ;
+ bool check_R = false ;
+ bool check_beta = false ;
+
+ bool ReadingStatus = false ;
+
+
+ while (!ConfigFile.eof())
+ {
+
+ getline(ConfigFile, LineBuffer);
+
+ // If line is a Start Up CharissaCHARISSA bloc, Reading toggle to true
+ if (LineBuffer.compare(0, 17, "CharissaTelescope")==0)
+ {
+ cout << "///" << endl ;
+ cout << "Telescope found: " << endl ;
+ ReadingStatus = true ;
+
+ }
+
+ // Else don't toggle to Reading Block Status
+ else ReadingStatus = false ;
+
+ // Reading Block
+ while(ReadingStatus)
+ {
+
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+ if(DataBuffer.compare(0, 1, "%") == 0) {
+ ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
+
+ }
+
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer=="CharissaTelescope") {
+ cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
+ ReadingStatus = false ;
+ }
+
+ // Position method
+
+ else if (DataBuffer=="X1_Y1=") {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+
+ A = TVector3(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+
+ }
+
+
+ else if (DataBuffer=="X16_Y1=") {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+
+ B = TVector3(Bx, By, Bz);
+ cout << "X16 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+
+ }
+
+
+ else if (DataBuffer=="X1_Y16=") {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+
+ C = TVector3(Cx, Cy, Cz);
+ cout << "X1 Y16 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
+
+ }
+
+ else if (DataBuffer=="X16_Y16=") {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+
+ D = TVector3(Dx, Dy, Dz);
+ cout << "X16 Y16 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+
+ }
+
+ // End Position Method
+
+ // Angle method
+ else if (DataBuffer=="THETA=") {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ cout << "Theta: " << Theta << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer=="PHI=") {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ cout << "Phi: " << Phi << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer=="R=") {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ cout << "R: " << R << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer=="BETA=") {
+ check_beta = true;
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ cout << "Beta: " << beta_u << " " << beta_v << " " << beta_w << endl ;
+
+ }
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) )
+ {
+ ReadingStatus = false;
+
+ ///Add The previously define telescope
+ //With position method
+ if ( check_A && check_B && check_C && check_D )
+ {
+ AddTelescope( A,
+ B,
+ C,
+ D) ;
+ }
+
+ //with angle method
+ else if ( check_Theta && check_Phi && check_R && check_beta )
+ {
+ AddTelescope( Theta,
+ Phi,
+ R,
+ beta_u,
+ beta_v,
+ beta_w);
+ }
+
+ check_A = false ;
+ check_B = false ;
+ check_C = false ;
+ check_D = false ;
+
+ check_Theta = false ;
+ check_Phi = false ;
+ check_R = false ;
+ check_beta = false ;
+ }
+ }
+ }
+
+ InitializeStandardParameter();
+ ReadAnalysisConfig();
+
+ cout << endl << "/////////////////////////////" << endl << endl;
+
+}
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::InitSpectra(){
+ m_Spectra = new TCharissaSpectra(m_NumberOfTelescope);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::FillSpectra(){
+ m_Spectra -> FillRawSpectra(m_EventData);
+ m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
+// m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
+}
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::CheckSpectra(){
+ m_Spectra->CheckSpectra();
+}
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::ClearSpectra(){
+ // To be done
+}
+///////////////////////////////////////////////////////////////////////////
+map< vector<TString> , TH1*> TCharissaPhysics::GetSpectra() {
+return m_Spectra->GetMapHisto();
+}
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::AddParameterToCalibrationManager()
+{
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ for(int i = 0 ; i < m_NumberOfTelescope ; ++i){
+
+ for( int j = 0 ; j < m_NumberOfStrip ; ++j){
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_Si_X"+itoa(j+1)+"_E","CHARISSA_T"+itoa(i+1)+"_Si_X"+itoa(j+1)+"_E") ;
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_Si_Y"+itoa(j+1)+"_E","CHARISSA_T"+itoa(i+1)+"_Si_Y"+itoa(j+1)+"_E") ;
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_Si_X"+itoa(j+1)+"_T","CHARISSA_T"+itoa(i+1)+"_Si_X"+itoa(j+1)+"_T") ;
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_Si_Y"+itoa(j+1)+"_T","CHARISSA_T"+itoa(i+1)+"_Si_Y"+itoa(j+1)+"_T") ;
+ }
+
+ for( int j = 0 ; j < m_NumberOfStrip ; ++j){
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_E","CHARISSA_T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_E") ;
+ Cal->AddParameter("CHARISSA", "T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_T","CHARISSA_T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_T") ;
+ }
+ }
+
+ return;
+
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::InitializeRootInputRaw()
+{
+ TChain* inputChain = RootInput::getInstance()->GetChain() ;
+ inputChain->SetBranchStatus( "CHARISSA" , true ) ;
+ inputChain->SetBranchStatus( "fMM_*" , true ) ;
+ inputChain->SetBranchAddress( "CHARISSA" , &m_EventData ) ;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::InitializeRootInputPhysics()
+{
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus( "CHARISSA" , true );
+ inputChain->SetBranchStatus( "EventMultiplicity" , true );
+ inputChain->SetBranchStatus( "EventType" , true );
+ inputChain->SetBranchStatus( "Layer1_TelescopeNumber" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_E" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_T" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_X" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_Y" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_EX" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_TX" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_EY" , true );
+ inputChain->SetBranchStatus( "Layer1_Si_TY" , true );
+ inputChain->SetBranchStatus( "Layer1_TelescopeNumber_X" , true );
+ inputChain->SetBranchStatus( "Layer1_TelescopeNumber_Y" , true );
+ inputChain->SetBranchStatus( "Layer2_TelescopeNumber" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_E" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_T" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_X" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_Y" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_EX" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_TX" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_EY" , true );
+ inputChain->SetBranchStatus( "Layer2_Si_TY" , true );
+ inputChain->SetBranchStatus( "Layer2_TelescopeNumber_X" , true );
+ inputChain->SetBranchStatus( "Layer2_TelescopeNumber_Y" , true );
+ inputChain->SetBranchStatus( "CsI_E" , true );
+ inputChain->SetBranchStatus( "CsI_T" , true );
+ inputChain->SetBranchStatus( "CsI_N" , true );
+ inputChain->SetBranchStatus( "TotalEnergy" , true );
+ inputChain->SetBranchAddress( "CHARISSA" , &m_EventPhysics);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::InitializeRootOutput()
+{
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch( "CHARISSA" , "TCharissaPhysics" , &m_EventPhysics );
+}
+
+
+///// Specific to CHARISSAArray ////
+/////////////////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::AddTelescope( TVector3 C_X1_Y1,
+ TVector3 C_X16_Y1,
+ TVector3 C_X1_Y16,
+ TVector3 C_X16_Y16)
+{
+ // To avoid warning
+ C_X16_Y16 *= 1;
+
+ m_NumberOfTelescope++;
+
+ // Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U = C_X16_Y1 - C_X1_Y1 ;
+ double Ushift = (U.Mag()-98)/2.;
+ U = U.Unit();
+ // Vector V on Telescope Face (parallele to X Strip)
+ TVector3 V = C_X1_Y16 - C_X1_Y1 ;
+ double Vshift = (V.Mag() -98)/2. ;
+ V = V.Unit() ;
+
+ // Position Vector of Strip Center
+ TVector3 StripCenter = TVector3(0,0,0);
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1;
+
+ // Geometry Parameter
+ double Face = 98; //mm
+ double NumberOfStrip = m_NumberOfStrip;
+ double StripPitch = Face/NumberOfStrip ; //mm
+ // Buffer object to fill Position Array
+ vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
+
+ vector< vector< double > > OneTelescopeStripPositionX ;
+ vector< vector< double > > OneTelescopeStripPositionY ;
+ vector< vector< double > > OneTelescopeStripPositionZ ;
+
+ // Moving StripCenter to 1.1 corner:
+ Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2.) ;
+ Strip_1_1+= U*Ushift+V*Vshift ;
+
+ for( int i = 0 ; i < m_NumberOfStrip ; ++i )
+ {
+ lineX.clear() ;
+ lineY.clear() ;
+ lineZ.clear() ;
+
+ for( int j = 0 ; j < m_NumberOfStrip ; ++j )
+ {
+ StripCenter = Strip_1_1 + StripPitch*( i*U + j*V );
+ lineX.push_back( StripCenter.X() );
+ lineY.push_back( StripCenter.Y() );
+ lineZ.push_back( StripCenter.Z() );
+ }
+
+ OneTelescopeStripPositionX.push_back(lineX);
+ OneTelescopeStripPositionY.push_back(lineY);
+ OneTelescopeStripPositionZ.push_back(lineZ);
+
+ }
+ m_StripPositionX.push_back( OneTelescopeStripPositionX ) ;
+ m_StripPositionY.push_back( OneTelescopeStripPositionY ) ;
+ m_StripPositionZ.push_back( OneTelescopeStripPositionZ ) ;
+
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::InitializeStandardParameter()
+{
+ // Enable all channel
+ vector< bool > ChannelStatus;
+ m_XChannelStatus.clear() ;
+ m_YChannelStatus.clear() ;
+ m_CsIChannelStatus.clear() ;
+
+ ChannelStatus.resize(m_NumberOfStrip,true);
+ for(int i = 0 ; i < m_NumberOfTelescope ; ++i)
+ {
+ m_XChannelStatus[i] = ChannelStatus;
+ m_YChannelStatus[i] = ChannelStatus;
+ }
+
+ ChannelStatus.resize(m_NumberOfStrip,true);
+ for(int i = 0 ; i < m_NumberOfTelescope ; ++i)
+ {
+ m_CsIChannelStatus[i] = ChannelStatus;
+ }
+
+
+ m_MaximumStripMultiplicityAllowed = m_NumberOfTelescope ;
+
+ return;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////
+void TCharissaPhysics::AddTelescope( double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w)
+{
+
+
+ m_NumberOfTelescope++;
+
+ double Pi = 3.141592654 ;
+
+ // convert from degree to radian:
+ theta = theta * Pi/180. ;
+ phi = phi * Pi/180. ;
+
+ //Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
+ TVector3 U ;
+ //Vector V on Telescope Face (parallele to X Strip)
+ TVector3 V ;
+ //Vector W normal to Telescope Face (pointing CsI)
+ TVector3 W ;
+ //Vector position of Telescope Face center
+ TVector3 C ;
+
+ C = TVector3 ( distance * sin(theta) * cos(phi) ,
+ distance * sin(theta) * sin(phi) ,
+ distance * cos(theta) );
+
+ TVector3 P = TVector3( cos(theta ) * cos(phi) ,
+ cos(theta ) * sin(phi) ,
+ -sin(theta) );
+
+ W = C.Unit() ;
+ U = W .Cross ( P ) ;
+ V = W .Cross ( U );
+
+ U = U.Unit();
+ V = V.Unit();
+
+ U.Rotate( beta_u * Pi/180. , U ) ;
+ V.Rotate( beta_u * Pi/180. , U ) ;
+
+ U.Rotate( beta_v * Pi/180. , V ) ;
+ V.Rotate( beta_v * Pi/180. , V ) ;
+
+ U.Rotate( beta_w * Pi/180. , W ) ;
+ V.Rotate( beta_w * Pi/180. , W ) ;
+
+ double Face = 98 ; //mm
+ double NumberOfStrip = m_NumberOfStrip ;
+ double StripPitch = Face/NumberOfStrip ; //mm
+
+ vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
+
+ vector< vector< double > > OneTelescopeStripPositionX ;
+ vector< vector< double > > OneTelescopeStripPositionY ;
+ vector< vector< double > > OneTelescopeStripPositionZ ;
+
+ double X , Y , Z ;
+
+ //Moving C to the 1.1 corner:
+ C.SetX( C.X() - ( Face/2 - StripPitch/2 ) * ( V.X() + U.X() ) ) ;
+ C.SetY( C.Y() - ( Face/2 - StripPitch/2 ) * ( V.Y() + U.Y() ) ) ;
+ C.SetZ( C.Z() - ( Face/2 - StripPitch/2 ) * ( V.Z() + U.Z() ) ) ;
+
+ for( int i = 0 ; i < m_NumberOfStrip ; ++i )
+ {
+
+ lineX.clear() ;
+ lineY.clear() ;
+ lineZ.clear() ;
+
+ for( int j = 0 ; j < m_NumberOfStrip ; ++j )
+ {
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j ) ;
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j ) ;
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j ) ;
+
+ lineX.push_back(X) ;
+ lineY.push_back(Y) ;
+ lineZ.push_back(Z) ;
+
+ }
+
+ OneTelescopeStripPositionX.push_back(lineX) ;
+ OneTelescopeStripPositionY.push_back(lineY) ;
+ OneTelescopeStripPositionZ.push_back(lineZ) ;
+
+ }
+ m_StripPositionX.push_back( OneTelescopeStripPositionX ) ;
+ m_StripPositionY.push_back( OneTelescopeStripPositionY ) ;
+ m_StripPositionZ.push_back( OneTelescopeStripPositionZ ) ;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////
+TVector3 TCharissaPhysics::GetPositionOfInteraction(const int i) const
+{
+ TVector3 Position = TVector3 ( GetStripPositionX( Layer1_TelescopeNumber[i] , Layer1_Si_X[i] , Layer1_Si_Y[i] ) ,
+ GetStripPositionY( Layer1_TelescopeNumber[i] , Layer1_Si_X[i] , Layer1_Si_Y[i] ) ,
+ GetStripPositionZ( Layer1_TelescopeNumber[i] , Layer1_Si_X[i] , Layer1_Si_Y[i] ) ) ;
+
+ return(Position) ;
+
+}
+
+TVector3 TCharissaPhysics::GetTelescopeNormal( const int i) const
+{
+ TVector3 U = TVector3 ( GetStripPositionX( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) ,
+ GetStripPositionY( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) ,
+ GetStripPositionZ( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) )
+
+ -TVector3 ( GetStripPositionX( Layer1_TelescopeNumber[i] , 1 , 1 ) ,
+ GetStripPositionY( Layer1_TelescopeNumber[i] , 1 , 1 ) ,
+ GetStripPositionZ( Layer1_TelescopeNumber[i] , 1 , 1 ) );
+
+ TVector3 V = TVector3 ( GetStripPositionX( Layer1_TelescopeNumber[i] , m_NumberOfStrip , m_NumberOfStrip ) ,
+ GetStripPositionY( Layer1_TelescopeNumber[i] , m_NumberOfStrip , m_NumberOfStrip ) ,
+ GetStripPositionZ( Layer1_TelescopeNumber[i] , m_NumberOfStrip , m_NumberOfStrip ) )
+
+ -TVector3 ( GetStripPositionX( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) ,
+ GetStripPositionY( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) ,
+ GetStripPositionZ( Layer1_TelescopeNumber[i] , m_NumberOfStrip , 1 ) );
+
+ TVector3 Normal = U.Cross(V);
+
+ return(Normal.Unit()) ;
+}
+
+///////////////////////////////////////////////////////////////////////////
+namespace CHARISSA_LOCAL
+{
+
+ // tranform an integer to a string
+ string itoa(int value)
+ {
+char buffer [33];
+sprintf(buffer,"%d",value);
+return buffer;
+ }
+
+ // DSSD
+ // X
+ double fSi_X_E(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaLayer1StripXEDetectorNbr(i) ) + "_Si_X" + itoa( m_EventData->GetCharissaLayer1StripXEStripNbr(i) ) + "_E",
+ m_EventData->GetCharissaLayer1StripXEEnergy(i) );
+ }
+
+ double fSi_X_T(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaLayer1StripXTDetectorNbr(i) ) + "_Si_X" + itoa( m_EventData->GetCharissaLayer1StripXTStripNbr(i) ) +"_T",
+ m_EventData->GetCharissaLayer1StripXTTime(i) );
+ }
+
+ // Y
+ double fSi_Y_E(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaLayer1StripYEDetectorNbr(i) ) + "_Si_Y" + itoa( m_EventData->GetCharissaLayer1StripYEStripNbr(i) ) +"_E",
+ m_EventData->GetCharissaLayer1StripYEEnergy(i) );
+ }
+
+ double fSi_Y_T(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaLayer1StripYTDetectorNbr(i) ) + "_Si_Y" + itoa( m_EventData->GetCharissaLayer1StripYTStripNbr(i) ) +"_T",
+ m_EventData->GetCharissaLayer1StripYTTime(i) );
+ }
+
+ // CsI
+ double fCsI_E(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaCsIEDetectorNbr(i) ) + "_CsI" + itoa( m_EventData->GetCharissaCsIECristalNbr(i) ) +"_E",
+ m_EventData->GetCharissaCsIEEnergy(i) );
+ }
+
+ double fCsI_T(const TCharissaData* m_EventData , const int i)
+ {
+return CalibrationManager::getInstance()->ApplyCalibration( "CHARISSA/T" + itoa( m_EventData->GetCharissaCsITDetectorNbr(i) ) + "_CsI" + itoa( m_EventData->GetCharissaCsITCristalNbr(i) ) +"_T",
+ m_EventData->GetCharissaCsITTime(i) );
+ }
+
+}
+
+
diff --git a/NPLib/Charissa/TCharissaPhysics.h b/NPLib/Charissa/TCharissaPhysics.h
new file mode 100755
index 0000000000000000000000000000000000000000..8e5bcb41c477f1c07ef3184422eddee1edac1314
--- /dev/null
+++ b/NPLib/Charissa/TCharissaPhysics.h
@@ -0,0 +1,310 @@
+#ifndef TCHARISSAPHYSICS_H
+#define TCHARISSAPHYSICS_H
+/*****************************************************************************
+ * 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 : December 2013 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold charissa treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+// STL
+#include <vector>
+#include <map>
+// NPL
+#include "TCharissaData.h"
+#include "TCharissaSpectra.h"
+#include "../include/CalibrationManager.h"
+#include "../include/VDetector.h"
+// ROOT
+#include "TVector2.h"
+#include "TVector3.h"
+#include "TObject.h"
+#include "TH1.h"
+
+using namespace std ;
+
+// Forward Declaration
+class TCharissaSpectra;
+
+class TCharissaPhysics : public TObject, public NPA::VDetector{
+ public:
+ TCharissaPhysics();
+ ~TCharissaPhysics();
+
+ public:
+ void Clear();
+ void Clear(const Option_t*) {};
+
+ public:
+ vector < TVector2 > Layer1_Match_X_Y();
+ vector < TVector2 > Layer2_Match_X_Y();
+ bool Match_Si_CsI(int X, int Y , int CristalNbr);
+ bool ResolvePseudoEvent();
+ int Layer1_CheckEvent();
+ int Layer2_CheckEvent();
+
+ public:
+
+ // Provide Physical Multiplicity
+ Int_t EventMultiplicity;
+
+ // Provide a Classification of Event
+ vector<int> EventType ;
+
+ // Telescope
+ vector<int> Layer1_TelescopeNumber ;
+ vector<int> Layer2_TelescopeNumber ;
+
+ // Si
+ vector<double> Layer1_Si_E ;
+ vector<double> Layer1_Si_T ;
+ vector<int> Layer1_Si_X ;
+ vector<int> Layer1_Si_Y ;
+
+ vector<double> Layer2_Si_E ;
+ vector<double> Layer2_Si_T ;
+ vector<int> Layer2_Si_X ;
+ vector<int> Layer2_Si_Y ;
+
+ // Use for checking purpose
+ vector<double> Layer1_Si_EX ;
+ vector<double> Layer1_Si_TX ;
+ vector<double> Layer1_Si_EY ;
+ vector<double> Layer1_Si_TY ;
+ vector<int> Layer1_TelescopeNumber_X ;
+ vector<int> Layer1_TelescopeNumber_Y ;
+
+ vector<double> Layer2_Si_EX ;
+ vector<double> Layer2_Si_TX ;
+ vector<double> Layer2_Si_EY ;
+ vector<double> Layer2_Si_TY ;
+ vector<int> Layer2_TelescopeNumber_X ;
+ vector<int> Layer2_TelescopeNumber_Y ;
+ // Si(Li)
+
+
+ // CsI
+ vector<double> CsI_E ;
+ vector<double> CsI_T ;
+ vector<int> CsI_N ;
+
+ // Physical Value
+ vector<double> TotalEnergy ;
+
+
+ public: // Innherited from VDetector Class
+
+ // 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() ;
+
+ // Same as above but for online analysis
+ void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
+
+ // Method related to the TSpectra classes, aimed at providing a framework for online applications
+ // Instantiate the Spectra class and the histogramm throught it
+ void InitSpectra();
+ // Fill the spectra hold by the spectra class
+ void FillSpectra();
+ // Used for Online mainly, perform check on the histo and for example change their color if issues are found
+ void CheckSpectra();
+ // Used for Online only, clear all the spectra hold by the Spectra class
+ void ClearSpectra();
+
+ public: // Specific to CHARISSA Array
+
+ // Clear The PreTeated object
+ void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+
+ // Remove bad channel, calibrate the data and apply threshold
+ void PreTreat();
+
+ // Return false if the channel is disabled by user
+ // Frist argument is either "X","Y","CsI"
+ bool IsValidChannel(const string DetectorType, const int telescope , const int channel);
+
+ // Initialize the standard parameter for analysis
+ // ie: all channel enable, maximum multiplicity for strip = number of telescope
+ void InitializeStandardParameter();
+
+ // Read the user configuration file; if no file found, load standard one
+ void ReadAnalysisConfig();
+
+ // Add a Telescope using Corner Coordinate information
+ void AddTelescope( TVector3 C_X1_Y1,
+ TVector3 C_X128_Y1,
+ TVector3 C_X1_Y128,
+ TVector3 C_X128_Y128);
+
+ // Add a Telescope using R Theta Phi of Si center information
+ void AddTelescope( double theta,
+ double phi,
+ double distance,
+ double beta_u,
+ double beta_v,
+ double beta_w);
+
+ // Use for reading Calibration Run, very simple methods; only apply calibration, no condition
+ void ReadCalibrationRun();
+
+ // Give and external TMustData object to TCharissaPhysics. Needed for online analysis for example.
+ void SetRawDataPointer(TCharissaData* rawDataPointer) {m_EventData = rawDataPointer;}
+ // Retrieve raw and pre-treated data
+ TCharissaData* GetRawData() const {return m_EventData;}
+ TCharissaData* GetPreTreatedData() const {return m_PreTreatedData;}
+
+ // Use to access the strip position
+ double GetStripPositionX( const int N , const int X , const int Y ) const{ return m_StripPositionX[N-1][X-1][Y-1] ; } ;
+ double GetStripPositionY( const int N , const int X , const int Y ) const{ return m_StripPositionY[N-1][X-1][Y-1] ; } ;
+ double GetStripPositionZ( const int N , const int X , const int Y ) const{ return m_StripPositionZ[N-1][X-1][Y-1] ; } ;
+
+ double GetNumberOfTelescope() const { return m_NumberOfTelescope ; };
+
+ // To be called after a build Physical Event
+ int GetEventMultiplicity() const { return EventMultiplicity; };
+
+ double GetEnergyDeposit(const int i) const{ return TotalEnergy[i] ;};
+
+ TVector3 GetPositionOfInteraction(const int i) const;
+ TVector3 GetTelescopeNormal(const int i) const;
+
+ private: // Parameter used in the analysis
+
+ // By default take EX and TY.
+ bool m_Take_E_Y;//!
+ bool m_Take_T_Y;//!
+
+ // Size Container to be used in the different loop of the analysis (value must be given locally)
+ unsigned int m_Layer1_StripXEMult;//!
+ unsigned int m_Layer1_StripYEMult;//!
+ unsigned int m_Layer1_StripXTMult;//!
+ unsigned int m_Layer1_StripYTMult;//!
+ unsigned int m_Layer2_StripXEMult;//!
+ unsigned int m_Layer2_StripYEMult;//!
+ unsigned int m_Layer2_StripXTMult;//!
+ unsigned int m_Layer2_StripYTMult;//!
+ unsigned int m_CsIEMult;//!
+ unsigned int m_CsITMult;//!
+
+
+ // Event over this value after pre-treatment are not treated / avoid long treatment time on spurious event
+ unsigned int m_MaximumStripMultiplicityAllowed ;//!
+ // Give the allowance in percent of the difference in energy between X and Y
+ double m_StripEnergyMatchingSigma ; //!
+ double m_StripEnergyMatchingNumberOfSigma ; //!
+
+ // Raw Threshold
+ int m_Si_X_E_RAW_Threshold ;//!
+ int m_Si_Y_E_RAW_Threshold ;//!
+ int m_CsI_E_RAW_Threshold ;//!
+
+ // Calibrated Threshold
+ double m_Si_X_E_Threshold ;//!
+ double m_Si_Y_E_Threshold ;//!
+ double m_CsI_E_Threshold ;//!
+
+ // Geometric Matching
+ // size in strip of a cristal
+ int m_CsI_Size;//!
+ // center position of the cristal on X
+ vector< int > m_CsI_MatchingX;//!
+ // center position of the cristal on X
+ vector< int > m_CsI_MatchingY;//!
+
+ // If set to true, all event that do not come in front of a cristal will be ignore all time (crossing or not),
+ // Warning, this option reduce statistic, however it help eliminating unrealevent event that cross the DSSD
+ // And go between pad or cristal.
+ bool m_Ignore_not_matching_CsI;//!
+
+ private: // Root Input and Output tree classes
+
+ TCharissaData* m_EventData;//!
+ TCharissaData* m_PreTreatedData;//!
+ TCharissaPhysics* m_EventPhysics;//!
+
+ private: // Map of activated channel
+ map< int, vector<bool> > m_XChannelStatus;//!
+ map< int, vector<bool> > m_YChannelStatus;//!
+ map< int, vector<bool> > m_CsIChannelStatus;//!
+
+ private: // Spatial Position of Strip Calculated on bases of detector position
+
+ unsigned int m_NumberOfTelescope;//!
+ unsigned int m_NumberOfStrip;//!
+
+ vector< vector < vector < double > > > m_StripPositionX;//!
+ vector< vector < vector < double > > > m_StripPositionY;//!
+ vector< vector < vector < double > > > m_StripPositionZ;//!
+
+ private: // Spectra Class
+ TCharissaSpectra* m_Spectra;//!
+
+ public: // Spectra Getter
+ map< vector<TString> , TH1*> GetSpectra();
+
+ ClassDef(TCharissaPhysics,1) // CharissaPhysics structure
+};
+
+namespace CHARISSA_LOCAL
+{
+
+ // tranform an integer to a string
+ string itoa(int value);
+ // DSSD
+ // X
+ double fSi_X_E(const TCharissaData* Data, const int i);
+ double fSi_X_T(const TCharissaData* Data, const int i);
+
+ // Y
+ double fSi_Y_E(const TCharissaData* Data, const int i);
+ double fSi_Y_T(const TCharissaData* Data, const int i);
+
+ // CsI
+ double fCsI_E(const TCharissaData* Data, const int i);
+ double fCsI_T(const TCharissaData* Data, const int i);
+
+}
+
+
+#endif
diff --git a/NPLib/Charissa/TCharissaSpectra.cxx b/NPLib/Charissa/TCharissaSpectra.cxx
index f0bf4a164a9e458215a8b21c92b8f97b1d396bc4..c09bfa8af976307451044cd989acf757d67b19b1 100644
--- a/NPLib/Charissa/TCharissaSpectra.cxx
+++ b/NPLib/Charissa/TCharissaSpectra.cxx
@@ -447,68 +447,59 @@ void TCharissaSpectra::FillPreTreatedSpectra(TCharissaData* PreTreatedData)
}
-/*
+
////////////////////////////////////////////////////////////////////////////////
void TCharissaSpectra::FillPhysicsSpectra(TCharissaPhysics* Physics)
{
+ cout << "TCharissaSpactra::FillPhysicsSpectra has to be implemented !" << endl;
+
TString name;
TString family= "CHARISSA/PHY";
// X-Y Impact Matrix
- for(unsigned int i = 0 ; i < Physics->Si_E.size(); i++){
+ for(unsigned int i = 0 ; i < Physics->Layer1_Si_E.size(); i++){
name = "CHA_IMPACT_MATRIX";
double x = Physics->GetPositionOfInteraction(i).x();
double y = Physics->GetPositionOfInteraction(i).y();
GetHisto(family,name)-> Fill(x,y);
- name = "CHA_THETA_E";
+ name = "L1_CHA_THETA_E";
double Theta = Physics->GetPositionOfInteraction(i).Angle(TVector3(0,0,1));
Theta = Theta/deg;
- GetHisto(family,name)-> Fill(Theta,Physics->Si_E[i]);
+ GetHisto(family,name)-> Fill(Theta,Physics->Layer1_Si_E[i]);
// STRX_E_CAL
- name = Form("CHA%d_XY_COR", Physics->TelescopeNumber[i]);
- GetHisto(family,name)-> Fill(Physics->Si_EX[i],Physics->Si_EY[i]);
-
-
- // Fill only for particle stopped in the first stage
- if(Physics->SiLi_E[i]<0 && Physics->CsI_E[i]<0){
+ name = Form("CHA%d_XY_COR", Physics->Layer1_TelescopeNumber[i]);
+ GetHisto(family,name)-> Fill(Physics->Layer1_Si_EX[i],Physics->Layer1_Si_EY[i]);
+
+
// E-TOF:
name = "CHA_E_TOF";
- GetHisto(family,name)->Fill(Physics->Si_E[i],Physics->Si_T[i]);
+ GetHisto(family,name)->Fill(Physics->Layer1_Si_E[i],Physics->Layer1_Si_T[i]);
- name = Form("CHA%d_E_TOF", Physics->TelescopeNumber[i]);
- GetHisto(family,name)->Fill(Physics->Si_E[i],Physics->Si_T[i]);
- }
+ name = Form("CHA%d_E_TOF", Physics->Layer1_TelescopeNumber[i]);
+ GetHisto(family,name)->Fill(Physics->Layer1_Si_E[i],Physics->Layer1_Si_T[i]);
- double Etot=0;
- if(Physics->SiLi_E[i]>0){
- name = "CHA_SILIE_E";
- Etot = Physics->SiLi_E[i];
- GetHisto(family,name)->Fill(Physics->SiLi_E[i],Physics->Si_E[i]);
-
- name = Form("CHA%d_SILIE_E", Physics->TelescopeNumber[i]);
- GetHisto(family,name)->Fill(Physics->SiLi_E[i],Physics->Si_E[i]);
- }
+ double Etot;
if(Physics->CsI_E[i]>0){
name = "CHA_CSIE_E";
- Etot += Physics->CsI_E[i];
- GetHisto(family,name)->Fill(Physics->CsI_E[i],Physics->Si_E[i]);
- name = Form("CHA%d_CSIE_E", Physics->TelescopeNumber[i]);
- GetHisto(family,name)->Fill(Physics->CsI_E[i],Physics->Si_E[i]);
+ Etot = Physics->Layer1_Si_E[i]+Physics->Layer2_Si_E[i]+Physics->CsI_E[i];
+ GetHisto(family,name)->Fill(Physics->CsI_E[i],Physics->Layer1_Si_E[i]+Physics->Layer2_Si_E[i]);
+ name = Form("CHA%d_CSIE_E", Physics->Layer1_TelescopeNumber[i]);
+ GetHisto(family,name)->Fill(Physics->CsI_E[i],Physics->Layer1_Si_E[i]+Physics->Layer2_Si_E[i]);
}
if(Etot>0){
name = "CHA_Etot_E";
- GetHisto(family,name)->Fill(Etot,Physics->Si_E[i]);
- name = Form("CHA%d_Etot_E", Physics->TelescopeNumber[i]);
- GetHisto(family,name)->Fill(Etot,Physics->Si_E[i]);
+ GetHisto(family,name)->Fill(Etot,Physics->Layer1_Si_E[i]);
+ name = Form("CHA%d_Etot_E", Physics->Layer1_TelescopeNumber[i]);
+ GetHisto(family,name)->Fill(Etot,Physics->Layer1_Si_E[i]);
}
}
-}*/
+}
diff --git a/NPLib/Charissa/TCharissaSpectra.h b/NPLib/Charissa/TCharissaSpectra.h
index 1002cbdcc36514878e4ccf82ff08f472ae0644c4..56bbd447754c61a06072f42897a9d5c9a094c271 100644
--- a/NPLib/Charissa/TCharissaSpectra.h
+++ b/NPLib/Charissa/TCharissaSpectra.h
@@ -1,5 +1,5 @@
-#ifndef TMUST2SPECTRA_H
-#define TMUST2SPECTRA_H
+#ifndef TCHARISSASPECTRA_H
+#define TCHARISSASPECTRA_H
/*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
@@ -34,11 +34,11 @@
// NPLib headers
#include "TCharissaData.h"
-//#include "TCharissaPhysics.h"
+#include "TCharissaPhysics.h"
using namespace std;
// ForwardDeclaration
-//class TCharissaPhysics;
+class TCharissaPhysics;
class TCharissaSpectra {
public:
@@ -62,7 +62,7 @@ class TCharissaSpectra {
// Filling methods
void FillRawSpectra(TCharissaData*);
void FillPreTreatedSpectra(TCharissaData*);
-// void FillPhysicsSpectra(TCharissaPhysics*);
+ void FillPhysicsSpectra(TCharissaPhysics*);
// Check the Spectra
void CheckSpectra();
@@ -72,7 +72,7 @@ class TCharissaSpectra {
TH1* GetHisto(TString family,TString name);
void WriteHisto(TString filename="VOID");
- private: // Information on MUST2
+ private: // Information on CHARISSA
unsigned int fNumberOfTelescope;
unsigned int fStripX;
unsigned int fStripY;
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index efa74ddfc3813f14b38c0e1ff8d31386b599a179..b26c96698e2d88912b9e3198bdd8e245fe7f5d12 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -65,6 +65,9 @@ Reaction::Reaction(){
// Need to be done before initializePrecomputeVariable
fKineLine3 = 0 ;
fKineLine4 = 0 ;
+ fLineBrho3 = 0 ;
+ fTheta3VsTheta4 = 0;
+ fAngleLine = 0;
//
fNuclei1 = new Beam();
@@ -127,6 +130,9 @@ Reaction::Reaction(string reaction){
fKineLine3 = 0 ;
fKineLine4 = 0 ;
+ fLineBrho3 = 0;
+ fTheta3VsTheta4 = 0;
+ fAngleLine = 0;
fNuclei1 = new Beam(A);
fNuclei2 = new Nucleus(b);
fNuclei3 = new Nucleus(c);
@@ -450,9 +456,6 @@ void Reaction::ReadConfigurationFile(string Path){
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::initializePrecomputeVariable(){
- // delete the previously calculated kinematical line:
- if(fKineLine3!=0) {delete fKineLine3 ; fKineLine3 = 0;}
- if(fKineLine4!=0) {delete fKineLine4 ; fKineLine4 = 0;}
m1 = fNuclei1->Mass();
m2 = fNuclei2->Mass();
@@ -463,13 +466,15 @@ void Reaction::initializePrecomputeVariable(){
fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
- ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
- ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
- ECM_3 = (s + m3*m3 - m4*m4)/(2*sqrt(s));
- ECM_4 = (s + m4*m4 - m3*m3)/(2*sqrt(s));
+ ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
+ ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
+ ECM_3 = (s + m3*m3 - m4*m4)/(2*sqrt(s));
+ ECM_4 = (s + m4*m4 - m3*m3)/(2*sqrt(s));
- pCM_3 = sqrt(ECM_3*ECM_3 - m3*m3);
- pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
+ pCM_1 = sqrt(ECM_1*ECM_1 - m1*m1);
+ pCM_2 = sqrt(ECM_2*ECM_2 - m2*m2);
+ pCM_3 = sqrt(ECM_3*ECM_3 - m3*m3);
+ pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
fImpulsionLab_1 = TVector3(0,0,sqrt(fBeamEnergy*fBeamEnergy + 2*fBeamEnergy*m1));
fImpulsionLab_2 = TVector3(0,0,0);
@@ -506,93 +511,103 @@ void Reaction::SetNuclei3(double EnergyLab, double ThetaLab){
////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetKinematicLine3(){
-
- if(fKineLine3==0){
- int size = 360;
- double x[size];
- double y[size];
- double theta3,E3,theta4,E4;
-
- for (int i = 0; i < size; ++i){
- SetThetaCM(((double)i)/2*deg);
- KineRelativistic(theta3, E3, theta4, E4);
- fNuclei3->SetKineticEnergy(E3);
-
- x[i] = theta3/deg;
- y[i] = E3;
- }
+TGraph* Reaction::GetKinematicLine3(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta3,E3,theta4,E4;
- fKineLine3 = new TGraph(size,x,y);
- // fKineLine3->SetTitle("Kinematic Line of particule 3");
- }
+ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei3->SetKineticEnergy(E3);
+
+ vx.push_back(theta3/deg);
+ vy.push_back(E3);
+ }
+ fKineLine3 = new TGraph(vx.size(),&vx[0],&vy[0]);
+
return(fKineLine3);
}
-
////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetKinematicLine4(){
- if(fKineLine4==0){
- int size = 360;
- double x[size];
- double y[size];
+TGraph* Reaction::GetKinematicLine4(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
double theta3,E3,theta4,E4;
- for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
- KineRelativistic(theta3, E3, theta4, E4);
- fNuclei4->SetKineticEnergy(E4);
-
- x[i] = theta4/deg;
- y[i] = E4;
- }
- fKineLine4= new TGraph(size,x,y);
- // fKineLine4->SetTitle("Kinematic Line of particule 4");
- }
+ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei4->SetKineticEnergy(E4);
+
+ vx.push_back(theta3/deg);
+ vy.push_back(E4);
+ }
+ fKineLine4= new TGraph(vx.size(),&vx[0],&vy[0]);
+
return(fKineLine4);
}
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetTheta3VsTheta4(double AngleStep_CM)
+{
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta3,E3,theta4,E4;
+
+ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+
+ vx.push_back(theta3/deg);
+ vy.push_back(theta4/deg);
+ }
+ fTheta3VsTheta4= new TGraph(vx.size(),&vx[0],&vy[0]);
+ return(fTheta3VsTheta4);
+}
+
////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetBrhoLine3(){
- int size = 360;
- double x[size];
- double y[size];
+TGraph* Reaction::GetBrhoLine3(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
double theta3,E3,theta4,E4;
double Brho;
- for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*deg);
KineRelativistic(theta3, E3, theta4, E4);
fNuclei3->SetKineticEnergy(E3);
Brho = fNuclei3->GetBrho();
- x[i] = theta3/deg;
- y[i] = Brho;
+ vx.push_back(theta3/deg);
+ vy.push_back(Brho);
}
- TGraph* LineBrho3= new TGraph(size,x,y);
- return(LineBrho3);
+ fLineBrho3= new TGraph(vx.size(),&vx[0],&vy[0]);
+ return(fLineBrho3);
}
////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetThetaLabVersusThetaCM(){
- int size = 360;
- double x[size];
- double y[size];
+TGraph* Reaction::GetThetaLabVersusThetaCM(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
double theta3,E3,theta4,E4;
- for (int i = 0; i < size; ++i){
- SetThetaCM(((double)i)/2*deg);
+ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*deg);
KineRelativistic(theta3, E3, theta4, E4);
- x[i] = fThetaCM/deg;
- y[i] = theta3/deg;
+ vx.push_back(fThetaCM/deg);
+ vy.push_back(theta3/deg);
}
- TGraph* AngleLine= new TGraph(size,x,y);
- return(AngleLine);
+ fAngleLine= new TGraph(vx.size(),&vx[0],&vy[0]);
+ return(fAngleLine);
}
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index f419e9595739f53fa83bcc77aa2905b8b1ace38a..9525e97a4a7412ddc554a9f2877c575bd9a0e2e6 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -74,8 +74,11 @@ namespace NPL{
bool fshoot4;
private: // use to display the kinematical line
- TGraph* fKineLine3 ;
- TGraph* fKineLine4 ;
+ TGraph* fKineLine3 ;
+ TGraph* fKineLine4 ;
+ TGraph* fTheta3VsTheta4;
+ TGraph* fLineBrho3;
+ TGraph* fAngleLine;
private:
Beam* fNuclei1; // Beam
Nucleus* fNuclei2; // Target
@@ -150,12 +153,14 @@ namespace NPL{
TVector3 fImpulsionCM_4;
// CM Energy composante & CM impulsion norme
- Double_t ECM_1;
- Double_t ECM_2;
- Double_t ECM_3;
- Double_t ECM_4;
- Double_t pCM_3;
- Double_t pCM_4;
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t ECM_3;
+ Double_t ECM_4;
+ Double_t pCM_1;
+ Double_t pCM_2;
+ Double_t pCM_3;
+ Double_t pCM_4;
// Mandelstam variable
Double_t s;
@@ -185,11 +190,20 @@ namespace NPL{
void SetNuclei3(double EnergyLab, double ThetaLab);
- TGraph* GetKinematicLine3();
- TGraph* GetKinematicLine4();
- TGraph* GetBrhoLine3();
- TGraph* GetThetaLabVersusThetaCM();
- void PrintKinematic();
+ TGraph* GetKinematicLine3(double AngleStep_CM=1);
+ TGraph* GetKinematicLine4(double AngleStep_CM=1);
+ TGraph* GetBrhoLine3(double AngleStep_CM=1);
+ TGraph* GetThetaLabVersusThetaCM(double AngleStep_CM=1);
+ TGraph* GetTheta3VsTheta4(double AngleStep_CM=1);
+ void PrintKinematic();
+ double GetP_CM_1() {return pCM_1;}
+ double GetP_CM_2() {return pCM_2;}
+ double GetP_CM_3() {return pCM_3;}
+ double GetP_CM_4() {return pCM_4;}
+ double GetE_CM_1() {return ECM_1;}
+ double GetE_CM_2() {return ECM_2;}
+ double GetE_CM_3() {return ECM_3;}
+ double GetE_CM_4() {return ECM_4;}
// Print private paremeter
void Print() const;
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 6b9518885a62dc4e5d672d286e594760bca4fea5..795f8badd4ad4a6a9038df828ce9201e2c308295 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -45,6 +45,7 @@
#include "Paris.h"
#include "TW1Physics.h"
#include "TS2Physics.h"
+#include "TCharissaPhysics.h"
#include "Shield.h"
#include "TSpegPhysics.h"
#include "TExlPhysics.h"
@@ -85,6 +86,7 @@ void DetectorManager::ReadConfigurationFile(string Path) {
Bool_t ScintillatorPlastic = false;
Bool_t IonisationChamber = false;
Bool_t Trifoil = false;
+ Bool_t Charissa = false;
Bool_t GeneralTarget = false;
Bool_t GPDTracker = false;
Bool_t HYD2Tracker = false;
@@ -211,6 +213,27 @@ void DetectorManager::ReadConfigurationFile(string Path) {
AddDetector("MUST2", myDetector);
#endif
}
+
+ ////////////////////////////////////////////
+ //////// Search for CHARISSA Array ////////
+ ////////////////////////////////////////////
+ else if (LineBuffer.compare(0, 13, "CharissaArray") == 0 && Charissa == false) {
+#ifdef INC_CHARISSA
+ Charissa = true;
+ cout << "//////// Charissa Array ////////" << endl << endl;
+
+ // Instantiate the new array as a VDetector Object
+ VDetector* myDetector = new TCharissaPhysics();
+
+ // Read Position of Telescope
+ ConfigFile.close();
+ myDetector->ReadConfiguration(Path);
+ ConfigFile.open(Path.c_str());
+
+ // Add array to the VDetector Vector
+ AddDetector("Charissa", myDetector);
+#endif
+ }
////////////////////////////////////////////
//////// Search for CATS Array ////////