PWGGA/PHOSTasks/PHOS_PbPb/macros/QA/ExtractQA.C

   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"
  15 using namespace std;
  16 #endif
  17
  18 //-----------------------------------------------------------------------------
  19 // Function declaration
  20 void QAFillEventSelection();
  21 void QAFillOccupancy();
  22 void QAFillClusters();
  23 void QAFillTracks();
  24 void QAFillNPi0();
  25
  26 void QAWriteEventSelection();
  27 void QAWriteOccupancy();
  28 void QAWriteClusters();
  29 void QAWriteTracks();
  30 void QAWriteNPi0();
  31
  32 TH1D *Pi0InvMass(TH2 *hMassPt, Float_t ptmin, Float_t ptmax);
  33 TH1F *AddWriteTH1F(const char *name, const char *title, Double_t* values, Double_t* errors);
  34 void 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);
  39 Double_t pi0massP1(Double_t *x, Double_t *par);
  40 Double_t pi0massP2(Double_t *x, Double_t *par);
  41 Double_t pi0mass  (Double_t *x, Double_t *par);
  42 Double_t bgP2     (Double_t *x, Double_t *par);
  43
  44 //-----------------------------------------------------------------------------
  45 // Global variabes
  46 const Int_t kNCents = 1;
  47 const Int_t kNPID = 8+4;
  48 const char* kPIDNames[kNPID] = {"All", "Allwou", "Disp", "Disp2", "Dispwou", "CPV", "CPV2", "Both",
  49                                 "Allcore", "Dispcore", "CPVcore", "Bothcore"};
  50
  51 Int_t runIndex;
  52 Int_t runNum[500];
  53 TList *listHist;
  54 Double_t run[500], erun[500];
  55 // Event selection:
  56 Double_t nTotal4[500];
  57 Double_t rVtx[500], rVtxZ10[500], rVtxZ10Cent[500];
  58 Double_t eVtx[500], eVtxZ10[500], eVtxZ10Cent[500];
  59 // Occupancy:
  60 Double_t nCellsM1[500], nCellsM2[500], nCellsM3[500];
  61 // Clusters:
  62 Double_t nCluEvent[500], eCluEvent[500];
  63 Double_t cluEnergy[500], ecluEnergy[500];
  64 Double_t nPhotPID[kNCents][kNPID][500], enPhotPID[kNCents][kNPID][500];
  65 Double_t mEnPID[kNCents][kNPID][500],   emEnPID[kNCents][kNPID][500];
  66 Double_t nPhotPIDHigh[kNCents][kNPID][500], enPhotPIDHigh[kNCents][kNPID][500];
  67 Double_t mEnPIDHigh[kNCents][kNPID][500],   emEnPIDHigh[kNCents][kNPID][500];
  68
  69
  70 // Tracks:
  71 Double_t nTracks0[500], eTracks0[500];
  72 Double_t nTracks1[500], eTracks1[500];
  73 // Pi0:
  74 Double_t mPi0[500], emPi0[500];
  75 Double_t wPi0[500], ewPi0[500];
  76 Double_t nPi0Event[500], ePi0Event[500];
  77
  78
  79 //-----------------------------------------------------------------------------
  80 void ExtractQA(const TString filelist="filelist.txt",
  81                const TString runlist="runlist.txt",
  82                const TString outputFileName = "outputQA.root")
  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();
 121     QAFillNPi0();
 122
 123     listHist->Clear();
 124     rootFile->Close();
 125     runIndex++;
 126   }
 127
 128
 129   TFile *fileQA = TFile::Open(outputFileName.Data(), "recreate");
 130   QAWriteEventSelection();
 131   QAWriteOccupancy();
 132   QAWriteClusters();
 133   QAWriteTracks();
 134   QAWriteNPi0();
 135   fileQA         ->Close();
 136
 137 }
 138
 139 //-----------------------------------------------------------------------------
 140 void 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
 152   nTotal4[runIndex] = hSelEvents->GetBinContent(4);
 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 //-----------------------------------------------------------------------------
 167 void 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 //-----------------------------------------------------------------------------
 189 void QAFillClusters()
 190 {
 191   TH1F *hTotSelEvents = (TH1F*)listHist->FindObject("hTotSelEvents");
 192   Double_t nEvents4 = hTotSelEvents->GetBinContent(4);
 193   if( ! nEvents4 ) return;
 194
 195   TH2F *hCluEvsClu = (TH2F*)listHist->FindObject("hCluEvsClu");
 196   TH1D *hClusterE = hCluEvsClu->ProjectionX("hClusterE",4,41);
 197
 198   cluEnergy[runIndex]  = hClusterE->GetMean();
 199   ecluEnergy[runIndex] = hClusterE->GetMeanError();
 200
 201   nCluEvent[runIndex]  = hClusterE->Integral() / nEvents4;
 202   eCluEvent[runIndex]  = TMath::Sqrt( hClusterE->Integral() )/nEvents4;
 203
 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
 209       hPhot->SetAxisRange(0., 100.);
 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();
 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
 223     }
 224   }
 225 }
 226
 227 //-----------------------------------------------------------------------------
 228 void 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 //-----------------------------------------------------------------------------
 246 void QAFillNPi0()
 247 {
 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
 257   TCanvas* canv = new TCanvas;
 258   canv->Divide(3,3);
 259   canv->cd(1);
 260
 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);
 273   //hReMass->SetAxisRange(0.1, 0.2);
 274   hReMass->DrawCopy("E");
 275   canv->cd(3);
 276
 277
 278   TH1* hReMiRatio = (TH1*)hReMass->Clone("hReMiRatio");
 279   TH1* hPi0SubBG  = (TH1*)hReMass->Clone("hPi0SubBG");
 280   hReMiRatio->Divide(hMiMass);
 281   hReMiRatio->DrawCopy("E");
 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
 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) ;
 311
 312   Double_t rangeMin=0.06 ;
 313   Double_t rangeMax=0.25 ;
 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");
 323   fitBG->SetParameters(fitM->GetParameter(3),
 324                        fitM->GetParameter(4),
 325                        fitM->GetParameter(5));
 326   hReMiRatio->SetAxisRange(rangeMin, rangeMax);
 327   hReMiRatio->DrawCopy();
 328   canv->cd(5);
 329
 330
 331   hMiMass->Multiply(fitBG) ;
 332   hPi0SubBG ->Add(hMiMass,-1.) ;
 333
 334
 335   fitG->SetParameters(100.,fitM->GetParameter(1),fitM->GetParameter(2)) ;
 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
 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
 358
 359   nPi0 = hPi0SubBG->IntegralAndError(iMin, iMax, ePi0);
 360   //ePi0 = TMath::Sqrt(nPi0);
 361   // Fit1Pi0(hMass,2,0.05,0.20,nPi0,ePi0);
 362
 363   mPi0     [runIndex] = pi0Peak;
 364   emPi0    [runIndex] = epi0Peak;
 365   wPi0     [runIndex] = pi0Sigma;
 366   ewPi0    [runIndex] = epi0Sigma;
 367   nPi0Event[runIndex] = nPi0/nEvents4;
 368   ePi0Event[runIndex] = ePi0/nEvents4;
 369 }
 370
 371 //-----------------------------------------------------------------------------
 372 void QAWriteEventSelection()
 373 {
 374   TH1F *hNSelected       = new TH1F("hNSelected","N_{selected}",runIndex,0,runIndex);
 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++) {
 380     hNSelected->SetBinContent(i+1, nTotal4[i]);
 381
 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
 399   grVtx          ->SetMarkerStyle(33);
 400   grVtxZ10       ->SetMarkerStyle(33);
 401   grVtxZ10Cent   ->SetMarkerStyle(33);
 402
 403   hNSelected     ->Write();
 404   grVtx          ->Write();
 405   grVtxZ10       ->Write();
 406   grVtxZ10Cent   ->Write();
 407 }
 408
 409 //-----------------------------------------------------------------------------
 410 void 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");
 439   grNCellsM1     ->SetMarkerStyle(33);
 440   grNCellsM2     ->SetMarkerStyle(33);
 441   grNCellsM3     ->SetMarkerStyle(33);
 442   grNCellsM1     ->Write();
 443   grNCellsM2     ->Write();
 444   grNCellsM3     ->Write();
 445 }
 446
 447 //-----------------------------------------------------------------------------
 448 void 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");
 458   grECluster ->SetMarkerStyle(33);
 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");
 468   grNCluster ->SetMarkerStyle(33);
 469   grNCluster ->Write();
 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) );
 474       AddWriteTH1F(Form("grNPhot%s_cen%d", kPIDNames[ipid], cent), Form("#LTN_{clusters}^{%s}#GT, c.bin=%d", kPIDNames[ipid], cent),
 475                    nPhotPID[cent][ipid], enPhotPID[cent][ipid]);
 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]);
 482     }
 483   }
 484 }
 485
 486 //-----------------------------------------------------------------------------
 487 void 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");
 496   grNTracks0 ->SetMarkerStyle(33);
 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");
 506 //   grNTracks1 ->SetMarkerStyle(33);
 507 //   grNTracks1 ->Write();
 508 }
 509
 510 //-----------------------------------------------------------------------------
 511 void 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");
 520   grNPi0     ->SetMarkerStyle(33);
 521   grNPi0->GetYaxis()->SetTitle("#LTN#GT");
 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");
 531   grMPi0     ->SetMarkerStyle(33);
 532   grMPi0->GetYaxis()->SetTitle("#LTM#GT");
 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");
 542   grWPi0     ->SetMarkerStyle(33);
 543   grWPi0->GetYaxis()->SetTitle("#LT#sigma#GT");
 544   grWPi0     ->Write();
 545 }
 546
 547 //-----------------------------------------------------------------------------
 548 TH1D * 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
 560   hMass->Sumw2();
 561
 562   return hMass;
 563
 564 }
 565
 566 //-----------------------------------------------------------------------------
 567 TH1F *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
 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
 586   return hist;
 587 }
 588
 589 //-----------------------------------------------------------------------------
 590 void 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 //-----------------------------------------------------------------------------
 649 Double_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 //-----------------------------------------------------------------------------
 658 Double_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 //-----------------------------------------------------------------------------
 667 Double_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 //-----------------------------------------------------------------------------
 675 Double_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 }