#include <fstream>
#include "TMath.h"
#include "TGeoBBox.h"
#include "PamVMCDetTrk.h"

using TMath::RadToDeg;

ClassImp(PamVMCDetTrk)

void  PamVMCDetTrk::DefineGeometry(){

  TGeoMaterial *m=0;
  TGeoMixture *mx =0;
  TString s = "ALUMINIUM"; 
  m = GetMat(s);
  if(!m){
    m = new TGeoMaterial(s,26.98,13,2.7);
    m->SetUniqueID(   9);
  }

  s = "IRON";
  m = GetMat(s);
  if(!m){
    m = new TGeoMaterial(s,55.85,26,7.87);
    m->SetUniqueID(  10);
  }


  s = "N2GAS";
  m = GetMat(s);
  if (!m){
    m = new TGeoMaterial(s,14.01,7,0.1250000E-02);
    m->SetUniqueID(  19);
  }

  s = "SILICON TR";
  m= GetMat(s);
  if(!m){
    m = new TGeoMaterial(s,28.09,14,2.33);
    m->SetUniqueID(  20);
  }

  s = "G10C";
  mx = (TGeoMixture*)GetMat(s);
  if(!mx){
    mx = new TGeoMixture(s,4,   1.70000    );
    mx->SetUniqueID(  27);
    mx->DefineElement(0,28.09,14,0.53);
    mx->DefineElement(1,16,8,0.3);
    mx->DefineElement(2,12.01,6,0.15);
    mx->DefineElement(3,1.01,1,0.2000000E-01);
  }

  s = "SBAR";
  mx = (TGeoMixture*)GetMat(s);
  if(!mx){
    mx = new TGeoMixture(s,3,   1.49000    );
    mx->SetUniqueID(  40);
    //Carbon fiber 90%, Resin 10%
    mx->DefineElement(0,12.01,6,0.845*0.9+0.545*0.1); //C
    mx->DefineElement(1,1.01,1,0.043*0.9+0.092*0.1); //H
    mx->DefineElement(2,16.0,8,0.112*0.9+0.363*0.1); //O
  }


  TGeoMedium * n = 0;

  s="ALUMINIUM";
  n = (GetMed(s))?GetMed(s): new TGeoMedium(s,1,9,-1,1,40,1,100,0.5000000E-01,0.136283,0.5000000E-03);

  s = "N2 GAS";
  n = (GetMed(s))?GetMed(s): new TGeoMedium(s,3,19,-1,1,40,1,0.5000000E-01,0.5000000E-01,0.1000000E-02,0.5000000E-03);

  s = "IRONTR";
  n = (GetMed(s))?GetMed(s): new TGeoMedium(s,5,10,-1,1,40,1,0.5000000E-01,0.5000000E-01,0.8000000E-01,0.5000000E-03); 

  s = "FIBER+RESIN";
  n = (GetMed(s))?GetMed(s): new TGeoMedium(s,29,40,-1,1,40,1,100,0.5000000E-01,0.5000000E-02,0.5000000E-03); 
 
  s = "SITRACKER";
  n = (GetMed(s))?GetMed(s): new TGeoMedium(s,12,20,1,1,40,1,0.1000000E-02,0.5000000E-01,0.5000000E-04,0.1569056E-01);

  s = "G10C";
  n = (GetMed(s))?GetMed(s):new TGeoMedium(s,17,27,-1,1,40,1,100,0.5000000E-01,0.1200000E-02,0.5000000E-03);

  TGeoVolume *v = 0;

  
  s = "SPEB";
  v = (GetVol(s))?GetVol(s):gGeoManager->MakeBox(s,GetMed("N2 GAS"),15.,14.4,22.615);   
  for(Int_t i = 1; i<=5; i++){
    s = "MGF";
    s += i;
    v = (GetVol(s))?GetVol(s):gGeoManager->MakeBox(s,GetMed("ALUMINIUM"),15.,14.4,4.45); 
  }
  s = "MGPL";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("IRONTR"),12.,11.4,4.); 
  s = "MGPI"; 
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("N2 GAS"),8.07,6.57,4.); 
  for(Int_t i = 1; i<=6; i++){
    s = "TRP";
    s += i;
    v = (GetVol(s))?GetVol(s):gGeoManager->MakeBox(s,GetMed("ALUMINIUM"),10.49,12.15,0.365);
  }
  for(Int_t i = 1; i<=6; i++){
    s = "TPA";
    s += i;
    v = (GetVol(s))?GetVol(s):gGeoManager->MakeBox(s,GetMed("N2 GAS"),8.2,9.9,0.365); 
  }
  s = "TPAI";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("ALUMINIUM"),8.144,1.25,0.5E-01); 
  s = "TRSL";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("SITRACKER"),2.6665,3.5,0.15E-01);
  s = "THBP";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("G10C"),2.6665,2.75,0.15E-01); 
  s = "TSPA";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("SITRACKER"),2.596,3.413,0.15E-01); 
  s = "TRCP";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("FIBER+RESIN"),0.35E-01,9.9,0.25); 
  s = "TBAL"; 
  v = (GetVol(s))?GetVol(s): gGeoManager->MakeBox(s,GetMed("ALUMINIUM"),1.49,2.25,0.365); 
  s = "MGPA";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("ALUMINIUM"),8.1,6.6,4.);
  s = "TPGA";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("N2 GAS"),10.29,10.65,0.4E-01); 
  s = "TPGI";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("N2 GAS"),8.07,6.57,0.15E-01); 
  s = "TPGU";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("N2 GAS"),10.49,12.15,0.1E-01); 
  s = "TPGD";
  v = (GetVol(s))?GetVol(s):  gGeoManager->MakeBox(s,GetMed("N2 GAS"),12.,11.4,0.5E-02); 


  ReadPaddlesPos();
   
  //Original, not shifted positions of the paddles TRSL respect to TPAs volumes.
  Double_t x14,x36,x25, y123, y456, za; 
  x36 = 5.405;
  x14 = -x36;
  x25 = 0.;
  y123 = 6.252;
  y456 = -0.75;
  za = 0.15E-01;
  Double_t xp[6] = {x14, x25, x36, x14, x25, x36};
  Double_t yp[6] = {y123, y123, y123, y456, y456, y456};
  Double_t x,y,z;
 

  GetVol("TRSL")->AddNode(GetVol("TSPA"),1,gGeoIdentity); //adding sensetive detector first
  GetVol("MGPL")->AddNode(GetVol("MGPA"),1,gGeoIdentity);
      GetVol("MGPA")->AddNode(GetVol("MGPI"),1,gGeoIdentity);

  Double_t zMGF =  ((TGeoBBox*)GetVol("SPEB")->GetShape())->GetDZ();
  for(Int_t i = 1; i<=6; i++){
    Int_t rot_id, j;
    TString ss="TRP";
    if( i<6 ){
      s = "MGF";
      s += i;
      zMGF -= ((TGeoBBox*)GetVol(s)->GetShape())->GetDZ();
      GetVol("SPEB")->AddNode(GetVol( s ),1,new TGeoTranslation(0.,0., zMGF )); 
      zMGF -= ((TGeoBBox*)GetVol(s)->GetShape())->GetDZ();
          GetVol( s )->AddNode(GetVol("TPGD"),1,new TGeoTranslation(0.,0.,-4.445));
	  GetVol( s )->AddNode(GetVol("TPGI"),1,new TGeoTranslation(0.,0.,-4.425));
	  GetVol( s )->AddNode(GetVol("MGPL"),1,new TGeoTranslation(0.,0.,-0.41)); 
	  GetVol( s )->AddNode(GetVol("TPGI"),2,new TGeoTranslation(0.,0.,3.605)); 
	  GetVol( s )->AddNode(GetVol("TPGA"),1,new TGeoTranslation(0.,-0.75,3.66));
	  ss += i;
	  GetVol( s )->AddNode(GetVol( ss ),1,new TGeoTranslation(dxflo[i-1], -2.25 + dyflo[i-1], 4.065));
    } else {
      s = "SPEB";
      ss += i;
      GetVol( s )->AddNode(GetVol( ss ),1,new TGeoCombiTrans(dxflo[i-1], 2.25 + dyflo[i-1], -22.25, GetRot("rot10")));
    }	
    s = "TPA";
    s += i;
    GetVol( ss )->AddNode(GetVol( s ),i,new TGeoTranslation(0.,-0.75,0.)); 
    GetVol( s )->AddNodeOverlap(GetVol("TRCP"),1,new TGeoTranslation(-8.1075,0.,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("TRCP"),2,new TGeoTranslation(-2.7025,0.,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("TRCP"),3,new TGeoTranslation(2.7025,0.,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("TRCP"),4,new TGeoTranslation(8.1075,0.,0.15E-01)); 
    //paddles//
    j = (i-1)*NP_TRK + 1;
    Int_t kk;
   
    for(Int_t pad = 0; pad<6; pad++){
      ss = "rot_";
      rot_id = 50 + j;
      ss += rot_id;
      if(i<6) kk = pad; 
      else{
	j = (i-1)*NP_TRK + pad + 1;
	if(pad<=2)
	  kk= 3+pad;
	else
	  kk = pad-3;	  
    
	ss = "rot_";
	rot_id = 50 + kk + (i-1)*NP_TRK + 1;
	ss += rot_id;
	(fctmap.find(rot_id - 50)->second)->fdy *=-1.;
	(fctmap.find(rot_id - 50)->second)->fdz *=-1.;
	//cout<<"kk+1="<<kk+1<<"j="<<j<<endl;
      }
      x = xp[pad] + (fctmap.find(rot_id - 50)->second)->fdx;
      y = yp[pad] + (fctmap.find(rot_id - 50)->second)->fdy;
      z = za  +  (fctmap.find(rot_id - 50)->second)->fdz;
      
      //cout<<"plane ="<<i<<" TRSL="<<j<<"x,y,z"<<x<<","<<y<<","<<z<<"RotName"<<ss<<endl;
      GetVol( s )->AddNode(GetVol("TRSL"),j ,new TGeoCombiTrans( x,y,z, GetRot( ss ))); 
      j++;
    }
      
    GetVol( s )->AddNodeOverlap(GetVol("THBP"),1,new TGeoTranslation(-5.405,-7.0015,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("THBP"),2,new TGeoTranslation(0.,-7.0015,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("THBP"),3,new TGeoTranslation(5.405,-7.0015,0.15E-01)); 
    GetVol( s )->AddNodeOverlap(GetVol("TPAI"),1,new TGeoTranslation(0.,-8.65,0.315)); 

    if(i<6){
      s = "MGF";
      s += i;
      GetVol( s )->AddNode(GetVol("TPGU"),1,new TGeoTranslation(0.,-2.25,4.44)); 
    }
  }
     
  GetVol("SPEB")->AddNode(GetVol("TBAL"),1,new TGeoTranslation(5.4,-12.15,-22.25));
  GetVol("SPEB")->AddNode(GetVol("TBAL"),2,new TGeoTranslation(-5.4,-12.15,-22.25));


  SetMotherProp(GetVol("SPEB"),1,new TGeoTranslation(0.,0.,49.274)); 
};

void  PamVMCDetTrk::DefineCuts(){

  TString  s ="ALUMINIUM"; // default GPAMELA parameters
  if (GetMed(s) && !GetCC(s))
  SetCC(s, new pCutControl(GetMedID(s), 0.0001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
			     0.001, 0.01, 1., 1., 1., 0., 1., 1., 1., 4., 1., 1., 1., 1.));
  s ="N2 GAS";
  if (GetMed(s) && !GetCC(s))
    SetCC(s, new pCutControl(GetMedID(s), 0.0001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
                             0.001, 0.01, 1., 1., 1., 0., 0., 1., 1., 4., 1., 1., 2., 1.));
  s = "IRONTR";  // default GPAMELA parameters
  if (GetMed(s) && !GetCC(s))
    SetCC(s, new pCutControl(GetMedID(s), 0.0001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
			     0.001, 0.01, 1., 1., 1., 0., 1., 1., 1., 4., 1., 1., 1., 1.));  
  s = "SITRACKER";
  if (GetMed(s) && !GetCC(s))
    SetCC(s, new pCutControl(GetMedID(s), 0.0001, 0.001, 0.00001, 0.00001, 0.0001, 0.0001, 0.00001, 0.00001,
                             0.001, 0.01, 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 2.));
  s = "FIBER+RESIN";
  if (GetMed(s) && !GetCC(s)) // default GPAMELA parameters
    SetCC(s, new pCutControl(GetMedID(s), 0.0001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
                             0.001, 0.01, 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 1.));

  s = "G10C";
  if (GetMed(s) && !GetCC(s))
    SetCC(s, new pCutControl(GetMedID(s), 0.00001, 0.00001, 0.00001, 0.00001, 0.001, 0.001, 0.01, 0.00001,
                             0.00001, 0.01, 1., 1., 1., 0., 1., 1., 1., 1., 1., 1., 1., 2.));
}

void  PamVMCDetTrk::ReadPaddlesPos(){
  //This method reads from file rotations and shifts of the paddles
  //The format is:  Omega, Beta, Gamma, dx, dy, dz
  
  cout<<"Loading POS files for TRK sensors..."<<endl;
  ifstream fcfile;
  TString PAM_VMC=gSystem->Getenv("PAM_VMC");
  TString PLATFORM=gSystem->Getenv("PLATFORM");
  TString PATH=PAM_VMC+"/lib/tgt_"+PLATFORM+"/parameters_l";
  Int_t ls, ss;

  TMap rotshiftmap;

  for(Int_t ik = 1; ik<=NP_TRK; ik++){
    TString filename = PATH;
    if( ik < 4 ){
      ls = ik;
      ss = 2;
    } else{
      ls = ik-3;
      ss = 1;
    }
    filename += ls;
    filename += "s";
    filename += ss;
    filename += ".dat";
    cout<<"  loading:  "<<filename<<endl;
    fcfile.open(filename);
    if(!fcfile){
      cout<<" !!! ERROR, file "<<filename<<" is missed...exiting"<<endl;
      exit(-1);
    }
    //reading from file, create structures and put them in array

    Double_t om, be, ga, dx, dy, dz;
    for (Int_t ih = 1; ih <= NP_TRK; ih++){
      fcfile >> om;
      fcfile >> be;
      fcfile >> ga;
      fcfile >> dx;
      fcfile >> dy;
      fcfile >> dz;

      Int_t ii = (6-ih)*6 +ik;
      fctmap[ii] = new pCombTrans( om, be, ga, dx, dy, dz );
    }
    fcfile.close();
  }

  for (Int_t ij = 1; ij<=6; ij++){
    dxflo[ij-1] = dyflo[ij-1] = 0.;
    for (Int_t il = 1; il<=6; il++){
      Int_t ii = (ij-1)*6+il;
      dxflo[ij-1] += (fctmap.find(ii)->second)->fdx;
      dyflo[ij-1] += (fctmap.find(ii)->second)->fdy;
    }
    dxflo[ij-1] /=6.;
    dyflo[ij-1] /=6.;
  }

  for (Int_t ij = 1; ij<=6; ij++){
    for (Int_t il = 1; il<=6; il++){
      Int_t ii = (ij-1)*6+il;
      (fctmap.find(ii)->second)->fdx = ( (fctmap.find(ii)->second)->fdx - dxflo[ij-1] )*1.E-4;
      (fctmap.find(ii)->second)->fdy = ( (fctmap.find(ii)->second)->fdy - dyflo[ij-1] )*1.E-4;
      (fctmap.find(ii)->second)->fdz = ( (fctmap.find(ii)->second)->fdz - 100. )*1.E-4;
    }
    dxflo[ij-1] *= 1.E-4;
    dyflo[ij-1] *= 1.E-4;
  }
  
  for (Int_t ij = 1; ij<=6; ij++){
    for (Int_t il = 1; il<=6; il++){
      Int_t ii = (ij-1)*6+il;
      
      Double_t om_1 = (fctmap.find(ii)->second)->fomega * 1.E-6;
      Double_t be_1 = (fctmap.find(ii)->second)->fbeta * 1.E-6;
      Double_t ga_1 = (fctmap.find(ii)->second)->fgamma * 1.E-6;
      
      if(ij==6){
	om_1 = -om_1;
	ga_1 = -ga_1;
      }

      TVector3 rot1(1.,om_1, -ga_1), rot2(-om_1,1.,be_1), rot3(ga_1,-be_1,1.);
      Int_t rot_no = 50 + ii;
      TString nmrot = "rot_";
      nmrot += rot_no;
   
      AddRot(nmrot.Data(), new TGeoRotation( nmrot.Data(), 
      				     rot1.Theta()*RadToDeg(),rot1.Phi()*RadToDeg(),
      				     rot2.Theta()*RadToDeg(),rot2.Phi()*RadToDeg(),
      				     rot3.Theta()*RadToDeg(),rot3.Phi()*RadToDeg() ));
      //cout<<"plno:"<<ij<<" padno:"<<il
      //  <<" t1:"<<rot1.Theta()*RadToDeg()<<" p1:"<<rot1.Phi()*RadToDeg()
      //  <<" t2:"<<rot2.Theta()*RadToDeg()<<" p2:"<<rot2.Phi()*RadToDeg()
      //  <<" t3:"<<rot3.Theta()*RadToDeg()<<" p3:"<<rot3.Phi()*RadToDeg() <<endl;
    }
  }
}