#define Muon_cxx
#include "Muon.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TRefArray.h"
#include "TString.h"
#include <iostream>
//#include <string>
using namespace std;
#include <fstream>
ofstream myout;



void   Muon::Loop(char* name, int museg, float mupt, int icol1)
{
//   In a ROOT session, you can do:
//      Root > .L T.C
//      Root > T trg
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch
   if (fChain == 0) return;
//   histogram booking

   TH1F* hsmuons = new TH1F("hsmuons","no. single muons",10,0.,10.);
   TH1F* hdmuons = new TH1F("hdmuons","no. di-muons",10,0.,10.);
   TH1F *h1 = new TH1F("h1","Di-muon mass ",40, 2, 4);
   TH1F *h2 = new TH1F("h2","Di-muon mass ",100, 1.5, 10);  
   TH1F *h3 = new TH1F("h3","Di-muon mass with trigger",40, 2, 4);  
   TH1F *h4 = new TH1F("h4","Di-muon mass with di trigger",40, 2, 4);
   TH1F *h5 = new TH1F("h5","Di-muon mass with di trigger",100, 1.5, 10);
   TH1F *h6 = new TH1F("h6","Di-muon mass with di trigger and mass cut",40, 2, 4);
   TH1F *h7 = new TH1F("h7","Di-muon mass with di trigger and central rank=1",40, 2, 4);
   TH1F *h8 = new TH1F("h8","central match for Psi mass ",7, -2, 5);
   TH1F *h9 = new TH1F("h9","Number of muons ",17, -2, 15);
   TH1F *h10 = new TH1F("h10","Di-trk mass ",40, 2, 4);
   TH1F *h20 = new TH1F("h20","Di-trk mass ",100, 1.5, 10);
   TH1F *h30 = new TH1F("h30","events with two rank = 1 muons ",4, -2, 2);
   TH1F *h40 = new TH1F("h40","Di-muon mass with at least 2 central rank=1",40, 2, 4);
   TH1F *h100 = new TH1F("h100","Di-muon mass with specific trigger",40, 2, 4);   
   TH1F *h101 = new TH1F("h101","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h102 = new TH1F("h102","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h103 = new TH1F("h103","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h104 = new TH1F("h104","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h105 = new TH1F("h105","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h106 = new TH1F("h106","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h107 = new TH1F("h107","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h108 = new TH1F("h108","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h109 = new TH1F("h109","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h110 = new TH1F("h110","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h111 = new TH1F("h111","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h112 = new TH1F("h112","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h113 = new TH1F("h113","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h114 = new TH1F("h114","Di-muon mass with specific trigger",40, 2, 4);
   TH1F *h115 = new TH1F("h115","Di-muon mass with specific trigger",40, 2, 4);

   Int_t nentries = Int_t(fChain->GetEntriesFast());
   Int_t ievent=0;
   Int_t itrgdi=0;
   Int_t ialltrigdi=0;
   int aa = 1000; int nn = 0;
   myout.open("run_event_number.txt");

   Int_t nbytes = 0, nb = 0;
   for (Int_t jentry=0; jentry<nentries;jentry++) {
     Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
     if (ientry < 0) break;
     nb = fChain->GetEntry(jentry);   nbytes += nb;
     // if (Cut(ientry) < 0) continue;
     
     ievent++;
     TMBMuon* muon = 0;
     TMBTrks* trks = 0;
     TMBTrig *  trig = 0;
     TMBGlob  * glob = 0;
     
      //
     const float MuonMass=0.10565;
     int smu=0;
     int dmu=0;
     //int hits, _wa, _wbc, _sa, _sbc;
     int q[50];
     int tq[500];
     int ntrks=0;
     float px[50],py[50],pz[50],pt[50],trketa[50];
     int rank[50];     
     float M,E1,E2;
     float tM,tE1,tE2;
     float tpx[500],tpy[500],tpz[500],tpt[500];
     int trigger_id[20];
     int rank_events;

     for(Int_t ii=0;ii<20;ii++) {trigger_id[ii]=0;}
     

     //Event information
     glob = fGlob;
     int eventnum=glob->evtno();
     int runnum= glob->runno();
     
     h9->Fill(fMuon->GetLast()+1);

     rank_events=0;

     for(Int_t i=0; i<fMuon->GetLast()+1; i++) {
       
       muon = (TMBMuon* ) fMuon->At(i);
       // hits = muon->nhit();
       // _wa  = (hits % 10);
       //_wbc = (hits % 1000)/10;
       //_sa  = (hits % 10000)/1000;
       //_sbc = (hits/10000);         
       
       //  if(_wa  > 1 && _wbc > 2 && _sa  > 0 && _sbc > 0) {

       if(muon->nseg()>0) {
	 
	 px[dmu]=muon->px();
         py[dmu]=muon->py();
	 pz[dmu]=muon->pz();        
         q[dmu]=muon->charge();
         pt[dmu]=muon->pT();
         rank[dmu]=muon->centralrank();
	 if(muon->centralrank()==1) rank_events++;

//           Int_t track_index=muon->TrkIndex();
//           cout << "track index " << track_index << endl;
//           cout << "number of tracks " << fTrks->GetLast() << endl;
//           trks = (TMBTrks* ) fTrks->At(track_index);
//           trketa[dmu]=trks->eta();
//           cout << "track eta " << trketa[dmu] << endl;
         dmu++; 
         if(muon->nseg() ==3 && muon->pT() > 2.5 && abs(muon->eta())<2.4) smu++;
       }
     }

     Int_t rank_value=-1;
     if(rank_events>1) {rank_value=1;}
     h30->Fill(rank_value);
    
     Bool_t trig_di = false;
     

     for(Int_t i=0; i<fTrig->GetLast()+1; i++) {
       trig =(TMBTrig*) fTrig->At(i);
       
       if (trig->getTrgName() == TString("2MU_A_L2M0")) {
	 trig_di = true;
	 trigger_id[0]=1;
       }
       if (trig->getTrgName() == TString("2MU_W_L2M0_2TRK2")) {
	 trig_di = true;
	 trigger_id[1]=1;
       }
       if (trig->getTrgName() == TString("2MU_A_L2ETAPHI")) {
	 trig_di = true;
	 trigger_id[2]=1;
       }
       if (trig->getTrgName() == TString("2MU_C_2L2_TRKLO")) {
	 trig_di = true;
	 trigger_id[3]=1;
       }
       if (trig->getTrgName() == TString("mu2ptxbtxx_fz")) {
	 trig_di = true;
	 trigger_id[4]=1;
       }
       if (trig->getTrgName() == TString("mu2ptxctxx_fz")) {
	 trig_di = true;
	 trigger_id[5]=1;
       }
       if (trig->getTrgName() == TString("mu2ptxatxx")) {
	 trig_di = true;
	 trigger_id[6]=1;
       }
       if (trig->getTrgName() == TString("2MU_A_2L2")) {
	 trig_di = true;
	 trigger_id[7]=1;
       }
       if (trig->getTrgName() == TString("mu2ptxwtxx_fz")) {
	 trig_di = true;
	 trigger_id[8]=1;
       }
       if (trig->getTrgName() == TString("mu2ptxatxx_fz")) {
	 trig_di = true;
	 trigger_id[9]=1;
       }
       if (trig->getTrgName() == TString("2MU_CA_2L2L_2TRK")) {
	 trig_di = true;
	 trigger_id[10]=1;
       }
       if (trig->getTrgName() == TString("2MU_CA_2L2_2TRK")) {
	 trig_di = true;
	 trigger_id[11]=1;
       }
       if (trig->getTrgName() == TString("2MU_CFT_2L2_2TRK")) {
	 trig_di = true;
	 trigger_id[12]=1;
       }
       if (trig->getTrgName() == TString("2MU_CA_2L2L")) {
	 trig_di = true;
	 trigger_id[13]=1;
       }
       
     }
     
//       if(!trig_di)  
//         {
//  	 for(Int_t i=0; i<fTrig->GetLast()+1; i++) {
//  	   trig =(TMBTrig*) fTrig->At(i);
//  	   cout << trig->getTrgName() << endl;
//  	 }
//  	 cout << " " << endl;
//         }
     

     if(trig_di) ialltrigdi++;

//      if( (dmu==2) && (q[0]*q[1]==-1) && pt[0]>1.5 && pt[1]>1.5 ) {
//         E1 = sqrt( pt[0]*pt[0] + pz[0]*pz[0] + MuonMass*MuonMass );
//         E2 = sqrt( pt[1]*pt[1] + pz[1]*pz[1] + MuonMass*MuonMass );
//         M = (E1+E2) * (E1+E2);
//         M -= (px[0]+px[1]) * (px[0]+px[1]);
//         M -= (py[0]+py[1]) * (py[0]+py[1]);
//         M -= (pz[0]+pz[1]) * (pz[0]+pz[1]);
//         if(M<0) {
//  	 cout << "Muon Mass negative" << endl;
//         }
//         M = sqrt(M);

     
     Bool_t write_event=false;
     for(Int_t ii=0;ii<dmu;ii++) {
       for(Int_t jj=0;jj<dmu;jj++) {
	 //	 if((q[ii]*q[jj]==-1) && pt[ii]>1.5 && pt[jj]>1.5 && q[ii]<0  && rank[ii]==1 && rank[jj]==1) {
           if((q[ii]*q[jj]==-1) && pt[ii]>1.5 && pt[jj]>1.5 && q[ii]<0 && rank_value==1) {
	   E1 = sqrt( pt[ii]*pt[ii] + pz[ii]*pz[ii] + MuonMass*MuonMass );
	   E2 = sqrt( pt[jj]*pt[jj] + pz[jj]*pz[jj] + MuonMass*MuonMass );
	   M = (E1+E2) * (E1+E2);
	   M -= (px[ii]+px[jj]) * (px[ii]+px[jj]);
	   M -= (py[ii]+py[jj]) * (py[ii]+py[jj]);
	   M -= (pz[ii]+pz[jj]) * (pz[ii]+pz[jj]);
	   if(M<0) {
	     cout << "Muon Mass negative" << endl;
	   }
	   if(M>0) {M = sqrt(M);}
	   
	   h1->Fill(M);
	   h2->Fill(M);
           if( (rank[ii]==1) && (rank[jj]==1)) {h7->Fill(M);}
           if(rank_value==1) {h40->Fill(M);}

           if( (M>3.0) && (M<3.2)) {
              h8->Fill(rank[ii]);
              h8->Fill(rank[jj]);
	   } 
	   if( (trig_di) && (M>2.0) && (M<4.0)) {h3->Fill(M);}
	   if(trig_di) {h4->Fill(M);}
	   if(trig_di) {h5->Fill(M);}
	   if( (trig_di) && (M>2.9) && (M<3.3)) {h6->Fill(M);}
	   if( (trig_di) && (M>2.0) && (M<4.0)) 
	     {
	       write_event=true;
	     }
	 }
       }
     }
      
     if(write_event) {myout << runnum << " " << eventnum << endl;}

       if(trig_di) {
	 
	 if(trigger_id[0]==1){ h100->Fill(M);}
         if(trigger_id[1]==1){ h101->Fill(M);}
	 if(trigger_id[2]==1){ h102->Fill(M);}
	 if(trigger_id[3]==1){ h103->Fill(M);}
	 if(trigger_id[4]==1){ h104->Fill(M);}
	 if(trigger_id[5]==1){ h105->Fill(M);}
	 if(trigger_id[6]==1){ h106->Fill(M);}
	 if(trigger_id[7]==1){ h107->Fill(M);}
	 if(trigger_id[8]==1){ h108->Fill(M);}
	 if(trigger_id[9]==1){ h109->Fill(M);}
	 if(trigger_id[10]==1){ h110->Fill(M);}
	 if(trigger_id[11]==1){ h111->Fill(M);}
	 if(trigger_id[12]==1){ h112->Fill(M);}
         if(trigger_id[13]==1){ h113->Fill(M);}
	 if(trigger_id[14]==1){ h114->Fill(M);} 
         if(trigger_id[15]==1){ h115->Fill(M);}

       }

       if(trig_di) itrgdi++;
       
//         for(Int_t i=0; i<fTrks->GetLast()+1; i++) {
	 
//  	 if(i<500 ) {
//  	   trks = (TMBTrks* ) fTrks->At(i);
//  	   tpx[ntrks]=trks->px();
//  	   tpy[ntrks]=trks->py();
//  	   tpz[ntrks]=trks->pz();        
//  	   tq[ntrks]=trks->charge();
//  	   tpt[ntrks]=trks->pT();
//  	   ntrks++;
//         }
       
//         }
     
//       if( (trig_di) && (dmu>1))  {
      
//         for(Int_t ii=0;ii<ntrks;ii++) {
//  	 for(Int_t jj=0;jj<ntrks;jj++) {
//  	   int good_pz=-1;
//  	   if(tpz[ii]<100 && tpz[ii]>-100 && tpz[jj]<100 && tpz[jj]>-100 ) {good_pz=1;} 
//  	   if((tq[ii]*tq[jj]==-1) && tpt[ii]>2.0 && tpt[jj]>2.0 && tq[ii]<0 && good_pz==1) {
//  	     tE1 = sqrt( tpt[ii]*tpt[ii] + tpz[ii]*tpz[ii] + MuonMass*MuonMass );
//  	     tE2 = sqrt( tpt[jj]*tpt[jj] + tpz[jj]*tpz[jj] + MuonMass*MuonMass );
//  	     tM = (tE1+tE2) * (tE1+tE2);
//  	     tM -= (tpx[ii]+tpx[jj]) * (tpx[ii]+tpx[jj]);
//  	     tM -= (tpy[ii]+tpy[jj]) * (tpy[ii]+tpy[jj]);
//  	     tM -= (tpz[ii]+tpz[jj]) * (tpz[ii]+tpz[jj]);
//  	     if(tM<0) {
//  	       cout << "Track Mass negative" << endl;
//  	     }
//  	     if(tM>0) {tM = sqrt(tM);}
//  	     h10->Fill(tM);
//  	     h20->Fill(tM);
//  	   }
//  	 }
//         }
       
//       }

      
     hsmuons->Fill((float) smu);
     hdmuons->Fill((float) dmu);
     
     if(ievent/aa>nn) cout << "Event #" << ievent << endl; nn = ievent/aa;
  
   } // jentry

TFile file(name,"recreate"); 
 hsmuons->Write(); hdmuons->Write(); h1->Write(); h2->Write(); h3->Write(); h4->Write(); h5->Write(); h6->Write(); h7->Write(); h8->Write(); h9->Write(); h10->Write(); h20->Write(); h30->Write(); h40->Write();
h100->Write(); h100->Write(); h101->Write(); h102->Write(); h103->Write(); h104->Write(); h105->Write(); h106->Write(); h107->Write();h108->Write();  h109->Write(); h110->Write(); h111->Write(); h112->Write(); h113->Write(); h114->Write(); h115->Write();
file.Close();
myout.close();

} // Loop