float fabs(float x) {
  if (x>0) return x; else return -x;
}

void effi_fake()
{ 
  Int_t neve = 0;  // number of events to loop over or 0 = all
  float pmin = 2.5; // momentum cut in GeV
  float pifake = 5.;     // (pi as K)/pi v %    FAKE rate
  float dLKH   = 100;    //obmocje deltaLKH

  f->SetBranchStatus("*",1);
  Stat_t nevent = f->GetEntries();
  if ((neve>0)&&(neve<nevent) ) nevent=neve;
  int count=0;

  float lkh[5];
  int ndiv  = 20;
  int nbins = 1000;
  float cut[20]; 
  float x[20],y[1000];

  TH2F *h100 = new TH2F("h100","L4-L3 vs. p, pi",ndiv,0,100,nbins,-dLKH,dLKH);
  h100->SetYTitle("logL(K)-logL([p])");
  h100->SetXTitle("p (GeV/c)");

  TH2F *h102 = new TH2F("h102","L4-L3 vs. p, K ",ndiv,0,100,nbins,-dLKH,dLKH);
  h102->SetYTitle("logL(K)-logL([p])");
  h102->SetXTitle("p (GeV/c)");


  TH1F *h111 = new TH1F("h111","test",nbins,-dLKH,dLKH);
  TH1F *h104 = new TH1F("h104","Cut for fake pi",ndiv,0,100);

  for( int i=0;i<nevent; event->Clear(),++i ) {
    if(i%10 == 0 ) {
      printf("Event:%d\n",i);
    }
    f->GetEntry(i);
 
    // Loops over all tracks. You get a pointer to the current track
    // by dereferencing the iterator.
    for(   GObjArray<RecoTrack*>::const_iterator track = event->tracks.begin();
         track != event->tracks.end();
         ++track ) {
      if ( ((*track))->GetFlr()==0 ) continue;
      lkh[2] = ((*track))->GetLrpi();
      lkh[3] = ((*track))->GetLrk();
      float dlkh = log(lkh[3])-log(lkh[2]);   
      int itip =((*track))->GetFlt();
      float p = fabs( ((*track))->p() );

      switch (itip) {
      case 3:
        h100->Fill(p,dlkh,1);
	break;
      case 4:
        h102->Fill(p,dlkh,1);
	break;
      }

    }
  }


  float bin= (2.*dLKH)/nbins;

  
  for (int i=0;i< nbins ; i++) y[i]=h111->GetBinCenter(i+1);


  for (int i=0;i< ndiv ; i++) {
    x[i]=h104->GetBinCenter(i+1); 
    char tit[100];
    sprintf(tit,"projpi%d",i+1);
    TH1D *hi1= h100->ProjectionY(tit,i+1,i+1);
    float h1sum= hi1->GetSumOfWeights();
    if (h1sum==0) h1sum=1;
    hi1->Scale(1/h1sum);

    sprintf(tit,"projK%d",i+1);
    TH1D *hi2= h102->ProjectionY(tit,i+1,i+1);
    float h2sum= hi2->GetSumOfWeights();
    if (h2sum==0) h2sum=1;
    hi2->Scale(1/h2sum);


    //    printf("%f %f \n",h1sum,h2sum);
    h1sum=0;
    cut[i]=0;
    for (int j=0;j< nbins ; j++) { 
      h1sum+=hi1->GetBinContent(j+1); 
      float hh=1-h1sum- pifake/100.;
      float hh1=hh+fabs(hh);
      if  (hh1<=0) { 
	cut[i]=y[j]; 
	//printf("%d cut[%d]= %f\n",j, i,cut[i]);
	break;
      }
    }         
    h104->SetBinContent(i+1,cut[i]); 


    /// neki s cutom je treba nardit sam ne vem kaj
  }
  
 
  // determine the efficiency
  TH1F *h101 = new TH1F("h101","pi as K",ndiv,0,100);
  TH1F *h103 = new TH1F("h103","K  as K",ndiv,0,100);
  TH1F *h201 = new TH1F("h201","all PI ",ndiv,0,100);
  TH1F *h203 = new TH1F("h203","all K  ",ndiv,0,100);

  for( int i=0;i<nevent; event->Clear(),++i ) {
    if(i%100 == 0 ) {
      printf("Event:%d\n",i);
    }
    f->GetEntry(i);
 
    // Loops over all tracks. You get a pointer to the current track
    // by dereferencing the iterator.
    for(   GObjArray<RecoTrack*>::const_iterator track = event->tracks.begin();
         track != event->tracks.end();
         ++track ) {
      if ( ((*track))->GetFlr()==0 ) continue;
      lkh[0] = ((*track))->GetLre();
      lkh[1] = ((*track))->GetLrmu();
      lkh[2] = ((*track))->GetLrpi();
      lkh[3] = ((*track))->GetLrk();
      lkh[4] = ((*track))->GetLrp();  
      float dlkh = log(lkh[3])-log(lkh[2]);   
      int itip =((*track))->GetFlt();

      float p =  fabs(((*track))->p()); 

      if (p<pmin) continue;
      int k=int(p/5.);
      if (k>19) k=19;
      
      if (itip==3) {
	h201->Fill(p);
	if (dlkh > cut[k] )  h101->Fill(p);
      } else   if (itip==4) {
	h203->Fill(p);
	if (dlkh > cut[k] )  h103->Fill(p);
      }
      count ++;     
    }
  }

  
  float kask  = h103->GetSumOfWeights();
  float kall  = h203->GetSumOfWeights();
  float piask = h101->GetSumOfWeights();
  if ( kall==0 ) kall = 0.0001;
  float kid = kask+piask;
  if ( kid == 0 )  kid = 0.0001;
  float efall = kask/kall;
  float eefall= sqrt((efall*(1-efall))/kall);
  float puri  = kask/ kid;
  h101->Divide(h201); 
  h101->SetMaximum(1);
  h101->SetXTitle("p (GeV/c)");
  h101->SetYTitle("pi as K / all K");
  
  
  h103->Divide(h203);
  h103->SetMaximum(1);
  h103->SetXTitle("p (GeV/c)");
  h103->SetYTitle("logL(K)-logL([p])");
  h103->SetYTitle("K as K / all K");

  
  for(int i=0;i<ndiv;i++) {
    float error=h103->GetBinContent(i+1);
    float div =  h203->GetBinContent(i+1);
    if (div>0) error= sqrt( error * (1-error) / div ); else error=0;
    h103->SetBinError(i+1,error);
  }
  
  for(int i=0;i<ndiv;i++) {
    float error=h101->GetBinContent(i+1);
    float div  =  h201->GetBinContent(i+1);
    if (div>0) error= sqrt( error * (1-error) / div ); else error=0;
    h101->SetBinError(i+1,error);
  }

  char title[200];
  sprintf(title,"p > %2f Gev/c",pmin);  
  printf("KasK %f Kall %f Effi %f \n",kask,kall,efall);
  printf("KasK %f piasK %f puri %f \n",kask,piask,puri);
  char tit[200];
  sprintf(tit,"Efficiency= %3.2f +- %3.2f ",efall,eefall);
  printf("%s\n",tit);
  char tit1[200];
  sprintf(tit1,"Purity= %3.2f +- %3.2f",puri,sqrt(kask)/ kid );
  printf("%s\n",tit1);

  // now draw the results

  TCanvas *c2= new TCanvas("c2","Kaon separation efficiency in RICH",
			   0 ,0,600,600);
  c2->Divide(2,2);
  c2->cd(1);   h102->Draw(); h104->Draw("SAME");
  c2->cd(2);   h100->Draw(); h104->Draw("SAME");
  c2->cd(3);   h103->Draw();
  c2->cd(4);   h101->Draw();
  gStyle->SetOptStat(0);
  TPaveText *pave = new TPaveText(50,0.8,100,0.9);  
  pave->AddText(tit);  
  pave->AddText(tit1); 
  c2->cd(3);  
  pave->Draw();
  c2->Update();  
  return;  
}