///////////////////////////////////////////////////////////
//
//  Puts some selected info from each run in a runlist
//  into a rootfile.  
//
//needs the following files: 
// 1. runlist.txt (or somefile with in runlist.txt format with the list of runs)
// 2. slopes/avgSlopes_hydrogen_dither_slug#.res
// 3. takes range0 as the argument
// 4. make sure to specify the correct cuts

#include <iostream.h>
#include <stdlib.h>
#include "../src/TaFileName.hh"
#include <string>

int mkSumTreeH(TString slugrangefname = "") {

  typedef UInt_t Long64_t;

  // Get Run list
  const int mxNRuns = 1001;
  UInt_t runlist[mxNRuns];
  UInt_t sluglist[mxNRuns];

  Int_t firstslug=0;
  Int_t lastslug=-1;
  
  ifstream sluglimitfile("./"+slugrangefname);
  if(!sluglimitfile.good()) {
    cout << "Failure to read slug limit file" << endl;
    // firstslug =0;
    // lastslug = 10000000;
    //    cout << "doing ALL slugs" << endl;
    return;
  }
  sluglimitfile >> firstslug >> lastslug;
  if (lastslug<firstslug) {
    cout << "problem with slug limit file" << endl;
    return;
  }

  cout << "first slug = " << firstslug << endl;
  cout << "last slug = " << lastslug << endl;
 
  // plug in some run list here...
  TString runlistStr = "./runlist.txt";

  TaRunlist dblist(runlistStr.Data());
  if(!dblist.GoodRunlist()) {
    cout << "Bad runlist" << endl;
    return;
  }

  // Fill the runlist with the runs
//   Int_t nruns = dblist.GetNumberOfRuns();
  
//   for (Int_t ir = 0; ir<nruns; ir++) {
//     runlist[ir] = dblist.GetRunByIndex(ir);
//     sluglist[ir] = dblist.GetSlug(runlist[ir]);
//     cout << runlist[ir] << "  " << sluglist[ir] << endl;
//   }
//   cout << "largest run number = " << runlist[nruns-1] << endl;

  // Fill the runlist with the runs in these slugs
  Int_t total_nruns = dblist.GetNumberOfRuns();
  Int_t nruns;
  
  Int_t thisrun;
  Int_t thisslug;
  for (Int_t ir = 0; ir<total_nruns; ir++) {
    thisrun = dblist.GetRunByIndex(ir);
    thisslug = dblist.GetSlug(thisrun);
    if (thisslug>=firstslug && thisslug<=lastslug) {
      runlist[nruns] = thisrun;
      sluglist[nruns] = thisslug;
      cout << runlist[nruns] << "  " << sluglist[nruns] << endl;
      nruns++;
    }
  }
  cout << "largest run number = " << runlist[nruns-1] << endl;

  //  gROOT->ProcessLine(".q");  
  //return;

  //**********************************************
  // standard procedure: 
  //
  // create tree and define leaves in the tree.
  //
  // loop over runs.
  //      make a histo, get appropriate mean, rms, etc info
  //      then delete histo.  Store data is word pointing to a tree leaf.
  //
  //      When done extracting data, fill tree.

  // end loop over runs
  // write and close tree
  //**********************************************

  //
  // Create tree and define desired data words
  //     

  // name a root file in compliance with 'pan' convention, 27Dec09, rupesh
  string slugnumstr = Form("%.2d%.2d",firstslug,lastslug);
  Int_t slugnum=atoi(slugnumstr.c_str());    // string to int
  cout << slugnumstr << "   " << slugnum << endl;
  TaFileName::Setup(slugnum,"runsummary");
   // TaFileName::Setup(1234,"runsummary");
  TString outrootfilestr = (TaFileName ("root")).Tstring();

  //  gROOT->ProcessLine(".qqqqqqqqqqqq");  
  //  return;

  //  TString outrootfilestr = Form("summaryrun/testsummary_slug%d_slug%d.root",
  //  				firstslug,lastslug);
  TFile *fRootFile = new TFile(outrootfilestr,"RECREATE");
  TTree *stats = new TTree("stats","Run statistics");
  stats->SetMarkerStyle(20);
  stats->SetMarkerSize(1.);
  stats->SetMarkerColor(2);

  const UInt_t ndet=3;
  TString det[ndet] = {"det1","det2","det_all"}; 
  //  Int_t icdet[ndet] = {1,1,2,2,1,2,3,3,3};
  Int_t icdet[ndet] = {1,2,3};

  const UInt_t nblumi=8; 
  TString blumi[nblumi] = {"blumi1","blumi2","blumi3","blumi4","blumi5","blumi6","blumi7","blumi8"}; 
  //  const UInt_t nflumi=2; 
  //  TString flumi[nflumi] = {"flumi1","flumi2"};
  const UInt_t nbpm=5; 
  TString bpm[nbpm] = {"4bx", "4by", "4ax", "4ay", "12x"}; 
  const UInt_t nbcm=4;
  TString bcm[nbcm] = {"1","2","5","10"};

  const UInt_t nEpics=4;
  TString strEpics[nEpics] = {"IGL1I00HALLAMODE","IGL1I00HALLBMODE",
			      "IGLdac2_G2Ch3Pos-IGLdac2_G2Ch4Neg",
			      "IGLdac3_ao_7"};
  TString nameEpics[nEpics] = {"Laser_Mode_A","Laser_Mode_B",
			      "Pita_OFfset",
			      "IA_Setpoint"};
  Double_t EpicsData[nEpics];

  Int_t runnumber;
  Int_t slowsign; //IHWP +1 for out, -1 for in
  UInt_t slug;
  Int_t ok_stats[ndet]; // Useful for runs without statistics 
  Double_t udet[ndet],udeterr[ndet],udetrms[ndet];
  Double_t normdet[ndet],normdeterr[ndet],normdetrms[ndet];
  Double_t rawbcm[nbcm], rawbcmerr[nbcm], rawbcmrms[nbcm];
  Double_t rawdet[ndet], rawdeterr[ndet], rawdetrms[ndet];
  Double_t regdet[ndet], regdeterr[ndet], regdetrms[ndet];
  Double_t ditdet[ndet], ditdeterr[ndet], ditdetrms[ndet];
  Double_t ublumi[nblumi],ublumierr[nblumi],ublumirms[nblumi];
  Double_t normblumi[nblumi],normblumierr[nblumi],normblumirms[nblumi];
  Double_t regblumi[nblumi], regblumierr[nblumi], regblumirms[nblumi];
  Double_t ditblumi[nblumi], ditblumierr[nblumi], ditblumirms[nblumi];  
  //  Double_t uflumi[nflumi],uflumierr[nflumi],uflumirms[nflumi];
  // Double_t normflumi[nflumi],normflumierr[nflumi],normflumirms[nflumi];
  // Double_t regflumi[nflumi], regflumierr[nflumi], regflumirms[nflumi];
  //Double_t ditflumi[nflumi], ditflumierr[nflumi], ditflumirms[nflumi];
  Double_t DD12,DD12rms,DD12err;
  Double_t bpmdat[nbpm],bpmrms[nbpm],bpmerr[nbpm];
  //  Double_t injbpmdat[ninj], injbpmerr[ninj],injbpmrms[ninj];
  Double_t bpmdiffdat[nbpm],bpmdiffrms[nbpm],bpmdifferr[nbpm];
  //  Double_t injbpmdiffdat[ninj], injbpmdifferr[ninj],injbpmdiffrms[ninj];
  Double_t bcmasymdat[nbcm],bcmasymrms[nbcm],bcmasymerr[nbcm];
  Int_t numevents;


  TH1F* hdet[ndet];
  TH1F* hdet_dit[ndet];
  TH1F* hdet_reg[ndet];
  TString hdetnam[ndet];
  TString hdetregnam[ndet];
  TString hdetditnam[ndet];

  TH2F* hdetlin_reg[ndet];
  TString hdetlinregnam[ndet];
  TH2F* hdetlin_dit[ndet];
  TString hdetlinditnam[ndet];
  TH1F* hDD;
  TString hDDnam;
  TH2F* hbcmlin;
  TString hbcmlinnam;


  stats->Branch("run",&runnumber,"run/I");
  stats->Branch("sign",&slowsign,"sign/I");
  stats->Branch("slug",&slug,"slug/I");
  stats->Branch("ok_statsL",&ok_stats[0],"ok_statsL/I");
  stats->Branch("ok_statsR",&ok_stats[1],"ok_statsR/I");
  stats->Branch("nevents",&numevents,"nevents/I");

  TString bname;

 
  for(UInt_t ibcm=0; ibcm<nbcm; ibcm++) { 
    bname = "bcm" +bcm[ibcm];
    stats->Branch(bname,&rawbcm[ibcm],bname+"/D");
    bname = "bcm" +bcm[ibcm]+"_err";
    stats->Branch(bname,&rawbcmerr[ibcm],bname+"/D");
    bname = "bcm" +bcm[ibcm]+"_rms";
    stats->Branch(bname,&rawbcmrms[ibcm],bname+"/D"); 
    
  } 
   

  for(UInt_t idet=0; idet<ndet; idet++) {
    bname = det[idet];
    stats->Branch(bname,&udet[idet],bname+"/D");
    bname = det[idet]+"_err";
    stats->Branch(bname,&udeterr[idet],bname+"/D");
    bname = det[idet]+"_rms";
    stats->Branch(bname,&udetrms[idet],bname+"/D");

    bname = det[idet]+ "_n";
    stats->Branch(bname,&normdet[idet],bname+"/D");
    bname = det[idet]+ "_n"+"_err";
    stats->Branch(bname,&normdeterr[idet],bname+"/D");
    bname = det[idet]+ "_n"+"_rms"; 
    stats->Branch(bname,&normdetrms[idet],bname+"/D");

    bname = "asym_n_" +det[idet]; 
    stats->Branch(bname,&rawdet[idet],bname+"/D");
    bname = "asym_n_" + det[idet] +"_err";
    stats->Branch(bname,&rawdeterr[idet],bname+"/D");
    bname = "asym_n_" + det[idet] +"_rms";
    stats->Branch(bname,&rawdetrms[idet],bname+"/D");
    
    bname = "reg_asym_n_" + det[idet];
    stats->Branch(bname,&regdet[idet],bname+"/D");
    bname = "reg_asym_n_" + det[idet]+ "_err";
    stats->Branch(bname,&regdeterr[idet],bname+"/D");
    bname = "reg_asym_n_" + det[idet]+ "_rms";
    stats->Branch(bname,&regdetrms[idet],bname+"/D");
 
    bname = "dit_asym_n_" + det[idet];
    stats->Branch(bname,&ditdet[idet],bname+"/D");
    bname = "dit_asym_n_" + det[idet]+ "_err";
    stats->Branch(bname,&ditdeterr[idet],bname+"/D");
    bname = "dit_asym_n_" + det[idet]+ "_rms";
    stats->Branch(bname,&ditdetrms[idet],bname+"/D");
  }



//   for(UInt_t iblumi=0;iblumi<nblumi;iblumi++) { 
//     bname = blumi[iblumi];
//     stats->Branch(bname,&ublumi[iblumi],bname+"/D");
//     bname = blumi[iblumi]+"_err";
//     stats->Branch(bname,&ublumierr[iblumi],bname+"/D");
//     bname = blumi[iblumi]+ "_rms";
//     stats->Branch(bname,&ublumirms[iblumi],bname+"/D");
 

//     bname = blumi[iblumi]+ "_n";
//     stats->Branch(bname,&normblumi[iblumi],bname+"/D");
//     bname = blumi[iblumi]+ "_n"+"_err";
//     stats->Branch(bname,&normblumierr[iblumi],bname+"/D");
//     bname = blumi[iblumi]+ "_n"+"_rms";
//     stats->Branch(bname,&normblumirms[iblumi],bname+"/D"); 
 
//     bname = "reg_asym_n_" + blumi[iblumi];
//     stats->Branch(bname,&regblumi[iblumi],bname+"/D");
//     bname = "reg_asym_n_" + blumi[iblumi]+ "_err";
//     stats->Branch(bname,&regblumierr[iblumi],bname+"/D");
//     bname = "reg_asym_n_" + blumi[iblumi]+ "_rms";
//     stats->Branch(bname,&regblumirms[iblumi],bname+"/D"); 

//     bname = "dit_asym_n_" + blumi[iblumi];
//     stats->Branch(bname,&ditblumi[iblumi],bname+"/D");
//     bname = "dit_asym_n_" + blumi[iblumi]+ "_err";
//     stats->Branch(bname,&ditblumierr[iblumi],bname+"/D");
//     bname = "dit_asym_n_" + blumi[iblumi]+ "_rms";
//     stats->Branch(bname,&ditblumirms[iblumi],bname+"/D");

//   }

//  for(UInt_t iflumi=0;iflumi<nflumi;iflumi++) {
//     bname = flumi[iflumi];
//     stats->Branch(bname,&uflumi[iflumi],bname+"/D");
//     bname = flumi[iflumi]+"_err";
//     stats->Branch(bname,&uflumierr[iflumi],bname+"/D");
//     bname = flumi[iflumi]+ "_rms";
//     stats->Branch(bname,&uflumirms[iflumi],bname+"/D");
                                                                                
                                                                                
//     bname = flumi[iflumi]+ "_n";
//     stats->Branch(bname,&normflumi[iflumi],bname+"/D");
//     bname = flumi[iflumi]+ "_n"+"_err";
//     stats->Branch(bname,&normflumierr[iflumi],bname+"/D");
//     bname = flumi[iflumi]+ "_n"+"_rms";
//     stats->Branch(bname,&normflumirms[iflumi],bname+"/D");
                                                                                
//     bname = "reg_asym_n_" + flumi[iflumi];
//     stats->Branch(bname,&regflumi[iflumi],bname+"/D");
//     bname = "reg_asym_n_" + flumi[iflumi]+ "_err";
//     stats->Branch(bname,&regflumierr[iflumi],bname+"/D");
//     bname = "reg_asym_n_" + flumi[iflumi]+ "_rms";
//     stats->Branch(bname,&regflumirms[iflumi],bname+"/D");
    
// //     bname = "dit_asym_n_" + flumi[iflumi];
// //     stats->Branch(bname,&ditflumi[iflumi],bname+"/D");
// //     bname = "dit_asym_n_" + flumi[iflumi]+ "_err";
// //     stats->Branch(bname,&ditflumierr[iflumi],bname+"/D");
// //     bname = "dit_asym_n_" + flumi[iflumi]+ "_rms";
// //     stats->Branch(bname,&ditflumirms[iflumi],bname+"/D");
                                                                                
//  }


// for(UInt_t iinj=0;iinj<ninj;iinj++) { 
//    bname = injbpm[iinj]; 
//   stats->Branch(bname,&injbpmdat[iinj],bname+"/D");
//   bname = injbpm[iinj]+"_err";
//   stats->Branch(bname,&injbpmerr[iinj],bname+"/D");
//   bname = injbpm[iinj]+"_rms";
//   stats->Branch(bname,&injbpmrms[iinj],bname+"/D"); 
 
//   bname = "diff_"+injbpm[iinj];
//   stats->Branch(bname,&injbpmdiffdat[iinj],bname+"/D");
//   bname = "diff"+injbpm[iinj]+"_err";
//   stats->Branch(bname,&injbpmdifferr[iinj],bname+"/D");
//   bname = "diff_"+injbpm[iinj]+"_rms";
//   stats->Branch(bname,&injbpmdiffrms[iinj],bname+"/D");

// }

for(UInt_t imon=0; imon<nbpm; imon++) { 
    bname = "bpm" +bpm[imon]; 
   stats->Branch(bname,&bpmdat[imon],bname+"/D");
   bname = "bpm" + bpm[imon]+ "_err";
   stats->Branch(bname,&bpmerr[imon],bname+"/D");
   bname = "bpm" + bpm[imon]+ "_rms";
   stats->Branch(bname,&bpmrms[imon],bname+"/D"); 

    bname = "diff_bpm" + bpm[imon];
    stats->Branch(bname,&bpmdiffdat[imon],bname+"/D");
    bname = "diff_bpm" + bpm[imon]+ "_err";
    stats->Branch(bname,&bpmdifferr[imon],bname+"/D");
    bname = "diff_bpm" + bpm[imon]+ "_rms";
    stats->Branch(bname,&bpmdiffrms[imon],bname+"/D");
  }

  for(UInt_t imon=0; imon<nbcm; imon++) {
    bname = "asym_bcm" + bcm[imon];
    stats->Branch(bname,&bcmasymdat[imon],bname+"/D"); 
    bname = "asym_bcm" + bcm[imon]+ "_err"; 
   stats->Branch(bname,&bcmasymerr[imon],bname+"/D");
    bname = "asym_bcm" + bcm[imon]+ "_rms";
   stats->Branch(bname,&bcmasymrms[imon],bname+"/D");
  
  }

  for(UInt_t iepi=0; iepi<nEpics; iepi++) { 
    bname = "E_" +nameEpics[iepi];
    stats->Branch(bname,&EpicsData[iepi],bname+"/D");
  }

  bname = "DD12";
  stats->Branch(bname,&DD12,bname+"/D"); 
  bname = "DD12"; 
  bname += "_err"; 
  stats->Branch(bname,&DD12err,bname+"/D"); 
  bname = "DD12";
  bname += "_rms";
  stats->Branch(bname,&DD12rms,bname+"/D");


  //
  // make histograms
  //  


  for(UInt_t idet=0; idet<ndet; idet++) {

    TString hname = "h_asym_";
    hname += det[idet];
    TString htit = "Asym_n_";
    htit += det[idet];
    hdetnam[idet] = hname;
    hdet[idet]=new TH1F(hname,htit,750,-15000,15000);

    hname = "h_asym_";
    hname += det[idet];
    hname += "_reg";
    htit = "Reg_Asym_n_";
    htit += det[idet];
    hdetregnam[idet] = hname;
    hdet_reg[idet]=new TH1F(hname,htit,500,-10000,10000);

    hname = "h_asym_";
    hname += det[idet];
    hname += "_dit";
    htit = "Dit_Asym_n_";
    htit += det[idet];
    hdetditnam[idet] = hname;
    hdet_dit[idet]=new TH1F(hname,htit,500,-10000,10000);


    hname = "h_";
    hname += det[idet];
    hname += "_reg_lin";
    htit = "Reg_";
    htit += det[idet];
    htit += "_lin_bcm1";
    hdetlinregnam[idet] = hname;
    hdetlin_reg[idet]=new TH2F(hname,htit,100,-4000,4000,100,-8000,8000);

    hname = "h_";
    hname += det[idet];
    hname += "_dit_lin";
    htit = "Dit_";
    htit += det[idet];
    htit += "_lin_bcm1";
    hdetlinditnam[idet] = hname;
    hdetlin_dit[idet]=new TH2F(hname,htit,100,-4000,4000,100,-8000,8000);

  }

//   for(UInt_t idet=0; idet<ndetcomb; idet++) {

//     TString hname;
//     hname = "h_asym_";
//     hname += detcomb[idet];
//     htit = "Asym_n_";
//     htit += detcomb[idet];
//     hdetcombnam[idet] = hname;
//     hdetcomb[idet]=new TH1F(hname,htit,750,-15000,15000);

//     hname = "h_asym_";
//     hname += detcomb[idet];
//     hname += "_reg";
//     htit = "Reg_Asym_n_";
//     htit += detcomb[idet];
//     hdetcombregnam[idet] = hname;
//     hdetcomb_reg[idet]=new TH1F(hname,htit,500,-10000,10000);

//     hname = "h_asym_";
//     hname += detcomb[idet];
//     hname += "_dit";
//     htit = "Dit_Asym_n_";
//     htit += detcomb[idet];
//     hdetcombditnam[idet] = hname;
//     hdetcomb_dit[idet]=new TH1F(hname,htit,500,-10000,10000);

//   }

  
  TString hname = "h_DD";
  htit = "DD12";
  hDDnam = hname;
  hDD = new TH1F(hname,htit,400,-2000,2000);
  hname="h_bcmlin";
  htit = "BCM_Lin";
  hbcmlinnam = hname;
  hbcmlin = new TH2F(hname,htit,100,-4000,4000,100,-1000,1000);

  //
  // loop over runs
  //     

  TString slopefile;
  for(UInt_t irun=0; irun<nruns; irun++) {
    
    runnumber = runlist[irun];
    
    slug = sluglist[irun];
    slowsign = dblist.GetSlowSign(runnumber);

    slopefile = 
      Form("../slopes/resultsslug/avgSlopes_hydrogen_dither_slug%d.res",slug);
    cout << "Getting dithering slopes..." << endl;

    TCut rawcut_base = "ev_num>1 && !cond_low_beam && !cond_beam_burp";
    rawcut_base += "!cond_mon_saturate";
    rawcut_base += "!cond_startup && !cond_pair_seq && !cond_evt_seq";
    rawcut_base += "!cond_adcx_bad && !cond_adcx_dac_burp";
    //    rawcut_base += "!cond_pos_burp";
    //    rawcut_base += "!cond_det_saturate && !cond_mon_saturate";
    //    rawcut_base += "!cond_c_burp";
    //    rawcut_base += "bmwcyc==0";

    TCut prcut_base = "ok_cut";
    //    prcut_base += "evt_bmwcyc[0]==0 && evt_bmwcyc[1]==0";
    
    Double_t minL = (Double_t)dblist.GetLeftDetLo(runnumber);
    Double_t maxL = (Double_t)dblist.GetLeftDetHi(runnumber);
    Double_t minR = (Double_t)dblist.GetRightDetLo(runnumber);
    Double_t maxR = (Double_t)dblist.GetRightDetHi(runnumber);

    TCut Lcut_raw;
    TCut Lcut_pr;
    if (maxL>minL && maxL>0) {
      Lcut_raw = Form("  ev_num>%d &&   ev_num<%d",minL,maxL);
      Lcut_pr  = Form("m_ev_num>%d && m_ev_num<%d",minL,maxL);
    } else if (maxL<0) {
      Lcut_raw = "0";
      Lcut_pr  = "0";
    } else {
      Lcut_raw = "1";
      Lcut_pr = "1";
    }

    TCut Rcut_raw;
    TCut Rcut_pr;
    if (maxR>minR && maxR>0) {
      Rcut_raw = Form("  ev_num>%d &&   ev_num<%d",minR,maxR);
      Rcut_pr  = Form("m_ev_num>%d && m_ev_num<%d",minR,maxR);
    } else if (maxR<0) {
      Rcut_raw = "0";
      Rcut_pr  = "0";
    } else {
      Rcut_raw = "1";
      Rcut_pr = "1";
    }

    TCut rawcut[4];
    TCut prcut[4];
    rawcut[0] = rawcut_base;
    prcut[0]  = prcut_base;
    rawcut[1] = rawcut_base + Lcut_raw;
    prcut[1]  = prcut_base  + Lcut_pr;
    rawcut[2] = rawcut_base + Rcut_raw;
    prcut[2]  = prcut_base  + Rcut_pr;
    rawcut[3] = rawcut_base + Rcut_raw + Lcut_raw;
    prcut[3]  = prcut_base  + Rcut_pr  + Lcut_pr;

    cout << "iteration     " << irun << endl;
    cout << "runnum    " << runlist[irun] << endl;

    Int_t fopened = open(runlist[irun],"regress"); 
    if(fopened==-1) continue;

    DitherAlias(asy,slopefile);    


////////////////////////////////////////////// 
   for(UInt_t ibcm=0; ibcm<nbcm; ibcm++) {
   TString drawcommand = "bcm"+bcm[ibcm];
   raw->Draw(drawcommand+">>htemp",rawcut[0],"goff");
    //    TH1* htemp = (TH1*) gROOT->FindObject("htemp");
    if (htemp && htemp->GetEntries()>0) {
      Int_t nev = htemp->GetEntries();
      rawbcm[ibcm]  = htemp->GetMean(); 
      rawbcmrms[ibcm] = htemp->GetRMS(); 
      rawbcmerr[ibcm] = rawbcmrms[ibcm] / sqrt(nev);
    } else {
      rawbcm[ibcm] = -1.E6; 
      rawbcmrms[ibcm] = -1.E6; 
      rawbcmerr[ibcm] =  -1.E6;
    }
    cerr << "***bcm: " << drawcommand << endl;
    cerr << "cut: " << rawcut[0] << endl;
    cerr << "Entries: " << htemp->GetEntries() << endl;
    cerr << "rawbcm: " << rawbcm[ibcm] << endl;
    cerr << "rawbcm_rms: " << rawbcmrms[ibcm] << endl;
    cerr << "rawbcm_err: " << rawbcmerr[ibcm] << endl;
   }
///////////////////////////////////////////////////////
 for(UInt_t ibpm=0; ibpm<nbpm; ibpm++) {
      TString drawcommand = "bpm"+bpm[ibpm];
      raw->Draw(drawcommand+">>htemp",rawcut[0],"goff");
      if (htemp && htemp->GetEntries()>0) {
      Int_t nev = htemp->GetEntries();
        bpmdat[ibpm] = htemp->GetMean();
       bpmrms[ibpm] = htemp->GetRMS();
       bpmerr[ibpm] = bpmdiffrms[ibpm] / sqrt(nev);
      } else {
        bpm[ibpm] = -1.E6;
       bpmrms[ibpm] = -1.E6;
       bpmerr[ibpm] =  -1.E6;
      }
    }

///////////////////////////////////////////////////////
   for(UInt_t idet=0; idet<ndet; idet++) {

      raw->Draw(det[idet]+ "/ bcm1>>htemp",rawcut[icdet[idet]],"goff");
      if (htemp && htemp->GetEntries()>0) {
        Int_t nev = htemp->GetEntries();
        normdet[idet] = htemp->GetMean();
        normdetrms[idet] = htemp->GetRMS();
        normdeterr[idet] = normdetrms[idet] / sqrt(nev);
      } else {
        normdet[idet] = -1.E6;
        normdetrms[idet] = -1.E6;
        normdeterr[idet] = -1.E6;
      }
      cerr << "***det: " << det[idet] << endl;
      cerr << "cut: " << rawcut[0] << endl;
      cerr << "Entries: " << htemp->GetEntries() << endl;
      cerr <<"normdet: "<< normdet[idet] << endl;
      cerr <<"normdet_rms: "<< normdetrms[idet] << endl;
      cerr <<"normdet_err: "<< normdeterr[idet] << endl;
     }
///////////////////////////////////////////////////////////////
  for(UInt_t idet=0; idet<ndet; idet++) {
      TString drawcommand = det[idet];
      raw->Draw(det[idet]+">>htemp",rawcut[icdet[idet]],"goff");
      if (htemp && htemp->GetEntries()>0) {
        Int_t nev = htemp->GetEntries();
        udet[idet] = htemp->GetMean();
        udetrms[idet] = htemp->GetRMS();
        udeterr[idet] = udetrms[idet] / sqrt(nev);
      } else {
        udet[idet] = -1.E6;
        udetrms[idet] = -1.E6;
        udeterr[idet] = -1.E6;
      }
     }

/////////////////////////////////////////////////////////////////
//   for(UInt_t iblumi=0; iblumi<nblumi; iblumi++) {
//     TString drawcommand = blumi[iblumi]+"n";
//     raw->Draw(blumi[iblumi]+"/ bcm"+bcm[2]+">>htemp",rawcut[0],"goff");
//     if (htemp && htemp->GetEntries()>0) {
//       Int_t nev = htemp->GetEntries();
//       normblumi[iblumi] = htemp->GetMean();
//       normblumirms[iblumi] = htemp->GetRMS();
//       normblumierr[iblumi] = normblumirms[iblumi] / sqrt(nev);
//     } else {
//       normblumi[iblumi] = -1.E6;
//       normblumirms[iblumi] = -1.E6;
//       normblumierr[iblumi] = -1.E6;
//     }
//   }
  
//   for(UInt_t iblumi=0; iblumi<nblumi; iblumi++) {
//     TString drawcommand = blumi[iblumi];
//     raw->Draw(blumi[iblumi]+">>htemp",rawcut[0],"goff");
//     if (htemp && htemp->GetEntries()>0) {
//       Int_t nev = htemp->GetEntries();
//       ublumi[iblumi] = htemp->GetMean();
//       ublumirms[iblumi] = htemp->GetRMS();
//       ublumierr[iblumi] = ublumirms[iblumi] / sqrt(nev);
//     } else {
//       ublumi[iblumi] = -1.E6;
//       ublumirms[iblumi] = -1.E6;
//       ublumierr[iblumi] = -1.E6;
//     }
//   }
////////////////////////////////////////////////////////
// for(UInt_t iflumi=0; iflumi<nflumi; iflumi++) {
//       TString drawcommand = flumi[iflumi]+"n";
//       raw->Draw(flumi[iflumi]+"/ bcm"+bcm[2]+">>htemp",rawcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//         Int_t nev = htemp->GetEntries();
//         normflumi[iflumi] = htemp->GetMean();
//         normflumirms[iflumi] = htemp->GetRMS();
//         normflumierr[iflumi] = normflumirms[iflumi] / sqrt(nev);
//       } else {
//         normflumi[iflumi] = -1.E6;
//         normflumirms[iflumi] = -1.E6;
//         normflumierr[iflumi] = -1.E6;
//       }
//      }
////////////////////////////////////////////////////
// for(UInt_t iflumi=0; iflumi<nflumi; iflumi++) {
//       TString drawcommand = flumi[iflumi];
//       raw->Draw(flumi[iflumi]+">>htemp",rawcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//         Int_t nev = htemp->GetEntries();
//         uflumi[iflumi] = htemp->GetMean();
//         uflumirms[iflumi] = htemp->GetRMS();
//         uflumierr[iflumi] = uflumirms[iflumi] / sqrt(nev);
//       } else {
//         uflumi[iflumi] = -1.E6;
//         uflumirms[iflumi] = -1.E6;
//         uflumierr[iflumi] = -1.E6;
//       }
//      }
/////////////////////////////////////////////////////
    
///////////////////////////////////////////////////////////////      
 
     // uncorrected detector asymmetries 
     for(UInt_t idet=0; idet<ndet; idet++) {

       // Setup cuts
        TString hrs;
        if(det[idet].Contains("1") || det[idet].Contains("2") ) hrs="L";
        else hrs="R";
	asy->Draw("asym_n_"+det[idet]+">>htemp",prcut[icdet[idet]],"goff");
	numevents=htemp->GetEntries();
	if (htemp && htemp->GetEntries()>0) {
	 Int_t nev = htemp->GetEntries();
	 rawdet[idet] = htemp->GetMean();
	 rawdetrms[idet] = htemp->GetRMS();
	 rawdeterr[idet] = rawdetrms[idet] / sqrt(nev);
  	 if (hrs.Contains("L")) {
  	   ok_stats[0] = 1; }
  	 else ok_stats[1] = 1;
	} else {
	  rawdet[idet] = -1.E6;
	  rawdetrms[idet] = -1.E6;
	  rawdeterr[idet] = -1.E6;
	  if (hrs.Contains("L")) {
	    ok_stats[0] = 0; }
	  else ok_stats[1] = 0;
	}
	TString prestr;
	if (slowsign == 1) prestr = "1*asym_n_";
	else prestr = "-1*asym_n_";
	fRootFile->cd();
     	asy->Draw(prestr+det[idet]+">>+"+hdetnam[idet],prcut[icdet[idet]],"goff");
	f_pan_standard->cd();
     }

     // uncorrected detector asymmetries 
//      for(UInt_t idet=0; idet<ndetcomb; idet++) {

//        // Setup cuts
//         TString hrs;
// 	asy->Draw("asym_n_"+detcomb[idet]+">>htemp",prcut[icdet[idet]],"goff");
// 	if (htemp && htemp->GetEntries()>0) {
// 	 Int_t nev = htemp->GetEntries();
// 	 rawdetcomb[idet] = htemp->GetMean();
// 	 rawdetcombrms[idet] = htemp->GetRMS();
// 	 rawdetcomberr[idet] = rawdetcombrms[idet] / sqrt(nev);
//        } else {
// 	 rawdetcomb[idet] = -1.E6;
// 	 rawdetcombrms[idet] = -1.E6;
// 	 rawdetcomberr[idet] = -1.E6;
//        }
// 	TString prestr;
// 	if (slowsign == 1) prestr = "1*asym_n_";
// 	else prestr = "-1*asym_n_";
// 	fRootFile->cd();
//      	asy->Draw(prestr+detcomb[idet]+">>+"+hdetcombnam[idet],prcut[icdet[idet]],"goff");
// 	f_pan_standard->cd();
//      }

  
////////////////////////////////////////////////////////////////    
    for(UInt_t ibpm=0; ibpm<nbpm; ibpm++) {
      TString drawcommand = "diff_bpm"+bpm[ibpm];
      asy->Draw(drawcommand+">>htemp",prcut[0],"goff");
      if (htemp && htemp->GetEntries()>0) {
      Int_t nev = htemp->GetEntries();
	bpmdiffdat[ibpm] = htemp->GetMean();
       bpmdiffrms[ibpm] = htemp->GetRMS();
       bpmdifferr[ibpm] = bpmdiffrms[ibpm] / sqrt(nev);
      } else {
	bpmdiffdat[ibpm] = -1.E6; 
       bpmdiffrms[ibpm] = -1.E6; 
       bpmdifferr[ibpm] =  -1.E6;
      }
    }   

// /////////////////////////////////////////////////////////////////
// for(UInt_t iinj=0; iinj<ninj; iinj++) {
//       TString drawcommand = injbpm[iinj];
//       raw->Draw(drawcommand+">>htemp","","goff");
//       if (htemp && htemp->GetEntries()>0) {
//       Int_t nev = htemp->GetEntries();
//        injbpmdat[iinj] = htemp->GetMean();
//        injbpmrms[iinj] = htemp->GetRMS();
//        injbpmerr[iinj] = injbpmrms[iinj] / sqrt(nev);
//       } else {
//        injbpmdat[iinj] = -1.E6;
//        injbpmrms[iinj] = -1.E6;
//        injbpmerr[iinj] =  -1.E6;
//       }
//     } 
//
// ///////////////////////////////////////////////////////////////// 
// for(UInt_t iinj=0; iinj<ninj; iinj++) {
//       TString drawcommand = "diff_"+injbpm[iinj];
//       asy->Draw(drawcommand+">>htemp",prcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//       Int_t nev = htemp->GetEntries();
//        injbpmdiffdat[iinj] = htemp->GetMean();
//        injbpmdiffrms[iinj] = htemp->GetRMS();
//        injbpmdifferr[iinj] = injbpmdiffrms[iinj] / sqrt(nev);
//       } else {
//        injbpmdiffdat[iinj] = -1.E6;
//        injbpmdiffrms[iinj] = -1.E6;
//        injbpmdifferr[iinj] =  -1.E6;
//       }
//     }
//
////////////////////////////////////////////////////////////////    
    for(UInt_t ibcm=0; ibcm<nbcm; ibcm++) {
      TString drawcommand = "asym_bcm"+bcm[ibcm];
      asy->Draw(drawcommand+">>htemp",prcut[0],"goff");
      if (htemp && htemp->GetEntries()>0) {
      Int_t nev = htemp->GetEntries(); 
	bcmasymdat[ibcm] = htemp->GetMean();
     bcmasymrms[ibcm] = htemp->GetRMS();
     bcmasymerr[ibcm] = bcmasymrms[ibcm] / sqrt(nev);
      } else {
	bcmasymdat[ibcm] = -1.E6;
        bcmasymrms[ibcm] = -1.E6;
        bcmasymerr[ibcm] =  -1.E6;
      }
    }
  
//////////////////////////////////////////////////////////////////    
  
    TString drawcommand = "asym_bcm1-asym_bcm2"; 
    asy->Draw(drawcommand+">>htemp",prcut[0],"goff");
    if (htemp && htemp->GetEntries()>0) {
     Int_t nev = htemp->GetEntries();
      DD12 = htemp->GetMean();
      DD12rms =htemp->GetRMS();
      DD12err = DD12rms / sqrt(nev);
    } else {
      DD12 = -1.E6;
      DD12rms = -1.E6;
      DD12err =  -1.E6;
    }
    fRootFile->cd();
    asy->Draw(drawcommand+">>+"+hDDnam,prcut[0],"goff");
    asy->Draw(drawcommand+":asym_bcm1>>+"+hbcmlinnam,prcut[0],"goff");
    f_pan_standard->cd();
    
//////////////////////////////////////////////////////////////// 
// regressed detector asymmetries

    for(UInt_t idet=0; idet<ndet; idet++) {

      asy->Draw("reg_asym_n_"+det[idet]+">>htemp",prcut[icdet[idet]],"goff");
      if (htemp && htemp->GetEntries()>0) {
	Int_t nev = htemp->GetEntries();
	regdet[idet] = htemp->GetMean();
	regdetrms[idet] = htemp->GetRMS();
	regdeterr[idet] = regdetrms[idet] / sqrt(nev);
      } else {
	regdet[idet] = -1.E6;
	regdetrms[idet] = -1.E6;
	regdeterr[idet] = -1.E6;
      }
      TString prestr;
      if (slowsign == 1) prestr = "1*reg_asym_n_";
      else prestr = "-1*reg_asym_n_";
      fRootFile->cd();
      asy->Draw(prestr+det[idet]+">>+"+hdetregnam[idet],prcut[icdet[idet]],"goff");
      asy->Draw(prestr+det[idet]+":asym_bcm1>>+"+hdetlinregnam[idet],prcut[icdet[idet]],"goff");
      f_pan_standard->cd();
    }  

// regressed detector combination asymmetries

//     for(UInt_t idet=0; idet<ndetcomb; idet++) {

//       asy->Draw("reg_asym_n_"+detcomb[idet]+">>htemp",prcut[icdet[idet]],"goff");
//       if (htemp && htemp->GetEntries()>0) {
// 	Int_t nev = htemp->GetEntries();
// 	regdetcomb[idet] = htemp->GetMean();
// 	regdetcombrms[idet] = htemp->GetRMS();
// 	regdetcomberr[idet] = regdetcombrms[idet] / sqrt(nev);
//       } else {
// 	regdetcomb[idet] = -1.E6;
// 	regdetcombrms[idet] = -1.E6;
// 	regdetcomberr[idet] = -1.E6;
//       }
//       TString prestr;
//       if (slowsign == 1) prestr = "1*reg_asym_n_";
//       else prestr = "-1*reg_asym_n_";
//       fRootFile->cd();
//       asy->Draw(prestr+detcomb[idet]+">>+"+hdetcombregnam[idet],prcut[icdet[idet]],"goff");
//       f_pan_standard->cd();
//     }  
    
 
// // regressed blumi asymmetries

//     for(UInt_t iblumi=0; iblumi<nblumi; iblumi++) {
//      asy->Draw("reg_asym_n_"+blumi[iblumi]+">>htemp",prcut[0],"goff");
//      if (htemp && htemp->GetEntries()>0) {
//        Int_t nev = htemp->GetEntries();
//        regblumi[iblumi] = htemp->GetMean();
//        regblumirms[iblumi] = htemp->GetRMS();
//       regblumierr[iblumi] = regblumirms[iblumi] / sqrt(nev);
//       } else {
//         regblumi[iblumi] = -1.E6;
//         regblumirms[iblumi] = -1.E6;
//         regblumierr[iblumi] = -1.E6;
//       }
//    }
////////////////////////////////////////////////////////
// // regressed flumi asymmetries
//
//     for(UInt_t iflumi=0; iflumi<nflumi; iflumi++) {
//       asy->Draw("reg_asym_n_"+flumi[iflumi]+">>htemp",prcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//         Int_t nev = htemp->GetEntries();
//         regflumi[iflumi] = htemp->GetMean();
//         regflumirms[iflumi] = htemp->GetRMS();
// 	regflumierr[iflumi] = regflumirms[iflumi] / sqrt(nev);
//       } else {
// 	regflumi[iflumi] = -1.E6;
// 	regflumirms[iflumi] = -1.E6;
// 	regflumierr[iflumi] = -1.E6;
//       }
//     }
/////////////////////////////////////////////////////////////////////
    // dithering detector asymmetries
    
    for(UInt_t idet=0; idet<ndet; idet++) {
      asy->Draw("dit_asym_n_"+det[idet]+">>htemp",prcut[icdet[idet]],"goff");
      if (htemp && htemp->GetEntries()>0) {
        Int_t nev = htemp->GetEntries();
        ditdet[idet] = htemp->GetMean();
        ditdetrms[idet] = htemp->GetRMS();
	ditdeterr[idet] = ditdetrms[idet] / sqrt(nev);
      } else {
	ditdet[idet] = -1.E6;
	ditdetrms[idet] = -1.E6;
	ditdeterr[idet] = -1.E6;
      }
      TString prestr;
      if (slowsign == 1) prestr = "1*(dit_asym_n_";
      else prestr = "-1*(dit_asym_n_";
      fRootFile->cd();
      asy->Draw(prestr+det[idet]+")>>+"+hdetditnam[idet],prcut[icdet[idet]],"goff");
      asy->Draw(prestr+det[idet]+"):asym_bcm1>>+"+hdetlinditnam[idet],prcut[icdet[idet]],"goff");
      f_pan_standard->cd();
    }
    

// ///////////////////////////////////////////////////////////////
// // dethering blumi asymmetries
//     // dithering detector asymmetries
    
//     for(UInt_t iblumi=0; iblumi<nblumi; iblumi++) {
//       asy->Draw("dit_asym_n_"+blumi[iblumi]+>>htemp",prcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//         Int_t nev = htemp->GetEntries();
//         ditblumi[iblumi] = htemp->GetMean();
//         ditblumirms[iblumi] = htemp->GetRMS();
// 	ditblumierr[iblumi] = ditblumirms[iblumi] / sqrt(nev);
//       } else {
// 	ditblumi[iblumi] = -1.E6;
// 	ditblumirms[iblumi] = -1.E6;
// 	ditblumierr[iblumi] = -1.E6;
//       }
// //       TString prestr;
// //       if (slowsign == 1) prestr = "1*(";
// //       else prestr = "-1*(";
// //       fRootFile->cd();
// //       asy->Draw(prestr+drawstr+")>>+"+hblumiditnam[iblumi],prcut[0],"goff");
// //       asy->Draw(prestr+drawstr+"):asym_bcm1>>+"+hblumilinditnam[iblumi],
// // 		prcut[0],"goff");
// //       f_pan_standard->cd();
//     }

//  //////////////////////////////////////////////////////////////////
//  // dethering flumi asymmetries
//      for(UInt_t iflumi=0; iflumi<nflumi; iflumi++) {
//       asy->Draw("dit_asym_n_"+flumi[iflumi]+">>htemp",prcut[0],"goff");
//       if (htemp && htemp->GetEntries()>0) {
//         Int_t nev = htemp->GetEntries();
//         ditflumi[iflumi] = htemp->GetMean();
//         ditflumirms[iflumi] = htemp->GetRMS();
//         ditflumierr[iflumi] = ditflumirms[iflumi] / sqrt(nev);
//        } else {
//          ditflumi[iflumi] = -1.E6;
//          ditflumirms[iflumi] = -1.E6;
//          ditflumierr[iflumi] = -1.E6;
//        }
//      }
///////////////////////////////////////////////////////////////////
//
    //
    // Epics words
    //  
    for(UInt_t iepi=0; iepi<nEpics; iepi++) {
      //      cout << strEpics[iepi] << endl;
      raw->Draw(strEpics[iepi]+">>htemp",rawcut[0],"goff");
      if (htemp && htemp->GetEntries()>0) {
        Int_t nev = htemp->GetEntries();
        EpicsData[iepi] = htemp->GetMean();
      } else {
	EpicsData[iepi] = -1.E6;
      }
      f_pan_standard->cd();
    }
    

    //    
    // FILL TREE
    //
    stats->Fill();

    if (f_pan_standard) f_pan_standard->Close();
    if (f_pan_regress) f_pan_regress->Close();
    if (f_pan_dither) f_pan_dither->Close();
 
  }
 
  // Write tree, root file and complete macro

  fRootFile->cd();
  stats->Write("",TObject::kOverwrite);

  for(UInt_t idet=0; idet<ndet; idet++) {
    hdet[idet]->Write("",TObject::kOverwrite);
    hdet_reg[idet]->Write("",TObject::kOverwrite);
    hdet_dit[idet]->Write("",TObject::kOverwrite);
    hdetlin_reg[idet]->Write("",TObject::kOverwrite);
    hdetlin_dit[idet]->Write("",TObject::kOverwrite);
  }
  hDD->Write("",TObject::kOverwrite);
  hbcmlin->Write("",TObject::kOverwrite);

  fRootFile->Write("",TObject::kOverwrite);

  cout<<"done"<< endl;
  gROOT->ProcessLine(".qqqqqqqqqq");  
  return 1;

}