-
Hugo Jacob authoredHugo Jacob authored
TCATSPhysics.cxx 36.07 KiB
/* Copyright (C) 2009-2016 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 CATS_LOCAL;
// STL
#include <cmath>
#include <algorithm>
#include <sstream>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <utility>
#include <vector>
using namespace std;
// NPL
#include "RootInput.h"
#include "NPDetectorFactory.h"
#include "RootOutput.h"
#include "NPOptionManager.h"
// ROOT
#include "TChain.h"
#include "TF1.h"
#include "TGraph.h"
ClassImp(TCATSPhysics)
///////////////////////////////////////////////////////////////////////////
TCATSPhysics::TCATSPhysics(){
m_EventData = new TCATSData ;
m_PreTreatedData = new TCATSData ;
m_EventPhysics = this ;
m_NumberOfCATS = 0 ;
m_Spectra = NULL ;
m_Zproj = 0 ;
}
///////////////////////////////////////////////////////////////////////////
TCATSPhysics::~TCATSPhysics(){
}
///////////////////////////////////////////////////////////////////////////
void TCATSPhysics::PreTreat(){
ClearPreTreatedData();
gRandom->SetSeed(0);
// X
unsigned int sizeX = m_EventData->GetCATSMultX();
for(unsigned int i = 0 ; i < sizeX ; 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);
unsigned int stripX = -1;
//Inversion X
if( *(m_CATSXInversion[m_EventData->GetCATSDetX(i)-1].begin() + m_EventData->GetCATSStripX(i)-1) != m_EventData->GetCATSStripX(i) ){
stripX = *(m_CATSXInversion[m_EventData->GetCATSDetX(i)-1].begin() + m_EventData->GetCATSStripX(i)-1);
}
else {
stripX = m_EventData->GetCATSStripX(i);
}
m_PreTreatedData->SetStripX( m_EventData->GetCATSDetX(i) , stripX, QX);
}
}
}
// Y
unsigned int sizeY = m_EventData->GetCATSMultY();
for(unsigned int i = 0 ; i < sizeY ; 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);
unsigned int stripY = -1;
//Inversion Y
if( *(m_CATSYInversion[m_EventData->GetCATSDetY(i)-1].begin() + m_EventData->GetCATSStripY(i)-1) != m_EventData->GetCATSStripY(i) ){
stripY = *(m_CATSYInversion[m_EventData->GetCATSDetY(i)-1].begin() + m_EventData->GetCATSStripY(i)-1);
}
else {
stripY = m_EventData->GetCATSStripY(i) ;
}
m_PreTreatedData->SetStripY( m_EventData->GetCATSDetY(i), stripY, QY );
}
}
}
return;
}
/////////////////////////////////////////////////////////////////////////////
void TCATSPhysics::BuildSimplePhysicalEvent(){
BuildPhysicalEvent();
}
//////////////////////////////////////////////////////////////////////////////
void TCATSPhysics::BuildPhysicalEvent(){
// std::cout << "test 1" << std::endl;
if (NPOptionManager::getInstance()->IsReader() == true) {
m_EventData = &(**r_ReaderEventData);
}
//m_EventData->Dump();
// std::cout << "test 2" << std::endl;
PreTreat();
sizeX = m_PreTreatedData->GetCATSMultX() ;
for( unsigned short i = 0 ; i < sizeX; i++ ){
// Insert detector number in the set, if the key already exist, do nothing
DetectorHitX.insert(m_PreTreatedData->GetCATSDetX(i));
}
// std::cout << "test 3" << std::endl;
// Correspond to CATS with both X and Y
sizeY = m_PreTreatedData->GetCATSMultY() ;
for( unsigned short i = 0 ; i < sizeY ; i++ ){
// Insert detector number in the set, if the key already exist, do nothing
// Only if the detector was hit on X as well if(DetectorHitX.find(m_PreTreatedData->GetCATSDetY(i))!=DetectorHitX.end())
DetectorHit.insert(m_PreTreatedData->GetCATSDetY(i));
}
// The number of CATS hit, i.e. the number of CATS that we are going to analyse
//std::cout << "StripX mult " << sizeX << std::endl;
for( unsigned short i = 0 ; i < sizeX; i++ ){
StrX = m_PreTreatedData->GetCATSStripX(i);
NX = m_PreTreatedData->GetCATSDetX(i);
CATS_X_Q = m_PreTreatedData->GetCATSChargeX(i) ;
if(DetectorHit.find(NX)!=DetectorHit.end()){
MapX[NX].push_back(std::make_pair(StrX,CATS_X_Q));
QSumX[NX]+= CATS_X_Q;
if(MaxQX.find(NX)==MaxQX.end()|| MaxQX[NX].second < CATS_X_Q ){
MaxQX[NX] = make_pair(StrX,CATS_X_Q);
}
}
}
for( unsigned short i = 0 ; i < sizeY; i++ ){
StrY = m_PreTreatedData->GetCATSStripY(i);
NY = m_PreTreatedData->GetCATSDetY(i);
CATS_Y_Q = m_PreTreatedData->GetCATSChargeY(i) ;
if(DetectorHit.find(NY)!=DetectorHit.end()){
MapY[NY].push_back(std::make_pair(StrY,CATS_Y_Q));
QSumY[NY]+= CATS_Y_Q;
if(MaxQY.find(NY)==MaxQY.end()|| MaxQY[NY].second < CATS_Y_Q ){
MaxQY[NY] = make_pair(StrY,CATS_Y_Q);
}
}
}
// std::cout << "test 5" << std::endl;
for(auto &DetN : DetectorHit){
// std::cout << "test 6" << std::endl;
// Return the position in strip unit
// Convention: the collected charge is atrributed to the center of the strip
// (histogram convention) so that a reconstructed position for a single strip
// goes from strip index -0.5 to strip index +0.5
double PosX = ReconstructionFunctionX[DetN](MaxQX[DetN],MapX[DetN], QSumX[DetN]);
// std::cout << "test 7" << std::endl;
double PosY = ReconstructionFunctionY[DetN](MaxQY[DetN],MapY[DetN], QSumY[DetN]);
//std::cout << "test " << std::reduce(QsumSample[DetN].begin(),QsumSample[DetN].end()) /(QsumSample[DetN]).size() << std::endl;
StripNumberX.push_back(PosX);
StripNumberY.push_back(PosY);
DetNumber.push_back(DetN);
ChargeX.push_back(QSumX[DetN]);
ChargeY.push_back(QSumY[DetN]);
// a shift - -1 is made to have PosX in between -0.5 and 27.5
// for the following calculation of the position in the lab.
PosX = PosX;
PosY = PosY;
// sx and sy are the X and Y strip number between which the PosX and PosY are
int sx0 = (int) PosX;
int sx1 = sx0+1;
int sy0 = (int) PosY;
int sy1 = sy0+1;
if(PosX>-1000 && PosY>-1000 && sx0 > -1 && sx1 < 28 && sy0 > -1 && sy1 < 28){
// px and py are the x and y coordinate of strip sx and sy
double px0 = StripPositionX[DetN][sx0][sy0];
double px1 = StripPositionX[DetN][sx1][sy1];
double py0 = StripPositionY[DetN][sx0][sy0];
double py1 = StripPositionY[DetN][sx1][sy1];
// Positon [Detector] = <PosZ,<PosX,PosY>>
Positions[DetN] = make_pair(StripPositionZ[DetN], make_pair(px0+(px1-px0)*(PosX-sx0),py0+(py1-py0)*(PosY-sy0)));
PositionX.push_back(px0+(px1-px0)*(PosX-sx0));
PositionY.push_back(py0+(py1-py0)*(PosY-sy0));
PositionZ.push_back(StripPositionZ[DetN]);
}
}
// Sorting Positions depending on Z
// std::sort(Positions.begin(), Positions.end(),
// [&](const auto& Pos1, const auto& Pos2){
// return Pos1.first < Pos2.first;
// });
// At least two CATS need to gave back position in order to reconstruct on Target
if(Positions.size()>1){
double t = (m_Zproj-Positions[2].first)/(m_Zproj-Positions[1].first);
PositionOnTargetX= (Positions[2].second.first - Positions[1].second.first*t)/(1 - t);
PositionOnTargetY= (Positions[2].second.second - Positions[1].second.second*t)/(1 - t);
// PositionOnTargetY= Positions[2].second.second + (Positions[2].second.second-3-Positions[1].second.second)*t;
if(Mask1_Z != 0 && Mask2_Z != 0)
{
double tmask1 = (Mask1_Z-Positions[2].first)/(Positions[2].first-Positions[1].first);
double tmask2 = (Mask2_Z-Positions[2].first)/(Positions[2].first-Positions[1].first);
PositionOnMask1X= Positions[2].second.first + (Positions[2].second.first-Positions[1].second.first)*tmask1;
PositionOnMask1Y= Positions[2].second.second + (Positions[2].second.second-Positions[1].second.second)*tmask1;
PositionOnMask2X= Positions[2].second.first + (Positions[2].second.first-Positions[1].second.first)*tmask2;
PositionOnMask2Y= Positions[2].second.second + (Positions[2].second.second-Positions[1].second.second)*tmask2;
}
else{
PositionOnMask1X= -1000;
PositionOnMask1Y= -1000;
PositionOnMask2X= -1000;
PositionOnMask2Y= -1000;
}
}
else{
BeamDirection = TVector3 (1,0,0);
PositionOnTargetX = -1000 ;
PositionOnTargetY = -1000 ;
PositionOnMask1X= -1000;
PositionOnMask1Y= -1000;
PositionOnMask2X= -1000;
PositionOnMask2Y= -1000;
}
BeamDirection = GetBeamDirection();
// Does not meet the conditions for target position and beam direction
return;
}
void TCATSPhysics::SetTreeReader(TTreeReader* TreeReader) {
TCATSPhysicsReader::r_SetTreeReader(TreeReader);
}
///////////////////////////////////////////////////////////////////////////
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
void TCATSPhysics::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("CATSDetector");
vector<NPL::InputBlock*> blocksMask = parser.GetAllBlocksWithToken("MASK");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
cout << "//// " << blocksMask.size() << " masks found " << endl;
vector<string> token = {"X1_Y1","X28_Y1","X1_Y28","X28_Y28","CATSNumber"};
vector<string> tokenMask = {"Z","MaskNumber"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(token)){
TVector3 A = blocks[i]->GetTVector3("X1_Y1","mm");
TVector3 B = blocks[i]->GetTVector3("X28_Y1","mm"); TVector3 C = blocks[i]->GetTVector3("X1_Y28","mm");
TVector3 D = blocks[i]->GetTVector3("X28_Y28","mm");
UShort_t N = blocks[i]->GetInt("CATSNumber");
AddCATS(A,B,C,D,N);
}
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
for(unsigned int i = 0 ; i < blocksMask.size() ; i++){
if(blocksMask[i]->HasTokenList(tokenMask)){
AddMask(blocksMask[i]->GetDouble("Z","mm"),blocksMask[i]->GetInt("MaskNumber"));
}
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
}
}
InitializeStandardParameter();
ReadAnalysisConfig();
}
/////////////////////////////////////////////////////////////////////
// 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() ;
// Option to use the nptreereader anaysis
if (NPOptionManager::getInstance()->IsReader() == true) {
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
inputTreeReader->SetTree(inputChain);
}
// Option to use the standard npanalysis
else{
std::cout << "////////////////////////////// TEST" << std::endl;
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();
// Option to use the nptreereader anaysis
if (NPOptionManager::getInstance()->IsReader() == true) {
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
inputTreeReader->SetTree(inputChain);
}
// Option to use the standard npanalysis
else{
inputChain->SetBranchStatus( "CATS" , 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( "DetMaxX" , true );
inputChain->SetBranchStatus( "DetMaxY" , true );
inputChain->SetBranchStatus( "PositionX" , true );
inputChain->SetBranchStatus( "PositionY" , true );
inputChain->SetBranchStatus( "PositionZ" , true ); inputChain->SetBranchStatus( "StripNumberX" , true );
inputChain->SetBranchStatus( "StripNumberY" , 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, UShort_t N){
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[N] = OneDetectorStripPositionX ; StripPositionY[N] = OneDetectorStripPositionY ;
StripPositionZ[N] = OneDetectorStripPositionZ ;
}
void TCATSPhysics::AddMask(Double_t Z, UShort_t MaskNumber){
if(MaskNumber == 1){
Mask1_Z = Z;
}
else if(MaskNumber == 2){
Mask2_Z = Z;
}
else{
std::cout << "Wrong Number for MASKS" << std::endl;
}
}
///////////////////////////////////////////////////////////////
void TCATSPhysics::Clear(){
DetNumber.clear();
ChargeX.clear();
PositionX.clear();
ChargeY.clear();
PositionY.clear();
PositionZ.clear();
StripNumberX.clear();
StripNumberY.clear();
QSumX.clear();
QSumY.clear();
Positions.clear();
MapX.clear();
MapY.clear();
MaxQX.clear();
MaxQY.clear();
DetectorHit.clear();
DetectorHitX.clear();
}
////////////////////////////////////////////////////////////////////////////
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;
SetReconstructionMethod(i+1, "X", "AGAUSS");
SetReconstructionMethod(i+1, "Y", "AGAUSS");
}
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 (whatToDo == "INVERSIONX") {
AnalysisConfigFile >> DataBuffer;
cout << whatToDo << " " << DataBuffer;
int Detector = atoi(DataBuffer.substr(4,1).c_str());
for(unsigned int strip = 0; strip < 28; strip ++){
*(m_CATSXInversion[Detector-1].begin()+strip) = 27-strip;
}
}
else if (whatToDo == "INVERSIONY") {
AnalysisConfigFile >> DataBuffer;
cout << whatToDo << " " << DataBuffer;
int Detector = atoi(DataBuffer.substr(4,1).c_str());
for(unsigned int strip = 0; strip < 28; strip ++){
*(m_CATSYInversion[Detector-1].begin()+strip) = 27-strip;
}
}
else if (whatToDo == "RECONSTRUCTION_METHOD") {
AnalysisConfigFile >> DataBuffer;
cout << whatToDo << " " << DataBuffer ;
// DataBuffer is of form CATSNX
// Look for the CATS Number removing the first 4 letters and the trailling letter
string Duplicate = DataBuffer.substr(4); // Duplicate is of form NX
Duplicate.resize(Duplicate.size()-1); // Duplicate is of form
unsigned int CATSNumber = atoi(Duplicate.c_str());
// Look for the X or Y part of the Detector, Basically the last character
string XorY(string(1,DataBuffer[DataBuffer.size()-1])) ;
// Get the Reconstruction Methods Name
AnalysisConfigFile >> DataBuffer;
cout << " " << DataBuffer << endl ;
// Set the Reconstruction Methods using above information
SetReconstructionMethod(CATSNumber,XorY,DataBuffer);
}
else if (whatToDo == "ZPROJ") {
AnalysisConfigFile >> DataBuffer;
cout << whatToDo << " " << DataBuffer << endl ;
m_Zproj = atoi(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< string , TH1*> TCATSPhysics::GetSpectra() {
if(m_Spectra)
return m_Spectra->GetMapHisto();
else{
map< string , TH1*> empty;
return empty;
}
}
///////////////////////////////////////////////////////////////////////////
void TCATSPhysics::WriteSpectra(){
if(m_Spectra)
m_Spectra->WriteSpectra();
}
/////////////////////////////////////////////////////////////////////
// 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"+NPL::itoa(i+1)+"_X"+NPL::itoa(j+1)+"_Q","CATS_D"+NPL::itoa(i+1)+"_X"+NPL::itoa(j+1)+"_Q") ;
Cal->AddParameter("CATS", "D"+NPL::itoa(i+1)+"_Y"+NPL::itoa(j+1)+"_Q","CATS_D"+NPL::itoa(i+1)+"_Y"+NPL::itoa(j+1)+"_Q") ;
Cal->AddParameter("CATS", "D"+NPL::itoa(i+1)+"_X"+NPL::itoa(j+1),"CATS_D"+NPL::itoa(i+1)+"_X"+NPL::itoa(j+1)) ;
Cal->AddParameter("CATS", "D"+NPL::itoa(i+1)+"_Y"+NPL::itoa(j+1),"CATS_D"+NPL::itoa(i+1)+"_Y"+NPL::itoa(j+1)) ;
}
}
return;
}
////////////////////////////////////////////////////////////////
void TCATSPhysics::SetReconstructionMethod(unsigned int CATSNumber, string XorY, string MethodName){
if(XorY=="X"){
if(MethodName=="ASECH") ReconstructionFunctionX[CATSNumber] = &(AnalyticHyperbolicSecant);
//else if(MethodName=="FSECH") ReconstructionFunctionX[CATSNumber] = &(FittedHyperbolicSecant);
//else if(MethodName=="AGAUSS") ReconstructionFunctionX[CATSNumber] = &(AnalyticGaussian);
//else if(MethodName=="CENTROIDE") ReconstructionFunctionX[CATSNumber] = &(Centroide);
else cout <<"WARNING: Wrong name for reconsctuction Method, using default AGAUSS" << endl;
}
if(XorY=="Y"){
if(MethodName=="ASECH") ReconstructionFunctionY[CATSNumber] = &(AnalyticHyperbolicSecant);
//else if(MethodName=="FSECH") ReconstructionFunctionY[CATSNumber] = &(FittedHyperbolicSecant);
//else if(MethodName=="AGAUSS") ReconstructionFunctionY[CATSNumber] = &(AnalyticGaussian);
//else if(MethodName=="CENTROIDE") ReconstructionFunctionY[CATSNumber] = &(Centroide);
else cout <<"WARNING: Wrong name for reconsctuction Method, using default AGAUSS" << endl;
}
}
///////////////////////////////////////////////////////////////
TVector3 TCATSPhysics::GetBeamDirection(){
TVector3 Direction;
if(Positions.size() <2)return Direction;
else{
Direction = TVector3 (Positions[2].second.first-Positions[1].second.first ,
Positions[2].second.second-Positions[1].second.second,
Positions[2].first-Positions[1].first );
Direction.Unit();
}
return(Direction) ;
}
///////////////////////////////////////////////////////////////
TVector3 TCATSPhysics::GetPositionOnTarget(){
double Pi = 3.14159265;
TVector3 Position = TVector3 (GetPositionOnTargetX() ,
GetPositionOnTargetY() ,
GetPositionOnTargetX()*tan(m_TargetAngle*Pi/180));
return(Position) ;
}
////////////////////////////////////////////////////////////////////////
namespace CATS_LOCAL{
////////////////////////////////////////////////////////////////////
/*double AnalyticGaussian(std::pair<UShort_t,UShort_t>& MaxQ,std::vector<std::pair<UShort_t,UShort_t>>& Map, Double_t QSum){
/*double gauss = -1000;
double Q[3];
double StripPos[3];
for(int j = 0; j<3 ; j++){
Q[j] = 0;
StripPos[j] = 0;
}
if(MaxQ.first> 3 && MaxQ.first< 26){
// central value taken using the Strip with Max charge
Q[0] = Buffer_Q[MaxQ.first-1] ;
// Look at the next strip on the left
if(Buffer_Q[MaxQ.first-2]!=-1){
Q[1] = Buffer_Q[MaxQ.first-2];
StripPos[1] = MaxQ.first-2;
}
// Look at the next next strip on the left
else if(Buffer_Q[MaxQ.first-3]!=-1){
Q[1] = Buffer_Q[MaxQ.first-3];
StripPos[1] = MaxQ.first-3;
}
// Look at the next next next strip on the left
else if(Buffer_Q[MaxQ.first-4]!=-1){
Q[1] = Buffer_Q[MaxQ.first-4];
StripPos[1] = MaxQ.first-4;
}
// Look at the next strip on the right
if(Buffer_Q[MaxQ.first]!=-1){
Q[2] = Buffer_Q[MaxQ.first];
StripPos[2] = MaxQ.first;
}
// Look at the next next strip on the right
else if(Buffer_Q[MaxQ.first+1]!=-1){
Q[2] = Buffer_Q[MaxQ.first+1];
StripPos[2] = MaxQ.first+1; }
// Look at the next next next strip on the right
else if(Buffer_Q[MaxQ.first+2]!=-1){
Q[2] = Buffer_Q[MaxQ.first+2];
StripPos[2] = MaxQ.first+2;
}
}
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;
}
return gauss;
return 0;
}
///////////////////////////////////////////////////////////////
double Centroide(std::pair<UShort_t,UShort_t>& MaxQ,std::vector<std::pair<UShort_t,UShort_t>>& Map, Double_t QSum){
double Centroide = 0 ;
double ChargeTotal = 0;
unsigned int sizeQ = Buffer_Q.size();
for(unsigned int i = 0 ; i < sizeQ ; i++){
if(Buffer_Q[i]>0){
Centroide += (i)*Buffer_Q[i-1] ;
ChargeTotal+=Buffer_Q[i-1];
}
}
if(ChargeTotal>0) Centroide = Centroide / ChargeTotal ;
else {
Centroide = -1000 ;
}
return Centroide ;
return 0;
}
*/
/////////////////////////////////////////////////////////////////////
double AnalyticHyperbolicSecant(std::pair<UShort_t,UShort_t>& MaxQ,std::vector<std::pair<UShort_t,UShort_t>>& Map,Double_t QSum){
double sech = -1000 ;
// std::cout << "test AnH 1" << std::endl;
if(MaxQ.second > 0 && MaxQ.first > 2 && MaxQ.first<27){
// if(Buffer_Q[MaxQ.first-1+1]==0||Buffer_Q[MaxQ.first-1-1]==0)
// return sech;
// std::cout << "test AnH 2" << std::endl;
double q2 = MaxQ.second;
double q1 = 0,q3 = 0;
for(auto &strip : Map){
if(strip.first == MaxQ.first - 1){
q1 = strip.second;
}
else if(strip.first == MaxQ.first + 1){
q3 = strip.second;
}
} //std::cout << "test q " << q1 << " " << q2 << " " << q3 << std::endl;
double vs[6];
// std::cout << "test AnH 3" << std::endl;
if(q1 > 0 && q3 > 0)
{
// QsumSample[DetNum].push_back(QSum);
vs[0] = sqrt(q2/q3);
vs[1] = sqrt(q2/q1);
vs[2] = 0.5*(vs[0] + vs[1]);
vs[3] = log( vs[2] + sqrt(vs[2]*vs[2]-1.0) );
vs[4] = abs((vs[0] - vs[1])/(2.0*sinh(vs[3])));
vs[5] = 0.5*log( (1.0+vs[4])/(1.0-vs[4]) ) ;
// std::cout << "test AnH 4" << std::endl;
if ( q3>q1 )
sech = MaxQ.first + vs[5]/vs[3] ;
else
sech = MaxQ.first - vs[5]/vs[3] ;
//std::cout << "test sech " << sech << std::endl;
}
}
return sech ;
}
/////////////////////////////////////////////////////////////////////
/*double FittedHyperbolicSecant(std::pair<UShort_t,UShort_t>& MaxQ,std::vector<std::pair<UShort_t,UShort_t>>& Map, Double_t QSum){
// Warning: should not delete static variable
static TF1* f = new TF1("sechs","[0]/(cosh(TMath::Pi()*(x-[1])/[2])*cosh(TMath::Pi()*(x-[1])/[2]))",1,28);
// Help the fit by computing the position of the maximum by analytic method
double StartingPoint = AnalyticHyperbolicSecant(Buffer_Q,MaxQ.first);
// if analytic method fails then the starting point in strip max
if(StartingPoint==-1000) StartingPoint = MaxQ.first;
// Maximum is close to charge max, Mean value is close to Analytic one, typical width is 3.8 strip
f->SetParameters(Buffer_Q[MaxQ.first-1],StartingPoint,3.8);
static vector<double> y ;
static vector<double> q ;
y.clear(); q.clear();
double final_size = 0 ;
unsigned int sizeQ = Buffer_Q.size();
for(unsigned int i = 0 ; i < sizeQ ; i++){
if(Buffer_Q[i] > Buffer_Q[MaxQ.first-1]*0.2){
q.push_back(Buffer_Q[i]);
y.push_back(i+1);
final_size++;
}
}
// requiered at least 3 point to perfom a fit
if(final_size<3){
return -1000 ;
}
TGraph* g = new TGraph(q.size(),&y[0],&q[0]);
g->Fit(f,"QN0");
delete g;
return f->GetParameter(1) ;
return 0;
}
*/
////////////////////////////////////////////////////////////////////////
double fCATS_X_Q(const TCATSData* m_EventData , const int& i){
static string name; name = "CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetX(i) ) ;
name+= "_X" ;
name+= NPL::itoa( m_EventData->GetCATSStripX(i) ) ;
name+= "_Q";
return CalibrationManager::getInstance()->ApplyCalibration( name,
m_EventData->GetCATSChargeX(i) + gRandom->Rndm() );
//return CalibrationManager::getInstance()->ApplyCalibration( name,
// m_EventData->GetCATSChargeX(i) + gRandom->Rndm() );
//m_EventData->GetCATSChargeX(i) + gRandom->Rndm() - fCATS_Ped_X(m_EventData, i) );
}
////////////////////////////////////////////////////////////////////////
double fCATS_Y_Q(const TCATSData* m_EventData , const int& i){
static string name;
name = "CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetY(i) ) ;
name+= "_Y" ;
name+= NPL::itoa( m_EventData->GetCATSStripY(i) ) ;
name+= "_Q";
return CalibrationManager::getInstance()->ApplyCalibration( name ,
m_EventData->GetCATSChargeY(i) + gRandom->Rndm() );
//return CalibrationManager::getInstance()->ApplyCalibration( name ,
// m_EventData->GetCATSChargeY(i) + gRandom->Rndm() );
//m_EventData->GetCATSChargeY(i) + gRandom->Rndm() - fCATS_Ped_Y(m_EventData, i) );
}
////////////////////////////////////////////////////////////////////////
bool fCATS_Threshold_X(const TCATSData* m_EventData , const int& i){
static string name;
name = "CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetX(i) ) ;
name+= "_X" ;
name+= NPL::itoa( m_EventData->GetCATSStripX(i) );
return CalibrationManager::getInstance()->ApplyThreshold(name,
m_EventData->GetCATSChargeX(i));
}
////////////////////////////////////////////////////////////////////////
bool fCATS_Threshold_Y(const TCATSData* m_EventData , const int& i){
static string name;
name ="CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetY(i) ) ;
name+= "_Y" ;
name+= NPL::itoa( m_EventData->GetCATSStripY(i) );
return CalibrationManager::getInstance()->ApplyThreshold( name,
m_EventData->GetCATSChargeY(i));
}
////////////////////////////////////////////////////////////////////////
double fCATS_Ped_X(const TCATSData* m_EventData, const int& i){
static string name;
name = "CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetX(i) ) ;
name+= "_X" ;
name+= NPL::itoa( m_EventData->GetCATSStripX(i) ) ;
return CalibrationManager::getInstance()->GetPedestal(name);
}
////////////////////////////////////////////////////////////////////////
double fCATS_Ped_Y(const TCATSData* m_EventData, const int& i){
static string name;
name = "CATS/D" ;
name+= NPL::itoa( m_EventData->GetCATSDetY(i) ) ;
name+= "_Y" ;
name+= NPL::itoa( m_EventData->GetCATSStripY(i) );
return CalibrationManager::getInstance()->GetPedestal( name );
}
}
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TCATSPhysics::Construct(){
return (NPL::VDetector*) new TCATSPhysics();}
NPL::VTreeReader* TCATSPhysics::ConstructReader() { return (NPL::VTreeReader*)new TCATSPhysicsReader(); }
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_cats{
public:
proxy_cats(){
NPL::DetectorFactory::getInstance()->AddToken("CATSDetector","CATS");
NPL::DetectorFactory::getInstance()->AddDetector("CATSDetector",TCATSPhysics::Construct);
NPL::DetectorFactory::getInstance()->AddDetectorReader("CATSDetector", TCATSPhysics::ConstructReader);
}
};
proxy_cats p;
}