From eb8d3857be126fbddce0a1f7bae23b7db5e5624e Mon Sep 17 00:00:00 2001 From: adrien-matta <a.matta@surrey.ac.uk> Date: Thu, 19 Dec 2013 16:44:19 +0100 Subject: [PATCH] * E-Theta for TiaraBarrel full detector --- NPLib/CATS/TCATSPhysics.cxx | 3361 ++++++++++++++------------- NPLib/CATS/TCATSPhysics.h | 544 ++--- NPLib/Tiara/TTiaraBarrelSpectra.cxx | 9 +- 3 files changed, 1966 insertions(+), 1948 deletions(-) diff --git a/NPLib/CATS/TCATSPhysics.cxx b/NPLib/CATS/TCATSPhysics.cxx index cd2a93e63..08f30e86a 100644 --- a/NPLib/CATS/TCATSPhysics.cxx +++ b/NPLib/CATS/TCATSPhysics.cxx @@ -1,1677 +1,1684 @@ -/***************************************************************************** - * Copyright (C) 2009-2013 this file is part of the NPTool Project * - * * - * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * - * For the list of contributors see $NPTOOL/Licence/Contributors * - *****************************************************************************/ - -/***************************************************************************** - * Original Author: Sandra Giron contact address: giron@ipno.in2p3.fr * - * * - * Creation Date : febuary 2010 * - * Last update : modification november 2011 by Pierre Morfouace * - * Contact adress : morfouac@ipno.in2p3.fr * - *---------------------------------------------------------------------------* - * Decription: * - * This class hold CATS treated data * - * * - *---------------------------------------------------------------------------* - * Comment: * - * * - *****************************************************************************/ - -#include "TCATSPhysics.h" -using namespace LOCAL_CATS; - -// STL -#include <cmath> -#include <algorithm> -#include <sstream> -#include <fstream> -#include <iostream> -#include <stdlib.h> -// NPL -#include "RootInput.h" -#include "RootOutput.h" -// ROOT -#include "TChain.h" -#include "TRandom.h" - -ClassImp(TCATSPhysics) - - /////////////////////////////////////////////////////////////////////////// -TCATSPhysics::TCATSPhysics() -{ - m_EventData = new TCATSData ; - m_PreTreatedData = new TCATSData ; - m_EventPhysics = this ; - m_Spectra = 0; - 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; -} - -///////////////////////////////////////////////////////////////////////////// -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; - } - } - - 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) -{ - ifstream ConfigFile; - ConfigFile.open(Path.c_str()); - string LineBuffer ; - string DataBuffer ; - - double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ; - TVector3 A , B , C , D ; - - bool check_A = false ; - bool check_B = false ; - bool check_C = false ; - bool check_D = false ; - - bool ReadingStatus = false ; - - - while (!ConfigFile.eof()) - { - getline(ConfigFile, LineBuffer); - - //If line is a Start Up CATS bloc, Reading toggle to true - if (LineBuffer.compare(0, 12, "CATSDetector") == 0) - { - cout << "CATS Detector found: " << endl ; - ReadingStatus = true ; - } - - // Else don't toggle to Reading Block Status - else ReadingStatus = false ; - - // Reading Block - while(ReadingStatus) - { - ConfigFile >> DataBuffer ; - // Comment Line - if(DataBuffer.compare(0, 1, "%") == 0) { - ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' ); - } - - // Finding another telescope (safety), toggle out - else if (DataBuffer.compare(0, 12, "CATSDetector") == 0) { - cout << "Warning: Another CATS is found before standard sequence of Token, Error may have occurred in CATS definition" << endl ; - ReadingStatus = false ; - } - - // Corner Position method - - else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) { - check_A = true; - ConfigFile >> DataBuffer ; - Ax = atof(DataBuffer.c_str()) ; - Ax = Ax ; - ConfigFile >> DataBuffer ; - Ay = atof(DataBuffer.c_str()) ; - Ay = Ay ; - ConfigFile >> DataBuffer ; - Az = atof(DataBuffer.c_str()) ; - Az = Az ; - - A = TVector3(Ax, Ay, Az); - cout << " X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl; - } - - else if (DataBuffer.compare(0, 7, "X28_Y1=") == 0) { - check_B = true; - ConfigFile >> DataBuffer ; - Bx = atof(DataBuffer.c_str()) ; - Bx = Bx ; - ConfigFile >> DataBuffer ; - By = atof(DataBuffer.c_str()) ; - By = By ; - ConfigFile >> DataBuffer ; - Bz = atof(DataBuffer.c_str()) ; - Bz = Bz ; - - B = TVector3(Bx, By, Bz); - cout << " X28 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl; - } - - else if (DataBuffer.compare(0, 7, "X1_Y28=") == 0) { - check_C = true; - ConfigFile >> DataBuffer ; - Cx = atof(DataBuffer.c_str()) ; - Cx = Cx ; - ConfigFile >> DataBuffer ; - Cy = atof(DataBuffer.c_str()) ; - Cy = Cy ; - ConfigFile >> DataBuffer ; - Cz = atof(DataBuffer.c_str()) ; - Cz = Cz ; - - C = TVector3(Cx, Cy, Cz); - cout << " X1 Y28 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl; - } - - else if (DataBuffer.compare(0, 8, "X28_Y28=") == 0) { - check_D = true; - ConfigFile >> DataBuffer ; - Dx = atof(DataBuffer.c_str()) ; - Dx = Dx ; - ConfigFile >> DataBuffer ; - Dy = atof(DataBuffer.c_str()) ; - Dy = Dy ; - ConfigFile >> DataBuffer ; - Dz = atof(DataBuffer.c_str()) ; - Dz = Dz ; - - D = TVector3(Dx, Dy, Dz); - cout << " X28 Y28 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl; - - } - - // End Corner Position Method - - ///////////////////////////////////////////////// - // If All necessary information there, toggle out - if (check_A && check_B && check_C && check_D) - { - ReadingStatus = false; - - ///Add The previously define telescope - - AddCATS( A , - B , - C , - D ); - - check_A = false; - check_B = false; - check_C = false; - check_D = false; - } - } - - } - 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); -} - -///////////////////////////////////////////////////////////////////// -// 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() -{ - TChain* inputChain = RootInput::getInstance()->GetChain() ; - inputChain->SetBranchStatus( "CATS" , true ) ; - inputChain->SetBranchStatus( "fCATS_*" , true ) ; - 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( "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" , &m_EventPhysics ) ; -} - -///////////////////////////////////////////////////////////////////// -void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28) -{ - m_NumberOfCATS++ ; - - // remove warning - C_X28_Y28 *= 1; - - // Vector U on Telescope Face (paralelle to Y Strip) - TVector3 U = C_X28_Y1 - C_X1_Y1 ; - U = U.Unit() ; - - // Vector V on Telescope Face (parallele to X Strip) - TVector3 V = C_X1_Y28 - C_X1_Y1 ; - V = V.Unit() ; - - // Position Vector of Strip Center - TVector3 StripCenter ; - // Position Vector of X=1 Y=1 Strip - TVector3 Strip_1_1 ; - - // Geometry Parameters - double Face = 71.12 ; //mm - double NumberOfStrip = 28 ; - double StripPitch = Face / NumberOfStrip ; //mm - - // Buffer object to fill Position Array - vector<double> lineX ; - vector<double> lineY ; - - vector< vector< double > > OneDetectorStripPositionX ; - vector< vector< double > > OneDetectorStripPositionY ; - double OneDetectorStripPositionZ ; - - // Moving StripCenter to 1.1 corner (strip center!) : - Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2) ; - - - for( int i = 0 ; i < 28 ; i++ ) - { - lineX.clear() ; - lineY.clear() ; - - for( int j = 0 ; j < 28 ; j++ ) - { - StripCenter = Strip_1_1 + StripPitch*( i*U + j*V ) ; - lineX.push_back( StripCenter.x() ) ; - lineY.push_back( StripCenter.y() ) ; - } - - OneDetectorStripPositionX.push_back(lineX); - OneDetectorStripPositionY.push_back(lineY); - } - - OneDetectorStripPositionZ = C_X1_Y1.Z(); - - StripPositionX.push_back(OneDetectorStripPositionX) ; - StripPositionY.push_back(OneDetectorStripPositionY) ; - StripPositionZ.push_back(OneDetectorStripPositionZ) ; - -} - -/////////////////////////////////////////////////////////////// -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() ; - - - ff = 0; - HitX = 0; - HitY = 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::ReadAnalysisConfig() -{ - bool ReadingStatus = false; - - // 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;} - - } - } -} - - -/////////////////////////////////////////////////////////////////////////// -void TCATSPhysics::InitSpectra(){ - m_Spectra = new TCATSSpectra(m_NumberOfCATS); -} - -/////////////////////////////////////////////////////////////////////////// -void TCATSPhysics::FillSpectra(){ - m_Spectra -> FillRawSpectra(m_EventData); - m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); - m_Spectra -> FillPhysicsSpectra(m_EventPhysics); -} -/////////////////////////////////////////////////////////////////////////// -void TCATSPhysics::CheckSpectra(){ - // To be done -} -/////////////////////////////////////////////////////////////////////////// -void TCATSPhysics::ClearSpectra(){ - // To be done -} -/////////////////////////////////////////////////////////////////////////// -map< vector<TString> , TH1*> TCATSPhysics::GetSpectra() { - return m_Spectra->GetMapHisto(); -} - -///////////////////////////////////////////////////////////////////// -// Add Parameter to the CalibrationManger -void TCATSPhysics::AddParameterToCalibrationManager() -{ - CalibrationManager* Cal = CalibrationManager::getInstance(); - for(int i = 0 ; i < m_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") ; - 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)) ; - } - } - - return; -} - -//////////////////////////////////////////////////////////////// -double TCATSPhysics::AnalyseX(int ff) -{ - double CalculatedStripX=0; - - ReconstructionMethodX[ff] = ChooseReconstruction(ff,"X"); - - if(ReconstructionMethodX[ff] == SECHS)CalculatedStripX = HyperbolicSequentMethodX(); - else if(ReconstructionMethodX[ff] == GAUSS)CalculatedStripX = GaussianMethodX(); - else if(ReconstructionMethodX[ff] == BAR3) CalculatedStripX = Barycentric3MethodX(); - else if(ReconstructionMethodX[ff] == BAR4) CalculatedStripX = Barycentric4MethodX(); - else if(ReconstructionMethodX[ff] == BAR5) CalculatedStripX = Barycentric5MethodX(); - else CalculatedStripX = Barycentric3MethodX(); - - return(CalculatedStripX); -} - -//////////////////////////////////////////////////////////////// -double TCATSPhysics::AnalyseY(int ff) -{ - double CalculatedStripY=0; - - ReconstructionMethodY[ff] = ChooseReconstruction(ff,"Y"); - - if(ReconstructionMethodY[ff] == SECHS)CalculatedStripY = HyperbolicSequentMethodY(); - else if(ReconstructionMethodY[ff] == GAUSS)CalculatedStripY = GaussianMethodY(); - else if(ReconstructionMethodY[ff] == BAR3) CalculatedStripY = Barycentric3MethodY(); - else if(ReconstructionMethodY[ff] == BAR4) CalculatedStripY = Barycentric4MethodY(); - else if(ReconstructionMethodY[ff] == BAR5) CalculatedStripY = Barycentric5MethodY(); - else CalculatedStripY = Barycentric3MethodY(); - - return(CalculatedStripY); -} - -//////////////////////////////////////////////////////////////////////// -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;} - } - - 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); -} - -////////////////////////////////////////////////////////////////////// -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{ - IStripX = (int) CalculatedStripX ; - - // Decimal Part - double DStripX = CalculatedStripX-IStripX ; - - 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; - } - - 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; -} - -///////////////////////////////////////////////////////////////////////// -double TCATSPhysics::CalculatePositionY(double CalculatedStripY, 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 ){ - 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); - } - } - - 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); - } - } - - else cout << "only 2CATS!! ff = " << ff << endl; - } - - else {positionY = -1000;} - } - - - 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); - } - } - - 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); - } - } - - return positionY; - -} - -//////////////////////////////////////////////////////////////////// -double TCATSPhysics:: GaussianMethodX() -{ - 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; - } - - - if(StripMaxX[ff]> 3 && StripMaxX[ff]< 26) - { - Q[0] = Buffer_X_Q[StripMax-1][ff] ; - StripPos[0] = 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]; - } - - 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(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]) ; - if(denom != 0){ - gauss = 0.5*num / denom; - } - else{gauss = -1000;} - - } - - else { - gauss = -1000; - } - - 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::Corrected4PointsX(double Position, double d) -{ - 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 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; - } - - int StripMax = StripMaxY[ff]; - - double gauss = -1000; - - if(StripMaxY[ff] > 2 && StripMaxY[ff]<27) - { - Q[0] = Buffer_Y_Q[StripMax-1][ff] ; - StripPos[0] = StripPositionY[ff][0][StripMax-1]; - - 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] ; - } - } - } - - - - Q0_Q1 = log(Q[0]/Q[1]); - Q0_Q2 = log(Q[0]/Q[2]); - - num = Q0_Q1 * (StripPos[2]*StripPos[2] - StripPos[0]*StripPos[0]) - Q0_Q2 * (StripPos[1]*StripPos[1] - StripPos[0]*StripPos[0]) ; - denom = Q0_Q1 * (StripPos[2] - StripPos[0]) - Q0_Q2 * (StripPos[1] - StripPos[0]) ; - if(denom != 0){ - gauss = 0.5*num / denom; - } - - } - - else { - gauss = -1000; - } - - 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 ; - double NumberOfPoint = 0 ; - double ChargeTotal =0; - - - for(int i = -2 ; i < 3 ; i++) - { - if(Buffer_X_Q[StripMax_+i][ff]!=-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 ; } - - } - else - { - Barycenter = -1000 ; - } - - return Barycenter ; -} - -/////////////////////////////////////////////////////////////// -double TCATSPhysics:: Barycentric5MethodY() -{ - double Barycenter = 0 ; - - if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) - { - int StripMax_ = StripMaxY[ff] -1 ; - double NumberOfPoint = 0 ; - double ChargeTotal =0; - - - for(int i = -2 ; i < 3 ; i++) - { - if(Buffer_Y_Q[StripMax_+i][ff]!=-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 -{ - Barycenter = -1000 ; -} - -return Barycenter ; -} - -/////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////// -double TCATSPhysics:: Barycentric3MethodX() -{ - double Barycenter = 0 ; - - if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27) - { - 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 -{ - Barycenter = -1000 ; -} - -return Barycenter ; -} - -/////////////////////////////////////////////////////////////// -double TCATSPhysics:: Barycentric3MethodY() -{ - double Barycenter = 0 ; - - if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) - { - int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ; - double NumberOfPoint = 0 ; - double ChargeTotal =0; - - for(int i = -1 ; i < 2 ; i++) - { - if(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 -{ - Barycenter = -1000 ; - // cout << "Strip max " << StripMax << endl; -} - -return Barycenter ; -} - -/////////////////////////////////////////////////////////////// -double TCATSPhysics:: Barycentric4MethodX() -{ - double Barycenter = 0 ; - - if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27) { - - int StripMax_ = StripMaxX[ff] -1 ; // Use because numerotation of array start at 0 ; - double NumberOfPoint = 0 ; - double ChargeTotal =0; - - if(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 ; - } - - } - - else - { - Barycenter = -1000 ; - // cout << "Strip max " << StripMax << endl; - } - - return Barycenter ; -} - -/////////////////////////////////////////////////////////////// -double TCATSPhysics:: Barycentric4MethodY() -{ - double Barycenter = 0 ; - - if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) { - - int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ; - double NumberOfPoint = 0 ; - double ChargeTotal =0; - - if(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 ; -} - -///////////////////////////////////////////////////////////////////// -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 ; -} - -////////////////////////////////////////////////////////////////// -double TCATSPhysics:: HyperbolicSequentMethodY() -{ - double sechs = -1000 ; - - if(StripMaxY[ff] > 2 && StripMaxY[ff]<27) - { - 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] ; - - 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_]; - - 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(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) / b + xmax; - - return Corrected_Position; -} - -/////////////////////////////////////////////////////////////// -double TCATSPhysics::CorrectedPositionY4(double Position, double b) -{ - double Corrected_Position = 0; - double xmax = 0; - int StripMax_ = StripMaxY[ff] -1; - - if(Buffer_Y_Q[StripMax_+1][ff] > Buffer_Y_Q[StripMax_-1][ff]) { - xmax = StripPositionY[ff][0][StripMax_] + 1.27 ; - } - - else{ - xmax = StripPositionY[ff][0][StripMax_-1] + 1.27; - } - - Corrected_Position = (Position - xmax) / b + xmax; - - return Corrected_Position; -} - - -/////////////////////////////////////////////////////////////// -TVector3 TCATSPhysics::GetBeamDirection() -{ - TVector3 Position = TVector3 (PositionX[1]-PositionX[0] , - PositionY[1]-PositionY[0] , - PositionZ[1]-PositionZ[0] ); - Position.Unit(); - return(Position) ; -} - -/////////////////////////////////////////////////////////////// -TVector3 TCATSPhysics::GetPositionOnTarget() -{ - double Pi = 3.14159265; - TVector3 Position = TVector3 (GetPositionOnTargetX() , - GetPositionOnTargetY() , - GetPositionOnTargetX()*tan(m_TargetAngle*Pi/180)); - Position.Unit(); - return(Position) ; -} - -//////////////////////////////////////////////////////////////////////// -namespace LOCAL_CATS -{ - // 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) ); - } - - 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)); - } - - 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) ) ); - } - - -} - - - - - - - - - - +/***************************************************************************** + * Copyright (C) 2009-2010 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Sandra Giron contact address: giron@ipno.in2p3.fr * + * * + * Creation Date : febuary 2010 * + * Last update : modification november 2011 by Pierre Morfouace * + * Contact adress : morfouac@ipno.in2p3.fr * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold CATS treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#include "TCATSPhysics.h" +using namespace LOCAL_CATS; + +// STL +#include <cmath> +#include <algorithm> +#include <sstream> +#include <fstream> +#include <iostream> +#include <stdlib.h> +// NPL +#include "RootInput.h" +#include "RootOutput.h" +// ROOT +#include "TChain.h" +#include "TRandom.h" + +ClassImp(TCATSPhysics) + + /////////////////////////////////////////////////////////////////////////// +TCATSPhysics::TCATSPhysics() +{ + 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; +} + +///////////////////////////////////////////////////////////////////////////// +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++; + DetectorID = m_PreTreatedData->GetCATSDetX(i); + } + 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) +{ + ifstream ConfigFile; + ConfigFile.open(Path.c_str()); + string LineBuffer ; + string DataBuffer ; + + double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ; + TVector3 A , B , C , D ; + + bool check_A = false ; + bool check_B = false ; + bool check_C = false ; + bool check_D = false ; + + bool ReadingStatus = false ; + + + while (!ConfigFile.eof()) + { + getline(ConfigFile, LineBuffer); + + //If line is a Start Up CATS bloc, Reading toggle to true + if (LineBuffer.compare(0, 12, "CATSDetector") == 0) + { + cout << "CATS Detector found: " << endl ; + ReadingStatus = true ; + } + + // Else don't toggle to Reading Block Status + else ReadingStatus = false ; + + // Reading Block + while(ReadingStatus) + { + ConfigFile >> DataBuffer ; + // Comment Line + if(DataBuffer.compare(0, 1, "%") == 0) { + ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' ); + } + + // Finding another telescope (safety), toggle out + else if (DataBuffer.compare(0, 12, "CATSDetector") == 0) { + cout << "Warning: Another CATS is found before standard sequence of Token, Error may have occurred in CATS definition" << endl ; + ReadingStatus = false ; + } + + // Corner Position method + + else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) { + check_A = true; + ConfigFile >> DataBuffer ; + Ax = atof(DataBuffer.c_str()) ; + Ax = Ax ; + ConfigFile >> DataBuffer ; + Ay = atof(DataBuffer.c_str()) ; + Ay = Ay ; + ConfigFile >> DataBuffer ; + Az = atof(DataBuffer.c_str()) ; + Az = Az ; + + A = TVector3(Ax, Ay, Az); + cout << " X1 Y1 corner position : (" << A.X() << ";" << A.Y() << ";" << A.Z() << ")" << endl; + } + + else if (DataBuffer.compare(0, 7, "X28_Y1=") == 0) { + check_B = true; + ConfigFile >> DataBuffer ; + Bx = atof(DataBuffer.c_str()) ; + Bx = Bx ; + ConfigFile >> DataBuffer ; + By = atof(DataBuffer.c_str()) ; + By = By ; + ConfigFile >> DataBuffer ; + Bz = atof(DataBuffer.c_str()) ; + Bz = Bz ; + + B = TVector3(Bx, By, Bz); + cout << " X28 Y1 corner position : (" << B.X() << ";" << B.Y() << ";" << B.Z() << ")" << endl; + } + + else if (DataBuffer.compare(0, 7, "X1_Y28=") == 0) { + check_C = true; + ConfigFile >> DataBuffer ; + Cx = atof(DataBuffer.c_str()) ; + Cx = Cx ; + ConfigFile >> DataBuffer ; + Cy = atof(DataBuffer.c_str()) ; + Cy = Cy ; + ConfigFile >> DataBuffer ; + Cz = atof(DataBuffer.c_str()) ; + Cz = Cz ; + + C = TVector3(Cx, Cy, Cz); + cout << " X1 Y28 corner position : (" << C.X() << ";" << C.Y() << ";" << C.Z() << ")" << endl; + } + + else if (DataBuffer.compare(0, 8, "X28_Y28=") == 0) { + check_D = true; + ConfigFile >> DataBuffer ; + Dx = atof(DataBuffer.c_str()) ; + Dx = Dx ; + ConfigFile >> DataBuffer ; + Dy = atof(DataBuffer.c_str()) ; + Dy = Dy ; + ConfigFile >> DataBuffer ; + Dz = atof(DataBuffer.c_str()) ; + Dz = Dz ; + + D = TVector3(Dx, Dy, Dz); + cout << " X28 Y28 corner position : (" << D.X() << ";" << D.Y() << ";" << D.Z() << ")" << endl; + + } + + // End Corner Position Method + + ///////////////////////////////////////////////// + // If All necessary information there, toggle out + if (check_A && check_B && check_C && check_D) + { + ReadingStatus = false; + + ///Add The previously define telescope + + AddCATS( A , + B , + C , + D ); + + check_A = false; + check_B = false; + check_C = false; + check_D = false; + } + } + + } + 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); +} + +///////////////////////////////////////////////////////////////////// +// 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() +{ + TChain* inputChain = RootInput::getInstance()->GetChain() ; + inputChain->SetBranchStatus( "CATS" , true ) ; + inputChain->SetBranchStatus( "fCATS_*" , true ) ; + 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( "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" , &m_EventPhysics ) ; +} + +///////////////////////////////////////////////////////////////////// +void TCATSPhysics::AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28) +{ + m_NumberOfCATS++ ; + + // remove warning + C_X28_Y28 *= 1; + + // Vector U on Telescope Face (paralelle to Y Strip) + TVector3 U = C_X28_Y1 - C_X1_Y1 ; + U = U.Unit() ; + + // Vector V on Telescope Face (parallele to X Strip) + TVector3 V = C_X1_Y28 - C_X1_Y1 ; + V = V.Unit() ; + + // Position Vector of Strip Center + TVector3 StripCenter ; + // Position Vector of X=1 Y=1 Strip + TVector3 Strip_1_1 ; + + // Geometry Parameters + double Face = 71.12 ; //mm + double NumberOfStrip = 28 ; + double StripPitch = Face / NumberOfStrip ; //mm + + // Buffer object to fill Position Array + vector<double> lineX ; + vector<double> lineY ; + + vector< vector< double > > OneDetectorStripPositionX ; + vector< vector< double > > OneDetectorStripPositionY ; + double OneDetectorStripPositionZ ; + + // Moving StripCenter to 1.1 corner (strip center!) : + Strip_1_1 = C_X1_Y1 + (U+V) * (StripPitch/2) ; + + + for( int i = 0 ; i < 28 ; i++ ) + { + lineX.clear() ; + lineY.clear() ; + + for( int j = 0 ; j < 28 ; j++ ) + { + StripCenter = Strip_1_1 + StripPitch*( i*U + j*V ) ; + lineX.push_back( StripCenter.x() ) ; + lineY.push_back( StripCenter.y() ) ; + } + + OneDetectorStripPositionX.push_back(lineX); + OneDetectorStripPositionY.push_back(lineY); + } + + OneDetectorStripPositionZ = C_X1_Y1.Z(); + + StripPositionX.push_back(OneDetectorStripPositionX) ; + StripPositionY.push_back(OneDetectorStripPositionY) ; + StripPositionZ.push_back(OneDetectorStripPositionZ) ; + +} + +/////////////////////////////////////////////////////////////// +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(); + + ff = 0; + HitX = 0; + HitY = 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::ReadAnalysisConfig() +{ + bool ReadingStatus = false; + + // 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;} + + } + } +} +/////////////////////////////////////////////////////////////////////////// +void TCATSPhysics::InitSpectra(){ + m_Spectra = new TCATSSpectra(m_NumberOfCATS); +} + +/////////////////////////////////////////////////////////////////////////// +void TCATSPhysics::FillSpectra(){ + m_Spectra -> FillRawSpectra(m_EventData); + m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData); + m_Spectra -> FillPhysicsSpectra(m_EventPhysics); +} +/////////////////////////////////////////////////////////////////////////// +void TCATSPhysics::CheckSpectra(){ + // To be done +} +/////////////////////////////////////////////////////////////////////////// +void TCATSPhysics::ClearSpectra(){ + // To be done +} +/////////////////////////////////////////////////////////////////////////// +map< vector<TString> , TH1*> TCATSPhysics::GetSpectra() { + return m_Spectra->GetMapHisto(); +} + +///////////////////////////////////////////////////////////////////// +// Add Parameter to the CalibrationManger +void TCATSPhysics::AddParameterToCalibrationManager() +{ + CalibrationManager* Cal = CalibrationManager::getInstance(); + for(int i = 0 ; i < m_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") ; + 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)) ; + } + } + + return; +} + +//////////////////////////////////////////////////////////////// +double TCATSPhysics::AnalyseX(int ff) +{ + double CalculatedStripX=0; + + ReconstructionMethodX[ff] = ChooseReconstruction(ff,"X"); + + 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); +} + +//////////////////////////////////////////////////////////////// +double TCATSPhysics::AnalyseY(int ff) +{ + double CalculatedStripY=0; + + ReconstructionMethodY[ff] = ChooseReconstruction(ff,"Y"); + + 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); +} + +//////////////////////////////////////////////////////////////////////// +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;} + } + + 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); +} + +////////////////////////////////////////////////////////////////////// +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{ + IStripX = (int) CalculatedStripX ; + + // Decimal Part + double DStripX = CalculatedStripX-IStripX ; + + 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; + } + + 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; +} + +///////////////////////////////////////////////////////////////////////// +double TCATSPhysics::CalculatePositionY(double CalculatedStripY, 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 ){ + 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); + } + } + + 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); + } + } + + else cout << "only 2CATS!! ff = " << ff << endl; + } + + else {positionY = -1000;} + } + + + 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); + } + } + + 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); + } + } + + return positionY; + +} + +//////////////////////////////////////////////////////////////////// +double TCATSPhysics:: GaussianMethodX() +{ + 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; + } + + + if(StripMaxX[ff]> 3 && StripMaxX[ff]< 26) + { + Q[0] = Buffer_X_Q[StripMax-1][ff] ; + StripPos[0] = 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]; + } + + 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(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]) ; + if(denom != 0){ + gauss = 0.5*num / denom; + } + else{gauss = -1000;} + + } + + else { + gauss = -1000; + } + + 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::Corrected4PointsX(double Position, double d) +{ + 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 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; + } + + int StripMax = StripMaxY[ff]; + + double gauss = -1000; + + if(StripMaxY[ff] > 2 && StripMaxY[ff]<27) + { + Q[0] = Buffer_Y_Q[StripMax-1][ff] ; + StripPos[0] = StripPositionY[ff][0][StripMax-1]; + + 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] ; + } + } + } + + + + Q0_Q1 = log(Q[0]/Q[1]); + Q0_Q2 = log(Q[0]/Q[2]); + + num = Q0_Q1 * (StripPos[2]*StripPos[2] - StripPos[0]*StripPos[0]) - Q0_Q2 * (StripPos[1]*StripPos[1] - StripPos[0]*StripPos[0]) ; + denom = Q0_Q1 * (StripPos[2] - StripPos[0]) - Q0_Q2 * (StripPos[1] - StripPos[0]) ; + if(denom != 0){ + gauss = 0.5*num / denom; + } + + } + + else { + gauss = -1000; + } + + 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 ; + double NumberOfPoint = 0 ; + double ChargeTotal =0; + + + for(int i = -2 ; i < 3 ; i++) + { + if(Buffer_X_Q[StripMax_+i][ff]!=-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 ; } + + } + else + { + Barycenter = -1000 ; + } + + return Barycenter ; +} + +/////////////////////////////////////////////////////////////// +double TCATSPhysics:: Barycentric5MethodY() +{ + double Barycenter = 0 ; + + if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) + { + int StripMax_ = StripMaxY[ff] -1 ; + double NumberOfPoint = 0 ; + double ChargeTotal =0; + + + for(int i = -2 ; i < 3 ; i++) + { + if(Buffer_Y_Q[StripMax_+i][ff]!=-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 +{ + Barycenter = -1000 ; +} + +return Barycenter ; +} + +/////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////// +double TCATSPhysics:: Barycentric3MethodX() +{ + double Barycenter = 0 ; + + if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27) + { + 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 +{ + Barycenter = -1000 ; +} + +return Barycenter ; +} + +/////////////////////////////////////////////////////////////// +double TCATSPhysics:: Barycentric3MethodY() +{ + double Barycenter = 0 ; + + if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) + { + int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ; + double NumberOfPoint = 0 ; + double ChargeTotal =0; + + for(int i = -1 ; i < 2 ; i++) + { + if(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 +{ + Barycenter = -1000 ; + // cout << "Strip max " << StripMax << endl; +} + +return Barycenter ; +} + +/////////////////////////////////////////////////////////////// +double TCATSPhysics:: Barycentric4MethodX() +{ + double Barycenter = 0 ; + + if(StripMaxX[ff] > 2 && StripMaxX[ff] < 27) { + + int StripMax_ = StripMaxX[ff] -1 ; // Use because numerotation of array start at 0 ; + double NumberOfPoint = 0 ; + double ChargeTotal =0; + + if(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 ; + } + + } + + else + { + Barycenter = -1000 ; + // cout << "Strip max " << StripMax << endl; + } + + return Barycenter ; +} + +/////////////////////////////////////////////////////////////// +double TCATSPhysics:: Barycentric4MethodY() +{ + double Barycenter = 0 ; + + if(StripMaxY[ff] > 2 && StripMaxY[ff] < 27) { + + int StripMax_ = StripMaxY[ff] -1 ; // Use because numerotation of array start at 0 ; + double NumberOfPoint = 0 ; + double ChargeTotal =0; + + if(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 ; +} + +///////////////////////////////////////////////////////////////////// +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 ; +} + +////////////////////////////////////////////////////////////////// +double TCATSPhysics:: HyperbolicSequentMethodY() +{ + double sechs = -1000 ; + + if(StripMaxY[ff] > 2 && StripMaxY[ff]<27) + { + 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] ; + + 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_]; + + 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(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) / b + xmax; + + return Corrected_Position; +} + +/////////////////////////////////////////////////////////////// +double TCATSPhysics::CorrectedPositionY4(double Position, double b) +{ + double Corrected_Position = 0; + double xmax = 0; + int StripMax_ = StripMaxY[ff] -1; + + if(Buffer_Y_Q[StripMax_+1][ff] > Buffer_Y_Q[StripMax_-1][ff]) { + xmax = StripPositionY[ff][0][StripMax_] + 1.27 ; + } + + else{ + xmax = StripPositionY[ff][0][StripMax_-1] + 1.27; + } + + Corrected_Position = (Position - xmax) / b + xmax; + + return Corrected_Position; +} + + +/////////////////////////////////////////////////////////////// +TVector3 TCATSPhysics::GetBeamDirection() +{ + TVector3 Position = TVector3 (PositionX[1]-PositionX[0] , + PositionY[1]-PositionY[0] , + PositionZ[1]-PositionZ[0] ); + Position.Unit(); + return(Position) ; +} + +/////////////////////////////////////////////////////////////// +TVector3 TCATSPhysics::GetPositionOnTarget() +{ + double Pi = 3.14159265; + TVector3 Position = TVector3 (GetPositionOnTargetX() , + GetPositionOnTargetY() , + GetPositionOnTargetX()*tan(m_TargetAngle*Pi/180)); + Position.Unit(); + return(Position) ; +} + +//////////////////////////////////////////////////////////////////////// +namespace LOCAL_CATS +{ + // 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) ); + } + + 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)); + } + + 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) ) ); + } + + +} + + + + + + + + + + diff --git a/NPLib/CATS/TCATSPhysics.h b/NPLib/CATS/TCATSPhysics.h index a7394fcac..a94109f9b 100644 --- a/NPLib/CATS/TCATSPhysics.h +++ b/NPLib/CATS/TCATSPhysics.h @@ -1,270 +1,274 @@ -#ifndef __CATSCAL__ -#define __CATSCAL__ -/***************************************************************************** - * Copyright (C) 2009-2013 this file is part of the NPTool Project * - * * - * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * - * For the list of contributors see $NPTOOL/Licence/Contributors * - *****************************************************************************/ - -/***************************************************************************** - * Original Author: Sandra Giron contact address: giron@ipno.in2p3.fr * - * * - * Creation Date : febuary 2010 * - * Last update : * - *---------------------------------------------------------------------------* - * Decription: * - * This class hold CATS treated data * - * * - *---------------------------------------------------------------------------* - * Comment: * - * * - *****************************************************************************/ - -// STL -#include <vector> -// ROOT -#include "TObject.h" -#include "TVector3.h" -#include <TRandom1.h> -#include <TRandom2.h> -#include <TRandom3.h> -// NPLib -#include "TCATSData.h" -#include "TCATSSpectra.h" -#include "../include/VDetector.h" -#include "../include/CalibrationManager.h" -#include "../include/DetectorManager.h" - -#define NBDETECTOR 2 -#define NBSTRIPS 28 - -// forward declaration -class TCATSSpectra; - -using namespace std ; -enum reconstruction{NO,SECHS,GAUSS,BAR3,BAR4,BAR5}; -enum correction{NOcor,cor}; - -class TCATSPhysics : public TObject, public NPA::VDetector -{ - - public: // Constructor and Destructor - TCATSPhysics(); - ~TCATSPhysics(); - - 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 - vector<int> DetNumberX; - vector<int> StripX; - vector<double> ChargeX; - - // Vector of dim = number of CATS - vector<int> StripMaxX; - - - // Vector of dim = multiplicity of event on all detector - vector<int> DetNumberY; - vector<int> StripY; - vector<double> ChargeY; - - // Vector of dim = number of CATS - vector<int> StripMaxY; - - // 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 Buffer_X_Q[NBSTRIPS][NBDETECTOR];//! - double Buffer_Y_Q[NBSTRIPS][NBDETECTOR];//! - - int HitX; //! - int HitY; //! - - vector<reconstruction> ReconstructionMethodX; - vector<reconstruction> ReconstructionMethodY; - - - 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;//! - - public: // Output data of interest - // for a CATS - void SetTargetAngle(double TargetAngle) {m_TargetAngle = TargetAngle;} - void SetTargetThickness(double TargetThickness) {m_TargetThickness = TargetThickness;} - - - // Remove bad channel, calibrate the data and apply threshold - void PreTreat(); - - // 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();} - - void BuildPhysicalEvent(); - - void BuildSimplePhysicalEvent(); - - // Same as above but for online analysis - void BuildOnlinePhysicalEvent() {BuildSimplePhysicalEvent();}; - - // Those two method all to clear the Event Physics or Data - void ClearEventPhysics() {Clear();} - void ClearEventData() {m_EventData->Clear();} - - // Method related to the TSpectra classes, aimed at providing a framework for online applications - // Instantiate the Spectra class and the histogramm throught it - void InitSpectra(); - // Fill the spectra hold by the spectra class - void FillSpectra(); - // Used for Online mainly, perform check on the histo and for example change their color if issues are found - void CheckSpectra(); - // Used for Online only, clear all the spectra hold by the Spectra class - void ClearSpectra(); - - 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;} - - // Return false if the channel is disabled by user - bool IsValidChannel(const string DetectorType, const int Detector , const int channel); - - void InitializeStandardParameter(); - - void AddParameterToCalibrationManager(); - - void ReadAnalysisConfig(); - - void ReadConfiguration(string); - - void AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28); - - double AnalyseX(int ff); - - double AnalyseY(int ff); - - double CalculatePositionX( double CalculatedStripX, correction method); - - double CalculatePositionY( double CalculatedStripY, correction method); - - reconstruction ChooseReconstruction(int ff, TString type); - - reconstruction StringToEnum(string type); - - 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); - - - // 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;} - - private: // Spectra Class - TCATSSpectra* m_Spectra;//! - - public: // Spectra Getter - map< vector<TString> , TH1*> GetSpectra(); - - - ClassDef(TCATSPhysics,1) // CATSPhysics structure -}; - - -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 +#ifndef __CATSCAL__ +#define __CATSCAL__ +/***************************************************************************** + * Copyright (C) 2009-2010 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Sandra Giron contact address: giron@ipno.in2p3.fr * + * * + * Creation Date : febuary 2010 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold CATS treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// STL +#include <vector> +// ROOT +#include "TObject.h" +#include "TVector3.h" +#include <TRandom1.h> +#include <TRandom2.h> +#include <TRandom3.h> +// NPLib +#include "TCATSData.h" +#include "TCATSSpectra.h" + +#include "../include/VDetector.h" +#include "../include/CalibrationManager.h" +#include "../include/DetectorManager.h" + +#define NBDETECTOR 2 +#define NBSTRIPS 28 + +// forward declaration +class TCATSSpectra; + + + +using namespace std ; +enum reconstruction{NO,SECHS,GAUSS,BAR3,BAR4,BAR5}; +enum correction{NOcor,cor}; + +class TCATSPhysics : public TObject, public NPA::VDetector +{ + + public: // Constructor and Destructor + TCATSPhysics(); + ~TCATSPhysics(); + + 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 + vector<int> DetNumberX; + vector<int> StripX; + vector<double> ChargeX; + + // Vector of dim = number of CATS + vector<int> StripMaxX; + + + // Vector of dim = multiplicity of event on all detector + vector<int> DetNumberY; + vector<int> StripY; + vector<double> ChargeY; + + // Vector of dim = number of CATS + vector<int> StripMaxY; + + // 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 Buffer_X_Q[NBSTRIPS][NBDETECTOR];//! + double Buffer_Y_Q[NBSTRIPS][NBDETECTOR];//! + + int HitX; //! + int HitY; //! + + vector<reconstruction> ReconstructionMethodX; + vector<reconstruction> ReconstructionMethodY; + + + 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;//! + + 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;} + + + // Remove bad channel, calibrate the data and apply threshold + void PreTreat(); + + // 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();} + + void BuildPhysicalEvent(); + + void BuildSimplePhysicalEvent(); + + // Method related to the TSpectra classes, aimed at providing a framework for online applications + // Instantiate the Spectra class and the histogramm throught it + void InitSpectra(); + // Fill the spectra hold by the spectra class + void FillSpectra(); + // Used for Online mainly, perform check on the histo and for example change their color if issues are found + void CheckSpectra(); + // Used for Online only, clear all the spectra hold by the Spectra class + void ClearSpectra(); + + + // Those two method all to clear the Event Physics or Data + void ClearEventPhysics() {Clear();} + void ClearEventData() {m_EventData->Clear();} + + 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;} + + // Return false if the channel is disabled by user + bool IsValidChannel(const string DetectorType, const int Detector , const int channel); + + void InitializeStandardParameter(); + + void AddParameterToCalibrationManager(); + + void ReadAnalysisConfig(); + + void ReadConfiguration(string); + + void AddCATS(TVector3 C_X1_Y1, TVector3 C_X28_Y1, TVector3 C_X1_Y28, TVector3 C_X28_Y28); + + double AnalyseX(int ff); + + double AnalyseY(int ff); + + double CalculatePositionX( double CalculatedStripX, correction method); + + double CalculatePositionY( double CalculatedStripY, correction method); + + reconstruction ChooseReconstruction(int ff, TString type); + + reconstruction StringToEnum(string type); + + 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); + + + // 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;} + + + private: // Spectra Class + TCATSSpectra* m_Spectra;//! + + public: // Spectra Getter + map< vector<TString> , TH1*> GetSpectra(); + + + + ClassDef(TCATSPhysics,1) // CATSPhysics structure +}; + + +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/Tiara/TTiaraBarrelSpectra.cxx b/NPLib/Tiara/TTiaraBarrelSpectra.cxx index b9e072906..2b5e2a23b 100644 --- a/NPLib/Tiara/TTiaraBarrelSpectra.cxx +++ b/NPLib/Tiara/TTiaraBarrelSpectra.cxx @@ -148,7 +148,10 @@ string name ; AddHisto2D(name, name,360,0,180,1000,0,30,BaseFamily); } } - + //// E Theta //// + // Inner Barrel + name = "IB_ETHETA_CAL"; + AddHisto2D(name, name,360,0,180,1000,0,30,BaseFamily); } @@ -276,6 +279,10 @@ string name ; GetHisto(family,name) ->Fill(Theta*rad/deg,Physics->Strip_E[i]); + + name = "IB_ETHETA_CAL"; + GetHisto(family,name) + ->Fill(Theta*rad/deg,Physics->Strip_E[i]); } } -- GitLab