#include "rootAnc.h"
#include "PISAEventHeader.h"
#include "MvdPISAHit.h"

void rootAncMvd(Int_t ancflag, PISAEventHeader *EventHeader)
{
//  26-Oct-2000 JPS modified
//  calculate r,phi,z,wfradr (correctly) where phi, theta were calculated 
//  (incorrectly) previously. "wfradr" is a variable used in the MVD code
//  to encode the address of the Si wafer which was hit.
//  Also removed code which did not save multiple hits from the same
//  track (this is normal and expected for the MVD).

//  For the barrel part of the MVD, r is just the "radius" of the barrel
//  in which the hit is seen, phi is at the center of the Si wafer,
//  z does not include loss of resolution from strip pitch.
//  For the pads, the loss of resolution in r, phi from the pad granularity
//  is not included. 

  static int icall = 0;
  static TFile *AncMvdFile = 0;
  static TNtuple *AncMvdNtuple = 0;

  if(ancflag == 1) {
    //
    // Check that the file has been opened
    //
    if(icall == 0 || AncMvdFile == 0){
      cerr << "\n rootAncMvd <E> bad call with ancflag = 1" << endl;
      exit(1);
    }  // safety check
    AncMvdFile->Write();
    AncMvdFile->Close();
    return;
  }  // check for close output file flag

  const int NTPL_PARAM = 24;
  float evt_ntpl[NTPL_PARAM];
  if(icall == 0){
    //
    // NTUPLE files
    //
    AncMvdFile = new TFile("ancmvd.root", "recreate", "MVD PISA NTUPLE");
    AncMvdNtuple = new TNtuple("AncMvd", "MVD Ancestors",
			       "TRACK:NFILE:PTOT:PTHETA:PPHI:"//
			       "R_VERTEX:Z_VERTEX:THET_VER:PHI_VER:"//
			       "ITPARENT:IDPARENT:IDPART:IHIT:NHIT:"// 
			       "ITORIGIN:IDORIGIN:HITCOUNT:R:PHI:Z:"//
			       "WFRADR:Z0_EVENT:B_IMPACT:EVENT");

  }  // initialization
  icall++;

  int verRows;
  int iRow;

  int true_track;
  int nfile;
  int error;
  float ptot;
  float ptheta;
  float pphi;
  float r_vertex;
  float z_vertex;
  float theta_vertex;
  float phi_vertex;
  int itparent;
  int idparent;
  int idpart;

  int itorigin;
  int idorigin;

  float phi;
  float xglobal, yglobal, zglobal;
  float rglobal;  // sqrt(xglobal*xglobal + yglobal*yglobal)
  int wfradr;     // wafer address
  int ntype;      // 0 = inner barrel, 1 = outer, 2 = north pads, 3 = south
  int itemp;      // dummy variable
  int nrow;       // phi index in barrel (0-5)
  int nwafer;     // z index in barrel (0-11)
  const float DEGRAD = 57.295779513;
  const float R_MVD_INNER_BARREL = 5.0;
  const float R_MVD_OUTER_BARREL = 7.5;
  const float Z_MVD_NORTH_PADS   = 34.915;
  const float Z_MVD_SOUTH_PADS   =-34.915;
  const float DZ_MVD_WAFERS      = 5.3;

  int nstore;
  int istore;
  const int MAXSTORE = 100000;
  int trueTimes[MAXSTORE];

  evt_ntpl[NTPL_PARAM - 1] = icall;
  nstore = 0;
  verRows = MvdPISAHit::GetMvdCount();
  MvdPISAHit *verghit = MvdPISAHit::GetMvdHitEvt();

//  The following loop looks at all tracks and saves them in an output
//  array. The original code (prior to 26 Oct 2000) stored a track only
//  the first time it hit the detector. Since it is normal for a track to
//  hit, for example, the inner and outer layer of the MVD, this cut
//  distorted the output. It was removed by commenting out the offending lines.

  for (iRow = 0; iRow < verRows ; iRow++) {

    true_track = verghit[iRow].GetMctrack();
    trueTimes[iRow] = 0; // start with zero because the
                         // loop below counting this track too
    if(true_track < 1)
      continue ;  // this is an error condition

    for (istore = 0; istore < verRows; istore++){
      if(true_track ==verghit[istore].GetMctrack() ){
        trueTimes[iRow]++;
      } // count previous hits caused by this track
    }

    dio_ptrkorigin(&true_track, &nfile, &error, &ptot, &ptheta, &pphi,
		  &r_vertex, &z_vertex, &theta_vertex, &phi_vertex, 
		  &itorigin, &idorigin, &idpart);  
/*
*  NOTE: idpart is primary particle id (=idorigin), not the input particle id.
*/
      
    dio_ptrkstack(&true_track, &nfile, &error, &ptot, &ptheta, &pphi,
		  &r_vertex, &z_vertex, &theta_vertex, &phi_vertex, 
		  &itparent, &idparent, &idpart);  
 
    evt_ntpl[NTPL_PARAM - 2] = EventHeader->GetImpactParameter();
    evt_ntpl[NTPL_PARAM - 3] = EventHeader->GetZvertex();
    evt_ntpl[0] = float(true_track);
    evt_ntpl[1] = float(nfile);
    evt_ntpl[2] = ptot;
    evt_ntpl[3] = ptheta;
    evt_ntpl[4] = pphi;
    evt_ntpl[5] = r_vertex;
    evt_ntpl[6] = z_vertex;
    evt_ntpl[7] = theta_vertex;
    evt_ntpl[8] = phi_vertex;
    evt_ntpl[9] = float(itparent);
    evt_ntpl[10] = float(idparent);
    evt_ntpl[11] = float(idpart);
    evt_ntpl[12] = iRow;
    evt_ntpl[13] = nstore;
    evt_ntpl[14] = itorigin;
    evt_ntpl[15] = idorigin;
    evt_ntpl[16] = trueTimes[iRow];

    wfradr  = verghit[iRow].GetWfr_Adr();
    
    ntype=wfradr/10000;
    ntype--;  // 0 or 1 = shell, 2 or 3 means N or S pads
  
    // The following if/else block converts the wafer address
    // and the x,y,z in local coordinates (stored in hits)
    // to r, phi, z in phenix coordinates. It should really
    // get some of the constants from someplace better than
    // the hardwired constants used here, but ...

    if(ntype <2) {

      // for barrel
      itemp  = wfradr-((ntype+1)*10000);
      nrow   = itemp/1000;          // 0-5 (phi)
      nwafer = itemp-nrow*1000;     // 1-12 (z)

      // converting from fortran to c index convention
      nwafer--; // now this is the panel number (0-11)
      nrow--;   // now this is the row number (0-5)

      // now convert from the (unfortunate) pisa convention in which
      // the first phi segment is usually at -120 deg to the official
      // offline convention in which the first phi segment is at phi=0.

      phi = -120. + (nrow*60.);   // get phi of this row
       
      zglobal = ((float)nwafer-5.5)*DZ_MVD_WAFERS + verghit[iRow].GetZin();

      if ( ntype == 0 ) rglobal = R_MVD_INNER_BARREL;   // inner shell
      if ( ntype == 1 ) rglobal = R_MVD_OUTER_BARREL;   // outer shell
    }
    else
    {
      // pads
      xglobal = verghit[iRow].GetXin();
      yglobal = verghit[iRow].GetYin();
      phi = DEGRAD*atan2(yglobal, xglobal);
      rglobal = sqrt(xglobal*xglobal + yglobal*yglobal);
      if ( ntype == 2 ) zglobal = Z_MVD_NORTH_PADS;   // north pads
      if ( ntype == 3 ) zglobal = Z_MVD_SOUTH_PADS;   // south pads
    }
    
    if( phi < -90.0 )
        phi = 360.0 + phi;

    evt_ntpl[17] = rglobal;
    evt_ntpl[18] = phi;
    evt_ntpl[19] = zglobal;
    evt_ntpl[20] = (float)wfradr;

    AncMvdNtuple->Fill(evt_ntpl);

  } // loop over iRow

  return;
}