// This macro is written to investigate the qualities of tracks 
// at the position of the detector the detector position are taken
// from phast

// The special purpose of this program is to determine the efficiency of the
// chambers

#include <iostream>
#include <cmath>
#include <map>
#include "TH1.h"
#include "TH2.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "Phast.h"
#include "PaSetup.h"
#include "PaEvent.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "G3part.h"

#include "W45P.h"

//#include "W45Plane.cc"

// If you want current event to be saved to output tree 
// (microDST) call at least once "e.TagToSave()" function

//w45 represents the full detector system it might become a class later on


void UserEvent12(PaEvent& e)
{

    
  static map <string, int> plane;
  static vector <W45P> w45;
  static double alpha, x, y, z, u, v, tt, t ,tsig ,eu, ev, ex, ey, d, dd;
  static int chn, dID, extra;
  //  static char dName[20];

  static bool first(true);
  // create histos here as static 
  static TTree* tree(NULL);

  if(first){ // histograms and Ntupes booking block
    
    tree = new TTree("nt","User Ntuple example"); // name (has to be unique) and title of the Ntuple
    tree->Branch("extra", &extra, "extra/I");
    tree->Branch("alpha",    &alpha, "alpha/D");  
    tree->Branch("x",    &x, "x/D");  
    tree->Branch("y",    &y, "y/D");  
    tree->Branch("z",    &z, "z/D");  
    tree->Branch("u",    &u, "u/D");  
    tree->Branch("v",    &v, "v/D");  
    tree->Branch("d",    &d, "d/D");  
    tree->Branch("dd",  &dd, "dd/D");  
    tree->Branch("ex",  &ex, "ex/D");  
    tree->Branch("ey",  &ey, "ey/D");  
    tree->Branch("eu",  &eu, "eu/D");  
    tree->Branch("ev",  &eu, "ev/D");  
    tree->Branch("t",    &t, "t/D");  
    tree->Branch("tt",  &tt, "tt/D");  
    tree->Branch("tsig",  &tsig, "tsig/D");  
    tree->Branch("chn", &chn, "chn/I");  
    tree->Branch("dID", &dID, "dID/I");  
    
    // end of histo booking 

    cout << "Detector Initialisation ..." << endl;
    
    const int nDet = PaSetup::Ref().NDetectors();    
    for ( int det = 0, index = 0 ; det < nDet ; det++ ){
      const PaDetect &d = PaSetup::Ref().Detector(det);
      // the detectors are ordered after their z position
      // i.e. we can start to put them in a list at position 
      string dn = d.Name();
      if ( dn.find("DW") != string::npos ){	
	w45.push_back(W45P(d.Nwires(),det,dn,d.Ca(),d.Sa(),
			  d.Uorig(),d.X(),d.Y(),d.Z(),d.XSiz(),d.YSiz()));
	plane[dn] = index++; 
      }
    }
    cout << "Data processing ..." << endl;
    first = false;
  } // end of histogram booking
  
  //loop over all tracks
  for ( int it = 0 ; it < e.NTrack() ; it++ ){
    const PaTrack &tr = e.vTrack(it);


    // cuts on track quality && after SM2  
    PaTPar helix; // buffer to store best helix for extrapolation
    //bool intrapol = true; //to mark tracks as extrapolated or
    //interpolated
    
    // search best helix for tracks to be extrapolated
    if  ( tr.NTPar() < 5 ){ 
      extra = 1;
      if ( tr.vTPar(0)(0) < 1200 && tr.vTPar(tr.NTPar()-1)(0) > 2000
	   && tr.vTPar(tr.NTPar()-1)(0) < 3000)
	helix = tr.vTPar(tr.NTPar()-1) ;
      else if ( tr.vTPar(0)(0) > 1900 && tr.vTPar(0)(0) < 2200 &&
		tr.vTPar(tr.NTPar()-1)(0) > 4000 )
	helix = tr.vTPar(0) ;
      else // bull_shit track ==> next track
	continue;
    }
    else
      extra = 0;

    //loop over all detectors and extrapolate track to the plane
    for ( unsigned int detId = 0 ; detId < w45.size() ; detId++ ){	 
      
      //track reconstructed before and after plane?
      //smoothing @ 2860, 3004.5, 3034, 3126
      //i.e. plane 0- 3 helix #1 = planeID/4+1
      //           4- 7 helix #2
      //           8-11 helix #3
      //          12-15 helix #4 
      
      // search best helix for tracks to be interpolated
      if ( tr.NTPar() > 4 )
	  helix = tr.vTPar(detId/4+1) ;

      PaTPar pos, rotpos;
      
      //extrapolate  result stored in pos
      helix.Extrapolate(w45[detId].getZ(), pos ,true);
      //rotate it to wire reference system (x = u)
      pos.Rotate(w45[detId].getCa(), w45[detId].getSa(), rotpos);  
	

      // track parameter
      
      alpha = rotpos(3);
      x = pos(1);
      y = pos(2);
      z = pos(0);
      u = rotpos(1);
      v = rotpos(2);
      eu = sqrt(rotpos(1,1));
      ev = sqrt(rotpos(2,2));
      tt = tr.MeanTime();
      tsig = w45[detId].getPropTime(w45[detId].getName(),pos);
      ex = sqrt(rotpos(1,1)); 
      ey = sqrt(rotpos(2,2)); 
      if ( eu > 0.15 ) continue;
      
      // hit parameter
      chn = w45[detId].getWire(rotpos(1));
      dID = detId;
      double RealDetID = w45[detId].getId();

      const vector< PaHit > &hit = e.Hits();
      for ( unsigned int nhit = 0 ; nhit < hit.size() ; nhit++ )
	{
	  if ( hit[nhit].iDet() != RealDetID ) continue;
	  const vector <PaDigit> &raw = hit[nhit].vDigits();
	  if ( chn - raw[0].IWire() != 0 ) continue;	
	  d = rotpos(1) - w45[detId].getPos(chn);
	  dd = d / sqrt(1 + rotpos(3) * rotpos(3) );
	  t = raw[0].vDigInfo()[0] ;
	  if ( t < -50 || t > 550 ) continue;
	  tree->Fill();
	} // end loop over hits

//       const vector< PaDigit > &dig = e.RawDigits();
//       for ( unsigned int ndig = 0 ; ndig < dig.size() ; ndig++ )
//  	{
//  	  if ( dig[ndig].DecodingMapName() != w45[detId].getName() ) continue;
//  	  if ( chn - dig[ndig].IWire() != 0 ) continue;	
//  	  d = rotpos(1) - w45[detId].getPos(chn);
//  	  dd = d / sqrt(1 + rotpos(3) * rotpos(3) );
// 	  t = dig[ndig].vDigInfo()[4] ;
// 	  tree->Fill();
// 	} // end loop over digits
      
    }
  }
}