//
// 23-apr-2002 for calculation radiation length with helping
// Nov 25, 2005; 
// Mar 04,2006 - 20% reduction 
// Mar 24,2006 - 10% reduction(pbth90); calculate(13, 1);

#define SCINT 20
#define EMCALMAT  21

void calculate(Int_t imat=13, int var=1);
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 

Int_t   lSC=77;
//Int_t   lSC=61; // 20% reduction
Float_t thicknessSC=0.16;
Int_t   lP=77;
//Int_t   lP=61;
Float_t thicknessP=0.16;
Float_t radlPB=0.562; // in cm
Float_t thicknessW = thicknessSC*lSC + thicknessP*lP;

class TGeant3TGeo; TGeant3TGeo *g3=0; // move from TGeant3TGeo to TGeant3TGeoTGeo

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, int var)
{
  if(imat<12 || imat>14) {
     printf("<E> wrong materials imat = %2i\n", imat);
     return;
  }
  initMaterials(imat);
  if(var == 0) { // 20% reduction; 1.28mm/1.92mm
    lSC = lP = 77;
    thicknessSC = 0.192;
    thicknessP  = 0.128;
  } else if(var == 1) { // 10% reduction; 1.44mm/1.76mm
    lSC = lP = 77;
    thicknessSC = 0.176;
    thicknessP  = 0.144;
  }
  thicknessW = thicknessSC*lSC + thicknessP*lP;

  // 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);
  // Nov 25,05 !!
  Float_t rhoTmp  = (rho[0]*thicknessSC*lSC + rho[1]*thicknessP*lP)/(thicknessSC*lSC + thicknessP*lP);

  thicknessEMCAL = thicknessSC*lSC + thicknessP*lP;
  Char_t ctmp[21];
  sprintf(ctmp, "%s + SC$", chW[imat-12]);

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

  g3->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(" number of radiation length %7.2f \n", thicknessEMCAL/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(!g3) {
    gSystem->Load("libgeant321");
    g3 = new TGeant3TGeo("C++ Interface to Geant3");
    g3->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
    g3->Gsmixt(SCINT, "Polystyrene$", aTmp, zTmp, dTmp, -2, wTmp);

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

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

// calculate(13)