void  mycode(Int_t specie = 0, 
			     Int_t centMin = 0, Int_t centMax = 10, 
			     Bool_t bAnalyticCorr= 1 /*kTRUE*/, 
			     Bool_t bMix = 0,
			    //Bool_t bMerged = 1, 
			     Bool_t bAll = 1){

  //Note: not merged for the centrality bin of 10-40% for some cases

  gROOT->LoadMacro("fitHisto.C");

  if(!bAnalyticCorr){
    if (!bMix) {
      cout<<"Can not proceed further if the option is non-mixed!"<<endl;
      return;
    }
  }
}


  //char fileName[256];

  //TF1 * fAccCorr = new TF1("fAccCorr", "[0]-[1]*TMath::Abs(x)", -1.6, 1.6);
  //fAccCorr->SetParameters(1., 0.625);

  TString specieName[] = {"K0s", "Lambda", "AntiLambda"};

  // defining bins for Z vertex 
   const Int_t nVtx  = 8;
   Double_t vertexBins[] = {-8., -6., -4., -2., 0., 2., 4., 6., 8.};

  // pt bins of associated particles for the analysis                          
  const Int_t nPtAsso = 4;
   const Double_t PtBins[5] = {1.0, 2.0, 3.0, 4.0, 8.0};

  //  const Int_t nPtAsso = 1;
  //const Double_t PtBins[2] = {1.0, 2.0};

  // pt bins of trigger particles for the analysis                             
  const Int_t nPtTrig = 5;//4
   const Double_t PtBinsV0Xi[5] = {1.0, 2.0, 4.0, 8.0, 15.0};

  //  Int_t rebin = 2;
  const Int_t nPhi = 20;  
  const Int_t nEta = 32;
  // defining bins for centrality                                              
  //const Int_t nCent  = 2; 
  //Double_t centBins[] = {0.,10.,20.,30.,40.,50.,60., 90.};

  Double_t massMean = 0.493677;
  Double_t massRes = 0.005;
    
  //if(specie == 1 || specie == 2){
   // massMean = 1.115683;
   // massRes = 0.0015;
  //}

  Double_t fitRangeLow = massMean-15.*massRes;
  Double_t fitRangeHigh = massMean+15.*massRes;

  //TH1F * hNtrig;

  TH3F * h3EtaPhiSB[nVtx][nPtTrig][nPtAsso];
  
  TH2F * h2EtaPhiSB[nVtx][nPtTrig][nPtAsso];
  TH2F * h2EtaPhiSBMix[nVtx][nPtTrig][nPtAsso];

  TH3F * h3EtaPhiSBCor[nVtx][nPtTrig][nPtAsso]; //corrected by acceptance
  TH3F * h3EtaPhiSBCorAll[nPtTrig][nPtAsso] = {{NULL}}; // sum up vtx bins
  
  TH2F * h2EtaPhiSBPeakCor[nVtx][nPtTrig][nPtAsso];
  TH2F * h2EtaPhiSBPeakCorAll[nPtTrig][nPtAsso] = {{NULL}};

  TH2F * h2EtaPhiS[nVtx][nPtTrig][nPtAsso]; //after background subtraction

  TH2F * h2EtaPhiSPerTrig[nPtTrig][nPtAsso] = {{NULL}}; // per trigger phi correlation
  
  TH1F * hMass[nVtx][nPtTrig][nPtAsso];

  THnSparseF *fHistTrigSibAllK0s;
  THnSparseF *fHistdPhidEtaSibK0s[nVtx];
  THnSparseF *fHistdPhidEtaMixK0s[nVtx];

   TFile * file = new TFile("ResultData18JUL.root", "READ");

    //gDirectory->cd(Form("Cent%d_%d", centMin, centMax));
    TIter nextkey(gDirectory->GetListOfKeys());
    TKey* key;
    while (Key = (TKey*)nextkey()){
      if(!strcmp( key->GetName(), Form("OutputCent%d_%d", centMin, centMax))){
	TList * V0Xi = (TList*)key->ReadObj();
     fHistTrigSibAllK0s = (THnSparseF *)V0Xi->FindObject("K0s")->FindObject("fHistTrigSibAllK0s");
      if(!fHistTrigSibAllK0s) cout<<"Does not work!"<<endl;
      
      fHistdPhidEtaSibK0s = (THnSparseF *)V0Xi->FindObject("K0s")->FindObject("fHistdPhidEtaSibK0s");
      if(!fHistdPhidEtaSibK0s) cout<<"Does not work!"<<endl;

      fHistdPhidEtaMixK0s = (THnSparseF *)V0Xi->FindObject("K0s")->FindObject("fHistdPhidEtaMixK0s");
      if(!fHistdPhidEtaMixK0s) cout<<"Does not work!"<<endl;

	  }
	}
 

  for(Int_t iVtx = 0; iVtx < nVtx; iVtx++)
   for(Int_t iVtx = 0; iVtx < nVtx; iVtx++)
    for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++)
      for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){
	if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break;

	Int_t ptTrigBinLow 
	  = fHistTrigSibAllK0s[iVtx]->GetAxis(0)->FindBin(PtBinsV0Xi[iPtT]+0.01);
	Int_t ptTrigBinHigh
	  = fHistTrigSibAllK0s[iVtx]->GetAxis(0)->FindBin(PtBinsV0Xi[iPtT+1]-0.01); //changwan change
/*
	Int_t ptAssoBinLow
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01);
	Int_t ptAssoBinHigh 
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01);
*/
	fHistTrigSibAllK0s[iVtx]->GetAxis(0)->SetRange(ptTrigBinLow, ptTrigBinHigh);
//	fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh);

	hMass[iVtx][iPtT][iPtA] 
	  = (TH1F*)fHistTrigSibAllK0s[iVtx]->Projection(0, "OE");
	hMass[iVtx][iPtT][iPtA]
	  ->SetName(Form("MassCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f",
			 centMin, centMax,
			 vertexBins[iVtx], vertexBins[iVtx+1],
			 PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			 PtBins[iPtA]));
	hMass[iVtx][iPtT][iPtA]
	  ->SetTitle(Form("Mass, Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f]; mass [GeV/c^{2}]", 
			  centMin, centMax,
			  vertexBins[iVtx], vertexBins[iVtx+1],
			  PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			  PtBins[iPtA], PtBins[iPtA+1]));
      }
  
  for(Int_t iVtx = 0; iVtx < nVtx; iVtx++)
    for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++)
      for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){
	if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break;
	
	Int_t ptTrigBinLow
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT]+0.01);
	Int_t ptTrigBinHigh
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT+1]-0.01); //changwan change
	
	Int_t ptAssoBinLow
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01);
	Int_t ptAssoBinHigh
	  = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01);
	

	fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptTrigBinLow, ptTrigBinHigh);//changwan change
//	fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh);

	//	Int_t binLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->FindBin(massMean-3.*massRes);
	//Int_t binHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->FindBin(massMean+3.*massRes);
	
	//fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->SetRange(binLow, binHigh);
	h3EtaPhiSB[iVtx][iPtT][iPtA] 
	  = (TH3F*)fHistdPhidEtaSibK0s[iVtx]->Projection(5, 2, 3, "OE"); //changwan change
	h3EtaPhiSB[iVtx][iPtT][iPtA]
	  ->SetName(Form("EtaPhiSBCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", 
			 centMin, centMax, 
			 vertexBins[iVtx], vertexBins[iVtx+1],
			 PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			 PtBins[iPtA]));
	h3EtaPhiSB[iVtx][iPtT][iPtA]
	  ->SetTitle(Form("S. R., Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta #phi; #Delta#eta; mass",
			  centMin, centMax,
			  vertexBins[iVtx], vertexBins[iVtx+1],
			  PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			  PtBins[iPtA], PtBins[iPtA+1]));
	//h3EtaPhiSB[iVtx][iPtT][iPtA]->RebinX(rebin);
      }

  for(Int_t iVtx = 0; iVtx < nVtx; iVtx++)
    for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++)
      for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){
        if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break;

	Int_t binLow = h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->FindBin(massMean-3.*massRes);                                                
        Int_t binHigh = h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->FindBin(massMean+3.*massRes);
	h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->SetRange(binLow, binHigh);


	h2EtaPhiSB[iVtx][iPtT][iPtA]
          = (TH2F*)h3EtaPhiSB[iVtx][iPtT][iPtA]->Project3D("xy");
        h2EtaPhiSB[iVtx][iPtT][iPtA]
          ->SetName(Form("h2EtaPhiSBCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f",
                         centMin, centMax,
                         vertexBins[iVtx], vertexBins[iVtx+1],
                         PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
                         PtBins[iPtA]));
        h2EtaPhiSB[iVtx][iPtT][iPtA]
          ->SetTitle(Form("S. R., Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi",
                          centMin, centMax,
                          vertexBins[iVtx], vertexBins[iVtx+1],
                          PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
                          PtBins[iPtA], PtBins[iPtA+1]));

      }
  
  if(bMix)
    for(Int_t iVtx = 0; iVtx < nVtx; iVtx++)
      for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++)
	for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){
	  if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break;

       Int_t ptTrigBinLow = fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT]+0.01);
	  Int_t ptTrigBinHigh = fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT+1]-0.01);  //changwan change

	  //Int_t ptAssoBinLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01);
	 // Int_t ptAssoBinHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01);
	  
      fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(ptTrigBinLow, ptTrigBinHigh);//changwan change
	//  fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh);
	  
	  Int_t binLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(massMean-3.*massRes);
	  Int_t binHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(massMean+3.*massRes);
	  
	  fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(binLow, binHigh);
	  h2EtaPhiSBMix[iVtx][iPtT][iPtA]
	    = (TH2F*)fHistdPhidEtaMixK0s[iVtx]->Projection(5, 2, 3, "OE");
	  h2EtaPhiSBMix[iVtx][iPtT][iPtA]
	    ->SetName(Form("EtaPhiSBMixCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f",
			   centMin, centMax,
			   vertexBins[iVtx], vertexBins[iVtx+1],
			   PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			   PtBins[iPtA]));
	  h2EtaPhiSBMix[iVtx][iPtT][iPtA]
	    ->SetTitle(Form("Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi",
			    centMin, centMax,
			    vertexBins[iVtx], vertexBins[iVtx+1],
			    PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1],
			    PtBins[iPtA], PtBins[iPtA+1]));
	  //h2EtaPhiSBMix[iVtx][iPtT][iPtA]->RebinY(rebin);
	}