#include <ppc/i3ppc.h>

#define ICETRAY
#include "ppc.cxx"

I3_MODULE(i3ppc);

i3ppc::i3ppc(const I3Context& ctx) : I3Module(ctx){
  AddOutBox("OutBox");
  verbose=false;
}

void i3ppc::Physics(I3FramePtr frame){
  log_trace("Entering i3ppc::Physics()");

  static bool ini=false;
  if (frame->GetStop() == I3Frame::Geometry || !ini){
    if(verbose) cout<<"Geometry frame"<<endl;
    I3GeometryConstPtr geo = frame->Get<I3GeometryConstPtr>();
    if(verbose){
      if(!geo) cout<<"    --> empty"<<endl;
      else cout<<"    --> contains "<<geo->omgeo.size()<<" entries"<<endl;
    }
    if(geo) if(geo->omgeo.size()>0){
      ppc::load(geo);
      ini=true;
    }
  }

  for (I3Frame::const_iterator iter=frame->begin(); iter!=frame->end(); iter++){
    I3MMCTrackListConstPtr mmclist = dynamic_pointer_cast<const I3MMCTrackList>(iter->second);
    if(mmclist){
      if(verbose) cout<<"  --> This is an I3MMCTrackList: "<<iter->first<<endl;
      for(I3MMCTrackList::const_iterator it=mmclist->begin(); it!=mmclist->end(); ++it){
	const I3MMCTrack& m=*it;
	const I3Particle& p=it->GetI3Particle();
	int minor = p.GetMinorID();
	unsigned long long major = p.GetMajorID();
	if(verbose) cout<<"      MMCTrack info added for particle with id: ("<<minor<<","<<major<<")"<<endl;
	i3mmctracks[make_pair(minor, major)]=&m;
      }
    }
  }

  ppc::hitnum=0;
  ppc::mchits->erase(ppc::mchits->begin(), ppc::mchits->end());
  bool trees=false, lists=false;
  for (I3Frame::const_iterator iter=frame->begin(); iter!=frame->end(); iter++){
    I3MCTreeConstPtr mctree = dynamic_pointer_cast<const I3MCTree>(iter->second);
    if(mctree){
      if(verbose) cout<<"  --> This is an I3MCTree: "<<iter->first<<endl;
      if(lists) continue; else trees=true;
    }
    else{
      I3MCListConstPtr mclist = dynamic_pointer_cast<const I3MCList>(iter->second);
      if(mclist){
	if(verbose) cout<<"  --> This is an I3MCList: "<<iter->first<<endl;
	if(trees) continue; else lists=true;
	mctree = I3MCTreeUtils::ListToTree(mclist);
	if(mctree) if(verbose) cout<<"  --> Converted to I3MCTree: "<<iter->first<<endl;
      }
    }
    if(mctree){
      int num=0;
      const vector<I3Particle>& primaries = I3MCTreeUtils::GetPrimaries(mctree);
      for(vector<I3Particle>::const_iterator it=primaries.begin(); it!=primaries.end(); ++it, num++){
	const I3Particle& p=*it;
	if(verbose) cout<<"extracting a particle "<<p.GetType()<<endl;
	pparticle(p, num, mctree);
      }
    }
  }

  i3mmctracks=map<pair<int,unsigned long long>,const I3MMCTrack *>();
  frame->Put("MCHitSeriesMap", ppc::mchits);

  log_debug("Pushing the frame");
  PushFrame(frame,"OutBox");
}

double i3ppc::res(double x, double mod){
  return x-mod*int(x/mod);
}

void i3ppc::pparticle(const I3Particle& p, int n, I3MCTreeConstPtr tree){
  int num=0;
  pair<int,unsigned long long> id = make_pair(p.GetMinorID(), p.GetMajorID());
  static int level=0;
  int type=p.GetType();
  if(verbose){
    for(int i=0; i<level; i++) cout<<"    ";
    cout<<"      --> particle # "<<n<<" of type "<<type<<" id: ("<<id.first<<","<<id.second<<")"<<endl;
  }

  const vector<I3Particle>& secondaries = I3MCTreeUtils::GetDaughters(tree, p);
  level++;
  for(vector<I3Particle>::const_iterator it=secondaries.begin(); it!=secondaries.end(); ++it, num++) pparticle(*it, num, tree);
  level--;

  if(type==5 || type==6){
    map<int,int> strings;
    {
      const double sol=0.299792458;
      double E0=p.GetEnergy()/I3Units::GeV;
      double t0=p.GetTime()/I3Units::ns;
      double ll=p.GetLength()/I3Units::m;

      double Ei, Ef, ti, tf;

      if(i3mmctracks.find(id)!=i3mmctracks.end()){
	const I3MMCTrack * i3mmctrack = i3mmctracks[id];
	ti=i3mmctrack->GetTi()/I3Units::ns;
	Ei=i3mmctrack->GetEi()/I3Units::GeV;
	tf=i3mmctrack->GetTf()/I3Units::ns;
	Ef=i3mmctrack->GetEf()/I3Units::GeV;
	if(Ei==0){
	  ti=t0;
	  Ei=E0;
	}
	if(Ef<0){
	  tf=t0+ll/sol;
	  Ef=0;
	}
	if(Ei>0){
	  if(verbose) cout<<"extracting a muon Ei="<<Ei<<" "<<secondaries.size()<<" secondaries"<<endl;

	  double th=180-p.GetZenith()/I3Units::degree;
	  double ph=res(180+p.GetAzimuth()/I3Units::degree, 360);

	  double sinth=sin(th*M_PI/180);
	  double costh=cos(th*M_PI/180);
	  double sinph=sin(ph*M_PI/180);
	  double cosph=cos(ph*M_PI/180);

	  double nx=sinth*cosph;
	  double ny=sinth*sinph;
	  double nz=costh; 

	  double totE=0;
	  for(vector<I3Particle>::const_iterator it=secondaries.begin(); it!=secondaries.end(); ++it){
	    double t=it->GetTime()/I3Units::ns;
	    if(t>=ti && t<=tf) totE+=it->GetEnergy()/I3Units::GeV;
	  }
	  double sl0=tf>ti?((Ef-Ei+totE)/(ti-tf)):0;
	  double slmax=(0.21+8.8e-3*log(Ei)/log(10.))*sol;  // for ice only at Ecut=500 MeV
	  double sl=min(sl0, slmax);

	  double x=p.GetX()/I3Units::m;
	  double y=p.GetY()/I3Units::m;
	  double z=p.GetZ()/I3Units::m+ppc::zoff;
	  double t=ti;
	  double l=sol*(tf-ti);
	  double E=Ei;
	  double Type=p.GetType();
	  double type=Type;
	  double new_t=tf;

	  for(vector<I3Particle>::const_iterator it=secondaries.begin(); it!=secondaries.end(); ++it){
	    new_t=it->GetTime()/I3Units::ns;
	    if(new_t>=ti && new_t<=tf){
	      l=sol*(new_t-t);

	      if(verbose) cout<<"muon at z="<<z<<" length="<<l<<" E="<<E<<" t="<<t<<endl;
	      if(Type==I3Particle::MuPlus || Type==I3Particle::MuMinus){
		ppc::mcpids=id;
		ppc::muon(x, y, z, nx, ny, nz, t, E, l);
	      }

	      x=it->GetX()/I3Units::m;
	      y=it->GetY()/I3Units::m;
	      z=it->GetZ()/I3Units::m+ppc::zoff;
	      double e=it->GetEnergy()/I3Units::GeV;

	      E-=e+sl*max(new_t-t, 0.);
	      t=new_t;
	      type=it->GetType();
	      ppc::mcpids=make_pair(it->GetMinorID(), it->GetMajorID());

	      if(type==I3Particle::DeltaE || type==I3Particle::Brems || type==I3Particle::PairProd ||
		 type==I3Particle::EPlus || type==I3Particle::EMinus){
		if(verbose) cout<<"e-m cascade at z="<<z<<" E="<<e<<" t="<<t<<endl;
		ppc::cascade(x, y, z, nx, ny, nz, t, e, 0);
	      }
	      else if(type==I3Particle::NuclInt || type==I3Particle::Hadrons){
		if(verbose) cout<<"hadr cascade at z="<<z<<" E="<<e<<" t="<<t<<endl;
		ppc::cascade(x, y, z, nx, ny, nz, t, e, 1);
	      }
	    }
	  }
	  l=sol*(tf-t);
	  if(verbose) cout<<"muon at z="<<z<<" length="<<l<<" E="<<E<<" t="<<t<<endl;
	  if(Type==I3Particle::MuPlus || Type==I3Particle::MuMinus){
	    ppc::mcpids=id;
	    ppc::muon(x, y, z, nx, ny, nz, t, E, l);
	  }
	}
      }
      else{
	ti=t0;
	Ei=E0;
	tf=t0+ll/sol;
	Ef=0;
      }
    }
  }
}