]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_PbPb/macros/QA/ExtractQA.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Improved labeling of QA histograms
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / macros / QA / ExtractQA.C
CommitLineData
a600ad2c 1#include <TCanvas.h>
2
3#if !defined(__CINT__) || defined(__MAKECINT__)
4#include <TGrid.h>
5#include <TFile.h>
6#include <TF1.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TList.h>
10#include <TMath.h>
11#include <TGraphErrors.h>
12#include <TString.h>
13#include "Riostream.h"
14#include "stdio.h"
15using namespace std;
16#endif
17
18//-----------------------------------------------------------------------------
19// Function declaration
20void QAFillEventSelection();
21void QAFillOccupancy();
22void QAFillClusters();
23void QAFillTracks();
24void QAFillNPi0();
25
26void QAWriteEventSelection();
27void QAWriteOccupancy();
28void QAWriteClusters();
29void QAWriteTracks();
30void QAWriteNPi0();
31
32TH1D *Pi0InvMass(TH2 *hMassPt, Float_t ptmin, Float_t ptmax);
b5bff3e2 33TH1F *AddWriteTH1F(const char *name, const char *title, Double_t* values, Double_t* errors);
a600ad2c 34void Fit1Pi0(TH1D *hMass,
35 const Int_t polN,
36 const Double_t mFitMin,
37 const Double_t mFitMax,
38 Double_t &nPi0, Double_t &ePi0);
39Double_t pi0massP1(Double_t *x, Double_t *par);
40Double_t pi0massP2(Double_t *x, Double_t *par);
41Double_t pi0mass (Double_t *x, Double_t *par);
42Double_t bgP2 (Double_t *x, Double_t *par);
43
44//-----------------------------------------------------------------------------
45// Global variabes
b5bff3e2 46const Int_t kNCents = 1;
47const Int_t kNPID = 8+4;
48const char* kPIDNames[kNPID] = {"All", "Allwou", "Disp", "Disp2", "Dispwou", "CPV", "CPV2", "Both",
49 "Allcore", "Dispcore", "CPVcore", "Bothcore"};
a600ad2c 50
51Int_t runIndex;
52Int_t runNum[500];
53TList *listHist;
54Double_t run[500], erun[500];
b5bff3e2 55// Event selection:
ff329095 56Double_t nTotal4[500];
a600ad2c 57Double_t rVtx[500], rVtxZ10[500], rVtxZ10Cent[500];
58Double_t eVtx[500], eVtxZ10[500], eVtxZ10Cent[500];
b5bff3e2 59// Occupancy:
a600ad2c 60Double_t nCellsM1[500], nCellsM2[500], nCellsM3[500];
b5bff3e2 61// Clusters:
a600ad2c 62Double_t nCluEvent[500], eCluEvent[500];
63Double_t cluEnergy[500], ecluEnergy[500];
b5bff3e2 64Double_t nPhotPID[kNCents][kNPID][500], enPhotPID[kNCents][kNPID][500];
65Double_t mEnPID[kNCents][kNPID][500], emEnPID[kNCents][kNPID][500];
ff329095 66Double_t nPhotPIDHigh[kNCents][kNPID][500], enPhotPIDHigh[kNCents][kNPID][500];
67Double_t mEnPIDHigh[kNCents][kNPID][500], emEnPIDHigh[kNCents][kNPID][500];
68
b5bff3e2 69
70// Tracks:
a600ad2c 71Double_t nTracks0[500], eTracks0[500];
72Double_t nTracks1[500], eTracks1[500];
b5bff3e2 73// Pi0:
74Double_t mPi0[500], emPi0[500];
75Double_t wPi0[500], ewPi0[500];
76Double_t nPi0Event[500], ePi0Event[500];
77
a600ad2c 78
79//-----------------------------------------------------------------------------
80void ExtractQA(const TString filelist="filelist.txt",
3b8b9fc1 81 const TString runlist="runlist.txt",
82 const TString outputFileName = "outputQA.root")
a600ad2c 83{
84 // Loop over per-run root files and run various QA functions
85
86 //TGrid::Connect("alien://");
87
88 ifstream in;
89 in.open(filelist.Data());
90
91 ifstream inRuns;
92 inRuns.open(runlist.Data());
93
94 char rootFileName[256];
95 TFile *rootFile;
96
97 Int_t runNumber=0;
98 runIndex = 0;
99 while (1) {
100 in >> rootFileName;
101 if (!in.good()) break;
102 inRuns >> runNumber;
103 if (!inRuns.good()) break;
104
105 if(169094 == runNumber )
106 continue;
107
108 runNum[runIndex] = runNumber;
109
110 printf("root file is %s, run # = %d\n",rootFileName,runNumber);
111 // char *runNum = strtok(rootFileName+35,".");
112 rootFile = TFile::Open(rootFileName,"read");
113 listHist = (TList*)rootFile->Get("PHOSPi0Flow/PHOSPi0FlowCoutput1");
114
115 run[runIndex] = runIndex+1;
116
117 QAFillEventSelection();
118 QAFillOccupancy();
119 QAFillClusters();
120 QAFillTracks();
3b8b9fc1 121 QAFillNPi0();
a600ad2c 122
123 listHist->Clear();
124 rootFile->Close();
125 runIndex++;
126 }
127
128
3b8b9fc1 129 TFile *fileQA = TFile::Open(outputFileName.Data(), "recreate");
a600ad2c 130 QAWriteEventSelection();
131 QAWriteOccupancy();
132 QAWriteClusters();
133 QAWriteTracks();
3b8b9fc1 134 QAWriteNPi0();
a600ad2c 135 fileQA ->Close();
136
137}
138
139//-----------------------------------------------------------------------------
140void QAFillEventSelection()
141{
142 TH1F *hSelEvents = (TH1F*)listHist->FindObject("hTotSelEvents");
143
144 Double_t nTotal = hSelEvents->GetBinContent(1);
145 Int_t nVtx = hSelEvents->GetBinContent(2);
146 Int_t nVtxZ10 = hSelEvents->GetBinContent(3);
147 Int_t nVtxZ10Cent = hSelEvents->GetBinContent(4);
148
149 if( ! nTotal )
150 return;
151
ff329095 152 nTotal4[runIndex] = hSelEvents->GetBinContent(4);
a600ad2c 153
154 rVtx[runIndex] = nVtx /nTotal;
155 rVtxZ10[runIndex] = nVtxZ10 /nTotal;
156 rVtxZ10Cent[runIndex] = nVtxZ10Cent /nTotal;
157
158
159 //Change in logic (plus clarification through rewriting): t(1+p)/n/n -> p(1-p)/n:
160 eVtx[runIndex] = TMath::Sqrt( rVtx[runIndex] * (1 - rVtx[runIndex]) /nTotal);
161 eVtxZ10[runIndex] = TMath::Sqrt( rVtxZ10[runIndex] * (1 - rVtxZ10[runIndex]) /nTotal);
162 eVtxZ10Cent[runIndex] = TMath::Sqrt(rVtxZ10Cent[runIndex] * (1 - rVtxZ10Cent[runIndex]) /nTotal);
163
164}
165
166//-----------------------------------------------------------------------------
167void QAFillOccupancy()
168{
169 TH2D *hCellNXZM1 = (TH2D*)listHist->FindObject("hCellNXZM1");
170 TH2D *hCellNXZM2 = (TH2D*)listHist->FindObject("hCellNXZM2");
171 TH2D *hCellNXZM3 = (TH2D*)listHist->FindObject("hCellNXZM3");
172
173 // Count cells in each module
174
175 Int_t nCells1=0, nCells2=0, nCells3=0;
176 for (Int_t cellX=1; cellX<=64; cellX++) {
177 for (Int_t cellZ=1; cellZ<=56; cellZ++) {
178 if (hCellNXZM1->GetBinContent(cellX,cellZ) > 0) nCells1++;
179 if (hCellNXZM2->GetBinContent(cellX,cellZ) > 0) nCells2++;
180 if (hCellNXZM3->GetBinContent(cellX,cellZ) > 0) nCells3++;
181 }
182 }
183 nCellsM1[runIndex] = nCells1;
184 nCellsM2[runIndex] = nCells2;
185 nCellsM3[runIndex] = nCells3;
186}
187
188//-----------------------------------------------------------------------------
189void QAFillClusters()
190{
191 TH1F *hTotSelEvents = (TH1F*)listHist->FindObject("hTotSelEvents");
b5bff3e2 192 Double_t nEvents4 = hTotSelEvents->GetBinContent(4);
193 if( ! nEvents4 ) return;
194
a600ad2c 195 TH2F *hCluEvsClu = (TH2F*)listHist->FindObject("hCluEvsClu");
196 TH1D *hClusterE = hCluEvsClu->ProjectionX("hClusterE",4,41);
a600ad2c 197
b5bff3e2 198 cluEnergy[runIndex] = hClusterE->GetMean();
199 ecluEnergy[runIndex] = hClusterE->GetMeanError();
a600ad2c 200
b5bff3e2 201 nCluEvent[runIndex] = hClusterE->Integral() / nEvents4;
202 eCluEvent[runIndex] = TMath::Sqrt( hClusterE->Integral() )/nEvents4;
a600ad2c 203
b5bff3e2 204
205 for(int cent = 0; cent < kNCents; ++cent) {
206 for(int ipid = 0; ipid < kNPID; ++ipid) {
207 TH1* hPhot = listHist->FindObject( Form("hPhot%s_cen%d", kPIDNames[ipid], cent) );
208
ff329095 209 hPhot->SetAxisRange(0., 100.);
b5bff3e2 210 double nPhot = hPhot->Integral() /nEvents4;
211 nPhotPID [cent][ipid][runIndex] = nPhot;
212 enPhotPID[cent][ipid][runIndex] = TMath::Sqrt( nPhot /nEvents4 );
213 mEnPID [cent][ipid][runIndex] = hPhot->GetMean();
214 emEnPID [cent][ipid][runIndex] = hPhot->GetMeanError();
ff329095 215
216 hPhot->SetAxisRange(1., 100.);
217 double nPhotHigh = hPhot->Integral() /nEvents4;
218 nPhotPIDHigh [cent][ipid][runIndex] = nPhotHigh;
219 enPhotPIDHigh[cent][ipid][runIndex] = TMath::Sqrt( nPhotHigh /nEvents4 );
220 mEnPIDHigh [cent][ipid][runIndex] = hPhot->GetMean();
221 emEnPIDHigh [cent][ipid][runIndex] = hPhot->GetMeanError();
222
b5bff3e2 223 }
224 }
a600ad2c 225}
226
227//-----------------------------------------------------------------------------
228void QAFillTracks()
229{
230 TH2F *hCenTrack = (TH2F*)listHist->FindObject("hCenTrack");
231 TH1D *hTrackMult = hCenTrack->ProjectionY("hTrackMult", 1, hCenTrack->GetNbinsY());
232
233 nTracks0[runIndex] = hTrackMult->GetMean();
234 eTracks0[runIndex] = hTrackMult->GetMeanError();
235
236// hTrackMult->DrawCopy();
237// hTrackMult->SetAxisRange(5., hTrackMult->GetXaxis()->GetXmax());
238// new TCanvas;
239// hTrackMult->DrawCopy();
240
241// nTracks1[runIndex] = hTrackMult->GetMean();
242// eTracks1[runIndex] = hTrackMult->GetMeanError();
243}
244
245//-----------------------------------------------------------------------------
246void QAFillNPi0()
247{
3b8b9fc1 248 TH1 *hTotSelEvents = (TH1*)listHist->FindObject("hTotSelEvents");
249 Int_t nEvents4 = hTotSelEvents->GetBinContent(4);
250
251 if( nEvents4 < 10000 ) {
252 Printf(" -> skipping due to small number of selected events: %d", nEvents4);
253 return;
254 }
255
256
a600ad2c 257 TCanvas* canv = new TCanvas;
258 canv->Divide(3,3);
259 canv->cd(1);
260
a600ad2c 261 TH2F *hReMassPt = (TH2F*)listHist->FindObject("hPi0All_cen0");
262 TH2 *hMiMassPt = (TH2*)listHist->FindObject("hMiPi0All_cen0");
263 TH1* hReMass = Pi0InvMass(hReMassPt,2.,10.);
264 TH1* hMiMass = Pi0InvMass(hMiMassPt,2.,10.);
265 hReMass->Rebin(2);
266 hMiMass->Rebin(2);
267 hReMass->SetAxisRange(0.,0.3);
268 hMiMass->SetAxisRange(0.,0.3);
269
270 // Draw
271 hReMassPt->DrawCopy("colz");
272 canv->cd(2);
3b8b9fc1 273 //hReMass->SetAxisRange(0.1, 0.2);
274 hReMass->DrawCopy("E");
a600ad2c 275 canv->cd(3);
276
277
278 TH1* hReMiRatio = (TH1*)hReMass->Clone("hReMiRatio");
279 TH1* hPi0SubBG = (TH1*)hReMass->Clone("hPi0SubBG");
280 hReMiRatio->Divide(hMiMass);
3b8b9fc1 281 hReMiRatio->DrawCopy("E");
a600ad2c 282 canv->cd(4);
283
284
285 TF1 * fitM = new TF1("fitM",pi0massP2,0.,1.,6) ;
286 fitM->SetLineColor(kRed);
287 fitM->SetParName(0,"A") ;
288 fitM->SetParName(1,"m_{0}") ;
289 fitM->SetParName(2,"#sigma") ;
290 fitM->SetParName(3,"a_{0}") ;
291 fitM->SetParName(4,"a_{1}") ;
292 fitM->SetParName(5,"a_{2}") ;
293
294 TF1 * fitG = new TF1("fitG",pi0mass,0.,1.,3) ;
295 fitG->SetLineColor(kRed);
296 fitG->SetParName(0,"A") ;
297 fitG->SetParName(1,"m_{0}") ;
298 fitG->SetParName(2,"#sigma") ;
299
300 TF1 * fitBG = new TF1("fitBG",bgP2,0.,1.,3) ;
301 fitBG->SetLineColor(kBlue);
302 fitBG->SetLineStyle(2);
303 fitBG->SetParName(0,"a_{0}") ;
304 fitBG->SetParName(1,"a_{1}") ;
305 fitBG->SetParName(2,"a_{2}") ;
306
3b8b9fc1 307 fitM->SetParameters(0.003, 0.134, 0.007, 0.14, -0.02, -0.01) ;
308 fitM->SetParLimits(0, 0., 1. ) ;
309 fitM->SetParLimits(1, 0.1, 0.2 ) ;
310 fitM->SetParLimits(2, 0., 0.05) ;
a600ad2c 311
312 Double_t rangeMin=0.06 ;
313 Double_t rangeMax=0.25 ;
3b8b9fc1 314 TFitResultPtr mrp = hReMiRatio->Fit(fitM,"Q","",rangeMin,rangeMax) ;
315 //mrp = hReMiRatio->Fit(fitM,"MQ","",rangeMin,rangeMax) ; // "M" always fail...
316 int error = mrp;
317 if( error % 1000) {
318 Printf(" -> fit of fitM to hReMiRatio failed with error code %d", error % 1000);
319 continue;
320 }
321 else if( error )
322 Printf("Warning: failure of 'improve result' of fit of fitM to hReMiRatio");
a600ad2c 323 fitBG->SetParameters(fitM->GetParameter(3),
324 fitM->GetParameter(4),
325 fitM->GetParameter(5));
3b8b9fc1 326 hReMiRatio->SetAxisRange(rangeMin, rangeMax);
327 hReMiRatio->DrawCopy();
328 canv->cd(5);
329
330
a600ad2c 331 hMiMass->Multiply(fitBG) ;
332 hPi0SubBG ->Add(hMiMass,-1.) ;
333
334
335 fitG->SetParameters(100.,fitM->GetParameter(1),fitM->GetParameter(2)) ;
3b8b9fc1 336 // fitG->SetParLimits(0,0.000,1.e+5) ;
337 fitG->SetParLimits(1, rangeMin, rangeMax) ;
338 fitG->SetParLimits(2, 0., rangeMax) ;
339 mrp = hPi0SubBG->Fit(fitG,"Q","",rangeMin,rangeMax);
340 if( (error=mrp) ) {
341 Printf(" -> fit of fitG to hPi0SubBG failed with error code %d, skipping", error );
342 continue;
343 }
344 hPi0SubBG->SetAxisRange(rangeMin, rangeMax);
345 hPi0SubBG->DrawCopy();
346
347
348
a600ad2c 349 Double_t pi0Peak = fitG->GetParameter(1);
350 Double_t pi0Sigma = fitG->GetParameter(2);
351 Double_t epi0Peak = fitG->GetParError(1);
352 Double_t epi0Sigma = fitG->GetParError(2);
353 Int_t iMin = hPi0SubBG->FindBin(pi0Peak - 3*pi0Sigma);
354 Int_t iMax = hPi0SubBG->FindBin(pi0Peak + 3*pi0Sigma);
355
356 Double_t nPi0,ePi0;
357
3b8b9fc1 358
359 nPi0 = hPi0SubBG->IntegralAndError(iMin, iMax, ePi0);
360 //ePi0 = TMath::Sqrt(nPi0);
a600ad2c 361 // Fit1Pi0(hMass,2,0.05,0.20,nPi0,ePi0);
362
a600ad2c 363 mPi0 [runIndex] = pi0Peak;
364 emPi0 [runIndex] = epi0Peak;
365 wPi0 [runIndex] = pi0Sigma;
366 ewPi0 [runIndex] = epi0Sigma;
b5bff3e2 367 nPi0Event[runIndex] = nPi0/nEvents4;
368 ePi0Event[runIndex] = ePi0/nEvents4;
a600ad2c 369}
370
371//-----------------------------------------------------------------------------
372void QAWriteEventSelection()
373{
ff329095 374 TH1F *hNSelected = new TH1F("hNSelected","N_{selected}",runIndex,0,runIndex);
a600ad2c 375 TH1F *grVtx = new TH1F("grVtx","N_{vertex exists} / N_{total}",runIndex,0,runIndex);
376 TH1F *grVtxZ10 = new TH1F("grVtxZ10","N_{vertex exists, |Z|<10 cm} / N_{total}",runIndex,0,runIndex);
377 TH1F *grVtxZ10Cent = new TH1F("grVtxZ10Cent","N_{vertex exists, |Z|<10 cm, centrality} / N_{total}",runIndex,0,runIndex);
378
379 for (Int_t i=0; i<runIndex; i++) {
ff329095 380 hNSelected->SetBinContent(i+1, nTotal4[i]);
381
a600ad2c 382 grVtx ->SetBinContent(i+1,rVtx[i]);
383 grVtxZ10 ->SetBinContent(i+1,rVtxZ10[i]);
384 grVtxZ10Cent ->SetBinContent(i+1,rVtxZ10Cent[i]);
385
386 grVtx ->SetBinError(i+1,eVtx[i]);
387 grVtxZ10 ->SetBinError(i+1,eVtxZ10[i]);
388 grVtxZ10Cent ->SetBinError(i+1,eVtxZ10Cent[i]);
389
390 grVtx ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
391 grVtxZ10 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
392 grVtxZ10Cent ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
393 }
394
395 grVtx ->LabelsOption("v");
396 grVtxZ10 ->LabelsOption("v");
397 grVtxZ10Cent ->LabelsOption("v");
398
b5bff3e2 399 grVtx ->SetMarkerStyle(33);
400 grVtxZ10 ->SetMarkerStyle(33);
401 grVtxZ10Cent ->SetMarkerStyle(33);
a600ad2c 402
ff329095 403 hNSelected ->Write();
a600ad2c 404 grVtx ->Write();
405 grVtxZ10 ->Write();
406 grVtxZ10Cent ->Write();
407}
408
409//-----------------------------------------------------------------------------
410void QAWriteOccupancy()
411{
412 // Number of fired cells in each module
413
414 TH1F *grNCellsM1 = new TH1F("grNCellsM1","N_{cells} in module 1",runIndex,0,runIndex);
415 TH1F *grNCellsM2 = new TH1F("grNCellsM2","N_{cells} in module 2",runIndex,0,runIndex);
416 TH1F *grNCellsM3 = new TH1F("grNCellsM3","N_{cells} in module 3",runIndex,0,runIndex);
417
418 for (Int_t i=0; i<runIndex; i++) {
419 grNCellsM1 ->SetBinContent(i+1,nCellsM1[i]);
420 grNCellsM2 ->SetBinContent(i+1,nCellsM2[i]);
421 grNCellsM3 ->SetBinContent(i+1,nCellsM3[i]);
422
423 grNCellsM1 ->SetBinError(i+1,0);
424 grNCellsM2 ->SetBinError(i+1,0);
425 grNCellsM3 ->SetBinError(i+1,0);
426
427 grNCellsM1 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
428 grNCellsM2 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
429 grNCellsM3 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
430
431 grNCellsM1->GetYaxis()->SetTitle("N_{cells}");
432 grNCellsM2->GetYaxis()->SetTitle("N_{cells}");
433 grNCellsM3->GetYaxis()->SetTitle("N_{cells}");
434 }
435
436 grNCellsM1 ->LabelsOption("v");
437 grNCellsM2 ->LabelsOption("v");
438 grNCellsM3 ->LabelsOption("v");
b5bff3e2 439 grNCellsM1 ->SetMarkerStyle(33);
440 grNCellsM2 ->SetMarkerStyle(33);
441 grNCellsM3 ->SetMarkerStyle(33);
a600ad2c 442 grNCellsM1 ->Write();
443 grNCellsM2 ->Write();
444 grNCellsM3 ->Write();
445}
446
447//-----------------------------------------------------------------------------
448void QAWriteClusters()
449{
450 TH1F *grECluster = new TH1F("grECluster","#LTE_{cluster}#GT",runIndex,0,runIndex);
451 for (Int_t i=0; i<runIndex; i++) {
452 grECluster ->SetBinContent(i+1,cluEnergy[i]);
453 grECluster ->SetBinError(i+1,ecluEnergy[i]);
454 grECluster ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
455 }
456
457 grECluster ->LabelsOption("v");
b5bff3e2 458 grECluster ->SetMarkerStyle(33);
a600ad2c 459 grECluster ->Write();
460
461 TH1F *grNCluster = new TH1F("grNCluster","#LTN_{clusters}#GT",runIndex,0,runIndex);
462 for (Int_t i=0; i<runIndex; i++) {
463 grNCluster ->SetBinContent(i+1,nCluEvent[i]);
464 grNCluster ->SetBinError(i+1,eCluEvent[i]);
465 grNCluster ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
466 }
467 grNCluster ->LabelsOption("v");
b5bff3e2 468 grNCluster ->SetMarkerStyle(33);
a600ad2c 469 grNCluster ->Write();
b5bff3e2 470
471 for(int cent=0; cent<kNCents; ++cent) {
472 for(int ipid = 0; ipid < kNPID; ++ipid) {
473 TH1* hPhot = listHist->FindObject( Form("hPhot%s_cen%d", kPIDNames[ipid], cent) );
ff329095 474 AddWriteTH1F(Form("grNPhot%s_cen%d", kPIDNames[ipid], cent), Form("#LTN_{clusters}^{%s}#GT, c.bin=%d", kPIDNames[ipid], cent),
b5bff3e2 475 nPhotPID[cent][ipid], enPhotPID[cent][ipid]);
ff329095 476 AddWriteTH1F(Form("grEn%s_cen%d", kPIDNames[ipid], cent), Form("#LTE_{clusters}^{%s}#GT, c.bin=%d", kPIDNames[ipid], cent),
477 mEnPID[cent][ipid], emEnPID[cent][ipid]);
478 AddWriteTH1F(Form("grNPhot%sHigh_cen%d", kPIDNames[ipid], cent), Form("#LTN_{clusters}^{%s}#GT, c.bin=%d", kPIDNames[ipid], cent),
479 nPhotPIDHigh[cent][ipid], enPhotPIDHigh[cent][ipid]);
480 AddWriteTH1F(Form("grEn%sHigh_cen%d", kPIDNames[ipid], cent), Form("#LTE_{clusters}^{%s}#GT, c.bin=%d", kPIDNames[ipid], cent),
481 mEnPIDHigh[cent][ipid], emEnPIDHigh[cent][ipid]);
b5bff3e2 482 }
483 }
a600ad2c 484}
485
486//-----------------------------------------------------------------------------
487void QAWriteTracks()
488{
489 TH1F *grNTracks0 = new TH1F("grNTracks0","#LTN_{tracks}#GT",runIndex,0,runIndex);
490 for (Int_t i=0; i<runIndex; i++) {
491 grNTracks0 ->SetBinContent(i+1,nTracks0[i]);
492 grNTracks0 ->SetBinError(i+1,eTracks0[i]);
493 grNTracks0 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
494 }
495 grNTracks0 ->LabelsOption("v");
b5bff3e2 496 grNTracks0 ->SetMarkerStyle(33);
a600ad2c 497 grNTracks0 ->Write();
498
499// TH1F *grNTracks1 = new TH1F("grNTracks1","#LTN_{tracks}#GT, N #geq 5",runIndex,0,runIndex);
500// for (Int_t i=0; i<runIndex; i++) {
501// grNTracks1 ->SetBinContent(i+1,nTracks0[i]);
502// grNTracks1 ->SetBinError(i+1,eTracks0[i]);
503// grNTracks1 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
504// }
505// grNTracks1 ->LabelsOption("v");
b5bff3e2 506// grNTracks1 ->SetMarkerStyle(33);
a600ad2c 507// grNTracks1 ->Write();
508}
509
510//-----------------------------------------------------------------------------
511void QAWriteNPi0()
512{
513 TH1F *grNPi0 = new TH1F("grNPi0","N(#pi^{0}) per event",runIndex,0,runIndex);
514 for (Int_t i=0; i<runIndex; i++) {
515 grNPi0 ->SetBinContent(i+1,nPi0Event[i]);
516 grNPi0 ->SetBinError(i+1,ePi0Event[i]);
517 grNPi0 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
518 }
519 grNPi0 ->LabelsOption("v");
b5bff3e2 520 grNPi0 ->SetMarkerStyle(33);
5f3e9080 521 grNPi0->GetYaxis()->SetTitle("#LTN#GT");
a600ad2c 522 grNPi0 ->Write();
523
524 TH1F *grMPi0 = new TH1F("grMPi0","M(#pi^{0})",runIndex,0,runIndex);
525 for (Int_t i=0; i<runIndex; i++) {
526 grMPi0 ->SetBinContent(i+1,mPi0[i]);
527 grMPi0 ->SetBinError(i+1,emPi0[i]);
528 grMPi0 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
529 }
530 grMPi0 ->LabelsOption("v");
b5bff3e2 531 grMPi0 ->SetMarkerStyle(33);
5f3e9080 532 grMPi0->GetYaxis()->SetTitle("#LTM#GT");
a600ad2c 533 grMPi0 ->Write();
534
535 TH1F *grWPi0 = new TH1F("grWPi0","#sigma(#pi^{0})",runIndex,0,runIndex);
536 for (Int_t i=0; i<runIndex; i++) {
537 grWPi0 ->SetBinContent(i+1,wPi0[i]);
538 grWPi0 ->SetBinError(i+1,ewPi0[i]);
539 grWPi0 ->GetXaxis()->SetBinLabel(i+1,Form("%d",runNum[i]));
540 }
541 grWPi0 ->LabelsOption("v");
b5bff3e2 542 grWPi0 ->SetMarkerStyle(33);
5f3e9080 543 grWPi0->GetYaxis()->SetTitle("#LT#sigma#GT");
a600ad2c 544 grWPi0 ->Write();
545}
546
547//-----------------------------------------------------------------------------
548TH1D * Pi0InvMass(TH2 *hMassPt, Float_t ptmin, Float_t ptmax)
549{
550 Int_t bmin = hMassPt->GetYaxis()->FindBin(ptmin);
551 Int_t bmax = hMassPt->GetYaxis()->FindBin(ptmax);
552
553 TString name = hMassPt->GetName();
554 name += "_M";
555 TH1D* hMass = hMassPt->ProjectionX(name,bmin,bmax);
556 char htitle[64];
557 sprintf(htitle,"%.1f<p_{T}<%.1f GeV/c",ptmin,ptmax);
558 hMass->SetTitle(htitle);
559
3b8b9fc1 560 hMass->Sumw2();
561
a600ad2c 562 return hMass;
563
564}
565
566//-----------------------------------------------------------------------------
b5bff3e2 567TH1F *AddWriteTH1F(const char *name, const char *title, Double_t* values, Double_t* errors)
568{
569 TH1F* hist = new TH1F(name, title, runIndex, 0, runIndex);
570 hist->SetContent(values);
571 hist->SetError(errors);
572
573 for(Int_t i=0; i<runIndex; i++) {
574 hist->GetXaxis()->SetBinLabel( i+1, Form("%d",runNum[i]) );
575 }
576
577 hist->LabelsOption("v");
578 hist->SetMarkerStyle(33);
579 hist->Write();
580
5f3e9080 581 if( TString(name).Contains("En") )
582 hist->GetYaxis()->SetTitle("#LTE#GT");
583 if( TString(name).Contains("NPhot") )
584 hist->GetYaxis()->SetTitle("#LTN#GT");
585
b5bff3e2 586 return hist;
587}
588
589//-----------------------------------------------------------------------------
a600ad2c 590void Fit1Pi0(TH1D *hMass,
591 const Int_t polN,
592 const Double_t mFitMin,
593 const Double_t mFitMax,
594 Double_t &nPi0, Double_t &ePi0)
595{
596 // This script takes a 2D histogram hMassPt with invariant mass vs
597 // pt of cluster pairs:
598 // - slices them along pt bins,
599 // - fits 1D invariant mass specrta by Gaussian+polynomial
600 // - calculates the number of pi0's as an integral of Gaussian
601 // - calculates background under pi0 as an integral of polynomial
602 // Author: Yuri Kharlov.
603
604 TF1 *fitfun = 0;
605 if (polN == 1)
606 fitfun = new TF1("fitfun",pi0massP1,0.1,0.7,5);
607 else if (polN == 2)
608 fitfun = new TF1("fitfun",pi0massP2,0.1,0.7,6);
609 fitfun->SetLineColor(kRed);
610 fitfun->SetLineWidth(2);
611
612 fitfun->SetParName(0,"A");
613 fitfun->SetParName(1,"m_{0}");
614 fitfun->SetParName(2,"#sigma");
615 fitfun->SetParName(3,"a_{0}");
616 fitfun->SetParName(4,"a_{1}");
617 if (polN == 2)
618 fitfun->SetParName(5,"a_{2}");
619
620 fitfun->SetParLimits(0, 1.000,1.e+5);
621 fitfun->SetParLimits(1, 0.09,0.18);
622 fitfun->SetParLimits(2, 0.003,0.020);
623
624 Int_t nM = hMass->GetNbinsX();
625 Float_t mMin = hMass->GetXaxis()->GetXmin();
626 Float_t mMax = hMass->GetXaxis()->GetXmax();
627 Float_t mStep = (mMax-mMin)/nM;
628
629 hMass->SetXTitle("M_{#gamma#gamma}, GeV/c");
630 hMass->SetAxisRange(0.01,1.0);
631 Float_t nPi0Max = hMass->Integral(hMass->FindBin(0.30),
632 hMass->FindBin(0.40));
633 for (Int_t iM=1; iM<nM; iM++)
634 if (hMass->GetBinContent(iM) == 0) hMass->SetBinError(iM,1.0);
635 hMass->SetMinimum(0.01);
636 if (polN == 1)
637 fitfun->SetParameters(nPi0Max/4,0.135,0.010,1.,0.);
638 else if (polN == 2)
639 fitfun->SetParameters(nPi0Max/4,0.135,0.010,1.,0.,0.);
640 hMass->Fit("fitfun","Q0","",mFitMin,mFitMax);
641
642 nPi0 = fitfun->GetParameter(0)*fitfun->GetParameter(2)/
643 mStep*TMath::Sqrt(TMath::TwoPi());
644 ePi0 = TMath::Sqrt(nPi0);
645}
646
647
648//-----------------------------------------------------------------------------
649Double_t pi0massP2(Double_t *x, Double_t *par)
650{
651 Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
652 (2*par[2]*par[2]) );
653 Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
654 return gaus+back;
655}
656
657//-----------------------------------------------------------------------------
658Double_t pi0massP1(Double_t *x, Double_t *par)
659{
660 Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
661 (2*par[2]*par[2]) );
662 Double_t back = par[3] + par[4]*x[0];
663 return gaus+back;
664}
665
666//-----------------------------------------------------------------------------
667Double_t pi0mass(Double_t *x, Double_t *par)
668{
669 Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
670 (2*par[2]*par[2]) );
671 return gaus;
672}
673
674//-----------------------------------------------------------------------------
675Double_t bgP2(Double_t *x, Double_t *par)
676{
677 Double_t back = par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
678 return back;
679}