//  -------------------------------------------------------------------------
//
//   ----- General Macro for getting information about sim-events
//
//         Author: Sebastian Kupny <sebastian.kupny@uj.edu.pl
//
//         Last Update: 11/09/12
//
//         Comments:
//            Macro prints out information about events
//  -------------------------------------------------------------------------
//
//   Usage inside ROOT interpreter: 
//        1) root -l "Kratta_get_interested_0903.C(\"asyeos_sim.root\")"
//        2)[root] .L Kratta_get_interested_0903.C
//          [root] Kratta_get_interested_0903("asyeos_sim.root")
//                         
//
//  -------------------------------------------------------------------------

#include <fstream>
#include <iostream>
#include <string>

using namespace std;

void Kratta_get_interested_0903 ( TString A_sim_file ) {
  TString dir = getenv("VMCWORKDIR");
  TString r3bdir = dir + "/macros";

  TString r3b_geomdir = dir + "/geometry";
  gSystem->Setenv("GEOMPATH",r3b_geomdir.Data());

  TString r3b_confdir = dir + "gconfig";
  gSystem->Setenv("CONFIG_DIR",r3b_confdir.Data());

// Output files
  //TString OutFile = "r3bsim.root";
  TString InFile = A_sim_file;
  


  // In general, the following parts need not be touched
  // ========================================================================

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();

  gSystem->Load("libGenVector");

  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  gSystem->Load("libField");
  gSystem->Load("libGen");

  //----  Load R3B specific libraries ---------------------------------------
  gSystem->Load("libR3Bbase");
  gSystem->Load("libR3BGen");
  gSystem->Load("libR3BPassive");
  gSystem->Load("libR3BData");
  gSystem->Load("libR3BCal");
  gSystem->Load("libR3BCalo");
  gSystem->Load("libR3BDch");
  gSystem->Load("libR3BGfi");
  gSystem->Load("libR3BLand");
  gSystem->Load("libR3BmTof");
  gSystem->Load("libR3BTof");
  gSystem->Load("libR3BATof");
  gSystem->Load("libR3BTra");
  gSystem->Load("libR3BChimera");
  gSystem->Load("libELILuMon");
  gSystem->Load("libData");
  gSystem->Load("libKratta"); 
  
  

  
  string stlInFile(A_sim_file);
  
  Double_t fEkin = atof ( (stlInFile.substr(20, 3)).c_str() );
  Double_t fMom = atof ( (stlInFile.substr(28, 3)).c_str() );
  //This is very bad done: Strongly depeend on the input argument and calling method in script ///TODO: Find better sollution
  cout <<"First substr " <<  stlInFile.substr(20, 3) << endl;
  cout <<"Second substr " <<  stlInFile.substr(28, 3) << endl;  
  
  cout <<"conv First substr " <<  fEkin << endl;
  cout <<"conv  Second substr " <<  fMom << endl;  
  
  Double_t fEkinInGeV = fEkin/1000.;

  
  TFile *f = new TFile(  InFile );
  TTree *t4 = (TTree*)f->Get("cbmsim");
  
  Int_t nevent = t4->GetEntries();
  cout << nevent << endl;

  TClonesArray * pp;
  t4->SetBranchAddress( "KrakowPoint", &pp);

  TH2F * hist1 = new TH2F("hist1", "Kratta deltaE-E", 1000, -0.02, 0.25, 1000, 0, 0.05);
  TH1F * hist2 = new TH1F("hist2", "Kratta. Energy sum\ndeltaE + E", 1000, -0.05, 0.30 );
  hist1->SetXTitle("E [GeV]");
  hist1->SetYTitle("deltaE [GeV]");
  Double_t fDeltaE;
  Double_t fE;
  
  Int_t fDetId;
  Int_t fTrackID;
  
  Int_t indexCut6 = 0;
  Int_t indexCut7 = 0;
  Int_t indexCut8 = 0;  
  Int_t indexCut9 = 0;    
  Int_t indexCut10 = 0;
  Int_t indexCut11 = 0;
  Int_t indexCut12 = 0;
  Int_t indexCut13 = 0;
  Int_t indexCut14 = 0;
  Int_t indexCut15 = 0;
  Int_t indexCut16 = 0;
  Int_t indexCut17 = 0;
  Int_t indexCut18 = 0;
  Int_t indexCut19 = 0;
  Int_t indexCut20 = 0;
  
  Int_t iFiredOnlyFirst = 0;
  Int_t iFiredBoth = 0;
  Int_t iFiredNoOne = 0;
  Int_t iTotalEnergyIsMax = 0;
  Int_t iDeltaEEnergyIsMax = 0; 
  Int_t iTotalEnergyIsNotMax = 0;
  Int_t iTotalEnergyIsLargerThanMax = 0; 
  
  
  for (Int_t i=0; i<nevent; i++) {

    t4-> GetEntry(i);  

    //pp->Print();
    Int_t nPoints = pp-> GetEntries();
    
    fDeltaE = 0;
    fE = 0;
    
    for (Int_t pi = 0; pi < nPoints; pi++ )
    {
      KrattaPoint *krattaPt = (KrattaPoint*) pp->ConstructedAt(pi);
      cout << "Event "<< i << ", point= " << pi << ", Eloss= " <<  krattaPt->GetEnergyLoss();
      krattaPt->Print();
      cout << ", DetectorID = " << krattaPt->GetDetectorID()  << endl;
      fDetId = krattaPt->GetDetectorID();
      fTrackID = krattaPt->GetTrackID();
      cout << "TrackID=" << fTrackID << endl;
      
      if ( fTrackID == 0 ){ ///Only primary particle
        if ( fDetId == 6 ) fDeltaE =  krattaPt->GetEnergyLoss();
        if ( fDetId == 7 ) fE = krattaPt->GetEnergyLoss();
      }
    }
    //fDeltaE *= 1000.;
    //fE *= 1000.;
    
    cout << "Pair (" << fE << "," << fDeltaE << ")" << endl;
    cout << "Total energy: " << fE + fDeltaE << endl;
    hist1->Fill( fE, fDeltaE );
    hist2->Fill( fE + fDeltaE );
    

    if ( fE == 0 && fDeltaE > 0 )
    {
      iFiredOnlyFirst++;
    }
    if ( fE > 0 && fDeltaE > 0 ) ){
      iFiredBoth++;
    }
    if ( fE == 0 && fDeltaE  == 0 ){
      iFiredNoOne++;
    }
    if ( ( 0.99*fEkinInGeV  <= ( fE  + fDeltaE)) &&  ( ( fE  + fDeltaE) < 1.01* fEkinInGeV ))){
      iTotalEnergyIsMax++; 
    }
    
    if ( ( 0.99*fEkinInGeV  <=   fDeltaE ) &&  ( fDeltaE < 1.01* fEkinInGeV )){
      iDeltaEEnergyIsMax++; 
    }
    if ( ( fE  + fDeltaE) < 0.99*fEkinInGeV ){
      iTotalEnergyIsNotMax++;
    }
    if (  1.01* fEkinInGeV <= ( fE  + fDeltaE) ){
      iTotalEnergyIsLargerThanMax++;
    }
    ///Cuty
    
    if ( fE > 0.228 && (fDeltaE > 0.02 && fDeltaE < 0.035) )
    {
      cout <<" Cut1: EventID: " <<  i << endl;
    }
    if ( fE < 0.005 && (fDeltaE > 0.0225 && fDeltaE < 0.03) )
    {
      cout <<" Cut2: EventID: " <<  i << endl;
    }
    if ( fE < 0.005 && fDeltaE < 0.005 )
    {
      cout <<" Cut3: EventID: " <<  i << endl;
    }
    if ( ( fE > 0.18 && fE < 0.22 ) &&  fDeltaE < 0.01 )
    {
      cout <<" Cut4: EventID: " <<  i << endl;
    }
    if ( fE < 0.1  )
    {
      cout <<" Cut5: EventID: " <<  i << endl;
    }
    
    ///2012-08-31
    if ( ( 0.185 < fE  && fE < 0.225 ) && ( 0.0225 < fDeltaE && fDeltaE < 0.03) )
    {
      cout <<" Cut6: EventID: " <<  i << endl;
      indexCut6++;
    }
    if ( ( fE < 0.185 ) && ( 0.0225 < fDeltaE && fDeltaE < 0.03) )
    {
      cout <<" Cut7: EventID: " <<  i << endl;
      indexCut7++;
    }
    ///Case, when particle is stopped in first crystal
    if (  fE == 0 )
    {
      cout <<" Cut8 (fE == 0): EventID: " <<  i << endl;
      indexCut8++;
    }
    ///Case, when particle is stopped in first crystal
    if (  fDeltaE == 0 )
    {
      cout <<" Cut9 (fDeltaE == 0): EventID: " <<  i << endl;
      indexCut9++;
    }    

    ///Na E    
    if (  fE < 0.01 )
    {
      cout <<" Cut10 for fE: EventID: " <<  i << endl;
      indexCut10++;
    }    
    if ( 0.01 < fE && fE < 0.05 )
    {
      cout <<" Cut11 for fE: EventID: " <<  i << endl;
      indexCut11++;
    }
    if ( 0.05 < fE && fE < 0.1  )
    {
      cout <<" Cut12 for fE: EventID: " <<  i << endl;
      indexCut12++;
    }
    if ( 0.1 < fE && fE < 0.15  )
    {
      cout <<" Cut13 for fE: EventID: " <<  i << endl;
      indexCut13++;
    }
    if ( 0.15 < fE && fE < 0.185  )
    {
      cout <<" Cut14 for fE: EventID: " <<  i << endl;
      indexCut14++;
    }
    if ( 0.185  < fE && fE < 0.225   )
    {
      cout <<" Cut15 for fE: EventID: " <<  i << endl;
      indexCut15++;
    }
    if ( 0.225 < fE    )
    {
      cout <<" Cut16 for fE: EventID: " <<  i << endl;
      indexCut16++;
    }    
    ///Na Delta E
    if (  fDeltaE < 0.0225 )
    {
      cout <<" Cut17 for fDeltaE: EventID: " <<  i << endl;
      indexCut17++;
    }
    
    if (  0.0225 <= fDeltaE  && fDeltaE < 0.03  )
    {
      cout <<" Cut18 forfDeltaE: EventID: " <<  i << endl;
      indexCut18++;
    }
    
    if (   0.03 <= fDeltaE   )
    {
      cout <<" Cut19 for fDeltaE: EventID: " <<  i << endl;
      indexCut19++;
    }
    if (  fDeltaE < 0.005 &&  fE < 0.01 )
    {
      cout <<" Cut20 for corner: EventID: " <<  i << endl;
      indexCut20++;
    }                  
    
  }

  cout << "Cut 6 index " << indexCut6 << endl;
  cout << "Cut 7 index " << indexCut7 << endl;
  cout << "Cut 8 index " << indexCut8 << endl;
  cout << "Cut 9 index " << indexCut9 << endl;

  cout << "Cut10 index " << indexCut10 << endl;
  cout << "Cut11 index " << indexCut11 << endl;
  cout << "Cut12 index " << indexCut12 << endl;
  cout << "Cut13 index " << indexCut13 << endl;
  cout << "Cut14 index " << indexCut14 << endl;
  cout << "Cut15 index " << indexCut15 << endl;
  cout << "Cut16 index " << indexCut16 << endl;
  cout << endl;
  cout << "Cut17 index " << indexCut17 << endl;
  cout << "Cut18 index " << indexCut18 << endl;
  cout << "Cut19 index " << indexCut19 << endl;
  cout << endl;
  cout << "Cut20 index " << indexCut20 << endl;


  cout <<  "fEkinInGeV \t";
  cout <<  "fMom  \t";
  cout <<  "iFiredOnlyFirst  \t";
  cout <<  "iFiredBoth   \t";
  cout <<  "iFiredNoOne   \t";
  cout <<  "iDeltaEEnergyIsMax   \t";
  
  cout <<  "iTotalEnergyIsMax   \t";
  cout <<  "iTotalEnergyIsNotMax \t";
  cout <<  "iTotalEnergyIsLargerThanMax \t";
  cout <<  endl;

  
  cout <<  fEkinInGeV  << "\t";
  cout <<  fMom  << "\t";
  cout <<  iFiredOnlyFirst  << "\t";
  cout <<  iFiredBoth   << "\t";
  cout <<  iFiredNoOne   << "\t";
  cout <<  iDeltaEEnergyIsMax   << "\t";
  
  cout <<  iTotalEnergyIsMax   << "\t";
  cout <<  iTotalEnergyIsNotMax   << "\t";
  cout <<  iTotalEnergyIsLargerThanMax   << "\t";
  cout <<  endl;
  
  cout << "Energy limitation " << 0.99*fEkinInGeV  << " " <<  1.01 *fEkinInGeV<<  endl;


  ofstream logfile;
  logfile.open ("Sim_output.dat", ios::out | ios::app ); 
  if ( logfile.is_open())
  {
     logfile <<  "Analysed file = " << stlInFile  << endl;;
    logfile <<  "fEkinInGeV \t";
    logfile <<  "fMomInGeV  \t";
    logfile <<  "iFiredOnlyFirst  \t";
    logfile <<  "iFiredBoth   \t";
    logfile <<  "iFiredNoOne   \t";
    logfile <<  "iDeltaEEnergyIsMax   \t";
    
    logfile <<  "iTotalEnergyIsMax   \t";
    logfile <<  "iTotalEnergyIsNotMax \t";
    logfile <<  "iTotalEnergyIsLargerThanMax \t";
    logfile <<  endl;


    logfile <<  fEkinInGeV  << "\t";
    logfile <<  fMom/1000.  << "\t";
    logfile <<  iFiredOnlyFirst  << "\t";
    logfile <<  iFiredBoth   << "\t";
    logfile <<  iFiredNoOne   << "\t";
    logfile <<  iDeltaEEnergyIsMax   << "\t";
    
    logfile <<  iTotalEnergyIsMax   << "\t";
    logfile <<  iTotalEnergyIsNotMax   << "\t";
    logfile <<  iTotalEnergyIsLargerThanMax   << "\t";
    logfile <<  endl;

    logfile << "Energy limitation " << 0.99*fEkinInGeV  << " " <<  1.01 *fEkinInGeV<<  endl;
    logfile.close();
  }


  gStyle->SetPalette(1);
  hist1->Draw("COLZ");
  hist1->SaveAs( Form ("./outputPict/Kratta_hist_DeltaE-E_Ekin_%3.0f_Mom_%3.0f.root",fEkin, fMom ) );
  
  hist2->Draw();
  hist2->SaveAs( Form ("./outputPict/Kratta_hist_TotalEnergySum_DeltaEE_Ekin_%3.0f_Mom_%3.0f.root",fEkin, fMom ) );
  
  cout << "END" << endl;

}