Merge pull request #29 from blackcathj/SPACALTower-NewEval-Tower

Build SPACAL towers as in engineering drawing

simulation/g4simulation/g4cemc/RawTower.h

@@ -20,6 +20,8 @@ class RawTower : public PHObject {
   virtual int get_bineta() const { PHOOL_VIRTUAL_WARN("get_ieta()"); return -1; }
   virtual int get_binphi() const { PHOOL_VIRTUAL_WARN("get_iphi()");; return -1; }
   virtual float get_energy() const { PHOOL_VIRTUAL_WARN("get_energy()"); return 0.0; }
+  virtual void set_light_yield(float l)  { PHOOL_VIRTUAL_WARN("set_light_yield()"); return ; }
+  virtual float get_light_yield() const { PHOOL_VIRTUAL_WARN("get_light_yield()"); return 0.0; }
 
   virtual bool is_adjacent(RawTower& tower) { PHOOL_VIRTUAL_WARNING; return false; }
 

simulation/g4simulation/g4cemc/RawTowerBuilder.cc

@@ -102,6 +102,7 @@ RawTowerBuilder::process_event(PHCompositeNode *topNode)
 
       if (verbosity > 2)
         {
+          std::cout << PHWHERE << " print out the cell:" << std::endl;
           cell->identify();
         }
 
@@ -113,6 +114,7 @@ RawTowerBuilder::process_event(PHCompositeNode *topNode)
 	   _towers->AddTower(cell->get_binz(), cell->get_binphi(), tower);
 	 }
        tower->add_ecell(cell->get_cell_id(), cell->get_edep());
+       tower->set_light_yield( tower->get_light_yield() +  cell->get_light_yield()  );
        rawtowergeom =  findNode::getClass<RawTowerGeom>(topNode, TowerGeomNodeName.c_str());
        if (verbosity > 2)
 	 {
@@ -224,6 +226,16 @@ RawTowerBuilder::CreateNodes(PHCompositeNode *topNode)
             {
               _netabins = cellgeo->get_zbins();// bin eta in the same number of z bins
             }
+          else if (_cell_binning == PHG4CylinderCellDefs::spacalbinning)
+            {
+              // use eta definiton for each row of towers
+              _netabins = cellgeo->get_etabins();
+            }
+          else
+            {
+              cout <<"RawTowerBuilder::CreateNodes::"<<Name()
+                  <<" - Fatal Error - unsupported cell binning method "<<_cell_binning<<endl;
+            }
 	  inner_radius = cellgeo->get_radius();
 	  first_layer = cellgeo->get_layer();
 	  ifirst = 0;
@@ -322,7 +334,9 @@ RawTowerBuilder::CreateNodes(PHCompositeNode *topNode)
     }
 
   if (_cell_binning == PHG4CylinderCellDefs::etaphibinning
-      or _cell_binning == PHG4CylinderCellDefs::etaslatbinning)
+      or _cell_binning == PHG4CylinderCellDefs::etaslatbinning
+      or _cell_binning == PHG4CylinderCellDefs::spacalbinning
+      )
     {
 //  rawtowergeom->set_etamin(_etamin);
 //  rawtowergeom->set_etastep(_etastep);

simulation/g4simulation/g4cemc/RawTowerGeomv2.cc

@@ -89,17 +89,48 @@ RawTowerGeomv2::get_etabin(const double eta) const
 
   int ibin = -1;
   int i = 0;
+//  for (bound_map_t::const_iterator iter = eta_bound_map.begin();
+//      iter != eta_bound_map.end(); ++iter)
+//    {
+//      if (eta >= iter->first && eta < iter->second)
+//        {
+//          ibin = i;
+//          break;
+//        }
+//
+//      i++;
+//    }
+
+  // switch to search for the closest bin
+  // since in a realistic calorimeter, there could be gaps
+  double min_deta = 10;
+
   for (bound_map_t::const_iterator iter = eta_bound_map.begin();
       iter != eta_bound_map.end(); ++iter)
     {
+      const double mean_eta = 0.5*(iter->first + iter->second);
+
       if (eta >= iter->first && eta < iter->second)
         {
+          // found the bin that the hit belong
+          min_deta = 0;
           ibin = i;
           break;
         }
+      else
+        {
+          const double deta = fabs(mean_eta - eta);
+          if (deta < min_deta)
+            {
+              min_deta = deta;
+              ibin = i;
+            }// keep searching
+        }
 
       i++;
     }
+
+
   if (ibin < 0 || ibin >= nphibins)
     {
       cout

simulation/g4simulation/g4cemc/RawTowerv1.cc

@@ -15,7 +15,8 @@ RawTowerv1::RawTowerv1() :
 
 RawTowerv1::RawTowerv1(const int ieta, const int iphi) :
   bineta(ieta),
-  binphi(iphi)
+  binphi(iphi),
+  light_yield(0)
 {}
 
 RawTowerv1::~RawTowerv1()

simulation/g4simulation/g4cemc/RawTowerv1.h

@@ -21,6 +21,9 @@ class RawTowerv1 : public RawTower {
   int get_binphi() const { return binphi; }
   float get_energy() const;
 
+  void set_light_yield(float l)  { light_yield = l; }
+  float get_light_yield() const { return light_yield; };
+
   bool is_adjacent(RawTower& tower);
 
   std::pair< std::map<unsigned int,float>::const_iterator, std::map<unsigned int,float>::const_iterator > get_g4cells()
@@ -30,10 +33,11 @@ class RawTowerv1 : public RawTower {
  protected:
   int bineta;
   int binphi;
+  float light_yield;
 
   std::map<unsigned int, float> ecells;
 
-  ClassDef(RawTowerv1,1)
+  ClassDef(RawTowerv1,2)
 };
 
 #endif /* RAWTOWERV1_H_ */

simulation/g4simulation/g4detectors/Makefile.am

@@ -33,6 +33,7 @@ pkginclude_HEADERS = \
   PHG4BlockGeomContainer.h \
   PHG4CylinderCell.h \
   PHG4CylinderCellv1.h \
+  PHG4CylinderCell_Spacalv1.h \
   PHG4CylinderCellv2.h \
   PHG4CylinderCellContainer.h \
   PHG4CylinderGeom.h \
@@ -42,6 +43,7 @@ pkginclude_HEADERS = \
   PHG4CylinderGeomContainer.h \
   PHG4CylinderCellDefs.h \
   PHG4CylinderCellGeom.h \
+  PHG4CylinderCellGeom_Spacalv1.h \
   PHG4CylinderCellGeomContainer.h
 
 libg4detectors_io_la_SOURCES = \
@@ -77,12 +79,16 @@ libg4detectors_io_la_SOURCES = \
   PHG4CylinderCell_Dict.cc \
   PHG4CylinderCellv1.cc \
   PHG4CylinderCellv1_Dict.cc \
+  PHG4CylinderCell_Spacalv1.cc \
+  PHG4CylinderCell_Spacalv1_Dict.cc \
   PHG4CylinderCellv2.cc \
   PHG4CylinderCellv2_Dict.cc \
   PHG4CylinderCellContainer.cc \
   PHG4CylinderCellContainer_Dict.cc \
   PHG4CylinderCellGeom.cc \
   PHG4CylinderCellGeom_Dict.cc \
+  PHG4CylinderCellGeom_Spacalv1.cc \
+  PHG4CylinderCellGeom_Spacalv1_Dict.cc \
   PHG4CylinderCellGeomContainer.cc \
   PHG4CylinderCellGeomContainer_Dict.cc
 
@@ -184,6 +190,8 @@ libg4detectors_la_SOURCES = \
   PHG4SpacalDetector.cc \
   PHG4ProjSpacalDetector.cc \
   PHG4FullProjSpacalDetector.cc \
+  PHG4FullProjSpacalCellReco.cc \
+  PHG4FullProjSpacalCellReco_Dict.cc \
   PHG4SpacalSteppingAction.cc \
   PHG4SpacalSubsystem.cc \
   PHG4SpacalSubsystem_Dict.cc

simulation/g4simulation/g4detectors/PHG4CylinderCell.h

@@ -28,19 +28,25 @@ class PHG4CylinderCell : public PHObject
   }
 
   virtual void add_edep(const PHG4HitDefs::keytype g4hitid, const float edep) {return;}
+  virtual void add_edep(const PHG4HitDefs::keytype g4hitid, const float edep, const float light_yield) {return;}
+
   virtual void set_cell_id(const PHG4CylinderCellDefs::keytype id) {return;}
   virtual void set_layer(const unsigned int i) {return;}
- 
+
   virtual double get_edep() const {return NAN;}
   virtual unsigned int get_layer() const {return 0xFFFFFFFF;}
   virtual PHG4CylinderCellDefs::keytype get_cell_id() const {return 0xFFFFFFFF;}
   virtual int get_binz() const {return -1;}
   virtual int get_binphi() const {return -1;}
   virtual int get_bineta() const {return -1;}
+  virtual float  get_light_yield() const  {    return NAN;  }
+  virtual int    get_fiber_ID() const {return -1;}
 
   virtual void set_phibin(const int i) {return;}
   virtual void set_zbin(const int i) {return;}
   virtual void set_etabin(const int i) {return;}
+  virtual void  set_light_yield(float lightYield)  {   return;  }
+  virtual void  set_fiber_ID(int fiberId) {   return;  }
 
   virtual void set_sensor_index(const std::string &si) {return;}
   virtual std::string get_sensor_index() const {return "";}

simulation/g4simulation/g4detectors/PHG4CylinderCellDefs.h

@@ -6,7 +6,7 @@ namespace PHG4CylinderCellDefs
   typedef unsigned int keytype;
   static int keybits = 8;
   static int cell_idbits = 32 - keybits;
-  enum {undefined = 0, sizebinning = 1, etaphibinning = 2, etaslatbinning = 3};
+  enum {undefined = 0, sizebinning = 1, etaphibinning = 2, etaslatbinning = 3, spacalbinning = 4};
 }
 
 #endif

simulation/g4simulation/g4detectors/PHG4CylinderCellGeom.cc

@@ -151,7 +151,7 @@ PHG4CylinderCellGeom::set_etabins(const int i)
 void
 PHG4CylinderCellGeom::identify(std::ostream& os) const
 {
-  os << "layer: " << layer
+  os << "PHG4CylinderCellGeom::identify - layer: " << layer
      << ", radius: " << radius
      << ", thickness: " << thickness;
   switch (binning)
@@ -168,8 +168,11 @@ PHG4CylinderCellGeom::identify(std::ostream& os) const
       break;
     case PHG4CylinderCellDefs::etaslatbinning:
       os << ", etabins: " << nzbins
-	 << ", etamin: " << zmin
-	 << ", etastepsize: " << zstep;
+   << ", etamin: " << zmin
+   << ", etastepsize: " << zstep;
+      break;
+    case PHG4CylinderCellDefs::spacalbinning:
+      os << ", etabins: " << nzbins<<" for Spacal";
       break;
     default:
       os << "no valid binning method: " << binning << endl;
@@ -317,6 +320,9 @@ PHG4CylinderCellGeom::methodname(const int i) const
     case PHG4CylinderCellDefs::etaslatbinning:
       return "Eta/numslat bins";
       break;
+    case PHG4CylinderCellDefs::spacalbinning:
+      return "SPACAL Tower bins";
+      break;
     default:
       break;
     }
@@ -340,14 +346,16 @@ void
 PHG4CylinderCellGeom::check_binning_method_eta(const std::string & src) const
 {
   if (binning != PHG4CylinderCellDefs::etaphibinning && 
-      binning != PHG4CylinderCellDefs::etaslatbinning)
+      binning != PHG4CylinderCellDefs::etaslatbinning&&
+      binning != PHG4CylinderCellDefs::spacalbinning)
     {
       if (src.size())
         cout << src<<" : ";
 
       cout << "different binning method used " << methodname(binning)
            << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning)
 	   << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning)
+     << " or " << methodname(PHG4CylinderCellDefs::spacalbinning)
            << endl;
       exit(1);
     }
@@ -359,14 +367,17 @@ PHG4CylinderCellGeom::check_binning_method_phi(const std::string & src) const
 {
   if (binning != PHG4CylinderCellDefs::etaphibinning && 
       binning != PHG4CylinderCellDefs::sizebinning &&
-      binning != PHG4CylinderCellDefs::etaslatbinning)
+      binning != PHG4CylinderCellDefs::etaslatbinning &&
+      binning != PHG4CylinderCellDefs::spacalbinning)
     {
       if (src.size())
         cout << src<<" : ";
 
       cout << "different binning method used " << methodname(binning)
            << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning)
-	   << " or " << methodname(PHG4CylinderCellDefs::sizebinning)
+           << " or " << methodname(PHG4CylinderCellDefs::sizebinning)
+           << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning)
+           << " or " << methodname(PHG4CylinderCellDefs::spacalbinning)
            << endl;
       exit(1);
     }

simulation/g4simulation/g4detectors/PHG4CylinderCellGeom.h

@@ -28,16 +28,16 @@ class PHG4CylinderCellGeom: public PHObject
   double get_etastep() const;
   double get_etamin() const;
 
-  std::pair<double, double> get_zbounds(const int ibin) const;
-  std::pair<double, double> get_phibounds(const int ibin) const;
-  std::pair<double, double> get_etabounds(const int ibin) const;
-  double get_etacenter(const int ibin) const;
-  double get_zcenter(const int ibin) const;
-  double get_phicenter(const int ibin) const;
+  virtual std::pair<double, double> get_zbounds(const int ibin) const;
+  virtual std::pair<double, double> get_phibounds(const int ibin) const;
+  virtual std::pair<double, double> get_etabounds(const int ibin) const;
+  virtual double get_etacenter(const int ibin) const;
+  virtual double get_zcenter(const int ibin) const;
+  virtual double get_phicenter(const int ibin) const;
 
-  int get_etabin(const double eta) const;
-  int get_zbin(const double z) const;
-  int get_phibin(const double phi) const;
+  virtual int get_etabin(const double eta) const;
+  virtual int get_zbin(const double z) const;
+  virtual int get_phibin(const double phi) const;
 
    void set_layer(const int i) {layer = i;}
    void set_binning(const int i) {binning = i;}