diff --git a/NPLib/CATS/Makefile b/NPLib/CATS/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..7895a0598a3d8058c719436419c3f936d82cccf7
--- /dev/null
+++ b/NPLib/CATS/Makefile
@@ -0,0 +1,43 @@
+include ../Makefile.arch
+
+#------------------------------------------------------------------------------
+SHARELIB = libCATSData.so libCATSPhysics.so
+
+all: $(SHARELIB)
+#------------------------------------------------------------------------------
+############### Detector ##############
+
+## CATS ##
+libCATSData.so: TCATSData.o TCATSDataDict.o
+ $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+
+TCATSDataDict.cxx: TCATSData.h
+ rootcint -f $@ -c $^
+
+libCATSPhysics.so: TCATSPhysics.o TCATSPhysicsDict.o
+ $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+
+TCATSPhysicsDict.cxx: TCATSPhysics.h
+ rootcint -f $@ -c $^
+
+# dependances
+TCATSData.o: TCATSData.cxx TCATSData.h
+TCATSDataDict.o: TCATSDataDict.cxx TCATSDataDict.h
+TCATSPhysics.o: TCATSPhysics.cxx TCATSPhysics.h
+#######################################
+
+############# Clean and More ##########
+clean:
+ @rm -f core *~ *.o *Dict*
+
+distclean:
+ make clean; rm -f *.so
+
+.SUFFIXES: .$(SrcSuf)
+
+###
+
+.$(SrcSuf).$(ObjSuf):
+ $(CXX) $(CXXFLAGS) $(INCLUDE) -c $<
+
+
diff --git a/NPLib/CATS/TCATSData.cxx b/NPLib/CATS/TCATSData.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3f0366121d7f2d61040d4062410b6c562f39f0e4
--- /dev/null
+++ b/NPLib/CATS/TCATSData.cxx
@@ -0,0 +1,87 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: Sandra Giron contact address: giron@ipno.in2p3.fr *
+ * *
+ * Creation Date : febuary 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold CATS raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#include <iostream>
+
+#include "TCATSData.h"
+
+
+ClassImp(TCATSData)
+
+TCATSData::TCATSData()
+{
+ // Default constructor
+
+ // X
+ fCATS_DetX.clear();
+ fCATS_StripX.clear();
+ fCATS_ChargeX.clear();
+ // Y
+ fCATS_DetY.clear();
+ fCATS_StripY.clear();
+ fCATS_ChargeY.clear();
+ // (Qfil)
+ fCATS_DetQ.clear();
+ fCATS_Charge.clear();
+}
+
+
+
+TCATSData::~TCATSData()
+{
+}
+
+
+
+void TCATSData::Clear()
+{
+ // X
+ fCATS_DetX.clear();
+ fCATS_StripX.clear();
+ fCATS_ChargeX.clear();
+ // Y
+ fCATS_DetY.clear();
+ fCATS_StripY.clear();
+ fCATS_ChargeY.clear();
+ // (Qfil)
+ fCATS_DetQ.clear();
+ fCATS_Charge.clear();
+}
+
+
+
+void TCATSData::Dump()
+{
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+
+ // X
+ cout << "CATS_MultX = " << fCATS_DetX.size() << endl;
+ for (UShort_t i = 0; i < fCATS_DetX.size(); i++)
+ cout << "DetX: " << fCATS_DetX[i] << " StripX: " << fCATS_StripX[i] << " ChargeX: " << fCATS_ChargeX[i] << endl;
+ // Y
+ cout << "CATS_MultY = " << fCATS_DetY.size() << endl;
+ for (UShort_t i = 0; i < fCATS_DetY.size(); i++)
+ cout << "DetY: " << fCATS_DetY[i] << " StripY: " << fCATS_StripY[i] << " EnergyY: " << fCATS_ChargeY[i] << endl;
+ // (Qfil)
+ cout << "MM_MultQ = " << fCATS_DetQ.size() << endl;
+ for (UShort_t i = 0; i < fCATS_DetQ.size(); i++)
+ cout << "DetQ: " << fCATS_DetQ[i] << " Charge: " << fCATS_Charge[i] << endl;
+}
diff --git a/NPLib/CATS/TCATSData.h b/NPLib/CATS/TCATSData.h
new file mode 100644
index 0000000000000000000000000000000000000000..cbe1358c2d4898a1774870a95108797145b9b2dc
--- /dev/null
+++ b/NPLib/CATS/TCATSData.h
@@ -0,0 +1,86 @@
+#ifndef __CATSDATA__
+#define __CATSDATA__
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: Sandra Giron contact address: giron@ipno.in2p3.fr *
+ * *
+ * Creation Date : febuary 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold CATS raw data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+#include <vector>
+
+#include "TObject.h"
+
+using namespace std;
+
+
+class TCATSData : public TObject {
+ protected:
+ // X strips
+ vector<UShort_t> fCATS_DetX;
+ vector<UShort_t> fCATS_StripX;
+ vector<UShort_t> fCATS_ChargeX;
+
+ // Y strips
+ vector<UShort_t> fCATS_DetY;
+ vector<UShort_t> fCATS_StripY;
+ vector<UShort_t> fCATS_ChargeY;
+
+ //Q Fil
+ vector<UShort_t> fCATS_DetQ;
+ vector<UShort_t> fCATS_Charge;
+
+ public:
+ TCATSData();
+ virtual ~TCATSData();
+
+ void Clear();
+ void Dump();
+
+ ///////////////////// SETTERS ////////////////////////
+ // X
+ void SetCATSDetX(UShort_t DetX) {fCATS_DetX.push_back(DetX);}
+ void SetCATSStripX(UShort_t StripX) {fCATS_StripX.push_back(StripX);}
+ void SetCATSChargeX(UShort_t ChargeX) {fCATS_ChargeX.push_back(ChargeX);}
+ // Y
+ void SetCATSDetY(UShort_t DetY) {fCATS_DetY.push_back(DetY);}
+ void SetCATSStripY(UShort_t StripY) {fCATS_StripY.push_back(StripY);}
+ void SetCATSChargeY(UShort_t ChargeY) {fCATS_ChargeY.push_back(ChargeY);}
+
+ //Q fil
+ void SetCATSDetQ(UShort_t DetQ) {fCATS_DetQ.push_back(DetQ);}
+ void SetCATSCharge(UShort_t Charge) {fCATS_Charge.push_back(Charge);}
+
+ ///////////////////// GETTERS ////////////////////////
+ // X
+ UShort_t GetCATSMultX() {return fCATS_DetX.size();}
+ UShort_t GetCATSDetX(Int_t i) {return fCATS_DetX.at(i);}
+ UShort_t GetCATSStripX(Int_t i) {return fCATS_StripX.at(i);}
+ UShort_t GetCATSChargeX(Int_t i) {return fCATS_ChargeX.at(i);}
+ // Y
+ UShort_t GetCATSMultY() {return fCATS_DetY.size();}
+ UShort_t GetCATSDetY(Int_t i) {return fCATS_DetY.at(i);}
+ UShort_t GetCATSStripY(Int_t i) {return fCATS_StripY.at(i);}
+ UShort_t GetCATSChargeY(Int_t i) {return fCATS_ChargeY.at(i);}
+ //Q fil
+ UShort_t GetCATSMultQ() {return fCATS_DetQ.size();}
+ UShort_t GetCATSDetQ(Int_t i) {return fCATS_DetQ.at(i);}
+ UShort_t GetCATSCharge(Int_t i) {return fCATS_Charge.at(i);}
+
+ ClassDef(TCATSData,1) // CATSData structure
+};
+
+#endif
diff --git a/NPLib/CATS/TCATSPhysics.cxx b/NPLib/CATS/TCATSPhysics.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..90e8e20a2d50fede385436cfd213f925e73b50c7
--- /dev/null
+++ b/NPLib/CATS/TCATSPhysics.cxx
@@ -0,0 +1,1752 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: Sandra Giron contact address: giron@ipno.in2p3.fr *
+ * *
+ * Creation Date : febuary 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold CATS treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+
+#include "TCATSPhysics.h"
+using namespace LOCAL_CATS;
+
+// STL
+#include <cmath>
+#include <algorithm>
+#include <sstream>
+#include <fstream>
+#include <iostream>
+#include <stdlib.h>
+
+// NPL
+#include "RootInput.h"
+#include "RootOutput.h"
+
+// ROOT
+#include "TChain.h"
+#include "TRandom.h"
+
+ClassImp(TCATSPhysics)
+///////////////////////////////////////////////////////////////////////////
+TCATSPhysics::TCATSPhysics()
+{
+ EventData = new TCATSData() ;
+ EventPhysics = this ;
+ NumberOfCATS = 0 ;
+}
+
+
+
+TCATSPhysics::~TCATSPhysics()
+{
+}
+
+
+// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+void TCATSPhysics::ReadConfiguration(string Path)
+{
+ ifstream ConfigFile;
+ ConfigFile.open(Path.c_str());
+ string LineBuffer ;
+ string DataBuffer ;
+
+ double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
+ TVector3 A , B , C , D ;
+
+ bool check_A = false ;
+ bool check_B = false ;
+ bool check_C = false ;
+ bool check_D = false ;
+
+ bool ReadingStatus = false ;
+
+
+ while (!ConfigFile.eof())
+ {
+ getline(ConfigFile, LineBuffer);
+
+ //If line is a Start Up CATS bloc, Reading toggle to true
+ if (LineBuffer.compare(0, 12, "CATSDetector") == 0)
+ {
+ cout << "///" << endl ;
+ cout << "CATS Detector 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.compare(0, 12, "CATSDetector") == 0) {
+ cout << "WARNING: Another CATS is found before standard sequence of Token, Error may occured in CATS definition" << endl ;
+ ReadingStatus = false ;
+ }
+
+ // Corner Position method
+
+ else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az ;
+
+ A = TVector3(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "X28_Y1=") == 0) {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz ;
+
+ B = TVector3(Bx, By, Bz);
+ cout << "X28 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 7, "X1_Y28=") == 0) {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz ;
+
+ C = TVector3(Cx, Cy, Cz);
+ cout << "X1 Y28 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 8, "X28_Y28=") == 0) {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz ;
+
+ D = TVector3(Dx, Dy, Dz);
+ cout << "X28 Y28 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+
+ }
+
+ // End Corner Position Method
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if (check_A && check_B && check_C && check_D)
+ {
+ ReadingStatus = false;
+
+ ///Add The previously define telescope
+
+ AddCATS( A ,
+ B ,
+ C ,
+ D );
+
+ check_A = false;
+ check_B = false;
+ check_C = false;
+ check_D = false;
+ }
+ }
+
+ }
+
+ // Should be called in Analysis...
+ ReadPedestal("../../../offline/calibrations/CATS/CATScoeff/Piedestaux_368.txt");
+ // path given from e530 directory
+
+ cout << endl << "/////////////////////////////" << endl<<endl;
+
+}
+
+// Add Parameter to the CalibrationManger
+void TCATSPhysics::AddParameterToCalibrationManager()
+{
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ for(int i = 0 ; i < NumberOfCATS ; i++)
+ {
+
+ for( int j = 0 ; j < 28 ; j++)
+ {
+ Cal->AddParameter("CATS", "D"+itoa(i+1)+"_X"+itoa(j+1)+"_Q","CATS_D"+itoa(i+1)+"_X"+itoa(j+1)+"_Q") ;
+ Cal->AddParameter("CATS", "D"+itoa(i+1)+"_Y"+itoa(j+1)+"_Q","CATS_D"+itoa(i+1)+"_Y"+itoa(j+1)+"_Q") ;
+ }
+ }
+}
+
+// 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 TCATSPhysics::InitializeRootInput()
+{
+ TChain* inputChain = RootInput::getInstance()->GetChain() ;
+ inputChain->SetBranchStatus( "CATS" , true ) ;
+ inputChain->SetBranchStatus( "fCATS_*" , true ) ;
+ inputChain->SetBranchAddress( "CATS" , &EventData ) ;
+}
+
+// Create associated branches and associated private member DetectorPhysics address
+void TCATSPhysics::InitializeRootOutput()
+{
+ TTree* outputTree = RootOutput::getInstance()->GetTree() ;
+ outputTree->Branch( "CATS" , "TCATSPhysics" , &EventPhysics ) ;
+}
+
+
+void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28)
+{
+ NumberOfCATS++ ;
+
+ // Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are along X axis)
+ TVector3 U = C_X28_Y1 - C_X1_Y1 ;
+ U = U.Unit() ;
+
+ // Vector V on Telescope Face (parallele to X Strip)
+ TVector3 V = C_X1_Y28 - C_X1_Y1 ;
+ V = V.Unit() ;
+
+ // Position Vector of Strip Center
+ TVector3 StripCenter ;
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1 ;
+
+ // Geometry Parameter
+
+ double Face = 71.12 ; //mm
+ double NumberOfStrip = 28 ;
+ double StripPitch = Face / NumberOfStrip ; //mm
+
+ // Buffer object to fill Position Array
+ vector<double> lineX ;
+ vector<double> lineY ;
+ //vector<double> lineZ ;
+
+ vector< vector< double > > OneDetectorStripPositionX ;
+ vector< vector< double > > OneDetectorStripPositionY ;
+ double OneDetectorStripPositionZ ;
+
+ // Moving StripCenter to 1.1 corner (strip center!) :
+ Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2) ;
+
+ //cout << "Strip_1_1X = " << Strip_1_1.X() << endl;
+ //cout << "Strip_1_1Y = " << Strip_1_1.Y() << endl;
+
+
+ for( int i = 0 ; i < 28 ; i++ )
+ {
+ lineX.clear() ;
+ lineY.clear() ;
+ // lineZ.clear() ;
+
+ for( int j = 0 ; j < 28 ; j++ )
+ {
+ StripCenter = Strip_1_1 + StripPitch*( i*U + j*V ) ;
+ //StripCenter += -TargetPosition ;
+ lineX.push_back( StripCenter.x() ) ;
+ lineY.push_back( StripCenter.y() ) ;
+ // lineZ.push_back( StripCenter.z() ) ;
+ }
+
+ OneDetectorStripPositionX.push_back(lineX);
+ OneDetectorStripPositionY.push_back(lineY);
+ }
+
+ OneDetectorStripPositionZ = C_X1_Y1.Z();
+
+ StripPositionX.push_back(OneDetectorStripPositionX) ;
+ StripPositionY.push_back(OneDetectorStripPositionY) ;
+ StripPositionZ.push_back(OneDetectorStripPositionZ) ;
+
+}
+
+
+// This method is called at each event read from the Input Tree. Aim is to build treat Raw dat in order to extract physical parameter.
+void TCATSPhysics::BuildPhysicalEvent()
+{
+ // voir les commentaire fait la ou la methode existe deja
+}
+
+// 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 TCATSPhysics::BuildSimplePhysicalEvent()
+{
+
+ // cout << EventData-> GetCATSMultX() << endl;
+ // pourquoi pas une methode avec qui prend que le strip max par exemple...
+}
+*/
+
+
+////////////////////////////////////////// LE RESTE DE LA CLASSE! ////////////////////////////////////////
+
+void TCATSPhysics::Clear()
+{
+ DetNumberX.clear() ;
+ StripX.clear() ;
+ ChargeX.clear() ;
+ ChargeSumX.clear() ;
+ MultOverThreshX.clear() ;
+ StripMaxX.clear() ;
+ DetNumberY.clear() ;
+ StripY.clear() ;
+ ChargeY.clear() ;
+ // ChargeY_test.clear() ;
+ ChargeSumY.clear() ;
+ MultOverThreshY.clear() ;
+ StripMaxY.clear() ;
+ // StripMaxY_test.clear() ;
+ DetNumberX_Position.clear() ;
+ DetNumberY_Position.clear() ;
+ DetNumberZ_Position.clear() ;
+ PositionX.clear() ;
+ PositionY.clear() ;
+ PositionZ.clear() ;
+
+ ReconstructionMethodX.clear() ;
+ ReconstructionMethodY.clear() ;
+
+ //FailedReconstructionX.clear() ;
+ FailedReconstructionY.clear() ;
+
+ ff = 0;
+ HitX = 0;
+ HitY = 0;
+
+ NumberOfCATS = 0;
+}
+
+
+
+void TCATSPhysics::Dump()
+{
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+
+ cout << "Number Of CATS : " << NumberOfCATS << endl;
+
+ for(unsigned int i= 0; i < DetNumberX.size() ; i++)
+ {
+ cout << "DetNumberX : " << DetNumberX.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < StripX.size() ; i++)
+ {
+ cout << "StripX : " << StripX.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < StripMaxX.size() ; i++)
+ {
+ cout << "StripMaxX : " << StripMaxX.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < ChargeX.size() ; i++)
+ {
+ cout << "ChargeX : " << ChargeX.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < ChargeSumX.size() ; i++)
+ {
+ cout << "ChargeSumX : " << ChargeSumX.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < MultOverThreshX.size() ; i++)
+ {
+ cout << "MultOverThreshX : " << MultOverThreshX.at(i) << endl;
+ }
+
+
+
+ for(unsigned int i= 0; i < DetNumberY.size() ; i++)
+ {
+ cout << "DetNumberY : " << DetNumberY.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < StripY.size() ; i++)
+ {
+ cout << "StripY : " << StripY.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < StripMaxY.size() ; i++)
+ {
+ cout << "StripMaxY : " << StripMaxY.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < ChargeY.size() ; i++)
+ {
+ cout << "ChargeY : " << ChargeY.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < ChargeSumY.size() ; i++)
+ {
+ cout << "ChargeSumY : " << ChargeSumY.at(i) << endl;
+ }
+ for(unsigned int i= 0; i < MultOverThreshY.size() ; i++)
+ {
+ cout << "MultOverThreshY : " << MultOverThreshY.at(i) << endl;
+ }
+
+
+ for(unsigned int i = 0; i < DetNumberX_Position.size() ; i++)
+ {
+ cout << "DetNumberX_Position : " << DetNumberX_Position.at(i) << endl;
+ }
+ for(unsigned int i = 0; i < DetNumberY_Position.size() ; i++)
+ {
+ cout << "DetNumberY_Position : " << DetNumberY_Position.at(i) << endl;
+ }
+ for(unsigned int i = 0; i < DetNumberZ_Position.size() ; i++)
+ {
+ cout << "DetNumberZ_Position : " << DetNumberZ_Position.at(i) << endl;
+ }
+
+ for(unsigned int i = 0; i < PositionX.size() ; i++)
+ {
+ cout << "PositionX : " << PositionX.at(i) << endl;
+ }
+ for(unsigned int i = 0; i < PositionY.size() ; i++)
+ {
+ cout << "PositionY : " << PositionY.at(i) << endl;
+ }
+ for(unsigned int i = 0; i < PositionZ.size() ; i++)
+ {
+ cout << "PositionZ : " << PositionZ.at(i) << endl;
+ }
+
+ for(unsigned int i = 0; i < ReconstructionMethodX.size() ; i++)
+ {
+ cout << "ReconstructionMethodX : " << ReconstructionMethodX.at(i) << endl;
+ }
+
+ for(unsigned int i = 0; i < ReconstructionMethodY.size() ; i++)
+ {
+ cout << "ReconstructionMethodY : " << ReconstructionMethodY.at(i) << endl;
+ }
+ /*
+ for(unsigned int i = 0; i < FailedReconstructionX.size() ; i++)
+ {
+ cout << "FailedReconstructionX : " << FailedReconstructionX.at(i) << endl;
+ }
+ */
+ for(unsigned int i = 0; i < FailedReconstructionY.size() ; i++)
+ {
+ cout << "FailedReconstructionY : " << FailedReconstructionY.at(i) << endl;
+ }
+
+
+}
+
+
+
+
+void TCATSPhysics::ReadPedestal(string PedestalPath)
+{
+ int NumberOfStrips = 28 ;
+ vector<double> line ;
+ line.resize(NumberOfStrips, 0);
+
+ Pedestal_X.resize(NumberOfCATS, line);
+ Pedestal_Y.resize(NumberOfCATS, line);
+
+ Threshold_X.resize(NumberOfCATS, line);
+ Threshold_Y.resize(NumberOfCATS, line);
+
+ string DataBuffer, XY;
+ int StripNumber = 0;
+
+ ifstream Pedestal_File;
+ Pedestal_File.open(PedestalPath.c_str());
+
+ if( Pedestal_File.is_open() )
+ cout << " Calibration File " << PedestalPath << " loading " << endl;
+ else {
+ cout << " Error, no calibration file" << PedestalPath << " found" << endl;
+ return;
+ }
+
+ while( !Pedestal_File.eof() ) {
+
+ Pedestal_File >> DataBuffer ;
+
+ for(int i = 0 ; i < NumberOfCATS ; i++)
+ {
+ std::ostringstream DetectorNumber ;
+ DetectorNumber << i+1;
+
+ if(DataBuffer == "CATS"+DetectorNumber.str()+"_X")
+ {
+ for( int k = 0 ; k < NumberOfStrips ; k++ )
+ {
+ Pedestal_File >> StripNumber;
+ Pedestal_File >> DataBuffer; Pedestal_X[i][k] = atof(DataBuffer.c_str()) ;
+
+ // definition du seuil : piedestal + 3 * sigma(piedestal)
+ Pedestal_File >> DataBuffer; Threshold_X[i][k] = Pedestal_X[i][k]+3*atof(DataBuffer.c_str()) ;
+
+ }
+ }
+ else if(DataBuffer == "CATS"+DetectorNumber.str()+"_Y")
+ {
+ for( int k = 0 ; k < NumberOfStrips ; k++ )
+ {
+ Pedestal_File >> StripNumber;
+ Pedestal_File >> DataBuffer; Pedestal_Y[i][k] = atof(DataBuffer.c_str()) ;
+
+ // definition du seuil : piedestal + 3 * sigma(piedestal)
+ Pedestal_File >> DataBuffer; Threshold_Y[i][k] = Pedestal_Y[i][k]+3*atof(DataBuffer.c_str()) ; // definition du seuil
+
+ }
+ }
+ }
+ }
+
+}
+
+
+void TCATSPhysics::BuildSimplePhysicalEvent()
+{
+ /*
+ int HitX = 0 ;
+ int HitY = 0 ;
+ */
+
+ //cout << "CATS Build... " << endl;
+
+ gRandom->SetSeed(0);
+
+ //EventData->Dump();
+ double CalculatedStripX, CalculatedStripY ;
+ double posX =0 , posY = 0;
+
+
+ // How many CATS?
+ int NumberOfCATSHit = 0 ;
+ int DetectorID = -1;
+
+ for( unsigned short i = 0 ; i < EventData->GetCATSMultX() ; i++ )
+ {
+ // if( NumberOfCATSHit < EventData->GetCATSDetX(i) ) NumberOfCATSHit = EventData->GetCATSDetX(i) ; //determination of the number of CATS detectors
+
+ if( EventData->GetCATSDetX(i) != DetectorID) {
+ NumberOfCATSHit++;
+ }
+
+ if(NumberOfCATSHit == 2) break;
+
+ }
+
+
+ // cout << "Nombre de CATS: " << NumberOfCATSHit << endl;
+
+
+ // INITIALISATION OF VECTORS
+ for(int ff1 = 0 ; ff1 < NumberOfCATSHit ; ff1++ )
+ {
+ MultOverThreshX.push_back(-1);
+ StripMaxX.push_back(-1);
+ ChargeSumX.push_back(-1);
+ ReconstructionMethodX.push_back(NO);
+
+ MultOverThreshY.push_back(-1);
+ StripMaxY.push_back(-1);
+ // StripMaxY_test.push_back(-1);
+ ChargeSumY.push_back(-1);
+ ReconstructionMethodY.push_back(NO);
+
+ // FailedReconstructionX.push_back(NO);
+ FailedReconstructionY.push_back(NO);
+
+ }
+ // cout << "Init OK" << endl;
+
+
+
+ for(int gg = 0 ; gg < NumberOfCATSHit ; gg++ )
+ {
+ //int ff = NumberOfCATSHit - gg -1 ;
+ ff = gg ;
+
+ CalculatedStripX = AnalyseX(ff, NumberOfCATSHit); // cout << "Analyse X = " << CalculatedStripX << endl;
+ CalculatedStripY = AnalyseY(ff, NumberOfCATSHit); // cout << "Analyse Y = " << CalculatedStripY << endl;
+
+ posX = CalculatePositionX(CalculatedStripX, cor); // cout << "Position X = " << posX << endl;
+ posY = CalculatePositionY(CalculatedStripY, cor); // cout << "Position Y = " << posY << endl;
+
+ DetNumberX_Position.push_back(ff+1);
+ DetNumberY_Position.push_back(ff+1);
+ DetNumberZ_Position.push_back(ff+1);
+
+ PositionX.push_back(posX) ;
+ PositionY.push_back(posY) ;
+ PositionZ.push_back(StripPositionZ[ff]) ;
+ }
+
+
+ if(NumberOfCATSHit>1)
+ {
+ /*
+ cout << "PositionX[0]= " <<PositionX[0] << " PositionX[1] = " << PositionX[1] << " PositionX[0]-PositionX[1]= " << PositionX[0]-PositionX[1] << endl;
+ cout << "PositionY[0]= " <<PositionY[0] << " PositionY[1] = " << PositionY[1] << " PositionY[0]-PositionY[1]= " << PositionY[0]-PositionY[1] << endl;
+ cout << "PositionZ[0] = "<< PositionZ[0] << " PositionZ[1] = "<< PositionZ[1] << endl;
+ */
+ if(PositionX[0] > -35 && PositionX[0] < 35 && PositionY[0] > -35 && PositionY[0] < 35 && PositionX[1] > -35 && PositionX[1] < 35 && PositionY[1] > -35 && PositionY[1] < 35 )
+ // if(PositionX[0] != -1000 && PositionY[0] != -1000 && PositionX[1] != -1000 && PositionY[1] != -1000)
+ {
+ BeamDirection = TVector3 (PositionX[1]-PositionX[0] ,
+ PositionY[1]-PositionY[0] ,
+ PositionZ[1]-PositionZ[0] );
+ // BeamDirection.Unit(); Adrien
+
+ double l = sqrt((PositionZ[0]-PositionZ[1])*(PositionZ[0]-PositionZ[1]));
+ double L = - PositionZ[1];
+
+ double t = (l+L) / l;
+
+ PositionOnTargetX = PositionX[0] + (PositionX[1] - PositionX[0]) * t ;
+ PositionOnTargetY = PositionY[0] + (PositionY[1] - PositionY[0]) * t ;
+
+ //Reconstruction in S1 plane
+ /*
+ PositionOnTargetX = PositionX[0] + (PositionX[1] - PositionX[0]) * (l+403.8)/l ;
+ PositionOnTargetY = PositionY[0] + (PositionY[1] - PositionY[0]) * (l+403.8)/l ;
+ */
+ //cout << PositionOnTargetX << " " << PositionOnTargetY << endl;
+
+ BeamDirection.Unit();
+ }
+
+ else
+ {
+ //cout << "One of the CATS position was not reconstructed ! Impossible to reconstruct position on target ..." << endl;
+
+ //cout << "PositionX[0] = " << PositionX[0] << " PositionY[0] = " << PositionY[0] << " PositionX[1] = " << PositionX[1] << " PositionY[1] = " << PositionY[1] << endl;
+
+
+ BeamDirection = TVector3 ( 1 ,
+ 0 ,
+ 0 );
+
+ PositionOnTargetX = 3000 ;
+ PositionOnTargetY = 3000 ;
+ }
+ }
+
+ else if(NumberOfCATSHit ==1)
+ {
+ cout << NumberOfCATSHit << endl;
+
+ BeamDirection = TVector3 ( 1 ,
+ 0 ,
+ 0 );
+ PositionOnTargetX = 5000 ;
+ PositionOnTargetY = 5000 ;
+ }
+
+ // cout << "CATS fine ! " << endl;
+
+ return;
+}
+
+
+double TCATSPhysics::AnalyseX(int ff,
+ int NumberOfCATSHit)
+{
+ // cout << "AnalyseX proccessing ..." << endl;
+
+ // double Chargex[28];
+ for(unsigned short z=0; z<28; z++) {
+ Chargex[z] = -1;
+ }
+
+ // int HitX = 0;
+ double ChargeSum_X = 0 ;
+ double ChargeX_Buffer = 0;
+ double CalculatedStripX=0;
+
+ for(UShort_t i =0; i<EventData->GetCATSMultX(); i++)
+ {
+ // cout << EventData->GetCATSDetX(i)<< endl;
+ //cout << EventData->GetCATSStripX(i)<< endl;
+ // cout << "ff+1 = " << ff+1 << endl;
+
+ if( EventData->GetCATSDetX(i) == ff+1 )
+ {
+ int NX = EventData->GetCATSDetX(i);
+ int StrX = EventData->GetCATSStripX(i) ; // cout << NX << " " << StrX << endl ;
+
+ if(NX > 2 || StrX > 28) cout << NX << " " << StrX << endl ;
+
+ double Q = EventData->GetCATSChargeX(i) + gRandom->Rndm() - Pedestal_X[NX-1][StrX-1] ;
+
+ ChargeX_Buffer = CalibrationManager::getInstance()->ApplyCalibration("CATS/D"+itoa(NX)+"_X"+itoa(StrX)+"_Q",Q);
+
+
+ if(EventData->GetCATSChargeX(i) > Threshold_X[NX-1][StrX-1] && NX <= NumberOfCATSHit && StrX < 28)
+ {
+ // cout << Threshold_X[NX-1][StrX-1] << endl;
+
+ // KNOWN INVERSION
+ if (StrX == 15 && NX == 1) StrX = 16 ;
+ else if (StrX == 16 && NX == 1) StrX = 15 ;
+
+ MultOverThreshX[ff]++;
+ ChargeSum_X += ChargeX_Buffer;
+ //ChargeSum += ChargeX_Buffer;
+ ChargeX.push_back( ChargeX_Buffer ); // cout << "ChargeX_Buffer = " << ChargeX_Buffer << endl;
+ Chargex[StrX-1] = ChargeX_Buffer ; //cout <<" Chargex[" << StrX-1 << "] " << Chargex[StrX-1] << endl;
+ StripX.push_back(StrX);
+ DetNumberX.push_back(NX) ;
+ HitX++;
+
+ if ( ChargeX[HitX-1] > Chargex[ StripMaxX[ff] -1 ] ) StripMaxX[ff] = StrX ;
+ //cout << "X " << ff+1 << " " << StrX << " " << StripMaxX[ff] << " " << HitX-1 << " " << ChargeX[HitX-1] << " " << StripMaxX[ff]-1 << " " << Chargex[StripMaxX[ff]-1] << endl;
+ // cout << "X " << ff+1 << " " << StrX << " " << StripMaxX[ff] << " " << HitX-1 << " " << ChargeX_Buffer << " " << StripMaxX[ff]-1 << " " << Chargex[StripMaxX[ff]-1] << endl;
+ }
+
+ }
+
+ ChargeSumX[ff] = ChargeSum_X;
+ }
+
+ //cout << "Charge x = " << Chargex << " StripMaxX = " << StripMaxX[ff] << endl;
+
+ ReconstructionMethodX[ff] = ChooseReconstruction("X");
+ //if (ff==1) cout << "StripMaxX = " << StripMaxX[ff] << endl;
+ //cout << ReconstructionMethodX[ff] << endl;
+
+ if(ReconstructionMethodX[ff] == SECHS)CalculatedStripX = HyperbolicSequentMethodX();
+ if(ReconstructionMethodX[ff] == GAUSS)CalculatedStripX = GaussianMethodX();
+ if(ReconstructionMethodX[ff] == BAR3) CalculatedStripX = Barycentric3MethodX();
+ if(ReconstructionMethodX[ff] == BAR4) CalculatedStripX = Barycentric4MethodX();
+ if(ReconstructionMethodX[ff] == BAR5) CalculatedStripX = Barycentric5MethodX(); // Chargex, StripMaxX[ff] );
+
+ if(CalculatedStripX < 35 && CalculatedStripX > -35) { } //FailedReconstructionX[ff] = NO ; }
+ // else { FailedReconstructionX[ff] = ReconstructionMethodX[ff] ; } // cout << CalculatedStripX << endl;}
+
+ //cout << "in AnalyseX : " << CalculatedStripX << endl;
+
+ // else cout << "Error in the choice of the method!" << endl;
+
+ // cout << "AnalyseX done!" << endl ;
+ return(CalculatedStripX);
+}
+
+double TCATSPhysics::AnalyseY(int ff,
+ int NumberOfCATSHit)
+{
+
+
+ // cout << "AnalyseY proccessing ..." << endl;
+
+ // double Chargey[28];
+ for(unsigned short z=0; z<28; z++) {
+ Chargey[z] = -1;
+ // Chargey_test[z] = -1;
+ }
+
+ // int HitY = 0;
+
+ double ChargeSum_Y = 0 ;
+ double ChargeY_Buffer = 0;
+
+ double CalculatedStripY=0;
+
+ for(UShort_t i =0; i<EventData->GetCATSMultY(); i++)
+ {
+ if( EventData->GetCATSDetY(i) == ff+1 )
+ {
+ int NY = EventData -> GetCATSDetY(i);
+ int StrY = EventData -> GetCATSStripY(i) ;// cout << NY << endl ; //" " << StrY << endl ;
+
+ if(NY > 2 || StrY > 32) cout << NY << " " << StrY << endl ;
+
+ if(StrY < 28)
+ {
+ double Q = EventData->GetCATSChargeY(i) + gRandom->Rndm() - Pedestal_Y[NY-1][StrY-1];
+
+ ChargeY_Buffer = CalibrationManager::getInstance()->ApplyCalibration("CATS/D"+itoa(NY)+"_Y"+itoa(StrY)+"_Q",Q);
+ }
+
+ else {ChargeY_Buffer = 0 ;}
+
+ if(EventData->GetCATSChargeY(i) > Threshold_Y[NY-1][StrY-1] && NY <= NumberOfCATSHit && StrY < 28)
+ {
+ // cout << Threshold_Y[NY-1][StrY-1] << endl;
+
+ // KNOWN INVERSION
+ if (StrY == 15 && NY == 1) StrY = 16 ;
+ else if (StrY == 16 && NY == 1) StrY = 15 ;
+ //else if (StrY == 4 && NY == 2) StrY = 14 ;
+ //else if (StrY == 14 && NY == 2) StrY = 4 ;
+
+ MultOverThreshY[ff]++;
+ //ChargeSum += ChargeY_Buffer;
+ ChargeSum_Y += ChargeY_Buffer;
+
+ // Normal treatment
+ /*
+ ChargeY.push_back( ChargeY_Buffer ); // cout << "ChargeY_Buffer = " << ChargeY_Buffer << endl;
+ Chargey[StrY-1] = ChargeY_Buffer ; // cout <<" Chargey[" << StrY-1 << "] " << Chargey[StrY-1] << endl;
+ */
+
+ //Specific treatment for e530 experiment ///////////////////////////////////////////////////////////////////
+ // pist15 of cats2 not working...
+
+ if(ff ==1 && StrY ==15)
+ {
+ ChargeY.push_back( -1 );
+ Chargey[StrY-1] = -1 ;
+ }
+
+
+ else{
+ ChargeY.push_back( ChargeY_Buffer ); // cout << "ChargeY_Buffer = " << ChargeY_Buffer << endl;
+ Chargey[StrY-1] = ChargeY_Buffer ; // cout <<" Chargey[" << StrY-1 << "] " << Chargey[StrY-1] << endl;
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ StripY.push_back(StrY);
+ DetNumberY.push_back(NY) ;
+ HitY++;
+
+ if (ChargeY[HitY-1] > Chargey[ StripMaxY[ff] -1 ]) StripMaxY[ff] = StrY; // stripmax si pas de piste supprimee
+
+ //cout << "Y " << ff+1 << " " << StrY << " " << StripMaxY[ff] << " " << HitY-1 << " " << ChargeY[HitY-1] << " " << StripMaxY[ff]-1 << " " << Chargey[StripMaxY[ff]-1] << endl;
+ //cout << "Y " << ff+1 << " " << StrY << " " << StripMaxY[ff] << " " << HitY-1 << " " << ChargeY_Buffer << " " << StripMaxY[ff]-1 << " " << Chargey[StripMaxY[ff]-1] << endl;
+ }
+ }
+ //ChargeSumY[ff] = ChargeSum;
+ ChargeSumY[ff] = ChargeSum_Y;
+ }
+
+ // cout << "Charge y = " << Chargey << " StripMaxY = " << StripMaxY[ff] << endl;
+
+ ReconstructionMethodY[ff] = ChooseReconstruction("Y");
+ //if (ff==1) cout << "StripMaxY = " << StripMaxY[ff] << endl;
+ // cout << ReconstructionMethodY[ff] << endl;
+
+ if(ReconstructionMethodY[ff] == SECHS)CalculatedStripY = HyperbolicSequentMethodY();
+ if(ReconstructionMethodY[ff] == GAUSS)CalculatedStripY = GaussianMethodY() ;
+ if(ReconstructionMethodY[ff] == BAR3) CalculatedStripY = Barycentric3MethodY();
+ if(ReconstructionMethodY[ff] == BAR4) CalculatedStripY = Barycentric4MethodY();
+ if(ReconstructionMethodY[ff] == BAR5) CalculatedStripY = Barycentric5MethodY(); // Chargey, StripMaxY[ff] );
+
+ if(CalculatedStripY < 35 && CalculatedStripY > -35) { FailedReconstructionY[ff] = NO ; }
+ else FailedReconstructionY[ff] = ReconstructionMethodY[ff] ;
+
+ // else cout << "Error in the choice of the method!" << endl;
+
+ // cout << "AnalyseY done!" << endl;
+ return(CalculatedStripY);
+}
+
+reconstruction TCATSPhysics::ChooseReconstruction(TString type)
+{
+ reconstruction method = NO;
+ if(ff ==0) { method = SECHS; } //cout << "sechs" << endl; }
+
+ else {
+ if(type == "Y") {
+ method = GAUSS;
+ }
+ else if(type == "X") {
+ method = GAUSS;
+ }
+ }
+ return(method);
+}
+
+double TCATSPhysics::CalculatePositionX( double CalculatedStripX,
+ correction method)
+{
+ double positionX=-1000;
+ int IStripX = 0;
+
+ if(ReconstructionMethodX[ff] == GAUSS) positionX = CalculatedStripX; // already in mm -> see gaussian method
+
+ else
+ {
+ // cout << "CalculatedStripX = " << CalculatedStripX << endl;
+ // Integer part
+ IStripX = (int) CalculatedStripX ;
+
+ // Decimal Part
+ double DStripX = CalculatedStripX-IStripX ; // cout << "DStripX = " << DStripX << endl;
+
+ if( DStripX > 0.5) {IStripX++; DStripX = DStripX-1 ;} else {DStripX = DStripX;}
+
+ // Calculate Geometrical Position
+ if( IStripX > 0 && IStripX < 29 ) {
+ //if( IStripX>0 && IStripY>0 && StripMaxX[ff] > 0 && StripMaxY[ff] > 0 && StripMaxX[ff] < 29 && StripMaxY[ff] < 29 )
+
+ // positionX = (DStripX)*2.54 + StripPositionX[ff][IStripX-1][0] ; // conversion en mm initiale
+
+ if(ff==0) //CATS1
+ {
+ // Warning : DStrip sign has to be determined carefully depending on the X direction
+
+ positionX = (DStripX)*2.54 + StripPositionX[ff][IStripX-1][0] ; //init avec le moins
+ // positionX = 2.54 * (15-CalculatedStripX) - 1.27 ;
+ //cout << "ecartX1 = " << positionX - (-(DStripX)*2.54 + StripPositionX[ff][IStripX-1][0]) << endl;
+
+ if(method == NOcor) positionX = positionX;
+ else if(method == cor){
+ if(ReconstructionMethodX[ff] == BAR3) positionX = CorrectedPositionX3(positionX, 0.6);
+ if(ReconstructionMethodX[ff] == BAR4) positionX = CorrectedPositionX4(positionX, 0.77);
+ }
+
+ }
+
+ else if(ff==1) //CATS2
+ {
+ // Warning : DStrip sign has to be determined carefully depending on the X direction
+
+ positionX = -(DStripX)*2.54 + StripPositionX[ff][IStripX-1][0] ;
+ // positionX = 2.54 * (CalculatedStripX-15) + 1.27 ;
+ //cout << "ecartX2 = " << positionX - ((DStripX)*2.54 + StripPositionX[ff][IStripX-1][0]) << endl;
+
+ if(method == NOcor) positionX = positionX;
+ else if(method == cor){
+ if(ReconstructionMethodX[ff] == BAR3) positionX = CorrectedPositionX3(positionX, 0.53);
+ if(ReconstructionMethodX[ff] == BAR4) positionX = CorrectedPositionX4(positionX, 0.67);
+ }
+ }
+ else cout << "only 2CATS!! ff = " << ff << endl;
+ }
+
+ else { positionX = -1000; } // cout << CalculatedStripX << " " << IStripX << " " << ff << endl; }
+
+ // cout << "positionX " << positionX << " ff " << ff << " IStripX " << IStripX <<endl;
+ }
+
+
+ return(positionX);
+}
+
+
+double TCATSPhysics::CalculatePositionY( //int StripMaxY,
+ double CalculatedStripY,
+ // int ff,
+ correction method)
+{
+ double positionY = -1000;
+
+ if(ReconstructionMethodY[ff] == GAUSS) positionY = CalculatedStripY; // already in mm -> see gaussian method
+
+ else
+ {
+
+ // Integer part
+ int IStripY = (int) CalculatedStripY ;
+
+ // Decimal Part
+ double DStripY = CalculatedStripY-IStripY ;
+
+ if( DStripY > 0.5) {IStripY++; DStripY = DStripY-1 ; } else {DStripY = DStripY; }
+
+ // Calculate Geometrical Position
+ if(IStripY > 0 && IStripY < 29 )
+ // if(IStripX>0 && IStripY>0 && StripMaxX[ff] > 0 && StripMaxY[ff] > 0 && StripMaxX[ff] < 29 && StripMaxY[ff] < 29 )
+ {
+ positionY = (DStripY)*2.54 + StripPositionY[ff][0][IStripY-1] ; // conversion en mm initiale
+
+ /* if(ff== 1) {
+ cout << CalculatedStripY << endl;
+ cout << positionY << endl;
+ cout << StripPositionY[ff][0][IStripY-1] << endl;
+ }*/
+ //positionY = 2.54 * ( CalculatedStripY-14) - 1.27;
+ // cout << "ecartY" << positionY - ((DStripY)*2.54 + StripPositionY[ff][0][IStripY-1]) << endl;
+
+ if(ff ==0){
+ if(method == NOcor) positionY = positionY;
+ else if(method == cor) {
+ if(ReconstructionMethodY[ff] == BAR3) positionY = CorrectedPositionY3(positionY, 0.6);
+ if(ReconstructionMethodY[ff] == BAR4) positionY = CorrectedPositionY4(positionY, 0.75);
+ }
+ }
+
+ else if(ff ==1){
+ if(method == NOcor) positionY = positionY;
+ else if(method == cor){
+ if(ReconstructionMethodY[ff] == BAR3) positionY = CorrectedPositionY3(positionY, 0.57);
+ if(ReconstructionMethodY[ff] == BAR4) positionY = CorrectedPositionY4(positionY, 0.7);
+ }
+ }
+
+ else cout << "only 2CATS!! ff = " << ff << endl;
+ }
+
+ else { positionY = -1000 ; } // cout << CalculatedStripY << " " << IStripY << " " << ff << endl; }
+ // cout << IStripX << " " << IStripY << endl;
+
+ }
+
+ return(positionY);
+
+}
+
+double TCATSPhysics:: HyperbolicSequentMethodX()
+{
+ double sechs = -1000 ;
+
+ if(StripMaxX[ff] > 2 && StripMaxX[ff]<27)
+ {
+ /*
+ cout << "Chargex[" << StripMaxX[ff]-1 << "] = " << Chargex[StripMaxX[ff]-1-1] << endl;
+ cout << "Chargex[" << StripMaxX[ff] << "] = " << Chargex[StripMaxX[ff]-1] << endl;
+ cout << "Chargex[" << StripMaxX[ff]+1 << "] = " << Chargex[StripMaxX[ff]-1+1] << endl;
+ */
+
+ double vs1 = sqrt( Chargex[StripMaxX[ff]-1]/Chargex[StripMaxX[ff]-1+1] ) ;
+ double vs2 = sqrt( Chargex[StripMaxX[ff]-1]/Chargex[StripMaxX[ff]-1-1] ) ;
+ double vs3 = 0.5*( vs1 + vs2 ) ;
+ double vs4 = log( vs3 + sqrt(vs3*vs3-1.0) ) ;
+ double vs5 = (vs1 - vs2)/(2.0*sinh(vs4)) ;
+
+ if(vs5<0) vs5=-vs5 ;
+
+ double vs6 = 0.5*log( (1.0+vs5)/(1.0-vs5) ) ;
+
+ if ( Chargex[StripMaxX[ff]-1+1]>Chargex[StripMaxX[ff]-1-1] )
+ { sechs = StripMaxX[ff] + vs6/vs4 ;}
+
+ else
+ { sechs = StripMaxX[ff] - vs6/vs4 ;}
+
+ // cout << "vs1 = " << vs1 << " vs2 = " << vs2 << " vs3 = " << vs3 << " vs4 = " << vs4 << " vs5 = " << vs5 << " vs6 = " << vs6 << endl;
+
+ }
+
+ else {
+ sechs = -1000; //cout << " StripMaxX[ff] = " << StripMaxX[ff] << " out of range for SECHS !" << endl;
+ }
+
+
+ // cout << "sechs = " << sechs << endl;
+
+ return sechs ;
+}
+
+double TCATSPhysics:: HyperbolicSequentMethodY()
+{
+ double sechs = -1000 ;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff]<27)
+ {
+ /*
+ cout << "Chargey[" << StripMaxY[ff]-1 << "] = " << Chargey[StripMaxY[ff]-1-1] << endl;
+ cout << "Chargey[" << StripMaxY[ff] << "] = " << Chargey[StripMaxY[ff]-1] << endl;
+ cout << "Chargey[" << StripMaxY[ff]+1 << "] = " << Chargey[StripMaxY[ff]-1+1] << endl;
+ */
+
+ double vs1 = sqrt( Chargey[StripMaxY[ff]-1]/Chargey[StripMaxY[ff]-1+1] ) ;
+ double vs2 = sqrt( Chargey[StripMaxY[ff]-1]/Chargey[StripMaxY[ff]-1-1] ) ;
+ double vs3 = 0.5*( vs1 + vs2 ) ;
+ double vs4 = log( vs3 + sqrt(vs3*vs3-1.0) ) ;
+ double vs5 = (vs1 - vs2)/(2.0*sinh(vs4)) ;
+
+ if(vs5<0) vs5=-vs5 ;
+
+ double vs6 = 0.5*log( (1.0+vs5)/(1.0-vs5) ) ;
+
+ if ( Chargey[StripMaxY[ff]-1+1]>Chargey[StripMaxY[ff]-1-1] )
+ { sechs = StripMaxY[ff] + vs6/vs4 ;}
+
+ else
+ { sechs = StripMaxY[ff] - vs6/vs4 ;}
+
+ // cout << "vs1 = " << vs1 << " vs2 = " << vs2 << " vs3 = " << vs3 << " vs4 = " << vs4 << " vs5 = " << vs5 << " vs6 = " << vs6 << endl;
+
+ }
+
+ else {
+ sechs = -1000; //cout << " StripMaxY[ff] = " << StripMaxY[ff] << " out of range for SECHS !" << endl;
+ }
+
+
+ // cout << "sechs = " << sechs << endl;
+
+ return sechs ;
+}
+
+
+double TCATSPhysics:: GaussianMethodX()
+//int StripMax)
+{
+ int StripMax_ = StripMaxX[ff]- 1;
+ double gauss = -1000;
+ double Q[3];
+ double StripPos[3];
+
+
+ for(int j = 0; j<3 ; j++)
+ {
+ Q[j] = 0;
+ StripPos[j] = 0;
+ }
+
+ // cout << "StripMaxX[ff]= " << StripMaxX[ff]<<endl;
+
+ if(StripMaxX[ff]> 3 && StripMaxX[ff]< 26)
+ {
+ Q[0] = Chargex[StripMax_] ;
+ StripPos[0] = StripPositionX[ff][StripMax_][0];
+
+ // cout << "Q[0] = " << Q[0] << " StripMaxX[ff]= " << StripMaxX[ff]<< " StripPos[0] =" << StripPos[0] << endl;
+
+ if(Chargex[StripMax_-1] !=-1) {
+ //cout << "no pb" << endl;
+ Q[1] = Chargex[StripMax_-1];
+ StripPos[1] = StripPositionX[ff][StripMax_-1][0];
+
+ // cout << "Q[1] = " << Q[1] << " Strip1 = " << StripMax-1 << " StripPos[1] =" << StripPos[1] << endl;
+ }
+ else {
+ if(Chargex[StripMax_- 2] !=-1) {
+ //cout << "coconutsX-1! " << endl;
+ Q[1] = Chargex[StripMax_-2];
+ StripPos[1] = StripPositionX[ff][StripMax_ - 2][0];
+ // cout << "Q[0] = " << Q[0] << " StripMaxX[ff]= " << StripMaxX[ff]<< " StripPos[0] =" << StripPos[0] << endl;
+ // cout << "Q[1] = " << Q[1] << " Strip1 = " << StripMaxX[ff]-2 << " StripPos[1] =" << StripPos[1] << endl;
+ }
+ else {
+ if(Chargex[StripMax_- 3] !=-1) {
+ //cout << "coconutsX-2! " << endl;
+ Q[1] = Chargex[StripMax_- 3];
+ StripPos[1] = StripPositionX[ff][StripMax_ - 3][0];
+ // cout << "Q[0] = " << Q[0] << " StripMaxX[ff]= " << StripMaxX[ff]<< " StripPos[0] =" << StripPos[0] << endl;
+ // cout << "Q[1] = " << Q[1] << " Strip1 = " << StripMaxX[ff]-3 << " StripPos[1] =" << StripPos[1] << endl;
+ }
+ //else cout <<"pb avec les charges X1" << endl;
+ }
+ }
+
+ if(Chargex[StripMax_+1] !=-1) {
+ //cout << "no pb" << endl;
+ Q[2] = Chargex[StripMax_+1];
+ StripPos[2] = StripPositionX[ff][StripMax_ + 1][0];
+
+ // cout << "Q[2] = " << Q[2] << " Strip2 = " << StripMaxX[ff]+1 << " StripPos[2] =" << StripPos[2] << endl;
+ }
+ else {
+ if(Chargex[StripMax_ +2] !=-1) {
+ //cout << "coconutsX+1! " << endl;
+ Q[2] = Chargex[StripMax_+2];
+ StripPos[2] = StripPositionX[ff][StripMax_ + 2][0];
+ // cout << "Q[0] = " << Q[0] << " StripMaxX[ff]= " << StripMaxX[ff]<< " StripPos[0] =" << StripPos[0] << endl;
+ // cout << "Q[2] = " << Q[2] << " Strip2 = " << StripMaxX[ff]+2 << " StripPos[2] =" << StripPos[2] << endl;
+ }
+ else {
+ if(Chargex[StripMax_ +3] !=-1) {
+ //cout << "coconutsX+2! " << endl;
+ Q[2] = Chargex[StripMax_+3];
+ StripPos[2] = StripPositionX[ff][StripMax_ + 3][0];
+ // cout << "Q[0] = " << Q[0] << " StripMaxX[ff]= " << StripMaxX[ff]<< " StripPos[0] =" << StripPos[0] << endl;
+ // cout << "Q[2] = " << Q[2] << " Strip2 = " << StripMaxX[ff]+3 << " StripPos[2] =" << StripPos[2] << endl;
+ }
+ //else cout << "pb avec les charges X2" << endl; //Q[1] = 1 ; Strip[1] = 0 ;
+ }
+ }
+
+ double Q0_Q1 = log(Q[0]/Q[1]);
+ double Q0_Q2 = log(Q[0]/Q[2]);
+
+ double num = Q0_Q1 * (StripPos[2]*StripPos[2] - StripPos[0]*StripPos[0]) - Q0_Q2 * (StripPos[1]*StripPos[1] - StripPos[0]*StripPos[0]) ;
+ double denom = Q0_Q1 * (StripPos[2] - StripPos[0]) - Q0_Q2 * (StripPos[1] - StripPos[0]) ;
+ gauss = 0.5*num / denom;
+
+ }
+
+ else {
+ gauss = -1000;
+ // cout << "Gaussian method X failed ! " << endl;
+
+ }
+
+ return gauss;
+
+}
+
+
+double TCATSPhysics:: GaussianMethodY()
+//int StripMax)
+{
+ double Q[3];
+ double StripPos[3];
+
+ double Q0_Q1, Q0_Q2;
+ double num, denom;
+
+ for(int j = 0; j<3 ; j++)
+ {
+ Q[j] = 0;
+ StripPos[j] = 0;
+ }
+
+ int StripMax_ = StripMaxY[ff] - 1;
+
+ double gauss = -1000;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff]<27)
+ {
+ Q[0] = Chargey[StripMax_] ;
+ StripPos[0] = StripPositionY[ff][0][StripMax_];
+
+ if(Chargey[StripMax_-1] !=-1) {
+ //cout << "no pb" << endl;
+ Q[1] = Chargey[StripMax_-1];
+ StripPos[1] = StripPositionY[ff][0][StripMax_-1];
+ }
+ else {
+ if(Chargey[StripMax_-2] !=-1) {
+ //cout << "coconutsY-1! " << endl;
+ Q[1] = Chargey[StripMax_-2];
+ StripPos[1] = StripPositionY[ff][0][StripMax_ - 2] ;
+ }
+ else {
+ if(Chargey[StripMax_-3] !=-1) {
+ //cout << "coconutsY-2! " << endl;
+ Q[1] = Chargey[StripMax_-3];
+ StripPos[1] = StripPositionY[ff][0][StripMax_ - 3] ;
+ }
+ //else cout << "pb avec les charges Y1" << endl; //Q[1] = 1 ; Strip[1] = 0
+ }
+ }
+
+ if(Chargey[StripMax_+1] !=-1) {
+ //cout << "no pb" << endl;
+ Q[2] = Chargey[StripMax_+1];
+ StripPos[2] = StripPositionY[ff][0][StripMax_ + 1];
+ }
+ else {
+ if(Chargey[StripMax_ +2] !=-1) {
+ //cout << "coconutsY+1! " << endl;
+ Q[2] = Chargey[StripMax_+2];
+ StripPos[2] = StripPositionY[ff][0][StripMax_ + 2] ;
+ }
+ else {
+ if(Chargey[StripMax_ +3] !=-1) {
+ //cout << "coconutsY+2! " << endl;
+ Q[2] = Chargey[StripMax_+3];
+ StripPos[2] = StripPositionY[ff][0][StripMax_ + 3] ;
+ }
+ //else cout << "pb avec les charges Y2" << endl; //Q[1] = 1 ; Strip[1] = 0
+ }
+ }
+
+ //else cout << "pb avec les charges Y2" << endl; //Q[1] = 1 ; Strip[1] = 0
+
+
+ Q0_Q1 = log(Q[0]/Q[1]);
+ Q0_Q2 = log(Q[0]/Q[2]);
+
+ num = Q0_Q1 * (StripPos[2]*StripPos[2] - StripPos[0]*StripPos[0]) - Q0_Q2 * (StripPos[1]*StripPos[1] - StripPos[0]*StripPos[0]) ;
+ denom = Q0_Q1 * (StripPos[2] - StripPos[0]) - Q0_Q2 * (StripPos[1] - StripPos[0]) ;
+ gauss = 0.5*num / denom;
+
+ }
+
+ else {
+ gauss = -1000;
+ //cout << "Gaussian method Y failed ! " << endl;
+ }
+
+ return gauss;
+
+}
+
+
+double TCATSPhysics:: Barycentric5MethodX()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27)
+ {
+ int StripMax_ = StripMaxX[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+
+ for(int i = -2 ; i < 3 ; i++)
+ {
+ if(Chargex[StripMax_+i]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxX[ff]+i)*Chargex[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargex[StripMax_+i];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ; } // cout << "Yo" << endl ;}
+
+ }
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMaxX[ff] << endl;
+ }
+
+ return Barycenter ;
+}
+
+double TCATSPhysics:: Barycentric5MethodY()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27)
+ {
+ int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+
+ for(int i = -2 ; i < 3 ; i++)
+ {
+ if(Chargey[StripMax_+i]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxY[ff]+i)*Chargey[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargey[StripMax_+i];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ; } // cout << "Yo" << endl ;}
+
+ }
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMaxY[ff] << endl;
+ }
+
+ return Barycenter ;
+}
+
+double TCATSPhysics:: Barycentric3MethodX()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27)
+ {
+ int StripMax_ = StripMaxX[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+ for(int i = -1 ; i < 2 ; i++)
+ {
+ if(Chargex[StripMax_+i]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxX[ff]+i)*Chargex[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargex[StripMax_+i];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ;} // cout << "Yo" << endl ;}
+ }
+
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
+ }
+
+ return Barycenter ;
+}
+
+double TCATSPhysics:: Barycentric3MethodY()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27)
+ {
+ int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+ for(int i = -1 ; i < 2 ; i++)
+ {
+ if(Chargey[StripMax_+i]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxY[ff]+i)*Chargey[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargey[StripMax_+i];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ;} // cout << "Yo" << endl ;}
+ }
+
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
+ }
+
+ return Barycenter ;
+}
+
+
+double TCATSPhysics:: Barycentric4MethodX()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27) {
+
+ int StripMax_ = StripMaxX[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+ if(Chargex[StripMax_+1] > Chargex[StripMax_-1]) {
+
+ // cout << "Barycentre droit" << endl;
+ for(int i = -1 ; i < 3 ; i++)
+ {
+ if(Chargex[StripMax_+i]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxX[ff]+i)*Chargex[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargex[StripMax_+i];
+ }
+ }
+ }
+
+ else {
+ // cout << "Barycentre gauche" << endl;
+ for(int i = -2 ; i < 2 ; i++)
+ {
+ if(Chargex[StripMax_+i]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxX[ff]+i)*Chargex[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargex[StripMax_+i];
+ }
+ }
+ }
+
+ if(ChargeTotal>0) {
+ Barycenter = Barycenter / ChargeTotal ;
+ }
+
+ }
+
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
+ }
+
+ return Barycenter ;
+}
+
+
+double TCATSPhysics:: Barycentric4MethodY()
+{
+ double Barycenter = 0 ;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) {
+
+ int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+ if(Chargey[StripMax_+1] > Chargey[StripMax_-1]) {
+
+ // cout << "Barycentre droit" << endl;
+ for(int i = -1 ; i < 3 ; i++)
+ {
+ if(Chargey[StripMax_+i]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxY[ff]+i)*Chargey[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargey[StripMax_+i];
+ }
+ }
+ }
+
+ else {
+ // cout << "Barycentre gauche" << endl;
+ for(int i = -2 ; i < 2 ; i++)
+ {
+ if(Chargey[StripMax_+i]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxY[ff]+i)*Chargey[StripMax_+i] ;
+ NumberOfPoint++;
+ ChargeTotal+=Chargey[StripMax_+i];
+ }
+ }
+ }
+
+ if(ChargeTotal>0) {
+ Barycenter = Barycenter / ChargeTotal ;
+ }
+
+ }
+
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
+ }
+
+ return Barycenter ;
+}
+
+
+
+double TCATSPhysics::CorrectedPositionX3(double Position, double a)
+{
+ double Corrected_Position = 0;
+ int StripMax_ = StripMaxX[ff] -1;
+ double xmax = StripPositionX[ff][StripMax_][0] ;
+ // cout << "xmax2 " << xmax << endl;
+ Corrected_Position = (Position - xmax) / a + xmax;
+
+ return Corrected_Position;
+}
+
+double TCATSPhysics::CorrectedPositionY3(double Position, double a)
+{
+ double Corrected_Position = 0;
+ int StripMax_ = StripMaxY[ff] -1;
+ double xmax = StripPositionY[ff][0][StripMax_];
+ // cout << "xmax2 " << xmax << endl;
+ Corrected_Position = (Position - xmax) / a + xmax;
+
+ return Corrected_Position;
+}
+
+double TCATSPhysics::CorrectedPositionX4(double Position, double b)
+{
+ double Corrected_Position = 0;
+ double xmax = 0;
+ int StripMax_ = StripMaxX[ff] -1;
+
+ if(Chargex[StripMax_+1] > Chargex[StripMax_-1]) {
+ if(ff ==0) xmax = StripPositionX[ff][StripMax_][0] - 1.27;
+ else xmax = StripPositionX[ff][StripMax_][0] + 1.27;
+ }
+
+ else{
+ if(ff ==0) xmax = StripPositionX[ff][StripMax_-1][0] - 1.27;
+ else xmax = StripPositionX[ff][StripMax_-1][0] + 1.27;
+ }
+
+ Corrected_Position = (Position - xmax) / b + xmax;
+
+ return Corrected_Position;
+}
+
+double TCATSPhysics::CorrectedPositionY4(double Position, double b)
+{
+ double Corrected_Position = 0;
+ double xmax = 0;
+ int StripMax_ = StripMaxY[ff] -1;
+
+ if(Chargey[StripMax_+1] > Chargey[StripMax_-1]) {
+ xmax = StripPositionY[ff][0][StripMax_] + 1.27 ;
+ }
+
+ else{
+ xmax = StripPositionY[ff][0][StripMax_-1] + 1.27;
+ }
+
+ Corrected_Position = (Position - xmax) / b + xmax;
+
+ return Corrected_Position;
+}
+
+TVector3 TCATSPhysics::GetPositionOnTarget()
+{
+ TVector3 Position = TVector3 (GetPositionOnTargetX() ,
+ GetPositionOnTargetY() ,
+ 3);
+ return(Position) ;
+}
+
+double TCATSPhysics::GetCATSChargeSumX(int i)
+{
+ double ChargeSum = 0;
+
+ for(unsigned int i = 0; i < ChargeSumX.size(); i++)
+ {
+ if (DetNumberX.at(i) == 1) ChargeSum = ChargeSumX.at(i);
+ else if(DetNumberX.at(i) == 2) ChargeSum = ChargeSumX.at(i);
+ }
+
+ return(ChargeSum);
+}
+
+int TCATSPhysics::GetCATSMultOverThreshX(int i)
+{
+ int MultOverThresh = 0;
+
+ for(unsigned int i = 0; i < MultOverThreshX.size(); i++)
+ {
+ if (DetNumberX.at(i) == 1) MultOverThresh = MultOverThreshX.at(i);
+ else if(DetNumberX.at(i) == 2) MultOverThresh = MultOverThreshX.at(i);
+ }
+
+ return MultOverThresh;
+}
+
+int TCATSPhysics::GetCATSStripMaxX(int i)
+{
+ int StripMax = 0;
+
+ for(unsigned int i = 0; i < StripMaxX.size(); i++)
+ {
+ if (DetNumberX.at(i) == 1) StripMax = StripMaxX.at(i);
+ else if(DetNumberX.at(i) == 2) StripMax = StripMaxX.at(i);
+ }
+
+ return StripMax;
+}
+
+//
+/*
+int TCATSPhysics::GetCATSDetNumberX_Position(int i)
+{
+ int DetNumber_Position = 0;
+
+ for(unsigned int j = 0; j < DetNumberX_Position.size(); j++)
+ {
+ if (DetNumberX.at(j) == 1) DetNumber_Position = DetNumberX_Position.at(j);
+ else if(DetNumberX.at(j) == 2) DetNumber_Position = DetNumberX_Position.at(j);
+ }
+
+ return DetNumber_Position;
+}
+*/
+
+double TCATSPhysics::GetCATSPositionX(int i)
+{
+ double Position = 0;
+
+ for(unsigned int j = 0; j < PositionX.size(); j++)
+ {
+ if (DetNumberX_Position.at(j) == i) { Position = PositionX.at(j);}
+
+ }
+
+ return Position;
+}
+
+
+double TCATSPhysics::GetCATSChargeSumY(int i)
+{
+ double ChargeSum = 0;
+
+ for(unsigned int i = 0; i < ChargeSumY.size(); i++)
+ {
+ if (DetNumberY.at(i) == 1) ChargeSum = ChargeSumY.at(i);
+ else if(DetNumberY.at(i) == 2) ChargeSum = ChargeSumY.at(i);
+ }
+
+ return(ChargeSum);
+}
+
+int TCATSPhysics::GetCATSMultOverThreshY(int i)
+{
+ int MultOverThresh = 0;
+
+ for(unsigned int i = 0; i < MultOverThreshY.size(); i++)
+ {
+ if (DetNumberY.at(i) == 1) MultOverThresh = MultOverThreshY.at(i);
+ else if(DetNumberY.at(i) == 2) MultOverThresh = MultOverThreshY.at(i);
+ }
+
+ return MultOverThresh;
+}
+
+int TCATSPhysics::GetCATSStripMaxY(int i)
+{
+ int StripMax = 0;
+
+ for(unsigned int i = 0; i < StripMaxY.size(); i++)
+ {
+ if (DetNumberY.at(i) == 1) StripMax = StripMaxY.at(i);
+ else if(DetNumberY.at(i) == 2) StripMax = StripMaxY.at(i);
+ }
+
+ return StripMax;
+}
+/*
+int TCATSPhysics::GetCATSDetNumberY_Position(int i)
+{
+ int DetNumber_Position = 0;
+
+ for(unsigned int j = 0; j < DetNumberY_Position.size(); j++)
+ {
+ if (DetNumberY.at(j) == 1) DetNumber_Position = DetNumberY_Position.at(j);
+ else if(DetNumberY.at(i) == 2) DetNumber_Position = DetNumberY_Position.at(j);
+ }
+
+ return DetNumber_Position;
+}
+*/
+double TCATSPhysics::GetCATSPositionY(int i)
+{
+ double Position = 0;
+
+ for(unsigned int j = 0; j < PositionY.size(); j++)
+ {
+ if (DetNumberY_Position.at(j) == i) { Position = PositionY.at(j);}
+ }
+
+ return Position;
+}
+
+
+
+namespace LOCAL_CATS
+{
+ // tranform an integer to a string
+ string itoa(int value)
+ {
+ std::ostringstream o;
+
+ if (!(o << value))
+ return "" ;
+
+ return o.str();
+ }
+}
diff --git a/NPLib/CATS/TCATSPhysics.h b/NPLib/CATS/TCATSPhysics.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1fd2dbab9ae39d3d9c00db6a4427001701f0c30
--- /dev/null
+++ b/NPLib/CATS/TCATSPhysics.h
@@ -0,0 +1,264 @@
+#ifndef __CATSCAL__
+#define __CATSCAL__
+/*****************************************************************************
+ * Copyright (C) 2009-2010 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: Sandra Giron contact address: giron@ipno.in2p3.fr *
+ * *
+ * Creation Date : febuary 2010 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold CATS treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// STL
+#include <vector>
+
+// ROOT
+#include "TObject.h"
+#include "TVector3.h"
+
+// NPLib
+#include "TCATSData.h"
+#include "../include/VDetector.h"
+#include "../include/CalibrationManager.h"
+
+using namespace std ;
+
+/* J'aime pas trop cette partie, je pense que deja ca pourrait etre mieux une variable interne te disant quel methode tu as utiliser
+et d'ailleur d'ecrire cette varaible dans l'arbre de sortie pour une question de tracabilite.
+Ensuite tu peux faire un Set et un Get de cette variable (je preconise un string avec un nom completement lisible... :p ).
+Ensuite dans ton build tu appelle une methode unique, qui elle appellera la methode correcte apres avoir fait les tests...
+
+Si apres tu veux vraiment ameliorer les performances le mieux est de definir un pointer de fonction que tu appelle a chaque event...
+mais c'est un peu plus complique,
+voila un liens si jamais ca t'interresse: http://www.newty.de/fpt/intro.html
+
+Ceci dit ce n'est que des points de detail.
+*/
+
+enum reconstruction{NO,SECHS,GAUSS,BAR3,BAR4,BAR5};
+//enum correction{BAR3cor,BAR4cor,NOcor,GAUSScor};
+enum correction{NOcor,cor};
+
+class TCATSPhysics : public TObject, public NPA::VDetector
+{
+
+ public: // Constructor and Destructor
+ TCATSPhysics();
+ ~TCATSPhysics();
+
+ public: // Output data of interest
+ //for a CATS
+
+ // marker of the cats used
+ int ff ;
+
+ // Vector of dim = multiplicity of event on all detector
+ vector<int> DetNumberX ;
+ vector<int> StripX ;
+ vector<double> ChargeX ;
+
+ // Vector of dim = number of CATS
+ vector<double> ChargeSumX ;
+ vector<int> MultOverThreshX ;
+ vector<int> StripMaxX ;
+
+
+ // Vector of dim = multiplicity of event on all detector
+ vector<int> DetNumberY ;
+ vector<int> StripY ;
+ vector<double> ChargeY ;
+ // vector<double> ChargeY_test ;
+
+ // Vector of dim = number of CATS
+ vector<double> ChargeSumY ;
+ vector<int> MultOverThreshY ;
+ vector<int> StripMaxY ;
+ // vector<int> StripMaxY_test;
+
+ // Calculate
+ vector<int> DetNumberX_Position ;
+ vector<int> DetNumberY_Position ;
+ vector<int> DetNumberZ_Position ;
+ vector<double> PositionX ;
+ vector<double> PositionY ;
+ vector<double> PositionZ ;
+ double PositionOnTargetX ;
+ double PositionOnTargetY ;
+
+ TVector3 BeamDirection ; //!
+
+ double Chargex[28]; //!
+ double Chargey[28]; //!
+
+ int HitX; //!
+ int HitY; //!
+
+ vector<reconstruction> ReconstructionMethodX;
+ vector<reconstruction> ReconstructionMethodY;
+
+
+ // vector<reconstruction> FailedReconstructionX;
+ vector<reconstruction> FailedReconstructionY;
+
+ private: // Root Input and Output tree classes
+
+ TCATSData* EventData ;//!
+ TCATSPhysics* EventPhysics ;//!
+
+ 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 InitializeRootInput() ;
+
+ // Create associated branches and associated private member DetectorPhysics address
+ void InitializeRootOutput() ;
+
+ // This method is called at each event read from the Input Tree. Aim 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() ;
+
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {EventData->Clear();}
+
+
+ private :
+
+ // redundant information : could be optimized in the future
+ vector< vector< vector<double> > > StripPositionX ; //!
+ vector< vector< vector<double> > > StripPositionY ; //!
+ vector<double> StripPositionZ ; //!
+
+ int NumberOfCATS ; //!
+
+ vector< vector <double> > Pedestal_X ; //!
+ vector< vector <double> > Pedestal_Y ; //!
+
+ vector< vector <double> > Threshold_X ; //!
+ vector< vector <double> > Threshold_Y ; //!
+
+
+ public : // Specific to CATS
+
+ void Clear();
+ void Dump();
+
+ void AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28);
+
+ void ReadPedestal(string PedestalPath);
+
+ double AnalyseX(int ff,
+ int NumberOfDetector);
+
+ double AnalyseY(int ff,
+ int NumberOfDetector);
+
+ double CalculatePositionX( double CalculatedStripX,
+ correction method);
+
+ double CalculatePositionY( double CalculatedStripY,
+ correction method);
+
+
+ reconstruction ChooseReconstruction(TString type);
+
+ // Calculate Strip touch using an array of Charge on Strip and Strip with Maximum Charge
+
+ double HyperbolicSequentMethodX();
+ double GaussianMethodX();
+ double Barycentric5MethodX();
+ double Barycentric4MethodX();
+ double Barycentric3MethodX();
+
+ double HyperbolicSequentMethodY();
+ double GaussianMethodY();
+
+ double Barycentric5MethodY();
+ double Barycentric4MethodY();
+ double Barycentric3MethodY();
+
+
+ double CorrectedPositionX3(double Position, double a) ;
+ double CorrectedPositionY3(double Position, double a) ;
+ double CorrectedPositionX4(double Position, double b);
+ double CorrectedPositionY4(double Position, double b);
+
+ // X
+
+ // Vector of dim = multiplicity of event on all detector
+ int GetCATSDetNumberX(int i) {return DetNumberX.at(i);}
+ int GetCATSStripX(int i) {return StripX.at(i);}
+ double GetCATSChargeX(int i) {return ChargeX.at(i);}
+
+ int GetCATSMultX() {return DetNumberX.size();}
+
+ // Vector of dim = number of CATS
+ double GetCATSChargeSumX(int i) ;
+ int GetCATSMultOverThreshX(int i) ;
+ int GetCATSStripMaxX(int i) ;
+ // int GetCATSDetNumberX_Position(int i) ;
+ double GetCATSPositionX(int i) ;
+
+ double GetPositionOnTargetX() {return PositionOnTargetX;}
+
+ // Y
+
+ // Vector of dim = multiplicity of event on all detector
+ int GetCATSDetNumberY(int i) {return DetNumberY.at(i);}
+ int GetCATSStripY(int i) {return StripY.at(i);}
+ double GetCATSChargeY(int i) {return ChargeY.at(i);}
+
+ int GetCATSMultY() {return DetNumberY.size();}
+
+ // Vector of dim = number of CATS
+ double GetCATSChargeSumY(int i) ;
+ int GetCATSMultOverThreshY(int i) ;
+ int GetCATSStripMaxY(int i) ;
+ // int GetCATSDetNumberY_Position(int i);
+ double GetCATSPositionY(int i) ;
+
+ double GetPositionOnTargetY() {return PositionOnTargetY;}
+
+ int GetCATSMult() {return PositionX.size();}
+
+ TVector3 GetPositionOnTarget();
+ TVector3 GetBeamDirection() {return BeamDirection;}
+
+ ClassDef(TCATSPhysics,1) // CATSPhysics structure
+ };
+
+
+
+
+
+namespace LOCAL_CATS
+{
+ // tranform an integer to a string
+ string itoa(int value);
+
+}
+
+#endif