PWG3/hfe/macros/CompareCorrectedSpectra.C

   1 void CompareCorrectedSpectra(const char *datafilea, const char *datafileb){
   2
   3   gROOT->SetStyle("Plain");
   4   gStyle->SetPalette(1);
   5   gStyle->SetOptStat(1111);
   6   gStyle->SetPadBorderMode(0);
   7   gStyle->SetCanvasColor(10);
   8   gStyle->SetPadLeftMargin(0.13);
   9   gStyle->SetPadRightMargin(0.13);
  10
  11   // Take files
  12
  13   TFile *filea = TFile::Open(datafilea);
  14   TGraphErrors *correctedSpectrum_a = (TGraphErrors *) filea->Get("AlltogetherSpectrum");
  15
  16   TH1D *histo = (TH1D*) filea->Get("RatioUnfoldingAlltogetherSpectrum");
  17   histo->SetName("historatio");
  18   TH1D *histoa = (TH1D*) histo->Clone();
  19   histoa->Sumw2();
  20   histoa->SetName("a");
  21   TH1D *histob = (TH1D*) histo->Clone();
  22   histob->Sumw2();
  23   histob->SetName("b");
  24
  25
  26   TFile *fileb = TFile::Open(datafileb);
  27   TGraphErrors *correctedSpectrum_b = (TGraphErrors *) fileb->Get("AlltogetherSpectrum");
  28
  29   // Style
  30
  31  correctedSpectrum_a->SetMarkerStyle(24);
  32  correctedSpectrum_a->SetMarkerColor(1);
  33  correctedSpectrum_a->SetLineColor(1);
  34
  35  correctedSpectrum_b->SetMarkerStyle(27);
  36  correctedSpectrum_b->SetMarkerColor(4);
  37  correctedSpectrum_b->SetLineColor(4);
  38
  39   //
  40
  41   TCanvas *c1 = new TCanvas("CorrectedSpectrum","CorrectedSpectrum",800,800);
  42   c1->cd(1);
  43   TH1D *total = new TH1D("total","",1,0.38,4.3);
  44   total->SetMaximum(1.0);
  45   total->SetMinimum(1.0e-09);
  46   total->SetXTitle("p_{T} [GeV/c]");
  47   total->SetYTitle("1/2#pip_{T} d^{2}N/dp_{T}dy [GeV/c]^{-2}, |#eta| < 0.8");
  48   total->SetTitleOffset(1.5,"Y");
  49   total->Draw();
  50   gPad->SetLeftMargin(0.13);
  51   gPad->SetLogy();
  52   gPad->SetTicks();
  53   gPad->SetFillColor(10);
  54   gPad->SetFrameFillColor(10);
  55   gStyle->SetOptStat(0);
  56   gStyle->SetOptFit(1111);
  57   gStyle->SetOptTitle(0);
  58   correctedSpectrum_a->Draw("Psame");
  59   correctedSpectrum_b->Draw("Psame");
  60   TLegend *leg = new TLegend(0.25,0.825,0.48,0.875);
  61   leg->AddEntry(correctedSpectrum_a,"Corrected spectrum a","lep");
  62   leg->AddEntry(correctedSpectrum_b,"Corrected spectrum b","lep");
  63   leg->SetFillStyle(0);
  64   leg->Draw("same");
  65
  66   // Ratio
  67
  68   Double_t x[300];
  69   Double_t ry[300];
  70   Double_t y[300];
  71   Double_t rex[300];
  72   Double_t rey[300];
  73
  74   double xa,ya,xb,yb,eya,exa,eyb,exb;
  75   Int_t npointsa = correctedSpectrum_a->GetN();
  76   Int_t npointsb = correctedSpectrum_b->GetN();
  77   if(npointsa != npointsb) {
  78     printf("Problem the two spectra have not the same number of points");
  79     return;
  80   }
  81   for(Int_t k = 0; k < npointsa; k++){
  82     correctedSpectrum_a->GetPoint(k,xa,ya);
  83     correctedSpectrum_b->GetPoint(k,xb,yb);
  84     //
  85     Double_t centerhisto = histoa->GetBinCenter(k+1);
  86     //printf("bin center %f and center %f\n",centerhisto,xa);
  87     histoa->SetBinContent(k+1,ya);
  88     histob->SetBinContent(k+1,yb);
  89     //
  90     if(TMath::Abs(xa-xb) > 0.0001) {
  91       printf("Problem the two spectra have not the same number of points");
  92       return;
  93     }
  94     eya = correctedSpectrum_a->GetErrorY(k);
  95     exa = correctedSpectrum_a->GetErrorX(k);
  96     eyb = correctedSpectrum_b->GetErrorY(k);
  97     exb = correctedSpectrum_b->GetErrorX(k);
  98     x[k] = xa;
  99     rex[k] = exa;
 100     if(yb > 0.0) y[k] = ya/yb;
 101     ry[k] = ya-yb;
 102     Double_t error = 0.0;
 103     if((yb > 0.0) && (ya > 0.0)) {
 104       error = TMath::Sqrt(TMath::Abs(eya*eya-eyb*eyb));
 105     }
 106     rey[k] = error;
 107     histoa->SetBinError(k+1,eya);
 108     histob->SetBinError(k+1,eyb);
 109   }
 110
 111   histo->Sumw2();
 112   histo->Divide(histoa,histob,1.0,1.0,"B");
 113
 114   TGraph *gratio = new TGraph(npointsa,&x[0],&y[0]);
 115   gratio->SetName("RatioOfSpectra");
 116
 117   TGraphErrors *gdiff = new TGraphErrors(npointsa,&x[0],&ry[0],&rex[0],&rey[0]);
 118   gdiff->SetName("Difference");
 119
 120   TCanvas *c2 = new TCanvas("ratio","ratio",800,800);
 121   c2->cd(1);
 122   TH1D *ratio = new TH1D("ratio","",1,0.38,10.0);
 123   ratio->SetMaximum(1.5);
 124   ratio->SetMinimum(0.5);
 125   ratio->SetXTitle("p_{T} [GeV/c]");
 126   ratio->SetYTitle("Ratio of corrected spectra a/b");
 127   ratio->SetTitleOffset(1.5,"Y");
 128   //ratio->Draw();
 129   histo->SetLineColor(2);
 130   histo->Draw();
 131   gPad->SetLeftMargin(0.13);
 132   gPad->SetTicks();
 133   //gPad->SetGridx();
 134   //gPad->SetGridy();
 135   gPad->SetFillColor(10);
 136   gPad->SetFrameFillColor(10);
 137   gStyle->SetOptStat(0);
 138   gStyle->SetOptFit(1111);
 139   gStyle->SetOptTitle(0);
 140   gratio->SetMarkerStyle(24);
 141   gratio->SetMarkerColor(4);
 142   gratio->SetLineColor(4);
 143   gratio->Draw("P");
 144
 145   TCanvas *c3 = new TCanvas("Difference","Difference",800,800);
 146   c3->cd(1);
 147   TH1D *diff = new TH1D("difference","",1,0.38,4.3);
 148   diff->SetMaximum(1.0e-04);
 149   diff->SetMinimum(-1.0e-04);
 150   diff->SetXTitle("p_{T} [GeV/c]");
 151   diff->SetYTitle("Difference of corrected spectra a-b");
 152   diff->SetTitleOffset(1.5,"Y");
 153   diff->Draw();
 154   gPad->SetLeftMargin(0.13);
 155   gPad->SetTicks();
 156   //gPad->SetGridx();
 157   //gPad->SetGridy();
 158   gPad->SetFillColor(10);
 159   gPad->SetFrameFillColor(10);
 160   gStyle->SetOptStat(0);
 161   gStyle->SetOptFit(1111);
 162   gStyle->SetOptTitle(0);
 163   gdiff->SetMarkerStyle(24);
 164   gdiff->SetMarkerColor(4);
 165   gdiff->SetLineColor(4);
 166   gdiff->Draw("P");
 167
 168
 169  return;
 170
 171 }
 172