add conditioning

src/omicron-x.cc

@@ -16,7 +16,7 @@ using namespace std;
 int main (int argc, char* argv[]){
 
   const string mapn = "map4";
-  const unsigned int overlap = 2;// s
+  const unsigned int overlap = 4;// s
   const double zero_lag_dt = 1.0;// s (+/-)
 
   if(argc<3) return -1;
@@ -25,17 +25,20 @@ int main (int argc, char* argv[]){
   string dir2 = (string)argv[2];
   bool savespectro=false;
   if(argc==4) savespectro=true;
-  TH2D *rho1, *rho2;
+  TH2D *rho1_raw, *rho2_raw, *rho1, *rho2;
   TH1D *xi, *xi_bg;
   TH1D *xi_bg_tot=NULL, *xi_bg_cum_tot=NULL, *xi_0_tot=NULL, *xi_0_cum_tot=NULL;
-  double window;
+  unsigned int nt, m_s;
+  double window, window_raw;
   double realpart, norm;
   double xi0_max=0;
   double dt;
   string gps1;
   fft *ft_1 = NULL;
   fft *ft_2 = NULL;
-
+  TF1 *fitfn = new TF1("fitfn", "gaus(0)");
+  double nsigma;
+  
   // get list of root files
   vector<string> file1;
   vector<string> file2;
@@ -72,36 +75,58 @@ int main (int argc, char* argv[]){
 
     // get rho1
     TFile *tf1 = new TFile(file1[f1].c_str(), "READ");
-    tf1->GetObject(mapn.c_str(), rho1);
-    rho1->SetDirectory(0);
+    tf1->GetObject(mapn.c_str(), rho1_raw);
+    rho1_raw->SetDirectory(0);
     tf1->Close();
 
     // get rho2
     TFile *tf2 = new TFile(fn2.c_str(), "READ");
-    tf2->GetObject(mapn.c_str(), rho2);
-    rho2->SetDirectory(0);
+    tf2->GetObject(mapn.c_str(), rho2_raw);
+    rho2_raw->SetDirectory(0);
     tf2->Close();
 
-    if(rho1 == 0){
-      cerr<<"rho1 == 0"<<endl;
+    if(rho1_raw == 0){
+      cerr<<"rho1_raw == 0"<<endl;
       return 1;
     }
-    if(rho2 == 0){
-      cerr<<"rho2 == 0"<<endl;
+    if(rho2_raw == 0){
+      cerr<<"rho2_raw == 0"<<endl;
       return 1;
     }
     
-    if(rho1->GetNbinsX()!=rho2->GetNbinsX()){
-      cerr<<"rho1->GetNbinsX()!=rho2->GetNbinsX()"<<endl;
+    if(rho1_raw->GetNbinsX()!=rho2_raw->GetNbinsX()){
+      cerr<<"rho1_raw->GetNbinsX()!=rho2_raw->GetNbinsX()"<<endl;
       return 1;
     }
     
-    if(rho1->GetNbinsY()!=rho2->GetNbinsY()){
-      cerr<<"rho1->GetNbinsY()!=rho2->GetNbinsY()"<<endl;
+    if(rho1_raw->GetNbinsY()!=rho2_raw->GetNbinsY()){
+      cerr<<"rho1_raw->GetNbinsY()!=rho2_raw->GetNbinsY()"<<endl;
       return 1;
     }
-    window = rho1->GetXaxis()->GetBinUpEdge(rho1->GetNbinsX()) - rho1->GetXaxis()->GetBinLowEdge(1);
-    cout<<"\trho map size: "<<rho1->GetNbinsX()<<" x "<<rho1->GetNbinsY()<<endl;
+
+    // conditioning (remove overlap)
+    window_raw = rho1_raw->GetXaxis()->GetBinUpEdge(rho1_raw->GetNbinsX()) - rho1_raw->GetXaxis()->GetBinLowEdge(1);
+    window = window_raw - (double)overlap;
+    nt = (unsigned int)(window/window_raw * (double)rho1_raw->GetNbinsX());
+    m_s = (unsigned int)((double)overlap/2.0/window_raw * (double)rho1_raw->GetNbinsX());
+    rho1 = new TH2D("rho1", "rho1", nt,
+                    rho1_raw->GetXaxis()->GetBinLowEdge(1)+(double)overlap/2.0,
+                    rho1_raw->GetXaxis()->GetBinLowEdge(1)+(double)overlap/2.0 + (double)nt*rho1_raw->GetXaxis()->GetBinWidth(1),
+                    rho1_raw->GetNbinsY(), rho1_raw->GetYaxis()->GetXbins()->GetArray());
+    rho2 = new TH2D("rho2", "rho2", nt,
+                    rho2_raw->GetXaxis()->GetBinLowEdge(1)+(double)overlap/2.0,
+                    rho2_raw->GetXaxis()->GetBinLowEdge(1)+(double)overlap/2.0 + (double)nt*rho2_raw->GetXaxis()->GetBinWidth(1),
+                    rho2_raw->GetNbinsY(), rho2_raw->GetYaxis()->GetXbins()->GetArray());
+    for(unsigned int l=0; l<rho1->GetNbinsY(); l++){
+      for(unsigned int m=0; m<rho1->GetNbinsX(); m++){
+        rho1->SetBinContent(m+1,l+1, rho1_raw->GetBinContent(m_s+m+1,l+1));
+        rho2->SetBinContent(m+1,l+1, rho2_raw->GetBinContent(m_s+m+1,l+1));
+      }
+    }
+    rho1->SetStats(0);
+    rho2->SetStats(0);
+    
+    cout<<"\trho map size: "<<rho1->GetNbinsX()<<" x "<<rho1->GetNbinsY()<<" "<<rho1_raw->GetYaxis()->GetXbins()->GetSize()<<endl;
     cout<<"\trho map time range: "<<window<<" s"<<endl;
 
     // Xi function
@@ -172,7 +197,13 @@ int main (int argc, char* argv[]){
     xi_0_tot->Fill(xi->GetBinContent(rho1->GetNbinsX()/2));
     if(xi->GetBinContent(rho1->GetNbinsX()/2)>xi0_max) xi0_max=xi->GetBinContent(rho1->GetNbinsX()/2);
     cout<<"Xi(0) = "<<xi->GetBinContent(rho1->GetNbinsX()/2)<<" (max = "<<xi0_max<<")"<<endl;
-  
+
+    // Gaussian fit
+    xi_bg->Fit(fitfn);
+    nsigma = (xi->GetBinContent(rho1->GetNbinsX()/2) - fitfn->GetParameter(1))/fitfn->GetParameter(2);
+    if(nsigma>2.0) cout<<">>>>>>>>>>>>>>>>>>> "<<nsigma<<" sigmas"<<endl;
+    else cout<<">> "<<nsigma<<" sigmas"<<endl;
+    
     outf->cd();
     xi->Write();
     xi_bg->Write();
@@ -183,6 +214,8 @@ int main (int argc, char* argv[]){
     
     delete xi;
     delete xi_bg;
+    delete rho1_raw;
+    delete rho2_raw;
     delete rho1;
     delete rho2;
   }
@@ -190,10 +223,16 @@ int main (int argc, char* argv[]){
   // cumulative distribution
   if(xi_bg_tot!=NULL){
     xi_bg_cum_tot = (TH1D*)xi_bg_tot->GetCumulative(kFALSE);
+    for(unsigned int m=0; m<xi_bg_cum_tot->GetNbinsX(); m++){
+      xi_bg_cum_tot->SetBinError(m+1, sqrt(xi_bg_cum_tot->GetBinContent(m+1)));
+    }
     xi_bg_cum_tot->Scale(1.0/(double)xi_bg_tot->GetEntries());
   }
   if(xi_0_tot!=NULL){
     xi_0_cum_tot = (TH1D*)xi_0_tot->GetCumulative(kFALSE);
+    for(unsigned int m=0; m<xi_0_cum_tot->GetNbinsX(); m++){
+      xi_0_cum_tot->SetBinError(m+1, sqrt(xi_0_cum_tot->GetBinContent(m+1)));
+    }
     xi_0_cum_tot->Scale(1.0/(double)xi_0_tot->GetEntries());
   }