//////////////////////////////////////////////////////////
//   This file has been automatically generated 
//     (Tue Mar 20 13:09:57 2018 by ROOT version5.34/10)
//   from TTree save/save
//   found on file: datatuple.root
//
//   DP data analysis
//
//////////////////////////////////////////////////////////
int ana_DP(const char* inFile = "datatuple.root",
	   const char* outFile = "ana_DP.root" ){


//Reset ROOT and connect tree file
   gROOT->Reset();

   //--- input file 
   TFile *inF = TFile::Open(inFile, "read");
   TTree * save = (TTree *) inF->Get("save");


  //--- analysis output file and analysis plots ---
  TFile *outF = new TFile(outFile,"recreate");

  
//Declaration of leaves types
   Int_t           runID;
   Int_t           spillID;
   Int_t           eventID;
   Int_t           triggerBits;
   Int_t           nTracks;
   Int_t           nHits1;
   Int_t           nHits2;
   Double_t        chisq1;
   Double_t        chisq2;
   Double_t        tx1_st1;
   Double_t        tx1;
   Double_t        ty1;
   Double_t        x1_st1;
   Double_t        x1;
   Double_t        y1;
   Double_t        tx2_st1;
   Double_t        tx2;
   Double_t        ty2;
   Double_t        x2_st1;
   Double_t        x2;
   Double_t        y2;
   Double_t        pz1;
   Double_t        pz2;

   // Set branch addresses.
   save->SetBranchAddress("runID",&runID);
   save->SetBranchAddress("spillID",&spillID);
   save->SetBranchAddress("eventID",&eventID);
   save->SetBranchAddress("triggerBits",&triggerBits);
   save->SetBranchAddress("nTracks",&nTracks);
   save->SetBranchAddress("nHits1",&nHits1);
   save->SetBranchAddress("nHits2",&nHits2);
   save->SetBranchAddress("chisq1",&chisq1);
   save->SetBranchAddress("chisq2",&chisq2);
   save->SetBranchAddress("tx1_st1",&tx1_st1);
   save->SetBranchAddress("tx1",&tx1);
   save->SetBranchAddress("ty1",&ty1);
   save->SetBranchAddress("x1_st1",&x1_st1);
   save->SetBranchAddress("x1",&x1);
   save->SetBranchAddress("y1",&y1);
   save->SetBranchAddress("tx2_st1",&tx2_st1);
   save->SetBranchAddress("tx2",&tx2);
   save->SetBranchAddress("ty2",&ty2);
   save->SetBranchAddress("x2_st1",&x2_st1);
   save->SetBranchAddress("x2",&x2);
   save->SetBranchAddress("y2",&y2);
   save->SetBranchAddress("pz1",&pz1);
   save->SetBranchAddress("pz2",&pz2);


   //--- canvas  & Histograms 
   
   TCanvas * c1 = new TCanvas("c1", "DP: Single Tracks", 10, 20, 600, 800);
   TCanvas * c2 = new TCanvas("c2", "DP: Dimuons", 500, 100, 600, 800);

   TH1F * h1Trig = new TH1F("h1Trig", "Trigger Bits", 130, 0, 130); 

   TH1F * h1Zvtx1  = new TH1F("h1Zvtx1", "Projected Z-vtx mu+", 100, -1000, 1000); // cm 
   TH1F * h1Zvtx2  = new TH1F("h1Zvtx2", "Projected Z-vtx mu-", 100, -1000, 1000); // cm 
   TH1F * h1ZvtxdZ  = new TH1F("h1ZvtxdZ", "Projected Z-vtx(mu+)-(mu-)", 100, -1000, 1000); // cm 

   
     // Variables
   const Float_t m_mu = 0.1057; // muon_mass = 105.6583745(24) MeV/c2
   const Float_t y_offset =  1.6; // beam y-off-set = 1.6 cm
   const Float_t PI   =  3.1415926;
   const Float_t PIm  = -3.1415926;
   const Float_t PI2 = 6.2831853;
  
   Double_t vZ1=-999, vZ2=-999, dZ=-999;
   Double_t mass=-999, mass_avgZ = -999; // raw dimuon mass, raw dimuon mass with recalcualted mean-Vz for opening angle
   

   
//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// save->SetBranchStatus("*",0);  // disable all branches
// TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname

   Long64_t nentries = save->GetEntries();
   cout << "nentries =" << nentries << endl;
   
   Long64_t nbytes = 0;
  for (Long64_t i=0; i<nentries;i++) {
      nbytes += save->GetEntry(i);

      //
      //quality cuts
      //

      if (nHits1 < 10) continue;
      if (nHits2 < 10) continue;
      if (chisq1 > 20) continue;
      if (chisq2 > 20) continue;

      if (fabs(x1_st1) >100 || fabs(x2_st1) >100) continue;
      if (fabs(tx1_st1) >100 || fabs(tx2_st1) >100) continue;

      if (fabs(y1) >100 || fabs(y2) >100) continue;
      if (fabs(ty1) >100 || fabs(ty2) >100) continue;

      vZ1 = (y_offset - y1)/ty1;
      vZ2 = (y_offset - y2)/ty2;
      dZ = vZ1 - vZ2;

      if (fabs(vZ1) >1000 || fabs(vZ2) >1000) continue;




      //loop events, fill histogram ---
      h1Trig->Fill(triggerBits);
      h1Zvtx1->Fill(vZ1);
      h1Zvtx2->Fill(vZ2);
      h1ZvtxdZ->Fill(dZ);
      
	
  } // end of loop 

  //--make plots

  c1->cd()->SetLogy();
  h1Trig->Draw();
  
  c2->cd();
  c2->Divide(2,2);

  c2->cd(1);
  h1Zvtx1->Draw();

  c2->cd(2);
  h1Zvtx2->Draw();

  c2->cd(3);
  h1ZvtxdZ->Draw();

  c2->cd(4)->SetLogy();
  h1Trig->Draw();

  //--- save output file
  outF->cd();

  h1Trig->Write();
  h1Zvtx1->Write();
  h1Zvtx2->Write();
  h1ZvtxdZ->Write();

  c1->Write();
  c2->Write();
  
  outFile->Close();
  
  return 0;

}