move SimPrimaryIonizationCol createAndPut to TrackHeedSimTool

Examples/options/tut_detsim_SDT_Heed.py

@@ -178,8 +178,8 @@
     #dedx_simtool.IonMobility_file ="/junofs/users/wxfang/MyGit/tmp/check_G4FastSim_20210121/CEPCSW/Digitisers/DigiGarfield/IonMobility_He+_He.txt" 
     dedx_simtool.gas_file         ="he_90_isobutane_10.gas"
     dedx_simtool.IonMobility_file ="IonMobility_He+_He.txt"
-    dedx_simtool.save_mc = True##IF this is False then ... 
-    dedx_simtool.debug = False 
+    dedx_simtool.save_mc = True
+    dedx_simtool.debug = False
     dedx_simtool.sim_pulse = True
     #dedx_simtool.model='/junofs/users/wxfang/MyGit/tmp/fork_cepcsw_20220418/CEPCSW/Digitisers/SimCurrentONNX/src/model_90He10C4H10_18mm.onnx'
     dedx_simtool.model='model_90He10C4H10_18mm.onnx'

Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp

@@ -87,8 +87,6 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) {
 
     // reset
     m_track2primary.clear();
-
-    auto SimPIonCol =  m_SimPrimaryIonizationCol.createAndPut();
 
 }
 
@@ -97,8 +95,6 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
 
     msg() << "mcCol size (after simulation) : " << mcCol->size() << endmsg;
     // save all data
-    auto SimPrimaryIonizationCol =  m_SimPrimaryIonizationCol.get();
-    //msg() << "SimPrimaryIonizationCol size ="<<SimPrimaryIonizationCol->size()<<endmsg;
     // create collections.
     auto trackercols = m_trackerCol.createAndPut();
     auto calorimetercols = m_calorimeterCol.createAndPut();

Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h

@@ -15,7 +15,6 @@
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
 #include "edm4hep/CaloHitContributionCollection.h"
-#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
 
 class Edm4hepWriterAnaElemTool: public extends<AlgTool, IAnaElemTool> {
 
@@ -130,8 +129,6 @@ class Edm4hepWriterAnaElemTool: public extends<AlgTool, IAnaElemTool> {
             "DriftChamberHitsCollection", 
             Gaudi::DataHandle::Writer, this};
 
-    // for ionized electron
-    DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>    m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this};
 
 private:
     // in order to associate the hit contribution with the primary track,

Simulation/DetSimDedx/src/TrackHeedSimTool.cpp

@@ -25,7 +25,10 @@ double TrackHeedSimTool::dndx(double betagamma) {return 0;}
 
 double TrackHeedSimTool::dedx(const G4Step* Step)
 {
-
+    if(m_beginEvt){
+        m_SimPrimaryIonizationCol =  m_SimPrimaryIonizationColWriter.createAndPut();
+        m_beginEvt = false;
+    }
     clock_t t0 = clock();
     double de = 0;
     float cm_to_mm = 10; 
@@ -43,15 +46,12 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
     if(g4Track->GetKineticEnergy() <=0) return 0;
     if(pdg_code == 11 && (tmp_str_pro=="phot" || tmp_str_pro=="hIoni" || tmp_str_pro=="eIoni" || tmp_str_pro=="muIoni" || tmp_str_pro=="ionIoni" ) ) return m_eps;//skip the electron produced by Ioni, because it is already simulated by TrackHeed
     if(m_particle_map.find(pdg_code) == m_particle_map.end() ) return m_eps;
-    edm4hep::SimPrimaryIonizationClusterCollection* SimPrimaryIonizationCol = nullptr;
     edm4hep::MCParticleCollection* mcCol = nullptr;
     try{
-        SimPrimaryIonizationCol =  const_cast<edm4hep::SimPrimaryIonizationClusterCollection*>(m_SimPrimaryIonizationCol.get());
         mcCol = const_cast<edm4hep::MCParticleCollection*>(m_mc_handle.get());
     }
     catch(...){
-        G4cout<<"Error! Can't find collection in event, please check it have been createAndPut() in Begin of event"<<G4endl;
-        G4cout<<"SimPrimaryIonizationCol="<<SimPrimaryIonizationCol<<",mcCol="<<mcCol<<G4endl;
+        G4cout<<"Error! Can't find MCParticle collection in event, please check it have been createAndPut() in Begin of event"<<G4endl;
         throw "stop here!";
     }
 
@@ -194,7 +194,8 @@ double TrackHeedSimTool::dedx(const G4Step* Step)
     clock_t t02 = clock();
     while (m_track->GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
         //auto chit = SimHitCol->create();
-        auto chit = SimPrimaryIonizationCol->create();
+        //auto chit = SimPrimaryIonizationCol->create();
+        auto chit = m_SimPrimaryIonizationCol->create();
         chit.setTime(tc);
         double cpos[3] = { cm_to_mm*( (xc - init_x)+position_x/CLHEP::cm) , cm_to_mm*((yc - init_y)+position_y/CLHEP::cm), cm_to_mm*(zc + position_z/CLHEP::cm)};
         chit.setPosition(edm4hep::Vector3d(cpos));
@@ -395,7 +396,8 @@ StatusCode TrackHeedSimTool::initialize()
   m_current_Parent_ID = -1;
   m_change_track = false;
   m_total_range = 0;
-  m_isFirst = false;
+  m_isFirst = true;
+  m_beginEvt = true;
   m_tot_edep = 0;
   m_tot_length = 0;
   m_pa_KE =0;
@@ -555,26 +557,16 @@ float* TrackHeedSimTool::NNPred(std::vector<float>& inputs)
 
 void TrackHeedSimTool::endOfEvent() {
     if(m_sim_pulse){
-        edm4hep::SimPrimaryIonizationClusterCollection* SimPrimaryIonizationCol = nullptr;
-        try{
-            SimPrimaryIonizationCol =  const_cast<edm4hep::SimPrimaryIonizationClusterCollection*>(m_SimPrimaryIonizationCol.get());
-        }
-        catch(...){
-            G4cout<<"Error! Can't find collection in event, please check it have been createAndPut() in Begin of event"<<G4endl;
-            G4cout<<"SimPrimaryIonizationCol="<<SimPrimaryIonizationCol<<G4endl;
-            throw "stop here!";
-        }
-        if(m_debug) G4cout<<"SimPrimaryIonizationCol size="<<SimPrimaryIonizationCol->size()<<G4endl;
+        if(m_debug) G4cout<<"SimPrimaryIonizationCol size="<<m_SimPrimaryIonizationCol->size()<<G4endl;
         clock_t t01 = clock();
         std::vector<float> inputs;
         std::vector<unsigned long> indexs_c;
         std::vector<unsigned long> indexs_i;
         std::map<unsigned long long, std::vector<std::pair<float, float> > > id_pulse_map;
-        for (unsigned long z=0; z<SimPrimaryIonizationCol->size(); z++) {
-            for (unsigned long k=0; k<SimPrimaryIonizationCol->at(z).electronCellID_size(); k++) {
-                //auto simIon = SimIonizationCol->at(k);
-                auto position = SimPrimaryIonizationCol->at(z).getElectronPosition(k);//mm
-                auto cellId = SimPrimaryIonizationCol->at(z).getElectronCellID(k);
+        for (unsigned long z=0; z<m_SimPrimaryIonizationCol->size(); z++) {
+            for (unsigned long k=0; k<m_SimPrimaryIonizationCol->at(z).electronCellID_size(); k++) {
+                auto position = m_SimPrimaryIonizationCol->at(z).getElectronPosition(k);//mm
+                auto cellId = m_SimPrimaryIonizationCol->at(z).getElectronCellID(k);
                 TVector3 Wstart(0,0,0);
                 TVector3 Wend  (0,0,0);
                 m_segmentation->cellposition(cellId, Wstart, Wend);
@@ -590,9 +582,9 @@ void TrackHeedSimTool::endOfEvent() {
                 float local_y = 0;
                 getLocal(wire_x, wire_y, position[0], position[1], local_x, local_y);
                 //std::cout<<"pos_z="<<pos_z<<",wire_x="<<wire_x<<",wire_y="<<wire_y<<",position[0]="<<position[0]<<",position[1]="<<position[1]<<",local_x="<<local_x<<",local_y="<<local_y<<",dr="<<sqrt(local_x*local_x+local_y*local_y)<<",dr1="<<sqrt( (wire_x-position[0])*(wire_x-position[0])+(wire_y-position[1])*(wire_y-position[1]) )<<std::endl;
-                float m_x_scale = 1;
-                float m_y_scale = 1;
-                local_x = local_x/m_x_scale;//FIXME, default is 18mm x 18mm, the real cell size maybe a bit different. need konw size of cell for each layer, and do normalization
+                //float ext_x_scale = 1;//TODO, default is 18mm x 18mm, the real cell size maybe a bit different. need konw size of cell for each layer, and do normalization using ext_x_scale
+                //float ext_y_scale = 1;
+                local_x = local_x/m_x_scale;
                 local_y = local_y/m_y_scale;
                 float noise = CLHEP::RandGauss::shoot(0,1);
                 inputs.push_back(local_x);
@@ -611,7 +603,7 @@ void TrackHeedSimTool::endOfEvent() {
                         //tmp_pluse.setValue(tmp_amp);
                         //tmp_pluse.setType(SimIonizationCol->at(tmp_index).getType());
                         //tmp_pluse.setSimIonization(SimIonizationCol->at(tmp_index));
-                        auto ion_time = SimPrimaryIonizationCol->at(indexs_c.at(i)).getElectronTime(indexs_i.at(i));
+                        auto ion_time = m_SimPrimaryIonizationCol->at(indexs_c.at(i)).getElectronTime(indexs_i.at(i));
                         id_pulse_map[indexs_c.at(i)].push_back(std::make_pair(tmp_time+ion_time,tmp_amp) );
                     }
                     inputs.clear();
@@ -632,7 +624,7 @@ void TrackHeedSimTool::endOfEvent() {
                 //tmp_pluse.setType(SimIonizationCol->at(tmp_index).getType());
                 //tmp_pluse.setSimIonization(SimIonizationCol->at(tmp_index));
                 //id_pulse_map[SimIonizationCol->at(tmp_index).getCellID()].push_back(std::make_pair(tmp_pluse.getTime(), tmp_pluse.getValue() ) );
-                auto ion_time = SimPrimaryIonizationCol->at(indexs_c.at(i)).getElectronTime(indexs_i.at(i));
+                auto ion_time = m_SimPrimaryIonizationCol->at(indexs_c.at(i)).getElectronTime(indexs_i.at(i));
                 id_pulse_map[indexs_c.at(i)].push_back(std::make_pair(tmp_time+ion_time,tmp_amp) );
             }
             inputs.clear();
@@ -641,21 +633,22 @@ void TrackHeedSimTool::endOfEvent() {
             delete [] res;
         }
         for(auto iter = id_pulse_map.begin(); iter != id_pulse_map.end(); iter++){
-            edm4hep::MutableSimPrimaryIonizationCluster dcIonCls = SimPrimaryIonizationCol->at(iter->first); 
+            edm4hep::MutableSimPrimaryIonizationCluster dcIonCls = m_SimPrimaryIonizationCol->at(iter->first); 
             for(unsigned int i=0; i< iter->second.size(); i++){
                 auto tmp_time = iter->second.at(i).first ;
                 auto tmp_amp  = iter->second.at(i).second;
                 dcIonCls.addToPulseTime(tmp_time);  
                 dcIonCls.addToPulseAmplitude(tmp_amp);  
             }
             if(dcIonCls.electronPosition_size() != dcIonCls.pulseTime_size()){
-                G4cout<<"Error ion size != pulse size"<<G4endl;
+                G4cout<<"Error: ion size != pulse size"<<G4endl;
                 throw "stop here!";
             }
         }
         clock_t t02 = clock();
         if(m_debug) std::cout<<"time for Pulse Simulation=" << (double)(t02 - t01) / CLOCKS_PER_SEC <<" seconds"<< std::endl;
     }
+    reset();
 }
 
 StatusCode TrackHeedSimTool::finalize()

Simulation/DetSimDedx/src/TrackHeedSimTool.h

@@ -53,11 +53,11 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
         double dndx(double betagamma) override;
         void getMom(float ee, float dx, float dy,float dz, float mom[3] );
         void reset(){ 
+            m_beginEvt = true;
             m_isFirst = true;
             m_previous_track_ID = 0;
             m_previous_KE = 0;
             m_tot_edep = 0;
-            //std::cout<<"m_tot_length="<<m_tot_length<<std::endl;
             m_tot_length = 0;
         }
         void endOfEvent();
@@ -89,7 +89,8 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
         Gaudi::Property<float> m_BField   {this, "BField", -3};
         Gaudi::Property<float> m_eps     { this, "eps"   , 1e-6  };//very small value, it is returned dedx for unsimulated step (may needed for SimTrackerHit)
         // Output collections
-        DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>    m_SimPrimaryIonizationCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this};
+        DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>    m_SimPrimaryIonizationColWriter{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Writer, this};
+        edm4hep::SimPrimaryIonizationClusterCollection* m_SimPrimaryIonizationCol;
         // In order to associate MCParticle with contribution, we need to access MC Particle.
         DataHandle<edm4hep::MCParticleCollection> m_mc_handle{"MCParticle", Gaudi::DataHandle::Writer, this};
 
@@ -102,19 +103,20 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
         Sensor* m_sensor;
         std::map<int, std::string> m_particle_map;
 
-        int m_previous_track_ID=0;
-        float m_previous_KE=0;
+        int m_previous_track_ID;
+        float m_previous_KE;
         int m_current_track_ID;
         int m_current_Parent_ID;
         int m_pdg_code;
         G4StepPoint* m_pre_point;
         G4StepPoint* m_post_point;
         G4double m_total_range;
-        bool m_isFirst=true;
+        bool m_isFirst;
+        bool m_beginEvt;
         bool m_change_track;
         edm4hep::MCParticle m_mc_paricle; 
-        float m_tot_edep=0;
-        float m_tot_length=0;
+        float m_tot_edep;
+        float m_tot_length;
         float m_pa_KE;
 
         G4double m_pre_x  ;
@@ -147,12 +149,12 @@ class TrackHeedSimTool: public extends<AlgTool, IDedxSimTool> {
         Gaudi::Property<bool> m_sim_pulse    { this, "sim_pulse"   , true  };
         Gaudi::Property<std::string> m_model_file{ this, "model", "model_test.onnx"};
         Gaudi::Property<int> m_batchsize     { this, "batchsize", 100};
-        Gaudi::Property<float> m_time_scale  { this, "time_scale", 99.0};
-        Gaudi::Property<float> m_time_shift  { this, "time_shift", 166.4};
-        Gaudi::Property<float> m_amp_scale   { this, "amp_scale" , 1e-2 };
-        Gaudi::Property<float> m_amp_shift   { this, "amp_shift" , 0    };
-        Gaudi::Property<float> m_x_scale     { this, "x_scale"   , 5.  };// in mm
-        Gaudi::Property<float> m_y_scale     { this, "y_scale"   , 5.  };// in mm
+        Gaudi::Property<float> m_time_scale  { this, "time_scale", 503.0};
+        Gaudi::Property<float> m_time_shift  { this, "time_shift", 814.0};
+        Gaudi::Property<float> m_amp_scale   { this, "amp_scale" , 1.15 };
+        Gaudi::Property<float> m_amp_shift   { this, "amp_shift" , 0.86 };
+        Gaudi::Property<float> m_x_scale     { this, "x_scale"   , 9.  };// in mm
+        Gaudi::Property<float> m_y_scale     { this, "y_scale"   , 9.  };// in mm
 
 
 

Simulation/DetSimInterface/include/DetSimInterface/IDedxSimTool.h

@@ -28,7 +28,6 @@ class IDedxSimTool: virtual public IAlgTool {
     virtual double dedx(const G4Step* aStep) = 0;
     virtual double dedx(const edm4hep::MCParticle& mc) = 0;
     virtual double dndx(double betagamma) = 0;
-    virtual void reset() {}
     virtual void endOfEvent() {}
 
 };

Simulation/DetSimSD/src/DriftChamberSensitiveDetector.cpp

@@ -72,5 +72,4 @@ DriftChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
 void
 DriftChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) {
     m_dedx_simtool->endOfEvent();
-    m_dedx_simtool->reset();
 }