StBbcSimulationMaker.cxx

#include "StBbcSimulationMaker.h"
#include "g2t/St_g2t_bbc_Module.h"
#include "TDataSetIter.h"
#include "TRandom.h"
#include <vector>
#include "StEvent.h"
#include "StTriggerDetectorCollection.h"
#include "StBbcTriggerDetector.h"
#include "StMessMgr.h"

TRandom BbcRndm = TRandom(0);

ClassImp(StBbcSimulationMaker)

  const float BbcTimingRMS = 900.E-12; 
const uint16_t NPMTsmall1 = 16;
const uint16_t NPMTlarge1 = 8;
const uint16_t NPMT1 = NPMTsmall1+NPMTlarge1;//# of PMTs on one side (East/West)
const uint16_t NPMT2 = 2*NPMT1; // ditto on both sides
const float dE1MIPper_gcm2 = 1.95E-3;  // in GeV/(g/cm**2), for polystyrene
const float PolystereneDensity = 1.032;   // in g/cm**3
const float TyleThickness = 1.; // in cm
const float dE_1MIP = dE1MIPper_gcm2*PolystereneDensity*TyleThickness;
const float NPhotoelectrons_1MIP = 15.; 
const float pC_per_Photoelectron = 0.3;
const short NADCbins = 256; /* bins are numbered from 0 to 255 */
const short NTDCbins = 256; /* ditto */
const float pC_perADCbin = 0.25; /* Central trigger barrel 
                                   Digitizer Boards (CDB) have such 
                                   conversion gain (see STAR NIM) */
const float ADC0 = 0.;  /* beginning of ADC range */
const float s_perTDCbin = .1E-9; /* so that 25.6 ns is full range (guess) */
const float TDC0 = 0.;  /* beginning of TDC range */
const float OuterFactor = 0.8; /* larger scintillators of the outer ring
                                  have less light output by this factor, 
                                  compared with the inner ring, per
                                  equal ionization */
const float SinglePhotoElectronResolution = 0.3; // according to Les Bland
/* Numbering: the real PMT numbering (used in the map) starts from 1.
ALL OTHERS start from 0 !
 */


// Note: BBC new factors are calibrated by 2017 p+p runs
const short ADCBin = 4096;

const float LightFactor[NPMT2] = {1., 1., 0.95, 0.92, 1., 1., 1., 1., 0.9, 0.96, 0.95, 1., 0.95, 0.97, 1., 1.,    // East small 
                                  0.33, 0.298, 0.25, 0.26, 0.22, 0.4, 0.24, 0.26,                                 // East large
                                  1.2, 0.9, 1., 0.6, 1., 1., 1.2, 1.2, 1.25, 0.4, 1., 1., 0.95, 0.95, 0.94, 0.95, // West small
                                  0., 0.243, 0., 0.2, 0.24, 0.18, 0.19, 0.2};                                     // West large

const float TileResolution[NPMT2] = {0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, // East small
                                     0.6, 0.6, 0.4, 0.5, 0.5, 0.6, 0.4, 0.5,                                         // East large
                                     0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.5, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, // West small
                                     0., 0.5, 0., 0.5, 0.6, 0.6, 0.4, 0.65};                                         // West large

const float pC_per_ADCBins[NPMT2] = {0.03, 0.03, 0.033, 0.033, 0.03, 0.03, 0.022, 0.03, 0.021, 0.022, 0.024, 0.03, 0.024, 0.024, 0.025, 0.027, // East small
                                     0.0226, 0.028, 0.021, 0.023, 0.014,  0.0125, 0.012, 0.0125,                                               // East large
                                     0.031, 0.03, 0.03, 0.018, 0.03, 0.03, 0.023, 0.024, 0.028, 0.005, 0.03, 0.03, 0.028, 0.023, 0.022, 0.028, // West small
                                     0., 0.027, 0., 0.06, 0.013, 0.014, 0.0129, 0.02};                                                         // West large

const int shift_ADC0[NPMT2] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   // East small
                               0, 0, 0, 0, -15, 0, -15, 0,                       // East large
                               0, 0, 0, 0, 0, 0, 0, 0, 0, -60, 0, 0, 0, 0, 0, 0, // West small
                               0, 0, 0, 0, -10, -7, -15, -10};                   // West large

//____________________________________________________________________________

bool IsSmall(short iPMT)
{
  if ( 0<=iPMT && iPMT<NPMTsmall1) return 1;
  if ( NPMT1<=iPMT && iPMT<NPMT1+NPMTsmall1) return 1;
  return 0;
}

/*------------------------------------------------------------------------*/
class BbcTOF 
{
private:
  vector<float> Times;

public: 
  BbcTOF():Times(vector<float>(NPMT2)){};  // NPMT1 !
  ~BbcTOF(){};
  void AddTOF(uint16_t ipmt, float time)
  {
    if (Times[ipmt]==0 || Times[ipmt]>time) {Times[ipmt]=time;}
  }
  float GetTOF(uint16_t ipmt)  
  {

    float tof = 0;
    if (Times[ipmt]!=0.){ tof = Times[ipmt]+BbcRndm.Gaus(0.,BbcTimingRMS); }
    return tof;
  }
  short GetTDC(uint16_t ipmt)
  { 
    float T = this->GetTOF(ipmt);

    if (T<TDC0) {return 0;}
    short N = (short)((T-TDC0)/s_perTDCbin);
    if (N>=NTDCbins) {return NTDCbins-1;}
    return N;
  }

};
/*------------------------------------------------------------------------*/
class BbcDE
{
private:
  vector<float> dE;

public:
  BbcDE():dE(vector<float>(NPMT2)){};
  ~BbcDE(){};
  void AddDE(uint16_t ipmt, float de)
  {
    de *= LightFactor[ipmt];
    dE[ipmt] += de;
  }
  float GetDE(uint16_t ipmt)
  {
    float PoissonMean = dE[ipmt]/dE_1MIP*NPhotoelectrons_1MIP;
    short NPhotoelectrons = BbcRndm.Poisson(PoissonMean);
    float Q = pC_per_Photoelectron*
      (1+BbcRndm.Gaus(0.,TileResolution[ipmt]))*NPhotoelectrons;
    return Q;
  }
  short GetADC(uint16_t ipmt)
  {
    float A = this->GetDE(ipmt);
    short N = (short)(A/pC_per_ADCBins[ipmt] + shift_ADC0[ipmt]);
    if (N>=ADCBin) {return ADCBin-1;}
    if (N<0){return 0;}
    return N;
  }
};
//_____________________________________________________________________________

StBbcSimulationMaker::StBbcSimulationMaker(const char *name):StMaker(name)
{
}

//_____________________________________________________________________________

StBbcSimulationMaker::~StBbcSimulationMaker()
{
  //
}


//_____________________________________________________________________________
Int_t StBbcSimulationMaker::Init()
{
  //   BBC tile           BBC PMT number      Comments/description
  // --------           --------------      --------------------
  for (short iew=1; iew<=2; iew++)
    {
      short EW = iew*1000; // for VID
      //      short EWshift = (iew-1)*NPMT1; if Akio numbered from West
      short EWshift = (2-iew)*NPMT1; // for PMT
  //      1                    1              inner small hex tile, 12:00 pos'n 90 degrees
  Geant2PMT[EW+113] = EWshift+1;
  //      2                    2              inner small hex tile,  2:00 pos'n 30 degrees
  Geant2PMT[EW+123] = EWshift+2;
  //      3                    3              inner small hex tile,  4:00 pos'n -30 degrees
  Geant2PMT[EW+133] = EWshift+3;
  //      4                    4              inner small hex tile,  6:00 pos'n -90 degrees
  Geant2PMT[EW+143] = EWshift+4;
  //      5                    5              inner small hex tile,  8:00 pos'n -150 degrees
  Geant2PMT[EW+153] = EWshift+5;
  //      6                    6              inner small hex tile, 10:00 pos'n 150 degrees
  Geant2PMT[EW+163] = EWshift+6;
  //      7                    7              outer small hex tile, 11:00 pos'n 120 degrees
  Geant2PMT[EW+111] = EWshift+7;
  //      8                    8              outer small hex tile, 12:00 pos'n 90 degrees
  Geant2PMT[EW+112] = EWshift+8;
  //      9                    7              outer small hex tile,  1:00 pos'n 60 degrees
  Geant2PMT[EW+121] = EWshift+7;
  //     10                    9              outer small hex tile,  2:00 pos'n 30 degrees
  Geant2PMT[EW+122] = EWshift+9;
  //     11                   10              outer small hex tile,  3:00 pos'n 0 degrees
  Geant2PMT[EW+131] = EWshift+10;
  //     12                   11              outer small hex tile,  4:00 pos'n -30 degrees
  Geant2PMT[EW+132] = EWshift+11;
  //     13                   12              outer small hex tile,  5:00 pos'n -60 degrees
  Geant2PMT[EW+141] = EWshift+12;
  //     14                   13              outer small hex tile,  6:00 pos'n -90 degrees
  Geant2PMT[EW+142] = EWshift+13;
  //     15                   12              outer small hex tile,  7:00 pos'n -120 degrees
  Geant2PMT[EW+151] = EWshift+12;
  //     16                   14              outer small hex tile,  8:00 pos'n -150 degrees
  Geant2PMT[EW+152] = EWshift+14;
  //     17                   15              outer small hex tile,  9:00 pos'n 180 degrees
  Geant2PMT[EW+161] = EWshift+15;
  //     18                   16              outer small hex tile, 10:00 pos'n 150 degrees
  Geant2PMT[EW+162] = EWshift+16;

  //     19                   17              inner large hex tile, 12:00 pos'n 90 degrees
  Geant2PMT[EW+213] = EWshift+17;
  //     20                   18              inner large hex tile,  2:00 pos'n 30 degrees
  Geant2PMT[EW+223] = EWshift+18;
  //     21                   18              inner large hex tile,  4:00 pos'n -30 degrees
  Geant2PMT[EW+233] = EWshift+18;
  //     22                   19              inner large hex tile,  6:00 pos'n -90 degrees
  Geant2PMT[EW+243] = EWshift+19;
  //     23                   20              inner large hex tile,  8:00 pos'n -150 degrees
  Geant2PMT[EW+253] = EWshift+20;
  //     24                   20              inner large hex tile, 10:00 pos'n 150 degrees
  Geant2PMT[EW+263] = EWshift+20;
  //     25                   21              outer large hex tile, 11:00 pos'n 120 degrees
  Geant2PMT[EW+211] = EWshift+21;
  //     26                   21              outer large hex tile, 12:00 pos'n 90 degrees
  Geant2PMT[EW+212] = EWshift+21;
  //     27                   21              outer large hex tile,  1:00 pos'n 60 degrees
  Geant2PMT[EW+221] = EWshift+21;
  //     28                   22              outer large hex tile,  2:00 pos'n 30 degrees
  Geant2PMT[EW+222] = EWshift+22;
  //     29                   22              outer large hex tile,  3:00 pos'n 0 degrees
  Geant2PMT[EW+231] = EWshift+22;
  //     30                   22              outer large hex tile,  4:00 pos'n -30 degrees
  Geant2PMT[EW+232] = EWshift+22;
  //     31                   23              outer large hex tile,  5:00 pos'n -60 degrees
  Geant2PMT[EW+241] = EWshift+23;
  //     32                   23              outer large hex tile,  6:00 pos'n -90 degrees
  Geant2PMT[EW+242] = EWshift+23;
  //     33                   23              outer large hex tile,  7:00 pos'n -120 degrees
  Geant2PMT[EW+251] = EWshift+23;
  //     34                   24              outer large hex tile,  8:00 pos'n -150 degrees
  Geant2PMT[EW+252] = EWshift+24;
  //     35                   24              outer large hex tile,  9:00 pos'n 180 degrees
  Geant2PMT[EW+261] = EWshift+24;
  //     36                   24              outer large hex tile, 10:00 pos'n 150 degrees
  Geant2PMT[EW+262] = EWshift+24;
    }

  // Create Histograms    
#ifdef BbcSimQa
  QaFile = new TFile("StBbcSimQa.root","recreate");
  QaBbcPmtdE = new TH1F("QaBbcPmtdE","BBC PMT",NPMT2,-0.5,NPMT2-0.5);
  QaBbcPmtTime = new TH1F("QaBbcPmtTime","BBC PMT",NPMT2,-0.5,NPMT2-0.5);
  QaBbcEastVid = 
    new TH1F("QaBbcEastVid","BBC PMT with East VID",256,-0.5,255.5);
  QaBbcWestVid = 
    new TH1F("QaBbcWestVid","BBC PMT with West VID",256,-0.5,255.5);
  QaBbcEastPmt =
    new TH1F("QaBbcEastPmt","BBC PMT with East #", 256,-0.5,255.5);
  QaBbcWestPmt =
    new TH1F("QaBbcWestPmt","BBC PMT with West #", 256,-0.5,255.5);
  // create a test map, reverse w.r.t. Geant2PMT 
  typedef map<short,short>::const_iterator CI;

  for (CI I=Geant2PMT.begin(); I!=Geant2PMT.end(); ++I)
    {
      PMT2Geant[(*I).second] = (*I).first;
    }
#endif  

   return StMaker::Init();
}

//_____________________________________________________________________________
Int_t StBbcSimulationMaker::Make()
{

 TDataSet* ds= GetInputDS("geant");
 assert(ds);
 StEvent* event = (StEvent*)GetInputDS("StEvent");
 assert(event);
 St_g2t_ctf_hit* g2t_bbc_hit = (St_g2t_ctf_hit*)ds->Find("g2t_bbc_hit");

 if (g2t_bbc_hit) 
   {
     short nBBChits = g2t_bbc_hit->GetNRows();
     BbcTOF TOFdata;
     BbcDE DEdata;
     g2t_ctf_hit_st *bbc_hit = g2t_bbc_hit->GetTable();
     for (short iBBChit=0; iBBChit<nBBChits; iBBChit++)
     {
       float De = bbc_hit[iBBChit].de;
       float TOF = bbc_hit[iBBChit].tof;
       short Vid =  bbc_hit[iBBChit].volume_id;

       short PMTid = Geant2PMT[Vid];
       if (PMTid == 0) {
         LOG_ERROR << "Cannot find a  PMTid in Geant2PMT for Vid = " << Vid << endm;
         continue;
       }
       PMTid -= 1;
       DEdata.AddDE(PMTid,De);
       TOFdata.AddTOF(PMTid,TOF);
     }

     StTriggerDetectorCollection* myTrig =event->triggerDetectorCollection();
     if (!myTrig) {
       Warning("Make"," NoStTriggerDetectorCollection, Make the new one\n");
       myTrig = new StTriggerDetectorCollection;
       event->setTriggerDetectorCollection(myTrig);
     }
     StBbcTriggerDetector& myBbc = myTrig->bbc();


     for (uint16_t iPMT = 0; iPMT<NPMT2; iPMT++)
       {          
          short ADC = DEdata.GetADC(iPMT);
#ifdef BbcSimQa
          //      QaBbcPmtAdc->Fill(iPMT,ADC);
          short Vid = PMT2Geant[iPMT+1];

          if (Vid<2000) {QaBbcWestVid->Fill(ADC);}
          if (Vid>2000) {QaBbcEastVid->Fill(ADC);}
          if (iPMT<NPMT1) {QaBbcEastPmt->Fill(ADC);}
          if (NPMT1<=iPMT && iPMT<NPMT2) {QaBbcWestPmt->Fill(ADC);}
          QaBbcPmtTime->Fill(iPMT,TOFdata.GetTOF(iPMT));
          QaBbcPmtdE->Fill(iPMT,DEdata.GetDE(iPMT));
#endif
          myBbc.setAdc(iPMT, ADC);
          if (IsSmall(iPMT))
            {
              short TDC = TOFdata.GetTDC(iPMT);
              myBbc.setTdc(iPMT, TDC);
            }
        }
   }
else
  {
    gMessMgr->Info
      ("MLK StBbcSimulationMaker::Make() could not inst g2t_bbc_hit\n");
  }

 return kStOK;
}
Int_t StBbcSimulationMaker::Finish()
{
#ifdef BbcSimQa
  QaFile->Write();
  QaFile->Close();
#endif
  return StMaker::Finish();
}