Merged auto testbeam geometry into master branch

Analysis/include/LHCf_reader.h

@@ -12,12 +12,14 @@
 
 #include "TTree.h"
 #include "TFile.h"
+#include "TKey.h"
 #include "TBranch.h"
 #include "TLeaf.h"
 #include "TCanvas.h"
 #include "TH1D.h"
 #include "TH2D.h"
 #include <vector>
+#include <string>
 
 #include "DataReader2021.h"
 #include "LHCf_Trigger1.h"
@@ -35,6 +37,8 @@ class LHCf_reader : public DataReader2021{
 
    void OpenFile( void );
    void OpenFileName( std::string _fileName );
+   void ConvertFile( void );
+   void GetEventWaveforms( int iEvent );
    void InitializeHistograms( void );
 
    void SetSkipPedestal ( bool _wantToSkip ) { m_skipPedestal = _wantToSkip; }
@@ -44,9 +48,9 @@ class LHCf_reader : public DataReader2021{
    void ShowEvent( int iEvent );
    void ReadEvents( int iStartEvent = -1, int iEndEvent = -1 );
 
-   std::vector < unsigned int > GetGroupSample(unsigned int dwrd_0,
-                                               unsigned int dwrd_1,
-                                               unsigned int dwrd_2);
+   std::vector < float > GetGroupSample(unsigned int dwrd_0,
+                                        unsigned int dwrd_1,
+                                        unsigned int dwrd_2);
 
    void VisualizeChannels ( void );
 
@@ -63,8 +67,8 @@ class LHCf_reader : public DataReader2021{
    std::vector < Channel* > m_v_Channel;
 
    //std::vector < Channel* > m_v_digitizedChannels;
-   std::vector < std::vector < unsigned int > > m_v_Channels;
-   std::vector < std::vector < unsigned int > > m_v_Channels_Corrected;
+   std::vector < std::vector < float > > m_v_Channels;
+   std::vector < std::vector < float > > m_v_Channels_Corrected;
 
    bool m_debug = true;
 
@@ -74,6 +78,7 @@ class LHCf_reader : public DataReader2021{
    bool m_skipBeam = false;
 
    unsigned int m_nGroups = 3;
+   unsigned int m_nChPerGroup= 8;
 
 };
 #endif

Analysis/include/caen_correction.h

@@ -15,7 +15,7 @@ public:
   virtual ~caen_correction () {};
 
   int init( unsigned int *index_cell, int nChips,
-            std::vector < std::vector < unsigned int > > m_v_Channels );
+            std::vector < std::vector < float > > m_v_Channels );
 
   float caen_corrected(const int sample, const int channel) const;
   float caen_time(const int sample, const int channel) const;

Analysis/src/Detector.cpp

@@ -139,7 +139,7 @@ Channel* Detector::GetElement(std::string _name){
 void Detector::SetBranches( TTree *_dataTree ){
 
     for( uint ch = 0; ch < m_Element.size(); ch++ ){
-      _dataTree->SetBranchAddress( ("Raw" + m_Element.at(ch)->name).c_str(), &m_Element.at(ch)->pWF );
+      if( m_Element[ch]->is_on) _dataTree->SetBranchAddress( ("Raw" + m_Element.at(ch)->name).c_str(), &m_Element.at(ch)->pWF );
     }
 
 }
@@ -183,6 +183,7 @@ void Detector::DeclareHistograms(){
  */
 void Detector::FillHistograms(){
     for( uint ch = 0; ch < m_Element.size(); ch++ ){
+        if (!m_Element.at(ch)->is_on) continue;
         m_Element.at(ch)->WF_histo->Reset();
         m_Element.at(ch)->PWF_histo->Reset();
         m_Element.at(ch)->FirstDerivative->Reset();

Analysis/src/LHCf_reader.cpp

@@ -74,7 +74,7 @@ void LHCf_reader::OpenFileName( std::string _fileName ){
 }
 
 // Get 3 32 bit dwords and return 8 Ch sample values of the group Ch 0,1,2…7
-std::vector<unsigned int>  LHCf_reader::GetGroupSample(unsigned int dwrd_0, unsigned int dwrd_1, unsigned int dwrd_2){
+std::vector<float>  LHCf_reader::GetGroupSample(unsigned int dwrd_0, unsigned int dwrd_1, unsigned int dwrd_2){
   /*
    Get 3 32 bit dwords and return 8 Ch sample values
    of the group Ch 0,1,2…7.
@@ -87,7 +87,7 @@ std::vector<unsigned int>  LHCf_reader::GetGroupSample(unsigned int dwrd_0, unsi
     x/(2^n)%2^m
     ===============================================
   */
-  std::vector<unsigned int> sample;
+  std::vector<float> sample;
   sample.push_back(dwrd_0%4096);                                 // CH0
   sample.push_back((dwrd_0/4096)%4096);                          // CH1
   sample.push_back((dwrd_0/16777216)%4096 ^ (dwrd_1%16)*256);    // CH2
@@ -104,6 +104,180 @@ void LHCf_reader::ShowEvent( int iEvent ){
   m_readTrigger1->Show( iEvent );
 }
 
+void LHCf_reader::ConvertFile(){
+   std::string outputDir =  Form("%s/ProcessedLHCf",std::getenv("JZCaPA"));
+   gSystem->Exec( Form("mkdir -p %s", outputDir.c_str() ) );
+   std::string base_filename = m_fNameIn.substr(m_fNameIn.find_last_of("/\\") + 1);
+   std::string of = Form("%s/%s",outputDir.c_str() ,base_filename.c_str() );
+   if (m_debug) printf(" output file:%s\n",of.c_str());
+   TFile* out_file = TFile::Open(of.c_str(),"RECREATE"); 
+   TObject *obj;
+   TKey *key;
+   TTree* tree[5];
+   TIter next( m_filePointer->GetListOfKeys());
+   int id = 0;
+   bool readCurrent = 0;
+   int RunNumber;
+   int evt;
+
+   // Get runNumber
+  std::string fname_str(base_filename.c_str());
+  std::string result = fname_str.substr(0, fname_str.find_last_of("."));
+  RunNumber =  atoi(result.substr(result.find_last_not_of("0123456789") + 1).c_str());
+
+    for ( int iCh = 0; iCh < m_nCh; iCh++){
+      std::vector < float > v_buffer, v_buffer2;
+      v_buffer.resize( m_nSamp );
+      v_buffer2.resize( m_nSamp );
+      m_v_Channels.push_back( v_buffer );
+      m_v_Channels_Corrected.push_back( v_buffer2 );
+    }
+
+
+   while ((key = (TKey *) next())) {
+      obj = m_filePointer->Get(key->GetName()); // copy object to memory
+      if (strcmp(key->GetClassName(),"TTree")) continue;
+      if (m_debug) printf(" found object:%s\n",key->GetName());
+      if (!strcmp(((TTree*)obj)->GetName(),"Trigger1")){
+        if (readCurrent) continue;
+        readCurrent = 1;
+        // Do Things with Trigger1
+        TTree* tree_to_read = (TTree*)obj;
+        if (m_debug) printf(" Nubmer of entries:%lld\n",tree_to_read->GetEntries());
+        // Deactivate the branch which needs to be processed
+        tree_to_read->SetBranchStatus("ZDCH",0);
+        // Clone skimmed
+        if (m_debug) printf(" About to clone LHCf branches\n");
+        tree[id] = tree_to_read->CloneTree(0); 
+        tree[id]->CopyEntries(tree_to_read);
+        tree[id]->SetBranchStatus("*",0);
+        if (m_debug) printf(" Successfully cloned LHCf branches\n"); 
+        // Reactivate ZDCH
+        tree_to_read->SetBranchStatus("*",0);
+
+        tree_to_read->SetBranchStatus("ZDCH",1);
+        //Create run number and event number branches
+
+        tree[id]->Branch("RunNumber", &RunNumber, "RunNumber/I" );
+        tree[id]->Branch("EventNumber", &evt, "EventNumber/I" );
+        // Generate new Branches for unpacking
+        if (m_debug) printf(" About to generate ZDCH WF branches\n");
+        for (int ch=0;ch<m_nChPerGroup*m_nGroups;ch++){
+          tree[id]->Branch( Form( "RawSignal%d", ch), "std::vector<float>", &m_v_Channels_Corrected.at(ch) );
+        }
+        if (m_debug) printf(" Finished generating ZDCH WF branches\n");
+        // Run over all events and unpack, as per the 2 event shift
+        int MaxEvents = tree_to_read->GetEntries();
+        for(evt=0;evt<MaxEvents;evt++){
+          // Get Samples from evt+2
+          if ((MaxEvents-evt) < 3) continue;// Fill Empty
+          GetEventWaveforms(evt+2);
+           
+          tree[id]->Fill();
+        }
+        tree[id]->SetEntries(MaxEvents);
+        tree[id]->SetName("tree"); // Temporary to see if solves JZCaPA issue
+      }else{
+        if (m_debug) printf(" Writing to TTree with ID: %d\n",id);
+        tree[id] = ((TTree*)obj)->CloneTree();
+      }
+      id++;
+  
+   }
+   if (m_filePointer) m_filePointer->Close();
+   
+   if (m_debug) printf(" Moving into new file.\n");
+   out_file->cd();
+   if (m_debug) printf(" Writing into new file.\n");
+   out_file->Write();
+   if (m_debug) printf(" Finished writing.\n");
+   out_file->Close();
+}
+
+
+void LHCf_reader::GetEventWaveforms(int iEvent){
+    m_v_Channels.clear();
+    m_v_Channels_Corrected.clear();
+
+    m_readTrigger1->GetEntry( iEvent );
+    bool trg_beam = false;
+    bool trg_pede = false;
+    trg_beam = (m_readTrigger1->GPI0_GPI0[0] >> 8) & 0x1;
+    trg_pede = (m_readTrigger1->GPI0_GPI0[0] >> 5) & 0x1;
+
+    //Checking trigger ADCs
+    double adc[2];
+    adc[0] = m_readTrigger1->ADC2_ADC2[0] - m_readTrigger1->ADC2_ADC2[32];// trg.scin 1
+    adc[1] = m_readTrigger1->ADC2_ADC2[2] - m_readTrigger1->ADC2_ADC2[34];// trg.scin 2
+
+
+    if(trg_pede) {
+      std::cout << "PEDESTAL TRIGGER - SKIP " << std::endl;
+     // return;
+    }
+    if(trg_beam){
+       std::cout << "BEAM TRIGGER " << std::endl;
+     }
+
+    for ( int iCh = 0; iCh < m_nCh; iCh++){
+      std::vector < float > v_buffer, v_buffer2;
+      v_buffer.resize( m_nSamp );
+      v_buffer2.resize( m_nSamp );
+      m_v_Channels.push_back( v_buffer );
+      m_v_Channels_Corrected.push_back( v_buffer2 );
+    }
+
+    //Retrieving calibration file from Util folder
+    std::string calibFile = std::getenv("JZCaPA") + std::string("/Utils/calib_0234_2G_0.dat");
+
+    caen_correction *cc;
+    cc = new caen_correction ( calibFile.c_str() );
+
+    unsigned int curser = 4;
+    unsigned int n_words_for_group_data_extraction;
+    unsigned int c_stop_flag;
+
+    //NSI - Noise subtraction implementation
+    unsigned int index_cell[m_nGroups];
+    // For group 1
+    for (int gid = 0; gid < m_nGroups; gid++){
+      n_words_for_group_data_extraction = m_readTrigger1->ZDCH_ZDCH[curser]%4096;
+      index_cell[gid] = (m_readTrigger1->ZDCH_ZDCH[curser]/1048576)%4096;
+      if(m_debug)
+        std::cout << "GID: " << gid << " n_words_for_group_data_extraction: " << n_words_for_group_data_extraction << std::endl;
+      //std::cout << "TR: " << (m_readTrigger1->ZDCH_ZDCH[curser]/4096)%2 << std::endl;
+      curser += 1;
+      c_stop_flag = curser + n_words_for_group_data_extraction;
+    // run over the dwords and extract each sample
+      int iBin = 0;
+      if(m_debug)
+        std::cout << "PRE --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
+
+      for (int dummy = 0; curser < c_stop_flag; curser+=3){
+        std::vector < float > sample_vec = GetGroupSample(m_readTrigger1->ZDCH_ZDCH[curser],m_readTrigger1->ZDCH_ZDCH[curser+1],m_readTrigger1->ZDCH_ZDCH[curser+2]);
+        // store samples
+        for ( int iSamp = 0; iSamp < m_nChPerGroup; iSamp++){
+          m_v_Channels.at( (gid * m_nChPerGroup) + iSamp ).at(iBin) = sample_vec.at(iSamp);
+        }
+        iBin++;
+      }
+      if(m_debug)
+        std::cout << "POST --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
+     // Move to next group data words
+     curser += 1;
+    }
+
+    // Perform signal correction using firmware values
+    cc->init(index_cell, m_nGroups, m_v_Channels);
+    for(int iCh = 0; iCh < m_nChPerGroup*m_nGroups; iCh++){
+      for(int iSamp = 0; iSamp < m_nSamp; iSamp++){
+        m_v_Channels_Corrected.at( iCh ).at( iSamp ) = cc->caen_corrected( iSamp, iCh );
+      }
+    }
+
+}
+
+
 void LHCf_reader::ReadEvents( int iStartEvent, int iEndEvent ){
 
   if( iStartEvent == -1 ) iStartEvent = 0;
@@ -138,7 +312,7 @@ void LHCf_reader::ReadEvents( int iStartEvent, int iEndEvent ){
      }
 
     for ( int iCh = 0; iCh < m_nCh; iCh++){
-      std::vector < unsigned int > v_buffer, v_buffer2;
+      std::vector < float > v_buffer, v_buffer2;
       v_buffer.resize( m_nSamp );
       v_buffer2.resize( m_nSamp );
       m_v_Channels.push_back( v_buffer );
@@ -173,7 +347,7 @@ void LHCf_reader::ReadEvents( int iStartEvent, int iEndEvent ){
         std::cout << "PRE --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
 
       for (int dummy = 0; curser < c_stop_flag; curser+=3){
-        std::vector < unsigned int > sample_vec = GetGroupSample(m_readTrigger1->ZDCH_ZDCH[curser],m_readTrigger1->ZDCH_ZDCH[curser+1],m_readTrigger1->ZDCH_ZDCH[curser+2]);
+        std::vector < float > sample_vec = GetGroupSample(m_readTrigger1->ZDCH_ZDCH[curser],m_readTrigger1->ZDCH_ZDCH[curser+1],m_readTrigger1->ZDCH_ZDCH[curser+2]);
         // store samples
         for ( int iSamp = 0; iSamp < nPerGroup; iSamp++){
           m_v_Channels.at( (gid * nPerGroup) + iSamp ).at(iBin) = sample_vec.at(iSamp);

Analysis/src/caen_correction.cpp

@@ -6,6 +6,11 @@
 
 using namespace std;
 
+/*
+Original code by Martin Purschke (sPhenix)
+Modified for JZCaPA LHCf data by Riccardo Longo, Yftach Moyal
+*/
+
 caen_correction::caen_correction ( const char *calibdata)
 {
 
@@ -71,7 +76,7 @@ caen_correction::caen_correction ( const char *calibdata)
 }
 
 int caen_correction::init ( unsigned int *index_cell, int nChips,
-                            std::vector < std::vector < unsigned int > > m_v_Channels )
+                            std::vector < std::vector < float > > m_v_Channels )
 {
   int chip,c,i,idx;
 

Analysis/userFunctions/convertLHCf.cpp

@@ -23,8 +23,9 @@ int main(int argc, char *argv[]){
   LHCf_reader* myReader = new LHCf_reader( fNameIn.c_str() );
   myReader->OpenFile();
 
-  myReader->ReadEvents( 215, 220  );
+  //myReader->ReadEvents( 215, 220  );
+  myReader->ConvertFile();
 
-  myReader->VisualizeChannels();
+  //myReader->VisualizeChannels();
   return 0;
 }

MonteCarlo/include/DetectorConstruction.hh

@@ -109,6 +109,12 @@ class Alignment {
     double x_table;
     /** Y position of the Desy Table **/
     double y_table;
+    /** Energy of the beam **/
+    double beam_energy;
+    /** Scan name of the beam **/
+    std::string beam_scan;
+    /** Type of the beam **/
+    std::string beam_type;
     /** First detector met by the beam **/
     std::string upstream_Det;
     /** Second detector met by the beam **/
@@ -131,6 +137,8 @@ class Alignment2021 {
 
     /** Run number being analyzed **/
     int runNumber;
+    /** Scan number **/
+    int scanNumber;
     /** Beam type - can be p or e in H2 2021 **/
     std::string beam_type;
     /** Beam energy **/
@@ -190,11 +198,12 @@ public:
   inline  void    SetPI0Flag            ( G4bool arg ){PI0 = arg;}
 
   virtual G4RotationMatrix QuaternionToRotationMatrix( TQuaternion Q );
-  virtual Survey*    GetSurvey             ( G4String name );
-  inline  Alignment  GetAlignment          ( ){ return *m_alignment;    }
-  inline  G4bool     GetOverlapsFlag       ( ){ return CHECK_OVERLAPS;  }
-  inline  G4bool     GetOpticalFlag        ( ){ return OPTICAL;         }
-  inline  G4bool     GetPI0Flag            ( ){ return PI0;             }
+  virtual Survey*       GetSurvey         ( G4String name );
+  inline  Alignment     GetAlignment      ( ){ return *m_alignment;    }
+  inline  Alignment2021 GetAlignment2021  ( ){ return *m_alignment2021;}
+  inline  G4bool        GetOverlapsFlag   ( ){ return CHECK_OVERLAPS;  }
+  inline  G4bool        GetOpticalFlag    ( ){ return OPTICAL;         }
+  inline  G4bool        GetPI0Flag        ( ){ return PI0;             }
   inline  std::vector< ModTypeZDC* >* GetZDCvec ( ){ return &m_ZDCvec; }
   inline  std::vector< ModTypeRPD* >* GetRPDvec ( ){ return &m_RPDvec; }
   inline  ModTypeTRI* GetTRI ( ){ return m_TRI; }

MonteCarlo/include/TestBeam2021PrimaryGenerator.hh

@@ -36,6 +36,7 @@
 
 #include "G4VPrimaryGenerator.hh"
 #include "DetectorConstruction.hh"
+#include "Randomize.hh"
 
 class G4Event;
 
@@ -50,7 +51,13 @@ class TestBeam2021PrimaryGenerator : public G4VPrimaryGenerator
   public:
     virtual void GeneratePrimaryVertex(G4Event*);
 
-    Alignment m_alignment;
+  private:
+    void GetWireChamberDist();
+
+    Alignment2021 m_alignment;
+    G4RandGeneral *m_RandX;
+    G4RandGeneral *m_RandY;
+    int m_WireChamberSize;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MonteCarlo/run1.mac

@@ -78,7 +78,7 @@
 # Descriptions of commands can be found in the UI command tree
 #################################
 
-/beam/type lhc
+/beam/type TB2021
 /beam/pos 0. 0. -140.
 /beam/projectBeam -80. mm
 /beam/nPrimaries 1
@@ -117,6 +117,6 @@
 ##############################################################
 ##############################################################
 # number of events
-/run/beamOn 1
+/run/beamOn 4
 ##############################################################
 ##############################################################

MonteCarlo/src/DetectorConstruction.cc

@@ -690,7 +690,7 @@ void DetectorConstruction::ConstructSDandField( ){
 
 void DetectorConstruction::LoadConfigurationFile( G4String _inFile  ){
 
-  if(_inFile = ""){
+  if(_inFile == ""){
     _inFile = std::getenv("JZCaPA");
     _inFile += "/Utils/Survey_2018.xml";
   }
@@ -762,6 +762,7 @@ void DetectorConstruction::LoadAlignmentFile2021( G4String _inFile ){
         m_XMLparser->getChildValue("Alignment",i,"run",run);
         if(run != m_runNumber) continue;
         std::cout << "Found Run Entry in Alignment file for run " << m_runNumber << std::endl;
+        m_XMLparser->getChildValue("Alignment",i,"beam_scan",m_alignment2021->scanNumber);
         m_XMLparser->getChildValue("Alignment",i,"beam_type",m_alignment2021->beam_type);
         m_XMLparser->getChildValue("Alignment",i,"beam_energy",m_alignment2021->beam_energy);
         m_XMLparser->getChildValue("Alignment",i,"x_table",m_alignment2021->x_det_table);
@@ -789,7 +790,7 @@ void DetectorConstruction::LoadAlignmentFile2021( G4String _inFile ){
 void DetectorConstruction::LoadAlignmentFile( G4String _inFile ){
   bool debug = false;
 
-  if( _inFile = ""){
+  if( _inFile == ""){
     _inFile = std::getenv("JZCaPA");
     _inFile += "/Utils/Alignment_2018.xml";
   }
@@ -819,6 +820,9 @@ void DetectorConstruction::LoadAlignmentFile( G4String _inFile ){
         if(debug){
           std::cout << "Found Run Entry in Alignment file for run " << m_runNumber << std::endl;
         }
+        m_XMLparser->getChildValue("Alignment",i,"beam_scan",m_alignment->beam_scan);
+        m_XMLparser->getChildValue("Alignment",i,"beam_type",m_alignment->beam_type);
+        m_XMLparser->getChildValue("Alignment",i,"beam_energy",m_alignment->beam_energy);
         m_XMLparser->getChildValue("Alignment",i,"x_table",m_alignment->x_table);
         m_XMLparser->getChildValue("Alignment",i,"y_table",m_alignment->y_table);
         m_XMLparser->getChildValue("Alignment",i,"upstream_Det",m_alignment->upstream_Det);

MonteCarlo/src/TestBeam2021PrimaryGenerator.cc

@@ -41,13 +41,13 @@
 #include "G4PrimaryVertex.hh"
 #include "G4PhysicalConstants.hh"
 #include "G4SystemOfUnits.hh"
-#include "Randomize.hh"
 #include "G4String.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 TestBeam2021PrimaryGenerator::TestBeam2021PrimaryGenerator()
-: G4VPrimaryGenerator()
+: G4VPrimaryGenerator(),
+  m_WireChamberSize(101)
 { }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -59,41 +59,61 @@ TestBeam2021PrimaryGenerator::~TestBeam2021PrimaryGenerator()
 
 void TestBeam2021PrimaryGenerator::GeneratePrimaryVertex(G4Event* event)
 {
-
   //Get the alignment from DetectorConstruction
-  DetectorConstruction* dC
-  = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction();
   if(event->GetEventID() == 0){
-    G4String _inFile = (std::string)std::getenv("JZCaPA") + "/Utils/Survey_2021.xml";
-    dC->LoadAlignmentFile( _inFile.c_str() );
-    m_alignment = dC->GetAlignment();
+    GetWireChamberDist();
   }
 
-  G4cout << "This doesn't do anything yet" << G4endl;
-
-  // // Define the position of the particle
-  // G4double x = yourRandomDistributionX; //x position
-  // G4double y = yourRandomDistributionY; //y position
-  // G4ThreeVector position(x,y,0); //Insert x and y into a vector. Beam will be at zero
-  // G4double time = 0*s; //Time is always zero for us
-  // //
-  // G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time); //Insert the position and time into a primary vertex
-  //
-  //
-  // // A note for Hongbo:
-  // // I think for now we can use the particle type as defined by the .xml file
-  // // but for future work it might be good for us to simulate an impure beam.
-  // // What I've written below should work fine. You just need to worry about x and y above
-  //
-  // G4ParticleDefinition* particleDefinition
-  //          = G4ParticleTable::GetParticleTable()->FindParticle(alignment.beam_type.c_str()); //get the particle definition for this beam
-  // G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition); //create a particle with the given definition
-  // particle->SetMomentumDirection(G4ThreeVector(0,0,1)); //set the momentum to be purely in Z
-  // particle->SetKineticEnergy(alignment.beam_energy); //get the beam energy from the alignment file
-  // vertex->SetPrimary(particle); //add the particle to the vertex created above
-  //
-  // event->AddPrimaryVertex(vertex); //add the vertex to this event for geant to simulate
+  G4double x = m_RandX->shoot()*(m_WireChamberSize-1)-(m_WireChamberSize-1)/2; //Map (0, 1) to (-50, 50), RandNum * 100 - 50
+  G4double y = m_RandY->shoot()*(m_WireChamberSize-1)-(m_WireChamberSize-1)/2;
+  G4ThreeVector position(x, y, 0);
+  G4double time = 0*s;
+  G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time);
+
+  G4cout << "Beam scan & type & energy ==> (" << m_alignment.scanNumber << m_alignment.beam_type << ", " << m_alignment.beam_energy << ")" << G4endl;
+  G4cout << "Particle Position ==> (" << position.x() << ", " << position.y() << ", " << position.z() << ")" << G4endl;
+
+  G4ParticleDefinition* particleDefinition
+           = G4ParticleTable::GetParticleTable()->FindParticle(m_alignment.beam_type.c_str()); // about the definition of the ParticleTable, I changed the Alignment_2021.xml file, "e" to "e-", "p" to "proton"
+  G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition);
+  particle->SetMomentumDirection(G4ThreeVector(0,0,1));
+  particle->SetKineticEnergy(m_alignment.beam_energy*GeV); //I'm not sure if I need units here
+  vertex->SetPrimary(particle);
+
+  event->AddPrimaryVertex(vertex);
 
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void TestBeam2021PrimaryGenerator::GetWireChamberDist(){
+  DetectorConstruction* dC
+  = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction();
+  m_alignment = dC->GetAlignment2021();
+
+  // Define the position of the particle
+  G4String _profile = (std::string)std::getenv("JZCaPA") + "/Utils/2021testbeam_scan/scan" + std::to_string(m_alignment.scanNumber) + ".txt"; //Select the corresponding scan file according to run number
+
+  // Attempt to open the file and check if it succeeded
+  std::ifstream _scanfile(_profile);
+  if ( !_scanfile.is_open() ){
+    G4cout << " ==== This Run does not belong to any processed Scan ==== " << G4endl; // You can check ("JZCaPA")/Utils/2021testbeam_scan/README.md
+    return;
+  }
+
+  G4double probListX[m_WireChamberSize],
+           probListY[m_WireChamberSize];
+  // Put the X and Y beam profile into two arrays
+  for(int j = 0; j < 2; ++j){
+    for(int i = 0; i <m_WireChamberSize; ++i){
+      if (!j){
+        _scanfile >> probListX[i];
+      }
+      _scanfile >> probListY[i];
+    }
+  }
+  _scanfile.close();
+
+  m_RandX = new G4RandGeneral(probListX,m_WireChamberSize); //Generate the random number, range (0, 1)
+  m_RandY = new G4RandGeneral(probListY,m_WireChamberSize);
+}

Utils/2021testbeam_scan/README.md

+# Each processed scan data of 2021 test beam has two column,
+# First column is the profile of X, the second is Y
+
+
+# The relation between Scans and Runs are as follows:
+scan1: run12-15
+scan2: run54-98
+scan3: run117-136
+scan4: run139-157
+scan5: run161-178
+scan6: run180-197
+scan7: run203-220
+scan8: run268-326
+scan9: run337-364
+scan10: run490-600
+scan11: run601-638
+scan12: run639-682
+scan13: run689-761
+(Scan up to 757, 758-761 additional statistics in the center (5,5))
\ No newline at end of file

Utils/2021testbeam_scan/scan1.txt

+    0     0
+    0     0
+    0     0
+    0     0
+    1     1
+    1     1
+    1     1
+    1     1
+    2     1
+    3     1
+    4     2
+    6     2
+    7     3
+    8     3
+   10     4
+   11     4
+   11     6
+   12     8
+   12    10