BesMucNoise.cc

Go to the documentation of this file.

//
//
//
//
 
#include "MucSim/BesMucNoise.hh"
#include "SimUtil/ReadBoostRoot.hh"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
 
#include "G4Box.hh"
#include "G4Trap.hh"
#include "G4VSolid.hh"
#include "Randomize.hh"
#include "math.h"
#include "stdlib.h"
#include "strstream"
#include <fstream>
#include <iostream>
#include <sstream>
 
using namespace std;
 
// const int m_kPart              = 3;
const int BesMucNoise::m_kSegment[m_kPart]  = { 4, 8, 4 };
const int BesMucNoise::m_kAbsorber[m_kPart] = { 9, 9, 9 };
const int BesMucNoise::m_kGap[m_kPart]      = { 8, 9, 8 };
const int BesMucNoise::m_kPanel[m_kPart]    = { 4, 3, 4 };
// const int m_kGasChamber        = 2;
 
BesMucNoise* BesMucNoise::fPointer = 0;
BesMucNoise* BesMucNoise::Instance( void ) {
  if ( !fPointer ) fPointer = new BesMucNoise();
  return fPointer;
}
 
BesMucNoise::BesMucNoise() {
  if ( fPointer )
  {
    G4cout << "BesMucNoise constructed twice." << G4endl;
    exit( -1 );
  }
  fPointer = this;
}
 
BesMucNoise::~BesMucNoise() {
  if ( m_ptrIdTr != NULL ) delete m_ptrIdTr;
}
 
void BesMucNoise::Initialize( G4String filename, G4LogicalVolume* logicalMuc, G4String temp ) {
  // m_noiseLevel = 2; // set to 2 temporary. In future, expected to get from cal servise,
  // like m_noiseLevel = m_ptrCalibSvc->getNoiseLevel();
  /*
    m_noiseLevel = m_ptrCalibSvc->getLevel();
 
       if(m_noiseLevel!=2 && m_noiseLevel!=3){
        G4cout<<"ERROR, In BesMucNoise, noise level should be 2 or 3 ..."<<G4endl;
        m_noiseLevel = 2;
       }
 
       for(int i = 0; i < m_kPart; i++) {
        for(int j = 0; j < m_kSegment[i]; j++) {
        for(int k = 0; k < m_kGap[i]; k++) {
          m_noise[i][j][k] = m_ptrCalibSvc->getBoxCnt(i,j,k);
          int stripNum = 0;
          if(i != 1) stripNum = 64; //endcap;
          else{
            if(k%2 == 0) stripNum = 48;
            if(k%2 == 1 && j != 2) stripNum = 96;
            if(k%2 == 1 && j == 2) stripNum = 112;
          }
 
          for(int strip = 0; strip < stripNum; strip++){
            m_noise_strip[i][j][k][strip] = m_ptrCalibSvc->getStripCnt(i,j,k,strip);
          }
        }
       }
     }
   */
}
 
void BesMucNoise::Initialize( G4String filename, G4LogicalVolume* logicalMuc ) {
  for ( G4int part = 0; part < 3; part++ )
  {
    for ( G4int seg = 0; seg < 8; seg++ )
    {
      for ( G4int gap = 0; gap < 9; gap++ )
      {
        m_noise[part][seg][gap] = 1; // init -1
      }
    }
  }
  // read noise from file
  G4cout << "filename: " << filename << G4endl;
  std::ifstream fin( filename );
  if ( !fin ) { G4cout << "error opening muc_noise data" << G4endl; }
  char buffer[200];
  fin.getline( buffer, 200, '\n' ); // get info whether add noise or not!
  std::istringstream stringBuf( buffer );
 
  // get strip area from logicalMuc!
  G4int tot_NoDaughter = logicalMuc->GetNoDaughters();
  // G4cout<<"---in noise::Init()---  "<<tot_NoDaughter<<G4endl;
  for ( G4int i = 0; i < tot_NoDaughter; i++ )
  {
    G4LogicalVolume* i_LogicalGap = logicalMuc->GetDaughter( i )->GetLogicalVolume();
    G4String         i_GapName    = i_LogicalGap->GetName();
 
    if ( i_GapName.find( "G" ) == 8 )
    {
      G4LogicalVolume* i_LogicalBox = i_LogicalGap->GetDaughter( 0 )->GetLogicalVolume();
      G4LogicalVolume* i_LogicalStripPlane =
          i_LogicalBox->GetDaughter( 0 )->GetLogicalVolume();
      // G4cout<<"---in noise::Init()---  "<<i<<"   "<<i_GapName<<"
      // "<<i_LogicalStripPlane->GetDaughter(1)->GetLogicalVolume()->GetName()<<G4endl;
      G4String strPart = i_GapName.substr( 5, 1 );
      G4String strSeg  = i_GapName.substr( 7, 1 );
      G4String strGap  = i_GapName.substr( 9, 1 );
 
      std::istrstream partBuf( strPart.c_str(), strlen( strPart.c_str() ) );
      std::istrstream segBuf( strSeg.c_str(), strlen( strSeg.c_str() ) );
      std::istrstream gapBuf( strGap.c_str(), strlen( strGap.c_str() ) );
 
      G4int part, seg, gap;
 
      partBuf >> part;
      segBuf >> seg;
      gapBuf >> gap;
      // G4cout<<"-------in noise::Init()----  "<<part<<" "<<seg<<"  "<<gap<<"
      // "<<i_LogicalStripPlane->GetNoDaughters()<<G4endl;
      G4int tot_NoStrip       = i_LogicalStripPlane->GetNoDaughters();
      area[part][seg][gap][0] = tot_NoStrip; // the first element is the total strip number
      G4float tot_Area        = 0;
      // get width between two strip!
      G4LogicalVolume* i_LogicalStrip1 =
          i_LogicalStripPlane->GetDaughter( 1 )->GetLogicalVolume();
      G4LogicalVolume* i_LogicalStrip2 =
          i_LogicalStripPlane->GetDaughter( 2 )->GetLogicalVolume();
      G4Box* temp1;
      G4Box* temp2;
 
      temp1          = (G4Box*)i_LogicalStrip1->GetSolid();
      temp2          = (G4Box*)i_LogicalStrip2->GetSolid();
      G4float Width1 = temp1->GetXHalfLength() * 2;
      G4float Width2 = temp2->GetXHalfLength() * 2;
      G4float pos1   = i_LogicalStripPlane->GetDaughter( 1 )->GetObjectTranslation().x();
      G4float pos2   = i_LogicalStripPlane->GetDaughter( 2 )->GetObjectTranslation().x();
      if ( ( part == 1 && gap % 2 != 0 ) || ( part != 1 && gap % 2 == 0 ) )
      {
        Width1 = temp1->GetYHalfLength() * 2;
        Width2 = temp2->GetYHalfLength() * 2;
        pos1   = i_LogicalStripPlane->GetDaughter( 1 )->GetObjectTranslation().y();
        pos2   = i_LogicalStripPlane->GetDaughter( 2 )->GetObjectTranslation().y();
      }
      G4float width_between_strip = pos2 - pos1 - Width1 / 2 - Width2 / 2;
      // G4cout<<"-----what I want---"<<width_between_strip<<"  "<<pos2<<"  "<<pos1<<"
      // "<<Width1/2<<"
      // "<<Width2/2<<G4endl;
 
      for ( G4int j = 0; j < tot_NoStrip; j++ )
      {
        G4LogicalVolume* i_LogicalStrip =
            i_LogicalStripPlane->GetDaughter( j )->GetLogicalVolume();
        G4Box* temp;
        temp           = (G4Box*)i_LogicalStrip->GetSolid();
        G4float Width  = temp->GetXHalfLength() * 2;
        G4float Length = temp->GetYHalfLength() * 2;
        if ( ( part == 1 && gap % 2 != 0 ) || ( part != 1 && gap % 2 == 0 ) )
        {
          Width  = temp->GetYHalfLength() * 2;
          Length = temp->GetXHalfLength() * 2;
        } // end if
        //       G4cout<<"----in noise::init()----  "<<i_LogicalStrip->GetName()<<" "<<Width<<"
        //       "<<Length<<G4endl;
        if ( j == 0 || j == ( tot_NoStrip - 1 ) ) Width = Width + width_between_strip / 2;
        else Width = Width + width_between_strip;
        G4float Strip_Area            = fabs( Width * Length );
        tot_Area                      = tot_Area + Strip_Area;
        area[part][seg][gap][j + 1]   = tot_Area;
        strip_area[part][seg][gap][j] = Strip_Area;
        // G4cout<<" id: ( "<<part<<" "<<seg<<" "<<gap<<") "<<Strip_Area<<" , w: "<<Width<<" L:
        // "<<Length<<G4endl;
      } // end for(j)
      // unitary
 
      box_area[part][seg][gap] = tot_Area;
 
      for ( G4int k = 1; k < tot_NoStrip + 1; k++ )
      {
        area[part][seg][gap][k] = area[part][seg][gap][k] / tot_Area;
        //       G4cout<<"----in noise::init---"<<area[part][seg][gap][k]<<G4endl;
      } // end for(k)
 
    } // end if(gap)
  }   // end for(i)
 
  //  G4cout<<"---in noise::init()---area--"<<G4endl;
  // for(G4int kk=0;kk<96;kk++){
  // G4cout<<area[1][0][0][kk+1]-area[1][0][0][kk]<<" ,";
  // }
 
  // init MucIdTransform and const pointer
  m_ptrIdTr = new MucIdTransform();
  CheckCalibSvc();
  InitProb();
}
 
void BesMucNoise::CheckCalibSvc() {
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode   sc         = svcLocator->service( "MucCalibConstSvc", m_ptrCalibSvc, true );
 
  if ( sc != StatusCode::SUCCESS ) { G4cout << "Can not use MucCalibConstSvc!" << G4endl; }
}
 
void BesMucNoise::InitProb() {
  m_HitMean = 0.08;
  int NHIT  = 20;
 
  for ( int i = 0; i < NHIT; i++ ) { m_Prob[i][0] = m_Prob[i][1] = 0.; }
 
  double sum = 0.;
  double EU  = TMath::Power( TMath::Exp( 1.0 ), -m_HitMean );
  for ( int i = 0; i < NHIT; i++ )
  {
    m_Prob[i][0] = EU * TMath::Power( m_HitMean, i ) / TMath::Factorial( i );
    sum += m_Prob[i][0];
    m_Prob[i][1] = sum;
    // G4cout << i << "\tProb:\t" << m_Prob[i][0] << "\tSumProb:\t" << m_Prob[i][1] << G4endl;
  }
}
 
G4int BesMucNoise::AddNoise( int model, BesMucHitsCollection* aMucHitCollection,
                             BesMucHitsCollection* aMucHitList ) {
  G4int noiseNum = 0;
  if ( model == 1 ) noiseNum = NoiseByCnt( aMucHitCollection, aMucHitList );
  else noiseNum = NoiseByNosRatio( aMucHitCollection, aMucHitList );
  // G4cout << "Noise hit(s) produced:\t" << noiseNum << "\tby model:\t" << model << G4endl;
  return noiseNum;
}
 
G4int BesMucNoise::NoiseByCnt( BesMucHitsCollection* MucHitCollection,
                               BesMucHitsCollection* MucHitList ) {
  G4int noiseNum = 0;
  m_noiseLevel   = m_ptrCalibSvc->getLevel();
 
  for ( int i = 0; i < m_kPart; i++ )
  {
    for ( int j = 0; j < m_kSegment[i]; j++ )
    {
      for ( int k = 0; k < m_kGap[i]; k++ )
      {
        if ( m_noiseLevel == 2 )
        {
          G4int hitNum = NoiseSampling( m_noiseLevel, i, j, k, 0 );
          for ( G4int ii = 0; ii < hitNum; ii++ )
          {
            G4int      strip         = GetStripNo( i, j, k );
            BesMucHit* noiseHit      = new BesMucHit( i, j, k, strip, -1, -1 );
            BesMucHit* noiseHit2     = new BesMucHit( i, j, k, strip, -1, -1 );
            bool       noiseHitExist = false;
            noiseHitExist            = IsExist( noiseHit, MucHitList );
            MucHitCollection->insert( noiseHit );
            if ( !noiseHitExist )
            {
              MucHitList->insert( noiseHit2 );
              noiseNum += 1;
              // G4cout<<"-------insert-------  "<<" "<<noiseHit2->GetStrip()<<" "<<i<<"
              // "<<j<<" "<<k<<" "<<l<<"
              // "<<m<<G4endl;
            }
            else delete noiseHit2;
          } // hitNum
        }   // box level
 
        if ( m_noiseLevel == 3 )
        {
          int stripNum = m_ptrIdTr->GetStripMax( i, j, k );
          for ( int strip = 0; strip < stripNum; strip++ )
          {
            G4int hitNum = NoiseSampling( m_noiseLevel, i, j, k, strip );
            if ( hitNum > 0 )
            {
              BesMucHit* noiseHit      = new BesMucHit( i, j, k, strip, -1, -1 );
              BesMucHit* noiseHit2     = new BesMucHit( i, j, k, strip, -1, -1 );
              bool       noiseHitExist = false;
              noiseHitExist            = IsExist( noiseHit, MucHitList );
              MucHitCollection->insert( noiseHit );
              if ( !noiseHitExist )
              {
                MucHitList->insert( noiseHit2 );
                noiseNum += 1;
                // G4cout<<"-------insert-------  "<<" "<<noiseHit2->GetStrip()<<" "<<i<<"
                // "<<j<<" "<<k<<" "<<l<<"
                // "<<m<<G4endl;
              }
              else delete noiseHit2;
            }
          } // stripNum
        }   // strip level
      }     // part
    }       // segment
  }         // layer
 
  return noiseNum;
}
 
G4int BesMucNoise::NoiseByNosRatio( BesMucHitsCollection* MucHitCollection,
                                    BesMucHitsCollection* MucHitList ) {
  G4int   noiseNum = 0;
  G4float random   = 0;
 
  // sample noise hit number;
  random = G4UniformRand();
  for ( int i = 0; i < 20; i++ )
  {
    if ( random < m_Prob[i][1] )
    {
      noiseNum = i;
      break;
    }
  }
  // G4cout << "Random for noiseNum:\t" << random << "\tProbSum:\t" << m_Prob[noiseNum][1] <<
  // "\t" << noiseNum << G4endl;
 
  int     prt, seg, lay, str, tmp_strip;
  G4float nosRatio = 0.;
 
  for ( int i = 0; i < noiseNum; i++ )
  {
    do {
      random    = G4UniformRand();
      tmp_strip = TMath::Nint( random * STRIP_MAX ); // 9152
      m_ptrIdTr->SetStripPos( tmp_strip, &prt, &seg, &lay, &str );
      nosRatio = m_ptrCalibSvc->getNosRatio( prt, seg, lay, str );
      random   = G4UniformRand();
      if ( random < nosRatio ) break;
    } while ( 1 );
 
    // G4cout<<"Strip:\t"<<prt<<"\t"<<seg<<"\t"<<lay<<"\t"<<str<<"\t"<<random<<"
    // "<<nosRatio<<G4endl;
 
    BesMucHit* noiseHit  = new BesMucHit( prt, seg, lay, str, -1, -1 );
    BesMucHit* noiseHit2 = new BesMucHit( prt, seg, lay, str, -1, -1 );
 
    bool noiseHitExist = false;
    noiseHitExist      = IsExist( noiseHit, MucHitList );
    MucHitCollection->insert( noiseHit );
    if ( !noiseHitExist ) MucHitList->insert( noiseHit2 );
    else delete noiseHit2;
  } // noiseNum;
 
  return noiseNum;
}
 
bool BesMucNoise::IsExist( BesMucHit* aNoiseHit, BesMucHitsCollection* aMucHitList ) {
  bool  isExist = false;
  G4int n_hit   = aMucHitList->entries();
  for ( G4int iNoise = 0; iNoise < n_hit; iNoise++ )
  {
    if ( aNoiseHit->GetTrackIndex() % 1000 ==
             ( *aMucHitList )[iNoise]->GetTrackIndex() % 1000 &&
         aNoiseHit->GetPart() == ( *aMucHitList )[iNoise]->GetPart() &&
         aNoiseHit->GetSeg() == ( *aMucHitList )[iNoise]->GetSeg() &&
         aNoiseHit->GetGap() == ( *aMucHitList )[iNoise]->GetGap() &&
         aNoiseHit->GetStrip() == ( *aMucHitList )[iNoise]->GetStrip() )
    {
      isExist = true;
      break;
    }
  }
 
  return isExist;
}
 
G4int BesMucNoise::NoiseSampling( int level, int prt, int seg, int lay, int strip ) {
  G4int    hitNum = 0;
  G4double lambda;
  // must get the time of this event
  G4double t = 800; // 800ns  temporary!
  G4double e = 2.71828182845904590;
 
  if ( level == 2 )
    lambda = m_ptrCalibSvc->getBoxCnt( prt, seg, lay ) * box_area[prt][seg][lay] * 1E-2 *
             1E-9; // 1E-2:mm2->cm2 1E-9:ns
  else if ( level == 3 )
    lambda = m_ptrCalibSvc->getStripCnt( prt, seg, lay, strip ) *
             strip_area[prt][seg][lay][strip] * 1E-2 * 1E-9; // 1E-2:mm2->cm2 1E-9:ns
  else lambda = 0.;
 
  // sample noise hit number;
  G4float  random = G4UniformRand(); //***use other random
  G4double prob   = 0;
  do {
    prob = prob + pow( e, -lambda * t ) * pow( lambda * t, hitNum ) / Factorial( hitNum );
    if ( random < prob ) break;
    hitNum++;
  } while ( 1 );
 
  // G4cout<<" hitNum: "<<hitNum<<G4endl;
  return hitNum;
}
 
G4float BesMucNoise::Factorial( G4int i ) {
  G4float fact = 1;
  if ( i == 0 || i == 1 ) return 1;
  else
  {
    for ( G4int ii = 2; ii <= i; ii++ ) { fact = fact * ii; }
    return fact;
  }
}
 
G4int BesMucNoise::GetStripNo( G4int part, G4int seg, G4int gap ) {
  G4int stripno;
  // G4float random=G4UniformRand();//***use other random
  G4float random = ( rand() % 100000 ) / 100000.0;
  if ( part == 1 )
  { // berrel
    G4float width = area[part][seg][gap][3] - area[part][seg][gap][2];
    stripno       = G4int( ( random - area[part][seg][gap][1] ) / width ) + 2;
    if ( stripno < 1 ) stripno = 1;
    if ( stripno > 111 ) stripno = 111;
    // G4cout<<"---in noise::GetStripNo()---stripno= "<<stripno<<"
    // "<<(random-area[part][seg][gap][1])/width<<G4endl;
    G4int step = IsNearestStrip( stripno, part, seg, gap, random );
    while ( step != 0 )
    {
      // G4cout<<"---in noise::GetStripNo()---while "<<G4endl;
      stripno += step;
      step = IsNearestStrip( stripno, part, seg, gap, random );
    }          // endl while
    stripno--; // the first element is total strip number
    return stripno;
  }
  else
  { // endcap
    G4int max, min, mid, pass;
    min = 1;
    max = area[part][seg][gap][0];
    mid = G4int( ( min + max ) / 2 );
    // G4cout<<"---in noise first---"<<min<<" "<<mid<<" "<<max<<G4endl;
    do {
      //  G4cout<<"-----in noise---"<<random<<" "<<area[part][seg][gap][min]<<"
      //  "<<area[part][seg][gap][mid]<<"
      //  "<<area[part][seg][gap][max]<<G4endl;
      pass = 0;
      if ( random > area[part][seg][gap][mid] )
      {
        min = mid;
        mid = G4int( ( min + max ) / 2 );
      }
      else if ( random < area[part][seg][gap][mid - 1] )
      {
        if ( random < area[part][seg][gap][1] )
        {
          pass = 1;
          mid  = 1;
          max  = 1;
        }
        else
        {
          max = mid - 1;
          mid = G4int( ( min + max ) / 2 );
        }
      }
      else { pass = 1; }
 
      if ( min == mid ) mid = max;
      // G4cout<<"-----in noise---"<<min<<" "<<mid<<" "<<max<<G4endl;
    } while ( pass == 0 && ( max > mid && mid > min ) );
    // G4cout<<"-----in noise---  "<<mid-1<<G4endl;
    return mid - 1;
  }
}
 
G4int BesMucNoise::IsNearestStrip( G4int stripno, G4int part, G4int seg, G4int gap,
                                   G4float random ) {
  if ( stripno == 1 ) { return 0; }
  else
  {
    if ( area[part][seg][gap][stripno] != 0 )
    {
      //   G4cout<<"--in noise::IsNearestStrip()--  "<<area[part][seg][gap][stripno]<<"
      //   "<<area[part][seg][gap][stripno-1]<<" "<<random<<G4endl;
      if ( random <= area[part][seg][gap][stripno] &&
           random > area[part][seg][gap][stripno - 1] )
        return 0;
      if ( random <= area[part][seg][gap][stripno - 1] ) return -1;
      if ( random > area[part][seg][gap][stripno] ) return 1;
    }
    else { return -1; }
  } // end else
  std::cerr << __FILE__ << ":" << __LINE__ << ": should not reach here" << std::endl;
  exit( 1 );
}