AliceO2Group · sawenzel · Feb 19, 2019 · Feb 18, 2019 · Feb 19, 2019 · Feb 19, 2019
@@ -6,10 +6,12 @@ set(SRCS
     src/SimConfig.cxx
     src/ConfigurableParam.cxx
     src/ConfigurableParamHelper.cxx
+    src/SimCutParams.cxx
    )
 
 set(HEADERS
     include/${MODULE_NAME}/SimConfig.h
+    include/${MODULE_NAME}/SimCutParams.h
     include/${MODULE_NAME}/ConfigurableParam.h
     include/${MODULE_NAME}/ConfigurableParamHelper.h
    )

@@ -0,0 +1,34 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef O2_SIMCONFIG_SIMCUTPARAM_H_
+#define O2_SIMCONFIG_SIMCUTPARAM_H_
+
+#include "SimConfig/ConfigurableParam.h"
+#include "SimConfig/ConfigurableParamHelper.h"
+
+namespace o2
+{
+namespace conf
+{
+// parameters to influence/set the tracking cuts in simulation
+// (mostly used in O2MCApplication stepping)
+struct SimCutParams : public o2::conf::ConfigurableParamHelper<SimCutParams> {
+  bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication
+
+  double maxRTracking = 1E20;    // max R tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
+  double maxAbsZTracking = 1E20; // max |Z| tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
+
+  O2ParamDef(SimCutParams, "SimCutParams");
+};
+} // namespace conf
+} // namespace o2
+
+#endif /* O2_SIMCONFIG_SIMCUTPARAM_H_ */
@@ -17,4 +17,6 @@
 #pragma link C++ class o2::conf::SimConfig+;
 #pragma link C++ class o2::conf::SimConfigData+;
 
+#pragma link C++ class o2::conf::SimCutParams + ;
+#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::conf::SimCutParams > +;
 #endif
@@ -0,0 +1,19 @@
+/*
+ * SimCutParams.cxx
+ *
+ *  Created on: Feb 18, 2019
+ *      Author: sandro
+ */
+
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "SimConfig/SimCutParams.h"
+O2ParamImpl(o2::conf::SimCutParams);
@@ -18,8 +18,10 @@
 #ifndef STEER_INCLUDE_STEER_O2MCAPPLICATIONBASE_H_
 #define STEER_INCLUDE_STEER_O2MCAPPLICATIONBASE_H_
 
-#include "FairMCApplication.h"
+#include <FairMCApplication.h>
 #include "Rtypes.h" // for Int_t, Bool_t, Double_t, etc
+#include <TVirtualMC.h>
+#include "SimConfig/SimCutParams.h"
 
 namespace o2
 {
@@ -32,14 +34,23 @@ namespace steer
 class O2MCApplicationBase : public FairMCApplication
 {
  public:
-  using FairMCApplication::FairMCApplication;
+  O2MCApplicationBase() : FairMCApplication(), mCutParams(o2::conf::SimCutParams::Instance()) {}
+  O2MCApplicationBase(const char* name, const char* title, TObjArray* ModList, const char* MatName) : FairMCApplication(name, title, ModList, MatName), mCutParams(o2::conf::SimCutParams::Instance())
+  {
+  }
+
   ~O2MCApplicationBase() override = default;
 
+  void Stepping() override;
+
   // specific implementation of our hard geometry limits
-  double TrackingRmax() const override { return 1E20; }
-  double TrackingZmax() const override { return 1E20; }
+  double TrackingRmax() const override { return mCutParams.maxRTracking; }
+  double TrackingZmax() const override { return mCutParams.maxAbsZTracking; }
+
+ protected:
+  o2::conf::SimCutParams const& mCutParams; // reference to parameter system
 
-  ClassDefOverride(O2MCApplicationBase, 1) //Interface to MonteCarlo application
+  ClassDefOverride(O2MCApplicationBase, 1)
 };
 
 } // end namespace steer

@@ -43,6 +43,14 @@ void TypedVectorAttach(const char* name, FairMQChannel& channel, FairMQParts& pa
   }
 }
 
+void O2MCApplicationBase::Stepping()
+{
+  // dispatch first to stepping function in FairRoot
+  FairMCApplication::Stepping();
+
+  // now the right place to put our specific step/filtering criteria:
+}
+
 void O2MCApplication::initLate()
 {
   o2::utils::ShmManager::Instance().occupySegment();

@@ -0,0 +1,212 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "TFile.h"
+#include "TTree.h"
+#include "ITSMFTSimulation/Hit.h"
+#include "TOFSimulation/Detector.h"
+#include "EMCALBase/Hit.h"
+#include "TRDSimulation/Detector.h"
+#endif
+
+template <typename Hit, typename Accumulator>
+Accumulator analyse(TTree* tr, const char* brname)
+{
+  Accumulator prop;
+  auto br = tr->GetBranch(brname);
+  if (!br) {
+    return prop;
+  }
+  auto entries = br->GetEntries();
+  std::vector<Hit>* hitvector = nullptr;
+  br->SetAddress(&hitvector);
+
+  for (int i = 0; i < entries; ++i) {
+    br->GetEntry(i);
+    for (auto& hit : *hitvector) {
+      prop.addHit(hit);
+    }
+  }
+  prop.normalize();
+  return prop;
+};
+
+template <typename T>
+struct HitStats {
+  int NHits = 0;
+  double XAvg = 0.; // avg 1st moment
+  double YAvg = 0.;
+  double ZAvg = 0.;
+  double X2Avg = 0.; // avg 2nd moment
+  double Y2Avg = 0.;
+  double Z2Avg = 0.;
+  double EAvg = 0.; // average total energy
+  double E2Avg = 0.;
+  double TAvg = 0.;  // average T
+  double T2Avg = 0.; // average T^2
+
+  void print() const
+  {
+    std::cout << NHits << " "
+              << XAvg << " "
+              << YAvg << " "
+              << ZAvg << " "
+              << X2Avg << " "
+              << Y2Avg << " "
+              << Z2Avg << " "
+              << EAvg << " "
+              << E2Avg << " "
+              << TAvg << " "
+              << T2Avg << "\n";
+  }
+
+  // adds a hit to the statistics
+  void addHit(T const& hit)
+  {
+    NHits++;
+    auto x = hit.GetX();
+    XAvg += x;
+    X2Avg += x * x;
+    auto y = hit.GetY();
+    YAvg += y;
+    Y2Avg += y * y;
+    auto z = hit.GetZ();
+    ZAvg += z;
+    Z2Avg += z * z;
+    auto e = hit.GetEnergyLoss();
+    EAvg += e;
+    E2Avg += e * e;
+    auto t = hit.GetTime();
+    TAvg += t;
+    T2Avg += t * t;
+  }
+
+  void normalize()
+  {
+    XAvg /= NHits;
+    YAvg /= NHits;
+    ZAvg /= NHits;
+    X2Avg /= NHits;
+    Y2Avg /= NHits;
+    Z2Avg /= NHits;
+    EAvg /= NHits;
+    E2Avg /= NHits;
+    TAvg /= NHits;
+    T2Avg /= NHits;
+  }
+};
+
+struct ITSHitStats {
+  int NHits = 0;
+  double XAvg = 0.; // avg 1st moment
+  double YAvg = 0.;
+  double ZAvg = 0.;
+  double X2Avg = 0.; // avg 2nd moment
+  double Y2Avg = 0.;
+  double Z2Avg = 0.;
+  double EAvg = 0.; // average total energy
+  double E2Avg = 0.;
+  double TAvg = 0.;  // average T
+  double T2Avg = 0.; // average T^2
+
+  void print() const
+  {
+    std::cout << NHits << " "
+              << XAvg << " "
+              << YAvg << " "
+              << ZAvg << " "
+              << X2Avg << " "
+              << Y2Avg << " "
+              << Z2Avg << " "
+              << EAvg << " "
+              << E2Avg << " "
+              << TAvg << " "
+              << T2Avg << "\n";
+  }
+
+  // adds a hit to the statistics
+  void addHit(o2::ITSMFT::Hit const& hit)
+  {
+    NHits++;
+    auto x = hit.GetStartX();
+    XAvg += x;
+    X2Avg += x * x;
+    auto y = hit.GetStartY();
+    YAvg += y;
+    Y2Avg += y * y;
+    auto z = hit.GetStartZ();
+    ZAvg += z;
+    Z2Avg += z * z;
+    auto e = hit.GetTotalEnergy();
+    EAvg += e;
+    E2Avg += e * e;
+    auto t = hit.GetTime();
+    TAvg += t;
+    T2Avg += t * t;
+  }
+
+  void normalize()
+  {
+    XAvg /= NHits;
+    YAvg /= NHits;
+    ZAvg /= NHits;
+    X2Avg /= NHits;
+    Y2Avg /= NHits;
+    Z2Avg /= NHits;
+    EAvg /= NHits;
+    E2Avg /= NHits;
+    TAvg /= NHits;
+    T2Avg /= NHits;
+  }
+}; // end struct
+
+// do comparison for ITS
+void checkITS(TTree* reftree, TTree* testtree)
+{
+  auto refresult = analyse<o2::ITSMFT::Hit, ITSHitStats>(reftree, "ITSHit");
+  refresult.print();
+  auto testresult = analyse<o2::ITSMFT::Hit, ITSHitStats>(testtree, "ITSHit");
+  testresult.print();
+}
+
+// do comparison for TOF
+void checkTOF(TTree* reftree, TTree* testtree)
+{
+  auto refresult = analyse<o2::tof::HitType, HitStats<o2::tof::HitType>>(reftree, "TOFHit");
+  refresult.print();
+  auto testresult = analyse<o2::tof::HitType, HitStats<o2::tof::HitType>>(testtree, "TOFHit");
+  testresult.print();
+}
+
+// do comparison for EMC
+void checkEMC(TTree* reftree, TTree* testtree)
+{
+  auto refresult = analyse<o2::EMCAL::Hit, HitStats<o2::EMCAL::Hit>>(reftree, "EMCHit");
+  refresult.print();
+  auto testresult = analyse<o2::EMCAL::Hit, HitStats<o2::EMCAL::Hit>>(testtree, "EMCHit");
+  testresult.print();
+}
+
+// do comparison for TRD
+void checkTRD(TTree* reftree, TTree* testtree)
+{
+  auto refresult = analyse<o2::trd::HitType, HitStats<o2::trd::HitType>>(reftree, "TRDHit");
+  refresult.print();
+  auto testresult = analyse<o2::trd::HitType, HitStats<o2::trd::HitType>>(testtree, "TRDHit");
+  testresult.print();
+}
+
+// Simple macro to compare properties of simulated hits
+// of two runs (for example reference run and test run).
+// Used for instance to validate MC optimizations or to study
+// the effect on the physics at the lowest level.
+void compareHits(const char* reffilename = "o2sim.root", const char* testfilename = "o2sim.root")
+{
+  TFile rf(reffilename, "OPEN");
+  TFile tf(testfilename, "OPEN");
+  auto reftree = (TTree*)rf.Get("o2sim");
+  auto testtree = (TTree*)tf.Get("o2sim");
+
+  checkITS(reftree, testtree);
+  checkTOF(reftree, testtree);
+  checkEMC(reftree, testtree);
+  checkTRD(reftree, testtree);
+}