#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TMath.h>
#include <TLegend.h>
#include <TImage.h>
#include <TStopwatch.h>
#include <TH1F.h>
#include <TH3F.h>
#include <TRandom3.h>
#include <TTree.h>
#include <TClonesArray.h>
#include <TBranch.h>
#include <TArrayD.h>
#include "VertexSimulator.h"
#include "PMultiplicity.h"
#include "Direction.h"
#endif

void simulation(){
    

//___________________________________________________________________________________
// definizione delle variabili
    
  Int_t numEvent=500;
  Double_t mu=0;
  Double_t sigmaXY=0.0001;
  Double_t sigmaZ=0.053;
  Double_t mom=0.6e9;
  Double_t beta=0.99;
  Double_t zetaPart=1;
  // beam pipe
  Double_t rayPipe=0.03;
  Double_t thickPipe=0.0008;
  Double_t X0Pipe=0.3528;
  // first layer
  Double_t rayLayer1=0.04;
  Double_t thickL1=0.0002;
  Double_t X0L1=0.0936;
  // second layer
  Double_t rayLayer2=0.07;


//___________________________________________________________________________________


  Int_t output3=1;			//una variabile di controllo per l'output


//___________________________________________________________________________________
// apro il file di input e lo leggo

  TFile *inf = TFile::Open("../histInput.root");
  if( ! inf ){
  cerr << "azz... no file" << endl;
  return;
  }
  
  TH1F *hInput;
    
  hInput = (TH1F*) inf->Get( "hR" );
  if( ! "hR" ){
  cerr << "azz... no mult histo"<< endl;
  }
  
  typedef struct {Double_t xVert,yVert,zVert;
  Int_t numParticle;} VERTEXMULT;
  static VERTEXMULT vertexMult;
  
  //TArrayD* vertexMult=new TArrayD(4);							//ho provato a vedere cosa viene fuori salvando un vettore al posto di una serie di variabili
  
  // Apertura del file di output
  TFile hfile("../simulation.root","RECREATE");
  // dichiarazione del TTree
  TTree *tree = new TTree("T","TTree con 2 branches");
  
  TClonesArray *hits1 = new TClonesArray("Hits1",numEvent);
  TClonesArray &hitsAr1 = *hits1;
  TClonesArray *hits2 = new TClonesArray("Hits2",numEvent);
  TClonesArray &hitsAr2 = *hits2;
  
  tree->Branch("VertMult",&vertexMult,"X/D:Y:Z:mult/I");
  tree->Branch("Hits1",&hits1);
  tree->Branch("Hits2",&hits2);
  

//___________________________________________________________________________________
// booking dei puntatori
  gRandom=new TRandom3();
  TH3F *hVertex = new TH3F("hVertex","Vertex",50,-0.6,0.6,50,-0.2,0.2,50,-0.2,0.2);
  TH1F *hXVertex = new TH1F("hXVertex","XVertex",100000,-6*sigmaXY,6*sigmaXY);
  TH1F *hYVertex = new TH1F("hYVertex","YVertex",100000,-6*sigmaXY,6*sigmaXY);
  TH1F *hZVertex = new TH1F("hZVertex","ZVertex",100000,-6*sigmaZ,6*sigmaZ);
  TH1F *hMultiplicity = new TH1F("hMultiplicity","Multiplicity",10000,-1,101);
  TH3F *hPipeInteractions = new TH3F("hPipeInteractions","PipeInteractions",50,-0.6,0.6,50,-0.2,0.2,50,-0.2,0.2);
  TH3F *hLayer1Interactions = new TH3F("hLayer1Interactions","Layer1Interactions",50,-0.6,0.6,50,-0.2,0.2,50,-0.2,0.2);
  TH3F *hLayer2Interactions = new TH3F("hLayer2Interactions","Layer2Interactions",50,-0.6,0.6,50,-0.2,0.2,50,-0.2,0.2);  
  TH1F *hTheta = new TH1F("hTheta","Theta",314,0,3.14);
  TH1F *hTheta1 = new TH1F("hTheta1","Theta1",314,0,3.14);
  TH1F *hTheta2 = new TH1F("hTheta2","Theta2",314,0,3.14);

//___________________________________________________________________________________
// definisco le variabili del vertice
  Double_t x=0;
  Double_t y=0;
  Double_t z=0;


//___________________________________________________________________________________
//inizio la simulazione  
  for (int i=0;i<numEvent;i++){
    Int_t output1=i+1;
    Int_t output2=numEvent/100;
    if(output1%output2==0){
    cout<<"Simulation... "<<output3<<"%\r";
    cout.flush();
    output3++;
    }


//___________________________________________________________________________________
// simulazione dei vertici
      
      VertexSimulator::simuVertex(x,y,z,mu,sigmaXY,sigmaZ);
      hVertex->Fill(z,x,y);
      hXVertex->Fill(x);
      hYVertex->Fill(y);
      hZVertex->Fill(z);


//___________________________________________________________________________________
//creazione della molteplicita`
      
      //multiplicity fron known distribution
      /*Int_t mult;
      molt=PMultiplicity::GetMultiplicity(1);
      hMultiplicity->Fill(mult);hXVeertex->Fill(x);
      */
      
      //multiplicity from hist
      Int_t mult=PMultiplicity::GetMultiplicity(hInput);
      hMultiplicity->Fill(mult);
      
      //straight multiplicity
      /*Int_t mult=PMultiplicity::GetMultiplicity();
      hMultiplicity->Fill(mult);*/
      
      vertexMult.xVert=x;
      vertexMult.yVert=y;
      vertexMult.zVert=z;
      vertexMult.numParticle=mult;
      
      /*vertexMult->AddAt(x,0);
      vertexMult->AddAt(y,1);
      vertexMult->AddAt(z,2);
      vertexMult->AddAt(mult,3);*/
      
//___________________________________________________________________________________
//trasporto le particelle
      Double_t cd[3];
    for(int j=0;j<mult;j++){
    
      Direction* particle=new Direction();
      
      Double_t theta=0;
      Double_t phi=0;
      Double_t thp=0;
      Double_t php=0;
      Double_t xDir=x;
      Double_t yDir=y;
      Double_t zDir=z;
        
      particle->GetDirection(theta,phi);
      Direction::GetInteraction(xDir,yDir,zDir,rayPipe,theta,phi);
      hPipeInteractions->Fill(zDir,xDir,yDir);
      hTheta->Fill(theta);
      //cout<<"theta="<<theta<<endl;
      //theta=particle->GetTheta();
      particle->Scattering(thp, php, mom, beta, zetaPart, thickPipe, X0Pipe);
      particle->Rotate(theta,phi,thp,php,cd);
      Direction::GetInteraction(xDir,yDir,zDir,(rayLayer1-rayPipe),theta,phi);
      new(hitsAr1[j]) Direction(xDir,yDir,zDir);
      hLayer1Interactions->Fill(zDir,xDir,yDir);
      hTheta1->Fill(theta);
      //cout<<"theta1="<<theta<<endl;
      //theta=particle->GetTheta();
      particle->Scattering(thp, php, mom, beta, zetaPart, thickL1, X0L1);
      particle->Rotate(theta,phi,thp,php,cd);
      Direction::GetInteraction(xDir,yDir,zDir,(rayLayer2-rayLayer1),theta,phi);
      new(hitsAr2[j]) Direction(xDir,yDir,zDir);
      hLayer2Interactions->Fill(zDir,xDir,yDir);
      hTheta2->Fill(theta);
      //cout<<"theta2="<<theta<<endl<<endl;
      delete particle;
    }

    tree->Fill();
    hits1->Clear();
    hits2->Clear();

}
    
//___________________________________________________________________________________
// drawings

    /*TCanvas* c=new TCanvas("Test Canvas","Test");
    c->cd();
    hVertex->SetXTitle("Z (m)");
    hVertex->SetYTitle("X (m)");
    hVertex->SetZTitle("Y (m)");
    hVertex->Draw();
    //hXVertex->Draw();
    //hYVertex->Draw();
    //hZVertex->Draw();
    //hMultiplicity->Draw();
    hPipeInteractions->SetMarkerColor(46);
    hPipeInteractions->Draw("SAME");
    hLayer1Interactions->SetMarkerColor(30);
    hLayer1Interactions->Draw("SAME");
    hLayer2Interactions->SetMarkerColor(40);
    hLayer2Interactions->Draw("SAME");
      
    c->Divide(2,2);
    c->cd(1);
    hTheta->Draw();
    c->cd(2);
    hTheta1->Draw();
    c->cd(3);
    hTheta2->Draw();*/
    

    // Save all objects in this file
    hfile.Write();

    // Close the file. 
    hfile.Close();
    
    cout<<"\nDone."<<endl;
	  
}