/////////////////////////////////////////////////////////////////////////////////
//
// Simple code to simulate the detector response of CDC
// Output root format is almost same as in Chen's stand-alone simulation
//
/////////////////////////////////////////////////////////////////////////////////
#include <memory>
#include <iostream>
#include <vector>

#include <cometEventLoop.hxx>

#include <TROOT.h>
#include <IG4Trajectory.hxx>
#include <IG4TrajectoryPoint.hxx>
#include <IG4PrimaryVertex.hxx>
//#include <IG4HitSegment.hxx>
#include <COMETHitSegment.hh>
#include <ICOMETEvent.hxx>
#include <IHitSelection.hxx>
#include <IMCHit.hxx>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <IHandle.hxx>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TVector2.h>
#include <TVector3.h>
#include <TArc.h>
#include <TF1.h>
#include <TGraph.h>
#include <TArrow.h>
#include <TRandom.h>
#include <TRandom3.h>

std::string ofileName;

using unit::mm;
using unit::cm;
using unit::eV;
using unit::keV;
using unit::ns;

namespace {
   // constants
   const double deg2rad        = TMath::DegToRad();
   const double rad2deg        = TMath::RadToDeg();
   const double two_pi         = TMath::TwoPi();
   const double half_pi        = TMath::Pi()/2.;
   // global geometry
   const double xoffset        = 641.00*cm; // cm
   const double zoffset        = 765.65*cm; // cm
   // gas simulation
   const double edep2ncls      = 1./(28.7*eV); // (#of clusters)/(energy deposit [eV]) from Garfield++
   // cdc geometry
   const double lenOffset      = 129.48*cm; // cm (Offset length of CDC)
   const double tiltAngle      = 10.*deg2rad;
   const double rminLayer      = 52.2*cm; // minimum radius of sense layer
   const double hCell          =  1.6*cm; // height of cell @EP
   const int    nLayer         =  18;
   const int    nWire0         = 396/2; // #of wires @ 0-th layer
   const int    nDeltaW        =  12/2; // step of the #of wires
   const int    innerskipHole  =  12/2; // #of skip holes @EP to get the stereo angle
   const int    outerskipHole  =  14/2; // #of skip holes @EP to get the stereo angle
   const double rinCDC         =  50.0*cm; // cm
   const double routCDC        =  81.5*cm; // cm
   const double driftV         = 0.025*mm/ns; // drift velocity
   const double tstart         =  200.*ns; // start time of default time window
   const double tend           = 1200.*ns; // end   time of default time window
   const double zmax           =  160.*cm;
}

// wire geometry function
double funcXYWire (double *x, double *par)
{
   // x -- x value
   // par[0] -- a_x
   // par[1] -- a_y
   // par[2] -- x_0
   // par[3] -- y_0
   double val = par[1]/par[0]*(x[0]-par[2]) + par[3];
   return val;
}

/// User code for the event loop.  This can be copied and used as a skeleton
/// for generic user loop programs.
class IUserLoop: public COMET::ICOMETEventLoopFunction {
public:
    /// Initialize any class specific variables, but most of the work can be
    /// done in Initialize.  Don't create histograms here!
    IUserLoop() {}
    virtual ~IUserLoop() {};

    /// Print a usage message.  This is generally called when there is a
    /// command line input error.
    void Usage(void) {
    }

    /// Set an option and return true if it is valid.  This is called by the
    /// event loop command line argument processing code for each "-O
    /// [name]=[value]" option found on the command line.  If the command line
    /// has "-O [name]" without a value, then the value string will be equal
    /// to "".  This must return false if the option was not correctly
    /// processed.
    virtual bool SetOption(std::string option,std::string value="") {
        return false;
    }

    /// Transform the position from global to local
    virtual bool transformPosition(TVector3 global, TVector3 &local, Bool_t trans = true) {
        local = global;
        if (trans) {
          local.SetXYZ(local.X()-xoffset, local.Y(), local.Z()-zoffset);
          local.RotateY(-half_pi);
        }
        return true;
    }

    /// Transform the momentum direction from global to local
    virtual bool transformMomentum(TVector3 global, TVector3 &local, Bool_t trans = true) {
        local = global;
        if (trans) {
          local = TVector3(-global.Z(), global.Y(), global.X());
        }
        return true;
    }

    /// Number of sense wire at i-th layer
    int GetNumberOfWirePerLayer(int ilayer) {
       return (nWire0+ilayer*nDeltaW);
    }

    /// Find the layer ID from hit position
    int FindLayer(TVector3 pos) { // return the number of sense wire (assuming almost square cell shape)
       int layerID = -1;
       int wirenum[nLayer];
       int rpos_wire[nLayer];
       double xy[2];
       for (int iLayer = 0; iLayer < nLayer; iLayer++) {
         if (iLayer>0) wirenum[iLayer] = wirenum[iLayer-1]+GetNumberOfWirePerLayer(iLayer-1);
         else wirenum[iLayer] = 0;
         GetWirePosition(wirenum[iLayer], xy, pos.Z());
         rpos_wire[iLayer] = hypot(xy[0], xy[1]);
       }
       for (int iLayer = 0; iLayer < nLayer; iLayer++) {
         if (pos.Perp()<rpos_wire[iLayer]) {
           layerID = iLayer;
           break;
         } else if (iLayer<nLayer-1 && fabs(pos.Perp()-rpos_wire[iLayer])<fabs(pos.Perp()-rpos_wire[iLayer+1])) {
           layerID = iLayer;
           break;
         } else if (iLayer==nLayer-1){// && fabs(pos.Perp()-rpos_wire[iLayer])<fabs(rpos_wire[iLayer]-rpos_wire[iLayer-1])/2.) {
           layerID = iLayer;
           break;
         }
       }
       return layerID;
    }

    /// Find the nearest wire from hit position
    int FindWire(int ilayer,TVector3 pos) {
       int wirenum = 0;
       for (int iLayer = 0; iLayer < nLayer; iLayer++) {
         if (iLayer>0) wirenum += GetNumberOfWirePerLayer(iLayer-1);
         if (iLayer==ilayer) break;
       }

       double xy0[2];
       GetWirePosition(wirenum, xy0, pos.Z());

       double phi_hit = pos.Phi();
       double phi0    = atan2(xy0[1], xy0[0]);
       double dphi    = two_pi/((double)GetNumberOfWirePerLayer(ilayer));
       double phi     = phi0;
       int wireID = -1;
       for (int iwire = wirenum; iwire < wirenum+GetNumberOfWirePerLayer(ilayer); iwire++) {
         if (fabs(phi_hit-phi)<dphi/2.) {
           wireID = iwire;
           break;
         } else if (fabs(phi_hit-phi-two_pi)<dphi/2.) {
           wireID = iwire;
           break;
         } else if (fabs(phi_hit-phi+two_pi)<dphi/2.) {
           wireID = iwire;
           break;
         }
         phi += dphi;
       }
       return wireID;
    }

    /// Calculate the wire position in XY from z position
    virtual void GetWirePosition (int iwire, double *xy, double zpos) {
       xy[0] = fXYWire[iwire]->GetParameter(0)*zpos + fXYWire[iwire]->GetParameter(2);
       xy[1] = fXYWire[iwire]->GetParameter(1)*zpos + fXYWire[iwire]->GetParameter(3);
    }

    COMET::IWritableMCHit *CreateMCHit(COMET::IChannelId channelId, std::vector < COMET::IG4VHit* > parents, double charge,double time){
       COMET::IWritableMCHit *mcHit = new COMET::IWritableMCHit();
       mcHit->SetCharge(charge);
       mcHit->SetTime(time);
       std::vector < COMET::IG4VHit* >::iterator parentIter;
       for (parentIter = parents.begin(); parentIter != parents.end(); ++parentIter) {
         if (*parentIter){
           mcHit->GetContributors().push_back( *parentIter );
         }
       }
       return mcHit;
    }

    int GetCellID(int iwire, int ilayer){
        int numberofwire = 0;
        int cellID = -1;
        for(int i=0;i<ilayer;i++){
            numberofwire += GetNumberOfWirePerLayer(i);
        }
        cellID = iwire - numberofwire;
        return cellID;
    }

    /// Called for each event inside the event loop, and returns true if the
    /// event should be saved to the output file.  If the remainder of the
    /// current file should be skipped, this should through the
    /// ENextEventLoopFile exception.  If the daughter class defines this
    /// operator, it should not define the Process() method.
    bool operator () (COMET::ICOMETEvent& event) {

      bool useOnlyFirstTurn = true;
      bool global2local     = true; // need to convert the coordinate...
      //bool global2local     = false; // need to convert the coordinate...

      //--- Added by hayashi
      //int EventNumber = event.GetEventId();
      //CdcCell_eventID = EventNumber;

      COMET::IHandle<COMET::IG4TrajectoryContainer> trajs=event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
      TLorentzVector momentum;
      TLorentzVector position;
      if (trajs) {
        //std::cout<<"Number of trajectories "<<trajs->size()<<std::endl;

        for (COMET::IG4TrajectoryContainer::const_iterator j=trajs->begin();j!=trajs->end(); j++) {
          //std::cout << " Primary ID "<< j->second.GetPrimaryId() << std::endl;
          if (j->second.GetPDGEncoding()==11 && j->second.GetTrackId()==1) {
            momentum=j->second.GetInitialMomentum();
            //position=j->second.GetInitialPosition();
            numSignal++;
            if (global2local) {
              momvtxx = -momentum.Z();
              momvtxy =  momentum.Y();
              momvtxz =  momentum.X();
            } else {
              momvtxx = momentum.X();
              momvtxy = momentum.Y();
              momvtxz = momentum.Z();
            }

            //posvtxx = position.X();
            //posvtxy = position.Y();
            //posvtxz = position.Z();
          }
        }
      }else{
        std::cout<<"Trajectories not added "<<std::endl;
        return false;
      }

      bool triggerHit = false;
      int  nhits_cth  = 0;
      int  nhits_cthI = 0;
      TVector3 posGlobal, posLocal;
      TVector3 momGlobal, momLocal;
      TVector2 mom2D, pos2D;
      COMET::IHandle<COMET::IG4HitContainer> cthg4Hits = event.Get<COMET::IG4HitContainer>("truth/g4Hits/CTH");
      if(cthg4Hits) {

        for (COMET::IG4HitContainer::const_iterator g4Hit = cthg4Hits->begin();
          g4Hit != cthg4Hits->end();++g4Hit) {

          COMET::IG4HitSegment* g4Seg = dynamic_cast<COMET::IG4HitSegment*>(*g4Hit);

          double edep = g4Seg->GetEnergyDeposit();

          if (!edep) continue;

          momGlobal = TVector3(g4Seg->GetMomentumX(), g4Seg->GetMomentumY(), g4Seg->GetMomentumZ());
          posGlobal = TVector3(g4Seg->GetStartX(),    g4Seg->GetStartY(),    g4Seg->GetStartZ());
          transformPosition(posGlobal, posLocal, global2local);
          transformMomentum(momGlobal, momLocal, global2local);

          CthHits_x.push_back(posLocal.X());
          CthHits_y.push_back(posLocal.Y());
          CthHits_z.push_back(posLocal.Z());
          CthHits_t.push_back(g4Seg->GetStartT());
          CthHits_px.push_back(momLocal.X());
          CthHits_py.push_back(momLocal.Y());
          CthHits_pz.push_back(momLocal.Z());
          CthHits_edep.push_back(edep);

          mom2D = momLocal.XYvector();
          pos2D = posLocal.XYvector();

          nhits_cth++;
          if (fabs(mom2D.DeltaPhi(pos2D))>=110.*deg2rad&&g4Seg->GetTrackId()==1) nhits_cthI++;
        }
        if (nhits_cth>0) triggerHit = true;
      }
      CthHits_nHits = nhits_cth;
      CthHits_nHitsI= nhits_cthI;
      //if (!triggerHit) std::cout << "no trigger hit" << std::endl;

      // Loop over CDC hits in the event and histogram position
      COMET::IHandle<COMET::IG4HitContainer> cdcg4Hits = event.Get<COMET::IG4HitContainer>("truth/g4Hits/CDC");

      if(cdcg4Hits) {

        int nhits = 0;
        // fill true hit information
        std::vector < TVector3 >  poshits(0); // position
        std::vector < TVector3 >  momhits(0); // momentum
        std::vector < double >   edephits(0); // energy deposit
        std::vector < double >    tofhits(0); // time of flight
        std::vector < int >      wirehits(0); // wire number
        std::vector < int >     layerhits(0); // layer number
        std::vector < int >       trkhits(0); // track id
        std::vector < int >      primhits(0); // primary id
        std::vector< std::vector < COMET::IG4VHit* > > parents(nWires);

        for (int iwire = 0; iwire < nWires; iwire++) parents[iwire].resize(0); // initialize

        double stepLen = 0;
        double edep    = 0;
        double time    = 0;
        //////// Merge true hits information
        for (COMET::IG4HitContainer::const_iterator g4Hit = cdcg4Hits->begin();
          g4Hit != cdcg4Hits->end();++g4Hit) {

          COMET::IG4HitSegment* g4Seg = dynamic_cast<COMET::IG4HitSegment*>(*g4Hit);

          edep = g4Seg->GetEnergyDeposit();

          if (!edep) continue;
          time = g4Seg->GetStartT();
          //std::cout << "time = " << time/ns << "[ns]" <<std::endl;
          TVector3 pre_global  = TVector3(g4Seg->GetStartX(), g4Seg->GetStartY(), g4Seg->GetStartZ());
          TVector3 post_global = TVector3(g4Seg->GetStopX(), g4Seg->GetStopY(), g4Seg->GetStopZ());
          const std::vector<TVector3>& path = g4Seg->GetHitSteps();
          int istep;
          for (istep = 0; istep < path.size(); istep++) {
            COMETLog("step position = " << path.at(istep).x() << "," << path.at(istep).y() << "," << path.at(istep).z());
            if (istep>0) {
              COMETLog("step length = " << sqrt( pow(path.at(istep).x()-path.at(istep-1).x(),2) +
                                                 pow(path.at(istep).y()-path.at(istep-1).y(),2) +
                                                 pow(path.at(istep).z()-path.at(istep-1).z(),2)) << " mm");
            }
            istep++;
          }

          momGlobal = TVector3(g4Seg->GetMomentumX(), g4Seg->GetMomentumY(), g4Seg->GetMomentumZ());
          transformMomentum(momGlobal, momLocal, global2local);

          //std::cout << "Pre position  : " << pre_global.X() << "  " << pre_global.Y() << "  " << pre_global.Z() << std::endl;
          //std::cout << "Post position : " << post_global.X() << "  " << post_global.Y() << "  " << post_global.Z() << std::endl;

          TVector3 pre_local, post_local, dir_local;
          transformPosition(pre_global,  pre_local,  global2local);
          transformPosition(post_global, post_local, global2local);
          /*CdcHits_nHits++;
          CdcHits_x.push_back(post_local.X());
          CdcHits_y.push_back(post_local.Y());
          CdcHits_z.push_back(post_local.Z());*/

          //std::cout << "Pre pos local  : " << pre_local.X() << "  " << pre_local.Y() << "  " << pre_local.Z() << std::endl;
          //std::cout << "Post pos local : " << post_local.X() << "  " << post_local.Y() << "  " << post_local.Z() << std::endl;
          //std::cout << "dE/dx = " << edep/keV << " [keV] / " << stepLen/cm << " [cm] " << std::endl;

          //if (fabs(momvtxz)>20.) std::cout << "pos = " << post_local.X() << " , " << post_local.Y() << " , " << post_local.Z() << std::endl;
          //int layer = FindLayer(post_local);
          //int wire  = FindWire(layer, post_local);

          ///////////////// generate clusters here !!!!!
          dir_local = post_local - pre_local;
          stepLen   = dir_local.Mag();
          dir_local = dir_local.Unit();
          int ncls = rnd->Poisson(edep * edep2ncls); // number of clusters
          double dedep = edep / ((double)ncls);
          double dstep = stepLen/((double)ncls);
          double dlen = 0.;
          TVector3 clsxyz;
          for (int icl = 0; icl < ncls; icl++) {
            dlen = rnd->Uniform(stepLen);
            // make cluster
            clsxyz = pre_local + dlen * dir_local;
            int layer = FindLayer(clsxyz);
            int wire  = FindWire(layer, clsxyz);

            if (layer<0||wire<0) {
              std::cout << "BAD WIRE OR LAYER POSITION ???" << std::endl;
              std::cout << layer << " " << wire << " " << post_local.X() << " " << post_local.Y() << " " << post_local.Z() << std::endl;
              std::cout << "  " << clsxyz.X() << " " << clsxyz.Y() << " " << clsxyz.Z() << std::endl;
            }

            poshits.push_back(clsxyz);
            momhits.push_back(momLocal);
            edephits.push_back(dedep);
            tofhits.push_back(time);
            wirehits.push_back(wire);
            layerhits.push_back(layer);
            trkhits.push_back(g4Seg->GetTrackId());
            primhits.push_back(g4Seg->GetPrimaryId());
            parents[wire].push_back(g4Seg);
            nhits++;
          }
          //std::cout << "wirehit = " << wire << std::endl;
        }
        std::vector < TVector3 > poswire(nWires);   // closest point
        std::vector < TVector3 > momwire(nWires);   // momentum @ closest point
        std::vector < double >   edepwire(nWires);  // total energy deposit
        std::vector < double >   timewire(nWires);  // drift time + tof
        std::vector < double >   ddistwire(nWires); // drift distance (minimum)
        std::vector < int >      nclusters(nWires); // number of ionized electrons
        std::vector < bool >     hitwire(nWires,false); //added by hayashi
        std::vector < double >   hitcell(nWires);   //added by hayashi

        int cdcCell_nHits = 0;
        std::vector < double > cdcCell_edep(0);
        std::vector < double > cdcCell_x(0);
        std::vector < double > cdcCell_y(0);
        std::vector < double > cdcCell_z(0);
        std::vector < double > cdcCell_wx(0);
        std::vector < double > cdcCell_wy(0);
        std::vector < double > cdcCell_wz(0);
        std::vector < double > cdcCell_t(0);
        std::vector < double > cdcCell_driftD(0);
        std::vector < double > cdcCell_driftDtrue(0);
        std::vector < int >    cdcCell_cellID(0);
        std::vector < int >    cdcCell_layerID(0);
        std::vector < double > cdcCell_px(0);
        std::vector < double > cdcCell_py(0);
        std::vector < double > cdcCell_pz(0);
        std::vector < int >    cdcCell_hittype(0);
        std::vector < int >    cdcCell_nclusters(0);
        std::vector < double > cdcCell_absmom(0);
        for (int iwire = 0; iwire < nWires; iwire++) {
          nclusters[iwire] = 0;
          if (!parents[iwire].size()) continue;
          hitwire[iwire] = true;
          double edepTot = 0;
          double tof     = 0;
          double ddist   = 1e9;
          double ddistsm = 1e9;
          double arriveT = 1e9;
          double tmpArrT = 1e9;
          double xy[2] = {0,0};
          int layerID = -1;
          int wireID  = -1;

          bool signalHits = false;
          for (int ihit = 0; ihit < nhits; ihit++) {
            if (wirehits[ihit]!=iwire) continue;
            edepTot += edephits[ihit];
            tof      = tofhits[ihit];
            //if (trkhits[ihit]==1) {
            if (ihit%20==0) {
              CdcHits_nHits++;
              CdcHits_x.push_back(poshits[ihit].X());
              CdcHits_y.push_back(poshits[ihit].Y());
              CdcHits_z.push_back(poshits[ihit].Z());
            }
            nclusters[iwire]++;
            GetWirePosition(iwire, xy, poshits[ihit].Z());
            //std::cout << "drifttime = " << hypot(xy[0]-poshits[ihit].X(), xy[1]-poshits[ihit].Y())/driftV << ", tof = " << tof/ns << " [ns]" << std::endl;
            tmpArrT = hypot(xy[0]-poshits[ihit].X(), xy[1]-poshits[ihit].Y())/driftV + tof;

            if (hypot(xy[0]-poshits[ihit].X(), xy[1]-poshits[ihit].Y())<ddist) {
              ddist   = hypot(xy[0]-poshits[ihit].X(), xy[1]-poshits[ihit].Y()); // closest point
              ddistsm = ddist + rnd->Gaus(0,resolution);
              arriveT = tmpArrT;
              poswire[iwire] = poshits[ihit];
              //std::cout << "track ID = " << trkhits[ihit] << std::endl;
              if (trkhits[ihit]==1 && momhits[ihit].Mag()>60.) {
                momwire[iwire] = momhits[ihit];
                signalHits     = true;
              }
            }
            layerID = layerhits[ihit];
            wireID  = wirehits[ihit];
          }
          if (wireID<0 || layerID<0) continue;

          edepwire[iwire] = edepTot;
          //std::cout << iwire << ", drift distance = " << ddist/mm << " [mm], arrival time = " << arriveT/ns
          //          << " [ns], energy deposit = " << edepTot/keV << " [keV]"  << std::endl;
          ddistwire[iwire] = ddist;
          timewire[iwire]  = arriveT;
          HitPosHist->Fill(poswire[iwire].X(), poswire[iwire].Y());
          EneDepoHist->Fill(edepTot/keV);
          DDistHist->Fill(ddist/mm);

          COMET::IWritableMCHit *mcHit;
          COMET::IChannelId channelId(iwire);
          mcHit = CreateMCHit(channelId, parents[iwire], edepTot, arriveT);

          cdcCell_nHits++;
          int cellID = GetCellID(wireID,layerID);
          cdcCell_edep.push_back(edepwire[iwire]);
          cdcCell_x.push_back(poswire[iwire].X()/10.);
          cdcCell_y.push_back(poswire[iwire].Y()/10.);
          cdcCell_z.push_back(poswire[iwire].Z()/10.);
          GetWirePosition(iwire,xy,poswire[iwire].Z());
          cdcCell_wx.push_back(xy[0]/10.);
          cdcCell_wy.push_back(xy[1]/10.);
          cdcCell_wz.push_back(poswire[iwire].Z()/10.);
          cdcCell_t.push_back(timewire[iwire]); // ???
          cdcCell_driftD.push_back(ddistsm/10.);
          cdcCell_driftDtrue.push_back(ddistwire[iwire]/10.);
          cdcCell_cellID.push_back(cellID);
          cdcCell_layerID.push_back(layerID);
          cdcCell_px.push_back(momwire[iwire].X());
          cdcCell_py.push_back(momwire[iwire].Y());
          cdcCell_pz.push_back(momwire[iwire].Z());
          if (signalHits) cdcCell_hittype.push_back(0); // only signal is added for the moment
          else cdcCell_hittype.push_back(1);
          cdcCell_nclusters.push_back(nclusters[iwire]);
          cdcCell_absmom.push_back(momwire[iwire].Mag());
        }

        std::vector<int> index(cdcCell_nHits);
        TMath::Sort(cdcCell_nHits, &cdcCell_absmom[0], &index[0]);

        int jhit  = 0, phit = -1;
        int nturn = 1;
        double dphi = 0, dz = 0.;
        //////// sort the signal hits in momentum ordering
        for (int ihit = 0; ihit < cdcCell_nHits; ihit++) {
          jhit = index[ihit];
          if (cdcCell_hittype[jhit]!=0) continue;
          if (ihit>0 && useOnlyFirstTurn) {
            phit = index[ihit-1]; // previous hit
            dphi = atan2(cdcCell_y[jhit], cdcCell_x[jhit]) - atan2(cdcCell_y[phit], cdcCell_x[phit]);
            dz   = fabs(cdcCell_z[jhit] - cdcCell_z[phit]);
            if (fabs(dphi)>5.*deg2rad && fabs(dphi-two_pi)>5.*deg2rad && fabs(dphi+two_pi)>5.*deg2rad) {
              //break;
              nturn++;
            }
          }
          CdcCell_nHits++;
          if (CdcCell_maxLayer<cdcCell_layerID[jhit]) CdcCell_maxLayer = cdcCell_layerID[jhit];
          if (nturn==1) CdcCell_nHits1st++;
          CdcCell_edep.push_back(cdcCell_edep[jhit]); 
          CdcCell_x.push_back(cdcCell_x[jhit]); 
          CdcCell_y.push_back(cdcCell_y[jhit]); 
          CdcCell_z.push_back(cdcCell_z[jhit]); 
          CdcCell_wx.push_back(cdcCell_wx[jhit]); 
          CdcCell_wy.push_back(cdcCell_wy[jhit]); 
          CdcCell_wz.push_back(cdcCell_wz[jhit]); 
          CdcCell_t.push_back(cdcCell_t[jhit]); 
          CdcCell_driftD.push_back(cdcCell_driftD[jhit]); 
          CdcCell_driftDtrue.push_back(cdcCell_driftDtrue[jhit]); 
          CdcCell_cellID.push_back(cdcCell_cellID[jhit]); 
          CdcCell_layerID.push_back(cdcCell_layerID[jhit]); 
          CdcCell_px.push_back(cdcCell_px[jhit]/1e3); // GeV
          CdcCell_py.push_back(cdcCell_py[jhit]/1e3); // GeV
          CdcCell_pz.push_back(cdcCell_pz[jhit]/1e3); // GeV
          CdcCell_hittype.push_back(cdcCell_hittype[jhit]); // only signal is added for the moment
          CdcCell_nclusters.push_back(cdcCell_nclusters[jhit]);
          CdcCell_nturn.push_back(nturn);
        }
        CdcCell_nTurns = nturn;

        // fill the noise hits....
        for (int ihit = 0; ihit < cdcCell_nHits; ihit++) {
          jhit = index[ihit];
          if (cdcCell_hittype[jhit]==0) continue;
          CdcCell_nHits++;
          CdcCell_edep.push_back(cdcCell_edep[jhit]); 
          CdcCell_x.push_back(cdcCell_x[jhit]); 
          CdcCell_y.push_back(cdcCell_y[jhit]); 
          CdcCell_z.push_back(cdcCell_z[jhit]); 
          CdcCell_wx.push_back(cdcCell_wx[jhit]); 
          CdcCell_wy.push_back(cdcCell_wy[jhit]); 
          CdcCell_wz.push_back(cdcCell_wz[jhit]); 
          CdcCell_t.push_back(cdcCell_t[jhit]); 
          CdcCell_driftD.push_back(cdcCell_driftD[jhit]); 
          CdcCell_driftDtrue.push_back(cdcCell_driftDtrue[jhit]); 
          CdcCell_cellID.push_back(cdcCell_cellID[jhit]); 
          CdcCell_layerID.push_back(cdcCell_layerID[jhit]); 
          CdcCell_px.push_back(cdcCell_px[jhit]/1e3); // GeV
          CdcCell_py.push_back(cdcCell_py[jhit]/1e3); // GeV
          CdcCell_pz.push_back(cdcCell_pz[jhit]/1e3); // GeV
          CdcCell_hittype.push_back(cdcCell_hittype[jhit]);
          CdcCell_nclusters.push_back(cdcCell_nclusters[jhit]);
          CdcCell_nturn.push_back(-1); // don't care about the #of turn of noise hit
        }

        CdcCell_mt      = 0;
        CdcCell_eventID = nGoodEvent;
        /*if (nGoodEvent<10) {
          TCanvas *tmpc = new TCanvas("tmpc", "tmpc", 800, 800);
          tmpc->cd();
          TH1F   *tmpfXY = gPad->DrawFrame(-routCDC/10.-5., -routCDC/10.-5., +routCDC/10.+5., +routCDC/10.+5.);
          TGraph *tmpgXY = new TGraph(CdcCell_nHits, &CdcCell_x[0], &CdcCell_y[0]);
          tmpgXY->SetMarkerColor(kBlue);
          tmpgXY->SetLineColor(kBlue);
          tmpgXY->Draw("LP same");
          TArrow *vmomXY = new TArrow(CdcCell_x[0], CdcCell_y[0], CdcCell_x[0]+100.*CdcCell_px[0], CdcCell_y[0]+100.*CdcCell_py[0], 0.01, "|>");
          vmomXY->SetLineWidth(2);
          vmomXY->SetLineColor(kCyan);
          vmomXY->Draw();
          tmpc->Print(Form("CyDetSignalEventDisplayXY%d.eps",nGoodEvent));
          tmpc->cd();
          TH1F   *tmpfXZ = gPad->DrawFrame(-routCDC/10.-5., -zmax/10., +routCDC/10.+5., +zmax/10.);
          TGraph *tmpgXZ = new TGraph(CdcCell_nHits, &CdcCell_x[0], &CdcCell_z[0]);
          tmpgXZ->SetMarkerColor(kBlue);
          tmpgXZ->SetLineColor(kBlue);
          tmpgXZ->Draw("LP same");
          TArrow *vmomXZ = new TArrow(CdcCell_x[0], CdcCell_z[0], CdcCell_x[0]+100.*CdcCell_px[0], CdcCell_z[0]+100.*CdcCell_pz[0], 0.01, "|>");
          vmomXZ->SetLineWidth(2);
          vmomXZ->SetLineColor(kCyan);
          vmomXZ->Draw();
          tmpc->Print(Form("CyDetSignalEventDisplayXZ%d.eps",nGoodEvent));

          delete tmpfXY;
          delete tmpfXZ;
          delete tmpgXY;
          delete tmpgXZ;
          delete vmomXY;
          delete vmomXZ;
          delete tmpc;
        }*/

        tree->Fill();
        nGoodEvent++;
      }

      CdcCell_eventID = 0;
      CdcCell_nHits   = 0;
      CdcCell_nHits1st= 0;
      CdcCell_nTurns  = 0;
      CdcCell_maxLayer= 0;
      CdcCell_t.clear();
      CdcCell_wx.clear();
      CdcCell_wy.clear();
      CdcCell_wz.clear();
      CdcCell_x.clear();
      CdcCell_y.clear();
      CdcCell_z.clear();
      CdcCell_px.clear();
      CdcCell_py.clear();
      CdcCell_pz.clear();
      CdcCell_driftD.clear();
      CdcCell_driftDtrue.clear();
      CdcCell_tof.clear();
      CdcCell_tstart.clear();
      CdcCell_tstop.clear();
      CdcCell_cellID.clear();
      CdcCell_layerID.clear();
      CdcCell_edep.clear();
      CdcCell_posflag.clear();
      CdcCell_hittype.clear();
      CdcCell_nclusters.clear();
      CdcCell_nturn.clear();
      CdcCell_mx = 0;
      CdcCell_my = 0;
      CdcCell_mz = 0;
      CdcCell_mt = 0;
      CdcHits_nHits = 0;
      CdcHits_x.clear();
      CdcHits_y.clear();
      CdcHits_z.clear();

      CthHits_nHits = 0;
      CthHits_nHitsI= 0;
      CthHits_x.clear();
      CthHits_y.clear();
      CthHits_z.clear();
      CthHits_t.clear();
      CthHits_px.clear();
      CthHits_py.clear();
      CthHits_pz.clear();
      CthHits_edep.clear();

      return true;
    }

    /// Called after the arguments are processes by before reading the first
    /// event.  The output file is open so any histograms will be added to the
    /// output file.
    virtual void Initialize(void) {
      nGoodEvent = 0;
      resolution = 0.2*mm;
      rnd = new TRandom3();

      innerWall = TArc(0,0,rinCDC);
      outerWall = TArc(0,0,routCDC);
      innerWall.SetLineStyle(2);
      innerWall.SetLineWidth(2);
      innerWall.SetFillStyle(0);
      outerWall.SetLineStyle(2);
      outerWall.SetLineWidth(2);
      outerWall.SetFillStyle(0);

      TString fileName = Form("test%s",ofileName.c_str());
      newFile = new TFile(fileName.Data(),"RECREATE");
      fc1= new TCanvas("c1","Electron");
      numSignal=0;

      tree = new TTree("tree", "cdc outputs");

      //--- added variables here (O.hayashi 14th Nov 2014)
      CdcCell_eventID = 0;
      CdcCell_nHits   = 0;
      CdcCell_nHits1st= 0;
      CdcCell_nTurns  = 0;
      CdcCell_maxLayer= 0;
      CdcCell_mx = 0;
      CdcCell_my = 0;
      CdcCell_mz = 0;
      CdcCell_mt = 0;

      momvtxx = 0; // momentum @ vertex
      momvtxy = 0; // momentum @ vertex
      momvtxz = 0; // momentum @ vertex
      posvtxx = 0; // position @ vertex
      posvtxy = 0; // position @ vertex
      posvtxz = 0; // position @ vertex

      tree->Branch("CdcCell_eventID",&CdcCell_eventID);
      tree->Branch("CdcCell_nHits",&CdcCell_nHits);
      tree->Branch("CdcCell_nHits1st",&CdcCell_nHits1st);
      tree->Branch("CdcCell_nTurns",&CdcCell_nTurns);
      tree->Branch("CdcCell_maxLayer",&CdcCell_maxLayer);
      tree->Branch("CdcCell_t",&CdcCell_t);
      tree->Branch("CdcCell_wx",&CdcCell_wx);
      tree->Branch("CdcCell_wy",&CdcCell_wy);
      tree->Branch("CdcCell_wz",&CdcCell_wz);
      tree->Branch("CdcCell_x",&CdcCell_x);
      tree->Branch("CdcCell_y",&CdcCell_y);
      tree->Branch("CdcCell_z",&CdcCell_z);
      tree->Branch("CdcCell_px",&CdcCell_px);
      tree->Branch("CdcCell_py",&CdcCell_py);
      tree->Branch("CdcCell_pz",&CdcCell_pz);
      tree->Branch("CdcCell_driftD",&CdcCell_driftD);
      tree->Branch("CdcCell_driftDtrue",&CdcCell_driftDtrue);
      tree->Branch("CdcCell_tof",&CdcCell_tof);
      tree->Branch("CdcCell_tstart",&CdcCell_tstart);
      tree->Branch("CdcCell_tstop",&CdcCell_tstop);
      tree->Branch("CdcCell_cellID",&CdcCell_cellID);
      tree->Branch("CdcCell_layerID",&CdcCell_layerID);
      tree->Branch("CdcCell_edep",&CdcCell_edep);
      tree->Branch("CdcCell_posflag",&CdcCell_posflag);
      tree->Branch("CdcCell_hittype",&CdcCell_hittype);
      tree->Branch("CdcCell_nclusters",&CdcCell_nclusters);
      tree->Branch("CdcCell_nturn",&CdcCell_nturn);
      tree->Branch("CdcCell_mx",&CdcCell_mx);
      tree->Branch("CdcCell_my",&CdcCell_my);
      tree->Branch("CdcCell_mz",&CdcCell_mz);
      tree->Branch("CdcCell_mt",&CdcCell_mt);
      //tree->Branch("CdcHits_nHits",&CdcHits_nHits);
      //tree->Branch("CdcHits_x",&CdcHits_x);
      //tree->Branch("CdcHits_y",&CdcHits_y);
      //tree->Branch("CdcHits_z",&CdcHits_z);
      tree->Branch("CthHits_nHits",&CthHits_nHits);
      tree->Branch("CthHits_nHitsI",&CthHits_nHitsI);
      tree->Branch("CthHits_x",&CthHits_x);
      tree->Branch("CthHits_y",&CthHits_y);
      tree->Branch("CthHits_z",&CthHits_z);
      tree->Branch("CthHits_t",&CthHits_t);
      tree->Branch("CthHits_px",&CthHits_px);
      tree->Branch("CthHits_py",&CthHits_py);
      tree->Branch("CthHits_pz",&CthHits_pz);
      tree->Branch("CthHits_edep",&CthHits_edep);
      tree->Branch("momvtxx",&momvtxx);
      tree->Branch("momvtxy",&momvtxy);
      tree->Branch("momvtxz",&momvtxz);
      //tree->Branch("posvtxx",&posvtxx);
      //tree->Branch("posvtxy",&posvtxy);
      //tree->Branch("posvtxz",&posvtxz);
      //-------------------------------------------------

      HitPosHist= new TH2F("2D Position of hits in the CDC","Hit positions",400,-routCDC-1.*cm,+routCDC+1.*cm,400,-routCDC-1.*cm,+routCDC+1.*cm);
      HitPosHist->SetXTitle("X (mm)");
      HitPosHist->SetYTitle("Y (mm)");

      EneDepoHist= new TH1F("Energy Deposit in one cell","Energy deposit",400,0,10);
      EneDepoHist->SetXTitle("Energy deposit (keV)");
      EneDepoHist->SetYTitle("Entries");

      DDistHist= new TH1F("Drift Distance in one cell","Drift Distance",400,0,12);
      DDistHist->SetXTitle("Drift distance (mm)");
      DDistHist->SetYTitle("Entries");

      nWires = 0;
      for (int iLayer = 0; iLayer < nLayer; iLayer++) {
        nWires += GetNumberOfWirePerLayer(iLayer);
      }
      fXYWire.resize(nWires);
      wirePos.resize(nWires);

      int wirenum = 0;
      int numWirePerLayer = 0;
      double zpos_wire; // @ EP
      double rpos_wire; // @ EP
      double phi_wire;
      double deltaPhi;
      double alphaAngle; // projection of stereo angle @ XY-plane
      double xEnd[2] = {0,0};
      double yEnd[2] = {0,0};
      double zEnd[2] = {0,0};
      for (int iLayer = 0; iLayer < nLayer; iLayer++) {
        std::cout << iLayer << "-th layer" << std::endl;
        if (iLayer>0) wirenum += numWirePerLayer;
        numWirePerLayer = GetNumberOfWirePerLayer(iLayer);
        deltaPhi   = two_pi/((double)numWirePerLayer);
        if(iLayer<14)  alphaAngle = (iLayer%2!=0)? deltaPhi*((double)innerskipHole) : (-deltaPhi*((double)innerskipHole));
        if(iLayer>=14) alphaAngle = (iLayer%2!=0)? deltaPhi*((double)outerskipHole) : (-deltaPhi*((double)outerskipHole));
        rpos_wire  = rminLayer + (((double)iLayer)+0.5)*hCell;
        zpos_wire  = rpos_wire*tan(tiltAngle) + lenOffset/2.;
        std::cout << " #of sense wire = " << numWirePerLayer << ", z length = " << 2.*zpos_wire << std::endl;

        int wireID;
        double a_x, a_y, x_0, y_0; // parameters for the wire function
        double minX, maxX;
        for (int icell = 0; icell < numWirePerLayer; icell++) {
          wireID   = wirenum + icell;
          double PhiOffset = 0;
          if(iLayer == 0 || iLayer == 1 || iLayer == 3 || iLayer == 6 || iLayer == 8 || iLayer == 10 || iLayer == 12 || iLayer == 16){
            PhiOffset = deltaPhi/2.;
          }
          phi_wire = ((double)icell) * deltaPhi + PhiOffset;
          if (icell==0) std::cout << " start phi = " << phi_wire*rad2deg << ", alpha angle = " << alphaAngle*rad2deg << std::endl;
          xEnd[0]  = rpos_wire * cos(phi_wire);
          yEnd[0]  = rpos_wire * sin(phi_wire);
          zEnd[0]  = zpos_wire;
          xEnd[1]  = rpos_wire * cos(phi_wire+alphaAngle);
          yEnd[1]  = rpos_wire * sin(phi_wire+alphaAngle);
          zEnd[1]  = -zpos_wire;
          wirePos[wireID] = TVector3((xEnd[0]+xEnd[1])/2., (yEnd[0]+yEnd[1])/2., 0);
          a_x      = (xEnd[1]-xEnd[0])/(zEnd[1]-zEnd[0]);
          a_y      = (yEnd[1]-yEnd[0])/(zEnd[1]-zEnd[0]);
          x_0      = wirePos[wireID].X();
          y_0      = wirePos[wireID].Y();
          minX     = std::min(xEnd[0],xEnd[1]);
          maxX     = std::max(xEnd[0],xEnd[1]);
          fXYWire[wireID] = new TF1(Form("fXYWire_%d",wireID), funcXYWire, minX, maxX, 4);
          fXYWire[wireID]->SetNpx(10000);
          fXYWire[wireID]->SetParameters(a_x, a_y, x_0, y_0);
        }
      }
    }


    /// Called after reading the last event.  The output file is still open,
    /// so you can add extra information.  Because of an ideosyncracy in the
    /// way root handles histograms, objects created in Initialize() will
    /// already be stored in the output file.
    virtual void Finalize(COMET::ICOMETOutput* output) {
      fc1->cd();
      HitPosHist->SetMarkerStyle(8);
      HitPosHist->SetMarkerSize(.5);
      HitPosHist->SetMarkerColor(kRed);
      HitPosHist->Draw();
      innerWall.Draw("only");
      outerWall.Draw("only");
      //fc1->Print("HitPos2D.png");
      fc1->cd();
      EneDepoHist->Draw();
      //fc1->Print("EnergyDeposit.png");
      fc1->cd();
      DDistHist->Draw();
      //fc1->Print("DriftDistance.png");

      newFile->cd();
      HitPosHist->Write();
      EneDepoHist->Write();
      DDistHist->Write();
      tree->Write();
      newFile->Close();

      std::cout<<"Number of electrons: " << numSignal <<std::endl;
    }

private:

    TArc innerWall;
    TArc outerWall;
    TFile    *newFile;
    TTree    *tree;
    TCanvas  *fc1;
    TRandom3 *rnd;
    TH1F *EneDepoHist;
    TH1F *DDistHist;
    TH2F *HitPosHist;
    int nWires; // number of total wires
    std::vector < TF1* >     fXYWire; // wire function
    std::vector < TVector3 > wirePos; // wire position @ z=0

    //--- Add variables here (14th Nov 2014 O.hayashi)
    int CdcCell_eventID;
    int CdcCell_nHits;
    int CdcCell_nHits1st;
    int CdcCell_nTurns;
    int CdcCell_maxLayer;
    std::vector < double > CdcCell_t;
    std::vector < double > CdcCell_wx;
    std::vector < double > CdcCell_wy;
    std::vector < double > CdcCell_wz;
    std::vector < double > CdcCell_x;
    std::vector < double > CdcCell_y;
    std::vector < double > CdcCell_z;
    std::vector < double > CdcCell_px;
    std::vector < double > CdcCell_py;
    std::vector < double > CdcCell_pz;
    std::vector < double > CdcCell_driftD;
    std::vector < double > CdcCell_driftDtrue;
    std::vector < double > CdcCell_tof;
    std::vector < double > CdcCell_tstart;
    std::vector < double > CdcCell_tstop;
    std::vector < int >    CdcCell_cellID;
    std::vector < int >    CdcCell_layerID;
    std::vector < double > CdcCell_edep;
    std::vector < int >    CdcCell_posflag;
    std::vector < int >    CdcCell_hittype;
    std::vector < int >    CdcCell_nclusters;
    std::vector < int >    CdcCell_nturn;
    double CdcCell_mx;
    double CdcCell_my;
    double CdcCell_mz;
    double CdcCell_mt;
    int CdcHits_nHits;
    std::vector < double > CdcHits_x; // low hits in CDC
    std::vector < double > CdcHits_y; // low hits in CDC
    std::vector < double > CdcHits_z; // low hits in CDC

    double momvtxx; // momentum @ vertex
    double momvtxy; // momentum @ vertex
    double momvtxz; // momentum @ vertex
    double posvtxx; // position @ vertex
    double posvtxy; // position @ vertex
    double posvtxz; // position @ vertex

    int CthHits_nHits;
    int CthHits_nHitsI; // indirect hit
    std::vector < double > CthHits_x;
    std::vector < double > CthHits_y;
    std::vector < double > CthHits_z;
    std::vector < double > CthHits_t;
    std::vector < double > CthHits_px;
    std::vector < double > CthHits_py;
    std::vector < double > CthHits_pz;
    std::vector < double > CthHits_edep;
    //----------------------------------------------------

    int nGoodEvent;
    double resolution;
    Long_t numSignal;
};

int main(int argc, char **argv) {
    IUserLoop userLoop;
    ofileName = argv[1];
    ofileName = ofileName.erase(0,26);
    std::cout << "output file name = " << "test" << ofileName << std::endl;
    cometEventLoop(argc,argv,userLoop);
}