Update analysis in Vendeta project

Projects/Vendeta/Analysis.cxx

@@ -46,63 +46,170 @@ void Analysis::Init(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::TreatEvent(){
+  
   ReInitValue();
 
+  
+  const double GPeak_Vendlg[4] = {45.68, 45.42, 36.54, 41.6};
+  const double GPeak_Vendhg[4] = {40.26, 43.84, 42.91, 44.21};
+  const double LoF_Vend[4] = {560.8, 552.4, 546.2, 565.5};
+  double ToF_Shift_Vendlg[4], ToF_Shift_Vendhg[4];
+  
+  for(int i = 0; i <=3; i++){
+    ToF_Shift_Vendlg[i] = LoF_Vend[i] / c_light - GPeak_Vendlg[i];
+    ToF_Shift_Vendhg[i] = LoF_Vend[i] / c_light - GPeak_Vendhg[i];
+  }  
+
+  double FC_Qmax = -100;
   unsigned int FC_mult = FC->AnodeNumber.size();
-  if(FC_mult==1){
-    int anode = FC->AnodeNumber[0];
-    double Time_FC = FC->Time[0];
-
-    Vendeta->SetAnodeNumber(anode);
-    Vendeta->BuildPhysicalEvent();
-
-    unsigned int Vendeta_mult = Vendeta->DetectorNumber.size();
-    for(unsigned int i=0; i<Vendeta_mult; i++){
-      int DetNbr          = Vendeta->DetectorNumber[i];
-      double Time_Vendeta = Vendeta->Time[i];
-      double Rdet         = Vendeta->GetDistanceFromTarget(DetNbr);
-      TVector3 DetPos     = Vendeta->GetVectorDetectorPosition(DetNbr);
-
-      double DT = Time_Vendeta - Time_FC;
-
-      double DeltaTheta = atan(63.5/Rdet);
-      double Theta_Vendeta = DetPos.Theta();
-      double Theta_random = ra.Uniform(Theta_Vendeta-DeltaTheta,Theta_Vendeta+DeltaTheta);
-      //cout << DT << " " << Rdet << endl;
-      //neutron->SetTimeOfFlight(DT*1e-9/(Rdet*1e-3));
-      neutron->SetTimeOfFlight(DT*1e-9/(0.55));
-      double En = neutron->GetEnergy();
-
-      // Filling output tree
-      Tof.push_back(DT);
-      ELab.push_back(En);
-      ThetaLab.push_back(Theta_random);
-      Q1.push_back(Vendeta->Q1[i]);
-      Q2.push_back(Vendeta->Q2[i]);
-      HG_status.push_back(Vendeta->isHG[i]);
-    }
-  }
+  if(FC_mult==1 ){
+    FC_Qmax = FC->Qmax[0];
+
+    if(FC_Qmax > 400){
+
+      int anode = FC->AnodeNumber[0];
+      double Time_FC = FC->Time[0];
+
+      Vendeta->SetAnodeNumber(anode);
+      Vendeta->BuildPhysicalEvent();
+      FC->BuildPhysicalEvent();
+
+      // VENDETA LG 
+      unsigned int Vendeta_LG_mult = Vendeta->LG_DetectorNumber.size();
+      for(unsigned int i=0; i<Vendeta_LG_mult; i++){
+        
+        int DetNbr          = Vendeta->LG_DetectorNumber[i];
+        double Time_Vendeta = Vendeta->LG_Time[i];
+        double Rdet         = Vendeta->GetDistanceFromTarget(DetNbr);
+        TVector3 DetPos     = Vendeta->GetVectorDetectorPosition(DetNbr);
+
+        double DT = Time_Vendeta - Time_FC + ToF_Shift_Vendlg[DetNbr-1];
+        
+        if(DT>0){
+          
+          double DeltaTheta = atan(63.5/Rdet);
+          double Theta_Vendeta = DetPos.Theta();
+          double Theta_random = ra.Uniform(Theta_Vendeta-DeltaTheta,Theta_Vendeta+DeltaTheta);
+          //cout << DT << " " << Rdet << endl;
+          //neutron->SetTimeOfFlight(DT*1e-9/(Rdet*1e-3));
+          //neutron->SetTimeOfFlight(DT*1e-9/(0.55));
+          neutron->SetBeta( (LoF_Vend[DetNbr-1] / DT )/ c_light); 
+
+          double En = neutron->GetEnergy();
+
+          // Filling output tree
+          LG_Tof.push_back(DT);
+          LG_ID.push_back(DetNbr);
+          LG_ELab.push_back(En);
+          LG_ThetaLab.push_back(Theta_random);
+          LG_Q1.push_back(Vendeta->LG_Q1[i]);
+          LG_Q2.push_back(Vendeta->LG_Q2[i]);
+          LG_Qmax.push_back(Vendeta->LG_Qmax[i]);
+        }
+      }
+
+      // VENDETA HG 
+      unsigned int Vendeta_HG_mult = Vendeta->HG_DetectorNumber.size();
+      for(unsigned int i=0; i<Vendeta_HG_mult; i++){
+        int DetNbr          = Vendeta->HG_DetectorNumber[i];
+        double Time_Vendeta = Vendeta->HG_Time[i];
+        double Rdet         = Vendeta->GetDistanceFromTarget(DetNbr);
+        TVector3 DetPos     = Vendeta->GetVectorDetectorPosition(DetNbr);
+
+        double DT = Time_Vendeta - Time_FC + ToF_Shift_Vendhg[DetNbr-1];
+
+        if(DT>0){
+          double DeltaTheta = atan(63.5/Rdet);
+          double Theta_Vendeta = DetPos.Theta();
+          double Theta_random = ra.Uniform(Theta_Vendeta-DeltaTheta,Theta_Vendeta+DeltaTheta);
+          //cout << DT << " " << Rdet << endl;
+          //neutron->SetTimeOfFlight(DT*1e-9/(Rdet*1e-3));
+          //neutron->SetTimeOfFlight(DT*1e-9/(0.55));
+          neutron->SetBeta( (LoF_Vend[DetNbr-1] / DT )/ c_light); 
+          double En = neutron->GetEnergy();
+
+          // Filling output tree
+          HG_ID.push_back(DetNbr);
+          HG_Tof.push_back(DT);
+          HG_ELab.push_back(En);
+          HG_ThetaLab.push_back(Theta_random);
+          HG_Q1.push_back(Vendeta->HG_Q1[i]);
+          HG_Q2.push_back(Vendeta->HG_Q2[i]);
+          HG_Qmax.push_back(Vendeta->HG_Qmax[i]);
+        }
+      }
+
+      //Process coincidences signals in VENDETA LG / HG
+
+      if(HG_Tof.size() > 0 && LG_Tof.size() > 0 ){
+        for(int j = 0; j < LG_Tof.size();j++){
+          for(int k = 0; k < HG_Tof.size(); k++){
+            if(abs(HG_Tof[k]-LG_Tof[j]) < 2 && HG_ID[k] == LG_ID[j]){
+              if( HG_Q2[k]>120000){
+                /* HG_ID[k] = */
+                HG_Tof[k] = - 100000;
+                HG_ELab[k] = - 100000;
+                HG_ThetaLab[k] = - 100000;
+                HG_Q1[k] = - 100000;
+                HG_Q2[k] = - 100000;
+                HG_Qmax[k] = - 100000;
+              }  
+              else if( HG_Q2[k]<120000){
+                /* HG_ID[k] = */
+                LG_Tof[k] = - 100000;
+                LG_ELab[k] = - 100000;
+                LG_ThetaLab[k] = - 100000;
+                LG_Q1[k] = - 100000;
+                LG_Q2[k] = - 100000;
+
+              }
+            }
+          }
+        }
+      } // if LG && HG
+    } // if FC_Qmax > 400
+  }// if FC = 1
 
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch(){
-  RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
-  RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
-  RootOutput::getInstance()->GetTree()->Branch("Tof",&Tof);
-  RootOutput::getInstance()->GetTree()->Branch("Q1",&Q1);
-  RootOutput::getInstance()->GetTree()->Branch("Q2",&Q2);
-  RootOutput::getInstance()->GetTree()->Branch("HG_status",&HG_status);
+  RootOutput::getInstance()->GetTree()->Branch("LG_ID",&LG_ID);
+  RootOutput::getInstance()->GetTree()->Branch("LG_ThetaLab",&LG_ThetaLab);
+  RootOutput::getInstance()->GetTree()->Branch("LG_ELab",&LG_ELab);
+  RootOutput::getInstance()->GetTree()->Branch("LG_Tof",&LG_Tof);
+  RootOutput::getInstance()->GetTree()->Branch("LG_Q1",&LG_Q1);
+  RootOutput::getInstance()->GetTree()->Branch("LG_Q2",&LG_Q2);
+  RootOutput::getInstance()->GetTree()->Branch("LG_Qmax",&LG_Qmax);
+
+  RootOutput::getInstance()->GetTree()->Branch("HG_ID",&HG_ID);
+  RootOutput::getInstance()->GetTree()->Branch("HG_ThetaLab",&HG_ThetaLab);
+  RootOutput::getInstance()->GetTree()->Branch("HG_ELab",&HG_ELab);
+  RootOutput::getInstance()->GetTree()->Branch("HG_Tof",&HG_Tof);
+  RootOutput::getInstance()->GetTree()->Branch("HG_Q1",&HG_Q1);
+  RootOutput::getInstance()->GetTree()->Branch("HG_Q2",&HG_Q2);
+  RootOutput::getInstance()->GetTree()->Branch("HG_Qmax",&HG_Qmax);
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::ReInitValue(){
-  ThetaLab.clear();
-  ELab.clear();
-  Tof.clear();
-  Q1.clear();
-  Q2.clear();
-  HG_status.clear();
+  LG_ThetaLab.clear();
+  LG_ELab.clear();
+  LG_Tof.clear();
+  LG_ID.clear();
+  LG_Q1.clear();
+  LG_Q2.clear();
+  LG_Qmax.clear();
+
+  HG_ThetaLab.clear();
+  HG_ELab.clear();
+  HG_Tof.clear();
+  HG_ID.clear();
+  HG_Q1.clear();
+  HG_Q2.clear();
+  HG_Qmax.clear();
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////

Projects/Vendeta/Analysis.h

@@ -43,12 +43,22 @@ class Analysis: public NPL::VAnalysis{
     static NPL::VAnalysis* Construct();
 
   private:
-    vector<double> ThetaLab;
-    vector<double> ELab;
-    vector<double> Tof;
-    vector<double> Q1;
-    vector<double> Q2;
-    vector<bool>   HG_status;
+    vector<double> LG_ID;
+    vector<double> LG_ThetaLab;
+    vector<double> LG_ELab;
+    vector<double> LG_Tof;
+    vector<double> LG_Q1;
+    vector<double> LG_Q2;
+    vector<double> LG_Qmax;
+
+    vector<double> HG_ID;
+    vector<double> HG_ThetaLab;
+    vector<double> HG_ELab;
+    vector<double> HG_Tof;
+    vector<double> HG_Q1;
+    vector<double> HG_Q2;
+    vector<double> HG_Qmax;
+    
 
   private:
     TVendetaPhysics* Vendeta;