* adding file for e748 analysis

NPLib/Core/CMakeLists.txt

@@ -2,6 +2,6 @@ add_custom_command(OUTPUT TAsciiFileDict.cxx COMMAND  ${CMAKE_BINARY_DIR}/script
 add_custom_command(OUTPUT NPDeltaSpectraDict.cxx COMMAND  ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPDeltaSpectra.h NPDeltaSpectraDict.cxx NPDeltaSpectra.rootmap libNPCore.so NPDeltaSpectraLinkdef.h)  
 
 add_custom_command(OUTPUT NPVDetectorDict.cxx COMMAND  ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPVDetector.h NPVDetectorDict.cxx NPVDetector.rootmap libNPCore.so NPCoreLinkdef.h)  
-add_library(NPCore SHARED NPCore.cxx NPRandom.cxx NPVAnalysis.cxx NPAnalysisFactory.cxx NPCalibrationManager.cxx NPOptionManager.cxx RootOutput.cxx RootInput.cxx TAsciiFile.cxx TAsciiFileDict.cxx NPDeltaSpectraDict.cxx NPDetectorManager.cxx NPVDetector.cxx NPVDetectorDict.cxx NPVSpectra.cxx NPDetectorFactory.cxx NPSpectraServer.cxx NPInputParser.cxx NPImage.cxx NPElog.cxx NPDeltaSpectra.cxx)
-target_link_libraries(NPCore ${ROOT_LIBRARIES}) 
-install(FILES NPCore.h NPVAnalysis.h NPAnalysisFactory.h NPRandom.h NPCalibrationManager.h NPOptionManager.h RootInput.h RootOutput.h TAsciiFile.h NPDetectorManager.h NPVDetector.h NPGlobalSystemOfUnits.h NPPhysicalConstants.h NPSystemOfUnits.h NPVSpectra.h NPDetectorFactory.h  NPSpectraServer.h NPInputParser.h NPImage.h NPElog.h NPDeltaSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+add_library(NPCore SHARED NPCore.cxx NPRandom.cxx NPVAnalysis.cxx NPAnalysisFactory.cxx NPCalibrationManager.cxx NPOptionManager.cxx RootOutput.cxx RootInput.cxx TAsciiFile.cxx TAsciiFileDict.cxx NPDeltaSpectraDict.cxx NPDetectorManager.cxx NPVDetector.cxx NPVDetectorDict.cxx NPVSpectra.cxx NPDetectorFactory.cxx NPSpectraServer.cxx NPInputParser.cxx NPImage.cxx NPElog.cxx NPDeltaSpectra.cxx NPXmlParser.cxx)
+target_link_libraries(NPCore ${ROOT_LIBRARIES} -lXMLIO) 
+install(FILES NPCore.h NPVAnalysis.h NPAnalysisFactory.h NPRandom.h NPCalibrationManager.h NPOptionManager.h RootInput.h RootOutput.h TAsciiFile.h NPDetectorManager.h NPVDetector.h NPGlobalSystemOfUnits.h NPPhysicalConstants.h NPSystemOfUnits.h NPVSpectra.h NPDetectorFactory.h  NPSpectraServer.h NPInputParser.h NPImage.h NPElog.h NPDeltaSpectra.h NPXmlParser.cxx DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})

NPLib/Core/NPXmlParser.cxx

+#include "NPXmlParser.h"
+#include <stdio.h>
+#include <iostream>
+using namespace NPL;
+
+////////////////////////////////////////////////////////////////////////////////
+block::block(){};
+block::~block(){};
+int block::AsInt(std::string name){
+  parameter p(name,"void");  
+  auto it = m_parameters.find(p);
+  if(it!=m_parameters.end()){
+    std::string s = it->GetValue(); 
+    return atoi(s.c_str());
+    }
+  else
+    return -1000;
+}
+
+double block::AsDouble(std::string name){
+  parameter p(name,"void");  
+  auto it = m_parameters.find(p);
+  if(it!=m_parameters.end()){
+    std::string s = it->GetValue(); 
+    return atof(s.c_str());
+    }
+  else
+    return -1000;
+}
+std::string block::AsString(std::string name){
+  parameter p(name,"void");  
+  auto it = m_parameters.find(p);
+  if(it!=m_parameters.end()){
+    std::string s = it->GetValue(); 
+    return s;
+    }
+  else
+    return "void";
+}
+ 
+
+////////////////////////////////////////////////////////////////////////////////
+parameter::parameter(){};
+parameter::~parameter(){};
+
+////////////////////////////////////////////////////////////////////////////////
+Channel::Channel(){};
+Channel::~Channel(){};
+
+////////////////////////////////////////////////////////////////////////////////
+void XmlParser::LoadFile(std::string file){
+     // First create engine
+   TXMLEngine* xml = new TXMLEngine;
+   // Now try to parse xml file
+   // Only file with restricted xml syntax are supported
+   XMLDocPointer_t xmldoc = xml->ParseFile(file.c_str());
+   if (xmldoc==0) {
+      delete xml;
+      return;
+   }
+   // take access to main node
+   XMLNodePointer_t mainnode = xml->DocGetRootElement(xmldoc);
+   // display recursively all nodes and subnodes
+   LoadNode(xml, mainnode, 1);
+
+   // Release memory before exit
+   xml->FreeDoc(xmldoc);
+   delete xml;
+  }
+
+//////////////////////////////////////////////////////////////////////////////////
+void XmlParser::DisplayNode(TXMLEngine* xml, XMLNodePointer_t node, Int_t level){
+   // this function display all accessible information about xml node and its children
+   printf("%*c node: %s\n",level,' ', xml->GetNodeName(node));
+   // display namespace
+   XMLNsPointer_t ns = xml->GetNS(node);
+   if (ns!=0)
+      printf("%*c namespace: %s refer: %s\n",level+2,' ', xml->GetNSName(ns), xml->GetNSReference(ns));
+   // display attributes
+   XMLAttrPointer_t attr = xml->GetFirstAttr(node);
+   while (attr!=0) {
+       printf("%*c attr: %s value: %s\n",level+2,' ', xml->GetAttrName(attr), xml->GetAttrValue(attr));
+       attr = xml->GetNextAttr(attr);
+   }
+   // display content (if exists)
+   const char* content = xml->GetNodeContent(node);
+   if (content!=0)
+      printf("%*c cont: %s\n",level+2,' ', content);
+   // display all child nodes
+   XMLNodePointer_t child = xml->GetChild(node);
+   while (child!=0) {
+      DisplayNode(xml, child, level+2);
+      child = xml->GetNext(child);
+   }
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+void XmlParser::LoadNode(TXMLEngine* xml, XMLNodePointer_t node, Int_t level){
+   // namespace
+   XMLNsPointer_t ns = xml->GetNS(node);
+   XMLNodePointer_t child = xml->GetChild(node);
+   if(xml->GetNodeName(child)=="dataroot"){// main node
+      std::cout <<" Loading XML file" << std::endl;
+     }
+   else{
+    while(child!=0) {
+      block b;
+      // getting attribute:
+      XMLNodePointer_t param = xml->GetChild(child);
+      while (param!=0) {
+        parameter p(xml->GetNodeName(param),xml->GetNodeContent(param));
+        b.AddParameter(p); 
+        param=xml->GetNext(param);
+      }
+      std::string name = xml->GetNodeName(child);
+      m_blocks[name].push_back(b);
+      child = xml->GetNext(child);
+    }  
+    std::cout << " XML file loaded for " <<m_blocks.size() << " detectors" << std::endl;
+  }
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+std::vector<NPL::block*> XmlParser::GetAllBlocksWithName(std::string name){
+   std::vector<NPL::block*> res;
+   auto it=m_blocks.find(name);
+
+    if(it!=m_blocks.end()){
+      unsigned int size = it->second.size();
+      for(unsigned int i = 0 ; i < size ; i++){
+          res.push_back(&(it->second[i])); 
+        }
+      }
+  
+   return res;
+  }
+/////////////////////////////////////////////////////////////////////////////////
+std::vector<std::string> XmlParser::GetAllBlocksName(){
+  std::vector<std::string> res;
+ for(auto it=m_blocks.begin(); it!= m_blocks.end();++it){
+     res.push_back(it->first);
+   } 
+   return res;
+  }

NPLib/Core/NPXmlParser.h

+#ifndef NPXMLPARSER_H
+#define NPXMLPARSER_H
+#include<string>
+#include<map>
+#include<set>
+#include<vector>
+#include "TXMLEngine.h"
+namespace NPL{
+  /////////////////////
+  class parameter{
+    public:
+      parameter();
+      parameter(std::string name, std::string value){
+        m_name=name;
+        m_value=value;
+        };
+
+     ~parameter();
+
+    private:
+      std::string m_name;
+      std::string m_value;
+    
+    public:
+      std::string GetName()  const {return m_name;};
+      std::string GetValue() const {return m_value;};
+
+    public:
+    bool operator<(const parameter p2){
+      return this->m_name<p2.m_name;
+      }
+    friend bool operator<(const parameter p1,const parameter p2){
+      return p1.m_name<p2.m_name;
+      }
+    friend bool operator==(const parameter p1,const parameter p2){
+      return p1.m_name==p2.m_name;
+      };
+  };
+
+  /////////////////////
+  class Channel{
+    public:
+      Channel();
+      ~Channel();
+      Channel(int device,int geo,int ch){
+        m_device=device;
+        m_geo=geo;
+        m_ch=ch;
+       };
+
+    private:
+      int m_device;
+      int m_geo;
+      int m_ch;
+    
+    public:
+      int norme() const {return (m_device*1000000+m_geo*1000+m_ch);} ;
+
+    public:
+      bool operator<(const Channel p2){
+        return this->norme()<p2.norme();
+        }
+      
+      friend bool operator<(const Channel p1,const Channel p2){
+        return p1.norme()<p2.norme();
+        }
+      
+      friend bool operator==(const Channel p1,const Channel p2){
+       return p1.norme()==p2.norme();
+       }
+  };
+  /////////////////////
+  class block{
+    public:   
+      block();
+      ~block();
+    
+    public:
+      int AsInt(std::string name);
+      double AsDouble(std::string name);
+      std::string AsString(std::string name); 
+      void AddParameter(parameter p){ m_parameters.insert(p);};
+
+    private:
+      std::string m_name;
+      std::set<parameter> m_parameters;
+  };
+  
+  /////////////////////
+  class XmlParser{
+    public:
+      XmlParser(){};
+      ~XmlParser(){};
+
+    public:
+      void CheckFile(std::string file, std::string DTD);
+      void LoadFile(std::string file);
+      void DisplayNode(TXMLEngine* xml, XMLNodePointer_t node, Int_t level);
+      void LoadNode(TXMLEngine* xml, XMLNodePointer_t node, Int_t level);
+    
+    public: // access by channel id
+      block* GetBlock(Channel c){
+        auto it = m_Channels.find(c);
+        if(it!=m_Channels.end())
+          return it->second;
+        else
+          return NULL;
+        };
+      void SetBlock(const Channel& c,block* b){
+        m_Channels[c]=b; 
+        } 
+ 
+    public:
+      std::vector<block*> GetAllBlocksWithName(std::string name) ;
+      std::vector<std::string> GetAllBlocksName();
+    private:
+      std::map<std::string,std::vector<block>> m_blocks;
+      std::map<Channel,block*> m_Channels;
+      
+  };
+
+}
+
+#endif

NPLib/Detectors/Samurai/TSamuraiFDC2Physics.cxx

@@ -33,8 +33,10 @@
 #include "TAsciiFile.h"
 #include "NPOptionManager.h"
 #include "NPDetectorFactory.h"
+#include "NPSystemOfUnits.h"
 //   ROOT
 #include "TChain.h"
+using namespace NPUNITS;
 ///////////////////////////////////////////////////////////////////////////
 
 ClassImp(TSamuraiFDC2Physics)
@@ -44,6 +46,7 @@ ClassImp(TSamuraiFDC2Physics)
     m_PreTreatedData    = new TSamuraiFDC2Data ;
     m_EventPhysics      = this ;
     //m_Spectra           = NULL;
+    ToTThreshold = 0;
   }
 
 ///////////////////////////////////////////////////////////////////////////
@@ -54,17 +57,35 @@ void TSamuraiFDC2Physics::BuildSimplePhysicalEvent(){
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::BuildPhysicalEvent(){
   PreTreat();
+  RemoveNoise();
+ 
+  // one map per detector
+  map<unsigned int, vector<double> > X ; 
+  map<unsigned int, vector<double> > Z ; 
+  map<unsigned int, vector<double> > R ; 
+
+  unsigned int size = Detector.size();
+  for(unsigned int i = 0 ; i < size ; i++){
+      
+    
+    }
+
   return;
 }
 
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiFDC2Physics::Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& dirX, double& dirZ,double& refX ){
+  
+  }
+
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::PreTreat(){
   ClearPreTreatedData();
   unsigned int size = m_EventData->Mult();
   for(unsigned int i = 0 ; i < size ; i++){
-    // EDGE=0 is the leading edge, IE, the real time.
-    // EDGE=1 is the trailing edge, so it helps build Tot
-    if(m_EventData->GetEdge(i)==0){
+    // EDGE=1 is the leading edge, IE, the real time.
+    // EDGE=0 is the trailing edge, so it helps build Tot
+    if(m_EventData->GetEdge(i)==1){
       int det   = m_EventData->GetDetectorNbr(i); 
       int layer = m_EventData->GetLayerNbr(i); 
       int wire  = m_EventData->GetWireNbr(i); 
@@ -72,7 +93,7 @@ void TSamuraiFDC2Physics::PreTreat(){
       double etime = 0;
       // look for matching trailing edge   
       for(unsigned int j = 0 ; j < size ; j++){
-        if(m_EventData->GetEdge(j)==1){
+        if(m_EventData->GetEdge(j)==0){
           int edet   = m_EventData->GetDetectorNbr(j); 
           int elayer = m_EventData->GetLayerNbr(j); 
           int ewire  = m_EventData->GetWireNbr(j); 
@@ -80,27 +101,66 @@ void TSamuraiFDC2Physics::PreTreat(){
           if(wire==ewire && layer==elayer && det==edet){
             etime = m_EventData->GetTime(j); 
           }    
-          
-        //if(etime<time)
-        //  break;
-        //else
-        //  etime=0;
         }
+        if(etime && etime>time)
+          break;
+        else
+          etime=0;
       }
       // a valid wire must have an edge
-      if(etime && time){
+      if(etime && time && etime-time>ToTThreshold){
         Detector.push_back(det);
         Layer.push_back(layer);       
         Wire.push_back(wire);
         Time.push_back(time);
-        ToT.push_back(time-etime);
-        }
+        ToT.push_back(etime-time);
+        SamuraiDCIndex idx(det,layer,wire);
+      }
     }
 
   }
   return;
 }
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiFDC2Physics::RemoveNoise(){
+  // Remove the noise by looking if a matching wire exist in the adjacent layer
+  // this done by looking at the closest plane with the same orientation
+  unsigned int size = Detector.size(); 
+  for(unsigned int i = 0 ; i < size ; i++){
+    bool match=false;
+    int det = Detector[i];
+    int layer = Layer[i];
+    int wire = Wire[i];
+    // look for matching adjacent wire   
+
+    for(unsigned int j = 0 ; j < size ; j++){
+      int adet = Detector[j];
+      int alayer = Layer[j];
+      int awire = Wire[j];
+      bool blayer = false;
+      if(layer%2==0){
+        if(layer+1==alayer)
+          blayer=true;
+      }
+
+      else{
+        if(layer-1==alayer)
+          blayer=true;
+      }
 
+      if(det==adet && blayer && abs(wire-awire)<=1){
+        match=true;
+        break;
+        }
+      }
+     
+     if(match)
+       Matched.push_back(true);
+     else
+       Matched.push_back(false);
+    }
+  return;
+}
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::Clear(){
   Detector.clear();
@@ -110,6 +170,7 @@ void TSamuraiFDC2Physics::Clear(){
   ToT.clear();
   ParticleDirection.clear();
   MiddlePosition.clear();
+  Matched.clear();
 }
 ///////////////////////////////////////////////////////////////////////////
 
@@ -117,7 +178,55 @@ void TSamuraiFDC2Physics::Clear(){
 
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::ReadConfiguration(NPL::InputParser parser){
+  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIFDC2");
+  if(NPOptionManager::getInstance()->GetVerboseLevel())
+    cout << "//// " << blocks.size() << " detector(s) found " << endl; 
+
+  vector<string> token= {"XML"};
+
+  for(unsigned int i = 0 ; i < blocks.size() ; i++){
+    cout << endl << "////  Samurai FDC2 (" << i+1 << ")" << endl;
+    string xmlpath = blocks[i]->GetString("XML");
+    NPL::XmlParser xml;
+    xml.LoadFile(xmlpath);
+    AddDC("SAMURAIFDC2",xml);
+  }
 }
+
+///////////////////////////////////////////////////////////////////////////
+void TSamuraiFDC2Physics::AddDC(string name, NPL::XmlParser& xml){
+  std::vector<NPL::block*> b = xml.GetAllBlocksWithName(name);  
+  // FDC2 case
+  if(name=="SAMURAIFDC2"){
+    unsigned int det=2;
+    unsigned int size = b.size();
+    for(unsigned int i = 0 ; i < size ; i++){
+      unsigned int layer = b[i]->AsInt("layer"); 
+      unsigned int wire  = b[i]->AsInt("wireid"); 
+      double X = b[i]->AsDouble("wirepos");  
+      double Z = b[i]->AsDouble("wirez");  
+      string sDir = b[i]->AsString("anodedir");
+      double T=0;
+      if(sDir=="X")
+        T=0*deg;
+      else if(sDir=="Y")
+        T=90*deg;
+      else if(sDir=="U")
+        T=30*deg;
+      else if(sDir=="V")
+        T=-30*deg;
+      else{
+        cout << "ERROR: Unknown layer orientation for Samurai FDC2"<< endl;
+        exit(1);
+      }
+     SamuraiDCIndex idx(det,layer,wire);
+     Wire_X[idx]=X;
+     Wire_Z[idx]=Z;
+     Wire_T[idx]=T;
+    }
+  }
+}
+
 ///////////////////////////////////////////////////////////////////////////
 void TSamuraiFDC2Physics::InitSpectra(){  
   //m_Spectra = new TSamuraiFDC2Spectra(m_NumberOfDetector);

NPLib/Detectors/Samurai/TSamuraiFDC2Physics.h

@@ -24,6 +24,7 @@
  *****************************************************************************/
 // STL
 #include <vector>
+#include <map>
 
 // NPL
 #include "TSamuraiFDC2Data.h"
@@ -31,6 +32,7 @@
 #include "NPCalibrationManager.h"
 #include "NPVDetector.h"
 #include "NPInputParser.h"
+#include "NPXmlParser.h"
 // ROOT 
 #include "TVector2.h" 
 #include "TVector3.h" 
@@ -42,6 +44,37 @@
 
 
 using namespace std ;
+// little class to index each of the DC wire
+class SamuraiDCIndex{
+  public:
+   SamuraiDCIndex(){};  
+   ~SamuraiDCIndex(){};  
+   SamuraiDCIndex(unsigned int det, unsigned int layer, unsigned int wire){
+     m_det=det;
+     m_layer=layer;
+     m_wire=wire;
+     m_norme=Norme();
+   };  
+  
+  unsigned int m_det;
+  unsigned int m_layer;
+  unsigned int m_wire;
+  unsigned int m_norme;
+    
+  int Norme() const {return (m_det*1000000000+m_layer*1000000+m_wire);} ;
+  bool operator<(const SamuraiDCIndex i2){
+    return this->Norme()<i2.Norme();
+    }
+  
+  friend bool operator<(const SamuraiDCIndex i1,const SamuraiDCIndex i2){
+    return i1.Norme()<i2.Norme();
+    }
+  
+  friend bool operator==(const SamuraiDCIndex i1,const SamuraiDCIndex i2){
+   return i1.Norme()==i2.Norme();
+   }
+  };
+
 
 class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
   public:
@@ -59,7 +92,7 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
     vector<int> Wire;
     vector<double> Time;
     vector<double> ToT;
-
+    vector<bool>   Matched;
     // Computed variable
     vector<TVector3> ParticleDirection;
     vector<TVector3> MiddlePosition;
@@ -68,6 +101,19 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
     // Projected position at given Z plan
     TVector3 ProjectedPosition(double Z);
 
+  private: // Charateristic of the DC 
+    void AddDC(string name, NPL::XmlParser&);//! take the XML file and fill in Wire_X and Layer_Angle
+    map<SamuraiDCIndex,double> Wire_X;//! X position of the wires
+    map<SamuraiDCIndex,double> Wire_Z;//! Z position of the wires
+    map<SamuraiDCIndex,double> Wire_T;//! Wire Angle (0 for X, 90 for Y, U and V are typically at +/-30)
+  
+  private: // Analysis
+    double ToTThreshold;//! a ToT threshold to remove noise
+    void RemoveNoise();
+    // Construct the 2D track and ref position at Z=0 based on X,Z and Radius provided
+    void Track2D(const vector<double>& X,const vector<double>& Z,const vector<double>& R,double& dirX, double& dirZ,double& refX )
+
+
   public: //   Innherited from VDetector Class
 
     // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token

Projects/e748/BeamID.cxx

+void BeamID(){
+  // Load the Main Tree
+  TFile* file = new TFile("../../Outputs/Analysis/e748_Physics_12Be.root");
+  TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+  TFile* fileR = new TFile("../../Outputs/Analysis/e748_12Be.root");
+  TTree* treeR = (TTree*) fileR->FindObjectAny("ResultTree");
+  tree->AddFriend(treeR);
+
+  // Load the IC chain
+  TChain* IC = new TChain("Numexo2");
+  IC->Add("/data/Transfert/e748/Merged/run_0315.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0316.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0317.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0318.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0320.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0321.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0323.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0325.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0326.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0327.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0328.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0329.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0330.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0331.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0339.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0341.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0342.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0346.root");
+  IC->Add("/data/Transfert/e748/Merged/run_0347.root");