#include "StHiCorrEta.h"
#include "TSystem.h"
#include "TClonesArray.h"
#include "TChain.h"
#include <TStopwatch.h>
#include <TFile.h>
#include <TTree.h>
#include <sstream>
#include "TF2.h"
#include "TProfile.h"
#include "TVector2.h"
#include "TString.h"
#include "TROOT.h"
#include "TObjString.h"
#include "../StHiMicroEvent/StHiMicroEvent.h"
#include "StarClassLibrary/StThreeVectorF.hh"


//define binning
const int StHiCorrEta::mNTrigBin=24;
const float StHiCorrEta::mMinPtTrig=3;
const float StHiCorrEta::mMaxPtTrig=15;
const int StHiCorrEta::mNdPhiBin=2*24;
// Phi range needs to be a multiple of 1 or 2 pi
const float StHiCorrEta::mMindPhi=-TMath::Pi();
const float StHiCorrEta::mMaxdPhi=TMath::Pi();

const int StHiCorrEta::mNdEtaBin=2*20; // 2*20 TPC only
const float StHiCorrEta::mMindEta=-2.0; // -2 TPC only
const float StHiCorrEta::mMaxdEta=2.0; // 2 TPC only

const int StHiCorrEta::mNVtxZBin=12;  // Vtx-z binning for mixed events
const float StHiCorrEta::mMinVtxZ=-30;
const float StHiCorrEta::mMaxVtxZ=30;

const float StHiCorrEta::minEtaFtpc=3;
const float StHiCorrEta::maxEtaFtpc=3.5;

//const double StHiCorrEta::mPtAssocThres[mNAssocThres+1]={2,3,4,5,10}; //last value is dummy for histogramming (standard)
const double StHiCorrEta::mPtAssocThres[mNAssocThres+1]={2,2.5,3,3.5,10};


const bool includeFtpc=false;

// Help Functions

int getHighPtCentBin(int centrality)
{

  Int_t centhighpt=-1;

  // analog zu highpt ?
  if (centrality==9) centhighpt=6;
  else if (centrality==8) centhighpt=5;
  else if (centrality==7) centhighpt=4;
  else if (centrality==6) centhighpt=3;
  else if (centrality==5) centhighpt=2;
  else if (centrality==4 || centrality==3) centhighpt=1;
  else if (centrality==2 || centrality==1) centhighpt=0;
  else centhighpt=-1;

  return centhighpt;
}

// ----------------------------------------------------------


StHiCorrEta::StHiCorrEta()
{
}
StHiCorrEta::~StHiCorrEta()
{
}

Float_t StHiCorrEta::dPhi(Float_t phi1, Float_t phi2) {
  static float twopi=2*TMath::Pi();

  Float_t d_phi=phi1-phi2;
  //if (d_phi < mMindPhi) d_phi += twopi - mMindPhi; 
  //if (d_phi > mMindPhi + twopi) d_phi -= twopi - mMindPhi; 
  //if (d_phi > mMaxdPhi) d_phi = mMindPhi + twopi - d_phi;

  //if (d_phi < mMindPhi) d_phi += twopi; 
  //if (d_phi > mMindPhi + twopi) d_phi -= twopi - mMindPhi; 
  //if (d_phi > mMaxdPhi) d_phi -= twopi;

  if (d_phi < mMindPhi) {d_phi = twopi + d_phi;} 
  if (d_phi > mMaxdPhi) {d_phi = d_phi - twopi;}


  // DEBUG:
  //if (fabs(d_phi)>TMath::Pi())
  //cout<<d_phi<<endl;

  return d_phi;
}

void StHiCorrEta::initHistos() {

  mult = new TH1D("mult","multiplicity",1000,0.,1000.);
  dca = new TH1D("dca","global dca",100,0.,3.);
  nhits = new TH1D("nhits","nfit hits",50,0.,50.);
  tphi = new TH1D("phi","phi angle",100,-3.14,3.14);
  eta = new TH1D("etaTrig","etaTrig",100,-5,5);
  etaAss = new TH1D("etaAss","etaAss",100,-5,5);
  ntrig_evt_h = new TH1D("n_trig_evt_h","Triggers per evt",10,0.5,10.5);
  hVertZ= new TH1D("vertexZ","Vertex z position",50,-50,50);
  hCentrality=new TH1D("centrality","Centrality class",10,0,10);

  //TProfile* v2prof = new TProfile("v2prof","v2 vs pt",240,0.,12.,-1000.,1000.);
  //TH1D* v2    = new TH1D("Flow_vPt_Sel2_Har2","v2 vs pt",240,0.,12.);
  
  mYieldPt = new TH1D("YieldPt", "YieldPt",30,0,15);
  mYieldPt->SetXTitle("P_{t} (GeV)");
  // mFlowPt = new TH1D("FlowPt","FlowPt",30,0,15);
//   mFlowPt->SetXTitle("P_{t} (GeV)");
//   mFlowPt->Sumw2();

  /*
  mPairYield2D = new TH2D("PairYield2D", "mHistPairYield2D",
				    mNTrigBin, mMinPtTrig, mMaxPtTrig, 
				    mNAssocBin, mMinPtAssoc, mMaxPtAssoc); 
  mPairYield2D->SetXTitle("P_{t,trig} (GeV)");
  mPairYield2D->SetYTitle("P_{t,assoc} (GeV)");
  */

  //mEtaPhiPtTrig=new TH3D("EtaPhiPtTrig","EtaPhiPtTrig",20,-1,1,48,-TMath::Pi(),TMath::Pi(),mNTrigBin,mMinPtTrig,mMaxPtTrig);
  mEtadPhiPtTrig=new TH3D("EtadPhiPtTrig","EtadPhiPtTrig",20,-1,1,mNdPhiBin,mMindPhi,mMaxdPhi,mNTrigBin,mMinPtTrig,mMaxPtTrig);

  for (Int_t iAssoc=0; iAssoc < mNAssocThres; iAssoc++) {
    Char_t hname[64];
    Char_t htitle[128];
    
    // Signal Histos
    sprintf(hname,"PtTrigdEtadPhi_%d",iAssoc);
    sprintf(htitle,"PtTrigdEtadPhi_%d: P_{t,assoc} > %4.2f GeV",
	    iAssoc,mPtAssocThres[iAssoc]);
    mPtTrigdEtadPhi[iAssoc] = new TH3D(hname,htitle,
				       mNTrigBin, mMinPtTrig, mMaxPtTrig, 
				       mNdEtaBin, mMindEta, mMaxdEta,
				       mNdPhiBin, mMindPhi, mMaxdPhi ); 

    mPtTrigdEtadPhi[iAssoc]->SetXTitle("P_{t,trig} (GeV)");
    mPtTrigdEtadPhi[iAssoc]->SetYTitle("#Delta #eta");
    mPtTrigdEtadPhi[iAssoc]->SetZTitle("#Delta #phi");

      
    // Mixed Histos
    sprintf(hname,"PtTrigdEtadPhiMix_%d",iAssoc);
    sprintf(htitle,"PtTrigdEtadPhiMix_%d: P_{t,assoc} > %4.2f GeV",
	    iAssoc,mPtAssocThres[iAssoc]);
    mPtTrigdEtadPhiMix[iAssoc] = new TH3D(hname,htitle,
					  mNTrigBin, mMinPtTrig, mMaxPtTrig, 
					  mNdEtaBin, mMindEta, mMaxdEta,
					  mNdPhiBin, mMindPhi, mMaxdPhi ); 

    mPtTrigdEtadPhiMix[iAssoc]->SetXTitle("P_{t,trig} (GeV)");
    mPtTrigdEtadPhiMix[iAssoc]->SetYTitle("#Delta #eta");
    mPtTrigdEtadPhiMix[iAssoc]->SetZTitle("#Delta #phi");

  }

  mNTrigPtTrig = new TH1D("NTrigPtTrig", "NTrigPtTrig",
			  mNTrigBin, mMinPtTrig, mMaxPtTrig); 
    
  mNTrigPtTrig->SetXTitle("P_{t,trig} (GeV)");
  
  // Efficiency histos
  mAssocEffPtTrigPtAssoc = new TH2D("AssocEffPtTrigPtAssoc", "AssocEffPtTrigPtAssoc",
				    mNTrigBin, mMinPtTrig, mMaxPtTrig, 
				    mNAssocThres, mPtAssocThres);
  mAssocEffPtTrigPtAssoc->Sumw2();
  mAssocEffPtTrigPtAssoc->SetXTitle("P_{t,trig} (GeV)");
  mAssocEffPtTrigPtAssoc->SetYTitle("P_{t,assoc} (GeV)");


  // mPtTrigPtAssocFlowAssoc = new TH2D("PtTrigPtAssocFlowAssoc", "PtTrigPtAssocFlowAssoc",
// 				      mNTrigBin, mMinPtTrig, mMaxPtTrig,
// 				     mNAssocThres, mPtAssocThres); 
//   mPtTrigPtAssocFlowAssoc->Sumw2();
//   mPtTrigPtAssocFlowAssoc->SetXTitle("P_{t,trig} (GeV)");
//   mPtTrigPtAssocFlowAssoc->SetYTitle("P_{t,assoc} (GeV)");

//   mPtTrigPtAssocFlowTrig = new TH2D("PtTrigPtAssocFlowTrig", "PtTrigPtAssocFlowTrig",
// 				      mNTrigBin, mMinPtTrig, mMaxPtTrig, 
// 				     mNAssocThres, mPtAssocThres); 
//   mPtTrigPtAssocFlowTrig->Sumw2();
//   mPtTrigPtAssocFlowTrig->SetXTitle("P_{t,trig} (GeV)");
//   mPtTrigPtAssocFlowTrig->SetYTitle("P_{t,assoc} (GeV)");

}

//-----------------------------------------------------------------------

void StHiCorrEta::histofill(Int_t maxevents, Int_t cent, Float_t etacut,Int_t mId)
{
//_____________________________

  //Float_t phi_bite=(Float_t)phi_bite_deg/180.*TMath::Pi();

  TF2*               mEfficiencyMap;
  TString            mEfficiencyFileName; 
  TString fName = "fEfficiency";

  /*
  if(cent==99)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.M1311ny.efficiencyfcn.hist.root";
  //if(cent==1)
  //mEfficiencyFileName="effic/P02gd.Highpt_minbias.01311ny.efficiencyfcn.hist.root";
  if(cent==0)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.01311ny.efficiencyfcn.hist.root";
  if(cent==1)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.11311ny.efficiencyfcn.hist.root";
  //if(cent==2)
  //mEfficiencyFileName="effic/P02gd.Highpt_minbias.11311ny.efficiencyfcn.hist.root";
  if(cent==2)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.21311ny.efficiencyfcn.hist.root";
  if(cent==3)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.31311ny.efficiencyfcn.hist.root";
  if(cent==4)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.41311ny.efficiencyfcn.hist.root";
  if(cent==5)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.51311ny.efficiencyfcn.hist.root";
  if(cent==6)
    mEfficiencyFileName="effic/P02gd.Highpt_minbias.61311ny.efficiencyfcn.hist.root";
  */

  /*
  if(cent==99)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.M1211ny.efficiencyfcn.hist.root";
  //if(cent==1)
  //mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.01311ny.efficiencyfcn.hist.root";
  if(cent==0)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.01211ny.efficiencyfcn.hist.root";
  if(cent==1)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.11211ny.efficiencyfcn.hist.root";
  //if(cent==2)
  //mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.11311ny.efficiencyfcn.hist.root";
  if(cent==2)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.21211ny.efficiencyfcn.hist.root";
  if(cent==3)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.31211ny.efficiencyfcn.hist.root";
  if(cent==4)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.41211ny.efficiencyfcn.hist.root";
  if(cent==5)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.51211ny.efficiencyfcn.hist.root";
  if(cent==6)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.61211ny.efficiencyfcn.hist.root";
  */

  
  if(cent==99)
    mEfficiencyFileName="effic/P02gd.HighpT_embedding_minbias.M1211ny.efficiencyfcn.hist.root";
  if(cent==100)
    mEfficiencyFileName="effic/Piplusminus.61211ny.efficiencyfcn.hist.root";
  //if(cent==1)
  //mEfficiencyFileName="effic/Piplusminus.01211ny.efficiencyfcn.hist.root";
  if(cent==0)
    mEfficiencyFileName="effic/Piplusminus.01211ny.efficiencyfcn.hist.root";
  if(cent==1)
    mEfficiencyFileName="effic/Piplusminus.11211ny.efficiencyfcn.hist.root";
  //if(cent==2)
  //mEfficiencyFileName="effic/Piplusminus.11211ny.efficiencyfcn.hist.root";
  if(cent==2)
    mEfficiencyFileName="effic/Piplusminus.21211ny.efficiencyfcn.hist.root";
  if(cent==3)
    mEfficiencyFileName="effic/Piplusminus.31211ny.efficiencyfcn.hist.root";
  if(cent==4)
    mEfficiencyFileName="effic/Piplusminus.41211ny.efficiencyfcn.hist.root";
  if(cent==5)
    mEfficiencyFileName="effic/Piplusminus.51211ny.efficiencyfcn.hist.root";
  if(cent==6)
    mEfficiencyFileName="effic/Piplusminus.61211ny.efficiencyfcn.hist.root";


  TFile* file=0; 
  if(mEfficiencyFileName==""){
    std::cout << "no efficiency file set " << std::endl;
  }
  else{
    file=new TFile(mEfficiencyFileName.Data());
    if(!file ||!file->IsOpen()){
      std::cout << "Cannot find efficiency file : " << mEfficiencyFileName << std::endl;
    }
    else{
        mEfficiencyMap = (TF2*)file->Get(fName.Data());
        delete file;
    }
  }
  if(!mEfficiencyMap){
    std::cout << "Cannot find efficency function : " << fName.Data() << std::endl;
    std::cout << "WARNING-USING DUMMY EFFIENCY" << std::endl;
    mEfficiencyMap = new TF2("fEfficiency","1");

  }
  
  
  std::cout << "\tefficiency file : " << mEfficiencyFileName << std::endl;
  std::cout << "\tefficiency fcn title : " << mEfficiencyMap->GetTitle() << std::endl;

  TStopwatch timer;
  timer.Start();
  Float_t mbytes =0;

  StHiMicroEvent* event = new StHiMicroEvent();

  TChain* pPicoChain = new TChain("StHiMicroTree");
  pPicoChain->SetBranchAddress("StHiMicroEvent", &event);

  for(int ii=0; ii<mInputFileList->GetLast()+1;ii++){
    //cout << "Adding file " << ((TObjString*) mInputFileList->At(ii))->GetString() << endl;
    pPicoChain->Add(((TObjString*) mInputFileList->At(ii))->GetString());
  }

  cout<<" ..."<<endl;
  //cout << "total events: " << pPicoChain->GetEntries() << endl;
  TString mLastFile=((TObjString*) mInputFileList->At(mInputFileList->GetLast()))->GetString();
  // DEBUG:
  //mLastFile="/eliza2/rnc/putschke/himicro_mid/st_physics_5033127_raw_4010003.himicro.root";
  cout<<"Last file in list : "<<mLastFile<<endl;

  cout << "requested number of events to process: " << maxevents << endl;

  // open output file
  TFile *hfile = (TFile*) gROOT->FindObject("histos.root");
  if (hfile) hfile->Close();
  char fileName[80];
  /*
  if(cent==21)
    sprintf(fileName, "root.flow.pttrig_%d-%d_GeV/ana%2d.root", int(pttriglow),int(pttrighigh),cent);
  else
  */
  //sprintf(fileName, "root.flow.corr/ana%2d.root",cent+10);
  //Char_t outDir[256];
  //sprintf(outDir,"deltaEta.eta_cut%02f",(etacut));
  //if (gSystem->OpenDirectory(outDir)==0)
  //gSystem->mkdir(outDir);
  if (cent<99)
    sprintf(fileName, "deltaEtaCent%1d.root",cent);
  else if (cent==99)
    sprintf(fileName, "deltaEtaMinBias.root");
  else if (cent==100)
    sprintf(fileName, "deltaEtaCentral_mId_%1d.root",mId);

  //sprintf(fileName, "root.flow/ana.central.%2d.root", cent+10);
  hfile = new TFile(fileName,"RECREATE","file with histograms");
  initHistos();

  //Float_t pi=TMath::Pi();
  //Float_t twopi=2*pi;

  Int_t nev=0;
  long int countTracks=0;
  int mixed_events_ready[mNVtxZBin];
  const int nmixed_events=50; // 50 events mixed
  Float_t mphi[mNVtxZBin][nmixed_events+1][1000];
  Float_t meta[mNVtxZBin][nmixed_events+1][1000];
  Float_t mpt[mNVtxZBin][nmixed_events+1][1000];
  Int_t n_associated_per_event[mNVtxZBin][nmixed_events+1];
  Int_t cur_mix_event[mNVtxZBin];

  //Int_t nres_events=0;
  //Int_t assoc_inplane=0;
  //Int_t assoc_outplane=0;
  Int_t assoc_simple=0;

  for(Int_t i_vtx_z=0;i_vtx_z < mNVtxZBin; i_vtx_z++){
    for(Int_t imixed_event=0;imixed_event<nmixed_events;imixed_event++){
      for(Int_t n=0;n<1000;n++) {
	mphi[i_vtx_z][imixed_event][n]=0.;
	meta[i_vtx_z][imixed_event][n]=0.;
	mpt[i_vtx_z][imixed_event][n]=0.;
      }
      n_associated_per_event[i_vtx_z][imixed_event]=0;
    }
    cur_mix_event[mNVtxZBin]=0;
    mixed_events_ready[mNVtxZBin]=0;
  }

  // fill correlation histos
  nev=0;
  long int itrig_simple=0;

  TString oldfile;
  TString newfile;

  // Histogram to keep track of number of unique triggers
  //TH1D *PtAssocHist = new TH1D("PtAssocHist","PtAssocHist",mNAssocBin,mMinPtAssoc,mMaxPtAssoc);

  //for(Int_t i=0; i< pPicoChain->GetEntries(); i++){
  for(Int_t i=0; i<maxevents;i++){
    pPicoChain->Clear();
  
    if(nev%10000 == 0) cout << "event " << nev << endl;
    
    if(nev>maxevents)break;
    mbytes += pPicoChain->GetEntry(i)*1.e-6;
    newfile=pPicoChain->GetFile()->GetName();
    if (mLastFile==newfile) break;

    if(newfile != oldfile)
      cout << "analyzing file " <<  pPicoChain->GetFile()->GetName() << endl;
    oldfile=newfile;
    
    nev++;

    // select only events with the trigger particle

    Bool_t goodevent=0;

    // centrality check

    if (cent<99)
      if(getHighPtCentBin(event->Centrality())!=cent) 
	continue;
    
    // check Vertex-z position

    if (fabs(event->VertexZ())>30) continue;

    for (Int_t nrt=0; nrt < event->NTrack(); nrt++) {
      
      StHiMicroTrack* track = (StHiMicroTrack*) event->tracks()->At(nrt);
      
      if(track->FitPts()<20)continue;
      if(fabs(track->DcaGl())>1.0)continue;
      if(fabs(track->EtaPr())>etacut)continue;
      if(track->PtPr()>mMaxPtTrig)continue;
      if(track->PtPr()<mMinPtTrig)continue;
      goodevent=1;
      //cout << " trigger pt=" << track->Pt() << endl;
    }
    
    if(goodevent==0)continue;

    // Fill Vertex Z histogram

    hVertZ->Fill(event->VertexZ());
    hCentrality->Fill(getHighPtCentBin(event->Centrality())); // Flow centrality definition -> HighPt def.

    Int_t vtx_z_bin=Int_t((event->VertexZ()-mMinVtxZ)/(mMaxVtxZ-mMinVtxZ)*mNVtxZBin);
    if (vtx_z_bin < 0)
      vtx_z_bin = 0 ;
    if (vtx_z_bin >= mNVtxZBin)
      vtx_z_bin = mNVtxZBin-1;
    
    Int_t cur_mix_evt=cur_mix_event[vtx_z_bin];
    n_associated_per_event[vtx_z_bin][cur_mix_evt]=0;
    for(Int_t n=0; n<event->NTrack(); n++){
      StHiMicroTrack* tr =(StHiMicroTrack*) event->tracks()->At(n);
      if (fabs(tr->EtaPr())<2)
	  {
	    // DEBUG:
	    //cout<<tr->nSigmaProton()<<" "<<tr->nSigmaPion()<<endl;

	    if(tr->FitPts()<20)continue;
	    if(fabs(tr->DcaGl())>1.0)continue;
	    if(fabs(tr->EtaPr())>etacut)continue;
	    if(tr->PtPr()<mPtAssocThres[0])continue;

	    //cout<<"Out: nPi = "<<(tr->nSigmaPion())<<endl;

	    if (!(tr->nSigmaPion()>0 && tr->nSigmaPion()<3) && mId==1) continue;
	    if (!(tr->nSigmaPion()>-5.0 && tr->nSigmaPion()<-3.0) && mId==2) continue;

	    //cout<<mId<<" nPi = "<<(tr->nSigmaPion())<<endl;
	  }
      //else if (includeFtpc)
	//{
	  //if(fabs(tr->DcaGl())>1.0)continue;
	  //if (fabs(tr->EtaPr())<minEtaFtpc || fabs(tr->EtaPr())>maxEtaFtpc) continue;
	  //if(tr->PtPr()<mPtAssocThres[0])continue;
	//}

      meta[vtx_z_bin][cur_mix_evt][n_associated_per_event[vtx_z_bin][cur_mix_evt]]= tr->EtaPr();
      mphi[vtx_z_bin][cur_mix_evt][n_associated_per_event[vtx_z_bin][cur_mix_evt]]= tr->PhiPr();
      mpt[vtx_z_bin][cur_mix_evt][n_associated_per_event[vtx_z_bin][cur_mix_evt]]= tr->PtPr();
      n_associated_per_event[vtx_z_bin][cur_mix_evt]++;
    }
 
    mult->Fill(event->NUncorrectedPrimaries());
    
    Int_t ntrig=-1;
    Int_t ntrigPerEvent=0;

    for (Int_t nrt=0; nrt < event->NTrack(); nrt++) {
      
      StHiMicroTrack* track = (StHiMicroTrack*) event->tracks()->At(nrt);

      if(track->FitPts()<20)continue;
      if(fabs(track->DcaGl())>1.0)continue;
      if(fabs(track->EtaPr())>etacut)continue;

      Float_t phi_i = track->PhiPr();; 
      //Float_t dpsi= dPhi(phi_i,psi_i); 

      mYieldPt->Fill(track->PtPr());
      //mFlowPt->Fill(track->Pt(),cos(2*dpsi));

      if(track->PtPr()>mMaxPtTrig)continue;
      if(track->PtPr()<mPtAssocThres[0])continue;

      countTracks++;
      tphi->Fill(track->PhiPr());
      dca->Fill(track->DcaGl());
      nhits->Fill(track->FitPts());
      
      //phi_psi->Fill(fmod((double)dPhi,(double)pi)); 
      //eta->Fill(track->EtaPr());
      
      Float_t pt1,eta1,phi1,pt2,eta2,phi2,dphi,mdphi;
      //Float_t dpsi_assoc, dpsi_trig;
		    
      //Int_t is_inplane=-1;

      // find trigger particle
      if(track->PtPr() > mMinPtTrig){
	pt1  = track->PtPr();
	eta1 = track->EtaPr();
	phi1 = track->PhiPr();
	ntrig= nrt;
	ntrigPerEvent++;

	eta->Fill(eta1);

	//if(ntrigPerEvent>1)continue; // one trigger per event
	//dpsi_trig=dpsi;
	itrig_simple++;
	mNTrigPtTrig->Fill(pt1);
      }
      else 
	continue;

      dphi=-99.;

      Int_t n_associated=0;
      //PtAssocHist->Reset();

      // find associated particle
      for (Int_t nrt2=0; nrt2 < event->NTrack(); nrt2++) {
      
	if(nrt2 == ntrig) continue;

	StHiMicroTrack* track2 = (StHiMicroTrack*) event->tracks()->At(nrt2);

	if(track2->FitPts()<20)continue;
	if(fabs(track2->DcaGl())>1.0)continue;
	if(fabs(track2->EtaPr())>etacut)continue;
	if(track2->PtPr()>pt1)continue;
	if(track2->PtPr()<mPtAssocThres[0])continue;
	//if (fabs(track2->nSigmaPion())<2) continue; 
	//if (fabs(track2->nSigmaProton())>2) continue;
	//if (fabs(track2->nSigmaPion())>2 && mId==1) continue;
	//if (fabs(track2->nSigmaProton())>2 && mId==2) continue; 
	if (!(track2->nSigmaPion()>0 && track2->nSigmaPion()<3) && mId==1) continue;
	if (!(track2->nSigmaPion()>-5.0 && track2->nSigmaPion()<-3.0) && mId==2) continue;

	pt2  = track2->PtPr();
	eta2 = track2->EtaPr();
	phi2 = track2->PhiPr();

	etaAss->Fill(eta2);

	n_associated++;
	// particle 1 is the trigger
	phi_i = track2->PhiPr(); 
	
	//if(phi_i<0)phi_i+=twopi; // ???? Genauer !!!
	
	//PtAssocHist->Fill(pt2);

	//dpsi_assoc=dPhi(track2->Phi(),psi_i);
	dphi=dPhi(phi2,phi1);

	Int_t iAssoc=0;
	while (iAssoc < mNAssocThres && pt2 > mPtAssocThres[iAssoc]) {
	  Float_t deta=eta2-eta1;
	  Float_t pt_assoc_hist=mAssocEffPtTrigPtAssoc->GetYaxis()->GetBinCenter(iAssoc+1); // pt-assoc value for filling histos
	  mPtTrigdEtadPhi[iAssoc]->Fill(pt1,deta,dphi,1./mEfficiencyMap->Eval(eta2,pt2));
	  //mPtTrigPtAssocdPhiEventAssoc->Fill(pt1,pt_assoc_hist,dpsi_assoc);
	  //mPtTrigPtAssocFlowAssoc->Fill(pt1,pt_assoc_hist,cos(2*dpsi_assoc));
	  //mPtTrigPtAssocdPhiEventTrig->Fill(pt1,pt_assoc_hist,dpsi_trig);
	  //mPtTrigPtAssocFlowTrig->Fill(pt1,pt_assoc_hist,cos(2*dpsi_trig));
	  mAssocEffPtTrigPtAssoc->Fill(pt1,pt_assoc_hist,1./mEfficiencyMap->Eval(eta2,pt2));
	  iAssoc++;
	  assoc_simple++; // ????
	}  // loop over associated thresholds
      }

      // eff. corr in mixed events !?

      //if(mixed_events_are_ready && ntrigPerEvent==1){
      if (mixed_events_ready[vtx_z_bin]) {
	for(Int_t imixed_event=0;imixed_event<nmixed_events+1;imixed_event++){
	  if (imixed_event == cur_mix_evt)
	    continue;
	  for(Int_t n=0;n<n_associated_per_event[vtx_z_bin][imixed_event];n++){
	    mdphi=dPhi(mphi[vtx_z_bin][imixed_event][n],phi1);
	    if (mpt[vtx_z_bin][imixed_event][n] < pt1) {
	      Float_t deta=meta[vtx_z_bin][imixed_event][n]-eta1;
	      Int_t iAssoc=0;
	      while (iAssoc < mNAssocThres && mpt[vtx_z_bin][imixed_event][n] > mPtAssocThres[iAssoc]) {
		mPtTrigdEtadPhiMix[iAssoc]->Fill(pt1,deta,mdphi,1./mEfficiencyMap->Eval(meta[vtx_z_bin][imixed_event][n],mpt[vtx_z_bin][imixed_event][n]));
		iAssoc++;
	      }
	    }
	  }
	}
      }

    }
    ntrig_evt_h->Fill(ntrigPerEvent);
    if ( ntrigPerEvent ) 
      cur_mix_event[vtx_z_bin]++;
    if (cur_mix_event[vtx_z_bin] == nmixed_events) {
      cur_mix_event[vtx_z_bin]=0;
      mixed_events_ready[vtx_z_bin]=1;
    }
  }
  //delete PtAssocHist;

  hfile->Write();
 
  cout << "Pt trig = [ " << mMinPtTrig << "; " << mMaxPtTrig << " ] GeV/c" << endl;
  cout << "Pt assoc therholds: ";
  for (Int_t iAssoc=0; iAssoc < mNAssocThres; iAssoc++) 
    cout << mPtAssocThres[iAssoc] << ", ";
  cout << endl;
  //cout << "Phi bite = " << phi_bite << endl;

  cout << "Total events scanned = " << nev << endl;
  cout << "Total tracks = " << countTracks << endl;

  cout << "    Total triggers=" << itrig_simple << endl;

  cout << "    Total associated=" << assoc_simple << endl;

  //  Stop timer and print results
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  printf("You have scanned %f Mbytes/Realtime seconds\n",mbytes/rtime);
  printf("You have scanned %f Mbytes/Cputime seconds\n",mbytes/ctime);
}