#include "bos.h"
#include "fpack.h"
#include "look.h"
#include "h1util.h"
#include "mdb.h"
#include "4vectors.h"
#include <iostream.h>
#include <iomanip.h>
#include <signal.h>
#include <time.h> 
#include <math.h> 
#include "banks.h"
#include "dstbanks.h"

// define GTRA bank structure for "good track" flags
struct gtrarow {
	int goodtrack;
	int goodcentral;
	int goodforward;
};

// define QTRA bank structure for track selection steering
struct qtrarow {
	float thetamin_cen;	
	float thetamax_cen;	
	float RSmax_cen;	
	float REmin_cen;	
	float RLmin_cen;	
	float ptmin_cen;	
	float pmin_cen;	
	float thetamin_fwd;	
	float thetamax_fwd;	
	float zSmax_fwd;	
	float zEmin_fwd;	
	float zLmin_fwd;	
	float ptmin_fwd;	
	float pmin_fwd;	
	float thetamin_cen2;	
	float thetamax_cen2;	
	float RSmax_cen2;	
	float REmin_cen2;	
	float RLmin_cen2;	
	float ptmin_cen2;	
	float pmin_cen2;

//	Constructor defines default values;
//	zero default values can be omitted.
	qtrarow () {		
	  thetamax_cen = 180;	
	  RSmax_cen    = 100;	
	  thetamax_fwd = 180;	
	  zSmax_fwd    = 250;	
	  thetamax_cen2= 180;	
	  RSmax_cen2 =   200;	
	}	
};

void trksel() {
/*	Simple track delection:
	Read parameters from QTRA bank,
	decide for each track if it is good or very good,
	the resolve ambiguities by going along list of equivalent tracks
	and keeping only the best (i.e. the one with the highest
	quality)
*/

        static const float raddeg = 180/3.14159265;
        static int first = 1;
        static qtrarow qtra;
	
	if (first) {
	  first = 0;
	  try {
	    StrBank<qtrarow> QTRA ("QTRA");
	    if (QTRA.is_open()) {
  	      qtra = *QTRA;
	      cout << "QTRA bank found. Parameters taken from there:\n";
	    } else {
	      cout << "QTRA bank not found. Default parameters taken:\n";
	    }	  
	  }
	  catch (StrBankLengthError) {
	    cout << "QTRA bank has wrong format. Default parameters taken:\n";
	  }
          cout << "qtra.thetamin_cen = "<< qtra.thetamin_cen << '\n';
          cout << "qtra.thetamax_cen = "<< qtra.thetamax_cen << '\n';
          cout << "qtra.RSmax_cen    = "<< qtra.RSmax_cen << '\n';
          cout << "qtra.REmin_cen    = "<< qtra.REmin_cen << '\n';
          cout << "qtra.RLmin_cen    = "<< qtra.RLmin_cen << '\n';
          cout << "qtra.ptmin_cen    = "<< qtra.ptmin_cen << '\n';
          cout << "qtra.pmin_cen     = "<< qtra.pmin_cen << '\n';
          cout << "qtra.thetamin_fwd = "<< qtra.thetamin_fwd << '\n';
          cout << "qtra.thetamax_fwd = "<< qtra.thetamax_fwd << '\n';
          cout << "qtra.zSmax_fwd    = "<< qtra.zSmax_fwd << '\n';
          cout << "qtra.zEmin_fwd    = "<< qtra.zEmin_fwd << '\n';
          cout << "qtra.zLmin_fwd    = "<< qtra.zLmin_fwd << '\n';
          cout << "qtra.ptmin_fwd    = "<< qtra.ptmin_fwd << '\n';
          cout << "qtra.pmin_fwd     = "<< qtra.pmin_fwd << '\n';
	}
	
	Table<dtrarow> dtra ("DTRA", 0); 
	if (dtra.is_open()) {
	  Table<gtrarow> gtra ("GTRA", 0, dtra.rows());
	  if (gtra.is_open()) {
	    for (int i = 1; i <= dtra.rows(); i++) {
	      dtrarow &track = dtra[i];
	      if (track.is_central()) {
	        float rs;
	        float re;
	        dtnvrow &dtnvtrack = *track.l_tcn_tr;
	        if (track.l_ver_tr->ityp_ve <= 1 &&
	            track.q_tr != 0 &&
	            track.thet_tr*raddeg > qtra.thetamin_cen &&
	            track.thet_tr*raddeg < qtra.thetamax_cen &&
	            (rs = dtnvtrack.rs_tn) < qtra.RSmax_cen &&
	            (re = dtnvtrack.getrend()) > qtra.REmin_cen &&
	            re - rs > qtra.RLmin_cen &&
	            1./track.ptinv_tr > qtra.ptmin_cen &&
	            track.getp() > qtra.pmin_cen) {
	          gtra[i].goodcentral = 1;
	          
	          if (track.l_ver_tr->ityp_ve <= 1 &&
	              track.q_tr != 0 &&
	              track.thet_tr*raddeg > qtra.thetamin_cen2 &&
	              track.thet_tr*raddeg < qtra.thetamax_cen2 &&
	              (rs = dtnvtrack.rs_tn) < qtra.RSmax_cen2 &&
	              (re = dtnvtrack.getrend()) > qtra.REmin_cen2 &&
	              re - rs > qtra.RLmin_cen2 &&
	              1./track.ptinv_tr > qtra.ptmin_cen2 &&
	                track.getp() > qtra.pmin_cen2) {
	            gtra[i].goodcentral = 2;
	            
	          }
	        }
	          
              } else {
	        gtra[i].goodcentral = 0;
	      }
	      if (track.is_forward()) {
	        float zs;
	        float ze;
	        dtnvrow &dtnvtrack = *track.l_tnf_tr;
	        gtra[i].goodforward =
	          track.l_ver_tr->ityp_ve <= 1 &&
	          track.q_tr != 0 &&
	          track.thet_tr*raddeg > qtra.thetamin_fwd &&
	          track.thet_tr*raddeg < qtra.thetamax_fwd &&
	          (zs = dtnvtrack.z0_tn) < qtra.zSmax_fwd &&
	          (ze = dtnvtrack.zend_tn) > qtra.zEmin_fwd &&
	          ze - zs > qtra.zLmin_fwd &&
	          1./track.ptinv_tr > qtra.ptmin_fwd &&
	          track.getp() > qtra.pmin_fwd;
	      } else {
	        gtra[i].goodforward = 0;
	      }
	      gtra[i].goodtrack = 
	        gtra[i].goodcentral > gtra[i].goodforward ?
	        gtra[i].goodcentral : gtra[i].goodforward; 
	    }
	    for (int i = 1; i <= dtra.rows(); i++) {
	      if (gtra[i].goodtrack) {
	        TablePointer<dtrarow> best = i;
	        for (TablePointer<dtrarow> P = dtra[i].l_nxn_tr; P; 
	             P = P->l_nxn_tr) {
	          if (gtra[P].goodtrack <=
	              gtra[best].goodtrack) {
	            gtra[P].goodtrack = 0;
	          } else {
	            gtra[best].goodtrack = 0;
	            best = P;
	          }
	        }
	      }
	    }
	    blist (bcs, "E+", "GTRA");
	  } else {
	    cout << "Error: Cannot open GTRA bank!\n";
	  }
	}
}


int main()
{
//---catch all exceptions in order to be able to clean up properly
  try {
        cout << "main() starts...\n";

//---synchronize output streams of Fortran/C and C++
	ios::sync_with_stdio ();

//---initialization of BOS is now done when defining an BOS array
	 
//---read i/o commands and steering banks
	fparmr (5);
	breadc ();
	 
//---check for input file
	char *inbos;

	if (ifrst ("BOSINPUT", 0)) 
	  inbos  = "BOSINPUT";
	else {
	  inbos = " ";
	  cout << "Error: No input file found!\n";
	  h1stop ();
	}
	 
//---signal handler for SGI
	signal (SIGTERM, (void (*)(...))fseqe);  // for kill[-15] pid
	signal (SIGUSR1, (void (*)(...))fseqe);  // for kill -16  pid

//---get job number from database
	int jobnr = mdb ("/JOB");

//---initialize random number generator
//   if the h1 data base is not written to (i.e. for private production)
//   then JOBNR=0 for all jobs. in this case the number of secomds since
//   some date in the past is used as seed. 
	h1rnin ((jobnr <= 10) ? time (0): jobnr + 1);
			 
	int iret;
	while (fseqr ("BOSINPUT", iret) >= 0) {
	  if (bosmdl.begjob) {
//	    Initialize histograms	  
	    bhs (1, 0, 50, 0., 5.);
	    stext (1, 1, "\\i{m\\_{\\g{pp}}} [GeV]");
	    stext (1, 4, "2-Prong Mass Spectrum, Assuming Pion Masses");
	    bhs (2, 0, 50, 0., 5.);
	    stext (2, 1, "\\i{m\\_{KK}} [GeV]");
	    stext (2, 4, "2-Prong Mass Spectrum, Assuming Kaon Masses");
	  }
	  if (bosmdl.revent) { 
	    try {
	  
//	      Track selection; produces bank GTRA	  
	      trksel ();
	             
	      Table<dtrarow> dtra ("DTRA", 0); 
	      Table<gtrarow> gtra ("GTRA", 0); 

/*	      Loop over DTRA-Tracks; count good tracks and very good
	      tracks; add up 4-vectors of very good tracks, usinng pion
	      (default) and kaon mass hypotheses.
*/	    
  	      fourvector vm;
	      fourvector vmK;
	      int ngoodtracks = 0;
	      int ngoodtracks2 = 0;
	      for (int i=1; i <= dtra.rows(); i++) {
	        if (gtra.is_open() && gtra[i].goodtrack) {
                  ngoodtracks++;
	          if (gtra[i].goodtrack >= 2) {
		    vm += dtra[i].getfourvector();
		    vmK += dtra[i].getfourvector(fourvector::mK);
		  }
	        }
	      }

//	      For 2-prongs, fill mass into histogram	    	    
              if (ngoodtracks == 2 && ngoodtracks2 == 2) {
	        shs (1, 0, vm.getm());
	        shs (2, 0, vmK.getm());
	      }
	    }
	    catch (BankError& BE) // catch by reference to keep info
	                          // on derived exceptions
	    { 
	      cout << "A bank error occured:\n";
	      BE.printmessage();
	      cout << "Event is ignored.\n";
            }
//          The following "catch" should not be necessary, because
//	    TableWidthError is a BankError. But the SUN g++ compiler 
//	    does not handle this correctly.
	    catch (TableWidthError& BE) {
	      cout << "A TableWidthError error occured:\n";
	      BE.printmessage();
	      cout << "Event is ignored.\n";
            }
	  }
	  if (bosmdl.endjob) { 
	    lstore (9, "2prong.look");
	  }
	}
	
  }
  catch (...) {
        cout << "main() finished due to an uncaught exception!\n";
  }

//--- clean up, even after uncaught exception

//	bosta ();
	 
	fparm ("CLOSE ALL");
	mdb ("/END");
	cout << "main() finished\n";
	return 0;
}