//**************************************************************************************
// multi-event run script of ROMEO
//
// Usage:
// root -b -q 'RomeoManyRun.C(....)'
//
// Arguments:
//
// nevts                ... number of events
// maxpe                ... input NPE per shot. 
// beamwidth            ... time width of a beam shot. By default
//                          time property of a beam is uniform within the beamwidth.
// weightHits           ... if True or 1, all PEs within a shot are merged into 1 MCHit with weight. 
// isAmasimMode         ... if True or 1 romeo treats input beam as PhotoElectrons 
//                          instead of Photons. IceSim default is True.
// isQuickMode          ... if True or 1 romeo skips photocathode simulation.
//                          IceSim default is True.
// doWaveformSaturation ... if true or 1, it calculates saturation at waveform.
//                          if false or 0, romeo activates CurrentBasedSaturationModule
//                          which uses average current within outputtimeresol
//                          as an input of saturation curve. (semi-current saturation mode)
// outputtimeresol      ... any hit within this timewindow will be merged into one
//                          charge object. It's also used to calculate average current 
//                          at semi-current saturation mode.
//                          Negative falue is replaced to 2ns. If you set 0, binning is skipped and 
//                          RMOPMTCurrentBasedSaturation is not activated.
// nbeamshots           ... number of beam shots. beam bunch is sequently shot every "beamwidth".
//
// Beaming examples : 
//
// A)  weightHits = kFALSE, maxpe=4, beamwidth=10ns, nbeamshots=4:
//           |||||||||||||||| weight 1 / hit
//           ----
//           10ns 
//           |______________| 4 beam shots
//
// B)  weightHits = kTRUE, maxpe=4, beamwidth=10ns, nbeamshots=4:
//           |   |   |   |    weight 4 / hit
//           ----
//           10ns 
//           |______________| 4 beam shots
//
//**************************************************************************************


#include <TRandom.h>
#include <iomanip.h>
#include <TMath.h>

//TFile foutp;

int RomeoManyRun(
         Int_t nevts = 100,   // number of events
         Int_t maxpe = 100,  // number of PE per beamshot per event
         Double_t beamwidth = 1e-8, // 10ns square beam
         Bool_t weightHit = kFALSE, // will you make weighted hit?
         Bool_t isAmasimMode=kFALSE, 
         Bool_t isQuickMode=kFALSE,
         Bool_t doWaveformSaturation=kFALSE,
         Double_t outputtimeresol = -1,   //-1: use default 
         Int_t nbeamshot = 1) // number of beam shot
{

  gROOT->Reset();

  gSystem->Load("libPhysics.so");
  gSystem->Load("libHist.so");
  gSystem->Load("libGraf3d.so");
  gSystem->Load("../lib/libRMONumericalRecipes.so");
  gSystem->Load("../lib/libRMOKernel.so");
  gSystem->Load("../lib/libRMOOMAnalog.so");
  gSystem->Load("../lib/libRMOGlassGel.so");
  gSystem->Load("../lib/libRMOPMT.so");
    
  //foutp = new TFile(outname,"RECREATE");

  gStyle->SetOptStat(1111);
  gStyle->SetPalette(1);
  gStyle->SetCanvasColor(33);
  gStyle->SetFrameFillColor(18);

  RMOParameters::SetVerboseLevel(0);
  
  //======================================================================
  // 1) define data dir
  // you must chose one data directory which includes your calibration data

  RMOParameters::SetDataDir("../../resources/data/");

  //======================================================================
  // 2) create modules
  // create your modules. 
  
  // * glass+gel
  RMOOMAnalog                 *omanalog  = new RMOOMAnalog("StandardDOMIce3");  
  RMOPhotonGenerator          *phgen     = new RMOPhotonGenerator("PhotonGeneratorIce3");
  RMOWaveLengthSimulator      *wlsim     = new RMOWaveLengthSimulator("WaveLengthSimIce3");
  RMOIsotropicAngleSimulator  *angsim    = new RMOIsotropicAngleSimulator("IsotropicAngleSimIce3");
  RMOQuickPhotonPropagator    *phprop    = new RMOQuickPhotonPropagator("QuickPropagatorIce3");
  // * photo-cathode
  RMOPMTPhotoCathode          *cathode   = new RMOPMTPhotoCathode("PhotoCathodeSimIce3");
  // old constructor
  // RMOPMTPhotoCathode          *cathode   = new RMOPMTPhotoCathode("TA1062.table",
  //                                                                  "qe_representative.data",
  //                                                                   "qe_representative.data",
  //                                                                   "PhotoCathodeIce3TA1062","");
  // this is alternative PE generator. 
  // If you use the generator, you should skip all modules listed above.
  RMOPMTQuickPEGenerator      *pegenerator  = new RMOPMTQuickPEGenerator("QuickPEGeneratorIce3");

  // * charge response and gain
  RMOPMTChargeResponse        *qresponse = new RMOPMTChargeResponse("ChargeResIce3");
  // there are two saturation sub-modules. 
  //RMOPMTChargeBasedSaturation *saturation = new RMOPMTChargeBasedSaturation("ChargeBasedSaturationIce3"); 
  RMOPMTCurrentBasedSaturation *saturation = new RMOPMTCurrentBasedSaturation("CurrentBasedSaturationIce3"); 

  // * waveform
  RMOPMTWaveForm              *waveform  = new RMOPMTWaveForm("WaveFormIce3");

  //======================================================================
  // 3) set module parameters.
  // all configureable parameters should be defined here

  //--------------------------
  // set verbose...
  // 0 : no output
  // 1 : default
  // 2 : function called
  // 3 : precise
  // 4 : debug
  // 5 : trace

  /*
  omanalog->SetVerboseLevel(2);
  phgen->SetVerboseLevel(2);
  wlsim->SetVerboseLevel(2);
  angsim->SetVerboseLevel(2);
  phprop->SetVerboseLevel(2);
  cathode->SetVerboseLevel(2);
  pegenerator->SetVerboseLevel(2);
  qresponse->SetVerboseLevel(2);
  saturation->SetVerboseLevel(2);
  waveform->SetVerboseLevel(2);
  */

  omanalog->SetVerboseLevel(0);
  phgen->SetVerboseLevel(0);
  wlsim->SetVerboseLevel(0);
  angsim->SetVerboseLevel(0);
  phprop->SetVerboseLevel(0);
  cathode->SetVerboseLevel(0);
  pegenerator->SetVerboseLevel(0);
  qresponse->SetVerboseLevel(0);
  saturation->SetVerboseLevel(0);
  waveform->SetVerboseLevel(0);
 
  //--------------------------
  // for glass+gel
  // choose proper wavelength acceptance with current mode.
  // Note that you must use "withoutQE" wavelength-acceptance
  // if you kill the AmasimMode.

  RMOWaveLengthAcceptance::SetInputFileName("wavelength_acceptance_withoutQE_13inch_0deg.data"); 
  //RMOWaveLengthAcceptance::SetInputFileName("wavelength_acceptance_withoutQE_8mm.data");

  RMOAngleAcceptance::SetInputFileName("angle_acceptance_withoutQE_13inch_420nm.data");
  //RMOAngleAcceptance::SetInputFileName("angle_acceptance_withoutQE_8mm_385nm.data");

  RMOQuickPhotonPropagator::SetCalibFileNames("QuickPhotonPropagatorMap_13inch_420nm.data");
  //RMOQuickPhotonPropagator::SetCalibFileNames("QuickPhotonPropagatorMap_8mm_385nm.data");

  //----------------------------
  // for photo-cathode
  // if you want to change QECE manualy set here. default is TA0003
  //cathode->Set2DCEFactorAt337nm(0.228897); // CE of TA1062 based on the Chiba calibration

  //----------------------------
  // for charge-response
  // if you want to change spe parameters manualy set here. default is TA0003
  //qresponse->SetSQ0_ov_Q0(3.25837e-01);    // Charge Response of TA1062
  //qresponse->SetQt_ov_Q0(1.609385e-01);
  //qresponse->SetPe(2.228744e-01);
  //qresponse->SetGain(5.8e4);  // if icetop low gain

  // Output time resolution:  
  // its not PMT time resolution but a kind of binning window in order to 
  // generate one charge datum from several PEs within same time window
  // default is 2ns
  // If you set 0, it doesn't any binning and
  // you must switch off RMOCurrentBasedSaturation module.
  if (outputtimeresol < 0) {
     cerr << "you put invalid (negative) output time resolution. "
          << "2ns will be set instead." << endl;
     outputtimeresol = 2.e-9;
  }
  RMOParameters::SetOutputTimeResolution(outputtimeresol);

  // PMT time resolution:  
  // Default is 4ns. However, recent version of hit-constructor (or hit-maker)
  // convolutes this value: we need to set 0 if you use romeo with icesim.
  // For calibration experiment, just leave it as default.

  RMOParameters::SetPMTTimeResolution(0); 

  //======================================================================
  // 4) connect all modules
  // if you run the script with argument "true", it skips all simulation 
  // of glass+gel and Photo-Cathode simulation.
  
  if (isQuickMode) {
     // add pe generator only
     omanalog->Add(pegenerator);
     
  } else {
     // simulate photo-cathode
     omanalog->Add(phgen);
     omanalog->Add(wlsim);
     //omanalog->Add(angsim);
     omanalog->Add(phprop);
     omanalog->Add(cathode);
  }
  omanalog->Add(qresponse);

  if (!doWaveformSaturation) {

     if (outputtimeresol > 0) {
        // add saturation module only when you have positive 
        // outputtimeresolution
        qresponse->Add(saturation);
     } 
     // switch off waveform saturation (default)
     waveform->DoSaturationAtWaveform(kFALSE);

  } else {

     // set saturation flag on
     waveform->DoSaturationAtWaveform(kTRUE);
  }

  omanalog->Add(waveform);

  //======================================================================
  // 5) Set amasim mode or kill it.
  //------------------------------------------------------------------
  // Set mode flag AFTER added all modules. Default mode is AmasimMode.
  // AmasimMode :
  //   All incoming photons must reach to the photocathode.
  //   The photon absorption at glass+gel and incoming angle efficiency are
  //   already took into account the number of arrival photons.
  //   ROMEO assigns wavelength and incoming angle for all photons
  //   following analytical model of wavelength spectrum and angle acceptance
  //   without killing any photon.
  // Other :
  //   All incoming photons are regarded photons at DOM surface.
  //   ROMEO does fast simulation of photon absorption and angle acception
  //   inside the glass_gel modules then some photons are killed by following
  //   wavelength acceptance and angle acceptance.
  //

  RMOAttModeFlag::SetAmasimModeAll(kFALSE);   // Photons are killed inside the ROMEO.

  //======================================================================
  // 6) Simulation start.
  // Each line contains "ExecuteTask" calls specified function for 
  // all module added to omanalog
  //------------------------------------------------------------------

  TH1D *hq = new TH1D("hq","npe no saturation",10000, 0, 10000);
  TH1D *hqw = new TH1D("hqw","npe at waveform",5000, 0, 5000);

  //----------------------------
  // Configure
  // this function should be called beforhand the Process call
  // and should be called only once per one run.

  omanalog->ExecuteTask("Configure");

  //----------------------------
  // Process
  // this is the main part of simulation. 
  // For multiple run, you must be careful about
  // following points:
  //  - Call SetOwner() for input buffer
  //  - Call ExecuteTask("ClearEvent") in the end of a run
  // otherwise you well see enormous memory leak!

  ofstream out("SumCharge.dat");

  Double_t t0 = 100.e-9;

  for(Int_t i_evt = 0; i_evt < nevts; i_evt++){

     // create input buffer
     Int_t ipe = 0;
     RMOMCPhotonHitBuf hitarray("MCPhotonTimeHit");
     hitarray.SetOwner();  // set owner to allow to delete stored components


     for (Int_t itbin=0; itbin<nbeamshot; itbin++) {

        if (weightHit) {

            // make only one hit with weight of maxpe
            Double_t time = t0 + itbin*beamwidth;
            Double_t dist = 50;      // 50m
            TVector3 dir(0, 0, -1);

            RMOMCPhotonHit *hit = new RMOMCPhotonHit();
            hit->SetTime(time);
            hit->SetDirection(dir);
            hit->SetDistanceFromSource(dist);
            hit->SetSourceType(10);  // 10 is SPE pulse
            hit->SetWeight(maxpe);
            hitarray.Add(hit);

        } else {
         
           // make 1hit per 1pe
           for (ipe=0; ipe<maxpe; ipe++) {

              Double_t time = t0 + itbin*beamwidth 
                              + gRandom->Uniform(0, beamwidth);
              Double_t dist = 50;      // 50m
              TVector3 dir(0, 0, -1);

              RMOMCPhotonHit *hit = new RMOMCPhotonHit();
              hit->SetTime(time);
              hit->SetDirection(dir);
              hit->SetDistanceFromSource(dist);
              hit->SetSourceType(10);  // 10 is SPE pulse
              hitarray.Add(hit);
           }
        }
     }

     // set input buffer to first module
     omanalog->SetBuffer(&hitarray);

     //****************
     // call Process. 
     //****************

     omanalog->ExecuteTask("Process");

     //****************
     // simulation has done!
     //****************

     // output information

     Double_t qtot = qresponse->GetTotalNpe();
     
     hq->Fill(qtot);

     RMOPMTWaveFormBuf &wfbuf = (RMOPMTWaveFormBuf&)waveform->GetBuffer();
     Double_t qwtot = wfbuf.Integral(wfbuf.GetFirstTime(), wfbuf.GetLastTime()) / 43.;
     qwtot /= -(1e7*1.60217733e-19); // charge to npe
     hqw->Fill(qwtot);

     // output sum charge of the event
     qresponse->OutputSumChargeOfEvent(out);

     // delete all event buffer
     omanalog->ExecuteTask("ClearEvent");    

     if(i_evt%100 == 0){
        cerr << "event = " << i_evt << endl;
     }
  }

  //----------------------------
  // Finish
  // this function should be called after the Process call
  // and should be called only once per one run.

  omanalog->ExecuteTask("Finish");
 
  TCanvas *c1 = new TCanvas("c1","c1", 10, 10, 800, 400);
  c1->Divide(2,1);

  c1->cd(1);
  TF1 *g1 = new TF1("g1","gaus",0, maxpe);
  hq->Fit(g1,"Q");
  Double_t qmean = g1->GetParameter(1);
  Double_t qsig = g1->GetParameter(2) / TMath::Sqrt(nevts);

  c1->cd(2);
  TF1 *g2 = new TF1("g2","gaus",0, maxpe);
  hqw->Fit(g2,"Q");
  Double_t qwmean = g2->GetParameter(1);
  Double_t qwsig = g2->GetParameter(2) / TMath::Sqrt(nevts);

  cout << "CHARGES:\t\t" << maxpe << "\t\t" 
                         << qmean << "\t\t" << qsig << "\t\t"
                         << qwmean << "\t\t" << qwsig << endl;

  // close output stream
  out.close();
  out.clear();

}