#include <TFile.h>
#include <TTree.h>
#include <Riostream.h>
#include <TMath.h>

void correlation_EEE(Double_t diffCut = 0.01)
{
//
// This macro correlates the events measured by two EEE telescopes according to their GPS time
// within a time window = diffCut
//
// Data are read from the two ROOT trees created for each telescope
//

//
// Open and read the configuration file
//
   ifstream config;
   config.open("C:\\root\\macros\\config_correlation_EEE.txt", ios::in);
   //Check the existence of the config file
   if(!config.is_open()){
                         cout << "Please check the config file!" << endl;
                         break;
   }
   
   TString tmp1, tmp2, tmp3, tmp4,tmp5,tmp6;
   config >> tmp1; // Read line#1 of the config file (telescope code 1)
   char *tel_code1 = new char[tmp1.Length() + 1];
   tel_code1 = tmp1.Data();
   config >> tmp2; // Read line#2 of the config file (telescope code 2)
   char *tel_code2 = new char[tmp2.Length() + 1];
   tel_code2 = tmp2.Data();
   config >> tmp3; // Read line#3 of the config file (date)
   char *date = new char[tmp3.Length() + 1];
   date = tmp3.Data();
   config >> tmp4; // Read line#4 of the config file (data path telescope#1)
   char *path1 = new char[tmp4.Length() + 1];
   path1 = tmp4.Data();
   config >> tmp5; // Read line#5 of the config file (data path telescope#2)
   char *path2 = new char[tmp5.Length() + 1];
   path2 = tmp5.Data();
   config >> tmp6; // Read line#6 of the config file (data path correlation tree)
   char *path = new char[tmp6.Length() + 1];
   path = tmp6.Data();
//
// Input files
//
   TFile *f1 = new TFile(Form("%s%s-%s.root",path1,tel_code1,date));
   TTree *t1 = (TTree*)f1->Get("chambers");
   TFile *f2 = new TFile(Form("%s%s-%s.root",path2,tel_code2,date));
   TTree *t2 = (TTree*)f2->Get("chambers");
//
// Define tree structure	
//
   Double_t rad=0.017453293;
   Double_t time1, time2;
   Double_t cx1, cx2, cy1, cy2, cz1, cz2,teta1,phi1,teta2,phi2,chi1,chi2,tof1,l1,tof2,l2;
   t1->SetBranchAddress("time", &time1);
   t1->SetBranchAddress("teta", &teta1);
   t1->SetBranchAddress("phi", &phi1);
   t1->SetBranchAddress("chi", &chi1);
   t1->SetBranchAddress("tof", &tof1);
   t1->SetBranchAddress("l", &l1);

   t2->SetBranchAddress("time", &time2);
   t2->SetBranchAddress("teta", &teta2);
   t2->SetBranchAddress("phi", &phi2);
   t2->SetBranchAddress("chi", &chi2);
   t2->SetBranchAddress("tof", &tof2);
   t2->SetBranchAddress("l", &l2);

   Int_t nent1 = t1->GetEntries();
   Int_t nent2 = t2->GetEntries();
//      
// Find time range
//
   Double_t t1min, t1max, t2min, t2max, range1, range2;
   t1->GetEntry(        0); t1min = time1;
   t1->GetEntry(nent1 - 1); t1max = time1;
   t2->GetEntry(        0); t2min = time2;
   t2->GetEntry(nent2 - 1); t2max = time2;
   range1 = t1max - t1min;
   range2 = t2max - t2min;
   //cout.precision(10);
   cout.setf(ios::fixed);
   cout <<"N.entry1 = "<< nent1<<"   N.entry2 = "<<nent2 << endl;
   cout << "Time range file 1: " << t1min << " --> " << t1max << ", range = " << range1 << endl;
   cout << "Time range file 2: " << t2min << " --> " << t2max << ", range = " << range2 << endl;

   Double_t tmin = TMath::Min(t1min, t2min);
   //Double_t tmax = TMath::Max(t1max, t2max);
   cout << "Common measure time interval = "<<(TMath::Min(t1max,t2max)-TMath::Max(t1min,t2min))<< " s"<<endl;
//
// Chain mesh: define starting time cell for both trees
//
   Int_t firstCell1 = (Int_t)((t1min - tmin) / diffCut);
   Int_t firstCell2 = (Int_t)((t2min - tmin) / diffCut);
   Int_t lastCell1  = (Int_t)((t1max - tmin) / diffCut);
   Int_t lastCell2  = (Int_t)((t2max - tmin) / diffCut);
   //cout << "Starting, ending cell for tree 1: " << firstCell1 << ", " << lastCell1 << endl;
   //cout << "Starting, ending cell for tree 2: " << firstCell2 << ", " << lastCell2 << endl;
//
// define complete cell range
//
   Int_t ncells = (Int_t)TMath::Max(lastCell1, lastCell2);
   cout << "#cells: " << ncells << endl;
   cout<<"Working..."<<endl;
//	
// define index collector for all cells
//
   TArrayI *cell = new TArrayI[ncells];
   for (Int_t i = 0; i < ncells; i++) cell[i].Set(0);
//	
// loop on TTree #2 and add each entry to corresponding cell
//
     Int_t cellIndex, size;
     for (Int_t i = 0; i < nent2; i++) {
	t2->GetEntry(i);
	cellIndex = (Int_t)((time2 - tmin) / diffCut);
	if (cellIndex >= 0 && cellIndex < ncells) {
	  size = cell[cellIndex].GetSize();
	  cell[cellIndex].Set(size+1);
	  cell[cellIndex][size] = i;
         }
     }
//	
// Define output correlation tree
//
   TFile *fileout = new TFile(Form("%s%s-%s-%s.root",path,tel_code1,tel_code2,date), "RECREATE");
   TTree *treeout = new TTree("tree", "Delta T");
   Int_t e1, e2;	
   Double_t diff,tetarel;
   treeout->Branch("time1", &time1, "time1/D");
   treeout->Branch("chi1", &chi1, "chi1/D");
   treeout->Branch("tof1", &tof1, "tof1/D");
   treeout->Branch("l1", &l1, "l1/D");
   treeout->Branch("teta1", &teta1, "teta1/D");
   treeout->Branch("phi1", &phi1, "phi1/D");
   treeout->Branch("time2", &time2, "time2/D");
   treeout->Branch("chi2", &chi2, "chi2/D");
   treeout->Branch("tof2", &tof2, "tof2/D");
   treeout->Branch("l2", &l2, "l2/D");
   treeout->Branch("teta2", &teta2, "teta2/D");
   treeout->Branch("phi2", &phi2, "phi2/D");
   treeout->Branch("diff", &diff, "diff/D");
   treeout->Branch("tetarel", &tetarel, "tetarel/D");


   for(e1 = 0; e1 < nent1; e1++) {
     if (!(e1 % 10000)) cout << "\rCorrelating entry #" << e1 << flush;
	t1->GetEntry(e1);
	cellIndex = (Int_t)((time1 - tmin) / diffCut);
	for (Int_t i = cellIndex - 1; i <= cellIndex + 1; i++) {
	  if (i < 0 || i >= ncells) continue;
	  for (Int_t j = 0; j < cell[i].GetSize(); j++) {
            size = cell[cellIndex].GetSize();
            e2 = cell[i].At(j);
	    t2->GetEntry(e2);      
	    diff=time1-time2;   
            tetarel=acos(cos(teta1*rad)*cos(teta2*rad)+sin(teta1*rad)*sin(teta2*rad)*cos(phi2*rad-phi1*rad))/rad;
            if(TMath::Abs(diff) <= diffCut) {    
	      treeout->Fill();
	    }
	  }
	}
      }
//
// Closing files
//
    f1->Close();
    f2->Close();
    cout << endl;
    fileout->cd();
    treeout->Write();
    fileout->Close();
    cout<<"Correlation tree completed"<<endl;
	
}