1  /* ======================================================================== *\


2  !


3  ! *


4  ! * This file is part of CheObs, the Modular Analysis and Reconstruction


5  ! * Software. It is distributed to you in the hope that it can be a useful


6  ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.


7  ! * It is distributed WITHOUT ANY WARRANTY.


8  ! *


9  ! * Permission to use, copy, modify and distribute this software and its


10  ! * documentation for any purpose is hereby granted without fee,


11  ! * provided that the above copyright notice appears in all copies and


12  ! * that both that copyright notice and this permission notice appear


13  ! * in supporting documentation. It is provided "as is" without express


14  ! * or implied warranty.


15  ! *


16  !


17  !


18  ! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uniwuerzburg.de>


19  !


20  ! Copyright: CheObs Software Development, 20002009


21  !


22  !


23  \* ======================================================================== */


24 


25  //////////////////////////////////////////////////////////////////////////////


26  //


27  // MSimRandomPhotons


28  //


29  // Simulate poissonian photons. Since the distribution of the arrival time


30  // differences of these photons is an exonential we can simulate them


31  // using exponentially distributed time differences between two consecutive


32  // photons.


33  //


34  // FIXME: We should add the wavelength distribution.


35  //


36  // The artificial night sky background rate is calculated as follows:


37  //


38  // * The photon detection efficiency vs. wavelength of the detector is obtained


39  // from "PhotonDetectionEfficiency" of type "MParSpline"


40  //


41  // * The angular acceptance of the light collectors is obtained


42  // from "ConesAngularAcceptance" of type "MParSpline"


43  //


44  // * The spectral acceptance of the light collectors is obtained


45  // from "ConesTransmission" of type "MParSpline"


46  //


47  // * The reflectivity of the mirrors vs wavelength is obtained


48  // from "MirrorReflectivity" of type "MParSpline"


49  //


50  // The rate is then calculated as


51  //


52  // R = R0 * Ai * f


53  //


54  // R0 is the night sky background rate as given in Eckart's paper (divided


55  // by the wavelength window). Ai the area of the cones acceptance window,


56  // f is given as:


57  //


58  // f = nm * Min(Ar, sr*d^2)


59  //


60  // with


61  //


62  // nm being the integral of the product of the mirror reflectivity, the cone


63  // transmission and the photon detection efficiency.


64  //


65  // d the distance of the focal plane to the mirror


66  //


67  // Ar is the total reflective area of the reflector


68  //


69  // sr is the effective solid angle corresponding to the integral of the


70  // cones angular acceptance


71  //


72  // Alternatively, the nightsky background rate can be calculated from


73  // a spectrum as given in Fig. 1 (but versus Nanometers) in


74  //


75  // Chris R. Benn & Sara L. Ellison La Palma NightSky Brightness


76  //


77  // After proper conversion of the units, the rate of the pixel 0


78  // is then calculated by


79  //


80  // rate = f * nsb


81  //


82  // With nsb


83  //


84  // nsb = Integral(nsb spectrum * combines efficiencies)


85  //


86  // and f can be either


87  //


88  // Eff. angular acceptance Cones (e.g. 20deg) * ConeArea (mm^2)


89  // f = sr * A0


90  //


91  // or


92  //


93  // MirrorArea * Field of view of cones (deg^2)


94  // f = Ar * A0;


95  //


96  //


97  // Input Containers:


98  // fNameGeomCam [MGeomCam]


99  // MPhotonEvent


100  // MPhotonStatistics


101  // MCorsikaEvtHeader


102  // [MCorsikaRunHeader]


103  //


104  // Output Containers:


105  // MPhotonEvent


106  // AccidentalPhotonRate [MPedestalCam]


107  //


108  //////////////////////////////////////////////////////////////////////////////


109  #include "MSimRandomPhotons.h"


110 


111  #include <TRandom.h>


112 


113  #include "MMath.h" // RndmExp


114 


115  #include "MLog.h"


116  #include "MLogManip.h"


117 


118  #include "MParList.h"


119 


120  #include "MGeomCam.h"


121  #include "MGeom.h"


122 


123  #include "MPhotonEvent.h"


124  #include "MPhotonData.h"


125 


126  #include "MPedestalCam.h"


127  #include "MPedestalPix.h"


128 


129  #include "MCorsikaRunHeader.h"


130 


131  #include "MSpline3.h"


132  #include "MParSpline.h"


133  #include "MReflector.h"


134 


135  ClassImp(MSimRandomPhotons);


136 


137  using namespace std;


138 


139  // 


140  //


141  // Default Constructor.


142  //


143  MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)


144  : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),


145  fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"),


146  fFileNameNSB("resmc/nightskylapalma.txt")


147  {


148  fName = name ? name : "MSimRandomPhotons";


149  fTitle = title ? title : "Simulate possonian photons (like NSB or dark current)";


150  }


151 


152  // 


153  //


154  // Check for the necessary containers


155  //


156  Int_t MSimRandomPhotons::PreProcess(MParList *pList)


157  {


158  fGeom = (MGeomCam*)pList>FindObject(fNameGeomCam, "MGeomCam");


159  if (!fGeom)


160  {


161  *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;


162 


163  fGeom = (MGeomCam*)pList>FindObject("MGeomCam");


164  if (!fGeom)


165  {


166  *fLog << err << "MGeomCam not found... aborting." << endl;


167  return kFALSE;


168  }


169  }


170 


171  fEvt = (MPhotonEvent*)pList>FindObject("MPhotonEvent");


172  if (!fEvt)


173  {


174  *fLog << err << "MPhotonEvent not found... aborting." << endl;


175  return kFALSE;


176  }


177 


178  fStat = (MPhotonStatistics*)pList>FindObject("MPhotonStatistics");


179  if (!fStat)


180  {


181  *fLog << err << "MPhotonStatistics not found... aborting." << endl;


182  return kFALSE;


183  }


184 


185  fRates = (MPedestalCam*)pList>FindCreateObj("MPedestalCam", "AccidentalPhotonRates");


186  if (!fRates)


187  return kFALSE;


188 


189  /*


190  fEvtHeader = (MCorsikaEvtHeader*)pList>FindObject("MCorsikaEvtHeader");


191  if (!fEvtHeader)


192  {


193  *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;


194  return kFALSE;


195  }*/


196 


197  fRunHeader = (MCorsikaRunHeader*)pList>FindObject("MCorsikaRunHeader");


198  if (fSimulateWavelength && !fRunHeader)


199  {


200  *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;


201  return kFALSE;


202  }


203 


204  MReflector *r = (MReflector*)pList>FindObject("Reflector", "MReflector");


205  if (!r)


206  {


207  *fLog << err << "Reflector [MReflector] not found... aborting." << endl;


208  return kFALSE;


209  }


210 


211  const MParSpline *s1 = (MParSpline*)pList>FindObject("PhotonDetectionEfficiency", "MParSpline");


212  const MParSpline *s2 = (MParSpline*)pList>FindObject("ConesTransmission", "MParSpline");


213  const MParSpline *s3 = (MParSpline*)pList>FindObject("MirrorReflectivity", "MParSpline");


214  const MParSpline *s4 = (MParSpline*)pList>FindObject("ConesAngularAcceptance", "MParSpline");


215 


216  // Ensure that all efficiencies are at least defined in the raneg of the


217  // photon detection efficiency. We assume that this is the limiting factor


218  // and has to be zero at both ends.


219  if (s2>GetXmin()>s1>GetXmin())


220  {


221  *fLog << err << "ERROR  ConeTransmission range must be defined down to " << s1>GetXmin() << "nm (PhotonDetectionEffciency)." << endl;


222  return kFALSE;


223  }


224  if (s2>GetXmax()<s1>GetXmax())


225  {


226  *fLog << err << "ERROR  ConeTransmission range must be defined up to " << s1>GetXmax() << "nm (PhotonDetectionEffciency)." << endl;


227  return kFALSE;


228  }


229  if (s3>GetXmin()>s1>GetXmin())


230  {


231  *fLog << err << "ERROR  MirrorReflectivity range must be defined down to " << s1>GetXmin() << "nm (PhotonDetectionEffciency)." << endl;


232  return kFALSE;


233  }


234  if (s3>GetXmax()<s1>GetXmax())


235  {


236  *fLog << err << "ERROR  MirrorReflectivity range must be defined up to " << s1>GetXmax() << "nm (PhotonDetectionEffciency)." << endl;


237  return kFALSE;


238  }


239 


240  // If the simulated wavelength range exists and is smaller, reduce the


241  // range to it. Later it is checked that at both edges the transmission


242  // is 0. This must be true in both cases: The simulated wavelength range


243  // exceed the PDE or the PDE range exceeds the simulated waveband.


244  const Float_t wmin = fRunHeader && fRunHeader>GetWavelengthMin()>s1>GetXmin() ? fRunHeader>GetWavelengthMin() : s1>GetXmin();


245  const Float_t wmax = fRunHeader && fRunHeader>GetWavelengthMax()<s1>GetXmax() ? fRunHeader>GetWavelengthMax() : s1>GetXmax();


246 


247  const Int_t min = TMath::FloorNint(wmin);


248  const Int_t max = TMath::CeilNint(wmax);


249 


250  // Multiply all relevant efficiencies to get the total transmission


251  MParSpline eff;


252  eff.SetFunction("1", maxmin, min, max);


253 


254  eff.Multiply(*s1>GetSpline());


255  eff.Multiply(*s2>GetSpline());


256  eff.Multiply(*s3>GetSpline());


257 


258  // Effectively transmitted wavelength band in the simulated range


259  const Double_t nm = eff.GetSpline()>Integral();


260 


261  // Angular acceptance of the cones


262  const Double_t sr = s4 && s4>GetSpline() ? s4>GetSpline()>IntegralSolidAngle() : 1;


263 


264  {


265  const Double_t d2 = fGeom>GetCameraDist()*fGeom>GetCameraDist();


266  const Double_t conv = fGeom>GetConvMm2Deg()*TMath::DegToRad();


267  const Double_t f1 = TMath::Min(r>GetA()/1e4, sr*d2) * conv*conv;


268 


269  // Rate in GHz / mm^2


270  fScale = fFreqNSB * nm * f1; // [GHz/mm^2] efficiency * m^2 *rad^2 *mm^2


271 


272  const Double_t freq0 = fScale*(*fGeom)[0].GetA()*1000;


273 


274  *fLog << inf << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", freq0) << "MHz" << endl;


275 


276  // FIXME: Scale with the number of pixels


277  if (freq0>1000)


278  {


279  *fLog << err << "ERROR  Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;


280  return kFALSE;


281  }


282  }


283 


284  if (fFileNameNSB.IsNull())


285  return kTRUE;


286 


287  // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList>FindObject("MMcRunHeader");


288  // Set NumPheFromDNSB


289 


290  // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and


291  // nsb_mean 0.20


292  // Magic pixel: 0.00885361 deg


293  // dnsbpix = 0.2*50/15


294  // ampl = MMcFadcHeader>GetAmplitud()


295  // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)


296 


297  // Conversion of the yaxis


298  // 


299  // Double_t ff = 1; // myJy / arcsec^2 per nm


300  // ff *= 1e6; // Jy / arcsec^2 per nm


301  // ff *= 3600*3600; // Jy / deg^2


302  // ff *= 1./TMath::DegToRad()/TMath::DegToRad(); // Jy/sr = 1e26J/s/m^2/Hz/sr


303  // ff *= 1e26; // J/s/m^2/Hz/sr per nm


304 


305  const Double_t arcsec2rad = TMath::DegToRad()/3600.;


306  const Double_t f = 1e32 / (arcsec2rad*arcsec2rad);


307 


308  // Read night sky background flux from file


309  MParSpline flux;


310  if (!flux.ReadFile(fFileNameNSB))


311  return kFALSE;


312 


313  if (flux.GetXmin()>wmin)


314  {


315  *fLog << err << "ERROR  NSB flux from " << fFileNameNSB << " must be defined down to " << wmin << "nm." << endl;


316  return kFALSE;


317  }


318  if (flux.GetXmax()<wmax)


319  {


320  *fLog << err << "ERROR  NSB flux from " << fFileNameNSB << " must be defined up to " << wmax << "nm." << endl;


321  return kFALSE;


322  }


323 


324  MParSpline nsb;


325 


326  // Normalization to our units,


327  // conversion from energy flux to photon flux


328  nsb.SetFunction(Form("%.12e/(x*TMath::H())", f), maxmin, min, max);


329 


330  // multiply night sky background flux with normalization


331  nsb.Multiply(*flux.GetSpline());


332 


333  // Multiply with the total transmission


334  nsb.Multiply(*eff.GetSpline());


335 


336  // Check if the photon flux is zero at both ends of the NSB


337  if (eff.GetSpline()>Eval(min)>1e5)


338  {


339  *fLog << warn << "WARNING  Total transmission efficiency at ";


340  *fLog << min << "nm is not zero, but " << eff.GetSpline()>Eval(min) << "... abort." << endl;


341  }


342  if (eff.GetSpline()>Eval(max)>1e5)


343  {


344  *fLog << warn << "WARNING  Total transmission efficiency at ";


345  *fLog << max << "nm is not zero, but " << eff.GetSpline()>Eval(max) << "... abort." << endl;


346  }


347 


348  // Check if the photon flux is zero at both ends of the simulated region


349  if (eff.GetSpline()>Eval(wmin)>1e5)


350  {


351  *fLog << err << "ERROR  Total transmission efficiency at ";


352  *fLog << wmin << "nm is not zero... abort." << endl;


353  *fLog << " PhotonDetectionEfficency: " << s1>GetSpline()>Eval(wmin) << endl;


354  *fLog << " ConeTransmission: " << s2>GetSpline()>Eval(wmin) << endl;


355  *fLog << " MirrorReflectivity: " << s3>GetSpline()>Eval(wmin) << endl;


356  *fLog << " TotalEfficiency: " << eff.GetSpline()>Eval(wmin) << endl;


357  return kFALSE;


358  }


359  if (eff.GetSpline()>Eval(wmax)>1e5)


360  {


361  *fLog << err << "ERROR  Total transmission efficiency at ";


362  *fLog << wmax << "nm is not zero... abort." << endl;


363  *fLog << " PhotonDetectionEfficency: " << s1>GetSpline()>Eval(wmax) << endl;


364  *fLog << " ConeTransmission: " << s2>GetSpline()>Eval(wmax) << endl;


365  *fLog << " MirrorReflectivity: " << s3>GetSpline()>Eval(wmax) << endl;


366  *fLog << " TotalEfficiency: " << eff.GetSpline()>Eval(wmax) << endl;


367  return kFALSE;


368  }


369 


370  // Conversion from m to radians


371  const Double_t conv = fGeom>GetConvMm2Deg()*TMath::DegToRad()*1e3;


372 


373  // Angular acceptance of the cones


374  //const Double_t sr = s5.GetSpline()>IntegralSolidAngle(); // sr


375  // Absolute reflector area


376  const Double_t Ar = r>GetA()/1e4; // m^2


377  // Size of the cone's entrance window


378  const Double_t A0 = (*fGeom)[0].GetA()*1e6; // m^2


379 


380  // Rate * m^2 * Solid Angle


381  // 


382 


383  // Angular acceptance Cones (e.g. 20deg) * ConeArea


384  const Double_t f1 = A0 * sr; // m^2 sr


385 


386  // MirrorArea * Field of view of cones (e.g. 0.1deg)


387  const Double_t f2 = Ar * A0*conv*conv; // m^2 sr


388 


389  // FIXME: Calculate the reflectivity of the bottom by replacing


390  // MirrorReflectivity by bottom reflectivity and reflect


391  // and use it to reflect the difference between f1 and f2


392  // if any.


393 


394  // Total NSB rate in MHz per m^2 and sr


395  const Double_t rate = nsb.GetSpline()>Integral() * 1e6;


396 


397  *fLog << inf;


398 


399  // Resulting rates as if Razmick's constant had been used


400  // *fLog << 1.75e6/(600300) * f1 * eff.GetSpline()>Integral() << " MHz" << endl;


401  // *fLog << 1.75e6/(600300) * f2 * eff.GetSpline()>Integral() << " MHz" << endl;


402 


403  *fLog << "Conversion factor Fnu: " << f << endl;


404  *fLog << "Total reflective area: " << Form("%.2f", Ar) << " m" << UTF8::kSquare << endl;


405  *fLog << "Acceptance area of cone 0: " << Form("%.2f", A0*1e6) << " mm" << UTF8::kSquare << " = ";


406  *fLog << A0*conv*conv << " sr" << endl;


407  *fLog << "Cones angular acceptance: " << sr << " sr" << endl;


408  *fLog << "ConeArea*ConeSolidAngle (f1): " << f1 << " m^2 sr" << endl;


409  *fLog << "MirrorArea*ConeSkyAngle (f2): " << f2 << " m^2 sr" << endl;


410  *fLog << "Effective. transmission: " << Form("%.1f", nm) << " nm" << endl;


411  *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f1): " << Form("%.2f", rate * f1) << " MHz" << endl;


412  *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f2): " << Form("%.2f", rate * f2) << " MHz" << endl;


413  *fLog << "Using f2." << endl;


414 


415  // Scale the rate per mm^2 and to GHz


416  fScale = rate * f2 / (*fGeom)[0].GetA() / 1000;


417 


418  // FIXME: Scale with the number of pixels


419  if (rate*f2>1000)


420  {


421  *fLog << err << "ERROR  Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;


422  return kFALSE;


423  }


424 


425  return kTRUE;


426  }


427 


428  Bool_t MSimRandomPhotons::ReInit(MParList *pList)


429  {


430  // Overwrite the default set by MGeomApply


431  fRates>Init(*fGeom);


432  return kTRUE;


433  }


434 


435  // 


436  //


437  // Check for the necessary containers


438  //


439  Int_t MSimRandomPhotons::Process()


440  {


441  // Get array from event container


442  // const Int_t num = fEvt>GetNumPhotons();


443  //


444  // Do not produce pure pedestal events!


445  // if (num==0)


446  // return kTRUE;


447 


448  // Get array from event container


449  // FIXME: Use statistics container instead


450  const UInt_t npix = fGeom>GetNumPixels();


451 


452  // This is the possible window in which the triggered digitization


453  // may take place.


454  const Double_t start = fStat>GetTimeFirst();


455  const Double_t end = fStat>GetTimeLast();


456 


457  // Loop over all pixels


458  for (UInt_t idx=0; idx<npix; idx++)


459  {


460  // Scale the rate with the pixel size.


461  const Double_t rate = fFreqFixed + fScale*(*fGeom)[idx].GetA();


462 


463  (*fRates)[idx].SetPedestal(rate);


464 


465  // Calculate the average distance between two consequtive photons


466  const Double_t avglen = 1./rate;


467 


468  // Start producing photons at time "start"


469  Double_t t = start;


470  while (1)


471  {


472  // Get a random time for the photon.


473  // The differences are exponentially distributed.


474  t += MMath::RndmExp(avglen);


475 


476  // Check if we reached the end of the useful time window


477  if (t>end)


478  break;


479 


480  // Add a new photon


481  // FIXME: SLOW!


482  MPhotonData &ph = fEvt>Add();


483 


484  // Set source to NightSky, time to t and tag to pixel index


485  ph.SetPrimary(MMcEvtBasic::kNightSky);


486  ph.SetWeight();


487  ph.SetTime(t);


488  ph.SetTag(idx);


489 


490  // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant) FIXME: Reset?


491 


492  if (fSimulateWavelength)


493  {


494  const Float_t wmin = fRunHeader>GetWavelengthMin();


495  const Float_t wmax = fRunHeader>GetWavelengthMax();


496 


497  ph.SetWavelength(TMath::Nint(gRandom>Uniform(wmin, wmax)));


498  }


499  }


500  }


501 


502  // Resort the photons by time!


503  fEvt>Sort(kTRUE);


504 


505  // Update maximum index


506  fStat>SetMaxIndex(npix1);


507 


508  // Shrink


509  return kTRUE;


510  }


511 


512  // 


513  //


514  // Read the parameters from the resource file.


515  //


516  // FrequencyFixed: 0.040


517  // FrequencyNSB: 5.8


518  //


519  // The fixed frequency is given in units fitting the units of the time.


520  // Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz.


521  //


522  // The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore


523  // 0.040 would mean 40MHz/cm^2


524  //


525  Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print)


526  {


527  Bool_t rc = kFALSE;


528  if (IsEnvDefined(env, prefix, "FrequencyFixed", print))


529  {


530  rc = kTRUE;


531  fFreqFixed = GetEnvValue(env, prefix, "FrequencyFixed", fFreqFixed);


532  }


533 


534  if (IsEnvDefined(env, prefix, "FrequencyNSB", print))


535  {


536  rc = kTRUE;


537  fFreqNSB = GetEnvValue(env, prefix, "FrequencyNSB", fFreqNSB);


538  }


539 


540  if (IsEnvDefined(env, prefix, "FileNameNSB", print))


541  {


542  rc = kTRUE;


543  fFileNameNSB = GetEnvValue(env, prefix, "FileNameNSB", fFileNameNSB);


544  }


545 


546  if (IsEnvDefined(env, prefix, "SimulateCherenkovSpectrum", print))


547  {


548  rc = kTRUE;


549  fSimulateWavelength = GetEnvValue(env, prefix, "SimulateCherenkovSpectrum", fSimulateWavelength);


550  }


551 


552  return rc;


553  }

