mollerClass.C

Go to the documentation of this file.

#define mollerClass_cxx
#include "mollerClass.h"
#include <TH3.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <iostream>
#include <TGaxis.h>
#include <TProfile2D.h>
#include <TF1.h>
#include <TPaveLabel.h>
#include "profile.h"

using namespace std;

Int_t eventtype=-1;  //This comes from 0=moller, 1=proton, 2=inelastic
Int_t num_thrown=0;
char inifilename[256];

void SetupHist(TH1F* histo, Int_t colorint) 
{
  //histo->Sumw2();
  if (colorint==0) colorint=kRed;
  else if (colorint==1) colorint=kGreen;
  else if (colorint==2) colorint=kBlue;
  histo->SetLineColor(colorint);
  //histo->SetMarkerColor(colorint);
  //histo->SetMarkerStyle(20);
  histo->StatOverflows(0);
}

void PrintHist(TH1F* histo[], char* histoname, Int_t numhistos=4, Bool_t dology=0, Bool_t debug=0)
{
  char my_name[256];
  char my_name2[256];
  snprintf(my_name2,200,"r_%i_0",eventtype);
  sprintf(my_name,"%s",histo[0]->GetName());
  int comp = strcmp( my_name, my_name2);

  if (debug) {
    Printf("Debug: %s\n",histoname);
    for  (int i=0; i<5; i++) {
      printf("  %f,",histo[i]->GetEntries());
    }
    printf("\n");
  }
  TCanvas printcanvas("printcanvas","canvas",640,480);
  Double_t minnum=0, maxnum=0;
  for (int i=0; i<numhistos; i++) {
    minnum = TMath::Min(minnum,histo[i]->GetMinimum(0));
    maxnum = TMath::Max(maxnum,histo[i]->GetMaximum());
    //printf("new min = %f\n",minnum);
  }
  printf("min = %f, max = %f\n",minnum, maxnum); 

  Double_t xmin;
  xmin = 0.75;
  if(!comp) xmin = 0.65;
  Double_t xmax=0.999;  
  Double_t ymax=0.999, deltay=0.045, ymin=ymax-deltay;
  char label[256][numhistos];
  TPaveLabel *pt[numhistos];
  TPaveLabel *pt_main;
  pt_main = new TPaveLabel(xmin,ymin,xmax,ymax,"Entries  Mean  RMS","NDC");
    if(!comp) pt_main = new TPaveLabel(xmin,ymin,xmax,ymax,"Entries  Mean  RMS  Int","NDC");
  pt_main->SetBorderSize(1);
  pt_main->SetFillColor(0);
  ymax=ymin;
  ymin=ymax-deltay;
  Bool_t first=1;
  for (Int_t j = 0; j < numhistos; j++) {
    Int_t i = numhistos - 1 - j;
    sprintf(label[i],"%.0f  %.4g  %.4g",histo[i]->GetEntries(),histo[i]->GetMean(),histo[i]->GetRMS());
    if(!comp) sprintf(label[i],"%.0f  %.4g  %.4g  %.4g",histo[i]->GetEntries(),histo[i]->GetMean(),histo[i]->GetRMS(),histo[i]->Integral());

    pt[i] = new TPaveLabel(xmin,ymin,xmax,ymax,label[i],"NDC");
    pt[i]->SetBorderSize(1);
    pt[i]->SetFillColor(0);
    if(dology==1&&minnum==0) minnum = 0.01;
    histo[i]->SetMinimum(minnum);
    histo[i]->SetMaximum(1.05*maxnum);
    histo[i]->GetXaxis()->SetTitleOffset(1.2);
    if (i==0) {
      pt[i]->SetTextColor(kRed);
      histo[i]->SetLineColor(kRed);
    }
    if (i==1) {
      pt[i]->SetTextColor(kGreen);
      histo[i]->SetLineColor(kGreen);
    }
    if (i==2) {
      pt[i]->SetTextColor(kBlue);
      histo[i]->SetLineColor(kBlue);
    }
    if (i==3) {
      pt[i]->SetTextColor(kBlack);
      histo[i]->SetLineColor(kBlack);
    }
    if (i==4) {
      pt[i]->SetTextColor(7);
      histo[i]->SetLineColor(7);
    }
    if (first) {
      histo[i]->Draw();
      pt_main->Draw();
      first=0;
    } else {
      histo[i]->Draw("same");
    }
    pt[i]->Draw();
    ymax=ymin;
    ymin=ymax-deltay;
  }

  //  pt4->SetTextColor(7);
  //  pt1->SetTextColor(kRed);
  //  pt2->SetTextColor(kBlue);
  //  pt3->SetTextColor(kGreen);


  //  char label[256][numhistos];//, label_1[256], label_2[256], label_3[256], label_4[256];
  //  //sprintf(formatstring,"%%.0f  %%e  %%e",);
  //  //sprintf(label,  "%.0f  %.3g  %.3g",sumhisto->GetEntries(),sumhisto->GetMean(),sumhisto->GetRMS());
  //  sprintf(label_4,"%.0f  %.3g  %.3g",histo[4]->GetEntries(),histo[4]->GetMean(),histo[4]->GetRMS());
  //  sprintf(label  ,"%.0f  %.3g  %.3g",histo[3]->GetEntries(),histo[3]->GetMean(),histo[3]->GetRMS());
  //  sprintf(label_1,"%.0f  %.3g  %.3g",histo[0]->GetEntries(),histo[0]->GetMean(),histo[0]->GetRMS());  
  //  sprintf(label_2,"%.0f  %.3g  %.3g",histo[2]->GetEntries(),histo[2]->GetMean(),histo[2]->GetRMS());
  //  sprintf(label_3,"%.0f  %.3g  %.3g",histo[1]->GetEntries(),histo[1]->GetMean(),histo[1]->GetRMS());
  //  Double_t xmin=0.75, xmax=0.999;  
  //  Double_t ymax=0.999, deltay=0.045, ymin=ymax-deltay;
  //  TPaveLabel *pt0 = new TPaveLabel(xmin,ymin,xmax,ymax,"Entries  Mean  RMS","NDC");
  //  pt0->SetBorderSize(1);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  TPaveLabel *pt4 = new TPaveLabel(xmin,ymin,xmax,ymax,label_4,"NDC");
  //  pt4->SetBorderSize(1);
  //  pt4->SetTextColor(7);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  TPaveLabel *pt = new TPaveLabel(xmin,ymin,xmax,ymax,label,"NDC");
  //  pt->SetBorderSize(1);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  TPaveLabel *pt1 = new TPaveLabel(xmin,ymin,xmax,ymax,label_1,"NDC");
  //  pt1->SetBorderSize(1);
  //  pt1->SetTextColor(kRed);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  TPaveLabel *pt2 = new TPaveLabel(xmin,ymin,xmax,ymax,label_2,"NDC");
  //  pt2->SetBorderSize(1);
  //  pt2->SetTextColor(kBlue);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  TPaveLabel *pt3 = new TPaveLabel(xmin,ymin,xmax,ymax,label_3,"NDC");
  //  pt3->SetBorderSize(1);
  //  pt3->SetTextColor(kGreen);
  //  ymax=ymin;
  //  ymin=ymax-deltay;
  //  pt0->SetFillColor(0); pt->SetFillColor(0);  pt1->SetFillColor(0);  
  //  pt2->SetFillColor(0);  pt3->SetFillColor(0);
  //  pt4->SetFillColor(0);
  //  pt0->Draw();
  //  pt->Draw();
  //  pt1->Draw();
  //  pt2->Draw();
  //  pt3->Draw();
  //  pt4->Draw();


  printcanvas.SetLogy(dology);
  printcanvas.Print(histoname);
}

void PrintHist2D(TH2F* histo[], char* histoname)
{
  TCanvas printcanvas("printcanvas","canvas",640,480);
  histo[0]->SetMarkerColor(kRed);
  histo[1]->SetMarkerColor(kGreen);
  histo[2]->SetMarkerColor(kBlue);
  histo[0]->SetLineColor(kRed);
  histo[1]->SetLineColor(kGreen);
  histo[2]->SetLineColor(kBlue);
  histo[3]->Draw("contourz");
  char histoname2[200];
  snprintf(histoname2,200,"%s_black.gif",histoname);
  printcanvas.Print(histoname2);
  histo[0]->Draw("contourz");
  snprintf(histoname2,200,"%s_red.gif",histoname);
  printcanvas.Print(histoname2);
  histo[1]->Draw("contourz");
  snprintf(histoname2,200,"%s_green.gif",histoname);
  printcanvas.Print(histoname2);
  histo[2]->Draw("contourz");
  snprintf(histoname2,200,"%s_blue.gif",histoname);
  printcanvas.Print(histoname2);
  histo[4]->Draw("contourz");
  //histo[0]->Draw("contour,same");
  //histo[1]->Draw("contour,same");
  //histo[2]->Draw("contour,same");
  printcanvas.Print(histoname);
}

void mollerClass::Loop()
{

  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntriesFast();

  TCanvas canvas("canvas","canvas",640,480);

  Double_t pi = TMath::Pi();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(kFALSE);
  gStyle->SetPalette(1,0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetCanvasBorderMode(0);

  Int_t num_detectors = 28;
  char plotname[200],outdir[200],namestem[200],plottype[10];


  snprintf(namestem,200,"%sdets%i_",outdir,num_detectors);
  GetPrivateProfileString ("General", "outdir", "output", outdir, 200, inifilename);
  GetPrivateProfileString ("General", "plottype", ".gif", plottype, 20, inifilename);


  TH1F* std_ave_asym=new TH1F("ave_asym","Average transverse asymmetry;detector number;#bar{A_{T}}   (ppm)",
                num_detectors/8,0,num_detectors); 
  TH1F* std_total_asym=new TH1F("total_asym","Total XSec and transverse asymmetry weighted events;detector number",
                  num_detectors/8,0,num_detectors); 
  TH1F* std_total_sample=new TH1F("total_sample","Total XSec weighted events;detector number",
                  num_detectors/8,0,num_detectors); 
  std_ave_asym->Sumw2();
  std_total_asym->Sumw2();
  std_total_sample->Sumw2();
  std_ave_asym->SetMarkerStyle(20);
  std_total_asym->SetMarkerStyle(20);
  std_total_sample->SetMarkerStyle(20);
  TH1F* std_red_ave_asym=new TH1F("std_red_ave_asym","Average transverse asymmetry;detector number;#bar{A_{T}}   (ppm)",
                  num_detectors,0,num_detectors); 
  TH1F* std_red_total_asym=new TH1F("std_red_total_asym","Total (XSec and transverse asymmetry weighted) events;detector number",
                    num_detectors,0,num_detectors); 
  TH1F* std_red_total_sample=new TH1F("std_red_total_sample","Total (XSec weighted) events;detector number",
                    num_detectors,0,num_detectors); 
  std_red_ave_asym->Sumw2();
  std_red_total_asym->Sumw2();
  std_red_total_sample->Sumw2();
  std_red_ave_asym->SetLineColor(kRed);
  std_red_total_asym->SetLineColor(kRed);
  std_red_total_sample->SetLineColor(kRed);
  std_red_ave_asym->SetMarkerColor(kRed);
  std_red_total_asym->SetMarkerColor(kRed);
  std_red_total_sample->SetMarkerColor(kRed);
  std_red_ave_asym->SetMarkerStyle(20);
  std_red_total_asym->SetMarkerStyle(20);
  std_red_total_sample->SetMarkerStyle(20);
  TH1F* std_green_ave_asym=new TH1F("std_green_ave_asym","Average transverse asymmetry;detector number;#bar{A_{T}}   (ppm)",
                    num_detectors,0,num_detectors); 
  TH1F* std_green_total_asym=new TH1F("std_green_total_asym","Total (XSec and transverse asymmetry weighted) events;detector number",
                    num_detectors,0,num_detectors); 
  TH1F* std_green_total_sample=new TH1F("std_green_total_sample","Total (XSec weighted) events;detector number",
                      num_detectors,0,num_detectors); 
  std_green_ave_asym->Sumw2();
  std_green_total_asym->Sumw2();
  std_green_total_sample->Sumw2();
  std_green_ave_asym->SetLineColor(kGreen);
  std_green_total_asym->SetLineColor(kGreen);
  std_green_total_sample->SetLineColor(kGreen);
  std_green_ave_asym->SetMarkerColor(kGreen);
  std_green_total_asym->SetMarkerColor(kGreen);
  std_green_total_sample->SetMarkerColor(kGreen);
  std_green_ave_asym->SetMarkerStyle(20);
  std_green_total_asym->SetMarkerStyle(20);
  std_green_total_sample->SetMarkerStyle(20);
  TH1F* std_blue_ave_asym=new TH1F("std_blue_ave_asym","Average transverse asymmetry;detector number;#bar{A_{T}}   (ppm)",
                   num_detectors,0,num_detectors); 
  TH1F* std_blue_total_asym=new TH1F("std_blue_total_asym","Total (XSec and transverse asymmetry weighted) events;detector number",
                     num_detectors,0,num_detectors); 
  TH1F* std_blue_total_sample=new TH1F("std_blue_total_sample","Total (XSec weighted) events;detector number",
                     num_detectors,0,num_detectors); 
  std_blue_ave_asym->Sumw2();
  std_blue_total_asym->Sumw2();
  std_blue_total_sample->Sumw2();
  std_blue_ave_asym->SetLineColor(kBlue);
  std_blue_total_asym->SetLineColor(kBlue);
  std_blue_total_sample->SetLineColor(kBlue);
  std_blue_ave_asym->SetMarkerColor(kBlue);
  std_blue_total_asym->SetMarkerColor(kBlue);
  std_blue_total_sample->SetMarkerColor(kBlue);
  std_blue_ave_asym->SetMarkerStyle(20);
  std_blue_total_asym->SetMarkerStyle(20);
  std_blue_total_sample->SetMarkerStyle(20);




  TH1F* eweight_ave_asym=new TH1F("eweight_ave_asym","Average transverse asymmetry, energy weighted detectors;detector number;#bar{A_{T}}   (ppm)",
                  num_detectors/8,0,num_detectors); 
  TH1F* eweight_total_asym=new TH1F("eweight_total_asym","Total (energy, XSec and transverse asymmetry weighted) events;detector number",
                    num_detectors/8,0,num_detectors); 
  TH1F* eweight_total_sample=new TH1F("eweight_total_sample","Total (energy and XSec weighted) events;detector number",
                    num_detectors/8,0,num_detectors); 
  eweight_ave_asym->Sumw2();
  eweight_total_asym->Sumw2();
  eweight_total_sample->Sumw2();
  eweight_ave_asym->SetMarkerStyle(20);
  eweight_total_asym->SetMarkerStyle(20);
  eweight_total_sample->SetMarkerStyle(20);
  TH1F* eweight_red_ave_asym=new TH1F("eweight_red_ave_asym","Average transverse asymmetry, energy weighted detectors;detector number;#bar{A_{T}}   (ppm)",
                    num_detectors,0,num_detectors); 
  TH1F* eweight_red_total_asym=new TH1F("eweight_red_total_asym","Total (energy, XSec and transverse asymmetry weighted) events;detector number",
                      num_detectors,0,num_detectors); 
  TH1F* eweight_red_total_sample=new TH1F("eweight_red_total_sample","Total (energy and XSec weighted) events;detector number",
                      num_detectors,0,num_detectors); 
  eweight_red_ave_asym->Sumw2();
  eweight_red_total_asym->Sumw2();
  eweight_red_total_sample->Sumw2();
  eweight_red_ave_asym->SetLineColor(kRed);
  eweight_red_total_asym->SetLineColor(kRed);
  eweight_red_total_sample->SetLineColor(kRed);
  eweight_red_ave_asym->SetMarkerColor(kRed);
  eweight_red_total_asym->SetMarkerColor(kRed);
  eweight_red_total_sample->SetMarkerColor(kRed);
  eweight_red_ave_asym->SetMarkerStyle(20);
  eweight_red_total_asym->SetMarkerStyle(20);
  eweight_red_total_sample->SetMarkerStyle(20);
  TH1F* eweight_green_ave_asym=new TH1F("eweight_green_ave_asym","Average transverse asymmetry, energy weighted detectors;detector number;#bar{A_{T}}   (ppm)",
                      num_detectors,0,num_detectors); 
  TH1F* eweight_green_total_asym=new TH1F("eweight_green_total_asym","Total (energy, XSec and transverse asymmetry weighted) events;detector number",
                      num_detectors,0,num_detectors); 
  TH1F* eweight_green_total_sample=new TH1F("eweight_green_total_sample","Total (energy and XSec weighted) events;detector number",
                        num_detectors,0,num_detectors); 
  eweight_green_ave_asym->Sumw2();
  eweight_green_total_asym->Sumw2();
  eweight_green_total_sample->Sumw2();
  eweight_green_ave_asym->SetLineColor(kGreen);
  eweight_green_total_asym->SetLineColor(kGreen);
  eweight_green_total_sample->SetLineColor(kGreen);
  eweight_green_ave_asym->SetMarkerColor(kGreen);
  eweight_green_total_asym->SetMarkerColor(kGreen);
  eweight_green_total_sample->SetMarkerColor(kGreen);
  eweight_green_ave_asym->SetMarkerStyle(20);
  eweight_green_total_asym->SetMarkerStyle(20);
  eweight_green_total_sample->SetMarkerStyle(20);
  TH1F* eweight_blue_ave_asym=new TH1F("eweight_blue_ave_asym","Average transverse asymmetry, energy weighted detectors;detector number;#bar{A_{T}}   (ppm)",
                     num_detectors,0,num_detectors); 
  TH1F* eweight_blue_total_asym=new TH1F("eweight_blue_total_asym","Total (energy, XSec and transverse asymmetry weighted) events;detector number",
                       num_detectors,0,num_detectors); 
  TH1F* eweight_blue_total_sample=new TH1F("eweight_blue_total_sample","Total (energy and XSec weighted) events;detector number",
                       num_detectors,0,num_detectors); 
  eweight_blue_ave_asym->Sumw2();
  eweight_blue_total_asym->Sumw2();
  eweight_blue_total_sample->Sumw2();
  eweight_blue_ave_asym->SetLineColor(kBlue);
  eweight_blue_total_asym->SetLineColor(kBlue);
  eweight_blue_total_sample->SetLineColor(kBlue);
  eweight_blue_ave_asym->SetMarkerColor(kBlue);
  eweight_blue_total_asym->SetMarkerColor(kBlue);
  eweight_blue_total_sample->SetMarkerColor(kBlue);
  eweight_blue_ave_asym->SetMarkerStyle(20);
  eweight_blue_total_asym->SetMarkerStyle(20);
  eweight_blue_total_sample->SetMarkerStyle(20);


  TH1F* rate_total=new TH1F("rate_total","Total rate for 85  #muA;detector number;rate   [Hz]",num_detectors,0,num_detectors);
  TH1F* num_events=new TH1F("num_events","Number of Moller e^{-}  in 1 hour at 85  #muA;detector number",num_detectors,0,num_detectors);
  TH1F* error=new TH1F("error","Uncertainty in 1 hour at 85  #muA;detector number;#Delta A_{T}   (ppm)",num_detectors,0,num_detectors);

  rate_total->SetMarkerStyle(20);
  num_events->SetMarkerStyle(20);
  error->SetMarkerStyle(20);
  rate_total->Sumw2();
  num_events->Sumw2();
  error->Sumw2();
  //  std_total_asym->SetMarkerStyle(20);
  //  std_total_sample->SetMarkerStyle(20);



  TH2F* test=new TH2F("test","test",num_detectors*2,-pi,pi,num_detectors*2,0,num_detectors);
  TH2F* paperfig=new TH2F("paperfig","A_{T}/sin#phi for E_{vert}=11 GeV;cos#theta_{c.m.};A_{T}/sin#phi   (ppm)",
              400,-1,1,400,-15,15);
  TH3F* paperfig2=new TH3F("paperfig2","A_{T};cos#theta_{c.m.};#phi_{c.m.};A_{T}   (ppm)",
               200,-1,1,200,-pi,pi,200,-15,15);
  TH2F* distro_red=new TH2F("distro_red","Distribution of red events;cos#theta_{c.m.};#phi_{c.m.}",
                400,-1,1,400,-pi,pi);
  distro_red->SetMarkerColor(kRed);
  TH2F* distro_green=new TH2F("distro_green","Distribution of green events;cos#theta_{c.m.};#phi_{c.m.}",
                400,-1,1,400,-pi,pi);
  distro_green->SetMarkerColor(kGreen);
  TH2F* distro_blue=new TH2F("distro_blue","Distribution of blue events;cos#theta_{c.m.};#phi_{c.m.}",
                 400,-1,1,400,-pi,pi);
  distro_blue->SetMarkerColor(kBlue);


  TH2F* distro_rphi=new TH2F("distro_rphi",";r   (m);#phi_{c.m.}", 400,0,1.3,400,-pi,pi);
  TH2F* distro_r_edet=new TH2F("distro_r_edet",";r   (m);E_{det}   (MeV)", 400,0,1.3,400,0,11000);
  TH2F* distro_r_evert=new TH2F("distro_r_evert",";r   (m);E_{vert}   (MeV)", 400,0,1.3,400,0,11000);
  TH1F* distro_r;
  if (eventtype==0) {
    distro_r=new TH1F("moller_distro_r",";r   (m);rate   (GHz)", 260,0,1.3);
  } else {
    if (eventtype==1) {
      distro_r=new TH1F("mott_distro_r",";r   (m);rate   (GHz)", 260,0,1.3);
    } else {
      if (eventtype==2) {
        distro_r=new TH1F("epinel_distro_r",";r   (m);rate   (GHz)", 260,0,1.3);
      }
    }
  }



  // **** Each quantity gets 5 histograms, within the r cut the three detectors [0],[1],[2] their sum [3] and 
  // **** outside the the r cut [4].
  Int_t numhistos = 5;
  TH1F *theta_cm[numhistos], *Eprime_vert_hist[numhistos], *Eprime_det_hist[numhistos], *theta_lab_hist[numhistos];
  TH1F *phi_wrap_hist[numhistos], *Q2_hist[numhistos], *W2_hist[numhistos], *W_hist[numhistos], *asym_hist[numhistos], *depol_hist[numhistos], *asymd_hist[numhistos], *asymEep_hist[numhistos],*asymIep_hist[numhistos],*z0_hist[numhistos];
  TH1F *Q2_eweight_hist[numhistos], *W2_Q2weight_hist[numhistos];
  TH1F *r_hist[numhistos], *r_asymtot[numhistos], *r_asymave[numhistos], *r_eweight[numhistos];
  TH2F *r_phi_hist[numhistos],*r_asym_hist[numhistos];
  char name[200];
  for (Int_t i=0; i<numhistos; i++) {
    // **** Setup the 1D histograms
    snprintf(name,200,"theta_cm_%i_%i",eventtype,i);
    theta_cm[i]=new TH1F(name,";#theta_{c.m.};rate   (GHz)/deg", 180,0,180);
    SetupHist(theta_cm[i],i);
    snprintf(name,200,"Eprime_vert_%i_%i",eventtype,i);
    Eprime_vert_hist[i]=new TH1F(name,";E'_{vert}   (GeV);rate   (GHz)/100MeV",110,0,11);
    SetupHist(Eprime_vert_hist[i],i);
    snprintf(name,200,"Eprime_det_%i_%i",eventtype,i);
    Eprime_det_hist[i]=new TH1F(name,";E'_{det}   (GeV);rate   (GHz)/100MeV",110,0,11);
    SetupHist(Eprime_det_hist[i],i);
    snprintf(name,200,"theta_lab_%i_%i",eventtype,i);
    theta_lab_hist[i]=new TH1F(name,";#theta_{lab}   (millirad);rate   (GHz)/0.3mrad", 100, 0,30);
    SetupHist(theta_lab_hist[i],i);
    snprintf(name,200,"phi_wrap_%i_%i",eventtype,i);
    phi_wrap_hist[i]=new TH1F(name,";#phi_{wrap};rate   (GHz)/0.402deg", 128,-25.7142857142857153,25.7142857142857153);
    SetupHist(phi_wrap_hist[i],i);
    snprintf(name,200,"Q2_%i_%i",eventtype,i);
    Q2_hist[i]=new TH1F(name,";Q^{2}   (GeV/c)^{2};rate   (GHz)/80(MeV/c)^{2}",200,0,0.016);
    SetupHist(Q2_hist[i],i);
    snprintf(name,200,"Q2_eweight_%i_%i",eventtype,i);
    Q2_eweight_hist[i]=new TH1F(name,";Q^{2}   (GeV/c)^{2};Energy weighted rate   (MeV*GHz)/80(MeV/c)^{2}",200,0,0.016);
    SetupHist(Q2_eweight_hist[i],i);
    snprintf(name,200,"W2_Q2weight_%i_%i",eventtype,i);
    W2_Q2weight_hist[i]=new TH1F(name,";W^{2}   (GeV^{2});Q^{2} weighted rate   ((GeV/c)^{2} *GHz)/0.1GeV^{2}",200,0,20);//0.880352,0.880357);
    SetupHist(W2_Q2weight_hist[i],i);
    snprintf(name,200,"W2_%i_%i",eventtype,i);
    W2_hist[i]=new TH1F(name,";W^{2}   (GeV^{2});rate   (GHz)/0.1GeV^{2}",200,0,20);//,0.880352,0.880357);
    SetupHist(W2_hist[i],i);
    snprintf(name,200,"W_%i_%i",eventtype,i);
    W_hist[i]=new TH1F(name,";W   (GeV);rate   (GHz)/30MeV",200,0,6);
    SetupHist(W_hist[i],i);
    snprintf(name,200,"asym_hist_%i_%i",eventtype,i);
    asym_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/0.2ppb",300,0,60);
    SetupHist(asym_hist[i],i);
    snprintf(name,200,"depol_hist_%i_%i",eventtype,i);
    depol_hist[i]=new TH1F(name,";Depolarization Factor  ;rate   (GHz)/1%",125,0,1.25);
    SetupHist(depol_hist[i],i);
    snprintf(name,200,"asymd_hist_%i_%i",eventtype,i);
    if(eventtype==0) asymd_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/0.2ppb",300,0,60);
    else if(eventtype==1) asymd_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/3ppb",200,0,600);
    else if(eventtype==2) asymd_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/5ppb",400,0,2000);
    SetupHist(asymd_hist[i],i);

    snprintf(name,200,"asymEep_hist_%i_%i",eventtype,i);
    asymEep_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/3ppb",200,0,600);
    SetupHist(asymEep_hist[i],i);
    snprintf(name,200,"asymIep_hist_%i_%i",eventtype,i);
    asymIep_hist[i]=new TH1F(name,";A_{PV}   (ppb);rate   (GHz)/5ppb",400,0,2000);
    SetupHist(asymIep_hist[i],i);

    snprintf(name,200,"z0_hist_%i_%i",eventtype,i);
    z0_hist[i]=new TH1F(name,";Target z   (mm);rate   (GHz)/cm",160,-800,800);
    SetupHist(z0_hist[i],i);
    // **** 1D histogrms for averages
    snprintf(name,200,"r_%i_%i",eventtype,i);
    r_hist[i]=new TH1F(name,";r   (m);rate   (GHz)/5mm",140,0.6,1.3);
    SetupHist(r_hist[i],i);
    snprintf(name,200,"r_eweight_%i_%i",eventtype,i);
    r_eweight[i]=new TH1F(name,";r   (m);Energy weighted rate   (MeV*GHz)/5mm",140,0.6,1.3);
    SetupHist(r_eweight[i],i);
    snprintf(name,200,"r_asymtot_%i_%i",eventtype,i);
    r_asymtot[i]=new TH1F(name,";r   (m);Asymmetry weighted rate (ppb*GHz)/5mm",140,0.6,1.3);
    SetupHist(r_asymtot[i],i);
    snprintf(name,200,"r_asymave_%i_%i",eventtype,i);
    r_asymave[i]=new TH1F(name,";r   (m);#bar{A_{PV}}   (ppb)/5mm",140,0.6,1.3);
    SetupHist(r_asymave[i],i);
    // **** Setup the 2D histograms
    snprintf(name,200,"r_phi_%i_%i",eventtype,i);
    r_phi_hist[i]=new TH2F(name,";#phi   (degrees);r   (m)",224,0,360,210,0.6,1.3);
    snprintf(name,200,"r_asym_%i_%i",eventtype,i);
    r_asym_hist[i]=new TH2F(name,";asym;r   (m)",42,0,42,70,0.6,1.3);
  }

  Long_t totalevents=0, numevcut_rcut=0, numevcut_thetacut=0, numevcut_volumecut=0,numevcut_kinecut=0,numevallpass=0;

  // *****  Here begins the loop

  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;

    totalevents++;

    Double_t r = sqrt(x*x+y*y);
//    Double_t phi_wrap  = TMath::ATan2(y,x)-pi/3.5*(TMath::Floor((TMath::ATan2(y,x)+pi/7)/(pi/3.5)));
    Double_t phi_wrap  = (TMath::ATan2(y,x)-pi)-pi/3.5*(TMath::Floor((TMath::ATan2(y,x)+pi/7-pi)/(pi/3.5)));
    //    Double_t phi_org_wrap = TMath::ATan2(y0,x0)-pi/3.5*(TMath::Floor((TMath::ATan2(y0,x0)+pi/7)/(pi/3.5)));
    Double_t theta_lab;
    if (eventtype!=2) {
      if(kineE2==kineE_org&&kineE1==kineE_org) {
        theta_lab = acos(pz2/sqrt(px2*px2+py2*py2+pz2*pz2)) + acos(pz1/sqrt(px1*px1+py1*py1+pz1*pz1));  
      }else if(kineE2==kineE_org&&kineE1!=kineE_org) {
        theta_lab = acos(pz2/sqrt(px2*px2+py2*py2+pz2*pz2));
      }else if(kineE2!=kineE_org&&kineE1==kineE_org) {
        theta_lab = acos(pz1/sqrt(px1*px1+py1*py1+pz1*pz1));  
      }else {
        theta_lab = 0;
      }
    } else {
      theta_lab = acos(pz1/sqrt(px1*px1+py1*py1+pz1*pz1));
    }
//    if (eventtype!=2) {
//      theta_lab = (kineE2==kineE_org)*acos(pz2/sqrt(px2*px2+py2*py2+pz2*pz2)) +
//        (kineE1==kineE_org)*acos(pz1/sqrt(px1*px1+py1*py1+pz1*pz1));  
//      // only if (kineE2==kineE_org) does px2 represent the detected particle
//    } else {
//      theta_lab = acos(pz1/sqrt(px1*px1+py1*py1+pz1*pz1));
//    }

    Double_t Q2_org = 4*kineE0/1000*(kineE_org/1000)*sin(theta_lab/2)*sin(theta_lab/2);
    Double_t mp=0.938272;
    Double_t nu=(kineE0-kineE_org)/1000;
    Double_t W_org = sqrt(-Q2_org + mp*mp + 2*mp*nu);
    Double_t y_org = 1-kineE_org/kineE0;
    Double_t asymMol_org = 10e8*0.0000116639/sqrt(2)/3.14159/0.0072973*(1-4*0.2388)*Q2_org*(1-y_org)/(1+pow(y_org,4)+pow((1-y_org),4)); 
    //10e8 to give ppb
    Double_t Din = 1.0 - (2./3.)*(11000. - kineE0)*(11000. - kineE0)/(kineE0*kineE0+11000.*11000. - (2./3.)*kineE0*11000.);

    Double_t asymEep_org = 1E9*0.0000116639/sqrt(2)/3.14159/0.0072973*(0.072/2.0)*Q2_org;
    Double_t asymIep_org = 1E9*(0.8E-4)*Q2_org;

//Set proper Asymetry for the generator eventtype
    Double_t asym_org;
    if(eventtype==0) asym_org = asymMol_org;
    else if(eventtype==1) asym_org = asymEep_org;
    else if(eventtype==2) asym_org = asymIep_org;

    Double_t asymd_org = Din*asym_org; //depolarization effect

    Double_t ecm = sqrt(2*0.000511*0.000511+2.*0.000511*kineE0/1000);
    Double_t gammap = kineE0/1000/ecm;
    Double_t betap = 1/gammap*sqrt(gammap*gammap-1);
    Double_t e1cm = gammap*0.000511;
    Double_t pcm2_org = (kineE_org/1000-gammap*e1cm)/gammap/betap;
    Double_t costheta_cms_org = pcm2_org/e1cm;

    Double_t theta_cms_org = acos(pcm2_org/e1cm);
    Double_t phi_org =  (kineE2==kineE_org)*TMath::ATan2(py2,px2)+(kineE1==kineE_org)*TMath::ATan2(py1,px1);
    Double_t phi_det = TMath::ATan2(y,x);
    Double_t phi_det_deg = (phi_det+TMath::Pi())*180/TMath::Pi();
    Double_t m_e = 0.000511;
    Double_t E_vert = kineE0/1000;

    Double_t s = 2*m_e*E_vert;
    Double_t t = -s/2*(1-costheta_cms_org);
    Double_t u = -s/2*(1+costheta_cms_org);
    Double_t alpha = 1/137.;
    Double_t dsigma_born = alpha*alpha/2/s*((t*t+t*u+u*u)/t/u)*((t*t+t*u+u*u)/t/u);
    Double_t bracket_term = 3*s*(t*(u-s)*TMath::Log(-t/s)-u*(t-s)*TMath::Log(-u/s))-2*(t-u)*t*u;
    Double_t proton_term = -alpha*alpha*alpha/8*m_e/sqrt(s)*sin(theta_cms_org)*sin(phi_org)/t/t/u/u;
    Double_t dsigma_phi = proton_term*bracket_term;
    Double_t asym_trans = dsigma_phi/dsigma_born*1000000;


    Double_t rate_GHz = 85*rate/(num_thrown)/1000000000;
    //    Double_t pi = TMath::Pi();

        //This may need to be changed depending on the field map:
    //Bool_t rcut = (r>880 && r<1000);// for proposal maps
    //Bool_t rcut = (r>900 && r<1100);// for default TOSCA 
    //Bool_t rcut = (r>700 && r<740);// 
    //Bool_t rcut = (r>740 && r<820);// 
    //Bool_t rcut = (r>820 && r<860);// 
    //Bool_t rcut = (r>860 && r<900);// 
    //Bool_t rcut = (r>900 && r<1060);// 
    //Bool_t rcut = (r>1060 && r<1160);// 
    Bool_t rcut = (r>0 && r<10000);// 

    //Bool_t thetacut = (theta_lab>0.0095);
    Bool_t thetacut = (theta_lab>0.004);
    Bool_t volumecut;
    if (eventtype==0) {
      volumecut = volume==0;
    } else {
      volumecut = ((volume==0)&&(track==1));
    }
//    Bool_t volumecut = volume==1;
    Bool_t kinecut = kineE_org>0;

    if (! rcut) numevcut_rcut++;
    if (! thetacut) numevcut_thetacut++;
    if (! volumecut) numevcut_volumecut++;
    if (! kinecut) numevcut_kinecut++;

    Bool_t maincut = volumecut && kinecut && thetacut;
    if (maincut && !thetacut) {
      printf("theta_lab=%f  kineE2=%f, kineE1=%f, kineE_org=%f   px1=%f, py1=%f, pz1=%f,  px2=%f, py2=%f, pz2=%f\n",
           theta_lab, kineE2, kineE1, kineE_org, px1, py1, pz1, px2, py2, pz2);
    }

    if (eventtype==2) {
      maincut = maincut && kineE2==0;
    }
    if (maincut) { // This is for stuff that may not make into the detectors

      distro_r->Fill(r/1000,rate_GHz);
      distro_rphi->Fill(r/1000,phi_org,rate_GHz);
      distro_r_evert->Fill(r/1000,kineE_org,rate_GHz);
      distro_r_edet->Fill(r/1000,kineE,rate_GHz);
      theta_cm[4]->Fill(acos(costheta_cms_org)*180/TMath::Pi(),rate_GHz);
      Eprime_vert_hist[4]->Fill(kineE_org/1000,rate_GHz);
      Eprime_det_hist[4]->Fill(kineE/1000,rate_GHz);
      theta_lab_hist[4]->Fill(theta_lab*1000,rate_GHz);
      phi_wrap_hist[4]->Fill(phi_wrap*180/TMath::Pi(),rate_GHz);
      Q2_hist[4]->Fill(Q2_org,rate_GHz);
      Q2_eweight_hist[4]->Fill(Q2_org,kineE*rate_GHz);
      W2_Q2weight_hist[4]->Fill(W_org*W_org,Q2_org*rate_GHz);
      W2_hist[4]->Fill(W_org*W_org,rate_GHz);
      W_hist[4]->Fill(W_org,rate_GHz);
      r_hist[4]->Fill(r/1000,rate_GHz);
      r_eweight[4]->Fill(r/1000,kineE*rate_GHz);
      r_asymtot[4]->Fill(r/1000,rate_GHz*asym_org);
      asym_hist[4]->Fill(asym_org,rate_GHz);
      depol_hist[4]->Fill(Din,rate_GHz);
      asymd_hist[4]->Fill(asymd_org,rate_GHz);
      asymEep_hist[4]->Fill(asymEep_org,rate_GHz);
      asymIep_hist[4]->Fill(asymIep_org,rate_GHz);
      z0_hist[4]->Fill(z0,rate_GHz);

      if (!(rcut)) {
        r_phi_hist[3]->Fill(phi_det_deg,r/1000,rate_GHz);
        r_asym_hist[3]->Fill(asym_org,r/1000,rate_GHz);
      }
      r_phi_hist[4]->Fill(phi_det_deg,r/1000,rate_GHz);
      r_asym_hist[4]->Fill(asym_org,r/1000,rate_GHz);
    }
    if (rcut and maincut) {
      numevallpass++;
      //if (1) printf("|");
      Double_t fillval;
      //fillval=num_detectors*((phi_det+pi)/(2.*pi));
      fillval=num_detectors*((phi_det+pi)/(2.*pi)-(1./56.));
      if (fillval<0) {
        //printf("%f  %f\n",fillval,fillval + num_detectors);
        fillval = fillval + num_detectors;
      }
      std_total_asym->Fill(fillval,asym_trans*totXS);
      std_total_sample->Fill(fillval,totXS);
      eweight_total_asym->Fill(fillval,kineE*asym_trans*totXS);
      eweight_total_sample->Fill(fillval,kineE*totXS);
      rate_total->Fill(fillval,85*rate/(num_thrown));
      
      num_events->Fill(fillval,60*60*85*rate/(num_thrown));


      if (TMath::Abs(phi_wrap)<pi/28.) {
        r_phi_hist[0]->Fill(phi_det_deg,r/1000,rate_GHz);
        r_asym_hist[0]->Fill(asym_org,r/1000,rate_GHz);
        //r_phi_hist[0]->Fill(r/1000,phi_det_deg);
      } else {
        if (TMath::Abs(phi_wrap)>pi/28.&&TMath::Abs(phi_wrap)<3*pi/28.) {
          r_phi_hist[1]->Fill(phi_det_deg,r/1000,rate_GHz);
          r_asym_hist[1]->Fill(asym_org,r/1000,rate_GHz);
          //r_phi_hist[1]->Fill(r/1000,phi_det_deg);
        } else {
          if (TMath::Abs(phi_wrap)>3*pi/28.) {
            //r_phi_hist[2]->Fill(r/1000,phi_det_deg);
            r_phi_hist[2]->Fill(phi_det_deg,r/1000,rate_GHz);
            r_asym_hist[2]->Fill(asym_org,r/1000,rate_GHz);
          }
        }
      }
      
      //if (Q2_org<0 || Q2_org>
      theta_cm[3]->Fill(acos(costheta_cms_org)*180/TMath::Pi(),rate_GHz);
      Eprime_vert_hist[3]->Fill(kineE_org/1000,rate_GHz);
      Eprime_det_hist[3]->Fill(kineE/1000,rate_GHz);
      theta_lab_hist[3]->Fill(theta_lab*1000,rate_GHz);
      phi_wrap_hist[3]->Fill(phi_wrap*180/TMath::Pi(),rate_GHz);
      Q2_hist[3]->Fill(Q2_org,rate_GHz);
      Q2_eweight_hist[3]->Fill(Q2_org,kineE*rate_GHz);
      W2_Q2weight_hist[3]->Fill(W_org*W_org,Q2_org*rate_GHz);
      W2_hist[3]->Fill(W_org*W_org,rate_GHz);
      W_hist[3]->Fill(W_org,rate_GHz);
      r_hist[3]->Fill(r/1000,rate_GHz);
      r_eweight[3]->Fill(r/1000,kineE*rate_GHz);
      r_asymtot[3]->Fill(r/1000,rate_GHz*asym_org);
      asym_hist[3]->Fill(asym_org,rate_GHz);
      depol_hist[3]->Fill(Din,rate_GHz);
      asymd_hist[3]->Fill(asymd_org,rate_GHz);
      asymEep_hist[3]->Fill(asymEep_org,rate_GHz);
      asymIep_hist[3]->Fill(asymIep_org,rate_GHz);
      z0_hist[3]->Fill(z0,rate_GHz);
      if (TMath::Abs(phi_wrap)<pi/28.) {
        eweight_red_total_asym->Fill(fillval,kineE*asym_trans*totXS);
        eweight_red_total_sample->Fill(fillval,kineE*totXS);
        std_red_total_asym->Fill(fillval,asym_trans*totXS);
        std_red_total_sample->Fill(fillval,totXS);
        distro_red->Fill(costheta_cms_org,phi_org, rate_GHz);
        theta_cm[0]->Fill(acos(costheta_cms_org)*180/TMath::Pi(),rate_GHz);
        Eprime_vert_hist[0]->Fill(kineE_org/1000,rate_GHz);
        Eprime_det_hist[0]->Fill(kineE/1000,rate_GHz);
        theta_lab_hist[0]->Fill(theta_lab*1000,rate_GHz);
        phi_wrap_hist[0]->Fill(phi_wrap*180/TMath::Pi(),rate_GHz);
        Q2_hist[0]->Fill(Q2_org,rate_GHz);
        Q2_eweight_hist[0]->Fill(Q2_org,kineE*rate_GHz);
        W2_Q2weight_hist[0]->Fill(W_org*W_org,Q2_org*rate_GHz);
        W2_hist[0]->Fill(W_org*W_org,rate_GHz);
        W_hist[0]->Fill(W_org,rate_GHz);
        r_hist[0]->Fill(r/1000,rate_GHz);
        r_eweight[0]->Fill(r/1000,kineE*rate_GHz);
        r_asymtot[0]->Fill(r/1000,rate_GHz*asym_org);
        asym_hist[0]->Fill(asym_org,rate_GHz);
        depol_hist[0]->Fill(Din,rate_GHz);
        asymd_hist[0]->Fill(asymd_org,rate_GHz);
        asymEep_hist[0]->Fill(asymEep_org,rate_GHz);
        asymIep_hist[0]->Fill(asymIep_org,rate_GHz);
        z0_hist[0]->Fill(z0,rate_GHz);
      } else {
        if (TMath::Abs(phi_wrap)>pi/28.&&TMath::Abs(phi_wrap)<3*pi/28.) {
          eweight_green_total_asym->Fill(fillval,kineE*asym_trans*totXS);
          eweight_green_total_sample->Fill(fillval,kineE*totXS);
          std_green_total_asym->Fill(fillval,asym_trans*totXS);
          std_green_total_sample->Fill(fillval,totXS);
          distro_green->Fill(costheta_cms_org,phi_org);
          theta_cm[1]->Fill(acos(costheta_cms_org)*180/TMath::Pi(),rate_GHz);
          Eprime_vert_hist[1]->Fill(kineE_org/1000,rate_GHz);
          Eprime_det_hist[1]->Fill(kineE/1000,rate_GHz);
          theta_lab_hist[1]->Fill(theta_lab*1000,rate_GHz);
          phi_wrap_hist[1]->Fill(phi_wrap*180/TMath::Pi(),rate_GHz);
          Q2_hist[1]->Fill(Q2_org,rate_GHz);
          Q2_eweight_hist[1]->Fill(Q2_org,kineE*rate_GHz);
          W2_Q2weight_hist[1]->Fill(W_org*W_org,Q2_org*rate_GHz);
          W2_hist[1]->Fill(W_org*W_org,rate_GHz);
          W_hist[1]->Fill(W_org,rate_GHz);
          r_hist[1]->Fill(r/1000,rate_GHz);
          r_eweight[1]->Fill(r/1000,kineE*rate_GHz);
          asym_hist[1]->Fill(asym_org,rate_GHz);
          depol_hist[1]->Fill(Din,rate_GHz);
          asymd_hist[1]->Fill(asymd_org,rate_GHz);
          asymEep_hist[1]->Fill(asymEep_org,rate_GHz);
          asymIep_hist[1]->Fill(asymIep_org,rate_GHz);
          r_asymtot[1]->Fill(r/1000,rate_GHz*asym_org);
          z0_hist[1]->Fill(z0,rate_GHz);
        } else {
          if (TMath::Abs(phi_wrap)>3*pi/28.) {
            eweight_blue_total_asym->Fill(fillval,kineE*asym_trans*totXS);
            eweight_blue_total_sample->Fill(fillval,kineE*totXS);
            std_blue_total_asym->Fill(fillval,asym_trans*totXS);
            std_blue_total_sample->Fill(fillval,totXS);
            distro_blue->Fill(costheta_cms_org,phi_org, rate_GHz);
            theta_cm[2]->Fill(acos(costheta_cms_org)*180/TMath::Pi(),rate_GHz);
            Eprime_vert_hist[2]->Fill(kineE_org/1000,rate_GHz);
            Eprime_det_hist[2]->Fill(kineE/1000,rate_GHz);
            theta_lab_hist[2]->Fill(theta_lab*1000,rate_GHz);
            phi_wrap_hist[2]->Fill(phi_wrap*180/TMath::Pi(),rate_GHz);
            Q2_hist[2]->Fill(Q2_org,rate_GHz);
            Q2_eweight_hist[2]->Fill(Q2_org,kineE*rate_GHz);
            W2_Q2weight_hist[2]->Fill(W_org*W_org,Q2_org*rate_GHz);
            W2_hist[2]->Fill(W_org*W_org,rate_GHz);
            W_hist[2]->Fill(W_org,rate_GHz);
            r_hist[2]->Fill(r/1000,rate_GHz);
            r_eweight[2]->Fill(r/1000,kineE*rate_GHz);
            r_asymtot[2]->Fill(r/1000,rate_GHz*asym_org);
            asym_hist[2]->Fill(asym_org,rate_GHz);
            depol_hist[2]->Fill(Din,rate_GHz);
            asymd_hist[2]->Fill(asymd_org,rate_GHz);
            asymEep_hist[2]->Fill(asymEep_org,rate_GHz);
            asymIep_hist[2]->Fill(asymIep_org,rate_GHz);
            z0_hist[2]->Fill(z0,rate_GHz);
          }
        }
      }
      test->Fill(phi_det,fillval);
      if (kineE0>10999.999) {
        paperfig->Fill(costheta_cms_org,asym_trans/sin(phi_org));
        paperfig2->Fill(costheta_cms_org,phi_org,asym_trans);
      }
    }
    
  }

  if (0) {
    rate_total->Draw();
    snprintf(plotname,200,"%sratetotal.gif",namestem);
    canvas.Print(plotname);
    num_events->Draw();
    snprintf(plotname,200,"%snumevents.gif",namestem);
    canvas.Print(plotname);

    std_red_ave_asym->Divide(std_red_total_asym,std_red_total_sample);
    std_green_ave_asym->Divide(std_green_total_asym,std_green_total_sample);
    std_blue_ave_asym->Divide(std_blue_total_asym,std_blue_total_sample);
    eweight_red_ave_asym->Divide(eweight_red_total_asym,eweight_red_total_sample);
    eweight_green_ave_asym->Divide(eweight_green_total_asym,eweight_green_total_sample);
    eweight_blue_ave_asym->Divide(eweight_blue_total_asym,eweight_blue_total_sample);
    Double_t errorfactor;
    // Set error bars expected for 1 hour of running at 85 uA
    for (Int_t i = 1; i<=num_detectors; i++) {
      errorfactor = 1000000/sqrt(num_events->GetBinContent(i));
      error->SetBinContent(i,errorfactor);
      if (std_red_ave_asym->GetBinContent(i)!=0) std_red_ave_asym->SetBinError(i,errorfactor);
      if (std_green_ave_asym->GetBinContent(i)!=0) std_green_ave_asym->SetBinError(i,errorfactor);
      if (std_blue_ave_asym->GetBinContent(i)!=0) std_blue_ave_asym->SetBinError(i,errorfactor);
      if (eweight_red_ave_asym->GetBinContent(i)!=0) eweight_red_ave_asym->SetBinError(i,errorfactor);
      if (eweight_green_ave_asym->GetBinContent(i)!=0) eweight_green_ave_asym->SetBinError(i,errorfactor);
      if (eweight_blue_ave_asym->GetBinContent(i)!=0) eweight_blue_ave_asym->SetBinError(i,errorfactor);
    }
    error->Draw();
    snprintf(plotname,200,"%serror.gif",namestem);
    canvas.Print(plotname);

    TF1 *fa1 = new TF1("fa1","1.5*sin(2*TMath::Pi()*(x+0.5)/28.)",0,28);
    fa1->SetLineWidth(2);
    snprintf(plotname,200,"%sstd_color_dettotal.gif",namestem);
    std_red_total_asym->Draw();
    std_green_total_asym->Draw("same");
    std_blue_total_asym->Draw("same");
    canvas.Print(plotname);
    snprintf(plotname,200,"%sstd_color_detsample.gif",namestem);
    std_red_total_sample->Draw();
    std_green_total_sample->Draw("same");
    std_blue_total_sample->Draw("same");
    canvas.Print(plotname);
    snprintf(plotname,200,"%sstd_color_detave.gif",namestem);
    std_blue_ave_asym->Draw();
    std_red_ave_asym->Draw("same");
    std_green_ave_asym->Draw("same");
    fa1->Draw("same");
    std_ave_asym->Divide(std_total_asym,std_total_sample);
    //  std_ave_asym->Draw("same");
    canvas.Print(plotname);



    TF1 *fa2 = new TF1("fa2","-2*sin(2*TMath::Pi()*(x+0.5)/28.)",0,28);
    fa2->SetLineWidth(2);

    snprintf(plotname,200,"%seweight_color_dettotal.gif",namestem);
    eweight_red_total_asym->Draw();
    eweight_green_total_asym->Draw("same");
    eweight_blue_total_asym->Draw("same");
    canvas.Print(plotname);
    snprintf(plotname,200,"%seweight_color_detsample.gif",namestem);
    eweight_red_total_sample->Draw();
    eweight_green_total_sample->Draw("same");
    eweight_blue_total_sample->Draw("same");
    canvas.Print(plotname);
    snprintf(plotname,200,"%seweight_color_detave.gif",namestem);
    eweight_blue_ave_asym->Draw();
    eweight_red_ave_asym->Draw("same");
    eweight_green_ave_asym->Draw("same");
    eweight_ave_asym->Divide(eweight_total_asym,eweight_total_sample);
    //  eweight_ave_asym->Draw("same");
    fa2->Draw("same");
    canvas.Print(plotname);


    test->Draw("colz");
    snprintf(plotname,200,"%stest.gif",namestem);
    canvas.Print(plotname);
  }

  // *** Draw the 1D histograms
  snprintf(plotname,200,"%stheta_cm%s",outdir,plottype);
  PrintHist(theta_cm,plotname,4,0);
  snprintf(plotname,200,"%sEprime_vert%s",outdir,plottype);
  PrintHist(Eprime_vert_hist,plotname,4,0);
  snprintf(plotname,200,"%sEprime_det%s",outdir,plottype);
  PrintHist(Eprime_det_hist,plotname,4,0);
  snprintf(plotname,200,"%stheta_lab%s",outdir,plottype);
  PrintHist(theta_lab_hist,plotname,4,0);
  snprintf(plotname,200,"%sphi_wrap%s",outdir,plottype);
  PrintHist(phi_wrap_hist,plotname);
  snprintf(plotname,200,"%sQ2%s",outdir,plottype);
  PrintHist(Q2_hist,plotname,4,0);
  snprintf(plotname,200,"%sQ2_eweight%s",outdir,plottype);
  PrintHist(Q2_eweight_hist,plotname,4,0);
  snprintf(plotname,200,"%sW2_Q2weight%s",outdir,plottype);
  PrintHist(W2_Q2weight_hist,plotname,4,0);
  snprintf(plotname,200,"%sW2%s",outdir,plottype);
  PrintHist(W2_hist,plotname,4,0);
  snprintf(plotname,200,"%sW%s",outdir,plottype);
  PrintHist(W_hist,plotname,4,0);
  snprintf(plotname,200,"%sr%s",outdir,plottype);
  PrintHist(r_hist,plotname,5,0);
  for (Int_t i=0; i<numhistos; i++) {
    r_asymave[i]->Divide(r_asymtot[i],r_hist[i]);
  }
  snprintf(plotname,200,"%sr_asymtot%s",outdir,plottype);
  PrintHist(r_asymtot,plotname,5,0);
  snprintf(plotname,200,"%sr_asymave%s",outdir,plottype);
  PrintHist(r_asymave,plotname,5,0);
  snprintf(plotname,200,"%sasym%s",outdir,plottype);
  PrintHist(asym_hist,plotname,4,0);
  snprintf(plotname,200,"%sdepol%s",outdir,plottype);
  PrintHist(depol_hist,plotname,5,1);
  snprintf(plotname,200,"%sasymd%s",outdir,plottype);
  PrintHist(asymd_hist,plotname,4,0);

  snprintf(plotname,200,"%sasymEep%s",outdir,plottype);
  PrintHist(asymEep_hist,plotname,4,0);
  snprintf(plotname,200,"%sasymIep%s",outdir,plottype);
  PrintHist(asymIep_hist,plotname,4,0);
  snprintf(plotname,200,"%sz0%s",outdir,plottype);
  PrintHist(z0_hist,plotname);
  // *** Draw the 2D histograms
  snprintf(plotname,200,"%s2Dr_phi%s",outdir,plottype);
  PrintHist2D(r_phi_hist,plotname);
  snprintf(plotname,200,"%s2D_r_asym%s",outdir,plottype);
  PrintHist2D(r_asym_hist,plotname);


  distro_r->Draw();
  snprintf(plotname,200,"%sdistro_r.gif",outdir);
  canvas.SetLogy(1);
  canvas.Print(plotname);
  if (0) {
    paperfig->Draw();
    snprintf(plotname,200,"%spaperfig.gif",outdir);
    canvas.Print(plotname);
    paperfig2->Project3DProfile("yx")->Draw("colz");
    snprintf(plotname,200,"%sasym_mag.gif",outdir);
    canvas.Print(plotname);
    distro_red->Draw("contourz");
    snprintf(plotname,200,"%sdistro_red.gif",outdir);
    canvas.Print(plotname);
    distro_green->Draw("contourz");
    snprintf(plotname,200,"%sdistro_green.gif",outdir);
    canvas.Print(plotname);
    distro_blue->Draw("contourz");
    snprintf(plotname,200,"%sdistro_blue.gif",outdir);
    canvas.Print(plotname);

    distro_rphi->Draw("contourz");
    snprintf(plotname,200,"%sdistro_rphi.gif",outdir);
    canvas.Print(plotname);
    distro_r_evert->Draw("contourz");
    snprintf(plotname,200,"%sdistro_r_evert.gif",outdir);
    canvas.Print(plotname);
    distro_r_edet->Draw("contourz");
    snprintf(plotname,200,"%sdistro_r_edet.gif",outdir);
    canvas.Print(plotname);
  }

  printf("Total number of events       = %8li\n", totalevents);
  printf("number not passing rcut      = %8li\n", numevcut_rcut);
  printf("number not passing thetacut  = %8li\n", numevcut_thetacut);
  printf("number not passing volumecut = %8li\n", numevcut_volumecut);
  printf("number not passing kinecut   = %8li\n", numevcut_kinecut);
  printf("number passing all cuts      = %8li\n", numevallpass);

  printf("Total rate is %5g GHz\n",distro_r->Integral());
  printf("Rate in detector zone is %5g GHz\n",distro_r->Integral(176,200));
  printf("Red   detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[0]->Integral(),r_eweight[0]->Integral());
  printf("Green detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[1]->Integral(),r_eweight[1]->Integral());
  printf("Blue  detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[2]->Integral(),r_eweight[2]->Integral());
  printf("Total detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[3]->Integral(),r_eweight[3]->Integral());

  printf("Rates between 0.66 and 0.82 m\n");
  printf("Red   detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[0]->Integral(12,44),r_eweight[0]->Integral(12,44));
  printf("Green detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[1]->Integral(12,44),r_eweight[1]->Integral(12,44));
  printf("Blue  detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[2]->Integral(12,44),r_eweight[2]->Integral(12,44));
  printf("Total detector rate: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",r_hist[3]->Integral(12,44),r_eweight[3]->Integral(12,44));
  printf("All detector plane rate between 0.66 and 0.82 m: %5g GHz,   energy weighted rate: %5g MeV*GHz\n",
       r_hist[4]->Integral(12,44),r_eweight[4]->Integral(12,44));

  printf("Red   ave asymmetry: %5g ppb\n",asym_hist[0]->GetMean());
  printf("Green ave asymmetry: %5g ppb\n",asym_hist[1]->GetMean());
  printf("Blue  ave asymmetry: %5g ppb\n",asym_hist[2]->GetMean());
  printf("Total ave asymmetry: %5g ppb\n",asym_hist[3]->GetMean());
  printf("{%5g ppb, %5g ppb, %5g ppb, %5g ppb}\n",
       asym_hist[0]->GetMean(), asym_hist[1]->GetMean(), asym_hist[2]->GetMean(),asym_hist[3]->GetMean());

  printf("Red   ave depol factor: %5.3f\n",depol_hist[0]->GetMean());
  printf("Green ave depol factor: %5.3f\n",depol_hist[1]->GetMean());
  printf("Blue  ave depol factor: %5.3f\n",depol_hist[2]->GetMean());
  printf("Total ave depol factor: %5.3f\n",depol_hist[3]->GetMean());
  printf("{%5.3f, %5.3f, %5.3f, %5.3f}\n",
       depol_hist[0]->GetMean(), depol_hist[1]->GetMean(), depol_hist[2]->GetMean(),depol_hist[3]->GetMean());

  printf("Red   ave asymmetry: %5g ppb\n",asymd_hist[0]->GetMean());
  printf("Green ave asymmetry: %5g ppb\n",asymd_hist[1]->GetMean());
  printf("Blue  ave asymmetry: %5g ppb\n",asymd_hist[2]->GetMean());
  printf("Total ave asymmetry: %5g ppb\n",asymd_hist[3]->GetMean());
  printf("{%5g ppb, %5g ppb, %5g ppb, %5g ppb}\n",
       asymd_hist[0]->GetMean(), asymd_hist[1]->GetMean(), asymd_hist[2]->GetMean(),asymd_hist[3]->GetMean());

  printf("Red   ave asymmetry: %5g ppb\n",asymEep_hist[0]->GetMean());
  printf("Green ave asymmetry: %5g ppb\n",asymEep_hist[1]->GetMean());
  printf("Blue  ave asymmetry: %5g ppb\n",asymEep_hist[2]->GetMean());
  printf("Total ave asymmetry: %5g ppb\n",asymEep_hist[3]->GetMean());
  printf("{%5g ppb, %5g ppb, %5g ppb, %5g ppb}\n",
       asymEep_hist[0]->GetMean(), asymEep_hist[1]->GetMean(), asymEep_hist[2]->GetMean(),asymEep_hist[3]->GetMean());

  printf("Red   ave asymmetry: %5g ppb\n",asymIep_hist[0]->GetMean());
  printf("Green ave asymmetry: %5g ppb\n",asymIep_hist[1]->GetMean());
  printf("Blue  ave asymmetry: %5g ppb\n",asymIep_hist[2]->GetMean());
  printf("Total ave asymmetry: %5g ppb\n",asymIep_hist[3]->GetMean());
  printf("{%5g ppb, %5g ppb, %5g ppb, %5g ppb}\n",
       asymIep_hist[0]->GetMean(), asymIep_hist[1]->GetMean(), asymIep_hist[2]->GetMean(),asymIep_hist[3]->GetMean());

  // Write some output to file
  char rootoutfilename[256];
  GetPrivateProfileString ("General", "rootoutfile", "defaultout.root", rootoutfilename, 256, inifilename);
  TFile *outfile=new TFile(rootoutfilename,"update");
  if (outfile->IsZombie()) {
    cerr << "\nError opening file " <<  rootoutfilename << endl << endl;
    exit(-1);
  }
  for (Int_t i=0; i<numhistos; i++) {
    r_hist[i]->Write("", TObject::kOverwrite);
    Eprime_det_hist[i]->Write("", TObject::kOverwrite);
    Eprime_vert_hist[i]->Write("", TObject::kOverwrite);
    theta_cm[i]->Write("", TObject::kOverwrite);
    theta_lab_hist[i]->Write("", TObject::kOverwrite);
    phi_wrap_hist[i]->Write("", TObject::kOverwrite);
    Q2_hist[i]->Write("", TObject::kOverwrite);
    Q2_eweight_hist[i]->Write("", TObject::kOverwrite);
    W_hist[i]->Write("", TObject::kOverwrite);
    W2_hist[i]->Write("", TObject::kOverwrite);
    W2_Q2weight_hist[i]->Write("", TObject::kOverwrite);
    r_asymave[i]->Write("", TObject::kOverwrite);
    asym_hist[i]->Write("", TObject::kOverwrite);
    asymd_hist[i]->Write("", TObject::kOverwrite);
    depol_hist[i]->Write("", TObject::kOverwrite);
    asymEep_hist[i]->Write("", TObject::kOverwrite);
    asymIep_hist[i]->Write("", TObject::kOverwrite);


  }
  outfile->Close();

}








int main(int argc,char **argv)
{

  if (argc < 2) {
    cerr << "\nToo few arguements!\nuseage: mollerClass inputfile\n\n";
    exit(-1);
    //snprintf(inifilename,100,"fitsignal.ini");
  }



  for (int i=0; i<argc; i++) {
    cout << argv[i] << " ";
  }
  cout << "\n";
  if (argc < 2) {
    cerr << "\nuseage:  mollerClass\n\n";
    exit(-1);
  } else {
    snprintf(inifilename,256,argv[1]);
    TGaxis::SetMaxDigits(3); 
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(kFALSE);
    gStyle->SetPalette(1,0);
    gStyle->SetTitleBorderSize(0);
    gStyle->SetCanvasBorderMode(0);




    TChain *geant = new TChain("geant");


    char filenamestring[256], lookupstring[256];
    //snprintf(filenamestring,256,"start");
    Int_t counter=1;
    snprintf(lookupstring,256,"file%i",counter);
//    printf("%s\n",lookupstring);
    GetPrivateProfileString ("Filenames", lookupstring, "0", filenamestring, 256, inifilename);
//    printf("%i %s  %s\n",counter,lookupstring,filenamestring);
//    snprintf(lookupstring,256,"file%i",2);
//    GetPrivateProfileString ("Filenames", lookupstring, "0", filenamestring, 256, inifilename);
//    printf("%i %s  %s\n",counter,lookupstring,filenamestring);
//    snprintf(lookupstring,256,"file%i",3);
//    GetPrivateProfileString ("Filenames", lookupstring, "0", filenamestring, 256, inifilename);
//    printf("%i %s  %s\n",counter,lookupstring,filenamestring);
//    snprintf(lookupstring,256,"file%i",4);
//    GetPrivateProfileString ("Filenames", lookupstring, "0", filenamestring, 256, inifilename);
//    printf("%i %s  %s\n",counter,lookupstring,filenamestring);

    while (strcmp(filenamestring,"0")!=0) {
      printf("%i  %s\n",counter,filenamestring);
      geant->Add(filenamestring);
//      printf("%i %s  %s\n",counter,lookupstring,filenamestring);
      counter++;
//      printf("%i\n",counter);
      snprintf(lookupstring,256,"file%i",counter);
//      printf("%i %s  %s\n",counter,lookupstring,filenamestring);
      GetPrivateProfileString ("Filenames", lookupstring, "0", filenamestring, 256, inifilename);
    }

    char num_thrown_string[50];
    GetPrivateProfileString ("General", "num_thrown", "0", num_thrown_string, 50, inifilename);
    num_thrown = atoi(num_thrown_string);
    if (num_thrown<=0) {
      printf("num_thrown cannot equal %i\n",num_thrown);
      exit(-1);
    }

    eventtype=GetPrivateProfileInt("General", "eventtype", -1, inifilename);
    switch (eventtype) {
    case 0:
      printf("Analysing MOLLER events.\n");
      num_thrown=2*num_thrown;
      break;
    case 1: 
      printf("Analysing ELASTIC ep events.\n");
      break;
    case 2: 
      printf("Analysing INelastic ep events.\n");
      break;
    default:
      printf("eventtype not set correctly\n");
      exit(-1);
    }
    printf("num_thrown = %i\n",num_thrown);


    mollerClass* mollerClassobj = new mollerClass(geant);
    mollerClassobj->Loop();

    return 1;
  }
}



/* emacs
 * Local Variables:
 * mode:C++
 * mode:font-lock
 * c-file-style: "stroustrup"
 * tab-width: 4
 * End:
 */