LBNFDeDxMap.cc

Go to the documentation of this file.

#include <sstream>
#include "LBNFDeDxMap.hh"
#include "LBNEVolumePlacements.hh"

LBNFDeDxMap::LBNFDeDxMap(const std::string &aName):
fName(aName),
fPhiKeyFact(10.0), // minimum step size ih = 2*PI/10.; 
fZKeyFact(1.0), // minimum step size in Z = 1 mm 
fRKeyFact(10.0), // minimum step size in R = 100 microns 
fPhiKeyMult(100000), // maximum number of Z step corresponds ~ 100 m., cell size 1 mm
fNumPhiSlices(6)
{

  if (aName == std::string("HornAConceptDesignJan2017")) {
    defineInnerCHornA();
  }  
  if (aName == std::string("HornBConceptDesignJan2017")) {
    defineInnerCHornB();
  }  
}
// based on the geometry defined in 
// LBNEVolumePlacements::PlaceFinalLBNFConceptHornA
// and 
void LBNFDeDxMap::defineInnerCHornA() {
  
  const double phiWidth = 2.0*M_PI/fNumPhiSlices;
  //
  // The IO transition section, upstream. First 
  //
  LBNFDeDxCell aCellT; 
  aCellT.G4Name = std::string("LBNFConceptHornAUpstrIOSect");
  aCellT.MARSName = std::string("hAuA_");
  aCellT.rWidth = 16.666666;
  aCellT.rCenter = 43. + 0.5*aCellT.rWidth; // defining boundaries.. 
  aCellT.zWidth = 250.; // oversized, for this volume, limits are in r,phi. 
                       // encompass the 6+ 3 cells for this region. 
  aCellT.zCenter = -125.; // crude, but O.K. 
  aCellT.phiWidth = phiWidth;
  aCellT.keyedByZ = false;
  aCellT.eDep = 0.; aCellT.eDepSq = 0.; aCellT.nEntries = 0;
  
  for (size_t kr=0; kr != 9; kr++) {
    if (kr == 3)  aCellT.MARSName = std::string("hAuB_");
    LBNFDeDxCell aCellR(aCellT);
    aCellR.rCenter = aCellT.rCenter + kr*aCellT.rWidth;
    std::ostringstream mNameStrStr; mNameStrStr << kr;
    aCellR.MARSName += std::string(mNameStrStr.str());                 
    for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellR);
      aCell.phiCenter = 0.5*aCellR.phiWidth + kPhi*aCellR.phiWidth; 
      unsigned int tag = this->getTagFromR(aCell.rCenter, aCellT.rCenter + 0.00001*aCellT.rWidth, 
                                                aCellT.rWidth,  aCell.phiCenter);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag 
//                << " " <<  aCell.MARSName << " r " << aCell.rCenter <<  std::endl;
      fMap.emplace(tag, aCell);
    }
  }
  rStarts.push_back(aCellT.rCenter); // Only approximate, two step sizes are used above... 
  rWidths.push_back(aCellT.rWidth);
  G4NamesInMapRTag.push_back(aCellT.G4Name);
  //
  // The Inner Conductor.. Note: this scheme prevents us to look at the outer conductor..
  // Which is O.K., I was asked to look at the heating of the inner conductor only...
  // Actually, we could use special cases, like for the water layer..  
  //
  aCellT.G4Name = std::string("LBNFConceptHornAICCyl");
  aCellT.MARSName = std::string("hAtz_");
  aCellT.rWidth = 1.001; // oversize by 1 microns.. 
  aCellT.rCenter = 44.; // non-defining boundaries.. 
  aCellT.zWidth = 19.833861; // exact... 
  aCellT.zCenter = 0.024 +  0.5*aCellT.zWidth;
  aCellT.phiWidth = phiWidth;
  aCellT.keyedByZ = true;
  aCellT.eDep = 0.;
  for (size_t kz=0; kz != 59; kz++) {
    LBNFDeDxCell aCellZ(aCellT);
    aCellZ.zCenter = aCellT.zCenter + kz*aCellT.zWidth;
    std::ostringstream mNameStrStr; mNameStrStr << kz;
    aCellZ.MARSName = std::string("hAtz_");
    aCellZ.MARSName += std::string(mNameStrStr.str());                 
    for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellZ);
      aCell.phiCenter = 0.5*aCellZ.phiWidth + kPhi*aCellZ.phiWidth; 
      unsigned int tag = this->getTagFromZ(aCell.zCenter, aCellT.zCenter + 0.00001*aCellT.zWidth, 
                                  aCellT.zWidth,  aCell.phiCenter);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << " " 
//                << aCell.MARSName << " z " << aCell.zCenter << std::endl;
      fMap.emplace(tag, aCell);
      // Add the Water layer (for extra-credits !);
      LBNFDeDxCell aCellW(aCell);
      aCellW.G4Name = std::string("LBNFConceptHornAICCylWater");
      aCellW.rCenter = 45.5;
      aCellW.rWidth = 0.5;
      tag += 1000000;
      fMap.emplace(tag, aCellW);
    }
  } 
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(std::string("LBNFConceptHornAICCyl"));
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(std::string("LBNFConceptHornAICCylWater"));
  // 
  // The tapered section.. 
  //
  aCellT.G4Name = std::string("LBNFConceptHornAICTaper");
  aCellT.MARSName = std::string("hAtz_");
  aCellT.rWidth = 1.2; // largely oversized.. don't matter
  aCellT.rCenter = 44.; // non-defining boundaries.. 
  aCellT.zWidth = 19.64469; // exact 
                       // ???? The length of the tapered section is 992.3 mm.. 
                       // perhaps a typo in the document.. 
  aCellT.zWidth = 19.84469;                    
  aCellT.zCenter = 1170.01 +  0.5*aCellT.zWidth;
  aCellT.phiWidth = phiWidth;
  aCellT.keyedByZ = true;
  aCellT.eDep = 0.;
  size_t nTapered = 50;
//
// We change the cell size, to check the integral.. And we found a serious bug in the cell key find algo!  
//
//  aCellT.zWidth *= 1.5152;
//  nTapered = 33; 
//  aCellT.zCenter = 1170.01 +  0.5*aCellT.zWidth;
  for (size_t kz=0; kz != nTapered; kz++) {
    LBNFDeDxCell aCellZ(aCellT);
    aCellZ.zCenter = aCellT.zCenter + kz*aCellT.zWidth;
    aCellZ.rCenter = aCellT.rCenter  - aCellZ.zWidth*kz*1.00776e-2; // the tapered slope..                     
    std::ostringstream mNameStrStr; mNameStrStr << kz;
    aCellZ.MARSName = std::string("hAds2z_");
    aCellZ.MARSName += std::string(mNameStrStr.str());
    for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellZ);
      aCell.phiCenter = 0.5*aCellZ.phiWidth + kPhi*aCellZ.phiWidth; 
      unsigned int tag = this->getTagFromZ(aCell.zCenter,aCellT.zCenter + 0.00001*aCellT.zWidth, 
                                                 aCellT.zWidth,  aCell.phiCenter);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << aCell.MARSName <<  std::endl;
      fMap.emplace(tag, aCell);
      // Add the Water layer (for extra-credits !);
      LBNFDeDxCell aCellW(aCell);
      aCellW.G4Name = std::string("LBNFConceptHornAICTaperWater");
      aCellW.rCenter = aCellZ.rCenter + 1.0 + 0.5;
      aCellW.rWidth = 0.5;
      tag += 1000000;
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << aCellW.MARSName <<  std::endl;
      fMap.emplace(tag, aCellW);
    }
  } 
 
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(std::string("LBNFConceptHornAICTaper"));
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(std::string("LBNFConceptHornAICTaperWater"));
  //
  // The downstream I/O section.. 
  //
  aCellT.G4Name = std::string("LBNFConceptHornADownstrIOSect");
  aCellT.MARSName = std::string("hA_dbB_");
  aCellT.rWidth = 19.11484;
  aCellT.rCenter = 33. + 0.5*aCellT.rWidth; // defining boundaries.. 
  aCellT.zWidth = 86.; // oversized, slightly... for this volume, limits are in r,phi. 
                       // encompass the 8 cells for this region. 
  aCellT.zCenter = 0.024 + 1170 + 992.8 + 0.5*aCellT.zWidth; // Fairly precise, to avoid overlap
  // with the last segment of the tapered section..  
  aCellT.phiWidth = phiWidth;
  aCellT.keyedByZ = false;
  aCellT.eDep = 0.;
  for (size_t kr=0; kr != 9; kr++) { 
    LBNFDeDxCell aCellR(aCellT);
    aCellR.rCenter = aCellT.rCenter + kr*aCellT.rWidth;
    std::ostringstream mNameStrStr; mNameStrStr << kr;
    aCellR.MARSName += std::string(mNameStrStr.str());                 
    for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellR);
      aCell.phiCenter = 0.5*aCellR.phiWidth + kPhi*aCellR.phiWidth; 
      unsigned int tag = this->getTagFromR(aCell.rCenter, aCellT.rCenter + 0.00001*aCellT.rWidth, 
                                           aCellT.rWidth, aCell.phiCenter);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << aCell.MARSName <<  std::endl;
      fMap.emplace(tag, aCell);
    }
  }
  rStarts.push_back(aCellT.rCenter);
  rWidths.push_back(aCellT.rWidth);
  G4NamesInMapRTag.push_back(aCellT.G4Name);
  //
  // Special volumes, Spider support ring holder. 
  //
  aCellT.G4Name = std::string("LBNFConceptHornARingSuppIC");
  aCellT.MARSName = std::string("hAcoll");
  aCellT.rWidth = 3.75;
  aCellT.rCenter = 45.01 + 0.5*aCellT.rWidth; // defining boundaries.. 
  aCellT.zWidth = 20.0;
  aCellT.zCenter = 1084.5; // defining boundaries.. 
  aCellT.keyedByZ = true;
  for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellT);
      aCell.phiCenter = 0.5*aCell.phiWidth + kPhi*aCell.phiWidth; 
      unsigned int tag = 2000000 + this->getTagFromZ(aCell.zCenter,aCellT.zCenter + 0.00001*aCellT.zWidth, 
                                                             aCellT.zWidth,  aCell.phiCenter);
      fMap.emplace(tag, aCell);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << aCell.MARSName <<  std::endl;
  }
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(aCellT.G4Name);
   //
  // Special volumes, Weld 
  //
  aCellT.G4Name = std::string("LBNFConceptHornAWeldIC0toIC1");
  aCellT.MARSName = std::string("hAcoll2");
  aCellT.rWidth = 0.625;
  aCellT.rCenter = 45.01 + 0.5*aCellT.rWidth; // defining boundaries.. 
  aCellT.zWidth = 23.3;
  aCellT.zCenter = 1140.065; // defining boundaries.. 
  aCellT.keyedByZ = true;
  for (size_t kPhi=0; kPhi != fNumPhiSlices; kPhi++) {
      LBNFDeDxCell aCell(aCellT);
      aCell.phiCenter = 0.5*aCell.phiWidth + kPhi*aCell.phiWidth; 
      unsigned int tag = 3000000 + this->getTagFromZ(aCell.zCenter,aCellT.zCenter + 0.00001*aCellT.zWidth, 
                                   aCellT.zWidth,  aCell.phiCenter);
      fMap.emplace(tag, aCell);
//      std::cerr << " LBNFDeDxMap::defineInnerCHornA, tag " << tag << aCell.MARSName <<  std::endl;
  }
  zStarts.push_back(aCellT.zCenter);
  zWidths.push_back(aCellT.zWidth);
  G4NamesInMapZTag.push_back(aCellT.G4Name);
  //
  // Done.... 
  //
  std::cerr << " LBNFDeDxMap::defineInnerCHornA with " << fMap.size() << " cells, keep going " << std::endl;
}
void LBNFDeDxMap::defineInnerCHornB() {

}  
void LBNFDeDxMap::fill(const G4Step* aStep) {
   
   if (aStep->GetStepLength() < 1.0e-6) return;
   G4StepPoint* prePtr = aStep->GetPreStepPoint();
   G4StepPoint* postPtr = aStep->GetPostStepPoint();
   if (prePtr == 0) return;
   if (postPtr == 0) return;
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   if (volPre == 0) return;
   std::string volNamePre(volPre->GetName());
   double zStart = std::nan("");
   double rStart = std::nan("");
   double zWidth = std::nan("");
   double rWidth = std::nan("");
   bool isRelevantRTag = false;
//   std::cerr << " LBNFDeDxMap::fill, number of RVols name " << G4NamesInMapRTag.size() << std::endl;
   for (size_t kN=0; kN != G4NamesInMapRTag.size(); kN++) {
     std::string volR = G4NamesInMapRTag[kN];
//     std::cerr << "  ... volR " << volR << " volNamePre " << volNamePre << std::endl;
     if (volNamePre.find(volR.c_str()) != std::string::npos) {
       isRelevantRTag = true;
       rStart = rStarts[kN];
       rWidth = rWidths[kN];
       break;
     }
   }
   bool isRelevantZTag = false;
//   std::cerr << " LBNFDeDxMap::fill, number of ZVols name " << G4NamesInMapZTag.size() << std::endl;
   for (size_t kN=0; kN != G4NamesInMapZTag.size(); kN++) {
     std::string volZ = G4NamesInMapZTag[kN];
//     std::cerr << "  ... volZ " << volZ << " volNamePre " << volNamePre << std::endl;
     if (volNamePre.find(volZ.c_str()) != std::string::npos) {
       isRelevantZTag = true; 
       zStart = zStarts[kN];
       zWidth = zWidths[kN];
       break;
     }
   }
   if ((!isRelevantRTag) && (!isRelevantZTag)) return;
//   std::cerr << " LBNFDeDxMap::fill, volume " << volNamePre << " is relevant.. " << std::endl;
//   if (isRelevantRTag) std::cerr << " R tag type... " << std::endl;
//   if (isRelevantZTag) std::cerr << " Z tag type... " << std::endl;
   
   // Assume the track is a straight line.  Accurate enough in the material of interest... 
   const double x = 0.5*(prePtr->GetPosition()[0] + postPtr->GetPosition()[0]);
   const double y = 0.5*(prePtr->GetPosition()[1] + postPtr->GetPosition()[1]);
   const double z = 0.5*(prePtr->GetPosition()[2] + postPtr->GetPosition()[2]);
   const double r = std::sqrt(x*x + y*y);
   double phi = std::atan2(y,x);
   if (phi < 0.) phi += 2.0*M_PI;
   unsigned int tag = 0;
   if (isRelevantRTag) tag = getTagFromR(r, rStart, rWidth, phi);
   if (isRelevantZTag) tag = getTagFromZ(z, zStart, zWidth, phi);
   // Special cases..
   if (volNamePre.find("Water") != std::string::npos) tag += 1000000;  
   if (volNamePre.find("RingSuppIC") != std::string::npos) tag += 2000000;  
   if (volNamePre.find("WeldIC0toIC1") != std::string::npos) tag += 3000000;  
   std::map<unsigned int, LBNFDeDxCell>::iterator it = fMap.find(tag);
   if (it == fMap.end()) {
       //
       // This sometimes occur at thedge of G4 volume. Ignore, small effect. 
//       std::cerr << " LBNFDeDxMap::fill, in vol  " << volNamePre 
//                 << " at r/z " << r << " / " << z  << " phi " << phi << " tag " << tag 
//               <<  " not a valid cell " << std::endl;
//      if (isRelevantRTag) std::cerr << " R-Map, rStart " << rStart << " rWidth " << rWidth << std::endl;       
//      if (isRelevantZTag) std::cerr << " Z-Map, zStart " << zStart << " zWidth " << zWidth << std::endl;       
      return; 
   }
   const double edd = aStep->GetTotalEnergyDeposit();
   const G4Track *aTr = aStep->GetTrack();
   it->second.eDep += edd;
   it->second.eDepSq += edd*edd;
   const int pid= aTr->GetParticleDefinition()->GetPDGEncoding();
   const double zPre = prePtr->GetPosition()[2];
   const double zPost = postPtr->GetPosition()[2];
   const double rPre = std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]);
   const double rPost = std::sqrt(postPtr->GetPosition()[0]*postPtr->GetPosition()[0] +
                                 postPtr->GetPosition()[1]*postPtr->GetPosition()[1]);
//   std::cerr << " DeDxMap, trkID   " << aTr->GetTrackID() <<  " length " << aStep->GetStepLength() 
//           << " at Pre, RZ  " << rPre << " / " << zPre << " Post " << rPost << " / " << zPost  << std::endl 
//             << " .............. E dep " << aStep->GetTotalEnergyDeposit() << std::endl;  
   if (std::abs(aTr->GetParticleDefinition()->GetPDGCharge()) > 1.0e-30) it->second.nEntries += 1;
   if ((tag == 2072) || (tag == 102072) || (tag == 202072) || (tag == 302072) || 
       (tag == 402072) || (tag == 402072)) {
     if ((fVectSpecialCell.size() < 1000000) && 
         (std::abs(aTr->GetParticleDefinition()->GetPDGCharge()) > 1.0e-30)) {
       LBNFDeDxEntryInCell aa;
       aa.tag = tag;
       aa.pid = pid;
       aa.length = aStep->GetStepLength();
       aa.rPre = rPre;
       aa.zPre = zPre;
       aa.rPost = rPost;
       aa.zPost = zPost;
       aa.eDep = edd;
       aa.eKin = aTr->GetKineticEnergy();
       fVectSpecialCell.push_back(aa);
     }                     
     
   }
   // 
   // What about the track coming from the inside of the horn, going back inside the IC ??
   //
   
}
void LBNFDeDxMap::dumpASCII(const std::string &aNameDir, int nEvts) {
  std::string aNameOut(aNameDir); 
  aNameOut += std::string("/OutDeDxMap_") + fName + std::string(".txt");
  std::ofstream fOut(aNameOut.c_str());
  fOut << " tag ZOrR G4Name MARSName zC zW rC rW  phiC phiW nEnt eDep sigEDep " << std::endl;
  for (std::map<unsigned int, LBNFDeDxCell>::const_iterator it=fMap.begin(); 
           it != fMap.end(); it++) {
           const double sEdd = it->second.eDep/nEvts;
           double sigEDep = 1.412*it->second.eDep/nEvts;
           if (it->second.nEntries > 1) { 
             const double mEdd = it->second.eDep/it->second.nEntries;
             const double w = ((it->second.eDepSq) - 
                                 (it->second.nEntries) * mEdd * mEdd)/(it->second.nEntries -1);
             sigEDep = std::sqrt(std::abs(w)) * it->second.nEntries/(nEvts*std::sqrt(it->second.nEntries));
           }             
           fOut << " " << it->first << " " <<it->second.keyedByZ 
                << " " << it->second.G4Name << " " << it->second.MARSName 
                << " " << it->second.zCenter << " " << it->second.zWidth 
                << " " << it->second.rCenter << " " << it->second.rWidth 
                << " " << it->second.phiCenter << " " << it->second.phiWidth 
                << " " << it->second.nEntries << " " << sEdd << " " <<  sigEDep << std::endl; 
  }        
  fOut.close();
  std::string aNameOut2(aNameDir); 
  aNameOut2 += std::string("/OutDeDxSpCell_") + fName + std::string(".txt");
  std::ofstream fOut2(aNameOut2.c_str());
  fOut2 << " k  tag pid l rPre zPre rPos zPos eDep eKin " << std::endl;  
  size_t kk = 0; 
  for (std::vector<LBNFDeDxEntryInCell>::const_iterator it=fVectSpecialCell.begin(); 
           it != fVectSpecialCell.end(); it++, kk++) {
           fOut2 << " " << kk << " " << it->tag << " " << it->pid  
                << " " << it->length << " " << it->rPre 
                << " " << it->zPre << " " << it->rPost 
                << " " << it->zPost << " " << it->eDep 
                << " " << it->eKin  << std::endl; 
  }        
  fOut2.close();
}