#define lambda_selector_cxx
// The class definition in lambda_selector.h has been generated automatically
// by the ROOT utility TTree::MakeSelector().
//
// This class is derived from the ROOT class TSelector.
// The following members functions are called by the TTree::Process() functions:
//    Begin():       called everytime a loop on the tree starts,
//                   a convenient place to create your histograms.
//    Notify():      this function is called at the first entry of a new Tree
//                   in a chain.
//    ProcessCut():  called at the beginning of each entry to return a flag,
//                   true if the entry must be analyzed.
//    ProcessFill(): called in the entry loop for all entries accepted
//                   by Select.
//    Terminate():   called at the end of a loop on the tree,
//                   a convenient place to draw/fit your histograms.
//
//   To use this file, try the following session on your Tree T
//
// Root > T.Process("lambda_selector.C")
// Root > T.Process("lambda_selector.C,"some options"")
// Root > T.Process("lambda_selector.C+")
//
#include "lambda_selector.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "iostream.h"
#include "TCut.h"
float Lambda_mass_min=1;
float Lambda_mass_max=1.5;
float Kshort_mass_min=0.3;
float Kshort_mass_max=0.9;
int i,j,k;

TH1F* hbestlambda_naive=new TH1F("naive lambda mass truth id","naive lambda mass truth id",100,1.08,1.2); 
TH1F* hbestlambda_corr=new TH1F("naive lambda massn truth id","naive lambda massn truth id",100,1.08,1.2);
TH1F* hbestlambda_tvx=new TH1F("lambda mass with true vertex","lambda mass with true vertex",100,1.08,1.2);

TH1F* hKshort_mass_background=new TH1F("Kshort mass background","Kshort mass background",100,Kshort_mass_min,Kshort_mass_max);
TH1F* hKshort_mass=new TH1F("Kshort mass from Kshort","Kshort mass from Kshort",100,Kshort_mass_min,Kshort_mass_max);
TH1F* hKshort_mass_lambda=new TH1F("Kshort mass from lambda","Kshort mass from lambda",100,Kshort_mass_min,Kshort_mass_max);

TH1F* hresolution_primvx=new TH1F("vx resolution primary","vx resolution primary",100,-0.1,0.1);
TH1F* hresolution_secondary1vx=new TH1F("vx resolution secondary1","vx resolution secondary1",100,-5,5);
TH1F* hresolution_secondary2vx=new TH1F("vx resolution secondary2","vx resolution secondary2",100,-10,10);
TH1F* hresolution_primvy=new TH1F("vy resolution primary","vy resolution secondary2",100,-0.1,0.1);
TH1F* hresolution_secondary1vy=new TH1F("vy resolution secondary1","vy resolution secondary1",100,-5,5);
TH1F* hresolution_secondary2vy=new TH1F("vy resolution secondary2","vy resolution secondary2",100,-10,10);
TH1F* hresolution_dvz_lambda=new TH1F("resolution in z for lambdas","resolution in z for lambdas",100,-200,200);
TH1F* hresolution_dvz=new TH1F("resolution in z for all","resolution in z for all",100,-200,200);

TH1F* hLxy_resolution=new TH1F("resolution of Lxy=0","resolution of Lxy=0",100,-10,10);
TH1F* hhhhLxy_res=new TH1F("resolution of Lxy>0","resolution of Lxy>0",100,-10,10);
TH1F* hLxy_true_all=new TH1F("true Lxy background","true Lxy background",100,0,100);
TH1F* hLxy_true_lambda=new TH1F("true Lxy lambda","true Lxy lambda",100,0,100);
TH1F* hLxy_recon_all=new TH1F("recon Lxy background","recon Lxy background",100,0,100);
TH1F* hLxy_recon_lambda=new TH1F("recon Lxy lambda","recon Lxy lambda",100,0,100);

TH2F* hdelta_Lxy_Lxy_true=new TH2F("Lxy vs delta Lxy","delta Lxy vs Lxy",50,-50,50,50,0,150);
TH2F* hdelta_Lxy_pt_recon=new TH2F("delta Lxy vs lowset pt","delta Lxy vs lowest pt",50,-50,50,100,0,10);
TH2F* hSvxhits_Lxy=new TH2F("Svx hits vs Lxy","Svx hits vs Lxy",50,0,10,50,0,150);
TH2F* hSvxhits_Lxy_recon=new TH2F("Svx hits vs recon. Lxy","Svx hits vs recon. Lxy",50,0,50,50,0,150);
TH2F* hCothits_Lxy=new TH2F("Cot hits vs Lxy","Cot hits vs Lxy",50,0,10,50,0,150);
TH2F* hCothits_Lxy_recon=new TH2F("Cot hits vs recon. Lxy","Cot hits vs recon. Lxy",50,0,50,50,0,150);




//________________________________________________________________________________________
               float cut_dvz=70;           //difference in z coordinate of vertex       //|
               float cut_Lxytrue=0;          //difference in z coordinate in truth      //|
               float cut_Lxy=100;                                                       //|
               float cut_D0_pion=0.1;                                                   //|
               float cut_D0_proton=0;                                                   //|
	       float cut_maxpt_pion=1;                                                  //|
	       float cut_minpt_proton=0.9;                                              //|
	       float cut_maxenergy_pion=1.2;                                            //|
	       float cut_minenergy_proton=1;                                            //|
	       //Mass windows                                                           //|
	       float lambda_mass_window_max =1.15;                                      //|
	       float lambda_mass_window_min =1.1;                                       //|
	       float jpsi_mass_window_max=3.12;                                         //|
	       float jpsi_mass_window_min=3.05;                                         //|
               float Kshort_mass_window_min=0.47;                                       //|
               float  Kshort_mass_window_max=0.53;                                      //|
//______________________________________________________________________________________//|


//____________________________________UNIVERSAL CUTS_____________________________________


//________________________________________________________________________________________
void lambda_selector::Begin(TTree *tree)
{
  Init(tree);


TString mass_cut[6]={"no cuts","cut on D0","cut on pt","cut on Lxy","cut on dvz","cut on Kshort mass"};
TString mass_type[3]={ "naive"," corrected"," true vertex"};
for(int i=0;i<6;i++){for(int j=0;j<3;j++){
    hlambda_mass[i][j] =new TH1F(mass_cut[i]+mass_type[j],mass_cut[i]+mass_type[j],100,Lambda_mass_min,Lambda_mass_max); 
}}


  strcpy(Option,GetOption());
  i=0;  while(Option[i]!='q'){
if(Option[i]=='1')Option_int[i]=1;
i++;
  }

    // Take info from tree here
TFile *f2 =new TFile("Lambda_analysisb.root","READ");
TTree *lambda_tree2 =(TTree*)f2.Get("Lambda_analysis");
number_of_entries=lambda_tree2->GetEntries();
f2->Close();   
}
//________________________________________________________________________________________
Bool_t lambda_selector::ProcessCut(Int_t entry)
{
  if(entry==number_of_entries-1){return kTRUE;}

  if(Option_int[0]==2){
    i=0;  while(Option[i]!='q'){
      Option_int[i]=1; i++;}}

  if(Option_int[1] || Option_int[2] || Option_int[3])return kTRUE;
  
   if(Option_int[4]){
    b_track_combinations->GetEntry(entry); 
    if(ffrom_lambda && fLxy!=-99 && fIs_truly_good)return kTRUE;}
  
  if(Option_int[5] || Option_int[6]){
      b_track_combinations->GetEntry(entry); 
      if(ffrom_lambda && fLxy!=-99 && fvalidity && fIs_truly_good)return kTRUE;}
  
  if(Option_int[7]){permit=1; return kTRUE;}
  
  if(Option_int[8]){
    b_track_combinations->GetEntry(entry);
    if(fvalidity && ftrack_GenpNumberi!=ftrack_GenpNumberj && fIs_truly_good)	return kTRUE;}
  
  if(Option_int[9]){b_track_combinations->GetEntry(entry);
  if(ftrack_GenpNumberi!=ftrack_GenpNumberj) return kTRUE;}
    
    if(Option_int[10]){
      b_track_combinations->GetEntry(entry);  
      if(fLxy!=-99 && fLxy!=0 && fvalidity && fIs_truly_good)     return kTRUE; }

    if(Option_int[11]){ b_track_combinations->GetEntry(entry);
    if(fvalidity && fIs_truly_good && fLxy>0) return kTRUE;}
    
    return kFALSE;
}
//________________________________________________________________________________________
void lambda_selector::ProcessFill(Int_t entry)
{
  if(Option_int[0]==2){
 strcpy(Option,GetOption());
  i=0;  while(Option[i]!='q'){
if(Option[i]=='1')Option_int[i]=1;
i++;
  }}

    if(Option_int[1] || Option_int[2] ||Option_int[3] || permit){
  bool break_flag[6];
  bool survivor0[3],survivor1[3];
    b_track_combinations->GetEntry(entry);      

    for(int i =1;i<4;i++){
      if(Option_int[i] || permit){
	if(Option_int[1]){lambda_mass0=flambda_mass0;lambda_mass1=flambda_mass1;}
	if(Option_int[2]){lambda_mass0=flambda_mass0_corr;lambda_mass1=flambda_mass1_corr;}
	if(Option_int[3]){lambda_mass0=flambda_mass0_tvx;lambda_mass1=flambda_mass1_tvx;}
	j=i-1;
	if(lambda_mass0>0){
 hlambda_mass[0][j]->Fill(lambda_mass0);
 hlambda_mass[0][j]->Fill(lambda_mass1);
 cut[0][j]++;
 if(ffrom_lambda){lambda_cut[0][j]++;}

  //proton pion
 if(fabs(fD0i)>cut_D0_proton && fabs(fD0j)>cut_D0_pion){
   hlambda_mass[1][j]->Fill(lambda_mass0);
 if(Option_int[0]){cut[1][j]++;if(ffrom_lambda)lambda_cut[1][j]++; break_flag[1]=true; }
  
 if(fPti>cut_minpt_proton && fPtj<cut_maxpt_pion){
 hlambda_mass[2][j]->Fill(lambda_mass0);
 if(Option_int[0]){cut[2][j]++;if(ffrom_lambda)lambda_cut[2][j]++; break_flag[2]=true; }
  
 if(fvalidity && fLxy_recon>cut_Lxy){
 hlambda_mass[3][j]->Fill(lambda_mass0);
 if(Option_int[0]){cut[3][j]++;if(ffrom_lambda)lambda_cut[3][j]++; break_flag[3]=true; }
  
 if(abs(fdvz)<cut_dvz){
 hlambda_mass[4][j]->Fill(lambda_mass0);
 if(Option_int[7] && i==1){ 
   if(ffrom_Kshort && fIs_truly_good){
     hKshort_mass->Fill(fKshort_mass_corr); 
     surviving_Kshort++;
}
   if(ffrom_lambda && fIs_truly_good){hKshort_mass_lambda->Fill(fKshort_mass_corr);}}
 
 if(Option_int[0]){cut[4][j]++;if(ffrom_lambda)lambda_cut[4][j]++; break_flag[4]=true; }
 if(fKshort_mass_corr>Kshort_mass_window_min && fKshort_mass_corr<Kshort_mass_window_max){
 hlambda_mass[5][j]->Fill(lambda_mass0);
 if(Option_int[0]){cut[5][j]++;if(ffrom_lambda)lambda_cut[5][j]++;}
 survivor0[j]=true;
 }}}}}}
 //other combination
	if(lambda_mass1>0){
 if(fabs(fD0j)>cut_D0_proton && fabs(fD0i)>cut_D0_pion){
 hlambda_mass[1][j]->Fill(lambda_mass1);
 if(Option_int[0]){cut[1][j]++;if(ffrom_lambda && break_flag[1]==false)lambda_cut[1][j]++;}
 if(break_flag[1])break_flag[1]=false;

 if(fPtj>cut_minpt_proton && fPti<cut_maxpt_pion){
 hlambda_mass[2][j]->Fill(lambda_mass1);
 if(Option_int[0]){cut[2][j]++;if(ffrom_lambda && break_flag[2]==false)lambda_cut[2][j]++;}
 if(break_flag[2])break_flag[2]=false;

 if(fvalidity && fLxy_recon>cut_Lxy){
 hlambda_mass[3][j]->Fill(lambda_mass1);
 if(Option_int[0]){cut[3][j]++;if(ffrom_lambda && break_flag[3]==false)lambda_cut[3][j]++;}
 if(break_flag[3])break_flag[3]=false;

 if(abs(fdvz)<cut_dvz){ 
 hlambda_mass[4][j]->Fill(lambda_mass1);
 if(Option_int[7] && i==1){ 
 if(ffrom_Kshort && fIs_truly_good){
   hKshort_mass->Fill(fKshort_mass_corr); 
   surviving_Kshort++;}
 if(ffrom_lambda && fIs_truly_good){
   hKshort_mass_lambda->Fill(fKshort_mass_corr);}}

 if(Option_int[0]){cut[4][j]++;if(ffrom_lambda && break_flag[4]==false)lambda_cut[4][j]++;}
 if(break_flag[4])break_flag[4]=false;

 if(fKshort_mass_corr>Kshort_mass_window_min && fKshort_mass_corr<Kshort_mass_window_max){
   hlambda_mass[5][j]->Fill(lambda_mass1);
 if(Option_int[0]){cut[5][j]++;if(ffrom_lambda && break_flag[5]==false)lambda_cut[5][j]++;}
 if(break_flag[5])break_flag[5]=false;
 survivor1[j]=true;
 }}}}}}}
    
      //look at background contributions after surviving
      if(survivor0[j]==true || survivor1[j]==true){
    if(abs(ftrack_idi)==211 && abs(ftrack_idj)==211)twopi_count[j]++;
    if(abs(ftrack_idi)==2212 && abs(ftrack_idj)==2212)twoproton_count[j]++;
    if(abs(ftrack_idi)==321 && abs(ftrack_idj)==321)twoproton_count[j]++;
    if(ftrack_idi==-1 && ftrack_idj==-1)twounknown_count[j]++;
    if((abs(ftrack_idi)==211 && abs(ftrack_idj)==2212) ||
       (abs(ftrack_idj)==211 && abs(ftrack_idi)==2212))piproton_count[j]++;
    //  if((abs(ftrack_idi)==211 && abs(ftrack_idj)==321) ||
    //   (abs(ftrack_idj)==211 && abs(ftrack_idi)==321)) pikaon_count[j]++;
     if((abs(ftrack_idi)==211 && ftrack_idj==-1) ||
         (abs(ftrack_idj)==211 && ftrack_idi==-1))piunknown_count[j]++;
     //      if((abs(ftrack_idi)==2212 && abs(ftrack_idj)==321) ||
     //  (abs(ftrack_idj)==2212 && abs(ftrack_idi)==321))protonkaon_count[j]++;
      if((abs(ftrack_idi)==2212 && ftrack_idj==-1) ||
         (abs(ftrack_idj)==2212 && ftrack_idi==-1))protonunknown_count[j]++;
     if((abs(ftrack_idi)==321 && ftrack_idj==-1) ||
         (abs(ftrack_idj)==321 && ftrack_idi==-1))kaonunknown_count[j]++;
    if(Option_int[i]){survivor0[k]=false;survivor1[k]=false;}  
      } 
    if(permit)break;
    }}




if(Option_int[4]){
    if(abs(ftrack_idi)==2212)hbestlambda_naive->Fill(flambda_mass0);
    if(abs(ftrack_idi)==211)hbestlambda_naive->Fill(flambda_mass1);}
  
if(Option_int[5]){
    if(abs(ftrack_idi)==2212)hbestlambda_corr->Fill(flambda_mass0_corr);
    if(abs(ftrack_idi)==211)hbestlambda_corr->Fill(flambda_mass1_corr);}

if(Option_int[6]){
    if(abs(ftrack_idi)==2212)hbestlambda_tvx->Fill(flambda_mass0_tvx);
    if(abs(ftrack_idi)==211)hbestlambda_tvx->Fill(flambda_mass1_tvx);}

 if(Option_int[7]){
if(entry ==number_of_entries-1)permit=0;
if(fIs_truly_good && ffrom_Kshort){
  max_reconstructed_Kshorts++;
   hKshort_mass_background->Fill(fKshort_mass_corr);}}

 if(Option_int[8]){
    b_track_combinations->GetEntry(entry);      
if(fLxy==0)hresolution_primvx->Fill(fvx);
if(fLxy==0)hresolution_primvy->Fill(fvy);
if(fLxy>0 && fLxy<5)hresolution_secondary1vx->Fill(fvx-fvx_truei);
if(fLxy>0 && fLxy<5)hresolution_secondary1vy->Fill(fvy-fvy_truei);
if(fLxy>50)hresolution_secondary2vx->Fill(fvx-fvx_truei);
if(fLxy>50)hresolution_secondary2vy->Fill(fvy-fvy_truei);
if(fLxy_recon>100 && ffrom_lambda)hresolution_dvz_lambda->Fill(fdvz);
if(fLxy_recon>100 && fLxy!=-99 && !ffrom_lambda && fIs_truly_good)hresolution_dvz->Fill(fdvz);
}
 
 if(Option_int[9]){
   b_track_combinations->GetEntry(entry);
if(fLxy==0 && fvalidity)hLxy_resolution->Fill(fLxy-fLxy_recon);
if(fLxy>0 && fvalidity)hhhhLxy_res->Fill(fLxy-fLxy_recon);
if(fvx_truei!=-99 && fvx_truej!=-99 && fLxy>0)hLxy_true_all->Fill(fLxy);
//if(ffrom_lambda)hLxy_true_lambda->Fill(fLxy);
//if(fvalidity)hLxy_recon_all->Fill(fLxy_recon);
 if(fvalidity && ffrom_lambda)hLxy_recon_lambda->Fill(fLxy_recon);}
 
 if(Option_int[10]){
   b_track_combinations->GetEntry(entry);      
   //   hdelta_Lxy_Lxy_true->Fill(fLxy-fLxy_recon,fLxy);
   if(fPti<fPtj)lowest_pt=fPti;
   else lowest_pt=fPtj;
   // hdelta_Lxy_pt_recon->Fill(fLxy-fLxy_recon,lowest_pt);
 }
 
    if(Option_int[11]){
      int martin= -9;
      //      b_track_combinations->GetEntry(entry);      
     //     hSvxhits_Lxy->Fill(fAxial_SVX_hitsi,fLxy);   
     //  hSvxhits_Lxy->Fill(fAxial_SVX_hitsj,fLxy);
     //     hSvxhits_Lxy_recon->Fill(fAxial_SVX_hitsi,fLxy_recon);
//     hSvxhits_Lxy_recon->Fill(fAxial_SVX_hitsj,fLxy_recon);
     //    hCothits_Lxy->Fill(fAxial_COT_hitsi,fLxy);   
//     hCothits_Lxy->Fill(fAxial_COT_hitsj,fLxy);
//     hCothits_Lxy_recon->Fill(fAxial_COT_hitsi,fLxy_recon);
//     hCothits_Lxy_recon->Fill(fAxial_COT_hitsj,fLxy_recon);
 }
 //__________________________FILL ADDITIONAL COUNTERS HERE____________________________
 if(Option_int[0]){
   if(fvalidity && fLxy!=-99)correct_goodvalidity_count++;
   if(!fvalidity && fLxy!=-99)correct_badvalidity_count++;
   if(fvalidity && fLxy!=-99)wrong_validity_count0++;
   if(!fvalidity && fLxy==-99)wrong_validity_count1++;
   if(fLxy>0)secondary_vertex_count++;
 }


}
//_______________________________________________________________________________________
void lambda_selector::Terminate()
{
  float purity_mass[6][3];
  float efficiency_mass[6][3];

  //lambda mass search
if(Option_int[1] || Option_int[2] ||Option_int[3]){
if(Option_int[1]){clambda_mass[0] =new TCanvas("lambda_mass","lambda_mass",0,0,900,500); 
clambda_mass[0]->Divide(3,2);}
if(Option_int[2]){clambda_mass[1] =new TCanvas("lambda_mass corr","lambda_mass corr",0,0,900,500);
clambda_mass[1]->Divide(3,2);}
if(Option_int[3]){clambda_mass[2] =
new TCanvas("lambda_masses with true vertex","lambda masses with true vertex",0,0,900,500);
clambda_mass[2]->Divide(3,2);}
 
    for(int i=1;i<4;i++){
      if(Option_int[i]){
	j=i-1;
 clambda_mass[j]->cd(1);   hlambda_mass[0][j]->Draw();
 clambda_mass[j]->cd(2);   hlambda_mass[1][j]->Draw();
 clambda_mass[j]->cd(3);   hlambda_mass[2][j]->Draw();
 clambda_mass[j]->cd(4);   hlambda_mass[3][j]->Draw();
 clambda_mass[j]->cd(5);   hlambda_mass[4][j]->Draw();
 clambda_mass[j]->cd(6);   hlambda_mass[5][j]->Draw();
      }}}

  //best lambda masses
 if(Option_int[4]){
 TCanvas *cbestlambda_naive =new TCanvas("naive lambda mass true id","naive lambda_mass true id",0,0,400,400);
 hbestlambda_naive->Draw();}

 if(Option_int[5]){TCanvas *cbestlambda_corr =new TCanvas("corr lambda mass true id","corr lambda_mass true id",0,0,400,400);
 hbestlambda_corr->Draw();}

 if(Option_int[6]){TCanvas *cbestlambda_truthvertex =new TCanvas("lambda_mass with true vertex","lambda_mass with true vertex",0,0,400,400);
 hbestlambda_tvx->Draw();}

 if(Option_int[7]){
TCanvas *cKshort_mass =new TCanvas("Kshort masses","Kshort masses",0,0,900,500);
cKshort_mass->Divide(2,1);
cKshort_mass->cd(1);
hKshort_mass_background->Draw();
 cKshort_mass->cd(2);
hKshort_mass_lambda->Draw();
hKshort_mass->SetLineColor(9);
 hKshort_mass->Draw("same");}
 

 if(Option_int[8]){
TCanvas *cresolution_vertex =new TCanvas("vertex resolution","vertex resolution",0,0,500,1000);
cresolution_vertex->Divide(2,4);
cresolution_vertex->cd(1);hresolution_primvx->Draw();
cresolution_vertex->cd(2);hresolution_primvy->Draw();
cresolution_vertex->cd(3);hresolution_secondary1vx->Draw();
cresolution_vertex->cd(4);hresolution_secondary1vy->Draw();
cresolution_vertex->cd(5);hresolution_secondary2vx->Draw();
cresolution_vertex->cd(6);hresolution_secondary2vy->Draw();
cresolution_vertex->cd(7);hresolution_dvz_lambda->Draw();
cresolution_vertex->cd(8);hresolution_dvz->Draw();}

 if(Option_int[9]){
TCanvas *cLxy =new TCanvas("Lxy plots","Lxy plots",0,0,400,800);
cLxy->Divide(1,3);
cLxy->cd(1);hLxy_resolution->Draw();
//cLxy->cd(2);hLxy_resolution->Draw();
cLxy->cd(2);hhhhLxy_res->Draw();
cLxy->cd(3);hLxy_true_all->Draw();
hLxy_true_lambda->Draw();
hLxy_recon_lambda->SetLineColor(8);
hLxy_recon_lambda->Draw("same"); 

}
 if(Option_int[10]){
TCanvas *cLxy_correlations =new TCanvas("Lxy correlations","Lxy correlations",0,0,400,800);
cLxy_correlations->Divide(1,2);
cLxy_correlations->cd(1);hdelta_Lxy_Lxy_true->Draw();
cLxy_correlations->cd(2);hdelta_Lxy_pt_recon->Draw();}

 if(Option_int[11]){
TCanvas *cSilicon_correlations =new TCanvas("Silicon correlations","Silicon correlations",0,0,400,800);
cSilicon_correlations->Divide(2,2);
cSilicon_correlations->cd(1);hSvxhits_Lxy->Draw();
cSilicon_correlations->cd(2);hSvxhits_Lxy_recon->Draw();
cSilicon_correlations->cd(3);hCothits_Lxy->Draw();
cSilicon_correlations->cd(4);hCothits_Lxy_recon->Draw();
 }


//_______________________________PRINT OUT STATS_______________________________

   if(Option_int[0]){
     cout<<"Total number of track pair combinations ="<<number_of_entries<<endl<<endl;;
     cout<<"vertex validity when vertex exists ="<<correct_goodvalidity_count<<endl;
     cout<<"vertex validity when vertex does notexists ="<<correct_badvalidity_count<<endl;
     cout<<"tags vertex when it does not exist ="<<wrong_validity_count0<<endl;
     cout<<"fails to tags vertex when vertex exists ="<<wrong_validity_count1<<endl;
     cout<<"Number of secondary vertices ="<<secondary_vertex_count<<endl<<endl;
if(Option_int[1] || Option_int[2] ||Option_int[3]){
  TString cut_list[6]={"no cuts","D0 cuts","pt cuts","Lxy cut","dvz cut","Ks cuts"};

  if(Option_int[1])cout<<"\nnaive mass statistics"<<endl;
  if(Option_int[2])cout<<"\ncorrected mass statistics"<<endl;
  if(Option_int[3])cout<<"\ntrue vertex mass statistics"<<endl;

  cout<<"cut type \tlambdas\ttotal\tpurity\t\tefficiency"<<endl;
  cout<<"__________________________________________________________________"<<endl;
    for(int i=1;i<4;i++){
      j=i-1;
      if(Option_int[i]){


  for(k=0;k<6;k++)
{ 
  if(lambda_cut[0][j]==0)lambda_cut[0][j]=0.01;
    purity_mass[k][j]=lambda_cut[k][j]/cut[k][j];
    efficiency_mass[k][j]=lambda_cut[k][j]/lambda_cut[0][j];
    cout<<cut_list[k]<<"\t\t"<<lambda_cut[k][j]<<"\t"<<cut[k][j]
	<<"\t"<<purity_mass[k][j]<<"\t"<<efficiency_mass[k][j]<<endl;
    
}
  cout<<endl;
  cout<<"pi proton ="<<piproton_count[j]<<endl;  
  cout<<"pi pi ="<<twopi_count[j]<<endl;  
  cout<<"proton proton ="<<twoproton_count[j]<<endl;  
  cout<<"unknown unknown ="<<twounknown_count[j]<<endl;  
  cout<<"k k ="<<twok_count[j]<<endl;  
  cout<<"pi k ="<<pik_count[j]<<endl;  
  cout<<"proton k ="<<protonk_count[j]<<endl;  
  cout<<"unknown pion ="<<piunknown_count[j]<<endl;
  cout<<"unknown proton ="<<protonunknown_count[j]<<endl;
  cout<<"unknown kaon ="<<kaonunknown_count[j]<<endl;
      }}}
 if(Option_int[7]){
   cout<<"maximum number of Kshorts ="<<max_reconstructed_Kshorts<<endl;
   cout<<"Number of Kshorts surving lambda mass cuts ="<<surviving_Kshort<<endl;

 }}
   cout<<"End of macro"<<endl;
}