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

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/Samurai/TSamuraiHodoscopePhysics.h

@@ -156,6 +156,11 @@ class TSamuraiHodoscopePhysics : public TObject, public NPL::VDetector {
     
     // clear the pre-treated object
     void ClearPreTreatedData()   {m_PreTreatedData->Clear();}
+    
+  
+    vector <double> GetCharge() {return Charge;}
+    vector <double> GetTime() {return Time;}
+    vector<int> GetID() {return ID;}
 
     // give and external TSamuraiHodoscopeData object to TSamuraiHodoscopePhysics. 
     // needed for online analysis for example

Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx

@@ -6,8 +6,10 @@ double devE(const double* parameter){
 
   //Beam spot offset
   TVector3 offset(parameter[0],parameter[1],parameter[2]);
-  unsigned int size = pos.size();
 
+  //Other variable initilizations
+  unsigned int size = pos.size();
+  TVector3 MugastNormal;
   double dE,Theta;
   TVector3 dir;
 
@@ -25,24 +27,53 @@ double devE(const double* parameter){
     //Particle path vector
     dir=*(pos[i])-offset;
 
-    //Detected energy, and angle of particle leaving target
-    double Theta= dir.Angle(TVector3(0,0,1));
+    //Define normal vector for the MG# of detection
+    //[[[[ UPDATE WITH NEW MG POSITIONS FROM SURVEY ]]]]
+    switch(detnum[i]){
+      case 1:
+        MugastNormal.SetXYZ(-0.453915, +0.455463, -0.765842);
+        break;
+      case 2:
+	MugastNormal.SetXYZ(-0.642828, +0.000000, -0.766010);
+        break;
+      case 3:
+	MugastNormal.SetXYZ(-0.454594, -0.450670, -0.768271);
+        break;
+      case 4:
+	MugastNormal.SetXYZ(-0.002437, -0.638751, -0.769409);
+        break;
+      case 5:
+	MugastNormal.SetXYZ(+0.452429, -0.454575, -0.767248);
+        break;
+      case 7:
+	MugastNormal.SetXYZ(+0.443072, +0.443265, -0.779232);
+        break;
+      default:
+	cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
+        return 1; // Exit code
+    }
+
+    //Angle leaving target, angle entering MUGAST & energy deposited in MUGAST
+    double ThetaTarget = dir.Angle(TVector3(0,0,1));
+    double ThetaMugast = dir.Angle(MugastNormal);
     double Energy = energy[i];
 
-    //NOTE!!! Not calucualting energy loss in Al???
+    //Energy loss in Al
+    Energy=Al.EvaluateInitialEnergy(
+		Energy,                      //energy after Al 
+		0.4*micrometer,              //thickness of Al
+		ThetaMugast);                //angle impinging on MUGAST
     //Energy loss in target
     Energy=CD2.EvaluateInitialEnergy(
-		    Energy,                      //energy after leaving target
-		    0.5*parameter[4]*micrometer, //pass through half target
-		    Theta);                      //angle leaving target
+		Energy,                      //energy after leaving target
+		0.5*parameter[3]*micrometer, //pass through half target
+		ThetaTarget);                //angle leaving target
 
     //Final value of Ex
-    double Ex = reaction.ReconstructRelativistic(Energy,Theta);
+    double Ex = reaction.ReconstructRelativistic(Energy,ThetaTarget);
     
-    //Fill histogram with ERROR in Ex!
-    //h->Fill(Ex-refE); 
-    h->Fill(Ex); 
-
+    //Fill histograms with Ex
+    h->Fill(Ex);
     switch(detnum[i]){
       case 1:
         h1->Fill(Ex); 
@@ -63,70 +94,76 @@ double devE(const double* parameter){
         h7->Fill(Ex); 
         break;
       default:
-        cout << "ERROR! Invalid detnum: " << detnum[i] << " @" << i << endl;
+        cout << "ERROR:: Invalid DetNum " << detnum[i] << " at event " << i << endl;
         return 1; // Exit code
     }
-
   }
   //End loop over events
 
   //Write vals to screen
   cout << "Mean: " << h->GetMean() 
        << "\t StdDev: " << h->GetStdDev() 
-       << "\t Thickness??: " << parameter[4] 
+       << "\t Thickness: " << parameter[3] << " um" 
        << endl;
 
   //Draw histogram(s)
   h->Draw();
   if(flagDraw){ InitiliseCanvas(); }
 
+  /*
+  cout << pow(h->GetMean()-refE,2)  << " + " 
+       << pow(0.1*h->GetStdDev(),2) << " = " 
+       << pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) 
+       << endl;
+  */
+
   //Adapt the metric as needed
   return sqrt( pow(h->GetMean()-refE,2) + pow(0.1*h->GetStdDev(),2) );
 }
 ////////////////////////////////////////////////////////////////////////////////
 void MinimizeBeamSpot(){
 
-  // Read data in
+  //Read data in
   LoadFile();
 
-  // Start with beam (0,0,0) and 4.7um 0.5mg/c2 target
-  double parameter[4] = {0.0, 0.0, 0.0, 4.7};   
+  //Start with beam (0,0,0) and 4.76um 0.5mg/c2 target
+  double parameter[4] = {0.0, 0.0, 0.0, 4.76};   
   devE(parameter);
 
-  // Function with 4 parameter XYZ and Target thickness
+  //Function with 4 parameter XYZ and Target thickness
   auto func = ROOT::Math::Functor(&devE,4);
  
-  // Minimizer
+  //Minimizer
   auto minim = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad"); 
 
   minim->SetPrintLevel(0);
   minim->SetPrecision(1e-10); 
 
-  // Set minimizer function
+  //Set minimizer function
   minim->SetFunction(func);
 
-  // Assign variable limits
+  //Assign variable limits
   minim->SetLimitedVariable(0,"X",parameter[0],0.01,-10,10);
   minim->SetLimitedVariable(1,"Y",parameter[1],0.01,-10,10);
   minim->SetLimitedVariable(2,"Z",parameter[2],0.01,-5,5);
-  minim->SetLimitedVariable(3,"T",parameter[3],0.01, 4.7-3.0, 4.7+3.0);
+  minim->SetLimitedVariable(3,"T",parameter[3],0.01,4.76*0.5,4.76*1.5);
 
-  // Don't draw iterations of minimizer
+  //Don't draw iterations of minimizer
   flagDraw = 0;
 
-  // Shrink it, babeyyy
+  //Shrink it, babeyyy
   minim->Minimize(); 
   
-  // Draw minimal value
+  //Draw minimal value
   flagDraw = 1;
 
-  // Pull values from minimizer
+  //Pull values from minimizer
   const double* x = minim->X();
   cout << "========================================" << endl;
   cout << "\t\tX =" << x[0] << endl;
   cout << "\t\tY =" << x[1] << endl;
   cout << "\t\tZ =" << x[2] << endl;
   cout << "\t\tT =" << x[3] << endl;
-  cout << "Minimum: " << devE(x) << endl;
+  devE(x);
   cout << "========================================" << endl;
 }

Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h

@@ -14,8 +14,8 @@ NPL::EnergyLoss Al("proton_Al.G4table","G4Table",100);
 using namespace std;
 
 bool flagDraw = 0;
-static auto h = new TH1D("h","h", 80,-1.,1.);
-static auto h1 = new TH1D("h1","h1", 40,-1.,1.);
+static auto h = new TH1D("h","All MG#'s", 60,-1.,1.);
+static auto h1 = new TH1D("h1","Individual MG#'s", 40,-1.,1.);
 static auto h2 = new TH1D("h2","h2", 40,-1.,1.);
 static auto h3 = new TH1D("h3","h3", 40,-1.,1.);
 static auto h4 = new TH1D("h4","h4", 40,-1.,1.);
@@ -125,6 +125,7 @@ void InitiliseCanvas(){
   // ----- ALL -----
   canv->cd(2);
   h->SetStats(0);
+  h->SetLineColor(kBlack);
   h->GetXaxis()->SetTitle("Ex [MeV]");
   h->GetYaxis()->SetTitle("Counts");
   h->Draw();