Merge branch 'NPTool.2.dev' of https://gitlab.in2p3.fr/np/nptool into NPTool.2.dev

NPLib/CMakeLists.txt

@@ -32,6 +32,9 @@ if(NPMULTITHREADING)
  message("Building application with no MutilThreading Support")
 endif()
 
+configure_file(scripts/build_dict.sh.in scripts/build_dict.sh @ONLY) 
+
+
 configure_file(Core/NPLibVersion.h.in Core/NPLibVersion.h @ONLY)
 set(DETLIST ${ETLIST})
 

NPLib/Core/NPDetectorManager.cxx

@@ -413,7 +413,7 @@ void NPL::DetectorManager::InitThreadPool(){
 ////////////////////////////////////////////////////////////////////////////////
 void NPL::DetectorManager::StartThread(NPL::VDetector* det,unsigned int id){ 
   // Let the main thread start 
-  std::this_thread::sleep_for(std::chrono::milliseconds(100));
+  std::this_thread::sleep_for(std::chrono::milliseconds(250));
   while(true){
     // Do the job if possible
     if(m_Ready[id]){

NPLib/Detectors/Minos/TMinosPhysics.h

@@ -160,6 +160,9 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
     double GetVertexX()  {return X_Vertex;} //!
     double GetVertexY()  {return Y_Vertex;} //!
     double GetVertexZ()  {return Z_Vertex;} //!
+    
+    double GetDeltaVertex()  {return Delta_Vertex;} //!
+    double GetTheta12()  {return Theta_12;} //!
   
     int GetNbrOfTracks(){return Tracks_P0.size();}
       

NPLib/Detectors/Nebula/TNebulaData.h

@@ -100,24 +100,44 @@ class TNebulaData : public TObject {
     fNebula_Td_ID.push_back(ID);
     fNebula_Td_Time.push_back(Time);
     };//!
-/*
-    //////////////////////    GETTERS    ////////////////////////
-    // Energy
-    inline UShort_t GetMultEnergy() const
-      {return fNebula_E_DetectorNbr.size();}
-    inline UShort_t GetE_DetectorNbr(const unsigned int &i) const 
-      {return fNebula_E_DetectorNbr[i];}//!
-    inline Double_t Get_Energy(const unsigned int &i) const 
-      {return fNebula_Energy[i];}//!
 
+    //////////////////////    GETTERS    ////////////////////////
+    // MULT //
+    // Charge 
+    inline unsigned int GetChargeUpMult() const
+      {return fNebula_Qu_ID.size();};
+    // Time
+    inline unsigned int GetTimeUpMult() const
+      {return fNebula_Tu_ID.size();};
+    // Charge
+    inline unsigned int GetChargeDownMult() const
+      {return fNebula_Qd_ID.size();};
     // Time
-    inline UShort_t GetMultTime() const
-      {return fNebula_T_DetectorNbr.size();}
-    inline UShort_t GetT_DetectorNbr(const unsigned int &i) const 
-      {return fNebula_T_DetectorNbr[i];}//!
-    inline Double_t Get_Time(const unsigned int &i) const 
-      {return fNebula_Time[i];}//!
-*/
+    inline unsigned int GetTimeDownMult() const
+      {return fNebula_Td_ID.size();};
+
+    // Value // 
+    // Charge 
+    inline UShort_t GetChargeUpID(unsigned int& i) const
+      {return fNebula_Qu_ID[i];};
+    inline double GetChargeUp(unsigned int& i) const
+      {return fNebula_Qu_Charge[i];};
+    // Time 
+    inline UShort_t GetTimeUpID(unsigned int& i) const
+      {return fNebula_Tu_ID[i];};
+    inline double GetTimeUp(unsigned int& i) const
+      {return fNebula_Tu_Time[i];};
+    // Charge 
+    inline UShort_t GetChargeDownID(unsigned int& i) const
+      {return fNebula_Qd_ID[i];};
+    inline double GetChargeDown(unsigned int& i) const
+      {return fNebula_Qd_Charge[i];};
+    // Time 
+    inline UShort_t GetTimeDownID(unsigned int& i) const
+      {return fNebula_Td_ID[i];};
+    inline double GetTimeDown(unsigned int& i) const
+      {return fNebula_Td_Time[i];};
+
 
   //////////////////////////////////////////////////////////////
   // Required for ROOT dictionnary

NPLib/Detectors/Nebula/TNebulaPhysics.cxx

@@ -35,28 +35,61 @@ using namespace std;
 #include "RootOutput.h"
 #include "NPDetectorFactory.h"
 #include "NPOptionManager.h"
-
+#include "NPSystemOfUnits.h"
 //   ROOT
 #include "TChain.h"
 
 ClassImp(TNebulaPhysics)
 
 
-///////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////////////////////////
 TNebulaPhysics::TNebulaPhysics()
-   : m_EventData(new TNebulaData),
-     m_PreTreatedData(new TNebulaData),
-     m_EventPhysics(this),
-     m_Spectra(0),
-     m_E_RAW_Threshold(0), // adc channels
-     m_E_Threshold(0),     // MeV
-     m_NumberOfBars(0) {
-}
+  : m_EventData(new TNebulaData),
+  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_NumberOfBars(0) {
+  }
 
 ///////////////////////////////////////////////////////////////////////////
 /// A usefull method to bundle all operation to add a detector
 void TNebulaPhysics::ReadXML(NPL::XmlParser xml){ 
-  m_NumberOfBars++;
+  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;;
 } 
 
 ///////////////////////////////////////////////////////////////////////////
@@ -69,115 +102,109 @@ void TNebulaPhysics::BuildSimplePhysicalEvent() {
 ///////////////////////////////////////////////////////////////////////////
 void TNebulaPhysics::BuildPhysicalEvent() {
   // apply thresholds and calibration
-  PreTreat();
-/*
-  // match energy and time together
-  unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
-  unsigned int mysizeT = m_PreTreatedData->GetMultTime();
-  for (UShort_t e = 0; e < mysizeE ; e++) {
-    for (UShort_t t = 0; t < mysizeT ; t++) {
-      if (m_PreTreatedData->GetE_DetectorNbr(e) == m_PreTreatedData->GetT_DetectorNbr(t)) {
-        DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
-        Energy.push_back(m_PreTreatedData->Get_Energy(e));
-        Time.push_back(m_PreTreatedData->Get_Time(t));
-      }
-    }
-  }*/
-}
-
-///////////////////////////////////////////////////////////////////////////
-void TNebulaPhysics::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();
-/*
-  // Energy
-  unsigned int mysize = m_EventData->GetMultEnergy();
-  for (UShort_t i = 0; i < mysize ; ++i) {
-    if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) {
-      Double_t Energy = Cal->ApplyCalibration("Nebula/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i));
-      if (Energy > m_E_Threshold) {
-        m_PreTreatedData->SetEnergy(m_EventData->GetE_DetectorNbr(i), Energy);
+  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 double threshold;
+  // loop on Qup
+  for (unsigned int qup = 0; qup < QUsize ; qup++) {
+
+    rawQup = m_EventData->GetChargeUp(qup);
+    rawTup=-1;
+    rawQdown=-1;
+    rawTdown=-1;
+    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;
+        } //if match ID 
+
+      }// Qdwown 
+
+      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){
+        calQup=aQu[ID]*(rawQup-bQu[ID]);
+        calQdown=aQd[ID]*(rawQdown-bQd[ID]);
+        // cal T 
+        calTup=aTu[ID]*rawTup+bTu[ID];
+        // slew correction
+        calTup -= slwTu[ID]/sqrt(rawQup-bQu[ID]);
+
+        // cal T
+        calTdown=aTd[ID]*rawTdown+bTd[ID];
+
+        // slew correction
+        calTdown -= slwTd[ID]/sqrt(rawQdown-bQd[ID]);
+
+        // average value of Up and Down
+        calQ=sqrt(calQup*calQdown); 
+        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));
+
+          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]);
+
+          if(ID<120)
+            IsVeto.push_back(0);
+          else
+            IsVeto.push_back(1);
+
+        }
       }
-    }
-  }
 
-  // Time 
-  mysize = m_EventData->GetMultTime();
-  for (UShort_t i = 0; i < mysize; ++i) {
-    Double_t Time= Cal->ApplyCalibration("Nebula/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
-    m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
-  }
-  */
+    }// if raw threshold up
+  } // Qup
 }
 
-
-
 ///////////////////////////////////////////////////////////////////////////
-void TNebulaPhysics::ReadAnalysisConfig() {
-  bool ReadingStatus = false;
-
-  // path to file
-  string FileName = "./configs/ConfigNebula.dat";
-
-  // open analysis config file
-  ifstream AnalysisConfigFile;
-  AnalysisConfigFile.open(FileName.c_str());
+void TNebulaPhysics::PreTreat() {
 
-  if (!AnalysisConfigFile.is_open()) {
-    cout << " No ConfigNebula.dat found: Default parameter loaded for Analayis " << FileName << endl;
-    return;
-  }
-  cout << " Loading user parameter for Analysis from ConfigNebula.dat " << endl;
-
-  // Save it in a TAsciiFile
-  TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
-  asciiConfig->AppendLine("%%% ConfigNebula.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 = "ConfigNebula";
-    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_RAW_THRESHOLD") {
-        AnalysisConfigFile >> DataBuffer;
-        m_E_RAW_Threshold = atof(DataBuffer.c_str());
-        cout << whatToDo << " " << m_E_RAW_Threshold << endl;
-      }
 
-      else if (whatToDo=="E_THRESHOLD") {
-        AnalysisConfigFile >> DataBuffer;
-        m_E_Threshold = atof(DataBuffer.c_str());
-        cout << whatToDo << " " << m_E_Threshold << endl;
-      }
 
-      else {
-        ReadingStatus = false;
-      }
-    }
-  }
+///////////////////////////////////////////////////////////////////////////
+void TNebulaPhysics::ReadAnalysisConfig() {
 }
 
 
@@ -185,8 +212,12 @@ void TNebulaPhysics::ReadAnalysisConfig() {
 ///////////////////////////////////////////////////////////////////////////
 void TNebulaPhysics::Clear() {
   DetectorNumber.clear();
-  Energy.clear();
-  Time.clear();
+  Charge.clear();
+  TOF.clear();
+  PosY.clear();
+  PosX.clear();
+  PosZ.clear();
+  IsVeto.clear();
 }
 
 
@@ -227,7 +258,6 @@ void TNebulaPhysics::InitSpectra() {
 ///////////////////////////////////////////////////////////////////////////
 void TNebulaPhysics::FillSpectra() {
   m_Spectra -> FillRawSpectra(m_EventData);
-  m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
   m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
 }
 
@@ -268,8 +298,8 @@ void TNebulaPhysics::WriteSpectra() {
 void TNebulaPhysics::AddParameterToCalibrationManager() {
   CalibrationManager* Cal = CalibrationManager::getInstance();
   for (int i = 0; i < m_NumberOfBars; ++i) {
-    Cal->AddParameter("NEBULA_ID"+ NPL::itoa(i+1)+"_T");
-    Cal->AddParameter("NEBULA_ID"+ NPL::itoa(i+1)+"_Y");
+    Cal->AddParameter("NEBULA_T_ID"+ NPL::itoa(i+1));
+    Cal->AddParameter("NEBULA_Y_ID"+ NPL::itoa(i+1));
   }
 }
 
@@ -313,14 +343,14 @@ NPL::VDetector* TNebulaPhysics::Construct() {
 //            Registering the construct method to the factory                 //
 ////////////////////////////////////////////////////////////////////////////////
 extern "C"{
-class proxy_Nebula{
-  public:
-    proxy_Nebula(){
-      NPL::DetectorFactory::getInstance()->AddToken("NEBULA","Nebula");
-      NPL::DetectorFactory::getInstance()->AddDetector("NEBULA",TNebulaPhysics::Construct);
-    }
-};
+  class proxy_Nebula{
+    public:
+      proxy_Nebula(){
+        NPL::DetectorFactory::getInstance()->AddToken("NEBULA","Nebula");
+        NPL::DetectorFactory::getInstance()->AddDetector("NEBULA",TNebulaPhysics::Construct);
+      }
+  };
 
-proxy_Nebula p_Nebula;
+  proxy_Nebula p_Nebula;
 }
 

NPLib/Detectors/Nebula/TNebulaPhysics.h

@@ -52,26 +52,64 @@ class TNebulaPhysics : public TObject, public NPL::VDetector {
     ~TNebulaPhysics() {};
 
 
-  //////////////////////////////////////////////////////////////
-  // Inherited from TObject and overriden to avoid warnings
+    //////////////////////////////////////////////////////////////
+    // 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
+    //////////////////////////////////////////////////////////////
+    // data obtained after BuildPhysicalEvent() and stored in
+    // output ROOT file
   public:
     vector<int>      DetectorNumber;
-    vector<double>   Energy;
-    vector<double>   Time;
+    vector<double>   Charge;
+    vector<double>   TOF;
+    vector<double>   PosY;
+    vector<double>   PosX;
+    vector<double>   PosZ;
+    vector<bool>     IsVeto;
 
-  /// A usefull method to bundle all operation to add a detector
-  void ReadXML(NPL::XmlParser); 
-  
-  //////////////////////////////////////////////////////////////
-  // methods inherited from the VDetector ABC class
+  public:
+    TVector3 GetPos(const unsigned int& i) const{
+      return TVector3(PosX[i],PosY[i],PosZ[i]);
+    }
+
+    // Return true if one veto fired
+    bool HasVeto(){
+      unsigned int size = IsVeto.size();
+      for(unsigned int i = 0 ; i < size ; i++){
+        if(IsVeto[i])
+          return true;
+      }
+      return false;
+    };
+
+    /////////// Get index of fastest neutron
+    int GetFirstHit(){
+      unsigned int size = TOF.size();
+      unsigned int index=0;
+
+      if(!size)
+        return -1;
+
+      double tof = TOF[0];
+      for(unsigned int i = 1 ; i < size ; i++){
+        if(tof<TOF[i]){
+          tof=TOF[i];
+          index=i;
+        }
+      }
+      return index;
+    };
+
+  public:
+    /// A usefull method to bundle all operation to add a detector
+    void ReadXML(NPL::XmlParser); 
+
+    //////////////////////////////////////////////////////////////
+    // methods inherited from the VDetector ABC class
   public:
     // read stream from ConfigFile to pick-up detector parameters
     void ReadConfiguration(NPL::InputParser);
@@ -126,40 +164,36 @@ class TNebulaPhysics : public TObject, public NPL::VDetector {
     void WriteSpectra();
 
 
-  //////////////////////////////////////////////////////////////
-  // specific methods to Nebula array
+    //////////////////////////////////////////////////////////////
+    // specific methods to Nebula 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 TNebulaData object to TNebulaPhysics. 
     // needed for online analysis for example
     void SetRawDataPointer(TNebulaData* rawDataPointer) {m_EventData = rawDataPointer;}
-    
-  // objects are not written in the TTree
+
+    // objects are not written in the TTree
   private:
     TNebulaData*         m_EventData;        //!
-    TNebulaData*         m_PreTreatedData;   //!
     TNebulaPhysics*      m_EventPhysics;     //!
 
-  // getters for raw and pre-treated data object
+    // getters for raw and pre-treated data object
   public:
     TNebulaData* GetRawData()        const {return m_EventData;}
-    TNebulaData* GetPreTreatedData() const {return m_PreTreatedData;}
 
-  // parameters used in the analysis
+    // parameters used in the analysis
   private:
     // thresholds
-    double m_E_RAW_Threshold; //!
-    double m_E_Threshold;     //!
+    double m_Q_RAW_Threshold; //!
+    double m_Q_Threshold;     //!
+    double m_V_Threshold;     //!
 
-  // number of detectors
+    // number of detectors
   private:
     int m_NumberOfBars;  //!
 
@@ -168,16 +202,43 @@ class TNebulaPhysics : public TObject, public NPL::VDetector {
     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
+    // spectra class
   private:
     TNebulaSpectra* m_Spectra; // !
 
-  // spectra getter
+    // spectra getter
   public:
     map<string, TH1*>   GetSpectra();