PWGPP/VZERO/checkThr.C

   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");
   7   Double_t nTotalEvts = (Double_t)h3->Integral(h3->FindBin(0.0),h3->FindBin(91.0))/0.9;
   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();
  26   Double_t nCentEvts = (Double_t)hCent->Integral(hCent->FindBin(0.0),hCent->FindBin(91.0));
  27   nentCentall = hCentAll->GetEntries();
  28   nentSemiCent = hSemiCent->GetEntries();
  29   Double_t nSemiCentEvts = (Double_t)hSemiCent->Integral(hSemiCent->FindBin(0.0),hSemiCent->FindBin(91.0));
  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");
  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   // }
  41
  42   ffCent = new TF1("ffCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
  43   ffCent->SetLineColor(kBlue);
  44   Double_t par0 = hCent->GetBinCenter(hCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.6 : 0.95));
  45   ffCent->SetParameter(0,par0);
  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);
  51   hCent->Fit(ffCent,"R+");
  52   hCent->Fit(ffCent,"R+");
  53   hCent->Fit(ffCent,"R+");
  54
  55   ffSemiCent = new TF1("ffSemiCent","(1.-1./(1.+TMath::Exp(-(x-[0])/[1])))*[2]",0,100);
  56   ffSemiCent->SetLineColor(kRed);
  57   par0 = hSemiCent->GetBinCenter(hSemiCent->FindLastBinAbove((strcmp(what,"Tr") == 0) ? 0.01 : 0.95));
  58   ffSemiCent->SetParameter(0,par0);
  59   ffSemiCent->SetParameter(1,0.3);
  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);
  65   hSemiCent->Fit(ffSemiCent,"R+");
  66   hSemiCent->Fit(ffSemiCent,"R+");
  67   hSemiCent->Fit(ffSemiCent,"R+");
  68
  69   //  new TCanvas;
  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
  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
  89   printf("\n\n\n");
  90
  91   printf("Purity: Centr      %.3f %.3f   Cuts ( 0.0-%.1f)%%   rate=%.1f Hz\n",
  92          nCentEvts/nCentEvts,
  93          nentCent/nentCentall,
  94          ffCent->GetX(0.99*ffCent->GetParameter(2)),
  95          30.*10./(nentCent/nentCentall));
  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,
  99          ffSemiCent->GetX(0.99*ffSemiCent->GetParameter(2)),
 100          30.*50./(ffSemiCent->Integral(0,50)/ffSemiCent->Integral(0,100)*nentSemiCent/nentSemiCentall));
 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);
 116
 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;
 122   }
 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   // }
 156 }