diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx index f8a1933c8f7c8712f0598ce13ce3dfc5042719de..0d3385c78ab1f0a26f92b00863484e2bbf2a670e 100644 --- a/NPLib/Core/NPDetectorManager.cxx +++ b/NPLib/Core/NPDetectorManager.cxx @@ -268,7 +268,7 @@ void NPL::DetectorManager::BuildPhysicalEvent(){ (it->second->*m_ClearEventPhysicsPtr)(); (it->second->*m_BuildPhysicalPtr)(); (it->second->*m_FillSpectra)(); - if(m_CheckSpectra) + if(m_CheckSpectra) { (it->second->*m_CheckSpectra)(); } } diff --git a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h index e63650cc71b1fbe501339488bfa1790c5463c32b..ae07f7afbbb134c6f5df82bd5d8bcdf1690e3697 100644 --- a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h +++ b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.h @@ -57,18 +57,18 @@ class TComptonTelescopePhysics : public TObject, public NPL::VDetector // DSSD int EventMultiplicity; vector<int> EventType; - //vector<int> TowerNumber; + vector<int> TowerNumber;// vector<int> DetectorNumber; vector<int> Strip_Front; vector<int> Strip_Back; - //vector<double> Strip_E; + vector<double> Strip_E;// vector<double> Strip_T; vector<double> Front_Energy; vector<double> Back_Energy; - //vector<double> Half_Energy; + vector<double> Half_Energy;// vector<double> Front_Time; vector<double> Back_Time; - //vector<bool> Same_FBTime; + vector<bool> Same_FBTime;// // Calorimeter vector<double> Calor_E; vector<double> Calor_T; @@ -240,8 +240,8 @@ class TComptonTelescopePhysics : public TObject, public NPL::VDetector Int_t m_CounterEvt[50]; //! Int_t m_CounterHit[50]; //! // physical events - vector<int> TowerNumber; - vector<double> Half_Energy; + //vector<int> TowerNumber; + //vector<double> Half_Energy; //vector<bool> Same_FBTime; }; diff --git a/NPLib/Detectors/Sofia/CMakeLists.txt b/NPLib/Detectors/Sofia/CMakeLists.txt index cf621e66f7216c94cc9e8ea5ec2962575e907bcd..65c96df0144859d82f0ae17a03d7f688b108d71c 100644 --- a/NPLib/Detectors/Sofia/CMakeLists.txt +++ b/NPLib/Detectors/Sofia/CMakeLists.txt @@ -5,6 +5,7 @@ add_custom_command(OUTPUT TSofTofWDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scrip add_custom_command(OUTPUT TSofTofWPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofTofWPhysics.h TSofTofWPhysicsDict.cxx TSofTofWPhysics.rootmap libNPSofia.dylib DEPENDS TSofTofWPhysics.h) add_custom_command(OUTPUT TSofMwpcDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofMwpcData.h TSofMwpcDataDict.cxx TSofMwpcData.rootmap libNPSofia.dylib DEPENDS TSofMwpcData.h) +add_custom_command(OUTPUT TSofMwpcPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofMwpcPhysics.h TSofMwpcPhysicsDict.cxx TSofMwpcPhysics.rootmap libNPSofia.dylib DEPENDS TSofMwpcPhysics.h) add_custom_command(OUTPUT TSofAtDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofAtData.h TSofAtDataDict.cxx TSofAtData.rootmap libNPSofia.dylib DEPENDS TSofAtData.h) add_custom_command(OUTPUT TSofAtPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofAtPhysics.h TSofAtPhysicsDict.cxx TSofAtPhysics.rootmap libNPSofia.dylib DEPENDS TSofAtPhysics.h) @@ -19,9 +20,9 @@ add_custom_command(OUTPUT TSofBeamIDDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts add_custom_command(OUTPUT TSofFissionFragmentDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofFissionFragment.h TSofFissionFragmentDict.cxx TSofFissionFragment.rootmap libNPSofia.dylib DEPENDS TSofFissionFragment.h) -add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofSciPhysics.cxx TSofSciPhysicsDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofAtPhysics.cxx TSofAtPhysicsDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTwimPhysics.cxx TSofTwimPhysicsDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx TSofBeamID.cxx TSofBeamIDDict.cxx TSofFissionFragment.cxx TSofFissionFragmentDict.cxx) +add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofSciPhysics.cxx TSofSciPhysicsDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofMwpcPhysics.cxx TSofMwpcPhysicsDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofAtPhysics.cxx TSofAtPhysicsDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTwimPhysics.cxx TSofTwimPhysicsDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx TSofBeamID.cxx TSofBeamIDDict.cxx TSofFissionFragment.cxx TSofFissionFragmentDict.cxx) target_link_libraries(NPSofia ${ROOT_LIBRARIES} NPCore NPPhysics) -install(FILES TSofSciData.h TSofSciPhysics.h TSofMwpcData.h TSofAtData.h TSofAtPhysics.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTwimPhysics.h TSofTofWData.h TSofTofWPhysics.h TSofBeamID.h TSofFissionFragment.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) +install(FILES TSofSciData.h TSofSciPhysics.h TSofMwpcData.h TSofMwpcPhysics.h TSofAtData.h TSofAtPhysics.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTwimPhysics.h TSofTofWData.h TSofTofWPhysics.h TSofBeamID.h TSofFissionFragment.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) diff --git a/NPLib/Detectors/Sofia/TSofAtPhysics.cxx b/NPLib/Detectors/Sofia/TSofAtPhysics.cxx index db82c5a7872e401fb7ca0388c39ea8e902779ea1..749e6744833fa8ea02f35e8f4b1d04628e127634 100644 --- a/NPLib/Detectors/Sofia/TSofAtPhysics.cxx +++ b/NPLib/Detectors/Sofia/TSofAtPhysics.cxx @@ -80,19 +80,10 @@ void TSofAtPhysics::BuildPhysicalEvent() { PreTreat(); unsigned int mysizeE = m_PreTreatedData->GetMultiplicity(); - if(mysizeE != 4) - return; - double E[4]={-1,-1,-1,-1}; for (UShort_t e = 0; e < mysizeE ; e++) { - E[m_PreTreatedData->GetAnodeNbr(e)-1] = m_PreTreatedData->GetEnergy(e); - } - - if(E[0]>0 && E[1]>0 && E[2]>0 && E[3]>0){ - for(int i=0; i<4; i++){ - AnodeNbr.push_back(i+1); - Energy.push_back(E[i]); - } + AnodeNbr.push_back(m_PreTreatedData->GetAnodeNbr(e)); + Energy.push_back(m_PreTreatedData->GetEnergy(e)); } } diff --git a/NPLib/Detectors/Sofia/TSofBeamID.cxx b/NPLib/Detectors/Sofia/TSofBeamID.cxx index f078b93b0458c121e1df9fca58fc28f7cc0e4e56..2129233354dc4c3f8b27960e21dc048136446d9e 100644 --- a/NPLib/Detectors/Sofia/TSofBeamID.cxx +++ b/NPLib/Detectors/Sofia/TSofBeamID.cxx @@ -52,6 +52,7 @@ void TSofBeamID::Clear() { fBeam_Brho = -1; fBeam_XS2 = -1000; fBeam_XCC = -1000; + fBeam_YCC = -1000; } @@ -70,5 +71,6 @@ void TSofBeamID::Dump() const { cout << "fBeam_Brho: " << fBeam_Brho << endl; cout << "fBeam_XS2: " << fBeam_XS2 << endl; cout << "fBeam_XCC: " << fBeam_XCC << endl; + cout << "fBeam_YCC: " << fBeam_YCC << endl; } diff --git a/NPLib/Detectors/Sofia/TSofBeamID.h b/NPLib/Detectors/Sofia/TSofBeamID.h index be129d80698ecc9cdcfe0fd7c958943fd0c5e532..ba810d5331781712b1ca4b787d61e33605792bf2 100644 --- a/NPLib/Detectors/Sofia/TSofBeamID.h +++ b/NPLib/Detectors/Sofia/TSofBeamID.h @@ -43,6 +43,7 @@ class TSofBeamID : public TObject { double fBeam_Brho; double fBeam_XS2; double fBeam_XCC; + double fBeam_YCC; ////////////////////////////////////////////////////////////// // Constructor and destructor @@ -75,6 +76,7 @@ class TSofBeamID : public TObject { inline void SetBrho(double val){fBeam_Brho = val;};//! inline void SetXS2(double val){fBeam_XS2 = val;};//! inline void SetXCC(double val){fBeam_XCC = val;};//! + inline void SetYCC(double val){fBeam_YCC = val;};//! ////////////////////// GETTERS //////////////////////// inline double GetZbeam() const {return fBeam_Z;}//! @@ -86,6 +88,7 @@ class TSofBeamID : public TObject { inline double GetBrho() const {return fBeam_Brho;}//! inline double GetXS2() const {return fBeam_XS2;}//! inline double GetXCC() const {return fBeam_XCC;}//! + inline double GetYCC() const {return fBeam_YCC;}//! ////////////////////////////////////////////////////////////// // Required for ROOT dictionnary diff --git a/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx new file mode 100644 index 0000000000000000000000000000000000000000..182332b55ed56692c1044605c8d6b3944bfbf8b5 --- /dev/null +++ b/NPLib/Detectors/Sofia/TSofMwpcPhysics.cxx @@ -0,0 +1,380 @@ +/***************************************************************************** + * Copyright (C) 2009-2020 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: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : November 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold SofMwpc Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "TSofMwpcPhysics.h" + +// STL +#include <sstream> +#include <iostream> +#include <cmath> +#include <stdlib.h> +#include <limits> +using namespace std; + +// NPL +#include "RootInput.h" +#include "RootOutput.h" +#include "NPDetectorFactory.h" +#include "NPOptionManager.h" + +// ROOT +#include "TChain.h" + +ClassImp(TSofMwpcPhysics) + + + /////////////////////////////////////////////////////////////////////////// +TSofMwpcPhysics::TSofMwpcPhysics() + : m_EventData(new TSofMwpcData), + m_PreTreatedData(new TSofMwpcData), + m_EventPhysics(this), + m_NumberOfDetectors(0){ + } + +/////////////////////////////////////////////////////////////////////////// +/// A usefull method to bundle all operation to add a detector +void TSofMwpcPhysics::AddDetector(TVector3 Pos){ + // In That simple case nothing is done + // Typically for more complex detector one would calculate the relevant + // positions (stripped silicon) or angles (gamma array) + DetPosX.push_back(Pos.X()); + DetPosY.push_back(Pos.Y()); + DetPosZ.push_back(Pos.Z()); + + m_NumberOfDetectors++; +} + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::AddDetector(double R, double Theta, double Phi){ + // Compute the TVector3 corresponding + TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta)); + // Call the cartesian method + AddDetector(Pos); +} + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::BuildSimplePhysicalEvent() { + BuildPhysicalEvent(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::BuildPhysicalEvent() { + + PreTreat(); + + unsigned int mysizeE = m_PreTreatedData->GetMultiplicity(); + + vector<double> ChargeArray; + ChargeArray.resize(288,0); + + int StripMaxX1[4]={0,0,0,0}; + int StripMaxX2[4]={0,0,0,0}; + int StripMaxY[4] ={0,0,0,0}; + double ChargeMaxX1[4]={0,0,0,0}; + double ChargeMaxX2[4]={0,0,0,0}; + double ChargeMaxY[4] ={0,0,0,0}; + + + for(int i=0; i<m_NumberOfDetectors; i++){ + Buffer_X1_Q.push_back(ChargeArray); + Buffer_X2_Q.push_back(ChargeArray); + Buffer_Y_Q.push_back(ChargeArray); + } + + for (UShort_t e = 0; e < mysizeE ; e++) { + int det = m_PreTreatedData->GetDetectorNbr(e); + int plane = m_PreTreatedData->GetPlane(e); + int strip = m_PreTreatedData->GetPad(e);// starts at 0 + double charge = m_PreTreatedData->GetCharge(e); + + // Xup for MWPC2 and 3 and X for MWPC1 and 4 + if(plane==1){ + Buffer_X1_Q[det-1][strip] = charge; + + if(charge>ChargeMaxX1[det-1]){ + ChargeMaxX1[det-1] = charge; + StripMaxX1[det-1] = strip; + } + } + // Xdown for MWPC2 and 3 and nothing for MWPC1 and 4 + else if(plane==2){ + Buffer_X2_Q[det-1][strip] = charge; + + if(charge>ChargeMaxX2[det-1]){ + ChargeMaxX2[det-1] = charge; + StripMaxX2[det-1] = strip; + } + } + // Y for all MWPCx + if(plane==3){ + Buffer_Y_Q[det-1][strip] = charge; + + if(charge>ChargeMaxY[det-1]){ + ChargeMaxY[det-1] = charge; + StripMaxY[det-1] = strip; + } + } + } + + double X1 = -100; + double X2 = -100; + double Y = -100; + for(int i=0; i<m_NumberOfDetectors; i++){ + double qleft = Buffer_X1_Q[i][StripMaxX1[i]-1]; + double qright = Buffer_X1_Q[i][StripMaxX1[i]+1]; + if(qleft>0 && qright>0){ + X1 = GetPositionX(ChargeMaxX1[i], StripMaxX1[i], qleft, qright); + X1 = X1 - DetPosX[i]; + } + qleft = Buffer_X2_Q[i][StripMaxX2[i]-1]; + qright = Buffer_X2_Q[i][StripMaxX2[i]+1]; + if(qleft>0 && qright>0){ + X2 = GetPositionX(ChargeMaxX2[i], StripMaxX2[i], qleft, qright); + X2 = X2 - DetPosX[i]; + } + double qdown = Buffer_Y_Q[i][StripMaxY[i]-1]; + double qup = Buffer_Y_Q[i][StripMaxY[i]+1]; + if(qdown>0 && qup>0){ + Y = GetPositionY(ChargeMaxY[i], StripMaxY[i], qdown, qup); + Y = Y - DetPosY[i]; + } + + DetectorNbr.push_back(i+1); + PositionX1.push_back(X1); + PositionX2.push_back(X2); + PositionY.push_back(Y); + } +} + +/////////////////////////////////////////////////////////////////////////// +double TSofMwpcPhysics::GetPositionX(double qmax, int padmax, double qleft, double qright){ + double fwx = 3.125; + double fSize = 200.0; + + double a3 = TMath::Pi() * fwx / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qleft) + TMath::Sqrt(qmax / qright)))); + + Double_t a2 = (a3 / TMath::Pi()) * TMath::ATanH((TMath::Sqrt(qmax / qleft) - TMath::Sqrt(qmax / qright)) / + (2 * TMath::SinH(TMath::Pi() * fwx / a3))); + + return (-1. * padmax * fwx + (fSize / 2) - (fwx / 2) - a2); // Left is positive and right negative +} + +/////////////////////////////////////////////////////////////////////////// +double TSofMwpcPhysics::GetPositionY(double qmax, int padmax, double qdown, double qup){ + double fwy = 3.125; + double fSize = 200.0; + double a3 = TMath::Pi() * fwy / (TMath::ACosH(0.5 * (TMath::Sqrt(qmax / qdown) + TMath::Sqrt(qmax / qup)))); + + Double_t a2 = (a3 / TMath::Pi()) * TMath::ATanH((TMath::Sqrt(qmax / qdown) - TMath::Sqrt(qmax / qup)) / + (2 * TMath::SinH(TMath::Pi() * fwy / a3))); + + return (padmax * fwy - (fSize / 2) + (fwy / 2) + a2); +} + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::PreTreat() { + // This method typically applies thresholds and calibrations + // Might test for disabled channels for more complex detector + + // clear pre-treated object + ClearPreTreatedData(); + + // instantiate CalibrationManager + static CalibrationManager* Cal = CalibrationManager::getInstance(); + + unsigned int mysize = m_EventData->GetMultiplicity(); + for (unsigned int i = 0; i < mysize ; ++i) { + m_PreTreatedData->SetDetectorNbr(m_EventData->GetDetectorNbr(i)); + m_PreTreatedData->SetPlane(m_EventData->GetPlane(i)); + m_PreTreatedData->SetPad(m_EventData->GetPad(i)); + m_PreTreatedData->SetCharge(m_EventData->GetCharge(i)); + } +} + + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::ReadAnalysisConfig() { + bool ReadingStatus = false; + + // path to file + string FileName = "./configs/ConfigSofMwpc.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigSofMwpc.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigSofMwpc.dat " << endl; + + // Save it in a TAsciiFile + TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig(); + asciiConfig->AppendLine("%%% ConfigSofMwpc.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" + string name = "ConfigSofMwpc"; + if (LineBuffer.compare(0, name.length(), name) == 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=="E_THRESHOLD") { + AnalysisConfigFile >> DataBuffer; + m_E_Threshold = atof(DataBuffer.c_str()); + cout << whatToDo << " " << m_E_Threshold << endl; + } + + else { + ReadingStatus = false; + } + } + } +} + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::Clear() { + DetectorNbr.clear(); + PositionX1.clear(); + PositionX2.clear(); + PositionY.clear(); + + Buffer_X1_Q.clear(); + Buffer_X2_Q.clear(); + Buffer_Y_Q.clear(); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::ReadConfiguration(NPL::InputParser parser) { + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SofMwpc"); + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> cart = {"POS"}; + vector<string> sphe = {"R","Theta","Phi"}; + + for(unsigned int i = 0 ; i < blocks.size() ; i++){ + if(blocks[i]->HasTokenList(cart)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// SofMwpc " << i+1 << endl; + + TVector3 Pos = blocks[i]->GetTVector3("POS","mm"); + AddDetector(Pos); + } + else if(blocks[i]->HasTokenList(sphe)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// SofMwpc " << i+1 << endl; + double R = blocks[i]->GetDouble("R","mm"); + double Theta = blocks[i]->GetDouble("Theta","deg"); + double Phi = blocks[i]->GetDouble("Phi","deg"); + AddDetector(R,Theta,Phi); + } + else{ + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } + + ReadAnalysisConfig(); +} + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::AddParameterToCalibrationManager() { + CalibrationManager* Cal = CalibrationManager::getInstance(); + + for(int sec = 0; sec < m_NumberOfDetectors; sec++){ + Cal->AddParameter("SofMwpc","SEC"+NPL::itoa(sec+1)+"_ALIGN","SofMwpc_SEC"+NPL::itoa(sec+1)+"_ALIGN"); + } +} + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::InitializeRootInputRaw() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchStatus("SofMwpc", true ); + inputChain->SetBranchAddress("SofMwpc", &m_EventData ); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::InitializeRootInputPhysics() { + TChain* inputChain = RootInput::getInstance()->GetChain(); + inputChain->SetBranchAddress("SofMwpc", &m_EventPhysics); +} + + + +/////////////////////////////////////////////////////////////////////////// +void TSofMwpcPhysics::InitializeRootOutput() { + TTree* outputTree = RootOutput::getInstance()->GetTree(); + outputTree->Branch("SofMwpc", "TSofMwpcPhysics", &m_EventPhysics); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPL::VDetector* TSofMwpcPhysics::Construct() { + return (NPL::VDetector*) new TSofMwpcPhysics(); +} + + + +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern "C"{ + class proxy_SofMwpc{ + public: + proxy_SofMwpc(){ + NPL::DetectorFactory::getInstance()->AddToken("SofMwpc","SofMwpc"); + NPL::DetectorFactory::getInstance()->AddDetector("SofMwpc",TSofMwpcPhysics::Construct); + } + }; + + proxy_SofMwpc p_SofMwpc; +} + diff --git a/NPLib/Detectors/Sofia/TSofMwpcPhysics.h b/NPLib/Detectors/Sofia/TSofMwpcPhysics.h new file mode 100644 index 0000000000000000000000000000000000000000..14935bd46a5c3700a0792f68bda9b1e301308c73 --- /dev/null +++ b/NPLib/Detectors/Sofia/TSofMwpcPhysics.h @@ -0,0 +1,162 @@ +#ifndef TSofMwpcPHYSICS_H +#define TSofMwpcPHYSICS_H +/***************************************************************************** + * Copyright (C) 2009-2020 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: Pierre Morfouace contact address: pierre.morfouace2@cea.fr * + * * + * Creation Date : November 2020 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class hold SofMwpc Treated data * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <vector> +#include <map> +#include <string> +using namespace std; + +// ROOT headers +#include "TObject.h" +#include "TH1.h" +#include "TVector3.h" +#include "TSpline.h" +// NPTool headers +#include "TSofMwpcData.h" +#include "NPCalibrationManager.h" +#include "NPVDetector.h" +#include "NPInputParser.h" + + + +class TSofMwpcPhysics : public TObject, public NPL::VDetector { + ////////////////////////////////////////////////////////////// + // constructor and destructor + public: + TSofMwpcPhysics(); + ~TSofMwpcPhysics() {}; + + + ////////////////////////////////////////////////////////////// + // Inherited from TObject and overriden to avoid warnings + public: + void Clear(); + void Clear(const Option_t*) {}; + + + ////////////////////////////////////////////////////////////// + // data obtained after BuildPhysicalEvent() and stored in + // output ROOT file + public: + vector<int> DetectorNbr; + vector<double> PositionX1; + vector<double> PositionX2; + vector<double> PositionY; + + private: + vector<double> DetPosX;//! + vector<double> DetPosY;//! + vector<double> DetPosZ;//! + vector<vector<double>> Buffer_X1_Q;//! + vector<vector<double>> Buffer_X2_Q;//! + vector<vector<double>> Buffer_Y_Q;//! + + /// A usefull method to bundle all operation to add a detector + void AddDetector(TVector3 Pos); + void AddDetector(double R, double Theta, double Phi); + + ////////////////////////////////////////////////////////////// + // methods inherited from the VDetector ABC class + public: + // read stream from ConfigFile to pick-up detector parameters + void ReadConfiguration(NPL::InputParser); + + // add parameters to the CalibrationManger + void AddParameterToCalibrationManager(); + + // method called event by event, aiming at extracting the + // physical information from detector + void BuildPhysicalEvent(); + + // same as BuildPhysicalEvent() method but with a simpler + // treatment + void BuildSimplePhysicalEvent(); + + // same as above but for online analysis + void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();}; + + // activate raw data object and branches from input TChain + // in this method mother branches (Detector) AND daughter leaves + // (fDetector_parameter) have to be activated + void InitializeRootInputRaw(); + + // activate physics data object and branches from input TChain + // in this method mother branches (Detector) AND daughter leaves + // (fDetector_parameter) have to be activated + void InitializeRootInputPhysics(); + + // create branches of output ROOT file + void InitializeRootOutput(); + + // clear the raw and physical data objects event by event + void ClearEventPhysics() {Clear();} + void ClearEventData() {m_EventData->Clear();} + + + ////////////////////////////////////////////////////////////// + // specific methods to SofMwpc array + public: + // remove bad channels, calibrate the data and apply thresholds + void PreTreat(); + + // clear the pre-treated object + void ClearPreTreatedData() {m_PreTreatedData->Clear();} + + // read the user configuration file. If no file is found, load standard one + void ReadAnalysisConfig(); + + // give and external TSofMwpcData object to TSofMwpcPhysics. + // needed for online analysis for example + void SetRawDataPointer(TSofMwpcData* rawDataPointer) {m_EventData = rawDataPointer;} + + double GetPositionX(double qmax, int padmax, double qleft, double qright); + double GetPositionY(double qmax, int padmax, double qdown, double qup); + + // objects are not written in the TTree + private: + TSofMwpcData* m_EventData; //! + TSofMwpcData* m_PreTreatedData; //! + TSofMwpcPhysics* m_EventPhysics; //! + + // getters for raw and pre-treated data object + public: + TSofMwpcData* GetRawData() const {return m_EventData;} + TSofMwpcData* GetPreTreatedData() const {return m_PreTreatedData;} + + // parameters used in the analysis + private: + double m_E_Threshold; //! + + // number of detectors + private: + int m_NumberOfDetectors; //! + + // Static constructor to be passed to the Detector Factory + public: + static NPL::VDetector* Construct(); + + ClassDef(TSofMwpcPhysics,1) // SofMwpcPhysics structure +}; +#endif diff --git a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx index 35a50d1722b53265f04e027b2d746c495571cba7..b2f606d3439d7a4a0059f1911f931b18d168117a 100644 --- a/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx +++ b/NPLib/Detectors/Sofia/TSofTrimPhysics.cxx @@ -48,7 +48,9 @@ TSofTrimPhysics::TSofTrimPhysics() m_PreTreatedData(new TSofTrimData), m_EventPhysics(this), m_NumberOfDetectors(0), - m_Beta(-1), + m_Beta(-1), + m_IsSplineSectionDriftTime(false), + m_IsSplineSectionAngle(false), m_BetaNorm(0.838), m_NumberOfSections(3), m_NumberOfAnodesPaired(3), @@ -194,7 +196,7 @@ void TSofTrimPhysics::BuildPhysicalEvent() { Ep3[i] = norm3 * Ep3[i] / (p0_3[i] + p1_3[i]*TMath::Power(m_Beta, -5./3.)); // Angle correction per pair: spline - Ep1[i] = Ep1[i] / fcorr_EvsA[i][0]->Eval(Ddt) * fcorr_EvsA[i][0]->Eval(0); + /*Ep1[i] = Ep1[i] / fcorr_EvsA[i][0]->Eval(Ddt) * fcorr_EvsA[i][0]->Eval(0); Ep2[i] = Ep2[i] / fcorr_EvsA[i][1]->Eval(Ddt) * fcorr_EvsA[i][1]->Eval(0); Ep3[i] = Ep3[i] / fcorr_EvsA[i][2]->Eval(Ddt) * fcorr_EvsA[i][2]->Eval(0); @@ -202,12 +204,17 @@ void TSofTrimPhysics::BuildPhysicalEvent() { Ep1[i] = Ep1[i] / fcorr_EvsDT[i][0]->Eval(DTp1[i]) * fcorr_EvsDT[i][0]->Eval(3000); Ep2[i] = Ep2[i] / fcorr_EvsDT[i][1]->Eval(DTp2[i]) * fcorr_EvsDT[i][1]->Eval(3000); Ep3[i] = Ep3[i] / fcorr_EvsDT[i][2]->Eval(DTp3[i]) * fcorr_EvsDT[i][2]->Eval(3000); - +*/ // Summing up Anode Energy per section Esec[i] = (Ep1[i] + Ep2[i] + Ep3[i])/3; + // Angle correction per section: spline + if(m_IsSplineSectionAngle) + Esec[i] = Esec[i] / fcorr_sec_angle[i]->Eval(Ddt) * fcorr_sec_angle[i]->Eval(0); + // 2nd DT correction per section: spline - Esec[i] = Esec[i] / fcorr_sec[i]->Eval(DTp2[i]) * fcorr_sec[i]->Eval(0); + if(m_IsSplineSectionDriftTime) + Esec[i] = Esec[i] / fcorr_sec_dt[i]->Eval(DTp2[i]) * fcorr_sec_dt[i]->Eval(0); // Section ALignement Esec[i] = Cal->ApplyCalibration("SofTrim/SEC"+NPL::itoa(i+1)+"_ALIGN",Esec[i]); @@ -308,17 +315,25 @@ void TSofTrimPhysics::ReadAnalysisConfig() { else if (whatToDo=="SPLINE_PAIR_DT_PATH") { AnalysisConfigFile >> DataBuffer; m_SPLINE_PAIR_DT_PATH = DataBuffer; - cout << "*** Loading Spline for Dritf Time correction per pair ***" << endl; + cout << "*** Loading Spline for Drtft Time correction per pair ***" << endl; LoadSplinePairDriftTime(); } else if (whatToDo=="SPLINE_SECTION_DT_PATH") { AnalysisConfigFile >> DataBuffer; m_SPLINE_SECTION_DT_PATH = DataBuffer; - cout << "*** Loading Spline for Dritf Time correction per section ***" << endl; + cout << "*** Loading Spline for Drift Time correction per section ***" << endl; LoadSplineSectionDriftTime(); } + else if (whatToDo=="SPLINE_SECTION_ANGLE_PATH") { + AnalysisConfigFile >> DataBuffer; + m_SPLINE_SECTION_ANGLE_PATH = DataBuffer; + cout << "*** Loading Spline for Angle correction per section ***" << endl; + LoadSplineSectionAngle(); + } + + else if (whatToDo=="E_THRESHOLD") { AnalysisConfigFile >> DataBuffer; @@ -389,22 +404,42 @@ void TSofTrimPhysics::LoadSplinePairDriftTime(){ /////////////////////////////////////////////////////////////////////////// void TSofTrimPhysics::LoadSplineSectionDriftTime(){ - TString filename = m_SPLINE_SECTION_DT_PATH; - TFile* ifile = new TFile(filename,"read"); + TString filename = m_SPLINE_SECTION_DT_PATH; + TFile* ifile = new TFile(filename,"read"); if(ifile->IsOpen()){ cout << "Loading splines..." << endl; for(int s=0; s<m_NumberOfSections; s++){ TString splinename = Form("spline_sec%i",s+1); - fcorr_sec[s] = (TSpline3*) ifile->FindObjectAny(splinename); - cout << fcorr_sec[s]->GetName() << endl; + fcorr_sec_dt[s] = (TSpline3*) ifile->FindObjectAny(splinename); + cout << fcorr_sec_dt[s]->GetName() << endl; } + m_IsSplineSectionDriftTime = true; } else cout << "File " << filename << " not found!" << endl; ifile->Close(); } +/////////////////////////////////////////////////////////////////////////// +void TSofTrimPhysics::LoadSplineSectionAngle(){ + TString filename = m_SPLINE_SECTION_ANGLE_PATH; + TFile* ifile = new TFile(filename,"read"); + + if(ifile->IsOpen()){ + cout << "Loading splines..." << endl; + for(int s=0; s<m_NumberOfSections; s++){ + TString splinename = Form("spline_sec%i",s+1); + + fcorr_sec_angle[s] = (TSpline3*) ifile->FindObjectAny(splinename); + cout << fcorr_sec_angle[s]->GetName() << endl; + } + m_IsSplineSectionAngle = true; + } + else + cout << "File " << filename << " not found!" << endl; + ifile->Close(); +} /////////////////////////////////////////////////////////////////////////// void TSofTrimPhysics::Clear() { diff --git a/NPLib/Detectors/Sofia/TSofTrimPhysics.h b/NPLib/Detectors/Sofia/TSofTrimPhysics.h index 923b46288ada001d9abb38d94f3c266a76eeaade..7053b20fd4188c2ddbdd6d471f44a4dbe4bbb10e 100644 --- a/NPLib/Detectors/Sofia/TSofTrimPhysics.h +++ b/NPLib/Detectors/Sofia/TSofTrimPhysics.h @@ -131,6 +131,7 @@ class TSofTrimPhysics : public TObject, public NPL::VDetector { void LoadSplinePairAngle(); void LoadSplinePairDriftTime(); void LoadSplineSectionDriftTime(); + void LoadSplineSectionAngle(); void SetBeta(double beta) {m_Beta = beta;} double GetBeta() {return m_Beta;} @@ -156,11 +157,15 @@ class TSofTrimPhysics : public TObject, public NPL::VDetector { string m_SPLINE_PAIR_ANGLE_PATH; //! string m_SPLINE_PAIR_DT_PATH; //! string m_SPLINE_SECTION_DT_PATH; //! + string m_SPLINE_SECTION_ANGLE_PATH; //! + bool m_IsSplineSectionDriftTime; //! + bool m_IsSplineSectionAngle; //! private: TSpline3* fcorr_EvsA[3][3]; //! TSpline3* fcorr_EvsDT[3][3]; //! - TSpline3* fcorr_sec[3]; //! + TSpline3* fcorr_sec_angle[3]; //! + TSpline3* fcorr_sec_dt[3]; //! // number of detectors private: diff --git a/Projects/ComptonTelescope/online/src/DecodeD.cpp b/Projects/ComptonTelescope/online/src/DecodeD.cpp index 3ef51c8149014f24bed107e60cabd8fb34956b30..3a09dd750da7436982bdd3c75c199e770ad85892 100644 --- a/Projects/ComptonTelescope/online/src/DecodeD.cpp +++ b/Projects/ComptonTelescope/online/src/DecodeD.cpp @@ -100,7 +100,7 @@ __EVENTTYPE__* DecodeD::getEvent() void DecodeD::decodeEvent() { - this -> Clear(); + /*this -> */Clear(); switch (datatype) { case D_ROOT: if (cursor < length) { diff --git a/Projects/ComptonTelescope/online/src/online_coinc.cpp b/Projects/ComptonTelescope/online/src/online_coinc.cpp index a8ad64d2676ed419b98101a0e2e2cd48aae51346..3b6cedac406382013936dbc79bb7685b8364181c 100644 --- a/Projects/ComptonTelescope/online/src/online_coinc.cpp +++ b/Projects/ComptonTelescope/online/src/online_coinc.cpp @@ -16,6 +16,10 @@ #include "DecodeD.h" #include "DecodeT.h" +#define __RUN__ 3 + +#define __1DET__ + #define __TEST_ZONE__ #undef __TEST_ZONE__ @@ -105,15 +109,47 @@ int main(int argc, char** argv) auto ccamPhys = (TComptonTelescopePhysics*) m_NPDetectorManager->GetDetector("ComptonTelescope"); ccamPhys->SetRawDataPointer(ccamData); +#ifdef __1DET__ + ifstream is; + is.open("/disk/proto-data/data/20210510_Bi207/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-05-10_07_41_50.raw"); + is.seekg(0, ios::end); + int len = is.tellg(); + is.seekg(0, ios::beg); + char* buff = new char[len]; + is.read(buff, len); + is.close(); + DecodeR* DR = new DecodeR(false, buff); + while(DR -> getCursor() < len) { + DR -> decodeBlobMFM(); + m_NPDetectorManager->ClearEventPhysics(); + m_NPDetectorManager->ClearEventData(); + ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), 64); + m_NPDetectorManager->BuildPhysicalEvent(); + m_OutputTree->Fill(); + m_NPDetectorManager->CheckSpectraServer(); + } + delete DR; + delete [] buff; + m_NPDetectorManager -> WriteSpectra(); +#else + // read data file/flux and fill ccamData object std::cout << "Reading data\n"; DecodeR* DR = new DecodeR(false); // Instantiates DecodeR object reading calorimeter data flux DecodeT* DT = new DecodeT(false); // Instantiates DecodeT object reading trigger data flux DecodeD* DD = new DecodeD(false); // Instantiates DecodeD object reading DSSSD(s) data flux // newframe_t* event; - //DD -> setTree("/disk/proto-data/data/20210304_run2/bb7_3309-7_cs137_20210304_14h35_conv.root"); + #if __RUN__ == 0 + DD -> setTree("../data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root"); + #elif __RUN__ == 1 + DD -> setTree("/disk/proto-data/data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root"); + #elif __RUN__ == 2 + DD -> setTree("/disk/proto-data/data/20210304_run2/bb7_3309-7_cs137_20210304_14h35_conv.root"); + #elif __RUN__ == 3 DD -> setTree("/disk/proto-data/data/20210305_run3/bb7_3309-7_cs137_20210305_14h53_conv.root"); - //DD -> setTree("../data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root"); + #elif __RUN__ == 4 + DD -> setTree("/disk/proto-data/data/20210407_run4/bb7_3309-7_cs137_20210407_14h53_conv.root"); + #endif int dlen = DD -> getLength(); int i = 0;// ROSMAP files loop counter @@ -126,8 +162,16 @@ int main(int argc, char** argv) ifstream iros, itrig; cout << "Loading data files " << std::flush; + #if __RUN__ == 0 + itrig.open("../data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary); + #elif __RUN__ == 1 + itrig.open("/disk/proto-data/data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary); + #elif __RUN__ == 2 + #elif __RUN__ == 3 itrig.open("/disk/proto-data/data/20210305_run3/mfm_trigger_20210305_run3.raw", ios::binary); - //itrig.open("../data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary); + #elif __RUN__ == 4 + itrig.open("/disk/proto-data/data/20210407_run4/mfm_trigger_20210407_run4.raw", ios::binary); + #endif itrig.seekg(0, ios::end); int tlen = itrig.tellg(); itrig.seekg(0, ios::beg); @@ -136,8 +180,16 @@ int main(int argc, char** argv) itrig.close(); cout << "... " << std::flush; + #if __RUN__ == 0 + iros.open("../data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary); + #elif __RUN__ == 1 + iros.open("/disk/proto-data/data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary); + #elif __RUN__ == 2 + #elif __RUN__ == 3 iros.open("/disk/proto-data/data/20210305_run3/mfm_rosmap_20210305_run3.raw", ios::binary); - //iros.open("../data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary); + #elif __RUN__ == 4 + iros.open("/disk/proto-data/data/20210407_run4/mfm_rosmap_20210407_run4.raw", ios::binary); + #endif iros.seekg(0, ios::end); int rlen = iros.tellg(); iros.seekg(0, ios::beg); @@ -157,6 +209,22 @@ int main(int argc, char** argv) resetCount = DT->getResetCount() - resetCount; resetCount = -resetCount; // T B C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! cout << "Found reset count: " << resetCount << endl; + + if (resetCount == 0) { + string ans; + cout << "reset count is 0. Continue ? y/[n]"; + cin >> ans; + if (ans != "y" and ans != "Y") { + DT -> setRaw(tbuff); + DT -> decodeBlobMFM(); + resetCount = DT -> getResetCount(); + while (not(DT->hasTrigged(2))) { + DT -> decodeBlobMFM(); + } + resetCount = DT->getResetCount() - resetCount; + cout << "Found reset count: " << resetCount << endl; + } + } int cr, cd, tr, td; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ @@ -261,7 +329,7 @@ int main(int argc, char** argv) DR -> decodeBlobMFM(); //#ifdef __TEST_ZONE__ - cout << "Entering test zone." << endl; +// cout << "Entering test zone." << endl; //#else DD -> rewind(); @@ -271,7 +339,7 @@ int main(int argc, char** argv) tr = DR -> getTime(); td = DD -> getTime(); - while(DR -> getCursor() < rlen and DD -> getCursor() < dlen) + while(DR -> getCursor() < rlen and DD -> getCursor() < dlen and c < 10000) { if (cr == cd) { #ifndef __TEST_ZONE__ @@ -283,7 +351,6 @@ int main(int argc, char** argv) //DR -> Dump(); //DD -> Dump(); c++; - cout << cc << " " << c /*<< "(" << cr << ", " << cd << ") : " << tr << " " << td*/ << "\n"; // Clear raw and physics data m_NPDetectorManager->ClearEventPhysics(); m_NPDetectorManager->ClearEventData(); @@ -307,6 +374,7 @@ int main(int argc, char** argv) m_OutputTree->Fill(); cc++; #endif + cout << cc << " " << c /*<< "(" << cr << ", " << cd << ") : " << tr << " " << td*/ << " |\t"; } // check spectra @@ -353,7 +421,7 @@ int main(int argc, char** argv) deltaT->Write(); } fout->Close(); - +#endif // Essential #if __cplusplus > 199711L && NPMULTITHREADING m_NPDetectorManager->StopThread(); diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx index f33a731ad288ab2dcdc832590078dfaa4fc379ed..dea167d58ae1d16503d5a7e7bcca23df9044d1ca 100644 --- a/Projects/s455/Analysis.cxx +++ b/Projects/s455/Analysis.cxx @@ -43,6 +43,7 @@ void Analysis::Init(){ SofTwim= (TSofTwimPhysics*) m_DetectorManager->GetDetector("SofTwim"); SofTofW= (TSofTofWPhysics*) m_DetectorManager->GetDetector("SofTofW"); SofAt= (TSofAtPhysics*) m_DetectorManager->GetDetector("SofAt"); + SofMwpc= (TSofMwpcPhysics*) m_DetectorManager->GetDetector("SofMwpc"); InitParameter(); InitOutputBranch(); @@ -74,7 +75,6 @@ void Analysis::TreatEvent(){ void Analysis::FissionFragmentAnalysis(){ unsigned int softofw_size = SofTofW->PlasticNbr.size(); unsigned int softwim_size = SofTwim->SectionNbr.size(); - unsigned int sofat_size = SofAt->Energy.size(); double TOF_CC[2]; double Plastic[2]; @@ -111,7 +111,7 @@ void Analysis::FissionFragmentAnalysis(){ } - if(softofw_size==2 && softwim_size==2 && sofat_size==4){ + if(softofw_size==2 && softwim_size==2){ for(unsigned int i=0; i< softofw_size; i++){ TOF_CC[i] = SofTofW->CalTof[i]; Plastic[i] = SofTofW->PlasticNbr[i]; @@ -205,7 +205,7 @@ void Analysis::FissionFragmentAnalysis(){ // Z calibration // double p0 = -1.1072; - double p1 = 0.27517; + double p1 = 0.27517*1.001; double Z1=-1; double Z2=-1; double Zsum=-1; @@ -247,6 +247,11 @@ void Analysis::FissionFragmentAnalysis(){ //////////////////////////////////////////////////////////////////////////////// void Analysis::BeamAnalysis(){ unsigned int sofsci_size = SofSci->DetectorNbr.size(); + double Z_p0 = -8.57894; + double Z_p1 = 0.560421; + double Y_p0 = 78.0711; + double Y_p1 = 0.022669; + if(sofsci_size==2){ double beta = SofSci->Beta[0]; //cout << "Set beta to " << beta << endl; @@ -262,10 +267,18 @@ void Analysis::BeamAnalysis(){ double velocity_mns = SofSci->VelocityMNs[0]; double Beta = SofSci->Beta[0]; double XS2 = SofSci->PosMm[0]; - double XCC = SofSci->PosMm[1]; + //double XCC = SofSci->PosMm[1]; + double XCC=0; + double YCC=0; + for(unsigned int i=0; i<SofMwpc->DetectorNbr.size(); i++){ + if(SofMwpc->DetectorNbr[i]==1){ + XCC = SofMwpc->PositionX1[i]; + YCC = SofMwpc->PositionY[i]; + } + } double LS2; - LS2 = fLS2_0*(1 + fK_LS2*Theta); + LS2 = fLS2_0;//*(1 + fK_LS2*Theta); velocity_mns = LS2/TofFromS2; Beta = velocity_mns * m/ns / NPUNITS::c_light; double Gamma = 1./(TMath::Sqrt(1 - TMath::Power(Beta,2))); @@ -273,6 +286,8 @@ void Analysis::BeamAnalysis(){ double AoQ = Brho / (3.10716*Gamma*Beta); double A = AoQ * Qmax; + Zbeam = Z_p0 + Z_p1*sqrt(Zbeam); + Zbeam = Zbeam/(Y_p0 + Y_p1*YCC)*Y_p0; // Filling Beam tree SofBeamID->SetZbeam(Zbeam); SofBeamID->SetQmax(rand.Gaus(Qmax,0.15)); @@ -283,6 +298,7 @@ void Analysis::BeamAnalysis(){ SofBeamID->SetBrho(Brho); SofBeamID->SetXS2(XS2); SofBeamID->SetXCC(XCC); + SofBeamID->SetYCC(YCC); } } } @@ -295,6 +311,11 @@ void Analysis::LoadCut(){ TString cutfile; TFile* file; for(int i=0; i<3; i++){ + // Q=77 + rootfile = Form("cutsec%iQ77.root",i+1); + cutfile = input_path + rootfile; + file = new TFile(cutfile); + cutQ77[i] = (TCutG*) file->Get(Form("cutsec%iQ77",i+1)); // Q=78 rootfile = Form("cutsec%iQ78.root",i+1); cutfile = input_path + rootfile; @@ -310,11 +331,6 @@ void Analysis::LoadCut(){ cutfile = input_path + rootfile; file = new TFile(cutfile); cutQ80[i] = (TCutG*) file->Get(Form("cutsec%iQ80",i+1)); - // Q=81 - rootfile = Form("cutsec%iQ81.root",i+1); - cutfile = input_path + rootfile; - file = new TFile(cutfile); - cutQ81[i] = (TCutG*) file->Get(Form("cutsec%iQ81",i+1)); } } @@ -350,14 +366,17 @@ int Analysis::DetermineQmax(){ double Theta = SofTrim->Theta[i]; double Esec = SofTrim->EnergySection[i]; - if(cutQ78[SecNbr-1]->IsInside(Theta,Esec)) + + if(cutQ77[SecNbr-1]->IsInside(Theta,Esec)) + Qsec[SecNbr-1] = 77; + else if(cutQ78[SecNbr-1]->IsInside(Theta,Esec)) Qsec[SecNbr-1] = 78; else if(cutQ79[SecNbr-1]->IsInside(Theta,Esec)) Qsec[SecNbr-1] = 79; else if(cutQ80[SecNbr-1]->IsInside(Theta,Esec)) Qsec[SecNbr-1] = 80; - else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec)) - Qsec[SecNbr-1] = 81; + //else if(cutQ81[SecNbr-1]->IsInside(Theta,Esec)) + //Qsec[SecNbr-1] = 81; } Qmax = max(Qsec[0],Qsec[1]); @@ -372,18 +391,24 @@ void Analysis::End(){ //////////////////////////////////////////////////////////////////////////////// void Analysis::InitParameter(){ - fLS2_0 = 136.3706933; - fDS2 = 9500; - //fDCC = -30000; + fLS2_0 = 135.614; + fDS2 = 8000; fDCC = -10000; - fK_LS2 = -2.5e-8; + fK_LS2 = -30e-8; + + + //fBrho0 = 14.1008; // 238U run 367 + //fBrho0 = 12.8719; // 238U run 368 //fBrho0 = 12.3255; // 238U run 369 - fBrho0 = 10.8183; // 182Hg - //fBrho0 = 10.6814; // 180Hg + //fBrho0 = 12.3352; // 238U run 419 + //fBrho0 = 10.9476; // 189Pb + //fBrho0 = 10.8183; // 182Hg + fBrho0 = 10.6814; // 180Hg //fBrho0 = 10.8138; // 187Pb //fBrho0 = 11.3418; // 216Th //fBrho0 = 11.2712; // 207Fr //fBrho0 = 10.6814; // 175Pt + //fBrho0 = 10.9970; // 199At } //////////////////////////////////////////////////////////////////////////////// diff --git a/Projects/s455/Analysis.h b/Projects/s455/Analysis.h index 747011bf95daf3f413cc443f5b4c5c2b114b4962..33244a7ad31e0cc96a218d26fd8027ba6d9a6f67 100644 --- a/Projects/s455/Analysis.h +++ b/Projects/s455/Analysis.h @@ -32,6 +32,7 @@ #include"TSofAtPhysics.h" #include"TSofTwimPhysics.h" #include"TSofSciPhysics.h" +#include"TSofMwpcPhysics.h" #include"TSofBeamID.h" #include"TSofFissionFragment.h" @@ -59,6 +60,7 @@ class Analysis: public NPL::VAnalysis{ private: TSofSciPhysics* SofSci; + TSofMwpcPhysics* SofMwpc; TSofTrimPhysics* SofTrim; TSofAtPhysics* SofAt; TSofTwimPhysics* SofTwim; @@ -74,6 +76,7 @@ class Analysis: public NPL::VAnalysis{ double fK_LS2; TRandom3 rand; + TCutG* cutQ77[3]; TCutG* cutQ78[3]; TCutG* cutQ79[3]; TCutG* cutQ80[3]; diff --git a/Projects/s455/configs/ConfigSofTrim.dat b/Projects/s455/configs/ConfigSofTrim.dat index 5815cd203cefd30d194adc398e65e50405a1ec52..caf5ae3a71bd4b56fdcdcc9fc59165cc076b8259 100644 --- a/Projects/s455/configs/ConfigSofTrim.dat +++ b/Projects/s455/configs/ConfigSofTrim.dat @@ -1,4 +1,5 @@ ConfigSofTrim SPLINE_PAIR_ANGLE_PATH ./calibration/SofTrim/spline/EvsA_spline.root SPLINE_PAIR_DT_PATH ./calibration/SofTrim/spline/EvsDT_spline.root -SPLINE_SECTION_DT_PATH ./calibration/SofTrim/spline/spline_section.root +SPLINE_SECTION_ANGLE_PATH ./calibration/SofTrim/spline/182Hg/spline_section_angle.root +SPLINE_SECTION_DT_PATH ./calibration/SofTrim/spline/182Hg/spline_section_dt.root diff --git a/Projects/s455/s455.detector b/Projects/s455/s455.detector index 58296af6c2815831c5eb2c688b231688d279941a..db54ef66d64def54964000304d63f988ba96e9c6 100644 --- a/Projects/s455/s455.detector +++ b/Projects/s455/s455.detector @@ -11,6 +11,9 @@ Target SofSci POS= 0 0 -136.3706933 m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +SofMwpc + POS= 7 -13 -883.5 mm +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SofSci POS= 0 0 0 m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -28,4 +31,11 @@ SofTofW THETA= -9.5 deg PHI= 0 deg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +SofMwpc + POS= 0 0 1.980 m +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +SofMwpc + POS= 0 0 2.576 m +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +SofMwpc + POS= 0 0 8 m