/*
PbWO4 : 0.7% 
Plastic : ~40% Factor 50, but loos in fibre, only 7% ->3.5%  (overall 5 times)->10 times

 */
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TTree.h"
#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TObject.h"
#include "TPostScript.h"
#include <TRandom3.h>
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Matrix/Matrix.h"

#define SMEARING
using namespace std;
using namespace CLHEP;

const int nhclwedge = 36;
const int nhcalLayer = 7; // 17;
const int nhcalEtaDiv = 34;


double thetval[34]={0.652481, 0.684299, 0.718216, 0.755979, 0.791376, 0.834511, 0.879688, 0.930175, 0.983599, 1.04162, 1.10367, 1.17053, 1.24089, 1.31505, 1.3919, 1.47074, 1.54993, 1.6316, 1.71117, 1.78886, 1.86516, 1.93752, 2.00573, 2.06932, 2.12961, 2.18582, 2.23783, 2.28572, 2.32959, 2.366, 2.40572, 2.44111, 2.47378, 2.50406};

const double pival = acos(-1.0);
const double twopi = 2*pival;
const double pibytwo = pival/2.;
const double pibyfour = pival/4.;
const double degtorad= pival/180.;
const double radtodeg= 180./pival;
static unsigned int mypow_2[32];

double sfdph(double phi1, double phi2) {
  double dphi = phi2 - phi1;
  if (dphi > M_PI) { dphi = 2*M_PI - dphi;}
  if (dphi < -M_PI){ dphi = 2*M_PI + dphi;}
  return dphi;
}

double dgap(Hep3Vector v1, Hep3Vector v2) {
  double dthe = v1.theta() - v2.theta();
  double dphi = sfdph(v1.phi(),  v2.phi());
  return sqrt (dthe*dthe + dphi*dphi);
}

TRandom3* gRandom3 = new TRandom3();

int main(int argc, char** argv) {

  if (argc<5) {
    cout <<"Give seedthreh, pedwidth and otherthrs in MeV unit"<<endl;
    return 0;
  }

  const double pedwidth = atof(argv[1])/1000.0; // ECAL Noise per crystal : 40 MeV threshold peak
  const double otherthrs = atof(argv[2])/1000.0; // Energy threshold for other crystal : 200 MeV
  const double hclwidth = atof(argv[3])/1000.0; // HCAL noise in each tower :  25 MeV 
  const double hclthrs = atof(argv[4])/1000.0; // threshold for Hcal energy : 100 MeV
  
  cout <<"threshol "<< argv[0]<<" "<<pedwidth<<" "<<otherthrs<<" "<<hclwidth<<" "<<hclthrs<<endl;

  for (int ij=0; ij<32; ij++) {
    mypow_2[ij] = (int)pow(float(2), ij);
  }
  
  gStyle->SetOptStat(111111);
  char datafile[100];

  unsigned   irun;                // Run number of these events
  unsigned   ievt;                //Event number
  unsigned   ngent;
  
  const unsigned int ngenmx=50;
  int   pidin[ngenmx]; 	  //PID of incident particle
  float momin[ngenmx]; 	  //Energy of incident particle
  float thein[ngenmx];	  //Initial polar angle of incident particle
  float phiin[ngenmx];     //Initial azimuthal angle of incident particle
  
  // For Simulation output of ECAL
  unsigned int nsimhtEC;
  static const unsigned int nsimhtmxEC=2000;
  unsigned int detidEC[nsimhtmxEC];
  float energyEC[nsimhtmxEC];
  float thetaEC[nsimhtmxEC];
  float phiEC[nsimhtmxEC]; 

  unsigned int nsimhtTk; // Looked only for early showering

  // For Simulation output of HCAL
  unsigned int nsimhtHL;
  static const unsigned int nsimhtmxHL=2000;
  unsigned int detidHL[nsimhtmxHL];
  unsigned int timeHL[nsimhtmxHL]; 
  
  char rootfiles[100];
  char outfile[100];
  char outfilx[100];
  ifstream file_db;

  sprintf(rootfiles, "test_pion_klong.log");
  
  int len = strlen(rootfiles);
  strncpy(outfilx, rootfiles, len-4);
  outfilx[len-4]='\0';
  sprintf (outfile,"%s_xwoecl_%s_%s_%s_%s.root",outfilx, argv[1], argv[2], argv[3], argv[4]);

  float totsimenr =0;
  float totdigienr =0;
  float hclsimenr =0;
  float hcldigienr =0;
  
  //  cout<<"1mean "<<endl;
  TFile* fileOut = new TFile(outfile, "recreate");  
  TTree* Tout = new TTree("T2", "test");;

  Tout->Branch("totsimenr",&totsimenr,"totsimenr/F");
  Tout->Branch("totdigienr",&totdigienr,"totdigienr/F");
  Tout->Branch("hclsimenr",&hclsimenr,"hclsimenr/F");
  Tout->Branch("hcldigienr",&hcldigienr,"hcldigienr/F");

  file_db.open(rootfiles);  
  while(!(file_db.eof())) {
    int nentrymx=-1;
    
    file_db >> datafile>>nentrymx;
    if (strstr(datafile,"#")) continue;
    
    cout <<"datafile = "<<datafile<<endl;
    
    if(file_db.eof()) break; 
    
    TFile* fileIn = new TFile(datafile, "read");
    TTree *Tin = (TTree*)fileIn->Get("T1");

    Tin->SetBranchAddress("irun",&irun);
    Tin->SetBranchAddress("ievt",&ievt);
    
    Tin->SetBranchAddress("ngent",&ngent);
    Tin->SetBranchAddress("pidin",pidin);
    Tin->SetBranchAddress("momin",momin);
    Tin->SetBranchAddress("thein",thein);
    Tin->SetBranchAddress("phiin",phiin);
    
    //Simulation output of ECAL
    Tin->SetBranchAddress("nsimhtEC", &nsimhtEC);
    Tin->SetBranchAddress("detidEC", detidEC);
    Tin->SetBranchAddress("energyEC", energyEC);
    Tin->SetBranchAddress("thetaEC", thetaEC);
    Tin->SetBranchAddress("phiEC", phiEC);

    //Simulation output of HCAL
    Tin->SetBranchAddress("nsimhtHL", &nsimhtHL);
    Tin->SetBranchAddress("detidHL", detidHL);
    Tin->SetBranchAddress("timeHL", timeHL);

    int nentries = Tin->GetEntries();
    cout <<"nentr "<< datafile<<" "<<nentries<<endl;;
    
    for (int iev=0; iev<min(nentries, nentrymx); iev++) {
      fileIn->cd();
      Tin->GetEntry(iev);
      if (nsimhtEC==0 && nsimhtHL==0) continue;
      
      totsimenr =0;
      totdigienr =0;
      hclsimenr =0;
      hcldigienr =0;
      
      //Smear simulated energy with the inclusion of different noises
      
      fileOut->cd();
      
    Tout->Fill();
    } //end of events loop
    
    fileIn->cd();
    delete Tin;
    delete fileIn;
    
  } //end of file loop
  fileOut->cd();
  fileOut->Write();
  
}