#include <TRandom.h>

#include "TGenVertexML.h"

//
// Author: Mikhail Zhabitsky <http://cern.ch/jabitski>
//________________________________________________________________________
//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TGenVertexML                                                         //
//                                                                      //
// Generates vertex position on a multilayer target                     //
// uniformly over z and according to 2 gaussians on x and y             //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

ClassImp(TGenVertexML)


 TGenVertexML::TGenVertexML( Int_t NL, double zth[], double zpos[],
			    double sigmax, double sigmay )
 : TGenVertex( 0.0, sigmax, sigmay )
{
  SetNameTitle( "TGenVertexML", "Multilayer target" );

  fNL = (NL>1) ? NL : 1;

  if ( (!zth || !zpos) && (fNL!=12) ) {
    fprintf( stderr, "Default parameters are valid only for NL=12!n");
    fNL=12;
  };
  fzth = new double[fNL];
  fzpos = new double[fNL];
  fzthint = new double[fNL];

  if (zth) 
    for (int i=0; i<fNL; i++)
      fzth[i] = zth[i];
  else {
    // Values adapted from GEANT-DIRAC 2.64_00 15/12/2005
    fzth[0] = 2*3.970E-4;
    fzth[1] = 2*3.970E-4;
    fzth[2] = 2*3.990E-4;
    fzth[3] = 2*3.970E-4;
    fzth[4] = 2*3.990E-4;
    fzth[5] = 2*3.980E-4;
    fzth[6] = 2*3.985E-4;
    fzth[7] = 2*3.985E-4;
    fzth[8] = 2*3.980E-4;
    fzth[9] = 2*3.995E-4;
    fzth[10] = 2*4.000E-4;
    fzth[11] = 2*4.000E-4;
  };

  fzthint[0] = fzth[0];
  for (int i=1; i<fNL; i++)
    fzthint[i] = fzthint[i-1] + fzth[i];
  fzthick = fzthint[fNL-1];

  if (zpos) 
    for (int i=0; i<fNL; i++)
      fzpos[i] = zpos[i];
  else {
    double zsum = 0.0;
    for (int i=0; i<fNL; i++) {
      fzpos[i] = zsum + 0.5*fzth[i];
      zsum = fzpos[i] + 0.1;
    };
  };
}


 TVector3 TGenVertexML::GenVertexTargetSR() const
{
  // in target system of reference
  double xtarget[3];
  for (int j=0; j<2; j++)
    xtarget[j] = gRandom->Gaus( 0, fbsigma[j] );

  double z = fzthint[fNL-1]*gRandom->Rndm();
  int l=0;
  for ( l=0; l<fNL; l++)
    if (z<fzthint[l])
      break;

  xtarget[2] = fzpos[l] + 0.5*fzth[l] + z - fzthint[l];

  TVector3 vect( xtarget );

  return(vect);
}


 void TGenVertexML::Print( Option_t *option ) const
{
  TNamed::Print(option);

  printf("Multilayer target: total thickness=%gn", fzthint[fNL-1] );
  for (int i=0; i<fNL; i++)
    printf( "Layer %i: [%.6g,%.6g]n",
	    i+1, fzpos[i]-0.5*fzth[i], fzpos[i]+0.5*fzth[i] );
  printf("Gaussians: sigmaX=%g,tsigmaY=%gn", fbsigma[0], fbsigma[1] );

  return;
}


 TGenVertexML::~TGenVertexML()
{
  delete[] fzth;
  delete[] fzpos;
  delete[] fzthint;
}