From c8795231674138e19e87ff8eef5c4ad64d2d5de8 Mon Sep 17 00:00:00 2001
From: Dmitry Shchukin <dmitry.shchukin@cern.ch>
Date: Tue, 11 Oct 2022 13:49:24 +0200
Subject: [PATCH 1/3] initial working state

---
 mkcode.cc |  30 +++-
 simu.h    | 506 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 532 insertions(+), 4 deletions(-)
 create mode 100644 simu.h

diff --git a/mkcode.cc b/mkcode.cc
index 011aef5..442596a 100644
--- a/mkcode.cc
+++ b/mkcode.cc
@@ -33,6 +33,9 @@ using namespace CS;
 #include "tracking_new.h"
 #endif
 
+//MC simpletracking
+#include "simu.h"
+
 int procev(RecoEvent &e);
 void endrun();
 
@@ -249,7 +252,8 @@ int main(int argc, char *argv[])
       return 1;
     }
   
-    const int Nevents = 100000000;
+    //const int Nevents = 100000000;
+    const int Nevents = 2000;
     int NevMCRead = 0;
 
     // event loop, read event by event
@@ -544,6 +548,25 @@ int procev(RecoEvent &e0)
 
       if(mymom < MomentumLowerCut) ipass = 0;
       if(mymom > MomentumUpperCut) ipass = 0;
+    } 
+    else {
+      if(MCFileFormatVersion >= 8) {
+        //Init_Geometry_MC2022B();
+        if(ievproc == 0) {
+          Reset_Calib();
+          Init_Calib_2022B_MM(1664928923);
+          LoadGeometryFromCalo();
+        }
+        DoDigitization(e);
+        man.Reset();
+        man.AddMicromega(e.MM1);    
+        man.AddMicromega(e.MM2);    
+        man.AddMicromega(e.MM3);
+        man.AddMicromega(e.MM4);
+        
+        track = simple_tracking_4MM_2016(e);
+        if(track) mymom = track.momentum;
+      }
     }
 #endif
 #ifdef TRACKINGTOOLS
@@ -976,7 +999,7 @@ int procev(RecoEvent &e0)
 
   double mychi2 = 0.;
   //mychi2 = calcShowerChi2(e, 2, 2, 3, 3);
-  mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
+  //mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
   if(ECALChi2UpperCut < 101. && mychi2 > ECALChi2UpperCut) ipass = 0;
   if(ECALMaxCell > -1) { // check max cell
     if(e.mc || ECALMaxCell <= 29) { // set needed cell numbers from cuts
@@ -1099,10 +1122,9 @@ int procev(RecoEvent &e0)
               << " Veto = " << veto
               << std::endl;
   }
-
+  
   momentum->Fill(mymom);   //track.momentum
 
-
   ehplot->Fill(ecal, hcal);
   srdplot->Fill(srd);
   srdplot0->Fill(e.SRD[0].energy/MeV);
diff --git a/simu.h b/simu.h
new file mode 100644
index 0000000..e687110
--- /dev/null
+++ b/simu.h
@@ -0,0 +1,506 @@
+#pragma once
+
+//root
+#include "TVector3.h"
+#include "TRandom.h"
+#include "TMath.h"
+
+//na64
+#include "p348reco.h"
+#include "mm.h"
+#include "conddb.h"
+
+//GLOBALS
+const double MCsignalCalibrationX = 0.05;
+const double MCsignalCalibrationY = 0.15;
+const int MCsizeX = 5;
+const int MCsizeY = 12;
+const bool simulatenoise = true;
+//method to compute charge spread
+enum DigiMethod {STD,GAUS};
+const DigiMethod MM_chargespread = GAUS;
+//Maximum threshold for MM
+const double MM_threshold = 22e-6; //22 keV, taken from https://www.sciencedirect.com/science/article/pii/S1875389212018019
+const double GEM_threshold = 22e-6; // same, waiting input for Michael Hoesgen
+//gem resolution in mm
+double GEMyres = 0.060;
+double GEMxres = 0.060;
+//maximal distance before the merging
+double MergingDistance = 1.75; //mm
+//decide if to replicate fake hits
+const bool ReplicateHits = true;
+
+//geometry definition
+void Init_Geometry_MC2017()
+{
+ Reset_Calib();
+ //imported from standard example1
+ ECAL_pos.pos = TVector3(     0,   0,     0);  
+ MAGD_pos.pos = TVector3(     0,   0, -13975.1);
+ MAGU_pos.pos = TVector3(     0,   0, -17975.1);
+ MAG_field = 1.7;
+ MM1_pos.pos  = TVector3(315.,0.,-19375);
+ MM2_pos.pos  = TVector3(315.,0.,-18575);
+ MM3_pos.pos  = TVector3(0.,0.,-1840);
+ MM4_pos.pos  = TVector3(0.,0.,-640); 
+ //all same angle for simplicity
+ MM1_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM3_pos.Angle = M_PI/4;
+ MM4_pos.Angle = M_PI/4;
+ //same pedestal calibration as in 2017
+ // timestamp 1506895200 is `2017 Oct 2 0:00 CEST` (last midnight before the end of beam time)
+ Init_Calib_2017_MM(1506895200);
+ NominalBeamEnergy = 100.;
+}
+
+void Init_Geometry_MC2018()
+{
+ Reset_Calib();
+ //imported from standard example1
+ ECAL_pos.pos = TVector3(     0,   0,     0);  
+ MAGU_pos.pos = TVector3(     0,   0, -19530.1);
+ //MAGD_pos.pos = TVector3(     0,   0, -15530.1);
+ //adjusted down position
+ MAGD_pos.pos = TVector3(     0,   0, -16600.1);
+ MAG_field = 1.85;
+ //position of MM
+ MM3_pos.pos  = TVector3(0.,      0., -17593.7);
+ MM4_pos.pos  = TVector3(0.,      0., -17226.7); 
+ MM5_pos.pos  = TVector3(-203.33, 0.,  287.35);
+ MM6_pos.pos  = TVector3(-206.667,0.,  387.35);
+ //all same angle for simplicity
+ MM1_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM3_pos.Angle = M_PI/4;
+ MM5_pos.Angle = M_PI/4;
+ MM6_pos.Angle = M_PI/4;
+ //same pedestal calibration as in 2018
+ // timestamp 1529445600 is `2018 June 20 0:00 CEST` (last midnight before the end of beam time)
+ Init_Calib_2018_MM(1529445600);
+ NominalBeamEnergy = 150.;
+}
+
+//placeholder for dummy geometry
+void Init_Geometry_MC2018_invis()
+{
+ Reset_Calib();
+ //imported from standard example1
+ ECAL_pos.pos = TVector3(     0,   0,     0);  
+ MAGU_pos.pos = TVector3(     0,   0, -19530.1);
+ //MAGD_pos.pos = TVector3(     0,   0, -15530.1);
+ //adjusted down position
+ MAGD_pos.pos = TVector3(     0,   0, -16600.1);
+ MAG_field = 1.85;
+ //position of MM
+ MM1_pos.pos  = TVector3(0.,      0., -16930);
+ MM2_pos.pos  = TVector3(0.,      0., -16130); 
+ MM3_pos.pos  = TVector3(-319.,      0.,  854);
+ MM4_pos.pos  = TVector3(-329.,      0.,  1048); 
+ MM5_pos.pos  = TVector3(-381, 0.,  3309);
+ MM6_pos.pos  = TVector3(-397.667,0.,  3521);
+ //all same angle for simplicity
+ MM1_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM3_pos.Angle = M_PI/4;
+ MM5_pos.Angle = M_PI/4;
+ MM6_pos.Angle = M_PI/4;
+ //same pedestal calibration as in 2018
+ // timestamp 1529445600 is `2018 June 20 0:00 CEST` (last midnight before the end of beam time)
+ Init_Calib_2018_MM(1529445600);
+ NominalBeamEnergy = 150.;
+}
+
+void Init_Geometry_MC2022B()
+{
+ //from conddb.h
+ ECAL_pos.pos      = TVector3(-379.9,  8.5,   24887.1);
+ MAGU_pos.pos = TVector3(     0,   0,   6350.+506.5);
+ MAGD_pos.pos = TVector3(     0,   0,   6350.+506.5+4000.); 
+ MAG_field = 1.99;
+ MM1_pos.pos  = TVector3(       1.2,     -1.3,     4377.4);
+ MM2_pos.pos  = TVector3(     -13.0,     -2.5,     5279.9);
+ MM3_pos.pos  = TVector3(    -313.0,     -3.4,    23181.7);
+ MM4_pos.pos  = TVector3(    -364.9,     -5.4,    24395.6);
+ MM1_pos.Angle = M_PI/4;
+ MM2_pos.Angle = M_PI/4;
+ MM3_pos.Angle = M_PI/4;
+ MM4_pos.Angle = M_PI/4;
+ NominalBeamEnergy = 100.;
+}
+
+//collection of function for Micromegas response simulation
+
+void MCRealHit(vector<MChit>& mchit,const double& EnergyThreshold)
+{
+ sort(mchit.begin(),mchit.end(),MChitByID);
+ vector<MChit> new_hit;
+ for(unsigned int i(0);i<mchit.size();)
+  {
+   int id = mchit[i].id;
+   double energy = mchit[i].energy;
+   int particle = mchit[i].particle;
+   double deposit(0),xpos(0),ypos(0),zpos(0);
+   while(mchit[i].id == id)
+    {
+     deposit += mchit[i].deposit;
+     xpos += mchit[i].deposit*mchit[i].xpos;
+     ypos += mchit[i].deposit*mchit[i].ypos;
+     zpos += mchit[i].deposit*mchit[i].zpos;
+     ++i;
+    }
+   //weighted mean of the hit
+   xpos /= deposit;
+   ypos /= deposit;
+   zpos /= deposit;
+   if(energy > EnergyThreshold) //deposit would be better, but would need Geant4 revision
+    new_hit.push_back(MChit (xpos,ypos,zpos,deposit,particle,id,energy));
+  }
+ //all computation is over, substitute vector with one ordered by energy
+ sort(new_hit.begin(),new_hit.end(),MChitByDeposit);
+ mchit = new_hit;
+  
+}//end MCRealHit
+
+void SimulateMultiplexing(MicromegaPlane& plane,unsigned int strip,const double& signal)
+{
+ //const unsigned int detectedsignal=round(gRandom->Poisson(signal)); 
+ const unsigned int detectedsignal=round(signal);
+ vector<unsigned int>& PlaneStrips = plane.strips;
+ const int reverse_index = plane.calib->MM_reverse_map.at(strip);
+ const vector<int> strips = plane.calib->MM_map.at(reverse_index);
+ for(unsigned int j(0);j<strips.size();++j)
+  {
+   //limit ADC count to 4095
+   if(PlaneStrips[strips[j]] + detectedsignal < 4095)
+    {
+     PlaneStrips[strips[j]] += detectedsignal;
+    }
+   else
+    {
+     PlaneStrips[strips[j]] = 4095;
+    }
+  }
+} ///end SimulateMultiplexing
+
+//method of computation of strip deposit in strip as function of distance from hit point
+double ComputeCharge(const double& signal,const double& dist2center,const unsigned int ClusSize)
+{
+ if( MM_chargespread == STD )
+  {
+   return signal / pow(2,dist2center+1);
+  }
+ else if( MM_chargespread == GAUS )
+  {
+   const double chargespread = ClusSize/4.;
+   return signal*TMath::Gaus(dist2center,0,chargespread);
+  }
+ else
+  {
+   cout << "WARNING: NO VALID METHOD FOR CHARGE SPREAD OVER STRIPS SELECTED" << endl;
+   exit(-1);
+  }
+}
+
+
+void SimulateHit(MicromegaPlane& plane,const double center,const unsigned int meansize, const double& signal)
+{
+ //define size of the cluster
+ unsigned int size(0); 
+ while (size < 2)
+  {
+   size = gRandom->Poisson(meansize);
+  }
+ //in case cluster is asymmetric, choose a direction
+ int direction(1);
+ if(size%2 == 0)
+  {
+   direction = 1-2*gRandom->Integer(2);
+  }   
+ //substract center
+ int strip = round(center);
+ if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size()))
+  return;
+ double dist2center = fabs(center - strip);
+ double stripsignal = ComputeCharge(signal,dist2center,size);
+ SimulateMultiplexing(plane,strip,stripsignal);
+ unsigned int counter = 1;
+ while(counter < size)
+  {
+   int distance = counter/2 +1;
+   strip = round(center+direction*distance);
+   dist2center = fabs(center - strip);
+   stripsignal = ComputeCharge(signal,dist2center,size);
+   if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size())){++counter;continue;}
+   //multiplex the strip
+   SimulateMultiplexing(plane,strip,stripsignal);
+   ++counter;
+   strip = round(center-direction*distance);
+   dist2center = fabs(center - strip);
+   stripsignal = ComputeCharge(signal,dist2center,size);
+   if((strip < 0) || (((unsigned int)strip + 1) > plane.strips.size()) || counter >= (size)){++counter;continue;}
+   ++counter;
+   SimulateMultiplexing(plane,strip,stripsignal);
+  }  
+}
+
+
+
+void MCDoHitAccumulation(Micromega& MM,vector<MChit> MMtruth, const int verbose = 0)
+{
+ //compute strip hit
+ const double angle = MM.geom->Angle; 
+ const double Kx = cos(angle);
+ const double Ky = sin(angle);
+ const double xcenter = MM.calib->xcalib.MM_Nstrips/2.0 - 0.5;
+ const double ycenter = MM.calib->ycalib.MM_Nstrips/2.0 - 0.5;
+ for(unsigned int i(0);i<MMtruth.size();++i)
+  {
+   double xpos = MMtruth[i].xpos - MM.geom->pos.X();
+   double ypos = MMtruth[i].ypos - MM.geom->pos.Y();
+   const double signal = MMtruth[i].deposit;
+   //add resolution and convert to strip size
+   static const double StripRes = 1/sqrt(12);
+   xpos /= MM.calib->xcalib.MM_strip_step;
+   ypos /= MM.calib->ycalib.MM_strip_step;
+   //xpos = gRandom->Gaus(xpos,StripRes);
+   //ypos = gRandom->Gaus(ypos,StripRes);
+   /*find best strips for the hit and signal 
+     invert the function in MMCalculatePos to account 45 degree conversion*/      
+   const double xstrip = (-Kx*xpos - Ky*ypos + xcenter);
+   const double ystrip = (Ky*xpos - Kx*ypos + ycenter);
+   const double signalX = signal*MCsignalCalibrationX;
+   const double signalY = signal*MCsignalCalibrationY;
+   if(verbose > 0)
+    cout << " hit with xcenter: " << xcenter << " and ycenter: " << ycenter << " and signalX: " << signalX << " and signalY: " << signalY << endl;
+   if(xstrip >= 0 && xstrip < MM.calib->xcalib.MM_Nstrips)
+    {
+     SimulateHit(MM.xplane,xstrip,MCsizeX,signalX);
+    }
+   if(ystrip >= 0 && ystrip < MM.calib->ycalib.MM_Nstrips)
+    {
+     SimulateHit(MM.yplane,ystrip,MCsizeY,signalY);
+    }
+  }
+} //end MCDoHitAccumulation
+
+//function to simulate noise
+void SimulatePlaneNoise(MicromegaPlane& plane,const vector<double>& sigma)
+{
+ //loop over channels
+ vector<unsigned int>& PlaneStrips = plane.strips;
+ for(unsigned int i(0);i<plane.calib->MM_Nchannels;++i)
+  {
+   const double noise = gRandom->Gaus(0,sigma[i]);
+   const vector<int> strips = plane.calib->MM_map.at(i);
+   const double currentsignal = PlaneStrips[strips[0]];
+   if(currentsignal + noise < 2*sigma[i]) //sigma used in data
+    {
+     continue;
+    }
+   else
+    {     
+     for(unsigned int j(0);j<strips.size();++j)
+      {
+       PlaneStrips[strips[j]] += noise;    
+      }
+    }
+  }
+}
+
+void SimulateNoise(Micromega& mm)
+{
+ SimulatePlaneNoise(mm.xplane,mm.calib->xcalib.sigma);
+ SimulatePlaneNoise(mm.yplane,mm.calib->ycalib.sigma);
+}
+
+void DoSmearing(vector<MChit>& mmhits,const double& xres,const double& yres)
+{
+ for(unsigned int i = 0; i < mmhits.size(); i++)
+  {
+   mmhits[i].xpos = gRandom->Gaus(mmhits[i].xpos,xres);
+   mmhits[i].ypos = gRandom->Gaus(mmhits[i].ypos,yres);
+  }
+}
+
+vector<MChit> DoMerging(vector<MChit> mmhits,const double maxD = MergingDistance)
+{
+ vector<MChit> mmhitsmerged;
+ for(auto p1 = mmhits.begin(); p1 != mmhits.end(); ++p1)
+  {
+   MChit newhit = *p1;
+   for(auto p2 = p1+1; p2 != mmhits.end(); )
+    {
+    if(abs(newhit.xpos-(*p2).xpos) < maxD && abs(newhit.ypos-(*p2).ypos) < maxD)
+     {
+      //DISTANCE IS TOO SHORT, DO MERGING
+      //TODO: half of distance is taken, use weighted mean instead?
+      newhit.xpos = 0.5*(newhit.xpos + (*p2).xpos);
+      newhit.ypos = 0.5*(newhit.ypos + (*p2).ypos);
+      newhit.zpos = 0.5*(newhit.zpos + (*p2).zpos);
+      newhit.deposit = newhit.deposit + (*p2).deposit;
+      //TODO: as convention the one with largest energy dominates, should it be deposit instead?
+      newhit.energy = newhit.energy > (*p2).energy ? newhit.energy : (*p2).energy;
+      newhit.id = newhit.energy > (*p2).energy ? newhit.id : (*p2).id;
+      //TODO: Dark photon always take precedence
+      if(abs(newhit.particle) == 100)
+       newhit.particle = newhit.particle;
+      else if(abs((*p2).particle) == 100)
+       newhit.particle = (*p2).particle;
+      else
+       newhit.particle = newhit.energy > (*p2).energy ? newhit.particle : (*p2).particle;
+      //erase merged hit
+      p2 = mmhits.erase(p2);
+     }// merging condition
+    else
+     {
+      ++p2;
+     }
+    }
+   mmhitsmerged.push_back(newhit);
+  }
+ return mmhitsmerged;
+} // end do merging function
+
+//CREATE SAME REDUNDANCY OBSERVED IN THE DATA DUE TO THE FAKE HIT MATCHING BETWEEN PLANES
+vector<MChit> CreateFakeHits(const vector<MChit>& mmhits)
+{
+ vector<MChit> mmhitsfake = mmhits;
+ for(auto p1 = mmhits.begin(); p1 != mmhits.end(); ++p1)
+  {
+   MChit newhit;
+   for(auto p2 = mmhits.begin(); p2 != mmhits.end(); ++p2)
+    {
+     //skip true combination
+     if(p1 == p2)continue;
+     //DISTANCE IS TOO SHORT, DO MERGING
+     newhit.xpos = (*p1).xpos;
+     newhit.ypos = (*p2).ypos;
+     newhit.zpos = (*p2).zpos;
+     newhit.deposit = 0.5*((*p1).deposit + (*p2).deposit);
+     //TODO: as convention the one with largest energy dominates, should it be deposit instead?
+     newhit.energy = (*p1).energy > (*p2).energy ? (*p1).energy : (*p2).energy;
+     //put Geant ID at -1 to see that it is a fake hits
+     newhit.id = -1;
+     //TODO: Dark photon always take precedence
+     if(abs((*p1).particle) == 100)
+      newhit.particle = (*p1).particle;
+     else if(abs((*p2).particle) == 100)
+      newhit.particle = (*p2).particle;     
+     else
+      newhit.particle = (*p1).energy > (*p2).energy ? (*p1).particle : (*p2).particle;
+     mmhitsfake.push_back(newhit);
+     /* cout<< "new hit created" << endl;      */
+     /* cout << newhit; */
+    }
+  }
+ return mmhitsfake;
+} // end do merging function
+
+
+void DoMCReco (Micromega& mm,vector<MChit>& hits,const MM_calib& calib,const DetGeo& geom)
+{
+ //eliminate redundant hits
+ MCRealHit(hits,MM_threshold);
+ //assign calib and geom
+ mm.Init(calib, geom);
+ //simulate multiplexing and real condition
+ //MCDoHitAccumulation(mm,hits,1);
+ MCDoHitAccumulation(mm,hits);
+ //simulate noise
+ //if(simulatenoise) 
+  //SimulateNoise(mm);
+ //do the rest of the reconstruction as for data
+ mm.DoMMReco();
+}
+
+void DoMCRecoGEM(vector<MChit>& hits)
+{
+ //take only physical hit
+ 
+ MCRealHit(hits,GEM_threshold);
+ 
+ //smear the hits
+ DoSmearing(hits,GEMxres,GEMyres);
+ 
+ //MERGED HITS
+ hits = DoMerging(hits);
+ 
+ //create fake hits
+ if(ReplicateHits)
+  {
+   hits = CreateFakeHits(hits);
+  }
+}
+
+//DIGITALIZE MC INFORMATION
+
+void DoDigitization(RecoEvent& e)
+{
+
+ //do reconstruction in mc for Micromegas
+ DoMCReco(e.MM1,e.mc->MM1truth,MM1_calib,MM1_pos);
+ DoMCReco(e.MM2,e.mc->MM2truth,MM2_calib,MM2_pos);
+ DoMCReco(e.MM3,e.mc->MM3truth,MM3_calib,MM3_pos);
+ DoMCReco(e.MM4,e.mc->MM4truth,MM4_calib,MM4_pos);
+ DoMCReco(e.MM5,e.mc->MM5truth,MM1_calib,MM5_pos);
+ DoMCReco(e.MM6,e.mc->MM6truth,MM2_calib,MM6_pos);
+ DoMCReco(e.MM7,e.mc->MM7truth,MM3_calib,MM7_pos);
+ DoMCReco(e.MM8,e.mc->MM8truth,MM4_calib,MM8_pos);
+ //merge gem hits
+ DoMCRecoGEM(e.mc->GEM1truth);
+ DoMCRecoGEM(e.mc->GEM2truth);
+ DoMCRecoGEM(e.mc->GEM3truth);
+ DoMCRecoGEM(e.mc->GEM4truth);
+
+}
+
+void LoadGeometryFromCalo()  
+{
+  // loop over MC file metadata
+  for (auto p : MCRunInfo)
+    {
+      const string key   = p.first;
+      const string value = p.second;
+      
+      if( key == "ECAL:ECALStart") ECAL_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "ECAL:ECALX") ECAL_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "General:BeamMomentum") NominalBeamEnergy = mcrunoption<double>(key, 100.);
+      else if( key == "MM:MM1X") MM1_pos.pos.SetX(mcrunoption<double>(key, 0.));
+      else if( key == "MM:MM1Y") MM1_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM1Z") MM1_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM2X") MM2_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM2Y") MM2_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM2Z") MM2_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM3X") MM3_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM3Y") MM3_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM3Z") MM3_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM4X") MM4_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM4Y") MM4_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM4Z") MM4_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM5X") MM5_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM5Y") MM5_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM5Z") MM5_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM6X") MM6_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM6Y") MM6_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM6Z") MM6_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM7X") MM7_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM7Y") MM7_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM7Z") MM7_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM8X") MM8_pos.pos.SetX(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM8Y") MM8_pos.pos.SetY(mcrunoption<double>(key, 0.));    
+      else if( key == "MM:MM8Z") MM8_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      //else if( key == "Magnet:Membrane0Z") MAGU_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      //else if( key == "Magnet:Membrane1Z") MAGD_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "Magnet:Magnet2End") MAGU_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      else if( key == "Magnet:Magnet1Start") MAGD_pos.pos.SetZ(mcrunoption<double>(key, 0.)-1000.);
+      else if( key == "Magnet:MagnetMagY") MAG_field = -1.*mcrunoption<double>(key, 0.);
+      else if( key == "SRD:ZSandwichBegin")SRD_pos.pos.SetZ(mcrunoption<double>(key, 0.));    
+      
+      // WCAL_pos is absent on conddb.h for now
+      //else if( key == "ECAL:ECAL1Start") WCAL_pos.pos.SetZ(mcrunoption<double>(key, 0.));  //ECAL1 is the name of WCAL in the simulation 
+      //else if( key == "ECAL:ECAL1X") WCAL_pos.pos.SetX(mcrunoption<double>(key, 0.)); //ECAL1 is the name of WCAL in the simulation
+    }
+}
-- 
GitLab


From 1313be5d9037910a9e68334fa2a028b15040e9f5 Mon Sep 17 00:00:00 2001
From: Dmitry Shchukin <dmitry.shchukin@cern.ch>
Date: Tue, 11 Oct 2022 14:09:03 +0200
Subject: [PATCH 2/3] clean up and helper script

---
 modifyMakefile.sh | 3 +++
 1 file changed, 3 insertions(+)
 create mode 100755 modifyMakefile.sh

diff --git a/modifyMakefile.sh b/modifyMakefile.sh
new file mode 100755
index 0000000..66c63c7
--- /dev/null
+++ b/modifyMakefile.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+#modify in ../Makefile to allow for simu.h usage
+sed -i.bak -r 's#na64\-kirsanov\/mkcode\.exe\: na64\-kirsanov\/\*\.inc#na64\-kirsanov\/mkcode\.exe\: na64\-kirsanov\/\*\.inc na64\-kirsanov\/simu\.h#g' ../Makefile
-- 
GitLab


From 5458f2a01419b18df477d048c5e4107f0e6190c7 Mon Sep 17 00:00:00 2001
From: Dmitry Shchukin <dmitry.shchukin@cern.ch>
Date: Tue, 11 Oct 2022 14:10:48 +0200
Subject: [PATCH 3/3] cleanup 2

---
 mkcode.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mkcode.cc b/mkcode.cc
index 442596a..9bb452d 100644
--- a/mkcode.cc
+++ b/mkcode.cc
@@ -551,19 +551,19 @@ int procev(RecoEvent &e0)
     } 
     else {
       if(MCFileFormatVersion >= 8) {
-        //Init_Geometry_MC2022B();
         if(ievproc == 0) {
           Reset_Calib();
           Init_Calib_2022B_MM(1664928923);
           LoadGeometryFromCalo();
         }
         DoDigitization(e);
+        
         man.Reset();
         man.AddMicromega(e.MM1);    
         man.AddMicromega(e.MM2);    
         man.AddMicromega(e.MM3);
         man.AddMicromega(e.MM4);
-        
+     
         track = simple_tracking_4MM_2016(e);
         if(track) mymom = track.momentum;
       }
-- 
GitLab