sEPDMappingScanner/src/CheckEpdMap.cc

0001 #include "CheckEpdMap.h"
0002
0003 #include <phool/PHCompositeNode.h>
0004 #include <phool/getClass.h>
0005 #include <phool/recoConsts.h>
0006 #include <phool/PHNodeIterator.h>
0007
0008 #include <cdbobjects/CDBTTree.h>
0009 #include <ffamodules/CDBInterface.h>
0010
0011 #include <calobase/TowerInfoDefs.h>
0012 #include <calobase/TowerInfoContainerv4.h>
0013 #include <epd/EpdGeomV2.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <TFile.h>
0018 #include <TTree.h>
0019 #include <TProfile2D.h>
0020 #include <TH1.h>
0021 #include <TH2.h>
0022 #include <TCanvas.h>
0023 #include <TLatex.h>
0024 #include <TColor.h>
0025 #include <TStyle.h>
0026
0027 #include <iostream>
0028 #include <iomanip>
0029 #include <cmath>
0030 #include <ffamodules/CDBInterface.h>
0031
0032 #include <g4main/PHG4HitContainer.h>
0033 #include <g4main/PHG4HitDefs.h>
0034 #include <g4main/PHG4Hit.h>
0035 #include <g4main/PHG4InEvent.h>
0036 #include <g4main/PHG4Particlev2.h>
0037 #include <epd/EPDDefs.h>
0038
0039
0040 //_____________________________________________________________________________
0041 CheckEpdMap::CheckEpdMap(const std::string& name)
0042   : SubsysReco(name)
0043 {}
0044
0045 //_____________________________________________________________________________
0046 int CheckEpdMap::InitRun(PHCompositeNode* topNode)
0047 {
0048   std::cout << "Initializing CheckEpdMap" << std::endl;
0049   geom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0050   if (!geom)
0051   {
0052     std::cerr << Name() << ": EpdGeom node TOWERGEOM_EPD not found\n";
0053     return Fun4AllReturnCodes::ABORTEVENT;
0054   }
0055
0056   std::cout << "Found TOWERGEOM_EPD" << std::endl;
0057
0058   auto [vkey, mapped_values] = build_key_vector();
0059   //auto vkey, mapped_values = build_key_vector();
0060   if (vkey.empty())
0061   {
0062     std::cerr << Name() << ": failed to build vkey vector\n";
0063     return Fun4AllReturnCodes::ABORTEVENT;
0064   }
0065
0066   check_map(vkey, mapped_values);
0067
0068   m_keyVec = std::move(vkey); //store key vector for later use
0069
0070   const char* title[2][2] = {
0071     { "South  old (key vector)", "North old (key vector)" },
0072     { "South  new (geom loop)" , "North  new (geom loop)" }
0073   };
0074
0075   for (int arm = 0; arm < 2; ++arm)
0076   {
0077     m_profOld[arm] = new TProfile2D(Form("profOld_%d",arm),  title[0][arm],
0078                                     NRING, -0.5, NRING-0.5,
0079                                     NPHI , -0.5, NPHI -0.5);
0080     m_profNew[arm] = new TProfile2D(Form("profNew_%d",arm),  title[1][arm],
0081                                     NRING, -0.5, NRING-0.5,
0082                                     NPHI , -0.5, NPHI -0.5);
0083   }
0084
0085   m_diffCount = new TH1I("diffCount",
0086                        "index where old (channel) =/= new (key);channel index",
0087                        744, -0.5, 743.5);
0088
0089
0090   m_diffADC   = new TH2F("diffADC",
0091                         "ADC(old) – ADC(new);ring;#phi sector",
0092                         NRING, -0.5, NRING-0.5,
0093                         NPHI , -0.5, NPHI -0.5);
0094
0095   m_hitCenterDist = new TH1F(
0096     "hitCenterDist",
0097     "Distance between true hit and tile center;#Delta r (cm)",
0098     400,   // bins
0099     0.0,   // min
0100     200.0  // max
0101   );
0102
0103   constexpr int NBINS = 1100;
0104   constexpr float RANGE = 110.;
0105
0106   for (unsigned arm = 0; arm < 2; ++arm) {
0107     for (unsigned r = 0; r < NRING; ++r)
0108     {
0109       const unsigned maxPhi = (r == 0 ? 12 : 24);
0110       for (unsigned phi = 0; phi < maxPhi; ++phi)
0111       {
0112           m_tileHist[arm][r][phi] = new TH2F(
0113             Form("tileHit_arm%u_ring%u_phi%u", arm, r, phi),
0114             Form("Tile hits: arm=%u ring=%u phi=%u;X–X_{ctr} (cm);Z–Z_{ctr} (cm)",
0115                 arm, r, phi),
0116             NBINS, -RANGE, RANGE,
0117             NBINS, -RANGE, RANGE);
0118
0119       }
0120
0121     }
0122   }
0123
0124
0125
0126   BuildSimMap();
0127   InitializeSimMatchMap();
0128
0129
0130
0131
0132   if (!m_simulationMode){
0133     towers = findNode::getClass<TowerInfoContainer>(topNode,
0134       "TOWERINFO_CALIB_SEPD");
0135     if (!towers) return Fun4AllReturnCodes::ABORTEVENT;
0136   }
0137
0138
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 //_____________________________________________________________________________
0143 int CheckEpdMap::process_event(PHCompositeNode* topNode)
0144 {
0145
0146   if (! m_simulationMode) {
0147     return processEventDataMode(topNode);
0148   }
0149   else {
0150     return processEventSimulationMode(topNode);
0151   }
0152
0153
0154
0155 }
0156
0157 //_____________________________________________________________________________
0158 std::pair<std::vector<unsigned>, std::vector<int>>
0159 CheckEpdMap::build_key_vector() const
0160 {
0161
0162   std::cout << "CheckEpdMap::build_key_vector called..." << std::endl;
0163
0164   static const std::string sEPDMapName   = "SEPD_CHANNELMAP";   // domain
0165   static const std::string sEPDFieldName = "epd_channel_map";   // branch
0166
0167   std::string calibFile = CDBInterface::instance()->getUrl(sEPDMapName);
0168   if (calibFile.empty())
0169   {
0170     std::cerr << "CheckEpdMap : no sEPD mapping file for domain "
0171               << sEPDMapName << " found – aborting\n";
0172     std::exit(1);
0173   }
0174
0175   std::cout << "Found sEPD mapping file" << std::endl;
0176
0177   std::unique_ptr<CDBTTree> cdbttree( new CDBTTree(calibFile) );
0178
0179   std::cout << "Created CDBTTree..." << std::endl;
0180
0181
0182   std::vector<unsigned> vkey;
0183   std::vector<int> mapped_vals;
0184   vkey.reserve(768);
0185
0186   std::cout << "Building key vector... " << std::endl;
0187
0188   for (int ch = 0; ch < 768; ++ch)
0189   {
0190     int keymap = cdbttree->GetIntValue(ch, sEPDFieldName);
0191
0192     if (keymap == 999)
0193       continue;
0194
0195     unsigned key = TowerInfoDefs::encode_epd(static_cast<unsigned>(keymap));
0196     vkey.push_back(key);
0197     mapped_vals.push_back(keymap);
0198   }
0199
0200   std::cout << "Key vector built with " << vkey.size() << " keys" << std::endl;
0201
0202   return { std::move(vkey), std::move(mapped_vals) };
0203 }
0204
0205 //_____________________________________________________________________________
0206 void CheckEpdMap::check_map(
0207   const std::vector<unsigned>& vkey, const std::vector<int> mapped_vals) const
0208 {
0209 std::cout << "CheckEpdMap::check_map called..." << std::endl;
0210
0211 //constexpr unsigned NRING = 16;
0212 //constexpr unsigned NPHI  = 24;
0213 int reverse[2][NRING][NPHI];
0214 std::fill_n(&reverse[0][0][0], 2*NRING*NPHI, -1);
0215
0216   std::cout << "tower index that maps to zeroth key: " << mapped_vals[0] << std::endl;
0217
0218 for (unsigned ch = 0; ch < vkey.size(); ++ch)
0219 {
0220   unsigned key   = vkey[ch];
0221   //unsigned key = TowerInfoDefs::encode_epd(ch);
0222   //int mapped_ch = mapped_vals[ch]; //this is the "channel number" that should be returned by decode_epd
0223
0224
0225
0226   int mapped_ch = ch;
0227   //we are passing mapped_ch to get_tower_at_ch, so this is what we should map to
0228
0229
0230 unsigned arm   = TowerInfoDefs::get_epd_arm   (key);
0231 unsigned rbin  = TowerInfoDefs::get_epd_rbin  (key);
0232 unsigned phibin= TowerInfoDefs::get_epd_phibin(key);
0233
0234 if (arm < 2 && rbin < NRING && phibin < NPHI)
0235 reverse[arm][rbin][phibin] = mapped_ch;
0236 else
0237 std::cerr << Name() << ": decoded bins out of range for ch=" << ch << "\n";
0238 }
0239
0240 std::memcpy(m_reverse,reverse, sizeof m_reverse); //store truth table for later
0241
0242 std::cout << "Truth table built!" << std::endl;
0243
0244
0245   std::size_t nBad = 0;
0246   for (unsigned arm = 0; arm < 2; ++arm) {
0247     for (unsigned r = 0; r < NRING; ++r)
0248     {
0249       unsigned maxPhi = (r == 0 ? 12 : 24);
0250       for (unsigned phi = 0; phi < maxPhi; ++phi)
0251       {
0252         int ch_truth = reverse[arm][r][phi];
0253
0254         //to generate map directly from (arm,r,phi) -> tower index/"channel" number
0255
0256         unsigned key_tmp = TowerInfoDefs::encode_epd(arm, r, phi);
0257         int ch_code      = TowerInfoDefs::decode_epd(key_tmp);
0258
0259         if (ch_truth != ch_code)
0260         {
0261           std::cout << Name() << ": mismatch arm="<<arm
0262                     << " r="<<r<<" phi="<<phi
0263                     << "  code="<<ch_code
0264                     << "  truth="<<ch_truth << '\n';
0265           ++nBad;
0266         }
0267       }
0268     }
0269
0270   }
0271
0272
0273   std::cout << Name() << ": total mismatches = " << nBad << std::endl;
0274
0275 }
0276
0277 int CheckEpdMap::End(PHCompositeNode*)
0278 {
0279
0280
0281   TH2F* h[2][2];   // [row][col] : row 0 = OLD, row 1 = NEW ; col 0 = SOUTH, col 1 = NORTH
0282   for (int arm = 0; arm < 2; ++arm)
0283   {
0284     h[0][arm] = (TH2F*)m_profOld[arm]->ProjectionXY(
0285                   Form("hOld_%d",arm));
0286     h[1][arm] = (TH2F*)m_profNew[arm]->ProjectionXY(
0287                   Form("hNew_%d",arm));
0288   }
0289
0290   TCanvas c("c","sEPD compare maps",1800,900);
0291   c.Divide(2,2);
0292
0293   TLatex lab; lab.SetTextAlign(22); lab.SetTextSize(0.025);
0294
0295   auto draw = [&](TH2F* hh, int pad, bool isNorth, bool isOld)
0296   {
0297   c.cd(pad);
0298   hh->SetMinimum(0);
0299   hh->Draw("colz");
0300
0301   const unsigned arm = isNorth ? 1 : 0;
0302
0303   for (int bx = 1; bx <= hh->GetNbinsX(); ++bx)
0304     for (int by = 1; by <= hh->GetNbinsY(); ++by)
0305     {
0306       unsigned ring = bx - 1;
0307       unsigned phi  = by - 1;
0308       if (ring == 0 && phi > 11) continue;   // dummy bins in innermost ring
0309
0310       int chan;
0311       if (isOld)           // ------- “old” (CDB‑based) mapping -------
0312       {
0313         chan = m_reverse[arm][ring][phi];
0314         //let's brute force a check here
0315         unsigned key  = TowerInfoDefs::encode_epd(arm, ring, phi);
0316         int geom_chan = TowerInfoDefs::decode_epd(key);
0317         if(chan != geom_chan){
0318           std::cout << "Mismatched channel found: (arm, r, phi) = (" << arm << ", " <<
0319           ring << ", " << phi << "). reverse channel mapping: " << chan << ", new channel mapping: " <<
0320           geom_chan << std::endl;
0321         }
0322       }
0323       else                 // ------- “new” (geometry) mapping -------
0324       {
0325         unsigned key  = TowerInfoDefs::encode_epd(arm, ring, phi);
0326         chan = TowerInfoDefs::decode_epd(key);
0327         //and another one
0328         int old_chan = m_reverse[arm][ring][phi];
0329         if(chan != old_chan){
0330           std::cout << "Mismatched channel found: (arm, r, phi) = (" << arm << ", " <<
0331           ring << ", " << phi << "). reverse channel mapping: " << old_chan << ", new channel mapping: " <<
0332           chan << std::endl;
0333         }
0334       }
0335
0336       lab.DrawLatex(hh->GetXaxis()->GetBinCenter(bx),
0337                     hh->GetYaxis()->GetBinCenter(by),
0338                     Form("%u", chan));
0339     }
0340   };
0341
0342
0343   //north: h[0][new/old]
0344
0345   draw(h[0][1], 1, true ,  true );   // north – old
0346
0347   draw(h[1][1], 2, true ,  false);   // north – new
0348
0349   draw(h[0][0], 3, false,  true );   // south – old
0350
0351   draw(h[1][0], 4, false,  false);   // south – new
0352
0353   c.SaveAs("sEPD_compare_channelMap.pdf");
0354
0355   TCanvas cDiffADC("cDiffADC","ADC(old) - ADC(new)",800,700);
0356   cDiffADC.SetRightMargin(0.12);
0357   m_diffADC->Draw("colz");
0358   cDiffADC.Print("diffADC.pdf");
0359
0360   TCanvas cDiffIdx("cDiffIdx","channel indices with ADC mismatch",800,600);
0361   m_diffCount->Draw();
0362   cDiffIdx.Print("diffCount.pdf");
0363
0364   PrintReverseTable();
0365
0366   TCanvas cDist("cDist","Δ(hit,center)",800,600);
0367   m_hitCenterDist->Draw();
0368   cDist.Print("hitCenterDistance.pdf");
0369
0370   std::cout << "Attempting to write output histsogram " << outputFile << std::endl;
0371
0372   if (! outputFile.empty())
0373   {
0374     TFile fout(outputFile.c_str(), "RECREATE");
0375     if (fout.IsOpen())
0376     {
0377       m_hitCenterDist->Write("hitCenterDist");
0378       fout.Close();
0379       std::cout << Name()
0380                 << ": wrote histogram to "
0381                 << outputFile << std::endl;
0382     }
0383     else
0384     {
0385       std::cerr << Name()
0386                 << ": ERROR opening output file "
0387                 << outputFile << std::endl;
0388     }
0389   }
0390
0391   //example: draw tile (0,10,7)
0392
0393   TCanvas cTile("cTile","example tile: (0,10,7)", 800,600);
0394   m_tileHist[0][10][7]->Draw("COLZ");
0395   cTile.Print("example_tile.pdf");
0396
0397
0398   Int_t palette[2] = { kRed, kGreen };
0399   gStyle->SetPalette(2, palette);
0400   gStyle->SetNumberContours(2);
0401
0402
0403   TCanvas cSimMatchPolar("cSimMatchPolar",
0404                         "sEPD Simulation Match (polar)",-1,0,
0405                         1400, 600);
0406   cSimMatchPolar.SetRightMargin(0.15);
0407   cSimMatchPolar.Divide(2,1);
0408
0409   cSimMatchPolar.cd(1);
0410   gPad->SetTicks(1,1);
0411
0412   gPad->DrawFrame(-3.8, -3.8, 3.8, 3.8);
0413   m_simMatchPolarS ->Draw("same COLZ pol AH");
0414   m_simMatchPolarS01->Draw("same COLZ pol AH");
0415   TLatex ts; ts.SetNDC(); ts.SetTextSize(0.05);
0416   ts.DrawLatex(0.30, 0.92, "sEPD South Sim-Match");
0417
0418
0419   cSimMatchPolar.cd(2);
0420   gPad->SetTicks(1,1);
0421   gPad->DrawFrame(-3.8, -3.8, 3.8, 3.8);
0422   m_simMatchPolarN ->Draw("same COLZ pol AH");
0423   m_simMatchPolarN01->Draw("same COLZ pol AH");
0424   TLatex tn; tn.SetNDC(); tn.SetTextSize(0.05);
0425   tn.DrawLatex(0.30, 0.92, "sEPD North Sim-Match");
0426
0427
0428   cSimMatchPolar.SaveAs("sEPD_sim_match_polar.pdf");
0429
0430
0431
0432   //c.Write();
0433
0434   return Fun4AllReturnCodes::EVENT_OK;
0435 }
0436
0437
0438 //_____________________________________________________________________________
0439 void CheckEpdMap::PrintReverseTable(std::ostream &os /* = std::cout */) const
0440 {
0441   using std::setw;
0442
0443   os << "\nstatic constexpr int epdchnlmap[2]"
0444         "[" << NRING << "][" << NPHI << "] = {\n";
0445
0446   for (unsigned arm = 0; arm < 2; ++arm)
0447   {
0448     os << "  { // arm " << arm << '\n';
0449     for (unsigned r = 0; r < NRING; ++r)
0450     {
0451       os << "    { ";
0452       for (unsigned phi = 0; phi < NPHI; ++phi)
0453       {
0454         int ch = m_reverse[arm][r][phi];
0455         if (ch < 0) ch = 999;          // mark unused slots
0456         os << setw(4) << ch;
0457         if (phi + 1 < NPHI) os << ", ";
0458       }
0459       os << " }";
0460       if (r + 1 < NRING) os << ',';
0461       os << '\n';
0462     }
0463     os << "  }";
0464     if (arm == 0) os << ',';
0465     os << '\n';
0466   }
0467   os << "};\n" << std::endl;
0468 }
0469
0470
0471
0472 int CheckEpdMap::processEventDataMode(PHCompositeNode * /*topNode*/)
0473 {
0474
0475   for (unsigned int ch = 0; ch < towers->size(); ++ch)
0476   {
0477     unsigned key = m_keyVec[ch];
0478     unsigned arm    = TowerInfoDefs::get_epd_arm   (key);
0479     unsigned rbin   = TowerInfoDefs::get_epd_rbin  (key);
0480     unsigned phibin = TowerInfoDefs::get_epd_phibin(key);
0481     auto *twr_ch = towers->get_tower_at_channel(ch);
0482     m_profOld[arm]->Fill(rbin, phibin, twr_ch->get_energy());
0483   }
0484
0485
0486   for (unsigned arm = 0; arm < 2; ++arm)
0487   {
0488     for (unsigned r = 0; r < NRING; ++r)
0489     {
0490       unsigned maxPhi = (r == 0 ? 12 : 24);
0491       for (unsigned phi = 0; phi < maxPhi; ++phi)
0492       {
0493         unsigned key = TowerInfoDefs::encode_epd(arm, r, phi);
0494         if (auto *twr_geom = towers->get_tower_at_key(key))
0495         {
0496           m_profNew[arm]->Fill(r, phi, twr_geom->get_energy());
0497         }
0498       }
0499     }
0500   }
0501
0502   for (unsigned arm  = 0; arm < 2;      ++arm)
0503   for (unsigned ring = 0; ring < NRING; ++ring)
0504   {
0505     unsigned maxPhi = (ring == 0 ? 12 : 24);
0506     for (unsigned phi = 0; phi < maxPhi; ++phi)
0507     {
0508       int ch_idx = m_reverse[arm][ring][phi];
0509       if (ch_idx < 0) continue;
0510       float adc_old = towers->get_tower_at_channel(ch_idx)->get_energy();
0511
0512       uint32_t key_new = TowerInfoDefs::encode_epd(arm, ring, phi);
0513       auto *tw_new     = towers->get_tower_at_key(key_new);
0514       float adc_new    = tw_new->get_energy();
0515
0516       if (adc_old != adc_new)
0517       {
0518         m_diffCount->Fill(ch_idx);
0519         m_diffADC  ->Fill(ring, phi, adc_old - adc_new);
0520       }
0521     }
0522   }
0523
0524   return Fun4AllReturnCodes::EVENT_OK;
0525 }
0526
0527 void CheckEpdMap::BuildSimMap() {
0528   //build map from (arm,sector,tile_idx) -> (arm,rbin,phibin)
0529   //taken from PHG4EPDModuleReco
0530   for (unsigned int arm = 0; arm < 2; arm++) //arm
0531   {
0532     for (int sec = 0; sec < 12; sec++) //sectors
0533     {
0534       for (int tile_id = 0; tile_id < 31; tile_id++) //tile_idx
0535       {
0536         unsigned int phi = Getphimap(tile_id) + 2 * sec;
0537         unsigned int r = Getrmap(tile_id);
0538
0539         if (r == 0)
0540         {
0541           phi = sec;
0542         }
0543
0544         m_tile_bin[arm][sec][tile_id] = std::make_tuple(arm, r, phi);
0545         m_module_coords[arm][r][phi] = std::make_tuple(arm,sec,tile_id);
0546
0547       }
0548     }
0549   }
0550 }
0551
0552
0553 int CheckEpdMap::Getphimap(int phiindex)
0554 {
0555   static const int phimap[31] = {0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
0556
0557   return phimap[phiindex];
0558 }
0559
0560 int CheckEpdMap::Getrmap(int rindex)
0561 {
0562   static const int rmap[31] = {0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15};
0563
0564   return rmap[rindex];
0565 }
0566
0567
0568
0569 int CheckEpdMap::processEventSimulationMode(PHCompositeNode *topNode)
0570 {
0571   auto *g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_EPD");
0572   if (!g4hit)
0573   {
0574     std::cerr << Name() << ": no G4HIT_EPD node found in sim‐mode\n";
0575     return Fun4AllReturnCodes::ABORTEVENT;
0576   }
0577
0578
0579
0580   auto* simTowers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_SIM_EPD");
0581   if (!simTowers)
0582   {
0583     std::cerr << Name() << ": no TOWERINFO_SIM_EPD node found in sim-mode\n";
0584     return Fun4AllReturnCodes::ABORTEVENT;
0585   }
0586
0587
0588
0589   auto range = g4hit->getHits();
0590
0591   // skip the event if there are no g4hits
0592   if (range.first == range.second) {
0593     return Fun4AllReturnCodes::EVENT_OK;
0594   }
0595
0596
0597   uint32_t arm   = 0;
0598   uint32_t ring  = 0;
0599   uint32_t phi   = 0;
0600   bool     found_match = false;
0601
0602
0603   for (auto iter = range.first; iter != range.second; ++iter) {
0604     PHG4Hit* hit        = iter->second;
0605     uint32_t sim_tileid = hit->get_hit_id() >> PHG4HitDefs::hit_idbits;
0606
0607     //get coords
0608     int this_tile   = EPDDefs::get_tileid(sim_tileid);
0609     int this_sector = EPDDefs::get_sector(sim_tileid);
0610     int this_arm    = EPDDefs::get_arm(sim_tileid);
0611
0612
0613     std::tie(arm, ring, phi)
0614       = m_tile_bin[this_arm][this_sector][this_tile];
0615
0616
0617     unsigned key = TowerInfoDefs::encode_epd(arm, ring, phi);
0618     auto*   twr = simTowers->get_tower_at_key(key);
0619     if (!twr) {
0620       std::cout << "WARNING: no tower found at location (arm,ring,phi) = (" <<
0621       arm << ", " << ring << ", " << phi << ")." << std::endl;
0622       //Fun4AllReturnCodes::ABORTEVENT;
0623     }
0624     float   e   = twr ? twr->get_energy() : 0;
0625
0626     if (e == -1.f) {
0627       std::cout << "Found match!" << std::endl;
0628       found_match = true;
0629       break;
0630     }
0631   }
0632
0633   std::cout << "Setting histogram bin (" << ring << ", " <<
0634   phi << ") to " << found_match << " with arm = " << arm << std::endl;
0635
0636   double w = ( found_match ? +1.0 : 0.1 );
0637
0638
0639
0640   int rbin = ring +1;
0641   int phibin = phi +1;
0642
0643   std::cout << "phibin" << phibin << " rbin " << rbin << " w = " << w << std::endl;
0644
0645   if (ring == 0)
0646   {
0647     if (arm == 0) m_simMatchPolarS01->SetBinContent(phibin, rbin, w);
0648     else          m_simMatchPolarN01->SetBinContent(phibin, rbin, w);
0649   }
0650   else
0651   {
0652     if (arm == 0) m_simMatchPolarS->SetBinContent(phibin, rbin, w);
0653     else          m_simMatchPolarN->SetBinContent(phibin, rbin, w);
0654   }
0655
0656   std::cout << "Filled histogram!" << std::endl;
0657
0658
0659
0660   return Fun4AllReturnCodes::EVENT_OK;
0661 }
0662
0663 void CheckEpdMap::InitializeSimMatchMap()
0664 {
0665   std::cout << "Booking sim histos!" << std::endl;
0666   if (!m_simMatchPolarS)
0667   {
0668     std::cout << "initializing m_simMatchPolarS" << std::endl;
0669
0670     m_simMatchPolarS = new TH2F(
0671       "h_simMatch_south", "Sim-match South;phi;ring",
0672       24, 0, 2*M_PI,
0673       16, 0.15, 3.5
0674     );
0675     m_simMatchPolarS->SetMaximum(1);
0676     m_simMatchPolarS->SetMinimum(0);
0677   }
0678   if (!m_simMatchPolarN)
0679   {
0680     m_simMatchPolarN = new TH2F(
0681       "h_simMatch_north", "Sim-match North;phi;ring",
0682       24, 0, 2*M_PI,
0683       16, 0.15, 3.5
0684     );
0685     m_simMatchPolarN->SetStats(0);
0686     m_simMatchPolarN->SetMaximum(1);
0687     m_simMatchPolarN->SetMinimum(0);
0688   }
0689
0690   if (!m_simMatchPolarS01)
0691   {
0692     m_simMatchPolarS01 = new TH2F(
0693       "h_simMatch_south01","Sim-match South inner;phi;ring",
0694       12, 0, 2*M_PI,
0695       16, 0.15, 3.5
0696     );
0697     m_simMatchPolarS01->SetStats(0);
0698     m_simMatchPolarS01->SetMaximum(1);
0699     m_simMatchPolarS01->SetMinimum(0);
0700   }
0701   if (!m_simMatchPolarN01)
0702   {
0703     m_simMatchPolarN01 = new TH2F(
0704       "h_simMatch_north01","Sim-match North inner;phi;ring",
0705       12, 0, 2*M_PI,
0706       16, 0.15, 3.5
0707     );
0708     m_simMatchPolarN01->SetStats(0);
0709     m_simMatchPolarN01->SetMaximum(1);
0710     m_simMatchPolarN01->SetMinimum(0);
0711   }
0712
0713   std::cout << "InitializeSimMatchMap: Done!" << std::endl;
0714
0715 }