//
// 23-apr-2202 for calculation radiation length with helping
// AliGeant3 (envelope for GEANT3)
//
//#include "AliGeant3.h"
#define SCINT 20
#define EMCALMAT  21

void calculate(Int_t imat=12);
void initMaterials(Int_t imat);

// See output from Gprint("MATE").
// imat = 12 (TUNGSTEN)
//        13 (LEAD)
//        14 (URANIUM)

// see Particle Physics, p.171
Double_t Es=21.0; // in mev 
Double_t Ec;      // will be calculate 

const Int_t   lSC=21;
const Float_t thicknessSC=0.5;
const Int_t   lP=20;
Float_t thicknessP=0.5;
Float_t radlPB=0.562; // in cm
Float_t thicknessW = thicknessSC*lSC + thicknessP*lP;

class AliGeant3; AliGeant3 *aliG3=0;

Float_t ubuf[100], nbuf;

Char_t *chW[3]={" W","Pb"," U"};
Float_t aSC, zSC, rhoSC, radlSC, abslSC;
Char_t  chSC[21];
Float_t  aP,  zP, rhoP,  radlP,  abslP; // passive material
Char_t  chP[21];
Float_t a[2], z[2], rho[2], w[2];

Float_t aEMCAL, zEMCAL, rhoEMCAL, radlEMCAL, abslEMCAL, mollRadEMCAL, thicknessEMCAL;

void calculate(Int_t imat)
{
  if(imat<12 || imat>14) {
     printf("<E> wrong materials imat = %2i\n", imat);
     return;
  }
  initMaterials(imat);

  // Calculate weights and density
  thicknessP = thicknessSC;
  if(imat!=13) {
    thicknessP = (radlP/radlPB) * thicknessSC;
  }
  Float_t wSC = rhoSC*thicknessSC*lSC;
  Float_t wP  = rhoP*thicknessP*lP;
  w[0] = wSC/(wSC+wP);
  w[1] = wP /(wSC+wP);
  Float_t rhoTmp  = rho[0]*w[0] + rho[1]*w[1];
  thicknessEMCAL = thicknessSC*lSC + thicknessP*lP;
  Char_t ctmp[21];
  sprintf(ctmp, "%s + SC$", chW[imat-12]);

  aliG3->Gsmixt(EMCALMAT, ctmp, a, z, rhoTmp, 2, w);
  aliG3->Gprint("MATE");

  aliG3->Gfmate(EMCALMAT, chSC, aEMCAL, zEMCAL, rhoEMCAL, radlEMCAL, abslEMCAL, ubuf, nbuf);

  
  Ec = 610./(zEMCAL+1.24); // in mev
  mollRadEMCAL = radlEMCAL*Es/Ec;

  printf(" #layer Sc %2i ; thickness %5.3f(cm)\n", lSC, thicknessSC); 
  printf(" #layer %s %2i ; thickness %5.3f(cm)\n", chW[imat-12], lP, thicknessP);
  printf("        thickness EMCAL %7.2f (%5.1f%%)\n",  
  thicknessEMCAL, 100*(thicknessEMCAL/thicknessW));
  printf("        density         %7.2f(g/cm**3) \n", rhoEMCAL);  
  printf("        abs length      %7.2f(cm) \n", abslEMCAL);  
  printf("        rad length      %7.2f(cm) \n", radlEMCAL);  
  printf("        Moliere R       %7.2f(cm) \n", mollRadEMCAL);
  printf(" Es %5.2f(MeV) : Ec %7.2f(MeV) \n", Es, Ec);
}

void initMaterials(Int_t imat)
{ 
  if(!aliG3) {
    aliG3  = new AliGeant3("C++ Interface to Geant3"); // See config.C
    aliG3->Gmate();         // define standard material 

  // see AliEMCAL.cxx
  // --- The polysterene scintillator (CH) ---
    Float_t aTmp[2] = {12.011, 1.00794} ;
    Float_t zTmp[2] = {6.0, 1.0} ;
    Float_t wTmp[2] = {1.0, 1.0} ;
    Float_t dTmp    = 1.032 ;     // density
    aliG3->Gsmixt(SCINT, "Polystyrene$", aTmp, zTmp, dTmp, -2, wTmp);

    aliG3->Gfmate(SCINT, chSC, aSC, zSC, rhoSC, radlSC, abslSC, ubuf, nbuf);
    a[0]   = aSC;
    z[0]   = zSC;
    rho[0] = rhoSC;
  }

  aliG3->Gfmate(imat, chP, aP, zP, rhoP, radlP, abslP, ubuf, nbuf);  
  a[1]   = aP;
  z[1]   = zP;
  rho[1] = rhoP;
}