PWGGA/PHOSTasks/PHOS_PbPb/macros/DrawPi0Flow.C

   1 DrawPi0Flow(const TString filename = "AnalysisResults.root",
   2             const TString eventType = "SemiCentral")
   3 {
   4   /* $Id$ */
   5
   6   gStyle->SetOptStat(0);
   7   gStyle->SetPadGridX(1);
   8   gStyle->SetPadGridY(1);
   9   TFile * f = new TFile(filename) ;
  10   TString listName;
  11   if (eventType.Contains("Central"))
  12     listName = Form("AliPHOSPi0Flow/PHOSPi0FlowCoutput1",
  13                     eventType.Data(),eventType.Data());
  14   else
  15     listName = Form("AliPHOSPi0Flow_%s/AliPHOSPi0Flow_%sCoutput1",
  16                     eventType.Data(),eventType.Data());
  17   cout << listName << endl;
  18   TList *histoList = (TList*)f->Get(listName.Data()); // lego train
  19
  20   TH1F *hev         = (TH1F*)histoList->FindObject("hTotSelEvents") ;
  21   TH2F* hZvertex    = (TH2F*)histoList->FindObject("hZvertex");
  22   TH2F* hCenTrack   = (TH2F*)histoList->FindObject("hCenTrack");
  23   TH2F *hCenPHOS    = (TH2F*)histoList->FindObject("hCenPHOS") ;
  24   TH2F* phiRP       = (TH2F*)histoList->FindObject("phiRP");
  25   TH2F* phiRPflat   = (TH2F*)histoList->FindObject("phiRPflat");
  26   TH2F* phiRPV0A    = (TH2F*)histoList->FindObject("phiRPV0A");
  27   TH2F* phiRPV0C    = (TH2F*)histoList->FindObject("phiRPV0C");
  28   TH2F* phiRPV0Aflat= (TH2F*)histoList->FindObject("phiRPV0Aflat");
  29   TH2F* phiRPV0Cflat= (TH2F*)histoList->FindObject("phiRPV0Cflat");
  30   TH2F* hCellNXZM1  = (TH2F*)histoList->FindObject("hCellNXZM1");
  31   TH2F* hCellNXZM2  = (TH2F*)histoList->FindObject("hCellNXZM2");
  32   TH2F* hCellNXZM3  = (TH2F*)histoList->FindObject("hCellNXZM3");
  33   TH2F *hPi0All_cen0     = (TH2F*)histoList->FindObject("hPi0All_cen0");
  34   TH2F *hPi0Allcore_cen0 = (TH2F*)histoList->FindObject("hPi0Allcore_cen0");
  35   TH2F *hPi0M11          = (TH2F*)histoList->FindObject("hPi0M11");
  36   TH2F *hPi0M22          = (TH2F*)histoList->FindObject("hPi0M22");
  37   TH2F *hPi0M33          = (TH2F*)histoList->FindObject("hPi0M33");
  38
  39   //-----------------------------------------------------------------------------
  40   TCanvas *c1 = new TCanvas("c1","Event selection");
  41   hev->GetXaxis()->SetRangeUser(0,6);
  42   hev->GetXaxis()->SetBinLabel(1,"Total");
  43   hev->GetXaxis()->SetBinLabel(2,"Has Vertex");
  44   hev->GetXaxis()->SetBinLabel(3,"abs(z_vertex) < 10.");
  45   hev->GetXaxis()->SetBinLabel(4,"Has Centrality");
  46   hev->GetXaxis()->SetBinLabel(5,"C. upper edge");
  47   hev->GetXaxis()->SetBinLabel(6,"C. lower edge");
  48   hev->GetXaxis()->SetBinLabel(7,"Has PClusters");
  49   hev->SetYTitle("N_{events}");
  50   hev->GetYaxis()->SetTitleOffset(1.2);
  51   hev->Draw();
  52   c1->Print("PHOS_EvSel.eps");
  53
  54   //-----------------------------------------------------------------------------
  55   TCanvas *c2 = new TCanvas("c2","Track multiplicity vs centrality");
  56   gPad->SetLogz();
  57   hCenTrack->SetXTitle("centrality (%)");
  58   hCenTrack->SetYTitle("Number of tracks");
  59   hCenTrack->Draw("colz");
  60   c2->Print("TrackMultCentrality.eps");
  61
  62   //-----------------------------------------------------------------------------
  63   TCanvas *c3 = new TCanvas("c3","PHOS multiplicity vs centrality");
  64   gPad->SetLogz();
  65   hCenPHOS->SetXTitle("centrality (%)");
  66   hCenPHOS->SetYTitle("Number of PHOS clusters");
  67   hCenPHOS->Draw("colz");
  68   c3->Print("PHOSMultCentrality.eps");
  69
  70   //-----------------------------------------------------------------------------
  71   TCanvas *c3cent = new TCanvas("c3cent","Centrality");
  72   hCenPHOS->SetXTitle("centrality (%)");
  73   hCenPHOS->SetYTitle("Number of PHOS clusters");
  74   TH1D *hCent = hCenPHOS->ProjectionX();
  75   hCent->SetTitle("Centrality");
  76   hCent->Draw();
  77   c3cent->Print("Centrality.eps");
  78
  79   //-----------------------------------------------------------------------------
  80   TCanvas *c4V0 = new TCanvas("c4V0","RP phi in VZERO");
  81   phiRPV0A1 = phiRPV0A->ProjectionX();
  82   phiRPV0C1 = phiRPV0C->ProjectionX();
  83   phiRPV0A1->SetLineColor(kRed);
  84   phiRPV0C1->SetLineColor(kBlue);
  85   phiRPV0A1->SetXTitle("VZERO RP #phi");
  86   phiRPV0C1->SetXTitle("VZERO RP #phi");
  87   phiRPV0A1->SetMinimum(phiRPV0A1->GetMinimum()*0.98);
  88   phiRPV0A1->SetMaximum(phiRPV0A1->GetMaximum()*1.02);
  89   phiRPV0A1->Draw();
  90   phiRPV0C1->Draw("same");
  91   leg = new TLegend(0.2,0.8,0.39,0.89);
  92   leg->SetFillColor(kWhite);
  93   leg->SetBorderSize(1);
  94   leg->AddEntry(phiRPV0A1,"V0A","l");
  95   leg->AddEntry(phiRPV0C1,"V0C","l");
  96   leg->Draw();
  97   c4V0->Print("V0RPphi.eps");
  98
  99   //-----------------------------------------------------------------------------
 100   TCanvas *c4V0flat = new TCanvas("c4V0flat","RP phi in VZERO flattened");
 101   phiRPV0Aflat1 = phiRPV0Aflat->ProjectionX();
 102   phiRPV0Cflat1 = phiRPV0Cflat->ProjectionX();
 103   phiRPV0Aflat1->SetLineColor(kRed);
 104   phiRPV0Cflat1->SetLineColor(kBlue);
 105   phiRPV0Aflat1->SetXTitle("VZERO RP #phi");
 106   phiRPV0Cflat1->SetXTitle("VZERO RP #phi");
 107   phiRPV0Aflat1->Draw();
 108   phiRPV0Cflat1->Draw("same");
 109   leg = new TLegend(0.2,0.8,0.39,0.89);
 110   leg->SetFillColor(kWhite);
 111   leg->SetBorderSize(1);
 112   leg->AddEntry(phiRPV0Aflat1,"V0A","l");
 113   leg->AddEntry(phiRPV0Cflat1,"V0C","l");
 114   leg->Draw();
 115   c4V0flat->Print("V0RPphiFlat.eps");
 116
 117   //-----------------------------------------------------------------------------
 118   TCanvas *c4TPC = new TCanvas("c4TPC","TPC RP phi");
 119   phiRP1 = phiRP->ProjectionX();
 120   phiRP1->SetXTitle("TPC RP #phi");
 121   phiRP1->Draw();
 122   c4TPC->Print("TPCRPphi.eps");
 123
 124   //-----------------------------------------------------------------------------
 125   TCanvas *c4TPCflat = new TCanvas("c4TPCflat","TPC RP phi flattened");
 126   phiRPflat1 = phiRPflat->ProjectionX();
 127   phiRPflat1->SetXTitle("TPC RP #phi");
 128   phiRPflat1->Draw();
 129   c4TPCflat->Print("TPCRPphiFlat.eps");
 130
 131   //-----------------------------------------------------------------------------
 132   TCanvas *c5 = new TCanvas("c5","XZ in M1");
 133   c5->SetLogz();
 134   hCellNXZM1->SetTitle("Cell occupancy in module 4");
 135   hCellNXZM1->SetXTitle("X (cells)");
 136   hCellNXZM1->SetYTitle("X (cells)");
 137   hCellNXZM1->Draw("colz");
 138   c5->Print("PHOS_XYM4.eps");
 139
 140   //-----------------------------------------------------------------------------
 141   TCanvas *c6 = new TCanvas("c6","XZ in M2");
 142   c6->SetLogz();
 143   hCellNXZM2->SetTitle("Cell occupancy in module 3");
 144   hCellNXZM2->SetXTitle("X (cells)");
 145   hCellNXZM2->SetYTitle("X (cells)");
 146   hCellNXZM2->Draw("colz");
 147   c6->Print("PHOS_XYM3.eps");
 148
 149   //-----------------------------------------------------------------------------
 150   TCanvas *c7 = new TCanvas("c7","XZ in M3");
 151   c7->SetLogz();
 152   hCellNXZM3->SetTitle("Cell occupancy in module 2");
 153   hCellNXZM3->SetXTitle("X (cells)");
 154   hCellNXZM3->SetYTitle("X (cells)");
 155   hCellNXZM3->Draw("colz");
 156   c7->Print("PHOS_XYM2.eps");
 157
 158   //-----------------------------------------------------------------------------
 159   TCanvas *c8 = new TCanvas("c8","gg mass vs pt, PID=All");
 160   c8->SetLogz();
 161   hPi0All_cen0->SetTitle("PID: All");
 162   hPi0All_cen0->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
 163   hPi0All_cen0->SetYTitle("p_{T} (GeV/c)");
 164   hPi0All_cen0->Draw("colz");
 165   c8->Print("PHOS_MggAll.eps");
 166
 167   //-----------------------------------------------------------------------------
 168   TCanvas *c9 = new TCanvas("c9","gg mass vs pt, PID=Allcore");
 169   c9->SetLogz();
 170   hPi0Allcore_cen0->SetTitle("PID: Allcore");
 171   hPi0Allcore_cen0->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
 172   hPi0Allcore_cen0->SetYTitle("p_{T} (GeV/c)");
 173   hPi0Allcore_cen0->Draw("colz");
 174   c9->Print("PHOS_MggAllcore.eps");
 175
 176   //-----------------------------------------------------------------------------
 177   TCanvas *c10 = new TCanvas("c10","gg mass in M1");
 178   hM1 = hPi0M11->ProjectionX("m1",31,400);
 179   hM1->SetAxisRange(0.,0.29.);
 180   hM1->SetTitle("#gamma#gamma in module 1");
 181   Fit1Pi0(hM1,2);
 182   hM1->Draw();
 183   TPaveText *txtM1 = new TPaveText(0.6,0.7,0.89,0.89,"NDC");
 184   txtM1->SetFillColor(kWhite);
 185   txtM1->SetBorderSize(0);
 186   txtM1->AddText("p_{T}>3 GeV/c");
 187   txtM1->AddText(Form("M_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 188                       hM1->GetFunction("fitfun")->GetParameter(1)*1000,
 189                       hM1->GetFunction("fitfun")->GetParError (1)*1000));
 190   txtM1->AddText(Form("#sigma_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 191                       hM1->GetFunction("fitfun")->GetParameter(2)*1000,
 192                       hM1->GetFunction("fitfun")->GetParError (2)*1000));
 193   txtM1->Draw();
 194   c10->Print("PHOS_MggM1.eps");
 195
 196   //-----------------------------------------------------------------------------
 197   TCanvas *c11 = new TCanvas("c11","gg mass in M2");
 198   hM2 = hPi0M22->ProjectionX("m2",31,400);
 199   hM2->SetAxisRange(0.,0.29.);
 200   hM2->SetTitle("#gamma#gamma in module 2");
 201   Fit1Pi0(hM2,2);
 202   hM2->Draw();
 203   TPaveText *txtM2 = new TPaveText(0.6,0.7,0.89,0.89,"NDC");
 204   txtM2->SetFillColor(kWhite);
 205   txtM2->SetBorderSize(0);
 206   txtM2->AddText("p_{T}>3 GeV/c");
 207   txtM2->AddText(Form("M_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 208                       hM2->GetFunction("fitfun")->GetParameter(1)*1000,
 209                       hM2->GetFunction("fitfun")->GetParError (1)*1000));
 210   txtM2->AddText(Form("#sigma_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 211                       hM2->GetFunction("fitfun")->GetParameter(2)*1000,
 212                       hM2->GetFunction("fitfun")->GetParError (2)*1000));
 213   txtM2->Draw();
 214   c11->Print("PHOS_MggM2.eps");
 215
 216   //-----------------------------------------------------------------------------
 217   TCanvas *c12 = new TCanvas("c12","gg mass in M3");
 218   hM3 = hPi0M33->ProjectionX("m3",31,400);
 219   hM3->SetAxisRange(0.,0.29.);
 220   hM3->SetTitle("#gamma#gamma in module 3");
 221   Fit1Pi0(hM3,2);
 222   hM3->Draw();
 223   TPaveText *txtM3 = new TPaveText(0.6,0.7,0.89,0.89,"NDC");
 224   txtM3->SetFillColor(kWhite);
 225   txtM3->SetBorderSize(0);
 226   txtM3->AddText("p_{T}>3 GeV/c");
 227   txtM3->AddText(Form("M_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 228                       hM3->GetFunction("fitfun")->GetParameter(1)*1000,
 229                       hM3->GetFunction("fitfun")->GetParError (1)*1000));
 230   txtM3->AddText(Form("#sigma_{#pi^{0}} = %.2f #pm %.2f MeV/c^{2}",
 231                       hM3->GetFunction("fitfun")->GetParameter(2)*1000,
 232                       hM3->GetFunction("fitfun")->GetParError (2)*1000));
 233   txtM3->Draw();
 234   c12->Print("PHOS_MggM3.eps");
 235 }
 236
 237 //=============================================================================
 238 Fit1Pi0(const TH1D *hMass = 0,
 239        const Int_t polN = 1)
 240 {
 241   // This script takes a 2D histogram hMassPt with invariant mass vs
 242   // pt of cluster pairs:
 243   // - slices them along pt bins,
 244   // - fits 1D invariant mass specrta by Gaussian+polynomial
 245   // - calculates the number of pi0's as an integral of Gaussian
 246   // - calculates background under pi0 as an integral of polynomial
 247
 248   TF1 *fitfun = 0;
 249   if      (polN == 1)
 250     fitfun = new TF1("fitfun",pi0massP1,0.1,0.7,5);
 251   else if (polN == 2)
 252     fitfun = new TF1("fitfun",pi0massP2,0.1,0.7,6);
 253   else if (polN == 101)
 254     fitfun = new TF1("fitfun",CombiBG,0.0,1.0,6);
 255   fitfun->SetLineColor(kRed);
 256   fitfun->SetLineWidth(2);
 257
 258   if (polN <= 2) {
 259     fitfun->SetParName(0,"A");
 260     fitfun->SetParName(1,"m_{0}");
 261     fitfun->SetParName(2,"#sigma");
 262     fitfun->SetParName(3,"a_{0}");
 263     fitfun->SetParName(4,"a_{1}");
 264     if (polN == 2)
 265       fitfun->SetParName(5,"a_{2}");
 266
 267     fitfun->SetParLimits(0,  1.000,1.e+5);
 268     fitfun->SetParLimits(1,  0.12,0.14);
 269     fitfun->SetParLimits(2,  0.003,0.020);
 270   }
 271   else if (polN == 101) {
 272     fitfun->SetParName(0,"a_{0}");
 273     fitfun->SetParName(1,"a_{1}");
 274     fitfun->SetParName(2,"a_{2}");
 275     fitfun->SetParName(3,"a_{3}");
 276     fitfun->SetParName(4,"a_{4}");
 277     fitfun->SetParName(5,"a_{5}");
 278   }
 279
 280   Int_t   nM     = hMass->GetNbinsX();
 281   Float_t mMin   = hMass->GetXaxis()->GetXmin();
 282   Float_t mMax   = hMass->GetXaxis()->GetXmax();
 283   Float_t mStep  = (mMax-mMin)/nM;
 284
 285   hMass->SetXTitle("M_{#gamma#gamma}, GeV/c");
 286   hMass->SetAxisRange(0.01,1.0);
 287   Float_t nPi0Max = hMass->Integral(hMass->FindBin(0.30),
 288                                     hMass->FindBin(0.40));
 289   for (Int_t iM=1; iM<nM; iM++)
 290     if (hMass->GetBinContent(iM) == 0) hMass->SetBinError(iM,1.0);
 291   hMass->SetMinimum(0.01);
 292   if      (polN == 1)
 293     fitfun->SetParameters(nPi0Max/4,0.135,0.010,1.,0.);
 294   else if (polN == 2)
 295     fitfun->SetParameters(nPi0Max/4,0.135,0.010,1.,0.,0.);
 296 //   else if (polN == 101)
 297 //     fitfun->SetParameters(1.,1.,100.,1.,1.,10.);
 298   Double_t mFitMin, mFitMax;
 299   mFitMin=0.10;
 300   mFitMax=0.18;
 301   hMass->Fit("fitfun","Q","",mFitMin,mFitMax);
 302   hMass->SetAxisRange(mFitMin,mFitMax);
 303   fitfun->SetNpx(1./mStep);
 304 //   TH1F *bgFun = fitfun->GetHistogram();
 305 //   hMass->Add(bgFun,-1.);
 306 }
 307
 308 //-----------------------------------------------------------------------------
 309 Double_t pi0massP2(Double_t *x, Double_t *par)
 310 {
 311   Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
 312                                        (2*par[2]*par[2]) );
 313   Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
 314   return gaus+back;
 315 }
 316
 317 //-----------------------------------------------------------------------------
 318 Double_t pi0massP1(Double_t *x, Double_t *par)
 319 {
 320   Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
 321                                        (2*par[2]*par[2]) );
 322   Double_t back = par[3] + par[4]*x[0];
 323   return gaus+back;
 324 }
 325
 326 //-----------------------------------------------------------------------------
 327 Double_t CombiBG(Double_t *x, Double_t *par)
 328 {
 329   if (x[0] > 0.12 && x[0] < 0.15) {
 330     TF1::RejectPoint();
 331     return 0;
 332   }
 333   Double_t back = par[0] + par[1]*x[0]
 334     + par[2]*TMath::Power(x[0],2)
 335     + par[3]*TMath::Power(x[0],3)
 336     + par[4]*TMath::Power(x[0],4)
 337     + par[5]*TMath::Power(x[0],5);
 338   return back;
 339 }
 340
 341 //-----------------------------------------------------------------------------
 342 Double_t CombiBGExp(Double_t *x, Double_t *par)
 343 {
 344   if (x[0] > 0.12 && x[0] < 0.15) {
 345     TF1::RejectPoint();
 346     return 0;
 347   }
 348   Double_t back = par[0] * TMath::Power(x[0],1.5) * TMath::Exp(-par[1]*x[0]);
 349   return back;
 350 }
 351