Adding analysis task for the centrality trigger calibration.
[u/mrichter/AliRoot.git] / PWG1 / VZERO / extractThr.C
CommitLineData
c014f07d 1void 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
41 Int_t nb = task->GetNThr();
42 new TCanvas;
43 h3->Sumw2();
44 TH1F **hh = new TH1F*[nb];
45 TH1F **hhall = new TH1F*[nb];
46 TF1 **ff = new TF1*[nb];
47 Double_t *nent = new Double_t[nb];
48 Double_t *nentall = new Double_t[nb];
49 for(Int_t i = 0; i < nb; ++i) {
50 hh[i] = (TH1F*)list->FindObject(Form("fV0PercentBins_%d",i));
51 hhall[i] = (TH1F*)list->FindObject(Form("fV0PercentBinsAll_%d",i));
52
53 nent[i] = hh[i]->GetEntries();
54 nentall[i] = hhall[i]->GetEntries();
55
56 hh[i]->Sumw2();
57 hh[i]->Divide(hh[i],h3,1,1,"B");
58 // TF1 *ff = new TF1(Form("ff_%d",i),"x > [0] ? TMath::Exp(-(x-[0])*(x-[0])/(2.*[1]*[1])) : 1.0",0,100);
59 ff[i] = new TF1(Form("ff_%d",i),"1.-1./(1.+TMath::Exp(-(x-[0])/[1]))",0,100);
60 ff[i]->SetLineColor(i%9+1);
61 printf("%d\n",i);
62 if (nent[i] > 1000) {
63 Double_t par0 = hh[i]->GetBinCenter(hh[i]->FindLastBinAbove(0.9));
64 printf("%f\n",par0);
65 ff[i]->SetParameter(0,par0);
66 ff[i]->SetParameter(1,0.5);
67 hh[i]->Fit(ff[i],"R+");
68 hh[i]->Fit(ff[i],"R+");
69 hh[i]->Fit(ff[i],"R+");
70 }
71 }
72
73 new TCanvas;
74 Double_t *eff99 = new Double_t[nb];
75 Double_t *eff05 = new Double_t[nb];
76 Double_t *purity = new Double_t[nb];
77 Double_t *purity2 = new Double_t[nb];
78 Bool_t first = kTRUE;
79 for(Int_t i = 0; i < nb; ++i) {
80 hh[i]->SetLineColor(i%9+1);
81 hh[i]->SetMarkerColor(i%9+1);
82 hh[i]->SetLineWidth(2.0);
83 hh[i]->SetStats(0);
84 if (first) {
85 first = kFALSE;
86 hh[i]->Draw("e");
87 hh[i]->GetXaxis()->SetTitle("Centrality percentile");
88 hh[i]->GetYaxis()->SetTitle("Efficiency");
89 }
90 else
91 hh[i]->Draw("e same");
92 eff99[i] = ff[i]->GetX(0.99);
93 eff05[i] = ff[i]->GetX(0.05)-eff99[i];
94 purity[i] = ff[i]->Integral(0,eff99[i])/ff[i]->Integral(0,100);
95 purity2[i] = purity[i]*nent[i]/nentall[i];
96 // printf("%d %d %d\n",i,nent[i],nentall[i]);
97 }
98
99 new TCanvas;
100 TGraph *gr99 = new TGraph(nb,eff99,eff05);
101 gr99->SetMarkerStyle(21);
102 gr99->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
103 gr99->GetHistogram()->GetYaxis()->SetTitle("#Delta(Centrality percentile @ Eff=5% - @ Eff=99%)");
104 gr99->GetHistogram()->GetYaxis()->SetRangeUser(-7,10);
105 gr99->SetTitle("");
106 gr99->Draw("AP");
107
108 TLegend *leg0 = new TLegend(0.75,0.55,0.89,0.89);
109 leg0->AddEntry(gr99,"A and C sides","P");
110 leg0->Draw();
111
112 new TCanvas;
113 TGraph *pur = new TGraph(nb,eff99,purity);
114 pur->SetMarkerStyle(21);
115 pur->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
116 pur->GetHistogram()->GetYaxis()->SetTitle("Purity");
117 pur->SetTitle("");
118 TF1 *fitPur = new TF1("fitPur","[0]-[1]*TMath::Exp(-x/[2])",6,50);
119 fitPur->SetParameter(0,0.9);
120 fitPur->SetParameter(1,0.4);
121 fitPur->SetParameter(2,9.0);
122 pur->Fit(fitPur,"R+");
123 pur->Fit(fitPur,"R+");
124 pur->Fit(fitPur,"R+");
125 pur->Draw("AP");
126
127 TGraph *pur2 = new TGraph(nb,eff99,purity2);
128 pur2->SetMarkerStyle(22);
129 pur2->GetHistogram()->GetXaxis()->SetTitle("Centrality percentile @ Eff=99%");
130 pur2->GetHistogram()->GetYaxis()->SetTitle("Purity (no phys & centr selection)");
131 pur2->SetTitle("");
132 TF1 *fitPur2 = new TF1("fitPur2","[0]-[1]*TMath::Exp(-x/[2])",6,50);
133 fitPur2->SetParameter(0,0.9);
134 fitPur2->SetParameter(1,0.4);
135 fitPur2->SetParameter(2,9.0);
136 pur2->Fit(fitPur2,"R+");
137 pur2->Fit(fitPur2,"R+");
138 pur2->Fit(fitPur2,"R+");
139 pur2->Draw("Psame");
140
141 TLegend *leg = new TLegend(0.75,0.55,0.89,0.89);
142 leg->AddEntry(pur,"A and C, Phys & centr selection","P");
143 leg->AddEntry(pur2,"A and C, No event selection","P");
144 leg->Draw();
145
146 new TCanvas;
147 Double_t *thrA = new Double_t[nb];
148 Double_t *thrC = new Double_t[nb];
149 Double_t ratio = task->GetRatio();
150 for(Int_t i = 0; i < nb; ++i) {
151 thrA[i] = task->GetThrA(i);
152 thrC[i] = task->GetThrC(i);
153 }
154 TGraph *thr = new TGraph(nb,eff99,thrA);
155 thr->SetMarkerStyle(21);
156 TF1 *fThr = new TF1("fThr","[0]*x*x*x+[1]*x*x+[2]*x+[3]",1.5,50);
157 thr->Fit(fThr,"R");
158 thr->Draw("AP");
159
160 TH2F *hall = (TH2F*)list->FindObject("fV0Charge2dAll");
161 h2->Draw("colz");
162 TLine *line0 = new TLine(fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio(),25000,fThr->Eval(centrValue)*task->GetRatio());
163 line0->Draw("same");
164 TLine *line1 = new TLine(fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio(),fThr->Eval(centrValue),25000);
165 line1->Draw("same");
166
167 TLine *line2 = new TLine(fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),25000,fThr->Eval(semiCentrValue)*task->GetRatio());
168 line2->Draw("same");
169 TLine *line3 = new TLine(fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),fThr->Eval(semiCentrValue),25000);
170 line3->Draw("same");
171
172 printf("Slope param: %f %f\n",f1->GetParameter(0),f2->GetParameter(0));
173 printf("Mean slope param: %f\n",0.5*(f1->GetParameter(0)+f2->GetParameter(0)));
174 printf("Delta slope param: %f %%\n",100.*(f1->GetParameter(0)-f2->GetParameter(0))/(0.5*(f1->GetParameter(0)+f2->GetParameter(0))));
175
176 printf("A&C @ %.1f %% ThrA = %f ThrC = %f\n",centrValue,fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio());
177 printf("A&C @ %.1f %% ThrA = %f ThrC = %f\n",semiCentrValue,fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio());
178
179 printf("\n\n");
180 FILE *fout=fopen("./trigger.txt","w");
181 if (!fout) {
182 printf("Failed to open local result file\n");
183 return;
184 }
185 printf("%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n\n",
186 task->GetMinThr(),
187 task->GetMaxThr(),
188 0.5*(f1->GetParameter(0)+f2->GetParameter(0)),
189 task->GetNThr(),
190 fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),
191 fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio());
192 fprintf(fout,"%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n",
193 task->GetMinThr(),
194 task->GetMaxThr(),
195 0.5*(f1->GetParameter(0)+f2->GetParameter(0)),
196 task->GetNThr(),
197 fThr->Eval(semiCentrValue),fThr->Eval(semiCentrValue)*task->GetRatio(),
198 fThr->Eval(centrValue),fThr->Eval(centrValue)*task->GetRatio());
199 fclose(fout);
200}