diff --git a/NPLib/Detectors/NebulaPlus/TNebulaPlusData.h b/NPLib/Detectors/NebulaPlus/TNebulaPlusData.h index f11345e2a2531945ec79af0f040bb643c26fe7e3..4ec8fef6c467904c8426ca601eb0289bd22f4945 100644 --- a/NPLib/Detectors/NebulaPlus/TNebulaPlusData.h +++ b/NPLib/Detectors/NebulaPlus/TNebulaPlusData.h @@ -89,30 +89,30 @@ class TNebulaPlusData : public TObject { ////////////////////// GETTERS //////////////////////// // MULT // - inline unsigned int GetUpMult() const + inline unsigned int GetMultUp() const {return fNebulaPlus_u_ID.size();}; - inline unsigned int GetDownMult() const + inline unsigned int GetMultDown() const {return fNebulaPlus_d_ID.size();}; // Value // // Up - inline UShort_t GetUpID(unsigned int& i) const + inline unsigned short GetIDUp(const int& i) const {return fNebulaPlus_u_ID[i];}; - inline double GetUpQ(unsigned int& i) const + inline double GetQUp(const int& i) const {return fNebulaPlus_u_Q[i];}; - inline double GetUpQ4(unsigned int& i) const + inline double GetQ4Up(const int& i) const {return fNebulaPlus_u_Q4[i];}; - inline double GetUpT(unsigned int& i) const + inline double GetTUp(const int& i) const {return fNebulaPlus_u_T[i];}; // Down - inline UShort_t GetDownID(unsigned int& i) const + inline unsigned short GetIDDown(const int& i) const {return fNebulaPlus_d_ID[i];}; - inline double GetDownQ(unsigned int& i) const + inline double GetQDown(const int& i) const {return fNebulaPlus_d_Q[i];}; - inline double GetDownQ4(unsigned int& i) const + inline double GetQ4Down(const int& i) const {return fNebulaPlus_d_Q4[i];}; - inline double GetDownT(unsigned int& i) const + inline double GetTDown(const int& i) const {return fNebulaPlus_d_T[i];}; ////////////////////////////////////////////////////////////// diff --git a/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.cxx b/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.cxx index 4c7dcaacf44635d4104c34260caf8c8e8a9164df..2cc21b541427acf6d5a3ad51dc1be33055f901a6 100644 --- a/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.cxx +++ b/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.cxx @@ -6,7 +6,7 @@ *****************************************************************************/ /***************************************************************************** - * Original Author: Adrien Matta contact address: matta@lpccaen.in2p3.fr * + * Original Author: Freddy Flavigny contact: flavigny@lpccaen.in2p3.fr * * * * Creation Date : December 2019 * * Last update : * @@ -21,6 +21,7 @@ *****************************************************************************/ #include "TNebulaPlusPhysics.h" +using namespace NEBULAPLUS_LOCAL; // STL #include <sstream> @@ -47,50 +48,56 @@ TNebulaPlusPhysics::TNebulaPlusPhysics() : m_EventData(new TNebulaPlusData), m_EventPhysics(this), m_Spectra(0), - m_Q_RAW_Threshold(0), // adc channels - m_Q_Threshold(7), // normal bars in MeV - m_V_Threshold(1), // veto bars in MeV + m_Q_RAW_Threshold(15000), // adc channels + m_Q_Threshold(0), // normal bars in MeV + m_V_Threshold(0), // veto bar in MeV + m_TotalNbrModules(0), + m_BarWidth(12), // bar width in cm + m_BarLength(180), // bar length in cm + m_TdiffRange(12), // +- range in (Td-Tu) in ns m_NumberOfBars(0) { } -/////////////////////////////////////////////////////////////////////////// -/// A usefull method to bundle all operation to add a detector -void TNebulaPlusPhysics::ReadXML(NPL::XmlParser xml){ - std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName("NEBULA"); - - for(unsigned int i = 0 ; i < b.size() ; i++){ - m_NumberOfBars++; - unsigned int id = b[i]->AsInt("ID"); - - // position - PositionX[id] = b[i]->AsDouble("PosX"); - PositionY[id] = b[i]->AsDouble("PosY"); - PositionZ[id] = b[i]->AsDouble("PosZ"); - // linear cal - aQu[id] = b[i]->AsDouble("QUCal"); - bQu[id] = b[i]->AsDouble("QUPed"); - aQd[id] = b[i]->AsDouble("QDCal"); - bQd[id] = b[i]->AsDouble("QDPed"); - aTu[id] = b[i]->AsDouble("TUCal"); - bTu[id] = b[i]->AsDouble("TUOff"); - aTd[id] = b[i]->AsDouble("TDCal"); - bTd[id] = b[i]->AsDouble("TDOff"); - - // T average offset - avgT0[id] = b[i]->AsDouble("TAveOff"); - - // slew correction T= tcal +slwT/sqrt(Qcal) - slwTu[id] = b[i]->AsDouble("TUSlw"); - slwTd[id] = b[i]->AsDouble("TDSlw"); - - // DT position cal - DTa[id] = b[i]->AsDouble("DTCal");//! - DTb[id] = b[i]->AsDouble("DTOff");//! - - - } - cout << " -> " << m_NumberOfBars << " bars found" << endl;; -} +///////////////////////////////////////////// +//// Innherited from VDetector Class //// +///////////////////////////////////////////// +void TNebulaPlusPhysics::ReadConfiguration(NPL::InputParser parser) { + + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("NEBULAPLUS"); + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " Nebula layers found " << endl; + + unsigned int det=0; + // Cartesian Case + vector<string> info + = {"POS","NumberOfModules", "Veto", "Frame"}; + string Type; + + for (unsigned int i = 0; i < blocks.size(); i++) { + + if (blocks[i]->HasTokenList(info)) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + + cout << endl << "//// Nebula Layer " << Type << " " << i + 1 << endl; + int NbrOfModulesInLayer = blocks[i]->GetInt("NumberOfModules"); + //det = i+1; + //m_DetectorNumberIndex[detectorNbr]=det; + TVector3 Pos = blocks[i]->GetTVector3("Pos", "cm"); + AddLayer(Pos, NbrOfModulesInLayer); + } + + else { + cout << "ERROR: Missing token for NebulaPlus, check your input " + "file" + << endl; + exit(1); + } + + } + + //InitializeStandardParameter(); + //ReadAnalysisConfig(); +} /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::BuildSimplePhysicalEvent() { @@ -101,103 +108,101 @@ void TNebulaPlusPhysics::BuildSimplePhysicalEvent() { /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::BuildPhysicalEvent() { + // apply thresholds and calibration -/* + //std::cout << "----------------" << std::endl; + //std::cout << "in BuildPhysical" << std::endl; + // instantiate CalibrationManager static CalibrationManager* Cal = CalibrationManager::getInstance(); + //Cal->Dump(); + //double offset = Cal->GetValue("NebulaPlus/Tdiff_Offset_58",0); + //double offset = Cal->GetValue("NebulaPlus/U3_ENERGY",1); + + static double rawQup,calQup,rawQdown,calQdown,rawTup,calTup,rawTdown,calTdown,calQ,calT,Y; static unsigned int ID; // All vector size - static unsigned int QUsize, QDsize, TUsize, TDsize ; - QUsize = m_EventData->GetChargeUpMult(); - QDsize = m_EventData->GetChargeDownMult(); - TUsize = m_EventData->GetTimeUpMult(); - TDsize = m_EventData->GetTimeDownMult(); + static unsigned int Usize, Dsize ; + Usize = m_EventData->GetMultUp(); + Dsize = m_EventData->GetMultDown(); static double threshold; - // loop on Qup - for (unsigned int qup = 0; qup < QUsize ; qup++) { + bool matchUD; + int unmatched=0; + + // loop on up + //std::cout << "size up" << Usize<< std::endl; + //std::cout << "size down" << Dsize<< std::endl; - rawQup = m_EventData->GetChargeUp(qup); - rawTup=-1; + for (unsigned int up = 0; up < Usize ; up++) { + + rawQup = m_EventData->GetQUp(up); + rawTup= m_EventData->GetTUp(up); rawQdown=-1; rawTdown=-1; + matchUD = false; if (rawQup > m_Q_RAW_Threshold) { - ID = m_EventData->GetChargeUpID(qup); - if(ID<121) - threshold=m_Q_Threshold; - else - threshold=m_V_Threshold; - - // look for associated Charge down - for(unsigned int qdown = 0 ; qdown < QDsize ; qdown++){ - if(m_EventData->GetChargeDownID(qdown)==ID){ - rawQdown=m_EventData->GetChargeDown(qdown); - if(rawQdown > m_Q_RAW_Threshold){ - // Look for the associate time - for(unsigned int tdown = 0 ; tdown < TDsize; tdown++){ - if(m_EventData->GetTimeDownID(qdown)==ID) { - rawTdown=m_EventData->GetTimeDown(qdown); - break; - } - }// TDown - }//if raw threshold down - break; + ID = m_EventData->GetIDUp(up); + + // look for associated down + for(unsigned int down = 0 ; down < Dsize ; down++){ + if(m_EventData->GetIDDown(down)==ID){ + + //std::cout << "-------" << std::endl; + //std::cout << "same ID:\t" << ID<< std::endl; + if( m_EventData->GetQDown(down) > m_Q_RAW_Threshold){ + rawQdown = m_EventData->GetQDown(down); + rawTdown = m_EventData->GetTDown(down); + calQup = fCalQUp(m_EventData,up); + calQdown = fCalQDown(m_EventData,down); + matchUD = true; + break; + }//if raw threshold down } //if match ID + }// down - }// Qdwown + // Got a Down-Up match, do the math + if(matchUD){ - if(rawTdown>0){ // Tdown is found, means Qdown as well - // look for Tup - for(unsigned int tup = 0 ; tup < TUsize ; tup++){ - if(m_EventData->GetTimeUpID(tup)==ID){ - rawTup = m_EventData->GetTimeUp(tup); - break; - } - } - } - // Got everything, do the math - if(rawTup>0){ // cal Q Up and Down - calQup=aQu[ID]*(rawQup-bQu[ID]); - calQdown=aQd[ID]*(rawQdown-bQd[ID]); + //calQup = aQu[ID] * (rawQup-bQu[ID]); + //calQdown= aQd[ID] * (rawQdown-bQd[ID]); // average value of Up and Down calQ=sqrt(calQup*calQdown); // cal T Up - calTup=aTu[ID]*rawTup+bTu[ID]; - // slew correction - calTup -= slwTu[ID]/sqrt(rawQup-bQu[ID]); + //calTup=aTu[ID]*rawTup+bTu[ID]; + calTup=rawTup; // cal T Down - calTdown=aTd[ID]*rawTdown+bTd[ID]; - // slew correction - calTdown -= slwTd[ID]/sqrt(rawQdown-bQd[ID]); + //calTdown=aTd[ID]*rawTdown+bTd[ID]; + calTdown=rawTdown; - if(calQ>threshold){ - calT= (calTdown+calTup)*0.5+avgT0[ID]+Cal->GetPedestal("NEBULA_T_ID"+NPL::itoa(ID)); - Y=(calTdown-calTup)*DTa[ID]+DTb[ID]+Cal->GetPedestal("NEBULA_Y_ID"+NPL::itoa(ID)); + calT= (calTdown+calTup)*0.5; + Y=(calTdown-calTup) - Cal->GetValue("NebulaPlus/Tdiff_Offset_"+NPL::itoa(ID),0); DetectorNumber.push_back(ID); Charge.push_back(calQ); TOF.push_back(calT); - PosY.push_back(Y+PositionY[ID]); - PosX.push_back(PositionX[ID]); - PosZ.push_back(PositionZ[ID]); - + PosX.push_back(m_BarPositionX[ID-1]); + PosY.push_back(Y/m_TdiffRange*m_BarLength); + PosZ.push_back(m_BarPositionZ[ID-1]); +/* if(ID<121) IsVeto.push_back(0); else IsVeto.push_back(1); - +*/ } + } else { + unmatched++;/*std::cout << "unmatched=" <<unmatched<< std::endl;*/ } - }// if raw threshold up } // Qup -*/ + } /////////////////////////////////////////////////////////////////////////// @@ -205,8 +210,17 @@ void TNebulaPlusPhysics::PreTreat() { } +/////////////////////////////////////////////////////////////////////////// +void TNebulaPlusPhysics::AddLayer(TVector3 Pos, int n_modules) { + for(int i=0; i<n_modules; i++){ + m_BarPositionX.push_back(Pos.X() + i * m_BarWidth); + m_BarPositionY.push_back(Pos.Y()); + m_BarPositionZ.push_back(Pos.Z()); + m_TotalNbrModules++; + } +} /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::ReadAnalysisConfig() { } @@ -227,8 +241,9 @@ void TNebulaPlusPhysics::Clear() { /////////////////////////////////////////////////////////////////////////// +/* void TNebulaPlusPhysics::ReadConfiguration(NPL::InputParser parser) { - vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("NEBULA"); + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("NEBULAPLUS"); if(NPOptionManager::getInstance()->GetVerboseLevel()) cout << "//// " << blocks.size() << " detector(s) found " << endl; @@ -251,7 +266,7 @@ void TNebulaPlusPhysics::ReadConfiguration(NPL::InputParser parser) { } } } - +*/ /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::InitSpectra() { m_Spectra = new TNebulaPlusSpectra(m_NumberOfBars); @@ -296,21 +311,18 @@ void TNebulaPlusPhysics::WriteSpectra() { m_Spectra->WriteSpectra(); } - - /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::AddParameterToCalibrationManager() { CalibrationManager* Cal = CalibrationManager::getInstance(); - vector<double> standardO={0}; - for (int i = 0; i < m_NumberOfBars; ++i) { - Cal->AddParameter("NEBULA_T_ID"+ NPL::itoa(i+1),standardO); - Cal->AddParameter("NEBULA_Y_ID"+ NPL::itoa(i+1),standardO); + vector<double> standard={8192, 6}; + for (int i = 0; i < m_TotalNbrModules; ++i) { + Cal->AddParameter("NebulaPlus","U"+ NPL::itoa(i+1)+"_ENERGY","NebulaPlus_U"+ NPL::itoa(i+1)+"_ENERGY",standard); + Cal->AddParameter("NebulaPlus","D"+ NPL::itoa(i+1)+"_ENERGY","NebulaPlus_D"+ NPL::itoa(i+1)+"_ENERGY",standard); + Cal->AddParameter("NebulaPlus","Tdiff_Offset_"+ NPL::itoa(i+1),"NebulaPlus_Tdiff_Offset_"+ NPL::itoa(i+1),standard); } } - - /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::InitializeRootInputRaw() { TChain* inputChain = RootInput::getInstance()->GetChain(); @@ -318,8 +330,6 @@ void TNebulaPlusPhysics::InitializeRootInputRaw() { inputChain->SetBranchAddress("NebulaPlus", &m_EventData ); } - - /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::InitializeRootInputPhysics() { TChain* inputChain = RootInput::getInstance()->GetChain(); @@ -327,14 +337,32 @@ void TNebulaPlusPhysics::InitializeRootInputPhysics() { } - /////////////////////////////////////////////////////////////////////////// void TNebulaPlusPhysics::InitializeRootOutput() { TTree* outputTree = RootOutput::getInstance()->GetTree("NebulaPlus"); outputTree->Branch("NebulaPlus", "TNebulaPlusPhysics", &m_EventPhysics); } +//////////////////////////////////////////////////////////////////////////////// +namespace NEBULAPLUS_LOCAL{ + +double fCalQUp(const TNebulaPlusData* m_EventData, const int& i){ + static string name; + name = "NebulaPlus/U"; + name += NPL::itoa(m_EventData->GetIDUp(i)); + name += "_ENERGY"; + return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetQUp(i), 1); +} + +double fCalQDown(const TNebulaPlusData* m_EventData, const int& i){ + static string name; + name = "NebulaPlus/D"; + name += NPL::itoa(m_EventData->GetIDDown(i)); + name += "_ENERGY"; + return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetQDown(i), 1); +} +} //////////////////////////////////////////////////////////////////////////////// // Construct Method to be pass to the DetectorFactory // diff --git a/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.h b/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.h index 32b1a9547f03740eca22cbcd4898283be3e732d4..4d4d6c55ce10e848e8021f48cb83faaf14334a7d 100644 --- a/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.h +++ b/NPLib/Detectors/NebulaPlus/TNebulaPlusPhysics.h @@ -76,7 +76,7 @@ class TNebulaPlusPhysics : public TObject, public NPL::VDetector { return TVector3(PosX[i],PosY[i],PosZ[i]); } - // Return true if one veto fired + // Return true if one veto fired (not used for now) bool HasVeto(){ unsigned int size = IsVeto.size(); for(unsigned int i = 0 ; i < size ; i++){ @@ -105,8 +105,8 @@ class TNebulaPlusPhysics : public TObject, public NPL::VDetector { }; public: - /// A usefull method to bundle all operation to add a detector - void ReadXML(NPL::XmlParser); + /// A useful method to bundle all operation to add a detector + void ReadCalibration(NPL::InputParser); ////////////////////////////////////////////////////////////// // methods inherited from the VDetector ABC class @@ -193,46 +193,33 @@ class TNebulaPlusPhysics : public TObject, public NPL::VDetector { double m_Q_Threshold; //! double m_V_Threshold; //! + private: //FF addition + int m_TotalNbrModules; //! + double m_BarWidth; //! + double m_BarLength; //! + double m_TdiffRange; //! + vector<double> m_BarPositionX; //! + vector<double> m_BarPositionY; //! + vector<double> m_BarPositionZ; //! + void AddLayer(TVector3 Pos, int n_modules); + + public: + double GetBarPositionX(const int N) {return m_BarPositionX[N - 1];} + double GetBarPositionY(const int N) {return m_BarPositionY[N - 1];} + double GetBarPositionZ(const int N) {return m_BarPositionZ[N - 1];} + public: void SetQThreshold(double t) {m_Q_Threshold=t;}; void SetVThreshold(double t) {m_V_Threshold=t;}; // number of detectors private: int m_NumberOfBars; //! - +/* private: // offset and inversion std::map<unsigned int, TVector3> m_offset;//! std::map<unsigned int, bool> m_invertX;//! std::map<unsigned int, bool> m_invertY;//! - - private: // xml calibration - // position - std::map<unsigned int , double > PositionX;//! - std::map<unsigned int , double > PositionY;//! - std::map<unsigned int , double > PositionZ;//! - - // linear cal - std::map<unsigned int , double > aQu;//! - std::map<unsigned int , double > bQu;//! - std::map<unsigned int , double > aQd;//! - std::map<unsigned int , double > bQd;//! - std::map<unsigned int , double > aTu;//! - std::map<unsigned int , double > bTu;//! - std::map<unsigned int , double > aTd;//! - std::map<unsigned int , double > bTd;//! - - // T average offset - std::map<unsigned int , double > avgT0;//! - - // slew correction T= tcal +slwT/sqrt(Qcal) - std::map<unsigned int , double > slwTu;//! - std::map<unsigned int , double > slwTd;//! - - // DT position cal - std::map<unsigned int , double > DTa;//! - std::map<unsigned int , double > DTb;//! - - +*/ // spectra class private: TNebulaPlusSpectra* m_Spectra; // ! @@ -247,4 +234,12 @@ class TNebulaPlusPhysics : public TObject, public NPL::VDetector { ClassDef(TNebulaPlusPhysics,1) // NebulaPlusPhysics structure }; + +namespace NEBULAPLUS_LOCAL{ +double fCalQUp(const TNebulaPlusData* Data, const int& i); +double fCalQDown(const TNebulaPlusData* Data, const int& i); + + + +} #endif