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

NPLib/Detectors/BigRIPS/CMakeLists.txt

@@ -10,8 +10,12 @@ add_custom_command(OUTPUT TBigRIPSICDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scr
 
 add_custom_command(OUTPUT TBigRIPSICPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSICPhysics.h TBigRIPSICPhysicsDict.cxx TBigRIPSICPhysics.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSICPhysics.h)
 
-add_library(NPBigRIPS SHARED TBigRIPSPPACData.cxx TBigRIPSPPACDataDict.cxx TBigRIPSPPACPhysics.cxx TBigRIPSPPACPhysicsDict.cxx TBigRIPSPlasticData.cxx TBigRIPSPlasticDataDict.cxx TBigRIPSPlasticPhysics.cxx TBigRIPSPlasticPhysicsDict.cxx TBigRIPSICData.cxx TBigRIPSICDataDict.cxx TBigRIPSICPhysics.cxx TBigRIPSICPhysicsDict.cxx)
+add_custom_command(OUTPUT TBigRIPSFocalDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSFocal.h TBigRIPSFocalDict.cxx TBigRIPSFocal.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSFocal.h)
+
+add_custom_command(OUTPUT TBigRIPSRecoDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSReco.h TBigRIPSRecoDict.cxx TBigRIPSReco.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSReco.h)
+
+add_library(NPBigRIPS SHARED TBigRIPSPPACData.cxx TBigRIPSPPACDataDict.cxx TBigRIPSPPACPhysics.cxx TBigRIPSPPACPhysicsDict.cxx TBigRIPSPlasticData.cxx TBigRIPSPlasticDataDict.cxx TBigRIPSPlasticPhysics.cxx TBigRIPSPlasticPhysicsDict.cxx TBigRIPSICData.cxx TBigRIPSICDataDict.cxx TBigRIPSICPhysics.cxx TBigRIPSICPhysicsDict.cxx TBigRIPSFocal.cxx TBigRIPSFocalDict.cxx TBigRIPSReco.cxx TBigRIPSRecoDict.cxx)
 
 target_link_libraries(NPBigRIPS ${ROOT_LIBRARIES} NPCore NPTrackReconstruction) 
-install(FILES TBigRIPSPPACData.h TBigRIPSPPACPhysics.h TBigRIPSPlasticData.h TBigRIPSPlasticPhysics.h TBigRIPSICData.h TBigRIPSICPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TBigRIPSPPACData.h TBigRIPSPPACPhysics.h TBigRIPSPlasticData.h TBigRIPSPlasticPhysics.h TBigRIPSICData.h TBigRIPSICPhysics.h TBigRIPSFocal.h TBigRIPSReco.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
 

NPLib/Detectors/BigRIPS/TBigRIPSFocal.cxx

+#include "TBigRIPSFocal.h"
+
+TBigRIPSFocal::TBigRIPSFocal(){};
+TBigRIPSFocal::~TBigRIPSFocal(){Clear();};
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSFocal::Clear(){
+    FPTrack.clear();
+}
+////////////////////////////////////////////////////////////////////////////////
+void  TBigRIPSFocal::Init(int N_FP){
+    std::vector<double> tmp;
+    for(int j=0;j<4;j++) tmp.push_back(-9999);
+    for(int i=0;i<N_FP;i++) FPTrack.push_back(tmp);
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSFocal::Print(int FPNbr){
+    std::cout << "----FP----:" << FPNbr <<std::endl;
+    std::cout << "X: " << FPTrack[FPNbr][0] << std::endl;
+    std::cout << "A: " << FPTrack[FPNbr][1] << std::endl;
+    std::cout << "Y: " << FPTrack[FPNbr][2] << std::endl;
+    std::cout << "B: " << FPTrack[FPNbr][3] << std::endl;
+
+}
+ClassImp(TBigRIPSFocal); 

NPLib/Detectors/BigRIPS/TBigRIPSFocal.h

+#ifndef TBigRIPSFocal_H
+#define TBigRIPSFocal_H
+#include "TObject.h"
+#include <vector>
+#include <iostream>
+class TBigRIPSFocal: public TObject{
+
+  public:
+    TBigRIPSFocal();
+    ~TBigRIPSFocal();
+
+  private:
+    std::vector<std::vector<double>> FPTrack  ;
+
+  public:
+    void Clear();
+    void Init(int);
+    void Print(int);
+  
+  public:
+    void SetFPTrack(int FPNbr, std::vector<double> track){
+            FPTrack[FPNbr]=track;
+    };
+    void SetFPTrack(int FPNbr, double x, double a, double y, double b){
+            FPTrack[FPNbr][0]=x;
+            FPTrack[FPNbr][1]=a;
+            FPTrack[FPNbr][2]=y;
+            FPTrack[FPNbr][3]=b;
+    };
+    std::vector<double> GetFPTrack(int FPNbr) {return FPTrack[FPNbr];}
+
+    ClassDef(TBigRIPSFocal,1); 
+};
+
+#endif

NPLib/Detectors/BigRIPS/TBigRIPSPPACData.cxx

@@ -4,30 +4,8 @@
 TBigRIPSPPACData::TBigRIPSPPACData(){};
 TBigRIPSPPACData::~TBigRIPSPPACData(){};
 
-/*
-void TBigRIPSPPACData::SetData(const int& FP, const int& ID, const int& Edge, const int& value,const int& VariableType){
-  fPPAC_FP.push_back(FP);
-  fPPAC_ID.push_back(ID);
-  fPPAC_Edge.push_back(Edge);
-  switch(VariableType){
-      case 0: fPPAC_TX1.push_back(value); break; 
-      case 1: fPPAC_TX2.push_back(value); break; 
-      case 2: fPPAC_TY1.push_back(value); break; 
-      case 3: fPPAC_TY2.push_back(value); break; 
-      case 4: fPPAC_TA.push_back(value); break; 
-      case 5: fPPAC_QX1.push_back(value); break; 
-      case 6: fPPAC_QX2.push_back(value); break; 
-      case 7: fPPAC_QY1.push_back(value); break; 
-      case 8: fPPAC_QY2.push_back(value); break; 
-      case 9: fPPAC_QA.push_back(value); break; 
-  }
-}
-*/
 ////////////////////////////////////////////////////////////////////////////////
 void TBigRIPSPPACData::Clear(){
- // fPPAC_FP.clear();
- // fPPAC_ID.clear();
- // fPPAC_Edge.clear();
   fPPAC_TX1.clear(); 
   fPPAC_TX2.clear(); 
   fPPAC_TY1.clear(); 
@@ -59,32 +37,4 @@ void TBigRIPSPPACData::Print(){
   //cout << "   - Multiplicity: " << Mult() << endl;
 
 }
-////////////////////////////////////////////////////////////////////////////////
-/*
-unsigned int TBigRIPSPPACData::MultLayer(unsigned int det , unsigned int layer, int edge){
-  unsigned int mult=0;
-  unsigned int size = fDC_DetectorNbr.size();
-  for(unsigned int i = 0 ; i< size ; i++ ){
-    if(fDC_DetectorNbr[i]==det)
-      if(fDC_LayerNbr[i]==layer)
-        if(fDC_Edge[i]==edge && edge!=-1) // edge type is specified (0 or 1)
-          mult++;
-        else if(edge==-1)// edge type is not specified
-          mult++;
-  }
-  return mult;
-
-}
-////////////////////////////////////////////////////////////////////////////////
-std::vector<int> TBigRIPSPPACData::GetWire(unsigned int det , unsigned int layer){
-  std::vector<int> wires;
-  unsigned int size = fDC_DetectorNbr.size();
-  for(unsigned int i = 0 ; i< size ; i++ ){
-    if(fDC_DetectorNbr[i]==det)
-      if(fDC_LayerNbr[i]==layer)
-        wires.push_back(fDC_WireNbr[i]);
-  }
-  return wires;
-}
-*/
 ClassImp(TBigRIPSPPACData); 

NPLib/Detectors/BigRIPS/TBigRIPSReco.cxx

+#include "TBigRIPSReco.h"
+
+TBigRIPSReco::TBigRIPSReco(){Init();}
+TBigRIPSReco::~TBigRIPSReco(){Clear();}
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Clear(){
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Init(){
+    aoq=-9999;
+    beta=-9999;
+    delta=-9999;
+    brho=-9999;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void TBigRIPSReco::RecBrho(std::vector<double> RecFPUpstream,std::vector<double> RecFPDownstream, std::vector<std::vector<double>> matrix,double BrhoCentral){
+
+    if(matrix.size()!=6) {
+        std::cout << "MATRIX used for Brho Rec HAS NOT 6 rows" << std::endl;
+        return;
+    }else{
+        for(int i=0;i<6;i++){
+            if(matrix[i].size()!=5){
+                 std::cout << "MATRIX used for Brho Rec. HAS NOT 5 columns" << std::endl;
+                 std::cout << "but" << matrix[i].size()<<std::endl;
+                 return;
+            }
+        }
+    }
+
+    if(RecFPUpstream[0]==-9999 || 
+       RecFPUpstream[1]==-9999 || 
+       RecFPDownstream[0]==-9999)
+        return;
+
+    double x_x = matrix[0][0];
+    double a_x = matrix[0][1];
+    double x_a = matrix[1][0];
+    double x_d = matrix[5][0];
+    double a_a = matrix[1][1];
+    double a_d = matrix[5][1];
+
+    double x_up = RecFPUpstream[0]; // downst. X
+    double a_up = RecFPUpstream[1]; // downst. X
+    double x_down = RecFPDownstream[0]; // downst. X
+    double a_down = RecFPDownstream[1]; // downst. A
+
+    delta = (a_a * x_down - x_a * a_down - x_up) / (x_d * a_a - x_a * a_d);
+    angle = (a_a * x_down - x_a * a_down - x_up) / (x_d * a_a - x_a * a_d);
+
+    brho = BrhoCentral*(1.0+delta*0.01);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::RecAoqOne(double tof, double length){
+    beta = length /(tof * c_mm_ns);
+    double gamma = 1./sqrt(1 - beta*beta);
+    aoq = (brho * c_mm_ns) / (mnucleon * beta * gamma);
+//std::cout << aoq << std::endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::RecZet(double dE,double ionpair, std::vector<double> coeff){
+   if(coeff.size()<2){
+       std::cout << "Zcoeff for Z reco not defined" << std::endl;
+       return;
+   }
+   double de_v = log(ionpair*beta*beta) - log((1-beta*beta)) - beta*beta;
+   z = coeff[0] * sqrt(dE/de_v)*beta + coeff[1];   
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Print(){
+    std::cout << "aoq: "  << aoq << std::endl;
+    std::cout << "z: "  << z << std::endl;
+    std::cout << "beta: " << beta << std::endl;
+    std::cout << "delta: "<< delta << std::endl;
+    std::cout << "brho: " << brho << std::endl;
+}
+
+ClassImp(TBigRIPSReco); 

NPLib/Detectors/BigRIPS/TBigRIPSReco.h

+#ifndef TBigRIPSReco_H
+#define TBigRIPSReco_H
+#include "TObject.h"
+#include <vector>
+#include <iostream>
+class TBigRIPSReco: public TObject{
+
+  public:
+    TBigRIPSReco();
+    ~TBigRIPSReco();
+
+  private:
+    double aoq;
+    double beta;
+    double delta;
+    double brho;
+    double angle;
+    double z;
+  public:
+    void Clear();
+    void Init();
+    void Print();
+    void RecBrho(std::vector<double>,std::vector<double>,std::vector<std::vector<double>>,double);
+    void RecAoqOne(double, double);
+    void RecZet(double,double, std::vector<double>);
+    const double c_mm_ns = 299.7792458; // in mm/ns
+    const double mnucleon = 931.49432;  // MeV
+  
+  public:
+    void SetAoq(double value){aoq=value;}
+    void SetBeta(double value){beta=value;}
+    void SetDelta(double value){delta=value;}
+    void SetBrho(double value){brho=value;}
+    void SetZ(double value){z=value;}
+
+    ClassDef(TBigRIPSReco,1); 
+};
+
+#endif

Projects/RIBF196/Analysis.cxx

@@ -40,6 +40,14 @@ Analysis::~Analysis(){
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::Init() {
 
+  NmaxFP=11;
+  //FP = new TBigRIPSFocal();
+  //FP->Init(NmaxFP+1);
+  Rec35 = new TBigRIPSReco();
+  Rec57 = new TBigRIPSReco();
+  Rec89 = new TBigRIPSReco();
+  Rec911 = new TBigRIPSReco();
+
   InitOutputBranch();
   InitInputBranch();
   // get PPAC, PL, and IC  objects
@@ -48,9 +56,25 @@ void Analysis::Init() {
   IC = (TBigRIPSICPhysics*)  m_DetectorManager -> GetDetector("BigRIPSIC");
 
   ReadXmls();
-
+/*
   matF35.ResizeTo(6,5); matF35 = RecReadTransferMatrix("mat1.mat");
   matF57.ResizeTo(6,5); matF57 = RecReadTransferMatrix("mat2.mat");
+  matF89.ResizeTo(6,5); matF89 = RecReadTransferMatrix("F8F9_LargeAccAchr.mat");
+  matF911.ResizeTo(6,5); matF911 = RecReadTransferMatrix("F9F11_LargeAccAchr.mat");
+*/
+  matrixF35  = RecReadTransferMatrix2("mat1.mat");
+  matrixF57  = RecReadTransferMatrix2("mat2.mat");
+  matrixF89  = RecReadTransferMatrix2("F8F9_LargeAccAchr.mat");
+  matrixF911 = RecReadTransferMatrix2("F9F11_LargeAccAchr.mat");
+
+  length37  = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) 
+             -(FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]); 
+  length35  =  FP_Z[FPtoID[5]] - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
+  length57  = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - FP_Z[FPtoID[5]];
+  length811 = (FP_Z[FPtoID[11]]+PL_NametoZ["F11pl-1"]) 
+             -(FP_Z[FPtoID[8]]+PL_NametoZ["F8pl"]); 
+  length89  =  FP_Z[FPtoID[9]] - (FP_Z[FPtoID[8]]+PL_NametoZ["F8pl"]);
+  length911 = (FP_Z[FPtoID[11]]+PL_NametoZ["F11pl-1"]) - FP_Z[FPtoID[9]];
 
   // initialize various parameters
   //Rand = TRandom3();
@@ -111,11 +135,15 @@ void Analysis::ReadXmls() {
       IC_IDtoName[ID] = bic[i]->AsString("NAME"); 
       string name = bic[i]->AsString("NAME");
       IC_NametoID[name] = ID; 
+      IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_0")); 
+      IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_1")); 
+      IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_2")); 
+      IC_Ionpair[name] = bic[i]->AsDouble("ionpair"); 
   }
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::TreatEvent() {
-    // Reinitiate calculated variable
+    // Reinitiate and clear variables
     ReInitValue();
     FP_PPACHitList.clear();
 
@@ -124,162 +152,106 @@ void Analysis::TreatEvent() {
     /////////////////////////////////////////////
 
     // Get all PPAC X,Y hits and store them in a HitList indexed by Focalplane
+    // For a given Focal plane, the HitList size depends on how many PPAC 
+    // are installed at this FP
     for(int i=0;i<PPAC->ID.size();i++){ 
         FP_PPACHitList[PPAC->FP[i]].AddPPACHit(PPAC->X[i],PPAC->Y[i],PPAC->ID[i]);
     }
 
-    // Perform FP reco. and store in map: index FP, content vector(x,a,y,b)
-    std::map<int,vector<double>> RecFP;
-    double NmaxFP=12;
-    for(int i=0;i<NmaxFP;i++){
-        for(int j=0;j<4;j++){
-            RecFP[i].push_back(-9999);}}
-
+    // For each FP, perform reco. and store in vector<vector>
+    // RecFP[FPNbr][0,1,2,3 for x,a,y,b]
+    for(auto it = FP_PPACHitList.begin();it!=FP_PPACHitList.end();++it){
+        if(it->first < NmaxFP){
+           RecFP[it->first] = RecFPposition(it->second.PPACHit,it->second.PPACID);
+        }
+    }
+    // Same as above but using an object from class TBigRIPSFocal
+    /*
     for(auto it = FP_PPACHitList.begin();it!=FP_PPACHitList.end();++it){
-        RecFP[it->first] = RecFPposition(it->second.PPACHit,it->second.PPACID);
+        if(it->first < NmaxFP){
+           FP->SetFPTrack(it->first, RecFPposition(it->second.PPACHit,it->second.PPACID) );
+           //FP->Print(it->first);
+        }
     }
+*/
 
-    fX=RecFP[8][0];
-    fA=RecFP[8][1];
-    fY=RecFP[8][2];
-    fB=RecFP[8][3];
-
-    delta35 = RecDeltaBrho(RecFP[3],RecFP[5],matF35);
-    Brho35 = BrhoD3 * (1.0 + delta35*0.01);
+    // Brho/delta reconstruction
+    Rec35->RecBrho(RecFP[3],RecFP[5],matrixF35,BrhoD3);
+    Rec57->RecBrho(RecFP[5],RecFP[7],matrixF57,BrhoD5);
+    Rec89->RecBrho(RecFP[8],RecFP[9],matrixF89,BrhoD7);
+    Rec911->RecBrho(RecFP[9],RecFP[11],matrixF911,BrhoD8);
 
-    delta57 = RecDeltaBrho(RecFP[5],RecFP[7],matF57);
-    Brho57 = BrhoD5 * (1.0 + delta57*0.01);
 
     /////////////////////////////////
     /// STEP2: TOF, flight length ///
     /////////////////////////////////
     
     std::map<string,double> PL_T; //index is plastic name in xml 
-    for(int i=0;i<PL->ID.size();i++) {PL_T[PL_IDtoName[PL->ID[i]]] = PL->T[i];
-        //if(PL->ID[i]==2) std::cout << "PL ID2 T::" << PL->T[i]  << std::endl;
-    }
+    for(int i=0;i<PL->ID.size();i++) {PL_T[PL_IDtoName[PL->ID[i]]] = PL->T[i];}
     tf3 = PL_T["F3pl"];
     tf7 = PL_T["F7pl"];
-    if(PL_T["F7pl"]!=-9999 && PL_T["F3pl"]!=-9999)  tof37 = PL_T["F7pl"] - PL_T["F3pl"];
+    if(PL_T["F7pl"]!=-9999 && PL_T["F3pl"]!=-9999) tof37 = PL_T["F7pl"] - PL_T["F3pl"];
     tof37 += tof_offset37;
-    
-    //double length35 = (FP_Z[5]+PL_NametoZ["F5pl"]) - (FP_Z[3]+PL_NametoZ["F3pl"]) ; 
-    //double length57 = (FP_Z[7]+PL_NametoZ["F7pl"]) - (FP_Z[5]+PL_NametoZ["F5pl"]) ; 
-    //double length37 = length35 + length57 ; 
-
-    double length37 =  (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]); 
-    double length35 =   FP_Z[FPtoID[5]] - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
-    double length57 =  (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - FP_Z[FPtoID[5]];
-
-    //std::cout << "length37:" << length37  << std::endl;
-    //std::cout << "Z FP7:" << FP_Z[FPtoID[7]] << std::endl;
-    //std::cout << "Zoff PlasticF7:" << PL_NametoZ["F7pl"] << std::endl;
-    //std::cout << "Z FP3:" << FP_Z[FPtoID[3]] << std::endl;
-    //std::cout << "Zoff PlasticF3:" << PL_NametoZ["F3pl"] << std::endl;
 
-    beta35 = length35 /(tof37/2. * clight);
-    beta57 = length57 /(tof37/2. * clight);
-    beta37 = length37 /(tof37 * clight);
+    tf8 = PL_T["F8pl"];
+    tf11 = PL_T["F11pl-1"];
+    if(PL_T["F11pl-1"]!=-9999 && PL_T["F8pl"]!=-9999) tof811 = PL_T["F11pl-1"] - PL_T["F8pl"];
+    tof811 += tof_offset811;
 
+    
     /////////////////////////////////
     /// STEP3: Beta, Gamma, AoQ   ///
     /////////////////////////////////
 
-    //   std::cout << "______________________" << std::endl;
-    //   std::cout << "Brho35:" << Brho35  << std::endl;
-    //   std::cout << "Tof37:" << tof37  << std::endl;
-    //   std::cout << "length35:" << length35  << std::endl;
-    aoq35 = RecAoqOneFold(Brho35, tof37, length37);
-    //   std::cout << "aoq35:" << aoq35  << std::endl;
-
-    //   std::cout << "______________________" << std::endl;
-    //   std::cout << "Brho57:" << Brho57  << std::endl;
-    //   std::cout << "Tof37:" << tof37  << std::endl;
-    //   std::cout << "length57:" << length57  << std::endl;
-    aoq57 = RecAoqOneFold(Brho57, tof37, length37);
-    //   std::cout << "aoq57:" << aoq57  << std::endl;
+    // Calculate Beta from length/tof and combine it with Brho to get AoQ
+    Rec35->RecAoqOne(tof37,length37);
+    Rec57->RecAoqOne(tof37,length37);
+    Rec89->RecAoqOne(tof811,length811);
+    Rec911->RecAoqOne(tof811,length811);
     
-    //////////////////
-    /// STEP4: Z   ///
-    //////////////////
+    ///////////////////////////////////////////////
+    /// STEP4: Z reco. from dE in IC and Beta   ///
+    ///////////////////////////////////////////////
     
     std::map<string,double> IC_dE; //index is IC name in xml 
-    for(int i=0;i<IC->ID.size();i++) {IC_dE[IC_IDtoName[IC->ID[i]]] = IC->CalSqSum[i];
+    for(int i=0;i<IC->ID.size();i++) {
+        IC_dE[IC_IDtoName[IC->ID[i]]] = IC->CalSqSum[i];
         //if(IC->ID[i]==1) std::cout << "IC ID1 CalSq:" << IC->CalSqSum[i]  << std::endl;
     }
     dE_ICF7 = IC_dE["F7IC"];
-    z_BR = RecZ(dE_ICF7, tof37, length37);
+    Rec35->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
+    Rec57->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
+
+    dE_ICF11 = IC_dE["F11IC"];
+    Rec89->RecZet(dE_ICF11,IC_Ionpair["F11IC"],IC_Zcoef["F11IC"]);
+    Rec911->RecZet(dE_ICF11,IC_Ionpair["F11IC"],IC_Zcoef["F11IC"]);
 
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-TMatrixD Analysis::RecReadTransferMatrix(string matfile){
+std::vector<std::vector<double>> Analysis::RecReadTransferMatrix2(string matfile){
   char buffer[256];
   std::ifstream matfin(matfile); 
   if(matfin.fail()) cout<<"Failed to open optical matrix file: "<<matfile<<endl;
   matfin.getline(buffer,256);
 
-  TMatrixD matrix(6,5);
+  std::vector<std::vector<double>> matrix;
+  //matrix.resize(5, vector<double>(6));
+
   double val[6];  
+  std::vector<double> row;
 
   for(Int_t i=0;i<6;i++){
+    row.clear();
     matfin>>val[0]>>val[1]>>val[2]>>val[3]>>val[4]>>val[5];
     //cout<<val[0]<<val[1]<<val[2]<<val[3]<<val[4]<<val[5]<<endl;
-    for(Int_t j=0;j<5;j++)
-      matrix(i,j) = val[j];
+    for(Int_t j=0;j<5;j++) row.push_back(val[j]);
+    matrix.push_back(row);
   }
   matfin.close();
-
-  matrix.Print();
   return matrix;
 }
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecDeltaBrho(std::vector<double> RecFPUpstream,std::vector<double> RecFPDownstream, TMatrixD matrix){
-
-    if(matrix.GetNrows()!=6 && matrix.GetNcols()!=5){
-        std::cout << "MATRIX used for Brho Rec. is NULL" << std::endl;
-        return -9999;
-    }
-
-    double x_a = matrix(0,0);
-
-    TMatrixD mat1(2,1);
-    mat1(0,0) = matrix(0,0); mat1(1,0) = matrix(0,1);
-    //mat1.Print();
 
-    TMatrixD mat2(2,2);
-    mat2(0,0) = matrix(1,0); mat2(0,1) = matrix(5,0);
-    mat2(1,0) = matrix(1,1); mat2(1,1) = matrix(5,1);
-    //mat2.Print();
-    mat2.Invert();
-
-    //TVectorD * fvec = (TVectorD *)fUpstreamFplArrayBuffer[i]->GetOptVector();
-    //TVectorD * bvec = (TVectorD *)fDownstreamFplArrayBuffer[i]->GetOptVector();
-
-    if(RecFPUpstream[0]==-9999 || 
-       RecFPUpstream[1]==-9999 || 
-       RecFPDownstream[0]==-9999)
-        return -9999;
-
-    TMatrixD xvec(2,1);
-    //xvec(0,0) = (*bvec)(0);
-    //xvec(1,0) = (*bvec)(1);
-    xvec(0,0) = RecFPDownstream[0]; // downst. X
-    xvec(1,0) = RecFPDownstream[1]; // downst. A
-
-    TMatrixD rvec = mat2 * (xvec - RecFPUpstream[0] * mat1);
-    Double_t angle = rvec(0,0); // change to mrad
-    Double_t delta = rvec(1,0);
-    // simple_delta = (*bvec)(0) / 33.; // only for fpl-5
-
-    //double BrhoRec = BrhoCentral*(1.0+delta*0.01);
-    //rips->SetDelta(delta);
-    //rips->SetAngle(angle);
-
-
-    //fReconstructed = true;
-    return delta;
-}
 ////////////////////////////////////////////////////////////////////////////////
 // Function transforming a list of PPAC Hit positions and IDs for a given FP
 // into  positions and angles at the FP position (FX,FY,FA,FB)
@@ -376,47 +348,29 @@ std::vector<double> Analysis::RecFPposition(std::vector<TVector2> HitList,std::v
 
 }
 
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecAoqOneFold(double Brho, double tof, double length){
-    double beta = length /(tof * clight);
-    double gamma = 1./sqrt(1 - beta*beta);
-    double aoq = (Brho * clight) / (mnucleon * beta * gamma);
-    return aoq;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecZ(double dE, double tof, double length){
- double ionpair = 4866.;
- double beta = length /(tof * clight);
- Double_t de_v = log(ionpair*beta*beta) - log((1-beta*beta)) - beta*beta;
- double z = 17.9143234533*sqrt(dE/de_v)*beta -13.0990490522;   
- return z;
-}
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::End(){
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch() {
 
-  RootOutput::getInstance()->GetTree()->Branch("EventNumber",&EventNumber,"EventNumber/I");
   RootOutput::getInstance()->GetTree()->Branch("RunNumber",&RunNumber,"RunNumber/I");
-  RootOutput::getInstance()->GetTree()->Branch("fX",&fX,"fX/D");
-  RootOutput::getInstance()->GetTree()->Branch("fY",&fY,"fY/D");
-  RootOutput::getInstance()->GetTree()->Branch("delta35",&delta35,"delta35/D");
-  RootOutput::getInstance()->GetTree()->Branch("delta57",&delta57,"delta57/D");
-  RootOutput::getInstance()->GetTree()->Branch("Brho35",&Brho35,"Brho35/D");
-  RootOutput::getInstance()->GetTree()->Branch("Brho57",&Brho57,"Brho57/D");
+  RootOutput::getInstance()->GetTree()->Branch("EventNumber",&EventNumber,"EventNumber/I");
+  RootOutput::getInstance()->GetTree()->Branch("Trigger",&Trigger,"Trigger/I");
+  RootOutput::getInstance()->GetTree()->Branch("TimeStamp",&TimeStamp,"TimeStamp/I");
+  RootOutput::getInstance()->GetTree()->Branch("RecFP",&RecFP);
+  //RootOutput::getInstance()->GetTree()->Branch("FP","TBigRIPSFocal",&FP);
+  RootOutput::getInstance()->GetTree()->Branch("Rec35","TBigRIPSReco",&Rec35);
+  RootOutput::getInstance()->GetTree()->Branch("Rec57","TBigRIPSReco",&Rec57);
+  RootOutput::getInstance()->GetTree()->Branch("Rec89","TBigRIPSReco",&Rec89);
+  RootOutput::getInstance()->GetTree()->Branch("Rec911","TBigRIPSReco",&Rec911);
   RootOutput::getInstance()->GetTree()->Branch("tof37",&tof37,"tof37/D");
-  RootOutput::getInstance()->GetTree()->Branch("aoq35",&aoq35,"aoq35/D");
-  RootOutput::getInstance()->GetTree()->Branch("aoq57",&aoq57,"aoq57/D");
-  RootOutput::getInstance()->GetTree()->Branch("aoq37",&aoq57,"aoq37/D");
-  RootOutput::getInstance()->GetTree()->Branch("tf3",&tf3,"tf3/D");
+  RootOutput::getInstance()->GetTree()->Branch("tof811",&tof811,"tof811/D");
   RootOutput::getInstance()->GetTree()->Branch("tf7",&tf7,"tf7/D");
-  RootOutput::getInstance()->GetTree()->Branch("beta35",&beta35,"beta35/D");
-  RootOutput::getInstance()->GetTree()->Branch("beta57",&beta57,"beta57/D");
-  RootOutput::getInstance()->GetTree()->Branch("beta37",&beta57,"beta37/D");
-  RootOutput::getInstance()->GetTree()->Branch("z_BR",&z_BR,"z_BR/D");
+  RootOutput::getInstance()->GetTree()->Branch("tf8",&tf8,"tf8/D");
+  RootOutput::getInstance()->GetTree()->Branch("tf11",&tf11,"tf11/D");
   RootOutput::getInstance()->GetTree()->Branch("dE_ICF7",&dE_ICF7,"dE_ICF7/D");
+  RootOutput::getInstance()->GetTree()->Branch("dE_ICF11",&dE_ICF11,"dE_ICF11/D");
 }
 
 
@@ -424,60 +378,29 @@ void Analysis::InitOutputBranch() {
 void Analysis::InitInputBranch(){
   RootInput:: getInstance()->GetChain()->SetBranchAddress("RunNumber",&RunNumber);
   RootInput:: getInstance()->GetChain()->SetBranchAddress("EventNumber",&EventNumber);
-/*
- *
-  if(!simulation){
-    // Vamos
-    RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",&LTS);
-    // Agata
-    RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("nbCores",&nbCores);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
-
-  }
-  else{
-
-    RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
-    RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
-    RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
-    RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
-    RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
-    RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
-
-  }
-*/
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::ReInitValue(){
-  fX = -9999 ;
-  fY = -9999 ;
-  fA = -9999 ;
-  fB = -9999 ;
-  delta35 = -9999 ;
-  delta57 = -9999 ;
-  Brho35 = -9999 ;
-  Brho57 = -9999 ;
+
+  RecFP.clear();
+  std::vector<double> reset;
+  for(int j=0;j<4;j++) reset.push_back(-9999);
+  for(int i=0;i<NmaxFP+1;i++) RecFP.push_back(reset);
+
+  //FP->Clear();
+  //FP->Init(NmaxFP+1);
+  Rec35->Init();
+  Rec57->Init();
+  Rec89->Init();
+  Rec911->Init();
   tf3 =-9999;
   tf7 =-9999;
   tof37 = -9999 ;
-  aoq35 = -9999 ;
-  aoq57 = -9999 ;
-  aoq37 = -9999 ;
-
-  beta35 = -9999 ;
-  beta57 = -9999 ;
-  beta37 = -9999 ;
   dE_ICF7 = -9999;
-  z_BR = -9999 ;
+  tf8 =-9999;
+  tf11 =-9999;
+  tof811 = -9999 ;
+  dE_ICF11 = -9999;
 
 }
 

Projects/RIBF196/Analysis.h

@@ -31,6 +31,8 @@
 #include "TBigRIPSPPACPhysics.h"
 #include "TBigRIPSPlasticPhysics.h"
 #include "TBigRIPSICPhysics.h"
+#include "TBigRIPSFocal.h"
+#include "TBigRIPSReco.h"
 #include <TMatrixD.h>
 #include <TRandom3.h>
 #include <TVector3.h>
@@ -77,11 +79,12 @@ class Analysis: public NPL::VAnalysis{
     void End();
     void ReadXmls();
     //TMatrixD* RecReadTransferMatrix(char*);
-    TMatrixD RecReadTransferMatrix(string);
+    //TMatrixD RecReadTransferMatrix(string);
+    std::vector<std::vector<double>> RecReadTransferMatrix2(string);
     std::vector<double> RecFPposition(std::vector<TVector2>,std::vector<int>);
-    double RecDeltaBrho(std::vector<double>,std::vector<double>, TMatrixD matrix);
-    double RecAoqOneFold(double, double, double);
-    double RecZ(double,double,double);
+    //double RecDeltaBrho(std::vector<double>,std::vector<double>, TMatrixD matrix);
+    //double RecAoqOneFold(double, double, double);
+    double RecZ(double,double,double, bool);
 
     void InitOutputBranch();
     void InitInputBranch();
@@ -90,28 +93,16 @@ class Analysis: public NPL::VAnalysis{
 
   private:
     //calculated variables
-    double aoq;
-    double zet;
-    double beta;
-    double delta35;
-    double delta57;
-    double Brho35;
-    double Brho57;
     double tof37;
-    double aoq35;
-    double aoq57;
-    double aoq37;
-    double beta35;
-    double beta57;
-    double beta37;
+    double tof811;
     double dE_ICF7;
+    double dE_ICF11;
     double z_BR;
+    double z_ZD;
     double tf3;
     double tf7;
-    double fX;
-    double fY;
-    double fA;
-    double fB;
+    double tf8;
+    double tf11;
 
     // List of Dipole Brho, manually set for the moment 
     // But method to read from root file should be implemented
@@ -124,17 +115,24 @@ class Analysis: public NPL::VAnalysis{
     double BrhoD7 = 5.9429;
     double BrhoD8 = 5.9381;
 
-    //double tof_offset37 = 0.;
-    //double tof_offset37 = 300.06;
-    double tof_offset37 = 299.968;
+    double length35;
+    double length57;
+    double length37;
+    double length89;
+    double length911;
+    double length811;
+
+    double tof_offset37 = 300.06;
+    double tof_offset811 = -134.503;
 
     const double clight = 299.7792458; // in mm/ns
     const double mnucleon = 931.49432;  // MeV
 
-
-    // intermediate variables
-    TMatrixD matF35;
-    TMatrixD matF57;
+    //5 columns/ 6 rows transfer matrices
+    vector<vector<double>> matrixF35;
+    vector<vector<double>> matrixF57;
+    vector<vector<double>> matrixF89;
+    vector<vector<double>> matrixF911;
 
 
     // Branches and detectors
@@ -144,6 +142,15 @@ class Analysis: public NPL::VAnalysis{
     map<int,BigRIPSPPACHitList> FP_PPACHitList  ;
     int EventNumber;
     int RunNumber;
+    int Trigger;
+    unsigned long long int TimeStamp;
+    //Focal* FP;
+    //TBigRIPSFocal* FP;
+    TBigRIPSReco* Rec35;
+    TBigRIPSReco* Rec57;
+    TBigRIPSReco* Rec89;
+    TBigRIPSReco* Rec911;
+    std::vector<std::vector<double>> RecFP;
 
     //maps containings infos from input xml config files
     // for PPACs
@@ -159,11 +166,14 @@ class Analysis: public NPL::VAnalysis{
     map<int,int>  IC_IDtoFP;
     map<int, string>  IC_IDtoName; 
     map<string,int>  IC_NametoID; 
+    map<string,vector<double>>  IC_Zcoef;
+    map<string,double>  IC_Ionpair;
     // for FocalPlanes
     map<int,int>  FP_IDtoFP;
     map<int,int>  FPtoID;
     map<int,double> FP_Z;
     map<int,double> FP_Zoffset;
+    int NmaxFP;
 
 };