#include "TKratTreat.h"
#include "TClonesArray.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"


// includes for modeling
#include "TGeoManager.h"
#include "TParticle.h"
#include "TVirtualMC.h"
#include "TGeoMatrix.h"
#include "TGeoMaterial.h"
#include "TGeoMedium.h"
#include "TGeoBBox.h"
#include "TGeoCompositeShape.h"
#include "TGeoShapeAssembly.h"


#include "TVector3.h"
#include "TMath.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TH2F.h"
#include <string>
#include <iostream>


#include "KrattaPoint.h"
#include "R3BMCTrack.h"		

using std::cout;
using std::endl;

		

TKratTreat::TKratTreat() :
  FairTask("R3B Land Digitization scheme ")
{
  fVerbose = 1;
}


TKratTreat::TKratTreat(Int_t verbose) :
  FairTask("R3B Land Digitization scheme ", verbose)
{
  fVerbose          = verbose;
  fKrakowPointName  = "KrakowPoint";
  fMCTrackName      = "MCTrack";
  fFSName           = "FASTSLOW";
  fSaveOutputToTree = kFALSE;
}


TKratTreat::~TKratTreat() {
}


void TKratTreat::SetParContainers() {

  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  if ( ! rtdb ) Fatal("SetParContainers", "No runtime database");

/*pc  fLandDigiPar = (R3BLandDigiPar*)(rtdb->getContainer("R3BLandDigiPar"));

  if ( fVerbose && fLandDigiPar ) {
      cout << "-I- TKratTreat::SetParContainers() "<< endl;
      cout << "-I- Container R3BLandDigiPar  loaded " << endl;
  }
pc*/
}




InitStatus TKratTreat::Init() {

  // Get input array 
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No FairRootManager");
  fKrattaPoints = (TClonesArray*) ioman->GetObject( fKrakowPointName );
  fKrattaMCTrack = (TClonesArray*) ioman->GetObject( fMCTrackName );
  cout << "linea 82 di TKratTreat::Init()" << endl;
  if ( ! gGeoManager )cout << "Init " << " No GeoManager"<<endl;
  getchar();

  hLight =  new TH1F("Light","Light",100,0.,500);
  hEnergy = new TH1F("Energy","Energy",100,0.,500);
  hLightEnergy =  new TH2F("LightEnergy","LightEnergy",500,0.,250,500,0.,250.);
  hFS =  new TH2F("FS","FS",250,0.,250,250,0.,250.);
  
  hLight2 =  new TH1F("Light2","Light2",100,0.,1000);
  hEnergy2 = new TH1F("Energy2","Energy2",100,0.,1000);
  hLightEnergy2 =  new TH2F("LightEnergy2","LightEnergy2",500,0.,1000,500,0.,1000.);
  hFS2 =  new TH2F("FS2","FS2",500,0.,1000,500,0.,1000.);

  hLight3 =  new TH1F("Light3","Light3",100,0.,1000);
  hEnergy3 = new TH1F("Energy3","Energy3",100,0.,1000);
  hLightEnergy3 =  new TH2F("LightEnergy3","LightEnergy3",500,0.,1000,500,0.,1000.);
  hFS3 =  new TH2F("FS3","FS3",500,0.,1000,500,0.,1000.);

  hLight4 =  new TH1F("Light4","Light4",100,0.,500);
  hEnergy4 = new TH1F("Energy4","Energy4",100,0.,500);
  hLightEnergy4 =  new TH2F("LightEnergy4","LightEnergy3",500,0.,500,500,0.,500.);
  hFS4 =  new TH2F("FS4","FS4",500,0.,500,500,0.,500.);

  hLight5 =  new TH1F("Light5","Light5",100,0.,500);
  hEnergy5 = new TH1F("Energy5","Energy5",100,0.,500);
  hLightEnergy5 =  new TH2F("LightEnergy5","LightEnergy5",500,0.,500,500,0.,500.);
  hFS5 =  new TH2F("FS5","FS5",500,0.,500,500,0.,500.);

  hFSTot =  new TH2F("FSTot","FSTot",1000,0.,2000,1000,0.,2000.);

  fs =new FS();
  ioman->Register( fFSName, "Sim Fast Slow", fs, fSaveOutputToTree );



  return kSUCCESS;

}


// -------------------------------------------------------------------------

// -----   Public method Exec   --------------------------------------------
void TKratTreat::Exec(Option_t* opt)
{
   cout<<"new event..............................."<<endl;
   
   double de[352],dlight[352],dslow[352];
   for (int i=0;i<352;i++){de[i]=0;dlight[i]=0;dslow[i]=0;} 
   Int_t nentries = fKrattaPoints->GetEntries();
   cout<<"   entries =" << nentries <<endl;
   
   
   Int_t TrackId=0;
   Double_t ratio=0;
   
   for (Int_t l=0;l<nentries;l++){

//       cout<<"loop over entries "<<l<<endl;
      // Get the Land Object in array
      KrattaPoint *fKratta_obj = (KrattaPoint*) fKrattaPoints->At(l);

//      Int_t paddle = int(land_obj->GetSector())-1; //note that paddle starts at 1 
//      Int_t scint = int(land_obj->GetPaddleNb());
      Double_t eloss = fKratta_obj->GetEnergyLoss()*1000.;  
      Double_t light = fKratta_obj->GetLightYield()*1000.;  
      Double_t slow = fKratta_obj->GetSlow()*1000.;  
      Double_t xIn = fKratta_obj->GetXIn();
      Double_t yIn = fKratta_obj->GetYIn();
      Double_t zIn = fKratta_obj->GetZIn();
      Double_t xOut = fKratta_obj->GetXOut();
      Double_t yOut = fKratta_obj->GetYOut();
      Double_t zOut = fKratta_obj->GetZOut();
      Double_t time = fKratta_obj->GetTime();
      Int_t detID = fKratta_obj->GetDetectorID();
      
//      Int_t media = int(land_obj->GetPaddleType());
      TrackId = fKratta_obj->GetTrackID();
      R3BMCTrack *aTrack = (R3BMCTrack*) fKrattaMCTrack->At(TrackId);      
      Int_t PID = aTrack->GetPdgCode();
      Int_t mother = aTrack->GetMotherId();
//      cout <<" TrackId " << TrackId << " PID = "<< PID << " mother = " << mother << 
//             " ELoss " << eloss << " light " << light <<" slow " << slow <<" detid = "<< detID;
      if(!gGeoManager) cout << " gGeoManager missing" << endl;	     
      gGeoManager->FindNode(xIn,yIn,zIn);
      TGeoNode* curNode = gGeoManager->GetCurrentNode();
      Int_t ntel=curNode->GetNumber();
      
      cout  << " " <<  curNode->GetNumber() <<endl;
      
      if(ntel<=495){
       de[ntel]=de[ntel]+eloss;
       Double_t lR=gRandom->Gaus(light,light*0.001/2.35);
       dlight[ntel]=dlight[ntel]+lR;
       Double_t sR=gRandom->Gaus(slow,slow*0.001/2.35);
       dslow[ntel]=dslow[ntel]+sR;
       
       
//       cout  << " " <<  ntel <<" "<<de[ntel-144] <<" "<<dlight[ntel-144]<<endl;
       if(PID>1000020040)ratio=0.8;
       if(PID==1000020040)ratio=0.8;
       if(PID==1000020030)ratio=0.8;
       if(PID==1000010030)ratio=0.8;
       if(PID==1000010020)ratio=0.8;
       if(PID==2212)ratio=0.8;
       if(PID==22)ratio=0.8;
//      dslow[ntel] =dslow[ntel] +5*(light-eloss*ratio);
//      cout << gMC->TrackCharge()<<endl;
       }
   }
   
   
   
 hLight->Fill(dlight[30]);
 hEnergy->Fill(de[30]);
 hLightEnergy->Fill(de[30],dlight[30],1.);
 hFS->Fill(dslow[30],dlight[30],1.);

 hLight2->Fill(dlight[35]);
 hEnergy2->Fill(de[35]);
 hLightEnergy2->Fill(de[35],dlight[35],1.);
 hFS2->Fill(dslow[35],dlight[35],1.);

 hLight3->Fill(dlight[0]);
 hEnergy3->Fill(de[0]);
 hLightEnergy3->Fill(de[0],dlight[0],1.);
 hFS3->Fill(dslow[0],dlight[0],1.);

 hLight4->Fill(dlight[5]);
 hEnergy4->Fill(de[5]);
 hLightEnergy4->Fill(de[5],dlight[5],1.);
 hFS4->Fill(dslow[5],dlight[5],1.);

 hLight5->Fill(dlight[10]);
 hEnergy5->Fill(de[10]);
 hLightEnergy5->Fill(de[10],dlight[10],1.);
 hFS5->Fill(dslow[10],dlight[10],1.);
 int ifired=0;
 for(int i=0;i<352;i++){
//  if(dslow[i]>800)dslow[i]=800;
//  if(dlight[i]>1100)dlight[i]=1100;
  if(dslow[i]>0 && dlight[i]>0){
   hFSTot->Fill(dslow[i],dlight[i],1.);
   fs->fast[ifired]=dlight[i];
   fs->slow[ifired]=dslow[i];
   fs->de[ifired]=de[i];
   fs->numtel[ifired]=i;
   ifired++;
  }
 }
 fs->multi=ifired;
 
}
// -------------------------------------------------------------------------

void TKratTreat::Reset(){
// Clear the structure
//   cout << " -I- Digit Reset() called " << endl;

   
//pc if (fLandDigi ) fLandDigi->Clear();
}   

void TKratTreat::Finish()
{
   hLight->Write();
   hEnergy->Write();
   hLightEnergy->Write();
   hFS->Write();
   hLight2->Write();
   hEnergy2->Write();
   hLightEnergy2->Write();
   hFS2->Write();
   hLight3->Write();
   hEnergy3->Write();
   hLightEnergy3->Write();
   hFS3->Write();
   hLight4->Write();
   hEnergy4->Write();
   hLightEnergy4->Write();
   hFS4->Write();
   hLight5->Write();
   hEnergy5->Write();
   hLightEnergy5->Write();
   hFS5->Write();hFSTot->Write();
// here event. write histos
//   cout << " -I- Digit Finish() called " << endl;
// Write control histograms
/*pc   hPMl->Write();
   hPMr->Write();
   hMult1->Write();
   hMult2->Write();
   hParticle->Write();
   hFirstEnergy->Write();
   hPaddleEnergy->Write();
   
   hTotalLight->Write();
   hTotalLight1->Write();
   hTotalEnergy->Write();
   hTotalEnergy1->Write();
   hTotalLightRel->Write();
   hTotalLightRel2->Write();
   hTotalEnergyRel->Write();

   hElossLight->Write();
pc*/
}

/*pcR3BLandDigi* TKratTreat::AddHit(Int_t paddleNr, Double_t tdcL, Double_t tdcR,
   Double_t tdc, Double_t qdcL, Double_t qdcR, Double_t qdc, 
   Double_t xx, Double_t yy, Double_t zz){   
  // It fills the R3BLandDigi array
  R3BLandDigi *digi =
    new((*fLandDigi)[fLandDigi->GetEntriesFast()]) R3BLandDigi(paddleNr,tdcL,tdcR,tdc,qdcL,qdcR,qdc,xx,yy,zz);
  return digi;
}pc*/
	      

ClassImp(TKratTreat)