[u/mrichter/AliRoot.git] / PWGPP / VZERO / checkThr.C

void checkThr(const char *filename, const char *what="")
{
  TFile *f = TFile::Open(filename);
  TList *list = (TList*)f->Get("coutput");

  TH1F *h3 = (TH1F*)list->FindObject("fV0Percent");
  Double_t nTotalEvts = (Double_t)h3->Integral(h3->FindBin(0.0),h3->FindBin(91.0))/0.9;

  TH1F *hCent;
  TH1F *hCentAll;
  TH1F *hSemiCent;
  TH1F *hSemiCentAll;
  Double_t nentCent;
  Double_t nentCentall;
  Double_t nentSemiCent;
  Double_t nentSemiCentall;
  TF1 *ffCent;
  TF1 *ffSemiCent;

  hCent = (TH1F*)list->FindObject(Form("fV0Cent%s",what));
  hCentAll = (TH1F*)list->FindObject(Form("fV0Cent%sAll",what));
  hSemiCent = (TH1F*)list->FindObject(Form("fV0SemiCent%s",what));
  hSemiCentAll = (TH1F*)list->FindObject(Form("fV0SemiCent%sAll",what));

  nentCent = hCent->GetEntries();
  Double_t nCentEvts = (Double_t)hCent->Integral(hCent->FindBin(0.0),hCent->FindBin(91.0));
  nentCentall = hCentAll->GetEntries();
  nentSemiCent = hSemiCent->GetEntries();
  Double_t nSemiCentEvts = (Double_t)hSemiCent->Integral(hSemiCent->FindBin(0.0),hSemiCent->FindBin(91.0));
  nentSemiCentall = hSemiCentAll->GetEntries();

  hCent->Sumw2();
  hCent->Divide(hCent,h3,1,1,"B");
  hSemiCent->Sumw2();
  hSemiCent->Divide(hSemiCent,h3,1,1,"B");
  // for(Int_t ibin = 1; ibin <= hSemiCent->GetXaxis()->GetNbins(); ++ibin) {
  //   printf("%f ",hSemiCent->GetBinError(ibin));
  //   if (hSemiCent->GetBinError(ibin) < 1e-6) hSemiCent->SetBinError(ibin,0.01);
  //   printf("%f ",hSemiCent->GetBinError(ibin));
  // }

  ffCent = new TF1("ffCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
  ffCent->SetLineColor(kBlue);
  Double_t par0 = hCent->GetBinCenter(hCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.6 : 0.95));
  ffCent->SetParameter(0,par0);
  ffCent->SetParameter(1,0.5);
  if (strcmp(what,"Tr") == 0)
    ffCent->SetParameter(2,0.6);
  else
    ffCent->SetParLimits(2,0.9999,1.0001);
  hCent->Fit(ffCent,"R+");
  hCent->Fit(ffCent,"R+");
  hCent->Fit(ffCent,"R+");

  ffSemiCent = new TF1("ffSemiCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
  ffSemiCent->SetLineColor(kRed);
  par0 = hSemiCent->GetBinCenter(hSemiCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.01 : 0.95));
  ffSemiCent->SetParameter(0,par0);
  ffSemiCent->SetParameter(1,0.3);
  printf("%d\n",par0);
  if (strcmp(what,"Tr") == 0)
    ffSemiCent->SetParameter(2,0.1);
  else
    ffSemiCent->SetParLimits(2,0.9999,1.0001);
  hSemiCent->Fit(ffSemiCent,"R+");
  hSemiCent->Fit(ffSemiCent,"R+");
  hSemiCent->Fit(ffSemiCent,"R+");

  //  new TCanvas;
  hSemiCent->SetLineColor(kBlue);
  hSemiCent->SetMarkerColor(kBlue);
  hCent->Draw("e");
  hSemiCent->SetLineColor(kRed);
  hSemiCent->SetMarkerColor(kRed);
  hSemiCent->Draw("e same");

  Float_t rawRatio = (strcmp(what,"Tr") == 0) ? (nentSemiCentall/ffSemiCent->GetParameter(2))/(nentCentall/ffCent->GetParameter(2)) : nentSemiCentall/nentCentall;
  Float_t rawRatioErr = rawRatio*TMath::Sqrt(1./nentCentall-1./nentSemiCentall);
  Float_t psRatio = (strcmp(what,"Tr") == 0) ? (nentSemiCent/ffSemiCent->GetParameter(2))/(nentCent/ffCent->GetParameter(2)) : nentSemiCent/nentCent;
  Float_t psRatioErr = psRatio*TMath::Sqrt(1./nentCent-1./nentSemiCent);
  Float_t psScRatio = (strcmp(what,"Tr") == 0) ? (nSemiCentEvts/ffSemiCent->GetParameter(2))/(nCentEvts/ffCent->GetParameter(2)) : nSemiCentEvts/nCentEvts;
  Float_t psScRatioErr = psScRatio*TMath::Sqrt(1./nCentEvts-1./nSemiCentEvts);

  Float_t centrInt = 100.*nCentEvts/ffCent->GetParameter(2)/nTotalEvts;
  Float_t centrIntErr = 100.*TMath::Sqrt(nCentEvts/ffCent->GetParameter(2)/nTotalEvts*(1.-nCentEvts/ffCent->GetParameter(2)/nTotalEvts)/nTotalEvts);
  Float_t centrIntFit = ffCent->GetX(0.5*ffCent->GetParameter(2));
  Float_t centrIntFitErr = ffCent->GetParError(0);

  printf("\n\n\n");

  printf("Purity: Centr      %.3f %.3f   Cuts ( 0.0-%.1f)%%   rate=%.1f Hz\n",
	 nCentEvts/nCentEvts,
	 nentCent/nentCentall,
	 ffCent->GetX(0.99*ffCent->GetParameter(2)),
	 30.*10./(nentCent/nentCentall));
  printf("Purity: Semi-Centr %.3f %.3f   Cuts ( 0.0-%.1f)%%   rate=%.1f Hz\n",
	 ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100),
	 ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100)*nentSemiCent/nentSemiCentall,
	 ffSemiCent->GetX(0.99*ffSemiCent->GetParameter(2)),
	 30.*50./(ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100)*nentSemiCent/nentSemiCentall));
  printf("Ratio of raw rates: %.2f +- %.2f\n",
	 rawRatio,
	 rawRatioErr);
  printf("Ratio of rates after PS: %.2f +- %.2f\n",
	 psRatio,
	 psRatioErr);
  printf("Ratio of rates after PS and SC: %.2f +- %.2f\n",
	 psScRatio,
	 psScRatioErr);
  printf("Integral of central trigger: %.2f %% +- %.2f %%\n",
	 centrInt,
	 centrIntErr);
  printf("Integral of central trigger (fit): %.2f %% +- %.2f %%\n",
	 centrIntFit,
	 centrIntFitErr);

  printf("\n\n");
  FILE *fout=fopen(Form("./check.%s.txt",filename),"w");
  if (!fout) {
    printf("Failed to open local result file\n");
    return;
  }
  printf("%.2f %.2f   %.2f %.2f   %.2f %.2f   %.2f %.2f   %.2f %.2f\n\n",
	 rawRatio,rawRatioErr,
	 psRatio,psRatioErr,
	 psScRatio,psScRatioErr,
	 centrInt,centrIntErr,
	 centrIntFit,centrIntFitErr);
  fprintf(fout,"%.2f %.2f   %.2f %.2f   %.2f %.2f   %.2f %.2f   %.2f %.2f\n\n",
	  rawRatio,rawRatioErr,
	  psRatio,psRatioErr,
	  psScRatio,psScRatioErr,
	  centrInt,centrIntErr,
	  centrIntFit,centrIntFitErr);
  fclose(fout);

  new TCanvas;
  h3->Rebin(4);
  TFitResultPtr fitRes = h3->Fit("pol0","S","",0,89);
  printf("Centrality flatness, Chi2=%.3f/%d\n",fitRes->Chi2(),fitRes->Ndf());

  // TH1F *hand = (TH1F*)list->FindObject("fV0Percent");
  // TH1F *handall = (TH1F*)list->FindObject("fV0PercentAll");
  // if (handall) {
  //   TH1F *hmb63 = (TH1F*)list->FindObject("fV0Percent63");
  //   TH1F *hmb63all = (TH1F*)list->FindObject("fV0Percent63All");
  //   printf("V0AND: purity=%f %%\n",hand->Integral(0,hand->FindBin(90.))/(Float_t)handall->GetEntries()*100.);
  //   printf(" MB63: purity=%f %% eff=%f %%\n",hmb63->Integral(0,hmb63->FindBin(90.))/hmb63all->GetEntries()*100.,hmb63->Integral(0,hmb63->FindBin(90.))/hand->Integral(0,hand->FindBin(90.))*100.);
  //   hmb63->Sumw2();
  //   hmb63->Divide(hmb63,hand,1,1,"B");
  //   new TCanvas;
  //   hmb63->SetStats(0);
  //   hmb63->SetTitle("MultA+MultC>=63 efficiency as a function of centrality percentile");
  //   hmb63->Draw();
  // }
}
Commit	Line	Data
c014f07d	1	void checkThr(const char filename, const char what="")
	2	{
	3	TFile *f = TFile::Open(filename);
	4	TList list = (TList)f->Get("coutput");
	5
	6	TH1F h3 = (TH1F)list->FindObject("fV0Percent");
8d939043	7	Double_t nTotalEvts = (Double_t)h3->Integral(h3->FindBin(0.0),h3->FindBin(91.0))/0.9;
c014f07d	8
	9	TH1F *hCent;
	10	TH1F *hCentAll;
	11	TH1F *hSemiCent;
	12	TH1F *hSemiCentAll;
	13	Double_t nentCent;
	14	Double_t nentCentall;
	15	Double_t nentSemiCent;
	16	Double_t nentSemiCentall;
	17	TF1 *ffCent;
	18	TF1 *ffSemiCent;
	19
	20	hCent = (TH1F*)list->FindObject(Form("fV0Cent%s",what));
	21	hCentAll = (TH1F*)list->FindObject(Form("fV0Cent%sAll",what));
	22	hSemiCent = (TH1F*)list->FindObject(Form("fV0SemiCent%s",what));
	23	hSemiCentAll = (TH1F*)list->FindObject(Form("fV0SemiCent%sAll",what));
	24
	25	nentCent = hCent->GetEntries();
35d3b2cb	26	Double_t nCentEvts = (Double_t)hCent->Integral(hCent->FindBin(0.0),hCent->FindBin(91.0));
c014f07d	27	nentCentall = hCentAll->GetEntries();
c014f07d	28	nentSemiCent = hSemiCent->GetEntries();
35d3b2cb	29	Double_t nSemiCentEvts = (Double_t)hSemiCent->Integral(hSemiCent->FindBin(0.0),hSemiCent->FindBin(91.0));
c014f07d	30	nentSemiCentall = hSemiCentAll->GetEntries();
	31
	32	hCent->Sumw2();
	33	hCent->Divide(hCent,h3,1,1,"B");
	34	hSemiCent->Sumw2();
	35	hSemiCent->Divide(hSemiCent,h3,1,1,"B");
8d939043	36	// for(Int_t ibin = 1; ibin <= hSemiCent->GetXaxis()->GetNbins(); ++ibin) {
	37	// printf("%f ",hSemiCent->GetBinError(ibin));
	38	// if (hSemiCent->GetBinError(ibin) < 1e-6) hSemiCent->SetBinError(ibin,0.01);
	39	// printf("%f ",hSemiCent->GetBinError(ibin));
	40	// }
c014f07d	41
8d939043	42	ffCent = new TF1("ffCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
c014f07d	43	ffCent->SetLineColor(kBlue);
8d939043	44	Double_t par0 = hCent->GetBinCenter(hCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.6 : 0.95));
c014f07d	45	ffCent->SetParameter(0,par0);
8d939043	46	ffCent->SetParameter(1,0.5);
	47	if (strcmp(what,"Tr") == 0)
	48	ffCent->SetParameter(2,0.6);
	49	else
	50	ffCent->SetParLimits(2,0.9999,1.0001);
c014f07d	51	hCent->Fit(ffCent,"R+");
	52	hCent->Fit(ffCent,"R+");
	53	hCent->Fit(ffCent,"R+");
	54
8d939043	55	ffSemiCent = new TF1("ffSemiCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
c014f07d	56	ffSemiCent->SetLineColor(kRed);
8d939043	57	par0 = hSemiCent->GetBinCenter(hSemiCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.01 : 0.95));
c014f07d	58	ffSemiCent->SetParameter(0,par0);
35d3b2cb	59	ffSemiCent->SetParameter(1,0.3);
8d939043	60	printf("%d\n",par0);
	61	if (strcmp(what,"Tr") == 0)
	62	ffSemiCent->SetParameter(2,0.1);
	63	else
	64	ffSemiCent->SetParLimits(2,0.9999,1.0001);
c014f07d	65	hSemiCent->Fit(ffSemiCent,"R+");
	66	hSemiCent->Fit(ffSemiCent,"R+");
	67	hSemiCent->Fit(ffSemiCent,"R+");
	68
8d939043	69	// new TCanvas;
c014f07d	70	hSemiCent->SetLineColor(kBlue);
	71	hSemiCent->SetMarkerColor(kBlue);
	72	hCent->Draw("e");
	73	hSemiCent->SetLineColor(kRed);
	74	hSemiCent->SetMarkerColor(kRed);
	75	hSemiCent->Draw("e same");
	76
8d939043	77	Float_t rawRatio = (strcmp(what,"Tr") == 0) ? (nentSemiCentall/ffSemiCent->GetParameter(2))/(nentCentall/ffCent->GetParameter(2)) : nentSemiCentall/nentCentall;
	78	Float_t rawRatioErr = rawRatio*TMath::Sqrt(1./nentCentall-1./nentSemiCentall);
	79	Float_t psRatio = (strcmp(what,"Tr") == 0) ? (nentSemiCent/ffSemiCent->GetParameter(2))/(nentCent/ffCent->GetParameter(2)) : nentSemiCent/nentCent;
	80	Float_t psRatioErr = psRatio*TMath::Sqrt(1./nentCent-1./nentSemiCent);
	81	Float_t psScRatio = (strcmp(what,"Tr") == 0) ? (nSemiCentEvts/ffSemiCent->GetParameter(2))/(nCentEvts/ffCent->GetParameter(2)) : nSemiCentEvts/nCentEvts;
	82	Float_t psScRatioErr = psScRatio*TMath::Sqrt(1./nCentEvts-1./nSemiCentEvts);
	83
	84	Float_t centrInt = 100.*nCentEvts/ffCent->GetParameter(2)/nTotalEvts;
	85	Float_t centrIntErr = 100.TMath::Sqrt(nCentEvts/ffCent->GetParameter(2)/nTotalEvts(1.-nCentEvts/ffCent->GetParameter(2)/nTotalEvts)/nTotalEvts);
	86	Float_t centrIntFit = ffCent->GetX(0.5*ffCent->GetParameter(2));
	87	Float_t centrIntFitErr = ffCent->GetParError(0);
	88
c014f07d	89	printf("\n\n\n");
	90
	91	printf("Purity: Centr %.3f %.3f Cuts ( 0.0-%.1f)%% rate=%.1f Hz\n",
8d939043	92	nCentEvts/nCentEvts,
	93	nentCent/nentCentall,
	94	ffCent->GetX(0.99*ffCent->GetParameter(2)),
	95	30.*10./(nentCent/nentCentall));
c014f07d	96	printf("Purity: Semi-Centr %.3f %.3f Cuts ( 0.0-%.1f)%% rate=%.1f Hz\n",
	97	ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100),
	98	ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100)*nentSemiCent/nentSemiCentall,
8d939043	99	ffSemiCent->GetX(0.99*ffSemiCent->GetParameter(2)),
35d3b2cb	100	30.50./(ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100)nentSemiCent/nentSemiCentall));
8d939043	101	printf("Ratio of raw rates: %.2f +- %.2f\n",
	102	rawRatio,
	103	rawRatioErr);
	104	printf("Ratio of rates after PS: %.2f +- %.2f\n",
	105	psRatio,
	106	psRatioErr);
	107	printf("Ratio of rates after PS and SC: %.2f +- %.2f\n",
	108	psScRatio,
	109	psScRatioErr);
	110	printf("Integral of central trigger: %.2f %% +- %.2f %%\n",
	111	centrInt,
	112	centrIntErr);
	113	printf("Integral of central trigger (fit): %.2f %% +- %.2f %%\n",
	114	centrIntFit,
	115	centrIntFitErr);
c014f07d	116
8d939043	117	printf("\n\n");
	118	FILE *fout=fopen(Form("./check.%s.txt",filename),"w");
	119	if (!fout) {
	120	printf("Failed to open local result file\n");
	121	return;
c014f07d	122	}
8d939043	123	printf("%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n\n",
	124	rawRatio,rawRatioErr,
	125	psRatio,psRatioErr,
	126	psScRatio,psScRatioErr,
	127	centrInt,centrIntErr,
	128	centrIntFit,centrIntFitErr);
	129	fprintf(fout,"%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n\n",
	130	rawRatio,rawRatioErr,
	131	psRatio,psRatioErr,
	132	psScRatio,psScRatioErr,
	133	centrInt,centrIntErr,
	134	centrIntFit,centrIntFitErr);
	135	fclose(fout);
	136
	137	new TCanvas;
	138	h3->Rebin(4);
	139	TFitResultPtr fitRes = h3->Fit("pol0","S","",0,89);
	140	printf("Centrality flatness, Chi2=%.3f/%d\n",fitRes->Chi2(),fitRes->Ndf());
	141
	142	// TH1F hand = (TH1F)list->FindObject("fV0Percent");
	143	// TH1F handall = (TH1F)list->FindObject("fV0PercentAll");
	144	// if (handall) {
	145	// TH1F hmb63 = (TH1F)list->FindObject("fV0Percent63");
	146	// TH1F hmb63all = (TH1F)list->FindObject("fV0Percent63All");
	147	// printf("V0AND: purity=%f %%\n",hand->Integral(0,hand->FindBin(90.))/(Float_t)handall->GetEntries()*100.);
	148	// printf(" MB63: purity=%f %% eff=%f %%\n",hmb63->Integral(0,hmb63->FindBin(90.))/hmb63all->GetEntries()100.,hmb63->Integral(0,hmb63->FindBin(90.))/hand->Integral(0,hand->FindBin(90.))100.);
	149	// hmb63->Sumw2();
	150	// hmb63->Divide(hmb63,hand,1,1,"B");
	151	// new TCanvas;
	152	// hmb63->SetStats(0);
	153	// hmb63->SetTitle("MultA+MultC>=63 efficiency as a function of centrality percentile");
	154	// hmb63->Draw();
	155	// }
c014f07d	156	}