//=====================================================//
//Macro that reads the output of the flow analysis and 
//retrieves the QC and the MH containers.
//The macro produces an output called outputMH.root 
//that contains the four histograms related to the 
//3 particle correlator and the two integrated v2{2,QC} 
//and v2{4,QC}
//Author: Panos.Christakoglou@cern.ch

const Int_t nCentralityBins = 9;
TString strCentralityBins[nCentralityBins] = {"0-5","5-10","10-20",
					      "20-30","30-40","40-50",
					      "50-60","60-70","70-80"};

void plotMHCentrality(const char* filename = "AnalysisResults.root",
		      const char* systemType = "PbPb",
		      Bool_t isPID = kTRUE,
		      Int_t charge = 1) {
  gStyle->SetPalette(1,0);

  //----------------------------------------------------------
  // >>>>>>>>>>> Load libraries <<<<<<<<<<<<<< 
  //----------------------------------------------------------
  gSystem->AddIncludePath("-I$ROOTSYS/include");
  gSystem->Load("libGeom");
  gSystem->Load("libVMC");
  gSystem->Load("libXMLIO");
  gSystem->Load("libPhysics");
  
  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libPWGflowBase.so");
  gSystem->Load("libPWGflowTasks.so");

  AliFlowTrackCuts::PIDsource sourcePID = AliFlowTrackCuts::kTOFbayesian;
  AliPID::EParticleType particleType=AliPID::kPion;

  //----------------------------------------------------------
  // >>>>>>>>>>> Open file - Get objects <<<<<<<<<<<<<< 
  //----------------------------------------------------------
  TFile *fInput = TFile::Open(filename);
  if(!fInput) {
    Printf("File not found!!!");
    break;
  }
  //fInput->ls();

  //Get the TList of the MH
  TList *listMH = GetMHResults(fInput,systemType,isPID,
			       sourcePID,particleType,charge);

  //Get the TList of the QC
  TList *listQC = GetQCResults(fInput,systemType,isPID,
			       sourcePID,particleType,charge);

  //Get the TList of the SP
  //TList *listSP = GetSPResults(fInput,systemType,isPID,
  //sourcePID,particleType,charge);

  TFile *fOutput = TFile::Open("outputMH.root","recreate");
  listMH->Write();
  listQC->Write();
  listSP->Write();
  fOutput->Close();
}

//____________________________________________________________//
TList *GetSPResults(TFile *fInput,
		    const char* systemType = 0x0,
		    Bool_t isPID = kTRUE,
		    AliFlowTrackCuts::PIDsource sourcePID = AliFlowTrackCuts::kTOFbayesian,
		    AliPID::EParticleType particleType=AliPID::kPion,
		    Int_t charge = 0) {
  //Function that reads the TDirectoryFile of the MH
  //and returns a TList with the relevant plots.
  TList *listOutput = new TList();
  listOutput->SetName("listSP");

  //______________________________________________________________//
  //Global variables
  AliFlowCommonHistResults *commonHistRes;
  TString gAliFlowCommonHistName;

  //______________________________________________________________//
  //Get the TDirectoryFile
  TString directoryNameSP = 0;
  directoryNameSP = "outputSPanalysis"; 
  TDirectoryFile *outputSPanalysis = dynamic_cast<TDirectoryFile *>(fInput->Get(directoryNameSP.Data()));
  if(!outputSPanalysis) {
    Printf("SP directory not found!!!");
    break;
  }
  //outputSPanalysis->ls();

  //______________________________________________________________//
  //Get the TList
  TString listNameSP = 0;
  TList *cobjSP;
  for(Int_t iCentralityBin = 0; iCentralityBin < nCentralityBins; iCentralityBin++) {
    //for(Int_t iCentralityBin = 0; iCentralityBin < 1; iCentralityBin++) {
    listNameSP = "cobjSP_"; 
    listNameSP += systemType;
    listNameSP += "_";
    listNameSP += strCentralityBins[iCentralityBin];
    if(isPID) {
      listNameSP += AliFlowTrackCuts::PIDsourceName(sourcePID);
      listNameSP += "_";
      listNameSP += AliPID::ParticleName(particleType);
    }
    else {
      listNameSP += "AllCharged";
    }
    if (charge < 0) listNameSP += "-";
    if (charge > 0) listNameSP += "+" ;
    cobjSP = dynamic_cast<TList *>(outputSPanalysis->Get(listNameSP.Data()));
    if(!cobjSP) {
      Printf("SP object list not found!!!");
      break;
    }
    
    //Common hist for SP
    Double_t nAnalyzedEvents = 0.;
    AliFlowCommonHist *commonHistSP = dynamic_cast<AliFlowCommonHist *>(cobjSP->FindObject("AliFlowCommonHistSP"));
    if(commonHistSP) {
      TList *clistSPCommonHist = dynamic_cast<TList *>(commonHistSP->GetHistList());
      if(clistSPCommonHist) {
	TH1F *gHistMultRPSP = dynamic_cast<TH1F *>(clistSPCommonHist->FindObject("Control_Flow_MultRP_AliFlowCommonHistSP"));
	if(gHistMultRPSP)
	  nAnalyzedEvents = gHistMultRPSP->GetEntries();
      }
    }

    //Integrated flow RP
    gAliFlowCommonHistName = "AliFlowCommonHistResultsSP";
    commonHistRes = dynamic_cast<AliFlowCommonHistResults *>(cobjSP->FindObject(gAliFlowCommonHistName.Data()));
    //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[0].Data());
    TH1D *histSP = commonHistRes->GetHistIntFlowRP();
    histSP->SetName("fHistIntegratedFlowRPSP");
    
    //______________________________________________________________//
    //Print the results
    Printf("\n============================================================");
    Printf("Analyzed events: %d",(Int_t)nAnalyzedEvents);
    Printf("v2(SP): %lf +- %lf",histSP->GetBinContent(1),
	   histSP->GetBinError(1));
    Printf("============================================================");

    TList *dirOutput = new TList();
    dirOutput->SetName(listNameSP.Data());
    dirOutput->Add(histSP);
    listOutput->Add(dirOutput);
  }//loop over centrality bins TLists

  return listOutput;
}

//____________________________________________________________//
TList *GetQCResults(TFile *fInput,
		    const char* systemType = 0x0,
		    Bool_t isPID = kTRUE,
		    AliFlowTrackCuts::PIDsource sourcePID = AliFlowTrackCuts::kTOFbayesian,
		    AliPID::EParticleType particleType=AliPID::kPion,
		    Int_t charge = 0) {
  //Function that reads the TDirectoryFile of the QC
  //and returns a TList with the relevant plots.
  TList *listOutput = new TList();
  listOutput->SetName("listQC");

  //______________________________________________________________//
  //Global variables
  const Int_t nQCMethods = 4;
  AliFlowCommonHistResults *commonHistRes[nQCMethods];
  TString methodIntegratedFlowQC[nQCMethods] = {"2ndOrderQC",
						"4thOrderQC",
						"6thOrderQC",
						"8thOrderQC"};
  TString gAliFlowCommonHistName[nQCMethods];

  //______________________________________________________________//
  //Get the TDirectoryFile
  TString directoryNameQC = 0;
  directoryNameQC = "outputQCanalysis"; 
  TDirectoryFile *outputQCanalysis = dynamic_cast<TDirectoryFile *>(fInput->Get(directoryNameQC.Data()));
  if(!outputQCanalysis) {
    Printf("QC directory not found!!!");
    break;
  }
  //outputQCanalysis->ls();

  //______________________________________________________________//
  //Get the TList
  TString listNameQC = 0;
  TList *cobjQC;
  for(Int_t iCentralityBin = 0; iCentralityBin < nCentralityBins; iCentralityBin++) {
    //for(Int_t iCentralityBin = 0; iCentralityBin < 1; iCentralityBin++) {
    listNameQC = "cobjQC_"; 
    listNameQC += systemType;
    listNameQC += "_";
    listNameQC += strCentralityBins[iCentralityBin];
    if(isPID) {
      listNameQC += AliFlowTrackCuts::PIDsourceName(sourcePID);
      listNameQC += "_";
      listNameQC += AliPID::ParticleName(particleType);
    }
    else {
      listNameQC += "AllCharged";
    }
    if (charge < 0) listNameQC += "-";
    if (charge > 0) listNameQC += "+";

    cobjQC = dynamic_cast<TList *>(outputQCanalysis->Get(listNameQC.Data()));
    if(!cobjQC) {
      Printf("QC object list not found!!!");
      break;
    }
    
    //Common hist for QC
    Double_t nAnalyzedEvents = 0.;
    AliFlowCommonHist *commonHistQC = dynamic_cast<AliFlowCommonHist *>(cobjQC->FindObject("AliFlowCommonHistQC"));
    if(commonHistQC) {
      TList *clistQCCommonHist = dynamic_cast<TList *>(commonHistQC->GetHistList());
      if(clistQCCommonHist) {
	TH1F *gHistMultRPQC = dynamic_cast<TH1F *>(clistQCCommonHist->FindObject("Control_Flow_MultRP_AliFlowCommonHistQC"));
	if(gHistMultRPQC)
	  nAnalyzedEvents = gHistMultRPQC->GetEntries();
      }
    }

    //2nd order QC
    gAliFlowCommonHistName[0] = "AliFlowCommonHistResults";
    gAliFlowCommonHistName[0] += methodIntegratedFlowQC[0].Data();
    commonHistRes[0] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[0].Data()));
    //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[0].Data());
    TH1D *histQC_2 = commonHistRes[0]->GetHistIntFlowRP();
    histQC_2->SetName("fHistIntegratedFlowRPQC_2");
    
    //4th order QC
    gAliFlowCommonHistName[1] = "AliFlowCommonHistResults";
    gAliFlowCommonHistName[1] += methodIntegratedFlowQC[1].Data();
    commonHistRes[1] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[1].Data()));
    //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[1].Data());
    TH1D *histQC_4 = commonHistRes[1]->GetHistIntFlowRP();
    histQC_4->SetName("fHistIntegratedFlowRPQC_4");
    
    //6th order QC
    gAliFlowCommonHistName[2] = "AliFlowCommonHistResults";
    gAliFlowCommonHistName[2] += methodIntegratedFlowQC[2].Data();
    commonHistRes[2] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[2].Data()));
    //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[2].Data());
    TH1D *histQC_6 = commonHistRes[2]->GetHistIntFlowRP();
    histQC_6->SetName("fHistIntegratedFlowRPQC_6");
    
    //8th order QC
    gAliFlowCommonHistName[3] = "AliFlowCommonHistResults";
    gAliFlowCommonHistName[3] += methodIntegratedFlowQC[3].Data();
    commonHistRes[3] = dynamic_cast<AliFlowCommonHistResults *>(cobjQC->FindObject(gAliFlowCommonHistName[3].Data()));
    //Printf("AliFlowCommonHist name %s",gAliFlowCommonHistName[3].Data());
    TH1D *histQC_8 = commonHistRes[3]->GetHistIntFlowRP();
    histQC_8->SetName("fHistIntegratedFlowRPQC_8");
    
    TString gHistEventsName = "gHistEvents";
    gHistEventsName += strCentralityBins[iCentralityBin];
    TH1D *gHistEvents = new TH1D(gHistEventsName.Data(),
				 ";;N_{analyzed events}",
				 1,0.5,1.5);
    gHistEvents->SetBinContent(1,nAnalyzedEvents);

    //______________________________________________________________//
    //Print the results
    Printf("\n============================================================");
    Printf("Analyzed events: %d",(Int_t)nAnalyzedEvents);
    Printf("v2(2,QC): %lf +- %lf",histQC_2->GetBinContent(1),
	   histQC_2->GetBinError(1));
    Printf("v2(4,QC): %lf +- %lf",histQC_4->GetBinContent(1),
	   histQC_4->GetBinError(1));
    Printf("v2(6,QC): %lf +- %lf",histQC_6->GetBinContent(1),
	   histQC_6->GetBinError(1));
    Printf("v2(8,QC): %lf +- %lf",histQC_8->GetBinContent(1),
	   histQC_8->GetBinError(1));
    Printf("============================================================");

    TList *dirOutput = new TList();
    dirOutput->SetName(listNameQC.Data());
    dirOutput->Add(histQC_2);
    dirOutput->Add(histQC_4);
    dirOutput->Add(histQC_6);
    dirOutput->Add(histQC_8);
    dirOutput->Add(gHistEvents);
    listOutput->Add(dirOutput);
  }//loop over centrality bins TLists

  return listOutput;
}

//____________________________________________________________//
TList *GetMHResults(TFile *fInput,
		    const char* systemType = 0x0,
		    Bool_t isPID = kTRUE,
		    AliFlowTrackCuts::PIDsource sourcePID = AliFlowTrackCuts::kTOFbayesian,
		    AliPID::EParticleType particleType=AliPID::kPion,
		    Int_t charge = 0) {
  //Function that reads the TDirectoryFile of the MH
  //and returns a TList with the relevant plots.
  TList *listOutput = new TList();
  listOutput->SetName("listMH");

  //______________________________________________________________//
  //Get the TDirectoryFile
  TString directoryNameMH = 0;
  directoryNameMH = "outputMHanalysis"; 
  TDirectoryFile *outputMHanalysis = dynamic_cast<TDirectoryFile *>(fInput->Get(directoryNameMH.Data()));
  if(!outputMHanalysis) {
    Printf("MH directory not found!!!");
    break;
  }
  //outputMHanalysis->ls();

  //______________________________________________________________//
  //Get the TList
  TString listNameMH = 0;
  TList *cobjMH;
  TCanvas *c[nCentralityBins];
  for(Int_t iCentralityBin = 0; iCentralityBin < nCentralityBins; iCentralityBin++) {
    //for(Int_t iCentralityBin = 0; iCentralityBin < 1; iCentralityBin++) {
    listNameMH = "cobjMH_"; 
    listNameMH += systemType;
    listNameMH += "_";
    listNameMH += strCentralityBins[iCentralityBin];
    if(isPID) {
      listNameMH += AliFlowTrackCuts::PIDsourceName(sourcePID);
      listNameMH += "_";
      listNameMH += AliPID::ParticleName(particleType);
    }
    else {
      listNameMH += "AllCharged";
    }
    if (charge < 0) listNameMH += "-";
    if (charge > 0) listNameMH += "+";

    //Printf("%s",listNameMH.Data());
    cobjMH = dynamic_cast<TList *>(outputMHanalysis->Get(listNameMH.Data()));
    if(!cobjMH) {
      Printf("MH object list not found!!!");
      break;
    }
    //cobjMH->ls();
    
    //______________________________________________________________//
    //Get the daughter TLists
    TList *listWeights = dynamic_cast<TList *>(cobjMH->FindObject("Weights"));
    //listWeights->ls();
    TList *listProfiles = dynamic_cast<TList *>(cobjMH->FindObject("Profiles"));
    //listProfiles->ls();
    TList *listResults = dynamic_cast<TList *>(cobjMH->FindObject("Results"));
    //listResults->ls();

    if((!listWeights)||(!listProfiles)||(!listResults)) {
      Printf("MH output lists not found!!!");
      break;
    }
    
    //______________________________________________________________//
    TProfile *fAnalysisSettings = dynamic_cast<TProfile *>(cobjMH->FindObject("fAnalysisSettings"));

    //______________________________________________________________//
    //Get the objects from the Results list
    TH1D *f3pCorrelatorHist = dynamic_cast<TH1D *>(listResults->FindObject("f3pCorrelatorHist"));
    TH1D *fDetectorBiasHist = dynamic_cast<TH1D *>(listResults->FindObject("fDetectorBiasHist"));
    TH1D *f3pCorrelatorVsMHist = dynamic_cast<TH1D *>(listResults->FindObject("f3pCorrelatorVsMHist"));
    TH1D *fDetectorBiasVsMHist = dynamic_cast<TH1D *>(listResults->FindObject("fDetectorBiasVsMHist"));
    
    //______________________________________________________________//
    //Get the objects from the Profile list
    TProfile *f3pCorrelatorPro = dynamic_cast<TProfile *>(listProfiles->FindObject("f3pCorrelatorPro"));
    TProfile *fNonIsotropicTermsPro = dynamic_cast<TProfile *>(listProfiles->FindObject("fNonIsotropicTermsPro"));
    TProfile *f3pCorrelatorVsMPro = dynamic_cast<TProfile *>(listProfiles->FindObject("f3pCorrelatorVsMPro"));
    TProfile2D *fNonIsotropicTermsVsMPro = dynamic_cast<TProfile2D *>(listProfiles->FindObject("fNonIsotropicTermsVsMPro"));
    TProfile *f3pCorrelatorVsPtSumPro = dynamic_cast<TProfile *>(listProfiles->FindObject("f3pCorrelatorVsPtSumPro"));
    TProfile *f3pCorrelatorVsPtDiffPro = dynamic_cast<TProfile *>(listProfiles->FindObject("f3pCorrelatorVsPtDiffPro"));
    
    Double_t g3pCorrelatorValue = 0., g3pCorrelatorError = 0.;
    GetCorrelatorAndError(f3pCorrelatorVsPtSumPro,
			  g3pCorrelatorValue,
			  g3pCorrelatorError,1,17);
    TString g3pHistName = "g3pHistName";
    g3pHistName += strCentralityBins[iCentralityBin];
    TH1D *g3pHist = new TH1D(g3pHistName,
			     ";;#LT#LTcos(#psi_{1}+#psi_{2}-2#phi_{3}#GT#GT",
			     1,0.5,1.5);
    g3pHist->SetBinContent(1,g3pCorrelatorValue);
    g3pHist->SetBinError(1,g3pCorrelatorError);

    //______________________________________________________________//
    //Draw the differential 3p correlator
    c[iCentralityBin] = new TCanvas(listNameMH.Data(),listNameMH.Data(),
				    iCentralityBin*50,
				    iCentralityBin*50,
				    500,500);
    c[iCentralityBin]->SetHighLightColor(10); 
    c[iCentralityBin]->SetFillColor(10);
    f3pCorrelatorVsPtDiffPro->DrawCopy("E");

    //if(iCentralityBin == 8) {
    //for(Int_t iBin = 1; iBin <= f3pCorrelatorVsPtDiffPro->GetNbinsX(); iBin++) {
	//if(f3pCorrelatorVsPtDiffPro->GetBinContent(iBin) != 0.) 
	//Printf("Entries: %lf - Error: %lf - Value: %lf",
	//f3pCorrelatorVsPtDiffPro->GetBinEntries(iBin),
	//f3pCorrelatorVsPtDiffPro->GetBinError(iBin),
	//f3pCorrelatorVsPtDiffPro->GetBinContent(iBin));
	//}
    //}

    //______________________________________________________________//
    //Print the results
    Printf("============================================================");
    Printf("=========================Bin: %s=========================",strCentralityBins[iCentralityBin].Data());
    cout<<"<cos(psi1 + psi2 - 2phi3)>: "<<
      g3pCorrelatorValue <<
      " +- " <<
      g3pCorrelatorError << endl;
    cout<<"<cos(phi1 + phi2 - 2phi3)>: "<<
      f3pCorrelatorHist->GetBinContent(1) <<
      " +- " <<
      f3pCorrelatorHist->GetBinError(1)<<endl;
    Printf("============================================================");

    TList *dirOutput = new TList();
    dirOutput->SetName(listNameMH.Data());
    dirOutput->Add(f3pCorrelatorHist);
    //dirOutput->Add(f3pCorrelatorVsMHist);
    dirOutput->Add(f3pCorrelatorVsPtSumPro);
    dirOutput->Add(f3pCorrelatorVsPtDiffPro);
    dirOutput->Add(g3pHist);
    listOutput->Add(dirOutput);
  }//loop over centrality bins TLists

  //listOutput->ls();
  return listOutput;
}

//____________________________________________________________//
void GetCorrelatorAndError(TProfile *g3pCorrelatorVsPt,
			   Double_t &g3pCorrelatorValue,
			   Double_t &g3pCorrelatorError,
			   Int_t iBinLow = 0,
			   Int_t iBinHigh = 0) {
  //Function to return the average value of the 3p correlator 
  //<cos(psi1 + psi2 - 2phi3)> and its error.
  //The first argument is one of the 3p TProfile objects vs pt.
  //The second and third argument give the integral and its error.
  //The fourth and fifth, if specified, indicate the lowest and 
  //highest bin the calculation should be performed for.
  Int_t gBinLow = 1, gBinHigh = g3pCorrelatorVsPt->GetNbinsX();
  if(iBinLow) gBinLow = iBinLow;
  if(iBinHigh) gBinHigh = iBinHigh;
  
  Double_t gSumXi = 0.;
  Double_t gSumYi = 0.;
  Double_t gSumXiYi = 0.;
  Double_t gSumXiYi2 = 0.;
  Double_t gSumXi2Yi2 = 0.;
  Double_t gSumDeltaXi2 = 0.;
  Double_t gSumYi2DeltaXi2 = 0.;
  Double_t dError = 0.; //Flow code driven error calculation

  Double_t kSumBi = 0., kSumBi2DeltaYi2 = 0.;
  Double_t kSumYiBi = 0., kSumYi2DeltaBi2 = 0.;
  Double_t kSumDeltaBi2 = 0.;

 for(Int_t iBin = gBinLow; iBin <= gBinHigh; iBin++) {
    gSumXi += g3pCorrelatorVsPt->GetBinEntries(iBin);
    gSumYi += g3pCorrelatorVsPt->GetBinContent(iBin);
    gSumXiYi += g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin);
    gSumXiYi2 += g3pCorrelatorVsPt->GetBinEntries(iBin)*TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2);
    gSumXi2Yi2 += TMath::Power(g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin),2);
    gSumDeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
    gSumYi2DeltaXi2 += TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2) + TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);

    dError += g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinError(iBin)*g3pCorrelatorVsPt->GetBinError(iBin);  

    //new error calculation
    kSumBi += g3pCorrelatorVsPt->GetBinEntries(iBin);
    kSumYiBi += g3pCorrelatorVsPt->GetBinEntries(iBin)*g3pCorrelatorVsPt->GetBinContent(iBin);
    kSumBi2DeltaYi2 += TMath::Power(g3pCorrelatorVsPt->GetBinEntries(iBin),2)*TMath::Power(g3pCorrelatorVsPt->GetBinError(iBin),2);
    kSumYi2DeltaBi2 += TMath::Power(g3pCorrelatorVsPt->GetBinContent(iBin),2)*TMath::Power(TMath::Sqrt(g3pCorrelatorVsPt->GetBinEntries(iBin)),2);  
    kSumDeltaBi2 += TMath::Power(TMath::Sqrt(g3pCorrelatorVsPt->GetBinEntries(iBin)),2);
  }
  
  g3pCorrelatorValue = -1000.;
  g3pCorrelatorError = 1000.;
  
  if(gSumXi != 0.)
    g3pCorrelatorValue = gSumXiYi/gSumXi;
  if((gSumXi != 0.)&&(gSumXiYi != 0.))
    g3pCorrelatorError = TMath::Abs((gSumXiYi/gSumXi))*TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumYi2DeltaXi2)/gSumXiYi),2) + TMath::Power((gSumDeltaXi2/gSumXi),2));
    //g3pCorrelatorError = TMath::Sqrt((1./TMath::Power(kSumBi,2))*(kSumBi2DeltaYi2 + kSumYi2DeltaBi2 + TMath::Power(kSumYiBi,2)*kSumDeltaBi2/TMath::Power(kSumBi,2)));
  if(gSumXi != 0.)
    dError /= TMath::Power(gSumXi,2);
  dError = TMath::Sqrt(dError);
  g3pCorrelatorError = dError;

  return;
}