Using AliMagWrapCheb

[u/mrichter/AliRoot.git] / PWG2 / drawProtonQAResults.C

void drawProtonQAResults(const char* filename1 = "Protons.QA.root",
			 const char* filename2 = "Protons.MC.QA.root",
			 const char* filename3 = "Protons.Efficiency.root") {
  //Macro to visualize the results of the proton QA task
  //TCanvas objects: 15
  gStyle->SetPalette(1,0);
 
  TFile *fQA = TFile::Open(filename1);
  TList *listGlobalQA = (TList *)fQA->Get("globalQAList");
  drawCutStatistics(listGlobalQA);
 
  TFile *fMC = TFile::Open(filename2);
  TList *listPDG = (TList *)fMC->Get("pdgCodeList");  
  TList *listMCProcesses = (TList *)fMC->Get("mcProcessList");  
  drawMCQA(listPDG,listMCProcesses);
  
  TFile *fEfficiency = TFile::Open(filename3);
  TList *listEfficiency = (TList *)fEfficiency->Get("efficiencyList");
  drawEfficiency(listEfficiency);

  fQA->Close();
  fMC->Close();
  fEfficiency->Close();
}

//________________________________________//
void drawCutStatistics(TList *list) {
  //Function to display the statistics from the cuts
  //The histogram shows the influence of each cut on the primary
  //and secondary (anti)protons
  const Int_t NQAHISTOSPERLIST = 26;
  
  Double_t gEntriesQA2DList[12] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  Double_t gEntriesQAPrimaryProtonsAcceptedList[NQAHISTOSPERLIST], gEntriesQAPrimaryProtonsRejectedList[NQAHISTOSPERLIST];
  Double_t gEntriesQASecondaryProtonsAcceptedList[NQAHISTOSPERLIST], gEntriesQASecondaryProtonsRejectedList[NQAHISTOSPERLIST];
  Double_t gEntriesQAPrimaryAntiProtonsAcceptedList[NQAHISTOSPERLIST], gEntriesQAPrimaryAntiProtonsRejectedList[NQAHISTOSPERLIST];
  Double_t gEntriesQASecondaryAntiProtonsAcceptedList[NQAHISTOSPERLIST], gEntriesQASecondaryAntiProtonsRejectedList[NQAHISTOSPERLIST];

  for(Int_t i = 0; i < NQAHISTOSPERLIST; i++) {
    gEntriesQAPrimaryProtonsAcceptedList[i] = 0.0;
    gEntriesQAPrimaryProtonsRejectedList[i] = 0.0;
    gEntriesQASecondaryProtonsAcceptedList[i] = 0.0;
    gEntriesQASecondaryProtonsRejectedList[i] = 0.0;
    gEntriesQAPrimaryAntiProtonsAcceptedList[i] = 0.0;
    gEntriesQAPrimaryAntiProtonsRejectedList[i] = 0.0;
    gEntriesQASecondaryAntiProtonsAcceptedList[i] = 0.0;
    gEntriesQASecondaryAntiProtonsRejectedList[i] = 0.0;
  }

  TList *fQA2DList = (TList *)list->At(0);
  GetQAEntries(fQA2DList,gEntriesQA2DList);
  TH3F *gHistYPtPDGProtonsPass = (TH3F *)fQA2DList->At(10);
  TH3F *gHistYPtPDGAntiProtonsPass = (TH3F *)fQA2DList->At(11);

  TList *fQAPrimaryProtonsAcceptedList = (TList *)list->At(1);
  GetQAEntries(fQAPrimaryProtonsAcceptedList,gEntriesQAPrimaryProtonsAcceptedList);

  TList *fQAPrimaryProtonsRejectedList = (TList *)list->At(2);
  GetQAEntries(fQAPrimaryProtonsRejectedList,gEntriesQAPrimaryProtonsRejectedList);

  TList *fQASecondaryProtonsAcceptedList = (TList *)list->At(3);
  GetQAEntries(fQASecondaryProtonsAcceptedList,gEntriesQASecondaryProtonsAcceptedList);

  TList *fQASecondaryProtonsRejectedList = (TList *)list->At(4);
  GetQAEntries(fQASecondaryProtonsRejectedList,gEntriesQASecondaryProtonsRejectedList);

  TList *fQAPrimaryAntiProtonsAcceptedList = (TList *)list->At(5);
  GetQAEntries(fQAPrimaryAntiProtonsAcceptedList,gEntriesQAPrimaryAntiProtonsAcceptedList);

  TList *fQAPrimaryAntiProtonsRejectedList = (TList *)list->At(6);
  GetQAEntries(fQAPrimaryAntiProtonsRejectedList,gEntriesQAPrimaryAntiProtonsRejectedList);

  TList *fQASecondaryAntiProtonsAcceptedList = (TList *)list->At(7);
  GetQAEntries(fQASecondaryAntiProtonsAcceptedList,gEntriesQASecondaryAntiProtonsAcceptedList);

  TList *fQASecondaryAntiProtonsRejectedList = (TList *)list->At(8);
  GetQAEntries(fQASecondaryAntiProtonsRejectedList,gEntriesQASecondaryAntiProtonsRejectedList);

  //_______________________________________________________//
  //Create the histograms
  const Int_t nx = 26;
  char *fCutName[nx] = {"Tracks","",
			"ITS Clusters",
			"#chi^{2}/N_{ITS-Clusters}",
			"TPC Clusters",
			"#chi^{2}/N_{TPC-Clusters}",
			"ExtCov11",
			"ExtCov22",
			"ExtCov33",
			"ExtCov44",
			"ExtCov55",
			"#sigma_{Vertex}",
			"#sigma_{Vertex-TPC}",
			"DCA_{xy}",
			"DCA_{xy}(TPC)",
			"DCA_{z}",
			"DCA_{z}(TPC)",
			"#chi^{2}(vertex)",
			"ITS refit",
			"TPC refit",
			"ESD pid",
			"TPC pid","",
			"N_{Secondaries}/N_{total}","",""};
  char *fCutITSName[6] = {"SPD_{1}","SPD_{2}",
			  "SDD_{1}","SDD_{2}",
			  "SSD_{1}","SSD_{2}"};

  //cut influence
  TH2F *hEmpty = new TH2F("hEmpty","",nx,0,nx,100,0,100); 
  hEmpty->SetStats(kFALSE); 
  hEmpty->SetMarkerStyle(kFullCircle);
  hEmpty->SetLineColor(2); hEmpty->SetMarkerColor(2);
  hEmpty->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= nx; i++) 
    hEmpty->GetXaxis()->SetBinLabel(i,fCutName[i-1]);
  hEmpty->GetXaxis()->SetLabelOffset(0.01); 
  hEmpty->GetXaxis()->SetLabelSize(0.045);

  //primary protons
  TH1F *hPrimaryProtons = new TH1F("hPrimaryProtons","",nx,0,nx); 
  hPrimaryProtons->SetStats(kFALSE); 
  hPrimaryProtons->SetMarkerStyle(kFullCircle);
  hPrimaryProtons->SetLineColor(4); hPrimaryProtons->SetLineWidth(2); 
  hPrimaryProtons->SetMarkerColor(4);
  hPrimaryProtons->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= nx; i++) {
    hPrimaryProtons->SetBinContent(i,-10.);
    hPrimaryProtons->GetXaxis()->SetBinLabel(i,fCutName[i-1]);
  }
  hPrimaryProtons->GetXaxis()->SetLabelOffset(0.01); 
  hPrimaryProtons->GetXaxis()->SetLabelSize(0.045);

  //secondary protons
  TH1F *hSecondaryProtons = new TH1F("hSecondaryProtons","",nx,0,nx); 
  hSecondaryProtons->SetStats(kFALSE); 
  hSecondaryProtons->SetMarkerStyle(22);
  hSecondaryProtons->SetMarkerSize(1.4);
  hSecondaryProtons->SetLineColor(2); hSecondaryProtons->SetLineWidth(2); 
  hSecondaryProtons->SetMarkerColor(2);
  hSecondaryProtons->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= nx; i++) {
    hSecondaryProtons->SetBinContent(i,-10.);
    hSecondaryProtons->GetXaxis()->SetBinLabel(i,fCutName[i-1]);
  }
  hSecondaryProtons->GetXaxis()->SetLabelOffset(0.01); 
  hSecondaryProtons->GetXaxis()->SetLabelSize(0.045);

  //primary antiprotons
  TH1F *hPrimaryAntiProtons = new TH1F("hPrimaryAntiProtons","",nx,0,nx); 
  hPrimaryAntiProtons->SetStats(kFALSE); 
  hPrimaryAntiProtons->SetMarkerStyle(kFullCircle);
  hPrimaryAntiProtons->SetLineColor(4); hPrimaryAntiProtons->SetLineWidth(2); 
  hPrimaryAntiProtons->SetMarkerColor(4);
  hPrimaryAntiProtons->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= nx; i++) {
    hPrimaryAntiProtons->SetBinContent(i,-10.);
    hPrimaryAntiProtons->GetXaxis()->SetBinLabel(i,fCutName[i-1]);
  }
  hPrimaryAntiProtons->GetXaxis()->SetLabelOffset(0.01); 
  hPrimaryAntiProtons->GetXaxis()->SetLabelSize(0.045);

  //secondary antiprotons
  TH1F *hSecondaryAntiProtons = new TH1F("hSecondaryAntiProtons","",nx,0,nx); 
  hSecondaryAntiProtons->SetStats(kFALSE); 
  hSecondaryAntiProtons->SetMarkerStyle(22);
  hSecondaryAntiProtons->SetMarkerSize(1.4);
  hSecondaryAntiProtons->SetLineColor(2); hSecondaryAntiProtons->SetLineWidth(2); 
  hSecondaryAntiProtons->SetMarkerColor(2);
  hSecondaryAntiProtons->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= nx; i++) {
    hSecondaryAntiProtons->SetBinContent(i,-10.);
    hSecondaryAntiProtons->GetXaxis()->SetBinLabel(i,fCutName[i-1]);
  }
  hSecondaryAntiProtons->GetXaxis()->SetLabelOffset(0.01); 
  hSecondaryAntiProtons->GetXaxis()->SetLabelSize(0.045);
  //_______________________________________________________//

  //1D for primary protons
  //cout<<"_____________________________________________________"<<endl;
  //cout<<"_______________PRIMARY PROTONS_______________________"<<endl;
  hPrimaryProtons->SetBinContent(1,GetPercentage(gEntriesQA2DList[0],gEntriesQA2DList[1]));

  for(Int_t i = 2; i < NQAHISTOSPERLIST-4; i++) 
    hPrimaryProtons->SetBinContent(i+1,GetPercentage(gEntriesQAPrimaryProtonsAcceptedList[i-2],
						     gEntriesQAPrimaryProtonsRejectedList[i-2]));
  
  //1D for secondary protons
  //cout<<"_____________________________________________________"<<endl;
  //cout<<"_______________SECONDARY PROTONS_____________________"<<endl;
  hSecondaryProtons->SetBinContent(1,GetPercentage(gEntriesQA2DList[2],gEntriesQA2DList[3]));
  
  for(Int_t i = 2; i < NQAHISTOSPERLIST-4; i++) 
    hSecondaryProtons->SetBinContent(i+1,GetPercentage(gEntriesQASecondaryProtonsAcceptedList[i-2],
						       gEntriesQASecondaryProtonsRejectedList[i-2]));
  hSecondaryProtons->SetBinContent(24,GetPercentage(gEntriesQA2DList[0],gEntriesQA2DList[2]));

  //1D for primary antiprotons
  //cout<<"_________________________________________________________"<<endl;
  //cout<<"_______________PRIMARY ANTIPROTONS_______________________"<<endl;
  hPrimaryAntiProtons->SetBinContent(1,GetPercentage(gEntriesQA2DList[4],gEntriesQA2DList[5]));
  
  for(Int_t i = 2; i < NQAHISTOSPERLIST-4; i++) 
    hPrimaryAntiProtons->SetBinContent(i+1,GetPercentage(gEntriesQAPrimaryAntiProtonsAcceptedList[i-2],
							 gEntriesQAPrimaryAntiProtonsRejectedList[i-2]));
  
  //1D for secondary antiprotons
  //cout<<"_________________________________________________________"<<endl;
  //cout<<"_______________SECONDARY ANTIPROTONS_____________________"<<endl;
  hSecondaryAntiProtons->SetBinContent(1,GetPercentage(gEntriesQA2DList[6],gEntriesQA2DList[7]));

  for(Int_t i = 2; i < NQAHISTOSPERLIST-4; i++) 
    hSecondaryAntiProtons->SetBinContent(i+1,GetPercentage(gEntriesQASecondaryAntiProtonsAcceptedList[i-2],
							   gEntriesQASecondaryAntiProtonsRejectedList[i-2]));
  hSecondaryAntiProtons->SetBinContent(24,GetPercentage(gEntriesQA2DList[4],gEntriesQA2DList[6]));

  TLatex *t1 = new TLatex();
  t1->SetTextSize(0.04);
  //_________________________________________________________//
  TCanvas *c1 = new TCanvas("c1","Cut Influence - Protons",0,0,700,400);
  c1->SetFillColor(10); c1->GetFrame()->SetFillColor(10);
  c1->SetHighLightColor(10); c1->SetBottomMargin(0.15);
  c1->SetGridx(); c1->SetGridy();

  for(Int_t i = 1; i <= nx; i++) {
    hPrimaryProtons->SetBinError(i,1.0);
    hSecondaryProtons->SetBinError(i,1.0);
  }
  hEmpty->DrawCopy();
  hPrimaryProtons->DrawCopy("EHISTSAME");
  hSecondaryProtons->DrawCopy("EHISTSAME");
  DrawMarker(20.5, 90, 20, 1.2, 4);
  t1->DrawLatex(21,88,"Primary p");
  DrawMarker(20.5, 80, 22, 1.2, 2);
  t1->DrawLatex(21,78,"Secondary p");

  c1->SaveAs("CutInfluence-Protons.gif");

  //_________________________________________________________//
  TCanvas *c2 = new TCanvas("c2","Cut Influence - AntiProtons",50,50,700,400);
  c2->SetFillColor(10); c2->GetFrame()->SetFillColor(10);
  c2->SetHighLightColor(10); c2->SetBottomMargin(0.15);
  c2->SetGridx(); c2->SetGridy();

  for(Int_t i = 1; i <= nx; i++) {
    hPrimaryAntiProtons->SetBinError(i,1.0);
    hSecondaryAntiProtons->SetBinError(i,1.0);
  }
  hEmpty->DrawCopy();
  hPrimaryAntiProtons->DrawCopy("EHISTSAME");
  hSecondaryAntiProtons->DrawCopy("EHISTSAME");
  DrawMarker(20.5, 90, 20, 1.2, 4);
  t1->DrawLatex(21,88,"Primary #bar{p}");
  DrawMarker(20.5, 80, 22, 1.2, 2);
  t1->DrawLatex(21,78,"Secondary #bar{p}");

  c2->SaveAs("CutInfluence-AntiProtons.gif");

  //_________________________________________________________//
  //ITS layers influence
  TH2F *hEmptyITS = new TH2F("hEmptyITS","",10,0,10,100,0,100); 
  hEmptyITS->SetStats(kFALSE); 
  hEmptyITS->SetMarkerStyle(kFullCircle);
  hEmptyITS->SetLineColor(2); hEmptyITS->SetMarkerColor(2);
  hEmptyITS->GetYaxis()->SetTitle("Influence of the ITS layers [%]");
  hEmptyITS->GetYaxis()->SetTitleOffset(1.3);
  for(Int_t i = 1; i <= 6; i++) 
    hEmptyITS->GetXaxis()->SetBinLabel(i,fCutITSName[i-1]);
  hEmptyITS->GetXaxis()->SetLabelOffset(0.01); 
  hEmptyITS->GetXaxis()->SetLabelSize(0.045);

  //_________________________________________________________//
  //primary protons
  TH1F *hPrimaryProtonsITS = new TH1F("hPrimaryProtonsITS","",7,0,7); 
  hPrimaryProtonsITS->SetStats(kFALSE); 
  hPrimaryProtonsITS->SetMarkerStyle(kFullCircle);
  hPrimaryProtonsITS->SetLineColor(4); hPrimaryProtonsITS->SetLineWidth(2); 
  hPrimaryProtonsITS->SetMarkerColor(4);
  hPrimaryProtonsITS->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= 6; i++) {
    hPrimaryProtonsITS->SetBinContent(i,-10.);
    hPrimaryProtonsITS->GetXaxis()->SetBinLabel(i,fCutITSName[i-1]);
  }
  hPrimaryProtonsITS->GetXaxis()->SetLabelOffset(0.01); 
  hPrimaryProtonsITS->GetXaxis()->SetLabelSize(0.045);

  //secondary protons
  TH1F *hSecondaryProtonsITS = new TH1F("hSecondaryProtonsITS","",7,0,7); 
  hSecondaryProtonsITS->SetStats(kFALSE); 
  hSecondaryProtonsITS->SetMarkerStyle(22);
  hSecondaryProtonsITS->SetMarkerSize(1.4);
  hSecondaryProtonsITS->SetLineColor(2); hSecondaryProtonsITS->SetLineWidth(2); 
  hSecondaryProtonsITS->SetMarkerColor(2);
  hSecondaryProtonsITS->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= 6; i++) {
    hSecondaryProtonsITS->SetBinContent(i,-10.);
    hSecondaryProtonsITS->GetXaxis()->SetBinLabel(i,fCutITSName[i-1]);
  }
  hSecondaryProtonsITS->GetXaxis()->SetLabelOffset(0.01); 
  hSecondaryProtonsITS->GetXaxis()->SetLabelSize(0.045);

  //primary antiprotons
  TH1F *hPrimaryAntiProtonsITS = new TH1F("hPrimaryAntiProtonsITS","",7,0,7); 
  hPrimaryAntiProtonsITS->SetStats(kFALSE); 
  hPrimaryAntiProtonsITS->SetMarkerStyle(kFullCircle);
  hPrimaryAntiProtonsITS->SetLineColor(4); hPrimaryAntiProtonsITS->SetLineWidth(2); 
  hPrimaryAntiProtonsITS->SetMarkerColor(4);
  hPrimaryAntiProtonsITS->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= 6; i++) {
    hPrimaryAntiProtonsITS->SetBinContent(i,-10.);
    hPrimaryAntiProtonsITS->GetXaxis()->SetBinLabel(i,fCutITSName[i-1]);
  }
  hPrimaryAntiProtonsITS->GetXaxis()->SetLabelOffset(0.01); 
  hPrimaryAntiProtonsITS->GetXaxis()->SetLabelSize(0.045);

  //secondary antiprotons
  TH1F *hSecondaryAntiProtonsITS = new TH1F("hSecondaryAntiProtonsITS","",7,0,7); 
  hSecondaryAntiProtonsITS->SetStats(kFALSE); 
  hSecondaryAntiProtonsITS->SetMarkerStyle(22);
  hSecondaryAntiProtonsITS->SetMarkerSize(1.4);
  hSecondaryAntiProtonsITS->SetLineColor(2); hSecondaryAntiProtonsITS->SetLineWidth(2); 
  hSecondaryAntiProtonsITS->SetMarkerColor(2);
  hSecondaryAntiProtonsITS->GetYaxis()->SetTitle("Influence of the cuts [%]");
  for(Int_t i = 1; i <= 6; i++) {
    hSecondaryAntiProtonsITS->SetBinContent(i,-10.);
    hSecondaryAntiProtonsITS->GetXaxis()->SetBinLabel(i,fCutITSName[i-1]);
  }
  hSecondaryAntiProtonsITS->GetXaxis()->SetLabelOffset(0.01); 
  hSecondaryAntiProtonsITS->GetXaxis()->SetLabelSize(0.045);

  //_______________________________________________________//
  TCanvas *c9 = new TCanvas("c9","ITS cluster map - (anti)protons",
			    100,100,700,400);
  c9->SetFillColor(10); c9->GetFrame()->SetFillColor(10);
  c9->SetHighLightColor(10); c9->Divide(2,1);
  for(Int_t i = 1; i <= 6; i++) {
    hPrimaryProtonsITS->SetBinError(i,1.0);
    hSecondaryProtonsITS->SetBinError(i,1.0);
    hPrimaryAntiProtonsITS->SetBinError(i,1.0);
    hSecondaryAntiProtonsITS->SetBinError(i,1.0);
  }
  c9->cd(1)->SetBottomMargin(0.15);
  c9->cd(1)->SetLeftMargin(0.15);
  c9->cd(1)->SetGridx(); c9->cd(1)->SetGridy();
  hEmptyITS->SetTitle("Protons");
  hEmptyITS->DrawCopy();

  for(Int_t i = 1; i < 7; i++) {
    hPrimaryProtonsITS->SetBinContent(i,GetPercentage(gEntriesQAPrimaryProtonsAcceptedList[19+i],
						      gEntriesQAPrimaryProtonsRejectedList[19+i]));
    hSecondaryProtonsITS->SetBinContent(i,GetPercentage(gEntriesQASecondaryProtonsAcceptedList[19+i],
							gEntriesQASecondaryProtonsRejectedList[19+i]));
  }
  hPrimaryProtonsITS->DrawCopy("EHISTSAME");
  hSecondaryProtonsITS->DrawCopy("EHISTSAME");
  DrawMarker(6.5, 90, 20, 1.2, 4);
  t1->DrawLatex(7,88,"Primary p");
  DrawMarker(6.5, 80, 22, 1.2, 2);
  t1->DrawLatex(7,78,"Secondary p");

  c9->cd(2)->SetBottomMargin(0.15);
  c9->cd(2)->SetLeftMargin(0.15);
  c9->cd(2)->SetGridx(); c9->cd(2)->SetGridy();
  hEmptyITS->SetTitle("Antiprotons");
  hEmptyITS->DrawCopy();
  for(Int_t i = 1; i < 7; i++) {
    hPrimaryAntiProtonsITS->SetBinContent(i,GetPercentage(gEntriesQAPrimaryAntiProtonsAcceptedList[19+i],
							  gEntriesQAPrimaryAntiProtonsRejectedList[19+i]));
    hSecondaryAntiProtonsITS->SetBinContent(i,GetPercentage(gEntriesQASecondaryAntiProtonsAcceptedList[19+i],
							    gEntriesQASecondaryAntiProtonsRejectedList[19+i]));
  }
  hPrimaryAntiProtonsITS->DrawCopy("EHISTSAME");
  hSecondaryAntiProtonsITS->DrawCopy("EHISTSAME");
  DrawMarker(6.5, 90, 20, 1.2, 4);
  t1->DrawLatex(7,88,"Primary #bar{p}");
  DrawMarker(6.5, 80, 22, 1.2, 2);
  t1->DrawLatex(7,78,"Secondary #bar{p}");

  c9->SaveAs("CutInfluence-ITS.gif");
  
  //Efficiency - Contamination plots
  DrawContamination(fQA2DList);
  DrawCutEfficiency(fQA2DList);
  DrawComposition(gHistYPtPDGProtonsPass,gHistYPtPDGAntiProtonsPass);
}

//________________________________________//
void DrawComposition(TH3F *gHistYPtPDGProtons,
		     TH3F *gHistYPtPDGAntiProtons) {
  //Function to display the composition of secondary (anti)protons
  //that survive the quality criteria
  Double_t nParticleCompositionProtonY[100], nParticleCompositionProtonPt[100];
  Double_t nParticleCompositionAntiProtonY[100], nParticleCompositionAntiProtonPt[100];
  Double_t gY[100], gPt[100];
  for(Int_t iBins = 0; iBins < 100; iBins++) {
    nParticleCompositionProtonY[iBins] = 0;
    nParticleCompositionProtonPt[iBins] = 0;
    nParticleCompositionAntiProtonY[iBins] = 0;
    nParticleCompositionAntiProtonPt[iBins] = 0;
    gY[iBins] = 0;
    gPt[iBins] = 0;
  }
  
  TGraph *gParticleProtonY[14];
  TGraph *gParticleProtonPt[14];
  TGraph *gParticleAntiProtonY[14];
  TGraph *gParticleAntiProtonPt[14];
  for(Int_t iParticle = 0; iParticle < 14; iParticle++) {
    GetComposition(iParticle,
		   gHistYPtPDGProtons,
		   nParticleCompositionProtonY,
		   gY, nParticleCompositionProtonPt, gPt);
    gParticleProtonY[iParticle] = new TGraph(gHistYPtPDGProtons->GetNbinsX(),
				       gY,nParticleCompositionProtonY);
    gParticleProtonY[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleProtonY[iParticle]->SetMarkerSize(1.2);

    gParticleProtonPt[iParticle] = new TGraph(gHistYPtPDGProtons->GetNbinsY(),
					gPt,nParticleCompositionProtonPt);
    gParticleProtonPt[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleProtonPt[iParticle]->SetMarkerSize(1.2);

    GetComposition(iParticle,
		   gHistYPtPDGAntiProtons,
		   nParticleCompositionAntiProtonY,
		   gY, nParticleCompositionAntiProtonPt, gPt);
    gParticleAntiProtonY[iParticle] = new TGraph(gHistYPtPDGAntiProtons->GetNbinsX(),
				       gY,nParticleCompositionAntiProtonY);
    gParticleAntiProtonY[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleAntiProtonY[iParticle]->SetMarkerSize(1.2);

    gParticleAntiProtonPt[iParticle] = new TGraph(gHistYPtPDGAntiProtons->GetNbinsY(),
					gPt,nParticleCompositionAntiProtonPt);
    gParticleAntiProtonPt[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleAntiProtonPt[iParticle]->SetMarkerSize(1.2);
  }

  //_________________________________________________________//
  char *fParticleName[14] = {"Primary","K_{L}","#pi","K_{S}","K",
			     "n","p","#Sigma^{-}","#Lambda","#Sigma^{+}",
			     "#Xi^{-}","#Xi^{0}","#Omega^{-}"};
  TLatex *t1 = new TLatex();
  t1->SetTextSize(0.04);

  TH2F *hEmptyY = new TH2F("hEmptyY","",100,-1.2,1.2,100,0,120); 
  hEmptyY->SetStats(kFALSE); 
  hEmptyY->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyY->GetXaxis()->SetTitle("y");

  TCanvas *c12 = new TCanvas("c12",
			     "Composition of accepted secondaries vs y",
			    550,550,700,400);
  c12->SetFillColor(10); c12->GetFrame()->SetFillColor(10); c12->Divide(2,1);
  c12->SetHighLightColor(10); c12->cd(1)->SetBottomMargin(0.15);
  c12->cd(1)->SetGridx(); c12->cd(1)->SetGridy();
  hEmptyY->SetTitle("Protons");
  hEmptyY->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    //if((iParticle == 0)||(iParticle == 2)||(iParticle == 5)||(iParticle == 6)||(iParticle == 8))
    if((iParticle == 0)||(iParticle == 2)||(iParticle == 6)||(iParticle == 8))
      gParticleProtonY[iParticle]->Draw("P");
    if(iParticle < 5) {
      DrawMarker(-1.1, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(-1.0,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(0.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
  }

  c12->SetHighLightColor(10); c12->cd(2)->SetBottomMargin(0.15);
  c12->cd(2)->SetGridx(); c12->cd(2)->SetGridy();
  hEmptyY->SetTitle("Antiprotons");
  hEmptyY->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if((iParticle == 0)||(iParticle == 6)||(iParticle == 8))
      gParticleAntiProtonY[iParticle]->Draw("P");
    if(iParticle < 5) {
      DrawMarker(-1.1, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(-1.0,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(0.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
  }
  c12->SaveAs("SurvivedSecondaries-Composition-Rapidity.gif");

  TH2F *hEmptyPt = new TH2F("hEmptyPt","",100,0.0,4.0,100,0,120); 
  hEmptyPt->SetStats(kFALSE); 
  hEmptyPt->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");

  TCanvas *c13 = new TCanvas("c13",
			     "Composition of accepted secondaries vs pT",
			     600,600,700,400);
  c13->SetFillColor(10); c13->GetFrame()->SetFillColor(10); c13->Divide(2,1);
  c13->SetHighLightColor(10); c13->cd(1)->SetBottomMargin(0.15);
  c13->cd(1)->SetGridx(); c13->cd(1)->SetGridy();
  hEmptyPt->GetXaxis()->SetRangeUser(gParticleProtonPt[0]->GetXaxis()->GetXmin()-0.2,gParticleProtonPt[0]->GetXaxis()->GetXmax()+0.2);
  hEmptyPt->SetTitle("Protons");
  hEmptyPt->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if(iParticle < 5) {
      DrawMarker(gParticleProtonPt[0]->GetXaxis()->GetXmin()+0.1, 
		 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(gParticleProtonPt[0]->GetXaxis()->GetXmin()+0.2,
		    113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(gParticleProtonPt[0]->GetXaxis()->GetXmax()*0.5, 
		 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(gParticleProtonPt[0]->GetXaxis()->GetXmax()*0.5+0.1,
		    113-5*(iParticle-5),fParticleName[iParticle]);
    }
    if((iParticle == 0)||(iParticle == 2)||(iParticle == 6)||(iParticle == 8))
      gParticleProtonPt[iParticle]->Draw("P");
  }

  c13->SetHighLightColor(10); c13->cd(2)->SetBottomMargin(0.15);
  c13->cd(2)->SetGridx(); c13->cd(2)->SetGridy();
  hEmptyPt->SetTitle("Antiprotons");
  hEmptyPt->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if(iParticle < 5) {
      DrawMarker(gParticleProtonPt[0]->GetXaxis()->GetXmin()+0.1, 
		 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(gParticleProtonPt[0]->GetXaxis()->GetXmin()+0.2,
		    113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(gParticleProtonPt[0]->GetXaxis()->GetXmax()*0.5, 
		 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(gParticleProtonPt[0]->GetXaxis()->GetXmax()*0.5+0.1,
		    113-5*(iParticle-5),fParticleName[iParticle]);
    }
    if((iParticle == 0)||(iParticle == 6)||(iParticle == 8))
      gParticleAntiProtonPt[iParticle]->Draw("P");
  }
  c13->SaveAs("SurvivedSecondaries-Composition-Pt.gif");
}

//________________________________________________//
void DrawContamination(TList *inputList) {
  //loops over the list entries and
  //draws the rapidity and pT dependence
  //of the percentage of primary and secondary
  //protons and antiprotons after the track cuts
  cout<<"Extracting the entries for the histograms in the list: "<<
    inputList->GetName()<<"..."<<endl;

  TLatex *t1 = new TLatex();
  t1->SetTextSize(0.04);

  TH2F *hPrimaryProtons = (TH2F *)inputList->At(0);
  TH2F *hSecondaryProtons = (TH2F *)inputList->At(2);
  TH2F *hPrimaryAntiProtons = (TH2F *)inputList->At(4);
  TH2F *hSecondaryAntiProtons = (TH2F *)inputList->At(6);

  //rapidity dependence
  //Protons
  TH1D *gYPrimaryProtons = (TH1D *)hPrimaryProtons->ProjectionX("gYPrimaryProtons",0,hPrimaryProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYSecondaryProtons = (TH1D *)hSecondaryProtons->ProjectionX("gYSecondaryProtons",0,hSecondaryProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYTotalProtons = (TH1D *)hPrimaryProtons->ProjectionX("gYTotalProtons",0,hPrimaryProtons->GetXaxis()->GetNbins(),"e");
  gYTotalProtons->Add(gYSecondaryProtons);

  TH1D *gYPrimaryProtonsPercentage = new TH1D("gYPrimaryProtonsPercentage",
					      "",
					      hPrimaryProtons->GetXaxis()->GetNbins(),
					      hPrimaryProtons->GetXaxis()->GetXmin(),
					      hPrimaryProtons->GetXaxis()->GetXmax());					      
  gYPrimaryProtonsPercentage->Divide(gYPrimaryProtons,
				     gYTotalProtons,1.,1.0);
  SetError(gYPrimaryProtonsPercentage,gYTotalProtons);
  gYPrimaryProtonsPercentage->Scale(100.);
  gYPrimaryProtonsPercentage->SetMarkerStyle(kFullCircle);

  TH1D *gYSecondaryProtonsPercentage = new TH1D("gYSecondaryProtonsPercentage",
						"",
						hSecondaryProtons->GetXaxis()->GetNbins(),
						hSecondaryProtons->GetXaxis()->GetXmin(),
						hSecondaryProtons->GetXaxis()->GetXmax());					      
  gYSecondaryProtonsPercentage->Divide(gYSecondaryProtons,
				       gYTotalProtons,1.,1.0);
  SetError(gYSecondaryProtonsPercentage,gYTotalProtons);
  gYSecondaryProtonsPercentage->Scale(100.);
  gYSecondaryProtonsPercentage->SetMarkerStyle(kOpenCircle);

  //Antiprotons
  TH1D *gYPrimaryAntiProtons = (TH1D *)hPrimaryAntiProtons->ProjectionX("gYPrimaryAntiProtons",0,hPrimaryAntiProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYSecondaryAntiProtons = (TH1D *)hSecondaryAntiProtons->ProjectionX("gYSecondaryAntiProtons",0,hSecondaryAntiProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYTotalAntiProtons = (TH1D *)hPrimaryAntiProtons->ProjectionX("gYTotalAntiProtons",0,hPrimaryAntiProtons->GetXaxis()->GetNbins(),"e");
  gYTotalAntiProtons->Add(gYSecondaryAntiProtons);
  
  TH1D *gYPrimaryAntiProtonsPercentage = new TH1D("gYPrimaryAntiProtonsPercentage",
						  "",
						  hPrimaryAntiProtons->GetXaxis()->GetNbins(),
						  hPrimaryAntiProtons->GetXaxis()->GetXmin(),
						  hPrimaryAntiProtons->GetXaxis()->GetXmax());					      
  gYPrimaryAntiProtonsPercentage->Divide(gYPrimaryAntiProtons,
					 gYTotalAntiProtons,1.,1.0);
  SetError(gYPrimaryAntiProtonsPercentage,gYTotalAntiProtons);
  gYPrimaryAntiProtonsPercentage->Scale(100.);
  gYPrimaryAntiProtonsPercentage->SetMarkerStyle(kFullCircle);
  
  TH1D *gYSecondaryAntiProtonsPercentage = new TH1D("gYSecondaryAntiProtonsPercentage",
						    "",
						    hSecondaryAntiProtons->GetXaxis()->GetNbins(),
						    hSecondaryAntiProtons->GetXaxis()->GetXmin(),
						    hSecondaryAntiProtons->GetXaxis()->GetXmax());					      
  gYSecondaryAntiProtonsPercentage->Divide(gYSecondaryAntiProtons,
					   gYTotalAntiProtons,1.,1.0);
  SetError(gYSecondaryAntiProtonsPercentage,gYTotalAntiProtons);
  gYSecondaryAntiProtonsPercentage->Scale(100.);
  gYSecondaryAntiProtonsPercentage->SetMarkerStyle(kOpenCircle);
  
  
  TH2F *hEmptyY = new TH2F("hEmptyCompositionY","",
			   100,-1.2,1.2,100,-10.0,130); 
  hEmptyY->SetStats(kFALSE); 
  hEmptyY->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyY->GetYaxis()->SetTitleOffset(1.3);
  hEmptyY->GetXaxis()->SetTitle("y");

  TCanvas *c7 = new TCanvas("c7","(Anti)Proton contamination vs y",
			    150,150,700,400);
  c7->SetFillColor(10); c7->GetFrame()->SetFillColor(10); 
  c7->SetHighLightColor(10); c7->Divide(2,1);

  c7->cd(1)->SetBottomMargin(0.15); 
  c7->cd(1)->SetLeftMargin(0.15); 
  c7->cd(1)->SetGridx(); c7->cd(1)->SetGridy();
  hEmptyY->SetTitle("Protons");
  hEmptyY->DrawCopy();
  gYPrimaryProtonsPercentage->DrawCopy("ESAME");
  gYSecondaryProtonsPercentage->DrawCopy("ESAME");

  DrawMarker(0, 55, kFullCircle, 1.2, 1);
  t1->DrawLatex(0.1,53,"Primaries");
  DrawMarker(0, 45, kOpenCircle, 1.2, 1);
  t1->DrawLatex(0.1,43,"Secondaries");

  c7->cd(2)->SetBottomMargin(0.15); 
  c7->cd(2)->SetLeftMargin(0.15); 
  c7->cd(2)->SetGridx(); c7->cd(2)->SetGridy();
  hEmptyY->SetTitle("Antiprotons");
  hEmptyY->DrawCopy();
  gYPrimaryAntiProtonsPercentage->DrawCopy("ESAME");
  gYSecondaryAntiProtonsPercentage->DrawCopy("ESAME");

  DrawMarker(0, 55, kFullCircle, 1.2, 1);
  t1->DrawLatex(0.1,53,"Primaries");
  DrawMarker(0, 45, kOpenCircle, 1.2, 1);
  t1->DrawLatex(0.1,41,"Secondaries");

  c7->SaveAs("Contamination-Protons-Rapidity.gif");

  //pT dependence
  //Protons
  TH1D *gPtPrimaryProtons = (TH1D *)hPrimaryProtons->ProjectionY("gPtPrimaryProtons",0,hPrimaryProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtSecondaryProtons = (TH1D *)hSecondaryProtons->ProjectionY("gPtSecondaryProtons",0,hSecondaryProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtTotalProtons = (TH1D *)hPrimaryProtons->ProjectionY("gPtTotalProtons",0,hPrimaryProtons->GetYaxis()->GetNbins(),"e");
  gPtTotalProtons->Add(gPtSecondaryProtons);

  TH1D *gPtPrimaryProtonsPercentage = new TH1D("gPtPrimaryProtonsPercentage",
					       "",
					       hPrimaryProtons->GetYaxis()->GetNbins(),
					       hPrimaryProtons->GetYaxis()->GetXmin(),
					       hPrimaryProtons->GetYaxis()->GetXmax());					      
  gPtPrimaryProtonsPercentage->Divide(gPtPrimaryProtons,
				      gPtTotalProtons,1.,1.0);
  SetError(gPtPrimaryProtonsPercentage,gPtTotalProtons);
  gPtPrimaryProtonsPercentage->Scale(100.);
  gPtPrimaryProtonsPercentage->SetMarkerStyle(kFullCircle);
  
  TH1D *gPtSecondaryProtonsPercentage = new TH1D("gPtSecondaryProtonsPercentage",
						 "",
						 hSecondaryProtons->GetYaxis()->GetNbins(),
						 hSecondaryProtons->GetYaxis()->GetXmin(),
						 hSecondaryProtons->GetYaxis()->GetXmax());					      
  gPtSecondaryProtonsPercentage->Divide(gPtSecondaryProtons,
					gPtTotalProtons,1.,1.0);
  SetError(gPtSecondaryProtonsPercentage,gPtTotalProtons);
  gPtSecondaryProtonsPercentage->Scale(100.);
  gPtSecondaryProtonsPercentage->SetMarkerStyle(kOpenCircle);


  //Antiprotons
  TH1D *gPtPrimaryAntiProtons = (TH1D *)hPrimaryAntiProtons->ProjectionY("gPtPrimaryAntiProtons",0,hPrimaryAntiProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtSecondaryAntiProtons = (TH1D *)hSecondaryAntiProtons->ProjectionY("gPtSecondaryAntiProtons",0,hSecondaryAntiProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtTotalAntiProtons = (TH1D *)hPrimaryAntiProtons->ProjectionY("gPtTotalAntiProtons",0,hPrimaryAntiProtons->GetYaxis()->GetNbins(),"e");
  gPtTotalAntiProtons->Add(gPtSecondaryAntiProtons);
  
  TH1D *gPtPrimaryAntiProtonsPercentage = new TH1D("gPtPrimaryAntiProtonsPercentage",
						   "",
						   hPrimaryAntiProtons->GetYaxis()->GetNbins(),
						   hPrimaryAntiProtons->GetYaxis()->GetXmin(),
						   hPrimaryAntiProtons->GetYaxis()->GetXmax());					      
  gPtPrimaryAntiProtonsPercentage->Divide(gPtPrimaryAntiProtons,
					  gPtTotalAntiProtons,1.,1.0);
  SetError(gPtPrimaryAntiProtonsPercentage,gPtTotalAntiProtons);
  gPtPrimaryAntiProtonsPercentage->Scale(100.);
  gPtPrimaryAntiProtonsPercentage->SetMarkerStyle(kFullCircle);
  
  TH1D *gPtSecondaryAntiProtonsPercentage = new TH1D("gPtSecondaryAntiProtonsPercentage",
						     "",
						     hSecondaryAntiProtons->GetYaxis()->GetNbins(),
						     hSecondaryAntiProtons->GetYaxis()->GetXmin(),
						     hSecondaryAntiProtons->GetYaxis()->GetXmax());					      
  gPtSecondaryAntiProtonsPercentage->Divide(gPtSecondaryAntiProtons,
					    gPtTotalAntiProtons,1.,1.0);
  SetError(gPtSecondaryAntiProtonsPercentage,gPtTotalAntiProtons);
  gPtSecondaryAntiProtonsPercentage->Scale(100.);
  gPtSecondaryAntiProtonsPercentage->SetMarkerStyle(kOpenCircle);
  
  TH2F *hEmptyPt = new TH2F("hEmptyCompositionPt","",
			   100,0.0,4.0,100,-10.0,130); 
  hEmptyPt->SetStats(kFALSE); 
  hEmptyPt->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyPt->GetYaxis()->SetTitleOffset(1.3);
  hEmptyPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");

  TCanvas *c8 = new TCanvas("c8","(Anti)Proton comtamination vs pT",
			    200,200,700,400);
  c8->SetFillColor(10); c8->GetFrame()->SetFillColor(10); 
  c8->SetHighLightColor(10); c8->Divide(2,1);

  c8->cd(1)->SetBottomMargin(0.15); 
  c8->cd(1)->SetLeftMargin(0.15); 
  c8->cd(1)->SetGridx(); c8->cd(1)->SetGridy();
  hEmptyPt->SetTitle("Protons");
  hEmptyPt->GetXaxis()->SetRangeUser(gPtPrimaryProtonsPercentage->GetXaxis()->GetXmin()-0.2,gPtPrimaryProtonsPercentage->GetXaxis()->GetXmax()+0.2);
  hEmptyPt->DrawCopy();
  gPtPrimaryProtonsPercentage->DrawCopy("ESAME");
  gPtSecondaryProtonsPercentage->DrawCopy("ESAME");

  DrawMarker(0.5, 55, kFullCircle, 1.2, 1);
  t1->DrawLatex(0.6,53,"Primaries");
  DrawMarker(0.5, 45, kOpenCircle, 1.2, 1);
  t1->DrawLatex(0.6,41,"Secondaries");

  c8->cd(2)->SetBottomMargin(0.15); 
  c8->cd(2)->SetLeftMargin(0.15); 
  c8->cd(2)->SetGridx(); c8->cd(2)->SetGridy();
  hEmptyPt->SetTitle("Antirotons");
  hEmptyPt->DrawCopy();
  gPtPrimaryAntiProtonsPercentage->DrawCopy("ESAME");
  gPtSecondaryAntiProtonsPercentage->DrawCopy("ESAME");

  DrawMarker(2.0, 55, kFullCircle, 1.2, 1);
  t1->DrawLatex(2.1,53,"Primaries");
  DrawMarker(2.0, 45, kOpenCircle, 1.2, 1);
  t1->DrawLatex(2.1,43,"Secondaries");

  c8->SaveAs("Contamination-Protons-Pt.gif");

  TFile *fout = TFile::Open("Contamination.root","recreate");
  gYPrimaryProtonsPercentage->Write();
  gYSecondaryProtonsPercentage->Write();
  gPtPrimaryProtonsPercentage->Write();
  gPtSecondaryProtonsPercentage->Write();
  gYPrimaryAntiProtonsPercentage->Write();
  gYSecondaryAntiProtonsPercentage->Write();
  gPtPrimaryAntiProtonsPercentage->Write();
  gPtSecondaryAntiProtonsPercentage->Write();
  fout->Close();
}

//________________________________________________//
void DrawCutEfficiency(TList *inputList) {
  //loops over the list entries and
  //draws the rapidity and pT dependence
  //of the percentage of primary and secondary
  //protons and antiprotons after the track cuts
  cout<<"Extracting the entries for the histograms in the list: "<<
    inputList->GetName()<<"..."<<endl;

  TLatex *t1 = new TLatex();
  t1->SetTextSize(0.04);

  TH2F *hPrimaryESDProtons = (TH2F *)inputList->At(0);
  TH2F *hPrimaryESDAntiProtons = (TH2F *)inputList->At(4);
  TH2F *hPrimaryMCProtons = (TH2F *)inputList->At(8);
  TH2F *hPrimaryMCAntiProtons = (TH2F *)inputList->At(9);

  //rapidity dependence
  //Protons
  TH1D *gYPrimaryESDProtons = (TH1D *)hPrimaryESDProtons->ProjectionX("gYPrimaryESDProtons",0,hPrimaryESDProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYPrimaryMCProtons = (TH1D *)hPrimaryMCProtons->ProjectionX("gYPrimaryMCProtons",0,hPrimaryMCProtons->GetXaxis()->GetNbins(),"e");
  gYPrimaryESDProtons->Divide(gYPrimaryMCProtons);
  SetError(gYPrimaryESDProtons,gYPrimaryMCProtons);
  gYPrimaryESDProtons->Scale(100.);
  gYPrimaryESDProtons->SetMarkerStyle(kFullCircle);

  //Antiprotons
  TH1D *gYPrimaryESDAntiProtons = (TH1D *)hPrimaryESDAntiProtons->ProjectionX("gYPrimaryESDAntiProtons",0,hPrimaryESDAntiProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYPrimaryMCAntiProtons = (TH1D *)hPrimaryMCAntiProtons->ProjectionX("gYPrimaryMCAntiProtons",0,hPrimaryMCAntiProtons->GetXaxis()->GetNbins(),"e");
  gYPrimaryESDAntiProtons->Divide(gYPrimaryMCAntiProtons);
  SetError(gYPrimaryESDAntiProtons,gYPrimaryMCAntiProtons);
  gYPrimaryESDAntiProtons->Scale(100.);
  gYPrimaryESDAntiProtons->SetMarkerStyle(kFullCircle);
  
  TH2F *hEmptyY = new TH2F("hEmptyEfficiencyY","",
			   100,-1.2,1.2,100,-10.0,130); 
  hEmptyY->SetStats(kFALSE); 
  hEmptyY->GetYaxis()->SetTitle("#epsilon [%]");
  hEmptyY->GetYaxis()->SetTitleOffset(1.3);
  hEmptyY->GetXaxis()->SetTitle("y");

  TCanvas *c10 = new TCanvas("c10","(Anti)Proton cut efficiency vs y",
			    250,250,700,400);
  c10->SetFillColor(10); c10->GetFrame()->SetFillColor(10); 
  c10->SetHighLightColor(10); c10->Divide(2,1);

  c10->cd(1)->SetBottomMargin(0.15); 
  c10->cd(1)->SetLeftMargin(0.15); 
  c10->cd(1)->SetGridx(); c10->cd(1)->SetGridy();
  hEmptyY->SetTitle("Protons");
  hEmptyY->GetXaxis()->SetRangeUser(gYPrimaryESDAntiProtons->GetXaxis()->GetXmin()-0.2,gYPrimaryESDAntiProtons->GetXaxis()->GetXmax()+0.2);
  hEmptyY->DrawCopy();
  gYPrimaryESDProtons->DrawCopy("ESAME");

  c10->cd(2)->SetBottomMargin(0.15); 
  c10->cd(2)->SetLeftMargin(0.15); 
  c10->cd(2)->SetGridx(); c10->cd(2)->SetGridy();
  hEmptyY->SetTitle("Antiprotons");
  hEmptyY->DrawCopy();
  gYPrimaryESDAntiProtons->DrawCopy("ESAME");

  c10->SaveAs("CutEfficiency-Protons-Rapidity.gif");

  //pT dependence
  //Protons
  TH1D *gPtPrimaryESDProtons = (TH1D *)hPrimaryESDProtons->ProjectionY("gPtPrimaryESDProtons",0,hPrimaryESDProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtPrimaryMCProtons = (TH1D *)hPrimaryMCProtons->ProjectionY("gPtPrimaryMCProtons",0,hPrimaryMCProtons->GetYaxis()->GetNbins(),"e");
  gPtPrimaryESDProtons->Divide(gPtPrimaryMCProtons);
  SetError(gPtPrimaryESDProtons,gPtPrimaryMCProtons);
  gPtPrimaryESDProtons->Scale(100.);
  gPtPrimaryESDProtons->SetMarkerStyle(kFullCircle);

  //Antiprotons
  TH1D *gPtPrimaryESDAntiProtons = (TH1D *)hPrimaryESDAntiProtons->ProjectionY("gPtPrimaryESDAntiProtons",0,hPrimaryESDAntiProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtPrimaryMCAntiProtons = (TH1D *)hPrimaryMCAntiProtons->ProjectionY("gPtPrimaryMCAntiProtons",0,hPrimaryMCAntiProtons->GetYaxis()->GetNbins(),"e");
  gPtPrimaryESDAntiProtons->Divide(gPtPrimaryMCAntiProtons);
  SetError(gPtPrimaryESDAntiProtons,gPtPrimaryMCAntiProtons);
  gPtPrimaryESDAntiProtons->Scale(100.);
  gPtPrimaryESDAntiProtons->SetMarkerStyle(kFullCircle);

  TH2F *hEmptyPt = new TH2F("hEmptyEfficiencyPt","",
			   100,0.0,4.0,100,-10.0,130); 
  hEmptyPt->SetStats(kFALSE); 
  hEmptyPt->GetYaxis()->SetTitle("#epsilon [%]");
  hEmptyPt->GetYaxis()->SetTitleOffset(1.3);
  hEmptyPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");

  TCanvas *c11 = new TCanvas("c11","(Anti)Proton cut efficiency vs pT",
			    300,300,700,400);
  c11->SetFillColor(10); c11->GetFrame()->SetFillColor(10); 
  c11->SetHighLightColor(10); c11->Divide(2,1);

  c11->cd(1)->SetBottomMargin(0.15); 
  c11->cd(1)->SetLeftMargin(0.15); 
  c11->cd(1)->SetGridx(); c11->cd(1)->SetGridy();
  hEmptyPt->SetTitle("Protons");
  hEmptyPt->GetXaxis()->SetRangeUser(gPtPrimaryESDAntiProtons->GetXaxis()->GetXmin()-0.2,gPtPrimaryESDAntiProtons->GetXaxis()->GetXmax()+0.2);
  hEmptyPt->DrawCopy();
  gPtPrimaryESDProtons->DrawCopy("ESAME");

  c11->cd(2)->SetBottomMargin(0.15); 
  c11->cd(2)->SetLeftMargin(0.15); 
  c11->cd(2)->SetGridx(); c11->cd(2)->SetGridy();
  hEmptyPt->SetTitle("Antirotons");
  hEmptyPt->DrawCopy();
  gPtPrimaryESDAntiProtons->DrawCopy("ESAME");

  c11->SaveAs("CutEfficiency-Protons-Pt.gif");

  TFile *fout = TFile::Open("CutEfficiency.root","recreate");
  gYPrimaryESDProtons->Write();
  gYPrimaryESDAntiProtons->Write();
  gPtPrimaryESDProtons->Write();
  gPtPrimaryESDAntiProtons->Write();
  fout->Close();
}

//________________________________________________//
void GetQAEntries(TList *inputList, Double_t *entries) {
  //loops over the list entries
  //extracts the entries for each histogram
  //cout<<"Extracting the entries for the histograms in the list: "<<
  //inputList->GetName()<<"..."<<endl;

  for(Int_t i = 0; i < inputList->GetEntries(); i++) {
    TH1F *gHist = (TH1F *)inputList->At(i);
    entries[i] = gHist->GetEntries();
    cout<<"Histogram: "<<gHist->GetName()<<
      " - Entries: "<<entries[i]<<endl;
    gHist = 0;
  }
}

//________________________________________________//
Double_t GetPercentage(Double_t nPassEntries,
		       Double_t nRejectEntries) {
  //returns the percentage of tracks that were rejected by a cut
  Int_t nTotalEntries = nPassEntries + nRejectEntries;

  if(nTotalEntries == 0)
    return -1;

  return 100.*nRejectEntries/nTotalEntries;
}

//________________________________________//
void drawMCQA(TList *listPDG, TList *listMCProcesses) {
  //Function to display the composition of secondary (anti)protons
  //The histogram shows the percentage of secondary (anti)protons
  //originating from each particle species.
  //The box summarizes the MC process that gave these secondary (anti)protons
  TDatabasePDG *db = TDatabasePDG::Instance();
  TParticlePDG *p = 0x0;
  
  TH3F *gHistYPtPDGProtons = (TH3F *)listPDG->At(0);
  TH3F *gHistYPtPDGAntiProtons = (TH3F *)listPDG->At(1);
  readProcesses(listMCProcesses);
  Double_t nParticleCompositionProtonY[100], nParticleCompositionProtonPt[100];
  Double_t nParticleCompositionAntiProtonY[100], nParticleCompositionAntiProtonPt[100];
  Double_t gY[100], gPt[100];
  for(Int_t iBins = 0; iBins < 100; iBins++) {
    nParticleCompositionProtonY[iBins] = 0;
    nParticleCompositionProtonPt[iBins] = 0;
    nParticleCompositionAntiProtonY[iBins] = 0;
    nParticleCompositionAntiProtonPt[iBins] = 0;
    gY[iBins] = 0;
    gPt[iBins] = 0;
  }
  
  TGraph *gParticleProtonY[14];
  TGraph *gParticleProtonPt[14];
  TGraph *gParticleAntiProtonY[14];
  TGraph *gParticleAntiProtonPt[14];
  for(Int_t iParticle = 0; iParticle < 14; iParticle++) {
    GetComposition(iParticle,
		   gHistYPtPDGProtons,
		   nParticleCompositionProtonY,
		   gY, nParticleCompositionProtonPt, gPt);
    gParticleProtonY[iParticle] = new TGraph(gHistYPtPDGProtons->GetNbinsX(),
				       gY,nParticleCompositionProtonY);
    gParticleProtonY[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleProtonY[iParticle]->SetMarkerSize(1.2);

    gParticleProtonPt[iParticle] = new TGraph(gHistYPtPDGProtons->GetNbinsY(),
					gPt,nParticleCompositionProtonPt);
    gParticleProtonPt[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleProtonPt[iParticle]->SetMarkerSize(1.2);

    GetComposition(iParticle,
		   gHistYPtPDGAntiProtons,
		   nParticleCompositionAntiProtonY,
		   gY, nParticleCompositionAntiProtonPt, gPt);
    gParticleAntiProtonY[iParticle] = new TGraph(gHistYPtPDGAntiProtons->GetNbinsX(),
				       gY,nParticleCompositionAntiProtonY);
    gParticleAntiProtonY[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleAntiProtonY[iParticle]->SetMarkerSize(1.2);

    gParticleAntiProtonPt[iParticle] = new TGraph(gHistYPtPDGAntiProtons->GetNbinsY(),
					gPt,nParticleCompositionAntiProtonPt);
    gParticleAntiProtonPt[iParticle]->SetMarkerStyle(iParticle+20);
    gParticleAntiProtonPt[iParticle]->SetMarkerSize(1.2);
  }

  //_________________________________________________________//
  char *fParticleName[14] = {"Primary","K_{L}","#pi","K_{S}","K",
			     "n","p","#Sigma^{-}","#Lambda","#Sigma^{+}",
			     "#Xi^{-}","#Xi^{0}","#Omega^{-}"};
  TLatex *t1 = new TLatex();
  t1->SetTextSize(0.04);

  TH2F *hEmptyY = new TH2F("hEmptyY","",100,-1.2,1.2,100,0,120); 
  hEmptyY->SetStats(kFALSE); 
  hEmptyY->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyY->GetXaxis()->SetTitle("y");

  TCanvas *c3 = new TCanvas("c3","MC secondary composition vs y - Protons",
			    350,350,700,400);
  c3->SetFillColor(10); c3->GetFrame()->SetFillColor(10);
  c3->SetHighLightColor(10); c3->SetBottomMargin(0.15);
  c3->SetGridx(); c3->SetGridy();
  hEmptyY->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    //if((iParticle == 0)||(iParticle == 2)||(iParticle == 5)||(iParticle == 6)||(iParticle == 8))
    if((iParticle == 0)||(iParticle == 2)||(iParticle == 6)||(iParticle == 8))
      gParticleProtonY[iParticle]->Draw("P");
    if(iParticle < 5) {
      DrawMarker(-1.1, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(-1.0,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(0.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
  }

  TCanvas *c5 = new TCanvas("c5","MC secondary composition vs y - antiProtons",
			    400,400,700,400);
  c5->SetFillColor(10); c5->GetFrame()->SetFillColor(10);
  c5->SetHighLightColor(10); c5->SetBottomMargin(0.15);
  c5->SetGridx(); c5->SetGridy();
  hEmptyY->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if((iParticle == 0)||(iParticle == 6)||(iParticle == 8))
      gParticleAntiProtonY[iParticle]->Draw("P");
    if(iParticle < 5) {
      DrawMarker(-1.1, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(-1.0,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(0.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
  }

  TH2F *hEmptyPt = new TH2F("hEmptyPt","",100,0.0,4.0,100,0,120); 
  hEmptyPt->SetStats(kFALSE); 
  hEmptyPt->GetYaxis()->SetTitle("Particle composition [%]");
  hEmptyPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");

  TCanvas *c4 = new TCanvas("c4","MC secondary composition vs pT - Protons",
			    450,450,700,400);
  c4->SetFillColor(10); c4->GetFrame()->SetFillColor(10);
  c4->SetHighLightColor(10); c4->SetBottomMargin(0.15);
  c4->SetGridx(); c4->SetGridy();
  hEmptyPt->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if(iParticle < 5) {
      DrawMarker(0.2, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(2.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(2.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
    if((iParticle == 0)||(iParticle == 2)||(iParticle == 6)||(iParticle == 8))
      gParticleProtonPt[iParticle]->Draw("P");
  }

  TCanvas *c6 = new TCanvas("c6",
			    "MC secondary composition vs pT - AntiProtons",
			    500,500,700,400);
  c6->SetFillColor(10); c6->GetFrame()->SetFillColor(10);
  c6->SetHighLightColor(10); c6->SetBottomMargin(0.15);
  c6->SetGridx(); c6->SetGridy();
  hEmptyPt->DrawCopy();
  for(Int_t iParticle = 0; iParticle < 10; iParticle++) {
    if(iParticle < 5) {
      DrawMarker(0.2, 115-5*iParticle, 20+iParticle, 1.2, 1);
      t1->DrawLatex(0.3,113-5*iParticle,fParticleName[iParticle]);
    }
    else {
      DrawMarker(2.2, 115-5*(iParticle-5), 20+iParticle, 1.2, 1);
      t1->DrawLatex(2.3,113-5*(iParticle-5),fParticleName[iParticle]);
    }
    if((iParticle == 0)||(iParticle == 6)||(iParticle == 8))
      gParticleAntiProtonPt[iParticle]->Draw("P");
  }
}

//________________________________________//
void GetComposition(Int_t iSpecies,
		    TH3F *gHist, 
		    Double_t *nParticleCompositionY,
		    Double_t *gY,
		    Double_t *nParticleCompositionPt,
		    Double_t *gPt) {
  //Returns the pT and y dependence of the MC composition
  Double_t ymin = gHist->GetXaxis()->GetXmin();
  Double_t ymax = gHist->GetXaxis()->GetXmax();
  Double_t nybins = gHist->GetNbinsX();
  Double_t ptmin = gHist->GetYaxis()->GetXmin();
  Double_t ptmax = gHist->GetYaxis()->GetXmax();
  Double_t nptbins = gHist->GetNbinsY();
  Double_t nTotalY[100], nTotalPt[100];
  for(Int_t iBins = 0; iBins < 100; iBins++) {
    nParticleCompositionY[iBins] = 0;
    nParticleCompositionPt[iBins] = 0;
    nTotalY[iBins] = 0.0;
    nTotalPt[iBins] = 0.0;
  }

  //rapidity dependence
  //cout<<"Ymin: "<<ymin<<" - Ymax: "<<ymax<<" - Ybins: "<<nybins<<endl;
  for(Int_t iXbins = 1; iXbins <= gHist->GetNbinsX(); iXbins++) {
    for(Int_t iZbins = 1; iZbins <= gHist->GetNbinsZ(); iZbins++) {
      for(Int_t iYbins = 1; iYbins <= gHist->GetNbinsY(); iYbins++) {
	nTotalY[iXbins-1] += gHist->GetBinContent(iXbins,iYbins,iZbins);
      }
    }
  }

  for(Int_t iXbins = 1; iXbins <= gHist->GetNbinsX(); iXbins++) {
    for(Int_t iYbins = 1; iYbins <= gHist->GetNbinsY(); iYbins++) {
      if(nTotalY[iXbins-1] > 0)
	nParticleCompositionY[iXbins-1] += 100.*gHist->GetBinContent(iXbins,iYbins,iSpecies+1)/nTotalY[iXbins-1];
      //if(nParticleCompositionY[iXbins-1] == 0) 
      //nParticleCompositionY[iXbins-1] = -10.0;
    }//pt loop
    gY[iXbins-1] = ymin + (iXbins-1)*(ymax - ymin)/nybins + 0.5*(ymax - ymin)/nybins;
    //cout<<"y: "<<gY[iXbins-1]<<
    //" - test: "<<ymin + (iXbins-1)*(ymax - ymin)/nybins + 0.5*(ymax - ymin)/nybins<<
    //" - Number of protons: "<<nY[iXbins-1]<<
    //" - Total: "<<nTotalY[iXbins-1]<<
    //" - Percentage: "<<nParticleCompositionY[iXbins-1]<<endl;
  }//y loop

  //pt dependence
  //cout<<"Ptmin: "<<ptmin<<" - Ptmax: "<<ptmax<<" - Ptbins: "<<nptbins<<endl;
  for(Int_t iYbins = 1; iYbins <= gHist->GetNbinsY(); iYbins++) {
    for(Int_t iZbins = 1; iZbins <= gHist->GetNbinsZ(); iZbins++) {
      for(Int_t iXbins = 1; iXbins <= gHist->GetNbinsX(); iXbins++) {
	nTotalPt[iYbins-1] += gHist->GetBinContent(iXbins,iYbins,iZbins);
      }
    }
  }

  for(Int_t iYbins = 1; iYbins <= gHist->GetNbinsY(); iYbins++) {
    for(Int_t iXbins = 1; iXbins <= gHist->GetNbinsX(); iXbins++) {
      if(nTotalPt[iYbins-1] > 0)
	nParticleCompositionPt[iYbins-1] += 100.*gHist->GetBinContent(iXbins,iYbins,iSpecies+1)/nTotalPt[iYbins-1];
      //if(nParticleCompositionPt[iYbins-1] == 0) 
      //nParticleCompositionPt[iYbins-1] = -10.0;
    }//pt loop
    gPt[iYbins-1] = ptmin + (iYbins-1)*(ptmax - ptmin)/nptbins + 0.5*(ptmax - ptmin)/nptbins;
    //cout<<"Pt: "<<gPt[iYbins-1]<<
    //" - test: "<<ptmin + (iYbins-1)*(ptmax - ptmin)/nptbins + 0.5*(ptmax - ptmin)/nptbins<<
    //" - Number of protons: "<<nY[iXbins-1]<<
    //" - Total: "<<nTotalPt[iYbins-1]<<
    //" - Percentage: "<<nParticleCompositionPt[iYbins-1]<<endl;
  }//pt loop
}

//________________________________________//
void readProcesses(TList *list) {
  char *fParticleProtonName[12] = {"K_{L}","#pi","K_{S}","K",
				   "n","p","#Sigma^{-}","#Lambda","#Sigma^{+}",
				   "#Xi^{-}","#Xi^{0}","#Omega^{-}"};
  char *fParticleAntiProtonName[8] = {"K_{L}","#pi","K_{S}","K",
				      "n","p","#Lambda","#Sigma^{+}"};
  Int_t iProtonCounter = 0, iAntiProtonCounter = 0;
  TH1F *gMCProcesses;
  for(Int_t iEntry = 0; iEntry < list->GetEntries(); iEntry++) {
    gMCProcesses = (TH1F *)list->At(iEntry);
    TString histName = gMCProcesses->GetName();
    if(histName.Contains("gHistProtons")) {
      cout<<"Protons coming from "<<fParticleProtonName[iProtonCounter]<<endl;

      iProtonCounter += 1;
    }
    if(histName.Contains("gHistAntiProtons")) {
      cout<<"Antiprotons coming from "<<fParticleAntiProtonName[iAntiProtonCounter]<<endl;

      iAntiProtonCounter += 1;
    }
    for(Int_t iBin = 1; iBin < gMCProcesses->GetNbinsX(); iBin++) {
      Double_t binContent = gMCProcesses->GetBinContent(iBin);
      if(binContent > 0) {
	Int_t processId = gMCProcesses->GetBinCenter(iBin);
	cout<<"\t Process ID: "<<processId<<" - "<<
	  gMCProcessName[processId]<<" - "<<
	  100.*binContent/gMCProcesses->GetEntries()<<"%"<<endl;
      }
    }
  }
}

//________________________________________________//
void SetError(TH1D *hEff, TH1D *hGen) {
  for(Int_t iBin = 1; iBin <= hEff->GetNbinsX(); iBin++) {
    Double_t error = 0.0;
    if((hEff->GetBinContent(iBin) <= 1  .)&&(hGen->GetBinContent(iBin) != 0))
      error = TMath::Sqrt(hEff->GetBinContent(iBin)*(1  . - hEff->GetBinContent(iBin))/hGen->GetBinContent(iBin));
    hEff->SetBinError(iBin,error);
  }
}

//________________________________________//
void drawEfficiency(TList *list) {
  //Function to display the reconstruction and PID efficiencies
  //for protons and antiprotons vs y and pT

  TH2F *hEmpty = new TH2F("hEmptyReconstructionEfficiency","",
			   100,-1.2,3.5,100,-10.0,130); 
  hEmpty->SetStats(kFALSE); 
  hEmpty->GetYaxis()->SetTitle("#epsilon [%]");
  hEmpty->GetYaxis()->SetTitleOffset(1.3);

  //Reconstruction efficiency
  TH2D *gHistPrimariesMCYPtProtons = (TH2D *)list->At(0);
  TH2D *gHistPrimariesMCYPtAntiProtons = (TH2D *)list->At(1);
  TH2D *gHistMCYPtProtonsFromWeak = (TH2D *)list->At(2);
  TH2D *gHistMCYPtAntiProtonsFromWeak = (TH2D *)list->At(3);
  TH2D *gHistMCYPtProtonsFromHadronic = (TH2D *)list->At(4);
  TH2D *gHistMCYPtAntiProtonsFromHadronic = (TH2D *)list->At(5);
  TH2D *gHistPrimariesESDYPtProtons = (TH2D *)list->At(6);
  TH2D *gHistPrimariesESDYPtAntiProtons = (TH2D *)list->At(7);
  TH2D *gHistESDYPtProtonsFromWeak = (TH2D *)list->At(8);
  TH2D *gHistESDYPtAntiProtonsFromWeak = (TH2D *)list->At(9);
  TH2D *gHistESDYPtProtonsFromHadronic = (TH2D *)list->At(10);
  TH2D *gHistESDYPtAntiProtonsFromHadronic = (TH2D *)list->At(11);

  //rapidity dependence
  TCanvas *c14 = new TCanvas("c14",
			     "(Anti)Proton reconstruction efficiency vs y",
			     600,600,700,400);
  c14->SetFillColor(10); c14->GetFrame()->SetFillColor(10); 
  c14->SetHighLightColor(10); c14->Divide(2,1);

  //Primary Protons
  TH1D *gYPrimariesESDProtons = (TH1D *)gHistPrimariesESDYPtProtons->ProjectionX("gYPrimariesESDProtons",0,gHistPrimariesESDYPtProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYPrimariesMCProtons = (TH1D *)gHistPrimariesMCYPtProtons->ProjectionX("gYPrimariesMCProtons",0,gHistPrimariesMCYPtProtons->GetXaxis()->GetNbins(),"e");
  gYPrimariesESDProtons->Divide(gYPrimariesMCProtons);
  SetError(gYPrimariesESDProtons,gYPrimariesMCProtons);
  gYPrimariesESDProtons->Scale(100.);
  gYPrimariesESDProtons->SetMarkerStyle(kFullCircle);

  //Primary AntiProtons
  TH1D *gYPrimariesESDAntiProtons = (TH1D *)gHistPrimariesESDYPtAntiProtons->ProjectionX("gYPrimariesESDAntiProtons",0,gHistPrimariesESDYPtAntiProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYPrimariesMCAntiProtons = (TH1D *)gHistPrimariesMCYPtAntiProtons->ProjectionX("gYPrimariesMCAntiProtons",0,gHistPrimariesMCYPtProtons->GetXaxis()->GetNbins(),"e");
  gYPrimariesESDAntiProtons->Divide(gYPrimariesMCAntiProtons);
  SetError(gYPrimariesESDAntiProtons,gYPrimariesMCAntiProtons);
  gYPrimariesESDAntiProtons->Scale(100.);
  gYPrimariesESDAntiProtons->SetMarkerStyle(kFullCircle);

  //Protons from weak decays
  TH1D *gYESDProtonsFromWeak = (TH1D *)gHistESDYPtProtonsFromWeak->ProjectionX("gYESDProtonsFromWeak",0,gHistESDYPtProtonsFromWeak->GetXaxis()->GetNbins(),"e");
  TH1D *gYMCProtonsFromWeak = (TH1D *)gHistMCYPtProtonsFromWeak->ProjectionX("gYMCProtonsFromWeak",0,gHistMCYPtProtonsFromWeak->GetXaxis()->GetNbins(),"e");
  gYESDProtonsFromWeak->Divide(gYMCProtonsFromWeak);
  SetError(gYESDProtonsFromWeak,gYMCProtonsFromWeak);
  gYESDProtonsFromWeak->Scale(100.);
  gYESDProtonsFromWeak->SetMarkerStyle(21);
  gYESDProtonsFromWeak->SetMarkerColor(2);

  //AntiProtons from weak decays
  TH1D *gYESDAntiProtonsFromWeak = (TH1D *)gHistESDYPtAntiProtonsFromWeak->ProjectionX("gYESDAntiProtonsFromWeak",0,gHistESDYPtAntiProtonsFromWeak->GetXaxis()->GetNbins(),"e");
  TH1D *gYMCAntiProtonsFromWeak = (TH1D *)gHistMCYPtAntiProtonsFromWeak->ProjectionX("gYMCAntiProtonsFromWeak",0,gHistMCYPtProtonsFromWeak->GetXaxis()->GetNbins(),"e");
  gYESDAntiProtonsFromWeak->Divide(gYMCAntiProtonsFromWeak);
  SetError(gYESDAntiProtonsFromWeak,gYMCAntiProtonsFromWeak);
  gYESDAntiProtonsFromWeak->Scale(100.);
  gYESDAntiProtonsFromWeak->SetMarkerStyle(21);
  gYESDAntiProtonsFromWeak->SetMarkerColor(2);

  //Protons from hadronic interactions
  TH1D *gYESDProtonsFromHadronic = (TH1D *)gHistESDYPtProtonsFromHadronic->ProjectionX("gYESDProtonsFromHadronic",0,gHistESDYPtProtonsFromHadronic->GetXaxis()->GetNbins(),"e");
  TH1D *gYMCProtonsFromHadronic = (TH1D *)gHistMCYPtProtonsFromHadronic->ProjectionX("gYMCProtonsFromHadronic",0,gHistMCYPtProtonsFromHadronic->GetXaxis()->GetNbins(),"e");
  gYESDProtonsFromHadronic->Divide(gYMCProtonsFromHadronic);
  SetError(gYESDProtonsFromHadronic,gYMCProtonsFromHadronic);
  gYESDProtonsFromHadronic->Scale(100.);
  gYESDProtonsFromHadronic->SetMarkerStyle(22);
  gYESDProtonsFromHadronic->SetMarkerColor(3);

  //AntiProtons from weak decays
  TH1D *gYESDAntiProtonsFromHadronic = (TH1D *)gHistESDYPtAntiProtonsFromHadronic->ProjectionX("gYESDAntiProtonsFromHadronic",0,gHistESDYPtAntiProtonsFromHadronic->GetXaxis()->GetNbins(),"e");
  TH1D *gYMCAntiProtonsFromHadronic = (TH1D *)gHistMCYPtAntiProtonsFromHadronic->ProjectionX("gYMCAntiProtonsFromHadronic",0,gHistMCYPtProtonsFromHadronic->GetXaxis()->GetNbins(),"e");
  gYESDAntiProtonsFromHadronic->Divide(gYMCAntiProtonsFromHadronic);
  SetError(gYESDAntiProtonsFromHadronic,gYMCAntiProtonsFromHadronic);
  gYESDAntiProtonsFromHadronic->Scale(100.);
  gYESDAntiProtonsFromHadronic->SetMarkerStyle(22);
  gYESDAntiProtonsFromHadronic->SetMarkerColor(3);

  c14->cd(1)->SetBottomMargin(0.15); 
  c14->cd(1)->SetLeftMargin(0.15); 
  c14->cd(1)->SetGridx(); c14->cd(1)->SetGridy();
  hEmpty->GetXaxis()->SetRangeUser(gYPrimariesESDAntiProtons->GetXaxis()->GetXmin()-0.2,
				   gYPrimariesESDAntiProtons->GetXaxis()->GetXmax()+0.2);
  hEmpty->GetXaxis()->SetTitle(gYPrimariesESDAntiProtons->GetXaxis()->GetTitle());
  hEmpty->SetTitle("Protons");
  hEmpty->DrawCopy();
  gYPrimariesESDProtons->DrawCopy("ESAME");
  gYESDProtonsFromWeak->DrawCopy("ESAME");
  gYESDProtonsFromHadronic->DrawCopy("ESAME");

  c14->cd(2)->SetBottomMargin(0.15); 
  c14->cd(2)->SetLeftMargin(0.15); 
  c14->cd(2)->SetGridx(); c14->cd(2)->SetGridy();
  hEmpty->SetTitle("Antiprotons");
  hEmpty->DrawCopy();
  gYPrimariesESDAntiProtons->DrawCopy("ESAME");
  gYESDAntiProtonsFromWeak->DrawCopy("ESAME");
  gYESDAntiProtonsFromHadronic->DrawCopy("ESAME");
  c14->SaveAs("ReconstructionEfficiency-Protons-Rapidity.gif");

  //pT dependence
  TCanvas *c15 = new TCanvas("c15",
			     "(Anti)Proton reconstruction efficiency vs pT",
			     650,650,700,400);
  c15->SetFillColor(10); c15->GetFrame()->SetFillColor(10); 
  c15->SetHighLightColor(10); c15->Divide(2,1);

  //Primary Protons
  TH1D *gPtPrimariesESDProtons = (TH1D *)gHistPrimariesESDYPtProtons->ProjectionY("gPtPrimariesESDProtons",0,gHistPrimariesESDYPtProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtPrimariesMCProtons = (TH1D *)gHistPrimariesMCYPtProtons->ProjectionY("gPtPrimariesMCProtons",0,gHistPrimariesMCYPtProtons->GetYaxis()->GetNbins(),"e");
  gPtPrimariesESDProtons->Divide(gPtPrimariesMCProtons);
  SetError(gPtPrimariesESDProtons,gPtPrimariesMCProtons);
  gPtPrimariesESDProtons->Scale(100.);
  gPtPrimariesESDProtons->SetMarkerStyle(kFullCircle);

  //Primary AntiProtons
  TH1D *gPtPrimariesESDAntiProtons = (TH1D *)gHistPrimariesESDYPtAntiProtons->ProjectionY("gPtPrimariesESDAntiProtons",0,gHistPrimariesESDYPtAntiProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtPrimariesMCAntiProtons = (TH1D *)gHistPrimariesMCYPtAntiProtons->ProjectionY("gPtPrimariesMCAntiProtons",0,gHistPrimariesMCYPtProtons->GetYaxis()->GetNbins(),"e");
  gPtPrimariesESDAntiProtons->Divide(gPtPrimariesMCAntiProtons);
  SetError(gPtPrimariesESDAntiProtons,gPtPrimariesMCAntiProtons);
  gPtPrimariesESDAntiProtons->Scale(100.);
  gPtPrimariesESDAntiProtons->SetMarkerStyle(kFullCircle);

  //Protons from weak decays
  TH1D *gPtESDProtonsFromWeak = (TH1D *)gHistESDYPtProtonsFromWeak->ProjectionY("gPtESDProtonsFromWeak",0,gHistESDYPtProtonsFromWeak->GetYaxis()->GetNbins(),"e");
  TH1D *gPtMCProtonsFromWeak = (TH1D *)gHistMCYPtProtonsFromWeak->ProjectionY("gPtMCProtonsFromWeak",0,gHistMCYPtProtonsFromWeak->GetYaxis()->GetNbins(),"e");
  gPtESDProtonsFromWeak->Divide(gPtMCProtonsFromWeak);
  SetError(gPtESDProtonsFromWeak,gPtMCProtonsFromWeak);
  gPtESDProtonsFromWeak->Scale(100.);
  gPtESDProtonsFromWeak->SetMarkerStyle(21);
  gPtESDProtonsFromWeak->SetMarkerColor(2);

  //AntiProtons from weak decays
  TH1D *gPtESDAntiProtonsFromWeak = (TH1D *)gHistESDYPtAntiProtonsFromWeak->ProjectionY("gPtESDAntiProtonsFromWeak",0,gHistESDYPtAntiProtonsFromWeak->GetYaxis()->GetNbins(),"e");
  TH1D *gPtMCAntiProtonsFromWeak = (TH1D *)gHistMCYPtAntiProtonsFromWeak->ProjectionY("gPtMCAntiProtonsFromWeak",0,gHistMCYPtProtonsFromWeak->GetYaxis()->GetNbins(),"e");
  gPtESDAntiProtonsFromWeak->Divide(gPtMCAntiProtonsFromWeak);
  SetError(gPtESDAntiProtonsFromWeak,gPtMCAntiProtonsFromWeak);
  gPtESDAntiProtonsFromWeak->Scale(100.);
  gPtESDAntiProtonsFromWeak->SetMarkerStyle(21);
  gPtESDAntiProtonsFromWeak->SetMarkerColor(2);

  //Protons from hadronic interactions
  TH1D *gPtESDProtonsFromHadronic = (TH1D *)gHistESDYPtProtonsFromHadronic->ProjectionY("gPtESDProtonsFromHadronic",0,gHistESDYPtProtonsFromHadronic->GetYaxis()->GetNbins(),"e");
  TH1D *gPtMCProtonsFromHadronic = (TH1D *)gHistMCYPtProtonsFromHadronic->ProjectionY("gPtMCProtonsFromHadronic",0,gHistMCYPtProtonsFromHadronic->GetYaxis()->GetNbins(),"e");
  gPtESDProtonsFromHadronic->Divide(gPtMCProtonsFromHadronic);
  SetError(gPtESDProtonsFromHadronic,gPtMCProtonsFromHadronic);
  gPtESDProtonsFromHadronic->Scale(100.);
  gPtESDProtonsFromHadronic->SetMarkerStyle(22);
  gPtESDProtonsFromHadronic->SetMarkerColor(3);

  //AntiProtons from hadronic interactions
  TH1D *gPtESDAntiProtonsFromHadronic = (TH1D *)gHistESDYPtAntiProtonsFromHadronic->ProjectionY("gPtESDAntiProtonsFromHadronic",0,gHistESDYPtAntiProtonsFromHadronic->GetYaxis()->GetNbins(),"e");
  TH1D *gPtMCAntiProtonsFromHadronic = (TH1D *)gHistMCYPtAntiProtonsFromHadronic->ProjectionY("gPtMCAntiProtonsFromHadronic",0,gHistMCYPtProtonsFromHadronic->GetYaxis()->GetNbins(),"e");
  gPtESDAntiProtonsFromHadronic->Divide(gPtMCAntiProtonsFromHadronic);
  SetError(gPtESDAntiProtonsFromHadronic,gPtMCAntiProtonsFromHadronic);
  gPtESDAntiProtonsFromHadronic->Scale(100.);
  gPtESDAntiProtonsFromHadronic->SetMarkerStyle(22);
  gPtESDAntiProtonsFromHadronic->SetMarkerColor(3);


  c15->cd(1)->SetBottomMargin(0.15); 
  c15->cd(1)->SetLeftMargin(0.15); 
  c15->cd(1)->SetGridx(); c15->cd(1)->SetGridy();
  hEmpty->GetXaxis()->SetRangeUser(gPtPrimariesESDAntiProtons->GetXaxis()->GetXmin()-0.2,
				   gPtPrimariesESDAntiProtons->GetXaxis()->GetXmax()+0.2);
  hEmpty->GetXaxis()->SetTitle("P_{T} [GeV/c]");
  hEmpty->SetTitle("Protons");
  hEmpty->DrawCopy();
  gPtPrimariesESDProtons->DrawCopy("ESAME");
  gPtESDProtonsFromWeak->DrawCopy("ESAME");
  gPtESDProtonsFromHadronic->DrawCopy("ESAME");

  c15->cd(2)->SetBottomMargin(0.15); 
  c15->cd(2)->SetLeftMargin(0.15); 
  c15->cd(2)->SetGridx(); c15->cd(2)->SetGridy();
  hEmpty->SetTitle("Antiprotons");
  hEmpty->DrawCopy();
  gPtPrimariesESDAntiProtons->DrawCopy("ESAME");
  gPtESDAntiProtonsFromWeak->DrawCopy("ESAME");
  gPtESDAntiProtonsFromHadronic->DrawCopy("ESAME");
  c15->SaveAs("ReconstructionEfficiency-Protons-Pt.gif");

  //______________//
  //PID efficiency//
  //______________//
  TH2D *gHistESDInitYPtProtons = (TH2D *)list->At(12);
  TH2D *gHistESDIdYPtProtons = (TH2D *)list->At(13);
  TH2D *gHistESDRecIdYPtProtons = (TH2D *)list->At(14);
  TH2D *gHistESDContamYPtProtons = (TH2D *)list->At(15);

  TCanvas *c16 = new TCanvas("c16",
			     "(Anti)Proton PID efficiency vs y and pT",
			     700,700,700,400);
  c16->SetFillColor(10); c16->GetFrame()->SetFillColor(10); 
  c16->SetHighLightColor(10); c16->Divide(2,1);

  //rapidity dependence
  //protons pid efficiency
  TH1D *gYESDIdProtons = (TH1D *)gHistESDIdYPtProtons->ProjectionX("gYESDIdProtons",0,gHistESDIdYPtProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYESDInitProtons = (TH1D *)gHistESDInitYPtProtons->ProjectionX("gYESDInitProtons",0,gHistESDInitYPtProtons->GetXaxis()->GetNbins(),"e");
  gYESDIdProtons->Divide(gYESDInitProtons);
  SetError(gYESDIdProtons,gYESDInitProtons);
  gYESDIdProtons->Scale(100.);
  gYESDIdProtons->SetMarkerStyle(kFullCircle);

  //protons pid contamination
  TH1D *gYESDContamProtons = (TH1D *)gHistESDContamYPtProtons->ProjectionX("gYESDContamProtons",0,gHistESDContamYPtProtons->GetXaxis()->GetNbins(),"e");
  TH1D *gYESDRecIdProtons = (TH1D *)gHistESDRecIdYPtProtons->ProjectionX("gYESDRecIdProtons",0,gHistESDRecIdYPtProtons->GetXaxis()->GetNbins(),"e");
  gYESDContamProtons->Divide(gYESDRecIdProtons);
  SetError(gYESDContamProtons,gYESDRecIdProtons);
  gYESDContamProtons->Scale(100.);
  gYESDContamProtons->SetMarkerStyle(kOpenCircle);

  c16->cd(1)->SetBottomMargin(0.15); 
  c16->cd(1)->SetLeftMargin(0.15); 
  c16->cd(1)->SetGridx(); c16->cd(1)->SetGridy();
  hEmpty->GetXaxis()->SetRangeUser(-1.0.,1.0);
  hEmpty->GetXaxis()->SetTitle(gYESDContamProtons->GetXaxis()->GetTitle());
  hEmpty->DrawCopy();
  gYESDIdProtons->DrawCopy("ESAME");
  gYESDContamProtons->DrawCopy("ESAME");

  //pT dependence
  //protons pid efficiency
  TH1D *gPtESDIdProtons = (TH1D *)gHistESDIdYPtProtons->ProjectionY("gPtESDIdProtons",0,gHistESDIdYPtProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtESDInitProtons = (TH1D *)gHistESDInitYPtProtons->ProjectionY("gPtESDInitProtons",0,gHistESDInitYPtProtons->GetYaxis()->GetNbins(),"e");
  gPtESDIdProtons->Divide(gPtESDInitProtons);
  SetError(gPtESDIdProtons,gPtESDInitProtons);
  gPtESDIdProtons->Scale(100.);
  gPtESDIdProtons->SetMarkerStyle(kFullCircle);

  //protons pid contamination
  TH1D *gPtESDContamProtons = (TH1D *)gHistESDContamYPtProtons->ProjectionY("gPtESDContamProtons",0,gHistESDContamYPtProtons->GetYaxis()->GetNbins(),"e");
  TH1D *gPtESDRecIdProtons = (TH1D *)gHistESDRecIdYPtProtons->ProjectionY("gPtESDRecIdProtons",0,gHistESDRecIdYPtProtons->GetYaxis()->GetNbins(),"e");
  gPtESDContamProtons->Divide(gPtESDRecIdProtons);
  SetError(gPtESDContamProtons,gPtESDRecIdProtons);
  gPtESDContamProtons->Scale(100.);
  gPtESDContamProtons->SetMarkerStyle(kOpenCircle);

  c16->cd(2)->SetBottomMargin(0.15); 
  c16->cd(2)->SetLeftMargin(0.15); 
  c16->cd(2)->SetGridx(); c16->cd(2)->SetGridy();
  hEmpty->GetXaxis()->SetRangeUser(0.0,1.2);
  hEmpty->GetXaxis()->SetTitle("P_{T} [GeV/c]");
  hEmpty->DrawCopy();
  gPtESDIdProtons->DrawCopy("ESAME");
  gPtESDContamProtons->DrawCopy("ESAME");

  c16->SaveAs("PIDEfficiency-Protons.gif");

  TFile *fout = TFile::Open("Reconstruction-PID-Efficiency.root","recreate");
  gYPrimariesESDProtons->Write();
  gYESDProtonsFromWeak->Write();
  gYESDProtonsFromHadronic->Write();
  gPtPrimariesESDProtons->Write();
  gPtESDProtonsFromWeak->Write();
  gPtESDProtonsFromHadronic->Write();
  gYPrimariesESDAntiProtons->Write();
  gYESDAntiProtonsFromWeak->Write();
  gYESDAntiProtonsFromHadronic->Write();
  gPtPrimariesESDAntiProtons->Write();
  gPtESDAntiProtonsFromWeak->Write();
  gPtESDAntiProtonsFromHadronic->Write();
  gYESDIdProtons->Write();
  gYESDContamProtons->Write();
  gPtESDIdProtons->Write();
  gPtESDContamProtons->Write();
  fout->Close();
}

//________________________________________________//
void DrawMarker(Double_t x, Double_t y, Int_t style, 
		Double_t size, Int_t color) {
  TMarker *m = new TMarker(x,y,style);
  m->SetMarkerSize(size);
  m->SetMarkerColor(color);
  m->Draw();
}

//________________________________________________//
const char * const gMCProcessName[45] = {
  "Primary particle emission",
  "Multiple scattering",
  "Energy loss",
  "Bending in magnetic field",
  "Decay",
  "Lepton pair production",
  "Compton scattering",
  "Photoelectric effect",
  "Bremstrahlung",
  "Delta ray",
  "Positron annihilation",
  "Positron annihilation at rest",
  "Positron annihilation in flight",
  "Hadronic interaction",
  "Nuclear evaporation",
  "Nuclear fission",
  "Nuclear absorbtion",
  "Antiproton annihilation",
  "Antineutron annihilation",
  "Neutron capture",
  "Hadronic elastic",
  "Hadronic incoherent elastic",
  "Hadronic coherent elastic",
  "Hadronic inelastic",
  "Photon inelastic",
  "Muon nuclear interaction",
  "Electron nuclear interaction",
  "Positron nuclear interaction",
  "Time of flight limit",
  "Nuclear photofission",
  "Rayleigh effect",
  "No active process",
  "Energy threshold",
  "Light absorption",
  "Light detection",
  "Light scattering",
  "Maximum allowed step",
  "Cerenkov production",
  "Cerenkov feed back photon",
  "Cerenkov photon reflection",
  "Cerenkov photon refraction",
  "Synchrotron radiation",
  "Scintillation",
  "Transportation",
  "Unknown process"
};