diff -pruN orig/caps.cpp new/caps.cpp
--- orig/caps.cpp	2012-12-15 17:13:23.000000000 +0200
+++ new/caps.cpp	2020-09-09 23:07:46.080566000 +0300
@@ -14,7 +14,7 @@
 #include <gsl/gsl_statistics.h>
 #include<sys/time.h>
 #include<iomanip>
-
+#include <bits/stdc++.h> 
 
 
 
@@ -69,6 +69,8 @@
 const gsl_rng_type * T;
 gsl_rng *r;
 
+vector<double> totaltempnew;
+double alphathresh = 0;
 int main(int argc, char *argv[]){
 
 
@@ -543,16 +545,27 @@ int main(int argc, char *argv[]){
 
 
 					print_splash(output);
+					OUTPUT << "\n\File1: " << files[i] << endl;
 					vec1.print_to_fasta(output.c_str());
+					OUTPUT << "\n\nFile2: " << files[j] << endl;
 					vec2.print_to_fasta(output.c_str());
 					int length1 = vec1.sequences[0].length();			
 					int length2 = vec2.sequences[0].length();
 
+					OUTPUT << "\n\nLength1: " << length1 << endl;
+					OUTPUT << "Length2: " << length2 << endl;
 
 
 					if(tree_in ==0){
 						tree1 = create_input_tree(vec1.names, vec1.sequences);
 						tree2 = create_input_tree(vec2.names, vec2.sequences);
+						
+						// Output the CAPS generated trees to the .out file of each pair
+						string temptre1 = TreeTemplateTools::treeToParenthesis(*tree1, true);
+						string temptre2 = TreeTemplateTools::treeToParenthesis(*tree2, true);
+						OUTPUT << "\n" << endl;
+						OUTPUT << "CAPS generated tree 1: " << temptre1 << endl;
+						OUTPUT << "CAPS generated tree 2: " << temptre2 << endl;
 					}/*else if(tree_in ==1 && variable==1){
 
 					   std::auto_ptr<DistanceMatrix> DS;
@@ -666,6 +679,7 @@ int main(int argc, char *argv[]){
 					int value = floor(((totaltemp.size())*(1-(threshval))))+1;
 
 					threshold = totaltemp[value];
+					totaltempnew = totaltemp;
 
 
 					/*=======================================================*/
@@ -870,6 +884,30 @@ int Chi_squared (int num_pairs, int num_
 
 }		/* -----  end of function Chi_squared  ----- */
 
+/* 
+ * ===  FUNCTION  ======================================================================
+ *         Name:  find_alpha
+ *  Description:  Find the index of an element in a vector totaltemp
+ *    Help from:  https://www.geeksforgeeks.org/how-to-find-index-of-a-given-element-in-a-vector-in-cpp/
+ *                https://stackoverflow.com/questions/8647635/elegant-way-to-find-closest-value-in-a-vector-from-above
+ *       Author:  Petar Petrov, University of Turku (Finland); pebope@utu.fi
+ * =====================================================================================
+ */
+double getIndex(std::vector<double> const& v, double K) 
+{ 
+    auto const it = std::lower_bound(v.begin(), v.end(), fabs(K));
+    //auto it = std::upper_bound(v.begin(), v.end(), fabs(K));
+  
+    if (it != v.end()) { 
+        int index = distance(v.begin(), it);
+	alphathresh = (((int)1+(double)v.size()-(int)index)/(double)v.size());
+	return alphathresh;
+        //cerr << index << "\t" << alphathresh << endl; 
+    } 
+    else { 
+        cerr << "ELEMENT NOT FOUND!" << endl; 
+    } 
+} 
 
 
 
@@ -890,9 +928,9 @@ int print_inter(vector<double>& Correl1,
 	output << endl << endl;
 
 	output << "Coevolving Pairs of amino acid sites\n";
-	output << "=============================================================================\n";
-	output << "Col1(real)\tCol2(real)\tDmean1\t\tDmean2\t\tCorrelation\tBootstrap value\n\n";
-	output << "=============================================================================\n";
+	output << "================================================================================================================================\n";
+	output << "Col1(real)\tCol2(real)\tDmean1\t\tDmean2\t\tCorrelation\tBootstrap value\tP-value1\tP-value2\tMean P-value\tCorrelation1\tCorrelation2\n\n";
+	output << "================================================================================================================================\n";
 
 	//double mean = average_vec<double>(Correl);
 	//double SD = SD_vf(Correl, mean);
@@ -951,9 +989,11 @@ int print_inter(vector<double>& Correl1,
 
 						//	}
 
+						double Alpha1 = getIndex(totaltempnew, Correl1[cor]);
+						double Alpha2 = getIndex(totaltempnew, Correl2[cor]);
 						//if(bootval>=bootcut && re1<=8 && re2<=8 ){
 						if(bootval>=bootcut){
-							output << i+1 << "(" << i-gaps1+1 << ")\t\t" << j+1 << "(" << (j+1)-gaps2 << ")\t\t" << averDi << "\t\t" << averDj << "\t\t" << (Correl1[cor]+Correl2[cor])/2 << "\t" << bootval << endl; 
+							output << i+1 << "(" << i-gaps1+1 << ")\t\t" << j+1 << "(" << (j+1)-gaps2 << ")\t\t" << averDi << "\t\t" << averDj << "\t" << (Correl1[cor]+Correl2[cor])/2 << "\t" << bootval << "\t" << Alpha1 << "\t" << Alpha2 << "\t" << (Alpha1+Alpha2)/2 << "\t" << Correl1[cor] << "\t" << Correl2[cor] << endl; 
 							signif.push_back(((Correl1[cor]+Correl2[cor])/2));
 							++pairs;
 							vector<int> tem;