diff --git a/NPLib/CATS/TCATSData.h b/NPLib/CATS/TCATSData.h
index ec451d2aaa2470bd0a0cf2b08c53971a2ea30991..8af5e07ac5d4dd3e6c83f4b52f0ecf3016388e21 100644
--- a/NPLib/CATS/TCATSData.h
+++ b/NPLib/CATS/TCATSData.h
@@ -53,35 +53,35 @@ class TCATSData : public TObject {
///////////////////// 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);}
+ 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);}
+ 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);}
+ 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);}
+ UShort_t GetCATSMultX() const {return fCATS_DetX.size();}
+ UShort_t GetCATSDetX(Int_t i) const {return fCATS_DetX.at(i);}
+ UShort_t GetCATSStripX(Int_t i) const {return fCATS_StripX.at(i);}
+ UShort_t GetCATSChargeX(Int_t i) const {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);}
+ UShort_t GetCATSMultY() const {return fCATS_DetY.size();}
+ UShort_t GetCATSDetY(Int_t i) const {return fCATS_DetY.at(i);}
+ UShort_t GetCATSStripY(Int_t i) const {return fCATS_StripY.at(i);}
+ UShort_t GetCATSChargeY(Int_t i) const {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);}
+ UShort_t GetCATSMultQ() const {return fCATS_DetQ.size();}
+ UShort_t GetCATSDetQ(Int_t i) const {return fCATS_DetQ.at(i);}
+ UShort_t GetCATSCharge(Int_t i) const {return fCATS_Charge.at(i);}
- ClassDef(TCATSData,1) // CATSData structure
+ ClassDef(TCATSData,2) // CATSData structure
};
#endif
diff --git a/NPLib/CATS/TCATSPhysics.cxx b/NPLib/CATS/TCATSPhysics.cxx
index e5ab70bbeaed2b4c13f586ea830a05bb86c6773c..9e3990d5b1a0f494c8afc670764f73d6ac38fb63 100644
--- a/NPLib/CATS/TCATSPhysics.cxx
+++ b/NPLib/CATS/TCATSPhysics.cxx
@@ -9,7 +9,8 @@
* Original Author: Sandra Giron contact address: giron@ipno.in2p3.fr *
* *
* Creation Date : febuary 2010 *
- * Last update : *
+ * Last update : modification november 2011 by Pierre Morfouace *
+ * Contact adress : morfouac@ipno.in2p3.fr *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold CATS treated data *
@@ -19,7 +20,6 @@
* *
*****************************************************************************/
-
#include "TCATSPhysics.h"
using namespace LOCAL_CATS;
@@ -30,31 +30,239 @@ using namespace LOCAL_CATS;
#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 ;
+ m_EventData = new TCATSData ;
+ m_PreTreatedData = new TCATSData ;
+ m_EventPhysics = this ;
+ m_NumberOfCATS = 0 ;
}
+///////////////////////////////////////////////////////////////////////////
+TCATSPhysics::~TCATSPhysics()
+{
+}
+///////////////////////////////////////////////////////////////////////////
+void TCATSPhysics::PreTreat()
+{
+ ClearPreTreatedData();
+ gRandom->SetSeed(0);
+ // X
+ for(int i = 0 ; i < m_EventData->GetCATSMultX() ; i++)
+ {
+ // Valid Channel X
+ if(IsValidChannel("X", m_EventData->GetCATSDetX(i), m_EventData->GetCATSStripX(i)) )
+ {
+ if( fCATS_Threshold_X(m_EventData , i) )
+ {
+ double QX = fCATS_X_Q(m_EventData , i);
+ m_PreTreatedData->SetCATSChargeX( QX );
+ //Inversion X
+ if( *(m_CATSXInversion[m_EventData->GetCATSDetX(i)-1].begin() + m_EventData->GetCATSStripX(i)-1) != m_EventData->GetCATSStripX(i) )
+ {
+ m_PreTreatedData->SetCATSStripX( *(m_CATSXInversion[m_EventData->GetCATSDetX(i)-1].begin() + m_EventData->GetCATSStripX(i)-1) );
+ }
+ else {m_PreTreatedData->SetCATSStripX( m_EventData->GetCATSStripX(i) );}
+ m_PreTreatedData->SetCATSDetX( m_EventData->GetCATSDetX(i) );
+ }
+ }
+ }
+
+ // Y
+ for(int i = 0 ; i < m_EventData->GetCATSMultY() ; i++)
+ {
+ // Valid Channel Y
+ if(IsValidChannel("Y", m_EventData->GetCATSDetY(i), m_EventData->GetCATSStripY(i)))
+ {
+ if( fCATS_Threshold_Y(m_EventData , i) )
+ {
+ double QY = fCATS_Y_Q(m_EventData , i);
+ m_PreTreatedData->SetCATSChargeY( QY );
+ //Inversion Y
+ if( *(m_CATSYInversion[m_EventData->GetCATSDetY(i)-1].begin() + m_EventData->GetCATSStripY(i)-1) != m_EventData->GetCATSStripY(i) )
+ {
+ m_PreTreatedData->SetCATSStripY( *(m_CATSYInversion[m_EventData->GetCATSDetY(i)-1].begin() + m_EventData->GetCATSStripY(i)-1) );
+ }
+ else {m_PreTreatedData->SetCATSStripY( m_EventData->GetCATSStripY(i) );}
+ m_PreTreatedData->SetCATSDetY( m_EventData->GetCATSDetY(i) );
+ }
+ }
+ }
+ return;
+ }
-TCATSPhysics::~TCATSPhysics()
+/////////////////////////////////////////////////////////////////////////////
+void TCATSPhysics::BuildSimplePhysicalEvent()
{
+ BuildPhysicalEvent();
}
+//////////////////////////////////////////////////////////////////////////////
+void TCATSPhysics::BuildPhysicalEvent()
+{
+ PreTreat();
+ double Pi = 3.14159265;
+
+ // How many CATS?
+ int NumberOfCATSHit = 0 ;
+ int DetectorID = -1;
+ double SumChargeX[2];
+ double SumChargeY[2];
+
+
+ for( unsigned short i = 0 ; i < m_PreTreatedData->GetCATSMultX() ; i++ )
+ {
+ if( m_PreTreatedData->GetCATSDetX(i) != DetectorID) {
+ NumberOfCATSHit++;
+ }
+ if(NumberOfCATSHit == m_NumberOfCATS) break;
+ }
+
+
+ // INITIALISATION OF VECTORS : DIM = NUMBER OF CATS
+ for(int k = 0 ; k < NumberOfCATSHit ; k++ )
+ {
+ // X
+ StripMaxX.push_back(-1);
+ ReconstructionMethodX.push_back(NO);
+ SumChargeX[k] = 0;
+
+ // Y
+ StripMaxY.push_back(-1);
+ ReconstructionMethodY.push_back(NO);
+ SumChargeY[k] = 0;
+ }
+
+ for(int p = 0 ; p < m_NumberOfCATS ; p++){
+ for(int z=0; z<28; z++) {
+ Buffer_X_Q[z][p] = -1;
+ Buffer_Y_Q[z][p] = -1;
+ }
+ }
+
+
+ for(unsigned int i = 0 ; i < m_PreTreatedData->GetCATSMultX() ; i++ )
+ {
+ int StrX = m_PreTreatedData->GetCATSStripX(i);
+ int NX = m_PreTreatedData->GetCATSDetX(i);
+ double CATS_X_Q = m_PreTreatedData->GetCATSChargeX(i) ;
+ Buffer_X_Q[StrX-1][NX-1] = CATS_X_Q;
+ SumChargeX[NX-1] += CATS_X_Q;
+
+ ChargeX.push_back(CATS_X_Q);
+ StripX.push_back(StrX);
+ DetNumberX.push_back(NX);
+ HitX++;
+ if(HitX==1) StripMaxX[NX-1] = StrX;
+ else if(ChargeX[HitX-1] > Buffer_X_Q[StripMaxX[NX-1] -1][NX-1] ) StripMaxX[NX-1] = StrX ;
+ }
+
+ for(unsigned int j = 0 ; j < m_PreTreatedData->GetCATSMultY() ; j++ )
+ {
+ int StrY = m_PreTreatedData->GetCATSStripY(j);
+ int NY = m_PreTreatedData->GetCATSDetY(j);
+ double CATS_Y_Q = m_PreTreatedData->GetCATSChargeY(j) ;
+ Buffer_Y_Q[StrY-1][NY-1] = CATS_Y_Q;
+ SumChargeY[NY-1] += CATS_Y_Q;
+
+ ChargeY.push_back(CATS_Y_Q);
+ StripY.push_back(StrY);
+ DetNumberY.push_back(NY);
+ HitY++;
+ if(HitY==1) StripMaxY[NY-1] = StrY;
+ else if(ChargeY[HitY-1] > Buffer_Y_Q[StripMaxY[NY-1] -1][NY-1] ) StripMaxY[NY-1] = StrY ;
+ }
+
+ double CalculatedStripX = 0, CalculatedStripY = 0;
+ double posX = 0 , posY = 0;
+
+ for(ff = 0 ; ff < NumberOfCATSHit ; ff++ )
+ {
+ CalculatedStripX = AnalyseX(ff);
+ CalculatedStripY = AnalyseY(ff);
+
+ posX = CalculatePositionX(CalculatedStripX, NOcor);
+ posY = CalculatePositionY(CalculatedStripY, NOcor);
+
+ 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]) ;
+
+ QsumX.push_back(SumChargeX[ff]);
+ QsumY.push_back(SumChargeY[ff]);
+ }
+
+ if(NumberOfCATSHit > 1){
+ if(PositionX[0] != -1000 && PositionY[0] != -1000 && PositionX[1] != -1000 && PositionY[1] != -1000)
+ {
+ double PositionOnTargetX_1;
+ double PositionOnTargetY_1;
+ double l = sqrt((PositionZ[0]-PositionZ[1])*(PositionZ[0]-PositionZ[1]));
+ double L = - PositionZ[1];
+ double t = (l+L) / l;
+
+ PositionOnTargetX_1 = PositionX[0] + (PositionX[1] - PositionX[0]) * t ;
+ PositionOnTargetY_1 = PositionY[0] + (PositionY[1] - PositionY[0]) * t ;
+
+ if(m_TargetAngle != 0)
+ {
+ double a = (PositionZ[1]-PositionZ[0])/(PositionX[1]-PositionX[0]);
+ double b = PositionZ[0] - a*PositionX[0];
+ PositionOnTargetX = b/(tan(m_TargetAngle*Pi/180.) - a);
+
+ double t_new = (l + L + PositionOnTargetX*tan(m_TargetAngle*Pi/180.)) / l;
+ PositionOnTargetY = PositionY[0] + (PositionY[1] - PositionY[0]) * t_new ;
+ }
+
+
+ else{
+ PositionOnTargetX = PositionOnTargetX_1;
+ PositionOnTargetY = PositionOnTargetY_1;
+ }
+ GetPositionOnTarget();
+ GetBeamDirection();
+ }
+
+ else
+ {
+ BeamDirection = TVector3 ( 1 ,
+ 0 ,
+ 0 );
+
+ PositionOnTargetX = -1000 ;
+ PositionOnTargetY = -1000 ;
+ }
+ }
+
+ else if(NumberOfCATSHit == 1)
+ {
+ BeamDirection = TVector3 ( 1 ,
+ 0 ,
+ 0 );
+ PositionOnTargetX = -1000 ;
+ PositionOnTargetY = -1000 ;
+ }
+
+ return;
+
+}
+///////////////////////////////////////////////////////////////////////////
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
void TCATSPhysics::ReadConfiguration(string Path)
{
@@ -81,9 +289,8 @@ void TCATSPhysics::ReadConfiguration(string Path)
//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 ;
+ ReadingStatus = true ;
}
// Else don't toggle to Reading Block Status
@@ -100,7 +307,7 @@ void TCATSPhysics::ReadConfiguration(string Path)
// 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 ;
+ cout << "Warning: Another CATS is found before standard sequence of Token, Error may have occurred in CATS definition" << endl ;
ReadingStatus = false ;
}
@@ -119,7 +326,7 @@ void TCATSPhysics::ReadConfiguration(string Path)
Az = Az ;
A = TVector3(Ax, Ay, Az);
- cout << "X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
+ cout << " X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl;
}
else if (DataBuffer.compare(0, 7, "X28_Y1=") == 0) {
@@ -135,7 +342,7 @@ void TCATSPhysics::ReadConfiguration(string Path)
Bz = Bz ;
B = TVector3(Bx, By, Bz);
- cout << "X28 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
+ cout << " X28 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl;
}
else if (DataBuffer.compare(0, 7, "X1_Y28=") == 0) {
@@ -151,7 +358,7 @@ void TCATSPhysics::ReadConfiguration(string Path)
Cz = Cz ;
C = TVector3(Cx, Cy, Cz);
- cout << "X1 Y28 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
+ cout << " X1 Y28 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl;
}
else if (DataBuffer.compare(0, 8, "X28_Y28=") == 0) {
@@ -167,7 +374,7 @@ void TCATSPhysics::ReadConfiguration(string Path)
Dz = Dz ;
D = TVector3(Dx, Dy, Dz);
- cout << "X28 Y28 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
+ cout << " X28 Y28 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl;
}
@@ -194,31 +401,15 @@ void TCATSPhysics::ReadConfiguration(string Path)
}
}
-
- // Should be called in Analysis...
- ReadPedestal("../../../offline/calibrations/CATS/CATScoeff/Piedestaux_368.txt");
- // path given from e530 directory
-
- cout << endl << "/////////////////////////////" << endl<<endl;
-
+ InitializeStandardParameter();
+ ReadAnalysisConfig();
+ m_method_CATS1X = StringToEnum(m_reconstruction_CATS1X);
+ m_method_CATS1Y = StringToEnum(m_reconstruction_CATS1Y);
+ m_method_CATS2X = StringToEnum(m_reconstruction_CATS2X);
+ m_method_CATS2Y = StringToEnum(m_reconstruction_CATS2Y);
}
-// 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::InitializeRootInputRaw()
@@ -226,55 +417,56 @@ void TCATSPhysics::InitializeRootInputRaw()
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus( "CATS" , true ) ;
inputChain->SetBranchStatus( "fCATS_*" , true ) ;
- inputChain->SetBranchAddress( "CATS" , &EventData ) ;
+ inputChain->SetBranchAddress( "CATS" , &m_EventData ) ;
}
+
+/////////////////////////////////////////////////////////////////////
// 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 TCATSPhysics::InitializeRootInputPhysics()
{
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus( "CATS" , true );
- inputChain->SetBranchStatus( "ff" , true );
- inputChain->SetBranchStatus( "DetNumberX" , true );
- inputChain->SetBranchStatus( "StripX" , true );
- inputChain->SetBranchStatus( "ChargeX" , true );
- inputChain->SetBranchStatus( "ChargeSumX" , true );
- inputChain->SetBranchStatus( "MultOverThreshX" , true );
- inputChain->SetBranchStatus( "StripMaxX" , true );
- inputChain->SetBranchStatus( "DetNumberY" , true );
- inputChain->SetBranchStatus( "StripY" , true );
- inputChain->SetBranchStatus( "ChargeY" , true );
- inputChain->SetBranchStatus( "ChargeSumY" , true );
- inputChain->SetBranchStatus( "MultOverThreshY" , true );
- inputChain->SetBranchStatus( "StripMaxY" , true );
- inputChain->SetBranchStatus( "DetNumberX_Position" , true );
- inputChain->SetBranchStatus( "DetNumberY_Position" , true );
- inputChain->SetBranchStatus( "DetNumberZ_Position" , true );
- inputChain->SetBranchStatus( "PositionX" , true );
- inputChain->SetBranchStatus( "PositionY" , true );
- inputChain->SetBranchStatus( "PositionZ" , true );
- inputChain->SetBranchStatus( "PositionOnTargetX" , true );
- inputChain->SetBranchStatus( "PositionOnTargetY" , true );
- inputChain->SetBranchAddress( "CATS" , &EventPhysics );
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus( "CATS" , true );
+ inputChain->SetBranchStatus( "ff" , true );
+ inputChain->SetBranchStatus( "DetNumberX" , true );
+ inputChain->SetBranchStatus( "StripX" , true );
+ inputChain->SetBranchStatus( "ChargeX" , true );
+ inputChain->SetBranchStatus( "StripMaxX" , true );
+ inputChain->SetBranchStatus( "DetNumberY" , true );
+ inputChain->SetBranchStatus( "StripY" , true );
+ inputChain->SetBranchStatus( "ChargeY" , true );
+ inputChain->SetBranchStatus( "StripMaxY" , true );
+ inputChain->SetBranchStatus( "DetNumberX_Position" , true );
+ inputChain->SetBranchStatus( "DetNumberY_Position" , true );
+ inputChain->SetBranchStatus( "DetNumberZ_Position" , true );
+ inputChain->SetBranchStatus( "PositionX" , true );
+ inputChain->SetBranchStatus( "PositionY" , true );
+ inputChain->SetBranchStatus( "PositionZ" , true );
+ inputChain->SetBranchStatus( "PositionOnTargetX" , true );
+ inputChain->SetBranchStatus( "PositionOnTargetY" , true );
+ inputChain->SetBranchStatus( "QsumX" , true );
+ inputChain->SetBranchStatus( "QsumY" , true );
+ inputChain->SetBranchAddress( "CATS" , &m_EventPhysics );
}
+/////////////////////////////////////////////////////////////////////
// Create associated branches and associated private member DetectorPhysics address
void TCATSPhysics::InitializeRootOutput()
{
TTree* outputTree = RootOutput::getInstance()->GetTree() ;
- outputTree->Branch( "CATS" , "TCATSPhysics" , &EventPhysics ) ;
+ outputTree->Branch( "CATS" , "TCATSPhysics" , &m_EventPhysics ) ;
}
-
+/////////////////////////////////////////////////////////////////////
void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28)
{
- NumberOfCATS++ ;
+ m_NumberOfCATS++ ;
// remove warning
C_X28_Y28 *= 1;
- // Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are along X axis)
+ // Vector U on Telescope Face (paralelle to Y Strip)
TVector3 U = C_X28_Y1 - C_X1_Y1 ;
U = U.Unit() ;
@@ -287,8 +479,7 @@ void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y2
// Position Vector of X=1 Y=1 Strip
TVector3 Strip_1_1 ;
- // Geometry Parameter
-
+ // Geometry Parameters
double Face = 71.12 ; //mm
double NumberOfStrip = 28 ;
double StripPitch = Face / NumberOfStrip ; //mm
@@ -296,7 +487,6 @@ void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y2
// Buffer object to fill Position Array
vector<double> lineX ;
vector<double> lineY ;
- //vector<double> lineZ ;
vector< vector< double > > OneDetectorStripPositionX ;
vector< vector< double > > OneDetectorStripPositionY ;
@@ -305,23 +495,17 @@ void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y2
// 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() ) ;
+ lineY.push_back( StripCenter.y() ) ;
}
OneDetectorStripPositionX.push_back(lineX);
@@ -336,976 +520,649 @@ void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y2
}
-
-// 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()
-{
+///////////////////////////////////////////////////////////////
+void TCATSPhysics::Clear()
+{
+ DetNumberX.clear() ;
+ StripX.clear() ;
+ ChargeX.clear() ;
+ StripMaxX.clear() ;
+ DetNumberY.clear() ;
+ StripY.clear() ;
+ ChargeY.clear() ;
+ StripMaxY.clear() ;
+ DetNumberX_Position.clear() ;
+ DetNumberY_Position.clear() ;
+ DetNumberZ_Position.clear() ;
+ PositionX.clear() ;
+ PositionY.clear() ;
+ PositionZ.clear() ;
+ QsumX.clear();
+ QsumY.clear();
+ ReconstructionMethodX.clear() ;
+ ReconstructionMethodY.clear() ;
- // cout << EventData-> GetCATSMultX() << endl;
- // pourquoi pas une methode avec qui prend que le strip max par exemple...
+
+ ff = 0;
+ HitX = 0;
+ HitY = 0;
}
-*/
-
-////////////////////////////////////////// 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;
-}
+////////////////////////////////////////////////////////////////////////////
+bool TCATSPhysics :: IsValidChannel(const string DetectorType, const int Detector , const int channel)
+ {
+ if(DetectorType == "X")
+ return *(m_XChannelStatus[Detector-1].begin()+channel-1);
+
+ else if(DetectorType == "Y")
+ return *(m_YChannelStatus[Detector-1].begin()+channel-1);
+
+ else return false;
+ }
+
+///////////////////////////////////////////////////////////////////////////////////
+void TCATSPhysics::InitializeStandardParameter()
+ {
+ // Enable all channel and no inversion
+ vector< bool > ChannelStatus;
+ vector< int > InversionStatus;
+ m_XChannelStatus.clear() ;
+ m_YChannelStatus.clear() ;
+ m_CATSXInversion.clear() ;
+ m_CATSYInversion.clear() ;
+
+ ChannelStatus.resize(28,true);
+ InversionStatus.resize(28);
+ for(unsigned int j = 0 ; j < InversionStatus.size() ; j++)
+ {
+ InversionStatus[j] = j+1;
+ }
+
+ for(int i = 0 ; i < m_NumberOfCATS ; ++i)
+ {
+ m_XChannelStatus[i] = ChannelStatus;
+ m_YChannelStatus[i] = ChannelStatus;
+ m_CATSXInversion[i] = InversionStatus;
+ m_CATSYInversion[i] = InversionStatus;
+ }
+ return;
+ }
-void TCATSPhysics::Dump() const
+///////////////////////////////////////////////////////////////////////////
+void TCATSPhysics::ReadAnalysisConfig()
{
- cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
+ bool ReadingStatus = false;
- 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;
- }
+ // path to file
+ string FileName = "./configs/ConfigCATS.dat";
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigCATS.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigCATS.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigCATS.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, 10, "ConfigCATS") == 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 == "DISABLE_CHANNEL") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ int Detector = atoi(DataBuffer.substr(4,1).c_str());
+ int channel = -1;
+ if (DataBuffer.compare(5,4,"STRX") == 0) {
+ channel = atoi(DataBuffer.substr(9).c_str());
+ *(m_XChannelStatus[Detector-1].begin()+channel-1) = false;
+ }
+
+ else if (DataBuffer.compare(5,4,"STRY") == 0) {
+ channel = atoi(DataBuffer.substr(9).c_str());
+ *(m_YChannelStatus[Detector-1].begin()+channel-1) = false;
+ }
+
+ else cout << "Warning: detector type for CATS unknown!" << endl;
+
+ }
+
+ else if (whatToDo == "INVERSION") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer;
+ int Detector = atoi(DataBuffer.substr(4,1).c_str());
+ int channel1 = -1;
+ int channel2 = -1;
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer;
+ if (DataBuffer.compare(0,4,"STRX") == 0) {
+ channel1 = atoi(DataBuffer.substr(4).c_str());
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ channel2 = atoi(DataBuffer.substr(4).c_str());
+ *(m_CATSXInversion[Detector-1].begin()+channel1-1) = channel2;
+ *(m_CATSXInversion[Detector-1].begin()+channel2-1) = channel1;
+ }
+
+ else if (DataBuffer.compare(0,4,"STRY") == 0) {
+ channel1 = atoi(DataBuffer.substr(4).c_str());
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ channel2 = atoi(DataBuffer.substr(4).c_str());
+ *(m_CATSYInversion[Detector-1].begin()+channel1-1) = channel2;
+ *(m_CATSYInversion[Detector-1].begin()+channel2-1) = channel1;
+ }
+ }
+
+ else if (whatToDo == "RECONSTRUCTION_METHOD") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer;
+ if (DataBuffer.compare(0,6,"CATS1X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_reconstruction_CATS1X = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS1Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_reconstruction_CATS1Y = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_reconstruction_CATS2X = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_reconstruction_CATS2Y = DataBuffer;
+ }
+ }
+
+ else if (whatToDo == "CORRECTION_METHOD") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer;
+ if (DataBuffer.compare(0,6,"CATS1X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_correction_CATS1X = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS1Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_correction_CATS1Y = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_correction_CATS2X = DataBuffer;
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_correction_CATS2Y = DataBuffer;
+ }
+ }
+ /*else if (whatToDo == "CORRECTION_METHOD") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl;
+ m_correction = DataBuffer;
+ }*/
+
+ else if (whatToDo == "CORRECTION_COEF") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer;
+ if (DataBuffer.compare(0,6,"CATS1X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_CorrectionCoef_CATS1X = atof(DataBuffer.c_str());
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS1Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_CorrectionCoef_CATS1Y = atof(DataBuffer.c_str());
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2X") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_CorrectionCoef_CATS2X = atof(DataBuffer.c_str());
+ }
+
+ else if (DataBuffer.compare(0,6,"CATS2Y") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ cout << " " << DataBuffer << endl;
+ m_CorrectionCoef_CATS2Y = atof(DataBuffer.c_str());
+ }
+ }
+
+
+ else {ReadingStatus = false;}
+
+ }
+ }
+}
- 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++)
+/////////////////////////////////////////////////////////////////////
+// Add Parameter to the CalibrationManger
+void TCATSPhysics::AddParameterToCalibrationManager()
+{
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ for(int i = 0 ; i < m_NumberOfCATS ; i++)
{
- cout << "MultOverThreshY : " << MultOverThreshY.at(i) << endl;
+
+ 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") ;
+ Cal->AddParameter("CATS", "D"+itoa(i+1)+"_X"+itoa(j+1),"CATS_D"+itoa(i+1)+"_X"+itoa(j+1)) ;
+ Cal->AddParameter("CATS", "D"+itoa(i+1)+"_Y"+itoa(j+1),"CATS_D"+itoa(i+1)+"_Y"+itoa(j+1)) ;
+ }
}
-
-
- 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;
- }
+
+ return;
+}
- 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;
- }
+////////////////////////////////////////////////////////////////
+double TCATSPhysics::AnalyseX(int ff)
+{
+ double CalculatedStripX=0;
- for(unsigned int i = 0; i < ReconstructionMethodX.size() ; i++)
- {
- cout << "ReconstructionMethodX : " << ReconstructionMethodX.at(i) << endl;
- }
+ ReconstructionMethodX[ff] = ChooseReconstruction(ff,"X");
- 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;
- }
-
+ 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();
+ return(CalculatedStripX);
}
-
-
-
-void TCATSPhysics::ReadPedestal(string PedestalPath)
+////////////////////////////////////////////////////////////////
+double TCATSPhysics::AnalyseY(int ff)
{
- 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() ) {
+ double CalculatedStripY=0;
- Pedestal_File >> DataBuffer ;
-
- for(int i = 0 ; i < NumberOfCATS ; i++)
- {
- std::ostringstream DetectorNumber ;
- DetectorNumber << i+1;
+ ReconstructionMethodY[ff] = ChooseReconstruction(ff,"Y");
- 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()) ;
+ 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();
+
+
+ return(CalculatedStripY);
+}
- // 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()) ;
+////////////////////////////////////////////////////////////////////////
+reconstruction TCATSPhysics::ChooseReconstruction(int ff, TString type)
+{
+ reconstruction method = NO;
+ if(ff==0){
+ if(type=="X"){method = m_method_CATS1X;}
+ else if(type=="Y"){method = m_method_CATS1Y;}
+ }
- // 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
+ if(ff==1){
+ if(type=="X"){method = m_method_CATS2X;}
+ else if(type=="Y"){method = m_method_CATS2Y;}
+ }
- }
- }
- }
- }
+ return(method);
+}
+/////////////////////////////////////////////////////////////////////////
+reconstruction TCATSPhysics::StringToEnum(string type)
+{
+ reconstruction method = NO;
+ if(type=="GAUSS"){method = GAUSS;}
+ if(type=="SECHS"){method = SECHS;}
+ if(type=="BAR3"){method = BAR3;}
+ if(type=="BAR4"){method = BAR4;}
+ if(type=="BAR5"){method = BAR5;}
+
+ return(method);
}
-
-void TCATSPhysics::BuildSimplePhysicalEvent()
+//////////////////////////////////////////////////////////////////////
+double TCATSPhysics::CalculatePositionX(double CalculatedStripX, correction method)
{
- /*
- int HitX = 0 ;
- int HitY = 0 ;
- */
-
- //cout << "CATS Build... " << endl;
+ double positionX=-1000;
+ int IStripX = 0;
- gRandom->SetSeed(0);
+ if(ReconstructionMethodX[ff] == GAUSS){positionX = CalculatedStripX;} // already in mm -> see gaussian method
- //EventData->Dump();
- double CalculatedStripX, CalculatedStripY ;
- double posX =0 , posY = 0;
+ else{
+ IStripX = (int) CalculatedStripX ;
-
- // How many CATS?
- int NumberOfCATSHit = 0 ;
- int DetectorID = -1;
+ // Decimal Part
+ double DStripX = CalculatedStripX-IStripX ;
- 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;
-
- }
-
+ if( DStripX > 0.5) {IStripX++; DStripX = DStripX-1 ;} else {DStripX = DStripX;}
+
+ // Calculate Geometrical Position
+ if(IStripX > 0 && IStripX < 29){
+ 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
+
+ if(method == NOcor) positionX = positionX;
+ else if(method == cor){
+ if(ReconstructionMethodX[ff] == BAR3) positionX = CorrectedPositionX3(positionX, 0.60);
+ 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] ;
+
+ if(method == NOcor) positionX = positionX;
+ else if(method == cor){
+ if(ReconstructionMethodX[ff] == BAR3) positionX = CorrectedPositionX3(positionX, 0.30);
+ if(ReconstructionMethodX[ff] == BAR4) positionX = CorrectedPositionX4(positionX, 0.67);
+ }
+ }
+ else cout << "only 2CATS!! ff = " << ff << endl;
+ }
- // cout << "Nombre de CATS: " << NumberOfCATSHit << endl;
-
+ else {positionX = -1000;}
+ }
+
+
+ if(ff==0 && CalculatedStripX != -1000){
+ if(m_correction_CATS1X == "Correction3Points"){
+ positionX = Corrected3PointsX(positionX,m_CorrectionCoef_CATS1X);
+ }
+ if(m_correction_CATS1X == "Correction4Points"){
+ positionX = Corrected4PointsX(positionX,m_CorrectionCoef_CATS1X);
+ }
+ }
+
+ if(ff==1 && CalculatedStripX != -1000){
+ if(m_correction_CATS2X == "Correction3Points"){
+ positionX = Corrected3PointsX(positionX,m_CorrectionCoef_CATS2X);
+ }
+ if(m_correction_CATS2X == "Correction4Points"){
+ positionX = Corrected4PointsX(positionX,m_CorrectionCoef_CATS2X);
+ }
+ }
+
+
+ return positionX;
+}
- // 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;
+/////////////////////////////////////////////////////////////////////////
+double TCATSPhysics::CalculatePositionY(double CalculatedStripY, correction method)
+{
+ double positionY = -1000;
+ if(ReconstructionMethodY[ff] == GAUSS){positionY = CalculatedStripY;} // already in mm -> see gaussian method
-
- for(int gg = 0 ; gg < NumberOfCATSHit ; gg++ )
- {
- //int ff = NumberOfCATSHit - gg -1 ;
- ff = gg ;
+ else{
+ // Integer part
+ int IStripY = (int) CalculatedStripY ;
- 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);
+ // Decimal Part
+ double DStripY = CalculatedStripY-IStripY ;
- 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;
+ if(DStripY > 0.5) {IStripY++; DStripY = DStripY-1;}
+ else {DStripY = DStripY; }
+
+ // Calculate Geometrical Position
+ if(IStripY > 0 && IStripY < 29 ){
+ positionY = (DStripY)*2.54 + StripPositionY[ff][0][IStripY-1] ; // conversion en mm initiale
+ 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);
+ }
+ }
- PositionOnTargetX = PositionX[0] + (PositionX[1] - PositionX[0]) * t ;
- PositionOnTargetY = PositionY[0] + (PositionY[1] - PositionY[0]) * t ;
+ else if(ff ==1){
+ if(method == NOcor) positionY = positionY;
+ else if(method == cor){
+ if(ReconstructionMethodY[ff] == BAR3) positionY = CorrectedPositionY3(positionY, 0.45);
+ if(ReconstructionMethodY[ff] == BAR4) positionY = CorrectedPositionY4(positionY, 0.7);
+ }
+ }
- //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 << "only 2CATS!! ff = " << ff << endl;
+ }
+
+ else {positionY = -1000;}
}
-
- 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 ;
+
+
+ if(ff==0 && CalculatedStripY != -1000){
+ if(m_correction_CATS1Y == "Correction3Points"){
+ positionY = Corrected3PointsY(positionY,m_CorrectionCoef_CATS1Y);
+ }
+ if(m_correction_CATS1Y == "Correction4Points"){
+ positionY = Corrected4PointsY(positionY,m_CorrectionCoef_CATS1Y);
+ }
}
- }
-
- 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 ;
+ if(ff==1 && CalculatedStripY != -1000){
+ if(m_correction_CATS2Y == "Correction3Points"){
+ positionY = Corrected3PointsY(positionY,m_CorrectionCoef_CATS2Y);
+ }
+ if(m_correction_CATS2Y == "Correction4Points"){
+ positionY = Corrected4PointsY(positionY,m_CorrectionCoef_CATS2Y);
}
-
-
- 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);
+ 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];
-
+ int StripMax = StripMaxX[ff];
+ double gauss = -1000;
+ double Q[3];
+ double StripPos[3];
- for(int j = 0; j<3 ; j++)
- {
- Q[j] = 0;
- StripPos[j] = 0;
- }
+ 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];
+ if(StripMaxX[ff]> 3 && StripMaxX[ff]< 26)
+ {
+ Q[0] = Buffer_X_Q[StripMax-1][ff] ;
+ StripPos[0] = StripPositionX[ff][StripMax-1][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];
+ if(Buffer_X_Q[StripMax-2][ff]!=-1){
+ Q[1] = Buffer_X_Q[StripMax-2][ff];
+ StripPos[1] = StripPositionX[ff][StripMax-2][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;
+ else {
+ if(Buffer_X_Q[StripMax-3][ff]!=-1){
+ Q[1] = Buffer_X_Q[StripMax-3][ff];
+ StripPos[1] = StripPositionX[ff][StripMax-3][0];
+ }
+ else {
+ if(Buffer_X_Q[StripMax-4][ff]!=-1){
+ Q[1] = Buffer_X_Q[StripMax-4][ff];
+ StripPos[1] = StripPositionX[ff][StripMax-4][0];
+ }
+ }
}
- }
-
- 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 ;
+
+ if(Buffer_X_Q[StripMax][ff]!=-1){
+ Q[2] = Buffer_X_Q[StripMax][ff];
+ StripPos[2] = StripPositionX[ff][StripMax][0];
+ }
+ else {
+ if(Buffer_X_Q[StripMax+1][ff]!=-1){
+ Q[2] = Buffer_X_Q[StripMax+1][ff];
+ StripPos[2] = StripPositionX[ff][StripMax+1][0];
+ }
+ else {
+ if(Buffer_X_Q[StripMax+2][ff]!=-1){
+ Q[2] = Buffer_X_Q[StripMax+2][ff];
+ StripPos[2] = StripPositionX[ff][StripMax+2][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;
+ if(denom != 0){
+ gauss = 0.5*num / denom;
+ }
+ else{gauss = -1000;}
}
else {
gauss = -1000;
- // cout << "Gaussian method X failed ! " << endl;
-
}
return gauss;
}
+/////////////////////////////////////////////////////////////////////////
+double TCATSPhysics::Corrected3PointsX(double Position, double c)
+{
+ double Corrected_Position = 0;
+ int StripMax_ = StripMaxX[ff] -1;
+ double xmax = StripPositionX[ff][StripMax_][0] ;
+
+ Corrected_Position = (Position - xmax) / c + xmax;
+
+ return Corrected_Position;
+}
-double TCATSPhysics:: GaussianMethodY()
-//int StripMax)
+/////////////////////////////////////////////////////////////////////////
+double TCATSPhysics::Corrected4PointsX(double Position, double d)
{
- double Q[3];
- double StripPos[3];
+ double Corrected_Position = 0;
+ double xmax = 0;
+ int StripMax_ = StripMaxX[ff] -1;
+
+ if(Buffer_X_Q[StripMax_+1][ff] > Buffer_X_Q[StripMax_-1][ff]) {
+ 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) / d + xmax;
+
+ return Corrected_Position;
+}
- double Q0_Q1, Q0_Q2;
- double num, denom;
+////////////////////////////////////////////////////////////////////////////
+double TCATSPhysics:: GaussianMethodY()
+{
+ 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;
- }
+ for(int j = 0; j<3 ; j++)
+ {
+ Q[j] = 0;
+ StripPos[j] = 0;
+ }
- int StripMax_ = StripMaxY[ff] - 1;
+ int StripMax = StripMaxY[ff];
- double gauss = -1000;
+ 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(StripMaxY[ff] > 2 && StripMaxY[ff]<27)
+ {
+ Q[0] = Buffer_Y_Q[StripMax-1][ff] ;
+ StripPos[0] = StripPositionY[ff][0][StripMax-1];
- 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
+ if(Buffer_Y_Q[StripMax-2][ff]!=-1){
+ Q[1] = Buffer_Y_Q[StripMax-2][ff];
+ StripPos[1] = StripPositionY[ff][0][StripMax-2];
+ }
+
+ else {
+ if(Buffer_Y_Q[StripMax-3][ff]!=-1){
+ Q[1] = Buffer_Y_Q[StripMax-3][ff];
+ StripPos[1] = StripPositionY[ff][0][StripMax-3] ;
+ }
+ else {
+ if(Buffer_Y_Q[StripMax-4][ff]!=-1){
+ Q[1] = Buffer_Y_Q[StripMax-4][ff];
+ StripPos[1] = StripPositionY[ff][0][StripMax-4] ;
+ }
+ }
}
- }
+
+ if(Buffer_Y_Q[StripMax][ff]!=-1){
+ Q[2] = Buffer_Y_Q[StripMax][ff];
+ StripPos[2] = StripPositionY[ff][0][StripMax];
+ }
+
+ else {
+ if(Buffer_Y_Q[StripMax+1][ff]!=-1){
+ Q[2] = Buffer_Y_Q[StripMax+1][ff];
+ StripPos[2] = StripPositionY[ff][0][StripMax+1] ;
+ }
+
+ else {
+ if(Buffer_Y_Q[StripMax+2][ff]!=-1){
+ Q[2] = Buffer_Y_Q[StripMax+2][ff];
+ StripPos[2] = StripPositionY[ff][0][StripMax+2] ;
+ }
+ }
+ }
- //else cout << "pb avec les charges Y2" << endl; //Q[1] = 1 ; Strip[1] = 0
Q0_Q1 = log(Q[0]/Q[1]);
@@ -1313,72 +1170,105 @@ double TCATSPhysics:: GaussianMethodY()
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;
+ if(denom != 0){
+ gauss = 0.5*num / denom;
+ }
}
else {
gauss = -1000;
- //cout << "Gaussian method Y failed ! " << endl;
}
return gauss;
}
+/////////////////////////////////////////////////////////////////////////
+double TCATSPhysics::Corrected3PointsY(double Position, double c)
+{
+ double Corrected_Position = 0;
+ int StripMax_ = StripMaxY[ff] -1;
+ double ymax = StripPositionY[ff][0][StripMax_] ;
+
+ Corrected_Position = (Position - ymax) / c + ymax;
+
+ return Corrected_Position;
+}
+
+/////////////////////////////////////////////////////////////////////////
+double TCATSPhysics::Corrected4PointsY(double Position, double d)
+{
+ double Corrected_Position = 0;
+ double ymax = 0;
+ int StripMax_ = StripMaxY[ff] -1;
+
+ if(Buffer_Y_Q[StripMax_+1][ff] > Buffer_Y_Q[StripMax_-1][ff]) {
+ ymax = StripPositionY[ff][0][StripMax_] + 1.27 ;
+ }
+
+ else{
+ ymax = StripPositionY[ff][0][StripMax_-1] + 1.27;
+ }
+
+ Corrected_Position = (Position - ymax) / d + ymax;
+
+ return Corrected_Position;
+}
+///////////////////////////////////////////////////////////////
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 ;
+ int StripMax_ = StripMaxX[ff] -1 ;
double NumberOfPoint = 0 ;
double ChargeTotal =0;
for(int i = -2 ; i < 3 ; i++)
{
- if(Chargex[StripMax_+i]!=-1) // Charge initialized to -1
+ if(Buffer_X_Q[StripMax_+i][ff]!=-1)
{
- Barycenter += (StripMaxX[ff]+i)*Chargex[StripMax_+i] ;
+ Barycenter += (StripMaxX[ff]+i)*Buffer_X_Q[StripMax_+i][ff] ;
NumberOfPoint++;
- ChargeTotal+=Chargex[StripMax_+i];
+ ChargeTotal+=Buffer_X_Q[StripMax_+i][ff];
}
}
if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
- else {Barycenter = -1000 ; } // cout << "Yo" << endl ;}
+ else {Barycenter = -1000 ; }
}
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 ;
+ int StripMax_ = StripMaxY[ff] -1 ;
double NumberOfPoint = 0 ;
double ChargeTotal =0;
for(int i = -2 ; i < 3 ; i++)
{
- if(Chargey[StripMax_+i]!=-1) // Charge initialized to -1
+ if(Buffer_Y_Q[StripMax_+i][ff]!=-1)
{
- Barycenter += (StripMaxY[ff]+i)*Chargey[StripMax_+i] ;
+ Barycenter += (StripMaxY[ff]+i)*Buffer_Y_Q[StripMax_+i][ff] ;
NumberOfPoint++;
- ChargeTotal+=Chargey[StripMax_+i];
+ ChargeTotal+=Buffer_Y_Q[StripMax_+i][ff];
}
}
@@ -1389,211 +1279,281 @@ double TCATSPhysics:: Barycentric5MethodY()
else
{
Barycenter = -1000 ;
- // cout << "Strip max " << StripMaxY[ff] << endl;
}
return Barycenter ;
}
+///////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////
double TCATSPhysics:: Barycentric3MethodX()
{
- double Barycenter = 0 ;
-
- if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27)
+ 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 ;}
+ int StripMax_ = StripMaxX[ff] -1;
+ double NumberOfPoint = 0 ;
+ double ChargeTotal =0;
+
+ for(int i = -1 ; i < 2 ; i++)
+ {
+ if(Buffer_X_Q[StripMax_+i][ff]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxX[ff]+i)*Buffer_X_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_X_Q[StripMax_+i][ff];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ;} // cout << "Yo" << endl ;}
}
-
- else
+
+
+ else
{
- Barycenter = -1000 ;
- // cout << "Strip max " << StripMax << endl;
+ Barycenter = -1000 ;
}
-
- return Barycenter ;
+
+ return Barycenter ;
}
+///////////////////////////////////////////////////////////////
double TCATSPhysics:: Barycentric3MethodY()
{
- double Barycenter = 0 ;
-
- if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27)
+ 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 ;}
+ 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(Buffer_Y_Q[StripMax_+i][ff]!=-1) // Charge initialized to -1
+ {
+ Barycenter += (StripMaxY[ff]+i)*Buffer_Y_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_Y_Q[StripMax_+i][ff];
+ }
+ }
+
+ if(ChargeTotal>0) Barycenter = Barycenter / ChargeTotal ;
+ else {Barycenter = -1000 ;} // cout << "Yo" << endl ;}
}
-
- else
+
+ else
{
- Barycenter = -1000 ;
- // cout << "Strip max " << StripMax << endl;
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
}
-
- return Barycenter ;
+
+ 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];
- }
+ 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(Buffer_X_Q[StripMax_+1][ff] > Buffer_X_Q[StripMax_-1][ff]) {
+
+ // cout << "Barycentre droit" << endl;
+ for(int i = -1 ; i < 3 ; i++)
+ {
+ if(Buffer_X_Q[StripMax_+i][ff]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxX[ff]+i)*Buffer_X_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_X_Q[StripMax_+i][ff];
+ }
+ }
+ }
+
+ else {
+ // cout << "Barycentre gauche" << endl;
+ for(int i = -2 ; i < 2 ; i++)
+ {
+ if(Buffer_X_Q[StripMax_+i][ff]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxX[ff]+i)*Buffer_X_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_X_Q[StripMax_+i][ff];
+ }
+ }
+ }
+
+ if(ChargeTotal>0) {
+ Barycenter = Barycenter / ChargeTotal ;
+ }
+
}
- }
-
- if(ChargeTotal>0) {
- Barycenter = Barycenter / ChargeTotal ;
- }
-
- }
-
- else
+
+ else
{
- Barycenter = -1000 ;
- // cout << "Strip max " << StripMax << endl;
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
}
-
- return Barycenter ;
+
+ 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];
- }
+ 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(Buffer_Y_Q[StripMax_+1][ff] > Buffer_Y_Q[StripMax_-1][ff]) {
+
+ // cout << "Barycentre droit" << endl;
+ for(int i = -1 ; i < 3 ; i++)
+ {
+ if(Buffer_Y_Q[StripMax_+i][ff]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxY[ff]+i)*Buffer_Y_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_Y_Q[StripMax_+i][ff];
+ }
+ }
+ }
+
+ else {
+ // cout << "Barycentre gauche" << endl;
+ for(int i = -2 ; i < 2 ; i++)
+ {
+ if(Buffer_Y_Q[StripMax_+i][ff]!=-1) { // Charge initialized to -1
+ Barycenter += (StripMaxY[ff]+i)*Buffer_Y_Q[StripMax_+i][ff] ;
+ NumberOfPoint++;
+ ChargeTotal+=Buffer_Y_Q[StripMax_+i][ff];
+ }
+ }
+ }
+
+ if(ChargeTotal>0) {
+ Barycenter = Barycenter / ChargeTotal ;
+ }
+
}
+
+ else
+ {
+ Barycenter = -1000 ;
+ // cout << "Strip max " << StripMax << endl;
}
+
+ return Barycenter ;
+}
- if(ChargeTotal>0) {
- Barycenter = Barycenter / ChargeTotal ;
+/////////////////////////////////////////////////////////////////////
+double TCATSPhysics:: HyperbolicSequentMethodX()
+{
+ double sechs = -1000 ;
+
+ if(StripMaxX[ff] > 2 && StripMaxX[ff]<27)
+ {
+ double vs1 = sqrt( Buffer_X_Q[StripMaxX[ff]-1][ff]/Buffer_X_Q[StripMaxX[ff]-1+1][ff] ) ;
+ double vs2 = sqrt( Buffer_X_Q[StripMaxX[ff]-1][ff]/Buffer_X_Q[StripMaxX[ff]-1-1][ff] ) ;
+ 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 ( Buffer_X_Q[StripMaxX[ff]-1+1][ff]>Buffer_X_Q[StripMaxX[ff]-1-1][ff] )
+ { sechs = StripMaxX[ff] + vs6/vs4 ;}
+
+ else
+ { sechs = StripMaxX[ff] - vs6/vs4 ;}
+
}
-
- }
+
+ else {
+ sechs = -1000;
+ }
+
+ return sechs ;
+}
- else
+//////////////////////////////////////////////////////////////////
+double TCATSPhysics:: HyperbolicSequentMethodY()
+{
+ double sechs = -1000 ;
+
+ if(StripMaxY[ff] > 2 && StripMaxY[ff]<27)
{
- Barycenter = -1000 ;
- // cout << "Strip max " << StripMax << endl;
- }
-
- return Barycenter ;
+ double vs1 = sqrt( Buffer_Y_Q[StripMaxY[ff]-1][ff]/Buffer_Y_Q[StripMaxY[ff]-1+1][ff] ) ;
+ double vs2 = sqrt( Buffer_Y_Q[StripMaxY[ff]-1][ff]/Buffer_Y_Q[StripMaxY[ff]-1-1][ff] ) ;
+ 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 ( Buffer_Y_Q[StripMaxY[ff]-1+1][ff]>Buffer_Y_Q[StripMaxY[ff]-1-1][ff] )
+ { sechs = StripMaxY[ff] + vs6/vs4 ;}
+
+ else
+ { sechs = StripMaxY[ff] - vs6/vs4 ;}
+
+ }
+
+ else {
+ sechs = -1000;
+ }
+
+ return sechs ;
}
+///////////////////////////////////////////////////////////////
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(Buffer_X_Q[StripMax_+1][ff] > Buffer_X_Q[StripMax_-1][ff]) {
if(ff ==0) xmax = StripPositionX[ff][StripMax_][0] - 1.27;
else xmax = StripPositionX[ff][StripMax_][0] + 1.27;
}
@@ -1608,13 +1568,14 @@ double TCATSPhysics::CorrectedPositionX4(double Position, double b)
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]) {
+ if(Buffer_Y_Q[StripMax_+1][ff] > Buffer_Y_Q[StripMax_-1][ff]) {
xmax = StripPositionY[ff][0][StripMax_] + 1.27 ;
}
@@ -1627,156 +1588,85 @@ double TCATSPhysics::CorrectedPositionY4(double Position, double b)
return Corrected_Position;
}
-TVector3 TCATSPhysics::GetPositionOnTarget()
-{
- TVector3 Position = TVector3 (GetPositionOnTargetX() ,
- GetPositionOnTargetY() ,
- 3);
- return(Position) ;
-}
-double TCATSPhysics::GetCATSChargeSumX(int i)
+///////////////////////////////////////////////////////////////
+TVector3 TCATSPhysics::GetBeamDirection()
{
- 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;
+ TVector3 Position = TVector3 (PositionX[1]-PositionX[0] ,
+ PositionY[1]-PositionY[0] ,
+ PositionZ[1]-PositionZ[0] );
+ Position.Unit();
+ return(Position) ;
}
-//
-/*
-int TCATSPhysics::GetCATSDetNumberX_Position(int i)
+///////////////////////////////////////////////////////////////
+TVector3 TCATSPhysics::GetPositionOnTarget()
{
- 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 Pi = 3.14159265;
+ TVector3 Position = TVector3 (GetPositionOnTargetX() ,
+ GetPositionOnTargetY() ,
+ GetPositionOnTargetX()*tan(m_TargetAngle*Pi/180));
+ Position.Unit();
+ return(Position) ;
}
-*/
-double TCATSPhysics::GetCATSPositionX(int i)
+////////////////////////////////////////////////////////////////////////
+namespace LOCAL_CATS
{
- double Position = 0;
-
- for(unsigned int j = 0; j < PositionX.size(); j++)
- {
- if (DetNumberX_Position.at(j) == i) { Position = PositionX.at(j);}
-
- }
-
- return Position;
-}
+ // transform an integer to a string
+ string itoa(int value)
+ {
+ std::ostringstream o;
+
+ if (!(o << value))
+ return "" ;
+
+ return o.str();
+ }
+ double fCATS_X_Q(const TCATSData* m_EventData , const int i)
+ {
+ return CalibrationManager::getInstance()->ApplyCalibration( "CATS/D" + itoa( m_EventData->GetCATSDetX(i) ) + "_X" + itoa( m_EventData->GetCATSStripX(i) ) + "_Q",
+ m_EventData->GetCATSChargeX(i) + gRandom->Rndm() - fCATS_Ped_X(m_EventData, i) );
+ }
+
+ double fCATS_Y_Q(const TCATSData* m_EventData , const int i)
+ {
+ return CalibrationManager::getInstance()->ApplyCalibration( "CATS/D" + itoa( m_EventData->GetCATSDetY(i) ) + "_Y" + itoa( m_EventData->GetCATSStripY(i) ) + "_Q",
+ m_EventData->GetCATSChargeY(i) + gRandom->Rndm() - fCATS_Ped_Y(m_EventData, i) );
+ }
-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);
- }
+ bool fCATS_Threshold_X(const TCATSData* m_EventData , const int i)
+ {
+ return CalibrationManager::getInstance()->ApplyThreshold( "CATS/D" + itoa( m_EventData->GetCATSDetX(i) ) + "_X" + itoa( m_EventData->GetCATSStripX(i) ),
+ m_EventData->GetCATSChargeX(i));
+ }
- return(ChargeSum);
+ bool fCATS_Threshold_Y(const TCATSData* m_EventData , const int i)
+ {
+ return CalibrationManager::getInstance()->ApplyThreshold( "CATS/D" + itoa( m_EventData->GetCATSDetY(i) ) + "_Y" + itoa( m_EventData->GetCATSStripY(i) ),
+ m_EventData->GetCATSChargeY(i));
+ }
+
+ double fCATS_Ped_X(const TCATSData* m_EventData, const int i)
+ {
+ return CalibrationManager::getInstance()->GetPedestal( "CATS/D" + itoa( m_EventData->GetCATSDetX(i) ) + "_X" + itoa( m_EventData->GetCATSStripX(i) ) );
+ }
+
+ double fCATS_Ped_Y(const TCATSData* m_EventData, const int i)
+ {
+ return CalibrationManager::getInstance()->GetPedestal( "CATS/D" + itoa( m_EventData->GetCATSDetY(i) ) + "_Y" + itoa( m_EventData->GetCATSStripY(i) ) );
+ }
+
+
}
-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)
- {
- char buffer [33];
- sprintf(buffer,"%d",value);
- return buffer;
- }
-}
diff --git a/NPLib/CATS/TCATSPhysics.h b/NPLib/CATS/TCATSPhysics.h
index 85ae1bc587d6688638d0b99f638946b6aa4f52aa..3093d666965de8aed3c341235e9076aa96a0914f 100644
--- a/NPLib/CATS/TCATSPhysics.h
+++ b/NPLib/CATS/TCATSPhysics.h
@@ -23,32 +23,24 @@
// STL
#include <vector>
-
// ROOT
#include "TObject.h"
#include "TVector3.h"
-
+#include <TRandom1.h>
+#include <TRandom2.h>
+#include <TRandom3.h>
// NPLib
#include "TCATSData.h"
#include "../include/VDetector.h"
#include "../include/CalibrationManager.h"
+#include "../include/DetectorManager.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
+#define NBDETECTOR 2
+#define NBSTRIPS 28
-Ceci dit ce n'est que des points de detail.
-*/
+using namespace std ;
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
@@ -57,11 +49,15 @@ 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
+ private: // Root Input and Output tree classes
+
+ TCATSData* m_EventData;//!
+ TCATSData* m_PreTreatedData;//!
+ TCATSPhysics* m_EventPhysics;//!
+
+ public :
+ // marker of the cats used
int ff ;
// Vector of dim = multiplicity of event on all detector
@@ -70,8 +66,6 @@ class TCATSPhysics : public TObject, public NPA::VDetector
vector<double> ChargeX;
// Vector of dim = number of CATS
- vector<double> ChargeSumX;
- vector<int> MultOverThreshX;
vector<int> StripMaxX;
@@ -79,28 +73,26 @@ class TCATSPhysics : public TObject, public NPA::VDetector
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 of dim = number of CATS
vector<int> StripMaxY;
- // vector<int> StripMaxY_test;
- // Calculate
+ // Vector of dim = number of CATS
vector<int> DetNumberX_Position;
vector<int> DetNumberY_Position;
vector<int> DetNumberZ_Position;
vector<double> PositionX;
vector<double> PositionY;
vector<double> PositionZ;
+ vector<double> QsumX;
+ vector<double> QsumY;
double PositionOnTargetX;
double PositionOnTargetY;
TVector3 BeamDirection ; //!
- double Chargex[28]; //!
- double Chargey[28]; //!
+ double Buffer_X_Q[NBSTRIPS][NBDETECTOR];//!
+ double Buffer_Y_Q[NBSTRIPS][NBDETECTOR];//!
int HitX; //!
int HitY; //!
@@ -108,156 +100,136 @@ class TCATSPhysics : public TObject, public NPA::VDetector
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 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. 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();
+
+ private :
+ vector< vector< vector<double> > > StripPositionX;//!
+ vector< vector< vector<double> > > StripPositionY;//!
+ vector<double> StripPositionZ;//!
+ int m_NumberOfCATS;
+ double m_TargetAngle;
+ double m_TargetThickness;
+ double m_CorrectionCoef_CATS1X;//!
+ double m_CorrectionCoef_CATS1Y;//!
+ double m_CorrectionCoef_CATS2X;//!
+ double m_CorrectionCoef_CATS2Y;//!
+
+
+ string m_correction_CATS1X;//!
+ string m_correction_CATS1Y;//!
+ string m_correction_CATS2X;//!
+ string m_correction_CATS2Y;//!
+
+ string m_reconstruction_CATS1X;//!
+ string m_reconstruction_CATS1Y;//!
+ string m_reconstruction_CATS2X;//!
+ string m_reconstruction_CATS2Y;//!
+ reconstruction m_method_CATS1X;//!
+ reconstruction m_method_CATS1Y;//!
+ reconstruction m_method_CATS2X;//!
+ reconstruction m_method_CATS2Y;//!
+
+ private :
+ // Map of activated channel
+ map< int, vector<bool> > m_XChannelStatus;//!
+ map< int, vector<bool> > m_YChannelStatus;//!
+ // Map of inverted channel
+ map< int, vector<int> > m_CATSXInversion;//!
+ map< int, vector<int> > m_CATSYInversion;//!
- // 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() {EventData->Clear();}
-
- // Give and external TMustData object to TMust2Physics. Needed for online analysis for example.
- void SetRawDataPointer(TCATSData* rawDataPointer) {EventData = rawDataPointer;}
-
-
- private :
+ public: // Output data of interest
+ // for a CATS
+ void SetTargetAngle(double m_TargetAngle) {m_TargetAngle = m_TargetAngle;}
+ void SetTargetThickness(double m_TargetThickness) {m_TargetThickness = m_TargetThickness;}
+
- // redundant information : could be optimized in the future
- vector< vector< vector<double> > > StripPositionX;//!
- vector< vector< vector<double> > > StripPositionY;//!
- vector<double> StripPositionZ;//!
- int NumberOfCATS;//!
+ // Remove bad channel, calibrate the data and apply threshold
+ void PreTreat();
- vector< vector <double> > Pedestal_X;//!
- vector< vector <double> > Pedestal_Y;//!
-
- vector< vector <double> > Threshold_X;//!
- vector< vector <double> > Threshold_Y;//!
+ // 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() ;
+
+ // Clear The PreTeated object
+ void ClearPreTreatedData() {m_PreTreatedData->Clear();}
- public : // Specific to CATS
+ void BuildPhysicalEvent();
- void Clear();
- void Clear(const Option_t*) {};
- void Dump() const;
-
- void AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28);
+ void BuildSimplePhysicalEvent();
- void ReadPedestal(string PedestalPath);
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
- double AnalyseX(int ff,int NumberOfDetector);
+ void Clear();
+ void Clear(const Option_t*) {};
+
+ // Give and external TCATSData object to TCATSPhysics, needed for online analysis
+ void SetRawDataPointer(TCATSData* rawDataPointer) {m_EventData = rawDataPointer;}
- double AnalyseY(int ff,int NumberOfDetector);
+ // Return false if the channel is disabled by user
+ bool IsValidChannel(const string DetectorType, const int Detector , const int channel);
- double CalculatePositionX( double CalculatedStripX, correction method);
+ void InitializeStandardParameter();
- double CalculatePositionY( double CalculatedStripY, correction method);
+ void AddParameterToCalibrationManager();
+ void ReadAnalysisConfig();
- reconstruction ChooseReconstruction(TString type);
+ void ReadConfiguration(string);
- // 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
+ void AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28);
- // 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 AnalyseX(int ff);
- double GetPositionOnTargetX() {return PositionOnTargetX;}
+ double AnalyseY(int ff);
- // Y
+ double CalculatePositionX( double CalculatedStripX, correction method);
- // 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);}
+ double CalculatePositionY( double CalculatedStripY, correction method);
- 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;}
+ reconstruction ChooseReconstruction(int ff, TString type);
+
+ reconstruction StringToEnum(string type);
- int GetCATSMult() {return PositionX.size();}
+ double CorrectedPositionX3(double Position, double a) ;
+ double CorrectedPositionY3(double Position, double a) ;
+ double CorrectedPositionX4(double Position, double b);
+ double CorrectedPositionY4(double Position, double b);
+ double Corrected3PointsX(double Position, double c);
+ double Corrected3PointsY(double Position, double c);
+ double Corrected4PointsX(double Position, double d);
+ double Corrected4PointsY(double Position, double d);
- TVector3 GetPositionOnTarget();
- TVector3 GetBeamDirection() {return BeamDirection;}
-
- ClassDef(TCATSPhysics,1) // CATSPhysics structure
- };
+
+ // Methode de reconstruction X
+ double HyperbolicSequentMethodX();
+ double GaussianMethodX();
+ double Barycentric5MethodX();
+ double Barycentric4MethodX();
+ double Barycentric3MethodX();
+ // Methode de Reconstruction Y
+ double HyperbolicSequentMethodY();
+ double GaussianMethodY();
+ double Barycentric5MethodY();
+ double Barycentric4MethodY();
+ double Barycentric3MethodY();
+ TVector3 GetBeamDirection();
+ TVector3 GetPositionOnTarget();
+
+ double GetPositionOnTargetX() {return PositionOnTargetX;}
+ double GetPositionOnTargetY() {return PositionOnTargetY;}
+ ClassDef(TCATSPhysics,1) // CATSPhysics structure
+};
namespace LOCAL_CATS
@@ -265,6 +237,12 @@ namespace LOCAL_CATS
// tranform an integer to a string
string itoa(int value);
+ double fCATS_X_Q(const TCATSData* Data, const int i);
+ double fCATS_Y_Q(const TCATSData* Data, const int i);
+ bool fCATS_Threshold_X(const TCATSData* Data, const int i);
+ bool fCATS_Threshold_Y(const TCATSData* Data, const int i);
+ double fCATS_Ped_X(const TCATSData* m_EventData, const int i);
+ double fCATS_Ped_Y(const TCATSData* m_EventData, const int i);
}
#endif
diff --git a/NPLib/Tools/CalibrationManager.cxx b/NPLib/Tools/CalibrationManager.cxx
index 68e399c9c2d13363acc06badecda23b61e2eb7d9..ad30e018ad29fa9531e73e855e1ce6b652599d1f 100644
--- a/NPLib/Tools/CalibrationManager.cxx
+++ b/NPLib/Tools/CalibrationManager.cxx
@@ -9,7 +9,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : october 2009 *
- * Last update : *
+ * Last update : december 2012 *
*---------------------------------------------------------------------------*
* Decription: *
* This class manage the calibration coefficient of the detector in NPA *
@@ -237,6 +237,73 @@ double CalibrationManager::ApplyCalibration(const string& ParameterPath , const
return CalibratedValue ;
}
+
//////////////////////////////////////////////////////////////////
+bool CalibrationManager::ApplyThreshold(const string& ParameterPath, const double& RawValue)
+{
+ map< string , vector<double> >::iterator it ;
+
+ // Find the good parameter in the Map
+ // Using Find method of stl is the fastest way
+ it = fCalibrationCoeff.find(ParameterPath) ;
+
+ // If the find methods return the end iterator it's mean the parameter was not found
+ if(it == fCalibrationCoeff.end() )
+ {
+ // cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
+ // cout << " ERROR: PARAMETER " << ParameterPath << " IS NOT FOUND IN THE CALIBRATION DATA BASE " << endl ;
+ // cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
+ return false;
+ }
+
+ // Else we take the second part of the element (first is index, ie: parameter path)
+ // Second is the vector of Coeff
+ vector<double> Coeff = it->second ;
+
+ // The vector size give the degree of calibration
+ // We just apply the coeff and returned the calibrated value
+
+ double ThresholdValue = 0 ;
+
+ if(Coeff.size()==2){ThresholdValue = Coeff[0] + 3*Coeff[1];}
+
+ if(RawValue > ThresholdValue)
+ {
+ return true;
+ }
+ else return false;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////
+double CalibrationManager::GetPedestal(const string& ParameterPath)
+{
+ map< string , vector<double> >::iterator it ;
+
+ // Find the good parameter in the Map
+ // Using Find method of stl is the fastest way
+ it = fCalibrationCoeff.find(ParameterPath) ;
+
+ // If the find methods return the end iterator it's mean the parameter was not found
+ if(it == fCalibrationCoeff.end() )
+ {
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
+ cout << " ERROR: PARAMETER " << ParameterPath << " IS NOT FOUND IN THE CALIBRATION DATA BASE " << endl ;
+ cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
+ }
+
+ // Else we take the second part of the element (first is index, ie: parameter path)
+ // Second is the vector of Coeff
+ vector<double> Coeff = it->second ;
+
+ // The vector size give the degree of calibration
+ // We just apply the coeff and returned the calibrated value
+
+ double PedestalValue = 0 ;
+
+ if(Coeff.size()==2){PedestalValue = Coeff[0];}
+
+
+ return PedestalValue;
+}
diff --git a/NPLib/Tools/CalibrationManager.h b/NPLib/Tools/CalibrationManager.h
index 2de05b5fc77a25ad14cf998ed52570025ee86137..0c620e587049bc3c729709fc10d33d76985f7376 100644
--- a/NPLib/Tools/CalibrationManager.h
+++ b/NPLib/Tools/CalibrationManager.h
@@ -11,7 +11,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : october 2009 *
- * Last update : *
+ * Last update : december 2012 *
*---------------------------------------------------------------------------*
* Decription: *
* This class manage the calibration coefficient of the detector in NPA *
@@ -59,6 +59,11 @@ class CalibrationManager
// call like : myCalibrationManager->ApplyCalibration( "MUST2/Telescope5_Si_X38_E" , RawEnergy )
// return the Calibrated value
double ApplyCalibration(const string& ParameterPath , const double& RawValue);
+
+ bool ApplyThreshold(const string& ParameterPath, const double& RawValue);
+
+ double GetPedestal(const string& ParameterPath);
+
public: // To be called after initialisation
@@ -68,7 +73,7 @@ class CalibrationManager
public: //Clear calibration if we have some calibration files to read
void ClearCalibration();
- public: //Get calibration coefficient vector
+ public: //Get correction coefficient vector
vector<double> GetCorrection(const string& ParameterPath);