#include <math.h>
#include <stdio.h>
#include <string.h>
#include <fstream>
#include <vector>
#include <TCanvas.h>
#include <TRandom.h>
#include <TF1Convolution.h>
#include <TF1.h>
#include <TH1F.h>
 
using namespace std;
 
void test()
{	 
	int i, j, k, k1, k2, Nruns, iRun, type[600], line;
	
    Float_t amp_direction1;
    Double_t tof_direction1;
    Int_t detn_direction1, BunchNumber_direction1, BunchNumber_direction_previous;
    
    Float_t amp_direction2;    
    Int_t detn_direction2;    			                
    Double_t tof_direction2;    
    
	int kclosest[16], kchosen;
	double closest_tofs[16], closest_amps[16];		
	double delta_tof_closest, delta_amp_closest;		
	
   	string condition1, condition2, condition3;	
   	int NBINS1, NBINS2;
   	double xmin1, xmax1;
   	double xmin2, xmax2;   		   	 		
	    	
   	std::vector<TH1D*> h1_vector_delta_tofs_B(32), h1_vector_delta_tofs_Cl(32), h1_vector_delta_tofs_AmOFF(32), h1_vector_delta_tofs_AmON(32);  
   	std::vector<TH1D*> h1_vector_delta_amps_B(32), h1_vector_delta_amps_Cl(32), h1_vector_delta_amps_AmOFF(32), h1_vector_delta_amps_AmON(32);     
   	
    std::vector<vector<TH1D*>> h1_matrix_delta_tofs_B(16, vector<TH1D*>(16)), h1_matrix_delta_tofs_Cl(16, vector<TH1D*>(16)), h1_matrix_delta_tofs_AmOFF(16, vector<TH1D*>(16)), h1_matrix_delta_tofs_AmON(16, vector<TH1D*>(16));  
    std::vector<vector<TH1D*>> h1_matrix_delta_amps_B(16, vector<TH1D*>(16)), h1_matrix_delta_amps_Cl(16, vector<TH1D*>(16)), h1_matrix_delta_amps_AmOFF(16, vector<TH1D*>(16)), h1_matrix_delta_amps_AmON(16, vector<TH1D*>(16));          		    		
		
    NBINS1=1e3, xmin1=-1e3, xmax1=1e3;  		
    NBINS2=2e5, xmin2=-1e5, xmax2=1e5;  	     
    for(i = 1; i <= 32; i++)
    {
        string position;
        if(i<=16) {position = "_V"; }
        else if(i>=17) {position = "_H"; }        
        string name1 = "h1_vector_delta_tofs_B_" + std::to_string(i) + position;             
        h1_vector_delta_tofs_B[i-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

        string name2 = "h1_vector_delta_amps_B_" + std::to_string(i) + position;              
        h1_vector_delta_amps_B[i-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2); 
        
        name1 = "h1_vector_delta_tofs_Cl_" + std::to_string(i) + position;             
        h1_vector_delta_tofs_Cl[i-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

        name2 = "h1_vector_delta_amps_Cl_" + std::to_string(i) + position;              
        h1_vector_delta_amps_Cl[i-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2); 
        
        name1 = "h1_vector_delta_tofs_AmOFF_" + std::to_string(i) + position;             
        h1_vector_delta_tofs_AmOFF[i-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

        name2 = "h1_vector_delta_amps_AmOFF_" + std::to_string(i) + position;              
        h1_vector_delta_amps_AmOFF[i-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2);   
        
        name1 = "h1_vector_delta_tofs_AmON_" + std::to_string(i) + position;             
        h1_vector_delta_tofs_AmON[i-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

        name2 = "h1_vector_delta_amps_AmON_" + std::to_string(i) + position;              
        h1_vector_delta_amps_AmON[i-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2);  
    }    
        
    for(i = 1; i <= 16; i++)
    {        
        for(j = 1; j <= 16; j++)
        {
            string name1 = "h1_matrix_delta_tofs_B_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";           
            h1_matrix_delta_tofs_B[i-1][j-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

            string name2 = "h1_matrix_delta_amps_B_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";            
            h1_matrix_delta_amps_B[i-1][j-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2); 
            
            name1 = "h1_matrix_delta_tofs_Cl_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";           
            h1_matrix_delta_tofs_Cl[i-1][j-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

            name2 = "h1_matrix_delta_amps_Cl_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";            
            h1_matrix_delta_amps_Cl[i-1][j-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2); 
            
            name1 = "h1_matrix_delta_tofs_AmOFF_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";           
            h1_matrix_delta_tofs_AmOFF[i-1][j-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

            name2 = "h1_matrix_delta_amps_AmOFF_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";            
            h1_matrix_delta_amps_AmOFF[i-1][j-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2);   
            
            name1 = "h1_matrix_delta_tofs_AmON_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";           
            h1_matrix_delta_tofs_AmON[i-1][j-1] = new TH1D(name1.c_str(), name1.c_str(), NBINS1, xmin1, xmax1);

            name2 = "h1_matrix_delta_amps_AmON_" + std::to_string(i) + "_V_" + std::to_string(j) + "_H";            
            h1_matrix_delta_amps_AmON[i-1][j-1] = new TH1D(name2.c_str(), name2.c_str(), NBINS2, xmin2, xmax2);         
        }                         
    } 		
				 
	for(i=408; i<=550; i++)
	{	        	
        string fileName = "./../ROOTfiles/ROOTfiles_original/run" + std::to_string(208000+i) + ".root";
	    ifstream file(fileName);
	    if (file.good()) 
	    {
	        //We read the Tree file.
		    std::cout << "File " << fileName << " does exist." << std::endl;
		    TFile *inputFile = new TFile(fileName.c_str(), "READ");			    
		    TTree *DSSSD_Tree = (TTree*) inputFile->Get("DSSD");
		    //We are only interested in the following branches, so we disable the rest not to be constantly copying them.
		    DSSSD_Tree->SetBranchStatus("*", 0);
            DSSSD_Tree->SetBranchStatus("detn", 1);
            DSSSD_Tree->SetBranchStatus("tof", 1); 
            DSSSD_Tree->SetBranchStatus("BunchNumber", 1); 
            DSSSD_Tree->SetBranchStatus("PSpulse", 1);   
            DSSSD_Tree->SetBranchStatus("amp", 1);  
            DSSSD_Tree->SetBranchStatus("area", 1);                                           
            gROOT->cd();                  			    
		    
		    //Trees for the vertical strips and a second one only to contain the specific events with associated bunch number.
           	TTree *Tree_vertical_strips;	
            TTree *Tree_horizontal_strips_whole;                	
            TTree *Tree_horizontal_strips;  
            
            //We create a Tree for vertical strips.             			    
            condition1 = "(detn >= 1) && (detn <= 16)";
		    Tree_vertical_strips = DSSSD_Tree->CopyTree(condition1.c_str());	
	        condition2 = "(detn >= 17) && (detn <= 32)";	        	 
            Tree_horizontal_strips_whole = DSSSD_Tree->CopyTree(condition2.c_str());  			    		        	 
	               
	        Tree_vertical_strips->SetBranchAddress("amp", &amp_direction1);		        
            Tree_vertical_strips->SetBranchAddress("detn", &detn_direction1);  		        
	        Tree_vertical_strips->SetBranchAddress("tof", &tof_direction1);            
	        Tree_vertical_strips->SetBranchAddress("BunchNumber", &BunchNumber_direction1);  		               
	               
	        //Now we look for coincidences of all the events in each strip.	        
	        for(k1=1; k1<=Tree_vertical_strips->GetEntries(); k1++)
	        {
	            Tree_vertical_strips->GetEntry(k1);
                if((k1/(Tree_vertical_strips->GetEntries()/10))*(Tree_vertical_strips->GetEntries()/10)==k1)
                {
                    cout << k1 << endl;
                }	   			            

                //We make this to avoid creating the same tree several times.
	            if((BunchNumber_direction1 != BunchNumber_direction_previous)||(k1 == 1))
	            {	  
		            condition3 = "(BunchNumber == " + std::to_string(BunchNumber_direction1) + ")";
	                Tree_horizontal_strips = Tree_horizontal_strips_whole->CopyTree(condition3.c_str());    

                    Tree_horizontal_strips->SetBranchAddress("detn", &detn_direction2);  			                    
                    Tree_horizontal_strips->SetBranchAddress("tof", &tof_direction2);   
                    Tree_horizontal_strips->SetBranchAddress("amp", &amp_direction2);  		                    		                    
                }	                
                
	            //Coincidences.
                for(k=1; k<=16; k++)
                { 
                    closest_tofs[k-1] = -1.0;
                }   		            
                for(k2=1; k2<=Tree_horizontal_strips->GetEntries(); k2++)
                {	    		            
                    Tree_horizontal_strips->GetEntry(k2);
                    if((abs(tof_direction1-tof_direction2)<abs(tof_direction1-closest_tofs[detn_direction2-17]))||(closest_tofs[detn_direction2-17] < 0))
                    {
                        closest_tofs[detn_direction2-17] = tof_direction2;					                                    
                        kclosest[detn_direction2-17] = k2;
                    }			                		                
                }                                      
                
                kchosen = 1;                    
                Tree_horizontal_strips->GetEntry(kclosest[kchosen-1]);
                delta_tof_closest = tof_direction1 - tof_direction2;                        
                for(k=1; k<=16; k++)
                {                        
                    Tree_horizontal_strips->GetEntry(kclosest[k-1]);

                    if( abs(tof_direction1-tof_direction2) <= abs(delta_tof_closest) )
                    {
                        kchosen = k;
                        delta_tof_closest = tof_direction1 - tof_direction2;
                        delta_amp_closest = amp_direction1 - amp_direction2;
                    }                     
                }		 
                   
                Tree_horizontal_strips->GetEntry(kclosest[kchosen-1]);  		                            
                BunchNumber_direction_previous = BunchNumber_direction1;                            	                
                
                h1_vector_delta_tofs_Cl[detn_direction1-1]->Fill(delta_tof_closest);
                h1_vector_delta_amps_Cl[detn_direction1-1]->Fill(delta_amp_closest);	
                
                h1_vector_delta_tofs_Cl[detn_direction2-1]->Fill(delta_tof_closest);
                h1_vector_delta_amps_Cl[detn_direction2-1]->Fill(delta_amp_closest);			                
                
                h1_matrix_delta_tofs_Cl[detn_direction1-1][detn_direction2-17]->Fill(delta_tof_closest);
                h1_matrix_delta_amps_Cl[detn_direction1-1][detn_direction2-17]->Fill(delta_amp_closest);	            
          
	        }		        	        		                             		        
	        DSSSD_Tree->ResetBranchAddresses();
		    delete DSSSD_Tree;	
		    inputFile->Close();  	
		    delete inputFile; 			    			    
		    
	        Tree_vertical_strips->ResetBranchAddresses();	
	        Tree_horizontal_strips->ResetBranchAddresses();	
	        Tree_horizontal_strips_whole->ResetBranchAddresses();			        		        		    
            delete Tree_vertical_strips;          
            delete Tree_horizontal_strips; 	            
            delete Tree_horizontal_strips_whole;  

            gDirectory->cd(); 
            gDirectory->Clear();    
            gDirectory->GetList()->Delete();                   
            gROOT->cd();      
            gROOT->GetListOfFiles()->Delete();                                                              
	    }			
	    file.close();	    			    		    	    			    
	}	    
	
	
	TFile* file = new TFile("./Intermediate_files/CC1_Distributions.root", "RECREATE");
    for(i = 0; i<32; i++) 
    {	
        h1_vector_delta_tofs_B[i]->Write();
        h1_vector_delta_tofs_AmOFF[i]->Write();
        h1_vector_delta_tofs_AmON[i]->Write();        
        h1_vector_delta_tofs_Cl[i]->Write(); 
        
        h1_vector_delta_amps_B[i]->Write();
        h1_vector_delta_amps_AmOFF[i]->Write();
        h1_vector_delta_amps_AmON[i]->Write();        
        h1_vector_delta_amps_Cl[i]->Write();                         
    }
    
    for(i = 0; i<16; i++) 
    {	
        for(j = 0; j<16; j++) 
        {
            h1_matrix_delta_tofs_B[i][j]->Write();
            h1_matrix_delta_tofs_AmOFF[i][j]->Write();
            h1_matrix_delta_tofs_AmON[i][j]->Write();        
            h1_matrix_delta_tofs_Cl[i][j]->Write(); 
            
            h1_matrix_delta_amps_B[i][j]->Write();
            h1_matrix_delta_amps_AmOFF[i][j]->Write();
            h1_matrix_delta_amps_AmON[i][j]->Write();        
            h1_matrix_delta_amps_Cl[i][j]->Write();           
        }                         
    }    
	file->Close();	
		
    
	return;
}