From f62a591d87f9470fe7c80663ac46256dc2927baa Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Thu, 23 Dec 2010 17:23:59 +0000
Subject: [PATCH] * Finish stripping scheme for GaspardTrackerAnnular    + bug
 was in the correspondance between strip number and coordinates    +
 excitation energy is now well reconstructed

* Add Destroy method in the singleton class NPOptionManager
   + it is used in the Gaspard analysis code

* Change default ROOT file name in general analysis macros
   + mySimul.root is now myResult.root

* Add DumpStrippingScheme method in GaspardTracker
   + it allows to dump the correspondance between strips and coordinates
     for any detector part of Gaspard
   + for the moment indexes are hardcoded for annular detector (should be
     improved)
---
 Inputs/EventGenerator/132Sndp.reaction      |  2 +-
 Inputs/EventGenerator/alpha.source          |  4 +--
 NPAnalysis/Gaspard/Makefile                 |  1 -
 NPAnalysis/Gaspard/src/Analysis.cc          | 22 +++++-------
 NPAnalysis/macros/ControlSimu.C             |  2 +-
 NPAnalysis/macros/GeometricalEfficiency.C   |  2 +-
 NPAnalysis/macros/TestAngleReconstruction.C |  2 +-
 NPLib/GASPARD/GaspardTracker.cxx            | 16 +++++++++
 NPLib/GASPARD/GaspardTracker.h              |  1 +
 NPLib/GASPARD/GaspardTrackerAnnular.cxx     | 38 ++++++++++++---------
 NPLib/GASPARD/GaspardTrackerAnnular.h       | 15 ++++----
 NPLib/Tools/NPOptionManager.cxx             | 12 +++++++
 NPSimulation/src/GaspardScorers.cc          | 16 ++++-----
 NPSimulation/src/GaspardTrackerAnnular.cc   |  2 --
 14 files changed, 78 insertions(+), 57 deletions(-)

diff --git a/Inputs/EventGenerator/132Sndp.reaction b/Inputs/EventGenerator/132Sndp.reaction
index ffe7e8e53..6507dc83e 100644
--- a/Inputs/EventGenerator/132Sndp.reaction
+++ b/Inputs/EventGenerator/132Sndp.reaction
@@ -7,7 +7,7 @@ Transfert
 	Target= 2H
 	Light= 1H
 	Heavy= 133Sn
-	ExcitationEnergy= 10.0
+	ExcitationEnergy= 0.0
 	BeamEnergy= 1320
 	BeamEnergySpread= 0
 	SigmaX= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index c8597743e..31cee696d 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -6,8 +6,8 @@
 Isotropic
 	EnergyLow= 0
 	EnergyHigh= 5
-	HalfOpenAngleMin= 0
-	HalfOpenAngleMax= 45
+	HalfOpenAngleMin= 158
+	HalfOpenAngleMax= 180 
 	x0= 0	
 	y0= 0	
 	z0= 0	
diff --git a/NPAnalysis/Gaspard/Makefile b/NPAnalysis/Gaspard/Makefile
index 39d6b7505..215794b28 100644
--- a/NPAnalysis/Gaspard/Makefile
+++ b/NPAnalysis/Gaspard/Makefile
@@ -8,4 +8,3 @@ clean:
 distclean:
 	make clean -C ./src
 	rm Analysis
-	rm AnalysisNew
diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc
index 196d2f3e4..0e4eb1ba6 100644
--- a/NPAnalysis/Gaspard/src/Analysis.cc
+++ b/NPAnalysis/Gaspard/src/Analysis.cc
@@ -5,13 +5,13 @@ using namespace std;
 
 int main(int argc,char** argv)
 {	
-
-  NPOptionManager* myOptionManager  = NPOptionManager::getInstance(argc,argv)   ;
-  string reactionfileName           = myOptionManager->GetReactionFilePath()    ;
-	string detectorfileName 		      = myOptionManager->GetDetectorFilePath()	  ;
-	string calibrationfileName        = myOptionManager->GetCalibrationFilePath()	;
-	string runToReadfileName          = myOptionManager->GetRunToReadFilePath()   ;
-	string OutputfileName             = myOptionManager->GetOutputFilePath()      ;
+   // command line parsing
+   NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+   string reactionfileName          = myOptionManager->GetReactionFilePath();
+	string detectorfileName          = myOptionManager->GetDetectorFilePath();
+	string calibrationfileName       = myOptionManager->GetCalibrationFilePath();
+	string runToReadfileName         = myOptionManager->GetRunToReadFilePath();
+	string OutputfileName            = myOptionManager->GetOutputFilePath();
 
    // Instantiate RootInput and RootOutput singleton classes
    RootInput:: getInstance(runToReadfileName);
@@ -89,7 +89,6 @@ int main(int argc,char** argv)
 
       // Get total energy
       double E = GPDTrack->GetEnergyDeposit();
-//      cout << i << "  " << E << endl;
 
       // if there is a hit in the detector array, treat it.
       double Theta, ThetaStrip, angle, ThetaCM;
@@ -101,15 +100,11 @@ int main(int argc,char** argv)
          ThetaCM = initCond->GetICEmittedAngleThetaCM(0) * deg;
 
          // Get exact scattering angle from TInteractionCoordinates object
-//         Theta = interCoord->GetDetectedAngleTheta(0) * deg;
-//         cout << interCoord << endl;
 //         interCoord->Dump();
 //         cout << i << " mult: " << interCoord->GetDetectedMultiplicity() << endl;
          DetecX = interCoord->GetDetectedPositionX(0);
          DetecY = interCoord->GetDetectedPositionY(0);
          DetecZ = interCoord->GetDetectedPositionZ(0);
-//         cout << "Detected position :" << endl;
-//         cout << "\t" << DetecX << "  " << DetecY << "  " << DetecZ << endl;
          TVector3 Detec(DetecX, DetecY, DetecZ);
 
          // Get interaction position in detector
@@ -119,11 +114,9 @@ int main(int argc,char** argv)
          // Get beam interaction coordinates on target (from initial condition)
          XTarget = initCond->GetICPositionX(0);
          YTarget = initCond->GetICPositionY(0);
-//         cout << XTarget << "  " << YTarget << endl;
          BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
          BeamPhi   = initCond->GetICIncidentAnglePhi(0)*deg;
          TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
-//         cout << BeamDirection.X() << "  " << BeamDirection.Y() << "  " << BeamDirection.Z() << endl;
 
          // Hit direction taking into account beam position on target
          TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
@@ -182,6 +175,7 @@ int main(int argc,char** argv)
    // delete singleton classes
    RootOutput::getInstance()->Destroy();
    RootInput::getInstance()->Destroy();
+   NPOptionManager::getInstance()->Destroy();
 
    return 0;
 }
diff --git a/NPAnalysis/macros/ControlSimu.C b/NPAnalysis/macros/ControlSimu.C
index eec454239..f8f9cdb4b 100644
--- a/NPAnalysis/macros/ControlSimu.C
+++ b/NPAnalysis/macros/ControlSimu.C
@@ -40,7 +40,7 @@
 #include "TInitialConditions.h"
 #include "TInteractionCoordinates.h"
 
-void ControlSimu(const char * fname = "mySimul")
+void ControlSimu(const char * fname = "myResult")
 {
    // Open output ROOT file from NPTool simulation run
    TString path = gSystem->Getenv("NPTOOL");
diff --git a/NPAnalysis/macros/GeometricalEfficiency.C b/NPAnalysis/macros/GeometricalEfficiency.C
index fed847d8c..8c8023369 100644
--- a/NPAnalysis/macros/GeometricalEfficiency.C
+++ b/NPAnalysis/macros/GeometricalEfficiency.C
@@ -42,7 +42,7 @@
 using namespace std ;
 
 
-void GeometricalEfficiency(const char * fname = "mySimul")
+void GeometricalEfficiency(const char * fname = "myResult")
 {
    // Open output ROOT file from NPTool simulation run
    TString path = gSystem->Getenv("NPTOOL");
diff --git a/NPAnalysis/macros/TestAngleReconstruction.C b/NPAnalysis/macros/TestAngleReconstruction.C
index b29fc482f..862995944 100644
--- a/NPAnalysis/macros/TestAngleReconstruction.C
+++ b/NPAnalysis/macros/TestAngleReconstruction.C
@@ -42,7 +42,7 @@
 #include "TInitialConditions.h"
 #include "TInteractionCoordinates.h"
 
-void TestAngleReconstruction(const char * fname = "mySimul")
+void TestAngleReconstruction(const char * fname = "myResult")
 {
    // Open output ROOT file from NPTool simulation run
    TString path = gSystem->Getenv("NPTOOL");
diff --git a/NPLib/GASPARD/GaspardTracker.cxx b/NPLib/GASPARD/GaspardTracker.cxx
index 66a35eddd..e0e0d5ac0 100644
--- a/NPLib/GASPARD/GaspardTracker.cxx
+++ b/NPLib/GASPARD/GaspardTracker.cxx
@@ -258,6 +258,22 @@ void GaspardTracker::DumpModulesMap()
 
 
 
+void GaspardTracker::DumpStrippingScheme(int moduleNumber)
+{
+   cout << "GaspardTracker::DumpStrippingScheme()" << endl;
+
+   for (int i = 1; i < 65; i++) {   // front part
+      for (int j = 1; j < 65; j++) {   // back part
+         cout << "strips X, Y: " << i << "  " << j << "\t--->\t (X,Y,Z) mm: "
+              << m_ModulesMap[moduleNumber]->GetStripPositionX(moduleNumber, i, j) << "\t"
+              << m_ModulesMap[moduleNumber]->GetStripPositionY(moduleNumber, i, j) << "\t"
+              << m_ModulesMap[moduleNumber]->GetStripPositionZ(moduleNumber, i, j) << endl;
+      }
+   }
+}
+
+
+
 double GaspardTracker::GetEnergyDeposit()		
 { 
    if (m_EventPhysics->GetEventMultiplicity() > 0) {
diff --git a/NPLib/GASPARD/GaspardTracker.h b/NPLib/GASPARD/GaspardTracker.h
index 49b1a3332..c617ef482 100644
--- a/NPLib/GASPARD/GaspardTracker.h
+++ b/NPLib/GASPARD/GaspardTracker.h
@@ -84,6 +84,7 @@ public:
 
 public:
    void DumpModulesMap();
+   void DumpStrippingScheme(int moduleNumber);
 
 private:
    map<int, GaspardTrackerModule*>	m_ModulesMap;
diff --git a/NPLib/GASPARD/GaspardTrackerAnnular.cxx b/NPLib/GASPARD/GaspardTrackerAnnular.cxx
index b16d8f90c..69a5f8cbc 100644
--- a/NPLib/GASPARD/GaspardTrackerAnnular.cxx
+++ b/NPLib/GASPARD/GaspardTrackerAnnular.cxx
@@ -20,7 +20,10 @@ GaspardTrackerAnnular::GaspardTrackerAnnular(map<int, GaspardTrackerModule*> &Mo
 	  m_EventPhysics(EventPhysics),
 	  m_EventData(0),
 	  m_PreTreatData(new TGaspardTrackerData),
-	  m_NumberOfModule(0)
+	  m_NumberOfModule(0),
+     m_NumberOfStripsTheta(16),
+     m_NumberOfStripsPhi(16),
+     m_NumberOfQuadrants(4)
 {
 }
 
@@ -247,14 +250,9 @@ void GaspardTrackerAnnular::AddModule(double zpos, double rmin, double rmax)
 {
    m_NumberOfModule++;
 
-   // Characteristics
-   int NbPhiStrips     = 16;
-   int NbThetaStrips   = 16;
-   int NbQuadrant = 4;
-
    // Theta & phi strips pitch
-   double thetaPitch = (rmax - rmin) / NbThetaStrips;
-   double phiPitch   = 2*M_PI / (NbQuadrant*NbThetaStrips);
+   double thetaPitch = (rmax - rmin) / m_NumberOfStripsTheta;
+   double phiPitch   = 2*M_PI / (m_NumberOfQuadrants*m_NumberOfStripsPhi);
 
    // Buffer object to fill Position Array
    vector<double> lineX;
@@ -266,25 +264,31 @@ void GaspardTrackerAnnular::AddModule(double zpos, double rmin, double rmax)
    vector< vector< double > >   OneModuleStripPositionZ;
 
    // loop on theta strips
-   for (int i = 0; i < NbThetaStrips*NbQuadrant; i++) {
+   for (int i = 0; i < m_NumberOfStripsTheta*m_NumberOfQuadrants; i++) {
       lineX.clear();
       lineY.clear();
       lineZ.clear();
 
       // center of theta strip
-      double r = rmin + thetaPitch/2 + thetaPitch*(i % NbThetaStrips);
-//      cout << i << "  " << i%NbThetaStrips << "   " << r << endl;
+      double r = rmin + thetaPitch/2 + thetaPitch*(i % m_NumberOfStripsTheta);
+         
+      // current quandrant
+      int quadrant = i / m_NumberOfStripsTheta;
 
       // loop on phi strips
-      for (int j = 0; j < NbPhiStrips*NbQuadrant; j++) {
+      for (int j = 0; j < m_NumberOfStripsPhi*m_NumberOfQuadrants; j++) {
+         // initialize x and y
+         double x = 0;
+         double y = 0;
          // center of phi strips
          double phi = phiPitch/2 + phiPitch*j;
+         if (phi > quadrant*M_PI/2 && phi < (quadrant+1)*M_PI/2) {
+            // calculate x and y projections
+            x = r * cos(phi);
+            y = r * sin(phi);
+            if (zpos < 0) y *= -1;
+         }
 
-         // calculate x and y projections
-         double x = r * cos(phi);
-         double y = r * sin(phi);
-
-//         cout << i << "  " << j << "  " << r << "  " << phi*180/M_PI << "   " << x << "   " << y << endl;
          // fill lineX,Y,Z
          lineX.push_back(x);
          lineY.push_back(y);
diff --git a/NPLib/GASPARD/GaspardTrackerAnnular.h b/NPLib/GASPARD/GaspardTrackerAnnular.h
index fa423e0f3..d3661d841 100644
--- a/NPLib/GASPARD/GaspardTrackerAnnular.h
+++ b/NPLib/GASPARD/GaspardTrackerAnnular.h
@@ -61,9 +61,9 @@ public:
    void AddModule(double zpos, double rmin, double rmax);
 
    // Getters to retrieve the (X,Y,Z) coordinates of a pixel defined by strips (X,Y)
-   double GetStripPositionX(int N ,int X ,int Y)        { return m_StripPositionX[N-1-m_index["Annular"]][X-1][Y-1]; }
-   double GetStripPositionY(int N ,int X ,int Y)        { return m_StripPositionY[N-1-m_index["Annular"]][X-1][Y-1]; }
-   double GetStripPositionZ(int N ,int X ,int Y)        { return m_StripPositionZ[N-1-m_index["Annular"]][X-1][Y-1]; }
+   double GetStripPositionX(int N ,int X ,int Y)        {return m_StripPositionX[N-1-m_index["Annular"]][X-1][Y-1]; }
+   double GetStripPositionY(int N ,int X ,int Y)        {return m_StripPositionY[N-1-m_index["Annular"]][X-1][Y-1]; }
+   double GetStripPositionZ(int N ,int X ,int Y)        {return m_StripPositionZ[N-1-m_index["Annular"]][X-1][Y-1]; }
    double GetNumberOfModule()                           { return m_NumberOfModule; }
 
 private:
@@ -77,12 +77,9 @@ private:
    //////////////////////////////
    // Geometry and stip number //
    //////////////////////////////
-   double m_FirstStageBaseLarge; // mm
-   double m_FirstStageHeight;    // mm
-   int    m_NumberOfStripsX;
-   int    m_NumberOfStripsY;
-   double m_StripPitchX;      // mm
-   double m_StripPitchY;      // mm
+   int m_NumberOfStripsTheta;
+   int m_NumberOfStripsPhi;
+   int m_NumberOfQuadrants;
 };
 
 #endif
diff --git a/NPLib/Tools/NPOptionManager.cxx b/NPLib/Tools/NPOptionManager.cxx
index dd36a146a..4d1262eca 100644
--- a/NPLib/Tools/NPOptionManager.cxx
+++ b/NPLib/Tools/NPOptionManager.cxx
@@ -88,3 +88,15 @@ void NPOptionManager::DisplayHelp()
     cout << "\t --run -R <arg>\t \t \t \t \t \t \t Set arg as the run to read file list" << endl  ;
     cout << endl << endl ;  
   }
+
+
+
+void NPOptionManager::Destroy()
+{
+   if (instance != 0) {
+      delete instance;
+      instance = 0;
+   }
+}
+
+
diff --git a/NPSimulation/src/GaspardScorers.cc b/NPSimulation/src/GaspardScorers.cc
index 4ebfe571f..cb5c4a575 100644
--- a/NPSimulation/src/GaspardScorers.cc
+++ b/NPSimulation/src/GaspardScorers.cc
@@ -679,12 +679,13 @@ G4bool GPDScorerFirstStageFrontStripAnnular::ProcessHits(G4Step* aStep, G4Toucha
    if (dummy < 0 && fabs(dummy) < 1e-6) dummy *= -1;
    G4double ThetaStripNumber = floor(dummy / ThetaStripPitch);
    ThetaStripNumber += PhiQuadrantNumber * GPDANNULAR::NbThetaStrips;
-
-//      G4cout << "POS: " << POS << G4endl;
-//      G4cout << "r, phi " << r << "  " << phi << G4endl;
-//      G4cout << "PhiWidth, PhiQuadrantNumber " << PhiWidth << "  " << PhiQuadrantNumber << G4endl;
-//      G4cout << "ThetaStripPitch, ThetaStripNumber, dummy " << ThetaStripPitch << "  " << ThetaStripNumber << "  " << dummy << G4endl;
-
+   ThetaStripNumber++;
+/*
+      G4cout << "POS: " << POS << G4endl;
+      G4cout << "r, phi " << r << "  " << phi << G4endl;
+      G4cout << "PhiWidth, PhiQuadrantNumber " << PhiWidth << "  " << PhiQuadrantNumber << G4endl;
+      G4cout << "ThetaStripPitch, ThetaStripNumber, dummy " << ThetaStripPitch << "  " << ThetaStripNumber << "  " << dummy << G4endl;
+*/
    if (ThetaStripNumber < 1e-6) {
     /*  G4cout << "POS: " << POS << G4endl;
       G4cout << "r, phi " << r << "  " << phi << G4endl;
@@ -695,7 +696,6 @@ G4bool GPDScorerFirstStageFrontStripAnnular::ProcessHits(G4Step* aStep, G4Toucha
    G4double edep = aStep->GetTotalEnergyDeposit();
    if (edep < 100*keV) return FALSE;
    G4int  index =  aStep->GetTrack()->GetTrackID();
-   ThetaStripNumber++;
    EvtMap->set(DetNbr + index, ThetaStripNumber);
 
    return TRUE;
@@ -763,6 +763,7 @@ G4bool GPDScorerFirstStageBackStripAnnular::ProcessHits(G4Step* aStep, G4Touchab
    // Interstrip should be taken into account here. To be done
    G4double PhiWidth = 360. / (GPDANNULAR::NbPhiStrips * GPDANNULAR::NbThetaQuadrant);
    G4double PhiStripNumber = floor(phi / PhiWidth);
+   PhiStripNumber++;
 
 //   G4cout << POS << G4endl;
 //   G4cout << "phi " << phi << " PhiWidth  " << PhiWidth << "  PhiStripNumber " << PhiStripNumber << G4endl;
@@ -770,7 +771,6 @@ G4bool GPDScorerFirstStageBackStripAnnular::ProcessHits(G4Step* aStep, G4Touchab
    G4double edep = aStep->GetTotalEnergyDeposit();
    if (edep < 100*keV) return FALSE;
    G4int  index =  aStep->GetTrack()->GetTrackID();
-   PhiStripNumber++;
    EvtMap->set(DetNbr + index, PhiStripNumber);
 
    return TRUE;
diff --git a/NPSimulation/src/GaspardTrackerAnnular.cc b/NPSimulation/src/GaspardTrackerAnnular.cc
index 8fc455a41..1bcb62bc0 100644
--- a/NPSimulation/src/GaspardTrackerAnnular.cc
+++ b/NPSimulation/src/GaspardTrackerAnnular.cc
@@ -246,8 +246,6 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber   ,
 
       // Silicon detector itself
       G4ThreeVector  positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);
-      G4cout << "position en z de l'annulaire " << MMpos.z() * mm << G4endl;
-      G4cout << "position en z de l'annulaire " << Silicon_PosZ * mm << G4endl;
 
       G4Tubs* solidSilicon = new G4Tubs("solidSilicon", 
                                          FirstStageRmin,
-- 
GitLab