#define  TP_cxx    
#include "TP.h"  
#include <TH2.h>   
#include <TStyle.h> 
#include <TCanvas.h> 
#include <string>
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>

void TP::Loop() 

{

  if (fChain == 0) return;

  TH1F * TagInvMass = new TH1F("TagInvMass","TagInvMass",100,2000.,4000.);
  TH1F * ProbeInvMass = new TH1F("ProbeInvMass","ProbeInvMass",100,2000.,4000.);

  TH2F * hex = new TH2F("hex","hex",50,-3.,3.,50,-3.2,3.2);
  TH2F * hall = new TH2F("hall","hall",50,-3.,3.,50,-3.2,3.2);
  TH1F * hptallbb = new TH1F("hptallbb","hptallbb",21,0.,20);
  TH1F * hptallbe = new TH1F("hptallbe","hptallbe",21,0.,20);
   TH1F * hdr  = new TH1F("hdr","dr",100,0.,10.); 
    TH1F * hdr2  = new TH1F("hdr2","dr2",100,0.,10.); 
  TH1F * hptdr   = new TH1F("hptdr","hptdr",21,0.,20);
  TH1F * hdrcut  = new TH1F("hpdrcut","All dr cuts",100,0.,1.); 
  TH1F * hptalle = new TH1F("hptalle","hptalle",21,0.,20);
  // 
  TH1F * stacoCombComb_jpsiInvmass_os_am_aa_cut3 =new TH1F("stacoCombComb_jpsiInvmass_os_am_aa_cut3","stacoCombComb_jpsiInvmass_os_am_aa_cut3", 20, 2.86 , 3.34);
  // 

 int nptbins2  = 11; 
  float ptbin2[12];
   ptbin2[0]=1.000;
 ptbin2[0]=1.000;  
 ptbin2[1]=5.000;  
 ptbin2[2]=7.000;  
 ptbin2[3]=10.000; 
 ptbin2[4]=14.000; 
 ptbin2[5]=18.000; 
 ptbin2[6]=25.000; 
 ptbin2[7]=30.000; 
 ptbin2[8]=35.000; 
 ptbin2[9]=40.000; 
 ptbin2[10]=45.000;
	ptbin2[11]=50.000;
 
  TH1F * hptall = new TH1F("hptall","hptall",nptbins2,ptbin2); 
  TH1F * hptselall = new TH1F("hptselall","hptselall",nptbins2,ptbin2); 
  TH1F * hptb   = new TH1F("hptb","hptb",nptbins2,ptbin2); 
  TH1F * hpte   = new TH1F("hpte","hpte",nptbins2,ptbin2); 
  TH1F * hptmsallqp = new TH1F("hptmsallqp","hptmsallqp",nptbins2,ptbin2); 
  TH1F * hptmsallqn = new TH1F("hptmsallqn","hptmsallqn",nptbins2,ptbin2); 
  TH1F * hptmsselbqp = new TH1F("hptmsselbqp","hptmsselbqp",nptbins2,ptbin2); 
  TH1F * hptmsselbqn = new TH1F("hptmsselbqn","hptmsselbqn",nptbins2,ptbin2); 
  TH1F * hptmsbqp = new TH1F("hptmsbqp","hptmsbqp",nptbins2,ptbin2); 
  TH1F * hptmsbqn = new TH1F("hptmsbqn","hptmsbqn",nptbins2,ptbin2); 
  TH1F * hptmsall     = new TH1F("hptmsall","hptmsall",nptbins2,ptbin2); 
  TH1F * hptmsallqnp = new TH1F("hptmsallqnp","hptmsallqnp",nptbins2,ptbin2); 
  TH1F * hptmsallqnn = new TH1F("hptmsallqnn","hptmsallqnn",nptbins2,ptbin2); 
  TH1F * hptmsallqpnp = new TH1F("hptmsallqpnp","hptmsallqpnp",nptbins2,ptbin2); 
  TH1F * hptmsallqpnn = new TH1F("hptmsallqpnn","hptmsallqpnn",nptbins2,ptbin2); 
  TH1F * hptmsallqnnp = new TH1F("hptmsallqnnp","hptmsallqnnp",nptbins2,ptbin2); 
  TH1F * hptmsallqnnn = new TH1F("hptmsallqnnn","hptmsallqnnn",nptbins2,ptbin2); 
  TH1F * hptmsb      = new TH1F("hptmsb","hptmsb",nptbins2,ptbin2); 
  TH1F * hptmsbqnp = new TH1F("hptmsbqnp","hptmsbqnp",nptbins2,ptbin2); 
  TH1F * hptmsbqnn = new TH1F("hptmsbqnn","hptmsbqnn",nptbins2,ptbin2); 
  TH1F * hptmsbqpnp      = new TH1F("hptmsbqpnp","hptmsbqpnp",nptbins2,ptbin2); 
  TH1F * hptmsbqpnn      = new TH1F("hptmsbqpnn","hptmsbqpnn",nptbins2,ptbin2); 
  TH1F * hptmsbqnnp      = new TH1F("hptmsbqnnp","hptmsbqnnp",nptbins2,ptbin2); 
  TH1F * hptmsbqnnn      = new TH1F("hptmsbqnnn","hptmsbqnnn",nptbins2,ptbin2); 
  TH1F * hptmse          = new TH1F("hptmse","hptmse",nptbins2,ptbin2); 
  TH1F * hptmseqnp = new TH1F("hptmseqnp","hptmseqnp",nptbins2,ptbin2); 
  TH1F * hptmseqnn = new TH1F("hptmseqnn","hptmseqnn",nptbins2,ptbin2); 
  TH1F * hptmseqpnp      = new TH1F("hptmseqpnp","hptmseqpnp",nptbins2,ptbin2); 
  TH1F * hptmseqpnn      = new TH1F("hptmseqpnn","hptmseqpnn",nptbins2,ptbin2); 
  TH1F * hptmseqnnp      = new TH1F("hptmseqnnp","hptmseqnnp",nptbins2,ptbin2); 
  TH1F * hptmseqnnn      = new TH1F("hptmseqnnn","hptmseqnnn",nptbins2,ptbin2); 
  TH1F * hptmsselall = new TH1F("hptmsselall","hptmsselall",nptbins2,ptbin2);
  TH1F * hptmsselallqnp = new TH1F("hptmselallqnp","hptmselallqnp",nptbins2,ptbin2); 
  TH1F * hptmsselallqnn = new TH1F("hptmselallqnn","hptmselallqnn",nptbins2,ptbin2); 
  TH1F * hptmsselallqpnp = new TH1F("hptmsselallqpnp","hptmsselallqpnp",nptbins2,ptbin2); 
  TH1F * hptmsselallqpnn = new TH1F("hptmsselallqpnn","hptmsselallqpnn",nptbins2,ptbin2); 
  TH1F * hptmsselallqnnp = new TH1F("hptmsselallqnnp","hptmsselallqnnp",nptbins2,ptbin2); 
  TH1F * hptmsselallqnnn = new TH1F("hptmsselallqnnn","hptmsselallqnnn",nptbins2,ptbin2); 
  TH1F * hptmsselb   = new TH1F("hptmsselb","hptmsselb",nptbins2,ptbin2);
  TH1F * hptmsselbqnp = new TH1F("hptmsselbqnp","hptmsselbqnp",nptbins2,ptbin2); 
  TH1F * hptmsselbqnn = new TH1F("hptmsselbqnn","hptmsselbqnn",nptbins2,ptbin2); 
  TH1F * hptmsselbqpnp = new TH1F("hptmsselbqpnp","hptmsselbqpnp",nptbins2,ptbin2); 
  TH1F * hptmsselbqpnn = new TH1F("hptmsselbqpnn","hptmsselbqpnn",nptbins2,ptbin2); 
  TH1F * hptmsselbqnnp = new TH1F("hptmsselbqnnp","hptmsselbqnnp",nptbins2,ptbin2); 
  TH1F * hptmsselbqnnn = new TH1F("hptmsselbqnnn","hptmsselbqnnn",nptbins2,ptbin2); 
  TH1F * hptmssele   = new TH1F("hptmssele","hptmssele",nptbins2,ptbin2);
  TH1F * hptmsseleqnp = new TH1F("hptmsseleqnp","hptmsseleqnp",nptbins2,ptbin2); 
  TH1F * hptmsseleqnn = new TH1F("hptmsseleqnn","hptmsseleqnn",nptbins2,ptbin2); 
  TH1F * hptmsseleqpnp = new TH1F("hptmsseleqpnp","hptmsseleqpnp",nptbins2,ptbin2); 
  TH1F * hptmsseleqpnn = new TH1F("hptmsseleqpnn","hptmsseleqpnn",nptbins2,ptbin2); 
  TH1F * hptmsseleqnnp = new TH1F("hptmsseleqnnp","hptmsseleqnnp",nptbins2,ptbin2); 
  TH1F * hptmsseleqnnn = new TH1F("hptmsseleqnnn","hptmsseleqnnn",nptbins2,ptbin2); 
  TH1F * hptselb   = new TH1F("hptselb","hptselb",nptbins2,ptbin2); 
  TH1F * hptsele   = new TH1F("hptsele","hptsele",nptbins2,ptbin2); 
  TH1F * hptselalldr = new TH1F("hptselalldr","hptselalldr",nptbins2,ptbin2); 
  TH1F * hptselbdr   = new TH1F("hptselbdr","hptselbdr",nptbins2,ptbin2); 
  TH1F * hptseledr   = new TH1F("hptseledr","hptseledr",nptbins2,ptbin2); 
  TH1F * hptmsselbdr   = new TH1F("hptmsselbdr","hptmsselbdr",nptbins2,ptbin2); 
  TH1F * hptmsseledr   = new TH1F("hptmsseledr","hptmsseledr",nptbins2,ptbin2); 
    

  int netabin = 60;

  TH1F * hetaall   = new TH1F("hetaall","hetaall",netabin,-3.05,3.05);
  TH1F * hetaallqnp   = new TH1F("hetaallqnp","hetaallqnp",netabin,-3.05,3.05);
  TH1F * hetaallqnn   = new TH1F("hetaallqnn","hetaallqnn",netabin,-3.05,3.05);
  TH1F * hetaallqpnp   = new TH1F("hetaallqpnp","hetaallqpnp",netabin,-3.05,3.05);
  TH1F * hetaallqpnn   = new TH1F("hetaallqpnn","hetaallqpnn",netabin,-3.05,3.05);
  TH1F * hetaallqnnp   = new TH1F("hetaallqnnp","hetaallqnnp",netabin,-3.05,3.05);
  TH1F * hetaallqnnn   = new TH1F("hetaallqnnn","hetaallqnnn",netabin,-3.05,3.05);
  TH1F * hetab     = new TH1F("hetab","hetab",netabin,-3.05,3.05);
  TH1F * hetabqp     = new TH1F("hetabqp","hetabqp",netabin,-3.05,3.05);
  TH1F * hetabqn     = new TH1F("hetabqn","hetabqn",netabin,-3.05,3.05);
  TH1F * hetaselbqp     = new TH1F("hetaselbqp","hetaselbqp",netabin,-3.05,3.05);
  TH1F * hetaselbqn     = new TH1F("hetaselbqn","hetaselbqn",netabin,-3.05,3.05);
  TH1F * hetabqnp   = new TH1F("hetabqnp","hetabqnp",netabin,-3.05,3.05);
  TH1F * hetabqnn   = new TH1F("hetabqnn","hetabqnn",netabin,-3.05,3.05);
  TH1F * hetabqpnp   = new TH1F("hetabqpnp","hetabqpnp",netabin,-3.05,3.05);
  TH1F * hetabqpnn   = new TH1F("hetabqpnn","hetabqpnn",netabin,-3.05,3.05);
  TH1F * hetabqnnp   = new TH1F("hetabqnnp","hetabqnnp",netabin,-3.05,3.05);
  TH1F * hetabqnnn   = new TH1F("hetabqnnn","hetabqnnn",netabin,-3.05,3.05);
  TH1F * hetae     = new TH1F("hetae","hetae",netabin,-3.05,3.05);
  TH1F * hetaeqnp   = new TH1F("hetaeqnp","hetaeqnp",netabin,-3.05,3.05);
  TH1F * hetaeqnn   = new TH1F("hetaeqnn","hetaeqnn",netabin,-3.05,3.05);
  TH1F * hetaeqpnp   = new TH1F("hetaeqpnp","hetaeqpnp",netabin,-3.05,3.05);
  TH1F * hetaeqpnn   = new TH1F("hetaeqpnn","hetaeqpnn",netabin,-3.05,3.05);
  TH1F * hetaeqnnp   = new TH1F("hetaeqnnp","hetaeqnnp",netabin,-3.05,3.05);
  TH1F * hetaeqnnn   = new TH1F("hetaeqnnn","hetaeqnnn",netabin,-3.05,3.05);
  TH1F * hetaselall= new TH1F("hetaselall","hetaselall",netabin,-3.05,3.05);
  TH1F * hetaselallqnp   = new TH1F("hetaselallqnp","hetaselallqnp",netabin,-3.05,3.05);
  TH1F * hetaselallqnn   = new TH1F("hetaselallqnn","hetaselallqnn",netabin,-3.05,3.05);
  TH1F * hetaselallqpnp   = new TH1F("hetaselallqpnp","hetaselallqpnp",netabin,-3.05,3.05);
  TH1F * hetaselallqpnn   = new TH1F("hetaselallqpnn","hetaselallqpnn",netabin,-3.05,3.05);
  TH1F * hetaselallqnnp   = new TH1F("hetaselallqnnp","hetaselallqnnp",netabin,-3.05,3.05);
  TH1F * hetaselallqnnn   = new TH1F("hetaselallqnnn","hetaselallqnnn",netabin,-3.05,3.05);
  TH1F * hetaselb  = new TH1F("hetaselb","hetaselb",netabin,-3.05,3.05);
  TH1F * hetaselbqnp   = new TH1F("hetaselbqnp","hetaselbqnp",netabin,-3.05,3.05);
  TH1F * hetaselbqnn   = new TH1F("hetaselbqnn","hetaselbqnn",netabin,-3.05,3.05);
  TH1F * hetaselbqpnp   = new TH1F("hetaselbqpnp","hetaselbqpnp",netabin,-3.05,3.05);
  TH1F * hetaselbqpnn   = new TH1F("hetaselbqpnn","hetaselbqpnn",netabin,-3.05,3.05);
  TH1F * hetaselbqnnp   = new TH1F("hetaselbqnnp","hetaselbqnnp",netabin,-3.05,3.05);
  TH1F * hetaselbqnnn   = new TH1F("hetaselbqnnn","hetaselbqnnn",netabin,-3.05,3.05);
  TH1F * hetasele  = new TH1F("hetasele","hetasele",netabin,-3.05,3.05);
  TH1F * hetaseleqnp   = new TH1F("hetaseleqnp","hetaseleqnp",netabin,-3.05,3.05);
  TH1F * hetaseleqnn   = new TH1F("hetaseleqnn","hetaseleqnn",netabin,-3.05,3.05);
  TH1F * hetaseleqpnp   = new TH1F("hetaseleqpnp","hetaseleqpnp",netabin,-3.05,3.05);
  TH1F * hetaseleqpnn   = new TH1F("hetaseleqpnn","hetaseleqpnn",netabin,-3.05,3.05);
  TH1F * hetaseleqnnp   = new TH1F("hetaseleqnnp","hetaseleqnnp",netabin,-3.05,3.05);
  TH1F * hetaseleqnnn   = new TH1F("hetaseleqnnn","hetaseleqnnn",netabin,-3.05,3.05);
  TH1F * hetaselalldr = new TH1F("hetaselalldr","hetaselalldr",netabin,-3.05,3.05);
  TH1F * hetaselbdr = new TH1F("hetaselbdr","hetaselbdr",netabin,-3.05,3.05);
  TH1F * hetaseledr = new TH1F("hetaseledr","hassled",netabin,-3.05,3.05);


  int nphibin = 20;

  TH1F * hphiall   = new TH1F("hphiall","hphiall",nphibin,-3.14,3.14);
  TH1F * hphiallqnp   = new TH1F("hphiallqnp","hphiallqnp",nphibin,-3.14,3.14);
  TH1F * hphiallqnn   = new TH1F("hphiallqnn","hphiallqnn",nphibin,-3.14,3.14);
  TH1F * hphiallqpnp   = new TH1F("hphiallqpnp","hphiallqpnp",nphibin,-3.14,3.14);
  TH1F * hphiallqpnn   = new TH1F("hphiallqpnn","hphiallqpnn",nphibin,-3.14,3.14);
  TH1F * hphiallqnnp   = new TH1F("hphiallqnnp","hphiallqnnp",nphibin,-3.14,3.14);
  TH1F * hphiallqnnn   = new TH1F("hphiallqnnn","hphiallqnnn",nphibin,-3.14,3.14);
  TH1F * hphie     = new TH1F("hphie","hphie",nphibin,-3.14,3.14);
  TH1F * hphieqnp   = new TH1F("hphieqnp","hphieqnp",nphibin,-3.14,3.14);
  TH1F * hphieqnn   = new TH1F("hphieqnn","hphieqnn",nphibin,-3.14,3.14);
  TH1F * hphieqpnp   = new TH1F("hphieqpnp","hphieqpnp",nphibin,-3.14,3.14);
  TH1F * hphieqpnn   = new TH1F("hphieqpnn","hphieqpnn",nphibin,-3.14,3.14);
  TH1F * hphieqnnp   = new TH1F("hphieqnnp","hphieqnnp",nphibin,-3.14,3.14);
  TH1F * hphieqnnn   = new TH1F("hphieqnnn","hphieqnnn",nphibin,-3.14,3.14);
  TH1F * hphib     = new TH1F("hphib","hphib",nphibin,-3.14,3.14);
  TH1F * hphibqp   = new TH1F("hphibqp","hphibqp",nphibin,-3.14,3.14);
  TH1F * hphibqn   = new TH1F("hphibqn","hphibqn",nphibin,-3.14,3.14);
  TH1F * hphiselbqp   = new TH1F("hphiselbqp","hphiselbqp",nphibin,-3.14,3.14);
  TH1F * hphiselbqn   = new TH1F("hphiselbqn","hphiselbqn",nphibin,-3.14,3.14);
  TH1F * hphibqnp   = new TH1F("hphibqnp","hphibqnp",nphibin,-3.14,3.14);
  TH1F * hphibqnn   = new TH1F("hphibqnn","hphibqnn",nphibin,-3.14,3.14);
  TH1F * hphibqpnp   = new TH1F("hphibqpnp","hphibqpnp",nphibin,-3.14,3.14);
  TH1F * hphibqpnn   = new TH1F("hphibqpnn","hphibqpnn",nphibin,-3.14,3.14);
  TH1F * hphibqnnp   = new TH1F("hphibqnnp","hphibqnnp",nphibin,-3.14,3.14);
  TH1F * hphibqnnn   = new TH1F("hphibqnnn","hphibqnnn",nphibin,-3.14,3.14);
  TH1F * hphiselall = new TH1F("hphiselall","hphiselall",nphibin,-3.14,3.14);
  TH1F * hphiselallqnp   = new TH1F("hphiselallqnp","hphiselallqnp",nphibin,-3.14,3.14);
  TH1F * hphiselallqnn   = new TH1F("hphiselallqnn","hphiselallqnn",nphibin,-3.14,3.14);
  TH1F * hphiselallqpnp   = new TH1F("hphiselallqpnp","hphiselallqpnp",nphibin,-3.14,3.14);
  TH1F * hphiselallqpnn   = new TH1F("hphiselallqpnn","hphiselallqpnn",nphibin,-3.14,3.14);
  TH1F * hphiselallqnnp   = new TH1F("hphiselallqnnp","hphiselallqnnp",nphibin,-3.14,3.14);
  TH1F * hphiselallqnnn   = new TH1F("hphiselallqnnn","hphiselallqnnn",nphibin,-3.14,3.14);
  TH1F * hphiselb  = new TH1F("hphiselb","hphiselb",nphibin,-3.14,3.14);
  TH1F * hphiselbqnp   = new TH1F("hphiselbqnp","hphiselbqnp",nphibin,-3.14,3.14);
  TH1F * hphiselbqnn   = new TH1F("hphiselbqnn","hphiselbqnn",nphibin,-3.14,3.14);
  TH1F * hphiselbqpnp   = new TH1F("hphiselbqpnp","hphiselbqpnp",nphibin,-3.14,3.14);
  TH1F * hphiselbqpnn   = new TH1F("hphiselbqpnn","hphiselbqpnn",nphibin,-3.14,3.14);
  TH1F * hphiselbqnnp   = new TH1F("hphiselbqnnp","hphiselbqnnp",nphibin,-3.14,3.14);
  TH1F * hphiselbqnnn   = new TH1F("hphiselbqnnn","hphiselbqnnn",nphibin,-3.14,3.14);
  TH1F * hphisele  = new TH1F("hphisele","hphisele",nphibin,-3.14,3.14);
  TH1F * hphiseleqnp   = new TH1F("hphiseleqnp","hphiseleqnp",nphibin,-3.14,3.14);
  TH1F * hphiseleqnn   = new TH1F("hphiseleqnn","hphiseleqnn",nphibin,-3.14,3.14);
  TH1F * hphiseleqpnp   = new TH1F("hphiseleqpnp","hphiseleqpnp",nphibin,-3.14,3.14);
  TH1F * hphiseleqpnn   = new TH1F("hphiseleqpnn","hphiseleqpnn",nphibin,-3.14,3.14);
  TH1F * hphiseleqnnp   = new TH1F("hphiseleqnnp","hphiseleqnnp",nphibin,-3.14,3.14);
  TH1F * hphiseleqnnn   = new TH1F("hphiseleqnnn","hphiseleqnnn",nphibin,-3.14,3.14);
  TH1F * hphiselalldr = new TH1F("hphiselalldr","hphiselalldr",nphibin,-3.14,3.14);
  TH1F * hphiselbdr  = new TH1F("hphiselbdr","hphiselbdr",nphibin,-3.14,3.14);
  TH1F * hphiseledr  = new TH1F("hphiseledr","hphiseledr",nphibin,-3.14,3.14);
  
  my_h1= new TFile ("TPout_EFMU4_EFMU4_2011_CombComb_staco.root", "RECREATE");
  
  int nphibins = 5;
 
  ///////////////////////////
  ///////////////////////////
  //MU0
  int nptbins  = 11; 
  float ptbin[12]={1,5,7,10,14,18,25,30,35,40,45,50};
  int netabins = 18; 
  float etabin[19]={-2.4,-2.0,-1.7,-1.5,-1.3,-1.05,-0.8,-0.5,-0.1,0,0.1,0.5,0.8,1.05,1.3,1.5,1.7,2.0,2.4}; 
  
  
  TH2F * hmapall    = new TH2F("hmapall","hmapall",netabins,etabin,nptbins,ptbin);                       
  TH2F * hmapselall = new TH2F("hmapselall","hmapselall",netabins,etabin,nptbins,ptbin);                 
  TH2F * hmapselalleff = new TH2F("hmapselalleff","hmapselalleff",netabins,etabin,nptbins,ptbin);        
  TH2F * hmapb      = new TH2F("hmapb","hmapb",netabins,-1.05,1.05,nphibins,-3.14,3.14);                 
  TH2F * hmape      = new TH2F("hmape","hmape",netabins,-1.05,1.05,nphibins,-3.14,3.14);                 
  TH2F * hmapselb   = new TH2F("hmapselb","hmapselb",netabins,-1.05,1.05,nphibins,-3.14,3.14);     	    
  TH2F * hmapsele   = new TH2F("hmapsele","hmapsele",10,-2.4165,2.4165,16,-3.14,3.14);		    
  TH2F * hmapselalldr = new TH2F("hmapselalldr","hmapselalldr",netabins,-1.05,1.05,nphibins,-3.14,3.14); 
  TH2F * hmapselbdr   = new TH2F("hmapselbdr","hmapselbdr",netabins,-1.05,1.05,nphibins,-3.14,3.14);     
  TH2F * hmapdr     = new TH2F("hmapdr","hmapdr",netabins,-1.05,1.05,nphibins,-3.14,3.14);       	    
  TH2F * hmapalle   = new TH2F("hmapalle","hmapalle",10,-2.4165,2.4165,16,-3.14,3.14);		    
  TH2F * hmapseledr = new TH2F("hmapseledr","hmapseledr",10,-2.4165,2.4165,16,-3.14,3.14);		    
  
  ofstream log;                                
  log.open("logTP.txt");		      
  bool DEBUG = false;    		      
						      
   Long64_t nentries = fChain->GetEntries();   
   Long64_t nbytes = 0, nb = 0;                 
   					       
   int ntotal = 0;			       
   int ntag   = 0;			       
   int nprobe = 0;                              
   int nprobetrigger = 0;		       
   int npass =  0;			       
   int nbarrel = 0;			       
   int n_1 = 0;         
   int n_2 = 0;	       
   int n_3 = 0;	       
   int n_barrel=0;      
   int n_endcap = 0;     
   int n_trig_barrel=0;  
   int n_trig_endcap=0;  
   int n_dr_barrel=0;    
   int n_dr_endcap=0;    
   int n_same = 0;       
			        
   int n_syst[20];       
   
   bool doPIcuts   = kFALSE;    
   bool IsMC = kTRUE; 
 
   bool DoCombComb = kTRUE ;
   bool DoCombAny  = kFALSE;
   bool doProbeSegmentTagged = kFALSE; 
   bool doProbeComb = kFALSE; 
   
   // Select kind of T&P (only one type a time!!) 
   bool DoL1 = kFALSE;  
   bool DoEF = kTRUE; 
   bool DoMix = kFALSE; 
   bool DoRev = kFALSE; 


   if ( (DoL1 + DoEF + DoMix + DoRev) > 1) 
     { 
       if (DEBUG) log << "More than 1 kind of TP: Exit !!"<<endl;
       return;
     }

   if (DEBUG) log << "number of entries = "<<nentries<<endl;

   for (int h=0; h<20; h++) n_syst[h]=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(runNumber==153565 || runNumber==153599) continue;

       bool IsCombComb = 0;
       if (mu1author ==6 && mu2author == 6) IsCombComb=1;
       bool IsCombAny = 0;
       if ( (mu1author ==6 && mu2author == 7) || (mu1author ==7  && mu2author == 6) || IsCombComb )  IsCombAny = 1; 
       if (DoCombComb && !IsCombComb) continue; // staco CombComb
       if (DoCombAny  && !IsCombAny)  continue; // staco CombAny
       
       
       // Jpsi event selection cuts
       bool IsJpsi  = 1;
       float deltaR = sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
       bool EtaCut  = 0;
       hdr->Fill(deltaR);
       float deltaR2=0.,dPhi,dEta;
      dPhi=fabs(mu1phi-mu2phi);
      dEta=fabs(mu1eta-mu2eta);
       if(dPhi>3.14159264) dPhi=2*3.14159264-dPhi;
       deltaR2=  sqrt(dPhi*dPhi+dEta*dEta);
       hdr2->Fill(deltaR2);

       if (fabs(mu1eta)<2.4 && fabs(mu2eta)<2.4) EtaCut = 1;
       
       if (! IsMC)
	  {
	    IsJpsi = 0;
	  if (invmass >2000 && invmass < 4000 && q == 0 && EtaCut && deltaR2>0.3 && deltaR2<2.5 ) // JpsiSelection 
	      {   
		IsJpsi = 1;
	      }
	  }

       // Here start Tag and Probe
       if (IsJpsi)
	 {
	   if (DEBUG) log << " ** Dimuon event found on event ** "<<jentry<<endl;
	   ntotal++; // new dimuon events

	   for (int muonindex= 1; muonindex <3; muonindex++)
	     {
	       if (DEBUG) log << " ** Muon Index = "<< muonindex<<endl;
	       
	       float tageta;
	       float tagphi;
	       float tagpt;   
	       float tagz0;   
	       float tagdr; 
	       float tagdrL1;
	       float tagdrEF;
	       float tagiso;
	       bool  isBadTagMuon= 0;
	     
	       if (DoL1 || DoRev)
		 {
		   if (muonindex == 1)
		     {
		       tagpt    = mu1pt;
		       tagz0    = mu1z0;
		       tagdr    = mu1EFMU4dr;
		       tagdrL1  = mu1EFMU4dr;
		       tagdrEF  = mu1EFMU20dr;
		       tagiso  = mu1_ptIsol40;
		       		       
		       if (fabs(mu1etas) < 100)
			 {
			   tageta   = mu1etas;
			   tagphi   = mu1phis; 
			 }
		       
		       if (fabs(mu1etas) >= 100 && fabs(mu1etae) < 100)
			 {
			   tageta   = mu1etae;
			   tagphi   = mu1phie; 
			 }
		       
		       if (fabs(mu1etas) >= 100 && fabs(mu1etae) >= 100)
			 {
			   isBadTagMuon == 1;
			   tageta   = mu1eta;
			   tagphi   = mu1phi; 
			 }
		     }
		   
		   if (muonindex == 2)
		     {
		       tagpt    = mu2pt;
		       tagz0    = mu2z0;
		       tagdr    = mu2EFMU4dr;
		       tagdrL1  = mu2EFMU4dr;
		       tagdrEF  = mu2EFMU20dr;
		       tagiso   = mu2_ptIsol40;

		       if (fabs(mu1etas) < 100)
			 {
			tageta   = mu2etas;
			tagphi   = mu2phis; 
			 }
		       
		       if (fabs(mu2etas) >= 100 && fabs(mu2etae) < 100)
			 {
			   tageta   = mu2etae;
			   tagphi   = mu2phie; 
			 }
		       
		       if (fabs(mu1etas) >= 100 && fabs(mu1etae) >= 100)
			 {
			   isBadTagMuon == 1;
			   tageta   = mu2eta;
			   tagphi   = mu2phi; 
			 }
		     }
	       	 } // end of (if ! DoEF)

	       if (DoEF || DoMix)
		 {
		   if (muonindex == 1)
		     {
		       tagpt    = mu1pt;
		       tagz0    = mu1z0;
		       tagdr    = mu1EFMU4dr;
		       tagdrL1  = mu1L1MU0dr;
		       tagdrEF  = mu1EFMU4dr;
		       tageta   = mu1eta;
		       tagphi   = mu1phi; 
		       tagiso   = mu1_ptIsol40;
		       //cout<<"mu1EFMU6dr "<<mu1EFMU6dr<<endl;
		     }
		   if (muonindex == 2)
		     {
		       tagpt    = mu2pt;
		       tagz0    = mu2z0;
		       tagdr    = mu2EFMU4dr;
		       tagdrL1  = mu2L1MU0dr;
		       tagdrEF  = mu2EFMU4dr;
		       tageta   = mu2eta;
		       tagphi   = mu2phi; 
		       tagiso   = mu2_ptIsol40;
		     }
		 } // end of (if ! DoEF)


	       bool IsTag = 0;
	       bool HasTagQuality = MuonQuality(muonindex);
	       float tagdrcut;
	      
	       ///deltaR cut 
	       if (  DoEF || DoMix) tagdrcut = 0.05; // EF DR cut
	       if (  DoL1 || DoRev) tagdrcut = 0.4;  // L1 DR cut
 
		 //tagdrcut = SearchDrCut(tagpt, tageta, tagphi, drmap);
	      
	       bool SameTriggerIndexL1 = false;
	       bool SameTriggerIndexEF = false;
	       bool CannotBeTagL1 = false;
	       bool CannotBeTagEF = false;

	       if (DEBUG) log<< " mu1 eta=" << mu1eta << " extr=" << mu1etae << " sp=" << mu1etas << endl;
	       if (DEBUG) log<< " mu2 eta=" << mu2eta << " extr=" << mu2etae << " sp=" << mu2etas << endl;
	       if (DEBUG) log << " mu1 EF trigger index = " << mu1EFMU4index <<" and DR is "<<mu1EFMU4index<<endl;
	       if (DEBUG) log << " mu2 EF trigger index = " << mu2EFMU4index <<" and DR is "<<mu2EFMU4index<<endl;
	       if (DEBUG) log << " mu1 L1 trigger index = " << mu1L1MU0index <<" and DR is "<<mu1L1MU0dr<<endl;
	       if (DEBUG) log << " mu2 L1 trigger index = " << mu2L1MU0index <<" and DR is "<<mu2L1MU0dr<<endl;
	       
	       if ( mu1L1MU0index == mu2L1MU0index ) SameTriggerIndexL1 = true;
	       if ( mu1EFMU4index == mu2EFMU4index ) SameTriggerIndexEF = true;
	   
	       if (SameTriggerIndexL1 == true) 
		 {
		   // If both muons inside the Dimuon could be matched with 
		   // the same Trigger Index Object then the muon (mu1 or mu2) 
		   // with the 
		   // higher DR is discarded as a TAG and is acceptes as a
		   // Probe but IS NOT MATCHED!! 

		   if (DEBUG) log << " This Event has same L1 trigger object index"<<endl;
		   float OtherDr ;

		   if (muonindex == 1 ) OtherDr = mu2L1MU0dr ;
		   if (muonindex == 2 ) OtherDr = mu1L1MU0dr ;

		   if (DEBUG) log << " tagdrL1 = "   << tagdrL1 <<endl;
		   if (DEBUG) log << " OtherDrL1 = " << OtherDr <<endl;

		   if (tagdrL1 > OtherDr ) CannotBeTagL1 = true;
		 }

	       if (CannotBeTagL1) if (DEBUG) log << " This tag for L1 has been rejected "<<endl;

	       if (SameTriggerIndexEF == true) 
	 {	 
		   if (DEBUG) log << " This Event has same EF trigger object index"<<endl;
		   float OtherDr ;

		   if (muonindex == 1 ) OtherDr = mu2EFMU4dr ;
		   if (muonindex == 2 ) OtherDr = mu1EFMU4dr ;

		   if (DEBUG) log << " tagdrEF = "   << tagdrEF<<endl;
		   if (DEBUG) log << " OtherDrEF = " << OtherDr<<endl;

		   if (tagdrEF > OtherDr ) CannotBeTagEF = true;
		 }

	       if (CannotBeTagEF) if (DEBUG) log << " This tag for EF has been rejected "<<endl;
	       
	       bool HasTriggerInEventL1 = false;
	       bool HasTriggerInEventEF = false;

	       if (   muEFMU4pass >0 ) HasTriggerInEventEF = true;
	       if (   muL1MU0pass >0 ) HasTriggerInEventL1 = true;

	       if (DEBUG) log << " L1 decision = "<< HasTriggerInEventL1<<endl;
	       if (DEBUG) log << " EF decision = "<< HasTriggerInEventEF<<endl;

	       if ( DoL1 || DoRev)
		 {
		   if ( HasTagQuality && HasTriggerInEventL1 
			&& fabs(tagdr) < tagdrcut && isBadTagMuon == 0 
                 && tagpt > 1000 && tagiso/tagpt < 0.2  // added
			&& !CannotBeTagL1 ) IsTag = 1;
		 }

	       if ( DoEF || DoMix) 
		 {
		   if ( HasTagQuality && HasTriggerInEventEF 
			&& fabs(tagdr) < tagdrcut && isBadTagMuon == 0 
                 && tagpt > 1000 && tagiso/tagpt < 0.2  // added
			&& !CannotBeTagEF && !CannotBeTagL1 ) IsTag = 1; 
		 } 

	       if (IsTag) 
		 {

		   if (DEBUG) log << " TAG is ok "<<endl; 
		   TagInvMass->Fill(invmass); 
		   ntag++;
		   float probeetaID; 
		   float probephiID; 
		   if (muonindex == 1) probeetaID = mu2eta; 
		   if (muonindex == 2) probeetaID = mu1eta; 
		   if (muonindex == 1) probephiID = mu2phi;
		   if (muonindex == 2) probephiID = mu1phi;
		   float probeeta ;
		   float probephi ;

		   float probept;   
		   float probez0;   
		   float probedr;
		   float probedrL1;
		   int   ProbeIndex;
		   bool isBadProbeMuon= 0;
		   int probeAuthor = -100;
		   int probeNtrt;
		   int probeq;
		   if ( DoL1 || DoMix )
		     {
		       if (muonindex == 2)
			 {
			   probept    = mu1pt;
			   probez0    = mu1z0;
			   probedr    = mu1L1MU0dr ;
			   probedrL1  = mu1L1MU0dr ;
			   ProbeIndex = 1;
			   probeAuthor= mu1author; 
			   probeNtrt  = mu1_nTRT;
			   probeq  = mu1q;
			   //probeNtrt  = mu1Ntrt;

			   if (fabs(mu1etas) < 100)
			     {
			       probeeta   = mu1etas;
			       probephi   = mu1phis; 
			     }
			   
			   if (fabs(mu1etas) >= 100 && fabs(mu1etae) < 100)
			     {
			       probeeta   = mu1etae;
			       probephi   = mu1phie; 
			     }
			   
			   if (fabs(mu1etas) >= 100 && fabs(mu1etae) >= 100)
			     {
			       // isBadProbeMuon= 1;
			       probeeta   = mu1eta;
			       probephi   = mu1phi; 
			     }
			 }
		       
		       if (muonindex == 1)
			 {
			   probept    = mu2pt;
			   probez0    = mu2z0; 
			   probedr    = mu2L1MU0dr; 
			   probedrL1  = mu2L1MU0dr; 
			   ProbeIndex = 2;
			   probeAuthor= mu2author;
			   probeNtrt  = mu2_nTRT;
			   probeq     = mu2q;
			   //probeNtrt  = mu2Ntrt;
			   
			   if (fabs(mu2etas) < 100)
			     {
			       probeeta   = mu2etas;
			       probephi   = mu2phis; 
			     }
			   
			   if (fabs(mu2etas) >= 100 && fabs(mu2etae) < 100)
			     {
			       probeeta   = mu2etae;
			       probephi   = mu2phie; 
			     }
			   
			   if (fabs(mu2etas) >= 100 && fabs(mu2etae) >= 100) 
			     {
			       //			       isBadProbeMuon= 1;
			       probeeta   = mu2eta;
			       probephi   = mu2phi; 
			     }  
			 }
		     } // end if !DoEF
		   
		   if (DoEF || DoRev) 
		    { 
		       if (muonindex == 1) 
			 {
		       probept    = mu2pt;       
		       probez0    = mu2z0;
		       probedr    = mu2EFMU4dr ; 
		       probedrL1  = mu2L1MU0dr ;
		       probeeta   = mu2eta;
		       probephi   = mu2phi;
		       probeAuthor= mu2author;
		       probeNtrt  = mu2_nTRT; 
		       probeq     = mu2q;
		       //probeNtrt=mu2Ntrt; 
		     }
		   if (muonindex == 2) 
		     {
		       probept    = mu1pt; 
		       probez0    = mu1z0; 
		       probedr    = mu1EFMU4dr ; 
		       probedrL1  = mu1L1MU0dr ; 
		       probeeta   = mu1eta;
		       probephi   = mu1phi;
		       probeAuthor= mu1author; 
		       probeNtrt  = mu1_nTRT;
		       probeq     = mu1q;
		       //probeNtrt=mu1Ntrt;
		     }

		 } // end of (if ! DoEF)

		
		   bool IsProbe = 0;   
		   bool HasProbeQuality = MuonQuality(ProbeIndex);
		   float probedrcut;
		   ///probe DeltaR cut	
		   float probedrcutL1=0.4;
		   float probedrcutEF=0.05;
		   // [-1.3 to -1.1], [-0.1 to 0.1] and [1.1 to 1.3] 

		   ///p cuts                      
		   float probetheta=0.;
		   float probep=0.;
		   probetheta=2*atan(exp(-1*probeetaID));
		   probep=probept/sin(probetheta);
		   bool pcut=false;
		   if(fabs(probeetaID)<=2.5 && probep>3000.) pcut = true;

		   ///trt cut 
		   bool trtCut=false;
		   if(   (fabs(probeetaID)<=2. && probeNtrt>10) ||
			 fabs(probeetaID)>2
			 ) trtCut=true;
		   
		   if (DEBUG) log<<"probetheta="<<probetheta<<" probeeta="<<probeetaID<<" probep="<<probep<<" probept="<<probept
		      << " probeAuthor=" << probeAuthor<<endl;
		  
		   ///runNmber cut
		   bool runNumberCut=true;
		   
		   ///Exclusion regions
		   bool excludeEtaRegion = false; 
		   if((probeetaID>-1.3 && probeetaID<-1.1)  
		      || (probeetaID>-0.1 && probeetaID<0.1) 
		      || (probeetaID>1.1 && probeetaID<1.3) ) excludeEtaRegion = true;
		   //if((probeetaID>-1.5 && probeetaID<-0.5) 
		   //   || (probeetaID>-0.1 && probeetaID<0.1) 
		   //   || (probeetaID>1.1 && probeetaID<1.3) ) excludeEtaRegion = true;
		   
		   if (DEBUG) log <<  " before  isBadProbeMuon=" << isBadProbeMuon << endl;
		   if(doProbeSegmentTagged && probeAuthor!=13) isBadProbeMuon = 1;
		   
		   if (DEBUG) log <<  " HasProbeQuality=" <<  HasProbeQuality << " isBadProbeMuon=" << isBadProbeMuon << " pcut="<< pcut << endl;

		   if ( HasProbeQuality && isBadProbeMuon == 0 && pcut && runNumberCut && trtCut) // && !excludeEtaRegion && probept> 8000) 
		     {
		       //cout<<"dentro if HasProbeQuality && isBadProbeMuon == 0 && pcut && runNumberCut"<<endl;
		       if (DEBUG) log << "good probe to be Tested " << endl;
		       ProbeInvMass->Fill(invmass);
		       stacoCombComb_jpsiInvmass_os_am_aa_cut3->Fill(invmass/1000.);

		         if (  DoEF || DoRev) probedrcut = 0.05; // EF DR cut
			 if (  DoL1 || DoMix) probedrcut = 0.4;  // L1 DR cut
		       
			 //			 if (probept > 10000)
			 if (true)
			   {
			     if (DEBUG) log << "Probe pt  = "<<probept<<endl;
			     if (DEBUG) log << "Probe eta = "<<probeetaID<<endl;
			     if (DEBUG) log << "Probe phi = "<<probephiID<<endl;
			   }

		       hptmsall->Fill(probept/1000.); // fill histo for barrel	 
		       hetaall->Fill(probeetaID);
		       hphiall->Fill(probephiID);
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmsallqnp->Fill(probept/1000.); 
                hetaallqnp->Fill(probeetaID); 
                hphiallqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmsallqnn->Fill(probept/1000.); 
                hetaallqnn->Fill(probeetaID); 
                hphiallqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmsallqpnp->Fill(probept/1000.); 
                hetaallqpnp->Fill(probeetaID); 
                hphiallqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmsallqpnn->Fill(probept/1000.); 
                hetaallqpnn->Fill(probeetaID); 
                hphiallqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmsallqnnp->Fill(probept/1000.); 
                hetaallqnnp->Fill(probeetaID); 
                hphiallqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmsallqnnn->Fill(probept/1000.); 
                hetaallqnnn->Fill(probeetaID); 
                hphiallqnnn->Fill(probephiID); } 
		       
		        hmapall->Fill(fabs(probeetaID),probept);
			hptall->Fill(probep/1000.);
		       if (fabs(probeetaID) < 1.05) 
			 { 
			   nprobe++;
			   hptmsb ->Fill(probept/1000.); // fill histo for barrel
			   hetab->Fill(probeetaID);
			   hphib->Fill(probephiID);
		       if (probeq>0)  {hptmsbqp->Fill(probept/1000.); 
		       hetabqp->Fill(probeetaID); 
                hphibqp->Fill(probephiID); } 
		       if (probeq<0)  {hptmsbqn->Fill(probept/1000.); 
		       hetabqn->Fill(probeetaID); 
                hphibqn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmsbqnp->Fill(probept/1000.); 
                hetabqnp->Fill(probeetaID); 
                hphibqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmsbqnn->Fill(probept/1000.); 
                hetabqnn->Fill(probeetaID); 
                hphibqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmsbqpnp->Fill(probept/1000.); 
                hetabqpnp->Fill(probeetaID); 
                hphibqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmsbqpnn->Fill(probept/1000.); 
                hetabqpnn->Fill(probeetaID); 
                hphibqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmsbqnnp->Fill(probept/1000.); 
                hetabqnnp->Fill(probeetaID); 
                hphibqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmsbqnnn->Fill(probept/1000.); 
                hetabqnnn->Fill(probeetaID); 
                hphibqnnn->Fill(probephiID); } 
		       hmapb->Fill(probeetaID,probephiID);
			 }
		       else
			 {	
			   nprobe++;
			   hptmse ->Fill(probept/1000.); // fill histo for Endcap
			   hetae->Fill(probeetaID);
			   hphie->Fill(probephiID);
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmseqnp->Fill(probept/1000.); 
                hetaeqnp->Fill(probeetaID); 
                hphieqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmseqnn->Fill(probept/1000.); 
                hetaeqnn->Fill(probeetaID); 
                hphieqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmseqpnp->Fill(probept/1000.); 
                hetaeqpnp->Fill(probeetaID); 
                hphieqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmseqpnn->Fill(probept/1000.); 
                hetaeqpnn->Fill(probeetaID); 
                hphieqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmseqnnp->Fill(probept/1000.); 
                hetaeqnnp->Fill(probeetaID); 
                hphieqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmseqnnn->Fill(probept/1000.); 
                hetaeqnnn->Fill(probeetaID); 
                hphieqnnn->Fill(probephiID); } 
			   hmape->Fill(probeetaID,probephiID);
			 }
		       
		       bool HasRequiredTriggerInEvent;
		       bool HasProbeTriggerL1 = false;
		       bool HasProbeTriggerEF = false;
		       bool HasProbeTrigger = false;

		       if (DoEF || DoRev)          HasRequiredTriggerInEvent = HasTriggerInEventEF;
		       if (DoL1 || DoMix) HasRequiredTriggerInEvent = HasTriggerInEventL1;
		       
		       
		       if (DEBUG) log << " probe dr l1 = "<<probedrL1<<endl; 
		       if (DoL1 || DoMix) if (DEBUG) log << " probe dr L1 = "<<probedr<<endl; 
		       if (DoEF || DoRev) if (DEBUG) log << " probe dr EF = "<<probedr<<endl; 

		       if (DEBUG) log << " CannotBeTagL1 = "<<CannotBeTagL1<<endl; 
		       if (DEBUG) log << " CannotBeTagEF = "<<CannotBeTagEF<<endl;
		       
		       if ( (DoL1 || DoMix) && SameTriggerIndexL1 == true) if (DEBUG) log << "This event must be rejected and non efficient as probe"<<endl;
		       
		       if ( fabs(probedrL1) < probedrcutL1 && SameTriggerIndexL1 != true) HasProbeTriggerL1 = true;

		       if (DEBUG) log << " HasProbeTriggerL1 = "<<HasProbeTriggerL1<<endl;

		       if ( DoEF  && SameTriggerIndexEF == true) if (DEBUG) log << "This event must be rejected in EF and non efficient as probe"<<endl;
		       if ( DoEF && fabs(probedr) <  probedrcut && HasProbeTriggerL1 && SameTriggerIndexEF != true)  HasProbeTriggerEF = true;
		       if ( DoRev && fabs(probedr) <  probedrcut && HasProbeTriggerL1 )  HasProbeTriggerEF = true;
		       if ( DoEF || DoRev ) if (DEBUG) log << " HasProbeTriggerEF = "<<HasProbeTriggerEF<<endl;

		       if ( (DoL1 || DoMix) && HasProbeTriggerL1) HasProbeTrigger = true;
		       if ( (DoEF || DoRev) && HasProbeTriggerEF) HasProbeTrigger = true;
		       
		       if (DEBUG) log << " Probe trigger matching decision is "<< HasProbeTrigger <<endl;
		       if ( HasProbeQuality && HasProbeTrigger)
			 {
			   
			   bool HasPassedL1 = false;

			   hptmsselall ->Fill(probept/1000.); 
			   hetaselall->Fill(probeetaID);
			   hphiselall->Fill(probephiID);
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmsselallqnp->Fill(probept/1000.); 
                hetaselallqnp->Fill(probeetaID); 
                hphiselallqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmsselallqnn->Fill(probept/1000.); 
                hetaselallqnn->Fill(probeetaID); 
                hphiselallqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmsselallqpnp->Fill(probept/1000.); 
                hetaselallqpnp->Fill(probeetaID); 
                hphiselallqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmsselallqpnn->Fill(probept/1000.); 
                hetaselallqpnn->Fill(probeetaID); 
                hphiselallqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmsselallqnnp->Fill(probept/1000.); 
                hetaselallqnnp->Fill(probeetaID); 
                hphiselallqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmsselallqnnn->Fill(probept/1000.); 
                hetaselallqnnn->Fill(probeetaID); 
                hphiselallqnnn->Fill(probephiID); } 
			   hmapselall->Fill(fabs(probeetaID),probept); 
			   //hmapselall->Fill(probeetaID,probept);
			   hptselall->Fill(probep/1000.);

			   if (fabs(probeetaID) < 1.05)
			     {   
			       nprobetrigger++; 
			       hptmsselb ->Fill(probept/1000.); // fill histo for barrel
			       hetaselb->Fill(probeetaID);
			       hphiselb->Fill(probephiID);
		       if (probeq>0)  {hptmsselbqp->Fill(probept/1000.); 
		       hetaselbqp->Fill(probeetaID); 
                hphiselbqp->Fill(probephiID); } 
		       if (probeq<0)  {hptmsselbqn->Fill(probept/1000.); 
		       hetaselbqn->Fill(probeetaID); 
                hphiselbqn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmsselbqnp->Fill(probept/1000.); 
                hetaselbqnp->Fill(probeetaID); 
                hphiselbqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmsselbqnn->Fill(probept/1000.); 
                hetaselbqnn->Fill(probeetaID); 
                hphiselbqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmsselbqpnp->Fill(probept/1000.); 
                hetaselbqpnp->Fill(probeetaID); 
                hphiselbqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmsselbqpnn->Fill(probept/1000.); 
                hetaselbqpnn->Fill(probeetaID); 
                hphiselbqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmsselbqnnp->Fill(probept/1000.); 
                hetaselbqnnp->Fill(probeetaID); 
                hphiselbqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmsselbqnnn->Fill(probept/1000.); 
                hetaselbqnnn->Fill(probeetaID); 
                hphiselbqnnn->Fill(probephiID); } 
			       hmapselb->Fill(probeetaID,probephiID);
			     }
			   if (fabs(probeetaID) >= 1.05)
			     {   
			       nprobetrigger++;
			       hptmssele ->Fill(probept/1000.); // fill histo for endcap
			       hetasele->Fill(probeetaID);
			       hphisele->Fill(probephiID);
		       if (probeq>0 && probeetaID >0 || probeq<0 && probeetaID <0)  { 
                hptmsseleqnp->Fill(probept/1000.); 
                hetaseleqnp->Fill(probeetaID); 
                hphiseleqnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0 || probeq<0 && probeetaID >0) { 
                hptmsseleqnn->Fill(probept/1000.); 
                hetaseleqnn->Fill(probeetaID); 
                hphiseleqnn->Fill(probephiID); } 
		       if (probeq>0 && probeetaID >0) { 
                hptmsseleqpnp->Fill(probept/1000.); 
                hetaseleqpnp->Fill(probeetaID); 
                hphiseleqpnp->Fill(probephiID); } 
		       if (probeq>0 && probeetaID <0) { 
                hptmsseleqpnn->Fill(probept/1000.); 
                hetaseleqpnn->Fill(probeetaID); 
                hphiseleqpnn->Fill(probephiID); } 
		       if (probeq<0 && probeetaID >0) { 
                hptmsseleqnnp->Fill(probept/1000.); 
                hetaseleqnnp->Fill(probeetaID); 
                hphiseleqnnp->Fill(probephiID); } 
		       if (probeq<0 && probeetaID <0) { 
                hptmsseleqnnn->Fill(probept/1000.); 
                hetaseleqnnn->Fill(probeetaID); 
                hphiseleqnnn->Fill(probephiID); } 
			       hmapsele->Fill(probeetaID,probephiID);
			     }
			 }
		     }// end HasProbequality
		 } // end IFisTag
	     } // if is Jpsi
	 } // iend of loop on jentry
     } // iend of loop on jentry

   hmapselall->Sumw2();
   hmapall->Sumw2();
   hmapselalleff->Divide(hmapselall,hmapall,1.,1.,"b") ;
   
   if (DEBUG) log << " ntotal Jpsi = " << ntotal<< endl;
   if (DEBUG) log << " ntag = " << ntag<< endl;
   if (DEBUG) log << " nprobe = " << nprobe<< endl;
   if (DEBUG) log << " nprobetrigger = " << nprobetrigger<< endl;
   
   stacoCombComb_jpsiInvmass_os_am_aa_cut3->Write();   
   hdr->Write();
   hdr2->Write(); 
   hptmsall->Write();
   hptmsallqnp->Write();
   hptmsallqnn->Write();
   hptmsallqpnp->Write();
   hptmsallqpnn->Write();
   hptmsallqnnp->Write();
   hptmsallqnnn->Write();
   hptmsselall->Write();
   hptmsselallqnp->Write();
   hptmsselallqnn->Write();
   hptmsselallqpnp->Write();
   hptmsselallqpnn->Write();
   hptmsselallqnnp->Write();
   hptmsselallqnnn->Write();
   hmapselalleff->Write();
   hetaall->Write();
   hetaallqnp->Write();
   hetaallqnn->Write();
   hetaallqpnp->Write();
   hetaallqpnn->Write();
   hetaallqnnp->Write();
   hetaallqnnn->Write();
   hetaselall->Write();
   hetaselallqnp->Write();
   hetaselallqnn->Write();
   hetaselallqpnp->Write();
   hetaselallqpnn->Write();
   hetaselallqnnp->Write();
   hetaselallqnnn->Write();
   hphiall->Write();
   hphiallqnp->Write();
   hphiallqnn->Write();
   hphiallqpnp->Write();
   hphiallqpnn->Write();
   hphiallqnnp->Write();
   hphiallqnnn->Write();
   hphiselall->Write();
   hphiselallqnp->Write();
   hphiselallqnn->Write();
   hphiselallqpnp->Write();
   hphiselallqpnn->Write();
   hphiselallqnnp->Write();
   hphiselallqnnn->Write();
   hmapall->Write();
   hmapselall->Write();
   hptall->Write();
   hptselall->Write();
   
   hptmsb->Write();
   hptmsbqp->Write();
   hptmsbqn->Write();
   hptmsselbqp->Write();
   hptmsselbqn->Write();
   hptmsbqnp->Write();
   hptmsbqnn->Write();
   hptmsbqpnp->Write();
   hptmsbqpnn->Write();
   hptmsbqnnp->Write();
   hptmsbqnnn->Write();
   hptmsselb->Write();
   hptmsselbqnp->Write();
   hptmsselbqnn->Write();
   hptmsselbqpnp->Write();
   hptmsselbqpnn->Write();
   hptmsselbqnnp->Write();
   hptmsselbqnnn->Write();
   hetab->Write();
   hetabqp->Write();
   hetabqn->Write();
   hetaselbqp->Write();
   hetaselbqn->Write();
   hetabqnp->Write();
   hetabqnn->Write();
   hetabqpnp->Write();
   hetabqpnn->Write();
   hetabqnnp->Write();
   hetabqnnn->Write();
   hetaselb->Write();
   hetaselbqnp->Write();
   hetaselbqnn->Write();
   hetaselbqpnp->Write();
   hetaselbqpnn->Write();
   hetaselbqnnp->Write();
   hetaselbqnnn->Write();
   hphib->Write();
   hphibqp->Write();
   hphibqn->Write();
   hphiselbqp->Write();
   hphiselbqn->Write();
   hphibqnp->Write();
   hphibqnn->Write();
   hphibqpnp->Write();
   hphibqpnn->Write();
   hphibqnnp->Write();
   hphibqnnn->Write();
   hphiselb->Write();
   hphiselbqnp->Write();
   hphiselbqnn->Write();
   hphiselbqpnp->Write();
   hphiselbqpnn->Write();
   hphiselbqnnp->Write();
   hphiselbqnnn->Write();
   hmapb->Write();
   hmapselb->Write();
   
   hptmse->Write();
   hptmseqnp->Write();
   hptmseqnn->Write();
   hptmseqpnp->Write();
   hptmseqpnn->Write();
   hptmseqnnp->Write();
   hptmseqnnn->Write();
   hptmssele->Write();
   hptmsseleqnp->Write();
   hptmsseleqnn->Write();
   hptmsseleqpnp->Write();
   hptmsseleqpnn->Write();
   hptmsseleqnnp->Write();
   hptmsseleqnnn->Write();
   hetae->Write();
   hetaeqnp->Write();
   hetaeqnn->Write();
   hetaeqpnp->Write();
   hetaeqpnn->Write();
   hetaeqnnp->Write();
   hetaeqnnn->Write();
   hetasele->Write();
   hetaseleqnp->Write();
   hetaseleqnn->Write();
   hetaseleqpnp->Write();
   hetaseleqpnn->Write();
   hetaseleqnnp->Write();
   hetaseleqnnn->Write();
   hphie->Write();
   hphieqnp->Write();
   hphieqnn->Write();
   hphieqpnp->Write();
   hphieqpnn->Write();
   hphieqnnp->Write();
   hphieqnnn->Write();
   hphisele->Write();
   hphiseleqnp->Write();
   hphiseleqnn->Write();
   hphiseleqpnp->Write();
   hphiseleqpnn->Write();
   hphiseleqnnp->Write();
   hphiseleqnnn->Write();
   hmape->Write();
   hmapsele->Write();
       
   TagInvMass->Write();
   ProbeInvMass->Write();
   my_h1->Write();
}

bool TP::MuonQuality(int index)
{
  bool mupass= 0;
  float mueta  ;
  float muphi  ; 
  float mupt   ;
  float muz0  ;
  if (index == 1)
    {
       mueta  = mu1eta;
       muphi  = mu1phi; 
       mupt   = mu1pt;
       muz0   = mu1z0;
    }
  else
    {
      mueta  = mu2eta;
      muphi  = mu2phi; 
      mupt   = mu2pt;
      muz0   = mu2z0;
    }

  float theta=2*atan(exp(-mueta));
  float mup=mupt/(sin(theta));
  
  /*
  if ( fabs(pvz) < 150 && fabs(mueta)<2.4 
       && mupt >1500 && mup>1500 
       && fabs(muz0-pvz) < 3 && msinfo == 1)
  */

  mupass=1;
  
  return mupass;
}

float TP::SearchDrCut(float pt, float muetas, float muphis, float t[7][32][42])
{

  //if (DEBUG) log << "test"<<endl;
  float dr_cut;
  int ptzone[100];
  //if (DEBUG) log << "Assign ptzone n2 "<<endl;
   
   for (int p=0; p<100; p++)
     {
	 if ( p< 4) ptzone[p]=0;
	 if ( p==4) ptzone[p]=1;
	 if ( p==5) ptzone[p]=2;
	 if ( p == 6 || p == 7) ptzone[p]=3;
	 if ( p == 8 || p == 9) ptzone[p]=3;
	 if ( p == 10 || p == 11) ptzone[p]=3;
	 if ( p == 12 || p == 13 || p== 14) ptzone[p]=3;
	 if ( p >= 15 ) ptzone[p]=3;
       }     
  
  float rpt = pt/1000.;
  int  i=(int)rpt;
  if (i >= 100) i = 99;
  int zone = ptzone[i]-1;  
  

  float infeta = -3.;
  float passoeta = 0.142; 
  float getetabin=fabs(muetas-infeta)/passoeta;
  Int_t ebin=((int)getetabin);
  
  float infphi = -3.14;
  float passophi = 0.196; 
  float getphibin=fabs(muphis-infphi)/passophi;
  Int_t pbin=((int)getphibin);
  

  if (zone < 0) dr_cut = 0.5; // flat 
  if (ebin > 41) ebin = 41;
  if (pbin > 31) pbin = 31;
  
  // if (DEBUG) log << " pt = "<<i<<" and zone = "<<zone<<endl;
  //if (DEBUG) log << " pbin ="<<pbin<<" and etabin= "<<ebin<<endl;
  
      
  if (zone >= 0) dr_cut = t[zone][pbin][ebin];
  if (dr_cut == 999) dr_cut = 0.5; // put flat dr cut when no entries in map
  
  
   
  
  return dr_cut;
}