/*****
*
* Neutron Detector Quick Look
* author  Marcelli-Malvezzi
* Version 1.00 - March 2006
*
* Not only new release of NeutronDetectorScan.c; all is changed and a new name is requested!
*
*
* Description: The aim of Neutron Detector QL software is to keep under control the time bahaviour of this detector. 
*              It creates three output canvases: the first and the second are relative to Background channels, 
*              while the third one is relative to trigger events ("events" channel) instead.
*              Each output relates to one PAMELA marshroot. See documentation for a more detailed description of the output. 
*
*
* Parameters:
* 	TSTring base          - the path to the root directory for the specific Pamela unpack session
*                               There is no default value, without this input the program will not run
*       TString outDir        - the path where to save the output image (Default = ./) 
*       TString format        - the format which will be used fo rsave the produced images (Default = "jpg") 
*       Float_t DeltaT        - the time interval in ms for ND records colection used for histograms ND_QL_1 plotting (Default = 1000 ms)
*       Float_t DeltaTevtime  - the time interval in minute for calculation of average ND count rate per minute, 
*                               see ND_QL_2 and ND_QL_3 plots (Default = 1 minute)
* Version 1.1 - June 2006
*
* Fixed bugs: for a large namber of events is not possible to initialize vectors, so all graphs have been converted in histograms
*
*/

#include <iostream>
#include <string>
#include <sstream>
#include <math.h>
#include "TStyle.h"
#include "TLatex.h"
#include "TFile.h"
#include "TList.h"
#include "TTree.h"
#include "TObjString.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TF1.h"
#include "TGaxis.h"
#include "TString.h"
#include "TPaveText.h"
#include "EventHeader.h"
#include "PscuHeader.h"
#include "TMultiGraph.h"
#include "physics/neutronDetector/NeutronEvent.h"
#include "physics/neutronDetector/NeutronRecord.h"


using namespace std;

void ND_QL(TString base, TString outDir, TString format, ULong_t DeltaT, ULong_t DeltaTevtime){ // DeltaT in ms, DeltaTevtime in minute
  
  //---------------- Variables initialization --------------------------//
  Int_t    tmpSize=0;
  Int_t    yTrig=0; 
  Int_t    yUpperBackground=0; 
  Int_t    yBottomBackground=0;
  ULong_t  lastime=0, firstime=0;
  Double_t scale= 1./DeltaTevtime;
  string title, title1;
  double obt;  
  Float_t     cut; 
  stringstream oss, noentries;
  stringstream oss1, oss2, oss3;

  //------- load root file --------------
  TFile  *file  = new TFile(base.Data()) ;

  TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
  filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();

  if ( outDir == "" ) outDir = ".";
  if (!file){
    printf("No such file in the directory has been found");
    return;
  }

  //------------------- Takes the tree of ND  ---------------------------//
  TTree *phystr = (TTree*)file->Get("Physics");
  TBranch *ndBr = phystr->GetBranch("Neutron");
  TBranch *headBr = phystr->GetBranch("Header");
  
  //PAMELA classes	
  pamela::neutron::NeutronEvent  *ne=0;
  pamela::neutron::NeutronRecord *nr=0; 
  pamela::EventHeader *eh=0;
  pamela::PscuHeader *ph=0;

  phystr->SetBranchAddress("Neutron", &ne);
  phystr->SetBranchAddress("Header", &eh);
  
  Long64_t nevents = phystr->GetEntries();
  const Int_t size = nevents;

  //--------------- output if there are 0 or <0  entries ---------------------------------------//
  if (nevents<=0){
    cout<<"WARNING: No Entries, RETURN"<<"\n";

    TCanvas *canvas4 = new TCanvas("No entries", "No entries ", 400, 200);
    canvas4->SetFillColor(10);
    canvas4->cd();

    TLatex *l = new TLatex();
    l->SetTextAlign(12);
    l->SetTextSize(0.15);
    l->SetTextColor(2);
    noentries.str("");
    noentries<< "ND_QL:";
    l->DrawLatex(0.05, 0.7, noentries.str().c_str()); 
    noentries.str("");
    noentries<< "No entries for this files";
    l->DrawLatex(0.05, 0.5, noentries.str().c_str()); 

    if (outDir == "./") {
      oss.str("");
      oss << filename.Data() << "ND_QL." << format.Data();
    } else {
      oss.str("");
      oss << outDir.Data() << filename.Data() << "_ND_QL." << format.Data();
    }
  
    canvas4->Update();
    canvas4->SaveAs(oss.str().c_str());
 
    return;
  }

  //------------ first and last  events ->  nbin ---------------------------------//
  phystr->GetEntry(0);
  ph = eh->GetPscuHeader();
  firstime = ph->GetOrbitalTime();  
  phystr->GetEntry(nevents-1);
  ph = eh->GetPscuHeader();
  lastime = ph->GetOrbitalTime();
  const ULong_t nint=((lastime-firstime)/(DeltaTevtime*60000));
  const ULong_t nint1=((lastime-firstime)/(DeltaT));
  const ULong_t nint2=lastime-firstime;
  int nbin=(int)nint;
  int nbin1=(int)nint1;
  int nbin2=(int)nint2;
  double obmin=firstime;
  double obmax=lastime; 

  //************************** HISTOGRAMS ***************************************************//
  //---------------------------- First histograms -----------------------------------------//
  oss.str("");
  oss <<"Bottom Background Neutrons with poissionian fit (DeltaT = " << DeltaT << " ms)";
  oss1.str("");
  oss1 <<"Upper Background Neutrons with poissionian fit (DeltaT = " << DeltaT << " ms)";
  TH1F *h1 = new TH1F (filename.Data(), oss.str().c_str(), nbin1,obmin,obmax);
  TH1F *h3 = new TH1F (filename.Data(), oss1.str().c_str(), nbin1,obmin,obmax);
  TH1F *h1bis = new TH1F (filename.Data(), oss.str().c_str(), 30, 0 ,30);
  TH1F *h3bis = new TH1F (filename.Data(), oss1.str().c_str(), 30, 0, 30);
  //---------------------------- Second histograms -----------------------------------------//
  title="";
  title=filename+":  Bottom Background & Upper Background: Time Evolution (DT="+DeltaTevtime+" min)";
  title1="";
  title1=filename+". Ratio: UpperBackground / BottomBackground: Time Evolution (DT="+DeltaTevtime+" min)";
  TH1F *histo1= new TH1F(title.c_str(),title.c_str(),nbin,obmin,obmax);
  TH1F *histo1bis= new TH1F(title1.c_str(),title1.c_str(),nbin,obmin,obmax);
  TH1F *histo2= new TH1F(title.c_str(),title.c_str(),nbin,obmin,obmax);
  //---------------------------- Third histograms -----------------------------------------//

  oss.str("");
  oss << filename.Data() << ":  " << "Trigger events: Time Evolution"; 
  oss1.str("");
  oss1 <<filename.Data()  << ":  " << "Trigger counting: Time Evolution (DT= "<<DeltaTevtime<<" min)";
  oss2.str("");
  oss2 << filename.Data() << ":  " << "Trigger counting (events with more than 10 n): Time Evolution  (DT= "<<DeltaTevtime<<" min)";

  TH1F *Trig = new TH1F(oss.str().c_str(), oss.str().c_str(),nbin2,obmin,obmax);
  TH1F *Trigger = new TH1F(oss1.str().c_str(), oss1.str().c_str(),nbin,obmin,obmax);
  TH1F *Triggercut = new TH1F(oss2.str().c_str(), oss2.str().c_str(),nbin,obmin,obmax); 

  //-------------------------- Fill histograms from root-ple -----------------------------//
  for (Int_t i = 0; i < nevents; i++){
    ndBr->GetEntry(i);
    headBr->GetEntry(i);
    tmpSize = ne->Records->GetEntries();
    if (ne->unpackError == 1) continue;    //Check if NeutronEvent is corrupted   
    ph = eh->GetPscuHeader();
    obt= ph->GetOrbitalTime();
    for (Int_t j = 0; j < tmpSize; j++){
      nr = (pamela::neutron::NeutronRecord*)ne->Records->At(j);
      yTrig             = yTrig            + (int)nr->trigPhysics;
      yUpperBackground  = yUpperBackground  + (int)nr->upperBack;
      yBottomBackground = yBottomBackground + (int)nr->bottomBack;
    }
    histo1->Fill(ph->GetOrbitalTime(), yUpperBackground);
    histo1bis->Fill(ph->GetOrbitalTime(), yUpperBackground);
    histo2->Fill(ph->GetOrbitalTime(), yBottomBackground);
    h1->Fill(ph->GetOrbitalTime(), yBottomBackground);
    h3->Fill(ph->GetOrbitalTime(), yUpperBackground);
    Trig->Fill(ph->GetOrbitalTime(), yTrig);
    Trigger->Fill(ph->GetOrbitalTime(), yTrig);
    if(yTrig >=10)    
      Triggercut->Fill(ph->GetOrbitalTime(), yTrig);
    yUpperBackground=0;
    yBottomBackground=0;
    yTrig=0;
  }

  for (Int_t i = 0; i < nevents; i++){ 
    h1bis->Fill(h1->GetBinContent(i));
    h3bis->Fill(h3->GetBinContent(i));
  }  

  //********************************* CANVASES ************************************************//
  //--------------------------------- First canvas --------------------------------------------// 
  TCanvas *histoCanv = new TCanvas("ND_QL_1", base.Data(), 1280, 1024);
  histoCanv->SetFillColor(10);
  histoCanv->Divide(1,2);

  histoCanv->cd(1);
  h1bis->SetLineColor(kRed);
  h1bis->SetFillStyle(3004);
  h1bis->SetFillColor(kRed);
  h1bis->GetXaxis()->SetTitle("Number of neutrons");
  h1bis->GetXaxis()->CenterTitle();
  h1bis->GetYaxis()->SetTitle("Number of events");
  h1bis->GetYaxis()->CenterTitle();
  gStyle->SetOptFit(1111);
  TF1 *h2 = new TF1 ("h2", "[0]*TMath::PoissonI(x, [1])" , -1., 30);
  h2->SetParName(0, "NEvents");
  h2->SetParName(1, "meanFreq");
  h2->SetParameter(0, size/2);
  h2->SetParameter(1, 1);
  h2->SetLineColor(kRed);
  h1bis->Fit("h2");
  h1bis->Draw();

  histoCanv->cd(2);
  h3bis->SetLineColor(kBlue);
  h3bis->SetFillStyle(3004);
  h3bis->SetFillColor(kBlue);
  h3bis->GetXaxis()->SetTitle("Number of neutrons");
  h3bis->GetXaxis()->CenterTitle();
  h3bis->GetYaxis()->SetTitle("Number of events");
  h3bis->GetYaxis()->CenterTitle();
  gStyle->SetOptFit(1111);
  TF1 *h4 = new TF1 ("h4", "[0]*TMath::PoissonI(x, [1])" , -1., 30);
  h4->SetParName(0, "NEvents");
  h4->SetParName(1, "meanFreq");
  h4->SetParameter(0, size/2);
  h4->SetParameter(1, 1);
  h4->SetLineColor(kBlue);
  h3bis->Fit("h4");
  h3bis->Draw();
  
  //---------------------------------- Second Canvas -----------------------------------------//
  TCanvas *finalCanv = new TCanvas("ND_QL_2", base.Data(), 1280, 1024);
  finalCanv->SetFillColor(10);
  finalCanv->Divide(1,2);

  finalCanv->cd(1);
  histo1->SetStats(kFALSE);
  histo1->GetYaxis()->SetTitle("Number of recorded Neutrons / DT   (min^-1) ");
  histo1->GetYaxis()->SetTitleSize(0.04);
  histo1->GetYaxis()->SetTitleOffset(0.5);
  histo1->GetYaxis()->CenterTitle(); 
  histo1->GetYaxis()->SetLabelSize(0.03);
  histo1->GetXaxis()->CenterTitle();
  histo1->GetXaxis()->SetTitleSize(0.04);
  histo1->GetXaxis()->SetTitleOffset(1);
  histo1->GetXaxis()->SetLabelSize(0.03);
  histo1->GetXaxis()->SetTitle("OBT   (ms)");
  histo1->SetMarkerSize(.5);
  histo1->SetMarkerStyle(21);
  histo1->SetMarkerColor(3);
  histo1->Scale(scale);
  histo1->SetMinimum(0.8);
  histo1->Draw("9p0");

  histo2->SetStats(kFALSE);
  histo2->SetMarkerSize(.5);
  histo2->SetMarkerStyle(21);
  histo2->SetMarkerColor(2);
  histo2->Scale(scale);
  histo2->SetMinimum(0.8);
  histo2->Draw("9p0same");

  TLegend *leg = new TLegend(.91,.65, .99, .81);
  leg->AddEntry(histo2,"Bottom");
  leg->AddEntry(histo1,"Upper");
  leg->Draw();

  finalCanv->cd(2);
  histo1bis->SetMarkerSize(.5);
  histo1bis->SetStats(kFALSE);
  histo1bis->SetMarkerStyle(21);
  histo1bis->SetMarkerColor(4);
  histo1bis->GetXaxis()->SetTitle("OBT   (ms)");
  histo1bis->GetXaxis()->SetTitleSize(0.04);
  histo1bis->GetXaxis()->CenterTitle();
  histo1bis->GetXaxis()->SetLabelSize(0.03);
  histo1bis->GetYaxis()->SetTitle("Upper / Bottom Background");
  histo1bis->GetYaxis()->CenterTitle();
  histo1bis->GetYaxis()->SetTitleSize(0.04);
  histo1bis->GetYaxis()->SetTitleOffset(0.5);
  histo1bis->GetYaxis()->SetLabelSize(0.03);
  histo1bis->Scale(scale);
  histo1bis->Divide(histo2);
  histo1bis->SetMinimum(0.);
  if(histo1bis->GetMaximum()<2)
    histo1bis->SetMaximum(2.0);
  histo1bis->Draw("9p0");
  TF1 *func1 = new TF1("func1", "1.4");
  func1->SetRange(firstime, lastime);
  func1->SetLineColor(3);
  func1->SetLineStyle(1);
  func1->SetLineWidth(3);
  func1->Draw("same");
  TF1 *func2 = new TF1("func2", "0.6");
  func2->SetRange(firstime, lastime);
  func2->SetLineColor(3);
  func2->SetLineStyle(1);
  func2->SetLineWidth(3);
  func2->Draw("same");
  
  //---------------------------------- Third Canvas -----------------------------------------//
  TCanvas *triggerCanv = new TCanvas("ND_QL_3", base.Data(), 1280, 1024);
  triggerCanv->SetFillColor(10);
  triggerCanv->Divide(1,3);

  triggerCanv->cd(1);
  Trig->SetStats(kFALSE);
  Trig->SetMarkerStyle(21);
  Trig->SetMarkerSize(.7);
  Trig->SetMarkerColor(2);
  Trig->GetXaxis()->SetTitle("OBT (ms)");
  Trig->GetXaxis()->CenterTitle();
  Trig->GetXaxis()->SetTitleSize(0.05);
  Trig->GetXaxis()->SetLabelSize(0.05);
  Trig->GetYaxis()->SetTitle("Number of Neutrons");
  Trig->GetYaxis()->CenterTitle();
  Trig->GetYaxis()->SetTitleSize(0.06);
  Trig->GetYaxis()->SetTitleOffset(0.5);
  Trig->GetYaxis()->SetLabelSize(0.05);
  Trig->Draw("9p0");

  triggerCanv->cd(2);
  Trigger->SetStats(kFALSE);
  Trigger->SetMarkerStyle(21);
  Trigger->SetMarkerSize(.7);
  Trigger->SetMarkerColor(4);
  Trigger->GetYaxis()->SetTitle("Number of Neutrons");
  Trigger->GetYaxis()->SetTitleSize(0.06);
  Trigger->GetYaxis()->SetTitleOffset(0.5);
  Trigger->GetYaxis()->CenterTitle(); 
  Trigger->GetYaxis()->SetLabelSize(0.05);
  Trigger->GetXaxis()->CenterTitle();
  Trigger->GetXaxis()->SetTitleSize(0.05);
  Trigger->GetXaxis()->SetLabelSize(0.05);
  Trigger->GetXaxis()->SetTitle("OBT (ms)");
  Trigger->SetMinimum(0.);
  Trigger->Draw("9p0");

  triggerCanv->cd(3);
  Triggercut->SetStats(kFALSE);
  Triggercut->SetMarkerStyle(21); 
  Triggercut->SetMarkerSize(.7);
  Triggercut->SetMarkerColor(4);
  Triggercut->GetYaxis()->SetTitle("Number of Neutrons");
  Triggercut->GetYaxis()->SetTitleSize(0.06);
  Triggercut->GetYaxis()->SetTitleOffset(.5);
  Triggercut->GetYaxis()->CenterTitle(); 
  Triggercut->GetYaxis()->SetLabelSize(0.05);
  Triggercut->GetXaxis()->CenterTitle();
  Triggercut->GetXaxis()->SetTitleSize(0.05);
  Triggercut->GetXaxis()->SetTitleOffset(1);
  Triggercut->GetXaxis()->SetLabelSize(0.05);
  Triggercut->GetXaxis()->SetTitle("OBT (ms)");
  Triggercut->Draw("9p");
 
  //************************** to save canvas ******************************************//

 if (outDir == "") {
    oss1.str("");
    oss2.str(""); 
    oss3.str("");
    oss1 << filename.Data() << "." << format.Data();
    oss2 << filename.Data() << "." << format.Data(); 
    oss3 << filename.Data() << "." << format.Data();
  } else {
    oss1.str("");
    oss2.str("");
    oss3.str("");
    oss1 << outDir.Data() << filename.Data() << "_ND_QL_1." << format.Data();
    oss2 << outDir.Data() << filename.Data() << "_ND_QL_2." << format.Data(); 
    oss3 << outDir.Data() << filename.Data() << "_ND_QL_3." << format.Data();
    }

  finalCanv->SaveAs(oss2.str().c_str()); 
  histoCanv->SaveAs(oss1.str().c_str());
  triggerCanv->SaveAs(oss3.str().c_str());

  file->Close();
 
}

int main(int argc, char* argv[]){
  TString path;
  TString outDir        ="./";
  TString format        ="jpg";
  ULong_t DeltaT        =1000;
  ULong_t DeltaTevtime  =5;
    
 if (argc < 2){
    printf("You have to insert at least the file to analyze \n");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  
  if (!strcmp(argv[1], "--help")){
        printf( "Usage: ND_QL FILE [OPTION] \n");
        printf( "\t --help                  Print this help and exit \n");	
        printf( "\t -outDir[path]           Path where to put the output [default ./] \n");
        printf( "\t -format[jpg|gif|ps]     Format for output files [default 'jpg'] \n");
	printf( "\t -DeltaT                 Time interval to bin histograms in ms [default = 1000ms] \n");
	printf( "\t -DeltaTevtime           Time interval to control time evolution of acquired data in minute [default = 1 min.] \n");
        exit(1);
  }
  path=argv[1];
   for (int i = 2; i < argc; i++){
     if (!strcmp(argv[i], "-outDir")){
       if (++i >= argc){
	 printf( "-outDir needs arguments. \n");
	 printf( "Try '--help' for more information. \n");
	 exit(1);
	} 
       else{
	 outDir = argv[i];
	 continue;
       }
     }  
     if (!strcmp(argv[i], "-format")){
       if (++i >= argc){
	 printf( "-format needs arguments. \n");
	 printf( "Try '--help' for more information. \n");
	 exit(1);
       } 
       else{
	 format = argv[i];
	 continue;
       }
     }
     if (!strcmp(argv[i], "-DeltaT")){
       if (++i >= argc){
	 printf( "-DeltaT needs arguments. \n");
	 printf( "Try '--help' for more information. \n");
	 exit(1);
       } 
       else{
	 DeltaT = atol(argv[i]);
	 continue;
       }
     } 
     if (!strcmp(argv[i], "-DeltaTevtime")){
       if (++i >= argc){
	 printf( "-DeltaTevtime needs arguments. \n");
	 printf( "Try '--help' for more information. \n");
	 exit(1);
       }
       else{
	 DeltaTevtime = atol(argv[i]);
	 continue;
       }
     }
   }
   
   ND_QL(argv[1], outDir, format, DeltaT, DeltaTevtime);
}