// -------------------------------------------------------------------------
// -----                  CbmLitEnvironment source file               -----
// -----                  Created 16/07/07  by A. Lebedev               -----
// -------------------------------------------------------------------------
 


#include "CbmLitEnvironment.h"

#include "CbmRunAna.h"
#include "CbmBaseParSet.h"
#include "CbmRuntimeDb.h"
#include "CbmDetector.h"
#include "CbmGeoStsPar.h"
#include "CbmGeoRichPar.h"
#include "CbmGeoTrdPar.h"
#include "CbmGeoPassivePar.h"
 

CbmLitEnvironment* CbmLitEnvironment::fInstance = NULL;

// constructors
CbmLitEnvironment::CbmLitEnvironment():
   fThickWall(5.),
   fField(NULL),
   fNofStations(0),
   fNofLayers(0)
{
   fvMaterials.clear();
   fNofLayersPerStation.clear();;
}

// destructor
CbmLitEnvironment::~CbmLitEnvironment() 
{
   if (fInstance != NULL) delete fInstance;
   
   for (std::vector<CbmLitMaterial*>::iterator i = fvMaterials.begin(); 
          i != fvMaterials.end(); i++)  delete (*i);
   fvMaterials.clear();
}

//
CbmLitEnvironment* CbmLitEnvironment::Instance()
{
   if (fInstance == NULL) {
      fInstance = new CbmLitEnvironment();
   }
   return fInstance;
}

// Initialization
void CbmLitEnvironment::Init()
{

}

// Field
CbmField* CbmLitEnvironment::GetField()
{
   if (fField == NULL) {
      CbmRunAna *Run = CbmRunAna::Instance();
         
      if(NULL == Run) {
         cout << "-E- CbmLitEnvironment::GetField : "
            << "Run Ana is not instantiated" << endl;
         return NULL;
      }
      
      cout << "-I- CbmLitEnvironment::GetField : "
         << "Reading Magnetic Field " << endl;
      
      fField = (CbmField*) Run->GetField();
   }
   return fField;
}

//
void CbmLitEnvironment::GetTrdLayout(Int_t &nofStations, 
                                     Int_t &nofLayers,
                                     std::vector<Int_t> &nofLayersPerStation)
{
   if (nofStations == 0 || nofLayers == 0 || nofLayersPerStation.empty()) {
   
      // Get the pointer to the CbmRunAna object
      CbmRunAna* ana = CbmRunAna::Instance();
      if(NULL == ana) {
         cout << "-E- CbmLitEnvironment::ReadTrdLayoutt :"
            << " no CbmRunAna object!" << endl;
         return;
      }
      
      // Get the pointer to run-time data base
      CbmRuntimeDb* rtdb = ana->GetRuntimeDb();
      if(NULL == rtdb) {
         cout << "-E- CbmLitEnvironment::ReadTrdLayout :"
            << " no runtime database!" << endl;
         return;
      }
      
      // Get the pointer to container of base parameters
      CbmBaseParSet* baseParSet =
         (CbmBaseParSet*) rtdb->getContainer("CbmBaseParSet");
      if(NULL == baseParSet) {
         cout << "-E- CbmLitEnvironment::ReadTrdLayout :"
            << " no container of base parameters!" << endl;
         return;
      }
      
      // Get the pointer to detector list
      TObjArray *detList = baseParSet->GetDetList();
      if(NULL == detList) {
         cout << "-E- CbmLitEnvironment::ReadTrdLayout :"
            << " no detector list!" << endl;
         return;
      }
      
      // Find TRD detector
      CbmDetector* trd = (CbmDetector*) detList->FindObject("TRD");
      if(NULL == trd) {
         cout << "-E- CbmLitEnvironment::ReadTrdLayout :"
            << " no TRD detector!" << endl;
         return;
      }
      // Determine the geometry version
      TString name = trd->GetGeometryFileName();
      if(name.Contains("9")) {
         cout << "-I- CbmLitEnvironment::ReadTrdLayout :"
            << " TRD geometry : 3x3." << endl;
         fNofStations = 3;
         fNofLayersPerStation.resize(3);
         for (Int_t i = 0; i < 3; i++) fNofLayersPerStation[i] = 3;
         fNofLayers = 9;
      } else if(name.Contains("12") || name.Contains("standard") || name.Contains("segmented")) {
         cout << "-I- CbmLitEnvironment::ReadTrdLayout :"
            << " TRD geometry : 3x4." << endl;
         fNofStations = 3;
         fNofLayersPerStation.resize(3);
         for (Int_t i = 0; i < 3; i++) fNofLayersPerStation[i] = 4;
         fNofLayers = 12;
      } else if(name.Contains("6x2")) {
         cout << "-I- CbmLitEnvironment::ReadTrdLayout :"
            << " TRD geometry : 6x2." << endl;
         fNofStations = 6;
         fNofLayersPerStation.resize(6);
         for (Int_t i = 0; i < 6; i++) fNofLayersPerStation[i] = 2;
         fNofLayers = 12;
      } else if(name.Contains("trd_special")) {
         cout << "-I- CbmLitEnvironment::ReadTrdLayout :"
              << " TRD geometry : 4-3-3." << endl;
         fNofStations = 3;
         fNofLayersPerStation.resize(3);
         fNofLayersPerStation[0] = 4;
         fNofLayersPerStation[1] = 3;
         fNofLayersPerStation[2] = 3;
         fNofLayers = 10;
      }
   }
   nofStations = fNofStations;
   nofLayers = fNofLayers;
   nofLayersPerStation = fNofLayersPerStation;
}

//
void CbmLitEnvironment::GetStsLayout()
{

}

//
const std::vector<CbmLitMaterial*> & CbmLitEnvironment::GetMaterials()
{

   if (fvMaterials.empty()) {
   
      fvMaterials.clear();
      
      ReadTrdMaterials();
      ReadStsMaterials();
      ReadRichMaterials();
      
      sort(fvMaterials.begin(), fvMaterials.end(), CbmLitMaterial::CmpUp);
   /*   
      cout << "-I- CbmLitEnvironment::ReadMaterials(): " << endl;
      cout << " Material  : [ Name, Z,   Thickness,   RadLength,   Density,   Zeff/Aeff,   Rmin,   Rmax" << endl;
      
      for (UInt_t i = 0; i < fvMaterials.size(); i++) {
         
         cout << " Material  : [ "
            << fvMaterials[i]->GetName() << ", "
            << fvMaterials[i]->GetZ() << ", "
            << fvMaterials[i]->GetThickness() << ", " 
            << fvMaterials[i]->GetRadLength() << ", "
            << fvMaterials[i]->GetDensity() << ", "
            << fvMaterials[i]->GetZeff() / fvMaterials[i]->GetAeff() << "] "
            << endl;
      }   
     */
      cout << "-I- CbmLitEnvironment::ReadMaterials(): " 
           << "materials initialized." << endl;
   }
  return fvMaterials;
}


void CbmLitEnvironment::ReadTrdMaterials(TObjArray *Nodes)
{
   for(int i = 0; i < Nodes->GetEntries(); i++){
      
      CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (Nodes->At(i));
      if (!node) continue;

      TString SName = node->getShapePointer()->GetName();
      
      TArrayD *P = node->getParameters();
      CbmGeoMedium *material = node->getMedium();
       
      Int_t NofComponents = material->getNComponents();
      //cout << "   NofComponents = " << NofComponents << endl;
      Double_t Zeff = 0.;
      Double_t Aeff = 0.;
      Double_t W = 0.;
      for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
         Double_t p[3];
         material->getComponent(iComp, p);
         W += p[2];        
      } 
      for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
         Double_t p[3];
         material->getComponent(iComp, p);
         Aeff += p[0] * p[2];
         Zeff += p[1] * p[2];
      } 
      Aeff /= W;
      Zeff /= W;
      
      CbmLitCircleMaterial *m = new CbmLitCircleMaterial();
      m->SetDensity(material->getDensity());
      m->SetAeff(Aeff);
      m->SetZeff(Zeff);
               
      material->calcRadiationLength();
      m->SetRadLength(material->getRadiationLength());
      
      if(SName.Contains("PGON")) { 
         m->SetInnerRadius(P->At(5));
         m->SetOuterRadius(P->At(6));
         m->SetThickness(-2 * P->At(4));
         m->SetXCenter(0.);
         m->SetYCenter(0.);         
      } else
      if(SName.Contains("BOX")) m->SetThickness(2 * P->At(2));
      else continue;
      
      CbmGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm
      CbmGeoVector centerV = node->getCenterPosition().getTranslation();
      Double_t z;
      
      if (m->GetThickness() < fThickWall) z = nodeV.Z() + centerV.Z();
      else z = nodeV.Z()+ centerV.Z() + m->GetThickness() / 2.0;
      m->SetZ(z);

      m->SetName(SName);
      
      if(m->GetRadLength() > 10e10) continue;
        
      fvMaterials.push_back(m);

   }
}

void CbmLitEnvironment::ReadTrdMaterials() {

   cout << "CbmLitEnvironment::ReadTrdMaterials: " << endl;
   
   cout << "   reading TRD materials..." << endl;

   CbmRunAna *Run = CbmRunAna::Instance();
   CbmRuntimeDb *RunDB = Run->GetRuntimeDb();

   CbmGeoTrdPar *TrdPar = (CbmGeoTrdPar*) (RunDB->getContainer("CbmGeoTrdPar"));

   TObjArray *NodesActive = TrdPar->GetGeoSensitiveNodes();
   TObjArray *NodesPassive = TrdPar->GetGeoPassiveNodes();

   cout << "   TRD Active Nodes:" << endl;
   ReadTrdMaterials(NodesActive);
   cout << "   TRD Passive Nodes:" << endl;
   ReadTrdMaterials(NodesPassive);

}

void CbmLitEnvironment::ReadStsMaterials(TObjArray *Nodes) 
{
   
   for (int i = 0; i < Nodes->GetEntries(); i++){
         
         CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (Nodes->At(i));
         if ( !node ) break;
               
         TArrayD *P = node->getParameters();
         CbmGeoMedium *material = node->getMedium();
            
         Int_t NofComponents = material->getNComponents();
         Double_t Zeff = 0.;
         Double_t Aeff = 0.;
         Double_t W = 0.;
         for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
            Double_t p[3];
            material->getComponent(iComp, p);
            W += p[2];    
         } 
         for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
            Double_t p[3];
            material->getComponent(iComp, p);
            Aeff += p[0] * p[2];
            Zeff += p[1] * p[2];
         } 
         Aeff /= W;
         Zeff /= W;
         
         CbmLitCircleMaterial *m = new CbmLitCircleMaterial();
         m->SetDensity(material->getDensity());
         m->SetAeff(Aeff);
         m->SetZeff(Zeff);
                          
         material->calcRadiationLength();
              
         m->SetRadLength(material->getRadiationLength());
         m->SetThickness(P->At(2));
         m->SetInnerRadius(P->At(0));
         m->SetOuterRadius(P->At(1));
         m->SetXCenter(0.);
         m->SetYCenter(0.);   
                     
         CbmGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm
         CbmGeoVector centerV = node->getCenterPosition().getTranslation();
         Double_t z;
         if (m->GetThickness() < fThickWall) z = nodeV.Z() + centerV.Z();
         else z = nodeV.Z()+ centerV.Z() + m->GetThickness() / 2.0;
         m->SetZ(z);
         
         if(m->GetRadLength() > 10e10) continue;
   
         fvMaterials.push_back(m);
   }
}

void CbmLitEnvironment::ReadStsMaterials() {

   std::cout << "CbmLitEnvironment::ReadStsMaterials: " << std::endl;
   
   cout << "   reading STS materials..." << endl;

   CbmRunAna *Run = CbmRunAna::Instance();
   CbmRuntimeDb *RunDB = Run->GetRuntimeDb();

     CbmGeoStsPar* StsPar = (CbmGeoStsPar*) (RunDB->getContainer("CbmGeoStsPar"));

   TObjArray *NodesActive = StsPar->GetGeoSensitiveNodes();
   TObjArray *NodesPassive = StsPar->GetGeoPassiveNodes();

   if( StsPar ){ //=== Sts stations ===
      cout << "   STS Active Nodes: " 
           << NodesActive->GetEntries() << endl;
      ReadStsMaterials(NodesActive);
      cout << "   STS Passive Nodes: " 
           << NodesPassive->GetEntries() << endl;
      ReadStsMaterials(NodesPassive);      
   }
}

void CbmLitEnvironment::ReadRichMaterials() 
{
   std::cout << "CbmLitEnvironment::ReadRichMaterials: " << std::endl;
   
   CbmRunAna *Run = CbmRunAna::Instance();
   CbmRuntimeDb *RunDB = Run->GetRuntimeDb();
   
   CbmGeoRichPar *RichPar = (CbmGeoRichPar*) (RunDB->getContainer("CbmGeoRichPar"));
    
  if( RichPar ) {
    cout << "CbmLitEnvironment::ReadRichMaterials" << endl;

    CbmGeoNode* rich1entrance = NULL;             // entrance window
    CbmGeoNode* rich1exit = NULL;                 // exit window
    CbmGeoNode* rich1gas = NULL;
    CbmGeoNode* rich1mgl = NULL;
    
    TObjArray *NodesActive = RichPar->GetGeoSensitiveNodes();
    TObjArray *NodesPassive = RichPar->GetGeoPassiveNodes();
    
    cout << " RICH active nodes:" << endl;
    
    for (Int_t i = 0; i < NodesActive->GetEntries(); i++) { 

      CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (NodesActive->At(i));
      if ( !node ) continue;  
      TString name = node->getName();
      cout << name << endl;
      
      if(name == "rich1mgl#1") { // RICH mirror (currently sensitive! - will change soon!)
        rich1mgl = node;
      }
      
    }
    
    cout << " RICH passive nodes:" << endl;
    
    for (Int_t i = 0; i < NodesPassive->GetEntries(); i++) {

      CbmGeoNode *node = dynamic_cast<CbmGeoNode*> (NodesPassive->At(i));  
      
      if ( !node ) continue;
      
      TString name = node->getName();
      cout << name << endl;
      
      if(name == "rich1entrance") {                  // RICH entrance window
         rich1entrance = node;
      } else if(name == "rich1mgl#1") {              // RICH mirror (see above)
        rich1mgl = node;
      } else if(name == "rich1gas1") {               // RICH gas
        rich1gas = node;
      } else if(name == "rich1exit") {               // RICH exit window
        rich1exit = node;
      }
    }

    if(rich1entrance && rich1gas && rich1mgl && rich1exit ) {
      
      CbmLitCircleMaterial* kapton1 = new CbmLitCircleMaterial();
      if (0 == ReadRichNode( kapton1, rich1entrance))  
         fvMaterials.push_back(kapton1);
      
      CbmLitCircleMaterial* kapton2 = new CbmLitCircleMaterial();
      if (0 == ReadRichNode( kapton2, rich1exit)) 
         fvMaterials.push_back(kapton2);
      
      CbmLitCircleMaterial* mirror = new CbmLitCircleMaterial();
      if (0 == ReadRichNode( mirror, rich1mgl))
         fvMaterials.push_back(mirror);;
      
      CbmLitCircleMaterial* gas1 = new CbmLitCircleMaterial();
      if ( 0 == ReadRichNode( gas1, rich1gas)) {
                
         CbmLitCircleMaterial *gas2 = new CbmLitCircleMaterial(*gas1);
   
         gas1->SetZ( 0.5*( (kapton1->GetZ() + kapton1->GetThickness()/2.) + (mirror->GetZ() - mirror->GetThickness()/2.) ) );
         gas1->SetThickness( (mirror->GetZ() - mirror->GetThickness()/2.) - (kapton1->GetZ() + kapton1->GetThickness()/2.));
   
         gas2->SetZ( 0.5*( (mirror->GetZ() + mirror->GetThickness()/2.) + (kapton2->GetZ() - kapton2->GetThickness()/2.) ) );
         gas2->SetThickness( (kapton2->GetZ() - kapton2->GetThickness()/2.) - (mirror->GetZ() + mirror->GetThickness()/2.) );
   
         
         fvMaterials.push_back(gas1);
         
         fvMaterials.push_back(gas2);
      }
         
    }  
    
   }
     
}

Int_t CbmLitEnvironment::ReadRichNode( CbmLitCircleMaterial *mat, CbmGeoNode *node)
{

  if ( !node ) return 1;   

  TString name = node->getName();
  TString Sname = node->getShapePointer()->GetName(); 
  CbmGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm
  CbmGeoVector centerV = node->getCenterPosition().getTranslation();
  TArrayD *P = node->getParameters();
  //Int_t NP = node->getShapePointer()->getNumParam();
  CbmGeoMedium *material = node->getMedium();
  material->calcRadiationLength();

  //tube.ID = node->getMCid();
  mat->SetRadLength(material->getRadiationLength());
  //tube.F = 1.;
  mat->SetZ(nodeV.Z()+centerV.Z());
  
  
        Int_t NofComponents = material->getNComponents();
      Double_t Zeff = 0.;
      Double_t Aeff = 0.;
      Double_t W = 0.;
      for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
         Double_t p[3];
         material->getComponent(iComp, p);
         W += p[2];        
      } 
      for (Int_t iComp = 0; iComp < NofComponents; iComp++) {
         Double_t p[3];
         material->getComponent(iComp, p);
         Aeff += p[0] * p[2];
         Zeff += p[1] * p[2];
      } 
      Aeff /= W;
      Zeff /= W;
           
      mat->SetDensity(material->getDensity());
      mat->SetAeff(Aeff);
      mat->SetZeff(Zeff);
            
      mat->SetName(Sname);
     
   if ( Sname=="TRAP" ) {
      mat->SetInnerRadius(0.);
      mat->SetOuterRadius(1000.);
      mat->SetXCenter(0.);
      mat->SetYCenter(0.);   
      mat->SetThickness(2. * P->At(0));
      cout << " shape is : " << Sname << endl;
      return 0;
   } else if ( Sname=="SPHE" ) {
      mat->SetInnerRadius(0.);
      mat->SetOuterRadius(1000.);
      mat->SetXCenter(0.);
      mat->SetYCenter(0.);   
      mat->SetZ(mat->GetZ() + 0.5 * (P->At(0) + P->At(1)));
      mat->SetThickness(P->At(1) - P->At(0));
      cout << " shape is : " << Sname << endl;
      return 0;
   } else {
      cout <<" unknown shape : "<<Sname<<endl;
      return 1;
   }

}

ClassImp(CbmLitEnvironment)