PWGPP/VZERO/extractThr.C

   1 void extractThr(const char *filename, const char *fparname, Float_t centrValue = 10.0, Float_t semiCentrValue = 50.0)
   2 {
   3   gROOT->LoadMacro(Form("%s/AliAnaVZEROTrigger.cxx+",
   4                         gSystem->pwd()));
   5
   6   AliAnaVZEROTrigger *task = new AliAnaVZEROTrigger;
   7   task->Setup(fparname);
   8
   9   TFile *f = TFile::Open(filename);
  10
  11   TList *list = (TList*)f->Get("coutput");
  12
  13   TH2F *h1 = (TH2F*)list->FindObject("fV0Charge2d");
  14   TH2F *h2 = (TH2F*)list->FindObject("fV0Charge2dPercent");
  15   h2->Divide(h1);
  16   TProfile *h1_pfx = h1->ProfileX("h1_pfx");
  17   TProfile *h1_pfy = h1->ProfileY("h1_pfy");
  18   new TCanvas;
  19   gStyle->SetPalette(1,0);
  20   h2->Draw("colz");
  21   new TCanvas;
  22   TF1 *f1 = new TF1("f1","[0]*x",0,task->GetCentCuts(0));
  23   f1->SetParameter(0,task->GetRatio());
  24   h1_pfx->Fit(f1,"WR+");
  25   h1_pfx->Fit(f1,"WR+");
  26   h1_pfx->Fit(f1,"WR+");
  27   h1_pfx->Draw();
  28
  29   new TCanvas;
  30   TF1 *f2 = new TF1("f2","x/[0]",0,task->GetCentCuts(0)*f1->GetParameter(0));
  31   f2->SetParameter(0,1./task->GetRatio());
  32   h1_pfy->Fit(f2,"WR+");
  33   h1_pfy->Fit(f2,"WR+");
  34   h1_pfy->Fit(f2,"WR+");
  35   h1_pfy->Draw();
  36
  37   TH1F *h3 = (TH1F*)list->FindObject("fV0Percent");
  38   new TCanvas;
  39   h3->Draw();
  40   Double_t nTotalEvts = (Double_t)h3->Integral(h3->FindBin(0.0),h3->FindBin(91.0))/0.9;
  41
  42   Int_t nb = task->GetNThr();
  43   new TCanvas;
  44   h3->Sumw2();
  45   TH1F **hh = new TH1F*[nb];
  46   TH1F **hhall = new TH1F*[nb];
  47   TF1 **ff = new TF1*[nb];
  48   Double_t *nent = new Double_t[nb];
  49   Double_t *nBinEvts = new Double_t[nb];
  50   Double_t *nentall = new Double_t[nb];
  51   for(Int_t i = 0; i < nb; ++i) {
  52     hh[i] = (TH1F*)list->FindObject(Form("fV0PercentBins_%d",i));
  53     hhall[i] = (TH1F*)list->FindObject(Form("fV0PercentBinsAll_%d",i));
  54
  55     nent[i] = hh[i]->GetEntries();
  56     nBinEvts[i] = (Double_t)hh[i]->Integral(hh[i]->FindBin(0.0),hh[i]->FindBin(91.0));
  57     nentall[i] = hhall[i]->GetEntries();
  58
  59     hh[i]->Sumw2();
  60     hh[i]->Divide(hh[i],h3,1,1,"B");
  61     //    TF1 *ff = new TF1(Form("ff_%d",i),"x > [0] ? TMath::Exp(-(x-[0])*(x-[0])/(2.*[1]*[1])) : 1.0",0,100);
  62     ff[i] = new TF1(Form("ff_%d",i),"1.-1./(1.+TMath::Exp(-(x-[0])/[1]))",0,100);
  63     ff[i]->SetLineColor(i%9+1);
  64     printf("%d\n",i);
  65     if (nent[i] > 100) {
  66       Double_t par0 = hh[i]->GetBinCenter(hh[i]->FindLastBinAbove(0.9));
  67       printf("%f\n",par0);
  68       ff[i]->SetParameter(0,par0);
  69       ff[i]->SetParameter(1,0.1);
  70       //      ff[i]->SetParLimits(1,0.01,2.);
  71       hh[i]->Fit(ff[i],"R+");
  72       hh[i]->Fit(ff[i],"R+");
  73       hh[i]->Fit(ff[i],"R+");
  74     }
  75   }
  76
  77   new TCanvas;
  78   Double_t *eff99 = new Double_t[nb];
  79   Double_t *eff05 = new Double_t[nb];
  80   Double_t *purity = new Double_t[nb];
  81   Double_t *purity2 = new Double_t[nb];
  82   Bool_t first = kTRUE;
  83   for(Int_t i = 0; i < nb; ++i) {
  84     hh[i]->SetLineColor(i%9+1);
  85     hh[i]->SetMarkerColor(i%9+1);
  86     hh[i]->SetLineWidth(2.0);
  87     hh[i]->SetStats(0);
  88     if (first) {
  89       first = kFALSE;
  90       hh[i]->Draw("e");
  91       hh[i]->GetXaxis()->SetTitle("Centrality percentile");
  92       hh[i]->GetYaxis()->SetTitle("Efficiency");
  93     }
  94     else
  95       hh[i]->Draw("e same");
  96     eff99[i] = ff[i]->GetX(0.99);
  97     eff05[i] = ff[i]->GetX(0.05)-eff99[i];
  98     purity[i] = (ff[i]->Integral(0,100)>0) ? ff[i]->Integral(0,eff99[i])/ff[i]->Integral(0,100) : 0;
  99     purity2[i] = (nentall[i]>0) ? purity[i]*nent[i]/nentall[i] : 0;
 100     //    printf("%d  %d %d\n",i,nent[i],nentall[i]);
 101   }
 102
 103   new TCanvas;
 104   TGraph *gr99 = new TGraph(nb,eff99,eff05);
 105   gr99->SetMarkerStyle(21);
 106   gr99->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
 107   gr99->GetHistogram()->GetYaxis()->SetTitle("#Delta(Centrality percentile @ Eff=5% - @ Eff=99%)");
 108   gr99->GetHistogram()->GetYaxis()->SetRangeUser(-7,10);
 109   gr99->SetTitle("");
 110   gr99->Draw("AP");
 111
 112   TLegend *leg0 = new TLegend(0.75,0.55,0.89,0.89);
 113   leg0->AddEntry(gr99,"A and C sides","P");
 114   leg0->Draw();
 115
 116   new TCanvas;
 117   TGraph *pur = new TGraph(nb,eff99,purity);
 118   pur->SetMarkerStyle(21);
 119   pur->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
 120   pur->GetHistogram()->GetYaxis()->SetTitle("Purity");
 121   pur->SetTitle("");
 122   TF1 *fitPur = new TF1("fitPur","[0]-[1]*TMath::Exp(-x/[2])",6,50);
 123   fitPur->SetParameter(0,0.9);
 124   fitPur->SetParameter(1,0.4);
 125   fitPur->SetParameter(2,9.0);
 126   pur->Fit(fitPur,"R+");
 127   pur->Fit(fitPur,"R+");
 128   pur->Fit(fitPur,"R+");
 129   pur->Draw("AP");
 130
 131   TGraph *pur2 = new TGraph(nb,eff99,purity2);
 132   pur2->SetMarkerStyle(22);
 133   pur2->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
 134   pur2->GetHistogram()->GetYaxis()->SetTitle("Purity (no phys & centr selection)");
 135   pur2->SetTitle("");
 136   TF1 *fitPur2 = new TF1("fitPur2","[0]-[1]*TMath::Exp(-x/[2])",6,50);
 137   fitPur2->SetParameter(0,0.9);
 138   fitPur2->SetParameter(1,0.4);
 139   fitPur2->SetParameter(2,9.0);
 140   pur2->Fit(fitPur2,"R+");
 141   pur2->Fit(fitPur2,"R+");
 142   pur2->Fit(fitPur2,"R+");
 143   pur2->Draw("Psame");
 144
 145   TLegend *leg = new TLegend(0.75,0.55,0.89,0.89);
 146   leg->AddEntry(pur,"A and C, Phys & centr selection","P");
 147   leg->AddEntry(pur2,"A and C, No event selection","P");
 148   leg->Draw();
 149
 150   new TCanvas;
 151   Double_t *thrA = new Double_t[nb];
 152   Double_t *thrC = new Double_t[nb];
 153   Double_t ratio = task->GetRatio();
 154   for(Int_t i = 0; i < nb; ++i) {
 155     thrA[i] = task->GetThrA(i);
 156     thrC[i] = task->GetThrC(i);
 157   }
 158   TGraph *thr = new TGraph(nb,eff99,thrA);
 159   thr->SetMarkerStyle(21);
 160   TF1 *fThr = new TF1("fThr","[0]*x*x*x*x+[1]*x*x*x+[2]*x*x+[3]*x+[4]",3.0,60);
 161   thr->Fit(fThr,"R");
 162   thr->Draw("AP");
 163
 164   new TCanvas;
 165   Double_t *intThr = new Double_t[nb];
 166   for(Int_t i = 0; i < nb; ++i) {
 167     intThr[i] = 100.*nBinEvts[i]/nTotalEvts;
 168   }
 169   TGraph *grInt = new TGraph(nb,intThr,thrA);
 170   grInt->SetMarkerStyle(21);
 171   TF1 *fIntThr = new TF1("fIntThr","[0]*x*x*x*x+[1]*x*x*x+[2]*x*x+[3]*x+[4]",3.0,60);
 172   grInt->Fit(fIntThr,"R");
 173   grInt->Draw("AP");
 174
 175   new TCanvas;
 176   TH2F *hall = (TH2F*)list->FindObject("fV0Charge2dAll");
 177   h2->Draw("colz");
 178   TLine *line0 = new TLine(fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio(),25000,fThr->Eval(centrValue)*task->GetRatio());
 179   line0->Draw("same");
 180   TLine *line1 = new TLine(fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio(),fThr->Eval(centrValue),25000);
 181   line1->Draw("same");
 182
 183   TLine *line01 = new TLine(fIntThr->Eval(centrValue),fIntThr->Eval(centrValue)*task->GetRatio(),25000,fIntThr->Eval(centrValue)*task->GetRatio());
 184   line01->Draw("same");
 185   TLine *line11 = new TLine(fIntThr->Eval(centrValue),fIntThr->Eval(centrValue)*task->GetRatio(),fIntThr->Eval(centrValue),25000);
 186   line11->Draw("same");
 187
 188   TLine *line2 = new TLine(fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),25000,fThr->Eval(semiCentrValue)*task->GetRatio());
 189   line2->Draw("same");
 190   TLine *line3 = new TLine(fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),fThr->Eval(semiCentrValue),25000);
 191   line3->Draw("same");
 192
 193   printf("Slope param: %f %f\n",f1->GetParameter(0),f2->GetParameter(0));
 194   printf("Mean slope param: %f\n",0.5*(f1->GetParameter(0)+f2->GetParameter(0)));
 195   printf("Delta slope param: %f %%\n",100.*(f1->GetParameter(0)-f2->GetParameter(0))/(0.5*(f1->GetParameter(0)+f2->GetParameter(0))));
 196
 197   printf("A&C @ %.1f %%   ThrA = %f   ThrC = %f\n",centrValue,fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio());
 198   printf("A&C @ %.1f %%   ThrA = %f   ThrC = %f\n",semiCentrValue,fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio());
 199
 200   printf("A&C (Integral) @ %.1f %%   ThrA = %f   ThrC = %f\n",centrValue,fIntThr->Eval(centrValue),fIntThr->Eval(centrValue)*task->GetRatio());
 201   printf("A&C (Intergal) @ %.1f %%   ThrA = %f   ThrC = %f\n",semiCentrValue,fIntThr->Eval(semiCentrValue),fIntThr->Eval(semiCentrValue)*task->GetRatio());
 202
 203   printf("\n\n");
 204   FILE *fout=fopen("./trigger.txt","w");
 205   if (!fout) {
 206     printf("Failed to open local result file\n");
 207     return;
 208   }
 209   printf("%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n\n",
 210          task->GetMinThr(),
 211          task->GetMaxThr(),
 212          0.5*(f1->GetParameter(0)+f2->GetParameter(0)),
 213          task->GetNThr(),
 214          fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),
 215          fIntThr->Eval(centrValue),fIntThr->Eval(centrValue)*task->GetRatio());
 216   fprintf(fout,"%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n",
 217           task->GetMinThr(),
 218           task->GetMaxThr(),
 219           0.5*(f1->GetParameter(0)+f2->GetParameter(0)),
 220           task->GetNThr(),
 221           fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),
 222           fIntThr->Eval(centrValue),fIntThr->Eval(centrValue)*task->GetRatio());
 223   fclose(fout);
 224 }