From 8aa31fd96f44a310721c2f3fecf39b42f6141969 Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Mon, 25 Feb 2019 11:21:11 +0100
Subject: [PATCH] * Progress on MUGAST   - Strip ordering is now working

---
 NPLib/Detectors/Mugast/MugastMap.h        | 513 +++++++++++-----------
 NPLib/Detectors/Mugast/TMugastPhysics.cxx |  90 ++--
 NPLib/Detectors/Mugast/TMugastSpectra.cxx |   6 +-
 3 files changed, 305 insertions(+), 304 deletions(-)

diff --git a/NPLib/Detectors/Mugast/MugastMap.h b/NPLib/Detectors/Mugast/MugastMap.h
index e7bfc3798..0999bbf95 100644
--- a/NPLib/Detectors/Mugast/MugastMap.h
+++ b/NPLib/Detectors/Mugast/MugastMap.h
@@ -538,136 +538,134 @@ namespace MUGAST_MAP{
 
 
       const int AnnularX[128]{
-      // unconnected strip //
-      65,//1
-      66,//2
-      67,//3
-      68,//4
-      69,//5
-      70,//6
-      71,//7
-      72,//8
-      73,//9
-      74,//10
-      75,//11
-      76,//12
-      77,//13
-      78,//14
-      79,//15
-      80,//16
-      81,//17
-      82,//18
-      83,//19
-      84,//20
-      85,//21
-      86,//22
-      87,//23
-      88,//24
-      89,//25
-      90,//26
-      91,//27
-      92,//28
-      93,//29
-      94,//30
-      95,//31
-      96,//32
-      97,//33
-      98,//34
-      99,//35
-      100,//36
-      101,//37
-      102,//38
-      103,//39
-      104,//40
-      105,//41
-      106,//42
-      107,//43
-      108,//44
-      109,//45
-      110,//46
-      111,//47
-      112,//48
-      113,//49
-      114,//50
-      115,//51
-      116,//52
-      117,//53
-      118,//54
-      119,//55
-      120,//56
-      121,//57
-      122,//58
-      123,//59
-      124,//60
-      125,//61
-      126,//62
-      127,//63
-      128,//64
-      // Connected strip
-      1,  // 65  
-      2,  // 66  
-      4,  // 67  
-      3,  // 68  
-      6,  // 69  
-      5,  // 70  
-      8,  // 71  
-      7,  // 72  
-      10, // 73  
-      9,  // 74  
-      12, // 75  
-      11, // 76  
-      14, // 77  
-      13, // 78  
-      16, // 79  
-      15, // 80  
-      18, // 81  
-      17, // 82  
-      20, // 83  
-      19, // 84  
-      22, // 85  
-      21, // 86  
-      24, // 87  
-      23, // 88  
-      26, // 89  
-      25, // 90  
-      28, // 91  
-      27, // 92  
-      30, // 93  
-      29, // 94  
-      31, // 95  
-      32, // 96  
-      64, // 97  
-      63, // 98  
-      62, // 99  
-      61, // 100 
-      60, // 101 
-      59, // 102 
-      58, // 103 
-      57, // 104 
-      56, // 105 
-      55, // 106 
-      54, // 107 
-      53, // 108 
-      52, // 109 
-      51, // 110 
-      49, // 112 
-      50, // 111 
-      48, // 113 
-      47, // 114 
-      46, // 115 
-      45, // 116 
-      44, // 117 
-      43, // 118 
-      42, // 119 
-      41, // 120 
-      40, // 121 
-      39, // 122 
-      38, // 123 
-      37, // 124 
-      36, // 125 
-      35, // 126 
-      34, // 127 
-      33 // 128 
+      128,//1
+      128,//2
+      128,//3
+      128,//4
+      128,//5
+      128,//6
+      128,//7
+      128,//8
+      128,//9
+      128,//10
+      128,//11
+      128,//12
+      128,//13
+      128,//14
+      128,//15
+      128,//16
+      128,//17
+      128,//18
+      128,//19
+      128,//20
+      128,//21
+      128,//22
+      128,//23
+      128,//24
+      128,//25
+      128,//26
+      128,//27
+      128,//28
+      128,//29
+      32,//30
+      128,//31
+      31,//32
+      128,//33
+      30,//34
+      128,//35
+      28,//36
+      29,//37
+      26,//38
+      27,//39
+      24,//40
+      25,//41
+      22,//42
+      23,//43
+      20,//44
+      21,//45
+      18,//46
+      19,//47
+      16,//48
+      17,//49
+      14,//50
+      15,//51
+      12,//52
+      13,//53
+      10,//54
+      11,//55
+      8,//56
+      9,//57
+      6,//58
+      7,//59
+      4,//60
+      5,//61
+      2,//62
+      3,//63
+      1,//64
+      33,  // 65  
+      34,  // 66  
+      34,  // 67  ////// WARNING ? ERROR IN MAPPING FILE?
+      128,  // 68  
+      36,  // 69  
+      128,  // 70  
+      35,  // 71  
+      128,  // 72  
+      38, // 73  
+      128,  // 74  
+      37, // 75  
+      128, // 76  
+      40, // 77  
+      128, // 78  
+      39, // 79  
+      128, // 80  
+      42, // 81  
+      128, // 82  
+      41, // 83  
+      128, // 84  
+      44, // 85  
+      128, // 86  
+      43, // 87  
+      128, // 88  
+      46, // 89  
+      128, // 90  
+      45, // 91  
+      128, // 92  
+      48, // 93  
+      128, // 94  
+      47, // 95  
+      128, // 96  
+      50, // 97  
+      128, // 98  
+      49, // 99  
+      128, // 100 
+      52, // 101 
+      128, // 102 
+      51, // 103 
+      128, // 104 
+      54, // 105 
+      128, // 106 
+      53, // 107 
+      128, // 108 
+      56, // 109 
+      128, // 110 
+      55, // 111 
+      128,// 112
+      58, // 113 
+      128, // 114 
+      57, // 115 
+      128, // 116 
+      60, // 117 
+      128, // 118 
+      59, // 119 
+      128, // 120 
+      62, // 121 
+      128, // 122 
+      61, // 123 
+      128, // 124 
+      63, // 125 
+      128, // 126 
+      64, // 127 
+      128 // 128 
       };
 
   ///////
@@ -677,135 +675,134 @@ namespace MUGAST_MAP{
    //16 wedges phi
    const int AnnularY[128]={
      // only strip connected values //
-     43,// 1
-     41,// 2
-     42,// 3
-     44,// 4
-     46,// 5
-     48,// 6
-     50,// 7
-     52,// 8
-     54,// 9
-     56,// 10
-     55,// 11
-     53,// 12
-     51,// 13
-     49,// 14
-     47,// 15
-     45,// 16
-     // end of connected strip //
-     81,//17
-     82,//18
-     83,//19
-     84,//20
-     85,//21
-     86,//22
-     87,//23
-     88,//24
-     89,//25
-     90,//26
-     91,//27
-     92,//28
-     93,//29
-     94,//30
-     95,//31
-     96,//32
-     97,//33
-     98,//34
-     99,//35
-     100,//36
-     101,//37
-     102,//38
-     103,//39
-     104,//40
-     105,//41
-     106,//42
-     107,//43
-     108,//44
-     109,//45
-     110,//46
-     111,//47
-     112,//48
-     113,//49
-     114,//50
-     115,//51
-     116,//52
-     117,//53
-     118,//54
-     119,//55
-     120,//56
-     121,//57
-     122,//58
-     123,//59
-     124,//60
-     125,//61
-     126,//62
-     127,//63
+     128,// 1
+     128,// 2
+     128,// 3
+     128,// 4
+     128,// 5
+     128,// 6
+     128,// 7
+     128,// 8
+     128,// 9
+     128,// 10
+     128,// 11
+     128,// 12
+     128,// 13
+     10,// 14
+     128,// 15
+     11,// 16
+     128,//17
+     9,//18
+     128,//19
+     12,//20
+     128,//21
+     8,//22
+     128,//23
+     13,//24
+     128,//25
+     7,//26
+     128,//27
+     14,//28
+     128,//29
+     128,//30
+     128,//31
+     128,//32
+     128,//33
+     128,//34
+     128,//35
+     6,//36
+     128,//37
+     15,//38
+     128,//39
+     5,//40
+     128,//41
+     16,//42
+     128,//43
+     44,//44
+     128,//45
+     1,//46
+     128,//47
+     3,//48
+     128,//49
+     2,//50
+     128,//51
+     128,//52
+     128,//53
+     128,//54
+     128,//55
+     128,//56
+     128,//57
+     128,//58
+     128,//59
+     128,//60
+     128,//61
+     128,//62
+     128,//63
      128,//64
-     1,//65
-     2,//66
-     3,//67
-     4,//68
-     5,//69
-     6,//70
-     7,//71
-     8,//72
-     9,//73
-     10,//74
-     11,//75
-     12,//76
-     13,//77
-     14,//78
-     15,//79
-     16,//80
-     17,//81
-     18,//82
-     19,//83
-     20,//84
-     21,//85
-     22,//86
-     23,//87
-     24,//88
-     25,//89
-     26,//90
-     27,//91
-     28,//92
-     29,//93
-     30,//94
-     31,//95
-     32,//96
-     33,//97
-     34,//98
-     35,//99
-     36,//100
-     37,//101
-     38,//102
-     39,//103
-     40,//104
-     57,  //105
-     58,  //106
-     59,  //107
-     60,  //108
-     61,  //109
-     62,  //110
-     63,  //111
-     64,  //112
-     65,  //113
-     66,  //114
-     67,  //115
-     68,  //116
-     69,  //117
-     70,  //118
-     71,  //119
-     72,  //120
-     73,  //121
-     74,  //122
-     75,  //123
-     76,  //124
-     77,  //125
-     78,  //126
-     79,  //127
-     80//128
+     128,//65
+     128,//66
+     128,//67
+     128,//68
+     128,//69
+     128,//70
+     128,//71
+     128,//72
+     128,//73
+     128,//74
+     128,//75
+     128,//76
+     128,//77
+     128,//78
+     128,//79
+     128,//80
+     128,//81
+     128,//82
+     128,//83
+     128,//84
+     128,//85
+     128,//86
+     128,//87
+     128,//88
+     128,//89
+     128,//90
+     128,//91
+     128,//92
+     128,//93
+     128,//94
+     128,//95
+     128,//96
+     128,//97
+     128,//98
+     128,//99
+     128,//100
+     128,//101
+     128,//102
+     128,//103
+     128,//104
+     128,//105
+     128,//106
+     128,//107
+     128,//108
+     128,//109
+     128,//110
+     128,//111
+     128,//112
+     128,//113
+     128,//114
+     128,//115
+     128,//116
+     128,//117
+     128,//118
+     128,//119
+     128,//120
+     128,//121
+     128,//122
+     128,//123
+     128,//124
+     128,//125
+     128,//126
+     128,//127
+     128 //128
    };
    
 }  
diff --git a/NPLib/Detectors/Mugast/TMugastPhysics.cxx b/NPLib/Detectors/Mugast/TMugastPhysics.cxx
index 9eff408b8..b1ee97b2b 100644
--- a/NPLib/Detectors/Mugast/TMugastPhysics.cxx
+++ b/NPLib/Detectors/Mugast/TMugastPhysics.cxx
@@ -90,17 +90,18 @@ void TMugastPhysics::PreTreat() {
           m_EventData->GetDSSDXEStripNbr(i))) {
       double EX = fDSSD_X_E(m_EventData, i);
       if (EX > m_DSSD_X_E_Threshold)
-        m_PreTreatedData->SetDSSDXE(MG_NOCHANGE,
-        m_EventData->GetDSSDXEDetectorNbr(i),
-        m_EventData->GetDSSDXEStripNbr(i), EX);
+      m_PreTreatedData->SetDSSDXE(MG_NOCHANGE,
+          m_EventData->GetDSSDXEDetectorNbr(i),
+          m_EventData->GetDSSDXEStripNbr(i), EX);
     }
   }
+
   //   T
   for (unsigned int i = 0; i < DSSDX_TMult; ++i) {
     if (IsValidChannel(0, m_EventData->GetDSSDXTDetectorNbr(i),
           m_EventData->GetDSSDXTStripNbr(i)))
       m_PreTreatedData->SetDSSDXT(MG_NOCHANGE,
-m_EventData->GetDSSDXTDetectorNbr(i),
+          m_EventData->GetDSSDXTDetectorNbr(i),
           m_EventData->GetDSSDXTStripNbr(i),
           fDSSD_X_T(m_EventData, i));
   }
@@ -114,7 +115,7 @@ m_EventData->GetDSSDXTDetectorNbr(i),
       double EY = fDSSD_Y_E(m_EventData, i);
       if (EY > m_DSSD_Y_E_Threshold)
         m_PreTreatedData->SetDSSDYE(MG_NOCHANGE,
-m_EventData->GetDSSDYEDetectorNbr(i),
+            m_EventData->GetDSSDYEDetectorNbr(i),
             m_EventData->GetDSSDYEStripNbr(i), EY);
     }
   }
@@ -124,7 +125,7 @@ m_EventData->GetDSSDYEDetectorNbr(i),
     if (IsValidChannel(1, m_EventData->GetDSSDYTDetectorNbr(i),
           m_EventData->GetDSSDYTStripNbr(i)))
       m_PreTreatedData->SetDSSDYT(MG_NOCHANGE,
-m_EventData->GetDSSDYTDetectorNbr(i),
+          m_EventData->GetDSSDYTDetectorNbr(i),
           m_EventData->GetDSSDYTStripNbr(i),
           fDSSD_Y_T(m_EventData, i));
   }
@@ -348,9 +349,9 @@ vector<TVector2> TMugastPhysics::Match_X_Y() {
 ////////////////////////////////////////////////////////////////////////////
 bool TMugastPhysics::IsValidChannel(const int& DetectorType,
     const int& telescope, const int& channel) {
-  if (DetectorType == 0)
+  if (DetectorType == 0){
     return *(m_XChannelStatus[m_DetectorNumberIndex[telescope] - 1].begin() + channel - 1);
-
+  }
   else if (DetectorType == 1)
     return *(m_YChannelStatus[m_DetectorNumberIndex[telescope] - 1].begin() + channel - 1);
 
@@ -468,18 +469,18 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
   string Type; 
 
   for (unsigned int i = 0; i < blocks.size(); i++) {
-   
+
     if (blocks[i]->HasTokenList(cart)) {
       if (NPOptionManager::getInstance()->GetVerboseLevel())
-      
-      Type = blocks[i]->GetMainValue();
+
+        Type = blocks[i]->GetMainValue();
       cout << endl << "////  Mugast Telescope " << Type << " " << i + 1 << endl;
-      
+
       int detectorNbr = blocks[i]->GetInt("DetectorNumber");
       if(Type=="Square") DetectorType[detectorNbr]=MG_SQUARE;
       else if(Type=="Trapezoid") DetectorType[detectorNbr]=MG_TRAPEZE;
       else if(Type=="Annular") DetectorType[detectorNbr]=MG_ANNULAR;
-   
+
       det = i+1;
       m_DetectorNumberIndex[detectorNbr]=det;
       TVector3 A = blocks[i]->GetTVector3("X1_Y1", "mm");
@@ -488,10 +489,10 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
       TVector3 D = blocks[i]->GetTVector3("X128_Y128", "mm");
       AddTelescope(A, B, C, D);
     }
-    
+
     else if (blocks[i]->HasTokenList(annular)) {
       if (NPOptionManager::getInstance()->GetVerboseLevel())
-      Type = blocks[i]->GetMainValue();
+        Type = blocks[i]->GetMainValue();
       cout << endl << "////  Mugast Telescope " << Type << " " << i + 1 << endl;
       int detectorNbr = blocks[i]->GetInt("DetectorNumber");
       if(Type=="Square") DetectorType[detectorNbr]=MG_SQUARE;
@@ -500,7 +501,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
       if(Type!="Annular"){
         cout << "ERROR: Using Mugast Annular Token for Square or Trapezoid detector " << endl;
         exit(1);
-        }
+      }
 
       det = i+1;
       m_DetectorNumberIndex[detectorNbr]=det;
@@ -511,13 +512,13 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
 
     else if (blocks[i]->HasTokenList(sphe)) {
       if (NPOptionManager::getInstance()->GetVerboseLevel())
-      Type = blocks[i]->GetMainValue();
+        Type = blocks[i]->GetMainValue();
       cout << endl << "////  Mugast Telescope " << Type << " " << i + 1 << endl;
       int detectorNbr = blocks[i]->GetInt("DetectorNumber");
       if(Type=="Square") DetectorType[detectorNbr]=MG_SQUARE;
       else if(Type=="Trapezoid") DetectorType[detectorNbr]=MG_TRAPEZE;
       else if(Type=="Annular") DetectorType[detectorNbr]=MG_ANNULAR;
-  
+
       det = i+1;
       m_DetectorNumberIndex[detectorNbr]=det;
       double         Theta = blocks[i]->GetDouble("THETA", "deg");
@@ -533,6 +534,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
         << endl;
       exit(1);
     }
+
   }
 
   InitializeStandardParameter();
@@ -541,7 +543,7 @@ void TMugastPhysics::ReadConfiguration(NPL::InputParser parser) {
   std::ofstream shapeFile(".MugastShape");
   for(auto& it:DetectorType){
     shapeFile << it.first << " " << it.second << endl;
-    }
+  }
   shapeFile.close();
   //ReadAnalysisConfig();
 }
@@ -571,8 +573,8 @@ void TMugastPhysics::WriteSpectra() {
 
 ///////////////////////////////////////////////////////////////////////////
 map<string, TH1*> TMugastPhysics::GetSpectra() {
- if(m_Spectra)
-  return m_Spectra->GetMapHisto();
+  if(m_Spectra)
+    return m_Spectra->GetMapHisto();
   else {
     map<string, TH1*> empty;
     return empty;
@@ -776,22 +778,22 @@ void TMugastPhysics::AddTelescope(TVector3 C_Center) {
   vector< vector< double > >   OneStripPositionX   ;
   vector< vector< double > >   OneStripPositionY   ;
   vector< vector< double > >   OneStripPositionZ   ;
-  
- /* The logic behind the strip numbering of S1 in NPTOOL: 
- The number of rings goes from 1->64, the number of sectors goes from 1->16 
- (4 per quadrant). There's a redundancy in the fact that 1->64 already contain 
- the information about the quadrant and the majority of these positions are 
- indeed not physical. Example: 
- A hit combining Ring 17 (first ring in Quadrant 2) and 
- Sector 4 (last sector in Quadrant 1) is not possible due to physical mismatch 
- of the detector frontside-backside layout. 
- The possible (allowed hits) are R(1-16)S(1-4), R(17-32)S(5-8), R(33-49)S(9-12),
- R(50-64)S(13-16). 
- The three loops however takes all the possible combintation that an analysis
- can produce. This works perfectly for cases where the detector does not have 
- "Quadrants", e.g. S3 type. For the S1 an extra condition is added to flag the
- non physical hit combinations. 
- */
+
+  /* The logic behind the strip numbering of S1 in NPTOOL: 
+     The number of rings goes from 1->64, the number of sectors goes from 1->16 
+     (4 per quadrant). There's a redundancy in the fact that 1->64 already contain 
+     the information about the quadrant and the majority of these positions are 
+     indeed not physical. Example: 
+     A hit combining Ring 17 (first ring in Quadrant 2) and 
+     Sector 4 (last sector in Quadrant 1) is not possible due to physical mismatch 
+     of the detector frontside-backside layout. 
+     The possible (allowed hits) are R(1-16)S(1-4), R(17-32)S(5-8), R(33-49)S(9-12),
+     R(50-64)S(13-16). 
+     The three loops however takes all the possible combintation that an analysis
+     can produce. This works perfectly for cases where the detector does not have 
+     "Quadrants", e.g. S3 type. For the S1 an extra condition is added to flag the
+     non physical hit combinations. 
+     */
   for(int iQuad = 0 ; iQuad < NumberOfQuadrant ; iQuad++){
     for(int iRing = 0 ; iRing < NumberofRing; iRing++){
       lineX.clear() ;
@@ -799,25 +801,27 @@ void TMugastPhysics::AddTelescope(TVector3 C_Center) {
       lineZ.clear() ;
       for(int iSector = 0 ; iSector < NumberofSector ; iSector++){
         //Build vector
-        StripCenter = Strip_1_1;
-        StripCenter = TVector3(R_Min+(iRing+0.5)*StripPitchRing,0, Z);
+        StripCenter = TVector3(C_Center.X()+R_Min+(iRing+0.5)*StripPitchRing,C_Center.Y(), Z);
         StripCenter.RotateZ( ( Phi_0 - (iSector+0.5)*StripPitchSector ) *M_PI/180.);
 
         // if the hit is not "allowed" (see comment above) use a default value
         if ( (iRing+(iQuad*NumberofRing))/NumberofSector != (iSector/NumberOfQuadrant) ) 
-          StripCenter.SetXYZ(-100,-100, Z-100);
-        
+          StripCenter.SetXYZ(-1000,-1000, Z-1000);
+
         lineX.push_back( StripCenter.X() );// these vectors will contain 16x4 = 64 elements
         lineY.push_back( StripCenter.Y() );
         lineZ.push_back( StripCenter.Z() );
       }
-
       OneStripPositionX.push_back(lineX);
       OneStripPositionY.push_back(lineY);
       OneStripPositionZ.push_back(lineZ);
-
     }
   }
+  vector<double> defaultLine;
+  defaultLine.resize(128,-1000);
+  OneStripPositionX.resize(128,defaultLine);
+  OneStripPositionY.resize(128,defaultLine);
+  OneStripPositionZ.resize(128,defaultLine);
   m_StripPositionX.push_back( OneStripPositionX ) ;
   m_StripPositionY.push_back( OneStripPositionY ) ;
   m_StripPositionZ.push_back( OneStripPositionZ ) ;
diff --git a/NPLib/Detectors/Mugast/TMugastSpectra.cxx b/NPLib/Detectors/Mugast/TMugastSpectra.cxx
index eb06f82ad..903f4fe5b 100644
--- a/NPLib/Detectors/Mugast/TMugastSpectra.cxx
+++ b/NPLib/Detectors/Mugast/TMugastSpectra.cxx
@@ -354,13 +354,13 @@ void TMugastSpectra::FillPhysicsSpectra(TMugastPhysics* Physics){
     Theta = Theta/deg;
     FillSpectra(family,name,Theta,Physics->DSSD_E[i]);
 
-//    // STRX_E_CAL
+    // STRX_E_CAL
 //    name = "MG"+NPL::itoa( Physics->TelescopeNumber[i])+"_XY_COR";
-//    FillSpectra(family,name,Physics->DSSD_EX[i],Physics->DSSD_EY[i]);
+ //   FillSpectra(family,name,Physics->DSSD_E[i],Physics->DSSD_EY[i]);
 
 
     // Fill only for particle stopped in the first stage
-    if(Physics->SecondLayer_E[i]<0 && Physics->SecondLayer_E[i]<0){
+    if(Physics->SecondLayer_E[i]<0 ){
       // E-TOF:
       name = "MG_E_TOF";
       FillSpectra(family,name,Physics->DSSD_E[i],Physics->DSSD_T[i]);
-- 
GitLab